gimx analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab qt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{\n",
      "  \"stdin_port\": 40042, \n",
      "  \"ip\": \"127.0.0.1\", \n",
      "  \"control_port\": 46206, \n",
      "  \"hb_port\": 60907, \n",
      "  \"signature_scheme\": \"hmac-sha256\", \n",
      "  \"key\": \"171997fc-394e-4788-838b-77e93811bdc2\", \n",
      "  \"shell_port\": 33663, \n",
      "  \"transport\": \"tcp\", \n",
      "  \"iopub_port\": 59033\n",
      "}\n",
      "\n",
      "Paste the above JSON into a file, and connect with:\n",
      "    $> ipython <app> --existing <file>\n",
      "or, if you are local, you can connect with just:\n",
      "    $> ipython <app> --existing /run/user/1000/jupyter/kernel-8d1525d7-a45a-41c1-b02a-9a99b1699861.json \n",
      "or even just:\n",
      "    $> ipython <app> --existing \n",
      "if this is the most recent IPython session you have started.\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "%connect_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try to mask out HUD/UI elements by detection regions that slowly change from frame to frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'images' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-3-9f6059782af9>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mimages\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimages\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mdiffs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mimages\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mimages\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mmask\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mstd\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdiffs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mthresh\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.3\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mmask\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcv2\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mGaussianBlur\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmask\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'images' is not defined"
     ]
    }
   ],
   "source": [
    "images = array(images)\n",
    "diffs = images[1:] - images[:-1]\n",
    "mask = std(diffs, 0)\n",
    "thresh = 0.3;\n",
    "mask = cv2.GaussianBlur(mask, (0,0), 2)\n",
    "mask[mask < mask.max() * thresh] = 0\n",
    "mask[mask >= mask.max() * thresh] = 1\n",
    "\n",
    "imshow(mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import cv2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fdc2bce3650>"
      ]
     },
     "execution_count": 190,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAADdCAYAAACmGfmGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8E2X+wPFPW25KuSk3RQTkEDkULYgUBcQVxIMVFY8F\n3J+r7KroKuLx63Y9wVXU37quuwK6iCgLioCKHFJEAZH7PqXcLYdIuaXH74+ZDEkzSSbJTGYm+b5f\nr7yaeTKZ+U6TfPPkmWeeB4QQQgghhBBCCCGEEEIIIYQQQgghhBBCCCEcoR+wBdgOjLI5FiGEEGFK\nAXYAGUB5YA3Qxs6AhBAiniVbsM2uKIk8DzgPfAwMtGA/QgghsCaRNwL2ei3vU8uEEEJYoJwF2ywN\nvUp6KRRYsGshhIhnzYDdSWVLrUjk+4EmXstNUGrlXgqAnl7LGerNrXKBLJtjMFsu8XVMucjxOF0u\n8XVMuUR/PHnqzWOR7lpWJPIVQEuUzHwAGAzc6b9algW7FkKIeJKBbyU3dom8CPgj8DVKD5bxwGYL\n9iOEEAJrEjnAV+otQWTYHYAFMuwOwGQZdgdgsgy7A7BAht0BmCwjZnuyotdKAsqwOwALZNgdgMky\n7A7AZBl2B2CBDLsDMFlGzPYkiVwIIVxOErkQQricJHIhhHA5SeRCCOFyksiFEMLlJJELIYTLSSIX\nQgiXk0QuhBAuZ9WVncKA2hxlCJMBmMAwTpJqc0QilkbwNikUU5NjWtkxarKeS1lIr5jHczXf8R1X\nx3y/Inp+wyHGSClkawuVOMsdfBzWBhpygPKc15YLSWMFl7OYHrrrX8cCmqjDpL/P73TXuYXP6MC6\nsOIw226a+cVXl8PcyBeUkMx/uNfvOe3ZQE8WUYcjUe33a67nIA0i3oYI7ne8D0AzdtsbSAinqMp6\nLuVrrvd7zHMMgT5Dd/ERFfhVW/6RK9hIOyvCTFA5oJO3bUvk2aHXichJUllHB+bRB4Ab+IpLWU9l\nzvisl4NvBAOYRWdWWRRV+DbQHoA0CmnKHp/Hysb+AO9Sn3xL9u9xjorMpr+p+4h3PVhMPQ5Rn/yo\nvmTtsoRuPp+jriz3eXwmN7GaTtryIKbRjo1h7aPs+6w9GygkjT00ZTq3Gd5OBX5lALP4kt9whsph\nxeAuCZLIjdrCJXzCYG15FGOoxFkbIzLubUZwhDracrby4sbESrqwmB4cp3rM9ukmyZRwHQvoxhK7\nQ7FcPvV5lwcAqMwZnmSsqds/R0Ve4SlD6/6GL7mCHwH/ik580U/kCXuy8xK22B1CxEbwtm377sJK\nHuUNbmKmbTE4USon6cVCnuP5hEjigM+vwIqcM337FTnHs7xgaF1PEk9UcrLT5foy15b9dmI1pSSx\njKs4TF1bYnCCWvzM5awgk6V2h2ILq38NbqKtpduPF5LIgSRKSabE7jBcpzOr6MwqtnAJC7jOp7kn\n3tXiZ3qwmI6ssTuUuPYpt1q+j9oc5Si1Ld+PlSSRA6UkUZK4rUxRu4QtXMIW8qnPFO6kkDS7Q7JM\nNU5wMzO4iJ/sDkWo6nHIZ7kR+9lPI8PP78M8dtOMpWSaHVrMRJu98oB1wGrQTmnXAuYB24C5QI0o\n9yFcoj75jGQcKRTbHYolUijmMV6XJB5DkbyXwvl1nUwJrdnq+l/k0SbyUpRZlDsBXdWyp1ASeStg\ngbosEsizvOCaHkBGVOYMT/Cq4RNvwjyjeTnkOp6+7R7DmMDF7DC0/bv5EIDezA87Nicxo2mlbFeY\nm4Ce6v0PgFwkmSecUYzhCHUoUt9inm5qbtKI/dzMjJj2Ac8li0XqxyeW3Uoj5Xldq3GCu/jI9O2/\nwLMh1yl7jQgo3RHf4uGQz23OLu1+H+Zp/ebdJtpEXgrMB4qBd4F/A+lAgfp4gbosEpB3Aswmh2PU\n5Ah1+Ii7bIwquE6spg2bacl2u0NhPMMZzni7wwgqn/oAnKWSLfsPdCFcTY7RhZWspEvQ5xeQTrqW\nrtwr2kTeHTgI1EVpTinbObtUvQlBTY5Rk2NR1TS9P5iN2RfwQ+i9XjN2U4cjAT/UXVgZcTyJbAWX\n2x0CVTgd8LGm7AmZyL/hWu5kitlhxVy0ifyg+vcw8BlKO3kBUB/IBxpAmVPKqlyv+xnE3/zZVuvJ\nIkBJZonEaNLVW08Stvv0YHHA8ZMA7mFSwMf20DTk9rfRKqK4YidPvQUXTSKvAqQAJ4CqQF+U60dn\nAvcBY9S/M/SenBXFjgVk+XwVilhbRWftvlVj9LRmqyXbdZNr+YZfqcBaLvNrvgk1wF18XNeQgW81\nd5HuWtEk8nSUWrhnO5NRuhuuAKYCw1G+Sm6PYh8igB1cDCjt0DX4xbY4vE/OeXiG5gXlhKHeySg3\n2sHFzKGf38UjViXyq/nOku26TT/m0I85AExmCABt2UQnVgd93m6aRbXfJ3iVV3kiqm3ESjSJfBfQ\nUaf8Z6B3FNsVBnje0H2Z67jLwz2xefRmPp1YHbQ904l+pQI7acFUqYvosuML2ruSYAbvseC7scSn\n14qb+pYn9JWdFTnHOSraHYareY8JH8h8ejOf3nRlOVnkOr6GfobKbKQdX3Bj0PXKcz6hB2tqx0am\nMcjuMAIqR5HW/TWQAczSLU+mhEqcJZkSV1z1ndCJ/G4+ZDzD7Q7D1brzPVfyAy/yTMh1l9OV5XSl\nHEU8w4sxiC48M7iZjbQL+eH3eIpXTK+1BTuxJ8JzGWtD9loJxG3jxyd0IhfmKEeR1k/cyEUYRZQj\nh2ySKNXa9x/mLavD9FFIms8sN8eoafi5D/OWz09yMzmt9uc9trcbLlDymE1/Q0n8M27hMV4HlIk0\n3CqhE7nUxs1Vk2M8w4vspYnulHRllZKkJVDvhBHJVGjel2n/h3sDJsTTVIlq2N3RvOwzlVki20pr\nu0MIyGhN/ATVLI4kNhI6kQvzlaOI5uzidqZGfJIwkt4GVs8KE6v5XHuyyK8XkFM5NQm+x/12hxBz\nksiFJdqwmWxytJ+rbh3Dog/zEmbGn0Tn/IuDAkvoRD6c8dK8YjFPEuzGEubS1zVjPmey1LbZl0R0\n2rMhrPHIPaLtd26nhE7kIrb6MpdL2MJCepHn0EEZ6pNPJktj0owirHEVyyghmRVcrp2DSaOQOhzh\nJy7yWTfDwOXvbiCJXMRUU/ZwHx9wjJpM57aIak5ma8BBynOeoUy0OxRhkm4s0W0S+ze/5wANteV4\n+dUliVzYoibHuJ/3tOU3eUR3vV8MTjBVjiJSOaktF5Lm13OlIudozwb6MzuCiEVZoS6YcqLf82/t\nvXYPk6jFz7rrPcg7sQwrapLIhSM8wptBHw80wUAmS7mOBbqPTWCYVvtyy+w+yZTG7VR5ThHqvZbk\nwpG3JZELV4gkEQ9jggWRWOsaFnFNgBHuRGw8yVi7QwhbQifydAp4gHcB5CIPIRKcm65cLSuhE3l5\nzgecKsrp3PymE+4h7zN3cNbADkIIVwp0nkLEhiRyIUTU2rPB7hASWkI3rQghIueZGcmqkSCFcZLI\nhRARCTQpQzy5mB2uGIPFSNPKBKAAWO9VVguYB2xDmafT+6qN0cB2YAvKhMxCCOFKdzLF7hAMMZLI\nJwL9ypQ9hZLIWwEL1GWAtsBg9W8/4B8G9yGEECJCRpLsYvBrBLsJ+EC9/wFws3p/IDAFOA/kATuA\nrlFHKYQQIqBIa8vpKM0tqH/T1fsNgX1e6+0DB4yKJIQQccyMk52l6i3Y435yve5nqDchhBDe8tRb\ncJEm8gKgPpAPNAAOqeX7gSZe6zVWy/xkRbhjIYSIled5zuYIMvCt5uqPwxNp08pM4D71/n3ADK/y\nO4AKQHOgJbA8wn0IIYQwwEiNfArQE6gD7AX+F3gFmAoMR6n3e2bZ3aSWbwKKgIcI3uwihBAiSkYS\n+Z0BynsHKH9JvQkhhIgB6eMthBAuJ4lcCCFcThK5EEK4nCRyIYRwOUnkQgjhcpLIhRDC5SSRCyGE\ny8nEEkKIuLKInuSqg4BkkEcGefQMcGl7vJAauRAibozlSS2JA+SRQS5ZvMsDEW3vOZ43KTJrSY1c\nCBEXnuc5SgLUTfOpH+NoYktq5EII1ztNlYBJ3GM8wyPadhqFET0vliSRCyHiUjY59Ga+thwq0btZ\n/B6ZECJh3c5UALrzvVZ2gIZ2hWM5SeRCiLiTTIluuf0TRVhDErkQIu58zB0AvMGjPuUlJPNX/pdf\nqWB4W4WkmRqbFSSRCyHiUg7ZHKe6X3kpSbzMaBsiso4kciGEcDkjiXwCymTL673K/gLsA1artxu8\nHhsNbAe2AH1NiVIIIUx0BT8aWu8HrrQ4EnMYSeQTgX5lykqB14FO6u0rtbwtMFj92w/4h8F9CCFE\nxKpwOuAJTj3X8o2h9ea6pC5qJMkuBo7plCfplA1Emaz5PMqkzDuArpEGJ4Qw3ybakkM2OWSznks5\nSyW7QzLFczxPCsWG1t1HY4ujia1oast/AtYC44EaallDlCYXj31Aoyj2IYQw2X/5rXb/U25lDKPC\n6sXhZM/ygqH1JjPE0Ho1deuwzhNpIn8HaA50BA4CrwVZt1SvMNfrlhdhEEKI8OnVWl9mNKW6P7Ld\npxH7TdvWH/m7aduKTB6+2VJfpIn8EEqCLgXe40LzyX6gidd6jdUyP1let4wIgxBChC9QclrGVTGO\nRISWgW+21BdpIm/gdf8WLvRomQncAVRAqbG3BJZHuA8hhJdPGKy1bb/OYxFvpwa/cD/v+ZW75cRe\nKLfwWdDHUyh2zfC0RhkZxnYK0BOoA+wFslG+Gjqi1Mh3gTbY7yZgqvq3CHiIAE0rQgjjpjGILVyi\nLZ+gGvPoQx/mRbS9hhygAr/6tY2vphOdWB1VrHarzVGyyfEpO0kqB9X6Z0u22xGWpYwk8jt1yiYE\nWf8l9SaEMMlOWviVtWVTxNtLolS3u94R6kS8TSdL5WRcJnAP6eMthAuMYox2vw2bGcaEqE/qeW9T\nuJvMECSESwxgFgDN2WVZt7gC0i3ZrpvpNUE5jSRyIVxgDKN8LtzpzCotsZtJrwkn0VXirOMTuTSt\nCOFCq+jMPPpEtY3D1DUpmvjmhpmFnB+hEMISEximW+702mesnSTV7hBCkkQuRAJ6kWcCjrHyBTfG\nOBoRLWkjF8KlVtE57H7k+dTnY+6gKMhHfx0dAP0La2bTn5+4CICHeSusfQvrSCIXwgXu4iO/ppCz\nVGIxPejB4pDP30prbfozI9bRgSRKuZkZPuUnSeUYNQ1vR8SGNK0I4QJN2Ktb/g3XspgeIZ8fThL3\nWMtlTOc2V5zsC+UADZnObczxm1ohtDQKLYjIXHYNd1aaHeMdrqGjdpa+DZsZz3BAGWDeSI1GCDst\npBffco2hdUcyzi/55OD7iUumhOtYYLjnyyjGUImzxoJ1mEncozUHeXRlOcvVsf4e4F3qkx90G2X/\nf/bJAZ28nTBNK58zULu/hG7a/VV0dkUi934jlR1HQsS/5VHOzzKMCWymDQB9mauVd2MJoDSl5FOf\npWTqPn8Mo3ic10jlZFRxxFouWX5JHHz/n4epGzKRV+e47kTOTuH+30wGVeOE3SEIERNbae1X1oS9\n9GWuTxL31oF19GUu2eTQjo1k6MwS8BqPs5F2ZofrCk5O4pBANfIOrON7uvuUVeQcj/CmTRGF53Fe\n85nZRYhAvuQ3pFBMZ1b5PbaHpkxkqF95TxZxET/RlD0MYhoA/+J/tBEDPaYxiAM0jHjUxVirzvGo\ntzGOkSZEYq2EaSMXws3KXqJvlHdF5U0eCes5wda/imVcz9dhx2OXv/FnTlFV97FQTZXjGEkhaVaE\nFQH9NnJJ5EK4hHNOuCm6sJL+zLY7DMP208jv18jdfKjbjOTNDYk8YZpWhHC72hzlKLXtDkOzki4U\nUc6vr7lTNWK/4cmZ3SZhTnYK4XZJDpxsS+/Eqoi9UIm8CbAQ2AhsAB5Wy2sB84BtwFyghtdzRgPb\ngS0QJ5MACiF0RdJuL8wXqmnlPDASWAOkAitREvhQ9e9YYBTwlHprCwxW/zYC5gOtQGdOKSESwGz6\na/fd1J4s3CVUIs9XbwAngc0oCfomlAmZAT4AclES+UCUyZrPA3nADqArsMzEmIVwvEX0JJcsn7KV\ndPFbL4tcn+UM8mjGbt1txut8miJ64ZzszAA6AT8A6UCBWl6gLgM0xDdp70NJ/ELErV+ooV09WEJy\nWMPAlk32ZVXhNL1Y6OihZccznOGMtzsMy4xknON6DJVlNJGnAtOBR8DvEslS9RaI7mO5Xvcz1JsQ\nbvIztfiBK6O+fD6Y01RxdBIH2Edju0OwXCXO2nQ+IE+9BWckkZdHSeKTQOtnVADUR2l2aQAcUsv3\no5wg9WislvnJMrBjIZyqkDT+jz/ZHYaIkTQKbUrkGfhWcxfprhWq10oSMB7YBLzhVT4TuE+9fx8X\nEvxM4A6gAtAcaAksNxyzEC7xDdfaHYKIkbNU4hD17A4jqFA18u7A3cA6YLVaNhp4BZgKDEep99+u\nPrZJLd8EFAEPEbzZRQjXmcrt2kiCIv6NYZTdIYQUKpF/R+Bae+8A5S+pNyHi0g4utjsER2nMPrtD\nSHhyib4QYZhNf85T3u4wHKMpexjKRLvDSHhyib4QYcinvt0hOMp1LLA7BMu5Yao3SeRCiIgFmlEo\nnoxknN0hhCSJXAgRsS1cYncIAknkQogoeSYyF/aRRC5EGPbLiBN+zlHR7hASniRyIQxy+gS8djlM\nXbtDSHiSyIUwaALD7A7BsYzOByqsIYlcCBE1+bViL0nkQhgUaBZ2AaUk8QLP2h2GZRrpj/3nGJLI\nhTDgY+6gmBS7w3C0YlJk6jebSCIXwoDttLQ7BFdwwwBT8UgSuRDCVPE4JZ3Tu51KIhfCRu3ZQBVO\n2x2Gqd5mBCeoZncYCUUSuRAGPMfzlmz3NqbzBK+SRa4rBmcyai2X2R1CQpFELoRBqZw0fZuFpAHQ\nk0WMZBzZ5JCuzWvuXgu4Tk58xpAkciEMOENlTpJq+nbHMdKv7A/8k2xyqMwZ0/cXS3KCOHZCJfIm\nwEJgI7ABeFgt/wuwD2X6t9XADV7PGQ1sB7YAfU2MVQjb/I0/W7btAzTULX+SsfyRv1u2XxE/QiXy\n88BIoB1wFTACaIMyD+frQCf19pW6fltgsPq3H/APA/sQIqF9wuCAj9XmKCN4O4bRmOdTbrU7hIQR\nKsnmA2vU+yeBzaD1w0nSWX8gMAXlCyAP2AF0jTpKIeJYIWn8g4cCPl6bo1Tg1xhGZJ54mXginq7s\nzECpfS9Tl/8ErAXGAzXUsobgMxPrPnB4B0whDLiVTy3d/mHq8h736z6WRCk3M8PS/VtlLn1Zz6V2\nhxH3jE6+nApMAx5BqZm/A/xVfex54DUIOLp8aTQBCnNsoi3LuIq9NNHKurGEPsyzMSr3aMdGZtPf\n0p4Y+2nEBIaRyVLasNnnsTZsZiCfs5ge/Ewty2KwwqfcShKltGeD3aHELSOJvDwwHfgQtGrBIa/H\n3wNmqff3g1emgMZqmZ9cr/sZ6k1Y57/81q9sCd1YQjdt+XJWcCNfxDIsV8lkKQvpZek+9tKEvTSh\nHofowWKf5NeRNXRkDevowGfcYmkcZpvObSRTQls2Rb0tveaaTJZGvV1nylNvwYVK5EkoTSebgDe8\nyhsAB9X7twDr1fszgY9QToQ2AloCy/U2nBUyNBFrK7icpuzhUu3lFN6u4Vu20Soml2sfoh7TuY0k\nSmnHRp/HOrCOb7mGo9S2PA4zLaInbdnERtoB0IptlOd8WNuYwDCfX5Ue33E1T/CqKXE6Swa+1dxF\numuFSuTdgbuBdSjdDAGeBu4EOqI0m+wCHlAf2wRMVf8WAQ8Rx00rh6nLGSqHXK8pe2IQTXCpnDTU\nD9pJyWEPTZnIUFI5ye1MpQl77Q6J+3mPHLJjtr9pDGIag7TloUykKXvoynK+8un163yHqOf3v8sm\nx/Dz/8X/cJAGYe3zONV5g0cBaMw+hjM+rOe7RahE/h36J0S/0inzeEm9xZ1ftHO6yrCmBaQbet6t\nfGp7LfdxXvNZfp3HdMfDWERPMllKRc7FKjRd33KN1oxxklQmMIxRjKESZwM+x/v1qcEvlsQV6IRk\nrExkKP2YQ2u2ui6RhzKHfvzAlbRgJwP5nGqc0B4rJC1oEi/V7UTn2wyzj8bmBeswRk92JrxD1OMd\nHozouZ7+tHYnc1De8CUkBx3UaBWdA7Y56o3JnUJxVDEVk+IzKcGzvKDbFj2GUT41uLLP8zaScaaP\nXVJMiiNGwZtDP+bQz+4wTDWb/qykCwA7acFUbvepPU/inqDPP0PlkL+UkqJoHHDC6x6MJPIE8ipP\nUEQ5fqVCWM87Qh0mMhRQPjBlaz/h/DzW8yLP+CwHm2nmKLWpzVGfD36slI1TRO8gDWignW6zltP7\ngkfDrilP/lLCtTRnl027D19FznERu1hLx4iev5k2ZAU4UWG1OfRjMkM4T3lDs9zspIUW61s8zLdc\nw3nKc57y6F0HdopUWrE9otie5zlKw7icYTldWUQWBwNc1u7RmwUkm3B65mdq8T6/40tuRP8aOBGN\nPDLoynLqk88yr2aQe/kPVbzGmunKj+ygJSfUQcbC1ZxdDOX9iJ47jpGco2JEzzXfIsC/5mTXO7MU\nsqnCadqx0adt06Mry7mYHRHv4CPuAmAXzSkq88OjpVfSuYuPwt623jZB6Tq2mB4Bn9eXuSG7SXni\nLiuSOD2sPjn3IO9Qz6dHqjGTuIefuCjq/Xv+Ny0j/DIpy/ManKaK439Sx4P7eS+s2rJnMK5dNDd0\n5Wik70+PcYzURqm0Xw7o5G1bE3koDTnI7/mX4Y3uozHfck1Eo65dxlp6Mz/ioUqN/tTvwko6soY1\nZWr2oZ5bg194hDfDjus97jc1GfVnts9yF1aGvQ2zmkXasonf8l+fsl+pwFyvsdrKxhvMKjqznK6G\nT2IL8zRlD3U5rC2H87oBFFFOdwz0SN6fZUkiD8xQIgeltmWkpvUhd7OTFtHGFVZ77wbas5RuHAiz\nS1QkIk3kZtXGezOf7nyvLS+iJ7leVwNkkUtPnaajt3iEYzq/uMzShL20YCegzHL/I1f4rdOabWSy\nhGbs9nvse7ozn96WxScik0YhIxlndxiAJPJgDCdygHuYxEX8pPtYHhl8wH1mxUUmS+nL3JDrmV3T\nDaYqp/gzf4vouZEk8noc4hD1dB8bwCxmMSDgc6twmj7Mo6M61los+1yHMoBZdGaVT9kYRskECA4V\n7Ul0MznnfayfyF3Ra2US9zCc8TT2GY9L6RJkZhIHpd9pJkt9+rDaLVgS38HF2v186rOA6wC4kym0\nYhuVOKslqovZ4bN+IIGSOBA0iYPSrvw5A/mcgWH/PLbaLAYwiwEMYXJU519E7HkuDitrCJO1+y3Y\nGVUXQzdzRY0c9K/KsvJb8gleDTopbixr5OkUMIhp1OGI32PB/gdlazSnSOVvPG56fEakU8Apqloy\ny04kqnGCx3gdULoV6p28Fvbyfv+a+VlPoZineYlkSgytf4Q6vM0I0/YfHf0auWsnffg7f7R0+yu4\n3NLth6OAdN5mBC/wrKEhAfScobJtSRyUk8mP85pjxtX2viCqPvk2RiICySGbM1RmDKNM3W4xKUGv\npThDZZ/bf7jX1P1bwbXVEKvHBFlIL67hW9O29wwvavcjvbCkmBTG8iTlKGIIk8kgj2d4MeD2Ckkj\njUJH1Djn05tMlnIjXzhu5L7hjOcD7iNPxuB0nLE8GZP9zKc3P3AlgO2flUjY1rTSl0yfbmKheP/M\nitUJqmDjdReSpjtxrkcKxYzgbWpyzO8xM4chrcJpMllKCsVh/T/t4HkNyx7/UCb6nJNYxlUsNzix\n1MO8RVVOcYqqfo+9pU0x6/+cZEqoznGf8lNUZQHXsZpOhvYt3CmFYp7hRZ/2dOeczAzFYSc7w0k6\n3qMHhvpZZKYNtA+YyD1tZrU5ShHlOE517bFUTvoNUuUt2EVD4TpNFRZwHVU5Zdo2y6rJMW7lU0pI\n9rk8fwnd2EYrQ9vw7krWgXV0YF3AdW/gK27wGpct0IfsCn7Uvij1mmyyyWE3zXzK9LogelTlFDcx\nkxSKHdW0Jvz9jve119L7Nf4vv9X9UvcI9dl0K1f8hih70rEkRk37gfqOjmOk9mVStomnFj/zJ/5P\n93nr6MBSMjlCHXMDhaBvXm+3MZ3p3OZTNoBZfomwDkdCth03YzeHqMcyrtKtxfZlLh1YZ9mXjN6v\nnbKCJe5otivs054NPp9N79e4ObvYSQvOUJnOrPIbBiTQLEUt2R7RhYROYVvTykPUDTrhrLdL2MJg\nPgGCj3hnhQ6s4xY+4ySpLCWTpWQGHDITAvd9nUcfn9l4rFSe82SpczB1YF1YV6vOow+gfDhasc2K\n8MIWqEZuVT/jAzTk3/zekm0Lc1XnOO3YaOiz1Y0lVOIsPVis+7g7XneHNa3U5TDZ5Di+beoIdVga\nRnv+GzxKTxbRSZuHQxFOEr+In0ingDZs5nu6s5XWAddtxm4ackBbTuUk3VhieF9lecfplESeTQ4l\nJPtcgWnkoq1Ief8/hbMdp7rhz5ZnvVyy6M18vqe79ku2NVtpxm4uZb0rJ4u2LZF/yW90L6d2mgM0\n5ECIkfa8Hac633KNXyKvznGtHb0Gv2iDBA1iWtDt3cHH/EwtvuFav8eu5AdHzJoTC8mUWJq8ReIo\nIdmvYraV1kErTE4XKpFXQhk3sSJQAfgcGA3UAj4BmqHMDHo7aFOyjAaGAcXAw6D/6TOaxFM5qTWr\nuNmjPlOehqcWP4dM+GZx0mXRTtaNJbRmK/tpxCo6W3LeQwijQiXys0Av4LS67nfA1cBNwDxgLDAK\neEq9tQUGq38bAfOBVmDwEiqUK+5uZkbAsVWcqOy0YpEMbiWcw8jY00voFrNzHkKEYqRpxdNlpALK\nRBTHUBJ5T7X8AyAXJZEPBKYA51Fq6juArsAyI8FYMT2X2VIo5lleiOi5nvMBUut1tljPPCREtIz0\n40sG1gAFwEJgI5CuLqP+9Qzg3BB8RrbaB8YGJHmM13WT+Bo68ipPGNlETDwdn/NKCyFczEiNvATo\nCFQHvgZyuI2wAAAJSUlEQVS/WXFL1VsgIYcje5B3dEcbXMHlfMGNgFKbbcl2BvMJyZTErC+5twd5\nx/BAO0KI6KVQzP2851Pmub4hn/p+5a/xuGMGZoulcHqtHAe+ALqg1MLrA/lAA9DmUdoPNPF6TmO1\nTEcuAFU4xWmdaZiW05WvuMGnbDstY9qH3Fs9DkU1XRRIk4oQ4Somhfn05m4+9HtM74K1x3mN7bQM\nOGWi++Spt+BCVWvrgDa9S2WgD7AamAnaQOD3ATPU+zOBO1Da05sDLYHl+pvOArIo5Vq/oYre4FG/\nJG63B3nH7hCESEg7acE0fmu4Z1BLtnM3H9Kf2Y65FiJyGXhyJV4zcpUVqkbeAOVkZrJ6mwQsQEnm\nU4HhXOh+CLBJLd8EFAEPEaJp5QyVeYmnGc3L2iA23uOWOIGbetAIEY820paNtDU8kbJn+j/PnJ3T\nGMRG2lkao51CJfL1QGed8p8h4ESHL6k3w85TnlKStETuffGME8hsMkI4wzs8qDtlXyiDmBbXidwx\nE0t4T0H2KG84pj25GifIZKndYQghVLMYQA7ZhqYt9FZ2qsh44phEPoU72UkLbXk2/W2MRlGF09p0\nYEIIZ5nMEHLIpoB0DlM35PpOm0PWTI4axvZD7gaI+IIbM/2GL7mCH+0OQ9igHEV2hyDC8E/+oN3/\nA/8kjUIqc8ZvvXTt0pf445gaubeXGW3r/m9mhiTxBNY1UEcr4Xj/5A+M5Und6Q9jNVm6HRxVI/co\nJoX+zKY73/uUF5LG+/zO8v1fxlrL9yGEsE4R5cghm5oc4xwVOU0Vu0OylCMTOcA2Wvn1Aa3JMbLJ\nYS9NmM5tpvdsuYXP4qDfqRDC4xg17Q4hJhybyE9QLeBjTdjLo7zBLpqzkF7s9bmYNHy3MZ02bCaF\n4qi2I+KD0XlIhXAKxyby2fTXOvMH0pxd2px8RoYU3Uwb9tEYQJtUOZrZdER8krHFhds4NpEDTOZu\nhuiMsaDHSEKWpC2EiEeO7LXisYMWTOIeu8MQCcYzAbUQbuHoRA7wExfxCYPtDkPEkT001W5lh0P+\nRRsjTgj3cHTTiscWLrE7BBEHzlGRV3jKr/weJmkDoy2mR6zDEiJqrkjkoMx8LZM6iGis5TLdcmm+\nE27n+KYVj4kMtTsE4WKv8oTjxrgXwiyuSeT7aBxHs36IWJlDP8byZNxf2ScSW5JN+y1FnVE+XEmU\ncjE7uIuPTA5JxJMZ3MxWWnOWSnaHIoSJckAnb7sukXtrxTaqcSKuh6cU4fuCG1nB5XaHIYQF9BN5\nqJOdlYBFQEWUeTg/B0YDfwHuBw6r6z0NfKXeHw0MA4qBh4G5UcUdhOdS6pV0ASBLndC5J4us2qVw\nsFNUZSmZksRFwjFSI68CnEZJ+t8BfwauA06A36wLbYGPgCuARsB8oBX4dTcxpUZuxABmAcocftU5\nHpN9iugdoQ57aAoo87rODzizoBCJJLIaOShJHJQaeQpwTF3W+xIYCEwBzqNMyrwD6AosCytWE81i\ngM/yECYDMg+nk+zgYubQj6PUtjsUIVzJSCJPBlYBLYB3gI3AIOBPwL3ACuBx4BegIb5Jex84azT3\nyQzR7qdTwC18FtczhzjRr1TgGDV9ZnYRQkTOSCIvAToC1YGvgSyUhP5X9fHngdeA4QGeXxpdiNYp\nIN0nmQxnPLU5CkB5zsuUXyYrohxr6MgX3Gh3KELElXCu7DwOfAFcDupZRcV7oDZEw37wGRy8sVqm\nw3sTGerNXuO9vouSKSGZEi5lPTcx08ao4sfLjPYb20QIEUyeegsu1MnOOkARSrNJZZQaeQ5K80q+\nus5IlJObd3HhZGdXLpzsvBj/WnnMTnaaJYVi0ijUlh/mLe1+Hk74GjJXHuYd09dczxo62tynO4/4\nepXyiK/jgfg7pjzMP57ITnY2AD5AaSdPBiYBC4D/oDS3lAK7gAfU9TcBU9W/RcBDOLhpJRzFpPhM\nG5WjfhHV5ihFLOZKqpDJUrvCM10e5rwFv+Z6lnGVCVuKVh6SJJwuj/g6pjxidTyhEvl6oLNO+b1B\nnvOSeksISk+LGswli2VcRUfW0IuFdoflGM5I4kLEN2mwNFEhafzERXaH4RhvM8LuEIRICHZdor8G\nAowpKoQQIpBFKD0HhRBCCCGEEEIIIYR79QO2ANuBUTbHYtQEoAClF49HLWAesA1lhEfvWXtHoxzf\nFqBvjGIMRxNgIcr1ABtQRqkE9x5TJeAHlHMvm4CX1XK3Ho9HCrCaCxfcuf148oB1KMe0XC1z8zHV\nAKYBm1Hed1fi7uMxLAVlIK0MoDzKB6+NnQEZ1APohG8iHws8qd4fBbyi3m+LclzlUY5zB87rHVQf\n5ToAgFRgK8rr4OZj8kwBVA5lvJ+rcffxADwGTAbt0mK3H88ulETnzc3H9AHKkN2gvO+q4+7jMSwT\nmOO1/JR6c4MMfBP5FiBdvV9fXQblW9f7l8YccHxn6hlAb+LjmKoAPwLtcPfxNEa5MroXF2rkbj4e\nUBJ52SEu3XpM1YGfdMptOZ5YfyM0AvZ6LTtudMQwpIM2bGIBF168hijH5eH0Y8xA+bXxA+4+pmSU\nGk8BF5qN3Hw844An8B3L383HA8pV3vNRRkz9vVrm1mNqjjKxzkSU0WH/DVTFpuOJdSKPi8v1dZQS\n/NicetypwHTgEZSJQry57Zg8o3Q2Bq5Bqcl6c9Px9AcOobQlB7rWw03H49EdpdJwAzACpcnSm5uO\nqRzKVe//UP+ewr91IWbHE+tEXnZ0xCb4fku5SQHKTydQxqQ5pN4PYwRIW5VHSeKTUJpWwP3HBBdG\n6eyCe4+nG3ATSlPEFOBalNfJrcfjcVD9exj4DGVwPbce0z719qO6PA0loefjzuMJSzlgJ8rP+Qq4\n52Qn+LeRj+VCm9dT+J/UqIDy82sn9l1BG0gSysBn48qUu/WY6nChd0Bl4FuU6QjdejzeenKhjdzN\nx1MFqKberwp8j9Jzw83H9C2oEwcr8xiPxd3HE5YbUHpJ7EA5AeAGU4ADwK8obfxDUc6+z0e/m9HT\nKMe3Bbg+ppEaczVKU8QalJ/vq1G6hbr1mC5Faadcg9K97Qm13K3H460nF3qtuPl4mqO8PmtQurx6\nPvtuPqbLUGrka4FPUU6Auvl4hBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBCx8P/ryFQt\nI27JTgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fdc3df5c550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mask2 = cv2.medianBlur(array(mask, dtype=np.float32), 5)\n",
    "mask2 = cv2.erode(mask2, ones((3,3)))\n",
    "imshow(mask2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f00c27df5d0>]"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imshow(images[4], cmap=cm.gray)\n",
    "plot(points[:, 0], points[:, 1], 'r.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 602,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f098ee01990>]"
      ]
     },
     "execution_count": 602,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGkNJREFUeJzt3X2QHOV94PHvoJVAQi8rIRu9mlVAGCTHSGhBwii3kzvA\n0voqkiulABU7YEgURznjq3sBBKnT+lJ35vCRgK8KUro4YGxLiZykODBIRnIxQOwDEYxAlizr5SKd\nV1cIZLTgtxgJ+v54ujXdq1mtdqd352W/n6qu7n66p6ef3Zn+TT8v/YAkSZIkSZIkSZIkSZIkSZJU\nt84BXgR2ALuBL8bpU4CtwF7gaaA19Zq1wD5gD3DdsJ2pJGnYjIvnLcALwFLgXuD2OP0O4J54eR4h\niIwG2oD9wFnDdaKSpOE1DngJmE+4Gzg/Tp8Wr0O4W7gj9ZotwJLhOkFJUv/y+LV+FuEu4AjwDLCL\nEBSOxNuPUA4SM4Du1Gu7gZk5nIMkKSctORzjfWABMAn4NvCbvbZH8dSX022TJA2zPAJD4m3gSWAR\n4S5hGvA6MB14I97nMDA79ZpZcVrGhRdeGB04cCDHU5OkEeEAcFG1B6m2KGkq5RZHY4FrgVeAx4Gb\n4vSbgMfi5ceBG4AxwBxgLrC990EPHDhAFEVNO61bt67m52DezJ/5a74JuLDKazpQ/R3DdOCrhABz\nFvA14DuE4LAJuBU4CPxOvP/uOH03cAJYg0VJklRXqg0MO4HLK6S/BVzTx2v+azxJkuqQfQhqoFgs\n1voUhkwz5w3MX6Nr9vzlpVDrE+hDFJeXSZLOUKFQgByu694xSJIyDAySpAwDgyQpw8AgScowMEiS\nMgwMGrTVq6FYhM5O6Omp9dlIyouBQYO2dy88+yxs3hyChKTmYGDQoI2Lh2hqb4f162t7LpLyYwc3\nDVpPT7hTWL8eWlv731/S0Mqrg5uBQZKahD2fJUlDwsAgScowMEiSMgwMkqSMagPDbOAZYBfwA+C2\nOH0KsBXYCzxNefhPgLXAPmAPcF2V7y9Jylm1tdfT4mkHMB54GVgJfAY4CtwL3AFMBu4E5gEbgCuA\nmcA24GLg/V7HtVWSJA1QvbRKep0QFAB+BvyQcMH/LcJY0MTzlfHyCmAjcJwwFvR+4Moqz0GSlKM8\n6xjagIXAi8D5wJE4/Ui8DjAD6E69ppsQSCRJdaIlp+OMB/4O+Dzw017bonjqS8VtXV1dJ5eLxaJj\ntUpSL6VSiVKplPtx8+j5PBr4FrAZuD9O2wMUCUVN0wkV1JcQ6hkA7onnW4B1hLuMNOsYJGmA6qWO\noQB8BdhNOSgAPA7cFC/fBDyWSr8BGAPMAeYC26s8B0lSjqqNLEuB54DXKBcJrSVc7DcBHyJUMv8O\nkDyx/y7gFuAEoejp2xWO6x2DJA2QD9GTJGXUS1GSJKnJGBgkSRkGBklShoFBkpRhYJAkZRgYJEkZ\nBgYNyurVUCxCZyf09PS7u6QGYmDQoOzdC88+C5s3hyAhqXkYGDQo48aFeXs7rF9f23ORlK+G7vm8\nenX45TpuHGzYALffnl1vbe33EBqknp7w91+/3r+zVC98JAahjPvZZ8PynDnw1lvw9tun7jdqFIwf\nD1dcAd/8phcySc3JR2KQLc6YMaNyUAB4772wbds2mDIF2tpCcDjnHPjAB6BQKE8LFlihKmlka+jA\nsGEDrFoFW7fCxIkhraWfoYeiCA4dCoHiV7+Co0ez2199NVSoTp5cDhabN4dAkg4gZ50FY8bAa68N\n7JxtzSPlz+9Vvhq6KCktKfP+0pfgU5+Cf/iHkL5kSbhL2Lw5BIVCIQSRt98OF/f33z+z448aFe48\nBiMJJK2tcNFF8GI8LNGqVbBp0+COKY0kq1fD178Ox4/DiRPl9AkTwvfyn/+5/F0eyd8r6xgG6LXX\nYPHicFGeNAmWLoUnn4R16+DNN+G73+37tU89BStWhA/lcLj6avjWt6wL0ciwejU88QS88Ub48Zb+\n6n/84+Gi/9prcOxY/8caPTocZ6R+dwwMOTt0KNRVJEVL48fDz34WgsLy5eGDedllw3pKA7ZoUahH\nSb4UvVtt1fLLUk/nosFL/o8HDsAFF4S779O1CLzkEvjRj8qvPysuvE5+3be0hO/Vyy8P/pySu/lR\no+D734ePfnTwx2p09RQY/gr4BPAG8Otx2hTgb4ALOHUEt7WEEdzeA24Dnq5wzLocqKenB26+GX7x\nC3jlFXjnnfDFuOgieOEFWLgwpFdT7FQLo0eHL/S554bWXRMnhkr5Q4dCenp5sBf1dAuykXyrX++S\nX++/+lX4fEfRmRe39jZhQviu9PddOF2R7qJFIWgsWBA+h7t2wf33ww03wFVXhbQvfAE+8YlQfHzB\nBYM712ZRT4HhN4CfAY9SDgz3Akfj+R3AZOBOYB6wAbgCmAlsAy4Gen8s6jIwDMShQ/Cxj8H8+fDL\nX5brPCBUbP/853D55eFDP1xFVMOhUAiBZuHCULeTBJVdu8LdWHt7aCxwpnc13mn0LX0RX7QotMx7\n+uns+je+Ud8/UgoFeO45uO++8J3YuRNmz4aXXgp361ddZX+ZgcgrMOSlDdiZWt8DnB8vT4vXIdwt\n3JHabwuwpMLxopHk4MEomjYtij74wSgqFpNSVqe+pkIhiqZMiaKzzgrrLS3h73f11VG0fHkUHTsW\nRX/wByFt8uQouuaaKPq934uijo4omjUr/J2T9GPHsv+LD384ikaPDsfs6AivGzs2rE+ZEv5XvSXv\nNWZMFE2aFEUzZmTPJb1fR0f2HEeNqv3fcyin+fOjaOrU8LeZOjWKOjuj6Pnnw/rZZ0fRq68O4Rdr\nBAKifC/t1WkjGxjS1USF1Pr/AH43te0vgd+ucLxa/31r6uDBcAFLLkIHD4aLzZIl4aI4dWoUbdoU\nLoyLFpUvLs1+kXGqzTR6dPjcLV4cgurEiafuM39+mE+cGD6XV18dRStWnBp4NbQgn8DQT6v/XPR3\nshW3dXV1nVwuFosUi8VcT6qW+iseueAC+PGPs+uHD596nFWr+n6PpPnun/xJqDxfsAB+8pNyU9nE\nqFFhXs/FDcrHuefCu++Wiy4nTgz1CIlkffFiOO+8sN+uXfC971l2X69KpRKlUqnWp9GnNk4tSpoW\nL0+nXJR0ZzwltgCLKxyv1oF3SHV0lH9prVo1fO977Fh4v4MHw/xMf80dOxaKAKZOLf96TM4/Kc5x\nGtx07rmhWCX9N620z+TJYXnevNMfb8qUKFq6NBxzyZJQXFMsln+9J58Bf8k3J3K6Y8hLG9nAkFQ6\nQwgE98TL84AdwBhgDnCAyhUltf77Dqnly8OXuL298he0d1l0o0pfhJLiseefL9endHaGMuYZM8LF\nKymLnjw5lOknF7tx48K80sWzd2BK9j3dlH5NsjxlSjiHD34wnMd552XPISk+mTQpij7wgSgaPz6s\nt7SEYpNk+9lnl4/b2hrOub095Ku9PaRv2hTy3NlZ/v8mxYWdneXAXSmADza4a2Qgp8CQR+31RqAD\nmAocAf4T8L+ATcCHOLW56l2E5qongM8D365wzDiPzam/J5PatFPSYNRTc9Wh0NSBoT+dneERHr2b\ndkrS6RgYmphjHUjVG4l9YAwMknQaI7FI1vEYJOk0HH528AwMDeqSS8KtcfLICUlZ6fFaRkIxUp4s\nSmpQra3lEetmzcp2iJM0MlmUNMKNHh3m48ZlH9AnDTVHS2t+BoYG9Y//GO4Udu/2cQUaXnv3hkrd\nzZtDkFD16i3YDsezkjQEej9PSRouVurmLwm2EIJErVtQGRgaWLqddh6D6UhnYsMG+9nkrd6CrZXP\nDSzdTvvss8MALRDGp37ssYF18Olr396DwXzzm14MpLzl1anVymdlfmWMHVtOL8Qfi4GUBfe17969\n8PrrYSD2bdvqv0y53spqpTPR2hqKj+rlR5eBoYGl22m3t4e0hQvh4YfD8kBuT5N9x48PQSC5qCbp\nybHr4Tb3dKwYlapnUVKTqHQrOpDb054euPhiePPNsJ48QqCnB26+OdyFPPxw/fyi6YsPINRI5rOS\nlLtmuKj6AEKNZAYG5c6LqtTYDAySpIxGb5W0jDAO9D7KQ4BKkupALe4YRgE/Aq4BDgMvATcCP0zt\nM+LvGHr3K7j99ubozDYSB0+Rhktedwy16Pl8JbCfMBY0wF8DK8gGhhGvdxf5N94or0+dCkePlrfV\nuvv8QNRb139Jp6pFUdJMIP2Un+44rW711Wkq785UyfFmz4YdO0Ja0ncg3SdhwYLycr33K+it3rr+\nDzU73OlM9XedGT06NBsvFGDp0qH9PNXijuGMyoi6urpOLheLRYrF4hCdTv/6+pWb96/f9PESR4/C\nypXhQ7FiBTzySPn9GrH10Eh6zk7ymUjGzfAOSadzJteZxHe/G/ZZs6ZEqVQa1vMcKkuALan1tZxa\nAR3Vk+XLowiiqL09io4d6z+92vcZNSrMIYrOO6+8vGpV9e+h4dPRUf7fTZ6cz2dEzau/60x6mj+/\n8ueJM/zhXY9agANAGzAG2AFc2muf4ftvnIFjx8JFufc/oq/0at8nuaAsXBhF11yTb/DR8EkH+o4O\n/386vf6uM88/H0XnnBNF117b92eJnAJDrfoxLAfuJ7RQ+grwxV7b4zyOTOmOZjByil6aTU8PzJ1b\nbiiQPGZEGip2cBshbN7Z2JrhMSNqHI3ewU1nyKeFNrb0E3ANCmoUjuBW50Za885mkzxnX2okFiXV\nOR9s15gsAlQtWMcg1bH0sKtWOmu4WMcwgth7tvFYBKhGZmBoAFZANx4rndXIrHxuAP76bDxWOquR\nWcfQAKyAlnQmrHyWJGVY+SzAimlJ+TMwNDgrpiXlzcDQ4KyYlpQ36xganBXTkhJWPqsiH8UgjVxW\nPqsi6xwkVauawLAK2AW8B1zea9taYB+wB7gulb4I2Blve6CK91YfrHOQVK1qAsNO4JPAc73S5wHX\nx/NlwIOUb20eAm4F5sbTsireXxX4KAZJ1aomMOwB9lZIXwFsBI4DB4H9wGJgOjAB2B7v9yiwsor3\nVwXJoxgMCpIGayjqGGYA3an1bmBmhfTDcbokqY709xC9rcC0Cul3AU/kfzqSpFrrLzBcO4hjHgZm\np9ZnEe4UDsfL6fTDfR2kq6vr5HKxWKRYLA7iVCSpeZVKJUqlUu7HzaMfwzPAfwBejtfnARuAKwlF\nRduAi4AIeBG4jVDP8CTwZWBLhWPaj0GSBqge+jF8EvgxsIRwkd8cp+8GNsXzzcAaQlAgXv5LQnPV\n/VQOCpKkGrLnsyQ1iXq4Y5AkNSEDgyQpw8AgScowMEiSMgwMkqQMA4MkKcPAIEnKMDBIkjIMDJKk\nDAODJCnDwCBJyjAwSJIyDAySpAwDgyQpw8AgScowMEiSMqoJDF8Cfgi8Cvw9MCm1bS1hlLY9wHWp\n9EXAznjbA1W8tyRpiFQTGJ4G5gOXAXsJwQDCmM/Xx/NlwIOURxR6CLgVmBtPy6p4f0nSEKgmMGwF\n3o+XXwRmxcsrgI3AceAgYWznxcB0YAKwPd7vUWBlFe8vSRoCedUx3AI8FS/PALpT27qBmRXSD8fp\nkqQ60tLP9q3AtArpdwFPxMt3A+8CG3I8L0lSjfQXGK7tZ/vNQCfwr1Jph4HZqfVZhDuFw5SLm5L0\nw30duKur6+RysVikWCz2cyqSNLKUSiVKpVLuxy30v0uflgH3AR3A0VT6PMLdw5WEoqJtwEVARKiL\nuI1Qz/Ak8GVgS4VjR1EUVXFqkjTyFAoFqO66TrUH2AeMAd6K1/83sCZevotQ73AC+Dzw7Th9EfAI\nMJZQJ3FbH8c2MEjSANVDYBhKBgZJGqC8AoM9nyVJGQYGSVKGgUGSlGFgkCRlGBgkSRkGBklShoFB\nkpRhYJAkZRgYJEkZBgZJUoaBQZKUYWCQJGUYGCRJGQYGSVKGgUGSlGFgkCRlVBMY/hR4FdgBfIfs\nOM9rCSO87QGuS6UvAnbG2x6o4r0lSUOkmpF+JgA/jZc/B1wG/D7lMZ+voDzm81zCmM/bgX8Tz5/C\nMZ8lKTf1MILbT1PL44Gj8fIKYCNwHDgI7AcWA9MJwWR7vN+jwMoq3l+SNARaqnz9fwE+DfwSuDJO\nmwG8kNqnm3DncDxeThyO0yVJdaS/wLAVmFYh/S7gCeDueLoTuB/4TF4n1tXVdXK5WCxSLBbzOrQk\nNYVSqUSpVMr9uFWXRcU+RKgz+AghSADcE8+3AOuAQ8AzwKVx+o1AB/DZCsezjkGSBqge6hjmppZX\nAK/Ey48DNwBjgDnxftuB14F3CPUNBUIR1GNVvL8kaQhUU8fwReDDwHvAAeCP4vTdwKZ4fgJYQ2iR\nRLz8CDCWcIdRqUWSJKmG8ipKyptFSZI0QPVQlCRJakIGBklShoFBkpRhYJAkZRgYJEkZBgZJUoaB\nQZKUYWCQJGUYGCRJGQYGSVKGgUGSlGFgkCRlGBgkSRkGBklShoFBkpRhYJAkZeQRGP498D4wJZW2\nFtgH7AGuS6UvAnbG2x7I4b0lSTmrNjDMBq4FDqXS5gHXx/NlwIOURxR6CLiVMA703Hi7JKmOVBsY\n/gy4vVfaCmAjcBw4COwHFgPTgQnA9ni/R4GVVb6/JCln1QSGFUA38Fqv9BlxeqIbmFkh/XCcLkmq\nIy39bN8KTKuQfjehHiFdf1D1ANRpXV1dJ5eLxSLFYjHPw0tSwyuVSpRKpdyPO9iL+UeA7wC/iNdn\nEe4AFgOfidPuiedbgHWEeohngEvj9BuBDuCzFY4fRVE0yFOTpJGpUChADj/SB1uU9APgfGBOPHUD\nlwNHgMeBG4Ax8ba5hHqF14F3CMGjAHwaeKyKc5ckDYH+ipLOVPrn/W5gUzw/AaxJbV8DPAKMBZ4i\n3E1IkupIrvUCObIoSZIGqNZFSZKkJmVgkCRlGBgkSRkGBklShoFBkpRhYJAkZRgYJEkZBgZJUoaB\nQZKUYWCQJGUYGCRJGQYGSVKGgUGSlGFgkCRlGBgkSRkGBklSRjWBoYswpOcr8bQ8tW0tsA/YA1yX\nSl8E7Iy3PVDFe0uShkg1gSEC/gxYGE+b4/R5wPXxfBnwIOURhR4CbiWMAz033i5JqiPVFiVVGkJu\nBbAROA4cBPYDi4HpwARge7zfo8DKKt9fkpSzagPD54BXga8ArXHaDEIRU6IbmFkh/XCcLkmqIy39\nbN8KTKuQfjehWOg/x+t/CtxHKCbKRVdX18nlYrFIsVjM69CS1BRKpRKlUin341YqChqMNuAJ4NeB\nO+O0e+L5FmAdcAh4Brg0Tr8R6AA+W+F4URRFOZ2aJI0MhUIBcriuV1OUND21/ElCayOAx4EbgDHA\nHEIl83bgdeAdQn1DAfg08FgV7y9JGgL9FSWdzn8DFhBaJ/0T8Idx+m5gUzw/AayJ9yFefgQYCzxF\nuJuQJNWRvIqS8mZRkiQNUD0UJUmSmpCBQZKUYWCQJGUYGCRJGQYGSVKGgUGSlGFgkCRlGBgkSRkG\nBklShoFBkpRhYJAkZRgYJEkZBgZJUoaBQZIGafVqKBahsxN6emp9NvkxMEjSIO3dC88+C5s3hyDR\nLKoNDJ8Dfgj8gDBwT2ItsA/YA1yXSl9EGOltH/BAle8tSTU1blyYt7fD+vW1PZc8VRMYfhP4LeCj\nwEeA/x6nzwOuj+fLgAcpDxzxEHArYbjPufH2EWcoBu+uF82cNzB/jS7v/G3YAKtWwdat0Nqa66Fr\nqprA8EfAF4Hj8fqb8XwFsDFOPwjsJ4zzPB2YQBj/GeBRYGUV79+wmvnL18x5A/PX6PLOX2srbNrU\nXEEBqgsMc4F/AbwAlID2OH0G0J3arxuYWSH9cJwuSaojLf1s3wpMq5B+d/zaycAS4ApgE/BruZ6d\nJKmhbAY6Uuv7ganAnfGU2EIoSppGqKhO3Aj8RR/H3g9ETk5OTk4DmvZTY38IfCFevhj4v/HyPGAH\nMAaYAxygXPn8IiFIFICnGKGVz5LUrEYDXyM0P30ZKKa23UWIXHuAj6fSk+aq+4EvD8tZSpIkSWoe\nywh3GfuAO2p8LoM1G3gG2EXo+HdbnD6FUJm/F3gaSDdw66tDYL0aBbwCPBGvN1PeWoG/JdSH7SYU\nfTZT/tYSPps7gQ3A2TR2/v4KOELIT2Iw+anXzreV8vclwufzVeDvgUmpbY2Wv36NIhQxtRGKqXYA\nl9byhAZpGrAgXh4P/IiQj3uB2+P0O4B74uWkTmY0Ie/7qf9Hlfw74BvA4/F6M+Xtq8At8XIL4UvX\nLPlrA/4PIRgA/A1wE42dv98AFpK9cA4kP0n953bgyni5nuo/K+XvWsr/h3to7Pz16ypCC6ZE79ZN\njeox4BpCBD8/TpsWr0OI8Om7oy2EJsD1ahawjdDzPbljaJa8TSJcOHtrlvxNIfxQmUwIek8QLjKN\nnr82shfOgeZnOtkWkzfQd4vJWmgjm7+0TwJfj5dzy189Rf+ZwI9T60nHuEbWRoj2LxI+qEfi9COU\nP7h9dQisV38O/Efg/VRas+RtDqEH/8PA94H/CZxL8+TvLeA+QgvC/wf0EIpcmiV/iYHmp5E7395C\nuAOAHPNXT4EhqvUJ5Gw88HfA54Gf9tqWtDnuS73+Lf418AahfqHQxz6NmjcIv6IvJzzf63Lg55x6\n19rI+bsQ+LeEHywzCJ/RT/Xap5HzV0l/+WlkdwPvEuqKclVPgeEwoeI2MZtslGskowlB4WuEoiQI\nv1ySXuTTCRdYODXfs+K0evQxwoMT/4nwPKx/SchjM+QNwuetG3gpXv9bQoB4nebIXzvwPeAnwAlC\nxeVVNE/+EgP5PHbH6bN6pdd7Pm8GOoHfTaU1U/5OaiF0hmsjdI5r1MrnAuEBgX/eK/1eyuV/d3Jq\nhVGlDoH1rINyHUMz5e05QodNgC5C3polf5cRWsqNJZznV4E/pvHz18aplc8DzU89d75tI5u/ZYSW\nZVN77deo+evXckLl2H5CRUojWkoof99BKHJ5hfBPmEKotK3UhK6vDoH1rINyq6RmyttlhDuGdFPA\nZsrf7ZSbq36VcHfbyPnbSKgveZdQR/kZBpefeu182zt/txCanB6ifH15MLV/o+VPkiRJkiRJkiRJ\nkiRJkiRJkiRJkqSy/w9ABQOFWnK/zwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f098eff8750>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(mv[:,0], mv[:,2], '.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1614,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.9389954384692558e-05, -0.035636072284646048, 29.733696485112709)"
      ]
     },
     "execution_count": 1614,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn import linear_model\n",
    "\n",
    "X = np.vstack((mv[:, 0]**2, mv[:, 0], ones(mv.shape[0]))).T\n",
    "y = mv[:, 2]\n",
    "\n",
    "# Robustly fit linear model with RANSAC algorithm\n",
    "model_ransac = linear_model.RANSACRegressor(linear_model.LinearRegression(fit_intercept=False))\n",
    "model_ransac.fit(X, y)\n",
    "inlier_mask = model_ransac.inlier_mask_\n",
    "outlier_mask = np.logical_not(inlier_mask)\n",
    "\n",
    "a, b, c = model_ransac.estimator_.coef_\n",
    "a, b, c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1615,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18.931283939195769"
      ]
     },
     "execution_count": 1615,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib qt\n",
    "\n",
    "X_in = mv[inlier_mask, 0]\n",
    "X_out = mv[outlier_mask, 0]\n",
    "y_in = mv[inlier_mask, 2]\n",
    "y_out = mv[outlier_mask, 2]\n",
    "plot(X_in, y_in, '.')\n",
    "plot(X_out, y_out, '.')\n",
    "t = linspace(X_in.min(), X_in.max(), 100)\n",
    "plot(t, a*t**2 + b*t + c)\n",
    "x_min = -b/(2*a)\n",
    "a*x_min**2 + b*x_min + c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1672,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f098b0f6650>]"
      ]
     },
     "execution_count": 1672,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "control = [10e3, 15e3, 20e3,\n",
    "           25e3, 26e3, 27e3, 28e3, 29e3, 29.5e3,\n",
    "           30e3]\n",
    "speed = [9.740644336872279/3, 36.030645952209511/3, 37.124312115088088, \n",
    "         101.34249310780383, 123.63380553013994, 80.687738743971806*2, 95.34390238190602*2, 110.64445509500386*2, 78.85052787064636*10,\n",
    "         162.92974488164242*10]\n",
    "plot(control, speed, '.-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 571,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: pylab import has clobbered these variables: ['diff']\n",
      "`%matplotlib` prevents importing * from pylab and numpy\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8VFX6+PHPnZKZTHpCQgKBUEMvQkAQlLj2iGuNuliw\nrNGfbnFh7SDIsjbAritYUL8aFcQGigorsICIChilhFAkJkAIkElIn3Z+f5xJJSglkMA8b173xcy9\nd+7cO5l57rnnnPscEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBCt2OvAHuDnJpaNB3xAdL15DwBb\ngGzg/OO+d0IIIZrdmcBpHBz4OwBfAL9QF/h7Az8CVqATsBUwnZC9FEIIcdh+LzAvB5xNzH8KuLfR\nvEuBdwE3sAMd+Ice4/4JIYRoZkdTIr8UyAd+ajS/nX9+jXyg/VHulxBCiOPEcoTrO4AHgfPqzTN+\nY311xHskhBDiuDrSwN8VXX+f5X+eCKwBTgd2ouv+qbds50Eb6NpVbdu27Yh3VAghAtw2oNuJerNO\nNN2rB5pu3A0COqN3sqmrAXUymzRpUkvvwjGR/W9Zsv8t52Ted6WUohlrUH6vjv9d4BsgGcgDbm4c\nxOs93gjM8f+/ELgTqeoRQohW5/eqev70O8u7NHr+qH8SQgjRSkk/+yOUmpra0rtwTGT/W5bsf8s5\nmfe9uf1Wj5zjxV9dJYQQ4nAZhgHNFLOlxC+EEAFGAr8QQgQYCfxCCBFgJPALIUSAkcAvhBABRgK/\nEEIEGAn8QggRYCTwCyFEgJHAL4QQAUYCvxBCBBgJ/EIIEWAk8AshRICRwC+EEAFGAr8QQgQYCfxC\nCBFgJPALIUSAkcAvhBABRgK/EEIEGAn8QggRYCTwCyFEgJHAL4QQAeb3Av/rwB7g53rzpgGbgCzg\nQyCi3rIHgC1ANnB+8+2mEEKI5vJ7gX82cGGjeV8BfYABQA462AP0Bq7x/38h8NJhbF8IIcQJ9nuB\neTngbDRvEeDzP14NJPofXwq8C7iBHcBWYGiz7KUQQohmc6wl8luAz/2P2wH59ZblA+0PZyMenw+l\nFO3Of4+ImEpiYyE3t255RgakpkJaGhQXH+MeCyFOGhkZ4HCAxQJ2OwwbJnGgOViO4bUPAS4g8zfW\nUU3NnDx5cu3j1NRUvk5KYld1NXtj5+O5/Fl491NGjowlL0+vk5MDy5bpxxkZMGfOMey1EOKkkZMD\nlZWA2YX3wjtYnX0ZrP5jQMSBpUuXsnTp0uOybeMw1ukEzAf61Zt3E3AbcA5Q5Z93v///x/3/fwFM\nQlcH1aeUang+KPV4uHbjRhYuUqjFH0LyO4R8/DnDk5OZOxfGjIGFCyElBRYtgsjIwz08IcTJLDQU\nyr3FcPWV4A6FDzJJGRASkHHAMAw4vJj9u46mqudC4B50nX5VvfmfAtcCQUBnoDvw3eFsMNRsZl7P\nJGJdwXDVaPh5IuVXn8niLcvIyIDMTEhPl6AvRKAxonbALSOgsB+xX39I+mWBGfSb2++dPd4FRgFt\n0N06J6F78QQBRf51VgF3+h8/iK739wB/B75sYpsHlfizLv8HJVe8i/3+m3lrSBAvXXc+6j9VMHgM\nLH6CURE38fHH8scWIpB8k/cNI5+/ErXsQRzr/8rGjZCU1NJ71XKas8TfLBs5QgcF/nWRb1By1k9w\n0xtYJ/2VstzObOjensjg//BI9x/4Ne86/hj6bz75WHqHCnGqy87IJndtLhvKNtDmmV48dNtFrFgR\n2EEfWr6qp9mZrF6Y/0fsj12He8rT2IblMGitlb3e27jKdwl0XMl3Xa6gzFXW0rsqhDiOfMrHptWb\nsK2xMWjzIAa/3pG8PAn6za1VlPjdufvIOeND7tj1MxuMt3AHl9IjMo7Ju56GkFg+vOU70j94mT2h\npXTa3gMX0fRcm0nH/lL3I8SposxVxg0f3cBF/76I5J+TCU0JZcCiAax/4Xnc24OxFZxB78zeWCOt\nLb2rLeKUq+oB8PnAZl6Fhx3oZoIqUvvEMGnb/1Fmgp8HVZP67U4cnlJ6MxUvLqqztnDPC0nk5Oi+\nvpmZ0g4gxMnoF+cvXPrepfzlo78wqHwQ1durOW3VQPaoJ/g161XUP6fAtm7EpsfSZ06flt7dFnHK\nVfUAmExg51J01/9lOOjK+LMHwQt/ITJ8O302V1Dt6YOTYeQwHjtu3KePrO3jv3Ch7uMvhDi5/Hf7\nfxn+2nD+POjPDKkeQtk3ZbiLD/Dj15fidH5N+NvvoLZ1Y0mnJDo+l9zSu3tKaDWBH2DQgF7AdQSx\nEC+rGPdCLPc96KH4yYcIb18AwK8dPHxz4UZ8wE5LGfcuu4AIijGbYcKEFt19IcQRUErxzLfPcN2H\n15F5ZSbnv3Y+5T+XQ8IueOEvuLKDODDsYfYvDmJauwEsiEyiwgjMap7m1qoC/ydLPyE9PZ38HXeS\nGPI4W3ia7wpHMeNvSTBxKlwxj455Hp655mJuHz+entXlpPEVMy034PXCxRe39BEIIQ5HhbuC6z+6\nnrey3uLbP3/LHzr/gcqcSjydV8MLf4EFo2HaPex1h3FXSX88Tg8rvjHRtm1L7/mpodXU8TfWoUMH\n8vPDgY+BBdiZiC20klev70nS7JvZ0KEzr9+xjzefmEho6T5SYt7h151/IiIChgyBuXOlvl+I1mhb\n0TaunHMl/dv2Z+bomQRbg9mUsYm97hfwXfIu/PthWDuQn4hgCr25yrKTp7YkcPejfycnJweHw0Fm\nZiaRAfYDPyUbdxsbOXIkK1euBCKBd4BgdNbnvRhAf2svxsX9P0yqF+13gYFiXXgFEw+kUY6VoCCd\n50O6gQnRsjIyqO2Acf3UT/nH17fx8FkPc+eQOzEMA4+njNUvX4Hb8itMegSjKoG7SjewQV2BwVii\nIr9n7Y9rGTt2LMv8SbvS09OZc6on62nklGzcbSw8PNz/qBi4BFgJfA8MQQFZ7k18UHE37QvcGFgA\nK6cdiGB8h9lgL8blgpEjW2jnhRC1cnJg2XIPC10PcNu8v/LptZ9y19C7MAyDiorNrF07DMPngL89\nh6u8PY8eSGabOgc4A8WXFBUXMXLkSBwOBwApKSnMmjWrZQ/qJNdqA39mZiZpaWn+s5wPmIjOArEA\nuB2ArTGK8qv/D0xeANwmL2fnted/VUvIozvm/FxSUyWFqxDHW+PU6RkZkJAANhss/3EX3HgOr+/6\nP/Z93R565LE84n/8b/RU1v5wJomJfyflhjlUnduJvxqDWU82FuMsYDsADoeDFStWkJmZSXp6OosW\nLQq4ap7m1mqremrk5ubSrVs3PB6Pf0539IiPP6BTBFXy6t296Pr6o6gDERj+Q7KwjzNIZzVDcRLD\ndWQycnQkpaX6kvOj2AxsuXIDgBDNITW1LnV6mzbg9YLTCXRZzMyQP5KcG4fVeyuO8jBKLD0g43UY\nuQIencyI5bfwxXIrt94KD4z38s60YazZ/wMAdrud7OxskqTONjCqemokJSVx3nnnAdDT3BMz24HT\n0UMJfAt0587nN8Os2zH6rq99XYU9GAyD4XxHGgt5mQwWLKjr87/1c7kBQIjm4q+FITQU9u0DZ4kH\nznkQLhtLty3JfJc2nB+ThlESlwDPjod2u+D2mXjX9+T/DSvmrrvg44/hH/eZiRsaC+gqnd27d0vQ\nPw5afYkfoLi4mIyMDCqXV5JTkEMOOf4lGcBU4C+AbugxGwYz1cu0t3Zld1IpFxbexfaqMZS6+2FX\n1VTwOTeYZvPL8DGEr5Qk/0I0h5rqHacT/rn4I2xYUYaXcV26Ysr4jms+iWSosR7GPYPx0TWoOVez\n12XjsdB+xAwJIfN9E7GxNdvSv/dZs2ZJlU49AdGrpylnjTyL5SuXA2BgoFDAaeigvxj4B1BFkMXg\ny7APwBnNt8OqOG19AbayTgDEsoT1bGJKu1lsGJEBs2ZJ0Beimbz+3ft0OT0ECAVgf4yPf90fzYsF\nD+BOWQdTJ8CGvmQPT2LiukT+Ot7EhClmTK2+7qHlBVRVT32h4frLlJKSQpGziMTERGAdMAjd7fNb\noAdmZaJg6q3caLuaf317JVeX3U0BBVgoppo2tOVP3LUrh8fnPkT7pBByc6FnTx3/G4/3K4T4bRkZ\n0LZjKeYrbuG2zAno4TpAoYhO2MUzUdfijiyDjFl4NvTjz+zlb6scdOrxAH/7Z6kE/RZwUpX4G18C\n1vX1r3Eb8G/gPmB2g9dGEcVcYx5m1fCQlxDLrPg+VFZCSYmel5hI7Xi/QohDy87IZsW8PcRXf0+o\n9T9kRwyl3f4/E15mggu+hDtehjdugk8uZTfB/Ive5JNFKVcA+wKyP/7RCtgSf2RkJHPmzKmt96vr\n61/jFSAVGIcePCyidolh8tUGfZ+hu3+6rNuZYY8DwOpPAeJwwIoVx+0QhDhluL1uNq3eRLciRWh5\nCiWO+5j4zDU4rEXw8BS45n34x9PwyWUspi13GoM4m0J68k9gH4P6D5L++C3kpCrxN1ZcXEyXLl1w\nOp2NltiBaegbv24AdLuACZjFTP57p42r3m1LoSoj0bODUPdculTv4leSKCWcGYMzCYqLlF6eQtST\nnZFN1vxKDlSbeCzexP3qbAZuvpcKNZRdieVE7bVh752Fcd+TsPxMeOU2yl0OnjW6syU4ksw54Hhp\nMwVrC5g1aBavvfOaNN4egYBt3G1KcXExN998M+Xl5SxZsqRef3+ANPRVwBvAZMANwABzb/7lfZww\nwmp2CXBjpopQNpPMozhxsIPOlBLOiyMyeXtBpJwERMCpSbewbRtMObCOzgf89aGmXTh8+9gfFYql\nuBMWqwvjz69C6jKYdg98P5SfiOAxepKCk7vb/co5O4e17MGc5CTwH0JaWhoLFy4kJSUF9xY3WSVZ\nQBzwKtAeuB7YBEBqbC8m7X2pye3EsoQ+TKl9/j7pzEufg1RFikBT/8asx6zfMMztAssW7MpDlbeX\nXpC8CR54ArZ3gWfuxl0awRt04gviGUcOZwY7GbJpCMFJwS12HKcCCfyHUNP4GxwczPxP5+Msrl8F\ndCvwmH96FvARYjPTw92HWF8CBRRgw0YkEZTzP8Ko5l0gBBNn2Nbxweb+h0z4Vj8JlVQPiZNVU9/j\ntDRY+KWbjbZ2xHoq2e0eT6F5OCavDcPigRtmw+gv4MW7mP71OrYQTB6P0o8wwhnPXjYQNyKOOQvm\nSLXOMWrOwN8S1PE2atQoha6/Ueemnlv7GLooWOafutTOj7AE1z0movZxOigFKpdElZ7+W++nlH/V\n31xPiNasqe/x0s1rFLcPVE6LWf0aG6sufvRRtdj4Wi3pPlMtebWLWjJ1uFoS/YFazFIVz7MKChWM\nVaNIVQMYUPdbkh/GMfN/ls3ipOrVc7jqZ/Gb+9Fc2rRp41+yHTgbneN/NfAXgrFjclTp1+EgCV2s\n70Ey0+mIAjqQzz1zBxJpFPOmLYOqYal12aiou109JUXfDybEyajme2wYMG9+OU6LjbY93mHNzLGs\nUZ/zU9n7jJt0GuZbX4Un7oM5V3PjhDzSih7nAsIoYygwmBD7B1TFVGLGDMDg0wZL751W5vcuG14H\nLgYKgX7+edHA+0ASsAO4Gp07GeAB9EjpXuBvwFdNbNN/8jp+Gvf3z83NpVOnTo3W6o4+PB+6Gmhr\n7RIHDnrRi0lMItK6h3D3PnozFQvlQN2H9n37y7in20dYrTpHyezZUs0jWr/6VTpbt8KOHfCCO4Nk\ncqjAwYqwv3B5aQ5uEnARBvjPCH1/hn9Ohx2dGD35B8pxAfeg75h/CHiFiNAISsp0A/Af2v2BqKFR\nvDr7VanmaQYnso7/TKAMeIu6wP8ksM///31AFHA/0BvIBIagW1IXA8noyFrfcQ/8TYmNjWXfvn2N\n5pqAvwIT0N0/Z6DPWVoqqUxiEgBRpv8xwDepwas/4VIu42MA0tORxl9xUqjfYGuxgMcDS0glFT3z\nK/tzBFX1q3tBSBncNgvOWKWHRfzfKP7APSheAXYDdwC/AmCz2aiuriYlJUXSJzezE9242wmYT13g\nzwZGAXuAeGAp0BNd2vcBT/jX+wLdh/LbRttrkcCfm5vLGWecQWFhYaMun6APcSbQBn3379raJc/x\nFH0ZyAG7jzB3NdHejdjYRxEpuAmmAichvE1wjIkeazKJSIqU1l7RamVkwAcf6GRqMTGwfz9gdrHY\nN48gFYaBm2qLCbMnAoUPI3UZ3PkifHs605/z8JXnO9w8CNwI/BN4u3bbQUFBfP/990ydOlUSrB0H\nLR34nehSfs3ri/zPn0cH+Xf8y14FFgLzGm2vRQJ/jZoTQFVVFUVFRY2W3oi+kMlED/xSXrvkcR7n\ndE73P1M0/OhcnM4NlOPhO04nynSA4T5/Kgn/pYCcC8TxdtB37N6Dv3T1S/sAdP2KjftHs6f4S2q+\n0wqF0W4X09v/g1Wb91N0oOai/ULgJfRoeOOAvbWbsRgWtv6yVVIoH0fNGfgtx/j632tpbnLZ5MmT\nax+npqaSmpp6jLtx+JKSkti5c2eDPv8//PCDf+lbwGfAdGADcDf4q3Lu536WRH8ARTEc/NkHsY7n\nOYNrSGMhu3zxAGyNSqGbv1ErJ6fuB5eRIdVCovkd9B0rPPhLV9OAS/QWuGA8xG4k/JVQCvB/q4Oq\nMcZkwmUfk3+bwx/02wHPAIPR1ToNm+5MmNj1yy5ik2KP/0EGkKVLl7J06dIWe/9OwM/1nmejq3gA\nEvzPQdfz319vvS+gtohcX8v1h6rH6XSq9PR05XQ6VVRUVL0unzXTKAUbFMxX0LmJ5agQQtQgBqn5\nfKqW8LVawn9VEb3URsarH3hGreUJVUg79RXnqlirU4FSERFK7djR0kcvTkUXXaS7YqakKOV0Hjzj\nttuUei94pVrMZ2oZH6p8U0d177l/V28nLFJf81+1ZMS/1JJ3EtSSSaNUB1OCAouCuxXsVTBFgb3J\n38Ga5Wta+tADAi3cnfNTYKz/8VhqisR6/rXonKyd0d1mvjvWHTxe6id8W7duHYmJiXz++efYbDYi\nIiKAZcBA9GXtd8AkdA6gOuWUs5a1zOApdHnJRBbPU0UipQyghKEU8mfOYzHPu/UoXyUlcM89J/JI\nRaDIzNQ1i7XjCtWb4QpzsKT8eUK8Tsw48BHF9+Gv8PWVPWkbtAPjyfvgz6/CjPHwyGT2+vqj27rS\ngJHAw0DVQe+ZlpbGoJGDTuhximP3e/VF76IbctugG3MfBj5Bj3zSkYO7cz6I7s7pQY+M/mUT2/Sf\nvFqvnj17sm3btnqNwIno6p+h6AatD2vXDcLKm7xFPPG6bpRfuZl72EsFVmAO1xNJH75hKCG8RTc2\n4LY66PFDJh37S0W/aH4ZGfD22+ByQUSkwtLvA6buu5Xh6zPYx8VQM4hReAnG2DeZnv0Zny1y+1/d\nAd3ONRxdj//hod6GXr168c0330gj7gkiKRuOs8jISEpqkvM3cDY63cM+dP3/Tw2WRqAHey+uPQ9C\nLDHM4QP/Y50DKJtx5HEaa+nHkhG9eW+BlXvvlcZfcWQaN+beey9EvZ5NgreSKkxMTSrkKfsY+m2e\nTjWJKEwYGGB1MT3xbj77ZVO9rTnQhZq/Ai+gg3/lb77/ZZddxkcffXTcjk80FLD5+E8Ua01y/oMs\noW6ox6/QmT/ja5eWUNIg6Jsw8RwvABDKJpKZAUAlHbCSyOk4OX1lTu0PWMZ+F0ei8XcmJwcSvJUM\npIRhOBlvzaLd/g6UWRMAM4ah4Oyv4Y2byFf5/q0Y6NTl2ehbcQYDj3CooB8WpjPapqSkMHv27CbX\nEa2flPibkJuby8iRI3n33Xe59ppr6eTpxMrClY3WikDfrXgL+irgKep3/wTdBXQoQ4G6jnJhrKeK\neNzEYaaMFG7FxH5SWMt6+hMWBhdfDLt3S+lfHKx+Kd/thsWLYSb+u25tLjZF3czggu7YjRyKQuJQ\n3nAclWCkrGH0hvsor6z/2zsbfeOiG333beNbbuqEhoZy1lln8dJLL3HPPfdIP/0WIFU9J1hqairL\n/N3i+tKX9ayvt7QzerjHUeiS0uvoJg7N5P9Xw46dZLrThlj2sAcbQUzibtqxiwSeYSgr+ZWk2h+z\n2+pg6Bb/jWEioDSuyim4N5tNr+7HrHyYyScCHz4iUFgJZSulJAMmFCYU8BQzWGX+HwcoxWaD8oqa\nLfcDHgd6oO+7nHtY+yPDJLYsqeo5wWqSvvWO6M1jPMb3EZmcWfvR/QKMAf6IbufeAFxDzd/Hhw9P\nvX9llLGWdaxmNVlk8R3fM41XKCaFfdzKckYCkEwOqSzjPPdCNoyUup9A1LgqpzKnkijlJhwvISTg\nIQ4foShslNIHsAJmDAxMyTnkd/4fRd5SPN6aoN8FfaftV+je1r04VNA3DIMRI0bUPh84cKAkWjuF\nSIn/MNQkfXtx2osU3lNI8rS2XDqsHwsLCppY+xx0zn8ruhPU/Ca3WXPlYMZMMMEkk8wjPEwooZgp\nxYcDE1WYjN30/PFKYvrrDKMZGTB/PlRXw+DBMHeuVAWdSuqPeFVWphPAvmnLoJ8th+LSGzBUVwC8\n1m1Y3IkobP7eZLqnzljrn9jl3YO3QYasRHQ+qqvQN2I9g07B1VDfPn35z8v/4bzzzmP16tV07NiR\nm2++GaUUb7zxhlTttDCp6mkFiouLufrqq1m8eDFNH8+lwBSgGn0PwMIGS2OIxoOXEup6D9VPCldf\nbHosfeb00eukNrzlXpLDnVoOSqkALDenMtK7DDchLHY8yA9tIngpNobn84JoUxitV+r7M1z/NqMf\n+Y7y2nbZ9ug8itcBs9D1+Q3TlNhtdj455xNmBc2SLJqtnAT+VqQm9UPTDHQp62F0L4kpwIJDbmsA\nA5jKvwitHQsYTFSgcBHEHqqBCr5gDLN5kns5zZHDwDMcWOdKC/Cp4uGQbKIrdHfMYoKINxcz1Psr\ndgqpIo69lhjCPWb/2gpj6Pfc+MvD7C5y4alNLNsBHfD/hG5zmobOrK4ZhoFSCofDwcaNGyW/zklC\nAn8rUlxcTNeuXZtI+FafAVyOTvxmAI8CH3BwxmrtTM5kAxtw46I7yTzCI4QSCkAMy+nNwxwggij/\n1cJaJrDP6E00BXQxzSBsZC+sH0sdUKtTr7X277GZZOVG4nDAR7EZ2Oa+jbvSxWrex0OM/wVewNxg\nE9OZziIW4cLVxBt0R2dNuQydI3E69ROpAcTFxbFo0SIuvvhiVqxYIUH/JCKNu61IZGQkp5+uUxLp\nVA9NUeg7IE9D39z8d3S/6Qwap4EAWM5yiiiilDJ/SogZDZaboTboAxjEEKQSKFOnsdt7K9Zli+Vm\ngNaoXmvtZQszahtut36eA5WVWPGiqH8PScOgT9gB8mPXNBH0U9CNtCvRefG7oUv8DYN+SkoKmzdv\npn///uTl5UnQD2AS+JtBZmYm6enpZGVlkZ6eXtsbwmjy5Pw5MAI96tdodK+gCeisGE1bwxoKKECh\n2M8IlvI1o7mYgfRnGIPZjc6KWHOTWBb96bFsFtPis/mv7X2yol7Bfd7ltUNFipbx8zbdOywnIoVn\ne89iIz3xYtB7X12lvonq2seqJidX522MDjqPc8ouJWtvTYcCAz043tfozOcr0F2LH0FnTm/o7LPP\nloFRRC2p6jkOanoBTZgwgZTBKbg97t9Yuzc6/cNV6Oqf56DBfQJ1etGLfeyjLW3JJhuP/36BEYxg\nKlMIZTNT+ICd5BKMmX9zC8qfIDWWJSSan8dx9rDaNoH6/cRjYyE3V/cmSUqC8HC5eexY1O+dk5Sk\nP9vS/GJmksHtzKIkdidVe/thQ5HNOAqNETypXiePPHazm7bE4oiuJrLLfr76of4JOwQ9bsTf0T1z\nZqBL+40HFwKz2YzXqyv+pQ/+ya85q3pawglNZdrSLrroIgWo0NDQJlPa1k2xCiYq2KngawWXKzD/\nzmv0NJKRaglL1BKWqAEMqJ1/NsPVEpao73lJuQjRKXpBqfR0pZRSo0bVzWrTpu5xo9XEUaj/2dZO\nEbmKEU8o7higGNdOOW129c4556g3O39y0N+uZoqwOPyPuyt4SsE+BR8oOOs3vxPt2rWr/e6lpKQo\np9PZ0h+JOEY0Y1rmltDSn98JVZP3/9xzzz2sIA5WBdcoWK4gT8EkBe0Pub4JkzIwGowRAKhggtUg\nBqkFvK0qiFM+ULeBOo1QZeFcBXp8gHFsUs+Z1qpZUT+qEFwqIqJRTvdAddttOnpfdNFvfhA1q4WE\nKBUerpTNplRsrFKGodRMblNl2JQHlA+UB1S12awSbTZlCQpSwdZwNT1xlupgJCqLYWri72vxFwC+\nUlCg4DEFHQ9az469QSHB4XCoHTt2NBhzQpz8kMB/8nE6neqyyy5Tl156qVq+fLkyDOOQwbxu6qfg\nRQX7FXyq4I/+YPDbr7PUWyeVVLWeh5UCNZL675muQKmnWVt7tfBK4nq1Y4cu6Qd8rKhfZP+NS5+D\nS/Y+Rfw6xahH1JIEfZXlMpvVF0OGqD+PH6/afPSRMtnrBjSJDbOpEEfjv2F3BY8q2KVgqYIxCoKa\n/Fu/G/SuOpB1QCml1I4dO1RiYqLaISP9nJJoxsAvdfwt5Lf7/zfmANLRA8F3RQ9r/AaHagsIJphK\nKulGN7rQhQIKsBOEBy9rWUcPevAUEzmbceTwD5wMI5gtjI8+hy2WLvxaGEQQbnwYpBr/YyUjCQuD\noUNPnjuFj3mM47Q03eUmJQUWLSLj3siD6ux37QKfD2aabuI6lUmwcpPDOCrpQDntqQirwKhqDz4z\nFi/MYBr5jhw2VG7Hc9BvIBzdznMzulfO2+gumZsPuYvLly9n5MiRR3hg4mQl/fhPAcXFxdx88824\nXC7cbjdZWVk4nU7c7t9qCAbdV3ssuoFvPzpAvAfsbHJtEyZ8/vsFRjKSrWylDW1w4GAiE/z3B7gI\nYyvVtMGHBYWVMHLow2S2cTsVdMCEi2hmEsMugnBzgEgqlq/hX28l1Q76ERkJa9bowHiiNQ70l11W\ndwds587QsePhnQRqtpPzXTHPVGaQwSzuDCqgm62SnqUrUUTiJRKFFRPl/lw5Furyr8JYxlJEERXo\nrGgmDMyPrFpUAAAgAElEQVRmhcvb+HsfBJwPXI8eyPxr9An9c5pqrK1Pgn7gkcB/iqq5CnCYHVR4\nK35nbQM4E51L/XJ0crj30fcLHJxDqIe5K9O9zzCBCWSRBRw6RUSNWJbgIooSBtY+78OU2uW7TImM\nOTOvQYqBxETIyzuMg21m9VMdpKfrPDc1BXabDVaurFv2W51bmkqZ8LTxPQNVTcptN3Co8Rq00VxM\nOYf6+1mAP6AT+l2G/rtlosd4OLgbZmMpKSnSLTNASeA/RdV0A91fuJ+vl319BK+0okuO16L7dq9H\n9+3+BD065u8LI4yZzCSBBP8c/TeazgzyycdGEE/xByAFRRAGVQSzBzMVtOdZ1tGHnbSjM7l0NbaR\nq5KoMq5ikxpBlTmI5yN6c8b51kOOM3CsXUvT0qDXwmy6WCqxhJqYbu/N9kIrhqGrY5QCqxUiIsDr\nheBgfSXwUG4GF7jm4yssxOy/MvIC64I7MKAyHxOK1byHm7aA158QzYJCMaP2s7Gx3viJCnXwmLSa\nHTgPfYL+I5CD7ro7B8g/xGsa+vzzz5k9e7bkwQ9gEvhPcfXvAzj77LNr00HYsVPVxIDXDQUB5wJX\nAJegS/+fAp+hB41vOk3EodRkfQSIIoq3eKs2fUQNE+UYKLw40CcML4O5nZ+Yhtt/Y1oMy+nLw/zI\nQ5TQlXD205fJWClHAT5MuLHyIwNxEs0YMilBB7iasQkqcfAn/3y7HYKCdJbSgQMhLAyuXL6OntX6\njuZScgljNoP5Hh8WQigjh3EUMQw3wVRSRD8eopyrqKQDJqrozVS2cTv7GY6bcBQGPsPAUCae4qna\nID+Bhwhrf4Ar995OketQJfsE9EDlF6NL+GvQJ+IPOZxgb8aMF90HPy0tjc8+++zw/2jilCSBP0Bl\npWXhXOjkwfAHWXVg1WG8wgQMQ98hfDE6GC0GvgQWAbuOeB8aVw+FkE0lifganQwMqjDhwks4oAN/\nPx5mHU8fsuqovvdJ51p0ncwSUkll2UHzG5thzWKQ28kmQunFA6SyuMHy+u9d8/5VRFPKAAAKk37B\nZUSTuKNR6g2zh7utd5JVtUV/BsODmHR3BJfc6KSsuqYu3gacAVzgnzqi894vQGdm/a1cTgeLi4uj\nsLBQqnZELQn8Acpd7CYnI4e4aXHceuetfPb5Z7UNt4cnEV0ldAF63IBCdIPiUmA5sOd3t2DDhg8f\nQQTRgx5MZjKhhDKDGeSRhx07E5hAGKG11US72EVb4vwNyhP92UcVUMp0/uMvSdt5gqsJoiMKKwoL\nPnLpxyQicGICshnHHs5DYcFCGYO5Hbs/6+RmxlFOJ6poxyDuJJfra0vy7UKfJce4A1fpWbVpNHz+\nHpGG/19N1U0ev7Ke9ai6hAmHEITOkXMWukQ/DNiIPql+CawGvId8tYEBBk2m9O7bty8LFiyQIQ5F\nAxL4BaDHBh42bBgul+t3soM2xQAGoE8Ao4CR6KReK4FV/mkjh1M1VL86qEYYYZgwNRhvAPQVwza2\nUUQRFiy0ox2b2NTkspd5mXja+k8gediwM5GJXMu1VFKJQtGNbjhx0pa2hBDCRCYSQggAM5jBKlbh\nxIlCYcZMOOG0pz1b2IIHN158mP1H+dvfylhgKLpUPwI9KHkOsAx94lwGjY71UGqqbur37KphtVpl\n0BPRJEnZIA5Sc3s+oMzmw0v10HAyKeiv4HYFbyrYrOCAgmUKZii4XkGfw7qB7FBTD3qo+cyvvbsY\nUFash1wWS+xBqQxSSVUmmrrLtW75Er5WS0I/VQOCehzlvsYrOF/B/QreV7Bd6Tudv1Iw2b8s/Ki2\nfdppp8mdtOKo0Epu4HoA3QHZB/yMvvMkBN2nMAndneRqoHFKSP8xiOZUU3pUSjFlyhTOO+88CgsL\nf/+FvykSXZ0xCJ1S+jT0IB9b0D2HsoFN6JLvNjhkF8ajZ3C87lM30MfSHT3oeC90wrx+6LaRLGAd\nsBbdMJtz2HuS0jcFCmDN/jUopegS1oXkYcnYQ+zMnj1bSvPiqLSGqp5O6MrhXuixBd9H33XSB9gH\nPIlOCB6FHhmiPgn8J0hN76Bp06YxZMgQ9u7V+dmbqpo5fHZ0gOwL9ER/BbqjB/J2otNM7wBy0b1X\ndqIbkff4p6YGEGluBhANtEU3aLdDt290QH91O6PLJkXok9hm9AlsE7oMs/uw38mBo/ZGrcTERH7+\n+efawF7z+Us9vWgOrSHwR6MrgYcBpcBH6HzCz6MrjPcA8ejKz56NXiuBvwXU3ByWkpJCZGQkixfr\nHi9Wq/Uw7hY+HAZ6jNdO9ab26IAb75/i0OWE/eiTRIl/KkVfLVT4l7v8k+7oqbdt9k9B6B40NvQF\nZggQhr46iQRi/P+XoruyFqADeZ5/2lFvOvorlJ7mnkzzTmOqaSqrfaul94047lpD4Ac9fNQM9GCy\nX6JvIXWiS/k12y6q97yGBP4WUL/0CdRWCxUUFLB69eoTuCdh6OAcDUT4pxB0PqIQ6gJ7EPorZEIH\nfx+6l4wLfXKoBsr90wHqTiL70F+7Q/eoOVoWi4UlS5bw3HPP8eyEZ/nl4l/o/Fln/j7171KqF8dd\nawj8XYH56JwBJeiRIOahS/z1A30R+hden5o0qa4feGpqKqmpqUe5G+JYHepKIBAEm4Kp9FUecnnN\nQCY1g5ODDGgiTpylS5eydOnS2uePPPIItHCvnmvQqQNr3AC8iK4kjffPS0C3/jXWMk3iokn1c7Y7\nnU4VHR2tABUWFqYWz1msok3RR9lLqOWniIiIJudHRkaqtLQ0dc6ocw5aZrVa1ejRo1VsbKwaMWKE\nateunRo1apQCGdBEtCyOVz+HIzAA3a0jGH0GehO4i7pGXdCNuo838dqW/vzEb2hq8I6asQTS0tLU\npZdeqnbs2KHS0tJUfHy8SktLU1lZWScsmEdHRqsUS4oaEdRPRYU1HdhrJrvdrrKydqigoHgVFxev\nsrKyGhyb0+lUaWlpKi4uTsXHx9fmsa8J9EDt+jKgiWhptILAD3AvOrXgz+jAb0VX6yxG9337Cmiq\n0rOlPz9xHNQf5i8yMvKwgnhTVxIDBhw8/GDNdhsH3t96T8MwVFZW1jEfiwR70VrQSgL/0Wrpz08c\nB/VLxTXDTDY1znBycrICVNeuXVVWVpZqG9dWhZnDFKDmzJlTu53wcH2DlM1mU2lpaU0G4KbeMzk5\n2V/SP7qg33i7QrQWNGPgl5QNotnVv3/g7rvvxuVy1aYiAA6rb3tubi4jR45kxYoVJB3GyC7SZ16c\n6lpDr55jIYFfCCGOUHMGflNzbEQIIcTJQwK/EEIEGAn8QggRYCTwCyFEgJHAL4QQAUYCvxBCBBgJ\n/EIIEWAk8AshRICRwC+EEAFGAr8QQgQYCfxCCBFgJPALIUSAkcAvhBABRgK/EEIEGAn8QggRYCTw\nCyFEgJHAL4QQAUYCvxBCBBgJ/EIIEWAk8AshRICRwC+EEAHmWAJ/JPABsAnYCJwORAOLgBzgK/86\nQgghWpFjCfzPAp8DvYD+QDZwPzrwJwP/9T8XQgjRihhH+boIYB3QpdH8bGAUsAeIB5YCPRuto5RS\nR/m2QggRmAzDgKOP2Q0cbYm/M7AXmA2sBV4BQoC26KCP//+2x7qDQgghmpflGF43CPgL8D3wDAdX\n6yj/dJDJkyfXPk5NTSU1NfUod0MIIU5NS5cuZenSpcdl20d72RAPrEKX/AFGAg+gq37OBgqABGAJ\nUtUjhBDHrDVU9RQAeehGXIBzgQ3AfGCsf95Y4ONj2jshhBDN7ljOHgOAV4EgYBtwM2AG5gAdgR3A\n1UBxo9dJiV8IIY5Qc5b4m2UjR0gCvxBCHKHWUNUjhBDiJCWBXwghAowEfiGECDAS+IUQIsBI4BdC\niAAjgV8IIQKMBH4hhAgwEviFECLASOAXQogAI4FfCCECjAR+IYQIMBL4hRAiwEjgF0KIACOBXwgh\nAowEfiGECDAS+IUQIsBI4BdCiAAjgV8IIQKMBH4hhAgwEviFECLASOAXQogAI4FfCCECjAR+IYQI\nMMca+M3AOmC+/3k0sAjIAb4CIo9x+0IIIZrZsQb+vwMbAeV/fj868CcD//U/F0II0YocS+BPBNKA\nVwHDP++PwJv+x28Clx3D9oUQQhwHxxL4nwbuAXz15rUF9vgf7/E/F0II0YpYjvJ1o4FCdP1+6iHW\nUdRVATUwefLk2sepqamkph5qE0IIEZiWLl3K0qVLj8u2jd9fpUmPAjcAHsAOhAMfAkPQJ4ICIAFY\nAvRs9FqlVJPnAyGEEIdgGAYcfcxuuK1m2MYo4J/AJcCTwH7gCXTDbiQHN/BK4BdCiCPUnIG/ufrx\n10Tyx4Hz0N05/+B/LoQQohVplrPHEZISvxBCHKHWWOIXQghxkpDAL4QQAUYCvxBCBBgJ/EIIEWAk\n8AshRICRwC+EEAFGAr8QQgQYCfxCCBFgJPALIUSAkcAvhBAB5mjTMjevefNgyxamq5U8VrWYIlNV\nk6tF2iP58fYfSYpMOsE7KIQ42fR8oSeb928+aP7nYz7nou4XtcAetR6tI1fPihUwbx4bP32NpLxS\nstvAqkT4NhFWJ8LW6Lo9TQxLJG9c3gnfaSHEySXy8UhKqkuILYPTd8KwfD09PQwWZJ58+cJaW1rm\nI3VQ4P/fBT0J27GbPaqUMban8OwdiIpqi+mhHIwIg7KXBuLdFQo2LxG3zqHyP8G4iruAtQKuHAPB\nJbXbGtB2AEtvWkqkXcZ5F+JU1POFnuTsz0E1Mc6TwwXh3d+irTOG7d2qGXVdGd+VnsN+qw1fwjpM\nhb1xtN/O6fc+y9wbXzmp4sQpF/g/PbMnu6OSGZKdTY/t2zH54EdrN5yeRMb43qXUFIHPF6RXHlUI\nTiv8FAVAWPB2TGG/Uh5agSe94UmgKSZMnJl0Jh9f+/FJ9UcXIhBkzM/g1bWvNhnU6T4OYoaDKQhK\ncwj6aTLDtpsp+2I2VcW9CPdYmKx+oio0gpjIKZSaCimPjWB9cieuXP41u30uri5bTsmB/tD7fdIn\nz2NO+pwTf5BH6ZQL/EPansM2XwzVcflUJzhw5ObSY1sk0SqKPRQQxTa2DriN/AOjYfcgqIoCDExB\nbnw2H5TaAIjst5VOt3/AWQunsimqnHUJsC+k6Z1I751+Uv3RhQgEqW+ksix3WcOZ3cdB1HDalMbQ\nbQskrCsga3Ey1tIOhCsPm9rspKzjRui/HXruh8Rwuni6szsugb6//MKgDVv5fx8sYMCeHN4nnWvj\n72Pg/f9kye0fnVSFv1Mu8Cd3n0zx2ans7xGCityAWv0cVDhhXjuoWg8cYGi7PlRHeqnavYc2VVXM\nqa4m1qcwUHixksIPbKA3kR2ctPnbBnb1qKZn3lY2fngh5NkwrNWM6HsWW9vuYkckqN/oz2QxLKy5\nfQ392/Y/vp+EEAEiY34GOftzDg7q9Z3xKVhCMXsg+F89MO8IJriyjLcqLyesupx2nlKsVBBmq2B9\nly6s7d6dVT278/7QPkSE2+lemo2x9j12VbrxFcZz0dxLMVeHENVjNoNZS2IpDL3FxIh2F7Dg5syT\nKujDKRj409LSOGfr/+OqbS/jdAykqqw91ZioNs3n2xereSbuPJxKwaap4CmH5HFg6oBlp4ekJ1/H\nrpzMiI7gjNxcQqt0j6AM9DBgDuAdYAcD2M9NbGYEZQQxkefwGFvAbMMR9U8qEzqi7l0K9hL9Pt7y\nwz6gCFsEWXdkSW8jEVAy5mcwf/N8iiqLUCg8Pk/DKpqUNyC4HWACQ5e0Zs6YwQ0LF2H1eSghnHKz\ngwhvKeMpZQv695qJHrN1M+PYywie4DXyyMeOnchLIiho9ys/td1cF7280PPtf3BGSTc6JD3BwMpf\nGbQbwqsNNhr92H7vT1S0h6pweHAreLFTOaHyRH5UzeKUCvwqI4P3DYPSzZvp/eMOHCV3U8JAAGJZ\nwno2cr35bWz2Crp3X8etdz/EY2H3sCs6mnFMZ7hvJUHeYnJKYfIGKHdD+M4Q3Jk+Kqv0H/cqYC5w\nBcPYTiU2bFSzmyx076BziGIcw7jeeIPE0zZw/+Qb2BnalkfUJMqn94dcM7g3whXXHLINQXobiUDT\noFom5Q2wxYDPDYYVLA7AYJwxg2+m30zGoo8Z6OqEIgYzVfRmKlbKUeggNIREfiBfb+s+wAb4oOub\nCWz7424IA4LA4tO1ObnnvIHLEcM487P8POVqdu7uTt6BLrRp8ytwgJTTxzA3q4SgvWFU7yvljtvh\n8sth+V5YXnkeX9341Yn/wI7RKRX4K889l9dKbiJu+FKK93XDvLUdXb/viNfsguBt3Jn6FTnl/aBg\nIBQMgOJOdOn8DY8+eglxcSUY9Y7ghwPdech1Jy6jqu7qQAHdxoGlA+PMmXTYs4aqch8vTYW8cuhJ\nd57iYVIZTwYlbEEBDsYxmDtMmTi6FrJzS1/9BqP2YH7wBwyfi+Bv76DUKDzig4+wRTCk3RDmXj33\npLvUFKe2jPkZvP3T27i8LpRS+PA1WG4xLHi6/U03sAbF1C2Ynkz7/CLC7CXceM8kBkZ8i+GtpNOj\nioLlZzKYNYRQQRZP1xbq3NE/seaCb3nGtYTS8ANUxFeAGX1Lqf83nVAKr34Ksb727ItO4JsJUbQP\n2kOCsYu//HkF+bu6QMxWaJsF8VkQ/yPkjIbv7yIh4SOu+NNVXHP1UqrLRmMxHwDgxwOh3HR+3kn5\n2zulAj9pabyWv5497Yr4pYOPC4ZEsHHbOXQ1SkkosfPUgmzy9lVir45ggnqYEIIxJz7CDckdyTtw\nBhT2gz19MVyhJCev5LGpl0NpJ7I+zqD7+ncZm/8VXrcNzDDt5TNJ6foNAMtXmnl4gpeuIXHkVzgx\nKwvVuPDiBSAdeIAHWGm6mAW+CRTwIxYsVHV8iynPjSU8fD9PzniNxauuRh2wgc2L0c1J0D8+Jvm9\nOCryojAHVZB/0Rgqgg6+SpDGZdHaNNmwWr8k/1RfHN9FU+HUnSsANtKLBHazm9uooAPVxOPDgodI\nQDGNp8gnn90U4MZFOWUEmc24LwEjzk51ZCUxyo3145lUFIymojIam62S9t2zeeCBW+lxQwjDXT8A\nsM0cT54jhIrgKsZcGEdJ4mYoi4c9A6BgIJa9fUjwxZGXfSZnjPiIDxd0pW1kf3788TyKixez1x3J\n2cN/pE3oyVkle2oF/uJiMq5P4LRzzsIRtYtX357CmvwziQ/fwW0T78YIdTOVCZSrYJjRA2ObBUUF\njt3nEVxWRLyvPYXsJY44gnEwznQ/4SYnIbxGhTeUMep93uADIrHiwctYLmW3EUFUVCEJCb9QWNiR\nNnH5+IKrufD+f9EpfAdVpblMtT1Pua0jeNzwzh1gKQQ3RMwLJ7QqlGKKsWDhWyzATVTQARMuuvEY\nLmANKaQzlxIiag/VMDyYzF4Sum1hV3gcvntXgyVXX5kEt4egKP0DW3MHuJq+mhjRYQQLxiw4KUss\n4sTJmJ/B6+texzv4tbrAveYO6P+kfm7WVTGNfTNjLCGr7OxydmKYWk0oZZjx1VbJbGYclXRgEvNZ\nyTqcFGNg4LvLp6tjbPU2Nn8mbB4NRgx0L+POR+6ke2g282eMZ+2q8ykrjcTj9nfTxtCTrQRiN0Hs\nBkb9NIefvd/jxk0pZXXbvQSIAdzAPODOh+CMwZjsZvqqX7nv+l6M+KoTSUO7AuB2F5OTk0Fy8iys\n1pP3d3NKBf5bMiqZHRfD8wtCSS604KkKYUybyZTsPx2juB2PP9aLxYsfYfGS61Eec+0uDxj1JddO\nelifFHx2uOt+2LwWgHAjlEpVjRcvoYQymMEUUYQNGxN4iL9ZU5jszqMdug3gZ8J5mH5UW7zExuZR\nWNgRMOjZ51se/Pd17AxpwyPf5lDuqqB+29W4wf3oYKqgymVm6h35lBdWADoPRhvAB1Rjw4qNBKLY\nRDUeugLh6CascKzWarp2+5HwMCf3TbyR3aFxlKowpvIg5TNSIM8OQW646Ewwb222P4Lctt769Hyh\nJ9ud21FK4VGe2vlhQWGUu8oPqnppoPs4cHQAb1XDzglnLgaTWT9Wiqix3TiwK47/+O7kMvPHRFME\nQI5nPHu4EF3fooAqdBQ3+Z8rpl/yNKtiVuGMcqKClf6iGzAuGTo4oMreiamWaZSb/dVAV50BRfXP\nBPV+PI790GYTMyse5Lr9q7ArLz/aHmJU9VOUc3DDqxkIGfcg1amn4baH4zOb0fW8BlaXQYd8uCVy\nOiEx+TiKbPR5YAojNl9w+B/+SeCUCvypqbBshYtHzszgsfv/RNKePdz8oo2u2RGEVfnw2nJxO55j\nbNIYSopSYG8fKElEh1X9hTZZXNC1HF+EggkbIcQD7z0F7m8gqgRKvFANlKO7C7jBPs+OrcrGy7xM\nPPH42M9GcxydvL8Syv9RgYkxZOIJ8fHee1249tpfqKgIQykTNnspQTYXT/7nD/RMWA9AsS+YbGcl\nUzdCuQvIA95H/36a8EfgI8CLiVuw8gPhhOPhAf5JPudSYViZZuvIvir9Ixo5ah7hofv5Na8nVns1\neQ/sZV94WzAMxqlpdPBtpcpnYuqafH2COpw/xKST77b1U1lNioFDqrmByRIOhg5843iKDuRRhZ2p\nxkRuVy/TyZtHL9N6DEPRYzrsXjIYZ3UEnXy/0I3ttT/6mv89JhNXqmF8oX7Ag4dQQim/pBxvshdC\naBglDDB3HYctqAMVb16F49tgXAf0SnZ7OUFBVZSWRuN228DsYqblBpJNG2lTeRHFnI6PIK43bmOX\n2g9AxCURFMUU1ZbeLVUmutGe/Dem4YoJJ8jt5f0Hbid12y5iP1rE/Z5n6OjLJyZoF75KF3ucJbzw\nyjNkF3XGHZlL4r5knnzxFlwxl1Gd9BDr588n5ZUw+vziwjyoH73n9scaaW32v92J0BoCfwfgLSAO\nfRqfBTwHRKPDXRKwA7gaKG702gaBPy0NFi6Ev/EMl5s/YqK9iI2Wjbhc0CMmhja0YbfTha06nOu6\n/JWOJcEEu0txRD1HSWUJ+91t6FTmosITwxjPPBKSvwVvB/ZUdsQd5sP8wEZKEj36SJVX/2BQsHMT\n3DEdc3ku58b3ZXdJJXZXOJHeNhRQgB07E3mAEKOA69sncO2e9gx3VxOEjxzCmEwfyqm7AunadS0x\nMXt4aMIYTCFVVLqrKaiCcg9M3QTlxYATiABKgGowzTMRWRVOFZVUUA1APPG0pS02bLSlilxy2cgQ\nPMxBn7X0+4WGFmG3l+HxBOHxWLHZKkhIyMXlMJET2Rl+iAaXATYPOH4CWxHcuxJiYw8uFR6FIFMQ\nEfYISqpKcPvcmAwTXuWtXR5tj6a4qrhBKTXKHsW6/8/eeUfHVV5r/3fK9KJRsaply5YtuYlmmV5v\nLiQBnBCCCbEBk0uwHDshIFrAEpaxQjGWiUOxZQIhAZxQHAjp1yQ4BWzcKK5yt7pktZnRjGbmtO+P\nc2ZGkg33fmvd71vAuttLS9a08553ztnv3s9+9vNWffC5oL0mC52KpuB3+rmi9Ape2f0KBgYCAleV\nXUU4HuZw/2F6Ij3EtPQKbxNtqLqKKIgEnAHiapxBZfBTjobp1PMuB0EGdRBuvQz63GAzYM02eGkC\nQrMLu6hzT/Ui/vnS9XS2lvDA0QYmKEc5U/8YF1EEQEdER+QwPyTGWBASuAqe467eDFrjMZy4qHHf\ngTcvyhXfWoBq10d4AlGHojC07j+GsTUfEiKUhZkR2EnTPy8wnXryDaIGGcdpVL/P1UPvkquG2S/e\nQy/n87j+LC200kEHPfSMoHr6v+KgXJKZIkYpHzDYdHk1Z8/6iAK5E0V1s2FvLWOO+3hiz1W0/2Ma\nITXAoxwibrTgweAJm5MO591IejFRj8g71+zh48l5nJjgJeDpYWZ4Nxcf7cWX0YpU2AZLHoZdpzFm\nzhimvzr9f+5C+f9onwXHn2/9fAh4gR3ANcB3gB5gBSYpKxP40aj3jnD8AwNQVvYkoYFJjLn0p6w8\n/g77c2C/WMj+rd/gQCJOlCNAE+ZaYjoSEYHG075Haf/Z0F4ImoyBhoCBSJwM9jOVOiQiGIgMEGAj\nV1BIOxM5QjPjCOFnLi8Tkrw4nEMUTdlPx9I+ol4nDBwAZ8Ckpiky7G+CZcsgMtJZ2rDhxo2GRjnl\n1HEnT1PMn7mIN/gZAWRchEjwFmHszGU9QTIoLf2QrOxOamrm4vUG2a2W8IBWSyQeAUcuxLtN7H9f\nPSgRrDsamkF8RUSPmfMgIVoIrIGExIvcw1QymEY9i4jQhBm0rcHHAc7mel4jSADQEUUNENB1M3MK\nBLooKjqMyxMmEOjmb2/PRdeHCbhODlGRuYMxgTZ2br+cgUQmTB5g0YMLmOw/TExVqVfvIGIfa86b\nGgFRNheaeFf6fKwFp3qyBRFo1uKoQfX5l1As9xIzbNTHqojYi6gWnzYzGsFD/eC3iThKzPkZPk/K\nADjzT72ofRIMMtr+G6+rvuirFAttxHBQ/1GQyOQHTK66IJkhkCBAQzm0uChxHibWlkGofwySTWHt\nmkrWv7yElpYyjneU4x3Ty+rue5mkHmHawEEk69rWEZBOIVnwLSo4gED77A7ysnPpyuwmKAVJiAno\nwAwqAqTF1q27OwXF+GcyIGRTyXb8hPjG7BNEI34+yQ1cRya/J0gc89S8OJjOWNoYpJUu80VJvD0T\nM6C5/i2wucHQqNYf5Tz5I/yEkFFTR5k1Hxy9oEXtfMAZ9JNl3RcBkOLsNaZSrh9DAN7nBep5hRba\nOMIRBi2s/xLPTJaefSVCcQuMa8YobobiVhqEalqVcQSCIZa8+TG2QzPhWAkEA3jPdHH6387634j/\nf+pDgDeBp6yfS4AuzIVhEzBl1GtH8vgNg3+v9rL1L99HHrgAnywz7/q55BblkJWRyW8OP8jZZzXi\n97cQR6afDAq04+j9zSwJ38FgbCpEQ1BfA1VVUFwMsTg/faWfif2n0WuTOT7rbSacOMrYv9+KoHiR\nUCiDlw8AACAASURBVJjJQlx0c5jxtFBCFHfKKZum4/f3Eh3yo+sCdleMorL9dNRZC4M5ejM6CzfB\nnmWmg24mBfH48TPEEHbsTKGch7iLfDpo4hxy0HiCxzlGG26gjvl43W+xwN9AW2gCNluCZ9bOIi+v\nGQEIGy629sQZY9cpcIIsgE+GmA4HtHHUxRcSUQXYXY/jGRW6DRIkRrgPCYnneJQKHMTIxc4JLuNu\nosSRgO1AhXlW7Odu3uVr9OKnnmlEGH6zaCRhNgBZjvPii2W89NISdrWcTbOzGGr3gdfEqV/gZrLp\nw0aCQ8YEftOwhOMtU8lw9jMms51t264gmvCSW3aIx+quJt/XwcqVjRxumQFOndrab7N+bS2trWWI\nDoXZtY/S4L2TCN7hV5XpdFeWweZsM9uZPABL3oNMD4i29OsMlWrhJxTTyvqGOnY0n4fuALIM6HKB\nQ4MHPoRny2HzGAjJ4NDwOEOMz9+P1xOitnYuixZtpre3AJstweV1awhvLKajdSKLP36ZMqagMQYX\nQc7iARyMjPiTRdIhCtCRARsqLkBm5ewGNpdtJuQK4VBtDMlxdGkUvj/szq2urKDYpVFAF6sbnmbX\n5stQFCdOxxCRUIZVQE1GDgKNVFHGgVQA1MAQf2Q3KhpePFTOPoNN2e+iK7pZPB0NV771Frjd5nwL\nQMMUaHFBhwtyY+DR+P4Di7hv3au0vF1JfyKHbsYwh9dxMYQAVGE2WB4GxgGGWIpXyKVd6EV0ONj9\n8JcxiopAi8P2R0GPgSRAh53x7wd56vZivD0l0FQOByez6qJKDv52CrETbnyKRp26jxd+EKMrV+bR\nJ7vJOWMSU16u+Nw6ffjsOf4S4O/ADEy3lznss/uG/Z20k4q7k4/dT55dJdTvpq7y1zz/y/lkRIsw\ndIX5JX7aB0+DvknIuspdd53B6ac3k5sLopguFynIHDImkCMO0E0e0fgQy7ceJqJiRkNRYAM4EjKK\nriEicu7sc/mw8ANigwmcGxyUxcp4iB/iElQ+sJewVJ5CJFI4aprSY3e5g7gzggyGsojLdlizA/Kt\n+EgJw+IFcGgkO0dCIoMMCinkGMdSEYyEhIhIOTn4CPAjFhJiJgF0FAQWMpNunICBz9fPunVnkZ/f\nPGJU+ykjy+inmzFE8NCf0CmIfUBM1alPLCLinGzeQMkehygmSHcxUIbpyzswu91ipqBdBRXUs5wC\njnAaS5CJjOiKFrmANoI4sfEGLeTRwwGqiTCOLVRyH2cTGbbtgyQl0LSkIwLTGQ3XzzDw+XoRRY1g\nMM96TLNeL476G2bMeJcfPzwbnzeIakgsvH0bR/acmf44UcPlijB58g7yc4/z9tvzUlmMKMYBCV23\nDftcc0Fr5DZO50ISZCOgMIH7KODgKW+YrbxAgmw03OnDEmHF7EZas1vZNXYXhmhgYOBUnDgUB2t+\ntYLbZt5B5PSh9GklqTNlVubhLQN0qoXVlAY3MSEDVj/RyN823oCiOsGjEnD3EO7JRtPsNLKAMg4Q\nxZMOYlx94G8DXxtbjrzHcn0nG3kbFRUPHsKzw2mGTLIGNhaYWg25l4OULs4WzS+gvzcfRbFTWvoR\nPn8/tTVz6bvuKi5P/A0HcXYyk3Ecw+BrDFGMSIxreJwWgmhoCAijBNhkYBws7oTsKAhwXymcXehD\nlkA2EoiawVMrnuJ421Sc7gg/WHgrwfsW0Hn2G/Q4Ozm7Moc/Nl3Nb474IfsguE/A71+H9dvAq/LN\njGxeP7PiFN/c58s+S47fi+n0l2NG/f2MdPR9mLj/cDOWLl2a+uPNNy/lrox/Mv7ff4c+oZ2l6/x0\nten45RM8eOl5eHomQFshRmsReq8fpWQV7oMLELQAOjYOyQZBUebZ3DhzbplHaWmQggJwOODRFc+z\ndd+XyBvYz1dOu572eAFH3/8aR4nQQxfc8mcosSKx3cDrIAsiF3vPpDc6hGy4uDOziuxILnZFQNIE\nfDQxnkfRBI1Zjr/SHDuD9J2rM750Lx2tE1FVG2Awfca7LFt+Hce8xdSyPB2lKlHo3g931EHvyfiv\nHz9TmEIttXjxspKVtNKCAye1LGEtjRynDQUPP2YR2WSho7NIOI17nr2S0tLdqSAYQDcEIoaHJsqo\nE+pMJhSGuRCED4Azz6STik7Q46AnYPswWmkYKK8GXzHYYrDfWjw0SOFJyXQ/Dp4NHpwxJzoKKhpl\nlFKAjSAfAOOoZwEZCEyj/qTFJNmyfypLOlkBhbNYiI0TSBg0cRdRitkiTWO7dhVHUDiCGU02W7/9\nwDIWYWMyotU9KhNBFwRiNhvf03X2qAXsWdyJ5BOxqU7WPruW/GA+mi3Iim8+QntGH53eXvroR5cs\nZovFbjnJykZBR2oElwKZMcgcgqwh+OdNfwVBQNTAFzIIDGhEassZHMhEVez8TFrAJLWJoBqwnHkm\njdxGPS/Ri4KIxNmSg16tBAcuAmSwnR0kUFJBxQhLUSHtULgU2+l+XAzySs792OwCGhJVxhoyHh+L\n2u7F7QjT3j6R/v48nojew9X8HgdxMukfccpf4wpa6MSBg4+FncSG1XuoXgLn5oIwBPuWI8WiFEUM\nbp8KRW7I9JkZrOpwIkgGkhinrw/a26G6OoZuj0LmYcg6wk8iT+LOVjhanMW6SVsZcPcy1gsOVzYH\nfjcfesuhdzIcuxQuPcE5jxzhPzMz8Xd2wiWXfMJV9dm0TZs2sWnTptTfy5Ytg8+A47cBvwf+BPzE\nemw/cCnQCRQA7/BfQD1XXgmT+i5gb998mvUKnB4F27xWjkx0EAz4ELQYhJvQ9zxE0ZszUc7Yw6Az\niE2Ftf3XkX+iHNoLMTry0aJuBMdBpuk1SFqMWUVP0Rw5j+o+lWIjhiIavFbexwMPX8/PfraClpYy\nnM4o8++8i9WJmwkaAgXv19PmiXDCBXS6IarBhgRcbiBkQ6bgoeDNMXgiGbQlTtDPADZsrGUN+eQj\noODlIDJhohQRJx8TiR/EyyHyWY2DMDs4axjeDuv4LlPEfQzqPr7rXkvelIPMq7sfxStTRBuZ9ONk\niBgOFBwsZA3d5FtQ037YtgzWREzHO8xERBw4GBpGkfPiJU4cO3ZcOMgnF40AeQj004sTNzXchpsi\n9lNACAf1TCMq6LjdEcqn7OCHS79LobcVCR0ViXc5jwAhCuigm1wiuKmnll/xbdxEETCI4URC4xCl\nhPGb9YBwK9WJByh2qmmsf6LpMKu9v6NYOmEyVqglYjghuMvE9Z35VDvXUxzfYWY0+8xCOmVp6d5q\nqZHzhC34CRHDQQwnXcYYIsHDFHndZMpxbIbCITWXsJBBvXYPETEAktcssgoCYHDXiWsp9nsZshew\n+pb1dLVOwDAkTrZkyG7+v4IbiHKYIEfJR8PHePzk0UkX7XSQxxj20jQiAjasWs35nMu7bEafraey\nMaFDwMi/E84+D+x2cNlAThZaDarVhymWOigQupFRea5hBa0tZdR1PMpXlb8gqQa6A2xhMBIC/WRS\nyXb+zFcppykFBh2kml7OQ8eOSBwXnYSZxOM8RQutNNFE3CIjgASzHZAtgeKFzCHkwABFYRg/AIOt\nrxE7kc+JjnHk5jbj8YTIyOjm7bfnWXMogJSAjONUK82cN2SQqcoYgsGgM0H7+A7Gt+bjdIYR8ju5\n5cgyehIhJMPJs9I9XDZjBR3LJ9L8yt2s2/8ypZ7NlEbCTOqD0n4oCQqc8Im0FfpY9qNzWX/drz63\nPTCfhYhfAH4B9AJ3Dnt8hfXYY5hF3QD/jeJuef3F9Gy5DF3Ngv6JOIaO4Nck9PgMYrFSYrFsNM2M\nUN0TgziX7KZ/XALQzJtGspuMneAu+OhBaI2AHyRZQHcYCIqAt91L42vryIvlWSegoXKU+ZnTaY/M\ngIQH88Y1ectud4gpU7ZRV3c9Hm8QBYmISXQjhoPD/XFq96pE4kCCERAJkIqsPIoH5wYnBbEC3Lip\nZQlevNYUasgMIhFGJRPBqj046WYAP9uYxfW8xgrutVJ4FzeJv+Qx7qPMOMhEjiJN7icUcNFSEwav\nhobIISaZjpWaUTi4ifGiabBrFzz44EnF6k8zEXEES8eHj1nM5EN2EiKKAwdTKONhvk0OEcKci4Hd\nKj0biAispMES3DK4i2IW8FJq8TPtVLG/werVF1JRsTn1KtWQkNBYtaqRjRvnkUi4P+G9Iy3ZaZrA\nziy20cxIhpEJl5TjwEOcyRhIGAhEEZEwsBPHQZTTWcSNlHGcPhy4qKWGRSymnXZ0NCSr7ekk7v1i\nzEanJLokAKXVMParFuMMRi4gUM1KimmlgA5kFHzCID9Z+QzNrVOwOWL8o+3LOPsNDBtMDh5iCY/y\nXX42moUJwH6quZePaKaTXexLj6+62qqPxaC+3rwuFmM2pAzPXPbWc88ZE/nXmgfobJuM4NC4bUk1\nZ/h28OSq1byd+i7M89jLNHO+BYlZxU/QnBfn3d2tPJx4hxa9H4fhppYavs23iTKUXgRny5BtB8WD\n8PrXcT08kUT5ZDSHH7sSJyMYwx3NJL9LIr/DoHLvAELUT1afxKInFS56o4oe4RiHMyE+DNb/PHfM\nfxYc/4XAP4CPSYPe9wNbgVcxM+tj/DfonABXvnwlbz91BUsmbCNrKMr4zk4mdHZS0tmJOxZjhz9A\nj5zN3KI7CA5Og/4SCI0FfXihxiBQPEDB93bRW6LRM0bAJKsYyadNlka4CfbWIaoRxJjApJcK8fcH\nqJHvwxsqQjhpSgwOCRqHyONiYx9ugkT5I3N5nuIZeygsPMz27VeQSDiYXLaD2rpv4fWG2WVM5cHE\nvUTshaSKwH3bIPPM1EYSbKszbzAPoALPgivoQkNBxsYQMZw4sWPjZRYxEYMY+QxRgJNuJCIpsSvN\nGncIH4pNxFAkZDRiOOn2ZFEw5SA76+wMep10UWBG5EYNVTRaPHCHGVUPXyjAhHz6DkH21DRuNDgI\n3/0udP/XWkUCAs/yLKWUcjM300Zbytkk6xqTKMWHn1pqmM8t9FlNRaMXGgkJF04cOCkgny66ySPX\nWlDvpYblfMTHqddnkkER+XTRi4rKgJUO+fAykzP4iD0kSODAQYgQCsqIY01nOgUU0Ekn+9mPhESU\nUT0So1kts18A31gQBbNb1gAkm8lyYZXFCHJSH19ERM6nWnqKy9mILKiAzsqVP+Odjd9GVW0YhoDD\nGSUec7PWWMjV/AEHcYZwkUUfO3jVqitozKSKNq5liGLqWc+/2E4/Ago6IoLF/EqfHYuBbC3tAc54\ngurA2xTTyq8f/RHvb/06DkVh2uT3WbZsDo2NK9i8+Wr6+3MxDJFGqria31sNYAZRwcGgKNEtLqZO\n2cZ7bKef/pFzlOy0jQHkwH19cNod4CpGdBSgO7JTzWbiY+X4tkpoio3S0h04nV3s2fNvDA7mm5/p\nHIDMoxA4hmAPMabn60y6cQ/dpRqHS8EwYoj/ugr9LR2pX0KTNM74/hm8U/XO/0b8/1Mf8n9pJzn+\niy5dw7/+ruOmCeGrvyUyrR0kNR1FXw5itkCFUsG9by/CU/Iahb3HKOnqJjscRsD0mwOih8qiZ2ge\nmgXBYlBG78Ji8EjDl5lx2hYOGmXUSnVE8Jj0u1Y3JY5DLFiyGMGnmzDEkrtge3fqwnUrLqZsmMKy\n2DILd2+glWPkcoBncPGi/DsyVQ8xZHJoZBJ7GBKcVHmfpHWwFMMAUdSoqHiX5fXfQPeoNFNCLt1k\n04NkxTu71Qk8EL6JyN4GE0cHE0tXgPb0nIwuyP4s9ixv8OawWkCNlV2YJhDibBbwMSsYohismfOz\niwoeRCaCimh1GyvIaKbfssYVx86HFv2ulMNM4hCi9VwfmcxkJ82MZ+rUzfyp+FL8nQqaE/bUQNwr\n8y7nkkmIGM4UfKVgYyFruJGXKaaFMg5gI4GIgYaIhI5gOS0zczBBs+3MJIdeCwZKZzbVNHAem7GT\noIkyllF30mJWTQOXsxE7CTREdjGDB6kf8bpqGtjf8HW2bL6KUCgLyZHgw8gs8ukkgY0bJq5j3rG3\nmKQfYQJHaWYc42hmwOcmazBMltGPw4JDkuMXMCmZh4lix8ujfJufsoGjhLBZC58XJw2spJVW3EAW\nTv7EB2h8iiXVLK3CcHX2exSLXcTCx2n4SGfc3xsJ9VZwomsCQSWAOCjickQoK9tOzpg2/r5pDk8m\n7rCySjfdjOEKNpJLF3OZwWGiOJGps1+DTTsPQ/MBAldwBepsNb3wuazrMQasmwoLD8LHh6DND06d\nRTWmZEPccNLYXcPrP1vMuO52VFnj9zc4WVzTyZPK3ZRxgPN5F9k663d945l9zuUEx3RBxnEIHAPB\ngP4JMDABumbAO1GQd8HVR2HGcWRDwzAC6M/+GuPE34GHuebaa3hjwxufNpOfafvCOf6Kinmcc852\nrv7KQbLtBp3BUt766HY2VUylO1dEsWlc9noVIf0YLX7o/s10jBnNkD0IigEbYLXvfirCF6LH7CAp\nCO5m8vxPMK43rdGvIjAgutnoPJvCeISoFmAuv7agBnMqMjI6eemlKVzzjS401Q4YVK2+jZcqrkk7\nBkMHQ7W421ZTmCCBqsKCBVBxLHUzOINOpsenU7dhKWtja2mlFQeOVNHWQONdsjiLg8OyiZ8TJIAk\nKZx//m/Zs+cCEgkHZWU7qKu7Hq/XjFwbGhppaSnD4YxSXXMrOd4u4thJGBILe26m215upufOQnBk\nmdz60CEYDMOd9fBgFeQWgxKDO+sReqMECGBgpAqDKmnpgGS/goRk9S+4CNOHjI2XWcwEBLq5FANT\nf0UiTCULcGJmBjpJ5y2xk+dSRdqZLOQ4N6ZwZR8HmEYdBgnAwD5sDJCmQorEmEo9R6giSjESMcby\nBJnW8ZLaMlHGoaFxM1/nShSKGeIy/sI06jlMVeqzksXeJNCSBI6OYHYsZgC/wgSQbrOe85AGlRYA\nv0Qi/gluekQGM5uTmVR3Y7Jpss8DWybVrOI83sPPADENYnIWEjrPNaxgy+arCIeSvAkT43c4ojid\nUUKhTBTFnWL5mLTN8YTw0S35eV07ysbZH6Fmq3gVLzM3zKQv1st+miwKsDEykzm9GrzFEI1D3jSQ\nPansr5qGVOfwWezg2tknzA53BLyePgYHs0mSHxq5zRrPIZrlHCrVvdhEnZttNvaJIqI2jgeoI8Pu\nRpGhOU9k8YzrUH0JUJywoRKkxyFWQZYWZRwtjKPZ+jnIePYzPudfjB8aIi5DyZ1gH/KQeEyksrKc\njRs3Egh8PqN9+AI6/ufPLiHnWAixdyo3Vl5CKDqBqUW1hFvvROmvxOuOU1U1l3HjgmRlgSRZA49D\n1OYkIdi5fWgF5b95kA/OvZFgVjEqw2iLoxmDzfDwn+7hnN4LERUXWE1ffpqYTl3q5geTJlrJNnZz\n+ilPxuPpJ5Fwo+tgc8axOWOomg1UkcKyJjrqek3evzp4ciPTxCpwFcOxGNTWE4jYGEvRMAjDk4ra\nzeyiGQf2VJ3Azy6msRSBOJVsZzdJytpIjDhtGqKkYegyoqjidEaYVL6TumVz8HhCyOhoiGiIRHDj\nJQoYxHDQFEywbGgBEdfkFNabanCa9cJIIbB4t5mCdWJGfxNI0/5bMD3lJ0hZjLbRcE/SzD0V4qm/\nffgAsGMnTJgEidRzLlzYsTOTM/mQj1NQ0ikticFb0NvoYjmQdooFwIyTmTsjbNYLVLteoVhoo4B2\n65vtptvIZjr7U+hZchFv75hIKJhDIuH6hAGay9F0jrCSGlyMRSLBT3iId2Yfp79sAEES8HZ4mfna\nTD6KfZQ+39mYjRoycPYL5ncm2MzehicqoNUNTt2UPfGq0FDO9JbDHGo6A0nSiUbNSN80I9UPcBH/\nRBzVcGYAIY+HjqwsOlqe45bZt9KS3TKC+UUEmGaNB+DCt0zxOENj0n9WMbn1GD+8CewWozSiOhiz\nMIcZR8xMsZlx9PvGMphZQof9r/z10AGOXwHHZ0CnFxx2F1tu2kL93fWsW7fuc+304Qvo+K9b9RAV\nm6byt7ZltHRkEBw4Sr5NxGP4eMDzA/y9JYiKE4EoXvZzx8wfs1O6jnjPlWj9FRAuRNAk7rqrkoqK\n3eTkgMu6b2KaA7uYwBDMhLspVsJvf9dJqxClwwdBQwBVhAEdCo20g2qG8/rPIRIYwqk4yIz46Qx0\nkq0d5rXXYmRajqtb8OEzYjhRUJH5F+dzLb9NFSz3MuUTi4myHMfjCaLpMmWTd3B33S38onE5bS2T\ncTojXF+zFM0rUE8NVfqTFA++Q8wxln4hhwJbxGTBdM8h8lYDuCOpm2pqfCorNjzG2lgjm9lMgoTV\nVbx0BPRjfR1k8h46XkRiTOIR4ogWDCIRx8EstrGEH6cKzPN4iaDgQxQFZsx4l+XLv4HXF0x+XMqJ\nOZ1Ramrm4vEGzT4f6/GOjokoip3BSAC7fQiltA9l2WEqntWg1c2BprMQJR3RpvLomn9nUv4eZDQG\n8VKlPM6N6hqKnQnWr1rOjuYL0OVBFt03n8nZncQE78ju3mRXb+4VlrZNA+epf8GOSlyDTgqI2Iuo\nFx5MZXQvMJ9XGmp5f/NVJBIOHI4h8vOP4fGEaGsrpb8/j2jURyMLKeMAiizx4KR72HrgS+i6bUQB\n+Qb8DNFCAOWkcvMCYBszrOLwj7iEu2jjWjq5isv4Uup1wmwBMVtEUzRLjdIqxPpjYBsAn9mxXB14\nh2Kx06zXGDVEVlWajVWWM3/BO4+a+f9JR/sEDERkKc4abTETtWNWA+Ov2My5TOQIMmpKngIMVGRs\nqKiSxIlAgKPK7TwSfZtmujigHUU1VLOOsBjzOpx68oIol1YjO4txxWNc9sd6ps6rZKbjY7xCGEkw\ns7pIBHp6oLGxke3b56HpdvB0sE5fwFXxLZzQFhGlmIStEFxBtMF8RHRwx2izhdhw4ZsMfnOQPT17\neO/W9z4X8iD/XfvCOf5fnyNT0KXT4jew3QRiHkQM+MVmuPhcyMqAIZ00ZQ+gXYYXDbPf37LVPMFp\nJ0XmBqYUWoh95DKAl1XuPLqjZpFIkhTWrq2ktHT3KR3W/d1V7HCfj+Yttj5Ohe0LYOy16Qv7w3qQ\n0lj8LetuYX73zaPGoTKTBXg59okToyDxLy7gWn5rMXmaqGQHmoW7z2ILA34/Hs8AHR2lnNwEZQIp\nCAYIInhVbl1bxYz8nRYWvoTIwDHY8wDVZ06k2BElpkvU93yFSNZFliy0YcJX/YcgHobf15uLyljr\nEDr8dN1qKrrTDTFmNnLquoIJsKkYqbAurW9khtUiK2mgn3fYSgYzOUEmOusB3ymkC8wzFEbkNCH8\n2FBGLFTNjE9BHcmI1AA+4HS6KGAFhxlAw04WD7AUFz4O4SWMzJkcw4WEiMJPeIhD6Dhw0kYbA/Qj\nD6PvrqSBzWwemUXMBq6qhmyLJfNxPegRcALnvgBuS+JBHTRpo6INbjkfeu2WNs92sxGwoTzVhTyt\nbAu1dTfwYmMdWzZfTSLhIKE4ECWdNcpiznBvpz1YkpI+EEWF3fppFNCJhzAixklROYAmivT7fOyL\nPMMK9WXrXPpHZj4l1VBQDMawrCaZ5Qk27lYeosge4s+r7qCjpRSnM0pbWymdnSWAgKYld1gxa0EZ\nBFPfXZNFHVXIYiWrUtdRgICpmSXaqCn6Dt4xQ1z+5fvQ7Gb2N1X3sSK3AG9uP2T2szJcx5/8F5iE\nDgEu8Pv5fUUFAdvnt1N3tH3hHP9bFxYx40Ccsb39oBu8a5/GdZ5nsDlDfOPrN1JUFCQ7G7KzITMT\nZBlC0QxszgSqIbJo6FF64gaqEUHbPwyCaPHwyInreOpvf2FgqJ+EruEWXGiGQTll1FFnOSkDiSAg\n4eMAY2lgJ1O4ntcJkt50AlRcrih2e5w1a2aRn9+MDgzhQsHOQp6h28gDbchMoSOHQA2fGgIYAg5D\nmW8yASWTmg01+GIjHeZMFuAbtlBYbn2E20y60SB+to7Q4kmaxkicC+yuCCXTPqah7krWrl2Rkka4\nvnZpSgohSSGMGQ7qjR8RCR6DPbUmPOUuBv90i34omF3KoSZw5YPD2lNg20LEYA8X/PJM3n/sLhLZ\nfiRFR1y4CFt3r9X78Az5FHIHd/IRH5mDS0aMApzefDrLX1lO4+WNbM7eTFgJY99gR42pKCjo6AgI\n+PClVFZNG43tmTZ/9nx6s3uxKTaIQDAQBAX8G7w44nby5GzckpMdsT1oo+GlEcwULzz2PcgdC0IM\nBkPgzoeMXJBkWH0WtHrBJUDtPqo9j6YVNKkhIvgQNXANgfpEGUabm0RTRuqbFUQFwxDBYs/M4yXc\nVh9GMjhoFouYo/8GtwV3GYLA+/yCHxsv00o7DhEyjGy6OMHuRfvQfXpacfnTGs4sho0gyPzQ+XPe\na7iFfe9dSCRsXlMORxRVtROPjyRONLKAebxMM9+jl/NR8Z10oNRdb0tAYAAy+yGrj5svl2gt84Aj\nAk3L0/tqu8zrrFpcTTEtKKrMn7Z+mY7AJWQO2Mnug+weKDsIvrCAPQHZvfDgQ7BvGswZM4ZXp48U\nZFuwYAEHDhzA7Xazfv36zxX884Vz/BsuOY/v3343nc89j+34MVyHuigmE4/mp0ZbigcPSsZxPPnP\nUBhsI7+vD0cicdLgdQR0RK4d18h/RF/j4ry3ETwae2tAtXxqFKeJDhs2EthYGFxM956nT3bMAC9i\nVvaALKszNln0dOKgjHKWUYcXt3V087xE4vhoYjL1xIHtVHL/xPuRc+IpqMBmS6QWj9TEALuYzhIe\nNpuVtJj5I0gg+yyqm2BG5IkBiHea495TD0fM5i0xIOJVvDRuWEt+LJ+VsxtozW7BoTip3XA/3pg/\ndTwb/bhoT9FCZUbOQQJ5RGFVRWCATCrZYUE/TcM0jgIpip+DOB/4Z3Cr8DPaI+NwOKIUljVxNLsI\n9R+FJrQmGVA6iLc3QU5uK9meLu6u+Q8yvT081PASHS2lHHeORQ/oFHX24HMOcE/Nd4h4nUTwEoUq\n+gAAIABJREFUUNvwKnktIXzOAa6vqaPBe4cJ1VgCaTh1CCSg05nGrT3KiKLkew23sH3jVaiqw5Ty\ndiWwD8IaZRHlNBHW/HQbeZTYjxKTXdR/vxrd6+JHv2lgXE8LEwabOK4XUxJsS7GgVGQieLARx46C\nzRInSxaKDUp4gDpcFFjfucxKVqUgOScO8snHjZs22mmb3TaSCvlD0gyeODBc/WKY2TS4/aJLGGvr\nJ4aTATK4lL9jJzECdtu9+3wMQzQXGgSqaeJLfEQO7XjppJ9Z6NjRUr0nIzO8GmpYPHsxvdm9KEoC\nd8xB5MuLUHMKqc74T4u+amegZx/52SXEBA8DBMinkxhOphu78QimquiKFY0cO1bG4cOVKIoTkNIy\nFIKdf9ke5EL1GG6CjBcfwaHG0e0yi7RxHBTjCHYHh7wtqOEB3LrEL79eTWnom0x5aRph3UZHB8yb\ndxe7dj0FJJgzZw6vvvr54fR/4Rx/1X1LyDhQzoeuPLrfmEjfg+30TdCJZKnIyiC6LCENHED/WwNa\nJAIb3LD2MTMFUBRYuBBvd4S1PEMBYxEQENAQhDAu2nEY4RGOTUUkiJ/nvucm+6x2MjIgIwOSWaFm\nwJKuZ7ks6w0qbe9hF1QO6BN55P12+hI9I89GxbwBRaADhNfgWfvjvFF1Bi2FEk6zbos35VN1vBzE\nxsCIMfXgx0cMMIjiIYKb45QQEV2EG5vILm0l/QmW5JYgsosKao3lJi3VwmPNWVbT6YGBKc2gxdOL\nxfDCdzK4TaIxGgjrBHzdPqKpxqwpo2oEEZJKqA6hC8UYg0IWybDSxy7cdBAjnzCm5oz4CQyeZEE9\naV0EOEQZ57DDOoY5wDSwJaAiY7d4902WNtAgdsaxChdBbCgM4eSYUExI9NItZHGF9g/yjRPDMiYB\n+RTwhwFokkSTUU2UYh4xXuafxh6CDKFajB0JeMZ/O/XhV2k3Tljv09N8+eFQiUx6Z6qRzb2mDIbV\naYwWg1gX1f4/pRzmq//Zxf5zlpNwZfOG7Vt4hXQfQVkDuFtAc0CV7wn2dZyBLEeZN28u0x74MT59\nHKY2jjN1UJEQIpIlCJfMF0003zwxlQZxBZvvuJi+WT4gBjvrQbgTzswBYYjqzE0Uix2pWkKVsC5F\noY3j4PmGR2lvmYTDEaW1tZTu7hI0TWZEFiYqNIrfoUxoIiqLfGloB3Z0NFHkXf23GHjMmZRUMzsI\nDDD/rCX0ZQSJDiVgg4ARy8FUhs/FTxZ5l31MtmuQLZk5sKEGOTGJgD6dAcGGP1MkPx+6uj6kt/fL\nVFaO+9yxfL5wjv/P7r9iHxIRRgxI4SyqSNDHTPlfNKvTGT1cj6ef6TPfJTyQBQ6Nu2v/IxUR1htL\niAhWlGKo0LcL+rthKB+iw7oTR1kZZQQIUMt9VqfkyEVDB0JksJ1ZzOEVfmL7AdNtewkqAWpL70X3\nG9TUzMXhHUQWzLYqw4DdTOcBHjGj0lATRI+CqyBd/FKGjWW4UJenDBI6JBSEhYuxdfciI2HDnsKZ\nR5oJEdloZSv1QBk6HkRiSMQwkAAbbg6Sx2o84gAfSqdzf+kD6H6DtrZS2tomk3QIoztmU18iMKUB\nXC2gW1z9fq+PBayhmwJLLfNkuCWlZrq3Lp1lDe8MHS2vPPw1SQcZPAD769KQwKexaj7NykY5XUEy\n8XZBplr8KcVCOwVGO7KgYDMUZBR0wYYNlUOUkkMvGQwgo/FEwxparfpQINDN9u1XWBuSyIiigWEI\nrDUWUsZBFL6BRh5OwkT5I91kMI+XkdBQsPEWs7mA9zjEOvSTdkIZaStnN/Cnij+Z6p0xzJ8katSM\nyaApBc6y5slZkCp6v+B/nGw5hoKNhxpeItESwOmM8qOam3B5B3l65Wq2brkyVeDu68vDMNKCe3uZ\nClzNEGMJcVrqmjmd79PFVy3sPr0p+0pW8gf+aP6xGDOTGX1qSbhJjfFCxuPk24JI1mI7+qUrVzbS\n3FxG/e4mMrCjonObdAklV32JeLQVWseQnaex7PGZVJzxa2y2AO+VvUf4aBh3hpvKHZW4xn8Se+qz\nZ184x//OLbcg/NvbrGyL8IfIgBkYxMD+rIAr4mRN2R3kxyZC0A/BDFNgasTgQ0zhDj52ZHO9+xlu\nDdm5XAsiYzCITBVnEXKqaJqN0tIP6ekZS5+Sg01TGV++mxuW1qShgtB+iB4znbJqRTvPREwFSyvl\n9mxwI8RMDRwNDRU1laL70XiTZnJGdXjGkZCAATI4V/oHR7QpgIAsx5k06UP8/v6UNv/y6INscl6K\nIY36egxL9VMfRQsd7uySjjAegw/qwYggdojYX7OTiCUQEVnLWkophVNEuyMthlmRNE0iDJhQlooP\nAxmJQWTCliYR+NhNGfV8wL0oFPAh01nKDJ4QFlFmpJ1eLm1M4yFsFnRmIBDGyzYqmcPrqd4KUVQI\nBE6g6xKbgxeRb3SjYKOSLZxwjUFVbbjdYTRNpqxsJ9nZbWzceHPq6jjzzD+ye/elKEpaOdPj6UdV\nJRIJD2uN76Ual75jX0t3YhyNLKCSfFTyGMTJTVxDO9n8mbdwWAqcfnbhoIs4BQxRYDm4dLYTodRy\n2mkbDpFkksm22dtIZCdwKHZ6M/tMyUOLRmoLSihoI5uzRn/HzgJzT4JRVtEQZ//Gi9FUG6KoMGnS\nR/j9/Sdr5KTqP6bDbhRvpYz9RCWJbPXfUI3Tre95dIHUYKX4OJuNzfQbQYzh9Y9i0tmN5cQlNcbP\nMx4n0zZk1cLWMLVBp7elmM6OEvJym+nuHoei2FEUG9GonzuNwxQT5SL+xnQewm5de+/yltWtnJwQ\nq9DviJuZQUYQMno5eOkTGBNCZOYkGDMG7PYCHI5xhA82wz8vgGcWYx9r5/yW80914X8m7Qvn+K+S\nj/N3l53o4BgMDHzeAX761EWc98sT9L0/iYFYgLmOF8koa6OqbhGSO8pD3E+V8BzF+gFigpv68PUj\nO10BcUjkgpZK9hQeICGplHeUU/faUrzDiqgnM941JAaJ0YoNDyo68/k6N8ltXKLuxU0Qw/0qPyxf\nwg3LliJ4FfoJUGBhlifp46QOpELXbvhBLZwfOaUM8mjzIuDHoASQKCGDMXTQQw7H+RU692LixnYk\n9rOTVipGfYKOxxPC6Yyk6IjJxSWEhwU8awq9GZYDSJLKDQ0S/SYs5CkFyXL+Shh2LDAjRvjkaPvT\novAL3jIbgADifRDrSFMuHfkmc2TvsPeUV4OzOP14clMa8/TMzVNsxWTG3BiblxP0tiEEBYwMA2fQ\nSSwvhmEzzPdEwH/EhzhRJCJGUDpUc6+44RI5QzBMXXn01gMjne7wzXJGMV1QIyB74TvnI7XbgDT1\ntcObQx/ZDBBg/fyf0N4+EU2TEQQzMzBt9O+RllycIBvJaj7bydMMUchI5pRpKfbR7L50Y5Zv5MuS\nG7ak9nuwQUxw00Q5y1hmFf3N2ki0JSeV3XR2lphCcD1/w24o1rSlYbStvECMQute06ysU+ZUtnL2\nSlrGNONEJKC42F60hz77oDnOYf0i1X1zOT/XgU3QOKhP4O+vDHDDNxVybSdY1bCGY8emYjg0blj6\nIH/oOMSLX60jPzCDjy9oRjvsRxQ9zNo7638j/v+PNsLxK1VVZO9qxug6itzRTpGRw6FYm9UdamcN\nT1MgeUAXULEj0kxMeo25ztUE4yWguFKnIYpx7E4dhy1EKJQDGIhuhbz7X2bykWrC8iD9Tuh3wsBp\n1ZBzXlo3Jwk9JGuZVtbx0+cfZ0bvWWBVDqxTINmHKliyBqYCZ4Q4v2AOz41i1kA1TRQzxCT2kCE+\nR0wU+Y7jWXKnHKZnaQdzvK8OY34s4WljEYV0gAC7jWnUBq8nsucxc4xJByVwyoVDQMCGDRl5hHIo\nmOn2ZjbTzwACAl48vMztvMljHETFDbwEaGTjJ4zNUnsZIINKdozoQ3C7g1ZTjwnpTJy4k5ycrtTi\n0tAALS3gdEIgAJ2dIDpsPFSr0Nhoiqupqh2vdwAqggwO5CA4NBx37wWfiiFA/MlyjDY3ODXISGDs\nzAJFBEUwD6sKMHEQwadiv2cP67JvZP0TdbS2lOF2mBumLF5sbpgi2xRufeIHXDRuI2tXPc6mt7+V\nOv6aNZXk5HewetVTbNl8NeFwFroOTmeU8vLt5OS009lZQkfHROJxG+FwrlUQNU2S4rhcAwwO5qau\nx6RzNshhGb+jmQ4ceKhlCS5k8thMkDNQ8SMQszIpPyCxkifYzBbChLFjx4kd7Wsaao5KuTiO7N97\n2XTdhySy1E9l6NzlXM+E6A50RSfPBWPStX1UBGKGg7WrVvPHP9yWejyt6+9mLi+bUtDCd7nYCKQ2\njzGwoeK1nPhISM/AYDaziQzLeqtXw+XTwC6Z5OogGXSoXnKlCDIqPiGCoFnQqCQhGjqiYJBIwOTH\nTUgxJsLr18IJBX73u0Z6+yaD3E/vUBm7Q1eTSy+9PpEXnevIVr3EDD+O6c9jFIQ54fcTWfh9Xrj2\nm4SbW0BwYZ9WwlU33cRvNm2C9evNi/QzbF8ox//dxxuIHDidUEAi5IOwDwZ9EPaa/x9yYUbLomRC\nHRiQiML222DijeZmEYIM4UFzB67ubjLwo6CmnN693EcfvSRQmEQpfjys4lIU/t2SF9AROE6E3zGX\nn6c49IosM5eX6FHzEEWdaRWbqV1+HW5vxNSYMZ6hW8g3pzEpfm9oJo69o2pkZOyyIuD9ZtTqanEx\n9VdTWRarO0VT1ckmMIifJuz0MEgxKmOw0cMgdtqZeooNW1LTDejIssbjay7itNJtI5avpPWRQQdj\nybX6hiN40tmLGjYXx111Jh/d+oAXZkG2AxQdFv4DuocXaSy8uXqWte2fAfX7IaJg4s4u0oXx5KZn\nLdZzAUYJev1XcyMgIeHBQ5RoWmztE7pwBQSM+4yR7Ji4+Xz1RVDshJgA9R9AZDTMgjm26ilQ7Bu2\nbaQVMKS2OUw+fs5bILupZhWX83ZKHyhIBh0Uml285DKdPawa1UdiW1yBt11lnNY6bKvQ9SOCiuTi\nEuFMdLwIxPFykKnU8TFPDtNkOrXdzHw6ZrejJjdWh5E1Ji1GtVI7Qjb76fOLKRZaPtV5lDfAwT9c\nRAQPpRyihOMIGLzLBchUoZ1UmwIwWDl71TAWWk2qCXFEn8R9pBHIZOtK4b0YJV8CWUKKR5n5y/vY\n6TTQXv0+RqwIyMQu55LQAmBkYgZuWdiIcxAYP2cOfMYZPl8ox/+l53/Jjh/WEIx1mFtqXfhNEMZB\nLHDyTzwDYhnm72GRPoDTFeLpdecRz3DRSw7T3NuRRZNl0Em+6ciU24n07oWsM0dG+jFrQ5HpJmxA\nIs6YO9cS6e3Gho2ZnEUvvcMalFyYyWwMEQEvTUxnGWHs/IUrKKSDKG5KOUQe3VbX7laaGY8oqkyt\n2ELX8m56vNmAAbrlNUQJDAHQTS48mgkbDJ6iH8BsILDwYBvJf6cu+I42E84ykNEtvFQkgkAMOyew\nESFGPioBdBwW8cdm1UtGLywnmyTHzaYdw8S8/f5edF0irrhRFRlDT+5yVZWKLG/zPsM16+5ndv4b\nOCy5BRWJBazlGKW8xWw8RBCAQdzsZ0pKgK3aWMlXhT+n9qyN4GYfU8iniyLaUsSmQdw0MYVecsin\nM7V3wJsN97Jt81dQEw5cjjC6JqOqJtasWx1BgqAhCDqCYO4iZkoPm9ffxIk7+dGRNUzmoBUlv0RY\n9HKnfoRv8VfsRJjKco4M0wUq4xEkwnzAz0c4Z4EQIgIicRR8pAFzE67ZyEZTjmK4lo6Xk9sWyqp5\nIf/XZAsDKMjsSEwgx6YyXdhnCd+l7VQCH4mESfm8/u3NTDKOpGi7K7iXebyMmyG28gJDFFuU4VYc\nioPaDbUsunUx7b52dElHjIvobh1D+JR6kgH3ToZZRVmpvROaKOcs4wNE6327BmdR7v4I1ASDgzBn\njpLeE0GKm2qdzj58vr2Mz/Ax98Fl/HjvZiIfTYDXz8DcE6ofmy2MonQwPKIYa7PR0t39vxH//2Mb\n4fir5m+i4tzHCEeP8MYbD7O/9SJ0l0jW4n3ExqhEPGbUb2bV1iWqq/CbO+H5iOm0CcEw3RaE5J6M\nw5twRLDPgMBPQS2E/9PemcfJUZV7/1tL792zZJbMJJlkkpAZspAgIYSwiIig3MuIoKACEhAlgKiX\nuLEkkpAIiAwRvIFE8AqvgooEEa/KJiogCSRkgex7MpOZySSZpWd6reW8f5yqXmaJ6Csy4e3f59Mz\n3VWnqp+qPvWcc57l96QikI70G0AA5i/6FGee9gdULBRFYKGxUUxhvuJW0HLlz7GJd7/N3NR3qCkd\nTlKJsLj9UmK/c6gUcsva5XCVTFkxhbuSd3IDX6GVFgQwhSksZtFRVgEyaifEXpScK4xSzGpO5lKe\ndJLOsggGu/EHeqmq2kt7+2iGOwUx3OzkG/khW8QUUFSepYEgcVmJiQfZy3hp2+9cDb5yadt2s017\ntsHmhdmkrncbWWMjn7vdQASUUoWIFqEn2EOGjz0nMjVzTJNsnykoLpzz7ELyAeXeMgGqUDlh/wks\n/tUiliWX8/vv/n7QHh9IBUikE3J/gKzPkz5yDORkdWbFs0pNirwqSeFlmzmSI517+Ui5hVeVIcLd\nBiieYjyKRRIfmjAoUXocg6HGTsYzlj3YaDzUuISL/riSz9u/wkuaLoo4RXmD/UptZjAC+hCxjSbA\nOfRwBgJPjqiCIjby3YbnaS5rpsVoxYgZmCUm9UY9bbE22krasAwLLaZhnWANlP+WxenPSv+FA7fP\nAEQJEyLB5EaToBvxdRv06pBIwI03bmZh231OvkeStepUPlO0hEXq9zlffZVkwEtLeBhXTvoi7aU+\neOTa/O/O1A2wYPHDqGe3EYlYFCU8FP/hJIrSIYrxUKq8Q/kp2+nwws9HAytgyoQpPPLoI3z43A/L\nAvWHIaiqbF63jjFTpx7lgocGPjCKXwjBl7/8ChvfmcP68jhaaRdxIeDp8jyFriLrdAbwEVR8WAJ6\nHAKxEBbTkfXN9/AhhnMILyn8SpL1+jS6AmFG2q2UTmjil/OnceFPNlDclGaLv55X5pVREXacstZc\nYrYX0h5ZSBpg40a49Vb4aKwfn3gRYHmPJ1ZxH1hVspCLEcr+t/s7r+bPv5izz35GXjuQFh40Raru\nbopRsPFgso06DlPODNbgJc12J8s4RjifGTTXpDR+LoRqZKz+lkWQiPUzcdzNXcxkJkCfxK55OQ5v\nNz9YcrwH2ItA62fjz7cDP9HPp9F3HhkMdhKPu21sioo66OkZlkkcyi+Qsor9jOUR9WrqlR14imN8\nvONlujNVPPsGw4s8J7YbIWKaHny+ONHoMAzDj6pa6Hoay/Jg2zqyH6qoqsmTJRdzVvVz7D44kc+a\nT9LcOxbL8qAoNhvtKVTTRhovL/ExRtDi0BdXcgm/xk8SC5W/cQb7GE0t+4gToJ0yLuFpmpD8MlEm\n4xZTby5rYlv1dlSh4DG9LHv4Iaq6K8j1JN/b0EhTWRPbq7ejCAVTMxl/qJYjkSMcCnXIW5DrLO/e\nCIkDMPxc2T9y7pPLoumucsaLnfiVVOYunjIbvEdAMUA1wZ2gWyioTpS/O8lI+P3Ecl5H9i3lggsu\nJVYWg5kOk2eJ4/wuHk9mBduxRtaj8JaSGVkazoC4Bios179EXWpXTp8q5gwm0UMrfmqZzy0U46OE\nNaSpoIc63EI/in8HQguDT8cOxzhwxgqiw9M8Xn8Fh/Uq1JRGUXw9Po/OT+5exH5PnM+3r2DKqR/l\nmWeCQ32yD3yAFD/ALX+6gGmff57k5FMZtf089M5KPCJFdXgJ0a7PYpnV6CSYxGI0YhhoJFU/iUlJ\nXrptDIGyHg6JSnqVEHda8+j0hrNX1bobxepEvH0HdCQhJWBFPr9PX0gKgDDLWcoIiikaIMHIxUpm\n8h/8kagS5snSi/hI9XMoIYu358G68FRKGqsY19TGcH8rf5lXSXc4x26OcPRrTEbN5NEfpKC1TeYZ\nOPkGCrAe8uJ2BPA3ZSYN2gq6reGOaaU/gsEuHnnkxEyWsA3YQnMGHeFct0xs28gJmdrAuZS7g0Yr\n5UmjZI8ROotX7yfWeVjWEOgBrpf2boQFb8yBaQuyrJ6db8kVRbhODm5xA965LuvHUXU50K2ZI+sA\nDzYLF0iHt2v+cCOncqN3epA2Yk/OMbnH20hm0ZEwtz7HZp++IctO6q5sBHDWn7IRUVYcxYgjFJW5\nvp9m7t8IDjCSFhQENgop/PhJogiRORQkVZJtybzEdBouv7yDeDxL54Bis1k4Va3wMCPyW9b1nE+L\ncj2dYibfb3iY5rImtpbuJHX8jTBM8gUFX16CVf4lUrMqwE4yt+h5anSnnkH3hcRC07IDiIvGekY2\nHSHi7+LO2y7kkR/fzXN/vDpTrB5kLL/BRaSpJobX8T35UDG4m0dYyao8ksBPciF5tM8GsMIP16Th\nNFmURTGT4B2FCA8HVUexwZeSFBeBBFwbaeS1ZZewetV5GIYPjyeOaQtMNc3wUW9w1rlfI2Em+eNf\n7icpfOjBLvRIF7FEEUVvfoE3mMF95+/itbJythgzuNB8gmd+OfQ1/wdK8b/6xhk89sR5BJonMumc\nz6HpNn6/ZNd0/6MpbE1P5c7kIhLxdoR1iKo3F5NUY8Q8EPNCQkeG/gVrIJ2ENYvhSEw6DnMVgBMB\nM5WpfI9FjOZtLIeZciw/IExHxsqwmkdJUYaKQeyH3yQ1tZ3tTgfOVq7qoxSFkHH2dgpWXwc72uVy\n5F3SEPe1cJwEVADf8J9A+pF9fOgXUQJNENkOWs4537wfek+QeZpWHwWevfOmNNPk/xzZbxaWEzGU\nlKsexXFa33gjbNkil9mznHqv27bBwoVHLd0o6ZA9RMmacHJXHbmpxVP4Kjv5NklG95HNZSRSM9vc\nBKE2zic/FjO7CpjKjWzi+1iZBCg33zkbC6oSRyOJjYpFBIU0Kb2V+82nOUBLprD9DdxIB0eIz70B\ntWYMVjIGixcNeO0ZGum+VafOBb78KPjLwDTgz9fDh78vBz5hgNErE9dUzQlvTYC3CDRf5vK+KRq5\n8r6VlDSleWfT6VxmP05UKXFCQLOP8hlMxGQ3q3IqiumoWH3qcPWT8RrkjL36XOkDU9z7mm/3yV09\n6BhEiPGl2RuY19SYVyXsEFcTd6KW3O/VSHI5X6GFVmxsggQwG0xSZen+JtEgMupzxqPgde7Tm9cR\niLdjO896MJ6i/tV7CBgdhHpVQu2VhFadyThGkQjGCJoG4XgPYQzCrCdc2UrQKygyDU79Ehg6fGrC\nJfzmsqHt2IUPmOJ/4IEwiqJRURkl4Jc6RfeCpuYLaKJwfedX2ZX0y+QqOyeuez95yjVEiAQJBDYa\nGr4Gi5jTwaetqOVHySuZSSMasUyXXr0EYgNQ7rucJh5/ihvnXU9POEhM+Fnc/UnmpBbJaAd/LYuV\nBcTUUmlq6W6D/TH0hXcTipnYWNQxgQUsoBidMNuYwsJ+q4gUKjqyOu0WvsXhopNRbA/Ub8e8fQEi\nLP0YrvNroOpSR7nt0PY6lB7vLLWFVDavOTPoo+QTTGEKnXQ6zJR6DiGaRTHvMIX56MSYiDyNF1iN\nnGQv41lOZh9eepnEYiDJ9Vj8liApBJMYwe9opZw4ApiDzE1IUcKVfJxP8Vd8pCgiiuYUmDTROEgF\nOjYeDOIE2EctD7KeLWhUEONhglRg0EMRZRxBAF2U5PAMSbv4rRxmP0kCKCwjQgUm+7mBBKO4nad4\nAZNTeZufoVFBms8Tposj7EaWWu4CZ+CQhHE6gghyUZHudyfJDp5FRVmOEIBNm6C5BaZXgpqUJHj6\nQPkgNs8qF2Zs6j1EOGd2D6EW+Ux1T4bjNrXRYVU6IcRxqkkSoZlq9jGFhbzO05mkLJ1OGmnkFTbT\nQVe/2rtz58Tkx7oTWGzfSszwQdADXrk6cAeB5xq/lmHmvHne5Sxd9kP++tIlmH3qCiznWi7gfxlO\nOzhO5u1OYZ3FPMF+2vAQ5m12E2tI5Q9MS+dCSY189jcuhtN/AVooc19QVBQBgc3NTL3l91w6503q\nTj+IN5TATEKsp4otr5xNaNdrtCtNLJ0BalEx+27ae0yUY/xAKf753/g0b4y+iF3BSpIBlZRXJeUT\npD0WY1+eg5LYS8IjZ/Qd0+Zihh0n4tuL0ZIpynaXEovE0A2dZSseoipZDc5MR3PI0nJNNQI4S3mB\n14TkOw+HO3nooZMYMWI/iuL6HQP4SCJQ0LAxnVllUnjYlq5koXKHVPLR7bB2ATTF5JPeJwxxJCM5\nzGE0VDzoLGMpIynJyKQRw0alkxL+m+cpRSOJyvcYjwilMU0vmmbh8aT6EbptYzzfZAlL+QplHMmU\nMGynirniB9T0vkzSSEnzhHeCrLK1SWbyZgqkDBAu6c5YdXSHo6eeBQ7VcjEbmcJ81rKUNGXYeNCI\nYRJGI0mYbfg4QpKqTFKRTowoQSJIIq5uIuhOHS8vaaIZxaxghmxeXSJI5QQlKUrWiuL90VzUAzUI\nX5L0tyUB0pQHoPJNUA3oHg/JMgi2Q6AN1MM6GoIuijlFeZ29oi7jS/CRwkfWxv0HPs5JrMdLOjOQ\nRCliPDsYziF8pOgmTAlRclcX63mYNGUIPKjEEehcyTXs55AzV1YpIUwH0YELmg+EBkXSOo+ZBXpE\nMqXaB6S5yXdXhme/1r+LMSXb6WwbQWvrOCpznPbLlt9Dc1MdLS3jOHx4JKBmlG41bTn9SNIiu7Ny\nlRRBdSeznjkeM+gFy5Jh0hdfLGWvrpb0uEEPdGwHDsPwj2RXJrl9yYJn1AZCSiyzklWQp0yn5euB\nHy2jpXU8bYdHUFK+iwMHJ7C0dz5npkcTZzS91CEpvTVQbS74zwtJlCUQpYJINExPZUz6anTwJcGf\nVBj2swl8uraICk89w6xSVvk28Monz+bwhA/jNVQCMQNv0qCzws/hSsH6k2cyNRIZ+LeBNPssAAAV\nYElEQVQYQvhAKf5P3zWPyaPepDzciukJcH9iPpcVLadabyKh+FkkbiOm5vwoPTF8X/4ORnsLdiYA\nW7qdIoRJksSHlwAqXaSwsAgT5DEm4VXDXKb/D/cYd1AndmYiIfrGR/v8UUZP3EL7rZ1ES/2I3JWu\nbciluKJBrxPVMlAUS5+qVJGOBAlfAkzwt/io+3WdU7s3BA4XiatYc1cCIvNfmnEUwIjAmuWQqsov\nenLbvMuIhLvzjnV/4F7hp8v2U0qXjLt/C9rH5GTYxjsh7NBUJDoh4vAIpbtAq5Lms5sW4z2SJp1j\nQuhrLtCSKsUNRcTKukgNEIsfwsMUxvIDPoeH8Xyv4QmaylrwGzq3r/gmVckDmdXQNuZykHMR6ChY\nhNnJwoanOVC2j4ChMTa2i20lVWhGBbetWEA4GQHsvFXIYHDLN7oEcsIZxAQ6Ebbzo4ZbeLxMkHKL\nnyRl/5rBDF7jb5iYA1YGy9wTV3G7qyuzF16dI30WI5ztjnMZ93kWloxIy3ksc+sDD/SwNjYuZ6XD\nz+8WjNm0aRY4ldQUZOZsmjIs/KiksPFlepaCCVgIgtzbcB8r61YS9UfBBNNrDh7dUzdX+l40SZ+S\nNSNmpayZXUVTUx1ZU1v+tYEAzQAtxVxrN/9pdWX2NjqJhj1EMdysyty+Vkq2qk2TF345EpbfCSPK\nAQXvxkP473+R3pZebvr6Wj52zh48HgvTzA46RzpKKR+TQvV6+PjMDQQCYxjK+EAp/qozz+Hw6+ei\nYrDkgTsoKTHxeuUq2OuVkwtdl7M+C+iyi2lVZOKLx+hAVwTbe+FQB1QPcxKF3IQaG9gPJb8qxkga\n6OiE8NFNzCmmsZRRhJjOHPZxRV7t1RQmcY+fZr2GEWYbXtKolmBbeBwdy/bgr4qSRMeLxZLGZZna\nt1+fN8eJFPKx3WG0zDhzbQtaNsKr86E8Jp2LSlbOfHNVkARJBAINjVmcSpSoUxCkiR5aSeJEFwEf\nAp5EPgvZ0M5f06uFsG0Z+xiJdDF9+ots2PAR0mkfE+re4hsLrmFu+B5J3XD0n00OYqr3XfzCtsw7\n0ItynJ5JKanql/4PRZfnVPT86CQ3RLJoMhmSt+h2WTM41S4pJNwoEbfM499DLhlbz3ZIHcoS5Lmk\ncH0pGHJDVEtOkuGrANFN8M7N2cH+wy9l5RSOs9zshc41UHaGoxjzI2uS+OkSxVQpB3micQGrV56P\nL51mYt2bVJY3s2bNeXR2VmYiniA/iqqdSs7jBXykWMtJXKL+nPvtPYwmjpKxyWcfbZeuIUoUPz6M\nBkPa1N3Q4mJgupPJ7nWKIcyeCUe8+D0xvvXQ5zi16q/4nc7ZQ4TrrOUcsarQP30SyXjWKSxyAgyW\ncy1X81NuwOJ/kZH0QWA6sA9oB3rJxpG50XtZ7eABvHD/DTByBFhJFKOLuWV/pUY9SNqCRcrtxEQp\nGkEMrwLC5qTHrqWsRcG3pgfvhE6enWJhWSCeknTVXm/2tXAh1NeD1zuK005r+vt96X3EB0rxl540\nk+/u1PCVVfPS4RfZEy/Frwa5w/tZSi0Ty6hEJ81x+l0ErG50N60beCswicAJO4h0qvj2hdmZriNK\nMZfxOFbYZvr0F+jsHIHfH6e4pJ3WtnF4/Cma5nU7yVN9JTMh1S2TqLwRSKSgYxv8bQFqPEF4RZjl\nyWVHSZASKCSQrrQ4oGHjQcFCxeAkrmMPs+nx1qAHYsQfWoxVJRVIfSOUrQQ1DT31sOl2WUPgoBhG\nWEniJ4Fwpl8bzTHMX9NGLBXP58V5yzHlyGRk+DEwiF5UnZpUfnx48bGMh/gV/80RVhHCphLYwkg8\nRBhGgG42AKNZzLXcz284xFuEsFnIDXiYkFesvP9dgU5KeZ7zuISnMiyUM1jFvqIaQqEuFrbdzQSx\nMy+UL5PUpFh5CkVR0gjh7fMNygDv3c8QCnXj8/dSXbWP9vbRVFbuZ1H7Ys6slDUb0iXgb4NdrZMp\nOpSgksOk8TJ74gOIfXNpT/Sg4+UqcR6zeTZTGEUAabzYKARIZXj3g5C5N7fzOw6wD51hTgJgJKOM\nDdJMoI4KynieF/v/UH0dsMk+++qAF5czMl1HJBDn0u9+j8bi7xAjAgge5SpG0Txg9S2A1OwzOLPp\nLYIkMnfOXQ0lqMZPOxoxEowgiTQZQZZwroVWuujKZks3KFAnyDBBuyUYw3UO86kGjcfLmgmtAahM\nQrtf/t8ZBlNFxUTYSj8m0Gy47+vs94yQ9j3zqyB2SlNT2Q/AHwCvAUYZul3CC4GPMeXATgIpg9sv\nWM6B8mqe+pzJsqXrGLNHx7Pwu6hFghkzNhdm/P+P+ATwQ2S4xSPA9/vsz1P8T5xyBmMPRRnRc4hh\nXd0ErAQ2Km/oJ7PLrHdioqUyMII2qgoT6tZy6bw7GFW0hxK1ExMP19kPcVCtImMcEWRnm0YUVl4L\ntVdATw0Ekyir7kSNxylaUYRICixMJlDH97mCGd7F+EQPumWBbWPhIUaIt5jOAUZSy96MTOmAwsSJ\nq1mw4BLCOWaWzMUCG+3judWaT6znAMqmhQgrNvidN51XC4M6W0E6XO/iTpaxjGaaM1nFRUCErXye\n79CKJKRYysNUMhY3Nj/EXmxUbFQ0LKJEWM0MLuXXrGRW5gHrpogKRwHKcoYjAQ1FsRFCRVUNhNAR\nAlTVzmOCdM1OuaaI+vq1DCs7QEvbeLb7x2OUKLJQSl8FELKyBb/fazTWw+8d2g0EPn8cy9JlUpXD\nAotiMnHqG9yx6CJKwh2A4MKGIyTiEYRQ8ft7SCZdxjMbny/GA6lvDJjn4NrZ3cgXHZNyDjvO51Je\n5Fwu5Lf8Fwm2I3PcRiNXdj8HilQVW1WxVJXd5tdJiFF8VHybOKnBDE8O3NA2Db6Sgpk3QeUsULyg\nenhW/Uwe13+UMNebP+Zy5XHOV3+fCfDKrZvc3V2GZXlRVRMQvG2elKnXqyAyjvo1uHFZxUxvOInX\nyv6GaZiEV4RIJ02SX0lkKZozeQl9HhABqg2alX2ppkAzbOb47qPG14whAryy9kZ2Fw8nqcss+LZq\nEF6DpSv+wIhAM72EeHOGTkuNyedWWHgD3YQ2zcGrjGbSE5PwlAzdUo1DWfFrwDbgY8ABZHDH54Et\nOW3yFP9PP3kpC//4Gw6ZAdIkmcA4KvVWqoixzTwOH2Fu5yZGsBWTcnSSTGQROjE6x0KiWGPdNy3s\nCEx7EAIHNPBZbLoV7CLJAqEoWVu5m7b/dyNihCkTYtbMh6WxTBIUgKfBg11mYxkWoRVB6pPHO5W4\n5Pnu5V5HGXu4i09j8hpj2MGd7GUHgnWcTpLfwABJT4piEYl0sWzZyRlnboa6QNTmx7LPeBQCI+Sh\nuRW5tiyG0XMgWgO9gzsSNTQs+uc01DGBTmcWF6OLEOAlRBUVBCmmDC9tDktoJTb7kDPcvsXEbSQ5\nRBovARI5AZnQQQnTWZdH+uYS2SVRWcwkYgMyOP4F+MgA2/tioJn/4N09d0b5EudwCU+hOwR8cfwE\nSeatKSxUtvMNbmYdr7MGA4OJjKSSEN1sYA82aWQEVgp5rw0MvHhJkcSPDx9+h9bbylR2AwaY5atI\nxa1Cgwl1ThlFHTmjrpglzWrxFPg15uoPMouVDHNCk13UN8KRlRMo7+lilzEh49vaw1hK6Oa1DOUx\naESxKM6rgzufeexkJydyIlc2zKaprknKYTq3tq++dk1s3mED3nvX9JWwfFTZB6gUbSi2TY8dYFe6\njgd6v0XcCtJz9YUZeo9x7MZLmkraEUBK9TLD8yc2eyex2jibMeo+bkrHedvr53BoOPf5axlrRAlH\nbc4OKnT2tmKLXiLV4yieu4CikSnu+qpOuKWEiksqmPzk5H5yDhUMZcU/C7gdOesHuNn5f3dOmzzF\n/0r5NLQiky987Rq6ykP06F3YgWrscLWUrusdeOe70BtjrkelplJSzTcuDlMTP47bPd8hQgTF8cDe\na95PsziAD19mWd0facLsRCfGSO5nHcczhv0ZXh03OzOteXjg1M/S1DMBny/O/PmX8dRT3Vx1lUx2\nus5Yyu4ln4D9GjQDRfvAF4VPXwb+7v53Nxu2LukHOslGAjnxy5qhMnnFZL6XvINiFEpYQ5RpmWpV\nk1iAZxCnZTdFvMXtbOFMevEOojwXOC/5GwSDUY6b+BaeYXF6D5YT9kcpKTnIm2s+gZH2YdSlsOdv\nQzw8HpoD4LNZOv8sRof35EUSHR0DKGFhORxFDh+R4smu0I6sAi0g7fCpdsk5rxdLe/mjj8FVs+UA\naDulz/SIDPFzKSQ+tLR/YpibdDVuDlSdn81steJgxuVnTyl52VTCli9VH+Q6YK5opEZpIomP79nz\niQs5+E9eksJqDrNr24momoXXicx67rn9XHml49cFvvjFzbS0jMO2dQKBaH6ylvudimCzmEw1rRQR\nzTPHqCTR6SJFFSop6tQ72SuuJSFGkaCKu3mMZuTzUEIxf+EVyfXjIpfzJ0CWLfmMZ0Fzku3+Ogel\n7mJExSxOfeIWPntVkMlszOP8mT17M0eOVOPxpDnxxJd49dVLchK95DUsF3O4hp+gOvJLx31+NFB2\n0uRjvnYzYTUEirT832v8iGbRQgsHMbCJNXTgH+ZljHE8zX++l1i6FOUjd6AU70U1fCgvL0Q95QFE\nyT7iigZPl6CGr2RUKEpHeQ8j113B7olBTjlnBz982c+0F6cVZvz/JD4DfBxwOV6vAGYCX81pk6f4\nf3DxJdz1qh8RuxBNszlhTBE7WnuIp0ziKYOAV6UnYQCCJSdfgyqkCV4YGnY6hPAlJQGUYiOEhkj5\nEZEotuK0A0TKBx5T1hZdf7JsM2EroqpVtnFXBCLnGFtHWCoiFUAUd+PSpL+6G04fnxNt47Z39ouc\nbYPtp2/73M99tuOeywLX5JnbfqDvzJzHyVkSAvjzhxDjW3gncZATJjltU5rc3xWGN45HGHreOfu+\nH2jbwO/y3w+8TTh/B+uCOa1n7IIRHaAKtq7VOL5lWp6sueh/tsG7uDKAItfRnPKKCkqODApuqGvO\ntuk7IJIES0VZVS+zgXLOp/T5D7CJfUyhdtD9ylGOVwbcr2SuQ8ncTQWmbUMZ3SqPEci+0F2EYugo\nayehmB44axVKKJHN0zI1lO4IytoTUEwdzlyFEkxk0966ilj1ZgmniQnyOydtRSk/jOLNRhwpa6eh\ndJRnPx+/BSUQRwkkJMkdAsVSURSBopmoQkOzdbxmCCul4Nk+leSklfKcig2WjqKZKJaKqlqZ61E6\nKlAi3SiawZpKnWdHCJZHZjI84JcDhWKzunoPUw6PwLOvltjPL+ML424iUH4tw7pPpPnkRcQfexEs\nHx8d3cVzG0JDWunDv1bxD/z0/PMY2IN0FMz+zcvcN6qWtsOy87btH4VqpuhMdAAKKdPtdgpbKkAV\nzidhoRBFdb5REaBgoQjZCVXhbgP8KecMFsr0t+D1maihGIqWbZN5ZY4zUVRQAmkUh/9NAYIGVMRz\nv5MMr4mSs22w/fyd9rnbj3aMi77nYID9AMqo9dAVJtULDdvdvZbTthtl+A5onuAcIzJHu+8H2pav\n5PKVZF8oA3SNd92D2+IQlfb+X/RaXFa+E1rHDdp88MFkoLZ/b7+S874PDsURccdEU7kdceC4QY/P\nDp4WF+cQCoo+7fp+T/5+WRv3aINyZn+sF5r6Tjyisk3NFsSeekR3GnpzJw0W0IUYsxmxqx4RS0Eq\nd2IRZWuJQU1njbOtF4w0wsg5R/0GxMrTsxMXTy8iFM2TQ2jZ9qZiYmKS9CURYRDjNyDCaYRuOW2k\n09hWrPxzFB1COFyM7aUG0ZHwtZ0beLLtawAoQuFD4nRsoeLtLsabCnP+1tN4KnwZh3tGMlnobLJ8\nTJuWZsVfSvAM/fytfyn+1TP+U5F2BNfUcwtyrpHr4N2JrAJaQAEFFFDAu8cuoP/sYghARwpXiwwq\nXA9MfD8FKqCAAgoo4L3H+cjInp3IGX8BBRRQQAEFFFBAAQUU8P8LPgFsBXYgK2cONdQAfwY2ARuB\nrznbhwEvIvNRXiA/XP0W5PVsBc77t0l6dGjAOuB3zudjSf4S4Clk7sdmZFTYsST/Lcj+8w4ytcHH\n0Jb/f4CDSHld/DPyTnfOsQO4/z2Uty8Gkv8HyP6zAXiabFVnGFryDyS7i2/gFgbOYijJ/q6hIc0/\ntUgSjqFo/68CTnTeh5Emq4nAPcC3ne3fIZuXMAl5HR7kde3k6EXr/l2YCzwOPOt8Ppbkfwz4ovNe\nRz60x4r8tchkWzc4/VfAbIa2/GciqZ5ylc8/Iq8bIPImcIrz/g9kAzzeawwk/7lk7+PdDF35B5Id\n5AT0OWAPWcU/1GR/15iFvBgXN5NN8BqqeAaZhbwVGO5sq3I+gxyBc1cuzyEjm95PjAJeAs4mO+M/\nVuQvRirOvjhW5B+GnCyUIget3yGV0FCXv5Z85fOPyltNfnb+54Bl74Wgg6CWgWfNABch2S5gaMpf\nS3/Zfw1MJV/x/0tl/3fOLkYi81VdNDvbhipqkaPxG8iH4KCz/SDZh2IE8jpcDIVrWgJ8i/xK88eK\n/GOBQ8BPgbXIisEhjh35O4BGJNdqC7JOy4scO/K7+Efl7bv9AEPjOkCuHv/gvD8W5L/QkeXtPtv/\npbL/OxX/P5zc9T4ijGRJ+TqyxEouMrkpg+D9vM4LkHyc6xg8R2Moy68jq00+6PyP0X9VOJTlHw/8\nF3LSMALZj67o02Yoyz8Q/p68Qxm3IQuhPfF+C/IuEQRuRdLeuPhX51oB/17FfwBpu3JRQ/5INVTg\nQSr9nyFNPSBnPS4hTTVZsuO+1zTK2fZ+4TTgk8gl4i+AjyKv41iRv9l5rXY+P4UcANo4NuQ/GXgd\nOIKkLnsaaeI8VuR38Y/0l2Zn+6g+29/v67gK+A/g8pxtQ13+8chJwwbkMzwKeAu54hrqsg+KYyG5\nSwH+D9Jckot7yNrXbqa/s8iLNFPs4j0aof8JnEXWxn8syf8KkmUeZBb4PRw78k9DRoMFHDkeA77C\n0Je/lv7O3X9U3jeQEVgK/34HYy358n8CGVlV3qfdUJS/lsH9EwM5d4eS7O8aQz256wykbXw90lyy\nDnkThyEdpgOFt92KvJ6tSIK6oYKzyEb1HEvyT0PO+HND8Y4l+b9NNpzzMeQKcijL/wukPyKN9MFd\nzT8nrxtSuBN44D2XOou+8n8RGda4j+wz/GBO+6Ekvyt7iuy9z8Vu8sM5h5LsBRRQQAEFFFBAAQUU\nUEABBRRQQAEFFFBAAQUUUEABBRRQQAEFFFBAAQUUUEABBRRQQAEFFFBAAQUUUEABBRRQQAEFDIT/\nC08lcvREjopCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff0ac72c410>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "from sklearn import linear_model\n",
    "x_mins = []\n",
    "y_mins = []\n",
    "for idx in range(len(M)):\n",
    "    mv = M[idx]\n",
    "    s = S[idx]\n",
    "    #plot(mv[:,0], mv[:,2], '.')\n",
    "    \n",
    "    X = np.vstack((mv[:, 0]**2, mv[:, 0], ones(mv.shape[0]))).T\n",
    "    y = mv[:, 2] / s\n",
    "\n",
    "    # Robustly fit linear model with RANSAC algorithm\n",
    "    model_ransac = linear_model.RANSACRegressor(linear_model.LinearRegression(fit_intercept=False))\n",
    "    model_ransac.fit(X, y)\n",
    "    inlier_mask = model_ransac.inlier_mask_\n",
    "    outlier_mask = np.logical_not(inlier_mask)\n",
    "\n",
    "    a, b, c = model_ransac.estimator_.coef_\n",
    "    a, b, c\n",
    "    \n",
    "    X_in = mv[inlier_mask, 0]\n",
    "    X_out = mv[outlier_mask, 0]\n",
    "    y_in = mv[inlier_mask, 2] / s\n",
    "    y_out = mv[outlier_mask, 2] / s\n",
    "    plot(X_in, y_in, '.')\n",
    "    #plot(X_out, y_out, '.')\n",
    "    t = linspace(X_in.min(), X_in.max(), 100)\n",
    "    plot(t, a*t**2 + b*t + c)\n",
    "    x_min = -b/(2*a)\n",
    "    x_mins.append(x_min)\n",
    "    y_mins.append(a*x_min**2 + b*x_min + c)\n",
    "    \n",
    "y_mins = array(y_mins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 572,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff06f517a90>]"
      ]
     },
     "execution_count": 572,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEmFJREFUeJzt3WuMHeV9x/HvYnPxenEWx+3aBHwpAQwqbSCE0kKbITIt\nuA24bUyIlGpDG5BKG6jaAqZq69MXVbgoCm1fVCKEyEU1wiSKa6qS2iCGqrSC0Jhws+NixS43LxC8\n3IICCtMXzxzv8dlzdmfOdffZ70canZnZmXNGs6PfPOc/c+YBSZIkSZIkSZIkSZIkSZKkWeFa4Cng\n6XwcYDGwA9gDbAeG+7NpkqQyfp4Q6McA8whBfhJwC3B9vswNwE192TpJUimfAe6omf5LQpjvBkby\neUvzaUnSDLca+AGh3DII/Bfw98DBmmUG6qYlSX0yf5q/7wZuJtTN3wGeAH5at0yWD5KkPpsu1AHu\nzAeAvwVeAMYIZZcDwDLglUYrnnTSSdnevXs7sJmSNKfsBT7ayopHFFjmZ/PX5cDvAJuBbcBoPn8U\n2Npwq/buJcsyhw4MGzdu7Ps2xDS4P92fM3kg3JDSkiIt9W8CHwbeB64G3iDc7bIF+ANgH3BZqxsg\nSeqcIqH+aw3mvQ6s6fC2SJLaVKT8ohkgSZJ+b0JU3J+d5f6cOQa6/P5ZXh+SJBU0MDAALeazLXVJ\nioihLkkRMdQlKSKGuiRFxFCXpIgY6pIUEUNdkiJiqEtSRLoe6mvXwvh4tz9FkgQ9+EUpZKxfD1u2\ndPmTJCkSM/oXpUNDcPCgrXVJ6oWetNQBW+uSVFA7LfWehPrZZ8OOHTA83OVPk6QIzOhQX78+Y8EC\n2L8fBgdh82bDXZKm0u1QvxH4PPAB8BRwBbAQuAdYwUTPR42q5lmWZSQJPPxwmGEZRpKm1s0LpSuB\nK4GzgDOAecDlwAZgB3AK8GA+3dTgYHhdsgReesnbHCWpW6YL9TcJfZMOErq+GwReAi4BNuXLbALW\nTfUmmzeHFvqpp8Ijj8D998NVV7W34ZKkyaYL9deBrwD/RwjzcUILfQQYy5cZy6ebGh4OJZdFi8L0\n2WfD7be3vtGSpMam63j6JOBPCGWYN4B7CfX1WhnV+xYbqFQqh8avvjphaCjh9tu9WCpJVWmakqZp\nR95rukL8Z4ELgS/m078HnAt8CrgAOAAsAx4CVjdY3z5KJamkbl4o3U0I8QX5B6wBngXuA0bzZUaB\nra18uCSps4qcCa4nBPcHwPcIrfZjgS3Acgrc0ihJKm5G//jIUJekcmb0A70kSb1jqEtSRAx1SYqI\noS5JETHUJSkihrokRcRQl6SIGOqSFBFDXZIiYqhLUkQMdUmKiKEuSREx1CUpIoa6JEXEUJekiBjq\nkhSRIqF+KrCzZngDuAZYDOwA9gDbAbuSlqQ+K9uzxhHAi8A5wJeA14BbgBuA44ANdcvb85EkldTL\nno/WAM8BzwOXAJvy+ZuAda1sgCSpc8qG+uXA3fn4CDCWj4/l05KkPppfYtmjgE8TSi31snyYpFKp\nHBpPkoQkSUp8pCTFL01T0jTtyHuVqdlcCvwhcFE+vRtIgAPAMuAhYHXdOtbUJamkXtXUP8dE6QVg\nGzCaj48CW1vZAElS5xQ9EywE9gOrgLfyeYuBLcByYB9wGTBet54tdUkqqZ2WeksrlWCoS1JJvbyl\nUZI0gxnqkhQRQ12SImKoS1JEDHVJioihLkkRMdQlKSKGuiRFxFCXpIgY6pIUEUNdkiJiqEtSRAx1\nSYqIoS5JETHUJSkihrokRaRoqA8D3wR2Ac8Cv0To+WgHsAfYni8jSeqjoqH+d8C/AacBv0DodHoD\nIdRPAR7MpyVJfVSku6QPATuBn6ubvxv4JDAGLAVSYHXdMnZnJ0kldbs7u1XAq8A3gO8BXyN0RD1C\nCHTy15FWNkCS1DnzCy5zFvDHwHeB25hcasnyYZJKpXJoPEkSkiRpYTMlKV5pmpKmaUfeq0jzfinw\n34QWO8D5wI2EcswFwAFgGfAQll8kqW3dLr8cAJ4nXBAFWAM8A9wHjObzRoGtrWyAJKlzip4JfhG4\nAzgK2AtcAcwDtgDLgX3AZcB43Xq21CWppHZa6i2tVIKhLkkldbv8IkmaJQx1SYqIoS5JETHUJSki\nhrokRcRQl6SIGOqSFBFDXZIiYqhLUkQMdUmKiKEuSREx1CUpIoa6JEXEUJekiBjqkhSRIn2UQugE\n403gp8D7wDnAYuAeYAXNO8mQJPVQ0ZZ6BiTAmYRAh9D59A5CN3cPMrkzaklSj5Upv9T3wnEJsCkf\n3wSs68gWSZJaVqal/gDwOHBlPm8EGMvHx/JpSVIfFa2pnwe8DPwMoeSyu+7vWT5IkvqoaKi/nL++\nCnybUFcfA5YCB4BlwCuNVqxUKofGkyQhSZLWtlSSIpWmKWmaduS9ivRWPQjMA94CFgLbgb8B1gA/\nAm4mXCQdZvLF0izLbMBLUhkDAwNQLJ8nr1tgmVWE1jmElv0/A18m3NK4BVhO81saDXVJKqnbod4O\nQ12SSmon1P1FqSRFxFCXpIgY6pIUEUNdkiJiqEtSRAx1SYqIoS5JETHUJSkihrokRcRQl6SIGOqS\nFBFDXZIiYqhLUkQMdUmKiKEuSREx1CUpIkVDfR6wE7gvn15M6IB6D6F7u+HOb5okqayioX4t8CxQ\n7cZoAyHUTwEeZHLfpJKkPigS6icAa4E7mOhe6RJgUz6+CVjX+U2TJJVVJNS/ClwHfFAzbwQYy8fH\n8mlJUp9NF+q/BbxCqKc36wQ1Y6IsI0nqo/nT/P1XCKWWtcAxwCLgLkLrfClwAFhGCP6GKpXKofEk\nSUiSpJ3tlaTopGlKmqYdea9mre9GPgn8OfBp4BbgR8DNhIukwzS+WJplmY14SSpjYGAAyuXzIWXv\nU68m9E3AhYRbGj+VT0uS+qylM0EJttQlqaRettQlSTOYoS5JETHUJSkihrokRcRQl6SIGOqSFBFD\nXZIiYqhLUkQMdUmKiKEuSREx1CUpIoa6JEXEUJekiBjqkhQRQ12SImKoS1JEpgv1Y4BHgSeAZ4Ev\n5/MXAzsIPR9tJ3RnJ0nqsyI9awwCPyZ0Uv2fhH5KLwFeI/RVegNwHPZRKkkd0e2ej36cvx4FzAMO\nEkJ9Uz5/E7CulQ+XJHVWkVA/glB+GQMeAp4BRvJp8teRrmydJKmU+QWW+QD4GPAh4N+BC+r+nuVD\nQ5VK5dB4kiQkSVJ2GyUpammakqZpR96rbM3mr4B3gS8CCXAAWEZowa9usLw1dUkqqZs19SVM3Nmy\nALgQ2AlsA0bz+aPA1lY+XJLUWdOdCc4gXAg9Ih/uAm4l3NK4BVgO7AMuA8YbrG9LXZJKaqel3tJK\nJRjqklRSt29plCTNEoa6JEXEUJekiBjqkhQRQ12SImKoS1JEDHVJioihLkkRMdQlKSKGuiRFxFCX\npIgY6pIUEUNdkiJiqEtSRAx1SYqIoS5JESkS6icS+iB9BngauCafvxjYAewBtjPR7Z0kqU+K9Kyx\nNB+eAIaA/wHWAVcArwG3ADcAxwEb6ta15yNJKqnbPR8dIAQ6wNvALuAjwCWE/kvJX9e1sgGSpM4p\nW1NfCZwJPAqMAGP5/LF8WpLUR/NLLDsEfAu4Fnir7m9ZPkxSqVQOjSdJwubNCXv2wN69sGIFLFoE\nmzfDsBV5SXNUmqakadqR9ypaszkS+FfgfuC2fN5uICGUZ5YRLqaurltvUk09SeDhhw9f6OijYXAQ\nPv5xuPdeA17S3NbtmvoA8HXgWSYCHWAbMJqPjwJbi3zg4GB4XbQovA4NwU9+AgcPwgMPwMknw9q1\nMD5e5N0kSbWKnAnOB/4DeJKJEsuNwGPAFmA5sA+4DKiP4kkt9fFxuOoquPVWuO66iTAHWLgQ3nkn\njC9ZAsccY4lG0tzTTku9pZVKmPaWxvFx+MIXYGAA3n47BPzQUBivtXQp7NplsEuK36wO9VrVVny1\n9b5oEbz55sTfbb1LmguiCfWq2hLNuefCgQONW+8LFsBZZxnwkuISXajXmq71XrVkCXziE4a7pNkv\n6lCvatR6rwZ8bSve2ruk2W5OhHqtqe6gAWvvkma3ORfq9cbH4bTTmtfebb1Lmk3mfKiDd85Iioeh\nXsM7ZyTNdoZ6E0XvnPHZM5JmEkN9GkXvnAFvjZTUf4Z6CUWfPWPrXVK/GOptKPLsGS+ySuolQ71D\n6mvwcHjrvcqLrJK6yVDvsEatdx9PIKlXDPUuKnqR1da7pE4x1HtkuscTVHmRVVI7uh3qdwK/CbwC\nnJHPWwzcA6ygea9HEFmo16t9PMFUt0h6kVVSGd3uo/QbwEV18zYAO4BTgAfz6TlneDg8U2b9enjy\nyfB67rkTf1+4EF57DV54AR55BO6/P5wE7H9VUrcUPROsBO5joqW+G/gkMAYsBVJgdYP1om6pN1Lk\nIqutd0lT6UVNfSWHh/pB4Lia93i9ZrrWnAv1WmWfQ7N/v0Evqb1Qn9+Bz8/yoaFKpXJoPEkSkiTp\nwEfODsPDsGVLGN+1q/lzaN59N5RnIJRqIJRpfFywNDekaUqaph15r3bKLwlwAFgGPITll0KmukXS\nMo0k6P6F0ka2AaP5+CiwtcX3mXOqrfcVKyZfZH3yydChB4QyTf1F1uOPh/PPhxNPDK9r13rRVdLh\nipwJ7iZcFF1CuDD618C/AFuA5czhWxq7oejjgqv80ZMUH398FKEyZZqqarnm/ffhvff88ZM0Wxnq\nkav/JWuR58JXWZeXZh9DfQ5q9siC2lb8VE+Y9PZJaeYy1HVYyF97bbEnTFZZl5dmFkNdDbVSl68+\njGzBAli1yqCX+sFQ17TaqcvXB72lG6m7DHW1ZKpHCR97LLz1VvOgr7JGL3Weoa6OqH0Y2W23NQ/6\nojX6/fsnbq+0nCMVZ6iraxoFfZkafSOWc6SpGerquSI1+tqgb6WcYytfc5WhrhmhUdBXb69sp5xT\nq1Er31/QKjaGumaNVso5RVv5MPlRCYa/ZiNDXbNeJ1r5jX5B20zt4xMalXo8CaifDHXNCc1a+VP9\ngrbd8K8q8g3A2r86xVCXaPyohCLh342TwFS1f78NaDqGulRCkVJPK98AytT+GynzbcBvCHHrZ6hf\nBNwGzAPuAG6u+7uhrlmv6DeAorX/Tn0baKTsNwRPEjNTv0J9HvADYA3wIvBd4HPArpplDPUOSdN0\nTnXa3W3d3J/T1f478W2gk98QptLsh2L1J4HFi1PGxxNPIB3STqi345eB79RMb8iHWpk6Y+PGjf3e\nhKjMtP158GCWrV+fZfv2Zdmll2bZunVhvMy8NWuyDMJw7LHhddGi1ucNDU3Mm37YWGLZxsPRR2fZ\nccdl2fHHZ9l552XZCSdk2chIf+etWRP+N70G9KU1/BngazXTnwf+oW6Z3u+NSM20EJrtYtyfBw+2\nfkIoe5KYfBLY2MMTSG+HJUvaP0mcd16WXXxx8RMEfQr138VQ75kYQ6if3J/Ta3aSaHQSOP30jTPq\nW0an5i1c2NkTxPr1xfY9bYR6OzWbc4EK4WIpwI3ABxx+sfQ54KQ2PkOS5qK9wEd7/aHz8w9eCRwF\nPAGc1uuNkCR1zsWEO2CeI7TUJUmSJM1kFwG7gf8FbujztsxW+4AngZ3AY/m8xcAOYA+wHfCu3ubu\nBMaAp2rmTbX/biQcr7uBX+/RNs4WjfZlBXiBcHzuJHxrr3JfTu1E4CHgGeBp4Jp8/ow9PucRyjEr\ngSOx1t6qHxL+ybVuAa7Px28AburpFs0uvwqcyeFB1Gz/nU44To8kHLfPAUf0ZCtnh0b7ciPwpw2W\ndV9ObynwsXx8iFDCPo0ZfHwW+VGSpvdD4MN183YDI/n40nxaza3k8CBqtv9u5PBvlN8h3N2lCSuZ\nHOp/1mA592V5Wwm/zO/I8dmNtP8I8HzN9Av5PJWTAQ8AjwNX5vNGCF+DyV9HGqyn5prtv+MJx2mV\nx2wxXwK+D3ydiVKB+7KclYRvQY/SoeOzG6Hel19CReg8wj/7YuCPCF+Ba/XtV2eRmG7/uW+n9o/A\nKkIZ4WXgK1Ms675sbAj4FnAt8Fbd31o+PrsR6i8SLgRUncjhZxkV83L++irwbeAcwtl7aT5/GfBK\nH7ZrNmu2/+qP2RPyeWruFSaC5w7C8Qnuy6KOJAT6XYTyC3To+OxGqD8OnMzEj5I+C2zrwufEbBA4\nNh9fSLja/RRhP47m80eZOBhUTLP9tw24nHC8riIcv49NWlu1ltWM/zYT9Xb35fQGCCWrZwmPLq+a\n0cenP0pqzyrC1e4nCLc8VffhYkKd3Vsap3c38BLwHuEazxVMvf/+gnC87gZ+o6dbOvPV78vfB/6J\ncMvt9wnhU3t9x305tfMJj1R5golbQi/C41OSJEmSJEmSJEmSJEmSJEmSJElSN/0/zlEPwR4ROIAA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff075637250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(y_mins, '.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 577,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "10\n",
      "20\n",
      "30\n",
      "40\n",
      "50\n",
      "60\n",
      "70\n",
      "80\n",
      "90\n",
      "100\n",
      "110\n",
      "120\n",
      "130\n",
      "140\n",
      "150\n",
      "160\n",
      "170\n",
      "180\n",
      "190\n",
      "eps_hat=10.1156770381, f_hat=853.105091424\n",
      "avg diff is 0.00018717084228\n",
      "std diff is 0.000215140241104\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEACAYAAAC57G0KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEdxJREFUeJzt3X2QVfV9x/H3urugK+Iusw5WQHCIieYPNRYVWzPeVBrR\ndOqkjXUwaWucFv+oqfRBAZsZ9o+2kdhaJ+OkilHHPqC1xunojGC18SaV1iitolZBwEB4EBoMi4gV\nQW7/+J2Fu5f7/LT3/Pb9mtnZex44+11m93N+53t+Zy9IkiRJkiRJkiRJkiRJkiSNWw8Cu4HXy+zz\nHWAjsA74XDuKkiQ17vOE0C4V8FcDTyevLwFebEdRkqTmmEXpgL8XuC5veT0wtdUFSZLKO6EJx5gG\nbMtb3g5Mb8JxJUkNaEbAA3QVLOeadFxJUp16mnCMHcCMvOXpybpRZs+endu8eXMTvpwkjSubgU/V\n8w+bMYJ/Evid5PVcYJgw62aUzZs3k8vlUvuxbNmyMa9hvNaf5tqtf+w/0l4/MLvecK5mBP8IcDkw\nSOi1LwN6k233EWbQXA1sAg4AX6+3GElS81QT8Auq2OfmRguRJDVXs26yRi+TyYx1CQ1Jc/1prh2s\nf6ylvf5GFM5+aaVc0k+SJFWpq6sL6sxqR/CSFCkDXpIiZcBLUqQMeEmKlAEvSZEy4CUpUm0N+NMW\nXc3W3cPt/JKSNG61NeD3DKzisuUL2/klJWncamvAd304yL4jOx3JS1IbtDXgJx38DPsH1jiSl6Q2\naGvAT2QyAH3Dc3hh8Yp2fmlJGnfa+rdotuzay2XLF/LC4hVc+Ve3sevQ2/TSx9qlK5k5tb+NpUhS\nOjTyt2jG7I+N9S/KsG/ghwBM33ct2+56rI2lSFI6pPKPjfXSB9iukaRWGbMR/Nbdw1y2fCETTjiJ\n9w5vtVUjSUWkskUzwlaNJJWWyhbNCFs1ktQaYx7wa5euZPq+a3lzybO2ZySpica8RSNJKi3VLRpJ\nUmsY8JIUKQNekiJlwEtSpAx4SYqUAS9JkTLgJSlSBrwkRcqAl6RIGfCSFCkDXpIiZcBLUqQMeEmK\nlAEvSZEy4CUpUtUE/HxgPbARWFxk+yCwGngVeAO4oVnFSZLqV+mPyHcDG4B5wA7gZWAB8FbePkPA\nRGApIew3AFOBwwXH8g0/JKlGrXzDj4uBTcAW4BDwKHBNwT7vApOT15OB9zg+3CVJbdZTYfs0YFve\n8nbgkoJ97gd+AOwETgF+q2nVSZLqVingq+mp3E7ov2eA2cCzwPnA/sIdh4aGjr7OZDJkMpnqqpSk\ncSKbzZLNZptyrEp9nbmEHvv8ZHkpcARYnrfP08BfAGuS5X8j3IxdW3Ase/CSVKNW9uDXAmcDs4AJ\nwHXAkwX7rCfchIVwc/UzwDv1FCNJap5KLZrDwM3AM4QZNQ8QZtDclGy/D/hL4CFgHeGEcRvw81YU\nK0mqXl3D/jrZopGkGrWyRSNJSikDXpIiZcBLUqQMeEmKlAEvSZEy4CUpUga8JEXKgJekSBnwkhQp\nA16SImXAS1KkDHhJipQBL0mRMuAlKVIGvCRFyoCXpEgZ8JIUKQNekiJlwEtSpAx4SYqUAS9JkTLg\nJSlSBrwkRcqAl6RIGfCSFCkDXpIiZcBLUqQMeEmKlAEvSZEy4CUpUga8JEXKgJekSBnwkhQpA16S\nImXAS1Kkqgn4+cB6YCOwuMQ+GeAV4A0g24zCJEmN6aqwvRvYAMwDdgAvAwuAt/L26QfWAFcC24FB\nYE+RY+VyuVyj9UrSuNLV1QWVs7qoSiP4i4FNwBbgEPAocE3BPtcD3yeEOxQPd0lSm1UK+GnAtrzl\n7cm6fGcDU4DngbXAbzetOklS3XoqbK+mp9ILXAhcAfQB/wm8SOjZS5LGSKWA3wHMyFuewbFWzIht\nhLbM/yUfPwLOp0jADw0NHX2dyWTIZDK11itJUctms2Sz2aYcq1Ljvodwk/UKYCfwEsffZD0HuIdw\nk3Ui8GPgOuDNgmN5k1WSatTITdZKI/jDwM3AM4QZNQ8Qwv2mZPt9hCmUq4HXgCPA/Rwf7pKkNqvr\nrFAnR/CSVKNWTpOUJKWUAS9JkTLgJSlSBrwkRcqAl6RIGfCSFCkDXpIiZcBLUqQMeEmKlAEvSZEy\n4CUpUga8JEXKgJekSBnwkhQpA16SImXAS1KkDHhJipQBL0mRMuAlKVIGvCRFyoCXpEgZ8JIUKQNe\nkiJlwEtSpAx4SYqUAS9JkTLgJSlSBrwkRcqAl6RIGfCSFCkDXpIiZcBLUqQMeEmKlAEvSZEy4CUp\nUga8JEWqmoCfD6wHNgKLy+x3EXAY+I0m1CVJalClgO8G7iGE/GeBBcC5JfZbDqwGuppZoCSpPpUC\n/mJgE7AFOAQ8ClxTZL9vAI8DP2tmcZKk+lUK+GnAtrzl7cm6wn2uAf42Wc41pzRJUiMqBXw1YX03\nsCTZtwtbNJLUEXoqbN8BzMhbnkEYxef7RULrBmAQuIrQznmy8GBDQ0NHX2cyGTKZTE3FSlLsstks\n2Wy2KceqNNruATYAVwA7gZcIN1rfKrH/Q8BTwBNFtuVyObs3klSLrq4uqLMzUmkEfxi4GXiGMFPm\nAUK435Rsv6+eLypJar129ssdwUtSjRoZwfskqyRFyoCXpEgZ8JIUKQNekiJlwEtSpAx4SYqUAS9J\nkTLgJSlSBrwkRcqAl6RIGfCSFCkDXpIiZcBLUqQMeEmKlAEvSZEy4CUpUga8JEXKgJekSBnwkhQp\nA16SImXAS1KkDHhJipQBL0mRMuAlKVIGvCRFyoCXpEgZ8JIUKQNekiJlwEtSpAx4SYqUAS9JkTLg\nJSlSBrwkRcqAl6RI9Yx1AfnOuXUhuw69zUd9m7lg1kymnDyZlb+5kv4T+8e6NElKna42fq1cLpcr\nu0P/ogz7Bn44al3Px4P0H7iItUtXMnOqQS9pfOnq6oI6s7raFs18YD2wEVhcZPtXgXXAa8Aa4Lx6\niumlL7z4aHL4fHAShyfsYc/AKmZ/+0L6F2U4bdHVbN09XM/hJWlcqeas0A1sAOYBO4CXgQXAW3n7\nXAq8CewjnAyGgLkFx6k4gt+6e5jLli9k5oY7WdN3K92n7OWTmc/RNzyH7txE9g+sCTt+OMgJR07k\n5EMzmchkR/eSotXICL6af3QpsIwQ3ABLks93lNh/AHgdmF6wvmLAjxgehoUL4Zt/PsyX7l3IC4tX\nMOdb17NnYBUcnAQTPxi1f/f7ZzHpkzPppc+wlxSVVgf8V4Argd9Plr8GXAJ8o8T+fwp8GlhYsL7q\ngC9mZHR/4Mhe9g48Bwcnw8T3i47uBw9exOSe03jv8FZDX1KqNRLw1cyiqSWVvwDcCPxysY1DQ0NH\nX2cyGTKZTNUHnjm1n213PXY06B+58U4WPHgrLywJo/v9EEb3fXvY07eKPR8Owil7AJh119kMHvRG\nraTOl81myWazTTlWNWeFuYSe+kiLZilwBFhesN95wBPJfpuKHKehEXw5haP7vuE5TMz1JyP9vJaO\no3tJKdPqFk0P4SbrFcBO4CWOv8l6JvADQvvmxRLHaVnAjxgJ+hcWrwAoaOmMDnr69hx97eheUqdq\ndcADXAXcTZhR8wDwLeCmZNt9wPeALwM/TdYdAi4uOEbLA76Yakf3Izdq9/dudnaOpI7RjoBvhjEJ\n+BHlRvfH3agdYUtH0hgz4OuUH/rHpmGG2Tm2dCR1AgO+CQpn59jSkdQJDPgWsKUjqRMY8G1iS0dS\nuxnwY8CWjqR2MOA7QCMtHYNeUikGfIeqtqXj6F5SKQZ8CpRr6XjDVlIpBnwKNXLD1tCXxg8DPuVq\nvWFbGPoTu0/0PWylSBnwkal0w7Zo6CcmHTqL7g988xMpFgb8OFD2L2V+NBlOfB92zKG7ayKfnOGb\nn0ixMODHqcL3sL1g+wq2XXw9701ZVbal4/vZSulhwI9zI+9hu2IF7DtY5ZufJJyiKXU2A17HKf/m\nJyXezzZh6Eudw4BXVY57P9siUzSdly91FgNedasm9J2XL40dA15N1ei8fG/iSs1jwKulGpmXbz9f\naowBr7Zr5CZuYWvH8JdKM+DVERrq5ycc8UujGfDqWNX18x3xS6UY8EqNYq2dZo74D/RuNfwVFQNe\n0WjKiD9RrN3jlE6ljQGvaNU+4q/c7ik1pdPwVycy4DVuFY74Kz2lW25KZ7mHuGwBaawY8FIRxWb1\nQOkpnRUf4ipQrv/v1YCaxYCXalRN+Bfr+1fb/6/1amDevpW8+5N++vpg5Uro93yghAEvNVG5vn+1\n/f9arwZO2HcWR/aeCQOb6T4wk4GTJpc9IXiFMH4Y8NIYKNf/h9quBnq7JrLv1DI3hAvVeIXgCSG9\nDHipQ1V7NfB7q67nua2r6D40mU96C68CGr9CqPeE8NEHfVywaSVT+vptHY0RA15KueGPhln41EK+\necmdfOnuWyueEEa2VXuF0NAJ4cAgHD6xZOso/4TgbKPmM+ClcaraK4SRbTWfED6eBBPKTCctty5R\n7oGzclcP/kmKwICXVJVaTwgzz9nLmnefK9s6qmu2UbVXD4lqpqTGekVhwEtqiWpaR7XONqq+nVTH\nlNRy6xJpO1m0OuDnA3cD3cD3gOVF9vkOcBXwIXAD8EqRfQx4aRwq98xBuauH2qekNnBFMaLGk8X0\nfdey7a7HmvZ/VUwjAV9JN7AJmAX0Aq8C5xbsczXwdPL6EuDFEsfKpdnzzz8/1iU0JM31p7n2XM76\nm2HLrr256X90be7fX98y6vOWXXuP21a4bvArl4/aNnjLVTmGyLF08qjPfYvm5AZumVdxXf62Lbv2\ntvx7B1o2Mr4UWJ23vCT5yHcvcF3e8npgapFjtfw/opWWLVs21iU0JM31p7n2XM76x1ph/Y2cLAq3\ntQMNBHxPhe3TgG15y9sJo/RK+0wHdtdblCS1ysyp/UfbKoWfq13X6rZMs5xQYXu1Z47C/pDNdkka\nY5Ua93OBIcKNVoClwBFG32i9F8gCjybL64HLOX4EvwmYXX+pkjQubQY+1YoD9yQHnwVMoPJN1rmU\nvskqSeowVwEbCCPwpcm6m5KPEfck29cBF7a1OkmSJEnNNZ/Ql98ILB7jWqoxA3ge+B/gDeAPk/VT\ngGeBt4F/BTr9WeduwgNnTyXLaaq/H3gceAt4kzBzK031LyX8/LwOrAQm0rn1P0i4X/Z63rpytS4l\n/C6vB77YphrLKVb/nYSfnXXAE8CpedvSUP+IPyHc85ySt66j6q/mQalOczpwQfJ6EqE9dS7wbeC2\nZP1i4I72l1aTPwb+EXgyWU5T/Q8DNyavewi/oGmpfxbwDiHUAf4J+F06t/7PA59jdMCUqvWzhN/h\nXsL3uYnKM/FarVj9v8qxuu4gffVDGGiuBn7CsYDvuPqreVCq0/0LMI/RD3Cdnix3qunAc8AXODaC\nT0v9pxICslBa6p9CGBQMEE5OTxECp5Prn8XogClV61JGX4WvJkysGGuzKD4CBvgy8A/J6zTV/8/A\neYwO+Jrrb3X6F3sIalqLv2YzzSKcXX9M+IEfmfq5m+JP63aKvwFuJVzejUhL/WcBPwMeAv4buB84\nmfTU/3Pgr4GfAjuBYUK7Iy31Q+lazyD8Do9Iw+/zjRyb5ZeW+q8h1PZawfqa6291wKf5gadJwPeB\nW4D9Bdta+vchGvRrwP8S+u+lnnPo5Pp7CDOxvpt8PsDxV32dXP9sYBFhcHAG4efoawX7dHL9hSrV\n2snfx58BHxPug5TSafX3AbcDy/LWlXteqWz9rQ74HYRe0ogZjD4DdapeQrj/PaFFA2Ekc3ry+hcI\nIdqJfgn4dcKl3SPArxC+j7TUvz35eDlZfpwQ9LtIR/1zgP8A3gMOE27yXUp66ofSPyuFv8/Tk3Wd\n6AbCMzpfzVuXhvpnEwYH6wi/w9OB/yJcRXVc/dU8KNVpuoC/I7Q58n2bY/2vJXTOTbJyLudYDz5N\n9f8I+HTyeohQe1rqP58w++okws/Sw8Af0Nn1z+L4m6zFah25yTeB0ErbTHvfU6KUWYyufz5hFtNg\nwX5pqT9fsZusHVV/sQelOtllhN71q4Q2xyuEH5gphBuXnTbNrZzLOTaLJk31n08YwedPc0tT/bdx\nbJrkw4Qrwk6t/xHCvYKPCffLvk75Wm8n/C6vB65sa6XFFdZ/I2Ea4VaO/f5+N2//Tq3/IMf+//O9\nw+hpkp1WvyRJkiRJkiRJkiRJkiRJkiRJkiSV9//TTMrXWc9+EQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff06f680990>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.optimize import minimize\n",
    "import sys\n",
    "\n",
    "def FrameDelta(X, rot_hz, epsilon, s=1, dt=1.0/60.0, w=1280, f=1280):\n",
    "    X = X - w/2 + epsilon\n",
    "    phi_dot = 2*pi * rot_hz\n",
    "\n",
    "    Theta1 = arctan(X / f)\n",
    "    Theta2 = Theta1 + s*dt*phi_dot\n",
    "    return f*(tan(Theta2) - tan(Theta1))\n",
    "\n",
    "def GetInliers(X, y):\n",
    "    X = np.vstack((X**2, X, ones(X.shape[0]))).T\n",
    "\n",
    "    # Robustly fit linear model with RANSAC algorithm\n",
    "    model_ransac = linear_model.RANSACRegressor(linear_model.LinearRegression(fit_intercept=False))\n",
    "    model_ransac.fit(X, y)\n",
    "    return model_ransac.inlier_mask_\n",
    "\n",
    "class MinFun (object):\n",
    "    def __init__(self, X, y, s):\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "        self.s = s\n",
    "        \n",
    "    def __call__(self, args):\n",
    "        rot_hz = args[0]\n",
    "        epsilon = args[1]\n",
    "        f = args[2]\n",
    "        y_hat = FrameDelta(self.X, rot_hz, epsilon, self.s, f=f)\n",
    "        diff = y_hat - self.y\n",
    "        #return diff.dot(diff)\n",
    "        return sum((abs(diff) + 1))\n",
    "\n",
    "tick = True\n",
    "\n",
    "thetas = []\n",
    "for idx in xrange(len(M)):\n",
    "#for idx in xrange(90,110):\n",
    "    #idx = 90\n",
    "    s = S[idx]\n",
    "    X = M[idx][:, 0]\n",
    "    y = M[idx][:, 2]\n",
    "    c = C[idx]\n",
    "    W = 1280\n",
    "    #inliers = GetInliers(X, y)\n",
    "    #X_in = X[inliers]\n",
    "    #y_in = y[inliers]\n",
    "    X_in = X\n",
    "    y_in = y\n",
    "\n",
    "    minfun = MinFun(X_in, y_in, s)\n",
    "    x0 = [sqrt(0.0577), eps_hat, f_hat]\n",
    "    #x0 = [sqrt(0.0577), 0, 1000]\n",
    "    if tick:\n",
    "        res = minimize(lambda x: minfun([x[0], eps_hat, f_hat]), x0, method='Powell')\n",
    "    else:\n",
    "        res = minimize(minfun, x0, method='Powell')\n",
    "\n",
    "    #print '{} rot_hz={:.4f}, eps={:.4f}, f={:.4f}'.format(res.success, res.x[0], res.x[1], res.x[2])\n",
    "    #print res\n",
    "    \n",
    "    if res.success:\n",
    "        thetas.append([idx, s, c, res.x[0], res.x[1], res.x[2]])\n",
    "        if idx % 10 == 0:\n",
    "            print idx            \n",
    "    else:\n",
    "        print 'Fail at', idx\n",
    "\n",
    "    #N = 200\n",
    "    #t = linspace(0, W, N)\n",
    "    #y_lin = FrameDelta(t, res.x[0], res.x[1], s, f=res.x[2])\n",
    "    \n",
    "    #plot(t, y_lin)\n",
    "    #plot(X_in, y_in, '.')\n",
    "\n",
    "thetas = array(thetas)\n",
    "\n",
    "if not tick:\n",
    "    eps_hat = mean(thetas[:120, 4])\n",
    "    f_hat = mean(thetas[:120, 5])\n",
    "    \n",
    "print 'eps_hat={}, f_hat={}'.format(eps_hat, f_hat)\n",
    "\n",
    "N = 130\n",
    "y_mins_norm = thetas[:N, 3].max()*y_mins/y_mins.max()\n",
    "\n",
    "plot(thetas[:N, 0], thetas[:N, 3], '.')\n",
    "plot(y_mins_norm[:N], '.')\n",
    "\n",
    "diff = abs(thetas[:N, 3] - y_mins_norm[:N])\n",
    "print 'avg diff is', sum(diff / N)\n",
    "print 'std diff is', std(diff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 653,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00926814497731\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff06d11f590>]"
      ]
     },
     "execution_count": 653,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEACAYAAABcXmojAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVPX+x/EXgopr4oILKKiY+4YLelXEPSs1r5kZtliZ\nlbbc7Ob1/lqsrrfs6s20Mk1Kr5pbuVduJSmKO26pJeaCCG6IgoICM78/zqhIKIsMZ2Z4Px+PeXDm\nnDnMx9Np3nyXcwZERERERERERERERERERERERERERETs7j7gEHAYGJ3N9lBgD7AX2AQ0y7TtmG19\nFLDNrlWKiIjduQPRgD9QHNgNNMzymvbAPbbl+4AtmbYdBSrat0QRESkoxXLY3hYjFI4BacB8oF+W\n10QCF23LWwHfLNvd7q5EEREpLDmFgg8Qk+n5Sdu623kG+CHTcyuwDtgBDMtPgSIiUng8cthuzcPv\n6gI8DXTItK4DEAdUAdZijE1szEuBIiJSeHIKhVigZqbnNTFaC1k1A77EGFO4kGl9nO3nWWAJRnfU\nLaFQt25d65EjR/JQsoiIAEeAgMJ+Uw/bG/sDJch+oLkWxrhDuyzrSwPlbMtlMGYm9czmPaxScN55\n5x2zS3ApOp4FR8eyYJG3npxcy6mlkA6MBFZjzEQKAw4Cw23bpwFvA17AVNu6NIwWQTVgcab3mQus\nKajCRUSk4OUUCgA/2h6ZTcu0/KztkdUfQIt81iUiIibIafaROJmQkBCzS3ApOp4FR8fSOTjCNQS2\n7jEREcktNzc3sMNnuFoKIiJyg0JBRERuUCiIiMgNCgUREblBoSAiIjcoFERE5AaFgoiI3KBQEBGR\nGxQKIiJyg0JBRMQBZGRA//5w/ry5dSgUREQcwLJlEBcHFU3+VnuFgoiIyaxWGD8e3ngD3Ey+I51C\nQUTEZBs3QmIi9OtndiUKBRER040fD6+/Du7uZleiW2eLiJhq3z7o2ROOHgVPz9zvp1tni4i4oAkT\n4OWX8xYI9qSWgoiISWJioHlzOHIEvLzytq9aCiIiLmbSJBg61AiElLQUJm+dTLol3dSaFAoiIia4\ncAFmzoS//Q02ndhE06lNiTgRQfK1ZFPr8jD13UVEiqgvvoA+faCGj4Vunz3N+13eZ3DTwWaXpVAQ\nESlsKSkweTKsXQsrf19JuRLleLTJo2aXBaj7SESk0M2YAUFB0KQJTIycyKj2o64PHJvOEarQ7CMR\nKTKuXoWAAFiyBDKqbWXgooEcefkIxd2L5+n32Gv2kbqPREQK0axZRguhVSsrXWaN5s3gN/McCPak\nloKISCFJS4P69WHOHEiovJLR60az5/k9eBTL+9/naimIiDi5b74Bf39oHXSNFl+8wX96/CdfgWBP\naimIiBSCjAxo1MiYihpR7H22xm5lxeAV+R5gVktBRMSJLVoElSuDd+MDTJ41majhUQ4z4ygzR6hI\nLQURcWkWCzRrBu98eIF3T3Ti5aCXea7Vc3f1O3XvIxERJ7VsGZQsk8LkM33pWbcnwwKHmV3Sbaml\nICJiR1YrtG4NjYd/SGK5zSx9dCnF3O7+73GNKYiIOKHlyyHNmsLapE9Y029NgQSCPSkURETsxGKB\nt96CDq/O4mSp1jSt2tTsknKkUBARsZNvv4USpVNYnfQRs7vPNrucXMlNO+Y+4BBwGBidzfZQYA+w\nF9gENMvDviIiLikjA8aOhWpPjiLIN4gOtTqYXVKu5NRScAc+BboDscB2YDlwMNNr/gCCgYsYITAd\naJfLfUVEXNK8eWCpv5gD11YR9UCU2eXkWk4thbZANHAMSAPmA/2yvCYSIxAAtgK+edhXRMTlpKXB\nmCk7iW8znPkPz+cez3vMLinXcgoFHyAm0/OTtnW38wzwQz73FRFxCZO+iuFMt7583X86bX3aml1O\nnuTUfZSXCwi6AE8D1zvOcr3v2LFjbyyHhIQQEhKSh7cVEXEc167B+z//h4d7h9K/Yf8C+73h4eGE\nh4cX2O+7nZwufGgHjMUYKwAYA1iA8Vle1wxYbHtddB731cVrIuIyPv4siTfi/Tj69734lvfNeYd8\nMus2FzuAeoA/UAIYhDFYnFktjEAYws1AyO2+IiIuIyUF3l02i2DfbnYNBHvKqfsoHRgJrMaYTRSG\nMXtouG37NOBtwAuYaluXhjHIfLt9RURc0qTJ6WS0+pS3e00zu5R8072PREQKQEIC1Bz8Ic0eWsvm\n59fZ/bbYuveRiIgDe238XixBE5k/eIdDfk9Cbjn2nZlERJxAxK9HmX2tPx92nYhfBT+zy7krCgUR\nkbtwJOEIveYH06PsKF4JecLscu6aI7RxNKYgIk6r0xf3E7W0M7ELRnNPIV64rDEFEREHszp6NVHH\noxnba2mhBoI9qaUgIpIP6ZZ06k1szpUV/+bE2n6ULFm476+WgoiIA5m2YzoJMVX57Lm+hR4I9qSW\ngohIHiWmJuI3oT7V163hwM/NKWbClB2zbnMhIiJZvPvzOCwH+jL9PXMCwZ7UfSQikgfxyfFM2x5G\nF8t+goPNrqbgKRRERPLg7dUfYdn9OFMm1jC7FLtQKIiI5FJ8cjyz9szk2fr7qVPH7GrsQ6EgIpIL\nRy8cpcfXfSke9RIfzHXNVgJooFlEJEfRCdF0+KoD1h3P8XHfsZQvb3ZF9qMpqSIid5B0NYn2Ye1p\n7zGC7Z+/wM6d4O5udlX2m5KqUBARuYOhy4ZizXDnl9e/5KswN7p0Mbsig65oFhEpZDEXY1h2aBnD\nU48R2NJxAsGeFAoiIrcxdcdU+vo/zpcjyrNrl9nVFA6FgohINlLSUpixawZNtm1i1CioVcvsigqH\nQkFEJBufb/8cP482xO6tx2tzzK6m8CgURESy2Ba7jfGbxlNyTiRhU3Cpu6DmRKEgIpLJpauXGPTt\nIEIuT8V6b1169jS7osKlKakiIpmM+H4E5xKv8tMrM4iKgpo1za4oe5qSKiJiZ5tObGLpb0tpsmE/\nr7/uuIFgT7rNhYgI8MeFPwhdHMqQilM4fsiL114zuyJzqKUgIkXekYQjdJnVhb+1/icTB/+Vb76B\nEiXMrsocGlMQkSIt3ZJO+7D2hDYN5fDsV0lLg+nTza4qZxpTEBGxgw8jPsTL04s2Ga/w0RL49Vez\nKzKXQkFEiqz9Z/bzydZP2DJ0F/1C3Jg0Cby8zK7KXOo+EpEiyWq1EjwzmMeaPEbCmheIjIQVK8DN\nET4Vc0HdRyIiBWjWnlmkpqcSUu45On0MO3c6TyDYk0JBRIqc44nHeWPtG/wYuooXQ915803w8zO7\nKsfgCLmo7iMRKTRpGWkEzwxmQMMBlNv3OmFhEBnpGN+mlhfqPhIRKQBvrX8LL08vBtR4jbb94Jdf\nnC8Q7EmhICJFxqroVczdN5edw3YR2r8Yo0ZBo0ZmV+VYcnObi/uAQ8BhYHQ22xsAkUAqMCrLtmPA\nXiAK2JbvKkVE7tKppFMMXTaUOf3nsGRuFS5ehNdfN7sqx5NTS8Ed+BToDsQC24HlwMFMrzkPvAQ8\nlM3+ViAESLjbQkVE8ivDkkHo4lBebP0i/m6defhNo9vIQ30lf5JTS6EtEI3xF38aMB/ol+U1Z4Ed\ntu3ZcYTBbBEpwsZtHEcxt2KM6fhPnn0WdRvdQU6h4APEZHp+0rYut6zAOozQGJa30kRE7t7+M/uZ\nsm0Ks/vPJmyGu7qNcpBT4+lu54p2AOKAKsBajLGJjXf5O0VEcsVitTB85XDeC3mPq+dq8OabEB6u\nbqM7yenQxAKZv2aiJkZrIbfibD/PAkswuqP+FApjx469sRwSEkJISEge3kJEJHvTdkwjw5LBsy2H\n0yUE/vEPaNzY7KryJzw8nPDwcLu/T079/R7Ab0A34BTGDKLB3DrQfN1YIAmYaHteGmOgOgkoA6wB\n3rX9zEwXr4lIgdt3eh9d/9eVDU9tYMmXDVm3Dtatg2Iu8tViZl28lg6MBFZjfMCHYQTCcNv2aUA1\njFlJ5QEL8ArQCPAGFmd6n7n8ORBERArclbQrPPLtI0zsOZHLxxsyaZJxbyNXCQR7coSZQWopiEiB\nenXVq5y9cpbpveYSGAjvvQeDBpldVcHSbS5ERHJh4/GNLDqwiH0v7OPvo6BtW9cLBHtSKIiIyzh5\n6SSPL3mcz+//nMifK/LDD7Bnj9lVOReFgoi4hDOXz9D9f90Z2XYk7bz60XIYzJ8P99xjdmXORWMK\nIuIS+i/oz70V7+WDbuO57z6j2+hf/zK7KvvRmIKIyG2sObKGfaf3MW/APMaPh9RUyHT5k+SBWgoi\n4tSuZVyj+RfN+aj7R3id7cPDD8OOHeDra3Zl9qWWgohINt78+U3qVaxH+0oPEtgDwsJcPxDsSaEg\nIk5rVfQq5u2fx67nonh6sBsDB8IDD5hdlXNTKIiIUzqWeIyhy4Yyb8A85n5Zmfh4+PZbs6tyfhpT\nEBGnc/naZTp81YGnWjxFe7dX6dMHtmyBOnXMrqzwaExBRASwWq08vfxpWlRrweDar9CmDXz5ZdEK\nBHtSKIiIU/kw4kOOJR7jpyG/0Pd+N4YMgX5Zvw9S8k3dRyLiNNYfXc+QJUPYPmw7k8fVYOdOWLUK\n3N3Nrqzw2av7SKEgIk4h+VoyzaY2Y0rvKVzd/wB/+5txPUKVKmZXZg6NKYhIkTZm3Rg6+XUiwPoA\nHYfD998X3UCwJ4WCiDi8hb8uZMXvK9gQuoveITBunHFvIyl46j4SEYd2/Ws1Vz22hn+NbEmlSsZs\nIzdH+PQykbqPRKTISb6WzMOLHubjXh+z9IuWnD1r3A67qAeCPTnCoVVLQUSy9eTSJ3F3c6dX6leM\nHg3btoG3t9lVOQa1FESkSPlkyyfsOLWDaa220f9BWLtWgVAYFAoi4nC+2fcNEyInsPjBCAZ0K8O0\nadCihdlVFQ3qPhIRh3Lg7AGCvw5m9eBwRj7ShN694e23za7K8ejiNRFxeemWdP4S9heeavE0G//7\nPBaLBpZvR2MKIuLy3vr5LSp4ViBmyXBOnICfflIgFDaFgog4hHEbxrHst2U8XSycaYvciIwET0+z\nqyp6FAoiYrqvo75m5p6ZvOu3gdee82bjRqhc2eyqiiZHaJhpTEGkCNsTv4fus7szrd0vPD+gEUuW\nQIcOZlfl+Ow1plCsoH+hiEhuXUy9yMBFA3mrzSReeawRn36qQDCbWgoiYgqr1crARQMp516FrW9N\n5ZlnYNQos6tyHpp9JCIuZfym8Rw5f5Tis+fQp48CwVGopSAihW7chnF8vXsmPut+pm7lmoSFaepp\nXqmlICIu4bNtnzFn7xyabt+Axa0606crEByJI/ynUEtBpIjYeHwjAxYOoHdsJH/srMuaNVCqlNlV\nOSe1FETEqZ2/cp5Hv3uU7smz2L2+Lr/8okBwRGopiEihGPzdYGIOVuP83I8JD4eqVc2uyLmppSAi\nTmvhrwv56cAuys7ZzcafFQiOLDcXr90HHAIOA6Oz2d4AiARSgayTynLaV0Rc3IbjG3jmuxG4L53L\n+jWl8PExuyK5k5yaHu7Ab0B3IBbYDgwGDmZ6TRXAD3gIuABMzMO+oO4jEZe1J34Pnb7sQcmV84ic\n242AALMrch1m3eaiLRANHAPSgPlAvyyvOQvssG3P674i4qLOXj5L97B+uK+ZwoaZCgRnkVMo+AAx\nmZ6ftK3LjbvZV0ScWNLVJDp9+jAp2wcT8cUgGjY0uyLJrZxC4W76ddQnJFIExVyModHEjsRE1WfL\nuH/RuLHZFUle5DT7KBaomel5TYy/+HMj1/uOHTv2xnJISAghISG5fAsRcSRxSXEETulMWuTz7Jz4\ndxo0cIRZ764hPDyc8PBwu79PTv/FPDAGi7sBp4BtZD9YDDAWSOLmQHNu99VAs4gLSExNpPGEYC5v\nHcSOj/9PYwh2ZtZ1CunASGA1xmyiMIwP9eG27dOAahgzi8oDFuAVoBGQfJt9RcTFXEm7QuCEPiTt\n7UrUpH9St67ZFUl+OULbTi0FESeWnpFBo/f6E3/8Hva+Owt/P313V2HQFc0i4nAsFmj/z7HEXU7i\n0Aff4VNdgeDsFAoiki9padBz5Ar2VpzJ/n/uwKd6cbNLkgKgWBeRPEtNhc5D17G5yjOsfmYR9Wro\nZkauQqEgInmSkACtHv2BXbUf48eh3xES0M7skqQAKRREJNeOH4fGT0znePNn+HnYcrrW7WR2SVLA\nNKYgIrkSFWUl+N038QxayO4RGwmoqAsRXJFCQURytHq1lYfCnqfmX/awacRmqpSpYnZJYicKBRG5\no6+/hpfm/4eaPbex66WNlC1R1uySxI4UCiKSrYwMGP0PK7MOfEbZbpP5edgWBUIRoFAQkT9JToZH\nh6Sw1TeUGg8c59tB6/Et72t2WVIIFAoicosTJ+DBvulc7DWYLq08mTMgkhLuJcwuSwqJpqSKyA2R\nkRDUzkqpgcNp2DSVOQP+p0AoYhQKIoLVClOnQr9+0OHtMbhV/ZXvHvlOgVAEqftIpIhLSYEXX4Rt\nO9J4cMo/iDz3AxFPRFCmRBmzSxMTqKUgUoQdPw6dOkHC1TPc80pX4tMPEjE0gkqlK5ldmphEoSBS\nRK1bB0FB0GXQAfYFtaNbnRBWPrZSgVDE6Ut2RIoYqxUmTIAJky/Tdez7rDn/JR/3+pgnmj9hdmmS\nB/qSHRG5a+fOwdChEH/aSqv3nyS1VAb7HtlHjXI1zC5NHIS6j0SKiA0boGVLaNAAHv9kKqevHWX+\ngPkKBLmFWgoiLi4jA8aNM6acfvUVeDZYz6Bvx7Lp6U2U9ChpdnniYBQKIi7s1CkIDQU3N9i5E46m\nb6L/gkEsHLiQepXqmV2eOCB1H4m4qB9+gFatoGtXePerzQxZ25W/LvwrMx+aSYh/iNnliYNSS0HE\nxSQnw+uvw6pVsGCBlV0lPuHhbz9gYs+JPNL4EV2lLHekUBBxIRER8OST0LkzLN9wjFG/DONi6kW2\nPLOF2l61zS5PnIC6j0RcwNWrMHo0DBwIEyZaaPPCVLotaEP32t3Z/MxmBYLkmloKIk5u9254/HGo\nVw82bk/k1Q1DOLfnHBue2kDDKg3NLk+cjFoKIk7q2jV47z3o0QPeeAPem/orvZe0oY5XHTYO3ahA\nkHxRKIg4ochICAyE7dth1y7wDFxEl/+F8GanN5ncezLF3YubXaI4KXUfiTiRpCT4v/+DRYtg0iRo\n1f0II9b8jd/P/86PoT/SukZrs0sUJ6eWgoiT+P57aNLEmHIavv0MqzyH0i4siHa+7djz/B4FghQI\ntRREHFx8PLz2GmzdCmFhkO6/ii4LnmZIsyFEvxxNBc8KZpcoLkQtBREHlZYGH39stA5q1oTVm2OZ\nfuERnl/5PHP+OoePenykQJACp1AQcUDh4cYdTX/80bggLeipxXSYHUj9SvU5MOIAXWt3NbtEcVHq\nPhJxILGxxi0qNm82WgkNOh7kjZ9Gc/DcQZY/upwg3yCzSxQXp5aCiANITYXx46F5c6hbF6b8sIap\nl3vQeVYwwX7B7H9hvwJBCoW+jlPERFYrLFgAY8ZAk2ZpPPTaTyyOncLh84d5p/M7DGg0AE8PT7PL\nFAdkr6/jzM0vvA+YBLgDM4Dx2bxmMtAbuAI8BUTZ1h8DLgEZQBrQNpt9FQpSJEVEwKhRxpfg9P/H\ncqbFjMSnvA+PN3ucZwOf1d1M5Y7M+o5md+BToDsQC2wHlgMHM73mfiAAqAcEAVOBdrZtViAESCiw\nikWcXHS0cfO67dth9Hun+aXMy8yM28Xs/rPp7N/Z7PKkiMtpTKEtEI3xF38aMB/ol+U1fYFZtuWt\nQAWgaqbtjtBFJWK6+Hh4+WVo1w4CW2Xw9wXTeO9cM2pX8Gfv83sVCOIQcgoFHyAm0/OTtnW5fY0V\nWAfsAIblv0wR53X+vNEyaNwY8Ejh7cUzWeDVkoWH5rAqdBXje4ynVPFSZpcpAuTcfZTbzv7btQY6\nAqeAKsBa4BCwMeuLxo4de2M5JCSEkJCQXL6tiOO6dMmYVjplCnQcvIlGH77D1+e20jGmIx92/5De\nAb2v9wuL5Cg8PJzw8HC7v09OZ2Q7YCzGYDPAGMDCrYPNXwDhGF1LYHzwdwZOZ/ld7wDJwMQs6zXQ\nLC7l8mX47DOYMAG69b6MV99/s+TYV0zoMYE+9ftQvmR5s0sUF2Cvgeacuo92YAwg+wMlgEEYA82Z\nLQeesC23AxIxAqE0UM62vgzQE9h31xWLOKiLF+Hf/4Y6dWBzVALPfjmZDc3vJcF6hKjhUYQ2C1Ug\niMPLqfsoHRgJrMaYiRSGMfNouG37NOAHjBlI0cBlYKhtWzVgcab3mQusKajCRRzFuXPwyScwdSp0\n6HeQTv+ZyLpT31I6/X6WDFpCW5/sZmKLOCZH6NBU95E4pbg4mDgRwr6y0CF0PSnNJ7E/cRsj2ozg\nxTYvUrl0ZbNLFBdm5sVr9qZQEKdy4AD897/w7ZqTNBw8k1NVv6JC6XK80PoFnmz+pGYSSaFQKIiY\nyGqF9euNweMdB89Q56n3+b3kNwxq8gjPBD5Dq+qtNJNICpVZVzSLFGlpabBwoREGVywXafDkJ2R0\n/oSgZkNYEXyIKmWqmF2iSIFSKIhkIz4epk+HadMteLfeRLVhc9mStIA2/g+wvct26njVMbtEEbtQ\nKIjYWK3G9xh8+ims2L+eKg9MIWXERjIqVKdj40GEtfyVGuVqmF2miF05QieoxhTEVJcvw7x5Rhgk\nlIyiUr8POO+5jXe7vEOPuj3wLe9rdokif6KBZpECZLUadymdMQPmb9xGpa6z8ai9iVT3s4xsM5KX\ngl6idPHSZpcpclsKBZECcO4czJkDM8KsXCgfTqle40gp/Tsjgp4n2C+YIJ8girsXN7tMkRwpFETy\nyWKBdeuMIPhxx37qPLCUK36LcStxmTEdxxDaLFRfaCNOR6Egkkd798LcuTBr7Q4sjedhvXcppcpk\nMKBRf/o16EenWp1wL+Zudpki+aJQEMmFEyfgm29g9rwUTpf/kbJdPuda2cMMa/MU/Rv2p3nV5rrI\nTFyCQkHkNhISYNEimLnwHHuZS8W2P3KhTCTtarXhieZPMLjJYI0TiMtRKIhkcu4cLFsGC79LIyLm\nFyr1mMF5r1X0a9iHgY3/Soh/CF6lvMwuU8RuFApS5MXHw5IlMGfFCXYl/UDFtqtJrLCe+lUCeDrw\nKUKbhioIpMhQKEiRdPw4LF5iYebabfzGCko1X4ml9Cnur38fD9TvRY86PahatqrZZYoUOoWCFAkZ\nGbB1KyxemcyiXWs5U2EF3Ps93mUq80iLPjzUsA/tfNtp1pAUeQoFcVmXLsHq1TD/xxOsPraCYvVX\ncrXqJppXCuKx1n3o2+BB3YBOJAuFgrgMi8W4huDH1el8F7mdvanf49lsBZaysfSqcz+PtuxDr4Be\n+j5jkTtQKIhTi4uD1WsyWLhhDxtP/ozVbz1p1SOoXtqP/k17M6BJH9r7tle3kEguKRTEqVy8CBER\nsHx9HCsPL+Ns+dXg/wuVPKvSvW5X+jXrQme/zvqSGpF8UiiIQ0tMhPUb0lgQsZWI45s5XWwnxf12\nYi11ns7V72dI0AN0q9OF6uWqm12qiEtQKIhDSUiA79efZ8mWXWw9tofTpcJx89tIFY+6dKjZifub\ntyaoVivqV6qvLiERO1AoiGksFtix7xIrNh9m84Gj7E7YRGKF9bhV/INaxQMJ9G3KgFbB9KjXhcql\nK5tdrkiRoFCQQnPpEqyJOMey7TvYFBvOCff1WCodoEJGAL7laxFcpx2DgkJoV7O17ikkYhKFgthF\nair8tO0UK3fuYuuJXURf3kVyuV0UK3WJGsVa0r5GMIPaduH+Zu3w9PA0u1wRsVEoyF07dymZ+RGR\nrDmwlcOnTxGbfIykslG4F0+nmiWQxpUCCakfyENBgdT3rk0xt2Jmlywit6FQkFyzWq38ejyeRVs3\nEf5HBDEX4om/doSUMgcpfakldUv8hYY1atG6ni8PBQUS4O2r7xgQcTIKBfmTtIw0Dp4+wk97DhEZ\nfYgDpw9xMvUQl0ocAosHFS4H0cCzM419a9H63po80rEVXuXUBSTiChQKRVS6JZ3jicc5eCaaLb9F\nszsmmujz0cRdO0yS+zG45EuZlAb4lGxAQ+8GtA9oSM/A+jSvVxn98S/iuhQKLiwxNZHjicc5lniM\ng/HH2Hk0mt/ORBNzJZpLnKDYlepknA2gfHpdfEvX497KAbSqHUC3lgG0aOKJp/74FylyFApOKsOS\nQVxyHCcuniDmYgwnLp7gj/MxHIo7ztHEY5xOPU6GJYMSV/xJP+9Hxnl/qpUIoG7FAFrUDKB9g9o0\nbVSSgAAoUcLsf42IOAqFggPKsGQQnxzPyUsnOXnpJLFJsZy8dJKYizEcTYjh2IUTnEuNpzRVKHm1\nJtbEmqTE1yLtXE2qefpRp5IfjX38aVLXi4AANxo0AB8fKKZJPyKSA4VCIUpNT+Xs5bOcuXyGM5fP\ncPryaeKS4ohNiuVU0iliL8Vy8uIpzlw5TVn3SpS1+FI8xRdLoi+pZ3xJPOGLe3JN/CrUoqFvDe6t\nW4KAAG48qlfXB7+I3B2Fwl3IsGSQkJJw40P+T48rtz5PTU+lYglvyrl745nhjcfVKlgv1SD1XA0u\nnfTh3NEalM6ogX+l6tTxK4GfH/j7c8vPChXs+k8SkSLOXqHgUdC/sDBYrVaSryXn+kM+ISWBCp4V\n8C7tjVcJb0rjTck0b9xSvEm/GMjV896knfbGGuMNR70pllyesj5u1PABX1+jS8evNviF3PzQL1vW\n7KMgIlLwcpMy9wGTAHdgBjA+m9dMBnoDV4CngKg87HvblsKBsweYtXvWnz7kz1w+g0cxD7zLeN/4\noC+DNyUzvHFP9caa7E1aojcp57xJivMmIbYSZ+I9OHsWypc3PuSvf9hnt1yhAprOKSIOzayWgjvw\nKdAdiAW2A8uBg5lecz8QANQDgoCpQLtc7ntHFxOLcTbGC/eU+lS/7E2lC97UPu9N8mlvzsWVJj4e\ndp2GUqWgalWoVs14XF+uWufWdd7erj+DJzw8nJCQELPLcBk6ngVHx9I55BQKbYFo4Jjt+XygH7d+\nsPcFZtlMg0bbAAADmElEQVSWtwIVgGpA7Vzse0dpcQ04Of8feHtDlSq2Rx3jw/36B33VqkYoiEH/\n4xUsHc+Co2PpHHIKBR8gJtPzkxitgZxe4wPUyMW+dxQcDGvW5GUPERG5GzlNjMzttCD1wIuIFAHt\ngFWZno8BRmd5zRfAo5meHwKq5nJfMLqYrHrooYceeuTpEY0JPIAjgD9QAtgNNMzymvuBH2zL7YAt\nedhXREScTG/gN4xUGmNbN9z2uO5T2/Y9QGAO+4qIiIiISFHxFXAa2JdpXUVgLfA7sAZjqup1Y4DD\nGOMPPTOtb2X7HYeBTzKtLwkssK3fAvgVbPkOJ7vjORZjBleU7dE70zYdz9urCawHfgX2Ay/b1uv8\nzJ/bHc+x6PzMD0+Mqfy7gQPAB7b1Tn9+dgJacuuH2EfAG7bl0cCHtuVGGAegOMZ4QzQ3Zy9tw7g2\nAoxxivtsyy8Cn9uWB2Fc8+DKsjue7wCvZfNaHc87qwa0sC2XxejObIjOz/y63fHU+Zl/pW0/PTA+\ntDviIuenP7d+iF2fhQTGiXTItpx1FtIqjAHq6tx6YdujGDObrr/m+jUOHsDZgiragfnz51AYlc3r\ndDzzZinGVfY6PwvG9eOp8/Pulca480NjTDw/7XkD56oYXSDYfl7/B9bAaGZel/lit8zrY23r4dYL\n5NKBixjNq6LmJYzB/DBuNid1PHPPH6MFthWdnwXBH+N4Xp9xqPMzf4ph/PV/mptdc6adn4V1V//r\n82ol/6Zi3DqkBRAHTDS3HKdTFvgOeAVIyrJN52felQW+xTieyej8vBsWjOPmCwQDXbJsL9Tz056h\ncBqj2QNG0+aMbTkWY7DqOl+MhIu1LWddf32fWrZlD+AeIKHgS3ZoZ7h5cszgZt+hjmfOimMEwmyM\n7g7Q+Xk3rh/POdw8njo/795F4HuMAWPTzk97hsJy4Enb8pPcPHmWY/R3lcD4y6IexgBJPHAJo+/L\nDXgcWJbN73oY+MmOdTuq6pmW+3NzvEHH887cMLozDmDcxv06nZ/5c7vjqfMzfypzs6utFNADY/aW\n05+f84BTwDWMvquhGH1W68h+StU/MUbNDwG9Mq2/PqUqGuM7Gq4rCSzk5pQqfzv8GxxJ1uP5NPA/\nYC9Gn+1SbvYxgo7nnXTEaJ7v5uZ0yfvQ+Zlf2R3P3uj8zK+mwC6M47kX+Lttvc5PERERERERERER\nERERERERERERERERERERERFxHf8POzLezqci4B8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff06d8e3a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot(thetas[:N, 2], thetas[:N, 3])\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def small_control(x, a, c, d):\n",
    "    return a*np.exp(c*x) + d\n",
    "\n",
    "popt, pcov = curve_fit(small_control, C[4:,0], y_mins_norm[4:], p0=(.001, 1e-3, 0))\n",
    "#print popt, pcov\n",
    "\n",
    "t = linspace(1000, 30000, 100)\n",
    "\n",
    "#plot(C[4:], y_mins_norm[4:])\n",
    "#plot(t, small_control(t, *popt))\n",
    "\n",
    "y_hat = small_control(C[4:], *popt)[:,0]\n",
    "\n",
    "print max(abs(y_hat - y_mins_norm[4:])) \n",
    "plot(C[4:], y_hat)\n",
    "plot(C[4:], y_mins_norm[4:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 693,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.00136871  0.00017554 -0.0057529 ]\n",
      "0.22262857715\n",
      "0.939737022742\n"
     ]
    }
   ],
   "source": [
    "def GetControl (phi, popt):\n",
    "    a, c, d = popt\n",
    "    inflection = y_mins_norm[4:].max()\n",
    "    fastest = y_mins_norm.max()\n",
    "    if phi < inflection:\n",
    "        return log((phi - d)/a)/c\n",
    "    else:\n",
    "        exp_control = log((inflection - d)/a)/c\n",
    "        max_control = 32000\n",
    "        coin = random.rand() * (fastest - inflection)\n",
    "        if coin + inflection < phi:\n",
    "            return max_control\n",
    "        else:\n",
    "            return exp_control\n",
    "\n",
    "print popt\n",
    "print y_mins_norm[4:].max()\n",
    "print y_mins_norm.max()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 698,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0043841921696713944"
      ]
     },
     "execution_count": 698,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "small_control(0, a, c, d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 697,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stored 'M' (list)\n",
      "Stored 'C' (ndarray)\n",
      "Stored 'S' (ndarray)\n"
     ]
    }
   ],
   "source": [
    "%store M\n",
    "%store C\n",
    "%store S"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}