{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# M-Estimators for Robust Linear Modeling" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from statsmodels.compat import lmap\n", "import numpy as np\n", "from scipy import stats\n", "import matplotlib.pyplot as plt\n", "\n", "import statsmodels.api as sm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* An M-estimator minimizes the function \n", "\n", "$$Q(e_i, \\rho) = \\sum_i~\\rho \\left (\\frac{e_i}{s}\\right )$$\n", "\n", "where $\\rho$ is a symmetric function of the residuals \n", "\n", "* The effect of $\\rho$ is to reduce the influence of outliers\n", "* $s$ is an estimate of scale. \n", "* The robust estimates $\\hat{\\beta}$ are computed by the iteratively re-weighted least squares algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* We have several choices available for the weighting functions to be used" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "norms = sm.robust.norms" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def plot_weights(support, weights_func, xlabels, xticks):\n", " fig = plt.figure(figsize=(12,8))\n", " ax = fig.add_subplot(111)\n", " ax.plot(support, weights_func(support))\n", " ax.set_xticks(xticks)\n", " ax.set_xticklabels(xlabels, fontsize=16)\n", " ax.set_ylim(-.1, 1.1)\n", " return ax" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Andrew's Wave" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function weights in module statsmodels.robust.norms:\n", "\n", "weights(self, z)\n", " Andrew's wave weighting function for the IRLS algorithm\n", " \n", " The psi function scaled by z\n", " \n", " Parameters\n", " ----------\n", " z : array_like\n", " 1d array\n", " \n", " Returns\n", " -------\n", " weights : ndarray\n", " weights(z) = sin(z/a)/(z/a) for \\|z\\| <= a*pi\n", " \n", " weights(z) = 0 for \\|z\\| > a*pi\n", "\n" ] } ], "source": [ "help(norms.AndrewWave.weights)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAHYCAYAAABHt465AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3iV5eH/8c99TiZZEEhYCYQQVgZLQHECLgRlqnXUWuuou9atgAtwV+us2lpHW7WKTAXFgYpFZcjIIiSEvZIQCCE7Oc/vD/32RxElQJL7jPfrurgk4SS8LxXy4fCc+zGO4wgAAAAINC7bAQAAAIANDGEAAAAEJIYwAAAAAhJDGAAAAAGJIQwAAICAxBAGAABAQAqy9RO3a9fOSUpKsvXTAwAAIECsWLGixHGcuIPfb20IJyUlafny5bZ+egAAAAQIY8ymQ72fSyMAAAAQkBjCAAAACEgMYQAAAAQkhjAAAAACEkMYAAAAAYkhDAAAgIDEEAYAAEBAYggDAAAgIDGEAQAAEJAYwgAAAAhIDGEAAAAEJIYwAAAAAhJDGAAAAAGJIQwAAICAxBAGAABAQGIIAwAAICAxhAEAABCQGMIAAAAISAxhAAAABCSGMAAAAAISQxgAAAABiSEMAACAgMQQBgAAQEBiCAMAACAgMYQBAAAQkBjCAAAACEgMYQAAAAQkhjAAAAACEkMYAAAAAYkhDAAAgIB02CFsjPm7MabIGJP1Mz9ujDHPGmMKjDFrjDEDmz4TAAAAaFqNeUb4dUkjf+HHz5HU48dv10j6y7FnAQAAAM3rsEPYcZyvJJX+wkPGSnrT+cG3klobYzo2VSAAAADQHIKa4HN0lrTlgLe3/vi+HU3wuQEgIFTXNai0ova/3/ZU/vjPilrVeZxDfkx4sFttIkIU2ypEbSKCFRsR8sO3ViEKcvMSEAA4nKYYwuYQ7zvk79rGmGv0w+UT6tKlSxP81ADgW0orapW/q1zriytUULRf64v3q6Bov7btrTrk442Rgl0/HbWOHNU1HHogB7uNuraNUEpcpLrHRyglPlLd4yKVEh+pViFN8ds+APiHpvgdcaukxAPeTpC0/VAPdBznFUmvSNKgQYMO/Ts4APgJx3G0oaRCyzaWatnGPVq2sVSbdlf+98fDgl3qHhep47q20YWDEtU+OvSHZ3gjQtSmVYjaRoQoOjxYbtehnm+Qaus92ltZq9L/Pntcp9LKWm3bU6X1xfu1rqhcn+TuUsOPzyi7jJTWKUaDktpoSFKsBiXFKi4qtEX+XQCAN2qKITxX0o3GmHckHS+pzHEcLosAEJBK9tfos9xd+iKvWMs2lqpkf60kKTYiRIO6ttElQ7qod8dodY+LUKeYcLl+ZuQ2RkiQS/HRYYqPDvvZx9TWe7S5tEIFRRXK3l6m5Rv36O2lm/XafzZKkrq1i9CQpFid3idep/aMU1iw+6h7AMDXGMf55SdmjTFvSxomqZ2kXZLulxQsSY7jvGSMMZKe1w8nS1RKusJxnOWH+4kHDRrkLF9+2IcBgNfbUlqpj7N3amH2Li3fVCqPI3VuHa7ju8VqcLdYDU6KVfe4CP3w26V9tfUeZW0v0/KNpVq6YY++27Bb5dX1Cg9269Se7XR2WgeN6B2v1q1CbKcCQJMwxqxwHGfQT95/uCHcXBjCAHxZUXm13l+xTXNWbdPaneWSpN4donR2WgedldZeqR2jvWb4Hk5dg0ffFZb+MOZzdmrXvhq5XUYnJMdqwoAEjcroqPAQnikG4LsYwgBwjBo8jr5aV6x3lm3WZ7lFqvc4Oq5rG52T3kFnpXZQl7atbCceM4/HUea2Mn2cvVPzM3do4+5KRYUFaVz/zvrV4ESld46xnQgAR4whDABHafveKr2zdLPeW7FVO8qq1TYiROcfl6ALByeqe1yk7bxm4ziOvttQqn8v26L5mTtUU+9ReudoXTS4i8YN6KzIUE6gAOAbGMIAcITW7SrXS1+u19xV29XgODqtZ5wuGpyoEb3bKyQosM7pLaus0+xV2/T20s1au7NcMeHBunxoV11+YpLaRnLyBADvxhAGgEZasWmP/vLFen2au0vhwW5dNCRRV57cTQltfP/Sh2PlOI5Wbtmrl79cr4+zdyks2KVfDUrUVackKzGWfz8AvBNDGAB+geM4WpxfoucXFWjphlK1bhWs356YpMuHJqlNBKcnHEpBUble/rJQs1dtk8eRxvTrpBuGd1dKfJTtNAD4HwxhAPgZa7bu1SPz1+qbwt3qFBOmq05J1kVDErkLWyPtKKvSq4s36K2lm1Vd16BfDU7ULWf0VPtfON8YAFoSQxgADrJpd4WeXLhO81ZvV2xEiG4ekaJLju8acNf/NpXd+2v03OcF+td3m+R2GV11crJ+f1qyosKCbacBCHAMYQD40YGDLcjl0lWndNM1pzLYmsrm3ZV6cmGe5v74B4ybRqToUv6AAcAihjCAgNfgcfSv7zbpiY/zVFnboAsHJeqPZ/T4xVsU4+hlbi3TIwtytWT9bnVrF6GpY9N1co92trMABCCGMICAlrWtTPfOytSarWU6OaWdHhiTyou6WoDjOPpiXbEemJutTbsrNbZ/J00a3UfxUfzhA0DL+bkhzCtBAPi18uo6/WnhOr35zUbFRoTqmYv6a0y/Tj5z+2NfZ4zR8F7xGnpLW734xXq99MV6fb62SHeO7K1LhnSR28V/BwD28IwwAL/kOI7mZ+7Ug/OyVby/Rr8+vqtuP7uXYsK5Dtim9cX7NWV2lpas361+ia01fVw6t20G0Oy4NAJAwNhTUatJszM1P3On0jpFa/r4DPVPbG07Cz9yHEdzVm3XtA9ztLeyTjeN6KEbhndXkJsX0wFoHlwaASAgfL52l+56P1N7K2t158heuuaUZAaWlzHGaNyAzhreK173zc3S05+u0+d5RXrqwn7qHhdpOw9AAOGrAwC/UFFTr3tmZup3ry9XbKsQzb7hJF0/LIUR7MViWgXrmYsG6PlLBmjT7gqNfnax3liyUR6Pnb+pBBB4eEYYgM9bvrFUt767Wlv2VOr3pyXr1jN7KjTIbTsLjXRu304anBSru95fo/vnZuuTnF164oK+6hgTbjsNgJ/jqRIAPqvB4+iZT/N14cvfyJGjf18zVPec04cR7IPaR4fptd8O1sPjM/T95j06++mv9EnOLttZAPwcQxiAT9q9v0a/fW2pnv50ncb276wFfzhVQ7rF2s7CMTDG6JLju2j+zaeoa9sIXf3mcj0yP1d1DR7baQD8FJdGAPA5yzaW6qa3Vqq0slaPTsjQrwYnci6wH0lqF6H3rh2q6R/m6uWvCvX95j167uKB6hDDTTgANC2eEQbgMxzH0StfrddFr3yr0GCXZl1/oi4a0oUR7IfCgt2aOi5dz1zUX9nb92nUs4u1OL/YdhYAP8MQBuATyqrqdPWbK/Tw/LU6K7W95t10stI6cSMGfze2f2fNvfFktYsM0W/+vlRPf7JODZwqAaCJMIQBeL3C4v0a/8J/9EVeke47N1UvXjpQ0WHcIS5QpMRHavYNJ2nCgAQ981m+fv+P5dpfU287C4AfYAgD8GqL84s17oX/aG9Vnd6+5gT97uRuXAoRgFqFBOnJC/rqobFpWpRXrIkvLtGW0krbWQB8HEMYgFdyHEdvLNmo3762TB1jwjXnhpM0OIlTIQKZMUa/GZqk168YrB1lVRr7wn+0dEOp7SwAPowhDMDr1DV4NGl2lu6fm63hveL1/vUnKjG2le0seIlTesRp9g0nqXV4sC7927f697LNtpMA+CiGMACvsqeiVpe9+p3e+m6zrhvWXa9cdpwiQznpEf8rOS5Ss64/SSckt9Vd72dq6gc5vIgOwBHjqwsAr7GltFKX/32ptu6t0tO/6qfxAxJsJ8GLxbQK1mu/HaxpH+bq1a83aHNppZ67eIDCgrmzIIDG4RlhAF4ha1uZJvxliXZX1Oqtq45nBKNRgtwuPTAmTQ+OSdOnubt06d++097KWttZAHwEQxiAdV/nl+iiV75VsMtoxrVDNYgXxeEIXX5ikl64ZKAyt5Zp4l+WaOseTpQAcHgMYQBWzVm1TVe8vlQJbcI18/qT1KN9lO0k+KhRGR315pVDVFReo4l/WaLcHftsJwHwcgxhANb89atC/eGdVRrYpY3+/fuh6hATZjsJPu6E5LZ679qhMjK68KVv9M363baTAHgxhjCAFuc4jqZ/mKPp83M1KqOD3vjdEMWEc6c4NI3eHaI18/oT1T4mTJf/fakWZO6wnQTASzGEAbQoj8fR5NlZ+uviDfrN0K567uKBvMofTa5T63DNuHao0jtH64a3vteslVttJwHwQgxhAC2mvsGj22es1r++26xrT+uuB8ekye3idsloHq1bhegfVx6v47u11a3vrtbbS7nxBoD/xRAG0CLqGjz6w79Xaeb323TrmT1118heMoYRjOYVERqk164YrNN6xumemZn6+9cbbCcB8CIMYQDNrrquQdf9c4U+XLND947qrZtP78EIRosJC3br5cuO09lp7fXQBzl68YsC20kAvARDGECzqqpt0NVvLtenuUWaOjZN15za3XYSAlBokFvPXzJQY/p10uMf5emphXlyHG7JDAQ6brEMoNlU1tbriteWadnGUj1+fl9dOCjRdhICWLDbpad/1V/hwW49+3mBauo9uvuc3vztBBDAGMIAmkV1XYOuemO5lm0s1dO/6q+x/TvbTgLkdhk9MiFDIUEuvfxVoYLcRrefxfXqQKBiCANoctV1P1wO8U3hbj11YT9GMLyKy2X04Jg01Xs8emHReoW43frDGT1sZwGwgCEMoEnV1Dfo+n99r8X5JXr8/L4aPyDBdhLwEy6X0fRxGaprcPT0p+sU5Da6YXiK7SwALYwhDKDJ1DV4dONbK/X52iI9PD6Da4Lh1Vwuo8cm9lV9g0dPfJynYLfhxZxAgGEIA2gS9Q0e/eGdlfokZ5ceHJOmS47vYjsJOCy3y+jJC/qpzuPo4flrFex26YqTutnOAtBCGMIAjlmDx9Gt767W/Mydmjy6jy4/Mcl2EtBoQW6X/vyr/mpocPTgvBwFu1369QldbWcBaAGcIwzgmDiOoylzsjR39XbdNbK3rjol2XYScMSC3S49e/EAnd47XlPmZGnOqm22kwC0AIYwgGPyp4Xr9NZ3m3XdsO66bhjXV8J3hQS59MKlAzUkKVa3vbtai/KKbCcBaGYMYQBH7W+LC/X8ogJdPCRRd57dy3YOcMzCgt366+WD1KtDlK775wot31hqOwlAM2IIAzgqM1Zs1bQPczUqo4OmjcvghgTwG9FhwXrjd0PUKSZcv3t9mXJ37LOdBKCZMIQBHLFPcnbprvfX6OSUdnr6V/3ldjGC4V/aRYbqzSuHqFVIkH7z96XatLvCdhKAZsAQBnBEvi3crRve+l7pnWP08mXHKTTIbTsJaBYJbVrpH1cOUV2DR5e9ulRF+6ptJwFoYgxhAI2Ws32frnpjubrEttLrvx2siFBOYIR/69E+Sq9fMUQl+2v0m78v1b7qOttJAJoQQxhAo2zfW6UrXl+qqLAg/ePKIWoTEWI7CWgR/RNb6+XLjlNB0X5d/8/vVVvvsZ0EoIkwhAEcVllVna54bZkqaxr02hWD1TEm3HYS0KJO6RGnRyf21dcFJbp75ho5jmM7CUAT4O81Afyi2nqPrvvnChWW7NcbVwxR7w7RtpMAK84/LkHb91bpqU/WKaFNK916Zk/bSQCOEUMYwM9yHEd3v79GS9bv1lMX9tOJKe1sJwFW3TQiRdv2VOnZz/KV0DpcFw5OtJ0E4BgwhAH8rKc+WaeZK7fp9rN6asLABNs5gHXGGE0bn64d+6p1z6xMtY8J02k942xnAThKXCMM4JDeWbpZz31eoIsGJ+qG4Sm2cwCvEex26cVLB6pX+yhd/88Vyt5eZjsJwFFiCAP4ia/WFWvS7Cyd1jNO08alc9c44CCRoUF67YrBigkP1hWvLdOOsirbSQCOAkMYwP8oKCrXDf/6Xj3bR+mFSwcqyM1vE8ChtI8O02tXDFFlbYOuemO5KmvrbScBOEJ8hQPwX3sqavW715crNNitv10+SJHcMAP4Rb06ROm5SwYod8c+3frv1fJ4OFYN8CUMYQCSfjgm7dp/rtDOfdV65TfHqXNrzgoGGmN4r3hNGp2qj7J36qlP1tnOAXAEeLoHgBzH0ZTZWfpuQ6meuai/BnZpYzsJ8Cm/OylJBUXlen5RgVLiIzVuQGfbSQAagWeEAejVrzfo38u36KYRKRrbny/gwJEyxujBMek6ITlWd76/Ris27bGdBKARGMJAgPssd5emz8/VqIwO+uMZ3CkLOFohQS795dLj1CkmTL//x3Jt3VNpOwnAYTRqCBtjRhpj8owxBcaYuw/x412MMYuMMSuNMWuMMaOaPhVAU8vbWa6b316p9E4x+tMF/eVycUwacCzaRITob5cPVk29R1e9sVwVNZwkAXizww5hY4xb0guSzpGUKuliY0zqQQ+bLOldx3EGSLpI0otNHQqgae2trNXVby5XRGiQ/vqbQQoPcdtOAvxCSnykXrx0oNbtKtft762W43CSBOCtGvOM8BBJBY7jFDqOUyvpHUljD3qMIyn6x+/HSNredIkAmlqDx9FNb6/UzrJqvXTZceoQE2Y7CfArp/SI072j+mhB1k69+MV62zkAfkZjTo3oLGnLAW9vlXT8QY95QNJCY8xNkiIkndEkdQCaxRMf52lxfokenZDBCRFAM7ny5G7K2lamJxfmKbVjtIb3jredBOAgjXlG+FAXDR789zwXS3rdcZwESaMk/cMY85PPbYy5xhiz3BizvLi4+MhrARyzeau366Uv1+vS47vooiFdbOcAfssYo0cm9FWfDtG6+Z2V2lBSYTsJwEEaM4S3Sko84O0E/fTShyslvStJjuN8IylMUruDP5HjOK84jjPIcZxBcXFxR1cM4KjlbN+nO2es0aCubXT/eWm2cwC/Fx7i1iu/OU7BbpeueXO59vPiOcCrNGYIL5PUwxjTzRgToh9eDDf3oMdslnS6JBlj+uiHIcxTvoAX2VNRq9//c7liwoP14q8HKiSI0xOBlpDQppWev2SACksqdOu/V3EbZsCLHPYroeM49ZJulPSxpFz9cDpEtjHmIWPMmB8fdpukq40xqyW9Lem3Di+TBbxGfYNHN729UrvKavSXXw9UfBQvjgNa0ond2+neUX20MGeXnl9UYDsHwI8adYtlx3HmS5p/0PvuO+D7OZJOato0AE3liY/z9HVBiR4/v68G8OI4wIrfnZSk7G1leuqTdUrvHK0RvdvbTgICHn83Cvi5j7J26uWvCvXrE7rowkGJh/8AAM3CGKOHJ2QorVO0bnlnlbaUcuc5wDaGMODHNpRU6I73VqtfYmtNOffg++AAaGlhwW795dLjJEnX/WuFqusaLBcBgY0hDPipqtoGXffPFQpyG7146UCFBnHnOMAbdGnbSk9d2F9Z2/bpwXnZtnOAgMYQBvyQ4ziaNDtTebvK9eeLBqhz63DbSQAOcEZqe10/rLveXrpF7y3fcvgPANAsGMKAH3p76RbN/H6b/nB6D53WkzO7AW9065k9dWL3tpo8O0s52/fZzgECEkMY8DNrtu7VA3OzdWrPON08ooftHAA/I8jt0rMXD1DrVsG67l8rVFZVZzsJCDgMYcCP7K2s1XX//F5xUaH686/6y+U61B3SAXiLdpGheuGSgdq2p0q3v7daHMEPtCyGMOAnPB5Hf/z3KhWX1+jFSwcqNiLEdhKARhiUFKt7R/XRJzm79MpXhbZzgIDCEAb8xF8XF2pRXrEmn9tH/RJb284BcASuOClJozI66PGP87RiU6ntHCBgMIQBP7Bi0x49/nGeRmV00GUndLWdA+AIGWP06MS+6tQ6TDe/vUp7K2ttJwEBgSEM+Li9lbW6+e2V6tQ6TI9O7CtjuC4Y8EXRYcF64ZKBKiqv1u3vreF6YaAFMIQBH+Y4jm5/b42Kyqv1wiUDFR0WbDsJwDHom9Ba95zTR5/m7tKrX2+wnQP4PYYw4MNe/XqDPs3dpXvO6aO+CVwXDPiDK05K0pmp7fXYR2u1aste2zmAX2MIAz5q1Za9euyjtToztb2uOCnJdg6AJmKM0RPn91V8VJhufOt7zhcGmhFDGPBBZVV1uvGt7xUfFaYnzue6YMDftG4VoucuGaCdZdW6awbXCwPNhSEM+BjHcXT3+2u0s6xaz10yQK1bcV4w4I8GdmmjO0f20kfZO/XmN5ts5wB+iSEM+Ji3lm7Wgqyduv3sXhrYpY3tHADN6KqTkzW8V5ymz89V7o59tnMAv8MQBnxI/q5yTf0gR6f0aKdrTkm2nQOgmblcRk9c0E/RYcG66e2VqqptsJ0E+BWGMOAjqusadNPbK9UqJEh/uqCfXC6uCwYCQbvIUD11YT8VFO3X1A9zbOcAfoUhDPiIRxes1dqd5Xrygr6Kjw6znQOgBZ3aM07XnJqst77brI+ydtrOAfwGQxjwAZ+v3aXXl2zUb09M0oje7W3nALDg9rN6KaNzjO6euUY7yqps5wB+gSEMeLmifT/cbrVPx2jdfU5v2zkALAkJcunZiweott6jW95ZpQYPR6oBx4ohDHgxj8fRre+uVmVtvZ67uL/Cgt22kwBY1K1dhB4ck6bvNpTqL18U2M4BfB5DGPBif/u6UF8XlOi+c9OUEh9lOweAFzj/uASd16+Tnv40X99v3mM7B/BpDGHAS2VtK9MTH+dpZFoHXTwk0XYOAC9hjNG0cenqEB2mP7yzUuXV3IIZOFoMYcALVdU26A/vrFRsRIgemZDBLZQB/I+Y8GA9c1F/bdtTpYfmcaQacLQYwoAXenRBrtYXV+jJC/qpTQS3UAbwU4OSYnX9sBS9t2KrFmTusJ0D+CSGMOBlvsgr0hvfbNLvTuqmU3rE2c4B4MX+cEYP9U2I0T2zMrVrX7XtHMDnMIQBL1JaUas7ZqxRz/aRunNkL9s5ALxcsNulp3/VX9V1Dbr9vdXycKQacEQYwoCXcBxHd7+/RmWVdfrzrwZwVBqARukeF6lJo1O1OL9Eb36z0XYO4FMYwoCXeG/5Vi3M2aXbz+6p1E7RtnMA+JBfH99Fw3vF6ZEFa7VuV7ntHMBnMIQBL7Bpd4UemJetocltddXJybZzAPgYY4weP7+fIkODdMs7q1RT32A7CfAJDGHAsvoGj2759yq5XUZ/urCfXC6OSgNw5OKiQvXoxL7K2bFPT32yznYO4BMYwoBlL325Xis379W0cenq1Drcdg4AH3ZmantdPCRRr3xVqKUbSm3nAF6PIQxYlLWtTH/+NF/n9u2osf07284B4Acmj05VYptWuu29VdpfU287B/BqDGHAkpr6Bt327mq1iQjR1LHptnMA+ImI0CA9eUE/bd1Tpekf5trOAbwaQxiw5KlP1ilvV7ken9iXu8cBaFJDusXq6lOS9fbSzVqUV2Q7B/BaDGHAguUbS/XKV4W6eEiihveOt50DwA/demZP9WwfqbtmrNHeylrbOYBXYggDLayipl63vbdaCW3CNWl0qu0cAH4qLNitpy7sr9KKWt03J9t2DuCVGMJAC3t4fq42l1bqyR/P/ASA5pLeOUY3n95Dc1dv1wdrttvOAbwOQxhoQV+uK9a/vtusq07upuOT29rOARAArh/WXf0SW2vy7CwV7au2nQN4FYYw0ELKKut014w16hEfqdvO6mU7B0CACHK79KcL+qmqtkF3z8yU4zi2kwCvwRAGWsiD87JVvL9GT13YX2HBbts5AAJISnyk7hrZW5+vLdJ7K7bazgG8BkMYaAGf5OzSzJXbdMPwFGUkxNjOARCAfntikoZ0i9XUeTnaUVZlOwfwCgxhoJntrazVvbMy1adjtG4cnmI7B0CAcrmMnjy/n+o9ju56n0skAIkhDDS7B+Zma09FrZ68oK9CgvglB8CeLm1b6Z5RvfXVumK9u3yL7RzAOr4qA83o4+ydmr1qu24ckaK0TlwSAcC+Xx/fVSckx2rqB7natpdLJBDYGMJAM9lTUatJs7KU2jFaN3BJBAAv4XIZPXF+P3kcR3e/v4ZLJBDQGMJAM7l/brbKqmr15AX9FOzmlxoA75EY20r3jOqjxfklemcZl0ggcPHVGWgGH2Xt0NzV23XTiB5K7RRtOwcAfuLSIV10Yve2mvZBjrbuqbSdA1jBEAaa2O79NZo0K0vpnaN13bDutnMA4JBcLqPHJvaVJN3FJRIIUAxhoIk9MC9H+6rruCQCgNdLjG2le0f30X8KduutpZtt5wAtjq/SQBP6OHun5q3erptH9FDvDlwSAcD7XTKki05KaatH5q/Vdk6RQIBhCANNpKyyTpNn/3BKxLVcEgHARxhj9OiEvmrwOLp3FjfaQGBhCANNZOqHOSqtqNXj5/flkggAPiUxtpXuGtlLX+QVa+b322znAC2Gr9ZAE/gir0gzVmzVdad1V3pnbpwBwPf8ZmiSBnVtowfnZatoX7XtHKBFMISBY1ReXad7ZmYqJT5SN53OjTMA+CaXy+jx8/uqpt6jybOzuEQCAYEhDByjRxas1a591Xri/L4KDXLbzgGAo5YcF6lbz+yphTm79MGaHbZzgGbHEAaOwZKCEr313WZdeXI3DejSxnYOAByzK0/upn4JMbp/brZ276+xnQM0K4YwcJQqa+t118w1SmrbSree2ct2DgA0iSC3S4+f30/l1XV6YF6O7RygWTGEgaP0xMd52lJapcfP76fwEC6JAOA/enWI0s0jemje6u1amL3Tdg7QbBjCwFH4fvMevb5koy47oauGdIu1nQMATe7aYd3Vp2O0pszJUllVne0coFkwhIEjVFPfoLtmrFHH6DDddU5v2zkA0CyC3S49NjFDxeU1enRBru0coFk0aggbY0YaY/KMMQXGmLt/5jEXGmNyjDHZxpi3mjYT8B4vLlqv/KL9mj4+Q5GhQbZzAKDZ9E1oratOSdbbS7fom/W7becATe6wQ9gY45b0gqRzJKVKutgYk3rQY3pIukfSSY7jpEm6pRlaAevydpbrxS8KNK5/Jw3vHW87BwCa3R/P6KmubVvpnplrVFXbYDsHaDJvAZQAACAASURBVFKNeUZ4iKQCx3EKHceplfSOpLEHPeZqSS84jrNHkhzHKWraTMC+Bo+ju95fo6iwYN13XprtHABoEeEhbj0yIUMbd1fqz5+us50DNKnGDOHOkrYc8PbWH993oJ6Sehpj/mOM+dYYM7KpAgFv8fqSjVq1Za/uPy9VsREhtnMAoMWc2L2dLhqcqL8uLlTm1jLbOUCTacwQNod438H3XQyS1EPSMEkXS/qbMab1Tz6RMdcYY5YbY5YXFxcfaStgzZbSSj35cZ5G9I7XmH6dbOcAQIu7Z1QftYsM1Z3vr1Fdg8d2DtAkGjOEt0pKPODtBEnbD/GYOY7j1DmOs0FSnn4Yxv/DcZxXHMcZ5DjOoLi4uKNtBlqU4zi6Z2am3C6jaePSZcyh/mwIAP4tJjxYU8elK3fHPr3yVaHtHKBJNGYIL5PUwxjTzRgTIukiSXMPesxsScMlyRjTTj9cKsGvEviFGSu26uuCEt11Tm91ah1uOwcArDk7rYNGZ3TUM5/la33xfts5wDE77BB2HKde0o2SPpaUK+ldx3GyjTEPGWPG/PiwjyXtNsbkSFok6Q7HcThnBT6vuLxG0z7M1eCkNrp0SBfbOQBg3QNj0hQe7NY972fK4zn4SknAtzTqHGHHceY7jtPTcZzujuNM//F99zmOM/fH7zuO49zqOE6q4zgZjuO805zRQEt56IMcVdU26JEJfeVycUkEAMRFhWrSqD5aurFU7yzbcvgPALwYd5YDfsaitUWat3q7bhieopT4SNs5AOA1LhiUoKHJbfXIglwV7au2nQMcNYYwcAgVNfWaPDtLPeIjdd2w7rZzAMCrGGP08IQM1dR7dP/cbNs5wFFjCAOH8OTCPG0vq9KjEzMUEsQvEwA4WLd2EfrD6T20IGunFmbvtJ0DHBW+wgMHWbVlr15fslG/Pr6rjusaazsHALzWNacmq3eHKN03J1vl1XW2c4AjxhAGDlDX4NHd769R+6gw3Tmyl+0cAPBqwW6XHp3YV7vKq/X4R3m2c4AjxhAGDvDKV4Vau7NcD45NU1RYsO0cAPB6/RNb6/KhSfrnd5u0YlOp7RzgiDCEgR9tKKnQM5/la2RaB52d1sF2DgD4jNvP7qWO0WG6+/1M1dQ32M4BGo0hDOiH2yjfOzNToUEuPTg2zXYOAPiUyNAgTRufrvyi/XrpC24sC9/BEAb0w22UvyncrbtG9lb76DDbOQDgc0b0bq9z+3bUC4sKuP0yfAZDGAFv9/4aTZ+fq0Fd2+gSbqMMAEftvvNSFRbs0qRZmXIcbr8M78cQRsCb/mGuKmrq9fCEDG6jDADHID4qTHef00ffFpbqvRVbbecAh8UQRkD7Or9EM1du0+9P7a6e7aNs5wCAz7tocKIGdW2jh+fnavf+Gts5wC9iCCNgVdc1aNLsTCW1baUbR6TYzgEAv+ByGT0yIUMVNfWa9mGu7RzgFzGEEbCe/7xAm3ZXavr4DIUFu23nAIDf6NE+Stee1l2zVm7T1/kltnOAn8UQRkBat6tcL325XhMGdtZJKe1s5wCA37lheIq6tYvQpNmZqq7jbGF4J4YwAo7H4+iemZmKCgvS5NGptnMAwC+FBbs1fVy6Nu2u1HOf59vOAQ6JIYyA8/ayzVqxaY8mjU5VbESI7RwA8FsnprTTxIEJevnLQuXtLLedA/wEQxgBpai8Wo8uWKsTu7fVxIGdbecAgN+bNLqPosKCdO+sTHk8nC0M78IQRkCZ+kGuauo9mjYuXcZwZjAANLfYiBBNHp2qFZv26J1lW2znAP+DIYyA8eW6Ys1bvV03DEtRclyk7RwACBgTBnbW0OS2enRBrorLOVsY3oMhjIBQVdugybMzlRwXoWuHJdvOAYCAYozRtPHpqq7zaNqHObZzgP9iCCMgPPd5vraUVmn6uAyFBnFmMAC0tO5xkbpuWHfNWbVdX60rtp0DSGIIIwCs21WuV74q1MSBCRrava3tHAAIWNcN665u7SI0ZU4WZwvDKzCE4dc8Hkf3/nhm8KTRfWznAEBAO/Bs4ec/L7CdAzCE4d/eXb5Fyzft0T2j+nBmMAB4gRNT2mnCgM56+av1yt/F2cKwiyEMv1Wyv0aPLFirId1idcFxCbZzAAA/und0H7UKCdKkWVmcLQyrGMLwW9M/zFVlbb0eHs+ZwQDgTdpFhureUb21dGOpZqzYajsHAYwhDL/0n4ISzVq5Tdee1l0p8VG2cwAAB7nguEQNTmqjhxfkavd+zhaGHQxh+J3qugZNnp2lrm1b6YbhKbZzAACH4HIZPTw+Q/ur6/XIgrW2cxCgGMLwOy99uV4bSio0dWy6woI5MxgAvFWP9lG65tRkzVixVd8W7radgwDEEIZf2VBSoRcXrdd5/Trp1J5xtnMAAIdx04geSmgTrkmzMlVb77GdgwDDEIbfcBxHU2ZnKTTIpSmcGQwAPiE8xK2pY9O1vrhCf11caDsHAYYhDL8xd/V2fV1QojtG9lJ8dJjtHABAIw3vHa9z0jvo2c/ytWl3he0cBBCGMPxCWWWdpn6Qo34JMbr0+K62cwAAR+j+89IU5DK6b062HIezhdEyGMLwC08sXKvSilpNH58ht4szgwHA13SICdNtZ/XSl+uKNT9zp+0cBAiGMHzeys179K/vNuvyE5OU3jnGdg4A4Cj9ZmhXpXWK1oPzslVeXWc7BwGAIQyfVt/g0aRZWYqPCtWtZ/a0nQMAOAZBbpemj89Q8f4a/WnhOts5CAAMYfi015dsVM6OfXrgvDRFhQXbzgEAHKP+ia112Qld9eY3G5W5tcx2DvwcQxg+a0dZlZ7+ZJ2G9YrTyPQOtnMAAE3k9rN7qW1kqCbNzlSDhxfOofkwhOGzHpqXo3qPo4fGpMsYXiAHAP4iOixYk0f30ZqtZfrXd5ts58CPMYThkxblFWlB1k7dNCJFXdq2sp0DAGhiY/p10kkpbfXER3kqKq+2nQM/xRCGz6mua9B9c7LUPS5CV5+abDsHANAMjDGaOjZdNfUeTf8w13YO/BRDGD7n+c8LtKW0SlPHpSs0yG07BwDQTJLjInXtsO6as2q7vs4vsZ0DP8QQhk8pKNqvl79arwkDOuvE7u1s5wAAmtn1w7qra9tWmjInS9V1DbZz4GcYwvAZjuNo8uxMhQe7de/oPrZzAAAtICzYralj07WhpEIvf1loOwd+hiEMnzF71TZ9W1iqu87prXaRobZzAAAt5NSecTq3b0e98EWBNpZU2M6BH2EIwyeUVdZp+oe56p/YWhcP7mI7BwDQwqacm6oQt0tT5mTJcThbGE2DIQyf8PjHa1VaUavp49PlcnFmMAAEmvbRYbr9rJ5anF+iDzN32M6Bn2AIw+ut2rJXby3drN+e2E1pnWJs5wAALLlsaJLSO0froXk5Kq+us50DP8AQhlerb/Bo0qxMxUeF6tazetrOAQBY5HYZTR+XoeL9NXrqk3W2c+AHGMLwav/4dpOyt+/TfeemKTI0yHYOAMCyfomtdenxXfTGko3K2lZmOwc+jiEMr7VrX7X+tHCdTu0Zp1EZHWznAAC8xB1n91ZsRIgmz86Sx8ML53D0GMLwWlM/yFFtg0cPjUmTMbxADgDwg5jwYE0a3UertuzV28s2286BD2MIwystzi/WB2t26IZhKUpqF2E7BwDgZcb176yhyW312IK1KtlfYzsHPoohDK9TXdegKbOz1K1dhK4dlmw7BwDghYwxmjouXVV1DXp4fq7tHPgohjC8zktfrtfG3ZWaOjZdoUFu2zkAAC+VEh+pa05N1szvt+nbwt22c+CDGMLwKhtLKvTiF+t1Xr9OOrlHO9s5AAAvd+PwHkpoE67Js7NUW++xnQMfwxCG13AcR1PmZCnU7dKU0X1s5wAAfEB4iFsPjU1TQdF+/e3rQts58DEMYXiNDzN3aHF+iW47q6fio8Ns5wAAfMSI3u11dlp7PftZvraUVtrOgQ9hCMMrlFfX6aF5OUrvHK3LhibZzgEA+Jj7z0uTyxg9OC/bdgp8CEMYXuGpT9apeH+Npo/LkNvFmcEAgCPTqXW4bjmjhz7NLdLC7J22c+AjGMKwLmtbmd5YslGXHt9F/RJb284BAPioK07qpl7to/TgvBxV1tbbzoEPYAjDKo/H0eTZWYqNCNEdZ/e2nQMA8GHBbpemj0/Xtr1VeuazfNs58AEMYVj19rLNWrVlryaN7qOY8GDbOQAAHzcoKVYXDkrQq4s3KG9nue0ceLlGDWFjzEhjTJ4xpsAYc/cvPO58Y4xjjBnUdInwVyX7a/TYgrU6ITlW4/p3tp0DAPATd5/TR5FhQZoyO0uO49jOgRc77BA2xrglvSDpHEmpki42xqQe4nFRkm6W9F1TR8I/PTJ/rarqGjRtXLqM4QVyAICmERsRortH9tbSjaV6//tttnPgxRrzjPAQSQWO4xQ6jlMr6R1JYw/xuKmSHpdU3YR98FPfFe7W+99v1dWnJCslPsp2DgDAz1w4KFEDu7TWw/Nztbey1nYOvFRjhnBnSVsOeHvrj+/7L2PMAEmJjuN80IRt8FO19R5Nnp2lhDbhumlED9s5AAA/5HIZTR+fobKqOj32UZ7tHHipxgzhQ/2d9X8vuDHGuCQ9Lem2w34iY64xxiw3xiwvLi5ufCX8yqtfb1B+0X49OCZN4SFu2zkAAD/Vp2O0rjgxSW8v3awVm/bYzoEXaswQ3iop8YC3EyRtP+DtKEnpkr4wxmyUdIKkuYd6wZzjOK84jjPIcZxBcXFxR18Nn7WltFLPfLZOZ6W21+l92tvOAQD4uVvO7KkO0WGaPDtL9Q0e2znwMo0Zwssk9TDGdDPGhEi6SNLc//tBx3HKHMdp5zhOkuM4SZK+lTTGcZzlzVIMn/bgvGy5jNH9Y9JspwAAAkBkaJAeGJOq3B379PqSjbZz4GUOO4Qdx6mXdKOkjyXlSnrXcZxsY8xDxpgxzR0I/7Ewe6c+zS3SLWf0UOfW4bZzAAAB4uy0DhreK05Pf7JOO8qqbOfAizTqHGHHceY7jtPTcZzujuNM//F99zmOM/cQjx3Gs8E4WGVtvR6cl6Ne7aN0xUndbOcAAAKIMUYPjklXvcfR1A9ybOfAi3BnObSIZz7L17a9VZo+Pl3Bbv63AwC0rC5tW+mmESman7lTX+QV2c6Bl2CRoNnl7SzXq4s36MJBCRqUFGs7BwAQoK4+NVnJcRG6b062qusabOfACzCE0aw8HkeTZ2cqMixId5/Tx3YOACCAhQa5NW1sujaXVuqFRQW2c+AFGMJoVjO+36plG/fonnN6KzYixHYOACDAnZjSTuP6d9JLX67X+uL9tnNgGUMYzWZPRa0eXbBWg7q20QXHJR7+AwAAaAGTRqcqLNitKbOz5DjO4T8AfoshjGbz2EdrVVZVp2nj0+VyHeoGhQAAtLy4qFDdObK3lqzfrbmrtx/+A+C3GMJoFss3luqdZVt01cnd1LtDtO0cAAD+xyVDuqhfYmtN/SBHZZV1tnNgCUMYTa6uwaNJs7LUuXW4/nBGD9s5AAD8hNtlNH1cukoravXEwrW2c2AJQxhN7u9fb1DernI9MCZNrUKCbOcAAHBI6Z1j9NsTu+lf323Wqi17befAAoYwmtS2vVX686f5OjO1vc5MbW87BwCAX3TrWT3VPipMk2Zlqr7BYzsHLYwhjCb1wNzsH/45Js1yCQAAhxcZGqT7zktV9vZ9evObTbZz0MIYwmgyC7N36pOcXbrljB7q3Drcdg4AAI1yTnoHDesVpz8tzNPOsmrbOWhBDGE0iYqaej0wN1u92kfpdyd3s50DAECjGWP00Jh01XscPfRBtu0ctCCGMJrEs5/la3tZtaaPT1ewm/+tAAC+pUvbVrr59B6an7lTi/KKbOeghbBYcMzW7tynV7/eoIsGJ2pQUqztHAAAjsrVpyQrJT5S983JUlVtg+0ctACGMI6Jx+No0qwsRYcH666RvW3nAABw1EKCXJo2Ll1bSqv0/KJ82zloAQxhHJN/L9+iFZv26N5RfdQmIsR2DgAAx+SE5LaaODBBr3xVqPxd5bZz0MwYwjhqJftr9OiCtTq+W6wmDuxsOwcAgCZx76jeiggN0qRZWfJ4HNs5aEYMYRy16R/mqrK2XtPHZ8gYYzsHAIAm0TYyVPec01tLN5ZqxoqttnPQjBjCOCpLCko0a+U2XXtad6XER9rOAQCgSV1wXKIGJ7XRwwtyVVpRazsHzYQhjCNWU9+gybOz1LVtK90wPMV2DgAATc7lMpo+PkP7q+v18Pxc2zloJgxhHLG/fLFehSUVmjo2XWHBbts5AAA0i57to3T1qcmasWKrvi3cbTsHzYAhjCOyoaRCLy5ar/P6ddKpPeNs5wAA0KxuHtFDibHhmjQrU7X1Hts5aGIMYTSa4ziaPDtTocEuTTm3j+0cAACaXXiIWw+NSdf64gq98tV62zloYgxhNNqcVdv1n4LduvPsXoqPCrOdAwBAixjeO16jMjrouc8LtLGkwnYOmhBDGI2yt7JWUz/IUf/E1rrk+K62cwAAaFH3n5emYLdLU+ZkyXE4W9hfMITRKI8uWKu9VXV6eHyG3C7ODAYABJb20WG6c2QvLc4v0dzV223noIkwhHFYSzeU6p1lW3Tlyd2U2inadg4AAFZcenxX9UuI0dQPcrS3krOF/QFDGL+ott6je2dlqnPrcN1yRg/bOQAAWON2GT08IUN7Kuv02EdrbeegCTCE8Yte+Wq9Cor266GxaWoVEmQ7BwAAq9I6xeh3JyXp7aVbtGxjqe0cHCOGMH7WxpIKPft5gUZldNDpfdrbzgEAwCvcckZPdW4drntncrawr2MI45Acx9GUOVkKcbt0/3lptnMAAPAaEaFBemhsmvKL9uuviwtt5+AYMIRxSHNXb9fi/BLdObKX2kdzZjAAAAc6vU97nZPeQc9+lq9Nuzlb2FcxhPET/3dmcL+EGF3KmcEAABzS/50tPHk2Zwv7KoYwfuLRBWu1p7JOD0/gzGAAAH5Oh5gw3X5WTy3OL9GcVZwt7IsYwvgf3xXu/u+ZwWmdYmznAADg1S4bmsTZwj6MIYz/qqlv0D2zMpXQhjODAQBoDLfL6JEJfX+4++r8XNs5OEIMYfzXi4vWq7C4QtPGpXNmMAAAjZTaKVpXndJN7y7fqm/W77adgyPAEIYkqaCoXH/5Yr3G9OukYb3ibecAAOBTbjm9pxJjwzVpVqaq6xps56CRGMKQx+Po3plZCgt2acq5qbZzAADwOeEhbk0fl6HCkgq9uKjAdg4aiSEMvbt8i5ZuLNWk0X0UFxVqOwcAAJ90as84jevfSX/5cr3yd5XbzkEjMIQDXHF5jR6en6vju8XqwkGJtnMAAPBpk89NVURokO6dlSmPh7OFvR1DOMA99EGOqus8enhChozhzGAAAI5Fu8hQ3Tuqj5Zt3KN3lm2xnYPDYAgHsEV5RZq3ertuGJ6i7nGRtnMAAPALFxyXoBOSY/XIglwV7au2nYNfwBAOUBU19Zo8K0vd4yJ07bBk2zkAAPgNY4weHp+hmnqPHpyXYzsHv4AhHKD+tHCdtu2t0mMT+yo0yG07BwAAv5IcF6k/nN5DH2bu0Cc5u2zn4GcwhAPQqi179fqSDfr1CV00KCnWdg4AAH7pmlOT1btDlKbMzlJ5dZ3tHBwCQzjA1DV4dPf7axQfFaY7R/a2nQMAgN8Kdrv06MS+2lVercc/yrOdg0NgCAeYV74q1Nqd5Zo6Ll3RYcG2cwAA8Gv9E1vrihO76R/fbtLyjaW2c3AQhnAAKSzer2c+y9eojA46M7W97RwAAALCbWf1VOfW4bp7ZqZq6rn9sjdhCAcIj8fRPTMzFRbk0gNj0mznAAAQMCJCgzRtfLoKivbrxUXrbefgAAzhAPHu8i36bsMPt1GOjwqznQMAQEAZ3ite4/p30otfFHD7ZS/CEA4ARfuqNX1+rk5I5jbKAADYMuXcVEWGBumu99dw+2UvwRAOAPfPzVZNvUePTOjLbZQBALCkbWSoppybqu8379U/vt1kOwdiCPu9j7J2aEHWTt1yRg91axdhOwcAgIA2fkBnndKjnR7/aK227a2ynRPwGMJ+rKyyTlPmZCutU7SuPoXbKAMAYNv/3X7ZkTRpVqYch0skbGII+7Hp83NUWlGrxyb2VbCb/9QAAHiDxNhWuvPsXvoir1izV22znRPQWEd+6uv8Er27fKt+f2qy0jvH2M4BAAAHuGxoko7r2kYPzstRyf4a2zkBiyHshypr63X3zDVKbhehm0/vYTsHAAAcxO0yemxihiprGvTA3GzbOQGLIeyHnvx4nbbuqdJj5/dVWLDbdg4AADiElPgo3Xx6ij5Ys0MLs3fazglIDGE/8/3mPXptyQb9ZmhXDU6KtZ0DAAB+we9P667eHaI0ZU6WyqrqbOcEHIawH6mpb9BdM9aoY3SY7hzZ23YOAAA4jGC3S0+c30/F5TV6dEGu7ZyAwxD2Iy8sWq/8ov2aPiFDkaFBtnMAAEAjZCTE6OpTk/X20i1aUlBiOyegNGoIG2NGGmPyjDEFxpi7D/Hjtxpjcowxa4wxnxljujZ9Kn5JzvZ9enFRgcYP6KzhveJt5wAAgCPwxzN6KqltK909M1OVtfW2cwLGYYewMcYt6QVJ50hKlXSxMSb1oIetlDTIcZy+kmZIerypQ/Hz6ho8umPGarVuFaL7zj34Pw0AAPB2YcFuPTaxrzaXVurxj/Js5wSMxjwjPERSgeM4hY7j1Ep6R9LYAx/gOM4ix3Eqf3zzW0kJTZuJX/LKV4XK3r5P08alqU1EiO0cAABwFI5PbqvLh3bVG99s1LKNpbZzAkJjhnBnSVsOeHvrj+/7OVdKWnAsUWi8dbvK9cyn+Rrdt6NGpne0nQMAAI7BnSN7K6FNuO6asUbVdQ22c/xeY4awOcT7DnljbGPMryUNkvTEz/z4NcaY5caY5cXFxY2vxCHVN3h0x4w1igwL0kNj0mznAACAYxQRGqRHJ/RVYUmFnvpkne0cv9eYIbxVUuIBbydI2n7wg4wxZ0iaJGmM4ziHvFeg4zivOI4zyHGcQXFxcUfTiwO8+vUGrd6yVw+MSVPbyFDbOQAAoAmclNJOFw/por8tLtTKzXts5/i1xgzhZZJ6GGO6GWNCJF0kae6BDzDGDJD0sn4YwUVNn4mDrS/erz99sk5npbbXeX25JAIAAH9y76je6hAdpjtmrFFNPZdINJfDDmHHceol3SjpY0m5kt51HCfbGPOQMWbMjw97QlKkpPeMMauMMXN/5tOhCTR4HN01Y43Cg92aNj5dxhzq6hUAAOCrosKC9cjEvioo2q9nP8u3neO3GnXXBcdx5kuaf9D77jvg+2c0cRd+wRtLNmr5pj166sJ+io8Ks50DAACawWk943TBcQl66ctCjUzrqIyEGNtJfoc7y/mYjSUVeuLjPI3oHa/xA37p8A4AAODrJp+bqrYRIbpjxmrV1nts5/gdhrAP8Xgc3TFjtYLcRg+Pz+CSCAAA/FxMeLAenZihtTvL9dznXCLR1BjCPuS1JRu1bOMePXBemjrEcEkEAACBYETv9jr/uAS9+MV6rdm613aOX2EI+4jC4v16/KO1OqNPvCYM5JIIAAACyZRzUxUXGarb3l3NKRJNiCHsAxo8jm5/b7XCgt1cEgEAQAD6v0sk8ov268+fcolEU2EI+4BXvy7U95v36qGxaYqP5pIIAAAC0bBe8bpocKJe/nI9N9poIgxhL1dQVK4nF67T2WntNaZfJ9s5AADAokmj+6hDdJhuf2+1quu4ROJYMYS9WH2DR7e9t0YRIW5NG8clEQAABLqosGA9fn4/rS+u0FOfrLOd4/MYwl7slcWFWr1lr6aOS1dcVKjtHAAA4AVO7tFOlx7fRX9dXKgVm0pt5/g0hrCXyttZrj9/kq/RGR11bl8uiQAAAP/fPaP6qHPrcN3+3hpV1tbbzvFZDGEvVFvv0R//vUrR4UF6aGya7RwAAOBlIkOD9MT5/bShpEKPLVhrO8dnMYS90LOf5Stnxz49MqGv2kZySQQAAPipod3b6sqTu+mNbzbp6/wS2zk+iSHsZb7fvEcvflGgC45L0Jmp7W3nAAAAL3bH2b2UEh+pO2asVllVne0cn8MQ9iJVtQ267d3V6hgTrvvOS7WdAwAAvFxYsFtPXdhPReU1enButu0cn8MQ9iKPLsjVhpIKPXlBP0WFBdvOAQAAPqBvQmvdNCJFM1du04LMHbZzfApD2Esszi/WG99s0u9O6qah3dvazgEAAD7khuEpyugco3tnZaqovNp2js9gCHuBsqo63fH/2rv36KrKO43jzy8XEiEhCCRG5WZIgEAhiVCHylIUvGC1yhQQENFxZqwgXlpbUMfLoNUqikWl1qrtzFgFb5QC3sZ7LS3jDWLCJYCAkauByCUQSELIO38kuEIaISEnec/J/n7WyuJkn32yn3DW2nmy8+73fSVfPZPbadqI3r7jAACACBMbHaVZY7NUWnFI/zF/uZxzviNFBIpwGLhn0Urt2FeuWWOzFR8b7TsOAACIQOkpibp1RB+9W7Bdryzd7DtORKAIe/bm8m2an7tFN5ybrgFdOviOAwAAItg1Z/bQ4LSOuvfVVdq0c7/vOGGPIuzR13vKdPuflyurS5JuGJbuOw4AAIhwUVGmmWOyZJJ+9tLnOlTFEImjoQh7UlXlNHVensoPVmnW2GzFRvNWAACAputyYlvdO7KfPvtql3734XrfccIa7cuTZ/+vUIu/KNadl2QqLTnBdxwAANCKjMw+VZcMOFmz3lmr5Zv3+I4TtijCHqwt2qsH3lyt4X1SdMUZ3XzHAQAArYyZ6f6R/ZWcGKebX8rVgYpDviOFJYpwCyuvPKSbX/xc7eNjNGP0AJmZ70gAAKAVSmobq0fGZGnDjlL96o0C33HCEkW4hf367bUq2FaiGaMGqHNCnO844/p1PAAADepJREFUAACgFTszvbOuPes0PffRV/pg9XbfccIORbgFLVlfrKcXb9CEf+qm4Zkn+Y4DAAAC4BcX9laf1ERNnZen4n3lvuOEFYpwC9mz/6B+/nKeenRqpzsuzvQdBwAABERcTLQeHZetkgOVuu1PrDpXG0W4BTjndMeC5dq+t1yPjs1W2zYxviMBAIAA6ZPaXtNG9Na7BUWa8/FG33HCBkW4BbyydLNey9+mW87vpayurB4HAABa3r8OOU1nZXTWfa+v0hdFe33HCQsU4Wa2Ycc+TV+0UoPTOmrS0J6+4wAAgICKijI9cnmW2rWJ0Y0v5KrsIFOqUYSbUXnlId30Yq7axETp0bE5io5iqjQAAOBPSmK8Zo7J0uqv9+rBN1f7juMdRbgZzXxrjVZsKdFDowYoNSnedxwAAACd2ydF1wzpof9ZUqj3Cop8x/GKItxM/rp2h55Z/KWuHNxNF/RL9R0HAADgW7dd1EeZJ7fX1Hn52l5S5juONxThZlC8r1y3vJynjJQE3XlxX99xAAAAjhAXE63Z47O1v6JSt7ycp6qqYE6pRhEOMeecpr6Sp5Kyg5p9RY7iY6N9RwIAAPgH6SmJuvuSfvrbumI9s3iD7zheUIRD7L//XqgP1uzQHT/MVJ/U9r7jAAAAfKfxZ3TViH6pevitNcrbtNt3nBZHEQ6h/M279cCbBTovM0VX/aC77zgAAABHZWZ6cFR/pSTG6cYXclVSdtB3pBZFEQ6RkrKDumFurpIT4jRzTJbMmCoNAACEvw5t22j2FTnasvuAbp8frCWYKcIh4JzT7fOXa8vuA3p8fI46tG3jOxIAAECDDezeUb+4oLdez9+muZ8EZwlminAIvPDJJr2ev00/v6CXBvXo6DsOAABAo113dprO7pWse15dpYJtJb7jtAiKcBMVbCvRPa+u1FkZnTXpbJZQBgAAkSkqyvTry7PU4YRYTZm7TKXllb4jNTuKcBOUlldqytxlan9CrGaNzVYUSygDAIAI1jkhTo+Ny1FhcanuWrCi1Y8Xpgg3wV0LV+jL4lI9Ni5bnRPifMcBAABosh/07KSbhmdofu4WzVu62XecZkURPk7zlm7W/GVbdNOwDJ3Zs7PvOAAAACFz47AMDU7rqLsXrtQXRXt9x2k2FOHjsPrrEt25YLkGp3XUTcMzfMcBAAAIqego02PjctS2TbQmz2m944Upwo20r7xS1z+/TInxsXp8fI6iGRcMAABaoZPax+vx8Tlav2Of7vhz65xfmCLcCM453fqnfBV+U6rZ43OUkhjvOxIAAECzGZLeWbec10sLPt+qOR+3vvmFKcKN8OySQr2ev01TL+yjwWmdfMcBAABodlPOTdc5vZN176urlL95t+84IUURbqBlG3fp/jcKdF7mSbru7DTfcQAAAFpEVJRp1uXZSk6M0+Tnl2n3/grfkUKGItwAO0srNGXOMqUmxeuRMVnMFwwAAALlxHZt9MSE07V9b5lueTlPVVWtY7wwRfgYDlU53fxirr4prdCTEwYqqW2s70gAAAAtLrtrB911SV+9v3q7nvxwve84IUERPobfvL9Oi78o1vQf9dP3Tk3yHQcAAMCbiYO760dZp+iRt9doyfpi33GajCJ8FB+s2a5H31urH59+qsaf0dV3HAAAAK/MTA/+uL/SkhN049xcbd19wHekJqEIf4fC4lLd/EKuMlPb6/6R/WXGuGAAAIB2cTH63ZUDVV5ZpcnPL1XZwUO+Ix03inA99ldUatLzSxUVZXpq4kCd0CbadyQAAICwkZ6SoEcuz1Le5j36z4UrI3axDYpwHc45TZuXr7VFezV7fI66dmzrOxIAAEDYubBfqm4clq6XPtukuZ9E5mIbFOE6fr/4S71Ws2jGWRnJvuMAAACErZ+e10vn9E7W9EUrtfSrXb7jNBpFuJYl64r1wJsF+mH/VE0ayqIZAAAARxMdZXpsbI5O6XCCrp+zVNv3lvmO1CgU4Rpbdh/QDS/kqmdygh4ancXNcQAAAA2Q1DZWT00cqJIDlZoyZ5kqKqt8R2owirCksoOHNOm5pTpYWaWnJg5UQlyM70gAAAARo09qe80YPUCfFu7Sfa+v8h2nwQLf+A7fHLdi6x79/qpBSktO8B0JAAAg4lyadYpWbNmjp/+6QZknt9f4M7r5jnRMgb8i/Nu/rNeivK2aemFvDc88yXccAACAiHXriD4a2itZdy1YoY83fOM7zjEFugi/s6pIM99eo8uyT9HkoT19xwEAAIho0VGmx8fnqFuntpo8Z5k27dzvO9JRBbYIr/l6r376Yq4GnJqkGaMGcHMcAABACCSdEKs/XP19VR6q0rV//Eyl5ZW+I32nBhVhMxthZmvMbJ2Z3VbP83Fm9lLN8x+bWY9QBw2lnaUV+vc/fqp2cTF6auIgxceychwAAEConNa5nZ6YcLrWFu3Vz176XFVV4bny3DGLsJlFS3pC0kWS+koab2Z96+z2b5J2OefSJc2SNCPUQUPl4KEqXT9nqYpKyvX0VYOUmhTvOxIAAECrc1ZGsu68uK/eXlWkWe+u9R2nXg25InyGpHXOuQ3OuQpJL0q6rM4+l0l6tubxPEnDLUzHGkxftFIfbdiph0YNUHbXDr7jAAAAtFrXDOmhsYO6avb76/Rq3lbfcf5BQ4rwqZI21fp8c822evdxzlVK2iOpUygChtJ7BUWa8/FGTRraUyNz6n4LAAAACCUz0y9Hfk/f73Gips3L14695b4jHaEh8wjXd2W37kCPhuwjM/uJpJ9IUrduLT+33Lm9UzRzTJb+mRIMAADQItrEROnJKwcqd+NuJSfG+Y5zhIZcEd4sqWutz7tIqntt+9t9zCxGUpKknXW/kHPuaefcIOfcoOTk5ONL3ARRUabRA7soOiosR20AAAC0Sp0T4nR+3/Bbr6EhRfhTSRlmdpqZtZE0TtKiOvssknR1zePRkt53zoXn7YEAAACAGjA0wjlXaWY3SHpLUrSk/3LOrTSzeyV95pxbJOkPkp4zs3WqvhI8rjlDAwAAAE3VkDHCcs69IemNOtvurvW4TNKY0EYDAAAAmk9gV5YDAABAsFGEAQAAEEgUYQAAAAQSRRgAAACBRBEGAABAIFGEAQAAEEgUYQAAAAQSRRgAAACBRBEGAABAIFGEAQAAEEgUYQAAAAQSRRgAAACBRBEGAABAIFGEAQAAEEgUYQAAAAQSRRgAAACBRBEGAABAIFGEAQAAEEgUYQAAAAQSRRgAAACBRBEGAABAIFGEAQAAEEgUYQAAAAQSRRgAAACBRBEGAABAIFGEAQAAEEgUYQAAAAQSRRgAAACBRBEGAABAIFGEAQAAEEgUYQAAAASSOef8HNhsh6SvvBxc6iyp2NOxAaApOH8BiFQ+z1/dnXPJdTd6K8I+mdlnzrlBvnMAQGNx/gIQqcLx/MXQCAAAAAQSRRgAAACBFNQi/LTvAABwnDh/AYhUYXf+CuQYYQAAACCoV4QBAAAQcBRhAAhzZtbVzOaZ2R4zKzGz+WbWzXcuAIh0FOEaZjbdzP7Fdw4AqM3M2kp6X1IfSVdLmigpQ9IHZtbOZzYAaKhw7VmBLsJmNsTMLq+zLdrMJplZb1+5AKCWayWlSRrpnFvgnFso6VJJ3SVd5zUZABxFJPSsQBdhSRslnW9m76j6astgSYsl9ZC03WMuADjsUkkfOefWHd7gnPtS0t8lXeYtFQAcW9j3rFZVhM1srpm5o3y8VXt/59wm59y1kh6WNFLSOElTnHO3Oed2fccxzjez18xsi5mVmdkmM5tpZrHN/x0CCKB+klbUs32lpL4tnAVAwDWma0VCz4ppji/q0a8k9Ze0S9K0mm1Jkv5X0oOSnqq9s5mdIukuSemSFkjaLekJM/tQ0kPf8SZlSXpP0m8llUo6XdIvJe2sOT4AhFJHVZ/T6top6cQWzgIADe5akdCzWlURds6tMLMukl5zzn0kSWY2tObp15xzhXVekibpL865yWY2XVKhpCmqHneXonp++DjnZh5+bGbRqv7z5DmShoTyewGAWuqb8N1aPAWAwGtk1wr7nhW2RdjMzpP0TgN2/dA5d07Na7pL6iApr9bz2ar+IbK87gudc3+rZ9shVf8WUl+mGElXqvoNzJDUqdbTLzcgKwA01i5VXxWu60TVf6UYAJpNY7pWJPSssC3CkpZIymzAfvtrPc6q+Te/1rYcSYXOuZKjfRHn3PQGHGuupIsk/UbSfZKKJcVL+kD1j+EDgKZaqepxwnX1lbSqhbMAwHF1rXDtWWFbhJ1z+yWtbuTLBkgqk7Sm1rYcHflby3Exs2xJYyRNcM7NrbV9tKr/RJnb1GMAQD0WSZppZmnOuQ2SZGY9VP1nwts85gIQTM3StXz1rFY1a4Sqf0tZWXPZXTV3GGYqBEVY0uFVnL5942sms7+v5tNlITgGANT1jKrH1S00s8vM7FJJCyVtUp0bgAGgBTRX1/LSs1pjEa79RqRKilX1nYZNlSupQtLDNVN7TJT0kaQ2krY757aG4BgAcATnXKmkYZLWSnpO0hxJX0oa5pzb5zMbgEBqrq7lpWe1miJcswxpTx05ZmWHqseUzKj5Dz1uzrlNkiaoejWnRaq+63GapHViWASAZuSc2+icG+Wca++cS3TOjaxnFhwAaFbN2bV89Sxzrr5ZeQAAAIDWrdVcEQYAAAAagyIMAACAQKIIAwAAIJAowgAAAAgkijAAAAACiSIMAACAQKIIAwAAIJAowgAAAAgkijAAAAAC6f8BdIK246R5/RUAAAAASUVORK5CYII=\n", "text/plain": [ "