{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Fitting censored data\n", "\n", "Experimental measurements are sometimes censored such that we only know partial information about a particular data point. For example, in measuring the lifespan of mice, a portion of them might live through the duration of the study, in which case we only know the lower bound.\n", "\n", "One of the ways we can deal with this is to use Maximum Likelihood Estimation ([MLE](http://en.wikipedia.org/wiki/Maximum_likelihood)). However, censoring often make analytical solutions difficult even for well known distributions.\n", "\n", "We can overcome this challenge by converting the MLE into a convex optimization problem and solving it using [CVXPY](http://www.cvxpy.org/en/latest/).\n", "\n", "This example is adapted from a homework problem from Boyd's [CVX 101: Convex Optimization Course](https://class.stanford.edu/courses/Engineering/CVX101/Winter2014/info).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "\n", "We will use similar notation here. Suppose we have a linear model:\n", "\n", "$$ y^{(i)} = c^Tx^{(i)} +\\epsilon^{(i)} $$\n", "\n", "where $y^{(i)} \\in \\mathbf{R}$, $c \\in \\mathbf{R}^n$, $k^{(i)} \\in \\mathbf{R}^n$, and $\\epsilon^{(i)}$ is the error and has a normal distribution $N(0, \\sigma^2)$ for $ i = 1,\\ldots,K$.\n", "\n", "Then the MLE estimator $c$ is the vector that minimizes the sum of squares of the errors $\\epsilon^{(i)}$, namely:\n", "\n", "$$\n", "\\begin{array}{ll}\n", " \\underset{c}{\\mbox{minimize}} & \\sum_{i=1}^K (y^{(i)} - c^T x^{(i)})^2\n", "\\end{array}\n", "$$\n", "\n", "In the case of right censored data, only $M$ observations are fully observed and all that is known for the remaining observations is that $y^{(i)} \\geq D$ for $i=\\mbox{M+1},\\ldots,K$ and some constant $D$.\n", "\n", "Now let's see how this would work in practice.\n", "\n", "\n", "## Data Generation" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "n = 30 # number of variables\n", "M = 50 # number of censored observations\n", "K = 200 # total number of observations\n", "\n", "np.random.seed(n*M*K)\n", "X = np.random.randn(K*n).reshape(K, n)\n", "c_true = np.random.rand(n)\n", "\n", "# generating the y variable\n", "y = X.dot(c_true) + .3*np.sqrt(n)*np.random.randn(K)\n", "\n", "# ordering them based on y\n", "order = np.argsort(y)\n", "y_ordered = y[order]\n", "X_ordered = X[order,:]\n", "\n", "#finding boundary\n", "D = (y_ordered[M-1] + y_ordered[M])/2. \n", "\n", "# applying censoring\n", "y_censored = np.concatenate((y_ordered[:M], np.ones(K-M)*D))\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "# Show plot inline in ipython.\n", "%matplotlib inline\n", "\n", "def plot_fit(fit, fit_label):\n", " plt.figure(figsize=(10,6))\n", " plt.grid()\n", " plt.plot(y_censored, 'bo', label = 'censored data')\n", " plt.plot(y_ordered, 'co', label = 'uncensored data')\n", " plt.plot(fit, 'ro', label=fit_label)\n", " plt.ylabel('y')\n", " plt.legend(loc=0)\n", " plt.xlabel('observations');\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Regular OLS\n", "\n", "Let's see what the OLS result looks like. We'll use the `np.linalg.lstsq` function to solve for our coefficients." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAF3CAYAAAD3rnzeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3Xt8VOWdP/DPMxMcDCCVUMYLTUL8ISJyK0HlJ2hQtFZbbK1WKLBYF7PieumF3W2brna7xq2V1upapZFta01W1FZr7ba/1kgGULFtqGiBFCpXqTJI2AohEEny/f0xM2EyOWfmnJlznfN5v17zyuTMmTPPc27zneeqRARERERE5G0htxNARERERLkxaCMiIiLyAQZtRERERD7AoI2IiIjIBxi0EREREfkAgzYiIiIiH2DQRkREROQDDNqIiIiIfIBBGxEREZEPMGgjIiIi8oEStxNgh5EjR0plZaWtn3HkyBEMGTLE1s/wsiDnP8h5B5h/5j+4+Q9y3gHm3878b9iw4YCIfDjXekUZtFVWVqK1tdXWz4jFYqipqbH1M7wsyPkPct4B5p/5D27+g5x3gPm3M/9Kqd1G1mP1KBEREZEPMGgjIiIi8gEGbUREREQ+UJRt2rQcP34ce/fuxbFjxyzZ3vDhw9HW1mbJtvzIS/kfPHgwRo8ejUGDBrmdFCIiItsEJmjbu3cvhg0bhsrKSiilCt7e4cOHMWzYMAtS5k9eyb+IoL29HXv37sWYMWPcTg4REZFtAlM9euzYMZSVlVkSsJF3KKVQVlZmWQkqERGRVwUmaAPAgK1I8bgSEVEQBCpoI+OGDh2ac51vfOMbWL58edZ1fv7zn2PLli1WJYuIiCiwGLTpaGoCKiuBUCjxt6nJ7RQVrqenx/HPZNBGRERkDQZtGpqagNpaYPduQCTxt7a28MDtJz/5CSZNmoTJkydj0aJFAID33nsPn/nMZzB9+nRMnz4dr7zyCoBEKdZNN92EmpoaVFVV4aGHHgKQmEbj6quvxuTJk3HeeefhqaeeAgC89NJLmDp1KiZOnIibbroJXV1dABKzQ3zzm9/EzJkz8cwzz2D79u248sorMW3aNMyaNQt//vOfAQA7d+7EjBkzMH36dPzrv/6rbh7q6+sxbtw4zJ07F1u3bu1b/thjj2H69OmYPHkyPvOZz6CzsxOvvvoqfvGLX+Cf/umfMGXKFGzfvl1zPSIiIi9pisdRuX49QrEYKtevR1M87naSEkSk6B7Tpk2TTFu2bBmwTE9FhUgiXOv/qKg4sc6hQ4cMb09EZNOmTXL22WfLe++9JyIi7e3tIiIyf/58WbdunYiI7N69W8455xwREbn77rtlxowZcuzYMXnvvfdkxIgR8sEHH8hPf/pTWbJkSd92//a3v8nRo0dl9OjRsnXrVhERWbRokTzwwAPJvFTIfffd17f+pZdeKtu2bRMRkddee01mz54tIiKf/OQn5fHHHxcRkYcffliGDBkyIA+tra1y3nnnyZEjR2Tv3r1y1llnyf333y8iIgcOHOhbr66uTh566CEREVm8eLE888wzfa/prVcoM8e3UC0tLY59lhcx/y1uJ8FVQc5/kPMuEpz8N+7bJ6Vr1ghaWvoepWvWSJ2N+QfQKgbiG5a0adizx9xyI1avXo3rrrsOI0eOBACMGDECANDc3IzbbrsNU6ZMwdy5c3Ho0CEcPnwYAHD11VcjEolg5MiRGDVqFOLxOCZOnIjm5mb8y7/8C9atW4fhw4dj69atGDNmDM4++2wAwOLFi7F27dq+z77hhhsAAB0dHXj11Vdx/fXXY8qUKfiHf/gHvPvuuwCAV155BfPnzweAvlLATOvWrcOnP/1plJaW4pRTTsHcuXP7Xtu0aRNmzZqFiRMnoqmpCZs3b9bchtH1iIiI3FC3Ywc6e3v7Levs7cVKl9KTLjDjtJlRXp6oEtVani8R0ezl2Nvbi/Xr1+Pkk08e8FokEul7Hg6H0d3djbPPPhsbNmzAr371K3z1q1/FFVdc0S940jJkyJC+z/rQhz6EjRs3aq5npBem3jo33ngjfv7zn2Py5Mn48Y9/jFgsVtB6REREbtiTbF6Uab/D6dDCkjYN9fVAaWn/ZaWlieX5uuyyy/D000+jvb0dAHDw4EEAwBVXXIGHH364bz29gCrlnXfeQWlpKRYuXIhly5bhj3/8I8455xzs2rULb731FgDgiSeewCWXXDLgvaeccgrGjBmDZ555BkAikHzjjTcAABdddBFWrVoFAGjSabx38cUX47nnnsPRo0dx+PBhvPDCC32vHT58GKeffjqOHz/e7/3Dhg3rKznMth4REZEXlKcVmKQb5XA6tDBo07BgAdDQAFRUAEol/jY0JJbna8KECairq8Mll1yCyZMn40tf+hIA4KGHHkJraysmTZqEc889FytWrMi6nT/96U84//zzMWXKFNTX1+PrX/86Bg8ejB/96Ee4/vrrMXHiRIRCIdxyyy2a729qasJ//dd/YfLkyZgwYQKef/55AMCDDz6I73//+5g+fTref/99zfd+9KMfxQ033IApU6Zg4cKFmDVrVt9r//7v/44LLrgAl19+Oc4555y+5fPmzcP999+PqVOnYvv27brrEREReUF9VRVKQ/3Do9JQCEtcSk86lWj/Vlyqq6ultbW137K2tjaMHz/ess/wyjRObvFa/q0+vtnEYjHU1NQ48llexPwz/0HNf5DzDgQr/03xOOp27MCeri6URyKor6rCmW1ttuVfKbVBRKpzrcc2bURERERpFkSjWBCN9lsWa2tzKTUnsHqUiIiIyAcYtBERERH5gCNBm1Lqh0qp/UqpTWnLRiilXlRK/SX591Sd9y5OrvMXpdRiJ9JLRERE5DVOlbT9GMCVGcu+AuAlERkL4KXk//0opUYAuBvABQDOB3C3XnBHREREVMwcCdpEZC2AgxmLrwHwePL54wA+pfHWjwF4UUQOisj/AngRA4M/IiIioqLnZpu2qIi8CwDJv1rj1p0J4O20//cml5HHxGIxfOITn8i5Xk1NDTKHY8n0ve99jxPJExERZfD6kB9acyZpDiynlKoFUAsA0Wh0wPRIw4cP7zcyfy5Pt7fj3955B3s/+ACjTzoJd59xBj5bVtb3ek9Pj6nt+VHfBLWhgbF9Zv47OzvR3d2dc5/09PTgyJEjWdd74IEH8KlPfQplafs7l2PHjjk2JVZHR0egp99i/pn/oOY/yHkHmH9P5N/IrPJWPABUAtiU9v9WAKcnn58OYKvGe+YD+EHa/z8AMD/XZ02bNk0ybdmyZcAyPY379knpmjWClpa+R+maNdK4b1/fOocOHTK8PRGRnTt3yoQJE/r+v//+++Xuu+8WEZFLLrlE/vmf/1mmT58uY8eOlbVr14qISHd3t3z5y1+W8847TyZOnCgPPfSQiIi0trbKxRdfLB/96EfliiuukHfeeSfrdjZt2iTTp0+XyZMny8SJE2Xbtm0iIvKd73xHJkyYIBMmTJAHHnigL53nnHOOLF26VKZMmSK7du2S3/zmN3LhhRfK1KlT5brrrpPDhw/LoUOH5Ne//rWMGzdOLrroIrn99tvl6quvHpDvzs5OueGGG2TixIny2c9+Vs4//3z5wx/+ICIit9xyi0ybNk3OPfdcueuuu0RE5MEHH5RBgwbJeeedJzU1NbrrZTJzfAvV0tLi2Gd5EfPf4nYSXBXk/Ac57yLMv535B9AqBmIpN6tHfwEg1Rt0MYDnNdb5DYArlFKnJjsgXJFcZqu6HTvQ2dvbb1lnby/qduyw7TO7u7vx+9//Ht/73vfwb//2bwCAhoYG7Ny5E6+//jrefPNNLFiwAMePH8ftt9+On/70p9iwYQNuuukm1NXVZd3OihUrcOedd2Ljxo1obW3F6NGjsWHDBvzoRz/C7373O7z22mt47LHH8PrrrwMAtm7dir/7u7/D66+/jiFDhuCee+5Bc3Mz/vjHP6K6uhrf/e53cezYMdx888144YUXsG7dOuzbt08zX48++ihKS0vx5ptvoq6uDhs2bOh7rb6+Hq2trXjzzTexZs0avPnmm7jjjjtwxhlnoKWlBS0tLbrrERERBY0j1aNKqScB1AAYqZTai0SP0G8BeFop9fcA9gC4PrluNYBbRGSJiBxUSv07gD8kN/VNEcns0GC5PV1dppZb4dprrwUATJs2Dbt27QIANDc345ZbbkFJSeIwjRgxAps2bcKmTZtw+eWXA0hUN55++ulZtzNjxgzU19dj7969uPbaazF27Fi8/PLL+PSnP40hQ4b0vW/dunWYO3cuKioqcOGFFwIAXnvtNWzZsgUXXXQRAOCDDz7AjBkzsG3bNowZMwZjx44FACxcuBANDQ0D8rV27VrccccdAIBJkyZh0qRJfa89/fTTaGhoQHd3N959911s2bKl3+tm1yMiIipmjgRtIjJf56XLNNZtBU7MyyoiPwTwQ5uSpqk8EsFujQCtPBLJe5slJSXoTSu9O3bsWL/XI8lth8NhdHd3A0hUXSvVv1mfiGDChAlYv3695udobedzn/scLrjgAvzP//wPPvaxj2HlypWp6mZNqUAu9XmXX345nnzyyX7rvPLKKwPSpkdrvZ07d2L58uX4wx/+gFNPPRU33njjgH1iZj0iIqJixxkRNNRXVaE0o/F9aSiE+qqqvLcZjUaxf/9+tLe3o6urC7/85S9zvueKK67AihUr+oKvgwcPYty4cXjvvff6grbjx49j8+bNWbezY8cOVFVV4Y477sDcuXPx5ptv4uKLL8bPf/5zdHZ24siRI3juuecwa9asAe+98MIL8corr+Ctt94CkOhwsG3bNpx99tnYuXMntm/fDgADgrqUiy++GE1NTQCATZs29VVtHjp0CEOGDMHw4cMRj8fx61//uu89w4YN6+uokG09IiKiIPF671FXpCaJrduxA3u6ulAeiaC+qmrA5LFmDBo0CHfddRcuuOACjBkzBuecc07O9yxZsgTbtm3DpEmTMGjQINx888247bbb8NOf/hR33HEH3n//fXR3d+MLX/gCJkyYoLudp556Co2NjRg0aBBOO+003HXXXRgxYgRuvPFGnH/++X2fNXXq1L4q1ZQPf/jD+PGPf4z58+ejK1n6eM8992D27NloaGjA1VdfjZEjR2LmzJnYtGlT5kdj6dKl+PznP49JkyZhypQpfZ83efJkTJ06FRMmTEBVVVVf9SsA1NbW4uMf/zhOP/10tLS06K5HREQUJCpbNZlfVVdXS+ZYYG1tbRg/frxln3H48GEMGzbMsu35jdfyb/XxzSYWi6GmpsaRz/Ii5p/5D2r+g5x3IBj5b4rHdQts7My/UmqDiFTnWo8lbURERBR4TfE4ardu7Rs9YndXF2q3bgWAgmrarMQ2bURERBR4bgz3ZRaDNiIiIgo8N4b7MotBGxEREQWe3rBehQz3ZTUGbURERBR4dgz3ZTUGbURERBR4C6JRNIwbh4pIBApARSSChnHjPNMJAWDQ5qi9e/fimmuuwdixY3HWWWfhzjvvxAcffAAg0ZX4E5/4xID3/PKXv8TUqVMxefJknHvuufjBD34wYJ2uri7MmTMHU6ZMwVNPPYUlS5Zgy5YtAIB7773X3kwREREViQXRKHbNmIHemhrsmjHDUwEbwKBNX1MTUFkJhEKJv8lR/fMlIrj22mvxqU99Cn/5y1+wbds2dHR09JvsPdPx48dRW1uLF154AW+88QZef/11zTFiXn/9dRw/fhwbN27EDTfcgJUrV+Lcc88FwKCNiIioWDBo09LUBNTWArt3AyKJv7W1BQVuq1evxuDBg/H5z38eQGJu0AceeAA//OEP0dnZqfmew4cPo7u7G2VlZQAS84qOGzeu3zr79+/HwoULsXHjRkyZMgXbt29HTU0NWltb8ZWvfAVHjx7FlClTsGDBgrzTTkRERO5j0Kalrg7IDKQ6OxPL87R582ZMmzat37JTTjkF5eXlffN6ZhoxYgTmzp2LiooKzJ8/H01NTf0mnQeAUaNGYeXKlZg1axY2btyIs846q++1b33rWzj55JOxcePGvvk/iYiIyJ8YtGnZs8fccgNEBEopw8tTVq5ciZdeegnnn38+li9fjptuuinvNBAREZF/MWjTUl5ubrkBEyZMQOZ8qIcOHcLbb7/dr3RMy8SJE/HFL34RL774In72s5/lnQYiIiLyLwZtWurrgdLS/stKSxPL83TZZZehs7MTP/nJTwAAPT09+PKXv4wbb7wRpZmfldTR0YFYLNb3/8aNG1FRUWHqcwcNGoTjx4/nnW4iIiLyBgZtWhYsABoagIoKQKnE34aGxPI8KaXw3HPP4ZlnnsHYsWNx9tlnY/Dgwf16d7700ksYPXp03+P111/Ht7/9bYwbNw5TpkzB3XffjR//+MemPre2thaTJk1iRwQiIiKfK3E7AZ61YEFBQZqWj3zkI3jhhRc0X6upqcHRo0cHLJ81a1bO7dbU1PQbCiS9dO6+++7DfffdZzqtREREQdAUj6Nuxw7s6epCeSSC+qoqz43PlsKgjYiIiAKpKR5H7dat6EyOzLC7qwu1W7cCgCcDN1aPEhERUSDV7djRF7CldPb2om7HDpdSlB2DNiIiIgqkPV1dppa7LVBBm4i4nQSyAY8rERHlozwSMbXcbYEJ2gYPHoz29nZ+wRcZEUF7ezsGDx7sdlKIiMhn6quqUBrqHwqVhkKor6pyKUXZBaYjwujRo7F371689957lmzv2LFjgQ4UvJT/wYMHY/To0W4ng4iIfCbV2YC9Rz1m0KBBGDNmjGXbi8VimDp1qmXb85ug55+IiIrDgmjUs0FapsBUjxIRERH5GYM2IiIiCpymeByV69cjFIuhcv16NMXjbicpp8BUjxIREREB/htUN8W1kjal1Dil1Ma0xyGl1Bcy1qlRSr2fts5dbqWXiIiIioPfBtVNca2kTUS2ApgCAEqpMIC/AnhOY9V1IvIJJ9NGRERExSc1z+hunw2qm+KV6tHLAGwXkd1uJ4SIiIiKT2aVqBavDqqb4pWgbR6AJ3Vem6GUegPAOwCWichm55JFREREfpardC3Fy4Pqpii3ZwhQSp2EREA2QUTiGa+dAqBXRDqUUlcBeFBExupspxZALQBEo9Fpq1atsjXdHR0dGDp0qK2f4WVBzn+Q8w4w/8x/cPMf5LwD/sx/M4DlAHJVekYBLAEwJ8s6duZ/9uzZG0SkOtd6XgjargHwjyJyhYF1dwGoFpED2darrq6W1tZWi1KoLRaLoaamxtbP8LIg5z/IeQeYf+Y/uPkPct4Bf+a/cv36nCVsFZEIds2YkXNbduZfKWUoaPPCOG3zoVM1qpQ6TSmlks/PRyK97Q6mjYiIiHwqV8cCP1SJpnO1TZtSqhTA5QD+IW3ZLQAgIisAXAdgqVKqG8BRAPPE7aJBIiIi8oXySES3pK3C4/OManE1aBORTgBlGctWpD1/GMDDTqeLiIiI/K++qmpAj9HSUAgN48b5KlhL8UL1KBEREZHlFkSjaBg3DhWRCBQSpWt+DdgA7wz5QURERGSZ1FAfe7q6UB6J4Inx430brKUwaCMiIqKi4te5RXNh9SgREREVFb/OLZoLgzYiIiIqCk3xeNax2bw+t2gurB4lIiIi3yuGuUVzYUkbERER+Z5WlWg6vw2kq4VBGxEREfletqpPvw/1kcLqUSIiIvI9vdkPjM4t6gcsaSMiIiLfq6+qQmmof1hTDFWi6Ri0ERERka+lBtLt7O1FOLmsWKpE07F6lIiIiHwrs9doD06UsBVTwAawpI2IiIh8rFgH0tXCoI2IiIh8S6/XqN8H0tXCoI2IiIh8S2/AXL8PpKuFbdqIiIjIN1KdDnZ3dSGMRBs2BUDS1im2XqMpDNqIiIjIF7Q6HQCJgC0VuFVEIkXZCQFg0EZEREQ+kW2qqlTAViwD6WphmzYiIiLyhVydC4qx80E6Bm1ERETkaU3xOCrXr+/Xbk1LMXY+SMfqUSIiIvKMVEeDPV1dGBEO41hvL45IrnCteDsfpGPQRkRERJ6Q2dGgvacn6/qp3qPF3PkgHYM2IiIicl1TPI7FbW3IHqadoAB019TYmCLvYZs2IiIiclWqhM1owAYUf/s1LQzaiIiIyFXZhvLQEoT2a1oYtBEREZGrzAzVUVZSgoZx44q+/ZoWtmkjIiIiV6R6iur1DVUARpSU4GB3N8oD0tkgGwZtRERE5LjMnqKZSkOhwJao6WHQRkRERLbTmuhdT1CG8DDL9aBNKbULwGEkjl+3iFRnvK4APAjgKgCdAG4UkT86nU4iIiIyJz1QS03oDmQP2BRQ1POHFsL1oC1ptogc0Hnt4wDGJh8XAHg0+ZeIiIg8KrP6M/ecBglBHMrDKD/0Hr0GwE8k4TUAH1JKne52ooiIiEhbaqBcM8N4AMEdysMoLwRtAuC3SqkNSqlajdfPBPB22v97k8uIiIjIY/IZKBdItGNjx4PslBiYhNXWBCh1hoi8o5QaBeBFALeLyNq01/8HwH+IyMvJ/18C8M8isiFjO7UAagEgGo1OW7Vqla3p7ujowNChQ239DC8Lcv6DnHeA+Wf+g5v/IOcdyJ3/ZgArAcRNbjcCYBmAOfknzRF2Hv/Zs2dvyGzTr8X1Nm0i8k7y736l1HMAzgewNm2VvQA+kvb/aADvaGynAUADAFRXV0uNzfORxWIx2P0ZXhbk/Ac57wDzz/wHN/9BzjvQP/+pDgZ7urowIhzGsd5eHDFQCJTqjODHid69cPxdDdqUUkMAhETkcPL5FQC+mbHaLwDcppRahUQHhPdF5F2Hk0pERBR4TfE47ty2De09Jyo/059nEwbw+PjxvgjQvMrtkrYogOcSo3qgBMB/i8j/U0rdAgAisgLAr5AY7uMtJIb8+LxLaSUiIgqEzJI0KIV2AKqtzXAv0HQcKNcargZtIrIDwGSN5SvSnguAf3QyXUREREGUqyQtn4DNT1WgXud2SRsRERG5SG8A3EKxdM16DNqIiIgCIlcHAqsCtrKSEjw4diwDNosxaCMiIgqAzBkKjHYgyGWIUhgcDuNgdzfKWRVqKwZtREREAVC3Y4fpGQq0pKpQ2VbNeQzaiIiIAmBPV1de72NJmncwaCMiIgqA8kgEuw0EbnolabFYDDUzZtibSMqKQRsREVEA1FdV9WvTpoUdCLyNQRsREVGRS/Ua7ezt7ZtCqiw5aC6rPf2DQRsREVER0ht/rQeJMdQePPtsBmk+w6CNiIioiGjNapA5/lpnby/qduxg0OYzDNqIiIiKROZYbNnk25uU3MOgjYiIyKfSq0BTbdWMKo9E7EoW2YRBGxERkUdpBWWpDgTt3d0D2qoZVRoKob6qyvL0kr0YtBEREXlItg4EALK2VcuGMxn4H4M2IiIij8hsk8YJ3CkdgzYiIiKPsGp+0BSWqhUXBm1EREQuSVWF7unqwohwuF/VZyFKQyE0jBvHYK3IMGgjIiJyQWZVaL4BW6qtWqqjAkvXiheDNiIiIoc1xeNY3NZmqMdnZlDG6aeCi0EbERGRg1IlbEbL1Z4YP55BGQFg0EZEROSI9KE8jKqIRBiwUR8GbURERDYzM71UCgfApUwM2oiIiCyW2Sv0f3t6kCtcUwBGlJSwrRrpYtBGRERUgMwA7VhvL47IiWFxjfQK5RAdZETI7QQQERH5TVM8jsr166FiMSxqa8Puri4IEgFaesBmREUkwoDNa5qagMpKIBRK/G1qcjtFABi0ERERmZJqn5bqUJDvVFOloRAax4/HrhkzGLB5SVMTUFsL7N4NiCT+1tZiVHOz2ylj0EZEVHQ8WkpQLKyYaioMsHTNq+rqgM7O/ss6O1G1cqU76UnDNm1ERMUkVUqQ+tJJlhIAABYscC9dRWSPiSE7tLD9msft2aO5OLJ/v8MJGYglbUREbrCrNEynlAB1ddZsP4BefuQR7D3tNFx86aV4+7TTMN9kNdkQpVBWUgIFtl/zhfJyzcVdo0Y5nJCBXCtpU0p9BMBPAJwGoBdAg4g8mLFODYDnAexMLnpWRL7pZDqJiCxnZ2mYTimB7nLqJ7Mn6Cd/8xs8vHw5hiRL1z4Sj6Nh+XIIgCfnzOl7X2qqKU4xVQTq6/tfnwBQWoodS5bgXPdSBcDd6tFuAF8WkT8qpYYB2KCUelFEtmSst05EPuFC+oiI7JGtNKzQoK28PBEEai0nAP1nJkifz1NrqI67V67sC9hShnR14d6VK/uCtjCAxznVVPFIXYN1dYkfO+XlQH099p95ZnCDNhF5F8C7yeeHlVJtAM4EkBm0EREVFztLw3RKCVBfX/i2fSw9UEuVigHom/9Tbyy1cp12TOnLewEGbMVmwYKBP6BiMVeSks4THRGUUpUApgL4ncbLM5RSbwB4B8AyEdnsYNKIiKxnZ2mYTilBsXdCyKzWhFJo7+7uK0lLD9TMDNGxZ9QoVMbjmstTyiORAlJOZJwSk4MAWp4ApYYCWAOgXkSezXjtFAC9ItKhlLoKwIMiMlZnO7UAagEgGo1OW7Vqla3p7ujowNChQ239DC8Lcv6DnHeA+bci/6OamzFu+XKE06rdeiIRbF22DPvT2kl5kRePfzOA5QAK69OpbX5zMx5La9MGAEciEdy8bBmenDMHEQDLAHj7qFnDi8feSXbmf/bs2RtEpDrXeq4GbUqpQQB+CeA3IvJdA+vvAlAtIgeyrVddXS2tra3WJFJHLBZDTU2NrZ/hZUHOf5DzDjD/luW/qcmXpWFeO/5N8TgWt7Uh90RR+Zvf3Ix7V65E+f792DtqFO6prcXKSy8tzo4GWc5Lrx17p9mZf6WUoaDNzd6jCsB/AWjTC9iUUqcBiIuIKKXOR2KIknYHk0lEZA+tNjNkSmpmArsCtiFKYXA4jFVz5uDVq6/Gwq4u3FNTgwYADTZ9pqs4xp/nudmm7SIAiwD8SSm1MbnsawDKAUBEVgC4DsBSpVQ3gKMA5onb9blEROQJhcxMkGrjlt57NNdQHTEPNES3lZ29mskSbvYefRmJ6ybbOg8DeNiZFBERkR+k9wQ1IxWoVRRjtaYVOMaf53mi9ygREVEuTfE47ty2TXd4jhQFYERJSb/eowzUDOAYf57HoI2IiDwv1X4tV3WEESyMAAAgAElEQVQo5/UsAMf48zwGbURE5FlmqkJZmlaggI7x5ycM2oiIyFVGp5XKpiISwa4ZM2xNZyCwV7OnMWgjIiJb6QVlqVkLzEwrpaU0FEJ9VZW1iSZv8em4hlZj0EZERLbJbIumFZQVMo5TWUkJHhw7llWixYzjx/Vh0EZERJbLd1gOo9h+LUA4flyfkNsJICKi4pIqXbMjYCsNhdA4fjx2zZhRHAFbUxNQWQmEQom/TU1up8h7OH5cHwZtRERkqXxnKpjf3Iyd8+ah59JLsXPePMxvbu73ellJSXEN55Gq9tu9GxA5Ue3HwK0/vXHiAjh+HIM2IiKyTFM8nlcJ2+eam/HY8uWojMcREkFlPI6Vy5fjc83NqIhE0Dh+PA7MnFk8ARugX+23cCFL3dLV1yfGi0unVCLILWQ/+bCUk0Gb1/nwpCKiYEpVi2YTTv4tC4dRVlIChUT7tB88/jiGZAR7pV1daGpsLJ6q0EzZqvdY6nbCggVAQwNQUZH4X6lEySSQ/37yaSknOyJ4GXvMEPlbkQ9ToDWUh56cMxX89a/ay4u53ZLetFEpAW1sryk1flxl5cB91tkJLF58Yj0j9Eo5zW7HYSxp87JsPWaIyBi3Sqt9+ks+XVM8jsr16xGKxTBy3TqMfPllzAZQEotBxWJY1NbWVxWaa1S1nG3RgtZuqakJ6OjIvV6uoDVotTF6+6Onx9z1ZdV2HMagzcvYY4aoMG4GTk796CrgS1srKFOx2ICgTJAYV629uxvAiQDN6PhqFZFI7upNrXZLRua99GPQkjov29tzr5staPXDDwOrj0+2/WHm+rJqOw5j0OZlQfvlGUR+/MLxEzdLq5340WXiSzszQBu6Zg0WWhSUZWN4toL0dktKJf42NGSvptLK/8KFwMiR3r6WtM5LIJHvdLmCVi/XxjQ1JY7DwoXWBpVawX06o9eXVdtxGIM2L8v3lyf5gx9+Jfudm6XVVv3oyhbY63xpd3zlKzkDNKNzehaiIhIxN0THggXArl1Ab2/ib652RXrBT3u7+9dStuOmd/6JmAtavVobk60ksdCgMhXch8Parxu9vqzajsMYtHlZPr88yT+8/CvZD4yUUrpZWm3Fj65cgb3Ol3PpX/9qKEDLNS5avhwbADdbcNLZCdx5pzfbM+qdfxUV5oJWr9TGpK5FpYCSkkTpmlYwnVJoULlgAfD444VXp9fVJY6LjwpHGLR5ndlfnkHnp+pGr/5K9gOj1WJullZb8aMrS2DfFI9j76hRmm/bo7M83XyNcdEeW77cdOCWqtBLlVdkLV1zsn0TkCjp0Qqc7L5P5PpBluu8NJo+L9TGpF+LQKIhfy5WBJVWVac//niix6hfCkdEJOsDwG0ATs21npce06ZNE7u1tLTY/hle5sn8NzaKlJaKJC7HxKO0NLHcQpblvaKif1pTj4oKa7ZvE08ce719p3XMGxsT6yslUlaWeCiVWJbHueFo/pXSzGOvUlK6Zo3Mr6uTjkik32sdkYjMr6sTtLQIWlpkfl2d7IxGpUcp2RmN9r22MxrV3PbOaLTvvWhpEZX8W7Z2rZStWydoaZFwclnFq69K4759xvJix/Wptc1cj7KyvNNh+NjrHDdRqn/aU+dl+rlodj9pbUdv2wVqaWkZuO2yMnP734Z7cj/Z8m7knpvl/XZe+wBaxUB8YyRouwfAWwCeBnAlAGVkw24+GLTZz5P5dygIsizvDgWZVvPEsdf7Usx2zC3a347mX+ecTg+s9IKy1Gt6QV2Pzj7sUSprUJZ3/s1en0YDj8ZG84FDnvcJw3kv5F5U6H3MxvvK5ro680FyZsBsd8CWLe+5gukc7/dF0JbYFhSAjwFYlQzg7gVwlpH3uvHwbNBm068fN3jiizuTkV+3FrA07z48Jzxx7LOVtOkdc4uCekfz39gox08+uV96M0vSsj30StN2R6OyR+e1XPsj7/ybuT7zCTwKLQEycJ8wnPdCAqdC72NW/XjVuDcd1TtnjATETtzbcuW9wNe9ELQZatOW3OC+5KMbwKkAfqqU+rZV1bRFjz0F7eeVRrlmWN1m0U9t+gqRq7u+1jH3QRvCAcNyjB6Nv/vSl7ArGkWvUtgVjeLmZcvw5Jw5hrZXvn+/7vKPfOc7zraHMnN95tNJJ/NaevBB7fyVlZlLXz4Kac9Y6H0s13lu5B6h830ViceNpSGltBRobHSuPbZe3lNzlO7enX1YFR/cI3IGbUqpO5RSGwB8G8ArACaKyFIA0wB8xub0FQ/2FLRfkAbn1OKVHwZ27M/MbQKJL0GtL2C9Y+6xoN7IuGlHRPDknDkYs2oVwqtXY8yqVTkDttRXUkUkgs4zz9Reqbzc+d7pZq5PK7489fKnF8wVEqxqnfP5/iArtHNBtvPc6D1C5/tKQjohQ1nZiXlBU0NouNGgP9u1nOooIXIicMtMo8fuEZpyFcUB+CaACp3XxhspznP64cnqUYeq7pziiSoyLWarG/OoxvBs3r3Qps/K9jSpY5m6TvS2aabtkwVp21xXZ7pKu3HfPql49VVRyQb9Q2IxQ9WbRh+6bdAKzbPGvi3o/Dd6rPTO5bIya5oT5NksQTPvdnWwyDef2dJj9B6h1wEmtS0r82olMx1T8mj36oXqUdcDLDsengzafNpTUI+pth1ebrOVx3HR7EHlhXx5oU2fle1pct18w2Hz+73Q49bYKN0ZDfuzfWk17tsnZWvXWhqgpT9K16zJ3Xsz3zzrfIFtrquz5/NyffagQSInnWR439tB89z34r09PUALh0+kR+9ayrxH6Kx7NBr15r0vXbbg1Mh90e+9R/348GTQ5tOegnoM5d8Pec4j0NHsQVVaKrJ0qbs3swK7sxuV9dhbFTgauem6cT7ppOvw6NH9StJSQ2MoGwK11DZNDbdhYV6PRqP677GjpDVXxwIzvSrtOPe9WouidSz00pq5D/MN2L0knx7mOXghaOPguk4J4uwGbrfjs2nE/KqVK7XztWKFu+3JjAzYaXebN6vahBhtu5Tv+ZRvuzuDMxCk5vAU8ynLqqykBE+MHw+pqSlstgEj+dfJa0SncwMAa6/5zDZhBw9qr2fkXLHz3PdqOyitYyFibH5Tne+r/QY7wHhCtv3v4RkPcmHQ5qSgzW6QrTGx3Y3/jd6k82j0q/ulJRlf0U53NMn1w8CJINqqEdrNfOGZ7dll5NzQOz910mVkBgKzhiiFspISKCQ6FTSOH48DM2cWHqgpBSxalPva0MlrV3peM/dTqrF3Jit63xUSHNl57nthVgIthc5vasf3lZOdvvR6mJeV+bvAxEhxnN8enqweLTKG8p+tMbHd1aZm2pmYrDYxPVaRE1WmRvJgUTVOzmNvV5smq6o5cp0bOlVD677/fbn97rtzzkCQz2NILCZl69aJsqva08j+NFhF9vY11+h3EDFa/WZVHuwe+yzjXNatHvRiOy8b7oH9rn0HOn0VzOLj4oXqUVeDKyRmWNiKxIC9X9F4PQLgqeTrvwNQaWS7DNrsV1CbtkLbphhhYzsTzTZt2T7P7puU0Zuh3k3cZIN+x879zBvu0qXW3PSzTA1V8eqrugPSHg+FZH5dXdYZCDwRoGnJt2G2xjEY0BHDyXNer4F9ru3n01lA47rqjkSMBSdeCOCM3heMrJfMU28h16IXO2yYFOigDYn5hbcDqAJwEoA3AJybsc6tAFYkn88D8JSRbTNos19BvUetCqiy3RxtvEFo9h7Vuollawhrxc08Vy8pIyUneXyxunruW/GFqLPPUlND6U3vJNAvVcsWyDnWcSCbXI2yjV4bRjuIZB4jKwMZE0FGQQG/RYGeq52vjOz3fEqe8ylR9WqHDROCHrTNAPCbtP+/CuCrGev8BsCM5PMSAAdgYO5TBm32Kyj/VgRUuW6ONt48dfOeeYPM9cVWSHqMVHfplZykSijy3P+2T+Fmd0mFxr5LD8b0Stoyg7ts83seiUTkc3V17gZq6XKdj4VWM+p92etVoxZy7udZvW26d3c+QYYfS5Ny5dNooF6M+yaDF4I2lVjXeUqp6wBcKSJLkv8vAnCBiNyWts6m5Dp7k/9vT65zQGN7tQBqASAajU5btWqVrenv6OjA0KFDbf0MLysk/6OamzFu+XKEu7r6lvVEIti6bJnh3kkXzpuHwRpTqhyLRvFa8tiPam5G1cqViOzfj65Ro7BjyRJLej8ZzbteGvXSa0Yh277k0kuhNK57UQprVq/O+dlmj72Z423FuWE0TVUrV+Kk/fuxZ9QofG3Jkr6ZBuY3N+Ox5csxJC0N6XqVQjhtP+2cNw+VOc5Ft2nt19QZ0BWNGr42cp13qWMFYMDnZcp3/+Q6f43cG4zIZzuFXltuyJVPvTxpybZvnLq27WTn9/7s2bM3iEh1zhWNRHZ2PABcD2Bl2v+LAPxnxjqbAYxO+387gLJc22ZJm/0Kzn+20hQHG9Xno6BJo61Kb64Sj2wlGQWO52a4pDH1HjO/sB3+Na43jtr8ujo5HgpppqU3WdqWKpnTrU71WrWPRR1EBrRpS+XfyDG3Yv/kOkcs6nSgVaWas02bH0uTctVK6OUpn9JTr7T3y5MXStocC9IGfDCrR33NtvwX2qjegZujqbybbXdmdJt6VZyZX556789WfZTjhmx6Kh8zX6IOBUCpaaWytUv7z2uuGVDtmf7oiETksU99Sje48/QXdQEMTeNlVRs6LfkGGfm0RcuoUjU0G4SX2rQZletHtBcHE3dB0IO2EgA7AIzBiY4IEzLW+Uf074jwtJFtM2izn235N3rDdfHmmHebLivSa1FnAkMNtXWOgempfFwoacuc6zPVW1Nv3k+9dmn/ec01sisaTcy5aKS0wU9f1HkqaLgfq/ZPPkGGBSVkhnvNF1swk9l71Mo8+Wh/BTpoS6QRVwHYlqz2rEsu+yaAucnngwE8g8SQH78HUGVkuwza7Gdb/s2UtORzsVtwg8g77zb2fMxrHk4j29U4Bqan8jHzJVpgcJvvXJ+6nQ9yVbnZcSw8Lu/hfrSqUe1i9lrLdnyNNA1wmkuBjuX591nJpBeCNldnRBCRX4nI2SJylojUJ5fdJSK/SD4/JiLXi8j/EZHzRWSHm+m1hZMjRHtdU1NiP2jRGvXczIjdTU3AyJHAwoXuTTVlxQjjeqOc9/bqb6+AKYv6EQEqKzGquXnga9lGq882U0Nm2gDD0701xeOoXL8eKhZDSSwGFYthUVsb2nt6cuclM5l6s1yk9ouZGRqyHQs/sOKepHXMn3gicQ45MRuM2Wst2/G14h5h5X3eienonOL2VIc+xGms3FRMF1+hUvtC6wu30ClhUttubx/4mt9uEGan8jF6jhkNSnbvxrjly81PB6b1JaqXNmDAunoB2u5kT7TUWSPGcjGA7jRUqf2ilb/MORwz3+NHVt6TrJoGyYkftnpTHgGF3yO09unChYkfkfnkpRgCndQxtXPasyLFoM1NxXDxAdbcVLX2BQCEw4XPE6e37RQv3CDS9+HIkYmH1v40M89hUxOweLGxc8xEUBLu6jrx/lS6Fy0CTj45Ma9frjkNU3Kc/+mBmpUBmpZ/u/lmdJ98cv+FmUFnsuRIUvm75RZvzjlZCK/dk5z6YZs6vnq07hFG73t695/2dnN5KZZAJ/2Y6vHzDx+bMWizSyFVUn65+ADrfkXmU+1X6LZT3L5BZO7D9vbEQ+tLKtek8Jnb1KsqzNwnetVZeqVJe/Zop/vo0cT7jJSs6BwX2bMHI9etw8K0QM3KAC1TWUkJ5tx2G0oeeyz7fk2WHK1ZvTqRv0ceMVyV6xteuyc5GUQuWJA4hloy7xFmgsls+85oXoop0Mn1I9rvP3xsxqDNDoVWSeW6+Iz+wnOiWsGqX5H57otCtg144waR6yaWeWM3Uu2Ua5tG2wjq7bsRI4yX4plJA4Ddo0bl1S5NzxClUFZSAgWgLBzue14RiaBx/HgcmDkTC6LR/KrzvFYFWOh27LwO8+F0EGm0JNtMMJlr3xnJSzEFOtnym/rhA7Cttw4GbXYwekGbqepKMRoQOlWtUOivyPQi/8xSnWzVfmYuaL32KmVl3igZMXLTNvsllW19Mzd4jX3XGw4Dhw8bL8XDiarOUCyGkevWYeTLL+NzCxfiSCTSb70jkQi+tmSJsbTpSJ1FqaCs45JLcGDmTPTW1ODArFl9z3fNmJEI1txm1bVa6HaamoCOjoHL3QwInA4ijZZkmwkms7WXA4zlxUig4/Z9zCi9/FZUJH74AGzrnY2RLqZ+e7g+5Iedw1YYHcfKysFnNdLYl/9cQ0VkGxQ1n2EB8u0ibmEXecu7fRsZbsPscbNyaJCMfdd1yilZ07ozGpVwcigNvXHRjEy0buSRmtEg9XlOzPXp2PG36pgb2Y7e+H9lZQPOF0eHvfDYkBA573t6+7qxMbEvM9c3eu8qcL5gq1hy7PMdHNmp4WKy8MKQH64HWHY8XA/a7Byt32hAaNXI8joXWN/I4LmmasqW53z2kwemiXFkrKJCv6Ts+LJLDbCZJWBLn3jd6kDNjQBNi+XH36prtZDtmLiuHB+rzEODr/bl3akfj1YNpm1Reiw79tk+t5Ap+mzGoM2mh+tBm52/DgstaTP7a0VnO0ej0f75zedXZD5fMtkuaIdu7LZcuOk3sbKyE/sz9Qs7tczhgYT7bSvHjAnHQ6EBAZiRqaGMBHpl69a5EqBpKcqSNhPXomcGmHVBv7w7EUzaNZi2SF7fU44ceztqHizCoM2mh+tBm4h9F7TRC83Al6xW1ccAOjfzXiereo28J5/Ji/PkyI3LjtK3QuS4kWoFXlpTQ/XovH9nNKpZquZWaVo2eU9jZuW0S3qfke92vFzSZie3SpqMsnMu3jzuv56491mV/zwwaLPp4YmgLZ3RG4Nd62U7+fOck69fSVu+8vmSydYOzuyvsTwD60D+2tQL3gHZpVPFqTs1lMajRynXqz2NMn38jZznVv3Iy3c7JtKoOf+kh6owDfNqSVM6LzS1SeNY/nN9d7GkrbgengraCikZs6o0JVcbgWwXQK42bYXK52af+Z58fo3l2t9Z0uXIjcvIXJdO/trU2c+ZJWTpjx4z83W6dBPOh+nj74F2mIbkWxrosc4Chnm1pCmdF5rapLE9/5nn4NKlnjq3GLTZ9PBU0OZUb89sN9xCenjqbNtTVSRWd2jIcaMstpK2dd//vrydbHeWKjUrW7tWytatE5XsAXqjRlVntrZo4WwlbQ5WZdvB9PG3sze5U7JdL/neu9zOq5dLmtK53dQmja3510vP0qWeuSYYtNn08FTQ5kRvTyOlRvn28NThqaAtn1+j2fZ3ji8hT7TrMBHoNO7bJxWvvtoXgJWtW9cXWGm1O9MLxoz2+lS5eth56CacD9tK2rxcYpXtesnn3uWFvHqxpMlpXmrT54MSaS8EbRxc125GB4csZBDJXIP5pgaMLCsb+F6zA2cmB7a95NJLvTNStdEBMdNl299emMonM09lZYbm9cwcxHbomjV900EJgPaeHrR3dwNIzOF578qVGJKcKiplSFcX7l25csC2n5wzB2NWrUJ49WqMWbUKT86Zo5n08tSAuXrH5ZFHrJlFwC/sGGXfadmul3zuXV7Iaz6Dmxcbq2b0sIIX7rs+wKDNapmj9V91lbEbQyE3EL2TevfuE4HVggXAgQNAY2P+cyWmjbiuRJwdqTrXLAhmbz7Z9rdXpvJJz9OBA4lHlvw1xeOo3bq1X4B2RCTrR5Tv329qeS6loRDqq6q08+Dkl4KZWTOsnu4tc3uAdaPsOzE1nZZs10s+9y4vfEHn82OP7FPofdeta8NpRorj/PZwonp0c13dwGLlQquD8m27UGgPUaPcKr62qypFb397oU2bSY379vX1vjTz0Gt3lq2DgdcGux24M0ycL1a367Fz2A23qxSt7D3qg6owLV689p3kSps2o53THLg2vFA96nqAZcfD9qCtsVG6M9oBSWmp9gCzbgU1RtJg9kZr55hB2bhxg0/fNxkD21rWczbfpCXbqKUHTCqPgE1vLDW/DXY7gJnzxep2TYWcq7m+eDwS6DgylZFHMWhrsfcDrC64sPja8ELQxurRfNTVIZzRDgidnUB7u/b6dhf5pxfz68lMg97k0rfeql/E7Fa1oRtVKalqvSeeAI4eTRzb5H4aX18PjBzpSvF7ehUokGiXBgDZK0H1PTlnDm5etgy7olH0KoXd0ShuXrYMv/3Yx1BWUgIFoCwc7nteEYmgDsCBmTO9Mdm6FjPni9XnViHby1Vd54UqRauwapK05NucopiujRwYtOXD7IlgZ1CTqsdftCjxv1ZnA6006DUEXrGifyC3aFHipmqmfZ7V3GxjprGfFJAI4vJtz6fR9iKzA8HIl1+GisVQEotBpS1b2NaGzt7evLMzRCmUlZQAAMLJZa9efTVeeeMNhHp7UbFvH/77nntwYNYsHJg5E701Nf2e75oxA9rdDzzEzPli9blV6PayfWl5pa2lVdLzWl+fuNb80B7Jq22nvJouJxTbtZEFg7Z86J0IZWXOBjVapWWHDgEnnZQ7DXqBp4j2/7t3A48/DixeDFRUQPLtyJDPTcXNXl7ZAvRsvd308qpxzLpvvhnNDz+s28MTGcuMUkC/ErLG8ePRccklODBzJqSmBt01NZBkIObZUrN8mDlfrD637DxXi7W3o16pvweDjlHNzd5Mq4/2oS2K9drQYqQO1W8P19q0pTojODUGlV49fllZ7jQYGbw1SxsBW6byyfV+N8b2ymdg4mx5zWNmgXwepWvW2NbmzBftesycL1aPVWXnlFRuD0grNhx/j7TVM+Ko3oDRbqe1CNp0FcyBa8MLbdpUYt3iUl1dLa2trbZ+xpavfx3nNjYmSmLKyxMRvdPtMUKhgSVjQKI6M1cVWuqXWXrVn1La29PYdiwWQ01NjfG0VlYmfv1lqqhIVJF4ldZ+SqeV/mx53bNHcx/3KoXw6tUFJVUh0batIhJBfVWVbaVnpo99kXEk/1rnXWmpJ9p9WZ7/Qu5jDpNQKDHcUSa302pkHzY1JWoGCvjO4rVvX/6VUhtEpDrXeqwezdP+OXPcH5SwkHp8rYbAt9wysIg5n21rybehqNvtNPIZmFgnT7JnD/aOGqX52h6d5bmk2qVVRCJ4Yvz44qzuDCIvDD7rFB+1R+rSu07dTmuufRj06tMiwqDNzwqtx89s9PzII/17oSqV/7Yz5XNj9sqNJmNg4pzt+XTytHvUKPzzkiU4kpoxIOlIJIKvLVliKkmloRAax48v3nZpQReg3nB+ao+0Y8kSb6Y11z4M0o+AIsegzc+s6jafXppVV5e40EUSw11Y1SU/nxuz1240ySB3zerVmqWrqR6gn1u4UDcwyxxiY1dyiI30KaFSHQiAEyVpmcNuNIwbxyCtmPmo9KlgPhr+Y/+cOd5Ma5CGiwm4ErcTQAVasKCwG0Zm25lUaZYV206X2o6ZNhV+uNEk24nInj2YNWoU/m8yMAMS83qW79+PPaNG9QVsQGJsNL15O1MOzJxpe9LJw+rrtdu0uV2iYxcr7zV282pas6WrvFy7nW0x/ggocgzagi5baZbVNyazNzsP32ia4nH87tFH8R/f+haGdHVBASiPx/HY8uUAjAVmesozSukogPL5kUOkJ2g/AooYq0f9xI5G+V4uzfJYW5dU9edsAIva2vClFSswJGNmjCFdXbh35cq8P2PAhOtB4XaHEy/Kd3R4okw+qoKm7FwpaVNK3Q/gkwA+ALAdwOdF5G8a6+0CcBiJMUa7jXSHLVq5qjHz5eHSLC+VNjTF47hp81Z8EEp0nxcA5fv3a66rtxwA0ItEo7VeJH4yHQpDKQUZ2g28F0HnY1W4c2MUdyIx6UI4DPT0nOi8evAgMGJE4rnW61Yv0/68Syzd9rXHmvDAkVoMwYlzu3NhLZYsBJ4OL7A5L/ls+xKHP89r54Gx/PsjL2a3k/3c93ZeFgBYgHYBwnuBnoVA2Z1mP++SIj2ns7+noiLx1XPmmXCdK+O0KaWuALBaRLqVUvcBgIj8i8Z6uwBUi8gBM9t3Ypw2x8ersWucszzHg8qZfwvGBPKKW5+N49HhbSd6BSTtnDcPlfH4gPV3RaMYs2rViQCtB4n3xiPAyirgJXYgSLcTlajEwHN7FyowBrucTxARUYbSUuCLX9yCe+4515bte3qcNhH5rYik5uR5DcBoN9LhK3ZVY9pRbO6VoTrykaymk1AIuz98Gj73ta/j0Q8NDNgA4Gt6w3f8/RJgXwS4dzxwaQ1weU3i7/wZDNg0lEP7HNZbTgnz0YSdqEQPQtiJSsyHD64vIp/q7ARWrnS/6YoX2rTdBODXOq8JgN8qpTYopWodTJP32DkEgNVtZ7w2VIdRacGmEkHFgTge++5yzF/drLm65vAdd/wTnvzd7QzQTNgD7XNYbzklArbHUItK7EYIgkrsxmOoZeBGZKP9+93vJGZb9ahSqhnAaRov1YnI88l16gBUA7hWNBKilDpDRN5RSo0C8CKA20Vkrc7n1QKoBYBoNDpt1apVFuVEW0dHB4YOHWrrZ6Qb1dyMccuXI5zW8L0nEsHWZcsSYwc5LFv+L7n0Us2pXkSpxBhnHnXhvHkYnK26MxOrPy2RCkD62rQBOIJS3IwGPAlzPyDmown3og7l2IM9KMfXUG96G37AKmUi5334w514+unf27Lt2bNnG6oedW1SdwCLAawHUGpw/W8AWGZkXdsnjBeXJs71wGTRKZvr6vTT4qMJoNP1KqWZ7h6lBk7K/mKL4LJ9WeeS58P4Yz4aZScqpAdKdqJC5qMxr210oLTfwg6U5rUtr+e1BzrnKpTr+eODj2J8lJaK1NVttu37BwYnjM+5gh0PAFcC2ALgw1nWGQJgWNrzVwFcaWT7vgjasgVgHgrONDU2SnckMvCMTqWzsTHxv97rHrP0Z/sk/MyrsnNUVPNq3RmN9g/Yfr1GTvtSV8gAABTvSURBVLpqn5SVJVYJhxN/y8oSD6VOPNd6vZBldm7b+Of1ejIvO1Ghefz2hCos/rzevPM3H41yRCOwvHlIo6k05MqrveeBsfz765w2uiz7ue+vvOTzedZc+97Ii/H3pL6G7Sys8XrQ9haAtwFsTD5WJJefAeBXyedVAN5IPjYjUa1qaPueD9qyBTV+CHiMlKR5PfBMWvqzfYJfrxG0tMj8ujrpyAhGOyIRmV9XJ3ipRbC6RcLPvCrXPNTqdrJd5UopsxE6JaWilKUfU1D+rSqFdvE+4dnj74Ag512E+fdC0ObKOG0i8n90lr8D4Krk8x0AJjuZLsfkaqjv1AwF+TLSk9UDU73c+mwcDb070DOiCzicNh6aoG+MNAzv6esZqjv91Ow5WPq38Xjk2kRbtVgs5kp+KAcvjzmYYlUvcA+NYUhEzuE0Vm7I58bthRkKUjz45dgvQEsFZR/Cif7Rw3sgmW8a3jNgOwOmnzoWwtL3x/UFbORhfpiqx8prxwM/jIjIWQza3JDrxu2xgGiA+nr0/P3f9+vJ6uSXY84ALUUV8CEChNsjqA1VMWDzCz+UPvkhsCQiz/LCOG3Bk21OTY/Nt6lpwQJsXbbMlXnsbn02jkdLt6JnZFfi7A0jEZxZeSYfC2Hp/45H93UzGLD5jdfn6+QckERUAJa0ucFIiYCXSwsA7J8zB+fec4/jn9vQuwMY3GvfB/QASztZHUo2YrUmEeWJQZtbst24eVPX1TOiK/dK+ToWYsBGRESexaCNfCHVjg1lJt+YmrXgkHbv0dSy8EG2XyMiIm9j0GaFpibPV2f6WaodW85q0VSA1gsgxI4ERERUXBi0FSo1yXiqN9ju3Yn/AQZuFrj12TgeHd7WN5ZaPwIGaEREFBgM2gqVbaBcBm2mGB7KI0UAuazGsfQRERG5iUFboawa4TzgdKtAs4y1Fj4YsTdRREREHsKgrVAenB3AT/pK18q6zA2GeyyE2lCVbekiIiLyGg6uWyg/DIbrQbc+G4d6fh0ePbUtMVCumYCNY6kREVEAsaStUH6YOschRiZo71uWra1aNhxLjYiIAopBmxU4GO7ANmk6E7QPWGZEcigP9hAlIqIgY9BGBcs6LEc+OJQHERHRAAzazEoOpHtJQKtCT1SBAngplntYDrNY/UlERKSJQZsZaQPpKiBwA+nmMyyHYQKowyW4pXssAzYiIiIN7D1qRraBdAOgoXdH7qmkzOgFIED4QARL/3c8eufOZMBGRESkgyVtZgR0IN1+Y6kZ1QOoIyWcoJ2IiMgiDNrMKNKBdHWH6jgcBgb1AqeK6YFv2S6NiIjIWqweNaMIB9JNtVPrGdmVOBuG90BO6e57jlKDAVt6VScDNiIiIsuxpM2MtIF0Zc8eqCLoPZp3OzUOy0FEROQoBm1mJQfSXROLoaamxu3UFKxnhIl2amlCB05Cz/X/1+LUEBERkR4GbQHRr91aWscADM1jY8dC+OS+iMUpJCIiomwYtBWxfoGa1gC4w3vMbTBtLLXPTmyzKplERERkAIO2IqU7EG4u6UN1HNYfoiMWY9BGRETkJAZtRSrvDgYK6J070/oEERERUUEYtBWJzLHWUGay6jMpfJBt1YiIiLyIQZvP3fpsHI+GtwGn9pwYT81sW7WUYyHUhqosSxsRERFZx5WgTSn1DQA3A3gvuehrIvIrjfWuBPAggDCAlSLyLccS6WE5Oxjo6UUisEuOr8ZppYiIiPzDzZK2B0Rkud6LSqkwgO8DuBzAXgB/UEr9QkS2OJVAL8q7g4EAS/82nkEZERGRT3l5GqvzAbwlIjtE5AMAqwBc43KaXJdvB4Nwe4QBGxERkY+5GbTdppR6Uyn1Q6XUqRqvnwng7bT/9yaXBdatz8bRU5bHDAZsq0ZEROR7SkTs2bBSzQBO03ipDsBrAA4gMTb/vwM4XURuynj/9QA+JiJLkv8vAnC+iNyu83m1AGoBIBqNTlu1apVVWdHU0dGBoUPzmU4gP9/70zA8P/ZI7lI2AXBUAcfDwLBuhNpPwif3RfCFiYctTY/T+feSIOcdYP6Z/+DmP8h5B5h/O/M/e/bsDSJSnWs929q0icgcI+sppR4D8EuNl/YC+Eja/6MBvJPl8xoANABAdXW12D0vaMzhuUfnHFivH7AlOxg4OXG70/n3kiDnHWD+mf/g5j/IeQeYfy/k363eo6eLyLvJfz8NYJPGan8AMFYpNQbAXwHMA/A5h5LoOboTu7ODARERUSC41abt20qpPyml3gQwG8AXAUApdYZS6lcAICLdAG4D8BsAbQCeFpHNLqXXdXqD3rKDARERUTC4UtImIot0lr8D4Kq0/38FYMD4bUHSb0y21PhqKexgQEREFBicEcGDsg6e60L7NSIiInIfgzaPyTl4bggIH4ig+7oZziaMiIiIXOXlwXUDycjgubqdEoiIiKhoMWjzGCMBmV6nBCIiIipeDNo8JmdAxs4HREREgcSgzWNqQ1XAsYzD0gtAEm3ZlnaOY+cDIiKiAGJHBI955Noo8CzQ0JHoPRo+yF6iRERExKDNkx65NopHwCCNiIiITmDQ5iHp47OxhI2IiIjSsU2bR6TGZ+sZ2QWEgJ6RXXi0dCtufTbudtKIiIjIAxi0eYTm+GyDexPLiYiIKPAYtHmE3vhsHEiXiIiIALZpc02/+UUFuuEzB9IlIiIigEGbK3LOL5rCgXSJiIgoiUGbC3LOLypAuJ29R4mIiOgEBm0O6qsSLcvRTk2A7utmOJMoIiIi8gUGbQ4xXCUKtmMjIiKigRi0OeDWZ+N4dHgbEDawMtuxERERkQYGbTbp1zv0Q9AfXEWQmBA+xHZsREREpI9Bmw1MVYW2R9h+jYiIiHJi0GYxVoUSERGRHRi0WcBwVWi6HmBp5zhWhRIREZEhDNoKZKYqtM+xEAM2IiIiMoVBW4FyDpSb0gtAsbMBERER5YdBW4EMTejeAyx9fzwDNSIiIsqbkdZXlEXOgXCPhRiwERERUcEYtBWoNlQFHMvYjb1IzB96IMK2a0RERGQJVo8W6JFro8CzQENHovdo+CDbrBEREZH1GLRZ4JFro3gEDNKIiIjIPq4EbUqppwCMS/77IQB/E5EpGuvtAnAYQA+AbhGpdiyRRERERB7iStAmIjekniulvgPg/SyrzxaRA/anioiIiMi7XK0eVUopAJ8FcKmb6SAiIiLyOrd7j84CEBeRv+i8LgB+q5TaoJSqdTBdRERERJ6iRMSeDSvVDOA0jZfqROT55DqPAnhLRL6js40zROQdpdQoAC8CuF1E1uqsWwugFgCi0ei0VatWWZENXR0dHRg6dKitn+FlQc5/kPMOMP/Mf3DzH+S8A8y/nfmfPXv2BiPt9m0L2nJ+sFIlAP4KYJqI7DWw/jcAdIjI8lzrVldXS2tra+GJzCIWi6GmpsbWz/CyIOc/yHkHmH/mP7j5D3LeAebfzvwrpQwFbW5Wj84B8Ge9gE0pNUQpNSz1HMAVADY5mD4iIiIiz3AzaJsH4Mn0BUqpM5RSv0r+GwXwslLqDQC/B/A/IvL/HE4jERERkSe41ntURG7UWPYOgKuSz3cAmOxwsoiIiIg8ye3eo0RERERkAIM2IiIiIh9g0EZERETkAwzaiIiIiHyAQRsRERGRDzBoIyIiIvIBBm1EREREPsCgjYiIiMgHGLQRERER+QCDNiIiIiIfYNBGRERE5AMM2oiIiIh8wLUJ4/3ue38ahjkH1qNnRBfCByOoDVXhkWujbieLiIiIihRL2vJw67NxPD/2CHpGdgEhoGdkFx4t3Ypbn427nTQiIiIqUgzaTLr12TgeHd4GDO7t/8LgXjT07nAnUURERFT0GLSZcOuzcTxauhUIa7/eM6LL2QQRERFRYDBoM6Ghd8fAErY04YMRB1NDREREQcKgzYSsJWnHQqgNVTmXGCIiIgoUBm0m6Jak9QBLO8ex9ygRERHZhkGbCbWhKuBYxi47FsLS98czYCMiIiJbMWgz4ZFro1jaOQ7hAxGgFwgfiLCEjYiIiBzBwXVNeuTaKB5BFLFYDDXXzXA7OURERBQQLGkjIiIi8gEGbUREREQ+wKCNiIiIyAcYtBERERH5AIM2IiIiIh9g0EZERETkAwzaiIiIiHyAQRsRERGRDzBoIyIiIvIBBm1EREREPqBExO00WE4p9R6A3TZ/zEgAB2z+DC8Lcv6DnHeA+Wf+g5v/IOcdYP7tzH+FiHw410pFGbQ5QSnVKiLVbqfDLUHOf5DzDjD/zH9w8x/kvAPMvxfyz+pRIiIiIh9g0EZERETkAwza8tfgdgJcFuT8BznvAPPP/AdXkPMOMP+u559t2oiIiIh8gCVtRERERD7AoM0kpdSVSqmtSqm3lFJfcTs9dlNKfUQp1aKUalNKbVZK3Zlc/g2l1F+VUhuTj6vcTqtdlFK7lFJ/SuazNblshFLqRaXUX5J/T3U7nXZQSo1LO8YblVKHlFJfKObjr5T6oVJqv1JqU9oyzeOtEh5K3g/eVEp91L2UF04n7/crpf6czN9zSqkPJZdXKqWOpp0DK9xLuTV08q97riulvpo89luVUh9zJ9XW0cn/U2l536WU2phcXlTHP8t3nbeufRHhw+ADQBjAdgBVAE4C8AaAc91Ol815Ph3AR5PPhwHYBuBcAN8AsMzt9Dm0D3YBGJmx7NsAvpJ8/hUA97mdTgf2QxjAPgAVxXz8AVwM4KMANuU63gCuAvBrAArAhQB+53b6bcj7FQBKks/vS8t7Zfp6xfDQyb/muZ68D74BIAJgTPK7Iex2HqzOf8br3wFwVzEe/yzfdZ669lnSZs75AN4SkR0i8gGAVQCucTlNthKRd0Xkj8nnhwG0ATjT3VR5wjUAHk8+fxzAp1xMi1MuA7BdROweuNpVIrIWwMGMxXrH+xoAP5GE1wB8SCl1ujMptZ5W3kXktyLSnfz3NQCjHU+YQ3SOvZ5rAKwSkS4R2QngLSS+I3wrW/6VUgrAZwE86WiiHJLlu85T1z6DNnPOBPB22v97EaAARilVCWAqgN8lF92WLBb+YbFWDyYJgN8qpTYopWqTy6Ii8i6QuNgBjHItdc6Zh/437KAcf0D/eAftnnATEqULKWOUUq8rpdYopWa5lSgHaJ3rQTv2swDEReQvacuK8vhnfNd56tpn0GaO0lgWiO63SqmhAH4G4AsicgjAowDOAjAFwLtIFJsXq4tE5KMAPg7gH5VSF7udIKcppU4CMBfAM8lFQTr+2QTmnqCUqgPQDaApuehdAOUiMhXAlwD8t1LqFLfSZyO9cz0wxz5pPvr/aCvK46/xXae7qsYy248/gzZz9gL4SNr/owG841JaHKOUGoTESdwkIs8CgIjERaRHRHoBPAafVwtkIyLvJP/uB/AcEnmNp4rCk3/3u5dCR3wcwB9FJA4E6/gn6R3vQNwTlFKLAXwCwAJJNuhJVgu2J59vQKJN19nupdIeWc71QBx7AFBKlQC4FsBTqWXFePy1vuvgsWufQZs5fwAwVik1JlnyMA/AL1xOk62S7Rj+C0CbiHw3bXl63f2nAWzKfG8xUEoNUUoNSz1HolH2JiSO++LkaosBPO9OCh3T71d2UI5/Gr3j/QsAf5fsSXYhgPdTVSnFQil1JYB/ATBXRDrTln9YKRVOPq8CMBbADndSaZ8s5/ovAMxTSkWUUmOQyP/vnU6fQ+YA+LOI7E0tKLbjr/ddB69d+2732PDbA4keI9v+f3v3G5pVGcZx/PvrD8T6MzIiVqARhJBhD82skUkvepNkCUUro1IikIhlvYhAqKgokaAUozCLqQVqWjIsyhjW0oEJNp2touhNSBBIkEh/4erFfQ1Oo80arcez/T7wsLPr3LvPfe88Z7uec57nXJRXFSuaPZ7/Yb7zKKd8DwED+VgAbAIGM94DtDV7rBM0/0sonxA7CHw+vM+B84Be4Ov8Oq3ZY53A30ELcBRorcQm7f6nJKffA79TXk3fN9r+plwieSn/HgwCc5o9/gmY+zeU9+4MH/+vZNtb85g4CBwAFjZ7/BM0/1Gf68CK3PdfATc2e/wTMf+MdwPLRrSdVPt/jP91J9Wx74oIZmZmZjXgy6NmZmZmNeCkzczMzKwGnLSZmZmZ1YCTNjMzM7MacNJmZmZmVgNO2sys1iRdLKnp94mT1JC0oPL9zZIea+aYzGxycdJmZjZC3gH+32pQ7usEQET0RMTK/25UZjbVOWkzs1qR9Iikw/lYnuHTJG3Iot7bJLVk25WShjL+fMbOl7Rd0v58XJvxJyWtk7QL2Chpn6RZle1+JKld0lxJ/Vkou1/SzKyQ8hTQKWlAUqekJZLW5s/OkNSb4+iVND3j3ZLWZD/fSrot422S+rKvw5OpGLeZjZ+TNjOrDUntwFLgauAa4H7gXGAmsC4iZgM/AQ9ImkYpOzQr489kN6uBFyLiKspd3ddXNtEO3BIRi4HNwO253Tbgwig1Fr8E5kcplP048GxE/JbLWyKiERFb+Ku1wMYcx5vAmsq6Nsrd2G8Chs/MLQY+iIgGcAXl7uxmNsWN5xKAmVmzzAPeiYjjAJLeBq4DvouIvdnmDaALeBH4BVgv6V1gZ66/AbislBoE4Jzh+rJAT0T8nMtbgQ+BJyjJ21sZbwU2SLqUUvbm9H8w7g5KwW0oZZFWVdbtiFKMfEjSBRnbD7yeBax3RISTNjPzmTYzqxWNEh9Zjy8i4g9gLrAdWAS8n+tOATryjFgjIi6KiGO57nilgyPAUUmzgU7KmTeAp4HdEXE5sBA4YxzzqI7318qyctt9wHzgCLBJ0j3j2IaZTTJO2sysTvqARZJaJJ1Jufz5CTBdUke2uRPYI+ksSpH794DllA8KAOwCHhzuUFKD0W0GHs1+BjPWSkmmAJZU2h4Dzubv9QN35PJdwJ6xJilpBvBDRLwKvAZcOVZ7M5sanLSZWW1ExAGgG/gU2Ed5P9qPwBfAvZIOAdOAlykJ1M6MfQw8nN10AXPyQwFDwLIxNrmNkmxtrcRWAc9J2gucWonvplx2HZDUOaKfLmBpjuVu4KETTPV6YEDSZ5T33a0+QXszmwIUMfKqgpmZmZmdbHymzczMzKwGnLSZmZmZ1YCTNjMzM7MacNJmZmZmVgNO2szMzMxqwEmbmZmZWQ04aTMzMzOrASdtZmZmZjXwJ7LN5/x5yIHPAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "c_ols = np.linalg.lstsq(X_ordered, y_censored, rcond=None)[0]\n", "fit_ols = X_ordered.dot(c_ols)\n", "plot_fit(fit_ols, 'OLS fit')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that we are systematically overestimating low values of $y$ and vice versa (red vs. cyan). This is caused by our use of censored (blue) observations, which are exerting a lot of leverage and pulling down the trendline to reduce the error between the red and blue points.\n", "\n", "## OLS using uncensored data\n", "\n", "A simple way to deal with this while maintaining analytical tractability is to simply ignore all censored observations. \n", "\n", "$$\n", "\\begin{array}{ll}\n", " \\underset{c}{\\mbox{minimize}} & \\sum_{i=1}^M (y^{(i)} - c^T x^{(i)})^2\n", "\\end{array}\n", "$$\n", "\n", "Give that our $M$ is much smaller than $K$, we are throwing away the majority of the dataset in order to accomplish this, let's see how this new regression does." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAF3CAYAAAD+RdykAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3Xt4FOXZP/DvZIOLCYISTPQtJiGWMwSQg6IgAePhFV9btdbEQKUaUqUeqo1WTbW+reG1mnpswYZ4qmzBQz3U1l+rK0RRsQqFUoQGNSSBWjcQFQiHhWye3x+TWfYws8eZ3ZnZ7+e6uEImk9mZ3ezuvffzPPctCSFAREREROaSle4TICIiIqJwDNKIiIiITIhBGhEREZEJMUgjIiIiMiEGaUREREQmxCCNiIiIyIQYpBERERGZEIM0IiIiIhNikEZERERkQgzSiIiIiEwoO90noIchQ4aI4uJiw29n//79yM3NNfx2zCiTrx3g9Wfy9WfytQO8fl5/5l6/kde+fv363UKIE6PtZ4sgrbi4GOvWrTP8dpqbm1FWVmb47ZhRJl87wOvP5OvP5GsHeP28/sy9fiOvXZKk9lj243AnERERkQkxSCMiIiIyIQZpRERERCZkizlpao4cOYKdO3fi0KFDuh1z0KBB2Lp1q27Hs5JMvnbA/tffv39/DB06FP369Uv3qRARUR/bBmk7d+7Ecccdh+LiYkiSpMsx9+3bh+OOO06XY1lNJl87YO/rF0Kgq6sLO3fuxLBhw9J9OkRE1Me2w52HDh1CXl6ebgEakV1JkoS8vDxds85ERJQ82wZpABigEcWIzxUiIvOxdZBGsRswYEDUfe655x40NDRE3OeVV17Bli1b9DotIiKijMUgrY/LBRQXA1lZ8leXK91nlDyfz5fy22SQRkREpA8GaZADspoaoL0dEEL+WlOTfKD2u9/9DqWlpZgwYQLmz58PANi1axcuu+wyTJ06FVOnTsV7770HQM5SXX311SgrK0NJSQkeffRRAHJbirlz52LChAkYN24cnnvuOQDAW2+9hUmTJmH8+PG4+uqr4fV6AcjdF37+859jxowZeOGFF/DZZ5/hggsuwOTJkzFz5kz861//AgBs374d06dPx9SpU3HXXXdpXkN9fT1GjhyJiy++GC0tLf7ty5Ytw9SpUzFhwgRcdtllOHDgAN5//3388Y9/xK233oqJEyfis88+U92PiIjIdEKyNflud7rPSF7ZZfV/kydPFqG2bNkStk1LUZEQcngW/K+oKHi/vXv3xnzMzZs3ixEjRohdu3YJIYTo6uoSQghRWVkp1qxZI4QQor29XYwaNUoIIcTPfvYzMX36dHHo0CGxa9cuMXjwYHH48GHx4osviurqav9xv/76a3Hw4EExdOhQ0dLSIoQQYv78+eKhhx7qu5Yi8ctf/tK//5w5c8S2bduEEEJ88MEHYvbs2UIIIf7nf/5HPPPMM0IIIX7961+L3NzcsGtYt26dGDdunNi/f7/YuXOnOPXUU8UDDzwghBBi9+7d/v3q6urEo48+KoQQ4qqrrhIvvPCC/2da+1lNPI+9VUV6zqxevTp1J2IymXztQvD6ef2r030KqbF8uRA5OUFBQI/TKW83AIB1Iob4hpk0AB0d8W2PxapVq/Cd73wHQ4YMAQAMHjwYAOB2u3H99ddj4sSJuPjii7F3717s27cPADB37lw4nU4MGTIE+fn58Hg8GD9+PNxuN37yk59gzZo1GDRoEFpaWjBs2DCMGDECAHDVVVfhnXfe8d/2FVdcAQDo7u7G+++/j8svvxwTJ07ED37wA/znP/8BALz33nuorKwEAH+WL9SaNWtwySWXICcnBwMHDsTFF1/s/9nmzZsxc+ZMjB8/Hi6XCx9//LHqMWLdj4iIKG3q6oCQkR6H1ytvTyPb1kmLR2GhPMSptj1RQgjVFXO9vb1Yu3Ytjj322LCfOZ1O//8dDgd6enowYsQIrF+/Hq+//jruuOMOnHfeeUHBkprc3Fz/bR1//PHYuHGj6n6xrOjT2mfBggV45ZVXMGHCBDz99NNobm5Oaj8iIqK0MSJbowNm0gDU1wM5OcHbcnLk7Yk655xz8Pzzz6OrqwsA8OWXXwIAzjvvPPz617/276cVQCk+//xz5OTkYN68eaitrcXf//53jBo1Cm1tbfj0008BAM8++yxmzZoV9rsDBw7EsGHD8MILLwCQA8d//OMfAICzzjoLK1euBAC4NCbfnX322Xj55Zdx8OBB7Nu3D6+99pr/Z/v27cPJJ5+MI0eOBP3+cccd588MRtqPiIjINLSyMslka3TAIA1AVRXQ2AgUFQGSJH9tbJS3J2rs2LGoq6vDrFmzMGHCBNxyyy0AgEcffRTr1q1DaWkpxowZg8cffzzicf75z39i2rRpmDhxIurr6/HTn/4U/fv3x1NPPYXLL78c48ePR1ZWFq699lrV33e5XHjiiScwYcIEjB07Fq+++ioA4JFHHsFvfvMbTJ06FXv27FH93dNOOw1XXHEFJk6ciHnz5mHmzJn+n/3iF7/A6aefjnPPPRejRo3yb6+oqMADDzyASZMm4bPPPtPcj4iIyDRUsjU+pzO5bI0OJHn+mrVNmTJFrFu3Lmjb1q1bMXr0aF1vx86tgaLJ5GsHMuP6Iz1nmpubUVZWltoTMolMvnaA18/rz6Drd7nkOWgdHUBhIbbMm4cx995ryE1JkrReCDEl2n7MpBERERFVVQFtbUBvL9DWhs7y8nSfEYM0IiIiIjNikEZERERkQgzSiIiIiEyIQRoRERGRCbGYLREREWU8l8eDutZWdHi9KHQ6MQ9AWZrPiZk0SkhzczMuuuiiqPuVlZUhtDxKqIcffpiN14mIKG1cHg9qWlrQ7vVCAGj3etHQtz2dGKT1cXk8KF67FlnNzSheuzbtD0w6CCHQ29ub8ttlkEZEROlU19qKAyHvf96+7enEIA3qEXRNS0tSgVpbWxvGjRvn/76hoQH33HMPADm79JOf/ATTpk3DiBEjsGbNGgCAz+dDbW0txo8fj9LSUjz22GMAgPXr12PWrFmYPHkyzj//fH+TdK3jfPzxx/4uBaWlpfjkk08AAA8++CDGjRuHcePG4eGHH/af5+jRo7Fo0SKcdtpp2LFjB9544w1Mnz4dp512Gi6//HJ0d3cDAP7yl79g1KhRmDFjBl566SXV6z548CAqKipQWlqKK664AgcPHvT/7LrrrsOUKVMwduxY/OxnPwMgd2D4/PPPMXv2bMyePVtzPyIiIqN0eL1xbU8VBmlQj6AP9PYaGkH39PTgww8/xMMPP4z//d//BQA0NjZi+/bt2LBhAzZt2oSqqiocOXIEN9xwA1588UWsX78eV199Nerq6iIe5/HHH8dNN92EjRs3Yt26dRg6dCjWr1+Pp556Cn/729/wwQcfYNmyZdiwYQMAoKWlBd/73vewYcMG5Obm4t5774Xb7cbf//53TJkyBQ8++CAOHTqEhQsX4rXXXsOaNWvwxRdfqF7X0qVLkZOTg02bNqGurg7r16/3/6y+vh7r1q3Dpk2b8Pbbb2PTpk248cYb8V//9V9YvXo1Vq9erbkfERGRUQqdzri2pwqDNKQngr700ksBAJMnT0ZbWxsAwO1249prr0V2tryeY/DgwWhpacHmzZtx7rnnYuLEibj33nuxc+fOiMeZPn06Fi9ejF/+8pdob2/Hsccei3fffReXXHIJcnNzMWDAAFx66aX+zFtRURHOOOMMAMAHH3yALVu24KyzzsLEiRPxzDPPoL29Hdu2bcOwYcMwfPhwSJKEefPmqV7XO++84/9ZaWkpSktL/T97/vnncdppp2HSpEn4+OOPsWXLFtVjxLofERFlGJcLKC4GsrLkry6XLoetLylBTlZwSOTs255OXN0JOVJuVwnIkomgs7Ozg+Z3HTp0KOjnzr5jOxwO9PT0AJDnhEmSFLSfEAJjx47F2rVrVW9H7ThXXnklTj/9dPz5z3/G+eefj6amJkTq0Zqbmxt0e+eeey5WrFgRtM97770Xdm5a1Pbbvn07Ghoa8NFHH+GEE07AggULwu6TePYjIqIM43IBNTWAMoe5vV3+HpBbOiWhqqAAAIJXd3q9/u3pwkwa1CPonKyspCLogoICdHZ2oqurC16vF3/605+i/s55552Hxx9/3B9sffnllxg5ciR27drlD9KOHDmCjz/+OOJxWltbUVJSghtvvBEXX3wxNm3ahLPPPhuvvPIKDhw4gP379+Pll1/GzJkzw373jDPOwHvvvYdPP/0UAHDgwAFs27YNI0aMwPbt2/HZZ58BQFgQpzj77LPh6vtks3nzZv9Q5d69e5Gbm4tBgwbB4/Hg//2//+f/neOOOw779u2Luh8REWWwurqjAZriwAF5uw6qCgrQNn06esvK0DZ9OtLfuZOZNADqEXR9SUlSEXS/fv1w99134/TTT8ewYcMwatSoqL9TXV2Nbdu2obS0FP369cPChQtx/fXX48UXX8SNN96IPXv2oKenBz/60Y8wduxYzeM899xzWL58Ofr164eTTjoJd999NwYPHowFCxZg2rRp/tuaNGmSf4hUceKJJ+Lpp59GZWUlvH3ZxXvvvRezZ89GY2Mj5s6diyFDhmDGjBnYvHlz2G1fd911+P73v4/S0lJMnDjRf3sTJkzApEmTMHbsWJSUlOCss87y/05NTQ3++7//GyeffDJWr16tuR8REWWwjo74ttuAFGkYzPAbl6QnAVwEoFMIMa5v22AAzwEoBtAG4LtCiK8iHWfKlCkitBbX1q1bMXr0aF3Pd9++fTjuuON0PaZVZPK1A5lx/ZGeM83NzSgrK0vtCZlEJl87wOvn9Zvo+ouL5SHOUEVFQEjCQQ9GXrskSeuFEFOi7Zfu4c6nAVwQsu12AG8JIYYDeKvveyIiIspk9fVATk7wtpwcebtNpTVIE0K8A+DLkM3fAvBM3/+fAfDtlJ4UERERmU9VFdDYKGfOJEn+2tiY2KIBg1aJ6s2Mc9IKhBD/AQAhxH8kScpP9wkRERGRCVRVJb2S08hVonpL65w0AJAkqRjAnwLmpH0thDg+4OdfCSFOUPm9GgA1AFBQUDB55cqVQT8fNGgQvvnNb+p6rj6fDw6HQ9djWkUmXzuQGdf/6aefYs+ePao/6+7uxoABA1J8RuaQydcO8Pp5/fa7/jMqKtBfpaPQoYICfBAQSxh57bNnz45pTpoZM2keSZJO7suinQygU20nIUQjgEZAXjgQOrlv69atuk/0zoTJ41oy+dqBzLj+/v37Y9KkSao/M9Xk4RTL5GsHeP28fhtef6dqWIH+nZ1B12qGa0/3wgE1fwRwVd//rwLwahrPhYiIiOyksDC+7WmU1iBNkqQVANYCGClJ0k5Jkq4BcB+AcyVJ+gTAuX3fW9LOnTvxrW99C8OHD8epp56Km266CYcPHwYgR+gXXXRR2O/86U9/wqRJkzBhwgSMGTMGv/3tb8P28Xq9KC8vx8SJE/Hcc8+hurra3zpp8eLFCZ+v1nHa2tpw+umnJ3xcuygrK0NoqZdQWo9roI0bN+L111/X89RUxXK+REQZx0KrRNO9urNSCHGyEKKfEGKoEOIJIUSXEOIcIcTwvq+hqz+NofNKDyEELr30Unz729/GJ598gm3btqG7uzuoOXqoI0eOoKamBq+99hr+8Y9/YMOGDaqp1g0bNuDIkSPYuHEjrrjiCjQ1NWHMmDEAkgvS9DqOlSjdHVIpVUEaERGp0HOVqMHMONyZespKj/Z2QIijKz2SCNRWrVqF/v374/vf/z4AubfmQw89hCeffBIHQtta9Nm3bx96enqQl5cHQO7LOXLkyKB9Ojs7MW/ePGzcuBETJ07EZ5995s+Y3H777Th48CAmTpyIqpA/tueffx633HILAOCRRx5BSV/Lq88++wwzZswAgIjH8fl8WLhwIcaOHYvzzjsPBw8eDDv/BQsW4MUXX/R/r0y4VMb1v/Od72DUqFGoqqry9xL96KOPcOaZZ2LChAmYNm0a9u3bB5/Ph1tvvRVTp05FaWmpP5sY6Ti33347xowZg9LSUtTW1gIA2tvbcc4556C0tBTnnHMOOvqqUi9YsAC33HILZs+ejZ/85CfYv38/rr76akydOhWTJk3Cq6/KI+wHDx5ERUUFSktLsWDBAtVrBoC//OUvGDVqFGbMmIGXXnrJv/3DDz/EmWeeiUmTJuHMM89ES0sLDh8+jLvvvhvPPfecPxOqtl8oIQRuvfVWjBs3DuPHj8dzzz0X9T5RPPHEE7j55pv93y9btsz/t0BElJGqquQCuL298lcTBmgA5Bd/q/+bPHmyCLVly5awbZqKioSQw7Pgf0VFQbvt3bs35kM+8sgj4kc/+lHY9okTJ4p//OMfYvXq1WLu3LlhP7/mmmvEiSeeKCoqKsTy5cuFz+cL2yf0d2fNmiU++ugjIYQQubm5qufzn//8R0yZMkUIIcRll10mpkyZInbu3Cmefvppcfvtt0c8zvbt24XD4RAbNmwQQghx+eWXi2effTbsNq666irxwgsv+L9XjrF69WoxcOBAsWPHDuHz+cQZZ5wh1qxZI7xerxg2bJj48MMPhRBC7NmzRxw5ckT89re/Fb/4xS+EEEIcOnRITJ48WbS2tmoep6urS4wYMUL09vYKIYT46quvhBBCXHTRReLpp58WQgjxxBNPiG9961v+85w7d67o6ekRQghxxx13+K/nq6++EsOHDxfd3d3iV7/6lfj+978vhBDi/fffFw6Hw3//KA4ePCiGDh0qtm3bJnp7e8Xll1/uf2yU6xFCiDfffFNceumlQgghnnrqKfHDH/7Qfwyt/QK9+OKLory8XPT09IgvvvhCnHLKKeLzzz/XvE8CH8/u7m5RUlIiDh8+LIQQYvr06WLTpk1htxHpObN69WrNn9ldJl+7ELx+013/8uXye5MkyV+XLzf05kx3/Slk5LUDWCdiiG/MuLoz9QzoByaEgCRJMW9XNDU14Z///CfcbjcaGhrw5ptv4umnn074PBQnnXQSuru7sW/fPuzYsQNXXnkl3nnnHaxZswaXXnpp1N8vKirCxIkTAQCTJ08O6/kZzbRp0zB06FAAwMSJE9HW1oZBgwbh5JNPxtSpUwEAAwcOBAC88cYb2LRpkz8rt2fPHnzyySc45phjVI9zxhlnoH///qiursbcuXP9c8LWrl3rz2zNnz8ft912m/98Lr/8cn9JjTfeeAN//OMf0dDQAAA4dOgQOjo68M477+DGG28EAIwbNw6lpaVh1/Wvf/0Lw4YNw/DhwwEA8+bNQ2Njo/+8r7rqKnzyySeQJAlHjhxRvW9i2e/dd99FZWUlHA4HCgoKMGvWLHz00UcYOHCg6n2iZEcBIDc3F3PmzMGf/vQnjB49GkeOHMH48eMjPVxEZEYWqu9lFS6PR9e+3XrjcCdgyEqPsWPHhk3a3rt3L3bs2IFTTz014u+OHz8eN998M95880384Q9/SPgcQk2fPh1PPfUURo4ciZkzZ2LNmjVYu3ZtTE3MnU6n//8Oh0N1Lld2djZ6e3sByMGoskhC6/cjBbKPPfYYNm7ciI0bN2L79u0477zzNI+TnZ2NDz/8EJdddhleeeUVXHBBaKcxWeBt5ebmBt3eH/7wB//tdXR0+HtYRgqo1Y4b6K677sLs2bOxefNmvPbaazh06FDC+4kI9QxjeWyqq6vx9NNP46mnnvIPwRORxdTVHQ3QFAcOyNspbi6PBzUtLWj3eiEAtHu9qGlpgUulhlq6MEgDDFnpcc455+DAgQP43e9+B0Ce0/XjH/8YCxYsQE7obfXp7u5Gc3Oz//uNGzeiqKgortvt16+fZsbm7LPPRkNDA84++2xMmjQJq1evhtPpxKBBg+I6jpbi4mKsX78eAPDqq69G/f1Ro0bh888/x0cffQTg6Jy8888/H0uXLvX//rZt27B//37N43R3d2PPnj248MIL8fDDD2Pjxo0AgDPPPBNKkWOXyxWUXQp0/vnn47HHHvMHQhs2bAAg31+uvnmJW7ZswaZNm1SvYfv27fjss88AACtWrPD/bM+ePfjGN74BAEHZ0OOOOw779u2Lul+gs88+G8899xx8Ph927dqFd955B9OmTdO8T0Kdfvrp2LFjB37/+9+jsrIy5t8jIhMxYNQnk9W1tuJAX2JBcaC3F3WtrWk6o3AM0gBDVnpIkoSXX34ZL7zwAoYPH44RI0agf//+Qasm33rrLQwdOtT/b8OGDbj//vsxcuRITJw4ET/72c/iHuqsqalBaWlp2MIBAJg5cyZ27NiBs88+Gw6HA6eccopm4BLpOFoWLlyIt99+G9OmTcPf/va3oGyVmmOOOQbPPfccbrjhBkyYMAHnnnsuDh06hOrqaowZMwannXYaxo0bhx/84AcRV2Hu27cPF110EUpLSzFr1iw89NBDAIBHH30UTz31FEpLS/Hss8/ikUceUf39u+66C0eOHEFpaSnGjRuHu+66CwBw3XXXobu7G6WlpXj44YdVg6L+/fujsbERc+fOxYwZM4KC6ttuuw133HEHzjrrLPh8Pv/22bNnY8uWLf6FA1r7BbrkkktQWlqKCRMmYM6cObj//vtx0kknRbx/Q333u9/FWWedhRNOCGvgQURWYKH6XlbQ4fXGtT0d0t4WSg9TpkwRoUOLW7du9Q9Z6SUTqs5ryeRrB+xx/RdddBFuvvlmnHPOOao/j/ScMUPl7XTJ5GsHeP2muv7QOWmAPOpjYPkIU12/zorXrkW7SkBW5HSibfp0Q69dkqSY2kIxk0Zkc19//TVGjBiBY489VjNAIyILsFB9L1MKqYe6fMMG5GQFh0E5WVmo7ytRZQZc3Ulkc8cffzy2bduW7tMgIj1UVTEoS4TKytgZt96Kvz7wAOZNmmTa1Z0M0oiIiMjeNFbGzrj//rhLSqWSrYc77TDfjigV+FwhIluz6MpY2wZp/fv3R1dXF998iKIQQqCrqwv9+/dP96kQERnDoitjbTvcOXToUOzcuRO7du3S7ZiHDh3K2DeyTL52wP7X379/f3/XAiIi26mvV18Zm0Q91FSwbZDWr18/DBs2TNdjNjc3Y9KkSboe0yoy+doBXj8RkaUpiy3q6uQhzsJCOUAz+SIM2wZpRERERH4WXBlr2zlpRERERFbGII2IiIjIhBikERERUUZxeTwoXrsWWc3NKF67Fi6PJ92npIpz0oiIiChjuDwe1LS04EBvLwCg3etFTUsLAJiq2wDATBoRERFlkLrWVn+ApjjQ24u61tY0nZE2BmlERESUMTq83ri2pxODNCIiIsoYhU5nXNvTiUEaERERZYz6khLkZAWHPzlZWagvKUnTGWljkEZEREQZo6qgAI0jR6LI6YQEoMjpROPIkaZbNABwdScRERFlmKqCAlMGZaGYSSMiIiLbs0pttEDMpBEREZGtWak2WiBm0oiIiMjWrFQbLRCDNCIiIrItl8eDdgvVRgvEII2IiIhsSRnm1GLG2miBOCeNiIiIbMXl8aCutVUzgwaYtzZaIGbSiIiIyDaU7FmkAA2AXBvN7QaKi4GsLPmry5WSc4wVgzQiIjrK5Qp608p3u9N9RvYVcl+bLUAwtQj3ndoigVBFTqccoNXUAO3tgBDy15oaUz0ODNKIiEjmcoW9aY1saDDVm5ZtqNzXaQ8QrBI0Rrnvoi0G8A9z1tUBBw4E//DAAXm7STBIIyIimcqblsPrNdWblm2YLUDQCHxMmUmNct+FLgaodLuxvaICvjlzsKOiAn/dvFmujdbRoX58re1pwCCNiIhkFnjTsg2z3dcagU9JU1N6zieSKPddYAP1SrcbyxoaUOzxIEsIDPV4MOPWW+WgtLBQ/Tha29OAQRoREcks8KZlG2a7rzUCH2dnZ4pPJAZR7rvABuqLm5qQGzr8qWTd6uuBnJzgn+XkyNtNgkEaERHJVN60fE6nqd60bMNsAYJG4OPNz0/xicQgwn2n9Oecv3UrAKBIK8js6ACqqoDGRqCoCJAk+Wtjo7zdJBikERHFwyqTqxOh8qbVUltrqjct2zBbgKAR+LRWV6fnfCLRuO9c5eX+0hsCcn/OHVpBphKUVlUBbW1Ab6/81WR/6wzSiIhiZcYVeXoLedPqLC9P9xnZl5kCBI3Ax7SPf8B95/rb31BcUoJ5W7eGld64vboaB0K7CphsSDMSBmlERLEy24o8Ij2ZKWiMUbTCtSvKy7GwttY8Gcs4sS0UEVGszLYijyjDxVK49r25c4F7703RGemLmTQioliZbUUeUYaLuXCtRTFIIyKKldlW5BFluNDCtYpKtxs7KirQXVaGqtNPt+y8UQZpRESxMtuKPKIMF1i4VrHgrbfwuwcfxFCPB5LFF/gwSCMiiocFJ1cT2VVg4VoJcuP0x55+GtkHDwbvaNEFPgzSiIiIKDkprh+oFK3Nam5GXWsr6ktK0FtWhrbp0zHg3/9W/yULLvBhkEZEpBc7F7ol0pLi+oGBZTeUorU1LS1weTzyDjZa4MMgjYhID5lQ6JZITYrrB6qV3TjQ24u61lb5Gxst8GGQRkSkBxa6JSvRM+ub4vqBWmU3/NtttMCHxWyJiPTAQrdkFUrWV/lQoWR9gcQCmcJC+Rhq23Xg8nhQ19qKDq8Xgx0OSACE2mkEluOoqrJkUBaKmTQiIj3YaB4M2ZzeWV8DhxdD5591+XxQ6y9g9aK1WhikERHpwUbzYMjm9M76Gji8GEvbJweAxpEjUVVQkPTtmQ2DNCIiPdhoHgzZnBFZX4PqB0Zr+wQAvYAtAzSAQRoRkX5Y6JaswCJZX5fHE1OQotUayg64cICIiCiTKB8e6urkIc7CQjlAC/xQ4XIBdXWYpfVznSmLA9q9XjgA+ADNBQKB7DoXTcEgjYiIKNNEWv0YsPpTApJf/RmFsjhAmXvm69uuFqBJAAZnZ+PLnh4UOp2oLymx7VAnwCCNiIiIAkVa/aljkBaYPYvH7hkzdDsHs2OQRkREREeloOZfaPYsVnaef6aGCweIiIjoqBTU/IultEYou88/U8MgjYiI0odN6c0nBas/YymtAchz0ACgyOm0bS20SDjcSURE6aF3eyLSR8DqT9HRAUnn1Z1KaQ2fxs+f2kc0AAAgAElEQVSV1Z1FGbAwIBoGaURElB4pmqBOCehb/fl2czPKysoSPkw8pTVysrIyMlsWCYc7iYgoPdiU3hwMGnIO7LsJRC6tYefWTskwbSZNkqQ2APsgP649Qogp6T0jIiLSVWGhPMSptp1Sw4Ah50RKa9i5tVMyzJ5Jmy2EmMgAjYhIJ2aaqG+R9kS2FmnIOQGh2bNYZVppjViZNpNGREQ6M9tE/VjaE5GxdBpyTrQwLZCZpTViZeZMmgDwhiRJ6yVJqkn3yRARWZ7OWRNdsCl9eiVZE83l8WDImjWYt3VrXAFappfWiJUkRLT2pekhSdJ/CSE+lyQpH8CbAG4QQrwT8PMaADUAUFBQMHnlypWGn1N3dzcGDBhg+O2YUSZfO8Drz+Trt9O1z5ozB5LKa76QJLy9apXq79jp+hNh9+vPd7sxsqEBjoAAy+d0oqW2Fp3l5arX7wbQBMATw/Er3W4sbmpCYWcnduTn447qaqwqL0c1gHI9L8QARj72s2fPXh/LVC7TBmmBJEm6B0C3EKJB7edTpkwR69atM/w8mpNcimxlmXztAK8/k6/fVtdeXKw+Ub+oSM5iqbDV9ScgI67f5dIccg68fpfHg5u2bUOXT6vCWbBKtxvLGhqQG5hhy8kBGhstkTE18rGXJCmmIM2Uw52SJOVKknSc8n8A5wHYnN6zIiKyOCMm6ptpIQIlRmXI2eXxoHjtWswGkN3cDKm5GfO3bo05QAOAxU1NwQEakP7hdYsx68KBAgAvS5IEyOf4eyHEX9J7SkREFqf3RH2zLUSghEQrOBupvpmWnKwsFHV2qv+QdfBiZsogTQjRCmBCus+DiMh2+irJ64IdAywrMDBLNiALlZedjUeGD5fbSbEOXlJMOdxJREQWwI4BxtJpKPndJUuw86ST0JuVhfaTTsKVP/0p5gesxkx2ZnrgSs3lo0dj94wZ8mpN1sFLGoM0IiKKTCtYSLJ8A0WgDCW3twNCHB1KjjNQe3fJEky65RYM9XiQJQSKPB4sa2hAhduty2nmZWfj2dGjIcrK0DZ9enApjaoqeZFAUREgSfJXiywaMAsGaUREgUICknyd3sx0lcLJ+vlut3awYHSmJJMXJehU06745z8Pm7yf6/VicVNTQqelmTXTwjp4SWGQRkSkUMlejGxoMFdwoFOGJVYlTU2R550ZlSlJ8XWajk5Dyf+lMXm/UGtSfwAlIHP0fS1yOrWzZsnK5IA8AgZpREQKleyFw+s1V8mAFHcNcEZboWdUpsSM3RFSSaeh5M/z81W3d4Rs1wrIVgPoKSszJjBTZHpAHgGDNCIihRUmwqf4HL0ab/KGzzuzwmNhJJ2Gktvuvhv7Q5qX73c6cWd1ddDQpZIhMzwgU5PpAXkEDNKI9Ma0vXVZYSJ8is+xtbo6PSv0rPBYGCnJoWSlGO3ZY8bg5ttuw46CAvRKEtoLCrCwthbvz51r3NBlvDI9II/AlHXSiCyLxT2trb4++PGD3MfQYaaSASrnaGTQ1FlejjGjR+tXADdWKb5OU4qxpp1S86zD68VghwOHenuxP6Dl47I5c+AqL/c3Mv+9keecCNZT08RMGpGemLa3NpXsRUttrbkC7HSUNUjHCj2Wb4iJy+NBTUsL2r1eCABdPl9QgKY40NuLutZWA05Ah5ED1lPTxEwakZ6Ytre+kOxFZ3MzxqTxdFTp2TXAzDLlOpNQ19qKA729Me3bEdpHM1l6jRzo3a7MRhikEemJaXsiMlDo0GY8Dc8LQxYQJE3PtmAMyFUxSCPSE+fREFECQoMvSBK6enr8Dc/zVOaaxROg5WRlob6kRN+T5siB4TgnjUhPnEdDkdh95a/dr88gavPKunp6ABxteK411ywWednZ/kUDusr0FbgpwCCNSG9sg2I/egQfdi/YaffrA4L/DoYMkf/pEJDWtbbiW2+8ge0VFfDNmYPtFRWoTKIdWV52NiTE0bopUZzwbzgGaURkPzpmdCL2royH3Vf+puP6Yn2cjQiyu7rkfzoEpGf9+c9Y1tCA4r4m6MUeD5bX18M3e3bcAVuR04ndM2agNxX1zzhyYDgGaUSUOqkYDtM5oxOxd2U87D5/J9XXF+vjrNffg1oQGiiOvwml0KzU3Izs5mbUNzWFNUHP6vtX7PFgWUODaqBW6XYHZd8WvPWW/vPOouHIgaEYpBFRaqRqOEznjE7U3pWxsvv8nVRfX6yPs15/D7E83hH2CQzMXn/sMTRfcgl8c+bg04oKFHo8EQ+b6/VicVNT0LZKtzss+7bsV79CVRLDpGQ+DNKIKDVSNRymc0ZHt96Vdp+/k6rrU7KxaqVugPDHWa+/h1geb419AhcGVLrdaAwJrmK6+c7OoLlmjc88E5Z9yz540D7D5wSAQRqR9VllRV2qhsN0zujo1rvS7vN3UnF9gdlYLaGPs15/D2pBaCCNvwmXx4Ortm71F5xdrDG0Ga0cbdbgwdg9bx5658xBW2UlBuzcqb6jXYbPCQCDNCJrs9KKulQNh+mc0eksL9cv+LD7/B2jry/avDC1x1mvv4fQIDQvT/4X4W9CyaAFVjMr1Bo+B9BWUCAHa5IU/IN+/YB9+4Kf56H7+G/AJsPnBIBBGpG1WWnFYKqGw4zI6Fg1uLJKljWU1nlHyhJpPc56/j0E/h3s3i3/U/mbUOafzQvIoCk6NIbPOwoKMPb557Hiiy+AZ58NPt+BA4HDh4N/QYjwQM1Ow+eBrPp3rAN2HCCyMiutGExlfz62mNGvr2KquFzy34aSJVIKtwaet1bbtaIiOVDSYuDfQ2ingNCuAKHurK7GsoaGoCHP/U4nHrz22qMFZ0PPN0sjnyKEfO127ndptb9jnTGTRhQrM36as9qKQatmpKzISlnW0LlmoUGOct5pXnyhZMiympsxZM0aDHj7bczbujWoU0C0rgArysvxg9padA8d6s+U5T7xBB695x7tmmZaz2clOA18PpnxdSoZVvo7NgCDNKIY6FbQVG92XzFoBLu9iSlCryvW1Y9mEG2uGSCfdxoXX6i1bkqkTVNOVhb++4YbMGDHjtg/rMT6PLfSHNVYWWm0wAAM0ohioFtBU73ZfcWg3uz4JgaoX5eVJpbH8oarnHeasrF1ra1h88viVeR0JtZDM9bnuR2zTlYbLdAZgzSiGOhW0NQIHEKMnR3fxAD167LSxPJob7gmOO+OkLIZ8cjJysLy0aOTa9MUy/M8WtbJilnkDB8tYJBGFAPdCppSetl16ETr/JWJ5QDgcBwNSM325qz2RqwEmKnMDgcEMd2nnIIb77kHswFkNzcj/oFNWV52dmLZs0REyjpZNYuc4aMFDNKIYqBbQVNKL7sOnUSaWK4EQL6+al1mfHNWeyN+9lk5mEhBdtjl8eDGe+7B/muu8QcxA3buxP/ddx8q3e6gOmeR5EpSUFeA5aNHY/eMGakJ0IDIWScrZ5EzeLSAJTiIYtBZXo4xo0enpnwEGae+Png5P2CPYDvSdUV6czbT328Ky2RAktDV0wMHAB8ACUDr44+HdQJQemauKC8PO2Ze33G+7OlBodOJ+pKS1AVjWiKVuZk/X/13rJ5Ftjlm0ohilcGf5mzDrkMnka7LrkO8oVTmW7k8HgxZsyasTEZXTw8A+DNkAtqdANS2SwB2z5yJ3TNmoLesLLm5ZnrTep3SM4ts5Nw2K86bMxCDNCIyN71etJXjKBmFZ5+1V7Cdijdnk1HqllX99Kc4EDBUifZ27L/mGrz+2GPo8sU2WKnZCUBle6HTmdR5p4VeE/CNnNtm1XlzBmKQRkTmpdeLdia/+Nt0dVxg3bL6pibkqAxV1jc1xXy8O6ursT8k+NrvdOLO6uqgbTlZWagvKYlycibMBumVRTZybpuV580ZhEEaEZmXXi/amfzib7MhXrW+mPEMVWpZUV6OhbW1cpNzSUJbQQFqamuxorwcjr59YqpzZuYPBHpM2TBy+DxThubjwCCNiMxLrxftTH/xt8l8ysDsWaB4hiq1SJADtW+uXAnHqlUoe/llXHjDDVgNoKesDCLWuWd2/0AQy/B5oplEGw/NJ4pBGhGZl14v2nzxtwWtqv+xDlUCR8tkAAjKkD07ejREWVl8AZkau38giDZ8nkwm0aZD88lgCQ4iMi+9SmbYtfSGzYWWztBaBKCUyFjc1ITCzk505OfjzupqrCgvhwR59WZRqspkFBaq9021yweCSGU+lO2JlnyJduwMxCCNiGQul/leHPV60eaLv+UoQ5tK5izaKs03zj8fb1xwQVD9s5QFZoEy4QNBpJp2yWYSDayXZ0UM0ojo6BCF8saiDFEA6X/B1OtFmy/+lhJrQ/OcrKzUtV2KRaZ/ILB7JjHFOCeNiKw32dmMJQ4oIcpqTam5GdnNzZCamzFkzZqwxQFaTBWgKWyyUCNMLM87zivTFTNpRGStyc5mzvpRRMocs3avN6glk9K8XBnQjLUAbZHTab4Aza5ifd5leiZRZ8ykEZG1Vj9aLetnJJNmFJXsWFZfVmzIu+9Cam7G/L72TEBwS6ZExFRUNhVM+hjoLp7nXToziTZ7PBikEZG1hiislPUzUjqKpsbwBhhYyyy0V2aiAZkiLzsbEmIsKpsKZi5cqzcrPO9s+HgwSCMia1Wlt1LWz0ipzijG8Abo8nhwVUAnAD0VOZ3ma2ieSVldKzzvbPh4MEgjyhTRsiBWmexspayfkVKd2dB6A5w3DyguxrtLlqCmpQWxzSaLj2mGNkNZIbukFzM/75TXNrVVpYClHw8GaUSZINFhgL4Xv1lz5phnfoeVsn5GSnVmI9IbXXs7TrvlFnzrjTcSOrTU91XpAJDncJhvaFONFbJLeknH8y6W+WWBr21aLPx4cHUnUSZIpAp4wGouCTDXKkrWPEt90VSt+lfKTXu9WNzU5K/+r0VZzZnWgrN6yYTCtYFS+byLdTWp2mtbIIs/HsykEWWCRIZlbDi/w1aMzmyEZlEvvBA9xx4b8VcKOzvDtkkInvCvW49MM2BW1zixvv5Eeg2zweMRNZMmSdL1AFxCiK9ScD5EZIREqoBn0nwbqzIgs+HyePC3pUvxf/fdh1yv159F3f/kk3jqggtw0QcfoMjj8Q9RBurIzw/63nTdAIzArK4xYn390XptKyqS59ZaXCzDnScB+EiSpL8DeBLAX4UQya6kJqJUSmRYhu1dbCm0aTkkKajfpQSg9fHHkRtS8T/X68VFH3yAYStXotLtxrKGhqB99juduLO62v+9pYcxKf1iff2x+ZBz1OFOIcRPAQwH8ASABQA+kSRpsSRJpxp8bkSkl0SGZcy8msvudC7IGdh6SSkoG1rDLLC4rNqwJQK2rygvx8LaWrQVFKBXktBWUICFtbX++WhFTmf6hjFtVsw03fLd7vTcn7G+/th8yDmmhQNCCCFJ0hcAvgDQA+AEAC9KkvSmEOI2I0+QiHQS77BMQHsX0dEByaj2Li4XW8gESrDtlVrLpTyHA4d6e7E/YPAjlmGQjvx8FHs8qtsVK8rLVRcJ5GRlYfmGDUBlZeofU7YM05fLhZENDYCSMU3l/RlPeykbDzlHzaRJknSjJEnrAdwP4D0A44UQ1wGYDOAyg8+PyNzs/qm9r3ba26tWGVM7zYYVwpMWx4INtQwZENwDc38Cs1PurK7GfqczaFvocKaavOxs/HXzZsy49db0PKZc7KKvujo4Qhvdp/L+tErtRgPFsrpzCIBLhRDnCyFeEEIcAQAhRC+Aiww9OyIzs0KAYfYgkm+q4TQmTIuOjqB+mAPefhvzAgIzPScKRxvOBIBcSQpatbl89GjsnjEDM+6/P32PKRe76Iv3Z9rFMiftbiGEanEcIcRW/U+JyCLMHmBYIYjkm0A4jYUZO/Lzg+aSJZIhi5UEOVD75sqVcKxahSkvvIA3LrggKCDrnjVLvU1TOh9TMxeXNfsHJjVG3Z9WvC/ShHXSiBJl9gDD7EEkYO431XRRmTC93+nE7VGGGuOhlM9QKvsDR6v9h9YyWw1g98yZsffN1HrshDD+Ddmsi12s8IFJTX09fCHD3knfn1a9L9KEQRpRosweYJg9iATM+6aaTlVVePeBB7AzwlBjMvKys/1BmBJ86VpcVu0xVRj9hmzWlX5W+MCkpqoKLbW1+t6fVr0v0oRBGlGi0hlgxDJcYPYgEjDvm6rRAh6/7lNOwY333BM012zmmDE4pW+ocdjKlTEHaNF6YCrzxgwtjRH4mKox+g3ZjJPNrfCBSUNnebm+96eF74t0YO9OokTFs0Q8GaElKi68EHjmmehlBqxS5NHGy+eB8OKxl775Jh66/35/IdgBO3fi/+67D7t7epLKluVlZ+OR4cPNUTxWeUyzsuQhrVCZ9obMwtBH8b6ICzNpRMkw+lO72vyNxx+PbbggU7NUJuHyeDBkzRr/Ckxlwv+djY2q1fwXNzXFfOzQfpgpyZAlwgrZ3FTgsP5RWsPh3d2cl6aCQRqRmanN39Ba1aeWnTDj0I9RTLBiLLRuWZfPF7ZPtGr+0eRkZeHZvqAspon86cTgRMYPTEcp90VeXvD2ri4uIFDBII3IzOIZFsq07ESgNKwYUwKyeOuWhTYhj7YdODrXrMjptFbDcisHJ31B/6w5c/QJ+jPpA1M0VVXAgAHh27mAIAyDNCIz0wq8JCn4+0zMTgRK8Yoxl8eDmpaWhOqWxVvNP3A1pqmzZlqsGJwEBP0Sy0ToIzTTrTYvDci8+YpRMEgjMjOt4aJrr7VmdsIoKV4xVtfaigO9vQn9rlo1/5q+EhtpWY1J4VgmQl9qme7QD5qKTB4RUGHa1Z2SJF0A4BHIq8mbhBD3pfmUiFIvVStIrS7BFWNaTckhSejq6Tm6bc0aQJLwZU8PBjscqnPNYpErSejvcGBFeTmeLy+HD3IwVl9SAhcDMfPQCu7b2+WAg8+/+GjNrZWk4Dm2mT4ioMKUQZokSQ4AvwFwLoCdAD6SJOmPQogt6T0zojSweYkKXSRQbkQZslQyYoFNyRVq2+IJ0CTIc9OUQIwZMYvQCvoB9XI3FJlW0CuEPBLAD6CazDrcOQ3Ap0KIViHEYQArAXwrzedERGYVxwR1ZcL/vK1bEx6yjMXCVauwd/58iDlz0FZZiSq327DbIp1F6prAYc/4aWW0i4qsN18xxcwapH0DwI6A73f2bSOidItU6iKdZTBimKAeOOFfb4FzydZs2YLGX/0KA3buZH9CK1KCfi2c3B4flmJJmCRiXJGUSpIkXQ7gfCFEdd/38wFME0LcELBPDYAaACgoKJi8cuVKw8+ru7sbA9SWDWeATL52IPnrz3e7UdLUBGdnJ7z5+WitrpbbrViEcv35bjdGNjTAERDk+JxOub8foPkzs1xrBQCPAcctgJzuV5xRUYH+nvBbOlRQgA9S8Fqlp0x+7tvpcUyUXo+/FV8Djfzbnz179nohxJRo+5k1SJsO4B4hxPl9398BAEKI/1Pbf8qUKWLdunWGn1dzczPKysoMvx0zyuRrB5K8fmVlU+h8KQutyPRfv9bSeaVPo9bP2toMPLvYuDwezNu6Vffj5mRlhdcu02qHJElyps9CbPfcD22zFmkelMsF3zXXBH3wsNpzN1m2e/zjYOS1S5IUU5Bm1uHOjwAMlyRpmCRJx0D+APzHNJ8TUWLstJw/UqmLdDVOVhliDaz8n93c7O8AEEloU3K1bRLCm5arFpdlOyRzirfocVWVnCWOtdyNCbpekL2YcnWnEKJHkqTrAfwV8uvkk0KIj9N8WkSJSVfwYgStVW9ZWcDxx8utXdR+xyguF3oWLkT2wYPy9+3t2H/NNXi9thbtfUMpylpMrTED1UxYgObmZpTNnBnfeVmluX2mifSBSSPw6iwvx5h7741+7NCMuRIAAhmTdSP9mTWTBiHE60KIEUKIU4UQfGUj67JTVkVr1ZvPB+zdCxxzTPB2HQOT0DZMQ959F20//vHRAK1PrteL+jialRvSZsnK7ZDszMgPTHbKmJNpmDZII7INO61sUoIPhyP8Z0eOAMcdp0tgEqkvptKGqaunJ+lm5UVOp3G1y6zUDilThumM/MBkp4w5mQaDNCKj2S2rUlWlPfn9yy+TCkxcHg+GrFkTFpBp9cVMpFm5IicrC/UlJXGdny2loTl92hj5gclOGXMyDQZpRKmQiqxKKrMhOr4hBU7yn791a1wV/eNtVq50C9Sc8J+JMmmYzsgPTHbKmJNpmHLhABHFKdWTlhOcGK/0yuzwejHY4cCh3t6gLFm8BYFW9C0OWNzUhMLOTnTk5+PO6mqsKC/3t2RS+m/G3Zqpr1TDLLu3rMm0YTqj2qyxzy4ZgEEakR0ksGoNQHw1owJpvCG5ystRt3atPwgLbEoeGpAl2qQ81IrycqwsL08uIAsVEPRKgL1X6iXYnJ5UsM8u6YxBGpEdJJINSTL75iovR11JSXBWLKAWWaJNyWORK0no73Dgy54eFBrRvDzRoNeKWC6EyLQ4J41ID2rzwcw+RyyJuUiBPTCjTe5PVq4kBRWPXT56NLpnzcLuGTPQW1aGtunT9Z9blklDgHZb2JKJMmV1bgZiJo0oWWoZqe9/X37DO3z46DazzRFLIhCpa23FAYPaGylzyZIeskxGpg0BWmGYLtGhebtjEV1bYyaNKFlqGakjR44GaAojV8yFZkPy8oBjjwXmz9f+ZJ3gCk2Xx4P2wF6GOsrLzsazo0dDqGXIUpkt4Eo9c8mkMiHxyqTVuRmIQRpRsuIZAjNyuEwp8/Hss8DBg3KLpkhvaHEEIm4gqEyGHtSGMXfPmKGeOUv1m3RA0Cs4BJh+RgciRn4AMPrDRSYNzWcgBmlEyYpnCCwVw2WxvqFpzEVylZeHVfuvB/zZs1hnngUGYaFNyeOeV5aObEFf0Pv2qlXm7xgQyI7zk4wIRJT7SZLkjLMRHwBS8eGCRXRtjXPSiJKlNh+sX7/gOWlA6obLYnhDU+qVtX/jG3A8/TR8kAOpSCs0o8nLzo5vtWU8c4yYLYiNXecn6T1HMPR+Cl30otdK3lSsEubqXFtjJo0oWWoZqaeeAp58MqUr5pTK/W0aLZF25OdjyLvv+ocslcyYEoYls0KzyOmMnBULze4sWhRfhoHZgtjYdX6S3nME1e6nUHp8AEjFhwuuzrU1BmlEelBr+5SiBtuh/S61WiX9pLoaXT09AOKv7B9J1B6YakM+jz8eXzCR6Ju0lecaJcKuGUe9A5FY7g89PgCk6sNFil5rKPUYpBFZmFKvLHBYckV5ORbW1qKtoAC9koS2ggIsrK31t1CKptLtxvaKCvjmzMH2igpUut1h+8TVA1Mta6GVsdN680zkTdrI+UBmXW1o54yjnoFItPtDr+FCrhKmJHFOGpHJhfa7hCShq6fH3wJJzYry8piDskCVbjeWNTQgt28otNjjwbKGBv8xAXnu2SPDh8devyyeLE6kN894a3kZOR/IrB0J4p2flKm1x9TuJ0mSA+6iIv3uB/bzpCQxSCMyMSVTphSODcyY6dtoSba4qckfoClyvV4sbmrC+3PnJlZcVmvSt/KmqNA7w2Dk0J9ZhxXjCQrsusggFqkMnqxQKJhMi8OdRCaiTP5Xyl98b+tWQyv7A3JTckBe3VnY2am6b1FnZ+Ltl7SGfK691tjJzkYO/Zl5WDHWYUG7LjKIld3ncZlxziTFjUEa2Z9JX6wCA7Kau+9Gx0knofLkk9F8ySWocLvR5fMh0fBMqVEGBAdhSq2yG5qbsXf+fIg5c9CzYAHEv/+N3TNnIksjyPBqrBiNidZ8siVLjH2TNHI+kB3mGqUzG2jS56RtmHXOJMWNw51kbyYd0gkcxqx0u/FQlHlgsYqp36XLBTzwgPp9ojGnqbW6GmPiOpMQ6RjyMXJIyw5zjdLVn9Skz0lbMeucSYobM2lkbyYd0glsUK41D+yRxx6LuspSkZOVheWjR8c2JBntBVwl69WZwCIEUzBySMvqw2Xpygaa9DlpK2adM0lxYyaN7M0kL1b+Cv9eb9iqTK15YEP27sWJe/cCOJpdkwD89YILglZ3xpQ9CxTtPlHLejU3x3Zsso50ZQNN8py0tXRlSUl3DNLI3tLwYqUWkEk4WkA2dFVmR34+ij2esONIId/ner347TPPYMC99yZ3gnwBJ0U6hqH592c8toqyDQ53knXFMvnY4CEdZfK/1NyM7OZmzZZLkSr8q3UI0Np/wL//nfQ5G3KfWHUiuFXP28rssOjC7NgqyjYYpJE1aa1eWrQo+E0XMOzFSpn8H09ApmZFeTluvu027OjrELCzoADeE05Q31mPbIPeL+BWXUlm1fO2OgYQqWH1OZMEgMOdZFVak48ff/xogVTlTbexUX6RSlKkeWXJKHI60fjznwM//zkAYCgQvgIO0DfboOcwl1VXkln1vO2ABV6JYsJMGlmT1iTj0J6Q8awaUxn6UoYzZwOqw5jJ0mxObqVsg1Unglv1vIkoYzCTRtakNflYTSxvuiq1m3oWLsRff/RjtJ93DoD4hzHD9EJeDdALIAvI2uWE9PsSzP9jAW4aLO/S1QU4HIDPB+TlVQGoQpcAHDsB3zwg7ya1/eLb9uWXwGDV24u0bZbmfp+KQhQj/LHYIRVi0pBEb8/Ia4ntvI/+7ixdbs/IazH29maZ/vqMve/C//atey2J3J72c9961xL7cQoLgXnz8lFWhrSSRGjmwYKmTJki1q1bZ/jtNDc3oyzdj1iamO7a1YYDQ3tBKoqKog93FherBn09WVn43h13xN+sXAnIfJBL/nucQFMJ8FYCbZVMrhIuLEMNcnH0sdiPHLx75nSc891mOIb44NvtwOPP1+CG95ek8UyDaZ33QjRiBUyYsSSilHI6fXjiCYchAxiSJK0XQkyJth8zaWRNajWeLrwQeOaZxOZxaWTbsnt7Y6v+L5ARAZkaJaBZjDoUogMdKETLmd/EedVvQepbtJp9og8/rF4KAKYJ1NTO+07UM0AjIgCA10vZO/QAAB83SURBVOtI+xRVZtLiYLpsUgpZ5tpdrtiKc4bs1/3lXgzY95XmYdsKCjBs5Ur1Hx7KAhpGZkxQFosjD2cj+8TwmXs9uxzo96OeNJwREVH8JEleIKv/cWPLpHHhANlLLMvOVUovZB/YC2+2dmLZ3xWgF3LWrKfv6xdOBmgqHEPUl1Y48mJbcvHYmYtw5OFs9C6XcOThbDx25iI9T4+IKCbprrHM4U7KPCqlF/r7fNg1cCBO6O5GtsrHpo4T8+WALIOGMZPh2+1QzaT5uhxRf/exMxfhh9VLTT1USkT253T6UF8f/TXLSMykkTmkovJ7320IjVWhefv24Xt33BFW/X+/04k7h9yNvOunI2+jHKA5+p63eXlHVwQluk2S9DmOsbcn4rq9x5+vgQjuGQ/hlbdHO861VzT6AzSF5ASu/W5jmu47kZLbM+/fgTD99Rl734X/7Vv3WhK5vfie++a+ltiPU1QE1Na2pL3qETNplH4q5S9QUyP/P8FnyKKXPGjsbYVvsBfY58CV77yJxt/cj1yvN6wnpqJXkn+ysLYWi5uaUNjZiY4h+Vh15d2o+fYY/L4soVOxhebmt+Ock7gEcAH4qhEY5AP2OCAdX4Pr31uC66P9qkt9SDQ7z4eeNExni//a7YXXz+vP1Otvbu4EMCat58AgjdJPx8rvi17yYKljG3CC72iH8kE+1Lsakev1RvxdZSXnwh/X4pu/fhk1WSVYcmkBroa8cILiVLUEQALDk1875Mcv1J70DjsQEaUagzRKvwQrvwdlywTkwfvjoTqI75/4H0Wu14vHGp/B73fdG9P+ZIDjawDvUiBwyNPbt52IKIMwSKP00+oeMHiwPIesowMdefm4Y2E1Vpx+PiRJghjQox6QaYxlduTno9jjiel08rr+Hdfpk86qwodKcXxNX2bOhlyLgK8bgeN9chbRztdKRHHhwgGKX+Ak/yFD5H/JTPivr5eLzgbq1w/Ytw9ob4ckBIp2e7DswQZUfvRXiIE98l9uHH+9d1ZXhy0I0Cx9k+411yQHKT/sAeYJ+atdgxbXIuDw0qPD8yf45O9dLDlCRAzSKF6hNca6uuR/ffXGUFMTf6Cm1kx84EDg8OGg3XK9Xixuagr79Uq3G9srKuCbMwfbKypQ6XaH7bOivBwLa2vRVlCAXklCW0EBlsz9NrzOY4N3jLVDAZEevm4MHtYF5O+/bkzH2RCRyXC4k+KjNsk/UIIT/lFVFfw7WeqfH0LnllW63VjW0OBfFFDs8QS3cRIADkqQehxYMaccK2aXA1mAo8uJmqwSOCvdsXUoIEpUpOHM4zWK+w6KregvEdkbgzSKT5TJ/DHvE43GPLWO/Pyg7xc3NYWt2sz1evHIo4/h/37bhMKuTkhK8HWxWvBVxaCMjOMfzuz7/gSfvCjCBTlQ40pWIoqAw50Un1jmayUxp2vRSx5kv7gWV145T72obHV18E1prNo8cd9eFO32QEpmGJYoVq5FwG+yAZfU97VvTlm04czja+SVq4HsvJJV634iIlUM0ig+apP8A+XkABdeGFP3ACUgk1Y1Q3p1DaTX38bSE7bCN8SLFeeFzyFbWFuLFbPLIe3NBnoBx24nvsr7RmznrQzDphvfpMwpmccl0uT/aMOZVUuAY64DvnLIK1m+csjfqy2UsPrfDhdJEMWNw51m4XJZY26Uck7KuQ4eLH//5ZfyeV94IfDMM1G7Byx6yYOlOS1A/741lipzcFaUl8vzygL1Ar0Xzzj6vfe+4G4FkegxDJuMaENflB7JPi5fNx79XYUTcgkRxDCcGUvRXzv87US8nyxyDUQpxkyaGYSumDT78FxVFdDWBvT2Art3y/96e+Vtr7+u3T2gz6KXPFg6aOvRAC0Oji9Dxo7UVoYqDdlCpbu0BlfymVOyj0ukbJlew5l2+NvhIgmiuDFIM4NIbZHilYpG5ZFoZKtER4d/aHPp8VuBROZFH8pCTVZJ+PbAoLGtDXjkkfAhWTOU1uCblDkl+7h8rfHHvMcR33CmkedoBpHuJyJSxeFOM0iwLVIYAxqVRxPammn7ifko7gyv7N9+Yj58QyL3zlTVC0A6WjJjyaUF0X8ndEjWLMPHXMlnTsk+LtHaWCXaw1TPczQDtvsiihszaWagNQwX6/Cckj2bN0+/jFwMlHllviFe+S/JAdy5MLyy/36nE3curFY/iBYBSHuzcd3XoyFml6HnO9NjC9AUodm1dAdoQOat5LOKZB8XvbJlRp6jGaTifiKyGWbSzKC+Pnzyu9bwXOgCg9CJ+mqSmTDfd3uzAjJSi44tl7Nned6wXpnKRP/FTU0o7OxER34+7qyuDl8AEMoHSPuzIQb0wPFlHFkzK8m0npRWocfjoke2DNAufGuXvx297ieiDMEgzQxiHZ5TG858/HF5sUEkiU6Yd7ngvWYhnN6DcizW3o7911yDr39cC9+52kGX6qrMSA5l4boDI+WgzCqrXBPFNylzMsPjEm0FpxnOkYhSisOdZhHL8JzaAoNoAVqkjFyUBQZdP7odTu/BoG25Xi8WPxHePzNuvQCEXOssKEBTW+W6aFF6F0MYxep1r0hfdljBaTZ8jpHFMZNmJfEOWxYVxZ6RU1lgcELXv1UPq1XlPyYCkPZl49qe4eHDmVqrXAOzhSlYDJESdqh7RfqywwpOM+FzjGyAmTQr0Rq2lEImhuXkAMuXx5eRO3AAuOkmf8aq68RT0HXccao3F9o/E4DcyNzX93WPw98VIPD/jt1OXPfVaPRePEN9vplWEBqaLTxwQF4kYeWsmsWyJt/88OHUZCQyOfOhVaLiADL3PkmGxZ5jRGqYSbMSrQUGV10lF5GNdR6XVi2zri5IXV0AgLzdO5URyaC1AWr9M4PmlCVDo6m6Jitn1ayUNXEtwjdOfPXoG55RGYlMz3yolag4Avn7AX1/F5l2nyTDSs8xIg3MpFmJWnX9xkZgyZK4yk10afS7DMnHIStgmwCwa+BAuX+msihAKZMRS4AWS5Fdtb6goVnCUGbpyRkvKxX2/LoRUioyEpme+VArUeGVgH4h+2XSfZIMKz3HiDQwSLOaZOt/uVzoPbQXocsNoiw/gARgf/9jsXLaBbENXQbcHoYMkYcnAxcEzJ8vB2CBAZtaEHrttZEbugPp78mZCCvVvUpVRoKZDzlQ+2EPME/IX3M1npmZdJ8kykrPMSINHO7MAEpXgO/+/c9Y9mADTvQefeUSAHYPHAgAOHHv3ojHKdzVGdzcPJrQBQqBtBYCKP8CnXWWnC3TGgpNd0/ORFip7lWqqt3boaq+3nifJM5KzzEiDcyk2dSilzxBvTJ9Q7xY/GQTcr3BHy0lAPuPPRY33XBDWKeAUF9pDJNqUlugoCbakKWSPVy+3Jw9ORMVmjWpWmLOifPH10CkIiPBzEc43ifJUXuOEVkIgzSbWfSSB9Kra7D0hK3+dk2Vq9zYXlGBIk94T01ALqmxorwcC2tr0VZQgF7IU2ICHerXH3kP3xffycQzDBnLvlpz8qy2aECLf+K8T46eT/DJ36c7UKtagn/v+pbx7XzYNigc7xOijMbhThtRemmi/9EQq9LtxrKGhrAMWqCOE/Mh7c3GijnlWDG7HMgCql5pxiO/fwJ5Xf8GCgvROm8exsQbDMWzWjMrSx4ejXYbasOhdvF149GVjQon5OGaNFea/3TajzC07BXjb4hV9cPxPiHKWAzSLEyZa+Yb7JUnl52AsCWai5vChzgD7Xc68dOKa8Pnms0uAx65x/9tZ3MzxsR7gmolQ7T4fNYtp6EXTpwnIqIAHO60GLW5ZsgC4EB4DQ1odwcQANoKCrDwxlsxcNZ1yZ+YWokNteHJ5cvlfw6Vic9WLaehF5YMICKiAKbLpEmSdA+AhQB29W26UwjxevrOyBwWveTBUse2o/OVoqh0u7G4qQmSRm/P9vwCfPPXL6MmqyT5IrTR2kypZcbmz1c/lhXLaehFrZgpJ4kTEWUss2bSHhJCTOz7Z50ALZaCrfEcR5KA7GwIScJt101A5Ud/jTlAW9bQgGKPR333nBwUP/gr9HxnenwBmtb1abWZipQV0yqbYcVyGnrhJHEiIgpgukyaZcXYtDzu4/h8kAAUd3qwvL4ey+vr0VFQgDurq49W/g8RcR6aVtP1KPLdbuChh9SvTyv7FSkrptXiyqrlNPTCSeJERNTHrJm06yVJ2iRJ0pOSJIWudzMXJbs0b1782SQ1EWqLZfX9K/Z4sKyhAZVut+p+WvPQIEmJdSkAUNLUpH19iWTF7F5Og4iIKEmS0JizZOiNSpIbwEkqP6oD8AGA3ZDntv8CwMlCiKtVjlEDoAYACgoKJq9cudK4E+7T3d2NAQMG+L/Pd7sxsqEBjgirJ4Uk4e1Vq2K+jVlz5mjOIwvVVlCAYb9fKQ+B+gA4gKzdx+A/P7gM+V99Ebb/oYICfJDg/aR1XkKSsPXOO8PuB5/TiZbaWnRqZPusJvSxzzSZfP2ZfO0Ar5/Xn7nXb+S1z549e70QYkrUHYUQpv0HoBjA5mj7TZ48WaTC6tWrgzcUFQkhNzjS/ldUFN+NxHLMvn8+QOweMlQISZJ/b/ly+RjLlwuRkxO8f07O0Z8n4GBBQeTrW75c/n/oudhE2GOfYTL5+jP52oXg9fP6V6f7FNLGyGsHsE7EEAeZbrhTkqSTA769BMDmdJ1LVNFWIiYyx6q+PnpD8T5ZkoS83TuPNi2vqdEufZHkUGJrdXXklkzJNn4nIiKiIKYL0gDcL0nSPyVJ2gRgNoCb031CmiLNuVICIyC+FZ99AVb7ELk9U09WlmqbJkjS0SblisA5cDoHTZ3l5ZxDRkRElEKmC9KEEPOFEOOFEKVCiIuFEP9J9zlpUst65eTIxVrb2uTva2rkLFdotiuCRceWo7jpBThWrUa/t96CY/VqzKurQ1t+AYQSIGnNWzOyzhizZURERCljuiDNUqINK8ZRPyy0kwAGBRStFcDKaRfg/qX/gKQESEVF6ueUyXXGjOBaBPwmG7P+PRv4TXb6m50TEVHGYJCWrEjZpRjrhymN0f0tnkIfFQnIOuwILjyrlcXL9DpjenItAg4vBU7wQZIgd3s4vJSBGhERpQSDNCPFWD+ssbcV6B826yyIb3BImQ89Fwfo1SnBbr5uDG7RBMjff92YjrMhIqIMwyDNSDFmu8ICMBWOL0OjBegzR0zpcBDnvLmMcLxPffsgje1EREQ6YpBmpBizXaoBWKBDWajJKjHmHBPpu5kpvnaob9+jsZ2IiEhHDNKMFkO2qyarBDgU8lD0AhCAY7cT1x0YGV8j9Hgk0nczUxxfA4QmOb1924mIiAzGBusmsOTSAuAloLG7Fb7BXji+dKImq8S4wCxQYaE8xKm2PdNVLQFcAL5qhBjkg7THIQdoVWyATkRExmMmLR0CJup3nXgK5t10j1x2A8B1X49Gz3empyZAA7hKNJqqJcAPe/D20NXAD3sYoBERUcowSEu1kIn6ebt34re/vQ+Vq9zwDfFiaU4LFr3kSd35GNBCioiIiJLHIC3VVCbq53q9WNzUJH/Tv1cuyZFK7CRARERkOgzSUk1jQn5hZ6f//7GU5CAiIiJ7Y5CWahoT8jvy8/3/j1qSg4iIiGyPQVqKKL05r7xyHvY7g4Ow/U4n7qyulr8xsiYaERERWQZLcKSA0psT/Xux4rxyIAtY3NSEws5OdOTn487qaqw4pxyO3SksvUFERESmxiDNYIte8mDpoK1AQJH6FeXlWFFeHrxjL9DznempPTkiIiIyLQ53GsifQYvQRajS7cb2igr4yuewuTkRERH5MZNmoMbeVqB/r+bPK91uLGtoQK63bzWn0twcYBkMIiKiDMcgTUeLXvKgsVdu7QQBIC/Czr3A4mVNRwM0hdLcnEEaERFRRuNwp06UoU3fEC8qV7mxvaoCvnPmYHtFBSrd7uCdfXL7p+JdneoHY3Nz0uJaBPwmG3BJfV8XpfuMiIjIIAzSdKIMbSpDmMUeD7KEQLHHg2UNDUcDtUNZuG7PaHkFp1YTczY3JzWuRcDhpcAJPkCC/PXwUgZqREQ2xSAtQUrdM2lVM6S3muHLk4ctFzeFD2Hmer1YvKwJjt1OXHdg5NESG2xuTvH4uhEIrXPs7NtORES2wzlpCQisexYqsL1T0PZdneElNpR5Z3V18hBnYaEcoHE+Gqk53qe+fZDGdiIisjQGaQmItGqzIz8fxR5P2Pav8r6hvo6gqopBGcXma4c8xBlqT4QaL0REZFkc7kxApAbod1ZXh7V98jqPRd7D9xl9WmR3x9cAoX963r7tRERkOwzSEhCpAfqK8nL84Ae3A0VFgCQBRUVwPrGM2TJKXtUS4JjrgK8cwP9v795j5CrLOI5/f7vYJSC3CtSqtKUGm4DRWrDSKGikUSAKWC8tbhSvDa0GKxpFmyDBO96A6BZXbABduSigjVeUq0hEBAsUSqViy8XaWryAYIt0H/847+p0ndnZHWbmHM75fZLJnH3O2Zn33XfOzLPvOXOeYbL7SUuyuJmZlY4Pd7Zgcc9MVmyrf04a23rY85VL4Jwzut4uq4D+AcBJmZlZFXgmrQUDC6aw5PFZ9G7ty2Y0dgDB/39708zMzKxFnklr0cCCKQwwgWRsaMjf4jQzM7Nx80zaeAwNwYwZvPLVLRZBHxrKanJu3AgR/6vR6WLqZmZm1oCTtGZqEiy1mmAtX57V5Kw1UqPTzMzMrA4nac20I8FqVIvTNTrNzMysASdpzbQjwXKNTjMzM5sgJ2nNtCPBco1OMzMzmyAnac20I8Hq74fBwZ0ucMvgoL/daWZmZg35EhzN1BRBj/vvR61ePsM1Os3MzGwCPJM2Hv39sGED119zDWzY4GSr6oaWwtd2gSGl+6V5t8jMzErISZrZRAwthSdWwD47QGT3T6xwomZmZm3nJM1sIv4+CH2jYn0pbmZm1kZO0swmYu8d9eN7NYibmZm1yEma2UT8vbd+/B8N4mZmZi1ykmY2EXsvhu2jYttT3MzMrI2cpLUqFV2np6e1ouv29NQ/AJOWwN96YZjsftKSLG5mZtZGvk5aK0aKro/U9Bwpug6+PEcV9A8ATsrMzKyzPJPWinYUXTczMzMbg5O0VrSj6LqZmZnZGJyktaIdRdfNzMzMxuAkrRXtKLpuZmZmNgYnaa3o74fBQZg+HaTsfnDQXxowMzOztnGSNgG3rryRjfs9m+GeHjYu+xArT/gIDA+76LqZmZm1nZO0cVq5bICTL/sU07dupieC6Vs3s/C8U1m5zJdiMDMzs/ZzkjYOS6/YzKsvPpPdt+98qfndt2/nqKEzc2qVmZmZlZmTtCaWXrGZFbutY9pfttRdf8DD9eNmZmZmT4WTtCYGh++DXYe5f//9665/4Fn142ZmZmZPhZO0JnZMzg5xfvw97+Gxvr6d1j3W18fV/afn0SwzMzMrOdfubKL3r33s2Hc7F8+fD8Bnzj+faVu2cP9++3PNiafzrrOX5txCMzMzKyMnaU0s7pnJim3rYNdhLp4/P0vWtvWw5PFZDCyYknfzzMzMrKR8uLOJgQVTWPL4LHq39sEw9G7tc4JmZmZmHeeZtHEYWDCFAaZw3XXX8ao3zcu7OWZmZlYBnkkzMzMzKyAnaWZmZmYF5CTNzMzMrICcpJmZmZkVkJM0MzMzswJykmZmZmZWQLkkaZLeLOkuScOSDhu17mOS1ktaJ+m1ebTPzMzMLG95XSdtDbAA+HptUNLBwCLgEOA5wC8kvSAidnS/iWZmZmb5yWUmLSLWRsS6OquOBy6JiO0R8UdgPTC3u60zMzMzy1/Rzkl7LvBAzc8PppiZmZlZpXTscKekXwDPrrNqeUT8oNGv1YlFg8dfDCxOP/5TUr2ZuXbbF9jahecpoir3Hdz/Kve/yn0H99/9r27/O9n36ePZqGNJWkTMb+HXHgQOqPn5ecCfGjz+IDDYwnO0TNJvI+Kw5luWT5X7Du5/lftf5b6D++/+V7f/Reh70Q53rgIWSeqTdCBwEPCbnNtkZmZm1nV5XYLjDZIeBOYBP5L0M4CIuAu4DLgb+CnwPn+z08zMzKool0twRMSVwJUN1n0a+HR3WzRuXT28WjBV7ju4/1Xuf5X7Du6/+19dufddEXXPyzczMzOzHBXtnDQzMzMzw0nauEg6OpWpWi/ptLzb02mSDpB0raS1qXzXB1L8DEkPSVqdbsfm3dZOkbRB0p2pn79NscmSfi7p3nS/T97tbDdJs2rGd7WkRyQtK/PYS1opaYukNTWxumOtzLnpveAOSXPya3l7NOj/FyTdk/p4paS9U3yGpH/VvA7Oy6/l7dGg/w1f72UqXdig75fW9HuDpNUpXsaxb/RZV5z9PyJ8G+MG9AJ/AGYCk4DbgYPzbleH+zwVmJOW9wB+DxwMnAF8OO/2delvsAHYd1TsLOC0tHwa8Pm829nhv0Ev8Gey6/mUduyBI4E5wJpmYw0cC/yE7JqOhwM3593+DvX/NcAuafnzNf2fUbtdGW4N+l/39Z7eB28H+oAD02dDb959aGffR63/EnB6ice+0WddYfZ/z6Q1NxdYHxH3RcQTwCVk5atKKyI2RcRtaflRYC2u/ADZuF+Yli8ETsixLd1wFPCHiNiYd0M6KSJuAP46KtxorI8HLorMr4G9JU3tTks7o17/I+KqiHgy/fhrsmtWllKD8W+kVKULx+q7JAFvAS7uaqO6aIzPusLs/07Smqt0qSpJM4CXADen0PvTNO/KMh7uqxHAVZJuVVbdAmBKRGyCbOcG9s+tdd2xiJ3foKsy9tB4rKv4fvAustmDEQdK+p2k6yUdkVejuqDe671K438EsDki7q2JlXbsR33WFWb/d5LW3LhLVZWNpGcClwPLIuIRYAXwfGA2sIlsKrysXh4Rc4BjgPdJOjLvBnWTpEnAccB3U6hKYz+WSr0fSFoOPAkMpdAmYFpEvAQ4FfiOpD3zal8HNXq9V2n8T2Tnf9JKO/Z1Pusablon1tHxd5LW3LhLVZWJpGeQvWiHIuIKgIjYHBE7ImIY+AZP42n+ZiLiT+l+C9k1/eYCm0emttP9lvxa2HHHALdFxGao1tgnjca6Mu8Hkk4CXgf0RzohJx3mezgt30p2TtYL8mtlZ4zxeq/E+EvaBVgAXDoSK+vY1/uso0D7v5O05m4BDpJ0YJpdWERWvqq00rkI3wTWRsSXa+K1x97fAKwZ/btlIGl3SXuMLJOdRL2GbNxPSpudBPwgnxZ2xU7/RVdl7Gs0GutVwNvTt7wOB/4xclikTCQdDXwUOC4iHq+J7yepNy3PJCvdd18+reycMV7vVSldOB+4JyIeHAmUcewbfdZRpP0/729XPB1uZN/o+D3Zfw7L825PF/r7CrIp3DuA1el2LPAt4M4UXwVMzbutHer/TLJvcN0O3DUy5sCzgKuBe9P95Lzb2qH+7wY8DOxVEyvt2JMlo5uAf5P9p/zuRmNNdrjja+m94E7gsLzb36H+ryc792Zk/z8vbfvGtE/cDtwGvD7v9neo/w1f78DyNP7rgGPybn+7+57iFwAnj9q2jGPf6LOuMPu/Kw6YmZmZFZAPd5qZmZkVkJM0MzMzswJykmZmZmZWQE7SzMzMzArISZqZmZlZATlJM7OnHUkzJOV+rTZJsyUdW/PzcZJOy7NNZlYeTtLMzPjvVdYnajbZdZUAiIhVEfG59rXKzKrMSZqZFZ6kUyWtSbdlKbyLpAtTEezvSdotbfs5SXen+BdTbD9Jl0u6Jd1enuJnSBqUdBVwkaSbJR1S87zXSTpU0lxJN6Xi0jdJmpUqkJwJLJS0WtJCSe+Q9NX0u9MlXZ3acbWkaSl+gaRz0+PcJ+lNKT5V0g3psdaUrYC1mU2ckzQzKzRJhwLvBF4GHA68F9gHmAUMRsSLgEeApZImk5XxOSTFP5Ue5hzgKxHxUrIrp59f8xSHAsdHxFuBS4C3pOedCjwnsjqF9wBHRlZc+nTgMxHxRFq+NCJmR8Sl7OyrwEWpHUPAuTXrppJd7fx1wMjM21uBn0XEbODFZFc/N7MKa2V638ysm14BXBkRjwFIugI4AnggIn6Vtvk2cApwNrANOF/Sj4AfpvXzgYOzUn0A7DlSnxVYFRH/SsuXAT8HPkGWrH03xfcCLpR0EFkZmWeMo93zyIpUQ1Zm6Kyadd+PrHj33ZKmpNgtwMpU8Pn7EeEkzaziPJNmZkWnBvHRNe0iIp4E5gKXAycAP03reoB5acZrdkQ8NyIeTeseq3mAh4CHJb0IWEg2swbwSeDaiHgh8Hpg1xb6Udve7TXLSs99A3Ak8BDwLUlvb+E5zKxEnKSZWdHdAJwgaTdJu5MdzvwlME3SvLTNicCNkp5JVhj+x8AyshP7Aa4C3j/ygJJm09glwEfS49yZYnuRJU8A76jZ9lFgD+q7CViUlvuBG8fqpKTpwJaI+AbwTWDOWNubWfk5STOzQouI24ALgN8AN5OdT/Y3YC1wkqQ7gMnACrKE6Ycpdj3wwfQwpwCHpZP47wZOHuMpv0eWXF1WEzsL+KykXwG9NfFryQ6jrpa0cNTjnAK8M7XlbcAHmnT1VcBqSb8jO2/unCbbm1nJKWL0EQMzMzMzy5tn0szMzMwKyEmamZmZWQE5STMzMzMrICdpZmZmZgXkJM3MzMysgJykmZmZmRWQkzQzMzOzAnKSZmZmZlZA/wGylaox2uc45AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "c_ols_uncensored = np.linalg.lstsq(X_ordered[:M], y_censored[:M], rcond=None)[0]\n", "fit_ols_uncensored = X_ordered.dot(c_ols_uncensored)\n", "plot_fit(fit_ols_uncensored, 'OLS fit with uncensored data only')\n", "\n", "bad_predictions = (fit_ols_uncensored<=D) & (np.arange(K)>=M)\n", "plt.plot(np.arange(K)[bad_predictions], fit_ols_uncensored[bad_predictions], color='orange', marker='o', lw=0);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the fit for the uncensored portion is now vastly improved. Even the fit for the censored data is now relatively unbiased i.e. the fitted values (red points) are now centered around the uncensored obsevations (cyan points).\n", "\n", "The one glaring issue with this arrangement is that we are now predicting many observations to be _below_ $D$ (orange) even though we are well aware that this is not the case. Let's try to fix this." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##Using constraints to take into account of censored data\n", "Instead of throwing away all censored observations, lets leverage these observations to enforce the additional information that we know, namely that $y$ is bounded from below. We can do this by setting additional constraints:\n", "\n", "$$\n", "\\begin{array}{ll}\n", " \\underset{c}{\\mbox{minimize}} & \\sum_{i=1}^M (y^{(i)} - c^T x^{(i)})^2 \\\\\n", " \\mbox{subject to} & c^T x^{(i)} \\geq D\\\\\n", " & \\mbox{for } i=\\mbox{M+1},\\ldots,K\n", "\\end{array}\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAF3CAYAAAD3rnzeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3Xt4G9WdP/73kZ0q6zgEcJC4BNsxCyGE3JoEyJIQO7gphS5dKJS4DgulxgtpuTawbdwNLK3z24LLrTSkxlAoUeMChZa29NetiBXCrUtSLg1xw4JjuwEiE7MlcUIMts73D12sy4w0kmY0o5n363n02B6NZs6RxtJH5/I5QkoJIiIiIrI2l9kFICIiIqLMGLQRERERFQEGbURERERFgEEbERERURFg0EZERERUBBi0ERERERUBBm1ERERERYBBGxEREVERYNBGREREVAQYtBEREREVgVKzC2CEyZMny+rqakPPceDAAUyYMMHQc1iZk+vv5LoDrD/r79z6O7nuAOtvZP23bdu2V0p5VKb9bBm0VVdXY+vWrYaeIxAIoLa21tBzWJmT6+/kugOsP+vv3Po7ue4A629k/YUQfVr2Y/coERERURFg0EZERERUBBi0ERERERUBW45pU/Lpp59i9+7dOHTokC7HmzRpErq7u3U5VjGyUv3Hjx+PKVOmYNy4cWYXhYiIyDCOCdp2796NiRMnorq6GkKIvI+3f/9+TJw4UYeSFSer1F9KicHBQezevRtTp041uzhERESGcUz36KFDh1BRUaFLwEbWIYRARUWFbi2oREREVuWYoA0AAzab4utKRERO4KigjbQrLy/PuM+tt96Ktra2tPv86le/wo4dO/QqFhERkWMxaFPh8wHV1YDLFf7p85ldovyNjo4W/JwM2oiIiPTBoE2Bzwc0NwN9fYCU4Z/NzfkHbj/72c8wa9YszJ49G5deeikA4IMPPsCXv/xlLFiwAAsWLMALL7wAINyKdcUVV6C2thY1NTW49957AYSX0TjvvPMwe/ZsnHrqqfjFL34BAHj22Wcxd+5czJw5E1dccQWGh4cBhFeHuO2227Bo0SI8/vjjeOedd3DOOedg3rx5WLx4Mf76178CAHbt2oWFCxdiwYIF+I//+A/VOrS2tmLatGk4//zzsXPnztj2Bx54AAsWLMDs2bPx5S9/GQcPHsSLL76Ip59+GjfddBPmzJmDd955R3E/IiIiK/EFg6h+6SW4AgFUv/QSfMGg2UUKk1La7jZv3jyZbMeOHSnb1FRVSRkO1xJvVVVj++zbt0/z8aSUcvv27fKkk06SH3zwgZRSysHBQSmllA0NDXLLli1SSin7+vrkySefLKWU8pZbbpELFy6Uhw4dkh988IE88sgj5SeffCKfeOIJ2dTUFDvu3//+d/nxxx/LKVOmyJ07d0oppbz00kvlXXfdFalLlfzBD34Q23/p0qXyrbfeklJK+fLLL8u6ujoppZT//M//LB955BEppZT33XefnDBhQkodtm7dKk899VR54MABuXv3bnnCCSfIO+64Q0op5d69e2P7tbS0yHvvvVdKKeVll10mH3/88dh9avvlK5vXN19dXV0FO5cVsf5dZhfBVE6uv5PrLqVz6r9hzx5ZtnmzRFdX7Fa2ebNsMbD+ALZKDfENW9oU9Pdnt12LTZs24aKLLsLkyZMBAEceeSQAwO/345vf/CbmzJmD888/H/v27cP+/fsBAOeddx7cbjcmT54Mj8eDYDCImTNnwu/349///d+xZcsWTJo0CTt37sTUqVNx0kknAQAuu+wyPPfcc7FzX3LJJQCAoaEhvPjii7j44osxZ84c/Nu//Rvef/99AMALL7yAhoYGAIi1AibbsmULLrjgApSVleGwww7D+eefH7tv+/btWLx4MWbOnAmfz4c333xT8Rha9yMiIjJDS08PDoZCCdsOhkLoMKk88RyTpy0blZXhLlGl7bmSUirOcgyFQnjppZfwD//wDyn3ud3u2O8lJSUYGRnBSSedhG3btuGZZ57Bd77zHSxbtiwheFIyYcKE2LkOP/xwvPbaa4r7aZmFqbbP5Zdfjl/96leYPXs2Hn74YQQCgbz2IyIiMkN/ZHhRsoECl0MJW9oUtLYCZWWJ28rKwttzdfbZZ+Oxxx7D4OAgAODDDz8EACxbtgz33XdfbD+1gCrqvffeQ1lZGVasWIFVq1bhz3/+M04++WT09vbi7bffBgA8+uijWLJkScpjDzvsMEydOhWPP/44gHAg+frrrwMAzjzzTHR2dgIAfCqD98466yw89dRT+Pjjj7F//3785je/id23f/9+HHPMMfj0008THj9x4sRYy2G6/YiIiKygMq7BJJ6nwOVQwqBNQWMj0N4OVFUBQoR/treHt+dqxowZaGlpwZIlSzB79mzceOONAIB7770XW7duxaxZs3DKKadg/fr1aY/zl7/8BaeddhrmzJmD1tZWfPe738X48ePx05/+FBdffDFmzpwJl8uFq666SvHxPp8PDz74IGbPno0ZM2bg17/+NQDgnnvuwY9//GMsWLAAH330keJjP/vZz+KSSy7BnDlzsGLFCixevDh23/e+9z2cfvrp+NznPoeTTz45tn358uW44447MHfuXLzzzjuq+xEREVlBa00NylyJ4VGZy4Umk8oTT4THv9nL/Pnz5datWxO2dXd3Y/r06bqdwyrLOJnFavXX+/VNJxAIoLa2tiDnsiLWn/V3av2dXHfAWfX3BYNo6elB//AwKt1utNbU4LjubsPqL4TYJqWcn2k/jmkjIiIiitPo9aLR603YFujuNqk0Y9g9SkRERFQEGLQRERERFYGCBG1CiIeEEANCiO1x244UQvxRCPG/kZ9HqDz2ssg+/yuEuKwQ5SUiIiKymkK1tD0M4Jykbd8G8KyU8kQAz0b+TiCEOBLALQBOB3AagFvUgjsiIiIiOytI0CalfA7Ah0mbvwTgkcjvjwD4F4WHfh7AH6WUH0op/w/AH5Ea/BERERHZnplj2rxSyvcBIPJTKW/dcQD+Fvf37sg2sphAIIAvfvGLGferra1FcjqWZHfffTcXkiciIkpi9ZQfSmsmKSaWE0I0A2gGAK/Xm7I80qRJkxIy82fy2OAg/vO997D7k08w5TOfwS3HHouvVFTE7h8dHc3qeMUotkCtKzW2T67/wYMHMTIykvE5GR0dxYEDB9Lud9ddd+Ff/uVfUBH3fGdy6NChgi2JNTQ05Ojlt1h/1t+p9Xdy3QHW3xL117KqvB43ANUAtsf9vRPAMZHfjwGwU+ExDQB+Evf3TwA0ZDrXvHnzZLIdO3akbFOzYc8eWbZ5s0RXV+xWtnmz3LBnT2yfffv2aT6elFLu2rVLzpgxI/b3HXfcIW+55RYppZRLliyRN998s1ywYIE88cQT5XPPPSellHJkZER+61vfkqeeeqqcOXOmvPfee6WUUm7dulWeddZZ8rOf/axctmyZfO+999IeZ/v27XLBggVy9uzZcubMmfKtt96SUkr5wx/+UM6YMUPOmDFD3nXXXbFynnzyyfLqq6+Wc+bMkb29vfIPf/iDPOOMM+TcuXPlRRddJPfv3y/37dsnf//738tp06bJM888U15zzTXyvPPOS6n3wYMH5SWXXCJnzpwpv/KVr8jTTjtNvvLKK1JKKa+66io5b948ecopp8g1a9ZIKaW855575Lhx4+Spp54qa2trVfdLls3rm6+urq6CncuKWP8us4tgKifX38l1l5L1N7L+ALZKDbGUmd2jTwOIzga9DMCvFfb5A4BlQogjIhMQlkW2GaqlpwcHQ6GEbQdDIbT09Bh2zpGREfzP//wP7r77bvznf/4nAKC9vR27du3Cq6++ijfeeAONjY349NNPcc011+CJJ57Atm3bcMUVV6ClpSXtcdavX4/rrrsOr732GrZu3YopU6Zg27Zt+OlPf4o//elPePnll/HAAw/g1VdfBQDs3LkT//qv/4pXX30VEyZMwPe//334/X78+c9/xvz583HnnXfi0KFDuPLKK/Gb3/wGW7ZswZ49exTrdf/996OsrAxvvPEGWlpasG3btth9ra2t2Lp1K9544w1s3rwZb7zxBq699loce+yx6OrqQldXl+p+RERETlOQ7lEhxEYAtQAmCyF2Izwj9L8APCaE+DqAfgAXR/adD+AqKWWTlPJDIcT3ALwSOdRtUsrkCQ266x8ezmq7Hi688EIAwLx589Db2wsA8Pv9uOqqq1BaGn6ZjjzySGzfvh3bt2/H5z73OQDh7sZjjjkm7XEWLlyI1tZW7N69GxdeeCFOPPFEPP/887jgggswYcKE2OO2bNmC888/H1VVVTjjjDMAAC+//DJ27NiBM888EwDwySefYOHChXjrrbcwdepUnHjiiQCAFStWoL29PaVezz33HK699loAwKxZszBr1qzYfY899hja29sxMjKC999/Hzt27Ei4P9v9iIiI7KwgQZuUskHlrrMV9t0KjK3LKqV8CMBDBhVNUaXbjT6FAK3S7c75mKWlpQjFtd4dOnQo4X535NglJSUYGRkBEO66FiJxWJ+UEjNmzMBLL72keB6l43z1q1/F6aefjt/97nf4/Oc/j46Ojmh3s6JoIBc93+c+9zls3LgxYZ8XXnghpWxqlPbbtWsX2tra8Morr+CII47A5ZdfnvKcZLMfERGR3XFFBAWtNTUoSxp8X+ZyobWmJudjer1eDAwMYHBwEMPDw/jtb3+b8THLli3D+vXrY8HXhx9+iGnTpuGDDz6IBW2ffvop3nzzzbTH6enpQU1NDa699lqcf/75eOONN3DWWWfhV7/6FQ4ePIgDBw7gqaeewuLFi1Mee8YZZ+CFF17A22+/DSA84eCtt97CSSedhF27duGdd94BgJSgLuqss86Cz+cDAGzfvj3Wtblv3z5MmDABkyZNQjAYxO9///vYYyZOnBibqJBuPyIiIiex+uxRU0QXiW3p6UH/8DAq3W601tSkLB6bjXHjxmHNmjU4/fTTMXXqVJx88skZH9PU1IS33noLs2bNwrhx43DllVfim9/8Jp544glce+21+OijjzAyMoLrr78eM2bMUD3OL37xC2zYsAHjxo3D0UcfjTVr1uDII4/E5ZdfjtNOOy12rrlz58a6VKOOOuooPPzww2hoaMBwpPXx+9//Purq6tDe3o7zzjsPkydPxqJFi7B9+/bkU+Pqq6/G1772NcyaNQtz5syJnW/27NmYO3cuZsyYgZqamlj3KwA0NzfjC1/4Ao455hh0dXWp7kdEROQkIl03WbGaP3++TM4F1t3djenTp+t2jv3792PixIm6Ha/YWK3+er++6QQCAdTW1hbkXFbE+rP+Tq2/k+sOOKP+vmBQtcHGyPoLIbZJKedn2o8tbUREROR4vmAQzTt3xrJH9A0Po3nnTgDIq6dNTxzTRkRERI5nRrqvbDFoIyIiIsczI91Xthi0ERERkeOppfXKJ92X3hi0ERERkeMZke5LbwzaiIiIyPEavV60T5uGKrcbAkCV2432adMsMwkBYNBWUHv27MHy5ctxwgkn4JRTTsG5554bWw5qZ2SGStT111+P22+/HU8++STOPnts4Yjnn38ec+bMiSXcjdfQ0IBZs2bhrrvuwpo1a+D3+wEAd999Nw4ePGhs5YiIiIpco9eL3oULEaqtRe/ChZYK2AAGbep8PqC6GnC5wj8jWf1zJaXEBRdcgNraWrzzzjvYsWMH1q5di2AwiOXLl6OzszO2bygUwhNPPIFLLrkEF154IcaPH4+f//znGBkZwcqVK7Fu3brYeqRRe/bswYsvvog33ngDN9xwA2677TbU19cDYNBGRESUFZ1jAL0wT5sSnw9obgaigU5fX/hvAGhszOmQXV1dGDduHK666qrYtjlz5gAAJk2ahEsuuQS33HILgPAi69XV1aiqqgIA/OhHP0J9fT3efPNNLFiwAP/0T/+Ucvxly5ZhYGAAc+bMwY9+9CM8+OCD+OIXv4j33nsP7733Hurq6jB58mR0dXXlVH4iIiJHUIkBPDfcAJicXJgtbUpaWsZerKiDB8Pbc7R9+3bMmzdP8b5Zs2bB5XLh9ddfBwB0dnaioaEhdn9NTQ0uueQS3HffffjBD36geIynn34aJ5xwAl577bWENUSvvfZaHHvssejq6mLARkREFKXWmqYSA9R0dBS6hCkYtCnp789uuw4aGhrQ2dmJkZER/PrXv8bFF18cuy8UCsHv96O8vBx9fX2GlYGIiMgRoq1pfX2AlGM9aj6f6me9e2CgwIVMxaBNSWVldts1mDFjBrZt26Z6f0NDAx577DH4/X7MmjULHo8ndt+Pf/xjnHrqqXjwwQfxjW98A3ZcL5aIiKhg0vWoqXzWD8d9LpuFQZuS1lagrCxxW1lZeHuOli5diuHhYTzwwAOxba+88go2b94MADjhhBNQUVGBb3/72wldo3v27MGdd96J22+/Heeccw6OO+44dGTZRDtx4kTs378/57ITERFZWrYTB9R6zvr6gKGh1O1lZehpasq3lHlj0KaksRFobweqqgAhwj/b23OehAAAQgg89dRT+OMf/4gTTjgBM2bMwK233opjjz02tk9DQwP++te/4oILLohtu/HGG3HzzTfjqKOOAhCeCdra2ooPP/xQ87mbm5vxhS98AXV1dTmXn4iIyJLSdXWqUes5EwIYHEzcVlEBtLdjIJKRwUycPaqmsTGvIE3Jsccei8cee0z1/htuuAE33HBDwraf//znCX8ff/zx6O3tTXlsdXU1tm/fHvv74Ycfjv1+zTXX4Jprrsmt0ERERFaWrqtT7XO8tTVxhiiAEACX0vCj8vLwcQIB3YqcK7a0ERERUfHKZfJgpEdtaMoUhIRAb7okugZOQswWW9qIiIioeFVWhrtElban09iIU2tq0Dc8DADYtXw5qoPB7I9TQGxpIyIiouKVx+TB/kjABgCrm5pwwO3O6TiF4qigjaky7ImvKxGRg+UxebAyLkjbWF+PK1etQq/Xi5BOkxD15pigbfz48RgcHOQHvM1IKTE4OIjx48ebXRQiIjJLYyPQ2wuEQuGfGgOt1poalLnGQqGN9fWY8dhj2Pj++1kdp1AcM6ZtypQp2L17Nz744ANdjnfo0CFHBwpWqv/48eMxZcoUs4tBRERFptHrRdUvf4nq227DsQMDeM/jQe+aNVh01llmF02RY4K2cePGYerUqbodLxAIYO7cubodr9g4vf5ERGQDPh8W3XRTLPXHlGAQU266CZg0yXKtbICDukeJiIiIEqTL8WZBDNqIiIjIcXzBIEK55HgzEYM2IiIichRfMIjmnTvRr7YIvIVys8UzLWgTQkwTQrwWd9snhLg+aZ9aIcRHcfusMau8REREZA8tPT04GAoVRW62eKZNRJBS7gQwBwCEECUA3gXwlMKuW6SUXyxk2YiIiMh+fMEgWnp6YqsgbIwsAr+2owOVAwPo93hQ/cMfWnISAmCd2aNnA3hHSqmwDgURERFRfqJdogdDoYTtG+vrY8FblduN3oULzSieJlYJ2pYD2Khy30IhxOsA3gOwSkr5ZuGKRURERMUsuXVNTZnLhdaamgKVKjfC7BUChBCfQTggmyGlDCbddxiAkJRySAhxLoB7pJQnqhynGUAzAHi93nmdnZ2GlntoaAjl5eWGnsPKnFx/J9cdYP1Zf+fW38l1B6xXf4/fj5qODrgHBjDs8aCnqQkDkRazKD+ANgDpwzXAC6AJQH2afYysf11d3TYp5fxM+1khaPsSgG9IKZdp2LcXwHwp5d50+82fP19u3bpVpxIqCwQCqK2tNfQcVubk+ju57gDrz/o7t/5Orjtgsfr7fEBzc2KOtbKylLVCq196KWMLm9YuUSPrL4TQFLRZIeVHA1S6RoUQRwshROT30xAu72ABy0ZERERWozEpbr8NukTjmTqmTQhRBuBzAP4tbttVACClXA/gIgBXCyFGAHwMYLk0u2mQiIiIzKUxKW6l263a0lbldqO1pgaNXq/epTOMqUGblPIggIqkbevjfr8PwH2FLhcRERFZWGUl0KeQcCIpKW5rTU3KjNEylwvt06YVVbAWZYXuUSIiIiLtWlvDY9jiKSTFbfR60T5tGqrcbgiEW9eKNWADrJPyg4iIiEib6GSDlpZwl2hlZThgi5uEEE310T88jEq3G49On160wVoUgzYiIiIqPo2NqisXJCfS7RseRvPOneGHFXHgxu5RIiIispXo2qLxDoZCaOnpMalE+mDQRkRERLbgCwbT5mbLlALE6hi0ERERkfF8PqC6GnC5wj99Pn0PH+kSTZdMt9Lt1vWchcYxbURERGSs5BUM+vrCfwOq49KypdQlGq/YEukqYUsbERERGUvjCgb5SNf1WeypPqLY0kZERETG0riCQT7UVj/QurZoMWBLGxERERkraaWCjNtz0FpTgzJXYlhjhy7ReAzaiIiIyFgaVzDIVTSR7sFQCCWRbXbpEo3HoI2IiIiM1dgItLcDVVWAEOGf7e26TEJInjU6irEWNjsFbACDNiIiIiqExkagtxcIhcI/DZw1aodEukoYtBEREVHRUps1WuyJdJUwaCMiIipWBiesLQZqCXOLPZGuEqb8ICIiKkYFSFhrRdFJB33DwyhBeAybACDj9rHbrNEotrQREREVowIkrLUapUkHQDhgE5Hf7ThrNIotbURERMWoAAlrrSbdUlUS9kqkq4QtbURERMWoAAlrrSbT5ALF+2007o9BGxERUTGJBiF9feGcZ/F0TFhrJb5gENUvvZQwbk1JyuSD6Li/vj5AyrFxf0UauDFoIyIiKhbxQQgQDkSigZuOCWvN9Py6ddh99NEIuVz429FH42vf/S5WdHcrrisaT3Hygc3G/XFMGxERUbFQCkKkDAdsvb2mFElPz69bh7k33ogJkQDt+GAQ97W1YRjAxvr6lP2js0er3G7lFRBsNu6PQRsREVGxsFkQEs8XDGLRbbfFAraoCcPDWNvRkRK0CQAjtbXpD1pZOdYqmby9CLF7lIiIqFjYdPJBNJXH8QMDivdXKmzXlDzX4IXqC41BGxFRrmw0K42KhM2CkKhoKo9+j0fx/uTtmpPnGrhQvRkYtBER5cJms9KoSNgsCImKpupY3dSEA0ktaAfcbqxuaor9XVFaml3yXIMWqjcDx7QREeUi3ay0Iv5QoCLQ2Gibayy6JFU0lUd03Nrajg5UDgyg3+PB2uZmdC5dqj7ZwEHY0kZElAsbDwhPx+P3s0uYdJG8JFXUxvp6TO3sxMRAALvXrEH7z36G0NKl6G1oQKPfb1JprYEtbUREubDZrDRNfD5Ma2sDoh+yDlmgnPShtNC7miq3GxtefRWLbrpprEWb15v5LW1CiF4hxF+EEK8JIbYq3C+EEPcKId4WQrwhhPisGeUkIkpg0wHhabW0oCQ5wWkRJyotKkU66SW6koEIBHBpXILcdAGbANC7cCEW3X67rRLj6sEqLW11Usq9Kvd9AcCJkdvpAO6P/CQiMk/0m35LS7hLtLIyHLDZuQXAoV3CpotOeimyFqdo92d0gfdMS1BFxVJ58HpLYXpLmwZfAvAzGfYygMOFEMeYXSgiIjvNStPEpjnCLK8Il2LyBYO4rLs7FrBplZDKg9dbCisEbRLAfwshtgkhmhXuPw7A3+L+3h3ZRkREhdTaitHkhKZ27xK2giJrcYq2sKXrAlVS5XYnpvJw4hCEDISUWhssDSqAEMdKKd8TQngA/BHANVLK5+Lu/x2A/09K+Xzk72cB3Cyl3JZ0nGYAzQDg9XrndXZ2GlruoaEhlJeXG3oOK3Ny/Z1cd4D1d3r9D/vtb3HKhg1wDwxg2ONBT1MTBhTWhLQjs177M5Yvx/hgMGX7Ia8XLxv8WRcvU/39ADoApJY0PTeAVQCUriKP34+ajg5LXG9Gvv51dXXbpJTzM+1netAWTwhxK4AhKWVb3LafAAhIKTdG/t4JoFZK+b7acebPny+3bk2Z06CrQCCA2kxrntmYk+vv5LoDrD/r79z6m1b35DFtQLjFqcBJdePrH50J2j88jCNLSnAoFMIBDfGEQLh7LeNC7xZk5OsvhNAUtJnaPSqEmCCEmBj9HcAyANuTdnsawL9GZpGeAeCjdAEbEZGjFOmsQkV2qoueLLQKgi8YxOQtW7AiMhNUAhgcHU0J2Br8fuxavhyjS5di1/LlaPD7UQLg0enTIWtrMVJbC1lbi96FC4siYLMKs2ePegE8JYSIluXnUsr/XwhxFQBIKdcDeAbAuQDeBnAQwNdMKisRkbUU6axCRXaqixEKvArC8+vWofq223DswADe9XjwveZmPLB0KUR3d8ZZoA1+Px5oa8OESHqP6mAQD7S1YeWxx2KRQ1tp9WJqS5uUskdKOTtymyGlbI1sXx8J2BCZNfoNKeUJUsqZUkpj+z2JiPRidMtREc4qVGWnuhQxXzCI5jVrMPfGGzElGIRLShwfDOKu229Hg9+vKW3H2o6OWMAWNWF4OJx3jfJidksbEZE9FaLlqMhmFaZlp7oUmfiVCgSAnvZ2xaBrbUdHbG3QdCoHBpTv4GuZNyuk/CAisp9CtBzZKY+VnepiYdEVClyBACZv2YLyzZtj49OA8CQBtaBLNRhL0u/xKN9R6NfShmMkGbQRERmhEC1HVsljpceHo1XqYmPxC7SrTSAA1IMute0ThEBFaSkEwrNBd69ZY/5rGW3p7usDpBxr6S7ywI1BGxGREQrRcmSFWYV6fThaoS5mKGBrUEtPj6YVClY3NeFAUhLlA243Vjc1AQin7QAii7pPn46hJUuwd9EihCKzQRetXGn+a2nTMZIc00ZEZITWVuXcWnq3NhR4VmGKdB+O2ZbL7LoUWoFnzPYnjVNTEx23trajA5UDA9jt8eD7zc3oXLpUe141s19Lm46RZEsbEZERnNJyZNMPx4IocGtQZfISZBHJOdW+6vdjY309ap96Chvffx+Ve/ag/bbbsAnIPq+aWePKbDpGkkEbEZFRnLCgvE0/HAuiwAFva00NylyJH/vRnGrVkfQe1cEgHrnzTsh3380/8a2Z48psOkaSQRsRkd0UsnXDph+OBVHAgDea1uNgKISSyLaKkhL8QCGnWunhMdE0AAAgAElEQVTHHye29kWupyVLl2Z3PZk5rsymLd0c00ZEZCeFXlkgesyWlnALUWVlOGAr8g/HgjB43GNy/rXoPNFRAGUuF+456SQcnymnWtz1JIDsriezu87NHldnALa0ERHZiRmtG07oBjaCQa1ByeuDAkhZyeBgKISWnp7MrX35XE/sOtcdgzYiIjsxu3WDsqNzwBvNxTY4Oppx3/7h4czd2/lcT+w61x27R4mI7KSyMtyFpbSdbCe+C7QE4a5PrSrd7szd2/lcT+w61x2DNiIiOylUfjgqCKWgrKKkBBACgyMjKWPVtCpzudBaUxP+I93Yr3yvJxuOKzMTu0eJiOxEbZwUYLt1GO0quj6oCARwady4tGhQNjg6isGREQCpY9XSiV/JoH3aNG3pPOKuJ2mHWZhFvh4pgzYiIrtJHicF2HIdRjuKXx8UyC4oS6eitBSPTp8OGVlqKqv8a5HrafOmTdafaJIuKLPBeqQM2ojspMi/RZJBCjmjlNdgXrSuD6pVdH3QvYsW5ZcotxhkCspssB4pgzYiu7DBt0gySKFmlPIazFq0K9QVCGDyli2xFrZ8lblc2DB9evpWNbsF2JmCMhvMrGbQRqQXs98AbfAt0nL0ek3NvjYKlS+L12BW4rtCJaApTYeS6Fi16EoHmsasGRFgm32dZwrKbJA3jkEbkR6s0MJgg2+RlqLymnr8fl2OU9Bro1D5ssy4Bs0OFHLkCwZxWXe3pq7Q5KCsoqQEFaWlEAgHaNGxaiO1tdrHrOUbYCc/7ytXmn+dZwrKbJA3jkEbkR6s0MJgg2+RlqLymtZ0dOhynIJeG4Vah9Hoa9DIQKGAwV+0hU1ru1pyULZ38WLsXbQIoVwmFUTlE2ArfRFZv9786zxTUGaD9UgZtBHpwQqtXDb4FmkpKq+dW22txiyPU/AW0EIsNWXkNWhkoJBNa2gewV10/NoKlRa2Br8fu5Yvx+jSpdi1fDka/H5Uud3GTCDIJ8BW+iIiVea5FvI61xKUFfmSawzaiPRghVYuG3yLtBSV127Y49HlOLZsATXyGjQyUNDYGurx+3Nu2UtO5ZGswe/HA21tqA4G4ZIS1cEgHmhrw4ZXX82uLlrlE2Bn8/xmc53r0dpZ5EFZJgzaiPRQyFaudG9sNn/DKiiV17SnqUmX49i2BdSoa9CoQCHdsZO213R0aG7ZS54V+q8Zxq+t7ejAhKSAbsLwMBbdfru2OmQrnwBb7fkVIvHvbK5zK4z9LAIM2oj0UKhWLr6xFY7KazpQX6/LcQwJqPUel2WlQf5GBAqZjp20Xa1rXPb3JwRo5Zs3Y0VkJYPorNBM0w0q1brdjexezDXAVvsictVVuV/nVhj7WQQYtBHppRCtXHZ7Y7NSUKBEr9e0ENeG3gG9XrNn9WJEoJDp2EnBn1rXeJ/HkxCgHVDrtlVR5Xbj4HHHKd/pcun7f6FXF6TSF5F163K/zq0y9tPiGLQRFRM7vbGx1VBfegf0es2e1YsRgUKmYycdq6epKSW4O+B2Y3W2XeYR8Qlwy//rv1IDRwAYHdXv/0LP/zm9v4hkau20+he8AmHQRlRM7DSo3W6thmbTO6DXa/asnoxssdRw7IH6eqC9Hbu9XoSEQK/XiytXrcLGbLvMEc65lpAANxo4lpSk7qzX/4WW/zmzgqN0rZ38ghfDoI2omNhpULudWg2tIJeAPt0HtF6zZ+2msRGVnZ1YsXo1AGDD2rWx9BxalblceGT69NRUHo2N4aBRiR7/F5n+58wMjtK1dvILXgyDNqJiYqe0HnZqNbSCbAP6TB/Qes2etYHoTNClACZv2aKankMtcJsgRMIKBmmXmDLy/yLTsc0OjtRaO/kFL8a0oE0IcbwQoksI0S2EeFMIcZ3CPrVCiI+EEK9FbmvMKCuRpdglrYedWg2tINuAPtMHtF6zZ3NlUjddcqoOpZmgrSrpOdZ2dKQsMbVh+nQMLVmifQUDI/8vMh3bqsERv+DFlJp47hEA35JS/lkIMRHANiHEH6WUO5L22yKl/KIJ5SMiI0WDiZaW8IdCZWX4w6PQQajPZ34Z9NLYqL3sWj6glY4XCORUtKxEWwGjQWW0FTBapnwPHwyipacHfcPDKAEwivB6nodCoYSZn2oLuKul56geGMDexYvzK5yR/xeZjl1ZGX6uk5kdHLW2Jl4PgGO/4JnW0ialfF9K+efI7/sBdANQmfNMRLZkdquhkwc4W7n1woBuumgLmggEcGmk5QxAbP3PbFJ19KuN68v03GltPTRrwoVVW7/tNCwkT5YY0yaEqAYwF8CfFO5eKIR4XQjxeyHEjIIWjIjszewxPGYy+gM6n+7NHLrpkrs1Jz//PEQggNJAICVQyy6LWqrVTU044HYnbhQiHPSr1bUYviBYOTgy+wueRQiZZRJA3QsgRDmAzQBapZRPJt13GICQlHJICHEugHuklCeqHKcZQDMAeL3eeZ2dnYaWe2hoCOXl5Yaew8qcXH8n1x2wV/2XLF0KofAeKIXA5k2bFB9jp/p7/H7UdHTAPTCAYY8HPU1NGcesaam/x+/HtLY2lMSN+xp1u7Fz1SpNY+LOWL4c44PBlO2HvF68rPDe7gfQBkB5VU9jNPj9WNvRgcpgEAJA/LoMSnXNtk5WZKdrPxdG1r+urm6blHJ+pv1MDdqEEOMA/BbAH6SUd2rYvxfAfCnl3nT7zZ8/X27dulWfQqoIBAKora019BxW5uT6O7nugM3qX12tPIanqir8bV6BreqfA031z+F5TZA8pg0ItwIqtPr4gkFc1t0N5dFn+psgBMaXlODDkRFUut3Y/pWvoHz37tQdk+vqcikvcC+EepoPi1F97e00LjQNI//3hRCagjYzZ48KAA8C6FYL2IQQR0f2gxDiNITLO1i4UhKRrVl1DE+xy3cWosZuOl8wiOadOw0L2JJTdbQAKTNBy999V/nByXW18hjCfBRDt6+NmDmm7UwAlwJYGpfS41whxFVCiKsi+1wEYLsQ4nUA9wJYLs3uzyUi+7DyGB6riYxRW7J0aeYxanoEKBrGMLX09OBgjq1U0e7M6PoDWlJ1KHbsaq2rXb8gmDEu1MFLWpmW8kNK+TwShwEo7XMfgPsKUyIicqRs0mQ4VVx3pQAyp+AwOEVDfMqObAiEJyFUud1oralJny9NKy11jXYfHjwYXqZqdDT8BUHPbkSzuigLndvN4HQwVmeJ2aNERGRh2bamqLVgAnm1kPiCQUzesiWW7FaNAFBRGm6TiLakVbndeHT6dEgtCW6zkam1Nr77EAgHbNGgTs+AzawuykJ3+zp5xjfMTa5LRETFIJfWlOQWzDxbSKLj1zJ1h5a5XOmXiTJCutbadEGGXkFbIc6hptCJb626akOBsKWNyInsOiZEqV52rWsh6dGakmMLSTT/2oru7owBW8Z1Pc1QiCDDzECm0ONC7TqhQyO2tBE5jV3HhCjV62tfC3+QfPLJ2DY71LXQ9GhNSRNYaF1WKp0qtxu9CxdqL0+hFGJpKLOXnyrkuFCHL2nFljYip8l1TEg2swfNoFSvTz8dC9iiHDT+RTdxrSkyh9YUXzCI3SpLP/V6PHkvK1XmcqG1pkbTvhnp3TJbiFmjdp2ZqsThM74ZtBE5TS5dKXEDnYVVczFl0xXkkPEvuoqk4Ni8aVNWywhFx6LdrLD00wG3G6ubmvJaVqqitFS/LlEjBvQXIshwWiDj4CWtGLQROY2WMSHJrQ3XXWf9GVvZdAU5ZPyLmZLHom2sr8eVq1ah1+tFSAj0er24ctUqbNSwrJWSaC61vYsWhQM2PVrIjJqZWIggw8GBjJNwTBuR02QaE6I0NkyNlVqslOo1blzimDbAvt1GFqI203NjfX3OQVqU4uxQvcZpOnxmIlkfW9qInCZTV4pSa4MaK7VYKdXrpz8FHnrIOd1GFpHPSgXpqHaF6tVC5vCZiWR9DNqI7Cpdd1G6rhStrQqFbLHS2vWlVC92GxWULxjMeqUCQNuyUrGu0GR6tZA5aUA/FSV2jxLZUT7dRWrpAyoqgPJyyP5+iEIuk2PXFCU2FO0WTSc+pQeEwIcjI6jMd1kpvVJexLc2F3o5KCINGLQR2VE+GdLVxrzdcw/Q2IjNgQBqa2t1L7IqM7O9U1pK+dXUGLpSgZ65u7gWLVkYgzYiO8qnu8hqrQ0cHG6aaFDWPzyMIyMtY4MASgIBjGJsAXYgfcAGwNiVCqx2zerJrIXgyZIYtBHZUb7dRVZqbTA727uNKQZlIyOxVrP4oGxwdCwsi/6mNb9aldtt/NJSVrpm9WLFoQEMIk3FiQhEdmSnAdV2qouJonnTXIEAJm/ZgvLNm7EishKBRDgoGxwZAZB9UJaOrqsVOI1ReeNy5PH79U8+TFlh0EZUaEYsYJ58TMA+GdKdlu1dB1oCNK1LROXDkgu4FxOLDQ2o6eiwVBDpROweJSokI7o71I7Z3h5OcWEHduz6MoAvGMR1b72V0JUZ/3uhGDrpwOr07D602NAA98CA8h0cX1owbGkjKiQjujss1oVChRFtTROBAEoDAYhAAJd2dxc8SEvOr2ar1rWkFmyP3595fz27Dy02NGDY41G+g+NLCyZj0CaE+KYQ4ohCFIbI9ozo7rBYFwoZL5oPLZrEVs8xaMmiQVk02S0wFqBdEwhg36WXQi5dipHLL4d89130LlxofsCmxxAEhQBsWltb+mPp/QXKYkMDepqaLBVEOpGWlrajAbwihHhMCHGOEEJkfAQRKTNimRwuveMYyYuw62mCEClBWZXbjUenT4esrcXexYuxd9EidAEYqa2FfPdd3HvHHSjfvdtag9L1au1SCMBKhofTB2BGfIGy0IoeA/X1lgoinShj0Cal/C6AEwE8COByAP8rhFgrhDjB4LIR2Y8R3R1W6UIxYoIFJXSDXhqZTJCvy/1+9C9fjtGlS/G35cuxZccODC1Zgr2LFkHW1oaDstra9K1mVu2W16tcuQRgTvgCZaEg0ok0jWmTUkoAeyK3EQBHAHhCCHG7gWUjsh8jujus0IWi91gePcpTyABSp/Olm/UJ5Nb9GW1Bi67huWXHDvz0rrtwfDAIl5SYEgxi0U03ZV/mQnXLZ/vc6lWuXAIwq3yBItvKOHtUCHEtgMsA7AXQAeAmKeWnQggXgP8FcLOxRSSyGSNmQpo9u9JKS00VOiFpjudLTmx7KBRKSMOR64SCaELcKrX1PBsa9HmtCjGzMZfnVq9yKSyNNep2oyRdAGbnlRmMwmS9WdHS0jYZwIVSys9LKR+XUn4KAFLKEIAvGlo6IioOVpoMUehuOw3nMzpvmtIYNNWuTb1eq0K0KuXyWmopl5bWO4UW7J2rVmUOKNh9qJ3VWuiLgJYxbWuklApfWwApZbf+RSKiomOlsTxagpI8ujM9fn/iY5VadSLn8wWDmLxli2GJbctcLmyYPl3bGLQovV6rQnTL5xJgZipXNoFCUgA2UF+fd5UIY/9/K1ZYc1ykhTFPG1EuOOg+Ub6tLno+n5mCkny+3ft84bQP8Y9VmVDf6/EYkjcterac86Hp2UJmdKtSrgFmunJZdQKFU8T//6lhuiJVDNqIssT19xTk0+pS6ISk+Xxot7SE0z7EkxLJyTcOuN1Y3dSke960itLSzN2fmVhh4opWRnTBWqkr34mU/v+S2Wm2rc4YtBFlqajX3zOyhTDXVpd8Wz6yXXc1nw/tNPv0er0ICYFerxdXrlqFjVl0pTX4/di1fDlG6+rw6dlnY7SuDn3Ll+PKTZtisz43TJ+OvYsW6ZO4tljGXeUbYCpd71bqyneiTP9nnG2bFtceJcpS0a6/V+hZlVrlE0Tlsu5qjrMLfcEglng8mBIMphbV68XUzs7M5VXQ4PfjgbY2TIi04LkiSXMrg0G0//CHaJ82zbpBVSHkOjNa7dq47DLgkUcSvygwUCgctf8/IByUc/ZoWqa2tEVWWNgphHhbCPFthfvdQohfRO7/kxCiuvClJEpUtOvvWXUsTz4tHzrOLnz+5psTZnhOfv75lHU9b25qwgG3O+Gh0a5QrZLzprU/8kgsYEthhdenWKldG888Uzzdw3ak1uW9YYO1W30twrSgTQhRAuDHAL4A4BQADUKIU5J2+zqA/5NS/iOAuwD8oLClJEpVtOvvWXUsTzbjlpK7u9LM3FSV1OU2NGUKmr/1LSw+5ZSEGZ6DIyMAEtf13FhfjytXrcqqKzR+4sCG6dNjKw+EIuPSyt99V72smepC6tJd74XoHuZkJWXFNKbSgsxsaTsNwNtSyh4p5ScAOgF8KWmfLwF4JPL7EwDO5tqnZDbD198z6s3eqmN5tL6JK01YUHs70Di70Pf++/D6fHhg6VLNxd1YX4+pnZ0o2bQJUzs7FQO2+Ja0jBMHMpW10K+PXYINM6935h9Lr1jGVFqQkDrlC8r6xEJcBOAcKWVT5O9LAZwupfxm3D7bI/vsjvz9TmSfvQrHawbQDABer3deZ47jS7QaGhpCeXm5oeewMifX38i6e/x+TGtrS5ihOOp2Y+eqVXnniNLr2Ga99mcsX47xCuPJJMZas4Ds6rQcQPSIDX4/1nZ0oHJgAP0eD1Y3NWU1mSDKi/A3UK2UXpcovV77fMqSXIZi+d834n9Ja93VrtVDXi9eNvizyUjF8tobxcj619XVbZNSzs+4o5TSlBuAiwF0xP19KYAfJe3zJoApcX+/A6Ai07HnzZsnjdbV1WX4OazMyfXPue4bNkhZVSWlEOGfGzak7lNVJWX4u3niraoq5/JmXYYMTHvthVB+bqLPTw51El1dEl1dsqGlRQ653QnHHHK7ZUNLi0TcPru8XjkqhNzl9SbcF72Vbd4sN+zZk33doq8LIGVJyVidcnh98qLh+iuq/30drvd4muuudq0Kkdf5zVZUr70BjKw/gK1SQ+xkZvfobgDHx/09BcB7avsIIUoBTALwYUFKR6Qnrd0lRo87K+ZuCbVuraqqrOsUXVYq2s+wtqMjZTLAhOFhrO3oADA2w7M6ssh6dTCIB9racOWmTYmTCnJJdguMvS5SAiMj4Z9mvD5WHfeYK7Oud7OHItili5tSmBm0vQLgRCHEVCHEZxDuqXg6aZ+nEV6sHgAuArApEpESFRetsxzNfrO3siwmLCSv9Tn5+ecV1/2MqlRJ41I5MIASqAd17T/7WcKkAl1yqOktmw9wp11/RgU3hViXVQ3H09maaUGblHIEwDcB/AFAN4DHpJRvCiFuE0KcH9ntQQAVQoi3AdwIICUtCFFR0NqCYeabvV6M+iDUMGFBba3PwZGRtOt+9qukcXFVVmKkthbVxZ6bT+sHuFHXnxVbfowMbsycIWnV1D6kC1PztEkpn5FSniSlPEFK2RrZtkZK+XTk90NSyoullP8opTxNStljZnmJcqa1BaPYp8Mb/S0/rrvL96c/obqmJiWXWi5rfa5WyL+WEKwUawtUNh/gPt/Y/iUl4W16XH9WbfkxOrgxq2vWbl3clIDLWBEVQjYtGMU87sygD8Jod2dygBbt4ozPpZaLjfX1+M63v60eLBdrC6jWD/DkRbxHR8fql+/1Z9WWH7sGN8X6BYM0YdBGVAjF3oKmlY4fhPGBmp4BmpIylwunX321erAc9/rJYnr9tH6AGxlYWTU4smtwU6xfMEgTBm1EhVLMLWha6fBB+Py6dfjb0Uej4ZhjELjgAjT4/boGaMkqSku1zfqMvH6bN23S5/UrxDgvrR/gRgZWVg2O7BrcOOULokMxaCMi/WT5QZjc7fnV734Xc2+8EccnpdZo8PtzLlL8Wp8VJSUJKTo2TJ+OvYsWFX7WZ6HGeWn9ADcysDIiONIj4DU7uDEyaHfCF0SHKjW7AERkI9EPh5aWcCtNZSXQ2gpffT1aXnoJ/cPDOLKkBBACgyMjEBjr5hxF+nxpWlcniB6zyu1Ga02NNdNwpOuO1PsDtrEx8zFbW8NBY3yZ9Gp1Urkmcq5nNOCNljUa8MafK5uymRHQ6FkHchS2tBGRLmKtZscdh9KHH4bYtAmTH30U5VOmKKbgAFLHpaXLl6YkunxVZK6jtrU+rcBq47yMbnXSs+XHqhMbslGoOlgx1QrlhUEbkZHs/qYZqZ90ubB49mz80+9+B2BssoBaXjQ1avnS+j2etAHaSG2t9QO1eNl0RxbqGiqWLjWrBby5KEQdrJpqhfLCoI3MZeegxu5vmnH1E1KiUofxZ0r50g643Vjb3Fy8AVoynw8YGkrdrtQdafdrKMppqzYUog52aJGkFAzayDxGfyCZHRDa/U1ToX7x63XmYmN9Pa5ctQq9Xi9CQmC314tX77wT7bfdVpwBWrLoNT84mLi9okK5O9Lu1xBgnVUbCqkQdbBDiySlYNBG5lH7QLruuvyDLSu0UNj4TdMXDCKkUg+18Wdqkrs9XzzvPLzw+utwhUKYsmcPFq1cmXtBrUbpmgeA8nLl7kgbX0Mx2QamZs/61EMh6mCHFklKwdmjZB61D57BwbGWiFxnVRVydp6aysqxDPPJ24uQLxhES08P+gCI7m6c6fGgOhhM2U9tXFrUBCEwvqQEH46MoNLKMzyNkG0QZrNrSFEugalZsz71ZHQdjJwRTKZh0EbmUftASpZLsGWFFooiedOMBWPDwyhBeBJBRSQtx4cjIziypASHQqGECQUS4fFnD7S1JaToOOB2Y3VTk+JxHBegKck2CCuSaygvTghMzaB3qhWyBAZtZB6lDyQ12QZbVvggsMibZjQoS86RFg2sknOlAUhYdF1tAfZo3rS1HR2oHBhAv8eD1U1N6KyvR6i21qDaFLlsgzCLXEOGckJgahY7tEhSAgZtZB6lD6ShodRB2kD2wZZVPggK/KaZHKAlt5DFB2B6rOG5sb4+JeltVdLsT0vx+cwNgHIJwuz+weuEwJRIJ5yIYGdmz57UIjk31D336DOryg6DlbPkCwbRvHNnQhLbbHKk6aHM5UJrTU1Bz6mZFSanAMWTD62Q+JwQacKgza6s8gGVLT2DLQd9EPiCQVzW3Y2DoVDBzqmU7FbTwut6yfZLCbPQE1GRY/eoXVlh9mSurNodFOlaW2JyF47SxIH4cWmFUFFaintOPNG8SQW5rN1YyCz0XFOSiAzAlja7ssLsSTtJyv5vVstlfBcooM+4NCC11ayipAQVpaUQSb9Xud1oAbB30SJzZ4Hm0mpWbFnolVrs2IpH5GhsabMrK8yetBODWy4zzfCsiNuWj2iOtPhjV2WZiiMQCORVBl3k8qWkEJNT9PqypNRi97WvhYcMfPLJ2Da24hE5Clva7MoOS71YiYEtl0oTCKLBWXwKjmwDNgEktJBtmD4dQ0uWYO+iRcWxhme6VqVcWs2KKQu90peETz8dC9ii7LaklZGc3krp9PrbBIM2u3La7Emj35AM7Fpr6enRfQJBmcuFR6dPx95FixB69130NjSg8ZhjiufNOtNEmly/lBg9OUWvL0vZfBko1iEPhQwiinVill6cXn8bYdBmZ3aaPZnuDb4Qb0gGtlz2x60okI/ouLSEWZzF+madaWyYVb+U6FWubL4MFOOQh0Jfl4WaOWxVTq+/jXBMG1lfphl5hZgpG5cAVPb3Q+Qwe1Rtuah8Zn5mHJdWrLOItXRHW3WWsR7lUhp/N25c4pg2oHiHPBT6unT6xCyn199G2NJG1pfpW2I2b0j5dMlEWi43b9qUseXSFwyi+qWXIAIBlAYCEIEALu3uTpn1OTg6ilw6RstcLmyYPj3zuLRifbMuxExPK1NqsfvpT4GHHrJe62IuCn1dOv16cnr9bYRBG1lfpjd4rW9Ial0yK1fqMrYmPlBTCtC0tqZFJxAAyik4skpiW6xv1pxIozy8wS5DHgp9XTr9enJ6/W2EQRtlJ9JStWTp0sINas/0Bq/1DUmtxW79+rzH1iTnT8s3b1ryDM+9ixeHJxVkO9uzWN+srTpmjfRR6OvS6deT0+tvIxzTRuqSF9c+91zgkUeAgwfDg94LlScqU34trQtOq7XYJa/PqWFsjdL4NL1U6rngejEvxm3VMWuUPzOuS6dfT06vv00waCNlSoP/16/PKcDJm5Y3eC1vSGoJh5UoBHixQA2A6O6OtabpGbAZsuA636zJinhdEmXNlO5RIcQdQoi/CiHeEEI8JYQ4XGW/XiHEX4QQrwkhtha6nI6m1JWYHLBFFWJQu9axPOkmGih1yahJ6pL1BYO44s3U7s8Gvx+7li/H6NKl2LV8ORr8/rTbEYo8eDT8s+HpTdj1leUYravDp2efjf1LluCsY09Hc7kPXxU+9IlqhIQLg2IyBl2TERIu/K2kOuX+6DYhgNLScA/I5MnhW6ZtzeU+/K2kGiEhMCJKERIidjyXK/5+F/pENdrr3k3Z9lXh03xspbrkUu5st7lcuTwnqfWrq1ui+/lyqYvW10Dv505r/Y14Xcw6X/R/7ay6pWmv2WKoSz7bli5dYpu6ZHO+Qo0E0kRKWfAbgGUASiO//wDAD1T26wUwOdvjz5s3Txqtq6vL8HMYbsMGKauqpBQi/HPDhrH7hJAyHKZlvlVVmVSBJBs2SFlWlli2srLEekXrnK4+SY+5+pd7JPxdEl2Jt4aWFjnkdic8dsjtlj/60pcUtzd84xaJs/fENjdggxxCmWIZPsY4eQifUS2j0v1DKJMN2KD5ZctUhiGUyR/h6pT7tZ473bHzLbeRN6Vys4zOvfG55g0Ifyy0tLxp2McXgK1SaoiftOxk5A3ABQB8KvcxaDNKpgBHLbBJDuaSgyIzqZVZKahMF5Ru2CCv/uUeWfL4ixLPdoVvXam3XV6v4uM/dbkUt+9CVcKmXVApbx635HNkumUqw6coyfnc2dQv23IbeVMrN8vozBufa96iN6/3Y8M+vrQGbVaYPXoFgN+r3CcB/LcQYpsQormAZbK/TLnP1GZ3XXUVUFUFqTQDyey17bLJ/aQ2I7WqCiv/oR73l+3E6OTh8AAClf+SyobFgGUAABymSURBVIEBxe0lKktSVaI/7d96yPaYmfYvyWLEXj71M+K5yJVaWVhGZ+JzTVEDAzpOEsuRYRMRhBB+AEcr3NUipfx1ZJ8WACMA1D7dz5RSvieE8AD4oxDir1LK51TO1wygGQC8Xi8CgUC+VUhraGjI8HMYaUl/f2zZo3iyvx+bAwHguOPgueEG1HR0wD0wgGGPBz1NTRiorwe+8hUMDQ2hvLw8/KBAAB6/H9Pa2lASXZKprw+jX/86dnZ3hx9TAGd4PBgfDKZsP+Tx4OWk18qzYkVieQGMut24+6wm3H9Yt6b/jP6jPKgeSD3fKEpQqhDs9KMy5e9qaJwYoVHyObTsn64ManXRcu5s6pdtuY2kVm6W0Zn4XFPU5MkfIxD4H3MLoaU5zogbgMsAvASgTOP+twJYpWVfdo9qkE1XooKU+ud5PF1oGdOWvH/cmL4Hr/uxbLjpu3KX1ytHhZC7vF7Z0NKi2DWKP3bJhlN+rDjWRWkcWLZjvjimzbxbMYxhKoYy2uXG55o3wOFj2gCcA2AHgKPS7DMBwMS4318EcI6W4zNo0yDbACdJSv3VxogJoX/Z00k3uUJFdPxaw+rUiQUhQA4cdlhi8Pb7zfIz5+6RFRXhN/ReVMlRCNnvqpJXTtgghZDyygkbZL8rvL0XVbIBG2RFhZQVFeFDl0SGio3tFx4/NgrEjhN/7L2okHtFRcJ5lM4df2yl82VbBrW65FM/pbrkUu5stwmRy3OiVL+Q7ufLpS5aXwP9nztt9TfidTHrfFr/14qhLvmdL2Sjumh/TPSjxMjPfasHbW8D+BuA1yK39ZHtxwJ4JvJ7DYDXI7c3Ee5W1XR8Bm0a5RDgRFmypS0HV/9yj8TvN6edWCARmf25ukWWPP6i/NK9W80utqlsce3ngfXvMrsIpnFy3aVk/a0QtJmSXFdK+Y8q298DcG7k9x4AswtZLsfRM7llplULTLDyySDaQz0YPXIY2F8CIQRk+QggEZ5csK8EmDQaW+BTbWIBAEwYHsbPfRvw89bvF/VYRiIiKl5WmD1KdmDU2nbpZqTG3Td41PFYcd2tEJsCEM8GILoCuP/w7rEZoJNGIQ8bCf9egvCq7HEBGwD0ezzpy9LXB1RXwxNNlEtERFRADNpIP1pXLdAqupRW0mLuD12/DiuuuxUHvv712H0Ve3fjJz/5LzRs8o8FZXFXt+oKBXFWNzXhQKZ1P/v6MK2tzULpsYmIyCkYtJF1qeSSW7rxNnz/F+sxIS5dBxDuwlzb0ZFymAa/Hw+0taE6GIRLSlQHg3igrS0lcNtYX48rV63CBxMPiy1TpaRkeHgsnx0REVGBMGgj85PiqlFJllv5wYDq+DOl7Ws7OjQHeBvr6nHLw29BbNgQ7uLNsmxERERGYdDmdCpdkJYI3FRWLej3eFTHnylt1xzgHXLh6o+mY92F3rGuXrXATW1FBSIiIoMwaHO6TMtZaWFUS53CUloH3G6sbmpSHH8WvS9BKLxygZL+yR6IfaVACCjZ68bVB6eFA7YMZRh1u02dFUtERM5kSsoPspBs1utUEm2piwZ+0ZY6IP+JCJHHD17/bRwx+C76PR6sbmrCxrhlsdZ2dKByYCB83xVN2Hh2PTAKwAWUDLrR7KpB9Z0/VExHUn3XDxE6f5GmMqClJfycVFZi54oVOEWvVClEREQasaXNKGaOE8vm3GrdfFq7/zK11OX5PKz8h3pM7nwUJZs2YWpnZ0LAtvHsekz1daLk2U34x/uewuHzr4Gsq4U8uxayrhYjFy0c6+rMJx1J0qzYQq2lSkREFI8tbUbIt/XJ50to2UFra+bHRR/T1xcOTKTUdu58k+Kqtcj19QGTJwODg4nb0pQlIRluNAHu4VD/aiEBeXattnLqmUiYiIjIBGxpM0I+48S0TAxIbr1auXLsMcBYwKbl3Pm2Qqm1yAmRGLBlKMvKJ4O4v2znWDJchVxryUo+zJBTjYiIyEYYtBkhn3FiWrobk4O6++9PfUw2584nKa7CQP2Elj4lkZUF4PNh5ZNBlD7xEu4/ohsYH9J+3kMuNLtqtO9PRERU5Bi0GSGfcWKZAj6loC6fMuVLqaUuXcAW1deHA1//Ov6+7Ufh1jWRxTlHoTzTk4iIyMYYtBlBqfVJ6zixTAFfLkldjV64PdJSt/KJ91HathG9Hm3B1IThYax9MDXBbVrxudSIiIgchEGbEfIZJ5Yp4NPaYiYiTVd6LdyeQfyYtNVXpuZQU2t7U0t8myAUPoBqLjUiIiIH4OxRo+Q6W1EhL1jC7FGl2Z7JSkqARx4p2GzJlU8Gcf+k7vDkASCWliM+h9qEjz/GUfv2pTxWcWUDiXCgFpdrjYEaERE5HYM2K0oX8MUHdcnpPYBwq5yBLWtjaTkAPBtQTcuxsb4+IadadNH2+DVAFVcwOORiaxoREZECdo8Wo+hsTymBRx/NPV1HlnJJyxG1sb4eV65ahV6vFyEh0Ov14spVq8YCOwmIfaUM2IiIiFSwpa3YZdsNm0vi3oj2UE92aTmSJLe+xcaqsQuUiIgoIwZtTpLjSg2xLtGKYdV9UowC4kApZPnI2OoG+0oghIAsH0HJhwzUiIiIssGgzUlUEvcOXv9teN014eWj9o8FVthfAowLAUfI7PKoaRmX5vMBLQ3ARdm3+BERETkRx7QVk3wXoVfJ8XbE4Ltj49QmjUIeNhL7HWUaA7Zs0nJoWaornXyfByIioiLElrZike8i9EC4RSu6PmkcxbQbmeSTliPdUl2Z6qLH80BERFSE2NJmddFWpRUrcl+EPkohca9i2g0NXHs/A3l2LWRdLUYuWpjd2DQj12YlIiKyKQZtVhbfjahG47JWK58MotRdg69edwN6PZG0G0d5ceWNqxJndGpxyIV/3uPOvJ8aI9dmJSIisil2jxZStuk2tCwOnybQGUuEOxxLgLtxWT02LssySIuSgNhfiqtGTsRXZnbndgxAeVWHbNZmVQpitS7vRUREVKTY0lYouQy+z9R6lCbQSUmEq/WVHg0nuUUIwEclsd9L9rpx9f9NR+j8Rfmn6TBybVYiIiKbYktboeQy+F6tVQkIBzppWupyToQrgND5i7J/XLaMWpuViIjIphi0FUouY7HUuhEVWqUSukL3lwAVozkVs+TDPMaqFUquAR8REVERY9BWKLmMxdLQqrTyySDuL3kLOGJ0LJ/apNwCNhxyodlVk9tjiYiIyFCmjGkTQtwqhHhXCPFa5Hauyn7nCCF2CiHeFkJ8u9Dl1JXaWKxzz02fKDa6OHwoFP7Z2BieCfrESxCbArj/8O5wkJZFAlyMRn4mj1njYu1ERESWZWZL211Syja1O4UQJQB+DOBzAHYDeEUI8bSUckehCqgrpVazc88FHnkkq0Sx0QkGWY9Xk8DVf5/OoIyIiKhIWXn26GkA3pZS9kgpPwHQCeBLJpcpP8mtZs88k3Wi2FwnGDT+KoB1N57OpZ+IiIiKlJlB2zeFEG8IIR4SQhyhcP9xAP4W9/fuyDZry2ZdzCwnJ6x8MojRiuGsi9TwzLN48Ce3577WJxEREZlOSCmNObAQfgBHK9zVAuBlAHsRHln1PQDHSCmvSHr8xQA+L6Vsivx9KYDTpJTXqJyvGUAzAHi93nmdnZ16VUXR0NAQysvLE7Z5/H5Ma2tDyfBYYDXqdmPnqlUYUFh14IzlyzE+GEzZfsjrxctJ5b/7LxPx6xMPZG5lkwA+FsCnJcDEEbgGP4P3/+3L8PzfHk3n8fj9qOnogHtgAMMeD3qamhTLrlR/p3By3QHWn/V3bv2dXHeA9Tey/nV1dduklPMz7iilNPUGoBrAdoXtCwH8Ie7v7wD4jpZjzps3Txqtq6srdWNVlZThtqzEW1WV8kE2bJCyrCxx37Ky8PYkJY+/KNHVpXx7tktiU5csefxFefUv96SeRwjlcgmRc3nebGkJ10uI8E+FfexK8bV3ENa/y+wimMrJ9Xdy3aVk/Y2sP4CtUkN8Y9bs0WPi/rwAwHaF3V4BcKIQYqoQ4jMAlgN4uhDly1m2udiyWBlg9EiVbtHIBIO0C7drXetTLQHwddcldvmuXIlpbW3sbiUiIiogs8a03S6E+IsQ4g0AdQBuAAAhxLFCiGcAQEo5AuCbAP4AoBvAY1LKN00qrza5LISukNJDiVrS25JBd+YZoVqXflILLgcHEwO09esTuoABZJxAQURERPkxJWiTUl4qpZwppZwlpTxfSvl+ZPt7Uspz4/Z7Rkp5kpTyBCml9ReXzDUXWxrRnGyjRw6H86zF05oMV2uLntZF19XGQWZaK5WIiIhyxhUR9KRjLrbYklSHIzG0DgEQ4Ra2ZleN9rxrWpZ+Ulo2Kxtagz4iIiLKGoM2vSUHR9XVWS0UnzF5riu8esHIRQv1K3OUUtA5NBTuHk0ikbQIg1J3KxEREenGysl17SHLyQlakueqTkrQQ/IYu3vuUezyffdLX9I0gYKIiIj0wZY2o2W5ULyWgExtUoIhVBatf/u44zCltrZw5SAiInI4trQZTevMzYiMAZnWyQd60jjDlYiIiIzDoM1oWeRiAxAOyA4lvSwhADI8lu3qg9O46DsREZEDsXu0ELTM3IxYd6EXeBJoHwrPHi35MMtZokRERGRLDNosaN2FXqwDgzQiIiIaw+5RC3no+nXoO+pohFwu9B11NB66fp3ZRSIiIiKLYNBmEQ9dvw6XrL8RVXuDcEmJqr1BXLL+RgZuREREBIBBm2Wc7bsNE5LW85wwPIyzfbeZVCIiIiKyEgZtOfL4/TmvJ6rk+MGBrLYTERGRs3AiQi58PkxrawOiLWMa1hNNlrC+qAR2HeVB9UAwZb+/VXhQpVe5iYiIqGixpS0XLS0oSerKjK0nqkF0fdHRycPhV6AEWH1lEw64ExPrHnC78WzjGp0KTURERMWMLW25yHI90WRK64turK8HAKzt6EDlwAD+VuHBs41rcMXdK/MqKhEREdkDW9pyobJuqOr2iJVPBlH6xEsYrVBeX3RjfT2mdnaixL8JVR/scV7A5vPpOk6QiIjIThi05aK1FaNJXZnp1hMFkrpERfrDF3RBeKvw+cLjAvv6ACnHxgkycCMiIgLAoC03jY3YuWqV5vVEVz4ZxP2TulO6RBWZsSC8FbS0hMcFxstinCAREZHdcUxbjgbq63HK97+ven/C7NDDkRIeN/j9sfFr/Ud5sLqpCY/NPc+564zmOU6QiIjI7tjSpoeksVgPXb8ucXaoQsD2QFsbqoPh1Q+qB4L4+d13YWS4x5kBG5DzOEEiIiKnYNCWL4WxWJesvxENz/+36kPWdnSkrH7g+K7A1tbwuMB4GcYJEhEROQmDtnwpjMWaMDyMtR0dqg+pHFBZ5cDJXYGNjeFxgRrHCRIRETkNx7TlSfb3K04GVQ3MAPSrrH7g+K7AxkYGaURERCrY0pan/gqP8nZP0vYQAAmU7HVjU8MadgUSERFRVhi05ek7KstPrW5qGtswClz99+mQdbUYuWhhOGkuuwKJiIgoC+wezdNjnz0PWDW2/FS/J5y+I7osFQ65cPXBaamzQtkVSERERFlg0JanZlcN7l+0bCxIA8a6Qgfdzs27RkRERLpi0JandRd6gSeB9qFwIt2SDxmoERERkf4YtOlg3YVerAODNCIiIjKOKUGbEOIXAKZF/jwcwN+llHMU9usFsB/AKIARKeX8ghWSiIiIyEJMCdqklJdEfxdC/BDAR2l2r5NS7jW+VERERETWZWr3qBBCAPgKgKVmloOIiIjI6szO07YYQFBK+b8q90sA/y2E2CaEaC5guYiIiIgsRUgpjTmwEH4ARyvc1SKl/HVkn/sBvC2l/KHKMY6VUr4nhPAA+COAa6SUz6ns2wygGQC8Xu+8zs5OPaqhamhoCOXl5Yaew8qcXH8n1x1g/Vl/59bfyXUHWH8j619XV7dNy7h9w4K2jCcWohTAuwDmSSl3a9j/VgBDUsq2TPvOnz9fbt26Nf9CphEIBFBbW2voOazMyfV3ct0B1p/1d279nVz3/9fevcdIVZ5xHP/+WBWziheKJXjh1lBSbRTRWolKjSWtEhVLjaCbirYpUWqUmqalIbGmra3amxoLBi8RLRW8oBK1FqMiVSteEAQFFCle6VK0UeMKVnj6x3k3GbY7C0N35uyZ+X2SyZ555uzM8+w7Z+bZ98yZA66/mvVL2qmmLc/do2OA1eUaNkl7SerTvgx8A1hZw/zMzMzMeow8m7aJwB2lAUkHSnooXe0PPClpOfAs8GBEPFzjHM3MzMx6hNyOHo2I8zqJvQuMTcvrgCNqnJaZmZlZj5T30aNmZmZmthPctJmZmZkVgJs2MzMzswJw02ZmZmZWAG7aKjVnDgwezNdOOgkGD86um5mZmVVZruceLZw5c2DyZGhrQwBvvJFdB2hpyTMzMzMzq3OeaavE9OnQ1rZ9rK0ti5uZmZlVkZu2Srz5ZmVxMzMzs27ipq0SAwdWFjczMzPrJm7aKnHFFdDcvH2suTmLm5mZmVWRm7ZKtLTArFkwaBAhwaBB2XUfhGBmZmZV5qatUi0tsH49Tzz2GKxf74bNzMzMasJf+bGLrlnRhzGb/s7Wvltoer83k3sNZcb4/nmnZWZmZnXKM227YMr8Vu4f9jFb+22BXrC13xZmNq9hyvzWvFMzMzOzOuWmrUJT5rcyc99VsOe27W/Ycxuztq3LJykzMzOre27aKjBlfiszm9dAU+e3b+27pbYJmZmZWcNw01aBWdvW/e8MW4mm93vXMBszMzNrJG7aKtDlTNrmXkzuNbR2yZiZmVlDcdNWgbIzaVvhwrbhPnrUzMzMqsZNWwUm9xoKmzv8yTb34sIPvuSGzczMzKrKTVsFZozvz4Vtw2na1Bu2QdOm3p5hMzMzs5rwl+tWaMb4/sygP4sWLeLEM0flnY6ZmZk1CM+0mZmZmRWAmzYzMzOzAnDTZmZmZlYAbtrMzMzMCsBNm5mZmVkBuGkzMzMzKwA3bWZmZmYF4KbNzMzMrADctJmZmZkVgJs2MzMzswJQROSdQ7eT9C/gjSo/TD9gU5Ufoydr5PobuXZw/a6/cetv5NrB9Vez/kERccCOVqrLpq0WJD0fEUfnnUdeGrn+Rq4dXL/rb9z6G7l2cP09oX7vHjUzMzMrADdtZmZmZgXgpm3Xzco7gZw1cv2NXDu4ftffuBq5dnD9udfvz7SZmZmZFYBn2szMzMwKwE1bhSSdLGmNpLWSpuWdT7VJOkTS45JWSXpZ0iUpfrmkdyQtS5exeedaLZLWS1qR6nw+xfpKekTSa+nn/nnnWQ2ShpeM8TJJH0qaWs/jL+kWSRslrSyJdTreylyXXg9ekjQyv8z/f2Vq/42k1am+eyXtl+KDJX1S8hy4Ib/Mu0eZ+ss+1yX9NI39GknfzCfr7lOm/nklta+XtCzF62r8u3iv61nbfkT4spMXoAl4HRgK7AEsBw7NO68q1zwAGJmW+wCvAocClwM/yju/Gv0N1gP9OsSuBqal5WnAVXnnWYO/QxPwT2BQPY8/MBoYCazc0XgDY4G/AAKOBZbknX8Vav8GsFtavqqk9sGl69XDpUz9nT7X0+vgcqA3MCS9NzTlXUN319/h9t8Bl9Xj+HfxXtejtn3PtFXmGGBtRKyLiE+BucC4nHOqqojYEBFL0/JHwCrgoHyz6hHGAbPT8mzgjBxzqZWvA69HRLW/uDpXEbEYeL9DuNx4jwNui8wzwH6SBtQm0+7XWe0RsTAiPktXnwEOrnliNVJm7MsZB8yNiC0R8Q9gLdl7RGF1Vb8kAWcBd9Q0qRrp4r2uR237btoqcxDwVsn1t2mgBkbSYOBIYEkKXZSmhW+p192DSQALJb0gaXKK9Y+IDZBt7MDnc8uudiay/Qt2o4w/lB/vRntN+C7Z7EK7IZJelPSEpBPySqoGOnuuN9rYnwC0RsRrJbG6HP8O73U9att301YZdRJriMNvJe0N3ANMjYgPgZnAF4ARwAayafN6dVxEjAROAX4gaXTeCdWapD2A04G7UqiRxr8rDfOaIGk68BkwJ4U2AAMj4kjgUuDPkvbJK78qKvdcb5ixT85m+3/a6nL8O3mvK7tqJ7Gqj7+btsq8DRxScv1g4N2ccqkZSbuTPYnnRMR8gIhojYitEbENuJGC7xboSkS8m35uBO4lq7W1fSo8/dyYX4Y1cQqwNCJaobHGPyk33g3xmiBpEnAq0BLpAz1pt+B7afkFss90fTG/LKuji+d6Q4w9gKTdgPHAvPZYPY5/Z+919LBt301bZZ4DhkkakmYeJgILcs6pqtLnGG4GVkXE70vipfvuvwWs7Pi79UDSXpL6tC+TfSh7Jdm4T0qrTQLuzyfDmtnuv+xGGf8S5cZ7AXBuOpLsWOCD9l0p9ULSycBPgNMjoq0kfoCkprQ8FBgGrMsny+rp4rm+AJgoqbekIWT1P1vr/GpkDLA6It5uD9Tb+Jd7r6Onbft5H7FRtAvZESOvkv1XMT3vfGpQ7/FkU74vAcvSZSxwO7AixRcAA/LOtUr1DyU7Qmw58HL7mAOfAx4FXks/++adaxX/Bs3Ae8C+JbG6HX+y5nQD8B+y/6a/V268yXaR/DG9HqwAjs47/yrUvpbsszvt2/8Nad1vp21iObAUOC3v/KtUf9nnOjA9jf0a4JS8869G/Sl+K3BBh3Xravy7eK/rUdu+z4hgZmZmVgDePWpmZmZWAG7azMzMzArATZuZmZlZAbhpMzMzMysAN21mZmZmBeCmzcwKTdJgSbl/T5ykEZLGllw/XdK0PHMys/rips3MrIP0DfCVGkH2vU4ARMSCiLiy+7Iys0bnps3MCkXSpZJWpsvUFN5N0ux0Uu+7JTWnda+U9EqK/zbFDpB0j6Tn0uW4FL9c0ixJC4HbJC2RdFjJ4y6SdJSkYyQ9nU6U/bSk4ekMKT8HJkhaJmmCpPMkXZ9+d5CkR1Mej0oamOK3Srou3c86SWem+ABJi9N9raynk3Gb2a5z02ZmhSHpKOB84KvAscD3gf2B4cCsiDgc+BCYIqkv2WmHDkvxX6a7uRb4Q0R8hexb3W8qeYijgHERcQ4wFzgrPe4A4MDIzrG4Ghgd2YmyLwN+FRGfpuV5ETEiIuaxveuB21Iec4DrSm4bQPZt7KcC7TNz5wB/jYgRwBFk385uZg1uV3YBmJnl5Xjg3oj4GEDSfOAE4K2IeCqt8yfgYuAaYDNwk6QHgQfS7WOAQ7NTDQKwT/v5ZYEFEfFJWr4TeAT4GVnzdleK7wvMljSM7LQ3u+9E3qPITrgN2WmRri657b7ITkb+iqT+KfYccEs6gfV9EeGmzcw802ZmhaIy8Y7n44uI+Aw4BrgHOAN4ON3WCxiVZsRGRMRBEfFRuu3jkjt4B3hP0uHABLKZN4BfAI9HxJeB04A9d6GO0ny3lCwrPfZiYDTwDnC7pHN34THMrM64aTOzIlkMnCGpWdJeZLs//wYMlDQqrXM28KSkvclOcv8QMJXsQAGAhcBF7XcoaQTlzQV+nO5nRYrtS9ZMAZxXsu5HQB869zQwMS23AE92VaSkQcDGiLgRuBkY2dX6ZtYY3LSZWWFExFLgVuBZYAnZ59H+DawCJkl6CegLzCRroB5IsSeAH6a7uRg4Oh0U8ApwQRcPeTdZs3VnSexq4NeSngKaSuKPk+12XSZpQof7uRg4P+XyHeCSHZR6IrBM0otkn7u7dgfrm1kDUETHvQpmZmZm1tN4ps3MzMysANy0mZmZmRWAmzYzMzOzAnDTZmZmZlYAbtrMzMzMCsBNm5mZmVkBuGkzMzMzKwA3bWZmZmYF8F8lMaye6lnW5wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import cvxpy as cp\n", "X_uncensored = X_ordered[:M, :]\n", "c = cp.Variable(shape=n)\n", "objective = cp.Minimize(cp.sum_squares(X_uncensored*c - y_ordered[:M]))\n", "constraints = [ X_ordered[M:,:]*c >= D]\n", "prob = cp.Problem(objective, constraints)\n", "result = prob.solve()\n", "\n", "c_cvx = np.array(c.value).flatten()\n", "fit_cvx = X_ordered.dot(c_cvx)\n", "plot_fit(fit_cvx, 'CVX fit')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Qualitatively, this already looks better than before as it no longer predicts inconsistent values with respect to the censored portion of the data. But does it do a good job of actually finding coefficients $c$ that are close to our original data?\n", "\n", "We'll use a simple Euclidean distance $\\|c_\\mbox{true} - c\\|_2$ to compare:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "norm(c_true - c_cvx): 1.49\n", "norm(c_true - c_ols_uncensored): 2.23\n" ] } ], "source": [ "print(\"norm(c_true - c_cvx): {:.2f}\".format(np.linalg.norm((c_true - c_cvx))))\n", "print(\"norm(c_true - c_ols_uncensored): {:.2f}\".format(np.linalg.norm((c_true - c_ols_uncensored))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusion\n", "Fitting censored data to a parametric distribution can be challenging as the MLE solution is often not analytically tractable. However, many MLEs can be converted into a convex optimization problems as show above. With the advent of simple-to-use and robust numerical packages, we can now solve these problems easily while taking into account the entirety of our information set by enforcing consistency conditions on various portions of the data." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" } }, "nbformat": 4, "nbformat_minor": 1 }