nmayorov/robust_regression.ipynb

## robust_regression.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Robust nonlinear regression in scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the main applications of nonlinear least squares is nonlinear regression or curve fitting. That is by given pairs $\\left\\{ (t_i, y_i) \\: i = 1, \\ldots, n \\right\\}$ estimate parameters $\\mathbf{x}$ defining a nonlinear function $\\varphi(t; \\mathbf{x})$, assuming the model:\n",
    "\\begin{equation}\n",
    "y_i = \\varphi(t_i; \\mathbf{x}) + \\epsilon_i\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where $\\epsilon_i$ is the measurement (observation) errors. In the least-squares estimation we search $x$ as the solution of the following optimization problem:\n",
    "\\begin{equation}\n",
    "\\frac{1}{2} \\sum_{i=1}^n (\\varphi(t_i; \\mathbf{x}) - y_i)^2 \\rightarrow \\min_x\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Such formulation is intuitive and convinient from mathematical point of view. From the probabilistic point of view the least-squares solution is known to be the [maximum likelihood](https://en.wikipedia.org/wiki/Maximum_likelihood) estimate, provided that all $\\epsilon_i$ are independent and normally distributed random variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So theoretically it is not optimal when $\\epsilon_i$ have distribution other than normal. Although in engineering practice it is usally not important, i.e. if errors behave as some reasonable random variables with zero mean a result of least-squares estimation will be satisfactory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The real problems start when data is contaminated by outliers (completly wrong measurements). In this case the least-squares solution can become significantly biased to avoid very high residuals on outliers. To qualitatively explain why it is happening, let's consider the simplest least-squares problem:\n",
    "\\begin{align}\n",
    "\\frac{1}{2} \\sum_{i=1}^n (a - x_i)^2 \\rightarrow \\min_a\n",
    "\\end{align}\n",
    "And the solution is the mean value:\n",
    "\\begin{align}\n",
    "a^* = \\frac{1}{n} \\sum_{i=1}^n x_i\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now imagine: $a = 1, n=4, x_1 = 0.9, x_2 = 1.05, x_3=0.95, x_4=10 \\text{ (outlier) }$. The solution $a^* = 3.225$, i.e. completely ruined by a single outlier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the well known robust estimators is l1-estimator, in which the sum of absolute values of the residuals is minimized. For demonstration, again consider the simplest problem:\n",
    "\\begin{equation}\n",
    "\\sum_{i=1}^n |a - x_i| \\rightarrow \\min_a\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the well-known solution is the [median](https://en.wikipedia.org/wiki/Median) of $\\{x_i\\}$, such that the small number of outliers won't affect the solution. In our toy problem we have $a^* = 1$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only disadvantage of l1-estimator is that arising optimization problem is hard, as the function is nondifferentiable everywhere, which is particularly troublesome for efficient nonlinear optimization. It means that we are better to stay with differentiable problems, but somehow incorporate robustness in estimation. To accomplish this we introduce a sublinear function $\\rho(z)$ (i.e. its growth should be slower than linear) and formulate a new least-squares-like optimization problem:\n",
    "\\begin{equation}\n",
    "\\frac{1}{2} \\sum_{i=1}^n \\rho \\left((\\varphi(t_i; x) - y_i)^2 \\right) \\rightarrow \\min_x\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Turns out that this problem can be reduced to standard nonlinear least squares by modifying a vector of residuals and Jacobian matrix on each iteration, such that computed gradient and Hessian appxoximation match the ones of the objective function. Refer to the [paper](http://link.springer.com/chapter/10.1007%2F3-540-44480-7_21) for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The implemented choices of $\\rho$ are the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Linear functior which gives a standard least squares: $\\rho(z) = z$.\n",
    "2. [Huber loss](https://en.wikipedia.org/wiki/Huber_loss): $\\rho(z) = \\begin{cases} \n",
    "z & z \\leq 1 \\\\\n",
    "\\sqrt{z} - 1 & z > 1\n",
    "\\end{cases}$\n",
    "3. Smooth approximation to absolute value loss, \"soft l1 loss\": $\\rho(z) = 2 (\\sqrt{1 + z} - 1)$\n",
    "4. Cauchy loss: $\\rho(z) = \\ln(1 + z)$.\n",
    "5. Loss by arctan: $\\rho(z) = \\arctan z$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The functions 2 and 3 are relatively mild and give approximately absolute value loss for large residuals. The last two functions are strongly sublinear and give significant attenuation for outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The loss functions above are written with the assumption that the soft threshold between inliners and outliers is equal to 1.0. To generalize, we introduce the scaling parameter $C$ and evaluate the loss as \n",
    "\\begin{equation}\n",
    "\\hat{\\rho}(r^2) = C^2 \\rho \\left( \\left(\\frac{r}{C} \\right)^2 \\right)\n",
    "\\end{equation}\n",
    "\n",
    "Note that if residuals are small, we have $\\hat{\\rho}(r^2) \\approx \\rho(r^2) \\approx r^2$ for any $\\rho(z)$ defined above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To illustrate the behaviour of the introduced functions we build a plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from matplotlib import rcParams\n",
    "rcParams['figure.figsize'] = (10, 6)\n",
    "rcParams['legend.fontsize'] = 16\n",
    "rcParams['axes.labelsize'] = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r = np.linspace(0, 5, 100)\n",
    "\n",
    "linear = r**2\n",
    "\n",
    "huber = r**2\n",
    "huber[huber > 1] = 2 * r[huber > 1] - 1\n",
    "\n",
    "soft_l1 = 2 * (np.sqrt(1 + r**2) - 1)\n",
    "\n",
    "cauchy = np.log1p(r**2)\n",
    "\n",
    "arctan = np.arctan(r**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1056cce10>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAGECAYAAACGdAwQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd0lNXaxuHfDoHQIl1ASihSFBGUltBEsCBgwQ4eBMux\nIOjBhqKfBjsqoFJEQQTsXUBRQBEJTQERG4Ii0nvvKbO/P97MZCYkISQzeWeS+1prVqbPE+eclZu9\n97O3sdYiIiIiIuEnyu0CRERERCRrCmoiIiIiYUpBTURERCRMKaiJiIiIhCkFNREREZEwpaAmIiIi\nEqbCKqgZY2oZY74zxvxujPnNGHN3+v2JxpiNxpjl6ZeubtcqIiIiEmomnPZRM8ZUA6pZa382xpQF\nlgFXANcCB6y1I1wtUERERKQARbtdgD9r7VZga/r1g8aYlUCN9IeNa4WJiIiIuCCspj79GWPqAOcA\ni9PvGmiMWWGMecMYU961wkREREQKSFgGtfRpz4+Be6y1B4FXgbpAc2ALMNzF8kREREQKRFitUQMw\nxhQHvgC+sta+lMXjdYDp1tqmWTwWXr+MiIiISA6stTku7QqrETVjjAHeAP7wD2nGmOp+T+sJ/Jrd\ne1hrdYnAy+OPP+56Dbro+yuKF313kX3R9xd5lylTLHXrWrZvz93YUlg1EwDtgP8AvxhjlqffNwTo\nZYxpDlhgLXC7S/WJiIiI5MnChXDfffDdd1ClSu5eE1ZBzVo7n6xH+b4q6FpEREREguXff+Hqq2Hy\nZGjSJPevC6upTym6OnXq5HYJkg/6/iKXvrvIpu8vMhw4AJdeCoMHwyWXnNxrw66ZID+MMbYw/T4i\nIiIS2dLS4Ior4LTTYNw4MH6tA8YYbCQ1E4iIiIgUJoMHw+HDMHp0YEjLrbBaoyYiIiJSWIwfD9On\nw6JFULx43t6jSE19mrxEWQk7hel/syIiUjh9+y3ccAMkJUGDBlk/JzdTn0VuRE1/5CObwraIiIS7\nP/+E3r3hww+zD2m5pTVqIiIiIkGycyf06AHPPQfnnZf/9ytyU5+F6fctivQdiohIuDp2DC68ENq2\ndYLaieRm6lNBTSKKvkMREQlH1sJNN8H+/fDxxxCVizlLrVETERERKQDPPgu//grz5uUupOWWgpqI\niIhIPnzwgbOZ7eLFUKZMcN9bzQQRKjExkSi/yB4VFcUTTzzhYkUiIiJFz6JFMGCAs1/aaacF//01\nohbB/LeqWLx4MTVr1nSxGhERkaJl7Vq48kqYNAmaNQvNZyioRTD/RfWtW7d2sZKcHTt2jJiYGLfL\nEBERCZq9e6F7dxgyxPkZKpr6LCSioqIYOnSo77Z3avTvv/+me/fuxMbGUqdOHZ588snjuiZ37NjB\nHXfcQc2aNSlZsiRnnHEG48ePD3jOzp07uf3222nUqBFlypShdu3a3HDDDWzevDnged7P/f3337n4\n4ouJjY3luuuuC90vLiIiUsBSUuDqq+GCC2DgwNB+lkbUCpGsdu3v2bMnN998M/fddx/Tpk3j8ccf\np1atWvTr1w+A/fv30759e44dO8bQoUOpW7cuX3/9NXfeeSfHjh1jwIABAOzevZuYmBiefvppqlat\nypYtW3jxxRdp164df/7553EjZpdffjm33norDz/8cMBaOhERkUhmLfTvDyVLwsiRof88BbVC7v77\n76dv374AdO7cmTlz5vDee+/5gtrLL7/M+vXr+e2336hfv77veXv37mXo0KH079+fqKgoGjZsyCuv\nvOJ737S0NBISEoiLi+Orr77iiiuuCPjce+65h4Gh/meGiIhIARs2DJYudbbhKFYs9J+noY5sGBO6\nS0HqnmnivEmTJqxfv953++uvvyY+Pp46deqQmprqu1x00UXs2rWLP/74w/fcV199lWbNmhEbG0vx\n4sWJi4sDYPXq1cd9bs+ePUP0G4mIiLjjgw9g7Fj44guIjS2Yz9SIWjYKy+b3FStWDLgdExPD0aNH\nfbe3b9/OmjVrKF68+HGvNcawa9cuAEaNGsU999zDfffdx8UXX0yFChVIS0sjPj4+4P28qlevHuTf\nRERExD0LFjjr0WbPhho1Cu5zFdSKuMqVK1OtWjVefvnlLB9v2LAhAO+//z4XXHABL7zwgu+xtWvX\nZvu+Wa2XExERiUR//+00D0yZErptOLKjoFbEde3alVGjRlGrVi2qVKmS7fOOHDlCuXLlAu578803\nQ12eiIiIq3btgm7dIDERunYt+M9XUCviBg0axAcffECHDh0YNGgQDRs25NChQ/z555/Mnz+fzz//\nHHAC3bBhw3j22Wdp1aoVc+bM4ZNPPnG5ehERkdA5ehSuuAJ69oTbb3enBgW1CGWMyXF6MbvHM99/\nyimnsHDhQp544gmGDRvGpk2bKF++PI0bN+aqq67yPe+xxx5j7969jBw5kqNHj9KpUydmzpxJvXr1\nTqouERGRSODxwE03QbVqzoHrbjGZNz+NZMYYm9PvY4w5brNXiSz6DkVEpCA8/LCzBcc330CpUqH5\njPS/aTmObmhETURERMTPa6/BJ5/AwoWhC2m5paAmIiIikm7GDKdxICkJKld2uxoFNREREREAfvoJ\n+vaF6dPh9NPdrsahkwlERESkyFu3Di691Jn2jI93u5oMCmoiIiJSpO3d6+yV9sADcOWVblcTSF2f\nElH0HYqISDAdO+ZsZNusGbz0UsF+dm66PhXUJKLoOxQRkWDxeOCGGyAlxTlwvVixgv18bc8hIiIi\nko2HH4YNG5yD1gs6pOWWgpqIiIgUOaNHw9SpsGCB+3ul5URBTURERIqUzz6DZ55xQlqlSm5XkzN1\nfUaoxMREoqKi8Hg8QXm/qKgo/u///i8o7yUiIhKuFi2C226DadOgbl23qzkxBbUIFuzDz3WYuoiI\nFGarVkHPnjB5MrRs6XY1uaOgFsEiqfvRWktKSorbZYiISBG1ZQtccokz5dmtm9vV5J6CWoT7559/\n6N69O7GxsdSpU4cnn3zSF+AmTZpEVFQU69evD3iNd9o0M4/Hw9NPP03NmjUpXbo05513HitWrDju\neZ9++inx8fGUKVOGChUqcO2117Jhw4aA59SpU4c+ffowceJEGjduTExMDDNmzAjiby4iIpI7+/c7\n4eymm+Dmm92u5uQoqEW4nj17csEFFzB16lSuuOIKHn/8cSZPnnzC12U1zTllyhS+/vprxo4dy6RJ\nk9i2bRtdunRhz549vueMGzeOq6++mrPOOotPPvmE1157jd9++43zzjuPgwcPBrz/d999x0svvcTQ\noUOZOXMmTZs2Dc4vLSIikkvJyXDVVdCmDTz6qNvVnDx1fUa4+++/n759+wLQuXNn5syZw3vvvUe/\nfv1yfF1W06ZHjx5l1qxZlErvU27Tpg0NGjRg5MiRPPHEExw8eJDBgwdz8803M2HCBN/rWrduTaNG\njXjjjTe45557fO+/d+9efvrpJ0499dQg/bYiIiK55/HALbdA6dLOdhyRuBRbQS0bZmjovk37ePDW\nlnXv3j3gdpMmTbKcrsyNbt26+UIaQFxcHPHx8SxatAiARYsWceDAAXr37k1qaqrveTVr1qRRo0bM\nmzfPF9QA4uPjFdJERMQ1Q4bAmjXwzTcQHaGJJ0LLDr1ghqlQqlixYsDtmJgYjh49mqf3qlq16nH3\nnXrqqaxcuRKA7du3A3DBBRdk+fpKfpvRGGOoXr16nuoQERHJr1Gj4PPPnb3SSpd2u5q8U1ArxEqW\nLAlAcnJywP27du3K8vnbtm3L8r4aNWoAGUFs8uTJNGnS5LjnxsbGBtzWdh8iIuKGDz+EYcNg/vzw\n39D2RBTUCrG4uDgAfv31V04//XQAUlNTmTVrVpYhasaMGRw+fJjS6f/0+Pfff1m8eDFDhgwBoG3b\ntsTGxvLXX3/Rp0+fAvotREREcm/OHBgwwJnurFPH7WryT0GtEPI2CrRq1Yr69evzwAMP4PF4KFGi\nBGPHjiU5OTnLZoKSJUty0UUX8cADD3D06FEef/xxypcvz6BBgwA45ZRTeOGFF7jrrrvYsWMHXbt2\npVy5cmzatInvv/+e888/n169egXUICIiUlB+/hmuv94ZUTv7bLerCQ4FtQhljMlyVMz//ujoaKZO\nncpdd91Fv379qFSpEv/73/+Ij4/niSeeOO61ffv2pXTp0gwYMICdO3fSunVrPvzwQ8qXL+97zm23\n3UatWrV44YUXePfdd0lNTaVGjRp07NiRc845J6AOERGRgrJ2LXTvDmPHQqdOblcTPKYwjXwYY2xO\nv48xRiM9EU7foYiIZLZjB7RrB//7H/Tv73Y1uZf+Ny3HkQ1teCsiIiIR6+BBZyTt2msjK6TllkbU\nJKLoOxQREa/kZOjRA2rXhvHjI29D29yMqCmoSUTRdygiIuCcOvCf/8CRI/DRR5G5oW1ugloE/loi\nIiJSlFnrrEfbuBFmzozMkJZbhfhXExERkcLomWfg+++di9/Jh4WSgpqIiIhEjPHj4Y03nKOh/HaP\nKrQU1ERERCQifP45PP64M5JWVI6TVlATERGRsPfdd3DbbfDVV9CggdvVFBztoyYiIiJhbdkyuO46\n52ioFi3crqZgKaiJiIhI2Fq1ytkr7fXXC9fRULmloCYiIiJhacMGuOgip8vziivcrsYdYRXUjDG1\njDHfGWN+N8b8Zoy5O/3+isaY2caY1caYWcaYItDnUXAmTpxIgwYNiImJoUKFCqxbt47ExETWrl17\n0u/VqVMnzj///ID7hgwZwkUXXUSlSpWIiopi8uTJwSpdREQKqZ07nZA2cCDcdJPb1bgnrIIakAIM\nstY2AeKBu4wxZwAPAbOttQ2Bb9NvSxBs3ryZ2267jfbt2/Pdd9/x7bffsnbtWp544ok8BTVjDCbT\nGR6jR4/m2LFjXHrppb7niIiIZMd7fufll8P997tdjbvCquvTWrsV2Jp+/aAxZiVQA7gMOC/9aZOB\nuSisBcVff/2Fx+PhxhtvpG3btgDMnTsXIE9HNVlrjwti+/fvB2DNmjVMmTIlfwWLiEihdvSoM815\n9tnw7LNuV+O+cBtR8zHG1AHOAX4Aqlprt6U/tA2o6lJZYWX16tX07NmTqlWrUqpUKeLi4rj22mtJ\nS0sDYNWqVfTs2ZMKFSpQunRpEhISmDlzpu/1/fr1801TdunShaioKOrWrUvnzp0BuPDCC4mKiiIq\nKop58+blu16d0SkiIjlJTYVevaBiRRg3LvIOWQ+FsBpR8zLGlAU+Ae6x1h7wH6Gx1lpjjP7iA927\nd6dSpUqMGzeOypUrs3HjRr766iustWzevJn27dtTrlw5xowZwymnnMKYMWPo3r07X3zxBV27duWx\nxx6jZcuW3H333YwdO5Zzzz2XUqVKMX/+fO666y5GjRpFq1atADjjjDNc/m1FRKQw83jgv/+Fw4dh\n2jQoVsztisJD2AU1Y0xxnJD2lrX28/S7txljqllrtxpjqgPbs3t9YmKi73qnTp3olNde3lDG+CCM\nLO3cuZM1a9YwcuRIevTo4bu/V69eAIwYMYK9e/fyww8/UK9ePQC6devGmWeeySOPPELXrl2pV68e\njRs3BuDMM8+kdevWAOzatQtwwpn3PhERkVCxFu67z9mKY/ZsiIlxu6LQmDt3rm95UW6FVVAzztDZ\nG8Af1tqX/B6aBvQFhqX//DyLlwOBQS1fwnyarnLlytSrV4/BgwezdetWzjvvPBr4bdU8b948EhIS\nfCENICoqiuuvv54nn3ySgwcPUrZsWTdKFxERCfDUUzBnDsydC2XKuF1N6GQeQBo6dOgJXxNua9Ta\nAf8BzjfGLE+/dAWeAy40xqwGOqffLvJmz55Ny5Ytefjhh2nUqBH169dn3LhxAOzevZvqWRyEVq1a\nNay17Nmzp6DLFREROc7o0TB5MsycCRUquF1N+AmrETVr7XyyD48XFGQtkaBu3bq+PclWrFjB6NGj\n6d+/P3Xq1KFSpUps2bLluNds3boVYwwV9P8GERFx2ZQpMGwYzJsH1aq5XU14CrcRNcmjZs2aMXz4\ncAB+//13zjvvPBYvXsy6det8z0lLS+ODDz7g3HPPzXHaMyZ9ccCRI0dCW7SIiBRZn30GgwfDrFlQ\nt67b1YSvsBpRk9z75ZdfuOeee7j++uupX78+aWlpTJo0ieLFi9O5c2eqVavGpEmTuPDCCxk6dCix\nsbGMHTuWv//+my+//DLH927YsCHR0dG88cYblC9fnpiYGBo3bpzrNW2Zt+H4/vvv2bFjB1u3bgVg\nyZIllC5dGoCrr746D7+9iIhEstmz4fbb4euvQZsK5ExBLUJVr16duLg4RowYwcaNGylZsiRnn302\nX3zxBeeccw4A8+fPZ/Dgwdx5550cO3aMc845hy+//JKLLroo4L0yb1BbqVIlRo8ezbBhw+jUqRMe\nj4fvvvuOjh07nrCurE4mSExM5Pvvv/c9PmbMGMaMGYMxxrfnm4iIFA0LFsANN8Cnn8K557pdTfgz\nhWkTUmOMzen3McZo09UIp+9QRCRyLV8OXbvCW28553gWdel/03LcD0xr1ERERCTk/vwTunWDV19V\nSDsZmvqUXLHW5jhNaYyhmLaRFhGRLKxd64SzYcPgyivdriayaERNcmXo0KGUKFEi24v/ZrsiIiJe\nGzdCly7w8MNw441uVxN5tEZNcmXLli1Z7svmFRMTQ5MmTUJeh75DEZHIsX07dOwIt94K99/vdjXh\nJzdr1BTUJKLoOxQRiQy7d8P550PPnhCs0x0LGwW14x/XH/kIp+9QRCT87d8PF14IHTrACy+AyTGK\nFF0Kasc/rj/yEU7foYhIeDt8GC65BM48E8aOVUjLiYLa8Y/rj3yE03coIhK+jh6Fyy5zzu2cNAmi\n1LKYIwW14x/XH/kIp+9QRCQ8JSfDVVdBqVLw7rsQrQ3ATkgb3oqIiEjIpaZC797OCNo77yikBZOC\nmgRVnTp16NOnj9tliIhIAUlLg3794OBB+PBDKF7c7YoKF2VeCaqsDmUXEZHCyeOBO+6ATZvgyy8h\nJsbtigofBTURERE5adbCPffAH3/AzJlQurTbFRVOmvqMcCtWrKBnz55UrlyZ0qVL07hxY5577jkA\nZs2aRbdu3TjttNMoU6YMTZs2ZcSIEXg8noD3iIqKYujQoQH3/fvvv0RFRTF58uSA+7///nsuvPBC\nypcvT9myZWnevDkTJ04MeI61lvfff58zzjiDsmXL0qpVKxYsWOB7fPjw4ZQsWZKdO3ce97p69erR\nq1evfP93ERGR0LHWOWlg8WKYMQPKlnW7osJLQS2C/fjjjyQkJLB27VpeeuklZsyYwb333sumTZsA\nWLt2LZ07d2bChAnMmDGDvn37kpiYyCOPPHLce2U3Xel//9SpU+nSpQupqam8/vrrTJs2jZtvvpn1\n69f7nmOtJSkpiZEjR/L000/zwQcfkJaWRo8ePdi3bx8AN998M1FRUbz55psBnzVr1iz+/fdf7rzz\nznz/txERkdCw1jm3c84cmDULypVzu6JCzlpbaC7Or5O9Ez0eaTp06GBr165tjxw5csLnejwem5KS\nYp966ilboUKFgMeMMXbo0KEB961du9YaY+zkyZN9r4+Li7OtWrXK8XPi4uJsxYoV7d69e333LV26\n1Bpj7Lvvvuu7r1+/fvb0008PeG3Pnj3tmWeemeP7F7bvUEQk0jz2mLVNm1q7Y4fblUS+9L9pOWYb\nrVHLhpk7N2TvbTt1yvd7HD58mIULF/Lggw9SsmTJLJ+zZcsWEhMT+frrr9myZQupqamAM0q2fft2\nTj311Fx/3qpVq1i/fj1Dhgw54XMTEhIo5/dPrLPOOguADRs2+O7r378/kydP5ttvv6VLly5s2bKF\nL774ghdeeCHXNYmISMF66in46COYOxcqV3a7mqJBQS0bwQhTobRnzx48Hg81a9bM8nGPx8Nll13G\n1q1bSUxMpHHjxpQqVYrPPvuMp59+mqNHj57U5+3atQsg28/zMsZQsWLFgPti0tuA/D+zVatWtGjR\ngnHjxtGlSxcmTJhA8eLF6du370nVJSIiBeP55+Gtt5yQdhL/zpd8UlCLUBUqVCAqKoqNGzdm+fia\nNWtYtmwZb7/9Nr179/bdP3Xq1OOeGxMTQ3JycsB93mDmVTn9n07ZfV5e3Hnnndxxxx1s3ryZCRMm\ncM0111C+fPmgvb+IiATHSy/B66/D999D9epuV1O0qJkgQpUuXZr27dvz9ttvZzk6dvjwYQCi/baH\nTklJ4Z133jmucSAuLo5ff/014L4vv/wy4HbDhg2pU6cOEyZMCNavQK9evShbtiy9evViw4YN3HHH\nHUF7bxERCY5XXnEuc+ZAjRpuV1P0aEQtgr344oucd955JCQkcN9991GjRg3++ecfVqxYwfDhw4mL\ni+ORRx6hWLFiREdHM3LkyCzPyrz++ut56qmneOaZZ2jTpg1JSUm8//77Ac8xxvDSSy9x5ZVX0rlz\nZ+644w4qV67MypUr2bFjB4mJiQAndQ5nqVKl6NevHy+99BJnn3028fHx+f5vIiIiwTNmDIwc6Ux3\n1q7tdjVFk0bUIljLli1ZsGABtWrVYuDAgXTv3p3hw4dTq1Ytihcvzueff061atW48cYbGThwIJ06\ndeKhhx46bkTt4YcfZsCAAYwePZqePXuyatUq3nrrreM+77LLLmP27NkA3HLLLVx++eVMmDCBunXr\n+p5zsqcSXH311QDcfvvtJ/vri4hICL32mrMubc4ciItzu5qiy5zMCEi4M8bYnH6frEaTxF2PPPII\no0aNYvPmzZTNxY6J+g5FREJv/Hh48kn47juoX9/tagqv9L9pOY5waOpTXLF8+XJWrVrFK6+8wu23\n356rkCYiIqE3cSI88YRCWrjQiJq4om7dumzbto2uXbvy1ltvUaZMmVy9Tt+hiEjoTJoEjz7qTHc2\nbOh2NYVfbkbUFNQkoug7FBEJDW9I+/ZbaNTI7WqKhtwENTUTiIiIFHEKaeFLQU1ERKQIU0gLbwpq\nIiIiRZRCWvhTUBMRESmCFNIig4KaiIhIEfPGGwppkUL7qImIiBQhr7/ubGarLTgig4KaiIhIEfHq\nq/Dss85mtqef7nY1khua+pRsff7554wcOdLtMkREJAjGjIHnnlNIizQKapKtzz//nBEjRrhdhoiI\n5NMrr8CLL8LcuToWKtIoqBVBycnJbpcgIiIFZMQIeOklZyStbl23q5GTpaAWof7++2/69OlDvXr1\nKF26NPXr16d///7s3bs34Hn9+vWjVq1aLFq0iLZt21K6dGkefPBBAHbs2EH//v2pVasWJUuWpHbt\n2tx4440kJyfTr18/pkyZwqZNm4iKiiIqKoq66f8PP3bsGIMGDaJp06bExsZSvXp1LrvsMlatWhXw\n2ZMmTSIqKooffviBG264gXLlylGjRg3uuecejh07VjD/oUREirDnnnPWpc2dC3XquF2N5IWaCSLU\nli1bqFmzJiNGjKBSpUr8888/PPPMM3Tr1o2FCxcGPHffvn306tWLBx54gOeee45SpUqxZ88e2rZt\ny969e3n00Uc5++yz2bZtG9OmTSM5OZnHHnuMnTt3smTJEqZPnw5ATEwM4AS1AwcOMGTIEGrUqMGe\nPXsYM2YMCQkJrFy5kqpVqwZ8fp8+fejduzefffYZCxcuJDExkQoVKpCYmFgg/61ERIqiJ56Ad991\nQlqNGm5XI3lmrS00F+fXyd6JHo9kKSkpNikpyRpj7PLly3339+3b1xpj7LRp0wKe/3//93+2WLFi\n9ueff872Pfv27Wtr1qx5ws9OS0uzhw4dsrGxsXbkyJG++998801rjLGJiYkBz+/Ro4dt2LBhbn+1\nAIX5OxQRCQaPx9pHH7X2zDOt3bLF7WokJ+l/03LMNhpRy8ZcMzdk793Jdsr3eyQnJ/Piiy8yZcoU\n1q9fz9GjR32PrV69mubNm/tulyhRgh49egS8ftasWbRu3ZpmzZrl6fM//PBDhg8fzurVq9m3b1/A\nZ2fWvXv3gNtnnXUW33zzTZ4+V0REsmctPPQQfP21M5JWpYrbFUl+KahlIxhhKpQefvhhRo8ezeOP\nP07btm2JjY1lw4YNXHnllQGhDaBKlSoYYwLu27VrF+ecc06ePnv69Olcf/319OvXj6FDh1K5cmWM\nMXTr1u24zwaoWLFiwO2YmBitURMRCTJr4d57Yd48ZzPbSpXcrkiCQUEtQr3//vv07duXIUOG+O7b\nv39/rl9fpUoVNm7cmOfPbtCgARMnTvTdl5KSwq5du/L0fiIikj8eD9x1Fyxf7hwLVb682xVJsKjr\nM0IdOXKE6OjAnP3mm29m+dzMo2kAF110ET/++CO//PJLtp8RExPDkSNHjrv/8OHDFCtWLOC+t956\nC4/Hk5vSRUQkiNLS4Oab4Y8/YPZshbTCRiNqEapr165MnjyZpk2bUr9+fT799FMWLVqU5XOd9YqB\nBg0axLvvvssFF1zAo48+yllnncXOnTuZNm0a48aNo2zZsjRp0oTx48czbtw4WrRoQcmSJWnatCmX\nXHIJU6dO5d5776V79+4sXbqU0aNHU758+Sw/S0REQiMlBfr0gV274KuvoHRptyuSYFNQi1CjRo3C\nWssjjzwCOAv233vvPVq3bh3wPGNMliNq5cqVY8GCBTz66KM899xz7Nq1i6pVq9KlSxdKlCgBwK23\n3srixYsZMmQIe/fupU6dOvzzzz/897//ZcOGDUycOJHXXnuN1q1bM336dHr27HncZ2X12dnVJCIi\nuXfsGFx3HaSmwvTpULKk2xVJKJjCNAJijLE5/T7GGI34RDh9hyIicOQIXHkllCnj7JWW/u9riTDp\nf9NyHLnQGjUREZEIcuAAdO8OFSvC++8rpBV2CmoiIiIRYs8euOgi52D1KVMgWguYCj0FNRERkQiw\nYwd07gzx8fD665Cp+V4KKQU1ERGRMLdpE3TsCJdeCiNGgPqxig4FNRERkTC2dq0T0m66yTloXSGt\naFFQExERCVN//umEtPvugwcfdLsacYOWIYqIiIShn35yujufew769nW7GnGLgpqIiEiYmT/f2Sdt\n3DjnpxS0sUzNAAAgAElEQVRdYRfUjDETge7Admtt0/T7EoFbgR3pT3vYWvt1Ht8/GGWKiIiExMyZ\n8J//OBvZXnih29WI28LuZAJjTAfgIDDFL6g9Dhyw1o44wWtzPJlAREQknH38Mdx1F3z2GbRt63Y1\nEmoReTKBtTYJ2JPFQxoKExGRQuvNN+Huu50RNYU08Qq7oJaDgcaYFcaYN4wx5d0uRkREJFhGjIDE\nRPjuO2je3O1qJJxESlB7FagLNAe2AMPdLUdERCT/rIVHH3VOGkhKgkaN3K5Iwk3YNRNkxVq73Xvd\nGDMBmJ7dcxMTE33XO3XqRKdOnUJZmoiISJ6kpcGAAbBkiRPSqlRxuyIJtblz5zJ37tyTek3YNRMA\nGGPqANP9mgmqW2u3pF8fBLSy1vbO4nVqJhARkbCXnAw33gjbtsHUqXDKKW5XJG7ITTNB2I2oGWPe\nA84DKhtjNgCPA52MMc0BC6wFbnexRBERkTw7fBiuugpKlICvvoKSJd2uSMJZWI6o5ZVG1EREJJzt\n2QM9ekCDBjBhAkSH3XCJFKSI3J5DRESkMNq0yTm3MyEBJk5USJPcUVATEREJsdWroX176NMHXngB\novTXV3JJeV5ERCSEli1zpjuffhpuvtntaiTSKKiJiIiEyJw5cP31MH48XH6529VIJNLgq4iISAh8\n/LET0j76SCFN8k4jaiIiIkE2Zgw8+yzMng3NmrldjUQyBTUREZEgsRb+7//gww+d0wbq1nW7Iol0\nCmoiIiJBkJoKt98Ov/4KCxboSCgJDgU1ERGRfDp8GK67DlJSnAaCsmXdrkgKCzUTiIiI5MOuXXDB\nBVC+PEyfrpAmwaWgJiIikkf//gvt2jmb2U6eDMWLu12RFDYKaiIiInmwfLkT0vr3h+ef12kDEhp5\nWqNmjCkNdAaaA9WAYsAO4B9glrV2c9AqFBERCTPffAO9e8PYsXD11W5XI4XZSeV/Y0xtY8x44Cfg\nWpyA9iewHDgMNAM+M8YsMMZ0DnaxIiIibnvnHbjhBmdDW4U0CTVjrc3dE425FrgMeNlau+QEz60M\n/A+oAtxjrT2a30Jzwxhjc/v7iIiInAxrnSnOMWPgq6+gSRO3K5JIZ4zBWmtyfE5ugo0xpgdQwVr7\n1kkWUB+401p7/8m8Lq8U1EREJBTS0mDgQJg/H2bMgJo13a5ICoNgBrUS1trkPBaR59fm4bMU1ERE\nJKgOH4ZeveDQIfjkEyhXzu2KpLDITVDL1Rq1/AStggppIiIiwbZ9O5x/vrNH2owZCmlS8NRMLCIi\nkoXVq6FtW7j4Ypg0CUqUcLsiKYp0hJSIiEgmCxfClVfCU0/Brbe6XY0UZQpqIiIifj78EAYMgClT\noGtXt6uRok5BTUREBGf7jRdegFGjYPZsaNbM7YpEFNRERERITXVG0RYtci7afkPCRY7NBMaYZsaY\nFcaYfcaY940xVdLvv8EY81XBlCgiIhI6Bw7ApZfCunWQlKSQJuHlRF2ficBjQFtgDvCeMaa6tfYd\noGWIaxMREQmpjRuhQweoXRumT4dTTnG7IpFAJwpqX1hrp1prf7fWvg5cDdxnjKlWALWJiIiEzLJl\nEB/vnNs5bhxEazGQhKETBTVrjGlqjBltjClvrd0LDAauAkqGvjwREZHgmzrV6eh85RV44AEwOe4N\nL+KeHP/9YK2daIzpCqwGDqbflwaMMcZsK4D6REREgsZaGDkShg93Thpo1crtikRylquzPrN8oTFP\nAJcAV1tr1wW1qjzSWZ8iIpKd1FTnYPUFC+CLL5x1aSJuys1Zn/mZkd8LzACUjEREJKzt3QvXXQdR\nUTB/vpoGJHLk56zPPcAwa+36YBUjIiISbP/845zZ2aiROjsl8uQnqC0BvjTG3GiM0a4zIiISdpKS\noF07ZzPbV15RZ6dEnvysUXsH2Ai0A1oB64FvgU+ttbOCVuHJ1aQ1aiIiAsDkyU5H51tvwcUXu12N\nyPFys0YtP0FtCDDNWvubMaYM0AHoDDS01l6RpzfNJwU1ERHxeODRR+GDD5ypzjPPdLsikayFNKil\nf8BFQGVr7bt5fpMgUlATESnaDh6EPn1g50749FOoUsXtikSyl5uglqs1asaYylndb62ddaKQZow5\nNTefISIikh/r1jnr0SpWhG+/VUiTwiG3zQR1jTE3n+ybG2PigTtP9nUiIiInY+FCSEiAfv1gwgQo\nUcLtikSCI1f9L9baJcYYjDEfA+8Dn6WfUJAlY0wz4G5gnbV2aHBKFREROd6UKXD//U7zwCWXuF2N\nSHCd1Bo1Y0xJ4B7gBmA3ztFSe4BUoCJQFWiKs3XH49bav4Jd8Anq0xo1EZEiIi0NhgyBTz6BadPU\nNCCRJ9Rdn42Bc3DCWQlgO7AWWGStTc7Tm+aTgpqISNGwbx/07g1HjsBHH0GlSm5XJHLyCqLr0wCV\ngN3WWk+e3yhIFNRERAq/1avhssvgwgthxAgoXtztikTyJmhdn9m8+Xk4m9xuB3YaY0YaY8rm9f1E\nREROZOZM6NAB7rsPRo1SSJPCLz+HadwI9AAM0BK4DlhsjOlgrd0TjOJEREQArIWRI+HFF501ae3b\nu12RSMHIzxq1u621r2S6ryvQzVp7dzCKy0NNmvoUESlkjhyBO+6AX36BqVOhdm23KxIJjpBOfQK1\njDFN/O+w1n6NMx0qIiKSbxs2QMeOkJICCxYopEnRk5+g9jzwtjFmhDGmhd/9WjEgIiL5lpQEbdrA\ntdfCO+9A6dJuVyRS8PIc1Ky1O4DzgWLAN8aYXcaYlUB5Y0z9YBUoIiJFi7Xw6qtw9dXw5pvwwANg\ncpwcEim88rU9h+9NjIkG2gAXAF3Sr28C3rfWDsn3B+S+Dq1RExGJYMeOwYABsGgRfP45nH662xWJ\nhE7I91HL4YPLAB2BM621w4P+Adl/roKaiEiE2rjRGUWrUQMmTYLYWLcrEgkt14KaWxTUREQi07x5\ncP31MHAgPPSQpjqlaMhNUMvPPmoiIiL5Yq2zce3TTzuHq198sdsViYQXBTUREXHFkSNw++2wYoWz\nJq1ePbcrEgk/+dmeQ0REJE/WroV27SA1FRYuVEgTyY6CmoiIFKivvoL4eOjb19kfrUwZtysSCV+a\n+hQRkQLh8cBTT8Frr8HHHzuHq4tIzhTUREQk5PbsgT59YN8+WLoUqld3uyKRyKCpTxERCamff4aW\nLaFBA5gzRyFN5GQoqImISMhMnAgXXuhsvzFyJBTXadAiJyXspj6NMROB7sB2a23T9PsqAh8AccC/\nwLXW2r2uFSkiIjk6ciTjKKh58+CMM9yuSCQyheOI2ptA10z3PQTMttY2BL5Nvy0iImFozRpo29YJ\naz/+qJAmkh9hF9SstUnAnkx3XwZMTr8+GbiiQIsSEZFcmToVEhLgllucrTfKlnW7IpHIFnZTn9mo\naq3dln59G1DVzWJERCRQSgoMGQIffgjTp0ObNm5XJFI4REpQ87HWWmOMTl4XEQkTGzfCdddB+fLw\n009QqZLbFYkUHpES1LYZY6pZa7caY6oD27N7YmJiou96p06d6NSpU+irExEpor7+Gvr1g//9Dx58\nEKLCbkGNSPiYO3cuc+fOPanXGGvDb3DKGFMHmO7X9fk8sMtaO8wY8xBQ3lp7XEOBMcaG4+8jIlLY\npKZCYiJMmuSsRTvvPLcrEok8xhistSbH54RbsDHGvAecB1TGWY/2GDAV+BCoTQ7bcyioiYiE3ubN\n0Ls3REc7Ia2qVg2L5ElEBrX8UFATEQmtmTOdqc4774RHHoFixdyuSCRy5SaoRcoaNRERcVFqKjz2\nGEyZAu+9B1r+K1IwFNRERCRHGzdCr15QurTT1XnqqW5XJFJ0qD9HRESy9cUXzoHq3brBV18ppIkU\nNI2oiYjIcY4dg4cegk8/hY8/hvbt3a5IpGhSUBMRkQB//QXXXw9xcbB8OVSs6HZFIkWXpj5FRMTn\n7bedA9VvuQU++UQhTcRtGlETEREOHICBA2HxYvjmG2jWzO2KRAQ0oiYiUuQtWQLnnutsYLtsmUKa\nSDjRiJqISBHl8cCLLzqXMWPgmmvcrkikaDiWeozlW5fn6rkKaiIiRdDmzXDjjXD0qDOiFhfndkUi\nhdeBYwdYtHER89bNI2l9Ess2L6NR5Ua5eq2OkBIRKWKmT4f//jfjGKho/ZNdJKh2Ht7J/PXzfcFs\n5Y6VnFv9XDrGdaRD7Q4k1ErglJhTdNaniIhkOHQI7rvPOa/z7behXTu3KxIpHDbs20DS+iSS1iUx\nb/08Nu7fSNtabelQuwMd4zrS8rSWlIwuedzrdNaniIgATpPADTdAq1bw889QrpzbFYlEJmstq3et\n9o2WJa1P4mDyQTrU7kCH2h34b4v/cnbVs4mOCk7E0oiaiEghlpbmNAsMHw4vv+yc2SkiuZfmSWPF\nthUkrUvyBbOS0SV905gdanegceXGGJPjwFiWNPUpIlKErV8Pffs63Z1TpqhhQCQ3jqUeY8nmJb5p\nzEUbFnFa7Gm+acwOcR2oXa52UD5LQU1EpAiyFt55BwYNgnvvhQcfhGLF3K5KJDx5OzK9I2ZLNy+l\nceXGvmDWvnZ7qpSpEpLPVlATESlidu+GO+6A3393GgbOOcftikTCi7cj0zti5u3I9AYzb0dmQVBQ\nExEpQmbOdM7ovOYaeOYZKFXK7YpE3Ldh34aAhf8b928koWaCL5i1qtEqy47MgqCgJiJSBBw+DIMH\nw9Sp8Oab0KWL2xWJuCNzR+a8dfM4lHIoYOF/s2rNgtaRmV/ankNEpJBbvNg5YaB1a1ixAipUcLsi\nkYKT5knjl22/+EJZ0vokSkWXokOcE8oebv9wnjsyQ+bQIfjhB1i0KFdPV1ATEYlAyckwdCi88QaM\nHg1XX+12RSKhdyz1GEs3L/WFsoUbFvo6Mns27smIi0cErSMzaHbsgPnzMy6//QbNmkH79rl6uaY+\nRUQizC+/OKNotWrB+PFQrZrbFYmEhn9H5rz181i2eVmBdWTmibWwZk1gMNu6FRISoEMHJ5y1auVb\nQKo1aiIihUhqasbmtcOGwU03QTjN6IjkV3ZnZLrRkZkrqanOmoP58yEpyflZrJgTyrzB7Kyzst0f\nR0FNRKSQ+PNP6NcPypRxpjvr1HG7IpH88+/InLduHpsObCKhZoJv8b+bHZlZ8q4v846WLV7sDG17\nQ1m7ds7/OXP5LygFNRGRCJeWBi+9BM8+C08+CbffDlFRblclcvK8HZneUObtyPSOlnWM6xjUMzKD\nYvt2WLAgY8Ts99+heXMnkHXoAG3bQqVKeX57BTURkQj211/OKFrx4jBxItSr53ZFIrnn7cj038PM\n/4zMjnEdaVSpUfh0ZHrXl3mnMOfPh23bMtaXdegALVsGdYNCBTURkQiUlgavvAJPPw2PPQYDBmgU\nTcKf/xmZ3o7M6rHV6Vi7o2+7jLjyYXTgbGoq/Pxz4ML/4sUzpjHbt4cmTUJ6/pqCmohIhPnzT7j5\nZufvxRtvwOmnu12RSNa8HZneETP/jkxvMAurjsyDB501Zd5Q9sMPznoybyhr3x7iCjZIKqiJiESI\n1FQYMQKefx4SE6F/f42iSXjJ7oxM71Rm2HVkbt0aOFq2ciWce27G+rKEBKhY0dUSFdRERCLAH384\nW22ULQsTJkDdum5XJOJ0ZPrv+O89IzMsOzKthVWrAoPZ7t3OYn/vVGaLFlAyTOpNp6AmIhLGkpPh\nuedg1KiMjs5wWVctRUvmMzKT1idxMPmg73zMjnEdw+qMTJKTYfnyjG7MBQucvWvat88IZmecEfbD\n0gpqIiJhaskSZy1aXBy8+qqzFZNIQUnzpLFi2wrfwv/MHZkdancIrzMy9+93zsb0jpYtWeIs4PSu\nLWvXLiL/T6SgJiISZg4fdjo5337bWZPWq5dG0ST0sjoj09uR2THO6coMqzMyN20KnMb86y9nawz/\n9WXlyrldZb4pqImIhJE5c+C226BNG2cT2yph1BAnhYv/GZlJ65NYunkpjSs39o2YhdUZmR6Ps9Df\nP5gdOJARytq3d5oASpRwu9KgU1ATEQkDu3bB/ffDt9/C2LHQo4fbFUlhE1EdmUePwtKlGTv+L1jg\n7O7vncJs3x4aNSoSQ80KaiIiLrIW3n8f7r0XrrnG2cA2NtbtqqQwyK4j07vwP6w6MnfvhoULM0bL\nfv7ZWejvDWbt2kH16m5X6QoFNRERl6xbB3feCRs2wPjxEB/vdkUSqay1rNq1yjeN6T0j03/hf9h0\nZFoL//6bEcoWLID16535fm8wi4939qIRBTURkYKWmgovv+wcoj5oEDzwQKFcWiMh5N+ROW/9POav\nn0+p6FK+3f7D6ozM1FT45ZeMUDZ/vrPmzLtNRrt20KwZRIdBiAxDCmoiIgXoxx+dvdAqVXK23GjQ\nwO2KJBJkdUbmabGn+UJZWHVkHjrkHL3kHTFbvBhq1gw8hqlu3SKxviwYFNRERArAvn3wyCPw8ccw\nfDj07q2/U5K97M7IDMuOzG3bMkbK5s+H33+H5s0zFv23bQuVK7tdZcRSUBMRCSFrnXA2aBBccgkM\nG+b60YEShrwdmd5gFrYdmd5jmPyD2c6dGQv+27WDVq2gVCm3Ky00FNRERELk779hwADYuNGZ5uzQ\nwe2KJFxs2Lch4CimsD0j89gx+OmnjPVlCxY4i/y9o2Xt28OZZ4b9MUyRTEFNRCTIjh6F55+HV16B\nwYPhf/+D4sXdrkrckvmMTG9Hpm99WTh1ZO7ZE3gM07Jl0LBh4DFMNWu6XWWRoqAmIhJEs2fDXXfB\nWWc5JwvUDpP13VJw0jxp/LLtl4ARs5LRJQOCWVickWmts0eM/6aya9dC69YZI2bx8XBKGEy5FmEK\naiIiQbBxo7Np7ZIlMGqUThYoSvw7Muetn8eiDYt8Z2R6t8uIKx/ndpnONhm//hq4TUZaWkYoa9fO\naQLQ8G9YUVATEcmH5GQYORJeeMEZSXvoIa2jLuyy68j0jpiFTUfmwYPONhneUPbDD1CjRsb5mO3a\nQb16aj8OcwpqIiJ59M03TrPA6ac7G9jWr+92RRIKEdORuWVLYDfmypVwzjkZo2Vt2zob+ElEUVAT\nETlJ69c7B6gvXeoEtEsvdbsiCaaszshsW6ut7yimsOjI9HicIOa/vmzv3owtMtq3hxYtoGQYdI5K\nviioiYjk0tGjzhTnyy87I2mDB2uaM9Jl7shMWp/EweSDvtGyjnEdObvq2e53ZB496iyA9AazhQuh\nQoXAbszGjbVNRiGkoCYicgLWwtSpTrPAuefCiy9CnTpuVyV5kdUZmSWjSwYEs7A4I3PnTieMeUfL\nfv7Z2a/Mf+F/9eru1igFQkFNRCQHK1c6+6Bt3OiMpF1wgdsVycnIqiPTe0Zm2HRkWuvsjuw/jbl5\ns7M1hjeUtWkDZcq4W6e4QkFNRCQLe/ZAYiK8+65zRuddd2nXgkjg7cj0BrOw7MhMTnZ2+/fu9L9g\nAcTEZISy9u2djfiKFXO3TgkLCmoiIn5SU+H112HoULjySnjiCagSBjstSNa8HZneYBaWHZl79jjT\nmN5QtmyZ0yrsDWbt2mlnZMmWgpqISLpvv3WmOStXdk4VaNbM7YokM29HpjeYec/I9I6Yud6Raa2z\nu783lM2f7+z+793tv107Z0qzXDn3apSIoqAmIkXeqlXwwAPw229Oo0DPntoDNBx4OzL9t8rw78gM\nizMyU1JgxYrAQ8shcNF/s2aaN5c8U1ATkSJr1y5navOdd5ytNgYO1LZTbvKekekfzMKuI3PfPufQ\ncm8oW7LEaQH2n8asU0dJX4JGQU1EipzkZBg7Fp55Bq65xmka0Dq0guffkZm0PomFGxaG1xmZ3kPL\n/UfL/vkHWrbMCGYJCVC+vHs1SqGnoCYiRYa18OmnznmcDRo405xnnul2VUVH2HdkpqY605j+22R4\nDy33TmWec46mMaVAFbqgZoz5F9gPpAEp1trWmR5XUBMpghYudI59OnzYOV3gwgvdrqjwy+6MTG8w\nc70jc98+WLw4cBqzdu2MYKZDyyUMFMagthZoYa3dnc3jCmoiRcjffzsjaD/8AE89Bf/5j7anCpUN\n+zb4Qtm8dfPYdGATCTUTfGvMXO3I9E5j+u9dtmaNM43pDWUJCc6xTCJhpLAGtZbW2l3ZPK6gJlIE\nbN8OTz4J770H993nbLuhczmDJ/MZmfPWzeNQyiHfaJnrZ2SmpjrHLvnvX5aaGrjov3lzKFHCnfpE\ncqkwBrV/gH04U5+vWWvHZ3pcQU2kEDtwAEaMgFdecUbPHn1UjQLB4H9Gpvfw8pLRJX3BrEPtDjSu\n3Ni9jszsujHbtYO2bTWNKRGrMAa16tbaLcaYKsBsYKC1NsnvcQU1kUIoJQXGj3dG0Tp3dn7Wq+d2\nVZErrM/ItBb+/TdwGnPt2uOnMdWNKYVAboKaizsJnjxr7Zb0nzuMMZ8BrYEk/+ckJib6rnfq1IlO\nnToVYIUiEkweD3zwATz2mBPMZsxwGvPk5Hg7Mr1Tmcs2L6NR5UZ0rN2R2869jSlXTHGvIzM5GZYv\nD5zGNCYjlN18szONqW5MKQTmzp3L3LlzT+o1ETOiZowpDRSz1h4wxpQBZgFDrbWz/J6jETWRQsBa\nJ5Q98ohznvUzz0CXLm5XFTkyd2T+seMPWlRvER4dmbt3O6HMG8yWLYP69QO7MePiNI0pRUKhmvo0\nxtQFPku/GQ28Y619NtNzFNREIlxSEjz8sHPW9VNPwRVX6G/2iYRtR6a18NdfTiDzBrONG52zMb1r\ny3Q2phRhhSqo5YaCmkjkWrbMaQ74808YOhRuuEFbbWTlRB2Zrp6RefSo80V6pzAXLoTSpZ1Q5g1m\nZ58N0RG16kYkZBTURCTs/fqrswbtxx+dqc5bbnGmO8WRXUem/+HlrnVkbtsWOI25YgWccUZgN2bN\nmgVfl0iEUFATkbC1apVzDuecOc6h6Xfeqb3QIOeOzI5xzjmZtcvVLvjCPB7444/Aacxdu5wOTG8o\na90aypQp+NpEIpSCmoiEnb//dtaeffklDBoEd98NZcu6XZV7sjsj0zti5toZmQcPOsOc3lC2eDFU\nrpwxWta2rXOYalRUwdcmUkgoqIlI2FizxgloX3wBAwc6Aa0oboXl7cj0BjPvGZneYOZKR6a1sGFD\nRihbuNAZ8mzePCOUtW0Lp55asHWJFHIKaiLiun/+cQLatGkwYIBz3FNRCmgb9m3wLfpPWp/Exv0b\nSaiZ4JvKdKUjMyXFOYLJG8oWLnSOYPJOYbZtC+eeq8WCIiGmoCYirvnrL2f/M/+AVtjPxM6uI9N/\n4b8rHZk7dzpHMHlD2bJlzg7C/ov+69bVPigiBUxBTUQK3MqV8PTTMHNm4Z/i9O/InLd+HvPXzw/o\nyOwY15FGlRoVbEemx+PsceINZQsXwpYt0KZNxhRmmzbau0wkDCioiUiB+fVXJ6DNmeM0Cdx1F5zi\n0ub3oRKWHZkHDzqHlHtD2aJFztClN5QlJEDTptqUTiQMKaiJSMj98IMzxfnjj05A69+/8HRxZteR\n6Q1mBd6RaS2sWxc4WuZd9J+QkHFgebVqBVeTiOSZgpqIhIS18N13TkD76y948EHn7OxI3wct8xmZ\n3o5M187IPHYs48By78XajFCmRf8iEU1BTUSCyuOB6dPhueecs7Ufesg56qlECbcry5vszsh0rSNz\n69aMRf+LFjmdmQ0bBk5j1qmjRf8ihYSCmogERXIyvPsuPP+8M2o2eDBcdVVkLXvK3JGZtD6Jg8kH\nfaGsY1xHzq56dsF1ZKamOgv7/Lsx9+51wph3t//WrQvPPLKIHEdBTUTy5eBBGD8eRo50BnYeegi6\ndImMAZ2szsgsFV2KDnEd3OnI3L3b2d3fO1q2ZIlzDqY3lCUkQOPG2ulfpAhRUBORPNm6FUaNgtdf\nh06dnBG0li3dripnYdWR6fE4+5T4T2Nu2gStWmUEs/h4qFixYOoRkbCkoCYiJ2XlShg+HD75BHr3\ndro4Tz/d7aqy5u3I9E5lutqRuW+f0/66aJFzWbwYKlUKXFt21lkQXcAb3YpIWFNQE5ETsha+/94J\naD/+6Ox/1r+/c/52OMnujMwC78j0eGD16oxQtmgRrF0LLVpkrC9LSNC5mCJyQgpqIpKt5GR4/31n\n/dmRI87o2Y03hs8WG94zMr3BbOP+jbSt1dZ3FFOBdWQeOOAkWO805uLFzq7+/qGsWTMoXjz0tYhI\noaKgJiLH2bkTXnsNxoyBJk2cgNa1q7tr2LPryPQ/iqlAOjKtdTaG8x8tW7MmY0NZ76V69dDWISJF\ngoKaiPj88ovTIPDxx9CzpxPQmjZ1p5Y0Txq/bPslIJi5ckbmwYMZo2XetWVlyjgL/b1ry5o3j9yN\n4kTEVWnWcigtjYPZXK6rWlVBTaQoS0uDadPglVecZVV33gm33Vbwy6ey6sisHludjrU7FlxHprXO\nf4TFi7MeLYuPd36edlpo6xCRsJRm7XFB6kBqauDtzI/nEMIOpqVx1OOhdFQUsdHRlC1WLOBSJiqK\nT5o2VVATKYp27YKJE53pzdNOg7vvhiuvLLiBobDoyNy/P2O0bPFi5xIbmzF9GR+v0TKRCJbs8XAg\nPUwd8AtOWd130O/+zAHLe/2Yx0OZ9BAV6/ezTPrPzGErNlPwKuv3XO/1UlFRROUwM6CpT5EiZulS\nJ5x99hlcdhkMGOBsbh9q2XVkeqcyQ96R6fHAn39mBDJvJ+a552aMlMXHa22ZiIu8I1beEOW97Pe/\nnemxnG5bcAKUX5CKzSJoeQNWwH3+gSv98dJRUQW3AXY6BTWRIuDoUfjoIxg9GrZtgzvugFtugSoh\nHLDydmR6R8w27Nvg68gskDMyd+929i3zBrMffnD2E4mPhzZt1IkpEiSe9DVW+/0C1f7UVPb7haz9\n6QFqf6bg5fuZ/viR9BErb1A6JT1Y+YesE932D1oxLgSrYFNQEynE/vrLOTlg8mQ45xxn/7Pu3YN/\n/jHmEYgAABgFSURBVGbmjsx56+ZxKOVQxo7/tTvQrFqz0HVkpqbCb79lhLLFi2HzZmeXf28oa9NG\n+5aJ+ElNnxbcn5bGvtRUX6Dy/7kv0+39WQSvQ2lplC5WjFP8wtUpmUKW975Y789M173PLVOsWI7T\ngEWRgppIIZOS4jQHjBsHK1bATTc5zQH16wfvM1zvyNyyJXC0bNmyjDMx4+OdS5MmkXUivEgu2fTp\nwX3pQWlfamrAdf/g5Q1a+7IIX8c8Hl9YKpcelMr5hSxvsCqXKXz53++dEiymcBUyCmoihcSaNTBh\nAkyaBI0aOdObPXtCTEz+3zunjkzvAeZx5ePy/0FZOXoUli8PnMLcvz8jkMXHO4vsypcPzeeLBJHH\nL2TtSw9Ne73X/e7L6faBtDRKRkVRLjrauaQHqHJ+ocv/uu9npue6sd5KTp6CmkgEO3bMaQoYPx5+\n/RX69IFbb4Uzzsjf+3o7Mr3BbNnmZTSq3CigI/PUMiGYRrTWSZze0bIffoDff4fGjTNCWZs20KAB\n6A+MuCDZ42FfaqovXO31Xk8PXP73H/czfbrQG7LKe4NWeoAql8vbpxQrRrSbu09LgVJQE4lAv/7q\nbK3xzjvOevj//hcuvzzvo2fejkzvVOYfO/6gRfUWoT8jc88eZ3uMH37IuJQu7YQx76VFC+c+kSBI\n9Xh8oWpPSoovXO3xC105XY55PL6A5f+zfKbb5YoVO+553sClkCUnQ0FNJELs2QPvvecEtG3boF8/\nZ/1ZvXon/16ZOzI37t9IQs0E3xqzkHRkpqQ4Rx/4h7JNm5wg1qZNxmiZNpOVHFhrOeLx+ILVnpSU\njOt+gSvz/d6fh9PSKBcdTYX0AFUhU9DyhqrM93vvK1OsmKYLpUApqImEsbQ0+OYbp2tzxgznvM2b\nboILLsj9OnlvR6Z/MAv5GZnWwrp1gaHs55+hbl1nPZk3lDVpAtEhPptTwtIxj4c9KSnsTg9Qu9OD\nVebr3sC12++6ASdkFS9OhfQAVSH9dnm/2/5BzPtYrLoKJcIoqImEod9/hylT4O23nQGmG2+EG26A\nihVP/FpvR6Z/MCsVXYr2tdvTMa5jaDoy9+6FJUucQOadyoyKCpzCbNkSTgnhhrZS4Gz6/lm708NV\nTj/3ZLovxVpfoKroF7j8r2cVxCpGR1NK3bxShCioiYSJ7dvh/fedgLZli9MY0KePM+iUE/+OzKT1\nSSzcsJDTYk/L2MMs2GdkJic7+378+GNGKNu0ydnhv00bZ8SsTRtnuwyNXEQEay2HPR52p6SwKyWF\nXelhald6sAr4mSmAlTDGF64qpQepin4/K2S6XVFTiCInRUFNxEWHDsHnnztNAQsXQo8ezuhZly7Z\nT21m1ZEZsjMyPR74+++MUPbjj04nw+mnZwSyVq00hRlG0qxlb2oqO72hyy985XQ7yhgqpoct/8BV\nqXhxKvldr+gXyCoUL06MFsaLhJSCmkgBS0lx1p29+y5Mnw5t28J//uN0bZYpc/zzM3dk+p+R2b52\ne9rWahu8jswtW5wpTG8oW7IEypVzQpn30qJF1oVK0KVZy56UFHZmuuxKD2LeMOb/c19qKqd4A1d0\nNJW9YcsvdGV1W9OJIuFJQU2kAHg8kJTkdG1+8omzDVivXnDddcefapRTR2bHuI60PK1lcDoy9+1z\ndvT3D2aHDzsjZN6jl1q1gqpV8/9Z4lvPtSMlhR1+oWtHphC2IznZd31vairl/MJW5Swu3vu9gatC\ndLS2fxApRBTURELEWmf51ocfwgcfOIHs+uudcFanjvc5gWdkhqwj8+hRp+tyyZKMYLZxIzRvnjFS\n1qqVs9eH1g3lirWWfamp7EhJYXt6wPJe3+l323dJTibKGKqkByvfzxIlAm+nh68qCl0igoKaSFB5\nPE44++gj+PhjKFsWrrnGCWhnnOF0ZK7YtsK3vmz++vnBPyMzNdVpG/WGsiVL4M8/nXOlvIGsdWs4\n80ytK/PjHfHanpLC9uRk309v+PJd9wtgpaKiqJIetqoUL86pfte9oauKXxgro+lFETlJCmoi+ZSW\nBosWOVOan3zihLNrr3UC2umNQnxGpscDq1c7YWzpUufnL79ArVoZU5itWjnHF5QqFbxfOkKkejzs\n9Ata2/xC2Lbk5IBAtj0lBYCqJUrw/+3dWWxc133H8e9/9uG+SdxEUbItL0m80JadxA69JKnrBEXS\nBgWaFEWKPvSpRQsULYrmoehb0KcG6EsfkhZIArhpkzbeijpoHFdUJduxrNWSLTmSLFIUJZJauM3M\nvTP39OHOUJcUKZOK5Jmhfh/g4Nxz7rnDQwwk/nDutjmZZHO5rgStSl8lkHUlk2QUvETkFlNQE7kB\nvg+vvx4Gs5/+NLyM62tfgy99dZbLzSvfkTk8OPzrvSPTOTh5MgxklfLOO9DZGYaxnTvD+uGHN/Tz\nyvwgYNL3mfA8zlfKsnYliF0qFmlPJJaEr+5y8OpOJsP+yD6teIlIrVFQE1mjuTl49VV44QV45ZXw\nhoDnvjZF96O7OeFde0fm8NbhG39HpnNw+nR4sf/bb4f1vn3hct3OnVfLww9DV9dN/10/biXnmC6H\nreXl/LL6SqlEVzK5GLS6Uyl6ynU0fHWXV710jZeI1DMFNZHrmJiAF18Mw9nICAw9Ncr2Z3ZR6h9h\n39Quzs6e5fGBxxevL7uhOzIrr1uKBrJ9+yCTCcPYI4+EZefOursDc75UYsLzOFcocK4ctM6tEMYm\nfZ+2RIKecuiqlO5kcjGEVerOZJK4bngQkduEgppIhHOwfz+8/DK89LLj/anj3PvsLjL3jHDa7SJX\nmg9PY24d5qltT63/jkzn4NSp8JTl8lBWCWSV0tt7637RX0Plbsdznsd4JIQtlkjbd47ecsjqLZfu\nct2TStGbTtNTPvWY1MqXyIbhnAMHuPI2LLZxy8ZE9l0zlkj/srFr+gxW+JzVPn+VYxZ/zkcdd71j\n1vGZS/YBrY+2KqjJ7W1uDl57DV56pcQLew/hBnfR/uAIUw0jNGczPLntBu/IrDzVvxLK3nknLI2N\nYRB7+OGaCmXOOWZLJc4WCox7HuPlwDW+vO15JM3oTaXoS6cXA1hvOXhF262JhF4TdAs553AlBwG4\n0tq2F/sCByVwQWR/pG/JuOixwbI6Omal2rH6vsiY5f3X9F1vjFu575r90bks73OR45dvR8esdfyy\neSzfd80817PvVrThxsZE+4jUFRYphKtD0TbG1f8jVhsb2XfN2LV8BlzzOat9/qrHsMbjrnfMej4z\nUj+671EFNbm9OAdHj8JL/1Xg3/f8kiMzIzR/aoT5jj30t/Txhbtu4B2ZxSIcOxYux1UC2YED4YX+\nQ0NXA9nQUFVOX+ZLJc553mIIi4ax6DZAfzls9afTi0Gsb1koa6qRx3o453BFh/MdgRfg/Gu3ne8I\n/KVtV4z0Fd01deAHS/uXl2h/aZXtYjn0LN8XqRf3l66WVftWaOOAGFjcFus1bccM4mCxpf3ElvVF\nxi22o58RHWMrjK38cV0+3q7+PCxyjEU+K/oZsRXmYKscXxkf3W9Lj108xpb9jEiIWDL/lbZX+nmR\nkHDNz79Z+25hG25szJI+IsfJTaFTn3JbuHgRXv7ZLM/v3suesREKvbsobd7HYOO9/OZ9w3zxrnW8\nI3NhIXzfZSWU7d8fJr+BgXCVrFIeegg6Om7p7+Wc42KxyNly4Bor12fLYaxSZkqlq+GrXEe3+8rb\nzR8RwIJiQJALCPIr1PmAoBDZrrTLfa7grrYLkbZX3vbK/Z5b0rfY9pe1vTAwEYdYMoYlDUtZuJ0y\nLGlX+5PL2omrfZZYoT+xdN/iz0jYqoU44Xa8fGw8sj++tE2cpX3xSF9k/JJxK7Vj+oMostEpqMmG\n5Hnw6sgUP/jf3ew6s4vJ7Ai26RiDqYf54t3D/M7Qkzy+dQ13ZE5Ph0HswIGw3r8/vMbs3nvDMDY0\nFJYHHwzvyLyJSs5xwfMYKwewsUgQW9z2PDKxGP2pFAPxFIOlFAPFJH3FJN1+nE1+nA4vTlPBcAsB\npfkSpYUSwXwQ1gsr1LlIOxcsllKuBA5i2RjxbJxYJkYsGwvr5SUdw9K2pF0plral7ZQt3U6Vx6Su\n7rNkpJ2MjEuWV2pERDYoBTXZEEol+O+9o/xwZBe7z4wwnhiB1jEG45/h83c8xTeeGOaJbY+ufkdm\nEITPKDtwICwHD4b1zEy4MjY0FNYPPRQ+0T+V+vXm6xwTnseZXI7xizkmpvJMTuW5NF1g5nKBucs+\n/pUiXbkY3YU4nfk47fkYLTloWDAyC47knCO2EODmAoqzRZzviDfFr5bGcmmKE2uMLbZjDeXthqvb\nsWwsbGdjYV+2vF0ulbYlTSs4IiIfIwU1qUulkuPlvcd5fs8u9pwd4WxiF5aaZ9Ce5Kltw3zzqSf5\n3I5V7sicmwtPXR48eLUcPhyepqyEsUrZto21vPuytFDCn/LxL/p40z4XJ/NMns9zeSrP3LRHftrH\nv+zjLpeIXwnIzDia5yC7AEEKis0xaIkTb02QbkuQbUvS1JYi1Zog0Zog3hwn0RLWlZJoLrebwnYs\nE1OIEhHZYBTUpC4UvBI/3n2Qn/xyN2+c28VEeoR4kGVbbHgxmA3ft+yOzCAIHxp78GD4WqVKOXs2\nXBV78MGr5f77F68nc85RvFLEP+/jXfDwL/h4kx7+pI8/6TN/ocD8hQKFKZ/SVJHYpRIucCy0Glda\n4WKTI9dq0BYn3pEg05GkqTNNa2eazk0Zursy9HY3kO1IkmhNhNdAiYiIrEBBTWrS1JUFfvDaW7x8\neISDF3cz3fAGqXw/dySG+cJdT/KHTw/z6N2ROzKnp8NVsWg5cgTa2+GBB+CBByjd+wBe76fwUn14\nUyW8CY/CuQLehId3zsM7Xy4XfEgbflecXEeMmXaYboWJ5hKjTSUutUFmU4qmzSnaN2fY1JOlvy3D\nQDbL1nSaLek0Wb2KSEREbgIFNam6IHDsPXaGH/3fXnad2suJ/B4WGo7SuHA/9zUM8+w9n+ObzzzB\nPQNdMD8f3mF55MhiGHOHj1CcNwp3fJpC34MUWndQSPdT8NsoTIE37lEYL1CaK5HoTRFsTpDfHGe2\n05jqcJxrc3zYWuSDJp+TzUXS3Sl6WzJszWQYSKcZLNdbMxm2ptO06dlgIiLyMVFQk4/d+PQs/zay\nj1fffZND028xkXgDZ0W6vc8y1PU4X37gs/z+Y5+kY/xDOHKE4OBR8vvGyL93icJ0nHzHfRSatlOw\nHvL5ZgpTcWKZOKktaehLkOtJMLPZuNAFYx0Bp9qKvNficyxTIIgZW8vha7AcvCoBbDCToS+V0rsh\nRUSkZiioyS11eS7Pf+45xKuH3+adiX2cKb1JIXuKpvkHuSv7GM90P8DX2zr4xMQ83v5x8scukxst\nkp/JUkgPkneb8f0s6U5HcluW4p1NzPUnmN4MZzfDqc4S77X7nIiHD2ztTCYXQ9fWTIbByLZWw0RE\npN4oqMlNM3FxjhfeOMQvjh3gnYl9jBb3kW88TmZ2B09e3slzM9t4bKGVnktxvLMl8tMpcoUufGsj\n3bhAottR2J5l5p4WJnc0c6YHTnSGq2Gn/QJXikW2rBLCBtNpBjIZ0loNExGRDURBTdYtCBxvvjfK\nq/sPs/fUIY5dPMBUcJi+eY/P/mqIRybvYMfsZjrnW4jNNZPPtxG3ArGWebx+48odWS7c3cSZHU0c\nH0zxbqvP6ZJHHFYNYYOZDN2pFDGthomIyG1EQU2u69iZSX62/yh7P3iXE6PvETs3zo7pOe6b6GP7\n9HY2z3eTyXXiii1Y8hJeR57Z/hgX7sxw5q5Gjt+Z5eDWBKdTJTYlk4vXhUVD2Nby9WKtNfL+SBER\nkVqhoCZ4fok9Rz9k97vHOX7kJIX3R+kZu8TAZaPvcjftcwOkvW4sSOI1XmF2k8/kYJLR7WmO35Hh\n0LYU53uMLc0NS+6UjF6o359Ok9JpSRERkXVRULtNFEsBb78/xttvnWRs30ni743Rcm6BzitpWvNd\nNOZ6yOTa8dMFZtrzTPTB6LYUH9yR4eSWOK4X2voyDLS2Lj4vbKB8XdhgOk2rLtIXERG56TZcUDOz\n54DvAHHgu865v1+2f0MGtSBwfPDhNO/s/oDzb5yg+P4E8akiDXMZGnMtNC200TjXSjERcKnT51wP\njA6kGN0Sw99UJNFnNA6m6eluZ6Cjg4Hyg1sHMhk2JZO6NkxERKQKNlRQM7M48D7wReAs8EvgG865\nY5ExdRfU/GKJ48fOc2LX+1w4NMr82RncxRLp2RTZXAMNC000zjXSMJ9itqXE+W4Y741zpd0n11Yg\naA9oHEzRdU8z3YOb6O/qoq98OrIrmayblbDXX3+dp59+utrTkBuk769+6burb/r+6ttaglo9XeH9\nGPCBc+40gJn9K/BV4Nj1DqqGy9PzvP/Wrzh9eIzpMxfJTeYoXSoRn4uRzKXI5DJkc1kaFtI0zaQI\nYjDXaRTa+8k3b2KhqUChr0DQViLdF9B7d5LB+7sZGOjlscZGNieTG+7BrfrPpr7p+6tf+u7qm76/\nja+eglo/MBppjwGfvpk/YG4mx+jJC5wfnWJ6/DKzU3MsXMqRm8njzxYJFgLIG7FCjLiXIFlIkSok\nSRdSpPNJMrkEDQtx0gVjpgW85hTxxi7iWQ8/65FvLOBvLmKtedKbAjoG43Tu7OPOT25jU2MjGb1D\nUkRERCLqKait6ZzmP33iRaw80pxhzoiVjFgQIxbYYkkUYyR8I+kbiaKR8iAWQD4DhbSjmHIk0lmy\niRSJVBNe2sdP+vipIsWUT7ExIN+1QKwJkq1JmjY1kNraTsc93XxiaDvNjdm6Oe0oIiIitamerlH7\nDPB3zrnnyu2/AYLoDQVmVh+/jIiIiAhsqJsJEoQ3E3wBGAfeYtnNBCIiIiIbSd2c+nTOFc3sT4FX\nCR/P8T2FNBEREdnI6mZFTUREROR2syGe8WBmz5nZe2Z2wsz+utrzkbUzs382s/Nmdrjac5H1MbMB\nM/uFmb1rZkfM7M+qPSdZOzPLmNmbZnbAzI6a2berPSdZHzOLm9l+M3up2nOR9TGz02Z2qPz9vXXd\nsfW+oraWB+FK7TKzYWAO+L5z7v5qz0fWzsx6gB7n3AEzawL2Ab+tf3v1w8wanHML5WuAdwN/6Zzb\nXe15ydqY2V8AjwDNzrmvVHs+snZmdgp4xDl38aPGboQVtcUH4TrnfKDyIFypA865EeBStech6+ec\nm3DOHShvzxE+fLqvurOS9XDOLZQ3U4TX/n7kHw2pDWa2Bfgy8F1Az4KqT2v63jZCUFvpQbj9VZqL\nyG3JzLYBQ8Cb1Z2JrIeZxczsAHAe+IVz7mi15yRr9g/AXwFBtSciN8QB/2Nmb5vZH19v4EYIavV9\n7lakzpVPe/4Y+PPyyprUCedc4Jx7CNgCPGlmT1d5SrIGZvZbwAXn3H60mlavnnDODQFfAv6kfBnQ\nijZCUDsLDETaA4SraiJyi5lZEvgJ8EPn3E+rPR+5Mc65K8ArwM5qz0XW5HHgK+XrnJ4HPm9m36/y\nnGQdnHPnyvUk8J+El3GtaCMEtbeBHWa2zcxSwO8BL1Z5TiIbnoXvSPsecNQ5951qz0fWx8y6zKyt\nvJ0FfgPYX91ZyVo4577lnBtwzm0Hvg685pz7ZrXnJWtjZg1m1lzebgSeBVZ98kHdBzXnXBGoPAj3\nKPAj3XVWP8zseWAPcLeZjZrZH1V7TrJmTwB/ADxTvsV8v5k9V+1JyZr1Aq+Vr1F7E3jJOffzKs9J\nbowuAaov3cBI5N/ey865n602uO4fzyEiIiKyUdX9ipqIiIjIRqWgJiIiIlKjFNREREREapSCmoiI\niEiNUlATERERqVEKaiIiIiI1SkFNREREpEYpqImIiIjUKAU1ERERkRqloCYiIiJSoxTURERERGpU\notoTEBGpRWb2LeCTwLeBR4EOoNU597dVnZiI3Fa0oiYisoyZPQu8CIwB/wj8CGgBvlTNeYnI7UdB\nTUTkWl3OuSPAE8C/OOcWgO8Bv1vdaYnI7cacc9Weg4hIzTGzRuAicKdzbqza8xGR25NW1EREVjYM\nnFFIE5FqUlATEVnZ54GfV3sSInJ7U1ATEVnZ3cB/VHsSInJ70zVqIiIiIjVKK2oiIiIiNUpBTURE\nRKRGKaiJiIiI1CgFNREREZEapaAmIiIiUqMU1ERERERqlIKaiIiISI1SUBMRERGpUQpqIiIiIjXq\n/wGurOISX0tH4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10565a5c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(r, linear, label='linear')\n",
    "plt.plot(r, huber, label='huber')\n",
    "plt.plot(r, soft_l1, label='soft_l1')\n",
    "plt.plot(r, cauchy, label='cauchy')\n",
    "plt.plot(r, arctan, label='arctan')\n",
    "plt.xlabel(\"$r$\")\n",
    "plt.ylabel(r\"$\\rho(r^2)$\")\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will show how robust loss functions work on a model example. We define the model function as \n",
    "\\begin{equation}\n",
    "f(t; A, \\sigma, \\omega) = A e^{-\\sigma t} \\sin (\\omega t)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which can model a observed displacement of a linear damped oscillator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define data generator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_data(t, A, sigma, omega, noise=0, n_outliers=0, random_state=0):\n",
    "    y = A * np.exp(-sigma * t) * np.sin(omega * t)\n",
    "    rnd = np.random.RandomState(random_state)\n",
    "    error = noise * rnd.randn(t.size)\n",
    "    outliers = rnd.randint(0, t.size, n_outliers)\n",
    "    error[outliers] *= 35\n",
    "    return y + error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define model parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = 2\n",
    "sigma = 0.1\n",
    "omega = 0.1 * 2 * np.pi\n",
    "x_true = np.array([A, sigma, omega])\n",
    "\n",
    "noise = 0.1\n",
    "\n",
    "t_min = 0\n",
    "t_max = 30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data for fitting the parameters will contain 3 outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_train = np.linspace(t_min, t_max, 30)\n",
    "y_train = generate_data(t_train, A, sigma, omega, noise=noise, n_outliers=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the function computing residuals for least-squares minimization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(x, t, y):\n",
    "    return x[0] * np.exp(-x[1] * t) * np.sin(x[2] * t) - y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use all ones as the initial estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x0 = np.ones(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import least_squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run standard least squares:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "res_lsq = least_squares(fun, x0, args=(t_train, y_train))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run robust least squares with `loss='soft_l1'`, set `loss_scale` to 0.1 which means that inlier residuals are approximately lower than 0.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "res_robust = least_squares(fun, x0, loss='soft_l1', loss_scale=0.1, args=(t_train, y_train))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define data to plot full curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_test = np.linspace(t_min, t_max, 300)\n",
    "y_test = generate_data(t_test, A, sigma, omega)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute predictions with found parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_lsq = generate_data(t_test, *res_lsq.x)\n",
    "y_robust = generate_data(t_test, *res_robust.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x106dc2e10>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAGECAYAAACcSOyeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4FVX+x/H33PRAChBIAoQWeijSi5SASBAUUBQsoKx1\nVwVkLT8VhODKuroqAnZUVEQRLLiIAqIGkI4ovUjooUMSSG/z+2MkEgmQcpM7ST6v57kPuXPnnvkm\nBPLJOWfOMUzTRERERERcz+HqAkRERETEomAmIiIiYhMKZiIiIiI2oWAmIiIiYhMKZiIiIiI2oWAm\nIiIiYhO2C2aGYTxlGMY2wzC2GIbxiWEYXq6uSURERKQ02CqYGYZRD7gPaGuaZkvADbjVlTWJiIiI\nlBZ3VxfwF2eBTMDXMIxswBeIc21JIiIiIqXDVj1mpmmeAV4GDgJHgATTNJe6tioRERGR0mGrYGYY\nRjjwCFAPqAlUNgzjDpcWJSIiIlJK7DaU2R5YZZrmaQDDML4EugKzz59gGIY29xQREZEywzRNo6Dn\n2qrHDNgJdDYMw8cwDAPoA2z/60mmaepRio+JEye6vIaK9tDXXF/zivDQ11xf84rwKCxbBTPTNDcB\nHwEbgM1/HH7HdRWJiIiIlB67DWVimuaLwIuurkNERESktNmqx0zsKTIy0tUlVDj6mpc+fc1Ln77m\npU9fc/szijL+6UqGYZhlrWYRERGpmAzDwCzDk/9FREREKiwFMxERERGbsN3kfxERETuxVm8SuTRn\nTrFSMBMREbkCzW2WS3F2cNdQpoiIiIhNKJiJiIiI2ISCmYiIiIhNKJiJiIiI2ISCmYiIiIhNKJiJ\niIgIYG3Z1KtXr0K9JyYmhkmTJunOVSdRMBMRERHAWvqhsMs/KJg5l4KZiIiIAMVbr03BzDkUzERE\nRJxg4cLlREWNJzIymqio8SxcuNzW7c6ZM4emTZvi7e1NixYt+Oqrr/K8np6eztixY2nZsiV+fn6E\nhoYycOBAdu3alXtOdHQ0zz77LAAeHh44HA4cjj+jxcSJE2nbti0BAQFUr16da665hrVr1zql/vJK\nK/+LiIgU08KFyxkzZjGxsZNzj8XGjgNgwIAetmt36dKl3H777dxwww1MmTKFEydO8Mgjj5CZmUnT\npk0BSEtL49y5czz99NPUqlWL+Ph4Xn/9dbp06cKOHTsIDg7mvvvuIy4ujvfee4+VK1fi5uaW5zpx\ncXE88sgj1K1bl+TkZGbNmkWPHj345ZdfaNGiRZHrL9dM0yxTD6tkERGR0lGQnzt9+44zwbzoERU1\nvljXLql2u3btakZEROQ5tmbNGtMwDLNXr175vic7O9tMTk42/fz8zClTpuQenzhxomkYhpmdnX3Z\na2ZlZZmZmZlmkyZNzDFjxhSrfju50vfHH68XOOdoKFNERKSY0tPzH4BKS3PL97gr283OzmbDhg3c\nfPPNeY536tSJevXq5Tk2d+5cOnXqRJUqVXB3d6dy5cokJSWxe/fuAl1r6dKl9OrVi6CgIDw8PPD0\n9GT37t0Ffn9FpGAmIiJSTF5eWfke9/bOtl27p06dIjMzk+Dg4Iteq1GjRu7HCxYs4NZbbyUiIoJP\nP/2UdevWsX79eqpXr05aWtoVr7Nx40b69++Pv78/77//PmvXrmX9+vW0bt26QO+vqDTHrBxbuHA5\n06YtIT3dHS+vLEaP7lusOQkiIpK/0aP7Ehs7Ls9csPDwpxk1qp/t2j3fe3X8+PGLXjt+/Dj169cH\nrJsDGjVqxPvvv5/7emZmJqdPny7Qdb744gs8PT358ssv88w9O3PmDFWqVCly/eWdglk5VVITRkVE\n5GLn/1+dPv0Z0tLc8PbOZtSofsX+/7Yk2nVzc6NDhw7MmzePiRMn5q5btnbtWg4cOJAbzFJSUi6a\nzD9r1ixycnLyHPPy8so9v3LlyrnHU1JS8tyhCfDjjz9y6NAhwsPDi1x/eadgVk5Nm7YkTygDiI2d\nzPTpzyiYiYiUgAEDepTI/68l0e6kSZPo27cvgwcP5v777+fkyZNER0cTEhKSux7Zddddx9dff80/\n//lPBgwYwIYNG3jttdcIDAzMs2ZZREQEAC+//DL9+vXDzc2N9u3bc9111zF16lRGjhzJyJEj2b17\nN8899xy1atXSmmeXoTlm5VRJTUQVEZGy75prrmH27Nns2rWLIUOG8PLLLzN16lSaNGmS24N23333\nMW7cOD777DMGDhzIokWLWLBgAQEBAXl2B7j++ut58MEHeeONN+jatSudOnUCoG/fvkybNo2VK1dy\nww038MEHHzBr1iwaNmxY6N0FKhKjrKVWwzDMslazK0RFjWfJkufyOf4Mixb9ywUViYiUTYZhqIdH\nLulK3x9/vF7gJKoes3Jq9Oi+hIePy3PMmjB6rYsqEhERkStRj1k5tnDhcqZP//6CCaPXan6ZiEgh\nqcdMLsfZPWYKZiIiIpehYCaXo6FMERERkXJKwUxERETEJhTMRERERGxCwUxERETEJhTMRERERGxC\nwUxERETEJhTMRERERGxCwUxERKSCmT9/PlOmTHF1GZIPBTMREZEKZv78+bzyyiuuLkPyoWAmIiIi\n+UpPT3d1CRWOgpmIiEgFMnLkSD766CPi4uJwOBw4HA7q16/PsmXLcDgcfPXVV9x3331Ur16d0NDQ\n3PfUr1//orYiIyPp1atXnmMnT57k73//O7Vr18bb25tmzZoxY8aMUvncygN3VxcgIiIipWfChAmc\nOnWK9evXs2DBAgC8vLxISEgAYNSoUfTv35/Zs2eTlpaW+z7DuHi7R8Mw8hw/e/Ys3bp1Iz09nUmT\nJlG/fn0WLVrEP/7xD9LT03n44YdL+LMr+xTMREREKpAGDRoQFBSEp6cnHTt2zD0eExMDQKdOnXjn\nnXcuel9+G3WbppknmE2dOpWDBw+ydetWwsPDAejduzcJCQlMmjSJBx98EIdDg3WXo2AmIiLiRMak\ni3uWnMWceHE4crYbb7yxyO9dtGgRnTt3pl69emRlZeUe79u3L++++y7bt2+nRYsWziiz3FIwExER\ncaLSCE8l6fy8sqI4ceIEsbGxeHh4XPSaYRicPn26OKVVCApmIiIikiu/uWTe3t5kZGRcdPz06dNU\nr14993lQUBAhISFMnTo137YbN27svELLKQUzERGRCsbLy4vU1NQCn1+3bl2OHz/OqVOnCAoKAiA2\nNpZdu3blCWb9+vVj+vTphIWF5TkuBacZeCIiIhVMREQEZ86c4a233mL9+vVs2bLlsucPHToUwzAY\nPnw4ixcvZvbs2QwePJjq1avnuSlg7Nix1KhRg+7du/P222/z008/8c033/DSSy8xePDgkv60ygX1\nmImIiFQw9957L2vWrOHpp58mISGBevXqMXPmzHyHMQHCw8P5/PPPGT9+PDfeeCNNmjRhypQpTJ48\nOc97/P39WbVqFc8++ywvvPACcXFxBAYG0rRpU4YMGVJan16ZZuR3+6udGYZhlrWaRUSk7DIMI9+l\nIkTgyt8ff7xe4Ft1NZQpIiIiYhMKZiIiIiI2oWAmIiIiYhMKZiIiIiI2oWAmIiIiYhMKZiIiIiI2\noWAmIiIiYhMKZiIiIiI2oWAmIiIiYhMKZiIiIiI2oWAmIiJSgURHR+Nw6Me/Xdnub8YwjEDDMD43\nDGOHYRjbDcPo7OqaREREypNLbVYurufu6gLyMRX41jTNmw3DcAcqubogERGR8kSbstuXrXrMDMMI\nALqbpvk+gGmaWaZpJrq4LBERkXJr6tSpNGvWDF9fX6pWrUqHDh2YP39+7uvZ2dmMHz+e0NBQKlWq\nRK9evdi2bRsOh4NJkya5sPLyyW49ZvWBk4ZhzARaA78AY0zTTHFtWSIiIuXP7Nmzeeyxx5g4cSLd\nu3cnNTWVTZs2ER8fn3tOdHQ0zz//PI8++ih9+/Zl/fr1DBw4ENCQaEmwWzBzB9oCD5umud4wjFeB\nJ4EJri1LRESk/Fm9ejWtWrVi/Pjxucf69euX+3F8fDxTpkzhgQce4MUXXwSgT58+uLm58eSTT5Z6\nvRWB3YLZYeCwaZrr/3j+OVYwyyM6Ojr348jISCIjI0ujNhERkSsryV4kJ88N69ixI2+++SajR49m\n4MCBdO3aFV9f39zXt2zZQkpKCkOHDs3zvltvvVXB7BJiYmKIiYkp8vttFcxM0zxmGMYhwzAam6a5\nG+gDbPvreRcGMxEREVspQxPr77zzTtLS0njvvfd444038PDwoH///rzyyivUrVuXo0ePAhAcHJzn\nfTVq1HBFuWXCXzuMCjsPz1aT//8wCphtGMYmoBXwbxfXIyIiUm7df//9rF27ltOnT/Phhx+ybt06\nhg0bBkBoaCgAx48fz/Oevz4X57FdMDNNc5Npmh1M02xtmuZNuitTRESk5AUEBDB06FBuueUWtm7d\nCkCrVq2oVKkSn332WZ5z58yZ44oSKwRbDWWKiIhI6bn//vvx9/enc+fO1KhRg927d/Pxxx8TFRUF\nQGBgIGPHjmXy5Mn4+flx7bXXsn79et5//30XV15+KZiJiIhUIIZh5C5z0a1bN2bOnMmsWbNITEyk\nZs2ajBgxIs+8qOjoaEzT5N133+W1116jc+fOLFiwgIiICFd9CuWaUdZW/zUMwyxrNYuISNllGIZW\nys+Hw+EgOjqaCRMq9opWV/r++OP1At+qa7s5ZiIiIiIVlYKZiIiIiE1ojpmIiIgUWk5OjqtLKJfU\nYyYiIiJiEwpmIiIiIjahYCYiIiJiEwpmIiIiIjahYCYiIiJiEwpmIiIiIjahYCYiIiJiEwpmIiIi\nUigffPABDoeDvXv3luo1Z86cWaBzY2JicDgcLF++vISrcj4FMxEREbG9Dz74gPfff9/VZZQ4BTMR\nEREhPT3d1SUICmYiIiIVTnR0NA6Hg23bthEVFYWfnx/Dhg0D4OjRo9x5551Ur14db29vWrduzezZ\ns/NtJy4ujsGDB+Pn50dQUBAPP/wwaWlpua9fakjx/FDowYMHc4998skntGnTBj8/PwICAmjVqhXv\nvPMOAJGRkSxfvpyVK1ficDhwOBz07t27UJ/z4sWL6dq1K4GBgfj5+dG0aVP+9a9/5Tlnzpw5NG3a\nFG9vb1q0aMFXX31FZGQkvXr1KtS1ikN7ZYqIiFRQgwYN4t577+Wpp57C4XCQnJxMz549SUxM5Pnn\nnycsLIxZs2YxYsQIUlJSuO+++/K8f/jw4QwbNoyHH36YtWvX8uyzz5KcnFzguWDn/fzzz4wYMYIx\nY8bw8ssvk5OTw44dO0hMTATgzTffZPjw4eTk5PD2228D4O/vX+D29+7dy8CBAxk6dCjR0dF4enqy\ne/du9u3bl3vO0qVLuf3227nhhhuYMmUKJ06c4JFHHiEzM5OmTZsW6vMpDgUzERGRCmrMmDGMGjUq\n9/lrr73Gnj17iImJoUePHgBERUVx/Phxxo8fz7333othGLnnDxgwgBdffBGAPn36YBgGEyZM4Omn\nn6ZRo0YFrmPNmjUEBgbyyiuv5B7r06dP7sfNmjXDz8+PnJwcOnbsWOjPc+PGjWRmZvLmm29SuXJl\nwOqFu9DEiRNp3rw5X3/9de6xpk2b0qVLFwUzERGRssqIiSmxts2/hIniuvHGG/M8X758ObVr184N\nZefdcccd/PDDD2zfvp2IiIjc40OHDs1z3rBhwxg/fjzr168vVDDr2LEj8fHxjBgxgmHDhtGtWzcC\nAwOL8Bnlr02bNnh4eDBs2DDuvvtuunfvTo0aNXJfz87OZsOGDTz11FN53tepUyfq1avntDoKQsFM\nRETEiZwdnkpSaGhonudnzpy56BhASEhI7usXCg4Ozvd5XFxcoero0aMH8+bNY/r06dx0000A9OzZ\nk1deeYWWLVsWqq38hIeHs3jxYl544QVGjBhBeno6HTt25IUXXqBHjx6cOnWKzMzMiz4fIE+AKw2a\n/C8iIlJBXTgsCVC1alWOHj160XnHjh3LfT2/4+cdP34cgFq1agHg7e0NQEZGRp7zTp8+fdE1hgwZ\nQkxMDAkJCXz11VccPXqUfv36FebTuazIyEi+++47EhMTWbp0Ke7u7gwYMIAzZ84QFBSEh4dHbv35\nfU6lRcFMREREACu8HD58mFWrVuU5/sknnxAcHEzz5s3zHJ87d26e53PmzMHhcNCpUycA6tatC8CW\nLVvynLdw4cKLQuF5vr6+DBgwgPvvv5+jR4/mhjgvLy9SUlKK/sn9wcPDg169evH444+TnJzMvn37\ncHNzo0OHDsybNw/TNHPPXbt2LQcOHCj2NQtDQ5kiIiICwMiRI5k6dSo33XQTkydPplatWsyePZul\nS5fyzjvvXBSmvvvuO5544gmuvfZa1q1bx7PPPstdd91FeHg4YA2V9uzZk+eff56goCCqV6/Oxx9/\nzL59+/IEoAkTJnDixAl69epFaGgohw8fZtq0abRp04Zq1aoBEBERwRtvvMHcuXNp0KAB/v7+NG7c\nuECf11tvvcWKFSvo378/tWvX5tSpUzz//PPUqlWLFi1aADBp0iT69u3L4MGDuf/++zl58iTR0dG5\nw7ilxjTNMvWwShYRESkd5fHnTnR0tOlwOMzs7OyLXjt69Kg5YsQIMygoyPTy8jJbt25tzp49O885\nM2fONB0Oh7lixQpz0KBBZuXKlc1q1aqZDz/8sJmWlpbn3MOHD5s33HCDGRgYaIaEhJjjxo0z3333\nXdPhcJgHDhwwTdM0Fy5caEZFRZmhoaGml5eXGRYWZt57773m0aNHc9s5duyY2b9/f9PPz880DMPs\n1avXJT+/n376yXQ4HOayZctM0zTN1atXm4MGDTLDwsJMLy8vMzQ01Bw6dKi5e/fuPO/79NNPzSZN\nmpheXl5mixYtzPnz55uRkZGXvdaVvj/+eL3AOccwL0isZYFhGGZZq1lERMouwzDQz52KKzIyEofD\nwY8//pjv61f6/vjj9fzHbfOhOWYiIiIil1GawVzBTEREROQSDMO45I0KJXK9stY9q6FMEREpTRrK\nlMvRUKaIiIhIOaVgJiIiImITCmYiIiIiNqFgJiIiImITCmYiIiIiNqEtmURERK6gNJdLkIpNwUxE\nROQytFSGlCYNZYqIiIjYhIKZiIiIiE0omImIiIjYhIKZiIiIiE0omImIiIjYhIKZiIiIiE0omImI\niIjYhIKZiIiIiE0omImIiIjYhIKZiIiIiE0omImIiIjYhIKZiIiIiE0omImIiIjYhIKZiIiIiE0o\nmImIiIjYhIKZiIiIiE0omImIiIjYhIKZiIiIiE0omImIiIjYhIKZiIiIiE0omImIiIjYhIKZiIiI\niE3YMpgZhuFmGMavhmEscHUtIiIiIqXFlsEMGANsB0xXFyIiIiJSWmwXzAzDqA30B94FDBeXIyIi\nIlJqbBfMgCnA40COqwsRERERKU22CmaGYVwPnDBN81fUWyYiIiIVjLurC/iLrsBAwzD6A96Av2EY\nH5mmeeeFJ0VHR+d+HBkZSWRkZGnWKCIiIpKvmJgYYmJiivx+wzTtOb/eMIyewGOmad7wl+OmXWsW\nERERuZBhGJimWeBRQFsNZeZDCUxEREQqDNv2mF2KesxERESkrChvPWYiIiIiFYaCmYiIiIhNKJiJ\niIiI2ISCmYiIiIhNKJiJiIiI2ISCmYiIiIhNKJiJiIiI2ISCmYiIiIhNKJiJiIiI2ISCmYiIiIhN\nKJiJiIiI2ISCmYiIiIhNKJiJiIiI2ISCmYiIiIhNKJiJiIiI2ISCmYiIiIhNKJiJiIiI2ISCmYiI\niIhNKJiJiIiI2ISCmYiIiIhNKJiJiIiI2ISCmYiIiIhNKJiJiIiI2ISCmYiIiIhNKJiJiIiI2ISC\nmYiIiIhNKJiJiIiI2IR7QU4yDGMRsBf4CYgxTfNkiVYlIiIiUgEZpmle+STD6AkMAnoArYHdwI9Y\nQe1b0zTTSrLIv9RiFqRmEREREVczDAPTNI0Cn1/YkGMYRgDQHbgVK6ylAfebpvlVoRoqIgUzERER\nKStKPJj95WKjgdXAq8AE0zR/KHJjBb+mgpmIiIiUCYUNZgWa/G8YRrRhGL/+8WeDC18zTXM90BPo\nV7hSRURERORCBb0r0w14AqgPbDEM43fDMH4BuvzxekNgXwnUJyIiIlJhFOiuTOAYgGmadxmG8TBw\nNeADfGsYRiCwBXirZEoUERERqRgKPMfMMIxuf5y/Ip/XmgBHTNM85+T68qtDc8xERESkTCjVyf+u\noGAmIiIiZUWJTP4XERERkZJX0DlmUlAnTsCGDbBtG5w7B6YJtWtDkybQuTN4e7u6QhEREbEp9Zg5\nQ3Y2zJ4NvXpB48bw6qtw5Ai4u1uPDRvgqaegRg0YPBgWLYKcHFdXLSIiIjajOWbF9e23MHYsVK8O\njz4K11136V6x+Hj48kt47TXIzISXXoJ+Wv5NRESkvNLk/9Jy7hyMGQMxMfDmm9C3LxgF/LqbJixY\nAI89Bk2bwjvvQEhIiZYrIiIipU+T/0vDoUPQrZsVsDZvhqiogocysM4dOBC2boVWreCqq+Cbb0qu\nXhERESkT1GNWWNu3W71jY8fCP/9ZuEB2KStXwrBh8OCD1lw0Z7QpIiIiLqehzJK0ezf07g3/+Q8M\nH+7cto8cgRtvhObNYcYM66YBERERKdMUzErKkSPQpQtMnAh3310y10hOhiFDwMcH5swBL6+SuY6I\niIiUCgWzkpCaCj17wqBBMG5cyV4rIwNuuw2ysuDzz8HDo2SvJyIiIiVGwczZTBNGjLDWKvvkk9KZ\n/5WRATffbC278emn4OZW8tcUERERp9Ndmc724Yfw22/w3nulNynf0xPmzoXTp60lOcpYeBYREZGi\nUY/Z5fz+O3TtCj/+CC1bls41L5SYCN27WzcaPPFE6V9fREREiqWwPWa69e9SsrPhzjvhmWdcE8oA\nAgKsnQU6d4YWLaB/f9fUISIiIqVCQ5mX8tZb1pIVDz/s2jpq14bPPoORIyE21rW1iIiISInSUGZ+\nDh+GNm1g+XJo1qxkr1VQ06db89xWrQJfX1dXIyIiIgWguzKdYdgwaw/LSZNK9jqFcf7uUIfDuiFB\nuwOIiIjYnoJZcf38M9x+O+zcab+eqZQUa5HbBx6wtm8SERERW1MwK46cHOjUyVqiwtlbLjnLnj1W\nOFu2zNq+SURERGxL65gVx6efWkOEt9/u6kourWFD+Pe/reCYkeHqakRERMSJ1GN2XmamNdF/xgzo\n1cv57TuTacLAgdCqFUye7OpqRERE5BLUY1ZUs2ZBnTr2D2Vg9eq9+y68/z6sXOnqakRERMRJbNVj\nZhhGGPARUAMwgXdM05z2l3Oc32OWkQFNmsDHH8PVVzu37ZL09dcwdixs2gR+fq6uRkRERP6iTE/+\nNwwjBAgxTfM3wzAqA78Ag03T3HHBOc4PZu++C/PmweLFzm23NNxzj7XZ+euvu7oSERER+YsyHcz+\nyjCM+cB00zR/uOCYc4NZdrZ1d+Pbb0NkpPPaLS3x8RARYQXLstTbJyIiUgGUmzlmhmHUA9oAa0v0\nQl9/DVWqQM+eJXqZElOlCkydCvfdB+nprq5GREREisGWweyPYczPgTGmaSaV2IVME154Af7v/8r2\nSvo332wto/Gf/7i6EhERESkGd1cX8FeGYXgAXwAfm6Y5P79zoqOjcz+OjIwksqhDkD//DAkJMGhQ\n0d5vF4ZhzTFr0waGDrXP/p4iIiIVTExMDDExMUV+v63mmBmGYQAfAqdN0xx7iXOcN8fslluseWUP\nPeSc9lzt9detRXKXL7f21BQRERGXKtOT/w3D6AYsBzZjLZcB8JRpmosuOMc5wezgQbjqKjhwoPws\nNZGTA127wv33w913u7oaERGRCq9MB7OCcFowe+opSE2FV18tflt28ssvMGAA7Nhh3RggIiIiLqNg\nVhBpaRAWBqtWQaNGzinMTv7xD3Bzg9dec3UlIiIiFVq5WS6jRH3xhTVRvjyGMrD2z5w3D377zdWV\niIiISCFUzGA2Y4a17ld5VbUq/Otf1k0NOTmurkZEREQKqOIFs127rPlXZX2JjCu55x5rD9BZs1xd\niYiIiBRQxZtj9vjj1rpfL77olHoWLlzOtGlLSE93x8sri9Gj+zJgQA+ntF1s69ZZAXTXLvD3d3U1\nIiIiFY4m/19OZibUrg0rVkDjxsWuZeHC5YwZs5jYuEnQ/CzUTSGw0ff07FOfJg3CqOXlRWMfHzr7\n+xPo4VHs6xXJXXdBzZrw/POuub6IiEgFpmB2OQsWWNsWrVxZ7DpSs7Pp8H9vsS30aiuU/e4HeyvB\nKS+a1F/MXfdGcTg9nZ0pKaw7d44IX1/uCA7mtho1CPL0LPb1CywuDlq1spbRqFev9K4rIiIiCmaX\ndfPNcO218MADRb5+ek4Ob8TF8d9Dh0jfepQzs7vDmmqQ4ZZ7Ts+e0cTEROd5z/KEBD48doxvz5xh\nRHAwT9SpQy0vr4vaL5Gh0UmTYPt2+Oyz4rUjIiIihVLYYGa7vTJLzJkz8P331h2ZRbQ8IYEHdu8m\n3Nubb1u25P8mfcuS5UMuOs/bOzvPcy+Hg2urVuXaqlU5lp7Ofw8dovX69TwaFsY/w8Lw+mP7pNyh\n0djJue+NjR0HULxw9thj0LSp1VN49dVFb0dERERKVMXpMXvrLfjxR5g7t9BvzcrJYeL+/Xxw7Biv\nNWrE4KAgDMO4IEg9B0G7IGgnQfXfZOCQhrSMaERwpWDqV6lP6+DW+Hj45Glzb2oqj+zZQ2xqKrOb\nNeMqPz+iosazZMlzF10/KuoZFi36V+E/5wvNmmUtOLt6tfbRFBERKSXqMbuUWbOsbZgKKT4zkyHb\ntuFuGPzavj01/pgflpWThdnoLLXHrOPg8Uo4sjzxT61BuybNCKjtxd74vaw+vJo9Z/aw69QuWoe0\nZkizIQyLGEYt/1o08PHh6xYtmH38ONdu3sy4OnVIS8//ryMtzS3f44Vyxx0wbZq1yfkddxS/PRER\nEXG6ihHMDh2CnTuhb99Cve1gWhr9Nm8mqmpVXgoPx80wyMjOYMYvM3hp9UuEVg5leOfhzG7yAbX8\na12yneSLDR5RAAAgAElEQVSMZFYcXMG8bfNotaIVfRr04YmuT9CuZjuGh4RwdUAAN23bxvGbm8Dq\n7Dzz1eDiodEicThgyhS4/Xa48Ubw9S1+myIiIuJUFWMoc8oU2LIF3n+/wG/Zl5pK5G+/Mbp2bR4N\nCwNg0Z5FjFk0hvAq4TzT4xm6hHUpXB3A2fSzvLfxPV5a/RLd63Tn+Wuep36V+qRkZ3PdTytYuz+B\n9EcHwFlreY3w8KeZOrWf89ZGGzoUWraEZ55xTnsiIiJySborMz9du8KECdCvX4FOP5CWRuRvv/FY\nWBgP1apFUkYSYxeN5Yd9P/Ba/9fo36h/ESrPKzkjmVdWv8LUtVN5qttTjOk8BofhxtDvl7HoXDIt\nP9xMQEYKo0Zd69wFa/ftg/btraBas6bz2hUREZGLKJj91cGD1oblx45BARZ5PZ2ZydUbN/JAzZqM\nDQtj16ldDJwzkK5hXZnWbxp+Xn7FqP5ie87s4b4F95GZncknQz6hTkAdXjp4kLePHmXZVVdRM58l\nNYrtySfh+HGYOdP5bYuIiEiuwgaz8n973uefw+DBBQpladnZDN66lRuCghgbFsbSvUvp8UEPnuj6\nBDMHzXR6KANoWLUhP9z5Azc0voEOMzrw3e/f8VidOtwTEkLv337jWHq606/JU0/Bd9/Bpk3Ob1tE\nRESKrPz3mHXpAtHREBV12dNM02Tkzp2k5uQwp3lzvtzxBQ99+xDzbplHj7qls/flyoMruWXeLTza\n5VH+2eWfTD5wgE9PnOCnq67KvRvUaV5/HebPhyVLrL1DRURExOk0lHmhgwehbVs4evSKPWavx8Xx\n9pEjrG7blq+2fcoT3z/Bd3d8R+uQ1k6ouuAOJh5k4KcD6VG3B6/2e5VJ+w8w/9QplrdpQ4C7E2+i\nzcyEFi1g6tQCz70TERGRwtFQ5oUKOIy5JjGRZ/fv56sWLVi8+2se//5xfrjzh1IPZQB1AuoQMzKG\njUc3cs//7mF8ndp0Dwhg6LZtZObkOO9CHh7w4ovWrgBZWc5rV0RERIqsfAezuXOt5SEu42xWFnfs\n2MFbjRuz/+gq/v7N31l4+0KaVW9WSkVeLNA7kMXDFxN3No7bvriNF+rXwc0wGPX77zi1h3PgQKhW\nDT74wHltioiISJGV36HMAwesZSGOHLlsj9mIHTvwdTgYW82kx8wezL1lLpH1Ip1XcDGkZ6Vz2xe3\nkZqVykdDPueaLdu5KyQkd101p1i/HgYNgt27oXJl57UrIiIiGsrM9cUXVxzGnHfiBOvPnmVczapc\n/8n1vNDnBduEMgAvdy/m3jKXSh6VeHDBSL6OaM4rhw4x/+RJ512kQwfo1Qteesl5bYqIiEiRlN8e\ns8hIePxxGDDgopcWLlzOSzN+YtU9nWn91W9ktfyCayJ68t++/3V+wU6QnpVO/0/607BKQ+7p8SID\ntm7l5zZtaOKsbZUOHLBuktCisyIiIk6luzIB4uOhbl1rEVUfnzwvLVy4nDFjFhN7++1w0gt2fox3\n05l8NmAmA6/vXYKVF8+59HP0/qg3UeFR1Gn6INMOH2Ztu3ZUcnPCBudgLTp76hS8+65z2hMWLlzO\ntGlLSE93x8sri9Gj+zp3FwcREbG9wgaz8rmJ+aJFVo/ZX0IZwLRpS4it+ig02Q2fHYPrZ5D2zi+8\nse91WwczPy8/vr39W7rP7M7ffarTLrAvD+zaxaxmzTCcsQ7ZU09B48aweTO0alX89iq43F8AYifn\nHouNHQegcCYiIpdUPueYLVgA11+f70upWR7w8O8wIxiuvwvmfwBJIaSlOannqQRVr1SdxcMX8+Kq\nFxjk2Mfm5GTeOnLEOY0HBFgbmz/+uHPaq+CmTVuSJ5QBxMZOZvr0711UkRTHwoXLiYoaT2RkNFFR\n41m4cLmrSxKRcqr89ZhlZVk9Zv/Nf77YsatrQpwPBD8GW4fB3msB8PbOLs0qi6xuYF0+H/o5g+cM\n5sPblnLn/v209/Ojg79/8Rt/4AGYNg0WL77iTglyeenp+f/TKgu/AEhe6v0UkdJU/nrMVq6E+vWh\nVq2LXopLT+d4r0ZU++UVCNoFP/wbgPDwpxk16trSrrTIuoZ15aW+L/HwlzfyQt0Qbtu+nXPOWCT2\nwkVns8tGULUrL6/8/z7Kyi8A8if1fopIaSp/wWzBArjhhnxfeiI2lntq1SCn7QI6He9Az27PExX1\nDFOn9itzv/ne2fpObm52Mx8tvYfuAf6M3rPHOQ0PGgRVqmjR2WIaPbov4eHj8hwra78AiEW9nyJS\nmsrfUOY338Ds2RcdXp6QwIrERHqemMnwNncw7elpLijOuf59zb8Z/NlgvPbP4MeqtzL3xAmG1qhR\nvEYNw1rT7MYbYdgwLTpbROeD/vTpz5CW5oa3dzajRpW9XwBEvZ8iUrrK13IZv/8OPXvC4cPg+LMz\nMMc0af/LL/T3OsusH/7G1n9sxc/Lr5QqLlkJaQm0e6cdd3d/gakpoaxv14663t7Fb/j226FJE5g4\nsfhtiZRh+c0xCw9/ukz2tItI6avY65hNmQI7dsA77+Q5/NmJE7x48AAJq25jWr+pDGh88aKzZdmv\nR3+l78d9GTHgOzaku/HTVVfhVtwlNPbvh3btYOtWCA11Sp0iZdXChcuZPv37C3o/r1UoE5ECqdjB\nrHdveOQRa3PuP2Tl5NB8/XraJC3HkbCRT4d8WkqVlq73Nr7HS6tfIajLhwwMqsHjdeoUv9EnnrAW\n650xo/htiYiIVEAVN5glJECdOnDsGFywVdG7R47wXtx+dv80kG0PbiWkckgpVlu67v76bk7leLCq\n+giWt2lD80qVitdgQoI1nLl0KbRs6ZwiRUREKpCKu4n5Dz/A1VfnCWVp2dk8e+AAnodm8UTXx8t1\nKAN4vf/rHDi+hn7uRxi5cydZOTnFazAwEMaPt3rOREREpMSVn2C2ZAn07Zvn0JtHjhDmSOPQocWM\n6TzGRYWVHh8PH+YMmcOiFQ/hnpPOi4cOFb/RBx6A2Fjr6ysiIiIlqnwEM9O8KJidy8riPwcPcmr7\ni7zQ5wW83Z1wp2IZ0Kx6M/7dezKJm8cx5dAhtiQlFa9BT0944QUtOisiIlIKykcwi42FjAxo3jz3\n0LS4OOqbZ6hBEjc3v9mFxZW++9reR1O/qrRJXc9dO3eSWdwhzcGDrb00P/zQOQWKiIhIvsrH5P83\n34S1a3NXq0/Ozqbe6tWw6RG+u/Et2tdsX/qFutiZ1DO0fusqQjrP5PrQcCbWq1e8BteuhZtugt27\nobg3FUihLVy4nGnTlpCe7o6XVxajR/e1zXINdq5NRMTVCjv5v3ys/L9kCdz8Z6/YjCNHCMo8Qsda\nERUylAFU9anK7Js+5uav/870tm9zQ7VqtPUrxqK6nTpBjx7w8sswYYLzCpUrcvYm2s4MUtrgW0TE\nucp+j1lWFlSvDjt3QnAw6Tk51F+9iuRfx7J5+OfUDazrumJtYMJPE/jyTBKOOrexoV07PB3FGL3e\ntw/at9eis6UsKmo8S5Y8l8/xZ1i06F+Faiv/VezHMXVqVJGClDNrExEpjyrechnr1kHduhAcDMCs\nY8fwSj/CiAadKnwoA5jQcwL+iWvJSonj+YMHi9dY/fpw993apqmUOXMT7WnTluQJZQCxsZOZPv17\nl9cmIiLlYSjz++9z78bMyslh8oF9nN41naeHf+biwuzB3eHOJzfNpu0HUbzqeIcbg4JoVZyNyZ9+\n2lp0dvRoaNHCeYXKJTlzE+3CBKnMnBz2pKayJzWV2NRUYtPSOJ2ZydmsLM5mZZKcnUXsPS2h/2+Q\n5g7J7nDCC455k1zHn6Pp6YR6eRW6RhGRiqzsB7MlSyA6GoDPT54kPeUo94V3oqZfTdfWZSP1Ausx\nrfczPLFlFiO9fVnXrh3uRR3SrFIFxo2zFp399lvnFir5Gj26L7Gx4y7aRHvUqH6FbutyIS81O5vl\niYn8nJjIj2dOsjEpGV8zDZ/MM5iph0k5u4fUlKOkZyTgRTZehklaVg54GeBdCQKrQ1hD3H0qsy+k\nDo1Wr8DHzYOO/gF0DAgkMjCQLv7+xRtOFxEp58r2HLPERKhdG06cwPT2pvnaVRzeFM3e4Z9QvVJ1\n1xZqM6ZpctuXt7MqcDAPNuzMk3WLMcybkQEREfDGG3Dttc4rUi7JWZtoXzTHzDuL4IGvU/+eCDZ7\nOaiUcYLMMxvIjP+VTn6V6RjcnPCq4TSo0oD6gfUJ8g3C18MXN8efPWzffLOMaa8vIik7kxz/M/Qa\nVI+g8EpsO7mdNaf2EZvpRkBQR4wq7Uhyr0L3gACuD6rBTdWrU1M9aiJSzlWsvTLnz7fCwZIlfHv6\nNCN+W8bfjY1M7n3xZGSBhLQEIt7rTVLLl1nTriPNirPsxZdfwqRJsHEjuGk+UVmycOFyJn+2mr1X\nhXIyoga+6QfIPPk9kX7u3BTem251utE0qCkOwzk9W+lZ6Ww+vpmY/TF8e2Alq5PS8Qu9lmS/VrSo\nVIk7Q8O4uXp1Qry8tPSGiJQ7FSuYPfSQNSH9sce4esNqNm96iYO3v0sVnyquLdLGlu1fxsDl79Go\n2UOsbd8RN6PA3yt5mSZ07w733gsjRzq1RikZOabJ5ydP8lzsdmJTk8mO+5oo30xGNh9I3/C+VPIs\nnfXp0rLSWHFgBXN3fMW8o/twBPchzb81zTNyiHtzG8f+9wjkWN+XxbljVETEDipWMGvcGObNY3N4\nOF3WLWessY7nekW7tL6y4P+WPsl72RE83TySf4aFFb2htWthyBDYtUuLztpYjmky78RxHtu9ldNJ\nh6lychGPNOzI3666kyDfIJfWlp2TzYqDK5i5+TM+jjtOTrWh4F4DvqwPX9eGJA8tvSEiZVrFCWaH\nDkG7dnDsGMO2/sr/Nr3F4SH/oZpvNVeXaHsZ2Rm0+fB6DjR4gt86dqWhr2/RG7vjDmu5kn//23kF\nitP87+RxHtqxiZNJR6ifGMNL7YbSv9F1GEXtKS1B3Xs/zc9nmkCfhdC6BVTrCt8H03HfQtbOf9LV\n5YmIFEnFWcfsp58gMpLjWVl8ffo0fwsOUigrIE83Tz6/YSrmgY+5betv5BQnnP/3v/DOO9ZWTWIb\nv6ek0HH1Um7eGEPA0S9Y2LQO22+dyYDG/W0ZygB8PRyw6S54eS48MRA++RJazmDDP1py3cbVrD97\n1tUlioiUuLLbY/a3v0HHjjx+TXem/foJ+/uPItRPq9EXxmvr3uDJ074836I3o8LqFL2hl1+21pP7\n7juw6Q/9iiIlO5vR29fz0YkzBJ5cxIy2/Rho0x6yv8pvV4K6EaPpMDqFZW5ZpNYeTGvTwSspWXQ8\nfty6Kzst7c+HaYKHR96HpycEBkK1ankf1atbr4kUk25YkSupGEOZOTlQrx5pixYRdOwg16f8wJwB\nL7q6tDLHNE16zbuLddWHs71zd+r5+BStocxMaN3aGs4cPNi5RUqBLTl1lGGbfyE1fhPP1AriiQ53\n4+Hm4eqyCi4tjRVvz2LNe99QK+E0DdJOEOGThd+Zk5gOB8dCApnSvT0fDriNFufO8vyho3TMzARv\nb/Dysn4pyMzM+0hPh4QEOH36z8eZM3DqFAQFWcPwFz4aN4ZmzaBmTf2SUY45K0w5e4szKZ8qRjCL\njYWrr+a1NSv4568L2NVzMPWr1Hd1aWXS8aTjNJw/gaYNb2Ndp55F71n58Ue45x7Ytg2KM2dNCi0l\nO5vbN/7IN/FJ9Mj8lXm9xpSNYf2DB2HZMlizxrqRZPt2aNjQWiMvIgKaN7ee16lj9XoBmdmZfLj5\nE57e/jOJwdfTprIfrzdvTzt//8JdOzsbjhyBAwf+fOzfbw3J79hh9cA1bWqFtIgIaNPGelQrA19X\nuSxnhim77xWr3jx7KGwwK5sr///0E2avXjy3bxc93I4rlBVDcOVgPu5wA0Nj9/HaoTqMqtOgaA31\n7g0dO8J//gPPPuvcIuWSYk7FMXDTBsyz2/ksoj1DGtv4a5+UZA15f/89LF0K8fEQGQldusDtt0Pb\ntnCFXlsPNw/ubXMXf2s9nI+3zuXJbYvpkpREN38/ZrTsTPgF77/sDyU3NwgLsx7dul18odOnYedO\nK6Rt3QrffAO//gpVq1p1XvgICXHiF0lK2qX3i32m0KHFznvF5hdAY2PHASic2VyZDGY/PPMqex4Y\nxqnUeKZ2GuHqcsq8QU2uZ3DsUzy+J4SbgmtRq6irsb/8Mlx1lfVDtmlT5xYpeZimycObfuKtk0n0\nNXfxxfVj8PWwYU/lyZPwv/9Zi0EvWwadO0NUFDzwALRsCUXcnsnN4cZdrW5jRMthfLhlLo9u/5Fm\n51K5sao/0yM6sv77NcX7oVStGlx9tfU4LycHYmOtRZV/+QVeecX6uHJl6NrVCphdu1r/BjzK0BBy\nBePMMOXMfWydzZkBVEpXmQxmTY7Gc0tQFcJObySixnBXl1MuzLz2Gep8+Qw3bDD4pWtU0YY0a9eG\nCRPg73+37prVHJ0SEZeaRI/VCzmUeo6ZDetwZ5MnXF1SXqmpVhibNQt+/hn69oXbbrOe/zEk6SwO\nw8HfWt3K8IghTNv4ERN/X878+HPU2Pw7h+Mm5Tm32D+UHA5o1Mh6DBtmHTNN+P13WL0aVq2C996D\nvXutnrTzQa1LF6hRo5ifqTjD2fQkUkLPQo9lUPMMBJ8DvwzwzWZDdTfqLnqLNNNBBg5MjD8ewB9/\nOjDxJAdPw8TLMMl6IAjP654l41RdiPeG05WoRgxdhrXnYPIZwnyruOzGGzv35lUU53vtC6tMBrOj\nQVWIb1SXBh9sBOUyp/D18OWbbrfRbdNWpu7bziMNIorW0EMPwccfw8yZcPfdzi1SmH94O0N37KRe\nRhwHeo4gtLKN9oT97Td4802YOxc6dIARI2DOHKtHqYR5uHnwaId7ePCqVP617l3+U+80fL4EZteC\nuS0h2+qZc/oPJcOwbhho3Bjuuss6lpgI69ZZQe3NN63jQUF/9sBdfbU1d02buZeIkynxLDu+i1Wn\nDrH5XAJ7MzI5keNJilsgpkcAjrt7YiTtw0w0INGEc574ZMfSvkowEdVrUtXDi0APb9wcBm4YGIaB\nAwOHYZCRk0VCRjrnstI5m5XBOf9M9pgGuz32k1Hfi2wfb1L8mzPZwyB69RoAvLIS8Ced6m451Pby\npHnlADpXDeOa0MYEeZXcwtx27s2rCPIOJU++4vkXKpOT/wf/fQLzW3vRY046y2ImXflNUmCjVkzn\n7Yz67Lv6Gmp5F/Euzd9+s4artmxRT4ETPfzLN7xxJod7Kyfydufh9lgCIz0dPvvMCiCHD8P991uB\nvFYtl5Z1Tf/H+dE/BIY4wL0hvNMYFjcmqu+E0p+UnZNj3diwcuWfj/h4qyetWzcrqHXocMX5dXKx\ng+eO8dWhTfx48jCbU9I4YlYiwysUz+wkqpFMHQ+D5pX96FAlhKuD6tPELwgvNzcWLlzO9Onfk5bm\nhrd3NqNGXev04T3TNDmYdJp1pw+wJfE4O5Pi2Z+WyuFMkzP4kO4RhFtOKgHZZ6npnkNjXx/aBwZz\nfc1mtAgILva/7/xvcniaqVP7aSjzMpx1w0TeG0PK+F2ZhmH0A14F3IB3TdN84S+vmz7/m0/qKyZR\nXr/Y4s6X8iQ7J5v6/5uMf5UWbOlxY9H/c3j8cTh61Oo9k2JJzkzj6mWfsC3Lh9lNwhlav6OrS7KC\nxVtvwfTp1lyxhx6C/v3B3R6d8Lk/lI6OhTvehQFBuKV787hvFZ6/foCry7P+baxaZQ31rlxp3c3c\nsmXeXrXgYFdXaSs5Zg4rj21n9oFNLEuIZ29OJTK8alI5O556bpm08/OnX0gD+oc2wd/D/mvUZWRn\nsurEHlac3M/GxJP8nprK/jSTZK8a4PDAPfEIdd2y6FW7Fr2r1+O6mk0J9PQu1DVKI4CWJ868Yzcy\nMpply6L/eFaGg5lhGG7ALqAPEAesB24zTXPHBeeYvrNfI2TiIaa92l/fZCVgd/x+mq/5mcnhjfi/\nxp2K1khyMrRoAW+/bc0xkiLZe+4E7X7+Bk/DZN3Vg6jr59q9LTl6FF54AT76CAYOhEcftQKFDV34\nQ8n0j+fsLR5srt6Wul7evNeyM72Carq6xD+lpFjDn+d71FavznsDQgUc/jRNkzXHtvPO3g38lHiW\nQ0ZV8KxCaM4ZOlb24eaajbmxdnN83MrHnKk/Q8FzUHcvXLUTn7Yr8GsRQqJPVdI9a+CVfY5gUmjm\n60m3qqEMrNmMlgE17NF7Xg44c/mT4vSY2ePX2z91BPaYprkfwDCMOcAgYMeFJ1XZuZtprw5RKCsh\njavUY0LwBsYdOMattROp6xtQ+EYqVbJ6VO6/HzZvhsKuMyX8cGQH123eRCtPk1WRd+Hp5sJ/rseP\nW4Hsgw9g5EhrmNrFw5VXMmBAj4v+j9h9Zi8j18yhz0Zo7mXy0VU9aRPg4rAL1tp/kZHWA/IOf/78\ns/W1j4+3biY4H9TK4fDn0aSTvLV7GfNPHGFHdiWyvEOpneNGr+p1uaNOC66pUR+3chpC8txFeSAc\nDoST+vUAekQ9w6JFj5CYnsR3cVv44cR+Np5N4OWzZ5lw+DSGw5sqOfGEe0In/6pEhTTkmuBwvAsS\nWHNyrOkI5x9ZWdYafzk5V/7TNK1fFPJ7GEb+x93c/tyRw9PT+tjd3TY3ijnzhonRo/sSGzvuojtj\nC8JuwawWcOiC54eBi7pstj79HIHefqVWVEU0oe3NfPrddPqsms/vfe4qWiNRUdCnjzWs+fbbzi2w\nnHtt1zLGHIzn1gCD2V3ucV0hJ09a+6G+9x4MH24NuYWW3a3PGldtwKr+T7Px+HZGrv+adus86Oyd\nycdtr6VBJRv98uBwWD3OLVpYS4vAn8OfK1da/6a2brV6K9u3h3btrEfz5rYZTi6oLaf28NLO5XyX\ncI6T3g0INLPoWDmEx8OacUvt5niVkx6xK7lsKEhLI+BUAree9eHW5BqQ4gNJSeScO87O+GOszDzH\nrx7ubPEL5OugExyrto+wU8doeeggbQ4dpv3BONofOEiNU6fzBrHMTCsgnd89w8PjzwB1pT/BCmc5\nOfk/8nstK+vPXTkyMqxHdvafIe3CwHa5j8/Xe+GjqMcueF7DSMSPs6ThTSYegBUYi3LDxPlfCqdP\nf4bFiwv3XrsNZQ4B+pmmed8fz4cDnUzTHHXBOaadai7PjibHU+fnxfxfSADPtb6uaI0kJlo/PN59\nV0OaBfSP9Z/zdoIHk0Mr8VSLPq4pIikJXnwRXn/dWuriqads30NWFMsOb+TuX5ey36c51/pm8kHb\nKEK8bbgeXH5SUmD9emtNtfOPQ4esf2/ng1qbNtaagjbqWTNNkx8P/8Kre9azLCmbJN8G1DITGFSt\nGo826UJ9XxsF5JKQnm4tYHzqVJ7Hh698RkLsVQRxKvdRjdMEux3Bx82w7uw9v9erv791t/MlHsfc\nc/iBRJaRzkaHwT5PPxIq18LDzCIk+yxNvdzoXDWYAXUiaFslxLW9kDk5F4e18x/nd+zCYHnhIy3t\n8s8LcE7GuSQyzqXiaZq4k0UGnmQ6wNPPBy9/vyuGvFRPd/ZU9uL3Sp7s8fVin683BytV4tvRE8v0\nHLPOQLRpmv3+eP4UkHPhDQCGYZgTJ07MfU9kZCSR57v/xemm71rGIwfOsKlDR1pUKeIP5iVL4L77\nrCHNgCIMi1YQ2TnZ9I75kFVZ/nzerDGDwlq5oIhs+PBDeOYZazeHyZOtLZHKuf/tW8U/tq7hmE9j\n+vik8W7bKMJ8ymCv/Nmz1l3R54Papk2wZ4/1d3i+By4iwvqzUaNSWwg3KyeLebEreHPfFtZleJLp\nU4eGJHBbaBijG3akqmcRF7V2tYyMP/dg/UvQyvfYqVNWCKhWzQpa58NWUBC/J6Tw2dKj7Dx9F6cI\n4jTVqFRnJk++NJB+N0cVe7gvOSOZpUe2svhYLBvOxrMnAxLcq4FHIIHZ8TTwNGldOYBOVWsSWSOc\nhpX8cNhkiLE0nZ+bmp7qoLJnBg/f242oyPakJZ9j+5nDbIs/xu7keA6kpRGXk8Nxw50zHpU46x1A\nqk8gldLO4btuDe6btuCfmUlgZiarv19epoOZO9bk/2uAI8A68pn8b6eaK4IuP85kX3o2R6LuxlHU\nycf33291bc+Y4dziyomkzFRa/zCT4/izpnNvWgS6YGL6jz9ak/krV7ZWte/QofRrcLEv965i7I51\nHPJqSA/PJN5vG0WDylVcXVbxZGRYi+Bu3WoNRW/daj0OHrQCW3i4tSdpePifH9evb/UGFENSehIz\ndv3Ih3F72JpTBcOrOq3dk7k7rBF/q9fafpP2MzP/3OD+UqHqr8dTUv7sxToftP4SuC465u9/yZBV\n2ndRZudk8+up3/kmbger4o+xOy2DkzlepHgEYXj4UTk7kRBHJg29vWjtX42rAkNoV7U29Xwq4V6M\nG1GcuYenM9pKzcrkt/jDbIo/wq5zp4lNOcvh9HROZJkk4kGqoxLZ7gE4spPxyU4iwMighrtBHS9v\nGlbyp7l/EG2q1KZFQAie+Xxfl/lNzA3DuI4/l8t4zzTN5//yuoJZKUvKTKPGD18w2N/gk663F62R\ns2ehVStraGyADZYrsJGD545x1bIv8PQMYHP3wdTwKfkFWfPYuxfGjrV+WL/4Itx0k20m47rKokMb\nGbX1Z2I9GtDRPYHXWkXSvmptp7Vvi82l09Jg3z5rm6k9e6w/z3988CBUqWINX9esefGfoaFWyKhS\nxbpx4Y/vl2NJJ5i6YylzTxxln6MmRo47fnv3U29bCs/eEMnAAT2LXG6Bv2ZZWZCQYIWsgjzOh61z\n56y9UAsSrs5/HBBQLv+tZGZnsuV0LMtPxLI+/jg7UpI4nAmJ+JDhEQgegXhlJ+NPKjXcTOp4eVHX\npxK1ff1oUKkKDSsHUa9SINU8PC7qdXPmkhSXamvK1L5079uF/cln2JV4gr3JZziQksjhtBROZGRw\nOodC47MAABfASURBVCuHRNNBCp6kO3zJcfPFkZWIT3YSgUYmNdwNanl50cDXj2Z+QbSuUpPWVWrj\n6160ZVjKfDC7EgUz1/juyA4GbPudeQ2DGVK/iEtoLFsGt95q7S9YhieQO9OaY1vptWEl4ZUC2dBj\nCN6leedlWpo1sX/qVHjsMSucFXWf1HJq+dHtPLw5hq2OMOqap5nYMIKRDYrXk+jMH0wlJivLuhM3\nLg6O/H979x1fdXX/cfx1bkLGTSAQICBDRpQQ9hYFjShDSiu21kFbq2JdVazS1lboT2qto66qtI+W\nVmrdoyiOQBnKtIgM2bKHyEoQSLjZyb3n98f3Kigz4Y5vct/Px+P7yOUm+ebc4/HmnfM93/PZ43w8\n+vHevV8Hm7L4eKaf35f/nH8es3r3o3Gxj0GfrSUzdzXeVedQRjKlJONNX8jwH3ShU89OTpg50VFV\ndWRdUfDY8tkW5s7cTFlhX1Ipoj4+MpJX0KWtl/R6HidUFRU5R0WFU/orPf3UR6NGR8JWw4YxtR1J\nTVUFqthRuIvlBz5nVUEeG4sL2FFaypd+i8/GUUICFR4vtl4axKcQ7y8h0ZaTQBWJBCjIP0TZoUZQ\nbqDUA6VxUOUho+lnXDiwMx5j8ECw2oJzoaUs4Kc04Kcs4Kc8EHAOa9m2cz8lgSaQ7IHEeEhMgIQk\niE8EfznGX0yivxgvFTTw+EmP85CRUI8WiV7O9tYnMyWdrAYZdGnUiuT48L33KZhJWEybtoBbV8xl\nT5ezuHDyRu69bWTNfoncfz8sXgwzZsT8m+Drm2bzky17uKRROjP6fzey6zlmzXI2he3aFZ5+OibW\nkZ2JHb58xqyYzn9LvXhtOTdnpPFQ9+EkxVV/jVYo90oKtdOZlQrYAPO+WMbELUuYW1TBYW8HWlT5\nGGkMdyY2IbuwmEf/73k2r/leMJIdOTqePYfvf6eX89v2RMfRd+EFj3+/Op91W0ZQTiJFpOKjPkWk\n0qnf6zz5j7HO5ff69Z2Pycl1chYr2qo7y1tcUcye4nw2H87n8+JDHKwopaCynJenLmJfYTdIrnKO\npADEQ0pqHu0ym2GBgAULX9cpTTCQ6PGQ6DEkGg9JcXEke+JYOGcd+Vv7QaEXClKhoAEcTOP87H+y\ncObDxHnccbm8usEMa22tOpwmSyTl5s63mZnjLAQsz71h+fPDtn3mfTY3d371T1ZZae0FF1j7+OOh\nb2gt8ocl/7KeGW/YW1bOs4FAIHI/eNcua6+6ytp27azNzY3cz60jSivL7S+Xv2MbTH/eema+ZS/4\n8Dn7wZ511TpHTs6E4yaSnJwJ4Wn0aTry//mRNmVmjrO5ufNtcUWxnbT2Pdv/v4/bhHcnWs8H023m\nh6/bX6+Zb3eXlhxzrlC/Rrf2Waw42diorqFDxx/3v+WwYb+L6rnCKZhbTjvnxPaUhZyWIxsfGrh3\nBJzbiW0XtWDixNnVP1l8PLzyirOWafnykLfV7fwBPz+Z/SAPFKTzx/YdmNQ9JzK7dgcCzl5yPXo4\nBbfXrtVavxpIik/giV4jKRx+A29nn0OZ9TBk7RYaTn+O6xe9ytbDe095DrcWl/7GBqcAjbeytXsH\nrtqVS/3Z/+HOfQbrbcuzXXPwXTyULZdcw2NdLqLFcWrqhvo1urXPYsUxYwPYuvWhGv0OuOuuoWRm\njv/Gc5mZ4xgzZkhUz+UmtWsnQomKb2x8eDAFHuoM9wf4/OVVNTth27ZOjcVRo2DZspipClBcUcyQ\n3PtYljaMF7I78eMW7SLzg7dscbYrKSmBuXOdrRLkjI1s3Y2RrbvhqyzngXUf8lr+fl5c8ilNK3Zy\nRXoDxmZdSMdGx14iPt6O4M4vk8si2fxj+PyVkPM2DNwKneKgaRYUpOLdmcE7fS5lcNOWp325PdSv\n0a19FitCuSP+0RuvHrn7tGaF1UN5LjfRGjM5peOuiblrASZnBWsvHEanph1rduKf/9xZRPzWW3V+\nvdnuw7sZOG0C+5pdyfQefRmUHoEyQH6/s37skUecDWLvvvvIjt0SFrtLfDy0fgFTDxxiX1xjUss+\np3+S5brWHbmybT9SElKA6BeXttay5eBWpny+lJn7d7O6LMChhLMhkAobymFRG5jbFQ4l1njtW6hf\nY7T7LJa5eV1kbaDF/xJyx7uLrP0546h4thclvo/Z8r3f0Si5Bns9VVRATg5873swblwIW+wuK/et\n5OIP/gKtrmJhnwF0TY3Adhjr1sHo0c5WBv/8p7M3lURUQUU5f9u6lDf2fc76qhQqrZ/08l309MYx\npElLRrTsQnbjTDzmzP4oOdWi7Ep/JRsPbGLOvg189OUuVhcXsz2QTGVKB7ymiux6fi5t0owW28t4\n5hdz2PatWalnnqn9MxByZo5/J7HGxulSMJOwON5fq70Hn8c5i+aQ9eV7fHLlROI9Nbgyvns39Ovn\n1GK8rO5dlnhvw/tcu+JD0lt9l//1HUibM9y085QqKuDRR51LxX/8o3MJs47PRtYG1lrW+Q7y4uer\nmHUgn22V8RTFNcCU7qKx/xAt61kyk710Tm1Ez0ZnkZnahIyUpqQnp1PvJHd+5ubOZ8yv32NH3h2Q\nuh+a7Set+1tkDWpJaZqXXVUeCuIagbcNSVTQ0lNBZ6+X7zRvx+XN29P8W9ujaFZKTkRjo+YUzCSi\n5h06wLBPP+ZHVYt5ftixU92nZeFC+OEP4eOPoX370DYwSqy1PLX4Gf5vz2EyW+Qwr88FNA53+Ztl\ny5xZstat4e9/dz6Ka5X6/Sw6uI8Z+7awzneI7eUV5Pk9+IyXKk8ipqoIW1GAx1+Chyo8NoAHP9aC\n38QRII6AiYf4BhCfBvFJUF4GBWWklueRk92GrvXTGdysPX3SGpNWy4qbi9QVCmYScU/t2Mx9G5bw\nYMND3Nv/zpqdZOJE+Mc/4KOPan09zapAFT//71hes1n0a9GH97v3wRvOtV2lpTBhglPj8qmn4Ec/\n0j5OtVxVIMCBqiryK8rJKyumpKqSIn85JVWV1PN48MYnkBqXyH1jX2HV/DugsB4UxUPwvT8n5/fM\nm/f76L4IEQGqH8z0J5ScsXvanMNSXyETts7l7LVvcm2Xq6t/kjvvhA0bnJmzadOcjSVroYKyAn7w\n9o182uQqRrbswr+zu5xRTblTWrAAfvYz6NkT1qyBjIzw/SyJmHiPh2YJCTRLSKBr6omLqT+9vxB2\neY95XttIiNReWnwiZ8wYw78796LDWf25adVc5u2YV5OTOKWBkpOddVG1cFZ084HN9H5hBCuaj+b2\nzH681Klr+EKZz+fs3D9qlLMn3BtvKJTFoLq6j5NILNOlTAmZveXl9FyyiJLNf+GjERPo1qxb9U9S\nXAyDBsHw4fDAA6FvZJjM3jqba2aMx3Z5iD9kZjOmVegKXh9jxgy49VYYPBieeMKp9ycxS4uyRdxN\na8wkqj4rLuaCZYuJ3/QnFl7xDNlNs6t/kvx8OP98uPdeJ4C4mLWWiUsmMmH1dMj6LZM6duLqcM1c\n7d/vFBr/3/+c9XhDNCsiIuJ21Q1mupQpIdUpJYX3u/emosNvyJlyE+v3r6/+STIyYOZMZ7uHyZND\n38gQqfBXcEvurTy6fRMJ2eOZ1r1neEKZtfDyy07B8ebNnXJKCmUiInWSZswkLKbk53PL+tXUW/NL\n5l3z+jdmzk61IebXNm+GSy5xLmmOHh3B1p/arsO7uHrKKPY2u5LkJv2Z3q07bZOPrRl4xnbsgNtu\ncyokTJ4MffqE/meIiEjY6K5McYUfZmSwr6IjD5onGPT61cy59g06Ne103B2kt251Fi8fE87OPRc+\n/NAJZx4P3HBDBF/Bic3eOpufvH8baT2eJKtRO97s3JkGod4jyu8/sknsL38Jv/oVhHsfNBERiTrN\nmElY/fmLL3h0+0b8K+/ircv/xsO3zq5+zbWNG+HSS516j3fcEeYWn5g/4OfhhQ8zcc3beHs+xeUZ\nrXgqMzP0d14uX+7UEU1OdtaSdegQ2vOLiEjEaMZMXOWe4O7zj3v+yg/evY2MtE7H/bqyspNswJqV\n5VQHGD4ctm93toeIcJmhnYU7+enUn7I/uQOBHs/ym7btub1ly9D+kAMHYPx4eOcdeOghuPFGlVMS\nEYkxeteXsLundWt+3eYckvtMYlvX1TDwUeCbs56n3BCzXTtYtAiWLoVrr4WysvA1+FteXfMqvf/R\nD9PuFnxn30hu126hDWV+v1NCKTvb2Vh3/Xq46SaFMhGRGKQZM4mIe1q3JtHjYfyASQSq7qfqrE/h\nvX9CeVpwQ8zTKGCeng6zZjlrzS69FN58E0I8a3X0jQmkFmKHrWZPQgntcv5DfGIDlmdn0zSUVQkW\nLYIxYyAlBWbPhu7dQ3duERGpdbTGTCJq2oED/GjVaup9/BLFvnfosfn7/G709dXbEDMQgEcfhWef\nde5UHDEiNG376saEbQ9Cj+dh8H14PZdTb+iPuad9O37Xpg1xoapBuXGjs2Zu6VLntai+pYhInaR9\nzMTVRjRuzPy+fUi+5GauvPIltvR/nw2NluIPVKO2n8cD48bBlCnOIvmxY6Gi4ozb9uyzs9jqGwU3\n5MB5L0DKVEq63EzHl5YxoW3b0ISyPXvg9tth4EDo3x82bYIf/1ihTEREAAUziYIe9euzuFcvNnqa\n0+vS95myeSYD/jWA1Xmrq3eigQNhxQrYtg369nV2xK+hnYU7Wd1uGlx/CQRuhov/BL4mcHNvkrYX\n1/i8X/viC6dQe5cu4PU6Bdvvvde581JERCRIwUyiomViIh/17EnH+uns7fAAA7reyeAXB3Nb7m3k\nFeWd/onS02HqVOey4DXXOOvP8k7/+/OK8hg7cyw9J/UkPrElpM+AS7LgsY7wZBaUxZ/6xoST2bjR\nKSvVvbsTyNavhyefhMaNa35OERGpsxTMJGoSPR6eOfdc/tahA28GMhkx4gM89dLI/ms2Y2eOZffh\n3ad3ImOcOzXXr4emTZ1ZqUcegcLCE37L0t1LuW7qdXT8a0d8/gCjr5jP4St+RfreT+FnfeFTpzC4\nc2NCNcsfBQIwbRpcdhlcdJFTYmrTJmebj2bNqncuERGJKVr8L65QUFnJuO3beXv/fu5u3pDdm5/j\n5VX/ZnD7wYzuMZrB7QdTL+74O99/u8TTuO93IOd/s2H6dGfbibvvhhYt2FGwg3c3vMtra18jrziP\nm/qMoazZd5iUd4DvpKfz+7Zt+WzOUiZOnE1ZWRxJSX7GjBly+jcmbNwIr7ziHA0bwl13ObN4SUkh\n7CkREalNqrv4X8FMXOVTn4/fbtvG1tJS7m6RQdXeGfxn7ctsOrCJoZlDGdJ+COe1Oo+sxlnEeeKO\nW+IpM3M8zzwzjPPbp+J75AEyps7iow6JvNjV4h1xBX07X8OWxCye37ePwY0acX/btmR5vdVvbCAA\nK1c6W3hMmeIs7L/2WucOy969taBfREQUzKRumF9QwGM7d7LM5+O6Zs0YnGL5Ys985uz4kGV7lrH7\n8G7aNGxD/pZyDu7pBVVJ4KmCBB+k5lGv6Wd4UxPodVYvLk7rznc3JbF+52Emd8hi7TnncJ3Pxy2t\nWpGVne2sUzsdBw/C6tXOsXgxfPCB871Dh8Lll8OgQRB3kgoGIiIScxTMpE7ZVFLCC/v28UpeHsYY\nLktP54IGDchOisNTlscNt/+JNZtGQnwpBOpBRSqUNKVXv1z+MOkXfOLzsbCwkOU+HxempXF9Sgoj\nly4lccECWLLEWfuVkODUo2ze3LnsmJQEiYng88H+/fDll7B3LxQVQbduztG7NwwZAmefHe0uEhER\nF1MwkzrJWsu64mJmHTrEJ4cPs6KoiJ1lZdiiMioONoQKD8Rb8FZBw0oSS0u5sM1Z9KtfnwFpaeQ0\nbEjK8WazrIX8fGd92JdfOqWevjrq13duJmjSxFm036qVLk+KiEi1KJhJzPBby6v/XcD4hz7ii7xf\nQJWBkjjapz/Is08Oq141ARERkTBQMJOYM23agprfSSkiIhJGCmYiIiIiLqFamSIiIiK1lIKZiIiI\niEsomImIiIi4hIKZiIiIiEsomImIiIi4hIKZiIiIiEsomImIiIi4hIKZiIiIiEsomImIiIi4hIKZ\niIiIiEsomImIiIi4hIKZiIiIiEsomImIiIi4hIKZiIiIiEsomImIiIi4hIKZiIiIiEsomImIiIi4\nhIKZiIiIiEsomImIiIi4hIKZiIiIiEsomImIiIi4hIKZiIiIiEsomImIiIi4hIKZiIiIiEsomImI\niIi4hGuCmTHmcWPMemPMKmPM28aYtGi3SURERCSSXBPMgFlAZ2ttd2ATcF+U2yNB8+bNi3YTYo76\nPPLU55GnPo889bn7uSaYWWtnW2sDwX9+ArSKZnvkCP2PHHnq88hTn0ee+jzy1Ofu55pg9i2jgenR\nboSIiIhIJMVH8ocZY2YDzY/zqXHW2veDXzMeqLDWvhrJtomIiIhEm7HWRrsNXzPG3ADcDFxqrS07\nwde4p8EiIiIip2CtNaf7tRGdMTsZY8xlwK+BnBOFMqjeixMRERGpTVwzY2aM2QwkAAeDT31srf15\nFJskIiIiElGuCWYiIiIisc6td2UewxhzmTFmgzFmszHmN9FuTywwxuwwxqw2xqwwxiyJdnvqImPM\nv4wxecaYNUc9l26MmW2M2WSMmWWMaRjNNtY1J+jz3xtjdgXH+org0goJEWNMa2PMXGPMOmPMWmPM\nXcHnNdbD5CR9rrEeJsaYJGPMJ8aYlcaYz4wxjwSfr9Y4rxUzZsaYOGAjMBjYDSwFRllr10e1YXWc\nMWY70Ntae/CUXyw1Yoy5ECgCXrTWdg0+9xjwpbX2seAfIY2stb+NZjvrkhP0+QTAZ619KqqNq6OM\nMc2B5tbalcaYVGA5cAVwIxrrYXGSPr8ajfWwMcZ4rbUlxph44CPgV8DlVGOc15YZs37AFmvtDmtt\nJfA6MDLKbYoVutkijKy1C4FD33r6cuCF4OMXcN5MJURO0OegsR421tp91tqVwcdFwHqgJRrrYXOS\nPgeN9bCx1pYEHyYAcTjvNdUa57UlmLUEvjjq37s4MsAkfCzwgTFmmTHm5mg3JoY0s9bmBR/nAc2i\n2ZgYMiZYq3eyLqmFjzGmLdATp8KLxnoEHNXni4NPaayHiTHGY4xZiTOe51pr11HNcV5bgpn7r7fW\nTQOstT2B4cAdwUtAEkHWWWug8R9+fwPaAT2AvcCT0W1O3RS8pPYW8Atrre/oz2msh0ewz6fg9HkR\nGuthZa0NWGt74JSVvMgYM+hbnz/lOK8twWw30Pqof7fGmTWTMLLW7g1+3A9MxbmkLOGXF1wfgjHm\nLCA/yu2p86y1+TYIeA6N9ZAzxtTDCWUvWWvfCT6tsR5GR/X5y1/1ucZ6ZFhrC4FpQG+qOc5rSzBb\nBpxrjGlrjEkArgHei3Kb6jRjjNcYUz/4OAUYCqw5+XdJiLwHXB98fD3wzkm+VkIg+Gb5le+jsR5S\nxhgDTAY+s9Y+fdSnNNbD5ER9rrEePsaYJl9dGjbGJANDgBVUc5zXirsyAYwxw4GncRbTTbbWPhLl\nJtVpxph2OLNk4FSIeEV9HnrGmNeAHKAJztqD+4F3gTeBs4EdwNXW2oJotbGuOU6fTwAuxrm0Y4Ht\nwK1HrQmRM2SMGQgsAFZz5DLOfcASNNbD4gR9Pg4YhcZ6WBhjuuIs7vcEj5estY8bY9KpxjivNcFM\nREREpK6rLZcyRUREROo8BTMRERERl1AwExEREXEJBTMRERERl1AwExEREXEJBTMRERERl1AwExER\nEXEJBTMRERERl1AwE5GYZ4yJM8ZsMMa0iHZbRCS2KZiJiDiFhtOttXui3RARiW0KZiIiMAiYE+1G\niIioVqaIxCxjzBU4Bc1HAcuAzcDfrbUbo9owEYlZCmYiEtOMMQnAQaCXtXZTtNsjIrFNlzJFJNYN\nAA4rlImIGyiYiUisGwLMi3YjRERAwUxEZDDBYGaMGWiMSYxuc0QklimYiUis6wx8EgxkA6y15dFu\nkIjELi3+F5GYZox5AqgE9gOTrLXFUW6SiMQwBTMRERERl9ClTBERERGXUDATERERcQkFMxERERGX\nUDATERERcQkFMxERERGXUDATERERcQkFMxERERGXUDATERERcQkFMxERERGX+H+LgWm4KdUEXQAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106d622e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t_train, y_train, 'o', label='data')\n",
    "plt.plot(t_test, y_test, label='true')\n",
    "plt.plot(t_test, y_lsq, label='lsq')\n",
    "plt.plot(t_test, y_robust, label='robust lsq')\n",
    "plt.xlabel('$t$')\n",
    "plt.ylabel('$y$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We clearly see that standard least squares react significantly on outliers and give very biased solution, whereas robust least squares (with soft l1-loss) give the solution very close to the true model."
   ]
  },
  {
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}