Fitting circles.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:7fa5cc3c072f9d4bf24ff9db314d8b791f4d6e0ee3906d18f98a790a918a2be6"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "import matplotlib.pyplot as plt\n",
      "from matplotlib.patches import Circle\n",
      "from matplotlib.collections import PatchCollection\n",
      "\n",
      "from scipy import optimize"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 49
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Your objective is to measure the radius of a circle as accurately as possible. You have `N` opportunities to make a measurement, where on the circumference of the circle would you place them?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "O = (0,0)\n",
      "layers = [2, 3, 4,5]\n",
      "curvature = 3.456*2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = plt.subplots(figsize=(8,8))\n",
      "\n",
      "intersect_x = []\n",
      "intersect_y = []\n",
      "\n",
      "for lay in layers:\n",
      "    circ = Circle(O, lay, edgecolor=\"black\", fill=False)\n",
      "    ax.add_patch(circ)\n",
      "    \n",
      "    # Find the \"upper\" points at which each\n",
      "    # layer intersects with the track\n",
      "    d = curvature\n",
      "    a = 1./d * np.sqrt((-lay)*(-2*d+lay)*(lay)*(2*d+lay))\n",
      "    y = a/2.\n",
      "    x = (d**2 - d**2 + lay**2) / (2*d)\n",
      "    intersect_x.append(x)\n",
      "    intersect_y.append(y)\n",
      "    \n",
      "track = Circle((curvature, 0), curvature, edgecolor=\"red\", fill=False)\n",
      "ax.add_patch(track)\n",
      "\n",
      "ax.plot(0,0, 'ko')\n",
      "\n",
      "ax.plot(intersect_x, intersect_y, 'rs')\n",
      "\n",
      "ax.set_ylim((-10,10))\n",
      "ax.set_xlim((-7,13))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 69,
       "text": [
        "(-7, 13)"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAHfCAYAAACBE6uXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xlczdn/B/BXm4xd0i4J2bc2ZSwpW6VkGQxmJswMEoYx\nwzAzYswMY/mOrSyDmLGNdUqEkGUYLdIi0hBabtkqhere7vn98fllGVvpc+793Ho/H4/7QPJ+H1uv\ne87nfM5HizHGQAghhBC10lb3AAghhBBCgUwIIYRIAgUyIYQQIgEUyIQQQogEUCATQgghEkCBTAgh\nhEhApQN53LhxMDY2RocOHZ5+7MGDB+jbty9sbGzQr18/5OXlVbYNIYQQUqVVOpDHjh2L8PDwFz62\naNEi9O3bF9euXYObmxsWLVpU2TaEEEJIlaYlxsEgN2/ehJeXFxITEwEArVu3xqlTp2BsbIzs7Gy4\nuLjg6tWrlR4sIYQQUlVxuYack5MDY2NjAICxsTFycnJ4tCGEEEKqDF3eDbS0tKClpfXSxzt37oz4\n+Hje7QkhhBDJ6NSpEy5duvTKn+MyQy5bqgYAmUwGIyOjlz4nPj4ejDGNes2bN0/tY6AX/f1Vxxf9\n3Wn2i/7+nr3eNBHlEsje3t7YsmULAGDLli3w8fHh0YYQQgipMiodyB9++CG6deuGlJQUNGnSBJs3\nb8bs2bNx7Ngx2NjY4MSJE5g9e7YYYyWEEEKqrEpfQ96xY8crPx4REVHZ0pLj4uKi7iGQSqC/P81F\nf3eajf7+ykeU257eqbGWFtTUmhBCCFGLN2UfHZ1JCCGESAAFMiGEECIBFMiEEEKIBFAgE0IIIRJA\ngUwIIYRIAAUyIYQQIgEUyIQQQogEUCATQgghEkCBTAghhEgABTIhhBAiARTIhBBCiARQIBNCCCES\nQIFMCCGESAAFMiGEECIBFMiEEEKIBFAgE0IIIRJAgUwIIYRIAAUyIYQQIgEUyIQQQogEUCATQggh\nEkCBTAghhEgABTIhhBAiARTIhBBCiARQIBNCCCESQIFMCCGESAAFMiGEECIBFMiEEEKIBFAgE0II\nIRJAgUwIIYRIAAUyIYQQIgEUyIQQQogEUCATQgghEkCBTAghhEgABTIhhBAiARTIhBBCiARQIBNC\nCCESQIFMCCGESAAFMiGEECIBFMiEEEKIBFAgE0IIIRJAgUwIIYRIAAUyIYQQIgEUyIQQQogEUCAT\nQgghEkCBTAghhEgABTIhhBAiARTIhBBCiARQIBNCCCESQIFMCCGESAAFMiGEECIBFMiEEEKIBFAg\nE0IIIRJAgUwIIYRIAAUyIYQQIgEUyIQQQogEUCATQgghEkCBTAghhEgABTIhhBAiARTIhBBCiARQ\nIBNCCCESQIFMCCGESAAFMiGEECIBuuoeACGEaByFAigoEF5FRYBSKbxKSwEtLUBb+9nrvfeAOnWA\nunUBXfqSS16P/nUQQkhJCZCVBWRmAhkZz77NyABkMiA/HygsfBbCcrkQsHXrAjVrAjo6zwIYeDGg\nnzwRfk1hIVCjxrNwrlsXqFcPMDMDLCyEl7n5s2/NzAA9PfX+uRCV0mKMMbU01tKCmloTQqqjx4+B\n+HggNhZITn4xdHNzAWPjF0Px+WCsX/9ZiNapI8x6tbQq1p+xZ+FcFtB5eS++EXj+zcCdO0CjRi+O\np21bwM4O6NRJeCNANM6bso8CmRBS9TwfvmWvf/8F2rQRAq1DhxdD19hYmOVKiUIB5OQ8C+j0dCAx\nUfi9pKQANjbC76XsRSGtESiQCSFV19vC195e+LZ9e0BfX92jFUdRkRDOMTHPfs8U0hqBApkQUnXk\n5QHh4cDRo0Ig/fvvs6XcsldVCt/yKgvpsoCOiXkW0vb2QP/+wqtePXWPtFqjQCaEaLa0NCAkBAgN\nBaKigF69AA8PwNGxeoZveZWF9IULwKFDwNmzgLMz4OUlvJo2VfcIqx0KZEKIZlEqgehoIYRDQoQN\nTmUh0qcPULu2ukeomQoLhZWFkBAgLEy4fu7tLbxsbZ/tEifcUCATQqTv8WMgIkIIi4MHAUPDZ2Hh\n6EhhIbbSUuCff5696Xn4EBg4UPjzdnUVdpIT0VEgE0Kk6e5d4K+/hECIjBSudXp7CzPh5s3VPbrq\n5do14ZJASAhw6ZIQymVviBo1UvfoqgwKZEKIdDAGnDsHBAYK1zX79wcGDQIGDAAaNlT36AgA3L8P\nHD4svFk6dkz4+/HzE1YqKnr/NXkBBTIhRP0KC4Ft24QgfvJE+AL/yScUwlJ37x6weTMQFCT8Xfn5\nAR9+CNSqpe6RaSQKZEKI+iQnC1/Mt20DXFyEL+iurnRNWNMolcCRI8IbqnPngI8/BiZNEm6rIuX2\npuyj/xGEEPHJ5cDu3UDv3oCbmzCzio8H9u0TdklTGGsebW3A3V24zhwbKxw40r070LcvsH+/cLIY\nqRSaIRNCxJORAWzYILxathRmw4MHCw9VIFVPcTGwZ48wa759G/j8c+DTTwFTU3WPTLJohkwI4ev8\neWDoUOGM6Hv3hHtdT50CRoygMK7K9PWB0aOBv/8WZs4ZGcKpaSNGCCeFkQqhGTIh5N1dvgzMnQtc\nvAjMmiVcV6xbV92jIuqUnw8EBwO//AK8/z6wcCFdZ34OzZAJIeK6dQsYO1a4Rtyzp3AP6+TJFMZE\neFTltGlAaqpw+tf77wMTJghPrSJvRIFMCCm/e/eAGTOEL7Tm5sIX3Rkz6IlC5GW1agGzZwsPuGjQ\nAOjYUfhxbq66RyZZFMiEkLcrLBSWHlu3FjbyJCUJP65fX90jI1JnYAAsXizssn/wQFi+XrxYOCqV\nvIACmRDyeiUlwJo1wo7p5GTh7OM1a2gXLak4Cwtg/XrgzBnhwSE2NsKP5XJ1j0wyKJAJIS9TKoHt\n24E2bYQHPRw6JPy4RQt1j4xoutathVul9u8Hdu0SHp+5e7dwpGo1R7usCSEvOnYM+Oor4ZaWxYuF\n07UI4YEx4d/b7NmAjg6wZEmV//dGR2cSQt4uL0/YoHXyJLBsmXCgBz1IgKiCUinMkr/6CvD0FG6Z\nqqI79um2J0LImx0+LBzqUbMmkJAADBlCYUxUR1tbOEwkMVHYt9ChA3D8uLpHpXI0QyakOnt+Vrxx\no/DQB0LULTxcOIazCs6WaYZMCHnZf2fFFMZEKgYMqJazZZohE1Ld0KyYaJIqNlumGTIhRECzYqJp\nqtFsmWbIhFQHNCsmVcHhw8JseeBAjZ0t0wyZkOqMZsWkqnB3r9KzZZohE1JVyeXAzJlASAjNiknV\nUzZbHjUK+Okn4WARDUAHgxBS3dy9CwwfDrz3nnDkZYMG6h4RIeK7d0+4f1lPD9ixA2jYUN0jeita\nsiakOomLAxwcAGdnIDSUwphUXYaGwJEjwpnrDg7A5cvqHlGlUCATUpXs2AH06ydseNGgZTxC3pmu\nLvC//wHffy+cg71/v7pH9M5oyZqQqqC0FPjmG+EpOgcOCA+DJ6S6iYkRjn0dOxaYN084klNi6Boy\nIVXZgwfAhx8KobxrF9CokbpHRIj65OQAw4YJ15P/+AOoV0/dI3oBXUMmpKpKSgIcHYVnyoaHUxgT\nYmws3A7VpAnQtSuQkqLuEZUbBTIhmurQIaB3byAgQHhcoq6uukdEiDTUqAGsWQN8+SXQo4fG3K9M\nS9aEaKKtW4GvvxY2sDg7q3s0hEjX6dPCEvbq1cKtgGr2puzj+pbaysoK9erVg46ODvT09BAVFcWz\nHSFVH2PA0qXCu/+TJ4XbPQghr9ezJxARAXh4CNeXp0xR94hei2sga2lpITIyEgYGBjzbEFI9KJXC\nyVtHjwJnzwIWFuoeESGaoWNH4f9M//5AdjawcCGgpaXuUb2E+zVkWpYmRAQlJcBHHwFRUcCZMxTG\nhFSUlZUQyseOAZ99BigU6h7RS7heQ7a2tkb9+vWho6ODCRMm4LPPPnvWmK4hE5EUFBRAJpNBJpMh\nKysLubm5kMvlUCgUL7y0tbWhq6sLXV1d6OnpQVdXF7Vq1YKJiQnMzMxgamqKxo0bQ1tq9y4+fizc\nW6mvD+zcKRyHSQh5N4WFwjXlmjWF2wT19VXaXm3XkP/++2+Ympri7t276Nu3L1q3bo0ePXo8/fmA\ngICn33dxcYGLiwvP4RANVFJSgsuXLyM5ORlZWVkvBG/Z95VK5dNANTMzg4GBwdPAff6lUCjw+PHj\npwEtl8vx6NGjF2rm5+fDyMgIpqamT+uVfb9Zs2awtbVFI1XeWvToEeDtDZiaAsHBtJOakMqqU0d4\n4Mro0cDgwcC+fUI4cxIZGYnIyMhyfa7KdlnPnz8fderUwZdffik0phky+Y+y8I2NjUVMTAxiY2Nx\n+fJlNGvWDO3bt4e5ufkLwVv2bd26daEl0vWgkpISZGdnvxT8WVlZ+PfffxEXFwcDAwPY2dnB3t4e\ndnZ2sLOz4xPShYWApydgbQ389hsdg0mImORy4TJQbq5wup2KVp7UclLX48ePUVpairp16+LRo0fo\n168f5s2bh379+r11UKTqUyqVSExMRHR09EvhWxZy9vb26Ny5M2rXrq3u4T6lVCqRmpqK2NjYp6+4\nuDg0bNjwhXE7OTmhXmVOCCosFHaF2tgA69dL8ghAQjSeQgH4+gobvUJDVRLKagnktLQ0DB48GACg\nUCgwevRofPPNN+UaFKmanjx5ghMnTiAkJAShoaGoU6cOnJycns40pRa+5aVUKvHvv/8+fWMRExOD\nixcvolu3bvD29oaXlxcsLS3LX/DRI2Fm3Lw5sGEDhTEhPJWWCjPlBw+EmTLH5WuAzrImapSTk4Ow\nsDCEhobixIkT6Ny589OQsrGxUffwuCkoKMDRo0cREhKCsLAwNGnS5Onv29bW9vUbx548AQYOFI79\n27SJwpgQVVAogFGjhDfD+/Zx3ehFgUxUhjGGK1euICQkBCEhIUhOTka/fv3g7e0Nd3d31W6IkgiF\nQoHz588/XRkoKCiAl5cXvL29cX7HDuikpwufyJhwNrWODjBgAAK2blXvwAmpTuRyYMQIQE9PeIwp\npzfDFMiEu8LCQmzbtg2BgYHIzc2Ft7c3vL290atXL+ir+LYCqUtJSUFoaChCQkJQ49w5RJSWvvQ5\nAb16IaCcOzMJISIpKgIGDBAOElmxgsvhIRTIhJvk5GQEBQVh27ZtcHFxgZ+fH1xdXaV3L69EzXZy\nwqILF176OAUyIWqSlwf06gWMHCk8Y1xkarsPmVRNcrkcBw4cQGBgIK5evYrPPvsM8fHxaNKkibqH\npnFqvmYDyc2bN5GZmQlzc3MVj4iQaq5BA+DwYeD99wEjI2D8eJW1pkAm5ZaZmYn169djw4YNsLGx\ngZ+fH3x8fFCjRg11D01z3bnzyg+XlJSgQ4cOcHV1hZ+fH3r37i3avdaEkLcwMwOOHBFmyo0bC4fz\nqACtK5K3unHjBsaMGYMOHTrg/v37OHbsGCIjIzF8+HAK48qIjARSU1/5UzY2Nrh16xbc3NwwdepU\ntG/fHnv27KHLPISoio2NcKLXp58C58+rpCVdQyavlZOTg4ULF2LHjh2YOnUqvvjii8oddkGeSUsD\nnJ0R0LGj8OCI/7KyQkBwMABh5/qRI0cwZ84c6OrqYtGiRXB1dVXteAmprg4dEh5GERUFiHAJiTZ1\nkQp5+PAhli5dijVr1uCjjz7C3Llz0bhxY3UPq+p49Ei4PjV2LDBtWrl/mVKpxJ9//olvv/0W1tbW\n+Pnnn2FnZ8dxoIQQAMDPPwuHhpw6VemDQ96UfbRkTZ4qLi7G//73P7Rs2RK3bt1CbGwsfv31Vwpj\nMTEmbBLp1AmYOrVCv1RbWxsjR47ElStXMHjwYHh5eWHEiBFIfc2yNyFEJLNnA5aWwOTJwv9hTiiQ\nCUpLS7Flyxa0atUKJ0+eREREBLZs2QIrKyt1D63qWbIEuH4dWLv2ne9x1NPTw6RJk5CamopOnTrB\n2dkZkyZNgkwmE3mwhBAAwv/VzZuB6GhgzRpubSiQq7m///4bnTt3xoYNG/DHH38gJCQEHTp0UPew\nqqbwcODXX4H9+0U5xL527dqYM2cOUlJSULt2bbRv3x4//PAD5HK5CIMlhLygTh1h2XrhQmFDJgd0\nDbmaevz4Mb777jvs2LEDK1euxNChQ+m2Gp7+/Ve4brx3L9C9O5cWt2/fxsSJEyGTyRAcHIxOnTpx\n6UNItRYRITyM4sIFYRm7gugaMnlB2aw4KysLCQkJGDZsGIUxTwUFwKBBwPz53MIYACwtLREWFoap\nU6eib9++WLBgAc2WCRFbnz7AV18BPj7A48eilqYZcjXy/Kx49erVGDJkiLqHVPUplcDQocLhAuvW\ncTkb91UyMjLw+eef02yZEB4YE2bJjAF//FGh/9c0QyZPZ8UymQwJCQkUxqqyeDGQkwOsWqWyMAYA\nCwsLmi0TwouWlvCs8qtXgZUrxStLM+Sq7flZ8Zo1azB48GB1D6n6SEgA3NyAixeF5xurCc2WCeHk\nxg3A0RH4+2+gVaty/RI6GKSaunTpEkaMGAE7OzusWrWq2jyL+MmTJ7h//z7kcjkUCsXTl7a2NnR1\ndaGrqws9PT3Url0bBgYGfK6fy+VA166Avz8wbpz49SuIMYbg4GDMmjULM2bMwKxZs2jfACFiWL0a\n2L4dOHNGeJb5W1AgV0O7d++Gn58fVq1ahZEjR6p7OKJ5+PAhEhMTkZGRgaysLMhkMshksqffz8rK\nQlFRERo1aoQaNWo8DWAdHR0wxqBQKJ4GdUFBAR4/fgwTExOYmZnB1NT06bdl32/VqhWsrKwqHl4L\nFgjn3x46pNKl6rfJyMjA0KFD0axZM2zatAm1atVS95AI0WxKJeDqCnh5AV9++dZPp0CuRpRKJQIC\nArBlyxYcOHAAXbp0UfeQ3tnDhw8RFxeHmJgYxMbGIjY2FhkZGWjfvj0sLS1fCs+ybxs2bFjuAH3y\n5Amys7NfCPSykM/MzERycjKKiopga2sLe3t72NnZwc7O7qWQDvD1BW7eFH5QWAjExwP29kCrVk/P\npJaKoqIifP7550hKSsJff/1Fj80kpLIqsHRNgVxNFBYW4qOPPsLdu3exd+9eGBsbq3tIFZKeno7Q\n0FCcPXv2afh27NjxhSBs06YNdHVV+9TQ7Ozsp28Iyt4clIW0o6Mj3N3dcWzuXMw/ffqlXxvQqxcC\nOB0iUBmMMSxfvhzLli3Dnj170K1bN3UPiRDNVs6lawrkaiAtLQ2DBg2Co6MjAgMDNeKxiIwxXLx4\nESEhIQgJCUF6ejo8PDzg6uqqtvAtr7KQPn/+PA4ePIjGly/jmELx0udJNZDLHD58GJ988gkWLVqE\ncRK41k2IxlIqhU2cAwe+cemaArmKi4yMxMiRIzF37lz4+/tLerNOUVERTpw4gZCQEBw8eBB16tSB\nl5cXvL294ezsLNkAfpvZTk5YdOHCSx//xskJP6voWarv6urVqxg0aBDc3d2xdOlSjf07IETtyrF0\nTfchV2FBQUEYMWIEtm3bhilTpkg2jJOTkzFlyhSYmZlh0aJFaNGiBY4fP46rV69iyZIl6NGjh0YH\nQc3XrEjExMSgd+/e2L17t2TvA27dujX++ecfXL16Fe7u7njw4IG6h0SIZrK2BgIChEerlpZW+JdT\nIGsopVKJKVOmYPXq1Th37hzc3NzUPaSXyOVy7N69G71790afPn1gYGCAhIQEnD59GjNnzkSrct63\npxFu337lh8uexLRmzRo0bdoU8+bNQ2ZmpooH93YNGzbEwYMH0alTJ3Tt2hXXr19X95AI0Ux+foC+\nvvAgmQqiJWsNVFpaivHjxyMtLQ2hoaGoV6+euof0goyMDGzYsAEbNmxAq1at4OfnBx8fH+jp6al7\naHzExyPA0RGwtRX+Iz7PyurpLuvLly8jKCgI27dvh6urK/z8/NC7d2/JrWqsXbsWP/74IyIiIqrW\nmyZCVOXGDeEcgrNnX1q6fmP2MTVRY2uNVlJSwkaMGMH69OnDHj16pO7hvCAhIYENHTqUNWzYkPn7\n+7OkpCR1D4m/khLGOndmbNOmcv+Shw8fsqCgINa+fXvWqlUr9ttvvzG5XM5xkBW3efNmZmZmxhIT\nE9U9FEI00+rVjDk7M6ZQvPDhN2UfBbIGKSoqYj4+PszDw4M9efJE3cN5Ki0tjX300UfMyMiILV++\nnD18+FDdQ1Kd+fMZc3dnTKms8C9VKpUsMjKS9ezZk7Vp04bt27ePKd+hDi87duxgxsbGLDY2Vt1D\nIUTzlJYy5uLC2NKlL3z4TdlHS9YaoqSkBEOGDEGNGjWwc+dOSdzWdOfOHfz444/4448/4O/vjy+/\n/FJyy+dc3bgBODgIh4BYWLxzGcYYwsPD8c0336BmzZpYtGgRXFxcxBtnJezfvx8TJ07EoUOHYGdn\np+7hEKJZyr5GJCQA5uYAaJe1xpPL5Rg5ciT09PSwa9cutYdxQUEBAgIC0KZNGzDGcOXKFcyfP796\nhTEAfP89MHVqpcIYEP6Duru74+LFi5g6dSrGjRsHd3d3xMXFiTTQdzd48GCsW7cOnp6eSEhIUPdw\nCNEs1tbAp58Kz0IvBwpkiSstLcXHH3+M4uJi7Ny5U60bo5RKJdatW4eWLVvi+vXriImJwcqVK2Fk\nZKS2MalNfDwQEQHMmCFaSW1tbYwaNQpXr16Fp6cnPDw8MHr0aGRnZ4vW4134+PhgxYoVGDBgAK5c\nuaLWsRCicWbNAvbvB1JS3vqpFMgSplQqMX78+KdHYer/dwevCqWlpaFPnz4IDg7GkSNH8Pvvv6NZ\ns2ZqG4/azZkjvOrWFb10jRo14O/vj9TUVFhZWaFTp07Ytm2bWi/xjBgxAosWLULfvn3x77//qm0c\nhGgcAwPh5K5vv33753K7oP0WamytMaZNm8Z69OjBCgsL1TaG0tJStmbNGmZoaMiWLFnCFP/ZMVgt\nnTrFmJUVY0VFKmkXHR3N2rdvzwYNGsRkMplKer7OunXrWNOmTdU+DkI0yqNHjJmZMRYV9cbsoxmy\nRP322284fPgwQkJCULt2bbWMoWxW/Pvvv+PMmTOYOXMmdMrxvM8qjTFhCeqHH16+55gTe3t7xMTE\noH379mqfLX/++ecYO3YshgwZguLiYrWMgRCNU6uWsOfkm2/e+Gm0y1qCzp49iyFDhuDMmTNqOZhB\nqVRi7dq1mDdvHmbNmoXp06dTEJf56y/gu++AuLhyPYxcbDExMRg7dixatGiBoKAgmJiYqHwMSqUS\nw4cPR7169bBx40bJHWxCiCTJ5UC7dtBKTaVd1pri9u3bGD58OH7//Xe1hPGdO3fQr18/mhW/Smmp\ncN3455/VEsbAs9lyu3bt0KlTJ4SEhKh8DNra2ggODkZsbCxWrlyp8v6EaCQ9PWDhwjd+Cs2QJeTR\no0fo3r07PvroI8wQcfdueV26dAk+Pj746KOPEBAQQEH8X5s3C69TpwAJzAr/+ecfDBs2DJMmTcKc\nOXNUPlO9efMmnJ2dsWXLFvTr10+lvQnRSEoltHR06PGLUscYw4gRI/Dee+8hODhY5V9cd+/eDT8/\nP6xZswbDhw9XaW+NUFQE2NgAO3cC3bqpezRPZWVlwcfHB9bW1ti0aRNq1aql0v6nT5/GBx98gLNn\nz6Jly5Yq7U2IJqKDQTTAjz/+iNu3b2PdunUqDWOlUol58+Zh5syZOHr0KIXx6wQGAl26SCqMAcDM\nzAynTp2Cnp4eunfvjvT0dJX279mzJ3744Qd4e3sjPz9fpb0JqWpohiwBBw4cwJQpUxAVFQVTU1OV\n9S0sLMTHH3+MO3fuYO/evTA2NlZZb6kL8PUFbt4UfqBQABcuAJ07A+3aPX16k5QwxrBs2TIsX74c\ne/bsQTcVv3Hw9/dHWloaQkJC6FIHIW/wpuyjQFazpKQkuLq6IiwsDA4ODirrm5aWhkGDBsHR0RFr\n1qxR66Ej74oxBplMhps3b0ImkyErK+uFb+/evQu5XA65XA6FQgEdHR3o6upCV1cX9erVg6mpKUxN\nTWFmZvb0W3Nzc7Ro0QI/9u2LgFOnXuoZ0KsXAiIjVf+bLadDhw7B19cXixYtwrhx41TWVy6Xo3//\n/nBwcMDixYtV1pcQTfOm7NNV8VjIc4qLizFy5Ej88ssvKg3jpKQk9O/fH7NmzcKUKVM04rYVxhiy\nsrIQGxv79BUTEwOFQoEWLVq8EKw9e/aEqakpGjduDH19fejq6kLn/zdSKBQKyOVy5OfnQyaTPQ3w\nxMREZGVlIT09HVlZWXDV1syrOR4eHjh9+jS8vLyQmZmJ7777TiV99fT0sHv3bnTp0gX9+vWDm5ub\nSvoSUpXQDFmNvv32WyQlJWH//v0qC8W4uDh4eHhg+fLl+PDDD1XS810VFhbi6NGjCAkJwdGjR1FS\nUgI7OzvY29vDzs4OdnZ2sLS0FP3PrqCgAN/16IFf4+Nf+rnx1tb45sgRtGjRQtSeYsvOzoabmxt8\nfHywcOFClf37Cg8Px8SJE5GYmIi6HI4VJUTT0ZK1BMXExMDT0xPx8fEqO9whKioKXl5eCAoKwpAh\nQ1TSs6IyMjIQGhqK0NBQnD17Fk5OTvD29oaHhweaNWumsmAJcHF55ZL1hyYmiATQsGFDeHt7w9vb\nG127dpXkddO7d++iX79+cHV1xdKlS1X2Zzd+/HjUqFEDQUFBKulHiCahXdYSU1xcDF9fXyxfvlyl\nYTxw4EBs3LhRcmFcUlKCP//8Ey4uLujUqRPOnTsHX19fpKen4+jRo/D394e1tbUkltZbtWqFzMxM\nBAcHQ1dXF5MmTUKTJk3w/fffIyMjQ93De0Hjxo1x4sQJnDlzBtOnT1fZG+Dly5cjLCwMx48fV0k/\nQqoKmiGrgaqXqi9duoT+/ftj48aNGDhwIPd+5ZWRkYH169djw4YNaN26Nfz8/ODj46PWR0yWCfD1\nBUJDAUND4Pmd71ZWL+2yvnz5MoKCgrB9+3a4uLjAz88Pbm5ukngDAQB5eXlwdXVFv3798PPPP6tk\nXLR0TcirvTH7xHiQxbtQY2u1io6OZkZGRip7Wk5SUhIzMTFhe/bsUUm/8jh37hwbPHgwa9iwIfP3\n92eXL18K1dmgAAAgAElEQVRW95BedvMmY40aCU9pKaeHDx+yoKAg1qFDB2ZjY8NWrlzJnjx5wnGQ\n5Xf37l3Wvn17FhAQoLKe48aNYxMmTFBZP0I0wZuyjwJZhYqKili7du3YH3/8oZJ+t27dYubm5mzb\ntm0q6fc2SUlJzNvbmzVp0oQFBgaygoICdQ/p9ebMYeyLL97plyqVSnbmzBnm5eXFLC0t2ebNmyXx\n2Mrs7GzWqlUrtmLFCpX0y8vLY02aNGHHjh1TST9CNAEFskTMnTuXDRo0iCmVSu69CgsLWefOndmy\nZcu493qbmzdvMl9fX9a4cWO2bNkyycwaX6uoiDFjY8auXq10qTNnzrD333+ftWvXjv31118q+bt/\nk7S0NGZiYsKOHDmikn6HDh1iTZs2ZQ8fPlRJP0KkjgJZAlS5VK1UKtkHH3zAPvnkE7UGQG5uLps+\nfTozMDBgc+fOZXl5eWobS4Vs385Ynz6ilVMqlSw0NJS1b9+edevWjZ0/f1602u/i1KlTzMjIiF27\ndk0l/WjpmpBnKJDVTKlUMkdHRxYcHKySfj/88APr2rWrWmeihw4dYhYWFuzTTz9lWVlZahvHO+ne\nnbF9+0Qvq1Ao2JYtW5iJiQmbMWMGe1SB69NiW7t2LWvdurVK3iTl5eUxU1NTFh0dzb0XIVJHgaxm\ne/bsYZ07d2alpaXce+3bt49ZWFioLQTz8vLYuHHjWNOmTTXz2mF8PGMWFozJ5dxa3L17l40cOZK1\nbNmSnT17lluft/Hz82Oenp4qub4dGBjI+oi46kCIpqJAViO5XM5atWrFDh8+zL1XQkICMzQ0ZFFR\nUdx7vcrhw4dZkyZN2IQJEzT3muHEiYwtWKCSVnv37lXrbLmkpIS5uLiwWbNmqaRXixYtNPNNGiEi\nokBWow0bNjAXFxfu13Lv3r3LmjVrppYd1cXFxWzixImaOysuk5/PWMOGjKlwdaFstmxjY8MSExNV\n1vf5/s2aNVPJzv+dO3cyOzs7tW9sI0SdKJDV5PHjx8zc3Jz9888/XPsolUrm4eHBvv76a659XiUn\nJ4d1796dDRo0iOXn56u8v6hWr2Zs+HC1tN66dStr3LgxO3DggMp7JyYmMkNDQ3blyhWufUpLS5mt\nrS37888/ufYhRMookNXkl19+YYMHD+beZ9OmTaxLly6spKSEe6/nxcXFMUtLS/btt9+q5Po4dx07\nMnbihNraX7hwgVlYWLCFCxeqfBa5evVq5uTkxP168pEjR1jLli1V/m+VEKmgQFaD3Nxc1rhxY5ac\nnMy1T3p6OmvcuDGLj4/n2ue//vzzT2ZoaFh1ZjvXrzNmZMSYmg/wyMzMZI6OjmzEiBEqva5cWlrK\nXFxc2JIlS7j2USqVzNXVla1bt45rH0Kk6k3ZR2dZc/LNN9/g7t27+O2337j1YIzB09MTzs7OKnvu\nLWMMCxYswObNm7F//3506dJFJX25W7ECSEgANm5U90hQVFSEzz//HElJSQgLC4Pp82dpc5SWlgYH\nBwecPXsWrVu35tYnOjoaPj4+SE1NRa1atbj1IUSK6GlPKpaVlYX169dj3rx5XPsEBwcjOzsbs2fP\n5tqnDGMMX375JQ4cOIALFy5UnTAGgJAQwNtb3aMAANSsWRNbtmzBkCFD0KtXL6Snp6ukb7NmzTB/\n/nyMHTsWpaWl3Po4ODjA2dkZq1at4taDEE1EM2QO/Pz8ULt2bSxZsoRbj4yMDNja2iIiIgIdO3bk\n1qeMUqnElClTEB0djSNHjqBhw4bce6pMXh5gaQnIZEDt2uoezQuWLVuGNWvW4Pjx42jWrBn3fkql\nEm5ubvD09MTMmTO59UlJSUH37t2RmpqKBg0acOtDiNS8KfsokEWWm5sLa2trpKSkwMjIiEsPVS9V\nM8bg7++PuLg4hIeHo169etx78hbg6wvcvCn84M4dICcH6NDhlY9XVLfVq1dj2bJlOHXqFCwtLbn3\nS0tLg6OjI86cOcN16XrUqFHo2rUrpk2bxq0HIVJDj19UoeXLl7PRo0dz7bF161aV7apWKpVsxowZ\nzMHBQXPOoi6Heb16MQa89JrXq5e6h/ZKy5cvZy1atGCZmZkq6bdmzRrWrVs3rru9z5w5w2xsbOi+\nZFKtvCn7dFX5zqCqUyqVCAoKQjDHGVZRURHmzp2LXbt2QU9Pj1ufMj/99BOOHz+OEydOoH79+tz7\nVUZpaSmysrIgk8mefpuTk4OSkhIoFAro6OhAV1cX9evXR05OjrqHWyHTp09HcXEx+vbti3PnznH/\nu5g4cSI2bNiA/fv3Y8iQIVx6vP/++9DX18eJEyfg5ubGpQchmoQCWUTHjx9HrVq14OzszK3HmjVr\nYGtry7VHmf3792Pt2rWIioqCgYEB934VdePGDZw7dw6xsbGIjY1FXFwc6tatC3Nzc5iamsLU1BTG\nxsaoXbs2dHV1oVQqIZfLkZWVhfv377+yZnp6Os6ePYuuXbuq5A1PRcyePRsZGRkYNWoUQkJCoKOj\nw62XtrY2fv75Z3zxxRfw9vaGrq74Xyq0tLTg5+eHwMBACmRCQNeQRTV48GC4u7vj888/51I/Pz8f\nLVu2xMmTJ9GuXTsuPcokJibC1dUVhw4dgoODA9de5VVaWooLFy4gJCQEISEhyM3NRc+ePWFnZwc7\nOzvY2tqWe7NZgIsLAk6deunjo83MkGJqiuvXr6N///7w9vaGu7u7ZDaxyeVy9O/fHw4ODli8eDHX\nXowxuLq6YsyYMRg/fjyXHgUFBWjatCkSEhJgYWHBpQchUkLXkFXg9u3bzMDAgBUUFHDrMXfuXObr\n68utfhl1nov9Krdv32bffvstMzExYR07dmRz585lFy5cqNTpYG+7hpyVlcU2bNjAvLy8WN26dZmP\njw87duyYJE4ku3fvnsrOnz5//jyzsLBgjx8/5tZj8uTJ7Pvvv+dWnxApeVP20QxZJN999x3y8/Ox\ncuVKLvVlMhnat2+PuLg4rjtt5XI5+vXrh65du2LRokXc+rwNYwwnT57EqlWrcOrUKYwZMwYTJ05E\n27ZtRakf4OsLXL4sHAbi5ARoaQk/8Ypd1oWFhdi2bRsCAwNRVFSESZMmYfz48ahbt64oY3kXSUlJ\n6N27t0pWMIYMGQJnZ2d89dVXXOpfvnwZffv2xa1btyR3mYAQsdEMmbPi4mJmYmLC9ZhMPz8/Nn36\ndG71y0ybNk1lz8h9nfPnzzMXFxfWqlUrtm7dOn6rDj/9xJi/f7k/XalUsrNnz7IRI0YwY2NjtmLF\nClZUVMRnbOVw4MABZmFhwe7evcu1T3JyMjM0NGS5ubncevTq1avqHMNKyBu8KfsokEWwc+dO1rt3\nb271U1NTWaNGjbh/4T1x4gQzNzdnDx484Nrnda5cucIGDx7MzM3N2YYNG5hcLufb0MmJsaNH3+mX\nxsfHMw8PD2ZlZcW2bt2qtqXs6dOns5EjR3LvM27cOPbNN99wq79r1y7m4uLCrT4hUkGBzFnPnj3Z\n7t27udUfN24cCwgI4FafMcYKCgpYs2bNWGhoKNc+ryKXy9mPP/7IDA0N2eLFi7ler3wqO5ux+vUZ\nKy6uVJlTp04xBwcH1rt3b3bjxg2RBld+jx49Yi1btmR79+7l2uf27dusQYMG3N6sla0yXb58mUt9\nQqTiTdlHZ1lXkkwmQ2JiIrw5nYOcm5uLffv2YdKkSVzql5k9ezZ69OiBgQMHcu3zX0lJSXByckJk\nZCRiY2Px9ddf47333uPf+O+/gfffB2rUqFSZnj174vz583B3d4eDgwMCAwOhVCpFGuTb1apVC5s3\nb8bkyZNx7949bn2aNGkCd3d3bNmyhUv9GjVqYNSoUdi9ezeX+oRoAgrkSjp48CD69++PGpX8wv46\nwcHB8PT05HYMJwCcPHkSBw4cwK+//sqtx38xxrB69Wr07t0bEyZMwJEjR1RyLORTsbGAvb0opXR0\ndPDVV1/h7Nmz+P3339G/f//X3ufMw/vvv49Ro0ZhypQpXPtMnjyZ6xsOb29vhISEcKlNiEZQ3UT9\nRWpsLaqBAwey7du3c6ldWlrKWrZsyf7++28u9RlTz1J1cXEx+/TTT1n79u3Z9evXVdb3Bf37M/bX\nX6KXVSgU7Ouvv2bW1tYsMTFR9Pqvo4qla6VSyTp27MiOHTvGpb5cLmcGBgYsPT2dS31CpOBN2Ucz\n5Ep49OgRTp06hQEDBnCpr4qTv5YsWQInJyeVLVXfuXMHbm5uuHfvHs6dOwdra2uV9H0BY8IM2c5O\n9NI6OjpYvHgx5s+fD1dXV/z111+i93iVWrVqYePGjZg+fTqKioq49Hj+ZC0edHV14eHhgdDQUC71\nCZE6CuRKiIiIgIODA7dTnAIDA+Hn5wetsntkRZaTk4PVq1fjp59+4lL/v9LS0uDs7AwXFxfs3btX\nfffxpqcDOjqAmRm3FmPGjMHBgwcxefJkBAUFcevzvB49eqBLly7cAhMARo8ejcjISGRkZHCpT8vW\npDqjg0EqYfz48ejYsSOXx8elp6ejc+fOuHXrFurUqSN6fQCYOnUqtLS0sGLFCi71n5eamgo3NzfM\nnj0bfn5+3Pu90b59wMaNQFgY91Y3btyAm5sbpk6diunTp3Pvd/nyZbi6uuLatWvcHkAxZcoUNGzY\nEAsWLBC99sOHD2FhYYHMzEy1HrxCCC9vyj6aIb8jpVKJsLAweHl5cam/fv16jB49mlsY37hxA9u2\nbcPcuXO51H9eWloa3Nzc8N1336k/jAFuy9WvYm1tjdOnT2P16tVYvXo1937t2rWDh4cHli5dyq3H\npEmTsGHDBpSUlIheu169enB2dsaxY8dEr02I1FEgv6OoqCgYGhpyuQaqVCqxceNGTJw4UfTaZb7/\n/ntMnTqV6+5tALh37x769u2LWbNm4bPPPuPaq9xUGMiAcMvQ8ePHsWTJEmzbto17v4CAAAQGBiI7\nO5tL/bZt28LGxgaHDh3iUt/Ly4uWrUm1RIH8jkJCQrjdexwdHY0GDRqIdm7zfyUkJCAiIgIzZszg\nUr+MXC7HsGHD8MEHH2Dy5Mlce5Vb2YYukW55Ki8rKyuEhYXhiy++QHR0NNdeTZs2xSeffIKFCxdy\n6zFs2DBuG9a8vLwQFhaG0tJSLvUJkSoK5HfEM5B51gaAX3/9FV988QX3a3TTpk1D3bp1uQZDhalg\nQ9frtG/fHr/99huGDBmCrKwsrr1mzZqFbdu2ITc3l0t9nqHZtGlTmJub4/z586LXJkTKKJDfQWZm\nJrKzs+Ho6MilfmhoKLdAfvDgAfbv38/t+bZl1q9fj8jISGzbtg06Ojpce1VITIywXM1p5/rbDBo0\nCBMnTsTgwYNRXFzMrY+xsTE8PDy4naxlZWUFU1NTXLhwgUt9Dw8PHDlyhEttQqSKAvkdxMTEwMHB\nAdra4v/xpaWlITs7G127dhW9NiCc/DVw4EA0btyYS31A2FE9d+5cHDhwAPXq1ePW552o+Prxq8yZ\nMwfGxsbcVw7K7hnmdbIWz2u9Xbt2RUxMDJfahEgVBfI7iI2NhR2nL+qhoaEYOHAgl1mlUqlEUFAQ\n153OpaWlGDt2LL777jvY2Nhw6/POJBDIWlpaWLduHdavX4/Y2Fhufbp164b33nsPJ06c4FKf5z3D\ndnZ2iI2N1fhbIwmpCArkd8AzkHleP46IiECdOnXg5OTEpT4ArFq1Ctra2vD39+fW451xPKGrokxN\nTbF8+XL4+vpyW7rmfbKWvb09cnNzkZqaKnptc3NzaGlpcTuAhBApokCuIMYYt0DOz89HVFQU+vbt\nK3ptgP/JX7du3cLChQuxadMmLsv5lZaRAWhrA+bm6h4JAGDUqFGwtrbG4sWLufUoO1krPT1d9Nra\n2trw8vLictSllpYW7O3tadmaVCsS/KopbZmZmSgtLUWTJk1Er33kyBH06NEDtWvXFr12fn4+jh8/\njpEjR4peu8z3338PPz8/tGjRgluPSklLA1q2VNuGrv/S0tLCr7/+ihUrVuDOnTtcetSpUweDBw/G\n3r17udTnFcjAs2VrQqoLCuQKio2Nhb29PZdZ5rlz59CrVy/R6wLPwp7XrU5JSUkIDw/HzJkzudQX\nRUYGYGGh7lG8oFmzZhgzZgzXDV6DBg3idq23Z8+eiI6O5nL7EwUyqW4okCuI5/Vj3pvFeN7bPGfO\nHMyePVt6u6qfl5kpuUAGgLlz52Lbtm24ceMGl/p9+vRBTEwMl3uS69evDzMzM1y9elX02rSxi1Q3\nFMgVxCs0S0tLcenSJdja2opeW6FQ4NChQ9wesRgTE4NLly5h0qRJXOqLJiNDMtePn2dkZAR/f39u\nT92qVasWXFxcEB4ezqW+nZ0dl2u95ubm0NbWpo1dpNqgQK4Anhu6UlJSYGxszOVRjn///TesrKxg\nwWl2uGbNGvj7+6NmzZpc6otGgkvWZSZPnoy9e/dyO1lLFbcoiU1LS4tb2BMiRRTIFcBzQ5em3kr1\n4MEDHDhwAOPGjeNSX1SZmZKcIQPCLNnT0xPBwcFc6g8cOBDh4eFcntDE81ovXUcm1QkFcgVcu3YN\nbdu25bKhi2cgR0REYMCAAVxqBwcHw8vLC4aGhlzqi0rCM2RAOFkrKCiIy8laJiYmsLa25jLbtLW1\nRXx8PJeNXe3atUNKSorodQmRIgrkCpDJZDDnNMPiFchPnjxBamoqOnfuLHptANiwYQPXx0SKRqEA\n7twBTE3VPZLXcnZ2xnvvvYfTp09zqe/g4MBltslzY5eZmRlkMpnodQmRIgrkCsjKyoIphy/oPDd0\nxcfHo3Xr1tDX1xe9dkpKCgoKCuDs7Cx6bdHl5ACNGgF6euoeyWtpaWnhgw8+4PZYQ01cWjY1NeX+\nZCxCpIICuQJkMhnMODy2Ly0tDY0aNeKyoUsV16Z5nfwlKokvV5fx9vbGX3/9xeVWH56B3LlzZyQk\nJIhe19TUFDKZjG59ItUCBXIF8JohZ2ZmctkoBgi3JPEMZC8vLy61RSfRe5D/q0OHDigtLUVycrLo\ntdu3b4/r16/j0aNHote2sLBAZmam6HVr166NGjVqIC8vT/TahEgNBXIF8Johy2QyLkEP8Jsh5+bm\nIiEhAb179xa9Nhd37gBGRuoexVtpaWnB29uby3GUNWrUQLt27RAfHy96bZ7XestmyYRUdRTIFcBr\nhpyVlcUl6IuLi5GamooOHTqIXjs6Ohq2trbSv/e4TEEBwOnYULF1794dFy5c4FLb1tYWFy9eFL0u\nz2u9ZmZmdB2ZVAsUyBWgaTPkrKwsGBsbcwlNntemuSgsBOrUUfcoyoXntd5mzZrh9u3botelGTIh\nlUeBXE4FBQUAwOXhDDwDmddSOM9r01xo0Ay5efPmePjwIe7evSt6bV6zzbp160KpVD79fyImCmRS\nXVAglxOvZWWetXlfm7a3t+dSmwsNCuSyIyN53UbEI9y0tLS4zZJpyZpUFxTI5cQz3HjOkHkEfUlJ\nCbKystC8eXPRa3OjQUvWANC6dWukpqaKXpfnbJNXbZohk+qCArmc8vLy0KBBAy61NW2GnJ2dDSMj\nI2hra9A/Hw2aIQP8QojnbJNX7YYNG9JtT6Ra4PYVNTw8HK1bt0bLli2xePFiXm1URqFQQI/DKU9K\npRIPHz7kEva8gp7n8j03GjZD5hluRUVFePLkiei1DQwM8ODBA9Hr6urqQqFQiF6XEKnhEsilpaXw\n9/dHeHg4kpOTsWPHDly5coVHK5VRKBTQ1dXlUldPT4/LaVf3799Ho0aNRK/Lc/meG7kcqFFD3aMo\nN57Xehs1aoR79+6JXltPT4/LAyYokEl1wSWQo6Ki0KJFC1hZWUFPTw8jR47kdj6vqsjlci4zZLlc\nziXoy2rX4BBCd+/eRePGjUWvy0tYWBj6X74MF39/9O/fH2FhYeoe0lsZGRlx2WUNCMHJI+B0dXUh\nl8tFr8trvIRIDZck+O9RkBYWFtwOOlAV3jNkHnjVlsvlXB5WwUNYWBimTZuG6/n5wKVLAIDr168D\nADw9PdU5tDeqUaMGl3AD+AayJtUlRGq4zJA14mEDFaRQKKCjo6MxdXnW5jlmsa1cufJpAJe5fv06\nVq1apaYRlY+Ojg63EOJVW9PqEiI1XGbI5ubmSE9Pf/rj9PR0WLziYP+AgICn33dxcYGLiwuP4YhC\nR0eHy4PjdXR0uFx3K6vNa8w86vJQXFz8yo8XFRWpeCQVo1Qqub3p4VVb0+oSogqRkZGIjIws1+dy\nCWR7e3ukpqbi5s2bMDMzw65du7Bjx46XPu/5QJY6PT09LkuIPJfjeF3T41WXh9ctrUv9DG7eewt4\nXX7hUZfnnwUhvP13sjl//vzXfi6XJWtdXV2sXr0a/fv3R9u2bTFixAi0adOGRyuV4RWcPDes8HoT\nYWBggPv374tel4epU6e+dIBJ8+bNMWXKFDWNqHzu378PAwMDLrU1bT8Er/ESIjXc/pW7u7vD3d2d\nV3mV08QNK/Xr1+dyoIImnZxUtnFr1UcfocjSEjVNTDBlyhRJb+gC+N5alpuby+W+d00LekKkht52\nlhOv4NTW1oa+vj4KCgpEf3AFr8MlNO1sYU9PT3g6OQH+/oCHh7qHUy68Dl8pKCgAY4zLQ1Ly8/O5\n1KUZMqkuNOjsQ/WqU6cOlyfZ8DyUn+fZwtnZ2WCMiV6bmzp1hOMzNQTv88153AnB801E7dq1Ra9L\niNRQIJeTJh7Kz2smW6tWLdSrVw+ZmZmi1+ambl2NCuTr16/DyspK9Lqa+JAUjTwZjpB3QIFcTrwP\n5dekGTIA2Nra4uLFi1xqc1G3rnCetYaIjY3l8rxpnuFGZ6cTUjkUyOXUsGFDPHnyhMuh/Kamptyu\n9fIKZDs7O8TExHCpzYUGLVlnZWWhpKQElpaWXGrzCLfHjx+juLiYy2YxmiGT6oICuZy0tLRgYmKi\nUTNZMzMzZGZmcjl4xN7eHrGxsaLX5UaDlqzLZsc8rvOmp6fD3Nxc9LoymQwmJiZcxkyBTKoLCuQK\n4DXj5LUcXqdOHZiZmSElJUX02g4ODoiKiuJ2ypjoNGjJ+p9//oGDgwOX2nFxcejcubPodWUyGbdl\nZVqyJtUFBXIF8Fpa5nmt187OjstM1tzcHKampoiKihK9NhcGBgCnpyeJ7eDBg1zu4VcqlYiLi4Ot\nra3otbOysjRusxghUkOBXAE8Z8i8dizzvNbr5eWFkJAQLrVFZ2EBaMCu8Js3b0Imk8HJyUn02teu\nXYOhoSGXE8B4BXJxcTEePnwIQ0ND0WsTIjUUyBXAa4bcvHlzpKenc9kwxvNar7e3NwWyyEJDQ+Hp\n6cnlYQqxsbGwt7cXvS4AJCYmom3btqLXzc7OhrGxMbS16UsVqfroX3kF8Joh6+vro3Xr1oiPjxe9\ntq2tLeLj47lc63VwcMCDBw+4XKMWnZkZIJMBEn9K1b59++Dt7c2lNq9bqXjW5nltmhCpoUCuAF4z\nZIDfTLZBgwYwNjbG1atXRa+tra0NX19frF+/XvTaoqtRA2jYELhzR90jea2UlBQkJyfDg9PxnrxC\n88mTJ7h27Ro6duwoem2e16YJkRoK5AqwtrbGtWvXuNTmea3XxcUFERERXGpPmDABW7ZswePHj7nU\nF5W5OZCRoe5RvNbatWsxfvz41z4ysjIePnyIuLg4dO3aVfTaCQkJaNWqFZdHWl67dg3W1tai1yVE\niiiQK6B58+bIy8vDvXv3RK/Nazc0wHfzlZWVFZycnLBr1y4u9UVlYSHZQH706BG2bt2KCRMmcKl/\n9OhRvP/++1we/qCJS+GESBEFcgVoa2ujS5cuXIKzQ4cO+Pfff7ls7OrTpw+io6ORm5srem0A8PPz\nw8qVK6X/sAlzc8lu7Prjjz/QrVs3NG3alEv9kJAQujZNiMRRIFcQr2u9PDd21a5dG7169UJ4eLjo\ntQFgwIAB0NHRwZ49e7jUF41EZ8hPnjzBwoUL8e2333Kpr1AocOjQIXh5eXGpzys0Hzx4gHv37sHG\nxkb02oRIEQVyBfFcWtbUW5S0tbXx888/Y+7cuZDL5Vx6iEKigbxq1So4Ojpyub4LAOfOnUPTpk1h\nYWEhem2eG7piY2PRpUsXuuWJVBv0L72CeAayo6Mjzp49y6X2wIEDceTIEW6B2bdvX1haWmLTpk1c\n6otCgkvWubm5WLJkCRYuXMitB8/l6ujoaLRt25bLhi5aribVDQVyBbVo0QK5ublcNna5u7tzC01T\nU1O0bdsWBw8eFL12mZ9//hnz58/HgwcPuPWoFEtL4MYNdY/iBfPnz4ePjw/atGnDpb5CocDOnTsx\nZMgQLvXLDjLhgQKZVDcUyBXEc2OXubk5mjdvzm2WPHHiRAQFBXGpDQgHhQwbNgxffPEFtx6V0rw5\ncP8+IJE3DGfOnMGff/6JRYsWcesREhICKysrdOjQgVt9TdwsRogUUSC/A0291jts2DBcunSJ273U\ngDBL/vvvvxEaGsqtxzvT0QG6dAEk8NjIx48fY9y4cQgMDESjRo249QkMDISfnx+X2ikpKSgsLOTy\nsIr79+/Thi5S7VAgvwPe9wyHhoZyuYWoZs2aGDduHNauXSt67TK1a9fGpk2bMGnSJGkuXdvZSSKQ\n586dC0dHR/j4+HDrcfXqVSQmJmLo0KFc6oeGhsLLy4vLM5AvXrxIG7pItUP/2t+BnZ0doqOjuYRm\np06dUFJSgitXroheGxCWrbdu3cr1ZK1evXphxIgRGDVqlPSelyyBQN6zZw/27t2LlStXcu3D8+Qv\ngO9ydUxMDC1Xk2qHAvkdtGzZEowxLudDa2lpcV22trKygrOzM3bu3MmlfpnFixejtLQUs2bN4tqn\nwtQcyJcuXcKkSZNw4MABrkvVjx49wu+//87t5K/79+8jPj4erq6uXOqHh4dzq02IVFEgvwMtLa2n\nS8s88H6s4dSpU/HLL79AoVBw66Grq4tdu3bhwIED2Lp1K7c+FWZjA9y7p5aNXXfu3IGPjw9Wr17N\n5cpkDPYAACAASURBVLrr89auXQsXFxduJ38dOnQIbm5uXG53un//PuLi4uDm5iZ6bUKkjAL5HfEM\nzV69eiE5OZnbk6X69OkDU1NTBAcHc6lfxsDAACEhIZg5cyaOHz/OtVe56egAnTsDFy+qtG1BQQEG\nDRqEMWPGYMSIEVx75efnY/Hixfjhhx+49di/fz+3k78OHz4MV1dXvPfee1zqEyJVFMjvqHfv3khM\nTMTdu3dFr62vr48RI0Zg48aNotcGhBn+okWLEBAQwOXs7Oe1bdsWe/bswciRI3HmzBmuvcpNxcvW\njx8/hre3N9q3b48FCxZw77dkyRIMHDgQbdu25VI/KysLJ0+e5HZvM89r04RIGQXyO9LX10ffvn0R\nFhbGpf6kSZOwbt06bsvKXbt2haOjI1avXs2l/vN69uyJHTt2YOjQoYiMjOTe763s7ABOj7r8r4KC\nAnh4eMDS0hJr167lvmtYJpMhKCgIAQEB3Hr89ttvGDlyJOrXry967eLiYhw9epTbYSOESBkFciXw\nXLbu3LkzmjZtyvV+3h9//BFLlixBXl4etx5l+vTpgz///BPDhw/Hvn37uPd7IxXNkLOzs9G3b1/Y\n2Nhg8+bN0NHR4d5z4cKF8PX1haWlJZf6crkc69evx6RJk7jUP3XqFNq2bQtjY2Mu9QmRMgrkSvDw\n8MDx48dRVFTEpb6fnx8CAwO51AaANm3awNvbm+tJUc9zcXHB4cOH8cUXX2DBggVQKpUq6fsSFWzs\nio2NhaOjI9zd3VUyMwaA1NRU/Pnnn5gzZw63HqGhoWjWrBmXh0mU1aflalJtMTVRY2tRde/enR06\ndIhL7aKiImZkZMSuXr3KpT5jjGVmZjJDQ0MWFxfHrcd/yWQy5uTkxIYNG8YKCwtV1vcFPXowduwY\nl9I7duxghoaGbM+ePVzqv0ppaSnr1asXW7p0Kdc+bm5ubPv27VxqK5VKZmlpyZKSkrjUJ0QK3pR9\nNEOuJJ7L1vr6+txP1jIzM8Mvv/wCX19flT060cTEBJGRkahTpw4cHBwQFRWlkr4v4LBsXVBQgIkT\nJ+Kbb75BREQEtxOyXiUoKAglJSVczxG/evUqkpKSuG3mSkhIgK6uLrfNaIRIHQVyJZUFMuNwahcA\nTJgwAVu3bsWjR4+41AcAX19fmJub46effuLW47/09fWxadMmzJs3D97e3pg9eza3pf9XEjmQIyIi\n0KFDBygUCly6dAmdOnUSrfbb3LhxA/PmzeN+nVpVJ3/xOIqTEI2guon6i9TYWnTt2rVjJ0+e5FZ/\n4MCBbN26ddzqM8ZYRkaGypeuy+Tk5LChQ4eyNm3asIiICNU0/fdfxkxMGCstrVSZ7Oxs9tlnn7Em\nTZqww4cPizS48itbql6yZAnXPvn5+czAwIDdvHmTS32lUslat27NTp8+zaU+IVLxpuyjGbIIeD/W\ncM6cOVi4cCHXGaS5ubnKl67LGBkZYffu3fjhhx8wceJE9OvXj9vDO55q3hxo2PCdb396+PAhvvvu\nO7Rt2xa1a9dGYmIiBgwYIPIg365sqXr69Olc+yxduhQDBw7kdvJXZGQkdHR00L17dy71CdEIKnxj\n8AI1thZdXl4ea9CgAcvKyuLWY9CgQWzZsmXc6jMmzFI8PT3ZjBkzuPZ5k5KSEhYYGMhMTU3Z8OHD\nWXR0NL9mX3/N2LffVuiX3L9/ny1evJgZGRmxjz/+mKWlpfEZWzlcunSJGRoact30x5iwCmBgYMD1\n9zps2DC2Zs0abvUJkYo3ZR8FskgmTJjAFixYwK1+UlISMzIyYnl5edx6MCYETosWLVhwcDDXPm9T\nWFjIfvnlF2ZlZcUcHBzY5s2b2ePHj8VtcvYsYx07lutTo6KimK+vL2vQoAEbM2YMS0hIEHcsFZST\nk8OaNm3Kdu7cyb2Xv78/mzZtGrf6mZmZrGHDhiw/P59bD0Kk4k3Zp/X/n6ByWlpa3DZCqUN8fDw8\nPT1x8+ZN6OrqcukxduxYmJubY+HChVzql0lOTkavXr0QGhoKJycnrr3eprS0FOHh4QgMDMSFCxfg\n5eUFb29v9O3bF3Xq1KlsccDUFIiKAqysXvgpxhji4uIQEhKCAwcOID8/H5MmTcLYsWPRuHHjyvWt\npJKSEvTp0wc9e/bk/m/hxo0bcHBwwNWrV7n9vufPn4+cnByu99wTIhVvyj4KZBF1794dX375JQYP\nHsyl/q1bt2Bra4vLly/DxMSES48yoaGhmDhxIqKiomBubs61V3ndvn0bISEhCAkJwT///IPu3buj\nV69esLOzg62tLQwMDCpedOxYwNYW8okTkZycjNjYWFy4cAGHDh1CzZo1MWjQIHh5eaF79+4qOWnr\nbRhjmDhxInJycrBv3z7uB46MGTMGLVu2xLx587jUl8vlsLKyQnh4ODp06MClByFSQoGsItu3b8fm\nzZtx7Ngxbj1mzJiB4uJirFmzhluPMosWLcLevXtx+vRpyT15Jz8/H0ePHsW5c+cQGxuLS5cuwdDQ\nEO3atYOZmRlMTU1hamoKY2Nj6OvrQ1dXF0qlEgqFArm5ucjKyoJMJoNFdDS6JybCVaGApaUl7Ozs\nYG9vjwEDBqBVq1aSuwUnMDAQgYGBOH/+POrWrcu1V3x8PPr374/U1FRuvfbu3YsVK1bg9OnTXOoT\nIjUUyCpSXFwMS0tLnD59Gq1ateLS4969e2jdujUuXLiA5s2bc+lRhjGG0aNHo7i4GLt27eK2FC8G\npVKJa9euISUlBTKZ7Gng5uTkoKSkBAqFAjo6OtDR0UGDBg2ehralgQEG+/vjSUoK6lpYqPu38UYH\nDx7Ep59+inPnzsHa2pp7P09PT/Tv3x9Tp07l1sPNzQ2fffYZRo4cya0HIVJCgaxCc+bMwZMnT/C/\n//2PW48ff/wRUVFROHDgAPcZXHFxMXx8fGBgYICtW7dKYtlWdB4egK8vMHy4ukfyWkePHsWYMWMQ\nFhYGBwcH7v0OHz4Mf39/JCcnczsI5MqVK+jduzdu376NGjVqcOlBiNS8KfvoPmSRqeJkrZkzZ+LG\njRvYvn07tx5l9PX1sW/fPshkMnz66acoLS3l3lPlvL0BTsefiuHEiRMYM2YM9u/f/3/t3XlYlOXe\nB/AvApmluQKCmuClpiQYZ8hccQVUZMg6WYq5b2CLWufVeuMVT+HSm3Zcosig0nJPYARDzYXclwkE\n0RNuJOKIehQlFJhh7veP54XEBQGZ53kGvp/rmotFmN9PR/jOfc/93LcsYZyXl4cpU6Zg5cqVFgtj\nQNr5a9KkSQxjov/HEbIFDB8+HH379rXovsJ6vR5Dhw5FamoqnJ2dLVanVEFBAYYNG4Y2bdogJiZG\n1dPXVXbxItC1K3D5MmBvr3Q35Wzbtg2jR4/Gpk2b0LdvX1lqTpgwAfXr17foZjdXrlxB586dkZKS\nYrGjIonUiFPWMktPT8egQYNw+vRpPPPMMxar89FHHyE9PV2WqWsAuH37NoYPH47GjRvjhx9+qF0j\nG29v4LPPgH79lO6kzJYtWzBx4kTExcWhZ8+estTcunUrpk+fjrS0NIsuGpsxYwZKSkqwfPlyi9Ug\nUiNOWcvMw8MD/v7+WLx4sUXrhIWFyTZ1DQBPPfUU4uPjYTQa4efnh6tXr8pSVxZaLbBli9JdAJAW\n0/3rX//ClClTkJiYKFsY5+XlYerUqYiOjrZoGGdlZWH16tX46KOPLFaDyBpxhGwhWVlZ0Gg0OHny\nJJycnCxWR+6pa0DarCMsLAxr165FfHy8xQ6rl1VKCvDaa8Dp04CClzoVFRUhJCQEer0e8fHxcL1n\nwxJLkmOqGgDGjBkDNzc3zJs3z6J1iNSIU9YKeffddyGEwLJlyyxaJywsDGlpabJNXZdau3Yt3nnn\nHURFRVnsjFzZCAG0bQskJAAKPcG4fPkyXnnlFbi4uOC77757/J3IqkCuqWq5Xs4hUqsKs+/xd+as\nHgVLyyY3N1c0a9ZMnD171qJ1ioqKhIeHh4iOjrZonQc5evSoaNOmjQgPDxclj3mUoeLmzhVi+nRF\nSh87dkyxf8fc3FzRunVrsXPnTovXGjZsmPj8888tXodIrSrKPgayhYWHh4vg4GCL18nIyBAtWrQQ\nBw8etHitexkMBtGzZ0/h6+sr/vjjD9nr15iLF4Vo2lSIW7dkK2kymcSnn34qWrRoIX766SfZ6pYq\nKioSvXv3Fh9++KHFa+3du1e0bdtWFBYWWrwWkVpVlH1c1GVhs2bNwi+//ILjx49btI67uztiYmLw\n6quvIicnx6K17tWyZUskJyejX79+0Gg0WLlypXW+HNGqFTBgAPDjj7KU+/e//43evXtj69atOHz4\nsOzT/kIIvPXWW2jWrBk+/vhji9eaM2cO5s2bZ9Frm4msmlzPCu6lYGnZLV26VAwdOlSWWgsWLBDe\n3t41f1RhJaWnpwuNRiP8/Pysc7S8c6cQHh5CmM0WK1E6Km7evLlYsWKFYlP9K1asEM8//7wsxx7q\ndDrRpUsXYTKZLF6LSM0qyj4GsgwKCwuFq6ur2L17t8Vrmc1mMXLkSDFq1ChhtmCoVMRoNIqIiAjR\nokUL8eWXXwqj0ahIH9ViNgvx3HNC7N1rkbtPS0sT3bt3F/369bP42oKK7Ny5Uzg6OoozZ85YvJbR\naBRdunQROp3O4rWI1I6BrAKxsbGiQ4cOoqCgwOK1CgoKhEajEQsXLrR4rYqkp6eLvn37is6dO4vN\nmzcr9gShypYuFWLkyBq9y/Pnz4s333xTODk5icjISEUXwJ09e1Y4OjrKsohLCCEiIiKEn5+f9Tz+\nRBbEQFaJkSNHilmzZslSKzs7W7i4uIiEhARZ6j2M2WwWW7duFV27dhUvvfSSLLMEj+3GDSGaNBHi\n8uXHvqvc3FzxzjvviGbNmom5c+eKWzIuGHuQW7duieeff14sX75clnrp6emiRYsW1vnyBZEFMJBV\n4urVq8LZ2Vns27dPlnoHDhwQDg4O4siRI7LUq0hJSYn48ccfhZubmxg8eLA4fPiw0i1VbPJkISIi\nqv3t165dE3PnzhXNmjUTb7/9tsjNza3B5qqnsLBQ+Pv7i8mTJ8syWi0uLhYajUasXLnS4rWIrAUD\nWUU2b94s29S1EELEx8cLR0dHkZKSIku9RykqKhLLly8Xrq6uwtvbW8TExMj2b1ElKSlCtGkjRBUW\nIZnNZnH48GExduxY0bhxYzFu3DhFXye+W3FxsQgKChKvvPKKKC4ulqUmp6qJ7sdAVpmRI0eKmTNn\nylZv48aNomXLluLEiROy1XwUk8kkEhMTRUBAgGjevLmYNWuWyMzMVLqt8nr0ECI+/pFfVlBQIKKj\no4VGoxFubm5i0aJF4urVqzI0WDlGo1G89tprIiAgQBQVFclSk1PVRA/GQFaZ0qnrvRZayfsgP/zw\ng3BxcREZGRmy1aysc+fOidmzZwsHBwfRt29fsWTJEllW/z7S6tVC+Ps/8I/y8/PF5s2bxbhx40Tz\n5s1FQECASExMVN1lPUajUQQHBwtfX19x584dWWpyqpro4SrKPu5lrZDY2FjMnj0bqampeOqpp2Sp\nuWrVKsyePRvbtm1T5YEQhYWF2LFjB3Q6HRISEtCsWTNotVpotVp069YNtra2cjcEPPsscOAA0L49\ncnJykJCQAJ1Oh71796J79+4IDAxEUFCQKs/0NRqNCA4ORl5eHuLi4mT7fzZ//nwkJycjKSlJ1r3V\niawBD5dQqVGjRqFly5ZYsmSJbDXXr1+Pd999F4mJidBoNLLVrSqz2Yxjx45Bp9NBp9MhJycH3t7e\n8Pb2hkajgUajwbPPPmuxX/h//vknUlJS8GR4OK4YDHhfCFy5cgVDhgyBVquFv78/GjdubJHaNaGo\nqAgjRoyA2WzGxo0b8eSTT8pS98SJE+jfvz/0er0qn6QQKY2BrFLXrl2Dp6cn1q9fjz59+shWNy4u\nDlOmTMGGDRvQr18/2eo+DoPBgGPHjkGv10Ov1+PYsWMwGo3QaDRo3749nJ2d4ezsDBcXl7L3HRwc\nUK/eg3eHzc/Ph8FgKLtdunQJBoMBFy9eRGpqKi5cuAAPDw/4deiAD+PicObnn9Gpe3fY2dnJ/Dev\nups3b2LEiBFo1KgR1qxZgyeeeEKWusXFxejZsyemTZuGSZMmyVKTyNowkFVs69atmDx5Mo4cOYJW\nrVrJVnfnzp0YNWoUwsPDERISIlvdmnTp0iXo9XpkZWXdF6wGgwH/+c9/YGtrCzs7O9jZ2cFsNsNk\nMsFkMqF+/frlwrv0fRcXF3h6esLd3R329vZSofHjpX2uP/lE2b9wJWRmZiIoKAj9+/fHsmXLZHsC\nIYTAtGnTkJubi9jYWE5VEz0EA1nlFi5ciJ9++gm//vorGjRoIFvdM2fOICgoCD4+Pli6dKlsIym5\nmM1mlJSUwGg0wmQylQvnevXqVT40LlwAvLyAEycAZ2fLNv0Ytm/fjtGjR+Pjjz/G1KlTZa0dGRmJ\nyMhIHDx40KLnKRNZOwayygkhEBwcDBsbG/zwww+yji5u3bqF4OBg5OfnY+PGjXBwcJCttlV57z3g\nzh0gMlLpTu4jhMDSpUuxcOFCbNiwAT4+PrLW37VrF0aNGoUDBw6gXbt2stYmsjYVZR+PX1QBGxsb\nREdH4/fff8enn34qa+1nnnkGcXFx6NmzJ7p164a0tDRZ61uNDz4ANmwAzpxRupNyioqKMGHCBHz7\n7bc4dOiQ7GF87tw5jBo1CmvWrGEYEz0mBrJKNGjQAHFxcVi2bBkSExNlrW1ra4v58+dj/vz5GDhw\nIDZv3ixrfavQogUwYwYQFqZ0J2UuX76M/v37Iz8/H/v374erq6us9fPz86HVahEWFoYBAwbIWpuo\nNmIgq0jr1q2xadMmjB8/HqdOnZK9/siRI5GUlISZM2ciJCQE+fn5svegajNnAsnJwG+/Kd0JNm/e\nDC8vL/j7+2PDhg1o2LChrPXNZjPefPNN9OrVC6GhobLWJrJa+/dX+McMZJXp0aMH/vd//xdarRY3\nbtyQvb5Go8Hx48dRXFwMDw8P7Ny5U/YeVOvpp4GPPpKmrxVy7do1jBw5EnPmzMGmTZswd+7ch17a\nZUlz587F9evXsXz5cq6oJqoMIYA5cyr8EgayCo0dOxaBgYF4/fXXYTKZZK/fpEkTREdH48svv8S4\nceM4Wr7bpEnS68i7dsleevPmzfDw8ICLiwtSU1PRq1cv2XsApM1lVq9ejU2bNtW6lflEFpOYCFy/\nXuGXMJBV6tNPP0W9evUQEhKi2Gr0IUOG4MSJEyguLoanpyd2KRBCqvPEE9L1yB98ID3jlcG9o+LF\nixfLtg3mvX799Ve89dZbiIuLg6OjoyI9EFmdkhLpd8b8+RV+GQNZpezs7LBx40ZkZGTg3XffVSyU\nGzdufN9o+datW4r0ohqvvw4UFwMWXvwmhMBPP/0ET09PxUfFAHDo0CH8/e9/x7p16/DCCy8o1geR\n1VmzBnjmGUCrrfjravQYiypQsLRVycvLExqNRvzjH/9Q/FzZvLw8MXnyZOHo6CiWLFki2+lBqvTz\nz0I895wQRqNF7n7//v2iT58+wt3dXezbt88iNapCr9cLR0dHsXXrVqVbIbIuhYVCuLoK8euvQoiK\ns48jZJVr3Lgxtm3bhm3btiE8PFzxXr7++mvs3LkTu3fvxnPPPYfvvvsOJSUlivalCH9/adeu776r\n0bvNyMhAUFAQ3njjDYwfPx5paWmKjooBID09HUOHDkVUVBSGDBmiaC9EVicqCnB3BypzXoFcTxLu\npWBpq5Sbmys6d+4s5s6dq/hIudS+fftE7969hbu7u4iLi1NNX7I5dEiIVq2EyM9/7LvKysoSY8eO\nFQ4ODmLx4sWqmX04fvy4cHZ2FuvWrVO6FSLrc+OGEE5OQhw/XvapirKPW2dakdzcXAwaNAgBAQFY\nsGCBKi43EUIgMTERH3zwARo1aoQFCxagb9++Srcln3HjpMuhvviiWt9+9epVzJ8/H6tWrUJoaCje\nf/991RzrqNfrERAQgGXLlmHEiBFKt0NkfcaPBxo0KLflboXZZ+HnBw+lYGmrdvXqVeHl5SXeeecd\nVY1ITSaTWL16tWjXrp3w9vYWMTExoqCgQOm2LO/6dWmUvGtXpb/FbDaLw4cPi7Fjx4omTZqI0NBQ\nYTAYLNhk1R04cEA4ODiI2NhYpVshsk6JiUK4ud03g1ZR9jGQrdCNGzfESy+9JKZMmSJKSkqUbqcc\nk8kkEhMTRUBAgGjevLmYNWuWyMzMVLoty0pIeOAP3r0KCgpEdHS00Gg0ws3NTXz66afi6tWrMjVZ\necnJycLBwYELuIiq68YNIVq3fuAT9Yqyj1PWVio/Px+BgYFo2rQpVq9eLfvWiZVx/vx5REVFISYm\nBl5eXggNDUVAQIBsZ/TKqoKp68zMTHz11VdYtWoVevTogdDQUPj7+yuyw9aj/Pjjj5gxYwbWrVuH\ngQMHKt0OkXUaPx546qkH/j7g8Yu1VFFREaZPn44jR44gPj4ebm5uSrf0QIWFhdi0aRMiIyORnZ2N\nV199FVqtFn369IG9vb3S7dWMGzcADw9g1SpgwACcPXsWW7ZsQVxcHE6ePImJEydiypQpqn2MSkpK\n8OGHH2LTpk2Ij49Hly5dlG6JyDolJgJvvw2kpQEPGCgxkGsxIQRWrFiBiIgIrFu3Dv369VO6pQqd\nOHECcXFx2LJlC06fPo3BgwcjMDAQQ4YMQZMmTZRur9pKSkqQ+fnncPrkE/g7OyP7xg0EBgYiMDAQ\nfn5+ePLJJ5Vu8aFu3ryJUaNG4c6dO9i4cSOaN2+udEtE1ikv768n5v37P/BLGMh1wC+//ILg4GCE\nh4cjJCRE6XYqxWAwICEhATqdDsnJyXjxxRcRGBiIAQMGwN3dXfVT2waDAQcPHkRCQgISExPh6OiI\n721s0KpjRzhs2KDKKel7nT59GlqtFgMHDsTnn39ee2YsiJRQwVR1KQZyHXHmzBkEBQXBx8cHy5Yt\ns6pfrgUFBfjll1+wZcsW7Nu3D9nZ2fDw8IC3tzc0Gg00Go2iIW0wGKDX66HX63Hs2DHo9XoUFRXh\nxRdfxNChQxEYGChNR5c+Q/7+e0DlZwTv2LEDo0ePxscff4wpU6Yo3Q6RdXvEVHUpBnIdcuvWLQQH\nByM/Px8bN26Eg4OD0i1VS35+PlJSUsoFYGlIe3l54dlnn4WLiwucnZ3L3jZr1qza12bfuXMHBoMB\nBoMBly5dKnubkZFRFr4ajabcE4S2bds+uF5iIvDWW0B6eoU/mEoRQmDp0qVYtGgR1q9fDx8fH6Vb\nIrJupWtIVq9+6FR1KQZyHVNSUoKwsDCsXbsW8fHx8PT0VLqlGlEa0sePH8fFixfLgrM0PO/cuYOW\nLVvCxcUFLVq0gL29Pezs7GBvbw9bW1uYzWaYTKayW35+/gO/19nZuezWqVOnisP3YR5zwxBLKSoq\nQkhICPR6PeLj4+Hq6qp0S0TWrxJT1aUYyHXU2rVr8fbbb2POnDmYOXMmbG1tlW7Jou4e5V67dq1c\n+BqNRtja2sLOzq4spJ9++umyEXbTpk1rduczFU5d6/V6jBs3Dp07d0ZMTIwqL5UjsjqVnKouxUCu\nw86dO4eJEyeisLAQ3377LTp16qR0S3XH1q3A9OmKT10XFRXhk08+wddff40lS5Zg1KhRqth2lcjq\nVWGqulRF2af+ZaD0WNq1a4edO3di9OjR6N27Nz777LO6eTqTEoYOlX5Ip00DFHryqdfr4e3tjfT0\ndBw/fhzBwcEMY6KaIAQweTLw8suVDuNH4Qi5DuFoWQG3bwO9egFjxgAzZ8pWlqNiIguLiAC2bAGS\nk4H69Sv9bRwhEwCOlhXx1FNAXBywaBHwyy+ylCwdFaelpSE1NZWjYqKalpAgneC0eXOVwvhROEKu\no86dO4cJEyagqKiIo2U57N4NjBwJHDwIWGj7zNJRcVRUFJYsWcIgJrKE338H+vQB4uOBHj2q/O0c\nIdN92rVrh127dmH06NHo1asX3nnnHVy5ckXptmqv/v2BDz+UXm8qKKjRuzabzVizZg3c3d3LXise\nPXo0w5iopt28CQQFAfPnVyuMH4WBXIfVq1cP06dPx6lTp2BjY4POnTsjPDwct27dUrq12unttwEv\nL+maxRqYHRJCICkpCRqNBkuXLsU333yDuLg4ODs710CzRFSO2Qy8+SYwcCAwaZJFSjCQCY6Ojli6\ndCmOHTuGs2fPomPHjli6dCmKioqUbq12sbEBvvoKyMqSXlN+DIcOHUL//v0xY8YMhIWFlX1MRBYS\nHi7tL/D55xYrwUCmMm5ubli9ejW2b9+OHTt24LnnnsOqVau48KsmPfmktBBk+XLg55+r/O2nTp3C\n8OHD8dprr+HNN9/EiRMn8Morr3B6msiSNm+WNvnZtAl44gmLlWEg0308PT2RkJCA1atX46uvvsIL\nL7yAhIQELsKrKa1bAxs2AGPHApmZlfqW7OxsTJgwAX379kWvXr2QmZmJiRMnqv5ELCKrd+IEMHWq\nFMqOjhYtxUCmh+rTpw/279+PiIgIzJkzB97e3oiJicHt27eVbs369eolXccYEAA8ZDGdEAKHDx/G\n2LFj0bVrV7Rs2RKZmZl4//330aBBA5kbJqqDcnKAYcOAf/0L0GgsXo6XPVGlmM1mbNu2DZGRkTh4\n8CDGjBmDadOmoWPHjkq3Zt3+53+kLTZ37wYaNQIA3L59G2vXrkVkZCTy8vIQEhKC8ePHo3nz5go3\nS1SH5OVJlzeNHg3Mnl1jd8u9rKlGZWVlISoqCjExMejatStCQ0MxbNgwTp9WhxDSdFhWFjKXLMGX\n0dFYtWoVevXqhdDQUPj5+aFePU5kEcnqzh3A318aFS9ZIi3IrCEMZLKIoqIibNq0CZGRkbhw4QKm\nTp2KSZMmoWXLlkq3ZjVMJhMS4uLgEBqK3Fu3oJ85E5OnTuWxiERKMZmA114DGjQAfvgBqOEniaPl\n1AAAEKtJREFUxAxksrjU1FR8+eWX2LBhA3x8fBAUFISAgAA4OTkp3ZrqGI1G7N27FzqdDj/99BPa\ntm2LtydNwmsrV6Jet241/oyciCpJCOkwmHPnpGMVLbCimoFMsrl58ya2bNkCnU6H7du3w93dHVqt\nFlqtFp07d66zl+fk5eUhKSkJOp0OSUlJaN++PbRaLV5++WV06dJF+qIbNwAfHyA4GJgzR9mGieqi\nsDDpcsS71nTUNAYyKaKoqAjJycllAW1vb4/AwEBotVr07t0b9vb2SrdoUefOnSv7ux89ehR9+/aF\nVqtFQEAAXFxcHvxNOTlA377SOcoyng5FVOdFREhT1MnJFr28SfZADg8PxzfffAMHBwcAwIIFCzB4\n8OBKN0W1jxACaWlp0Ol00Ol0OHv2LPz8/NCjRw9oNBq88MILaNiwodJtVpvZbMbp06eh1+tx7Ngx\n7NixA1euXMGwYcOg1WoxaNAgPP3005W7s+xsKZRnzpS22yQiy1q0CIiJAfbsASy89azsgTxv3jw0\natQIs2bNqlZTVPvl5ORg27ZtOHr0KPR6PTIyMtC2bVt4e3tDo9GoOqTvDV+9Xo+UlBQ0b94cGo0G\n3t7e8PHxwUsvvQRbW9vqFfnjD6BfP+D996XRMhFZxmefAVFRUhi3amXxcooEcsOGDfHee+9Vqymq\ne4xGIzIyMsqFXEZGBlxdXaHRaODh4QEXFxc4OzuXvX3mmWcs9pp0cXExLl++DIPBgEuXLsFgMODM\nmTMPDF+NRoO//e1vNX+d8PnzwIABwIwZwLvv1ux9ExGwYIE0Mt69W9pBTwaKBPK3336Lxo0bw9vb\nG4sXL0aTJk0q3RQRUD6kMzIyYDAYygWk2WyGs7NzuZB2cXFBs2bNYG9vDzs7u3I3s9kMk8lUdjMa\njbh9+3a5+yx9e/PmTTg6Opa779InBxYJ34f54w8plENDgQqe4BJRFX38sfSa8e7dwMPWdFiARQLZ\n19cXly9fvu/zERER6N69e9nrx2FhYTAYDIiOjr6vqblz55Z93K9fP/Tr1686rVAdlZ+ff19IX7p0\nCdevX0dJSUm54DWZTLC1tb0vpBs0aHBfoDs7O6NFixbVn26uadnZ0pFvb7wBzJvHS6KIHofZLO28\nlZgI7Nxp8deM9+zZgz179pR9PG/ePOVWWWdlZSEwMBDp6enlC3OETFR5ubnA0KHSzkGRkQB3RSOq\nOqMRmDgROH0aSEgAFNiOtqLss8iefAaDoez92NhYeHh4WKIMUd3h5CQtOjl/XtpF6M4dpTsisi4F\nBUBQEHD9ujQyVuHe8BYZIY8ZMwapqamwsbGBm5sboqKi7tuxiSNkomooLpaObczJAXQ64J61GUT0\nANeuSSerubsDX38NKLgHAjcGIapNzGZg1izpWf62bbIuSCGyOn/8IR0U8cor0uYfCq/BkH3Kmogs\nqF494PPPpS02e/QA9HqlOyJSp4MHpbPHQ0KA+fMVD+NHYSATWSMbG2m/6yVLgMGDgR9/VLojInX5\n5hvpNeOvvrKa6/g5ZU1k7dLTgZdflm6LFnEFNtVtxcXStrM7dwLx8cBzzyndUTmcsiaqzTw8gKNH\npWAeMgT4z3+U7ohIGbm5wKBBwIULwOHDqgvjR2EgE9UGzZoBW7cCXl7Aiy8CaWlKd0QkL70e6NZN\n2gM+Ph5o3FjpjqqMU9ZEtc2aNdJrZpGR0jXLRLXd6tXSlQdRUdJqahXjZU9Edc1vv0m/mIKDgX/+\nE1DLNqBENclkAv7rv6Rr8uPigC5dlO7okRjIRHXR1avAiBHSIq+YGKBNG6U7Iqo5588D48YBDRoA\na9cCTZsq3VGlcFEXUV3k4ADs2CG9pva3vwHR0QCfBJO1M5uBL76Q1koMGyYdEmElYfwoHCET1QXp\n6dJowsEBWLmSo2WyTufPAxMmAIWFwLffAp06Kd1RlXGETFTXeXgAhw4BffpwtEzW5+5R8dChwL59\nVhnGj8IRMlFdw9EyWZNaMCq+G0fIRPSX0tFy794cLZN61ZFR8d04Qiaqy9LSgPHjOVomdTl/Hpg4\nUTr3uxaMiu/GETIRPZinJ0fLpB5ms7ShzYsvStvA1oFR8d04QiYiyd2j5RUrgPbtle6I6pJ//xsI\nDa2Vo+K7cYRMRI9WOlru3x/o3l06Q9ZgULorqu0uXgQmTZKuABg2rM6Niu/GQCaiv9jbA7NnA7//\nDjz9tLQV4X//N3DzptKdUW1z/Trwj38AXbtKszKZmdJ+1HV4m1cGMhHdr3lz4LPPgJQUaZTcoYP0\ncWGh0p2RtSsoABYskI5GzM+XLsNbsKDW7Lb1OBjIRPRwzz4r7YO9Zw+wfz/QsaP0scmkdGdkbYxG\n4KuvpP9DqanAgQPSxy4uSnemGlzURUSVd+gQMGcOcOUKEBEBvPwyYGOjdFekZmYzsHEj8NFHgJub\nNBrWaJTuSjE87YmIao4QQFIS8MEHwJNPAgsXSgdYEN1NCOlwkw8+AOrVk/6fDByodFeKYyATUc0z\nm4F164CwMGkaMjwc6NaNI+a6Tgjp5Y25c6UV1BERwKuv8v/F/2MgE5HlFBdLu3wtXgw0ayZdS/rG\nG8BTTyndGcnpzz+BNWukjT1u3wbef1/ag9rOTunOVIWBTESWZzYD27ZJv5APHgTGjAGmTZNGz1R7\nnToFfPkl8OOPgI+P9IRs4EBpmpruw41BiMjy6tWTtjvcsgU4ehSoX1/aktPPD4iL48rs2sRolBZq\nDRgg3Ro3llZOx8YCvr4M42riCJmILKeoCNi0SRo1X7gATJ0q7crUsqXSnVF15OQAX38tvUTRoYM0\nGh4+HHjiCaU7sxocIRORMurXB4KDpUU+W7YA2dlA587Sa8y//sqDLKyBEMCuXcDf/y4d3XntGrB9\nO5CcDLz+OsO4BnGETETyunkTWLVKGjXb2kqvNQcFSTs3kToIAZw8Ceh0wPffS1uqhoYCo0cDjRop\n3Z1V46IuIlIfIaQdwDZskH7xP/00oNVKt549uTpXbkajdLCDTifdjEbpsXj9dWktAC9bqhEMZCJS\nNyGA3377Kwyys4GhQ6VA8PfnqMxSbt6UNnnR6YCffwbatfvrSVHXrgxhC2AgE5F1uXABSEiQgmL/\nfqBXLykkAgOBNm2U7s66nT8vvZ6v0wFHjkjHHmq10tGHrVop3V2tx0AmIut165a0iEinA7ZulQJZ\nq5VG0C+8IC0co4crLJRO7UpMlP4NL18GAgKkf0NfX6BhQ6U7rFMYyERUO5hM0qYjOp20CcmZM9Jh\n9t7e0oEFGo20EriuhnRhIZCWBuj1wLFj0tvMTGnBnL+/FMIvvVSnzxxWGgOZiGqn27elACoNH73+\nr5DWaP4K6toY0qXhe/ffPTNT2hnt7iconp7SISCkCgxkIqo77twBjh//K6T0euD06b9CujSg27QB\nnJ2lS3rUrLgYuHRJOqihdPR798i39O/E8LUKDGQiqtvuDemTJ6Vdp3JzgebNpcVMrVtLt9L3S9+6\nuEiXZNX0imMhpAMZSsM2J0d6e/f7OTnA9evSzmatWgHPP8/wtXIMZCKiBzGZpFC+OwDvDcdLl6RR\nasOG0q1RI+lW+n6DBtLezaU3QDpoo/R2+7YUvPn50q30/T//lALVxeXhTwZatwYcHfmaby3CQCYi\nehxGI1BQcH+o5udLo28hpPAtKZFG0qXhbGMjHUP5oCBv2JCbn9RBDGQiIiIV4OESREREKsdAJiIi\nUgEGMhERkQowkImIiFSAgUxERKQCDGQiIiIVYCATERGpAAOZiIhIBRjIREREKsBAJiIiUgEGMhER\nkQowkImIiFSAgUxERKQCDGQiIiIVYCATERGpAAOZiIhIBRjIREREKsBAJiIiUgEGMhERkQowkImI\niFSAgUxERKQCDGQiIiIVYCATERGpAAOZiIhIBRjIREREKsBAJiIiUgEGMhERkQowkImIiFSAgUxE\nRKQCDGQiIiIVYCATERGpAAOZiIhIBRjIREREKsBAJiIiUgEGMhERkQowkImIiFSAgUxERKQCDGQi\nIiIVYCATERGpAAOZiIhIBRjIREREKsBAJiIiUgEGMhERkQowkImIiFSAgUxERKQCDGQiIiIVYCAT\nERGpAAOZiIhIBRjIREREKsBAJiIiUgEGMhERkQowkImIiFSAgUxERKQCDGQiIiIVYCATERGpAAOZ\niIhIBRjIREREKsBAJiIiUgEGMhERkQowkImIiFSAgUxERKQCDGQiIiIVYCATERGpQLUDeePGjXj+\n+edha2uL3377rdyfLViwAB06dECnTp2wffv2x25SLfbs2aN0C/QY+PhZLz521o2PX+VUO5A9PDwQ\nGxsLHx+fcp8/efIk1q9fj5MnTyIpKQmhoaEwm82P3aga8D+VdePjZ7342Fk3Pn6VU+1A7tSpEzp2\n7Hjf5+Pj4zFy5EjY29vD1dUV7du3x5EjRx6rSSIiotquxl9DvnTpElq3bl32cevWrZGTk1PTZYiI\niGoVu4r+0NfXF5cvX77v8/Pnz0dgYGCli9jY2Nz3ua5duz7w82o3b948pVugx8DHz3rxsbNufPwk\nXbt2feifVRjIO3bsqHKxVq1aITs7u+zjixcvolWrVvd9XWpqapXvm4iIqLaqkSlrIUTZ+1qtFuvW\nrUNxcTHOnz+P06dPo1u3bjVRhoiIqNaqdiDHxsaiTZs2OHToEAICAjBkyBAAgLu7O0aMGAF3d3cM\nGTIEkZGRVjk1TUREJCcbcffwloiIiBTBnbqqKDw8HK1bt4aXlxe8vLyQlJSkdEv0CElJSejUqRM6\ndOiARYsWKd0OVZGrqys8PT3h5eXFl79UbsKECXBycoKHh0fZ565fvw5fX1907NgRfn5+yMvLU7BD\ndWMgV5GNjQ1mzZqFlJQUpKSkYPDgwUq3RBUoKSnBW2+9haSkJJw8eRJr167FqVOnlG6LqsDGxgZ7\n9uxBSkoK9zRQufHjx983SFm4cCF8fX2RmZmJgQMHYuHChQp1p34M5GrgLL/1OHLkCNq3bw9XV1fY\n29vjjTfeQHx8vNJtURXxZ8469OnTB02bNi33OZ1Oh7FjxwIAxo4di7i4OCVaswoM5GpYvnw5unbt\niokTJ3L6ReVycnLQpk2bso+5UY31sbGxwaBBg+Dt7Y2VK1cq3Q5VUW5uLpycnAAATk5OyM3NVbgj\n9WIgP4Cvry88PDzuu+l0OoSEhOD8+fNITU2Fs7Mz3nvvPaXbpQpwhb/1279/P1JSUvDzzz/jiy++\nwN69e5VuiarJxsaGP5MVqHBjkLqqshuiTJo0qUo7lpH87t2oJjs7u9zWrqR+zs7OAAAHBwcMHz4c\nR44cQZ8+fRTuiirLyckJly9fRsuWLWEwGODo6Kh0S6rFEXIVGQyGsvdjY2PLrSYk9fH29sbp06eR\nlZWF4uJirF+/HlqtVum2qJJu376N/Px8AEBBQQG2b9/Onzkro9Vq8f333wMAvv/+e7z88ssKd6Re\nHCFX0ezZs5GamgobGxu4ubkhKipK6ZaoAnZ2dlixYgX8/f1RUlKCiRMnonPnzkq3RZWUm5uL4cOH\nAwBMJhOCg4Ph5+encFf0MCNHjkRycjKuXbuGNm3a4J///CfmzJmDESNGIDo6Gq6urtiwYYPSbaoW\nNwYhIiJSAU5ZExERqQADmYiISAUYyERERCrAQCYiIlIBBjIREZEKMJCJiIhUgIFMRESkAv8H5pc+\nC9yTE8oAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10fcc9f10>"
       ]
      }
     ],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def calc_R((xc, yc), x, y):\n",
      "    \"\"\"Euclidean distance of each 2D point from the center (xc, yc)\"\"\"\n",
      "    return np.sqrt((x-xc)**2 + (y-yc)**2)\n",
      "\n",
      "def f(center, x,y):\n",
      "    \"\"\"Calculate the distance between the\n",
      "    data points and the circle centered at center=(xc, yc)\n",
      "    \"\"\"\n",
      "    Ri = calc_R(center, x,y)\n",
      "    return Ri - Ri.mean()\n",
      "\n",
      "def estimate_center(x,y):\n",
      "    return x.mean(), y.mean()\n",
      "\n",
      "Xs = np.array(intersect_x)\n",
      "Xs += np.random.normal(scale=0.01, size=Xs.shape)\n",
      "Ys = np.array(intersect_y)\n",
      "Ys += np.random.normal(scale=0.01, size=Ys.shape)\n",
      "\n",
      "# Use mean in x and y as starting point for the fit\n",
      "center_estimate = estimate_center(Xs, Ys)\n",
      "#center_fit, ier \n",
      "res = optimize.leastsq(f, center_estimate, args=(Xs, Ys), full_output=True)\n",
      "center_fit = res[0]\n",
      "pcov = res[1]\n",
      "\n",
      "# Gotta rescale the covariance matrix by this factor\n",
      "# not sure why they divide by it in the first place\n",
      "s_sq = (f(center_fit, Xs,Ys)**2).sum()/(len(Xs)-len(center_fit))\n",
      "pcov = pcov * s_sq\n",
      "\n",
      "R_fit = calc_R(center_fit, Xs,Ys).mean()\n",
      "\n",
      "#center_fit, ier = optimize.leastsq(full, (3.,0.,3), args=(Xs, Ys))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 166
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f(center_estimate, Xs, Ys)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 167,
       "text": [
        "array([ 0.51426879, -0.51310497, -0.52524274,  0.52407892])"
       ]
      }
     ],
     "prompt_number": 167
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Xs_=[]\n",
      "Ys_=[]\n",
      "\n",
      "err_xs = []\n",
      "err_ys = []\n",
      "\n",
      "for n in xrange(1000):\n",
      "    Xs = np.array(intersect_x)\n",
      "    Xs += np.random.normal(scale=0.01, size=Xs.shape)\n",
      "    Ys = np.array(intersect_y)\n",
      "    Ys += np.random.normal(scale=0.01, size=Ys.shape)\n",
      "\n",
      "    # Use mean in x and y as starting point for the fit\n",
      "    center_estimate = estimate_center(Xs, Ys)\n",
      "    \n",
      "    res = optimize.leastsq(f, center_estimate, args=(Xs, Ys), full_output=True)\n",
      "    center_fit = res[0]\n",
      "    pcov = res[1]\n",
      "\n",
      "    s_sq = (f(center_fit, Xs,Ys)**2).sum()/(len(Xs)-len(center_fit))\n",
      "    pcov = pcov * s_sq\n",
      "    p = np.sqrt(pcov.diagonal())\n",
      "    err_x = p[0]\n",
      "    err_y = p[1]\n",
      "    Xs_.append(center_fit[0])\n",
      "    err_xs.append(err_x)\n",
      "    \n",
      "Xs_ = np.array(Xs_)\n",
      "err_xs = np.array(err_xs)\n",
      "fig, ax = plt.subplots(ncols=2, figsize=(2*4,4))\n",
      "ax[0].hist(Xs_-curvature)\n",
      "ax[1].hist((Xs_-curvature)/err_xs, range=(-4,4))\n",
      "# this doesn't look right ... can't do pulls?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 168,
       "text": [
        "(array([  36.,   39.,   57.,   90.,  147.,  161.,  114.,   54.,   37.,   17.]),\n",
        " array([-4. , -3.2, -2.4, -1.6, -0.8,  0. ,  0.8,  1.6,  2.4,  3.2,  4. ]),\n",
        " <a list of 10 Patch objects>)"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAEACAYAAACNuW5TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtQVPf5P/D3ErA2X0HRyEJY081wCSzgssZAkiktFher\nqZSoQ4O5bAX7h6YXY2ow6aQRp8pGm7GaxOZSTUgzFe1MCzRFBnNZtCZmTcAkFVOIXSLX/cUgihqD\nwOf3B3ELclt2z+5+YN+vGWfw7Nnn+Rx2n304Z8/5HJUQQoCIiIh8LsDXAyAiIqJ+bMpERESSYFMm\nIiKSBJsyERGRJNiUiYiIJMGmTEREJIlRm/KVK1eQmpqK5ORk6HQ6PP744wCAjo4OGI1GxMbGIjMz\nE52dnY7nFBUVISYmBnFxcaiqqvLs6InIaXl5eVCr1UhKSnIss1qtSElJgcFgwB133IHjx487HmMt\nE/mAGMOlS5eEEEJcvXpVpKamiiNHjogNGzaIp59+WgghhNlsFgUFBUIIIU6ePCn0er3o7u4WNptN\nREVFid7e3rFSEJEXHD58WNTU1IjExETHsu9///uisrJSCCFERUWFSE9PF0Kwlol8ZczD1zfeeCMA\noLu7G729vQgNDUV5eTlMJhMAwGQyobS0FABQVlaG3NxcBAUFQavVIjo6Glar1YN/UhCRs9LS0hAa\nGjpoWUREBM6fPw8A6OzsRGRkJADWMpGvBI61Ql9fH+bNm4fTp09jzZo1SEhIgN1uh1qtBgCo1WrY\n7XYAQGtrK+68807HczUaDVpaWjw0dCJyl9lsxne/+138+te/Rl9fH9577z0ArGUiXxlzTzkgIAAn\nTpxAc3MzDh8+jHfeeWfQ4yqVCiqVasTnj/YYEflWfn4+du3ahTNnzmDHjh3Iy8sbcV3WMpHnjbmn\nfM306dNxzz334MMPP4RarUZ7ezvCw8PR1taGsLAwAEBkZCSampocz2lubnYcDhsoOjoap0+fVmD4\nRJNbVFQUPvvsM4/Ft1qtePPNNwEAK1aswOrVqwGwlomU5mwtj7qnfPbsWceZ1V999RUOHToEg8GA\nrKwsFBcXAwCKi4uRnZ0NAMjKykJJSQm6u7ths9nQ0NCAlJSUIXFPnz4NIYRX/z311FPMyZwTLqen\nG150dDSqq6sBAG+//TZiY2OlrGVv/O4nSw5ui5w5nK3lUfeU29raYDKZ0NfXh76+Pjz44IPIyMiA\nwWBATk4O9uzZA61WiwMHDgAAdDodcnJyoNPpEBgYiN27d/OQF5EkcnNzUV1djbNnz2LOnDnYvHkz\nXnrpJTz88MP4+uuv8e1vfxsvvfQSANYyka+M2pSTkpJQU1MzZPnMmTMdh7yu98QTT+CJJ55QZnRE\npJh9+/YNu/z9998fdjlrmcj7/GZGr/T0dOZkzgmXk/p543c/WXJ4Kw+3xTNUQgjh9aQqFXyQlmjC\nkb1WZB8fkSycrRW/2VMmIiKSHZsyERGRJNiUiYiIJMGmTEREJAk2ZSIiIkmwKRMREUmCTZmIiEgS\nbMpERESSYFMmIiKSBJsyERGRJNiUiYiIJMGmTEREJAk2ZSIiIkmwKRMREUmCTVkyISEzoVKp3PoX\nEjLT15tBEsrLy4NarUZSUtKg5c8++yzi4+ORmJiIgoICx/KioiLExMQgLi4OVVVV3h4ukV8K9PUA\naLCurnMA3Ls/bVeXSpnB0KSyatUq/OIXv8BDDz3kWPbOO++gvLwcH3/8MYKCgvDFF18AAOrq6rB/\n/37U1dWhpaUFCxcuRH19PQIC+Hc8kSexwoj8RFpaGkJDQwct++Mf/4jHH38cQUFBAIDZs2cDAMrK\nypCbm4ugoCBotVpER0fDarV6fcxE/oZNmciPNTQ04PDhw7jzzjuRnp6ODz74AADQ2toKjUbjWE+j\n0aClpcVXw6RxUuJrMH4l5hs8fE3kx3p6enDu3DkcO3YMx48fR05ODv773/8Ou65KNfzXIps2bXL8\nnJ6ejvT0dA+MlMZDia/BBsfjV2LjZbFYYLFYxv08NmUiP6bRaLBs2TIAwB133IGAgACcPXsWkZGR\naGpqcqzX3NyMyMjIYWMMbMpE1O/6P1ALCwudeh4PXxP5sezsbLz99tsAgPr6enR3d+Omm25CVlYW\nSkpK0N3dDZvNhoaGBqSkpPh4tESTH/eUifxEbm4uqqur8eWXX2LOnDnYvHkz8vLykJeXh6SkJEyZ\nMgWvvfYaAECn0yEnJwc6nQ6BgYHYvXv3iIeviUg5KiGEcl88OJtUpYIP0k4I/R987v5u+PudLGSv\nFdnH56+U+RwZFJGvs5ucrRUeviYiIpIEmzIREZEk2JSJiIgkMWpTbmpqwoIFC5CQkIDExETs2rUL\nQP8lEBqNBgaDAQaDAQcPHnQ8h/PlEhERuWbUE73a29vR3t6O5ORkXLx4EbfffjtKS0tx4MABBAcH\nY/369YPWr6urw8qVK3H8+PFR58vlySEj44leNJDstSL7+PwVT/SSjyIneoWHhyM5ORkAMG3aNMTH\nxzum2hsuOOfLJSIicp3T3yk3NjaitrYWd955J4D+273p9Xrk5+ejs7MTAOfLJSIicodTTfnixYtY\nsWIFdu7ciWnTpmHNmjWw2Ww4ceIEIiIi8Oijj474XE44QERE5JwxZ/S6evUqli9fjgceeADZ2dkA\ngLCwMMfjq1evxtKlSwHA5flyOYk9UT9XJ7Enoslh1BO9hBAwmUyYNWsWduzY4Vje1taGiIgIAMCO\nHTtw/Phx/OUvf3Gc6GW1Wh0nen322WdD9pZ5csjIlDlBIwhAj1sRgoNDceFCh5vjIHfJXiuyj89f\n8UQv+ThbK6PuKR89ehSvv/465s6dC4PBAADYunUr9u3bhxMnTkClUuHWW2/Fiy++CIDz5cqjB+4W\nJG/VRkTkfZz7WjJKXRLFy6omB9lrRfbx+SvuKctHkT1lIiLyrJCQmejqOufrYZAk2JSJiHyovyEr\nvRfKr58mKs59TUREJAk2ZSIiIkmwKRP5iby8PKjVaiQlJQ157JlnnkFAQAA6Ov53GRxvLkPkfWzK\nRH5i1apVqKysHLK8qakJhw4dwne+8x3Hsrq6Ouzfvx91dXWorKzE2rVr0dfX583hEvklNmUiP5GW\nlobQ0NAhy9evX49t27YNWsabyxD5BpsykR8rKyuDRqPB3LlzBy3nzWWIfIOXRBH5qcuXL2Pr1q04\ndOiQY9lokxuMNDsf57EnGsrVeezZlIn81OnTp9HY2Ai9Xg+g/wYyt99+O95//32Xby5DRP2u/wO1\nsLDQqefx8DWRn0pKSoLdbofNZoPNZoNGo0FNTQ3UajWysrJQUlKC7u5u2Gw2NDQ0ICUlxddDJpr0\n2JSJ/ERubi7uvvtu1NfXY86cOXjllVcGPT7w8PTAm8ssXryYN5ch8hLekEIyvCEFDSR7rcg+volA\n+ZtHAMp8BgyOx9fZPc7WCveUiYiIJMGmTEREJAk2ZSIiIkmwKRMREUmCTZmIiEgSbMpERESSYFMm\nIiKSBJsyERGRJNiUiYiIJMGmTEREJAk2ZSIiIkmwKRMREUmCTZmIiEgSbMpERESSYFMmIiKSBJsy\nERGRJEZtyk1NTViwYAESEhKQmJiIXbt2AQA6OjpgNBoRGxuLzMxMdHZ2Op5TVFSEmJgYxMXFoaqq\nyrOjJyKn5eXlQa1WIykpybFsw4YNiI+Ph16vx7Jly3D+/HnHY6xlIu8btSkHBQVhx44dOHnyJI4d\nO4bnn38ep06dgtlshtFoRH19PTIyMmA2mwEAdXV12L9/P+rq6lBZWYm1a9eir6/PKxtCRKNbtWoV\nKisrBy3LzMzEyZMn8dFHHyE2NhZFRUUAWMtEvjJqUw4PD0dycjIAYNq0aYiPj0dLSwvKy8thMpkA\nACaTCaWlpQCAsrIy5ObmIigoCFqtFtHR0bBarR7eBCJyRlpaGkJDQwctMxqNCAjo/xhITU1Fc3Mz\nANYyka84/Z1yY2MjamtrkZqaCrvdDrVaDQBQq9Ww2+0AgNbWVmg0GsdzNBoNWlpaFB4yEXnC3r17\nsWTJEgCsZSJfCXRmpYsXL2L58uXYuXMngoODBz2mUqmgUqlGfO5Ij23atMnxc3p6OtLT050ZCtGk\nZrFYYLFYvJ53y5YtmDJlClauXDniOqxlIue5WstjNuWrV69i+fLlePDBB5GdnQ2gf++4vb0d4eHh\naGtrQ1hYGAAgMjISTU1Njuc2NzcjMjJy2LgDC5mI+l3f1AoLCz2e89VXX0VFRQXeeustxzLWMpF7\nXK3lUQ9fCyGQn58PnU6HdevWOZZnZWWhuLgYAFBcXOxo1llZWSgpKUF3dzdsNhsaGhqQkpIy3m0h\nIi+prKzE9u3bUVZWhqlTpzqWs5aJfGPUPeWjR4/i9ddfx9y5c2EwGAD0XyaxceNG5OTkYM+ePdBq\ntThw4AAAQKfTIScnBzqdDoGBgdi9e/eoh7aJyHtyc3NRXV2Ns2fPYs6cOSgsLERRURG6u7thNBoB\nAHfddRd2797NWibyEZUQQng9qUoFH6SdEPo/+Nz93SgTg6+R78leK7KPbyJQpuaHRFU4Jl9ndzlb\nK5zRi4iISBJsykRERJJw6pIoIiLyZ4GKnlMQHByKCxc6FIs3mbApExHRGHqg5HfUXV08aXAkPHxN\nREQkCTZlIiIiSfDwtcJCQmaiq+ucr4dBREQTEJuywvobsjvfvfC7FiIif8XD10RERJJgUyYiIpIE\nmzIREZEk2JSJiIgkwaZMREQkCTZlIiIiSbApExERSYJNmchP5OXlQa1WIykpybGso6MDRqMRsbGx\nyMzMRGdnp+OxoqIixMTEIC4uDlVVVb4YMpHfYVMm8hOrVq1CZWXloGVmsxlGoxH19fXIyMiA2WwG\nANTV1WH//v2oq6tDZWUl1q5di76+Pl8Mm8ivsCkT+Ym0tDSEhoYOWlZeXg6TyQQAMJlMKC0tBQCU\nlZUhNzcXQUFB0Gq1iI6OhtVq9fqYifwNmzKRH7Pb7VCr1QAAtVoNu90OAGhtbYVGo3Gsp9Fo0NLS\n4pMxEvkTNmUiAgCoVKpRb2Sv5E3uiWh4vCEFkR9Tq9Vob29HeHg42traEBYWBgCIjIxEU1OTY73m\n5mZERkYOG2PTpk2On9PT05Genu7JIRNNCBaLBRaLZdzPUwkh3LmlkUtUKhV8kNYr+vcm3L1LlLu/\nG2ViTNbXaCJRulYaGxuxdOlSfPLJJwCAxx57DLNmzUJBQQHMZjM6OzthNptRV1eHlStXwmq1oqWl\nBQsXLsRnn302ZG95Mteyt7j/mTFsVIVjKh/P3943ztYK95SJ/ERubi6qq6tx9uxZzJkzB5s3b8bG\njRuRk5ODPXv2QKvV4sCBAwAAnU6HnJwc6HQ6BAYGYvfu3Tx8TeQF3FNWGPeUSUmy14rs45sIuKfs\nH5ytFZ7oRUREJAk2ZSIiIkmwKRMREUmCTZmIiEgSYzbl4Sax37RpEzQaDQwGAwwGAw4ePOh4jJPY\nExERuWbMs6+PHDmCadOm4aGHHnJc21hYWIjg4GCsX79+0LrXrm08fvy449rG+vp6BAQM7v2T+YxN\nnn1NSpK9VmQf30TAs6/9g2JnXw83iT2AYYNzEnsiIiLXufyd8rPPPgu9Xo/8/HzHPVg5iT0REZHr\nXJrRa82aNfjtb38LAHjyySfx6KOPYs+ePcOuO9IsQJwvl2goV+fLJe8JCZmJrq5zvh4GTVIuNeVr\nk9YDwOrVq7F06VIArk9iT0T9rv8DtbCw0HeDoWH1N2Slv68l6ufS4eu2tjbHz3//+98dZ2ZnZWWh\npKQE3d3dsNlsaGhoQEpKijIjJSIimuTG3FO+fhL7wsJCWCwWnDhxAiqVCrfeeitefPFFAJzEnoiI\nyB28IYXCeEkUKUn2WpF9fJ6g/CVMvCTKH/CGFERERBMM76dMIwh0+6uH4OBQXLjQodB4iIgmPzZl\nGkEP3D1c1dXF8wmIiMaDh6+JCEVFRUhISEBSUhJWrlyJr7/+Gh0dHTAajYiNjUVmZqZjkiAi8hw2\nZSI/19jYiJdffhk1NTX45JNP0Nvbi5KSEpjNZhiNRtTX1yMjIwNms9nXQyWa9NiUifxcSEgIgoKC\ncPnyZfT09ODy5cu4+eabUV5eDpPJBAAwmUwoLS318UiJJj82ZSI/N3PmTDz66KO45ZZbcPPNN2PG\njBkwGo2w2+1Qq9UAALVaDbvd7uOREk1+bMpEfu706dP4wx/+gMbGRrS2tuLixYt4/fXXB62jUqk4\nERCRF/DsayI/98EHH+Duu+/GrFmzAADLli3De++9h/DwcLS3tyM8PBxtbW2D5rwfiDeXIRrK1ZvL\ncEYvhU2mGb04K5jveaNWPvroI9x///04fvw4pk6dip/+9KdISUnB559/jlmzZqGgoABmsxmdnZ1D\nTvaazLU8Es7opUw8f3zfOLPNbMoKY1MeHGOyvs7e4q1a2bZtG4qLixEQEIB58+bhT3/6E7q6upCT\nk4MzZ85Aq9XiwIEDmDFjhk/GJxM2ZWXi+eP7hk3ZB9iUB8eYrK+zt8heK7KPzxPYlJWJ54/vG859\nTURENIGwKRMREUmCTZmIiEgSbMpERESSYFMmIiKSBJsyERGRJNiUiYiIJMGmTEREJAk2ZSIiIkmw\nKRMREUmCTZmIiEgSbMpERESSYFMmIiKSBJsyERGRJNiUiYiIJMGmTEREJIkxm3JeXh7UajWSkpIc\nyzo6OmA0GhEbG4vMzEx0dnY6HisqKkJMTAzi4uJQVVXlmVETkaI6OzuxYsUKxMfHQ6fT4f333x+1\nzonIM8ZsyqtWrUJlZeWgZWazGUajEfX19cjIyIDZbAYA1NXVYf/+/airq0NlZSXWrl2Lvr4+z4yc\niBTzq1/9CkuWLMGpU6fw8ccfIy4ubsQ6JyLPGbMpp6WlITQ0dNCy8vJymEwmAIDJZEJpaSkAoKys\nDLm5uQgKCoJWq0V0dDSsVqsHhk1ESjl//jyOHDmCvLw8AEBgYCCmT58+Yp0Tkee49J2y3W6HWq0G\nAKjVatjtdgBAa2srNBqNYz2NRoOWlhYFhklEnmKz2TB79mysWrUK8+bNw89+9jNcunRpxDonIs8J\ndDeASqWCSqUa9fHhbNq0yfFzeno60tPT3R0K0YRnsVhgsVi8mrOnpwc1NTV47rnncMcdd2DdunVD\nDlWPVuesZaKhXK5l4QSbzSYSExMd/7/ttttEW1ubEEKI1tZWcdtttwkhhCgqKhJFRUWO9RYtWiSO\nHTs2JJ6TaSckAAIQbvxz9/lyxSD3eON32NbWJrRareP/R44cEUuWLBFxcXHD1rm3xycbZWpD6Vrz\ndEzl4/kbZ7fZpcPXWVlZKC4uBgAUFxcjOzvbsbykpATd3d2w2WxoaGhASkqKKymIyEvCw8MxZ84c\n1NfXAwDefPNNJCQkYOnSpcPWORF5juqbDj6i3NxcVFdX4+zZs1Cr1di8eTN+/OMfIycnB2fOnIFW\nq8WBAwcwY8YMAMDWrVuxd+9eBAYGYufOnVi0aNHQpCoVxkg7YfUf4nNn29x9vlwxJuvr7C3eqpWP\nPvoIq1evRnd3N6KiovDKK6+gt7d3xDr39vhk4n6ND4mocDxPxFQ+nj++b5zZ5jGbsidM5kJmUx4c\nY7K+zt4ie63IPj4ACAmZia6ucwpHlbnheSImm7K7nK0Vt0/0IiKSWX9DVrpBEXkGp9kkIiKSBJsy\nERGRJHj4moiIvCxw1PktXBEcHIoLFzoUjekLbMpERORlPVD65LaursnxXT8PXxMREUmCTZmIiEgS\nbMpERESSYFMmIiKSBJsyERGRJNiUiYiIJMGmTEREJAk2ZSIiIkmwKRMREUmCTZmIiEgSbMpERESS\nYFMmIgBAb28vDAYDli5dCgDo6OiA0WhEbGwsMjMz0dnZ6eMREk1+bMpEBADYuXMndDqd4+49ZrMZ\nRqMR9fX1yMjIgNls9vEIiSY/NmUiQnNzMyoqKrB69WoI0X/3nvLycphMJgCAyWRCaWmpL4dI5BfY\nlIkIjzzyCLZv346AgP99JNjtdqjVagCAWq2G3W731fCI/Abvp0zk59544w2EhYXBYDDAYrEMu45K\npRrxpvSbNm1y/Jyeno709HTlB0k0wVgslhHraTQqce1YlRepVCr4IK1X9H9wubNt7j5fphhB6L+Z\nueuCg0Nx4UKHm+OYuLxRK0888QT+/Oc/IzAwEFeuXMGFCxewbNkyHD9+HBaLBeHh4Whra8OCBQvw\n6aefenx8ISEz0dV1TtGY7r+XB1KiNjwZzxMxZY/XH1PmvuJsrbApK4xNWfkYk/W94gxv10p1dTV+\n//vf4x//+Acee+wxzJo1CwUFBTCbzejs7Bxyspcnxud+DQ2J6GfxPBFT9nj9MWX+rHC2VvidMhEN\ncu0w9caNG3Ho0CHExsbi7bffxsaNG308MqLJj3vKCuOesvIxJut7xRmy14pKpUJJSYli8b71rW/h\n3nvvhdx7ZbLH80RM2eMBSnxdNpDSX53x8LWPsCkrH2OyvlecIXutqFQqBAfnKBbvq6/+gZ6eryB3\nA5A9nidiyh7PEzGVrT02ZR9hU1Y+xmR9rzhD9lpR+vvf//u/7+DSpTOKxpS/oUyMBiV3PE/E9E1T\n5iVRA3jmrE8iIiLnuNWUtVotQkJCcMMNNyAoKAhWqxUdHR34yU9+gs8//xxarRYHDhzAjBkzlBqv\nR/U3ZCX2DomIiMbPrbOvVSoVLBYLamtrYbVaAXC+XCIiIle5fUnU9cfIOV8uERGRa9zeU164cCHm\nz5+Pl19+GQDnyyUiInKVW98pHz16FBEREfjiiy9gNBoRFxc36PHR5sslIiKiwdxqyhEREQCA2bNn\n495774XVaoVarUZ7e7tjvtywsLBhn8tJ7ImGcnUSeyKaHFy+Tvny5cvo7e1FcHAwLl26hMzMTDz1\n1FN48803fTJfrhKUueaS1ykrHUPG94q3yFor1/A6ZRnjeSKm7PE8EXOCXadst9u/mQ4P6Onpwf33\n34/MzEzMnz8fOTk52LNnj+OSKCIiIhobZ/QagHvKcsaQ8b3iLbLWyjXcU5Yxnidiyh7PEzF9s6fM\nu0QRERFJgk2ZiIhIEmzKREREkmBTJiIikgSbMpGfa2pqwoIFC5CQkIDExETs2rULANDR0QGj0YjY\n2FhkZmais7PTxyMlmvzYlIn8XFBQEHbs2IGTJ0/i2LFjeP7553Hq1CneXIbIB9iUifxceHg4kpOT\nAQDTpk1DfHw8WlpaeHMZIh9gUyYih8bGRtTW1iI1NZU3lyHyAbfmviaiyePixYtYvnw5du7cieDg\n4EGPjX5zmU0Dfk7/5h+Rf3N1HnvO6DUAZ/SSM4aM7xVv8VatXL16FT/60Y+wePFirFu3DgAQFxcH\ni8XiuLnMggUL8Omnnw4ZH2f0ki2eJ2LKHs8TMTmjFxH5gBAC+fn50Ol0joYMAFlZWSguLgYAFBcX\nIzs721dDJPIb3FMegHvKcsaQ8b3iLd6olX/961/43ve+h7lz5zoOURcVFSElJQU5OTk4c+aM4+Yy\nM2bMGDI+7inLFs8TMWWP54mYvtlTZlMegE1Zzhgyvle8RdZauYZNWcZ4nogpezxPxAwC0KNgPHj2\n1o1ERESTVw+U/0NkbPxOmYiISBJsykRERJJgUyYiIpIEmzIREZEk2JSJiIgkwaZMREQkCTZlIiIi\nSbApExERSYJNmSQX6LhDkSv/QkJm+noDiIicxhm9SHLuzarT1eXcLDpERDLgnjIREZEk2JSJiIgk\nMSkOX588eRKPP74Fvb2uH+YM4J8nRETkY5OiKdfW1uLQoTO4cuVhl2NMmfInBUdEREQ0fh5pypWV\nlVi3bh16e3uxevVqFBQUeCLNIIGBWgC5Lj//hhuqAbyt1HCIJgVf1DKRP1P8oG1vby9+/vOfo7Ky\nEnV1ddi3bx9OnTqldBoXWJiTOd3LaPF+Tl+Sq5YtzCEdyyTK440czlG8KVutVkRHR0Or1SIoKAj3\n3XcfysrKlE7jAgtz+mVO965zHnits781Zblq2cIc0rFMojzeyOEcxZtyS0sL5syZ4/i/RqNBS0uL\n0mmInHTtOmfX/3V1dUGlUqGwsNCvJjFhLRN5n+LfKatU3p+sISAgAD097yAkZOmI61y58h9Mnfrh\niI93d3/siaHRpHCtsW/65t/4TcRJTJyt5dHqbrwuX/5/isUimogUb8qRkZFoampy/L+pqQkajWbQ\nOlFRUR5p3leutI76eHd3gxNRlBjXwBiFbj5fiTF4I8ZI2+nJcTj7u1VyDK68nt9EcOE9HxUV5XI+\ndzlby6dPv+GB7MP9rlz/3Tv/HvDGe2q4HJ74o03pmL56TcaTR97XxdlaVgkhXL+4dxg9PT247bbb\n8NZbb+Hmm29GSkoK9u3bh/j4eCXTEJGHsZaJvE/xPeXAwEA899xzWLRoEXp7e5Gfn88iJpqAWMtE\n3qf4njIRERG5xiuTS/71r39FQkICbrjhBtTU1Iy4nlarxdy5c2EwGJCSkuKVnJWVlYiLi0NMTAye\nfvppt3J2dHTAaDQiNjYWmZmZ6OzsHHY9d7fTmTH/8pe/RExMDPR6PWpra8edw5W8FosF06dPh8Fg\ngMFgwO9+9zu38uXl5UGtViMpKWnEdZTezrFyKr2NTU1NWLBgARISEpCYmIhdu3YNu54nXk+lPfPM\nMwgICEBHR4fisZ988kno9XokJycjIyNj0HfdStqwYQPi4+Oh1+uxbNkynD9/XvEczn42uULJz7OR\nOFOX7nK2Ltxx5coVpKamIjk5GTqdDo8//rjiOQbq7e2FwWDA0qVjnBgpvODUqVPiP//5j0hPTxcf\nfvjhiOtptVrx5Zdfei1nT0+PiIqKEjabTXR3dwu9Xi/q6upczrlhwwbx9NNPCyGEMJvNoqCgYNj1\n3NlOZ8b8z3/+UyxevFgIIcSxY8dEamqqS7nGm/edd94RS5cudTvXNYcPHxY1NTUiMTFx2Mc9sZ1j\n5VR6G9va2kRtba0QQoiuri4RGxvrlddTaWfOnBGLFi1StIYHunDhguPnXbt2ifz8fMVzCCFEVVWV\n6O3tFUIGou0lAAAF1UlEQVQIUVBQMGINu8PZz8PxUvrzbCRj1YgSnKkLJVy6dEkIIcTVq1dFamqq\nOHLkiOI5rnnmmWfEypUrx/z88MqeclxcHGJjY51aVyh0NN2ZnEpPjlBeXg6TyQQAMJlMKC0tHXFd\nV7fTmTEPHEdqaio6Oztht9tdyjeevIByrx8ApKWlITQ0dMTHPbGdY+UElN3G8PBwJCcnAwCmTZuG\n+Ph4tLYOvorAE9uptPXr12Pbtm0eix8cHOz4+eLFi7jppps8ksdoNCLgm7vTpKamorm5WfEc4/k8\nHA9vTfbiTI24y5m6UMKNN94IAOju7kZvby9mzvTMfALNzc2oqKjA6tWrx/z8kOreSCqVCgsXLsT8\n+fPx8ssvezyf0pMj2O12qNVqAIBarR7xg9Od7XRmzMOt4+6HizN5VSoV3n33Xej1eixZsgR1dXVu\n5XRlTJ74EB3Ik9vY2NiI2tpapKamDlrui+0cj7KyMmg0GsydO9ejeX7zm9/glltuQXFxMTZu3OjR\nXACwd+9eLFmyxON5lDJZJ3sZqS6U0NfXh+TkZKjVaixYsAA6nU7xHADwyCOPYPv27Y4/+Eaj2NnX\nRqMR7e3tQ5Zv3bp17GPo3zh69CgiIiLwxRdfwGg0Ii4uDmlpaR7L6cp1oyPl3LJly5DYI8Uf73a6\nMubr/xpz97pwZ54/b948NDU14cYbb8TBgweRnZ2N+vp6t/KORentHIuntvHixYtYsWIFdu7ciWnT\npg153Nvbeb3R3vdFRUWoqqpyLHP1SMJY9bxlyxZs2bIFZrMZjzzyCF555RWP5AH6t2vKlClYuXKl\nx3IozRcTN3naWHXhroCAAJw4cQLnz5/HokWLYLFYkJ6ermiON954A2FhYTAYDE5N1atYUz506JDb\nMSIiIgAAs2fPxr333gur1Tpqs3I3pzOTI4wnp1qtRnt7O8LDw9HW1oawsLBh1xvvdo53zNev09zc\njMjISKfiu5N34CHGxYsXY+3atejo6PDYISFPbOdYPLGNV69exfLly/HAAw8gOzt7yOO+2M7rjfS+\n//e//w2bzQa9Xg+gf2y33347rFbriO//8ea43sqVK93agx0rz6uvvoqKigq89dZbHsvhCa58nsls\nrLpQ0vTp03HPPffggw8+ULwpv/vuuygvL0dFRQWuXLmCCxcu4KGHHsJrr7027PpeP3w90l/Rly9f\nRldXFwDg0qVLqKqqUuzsvpFyzp8/Hw0NDWhsbER3dzf279+PrKwsl/NkZWWhuLgYAFBcXDzsG8nd\n7XRmzFlZWY4X/NixY5gxY4bjsLqrnMlrt9sdv2ur1QohhMcaMuCZ7RyL0tsohEB+fj50Oh3WrVs3\n7Dq+2E5nJSYmwm63w2azwWazQaPRoKamZtwNeSwNDf+bja+srAwGg0HR+NdUVlZi+/btKCsrw9Sp\nUz2SYyAlz09Q+vPMl5ypC3edPXvWcYXMV199hUOHDnnkfbV161Y0NTXBZrOhpKQEP/jBD0ZsyAC8\nc/b13/72N6HRaMTUqVOFWq0WP/zhD4UQQrS0tIglS5YIIYQ4ffq00Ov1Qq/Xi4SEBLF161aP5xRC\niIqKChEbGyuioqLczvnll1+KjIwMERMTI4xGozh37tyQnEps53BjfuGFF8QLL7zgWOfhhx8WUVFR\nYu7cuYqd4TlW3ueee04kJCQIvV4v7rrrLvHee++5le++++4TERERIigoSGg0GrFnzx6Pb+dYOZXe\nxiNHjgiVSiX0er1ITk4WycnJoqKiwiuvpyfceuutHjn7evny5SIxMVHo9XqxbNkyYbfbFc8hhBDR\n0dHilltucbwWa9asUTzHSJ9NSlDy82wk12pkypQpQqPRiL179yqeY7i6OHjwoKI5Pv74Y2EwGIRe\nrxdJSUli27ZtisYfjsViGfPsa04eQkREJAmpzr4mIiLyZ2zKREREkmBTJiIikgSbMhERkSTYlImI\niCTBpkxERCQJNmUiIiJJsCkTERFJ4v8DZtaUH9SMhHoAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1119ecbd0>"
       ]
      }
     ],
     "prompt_number": 168
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Ms = []\n",
      "Bs = []\n",
      "\n",
      "err_Ms = []\n",
      "err_Bs = []\n",
      "\n",
      "X = np.linspace(0,10,60)\n",
      "M = 4.\n",
      "B = 2.5\n",
      "uncertainty = 0.1\n",
      "\n",
      "for n in xrange(1000):\n",
      "    Ys = M*X + B\n",
      "    Ys += np.hstack([np.random.normal(scale=uncertainty, size=len(X)/3),\n",
      "                     np.random.normal(scale=uncertainty*5, size=len(X)/3),\n",
      "                     np.random.normal(scale=uncertainty, size=len(X)/3)])\n",
      "\n",
      "    uncertainties = np.hstack([np.ones(len(X)/3)*uncertainty,\n",
      "                               np.ones(len(X)/3)*5*uncertainty,\n",
      "                               np.ones(len(X)/3)*uncertainty])\n",
      "    errfunc = lambda p, x, y: ((p[0]*x + p[1]) - y) / uncertainties\n",
      "    \n",
      "    res = optimize.leastsq(errfunc, (0.,0.), args=(X, Ys), full_output=True)\n",
      "    m,b = res[0]\n",
      "    Ms.append(m)\n",
      "    Bs.append(b)\n",
      "    \n",
      "    pcov = res[1]\n",
      "\n",
      "    s_sq = (errfunc((m,b), X,Ys)**2).sum()/(len(X)-2)\n",
      "    pcov = pcov * s_sq\n",
      "    p = np.sqrt(pcov.diagonal())\n",
      "    err_M = p[0]\n",
      "    err_B = p[1]\n",
      "    err_Ms.append(err_M)\n",
      "    err_Bs.append(err_B)\n",
      "    \n",
      "Ms = np.array(Ms)\n",
      "err_Ms = np.array(err_Ms)\n",
      "Bs = np.array(Bs)\n",
      "err_Bs = np.array(err_Bs)\n",
      "\n",
      "\n",
      "fig, ax = plt.subplots(ncols=3, figsize=(3*4,4))\n",
      "ax[0].hist(Ms - M)\n",
      "ax[0].set_title(\"Bias of M\")\n",
      "ax[1].hist((Ms-M)/err_Ms, range=(-4,4))\n",
      "ax[1].set_title(\"Pull of M\")\n",
      "ax[2].hist((Bs-B)/err_Bs, range=(-4,4))\n",
      "ax[2].set_title(\"Pull of B\")\n",
      "\n",
      "print np.mean(Ms-M), np.std(Ms-M)\n",
      "print np.std((Ms-M)/err_Ms), np.std((Bs-B)/err_Bs)\n",
      "# conclusion: I do know how to do toy studies to get pulls etc"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.000220592860268 0.00452046793266\n",
        "1.0216160568 1.0375277177\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAsQAAAEKCAYAAAAVX0PaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9YlHW+//HXILj9EBIyBmX0YIAhikD11dpdN1zFaguy\n7KJwU47h7rns7OlYnuzH2Qo6Jbid6lK7rK5Wi9qrzN09gaeTHLOrcautg5W2bbRhLSa/y4AELVG5\nv3+wzor8cJi5Z+aG+/m4Lq4LZub+vN/33Pd75s099/0Zh2EYhgAAAACbCgt1AgAAAEAo0RADAADA\n1miIAQAAYGs0xAAAALA1GmIAAADYGg0xAAAAbI2GeBhavny5HnzwwaDGfOKJJ+R0OhUVFaW2trag\nxgZGumeffVazZ8/2/B0WFqa//vWvPo21dOlSxcTE6JJLLjErPQCnoGZHHhpiC0pISNBZZ52lyMhI\nxcTE6Oqrr1Z9fb3n/ieeeEK//OUvg5bP0aNHtXLlSr3++us6ePCgoqOje92/b98+hYWF6cILL+x1\n+4EDBzR69GhNnjw5aLkCoXRy7cbFxWnp0qU6dOhQ0OK/+eab2rFjhxobG/Xuu+/2uf/ZZ59VWFiY\nbr/99l63V1RUKCwsTEuXLg1WqoAlDIeaHTVqlCIjIxUZGanExEQ9+eSTQcvPTmiILcjhcOiVV15R\nR0eHmpqa5HQ69S//8i8hy6e5uVnfffedpk6dOujjvv32W3388ceev1944QWdf/75cjgcgU4RsIST\na/eDDz7Qe++9F9RPc7744gslJCTojDPOGDC/xMRE/fa3v9Xx48c9t5eVlWnKlCnUKmzH6jUrST/4\nwQ/U0dGhjo4O/f73v9eqVau0Z8+eoOVoFzTEFve9731PCxcuVHV1tee2f/zHf9S9994rSWpra9PV\nV1+t2NhYxcTEKCcnRw0NDZ7HPvvss0pMTFRUVJTOP/98vfDCC/3GOXLkiFasWKH4+HjFx8frtttu\nU1dXl2pqajyN8NixYzVv3rwBc128eLHKyso8fz///PNasmSJ+DJE2NGECRN0xRVX6OOPP9YXX3yh\nsLAwdXd3e+7PysrSxo0bhzxuY2OjcnNzde655yo5OVm//vWvJUkbN27Uz372M73zzjuKjIxUcXFx\nv8vHxcUpLS1N//u//ytJam1t1TvvvKPc3FxqFbZm1Zo9uS4zMjI0depU/eUvfxlyHhgcDbFFnSiA\nw4cP66WXXtKll17quc/hcHiO5BiGocLCQu3fv1/79+/XmWeeqV/84heSpEOHDulf//VfVVlZqYMH\nD+qdd95RRkZGv/EeeughVVVV6cMPP9SHH36oqqoqPfjgg5oyZYrnqO8333yjHTt2DJjzT3/6U23e\nvFmGYai6ulqdnZ2aNWuWKc8HMFycqN26ujpt27ZNmZmZ/TaaJ9fxUNx4442aNGmSmpqa9Lvf/U73\n3HOP3njjDRUWFurJJ5/UpZdeqo6ODt1///0D5rZ48WI999xzkqTNmzfrmmuu0fe+970h5wKMBFau\n2VNVVVWppqZGF1988ZDzwOBoiC3IMAwtWLBA0dHRGjt2rF5//XX927/9W5/HSFJMTIyuvfZanXHG\nGRozZozuuece7dy50/O4sLAwffTRR/r222/ldDqVmprab8wXXnhB9913n8aNG6dx48bp/vvv1/PP\nP98r1um4XC5dcMEFeu211/Tcc89pyZIlvqw+MGydXLuzZ89WVlaW7rnnHtPGr6ur0x//+EetWbNG\no0ePVnp6upYtW+Zpbr2t1WuvvVZut1sHDx7U888/r4KCAtNyBIaT4VCz7777rqKjoxUVFaVLLrlE\nS5YsUVJSkmk5ogcNsQU5HA5VVFSora1NR44c0fr163XZZZfpyy+/7PPYw4cP65/+6Z+UkJCgc845\nR5dddpm++eYbGYahs88+Wy+99JKefPJJTZgwQVdffbU+/fTTfmM2NjbqH/7hHzx/T5o0SY2NjUPO\ne8mSJXrmmWe0efNmLV68mI9gYSsn1+6+ffv0+OOPm3rktbGxUTExMTr77LM9t02aNKnXaVLeOOOM\nM3TVVVfpP/7jP9Ta2qpLL72UWoUtDYeaveSSS9TW1qaDBw+qublZf/7zn01t2tGDhtjiHA6Hrr32\nWo0aNUpvvfVWr9sl6ZFHHlFNTY2qqqr0zTffaOfOnTIMw/PmNn/+fG3fvl3Nzc1KSUnRz372s37j\nTJgwQfv27fP8vX//fk2YMGHI+V533XV69dVXlZiYKJfLNeTlgZHoxJvh4cOHPbc1NzcPeZwJEyao\ntbVVnZ2dntv279/vU60tWbJEjz76qG666aYhLwuMdFasWUmKjY3Vddddp//+7//2aXkMjIbYok40\ntIZheP57PXFx28kNb2dnp84880ydc845am1t7XVS/pdffqmKigodOnRIEREROvvsszVq1Kh+4+Xn\n5+vBBx/UgQMHdODAAT3wwANavHjxkPM+++yz9cYbb3guGgAgnXfeeYqPj9fzzz+v48ePa9OmTfr8\n88+HPM7EiRP1/e9/X3fffbeOHDmiP/3pT9q0aZNPTe1ll12mHTt2hHQGG8CqrFizkvT111/r5Zdf\n1vTp031aHgOjIbaonJwcRUZG6pxzztG9996r5557ztMQn3xi/4oVK/Ttt99q3Lhx+v73v68rr7zS\nc193d7cee+wxxcfH69xzz9Wbb76pJ554ot94v/zlL3XxxRdrxowZmjFjhi6++OJecx2f7kKCk++/\n8MILe809zFROgPT000/r4Ycf1rhx41RdXa0f/OAHnvtOvVhnsJp58cUXtW/fPk2YMEHXXXedHnjg\nAf34xz/ud5xTnXr/nDlzNHbsWK+WBezGKjV7YhaKyMhIpaamyul0av369SasIU7mMAY5cayurk5L\nlizRl19+KYfDoZ///Oe69dZbVVRUpF//+tc677zzJEmrV6/WlVdeKUkqKSnRpk2bNGrUKK1bt07z\n588PzpoA0HfffafLLrtMR44cUVdXl6655hqVlJSotbVVN9xwg2fOyy1btngaIWoWCA3qFbCOQRvi\n5uZmNTc3KyMjQ52dnbroootUXl6uLVu2KDIyss+3HVVXV2vRokXatWuXGhoaNG/ePNXU1CgsjAPR\nQLAcPnxYZ511lo4dO6Yf/vCH+s///E9t3bpV48aN06pVq7RmzRq1tbWptLSUmgVCjHoFrGHQKoqL\ni/PMWztmzBhNnTrVc2Vkf310RUWF8vPzFRERoYSEBCUlJamqqioAaQMYyFlnnSVJ6urq0vHjxxUd\nHa2tW7d6ptYqKChQeXm5JGoWCDXqFbAGr/+t3Ldvn3bv3q1LLrlEkrR+/Xqlp6ersLBQ7e3tknqm\nFzn5ykmXyzXk6YAA+Ke7u1sZGRlyOp2aM2eOpk2bppaWFjmdTkmS0+lUS0uLJGoWCDXqFbAGrxri\nzs5OXX/99Vq7dq3GjBmj5cuXq7a2Vnv27NH48eO1cuXKAZflIg0guMLCwrRnzx7V19frD3/4g954\n441e93tzEQeA4KBeAWsIP90Djh49qoULF+qmm27SggULJPXMg3fCsmXLlJOTI0mKj49XXV2d5776\n+nrFx8f3GTMpKcmn6UuAkSoxMVGfffaZqWOec845uuqqq/T+++/L6XSqublZcXFxampq8tSwNzVL\nvQK9Ua/A8OFtvQ56hNgwDBUWFio1NVUrVqzw3N7U1OT5/eWXX1ZaWpokKTc3V5s3b1ZXV5dqa2u1\nd+9ezZw5s8+4n3/+uWcu3WD83H///cQjnqVjmvUGduDAAc8pTN9++61ee+01ZWZmKjc3V2VlZZKk\nsrIyzz+33tRsMOo10M93MLbnSIjBOnj3Q72yn9glxkhYB2/rddAjxG+//bZ+85vfaMaMGcrMzJTU\nM8Xaiy++qD179sjhcGjy5Ml66qmnJEmpqanKy8tTamqqwsPDtWHDBj7OAYKoqalJBQUF6u7uVnd3\ntxYvXqy5c+cqMzNTeXl52rhxo2caJ4maBUKJegWsY9CG+Ic//KG6u7v73H5izuH+3HPPPXzHNhAi\naWlp+uCDD/rcHhMTox07dvS7DDULhAb1CliHLSYvzMrKIh7xLB/TzgL9fAdje46EGKwDvMF+Yp8Y\nI2EdvDXoF3MELKjDoRCEBSzLyjVh5dyAULByTVg5NyAUvK0JWxwhBgAAAAZCQwwAAABboyEGAACA\nrdEQAwAAwNZoiAEAAGBrNMTDRFRUjOc77X35iYqKCfUqAAAAWBLTrg0TPd9G5M9zxnNuZVauCSvn\nBoSClWvCyrkBocC0awAAACHi7ye7fNIbXIN+dTMAwDqiomLU0dFmyliRkdE6eLDVlLGAkcDM+vo7\n847Wd3Q4TBsLfXHKxDDBKRMjm5Vrwsq52Y3/rwO9RmO7+sjKNWHl3KzO3PqSJPPHY9sOHadMAAAA\nAF6gIQYAAICt0RADAADA1miIAQAAYGs0xAAAALA1GmIAAADYGg0xAAAAbI2GGAAAALbGN9UBQIAE\n5puvAABmoyEGgADpaYbN/uYrAIDZOGUCAAAAtkZDDAAAAFujIQYAAICt0RADAADA1miIAQAAYGs0\nxAAAALA1GmJgBKmrq9OcOXM0bdo0TZ8+XevWrZMkFRUVyeVyKTMzU5mZmdq2bZtnmZKSEiUnJysl\nJUXbt28PVeqA7VCvgHU4DMMwc5JM74I6HApB2GHN4XDIv/lMec6tzKyaaG5uVnNzszIyMtTZ2amL\nLrpI5eXl2rJliyIjI3X77bf3enx1dbUWLVqkXbt2qaGhQfPmzVNNTY3Cwv7+vzL16jv/67bPiCaO\nx3b1FfU6Mlm7XnvGY9sOnbc1wRFiYASJi4tTRkaGJGnMmDGaOnWqGhoaJKnfF4SKigrl5+crIiJC\nCQkJSkpKUlVVVVBzBuyKegWsg4YYGKH27dun3bt365JLLpEkrV+/Xunp6SosLFR7e7skqbGxUS6X\ny7OMy+XyvCEDCB7qFQgtGmJgBOrs7NT111+vtWvXasyYMVq+fLlqa2u1Z88ejR8/XitXrhxw2Z6P\nDQEEC/UKhF54qBMAYK6jR49q4cKFuummm7RgwQJJUmxsrOf+ZcuWKScnR5IUHx+vuro6z3319fWK\nj4/vM2ZRUZHn96ysLGVlZQUmecCC3G633G53QMamXgFz+VqvXFQ3THBR3chmVk0YhqGCggKde+65\neuyxxzy3NzU1afz48ZKkxx57TLt27dILL7zguUinqqrKc5HOZ5991uuoE/XqO2tfpMN29RX1OjJZ\nu157xmPbDp23NcERYmAEefvtt/Wb3/xGM2bMUGZmpiRp9erVevHFF7Vnzx45HA5NnjxZTz31lCQp\nNTVVeXl5Sk1NVXh4uDZs2MBHsECQUK+AdQx6hLiurk5LlizRl19+KYfDoZ///Oe69dZb1draqhtu\nuEFffPGFEhIStGXLFo0dO1ZSzxyJmzZt0qhRo7Ru3TrNnz+/b1D+gx0yjhCPbFauCSvnZnXWPuLE\ndvWVlWvCyrlZnbXrtWc8tu3QeVsTgzbEA82R+Mwzz2jcuHFatWqV1qxZo7a2NpWWlno1R+JQksPf\n0RCPbFauCSvnZnXWfoNlu/rKyjVh5dysztr12jMe23boTJmHeKA5Erdu3aqCggJJUkFBgcrLyyUx\nR6K1hcvhcPj1ExUVE+qVAAAAMJ3X066dmCNx1qxZamlpkdPplCQ5nU61tLRIYo5Eazumnv9Uff/p\n6GgLftoAAAAB5tVFdZ2dnVq4cKHWrl2ryMjIXvedOHo4kIHus9O0MFFRMTST6CWQ0zgBAIChOW1D\nfGKOxMWLF3vmSHQ6nWpublZcXJyampo8cyZ6O0ei1LshHul6mmF/z/vhSuKR5NR/AouLi0OXDAAA\nNjfoKROGYaiwsFCpqalasWKF5/bc3FyVlZVJksrKyjyNcm5urjZv3qyuri7V1tZq7969mjlzZgDT\nBwAAAPwz6CwTb731ln70ox9pxowZnlMfSkpKNHPmTOXl5Wn//v19pl1bvXq1Nm3apPDwcK1du1aX\nX35536A2uwrWnCtX/Z9lwowc7LTdgsnKNWHl3KzO2lets119ZeWasHJuVmfteu0Zj207dKZMuxYo\nditYGmKcjpVrwsq5WZ2132DZrr6yck1YOTers3a99ozHth06U6ZdAwAAAEY6GmIAAADYGg0xAAAA\nbI2GGAAAALZGQwwAAABboyEGAACArdEQAwAAwNZoiAEAAGBrNMQAAACwNRpiAAAA2BoNMQAAAGyN\nhhgAAAC2RkMMAAAAW6MhBgAAgK3REAMAAMDWaIgBAABgazTEAAAAsDUaYgAAANgaDTEAAABsjYYY\nAAAAtkZDDAAAAFujIQZGkLq6Os2ZM0fTpk3T9OnTtW7dOklSa2ursrOzNWXKFM2fP1/t7e2eZUpK\nSpScnKyUlBRt3749VKkDtkO9AtbhMAzDCHpQh0MhCBsyDodDkr/r6+8Y5uRgp+0WTGbVRHNzs5qb\nm5WRkaHOzk5ddNFFKi8v1zPPPKNx48Zp1apVWrNmjdra2lRaWqrq6motWrRIu3btUkNDg+bNm6ea\nmhqFhf39f2W71auZzKn9XiOaOB7b1VfU68hk7XrtGY9tO3Te1gRHiIERJC4uThkZGZKkMWPGaOrU\nqWpoaNDWrVtVUFAgSSooKFB5ebkkqaKiQvn5+YqIiFBCQoKSkpJUVVUVsvwBO6FeAeugIQZGqH37\n9mn37t2aNWuWWlpa5HQ6JUlOp1MtLS2SpMbGRrlcLs8yLpdLDQ0NIckXsDPqFQit8FAnAMB8nZ2d\nWrhwodauXavIyMhe9zkcjr99NNi//u4rKiry/J6VlaWsrCyzUgUsz+12y+12B2x86hUwj6/1SkMM\njDBHjx7VwoULtXjxYi1YsEBSz1Gm5uZmxcXFqampSbGxsZKk+Ph41dXVeZatr69XfHx8nzFPfoMF\n7ObUprK4uNi0salXwFy+1iunTAAjiGEYKiwsVGpqqlasWOG5PTc3V2VlZZKksrIyzxtvbm6uNm/e\nrK6uLtXW1mrv3r2aOXNmSHIH7IZ6BayDWSaCgFkmcDpm1cRbb72lH/3oR5oxY4bno9SSkhLNnDlT\neXl52r9/vxISErRlyxaNHTtWkrR69Wpt2rRJ4eHhWrt2rS6//PKA5GZH1r5qne3qK+p1ZLJ2vfaM\nx7YdOm9rgoY4CGiIcTpWrgkr52Z11n6DZbv6yso1YeXcrM7a9dozHtt26Jh2DQAAAPACDTEAAABs\njYYYAAAAtkZDDAAAAFujIQYAAICtnbYhvvnmm+V0OpWWlua5raioSC6XS5mZmcrMzNS2bds895WU\nlCg5OVkpKSnavn17YLIGAAAATHLahnjp0qWqrKzsdZvD4dDtt9+u3bt3a/fu3bryyislSdXV1Xrp\npZdUXV2tyspK3XLLLeru7g5M5gAAAIAJTtsQz549W9HR0X1u729Ot4qKCuXn5ysiIkIJCQlKSkpS\nVVWVOZkCAAAAAeDzOcTr169Xenq6CgsL1d7eLklqbGyUy+XyPMblcqmhocH/LAEAAIAA8akhXr58\nuWpra7Vnzx6NHz9eK1euHPCxJ76OEgAAALCicF8Wio2N9fy+bNky5eTkSJLi4+NVV1fnua++vl7x\n8fH9jlFUVOT5PSsrS1lZWb6kAgxLbrdbbrc71GkAAABJDsOLL3jet2+fcnJy9NFHH0mSmpqaNH78\neEnSY489pl27dumFF15QdXW1Fi1apKqqKjU0NGjevHn67LPP+hwlttt3rZvz/ej+jmFODnbabsFk\n5Zqwcm5WZ07t9xrRxPHYrr6yck1YOTers3a99ozHth06b2vitEeI8/PztXPnTh04cEATJ05UcXGx\n3G639uzZI4fDocmTJ+upp56SJKWmpiovL0+pqakKDw/Xhg0bOGUCAAAAlubVEWLTg9rsP1iOEON0\nrFwTVs7N6qx9xInt6isr14SVc7M6a9drz3hs26Hztib4pjoAAADYmk8X1QEAhrtwU09pi4yM1sGD\nraaNBwDBREMMALZ0TGZ+nNvRwfUiAIYvTpkAAACArdEQAwAAwNZoiAEAAGBrnEMMAH8TFRWjjo62\nUKcBAAgyGmIA+JueZtjseUgBAFbHKRMAAACwNY4QAwCAYYnTnGAWGmIMgf8T+TN5PwDALOae5sQp\nTnZGQ4wh8H8ifybvBwAAVsM5xAAAALA1GmIAAADYGg0xMILcfPPNcjqdSktL89xWVFQkl8ulzMxM\nZWZmatu2bZ77SkpKlJycrJSUFG3fvj0UKQO2Rb0C1uEwDMPMSTe9C+pwKARhfWbOVaz+rq/DzzH8\nXd68MYbTtg8Ws2rizTff1JgxY7RkyRJ99NFHkqTi4mJFRkbq9ttv7/XY6upqLVq0SLt27VJDQ4Pm\nzZunmpoahYX1/j95uNWrP3ouGjV7HmKrjmd+bnbaT6hXazC3Zq1crz3j2WnbmsXbmuAIsRf+fhWr\nrz9AcMyePVvR0dF9bu/vxaCiokL5+fmKiIhQQkKCkpKSVFVVFYw0AYh6BayEhhiwgfXr1ys9PV2F\nhYVqb2+XJDU2Nsrlcnke43K51NDQEKoUAfwN9QoEH9OuASPc8uXLdd9990mS7r33Xq1cuVIbN27s\n97EDzTNdVFTk+T0rK0tZWVlmpwlYltvtltvtDkos6hXwj6/1SkMMjHCxsbGe35ctW6acnBxJUnx8\nvOrq6jz31dfXKz4+vt8xTn6DBezm1KayuLg4YLGoV8A/vtYrp0wAI1xTU5Pn95dfftlzRXtubq42\nb96srq4u1dbWau/evZo5c2ao0gQg6hUIFY4QAyNIfn6+du7cqQMHDmjixIkqLi6W2+3Wnj175HA4\nNHnyZD311FOSpNTUVOXl5Sk1NVXh4eHasGGD31/NDcB71CtgHUy75gX/p3WxwpRnVsihZ4zhtO2D\nxco1YeXczMa0a/6NZ6f9xKrrauXcAoFp13A6TLsGAAAAeIFTJgAAACwv3LTTZCIjo3XwYKspY40U\nNMQAAACWd0xmnYLR0cH556filAkAAADYGg0xAAAAbI2GGAAAALZGQwwAAABboyEGAACArdEQAwAA\nwNZoiAEAAGBrNMQAAACwNRpiAAAA2BoNMQAAAGzttA3xzTffLKfTqbS0NM9tra2tys7O1pQpUzR/\n/ny1t7d77ispKVFycrJSUlK0ffv2wGQNAAAAmOS0DfHSpUtVWVnZ67bS0lJlZ2erpqZGc+fOVWlp\nqSSpurpaL730kqqrq1VZWalbbrlF3d3dgckcAAAAMMFpG+LZs2crOjq6121bt25VQUGBJKmgoEDl\n5eWSpIqKCuXn5ysiIkIJCQlKSkpSVVVVANIGAAAAzOHTOcQtLS1yOp2SJKfTqZaWFklSY2OjXC6X\n53Eul0sNDQ0mpAkAAAAEht8X1TkcDjkcjkHvBwAAAKwq3JeFnE6nmpubFRcXp6amJsXGxkqS4uPj\nVVdX53lcfX294uPj+x2jqKjI83tWVpaysrJ8SQUYltxut9xud6jTAAAAkhyGYRine9C+ffuUk5Oj\njz76SJK0atUqnXvuubrzzjtVWlqq9vZ2lZaWqrq6WosWLVJVVZUaGho0b948ffbZZ32OEjscDnkR\n1jJ68vcnX3+XN2MMK+TQM8Zw2vbBYuWasHJuZvO/1vuMaOHxzM/NTvuJVdfVyrkFgrk1a+V6NXs8\n++wn3tbEaY8Q5+fna+fOnTpw4IAmTpyoBx54QHfddZfy8vK0ceNGJSQkaMuWLZKk1NRU5eXlKTU1\nVeHh4dqwYQOnTAAAAMDSvDpCbHrQYfYfLEeIzR1jOG37YLFyTVg5N7NxhNi/8ey0n1h1Xa2cWyBw\nhNj3seyyn3hbE3xTHQAAAGyNhhgAAAC2RkMMAAAAW6MhBgAAgK3REAMAAMDWaIgBAABgazTEwAhy\n8803y+l0Ki0tzXNba2ursrOzNWXKFM2fP1/t7e2e+0pKSpScnKyUlBRt3749FCkDtkW9AtZBQwyM\nIEuXLlVlZWWv20pLS5Wdna2amhrNnTtXpaWlkqTq6mq99NJLqq6uVmVlpW655RZ1d3eHIm3AlqhX\nwDpoiIERZPbs2YqOju5129atW1VQUCBJKigoUHl5uSSpoqJC+fn5ioiIUEJCgpKSklRVVRX0nAG7\nol4B66AhBka4lpYWOZ1OSZLT6VRLS4skqbGxUS6Xy/M4l8ulhoaGkOQIoAf1CoRGeKgTABA8Dofj\nb191OvD9/SkqKvL8npWVpaysLJMzA6zL7XbL7XYHPS71Cgydr/VKQwyMcE6nU83NzYqLi1NTU5Ni\nY2MlSfHx8aqrq/M8rr6+XvHx8f2OcfIbLGA3pzaVxcXFAYtFvQL+8bVeOWUCGOFyc3NVVlYmSSor\nK9OCBQs8t2/evFldXV2qra3V3r17NXPmzFCmCtge9QqEBkeIgREkPz9fO3fu1IEDBzRx4kQ98MAD\nuuuuu5SXl6eNGzcqISFBW7ZskSSlpqYqLy9PqampCg8P14YNGwb9eBaAuahXwDochmEYQQ/qcCgE\nYX3W86LjT77+Lm/GGFbIoWeM4bTtg8XKNWHl3Mzmf633GdHC45mfm532E6uuq5VzCwRza9bK9Wr2\nePbZT7ytCU6ZAAAAgK3REAMAAMDWaIgBAABgazTEAAAAsDUaYgAAANgaDTEAAABsjYYYAAAAtkZD\nDAAAAFujIQYAAICt0RADAADA1miIAQAAYGs0xAAAALA1GmIAAADYGg0xAAAAbI2GGAAAALZGQwwA\nAABboyEGAACArdEQAwAAwNZoiAEAAGBrNMQAAACwNRpiAAAA2Fq4PwsnJCQoKipKo0aNUkREhKqq\nqtTa2qobbrhBX3zxhRISErRlyxaNHTvWrHwBAAAAU/l1hNjhcMjtdmv37t2qqqqSJJWWlio7O1s1\nNTWaO3euSktLTUkUAAAACAS/T5kwDKPX31u3blVBQYEkqaCgQOXl5f6GAAAAAALG7yPE8+bN08UX\nX6ynn35aktTS0iKn0ylJcjqdamlp8T9LAAAAIED8Oof47bff1vjx4/XVV18pOztbKSkpve53OBxy\nOBz9LlsbPJcrAAAPmUlEQVRUVOT5PSsrS1lZWf6kAgwrbrdbbrc71GkAAABJDuPUcx58VFxcrDFj\nxujpp5+W2+1WXFycmpqaNGfOHP3lL3/pHdTh6HOqhZX1NPX+5Ovv8maMYYUcesYYTts+WKxcE1bO\nzWz+13qfES08nvm52Wk/seq6Wjm3QDC3Zq1cr2aPZ5/9xNua8PmUicOHD6ujo0OSdOjQIW3fvl1p\naWnKzc1VWVmZJKmsrEwLFizwNQQAABhBoqJiPJ8em/EDX4Wbuh2iomJCvUJ+8/kIcW1tra699lpJ\n0rFjx/TTn/5Ud999t1pbW5WXl6f9+/cPOO3acPsPliPE5o4xnLZ9sASjJnydJnG41as/OELs33h2\n2k+oV99Yu8asnJvZ49mn/r2tCdNOmRgKqxfsqWiIzR1jOG37YAlGTUyePFnvv/++YmL+/p/8qlWr\nNG7cOK1atUpr1qxRW1tbn6kSh1u9+sPab9Zmj2d2bhGSjpk2WmRktA4ebDVtPDNRr76zdo1ZOTez\nx6MhPhXfVAfYCNMkInCOqecN1pyfjo62IOdvPdQrEDw0xIBNME0iMHxQr0Bw+TXtGoDhw59pEgEE\nF/UKBBcNMWAT48ePlySdd955uvbaa1VVVSWn06nm5mbPNImxsbH9Lsu84bCzUMwbTr0CvvG1Xrmo\nzgtcVGfuGMNp2wdLoGvi8OHDOn78uCIjI3Xo0CHNnz9f999/v3bs2KFzzz1Xd955p0pLS9Xe3j7s\nLtIxk7Uv+DF7PCvn1jOeVfc76tV31q4xK+dm9njUa5/HjfSGOCoqxqSLM4Z7M2qFHCR/r0S38pXn\n/gh0TdhpmkR/WPvN2uzxrJxbz3hW3e+oV99Zu8asnJvZ41GvfR430htic4pvJDSjVsjBjDGsW3T+\nsPKbmJVzM5u136zNHs/KufWMZ9X9zso1YeXcJKvXmJVzM3s86vVUzDIBAAAAW+OiOgDDmnmnRQEA\n7IqGGMCw1tMMm/kxIgDAbjhlAgAAALZGQwwAAABboyEGAACArdEQAwAAwNZoiAEAAGBrNMQAAACw\nNRpiAAAA2BoNMQAAAGyNhhgAAAC2RkMMAAAAW6MhBgAAgK3REAMAAMDWaIgBAABgazTEAAAAsDUa\nYgAAANgaDTEAAABsjYYYAAAAtkZDDAAAAFujIQYAAICt0RADAADA1miIAQAAYGvhoU5gMI8//pSe\nffa3oU4DAAAAI5ilG+Lt2/+g99//f5Lm+jhCi6TXTcwIoRcuh8Ph1wiRkdE6eLDVpHwAYOSKiopR\nR0dbqNMAAs7SDXGP6ZLm+bjsfjMTgSUck2T4NUJHh38NNQDYRU8z7N9rbm+8/sKaOIcYAAAAthaQ\nhriyslIpKSlKTk7WmjVrAhECgEmCXa9RUTFyOBym/WCkCjdtH4mKign1ypiG91dYk3n1GqqaNb0h\nPn78uH7xi1+osrJS1dXVevHFF/XJJ5+YHWaI3MQj3tAiuoMfMxRCUa9//wj25J83+rnN2x9vuE3L\nP7QxAs1toRgnTo8a6k/ffWmknANrnfdX9wiIEejxg8VtkfF9rVfr1KzpDXFVVZWSkpKUkJCgiIgI\n3XjjjaqoqDA7zBC5iUe8oUW0SUNsnXp1D/PxgxUj0NwjIEagxw8d+9RrMGIEevxgcQ/z8YMV4/RM\nv6iuoaFBEydO9Pztcrn0f//3f2aHAWACb+v10UcfNSUepzgAvgt2vQKh4/+MUkOPaDIzVyA8PExn\nnvmfiojY7NPyhnFYHR2mpYMRw7tCKy4u7vf2kTRtm7f1unLlygBnAgRS8N9cA4F6hX34P6PU33lX\nN6Y3xPHx8aqrq/P8XVdXJ5fL1esxiYmJQ3px+vbbPX5mVfy3H38M9cX01HhmvBgPNoY362dmDv48\nn77kcXK80L6xdXS0mf7mmpiYaOp43vK2Xj///HOTI/f3/AV6n/J2fH+2bX8xzN5fzRwvFNthKDF8\nXVd/X+8HR71KI3c/sXK9DjSer9vCKtthoBjmPHfe1qvDMAwzJxjUsWPHdMEFF+j111/XhAkTNHPm\nTL344ouaOnWqmWEAmIB6BYYP6hUIHNOPEIeHh+vxxx/X5ZdfruPHj6uwsJBiBSyKegWGD+oVCBzT\njxADAAAAw4lp0661trYqOztbU6ZM0fz589Xe3t7v4waaVPy3v/2tpk2bplGjRumDDz7w3L5v3z6d\neeaZyszMVGZmpm655ZaAxpOkkpISJScnKyUlRdu3bzdtHQda/tR1vPrqq0878fqtt96q5ORkpaen\na/fu3T7HHmwZf+IVFRXJ5XJ51qmystKUeDfffLOcTqfS0tJ6PT5Q6zdQvECsX11dnebMmaNp06Zp\n+vTpWrdunVfrF0yPPPKIwsLC1Npq7kWF9957r9LT05WRkaG5c+f2Ok/SLHfccYemTp2q9PR0XXfd\ndfrmm29MjzHY64o/Av1lDAPt52YabP82w3fffadZs2YpIyNDqampuvvuu00d/4Tjx48rMzNTOTk5\nARnfTIGqVynwNUu9Dox6HRqva9YwyR133GGsWbPGMAzDKC0tNe68884+jzl27JiRmJho1NbWGl1d\nXUZ6erpRXV1tGIZhfPLJJ8ann35qZGVlGe+//75nmdraWmP69OlBi/fxxx8b6enpRldXl1FbW2sk\nJiYax48fNyXmQMufvI6DLX/C//zP/xhXXnmlYRiG8e677xqzZs3yOXag4hUVFRmPPPLIkJ6f08Uz\nDMP4wx/+YHzwwQd99olArN9g8QKxfk1NTcbu3bsNwzCMjo4OY8qUKcYnn3wy6PoF0/79+43LL7/c\nSEhIML7++mtTxz548KDn93Xr1hmFhYWmjm8YhrF9+3ZPLd95550BeQ4Hel3xhzf7lL8G2s/N1N/+\nbfZ6HDp0yDAMwzh69Kgxa9Ys48033zR1fMMwjEceecRYtGiRkZOTY/rYZgpkvRpG4GuWeh0Y9To0\n3tasaUeIt27dqoKCAklSQUGBysvL+zxmsEnFU1JSNGXKlJDHq6ioUH5+viIiIpSQkKCkpCRVVVWZ\nEtPf5ftb91mzZqm9vV3Nzc0+xQ5UPEky+jkbx594kjR79mxFR0f3GTcQ6zdYPLPXr6WlRXFxccrI\nyJAkjRkzRlOnTlVDQ8Og6xdMt99+u371q18FZOzIyEjP752dnRo3bpzpMbKzsxUW1vOSN2vWLNXX\n15seY6ivY94IxpcxDLafm6W//buxsdHUGGeddZYkqaurS8ePH1dMjLlf/1pfX69XX31Vy5Yt67f+\nrSSQ9SoFvmap14FRr94bSs2a1hC3tLTI6XRKkpxOp1paWvo8pr9JxU+84Q+mtrZWmZmZysrK0ltv\nvRXQeI2Njb2msTl5GX9jDrb8iXVctmyZIiIiBs15oBiNjY1Dju3Nc+RLPElav3690tPTVVhY6PmI\n3594gwnE+p2Omet36ov9vn37tHv3bs2aNWvQ9QuWiooKuVwuzZgxI2Ax/v3f/12TJk1SWVmZ7rrr\nroDFkaRNmzbpJz/5SUBjmMXX/dPKTt2/zdLd3a2MjAw5nU7NmTNHqamppo5/22236eGHH/Y0alYV\njHqVglez1GtoDdd6lYZWs0OaZSI7O9tz5OxkDz30UK+/HQ5Hv/O0DnXu1uzsbDU1NWnSpEk6duyY\n6uvr9eMf/1jPP/98QOJJPS8k5eXlKi0tldSzc7/22mu9mtShxDQMY8DHnbh9woQJqqurU3R0tB5+\n+GEVFxero6Oj13/gp/Lm6IQ3sb19joZ6NGT58uW67777JPWca7Zy5Upt3LjR53hD2ZZmrN/plgvk\n+nV2dur666/X2rVrNWbMmD7LDrTv+Wuw+i4pKel1Pr0vR8cGGn/16tXKycnRQw89pIceekilpaW6\n7bbb9Mwzz5geQ+pZn9GjR2vRokVDHt/bGGYaCV8ocbLT7d/+CAsL0549e/TNN9/o8ssvl9vtVlZW\nliljv/LKK4qNjVVmZqYlvto90PU6WAyzapZ6tb7hWq/S0Gt2SA3xa6+9NuB9TqdTzc3NiouLU1NT\nk2JjY/s8xptJxU8Xb86cOUpOTg5IPEmei/ZO/Ld7xRVXqLi4WLNmzfIpZn19veLj4yUN/ByNHj1a\no0ePltTzUcjo0aO1d+9eXXjhhf3m3F8Ml8ulo0ePDjm2N8+Rt/FOXvbk52bZsmWeFx5f451Yj4GY\nvX6nixeo9Tt69KgWLlyom266SQsWLDjt+plpoPr+85//rNraWqWnp3vyveiii1RVVTWkPAZ7/TjZ\nokWLfD4adLoYzz77rF599VW9/vrrPo3vTQyz+fI6ZlUD7d9mO+ecc3TVVVfpvffeM+0N9o9//KO2\nbt2qV199Vd99950OHjyoJUuW6LnnnjNl/KEKdL0OFuNUvtYs9Wptw7leJR9q1qyTlu+44w6jtLTU\nMAzDKCkp6fcE+KNHjxrnn3++UVtbaxw5cqTfk82zsrKM9957z/P3V199ZRw7dswwDMP4/PPPjfj4\neKOtrS1g8U5cVHfkyBHjr3/9q3H++ecb3d3dpqzjQMufvI6ffvqpMWrUKOPDDz8cMOeTL8p65513\nPBdl+RLbm+fIl3iNjY2e5R999FEjPz/f73gn9HehZSDWb7B4gVi/7u5uY/HixcaKFSuMU3mz7wVL\nIC7Sqamp8fy+bt0646abbjJ1fMMwjG3bthmpqanGV199ZfrYpzr1dcUf3uxTZhjoAmazDLZ/m+Gr\nr74y2traDMMwjMOHDxuzZ882duzYEZBYbrfbuPrqqwMyttkCdVFdoGuWeh0c9To03tSsaQ3x119/\nbcydO9dITk42srOzPSva0NBg/OQnP/E87tVXXzWmTJliJCYmGqtXr/bc/l//9V+Gy+UyzjjjDMPp\ndBpXXHGFYRiG8bvf/c6YNm2akZGRYVx44YXGK6+8EtB4hmEYDz30kJGYmGhccMEFRmVlpWnrONDy\nv//973utY1FRUZ/ln3zySePJJ5/0jPXP//zPRmJiojFjxoxeV8cONfZAy/gbb/HixUZaWpoxY8YM\n45prrjGam5tNiXfjjTca48ePN0aPHm24XC5j06ZNAV2/geIFYv3efPNNw+FwGOnp6UZGRoaRkZFh\nbNu27bTrF2yTJ082/Q124cKFxvTp04309HTjuuuuM1paWkwd3zAMIykpyZg0aZLnuV2+fLnpMQZ7\nXfHHQHVmloH2czMNtn+b4U9/+pORmZlppKenG2lpacavfvUr08Y+ldvttvwsEycEol4NI/A1S70O\njHodOm9qli/mAAAAgK1Z+1JZAAAAIMBoiAEAAGBrNMQAAACwNRpiAAAA2BoNMQAAAGyNhhgAAAC2\nRkMMAAAAW6MhBgAAgK39fwIdbgvQN9baAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x112f50050>"
       ]
      }
     ],
     "prompt_number": 164
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}