{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Least-squares regression with MOSEK\n",
    "\n",
    "Regression belongs to the basic toolbox of statistics, finance, machine learning, biology and many other disciplines that involves constructing approximate models of data. In this notebook we show how to quickly construct some and solve some standard regression models using the *Fusion* interface for Python.\n",
    "\n",
    "Table of contents:\n",
    "\n",
    "* least-squares (LSE, mean) regression\n",
    "* regularization (ridge and lasso)\n",
    "* robust regression with Huber loss function\n",
    "\n",
    "Other related material:\n",
    "\n",
    "* [Regression techniques for portfolio optimization using MOSEK](https://arxiv.org/abs/1310.3397)\n",
    "* Another MOSEK notebook on regression under Lp-norms\n",
    "* Future MOSEK notebooks will discuss LAD (least absolute deviation, $L_1$ or median) regression, quantile regression and mixed-integer piecewise-linear regression models.\n",
    "\n",
    "A general linear regression model approximates the relation between a vector of $d$ independent variables (features)  $x=(x_1,\\ldots,x_d)$ and a dependent variable $y$. (To allow for intercept we are often going to assume $x_1=1$). Suppose we have $n$ input vectors $x$ arranged in the rows of a matrix $\\mathbf{X}\\in\\mathbb{R}^{n\\times d}$ and $n$ corresponding observations $y$ arranged in a vector $\\mathbf{y}\\in\\mathbb{R}^n$. The goal is to find a vector of weights $\\mathbf{w}\\in\\mathbb{R}^{d}$ so that in some approximate sense we can write\n",
    "\n",
    "$$\\mathbf{y}\\approx \\mathbf{X}\\mathbf{w}.$$\n",
    "\n",
    "We call $\\mathbf{r}=\\mathbf{y}-\\mathbf{X}\\mathbf{w}$ the *residual* and we aim to minimize some penalty function $\\phi(\\mathbf{r})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mosek.fusion import *\n",
    "import numpy, sys\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Least-squares regression\n",
    "\n",
    "In the simplest *unconstrained least squares problem* (LSE, mean regression) the penalty function is\n",
    "\n",
    "$$\\phi(\\mathbf{r})=\\|\\mathbf{r}\\|_2=\\sqrt{\\sum_{i=1}^n r_i^2}$$\n",
    "\n",
    "and solving the least-squares regression problem is equivalent to the minimization problem\n",
    "\n",
    "$$\\mathrm{min}_{\\mathbf{w}\\in\\mathbb{R}^d} \\|\\mathbf{y}-\\mathbf{X}\\mathbf{w}\\|_2.$$\n",
    "\n",
    "A conic formulation of this minimization problem is:\n",
    "\n",
    "$$\n",
    "\\begin{array}{rl}\n",
    "\\mathrm{minimize} &            t \\\\\n",
    "\\mathrm{s.t.}     &            (t, \\mathbf{y}-\\mathbf{X}\\mathbf{w}) \\in \\mathcal{Q} \\\\\n",
    "                  &            \\mathbf{w}\\in\\mathbb{R}^{d} \\\\\n",
    "                  &            t\\in\\mathbb{R} \n",
    "\\end{array}\n",
    "$$\n",
    "This formulation can be translated directly into a *Fusion* model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Least squares regresion\n",
    "def lse(X, y):\n",
    "    n, d = len(X), len(X[0])\n",
    "    M = Model(\"LSE\")\n",
    "    \n",
    "    # The regression coefficients\n",
    "    w = M.variable(\"w\", d)\n",
    "    \n",
    "    # The bound on the norm of the residual\n",
    "    t = M.variable(\"t\")\n",
    "    r = Expr.sub(y, Expr.mul(X,w))\n",
    "    # t \\geq |r|_2\n",
    "    M.constraint(Expr.vstack(t, r), Domain.inQCone())\n",
    "\n",
    "    M.objective(ObjectiveSense.Minimize, t)\n",
    "    return M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Least squares example\n",
    "As an example we use the basic model to estimate a simple polynomial regression problem with synthetic data points sampled from a degree $3$ planar curve with Gaussian error. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd8zdf/wPHXSSQae68YoaV2hdh7bwmqRiml6NBftaWlOtRo1SittlSHUq1RRRSltX2pESOorahZhKgQss7vj3ujuXFvcm/uTt7PxyMPd3zu/bw/V3Lfn88573OO0lojhBBCJPNxdwBCCCE8iyQGIYQQJiQxCCGEMCGJQQghhAlJDEIIIUxIYhBCCGFCEoMQQggTkhiEEEKYkMQghBDCRDZ3B5ARhQoV0kFBQe4OQwghvMrevXuva60Lp7edVyaGoKAgIiIi3B2GEEJ4FaXUOWu2k6YkIYQQJiQxCCGEMCGJQQghhAlJDEIIIUxIYhBCCGFCEoMQQggTkhiEEEKYkMQghBDChCQGIYQQJiQxCCGEMCGJQQghhAlJDEIIIUxIYhBCCGFCEoMQQggTdicGpVQBpdTvSqmTxn/zW9juI6XUYeNPzxSPl1VK7TK+frFSyt/emIQQQmScI64YRgEbtNblgQ3G+yaUUh2BmkANoC4wUimVx/j0R8B04+tvAoMcEJNwoBX7L9Jw0kbKjlpNw0kbWbH/ortDEkI4kSMSQygwz3h7HhBmZpvKwBatdYLW+g4QCbRTSimgBbA0ndcLN1mx/yKjlx3iYnQsGrgYHcvoZYckOQiRiTkiMRTVWl8GMP5bxMw2kUB7pVQOpVQhoDlQCigIRGutE4zbXQACHRCTcJAp644TG59o8lhsfCJT1h13U0RCCGezamlPpdR6oJiZp8ZY83qt9W9KqdrADuAa8AeQAChzm1uIYQgwBKB06dLW7FY4wKXoWJseF0J4P6sSg9a6laXnlFL/KKWKa60vK6WKA1ctvMdEYKLxNT8CJ4HrQD6lVDbjVUNJ4JKF188B5gCEhISYTR7C8UrkC+CimSRQIl+AG6IRQriCI5qSVgL9jbf7A+GpN1BK+SqlChpvVweqA79prTWwCXgyrdcL9xnZ9nEC/HxNHgvw82Vk28fdFJEQwtkckRgmAa2VUieB1sb7KKVClFJfG7fxA7YppY5gOOvvm6Jf4U3gNaXUKQx9Dt84ICbhIGHBgXzYrRqB+QJQQGC+AD7sVo2wYOkKEiKzUoaTdu8SEhKiIyIi3B2GEEJ4FaXUXq11SHrbychnIYQQJiQxCCGEMCGJQQghhAlJDEIIIUxIYhBCCGFCEoMQQggTkhiEEEKYsGpKDGG/FfsvMmXdcS5Fx1IiXwAj2z4ug8SEEB5JEoMLJE9dnTxLafLU1YAkByGEx5GmJBeQqauFEN5EEoMLyNTVQghvIonBBSxNUS1TVwshPJEkBheQqauFEN5EOp9dILmDWaqShBDeQBKDi4QFB0oiEEJ4BWlKEkIIYUISgxBCCBOSGIQQQpiQxCCEEMKEdD4Ll5C5ooTwHpIYhNPJXFFCeBdpShJOJ3NFCeFdJDEIp5O5ooTwLpIYhNPJXFFCeBdJDMLpMjJX1Ir9F2k4aSNlR62m4aSNrNh/0dlhCiGMpPNZOJ2tc0VJZ7UQ7mVXYlBKFQAWA0HAWeAprfVNM9t9BHQ03h2vtV5sfPw7oClwy/jcAK31AXtiEp7Jlrmi0uqslsQghPPZ25Q0CtigtS4PbDDeN6GU6gjUBGoAdYGRSqk8KTYZqbWuYfyRpCCks1oIN7M3MYQC84y35wFhZrapDGzRWidore8AkUA7O/cr0uHNbfTSWS2Ee9mbGIpqrS8DGP8tYmabSKC9UiqHUqoQ0BwoleL5iUqpg0qp6Uqp7JZ2pJQaopSKUEpFXLt2zc6wM7fkNvqL0bFo/muj95bkIAsbCeFe6SYGpdR6pdRhMz+h1uxAa/0bsAbYASwE/gASjE+PBioCtYECwJtpvM8crXWI1jqkcOHC1uw6y/L2AWVhwYF82K0agfkCUEBgvgA+7FZN+heEcJF0O5+11q0sPaeU+kcpVVxrfVkpVRy4auE9JgITja/5EThpfPyycZP7Sqm5wAgb4xdmZIY2elnYSAj3sbcpaSXQ33i7PxCeegOllK9SqqDxdnWgOvCb8X5x478KQ//EYTvj8RrO7AOQNnohhD3sHccwCViilBoE/A30AFBKhQDPa62fA/yAbYbvfv4F+mqtk5uSflBKFQYUcAB43s54vIKz6/RHtn3c5P3BcW30MkuqEJmf0lq7OwabhYSE6IiICHeHkWENJ23koplmncB8AWwf1cIh+3DGF3jqhAaGhCPt/0J4B6XUXq11SHrbychnN3BFH4Az2uhl4JkQWUOWTQzubBIpkS/A7BWDp/cBZIZObSE8kac10WbJSfTcXefvrXX6lhKXBh4dvYYgLxxMJ4S7ufv7yJwsecXgyCaRjGR6ayaV01qTGBdH/LFjJEZGku3CBfyiovC9cgWio+HePcOP1hAQYPjJmxeKFzf8lC4NlSpB+fKQ3eK4QZuY69ROlmjsq5IJ74SwjSc20WbJxOCoJhF7qovCggNpW7EABw4c4M8//2T7j7+y6K3DFDt5kopXr1IzJoZqWpPyHD0WuATczpaNRH9/eOQRHnnkEfL6+5PT15eccXH4Xb+Oik1xHD4+UKEC1K0L9eqxscBjvHNacenWPZsvWVMmNHNNYQ/ilH4HIazmiU20WTIxOKqN39ZMHxsby5YtW1i7di3bt2/nwIEDFEhIIAx4UimaaY0fEO/jw/nAQCJLl+ZGYCDRJUtys0ABbmpN7L173LhxgytXrnD58mXOnj3LP2fPPthHvrx5aVK7Nm0rVKBp0aJUSErC7/BhWLMG5s2jBVA1Z362lQ1mS9laTLhxE6hHWHCgVVc/yZ3aZUetJq16ttS/1J7WhiqEp/DEPscsmRgcVedvTaa/c+cO4eHhLFy4kA0bNhAbG0uO7NkZVqEC35UtS6XTp/FJSkI/+iiqe3do1w6/unUpFxBAOSvjuHXrFidOnCAyMpKIiAj27NnD/82dS2JiItmzZ6dRo0Z0GTOG9Uc0j58/S+Oz+2l+OoLuhzdyf0029qyqTUSfXoy7XZobyg9I/+rH0i9zyueTyfoKQljmzHFHGZVlxzE44gw2rfEI01vk4vPPP+fnn3/mzp07lCpViic7dWKQjw+VVq3C59w5KFIEnn0Wnn4aqlYFwyBAh/j333/Ztm0bGzZsYO3atRw9ehQA/6KPkqNSU3JVakSdf6/T4fh22h/fTonb17ntH8AvlZrwU7XW7C/xOChlcWyFuTENyVKPbXDFuA0hvFnK76O8AX4oBdF34x1+dW3tOIYsmxgcIfWXo05KJOHUDnKd/I0Th/eTO3duevXqRb8ePWh44AA+06bBP/9A/frw6qsQGgr+/i6J9eTJk7R+ZRqXD2wi7vIJQPFImWrkqt6WXOXrUefyCXocWk+H4/8jR/x9DhZ7jLm1urCmYmOOTzU3m/p/v8wXo2PxVYpErQk084tsqdlJAWcmdTTzjBBZk7MHkUpicJEV+y8y+dcjnNz5O3d2LiL22nnKly/Pyy+/zID+/cm9ejWMHg3nzkGrVjBmDDRt6tCrA2tiTNlhHH/jIneObCHm8AYSb/2DT4685K7emlzBHcmbPSddj2xmQMRKHrtxgeu5C1Do7Tfh+echT56H3tOaKy65YhDCOs7+W7E2MWTJcQyO9Mg/h7j1w3Cu/zKFR4vmZenSpRw7dowqRSpyrmo96NOHE3F+/O/LxfD779CsmcuTQnKNdDK/AoHka9SHwKFfUaTH+2QvWZlbu5ZxcfYgzq77nG+Klaf1c18wqPcEkqpUhTffhDJl4L33IDra5rprbx23IYSreUqFkiQGI1tnOz116hRdunShffv2xMfHs2TJEiIjI+nesSMnBw+ncZ/2FI26zIgOw2nXdxqDz+dxy4AVc5VTyZTyIaBcLYp0HUPNkQsI7TuY+39FcGX+q9xc9j6PtqxEkT+2wO7dhoQ2bhyUK8elt8ai794xeS9L6z0kX1nExifia0yIsr6CEOZ5yszIkhiwbeRhfHw8EydOpEqVKmzatInJkydz+PBhevTogc+BAxAczOPfziS8cjNaP/cFS6u1IsnH120L5VhzphHg58s7vZqwYv5sbly9xEcffUT2f88z5rnuNG/enO1xcbB8OezfD/Xr8+Lar9n65WB6Rq7DJ+m/pGOuRDXl1Uqi1g+uFFInBW9eilQIR/GUq2tJDFi/4tnBgwepV68eb7/9NqGhoZw4cYKRI0fi7+cHn3xi6FS+fZv+Pd5nRMdXiQ7IY/J6dwxYsXSm4auU2dXR8uTJwxtvvMGZM2eYMWMGx44do1GjRnTq1IlIpWD1al4YOoO/8xXjo7UzWTVvOPXPRZrdl7WfqydOCSCEO3jK6oVZchxDaum162mt+fjjjxk9ejT58+dn6dKl+Jarx5PzjhJzZSuf//YJjY7+AZ06wdy5RH6+D2LjH3q/vAF+FmNw1gAwSzXS6f2y5ciRg1deeYXnnnuOmTNn8tFHHxEcHEz//v1p3HsY/QpVoOWhzYza/B0LF43h10qNUR9PMzkWS+McUn/eGZkSQAbMiczKE1YvlCsG0m7Xi4qKokuXLowYMYKOHTvy559/4luuHqOXHSL76ZMs+/516hzfwwdthrLi/VlQqJDFvmVLjzvzjNneM5CcOXMyatQo/vrrL0aMGMGPP/7Iy12bUCNqAxG1mtLquVl81WoAbf7aTbsnm3P4tXd556f9Vg9+A9s73OQKQwjnksSA5Xa90MC71KhRg99++41PP/2UZcuWUahQIaasO07t47tZ8f3r5L0XQ5/eE5kT3Jkpv50ADANTzLH0uLVNLhkVFhzI9lEtODOpI9tHtcjQ2Uj+/PmZPHkyR48epUOHDiyaPY1/5r7E5KbZeO63b/E9cgSaNaPq9PEs/HY4Vf45bfZ9zLWX2trh5uzPS4isThID5s+q2wec4t3BT5ItWzZ27NjByy+/jHF5UppuXs7cpe9zIW8RQp+ZTkTJKsB/Z7i2ftF5SomaNcqVK8dPP/3Epk2byJEjB6GhoXTu3JkzSsEvv/B82FsUiblB+LxXeXPzd2SPv//gtZauVmztcPOmz0sIbySJwSj5rPr0B+1pHrOR6WNepm7duuzevZtatWoZNtIaJk7kg3WfsaVsTZ58ejIX8xZ58B7JX/y2ftF5SomaLZo1a8aBAweYOnUqW7ZsoWrVqkz7+GMOhjSj1XOz+LlqS17YtZRV84ZT7fLJBwN0zF2tJCfmfCn6YO4nJDJ88QGzFUre+HkJ4U0kMaQQFxdH3759mThxIs899xy///47hQsXNjyZlATDh8Pbb/N3x+680vM97vr/90WU8ovf1nZ9TylRs5Wfnx+vv/46R44coWXLlowYMYJ/FrzG7TvXebPDK/R7ahy57t9l+fev881fv0C8+aa0ZPcTkh7cTjIOyDfXf+Ctn5cQ3kKmxDCKiYmhe/fu/Pbbb0yaNIk33njjQdMRSUkwdCh8/TW89hpMmcKKyMsOrYrxhCobe2LQWrNs2TKGDRvGtevXKd60D77BXXk8hw9f7fue0quWQu3a8OOP8NhjD73e0lQAyVJPCeAJn5cQ5njy76bMlWSDqKgoOnbsyJ49e/jqq68YOHDgf08mJcELL8CcOfD224bRvy6c0sJVHDV5140bN3jllVdYsGAB1atXZ/78+TzxxBPw888weDDExcHMmTBggMnnmN76DjLhnvAGzp4Ez14yV5KVoqKiaNmyJQcOHGDZsmWmSUFreOklQ1J4661MmxTAcZU+BQoU4Pvvvyc8PJyrV69Su3ZtJk+eTGJYGBw8CHXqwMCB0KcP/Pvvg9el1z8g/QfCG2SWirksnRiSk8Lx48dZuXIloaGhphu88QbMng2jRsGECZk2KYDjK326dOnCoUOH6NKlC2+++SbNmzfnbEKCYSLBDz6An36CWrUM02xgvt8gmfQfCG+RWSrmsmxiiIqKolWrVhw7dozw8HDatGljusG0aTB1KgwbZvgiy8RJAZxT6VOoUCF++ukn5s+fT2RkJDVq1GDx0qWGacg3b4bYWKhXD2bNIqxGiQcd9oBMuCe8UmapmLO7j0Ep1QMYC1QC6mitzTb+K6XaAZ8AvsDXWutJxsfLAouAAsA+oJ/WOi6tfdrbxxATE0PLli2JjIxk5cqVDyeFBQugXz/o0QMWLgRf82eymYmz20bPnDnD008/zR9//MGAAQOYOXMmue7dM3zOa9dC377w5ZeQI0eaMXpqp54Q4KS/o9u3IXduh8Tnyj6Gw0A3YGsawfgCnwPtgcpAb6VUZePTHwHTtdblgZvAIAfEZFFcXBzdu3dn7969LFmy5OGksH69YbnN5s3h+++zRFIA50/eVbZsWbZu3co777zD/PnzCQkJ4eClS6yYMIevW/UnacEPnHq0Kr+HbzP7epkGQ3gDp/wddewIPXs6LEZrOKwqSSm1GRhh7opBKVUfGKu1bmu8P9r41CTgGlBMa52QejtLMnrFkJSUxNNPP82iRYt47MmRJDza1PTM8/hxqFsXSpWC//0P8ua1eR/myJmu6WeQ88ZxLi2fzJ3b0RRoOQT/qm1oemYfn/wyFV+dxJ9TvqD+y8+YvF5WgRNZ0qFDUL06evJk1MiRdr+dp1UlBQLnU9y/YHysIBCttU5I9bhTjB49mkWLFlG45SDiH21qcua5esuf0KUL+PnBL784NClk9TPd1J9BTIHHKdhvBtlLVuHKmplcXzWVzSWr0HnADM7nLUrdVwYY+nVSnLRklk49IWzy+efcV4pnNm926W6tSgxKqfVKqcNmfkLTf7XhLcw8ptN43FwMQ5RSEUqpiGvXrlm5W1N9+vShdJuB5AjpavJ43P04ijz3DJw5Y1iQJigoQ+9vTmYpX7OHuc8g3j83+buPJW/jvtw9spUr81/jTEIc3ftO5peKTQxrYz/1FNwxrBSXWTr1hLBadDSJ8+axQGvqtm/v0l1blRi01q201lXN/IRbuZ8LQKkU90sCl4DrQD6lVLZUj5uLYY7WOkRrHfJgmgobPfHEE/gEd3vo8VGb51L71D5DaWqjRhl6b0vkTNfysSrlQ74GvSjSczyJsbe4PP81ok7tYXK/d2DKFMOguMaN4fx5mQZDZD3z5uF77x7f587Ns88+69Jdu6opaQ9QXilVVinlD/QCVmpDB8cm4Enjdv0Ba5NNhqQ+w2x7fAeD96zg53qhhoFXTt5feo9nBqmX6cyXw/wCRfkC/Ajw8yUgqAbF+3+Cf6EyXF/5EYWP/kTC8OGwahWcPg21axN272+PWNlKCJdISiJuxgx2AA2GDSNnzpwu3b3diUEp1VUpdQGoD6xWSq0zPl5CKbUGwNiHMAxYBxwFlmit/zS+xZvAa0qpUxj6HL6xN6a0pDzzDLpxkSlrZhBZ4nH8Znzs9P0ly8xnuub6VGLuJeDna9pqGODny9guVR582WfLU4gST08iV3BHwr//kpoNmnGtdm3YuRNy5YLmzQk7ttXudSWE8Arr1+N/9iyzfX15+eWXXb57u5f21FovB5abefwS0CHF/TXAGjPb/QXUsTcOayV/mXyy6iBffPsh2teXK3Pm0aVuOafuL6tUJZntT0jS5AvwI2f2bBY/A0PtNxRs8wLZi5fn8LrPqfJETX5b8ws1du2Cbt0M02icOAHvvpvpBxyKrC1u+nSigex9+lC8eHGX7z/rTqL3/POGAVVr1oCLO3YyM0uT4aU1CZ65UtT7V05xY8VEfONimDt3Lj27doUhQ2DePEOC+PZbyJ7dYhxSIiy81qlT6AoVGK81XQ8epFq1ag57a08rV/UsK1caksLIkZIUHCwjfSrmOqezF3uMov0+pmbNmvTq1Yu3xo4l6ZtvYOJEw9TdrVtDVJTZ95MSYeHNEqZNI05rjjRt6tCkYIuslxguX4ZBgyA4GMaPd3c0mU5G+lQsJY1SgSXYuHEjQ4YM4cMPP6Rrt27cfvllwzQlu3ZBgwaGzulUpERYeK2bN9Fz5/IjMPS999wWRtZKDElJhukuYmLghx/SbIoQGZORKQHSSib+/v7Mnj2bmTNnsnr1aho0aMCZunVhwwa4fh3q14fdu01eKyXCwlslfvklfvfvs7FaNZo1a+a2OOzufPYqn30G69bBF19ApUrujsZjOLo9Piw40KbXp9dBr5Ri2LBhVKxYkR49elCnTh2WL19Ooz/+gHbtoFkzWLIEOnUCDFcg5qbPyMwlwiITiI/n3pQp7ASemjjxvxUk3SBrdT6HhxtGNs+dm+mrWqz9svf0FadSO3HiBJ07d+bs2bN89dVXPNO2rSEh7NsHs2bBkCFed0xCACT98AM+ffvyYpkyfPbXX/j4OL5BRzqfzQkNhe++yxJJwdrOV29rjz9yJyf5ek3Gp3hF+vfvT/fR40nasMFw5TB0KIwda7K2g6cPhks9GFA6yLMorbk1dizHgEYTJjglKdgiazUl2cCbyx3T+rJPfQze1B7/35VANor0GMeN32ezbO7nNDp/lg0/LyRg+HB4/324dImwL77w+P+v1Fc2yQkc8PjYhWPpTZvIf+oU0woVYmyvXu4ORxKDOd7+B2vLl703tcenTHjKNxsF2r6EX4ES/LFhLs1at2ZleDhFS5QwlLReuQKLF0OA5x1HMlsSuMjcokaOJB4oO3Ys2bK5/2s5azUlWcnbmldSs2UsgTdN2ZE6sSmlyFOnG4XD3uLw4cPUrVePI336wOefG+ZZatMGbt50U7Tp86arNeE8eu9eCu3bx7x8+eg3eLC7wwEkMZjl7X+wtnzZO3vlNluk195uKeGVr9OCrVu3cv/+fRo0aMDGihUNVwu7d0OTJnDRM9vts+IEi+JhV159lVtAsbFj8ff3d3c4gCQGs7z9D9bWL/uw4EC3T05nTYd5WgmvVq1a7Ny5k5IlS9K2bVu+u3MHfv0Vzp6Fhg0Ncyx5GG+6WhPOoU+coMi2bfyQJw99XnjB3eE8IInBjMzwB+sJX/a2sKb5Lr2EV6ZMGbZv306zZs149tlnGbt1K3rzZrh717DOxr59Ljyi9HnS1Zpwj/PDhxMP5H33XY+5WgDpfDYrq82ImhGOrtqytvkuvcFzm/6KIabZSHJdgffff5+t+46wduNG/Dt2NAyEW7nS8G8GuXswoMg89PnzFFu7lp9y5aKHG6bWToskBgvkD9YyZ1RtOaI6KmVcBTq8gm++omz65QdCrlzif7/+Sp4ePQzjHRYvNoxpsZG3V6sJz3Lm+ecppTXZ3n7bo64WQJqSRAY4o2rLEc13JuWsSpGvYW8KdnyVw3t30rBnTy4uXAhPPAHdu8P8+TbH6O3VasJzJJ4/T+Cvv7Iib166v/66u8N5iCQGYTNnVG05or3d3P5zVW1JkR7vc+7cOep26MDhTz4xNCX17w8zZtgUo7dXqwnPcXzQIHy1JvcHH3jEuIXUPC8i4fGcNSjO3uY7S3E9WqM+s976Hx06dKBBmzYsX7SIlvnywauvGsY5jB1r1TQp3jQYUNjHmTMf3D97lnLr17OmYEE6e1AlUkpyxSBs5qlVW2nFVb16dXbu3ElQUBDtQkNZ0KkTDBwI48bBK68YpmS34/1F5uHshZ7+HDCAbFpT8OOP3TqDalokMQibeWqZZXpxlSxZkm3bttG4cWP6PfssE8uWRb/6KsycaVinIyHBrvcXmYMz+5JunzhBxS1b2FisGA2feSbd7d01yWLWmnZbCCAuLo6BAwfyww8/MHTIED4vWRLfd9+FsDDD6nCPPOLuEIUbZWTdcmvtrFmTWvv3c2z5cqqFhaW5rTOmj5dptx1ApkTOnPz9/Zk/fz6jRo3iyzlz6LpnD3EffwwrVhjWdoiJcXeIwo2cNfPBha1bqbV/P1srVEg3KYB7q+AkMVggC8pnbj4+Pnz44Yd89tlnrFq1iiaLF/PvZ5/B5s3QqhXcuOHuEIWbOKsv6cwzzxAPVFywwKrt3VkFJ4nBAqlZzxpeeuklli1bRmRkJLVmzODyZ5/B/v2GktYrV9wdnnADZ/QlRS5YQMNz59hbvz6BtWtb9Rp3ztkmfQwWOLOdUbhf6nLEzsVimDpiID4+Pmx97z0qvvkmFCsG69dDUJC7wxVeTGvNtoIFqREdjc+ZM+QqU8aq13ltH4NSqodS6k+lVJJSyuLOlFLtlFLHlVKnlFKjUjz+nVLqjFLqgPGnhj3xOJK3z7AqLDPXTDjvtD/vf72cXLlyUXPkSLa99x5ERRlmZj1yxN0hCy+27p13aHLzJidDQ61OCuDeKji7rhiUUpWAJOBLYITW+qHTeKWUL3ACaA1cAPYAvbXWR5RS3wGrtNZLbdmvK64YZEH5zKvhpI1mB6oF5gtg2bNV6NixI/v372fRmDH0+OoriI+HtWshJN0TLSFM3Lp5k7+KFCFQKQpFReGTO7db43HJFYPW+qjWOr1G9zrAKa31X1rrOGARYPsMZi4mNeuZV1qdekWLFmXz5s20a9eOp8aPZ1pYGDp3bmjRArZscXGkwtut7NOH4IQE7rz1ltuTgi1cMSVGIHA+xf0LQN0U9ycqpd4FNgCjtNb3XRCTVWSG1cwpvaktcuXKRXh4OC+++CIjZs/mXLdufHLkCKpdO1iyBDp3dnXIwgsd3rOHJmvX8nfBgpR99113h2OTdK8YlFLrlVKHzfxYe9Zvbsx3cvvVaKAiUBsoALyZRhxDlFIRSqmIa9euWblrIR5mTTlitmzZ+PLLL5kwYQIzly2jW6FCRJV5jISwMIZ3HiHjWkSatNbs6NmTMkDer78GH+8qAE33ikFr3crOfVwASqW4XxK4ZHzvy8bH7iul5gIj0ohjDjAHDH0MdsYksjBrF2JSSjFmzBhKly7NswMHUjFvCcKLVeDjVR/z7v27jL7TxeT9Ukpd9dS8YmE2HbsmCz9lEUs+/ZReZ85wtnp1gqwYzOZpXNGUtAcor5QqC1wEegF9AJRSxbXWl5VhJqkw4LAL4hHCpmbCfv36MWnLPxz9/j1a38vOslJVmPD7LPLdu82UHAON8YYBAAAY/0lEQVQeeh9zC/os2Pn3g+dlgZ/M7erVq+g33iBAKUovXuzucDLE3nLVrkqpC0B9YLVSap3x8RJKqTUAWusEYBiwDjgKLNFa/2l8ix+UUoeAQ0AhYII98QjhLHcLVaJY38nc98lG5yunWFzmCUZsW8DAZTMfmpnV3ODI1GSwZOZgbtqcWX370isujluDBuFTsaK7Q8wQGeAmhBWSS1wTbkdxden7JFw9w+dlqvPCuUh45hn4+mvw8wMsD45MTQZLejdzJe2c3cuKxe9RNk8e8l66BDlzui9AM2QSPSEcKLnDOlvughTrM4ns5Wry4rlIfgiuY1gmtFs3uHsXsH4QpAyW9G6prwyT4mLpEj6VGkCOWbM8LinYQhKDEFZIOa7FN3sOggd+SNsn+9F3/26+qlkTvXo1tG0L0dFmq55SkwV+vF/q8TDZfv+S9+/dZmvRx/Dr3dtNUTmGLO0phJVSd1hr3Ypp06oz9I03OPfYY4zftQvVpAlh69ZBt2pSlZTJpRwPE3tmPzMOryeH8mFWnzE08dCV2awlfQxC2GnZsmX07duX7nnzMu/WLXyKFoV166BCBXeHJhwsZRly3gA/7sQlcP9uDDW+HMzK2H+Z1qAXMxv3JdBDk7/0MQjhIt26dWPLli2sB1r4+BB386Zh8r3du90dWpbm6IW2Uk++GB0bDxrUhtl8Fvsvx/IWY3aDnoD3r98iiUEIB6hduza7du3iRrlyVL99m3+1hubNDZPvCZdzxkJb5sqQbx3dzuuHNlESeKvz68T7+j14zptLkiUxCOEgpUuXZvv27TzWvj0VoqK4kCMHunNnmDfP3aFlOc5YaCt1Z3PCv1d5fPXHDAPmB3dgX2CldF/jLSQxCOFAuXPnJjw8nKdfe43K169zIE8eGDAAPvgAvLA/z1s5Y1nMlOXFOimR2PDJzI2/z6W8hfkh9AWzr9HglfNqSWIQwsF8fX2ZNm0a07/+mkb//ssvuXPDmDHw4ouQkODu8LIEZyy0lbIM+daORUy4dIxyaM5Mm81LXYItlih7Y3+DJAYhnGTQoEGs3biRQf7+zPD3h9mzISwMYmLcHZpNHN2J6wrWzKBrq+SxLLmijlF3+yJeAE73G0rjQd1MxrmY4239DVKuKoSTnT17ltDQUBodOsRMQNWsiVq1yrCmtIfz5pUMU5eWKgXRd+PtGkdy4cIFWtSowbZbtyhUvjy++/bBI4+YbGPPevGpZ+V1dMmrlKsK4SGCgoLYsWMH1558ki5ac//AAZLq1IHDnj+ZsDM6cV0lLDiQ7aNaML1nDe4nJHHzbrxdFUr379/nye7dmRYdTRGl8F2w4KGkABlvxnJGJVVGSWIQwgVy5szJ4sWLaTBxIo0SE4m6fJmkBg3g99/dEo+1zUPO6MR1NUclt9dee43au3fTOTER9dFHULOm2e0y2ozlSUlYEoMQLqKU4q233mLc6tU0Dwjg6N27JLVrB3PmuDQOW85MndGJ62qOSG7ffPMNf3zxBdN9faFTJxg+3OK2GV0v3pOSsMyVJIQTWGorXrH/IlMOPkJ072m0XD6RudfO0X7oUPSff6KmTYNszv+TTOvMNPWX18i2j5vtY/CkCQDfXnGIhbvOk6g1vkrRu24pJoRVe/B8emt8p2fTpk2MGDqUwzly4Js/P8ydC+nMhZSR9eLtjdOR5IpBCAezdEb+9opDDx7Plr8E/n2n0adSE6YD6tNPiW/fHm7dsnvf6TUR2XJmmtGz34yytQLq7RWHWLDzbxKNRTSJWrNg59+8veLQg23sqVA6ceIET3brxo/ZAygae49eTV+m4dcHndLu74xKqoySqiQhHCx5UZ/UfJV68AWWTGvN7X2reGrDV8zUSSSULUvAunVQvrzN+7W2gshSfIH5Atg+qoXN+3WUjFRAPTp6zUOfKRg+69MfdjB5b1urfaKioqhfvz79zl/gnXuxjGsxmG9rh1oVV0Z5SlWSNCUJr+TsPyB7WDojN/cFppQiT63OhBcrz8ll41l85gzqiSfIvnw5qm1bm/ZrbRORpzYP2dLElczcZ2rucVubdu7cuUOnTp2ofOYMYxISWV65Gd+GdLE6rozKSBOUM0hTkvA6nlTWZ46lNmHfNNqlswdW5MRzs2hV9glOxMaS1K4dsePH2zSNhrVNRK5uHrJWRjpfLX2maX3W6YmPj6dHjx5E79rFEn9/jhQtx+h2wx7qV7gYHes1A/5sJYlBeB1PKuszx1Jbce+6pdJ8nW9AHm72GM+mCRNYoRQB777L9VatrB4pbUsFUXKN/5lJHdk+qoXbkwJkrALK0mea3mdtSVJSEgMHDiTi11/ZVbAg/jlz8m7/8dzze3i8AnjeSYmjSGIQXseTyvrMsXRGHlKmQHrFLPgoH8p0GEDJHTv4KH9+8m/cyD+PPkrisWPp7teTOi8zIiPxTwirRt96pR9cIfgqRd96pU2qkqyVlJTESy+9xNIFC9hfqhR5YmLgl194pleTNJdq9aSTEkeRPgbhdTyprM+S1G3Fyc1f6bUMaQxXRNtHtaDS2bNM7tqVwRs3cq9qVWI+/ZSiL76Y5j4xvt4T+17Sk9H4J4RVy1AiSCkpKYlhw4YxZ/ZsIitUoMTJk7B0KdStS5hxmynrjpv9vQPPOSlxFKlKEl7HG+fvsVQJZE7qOXWWzZhBqREjqJ2YyJ+tWlF51SpU9uxOijTrSU4Ks2bNYkdwMPX374ePP4ZXX31oW0+t6LKWzJUkMi1P7TxNiy1nlKmvfLoNH06Ro0f5OTCQKuvXc7JYMa7u2uXoELOkhIQEBg8ezKxZs9hUt64hKYwcaXFks7c311lLmpKEV/KUsj5rWWr+UmAyE6elL5ky5ctT6u+/+eXZZ2k8fz4+9eqx7aWXaDRzJsqOCpysLDY2lt69exMeHs765s1ptmkTvPACfPSRxZHN3t5cZy1pShLCBSw1f3WvFcimY9ds+pI5vWEDsWFhVI2J4bdSpaj822+UrFjR2YeQqURHRxMaGsq2bdvYFhpKwxUroF8/+O478Mm8DSkuGeCmlOoBjAUqAXW01ma/rZVS3wKdgKta66opHi8ALAaCgLPAU1rrm/bEJIQncuSZ5qMtW5J47Rq7O3Wi1YYNnKlcmSWvvEL3qVPx9bVcPSMMTp06RefOnTl96hSR3bpR7eefoWdP+PbbTJ0UbKK1zvAPhoTwOLAZCEljuyZATeBwqscnA6OMt0cBH1mz31q1amkhhNaXFi3SVx55RCeA/rp4cb1n+3Z3h+TR3p+9WGcLyK19sufS31RvpjVoPWCA1gkJ7g7NJYAIbcV3rF3pUWt9VGudbgGv1norcMPMU6HAPOPtefCgMkwIYYXiPXtS5PJl/m7cmEGXL+PfsCETu3cnKirK3aF5FK01g0a+z3sv9sEvIB/zywYz8OBmfqzViRXDxoFcaZlw93VTUa31ZQDjv0UsbaiUGqKUilBKRVy7ds1lAQrh6VS+fJTdupW7ixYRlCMHby5bxnclS/L5tGnEx8e7Ozy72bvm9M2bN+nWrRvfTh1LkaAa/JqrAE8f28bM+j15q+VQpvx+0kmRe690E4NSar1S6rCZn1BXBJhMaz1Hax2itQ4pXLiwK3cthFfI0bMnec6f53bnzrx+7x6tR4ygT6kyVHpmAkFvrvLKeX3snRdrx44dBAcHs2rVKqo26MXm21E0uHCYke3/j2lN+oFSmW5wmiOkmxi01q201lXN/IQ7YP//KKWKAxj/veqA9xQi6ypQgPwrV6LXraNQwcL89M9lRn3/DknzX+P0wd0OndfH3jN5a1iaF+v1JZFp7i82NpYRI0bQqFEjlFIcnDaNbQfWEHjrKgN6vM9P1ds82FaDVyZNZ3J3U9JKoL/xdn/AEclGiCxPtWlD91fm8Vnd7vRWPhy6coreC0dz4fvRvDNnmd3v76oZbtOawtzS/rZu3UqNGjWYNm0azw8ZwrG+fak0fDiqRHF6DvqE7UE1HnpNZp0ML6PsSgxKqa5KqQtAfWC1Umqd8fESSqk1KbZbCPwBPK6UuqCUGmR8ahLQWil1EmhtvC+EcICzd5KY2uxZ2g36gn1BNZgO7Pr7IEGz/4+WLVqwadOm5OpAm7lqhtu05r9Kvb/z58/Tu3dvmjZtyv3799m2cCFf/PUX2SdMgN69yXtgL0MHtyfQwntmxsnwMkoGuAmRSZnM66M1rU7t5q2NX1Mu+jI7/PwYHh9PUq1avPbaa/To0QM/Pz+r37vsqNWY++ZIPc+TvcwNDEy9v/2jGjJ9+nSmTp1KUlISb77xBqNLlyb7669DfLxh3qMhQ0xGM7sqfk8jcyUJkcWZzOujFOvL1yX0hTlEjppI/Xz52A1MPn6cqU8/TVBQEO+99x7nz5+36r0zsnZCRiTPi2Vu4Z2k+3dI3LuUsmXLMm7cODp27MjJ339n7IEDZH/uOahSBSIjYejQh6a4cFX83koSgxCZlLnJBsf3COaJD99CnToF48fTPFs29gFLExL4fdw4goKC6NSpE0uWLCE21nK1jisnkwsLDmTaU0882F/8zUvcWP8lF78YwPn139G0aVMO7NzJkmrVKNm6NWzYAFOmwNat8Nhjbo/fG0lTkhBZ2a1b8MknMGMG3LzJuZIlGXvnDvNv3iRn7tx07dqVrl270qZNG3LkyGHyUleuux0TE8M7n8zlu3nfE30yAuXrS5O2Xfj43TeoeegQjBsH58/DU0/BtGlQsmS67+nJ64Y7i7VNSZIYhBCG5UO//dbQHn/uHPeKFGFNyZKMPnWKE//+yyOPPEKrVq1o3bo1LVq0oEqVKhme1dXaL+SrV6+ydu1aVq9ezapVq7h79y5lypShf//+PD9gAMU3bYIPPoDTp6FuXcPtFp6/JoI7SWIQIpNxyRluQgKEh8OsWbBhA9rXl6iQEFblz8+UY8c4cvYsAIULF6Zu3brUrl2bkJAQKleuTKlSpdKdxM/SLLPjOlWgcq5Y9u7dy86dO/njjz84cOAAAEWLFiU0NJS+ffvSMCgIn6++gtmz4do1qFEDxo+Hjh0tTpUt/iOJQYhMxC2r1p04Ad98Az/8ABcvQs6c3GncmL2BgSyOiWHzoUMcPXr0Qclr9uzZeeyxxyhZsiTFihWjWLFi5M6dm1M34th6+hbRd+5BUiJJCXEk3Ysh8c5NEu9EkxB9hYRb/4BOAiB37tzUrVuXJk2a0LFjR2oEBeETHm6IY+NGQ2ydOhkW02neXBKCDSQxCJGJuHVJyaQkQ0fuokWwciVcvmz4Mg4OJq5BA04GBnIgWzYO/PMPJ0+e5NKlS1y5coUrV65YnKtJZfPHJ0c+fHPkJVveovgVCOSzFzsRHBxMxaAgfA8fNnQir10LO3ZAYiKUKwdPPw3PPGOxU1mkzSXrMQghXMPSCGCXzPPj48OKvOWZEvQUl/t1pHnMeUYmnKLi8X34z5lDlbg4qgBPFy0KwcFQvz6ULQtlyzJ4zRlO3oNoXz8SsvmDrx/4+JJDJ5H3/h3y3o+hxL/XqXI/ir4bNxr6OA4dMjRpAdSsCW+8AV26GPoRLFwdZMWOZGeSxCCEF7C0NKgr6u5NmrGUDxtyl2GHXzk+fP51wioWgIgI2LfP8BMZCf/7n6EzG/jKlh0VKwZVqxrWXK5dGxo0gKJFbYuP/6a3ACQ5ZJAkBiG8wMi2j5sdAXw3LoEV+y869QswrekvwoJbQOPGhp9kWsP163D2LK9/sZ74qBvkvReDX2IiCg1acy97ALf8c+JXsACdOtamZfu6kKoc1jHxSWLICEkMQniB5C+4sSv/JDr2v3b7m3fjnX52bHMzllJQuDAULkzj/yvh9E5zRzezSbOUjHwWwmuEBQeSM/vD53LOnvzNnukjzI2+dnQllSOnt3DVrLGeTq4YhPAi7uiENteMZcv0EWHBgU4947Y3vpSkWcpAEoMQXsQdndDJX4ie2rziyPjcWv3lQSQxCOFFHHl2bAtnn/Xby1HxubP6y5NIH4MQXsQVbfaZia3Lj8qsqwZyxSCEl/H0s3dPkZHxDZ7ebOYqkhiEEB7LntLRjHYkS+KVxCCEsIMza/7tHdEsHckZJ30MQogMcXbNf1pn/NaQ5TszThKDECJD7P3iTo+9Z/zmOpIVhgRmTUd0ViaJQQiRIc5uqrH3jD9lBRcYkkLyIgNZdUSztSQxCCEyxNlNNY4oHQ0LDmT7qBYE5gsg9cozzp5KxJtJYhBCZIiza/4dOWZDOqJtI1VJQogMcUXNv4xodg+7EoNSqgcwFqgE1NFam11vUyn1LdAJuKq1rpri8bHAYOCa8aG3tNZr7IlJCOE61n5xu3sqa3dNJeKt7L1iOAx0A75MZ7vvgM+A+Waem661nmpnHEIID5IyEeQN8ONOXALxiYZWfmvGIzg6kciIZtvYlRi01kcBlIV1WFNst1UpFWTPvoQQ3iH1wLSUCwslS2sEsrOW6pQRzdbzhM7nYUqpg0qpb5VS+d0djBDCPubGN5hjqePX2eMjRPrSTQxKqfVKqcNmfkIdsP9ZwKNADeAyMC2NOIYopSKUUhHXrl2ztJkQws2srfSx1PErFUTul25Tkta6lbN2rrX+J/m2UuorYFUa284B5gCEhISkLkkWQngISxVAKaXV8SsVRO7n1qYkpVTxFHe7YujMFkJ4MXPjG/x8FPlz+Fk1HkHWRHA/e8tVuwIzgcLAaqXUAa11W6VUCeBrrXUH43YLgWZAIaXUBeA9rfU3wGSlVA0MI9XPAkPtiUcI4V7J1USx8Yn4KkWi1gTaWAEUFhxIxLkbLNx1nkSt8VWK7rWk49iV7K1KWg4sN/P4JaBDivu9Lby+nz37F0J4jtTVRIlaPzjTt+VLfcX+i/y89yKJWj94n5/3XiSkTAFJDi7iCVVJQohMwFHVRFKV5H6SGIQQDuGoaiKpSnI/SQxCCIdw1GyrssCO+0liEEI4hKOqiaQqyf1kdlUhhEM4aj4imdfI/ZTW3jdWLCQkREdEmJ3IVQghhAVKqb1a65D0tpOmJCGEECYkMQghhDAhiUEIIYQJSQxCCCFMSGIQQghhQhKDEEIIE5IYhBBCmJDEIIQQwoQkBiGEECYkMQghhDAhiUEIIYQJSQxCCCFMSGIQQghhQhKDEEIIE5IYhBBCmPDK9RiUUteAcxl8eSHgugPD8QZZ7ZjleDO/rHbMjjreMlrrwult5JWJwR5KqQhrFqrITLLaMcvxZn5Z7ZhdfbzSlCSEEMKEJAYhhBAmsmJimOPuANwgqx2zHG/ml9WO2aXHm+X6GIQQQqQtK14xCCGESEOmTQxKqXZKqeNKqVNKqVFmns+ulFpsfH6XUirI9VE6jhXH+5pS6ohS6qBSaoNSqow74nSk9I45xXZPKqW0Usqrq1isOV6l1FPG/+c/lVI/ujpGR7Lid7q0UmqTUmq/8fe6gzvidBSl1LdKqatKqcMWnldKqU+Nn8dBpVRNpwWjtc50P4AvcBooB/gDkUDlVNu8CMw23u4FLHZ33E4+3uZADuPtF7z5eK09ZuN2uYGtwE4gxN1xO/n/uDywH8hvvF/E3XE7+XjnAC8Yb1cGzro7bjuPuQlQEzhs4fkOwK+AAuoBu5wVS2a9YqgDnNJa/6W1jgMWAaGptgkF5hlvLwVaKqWUC2N0pHSPV2u9SWt913h3J1DSxTE6mjX/xwDjgcnAPVcG5wTWHO9g4HOt9U0ArfVVF8foSNYcrwbyGG/nBS65MD6H01pvBW6ksUkoMF8b7ATyKaWKOyOWzJoYAoHzKe5fMD5mdhutdQJwCyjokugcz5rjTWkQhjMPb5buMSulgoFSWutVrgzMSaz5P64AVFBKbVdK7VRKtXNZdI5nzfGOBfoqpS4Aa4CXXROa29j6d55h2Zzxph7A3Jl/6vIra7bxFlYfi1KqLxACNHVqRM6X5jErpXyA6cAAVwXkZNb8H2fD0JzUDMMV4TalVFWtdbSTY3MGa463N/Cd1nqaUqo+8L3xeJOcH55buOw7K7NeMVwASqW4X5KHLzMfbKOUyobhUjStyzhPZs3xopRqBYwBumit77soNmdJ75hzA1WBzUqpsxjaZFd6cQe0tb/T4VrreK31GeA4hkThjaw53kHAEgCt9R/AIxjmFMqsrPo7d4TMmhj2AOWVUmWVUv4YOpdXptpmJdDfePtJYKM29vB4oXSP19is8iWGpODNbc/J0jxmrfUtrXUhrXWQ1joIQ79KF611hHvCtZs1v9MrMBQZoJQqhKFp6S+XRuk41hzv30BLAKVUJQyJ4ZpLo3StlcAzxuqkesAtrfVlZ+woUzYlaa0TlFLDgHUYqhu+1Vr/qZQaB0RorVcC32C49DyF4Uqhl/sito+VxzsFyAX8ZOxj/1tr3cVtQdvJymPONKw83nVAG6XUESARGKm1jnJf1Bln5fG+DnyllHoVQ5PKAC8+uUMptRBDM2AhY7/Je4AfgNZ6NoZ+lA7AKeAu8KzTYvHiz1EIIYQTZNamJCGEEBkkiUEIIYQJSQxCCCFMSGIQQghhQhKDEEIIE5IYhBBCmJDEIIQQwoQkBiGEECb+H1SEjM+MmvXPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7da8842d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Prepare the data\n",
    "N = 100\n",
    "sigma = 0.03\n",
    "x = numpy.sort(numpy.random.uniform(0.0, 1.0, N))\n",
    "def f(t):\n",
    "    return 1.5*t**3 - 2*t**2 + 0.5*t - 1\n",
    "y = f(x) + numpy.random.normal(0.0, sigma, N)\n",
    "X = [[1,t,t**2,t**3] for t in x]\n",
    "\n",
    "# Solve the model\n",
    "M = lse(X, y)\n",
    "M.solve()\n",
    "w = M.getVariable(\"w\").level()\n",
    "\n",
    "# Plot the data points, original function (black) and regression function (red)\n",
    "plt.scatter(x,y)\n",
    "ticks = numpy.linspace(0, 1, 100)\n",
    "plt.plot(ticks, f(ticks), 'k', color='black')\n",
    "plt.plot(ticks, sum([w[i]*ticks**i for i in range(4)]), 'k', color='red')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the unconstrained model leaves open numerous possibilities. It wold be just as easy to phrase the problem without intercepts, add additional linear or conic constraints on $\\mathbf{w}$, compute various statistics of the solution and so on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regularisation\n",
    "\n",
    "Regularisation is a technique which helps stabilize computational results, which tend to be sensitive to small perturbations when the matrix $\\mathbf{X}$ has nearly dependent columns. In machine learning regularisation is a major technique that  helps avoid overfitting. A general regularized least squares regression problem is usually written in the form\n",
    "\n",
    "$$\\mathrm{min}_{\\mathbf{w}\\in\\mathbb{R}^d} \\|\\mathbf{y}-\\mathbf{X}\\mathbf{w}\\|_2^2 + \\lambda\\phi'(T\\mathbf{w})$$\n",
    "\n",
    "where $T$ is a linear operator, for example $T\\mathbf{w}=\\mathbf{w}$ or $T\\mathbf{w}=\\mathbf{w}-\\mathbf{w}_0$, $\\phi'$ is an additional penalty function and $\\lambda$ controls the tradeoff between regularisation and fit. The two main regularisation terms occurring in practice are:\n",
    "\n",
    "* quadratic, also known as *ridge regression*, leading to the problem $$\\mathrm{min}_{\\mathbf{w}\\in\\mathbb{R}^d} \\|\\mathbf{y}-\\mathbf{X}\\mathbf{w}\\|_2^2 + \\lambda\\|T\\mathbf{w}\\|_2^2$$\n",
    "* linear, also known as *lasso (least absolute shrinkage and selection operator)*, leading to the problem $$\\mathrm{min}_{\\mathbf{w}\\in\\mathbb{R}^d} \\|\\mathbf{y}-\\mathbf{X}\\mathbf{w}\\|_2^2 + \\lambda\\|T\\mathbf{w}\\|_1.$$ When $T=\\mathrm{id}$ lasso regression tends to give preference to sparser solutions.\n",
    "* *Elastic net*, a combination of quadratic and linear regularisation terms:\n",
    "$$\\mathrm{min}_{\\mathbf{w}\\in\\mathbb{R}^d} \\|\\mathbf{y}-\\mathbf{X}\\mathbf{w}\\|_2^2 + \\lambda_1\\|T_1\\mathbf{w}\\|_1 + \\lambda_2\\|T_2\\mathbf{w}\\|_2^2.$$\n",
    "\n",
    "Below is the most general version formulated as a conic model, for simplicity with both $T_1$ and $T_2$ being the identity.\n",
    "\n",
    "$$\n",
    "\\begin{array}{rl}\n",
    "\\mbox{minimize}      &   t + \\lambda_1 p_{\\mathrm{lasso}} + \\lambda_2 p_{\\mathrm{ridge}} \\\\\n",
    "\\mbox{subject to: residual} &  (0.5, t, \\mathbf{y}-\\mathbf{X}\\mathbf{w}) \\in \\mathcal{Q}_{\\mathrm{rot}} \\\\\n",
    "                       &  \\mathbf{w}\\in\\mathbb{R}^{d} \\\\\n",
    "\\mbox{lasso}    &  p_{\\mathrm{lasso}} \\geq \\mathbf{p}_1+\\cdots+\\mathbf{p}_d \\\\\n",
    "                       &  -\\mathbf{p} \\leq \\mathbf{w} \\leq \\mathbf{p} \\\\\n",
    "                       &  \\mathbf{p} \\in \\mathbb{R}^d \\\\\n",
    "\\mbox{ridge}    &  (0.5, p_{\\mathrm{ridge}}, \\mathbf{w}) \\in \\mathcal{Q}_{\\mathrm{rot}}\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "Again, we have a straightforward translation to a *Fusion* model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Implement the lasso part of the constraints, and return the bounding variable\n",
    "def lassoVar(M, w, d):\n",
    "    p = M.variable(\"p\", d)\n",
    "    plasso = M.variable(\"plasso\")\n",
    "    \n",
    "    # plasso >= sum(p_i)\n",
    "    M.constraint(Expr.sub(plasso, Expr.sum(p)), Domain.greaterThan(0.0))\n",
    "    \n",
    "    # p_i = |w_i|\n",
    "    M.constraint(Expr.add(p,w), Domain.greaterThan(0.0))\n",
    "    M.constraint(Expr.sub(p,w), Domain.greaterThan(0.0))\n",
    "    \n",
    "    return plasso\n",
    "\n",
    "# Implement the ridge part of the constraints, and return the bounding variable\n",
    "def ridgeVar(M, w):\n",
    "    pridge = M.variable(\"pridge\")\n",
    "    M.constraint(Expr.vstack(0.5, pridge, w), Domain.inRotatedQCone())\n",
    "    return pridge\n",
    "    \n",
    "# Regularized least-squares regression\n",
    "def lseReg(X, y, lambda1 = 0, lambda2 = 0):\n",
    "    n, d = len(X), len(X[0])\n",
    "    M = Model(\"LSE-REG\")\n",
    "    \n",
    "    # The regression coefficients\n",
    "    w = M.variable(\"w\", d)\n",
    "    \n",
    "    # The bound on the norm of the residual\n",
    "    t = M.variable(\"t\")\n",
    "    r = Expr.sub(y, Expr.mul(X,w))\n",
    "    # t \\geq |r|_2^2\n",
    "    M.constraint(Expr.vstack(0.5, t, r), Domain.inRotatedQCone())\n",
    "    \n",
    "    # The objective, add regularization terms as required\n",
    "    objExpr = t.asExpr()\n",
    "    if lambda1 != 0: objExpr = Expr.add(objExpr, Expr.mul(lambda1, lassoVar(M,w,d)))\n",
    "    if lambda2 != 0: objExpr = Expr.add(objExpr, Expr.mul(lambda2, ridgeVar(M,w)))\n",
    "    M.objective(ObjectiveSense.Minimize, objExpr)\n",
    "    return M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regularized least squares example\n",
    "We consider the same dataset we had previously, but this time we try to fit a polynomial of very high degree. This will lead to overfitting, which can then be reduced by varying the regularization parameters $\\lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sparsity report:\n",
      "lambda1  nonzeros in w\n",
      "\n",
      "0.1      4 \n",
      "0.01     4 \n",
      "0.001    6 \n",
      "0.0001   10\n",
      "1e-05    8 \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD8CAYAAABpcuN4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd4VFX6wPHvmZn03ishFQg1QIBQlSYIiugKYkFcVBTbqqs/dXd11V0V19W1ocKKoq6IigpYQJrSIbRQE0gjnZRJ71PO748ZQgJJCCQQyvk8z32S3Dn3zpkR553T3iOklCiKoihKSzSdXQFFURTl0qYChaIoitIqFSgURVGUVqlAoSiKorRKBQpFURSlVSpQKIqiKK1qd6AQQngKIdYKIZKtPz1aKPcvIcRhIUSiEOJdIYSwnh8ohDgohEhpfF5RFEW5NHREi+JZYL2UMgpYb/27CSHEMGA40BfoDQwCrrE+/CEwB4iyHhM7oE6KoihKB+mIQHET8Jn198+Aqc2UkYA9YAvYATZAvhAiAHCVUm6XlpV/n7dwvaIoitJJdB1wDz8pZR6AlDJPCOF7egEp5XYhxG9AHiCA96WUiUKIWCC7UdFsIKi5JxFCzMHS8sDJyWlgjx49OqDqiqIoV489e/YUSSl9zvW6NgUKIcQ6wL+Zh/7axusjgWgg2HpqrRBiFFDTTPFmc4pIKRcCCwFiY2Pl7t272/LUiqIoipUQIuN8rmtToJBSjmvlifOFEAHW1kQAUNBMsZuBHVLKSus1q4A44AtOBQ+sv+e2tfKKoijKhdcRYxQrgVnW32cBK5opkwlcI4TQCSFssAxkJ1q7rCqEEHHW2U53t3C9oiiK0kk6IlDMA8YLIZKB8da/EULECiE+tpZZBqQCB4H9wH4p5Y/Wx+YCHwMp1jKrOqBOiqIoSgcRl2OacTVGoSgKgMFgIDs7m9ra2s6uyiXF3t6e4OBgbGxsmpwXQuyRUsae6/06YtaToihKp8jOzsbFxYXQ0FDUWl0LKSV6vZ7s7GzCwsI65J4qhYeiKJet2tpavLy8VJBoRAiBl5dXh7ayVKBQFOWypoLEmTr6PVGBQlEURWmVChSKoijttHr1arp3705kZCTz5s074/FNmzYxYMAAdDody5Yt64Qato8KFIqiKO1gMpl4+OGHWbVqFUeOHOGrr77iyJEjTcqEhISwePFi7rjjjk6qZfuoWU+KoijtEB8fT2RkJOHh4QDMmDGDFStW0LNnz4YyoaGhAGg0l+d3cxUoFEW5Irz042GO5JZ36D17Brry9xt7tVomJyeHLl26NPwdHBzMzp07O7Qene3yDG+KoiiXiOYWLV9pM7FUi0JRlCvC2b75XyjBwcFkZWU1/J2dnU1gYGCn1OVCUS0KRVGUdhg0aBDJycmkp6dTX1/P0qVLmTJlSmdXq0OpQKEoitIOOp2O999/nwkTJhAdHc306dPp1asXL7zwAitXrgRg165dBAcH8+233/LAAw/Qq1fntH7Ol0oKqCjKZSsxMZHo6OjOrsYlqbn35nyTAqoWhaIoitIqFSgURVGUVrU7UAghPIUQa4UQydafHi2U+5cQ4rAQIlEI8a51RzuEEL8LIY4KIRKsh29766QoiqJ0nI5oUTwLrJdSRgHrrX83IYQYBgwH+gK9gUFYtkM96U4pZYz1aG7PbUVRFKWTdESguAn4zPr7Z8DUZspIwB6wBewAGyC/A55bURRFucA6IlD4SSnzAKw/z+g6klJuB34D8qzHr1LKxEZFPrV2Oz0vrrQljYqiKJe5NgUKIcQ6IcShZo6b2nh9JBANBANBwBghxCjrw3dKKfsAI63HzBbuMUcIsVsIsbuwsLAtT6soinJRnC3NeF1dHbfddhuRkZEMGTKE48ePA6DX6xk9ejTOzs488sgjF7nWbdemQCGlHCel7N3MsQLIF0IEAFh/NjfGcDOwQ0pZKaWsBFYBcdZ751h/VgBLgMEt1GGhlDJWShnr4+Nzrq9TURTlgmhLmvFFixbh4eFBSkoKTzzxBM888wwA9vb2/OMf/+Df//53Z1S9zTqi62klMMv6+yxgRTNlMoFrhBA6IYQNloHsROvf3gDW8zcAhzqgToqiKBdF4zTjtra2DWnGG1uxYgWzZlk+Jm+99VbWr1+PlBInJydGjBiBvb19Z1S9zToiKeA84BshxL1YAsI0ACFELPCglPI+YBkwBjiIZWB7tZTyRyGEE/CrNUhogXXAfzugToqiXG1WPQsnDnbsPf37wPVndiU11pY0443L6HQ63Nzc0Ov1eHt7d2x9L5B2BwoppR4Y28z53cB91t9NwAPNlKkCBra3DoqiKJ2lLWnGL/dU5CrNuKIoV4azfPO/UNqSZvxkmeDgYIxGI2VlZXh6el7sqp43lcJDURSlHdqSZnzKlCl89plludmyZcsYM2aMalEoiqJcLRqnGTeZTMyePbshzXhsbCxTpkzh3nvvZebMmURGRuLp6cnSpUsbrg8NDaW8vJz6+nqWL1/OmjVrmuy3fSlQacYVRblsqTTjLVNpxhVFUZSLRgUKRVEUpVUqUCiKoiitUoFCURRFaZUKFIqiKEqrVKBQFEVRWqUChaIoSjudb5pxgNdee43IyEi6d+/Or7/+2nB+9uzZ+Pr60rt374vxElqlAoWiKEo7tCfN+JEjR1i6dCmHDx9m9erVPPTQQ5hMJgDuueceVq9efdFfT3NUoFAURWmH9qQZX7FiBTNmzMDOzo6wsDAiIyOJj48HYNSoUZdMPiiVwkNRlCvC6/Gvk1Sc1KH37OHZg2cGP9NqmfakGc/JySEuLq7JtTk5OR34CjqGalEoiqK0Q3vSjF8u6cdVi0JRlCvC2b75XyjtSTPelmsvBe1uUQghPIUQa4UQydafHi2Ue10Icch63NbofJgQYqf1+q+FELbtrZOiKMrF0p4041OmTGHp0qXU1dWRnp5OcnIygwcP7oyX0aqO6Hp6FlgvpYwC1lv/bkIIMRkYAMQAQ4CnhRCu1odfB/5jvb4EuLcD6qQoinJRNE4zHh0dzfTp0xvSjK9cuRKAe++9F71eT2RkJG+99VbDFNpevXoxffp0evbsycSJE5k/fz5arRaA22+/naFDh3L06FGCg4NZtGhRp73GdqcZF0IcBa6VUuYJIQKA36WU3U8r8zRgJ6X8p/XvRcCvwLdAIeAvpTQKIYYCL0opJ7T2nCrNuKIooNKMt+ZSSzPuJ6XMA7D+9G2mzH7geiGEoxDCGxgNdAG8gFIppdFaLhsIau5JhBBzhBC7hRC7CwsLO6DaiqIoSlu0aTBbCLEO8G/mob+25Xop5RohxCBgG5YWxHbACDQ3vN9sE0dKuRBYCJYWRVueV1EURWm/NgUKKeW4lh4TQuQLIQIadT0VtHCPV4BXrNcsAZKBIsBdCKGztiqCgdxzfA2KoijKBdQRXU8rgVnW32cBK04vIITQCiG8rL/3BfoCa6RlgOQ34NbWrlcURVE6T0eso5gHfCOEuBfIBKYBCCFigQellPcBNsBm60KScuCuRuMSzwBLhRD/BPYBnTe0f5Es35fDG78eJbe0hkB3B56e0J2p/ZsdmlEURel07Q4UUko9MLaZ87uB+6y/1wI9W7g+Dbj0Jg5fIMv35fDc9wepMVgSf+WU1vDc9wcBVLBQFOWSpFJ4XGRv/Hq0IUicVGMw8cavRzupRoqitNeFSDPe0j3ff/99IiMjEUJQVFR0QV/XSSpQXGS5pTXndF5RlEvbhUgz3to9hw8fzrp16+jatetFe40qUFxkge4O53ReUZRL24VIM97aPfv3709oaOhFfY0qKeBF9vSE7k3GKAAcbLQ8PaF7K1cpinI2J159lbrEjk0zbhfdA/+//KXVMhcqzfjZ7nkxqUBxkZ0csFaznhTlynAh0oybzeaz3vNiUoGiE0ztH6QCg6J0sLN9879QLlSa8Usp/bgao1AURWmHC5FmvC33vJhUi0K5bKiFisqlqHGacZPJxOzZsxvSjMfGxjJlyhTuvfdeZs6cSWRkJJ6enixduhRommZcp9M1STPe3D0B3n33Xf71r39x4sQJ+vbty6RJk/j4448v6Gtsd5rxzqDSjF99Tl+oCJZJAK/d0kcFi6uYSjPeskstzbiiXHBqoaKidB4VKJTLglqoqCidRwUK5bKgFioqSudRgaKDLN+Xw/B5Gwh79meGz9vA8n05nV2lK8rTE7rjYKNtck4tVFSUi0PNeuoAKiPshacWKipK51GBogO0NtCqPsg6jlqoqCido11dT0IITyHEWiFEsvWnRwvlXhdCHLIetzU6v1gIkS6ESLAeMe2pT2dRA62KcvWaPXs2vr6+9O7d+5yv3bNnD3369CEyMpLHHnusIaXHiy++SFBQEDExMcTExPDLL790dLXPSXvHKJ4F1kspo4D11r+bEEJMBgYAMcAQ4GkhhGujIk9LKWOsR0I769Mp1ECroly97rnnHlavXn1e186dO5eFCxeSnJxMcnJyk/s88cQTJCQkkJCQwKRJkzqquuelvYHiJuAz6++fAVObKdMT2CilNEopq4D9wMT2PGlaQRnzfklka0oRtad1+XSGiznQajKZqSypQ59bSV5KKccPFpG6r6DJkXlET356OaX51dRWGZpNPKYoSscYNWoUnp6eTc6lpqYyceJEBg4cyMiRI0lKOjOrbV5eHuXl5QwdOhQhBHfffTfLly+/WNU+J+0do/CTUuYBSCnzhBC+zZTZD/xdCPEW4AiMBhrv6vGKEOIFrC0SKWVdc08khJgDzAEYGKDh8Z0jyd/hwUHhicnRD0evYHyCwvALCkXjGgAuAeAWDDq7dr7Es+vogVazWVJeWENxbhX63EpKTlRTWVxLRXEtVaV1nOvnvq29FhcvB1y97fEMdMKniws+IS64eNl3akZKRelIm785RlFWZYfe07uLMyOndzvn6+bMmcNHH31EVFQUO3fu5KGHHmLDhg1NyuTk5BAcHNzwd+MU42BJ4fH5558TGxvLm2++iYdHsz37F8VZA4UQYh3g38xDf23LE0gp1wghBgHbgEJgO2C0PvwccAKwBRYCzwAvt3CfhdYyxEaHSE3cTOxOZOCjz0ZXnYxn1g4cs5uJMS4B4B5iPbo2+j0E3LqAzrYtL6NZ5upqDPn5GPMLGF2Qz0ifIkwOFZgrqzB/s47sTyqRRiMNn+xSglaLxtGxyVHn5E0JXuhrHCksFhQXGjAarGmGBbh42uPqZU9Qdw9cPO1xcrfDzlGHnYMOW0cdOhuNpaDlSTDUmqirNlJXY6S6vJ4KfS0V+hpK86s5flCPNFvqY+9sQ1A3D4J7eBDc3QM3XwcVOBSlnSorK9m2bRvTpk1rOFdXd+ZnU2vpyefOncvzzz+PEILnn3+eP//5z3zyyScXrtJncdZAIaUc19JjQoh8IUSAtTURABS0cI9XgFes1ywBkq3n86xF6oQQnwJPtanWTr7YXv9Kk+hVUFbDuqTjHDl2lKzjadjX5BMkioiuKaG7LMG/eDsO1d8jZOOuKgGugeAVAV6RTQ/3rqC1vD2mykpqDx2i7lgydWmp1KelU5eWhqm5/WqFQOPsbDmcHBE2tiBAWD/IpdlMTR0UaQMpsg+l1NmfWnsvADSmSlwqMgmoyMRVW4GHuwaPLi44du2CXUQ49tHR6AIC2vVhbqw3oc+tojCzgvy0MrKPlpC61/Kfzd3PkYj+PkQM8MW7i7MKGspl5Xy++V8IZrMZd3d3EhKaDrmaTCYGDhwIWLLJzp07l+zs7IbHG6cS9/Pzazh///33c8MNN1yEmresvV1PK4FZwDzrzxWnFxBCaAF3KaVeCNEX6AussT52MsgILOMbh863Ir5uDkwZEs2UIdFIKUkuqGRLchHfphSxI01Pdb0JG2FidKCRsf51xLqVE6rVoy07DvpUOPQd1JYhJdRXaKkudKCmwpOaIi31+jqwBn+Niwu2ERE4jxqFbUgIOn8/bPz80Pn5ofPxQeN85gesyWQmP62czMN6Mo8UU5hZAYCDiw1BUe74Bdjg7WrAlTJkiS/GE2bqMzKpz8ykduNuqhoFJK2bG3bR0dj36IFD3z44DByITaN/VGejs9XiF+qKX6grvUcFIaWkrKCGrMRi0hIK2bsmkz2rM3DzcSB6eAA9hgbg5Hbhu+8U5Urh6upKWFgY3377LdOmTUNKyYEDB+jXr98ZwcPFxYUdO3YwZMgQPv/8cx599FHAMn4REBAAwA8//HBeM6o6UruyxwohvIBvgBAgE5gmpSwWQsQCD0op7xNC2AN7rZeUW88nWK/fAPhg6TdJsD521k7Gxtlj25J6uqq+js1p6WxOTWd3djbpJflITSW2tjX4uhsJ1NXTN72I0GOFdD1WhluZpWes3BHSAgSpAZASIEj3g2JnAcLSPrBBg06jQ6exQaezQ6fR4aBzwMnGCU+DH95FXXErCML+hCfCoAMhsQuSeHezo0tPT0LD/fFw8EAjWp9TYKqspD4lhdrERGqPJFKblETdsWNIa3PWJigIh4EDcBwwEMchg7ENDT3v1kBNZT3p+4s4uuMEucmlCI0gtI8XfUcHE9TdQ7UylEvKpZA99vbbb+f333+nqKgIPz8/XnrpJcaMGcPcuXPJy8vDYDAwY8YMXnjhhTOu3b17N/fccw81NTVcf/31vPfeewghmDlzJgkJCQghCA0NZcGCBQ2Bo606MnvsZZ1m/NSK6Do0tsUI2yLs7Mu5pqcOZ6dK8qryyKvMo7CmEEnT1+lWKRmcLBicBL0yjejMUGmn4ViYG3k9fKkZEIo2xAetRodWaNDWV6GtLkZbXYymqghjVSGGmmKMmDEiMGKHpj4GbXU/7EqjsKu2dCdV25eR65FMuutBslyTqNfVNqmHVmjxdvAmyDmIAOcAAp0CCXQOJNApkCCXIIKcg9Bpzmz4SYOB2qSj1OzdQ/WevVTv3dvQFWYTFITTiBE4jRiOU1wcWheX83qfS05Ukbgtj6TtedRUGPAJcaH/+BAiBvig0arsL0rnuxQCxaXqqg8UEb0i5IOLH2TJ3t3Ua/IRNiUI0WiPWakjxDWQAOcAApwsh7eDN94447UzGfs12zHvOQBSYtM1BDFqDEfC+7PO7MWWtGL0VfUARPk6MyLKmxGR3gwJ98LZ7tQHtpSS4uwysnankJWoJydHg8mkRSsMBNkcpItdAiF2+/Dw1CAC+2L270eFTzf07oHotQJ9jR59rR59jZ786nxyK3PJq8rjRNUJTI3GUXQaHSEuIYS7hRPmFtZwhLuF42jj2KQ+howMqrZvp3LLVqq3b8dcXQ1aLY4DBuAyfhwu48Zhcx7bKRoNJo7uOEHCuixK86tx9bZn0OQwug3xR6NRLQyl86hA0bKrPlB08Y2SNzx6B3vdDlMmbDHXe2Ou88Fc7400eIDJifR5NwKWD9CahARKv/uOilWrMVdVYdOlC25TpuAy4TrsoqKadKeYzZLEE+VsTSlic3IR8enF1BnN6DSCwYFuxDk5418jqc6soqrU0vXj4e9Il56ehPTyIjDKHRtjGeTtb3oUp556Ac7+EBxrPQZBYH+wdQLAaDZSWF1ITmUOOZU5pJelW47ydLLKszBKS7eYQNDVtSvRntH08OpBD88eRHtG42FvmUInDQZqEhKo3LyFyt9/p+7YMQDse/fGZfx4XMaPxy487Jzed2mWpB8oYvcvxynMrMDD35EhU8IJ7++juqSUTqECRcuu+kAR2TVaPjl5PmYzZOpMJNmYSLExUW3tDQlyd2DLkyMoX72a4s+/oPbQIYSjI64TJ+J+81QcYmPP+sEmzZKywhoyj5Vw5EABhekVaCotH9K1QpJjK7EJdqRHP19G9fcn3Nup9XvWlsOJg3DiAOTug+zdp4KH0IJfT0vQOHl4RoCmafeOwWwguyKbtLI0jpUcI0mfRFJxErlVuQ1l/J386ePdhxifGGJ8Y4j2jMZGa0P98eNUrFtH+dq11O4/AIBddDRuN96I6+TJ2Pg1twSmhfdGStL2FbJzZRolJ6rx7erCyNu64R/u1uZ7KEpHSExMpEePHuqLymmklCQlJV3dgSI2Nlb+tmYLP3x3lMyEIjxMln8kRRozhToD42xzCTywDvvcozh28cdz5l2433QTGienM+5lNJioLK6jXF9DeVEt+pxK9NmVFGVXYqizdAHZOeoIiHQnIMIN52Anko31bEnVsyWlkKxiSz6nQDd7hkd6MyLKm+GR3ng7t2GmUJUecvZA9i7LkbMH6sotj9m7W1ocXYdByDAIGtDi4sHS2lKSSpJI0idxRH+EA0UHyKm0LNyx1djSy7tXQ+CI9Y/FobiaijVrKPvpZ2oPHAAhcBoah+uUKbiMG4/W+cz3qTlms+TojhPsWJFKdVk9PeL8ibs5Qs2SUi6a9PR0XFxc8PLyUsHCSkqJXq+noqKCsLCmvQZXXaA4Oevph73ZLPrxGJ6FVcTWVGGjccOsPfVBpdEJHF1ssXXQIYRAaMBsktTXWBakGWqbpgCxtdfiFeyMd5Az3l1c8At3xdPfCdFCX3ymvprNKYVsTSlia4qeshoDANEBroy0Bo3BoZ442Gqbvb4JsxmKjp0KHFnxUJhoeUxnD0HWwNF1GHQZ3NBd1ZyC6gL2F+5nf8F+EgoTOKI/gsFsQCM0RHtGMyRgCEP8h9Cr2oP6VWspW/kjhuxshL29peU1fToO/WPa9D9ffa2RPasySFiXidZGw+Abwug7posav1AuOIPBQHZ2NrW1tWcvfBWxt7cnODgYGxubJuev2kBhLCmheNEiipd8haypwXnCBOxuv58qex/KCmuoLq+nprye+joT0iyRZonQCOvKZhvsnXW4eDng4mWPi6c9zh525/3NxGSWHM4tY3NyEVuSi9iTUUK9yYytVsPArh6MiPJmZJQ3vQLd0Lb1Q7S6GDK3Q8Y2yNhqGe+QZtDoICDGEjRCR1h+2rU8u6neVM+BwgPEn4hnZ95ODhQdwGg2otPo6OfTj6H+cYwq88d9/V7Kf/4Fc1UVdlFRuE+fjtuUG9G6nb1bqTS/mi3fJpNxSI9vVxfG3B2NV5BzW98+RVEuJCkRGs3VFSjit2yh+PMv0C9ciLm6GtdJk/B+8AHsIiM7u3oNaupNxB8vbhgYT8yzdCu5O9owLMKL4ZHejIz0IcTL8Sx3aqS2HLLjrYFjm6W7ylRvCRxBsRB+reUIjgWtTYu3qTZUs69gHztP7GRn3k6O6C3pt3wdfRntFcfYZHt81+6n/tBhhJ0drhMn4D5jBg4xrbcypJSk7Clg89fHqKsyMuD6rsRODEVro6bTKkqnMZvgv6MRD26+egJF/8hI+U1wF4wnTuA8Zgy+f34Su4iIM8q1ZTHexVRUWcfWFEtrY0tKEXllluZyiKejJWhEeTMswgt3x3PIP2WosXRRpf0O6RstA+XSDLbO0HW4NXBcA749oZUP+KKaIjZnb2Zzzma25W6jylCFjcaGyYZoJhzU4r3pMLKqGvs+ffC8eyauEyYgbM+s58n3vKS4hkkmB8KqwDPQiXF/7IlPl/Nbz6EoSjsdWQHf3I14qfzqCRS97R3kjzfeiO/TT+E0eHCzZU7fnhQsqb9fu6XPJbFLmpSStKKqhqCxI1VPRZ0RIaBPkJu1teHNwFAP7HRtGN84qaYEjm+xBI6030GfYjnv5AsRoyFyPESOBUfPFm9hMBnYW7CXTdmb2Ji9kYzyDOzr4a6MYEZtr8A+R4/WxxuPGTPwmDEDnZdlcWFz73kPs46bjPZQbyZuSgQx47q0ON6jKMoF8vF4qCpAPH7g6gkU/aOi5N6jRxGalrszhs/bQE4zO8wFuTuw9dkxF7J658VoMrM/u4wtyUVsTSlib2YJRrPE3kbDoFBPRkZ5MyLShx7+Luc2SFyWDWkbIe03SN0A1XoQGks3VdR4y+Hf74ypuCdJKUkpTWFtxlrWZqwltSSZvumS2/Y7E5lUDjY2uE2ejOc9sxi7PLfZ9zzcxYGnvXxISygkqLsH4+6JxtnD/nzfKkVRzkXmTvjkOrj+DUTcA1dPoGg8mN2SsGd/prlXJoD0eZMvSL06UmWdkZ1perZYu6qSCywpsLycbC3TcK1Tcc9pFz2zCXITIHmN5cjdB0hLayNynCVoRIwGh5bz3qeVpbEuYx1rM9ZSlpzI9bslYw6Bbb2Z3X5RfBM5loPeEU26uQSQ9tokErflsfmbZLRaweiZPYjo3/a1G4qinKev74L0zfDkEYSdswoUjV1uLYqzyS+vbeim2pJSRGGFZVV4uLdTQ5qRuAgvXO1bHsA+Q2UhpK63BI2U9VBbaln8FxIHPSZD90ng2fLq7YzyDH5O+5n1h1fQa1M21+8y414NR718+CZiItsD+iCFpsl7XppfzdpPj1BwvJy+Y4IZdkskWp0a6FaUC0KfCu8NhJFPwtgXrt7psS251Mco2kNKybH8Smtro5Cd6cVU15vQagT9gt2srQ0f+oe4Y9PW5H0mo2UGVfKvcHQ1FBy2nPftZQkaPSZDQL9mB8SllOwv3M+HW5fguGEtU3bV4l8KOW4OrOh+LeMeeZSbBoc1eioz275P4cCGbPzCXJlwf29cPFVXlKI0p12TclY+BvuXwuMHwMVfBYrmXGqzni6UeqOZfZklbLFOwz2QXYpZgpOtliHhXg3dVFG+p/bKOOt7U5wGSb/A0V8s6zikGVyDocckS9DoOrzZ6bfL9hznjd+X0ztjNbfszyQ8X1LhakPNtOuIfeAvOLmeGkRP2VPAhi8S0WgF4//Yi669vS74e6Uol5N2feEty4F3+sGAu+GGt4CreMHd5ehCB7CyGgPbU/WWqbgpRaQXVQHg52rH8EhvHGy0fLcnm1rjqYy7rf7jqyqCY79C0s+WripjLdi7QbeJ0OtmiBjTbHqRouoiNi5/D5svfyQqtYZyR0HGDf3oc//TRHcZAFi6olYvPIQ+p5LYSaEMviFMzYpSFKvGXehBRg2RBksPQbGnDT/+/Sxd6KuehfiF8Ng+8OgKdGKgEEJMA14EooHBUspmP8GFEBOBdwAt8LGUcp71fBiwFPDEssHRTCllfWvPeTkHis7oEssuqW5Y9LctVU9xVfNvb5vGb+qrIPUf91maAAAgAElEQVQ3S9A4+otlXMPO1TKe0Wtqs0FDSsm+dV9R/NECgg4XUO4Au64NIHT2XK7reRMak5ZNS4+RuC2Prn28GD+7F3YO7d18UVEuf2HP/oxGwoRqG3oZdBitU3S0wNCpEQycGNr8hZWF8HYf6H0LTP2g4XRnBopowAwsAJ5qLlBYt0M9BowHsoFdwO1SyiNCiG+A76WUS4UQHwH7pZQftvacl3Og6OxBdrNZEv6XX1p8/Mnx3RgR5U3fIDd0ZxvfMNZbFvkdXg5JP0JtWaOgcbNlBtVpQaNw11ZS356H254UKu1hw1AnXO+cwc0xd1G0x8jmr5Nx9XFg0tw+ePi3LTmholypRry2gSE5JsKNWrbZGdhpb0QAU40OhFbB2Hui6RHXzM53616ELW/DI7vAO6rhdKd3PQkhfqflQDEUeFFKOcH693PWh+YBhYC/lNJ4ermWXM6B4lKYtttSsLLRCoxmiZTgYq9jaLhXQ2LDsLOlUW82aLhZxjROdk81GtOoOXSIY/95BdutCVTZwepBWqpvGcN1XneR9m0dJqOZ8bN7kWCquyrGmRSlOYvf20vV4VJ+dajngJ2lF8LBRsurU3sjf8unOLeKO16Mw9G1UZaEmhL4Tx/LdPdpnzacrjfVY6ezO69AcbHa90FAVqO/s4EhgBdQKqV1Nx7L+WY/BYQQc4A5ACEhIReuphdYoLtDsx/S57Qeop2entC9xe6vUd182G5Nob45uYg1R/IBS4tneKQXI6J8GB7hhdfpadR1tqcW8Bn/Yw0aP0DST7D/K3D0sgSMPtOhy2Acevem36KvqE1MJPO9t/jDhi1U71rLT4PWkzkqliGpd/HzBwfY4Wgkx8YAAnJKa3ji6wQe/zqBIBU0lCvc8YNFVB0uxamPB/qKEsRpX5ZKgjz46uV49q7JYMStp1oN7PgI6itg5J8bTiWXJPN/m/7vvOvSpkAhhFgH+Dfz0F+llCvacotmzslWzp95UsqFwEKwtCja8JyXpJY+pJ+e0P2i1eHkh2tL39Qn9w1gct8ApJRkFlez2bpafPWhE3yzOxuAno3TqId5Ym/TKM1Ik6DxtmVF+IGvYd//YNfH4N4V+k6HPtOxj46m2wf/pfboMfLfe4fp6zZQvWc3PwzZT7XHLOJK+uFhY2aVowmDOPWPI6e0hue+P9jk9SjKlaKu2sBv/0vCK8iJaXP6cU8zSTU9/J3oNsiPw5tyiJ0Yir2zjWWPm+3zIXoK+PcGYNmxZcyLn4eTzfl35aqup05wuU7bNZklB3PKrAPjhezJKMFgktjqNAwK9WjIhtsr0LX5NCN1FZD4kyVopG+0TLkN6GdpZfS5FVz8qTl8mMJ336Nq40ZKHbRsiBuNq7yJIrtyfnA0UC6bpi2/XBdQKkprtixLZv/6LKY/NwifkJaTaRZlV/D1P3cx8rYo+o7uAmv+Btveh4d2YPCO4PX41/n66NcMCxzGKyNewcfR55Ieo9BhGcweC+RgGcy+Q0p5WAjxLfBdo8HsA1LKD06/R2OXe6C4UlTXG4lPL25YMZ50ogKwpFEfHuHdsGK8i2czadQrTsCh7+HgN5ZUIkIDYaOg3+0QPYXqw0dZ9eTL9MxLIj2wL6nd7qbWxsBPQZvJqu6Fud4PuHxSsihKW5XmV/PVyzvpEefP6Jln3w/861fiEUIw/ZFgeDcGet1MyfWv8eTvT7I7fzd/7P1H/tT/T2g12k6d9XQz8B7gA5QCCVLKCUKIQCzTYCdZy00C3sYys+sTKeUr1vPhnJoeuw+4S0pZ19pzqkBxaSqsqGNbalHDxk0nyi1p1Lt6OTZkwx0W4Y2b42kL9YqS4cA3lpZGaQbYukDvm9nkNIF3vk3ntoOrCKupZV+/udTZufBbxJccsTdSXzSGAIcI1aJQriirPjpIZmIxd70c16Zthfevz2LLt8ncOfp33I/OJ2f2z8yJ/wcnqk7w4rAXuTHixoaynT7r6WJSgeLSJ6UktbCKLcmFbEnRsyNNT2WdEY01jfrJvcUHdm2URt1shsxtkLDEMnvKUEWFcxj/qxnGrnR/Jh/bQ03wZMrcIymw+ZUf+v9MN/c4/j7yT/Tx6dO5L1hROkBhZgXfvLqLwTeGMWhyy3nWGisvquGLv21nuOtnOA2z5cH6NGpNtcwfO58Y35gmZVWgUC5pBpOZ/VmlDdlw92WVYjJLHGy0DA7zbEgz0sPfxTINt67CstnKvi8twUNoyPMexvKDITjqu1Hp1Q/3sgPEh3/Nb1EVDA0azpy+cxjoN7CzX6qinLdVCw6SnVTC3a8OO6dFp1899T0mYy4fDv4SO50DC8YvIMoj6oxyKlAol5WKWgM704obsuGmWNOoeztb0qif3PEvwM3BkgEzYYllmm15DmZbd3aWz2VvySBcKjKIKPuOr4bo2RhSSaz/IB7o9wBD/Iec997nitIZ9DmVLP1HPLGTQhkyJbxN18j6esjfzw/vf0FO+Y2sG/0B869/j2CX4GbLq0ChXNbyymrYmqJv6KoqqrQMU0X4ODVkw40LdcMldyskfAmJP5FW2Y+1ZU+iNdTSd/8HiEDJJ0Nr2BJYzkC/gTza/1HVwlAuG2sWHeb4gSLufmWYZaprM6SUVO/YQekPP1C9fQfGwkKkgKTQ7uR1fYxR9wTQJ67lAXAVKJQrhpSSo/kVDbOpdqYVU2OwpFGP6eLOiEhvrgnR0U+/muJta/gl7XZqTK70Tl+CV+ZuKnp3ZcGQCuJ9yxkWOIxHYh5RYxjKJa2ssJovX9hBzPgQht0S2WyZupQUTrz0MtW7dqFxc8N51CgqbYv4KXs73XMdyQybR1jpDq79y404xMQ0ew8VKJQrVp3RxN6MUsv6jZQiDjZKox4X5skU1xxqdxgoqvCmb+3X+B7YgrnaTGlMKO8PLuWAVyXXdrmWR2IeobvnxVvYqChttfmbYxzamMPdrwzDyf3MpJol//uSgjfeQOPkhPfDD+M+fRpZFceZtfIPaBF8dstPbP13GqbMDAbsf5vg+e/jPHz4Gc+jAoVy1SirNrA9rahhYPy4vhqdhFsNOrpU2xDmtJvY7A8oTbTHVC8oHBDMe4PLSfKo5rqu1/FwzMOEu7etD1hRLrT6GiOLn9tKWF9vxs/u1eQxaTBw4uV/UPrttzhfcw0Br/wTnbc3eZV5zFpxM7V15Swe9DfC+9zB1u9SOLAhi3HZ72POOk7o10uxi2o6oK0ChXLVyiqublgtXp1QTGyFljJtFQM8vqR/0jZKjzpgNmjIG+DD23HVZLgbmBw2mYf7P0yQ86W/Il65sp1cBzHtuVh8u7o2nDfX1ZH96KNUbdqM1wMP4POnxxAaDUU1Rdzz810UV2SxyDaC6DtXApC6r4DVCw4x9f4wqh6/G62TM2HLf0Bjf2r3yPMNFGqzYuWy18XTkRmDQ5h/50AWvT6OnreG44YziWUPMCvibZZNGI0mGgL35/P6h5W8usqePftXceMPN/KvXf+ipLaks1+CcpUymyUHfssiIMKtaZCoryf7sceo2rQZ/xdfxPeJxxEaDWV1ZcxZO4eCqjw+KCwjetK7DdecTPVRUmlD0L/+Rf3x4xTNn98h9VSBQrmiaDSC0eNCue25QXi52jG1yoO+I//M9/d/x9czZlMe6UzE/jLe/qCGZ34y8lP8F0z87nr+e+C/1BjPzOqrKBfS8QNFlBfV0ndMl4Zz0mAg5/EnqNq4Cf+XX8Jjxm0AVBmqeHDtgxwvTePdvDxi4h4Hj9CG61w87bFz0lGYWYHT0KG4/eEW9J98Sl16ervrqQKFckXyDnZm2nOx+Ie7kbkqi1E1tjz/1/8j8qut5L7xHlXRwfQ9UM+CDww89HMZn259h1FfjueFDR9TUFHd2dVXrhIHNmTh7GFHeIw3YBm4znvpJSo3bMDvhefxmD4dgBpjDY+sf4TE4kTerDAT59wVhj3W5F5CCHxDXCjItORc833iCYStLUXvvQdAaXWrG4e2Su03qVyxHJxtufFPMWz9JpmEtZkU51Zy3b29GDN5PEweT11aCtn/+AuDtx9k0EEDmwcUsWjY2/yY8hkepumM6zqWUVG+xIZ6NE2jrigdoDCrgpxjpQy9JQKNdTdJ/X8/pmzZd3g/NBfPO+4AwGAy8OTvT7Infw/zPOMYnfY1zPrRks7/NCV2grLESiKf+Rk/Dwden3AzcvkSdHffy02/nDjvuqoWhXJF02o1jLq9O9fe2Z3sxBKWvb6HkhNVANiFRxLx6TeE//QjrsOGMWqn4NMPjMyML6Rct4BfMucye8mX9HtpDXd9vJOPNqZyKKcMs/nymwCiXDqW78th+LwNPP/GdowCjrtaMgiUr1pF4Vtv4XrDDXg/+igARrORZzY/w5acLbzQczaT9n0HfWdYMi03c99lqQVoAU+TsGz0Ze6J2daOda+8S3E7WhRq1pNy1chNLmX1woOYDGbGzupJeH+fJo/XHjlC4dtvUblpK0YHWB4HP8RqidYEUFF6JwfyLeU9HG0YZs2GO7ylNOqK0ozl+3J47vuDiDoTD5Tbc9DWxFY3M2/F2BH20uPY9+lDyCeL0NjZYZZmnt/6PCtTV/L0wCe5e/MiqCqAh3aAo+cZ9x4+bwO1+lpmV9jzk2M9ibaWzdGePPAdo47vIv2DpdwyuvclvRWqonS6wCh3pj03iNULDrJqwUFixocQNzUcrbXZb9+zJ10Wfkz1vn0UvvMut/62gxt3S5YMy2FD3zeY49eF3v6PsbHIjy0phfx8IA+AUC/Hhr03hkZ44+bQfPoFRXnj16PUGEzE1evQIdhrZ0RbVYXzq/9A6+1F8PvvobGzQ0rJvPh5rExdyUMxD3F3UT7kH4QZXzUbJAByS2sQGjAj8TKdynO2LHQ449O2Mypl+3nXWwUK5ari4mnPLU8NZOsyy7hFfnoZE+7r3WQ1rGP//nRd/ClVO3ZS+M47/PHXfUzbqeGz4Zmsiv4zD9j48drNT5HpOZLNKZaNm37Ym8P/dmRa0qgHuzPSmg23f4j7qTTqylUvt7QGjYSYOh3HdSZKNCZeil+CW3UZwZ8uQefhAcC7+97lq6SvmNVzFg/6DocVoy07QfaY1OK9A90dyCmtoUQj8TKfGlU44RGIfUwMpcuXn3e92zVGIYSYJoQ4LIQwCyFabM4IISYKIY4KIVKEEM82Or9YCJEuhEiwHs0nKFGUDqS1sYxbjJ/dk8KsSr5+JZ6spOIzyjnFDaHrki/psnAB3oHdefhnyWuLdWxKLOQPW54i47truEf7K4tu70HC36/j2weH8uiYKHQawYcbU5mxcAcxL63lnk/j+XhzGkknyrkcu3qVjhPo7kCUQYuLtLQmbju2gUH5SXw95FYc+vYF4OODH/PxwY+Z1m0af+73MGL5XHDwhOtfb/XeT0/ojoONFr32VItCAH+bHI371KnUp6Sed73bNUYhhIgGzMACWt4GVYtlG9TxQDaWbVBvl1IeEUIsBn6SUi47l+dVYxRKRynOq2L1goOU5lcz6IYwBl4f2ux+31JKKtato/Cdd6lPSSHXX8cXI8yIoHqerjLQs/cdMGROw7z28pNp1JML2ZJSRGqhZQDd29mOEZFejIjyYUSkN/5u9mc8l3LlWr4vh12LEnE0w+7ag/xz23/ZHDIA/3nzmDogmCWJS3gt/jUmh0/m1RGvoln1LMQvgDu+hW7Xten+q79Kolc5vO1ey/9d34MHronAVFZG8shRRB880HkpPM6yX/ZQ4EUp5QTr388BSClfU4FCuRTU1xrZuOQox+LzCYh0Y9w9PXH1dmi2rDSZKP9lFYXvvYchM5O0IC1LRkpCvWt5tLQM/8gJMORBCB0BjfbDyC2tYUtKEVutR1GlZQZKpK+zJY16pDdxEV4426ne4CvZyR3sEp2quWvti1Q6ulL1n4XcNDSS5SnLeX7r84zuMpo3r30Tm+R18NUMiHsIJr7WpvtX1Rl56t/b6J1lpM99PRgVG9jwWPXevTgNHHjJBopbgYlSyvusf88EhkgpH7EGiqFAHbAeeLal/bKFEHOAOQAhISEDMzIy2l1vRWns6M4TbPrqKACjbu9O9yH+LZaVBgNlK1ZQMH8+prwTJIZoWD5Cy2DPOmbr83H07Q1xc6H3H8CmaavBbD4tjXq6nlqDGd3JNOpRlk2b+ga7Y6NVM9ivJOs/TyRldz6jCxYhkxMJ/fZb7MLDWHN8DU9vepoh/kN4f+z72FYVwYfDwS0I7lsPurPvnV1db+SPn+4iM7WUu8rtmDinNxEDfJuUuWBJAYUQ64Dm/o/5q5RyhbXM77QcKKYBE04LFIOllI8KIQKAE4AtsBBIlVK+fLZKqxaFcqGUF9Ww7tMj5KWWETXIj2tu74adY8uzmMz19ZR+8y0FCz5EFupJDIa11zhynS9MyUtB6+QDsfdC7Gxw8Wv2HnVGE3sySthqzYZ7IKcMKcHZTkdcuFdDV1WEj5Pate8yVlNZz2fPbiPELpewn/9B0Ntv4zpxAhuzNvL4b4/Tx6cPH437CEeNLXwxFXL2wAObwPvMLU1PV2swMXvxLnak6Xlzah9yP01h6M0RDJjQtUm58w0UZ23nSinHnetNT5MNdGn0dzCQa713nvVcnRDiU+Cpdj6XorSLq7cDU5/sz57VGez6+Ti5yaVce0d3Qvt6N1teY2uL51134j7tVkq/+w754XyivyzmaBD8dWwMU/2diNs4D7a8Bb1vhbgHIaBfk3vY6bQMi/BmWIQ3T0+wpFrYnqpns7Wbal1iPgABbvYNW8QOi/DGx+Xs3zKVS8eRLbmYjGZ8tn2M56xZuE6cwPbc7Tz5+5N09+zO/LHzcbRxhF//Csc3w80L2hQkKuuMPPjFHran6XlzWj9uHhDMom8yKCvquNxlF6PrSYdlMHsskINlMPsOKeVhIUSAlDJPWL4m/QeolVI+e/o9TqdaFMrFkH+8nN++SESfU0XUID9GTo/CweXMtAmNmevrKf3+e3I+eBddQQnJAXB0cnduifIn4tBKMFRB1+GWbqnuk0Bz9qmzWcXVbE62jm+kFlFabQCgm58zQ8O9GBrhTVy4J+6OrddN6Txmk5nPn92MbXYiw8Rmun62mL3FB5i7bi5dXLrwyYRPcLNzg4PL4Lt7YfAcmPTGWe+rr6xj9uJdHMot519/6MsfBlr2yv523m5s7bXc9Hj/JuU7ZT8KIcTNwHuAD1AKJEgpJwghAoGPpZSTrOUmAW8DWuATKeUr1vMbrNcKIAF4UEpZebbnVYFCuVhMRjN7f81g9y/HsXXQMXJ6FFGD/M7aBSTr6yn6fhnZH7yLfUEZqQGCgluHM7VPFJ57voCyTHAPsQx8978L7N0arl2+L4c3fj1KbmkNge4OPD2hO1P7W/bNMJklh3PL2JaqZ3uqnvh0yzaxQkDPAFeGRXgxLMKbQWGeamD8EpK8LYs1nyfT7/gShnz2KkmafO5fez++jr58OuFTvBy84MRB+Hg8BMZYcjlpW1+4mV1Szd2fxJNTUsP8OwYwrueprs01iw6Tn17GzH8Oa3KN2rhIUTpIcx/UI33d2PBFEgXHywmMcmfkbd3wDnY+672kwUDOsiXkffA+zoWVZPhrqZt5IxMHDcRu9yLI3Aa2zpZgMXgOyzPtee77g9QYTA33cLDR8totfRqCRWP1RjMHsksbAseezBLqjWa0GkHfYLeGwDGwq0ps2FmklCx9eBnV1TB9Thdyergz+9fZuNm6sXjiYvyc/KCyAD4eCyYDzNnY4njWSfHpxcz93x7qTWY+uWcQg0KbrtbeuTKNPaszeOC9axoyD4AKFIrSIU7m4mnug3pKv0CObMllx4pU6quN9L4mmME3hmHvdPaUHdJoJHnpIvQLFuBeWEOOvw12993F8GvGoNn1Xzj0HZiNbNbE8kHtdWw398TS0LYIcndg67Njzvo8tQYTezNKLIEjTc/+rFKMZomtVkP/EHfLWEikF/2C3bHVqRlVF8Ox+UtZe9CX/l2K8HmoP/f+ei92OjsWT1xs2WGxvgoW3wAFifDHnyFoYKv3+3JnBn9fcZgQT0cW3h1LpO+ZX1gSt+Wy4fMk7nw5DnffU7nIVKBQlA4wfN4GckrPHARs/EFdW2UgfmUahzblYOdow4CJXelzTRA627N/Y5dGI3v+9w7Viz7Hp7CeggAHvObcT6/rpyD2Lka/8UO8RAWJ5hA+MU1kpWkYddgigPR5k8/59VTWGdl1vJjt1hbHoVzLjCoHGy2xoR7WQXQvegW6olNTcTtc9d59/PTqb5R692TYCyE8tOlBbDQ2fDLxE7q6dgWzCb6+C46thtu+bDVFR1Wdkb+vPMyyPdmM7u7D2zP6t5hXLDe5hB/e3MeNj/YjpJdXw3kVKBSlA4Q9+zPN/R/R3Ad1UXYl275LJiuxBCc3W2InhxE9PKBJU78lRkM9Gxb/E93n3xNQaKLc1wm/OQ9y/3F/Ymo2Mlu7imhNFkXSlS9NY1nvdAMrn7u13a+vrNrAjnR9Q+A4mm/Z5MbFTseQcE/iwr0YHOZJzwAVONrLqNdzePq9bIl8lK7D7Jhn9yz2Wns+mfAJIa4hICWs+j+IXwjXv2FZ2d+C/Vml/GnpPjKLq3lkdCR/GtcNbTMZBE6qLKnjs+e2MmpGN/pcG9xwXgUKRekAbWlRnC7naAk7VqRxIq0MV297+l/XlR5x/m1qYVTWVfDj4udxXbqW8DwzlW6OfB02hh+DhzLANoXZ2lWM1ewDjRZNn1sts6UCOy4lWmFFHTvSLN1U21P1pBdZUo042WoZ0NWDIWGeDA7zom+wmxrjOAfSZCLz3vtIKA0jJ3A4ywa/js4ZFk1YRBcX62qBDf+ETW/A0EdgwivN3qfOaOKD31KZ/1sKvi52/Oe2GIaEezVbtsnzmyUfPfY7/cZ0YdgtkQ3nVaBQlA7Q2hhFc4PJJ0kpyTikZ9dP6RRkVODgYkPf0cH0via4TWMYuRW5LPvybwR+v51emVBub8cPYdeyO2Ycz4zxYHzFCtj3P8v02pBhloDRY3Kbpteei4LyWuKPFxOfbjmSTlhaHLY6DTHB7gwO82RwmCcDunqoWVWtKHjzLfIWL2HbyHkc89lHYt/1LJqwyDImAbDp37DhHzBgFtz4TpN0LyftSNPzlx8OklZYxU0xgbw8pTdurSz+PN0Xf9uGX5gb193bq+GcChSK0kFam556NlJKco+Vsm9tJhmH9OhsNEQM9KXniEACItzOOq32YOFBlnzzAr1/OcqAVIl0csD7zpl4zrobnaPWEix2Ljg1vXbwAzBgZpPptR2ptLqeXcdLiE+3TMU9lFuOySzRagS9Al0ZHGoJHINCPfFwUus4AMp++pncp54i+bq5ZNX3ZuOwT3nnln8R4BxgKbB9Pvz6F+h7G0z98Ixgn1VczZtrjrI8IZcung78c2ofrunm08wzte6HN/cizZJbnj41OK4ChaJcYvQ5lRzcmMOx+BMYak14+DvSY2gAkQN9W0w6CJZgsy5zHUtXvsawDScYmiQRtrZ4TJuO1+w/YuPnC0d/gR0fnppeG3Mna11u4sWtdecV4Nqqqs7I3swS4tOL2ZleTEJWKfVGc8PjjjZaJvcN4I/Dw+ju79JqP/pJ7QnMl5qag4fIuOsuKrqFsdn3Xsq8TvDws1PwdbTmXNr8Jqx/GaKnwK2fgvZUq6ykqp75v6Xw+fYMhID7RobxyOgoHNrQhdmcdZ8eISe5hFmvDm84pwKFolyiDHUmUvbkc2RLLifSygHw7epCxABfwvp54+7n2GxLo95Uz1dJX7F8/Qdct6WSUYclGqHB7aab8LrvPuzCwiA3AXZ+hPnAt2A2sd7cn09NE9lm7oWDje6sXWbt9e3uLP76wyHqTeYzHnOy1dKvizsDQjzoH+JO/xAPPE9rdZxvV9+lyFBQwPFp06k21zJ/7DD6nbiB65+KJjwywDJwve5F2Pq2ZQOiqR80LKjLKa3h481pLI3PotZoYtrAYJ4Y340At5a/TLTFjuWp7F2TyYPvX9uQOl8FCkW5DJQX1ZCyt4DUPQUUZFj6/5097AiO9iQk2pOASDecPZpmmy2pLeHD/R+yPv5rpsbD6AQTWqMZl4kT8Lr3Phx692LKa8sYW/Uzd2rX4S3KSTJ34TPTdexyHsu65859Wm1btTT47+Foww19A9mXVUJiXgUms+VzJszbif5d3Okd5EafYDce+2ofeWW1Z1zf1nUjlwpzXR0ZM2dSdTSJv9xhw9icfxLWzZ8pj/QHkxF+eQr2fGpJEDnp30gh2JlezFfxmQ1b6k6JCWTuNRFE+bl0SJ0Obcph45KjzHptOM4elrxgKlAoymWmXF9D1pFishKLyU4qoa7aCP/f3p3HR13eix7/fGef7PtOSFiCIKuAgrih4IIti0trq9VWrbWnvbXW1tve/qG3p+e0R6ve29N66lKv2lrFU7WAqBSsgpVFkDUQEtaEhJB9X2YyM8/9Y0aIkExiQhbC9/16zSszv9/zm9+TJ8l88+xAZKyDlJwYUnNjSMyIIi4tgpgkN0ebjvC7Hb/jk31/56btDq791I+1zUvE3Dk8KFPYnpyHUzpYbN3I3db3mGgpocm4iZ799eDqtWmTz/r30JvhxK1eH7tLG9hRUs/2kjp2lNRT3dzlbgJdXj/cmUCAYz96kJb3/s7jN1kYn/0AMTvHcevPZpGSCvz1bji4Fi57kGMzfsLbe07w+rZjHKluIdpl45aZWdx7+Rgy4/pXgzjd0T3VrP79bm5+eCZpY4J9WBoolDqHBQKGqpImKo40UHG0kYojjTRUnvpP3WIT4lKCAcPrbmFb02YONe9mZrmXeXubcdVUU+qO4c2xl/JR5jT8FgsXyQG+7f6QG2QT+D0w6pJgwJi09Iw9MvqqL8OJASoa28kva+CHy3fS1O4747zVIiyYmEJeavTJR25S5LCcTV78y0dp/fNy/jTfQtbdP8D52ula5wgAABn9SURBVCRSc2P40u1xmL98FaoP8OG4h/lNzaXsPR5sepydE89ts7NZNCW9z30QPakubWb5Lz/h2nsvZPys4JIgGiiUGmE8rR3UnWil7kQLdeXBr0217TTXeU7WPk4nxo/F106bRWh0OMlIjSIzwY6jtQRHfT4OzwkcjgCO0VNxjJ+HIyUTh8uGw23D4bKe/Grp5WS7/vYxdHW9zSJMyoih2eOjuKb1ZLOVzSLkJEWSlxrFuJRoRidEMDoxguzECJKjnEOyV0fR00/g/+3zrJltY8L//g+id+Wxa00J2ZfVcNWRH2P8HdzvfYCNgclclB3H9ZPTuP7CdLITI3p+837ytHbw/I8+4tKbxzFjYTYwgPtRKKWGhjPCTtqY2JPNBp11eP201Hlorm9na/F2/l70Pk2NLeRUJDGm3EVSh6HDGYXFkU1DIJEObxre9iS8bT5MC1AH7KwBarq8t81pxeGy4oq0E53oIibJTWySm5gkF/FpkcQmuxGLnAwGfR211NP17R1+Dle1cKCyiaKKJooqmtl3vJF380/Q+X/cCIeV7IQIshMiyIqPIDXGSUqMk9RoFykxTlJiXEQ7bWclmAQChrpWLx/9v/9gwu9f5ZMLHOy5/ud88EYk8w8VY3Ue5cYDD1FsyeL1cb/i5okzeGp8Eqkxg7s/usNtw+600lx3Zh/QF6U1CqVGAF/Ax6pDq3h619OcaDnBjS3j+Np2N46PdyIOB7HLlpL4rW9hz87G5w3grS7Hu/0tvPnv4W1sxGtNxpsxD2/qHLzOdLztfrxtPtpbfDTWtNFY1Ya3/dR//XanlcTMKJJGRZGaG0PG+DhiEs9uG3s4Hp+f0ro2SmpaKa5poaS2jZLaFoprWimrb6PV6z/jGofVQrTLRozbHvzqshPltOGwWbBZBKtFsFkFiwg+v8Hj8+PxBfD4ArR4fNS2eKlt8VLX2sZlzX/m4XV72J/h5H9d/EMSokexpNVPbLVwR9K/YJ0yn4ilTyLOnlcYHkh/eXQz8emR3PCdKYA2PSmlCA6pfevAWzyf/zwnWk5wVSCPu/ck4l67GePzEX3ttSTe/S3c00K77AUCUPxP2PEK7FsBvjZImgDTvw7TboPo4C7Ixhg8LT4aqtqoOd5MdWkz1ceaqC5tpiMUQKITXGTmxZE9OZHsCxNxuoeuwaLZ46OysZ2KRg+VTe1UNnqobvHQ1O6jqd1HY1sHTe0dNLX78AUMvkAAv9/gCxj8AYPdasFpt+C0WXDarLjtVhIiHUREtOLP/3fuXX6UqlHxpD29nDGj0mlZ/zrL30hhatQaLrt9Osy4fci+985W/J8d+Lx+bn44GBuGauOiW4FHgYkE98Hu8tNbRF4AvgRUGmMmdzqeACwHcoCjwFeMMXU93VcDhVLhef1eVhxawfO7n+d4y3HmOC7g/qIsolZvJNDYiHv6dBLuupPohQsRW+gDvb0R9r4FO/8CxzaDWGDcguC4/wk3QBf/HZuAoeZ4M2VF9Rw/UM/xonraWzqwWIXMvDhypyUzbmZKjzsDngs+KPmA11/+Gd99tQH/mCym/uVNrN5KzNsP8da2BdSZXG7/6QW4ssYOdVZPWvvCXsoPNXDnvwU3MBqqQDERCADP0M1WqKF0VwDNwMunBYrHgFpjzK9F5KdAvDHmf/Z0Xw0USp2pqxnON05NYeWhlTy35znKmsuYFpHH/WUTSFu9jY5jx7Clp5Nw+9eJu/VWrLGd+kKqD8LOV2D3cmgsA5sbJlwPk2+GcQu7HTUVCBgqDjdwZHc1R3ZVU1/RisUiZE9O5II5aeRMScJqH34jl8Jp7Wjlsa2PUbrqv3lgpcE2Jpe8F57FuvMZ2PQ0Be0L+Eftfcy/fQKTLh9ekwQ/fuMgez4o5Tv/eSUiMrRNT+H2zO6UJgd4+7RAUQhcFdo3Ox340Bgzoaf7aaBQ6vN6Gn3UEejg7UNv80L+CxxtPEpWRAbfa5vDhe8foX3LVsTtJm7ZUuLv+AbOMbmn3jgQCNYu8t+AvX+D1mpwRMPELwWDRu4VYHN2m6+asmYKN5+g8JMTtDZ4cUbauPCyDCZfmUV0wuB27vbF1hNbeWTjI+R9VMJ97/lxT5tG9v2X07HpSVyeal7vWEhZ3f1EJUdw7yNzkV4sWTKYdq4r4eO/HuTeJy/HGWE/ZwNFvTEmrtPrOmNMfDfX3wfcB5CdnT2zuLi43/lWaqTo7XyGgAmw/th6Xtz7ItsrtxPjiOEe9wLmb2rB8+5ajNdL5JVXkHDnnUReeunnRwn5fXB0QzBo7FsFnoZg0Bh3TXAl2/ELwd3lny+BgKF0fy17PzrOkZ1VAOROT2bq/CwyxscNydDWcOra63hi2xOsOPg3vrktikXrGqjISGHc3GNk2Wv4NJDHv3bcQWbzRC7osLI83sePvzJ52C07UvTJCda+sI+vPXIJCemRAxcoRGQdkNbFqZ8bY1aE0nzIAAeKzrRGodTnfZENlz6zq2oXL+a/yPsl72Oz2FiccCU3743CufID/NU1OMaNJf62rxG7ZDHW6NOWlfB54PCHsH81FL4LLZUgVhh9aTBojFsAieO6XD67saaNvRvK2PvP43hafKSNiWHmDTmMnpw45AHDH/Cz8tBKnvz0SbwtTfxqQzrpW47iHO0j9+JKdso4nvLdzIbAVPK8Vpa0OvnY1cFGl29YLjtSWljHiqd2sOSH08m6IGHg5lEYYxb0LYu9UiEi6Z2anioH8F5KjVgZce4uaxQZYZaFmJY8jafmP0VxYzGvFLzCqkOreCO1mQsfyuOeiisZs3Y/Fb/8JZVPPEHMjYuIv+1ruCeH9jawOSHvuuAjEIDj20NB4x1476fBNDFZMOaq0ONKiAquoBqT6GbusnHMvjGXgo3lbP97Mat/v5ukUVHMvD6HsTOSh6QJZ9PxTTz56ZPsr93PZb5MfvBKCxw/QvLUJvbkjeFR/7fYEJgKCJEBuLbNQbk1wGZncPLj8S7Kf6hFxgYHEbQ0ePv1PkPd9PQ4UNOpMzvBGPNwT/fTGoVSn3c2VmFt7Whl9ZHVvLb/NYrqioiyR3GHzOWq7V6s6zZi2tpwTZ5M/G1fJWbRIiwR3cwurj0Chz8I1jgOr4f2+uDxlEnBZURGXQxZF0PiWBDB7w9QtOUEn75XTENlG/HpkVyyOJcx05MHpYaxt3ovv9v+f/ln+SYyLS5+tq2BtA8EAkLGXZewpPZyDpv0k+ktBr7a7CDFb+HlaA911uBn6HCsUZw+O3uoRj0tA/4TSAbqgZ3GmOtEJAN43hizKJTuVeAqIAmoAB4xxvxRRBKB14FsoAS41RhT29N9NVAodaazta+DMYZdVbt4df+rvF/yPh6/h4mObO48lsP4Dw8ROFyMJTqa2C9/mdhlS3FNntz9B3rAD+W7gkHj6EdQug08wfWOcCdA1mxInwopkwgkT+LQ0Rg+eaeE+opWkrOjuWTJGLInJZz9gOFpZtu+5Tx34HU2th0nxh/g/qoGrv7YQWORE9fYLDJ/+zSOsePP6P+5ptXORV4bqyK87HcEA/NwXRrdGMMzP1jPlCszmXfLeJ1wp5Q6+xq9jaw9upaVh1ayvXI7YmBp6wVctxMSNheCtwPH2LHELl1C7OLF2FNTw79hIADVhXBsCxzbCqVboeYAmNB+FlYngaQLKPRczdaS2TS1uknPNMxZGEfGpPRgcOntgobeFmiugOZKqC+B6gN4qvazpqGA5aaR3S4nCX4/d/mcLG6YSt3KA/iq6ki4525SHngAcQSbbTrX1qZ5rFzb5uBTl48dCVDf2jHsN1vqvCWqBgql1IA61niMVYdX8e6RdznaeJSIdsOtpZlcvsdPzP4ysFiInDuX2KVLiJo/H2tUL5ev6GgPBo+KfVC5Fyr3Q30x/roy9jVdzrbmW2kNJDDKsYM50a+QEnEiuPWrzQFWB1jsYLGAzxucWe7zgKc5uL84YIACh513oqJYER1NvQVybDHcNuoavpzwJRqf+TNNa9fiGDeW9F/8goiLLjoji3/bUcbyNwuZV2U47hZm35HHsplZZ7F0B84bj23Darey9MEZGiiUUoPncP1h1pWsY13xOgpqC0itNdxYGMFle/xE1bSC3U7UvHlEX3cd0VfP//xkvt4KBKClko7Ko+RvqmP7VhvtHiu56ZXMzt1NclQ1+L3g7wg2cdmcYHOB3YXfFsF+u4V/dFSzpvEAxe1V2MTK/Oyr+eqErzKDbGr+8Az1f/0r4nSSdN+3Sbz77pO1iNMd/LSStS/sJTU3hsU/mI5tgJYGHwjvPrOHuvIWvv7oHA0USqmhUdZcxsdlH7O5fDNbjm8i7UgTlxQGuKzISny9D2O1wMypJF1zLbFXXIUjJ6dPfQ7eNh+7/nGMneuO4W3zkTMlkVmLcknNjcEX8HGo/hB7qvewpXwLm8s3U++pxyIWZqfN5oacG7gm+xqcB0qpfeklGt97D0SI/8pXSPru/diSkrq9b/76Uta/VkT6mFhu/N5UnBH2/hTXoNvwaiFFWyu498krNFAoNdjOVufxSOIP+NlXs48tJ7awu3IXDbu2M3F3HRcXGtJDq7g1Jrqon5aDmTWViBkzSMrOIyUihRhHDHaLvcsgYoyhzddGnaeO6vpaCtafoOYTwGOlPrmUjWmrKIncDwLJ7mTmZszl0oxLmZsxl+g6D43vvEPDylV4CguxREYSd8stxH/jGziyuv95dXj9fPRaEQUbyxk9JZHrvj0Z+zlUk/jM1tVH+GTVEe7/3VXY7FYNFEoNlrMxHPV8YIyhvKWc3dW7Kd3/KWbLdhJ2lZB7sBl3aGh/VQwUZQoHM4SyFCs1qRG0x7sxEgw8voCPdn87HYGOz7233e9kdvW1TCydh93rxpEcYMLlqczMiydwqJDWTz6hZeNGPAcOAuCaNpXYLy8mdumSHvtPju2vZcOrRdRXtDLzhtFc/OUxWIbZ8hy9lb++lPWvFvHNX88jKt6lGxcpNVgeX1P4uSAB0Nbh5/E1hRooOhERMqIyyIjKgJzr4frg8Y72Nip3bqbu0y1E7d7D7L0HmVfQSHCN0Qa87laaU6Noj4vAEx+PLz4Ke1QMrogYXJExxNljiLdEYE8wtJsCjpbZOHQikz1vWijwFZNctYvUut1kjkshZckSohcuxDF6dNi8GmMoP9jA9jXFFOfXEJPkYvEPpjNqUsKAl9NAckUF+13amvs+6U4DhVJ90N0s3OE4O3c4srvcZM6ZT+ac+SeP+aqr8Rw6jOfQQbwHDxF/7Bi+qip8uyrx1xad8R6f7UcgDgepaWlkpqXRlDGFEsdEytxzKU+fyz63jfT2WNILDImN1cQmu3FHO7DaLfg7ArQ2eqk93sKJIw0c2VlFY3U7rig7lywZw/RrRp1TndbdcUcH+1Tamjp6SNk9DRRK9UFflsxQ4dmSkrAlJRF5ycVnnDM+H6u2HOa//r6P2tomkuMi+M6Cidx4cS7icn2uX2MK4O8IUFJQy9E91ZQfqKc4v+stXz9jsQlZefHMWpTDuFmp52RfRHfcWqNQamj85LoJXfZR/OS6HlfJV32wYk8FP3vvMG0ddohMoKIDHl5Xgj8mtsumPqvdQu7UJHKnBkcztTd3UFfRSmNVK+0tPvy+ABarEBHjIC41gsSMqHNun4ze0hqFUkPksw8nHfU0OPrSJ9TlqLQ559/PxxlpBwkGy77SQKFUHy2dkamBYZB80T6h00elldW38bM39wCngvz5MrzZYhFckXbamvre9DQy61pKqRGlu76f7o6Hq4HAqUBSVt+G4VQg+duOsrOa7+HCHWWnrR81Cg0USqlh7yfXTcBt/3wHc7g+oZ5qID0FkpHGHe3QGoVSamRbOiOTX900hcw4N0Jw74dwkxt7qoGcb8Ob3VF27aNQSo18X6RPqKdRaefb8GZXtIO2ovo+X9+vGoWI3Coie0UkICLdTgsXkRdEpFJE8k87/qiIlInIztBjUX/yo5RS0HMN5Is2ZZ3r3FF22luHrkaRD9wEPNNDuheB3wEvd3HuKWPMb/qZD6WU+pxwNZDzbXizO9oe3Jijj/oVKIwxBUCPSwYbYzaE9sxWSqlh4Xwa3vzZ7Oy+Gg6d2d8Xkd2h5qn47hKJyH0isk1EtlVVVQ1m/pRS6pz22ezsvuoxUIjIOhHJ7+KxpF93DvovYCwwHSgHnuguoTHmWWPMLGPMrOTk5LNwa6WUOj+4ovoXKHpsejLGLOjXHcK/d8Vnz0XkOeDtgbqXUkqdr1yRA1yjGEgikt7p5TKCneNKKaXOoiENFCKyTERKgbnAahFZEzqeISLvdEr3KrAJmCAipSJyT+jUYyKyR0R2A/OBB/uTH6WUUmeyOazkXZLa5+t1K1SllDpPiEiftkIdDqOelFJKDWMaKJRSSoWlgUIppVRYGiiUUkqFpYFCKaVUWBoolFJKhaWBQimlVFgaKJRSSoWlgUIppVRYGiiUUkqFpYFCKaVUWBoolFJKhaWBQimlVFgaKJRSSoWlgUIppVRY/d246FYR2SsiARHpco1zERklIh+ISEEo7QOdziWIyFoRORD6Gt+f/CillDr7+lujyAduAjaESeMDHjLGTATmAN8TkUmhcz8F3jfGjAfeD71WSik1jPQrUBhjCowxhT2kKTfGbA89bwIKgMzQ6SXAS6HnLwFL+5MfpZRSZ9+g9lGISA4wA9gSOpRqjCmHYEABUgYzP0oppXpm6ymBiKwD0ro49XNjzIre3khEooA3gB8aYxp7n8WT198H3AeQnZ39RS9XSinVRz0GCmPMgv7eRETsBIPEK8aYNzudqhCRdGNMuYikA5Vh8vEs8CzArFmzTH/zpJRSqncGvOlJRAT4I1BgjHnytNMrgbtCz+8Cel1DUUopNTj6Ozx2mYiUAnOB1SKyJnQ8Q0TeCSWbB3wDuFpEdoYei0Lnfg0sFJEDwMLQa6WUUsNIj01P4Rhj3gLe6uL4cWBR6Pk/Aenm+hrgmv7kQSml1MDSmdlKKaXC0kChlFIqLA0USimlwtJAoZRSKiwNFEoppcLSQKGUUiosDRRKKaXC0kChlFIqLA0USimlwtJAoZRSKiwNFEoppcLSQKGUUiosDRRKKaXC0kChlFIqLA0USimlwtJAoZRSKqz+7nB3q4jsFZGAiMzqJs0oEflARApCaR/odO5RESnrYuc7pZRSw0S/drgD8oGbgGfCpPEBDxljtotINPCpiKw1xuwLnX/KGPObfuZDKaXUAOnvVqgFACJd7nT6WZpyoDz0vElECoBMYF+3FymllBo2+luj+EJEJAeYAWzpdPj7InInsI1gzaOum2vvA+4LvfSISP4AZvVckgRUD3Umhgkti1O0LE7RsjhlQl8uEmNM+AQi64C0Lk793BizIpTmQ+DHxphtYd4nClgP/Jsx5s3QsVSCP0AD/CuQboy5u8dMi2wzxnTZJ3K+0bI4RcviFC2LU7QsTulrWfRYozDGLOhblk4RETvwBvDKZ0Ei9N4VndI8B7zd33sppZQ6uwZ8eKwEOzD+CBQYY5487Vx6p5fLCHaOK6WUGkb6Ozx2mYiUAnOB1SKyJnQ8Q0TeCSWbB3wDuLqLYbCPicgeEdkNzAce7OWtn+1PvkcYLYtTtCxO0bI4RcvilD6VRY99FEoppc5vOjNbKaVUWBoolFJKhTWsA4WIXC8ihSJyUER+2sV5p4gsD53fEpqnMSL1oix+JCL7RGS3iLwvIqOHIp+Doaey6JTuFhEx3S0vMxL0pixE5Cuh3429IvKXwc7jYOnF30h2aDmhHaG/kxG5ZJCIvCAild3NNZOg34bKabeIXNTjmxpjhuUDsAKHgDGAA9gFTDotzb8Afwg9vw1YPtT5HsKymA9EhJ5/93wui1C6aGADsBmYNdT5HsLfi/HADiA+9DplqPM9hGXxLPDd0PNJwNGhzvcAlcUVwEVAfjfnFwHvAgLMAbb09J7DuUZxMXDQGHPYGOMFXgOWnJZmCfBS6PlfgWsk3Hoi564ey8IY84ExpjX0cjOQNch5HCy9+b2A4ATOx4D2wczcIOtNWXwb+L0JrXhgjKkc5DwOlt6UhQFiQs9jgeODmL9BY4zZANSGSbIEeNkEbQbiTpuqcIbhHCgygWOdXpeGjnWZxhjjAxqAxEHJ3eDqTVl0dg/B/xhGoh7LQkRmAKOMMSN9Amdvfi/ygDwR+VhENovI9YOWu8HVm7J4FLgjNKT/HeB/DE7Whp0v+nkyuGs9fUFd1QxOH8vbmzQjQa+/TxG5A5gFXDmgORo6YctCRCzAU8A3BytDQ6g3vxc2gs1PVxGsZX4kIpONMfUDnLfB1puy+BrwojHmCRGZC/wpVBaBgc/esPKFPzeHc42iFBjV6XUWZ1YVT6YRERvB6mS4Kte5qjdlgYgsAH4OLDbGeAYpb4Otp7KIBiYDH4rIUYJtsCtHaId2b/9GVhhjOowxR4BCgoFjpOlNWdwDvA5gjNkEuAguGHi+6dXnSWfDOVBsBcaLSK6IOAh2Vq88Lc1K4K7Q81uAf5hQb80I02NZhJpbniEYJEZqOzT0UBbGmAZjTJIxJscYk0Owv2axCbNg5TmsN38jfyM40AERSSLYFHV4UHM5OHpTFiXANQAiMpFgoKga1FwODyuBO0Ojn+YADSa4HUS3hm3TkzHGJyLfB9YQHNHwgjFmr4j8AthmjFlJcA2pP4nIQYI1iduGLscDp5dl8TgQBfx3qD+/xBizeMgyPUB6WRbnhV6WxRrgWhHZB/iBnxhjaoYu1wOjl2XxEPCciDxIsKnlmyPxH0sReZVgU2NSqD/mEcAOYIz5A8H+mUXAQaAV+FaP7zkCy0kppdRZNJybnpRSSg0DGiiUUkqFpYFCKaVUWBoolFJKhaWBQimlVFgaKJRSSoWlgUIppVRY/x/WUXtLPsuwQgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7d90152d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Prepare the data\n",
    "N = 25\n",
    "degree = 15\n",
    "sigma = 0.03\n",
    "x = numpy.sort(numpy.random.uniform(0.0, 1.0, N))\n",
    "def f(t):\n",
    "    return 1.5*t**3 - 2*t**2 + 0.5*t - 1\n",
    "y = f(x) + numpy.random.normal(0.0, sigma, N)\n",
    "X = [[t**i for i in range(degree)] for t in x]\n",
    "\n",
    "# Solve a number of lasso-regularized models\n",
    "print( \"Sparsity report:\\nlambda1  nonzeros in w\\n\")\n",
    "plt.axis([0.0, 1.0, -1.25, -0.80])\n",
    "plt.scatter(x,y)\n",
    "ticks = numpy.linspace(0.0, 1.0, 10000)\n",
    "lambdas = [10**(-i) for i in range(1,6)]\n",
    "\n",
    "for lambda1 in lambdas:\n",
    "    M = lseReg(X, y, lambda1, 0)\n",
    "    M.solve()\n",
    "    w = M.getVariable(\"w\").level()\n",
    "    print(\"{0: <8} {1: <2}\".format(lambda1, sum([0 if abs(val)<1e-5 else 1 for val in w])))\n",
    "\n",
    "    # Plot the data points and regression function\n",
    "    plt.plot(ticks, sum([w[i]*ticks**i for i in range(degree)]))\n",
    "\n",
    "plt.legend(lambdas, loc='upper right')\n",
    "plt.show()    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Robust regression with Huber loss function\n",
    "\n",
    "The penalty function $\\phi(x)=x^2$ grows quite rapidly for large values of $x$, making the least-squares approximation overly sensitive to outliers. One solution to this issue in the realm of convex optimization is to replace it with the penalty defined as:\n",
    "$$\n",
    "\\phi_{\\mathrm{Huber}}(x) = \n",
    "\\begin{cases}\n",
    "x^2       &  |x|\\leq M, \\\\\n",
    "M(2|x|-M) &  |x|\\geq M,\n",
    "\\end{cases}\n",
    "$$\n",
    "for some value of $M$. This *Huber loss function* agrees with the quadratic loss for small values of $x$ (small residuals) and otherwise it is the slowest-growing function preserving that property while remaining convex. Thus it assignes much lower loss to distant outliers than pure quadratic regression, maing it more robust. The *robust regression with Huber loss function* is now the problem\n",
    "\n",
    "$$\\mathrm{min}_{\\mathbf{w}\\in\\mathbb{R}^d} \\phi_{\\mathrm{Huber}}(\\mathbf{y}-\\mathbf{X}\\mathbf{w}).$$\n",
    "\n",
    "It is left as an exercise to prove that an equivalent formulation of this optimization problem is\n",
    "\n",
    "$$\n",
    "\\begin{array}{rl}\n",
    "\\mathrm{minimize} &            \\|\\mathbf{u}\\|_2^2+2M\\|\\mathbf{v}\\|_1 \\\\\n",
    "\\mathrm{s.t.}     &            -(\\mathbf{u}+\\mathbf{v})\\leq \\mathbf{y}-\\mathbf{X}\\mathbf{w} \\leq \\mathbf{u}+\\mathbf{v}\\\\\n",
    "                  &            0\\leq \\mathbf{u}\\leq M\\\\\n",
    "                  &            0\\leq \\mathbf{v} \\\\\n",
    "                  &            \\mathbf{u}, \\mathbf{v}\\in \\mathbb{R}^n, \\mathbf{w}\\in\\mathbb{R}^d\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "We proceed directly to a *Fusion* implementation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Robust regression with Huber loss function\n",
    "def huber(X, y, MConst):\n",
    "    n, d = len(X), len(X[0])\n",
    "    M = Model(\"LSE-HUBER\")\n",
    "    \n",
    "    # The regression coefficients and other variables\n",
    "    w = M.variable(\"w\", d)\n",
    "    u = M.variable(n, Domain.inRange(0.0, MConst))\n",
    "    v = M.variable(n, Domain.greaterThan(0.0))\n",
    "    \n",
    "    # The residual and bounds on its absolute value (coordinatewise)\n",
    "    r = Expr.sub(y, Expr.mul(X,w))\n",
    "    ab= Expr.add(u,v)\n",
    "    M.constraint(Expr.add(ab,r), Domain.greaterThan(0.0))\n",
    "    M.constraint(Expr.sub(ab,r), Domain.greaterThan(0.0))\n",
    "    \n",
    "    # t >= |u|_2^2 and the objective\n",
    "    t = M.variable()\n",
    "    M.constraint(Expr.vstack(0.5, t, u), Domain.inRotatedQCone())\n",
    "    M.objective(ObjectiveSense.Minimize, Expr.add(t, Expr.mul(2*MConst, Expr.sum(v))))\n",
    "    return M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Huber regression example\n",
    "\n",
    "We demonstrate the Huber regression compared to standard least-squares regression on an example with distant outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD8CAYAAABkbJM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmczWX7wPHPNTOWIYwtMYjKvmSZlFRKQjwxRNEjLaRNafNE+9PvqRTtUdGiXSmGoqayp5Q9Y8ugYgzZhjAMM9fvj/s7GpwxhzMzZ87M9X69zmvO+X7vc851Rp1rvvdy3aKqGGOMMTkJC3YAxhhjQoMlDGOMMX6xhGGMMcYvljCMMcb4xRKGMcYYv1jCMMYY4xdLGMYYY/xiCcMYY4xfLGEYY4zxS0SwA8hNlSpV0lq1agU7DGOMCSmLFi3arqqVc2pXqBJGrVq1WLhwYbDDMMaYkCIif/jTzrqkjDHG+MUShjHGGL9YwjDGGOMXSxjGGGP8YgnDGGOMXwJOGCJSQUS+E5G13s/y2bR7VkQSvNu1WY6PE5ENIrLUuzXzjouIvCIiiSLyq4i0CDRWY4wxpy43rjCGAtNVtQ4w3Xt8FBHpArQAmgHnA0NEpGyWJkNUtZl3W+oduxKo490GAq/nQqzGGGNOUW4kjG7Ae97994BYH20aArNV9bCq7gOWAZ38eN331ZkPRIlI1VyI1xhjzCnIjYRRRVWTAbyfp/toswy4UkRKiUgl4DKgRpbzT3ndTi+KSAnvWDSwMUubTd4xY4wxQeDXSm8R+R44w8eph/15vqp+KyLnAT8C24CfgMPe6WHAFqA4MAZ4EHgSEF8v5SO2gbguK2rWrOlPOMYYY06BXwlDVdtnd05EtopIVVVN9rqM/srmNZ4CnvKe8zGw1jue7DU5KCLvAg94jzdx9FVIdWCzj9cdg0s0xMTEHJdQjDHG5I7c6JKaAtzg3b8BmHxsAxEJF5GK3v2mQFPgW+9xVe+n4MY/ErK8bj9vttQFwO4sycUYY0w+y43ig8OBz0SkP/An0AtARGKA21R1AFAMmOtyAnuAvqqa2SX1kYhUxnVBLQVu845PAzoDicB+4KZciNUYY8wpEtXC04sTExOjVq3WGGNOjogsUtWYnNrZSm9jjDF+sYRhjDHGL5YwjDHG+MUShjHGGL9YwjDGGOOXQrWnd2EXtySJEfFr2JySSrWoSIZ0rEdsc6uWYozJH5YwQkTckiSGTVxO6qF0AJJSUhk2cTnACZOGJRljTG6xLqkQMSJ+zZFkkSn1UDoj4tdk+5zMJJOUkoryT5KJW5KUx9EaYwojSxghYnNK6kkdh1NLMsYYkx1LGCEiqlQxn8erRUVm+5xTSTLGGJMdG8MIAXFLkji4L5Uqe1OouH8P5Q78TZgqYeFhDDynDixfDmecARUrQtg/fwNUi4okyUdyOFGSMcaY7FjCKEAyB6h3bttFp13ruC0siXrbfifmh4Uk7Ewm7PjtQODjLPcjIuDss6FxY2jShBEVz+auXaXZIcWPNIksFs6QjvXy/sMYYwodSxgFxJffL2fp02N5bvUPxGxNoETGIQ5LODurn83SCufwed1L+atceXZGliUlsgzpEka4KuMHtIIdO2DLFti8GdasgWXLYOJELlRlQXg4K6vV5dszW7Aoph29ru1os6SMMafEEkY+S0uDBQvcbeUKpcyP8XRIHE2ntK+5isMspzGvMJjvac8PehH7N5Z2T0yEsFIHiSi3n4io/RSv/DdVzz7AX42bcbqvTXH37YMffyRs1iwaz5hB43kfww8fQXwDuO46uPlmqFYtXz+7MSa0WXnzfJCYCF98AfHx8NNPcPjAIXoznqHhI2iUvpyUUlV5L6oLE2peyZrKtSBckTD376JpEXRvVItJ87dxYFcJDu+O5FBKadJ3lzry+nXrQrt27nbFFRAV5SOIzZth0iT47DOYMwfCw+Gqq+COO6B9exBfO+IaY4oCf8ubW8LII9u2wbhx8NFHrocI4NymyqBaX3Htwgcos/k3aNQI/vMf6N2bNi/84HOAOjoqknlD2x23AO+OC+sTnVGNBQtg9mx327sXihWDyy+HHj3crWJFH8ElJsLYsfDuuy7Qli1h6FDo3t0lEmNMkWIJI0h+/hleeQU+/9x1P7VuDb16wbWNEqj27GCYMQPq1YNnn4WuXY/8ZX/sSm5wA9TP9Gji15jDoUPwyy8webK7mlm/HooXd29x883QoYOPXHDwIHzwATz3HKxdCw0awNNPQ7dudsVhTBGSLxsoiUgFEflORNZ6P8tn0+5ZEUnwbtdmOT5XRJZ6t80iEucdv1REdmc591ggceaHefPcl/IFF8BXX8Gtt8KKFfDjnMPcu/8pqv2rBWmLFvPCVXdR56pnabOqDHFLNx95fmzzaJ7p0YToqEgEd2Xhb7IAd2XRpo377k9MhMWL4fbbYeZM6NwZzjkHXngBdu/O8qQSJWDAAFi1Cj79FFTdVUabNvDDD7n7CzLGhD5VPeUb8Bww1Ls/FHjWR5suwHe4AfbSwEKgrI92XwD9vPuXAl+dbDwtW7bU/LZihWrHjqqgWrmy6nPPqf79d5aTLVuqgm7s0FUvuG+8nvngV0du9R/5Wict3pSn8R08qDphguoll7gYTztN9b77VJOTfTQ+dEh17FjVatVc4759s2lojClMgIXqx3dsoCu9uwHvefffA2J9tGkIzFbVw6q6D1gGdMraQETKAO2AuADjyTcpKXDPPdC0qeuGGjECNmyAIUPgtNNwgxfnnQd//gkTJnDtZYNJLnbaUa+RH2U6iheHnj3dGMfCha636eWXoXZtGDzYjYUfERHhrjjWroVHHnED5PXquSekp2f7HsaYoiHQhFFFVZMBvJ++JnguA64UkVIiUgm4DKhxTJvuwHRV3ZPlWGsRWSYiX4tIowDjzFWTJ0P9+m6sYsAA+O03eOABKF0aN3AxaBD07QsxMfDrr9CzZ4Eo09GyJXz4Iaxe7WbWjhrluqoeeQT2ZP3NlyoF//d/kJDgBmHuuQcuvtit8TDGFFk5JgwR+T7L+EPWWzd/3kBVvwWmAT8CnwA/AYePadbHO5dpMXCmqp4LvMoJrjxEZKCILBSRhdu2bfMnpFO2axdcfz3ExkLVqu4v9jfegMqVvQY7drgpSqNGuQzy/feuZAfZl+MIRpmOc86Bt992ia57d3jqKbdAfNQoOJz1X6ZOHfj6azcwvno1NGsGI0dCRka+x2yMCb6AZkmJyBrgUlVNFpGqwCxVPWHdCRH5GPhQVad5jysCvwHRqnogm+f8DsSo6vYTvXZezpKaNw+uvRa2boWHH4aHHoJpK/6Z6tpCdzPus8cok7wJ3nvPNc4i0FlQeWnRIteVNnMmnHuuSxxt2hzTKDnZjaJPnuwWfHzwgS38M6aQyJdZUsAU4Abv/g3AZB+BhHtJARFpCjQFvs3SpBdugPtAluecIeLmdYpIKy/OHQHGekpU4cUX4dJLoWRJmD8fnnjCJYvMvSYabl3H66MGkbFlK3NHfXxcsoDAZ0HlpZYtYfp0mDABdu6Eiy6CG290F0xHVK3qFv699Zb7JTRtClOmBCtkY0ww+DMynt0NqAhMB9Z6Pyt4x2OAt7z7JYGV3m0+0OyY15gFdDrm2CBgBW78Yz5woT/x5PYsqb17VXv2dBOGYmNVU1L+OXfhM9P1zAe/0qv6vaC7S5TWTWUq6+X9R+uFz0zP1Rjy2969qsOGqUZEuFlfn3yimpFxTKPVq1WbNXO/mAcecLOrjDEhCz9nSdnCvWwkJ7vKGUuWuDV2999/9Fq22kOn0iT5Nz789FF2RZahT5+n2Vz2dATYMLxLrsQQTL/+6gb0Fyxwv4cxY44MxzgHD8J998Ho0dC2LYwff0wDY0yoyK8uqUIpIcEtwFu92nXZP/DA8Quf2/39Bx9++igpJU87kiyg8Ow10bSpq3s1ciR8952rmP7FF1kalCjhBjs++MAtMW/Rws0vNsYUWpYwjvHDD27A99AhV6PvX//y0WjlSl7/8CF2R5ah93XPHEkWhW2vifBwd2W1eDHUquXWc9xwA/z9d5ZGffu6RBEZ6a40PvwwWOEaY/KYJYwspk+Hjh3d+O7PP7s/mo+TlASdOlG8VCSrPpiI1DyzwA1i57YGDdzVxqOPunzQooWbWXVEkybuF3bBBW7e8bBhNvXWmELIxjA806a56q5167oumCpVfDRKSYFLLoHff3eXH82aBRRvKJo71y3627rV1a0aPDhLd11aGtx1lxvw6NnTdVeVLBnUeI0xObMxjJMwd65bjNeokVuL4DNZpKW5jLJ6NUycWCSTBbgF30uXwpVXwr33ukq8R1aJFy/uVjKOHOnK9bZvf8zcXGNMKLOEgSv5dPfdrkvK5/4R4L4dZ850S6Tbt8/X+AqaihUhLs7lhbg49/tbscI7KeIGPj791E2xatPGFdkyxoQ8Sxi4XpORI7PZqQ7cYrXRo910qeuvz9fYCqrMvDB9uiuZ3qqVW/h3xDXXuNIof/3lVgIeySjGmFBlCSMnP/7otjHt0AGGDw92NAVO27Zurcq557oc8eijWca7L77YjfWourEfm3ZrsohbkkSb4TOoPXQqbYbPIG5JUrBDMjmwhHEiW7fC1VdDzZpuYZptX+pT1aqut+7mm+F//3MFDY9MvW3c2M1VLl/eFWacPj2osZqCIbO2WlJKKgokpaQybOJySxoFnCWM7GRkQL9+bmbUpEnuC89kq0QJ13P3yiswdarrhdq40Tt51lkuaZx1FnTp4qakmSJtRPyaowpxQv7sD2MCYwkjOyNHwrffus2DmjQJdjQhQcTNqp061c08btXKlYAHXNmQmTPdVLTYWLeE3hRZBWF/GHPyLGH4Mn++q2Heqxfcckuwowk5HTu6oZ8SJdzQRVzmbiYVK7ouqRYtXFffUaPkpigpSPvDGP9ZwjjWnj3Qpw9Ur+4WoB1bRMr4pVEjN8bdpInLDaNHeyeiotzKyNat3e/5qAJVpqgY0rEekcWOHhMsbKV1CiNLGMe6/363D/fHH59gnq3xR5UqMGMGdO4Md94JQ4d6M6jKlHHjGOefD717u4WQpkgpyPvDmOxZaZCsvvnGLWEeOhSeeSb3AiviDh9225y/+aZbxvL221CsGO5qrlMnt8Dv88+hm1+7/hpjcpmVBjlZKSluA4hGjdyWeibXRETA66/Dk0+68lI9esD+/UDZsi5Jt2zpFnHExwc7VGPMCVjCyHTvvbBlC4wb50ZrTa4ScYv6Ro92s6g6doRdu3BJ4+uvoWFDt4BjzpyjnmeLu4wpOAJOGCLSS0RWiEiGiGR7SSMinURkjYgkisjQLMdri8jPIrJWRD4VkeLe8RLe40TvfK1AY83WzJkuUQwbBjE5XpWZANx+u1sD+fPPcNllrnII5cu7Kcy1arl1Gt6KcFvcZUzBkhtXGAlAD2BOdg1EJBwYBVwJNAT6iEhD7/SzwIuqWgfYBfT3jvcHdqnqOcCLXru8cfHFrs/k0Ufz7C3MP665Br76Cn77zU273bgRqFzZ1Z46/XQ3jrRihS3uMqaACThhqOoqVc3p/+BWQKKqrlfVNGA80E1EBGgHfO61ew+I9e538x7jnb/ca5/7IiLgtttceW6TLzp0cBcVyckuXycmAtWquaRRsiR06ID88bvP59riLmOCI7/GMKKBjVkeb/KOVQRSVPXwMcePeo53frfX/igiMlBEForIwm3btuVR+CYvXHSR6w3cu9ddaaxeDdSu7TJJaiqfTHiMSvt2Hfc8W9xlTHD4lTBE5HsRSfBx83cepK8rAz3B8RM95+gDqmNUNUZVYypXruxnOKagaNECZs1y6zPatoWEBFzBwqlTqbpvJ+9PeJzTDu4/0t4Wdxnjww8/wObNef42fiUMVW2vqo193PwtCLQJqJHlcXVgM7AdiBKRiGOOH/Uc73w5YKef72dCSOPGLmlERMCll7py6bRuTcSkidTf8SfvfjmcEocP2eIuY3xRhX//23Wr57H86pJaANTxZkQVB3oDU9StGpwJ9PTa3QBkJqEp3mO88zO0MK0yNEepXx9mz4ZSpVwV9CVLgE6dCHvnHc5bt5g1SZ8y7z+XWrIw5liLF7vqFD165Plb5ca02u4isgloDUwVkXjveDURmQZHxiAGAfHAKuAzVc3cgu1B4D4RScSNUbztHX8bqOgdvw84MhXXFE7nnOOSRpkyLmksXoxbGv7ss24u7gMPBDtEYwqeiRPdXj1XXZXnb2WlQUyBs2GDW6OxZ4+bNNWiucI997jNNl56CQYPDnaIxhQcDRu67QNmzDjll7DSICZk1a7txjTKloX27WHpMoEXXnArwe+914oVGpNpzRpYtSpfuqPAEoYpoGrVclNuTzvNJY3lK8Phww9dhdt//xt++inYIRoTfJMmuZ+xsSdul0ssYZgCq3ZtlzRKlnRjGit/LwVTpkB0NHTtCuvWBTtEY4Jr4kQ47zy3f08+sIRhCrSzz3ZdsxER0K4d/LarsitWmJHh6k7tOn5hnzFFwqZNbmuA7t2Jj/eqJeQxSximwKtb1+3smpHhrjQ2RNRx+75u2OD6btPSgh2iMfnP2/v40FU9uOEGt/dbXrOEYUJCgwZuZ9d9+1zS2FT7YrcT06xZcOutbvGSMUXJpEnQoAHfbKjH1q3Qv3/OTwmUJQwTMs4915WZ2rHDJY1xZ1zO25f3g3HjGN15oJU9N0XHtm3uj6UePXjnnX+KPOc1SxgmpMTEuO3A//gzg9v6luaJRn2Y0uASbvvmbb5/5k1LGqZomDgRMjLYeXkvvvoK+vXztj3OY5YwTMhp0wbOvm4ZB7eX5q/Pz+f+9vexrGpdnot7jrh3vwp2eMbkvQkToG5dxi1uyuHDcNNN+fO2ljBMSNpfeTOVuy4hLbkcm6a04ZZuj7KrZFmeefdht8mGMYXVX3/BzJlor2t4+x3h/PPdYu/8YAnDhKRqUZGUqruVil1+5cAflVg5/XL693iMqIP73IrwAweCHaIxecPrjkpo0IuVK+Hmm/PvrS1hmJA0pGM9IouFc1qjJCpckUDq2jP4aXF3lj31itsT/JZbbOaUKZw++wzq1WP03CZERsK11+bfW1vCMCEptnk0z/RoQnRUJGVb/EGNK9azZ3k0H2+8Bf3vk66MyIgRwQ7TmNy1dSvMnk1a7DV8/Ilw9dVQrlz+vX1Ezk2MKZhim0cf2R9DFYYMgeefh4qPP8Lj1yTA0KHQqJFbEW5MYeB1R00r3Ys9e9yFdH6y8uam0FB1i5fefRdef34/t314kas3NX++W/lnTKi77DLYsoULyq5kz9/CihUgvjazPklW3twUOSIwZgx06wa331+KuBvjXOXCbt2s5pQJfZs3w+zZbLnkGn7+RRg4MHeSxckIKGGISC8RWSEiGSKSbXYSkU4iskZEEkVkaJbjH3nHE0TkHREp5h2/VER2i8hS7/ZYIHGaoiMiwm3O17Yt9Lq/JvOHfAG//w59+kB6erDDM+bUffopqPLmnj6UKOEW6+W3QK8wEoAewJzsGohIODAKuBJoCPQRkcxZwx8B9YEmQCQwIMtT56pqM+/2ZIBxmiKkZEmYPBmaNIHLH7+I9fe9BvHx8NBDwQ7NGL/FLUmizfAZ1B46lTbDZ7DrrfdIb9aC56fW55proEKF/I8poIShqqtUdU0OzVoBiaq6XlXTgPFAN+/509QD/ALkT1F3U+iVK+eqoFetCueNHcjOa2+D555zlx/GFHBxS5IYNnE5SSmpKFB8fSLlVy5jcuXu/P03DBwYnLjyYwwjGtiY5fEm79gRXlfU9cA3WQ63FpFlIvK1iDTK+zBNYVOliitWWLw4xPzwMgfOu8itclq6NNihGXNCI+LXkHrony7UbitnkYEwNKEnDRu68jjBkGPCEJHvvTGGY2/d/HwPX8Myx07NGg3MUdW53uPFwJmqei7wKhB3gvgGishCEVm4bds2P0MyRcVZZ7neqJ17i9N+1+ekR1Vw21lu3x7s0IzJ1uaU1H8eqNJ15WzmndGctcn1ue22/B/szpRjwlDV9qra2Mdtsp/vsQmokeVxdWBz5gMReRyoDNyX5T33qOpe7/40oJiIVMomvjGqGqOqMZUrV/YzJFOUNG3qdnZduLEK/SvEkZ68hYUXdOCc/0yhzfAZVuHWFDjVoiKP3G+yJZGzdm3mk7A+hBU/zA03BC+u/OiSWgDUEZHaIlIc6A1MARCRAUBHoI+qZmQ+QUTOEHE5VERaeXHuyIdYTSF1ySVu+OKDVS25PWoEMeuW8OCsd0lKSWXYxOWWNEyBkln6BiB25SwOhhXjky030zH2AGXLBi+uQKfVdheRTUBrYKqIxHvHq4nINABVPQwMAuKBVcBnqrrCe4k3gCrAT8dMn+0JJIjIMuAVoLcWphWGJihiY6FWtzWM/esuXq9wI7csiKPrylmkHkpnRHxOczeMyT+ZpW9qlC3OVavmMLPSZaRkVOD5J04LalwBlQZR1UnAJB/HNwOdszyeBkzz0c7n+6vqa8BrgcRmjC8ZdddR7iLh7h/G0KTMCp79+lUSK9ZkFWcFOzRjjhLbPJrYvxJg3y6+KHkr7dsHv2CBrfQ2RUq1qEjKXZhIyWab6fH3V+wMj+LNSU9Rv1hasEMz5njjxnHwtAq8v6MLgwYFOxhLGKaIGdKxHqWKh1PhigT21U2n+8EpnP73TsZ9/7KtBDcFS0oKxMUxtdx1VD2zBP/6V7ADsoRhipjMvuHqFSKpdNVS1tSuw92Mosr82fDII8EOz5h/fPYZHDjAU0k3cscdEB4e7ICsWq0p4lJS4OKL4d41t3LzoTFur+SePYMdljFw4YVsWrWH+mnL2bhJKF8+797K32q1th+GKdKiouCbb+DS1q/QdMuvtLjxJsIaNoSGDYlbksSI+DVsTkmlWlQkQzrWO7L/hjF56rff4KefeDVsBDffmbfJ4mRYl5Qp8qKj4ctvS3Bj6c/ZeaAUh7t256s5q46q5WPrNUy+eu89MiSMD/Xf3HNPsIP5hyUMY4D69WHstGj6hE+A9espM/AWDqQdOqqNrdcw+SI9nYz33uf78E607lGVswrQjG9LGMZ4WreGuz+/hAd4nrZr5nHnvAnHtTmqxo8xeeH77wlL2sTYwzdy3305N89PljCMyeKqq6DJmLv4kH9z37wPuWTdoqPOZ63xY0xe0DFj2RlWia2tunLhhcGO5miWMIw5Rv8BwvRbRvIrTXl50gvUSNkCQGSxcIZ0rBfk6EyhtnUrOnkyb2fcyN1DSgQ7muNYwjDGh3fePIOXr/wA0oXRH43grFLCMz2a2Cwpk6f03XGEpR9meq0BdO8e7GiOZwnDGB9E4K0vmzDqwo9ptPc33vn2Y2KbVQt2WKYwy8gg9ZWxzOYSrnm0XoFYqHcsSxjGZCM8HB6YfiXvnPkkteZ9xJq7rB6myUOzZlEqeR1flL+Fvn2DHYxvljCMOYGSJaHn4oeYVeYqzhp1H2venpvzk4w5BdueHstOylN32NUULx7saHyzhGFMDqIqhFHvlw/YGFGb8gN7sWHe5pyfZMzJ2L6dqJkTmVDyem66o+DOxLOEYYwfqtYvh0ycSKmMvey6vCdb/rRy6Cb3bP6/tymWkYYOGEjp0sGOJnuWMIzxU+2rGpP81Lu0OPgTc1rcw+7dwY7IFArp6USMfZ254ZfS+/8aBTuaEwo4YYhILxFZISIZIpJttUMR6SQia0QkUUSGZjk+TkQ2eFu0LhWRZt5xEZFXvPa/ikiLQGM1JlB1HurFhp5DuGbH67xx/rscOBDsiEyoW/fqVE5P/YPkqwcRFRXsaE4sN64wEoAewJzsGohIODAKuBJoCPQRkYZZmgxR1Wbebal37EqgjncbCLyeC7EaE7DanzxNcqPLGbzmdh7rvND2XTIB2fP0KJIkmo6juwU7lBwFnDBUdZWq5lSRrRWQqKrrVTUNGA/k9NvpBryvznwgSkSqBhqvMQGLiKDqrPGklT+DQTN78J+btlGItpUx+WjFxDU03/Ytv112G+UqFvzdJvJrDCMa2Jjl8SbvWKanvG6nF0WkhJ/PMSZ4KlWi7PcTqRqxjS4fXEvtyxKpPXQqbYbPsBLoxm/rHhhNGsVo+cYtwQ7FL34lDBH5XkQSfNz8vYYSH8cy/yYbBtQHzgMqAA/68ZyssQ0UkYUisnDbtm1+hmNMLmjRgmWPDacdM7lr9uvsXljL9s0wOYpbkkSb4TM457pvabthHAvqxFK2TpVgh+UXvxKGqrZX1cY+bpP9fJ9NQI0sj6sDm73XTva6nQ4C7+K6r074nGNiG6OqMaoaU7lyZT/DMSZ3DC7WhHdaXMX9vECn6cvYmxBt+2aYbMUtSTqyMde/4hdTjj28fF6bkPkDI7+6pBYAdUSktogUB3oDUwAyxyVERIBY3CA63vl+3mypC4DdqpqcT/Ea45fNKak83a4/86s35i0ZQPWpGexPPN32zTA+jYhfQ+qhdA6uK8+dO99mYdmm/FL97JD5AyM3ptV2F5FNQGtgqojEe8ericg0AFU9DAwC4oFVwGequsJ7iY9EZDmwHKgE/M87Pg1YDyQCY4E7Ao3VmNxWLSqSw+ERDIp9kJTSpYkLj0XjanJwY0Ub0zDH2ZySiipc9N0G6pDI25d0AZGQ+QMjN2ZJTVLV6qpaQlWrqGpH7/hmVe2cpd00Va2rqmer6lNZjrdT1SZeF1dfVd3rHVdVvdNr30RVFwYaqzG5bUjHekQWC2d76fLc2uNhzmALn8q1bJtwLge2lLUxDXOUalGR7F9dlbt3v8EfkdX4tsEFR46HAlvpbUwAYptH80yPJkRHRZJQtQ6PdLqLdodnM0L+w1+fteLQ9tNsTMMccU+7etSZ+TcX8wPjLuhCelh4SG3MVfAn/hpTwMU2jz6ysVKtodBoy1ruWfQaS4s348NPr+OMvj+RRGh0OZi8tfWXaO7++372FDuNCedeQXRUJEM61guZjbksYRiTi8JFeOqy/tTb9gdvJN3BqoMNWTz+fKr1nR/s0EyQ7dwJYx/5g8XyOWGD7yNhRK9gh3TSrEvKmFyUrsrh8Aju7PYg204rz6Ri3ai0dxfJ41uxfXswxJ4hAAAbzUlEQVSwozPB9OSTMGD384SFC3L3XcEO55RYwjAmF0V7g5e7SpVjYI+HiTq0h8lRnQnfFU7HjliF2yJq9Wr47LW/uDV8LNKvH9SokfOTCiBLGMbkosxZUwCrTj+LBzrfS6vty/nhvFtY/qvSuTPs3RvkIE2+u/9+uD/8JYplHIQHH8z5CQWUJQxjclHWWVMCLLvgClbfMpiW8z/ml+tfZf586NoVUm0MvMj45huYNy2FQWGjkF69oG7dYId0ymzQ25hclnXWFAAZl8K2P2j23r18O6QBVzx3Bd27w+TJUKJEti9jCoEDB2DQIHi84ihK7NjDDWe0Z87QqVQLsdlRmewKw5g8FrcsmfbNBrC6Qg3Oe7kHT9yxgPh4uOYaSLOdXgu14cMhed0+bj3wArPOacXsyGoohOyCTksYxuShzGJziakw4OpHOSRhXP9Jd+6553emTIHeveHQoWBHafLC2rXwzDPwZvM3KbVvJ69ccPQ02lBc0GkJw5g8lFlsDmBT1Bnc1v0holO2ct1XV/PKC4eZNAmuu86SRmGjCnfeCRVL7KXPn8P54cxmLI5ucFy7UKkhlckShjF56NgvhAU1GvNwxzs4L3Exd62/lxdegM8/h6ott1LrP9OsWGEh8dln8N138PmlrxG+Yxvvd+7vs12o1JDKZAnDmDzk6wthQtMOfHxxL3jtNdptfYrKl69mx/IqbPuyGZt2HODeT5dSyyrdhqxt2+Cuu6Bts920/uE56NKFzrd0PzLdOlMo1ZDKZAnDmFyQuYvasSXNs67LyBRZLJxSL4yEq66i8XOPcWWFzyh/2Ur2r67G9i+bk5HuNpsM1YHRom7wYEhJgU8vfAnZtQuefPK46dbRUZE806NJyM2Ssmm1xgQoc2A7c6wi84seOPKFMCJ+DZtTUo+eTvnxx6w5uymvTR5Oct+KLBTYNaMh2xUqdV2ChOuRgdFQ+2IpqiZPhk8+geeG7qTK6Bege3do0QLwMd06BInqcdtkh6yYmBhduNC2zTD5q83wGST5GLyMjopk3tB2J3xu96HjeWPUnRwKiyC23/OsX9WcXdMbEXnOVip3W4xEZCDAhuFd8ih6k1t27YKGDaFKFVjUYRjhI5+FZcugSZNgh5YjEVmkqjE5tbMuKWMClN1MF39mwNxw7cXcec3jVNy/m7e+eJLTz11NhQ7LSU2swl9fxJBxKCzkBkaLqrvucuMXHz6zkfBXX4J//zskksXJCChhiEgvEVkhIhkikm12EpFOIrJGRBJFZGiW43NFZKl32ywicd7xS0Vkd5ZzjwUSpzF5KbsvdH++6GObR9N30NU80edhmiYn8vKXIyl37gYqXrmMA79XYvvn5zPoovq5HbLJZZ98Ah99BI8+Co3HP+Lm1f7vfzk/McQEeoWRAPQA5mTXQETCgVHAlUBDoI+INARQ1YtVtZmqNgN+AiZmeerczHOq+mSAcRqTZ7Ib2PZ3Bkxs82ieff9Rwl5+iY5r5/P0j+9Tpukm6ly7koNJ5Xnp/mrs2JEXkZvc8McfcPvt0Lo1PNxlKXzwgRv5PvPMYIeW6wJKGKq6SlVzWqrYCkhU1fWqmgaMB7plbSAiZYB2QFwg8RgTDLk2A+buu2HwYHrP+4INVRP5bXwjJk0Uli2Dtm1h8+Y8Cd8EID0drr8eMjLgww+UiKEPQIUKMGxYsEPLE/kxSyoa2Jjl8Sbg/GPadAemq+qeLMdai8gyYDPwgKqu8PXiIjIQGAhQs2bNXAvamJORazNgnn8e/vwT7r0Xqlen69VX8/XXrsLtxRe7xWBnnRX425jcMXw4zJ0L770HZ/32DUyfDi+9BFFRwQ4tT+R4hSEi34tIgo9bt5yem/kSPo4dOzWrD/BJlseLgTNV9VzgVU5w5aGqY1Q1RlVjKleu7GdIxhRQ4eGuM7x1azdoOncul13mvodSUqBNG1i6NNhBGoCZM+Gxx1w9sOuvTYP77oNzznH9U4VUjglDVduramMft8l+vscmIOv2UtVxVw0AiEhFXLfV1CzvuUdV93r3pwHFRKSSn+9nTGiLjIQpU6BWLejWDVaupFUr+OEHKFbMdU/NmhXsIIu2zZtdoqhbF8aMAXn1Fbet3ssvQ/HiwQ4vz+THtNoFQB0RqS0ixYHewJQs53sBX6nqgcwDInKGiIh3v5UXpw37maKjYkX4+mu3YUbHjrBxIw0awI8/QvXq7tCECcEOsmg6dMgli717XR2wMn9vhv/+F666Cjp3DnZ4eSrQabXdRWQT0BqYKiLx3vFqIjINQFUPA4OAeGAV8Nkx4xG9Obo7CqAnkOCNYbwC9NbCtMLQGH/Uru22a9uzBzp0gB07qF7d9Zmfd57bT+P5590MTpN/hg1z/wZjxkCjRsCQIS6LvPhisEPLe6paaG4tW7ZUYwqd2bNVS5RQPf981b17VVU1NVW1Vy9VUB00SPXw4SDHWES8+677nd95p3dg9mx34NFHgxlWwICF6sd3rK30Nqagu+QSGD8eFixwtYkOHqRkSXfogQfgtdcgNhb+/jvYgRZuc+fCwIHQvr13MXHwINx2G9SsCUOH5vj8wsAShjGhIDYW3nrLzau97jo4fJiwMBgxAkaNcsMdbdq4RWQm961f73J17dpur4tixYCnn4ZVq+D116FUqWCHmC8sYRgTKm66yc3xnzgRBgxwq8WAO+5wCePPP6FVK/jppyDHWcjs2AH/+pf7dX/1FZQvDyQkuP1Xr7uu0A90Z2UJw5hQMniwm5Hz3ntuZbg34n3FFTB/PpQp46bdvvVWkOMsJPbtgy5d3BXGpElQpw5uefeAAVC2rEvgRYjth2FMqHn0UTdgMXKkm3Y7ciSIUL8+/PIL9OkDt9wCixYV+mUBeSotDa6+2g0dffGFS8QAvPoq/PwzfPghFLHFwpYwjAk1IvDcc27Q9YUXXNJ46ikQoUIFmDYNHnrINVm2zPW5V68e7KBDS3o63HgjxMe7q7XYWO/EqlVuXm3nzq47qoixhGFMKBJxlw8HD7q+9GLFiOt2CyO+/c3t7Fc+kgeGN+ON/1WgeXNXbaRDh2AHHRoOH4YbbnAly4cPh/79vRNpaa5cy2mnwdtvu3+DIsbGMIwJVSJuhs7NN8OTT5I8+D8k7dqP4raJ/WrvL/Qdvob94Xvp2FGpcfkGvlhg+4OfyOHD0K8ffPyxy8MPPpjl5BNPwJIlMHYsnHFGsEIMKksYxoSysDAYO5Yp53Xm9h8+Ycic948MhKceSuebjYlU/PdcSjfexKYZtbm+RynGTN0S5KALpswLiE8+gWefPWZpxdy57nJjwIAs/VNFj3VJGRPqwsK457Lb+PtQBnfOn0BERjrPXHoTiKBAWLEMKnX5lcja29kR35g7ep5G2XddPSTj7NkDPXq4qsAjR8L992c5uWMH9O3r6soXhfIfJ2AJw5hCoGr50jzS4Q4OhUVw6y8TKZ2WyqMdbkfln06E0g03U7zaLnZ82Zw+fcozaZJb9FepiNeBTk52Y9gJCW62cr9+WU5mZLhksWULzJvnxi+KMOuSMqYQGNKxHiWLF+OJ9rcy+oKe9F36NS9MfZHwjPSj2hWLSqXFHUt5+mm3rqBxY5js70YFhdDSpW7rkbVr3aK8o5IFuH25v/kGXnkFYmKCEmNBYgnDmELgyDax5Usxou2NvNGxP91XzGTM5GcocejgkXaRxcL5T+e6DBsGCxdC1aquS75nz6K3BezHH8OFF7qB7tmzXcn4o3z7rRvovv56V0TKWLVaYwqt117TDBFdWquxNh08Xi98ZrpOWrzpqCZpaapPP61asqRq2bKqo0cX/sq3Bw+q3nuvKzJ7ySWqW7b4aLR2rWqFCqqNG6vu25fvMeY3/KxWG/Qv+dy8WcIw5hiffaZavLhqo0aqGzdm22ztWtXLL3ffCM2auardhdGqVaotWuiRsvBpaT4a7dihWreuasWKqomJ+R5jMPibMKxLypjCrFcv1wf/559wwQXZbgh+zjmuEO6nn7pJQW3bwrXXQmJiPsebRzIyYPRoaNHCVfSdONFV+ChW7JiGaWmuf27DBjfIc/bZQYm3oAo4YYhILxFZISIZIpLtqJCIvCMif4lIwjHHK4jIdyKy1vtZ3jsuIvKKiCSKyK8i0iLQWI0pki67zG0ILgIXXQRTp/psJuJ28Vu9Gh5/HL78Eho0cNVwk5PzOeZc9OuvbkuRO+90P5cvd6XKj6Pq9reYOdPVA7n44nyPtaDLjSuMBKAHMCeHduOATj6ODwWmq2odYLr3GOBKoI53Gwi8nguxGlM0NW3qKhPWrw9du7pZP+p7b9dSpdxY77p1bqx37Fi3BOHuu2HjxvwNOxApKW49RYsWLgm+8w7c+nQSPd+bQe2hU2kzfAZxS7yV76puWfe778Ijj/iYLmUgFxKGqq5S1TV+tJsD7PRxqhvwnnf/PSA2y/H3vS62+UCUiFQNNF5jiqyqVd10oK5dXZn0/v3hwIETNh81yn3Z9u7tqpCcfbarRPLrr/kY90nat88tyq5d29VmvPlmWLMGyjdL4qFJy0lKST1SPmXYxOUuaTz9NIwYwfpe/WhT6tLjE4oBCsa02iqqmgzg/TzdOx4NZP17ZpN3zBhzqkqXdrW6H3/c/TXdti0knfhL8eyzXdN161yPzfjxcO65rnvn009d/cOCYPt2971/9tmuoGybNq7005gxULEijIhfQ+qho9elpB5KZ8Pjw+GRR/izy9X8q+61JO0+cHxCMYCfCUNEvheRBB+3bnkYm69SkMddQ4vIQBFZKCILt23blofhGFNIhIW5PqdJk2DlStdn8/33OT6tZk3Xk7VpkyufkZTkrjyqVnWJZN68I5sA5htVtzXFLbdAjRrw8MMumc2b5xbiNWv2T9vNKanHPf/mBZO590u3KXrf1rey//DRXzGph9IZEZ9jB0qR4VfCUNX2qtrYxy031ohuzexq8n7+5R3fBNTI0q46cNzSIlUdo6oxqhpTuYhtZmJMQGJj3bhGpUqu9vnjj7uNIHJQoYIbG1i71k3A6twZPvjAjadHR7v6fJMnw+7deRN2errb1OjRR93srgsucHsZ9evnynvEx7sFeceqFhX5zwNV7p73CY/NGMuMxu5S6c+/03y+n69EU1QVhFpSU4AbgOHez8lZjg8SkfHA+cDuzK4rY0wuadDAJY1Bg+DJJ90Yx/vvu8uJE4hbksSI+DVu743Gkbx1W30y/qjGlCkwYYLbLiIszP21f9FFrqpGkybu7UqWPLkQd+1ys4GXLnVXDjNmuGNhYXD55S5xdO8O5cqd+HWGdKzHsInLOZB2iIdmvsMtC+KY2PQKwt8aS9yKbQg+ujA4JtEUcaLZzJTw+wVEugOvApWBFGCpqnYUkWrAW6ra2Wv3CXApUAnYCjyuqm+LSEXgM6Am8CfQS1V3iogAr+FmVu0HblLVhSeKJSYmRhcuPGETY0x2xo1ziSMiAl57zdX69rFJUNySJIZNXH7UeEBksXCe6dGE2ObRpKW5L/bZs2HOHLfXeKr3R3pYGFSr5q5EoqOhfHmIjHQzs1Rdu/373XjEpk1uVlbWnuYaNaB9e7eHefv2J79D6pc/ruW0W/tzWcJcJrTuTrFXXya2ZQ3aDJ9Bko8rCQFevLYZsc0L9/CpiCxS1RyLZQWcMAoSSxjGBGjdOrfd3Lx5bgHba69BlSpHNcnuyzU6KpJ5Q9sdd/zwYdd9lZDgbn/84cY/kpLc1NfUVHcTcYmjVCmIinLJoUYNN6W3eXM3HnH66ce9vP82boRu3dy+tSNHwj33HEmItYdO9Xl1AfD78C4BvGlo8DdhFIQuKWNMEBzVrRQVyZCO9Yhtfra7NBgxwo1pTJ/u7t9885Ev1+z69LM7HhHhuqIaNHALz4Ni+nR3xZSa6kbDr7zyqNPVoiKzTYLmHwVhWq0xJp9ldiv5XJMQHu62m1u2zNU/HzDArRZfvhzIvk+/QPb1p6W5z3LFFa7/66efjksW4MY3IouFH3Usslg4QzrWy69IQ4IlDGOKoOzWJBw1hbR+fZg1yy1kWL7c9QndfjsPt6oUGl+uCQluxP3ZZ92820WLoGFDn02PlIePikRwVxaZYzLmHzaGYUwRlF2fvQAbfPXZ79zpZlGNGgWlSrGqzwAGn9GWtQfCs3RnFZAv13374L//dduplivnEl6PHsGOqkDzdwzDrjCMKYJOulupQgV46SV3pdGuHQ3efIFvX76BDaUXM29gAZlFlJ7upgQ3bOjGXfr1c3VNLFnkGksYxhRBp9xnX7++WyG+eLGrDfLYY24q0x13uC/nYEhPd/VKGjVyM7wqVIC5c91ikKK+YXkus4RhTBEUcJ998+ZuOffSpW7jjLffdtOg2rZ1pcF9LPOOW5JEm+E+KsWeqq1b/yke1aePm471+edurOKiiwJ7beOTjWEYYwK3datLGu+/70rDlizplmF36QKdOxO3M+KEi/38tmWL26gjLs7t+HToELRr565wYmPdDC9z0mzhnjEm/6m6Qk8ffeTWO6xfD0ByVBUWnVGXJVXr8lvlM/kz6gw2l63M6RXL+lzshyr89Rf8/rub7fTzz27JeEKCO1e7Nlx9tZvyW6+Azc4KQZYwjDHBpQq//QZff83UN77g3OQ1VN/zT52PdAkjpeRplKwQRemKUW5hYOay7+3b/6knAm7p9/nnu66mbt3c+hAfZUvMqbGV3saY4BJxf/3Xq8fTB5qSlJJKpX27qL0ziTNTtlBzVzLlD/xNuUMHOK9CMaqWLeEKS0VGuoHr2rWhVi2oU8fdwmzINdgsYRhj8lxmpdjtpcuzvXR5FtRofNT57OpQmYLFEoYxJs9lDmzf8+lSn+dtz4nQYNd4xph8Eds8OttifgWyDpU5jiUMY0y+sSJ/oc26pIwx+Saza+r4suoFoLSIyZElDGNMvoptHm0JIkQF1CUlIr1EZIWIZIhItnN4ReQdEflLRBKOOT5CRFaLyK8iMklEorzjtUQkVUSWerc3AonTGGNM4AIdw0gAegBzcmg3Drc397G+AxqralPgN2BYlnPrVLWZd7stwDiNMcYEKKCEoaqrVHWNH+3mADt9HP9WVQ97D+cD1QOJxxhjTN4pSLOkbga+zvK4togsEZHZInJxsIIyxhjj5DjoLSLfA2f4OPWwqk7OjSBE5GHgMPCRdygZqKmqO0SkJRAnIo1UdY+P5w4EBgLUrFkzN8IxxuSCuCVJNhuqkMkxYahq+7wMQERuAP4FXK5eJURVPQgc9O4vEpF1QF3guMqCqjoGGAOu+GBexmqM8U/ckqSjypknpaQybOJyAEsaISyoXVIi0gl4EOiqqvuzHK8sIuHe/bOAOsD64ERpjDlZI+LXHLX3BUDqoXRGxOc45GkKsECn1XYXkU1Aa2CqiMR7x6uJyLQs7T4BfgLqicgmEenvnXoNKAN8d8z02UuAX0VkGfA5cJuqHjdobowpmLKrDWU1o0JbQAv3VHUSMMnH8c1A5yyP+2Tz/HOyOf4F8EUgsRljgqdaVCRJPpKD1YwKbQVplpQxppCwmlGFk5UGMcbkOqsZVThZwjDG5AmrGVX4WJeUMcYYv1jCMMYY4xdLGMYYY/xiCcMYY4xfLGEYY4zxiyUMY4wxfrGEYYwxxi+WMIwxxvjFEoYxxhi/WMIwxhjjF0sYxhhj/GIJwxhjjF8sYRhjjPGLJQxjjDF+CXSL1l4iskJEMkQk5gTt3hGRv0Qk4ZjjT4hIkrc961IR6Zzl3DARSRSRNSLSMZA4jTHGBC7QK4wEoAcwJ4d244BO2Zx7UVWbebdpACLSEOgNNPKeN1pEwrN5vjHGmHwQUMJQ1VWqusaPdnOAnSfx0t2A8ap6UFU3AIlAq1MM0xhjTC4oCGMYg0TkV6/bqrx3LBrYmKXNJu/YcURkoIgsFJGF27Zty+tYjTGmyMoxYYjI9yKS4OPWLRfe/3XgbKAZkAw8n/m2PtqqrxdQ1TGqGqOqMZUrV86FkIwxxviS457eqto+r95cVbdm3heRscBX3sNNQI0sTasDm3N6vUWLFm0XkT8CCKkSsD2A54eaovZ5wT5zUWGf+eSc6U+jHBNGXhKRqqqa7D3sjhtEB5gCfCwiLwDVgDrALzm9nqoGdIkhIgtVNdvZXoVNUfu8YJ+5qLDPnDcCnVbbXUQ2Aa2BqSIS7x2vJiLTsrT7BPgJqCcim0Skv3fqORFZLiK/ApcB9wKo6grgM2Al8A1wp6qmBxKrMcaYwIiqz6GBIqmo/VVS1D4v2GcuKuwz542CMEuqIBkT7ADyWVH7vGCfuaiwz5wH7ArDGGOMX+wKwxhjjF+KXMIQkU5efapEERnq43wJEfnUO/+ziNTK/yhzlx+f+T4RWektoJwuIn5NsSvIcvrMWdr1FBE9US20UOHPZxaRa7x/6xUi8nF+x5jb/Phvu6aIzBSRJd5/3519vU6oyK4uX5bzIiKveL+PX0WkRa4GoKpF5gaEA+uAs4DiwDKg4TFt7gDe8O73Bj4Ndtz58JkvA0p5928vCp/Za1cGVwdtPhAT7Ljz4d+5DrAEKO89Pj3YcefDZx4D3O7dbwj8Huy4A/zMlwAtgIRszncGvsYtfr4A+Dk337+oXWG0AhJVdb2qpgHjcXWrsuoGvOfd/xy4XER8rTwPFTl+ZlWdqar7vYfzcQslQ5k//84A/wc8BxzIz+DyiD+f+RZglKruAlDVv/I5xtzmz2dWoKx3vxx+LAAuyDTnunzdgPfVmQ9EiUjV3Hr/opYw/KlRdaSNqh4GdgMV8yW6vOF3XS5Pf9xfKKEsx88sIs2BGqr6FYWDP//OdYG6IjJPROaLSHYVpEOFP5/5CaCvt15sGnBX/oQWNCf7//tJCepK7yDwp0aV33WsQoTfn0dE+gIxQNs8jSjvnfAzi0gY8CJwY34FlA/8+XeOwHVLXYq7ipwrIo1VNSWPY8sr/nzmPsA4VX1eRFoDH3ifOSPvwwuKPP3+KmpXGP7UqDrSRkQicJexJ1OavaDxqy6XiLQHHga6qurBfIotr+T0mcsAjYFZIvI7rq93SogPfPv73/ZkVT2kbtuANbgEEqr8+cz9cVUjUNWfgJK4mkuF1SnV4fNXUUsYC4A6IlJbRIrjBrWnHNNmCnCDd78nMEO90aQQleNn9rpn3sQli1Dv14YcPrOq7lbVSqpaS1Vr4cZtuqrqwuCEmyv8+W87DjfBARGphOuiWp+vUeYufz7zn8DlACLSAJcwCvM+CFOAft5sqQuA3fpPvb6AFakuKVU9LCKDgHjcDIt3VHWFiDwJLFTVKcDbuMvWRNyVRe/gRRw4Pz/zCOA0YII3vv+nqnYNWtAB8vMzFyp+fuZ4oIOIrATSgSGquiN4UQfGz898PzBWRO7Fdc3cGMp/AHp1+S4FKnnjMo8DxQBU9Q3cOE1n3KZz+4GbcvX9Q/h3Z4wxJh8VtS4pY4wxp8gShjHGGL9YwjDGGOMXSxjGGGP8YgnDGGOMXyxhGGOM8YslDGOMMX6xhGGMMcYv/w+rHqTqriWkVgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7d901724e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Prepare the data\n",
    "N = 30\n",
    "sigma = 0.005\n",
    "x = numpy.sort(numpy.random.uniform(0.0, 1.0, N))\n",
    "def f(t):\n",
    "    return 1.5*t**3 - 2*t**2 + 0.5*t - 1\n",
    "y = f(x) + numpy.random.normal(0.0, sigma, N)\n",
    "# Append an outlier\n",
    "x = numpy.append(x, 0.73)\n",
    "y = numpy.append(y, -0.99)\n",
    "X = [[t**i for i in range(4)] for t in x]\n",
    "\n",
    "# Solve the models\n",
    "LSE, HUB = lse(X, y), huber(X, y, 5*sigma)\n",
    "LSE.solve()\n",
    "HUB.solve()\n",
    "wLSE = LSE.getVariable(\"w\").level()\n",
    "wHUB = HUB.getVariable(\"w\").level()\n",
    "\n",
    "# Plot the data points, LSE approximation (blue) and Huber approximation (red)\n",
    "plt.scatter(x,y)\n",
    "ticks = numpy.linspace(0, 1, 100)\n",
    "plt.plot(ticks, sum([wLSE[i]*ticks**i for i in range(4)]), 'k', color='blue')\n",
    "plt.plot(ticks, sum([wHUB[i]*ticks**i for i in range(4)]), 'k', color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The red curve is clearly less affected by the outlier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/80x15.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">Creative Commons Attribution 4.0 International License</a>. The **MOSEK** logo and name are trademarks of <a href=\"http://mosek.com\">Mosek ApS</a>. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API` are not guaranteed. For more information contact our [support](mailto:support@mosek.com). "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}