coordinate_descent_ru.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ограниченная линейная регрессия с помощью координатного спуска\n",
    "*Строим линейную регрессию с неотрицательными коэффициентами, используя только *```numpy```*!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Этот блокнот мотивирован задачей [неотрицательной линейной регрессии](http://en.wikipedia.org/wiki/Non-negative_least_squares). \n",
    "\n",
    "Обычно, когда вы строите линейную регрессию, вы позволяете коэффициентам быть какими угодно. Но иногда вы точно заранее знаете, что коэффициенты должны быть неотрицательны, есть тому физические или экономические причины. \n",
    "\n",
    "Например, в Альфа-Банке я строил модель вероятности дефолта на базе созданных кем-то экспертных правил (\"fraud-триггеров\"), которые были заранее размечены как позитивные и негативные. Все коэффициенты регрессии, которые таким ожиданиям не соответствуют, хочется занулять, а остальные - пересчитывать.  \n",
    "\n",
    "В таких случаях можно запариться и оценивать [байесовскую регрессию](http://en.wikipedia.org/wiki/Bayesian_linear_regression). Но если по каким-то причинам лень, можно решать обычную задачу МНК, потребовав, чтобы коэффифиценты были неотрицательные. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Модель ```Lasso``` из ```scikit-learn``` уже может требовать неотрицательности коэфффициентов. Основана она на алгоритме [координатного спуска](https://en.wikipedia.org/wiki/Coordinate_descent). Мы можем реализовать его самостоятельно, используя покоординатный [метод Ньютона](http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization). \n",
    "\n",
    "Работает метод Ньютона очень просто: для параболы $y=f(x)=ax^2+bx+c$ вершина находится в $\\frac{-b}{2a}$, и в эту вершину можно сразу шагать, чтобы прийти к оптимуму. \n",
    "\n",
    "Если функция $f(x)$ нелинейная, её можно разложить в ряд Тейлора $f(x)\\approx \\frac{f^{''}(x_0)}{2}(x-x_0)^2+f^{'}(x_0)(x-x_0)+f(x_0)$ в точке $x_0$, и попробовать шагнуть в точку $x_1 = x_0 - \\frac{f^{'}(x_0)}{f^{''}(x_0)}$. Там будет минимум параболы, но не факт, что минимум самой $f$. Не страшно! Повторим попытку из $x_1$, и рано или поздно мы доберёмся до максимума. Ну, а в задаче МНК $f$ действительно парабола, так что к оптимуму мы придём сразу. \n",
    "\n",
    "Если $x$ многомерная, то вместо $\\frac{f^{'}(x_0)}{f^{''}(x_0)}$ используется вектор $g_f(x_0)H^{-1}_f(x_0)$, где $g$ -вектор первых производных $f$, а $H$ - матрица вторых производных. \n",
    "\n",
    "Если бы в задаче не было ограничений, то $gH^{-1}$ мне было бы достаточно для счастья. Но делая такой шаг, можно выйти за ограничение. Поэтому я буду брать по одной координате, и идти в направлении $gH^{-1}$, пока не дойду до минимума по этой координате или до ограничения. Пройдя так через все координаты по нескольку раз, я обезательно дойду до оптимума. \n",
    "\n",
    "Вот пример модельки. Для простоты, я не включал в модель свободный член:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "# генерируем данные - центрированную линейную регрессию\n",
    "np.random.seed(1)\n",
    "n = 1000\n",
    "X = np.random.normal(size=(n, 10))\n",
    "y = np.dot(X, [0,0,0,0,0,1,2,3,4,5]) + np.random.normal(size=n, scale=3) # я точно знаю, что истинные координаты неотрицательные"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Начнём с решения обычного МНК, ограничим его, и найдём оптимальные коэффициенты координатным спуском. \n",
    "\n",
    "Всё отрицательное стало нулевым, а положительное немножко под это подстроилось. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.029 -0.043  0.05  -0.042 -0.09   0.834  1.947  3.177  4.091  5.138]\n",
      "[ 0.     0.     0.05   0.     0.     0.838  1.944  3.181  4.089  5.136]\n",
      "Wall time: 13.1 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# начальное решение неограниченным МНК (часть координат отрицательные)\n",
    "beta = np.dot(np.linalg.inv(np.dot(X.transpose(),X)), np.dot(X.transpose(), y))\n",
    "print(beta.round(3))\n",
    "# запоминаем, где мы были на предыдущем шаге\n",
    "prev_beta = beta.copy()\n",
    "# ограничим решение нулём снизу\n",
    "beta = np.maximum(beta, 0)\n",
    "# запустим покоординатный метод Ньютона\n",
    "hessian = np.dot(X.transpose(), X) # это матрица всех вторых производных\n",
    "# крутим алгоритм, пока бета не перестанет меняться\n",
    "while not (np.abs(prev_beta - beta)<1e-20).all():\n",
    "    prev_beta = beta.copy()\n",
    "    for i in range(len(beta)):\n",
    "        gradient = np.dot(np.dot(X,beta)-y, X) # а это вектор всех первых производных\n",
    "        beta[i] = np.maximum(0, beta[i] - gradient[i] / hessian[i,i]) # обновление\n",
    "print(beta.round(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Можно попытаться увеличить скорость алгоритма, рассчитывая не весь градиент, а только одну его координату за раз. В этом случае оказывается, что координатный спуск часто в несколько раз быстрее обычного градиентного спуска, которым обучают нейросети!\n",
    "\n",
    "С 10 координатами это сокращает время примерно на треть, но в более многомерном случае экономия была бы более существенной. \n",
    "\n",
    "Пробуем:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.029 -0.043  0.05  -0.042 -0.09   0.834  1.947  3.177  4.091  5.138]\n",
      "[ 0.     0.     0.05   0.     0.     0.838  1.944  3.181  4.089  5.136]\n",
      "Wall time: 8 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "beta = np.dot(np.linalg.inv(np.dot(X.transpose(),X)), np.dot(X.transpose(), y))\n",
    "print(beta.round(3))\n",
    "prev_beta = beta.copy()\n",
    "beta = np.maximum(beta, 0)\n",
    "hessian = np.dot(X.transpose(), X)\n",
    "resid = y - np.dot(X, beta)\n",
    "while not (np.abs(prev_beta - beta)<1e-20).all():\n",
    "    prev_beta = beta.copy()\n",
    "    for i in range(len(beta)):\n",
    "        derivative = np.dot(X[:,i], resid) / hessian[i,i]\n",
    "        prev = beta[i]\n",
    "        beta[i] = np.maximum(0, beta[i] + derivative)\n",
    "        resid += X[:,i] * (prev-beta[i])\n",
    "print(beta.round(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Больше фактов, алгоритмов и ссылок есть в [презентации](http://www.cs.cmu.edu/~ggordon/10725-F12/slides/25-coord-desc.pdf) Тибширани, создателя Лассо-регрессии.\n",
    "\n",
    "Здесь же хочется нарисовать картинку для двумерного случая. Начнём со случая, когда ограничения не активны:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.019  4.908]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAGPCAYAAAAUUwckAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xm4HGWZ/vHvE/aQEJWALBpiiCIo\niGxBIRhnxHFfRsYdRR3jhhuDouOGiDoKoiAyGhwFQdT5MYqIK6gxrIGgDvsSIMIAEsISEhLWPL8/\n3qqTOn2qu6u7q7q6qu7PdZ0rOd3VXe/p7tPn7ud96i1zd0RERERkuCaVPQARERGRJlIIExERESmB\nQpiIiIhICRTCREREREqgECYiIiJSAoUwERGRGrPgiWa2k5lZ2eOR9RTCpHbM7Flm9rrE9zua2cvL\nHJOIyDCZ2RQzO9LMrgEeBu4FrgOeX+7IymNmm5nZ7Wb27rLHEmtUCDOzQ8zMzWzflOsOMrPHzewC\nM9usjPFJbrYAvm1m+5nZ04CvAbuWPCYRkaEws62BS4DDgZ8D/wwcCLwAuLzEoZXtg4ADp5Y9kNiG\nZQ9gFJjZi4AfAlcBr3D3tSUPSQbg7heb2SLgguiiG4B3lTgkEZFhOgGYARzg7n8uezCjwMw2BD4K\nLHD3R8oeT6zxIczM9gF+Bvwf8E/ufn/JQ5IcuPvrzGxnYDJw5Sj90omIFMXMdgDeABymADbOS4Ft\ngB+XPZCkRk1Htor+SP8KWA282N3/nrLN/mb2BzNbbWarzOxcM5uTst0MM/uhma0ws4fM7K9m9raW\nbWZG06GfMLMPmtktZrbWzC5Ju8/oNr+JbtP69e2W7ZaZ2Sktl82Jt2+5fMLto8tPMbOHWi47xMzO\nM7O/m9nDZnaDmR1hZpMS2yxsM0Zv3X+P+07dNnH9vGibeS1jWQjg7te6++XAttFz4mY2s939Je5j\nWzP7jpn9X/QzLzOzk81sauIx6fTzfjva7unR9x9N2ccu0XWHJi6bZmbHmNnN0X5vN7MzzGz7Dj/v\nk81sqZldbWZP7PJzjT02icv+LbrP97Rcnul1H2377TaPw28y7H+H1ucmerw7Pb7LErc3M/uAhd+3\ntRZ+/841s7mJbTY3s6+Y2U1m9kjLfU1oTWgZ35FdxvKJlG13sfBecL+Z3Wdm3zWzLVLu+11mdkX0\n8y83sx+Y2XYp2z3DzH4UbfOQmd1oZt/IOL7W18urzGyxma2JxvYzM3tmp8egy+Mz18z+28z+Fr1m\n77Dwu/OEjLfv+PyZ2TvavD7fG13+9sRlkyy8N91g639/vmlm09rsO37dtX69MbFN/J59SMttj4gu\nXxh9P6/NfSW/jkzcvuN7TLTNc8zsl2Z2V8v9XNflYY1f025mF0eP6/1mdlbrc514/WzT7s4s/fd2\nYzP7tJldF43/LjP7uYV+3JnW/bE4Jbqf+L30hWZ2YvT8rzazM81s25SxDPL6fS1wh7tfm3K/zzSz\nGRnvJ1eNrYRFD/jvCI/BP7j7TSnbHACcS6iSHU0Ire8F/mRmL3D3xdF204ELgScB3wTuAN4InGpm\nT3T341vu+s3AE4GTovv8AHCeme3p7jekDHcFoYwaOy3jj/m1jNt1cihwPfAbYA2hr+A/gGnAv0fb\nfBH4bvT/6cDXgQXA+TnsPw9fAjbJsmH0ZnQpsBVwMmGKelvCL/CWwKrE5kcCra+bsefG3W80s4uA\ngwmPSdLBwKNEn8rMbHPgT4TetVOBywivp5cBs4HbU8Y6lfAhYmNgnrvfl+VnTNz+rcAxwGfd/TuJ\nyzO97lMcnPh/68/bTtpz8xFgSvT/nQmvsy8B8Zvn6sS2C4B/Bc4DTgGM0Hg8l/Wvv2OA9wE/IDzG\nj0TXz884Rgi/BysT38ev8zQ/JjxfnwJ2j8Y3A3hxvIGF8PZlYBHwsej6Q4G5ZvbcuCJvZs8ivLcA\nfIfweptJqHR8BPgpsDSx768THqcFicuuje7rTYS2i/8FPk34Hf4gcFH03nNLpkdivNcTXqffBf4O\nPCf6eZ8N7Jfh9h2fP3f/vpm9BjjWzM5195vNbBZwLPBzd0/29pwEvAc4Gzg+GsP7gTlmtp+7P5qy\n//NZ/1jFr7WOovf7T7ZcfC3jX//zo/tLvm9fEd2+63tMFBx/A2wU/SzLCL1Mn4oeo06eFP37deDK\n6DZbMPhzTTT+SYTH+J+A/wFOJMw4vBDYM7os+Vj8c/SzJS9rfd88HngQOArYIRrrTtFYH4n2O+jr\nd39gSZvrriW8N8zrch/5c/fGfAGHEF7IryQcJeKEviFrs/0S4B5gq8Rl2xP+EF+QuOzY6L7+MXHZ\nxsBiwh+MadFlM6PtHgZmJrZ9BuEP8hkpY/g9cHPLZQ58u+WyZcApie9fB6wDfhme5s63jy4/BXio\n5bLJKdt9N/q5Nkm5Lv4ZD2nzmPay79RtE9fPi7aZl7hsIbAw8f2e0eNwTrTtzHb3lxjHOmBuynXW\n8jrat9vPR3gzdmCX5P0AfwPOTlx2ZLTdWzrsd+znjV5f5xIC+jMzvv7HHhvCG+gjwIn9vu4T1/0X\nsC7l9fibdvvP+tykPccp1y1o95hF/78D+H3L9W2fw5bt4udlmzav80+kbPvrlv0fFV3+4uj76cBD\nwB+BDRPbvTra7uiWx2wN8PSW/U9qM95lJN4HEpdvBNxJeN/bPHH5c4HHgdOzvIZS7jft/eGt0c+x\nX5fbZn3+ngzcTQhMG0b/3g1sndjm2dF9ndZyPx+OLv/XlH08BnwvZTxvTHmeD0lc9k1C4FySfD23\n3PcpwLIO13V7j3lxtN+DU36HruvyuL4ruu3NwNTE5bu1PtftXt8p+0z+3h4S3eZTnZ631n20ue/4\nvq4g8fcEeGd0+fw8Xr/R62Yd8I0213u757Lor6ZOR34PmAWcTvi0NqFpO/q0sidwqrvfHV/u7rcD\nZwDPt/XTP68A/uLuv09s9wjwDWBzJqbrs919WWLbG4DfEuasW21KeMPOzMw2InzK/gmhopJmUzOb\nnvwipVrk7mui+9zAwjoz0wm/lJsDO/Uyrl733bLtlpaYAu3B1whHCZ3ZbcPo/l9LCA8Tqnge/bb2\n6CeE5y/5KfAFhMpHsqJ5EHCNu/8ww36NUNWZC7zc3btNT4y/sdnehE+rvwA+1HJdL6/7WM+v0Ujm\n56aNg6J/P916RctjNpUQVoflxJb9nxD9+7Lo3xcRXu9fd/fH4o3c/eeEqvPLAcxsK8Jr5RR3vzG5\nA3df1+OY9iT0w5zk7g8m7ucvhCrUy8x6Xz8q8f5gZrZFYlYg3mcnmZ4/d7+LUOHan/CheX/gPe6+\nPHGTV0T/HtNyV98GHiB6TGPRe+QG9P7e+vRoLJ9jfEU26+2zvsfE05L9vG4fjv79truPVe7d/QrC\nB7e05/pJ0fvs5hnu/yDgfkLxYZw+3yMBvuPuDye+/0G0j/h3ZtDX75MI75upswXubu4+r8+xD6Sp\nIWxLQgJ/O3Ax8FUze3LLNjOjf9P+wF1DeEJnJLZtt13yvmI3MtENwBNsYv/CVoyfBsnifdHYWkvm\nSW8nfJpMfr2xdSMLvUGLCJ/G7422i8NDpr6Pfvfdsu0K4CEzu8jMMq1zY2avIvwROzzjuLYilO2v\nzLh9V+6+knCI+JsTbxJvJTynv0hsuiNhWiKLLxGmozYhhOFebE+ojm5O+D1ofdOcGf2b5XUf6/k1\n2sdzk2ZHYHnLH+M0FwKvNrM3mdl2UVCY0uU2g2gNTCsIb/47RBfNjP5Ne4yvTVw/K/o3j9djp31e\nQ2iPmNC31o2ZPdXMfkx4/lcSfldvjq7u9v6Q9fnD3X9K+OAwB/hR9H3STMJr+fqW2z0cjWdmy/Zb\nRf/2+t76FcL073e7bdhG1veYywiV6s+b2V6JD6sbZdhHHGbavb7SnuurCc/dajO720Jvart97Qjc\n0BKaBtX6O/MYcAvZfmd6ef2O3EK1TQ1hn3T3M6JPk/MJb8gndLlNUtYnst12aZ8WJmxrZhsATyFM\nMWTbYQhxnyF8Gu90u3MI/V3Jr9+13NcswienJxD6T14RbXdEtEm/r5+u+07Z9sWEnpmnAOdEb0ht\nWTgc+SvA/7j7RRnHFT8H/X6aa+dUosPFzWwTwifJ/+furZ/Cs+53DqFcfwmwwMwm9zCWuL/sTYQQ\n1EtfVLvX80x6e43289y0G0+Wx2w+4Y/aGYSf/W7ClFJRMv1+t2Ep/8/79dhpn9lvFKo6vyNU9r5M\nqPC8GHhJtEm394eszx8WGv3jA0N26hAQ2u2n1czo32WZ78Rsf8LP+HF3f7yH/aeNpePP7e63Am8j\njPMy1n9YzfIBdFX3TSZ4PeF99tWEfq/DCX9H0mR+3nqQ1+9MO/dE++h48FIZmhrC/hT/x92vAr4K\nvN7Gr6q+LPo37ciLZxKe0FsT27bbLnlfsWekbPt04P6ochJ7NrAZ7acU03yK8Lx+sct2t7v7eckv\nwpx70qsIU02vdPf/dPdfRtv11ADe575btz3X3RcQGl2fCOzTZR/zCZ/YPtFlu6TlhKmLvBd2/R3h\n5zuY8JhOY+LBFTf1sN8j3f37hGn0pxJ6jrK6E3ipu/+Y0OT6VYuOvowsi/7N8ronmprckd5eo/08\nN2mWAk+2sDBlW9EftDcTntsfEP7YtE5b5Wnc73f0geEJhD5A6P4Yx9fHDfd5vB677fM+wuPTi92i\n2/6bu3/Z3c9y93NZXwnrJtPzF/kmsDUhHOzBxCnMZYQ/xuNaJMxsY+BpTHwP3jv6t5fX7bHAH939\nnB5u0yrze4y7/4TwPv44YebmQKLm/i7iqlIvz/X50fvs2e7+LkJV/iUTbh0sBZ4RPbZ5af2d2ZAQ\nQLP+znR8/Uah+UbCa2GkNDWEtTqa8ASdZGZTADwsV3E58LZk1cXCIeRvAS7y9UejnQM818xemNhu\nI0JT6BpCD1XSKy2xVIKZPYPQKP2blu3eQPij9/OMP8cOhB6fo73HI+XaiD/tjX3SiCo5h6ZvPhTx\nWDp9Ep1K6Nn4T3df2mG7caLK6M+Al6ZNefbTMxPd7+OEwHMQITgtY+KRo2cCu1ji8PgO+10U3e81\nhDfpj0R9Xlnc4OuXYvkwYeripMRYe3ndA/wL4X3krIz77+u5aSPuJZsQQpOPWVRR/iHhaM/3RqH/\nmtbb5OjQlucs7rv7dfTveYTH/cPR2OJxvpIQIn4JY9OYfwIOMbMdkzvo47V4OaGZ/H2WOCOImT2H\n8Mf9V33080x4f4hknWLO+vy9hjCF/1l3/xrhKNF/N7O9Ejf5ZfTvYS139R7CNNUvWy5/PeH38K8Z\nx/oGwge/j2XcPlUv7zFmthvhSPSj3P3UrB+APRzpfy3wnmSPl5k9m1CpzPJcG+3fY88kfBBufaz7\nfo+MxprsC34b4YNL/DuTx+v3QmCvtCtMS1SUy90fMrP5hKOVjiZMvUF4kZ0HXGJmJxNemO8jzMsn\n32j+g9DT9AszO4FQbXgDYb2Wj7ZUtyD0f51vZt8i/AE7lPCm/HkYW67gA9H+bwKeZ2bPa7mPp5vZ\ngdEnz9g8wqfQb/XzOKT4LaEv4Rwz+w6hB+lgwlEmw/JUM3sJ4bHfgXAI+R2EqbjntrnNHoRej14q\nRLFPEn6pf29mCwi9Ek9m/WHWy/q4TwhTkocTwvbRKW8YxxCOaP2hmR1I+IT+BMLBGp8lUb1t8WVC\nuPsvC4dppx2Gn8rd7zGzDwNnmNkbok/ekPF1b2YfAL5AaKDdzsKSF7HNCeuzvcbdkwFtkOemdfwL\nLaw39J5o6jz+Q/s8QsXgS9H3nyVUPvbx4ZwNY3vgV2Z2DmHJhncTjs78bTTuFRbWjPoyYWmanxIq\nmh8kvL6SDc8fJDSjL4l+B28iTG2/kVA9z8TdHzWzwwhh9EIzO431h/ivpGXqycLafn/q0qx8HeHD\n69fM7CmEntGXEloGsoyp6/Nn4eCE7xB6d+PH5XDC7+ip0Wv+IXe/Mnp83hO1ZJxLmEl4D+F36dTo\n55pFCFLPJ/RkvjmRG3aO/t3fzK6JGtljLyEcgZfHKX+6vsdELQY/jsbebVYjzeGEn++i6DGOl6iY\n8FxHXmhm9xFmXl4JPAv4eJv7Po0Qir9sZnsQ3ps2JSxR8ROyL6HU6o9m9iNCBeyDhA9Kp0Dvr982\nfg68w8ye5e5Xt1ynJSqG8UWXw9IJh9o/DuyduGwuIZw9SDga5ry02xPeGM8gzD0/RFjL5G0t28yM\n9v8JwovnlmjbxcDzU7br9rUwcZtl0WWvb9nnkQy2RMVLgb8Aa4HbCG8IB9J+2YB47Ie0eYx72Xfy\nZ11H+CR0FvDs6Pp5reMgVB2d0LeR9tzPzPA6eQrwfeAuQji+hfCHYEq311G7ny+67vLo+p3aXP9E\nwno5txHC7/8R3nS2a/fzRpfvQzjc/jNdfq6FpByGTegBuQvYspfXfcbX6LJBnpt2P3Pi+kmEKeqr\no+dqBWH6d7/Ez/EYcGgv7wWtvz/0tkTFswjvBSsJAfV7RMvUtNzHvxIatONx/wDYPmW7nQnVh3sJ\nv4c3AMe1Ge8yUpaoSFz/asIaVWujsZ1FyxInhB5ZJzTAd/td2YlQwV8Zje8MQqBwwrR5t9t3e/7O\njF6DrUt0zCW8Vx/bcl9HEILhI4QPaycmH/vE897t68iW53ktMCPL71N03Sm0WaIi43vMydHj+dSU\nfXZcoiKx7UsI4XVt9PycDezc5vUdf60lhJ/DiZZBSfs5CaHrKMLU5COsf2/eJWUcR9J9iYoXEooH\nKwjvNz8let/r9fXb4fHYMBrnF1KuG/f3dJhf8ZokMgTRFOQthAMD/iPDdk/zNs31ZvZewno283Ie\nphTEzC4mvLGlrjxfNVG15IXuvrDN9S8hBNKZwxxXWaLq1ueAbT3l7BtVYWYvI7RYPMfdcztSeBRY\nWPn+yE6vSQtHe17n7kcOaViNFT0f3wee5+6XDGF/RxDaA2Z5vkd39k09YSJDEPVj7Es0LSIywl4I\n/LhuAUyE9asgvL3UUSSoJ2w0rSZMQ3VaDPAGQt+DjLAofO1JaIKPp5zq4oeE6ZR27iA0IUuFuPtA\nzecj7ia6vybPp/3R2lJhHnpCt++64RAphI0gD0dFvbXLNn8A/jCcEckADiI0ht8IvMnde15le1S5\ne7fX6BWMP3eeSKk8rFLf8Zy27p7XgU0iXaknTERERKQE6gkTERERKUElpiO33GKqP3XrrbpvOAIe\nW/0gbNjuXNQyiEdWP8ikzXs9VWLzPLw6nN92oyl6rKRYj65+cNz3m+g119a6B9c/VhvrccrfY+Fg\nxw1H5LH935tuWeHuXYNLJULYU7feinOPPbrsYWSy4sLoLBhbzy53IDV06/mLAZgyJ3XRY4ncvGgx\n2+7f7axOIvm684JLx/4/64BarMJSiNWLlwAwY64eo9wtDyfhmL5f1hOIFGfr177lb9230nRk7kbh\nya8rvWl1d/OixWUPQRpq2/33GQv/Ny9arNdiG1Pm7MWUOXtx6/mLxz5YSk62ng1bz2bFhZetL4iM\nuEpUwkQkO1XBpEzJ118yiKk6Nt6UOXuxevGSsSCmD5k52no2LF/KigsvG/nCiCphBZi+395jZVER\nkaZSdayzuCoGqCqWt6glaNSrYgphIjWhP3AyqhTGOtMUZUGi6UlgZIOYQlhBVA0rTtzYKhNpKlJG\nmcJYZ8mqmMJYjkY4iCmESaWob0Kk+uIwtu3++yiMtdAUZUFGdHpSIaxAqoaJiHSmMJZOVbECjOD0\npELYMCiI5U5TkuPpj5dUXWsY02taVbHCjFAQUwgr2KgfHltFmpJMp34wqQP1jU2kxv0CjMiaYgph\nw6JqmIhIZgpjE6kqVoCSq2JarHUIpu+390iUPetm9eIlOoWR5GrpwksKvf/Z8/Yt9P7rKA5id15w\nKTcvWtz4RV9bg5hmBnJQ4uKuCmHDtHypzimZkxlz5+jToAykXeCaWdAf+WWLFncMeQponSWrYqAV\n+JMr7iuI5aClIjasMKYQNiSqhomUqzUAFRW22um0v3YBTcFsIoWx9VQVK8CQq2IKYcOmaliuNCUp\nnZQdvLJqNy4Fs/YUxtZTVSxnQwxi5u6F7iAPu8+e5ecee3TZw8jFigsvUwjL0a3nL1YII/wh0tGR\n6yXDy6gGr34sS2lOVygL/WKxJocxWL98j8JYDhIH1PUaxrZ+7Vsud/euf5wUwkqgIJafuAzf9CCm\nEBbUNXy1o1A2nsJYoCCWsyiM9RLEsoYwTUeWRdOSuVCDvsTiANaE8BVr/Vlbe8uaFsh0JGWgXrGc\nFTg9qXXCSqAFXEXy1cQAlmbmAXPGviA8LvFXk+hUSIHWFctRQeuJaTqyJGNPpKphuVBvWHOnJJcu\nvKTx4aub1mnLJlXINEU5/jRvqorlIMP0ZNbpSFXCSjL25GklfZG+KYBl0+QKmVbe1zkoc5djVUwh\nrESalszPjLlzdFLvhlEA609TA1lrGGsiBbEc5RTE1Jg/CtSkL9KTugeGYUmG2GRTf52nK5u+vpia\n9nOUQ8O+esJGgPrD8tPk3rD4j0rd+8LUhF+8ZA9ZnQNZ0/vFtJRFTlLWE1NPWIVoWjJfTZ2WbNIf\nEQWwYqVNV9ZR06coNT2Zk61n9z09qRA2StSkPzB9oqu3uoaBUdWkMAbNbNyPm/ZvPX+xwtig+ghi\nCmEjQkdL5qup1bAmUBVs+JLN/HVt5FdVTFWxXPTYVqQQNkI0LZmPJlfDZh0wZ1yfi0je6l4da3pV\nDBTEBtZDEFMIG0WqhuVC1TCR4tQ5jDW5KqbpyeFSCBsxmpbMR5OrYUAtq2FaF2w01T2MgapiUhyF\nsBGkacl8NHUB1yYdJSmjo65hrOlVMVAQK5JC2Iiavt/eqoblpIlBTKQsdQ5j0LyqmIJYsRTCRp2C\n2ECaOi2pBn0pWx3DWFOrYuoTK45C2AhTf1h+VA0TKUddwxg0K4iBqmJFUAgbceoPG1xcDWtaEFM1\nTEZJaxirurgqpulJGYRCWAWoP2xwTZ2WrJPZ8/Ydd05DqabWRV+rrolVMQWx/CiEVYmC2MBUDRMZ\nDXWaomxqEFOf2OAUwipC/WGDa+q0JNRz3TCpvjpNUWp6sjk/c54UwipE/WGDa+K0pNYNk1FXpynK\nplbFQEGsHwphFaP+sME1dRHXulTD1BdWX3WqioGCmHSnEFZVCmIDa1IQq0s1bPa8fcseghSsLlWx\nJk5PKoj1TiGsgtQfNrgm9oepSV+qRFWxalIQ6425e9lj6Gr32bP83GOPLnsYI2fFhZeF/2w9u9yB\nVNit5y8ee9NogvgPQfyHoarqfjLvGxZe3NftnjHveTmPZDTEU9BVroTGH4DqUpXuJv6A28Q+XIBd\nPn7o5e7e9Y/LhsMYjBRj+n57hyC2fKmCWJ9mzJ3TqCA264A5jflEXjXJ4LXjC/oLG2nhrQ7BbOYB\nc1i2aDFLF15S2SC27f77cOcFl3LzosWNCGKtFbGmhrFuVAmrgRUXXqYQNoD4TaIpQQxCRUzVsNER\nh6d+w1c7N/1p/FRe1QNZ8qCMqoaxplXEoJlVsayVMPWE1YCOmBxME/vDoD5HS1ZdUQEsvs/4K95X\n/FVFdVhXTA37kqQQVicKYn1rWhCLP4VXOYjV4TRGRQawVu0CWRVVPYiBGvYlUAirCR0xObgmlcqh\nPtMhVQ9iwwhgafusehhTEKseBbGJFMJqREFscE1cyLXq1TDpX9XDWB3WFFMQa7ahhTAzO9LMvOXr\n78Paf1MoiOWjKUGsDtOSUP1qWNnqEMagulUxBbHmGnYl7Hpg28TXrkPefyMoiA1G/WHVElfDFMQG\nV+UwpiBWLQpiwbBD2GPu/vfE191D3n9j6GTfg1EQqxYFsXy1hrGqUBCrFgWx4YewWWZ2u5ndYmY/\nNrNZQ95/o2jpisE0NYhVVRX7w54x73kT1vIaJVWsiimIVUvTg9gwQ9hi4BDgpcC7gW2Ai8xsy7SN\nzWy+mS0xsyX3PLBqeKOsIwWxvjUxiFW1GhZTNSxfVayKKYhVS5OD2NBCmLv/2t3/292vcPfzgFdE\n+397m+0XuPte7r7XlltMHdYwa0f9YYNr2tIVoGnJYRr1alisikEsPnKyihTEmqG0JSrcfTVwNfD0\nssbQFApig2vS0hXqDyuHglgxFMSqo4lBrLQQZmabAs8E7ixrDE2iIDY4BbHqqFoQi8/pWJUgtuML\n9lUQGxIFsXob5jphx5rZC8zsaWY2BzgT2Bw4dVhjaDoFsXwoiFWDgth6j69bxz2rVrJ67Zrc7lNB\nbHgUxOprmJWwpwA/IqwV9lPgYWBfd//bEMfQeApig2lioz4oiA1LUUHs/gdX8fbjPs+bj/kst624\ni3tWrcwtlFUpiIGa9auiKUHM3L3sMXS1++xZfu6xR5c9jFpZceFl4T9bzy53IBUVvzHEbxR1F7/x\nx38Iqib+wzuzIstw5H1i78fXreOOe+/mAycdwzpfN3b5JJvEGR87iimbTe77vuPAGAfIURcH8iou\naQLrPxBVfUmZrOIPvFU7QGqXjx96ubt3/QOhc0c2lCpig1FFrFqaXhHbYNIknjr9yZzxsaM49bDP\ncephn+Mb8w9jna/j4cceHei+q9asr+UrqqXuFTGFsAZTEBuMgli1ND2IAUzZbDJbTp3GllOn8cQp\n+S39k1fFblgUxKqlzjMOCmENpyA2GAWxaqliEIvXESvqyMn7Vq8a6xHL8tWpj6wq1TCoTxBriilz\n9qplNUw9YQKoR2xQ6hGrlqr1iEH+fWKr167hzcd8dlyPWBad+shu+tMllekNiy1btLiy/WEQPhCp\nP2z0ZO0JUwiTMQpig1EQq5YqBjHIN4ytXrump56w+1av4iMLjuPUwz7HllOnTbi+qiEM1KhfFVUJ\nYmrMl55panIwmpqslqpNTcaSvWKDTlEme8SyfOXZRzYq6jIt2bT+sLpMTSqEyTgKYoNREKuWZBCr\nUhiLe8UgnzDWdFWrhrZSEKsuhTCZQEFsMApi1TJ73r6VroopjOWjyivqQ3XbAvpVlyCmECapFMQG\noyBWPVUNYjAaYawuAbDqQawp1TCoR/+tQpi0pSA2mCYGsVkHzOHOCy6tbBirchCD9WEsuazFMMNR\n1ZryW1V9WjLWpCAG1a6GKYRu27kTAAAgAElEQVRJRwpig5kxdw4z5s5pTBCD6lfF4unJqvWJtUqr\njhURyOL7rXoAS6p6NQyaE8SqPi2pECZdKYgNLg5iTQljVQ9iUP2qWCxZHYPxgWyQUJa8fZ0CWB2q\nYeoPqw6tEyaZaR2xwWktsepJVkXq8Ae6VS+r3K98eC2fvfgXHPW8VzJtk81qFb6Sqr52GDRv/TAI\nrR+jsn6YFmuVQiiIDU5BrJqqurhrnu55YCVvPe5ITj/sSLbcYuJirXVS9ZX0QUGsTFqsVQqhqcnB\nqWG/muoyPSnNUfUPPv2q0rSkQpj0TEFscE0LYlCfPrE6NO1Ld1VfNyypKU36UL3+sA3LHoBU0/T9\n9g5Tk8uXamqyTzPmzuHW8xezevGSxkxNzjpgDjcvWsydF1xa6U/ps+fty9KFl4wFsaZNUd63+oGy\nh9DVxhttxNSUk4w3zbb771PpDz79mDJnr8p8wFVPmAxEPWKDa1qPGIz/ZF7lMAbN6hVbtXYNb/zq\nZ1jn68oeSleTbBI//vgXBgpidegLi915waWN6g2DcvvD1JgvQ6UwNrgmh7G6BDGofxhbtXYNjzz6\naNnD6Oi+1Q/wwQXHDXwAQR2Okow1tUkfKCWIZQ1hmo6UXGh6cnCanqxuEIv/SDdhinLqZpNhs7JH\nMRwzD5hTm94/TUuOJjXmS27UsD+4pjfsV/2PRJVPBi7N0KQmfQhBbJSb9BXCJFcKYoNrahCrw9GT\nMR1FKaOoytXmQY1qEFMIk9wpiA0uGcSaFsagHkEMxq8tpjAmo6KJ1bBRlSmEmdkBZqb+MclMQWxw\n8cm/oZlVsTpMT8LEKUqFMSlTU6thozotmbUS9kfgSUUOROpHQSwfTQxiUM+qmMKYSLlGLYhlDWFW\n6Ciktqbvt3cIY8uXKowNQEGsHlUxUBiT8m27/z6Nm5KE0ZyW7KUn7ClmNiPtq7DRSW2oKja4eHqy\niX1idauKgcKYSFlGqRrWSwi7DLil5WtZ9K9IVwpi+VBVrD5VMVAYG2V1WjE/japh5eslhM0BZrV8\nPS36VyQTBbF8NDmI1bEqBgpjMlxNbdCH0WrSzxrCHLjV3f+W9lXkAKV+FMTy0dRlLKC+VTFYH8aS\n64wpkInkbxSCmBrzpRQKYvlo6jIWUO+qWEzVsfLUfSoy1sQpSRidacmsIexpwN1FDkSaR0dO5qep\nQQzqt65YmrQwpkBWnKY8tk2ekoyVXQ3LGsJ2AOZGi7ZO+CpygFJ/qorlo8nTk1DvKcpYcqoSVB0r\nUhOqYE03CtWwrCFsIWHB1uS/CxPfiwxEQSwfTZ6ehGZMUcZUHSuGHsPmKbMaljWEbQVsHX2tAfaJ\nLosvL9SkDTcoehcyAhTE8tPkIAbNmKKMtauOKUz0Ln7MmlYFa2pfGJRfDct0Pkh3vyf+v5kB3Ju8\nbBg222oaa+9eOcxdSgniILbiwsvCBVvPLnE01RYHsfhTXtlvNmWYdcAcbl60eCyI1b0HJhkeli68\nZFwQmxlVCCVdUwPYtvvvU/sPKlncev7isffMYeplnbDSTNp0MyAEMWkGVcXyo6pYc6Yok1Qhy66p\nAUyCMj+gZqqEtfDoa6gm77Rr9L8rAVQVa4Dp++0dKmLLl6oiNqAZc+dw6/mLx4JYU6tisH7qpe5V\nsSRVyNIpfElSGdWwTCHMzFaxPnhNAa4ws7Eg5u5bFDC2VJN32pU111+p6cmGUBDLT3J6cvXiJY0M\nYtC8KcpWrYFj6cJLxn3fhFCWDKEKYALhg2kZswVZK2GHFjqKHimINYv6xPKlqtj4qlhTw1isU5UM\n6hXKFL5k1GRtzD+16IH0SkGseVQVy4+qYoHC2HhpVbI6hDKFL8lq2FOSmXvCzOzJwMHAjsBn3H2F\nme0H3OHutxQ1wE4UxJpHQSxfqooFCmPp0gJL6/RlbJTC2W2L/zL276pNNlPwkkzKmJLM2hO2J/B7\n4BbgWcAxwArgQOAZwJuLGmA3CmLNo+nJfKkqtl5rGFMQm6hdMOt01GWRAa3jfp+3B9O30FH10pth\nVsOyVsKOBY53989FTfqx3wLv6GfHZvbvwBeBb7n7QD1nOnKymVQVy5eqYus1+UjKfnSqNHULaEXs\ne8UDK+HCnxe2T6mvYVfDsoawPYF3pVx+J/DkXndqZvsC7wau6PW2nagq1jwKYvlSVWw8hbHBaSpQ\npL2si7WuBZ6YcvkzgeW97NDMpgE/JIS6+3q5bRZxVUwLuzbHuIVdtbhrLpp+MvBWTToNkogM73yS\nWUPYz4HPmdkm0fduZjOBrwD/0+M+FwBnuvsferxdZgpizTN9v721yn7Omn4y8DQKYyL1N8wZgKwh\n7HDgScDdwGTgAmApcD/w6aw7M7N3A7OBz2TYdr6ZLTGzJXffc2/WXYxREGsmBbH8qSo2kcKYiOQh\n6zphDwD7m9k/AHsQwtuf3f28rDsys52ALwFz3f2RDPtcQKiasdfuu/V1miQ17DeTjp7Mn3rF0rX2\njIH6xqRa7rzg0rHXsQxfTyfwdvc/uPux7v7VXgJY5HnAdOAqM3vMzB4DXgC8P/p+k84375+qYs2k\nqlj+4ilKVcXGaz1JuKpjItU2Zc5eQ+kLy7pO2GGdrnf34zLczVlA67v294EbCRWyrtWxQejIyWbS\n0ZPF0HIW6bToq4j0opd1wtYQjoS0lusc6BrC3P1+Qg/ZGDN7ELjX3a/KOI6BKIg1k6Yni6EpyvYU\nxqQKVLEtX9bpyC8B64DzgH3d/WmJr1nFDS9/k3falck77cpmW03T9GTDaHqyGGrcb0/TlDLq1A9W\nrkwhzN0/TVgTbCPgejP7VB49XO4+b9DV8vulPrFm0ppixWhdzkJhbLw4jCWPqFQgkzLp9dfdMPrC\nMjfmu/sd7v4OYF70tdTM3lbQuIZCQayZtKZYcbS2WHeqjsmoUBWsfJlCmJntEX9FtzkCOBk40cwu\nL3KARVMQay4FseKoKtadqmNSFr3ORkfWxvwlhAb8uCk/+f/d8x7UsGk9seZS035xWhv3QUdRtqNG\nfhk2VcFGQ9bpyKcBs6J/0/5fC6qKNZeqYsXRFGV2aVOVqlpInvR6Gi1ZV8z/W9EDGRVaxqK5VBUr\nVrIyBqqKdZKsUqg6JnmJX0eqgo2OnlbMb2VmG5vZ49HXmrwGVTZVxJpNVbFiqV+sN+odkzwpgPWu\nyCMks66Yv47QB5bK3TfIbUQjQn1izTZupX1QVSxn6hfrT1rvGKhCJt3pHJH9mTJnr0I/LGZtzP+X\nNpdvDJye01hGkqYnm2vc9KROe1QIhbH+tJuuBAWypPtWP1D2EEbC8kv/wg777sGKB/Q3rJ2NN9qI\nLTabPPT9mnvbAlf3G4cFW9cUXQnba/fd/NLzzilyF12tuV4VsSYb6xMDhbECqV9sMDcvGj9t0tRA\ntmrtGg76yqdZt25d2UORithg0iR+8ckvpgax1YuXjH1gzGqXjx96ubt3fSPLWglrPE1PNpuqYsOh\n5v3BqEIWTN1sMmcecTQPP/po2UMpzfJL/wLADvvuUfJIRt89qx/gXScdyyOPPgqbDXffCmE90vRk\ns6lXbDgUxgbX9EA2dbPJTB3yH9RRcecFl/LETSerB6wCsjbmn93mqoGOrqwqBbFmU1VseGbMnaN+\nsRw0PZA1iZahqJaslbB7Olz3gzwGUjXJIAaanmwirSs2HGrez5cCWT0ln0cFsOrIuljrO4oeSBXF\nfWKqijXbuClKBbHCKIzlr/WPtRr7q0nVr+pST1gOND0pqooNj8JYcTpVyUChbNQofFVf1p6wKzpd\n7+675TOc6tL0pIAa94dJYaxYaVUyhbLRoKnH4Sn6rB5ZK2E7A2uA7wJa/a4NTU8KqHF/2BTGhkOh\nrHwKX+XodY2wXmQNYc8GjgEOBj4PfNvdHy9sVBWn6UkBVcWGTWFsuBTKhkPBq96yNuZfD7zKzOYR\nwtihZvZxd/9FkYOrMk1PCqgqVgaFsXKkBYTWRn9QMMuiNcwqfNVXT4357r4Q2NvM3gqcaGYfBf7N\n3f9SxOCqTtOTElPj/vApjJUvS7UMFMzSHhMFr2bI2ph/QsrFvwPeCFwKbJTnoOpG05MS0xTl8CmM\njY52wSKtYharY0BT6JJY1krYrm0uL/awgRrRuSclpinKciiMja5O4SwtsMRGPaB1GrtC1+gr+shI\nyN4T9sKiB9IUqopJTFWxciSPdNK5KUdbt6DSqYKWJs/Q1ilgJSlsVVuRR0ZCjz1hZrYh8FTAgNvc\nvbmnqB+AmvYlpqpYuVpPFA4KZFXSS8DpVlUrct8i7WTtCZsKfBV4O7AJIYQ9YmanAR9z9/uLG2I9\nqWlfktS4Xy5NVdafQpP0YhhTkQCTMm53OnAAIYTtRugRexvwfOC0YobWDHEYi6ti0mxxGGP50vXT\nlDI0M+bOGQtkqxcvGdobsYiMnqKnIiH7dOSBwP7u/ufEZVeb2c3AovyH1Sxq2pckTVGWT5UxERmG\nrCHs98BM4M8tl+8AnJfngJpMTfuSpCnK8qU18YMCmUidDbMC3jGEmdn3ov8+CJxqZq8BHkvc9jXA\nLxPb4e7vLGKgTaGmfWmloyhHg6pjIs0xjKlI6F4J2yqx3UbR94+1XLYFMKWQ0TWUmval1YQpSlAY\nK4nCmEh9DbsPtGMIc/dXApjZGcD73P37yevN7J3AP7r7W4obYnNpelJaqV9sdGiqUqSehlUFg+w9\nYa8Gvpxy+eXAifkNR1qpaV/SqF9stKg6JlJ9D/7lf9l008lD3WfWJSouAr5nZgea2VOjrxcD3wUu\nKW54EtNSFpJGS1qMlrQlLrTMhUh1DLMKBtkrYW8mVLx+xfrg5sBZwKEFjEtSqComadQvNnpapypV\nHRORNFnPHXk38AYzmwY8Lbp4mVbKL4d6xSSN+sVGU9pUJSiQiYyKB//yvwBsv89zh77vns4dSViq\nYiXhtEUP5j8cyUpLWUg76hcbTaqOiUirTD1hZjbVzP4TeABYCtwIrDKzk83sCUUOUNqbvNOu6hWT\nttQvNrrUOyYyGsr+vev33JG7oXNHjoxkEFMYk6Tp++2tMDbC4jCm81WKDF/8u1bGNGRM546sCS3w\nKp2oeX/0ad0xkeGJA9iMuXNYvrK89vaslbD43JGtdO7IEaOqmHSiylg1aLpSpHjDXo4ijc4dWUOq\nikk3OpKyGto184MqZCL9Wr14yUgEMNC5I2tNR1BKNzqSsjoUyEQGN0oBDHTuyNpTVUyyUBirFgUy\nkd6NWgADnTuyMVQVkyym77e3mvcrRoFMpLNkE/6oyRrC4nNH/jtwXXTZzsAXyXjuSDP7APAe1jf4\nXw0c7e6/zDxaGYiqYpKFjqSsLgUykfFGOYDBcM8d+X/AEYSFXicR1hw7y8z2dPcrMo9YBqaqmGSh\nMFZtCmTSdKMewGCI545095+3XPQpM3sf8DxAIWzIVBWTrBTGqk+BTJqmCgEMejx3pLuvNLPbgLXu\nvqbfnZrZBsC/EI6qvKjNNvOB+QAznrJ9v7uSLuIwBlcCqopJewpj9aBAJnVWlfAVyxzCzOxdwGeB\np0Tf3wIc6e6n93AfuwIXA5sCq4HXuvuVadu6+wJgAcBeu+/mWfch/UlOUSqISScKY/XR+ocquVI/\nKJRJtVQtgEHGEGZmBwPHAV8AjgI+Ajyb0Ky/qbt/N+P+rgd2B54AvI6wAOw8d7+q55FL7lQVk14o\njNWPqmRSRcnXaZUCGGSvhB0OfMTdv29mnwPOc/cFZnYrcBiQKYS5+yNAfJ6UJWa2N/BR4F09jlsK\npMZ96YXCWD11CmSgUCajoYrVr6SsIWwn4E8pl58FHD3A/icBmwxweymIGvelVwpj9ZU2bakqmZSp\n6uErljWE3Q9slnL5VOC+LHdgZv8B/BK4Lbrdm4F5wMszjkFKoKqY9EphrP7USyZlqfLUY5qsIezX\nwFzCAqu7ALdHl784ui6LbYDTo39XEpaleKm7/zbzaKUUqopJPxTGmkNTl1K0uoWvWNZ1wt6R+P9t\nif9/JeuO3P2QnkYmI0eN+9IPhbFm6TZ1CQplkl1dw1esp3XCzOwfCJUwB65x9z8WMioZaZqilH4o\njDWTQpn0qu7BKynrEhXbAz8D9gTuiC7ezsyWENb6uqPtjaWWNEUp/VIYa7a0P6qtPWWgYNZETQpf\nsayVsBOAx4HZ7n4LgJnNIvR4nQAcVMzwZNSpKib9UhiTWJZqGSiY1VETg1dS1hB2IDAvDmAA7n6z\nmX0I+H0hI5PKUFVMBhGHMYgCGSiMNVy7apmCWT00PXgl9dQTlmJdLqOQWlDjvgxqXHUspkAmKJhV\nWetz1PTglZQ1hP0eOMHM3hQfHWlmM4DjUSVMWmiKUgalqUrJopdgBgpnw5L2+Ct4pcsawj4E/By4\n2czuIBwduT1hra8PFTQ2qTBNUUoeFMakV+3+2CucFaPdY6rQlU3WdcJuA/YwswOBZwJGWKLivCIH\nJ9WnKUrJg8KYDKpTKEg7OjNJIS1Q4MpfTz1h7n4ucG7yMjPbwt0fyHVUUjuaopQ8KIxJEboFtHbh\nI6kOQa3bz6mwlb+s64R9xN2/kXL5S4AFwIy8Byb1oylKyUvqEZWgQCa5yxo8ulXTOikqwGUJj60U\ntIYrayXs01HF6ygAM5sCfAN4I/CZogYn9aQpSsmTqmMyCgYJL4MEuE4UqEZf1hA2D/idmW0B/Bb4\nLuEk3nu6+/UFjU1qTlOUkieFMakqhaXmytqYf5WZHQCcB3wU+DTwFXfXOmEyEE1RSt40VSkiVZG5\nMd/dl5rZ/oTG/GcXNyRpIk1RShFUHRORUZa1Mf9KwtpgAJOBNwHPN7NVAO6+WzHDk6bRFKUUQWFM\nREZR1krYmYWOQiRBU5RSFE1VisgoydoT9vmiByLSSlOUUiRVx0SkbJPKHoBIN3EY22yraWPTlCJ5\nmb7f3usrZMuXrv8SESlYTyvmi5SldYoSVBmTfE2YqlR1TEQKpkqYVMrknXYdVxkTKUJqdUxEJGeq\nhEklqV9MhkGN/CJSJFXCpNLULybDot4xEclb1nXCDut0vbsfl89wRHqnfjEZJvWOiUhe2oYwM9vZ\n3a+Nvj0WWAMsB6xlUwcUwqR0CmMybKnLXIACmYhk0qkS9h0zu9Pd3wB8CfgQ4dyRn3H3u4YyOpE+\naLFXGTZVx0SkH51C2KuAu81sS3f/tJmdBHwRuN7MjgGOdfeHhzJKkT6oeV/KoGZ+EcmqU2P+U6J/\nHwJw9zvc/R3AvOhrqZm9rdDRieRAzftSFjXzi0gnnSph5wIfd/cHzWyPluuOAF4BnGhmH3b3PQsb\noUgO1C8mZdJ0pYikaRvC3H3bxLdLCA34cVN+8v+7FzM0kfwpjEnZNF0pIrGsi7U+rdBRiAyZmvdl\nFOjoSpFmyxTC3P1vRQ9EpAxq3pdR0Ha6EhTIRGos82mLzGw34HBgF8J05DWEIySvLGhsIkOjMCaj\nQoFMpDkynbbIzF4F/Bl4KvBr4DfADODPZvbK4oYnMlw6klJGSXx0pU4mLlJPWSthRwNfdPfPJS80\ns6Oi636R98BEyqLmfRlFaugXqZ+sIewZwGkpl58GfDy/4YiMDoUxGVUKZCL1kDWELQf2BFrr4HsC\nhZ/C6P6H1xW9C5G2FMZklCmQiVRX1hB2MuFckrOBiwiN+fsTGvWPKWhs603agLNvWgPAq3acXPju\nRNIojMmoUyATqZZeesJWA/8GfCG67A7gc8AJBYxrAps8DV+zUmFMSqc1xqQK2gYyUCgTGRFZ1wlz\n4OvA181sanTZqiIHlsYmh+pDHMYUxKRMWtZCqkLLXoiMpszrhAGY2Y7AztH/r3H3mwsZVbdxqCom\nI0RhTKpEgUxkdGQKYWa2JfBfwKuAdesvtnOAd7r7PQWNr/2YVBWTEaMwJlWjQCZSrqyVsO8Cs4G5\nwOLosjnAfxKa9v85/6FlE4exs28Kf/AUxqRsCmNSRclABuojExmGrCHsn4B/dPeLE5ddaGbvAc7L\nf1i90xSljJrJO+2qIymlslQlEyleptMWAXcDD6ZcvgbINBVpZp80s8vM7AEzu9vMfmFmz8460Ez7\nmDwtURlbk+ddi/Rl8k676lRIUnltT5+kUyiJDCRrJewo4BtmdrC73w5gZtsDX4uuy2IecBJwGWDR\n7c4zs13c/d6eRt2FqmIyarTGmNSFpi1F8pM1hH0EmAksM7Pbo8u2Bx4CtjazD8UbuvtuaXfg7v+U\n/N7MDgZWAvtRwLkn1bgvo0hhTOpG05Yi/csaws4sYN9TCdOh9xVw32PUuC+jSGFM6khVMpHeZF2s\n9fMF7Pt44K/AxWlXmtl8YD7AVttuP/DONEUpo0hhTOqsY5UMFMqk8XpdrPUfgF0I54682t0X9rNT\nMzuOcO7J/d398bRt3H0BsABg9rOe4/3sZ8J+NUUpI0phTOoutUqmUCYNl3Wx1u2BnwF7Es4ZCbCd\nmS0BXuvud7S98cT7+jrwRuCFZa64D5qilNGjMCZNoVAmkr0SdgLwODDb3W8BMLNZwOnRdQdluRMz\nO54QwOa5+3W9DzdfmqKUUaUwJk2jUCZNlDWEHUgITrfEF7j7zdFRkb/Pcgdm9i3gYOA1wH1mtk10\n1Wp3X93DmHOlKUoZZQpj0lQKZdIEPfWEpVjXfZMx74/+bQ1tnweOHHAcA1NVTEbZ+lMhgU6HJE2k\nUCZ1lDWE/R44wcze5O63AZjZDMIRjpkqYe5u/Q1xeFqrYqAwJqNH56YUmRjKIGVJDFAwk5GWNYR9\nCPg5cLOZ3UE4OnJ74IroulrRFKVUQWsYAwUyaTZVy6Rqsq4Tdhuwh5kdCDyTcNqha9x9JE7eXRRN\nUUoVqG9MJJ2qZTLqsi5R8TbgJ+5+LnBusUMaLZqilKpQGBPpLlO1DBTMZCiyTkd+H/gNsLzAsYw0\nTVFKVSiMiWTXtlqmYCZDkDWEjXxT/bBoilKqQmFMpD9pwQzaTGWCwpn0rZclKp5nZqkn23b3RTmN\npxI0RSlVojAmko+eqmagcCZd9RLCftbmcgc2yGEslaMpSqkShTGR/HWsmimcSRe9hLBt3L2xPWGd\n6FyUUiVa+FWkeApnkkXWEOaFjqIm1C8mVaO1xkSGq104gw49ZzGFtNHTLlBnpMb8nGmKUqpIU5Ui\n5esa0Lr9wVdIK06Hxz71eTs2291mXax1Ura7k5imKKWKFMZERlOngBbrWkmLKayNd8+y9f8+PCV1\nkyyPfz+yLta6ibs/HP1/e2A+MBk4293PL2RkNaEpSqkihTGR6skaFDKHtVZVCW99ThE+ae/dmf6k\nJ+Y8mM46hjAz2wn4KfBMM7sCeAthxfwtgHXAR83sIHc/q/CRVpiWtJCqSmviBwUykSrrt6rTd3gr\nQS8/42P33gffLnAwHXSrhB0L3Al8HHgj8Cvgt8C/Rtd/E/gEoBCWgfrFpMpUHRNptqKm5JqsWwjb\nFzjQ3f9qZouAlcBJ7r4OwMy+CVxS8BhrR/1iUmUKYyIi+ejWcL8lcAeAu68CHgTuTVx/HzC1mKHV\n3/owtmZsmlKkKibvtOtYINtsq2ljgUxERLLJ0pjfukaY1gzLkfrFpOrUNyYi0p8sIex0M3s4+v+m\nwMlmFpdtNilmWM2jfjGpA01Viohk1y2Endry/ekp2/wgp7EI6heTelAYExHprmMIc/d3DGsgMp7W\nF5M60FSliEh7Wgl/hNnkaWrel9pQI7+IyHgKYRWgMCZ1ojAmIhJkPYG3jAA170udaKpSRJpOIayC\n1LwvdaNGfhFpIoWwClPzvtRNWhgDBTIRqSeFsIrTYq9SR8mpSlXHRKSu1JhfE2rel7pSI7+I1JUq\nYTWj5n2pKzXyi0jdKITVlJr3pc7UyC8idaAQVnNq3pc6U3VMRKpMPWENoH4xaQL1jolI1agS1iA6\nklKaQNUxEakKVcIaSJUxaQpVx0RklKkS1mA6klKaQtUxERlFCmGiIymlUbQqv4iMCoUwGaMwJk2i\nVflFpGwKYTKBlrWQplF1TETKoBAmqXQkpTSRqmMiMkwKYdKRwpg0lZr5RaRoCmGSicKYNJmmK0Wk\nCAph0hOFMWmydtOVoEAmIr1TCJO+KIxJ06l/TEQGpRAmA1EYE1H/mIj0RyFMcqEwJhKof0xEshrq\nuSPN7AAzO9vMbjczN7NDhrl/KV7reSlFmio+b2XruSt1/koRiQ27EjYFuAr4QfQlNaXV90XWU/+Y\niKQZaghz918BvwIws1OGuW8px7pNprDynrs55ZLlALx45mYlj0ikZNO2Hvvv2hX3jP3/oXtWlTEa\nkcZbft/9pe17ZHvCzGw+MB9gq223L3k00q+V99zNO1+0R9nDEBERaWuDSZPYdOONh77fkQ1h7r4A\nWAAw+1nP8ZKHI32atuVWfO+8P0+43NeGT/2qjImMt/ama8d9rwqZSPE23XhjnjBl86Hvd2RDmNTD\nBhtswJO23iblmnDZkgfVMyYyzjZPHvvvmuuvhOlbjn2vHjKRelEIk1KpgV+kPa3QL1JvCmEyEmzy\nNK0xJtKBAplI/Qw1hJnZFGB29O0kYIaZ7Q7c6+63DnMsMnq04KtINgpkIvUw7ErYXsAfE99/Pvo6\nFThkyGOREaUwJpKdAplIdQ17nbCFgA1zn1JdCmMivVEgE6kW9YTJyIvDGKiBXyQrBTKR0acQJpWi\noylFepcMZMGVY/9TIBMpj0KYVJLCmEj/xoeyK8ddp1AmMjwKYVJpCmMig9G0pUh5FMKkFhTGRAan\nQCYyXAphUiutYQwUyET6oUAmUjyFMKklLW8hkp9Ojf2gUCbSL4UwqTWFMZH8qUomkg+FMGkEhTGR\nYiiQifRPIUwaRWFMpDiathTpjUKYNJLCmEjxVCUT6UwhTBpNYUxkOFQlE5lIIUwEhTGRYetUJQOF\nMmkGhTCRBIUxkeFrrZJp6lKaQiFMJEUcxkCr8IsMm6YupSkUwkS60Cr8IuXS1KXUlUKYSEaaqhQp\nX7epS1Aok+pQCBPpkSKj3+kAAAv1SURBVMKYyOhQP5lUmUKYSJ8UxkRGj/rJpEoUwkQGpDAmMroU\nymSUKYSJ5ERHVIqMvm6hDBTMZHgUwkQKoCMqRapBjf5SJoUwkQJpqlKkWhTKZJgUwkSGQGFMpJoU\nyqRICmEiQ6QwJlJtE3vKQH1l0i+FMJESpDXxgwKZSBWpWib9UggTKZmqYyL1omqZZKUQJjIiFMZE\n6kvBTNIohImMGIUxkWbIMo0JCmZ1phAmMqLUNybSLGnVMgWzelMIE6kAVcdEmknBrN4UwkQqRGFM\nRNL7yyCtxwwUzkaZQphIBSmMiUgrVc2qRyFMpMLUNyYinahqNtoUwkRqQtUxEcmq13AGCmhFUAgT\nqRmFMRHpV7tw1m5aExTOBqEQJlJTmqoUkby0r5yBqmf9UwgTaQBVx0SkKApo/VMIE2mQtDAGCmQi\nUox+A1qs7kFNIUykgZJTlaqOiUgZOge0zn1osaqHNIUwkYbTVKWIjKJuIS3oXk2D0Q1rCmEiAqiR\nX0SqJ1tQg6xhLWkYwU0hTEQmUHVMROoke1gLskyF5mHoIczM3g98DNgWuBr4iLufP+xxiEh3auQX\nkSbqNbT1a9JQ9hIxszcAxwNfAp4LXAT82sxmDHMcItIbmzxt7Avg7JvWjAtlIiLSu6GGMOAw4BR3\nP9ndr3X3DwJ3Au8b8jhEpE9pYUyBTESkd0MLYWa2MbAn8LuWq34HPH9Y4xCRfKRVx0REJLth9oRN\nBzYA7mq5/C7gRa0bm9l8YH707cOv2W27q4odnhRoOrCi7EFI3/T8VZuev+rSc1ddO2TZqIyjI73l\ne0u5DHdfACwAMLMl7r7XEMYmBdDzV216/qpNz1916bmrv2H2hK0AHge2abl8ayZWx0RERERqbWgh\nzN0fAS4HDmy56kDCUZIiIiIijTHs6cjjgNPM7FLgQuC9wHbAt7vcbkHRA5NC6fmrNj1/1abnr7r0\n3NWcuU9oxyp2h2Gx1o8TFmu9Cviouy8a6iBERERESjb0ECYiIiIiw1+sVURERERQCBMREREpRWVC\nmJk9ycy+aWbXmdlaM7vNzP7TzLYse2ySjZnNN7M/mtn9ZuZmNrPsMUl7ZvZ+M7vFzB4ys8vNbG7Z\nY5JszOwAMzvbzG6PftcOKXtMko2ZfdLMLjOzB8zsbjP7hZk9u+xxSTEqE8IIR1FuT2jq3xV4K3AA\n8KMyByU9mUw4TdWRJY9DujCzNwDHA18CnktYRubXZjaj1IFJVlMIBz59GFhb8likN/OAkwin8/sH\n4DHgPDN7UpmDkmJUujHfzF4GnAM8wd0fKHs8ko2Z7QVcBjzN3ZeVPBxJYWaLgSvc/d2Jy24EznT3\nT5Y3MumVma0GDnX3U8oei/TOzKYAK4HXuPsvyh6P5KtKlbA0WwAPAzpzsEhOzGxjYE9C1TLpd4RP\n5yIyPFMJf6vvK3sgkr/KhjAzewLwBeBkd3+s7PGI1Mh0YAMmnk7sLiaedkxEinU88Ffg4rIHIvkr\nPYSZ2dFR42inr3ktt9kc+AVwO6FHTErSz/MnldHaq2Apl4lIQczsOGB/4HXu/njZ45H8Dfu0RWm+\nAZzeZZtb4/9E8+O/ir59hbs/VNTAJJOenj+phBXA40ysem3NxOqYiBTAzL4OvBF4obvfXPZ4pBil\nhzB3X0F40+/KzKYCvyZ8In+Ju68ucmzSXS/Pn1SDuz9iZpcDBwL/L3HVgcD/lDMqkeYws+MJAWye\nu19X9nikOKWHsKyiAPY7QjP+a4DNo2lJgHvd/ZHSBieZmNk2hOrKM6KLdol6+25193vLG5mkOA44\nzcwuBS4E3ktYJubbpY5KMolmDGZH304CZpjZ7oT3SlWmR5iZfQs4mPB37r7ofRNgtQoP9VOZJSqi\nvqI/trn6he6+cHijkX6Y2ZHA51KueocOnx89ZvZ+Qs/ltoQ1pz7q7ovKHZVk0eH98lR3P2S4o5Fe\nmFm7P8qfd/cjhzkWKV5lQpiIiIhInZR+dKSIiIhIEymEiYiIiJRAIUxERESkBAphIiIiIiVQCBMR\nEREpgUKYiIiISAkUwkSkFGa2q5ndaGZTzWxHM7ux7DGJiAyTQphIRZnZKV1OnD6z7DF2cRVwJXA/\ncD3wrXKHk87MVpvZIUVtLyLNVZnTFolIqvMIpzhJeg7wmxLG0hMPK0X/s5ltBTzk7qvKHpOIyDCp\nEiZSbQ+7+9+TX8A9rRuZ2S5m9kszW2Vmy83sR4lz0sVVtXPM7NNmdldUzfm+mW2W2GZhosr2kJld\nZWavS1x/pJldlTZIM5sX3W569P0hZrYawN3vdvdVZnZ+tM1e7X7YaAwnJr5/S/QzzU1cdoCZLY7G\neJeZfd3MNo6u287Mfm9mK83sQTP7g5nt27KPfc3sWjN7wMz+Lbr4GWZ2dbSvrwy4vZvZQYnv/zW6\nLPlzLTOzw1tud6KZLUx8P/YYJi6b8Bhmfe5b7ueg+PQ50X7aVlwRkb4phInUnJltCywiTP/tA7wI\nmAKcbWbJ94AXEKpo/wi8Dngx8JXx98b3CeeS3BO4CDjdzDbKYYz/DOze421eDSwADnL386PLtgd+\nDfwFeC7wLuBNwJfjmwGnAHsDc4BLgT+a2a7R7TcC/h9wSbTNLGAz4F+AdxCqju83szf1s33Kz7A5\ncBQw8ImZ0x7DHp77Tn5CeM63BT4C/F/i+20HHbdIkymEidTf+4D/dfcj3P1ad78CeBshNCSrTo8T\nTqZ+lbv/FjgCeE8UFGJromrbdcAdwAPR7foWBZn/YGLg63SbFwFnAG+Lxhp7P3An8P7oZz0H+ARw\nqJlNdvfb3f00d78h+jk/AfyBcKJygJcAWwMfdvfrgQ8CDwNfcfdL3f0s4DRgfp/bt/oYcA1wedaf\nvc3j0e4xzPrct+XuaxNV1pXA4y2VVxHpk0KYSP3tCRwQTTGujqawbouu2zGx3RXunqzIXAxs3LLN\n/Oj2DxGqIm9193WJ63eO9rHSzK6Lpiity/jeTwhzP+zh5zkLeAS4sOW6nYGLW8Z0QfRzzIaxKczk\nY3EgsEe07WzgDnd/ACC6n3XAY4n7uya+rz62HxNVqQ4DDm+9rg/tHsOsz/1LWrY5LYcxiUgXaswX\nqb9JwC9J/2N/V4/39RPg88AmwBuAH5nZLu6+PLr+JuBl0T73BL5H+KN/U9qdmdkTgM8Qpj+z9hft\nS6g4vRz4DvDq5F12uJ/48rOBxYnL/x3YLeO+k/sZdPujgTPd/a/dc2qHO+78GGZ97hcxvlr3EuCb\nfQ9KRDJRCBOpvz8Drwf+5u6PdthuVzPb3N0fjL7fl1BtSgaole6+FEIjPiHAHACcGV3/SHw9cIOZ\nzSf0ZqWGMOBTwAXu/ifLvqTGj9z9RDP7GXC1mR3s7nHl5hrg9WY2KVEN2z/5c0RHYa6KfoZJhCrY\nn6Ntbwa2M7Mt3P2B6PpJjH+v3CXx8/S6fWw3Qt/YMzP+zJ10egyzPvdrEs8bZqZpRpEh0HSkSP19\nC5gG/MTM5pjZLDN7kZktMLOpie02BL5nZs8yswMJPUYnJ0IZwGQz28bMZhCm0iYR1viKmZltamaT\nLRyxuAehKTzNJsB7Wd+PldW9AO5+O/Bh4Hgz2y667iRgO+AkM9vZzF4e/RwnuvsaM9vBzN5nZs80\ns90IBxrsHG0DYWmPe6L73Ak4IRrnEWa2t5m9htBs/199bh87DPi6u9/R4efcMHosNzWzTYENgEkW\nHekZ6fYYZn3uRaQECmEiNRf9od+P0Kv0G+Bqwh/nh6Ov2J+i6/4I/IzxDeuxdxAa328kHHn4Tne/\nMnH9LsBaQqXpR4SjF09uM7QNgVPc/YYBfrZTCT1fC6LvbwdeSqi+/ZUwHfojQsUOQkXsIMLRjIsJ\nfVEHxmNw94cJlaP9gCWEStda4L8JR1WeHu3rtH62T1gFHNPlx/tydF/x13uBucDvEtt0fAx7eO5F\npAQW1ksUkSYzs1OA6e7+irLHMmqiRvVD3f2UIrbvcSzzgCPdfV7e9y0iw6dKmIhIdawjVPNEpAbU\nmC8iUhHuvoiwiK6I1ICmI0VERERKoOlIERERkRIohImIiIiUQCFMREREpAQKYSIiIiIlUAgTERER\nKcH/B70H/28qfdsMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cfff410208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "matplotlib.rcParams.update({'font.size': 14})\n",
    "\n",
    "np.random.seed(1)\n",
    "n = 1000\n",
    "X = np.random.normal(size=(n, 2))\n",
    "X[:,0]+=X[:,1] * 0.7 # correlated X for more fun!\n",
    "y = np.dot(X, [1, 5]) + np.random.normal(size=n, scale=3)\n",
    "\n",
    "beta = np.array([-1, 1]).astype(float)\n",
    "prev_beta = beta + 1\n",
    "# improve the solution by restricted coordinate-wise Newton descent\n",
    "hessian = np.dot(X.transpose(), X)\n",
    "betas = [beta.copy()]\n",
    "while not (np.abs(prev_beta - beta)<1e-20).all():\n",
    "    prev_beta = beta.copy()\n",
    "    for i in range(len(beta)):\n",
    "        grad = np.dot(np.dot(X,beta) - y, X)\n",
    "        beta[i] = np.maximum(0, beta[i]-grad[i] / hessian[i,i])\n",
    "        betas.append(beta.copy())\n",
    "print(beta.round(3))\n",
    "betas_history = np.vstack(betas)\n",
    "\n",
    "h = 0.01\n",
    "xx, yy = np.meshgrid(np.arange(-2, 3, h), np.arange(0, 6, h))\n",
    "Z = np.mean((y[:, np.newaxis] - np.dot(X, np.vstack([xx.ravel(), yy.ravel()])))**2, axis=0)\n",
    "Z = Z.reshape(xx.shape)\n",
    "cm = plt.cm.RdBu\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.contourf(xx, yy, Z, cmap=cm, alpha=0.5, levels=[8, 8.55, 9, 10, 12, 16, 25, 40, 70])\n",
    "\n",
    "plt.plot(betas_history[:,0], betas_history[:,1], color='k', lw=1.2)\n",
    "plt.xlabel('Первый коэффициент')\n",
    "plt.ylabel('Второй коэффициент')\n",
    "plt.title('Координатный спуск кодится просто, а сходится быстро (:');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "А теперь посмотрим, как наш алгоритм обходит ограничения:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.     1.959]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAGPCAYAAAAUUwckAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xe4JGWZ9/HvfYYJMAlkhiwMQZLI\ni4KCShhcWTEnzKLguqisi1lXRUXXxCsGzGJCwcC+uiJmSSNBGMRIEHCAAQRkGECYnM79/vFUMzV9\nOlR3V67f57rONXOqq6uerq7T53fu56mnzN0RERERkXyNFd0AERERkSZSCBMREREpgEKYiIiISAEU\nwkREREQKoBAmIiIiUgCFMBEREZECKISJiIiIFEAhrKLM7Dgz8x5fzyq6jSIiItLdZkU3QEZ2CnBz\nh+V/zLkdIiIiMgCFsOr7lbtfWXQjREREZDDqjmwIM3upmV1pZivM7J9mdpmZPbdtnfndujdj68yL\nlh3X9tx3RcsXJGzPUWZ2sZkti9p0sZk9qcN625vZV8zs72a2xswWm9lXzWxmgi7Zh9tpZmea2eoh\nDl2rHTub2XfMbKmZrTazP5nZq9rWmdejHW+M1mkd41eY2QfN7C4zW2lm55vZXm3b2z9q983RPu81\ns++Z2SPb1msdh0M6tNvN7Msd1p3Xtt4vouWnxJYtaH8/e7z/jzKz70dtXGNm15jZ8W3rtF77Szu0\nc7GZ/bLDuvPb1vtStPzMtuWzzOy0aDtrzOx2M/u4mU1t31dSZjbFzE42sxuibd5jZj82s0e3HYtu\nXzdE641F7flxh31MMrO7zewHsWVmZv8RnWOronPufDM7LLZO+3s1yczOtfDzdFBs+XPM7Kdmdmf0\nGm4zs1M7HZcB30OP7yd6bBszW9vl/NgzOneXROfy38zsM23rHBGdc/e1Hcdf0kfb+huiY/pdM9up\nw7ov7fJ+Tfh8MLP9zOw8C5+ZK83sCjN7epc2nNllu2fG1jklWrZdv9fUYfuJf36iZVPM7H1mdmP0\nfv7DzL5sZlt2eG6ndsfPyTELn/E3Rdu608w+Z2azB30dsilVwqpvtpnNaV/o7ktb/zez9wIfBn4H\nfAhYCRwIPA2Y8IsB+Cjw1+j/LwCe36sB0f7fnbTBZvYy4GzgYuBkwh8DrwEuMrMj3H1htN52wFXA\nXOCrwLXA9lF7tgYuAY6Nbfq90b8fiS37bdJ29WjvHOBy4BHA54C7gJcC3zKzrdz99LannAP8tG3Z\nVW3fvwuYBJwGbAW8CbjYzB7j7vdF6xwF7E04Vn8H9gBeDzw+Wm/VqK8ten1PBY4e4fn7EI7PUuBT\nwD+BZwLfMLPZ7v6ZXs8fYD97A6/tsHxzYAGwK3AGoXv+AOBthOP3vCH2NQacR/gZ+SHweWAL4EjC\nz851sdU7vd/vBQzA3cfN7CzgHWa2dez9BXgqsB1wVmzZGdHrvAA4M9rOk4DDgEu7NPlLwNOBZ7r7\n1bHlrwHWE87bB6LtvB14JPDy2Osd9D1cDRwPxPf1KmADMDm+ooXQenn07VcI78884CXAm6N15gE/\nj/b/MeAf0fqf7vJ6O/kR8L+En6v/A/wnsBNweJf1459zJwBPaGv3nlG710TtWE54zT81s2Pc/Ucd\ntrkUeEvs+0HanxozM8LxeArwNeAaYHfgjcCBZvYkd18Xe8pfgE+0bea22P+/CLyO8DNxOrAfcCJw\nsJk9uW1bMgh311cFv4DjAO/xNSNab3fCB+PPgc3atmFt3z81eu7hsWWnhNPk4e/nRescF1v2OcKH\n5tXAgj7tng7cD5zZtnwL4FbgwtiyM4Fx4LAO27EOyxZ023+0rdVDHuvTotf8L7FlU4CFhA/m2W3H\n5r96bGt+tM49wJax5U+Jln80fkw6PP/QaL1XdDgXDumwvgNf7rDuvOj7MeBPhBDhwCmxdS8CftO2\nvU7v/6+B69vbSwgnDwHT2177Szu0czHwyw7HaX5s2U8Ivyxuj58/wHsIf1js07bNE6NtTDh/Bvj5\nem+3c6/X+x2dizfEvt8zWvfEtvXOIvzintz2us/odc7H3yvgvwk/4y/p8JxO59DJhJ+rnUZ4D79L\nCHXTYuteB3ynw/mxIHp/HtW27bHY/0/o9F61nxc93q9Nzt1o2f8AyzuseywTP+fOpO3zAfgBsA7Y\nO7ZsFiGc3BZvf/TY2cAtHdofP1dPifa93RDnZOvYJ/n5eXn0Hh/Ztt4zom0cm/QYEwKXA2e1LX9T\ntPy1g74WfW38Undk9Z1EqJi0f7WqJM8n/KI9xd3Xx5/o0U9SzJTo3zVJd25mjyL8hfQBQiDp5yhC\n5ec7Zjan9UUIYRcAh5nZ5KgS8XzCh8OEv/47tD1pe1v73HyApz0L+KO7Xxjb/1rgM4RQOX+Ipnzb\n3f8Z295FhF9iz4gtWxlr9wwz2xq4gVClOLDDNmfHj2mnCmkHryJ8yL6zw2NLCJWErsxsK0J4/x9g\ni7Z9/wKYCRzU9rSZHdrZ87PIQrfks4B3EH65xL2YULG4t22bF0SPP6XXtrs4hnCcT2t/YJhzz91v\nAq4kVrk1s+mEc/wc31hJOCb69+Qk+zWzN0TrvtXdz+nwnJXRemNm1qqaX0Korj0uemyY9/DnhM+J\n50XbOIRQqf1uW/vmAkcQgsjf2toWfx9nRv8uZXittm9rZkcTfi4v6rDetOjfrsMTzGwSoTr8E3e/\nIdbmh4AvAzsDj2l72hSSf3Y+Imrr9ITrxyX5+XkxcBNwTdt6VxE+pwf5mWhdad9eKfsyIaA/c/CX\nIC3qjqy+33nvgfm7R/9ek2BbrbECScJUy6nAIkLJ+2UJ1t8z+vfXfdoxRvirM0m7k5oK3Nv6xszu\nBD7n7qf2ed48QjdHu+tjjw/qbx2W3UTo7mq1byvg44RfzI9oW3dLJuo7diYuCqIfBr7u7teHHoxN\n/BZ4iZm9Gfg+oVtrq7Z1HkX4hf6B6KuTbdq+PyP6and9h2WtrpXTgPPd/Vcd2rknofvp3vYHuuw/\nid2Bm9w98R8kCXwL+JKZ7e7uNxMC2HQ27YrcHVji7ksSbO/pbAxH7ecH8HBX4CcIgaT9D4/WOTTM\ne7g+avdrCOfGawhdVfe1rbdb9G+/n+NWd+VpZvYuNnZHDlIoeEf01XIFofuw3dzo3wd7bGsu4b25\nocNj8Z/7P8eWb0nyz86Hu7PNbCmhEvceT9atl+TnZ09gL9L5mZhHqHjdGF/o7mvM7BaG+/yTiEJY\n/RnhByiJ7aN//9FzrdaGzQ4l/CJ5trtv6PDLsZPWh+pxwJ1d1nmQjb9Uhqp4dbGOjZWmmVEbPm5m\n97r7N4bYXqIX3EWn19W+ve8TxrOcRphyZFn0vO/T+ZfTSWwc49Jyfo82vA2YTfdfvGcQxkR9mu5j\nW1rt+DShOtJJ+y/gjxC6qOLO7tHOVwCPjb66teEiwliiTv7eY9vdDPJzk9Q5hOrpK4EPRv8uavsj\napD9PgH4OqF7/91m9gN3f/hYR4OmLyZ0Bb6HMBZrFbAj4Zd+670b5j0E+AZwbTRW7yV0/iOsdU73\nfE3ufqWZvYVwXNr31TGcd3AW8O1on7sSKoS/iMY/xXsB5kXtuT3hdtt1+7nfntC1l8SLCd25WwDP\nJYzTWwW8P8Fzk/z8jBGO25u6bGOUimNcFj8njaIQVn+LCD8o+xEG5vfyaOBe33TgcC+nARe7e/ug\n5F5ac5rd6+4XdFvJzJYQSt3tJf9RjMf3aWY/IXS5HU34hdLNYsIA73Z7xx4f1J4dlj2KaDBsdAXT\nvxK6kT/YWsHMpjGxGtUyoSraIxhvQ+iC/IS7dwzd7r4aeGY0QPmRhPNoWzb9wL8l+nd9r/ezzbXt\n61r3K1enEX7pfMvd/9JlnZuBmQPsP4lFwJPMbErU9Twyd3/AzM4DXmlmXyF0AX6ow36fZmbbJKiG\nnU8YCjCF0C34NTN7Yqyb70hCRWe+u/+m9SQzO6ptO8O8h7j7X83sKkK4XA78Cnh8h9cDCX6O3f0z\nZrYL4aKEY6Nt9grn7W5p+/leSrio4ilsWnl/POEc7HVhy73AChL+3JvZFEJ37C8StvXS2M/deWb2\nBMLnUJIQluTn52bCkIWL2rp9h7GY8LO/F2FMZmufUwhht1OXrySkMWH19yPCGJoPROMcHmax39DR\nX83PJPkP1EsIf4m/o9+KbX5JGGtzsnW+TH4uPDxe5EfA063z1BWjVKEe3kz0taHPej8FHmtm8a7C\nyYS/Mlcy8a/SJF5lsUvFzewphBDc+hBvfXC2v863kM7P7QcJv+Q+2W9Fd7/J3S+MPvgvb3tsCeGc\nOcHaps6Aje/nCN4EzAHe12Od7xOuGH1Oh/1vbmYzhtjvDwhh960dtjnKufctwi/r0wlX8Z3V9nhr\nWoD2cNZpv7919w1RmPh3Qrh4c+zx1nkd/zkfo+01jfgefgPYnzDGccLPkYertH8DHGdmu8cfa389\nZvY0wvv9H+5+bnS+DT2tDBvHfrXGumJhGpgDgHN7PTF6Lb8EnhX9EdJ6/kxC8L2dTSt2L4j2N2wg\nSfI5NIjvE/5geuOEHZltFg11SOpn0b/tPwuvIwwZ+RkyNFXCqu9pZrZHh+W/d/e/uvvNZvYhwlU5\nl5vZDwnB4XGED7j/MLOXE/4Cm0X3Lp12RwNnu/vvB2msuy8zsxOA7wF/NrPvELo/dyL85b6CMNYF\nwrQXRwEXmtkZhHEU27Jx2ozFg+wbGIsG7ELojjyWMI6j03ivuI8TpqT4iZl9FribEEIPAd7i7r3G\nlnRzD+H9+HrUhjdHyz4JYQCwhTm63hn9xXkb4crII5g47mYYRxOualqRwrbeQAhnfzGzrxLGu80h\nnGNPI5xXwzoa+LC7d+u6hlCRfSbwv2Z2NmHw8VTCX+4vJnRBXwlhniZC9+uR7r6gxzbPInQXfszM\nHkcIEtMI5+g5TAxPSf2K8D6/GLjM3W+NP+juCyzMK/U6M9uNjb/gnkioQny000bd/TfRsf9vM/tR\ntN3LCefKt8zsc4Tu+GOATqF02PfwG4Q/lh7q8Zr/E7gMuDqqAN5MGNj+UkL1FzPbltCV+D13/3aP\nbfWyv5m9ko3dkScR/ZxF+3gaGweXL4vWbdmN8PnwSuBH0c/FewmfP5ea2RfYOEXFzsCLPEw9MpPw\n2XoiYQxlr7GucUea2QOEcXrPJvwB1unimGF9h/A5ebqZHUE4f53wB8AxhKEI30+yIXe/JnrfXhf9\nsX4+oWfldYTelW+l2O7mKfryTH0N90X/KSre3rb+qwhTSKwijEW4DHhO9NjPoq/Hd9jPKXSeomIV\nsHPbugvoM0VFbN0nESpM9xPC4K2ED4Wntq23E/BNwofpmmi9rxBNwZF0/4QxMPHjs4wwzur4hO3d\nmXDl131Re/8MvKptndaxSTJFxSsI1Y67omN5AROnWNiecMXafYRfcj8lfIguZtPL3lvnwiBTVFwD\nTOqw7il9jkPrNR7X4fh8jTDOb230ui4ATujw2geZouIeQldj+7pnti3bIjpXb4jOk/sIvyDez6ZT\ngZxGqDLu3et1RutOi96jRdFr+gehgrJvv/ebtikq2h77ZPS8E7o8PkaoeF4XvZalhF/uT+71XhGC\n0t+BC2LLDib8rK+IjuUXCV2Dqb+H0eOHdNn2PoQq3/2E8/0m4FPRY0YIp7cAs3qdFz3eq/jP9njs\nvdqv7T3p9ZnZ+poXe85+hKlRHozafQXwjNjjuxP+QPoUbedpp3OVjVNUtL5WEcZuvZ22KS96fHb0\n/fmJlk0i/HH3p2g//yR8bp0K7DjIMY7OyXcRwnnr3Pg80fQ8+hr+qzXfjYjkwMJUCxcDL3P3RH+J\nSnqiMUy3ufuLCmzDxwgha3t3f6CodjRNVFle4O6ndHl8GiGs7Orui/NrmTSZuiNFpBHMbBZhKotX\nF9iGKdH+z1MAExGFMBFpBA8TbQ59L8lRmNk2hKshn0voZv5UEe1ouPOZOIVL3AbCWKpB5kkUGYlC\nmIhI9vYl/IK/lzC7fa8JliUD7v6RPo+vI1yMIZIbjQkTERERKYDmCRMREREpQCW6I2dt9QjfZocJ\ncwjmYzzMn7flVOXVKmhVdjvNpzm+OkyQPb4+zTkRpcnWL2+bZm2zQoacVcf6cCvOzWYMc99qke7G\nNpvE2LT226MW5/d/vmapu/edrLoSIWybHR7JJ78/0L2JU+Erwxycz9l9i9z3LelZeePGia1X3TvM\nvKoiGy29vO3uX9t0mitZuloS7mQ058ntdzgSGc7mc2cDsMVead7lbjST5u5yW5L1KhHC8qbwVV1X\nXxMufjroMfsAGwOYwpeMQsErRdvsAUsWsfTy3ymISWrKFMAGoRDWRgGs2q5ftBgIIUwBTEah4JUh\nBTFJSasKVlUKYTEKYPWwYeVyVt54jcKXDGWT8KXglR0FMUlJVatgoBAGbAxfoABWdRtWhnkWFcBk\nEApeBVEQkxFsPnd2pQMYKISp+lUTra7H8XW68lGSUXdjSURBTGQQVe+GbGl0CFMAqweN/ZJBqOpV\nTqqGyaCqXgWDBocwBbDq6zT1xPP22Kmo5kiJKXiVnKphMoC6VMGggSFM4aseVP2SJBS+qkXVMEmq\nDlUwaFgIUwCrvn4Tr/5xyQMAPHabrXJrk5SLgldFqRomCdRhMH5cY0KYAlj1Jal+3fZQuI2MQljz\nKHzVg6ph0k2duiFbGhHCFMCqT92P0omCV82oGiZ91KkKBjUPYQpf1afwJZ0ofIk0Sx2rYFDjEKYA\nVn0KYNJO4asZ1CUpndStCgY1DWEKYNXWb/B9L5uNWdrNkYIpeDWMuiSlTd0G48fVKoTp9kPVN2r1\n61m77Zhmc6RACl8iUtduyJbahDBVv6pP3Y8CCl8isqm6VsGgJiFMAaza0gxfV99zHwAHbbv1yNuS\nfCl8iUhc3atgUIMQpgBWbWlXv/6+bBUAB22byuYkBw+HLwUvadlmDw3OF6DeVTCocAhT+Kq2UQbf\nS/Wp6iUivdR5MH5cJUOYAli1aexXcyl8iUg/TeiGbKlcCFMAqy5Vv5pL4UtEBtGEKhhULIQpgFVX\nXtWvaZuNZbp9GYzCl4gMoklVMKhKCBvfgK98UOGrovLsfjx63g6Z70P6U/gSkUG1AlhTqmBQlRCG\nql9VpLFfzaPwJSKjaFIAg4qEsC2nqoupaooKYFfcvRSAJ24/J9f9Np3Cl4iMomndkC2VCGFSLUVW\nwO5ZsTr3fTaZwpeIpKVpVTBQCJMUqfuxORS+RCQtTZkTrBOFMEmFAlgzKHyJSJqa2g3ZohAmI9Hc\nX82g8CUiWWlqFQwUwmQEZax+TZ+sUzpNCl8ikpUmd0O26DeWDKWMAQzgqF22K7oJtaDwJSJZano3\nZItCmAykrOFL0qHwJSJ5aXoVDBTCZABVCGCX3XkvAIfuOLfgllSLwpeI5EVVsI0UwqSvKg2+X7pq\nTdFNqBSFLymlJYuKboFkpIm3JupFIUx6qkL1Swan8CVlN+fJjy+6CZIRBbCNFMKkKwWweno4gCl8\niUiO1A05kUKYTKDwVU8KXyJSFHVDdqYQJpuoegDbcurkoptQOgpfIlIGCmATKYTJw6oewADmP3Lb\noptQGhr3JSJloG7I7hTCpBbhSzZS+JJK05WRtaJuyN4UwhqubgFswR33AM2siCl8SV3oysh6UQDr\nTiGsoao099cg/rlmXdFNKITGfYlI2agbsj+FsAaqW/WryRS+RKSM1A2ZjEJYwyiA1YPCl4iUnQJY\nfwphDVHX7sem0bgvqbUlizQerAbUDZmcQlgDNKn6NWfzqUU3IRMKXyJSBeqGHIxCWM01KYABHLrj\n3KKbkDp1PYpIlSiAJacQVlPqfqw+hS8RqRJ1Qw4utxBmZqcAH2hbfI+7b5dXG5qiadWvuPNv+wcA\nR+1S3dNKXY/SSBoPVmnqhhxO3pWwG4H5se835Lz/2mtyAANYsW590U0YmsKXiFSZAtjg8g5h6939\nHznvsxGaHr6qTl2PIlJV6oYcXt4hbDczuxNYCywE3uPut+TchtpRAKsuha9yeHDlStasW1t0M5rr\nvsUArL//gWLbkdDc2bOZNGms6GaUgrohR5NnCFsIHAfcAGwDnAz81swe7e73ta9sZicAJwDsvNOO\nOTazWhTAqkldj+Xx4MqVHPqh/2LD+HjRTZEvF92AZP7y9c+z3SO2KroZpaEANrzcQpi7/yL+vZld\nCdwCvBr4VIf1zwDOADjogP09jzZWicJXZ9tOn1Z0E3pS+CqfNevWsmF8nB+c9C7mzJxZdHOa6b7F\nPOLxBxTdisTmzlb3G4QqmALYaAqbosLdl5vZdcCjimpDVSmAdffE7ecU3YSu1PVYbnNmzmSb2VsW\n3YzmWbIIZsxgjipLlaJxYOkoLISZ2TRgb+DiotpQRQpg1aPql0hvmpqimlQFG12e84SdBvwEuJ0w\nJux9wHTgW3m1ocoUvpL55eK7ADh63g4Ft0ThS0TqSd2Q6cmzErYT8D1gDnAvcCVwiLvflmMbKkkB\nLLnV68sxuFpdjyIJLFlUdAtkQOqGTFeeA/Nfmte+6kQBrFoUvkQGo67I6tB0FOnTvSNLSuGrWtT1\nKCJNoACWLoWwElIAqxZVv0SGoK7ISlE3ZDYUwkpGAWw0O83cPLd9KXyJjEZdkdWgbsjsKISVhMJX\nOg7aduvM96GuRxFpGgWwbCiElYACWHWo+iWSgiWLVAWrCHVDZkshrGAKYOn66S13AvCs3dK936iq\nXyLSNOqGzJ5CWEEUvrKxfjz924yq+iWSIg3IrxQFsGwphBVAAawaFL5EsqGuyPLTrPj5UAjLmQJY\n+anrUUSaTOPA8qMQlhOFr2pQ9UskQxqQX3oaB5YvhbAcKIDlZ5dZ04d6nsKXiEigAJYfhbCMKYDl\n67HbbDXQ+up6FMmJBuSXnroh86cQlpFW+AIFsLJS9UskX+qKLC91QxZDISwDqn4V59xFfwfgeXvs\n1HUdVb9EcqYqWCUogOVPISxlCmDlpuqXSDFUBSsvTUdRHIWwFCmAlZfCl4jIRBoHViyFsBQofJWb\nAphIgTQtRWlpHFjxFMJGpABWXgpfIiK9KYAVSyFsBApg5bPHljNYfvNtLL3n7rBAAUykOKqClZa6\nIctBIWwICl/ltd2Nfwv/UfgSEelI3ZDloRA2IAWwcmp1Pa4HmLu7TmyRomlailJSACuXsaIbUCUK\nYOUUH/t16aw5XLpG749IGagrspwUwMpDBYMEFL7KSZOuipSUqmClpHFg5ZOoEmZmh5tZIwObAlg5\nbXLlowKYSOmoClYu6oYsp6TB6mJge2BJhm0pHQWw8tG0EyIlpypYaSmAlU/SMWGWaStKSAGsfBTA\nRKpBVbBy0W2JymuQLsadzGxapwfc/faU2lM4ha/yGSR8zdus4ykqInlQFax0NA6s3AYJYb/rsMwA\nByal05xiKYCVz6DVr10VwkQKpSpYeWgcWPkNEsIOBu7NqiFFUwArl2G7Htf4OABTTbOviORKVbBS\nUgArt6QhzIHb3b12A/MVvspl1GknfrvmIQCOnLZlWk0SkYRUBSsPjQOrhqQhrJYD8xXAykUD70Uq\nSlWwUtE4sOpIGsJ2pWZdkQpg5aFJV0WqT1WwctA4sGpJGsJ2AXYx61wQc/dLUmtRDhTAykPVL5GK\nUxWsdBTAqiNpCFtAGBdmsX9bKnN1pMJXeaj6JVIfqoKVg8aBVU/SEDY3+teAW4Ejo38rQwGsPLKs\nfu2+2eapb1NEulAVrDQ0DqyaEoUwd7+v9f+oS/L++LKyUwArhzyqXztvNjWT7YpIZ6qCFU/jwKqr\n1jflVvgqj7zGfq0c3wDAFmOV6CEXqa4lixTASkQBrJqGCWEefZWaAlh55Dn4fuHaZYDmCcvC7Zcu\nTGU7Ox92cCrbkQKpG7I0NA6s2hKFMDNbxsbgNQP4i5k9HMTcfVYGbRuaAlg56MrHauoVtmYcfNBI\n216+8OoJ21coqyZVwYqncWDVl7QS9sZMW5EiBbByUACrhm6Ba9Sw1U37duOhTGGsIlQFKwWNA6uH\npAPzv5V1Q0al8FUOCl/l1il0ZRW4kojvO942BbJyUxWsHBTAyum8m1cmXjfxmDAz2xY4FtgdeJ+7\nLzWzJwN3uXuh01UogJWDAlg5tQevIkNXL612tapjCmIlpCpYKWgcWHkNEsAg+ZiwA4ELCXODPRr4\nBLAUOArYE3j5QHtNkQJY8coUvvbSPGGVCV3dzDj4IAWxElMVrFgaB1ZerQBmWyR/j5JWwk4DTnf3\nD0SD9Ft+BRyfeG8xZvYe4CPAF9x9qDFnCmDFK1MAA9ihofOEVT14tWu1X+PFSkRVsMJpHFh5DRPA\nIHkIOxD4tw7L7wa2HWiPgJkdAvw78JdBnwsKX2VQ1tsOPTS+HoBZY7WeAg+oX/DqRFWxclEVrHgK\nYOUzbACD5CFsFbBVh+V7A0sG2aGZzQa+Qwh17x/kuaAAVgZlq37F/X7tcqC+84QVHbweWrWStevW\n5bpP9nkUK/74Z6Y9+M9MNr902bL+KzWdJmYtnMaBldswAQySh7AfAx8wsxdF37uZzQNOBX444D7P\nAH7g7heZ2UAhTAGsWGWtfjVBPHwVVfF6aNVKnv2x97JhfLyQ/XP+oB81yU0aG2Pq5CmZbb/S1A1Z\nOI0DK6/zbl45dACD5CHs7cDPgXuBLYDLCN2QlwMnJ92Zmf07sAfhKst+654AnADwyG3nKoAVrMzV\nr7oqQ/CKW7tuHRvGx/n6iW9n6xn5z8+84o9/BmDHJzw29W1PnTyF2Vtskfp260JVsOJoHFh5DXol\nZCdJ5wl7CDjUzJ4CPA4YA/7g7hck3ZGZ7QV8FDjM3dcm2OcZhKoZj9v7Ua7wVQxVv/JVdHdjElvP\nmMWcWfn/ZT7niMNZvvBqtpldz67mUlIVrFAKYOU1yjiwuIFGL7v7RcBFQ+7ricAc4Fozay2bBBxu\nZq8Hprv7mk5PHF+/YchdyihU/cpP2apeZaZB+vlSFaxYCmDlk1YAg+TzhL211+Pu/qkEmzkXuLpt\n2TeBvxEqZH2rY5KfqgawfSdXp0upClWvsmldLSk50GD8QmkcWDmlGcBgsHnCVhKuhLS2xxzoG8Lc\n/Z/AJpc3mdkK4H53vzZhOyQMbU1FAAAgAElEQVRjVQ1fLdtOKv/galW9pPTUDVkodUOWU9oBDJKH\nsI8CJwEXEG5ZdE9qLZDSqHoAA3ggmidsqxLOE6bwlY4ZBx+kLskcqApWDAWwckszgEHygfknm9kX\nCTPc32hmnwBO6zaGKyl3nz/K8yUddRp8/6eSzROm4CWVoypY4RTAymfUqSi6SVwucPe7gOPN7HTC\nvSNfb2bvdfdvp94qyU0dql9lpPAlVaYqWDE0IWs5pTEVRTdJB+Y/rm3Ru4BnAZ83sze5+4Gpt0wy\nVafqV5kofEmlqQpWGA3EL6csxoHFJa2EXU0YgN8alB///wFpN0qypepX+hS+8qdxYSmLApiqYPnT\nOLByyjqAQfIQtmtmLZBcKYClS+GrGJqqIhsKYPlTACunPAIYJB+Yf1umrZDMNSV8PWby9Fz20wpf\nCl5SC+qGLJQCWLnkFcBgwBnz25nZFGBV9O0ad6/OTJkN0pQABjBn0uRMt6/wJXWlKlj+NA6svPII\nYJB8YP44YRxYR+4+KbUWSWqaOPh+6YZ1QLphTF2OUmuaGb8Q6oYsp6ymougmaSXsRV2WTwHOTqkt\nkqImVb/irlm3AoAjJ40+T5jCl9SeuiELoQBWTllORdFN0jFhP+y03MymptscSUNTA1haFL6kSVQF\nK4YCWLnkOQ4srnz3dpGhKXyNRuFLGkVVsEJoQtbyKSqAgUJYbSiAjUYD7qVRNCdYITQQv3yKDGCQ\nfGD+eV0eGkuxLTKEJg6+T5PClzSVAli+NA6sfIoOYJC8EnZfj8d078iCqPo10QFTZiRaT+Gr2jRR\n6wjUDZk7BbDyKjKAQfKB+cdn3RAZjAJYZ1uN9T6lFb7qQ7csGoK6IQujAFYueU9F0Y3GhFWMwldv\n92xYC8C2k6ZsslyD7kUCBbB8aSB++RQxFUU3SceE/aXX4+6+fzrNkV4UwPq7fl344WqFMIUvkYi6\nIXOngfjlU4ZxYHFJK2H7ACuBrwEPZdcc6UYBbDAKXyITqQqWH40DK5+yBTBIHsL2Az4BHAt8EPiy\nu2/IrFXyMIWvwax+cNnD/1f4Eono1kS5UgArnzIGMEg4xYS73+juzwFeDBwHXGtmz86yYaIANqhW\n9WvSzJkKYDWmKyMHpG7IQiiAlUdZAxgMODDf3RcAjzezVwKfN7O3AG9z9z9m0bim0txfg4lf8Thp\n/bI+a0sd6MrIwagKlh8NxC+nMgYwSD4w/7MdFv8aeClwFTA5zUY1mapfyXUa97XPpC2Kao5I+agb\nMlcaiF8+ZZmKopuklbBusV79AilSAEum16D76TYp7+aIlJO6IXOlcWDlU6apKLpJOlnrkVk3pMkU\nvpJJcsXjveNhnrC5Y1M6Pi7VpvFgg1EVLB8KYOVT5nFgcQONCTOzzYBHAgbc4e7rMmlVgyiAJZN0\npvvbxtcACmF1pvFgCagbMncKYOVRlQAGyceEzQT+L/BqYCohhK01s7OAd7j7P7NrYn0pgPWn+b5E\nBqRuyFxpIH65VCmAQfJK2NnAHoQQdn207NHAB4CzAE1XMQCFr/4UvkSGpypYPjQQv1yqFsAgeQg7\nCjjU3f8QW3admd0CXJJ+s+pLAaw/3WRbOlm+8Gp1RfajbsjcaBxYuVQxgEHyEHYhMA/4Q9vyXYAL\n0mxQnSmA9abwJTICdUPmRgGsnKoWwKBPCDOzb0T/XQF8y8yeB6yPPfd5wM9i6+Hur8mioVWm8NVb\nml2P+02aPmpzRKonCmCqguVHAaw8yj4XWC/9KmFzY+tNjr5f37ZsFjAjk9bVgAJYd1mM+5pmie7E\nJRWjrsj+FMDyoYH45VKFucB66RnC3P3ZAGb2XeAN7v7N+ONm9hrgX9z9Fdk1sboUwLrLquvxH9E8\nYdtpigppCnVD5kYD8culquPA4pKOCXsu8LEOy38PfD695tSDwld3WV/1+PdonjCFsPrQBK09qBsy\nNxoHVi51CGCQPIT9FviGmb0HuCFatg/wEeDKLBpWVQpg3WngvQxLXZHdKYBlTwGsXOoSwCB5CHs5\noeL1c6A16MaBc4E3ZtCuSlIA60zhSyQD6obMhQJYudQpgEHye0feC7zEzGYDu0aLF2um/EDhqzNN\nuNrfQ6tWsnZdNe7+dd/yh3Lfp7oiu1A3ZK4UwMqhbgEMBrx3JGGqigcJty1akX5zqkcBrDNVv/p7\naNVKnv2x97JhfLzopiQ2aWyMKZMn57pPdUV2pgCWPV0JWT51CmAw/L0jAdY1/d6RCmATFR2+9q/Q\nPGFr161jw/g4Xz/x7Ww9Y1bRzUlkyuTJzNp8i1z2pSpYF5oVPxe6ErJcqjwXWC+6d+QQFL4mKkvX\n45QKzhO29YxZzJlVvw+XNKgK1kbjwHKhcWDlUvW5wHrRvSMHpAA2UdHVr7i7oikqdhib2mdNkYrR\nOLBcKICVSx3HgcXp3pEDUADbVFmqX3F3RZO1KoRVm7oiO1MAy5YCWLnUPYCB7h2ZiMLXRGWqfkk9\nqSsyRt2QuVEAK4cmBDDQvSP7UgDblMKXZE1VsDbqhsyFroQsj6YEMNC9I3tSANuojF2PUl+qgm1K\nASxbuhKyPJoUwED3juxI4WtTqn5JXlQFa6PpKDKncWDl0bQABjneO9LM/gN4HWGAP8B1wIfd/WeJ\nW5sDBbCNqlj9euyk2vaMN4aqYBGNA8ucAlj5NCmAQb73jvw78C7gb9E2Xg2ca2YHuvtfErc4Qwpg\nG1W1+jXJrOgmyJBUBYvROLDMKYCVS10nY+0nt3tHuvuP2xa918zeADwRKDSEPRy+oPEBrIrVr7g7\nxlcD8MixaQW3RIahKthGCmDZUwArhzpPxtrPQPeOdPcHzewOYJW7D33UzGwS8CLCVZW/7bLOCcAJ\nADvNnTPsrvpS9Wujqla/4u4ZDzfDVgirFlXBYtQNmTldCVkeTRwHFpf4Hi9m9m9mdhuwBFhmZovM\n7JWD7MzMHmNmy4E1wJeB57v7NZ3Wdfcz3P0gdz9o61kzB9lNYgpgwe2XLqxFAJNqUxUMdUPmQFdC\nlkfTAxgkDGFmdizwKeBzhAD1BsL4sG+Y2WsH2N+NwAHAIcCXCBPA7jdQi1OiABbEw5cCmBRBVbBN\nKYBlR+PAykMBLEjaHfl24M3u/k0z+wBwgbufYWa3A28FvpZkI+6+FmjV2q82s8cDbwH+bcB2D03h\nK1DlS8pEVTA0HUXGFMDKQwFso6TdkXsBv+mw/FxgtxH3n9tN/hTAAgUwKQtVwSIaB5YpBbDyUADb\nVNJK2D+BzTssnwk8kGQDZvZx4GfAHdHzXg7MB56ZsA0jUQBrRvg6aLNsxg9KdhpfBdM4sFwogBVP\nAWyipCHsF8BhhAlW9wXujJb/a/RYEtsBZ0f/PkiYluLp7v6rxK0dgqafCJoQwKRaVAXbSAEsO7oS\nslwUwDaVdJ6w42P/vyP2/1OT7sjdjxuoZSlQ9at54eu2aJ6wXTRFRSWoCqZxYFnSlZDl0dTJWPsZ\naJ4wM3sKoRLmwPXufnEmrUqBAljzAhjAvdE8YQph5aYqGBoHljGNAyuPJk/G2k+iEGZmOwI/Ag4E\n7ooW72BmVxPm+rqr65ML0PQA1sTwJdXT6CqYxoFlSgGsPDQOrLeklbDPAhuAPdz9VgAz240wxuuz\nwDHZNG8wTQ9foAAm5acqWKAAlg0FsPJQAOsvaQg7CpjfCmAA7n6LmZ0EXJhJywbU9ABW9Xs+SrM0\nvQqmAJYNBbDyUABLZqAxYR2Mp9KKESmAqfrVMsmKboH0snzh1Y0PYJItBbDiKYAllzSEXQh81sxe\n1ro60sx2Bk6n4EpYkwOYql8TPXaS5gkrq8Z3Q2ocWKY0FUU5KIANJmkIOwn4MXCLmd1FuDpyR8Jc\nXydl1Laemhy+QNUvqaZGV8FQAMuKpqIoBwWwwSWdJ+wO4HFmdhSwN2CEKSouyLJx3TQ5gKn61dst\n46sA2G2s0w0eYutdsrDn4+12O7zZ4WFUqoJpHFhWNA6sHBTAhjPQmDB3Px84P77MzGa5+0OptqoH\nBTCFr17uH18PwG7RXVG7ha3tD31C4m3efdlVE7ajUDa4xlbBNA4sMwpg5aIANrik84S92d0/02H5\n0cAZwM5pN6xd028/pACWzOoHw98Dt1x3/cPLBglcnbQ/Px7KFMb6a3QVTOPAMqMAVh6aDX94SSth\nJ0cVrw8BmNkM4DPAS4H3ZdW4lvXLV4T/KHxJm/YK1ZQn7QOMHrx6aW27FcYUxPprZBVMASxzCmDF\n02z4o0kawuYDvzazWcCvgK8RbuJ9oLvfmFHbNqUAJpH24BUPXEvHl+XWju0PfYKCWB9Nn5JCASwb\nuhKyHDQObHRJB+Zfa2aHAxcAbwFOBk5193zmCdtsai67KQuFr4l6Ba+4yYzl0ZwJ7VAQm0jdkJIF\nBbByUABLR+KB+e6+yMwOJQzM3y+7JjWbAtim4uErSRfj3mPTs2xOTwpiEzWyCqZuyMxoKopyUABL\nT9KB+dcQ5gYD2AJ4GfAkM1sG4O77Z9O85tDUExsNGrwGsWzVStasW5fqNgEm778XS676I7//9cXs\ncsjjEj3nvuW5XVScu8ZWwRTAMqOB+OWgAJaupJWwH2TaioZT9St5d2M/i6N5wuZ1mCds2aqVHHPq\nyYyPZ9yLfsm5iVedNDbGlMmTM2xMcRpZBUMBLAsKYOWgAJa+pGPCPph1Q5qq6QEs7arXMtZ3fWzN\nunWMj4/zxde/ja1mzBp5X50sueqPiSthAFMmT2bW5ltk0paiNHYwviZkzYQCWDkogGVj1Bt4y5AU\nvrLrcuxnqxmzmDMrmw+SOU+dz92XXdXYsWFN74aUdCmAlYMCWHYUwgrQ1ABWZPCS/DSuCqZxYJlS\nACsHBbBsKITlqKmD7xW+mqGRVTAFsMxoKopy0Gz42VIIy0kTq19FhK+pOc8TJkErgDWuCoYCWBY0\nFUU5KIBlTyEsB00LYEVWvvYscJ6wpmtcANM4sExoHFg56HZE+Ug6T9hbez3u7p9Kpzn1ovAlTaBu\nSEmLAlg5aCB+frqGMDPbx93/Gn17GrASWAJY26oOKIS1aVIAK1P4umU8fHjsNlavaR/KqpHdkApg\nmVAAKwcFsHz1qoR9xczudveXAB8FTiLcO/J97n5PLq2rqKYEsDKFr5YVbCi6CY2jACajUgArBwWw\n/PUKYc8B7jWzrd39ZDP7IvAR4EYz+wRwmruvyaWVFaHwJU3SyG5IFMDSpgBWDgpgxeh1KdlO0b+r\nAdz9Lnc/HpgffS0ys1dl2roKaUIAu+WShQ8HsO0PfYICWIM1uRtS0qcAViwFsOL0qoSdD7zT3VeY\nWft9WN4FPAv4vJm9yd0PzKyFJdeEub9U+ZJOmhjAVAVLl+YCK54CWLG6hjB33z727dWEAfitQfnx\n/x+QTdPKr+7VryqGr82ZVHQTaq9x3ZAKYJnQXGDFUwArXtJ5wnbNtBUV1JQAVpXw1bKHrorMVOO6\nIRXAMqFxYMVTACuHRCHM3W/LuiFVUffux6qGr7K4+7Krim5C5hoTwCIKYOlSACsPBbDiJZ4x38z2\nB94O7EvojryecIXkNRm1rXTqXP2qS/haFM0TVmRFbLfD6xlSmtoNKelRACsH3Y6oPJLOmP8c4H+B\nS4FfRIsPBf5gZi9w959k1L7SqGsAq+K4r15WaZ6wTKgbUkalAFYOCmDlkrQS9mHgI+7+gfhCM/tQ\n9FhtQ5jCl0igACbDUgArB90PsnyShrA9gbM6LD8LeGd6zSmXugcwha901XU8WKO6IRXAMqMAViwN\nxC+npCFsCXAg0D5I4kCglrcwqmMAU/jKXt3GgzWuGxIFsLRpLrDiKYCVV9IQ9lXCvST3AH5LGJh/\nKGGg/icyalsh6hy+oP4BbLrmCUtdYwLYkkUKYClTACueAli5DTImbDnwNuC/o2V3AR8APptBuwpR\n5wBW9/DVsltBV0XWsSty+cKrGxXAJF2ajLV4CmDll3SeMAc+DXzazGZGy5Zl2bC81S2ANS18lUGd\nuiI1DkxGoYH4xVMAq4bE84QBmNnuwD7R/69391syaVWO6jb5apO6Hju5aXwFAHuOTS+4JdXVqHFg\nCmCpUwArngJYdSSdJ2xr4OvAc4DxjYvtp8Br3P2+jNqXKVW/6mfNw6dnfurYFakAJsNQACueAli1\nJK2EfQ3YAzgMaJVaDga+RBi0/4L0m5atOgUwha/i1aUrslHdkCiApUkBrHgKYNWTNIQ9DfgXd78i\ntuxyM3sdcEH6zcpOnbofm9712M3KVatYumFiReyB5Q+lvq86VcGa1g2pAJYeBbDiKYBVU9IQdi+w\nosPylUCirkgzezehYrYXsAa4Eni3u1+bsA0jU/Wr/lauWsX7/+/HGffO3ZJjY2NMnTw51X3WpQoG\nzQlgkj4FsOIogFVX0hD2IeAzZnasu98JYGY7Ap+MHktiPvBF4HeARc+7wMz2dff7B2r1EOoSwFT9\n6m3qOmfcx/ni69/GVjNmTXx88mRmbl7czb3LqjHTUWgcWOo0F1g5KIBVU9IQ9mZgHrDYzO6Mlu0I\nrAa2MbOTWiu6+/6dNuDuT4t/b2bHAg8CTybje0/WLYApfHW309g0ALaaMYs5s7L9ULr7sqtqUQVr\nzDgwBbDUKYAVTzfkrrakIewHGex7JjAGPJDBtgGFL5F+GjMOTAEsdZqMtXgKYNWXdLLWD2aw79OB\nPwFXdHrQzE4ATgDYfsutBt54HQKYuh4Hd/P4ylz2U6cB+bUPYBEFsPRoIH7xWuPApNoGnaz1KcC+\nhHtHXufuC4bZqZl9inDvyUPdfUOnddz9DOAMgP122tkH2X6dApjC12DW5ThPWNW7Ips0DkwBLD0K\nYMXTQPz6SDpZ647Aj4ADCfeMBNjBzK4Gnu/ud3V98sRtfRp4KXBk2jPu12H6CVW/yq8OVbCmjQOT\ndCiAFU8BrF6SVsI+C2wA9nD3WwHMbDfg7OixY5JsxMxOJwSw+e5+w+DN7U7VLwGYOWMGH3jrO9hq\nxsxM91PlKpjGgckwFMCKpwBWP0lD2FGE4HRra4G73xJdFXlhkg2Y2ReAY4HnAQ+Y2XbRQ8vdffkA\nbZ6g6gFM1a/0jI2NMXvWLCaNjWWy/TpUwUABTAajAFY8BbB6GmhMWAeDDMA5Mfq3PbR9EDhl2AbU\nJYApfKVjS0t3ItZOql4FUwCTQSiAFU8BrL6ShrALgc+a2cvc/Q4AM9uZcIVjokqYu9twTeysLuEL\nFMDS9Eibltm2q14Fa8w4MBTA0qYAVhwFsHpLGsJOAn4M3GJmdxGujtwR+Ev0WK7qEsAUvqqnqlWw\nJo0DUwBLjyZjLZYCWP0lnSfsDuBxZnYUsDfhtkPXu3vuN++ucgBT9St7142H4YWPHpuR6nbrMDt+\nEwKYpEcBrFgKYM2QdIqKVwHnuPv5wPnZNqm7OgQwha9sjTPQlHKJ1KEbsikBTFWwdGg2/GIpgDVH\n0u7IbwK/BJZk2Jau1i5fwe2XLqx0+AIFsCqrahWsEePAFMBSpYH4xVIAa5akISzVQfXDqHIAU/iq\nripXwRoxDkwBLFUKYMVSAGueQaaoeKKZdbzZtrtfklJ7OhqbPj3LzWdCAaw+qloFAwUwSU4BrBwU\nwJplkBD2oy7LHZiUQltqQeGrWI9IcZ6wqlfBah3AIgpg6VAAK955N69UAGugQULYdu5eyJiwqlAA\nK96OKc8TVsUqWFPGgSmApUMBrHgKYM2VNISlf8lZjSh81U9Vp6Ro0jgwGZ0CWPEUwJot6Q32Ch+Y\nX1YKYOVyzfgyrhlfNtI2qtoN2aQApirY6BTAitcaiC/NlXSy1mzuhlxhmnqi3qpYBQMFMElGAax4\nuhJSIPlkrVPdfU30/x2BE4AtgPPc/dIM21dKqn7VV5WrYApgMggFsOIogElLzwqXme1lZtcBK83s\nj2a2L3AV8FZCELvYzJ6XQztLQwGsvloBrGpVsNoPxFcAS5VuR1QsBTCJ61cJOw24G3gn8FLg58Cv\ngNdGj38O+C/g3KwaWBYKX81Q1QBW2yqYAliqFMCKpQAm7fqFsEOAo9z9T2Z2CfAg8EV3Hwcws88B\nV2bcxsIpgFXHHJsy1POq2A1Z+wAWUQBLh+4HWSwFMOmkXwjbGrgLwN2XmdkK4P7Y4w8AMzNqW+E0\n+L56trepQz+3alUwqHkA01xgqdFA/GIpgEk3SQbmt88R1og5w5pU/Vq0oHMxc4/5h+TcktFt8HB6\nTrLks6pUcU6wpgzEl9EpgBVLAUx6SRLCzjazNdH/pwFfNbPW5CbDlx1KrK4BrFvYmtclgMTXr0og\nu96XA/AYS1agrXI3ZG1pHFhqFMCKpQAm/fQLYd9q+/7sDut8O6W2FK6O3Y/twatb4Ookvm5rO1UJ\nY4OoUhWs9uPAFMBSowBWLAUwSaJnCHP34/NqSNHqVP2KB69BQlcv8w4/mMWXLGTRgitrE8Sq1g2p\nACZJKYAVSwFMkhrkBt61VYcAlkXwalenIFbFbkhQAJP+FMCKpQAmg2h0CFP4GlwriNVB1apgCmDS\njwJYsRTAZFCNDWFVDmB5B68q2SbBPGFV7YasMwWw0SmAFUsBTIbRyBBW1QBWlvA17/CDS9sluW2f\necKq1g3ZhHFgCmCjUwArlgKYDKtRIUzhq/7WhZs5MNkm3ha1aveGVACTJBTAiqUAJqNoTAirYgBT\n+BrcDb4C6D5PmAJYSWgy1lQpgBVLAUyG1YgQVrUApvCVvqp1Q0L9A5iqYKPTDbmLdd7NKxXAZCS1\nDmEKXwLV7IZUAJN+FMCKpQAmaahtCKtSAFP4yl6VAlhtKYClpjUOTIqhACZpqWUIq0oAq2r4WnzJ\nwlJeGdlJlaajqPU4MAWw1GggfrEUwCRNtQphVbn3Y1XDVxVsF5uiokrjwBTAJAkFsGK1roQUSUtt\nQljVql8KX9mYG03WWqVxYApgkoQCWLE0FYUk9dzHbJ943VqEsCoEMIWvfKyJ5gmDagSwlloGsIgC\n2OgUwIqlACZJDRLAoAYhrOwBrG5dj2UfD3aTr2Dtg8s4tCLHuu5XQiqAjU4BrFgKYJLUoAEMKhzC\nyh6+QNWvIqx9cFnRTUhMAUz6UQArlgKYJDVMAIOKhrCyBzCFr2LcfdlV8OidmDZ7VtFN6asJU1HI\naBTAiqUAJkkNG8CggiGszAGsbl2P7creFQlUKoDVsgqmgfipUAArlgKYJDVKAIOKhbAqBLA6hq8q\naM0Hdv/6cndHKoBJPwpgxVIAk6RGDWBQkRC2ZvkKbrlkocJXgRbH5mArm/h8YLuMTe2xZrEUwKQf\nBbBiKYBJUmkEMKhICANVv8qgjF2R7fOBzR2bUmRzulIAk34UwIqlACZJpRXAoCIhbPKM6UU3YRNN\nC19lroLBpvOBrfANAEy3SUU1pysFMOlGAaxYCmCSVJoBDCoSwsqi7gPveylrFax9Qta/bggfpgdt\nNrOIJnVU26koFMBSoQBWLAUwSSrtAAYKYYk1rfrVUtYqWFXuC1nbqSgUwFKlAFYMBTBJKosABjCW\nyVa7MLPDzew8M7vTzNzMjstz/8NYtODKxgawlrJVwapyX8jajgNTAEvN5nNnK4AVRAFMksoqgEH+\nlbAZwLXAt6OvUmt6+CpjFUwBrBwUwEanAFYcBTBJKssABjmHMHf/OfBzADM7M899D6rpAaylbFUw\nUAArlG5HlAoFsOIogElSWQcw0JiwCRS+grJWwfoFsF3HpuXUms4UwKQfBbDiKIBJUnkEMChxCDOz\nE4ATALaZvVUu+1QAC1oBrExVsKQD8bcem5xxS7qrewCT0bWuhJT8KYBJUnkFMChxCHP3M4AzAPbc\ncWfPcl8KXxOVMYAl6YZc5usBmGnFnNp1DmCqgo1GU1EURwFMksozgEHOV0eWkQLYpsrWDTnoQPwb\nN6zixg2rsmxSR5oLTHpRACuOApgklXcAgxJXwrLW5IlX+ylLFaxqV0LWjgJYKhTAiqMAJkkVEcAg\n5xBmZjOAPaJvx4CdzewA4H53vz2vdqj61VnZqmBQnQBWuyqYAlgqFMCKowAmSRUVwCD/7siDgD9G\nX5sDH4z+/6G8GqAA1lnZBuMnuRKyaApg0osCWHEUwCSpIgMY5D9P2ALA8txni8JXf2UKYGWnACa9\nKIAVRwFMkio6gEFDxoQpgPVWpm7IUceB7TFp8zSb05ECmPSiAFYcBTBJqgwBDBoQwhTAeitTN2Qa\nA/G3zHhqCgUw6UUBrDgKYJJUWQIY1DiEKXwlV5cABvDPaJ6wLMKYApj0ogBWPAUw6adMAQxqOk+Y\nAlgyZeqGhHSuhFy0YRWLMpgnTAFMelEAK9Z5N69UAJO+yhbAoIYhTAEsmbJ1Q5b9SkhQAJPOFMCK\npQAmSZQxgEGNuiMVvgZXlgBWdrWdDR8FsFEpgBVLAUySKGsAg5pUwhTABrP4koWlCmBlroLVNoAt\nWaQANiIFsGIpgEkSZQ5gUIMQpgA2mLKMA6tKAKslBbCRKYAVSwFMkih7AIMKd0cqfA2uLOPAsgxg\ne6U0T1idB+IrgI1GAaxYCmCSRBUCGFQ0hCmADa/OAQxgZgpTUyiASTcKYMVSAJMkqhLAoILdkQpg\nwynLODDItgvyvvF13De+bujn1zmAyWgUwIqlACZJVCmAQcUqYQpgwynTOLCsx4DdOr4agK3HJg/8\n3LoHMFXBhqcAViwFMEmiagEMKhLC1ixbzqIFVyp8DaFM48CqMAhfAUzaKYAVSwFMkqhiAIMKdUcq\ngA2uTAGszBTApBsFsGIpgEkSVQ1gUJEQNnXG9KKbUFllCWBlrYIpgEk3CmDFat2QW6SXKgcwqEgI\nk8GVYRyYAlhBFMBGpgBWrFYAUxVMeql6AIOKjAmTwZShG7KoALbPpC0SracAJt0ogBVLAUySqEMA\nA4Ww2mlyAAOYbpP6rqMAJt0ogBVLAUySqEsAA3VH1krTAxjAveNruXd8bd/1FMCknQJYsRTAJIk6\nBTBQJax2mhzAAG4bX21jpUMAABEjSURBVAPA3LEpHR+v5Q25FcBGpgBWLAUwSaJuAQxUCauNssyI\nX9ZB+KAAJr0pgBVDAUySqGMAA1XCJrhpwRU9H99z/hNzaklyZbkSUgEsZwpgqdh87mwFsIIogEkS\ndQ1goBAGTAxeux/RuaJ082+u5KYFV5QqiJVlHJgCWM4UwFKhAFYcBTBJos4BDBoewuLhq1vwitv9\niEO4+TdXZtmkgZQlgJVZ60rIWlEAS4UCWHEUwCSJugcwaGgIGzR8lVGZAliZqmD7Tdp4d4VaTkWh\nAJYKBbDiKIBJEk0IYNDAENYKYFUNX3EKYBNNs3CtiQKYdKMAVhwFMEmiKQEMGhTC0ghfZemKLPpK\nyLIGMIB/jK9l9aKbmYMCmEykAFYcBTBJokkBDBoQwtLueix6UL4CWG+3PbgU5s7mcfvtV3RT0qMA\nlgoFsOIogEkSTQtgUPMQlmbXYxmqYEVPRVH2ALZ84dWw5w5Mmz2z6KakRwEsFQpgxVEAkySaGMCg\npiEsq3FfRVbBih6IX4kABgpgMoECWHEUwCSJpgYwqOGM+VkEsKKrYApgvWkQvnSjAFYcBTBJoskB\nDGpUCcv6qseiqmAKYL0pgEk3CmDFUQCTJJoewKAmISzLAFZkFUwBrLdOAexJU2cV1Zx0KIClQgGs\nOApgkoQCWFDpEJbXpKtFVMEUwHrrVgGbahXuYVcAS4UCWHEUwCQJBbCNKhvC8ph09ebfXKkAVkK9\nuiBvXb8agF03m5Zrm0amAJYKBbDiKIBJEgpgm6pk2SCvAFYkBbDO+o0BW7x+NYujIFYZCmCpUAAr\njgKYJKEANlHlQlieAayoKpgCWGcahC/dKIAVRwFMklAA66wy3ZF53/NRAaxcFMCkGwWw4iiASRIK\nYN1VohK2etkKIJ8AVuQ4MAWwzhTApBsFsOIogEkSCmC9VSKEQX4BrAhF344IFMBypQCWCgWw4iiA\nSRIKYP1Vojty6szpme+jqHFgZbgSsk4B7LCpJf+loACWCgWw4iiASRIKYMlUIoRlTQGsfIatgG1m\nlkVz0qEAlgoFsOIogEkSCmDJVaY7MisKYOUzShfkonWrWLRuVdpNGp0CWCoUwIqjACZJKIANptEh\nTAGsfEYdA3bHhjXcsWFNmk0anQJYKhTAiqMAJkkogA0u9xBmZiea2a1mttrMfm9mh+XdhjgFsPLQ\nIHzpRgGsOApgkoQC2HByDWFm9hLgdOCjwGOB3wK/MLOd82wHFDMVhQJYdwpg0o0CWHEUwCQJBbDh\n5V0Jeytwprt/1d3/6u7/CdwNvCHPRiiAlYsCmHSjAFYcBTBJQgFsNLmFMDObAhwI/LrtoV8DT8qr\nHUXMBaYA1l0tA1hEAWw0CmDFUQCTJBTARmfuns+OzHYA7gSOcPdLYsvfD7zC3fdqW/8E4ITo2/2A\na3NpaLXMAZYW3YgS0nHpTMelMx2XznRcOtNx6UzHZVO7uPvcfisVMU9Ye+qzDstw9zOAMwDM7Gp3\nPyiHtlWKjktnOi6d6bh0puPSmY5LZzounem4DCfPMWFLgQ3Adm3LtwHuybEdIiIiIoXLLYS5+1rg\n98BRbQ8dRbhKUkRERKQx8u6O/BRwlpldBVwOvB7YAfhyn+edkXXDKkrHpTMdl850XDrTcelMx6Uz\nHZfOdFyGkNvA/Id3aHYi8E5ge8Jg+7fEB+qLiIiINEHuIUxEREREGn7vSBEREZGilDKEmdkJZnax\nmf3TzNzM5iV83gvN7HozWxP9+/xsW5ovM5tqZp8zs6VmtsLMzjOznfo855ToGMa//pFXm7Mw6P1H\nzeyIaL3VZnaLmb0+r7bmaZDjYmbzO5wXbmZ759nmLJnZ4dHPyJ3RazsuwXMeY2a/MbNV0fPeb2aW\nQ3NzM+hxMbN5Xc6Vo3Nqci7M7N1m9jsze8jM7jWzn5jZfgmeV9tzZphj0pTzJS2lDGHAFoSZ9E9J\n+gQzeyJwDvAd4IDo3/9nZnWaiv0zwAuBlwGHAbOAn5rZpD7Pu5EwBq/1VdlpyAe9/6iZ7Qr8PFrv\nscDHgM+Z2QvzaXE+Rrgv66PZ9Nz4W5btzNkMwrjTNwGr+q1sZrOA8wlT5jweOAl4B+F2a3Uy0HGJ\nOZpNz5WL0m9aoeYDXyTcweUpwHrgAjN7RLcnNOCcmc+AxySm7udLOty9tF/AQYSJXOclWPcc4Py2\nZRcA3yv6daR0LGYDawl3F2gteyQwDjytx/NOAa4tuv0pHoeFwFfblv0N+FiX9U8F/ta27GvAFUW/\nloKPy/zoZ2tO0W3P6fgsB47rs84bgIeAzWPLTibc6cOKfg0FHpd50blyUNHtzfnYzCDMbflsnTMD\nHZNGni/DfpW1EjaMJzLxvpS/Isf7UmbsQGAysdfo7ncAf6X/a9wtKpPfambfN7PdMmxnZoa8/2i3\n8+IgM5ucbguLMeRxabnazO42swvN7MhMGlgdTwQudfd4dehXhGl05hXSonL5XzNbYmaXm9kxRTcm\nBzMJvUUP9FinaedMkmPS0rTzZSh1CmHbMXHm/XuYOEN/VW1H+Auk/d5c/V7jQuA44OnAv0fr/tbM\nts6gjVmbA0xisPe523mxWbS9OhjmuNxN+Cv+hcALCF3WF5rZ4Vk1sgK6nSutx5pqOfB24MXAM4AL\ngXPM7JWFtip7pwN/Aq7osU7Tzpkkx6Sp58tQcpus1cw+DLy3z2pHuvuCEXaT6L6UZZL0uPTaBD1e\no7v/om1/VwK3AK8mTJ5bRYO+z53W77S86hIfF3e/kRC8Wq6ILoB5O9Dkefuacq4k5u5LgU/GFl1t\nZnMI8z2eXUyrsmVmnwIOBQ519w19Vm/EOZP0mDTxfBlFnjPmf4b+b8DtI2z/H1TzvpRJj8shhGrH\nHODe2GPbMMAvTXdfbmbXAY8asJ1lMMz9R7udF+uB+1JtXXHSui/rQuClaTWqgrqdK1D+z5G8LQSO\nL7oRWTCzTxN+Do5091v6rN6Ic2bAY9JJbc+XUeUWwqJ03N6VlqYrCPeh/ERsWenvS5n0uJjZ74F1\nhNf03WjZTsA+DPAazWwasDdw8TDtLZK7r42Ow1HA/4s9dBTwwy5PuwJ4Xtuyo4Cr3X1d+q3M35DH\npZMDCN2UTXUFcKqZTXP31dGyo4C7gMWFtaqcanmumNnphLAx391vSPCU2p8zQxyTTmp5vqSi6CsD\nOn0R/rI4AHg5oaT7jOj7R8TWuZDYlV+EAcjrgXcTQsa7CaHl4KJfT4rH5UuEq26eSpiG4GJC//yk\n2Do3AG+MfX8acASwK3Aw8FPC1Ty7FP16hjwGLyFcJfpaQgA9nTAGYZfo8W8D346tvyuwglBx3Cd6\n3lrghUW/loKPy5sJ4fRRhGkqPhb9rL2g6NeS4jGZEX1uHACsBN4f/X/n6PGPARfG1p9NqGx8H9iP\nMFbuIeBtRb+Wgo/Lq6PP4n2AvQhd1msJt5wr/PWkeFy+EL3fT4l+B7W+ZsTWadQ5M+QxacT5ktox\nLroBXd74U6JfCO1fx8XWWQyc2fa8YwghZC3hqsHa/EKJXt804HOEbrSVwE+AR7at48Apse+/T/ir\nbC0hwP0Q2Lfo1zLicTgxev/XAL8HDo89tgBY0Lb+EcAfovVvBV5f9Gso+rgQxmcsIswTdT9wKfCM\nol9DysdjfpfPkTOjx88EFrc95zGE7v3VhL/cP0DNphoY9LhEv1SvJ/wx8xBwNfDKol9HBsel0zFp\n/zxt1DkzzDFpyvmS1pfuHSkiIiJSgDpNUSEiIiJSGQphIiIiIgVQCBMREREpgEKYiIiISAEUwkRE\nREQKoBAmIiIiUgCFMBEphJk9xsz+ZmYzzWx3M/tb0W0SEcmTQphIRZnZmWbmPb7mFd3GPq4FrgH+\nSbiZ+BeKbU5nZrbczI7Lan0Raa48b+AtIum7ADi2bdn/AX5ZQFsG4mGm6BeY2VxgtbsvK7pNIiJ5\nUiVMpNrWuPs/4l+E21ptwsz2NbOfmdkyM1tiZt8zs+1ij59pZj81s5PN7J6omvNNM9s8ts6CWJVt\ntZlda2YvjD1+ipld26mRZjY/et6c6PvjzGw5gLvf6+7LzOzSaJ2Dur3YqA2fj33/iug1HRZbdriZ\nLYzaeI+ZfdrMpkSP7WBmF5rZg2a2wswuMrND2vZxiJn91cweMrO3RYv3NLPron2dOuL6bmbHxL5/\nbbQs/roWm9nb2573eTNbEPv+4WMYWzbhGCZ979u2c4yZeWw/XSuuiMjQFMJEas7Mtifc2+5a4AmE\nG8DPAM4zs/hnwBGEKtq/AC8E/hU4ddOt8U1ge+BA4LfA2WY2OYU2voBwE+lBnvNc4AzgGHe/NFq2\nI/AL4I+Em9z/G/Aywk2GAYxwr7vHE25ofxVwsZk9Jnr+ZOD/AVdG6+wGbA68CDieUHU80cxeNsz6\nHV7DdOBDhJutj6TTMRzgve/lHMJ7vj3hxu9/j32//ajtFmkyhTCR+nsD8Gd3f5e7/9Xd/wK8ihAa\n4lWnDcDx7n6tu/8KeBfwuigotKyMqm03EG4M/1D0vKFFQebjTAx8vZ7zVOC7wKuitracSLiJ8onR\na/0p8F/AG81sC3e/093Pcvebotf5X8BFhBuaAxwNbAO8yd1vBP6TcEP0U939Knc/FzgLOGHI9du9\ng3Cz498nfe1djke3Y5j0ve/K3VfFqqwPAhvaKq8iMiSFMJH6OxA4POpiXB51Yd0RPbZ7bL2/uHu8\nInMFMKVtnROi568mVEVe6e7jscf3ifbxoJndEHVRWp/2nUgIc98Z4PWcC6wFLm97bB/girY2XRa9\njj3g4S7M+LE4CnhctO4ewF3u/hBAtJ1xYH1se9e3tjXE+g+LqlRvBd7e/tgQuh3DpO/90W3rnJVC\nm0SkDw3MF6m/MeBndP5lf8+A2zoH+CAwFXgJ8D0z29fdl0SP3ww8I9rngcA3CL/0b+60MTPbEngf\nofsz6fiiQwgVp2cCXwGeG99kj+20lp8HLIwtfw+wf8J9x/cz6vofBn7g7n/qn1N7bLj3MUz63l/C\nptW6o4HPDd0oEUlEIUyk/v4AvBi4zd3X9VjvMWY23d1XRN8fQqg2xQPUg+6+CMJAfEKAORz4QfT4\n2tbjwE1mdgJhbFbHEAa8F7jM3X9jyafU+J67f97MfgRcZ2bHunurcnM98GIzG4tVww6Nv47oKsxl\n0WsYI1TB/hCtewuwg5nNcveHosfH2PSzct/Y6xl0/Zb9CePG9k74mnvpdQyTvvcrY+8bZqZuRpEc\nqDtSpP6+AMwGzjGzg81sNzN7qpmdYWYzY+ttBnzDzB5tZkcRxhh9NRbKALYws+3MbGdCV9oYYY6v\nFjOzaWa2hYUrFh9HGBTeyVTg9Wwcj5XU/QDufifwJuB0M9sheuyLwA7AF81sHzN7ZvQ6Pu/uK81s\nFzN7g5ntbWb7Ey402CdaB8LUHvdF29wL+GzUzneZ2ePN7HmEwfZfH3L9lrcCn3b3u3q8zs2iYznN\nzKYBk4Axi670jPQ7hknfexEpgEKYSM1Fv+ifTBir9EvgOsIv5zXRV8tvoscuBn7EpgPWW44nDHz/\nG+HKw9e4+zWxx/cFVhEqTd8jXL341S5N2ww4091vGuG1fYsw5uuM6Ps7gacTqm9/InSHfo9QsYNQ\nETuGcDXjQsK4qKNabXD3NYTK0ZOBqwmVrlXA/xCuqjw72tdZw6wfswz4RJ+X97FoW62v1wOHAb+O\nrdPzGA7w3otIASzMlygiTWZmZwJz3P1ZRbelbKKB6m909zOzWH/AtswHTnH3+WlvW0Typ0qYiEh1\njBOqeSJSAxqYLyJSEe5+CWESXRGpAXVHioiIiBRA3ZEiIiLy/9utYwEAAACAQf7Ww9hTFDGQMACA\ngYQBAAwkDABgIGEAAAMJAwAYBBm1854IOsSMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1cfff431240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(1)\n",
    "n = 1000\n",
    "X = np.random.normal(size=(n, 2))\n",
    "X[:,0]-=X[:,1] * 2 # correlated X for more fun!\n",
    "y = np.dot(X, [-1, 0.001]) + np.random.normal(size=n, scale=3)\n",
    "\n",
    "beta = np.array([2, 5]).astype(float)\n",
    "prev_beta = beta + 1\n",
    "# improve the solution by restricted coordinate-wise Newton descent\n",
    "hessian = np.dot(X.transpose(), X)\n",
    "betas = [beta.copy()]\n",
    "while not (np.abs(prev_beta - beta)<1e-20).all():\n",
    "    prev_beta = beta.copy()\n",
    "    for i in range(len(beta)):\n",
    "        grad = np.dot(np.dot(X,beta) - y, X)\n",
    "        beta[i] = np.maximum(0, beta[i]-grad[i] / hessian[i,i])\n",
    "        betas.append(beta.copy())\n",
    "print(beta.round(3))\n",
    "betas_history = np.vstack(betas)\n",
    "\n",
    "h = 0.01\n",
    "xx, yy = np.meshgrid(np.arange(-1, 3, h), np.arange(0, 6, h))\n",
    "Z = np.mean((y[:, np.newaxis] - np.dot(X, np.vstack([xx.ravel(), yy.ravel()])))**2, axis=0)\n",
    "Z = Z.reshape(xx.shape)\n",
    "cm = plt.cm.RdBu\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.contourf(xx, yy, Z, cmap=cm, alpha=0.5, levels=[8, 8.55, 9, 10, 12, 16, 25, 40, 70])\n",
    "\n",
    "yl = plt.ylim()\n",
    "plt.vlines(0, *yl, linestyle='--', color='grey')\n",
    "plt.ylim(yl)\n",
    "\n",
    "plt.plot(betas_history[:,0], betas_history[:,1], color='k', lw=1.2)\n",
    "plt.xlabel('Первый коэффициент')\n",
    "plt.ylabel('Второй коэффициент')\n",
    "plt.title('Если есть ограничение, спускаемся вдоль него');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}