aphearin/random_subvol_mockpop.ipynb

## random_subvol_mockpop.ipynb
{
 "metadata": {
  "name": "",
  "signature": "sha256:1a755ff5836465fe86b364387172d1ae969fcc86d63587e9149b096bfb74e2e2"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Demo of how to generate Halotools predictions using random subvolumes"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In principle, *all* of this can get convenience wrapper functions put into the Halotools code base. However, I'm not sure how much of that is appropriate. On the one hand, the science case for being able to generate subvolume-based predictions is very clear. On the other hand, it has a fairly specific use-case, and it's not all that much code to implement by-hand, so it may be more appropriate to just create a demo page on `halotools.readthedocs.org` and put links to it at various places in the docstrings. There is also room for in-between scenarios, where some but not all of the functions in this notebook get incorporated into the repo. \n",
      "\n",
      "The code below shows the rigamarole that is necessary to do the following: \n",
      "\n",
      "1. Generate a masking function a random subvolume of a periodic box\n",
      "2. Populate `multidark` using the `masking_function` feature\n",
      "3. Compute the two-point function of the points in the subvolume\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Define the masking function"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def subvolume_masking_function(t, xmin, xmax, ymin, ymax, zmin, zmax):\n",
      "    \n",
      "    mask = np.ones(len(t), dtype = bool)\n",
      "    \n",
      "    mask *= t['halo_x'] >= xmin\n",
      "    mask *= t['halo_x'] < xmax\n",
      "    mask *= t['halo_y'] >= ymin\n",
      "    mask *= t['halo_y'] < ymax\n",
      "    mask *= t['halo_z'] >= zmin\n",
      "    mask *= t['halo_z'] < zmax\n",
      "    \n",
      "    return mask"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def get_random_tuple(num_divs):\n",
      "    return np.random.random_integers(0, num_divs-1, 3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def get_bin_boundaries(idx, Lbox, num_divs, verbose=False):\n",
      "    binsize = Lbox/float(num_divs)\n",
      "    lower_bound = idx*binsize\n",
      "    upper_bound = lower_bound + binsize\n",
      "    if verbose is True: print(lower_bound, upper_bound)\n",
      "    return lower_bound, upper_bound"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from functools import partial"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def get_random_masking_function(num_divs, Lbox, verbose=False):\n",
      "    ix, iy, iz = get_random_tuple(num_divs)\n",
      "    xmin, xmax = get_bin_boundaries(ix, Lbox, num_divs, verbose=verbose)\n",
      "    ymin, ymax = get_bin_boundaries(iy, Lbox, num_divs, verbose=verbose)\n",
      "    zmin, zmax = get_bin_boundaries(iz, Lbox, num_divs, verbose=verbose)\n",
      "    \n",
      "    f = partial(subvolume_masking_function, \n",
      "                xmin = xmin, xmax = xmax, \n",
      "                ymin = ymin, ymax = ymax, \n",
      "                zmin = zmin, zmax = zmax)\n",
      "    return f, xmin, xmax, ymin, ymax, zmin, zmax"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Populate a randomly selected subvolume"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from halotools.empirical_models import PrebuiltHodModelFactory\n",
      "model = PrebuiltHodModelFactory('zheng07')\n",
      "\n",
      "from halotools.sim_manager import CachedHaloCatalog\n",
      "halocat = CachedHaloCatalog(simname = 'multidark', redshift = 0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f, xmin, xmax, ymin, ymax, zmin, zmax = get_random_masking_function(num_divs=5, Lbox = halocat.Lbox, verbose=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(600.0, 800.0)\n",
        "(400.0, 600.0)\n",
        "(600.0, 800.0)\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "model.populate_mock(halocat = halocat, masking_function = f, enforce_PBC = False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Sanity checks on our results \n",
      "\n",
      "Let's just inspect that everything works as advertised"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(model.mock.galaxy_table['halo_x'].min(), model.mock.galaxy_table['halo_x'].max())\n",
      "print(model.mock.galaxy_table['halo_y'].min(), model.mock.galaxy_table['halo_y'].max())\n",
      "print(model.mock.galaxy_table['halo_z'].min(), model.mock.galaxy_table['halo_z'].max())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(600.00793, 799.99335)\n",
        "(400.02014, 599.99426)\n",
        "(600.01825, 799.99786)\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(model.mock.galaxy_table['x'].min(), model.mock.galaxy_table['x'].max())\n",
      "print(model.mock.galaxy_table['y'].min(), model.mock.galaxy_table['y'].max())\n",
      "print(model.mock.galaxy_table['z'].min(), model.mock.galaxy_table['z'].max())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(599.99120311265528, 800.34901486195088)\n",
        "(399.84426441111742, 600.20956488958211)\n",
        "(598.90828195885126, 800.11960622230413)\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "centrals_mask = model.mock.galaxy_table['gal_type'] == 'centrals'\n",
      "cens = model.mock.galaxy_table[centrals_mask]\n",
      "\n",
      "print(cens['x'].min(), cens['x'].max())\n",
      "print(cens['y'].min(), cens['y'].max())\n",
      "print(cens['z'].min(), cens['z'].max())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(600.0079345703125, 799.99334716796875)\n",
        "(400.0201416015625, 599.9942626953125)\n",
        "(600.01824951171875, 799.99786376953125)\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Looks good. Every selected halo lies within the subvolume. There do exist galaxies that lie slightly outside the subvolume. All such galaxies are satellites, as expected. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Second verse, same as the first"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f, xmin, xmax, ymin, ymax, zmin, zmax = get_random_masking_function(num_divs=5, Lbox = halocat.Lbox, verbose=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(400.0, 600.0)\n",
        "(200.0, 400.0)\n",
        "(800.0, 1000.0)\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "model.populate_mock(halocat = halocat, masking_function = f, enforce_PBC = False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(model.mock.galaxy_table['halo_x'].min(), model.mock.galaxy_table['halo_x'].max())\n",
      "print(model.mock.galaxy_table['halo_y'].min(), model.mock.galaxy_table['halo_y'].max())\n",
      "print(model.mock.galaxy_table['halo_z'].min(), model.mock.galaxy_table['halo_z'].max())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(400.00012, 599.99329)\n",
        "(200.00516, 399.98618)\n",
        "(800.00836, 999.99298)\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(model.mock.galaxy_table['x'].min(), model.mock.galaxy_table['x'].max())\n",
      "print(model.mock.galaxy_table['y'].min(), model.mock.galaxy_table['y'].max())\n",
      "print(model.mock.galaxy_table['z'].min(), model.mock.galaxy_table['z'].max())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(399.71707722257315, 600.25456444204667)\n",
        "(199.85509237951328, 399.98617553710938)\n",
        "(799.51701910371025, 1001.2407316564478)\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Computing observables"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from halotools.mock_observables import tpcf\n",
      "from halotools.empirical_models.factories.mock_helpers import three_dim_pos_bundle"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rbins = np.logspace(-1, 1.25, 15)\n",
      "rmax = rbins.max()\n",
      "rbin_centers = (rbins[1:] + rbins[0:-1])/2."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pos = three_dim_pos_bundle(model.mock.galaxy_table, 'x', 'y', 'z')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Since we have PBCs turned off, we must use `randoms`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "num_randoms = 1e5\n",
      "xran = np.random.uniform(xmin, xmax, num_randoms)\n",
      "yran = np.random.uniform(ymin, ymax, num_randoms)\n",
      "zran = np.random.uniform(zmin, zmax, num_randoms)\n",
      "randoms = np.vstack((xran, yran, zran)).T"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "clustering = tpcf(pos, rbins, \n",
      "                period=None, \n",
      "                randoms = randoms, \n",
      "                approx_cell1_size = [rmax, rmax, rmax])\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's compare our results to the full simulation volume case with PBCs turned on (a highly non-trivial sanity check)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "r2, xi2 = model.compute_average_galaxy_clustering(simname = 'multidark', \n",
      "                    approx_cell1_size = [40, 40, 40], num_iterations=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(r2, r2*r2*xi2, label='computed from subvolume')\n",
      "\n",
      "plt.plot(rbin_centers, rbin_centers*rbin_centers*clustering, \n",
      "         '--', color='red', label='computed from entire box')\n",
      "\n",
      "plt.loglog()\n",
      "plt.ylim(ymin = 25, ymax = 150)\n",
      "plt.xlim(xmin = 0.1, xmax = 30)\n",
      "plt.legend(frameon=False, loc='best')\n",
      "plt.xticks(size=15)\n",
      "plt.yticks(size=15)\n",
      "plt.xlabel('R', fontsize=15)\n",
      "plt.ylabel(r'$R^{2}\\xi(R)$', fontsize=25)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 36,
       "text": [
        "<matplotlib.text.Text at 0x12f187910>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEZCAYAAAA+MZraAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VFX6x/HPoWnQgKH3Liuo6C4uRWniTyzrKoq6WBAE\n7Looq2sDiWV1ZVfBtiorC6yIKCJYFhakBAELKmJBUJReXFBqgJD2/P44SUwgZSaZyc0k3/frNS8y\nM3fuPbkZ5plz7jnP48wMERGRIFQKugEiIlJxKQiJiEhgFIRERCQwCkIiIhIYBSEREQlMlaAbEGuc\nc5pOKCISJjNz+T2unlAxmFmZvY0aNarM7784+wjnNaFsW9Q2hT1f0HPRPvdl4W+n90b5fG9Eu42F\nURAqZ3r16lXm91+cfYTzmlC2LWqbwp6P9jmOFr039N4oTFBtd0VFKcnLOWc6Z5KfxMREEhMTg26G\nlEEV/b3hnMM0HCcSXbH8LViiS++NgqknFCb1hEREwqOekIiIlEkKQiIiEhgFIRERCYyCkIiIBEZB\nSEREAqMgJCIigVEQEpEyq1evXowfP77A56+99lpq1apFly5dSrFVJTNo0CBGjhwZM/uNNiUwFZGo\nWL9+Pa1atSI9PZ1KlYr3fdc5h3P5Li9h8eLFzJs3j61bt3L00UeXpKmlqrDfqSzuN9rUExKRqIrW\n4u4NGzbQokWLAgNQenp6VI4bCdE6J7G4kF5BSCRGbdq0iUsuuYR69epRp04dbrvtNgAyMzN55JFH\naNGiBfXr12fgwIHs3bsX8L2TSpUqMXHiRJo1a0bt2rV54YUX+OSTT+jQoQMJCQk5+wGYOHEiZ5xx\nBrfddhvHHXcc7dq1Y8GCBTnPt2jRgvnz5+fcT0xMZMCAAQD06NEDgOOOO474+Hg+/vhjAP71r3/R\nvn17atWqxbnnnsvGjRtzXv/ee+9xwgkncNxxx3HbbbcVmIV5/PjxXHfddXz44YfEx8fz4IMPkpSU\nRJMmTRg9ejQNGzZkyJAhpKamcvvtt9O4cWMaN27MHXfcQWpqKkDO9n/729+oV68ejRo1YubMmcya\nNYu2bdtSu3Zt/vrXvxZ4/mfNmsWJJ55IjRo1aNKkCU888UTOOevevXuebStVqsTatWtz7v/000/0\n6dOHGjVq0KtXr5xzcNNNN3HXXXflee1FF13EmDFjAFi1ahW9evUiISGBk046iXfeeSfPttk9oaLa\nMGjQIG6++WbOP/984uPj6d69Oz/++CPDhg0jISGBdu3asWLFipzXbt26lX79+lGvXj1atWrFM888\nU+B5CVvQ6cNj7eZPmUiw0tPTrUOHDjZ8+HA7cOCApaSk2NKlS83MbPz48damTRtbt26dJScn2yWX\nXGIDBgwwM7N169aZc85uuukmO3TokM2dO9eqVatmffv2tR07dtiWLVusXr16tmjRIjMzmzBhglWp\nUsXGjh1r6enp9tprr1nNmjVt165dZmbWokULmz9/fk67EhMT7eqrrzYzs/Xr15tzzjIyMnKenzlz\nprVp08ZWr15tGRkZ9sgjj9jpp59uZmY7duyw+Ph4mz59uqWnp9uYMWOsSpUqNn78+HzPwcSJE61b\nt2459xcuXGhVqlSxe+65x1JTU+3gwYM2cuRI69q1q+3YscN27Nhhp59+uo0cOTLP9g8//LClp6fb\nP//5T6tdu7ZdeeWVlpycbCtXrrS4uDhbv359vsdv0KCBLVmyxMzMdu/ebcuXL885Z7nbZWbmnLMf\nfvjBzMwGDhxo8fHxtnjxYjt06JANGzYsZ/v333/fmjZtmvO6nTt3WlxcnG3bts1SU1OtdevW9thj\nj1laWpotWLDA4uPj7dtvvzUzs0GDBuX8bqG0oU6dOrZ8+XJLSUmx3r17W/Pmze3ll1+2zMxMGzFi\nhJ155plmZpaRkWG/+c1v7OGHH7a0tDRbu3attWrVyubMmZPveclP1udm/p+pBT2hm4KQFA0icwvX\nBx98YHXr1s3zAZ+td+/e9vzzz+fc//bbb61q1aqWkZGRE4S2bt2a83zt2rXt9ddfz7nfr18/Gzt2\nrJn5D7NGjRrl2X+nTp1s8uTJZnZkEBo1alROEMo+Vu42nnvuuXmCSkZGhlWvXt02bNhgkyZNsq5d\nu+Y5VpMmTQoMQod/0C5cuNCqVatmhw4dynmsdevWNnv27Jz7c+bMsRYtWuRsHxcXZ5mZmWZmtnfv\nXnPO2bJly3K279ixo82cOTPf4zdr1sxefPFF27NnT6HtMjsyAFxxxRU5zyUnJ1vlypVt8+bNlpmZ\nac2aNbP333/fzMzGjRtnZ511lpn5ANWgQYM8+73iiissMTHRzMILQoMGDbLrr78+57lnnnnG2rdv\nn3P/yy+/tOOOO87MzD766CNr1qxZnn09+uijdu211+Z7XvJTWBDScJxICUQqDIVr06ZNNG/ePN8L\n/tu2baN58+Y595s1a0Z6ejr/+9//ch6rX79+zs9xcXFH3N+/f3/O/caNG+fZf/Pmzdm6dWv4jcZf\nx8ke8klISKB27doAbNmyhW3bttGkSZM82zdt2jSs/detW5dq1arl3N+6desR5yJ322vXrp0zhBUX\nFwcceW5yn4vcpk+fzqxZs2jRogW9evXio48+CqmNzrk8v+cxxxxDrVq12Lp1K845+vfvz6uvvgrA\nlClTuOqqq3J+l8PPR0n+FvXq1cv5+eijj85zPy4ujuTkZMD/zbZu3ZrzN0tISOCxxx5j+/btxTru\n4RSERGJQ06ZN2bhxIxkZGUc816hRI9avX59zf+PGjVSpUiXPh2s4tmzZkuf+hg0baNSoEeA/QHN/\nSP/44485P+c3U6tZs2aMGzeOXbt25dz2799P165dadiwIZs2bcrZ1szy3A/F4cfM71xkt72kTjvt\nNGbOnMmOHTvo27cvl19+OeDPyYEDB3K2y31O4MjfKzk5mZ07d+a064orruCNN95gw4YNLFu2jH79\n+uX8Lps2bcoekQH83+LwLwmhtCEcTZs2pWXLlnn+Znv37uXdd98t9j5zUxASiUGdO3emYcOG3HPP\nPRw4cICUlBQ++OADwH+IjRkzhvXr15OcnMx9991H//79w5omnfuDbvv27Tz99NOkpaUxbdo0Vq9e\nzfnnnw/AqaeeytSpU0lPT+fTTz9l+vTpOYGgbt26VKpUiR9++CFnXzfeeCOPPvoo33zzDQB79uxh\n2rRpAJx//vmsXLmSGTNmkJ6eztNPP12iD8/sc/HII4/w008/8dNPP/HQQw/lTJwoibS0NF555RX2\n7NlD5cqViY+Pp3LlygCccsoprFy5ki+++IKUlJR8i9nNmjWLpUuXkpqaysiRI+natWtOMDn11FOp\nU6cOQ4cO5dxzz6VGjRoAdOnSherVqzN69GjS0tJISkri3XffpX///sAvl1ZCaYOF0f3u1KkT8fHx\njB49moMHD5KRkcHXX3/Np59+Gu5py5eCkEgMqlSpEu+88w7ff/89zZo1o2nTprz++usADB48mAED\nBtCjRw9atWpF9erV88xmCmUtSe5tOnfuzJo1a6hbty4jR45k+vTpJCQkAPDwww/zww8/kJCQQGJi\nYs7QEUD16tW5//77OeOMM0hISGDZsmX07duXu+++m/79+1OzZk1OPvlk5syZA0CdOnWYNm0a99xz\nD3Xq1OH777+nW7duhbbx8N/l8PsjRozgtNNOo0OHDnTo0IHTTjuNESNGFLh9OOtsJk+eTMuWLalZ\nsybjxo3jlVdeAaBt27Y88MAD/N///R+/+tWv6N69e579Oue46qqrePDBB6lduzaff/45kydPzrPv\nK6+8kgULFnDllVfmPFa1alXeeecdZs+eTd26dbn11lt5+eWXadu27RHnI5Q2FHY/97moXLky7777\nLitWrKBVq1bUrVuX66+/PmfGZUmpqF2YVNROKpKJEycyfvx4Fi9eHHRTJIapqJ2IiJRJCkIiUqBY\nTQUjsUPDcWHScJyISHg0HCciImWSgpCIiARGQUhERAKjICQiIoFREBIRkcAoCIlImVUey3sX5KST\nTuL999+PyL5y13Uq61TeW0SiQuW9CzZo0CCaNm3Kww8/nPPY119/HbH9x9LaLvWERCSqorWuLpbL\ne4crv2zphYmltYwKQiIxSuW9gy3vfejQIe68806aN29OgwYNuOmmm0hJScmz7yeffJL69evTqFEj\nJk6cCMC4ceOYMmUKo0ePJj4+nosuuijnXGaf28TERC699FIGDBhAzZo1mTRpEnv27GHIkCE0atSI\nJk2aMHLkSDIzM/Ntm3OOlJQU+vfvT40aNejYsSNffvllzvMFlQn/4IMPqFu3Lps3bwbgiy++oFat\nWnz33XcFnocSK6janW6qrCpll8p7B1/e+/bbb7eLLrrIdu3aZfv27bPf//73du+99+bZ96hRoyw9\nPd1mzZpl1atXt927d5tZ3iqo2XKfy1GjRlnVqlXtrbfeMjOzgwcPWt++fe3GG2+0AwcO2Pbt261T\np0724osv5tu27Ndnn8u///3v1rJlS0tPTy+yTPj9999vvXv3tgMHDthJJ51kzz33XL7HCAcq760g\nJFEyapTlWyx11KjQty9o20KovHew5b0zMzPtmGOOySmXbeb/Ji1btsyz79y/e7169ezjjz82Mx+E\nRowYkWefhwehnj175jz3448/2lFHHWUHDx7MeWzKlCl25pln5ntuRo0aledcZmZmWsOGDW3x4sVF\nlglPS0uzjh072kknnWTnnXdevvsPV2FBSBMTREoiMdHforV9AWK9vPef/vSnPI/HWnnvHTt2cODA\nATp27JjzmJnlGR6rXbt2nr9P9erVc0pmhyL3udiwYQNpaWk0bNgw57HMzEyaNWsW0uuzS4pn/+75\nlQnPrqBbpUoVBg4cyLBhwxg7dmzI7S0uXRMSiUEq752/0irvXadOHeLi4vjmm29yfo/du3eHXOgt\n3MKCTZs25aijjuLnn3/OOd6ePXv46quvCnx97nOXmZnJ5s2bady4cYFlwrOD1pYtW3jooYcYPHgw\nw4cPz7mGFi0KQiIxSOW9QxOt8t6VKlXiuuuu4/bbb2fHjh2A//CeO3duSK+vX78+a9euDfl4DRs2\npE+fPgwfPpx9+/aRmZnJDz/8UOi6os8++yznXI4dO5ajjz6aLl260KlTpwLLhJsZgwYNYujQobz0\n0ks0bNiQkSNHhtzO4lAQEolBKu8dfHnvxx9/nDZt2tClSxdq1qzJ2WefnWcWWWH7GjJkCN988w0J\nCQlccsklIf1u//73v0lNTc2ZWXjZZZcVGKSdc/Tt25fXXnuNWrVq8corr/Dmm29SuXJlqlWrVmCZ\n8KeffpqffvopZ/3ShAkTmDBhAkuXLg35vIRL9YTCpHpCUpGovLdEguoJFcE519Q5N985941z7mvn\n3ONBt0lEpCJQEPLSgLvMrD3wa6Czc+7IPrJIBaPy3hJtGo7Lh3PuaeB7M3s6n+c0HCciEobChuO0\nTugwzrnaQF/g7KDbIiJS3sX0cJxzro1z7kXn3JfOuQzn3MICtmufdc1nv3Nui3PuQefcEb+7c+4o\n4A1gjJl9G+32i4hUdLHeE2oPnAd8iP9djhgnc84lAPOAr4ELgTbAE/gAPDLXdpWBV4DPzGxM1Fsu\nIiKxfU3I5bpA45x7A6hlZr0P2+Ze4E6guZklZz12F5AINDCzfVmPvQRUMrPBRRxT14RERMJQbqdo\nhxgNzgPmZAegLK8BcUAPAOfcGcBgoKNz7vOs260Rb7CIiOQR68NxofgVfjguh5ltdM4dAE4A/mNm\nS4nxgCwiEosqQhBKAHbn8/iurOfClpgrC3KvXr3o1atXcXYjIlIuJSUlkZSUFNK2MX1NKLdCrgml\nAncevubHObcJmGRmIwiDrgmJiISn3F4TCtEuoGY+jydkPSciIgGpCEFoNdAu9wPOuaZA9aznREQk\nIBUhCM0GznHOHZvrsT8AB4BFwTRJREQgxoOQcy7OOXepc+5SoDFQL/u+cy4ua7MXgEPAm865s5xz\n1wOjgCcPm7YdO95+G37+OehWiIiUWExPTHDOtQCyyxNm/yIu6+eWZrYxa7t2wLNAV/x1oJeAxOLM\nMCgTExOuugo6dYJhw4Jth4hICAqbmBDTQSgIZSIILVzoA9AXX4DS7ItIGVfRZ8eVDwsXwv79/uee\nPf3Pn34abJtEREpIQSgWbNsGl14Ku7JmlFeqBIMHw0svBdsuEZES0nBcmAIZjrv1VjjqKHjiiV8e\n27IFTj4ZNm2CY44p3faIiIRBRe1i2bp18OqrsPqwJU2NG8OMGVC1ajDtEhGJAPWEwlTqPaFBg6B5\nc3jwwdI7pohIBKknFGGJiYmlk7h0/XqYNQvWrInucUREoiCURKbqCYWp1HtCW7b4oTcRkRildUIR\nVCbWCYmIxBCtEyrv9u2DDRuCboWISNgUhMqDadPgj38s0S5SUyEjI0LtEREJkYbjwlQmh+OSk6Fp\nU/jmG2jYMKSXmMHKlTB3rr8tWeKDUOvW0LYtHH+8/zf7Vq+eMgSJSPHomlAERT0IzZ7tZ8OF27MZ\nOtRHjrvvLnCTn36CefNgzhwfeKpVg3POgT59oHdvf//77+G77468paXlDUrZt+OPhxo1Svg7i0i5\npiAUQVENQpmZ8Otfw0MPwUUXhffaDz/0a4pWr87psqSm+ofnzvWBZ80a6NXLB50+faBNm9B7Nz//\n7F9/eHBas8YHofwCVKtWPtGDiETJuHG+tMs775TpoQqtE4oVr78ORx8NF14Y/mu7dMEqV2bzq0t4\na2d35s6FRYt8MDjnHHjySejSxfd2iqN2bX/r0iXv45mZsHVr3sC0aJH/d+NGP0r4+9/DH/7gq0+U\n4f8nIrEnIwMefjim/2OpJxSmqPWE0tKgfXt44QU466yQX7Z7N8yf73s7aTP/wyaa0uyCDvTp43dT\np07kmxqqtFRj1WrHm2/6zEOpqdC/PwxpMZ/Ws5/BJSRArVqQ/W+HDtCtW3ANFpGoKNXhOOdcZaAJ\nUBM4FkjBF5LbZGbpET1YAKIWhP75T3jtNX/RphBpab6CQ/Z1na+/hjPO+OXaTrt2ZehL0XXX+W7Q\nhRdi5ssfTZ0KSa9s4bTMZZzdcSddfrWL+tV2wc6dfijy+uuP3M8rr8DIkXkDVvPmcM01cNJJpf97\niZS2lSvhxBODbkWxRTUIOeeqAmcBffGVS08A8suqmQmsApYCM4H3zCzmJgVHLQjdfjtccQV07pzn\n4T17/HWdpUv97ZNP/Ay2Pn184DnjDD+CV+Z89pkPQKtXHzFzwQyWLfMB6fXXfW+tf38/ZNeqVT77\nOnDAl7PYlRWsdu70/ylbtvQlLUTKq88/h/vug2+/9T/XrBl0i4olKkHIOXcMcAvwJ6AG8CmwHFgH\nbAGSgUPAUfgeUWOgNfBr4DTgZ+AZ4GkzO1CsRgQgmhMTzPya0+yAs3Qp/PADnHaaDzbdukHXrnDc\ncVE5fOSY+cJ711zjZ+0VIjPTTw+fOhXeeANatPAB6fLLoUmTEhy/zHQHRYphzRrf+1+0CEaM8KMK\nxb2gWwZEPAg5584CJgFfAeOAuWa2P4zXHwucC1wP/Aq40cxmh92QAEQyCKWn+yGqpUv9B/HSpf46\nY7duPuiccYYfoYq5ag1vvAGPPOJ7Q5Urh/yy9HRfQHbqVJg5018i69/f1/OrXz/EnZjBb34DHTv6\nINitmy8CKBIr3nrLf3m74w4YNqxc1AuLaBByzt0PnA/cbmafRKBxnYExwCwze6Sk+4u2kgShvXvz\nDq0tW+YvbWQHnDPO8MNREfkSf/AgxMVFYEdhSknx0WP8eDjzzGLvJjXVX/OaOhXefRd++1s/XHfJ\nJf6SUKG2bvXXkSZN8kN5Awb4gNS6dbHbI1Jq9u7138iKfKPHjogFIefc3fgJByMjeT0nazLDQ8DB\nsh6InHM2atSoIks5mPkpyrmH1r7/3n9Bzw44XbtG6X22YweccooviFfaC3V27YKXXy5xGqHcDh70\nFS2mTvWBqXt3uPpq30OqUtgiAzM/jj5pkj8Xb78dsTaJSNGySzk8+OCDEQtCl5vZ6xFr4ZH772dm\n06O1/0gItSe0fbuvvp0dcLp180NreYZ1d+3ys72i4ayz4IYb/MWVcmTfPr8u78UXYfNmnyBi4EAt\nipUYk5rqF5q2b+/TlZRgN888AxdfXMCknjIiYlm0oxmAsvZfpgNQOOrVgx9/hDffhD/9yU96yxOA\ntm6FE07wXe9oGDLED4mVM/HxcOWV/nrtxIn+/LZpA08/7UfewvbQQ/4P9OWXkW6qyJEyMvxIwQkn\n+O59vXrF3tWcOX5p3fz5sX3Zs0RNd86d6pwb7py72jkXH8L2/Zxzzzvn6pbkuLGi0Gs7jzzir1NE\nK/HaxRf7iQHluMRD9+7w3//CjBl+QkOrVvD44763FLL+/f0c9wsugFNP9ZMqRKLhs8/8e+yFF/w3\nqFmzirXObe1a6NsXbrkF/v53v5sWLSLe2tJjZsW64admZ+a6bQO6hPC63wDvATWLe+wgb/6UldD3\n35vVrm22Y0fJ91WYW24xS0yM7jHKkC+/NLviCrM6dfyv/fPPYbw4I8Ns7lyzFi3M7r7bLDMzau2U\nCqpHD7NJk4r93kpONhsxwn90PPqoWUpKhNsXRVmfm/l+pharJ+ScOwV4HL8OaEfWw/WBWc65pkUE\nveXAZOCK4hy7XEhM9Bfuo51TZ+jQ0lkvs2QJfPVV9I9ThJNPhilT/CSQDRt8hu977vHX54pUqRKc\nfbZfDXziiVpnJJE3Z44f/QjzvWXmF3W3a+fXDa5YAffeW36ugxZ3ndCLQDxwrZkdyhpeuwp4AEgy\ns0uytmuN7/n8DKwHfgRSgeFAqpmNjcQvUZpKvE7o66/9pIHvv/cXOGJdaqr/0H72WZ/CoQzZsAFG\nj/Z56wYMgLvuKsECWJEAfPWV/766a5efgNC9e9AtKp5olPfuBNxqZocAzGxHVkDpDHR3zh2fNZ17\nNfAafvhtDbAP33t6DPi8mMeObXXq+K/r5SEAgQ8+2am6y5jmzeG553yGn6pV/UXcG27wY+oiZdmu\nXT74nHUWXHaZzxcZqwGoKMWemGBmO/N5bA3wCDAaH2g2AF/irxm5XLdPzGxRcY8d0xo0CCtLdpn2\n00/w2GP+6mgZ1rChb+J330Hdur6kxDXXwKpVYezkhx9U/1yiLiPD5zJu184nK161Cm6+uYj1cDGu\nuMNxi80s37jsnGsAbAUeMrPErMfigR7AicBuYIqZJRe30UEqk+W9g3LLLT4tz9NPB92SsOze7XtI\nTz/tU9zdd5+ftFSoP/zBT7ubMiUGkvdJ4Mx81o7LLgv54s2HH8Jtt/nJms8849cVlhfRyB230MwK\nzMninNsKNLNyULrhcApCWfbu9ZlVP/zQV7uLQcnJfr3g3//uM1ncf/+RRftypKX5i0qzZv2S2E6k\nIBMmwBNP+Nxc1asXuum2bX4Czfz5fonBlVeWv3kx0bgmVFRKzTXlMQDFtIED/cyvSKlRA775JmYD\nEMCxx8Lw4f4a0Xnn+c7Otdf6tF1HqFoVxo71kapnT784SSQ/q1bBn//s64MVEoBSU/0XoJNP9kPG\nq1bBVVeVvwBUlOIGoY7OuUedc72cc/nlF88sSaPKnSVLfA6zILVtG/kMCuVkoProo/24+6pV8L//\n+Zx0KSkFbDxwoO8NDR/urxOJ5HbwoE+V9dhjhRahy852sHChH0z461/Lz1ylcBV3OC53kEkBPgYW\nAklZP/+3iOG6U81sRdgHLgPCHo7LyPAXHB591Bd5C8qWLf4r1+bNRQ4PVGSpqT456s6dftTt2GML\n2DAlpYxWE5RA3Xijr0Q5ZUq+XZq1a/33l6+/9h3rCy4IoI0BiMZw3H7gVWAjcDTQE0jEB6FdwK+d\nc6Odc+c75/LLS/NUMY8be6ZO9V9xgn63NW7s03YrLU2hqlXz64patPBrV3ftKmBDBSA5XFqaHx14\n8cUjApCZv0TUqZPPI7lyZfAfCWVFcXtCH5lZl6yfm+ODUK+sf1setnkG8AWwCB+kPgAWm1lMXtkN\nqyeUluYTFY4fD4WUfSg1b74JTz3ls38WR0qK/5SO5WyJITKDO++EefP80EmDBkG3SGJVcrKvQr9+\nPUyb5tevVTTR6Al9m/2DmW0ws3+b2WAzaw00A64BxgPfA5XxWRPuAN4CtuOrqcasxMREkpKSit7w\nX//yhdTKQgAC/9UrOdkPFxTHvff6FAQVgHP+onG/ftCjR4h5YBcs8Fm5M3VJVLzvvvM9nxo14P33\nK14ASkpKIjExsdBtitUTCodzrhG+h5TdW2qLT2YXet3nMiTknpCZLyw3frwvCxrrvvvOF0b65hu/\n4rMCeeopePJJX1DvV4V9ffrxRx+16taFf/87ehnSJSa89RZcdx385S/+34os4uuEStiYlsByM4tS\nNbfoCms47sCB8jMJ4MILfd6Qu+4KuiWBmDDBz86eNauIha2pqTBsGCQl+U+htm1Lq4lSRmRk+BzF\nkyb54bfOnYNuUfCiMRxXbGa2Dp9HrvwrLwFo3jx/JTWCJbtjzbXX+gwLffrABx8UsmG1avD8834K\nVLduMHt2qbVRStnWrf7LWa6FZTt3+lHvJUt8vjcFoKIFdYV5Y0DHlXClp8Mdd8Df/lZ+cscX06WX\n+qKYffvCe+8VsfF11/meUHn5IiJ5ZWT4laW//W3OerkVK/zdE0/0748SFE2tUEp9OC7WVbi0PRkZ\nPjtAv34Vbyl3AZYsgUsu8TNxL7446NZIIB56yM8ynTsXKldm8mT/Xe3ZZ33mDcmrTF0TinXlJgjN\nm+eTVg0YEHRLYtLy5fC73/mV7gMHBt0aKVWLFvmy8MuXk1q7IXfe6UddZ8woVrXuCiFi14Scc1FN\nHxzt/UsuxxwDjzziZ/FJ2H7zG59yZeRI/+1XKojdu31KjQkT2EZDeveGdet8WkYFoOIJ95rQS865\nqJxq51xHYGI09i356NLFl2FYsiTolsSsE07waz/GjvXTcEOK52+/XUhiOinzataEKVNYGn8uv/2t\nn6jy1lv2Bth7AAAVX0lEQVSq7lES4Qah64F/OOduiFQDnHOVnHN/BJ4Fro3UfqUIzsGQIZFPalrB\ntGgBixf77Ex//nMRgcgMJk/2RWMkJhmO577sziWX+DIgDzxQIRKIRFXY14Scc8fgsyG0BR4G3jaz\nsEtOOueqAJcB9wErgSFmtj/c/ZS2cnNNCGD7dr+OZePGvAsrly71CU8vvzy4tsWYnTt9OYhTTvEz\ntCsXtBR73z4/b/dPf/JfAiRmHDzo85OuWOEzYLVuHXSLYkdE1wmZ2X4z649PWPoosNU5N945N8A5\n16GAhKU45+Kdc6c45wY7514GfgRGAPeZWf9YCEDlTr16vtT4O+/88lhGhv+mrtQzYalVy8/1+P57\nP3M3La2ADePj/SfYPffAZ5+Vahul+Natg9NP9ysWPvhAASiSSjQ7zjlXCfg9MAA4F8heFGH4Mt7J\nQDXgWOCYrOf2A/8BXjaz/xT74AEpVz0h8BVS4+N/mX49YYIvcr90qaZkF0NKiu9AZmb61fJxcQVs\n+MYbPvvEp5/GdGHAcs+Mha/voP8f63HffX69tv5bhK9Upmg756oCvwY6AC2ABHyZhxRgJ7AO+BJY\nEctVV8tdEMotOdknR5sxw+ecl2JJS4NBg/yI5ttvF5JC7okn/PL6QhPSSVAyM2F2v5eoPevfpL73\nPj16BN2i2KV1QhFUroPQiBE+3/zkyUG3JOZlZsItt/iOzn//q85OrNm7F0ZcvJIH3+9F6nvvU79X\nu6CbFNOikjvOOXdsAaW9JRZlZsKyZb4ssZRYpUrwj3/4S249esCOHUG3SEK1bh30OO0A93x+Occ+\nN1oBKMqK3RPKmlyw38xujGyTyrZy3ROSqLjrLl8B4913dT2hrFuzxn9x+G+TobRvleKTBeqPVmLR\nyqJdBz9NO78DdnbOFZbwPnu7Qc65Mc65M0vQDpEy7dFH4aeffBbuImWEvdpBImTVKjjzTHjy+tW0\n3/exn2uvABR1JekJzQHizOyIy3XOuclADTO7sJDXXwW8DKQC6cAdZvbPYjWmFKknJMWxdq1fHvTe\ne4XUI5o40dchmjBBH36l7Kuv4Jxz4PHHs9IppqfnZMeWkotoT8g593/OuceBk4CCiqVU5pcp2QW5\nDcjEl/5uDgwNty1BCbm8t0iWVq18ep/+/WF/QSviLrvMrx0aN65U21bRff45nH02jBmTK5+vAlBE\nRKW8t3MuDjgVOAu4BOhpZvsO2+YtoJ6ZdS1gH+2Br4F5ZtYn67H38+tVlTXqCUlJXHONL8v0z4L6\n/GvW+DLq77yjimjR9vnnbJi8mE6T/8jzz/vyHBIdkc6YcNDMPgTGAE2Aj3Nf08lawAqwtpDdZOcr\n+VfWa1oVsb1IufDcc37E7Y03Ctjg+ON9hLrsMk2pi5YdO+CGG0j9v/MYO64648crAAWp2BMTstLs\nXIGfnDDPObfJObcU2ARcALTJWsCaR1bOuKvxGRVmZD18LjC2uG0RiRXx8TBlCtx8s0/Zl6+LLvK5\nf+6/v1TbVu6lpcFTT0H79mz+OY4T3SrOfWMoF1wQdMMqthIvVnXO/R54EWiQ9VAmcHfWz6cD/c0s\nLdf2twDPAP8ws1tLdPAAaDhOImH0aD/itnBhAZcf0tN9xsz4+FJvW7k1YgQsW8aHfxjLhfe057XX\noHfvoBtVMUQ9Y4Jz7iigB1AL+NjM1mc9PgboDfwZWAp0BmYCVYB22dvFEgUhiYTMTD8bq1s3GDUq\n6NZUEIcOMXt+NQYOckyfDt27B92giiPQtD3OuWuA0UC9rIf2A4PMbHpUDxwlCkISKdu2+Qqtr7+u\nD8TS8PbbMHSoL0LXNd8pUxItgeeOy7o21B2f1HSJmf0v6geNEgUhiaT//MdfH1qxAhISgm5NOWDm\nL7q1a+cjfJbp030uv3ffhdNOC7B9FVS0MiaEzMzSzGyBmU2P5QAkEmm/+x307QvXX19EVdaDB/2s\nOX0BKthnn/nxzTFj8jz86qtw660+kawCUNmjwrQiAXv8cb88qMhK6//4Bzz7bKm0KaZs3+7H2S64\nAAYP9ol4s3pBkybBnXcWkalCAqVSDmHScJxEw6pVPtv2++/7kaR8rV3rL2bMmOHLfIqfRXjiiT4A\nPfAA1KyZ89RLL8GDD/oAdMIJAbZRgr8mVJ4oCEm0jBvnOzsffQRHH13ARv/5D9xwgx96ql+/VNtX\nZu3fD8fkzRL23HN+Gvz8+dCmTUDtkhwKQhGkICTRYuYTJTRp4vPMFeiBB3yXad485TjLx5gx8Mwz\nsGABtGgRdGsEIjAxwTnXxTn3lnNupXPubefcb3M919o597Bz7hnn3J2RarRIReOcn3swYwbMmlXI\nhqNGwfnnQ2pqqbUtUPv3w9/+BuedV+TEjL/+1fcmFy1SAIoVRfaEnHNnAAuAn4GNQEv8otR7zezv\nuba7H3jIzCpHr7nBU09Iom3xYrj8cli+HBo2DLo1AUpO9uNqTz4JvXrByJFw0kn5bmoGDz/sZ8LN\nnw+NGpVuU6VwJe0JPQQMN7NGZtbFzOoDFwKDnHO359ouA1ARFJES6t7dT9keONBnVqiQxo3z9S9W\nrPDjaq+9VmgAGjECpk3zyWEVgGJLKEHooJk9l/sBM5sNnAZ0yCpOJyIRNHIkHDjgOwEhK0899DZt\n/Jjaq6/62W8FMPPl02fN8nn4NFcj9oQShNbk96CZpZjZYKC5c+4yoBz9DxAJVpUq8MorfobXp5+G\n8AIzuPBCn5OmPOjdu5C56p4ZDBvmY9X8+VCnTim1TSIqlCD0P+dcc+dconPur4c/aWaPAnHA2SgQ\niURM8+b+ksgVV8C+fUVs7JzvPt1yC/zlL7HRK9q1y084KMYECzO47Ta/LnXePKhVKwrtk1IRShAa\nAwzGl2e4PL8NzOzfwJPAgcg1TUQuuwx69vQfuEXq1Ml/Kr/1lo9cB8rof8edO33APP54v0q3wHrn\n+cvM9LF2+XKYOzfP+lSJQUUGITM7ZGajgPpAniuDzrlTnXPDnXNXA4uBFoXtyznXzzn3vHOubgna\nLFKhPPWUX8D66qshbNyokR+fqlLFz3DYtSvq7QvZzz/7Qn3HH+9TiC9bBv/6V1iZW7MD0Bdf+Fxw\nNWpEsb1SKkJe6WZme3Pfd879Cfhbrof+B1yMn8pd0D6mO+fWAVOcc5ea2Z4w2ytS4RxzjA9A55wD\nnTv7SWOFiouDl1+GmTPhuONKpY0h+eADX1r7s8+KtYgnM9NnHP/6ax+AVO+vfChWxgTn3CnAZ0Aa\nsBfI7tnsBk4xs01FvH4gEGdmL4R98IBpnZAEZexYmDrVryOqWjXo1pSuzEy48UY/ejdrlgJQrIlG\nKYebgdeB47LWDdUHhmc991SuA7d2zl3mnOvtnGvlnKvunKuStX1B2bHKvMTERJKSkoJuhlQww4b5\nC/CJiUG3pAhffAF7IjfIkZnp0+WtXq0AFGuSkpJILOINW9ye0OfAWWa287DHjwc+AE4HLgEeASqT\n/6y53ma2KOyDB0w9IQnS9u3w61/D5Mlw5pnF2MHmzXDssdEZplu82OfN+fxzn3uoc+cS7zIz0y/c\nXbPG52499tgItFNKXVSK2h0egLIeW4MPPKOBx4ANwJdAJj6bQvbtk1gMQCJBq1cPJkyAa67x1/nD\n9vrr0KULfPddZBpk5suVdusG114LF13kS05EKAANHaoAVN4Vtye02My6F/BcA2ArPo9cYtZj8UAP\n4ET8daMpZpZc3EYHST0hKQvuvNPHkRkzoHK42Rr/+U8/S23yZOjTp2QN+fRTuO46uPde6NevGI3J\nX0aGD0Dr1/sYd1ilBokxES/l4JxbaGYFDgY457YCzcwsPeydl3EKQlIWpKb6RNr16/vqoWFXdMjO\nknr33f5ikytB2kezkr3+MBkZMGQIbNwI77yjAFQeRGM4rqi5OWvKYwASKSuqVfMf0Dt3Qv/+xUg6\n0L27X3w0YQLMmVP09nv2FDz+F+EANHgwbNqkHlBFUdwg1NE596hzrpdzrlo+z1fU3L8ipSYuzi8F\nSk/3I2EpKWHuoHlzH4jOOafgbX78Ee65xy9Omj27RO0tSkaGv6y0ZYsPsNWrR/VwUkYUNwgdBdyD\nrzO02zm30Dn3gHOuh3PuqKJe7Jw7tZjHFZFcjjrKlzCoXt3nLw07U09cXP49mR9+gJtugvbtfVqd\n5cvh6qsj0ub8ZGT40hXbtsHbbysAVSTFvSa0D3gbOANodtjTKUAqMA5IApbkk21hkZn1LE6Dg6Zr\nQlIWpaf7YayNG/0wVolmkm3fDief7Ccc/PGPfkpeFKWn+wC0Y4dPexcXF9XDSQCiMTHhIzPrkvVz\nc6An0Cvr35aHbZ4BfAEswgelD4DFZtY+7AOXAQpCUlZlL+pcudKPnJUoseehQ76bFWXp6b9MN585\nUwGovIpGEJpkZgMLeK4JvwSknkCbwzYxgFgtA64gJGVZZqaf7PbRR36+QVkucZCeDgMG+ByrM2fC\n0TGbQ0WKEvEgFObBG/FLQOoFtAVMQUgkOszgz3/2ZQ7mzYO6ZTBnfXq6v8S0Z49f66QAVL4FGoTy\naUxLYLmZhZ6/vQxREJJYYAajRsEbb/iqow0bBt2iX6SlwVVXQXIyvPmmAlBFUFgQCneJW4mZ2Trn\nXL4lw0UkMpyDhx7yH/A9e/pA1LRp0K3yAejKK/0svhkzSuWyk5RxpR6EsmwM6LgiFcp99/lA1KMH\nLFgALQ+fNlSK0tJ8wdeUFN8DUgASCCgImdmlQRxXpCIaPvyXHtG8edC2bem3IS3NZ3ZIS4Pp0xWA\n5BfFzqItIrHj5pt9HaIzz4Rvvim94yYnw/jx0LWrn4wwbZoCkOQV1HCciJSywYN9ADjrLF8e+5RT\nonMcMz9FfPx43+vp2dNPkjj//Igl2ZZyREFIpAK56iofiPr08ZkVfvvbyO17xw54+WV46SXf6xky\nxJfjbtAgcseQ8kdBSKSCufRSH4h+9zs/Q+2MM4q/r4wMeO893+t57z1f0+6FF3yS7ggm15ZyrNTX\nCcU6rROS8mLOHL9gdNo06NUrvNeuX++rQEyY4GsaDRniZ76VKFWQlFtlarFqrFMQkvJk4UL4wx9C\nK7J66JBPrzN+vE+qfeWVPvhE69qSlB8KQhGkICTlzdKlcPHFPrj8/vdHPv/VV/65V16BDh182e2L\nL1amAwmdglAEKQhJefTJJ3DBBfDcc/6a0d69MHWqDz5btvhic9de62vbiYRLQSiCFISkvFqxAs47\nz6/pWbAAevf2vZ5zztHUaikZBaEIUhCS8mz1ap9n7tJL/YQDkUhQEIogBSERkfAUFoSUtkdERAKj\nIFQMiYmJJCUlBd0MEZEyLSkpicTExEK30XBcmDQcJyISHg3HiYhImaQgJCIigVEQEhGRwCgIiYhI\nYBSEREQkMApCIiISGAUhEREJjIKQiIgERkFIREQCoyAkIiKBURASEZHAKAiJiEhgFIRERCQwCkIi\nIhIYBSEREQmMgpCIiARGQUhERAKjICQiIoFREBIRkcAoCImISGAUhEREJDAKQiIiEhgFIRERCYyC\nkIiIBEZBSEREAqMgJCIigVEQEhGRwCgIiYhIYBSEREQkMApCIiISGAUhEREJjIKQiIgERkFIREQC\noyAkIiKBURAqhsTERJKSkoJuhohImZaUlERiYmKh2zgzK53WlBPOOdM5ExEJnXMOM3P5PaeekIiI\nBEZBSEREAqMgJCIigVEQEhGRwCgIiYhIYBSEREQkMApCIiISGAUhEREJjIKQiIgERkFIREQCoyAk\nIiKBURASEZHAKAiJiEhgFIRERCQwCkIiIhIYBSEREQmMgpCIiARGQUgkQlTyXQqi90bBFIREIkQf\nNFIQvTcKpiBUzkT7zR6J/RdnH+G8JpRti9qmsOdj9QNF7w29NwoTVNsVhMoZfdDog6Ygem/ovVGY\noNruzCyQA8cq55xOmIhImMzM5fe4gpCIiARGw3EiIhIYBSEREQmMgpCIiARGQSjCnHPPO+c2O+cy\ng26LlB3OuZOcc8udc985595yzh0bdJukbKjonxkKQpH3CvCboBshZc4LwH1m1hZYDfw54PZI2VGh\nPzM0Oy5KnHOZZqYgLzjn6gOfmVmTrPttgRlmdmKwLZOypKJ+ZlS4X1gkAE2AzbnubwKaBtQWkTKl\nwgYh51wb59yLzrkvnXMZzrmFBWzX3jk33zm33zm3xTn3oHOuwp63iiSC75F8F+lJ7NLnR+RUCboB\nAWoPnAd8iD8PR4xLOucSgHnA18CFQBvgCXzwHpm1zRDg1qyX3GxmH0a95VJaIvIewfeCmuR6WTPy\n9owk9kTqvVHhVdhrQs45Z1m/vHPuDaCWmfU+bJt7gTuB5maWnPXYXUAi0MDM9hW0byCjIo7vlieR\nfI8455YAfzGz2c650cAhM9MHUYyK9OdHRf7MqHC/cDYLLfqeB8zJfgNleQ2IA3rm9wLn3EvARsCc\nc5ucc+NK3FgJRITfIzcBf3HOfQecAIyOWEOl1EXyvVHRPzMq8nBcKH6F707nMLONzrkDWc+9e/gL\nzGxoKbVNyoaQ3iNm9hUVeBpuBRXqe6NCf2ZU2J5QiBKA3fk8vivrORG9R6Qgem+EQEFIREQCoyBU\nuF1AzXweT8h6TkTvESmI3hshUBAq3GqgXe4HnHNNgepZz4noPSIF0XsjBApChZsNnHNYssk/AAeA\nRcE0ScoYvUekIHpvhKDCzo5zzsUBv8u62xiId85dmnX/P2Z2EJ908o/Am865x4HWwCjgycOmXUo5\npPeIFETvjcipyItVWwBrs+5mnwSX9XNLM9uYtV074FmgK34c9yUgMcR1AhLD9B6Rgui9ETkVNgiJ\niEjwdE1IREQCoyAkIiKBURASEZHAKAiJiEhgFIRERCQwCkIiIhIYBSEREQmMgpCIiARGQUgkRjnn\nEp1zmbluW51zM5xzbYNum0ioFIREYtseoEvW7Q6gPTDPOXdMoK0SCVGFTWAqUk6km9myrJ+XOefW\nAx8C5wFvBNYqkRCpJyRSvnyZ9W+LIBshEioFIZHypVnWv6rcKTFBw3EiMc45VxlfRqA5vmzAXuDt\nQBslEiIFIZHYVhtIy3V/L3Ceme0IqD0iYdFwnEhs2wOcBnQGbsD3iIYG2iKRMKgnJBLb0s1sedbP\nnzjnDgL/ds5NMbP5QTZMJBTqCYmUI2Y2GVgJjAq6LSKhUBASKX8eBbo557oF3RCRoigIicQuy7od\n7jVgDXBv6TZHJHzOLL/3sIiISPSpJyQiIoFREBIRkcAoCImISGAUhEREJDAKQiIiEhgFIRERCYyC\nkIiIBEZBSEREAvP/hlWckLOT/rYAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x12f19b7d0>"
       ]
      }
     ],
     "prompt_number": 36
    }
   ],
   "metadata": {}
  }
 ]
}