Ensemble Learning - Boosting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boosting\n",
    "\n",
    "Boosting is an ensemble learning method that uses a group of weak learners to make a strong learner. Models are trained sequentially and the errors of earlier learners indicate hard examples which the later learners will focus on. Thus boosting algorithms are said to train models in a forward stage-wise manner.\n",
    "\n",
    "## Adaptive Boosting (AdaBoost)\n",
    "\n",
    "AdaBoost typically uses decision stumps (decision trees with a single split) as the weak classifiers although the idea can be applied to any type of weak classifier. Whereas in a random forest, all decision trees get an equal vote in the final prediction, the decision stumps' votes are weighted by how well they generally perform. Each stump is made by taking the previous stump's mistakes into account. It gives more weight to the data points where the previous classifier made mistakes (hence \"boosting\" these data points).\n",
    "\n",
    "### Binary AdaBoost\n",
    "\n",
    "The general formula for AdaBoost for binary classification (where $Y \\in \\{-1,+1\\}$) is:\n",
    "\n",
    "$$H(x) = sign\\left(\\displaystyle{\\sum_{t=1}^T\\alpha_th_t(x)}\\right)$$\n",
    "\n",
    "where:\n",
    "- $H$ is the AdaBoost classifier\n",
    "- $T$ is the number of weak classifiers\n",
    "- $h_t$ are the weak classifiers (typically decision stumps)\n",
    "- $\\alpha_t$ are the coefficients that correspond to the importance of the classifiers' votes\n",
    "\n",
    "$$\\alpha_t = \\displaystyle{\\frac{1}{2}}\\ln\\left(\\frac{1-\\epsilon_t}{\\epsilon_t}\\right)$$<br>\n",
    "$$\\epsilon_t = \\displaystyle{\\sum_{i\\text{ where }h_t(x)\\text{ is wrong}}w_i}$$<br>\n",
    "$$\\displaystyle{\\sum_{i=1}^Nw_i}=1$$\n",
    "\n",
    "where:\n",
    "- $w_i$ are the weights of the data points\n",
    "    - initially $w_i = \\frac{1}{N}$, meaning all data points are equally important\n",
    "    - sample weights adjust as new classifiers are trained:\n",
    "        - if $h_{t-1}$ got data point $i$ wrong, $w_i$ increases in $h_t$\n",
    "        - if $h_{t-1}$ got data point $i$ right, $w_i$ decreases in $h_t$\n",
    "        - sample weights should always add to $1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEXCAYAAACtTzM+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3Xd4XOWZ9/HvrS5Z3ZZlVeRewTYWxSRUmxow7yakB8LuBhOy2YQN2eQlhIR3U3dTls2SgpOwmBZgCdkEQkkoDqHYYGOvMe5NttwkN1ku6vf7x4yEbKuMbWnOSPP7XNdc1sw5c577yDPz0/M8Z84xd0dERCQh6AJERCQ2KBBERARQIIiISJgCQUREAAWCiIiEKRBERARQIIiISJgCYZAzs81mNjvoOtqZ2XgzW2pm9Wb2hVjbXoRtvmtmF3XXfhA1Rctgfz3FOwVCjDOz583sX7p4/Foz22lmSUHUdQq+Aixw9yx3/0msbS/8gXck/AGz38xeN7PPmlnHe8XdJ7v7gh7a7+t9PNH6u/3A1utJeqJAiH33A9ebmR3z+PXAw+7eEv2STslpwLtBb6+XD75r3D0rvO3vA18Ffn0C7fdHTX3lfvR6ku64u24xfAPSgTrggk6P5QENwNTw/f8LbADqgZXA33RadzMwO/yzA2M6Lbsf+Han+8XAb4FaYBPwhU7LvgpsC7exBpjVTb0TgQXAfkJv1Dmdlr0EtIZrPwiM6+L5ScCd4br3AJ8g9FfgHV2se9z2eml/c3g/lgONQFIX2+z4fXV67GygDZjSeZ1u2u/qsZ5+r8fV1NP6nZ7z5fBz6oDHgDTgwXCdR8Jtf6U/X0+9vaYi2I9Tek119bvu4rk5wG+A3eF2lgEJQb+vY/UWeAG6RfCfBL8EftXp/s3Ask73Pxx+8yUAHwUOAUXhZR1v4F7evAnAEuAbQAowCtgIXA6MB7YCxeF1K4DRXdSZDKwHvhbexiXhN+H4TussAD7Tw75+H3gh/Ea+BlgV/lDK6mb9ju311n74d7EMKAPSu9neUR94nR7fAtzSxe/0uP05pqZuf69d1dTb+p2e82b4/zw//Dv6bE/198frqafXVAT73Sevqa5+/8c8/xfAQ8CQcE1Tgn4/x/JNQ0YDw3zgw2aWHr5/Q/gxANz9v919u7u3uftjwDpCf9WeiLOAAnf/F3dvcveNhD44Pkbor7BUYJKZJbv7Znff0MU2zgUyge+Ht/ES8DTw8UgKMLNs4FZgrrvXAYuACcBD7l4fwSYiaf8n7r7V3Y9EUlMn2wl9+J6onn6vXdUUyfrtz9nu7nuBp4BpJ1BT0K8niNJrCmgO1384vD8rTnA/4ooCYQBw91cJdbuvNbNRhN5sj7QvN7MbzGxZeBJ0PzAFGHaCzZwGFLdvI7ydrwGF7r6e0Af1XUCNmT1qZsVdbKMY2OrubZ0eqwJKIqzhEmBt+MMDQn8R1gH/GeHzI2l/a4TbOlYJsPckntft77WbmiJZH2Bnp58PE/rQjEjQr6dwDdF6Ta0mNLx2yMxuPsF9iDsD7YiCePYAob/kxgN/cvddAGZ2GqG/vGYBb7h7q5ktA46dNITQB0dGp/sjgOrwz1uBTe4+tqvG3f0R4JHwX/H3Av9KaCKys+1AmZkldHoDlwNrI9zH4vA22s0FtkXYO4i0/RM+37uZnUXoA+jVE30uvfxeu6gpkvV7Eun+9cXrCbp/TfW6H/39mjKzOcDngTPdfV1v64t6CAPJA4QmMm+iU/ee0NioE/qLDzP7W0J/0XVlGfAJM0s0syuACzstexM4YGZfNbP08DpTzOys8LHel5hZKqEJvCOEuvzHWkRovPkrZpYcPlb/GuDRCPexGphmZkVmdg6hD4fhZpYS4fNPtf2jmFm2mV0dfv5D7v7OSWym299rH61/rF2Exut70xevJ+j+NdXjfkTpNTWF0KR1e9iVm1leBM+LWwqEAcLdNwOvE3rD/qHT4yuBHwFvEHrhnw681s1mvkjozbQf+CTwP5220xpeNo3QESG7gV8RmtxNJTTZu5vQUMVwQt3/Y2tsAuYAV4bX/Rlwg7uvjnA3nwP+RGiS9DfABwl94LwUyZP7oP12T5lZPaG/cu8Afgz87Qluo72mnn6vp7x+F74HfD08TPPlHurazKm/nqCb11QE+xGN19R/AU3AVjOrA55EoyI9MnddMU1ERNRDEBGRsEC7T2a2mdAxxa1Ai7tXBlmPiEg8i4XxtIvdfXfQRYiIxDsNGYmICBDwpLKZbQL2ETrM7V53n9fT+sOGDfOKiopolCYiMmgsWbJkt7sX9LZe0ENG73P37WY2HPizma1291c6r2Bmcwl9QYny8nIWL14cRJ0iIgOWmVVFsl6gQ0buvj38bw3wO7o4X4q7z3P3SnevLCjoNeBEROQkBRYIZjbEzLLafwYuA3TiKRGRgAQ5ZFQI/C58nY4k4BF3fy7AekRE4lpggRA+o+XUoNoXEZGj6bBTEREBFAgiIhKmQBARESBOAuHFVbv42YL1QZchIhLT4iIQFqyp5ZevbOx9RRGROBYXgZCYYLTpsg8iIj2Ki0AwgzYlgohIj+IiEBLNaNWV4UREehQfgZBgtCkQRER6FBeBYGa0tQVdhYhIbIuLQEhMQENGIiK9iI9AMA0ZiYj0Ji4CwcxwhyCvDiciEuviIhASEwyAVh16KiLSrbgKBOWBiEj34iIQQtfgQfMIIiI9CDwQzCzRzJaa2dP91UaiachIRKQ3gQcC8EVgVX828N6QkQJBRKQ7gQaCmZUCHwB+1c/tAOjLaSIiPQi6h3A38BWg249qM5trZovNbHFtbe1JNZIYnkPQl9NERLoXWCCY2dVAjbsv6Wk9d5/n7pXuXllQUHBSbWnISESkd0H2EN4HzDGzzcCjwCVm9lB/NPTekJECQUSkO4EFgrvf7u6l7l4BfAx4yd0/1R9tdXwxTT0EEZFuBT2HEBXth52qgyAi0r2koAsAcPcFwIL+2n7HF9OUCCIi3YqPHoLOZSQi0qu4CgQdZSQi0r24CISOo4wUCCIi3YqLQHjvXEYBFyIiEsPiIxDCe6kegohI9+IiEExnOxUR6VVcBEKi5hBERHoVH4GgK6aJiPQqLgIhoeN7CJpVFhHpTlwEQlpSaDcbmhUIIiLdiYtASE9JBOBIU2vAlYiIxK74CITkcCA0KxBERLoTF4GQpkAQEelVXARC+5BRgwJBRKRb8REIyQoEEZHeBHlN5TQze9PM/tfM3jWz/9dfbXUMGTXpKCMRke4EeYGcRuASdz9oZsnAq2b2rLsv7OuGEhOM1KQEDjW19PWmRUQGjcACwd0dOBi+mxy+9dt3iXMzktl/uKm/Ni8iMuAFOodgZolmtgyoAf7s7ou6WGeumS02s8W1tbUn3VZeRgr7DjefQrUiIoNboIHg7q3uPg0oBc42syldrDPP3SvdvbKgoOCk21IPQUSkZzFxlJG77wcWAFf0VxvqIYiI9CzIo4wKzCw3/HM6MBtY3V/t5WakqIcgItKDII8yKgLmm1kioWB63N2f7q/G8ocks+9wM+7eccEcERF5T5BHGS0HpkervbyMFFrbnAMNLeSkJ0erWRGRASMm5hCiIS8jBYC9hzRsJCLSlbgJhOLcdAC27z8ScCUiIrEpbgKhNC8UCNX7DgdciYhIbIqbQCjKSSMxwajepx6CiEhX4iYQkhITGJGdpkAQEelG3AQCQFl+Opt2Hwq6DBGRmBRXgTBhRDZrd9XT2tZv59ATERmw4ioQJhVlc7iplao96iWIiBwrvgKhOBuAlTsOBFyJiEjsiatAGDM8k8QEY5UCQUTkOHEVCGnJiYwdnsmyrfuDLkVEJObEVSAAzBw9lMWb99HQ3Bp0KSIiMSXuAuF9o4fR2NLG21v2BV2KiEhMibtAOGdUPokJxmvrdwddiohITIm7QMhKS+asijyeW7ETd30fQUSkXdwFAsAHzihmQ+0h1u46GHQpIiIxI8hLaJaZ2ctmtsrM3jWzL0ar7SsmjyDB4H+WbYtWkyIiMS/IHkILcJu7TwTOBf7BzCZFo+GCrFRmTSzk8be20tiio41ERCDAQHD3He7+dvjnemAVUBKt9m+YeRp7DjXxzDs7otWkiEhMi4k5BDOrIHR95UVdLJtrZovNbHFtbW2ftfm+0cMYVTCEX76yiTad7E5EJPhAMLNM4LfAre5+3Dkl3H2eu1e6e2VBQUGftZuQYHz+4jGs3HGA59/d2WfbFREZqAINBDNLJhQGD7v7k9Fu/9ppJYwuGMIP/7SG5ta2aDcvIhJTgjzKyIBfA6vc/cdB1JCYYNx+5UQ21B7iV3/dFEQJIiIxI8gewvuA64FLzGxZ+HZVtIuYPamQyycX8h8vrmXLnsPRbl5EJGYEeZTRq+5u7n6Gu08L354Jopa75kwmKSGBf3p8mYaORCRuBT6pHAuKctL57gdPZ0nVPn7857VBlyMiEggFQticqcV8/Owyfr5gA7/XN5hFJA4pEDq5a85kzh6Zzz//93IWbdwTdDkiIlGlQOgkNSmRedfPoCw/nc/MX8ySKl0zQUTihwLhGLkZKTz0mXMYmpnCDb9exFub9wZdkohIVCgQulCUk85jN8+kMCeN63+9iGd1viMRiQMKhG4UZqfx+M0zmVSUzS0Pv81PX16vC+qIyKCmQOjBsMxUHrnpXK6dVswPnl/DzQ8uYf/hpqDLEhHpFwqEXqQlJ3L3R6dxx1UTeXlNDVf+x19ZqCOQRGQQUiBEwMy46YJRPHnL+0hNSuDjv1zIN3+/gvqG5qBLExHpMwqEE3B6aQ5Pf+F8Pj2zggcWVnHpj1/h+Xd3am5BRAYFBcIJykxN4q45k3nylvPIzUjm5geX8KlfL2LFtrqgSxMROSUKhJM0vTyPp/7x/Xzzmkms3H6Aa+55lS89voyte3XGVBEZmGwgDXdUVlb64sWLgy7jOHVHmvnZgvX812ubaWtz/mZ6CZ+7eAwjhw0JujQREcxsibtX9rqeAqHv7Kg7wr1/2chv3txCc2sb10wt5uYLRjOpODvo0kQkjg2IQDCz+4CrgRp3n9Lb+rEeCO1q6hv49V838eDCKg43tXJ2RT43vq+CyyYVkpSoUToRia6BEggXAAeBBwZTILSrO9zM44u3Mv+NzVTvO0JRThqfPKec62aUMSInLejyRCRODIhAADCzCuDpwRgI7VrbnJdW1zD/9c28un43CQbnjy3guhmlXDqpkLTkxKBLFJFBbNAEgpnNBeYClJeXz6iqqopecf1g8+5D/Pbtan67pJrtdQ3kpCdzzdQi5kwtofK0PBISLOgSRWSQGTSB0NlA7SF0pbXNeX3Dbp5YUs1zK3bS2NJGYXYqV04p4uozijizXOEgIn0j0kBIikYxcrzEBOP8sQWcP7aAg40tvLhqF39cvoNH3tzC/a9vZkR2GldMGcGlkwo5qyKflCRNRotI/1IgxIDM1CSunVbCtdNKqG9o5qXVNTzdKRyyUpO4YFwBsyYO5+Lxw8kbkhJ0ySIyCAV9lNFvgIuAYcAu4Jvu/uvu1h9MQ0aRONzUwmvr9/Diql28uLqG2vpGEgzOLM/j4gnDef+YYUwpySFRQ0si0oMBM4dwIuItEDpra3NWbK/jhVU1vLR6Fyu2HQAgJz2Z80YP5f1jh3H+mALKh2YEXKmIxJpIA0ED0wNEQoJxRmkuX7p0HE//4/ks/vps/uNj07hsUiHLtu7njt+t4IIfvMz5//YStz/5Dr9fto2ddQ1Blz3gHDlyhAsvvJDW1lbmz5/P2LFjGTt2LPPnz+9y/bvuuouSkhKmTZvGtGnTeOaZZ7pcr7ttzZ49m3379vXLvoicKPUQBgF3Z+PuQ7y6bjevrt/Nwg17qG9sAaA8P4OzR+Zz9sh8zhmZT3l+BmYaYurOT3/6U1paWrj++uuprKxk8eLFmBkzZsxgyZIl5OXlHbX+XXfdRWZmJl/+8pe73ebevXu73db8+fOprq7mjjvu6O9dkzimo4ziiJkxuiCT0QWZfPq8Clpa21i9s55Fm/by5qbQHMQTS6oBKMxO5ayKUDhML89jwogsnU6jk4cffphHHnmE559/nksvvZT8/HwALr30Up577jk+/vGPn/A2e9rWnDlzOP/88xUIEhN6DQQzm+HuS6JRjPSNpMQEppTkMKUkh79//0jcnfU1B8MBEbo9vXwHAOnJiZxemsP08lyml+VxZnkuw7Pj87QaTU1NbNy4kYqKCp544gnKyso6lpWWlrJt27Yun3fPPffwwAMPUFlZyY9+9KPjehHbtm3rdlt5eXk0NjayZ88ehg4d2g97JRK5SHoInzaz24EfuvtCADP7sbt/qX9Lk75iZowtzGJsYRafOvc03J3qfUdYunU/S7fsY+mW/dz36iaaWzcCUJKbzrTyXKaX5TK9PI/JxdlxcXqN3bt3k5ubC9DlVfC6Gmq75ZZbuPPOOzEz7rzzTm677Tbuu+++o9bpbVvDhw9n+/btCgQJXCSBUAPMAZ40s3ogBXijX6uSfmVmlOVnUJafwZypxQA0NLeycscBlm55LyT+GO5FJCYYY4dncnpJDqeXhnoek4oGX0ikp6fT0BCaiC8tLWXBggUdy6qrq7nooouOe05hYWHHzzfddBNXX331cev0tq2GhgbS09NPuX6RU9XrpLKZrQamunujmRUD3wOWuvvd0SiwM00qR1fNgQaWbt3PO9V1vLOtjhXb6thzqAl4LySmlORwesl7IZGeMrBDoqysjHXr1nH48GFmzJjB22+/DcCZZ57JkiVLOuYB2u3YsYOioiIA/v3f/51Fixbx6KOPsm3bNm644QZefPFF9u7d2+223J3S0lKqqqpIStKUnvSPvpxU3gqMBFa7+3ZCQ0irgKgHgkTX8Ow0Lp88gssnjwBCQx876ho6wmF5dR0vr67pmLBOTDDGFGQyuSSbSUXZTAzf8gfQN6svu+wyXn31VWbPns2dd97JWWedBcA3vvGNjjD4zGc+w2c/+1kqKyv5yle+wrJlyzAzKioquPfee4FQULR/wOfn53e7rSVLlnDuuecqDCQmRNJDmAT8N/B2+FYKXBhJ2vQ19RBiz7Eh8c62OlZuP0BNfWPHOoXZqUwY0R4QWUwqymbksCExeXTT0qVL+fGPf8yDDz54Stu55557KC8vZ86cOT2u98UvfpE5c+Ywa9asU2pPpCd91kNw95VmdiYwG5gO7ASuPfUSZTAwM4pz0ynOTe/oSQDsOdjIqh31rNpxgFU7DrByxwFe37Cb5tbQHyCpSQmMK8xiYlFWR09i4ohscjKSg9oVAKZPn87FF19Ma2sriYknP/z1+c9/PqL1pkyZojCQmKEvpknUNLW0sb7mYEdIrNp5gFU76tkbnpcAKMpJY2xhFuMLMxlXmMX4EVmMGZ5JRoqGVEROlr6YJjEnJSmBScXZTCrO7njM3amtb2TljlA4rNtVz5pd9czfuIemljYAzKAsL4NxhVmMK8xk/IgsxhVmMapgCKlJA3sSWySWKBAkUGbG8Ow0hmencdH44R2Pt7Y5VXsOsXbXQdaGQ2LdrnoWrKmhpS3Uq01MMCqGZnQERPutYmhGTM5PiMQ6BYLEpMQEY1RBJqMKMrliyntzE00tbWzafYi1u+o7bqt21PPsip20j36mJCYwqmAIY4ZnMnZ4aMhpzPBMKoZlqEch0gMFggwoKUkJjB8RmlvorKG5lfU1nXsTB1leXccf39nRERSJCcZp+RkdATG2MJMxBVmMHj5EcxQiBBwIZnYF8B9AIvArd/9+kPXIwJWWnNhx/qbOGppb2VB7kPU1793W1RzkpdXvDT1B6HQdoR5F5nuBMTyr16OeHn7nYe548Q621G2hPKec78z6Dp88/ZP9so8i/S2wQDCzROCnwKVANfCWmf3B3VcGVZMMPmnJiUwuzmFy8dFB0dzaRtWeQ0eFROgEgHtoaG7rWG9YZmpHSIR6FKGfC7JSeWTFI8x9ai6Hmw8DUFVXxdyn5gIoFGRACuywUzObCdzl7peH798O4O7f6+45WVlZPmPGjKMe+8hHPsLnPvc5Dh8+zFVXXXXcc2688UZuvPFGdu/ezXXXXXfc8ltuuYWPfvSjbN26leuvv/645bfddhvXXHMNa9as4eabbz5u+de//nVmz57NsmXLuPXWW49b/t3vfpfzzjuP119/na997WvHLb/77ruZNm0aL7zwAt/+9rePW37vvfcyfvx4nnrqKX70ox8dt/zBBx+krKyMxx57jJ///OfHLX/iiScYNmwY999/P/fff/9xy5955hkyMjL42c9+xuOPP37c8vZz8Pzwhz/k6aefPmpZeno6zz77LADf+ta3ePHFF49aPnToUH77298CcPvtt/PGG0efAqu0tJSHHnoIgFtvvZVly5YdtXzcuHHMmzcPgLlz57J27dqjlk+bNo277w59Yf5Tn/oU1dXVRy2fOXMm3/te6OX0oQ99iD179hy1fNasWdx5550AXHnllRw5cgSAxpY2jjS1MKbyQkZd8jHW1RzkT/92C62dehRJCUbqxAlw5giaWjZz6LGXcRqBZlKTUjm39Fy99vTaA07stdfu6quv7rjGRlfn0DrRz72//OUvMX/YaQmh02K0qwbOOXYlM5sLzAVITU2NTmUS11KTEkhNSuHcUUP58nVTAbjo8XyaWts40tQaujW3sq+tkZzWsxnSciHNbZvCz26jramR9TUHeXHVLorf3cnQpEbcQ4fPisSyIHsIHwYud/fPhO9fD5zt7v/Y3XP0xTSJJRV3V1BVV0WCZ5HcVkayl5HcVkZ24jiK089ke6dLmHY+8ql9fmJsYSYVQ4eQkqRDZKV/DYQvplUDZZ3ulwLbA6pF5IR9Z9Z3wnMI9TQmrqSRlWQkZ/Cda+bxydNnUd/QzIbaQ6zbVc/62oOs7+7Ip6EZjO0UEu1XvxvoZ46VgSfIQHgLGGtmI4FtwMeATwRYj8gJaZ847u4oo6y0ZKaV5TKtLPeo5x1pCh35tKH2IOt2HWRdTT3rag7ywqqajnmK9m9nH3XUU2HoOxWZqTpEVvpHoOcyMrOrCJ1GOxG4z92/09P6GjKSwayppY3New51hET7EVAbaw/R1PrekU9FOWlHDTu1Hy6bmzFwTjMu0RXpkJFObicS41pa29iy93DHobGhw2Tr2VBziCPNrR3rHXuI7LjCLCaMyFJQyICYQxCRCCQlJnScxuPyye893tbmbNt/pCMg1u06yPrag/zP0m3UN7Z0rFeYncr4EdlMGJHF+E5nkB1sl0CVU6dAEBmgEhLeuzb2xRPeOzGgu7PzQAOrd9azJnxbvbOehRv2dAw9tZ8YcMKI7I5TgUwYkUVZXgYJCTo+Nl4pEEQGGTOjKCedopx0Lu50BtmW1tAcxepOIfHOttBRT+0yUhIZW5jFhHBPYkL4CncadooPmkMQiXOHGltCJwUMh8SanaETBHa+cFFJbjqTw9eyCJ0KJJuinDRM37YbEDSHICIRGZKaxPTyPKaX53U85u7sPtjE6p0HeHd7+62OP6/a1fEdiryM5I5wCAVFNiOHZZKoIacBS4EgIscxMwqyUinIKuD8sQUdjx9qbGH1znpWbq/rCIr/em1zx9xEenIiE4qymFyczRmloe9gjC5QSAwUGjISkVPS3Bq6VvbKTj2JldsPdBzpNCQldGryqWW5TC3NZWpZDiW56RpuiiINGYlIVCQnJjCxKJuJRdl8KHwy4rY2Z9OeQ/zv1v2hW3Ud93fqSQwdksLUslzOKA0FxbTSXPKGaOI6aAoEEelzCQnWcU6mD55ZCoS+ib1mZz3LqkMhsbx6Py+vqemYkxhdMITK0/KprMijsiKfiqEZ6kVEmYaMRCQwBxtbeKe6jre37GNJVehWd6QZgGGZKZxZnkdlRR4zTsvn9JIcnRn2JGnISERiXmZqEjNHD2Xm6KFAaKhpQ+1B3tq8j8VVe1lStY8/rdwFhK5TcWZ5HueNHsp5Y4ZyRmkuyYkKiL6kHoKIxLSa+gbertrHW5v3sXDjHlbuOIB76Et0Z1XkhwJi9DAmFWfraKZu6OR2IjIo7TvUxKJNe3h9wx7e2LCHdTUHAchOS+LcUUO5aPxwLhpfQHFuesCVxg4NGYnIoJQ3JIUrphRxxZQiAGoONPDGxlA4/HXd7o4hpgkjsjrCYcZpeRpeioB6CCIyaLg762sOsmBNLS+vqeHNTXtpaXOyUpM4f9wwLps0gksmDic7LTnoUqNKQ0YiEvfqG5p5bf0eFqyp4eU1New60EhyovH+McO4ckoRl04qjIvvP8R0IJjZh4G7gInA2e4e0ae8AkFETlZbm7N0636eW7GDZ1fspHrfERITjHNH5fOB04v5wOlF5GQMzp5DrAfCRKANuBf4sgJBRKLJ3Xl3+wGeDYfDxtpDpCQmMHvScD44vZQLxxcMqjmHmJ5UdvdVgL6FKCKBMDOmlOQwpSSHL182nhXbDvDk0mr+sGw7z7yzk/whKcyZWswnzilnXGFW0OVGTaBzCGa2gF56CGY2F5gLUF5ePqOqqipK1YlIvGlubeOVtbU8+fY2/rxyF02tbZwzMp/rZ57G5ZNHDNheQ+BDRmb2AjCii0V3uPvvw+ssQENGIhKD9h5q4vHFW3loYRXV+44wPCuVj59dzg0zT2NoZmrQ5Z2QwAMhEgoEEYl1rW3OX9bW8OAbVby8ppa05AQ+dlY5N10wipIB8uW3mJ5DEBEZKBITjEsmFHLJhELW19Tz8wUbeWhhFQ8trOJvppfwxdljKc3LCLrMPhHUUUZ/A/wnUADsB5a5++W9PU89BBGJBdv2H+GXr2zkkTe3gMMNM0/jHy4eE7PfaRgQQ0YnSoEgIrFk+/4j3P3CWp5YUs2QlCQ+f8kY/u79I2Nu8jnSQIitqkVEBpDi3HT+7bqpPH/rBZwzKp/vPbuaD/zkr7y5aW/QpZ0UBYKIyCkaW5jFrz59Fr+6oZJDja185N43uP3J5RwKX1d6oFAgiIj0kdmTCvnzly7g5gtG8ehbW7nqJ3/l7S37gi4rYgoEEZE+lJGSxO1XTeSxuTNpaXU+/Is3+MVfNjAQ5msVCCIi/eDskfk8e+v5XDF5BN9/djVffHQZR5qG1TeCAAAJhklEQVRagy6rRwoEEZF+kp2WzD2fmM4/Xz6ep5Zv56Pz3mDPwcagy+qWAkFEpB+ZGf9w8RjmXV/Jmp31fOTeN9hZ1xB0WV1SIIiIRMGlkwqZ/3dns+tAIx++9/WYDAUFgohIlJw7aigPfeYc9h1q5ob7FrH/cFPQJR1FgSAiEkXTynKZd/0MNu8+zE0PLKappS3okjooEEREouy8McP44Uem8tbmfXzv2VVBl9NBgSAiEoA5U4v52/dV8F+vbeaPy3cEXQ6gQBARCcztV05kenkuX/vdO9TWB384qgJBRCQgKUkJ/PDDUznS1MpdT70bdDkKBBGRII0uyOQLs8bwx+U7+Mva2kBrCSQQzOwHZrbazJab2e/MLDeIOkREYsHcC0ZTlp/Ovz67mra24M55FFQP4c/AFHc/A1gL3B5QHSIigUtJSuC2S8ezcscBnn4nuAnmQALB3f/k7u0nCl8IlAZRh4hIrJgztZjRBUOY90pwZ0aNhTmEvwOe7W6hmc01s8Vmtri2NtjxNRGR/pKQYPzd+0eyYtsBFlcFcw2FfgsEM3vBzFZ0cbu20zp3AC3Aw91tx93nuXulu1cWFBT0V7kiIoH74PRScjOSuf+1zYG0n9RfG3b32T0tN7NPA1cDs3wgXDlCRKSfpack8n+mlfDIm1s40NBMdlpyVNsP6iijK4CvAnPc/XAQNYiIxKJrpxXT1NLGcyt2Rr3toOYQ7gGygD+b2TIz+0VAdYiIxJRpZbmU5acHEgj9NmTUE3cfE0S7IiKxzsy4ZPxwHl9cTWNLK6lJiVFrOxaOMhIRkU4uGFfAkeZWFm+O7tFGCgQRkRgzc/RQkhONV9ZF91B7BYKISIzJSElicnEOS6v2R7VdBYKISAw6szyP5dv209wavSuqKRBERGLQ1LIcGprbWLfrYNTaVCCIiMSgCSOyAVhXUx+1NhUIIiIxaOSwISQlGGt2KhBEROJaSlICI4cNYa2GjEREZOSwIWzZeyhq7SkQRERiVEleOtv2HYna9REUCCIiMaokN51DTa3UHWmOSnsKBBGRGFWSmw5A9b4jUWlPgSAiEqNK8kKBsG2/AkFEJK6NyE4DoKa+MSrtKRBERGJU3pAUAPYcHMSBYGbfMrPl4Yvj/MnMioOoQ0QkliUnJpCbkczeQ01RaS+oHsIP3P0Md58GPA18I6A6RERiWv6QFPYcHMSB4O4HOt0dAkTnIFsRkQEmOy2ZAw3ROew0kEtoApjZd4AbgDrg4qDqEBGJZVlpSRxoaIlKW/3WQzCzF8xsRRe3awHc/Q53LwMeBj7fw3bmmtliM1tcWxvdqweJiAQtOy2Z+oHeQ3D32RGu+gjwR+Cb3WxnHjAPoLKyUkNLIhJXstKSqB/oPYSemNnYTnfnAKuDqENEJNaFAmGA9xB68X0zGw+0AVXAZwOqQ0QkpqWnJNHQ3Ia7Y2b92lYggeDuHwqiXRGRgSYtOTSQ09jSRlpyYr+2pW8qi4jEsNSkUAg0Nrf1e1sKBBGRGNbeQ2hoae33thQIIiIxTD0EEREB1EMQEZGwNPUQREQEIFU9BBERAToONVUPQUQkzhVkpnLV6SPIzUju97YCO9upiIj0rmLYEH72yRlRaUs9BBERARQIIiISpkAQERFAgSAiImEKBBERARQIIiISpkAQERFAgSAiImHmPnCuW29mtYQuuXkyhgG7+7CcgUD7HB+0z/HhVPb5NHcv6G2lARUIp8LMFrt7ZdB1RJP2OT5on+NDNPZZQ0YiIgIoEEREJCyeAmFe0AUEQPscH7TP8aHf9zlu5hBERKRn8dRDEBGRHigQREQEGISBYGZXmNkaM1tvZv+3i+WpZvZYePkiM6uIfpV9K4J9/pKZrTSz5Wb2opmdFkSdfam3fe603nVm5mY24A9RjGSfzewj4f/rd83skWjX2NcieG2Xm9nLZrY0/Pq+Kog6+4qZ3WdmNWa2opvlZmY/Cf8+lpvZmX1agLsPmhuQCGwARgEpwP8Ck45Z53PAL8I/fwx4LOi6o7DPFwMZ4Z9viYd9Dq+XBbwCLAQqg647Cv/PY4GlQF74/vCg647CPs8Dbgn/PAnYHHTdp7jPFwBnAiu6WX4V8CxgwLnAor5sf7D1EM4G1rv7RndvAh4Frj1mnWuB+eGfnwBmmZlFsca+1us+u/vL7n44fHchUBrlGvtaJP/PAN8C/g1oiGZx/SSSfb4J+Km77wNw95oo19jXItlnB7LDP+cA26NYX59z91eAvT2sci3wgIcsBHLNrKiv2h9sgVACbO10vzr8WJfruHsLUAcMjUp1/SOSfe7s7wn9hTGQ9brPZjYdKHP3p6NZWD+K5P95HDDOzF4zs4VmdkXUqusfkezzXcCnzKwaeAb4x+iUFpgTfb+fkKS+2lCM6Oov/WOPq41knYEk4v0xs08BlcCF/VpR/+txn80sAfh34MZoFRQFkfw/JxEaNrqIUC/wr2Y2xd3393Nt/SWSff44cL+7/8jMZgIPhve5rf/LC0S/fn4Nth5CNVDW6X4px3chO9YxsyRC3cyeumixLpJ9xsxmA3cAc9y9MUq19Zfe9jkLmAIsMLPNhMZa/zDAJ5YjfW3/3t2b3X0TsIZQQAxUkezz3wOPA7j7G0AaoZPADVYRvd9P1mALhLeAsWY20sxSCE0a/+GYdf4AfDr883XASx6erRmget3n8PDJvYTCYKCPK0Mv++zude4+zN0r3L2C0LzJHHdfHEy5fSKS1/b/EDqAADMbRmgIaWNUq+xbkezzFmAWgJlNJBQItVGtMrr+ANwQPtroXKDO3Xf01cYH1ZCRu7eY2eeB5wkdoXCfu79rZv8CLHb3PwC/JtStXE+oZ/Cx4Co+dRHu8w+ATOC/w/PnW9x9TmBFn6II93lQiXCfnwcuM7OVQCvwz+6+J7iqT02E+3wb8Esz+ydCQyc3DuQ/8MzsN4SG/IaF50W+CSQDuPsvCM2TXAWsBw4Df9un7Q/g352IiPShwTZkJCIiJ0mBICIigAJBRETCFAgiIgIoEEREJEyBICIigAJBRETCFAgipyh8zv5l4dui8LmURAYcfTFN5BSZ2TrgfHffGXQtIqdCf8mInLpngHfM7O6gCxE5FYPqXEYi0WZm5xE6JXFR+PoaIgOWeggip+bDwNrwidjMzLJ7fYZIjNIcgsgpMLOzCZ1B14EjwOfcfUmwVYmcHAWCiIgAGjISEZEwBYKIiAAKBBERCVMgiIgIoEAQEZEwBYKIiAAKBBERCfv/9we+CwZXoz4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "eps = np.arange(1e-4,1-1e-4,0.001)\n",
    "y = .5*np.log((1-eps)/eps)\n",
    "\n",
    "plt.plot(eps,y)\n",
    "plt.axhline(y=0,color='k',linestyle='--')\n",
    "plt.scatter(0.5,0,color='g')\n",
    "plt.annotate('(0.5,0)',(0.5,0.5),horizontalalignment='center')\n",
    "plt.title('Values of $\\\\alpha$ for Different Values of $\\epsilon$')\n",
    "plt.xlabel('$\\\\epsilon$')\n",
    "plt.ylabel('$\\\\alpha$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the error rate is half, then the weight for that classifier is $0$ because it tells us nothing of value. As the error approaches $0$, the weight for classifier increases. As the error approaches $1$, the weight becomes negative (whatever the classifier predicts is almost always the opposite of the truth). We generally do not want negative alphas because this breaks AdaBoost. As we will see later, this will have the effect of boosting the correctly classified data points instead of the misclassified ones. Also, weak learners by definition are still supposed to be better than random ($\\epsilon<0.5$ in the binary case).\n",
    "\n",
    "But how do we calculate the new weights?\n",
    "\n",
    "$$w_{t+1}=\\displaystyle{\\frac{w_te^{-\\alpha_th_t(x)y(x)}}{z}}$$\n",
    "\n",
    "where:\n",
    "- $y(x)$ is the true label ($\\{-1,+1\\}$)\n",
    "- $h_t(x)$ is the prediction ($\\{-1,+1\\}$)\n",
    "- z is the normalizing factor (simply the sum of all the weights at $t$)\n",
    "\n",
    "If $h(x) = y(x)$, the product of the two is $+1$. If $h(x)\\neq y(x)$, the product is $-1$.\n",
    "\n",
    "#### A Faster Way to Calculate the New Weights\n",
    "The formula for the new weights turns out to have a nice property that allows us to skip all the exponentials and logarithms.\n",
    "\n",
    "If we substitute $\\alpha$ into the equation, there are two possibilities.\n",
    "\n",
    "##### Correct\n",
    "\n",
    "$\\begin{align*}\n",
    "w_{t+1} &= \\displaystyle{\\frac{w_t}{z}}e^{-\\frac{1}{2}\\ln\\left(\\frac{1-\\epsilon_t}{\\epsilon_t}\\right)} \\\\ \n",
    "&= \\displaystyle{\\frac{w_t}{z}}e^{\\ln\\left(\\sqrt{\\frac{\\epsilon_t}{1-\\epsilon_t}}\\right)} \\\\\n",
    "&= \\displaystyle{\\frac{w_t}{z}}\\sqrt{\\frac{\\epsilon_t}{1-\\epsilon_t}}\n",
    "\\end{align*}$\n",
    "\n",
    "##### Wrong\n",
    "\n",
    "$\\begin{align*}\n",
    "w_{t+1} &= \\displaystyle{\\frac{w_t}{z}e^{\\frac{1}{2}\\ln\\left(\\frac{1-\\epsilon_t}{\\epsilon_t}\\right)}} \\\\\n",
    "&= \\displaystyle{\\frac{w_t}{z}}e^{\\ln\\left(\\sqrt{\\frac{1-\\epsilon_t}{\\epsilon_t}}\\right)} \\\\\n",
    "&= \\displaystyle{\\frac{w_t}{z}}\\sqrt{\\frac{1-\\epsilon_t}{\\epsilon_t}}\n",
    "\\end{align*}$\n",
    "\n",
    "The sum of all $w_{t+1}$ should add up to 1, which means $z$ as a normalizing factor should be the sum of all the $w_t$ and their scaling factors, i.e.:\n",
    "\n",
    "$\\begin{align*}\n",
    "z &= \\sqrt{\\frac{\\epsilon_t}{1-\\epsilon_t}}\\displaystyle{\\sum_{Correct}w_t} + \\sqrt{\\frac{1-\\epsilon_t}{\\epsilon_t}}\\displaystyle{\\sum_{Wrong}}w_t\\\\\n",
    "  &= \\sqrt{\\displaystyle{\\frac{\\epsilon_t}{1-\\epsilon_t}}}(1-\\epsilon_t) + \\sqrt{\\displaystyle{\\frac{1-\\epsilon_t}{\\epsilon_t}}}(\\epsilon_t)\\\\\n",
    "  &= \\sqrt{\\epsilon_t(1-\\epsilon_t)} + \\sqrt{\\epsilon_t(1-\\epsilon_t)}\\\\\n",
    "  &= 2\\sqrt{\\epsilon_t(1-\\epsilon_t)}\n",
    "\\end{align*}$\n",
    "\n",
    "Substituting $z$ into the earlier formulas:\n",
    "\n",
    "##### Correct\n",
    "\n",
    "$\\begin{align*}\n",
    " w_{t+1} &= \\displaystyle{\\frac{w_t}{2}}\\left(\\frac{1}{\\sqrt{\\epsilon_t(1-\\epsilon_t)}}\\right)\\left(\\sqrt{\\frac{\\epsilon_t}{1-\\epsilon_t}}\\right)\\\\\n",
    "         &= \\displaystyle{\\frac{w_t}{2}}\\left(\\frac{1}{1-\\epsilon_t}\\right)\n",
    "\\end{align*}$\n",
    "\n",
    "##### Wrong\n",
    "\n",
    "$\\begin{align*}\n",
    " w_{t+1} &= \\displaystyle{\\frac{w_t}{2}}\\left(\\frac{1}{\\sqrt{\\epsilon_t(1-\\epsilon_t)}}\\right)\\left(\\sqrt{\\frac{1-\\epsilon_t}{\\epsilon_t}}\\right)\\\\\n",
    "         &= \\displaystyle{\\frac{w_t}{2}}\\left(\\frac{1}{\\epsilon_t}\\right)\n",
    "\\end{align*}$\n",
    "\n",
    "Now, taking the sum of the $w_t$ gives\n",
    "\n",
    "##### Correct\n",
    "\n",
    "$\\begin{align*}\n",
    " w_{t+1} &= \\displaystyle{\\frac{1}{2(1-\\epsilon_t)}}\\displaystyle{\\sum_{Correct}w_t}\\\\\n",
    "         &= \\displaystyle{\\frac{1}{2(1-\\epsilon_t)}}(1-\\epsilon_t)\\\\\n",
    "         &= \\displaystyle{\\frac{1}{2}}\n",
    "\\end{align*}$\n",
    "\n",
    "##### Wrong\n",
    "\n",
    "$\\begin{align*}\n",
    " w_{t+1} &= \\displaystyle{\\frac{1}{2(\\epsilon_t)}}\\displaystyle{\\sum_{Correct}w_t}\\\\\n",
    "         &= \\displaystyle{\\frac{1}{2(\\epsilon_t)}}(\\epsilon_t)\\\\\n",
    "         &= \\displaystyle{\\frac{1}{2}}\n",
    "\\end{align*}$\n",
    "\n",
    "What this all means is that when calculating the new weights, we only have to rescale the weights of the data points with correct and wrong predictions so that they both add up to $\\frac{1}{2}$ each.\n",
    "\n",
    "#### Example\n",
    "For example, when starting out, all weights are $\\frac{1}{N}$. But if $75\\%$ of the predictions are correct, then the new weights for the correctly predicted data points have to scale in such a way that they only add up to $\\frac{1}{2}$. This means that the sum $\\frac{3}{4}$ has to shrink to $\\frac{1}{2}$ and the new weights are $\\frac{2}{3}$ of the old weights (because $\\frac{3}{4}\\cdot\\frac{2}{3}=\\frac{1}{2}$). The set of new weights for the correctly predicted data points are $\\frac{1}{N}\\cdot\\frac{2}{3}=\\frac{2}{3N}$ for each data point. For the $25\\%$ of the data on which the classifier predicted incorrectly, the new weights are twice as much as the old weights ($\\frac{1}{4}\\cdot 2=\\frac{1}{2}$), which means the set of new weights for the incorrectly predicted data points is $\\frac{1}{N}\\cdot 2=\\frac{2}{N}$.\n",
    "\n",
    "Note that nothing about the $75-25$ split says either of these should be the correctly classified half. If the error rate were $75\\%$ instead, the same scaling factor of $\\frac{2}{3}$ would be used on the misclassified data points, meaning they become less important in the next round of training and the correctly classified $25\\%$ of the data would be boosted instead, becoming twice as important. **This is why we don't want $\\epsilon>0.5$, and as a consequence, we don't want $\\alpha < 0$.**\n",
    "\n",
    "#### Basic Algorithm for Binary AdaBoost\n",
    "\n",
    "1. Start by giving all data points weights $\\frac{1}{N}$, where $N$ is the number of data points\n",
    "2. For $t=1,\\dots,T$:\n",
    "    1. Train a weak learner $h_t(x)$ on the weighted data\n",
    "        - Use a metric that allows sample weights (e.g. the weighted Gini impurity) to prioritize more \"important\" data points when determining the decision boundary\n",
    "        - Aim for a learner with low error $\\epsilon_t$ (ideally $\\epsilon_t<0.5$).\n",
    "    2. Set $\\alpha_t = \\displaystyle{\\frac{1}{2}}\\ln\\left(\\frac{1-\\epsilon_t}{\\epsilon_t}\\right)$\n",
    "    3. Scale the weights of all data points such that $\\displaystyle{\\sum_{Wrong}} w_{t+1} = \\frac{1}{2}$ and $\\displaystyle{\\sum_{Correct}} w_{t+1} = \\frac{1}{2}$\n",
    "3. Output final hypothesis $H(x) = sign\\left(\\displaystyle{\\sum_{t=1}^T\\alpha_th_t(x)}\\right)$\n",
    "\n",
    "### Multi-Class AdaBoost\n",
    "\n",
    "The algorithm described above works in the binary case, but does not work well in the one-versus-rest multi-class case. Stagewise Additive Modeling using a Multi-class Exponential loss function (SAMME) is an algorithm that tweaks the basic AdaBoost algorithm a little bit to make it work in the multi-class case.\n",
    "\n",
    "#### Basic Algorithm for Multi-Class AdaBoost (SAMME)\n",
    "1. Start by giving all data points weights $\\frac{1}{N}$, where $N$ is the number of data points\n",
    "2. For $t=1,\\dots,T$:\n",
    "    1. Train a weak learner $h_t(x)$ on the weighted data\n",
    "        - Use a metric that allows sample weights (e.g. the weighted Gini impurity) to prioritize more \"important\" data points when determining the decision boundary\n",
    "        - Aim for a learner with low error $\\epsilon_t$ (ideally $\\epsilon_t<1-\\frac{1}{K}$, where $K$ is the number of classes).\n",
    "    2. Set $\\alpha_t = \\ln\\left(\\frac{1-\\epsilon_t}{\\epsilon_t}\\right)+\\ln(K-1)$.\n",
    "    3. Set $w_{t+1}=\\displaystyle{\\frac{w_te^{\\alpha_tc}}{z}}$, where $c=\\begin{cases}1&\\text{if }h(x)=y(x)\\\\0&\\text{if }h(x)\\neq y(x)\\end{cases}$ and $z$ is the normalizing factor\n",
    "3. Output final hypothesis $H(x) = \\displaystyle{\\text{arg}\\max_{k\\in K}}\\left(\\displaystyle{\\sum_{t=1}^T\\alpha_th_t(x)}\\right)$, where $h(x)\\in \\{0,1\\}^K$ and $h(x)_k=\\begin{cases}1&\\text{if }k=\\text{predicted class}\\\\0&\\text{otherwise}\\end{cases}$\n",
    "\n",
    "There are a few minor differences between SAMME and the binary AdaBoost:\n",
    "- The definition of $\\alpha$\n",
    "    - Since $\\alpha$ is no longer multiplied by $\\frac{1}{2}$, the sample-weights-summing-to-one-half-trick can no longer be used\n",
    "- In the multi-class case, we need $\\epsilon < 1 - \\frac{1}{K}$ to get a positive $\\alpha$.\n",
    "    - In the binary AdaBoost, we need $\\epsilon < 0.5$ to get a positive $\\alpha$. The generalized expression above still applies.\n",
    "    - The multi-class case still requires that a weak classifier be better than random guessing.\n",
    "        - Example: with three classes, randomly guessing means getting $\\frac{2}{3}$ of the guesses wrong, so the error rates must always be less than $\\frac{2}{3}$.\n",
    "- Obtaining the new sample weights\n",
    "    - When $h(x)\\neq y(x)$, the numerator remains unchanged because $c=0$ and $e^c=1$. The normalizing factor will still change these weights however as the distribution changes.\n",
    "    - In the binary algorithm, the $\\alpha$ is either multiplied by $-1$ or $+1$ depending on whether the prediction is correct or not.\n",
    "- The definition of $H(x)$\n",
    "    - Since the sign trick can no longer be used in a multi-class problem, the output of each $h_t(x)$ is simply a $K$-dimensional one-hot vector. These one-hot vectors are multiplied by the $\\alpha$s and summed over all $T$. Whichever class maximizes this $K$-dimensional vector is the final prediction.\n",
    "        - Example:\n",
    "            - given: $K=3$, $T=2$, and a single data point $x$, if:\n",
    "                - $\\alpha_1h_1(x) = 0.8\\begin{bmatrix}0&1&0\\end{bmatrix}$\n",
    "                - $\\alpha_2h_2(x) = 0.7\\begin{bmatrix}1&0&0\\end{bmatrix}$\n",
    "                - $H(x) = \\begin{bmatrix}0&0.8&0\\end{bmatrix} + \\begin{bmatrix}0.7&0&0\\end{bmatrix} = \\begin{bmatrix}0.7&0.8&0\\end{bmatrix}$\n",
    "                - then $\\text{arg}\\max_kH(x) = 2$, so the predicted class is $2$\n",
    "\n",
    "However, despite the minor differences, setting $K=2$ when running SAMME is pretty much equivalent to running the binary AdaBoost.\n",
    "\n",
    "## Gradient Boosting\n",
    "\n",
    "AdaBoost uses a specific loss function (exponential loss), but boosting could be generalized to any differentiable loss function. This is where gradient boosting comes in. Gradient boosting uses the gradients to determine where errors were made, whereas AdaBoost uses the sample weights. The generalized gradient boosting algorithm is usually defined as such:\n",
    "\n",
    "Given $N$ data points ($X$,$Y$) and a differentiable loss function $L(y,\\gamma)$,\n",
    "\n",
    "1. Set $f_0(x)=\\displaystyle{\\text{arg}\\min_\\gamma\\sum_{i=1}^NL(y_i,\\gamma)}$.\n",
    "    - Initialize model with a constant (single prediction).\n",
    "    - This is the **population minimizer**, the constant that minimizes the loss function across *all* data points because they are yet to be split.\n",
    "2. For $m=1,\\dots,M$:\n",
    "    1. Compute $r_{im} = \\displaystyle{-\\left[\\frac{\\partial L(y_i,f_{m-1}(x_i))}{\\partial f_{m-1}(x_i)}\\right]} \\text{ for }i=1,\\dots,N$.\n",
    "        - Compute pseudo-residuals (negative gradients). This is $-1\\times$ the first derivative of loss function.\n",
    "    2. Fit base learner on the pseudo-residuals $r_{im}$ with terminal regions $R_{jm}$, $j=1,\\dots,J_m$.\n",
    "        - Instead of fitting the tree on the $y_i$s, fit them on the error of the predictions thus far.\n",
    "        - Hastie et al. suggest that the number of leaves should be equal for all trees ($J_m=J$ $\\forall m$) to prevent the problem of having large trees at the beginning and smaller ones towards the end.\n",
    "            - When $J = 2$ (decision stumps), the algorithm only considers main effects.\n",
    "            - When $J = 3$, the algorithm allows interaction between two variables.\n",
    "            - $J$ then controls how many variables are allowed to interact.\n",
    "            - In practice, $J=2$ would be insufficient in most cases. But it is also unlikely that $J>10$ would be required. Hastie et al. suggest that $4\\leq J\\leq 8$ would suffice for most cases.\n",
    "    3. Compute $\\gamma_{jm} = \\displaystyle{\\text{arg}\\min_\\gamma \\sum_{x_i \\in R_{jm}}L(y_i,f_{m-1}(x_i)+\\gamma)}$, $\\text{ for }j=1,\\dots,J_m$.\n",
    "        - For each region $R_j$, calculate the constant $\\gamma$ that minimizes the loss function across all the data points in said region when added to their previous predictions.\n",
    "        - This is the population minimizer in each region.\n",
    "    4. Update $f_m(x) = f_{m-1}(x)+\\nu\\gamma_{jm}$, where $0<\\nu\\leq 1$ and $j$ is the index of the leaf node $x$ ends up in.\n",
    "        - The new prediction is the previous prediction plus the \"step\" $\\gamma$ scaled by a learning rate $\\nu$.\n",
    "            - The use of learning rates is a regularization method called shrinkage.\n",
    "            - In practice, gradient boosting models with small learning rates ($\\nu < 0.1$) have been found to generalize better than models without shrinkage ($\\nu=1$).\n",
    "        - This is equivalent to gradient descent. By fitting the tree on the negative gradients and setting the $\\gamma$ in each leaf to the value that minimizes the loss function across all the data points in that region, the direction of steepest descent in each region (that could be known given the data) is captured.\n",
    "3. Output $f_M(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Gradient Boosting](img/boosting/gradient_boosting.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image above shows what gradient boosting does at a very high level. The predictions by the regression trees are just added across stages to get the final prediction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient Boosting in Action: Regression with (Half) Squared Error Loss\n",
    "\n",
    "The above is easiest to explain with a simple loss function that has an easy-to-deal-with first derivative.\n",
    "\n",
    "Let $L(y,\\gamma) = \\displaystyle{\\frac{1}{2}}(y-\\gamma)^2$, the half of the squared error loss. The first derivative $\\displaystyle{\\frac{\\partial L(y,\\gamma)}{\\partial\\gamma}} = \\gamma-y$ and the negative gradient $-\\displaystyle{\\frac{\\partial L(y,\\gamma)}{\\partial\\gamma}} = y-\\gamma$. This is simply the residual. Plugging in this loss function into the generalized algorithm, we get:\n",
    "\n",
    "1. Set $f_0 = \\displaystyle{\\frac{1}{N}\\sum_{i=1}^N}x_i$.\n",
    "    - This is simply the mean of all data points (because the population minimizer of the squared loss is the mean).\n",
    "2. For $m=1,\\dots,M$:\n",
    "    1. Compute $r_{im}=y_i-\\gamma_{im-1}$.\n",
    "        - This is the previous prediction minus the true value (the negative gradient).\n",
    "    2. Fit a regression tree on these negative gradients.\n",
    "    3. Compute $\\gamma_{jm}=\\displaystyle{\\frac{1}{|R_{jm}|}\\sum_{r_{im}\\in R_{jm}}r_{im}} \\text{ for }j=1,\\dots,J_m$.\n",
    "        - For each leaf node, get the average negative gradient. This is the minimizer of the loss function in the region.\n",
    "    4. Update $f_m(x) = f_{m-1}(x)+\\nu\\gamma_{jm}$, where $0<\\nu\\leq 1$ and $j$ is the index of the leaf node $x$ ends up in.\n",
    "3. Output $f_M(x)$.\n",
    "\n",
    "### Other Loss Functions\n",
    "\n",
    "For regression problems, half squared error loss is easy to calculate but not robust to outliers. Because it pays too much attention to outliers, the performance of the model is affected negatively. In practice, other loss functions like absolute loss and Huber loss are preferred.\n",
    "\n",
    "- Absolute loss\n",
    "    - $L(y,\\gamma)=|y-\\gamma|$\n",
    "    - $-\\displaystyle{\\frac{\\partial L(y,\\gamma)}{\\partial\\gamma}}=\\text{sign}(y-\\gamma)$\n",
    "    - The population minimizer is the median. This makes procedure (1) and (2C) easy.\n",
    "- Huber Loss\n",
    "    - $L(y,\\gamma)\\begin{cases}\\frac{1}{2}(y-\\gamma)^2 & \\text{if }|y-\\gamma|\\leq\\delta\\\\\n",
    "\\delta|y-\\gamma|-\\frac{1}{2}\\delta^2 & \\text{otherwise}\\end{cases}$\n",
    "    - $-\\displaystyle{\\frac{\\partial L(y,\\gamma)}{\\partial\\gamma}}\\begin{cases}y-\\gamma&\\text{if }|y-\\gamma|\\leq\\delta\\\\\\text{sign}(y-\\gamma)&\\text{otherwise}\\end{cases}$\n",
    "    - There is no closed form for the population minimizer, so it is obtained iteratively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEdCAYAAADjFntmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3Xl4TNf/wPH3SSQSQuw7iZ1IIogl1rTWWmqJbt+qrbVUN7rpooJqi1a1fkWr1VJV1dqqFKW2IJaESAixxk42JCF7zu+PO9Ige2bmTpLzep55npm5d+79TJb7uffccz5HSClRFEVRFAArvQNQFEVRLIdKCoqiKEoGlRQURVGUDCopKIqiKBlUUlAURVEyqKSgKIqiZFBJQSkShBDPCyH+0TuO+4QQ9kKIv4QQd4QQf+gdj6URQmwWQozQOw4l/4Qap1CyCCH+B7wJNAPigCDgEynlXl0DK2KEEC8ArwEdpZSpWSyfBjSSUg7TITYJ3APu/3OnSikrmHB/09DpuyrGp64UShAhxJvAV8CnQHWgHrAQGKhnXLkRQpTSO4YsOAGns0oIFqKllNLB8DBZQlCKISmlepSAB+AIxANP5bBOabSkcc3w+AoobVjmDVwB3gUigOvAIKAvcBqIAT7ItK1pwGpgFdoVyRG0A9X95e8B5wzLQoHBmZaNBPYB8wzbnWl4b69huTAsiwDuAMGAa6bv+TMQCVwEpgBWmba7F/gCuAVcAJ7I4efRHNgF3AZOAE8a3p8OJAMphp/pi1l8dhrwS362a1jW1/DziAOuAm8b3q8CbDR8Jgbwu/+9sti+RDtzf/j9jJ9hVusCS4EFwCbD/g8CDTOt2wLYZtj/TeADoM9DP4tjhnV3AS8ZnlsZfg8XDb+znwFHwzJnQwwjgEtAFPBhpn22AwKAWMM+v9T7f6m4P3QPQD3M9IvW/nlTgVI5rDMDOABUA6oC+4GPDcu8DZ+fCtgAYwwH3l+BcoYDRiLQwLD+NMOBYqhh/bcNB2Ebw/KngFqGA8YzwF2gpmHZSMO+XgNKAfY8mBR6A4FABbQE0TzTZ38G/jTE5IyWsF7MtN0UQ+zWwMtoyU9k8bOwAc4aDny2wOOGA2XTTN8vy4N+TsvzsN3rQBfD84pAa8Pzz4BvDZ+3AbpkFbdh3cIkhRi0A3EpYAXwm2FZOUNsbwF2htfts/uuPJgURhu+cwPAAVgLLDcsczbE8L3h99wSSAKaG5b7Ay8YnjsAHfT+XyruD9V8VHJUBqJkzs0dzwMzpJQRUspItDPiFzItT0G7/5AC/IZ29vq1lDJOSnkC7azXPdP6gVLK1Yb1v0Q7mHQAkFL+IaW8JqVMl1KuAs6gHYzuuyal/D8pZaqUMuGhOFPQDkrN0A6MJ6WU14UQ1mgJ5n1DTOHA3Ie+w0Up5fdSyjRgGVATrSntYR3QDkKzpJTJUsodaGfqz+Xw88uL3LabArgIIcpLKW9JKY9ker8m4CSlTJFS+knDkTIbR4QQtw2P+fmIb62U8pDh72QF4GF4vz9wQ0o5V0qZaPj5HszjNp9HO8M/L6WMB94Hnn2oWXC6lDJBSnkMOIaWHED73o2EEFWklPFSygP5+C5KAaikUHJEA1VyaZ+vhXaJf99Fw3sZ2zAcTAHuH6hvZlqegHbAu+/y/SdSynS05qdaAEKI4UKIoPsHLsAVLck88tmHGQ6k36A1ddwUQiwWQpQ3fN42i+9QO9PrG5m2c8/wNHPM99UCLhvizm5bBZHbdn3QmpAuCiF2CyG8DO9/jna2/Y8Q4rwQ4r1c9tNaSlnB8Hg9H/HdyPT8Hv/9bOqiNfcVRFZ/V6V4MBlnt98XgSbAKSHEYSFE/wLGoOSRSgolhz9a886gHNa5hnYD9b56hvcKqu79J0IIK6AOcE0I4YTWXPAqUFlqN0KPozUF3Zdjtzgp5XwpZRu0ZqsmwDto7dEpWXyHqwWI/RpQ1xB3YbeV5+1KKQ9LKQeiNeGtB343vB8npXxLStkAGAC8KYTons993wXK3H8hhKiRj89eBhpmsyy3LoxZ/V2l8uAJRdYblvKMlPI5tJ/HbGC1EKJs7uEqBaWSQgkhpbyDdj9ggRBikBCijBDCRgjxhBBijmG1lcAUIURVIUQVw/q/FGK3bYQQQwxXJxPR2ooPAGXRDiSRAEKIUWhXCnkihGgrhGgvhLBBO9AlAmmGq5jfgU+EEOUMyefNAn6Hg4Ztv2v4OXmjHYx/y8c2rIQQdpkepXParhDC1jAew9HQ5BYLpBm+c38hRCMhhMj0flqWe83eMaCFEMJDCGGHdi8grzYCNYQQE4UQpQ0/3/aGZTcB54cSXWYrgUlCiPpCCAe03m+rcmnKBEAIMUwIUdVwZXXb8HZ+v7eSDyoplCBSyi/RDpJT0A7Il9HO1tcbVpmJ1tMjGAhB6zE0sxC7/BOtjf8WWrv+EEN7eChaW78/2gHFDa23UV6VR7vSuIXWFBGN1qMItJvTd4HzaD2NfgV+zG/gUspk4EngCbQrkIXAcCnlqXxs5jm0JrX7j3N52O4LQLgQIhYYD9zv+98Y2I7Ww8cfWCil3JXP73QarTPBdrR7OHkemyKljAN6oiWwG4bPP2ZYfH/wXrQQ4kgWH/8RWA7sQetskIj2e8qLPsAJIUQ88DXwrJQyMa9xK/mnBq8pJqEGNClK0aSuFBRFUZQMKikoiqIoGVTzkaIoipJBXSkoiqIoGVRSUBRFUTJYYvXJHFWpUkU6OzvrHYaiKEqREhgYGCWlrJrbekUuKTg7OxMQEKB3GIqiKEWKEOJi7mup5iNFURQlE5UUFEVRlAwqKSiKoigZitw9BaV4SElJ4cqVKyQmqjI2lszOzo46depgY2OjdyiKmZg8KRgmPgkArkop+z+0rDTaTFlt0IqaPWOYGEUp5q5cuUK5cuVwdnZGK/ypWBopJdHR0Vy5coX69evrHY5iJuZoPnoDOJnNsheBW1LKRmhz7s42QzyKBUhMTKRy5coqIVgwIQSVK1dWV3MljEmTghCiDtAP+CGbVQaiTYkI2iTv3YU6SpQY6ldt+dTvqOQx9ZXCV8C7QHo2y2tjmHbRMOHGHbS5hI3v3E5Y1BkSbplk80rRtG7dOoQQnDr13zQJu3bton//ws/6OHLkSFavXp3jOrt27WL//v352u6uXbtwdHTEw8Mj47F9+/bChKpYOilhaX8IyPfUIPlmsqRgmEs1QkoZmNNqWbz3SIU+IcRYIUSAECIgMjKyYAGVqQQ3QyD494J9XimWVq5cSefOnfntt/xMqGY8BUkKAF26dCEoKCjj0aNHjweWSylJT3/wXCwtLW8TluV1PcWMLh+EcD+wMn3fIFNeKXQCnhRChKNNYfi4EOLhaRGvYJjH1zBloyMQ8/CGpJSLpZSeUkrPqlVzHaWdtZotoaYHBC7Tsq5S4sXHx7Nv3z6WLFnySFKIjY1l8ODBuLi4MH78eNLT00lLS2PkyJG4urri5ubGvHnzAAgKCqJDhw64u7szePBgbt169GrU2dmZqKgoAAICAvD29iY8PJxvv/2WefPm4eHhgZ+fH5GRkfj4+NC2bVvatm3Lvn15n5AuPDyc5s2bM2HCBFq3bs3ly5dxcHBg6tSptG/fHn9/f/79919atWqFm5sbo0ePJikpKSO+GTNm0LlzZ/74449c9qSYXeAysHWAFkNMviuTJQUp5ftSyjpSSmfgWWBHFrNwbQBGGJ4PNaxjuiN2mxEQcQKu5nTxopQU69evp0+fPjRp0oRKlSpx5Mh/M0keOnSIuXPnEhISwrlz51i7di1BQUFcvXqV48ePExISwqhRowAYPnw4s2fPJjg4GDc3N6ZPn56n/Ts7OzN+/HgmTZpEUFAQXbp04Y033mDSpEkcPnyYNWvW8NJLL2X5WT8/vweaj86dOwdAWFgYw4cP5+jRozg5OXH37l1cXV05ePAgnp6ejBw5klWrVhESEkJqaiqLFi3K2KadnR179+7l2WefLeiPVDGFhNtwYh24DYXSDibfndnHKQghZgABUsoNwBJguRDiLNoVgmn/Gl2HwtYPIXAp1PE06a6UvJv+1wlCr8UadZsutcrjO6BFjuusXLmSiRMnAvDss8+ycuVKWrduDUC7du1o0KABAM899xx79+6le/funD9/ntdee41+/frRq1cv7ty5w+3bt+nWrRsAI0aM4Kmnnipw3Nu3byc0NDTjdWxsLHFxcZQrV+6B9bp06cLGjRsfeC88PBwnJyc6dOiQ8Z61tTU+Pj6AljDq169PkyZNMmJdsGBBxs/gmWeeKXDcigmF/AGpCdB6RO7rGoFZkoJhgvFdhudTM72fCBT8Pyi/7MqD6xA4vhb6fAaly+X+GaVYio6OZseOHRw/fhwhBGlpaQghmDNnDvBorxshBBUrVuTYsWNs3bqVBQsW8Pvvv2c0IeWmVKlSGW38OXXxTE9Px9/fH3t7+wJ9r7Jlyz7w2s7ODmtra0C7z5CfzyoWQEo4sgxquEGtVmbZZckb0dx6JBz9BUJWg+covaNRINczelNYvXo1w4cP57vvvst4r1u3buzduxfQmo8uXLiAk5MTq1atYuzYsURFRWFra4uPjw8NGzZk5MiRODo6UrFiRfz8/OjSpQvLly/PuGrIzNnZmcDAQJ544gnWrFmT8X65cuWIjf3vKqlXr1588803vPPOO4B2v8LDw8Mo37lZs2aEh4dz9uxZGjVqlG2sigW5dhRuhEDfL8BM3YNLXu2jOp5QzUXLvkqJtXLlSgYPHvzAez4+Pvz6668AeHl58d577+Hq6kr9+vUZPHgwV69exdvbGw8PD0aOHMlnn30GwLJly3jnnXdwd3cnKCiIqVOnPrI/X19f3njjDbp06ZJx5g4wYMAA1q1bl3Gjef78+QQEBODu7o6LiwvffvttlvE/fE8ht66voF01/PTTTzz11FO4ublhZWXF+PHj8/wzU3RwZBmUsgf3p822yyI3R7Onp6cs9HwKB76FLZNhnB/UdDdOYEq+nDx5kubNm+sdhpIH6nelk6R4mNsUmj8Jgxflvn4uhBCBUspcb6aWvCsF0LJuKTt1taAoiuU6sRaS46HNSLPutmQmhTKVwGUgBP8Byff0jkZRFOVRgcugajOo286suy2ZSQG07l1JdyB0vd6RKIqiPOjmCbgaoB2nzFx/quQmBaeOULmxlo0VRVEsSeAysLaFluYfSFhyk4IQ0Ho4XD4AEdlV9lYURTGzlAQI/k27wVymktl3X3KTAoDH/8DKBo78rHckiqIomtA/IfGOVpZHByU7KZStAs36wbGVkKImEilpPvnkE1q0aIG7uzseHh4cPHhQ75AArVyFq6trlu/b29s/MD7h55/VCU2xE7gMKjUA5y667L7kjWh+WJsR2s3mUxu1glNKieDv78/GjRs5cuQIpUuXJioqiuTkZJPuMy0t7YGBawXRsGFDgoKC8rWfvO43NTWVUqXUIUFXkafh0n7oMc3sN5jvK9lXCgD1vaGCk1YkTykxrl+/TpUqVShdujQAVapUoVatWgBs2bKFZs2a0blzZ15//fWMCXemTZvGF198kbENV1dXwsPDARg0aBBt2rShRYsWLF68OGOdh0tXBwYG0q1bN9q0aUPv3r25fv06AIGBgbRs2RIvLy8WLFiQ7+/z8H4eLoWdXXlvb29vPvjgA7p168bXX3+d/x+kYlxHlmlzJng8r1sIKilYWUHrF7QJLKLP6R2NYia9evXi8uXLNGnShAkTJrB7925AK1Y3ZswY/vrrL/z8/Lhx40aetvfjjz8SGBhIQEAA8+fPJzo6GuCB0tXt27fntddeY/Xq1QQGBjJ69Gg+/PBDAEaNGsX8+fPx9/fPcT/nzp17oPnIz8/vkf107twZeLAUdk7lvW/fvs3u3bt566238vdDVIwrNUlrym76BDhU0y0Mda0I4DEMdn6m3XDumbda+IoRbX5PK/plTDXc4IlZ2S52cHAgMDAQPz8/du7cyTPPPMOsWbPw8PCgfv36NG7cGIBhw4Y9cOafnfnz57Nu3ToALl++zJkzZ6hcufIjpauPHz9Oz549Aa1Zp2bNmo+U337hhRfYvHlzlvvJrvko837uu18KO7fy3qpktoU4tQnuRWtFO3WkkgJA+ZrQpI9WPfWxD6BUab0jUszA2toab29vvL29cXNzY9myZXh4eGQ7WX3m8tfwXwnsXbt2sX37dvz9/SlTpgze3t4Zyx4uXd2iRYtHrgZu376d7T7zKvN+7strKWxVMttCBPwIjvWg4WO6hqGSwn1tX4SwTRC6AdzNN8WDQo5n9KYSFhaGlZVVxhVBUFAQTk5ONGvWjAsXLnDu3DkaNmzIypUrMz7j7OycMbHNkSNHuHDhAqCdiVesWJEyZcpw6tQpDhw4kOU+mzZtSmRkJP7+/nh5eZGSksLp06dp0aIFjo6O7N27l86dO7NixQqjfte8lvdWdBQZpjVhd/cFq8J1Rigsk91TEELYCSEOCSGOCSFOCCEeaZcRQowUQkQKIYIMj6znHjSHBo9p3cAOf69bCIr5xMfHM2LECFxcXHB3dyc0NJRp06ZhZ2fH4sWL6devH507d8bJySnjMz4+PsTExODh4cGiRYsyZjDr06cPqampuLu789FHHz0w81lmtra2rF69msmTJ9OyZUs8PDzYv38/AD/99BOvvPIKXl5eOU6w8/A9hfnz5+fp++alvLeio8M/aCOYWw/XOxLTlc4W2vVwWSllvBDCBtgLvCGlPJBpnZGAp5Ty1bxu1yils7PjvwC2fqBKaptBUSnHvGvXLr744otHpr4sSYrK76rISoqHuc2gWV8Ykvv9q4LSvXS21MQbXtoYHpY9eYPH/7QJLQ7/oHckiqKUFMGrIDkO2urXUJKZSbukCiGshRBBQASwTUqZ1ZBRHyFEsBBitRCirinjyZV9RW0AW8gfkHBb11AUy+Dt7V2irxIUE5MSDi+BGu5Qp63e0QAmTgpSyjQppQdQB2gnhHh47P5fgLOU0h3YDmRZslQIMVYIESCECIiMjDRlyFq2Trmn9RdWFEUxpUv+EHFCO+7oNIL5YWYZvCalvA3sAvo89H60lDLJ8PJ7oE02n18spfSUUnpWrVrVpLFSy0PL2Id/0LK4oiiKqRz+AUo7gpvl9Hg0Ze+jqkKICobn9kAP4NRD69TM9PJJwDJqWLd9CaLPwvldekeiKEpxFXdT6wLf6nmwLaN3NBlMeaVQE9gphAgGDqPdU9gohJghhHjSsM7rhu6qx4DXgZEmjCfvXAZBmcrqhrOiKKZz5GdITwHPF/WO5AEmG7wmpQwGWmXx/tRMz98H3jdVDAVmYwetXoD98+HOFXCso3dEiqIUJ2mpEPiTNj6qSiO9o3mAKoiXHc/R2j0FVT1VURRjO70ZYq9CuzF6R/IIlRSyU9EJmvTWJrxINW2dfUU/Dg4OD7xeunQpr76a/VjK7CbAKYyEhAS6detGWloaoJXNGDx4MG3atMHNzY0ffihYM+bo0aOpVq3aI/EmJyfTtWtXUlNTCx27UkCHvofydaBxb70jeYRKCjlpOwbuRsDJDXpHohQTUsoHiuqBVnZ7yJAhGQXt1qxZQ7ly5QgMDCQkJITnny9Ybf2RI0eyZcuWR963tbWle/furFq1qkDbVQop8jRc2A2eo8Da8srPqaSQk4aPQ8X66oZzCfTwFcEXX3zBtGnTAG2GshEjRuDu7s7QoUO5d+9exnq//PIL7dq1w8PDg3HjxpGWlkZ4eDjNmzdnwoQJtG7dmsuXLz+wrxUrVjBw4MCM161bt2b37t14enri6+ubMRFQfnXt2pVKlbKe+H3QoEFGL7yn5FHAEm1ueAuoc5QVy0tTlsTKSque+s8UuHEcahi32UDRzD40m1Mxp3JfMR+aVWrG5HaTc10vISEBDw+PjNcxMTE8+eSTOXxCq7C6ZMkSOnXqxOjRo1m4cCFvv/02J0+eZNWqVezbtw8bGxsmTJjAihUr6Nq1K2FhYfz0008sXLjwgW0lJydz/vx5nJ2dAa3p6N133yU4OJiyZcvy+OOP4+HhweDBgzM+06VLF+Li4h6J64svvqBHjx65fmfQZo07fPhwntZVjCj5LgT9Ci0G6TqRTk5UUsiNx/OwY6Z2tTDgK72jUYzM3t7+gUlrli5dSm4FF+vWrUunTp0AbRKe+fPn8/bbb/Pvv/8SGBhI27ZauYKEhASqVatG165dcXJyyrJ6alRUFBUqVMh4/d1339G7d28cHR0B8PLyemT2t/uzrRWGtbU1tra2xMXFUa5cuUJvT8mj4N8hKdZi6hxlRSWF3JSpBK5DtV9mz+lg56h3RMVOXs7ozS27CXWARybEuf9aSsmIESP47LPPHlgeHh6e7UQ29vb2D2z76NGjjBgx4oHXAwYMeOAzxrhSAEhKSsLOzi7P6yuFJKV2clndDeq21zuabKl7CnnR9kVIMVz2KSVC9erViYiIIDo6mqSkpAeK4l26dClj9rSVK1dmzIncvXt3Vq9eTUREBKA1RV28eDHH/VSsWJG0tLSMxFCxYkWOHj0KwKZNm4iNjaVjx44PfMbPz4+goKBHHvlJCNHR0VStWhUbG5s8f0YppEsH4OZx7XhiIXWOsqKSQl7Ubq1l9oPfQnqa3tEoZmBjY8PUqVNp3749/fv3p1mzZhnLmjdvzrJly3B3dycmJoaXX34ZABcXF2bOnEmvXr1wd3enZ8+eXL9+Pdd99erVi7179wLwzjvvsG7dOlq2bMn333/P2rVrsbIq2L/pc889h5eXF2FhYdSpU4clS5ZkLNu5cyd9+/Yt0HaVAjqwAOwqgLuFz4ktpSxSjzZt2khdHF8rpW95KU9u1Gf/xUxoaKjeIViMI0eOyGHDhpl1n4MHD5anTp3K07rqd2UEMReknFZBym2+uoUABMg8HGPVlUJeNRugTartvzD3dRUlH1q1asVjjz2WMXjN1JKTkxk0aBBNmzY1y/4U4OBiEFbQbqzekeRKJYW8si4F7cfCxb1w/Zje0SjFzOjRozMGr5mara0tw4dbZh/5YikxVit+12IwlK+ldzS5UkkhP1q9ADZl4cAivSNRFKWoCFqhTbfZ4WW9I8kTlRTyw74CtBoGIash7kbu6yuKUrKlp2knkXU7QO0s5xCzOCop5Ff7cZCeqkpfGIFUM9tZPPU7KqSwv+H2RfCaoHckeaaSQn5VbghNn4CAHyElQe9oiiw7Ozuio6PVQceCSSmJjo5WA9wK48AirYNK0356R5JnJhvRLISwA/YApQ37WS2l9H1ondLAz2hzM0cDz0gpw00Vk9F0mKCdAQT/Dm1G5L6+8og6depw5coVIiMj9Q5FyYGdnR116qhJpgrkWhBc3Ae9PrHIaqjZMWWkScDjUsp4IYQNsFcIsVlKeSDTOi8Ct6SUjYQQzwKzAQsf2QE4d4YabtpZQOvhFj060VLZ2NhQv359vcNQFNM5sBBsHaD1C3pHki8maz4yjJeIN7y0MTwebisYCCwzPF8NdBcPF5axREJoVwuRJ+H8Tr2jURTF0sReh+NrtY4pRaxemknvKQghrIUQQUAEsE1KefChVWoDlwGklKnAHaByFtsZK4QIEEIEWExzg6sPlK2mBrMpivKowz9oHVLaj9M7knwzaVKQUqZJKT2AOkA7IcTDExJkdVXwyJ1HKeViKaWnlNKzatWqpgg1/0qV1uZXPbsNIsP0jkZRFEuRkqB1RGnWDyo10DuafDNL7yMp5W1gF9DnoUVXgLoAQohSgCMQY46YjKLNKLAurRXKUxRFAQheBQkxRWaw2sNMlhSEEFWFEBUMz+2BHsDD02ttAO533xkK7JBFqY+iQ1VwfxqCVsK9opPLFEUxESm1Dig13MGpk97RFIgprxRqAjuFEMHAYbR7ChuFEDOEEPfnO1wCVBZCnAXeBN4zYTym0WECpCZA4E96R6Ioit7O/QuRp8DrlSLbK9FkXVKllMFAqyzen5rpeSLwlKliMIvqLtDAGw59D16vQSlbvSNSFEUv/gvBoTq0GKJ3JAWmRjQbg9drEHcdQn7XOxJFUfRyPVi7Umg3tkifHKqkYAyNumvzru77GjLN66soSgmy72uwLQdtX9I7kkJRScEYhIDOEyHqNJzerHc0iqKYW8wFOLEWPEdq1ZSLMJUUjMVlEFRwgr3ztB4IiqKUHP7fgFUp6PCK3pEUmkoKxmJdCjq+BlcOw8X9ekejKIq5xEfC0V/A/RkoX1PvaApNJQVjajUMylTRrhYURSkZDn4LqUnQ6Q29IzEKlRSMycYeOozXSl/cOK53NIqimFpSHBz+Hpr3hyqN9Y7GKFRSMLa2L2nlcvd9pXckiqKYWuBSSLwDnSbpHYnRqKRgbPYVoc1IrWzurYt6R6MoiqmkJmuD1Zy7QJ2iMf9yXqikYAper4Cw0nokKIpSPIX8DnHXtO7oxYhKCqZQvha0fAaOLIe7UXpHoyiKsaWnw96vtBkYG3bXOxqjUknBVDq+AamJcPA7vSNRFMXYwv6G6DPQaWKRLXyXHZUUTKVqE22SjUOLISk+9/UVRSkapNS6nVd01gatFjMqKZhS50mQeBuOLMt9XUVRioaL++BqgDZY1dpkhaZ1o5KCKdXx1Hom7P9G66mgKErRt3celK0KHs/rHYlJmHLmtbpCiJ1CiJNCiBNCiEeG+wkhvIUQd4QQQYbH1Ky2VaR1mqj1UFBltRWl6LsRAme3Q/vx2mDVYsiU1z6pwFtSyiNCiHJAoBBim5Qy9KH1/KSU/U0Yh74addd6KPh9Ce7PFsvLTUUpMfZ8DqXLQ9sX9Y7EZEx2pSClvC6lPGJ4HgecBGqban8WSwjoNhlizsHxNXpHoyhKQd0MhdA/of04bZBqMWWWewpCCGe0qTkPZrHYSwhxTAixWQjRwhzxmF3TflDdVTvLSE/TOxpFUQpizxxtEp0OE/SOxKRMnhSEEA7AGmCilDL2ocVHACcpZUvg/4D12WxjrBAiQAgREBkZadqATcHKCrq9q/VrPr5W72gURcmviFNwYj20HwtlKukdjUmZNCkIIWzQEsIKKeUjR0MpZayUMt7w/G/ARghRJYv1FkspPaWUnlWrVjVlyKbTbABUc9HONtTVgqLRgOORAAAgAElEQVQULXvmgG1Z8HpV70hMzpS9jwSwBDgppfwym3VqGNZDCNHOEE+0qWLS1f2rhajTcGKd3tEoipJXkWHaFX67McX+KgFM2/uoE/ACECKECDK89wFQD0BK+S0wFHhZCJEKJADPSlmM57JsPhCqNofdc6DFYLCy1jsiRVFys3sO2JQBr9f0jsQsTJYUpJR7gRyLgkgpvwFKTinR+1cLq0dB6Hpw9dE7IkVRchJ5Wus12OkNKFtZ72jMQo1oNjeXQVC1Gez+XKu0qCiK5drzuTZIrWPJuEoAlRTMz8oKur4DkSfh5J96R6MoSnaizsLx1dpsimUf6f9SbKmkoIcWg6FKE62tUl0tKIpl2vM5lLKDjq/rHYlZqaSgBytr6PouRITCqb/0jkZRlIdFn9PqlXmOBoci2g2+gFRS0IvrEKjcWF0tKIol2vMFWJfWbjCXMCop6MXKWuuJdPM4hG3SOxpFUe6LOQ/Bq7Sidw7V9I7G7FRS0JOrD1RuBLtna7M5KYqivz1zwdqmxN1LuE8lBT1ZWWs9kW6EwMkNekejKEr0OTi2UruXUK663tHoQiUFvbk9pY1b2DET0lL1jkZRSrYdM7UeR50n6R2JblRS0JuVNTz+kVYT6dhKvaNRlJLrWhCcWAter5TIewn3qaRgCZr1g9qesOszSEnUOxpFKZn+nQH2laBj8a+EmpM8JQUhREMhRGnDc28hxOtCiAqmDa0EEQJ6TIPYq3D4B72jUZSS58IeOPcvdHkL7Bz1jkZXeb1SWAOkCSEaoZXDrg/8arKoSqL6XaBhd/CbC4l39I5GUUoOKWH7dChfWytpUcLlNSmkSylTgcHAV1LKSUBN04VVQnWfCgkxsL/kFI5VFN2d2gRXA8D7fbCx0zsa3eU1KaQIIZ4DRgAbDe/ZmCakEqyWB7QYAv4LID5C72gUpfhLT9PuJVRpAi2f0zsai5DXpDAK8AI+kVJeEELUB34xXVgl2ONTIDVRG2avKIppHfsNosK0/ztrU845VnTkKSlIKUOllK9LKVcKISoC5aSUs3L6jBCirhBipxDipBDihBDikSIiQjNfCHFWCBEshGhdwO9RfFRuCK1fgIAf4Va43tEoSvGVkqj1+KvVCpo/qXc0FiOvvY92CSHKCyEqAceAn4QQWc67nEkq8JaUsjnQAXhFCOHy0DpPAI0Nj7HAonxFX1x1m6yNX9j5md6RKErxFfAj3Lms9fwTOU4SWaLktfnIUUoZCwwBfpJStgF65PQBKeV1KeURw/M44CRQ+6HVBgI/S80BoIIQwmQ3sOOTisiI4fK1oP04rSjXzRN6R6MoxU9iLPh9AQ28tUcRkJiSRmqa6Ssq5zUplDIcrJ/mvxvNeSaEcAZaAQcfWlQbuJzp9RUeTRxG4Xcmki6zd7D2yBVkUSg+12kilC4P/36sdySKUvz4L4B70VqPvyLg4Plonvjajx/2XjD5vvKaFGYAW4FzUsrDQogGwJm8fFAI4YA2zmGi4WrjgcVZfOSRI7YQYqwQIkAIERAZGZnHkB9U09Ge+lXK8ubvxxi19DBXbycUaDtmU6YSdH4DTm+GSwf0jkZRio/4SPD/BlwGQu02ekeTo7jEFKasD+GZxQdITU/HrbbpB9YJU541CyFs0K4stkopH7kHIYT4DtglpVxpeB0GeEspr2e3TU9PTxkQEFCgeNLSJcv9w5mzNQwBvPdEM55v74SVlYW2JybfhfmtoGJ9GL1FtXsqijH8/a5WOeCVg1Clsd7RZGvnqQg+WBfCjdhERneqz1u9mlDGtuA9pIQQgVJKz9zWy+uN5jpCiHVCiAghxE0hxBohRJ1cPiPQRj+fzCohGGwAhht6IXUA7uSUEArL2kowslN9tk7sSmuninz05wmeXXyA85Hxptpl4diW1QbUXD4Aoev1jkZRir7I0xCwROvhZ6EJIeZuMhN/O8qopYdxKF2KNS935KP+LoVKCPmRpysFIcQ2tLIWyw1vDQOel1L2zOEznQE/IAS4f3fkA6AegJTyW0Pi+AboA9wDRkkpc7wMKMyVQmZSSlYHXuHjjaEkpqYzqUcTxnSpTylrC6sRmJ4G33XVboy9eliNuFSUwljxlNYc+9oRi5t7WUrJxuDrTNtwgjsJKUx4rBGvPNaQ0qWsjbL9vF4p5DUpBEkpPXJ7zxyMlRTui4hN5KM/j7P1xE1ca5dnjk9LXGqVN9r2jeL8bvj5Se2mWJe39I5GUYqmM9thhQ/0/Bg6WdasajdjE/lw3XG2n7yJex1HZvu407ymcY9DRm0+AqKEEMOEENaGxzAgunAhWoZq5e347gVPFj7fmht3knjym718sTWMxJQ0vUP7T4Nu0LQf+H0JcTf0jkZRip60VNj6gXZ/rv04vaPJIKVk1eFL9PhyN35nIvmgbzPWvtzR6AkhP/KaFEajdUe9AVwHhqKVvig2+rrVZPubXRnoUZtvdp6l33w/Ai/G6B3Wf3p9DKlJsEN1UVWUfAv8SStn0WsmlCqtdzQAXIq+x/M/HGTymhBcapZn68SujO3aUPcm7LyWubgkpXxSSllVSllNSjkIbSBbsVKhjC1zn27JstHtSExJZ+i3/kzbcIK7ljDorXJD7Qzn6ApthihFUfIm4Rbs/AScu2gTWuksLV3yg995en+1h+Ard/hksCsrx3TAuUpZvUMDCjfz2ptGi8LCdGtSla2TujK8gxPL/MPpNW8Pe04XbHyEUXV9Rxu/sPUDrQa8oii52z0HEm5Dn89079Z9+mYcPov2M3PTSbwaVmbbm10trlt8YZKC5XwLE3AoXYrpA135fZwXpW2sGP7jId7+4xh37qXoF5R9BXjsQ7i4D05u0C8ORSkqos7AocXQejjUcNMtjOTUdL7efoZ+8/24FHOPr5/1YMkIT2o62usWU3YKkxRKxKlqW+dK/P16FyZ4N2Td0av0mLebLcdNNpQid61HQDUX+OcjNZ+zouTmnylQyl4rja2TY5dv8+Q3e5m3/TRPuNZk2yTt3qWw0MGoOSYFIUScECI2i0ccUMtMMerOzsaad/s0489XOlHVoTTjfznCy78EEhGnw0HZuhT0/gRuX4SDqqisomTr7L9wegt0fRscqpl99wnJaXz690kGL9zH7Xsp/DDck/nPtaKyg2Xc6M6OSctcmIKxxynkV0paOt/7neer7Wewt7Hmo/4u+LTWIev/+iyE74XXAqFcdfPuW1EsXVoqfNsZUhPglUNm73Hkfy6a99YGczH6Hs+1q8f7fZtR3k7fySqNPU5BMbCxtmKCdyM2v9GFJtUdePuPYwz/8RCXY+6ZN5BeM7U/+J0zzbtfRSkKjiyFyJPaQDUzJoTYxBTeXxvCc99rRSx/HdOez4a46Z4Q8kMlhQJqWNWBVWO9mDGwBUcu3qL3V3tYuu8C6elmuvKq0gjajYMjy+H6MfPsU1GKgoTbsOMTcOoMzQeYbbf/nrxJry/3sOrwJcZ0qc+WN7rSsWEVs+3fWFRSKAQrK8FwL2e2TuqKp3Mlpv0VytPf+XM2wkwF9roZuqhuehvSTT/5hqIUCTs+hsTb0OdTs3RBjY5P4vWVR3lxWQCO9jasndCJD/u5YG9rnJpF5qaSghHUqViGZaPaMveplpyNjKfv134s2HmWFFPPkmRfUWtGunIIjv5s2n0pSlFwJRAOL4F2Y6FmS5PuSkrJn0FX6TlvD5uPX2dSjyb89VpnPOpWMOl+TU3daDayyLgkpm04waaQ67jULM+coe64mnJiDClhaX+4eRxeDbC4yo+KYjZpqfD9YxAfoVUUtjNd/aDrdxKYsu44/56KwKNuBeYMdadJ9XIm258xqBvNOqlarjQLnm/Nt8PaEBmfxMAF+5i95ZTpCuwJAf3mQnI8bCsaUwsqikkc/gFuBGsjl02UENLTJSsOXqTnl3vYdy6KKf2as+bljhafEPJDJQUT6eNag+2TuuHTujaLdp2j79d+HA43UYG9as2g4+tw7Fe44GeafSiKJYu9BjtmQsPu0GKwSXYRHnWX//1wgA/XHce9jiP/TOzGS10aYG1BJSqMQSUFE3IsY8OcoS355cX2JKel89S3/kz98zjxpiiw1/UdqFAPNr0FqcnG376iWLIt70NaMvT7wug3l1PT0lm85xy9v9rDiauxzBrixoqX2lOvchmj7sdSmCwpCCF+NEzfeTyb5d5CiDtCiCDDo9i2fXRuXIV/JnVldKf6LD9wkd7z9rArLMK4O7EtA33nauWB/f/PuNtWFEt2Zrs2XW3Xd6BSA6Nu+tSNWHwW7efTv0/RpXFVtr3ZjWfb1bPYEhXGYLIbzUKIrkA88LOU0jWL5d7A21LK/vnZrqXfaM5N4MVbTF4TzNmIeIa0rs1H/VyoWNbWeDtYNQzObIMJB6BSfeNtV1EsUUoCLOwAVjbw8j6jDVRLSk1jwc5zLNx5Fkd7G6YPbEE/t5pFOhnofqNZSrkHsJhZalLTU/G/5q93GLRxqsim1zvz+uON2BB0jZ7zdrMp+DpGS859ZoNVKdj8riqvrRR/fnPhVrjW2cJICeHopVsM+L+9zP/3DANa1mL7m93o715L94QQGh1KdILpJ7zU+56ClxDimBBisxCihSl39OfZPxm7bSzv7H7HLD/YnJQuZc2bvZqy4dXO1HS055VfjzBueSARsUYosOdYGx77AM78Ayf/Kvz2FMVSRZ6GvV+B+zPalLWFdC85lY83hjJk0X7iElP5aWRb5j3jYdwr+QJITE3kq8Cv+N+m/7EwaKHJ92fScQpCCGdgYzbNR+WBdCllvBCiL/C1lLJxNtsZC4wFqFevXpuLFy/mO5aUtBR+PP4j3wV/RxmbMkxuO5n+Dfrrnv1T09JZsvcCX247jW0pKz7q58JTnnUKF1daKiz2hnvR8OohKF18usspCqBdBS8boHVBfTWg0FVQ952N4r21wVyOSWBYh3pM7tOMchZQr+jIzSP47vclPDacwY0G85bnWziWLti4p7w2H+mWFLJYNxzwlFJG5bReYe8pnLt9jqn7pxIcGUyX2l2Y6jWVGmVrFHh7xnIh6i6T1wRz6EIMnRtV4dPBboXr3XD5MCzpCV6vaKW2FaU4Cf4d1o6Bfl9C2xcLvJk7CSl89vdJfjt8mfpVyjJriBvtG1Q2YqAFczflLl8FfsVvYb9R26E2vl6+eNXyKtQ2LT4pCCFqADellFII0Q5YDTjJXAIyxo3mtPQ0Vp5ayfyj87ESVkxqPYmnmj6FldC3NS09XfLroUvM2nyKtHTJ272bMrKjc8H7Qf81EY4sg5e2Q+02xg1WUfRyNwoWtIeKTvDidrAq2P/tPyduMGX9caLikxjTtQGTejTBzkb/ekV7r+5lhv8Mbty9wfPNn+e1Vq9Rxqbw3V91TwpCiJWAN1AFuAn4AjYAUspvhRCvAi8DqUAC8KaUcn9u2zVm76MrcVeY7j+dA9cP0KZ6G6Z3nI5TeSejbLswrt1O4MN1IewMi6RVvQrM8XGncUFGTCbchoVe2ujOsbvBxs74wSqKuf0+AsL+1v6mq7vk++NR8Un4bjjBpuDrNKtRjjlD3XGvo3+9otuJt/k84HM2nNtAA8cGTO84HY9qHkbbvu5JwVSM3SVVSsn6s+v5/PDnJKcnM8FjAsNdhlPKqpTR9lHQuP4Musb0v05wNymNVx9vxPhuDbEtlc+zojPbYYUPdJoIPaebJlhFMZfja2H1KOg+Fbq8la+PSilZH3SV6X+Fci8pjde7N2Jct4bYWOvbQiClZNvFbXxy8BNik2IZ7Taace7jsLU27g1ulRTyKeJeBJ8c+IQdl3fgUtmFGR1n0LRSU6PvJ7+i4pOY/lcofx27VvCzmj9fhaAV8OI2qJPr34SiWKb4SFjYHio4aX/L1nk/cbtquPreFRZJ63paAbtG1fTvgBF5L5JPDn7Cv5f+NflxRyWFApBS8s/Ff/j04KcmzdgFsS30JlPWhxAZl8SYLg2Y2KNJ3uu1J97RmpFsy8I4P9WMpBQ9UsLvw7U5l8f5afW+8uB+AbtZm08hgXd7N+UFr0LcpzOSjBaKgM9JTkvmFY9XeMHlBZO2UKikUAimbtsrqNhErafEykOXca5chlk+7nTIa0+Js//CL0O0wnm9PjZtoIpibMfXwOrR0GMadJ6Up4+cj4znvTUhHArXevR9NsSNupX0r1f08L3MaV7TcHZ0Nvl+VVIwAlP1Aiis/WejeG9tCJdi7vF8+3q890Qe+1RveB2OLofR/0DdtqYPVFGMIT5C621Uqb72t5tLs1FqWjrf+11g3vbT2JWyYkp/F55qU8ixP0aQudejQPBmmzfN2utRJQUjebi/8FSvqXSs1dFs+89OQnIac/8J48d9F6he3o5PBrvyeLPqOX8oMRYWdYRSdjDeD2zszROsohSUlP/V8xrvB1Vzbm8PvRbLu2uOcfxqLL1bVOfjga5UK69/c+m52+fw3e/LschjdK7dGV8vX7OPj1JJwcgCbwYybf80o4wsNKagy7eZvDqYsJtxDPSoxdT+LlR2yKEGzLmdsHwQeL2qBrUpli9kNax5EXrOgE5vZLtaYkoa3+w4y7e7z1GhjA3Tn3Slr1sN3a8OUtJT+DFEq6RQ1qYsk9tNpl/9frrEpZKCCSSlJbEoaBFLTyylol1FprSfQnen7rrEkllyajoLd51lwc6zlLOzYdqTLRjgnkNFx78mQuBSGL0V6rU3a6yKkmdxN7XeRpUbaX+rVll3rAi8GMO7q4M5F3kXn9Z1+Kh/cyqU0b9zyImoE0zdP5XTt07zhPMTTG43mcr2+o2WVknBhEKjQ/Hd78upmFP0curF++3fp4p9FV1jAgi7Ece7a4I5dvk2PZpXY+YgN2o4ZnHpnBQHCzuCtQ2M36vNxaAolkRK+O15OPev9jda5dGyaHeTUvl8axjL/MOp5WjPJ4Nd8W5auBpIxpCYmsjCoIUsC11GFbsqfNjhQx6v97jeYamkYGop6SksPb6URccWYV/KnsntJjOgwQDdL1fT0iU/7bvAF/+EYWNlxft9m/Ncu7qPxnV+N/z8JLQfD0/M1idYRcnOsd9g3TjoNRM6vvbIYr8zkby3JoRrdxIY3sGJd/o0w6G0vgNOAQJuBOC735dLcZfwaezDm55vUt7WNPNF55dKCmZy/s55fPf5EhQZRKfanfDt4EtNh5p6h8XF6Lu8tyYE//PReDWozCwfN5wql31wpb/fhUPfwXOroGkffQJVlIdFn4PvukINdxi58YFmozv3Uvh4UyirA6/QoGpZZvu409a5ko7BauKT45kXOI/fT/9OHYc6TOs4jfY1LatpViUFM0qX6fx26je+OvIVAsGkNpN4uunTuhfYk1Ky6vBlPtl0kpT0dN7q2ZTRnev/N3AnJRF+6AGxV7VZq8rX0jVeRSE1Savue/sSjN+nzQ9isOX4dT768wQxd5MZ360Brz3e2CIK2O25socZ/jOITIhkWPNhvNrqVexLWV7PPpUUdHA1/ioz/Gew/9p+WldrzbSO06jvqP+UmDfuJDJlfQjbT0bQso4jc4a2pGkNwxD/qDPwXTeo1QpGbMj2Zp6imMWW9+HAQnh2JTTrC0BEXCK+f55g8/EbtKhVnjlD3WlRS/+ef7cSbzH78Gw2nd9EowqNmN5xOu5V3fUOK1sqKehESsmGcxuYc3gOiamJvOzxMiNbjLSIAnsbg68zbcMJYhNTmODdiFcea6QV2Av6Fda/DI99CN3e1TVOpQQL2wIrn8m4zyWlZHXgFWZuOklCShoTezRmTJcGFlHAbkv4Fj47+BlxKXGMcRvDGLcx2FjrPylPTlRS0FlUQhSfHvyUbRe30bxSc2Z0mkGzSnmr12JKMXeTmfHXCdYHXaNJdQdm+7jTqm4FWDsWjq+GkZvASf/BeUoJE3sNFnXSmote+pcrcWm8vzYEvzNReDpVZPZQdxpWddA7Sm7evcnMgzPZdXkXrpVdmd5pOk0qNtE7rDxRScFCbL+4nZkHZnI76TajXUczruU4SlsbZ4Lxwthx6iYfrjvOzdhERneqz1vdamH/02Nam+74vVBG/5t3SgmRngbLnoRrR0kfs4vlZ22ZveUUAJP7NOOFDk5YWUABuzVn1jA3YC6p6am82upVhjUfhnURam7VPSkIIX4E+gMR2cy8JoCvgb7APWCklPJIbtstakkB4E7SHb4I+IL1Z9dT37E+MzrOsIgCe3GJKczecopfDlyiXqUy/J83tNzyFDTuBc+uAJ271yolxK7ZsOtTbj7+Fa+caErAxVt0bVKVTwe7Uqei/mNoLsdeZpr/NA7dOETbGm2Z5jWNeuXr6R1WvuU1KZiycW4pkFM/xyeAxobHWGCRCWPRlWNpRz7u9DHf9fiOpNQkhm8ezmcHP+Neyj1d4ypnZ8PMQW78NrYDVgIGrr3HxurjIGwTHPpe19iUEiJ8H3L3LMKq9aXL1hqciYhn7lMtWTaqre4JIS09jWUnljFkwxBtwKqXL0t6LSmSCSE/9Jyj+Ttgl5RypeF1GOAtpbye0zaL4pVCZvdS7jH/6Hx+PfkrNcvWxNfLl4619W/DT0hO46vtp/ne7xw/231JRxGM1ZgdUNNye1MoRdy9GFIWeBGRIOh1bybd3Ooz/UlXqpbTv3n1zK0z+O73JSQqBO863kzpMIXqZXMpOGnhLOFKITe1gcuZXl8xvFeslbEpw3vt3uPnJ36mdKnSjNs+jil7p3An6Y6ucdnbWvN+3+asf6Uz/1duEhFpDtz88X9Ex0TrGpdSPCUmp3Lm++HI+EgmW73J3GGdWPh8G90TQkpaCouCFvH0xqe5EneF2V1mM//x+UU+IeSHnkkhqwbrLC9bhBBjhRABQoiAyMhIE4dlHh7VPPhjwB+McRvDxvMbGbh+INsubtM7LNzrVGD56/3Y5/4ZVZKvcHT+c/x59DJFrUOCYrkOh8ew4ovXaXzLj801J7DgzVH0cdW/CkBIZAhPb3yahccW0tOpJ+sHradvg766l64xN9V8ZAFOxZxi6r6pnIw5SU+nnnzQ/gOLKLAX8c9cqu2fwdyUoZxoPJ6Zg1ypVcHyRmoqRUN8Uipztpzi6sG1LLGdS4Tzk1Qb8bPuHRoSUhNYcHQBy08up4p9FT7q8BHedb11jckUikLz0QZguNB0AO7klhCKq2aVmvFrv1+Z2Hoiuy/vZuD6gaw/u173s/NqPd8k3e1p3rJZjf25LfSat4cVBy+Snq6uGpT82RUWQe95e9h/cD8L7RaRVqMl1Z5frHtCOHT9ED4bfFgWugyfxj6sH7i+WCaE/DBll9SVgDdQBbgJ+AI2AFLKbw1dUr9B66F0Dxglpcz1EqA4XilkduHOBabtn8aRiCN41fTCt6MvtR10vNWSkgA/PUF65GneqzSP3y860L5+JWb5uFO/StncP6+UaLfuJvPxplDWHrmKexXJ71YfYpd2F8buAsc6usUVlxzHl4Ffsvr0auqWq8v0jtNpW6N4T1Gr+zgFUynuSQG0Anu/h/3OvMB5SCRvtH6D55o9p1+BvTtXYbE3srQD6z2XM/WfqySnpvNmzya82Lk+pXQuO6BYHiklm4/fYOqfx7l9L4WXuzozMWIK1uF7YMRf4OSlW2y7L+9mxoEZRCVEMdxlOBM8JlhkATtjU0mhGLgef53pB6az7+o+PKp6ML3jdBpUaKBPMJcOwNL+UL8LNwf8wpQNJ9kWehP3Oo7M9nGneU3LqBmv6C8iNpGP/jzO1hM3ca1dntk+7rQ4MRf2fQ39vwLPUbrEFZMYw6xDs9h8YTONKjTi404f41rlkdudxZZKCsWElJKN5zcy+/Bs7qXc4+WWLzPSdSQ2VjoU3wpcCn+9AR1fR/acwd8hN/DdoJ0JTvBuyCuPN6J0qaIz7F8xLiklfwReYebGUJJS05nUswkvda5PqdC12jzLnqOh/zxd4vr7wt/MOjSL+JR4xrmP40XXFy2+gJ2xqaRQzEQlRDHr0Cy2hm+lacWmzOg0A5fKLuYPZNNbcPgHGPI9uD+ttRlvDGXt0as0qqYV2GvjVNH8cSm6uhxzj/fXhrD3bBTtnCsxy8eNBlUd4FoQ/NgbarWG4X9CKfPOnXzj7g1mHpjJ7iu7ca/izvSO02lUsZFZY7AUKikUU/9e+peZB2ZyK/EWI1qM4OWWL2NXKot5mE0lLQV+HghXA2H0Fm0eBmBnWAQfrg3hemwiIzs6807vppSx1X96RMW00tIly/aH8/nWMKwEvNe3Oc+3q6cVsIuPhMXe2opjd4FDVbPFlS7TWX16NV8Gfkm6TOe1Vq/xv2b/K1IF7IxNJYVi7E7SHb4M/JK1Z9biXN6ZaR2n0aZ6G/MFEB8J3z8G6anw4jaoUBfQCuzN2RLG8gMXqVPRnllD3OncWP/xFoppnLkZx+Q1wRy5dBvvplX5ZLAbte+PY0m+p5083AiG0VuhlvkKQF6KvYTvfl8CbgbQvmZ7fL18qVuurtn2b6lUUigB/K/5M91/Olfjr/JM02eY1GYSZW3M1E30xnH4qS+Uq67902cqtX3oQgyT1wRzIeouT3vW4cN+Ljjal6z22+IsJS2db3ed4/92nKVMaWt8B7gwyKP2fyN/01Lgt+fhzD/w9DJwGWiWuFLTU/kl9Be+CfoGWytb3m77NoMbDS5xI5Kzo5JCCXEv5R7/d/T/WHFyBdXLVmdqh6l0qdPFPDsP3wvLh2hF84b/Cbb/JaTElDS+/vcMi/ecp3JZWz4e5ErvFjXME5diMiFX7vDO6mOcuhFHP/eaTH+yBVUcMtUrkhL+fBWCfoF+X0LbF80SV1hMGL77fTkRfYLH6j7GlA5TqFammln2XVSopFDCHIs8hu8+X87dOceABgN4t+27VLCrYPodh26AP0ZAo57aHAwP9eg4fvUO764OJvR6LP3cajLtyRa6Fz1T8i8xJY1520/z/Z7zVHEonX2S3z4d9n4J3d6Dx943eU76yNMAABwXSURBVFzJacksDl7MkpAllC9dng/af0Avp17q6iALKimUQA//g7zf/n16O/U2/T9IwI+wcRK0/B8MWvhI6YKUtHQW7znP19vPYG9rzdT+LgxpXVv94xYRB85H896aYMKj7/Fs27q837d51s2BBxbBlvegzSit66mJf7+6nQgVUSoplGBhMWFM3T+V0OhQHq/7OB92+ND0l9KG2bPoNBF6Ts9ylbMR8UxeE0zgxVt0a1KVT4dkujGpWJy4xBRmbT7FioOXqFtJ6zjQqVE2HQdCVmtjEZoPgKeWgQl7+ejaZFqEqaRQwqWmp7I8dDkLghaY56ablNoYhoAl0Psz8JqQ5Wrp6ZKf/cOZszUMAUx+ohnD2us/B6/yoJ2nIvhgXQg3YhMZ1bE+b/dukn0X43M7YMXTULcdDFsLNqbrIn3g+gGm7Z+W0bliYuuJONg6mGx/xYlKCgpg5u556Wnwx0g4uQF8loDb0GxXvRxzjw/WheB3Joq2zhWZ5eNOw6rqn1tvMYbBiOuOXqVxNQfmDHWnVb0cBiNeO6qVP6noDCM3gb1pmm9ik2OZGzCXtWfW4lTeiekdp5u3G3YxoJKCksGsA3lSEmHFUK1W0nMroXHPbFeVUrLmyFU+3hhKQkoaE3s0ZkyXBtioAntmJ6VkY/B1pm04wZ2EFF55rBETHmuYc9mSyNPw0xNgWwZG/wPlTTNRzo5LO5h5YCYxiTGMbDGS8S3Hm3fAZjGhkoLyiBt3b/DxgY/Zc2WPaYf8J96BZQMg4iQ8/TM0/f/27jw6qipb4PBvkxkIYUiCDIGQMMkQRlFAUAQBQZltgfcU7eah9sJpaSu0M9oNOD1t235OONPAU0CjhkZBERCZRBIIgQAJZCCQhJCQeTzvj1um82ISqpJUpYD9rZVFVerWqV03Re17zzl3n5vq3Dw9t4inv4hlw8HT9O1oFVDr1ymg8eNSNTpzvojH1x9kU5xV4PCFWRH0vuICBQ7PHIKPpgACd22AwMb/HFUv7fLsyGfp265vo7/O5UKTgqqRy4qDFZ6DT2ZCWjTMes+uC5g2HEjjyS9iOVdQwj3XhXHfDT3w9bp8yxI4mzGGNXuS+UtUHCVlFTwyvhd3jQy9cCn0tBjramVPH6sMdmCPRo/LbYpAXkI0Kag6ZRVlsWzXMjac2ECPNj1YMmJJ45cRLjpvdSWl7IUZb9c5xvCr7IISnvsqjrX7UggLasELMyMYGtr2gs9Tjkk6W8CidTHsOH6Wa8LasmxGBKH2LJqUug8+ng7eLWFeJLQLb9S40vLSWLJzCdtTtzMgaABLRixpunLxlxi3SAoiMhF4DfAA3jXGLKv2+J3Ai0Cq7Vd/N8a8W1ebmhQa15bkLTy38zkyCzOZ12ce9w68t3EXHCnOg3/eBkk7YOo/YOAcu562NT6DxesOcCqnkHnDrQJ7LXy0wF5DlVcY3v8xkZe+OYJXs2YsnnQls68KsW/2V/Ju6+zPr411htCma6PFVdPCUrN7zb6sC9g1tiZPCiLiAcQDNwIpwB5gjjHmUJVt7gSGGmMW2tuuJoXGV3Vpwi7+XXhmxDONuzRhSQGsngMJP8CUv8HgO+x6Wn5xGS9uPMKHP52gY4Aff53Rn+t6uq7S5qXmyOlcHl0bQ3RyNmN7B/P89H50CLDzAODkDlh5K7QMthJCIy6leSLnBE/veNp9lqC9RNmbFJw5zWMYcMwYk2CMKQFWA66pjKUc4u/tz9PDn2bF+BUYDL/f+HuW/LSE3JLcxnkB7+YwZzV0HwuR98Hud+x6WgsfT56Z0pdP7x6Oj1cz5r23m4f/N5rsgpLGiesyUVJWwaub4rn59W0kZxXw2uyBvDtvqP0JIeEH6wyhVUe4M6rREkJZRRkrDqxgZuRMjmYf5bmRz/HWjW9pQmhizkwKnYDkKvdTbL+rbqaIxIjIZyKi9W2b0LAOw1g7ZS139r2TtUfXMu2LaWxO2kyjnE16+cHsf0KvSRD1CPz0ht1PHRralqj7R7FwTHc+35/KuFe2EnUgreExXQaik7O55fXtvLrpKJP6d+Dbh0YzdaADJUaObYJ//u7f1yE00rTT2MxY5n49l1f3vcrozqOJnBbJtO7TtPSJG3Bm99GtwARjzHzb/duBYcaY+6ps0w7IM8YUi8g9wO+MMTfU0NYCYAFAly5dhpw8edIpMat/O5h5kKd3PE38uXiu73w9i69eTMeWHRvecFmJVQ4hLhKGL4QblzhUEiH2lFVgL/bUeSb2vYIlU/sS3ErnrFdXWFLOK98eYcX2RIL9ffnL9H6MvbK9Y43s+8iqaRV8Jdz+BbRo1+C4cktyef2X11l9eDWBfoEsGraI8aHjG9yuujB3GFMYDjxjjJlgu78YwBiztJbtPYAsY0ydE9R1TMF1yirKWBm3kjf2W0f19wy4h9v73N7wqYHlZbBxMex+G3reBDPfBR/7r2YuK6/gnW2J/PemeHw9m/HkzX2YNaSzHmXa7DieyaK1B0jKKmDu1V1YdFNvWvk68DerKIdNz8COv0HYGLj1gwZfqWyM4V8n/sULe17gbOFZ5vSew8JBC/H39m9Qu8p+7pAUPLEGmsdizS7aA8w1xsRW2aaDMSbNdns68Jgx5pq62tWk4HppeWks272M75K/o3vr7jx5zZMMbj+44Q3vfgc2PAbBfWDuaof7qo9n5LFobQx7TpxjVI9A/jq9PyFtmzc8rovU+aJSlkYdZtXuJLq2a86yGREMD3fw6L44D9YtgCNfw1XzYeJy8GjYrK+k80k8v/N5fkr7iT7t+vDUNU/RN1AvQnO1Jk8KtiAmAa9iTUl9zxjzFxFZAuw1xkSKyFJgClAGZAH3GmMO19WmJoWm833S9yzdvZS0/DSmd5/OQ0Meoo1vHXVx7HFsE3x6l23MYRV0dqyeTUWFYeWukyzbcBgDPDqhF3cMD73sCuxtjjvD4+sPkp5bxPxRYTw0rid+3g5O58xJhVW3wZlYmLgMrr67QTGVlJew4uAK3o15Fy8PL+4fdD+39bpNp5k2EbdICs6gSaFpFZQW8GbMm3wc+zHNvZqzIGIBc3rPwdvDu/6NpsdZg5l56TDtf6DfDIebSM0u5M/rDvBDfAZDurZh+cz+dA++9LsmzuYV8+yXh4iMPkXvK/xZPjOCASH16OpJ3Qer5kBJPtz6fp01qy7EGMPGkxt57efXSMlLYWLoRP501Z90JbQmpklBOdXRc0d5ae9L7Di1g04tO/Hg4AeZENqABX3yM2H1XEjeBWOegNGPOLxIizGG9b+ksuSrQxQUl/PAuB4sGH1pFtgzxhAZfYpnImPJKy5j4Zge3Ht9ON6e9XivsZ/D+nugRRDMXQPt+9Q7rn1n9vHy3peJyYyhR5sePDL0EUZ0HFHv9lTj0aSgXGJH6g5e/vll4s/F0z+wPw8Pfbj+JY1Li+DL+yFmDfSZBre8Vq8BzozcYp75MpavY9K4skMrXpx1aRXYS8sp5In1B9l8OJ2BIa15YVYEPdvX46yovBS2LINtL0Hnq6wpwy3rdzR/IucEr+57lc1Jmwn2C2bhoIVMCZ+iXUVuRJOCcpnyinK+TPiS1/e9TnphOjeE3MCDQx6kW0A3xxszBn58FTY/Z10sNeNt6Fq/I82Nsad54vODZOWX8F+jwnhw3MVdYK+iwrBqTxJLow5TXmF4ZEIv7hwRikd9xk+yEmDtf0HqXhj4nzD55XotjpNVlMWb0W/y6ZFP8fbw5q5+d3FHnzto7nX5Dvi7K00KyuUKywr5+NDHrDiwguLyYqZ2n8r8fvMJaVWPaxJT9lrXM2QnwahH4LrH6jULJqeglL9GxbFmbzJhgS1YNjOCYd0uvgJ7JzLzWbQuhp0JWYzs3o6l0yPo0q4eX7zGQPQqiPoTiAfc8mq9xnCyi7L5JO4TVsatpLCskJk9ZnLvwHsJ9KtluU7V5DQpqCaTWZjJ2zFvszZ+LeWmnEndJjE/Yj5hAQ5WuyzOhahHIfqfVvfGjLehbf0qZm4/msni9TEkZxVy+zVdeXRiL/wdmbvfRMrKK3jvx0Re/iYeb89mPDH5Sn43NKR+YzeF2dbFaLHroOtImP4WtHYsYWcWZvJR7EesPrKawrJCxnUZx32D7tNKphcBTQqqyWUUZPBB7Ad8Gv8pRWVFjA8dz4KIBfRs09Oxhg6uhS8fAlMOk16CAbMdHoQGKCgp46WN8by/I5EOrXz5y4z+jOnlvjNi4tLO89jaGGJScrixT3uen9aP9vW9evvkDuv6g/OnYMyf4dqHHLqS/Ez+GT6I/YDP4j+jpKKECaETWNB/gXMWaVJOoUlBuY2soiw+PvQxqw6vIr80nzEhY7g74m7HLmDKToJ1d1sluPvNhJterHfZhX1J53jssxiOpucxY1Annry5D21aNGBKbSMrLivnje+P84/vj9G6uRfPTunHpP5X1O/soLQItr4I21+B1l2ttbMduBYkNS+V9w68x/pj66kwFdwcdjPz+88nNCDU8VhUk9KkoNxOTnEOK+NW8kncJ+SW5DIoeBBzes9hXJdx9q38VlFufbl9v9QqizHmCRj6+3qNNRSXlfPGd8f4x5bjBPh58ezUvkzu36HJS2U0WsIyBo5ssMqJnDsBA/8DbloOPheepWSMYdfpXayKW8WWlC00k2ZM6z6NP/T7A539G69ktnItTQrKbeWV5LH26FrWHFlDcm4ygX6B3NrzVmb1nGXfBU7pcVZ5jMQfrBIZNy2HbqPrFUvVLprxfdrzXEO6aBqgUbu2Mo7AvxbB8e8gqLd1dXL4mAs+Lb80n8jjkaw+vJqEnATa+LRhZs+Z3NbrNq5ocUX9YlFuQ5OCcnsVpoIfU39k1eFVbE/djod4MK7rOOb0nsOg4EF1H7UbA3FfwjePW11LfabC+OehdReH42jUwdx6+PFYJovWWYPgdwzvyqMTe9OyPqvMFeXAluWw+y3wagFjFlv1iy5wFpaQk8Dqw6uJPB5Jfmk+fdv1Ze6Vc5kQOgEfD596vivlbjQpqItK0vkk1hxZw/pj68ktySU8IJzJYZOZFDap7kVXSgthx+uw7RXAWAOoIx+waik5KDEzn0VrY9iVmMWI8HYsm1HPaZ92yiks5a9fW9NluwW2YHl9p8tWVMD+lbD5WevK8MF3wNinoEXt00Ozi7L55uQ3fJ3wNfvS9+HVzIuJoROZ03sO/YP6N+BdKXelSUFdlApKC4hKjCLyeCS/pP8CwODgwUwOm8yE0AkE+NRyZXJ2Mnz7lDXd0r8DXHMvDLkTfB27krlRLxCrw8bY0zz5+UHO5pewYHQYD4ytx4V1ZSXWzKwdr0N6LIRcbXWldRxU4+ZFZUX8kPIDXyV8xfbU7ZRVlBEWEMYt4bcwvft02vk1fL0E5b40KaiLXkpuClGJUXyd8DUJOQl4NvNkVKdRTA6bzKhOo2q+avbEj/DDcmu8wdsfhsyzEoSDZbnTcgp5fP1BvmtoKYlqMnKLeSYylq8PpNGnQyteqE8JjqIc2Ps+7HoTctOscZVRD1uzsqp1eZWWl7L3zF6iEqPYdHITeaV5BPkFcVO3m7g57GZ6t+3d5IPryjU0KahLhjGGw1mH+SrhKzYkbiCjMAOvZl4MbT+U0Z1HM7rzaLq0qjaWcGo//PR3OLjO+qLsN9Na6a1DhEOv21hF56oX67t/bHfuvi7csWJ9OSmw83/g5w+hJBe6XQcj7rfWvq7yxZ5ZmMm2lG1sS93GjlM7yC/Np4VXC8Z1GcfksMkMu2KY1iS6DGlSUJek8opyfj7zM1tTtrI1dSuJOYkAhLYK5dpO1zK682iGth/67ymu2Umw803Y9yGU5EHY9TB4HvQYb/dqb5m28tRf1rM8dWp2IY+vP8CWI7+W9Y6ge7CdK82Vl0LiVohebXWNGWOVpRhxH3QYAFgD9ofOHrL2ScpWYs9a61gF+wUzqvMoRnUexciOI/H11GVLL2eaFNRlITk3mW0p29iaupU9aXsoqSjB18OXvoF9GRg0kEHBgxgQNIDWBvj5fdj1ltXl4ukL4WOhzxToOdGuaqzfHjrDE58fICO32K6FbOq9AFBpESR8D4ci4UgUFGVbXWGD74Br7qWgRTsOZh5kf8Z+9qfvJzojmvMl5xGEiKCIyrOnXm16adeQquQWSUFEJgKvYa289q4xZlm1x32Aj4AhwFngNmPMibra1KSgalNQWsDu07vZfXo3+9P3E3c2jjJTBlhnEgODBzIwMILeJSV0S/qZ5oc3QO4paOYFYdfBlVOg9+Q6Z+1YS17GsWp3MqHtmrO0liUvEzLyWLT2ALtPZNm3VGhJvrUK3aFIiN9odQ/5BFDacwInuw0n3r8t0Vlx/JL+C/Hn4ik35QB0b92dAUEDGNJ+CNd2urbhK+GpS1aTJwUR8cBao/lGIAVrjeY5xphDVbb5IxBhjLlHRGYD040xt9XVriYFZa+isiJiz8ayP906ot6fsZ/s4uzKxzu26Eg3n7aEFxcSln6c8JzThJaWERDQBekw0Oqe+fWnWqLYcSyTResOkJRVwNyru7Dopt608vWirLyCd7Yl8t+b4vH1bMaTN/dh1pDO//+IvTgXTh+AtGhIi4G0aPIzj3DSUzjesg2J7XtyvHkACWXnSc5NqUwAfp5+RARGMCB4AAODBhIRFFH7bCylqnGHpDAceMYYM8F2fzGAMWZplW022rb5SUQ8gdNAkKkjKE0Kqr6MMZw8f5Jj2cdIyEngePZxEnISSMxJpLi8uHI7b4SgckNgaTFB5eUElpcT5NmCIP8QWvl3xNu7Fb4+rWjmHcC2xBI2Hy/Cy7cd0666ks0HE0g/e4arOgsTe/nhZfIpKc6huPg8+fnpZJxPIqM4mwyPZmR4epDp6U2GpweF/Psj7ymehLQKITwgnG4B3QhvHU731t0Jbx2OZ7N6XNSmFPYnBWd+wjoByVXupwBX17aNMaZMRHKAdkCmE+NSlykRITQg9DfF3MoryjmVf4rEnERO5JwgszCTjMIMMvLSSMxLZXfROc5XFEN5KmSn/rZh21pCL6YAra2fFGD9iRqC8IXmfm0I8mlNYMsO9GnZiUC/QIKaBxHiH0JYQBhd/LvYVwtKKSdwZlKoaYSr+hmAPdsgIguABQBdujhexkCpung08yDEP4QQ/xBGd665hlJxeTGZhZnkleRRVJJHSWEWxYVZFBdlU1yUQ2FRDiczMujQpjX+Ldrg6xuAt09rfPza4OPXFh+fVjT3ak6QX5CuSqbcmjOTQgpQdQWPzsCpWrZJsXUfBQBZ1RsyxrwNvA1W95FTolWqDj4ePnWX21DqEuH4VTj22wP0EJFuIuINzAYiq20TCcyz3Z4FfFfXeIJSSinnctqZgm2MYCGwEWtK6nvGmFgRWQLsNcZEAiuAj0XkGNYZwmxnxaOUUurCnDqVwRgTBURV+91TVW4XAbc6MwallFL2c2b3kVJKqYuMJgWllFKVNCkopZSqpElBKaVUJU0KSimlKl10pbNFJAM4Wc+nB+KeJTTcNS5w39g0LsdoXI65FOPqaowJutBGF11SaAgR2WtPQShXc9e4wH1j07gco3E55nKOS7uPlFJKVdKkoJRSqtLllhTebuoAauGucYH7xqZxOUbjcsxlG9dlNaaglFKqbpfbmYJSSqk6XNJJQUReFJHDIhIjIutFpHUt200UkSMickxEFrkgrltFJFZEKkSk1pkEInJCRA6IyH4RcfoapA7E5dL9ZXvNtiLyrYgctf1b4wr1IlJu21/7RaR6qfbGiqXO9y8iPiKyxvb4LhEJdUYc9YjrThHJqLJ/5rsorvdEJF1EDtbyuIjI32xxx4jIYDeJ63oRyamyv56qaTsnxBUiIt+LSJzt/+MDNWzjvH1mjLlkf4DxgKft9nJgeQ3beADHgTDAG4gG+jg5riuBXsAWYGgd250AAl24vy4YV1PsL9vrvgAsst1eVNPf0vZYnpPjuOD7B/4IvGm7PRtY44L9Y09cdwJ/d9XnqcrrjgYGAwdreXwSsAFrJcZrgF1uEtf1wFdNsL86AINtt/2B+Br+lk7bZ5f0mYIx5htjTJnt7k6s1d+qGwYcM8YkGGNKgNXAVCfHFWeMOeLM16gPO+Ny+f6ymQp8aLv9ITDNBa9ZE3vef9VYPwPGikhNS8+6Oq4mYYzZSg0rKlYxFfjIWHYCrUWkgxvE1SSMMWnGmH2227lAHNZ69lU5bZ9d0kmhmt9jZdbqOgHJVe6n8Ns/QFMxwDci8rNtnWp30FT7q70xJg2s/zRAcC3b+YrIXhHZKSLOSBz2vP/KbWwHJTlAOyfE4mhcADNt3Q2fiUhIDY83BXf+PzhcRKJFZIOI9HX1i9u6HgcBu6o95LR95tRFdlxBRDYBV9Tw0OPGmC9s2zwOlAEra2qiht81eEqWPXHZYaQx5pSIBAPfishh29FNU8bllP0FdcfmQDNdbPssDPhORA4YY443Rnw29rx/p+2jOtjzml8Cq4wxxSJyD9bZzA1OjsseTbG/7LEPqzREnohMAj4HerjqxUWkJbAWeNAYc776wzU8pVH22UWfFIwx4+p6XETmATcDY42tM66aFKDqEVNn4JSz47KzjVO2f9NFZD1WF0GDkkIjxOWU/QV1xyYiZ0SkgzEmzXaanF5LG7/uswQR2YJ1lNWYScGe9//rNiki4gkE4PxuigvGZYw5W+XuO1jjbO7AaZ+phqj6RWyMiRKRf4hIoDHG6TWRRMQLKyGsNMasq2ETp+2zS7r7SEQmAo8BU4wxBbVstgfoISLdRMQba2DQKbNWHCEiLUTE/9fbWIPmNc6ScLGm2l+RwDzb7XnAb85qRKSNiPjYbgcCI4FDjRyHPe+/aqyzgO9qOSBxaVzV+pynYPVVu4NI4A7bjJprgJxfuwqbkohc8etYkIgMw/q+PFv3sxrldQVr/fo4Y8wrtWzmvH3m6pF1V/4Ax7D63fbbfn6dEdIRiKqy3SSsEf7jWN0ozo5rOlamLwbOABurx4U1iyTa9hPrLnE1xf6yvWY7YDNw1PZvW9vvhwLv2m6PAA7Y9tkB4A9OiuU37x9YgnXwAeALfGr7/O0Gwly0jy4U11LbZyka+B7o7aK4VgFpQKnt8/UH4B7gHtvjArxhi/sAdczIc3FcC6vsr53ACBfFdS1WV1BMle+uSa7aZ3pFs1JKqUqXdPeRUkopx2hSUEopVUmTglJKqUqaFJRSSlXSpKCUUqqSJgWllFKVNCko5SIi0l9ETotIv6aORanaaFJQynX+jHWB3Z+bOhClaqMXrymllKqkZwpKKaUqaVJQqha2MYAfq9wfLCLfNVU7SrmCdh8pVQsRaYZVjriTMaZcRL4HHja2VbFc3Y5SrnDRr6eglLMYYypEJBboKyI9gKTqX+T2LFpkTztKuQtNCkrVbSfWugx/BCZWf9DYv2hRne0o5S40KShVt53AB8AbxphUN2hHKafSMQWl6mDr7vkB6GGMyW/qdpRyNp19pFTdHgAWN8IXeWO1o5RTaVJQqgYiEi4ihwE/Y8yHTd2OUq6i3UdKKaUq6ZmCUkqpSpoUlFJKVdKkoJRSqpImBaWUUpU0KSillKqkSUEppVQlTQpKKaUqaVJQSilVSZOCUkqpSv8HQV3RUwL/m4AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def huber_loss(actual, pred, delta):\n",
    "    absolute_error = np.abs(actual - pred)\n",
    "    huber = np.zeros(len(actual))\n",
    "    mask = absolute_error<=delta\n",
    "    \n",
    "    huber[mask] = .5*(actual[mask]-pred[mask])**2\n",
    "    huber[~mask] = delta*absolute_error[~mask]-.5*delta**2\n",
    "    \n",
    "    return huber\n",
    "\n",
    "pred = np.arange(-2,2.1,.1)\n",
    "actual = np.zeros(len(pred))\n",
    "\n",
    "absolute_error = np.abs(actual-pred)\n",
    "squared_error = (actual-pred)**2\n",
    "huber_d1 = huber_loss(actual,pred,1)\n",
    "\n",
    "plt.plot(pred,absolute_error,label='Absolute Error')\n",
    "plt.plot(pred,squared_error,label='Squared Error')\n",
    "plt.plot(pred,huber_d1,label='Huber ($\\delta=1$)')\n",
    "plt.title('Comparison of Loss Functions')\n",
    "plt.xlabel('$y-\\hat{y}$')\n",
    "plt.ylabel('Loss')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The squared error tends to be dominated by outliers. For the absolute loss, however, the gradient is the same whether the prediction is close or far away from the true value. The Huber loss combines the good properties of both so that within some neighborhood $\\delta$, the loss is measured by squared error, resulting in smaller gradients as the predictions get closer to the actual values. To mitigate the problem of squared error being influenced too much by outliers, the loss outside that neighborhood is measured by absolute loss instead. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For classification problems, loss functions include exponential loss and deviance. Exponential loss in the binary classification case gives rise to AdaBoost. However, in practice, multinomial deviance tends to be a more robust loss function. In a $K$-class classification problem, $K$ trees are fit at each stage $m$ (in procedures 2A-2D), one for each class. Each tree restates the problem as a one-versus-rest problem. In the binary case, only one tree is fit at each stage.\n",
    "\n",
    "- Multinomial Deviance\n",
    "    - $L(y,f(x))=-\\displaystyle{\\sum_{k=1}^K}\\mathbb{1}(y=\\mathcal{G}_k)f_k(x)+\\log\\left(\\sum_{k=1}^Ke^{f_k(x)}\\right)$\n",
    "        - $\\mathcal{G}_k$ is one of the $K$ possible unordered classes $\\{\\mathcal{G}_1,\\mathcal{G}_2,\\dots,\\mathcal{G}_K\\}$\n",
    "        - $f_k(x)$ is the $k^\\text{th}$ tree's prediction for data point $x$\n",
    "        - $\\mathbb{1}$ is the indicator function. This means the first term will only take one tree's prediction for $x$. This is the prediction of the $k^\\text{th}$ tree where $k$ corresponds to the data point's true class.\n",
    "        - Example 1: $x$ belongs to class $\\mathcal{G}_2$ in a three-class problem ($\\{\\mathcal{G}_1,\\mathcal{G}_2,\\mathcal{G}_3\\}$)\n",
    "            - $f_1(x)=0.3$\n",
    "            - $f_2(x)=0.7$\n",
    "            - $f_3(x)=0.3$\n",
    "            - Then:<br>$\\begin{align*}L(y,f(x))&=-0(f_1(x))-1(f_2(x))-0(f_3(x))+\\log(e^{f_1(x)}+e^{f_2(x)}+e^{f_3(x)})\\\\&=-0.7+\\log(e^{0.3}+e^{0.7}+e^{0.3})\\\\&\\approx 0.85\\end{align*}$\n",
    "        - Example 2: $x$ belongs to class $\\mathcal{G}_2$ in a three-class problem ($\\{\\mathcal{G}_1,\\mathcal{G}_2,\\mathcal{G}_3\\}$)\n",
    "            - $f_1(x)=0.7$\n",
    "            - $f_2(x)=0.3$\n",
    "            - $f_3(x)=0.3$\n",
    "            - Then:<br>$\\begin{align*}L(y,f(x))&=-0(f_1(x))-1(f_2(x))-0(f_3(x))+\\log(e^{f_1(x)}+e^{f_2(x)}+e^{f_3(x)})\\\\&=-0.3+\\log(e^{0.7}+e^{0.3}+e^{0.3})\\\\&\\approx 1.25\\end{align*}$\n",
    "            - The loss is higher in this example because $x$ is misclassified.\n",
    "    - $-\\displaystyle{\\frac{\\partial L(y,f(x))}{\\partial f(x)}}=\\mathbb{1}(y_i=\\mathcal{G}_k)-\\frac{e^{f_k(x)}}{\\sum_{k=1}^Ke^{f_k(x)}}\\text{ for each tree }k$\n",
    "        - These negative gradients are calculated for each of the $K$ trees created (recall that procedures 2A-2D are performed $K$ times for each $m$) and the first term is either $0$ or $1$, depending on whether this tree corresponds to the true class.\n",
    "        - Example: Using the given from the previous example (Example 2):\n",
    "            - Negative gradient for x in tree 1 ($\\mathcal{G}_1$ vs. not-$\\mathcal{G}_1$)<br>$\\begin{align*}&=0-\\frac{e^{0.7}}{e^{0.7}+e^{0.3}+e^{0.3}}\\\\&\\approx -0.43\\end{align*}$\n",
    "            - Negative gradient for x in tree 2 ($\\mathcal{G}_2$ vs. not-$\\mathcal{G}_2$)<br>$\\begin{align*}&=1-\\frac{e^{0.3}}{e^{0.7}+e^{0.3}+e^{0.3}}\\\\&\\approx 0.71\\end{align*}$\n",
    "            - Negative gradient for x in tree 3 ($\\mathcal{G}_3$ vs. not-$\\mathcal{G}_3$)<br>$\\begin{align*}&=0-\\frac{e^{0.3}}{e^{0.7}+e^{0.3}+e^{0.3}}\\\\&\\approx -0.29\\end{align*}$\n",
    "    - Line search is performed to find the population minimizer required in procedure (2C).\n",
    "    - Because this is a multi-class problem, predictions are stored in a row vector.\n",
    "        - Example: Predictions in the previous example would be stored like so: $\\begin{bmatrix}0.7&0.3&0.3\\end{bmatrix}$\n",
    "            - Each tree $f_k$ in stage $m$ updates the vector at index $k$. For example, if after training tree 1 the calculation at procedure (2C) returns $-0.5$ and we're using a learning rate of $0.1$, then the prediction vector gets updated like so $\\begin{bmatrix}0.65&0.3&0.3\\end{bmatrix}$.\n",
    "        - A softmax function is used to turn the predictions into a probability distribution.\n",
    "        - The final prediction is the class with the largest probability.\n",
    "\n",
    "### Regularization\n",
    "\n",
    "One form of regularization is shrinkage, which is the learning rate $\\nu$. Another method combines boosting with bagging. With *stochastic gradient boosting*, only a fraction $\\eta$ of the data is used for training at each step. A typical value for $\\eta$ is $\\frac{1}{2}$ though this choice may depend on the size of the data set. The boostrapping is done by taking random samples without replacement. The rest of the algorithm proceeds as usual. To further prevent over-fitting, column subsampling may also be performed (as in random forests).\n",
    "\n",
    "In the gradient tree boosting algorithm described above, the regression trees that are trained on the negative gradients are ordinary regression trees. The other regularization methods mentioned above do not affect the growth of the tree. Other regularization methods penalize complicated trees.\n",
    "\n",
    "#### xgboost\n",
    "\n",
    "`xgboost` is a popular open-source implementation of gradient boosting. It is highly optimized for large data sets and even sparse matrices. It also support parallelized computations. One of the things that sets `xgboost` apart from the ordinary gradient boosting methods is its split-finding algorithm. It uses a different metric for determining optimal splits than normal regression trees.\n",
    "\n",
    "$$Obj(\\theta) = L(\\theta) + \\Omega(\\theta)$$\n",
    "\n",
    "where:\n",
    "- $L(\\theta)$ is the loss function\n",
    "- $\\Omega(\\theta)$ is the regularization function\n",
    "\n",
    "Optimizing this objective function achieves two things:\n",
    "- optimizing the loss encourages predictive models\n",
    "- optimizing regularization encourages simple models\n",
    "\n",
    "The objective function across the entire model is given by:\n",
    "\n",
    "![Objective Function Across Model](img/boosting/obj_model.png)\n",
    "\n",
    "where:\n",
    "- $\\Gamma$ and $\\lambda$ are regularization parameters\n",
    "- $w_{jm}$ is the value the regression tree at stage $m$ returns at its leaf $j$\n",
    "\n",
    "The regularization penalizes the number of leaves $J_m$ in each tree and the L2-norm of the leaf values $w_{jm}$.\n",
    "\n",
    "The objective function at a stage $m$ is given by:\n",
    "\n",
    "![Objective Function per Stage](img/boosting/obj_stage.png)\n",
    "\n",
    "where:\n",
    "- $T_m$ is the current tree\n",
    "    - this tree maps the input $x_i$ to the weight $w_j$ at the leaf $x_i$ belongs to\n",
    "    - $T_m(x_i)$ in other words is simply the current tree's prediction for $x_i$\n",
    "    \n",
    "The first term in this expression looks a lot like procedure (2C) in the gradient boosting algorithm except with the regularization term added to it. (I don't use $\\gamma$ to represent the current prediction because the weights used in the regularized trees grown by `xgboost` are different as will be shown later.)\n",
    "\n",
    "Second-order approximation can be used to optimize the objective function in the general setting:\n",
    "\n",
    "$$Obj\\simeq\\displaystyle{\\sum_{i=1}^N}\\left[L(y_i,f_{m-1}(x_i))+g_iT_m(x_i)+\\displaystyle{\\frac{1}{2}}h_i(T_m(x_i))^2\\right]+\\Gamma J_m+\\displaystyle{\\frac{\\lambda}{2}}\\sum_{j=1}^{J_m}w_{jm}^2$$\n",
    "\n",
    "where:\n",
    "- $g = \\partial_{f_{m-1}(x)}L(y,f_{m-1}(x))$, the first derivative of the loss function\n",
    "- $h = \\partial_{f_{m-1}(x)}^2L(y,f_{m-1}(x))$, the second derivative of the loss function\n",
    "\n",
    "**Recall: Taylor Series**\n",
    "\n",
    "$$f(x)=f(a)+\\displaystyle{\\frac{f'(a)}{1!}}(x-a)+\\displaystyle{\\frac{f''(a)}{2!}}(x-a)^2+\\dots$$\n",
    "\n",
    "where:\n",
    "- $f(x)$ is a function that is infinitely differentiable at $a$\n",
    "\n",
    "In the case of the objective function, $f(a)$ is the loss function evaluated at the previous prediction $f_{m-1}(x)$ (i.e. $a=f_{m-1}(x)$). Note that these two $f$s are different. The function being approximated by the second-order approximation is the loss function. The function that results in the previous prediction is a sum of tree predictions. $(x-a)$ in this case is given by $T_m(x)$ because the current tree, which has been trained on the negative gradients of the previous predictions, measures the direction in which to move away from the previous prediction $a$.\n",
    "\n",
    "The objective function can be simplified by removing constant terms:\n",
    "\n",
    "$$\\tilde{Obj}=\\displaystyle{\\sum_{i=1}^N}\\left[g_iT_m(x_i)+\\displaystyle{\\frac{1}{2}}h_i(T_m(x_i))^2\\right]+\\Gamma J_m+\\displaystyle{\\frac{\\lambda}{2}}\\sum_{j=1}^{J_m}w_{jm}^2$$\n",
    "\n",
    "Then the objective function can be regrouped by leaf:\n",
    "\n",
    "$$\\tilde{Obj}=\\displaystyle{\\sum_{j=1}^{J_m}}\\left[\\left(\\displaystyle{\\sum_{i\\in R_{jm}}}g_i\\right)w_{jm}+\\displaystyle{\\frac{1}{2}}\\left(\\displaystyle{\\sum_{i\\in R_{jm}}}h_i+\\lambda\\right)w_{jm}^2\\right]+\\Gamma J_m$$\n",
    "\n",
    "We can simplify notation by defining $G_{jm} = \\displaystyle{\\sum_{i\\in R_{jm}}g_i}$ and $H_{jm} = \\displaystyle{\\sum_{i\\in R_{jm}}h_i}$:\n",
    "\n",
    "$$\\tilde{Obj}=\\displaystyle{\\sum_{j=1}^{J_m}}\\left[G_{jm}w_{jm}+\\displaystyle{\\frac{1}{2}}(H_{jm}+\\lambda)w_{jm}^2\\right]+\\Gamma J_m$$\n",
    "\n",
    "This looks like the optimization problem:\n",
    "\n",
    "$$\\text{minimize}\\quad Gx+\\displaystyle{\\frac{1}{2}}(H+\\lambda)x^2$$\n",
    "\n",
    "Getting the derivative and setting to $0$ will give the optimizer $x^*$:\n",
    "\n",
    "$$G+(H+\\lambda)x=0$$<br>\n",
    "$$x^*=-\\displaystyle{\\frac{G}{H+\\lambda}}$$<br>\n",
    "\n",
    "This means that:\n",
    "\n",
    "$$w^*_{jm}=-\\displaystyle{\\frac{G_{jm}}{H_{jm}+\\lambda}}$$\n",
    "\n",
    "Plugging this into the objective function:\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\tilde{Obj}&=\\displaystyle{\\sum_{j=1}^{J_m}}\\left[G_{jm}\\left(-\\displaystyle{\\frac{G_{jm}}{H_{jm}+\\lambda}}\\right)+\\displaystyle{\\frac{1}{2}}(H_{jm}+\\lambda)\\left(-\\displaystyle{\\frac{G_{jm}}{H_{jm}+\\lambda}}\\right)^2\\right]+\\Gamma J_m\\\\\n",
    "&=-\\displaystyle{\\frac{1}{2}\\sum_{j=1}^{J_m}\\frac{G_{jm}^2}{H_{jm}+\\lambda}}+\\Gamma J_m\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Objective Function](img/boosting/xgboost_objective_func.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above equation would be used to evaluate the quality of a tree structure, but since there are so many possible trees to grow, in practice, the tree is allowed to grow greedily. The equation can be used to evaluate the quality of the tree so it can also evaluate the tree at different points in its growth. Gain can be defined as such:\n",
    "\n",
    "![Gain](img/boosting/gain.png)\n",
    "\n",
    "The split-finding can be exact (the algorithm evaluates the above for all possible splits) or approximated (the algorithm only selects a subset of splits). In `xgboost`, the approximate algorithm for split-finding splits features by percentiles then only gets the best split from there. In the case of categorical variables, the column is represented with several one-hot encoded columns instead. This may result in sparse matrices, but `xgboost` is highly optimized for them anyway.\n",
    "\n",
    "Two ways to grow the tree:\n",
    "- Pre-stopping\n",
    "    - Stop growing the tree when splitting results in negative gain\n",
    "- Post-pruning\n",
    "    - Since it is possible that a split that results in a negative gain now can result in nicer splits later on, grow the tree to maximum depth.\n",
    "    - Recursively prune all leaves with negative gain"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}