{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Signal and noise in DWI data\n",
      "\n",
      "In this notebook we will examine the distributions of SNR and noise in DWI data at different b values. \n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.stats as stats\n",
      "\n",
      "import nibabel as ni\n",
      "\n",
      "import osmosis.viz.mpl as viz\n",
      "import osmosis.utils as ozu\n",
      "import osmosis.model.sparse_deconvolution as ssd\n",
      "import osmosis.model.dti as dti\n",
      "import osmosis.model.analysis as oza"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import osmosis as oz\n",
      "import osmosis.io as oio\n",
      "oio.data_path = os.path.join(oz.__path__[0], 'data')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "subject = 'SUB1'\n",
      "data_1k_1, data_1k_2 = oio.get_dwi_data(1000, subject)\n",
      "data_2k_1, data_2k_2 = oio.get_dwi_data(2000, subject)\n",
      "data_4k_1, data_4k_2 = oio.get_dwi_data(4000, subject)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mask = np.zeros(ni.load(data_1k_1[0]).shape[:3])\n",
      "wm_nifti = ni.load(oio.data_path + '/%s/%s_wm_mask.nii.gz'%(subject, subject)).get_data()\n",
      "mask[np.where(wm_nifti==1)] = 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "TM_1k_1 = dti.TensorModel(*data_1k_1, mask=mask, params_file='temp')\n",
      "TM_1k_2 = dti.TensorModel(*data_1k_2, mask=mask, params_file='temp')\n",
      "TM_2k_1 = dti.TensorModel(*data_2k_1, mask=mask, params_file='temp')\n",
      "TM_2k_2 = dti.TensorModel(*data_2k_2, mask=mask, params_file='temp')\n",
      "TM_4k_1 = dti.TensorModel(*data_4k_1, mask=mask, params_file='temp')\n",
      "TM_4k_2 = dti.TensorModel(*data_4k_2, mask=mask, params_file='temp')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.bvecs\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.bvecs\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.bvecs\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.bvecs\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.bvecs\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.bvals\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.bvecs\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rmse1k = oza.rmse(TM_1k_1, TM_1k_2)\n",
      "rmse2k = oza.rmse(TM_2k_1, TM_2k_2)\n",
      "rmse4k = oza.rmse(TM_4k_1, TM_4k_2)\n",
      "\n",
      "rmse1k = rmse1k[np.isfinite(rmse1k)]\n",
      "rmse2k = rmse2k[np.isfinite(rmse2k)]\n",
      "rmse4k = rmse4k[np.isfinite(rmse4k)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_1.nii.gz\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b1000_2.nii.gz"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_1.nii.gz"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b2000_2.nii.gz"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_1.nii.gz"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Loading from file: /Users/arokem/usr/lib/python2.7/site-packages/osmosis/data/SUB1/SUB1_b4000_2.nii.gz"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = viz.probability_hist(rmse1k, label='b=1000', color=[0.8, 0.8, 0.8])\n",
      "fig = viz.probability_hist(rmse2k, fig=fig, label='b=2000', color=[0.59, 0.59, 0.59])\n",
      "fig = viz.probability_hist(rmse4k, fig=fig, label='b=4000', color=[0.32, 0.32, 0.32])\n",
      "ax = fig.axes[0]\n",
      "ax.set_xlabel('RMSE')\n",
      "ax.set_ylabel('P(RMSE)')\n",
      "plt.legend()\n",
      "ax.text(80, 0.05, 'median(1000): %2.2f \\nmedian(2000): %2.2f \\nmedian(4000): %2.2f'%(np.median(rmse1k), np.median(rmse2k), np.median(rmse4k)))\n",
      "fig.savefig('figures/rmse_distributions.svg')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlclXX6P/7X2dgOOyIgRwEF2QUUt5m+haWR1pCpmUum\n5jSkY06l5TRNpZ/GbTRLs0ZrKs1KKWvSKSK3UStDLHAZMUUCZREUkOVwDmfj/fuD333PuTnn4AHO\nhl7Px8OH3CsX27nOe7mvt4gxxkAIIYRYQezsAAghhPQdlDQIIYRYjZIGIYQQq1HSIIQQYjVKGoQQ\nQqxGSYMQQojV7J408vLyEBcXh5iYGKxbt87k+C+//IKxY8fCw8MDr732Gr+/oqIC48aNQ2JiIpKS\nkrB582Z7h0oIIeQmRPZ8TsNgMCA2NhYHDx5EeHg4Ro4ciV27diE+Pp4/5/r167h8+TK+/PJLBAQE\nYOnSpQCAmpoa1NTUIDU1FUqlEiNGjMCXX34puJYQQohj2bWlUVBQgOjoaERGRkImk2HGjBnYu3ev\n4Jzg4GCkp6dDJpMJ9oeGhiI1NRUA4O3tjfj4eFRXV9szXEIIITdh16RRVVWFgQMH8tsKhQJVVVXd\nvk95eTmKioowevRoW4ZHCCGkm6T2vLlIJOr1PZRKJaZNm4ZNmzbB29vb5vcnhJDbUU9HJuza0ggP\nD0dFRQW/XVFRAYVCYfX1Op0OU6dOxaOPPorJkyebPYcx5nL/XnnlFafHQDFRTLdjXBSTdf96w65J\nIz09HSUlJSgvL4dWq0VOTg6ysrLMntv5C2GMYcGCBUhISMDTTz9tzzAJIYRYya7dU1KpFFu2bEFm\nZiYMBgMWLFiA+Ph4bNu2DQCQnZ2NmpoajBw5Es3NzRCLxdi0aROKi4tx6tQpfPTRRxg2bBjS0tIA\nAGvWrMF9991nz5AJIYR0wa5JAwAmTpyIiRMnCvZlZ2fzH4eGhgq6sDh33HEH2tvb7R2eXWRkZDg7\nBBMUk3UoJuu5YlwUk/3Z9TkNexOJRL3unyOEkNtNb1477d7SIIQQcwIDA3Hjxg1nh3FLCwgIQEND\ng03vSS0NQohT0N+v/Vn6Hvfme08FCwkhhFiNkgYhhBCrUdIghBBiNUoahBBiJDIyEocOHXJ2GC6L\nkgYhhBgRiUS9qmun0+kwbdo0REVFQSwW4+jRoybnLF++HP369UO/fv3w5z//WXCsvLwc48aNg1wu\nR3x8vEkC++STTxAREQFvb2889NBDDp+BRkmDEEJs7M4778RHH32E0NBQkwS0bds27N27F2fOnMGZ\nM2fw73//m6+SAQAzZ87EiBEj0NDQgFWrVmHatGmoq6sDAJw7dw5PPvkkPv74Y9TW1sLLywuLFi1y\n6NcG1of18fAJua256t9vZGQkW7NmDUtISGABAQFs/vz5rK2trUf3UigU7OjRo4J9Y8eOZe+++y6/\n/f7777MxY8Ywxhi7cOECc3d3Z0qlkj9+5513sq1btzLGGHvhhRfY7Nmz+WOlpaXMzc1NcL4xS9/j\n3nzvqaVBCCFGGGP45JNPsH//fpSWluLixYv429/+hoqKCvj7+yMgIMDsv927d1t1/+LiYqSkpPDb\nw4YNw7lz5wB0tCQGDx4MuVzOH09JSREcN7528ODBcHd3x8WLF23xpVuFnggnhLikn3/+2Sb3GTFi\nRLfOF4lEWLx4McLDwwEAL774Ip566im8+uqraGxs7HU8SqUSfn5+/Lavry+USqXZYwDg4+ODq1ev\nAgBaW1tNjvv6+qKlpaXXcVmLkgYhxCV198XeloxXHB00aJBNl5r29vZGc3Mzv93U1MQvMNf5GHfc\nx8eHP97U1GTxuCNQ9xQhhHRy5coVwcfcgnLe3t7w8fEx+2/Xrl1W3TsxMRGnTp3it0+fPo2kpCT+\n2K+//sq3PLjjiYmJ/PHTp0/zx0pLS6HVajF06NBefb3d0uPREBfQx8PvEa1Wy86fP+/sMAjpNVf9\n+42IiGDJycmssrKS1dfXs9/+9rfsxRdf7NY92tramFqtZgqFgu3fv5+p1Wr+2NatW1l8fDyrqqpi\nlZWVLCEhgW3bto0/PmbMGLZs2TKmVqvZ559/zvz9/VldXR1jjLFz584xX19f9t133zGlUslmzpzJ\nZs6caTEOS9/j3nzvXfOnZiVX/aWzp8uXL7Nt27YxvV7v7FAI6RVX/fuNjIxka9euZQkJCczf35/N\nmzdP8KJvjYiICCYSiZhYLOb/v3z5Mn/8+eefZ4GBgSwwMJAtX75ccG15eTnLyMhgnp6eLC4ujh06\ndEhw/JNPPmGDBg1icrmcTZ48md24ccNiHPZIGlTlto/56aefcPLkSUybNg3BwcHODoeQHrsd/34d\njarcEly7dg2enp6ora11diiEkNsQJY0+hDGGa9euIT4+npIGIcQpKGn0Ia2trQCAIUOG4Nq1a06O\nhhByO6Kk0Ydcu3YN/fv3R2BgIFpbW6HRaJwdEiHkNkNJow/hkoZYLEZwcDC1NgghDkdJow/hkgYA\n9O/fn8Y1CCEOR0mjD6mrq0O/fv0AAIGBgQ6vo08IIZQ0+gi9Xg+dTgdPT08AgJeXF9RqtZOjIoTc\nbihp9BEqlQpyuZxf0MXT0xNtbW1OjoqQWw8t99o1Shp9RGtrK7y8vPhtDw8PamkQYge9Xe41Pz8f\nEyZMQFBQEPr374/p06ejpqZGcA4t90rsjmtpcLiWBpVhIMS1NDY24sknn8Tly5dx+fJl+Pj4YP78\n+fxxWu61C9988w2LjY1l0dHRbO3atSbHz58/z8aMGcPc3d3Zhg0bunUtY65b8MweTp8+zY4dOybY\n989//rPbhdQIcRWu+vdry+VeGWPs559/Zj4+Pvw2LfdqgcFgwOLFi5GXl4fi4mLs2rUL58+fF5wT\nFBSEN998E8uWLev2tbcblUol6J4COlob1EVFiG0xGy/3euzYMX69DICWe7WooKAA0dHRiIyMBADM\nmDEDe/fuRXx8PH9OcHAwgoOD8fXXX3f72ttNa2srv/wkh0saAQEBToqKEPv5xz/+YZP7LFy4sFvn\n23K51zNnzuDVV1/Fvn37+H203KsFVVVVgiUTFQoFTpw4Yfdrb1XU0iC3m+6+2NuSLZZ7vXTpEiZN\nmoTNmzfjt7/9Lb+/ry/3arek0ZvZB925dsWKFfzHGRkZyMjI6PHndWWdB8IBmnZLiL1YWu41Pj7e\n4uvTO++8g5kzZwIALl++jAkTJuDll1/G7NmzBedxy72mp6cDsLzcK5dITp8+jTlz5vDHe7Lc65Ej\nR3DkyJFufAcss1vS4L7JnIqKCigUCptfa5w0bmWdp9wCNO2WEHtgjOGtt97CAw88AE9PT6xatQqP\nPPIIBg4cKFi725KqqircfffdWLx4Mf7whz+YHH/sscewceNGTJo0CYwxbNy4EX/6058AAEOHDkVq\naipWrlyJV199Fbm5ufjvf/+LqVOnAgBmz56NsWPH4vvvv0daWhpeeuklTJ061eQNZWed31CvXLmy\nG9+RTno8hH4TOp2ODR48mJWVlTGNRsNSUlJYcXGx2XNfeeUVwewpa6+1Y/guRafTsa1bt7L29nbB\nfnMzqgjpK1z177e3y72uWLGCiUQi5u3tzf8znj3FGC33atE333yDp59+GgaDAQsWLMALL7zAz0fO\nzs5GTU0NRo4ciebmZojFYvj4+KC4uBje3t5mr+3sdlkusrm5GXv37uWbqJySkhKUlZXh3nvvdVJk\nhPTc7fL360z2WO6V1gjvA2pqavDDDz/wTVROZWUlCgsLkZWV5aTICOm52+Xv15lojfDbVGtrq9k+\nS5o9RQhxNEoafYC56bYAJQ1CiONR0ugDzE23BTpmT2k0GrS3tzshKkLI7YiSRh9gbrotAIjFYshk\nMlornBDiMJQ0+gC1Ws0vvtQZdVERQhyJkkYfoNFo4OHhYfYYJQ1CiCNR0nBBR44cwWeffcZvazQa\nuLm5mT2XSokQQhyJkoYLKisrw6FDh/h51BqNBu7u7mbPpZYGIbZFy712jZKGC6qvr8f169fxzjvv\n4NChQ2hra7OYNNzd3amlQYgN9Xa5V2P/93//B7FYjMOHDwv203KvxKYaGhowbtw4SKVS5OTkoKWl\nBVKp+dqS1D1FiGsqLS3Fnj17MGDAAMH+vr7cKyUNF1RfX4/p06djwYIFiI6O7nJKrYeHByUNQmys\noKAAiYmJCAwMxOOPP96jae2LFy/GunXrIJPJBPt37NiBZcuWYcCAARgwYACWLVuG7du3AwAuXryI\noqIirFy5Eu7u7pgyZQqGDRuGzz//HADw8ccfIysrC3fccQfkcjleffVVfPHFF2htbe3112wtShou\nxmAwoLGxkV+Nr1+/flAqlYJ3IsYoaRBiW8wGy71+9tln8PDwwMSJE03uT8u9Eptqbm6GXC7n350E\nBgbiypUruHz5MrKzs03Op6RBblWdC3T2FPcu3Vq9Xe61paUFL774Ig4ePGj2OC33Smyqvr4egYGB\n/La/vz/fNNbpdCZNXUoa5FbV3Rd7W+rNcq8rVqzAnDlzMGjQIH6fcUXZvr7cK3VPuZj6+noEBQXx\n28bvKsz1W1LSIMT2LC336u3tDR8fH7P/du3aBQA4fPgwNm/ejLCwMISFhaGiogLTp0/H+vXrAfxv\nuVeOpeVejY8nJibyx3uy3KstUUvDxVy9ehVhYWH8tkQiQVpaGqqqqqBUKuHv7y84XyaTwWAwQK/X\nW5xhRQixHuvlcq+HDh2CXq/n7zVy5Ei8/vrr/PiGM5Z7tSV6lXEx1dXVGDJkCL+t1Wpx77334osv\nvoBKpTI5XyQS8a0NrolLCOk5kUiE2bNn495770V1dTUmT56Mv/71r1Zfb9y9DHS88QsICOCLjmZn\nZ+PXX39FcnIyAOCJJ54QrCW+e/duzJs3D4GBgYiIiMDnn3/O9z4kJCRg69atmD17Nurr6zFhwgR8\n8MEHvf2Su4VW7nMxf/3rXzFjxgy+ufrDDz/A29sb+/btw/3334/hw4ebXJOTk4Px48cLurUIcXW3\n4t+vq6GV+24D1dXVgoeBuLpT3t7eFudie3h4UCkRQohDUNJwIa2trdBoNPwzGsD/6k55eXl1mTRo\nMJwQ4giUNFxIbW0tQkNDBXVvuKRxs5YGJQ1CiCPQQLgLUSqVJoPZXNKQy+Um87c5lDRIXxQQEGCz\nwoDEPONeC1uhpOFCVCqVybKuxkmDeyq0Mw8PD4c+EUqILTQ0NDg7BNID1D3lQswlDa1WyycNS3PE\nqaVBCHEUShouRKVSCR7S4R7ak8lkkMvlFsc0qDw6IcRRKGm4kM4tDa1WCzc3N4hEIvj5+VkslkYL\nMRFCHIWShgvpnDSMl3kNCwtDbW0tDAaDyXXU0iCEOAolDRfS2tpqMWm4u7vDz88P169fN7mOxjQI\nIY5CScOFdNXSAIABAwagqqrK5DrjooWEEGJPdk0aeXl5iIuLQ0xMDNatW2f2nCVLliAmJgYpKSko\nKiri969ZswaJiYlITk7GrFmzerTcYl9zs6ShUChQVVWF9vZ2Qd0Y46KFhBBiT3ZLGgaDAYsXL0Ze\nXh6Ki4uxa9cunD9/XnBObm4uLl26hJKSErzzzjtYuHAhAKC8vBzvvvsuCgsLcfbsWRgMBsFSircq\ncwPhxkkjPDwcV65cQXFxMWpqagTXuru73xaJlRDiXHZLGgUFBYiOjkZkZCRkMhlmzJiBvXv3Cs7Z\nt28f5s6dCwAYPXo0GhsbUVtbC19fX8hkMqhUKuj1eqhUKn7pxVuZuZaGm5sbv+3j44Py8nL4+fmZ\nPBjl5uYGrVbrsFgJIbcnuz0RXlVVJVgyUaFQ4MSJEzc9p6qqCsOHD8fSpUsxaNAgeHp6IjMzE+PH\njzf7eVasWMF/nJGRgYyMDJt+HY7UOWkYL++qUqmgVquhUqmgUChw9uxZqNVqeHp6AqCkQQix7MiR\nIzhy5IhN7mW3pGFtTRlzNd1LS0vxxhtv8O+qH374YXz88ceYPXu2ybnGSaOvM5c0uKSgVCoRERGB\nhoYGMMYQEBCAGzduUNIghNxU5zfUK1eu7PG97NY9xa2py6moqIBCoejynMrKSoSHh+Onn37Cb37z\nGwQFBUEqlWLKlCk4fvy4vUJ1CTqdDnq9XjCGodVq+ZaGWq2Gn58fvLy80NTUxCcNLulS0iCEOILd\nkkZ6ejpKSkpQXl4OrVaLnJwcZGVlCc7JysrChx9+CADIz8+Hv78/QkJCEBsbi/z8fKjVajDGcPDg\nQSQkJNgrVJfAlRAxbqHpdDp+TEOtVsPDwwPBwcG4fv065HI5DAYDP2OKkgYhxBHs1j0llUqxZcsW\nZGZmwmAwYMGCBYiPj8e2bdsAdKyTO2nSJOTm5iI6OhpyuZxf6zY1NRWPPfYY0tPTIRaLMXz4cMEa\nurei5uZm+Pr6CvZxYxqMMX78gksaQ4cOFXRRubu7U9IghNidXUujT5w4ERMnThTsy87OFmxv2bLF\n7LXPP/88nn/+ebvF5mrMJQ2ue0qr1UIikUAqlfJJA+iolX/58mWEhYXBzc3NYkFDQgixFXoi3EU0\nNzfDx8dHsI/rnjKeJRUUFIT6+noA4LuouKm51NIghNgbJQ0X0VX3lHHS8Pf356vdikQiyOVyqFQq\nShqEEIegpOEiWlpa4OfnJ9jHdU9ZShoA+HU2KGkQQhyBkoaLsLZ7qnPS8PLyopYGIcRhKGm4iM7d\nU9y6GSKRCBqNBh4eHgAoaRBCnIuShovonDS4rilukFss7vhRyeVyaLVavjihVCqFTCZDe3s7JQ1C\niN1R0nARnZOGua4poKPl4e/vj6amJn6fl5cX9Ho9tFqt2bIshBBiK5Q0XIS5pNF5EJwTEBBg0kXV\n1tYGiURCCzERQuyKkoaLaGlpEQyEm5s5xfHz88ONGzf4beMZVLSmBiHEnihpuACdTof29nZBsUJL\n3VMA+PIhHC8vL6jVahoMJ4TYHSUNF9Da2govLy+TYoVcCRHjZAIA/fv3x7Vr1/htiUQCmUwGqVRK\nSYMQYleUNFyAUqmEt7e3YJ9Wq4VUKoVEIjFZmyQkJAS1tbWCfXK5HBKJhJIGIcSuKGm4AK4sujGd\nTscXKezMXNLgWiqUNAgh9tRlldvCwkLs2rULx44dQ3l5OUQiESIiInDnnXdi1qxZSEtLc1SctzSl\nUtmtpBEaGoqamhowxvhWCCUNQogjWEwakyZNQkBAALKysrBo0SKEhYWBMYarV6+ioKAAGzZsQGNj\nI77++mtHxntLam1tNUkaxuXQO/P29oZYLEZLSws/Tdfd3Z2SBiHE7iwmjQ8++AAhISEm+wcPHozB\ngwdjxowZgsFY0nPmkgY3e8pc0gA6uqhqamr4pCGTyfiSI4QQYi8WxzSMp3RyS4py8vPzAXTM4iG9\nZ657SqvVQiwWW0waoaGhgnENkUgENzc3k58VIYTYksWkMXPmTP7j3/zmN4JjCxcutF9EtyGVSmUy\ne0qn00EkEnXZ0ug8GE5JgxBib1bNnqJ6RvZlaSD8ZkmjpqZGsM/d3Z26p1xMRkYGCgsLAQD3338/\nmpube3yvLVu2YPv27QCAzz77DImJiZBIJPz9OWvWrEFMTAzi4uKwf/9+fv/PP/+M5ORkxMTE4E9/\n+hO/X6PR4JFHHkFMTAzGjBmDy5cvAwBqa2sxadIkq2K77777kJqaisTERCxYsAA6nQ4AsHXrVgwb\nNgxpaWkYO3YsTp8+bfb6Dz74AMnJyUhJScHEiRP51SmtvZ44Dk25dQGWBsK7ShrcDCpjnp6eNBDu\nYoyfsfn6669NVme0FmMM7733Hh599FEAQHJyMv71r3/hzjvvFJxXXFyMnJwcFBcXIy8vD4sWLeLf\n9C1cuBDvvfceSkpKUFJSgry8PADAe++9h6CgIJSUlOCZZ57B8uXLAXS8MQkICDBJSubs2bMHp06d\nwrlz59DU1IScnBwAwOzZs3HmzBkUFRXhL3/5C5YuXWpyrVarxbJly3D06FGcPn0aw4YNw5YtW6y+\nnjiWxaRRWVmJJUuW4KmnnkJVVRX/MbdNbMfSQDgAq8c0AMDDw4OSRi+Vl5cjLi4O8+fPR2xsLGbP\nno39+/fjt7/9LYYOHYqTJ08C6PiZPf744xg9ejSGDx+Offv2AQDUajVmzJiBhIQETJkyBWq1mr93\nZGQkGhoaAAAPPfQQ0tPTkZSUhHfffZc/x9vbG3/961+RmpqKsWPH8pNNfvjhB8TFxfG/D3FxcRg6\ndKhJ/Hv37sXMmTMhk8kQGRmJ6OhonDhxAlevXkVLSwtGjRoFAHjsscfw5ZdfAgD27duHuXPnAgCm\nTp2KQ4cO8ffLysrCrl27bvp947pXdTodtFot+vXrBwCCempKpZLfb0wqlSIgIABKpRKMMTQ3NyM8\nPNzq64ljWZw9tX79ev5d0ogRIwTH0tPT7RvVbaYnSSMwMBAtLS3QaDR8mRFqadhGaWkpPv/8cyQk\nJGDkyJHIycnBDz/8gH379mH16tX417/+hVWrVuGee+7B+++/j8bGRowePRrjx4/H1q1b4e3tjeLi\nYpw9exbDhw/n72vc6nj//fcREBAAtVqNUaNGYdq0aQgICIBKpcLYsWPxt7/9DcuXL8e7776LF198\nEd9//z1Gjhx509irq6sxZswYfluhUKCqqgoymQwKhYLfHx4ezr/5q6qqwsCBAwF0/L75+fmhoaEB\ngYGBGDVqFDZu3Mjf+4knnrA4zT4zMxMnT57EhAkTcN999/H73377bWzcuBGtra04fvy4yXVisRib\nNm1CUlISvL29MXToULz11ltWX08cy2JLY968eZg7dy7mzp2LefPmYd68ecjKyuL3EdtRKpUmS71y\na2NYShoSiQTBwcGCac/cuhqkd6KiopCYmAiRSITExESMHz8eAJCUlITy8nIAwP79+7F27VqkpaVh\n3Lhx0Gg0uHLlCr777jtBF9KwYcPMfo5NmzbxrYmKigqUlJQA6JjMcP/99wPoeLPGfb4rV64gNDTU\njl+1eWFhYXwMAwYM6PK5rG+//RZXr16FRqPBjh07+P2LFi3CpUuXsHHjRjz++OMm1zU3N2PJkiU4\nffo0qqurkZycjDVr1lh9PXEsi0lj5cqVOH/+PICOgbJx48YhOjoaISEhOHDggMMCvB20tLQIZk8x\nxvgXf4lEYvE6S0mDJi70jnGBSLFYDDc3N/5j46T8xRdfoKioCEVFRXy3FnDziSNHjhzBoUOHkJ+f\nj1OnTiEtLY2f9SaTyQSf2/jzWfNzDQ8PR0VFBb9dWVkJhUKB8PBwVFZWmuznrrly5QoAQK/Xo6mp\nCYGBgfzn7Fz7rCvu7u6YOnUq341n7JFHHjE7PnL+/HlERUUhKioKAPDwww+bbVFYup44lsWkkZOT\nw/8R7NixA4wxXL9+HUePHsVf/vIXhwV4q9Pr9dBoNPDy8hLs454G7+oPNiQkRJA0PD09IRKJqLXh\nAJmZmdi8eTO/XVRUBAC488478cknnwAA/vvf/+LMmTMm1zY3NyMgIAAeHh745Zdf+OeeuhIREWEy\n8YFjnEyysrKwe/duaLValJWVoaSkBKNGjUJoaCh8fX1x4sQJMMawc+dOPPjgg/w1XMtgz549uOee\ne/j7Xb16FREREV3G1traiqtXrwLo+N396quv+BJDXAsK6JgIYK7lNXjwYPzyyy+oq6sDABw4cAAJ\nCQlWX08cy+KYBleWAgDy8vIwY8YMSCQSxMfH04uSDXEVbo2TA7cAk6WuKU7//v1x/fp1flsikUAi\nkUCtVgvesZLu6Zyojbe5j1966SU8/fTTGDZsGNrb2zF48GDs27cPCxcuxPz585GQkID4+Hiz43/3\n3Xcftm7dioSEBMTGxmLs2LEWPxe3fccdd/AzigDgX//6F5YsWYK6ujrcf//9SEtLwzfffIOEhARM\nnz4dCQkJkEqlePvtt/l7vP3225g3bx7UajUmTZrEjzssWLAAc+bMQUxMDIKCgrB7927+8xQUFPAz\ntCyNabS2tuLBBx+ERqMBYwyZmZl8N9Jbb72FgwcPQiaTITg4GB988AF/XVpaGoqKihAcHIzVq1dj\n3LhxEIvFiIyM5KcWd3U9cRJmwejRo9mZM2fYtWvXWEBAACstLeWPDR061NJlDtVF+H1GRUUFe+qp\npwT7bty4wXbu3MnOnz/f5bXff/89+/vf/y7Yt337dlZdXW3zOIlztbe3s9TUVKbRaBz6eWfNmsUK\nCwsd+jmJ/fXmtdNi99Qbb7yBadOmITY2Fs888wwGDx4MoKOJaDwjhPRO5/EMoGPmlFQq7XZLA+jo\nE29tbbV5nMS5RCIRnnjiCXz88ccO+5zXrl1DY2MjVbMmAhaTxpgxY3DhwgU0NDTgpZde4vfff//9\nVs3bBjq6teLi4hATE4N169aZPWfJkiWIiYlBSkoK3y8MAI2NjZg2bRri4+ORkJBgVb9vX2Rp5pS1\nScNcKRHjZwPIrWPRokWYP3++wz5f//79qYo1MWHxVem1116DSCQyO2NDJBLh2Wef7fLGBoMBixcv\nxsGDBxEeHo6RI0ciKysL8fHx/Dm5ubm4dOkSSkpKcOLECSxcuJBPDn/6058wadIk7NmzB3q9/pZ9\n92yppWGpLLoxX19faLVawTrilDQIIfZksaXx3HPPYefOnaivr4dSqRT8a2lpuemNCwoKEB0djcjI\nSMhkMsyYMQN79+4VnGP8JOro0aPR2NiI2tpaNDU14bvvvuMH07gHjm5Flloa1iQNkUiEoKAgQUVi\nd3d3KlroQuxVe+q5555DfHw8UlJSMGXKFDQ1NfHn9cXaUy+++CIGDRpk8rcAAJ9++ikSExORlJSE\n2bNnWxUPsR+Lr0rcqn25ubkYPnw4Zs6ciXvuuQdisXXlqoyfMgU6nkw9ceLETc+prKzkH1ybP38+\nTp8+jREjRmDTpk2CaamcFStW8B9nZGQgIyPDqvhchbmWhrVJAwACAgLQ0NCAAQMGAOgoJUJJw3V0\nrj3VU+z/rz3FPf9w7733Yt26dRCLxfjzn/+MNWvWYO3atYLaU1VVVRg/fjxKSkogEon42lOjRo3C\npEmTkJdDYcCwAAAgAElEQVSXh/vuu09QeyonJwfLly/H7t27BbWnbjaOuWfPHv73eNq0acjJycGj\njz6K2bNn48knnwQA/Pvf/8bSpUtx8OBBk+sffPBBPPXUU4iJiRHsLykpwdq1a3H8+HH4+fnx03JJ\n9xw5cgRHjhyxyb0sZoDU1FSsW7cORUVF+P3vf499+/YhMTGRr7FzM9Y+ENS5+4t7zqCwsBCLFi1C\nYWEh5HI51q5da/b6FStW8P/6WsIALLc0ulpLw1hgYCBfzwjoeFaDKt32XF+pPTVhwgT+Ddzo0aP5\nB/f6Yu0pAPyzJJ29++67WLx4Md/TQLWneiYjI0PwWtkbN202XL9+HUVFRThz5gwUCgWCg4OtunHn\nJ1MrKioEtW/MnVNZWYnw8HAoFAooFAq+1s60adNu2SdBLY1p9CZpUP2p3iktLcWyZcvwyy+/4MKF\nC3ztqQ0bNmD16tUAwNeeOnHiBA4fPoznnnsOKpUK//jHP/jaUytXrsTPP//M37dz7amffvoJJ0+e\nxObNm/kuRq721KlTp3DnnXfyCaWr2lPvv/8+341UXV0t+Dvjak913m9N7Smg48X82LFj/L25Eifm\nZGZmIiQkBJ6enia1p6Kjo/Hss88KyoNYo6SkBBcuXMAdd9yBsWPH4ttvv+3W9cT2LCaN9957D5mZ\nmZg+fToYY/j0009x4MABwYNIXUlPT0dJSQnKy8uh1WqRk5ODrKwswTlZWVn48MMPAXSsBujv74+Q\nkBCEhoZi4MCBuHjxIgDg4MGDSExM7OnX6NJ6UhbdWOek4eXlBZ1OR6VEeqEv1Z5atWoV3NzcMGvW\nLFt+C3iOqD3VFb1ej0uXLuHo0aPYtWsXnnjiCcH4DXE8i69KTzzxBJKSkhAREYFvv/1WkOFFItFN\nu6mkUim2bNmCzMxMGAwGLFiwAPHx8di2bRsAIDs7G5MmTUJubi6io6Mhl8sFT3u++eabmD17NrRa\nLYYMGXLLPgmqVqtNxmq4loK1SePChQv8tqenJ/R6Pdrb27usW0Us607tqc598ED3ak95eHhg3Lhx\nPao9tX37duTm5gq6k3pTe2rAgAE2qz114sQJk8KmjzzyCD++YS2FQoHRo0dDIpEgMjISQ4cOxaVL\nl0wqbxPHsfiqdPjwYQAwO+3W2l+iiRMnYuLEiYJ92dnZgm3j0gjGUlJSzBY9u9WoVCqzSUMqlVr1\not+5peHm5gaDwQCtVstPwyW2x9WeevPNNwF01J5KS0vja0+NGzfOrrWn8vLysH79ehw9ehQeHh78\n/qysLMyaNQvPPvssqqqq+NpTIpGIrz01atQo7Ny5E0uWLOGv2bFjB8aMGdPj2lPNzc0ICwvja0/d\ne++9ADq6l7jE2pPaUZMnT8auXbswb9481NXV4eLFi/yDxsQ5LCaNrgaVCwoK7BHLbUmlUpm8uGu1\nWnh5eVmVnDsnDXd3dxgMBuh0OkoaPdQXak899dRT0Gq1mDBhAgBg7NixePvtt/tk7SkAeP7557Fr\n1y6o1WoMHDgQTzzxBF5++WVkZmZi//79/NK2GzZsQEBAQBc/PWJ3luqLGAwGtmfPHrZu3Tr29ddf\nM8YYO3nyJJswYQJLSUnpcd0SW+oi/D5j1qxZTKVSCfbl5OSw48ePW3V9S0sLe/TRR/ltjUbDtm3b\nxurq6mwaJ3Euqj1FbKk3r50WB8L/8Ic/4O2338aNGzfwt7/9DVOnTsXcuXOxaNEinDp1ynFZ7RbG\ndSMZdy8AHbOnjPvVu+Lp6Ym2tja+C1Emk6G9vZ2m3d5iqPYUcRUWu6fy8/Nx5swZiMVitLW1ITQ0\nFKWlpQgKCnJkfLc0tVoNDw8Pk+4QrVZrddLgHgLUaDT8vaho4a1p0aJFDv18VHuKmGOxpSGTyfiH\nhzw8PBAVFUUJw8bMzZwCutfSADqm2Ro/RObm5kZJgxBiFxaTxi+//ILk5GT+34ULF/iPafUs22ht\nbTVJGu3t7Whvb+9W0vD09BQkDXd3d6hUKpvFSXrOXrWnOK+99hrEYrFgMoQr1Z46duwYhg8fDplM\nhs8//9zi9ZZiA6j2lKux2D3FrQ9O7Me4Oi2HW0ujOyvveXp6CpIEV0qkvb3d6lphxD7sVXsK6Kiy\ncODAAcGUWFerPRUREYEdO3Zgw4YNXV5vKTaqPeV6LL6iREREIDIy0uI/wLqF7ollKpXK7NPg1hYr\n5HCD4RzuYTQaDO++vlJ7CgCeffZZ/P3vfxfE72q1pyIiIpCcnNzlm5euYqPaU67H4k8yIyMD69ev\n50t5GLtw4QLWrVuHu+66y67B3eosPaPRk6Rh3NLg1nenpNEzfaH21N69e6FQKEy6il2x9tTNVFVV\nWYyNak+5HotJY//+/QgKCsIf//hHhIWFYejQoYiJiUFYWBgWL16MkJAQsyWOifXMdU/1JGl0Hgh3\nd3cHY4xKpPeQq9eeUqlUWL16NVauXMnfz16tflvUnuoNnU5HtadcjMVXJnd3dzz++ON4/PHHYTAY\n+L7Efv36UU0jGzHXPWXtqn3GPDw8TJKGWq2mlkYPuXrtqdLSUpSXlyMlJQVARx2pESNG4MSJEy5d\ne8rSfbqKbeDAgVR7ysVYbGmo1Wq8/vrr+OMf/4h//vOfCAoKQkhICCUMG7LUPSUWi7v1fTY35dZg\nMFBLw4642lMcrhwGV3sKgN1qTyUnJ6O2thZlZWUoKyuDQqFAYWEhQkJCkJWVhd27d0Or1aKsrIyv\nPRUaGsrXnmKMYefOnXjwwQcB/K/2FIAe1566evUqAPC1pzo/EMgYs5hMw8LCLMY2efJkfvEgqj3l\nGiwmjblz5/LT4HJzc7F06VJHxnVbMNc9pdForF5Lg2Nuyq3BYKCWRg9ZW3tKp9Nh2LBhSEpKwiuv\nvAKgYxaQUqlEQkICXnnlFYu1p/R6PRISEvDCCy9YXXvqp59+umm8xrWnJk6caFJ76ve//z1iYmIQ\nHR0tqD1VX1+PmJgYvPHGG4IFzzrXnjI3psHVnkpJScHw4cMxaNAgvvbUyZMnMXDgQOzZswfZ2dlI\nTk7mrzNOLJZiy8zMRFBQEBITE3H33XdT7SkXIGIW0n9ycjLOnj0LoOPdw8iRI/l3U67CXAXevuSt\nt95CbGws32cOdCyze+XKFUyePNnq+3z99deoqanBggULAHT0f58+fRrh4eFITU2labe3AMYYhg8f\njhMnTvDdZY4we/ZsLFu2jEqJ3GJ689pp8dXE+J1ud971EuuZKxei1Wq79YwGYDp7ysPDAxqNBu7u\n7tRFdYug2lPEVVjMBmfOnBGs76tWq/ltkUjUqydbSYe2tjaTpNHW1tajpNG5e4pLGhqNxmypEtL3\nUO0p4gosJg2DweDIOG5L3Au7Ma1W2+3uh85Jw8PDA21tbfz/hBBiK9TZ7UTmkoZGo+l20ggICBCU\nV3Bzc4NOp4ObmxsNhhNCbIqShhOZSxrci313KBQKNDQ0QKlUAujoPnRzc4NYLKakQQixKUoaTmSp\ne6o7FW6BjjU1Bg8ejEuXLvH76KlwQog9UNJwIm7hJGM6nc5knzViYmL4UhRAx7iGwWBAe3s7jU8R\nQmyGkoYTmWtp6PX6brc0gI56SVeuXOG33d3d+VYLdVERQmyFkoYTWUoanZ8St4afn5+gkBv3jAYN\nhhNCbImShpNwK5sZPzjJGINer+9R95Sfn5/g2RnjB/y0Wm3vAyaEEFDScBpLrQzjqqrd4evra9LS\n4KbvUkuDEGIrlDScxNwgOLeWRk8qCfv6+kKpVPKD3sbdU9TSIITYCiUNJ7E03banSUMikcDLy4t/\nVoO6pwgh9kBJw0nMJQ21Wg2JRNKtRW+MGQ+Gd+6e6svVgAkhrsOuSSMvLw9xcXGIiYnBunXrzJ6z\nZMkSxMTEICUlxaT0usFgQFpaGn73u9/ZM0ynsJQ0elNR2NfXlx8M57qnuFYLPatBCLEFuyUNg8GA\nxYsXIy8vD8XFxdi1axfOnz8vOCc3NxeXLl1CSUkJ3nnnHSxcuFBwfNOmTUhISOjxO29XZqnuVHcr\n3Bozbmlw3VMikYi6qAghNmO3pFFQUIDo6GhERkZCJpNhxowZ2Lt3r+Ccffv28esIjx49Go2Njait\nrQXQsU5wbm4ufv/739+SXStcFdrO+3qTNDq3NLhZUzQYTgixFbutrlRVVYWBAwfy2wqFAidOnLjp\nOVVVVQgJCcEzzzyD9evX33TdjhUrVvAfZ2RkICMjwybx25u5Eui2aGk0NjYC+F/SYIzRtFtCbnNH\njhzh11rvLbslDWu7lDq3Ihhj+Oqrr9C/f3+kpaXd9As1Thp9ibmWRk/W0jAWFhaGwsJCAB2zqcRi\nMXQ6HXVPEXKb6/yGeuXKlT2+l926p8LDw1FRUcFvV1RUQKFQdHlOZWUlwsPDcfz4cezbtw9RUVGY\nOXMmDh8+jMcee8xeoTqFrRZgMjZgwABUVVXx29y4BnVPEUJsxW5JIz09HSUlJSgvL4dWq0VOTg6y\nsrIE52RlZeHDDz8EAOTn58Pf3x+hoaFYvXo1KioqUFZWht27d+Puu+/mz7tVKJVKyOVywT5bJI3q\n6mq0t7cDoPpThBDbs1v3lFQqxZYtW5CZmQmDwYAFCxYgPj4e27ZtAwBkZ2dj0qRJyM3NRXR0NORy\nOT744AOz97oVZ0+1trYiJCREsI/rSuopuVwOT09PNDQ0oF+/fnxLIyAggFoahBCbsFvSAICJEydi\n4sSJgn3Z2dmC7S1btnR5j7vuugt33XWXzWNzNqVSiSFDhgj29TZpAB1dftXV1ejXrx8/GC6RSMAY\ng8Fg6NHT5oQQwqEnwp2ktbXVpHuqpwswGQsLC0NNTQ2A/3VPccu/UhcVIaS3KGk4iVKphI+Pj2Bf\nT9fSMObv7y94wI9b7pUGwwkhtkBJw0nMDYT3dNU+Y8bPanh6ekKtVgMATbslhNgEJQ0naW1thbe3\nN7/NreXd2+4p48WYjJMGdU8RQmyBkoYTMMZMWhp6vR7t7e29mnILCOtPeXl5CZIGtTQIIb1FScMJ\nNBqNyQp9Go0G7e3tvapyCwhX8PP09IRKpQJA3VOEENugpOEE5mZOaTQaSKXSXj+TYtzSoO4pQoit\nUdJwAqVSKRjPAHq/lgbH29sbKpUKer2ef7jPuAWj0+l6/TkIIbcvShpOYC5p9LYsOkcikcDHxwct\nLS0Qi8Xw8PCAWq2GSCQSTMElhJCeoKThBOam23LdU7ZgqYuKkgYhpLcoaTiBSqWCl5eXYB9XjdYW\nLD2rYfwxIYT0BCUNJ7CUNGzRPQUAgYGBqK+vB9Ax7ZabQUUtDUJIb1HScAKVSmW2e6q3T4NzgoKC\n+KTRuaVBSYMQ0huUNJzAXEtDp9PZrHsqKCgIDQ0NAIRJQyaTwWAwQK/X2+TzEEJuP5Q0nMBShVtb\ntTQsdU/RDCpCSG9R0nACSy2N3tad4ljqngJoXIMQ0juUNJygc9JgjNm0pWGpe4rbpqRBCOkpShpO\n0HkgXK/XgzFmszENX19ftLW1QaPRCIoWAuAf9iOEkJ6gpOEEnVsatqpwyxGJRPy4BtfSYIwBoJYG\nIaR3KGk4QWtrq0nSsMUCTMZCQ0NRU1MDiUQCqVTKFyt0c3ODTqeDwWCw2ecihNw+KGk4gbmWBldg\n0FbCwsJw9epVAMJxDZpBRQjpDUoaDqbX62EwGAStCr1eb9PZU8D/WhpAx7Tb1tZW/hglDUJIT1HS\ncDCVSgVPT0/BuhkajQYGg8FmYxpAR0uDSxre3t6CpEHjGoSQnqKk4WDmHuxTqVRwc3Pr9QJMxkJD\nQ/nuKblcbpI0uAf+CCGkOyhpOJi5B/vUarVNu6YAICQkBHV1dTAYDPD29oZSqeSPcU+JczOqCCHE\nWpQ0HKylpQU+Pj6CfWq12qYzp4COOlN+fn6oq6sz6Z7iqunSKn6EkO6ipOFg5pKGRqOBp6enzT9X\nSEgIamtrIZfLBS0NkUgkqElFCCHWoqThYEql0uFJo3NLAzAd5yCEEGvYPWnk5eUhLi4OMTExWLdu\nndlzlixZgpiYGKSkpKCoqAgAUFFRgXHjxiExMRFJSUnYvHmzvUN1iJaWFpP1wblyH7YWEhKCa9eu\nwdPTk5+hxaGWBiGkJ+yaNAwGAxYvXoy8vDwUFxdj165dOH/+vOCc3NxcXLp0CSUlJXjnnXewcOFC\nAB397q+//jrOnTuH/Px8vPXWWybX9kWdWxqMMWi1Wru2NEQikUnLggbDCSE9YdekUVBQgOjoaERG\nRkImk2HGjBnYu3ev4Jx9+/Zh7ty5AIDRo0ejsbERtbW1CA0NRWpqKoCO5wzi4+NRXV1tz3AdovOY\nBvewnz2TBgCTcQ1uiq9Wq7X55yWE3LrsmjSqqqowcOBAfluhUKCqquqm51RWVgrOKS8vR1FREUaP\nHm3PcB1CqVQKuqfsmTTCwsJQXV0NxpjZcQ3qoiKEdJfUnje39mG1zl0kxtcplUpMmzYNmzZtMhkL\nAIAVK1bwH2dkZCAjI6NHsTqKuZaGrUuIcHx8fODm5oaGhgazA99c0ggICLD55yaEuI4jR47gyJEj\nNrmXXZNGeHg4Kioq+O2KigooFIouz6msrER4eDiAjucIpk6dikcffRSTJ082+zmMk0Zf0Hkg3B7F\nCo0NHDgQV65cgbe3N5qbmwXH5HI5331FCLl1dX5DvXLlyh7fy67dU+np6SgpKUF5eTm0Wi1ycnKQ\nlZUlOCcrKwsffvghACA/Px/+/v4ICQkBYwwLFixAQkICnn76aXuG6VCdWxo6nQ5ardZuSWPQoEG4\ncuVKly0NGgwnhFjLrklDKpViy5YtyMzMREJCAh555BHEx8dj27Zt2LZtGwBg0qRJGDx4MKKjo5Gd\nnY23334bAPDDDz/go48+wn/+8x+kpaUhLS0NeXl59gzX7gwGA9RqtWB6bVtbG8RiMaRS+zT6uKTR\nuZQI0DFDTSwW02A4IcRqdu2eAoCJEydi4sSJgn3Z2dmC7S1btphcd8cdd6C9vd2usTmaUqmEl5cX\nJBIJv0+lUtm8hIixQYMGIS8vz+LDfNx+e8ZACLl10BPhDtTU1AR/f3/BPnvUnTIWERGBqqoqyGQy\ntLW1Qa/XC46ba4EQQogllDQcqKmpCX5+foJ9bW1tdhvPAAB3d3eEhISgqqoKvr6+aGpqEhynpEEI\n6Q5KGg5kKWnY4xkNY4MHD0ZZWRkCAgLQ0NAgOObl5QWNRmPSAiGEEHMoaTiQuaRhr2KFxqKiovDr\nr78iICAAN27cEBzjSoxQa4MQYg1KGg5kKWnYo1ihsQEDBuDq1atmkwbQ8RAgJQ1CiDUoaThQ54Fw\nxhjUarVJIrE1rtptYGCg2aTh7e2NlpYWu8ZACLk1UNJwoKamJvj6+vLb3IN9xvvsITg4GHV1dfDx\n8UFzc7OgRDrQMe22ra3NZD8hhHRGScOBOndPaTQaaLVak0WZbM3d3R1yuZwvYdK5nIhYLIaXlxct\nykQIuSlKGg7UuXuqpaXFbhVuO+PKpJubQQVQFxUhxDqUNBykvb0dN27cECSNxsZGyOVyq6sB90b/\n/v1x7do1BAQEoLGx0eQ4DYYTQqxBScNBuPLkxg/yNTY2mi33bg+hoaGoqqqy2NKQy+VQqVS3XOkW\nQohtUdJwkKtXryIsLEywr6Wlxe4zpzhJSUk4e/YsAgMDUVdXZ3JcIpHA09OTuqgIIV2ipOEgnZMG\nYwwqlcqkFpW9xMbGoqqqCm5ublAqldBoNCbn+Pn5mQySE0KIMUoaDtI5aXAzp+w93ZYjk8mQkJCA\n//73vwgODsa1a9dMzjFXm4oQQoxR0nCQmpoahIaG8tttbW0OTRoAEBMTg9LSUn4mVWdeXl4wGAxm\nWyGEEAJQ0nCYzi2NtrY2tLW12f0ZDWODBw/Gr7/+ajFpiEQi+Pr6UhcVIcQiShoOYDAYUFtbK0ga\n3PRWRy5+xFW75abfmlvm1c/Pj7qoCCEWUdJwgPr6evj4+AgSBDfd1hHPaHD8/f0hk8mgUqkglUrN\ntih8fX3R0tJCU28JIWZR0nCA6upqk+m25lbxc4Rhw4bh+++/t9hFJZVK4enpSQ/6EULMoqThAFev\nXsWAAQP4bYPBgObmZoSEhDg8lsmTJ+Orr75CUFAQLl++bPYcGtcghFhCScMBzE23ValUTkkagwYN\nQmpqKs6ePYsrV65ArVabnBMQEID6+npotVqHx0cIcW2UNBygtLQUgwYN4rdVKhVaW1sRHBzslHjm\nzJmDo0ePwsfHBxcuXDA57unpif79+6O8vNzsYDkh5PZFScPOmpqacOXKFSQmJvL7rl27Bi8vL7i5\nuTklpoCAAMyZMwcFBQUoLi42mxhCQ0PBGDM77kEIuX1R0rCzn376CSkpKZDJZPy+uro6BAUFOTEq\n4K677kJtbS0MBgMqKipMjotEIkRFRaG2tpbW2SCE8Chp2NnZs2eRlpbGbzPGUFdX55TxDGMSiQSx\nsbHw9fVFfn6+2daGm5sbIiIiUFpaCp1O54QoCSGuhpKGnV24cAFxcXH8dktLC5RKJRQKhROj6pCQ\nkID6+npIpVJcvHjR7Dn+/v4IDg5GaWkpPbtBCKGkYU83btyASqUSTLe9fPkyNBqN07unACA9PR3f\nffcd/P39UVBQYHG2VGhoKNzc3HD58mUaGCfkNkdJw47OnDmD2NhY/qlvlUqFc+fOITU1FRKJxMnR\ndUy/ff7557Fz5040NTXhP//5j9mkIBKJEBERAZ1OhzNnzqCiooLGOQi5Tdk1aeTl5SEuLg4xMTFY\nt26d2XOWLFmCmJgYpKSkoKioqFvXuqrDhw9jz5492L59Ox544AF+/9mzZ6HT6TB8+HCHx3TkyBGz\n++Pj4/HSSy+hoqIChw8fxpkzZ8x2Q0kkEgwdOhSxsbGQSCQoLS3FpUuXoFKpbB6TM1FM1nPFuCgm\n+7Nb0jAYDFi8eDHy8vJQXFyMXbt24fz584JzcnNzcenSJZSUlOCdd97BwoULrb7WlW3fvh0nT57E\nqlWrMGzYMABAbW0tzp49i4yMDKe0Mrr6xR08eDBefvlltLS0YOPGjXj11VeRl5dntqihh4cHBgwY\ngKSkJPj6+qKkpAQXL15EbW1tt5eLdcU/JorJeq4YF8Vkf1J73bigoADR0dGIjIwEAMyYMQN79+5F\nfHw8f86+ffswd+5cAMDo0aPR2NiImpoalJWV3fRaV6HT6aBWq+Ht7Y0ff/wRpaWlOHfuHF544QUM\nGDAAOp0Oly5dwvHjx5GYmIiBAwc6O2Sz5HI53nzzTVy6dAk5OTnYuXMn9u/fD6lUCoVCgeHDh0Mu\nl6O+vh79+vVDSEgIvLy8EBERAb1eD5VKhevXr/NrcXh4eEAmk0EqlUIkEsHd3R1SqRRisRhisRgi\nkQgajQatra2QSCQQiUT8fg63nxDiOuyWNKqqqgQvkAqFAidOnLjpOVVVVaiurr7ptc6k1+vx3HPP\nwWAwoKGhASKRCFqtFn5+flAoFAgICMDJkyeRn58PrVYLuVyO9PR0JCcnOzv0LnHTcF9++WXcuHED\nP/74I5qbm3H27Fls2bIFjDF4enry4xlSqRReXl5Qq9Xw8vICIFxcSq/XQyaTgTEGqbTjV02n08Fg\nMEAqleLkyZOora3ln2GRSCQwGAwYNGgQ+vfvD7FYjPb2djDG4O7uDjc3N7S3t0Ov10Or1UIikfD3\nZ4yZJBguEXHJyPg4dw13HXfs3Llz+PTTT/lzLDE+JhaL+c9nfG/jfZaSn/G5xnEbn19cXIzPPvvM\n7L17m1S5+0VHR/Pr1XP37+rzMMbQ0NCAkpKSLr8+4/twzN2rq+Pd0dDQgEuXLlkcm3MGLqbOuBj9\n/PzQv39/R4fVY3ZLGtb+gHo7G8dV34keOnTI2SGYWLlypbNDMHHmzBlnh2Biz549zg7BBJc0XM2b\nb77p7BBMUEz2ZbekER4eLnjSuKKiwuTZhM7nVFZWQqFQQKfT3fRaoPcJhxBCSPfYbSA8PT0dJSUl\nKC8vh1arRU5ODrKysgTnZGVl4cMPPwQA5Ofnw9/fHyEhIVZdSwghxPHs1tKQSqXYsmULMjMzYTAY\nsGDBAsTHx2Pbtm0AgOzsbEyaNAm5ubmIjo6GXC7HBx980OW1hBBCnIz1Ud988w2LjY1l0dHRbO3a\ntU6J4cqVKywjI4MlJCSwxMREtmnTJsYYY/X19Wz8+PEsJiaGTZgwgd24ccPhsen1epaamsoeeOAB\nl4jpxo0bbOrUqSwuLo7Fx8ez/Px8p8e0evVqlpCQwJKSktjMmTNZW1ubU2KaP38+69+/P0tKSuL3\ndRXH6tWrWXR0NIuNjWXffvutw2JatmwZi4uLY8OGDWMPPfQQa2xsdHpMnA0bNjCRSMTq6+tdIqbN\nmzezuLg4lpiYyJ5//nmHxmQprhMnTrCRI0ey1NRUlp6ezgoKCnoUV59MGnq9ng0ZMoSVlZUxrVbL\nUlJSWHFxscPjuHr1KisqKmKMMdbS0sKGDh3KiouL2XPPPcfWrVvHGGNs7dq1bPny5Q6P7bXXXmOz\nZs1iv/vd7xhjzOkxPfbYY+y9995jjDGm0+lYY2OjU2MqKytjUVFRrK2tjTHG2PTp09n27dudEtOx\nY8dYYWGh4A/cUhznzp1jKSkpTKvVsrKyMjZkyBBmMBgcEtP+/fv5z7V8+XKXiImxjjdvmZmZLDIy\nkk8azozp8OHDbPz48Uyr1TLGGLt27ZpDY7IU11133cXy8vIYY4zl5uayjIyMHsXVJ5PG8ePHWWZm\nJr+9Zs0atmbNGidG1OHBBx9kBw4cYLGxsaympoYx1pFYYmNjHRpHRUUFu+eee9jhw4f5loYzY2ps\nbGRRUVEm+50ZU319PRs6dChraGhgOp2OPfDAA2z//v1Oi6msrEzwB24pjtWrVwta1pmZmezHH390\nSEzGvvjiCzZ79myXiGnatGns9OnTgqThzJgefvhhdujQIZPzHBmTubhmzJjBcnJyGGOMffLJJz3+\n+Q2mf2MAAAZTSURBVPXJ2lOWnu9wpvLychQVFWH06NGora3lS5+HhIQ4fCGjZ555BuvXr+efHwDg\n1JjKysoQHByM+fPnY/jw4XjiiSfQ2trq1JgCAwOxdOlSDBo0CAMGDIC/vz8mTJjg9J8dx1Ic1dXV\ngpmEzvrdf//99zFp0iSnx7R3714oFAq+8gLHmTGVlJTg2LFjGDNmDDIyMvDTTz85PSYAWLt2Lf87\n/9xzz2HNmjU9iqtPJg1XezZDqVRi6tSp2LRpE3x8fATHbPEAVnd89dVX6N+/P9LS0ixOSXZ0THq9\nHoWFhVi0aBEKCwshl8uxdu1ap8ZUWlqKN954A+Xl5aiuroZSqcRHH33k1JgsuVkcjo5x1apVcHNz\nw6xZsyye44iYVCoVVq9eLXj+yNLvvKNiAjp+32/cuIH8/HysX78e06dPd3pMALBgwQJs3rwZV65c\nweuvv47HH3+8R3H1yaRhzTMgjqLT6TB16lTMmTMHkydPBtDxzrCmpgYAcPXqVYc+7Xn8+HHs27cP\nUVFRmDlzJg4fPow5c+Y4NSaFQgGFQoGRI0cCAKZNm4bCwkKEhoY6LaaffvoJv/nNbxAUFASpVIop\nU6bgxx9/dGpMxiz9vMw92xQeHu6wuLZv347c3Fx8/PHH/D5nxVRaWory8nKkpKQgKioKlZWVGDFi\nBGpra536fVIoFJgyZQoAYOTIkRCLxairq3P6z66goAAPPfQQgI6/wYKCAgDd//n1yaThKs9xMMaw\nYMECJCQk4Omnn+b3Z2VlYceOHQCAHTt28MnEEVavXo2KigqUlZVh9+7duPvuu7Fz506nxhQaGoqB\nAwfyCz0dPHgQiYmJ+N3vfue0mOLi4pCfnw+1Wg3GGA4ePIiEhASnxmTM0s8rKysLu3fvhlarRVlZ\nGUpKSjBq1CiHxJSXl4f169dj79698PDwEMTqjJiSk5NRW1uLsrIylJWVQaFQoLCwECEhIU79Pk2e\nPBmHDx8GAFy8eBFarRb9+vVzakwAEB0djaNHjwLoqMQ9dOhQAD34+dlh/MUhcnNz2dChQ9mQIUPY\n6tWrnRLDd999x0QiEUtJSWGpqaksNTWVffPNN6y+vp7dc889Tp1yyxhjR44c4WdPOTumU6dOsfT0\ndMF0TWfHtG7dOn7K7WOPPca0Wq1TYpoxYwYLCwtjMpmMKRQK9v7773cZx6pVq9iQIUNYbGwsPxvG\n3jG99957LDo6mg0aNIj/XV+4cKFTYnJzc+O/T8aioqIEU26dFZNWq2WPPvooS0pKYsOHD2f/+c9/\nHBqTcVzGv1MnT55ko0aNYikpKWzMmDGssLCwR3GJGKNaHIQQQqzTJ7unCCGEOAclDUIIIVajpEEI\nIcRqlDQIIYRYjZIGIVaQSCRIS0vDsGHDMGXKFCiVSgAdlQDEYjFeeukl/ty6ujrIZDI89dRTAIAL\nFy4gIyMDaWlpSEhIQHZ2NoCOtaP9/PyQlpbG/+OmahLiqihpEGIFLy8vFBUV4cyZM/D19eVL/ANA\nVFQUcnNz+e3PPvsMSUlJ/FO1S5YswdKlS1FUVITi4mI+mQDAXXfdhaKiIv7f3Xff7bgvipAeoKRB\nSDeNHTsWpaWl/LaXlxfi4+Px888/AwA+/fRTTJ8+nS9pUVNTI3jCNikpif+YZryTvoaSBiHdYDAY\nsH//fsELPwDMmDEDu3fvRmVlJSQSCQYMGMAfe+aZZ3D33Xdj0qRJeOONN9DU1MQf++677wTdU2Vl\nZQ77WgjpCUoahFhBrVYjLS0NYWFhqKiowJNPPik4npmZiQMHDmD37t145JFHBMfmzZuH8+fP4+GH\nH8aRI0cwZswYaLVaAMD/+3//T9A9FRUV5bCviZCeoKRBiBU8PT1RVFSEy5cvw8PDA3v37hUcl8lk\nGDFiBDZu3IiHH37YpNspLCwM8+fPx5dffgmpVIpz5845MnxCbIaSBiHd4Onpic2bN+PFF180SQxL\nly7FunXr4O/vL9ifl5cHnU4HoGN8o76+3qHVTQmxJUoahFjBeH2B1NRUREdH49NPPxWsdZGQkIA5\nc+bw53P7Dxw4gOTkZKSmpuK+++7Dhg0b+FLnncc0vvjiCwd/ZYR0DxUsJIQQYjVqaRBCCLEaJQ1C\nCCFWo6RBCCHEapQ0CCGEWI2SBiGEEKtR0iCEEGK1/w8vg5XP9EXWqwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x100741ad0>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = plt.subplots(1)\n",
      "\n",
      "rmse1k = rmse1k[np.isfinite(rmse1k)]\n",
      "rmse2k = rmse2k[np.isfinite(rmse2k)]\n",
      "rmse4k = rmse4k[np.isfinite(rmse4k)]\n",
      "\n",
      "m = [np.median(x) for x in [rmse1k, rmse2k, rmse4k]]\n",
      "e_up = [np.abs(np.median(x)-stats.scoreatpercentile(x,0.95)) for x in [rmse1k, rmse2k, rmse4k]]\n",
      "e_down = [np.abs(np.median(x)-stats.scoreatpercentile(x,0.05)) for x in [rmse1k, rmse2k, rmse4k]]\n",
      "            \n",
      "ax.errorbar([1000,2000,4000], m, yerr=np.vstack([e_up, e_down]), \n",
      "            fmt = 'o', color='k')\n",
      "\n",
      "ax.set_ylabel('RMSE')\n",
      "ax.set_xlabel('b value')\n",
      "ax.set_xlim([900, 4500])\n",
      "ax.set_ylim([10,100])\n",
      "ax.set_xticks([1000,2000,4000])\n",
      "ax.loglog()\n",
      "fig.savefig('figures/rmse_bars.svg')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEUCAYAAAAiMOHqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD5JJREFUeJzt3X9oVfUfx/HX/W6SQZmWZHo3mjlxW5tZaaYwnEQskQxk\nvySRXIaZ2yQqoj/EmQ0yCGqaFCgKoraNrGWry9J1wYgcWvmjaf5oq2GFf0iauR867/ePuhev27FZ\n95x7797PB1zw3uu29x16n/ece+7n+EKhUEgAAHP+F+8BAADxQQAAwCgCAABGEQAAMIoAAIBRBAAA\njCIAAGAUAQAAo1LjPYAkNTY2qqmpSefPn9fTTz+tRx99NN4jAcCQ50ukTwL//vvvevHFF7Vx48Z4\njwIAQ55ru4DKy8s1ZswY5eXlRd0eCASUlZWliRMnau3atVH3vfbaa6qoqHBrJADAVVwLwOLFixUI\nBKJu6+vrU0VFhQKBgNra2rRjxw4dPXpUoVBIL7/8subMmaMpU6a4NRIA4CquvQeQn5+vjo6OqNta\nW1uVmZmpjIwMSVJZWZkaGxu1e/du7dmzR+fPn9fJkye1dOlSt8YCAPzN0zeBT58+rfT09Mj1tLQ0\n7du3T+vWrVNlZeU/fr3P53NzPAAYsgZ6u9fTw0Bj8QQeCoViflm1apUr35eLjQv/fvh9JvrjdeJp\nAPx+vzo7OyPXOzs7lZaW5uUIQMwVFBTEe4QhxdrvM56P19MATJ06VSdOnFBHR4d6e3tVV1enefPm\neTkCEHPWnrDcZu33OSQDsGDBAs2cOVPHjx9Xenq6Nm/erNTUVK1fv16FhYXKyclRaWmpsrOz3Rph\n0Kz9gwMAKcE+CPZPfD7fdfdnAQD6c3ruZC0gADCKAACAUQQAAIwiAABgVNIFoLq6WsFgMN5jAEDC\nCwaDqq6udryfo4AAYIjjKCAAQBQCAABGEQAAMIoAAIBRBAAAjCIAAGAUAQAAowgAABhFAADAqKQL\nAEtBAMDgsBQEABjHUhAAgCgEAACMIgAAYBQBAACjCAAAGEUAAMAoAgAARhEAADCKAACAUQQAAIwi\nAABgVNIFgMXgAGBwWAwOAIxjMTgAQBQCAABGEQAAMIoAAIBRBAAAjCIAAGAUAQAAowgAABhFAADA\nKAIAAEYRAAAwigAAgFEEAACMSroAsBw0AAwOy0EDgHEsBw0AiEIAAMAoAgAARhEAADCKAACAUQQA\nAIwiAABgFAEAAKMIAAAYRQAAwCgCAABGEQAAMIoAAIBRBAAAjCIAAGAUAQAAo5IuAJwRDAAGhzOC\nAYBxnBEMABCFAACAUQQAAIwiAABgFAEAAKMIAAAYRQAAwCgCAABGmQ/ARx99pD/++CPeYwCA51Lj\nPUC8BINBBYNB1dbWauHChbr99tslSQUFBSooKIjvcADgAfNLQWRmZioQCCgzMzOm3xcAEgVLQQAu\n+Oyzz3Ty5Ml4jwH8K2Z3AQH/RXgXYn19vXJzc5WTkyOJXYhILgQA+BfCT/RHjhxRSUmJioqK4j0S\nkkj4BcTx48fV19en7OxsSd6/gCAAAOCx8BN9TU2NLl68eN01+93EewAAYBQBAACjCAAAGEUAAMAo\nAgAARiVdAKqrqxUMBuM9BgAkvGAweN0jjJLuMNB4HS4FAMkmfLjp6tWrB7w/6bYAAACxQQAAwCgC\nAABGEQAAMIoAAIBRBAAAjCIAAGAUAQAAoxwD0NLSEvlze3t71H07d+50byIAgCccA/DCCy9E/jx/\n/vyo+9asWePeRAAAT7ALCACMMh2ApqYm/fLLL3ryySdVWFiopqameI8EAJ5xXAzuxx9/1Lx58xQK\nhdTe3q7HH388ct+17wkko6amJq1YsUJdXV1qbW2VJJ06dUqSNHfu3HiOBgCe8IVCodBAd/zTkste\nnrk+zOfzyWHcG1ZYWKjm5uYBbw8EAjH5GRjampqaVF5ertGjRystLU1VVVW8eMANCZ8UvqamxtWf\n4/Tc6bgFcO0TfG9vr77//nv5/X7deeedMR/Qaz09PQPe3t3d7fEkSEbhLcgzZ87ozJkzamtrYwsS\nScfxPYClS5fqyJEjkqRz587pvvvu06JFizRlyhRt377dswHdctNNNw14+/Dhwz2eBMmotrY28oQf\ndurUKa1bty5OEwE3zjEAe/fuVW5uriRp8+bNmjRpkg4fPqxvvvlGb7zxhmcDuqWqqkoTJkyIum3C\nhAmqrKyM00RIJmxBYihw3AV09Svk5uZmFRcXS5Luuusu96fyQHgzvbi4WHl5eRo1apQqKyvZfMeg\nsAWJocAxALfddpt27dolv9+vr776Sps2bZIkXbp0aci8ypk7d67GjRunbdu2KTMzM97jIIlUVVXp\n1KlTUbuB2IJEsnEMwHvvvaeqqir99ttveuuttzR27FhJ0p49e3iVDPPC/wfCRwGlp6ezBYkb0tTU\npC1btqi3t1f79++Py1FkjoeBJqJYHgYalpmZqUAgwBYA/pWioiKVlZWpqKgo3qMgiYSPIrt2C/Lt\nt992JQI3fBhoZWWl4xf5fD7V1tbGdkIAMOJ6R5F5uRXgGIB3331Xubm5Kikp0bhx4yQpEgOfz+fN\ndAAwBCXKUWSOAfj111/V0NCg+vp6paSkqLS0VMXFxRo5cqSX8wHAkJMoR5E5fg5g9OjRWrZsmb74\n4gtt2bJF586dU05OjrZu3erlfAAw5CTK55ActwDCDhw4oPfff1+ff/655syZowcffNCLuRxVV1er\noKAgLmsRAUAshPfzP//88+rp6VF2drYrR5EFg8HrruvmeBTQypUr9emnnyo7O1tlZWUqLCzUsGHD\nYjrcjeIoICQajgLCf5Gwi8HV1NRo/PjxOnjwoA4ePKhXXnkl6psdOnTInUkBAJ647vkAnHAUEAAk\nP8cAZGRkDHh7KBRSfX297r77brdmAgB4wPEooAsXLujNN9/Uc889pw0bNujKlSv68MMPde+992rb\ntm1ezggAcIHjFsCiRYs0YsQIzZgxQ83NzdqyZYuGDx+u7du3a8qUKV7OCABwgWMATp48GXmjd8mS\nJRo7dqx++ukn3XzzzZ4NBwBwj+MuoJSUlKg/+/1+nvwBYAhx3AI4dOiQbr311sj1rq6uyHWfz6fz\n58+7Px0AwDWOAejr6/NyDgCAxxx3AQEAhjYCAABGEQAAMIoAAIBRBAAAjPrH8wEA6C+8znpWVpYO\nHTqkI0eOSBLnqkBScTwfQCKK5fkAwv+BL1++rJSUlMgKp/wHBuCVhD0fwFDHEz0A63gPAACMIgAA\nYBQBAACjCAAAGEUAAMAoAgAARhEAADDK7AfBACBewh9EvZZbn09yeu4kAAAwxDk9d7ILCACMIgAA\nYBQBAACjCAAAGEUAAMAoAgAARhEAADCKAACAUQQAAIwiAABgFAEAAKMIAAAYlXQBqK6uHnAVPQBA\ntGAwqOrqasf7WQ0UAIY4VgMFAEQhAABgFAEAAKMIAAAYRQAAwCgCAABGEQAAMIoAAIBRBAAAjCIA\nAGAUAQAAowgAABhFAADAKAIAAEYRAAAwigAAgFEEAACMIgAAYBQBAACjCAAAGEUAAMAoAgAARhEA\nADCKAACAUQQAAIwiAABgFAEAAKMIAAAYRQAAwCgCAABGEQAAMIoAAIBRBAAAjCIAAGAUAQAAowgA\nABhFAADAKAIAAEYRAAAwigAAgFEEAACMIgAAYBQBAACjCAAAGEUAAMAoAgAARhEAADCKAACAUQQA\nAIwiAABgFAEAAKMIAAAYRQAAwCgCAABGEQAAMIoAAIBRBAAAjCIAAGAUAQAAowgAABiVEAFob2/X\nkiVLVFxcHO9RAMCMhAjA+PHjtXHjxrj9/GAwGLefDQDx4loAysvLNWbMGOXl5UXdHggElJWVpYkT\nJ2rt2rVu/fgbQgAAWORaABYvXqxAIBB1W19fnyoqKhQIBNTW1qYdO3bo6NGjbo0AeIIXELFl7fcZ\nz8frWgDy8/M1atSoqNtaW1uVmZmpjIwMDRs2TGVlZWpsbNTZs2f17LPP6rvvvkuYrQJgsKw9YbnN\n2u8zro835KL29vZQbm5u5HpDQ0NoyZIlketbt24NVVRUDPr7SeLChQsXLv/iMpBUecjn8/2nr/+r\nAQCAWPD0KCC/36/Ozs7I9c7OTqWlpXk5AgDgb54GYOrUqTpx4oQ6OjrU29ururo6zZs3z8sRAAB/\ncy0ACxYs0MyZM3X8+HGlp6dr8+bNSk1N1fr161VYWKicnByVlpYqOzvbrREAANfhC7FjHQBMSohP\nAgMAvEcArnHs2DEtW7ZMJSUl2rRpU7zHAQDXsAvIwZUrV1RWVqb6+vp4jwIArjCxBXCj6xLt2rVL\nc+fOVVlZmdejAoBnTGwB7N27V7fccosWLVqkw4cPS/prXaJJkyZp9+7d8vv9mjZtmnbs2BF1VNIT\nTzyhxsbGeI0NAK7y9JPA8ZKfn6+Ojo6o265el0hSZF2iM2fOaOfOneru7tbs2bO9HxYAPGIiAAM5\nffq00tPTI9fT0tK0b98+zZo1S7NmzYrjZADgDRPvAQzkv65LBADJzmwAWJcIgHVmA8C6RACsMxEA\n1iUCgP5MHAYKAOjPxBYAAKA/AgAARhEAADCKAACAUQQAAIwiAABgFAEAAKMIAEzr6Ojod56IWHjq\nqaf0wQcfxPz7ArFEAAAX+Hw+FhxEwiMAMO/y5ctauHChcnJyVFxcrK6urqj7jx07punTp0eud3R0\naPLkyZKkV199VQ899JDy8vK0dOnSAb9/RkaGzp49K0nav39/5DwTf/75p8rLyzV9+nQ98MAD+vjj\nj914eIAjAgDzfvjhBy1fvlxtbW0aMWKENmzYEHV/VlaWent7IycVqquri5wutLKyUq2trTp8+LC6\nurr0ySef9Pv+TlsCNTU1euSRR7Rv3z61tLTopZde0sWLF2P74IDrIAAwLz09XTNmzJAkLVy4UF9+\n+WW/v1NSUqK6ujpJUn19vUpLSyVJLS0tevjhhzV58mS1tLSora1t0D+3ublZr7/+uu6//37Nnj1b\nPT09UUuUA24ze0YwIOzqV+ihUGjAV+ylpaUqLi7W/Pnz5fP5NGHCBHV3d2v58uU6cOCA/H6/Vq9e\nre7u7n5fm5qaqitXrkhSv/t37typiRMnxvgRAYPDFgDM+/nnn/X1119LkrZv3678/Px+f+eee+5R\nSkqK1qxZE9n9E34yv+OOO3ThwgU1NDQM+P0zMjK0f/9+SYo6MqiwsFC1tbWR699++21sHhAwSAQA\npvl8Pk2aNEnvvPOOcnJydO7cOS1btmzAv1taWqpt27appKREkjRy5Eg988wzys3N1WOPPRb1RvHV\nVq1apRUrVmjatGlKTU2NbGGsXLlSly5d0uTJk5Wbm6tVq1a58yABB5wPAACMYgsAAIwiAABgFAEA\nAKMIAAAYRQAAwKj/A1b4eXvspXnmAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x108d8f610>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def calc_snr(model):\n",
      "    b0 = model.data[..., model.b0_idx][model.mask]\n",
      "    b_data = model.data[..., model.b_idx][model.mask]\n",
      "    s = np.mean(b_data, -1)\n",
      "    sigma = np.std(b0, -1)\n",
      "    # Correct for small sample (according to http://nbviewer.ipython.org/4287207)\n",
      "    nb0 = len(model.b0_idx)\n",
      "    bias=sigma*(1-np.sqrt(2/(nb0-1))*(gamma(nb0/2)/gamma((nb0-1)/2)))\n",
      "    n = sigma + bias\n",
      "    return s/n\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "snr1k = calc_snr(TM_1k_1)\n",
      "snr2k = calc_snr(TM_2k_1)\n",
      "snr4k = calc_snr(TM_4k_1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = viz.probability_hist(snr1k, label='b=1000', color=[0.8, 0.8, 0.8])\n",
      "fig = viz.probability_hist(snr2k, fig=fig, label='b=2000', color=[0.59, 0.59, 0.59])\n",
      "fig = viz.probability_hist(snr4k, fig=fig, label='b=4000', color=[0.32, 0.32, 0.32])\n",
      "ax = fig.axes[0]\n",
      "ax.set_xlabel('SNR')\n",
      "ax.set_ylabel('P(SNR)')\n",
      "plt.legend()\n",
      "ax.text(20, 0.20, 'median(1000): %2.2f \\nmedian(2000): %2.2f \\nmedian(4000): %2.2f'%(np.median(snr1k), np.median(snr2k), np.median(snr4k)))\n",
      "fig.savefig('figures/snr_distributions.svg')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8U1X6P/DPTZo2zdJ9X2gp6d5SyiqDYisyWBkrCo4s\ngxsiOgKOM1XHlzjCKCjuCsPId8QRBbE/BS0M0EHBDotCFVtZWqClLd23dG/TZun9/dHJnYYmbZIm\nzcLzfr36ssndHq6QJ+ece57DsCzLghBCyA2PZ+sACCGE2AdKCIQQQgBQQiCEEPJflBAIIYQAoIRA\nCCHkvyghEEIIAWDlhJCbm4u4uDhER0dj8+bNQ7bn5eXB09MTqampSE1NxSuvvGLNcAghhAzDxVon\n1mg0WL16Nb799luEhoZi2rRpyMzMRHx8vM5+t956K/bv32+tMAghhBjJai2E/Px8yGQyREZGQiAQ\nYPHixcjJyRmyH82LI4QQ+2C1hFBTU4Pw8HDudVhYGGpqanT2YRgG33//PVJSUnDnnXeiqKjIWuEQ\nQggZgdW6jBiGGXGfyZMno6qqCiKRCIcPH8aCBQtw5coVs85FCCFkKFN6YazWQggNDUVVVRX3uqqq\nCmFhYTr7SKVSiEQiAEBGRgZUKhVaWlr0no9lWbv6eemll2wegyPEZK9xUUwU040Ql6mslhCmTp2K\nkpISVFRUQKlUIjs7G5mZmTr7NDQ0cEHn5+eDZVn4+PhYKyRCCCHDsFqXkYuLC7Zu3Yp58+ZBo9Fg\nxYoViI+Px/bt2wEAq1atwpdffom///3vcHFxgUgkwueff26tcAghhIzAagkBGOgGysjI0Hlv1apV\n3O9PPvkknnzySWuGYDVpaWm2DmEIe4wJsM+4KCbjUEzGs9e4TMGw5nQ0jTGGYUzuD5PL5Th48CDS\n0tIwbtw4K0VGCCH2y9TPTqu2EGwpJycH3333HRiGwfLly20dDiE3JB8fH7S2tto6DKfn7e1t8IEc\nUzhlQlCpVMjLy8P8+fPpLyMhNtTa2mrW0y7ENJZ6NN8pi9uVlZUhMDAQ48aNQ1dXl63DIYQQh+CU\nCeHy5cuIjY2FRCJBd3e3rcMhhBCH4PQJgVoIhBBiHKdMCA0NDQgNDYVEIkFnZ6etwyGE2KHIyEgc\nPXrU1mHYFadMCJ2dnZBKpRCLxdRlRAjRi2GYUQ3GqlQqLFq0COPHjwePx8N//vOfIfs899xz8PPz\ng5+fH/785z/rbKuoqEB6ejrEYjHi4+OHJKfPPvsMERERkEgkuOeee8bkARmnTQgSiQQikQh9fX1Q\nq9W2DokQ4oRmz56NXbt2ISgoaEhy2b59O3JycnDu3DmcO3cOBw4c4Co1AMCSJUswZcoUtLS0YOPG\njVi0aBGam5sBABcvXsTjjz+O3bt3o6GhASKRCL///e+t/udxuoSgVCqh0WggFArBMAy1EgghBuXn\n5yMxMRE+Pj545JFH0NfXZ/SxAoEAa9euxaxZs8Dn84ds37lzJ7KyshASEoKQkBBkZWXh448/BgBc\nuXIFBQUF2LBhA9zc3HDvvfdi4sSJ2Lt3LwBg9+7dyMzMxM033wyxWIyXX34Z+/bts/pnmdMlBG13\nkTZb08AyIUQflmXx2Wef4ciRI7h69SquXLmCV155BVVVVfDy8oK3t7feH2NrrhUVFSElJYV7PXHi\nRFy8eBHAQAsgKioKYrGY256SkqKzffCxUVFRcHNz07s8gCU53cQ0bULQooFlQuzb2bNnLXKeKVOm\nmLQ/wzBYvXo1QkNDAQAvvPAC1qxZg5dffhltbW2jjqerqwuenp7caw8PD+7L6fXbgIHlAOrq6gAA\n3d3dQ7Z7eHhY/bPM6RNCcHAwqqqqEBcXZ8OoCCGGmPpBbkmDV3UcN24camtrLXZuiUSCjo4O7nV7\nezskEonebdrt2s8uiUSC9vZ2g9utxWm7jLSSkpJw4cIFG0ZECLFXlZWVOr9rF/aSSCSQSqV6f/bs\n2WPUuRMTE1FYWMi9/uWXX5CUlMRtKysr0+nO/uWXX5CYmMht/+WXX7htV69ehVKpRExMzKj+vCNx\nyoSgzcIAkJycjPPnz1M9FUKIDpZl8be//Q01NTXckz73338/wsPD0dXVhc7OTr0/S5Ys4c7R19eH\n3t7eIb8DwAMPPIC3334btbW1qKmpwdtvv42HHnoIABATE4NJkyZhw4YN6O3txb59+3DhwgUsXLgQ\nALBs2TIcOHAAJ0+eRHd3N1588UUsXLhQZ8zBGpy+yyggIAC9vb3o7e2Fu7u7DSMjhNgThmGwbNky\n/PrXv0ZtbS0WLFiAdevWmXSO2NhYVFZWgmEYzJs3DwzDoLy8HOPGjcOqVatQVlaG5ORkAMDKlSvx\n2GOPccd+/vnneOihh+Dj44OIiAjs3bsXvr6+AICEhAR88MEHWLZsGeRyOebOnYt//vOflvvDG+B0\n6yF88sknkEqluOeee7j3Vq5ciVdffRV+fn7WCpEQooc5a5kQ0xm6z6bef6frMurp6YFIJNJ5j+Yi\nEELIyJwuISgUiiFdQ5QQCCFkZDdEQqAy2IQQMjKnSwg9PT16Wwg0W5kQQobndAmht7eXxhAIIcQM\nTpcQenp6IBQKdd6jhEAIISNzuoSgUCiohUAIIWZwyoRAg8qEEGI6p0oIKpUKGo0Grq6uOu+LRCIa\nVCaE6KAlNIdyqoSgLU9x/cpF1EIghFxvtEtonj59GnPnzoWvry8CAgLw29/+FvX19Tr70BKaNqSv\nuwgAhEKhTtEpQggZrba2Njz++OO4du0arl27BqlUiocffpjb7ohLaIJ1AMaGWV5ezv7hD38Y8n51\ndTX75JNPWjosQsgI7PkjJjIykn311VfZhIQE1tvbm3344YfZ3t5es8939uxZViqVcq9nzpzJ/uMf\n/+Bef/TRR+xNN93EsizLXr58mXVzc2O7urq47bNnz2Y/+OADlmVZ9vnnn2eXLVvGbbt69Srr6uqq\ns/9ghu6zqfffqVoIhiqaCoVCk9ZKJYQ4P9bCS2geP36cW+8AoCU0bU7fLGWAuowIsWd///vfLXKe\nJ554wqT9LbmE5rlz5/Dyyy9j//793Hu0hKaNdXd36yyOo+Xm5obe3l6wLDuqQSRCiOWZ+kFuSZZY\nQrO0tBR33nkn3n//fcyaNYt7n5bQtLGurq4hk9IAwMXFBTweDyqVygZREULs1WiX0Lx27Rrmzp2L\nv/zlL1i2bJnOuR1xCU37HfEZxNgwv/jiC3bXrl16tz3wwANsR0eHJcMihIzAnj9iIiIi2OTkZLa6\nupqVy+XsrFmz2BdeeMHo46urq9moqCj2zTff1Lv9gw8+YOPj49mamhq2urqaTUhIYLdv385tv+mm\nm9isrCxWoVCwe/fuZb28vNjm5maWZVn24sWLrIeHB3vixAm2q6uLXbJkCbtkyRKDsRi6z6bef6dq\nIXR3dxtcc1TbbUQIIYDuEpoTJkxAdHS0SUtofvjhhygvL8f69eu51oOHhwe3fdWqVbjrrruQnJyM\niRMn4q677hqyhOZPP/0EHx8fvPDCCwaX0AwMDIRCocC2bdss94c3wKmW0Ny2bRuio6Mxd+7cIdvW\nrl2LZ555RqfPkBBiXbSE5thwiCU0c3NzERcXh+joaGzevNngfj/++CNcXFywb9++UV3P0KAyQE8a\nEULISKyWEDQaDVavXo3c3FwUFRVhz549KC4u1rvfc889hzvuuGPU3ySG6zKiuQiEEDI8qyWE/Px8\nyGQyREZGQiAQYPHixcjJyRmy35YtW7Bo0SL4+/uP+ppdXV00hkAIIWay2jyEmpoanf76sLAwnDlz\nZsg+OTk5OHbsGH788cdh5wisX7+e+z0tLQ1paWlD9qEuI0LIjSwvLw95eXlmH2+1hGDMBLA//OEP\neO2117iBj+G6jAYnBEO6u7v1zkMAqMuIEOL8rv+yvGHDBpOOt1pC0E7w0KqqqkJYWJjOPmfPnsXi\nxYsBAM3NzTh8+DAEAgEyMzNNvl5/f7/e1dK0hEIhFAqFyeclhJAbhdUSwtSpU1FSUoKKigqEhIQg\nOztbZ4YfAJSVlXG/P/zww7jrrrvMSgYAoFQqIRAIwOfz9W53c3OjFgIhhAzDagnBxcUFW7duxbx5\n86DRaLBixQrEx8dz9cBXrVpl0eupVKohK6UN5u7uTi0EQggZhlWL22VkZCAjI0PnPUOJ4J///Oeo\nrqVtIRji5eVlVuEqQohzioyMxI4dOzBnzhxbh2I3nKZ0hUqlGjYh+Pj4jMkSdIQQxzDaJTQH++tf\n/woej4djx47pvE9LaNqIUqkctsvI29sbLS0tYxgRIeRGcPXqVXz55ZcICQnRed8Rl9B0moRALQRC\niKny8/ORmJgIHx8fPPLII2Y9eLJ69Wps3rx5yOfPzp07kZWVhZCQEISEhCArKwsff/wxAODKlSso\nKCjAhg0b4ObmhnvvvRcTJ07E3r17AQC7d+9GZmYmbr75ZojFYrz88svYt28furu7R/1nHs4NkxCk\nUil6e3uhVCrHMCpCiL1iLbCE5hdffAGhUDhkrBSgJTRtaqSEwDAMvL290draisDAwDGMjBAynIUL\nF1rkPNpv18Ya7RKanZ2deOGFF/Dtt9/q3U5LaNrQSGMIwEC3UUtLCyUEQuyIqR/kljSaJTTXr1+P\n5cuXY9y4cdx7g6st0BKaNjRSCwEYuMmDl6wjhNzYRrOE5rFjx/D+++8jODgYwcHBqKqqwm9/+1u8\n8cYbABxzCU2naSEYkxBotjIhRItlWfztb3/Db37zG7i7u2Pjxo24//77ER4ebtQXx6NHj0KtVnPn\nmjZtGt555x1uPOGBBx7A22+/jTvvvBMsy+Ltt9/GU089BQCIiYnBpEmTsGHDBrz88ss4dOgQLly4\nwHWfLVu2DDNnzsTJkyeRmpqKF198EQsXLjRYzdlSbqiEQAXuCCFag5fQrK2txYIFC0xaQtPHx0fn\nNZ/Ph7e3N1dPbdWqVSgrK0NycjIAYOXKlUOW0HzooYfg4+ODiIgIg0toyuVyzJ07d9STd43hNEto\n/vvf/0ZFRcWwJTF27NiBoKAgzJ8/39IhEkL0oCU0x4ZDLKE5lowZVKZFcgghxDCnSQjUZUQIIaPj\nNAlhpOJ2ALUQCCFkOE4zqKxWqw0ujqNFLQRCxpa3t7fFCsgRw7y9vS1yHqdJCEqlEl5eXsPuQy0E\nQsYWFZR0LE7TZUTzEAghZHScJiEY85SRUCikFgIhhBjgNAlBpVLBxWX4HjBqIRBCiGFOlRCMaSFQ\nQiCEEP2cKiHQY6eEEGI+p0kIxo4hUAuBEEL0c5qEYOwYArUQCCFEP6dJCAqFAkKhcNh93NzcoFQq\nqdgWIYTo4TQJoaenZ8Ra4Xw+H3w+n9ZVJoQQPZwqIYxUugIA3N3doVAoxiAiQghxLE6REDQaDfr6\n+uDu7j7iviKRiBICIYTo4RQJQTt+wOON/MehFgIhhOjnFAnB2O4iYCAh9PT0WDkiQghxPE6TEIxd\nfJq6jAghRD+nSQjGthCEQiElBEII0cMpEkJ3d7fRCYFaCIQQop9TJAQaQyCEkNFzmoRAYwiEEDI6\nTpMQaAyBEEJG54ZLCCKRiLqMCCFED6smhNzcXMTFxSE6OhqbN28esj0nJwcpKSlITU3FlClTcOzY\nMbOu09vba9QsZYAmphFCiCHD14seBY1Gg9WrV+Pbb79FaGgopk2bhszMTMTHx3P73H777bj77rsB\nAOfPn8c999yD0tJSk6/V29sLNzc3o/alFgIhhOhntRZCfn4+ZDIZIiMjIRAIsHjxYuTk5OjsM3gg\nuKurC35+fmZdq7e3d8TS11rUQiCEEP2slhBqamoQHh7OvQ4LC0NNTc2Q/b7++mvEx8cjIyMD77//\nvlnXUigURncZSSQSdHV1mXUdQghxZlbrMmIYxqj9FixYgAULFuDEiRNYvnw5Ll++rHe/9evXc7+n\npaUhLS2Ne21KC0EqlVJCIIQ4pby8POTl5Zl9vNUSQmhoKKqqqrjXVVVVCAsLM7j/LbfcArVaDblc\nDl9f3yHbByeE65mSELRjCBqNBnw+36hjCCHEEVz/ZXnDhg0mHW+1LqOpU6eipKQEFRUVUCqVyM7O\nRmZmps4+V69e5Zaz/PnnnwFAbzIYiSkJgc/n08AyIYToYbUWgouLC7Zu3Yp58+ZBo9FgxYoViI+P\nx/bt2wEAq1atwt69e/HJJ59AIBBAIpHg888/N+tapjxlBAx0G3V2dkIqlZp1PUIIcUYM6wArzjMM\ng+HCXLp0KXbs2GH0wPKf//xnPPLII4iJibFUiIQQYndG+uy8nsPPVO7v74dKpTKrhUAIIeR/HD4h\n9PX1wdXV1ajlM7Xo0VNCCBlqxDGEixcv4vjx46ioqADDMIiMjMQtt9yCxMTEsYhvRKYMKGtRC4EQ\nQoYy+LX6008/xfTp05GVlYX6+npERUUhMjISdXV1yMrKwrRp07Br166xjFUvcxKCRCKhhEAIIdcx\n2EJobW3F0aNHDT6J09HRgY8//thacRnNnIQgFotRV1dnpYgIIcQxGUwIa9euNXhQT08PPDw8ht1n\nrJj6yClAi+QQQog+w47ENjY24scff0RLSwuAgQJ0mzZtQmxs7JgEZwxzWgg0MY0QQoYymBC2bduG\n5ORkrF27FvHx8diyZQtSUlLQ2NiI/Pz8sYxxWH19fSa3EKjiKSGEDGUwIWzduhXFxcX44YcfkJ+f\nj6ysLPy///f/8O677yI4OHgsYxyWUqmEq6urScdQC4FYWlpaGld+Zf78+ejo6DD7XFu3buXG5774\n4gskJiaCz+dz59d69dVXER0djbi4OBw5coR7/+zZs0hOTkZ0dDSeeuop7v2+vj7cf//9iI6Oxk03\n3YRr164BABoaGnDnnXeOGFdnZydSU1O5H39/fzz99NMG96+srIREIsFbb70FYKCref78+YiPj0dS\nUhKef/55o+8JGRsGE4Kbmxt8fHwAABEREYiLi8OUKVPGLDBjmZMQ3N3dKSEQixpc3ffgwYPw8PAw\n6zwsy2LHjh343e9+BwBITk7GV199hdmzZ+vsV1RUhOzsbBQVFSE3Nxe///3vuRmpTzzxBHbs2IGS\nkhKUlJQgNzcXALBjxw74+vqipKQETz/9NJ577jkAQGBgILy9vYcknOtJpVIUFBRwPxEREVi4cKHB\n/f/4xz9i/vz53GuGYfDss8+iuLgYBQUFOHXqFBcbsQ8GE0J1dTXWrl2LNWvWYM2aNairq+Ne28Ng\nspZKpTKrhUBdRje2iooKxMXF4eGHH0ZsbCyWLVuGI0eOYNasWYiJicGPP/4IAOju7sYjjzyCGTNm\nYPLkydi/fz+AgTU4Fi9ejISEBNx77706f58iIyO5cbd77rkHU6dORVJSEv7xj39w+0gkEqxbtw6T\nJk3CzJkz0djYCAA4deoU4uLi4OIy8LxHXFyc3hIrOTk5WLJkCQQCASIjIyGTyXDmzBnU1dWhs7MT\n06dPBwA88MAD+PrrrwEA+/fvx4MPPggAWLhwIY4ePcqdLzMzE3v27DH6/l25cgWNjY24+eab9W7/\n+uuvERUVhYSEBO49d3d33HrrrQAAgUCAyZMn610jhdiOwYTwxhtvYMqUKdzP66+/jsmTJ3Ov7YV2\nprIpqMuIAAPVdrOysnDp0iVcvnwZ2dnZOHXqFN58801s2rQJALBx40bMmTMHZ86cwbFjx/DMM8+g\np6cHf//73yGRSFBUVIQNGzbg7Nmz3HkHtxY++ugj/PTTT/jxxx/x/vvvo7W1FcBA98nMmTNRWFiI\n2bNnc8ni5MmTmDZt2oix19bW6pST1y5Adf37oaGh3Ifu4EWrXFxc4OnpySWu6dOn4/jx49y5B3+z\n1+fzzz/H4sWL9W7r6urC66+/PmzJ+ra2Nhw4cABz5swZ8c9Kxo7Bx04feuihMQzDfOZ0Gbm6ukKt\nVkOlUkEgEFgpMmLvxo8fz824T0xMxO233w4ASEpKQkVFBQDgyJEjOHDgAN58800AA19AKisrceLE\nCa5/Pjk5GRMnTtR7jffee4/7hl5VVYWSkhJMnz4drq6u3IfulClT8M033wAY6Hc39K3bmoKDg7k/\nc0hICA4ePDjs/tnZ2QYnpq5fvx5PP/00RCKR3sJqarUaS5YswVNPPYXIyMjRhk4syGBCePjhh/W+\nr/3289FHH1knIhMplUqjq5xqMQwDkUiE3t5eSgg3sMFPp/F4PO6LBY/Hg1qt5rbt27cP0dHRQ44f\nqYpkXl4ejh49itOnT0MoFCI9PR29vb0AoPP37vrrGVOd8voFqKqrqxEWFobQ0FBUV1cPeV97TGVl\nJUJCQqBWq9He3s6NE7Isa/Qqh7/88gvUajVSU1P1bs/Pz8fevXvx7LPPoq2tDTweD+7u7vj9738P\nAHjssccQGxtrV13PZIDBhDB//nyd0qkMw6Cqqgpvv/02NBrNmAU4EpVKZdYAnrbbiNZEIMOZN28e\n3n//fWzZsgUAUFBQgNTUVMyePRufffYZ0tPTceHCBZw7d27IsR0dHfD29oZQKMSlS5dw+vTpEa8X\nERGB+vp6vdsGJ4rMzEwsXboUf/zjH1FTU8O1PBiGgYeHB86cOYPp06fj008/5T54MzMzsXPnTtx0\n00348ssvdbpr6urqEBERYdQ92bNnD5YuXWpwu7brCRhYsUsqlXLJYN26dejo6MCOHTuMuhYZWwbH\nEBYtWoSFCxdi0aJFSE1NxeHDh7Ft2zY8//zzKC8vH8sYh6VUKk2ehwDQk0Zk6Lrfg19rf3/xxReh\nUqkwceJEJCUl4aWXXgIw8CRPV1cXEhIS8NJLL2Hq1KlDzn/HHXdArVYjISEBzz//PGbOnGnwWtrX\nN998M3766Sdu21dffYXw8HCcPn0a8+fPR0ZGBgAgISEBv/3tb5GQkICMjAxs27aNO8e2bdvw6KOP\nIjo6GjKZDHfccQcAYMWKFZDL5YiOjsa7776L1157jbtOfn4+9yTTSGMIX3zxBZYsWaLz3oEDB7h7\nY0h1dTU2bdqE4uJiTJ48GampqXbT00AGDLtATnFxMTZu3Iiff/4ZzzzzDJYvX849/TCWhlvkYcuW\nLUhKSkJ6erpJ51y3bh2WLl2q8xQEIbbGsiwmT56MM2fOmDw2NhrLli1DVlaWwW4g4pgstkDOokWL\nMH/+fMycORN5eXnIzMxER0cHWlpauCcT7IE5g8rAQAuhu7vbChERYj6GYbBy5Urs3r17zK7Z2NiI\ntrY2SgbEcAtBO/qvb6CJYRiUlZVZNbDrr2coy7366quYM2cO99y1sd555x1MmTJlyIQfQghxFqa2\nEAz2/2gfQbN35rYQaC4CIYToMthldO3aNbS1tXGvjx07hrVr1+Ltt9+GUqkck+CMQQmB2ANr1TJ6\n5plnEB8fj5SUFNx7771ob2/n9hvrWkaDZWZmIjk5We+2lpYWpKenQyqVYs2aNTrbDMVG7IPBhHDf\nffdxH5iFhYW47777EBERgcLCQu4RMntACYHYA2vVMvr1r3+Nixcv4pdffkFMTAxeffVVALapZaS1\nb98+SKVSg/MWhEIhXnnlFW4y32CGYiP2wWBC6O3tRUhICABg165dWLFiBf70pz/h448/xpkzZ8Ys\nwJGMJiFQPaMbl6PUMpo7dy54vIF/pjNmzOAmndmqllFXVxfeeecdrFu3zmDftEgkwqxZs4Y8Dj5c\nbMQ+GEwIg/9nHz16FLfddtvAAbxh19QZc9ZoIVRWVqK0tBT9/f2jDY/YMUerZfTRRx9xXTu2qmX0\n4osvIisrCyKRaMT7e30LoqamxmBsxD4Y/HRPT0/Hfffdh7Vr16KtrY1LCLW1tWZNBLMWSycEpVKJ\nY8eOobCwENnZ2dw3N+J8tLWMGIYZtpbRa6+9htTUVKSnp+vUMhpconq4WkbaVoC2lhGAIbWMtNer\nrKxEUFDQkPNs3LgRrq6uw84QHg1jahkVFhairKwMd999t0lPrhDHYfApo3fffRfZ2dmor6/HyZMn\nuQ/dhoYGbNy4ccwCHMlo5iHoSwi//PILwsPDcdttt6GsrAyHDh1CRkYGAgMDLREusSOOUsvo448/\nxqFDh3S6eGxRy+j06dP46aefMH78eKjVajQ2NuK2227DsWPHhj1ucMzXxxYaGmrUsWRsGGwhMAyD\nJUuW4Omnn9b5n5aamop58+YBMK4Il7VZsoXQ09OD8+fPY9q0aWAYBhMmTEB6ejoOHTrE/UMmNxZt\nLSOtgoICAOBqGQGwai2j3NxcvPHGG8jJydFZOzwzMxOff/45lEolysvLuVpGQUFBXC0jlmXx6aef\n4u677+aO2blzJwCYVcvo8ccfR01NDcrLy3Hy5EnExMQMmwyu/3wIDg4eEtuCBQtGvCdk7BhMCGlp\naXjjjTdw5cqVIdsuX76MzZs3c4td2Ep/fz/UarVZFUvFYjGXEFiWxaVLl/DVV18hISFB5wmRiIgI\n+Pn5GSw4RhyXI9QyWrNmDbq6ujB37lykpqZyT/jZspYRMLRFcX0to8jISO4hlPDwcFy6dGnY2Ih9\nMDhTua+vD7t378aePXtw4cIFSKVSsCyLrq4uJCUlYdmyZVi6dOmY1FsxNNuur68PDz30kEkrPWl1\ndHRgzZo12LlzJ0pKSnD27FnMnj2be7JqsPz8fAAweTY0IaaiWkbEkiw2U9nNzQ2PPPIIHnnkEWg0\nGjQ3NwMA/Pz8wOfzRx+pBZjbXQT877FTlmVRVVWFpKQkvckAAAICAnD+/PnRhEqIUQbXMjK0Joml\nUS0jomUwISgUCnzwwQcoLS3FxIkTsWLFCptUOh2OUqk0e4EbFxcX8Pl89PX1oaamBpMnTza4b2Bg\nII4ePWrSIiKEmGusJ34GBASMuEIauTEYHEN48MEHuWnmhw4dwp/+9KexjMsoo2khAANPGjU0NIBl\nWXh6eg67n1Ao1CnlQQghzsZgQiguLsauXbvw+OOPY+/evTqrINkLSySEa9euITQ0dMRv/oGBgWho\naDD7WsS5WauWkdZbb70FHo+nU3p+rGsZKRQKzJ8/H/Hx8UhKSsLzzz9vcF9DsaWlpSEuLg6pqalI\nTU3luqLdtvdMAAAgAElEQVSJfTCYEAZ3D9lbV5GWSqUaVUIQCoWoqqrSmT1pSEBAACUEYpC1ahkB\nQFVVFb755hudx0JtVcvo2WefRXFxMQoKCnDq1Cm9tYiGi41hGHz22WcoKChAQUEB/Pz8zLpPxDoM\nJoRz585BKpVyP+fPn+d+N/cvu6WNtoUgFApRV1dn1OQYaiE4F0epZQQAf/zjH/H666/rxG+LWkbu\n7u7co+YCgQCTJ0/WW3rCUGxa9jB/iehnMCFoNBp0dnZyP2q1mvt9NM1hSxptQuDz+eDxeJBIJCPu\n6+fnh/b2dqhUKrOvR+yLI9QyysnJQVhY2JDSGLaqZaTV1taGAwcO6ExuGy622tpa7vWDDz6I1NRU\nvPLKK8Neg4w9q1eqy83NRVxcHKKjo7F58+Yh23fv3o2UlBRMnDgRs2bN0jvj05DRPGUEDHQ5Gdtk\n5fP58Pb2hlwuN/t6xL7Yey2jnp4ebNq0CRs2bODOZ61v18bUMtJSq9VYsmQJnnrqKW5lRWPt3r0b\nFy5cwIkTJ3DixAl8+umno4iaWJpVBwc0Gg1Wr16Nb7/9FqGhoZg2bRoyMzMRHx/P7RMVFYXjx4/D\n09MTubm5eOyxx4ya4g+MroXQ39+P3t5eeHl5GX1MQEAAGhsb9RYfI47H3msZXb16FRUVFUhJSQEw\nUPtnypQpOHPmjE1qGWk99thjiI2Nxdq1a/Vu1xebtltWO9dHIpFg6dKlyM/Px/Lly426LrE+q7YQ\n8vPzIZPJEBkZCYFAgMWLFyMnJ0dnn5kzZ3KPfA6u926M0Qwq19bWQigUmjSvwN/fH01NTWZdjzgm\nW9YySk5ORkNDA8rLy1FeXo6wsDD8/PPPCAwMtEktIwBYt24dOjo68M477xjcx1Bsgye4qlQqHDhw\nwOCqa8Q2rJoQBvdZAv/r5zRkx44dJi3lN5oWQmlpKQIDA01aJEfbQiDOwRFqGRmKzxa1jKqrq7Fp\n0yYUFxdj8uTJSE1NxUcffQRAt5aRodh6e3txxx13ICUlBampqQgPD8fKlSv1/lmJbRisZWQJe/fu\nRW5uLjdgtmvXLpw5cwZbtmwZsu93332HJ598EqdOnYK3t7dukAyjUzgrLS0NaWlp2L9/P+RyuclT\n/Pv7+/Hxxx9DLBajr68PDzzwgNHH7dixAw8++OCY1pkhNw6qZURGIy8vD3l5edzrDRs2WKaWkSVc\n35do6Jn/c+fOYeXKlcjNzR2SDLTWr18/5D1zu4za29vh6uoKT09PbnKOMXg8Hvz8/NDU1ER13IlV\nUC0jMhraL8tagx9IMIZVE8LUqVNRUlKCiooKhISEIDs7e8izzpWVlbj33nuxa9cuyGQyk85vbpeR\nXC6Hn58f3N3dTV7nQDuOQAmBWAvVMiK2YtWE4OLigq1bt2LevHnQaDRYsWIF4uPjsX37dgDAqlWr\n8Ne//hWtra144oknAAw8faEtNz0SpVIJd3d3k+OSy+Xw9fVFf3+/SWMIwMA/Hu3jeYQQ4kysXpMi\nIyMDGRkZOu+tWrWK+/3DDz/Ehx9+aNa5lUqlWes7Nzc3IyEhAZ2dnSYnhKCgIHz//ffo7+8Hj2f1\naRyEEDJmHPoTbTRdRr6+vhAKhSYnBA8PD4jFYtTV1Zl8XUIIsWcOnRBUKpXJM5UVCgVUKhWkUimC\ng4NRW1sLjUZj0jkmTJiAq1evmnQMIYTYO4dOCOa0ELStA4ZhIJVKERAQgLKyMpPOMWHCBJSVlaG/\nv9+k4wghxJ7dsAlBKyEhARcvXjTpHJ6entRtRAhxOjdcQmhubtZJCOHh4VypAFNERUWZ3LIghBB7\ndsMlBO0cBC2RSISenh6Trx0SEkJlLAghTuWGSggqlUqnwiNgfkLw9fVFS0sLjSMQQpyGQyeE7u5u\niMVio/dvbm6Gt7e3zmpUYrEY3d3dJl/b1dUV7u7udrNYECGEjNYNlRAaGxsRGBio8565LQRgYBU1\nWiScEOIsHDYhaDQa9PX1QSgUGn1MY2MjAgICdN4Ti8VmJwRfX19aQY0Q4jQcNiH09PRAJBKZVD6i\noaFhSEIYTQuBEgIhxJk4fEIwlkKhQF9f35AlM4VCIZRKpc4ShsaiLiNCiDNx2IRgzvhBQECA3lWy\nRCKRyTWNAEAqlUKpVJpcQpsQQuzRDZcQ9DG324hhGOo2IoQ4jRsqIfj7++vdJhKJzHr0FKBuI0KI\n87hhEkJbW5vOhLTBRjOwHBwcjJqaGrOOJYQQe3JDJASNRoPu7m5IpVK9282dnAYM1EKqra01a1Ca\nEELsiUMnBGOfMmpvb4dEIgGfz9e73cPDw+wZx25ubvDz80Ntba1ZxxNCiL1w2ISgUCiMTghtbW1D\nHjcdzMvLC21tbWbHEh4ejsrKSrOPJ4QQe+CwCaGvr8/ownYjJQRPT89RJYSIiAiDCaG/vx8KhQLt\n7e1ob28Hy7JmX4cQQqzJZeRd7JMpy2e2t7cPqWE0mJeXF4qKisyOxdfXl6uk6unpCQDo7e1FTU0N\n2tvb4erqCjc3N6jVatTU1CAsLAweHh5mX48QQqzBoROCKS2E2NhYg9u9vLzQ3t5udiwMwyAsLAw1\nNTUQi8Wora1FS0sLgoKCMH78eK68BsuyaGtrQ2VlJdzc3BAVFWVwXIMQQsaawyYEpVJpdAvB2l1G\nABAQEICmpiaIxWJoNBokJSXplNkGBhKHt7c3vLy8cO3aNVy9ehUymcykekyEEGItDvtJZGwLobe3\nF/39/XB3dze4z2hbCADg7++PxsZGtLa2Ijw8fEgyGIxhGERERIDP56OiooLGFQghdsFhE4Kxq6Vp\nWwfX1zAaTCwWQ6lUoq+vz+x4fH190draCk9PT6NaLgzDYPz48VCpVKisrIRGozH72oQQYgkOmxCM\nHVQeqbsIGPhwDgwMRH19vdnx8Pl8uLm5mbSkJ4/Hg0wmg0qlwrlz51BcXIyamhpKDoQQm7ghEoL2\nyZ/hRERE4Nq1a2bHI5fL4enpaXLXE5/Ph0wmQ0pKCsLCwtDT0zOqxEQIIeZy2IRgbJdRe3v7iC0E\nYCAhVFRUmBULy7JoaGhAaGgompqazDoHj8eDVCrFuHHj0NTURKUwCCFjzmETgrEthMFzA4YzmhZC\ne3s7+Hw+wsLCzE4IWm5ubvD29kZDQ8OozkMIIaZy2IRgTAuBZVmjE0JoaCjq6urMiqWhoQGBgYHw\n8/NDa2vrqMcAgoKCqJVACBlzDpsQjGkh9PT0QCAQGNW15OnpaVaBu+7ubiiVSnh7e0MgEMDLywvl\n5eUmn2cwaiUQQmzBqROCsa0DYGBNBJVKBZVKZVIcDQ0NOktzzp49GydPnhz1RDdqJRBCxppDJgSN\nRgONRjPs5C/AtITAMAykUqlJrYS+vj50dHTAz8+Pey8wMBDTp0/Hv//9b/T39xt9ruu5ubnBx8cH\nVVVVZp+DEEJM4ZAJQa1WQyAQDDvZDDAtIQAD6yKY8thoY2Mj/Pz8htQjio+Ph0qlMnuNBa3Q0FB0\nd3fTms2EkDHhkAnBlEdOTUkIpowjaDQayOVyBAQEDNnGMAx8fHzQ0tJi9LX14fP5iIqKQnV1NXp7\ne0d1LkIIGYnDJgRLPnKqZUoLobm5GR4eHgYTk6+v76gTAjAwthEcHIzy8vJRdUERQshIrJoQcnNz\nERcXh+joaGzevHnI9kuXLmHmzJkQCoV46623jD6vMQPKpjxyqmVsC0E7mzg4ONjgPpZoIWj5+/vD\n1dXV7MdiCSHEGFZLCBqNBqtXr0Zubi6KioqwZ88eFBcX6+zj6+uLLVu2ICsry6RzG9NlZMojp1rG\nrK2sUqlw9epVhIeHD1tB1ZIJgWEYhIeHo6mpaVQF+AghZDhWSwj5+fmQyWSIjIyEQCDA4sWLkZOT\no7OPv78/pk6davS6BlqWfuRUa6R1Efr7+1FaWgo/Pz/4+PgMey4vLy90dHRYrFCdq6srAgMDUV1d\nbZHzEULI9ay2QE5NTQ3Cw8O512FhYThz5ozZ51u/fj33e1RU1Ijf/I0tajeYn58ffvjhB4Pbr127\nBqFQiKCgoBHPxefz4eHhgba2Nvj6+poUhyGBgYG4ePEiOjo6aAlOQsgQeXl5yMvLM/t4qyWEkR4J\nNdXghHDu3DmUlZUNu79cLjf5g9jPzw/Nzc16t6nVarS3tyM5OdnoP5u228hSCYHH4yEsLAxVVVVI\nSEiw+D0mhDi2tLQ0pKWlca83bNhg0vFW6zIKDQ3VmVRVVVWFsLAwi5zbmC6j5uZmnQljxvDz80NT\nU5PeFcza2toglUpNWgPZkuMIWl5eXhAIBGhsbLToeQkhxGoJYerUqSgpKUFFRQWUSiWys7ORmZmp\nd19Tl5Ac6bHT/v5+yOVykxOCu7s7BAIBurq6hmxra2uDt7e3Sefz8fGx+KQyhmEwbtw41NfXo7W1\n1aLnJoTc2KzWZeTi4oKtW7di3rx50Gg0WLFiBeLj47F9+3YAwKpVq1BfX49p06aho6MDPB4P7733\nHoqKiiCRSIY990hPGXV0dEAoFMLNzc3kuLWtBKlUyr2n0WjQ2dmJ8ePHm3QuX19fNDc3g2VZi3bv\nCIVCREdHo6SkBAzDGLXeAyGEjMRqCQEAMjIykJGRofPeqlWruN+DgoLMqtUzUpdRU1OTya0DLe04\nQlRUFPdee3u7yd1FwMBjrO7u7qiurtYZYLcEkUgEmUyG0tJSMAxj8gA6IYRcz2FnKg/XQjBn/EAr\nNDR0yKOdra2tZn0LZxgGCQkJKCoqMiuWkYjFYkyYMAEVFRU0P4EQMmpOmRDMGT/QioyM1FnPoL+/\nHx0dHWZ3y0RHR6O6uhrd3d1mHT8SiUSC4OBglJWVUWkLQsioOGRCUKlUBhMCy7KjaiGMHz9eZ23l\n9vZ2iMXiEUttG+Lq6gqZTIZLly6ZdbwxtKUtampqrHYNQojzc8iEMNxTRtpv4mKx2Kxzh4aGQi6X\nQ6FQABh4umi0g7bx8fG4dOmSyU9TGYthGERERKCtrW3UC/MQQm5cDpkQhmshtLS0wMfHx+ynevh8\nPvz9/dHU1IT+/n60t7ePOiH4+/uDYRirzh1wcXFBVFQUrl27ZrXuKUKIc3PIhNDX12cwIXR0dIz6\niRtt6eqOjg64u7ubVCBPH4ZhuMdErUksFiMiIgKlpaWUFAghJnPIhDBcC6G9vX3UdX58fX0hl8vR\n2tpq8mQ0Q6Kjo1FaWmr1gV8vLy8uKfT09Fj1WoQQ5+KQCWG4MQRLFH7TTihrb2+3WELw8vKCRCJB\nbW2tRc430rXGjRuHkpISqNVqq1+PEOIcHDIhDNdCsFSXUV1dHUQikcmluYcTGxuLn3/+eUweD/X2\n9oanpyfVPCKEGM0hE4KhMQSWZS3SQvDx8UFjY6PFWgdaiYmJYFkWZ8+eteh5DQkKCkJTU5PF1mQg\nhDg3h0wIhkpXmLNKmj7e3t6Qy+UWTwg8Hg9z585FcXHxmMwZEAqF8PDwoFYCIcQoDpsQ9H3oW2rh\nGIlEgs7OTm4ugiWJRCJMnToVFy9etPi59QkKCkJjYyO1EgghI3LIhGCodIWlEkJnZydkMpnVahCN\nHz8eVVVVUKlUVjn/YO7u7pBKpQYX/iGEEC2nSgiWeORUrVajo6MDKSkpuHDhwqjOZYi7uzsCAgJQ\nWVlplfNfLygoCPX19WOSgAghjsthE4K+MQRLPGGkHUyOj4+36kSyqKioEZcBtRSRSAR/f38UFRWh\npaXFaiU0CCGOzSETgqExhNG2ENRqNRobGxEUFITx48ejsrLSas/xW/v81wsJCYFMJkNdXR3Kyspo\nfgIhZAiHTAjWGkNobGyEp6cnhEIh3N3d4efnN2RtBEvRfmsfXGrb2sRiMeLj4+Hq6oorV65QFxIh\nRIfDJoTru4z6+vqgVqshEonMOqdGo0FTUxOCg4O59+Lj43H8+PFRxTqciRMnorCwcEy7cHg8HsLC\nwuDl5YUrV65AqVSO2bUJIfbN4RKCtqvj+vUJ5HL5qKqcNjY2wsPDA0KhkHtv6dKl+O6779DQ0GB+\nwMOIiIiARqOxWivEEIZhEBISAh8fH1y5coVWWyOEAHDAhGBo/KCpqQkBAQFmnVOj0aCxsVGndQAA\nnp6eSEpKstriNgzDIDU1FYWFhVY5/0iCg4MREBCAS5cuobm5mQabCbnBOVxCMFS2orGx0eyEoK91\noBUbG2vV1c5kMhk6OjrGvOtIKyAgADExMWhsbERpaSl1IRFyA3O4hNDb26v3g7upqQn+/v4mn89Q\n60ArNjYW58+ft9oALJ/Px913342SkhL85z//sUlScHd3R3x8PMRiMYqLiyGXy6m1QMgNyOESgkKh\ngLu7u857fX196OnpMWtls+FaBwAwYcIEhIWFITs726x4jSGRSLBgwQI0NzejuLjYatcZjnZcITo6\nGg0NDdRaIOQG5BQJoampCX5+fuDxTPvjjNQ6AAaeylm6dClOnDhh1W/NAoEA6enpOHPmDDo7O612\nnZGIRCLExcVBLBbj0qVLtPIaITcQp0gI5o4fjNQ60AoPDwcAq1co9fX1xcSJE3Hs2DGbThzj8XgI\nCQnBuHHjUFpaivb2dpvFQggZOw6XEPSNIZgzfmBM60CLYRhMnDjRarWNBktNTYVYLMa//vUvmz8O\n6uXlhQkTJqCiogINDQ00kY0QJ+dwCcFQC8HUhNDQ0GBU60ArIiICVVVVJl3DHDweD3PmzIGvry8O\nHjxo87LVEokEsbGx6OzsxIULF1BcXIza2loqfUGIE3K4hNDb26uTELq6uqBWq00qaqdUKtHY2IjQ\n0FCjjwkPDx+ThAAMtEhuvvlmuLu744cffhiTaw5HKBRCJpMhJSUFYWFhUCqVuHjxIs1dIMTJOFxC\nUCgUOt/q6+rqEBwcbNIM5ZqaGgQEBJi0slpYWBgqKyvHrGQ1wzC47bbbcO3aNRQWFtq8pQAMtF6k\nUikiIyMhk8nQ3NyMS5cuoaWlZUzWiSaEWJdDJoTBLQRtQjBWV1cXOjs7ERgYaNJ1fXx80NnZiaef\nfnrMBlnd3Nwwf/58VFZWIjs7266WwhSLxYiNjUVQUBCam5tx/vx51NTU2HzcgxBiPodPCLW1tQgJ\nCTHqWJZlUVVVhdDQUPD5fJOuyzAMNmzYgF/96lfYsmULzp49a9Lx5vLy8sJdd92F6dOn4/Dhw3b1\nxA/DMPD29kZMTAxiY2PR39+PS5cu4fLly5DL5XbRqiGEGM9l5F3sy+AxhN7eXnR3d8PX19eoY5ua\nmsAwDHx8fMy6dlJSEry8vHD8+HF8+OGHOHr0KJ555hmzC+oZi2EYyGQy9PX14eDBg7j99tvh7+9v\n9euaQigUIjw8HKGhoWhvb4dcLkdlZSUkEgm8vLzg5eWld1EjQoj9cLiEMHgMoa6uDoGBgUZNSGtq\nakJ9fT1iYmJG9UEaFhaGpUuX4r777sPq1atRXV3NzVOwtsTERLAsiyNHjkAgECAgIAC+vr6Ijo4e\n8uSVrfB4PHh7e8Pb2xsajQbt7e1oa2tDTU0NRCIRfHx84OnpCRcXF7tKaIQQB0wIg1sIxo4fNDY2\noqGhAbGxsXBzc7NIHAKBAFOnTkV+fv6YJQRgoJWSmJiIhoYGtLS0oL6+Hj/++CNkMhkmT54MqVQ6\nZrGMhM/nw8fHBz4+Pujv70d7eztaWlq4ct/ahYjEYjFEIhHc3d0pSRBiQw6XELQtBJZlUVlZiVtv\nvdXgvizLoq6uDi0tLYiJibFYMtBKT0/Hxo0b4enpiVtvvXXMukQYhkFQUBCCgoKQkJAAhUKBc+fO\n4YsvvkBUVBSSk5ON7kYbK4NbDizLQq1Wo7e3Fz09Pejs7ERDQwOUSiU8PT3h5eUFDw8P9Pf3Q61W\nQ61WQywWmzzuQwgxjVUHlXNzcxEXF4fo6Ghs3rxZ7z5r165FdHQ0UlJSUFBQMOz5WJaFXC6Hl5cX\nSkpK4OrqiqCgIL37KpVKXLlyBZ2dnRZtGWjl5eVBJpPhhRdewPfff4/Vq1fj0qVLKCsrG/MZve7u\n7pgxYwZCQ0MhkUhw8OBBfPbZZzh8+DBOnTqFwsJCtLW1jWlMg+Xl5em8ZhgGAoEAUqkUgYGBGD9+\nPBITE5GUlAQPDw+0tLTg/PnzKC4uRkVFBWpra3Hu3DmUlpZCLpejs7MT3d3d6OnpgVKpNGsuxPUx\n2QOKyTj2GBNgv3GZwmotBI1Gg9WrV+Pbb79FaGgopk2bhszMTMTHx3P7HDp0CKWlpSgpKcGZM2fw\nxBNP4PTp0wbPWVNTAxcXF/j4+ODIkSO47bbbhnQxsCyL1tZWVFVVISAgAEFBQVbphsjLy0NaWhpk\nMhn+8pe/4Pvvv8e7776Lvr4+SKVSeHh4YOXKlRg3btyYdYP88MMPWL9+PSZPnsx1z3R2dqKjowM5\nOTkQi8Xw8PCASCSCr68vvL29uS6b61egsyTtvRqJQCCAn58f/Pz8hmxTq9XceIRKpQLLsmBZFiqV\nCv39/XBzc4OrqysEAgFcXFwgEAjg5uYGNzc3CAQC8Hg8nf8PxsY0ligm49hjTID9xmUKq30K5Ofn\nQyaTITIyEgCwePFi5OTk6CSE/fv348EHHwQAzJgxA21tbWhoaDA4R6CwsBDJyckoLi6Gt7c397hp\nf38/lEol5HI55HI5XF1dMWHCBEgkEmv98Yb41a9+hV/96ldQqVTIz89HSUkJ1q9fj6SkJFy7do37\nkBKJRLj33nsBDMx+dnV1hUgkgkKhADDw7dnYchqGDO6e0br55ptRX1+P7u5udHd3o7a2FkVFRVAo\nFFAoFPDw8OC6ZViWRX9/P3x8fODv7w83Nze4uLjA3d0dQqGQ+4Dt7+8HwzDg8XjQaDRQqVTg8Xjg\n8/kWTzAuLi7w9fXV2xWmVqvR19cHpVLJdTEpFAq0t7ejr6+PSxrae+Pi4oKmpiZcvnyZeyCBYRjw\n+XzuR5tAtO+7uLhw2wbr7+/n7pf2h2VZ7v+39l5pr6H97+DfAXAJTqPRgGVZ8Pl8Gk8hY85qCaGm\npkZnsDUsLAxnzpwZcZ/q6mq9CaGwsBB79+7Fs88+CxcXF8yYMQMXL17k/rELBAJ4eXnZ/IkbgUCA\nWbNmYdasWbj77ruxbds2PPbYY9wHS1lZGfbs2QO1Wo3q6mq4uLjA09MTra2t4PF46O3tha+vL1iW\nhbe3N5RKJbcugaenJ/r7++Hu7o62tjYwDAO1Wg1fX1/w+Xx88803UKlUEAgEmDNnDm655RYuLm0F\nU300Gg1aWlqgUCigVqu5DzC5XM51galUKigUCvT29nLf0BmG4T7IeDweXF1duX5/YGDQWFtGe8+e\nPdy8BO0HZV9fH7cCnqurK/dheP0HpvbDWZtseDwed13tB7BGo+G2X3+sdh/tcQDQ2tqKy5cvc+cA\n/vehrI1Ba/C1NBoNd17tNbT/vf4Y7XW1r0dy4cIFZGdnIzg4GB4eHkP+7IYShL57Zuh37Z/l+tgH\nxzn4OnK5XO+KgfqS2uB4rt9X3/v67sn191GflpYWlJaWGjzP4NfGnM+YmIwxOC5D/7+HO9fgY2Qy\nmW2+ELBW8uWXX7KPPvoo9/rTTz9lV69erbPPb37zG/bkyZPc6zlz5rBnz54dci4A9EM/9EM/9GPG\njyms1kIIDQ3VKQZXVVWFsLCwYfeprq7WW3COpQJqhBBidVZ7ymjq1KkoKSlBRUUFlEolsrOzkZmZ\nqbNPZmYmPvnkEwDA6dOn4eXlZXKNIUIIIZZhtRaCi4sLtm7dinnz5kGj0WDFihWIj4/H9u3bAQCr\nVq3CnXfeiUOHDkEmk0EsFuOf//yntcIhhBAyEtNGBsbW4cOH2djYWFYmk7GvvfaarcPhREREsMnJ\nyeykSZPYadOm2SSGhx9+mA0ICGCTkpK49+RyOXv77bez0dHR7Ny5c9nW1labx/TSSy+xoaGh7KRJ\nk9hJkyaxhw8fHtOYKisr2bS0NDYhIYFNTExk33vvPZZlbXuvDMVk63ulUCjY6dOnsykpKWx8fDz7\n5z//mWVZ294rQzHZ+l6xLMuq1Wp20qRJ7G9+8xuWZW3/709fTKbeJ7tNCGq1mp0wYQJbXl7OKpVK\nNiUlhS0qKrJ1WCzLsmxkZCQrl8ttGsPx48fZn3/+WefD95lnnmE3b97MsizLvvbaa+xzzz1n85jW\nr1/PvvXWW2Max2B1dXVsQUEBy7Is29nZycbExLBFRUU2vVeGYrL1vWJZlu3u7mZZlmVVKhU7Y8YM\n9sSJEzb/e6UvJnu4V2+99Ra7dOlS9q677mJZ1vb//vTFZOp9stvy14PnMQgEAm4eg71gbTzQfcst\nt+jMMwB053U8+OCD+Prrr20eE2DbexUUFIRJkyYBGFgOND4+HjU1NTa9V4ZiAmz/90okEgEYmOmv\n0Wjg7e1t879X+mICbHuvqqurcejQITz66KNcHLa+T/piYgc9ZmwMu00I+uYoaP/R2BrDMLj99tsx\ndepU/OMf/7B1OJzBk/oCAwPR0NBg44gGbNmyBSkpKVixYoVNS2hUVFSgoKAAM2bMsJt7pY3ppptu\nAmD7e9Xf349JkyYhMDAQ6enpXCFFW94rfTEBtr1XTz/9NN544w2dSsu2vk/6YmIYxqT7ZLcJwZ5n\naZ46dQoFBQU4fPgw/va3v+HEiRO2DmkIUyfkWMsTTzyB8vJyFBYWIjg4GH/6059sEkdXVxcWLlyI\n9957b0hFWFvdq66uLixatAjvvfceJBKJXdwrHo+HwsJCVFdX4/jx4/juu+90ttviXl0fU15enk3v\n1b/+9S8EBAQgNTV12AloY3mfDMVk6n2y24RgzDwGW9GW3Pb398c999yD/Px8G0c0IDAwEPX19QAG\nSoMHBATYOCIgICCA+8fx6KOP2uReqVQqLFy4EMuXL8eCBQsA2P5eaWP63e9+x8VkD/dKy9PTE/Pn\nz8fZs2dtfq+uj+mnn36y6b36/vvvsX//fowfPx5LlizBsWPHsHz5cpveJ30xPfDAAybfJ7tNCMbM\nY431kx8AAANhSURBVLAFbblmAOju7saRI0eQnJxs46gGZGZmYufOnQCAnTt3ch80tlRXV8f9/tVX\nX435vWJZFitWrEBCQgL+8Ic/cO/b8l4ZisnW96q5uZnrUlAoFPjmm2+Qmppq03tlKCbtBy8w9vdq\n06ZNqKqqQnl5OT7//HPcdttt+PTTT216n/TF9Mknn5j+d8qSI9yWdujQITYmJoadMGECu2nTJluH\nw7Isy5aVlbEpKSlsSkoKm5iYaLO4Fi9ezAYHB7MCgYANCwtjP/roI1Yul7Nz5syx2WNv18e0Y8cO\ndvny5WxycjI7ceJE9u6772br6+vHNKYTJ06wDMOwKSkpOo/e2fJe6Yvp0KFDNr9X586dY1NTU9mU\nlBQ2OTmZff3111mWZW16rwzFZOt7pZWXl8c90WPrf39a3333HRfT7373O5PuE8OyVBeCEEKIHXcZ\nEUIIGVuUEAghhACghEAIIeS/KCEQQggBQAmBEIM2btyIpKQkpKSkIDU1Ffn5+UhLS8O0adO4fX76\n6Sekp6cDGFhT19PTE6mpqUhISMC6detsFTohZrHeyuqEOLAffvgBBw8eREFBAQQCAVpaWtDX1weG\nYdDU1ITc3FzccccdQ46bPXs2Dhw4gN7eXqSmpuKee+7BlClTbPAnIMR01EIgRI/6+nr4+flBIBAA\nAHx8fLgZ6llZWdi4ceOwxwuFQkyaNAllZWVWj5UQS6GEQIgev/71r1FVVYXY2Fg8+eSTOH78OLdt\n5syZcHV1RV5ensF6NS0tLcjPz0dCQsJYhUzIqFFCIEQPsViMs2fP4v/+7//g7++P+++/nytLAADr\n1q3DK6+8MuS4EydOYNKkSQgPD8eCBQu4ypyEOAJKCIQYwOPxcOutt2L9+vXYunUr9u7dC2CgkmV6\nejoUCgVOnz6tc8wtt9yCwsJCXLx4Efv27dMp0EiIvaOEQIgeV65cQUlJCfe6oKAAERERAP63MMu6\ndeuwefNmvd1GkZGReOqpp/Dyyy+PTcCEWAA9ZUSIHl1dXVizZg3a2trg4uKC6OhobN++HYsWLeIS\nQEZGhk6J4+tr4D/++OOIiYlBdXW13ZRuJ2Q4VNyOEEIIAOoyIoQQ8l+UEAghhACghEAIIeS/KCEQ\nQggBQAmBEELIf1FCIIQQAgD4/30oawe/0i41AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10075ff50>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = plt.subplots(1)\n",
      "\n",
      "m = [np.median(x) for x in [snr1k, snr2k, snr4k]]\n",
      "e_up = [np.abs(np.median(x)-stats.scoreatpercentile(x,0.95)) for x in [snr1k, snr2k, snr4k]]\n",
      "e_down = [np.abs(np.median(x)-stats.scoreatpercentile(x,0.05)) for x in [snr1k, snr2k, snr4k]]\n",
      "            \n",
      "ax.errorbar([1000,2000,4000], m, yerr=np.vstack([e_up, e_down]), \n",
      "            fmt = 'o', color='k')\n",
      "\n",
      "ax.set_ylabel('SNR')\n",
      "ax.set_xlabel('b value')\n",
      "ax.set_xlim([900, 4500])\n",
      "ax.set_ylim([0, 40])\n",
      "ax.set_xticks([1000,2000,4000])\n",
      "\n",
      "\n",
      "ax.loglog()\n",
      "fig.savefig('figures/snr_bars.svg')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAENCAYAAADt3gm6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEIhJREFUeJzt3X9MleX/x/HXSRBT1kep5YpzFoYoklg5zGQjsbZOSlHT\nIdDsB2TDHyBr1WbTJmVu0dqaaM02yzVnBAYb1plnztFZ+YcwTDcXGmcok3CuNra2jB8K5/tHcb4c\nOceOdM59n7vzfGxn474P5+Z9zs6uF9d93dd123w+n08AgLh2m9kFAADMRxgAAAgDAABhAAAQYQAA\nEGEAAJCUYHYBk2Gz2cwuAQAsKdRsAkv2DHbs2CGfzxfxR7SOyyM+Hnx/+Eyt8H5DsWQY5Ofnm10C\nMAHfy8iLt8/UzPdLGAARwvcy8uLtMyUMYkS8ffEAYIzNd7OTSACAuEDPAABAGAAACAMAgAgDAIAs\nGgYej0c1NTXyeDxmlwIAMS+cNpOriQAA1uwZAAAiizAAABAGAADCAAAgwgAAIMIAACDCAAAgwgAA\nIMIAACDCAAAgwgAAIMIAACDCAAAgwgAAIIuGAfczAIDwcT8DAEBYLNkzAABEFmEAACAMAACEAQBA\nhAEAQIQBAECEAQBAhAEAQIQBAECEAQBAhAEAQIQBAEBSgtkFmM3j8cjj8airq0tdXV16+umnJUn5\n+fnKz883tzgAMAirlv6tsbFRX3/9tRobG80uBQAMx2kiAABhAAAgDAAAIgwAAIrBMLh48aLWr1+v\noqIis0sBgLgRc2EwZ84c7d+/3+wyACCuGBIG5eXlmj17trKzswP2u91uZWZmKiMjQ7W1tUaUAgAI\nwpAwKCsrk9vtDtg3MjKiyspKud1udXZ2qr6+XufOnTOiHADADQwJg7y8PM2aNStgX3t7u+bOnau0\ntDQlJiaqpKRELS0t6u/v14YNG3TmzBl6CwBgENOWo+jr65PD4fBv2+12tbW1KSUlRfv27bvpa8eW\nkBjD0hEAMNGttJWmhYHNZpv0a2n8AeCf3UpbadrVRKmpqert7fVv9/b2ym63m1UOAMQ108IgJydH\nXq9XPT09Gh4eVkNDgwoLC80qBwDimiFhUFpaqtzcXHV1dcnhcOjAgQNKSEjQ3r175XQ6lZWVpeLi\nYi1YsCCs43k8HtXU1AScCwMABBdOm8kS1n9jCWsA8SzmZiADAIxHGAAACAMAgEXDINIDyC6XS7t2\n7ZLH45HT6ZTL5YrIcQEgFjCAHAaXy6Xq6mp1d3f796Wnp2v37t0qKCgwsTJYic/n0+joqKZMmWJ2\nKcCkxH0YOJ1OHTt2LOj+GxfXA240Nt2/v79fhw4dUlVVlSRmycN6TFuOIlYMDQ0F3T84OGhwJbCi\nsUbf6/Xq6NGjqqmpMbskYFIsOWYQSUlJSUH3T5s2zeBKAMA8lgyDSA4gb9myRenp6QH70tPT/d19\nALA6BpDD5HK5tG3bNl2+fFmLFy9WVVUVg8e4JV6vV6tWrZLX6zW7FGBS4n7MQJIKCgp09epVlqMA\nELcseZoIABBZhAEAwJphwBLWABC+cNpMS44ZMKEHAMIXTptpyZ4BACCyCAMAAGEAACAMAAAiDAAA\nsmgYcGkpAISPS0sBAFxaCgAID2EAACAMAACEAfCvuVwulZeXq6+vT06nUy6Xy+ySgFtmyQFkIFa4\nXC5VV1eru7tbknTs2DH/z9wgCVZCzwD4F+rq6vyN/5ju7m7t2bPHpIqAySEMgH9haGgo6P7BwUGD\nKwH+HUuGAZPOECuSkpKC7p82bZrBlQChMekMiLItW7aou7s74FRRenq6qqqqTKwKCBROm2nJMABi\nxdggcW1trTo6OvTYY4+pqqqKwWNYDmEA/EsFBQWaN2+eVq1aJbfbbXY5sBiPxyOPx6N9+/bplVde\nUWJioiTjz4AQBgBgorFG/8MPP9Rbb72l5ORkU+qw5AAyACCyCAMAAGEAACAMAAAiDAAAsmgYMAMZ\nAMLHDGQAALe9BACEhzAAABAGAADCAAAgwgAAoEmGwS+//BLpOgAAJrrppaWnTp3ShQsXlJWVpQce\neEC9vb3auXOn3G63Ll26ZFSNUTW2fGx/f7+mTp2qmpoaSVy+CiC+hAyD7du3q6mpSQ899JC2bt2q\n5557Ts3NzaqurlZdXZ2RNUYVjT4A3CQMmpubdfr0aU2bNk39/f1yOBz66aeflJaWZmB5AAAjhBwz\nSEpK8t/UOyUlRRkZGQQBAESBy+XSwMCAVq5cKafTKZfLZXgNIXsGFy5c0DPPPOPf7unp8W/bbDYd\nOXIk+tUBwH+cy+VSdXW1RkdHdeLECUlSd3e3JBl6L+2QYdDS0hKw/frrr/t/ttls0asIAOJIXV2d\nv/Ef093drT179sRGGDCoCgDRNzQ0FHT/4OCgoXWEDIMVK1YE3T/WK2htbY1ORWEYuxyUK4EAWF1S\nUlLQ/WNjtpEQTptp8/l8vmBPdHR0/P8v/R0AJ0+eVG1tre6+++6A54F45/V6tWrVKnm9XrNLgcWM\njRmMP1WUnp6u3bt3x8ZpopycHP/PHo9H7733ngYGBvTpp59q5cqVhhQHAP91Yw1+YWGhcnNzNWPG\nDFVVVRkaBNI/zEB2u93atWuXpk6dqu3bt4c8dQQAmLyCggLdfvvtOnr0qJKTk02pIWQYLFmyRL/9\n9pveeOMNLVu2TJL0448/yufzyWazafHixYYVCQCIrpBhMGPGDCUnJ6upqUlNTU0Tnv/uu++iWhgA\nwDghw+CDDz6Qw+HQPffcI0n64osv1NTUpPvuu8+/mBsA4L8h5HIUFRUV/kuevv/+e23dulUvvfSS\n/ve//6miosKwAgEA0ReyZzA6OqqUlBRJUkNDgyoqKrRmzRqtWbNGDz74oGEFAgCiL2TPYGRkRNeu\nXZMkHT9+POBKouvXr0e/MgCAYUL2DEpLS7V8+XLdddddmj59uvLy8iT9Nblm5syZhhUIAIi+kGGw\nbds2Pf7447py5YqefPJJ3XbbX50In8+nPXv2GFYgACD6bjrpbGx+wXjz5s2LWjEAAHOEHDMAAMQP\nwgAAQBgAAAgDAIAIAwCACAMAgAgDAIAIAyAiZs6cqbKyMrPLACYt5D2QzXL16lVt2rRJSUlJys/P\n1/PPP292SUBIYzcav9HNbjwOBJOcnKwrV66YdqezmAuDgwcPKiUlRQUFBSopKdFXX31ldkkAEHVm\nh4Ehp4nKy8s1e/ZsZWdnB+x3u93KzMxURkaGamtrJUl9fX1yOBySpClTphhRHgDEPUPCoKysTG63\nO2DfyMiIKisr5Xa71dnZqfr6ep07d052u129vb2S/rqnAgAg+gwJg7y8PM2aNStgX3t7u+bOnau0\ntDQlJiaqpKRELS0tWr16tZqamrRp0yYVFhYaUR4AxL2brloaTeNPB0mS3W5XW1ubpk+frs8///ym\nr71x0I7BOgCY6FbaStPCwGazTfq1NP4A8M9upa00bZ5Bamqqf2xAknp7e2W3280qBwDimmlhkJOT\nI6/Xq56eHg0PD6uhoYExAgAwiSFhUFpaqtzcXHV1dcnhcOjAgQNKSEjQ3r175XQ6lZWVpeLiYi1Y\nsCCs43k8HtXU1ASd7AMACBROmxlzk84AIB7FxaQzAEBsIwwAANYMA8YMACB8jBkAgEUwZgAAMB1h\nAAAgDAAAFg0DBpABIHwMIAOARTCADABQUVGREhJMW0iangEAmOnGew6MMXqpfsIAAMBpIgCARcOA\nq4kAIHxcTQQACIslewYAgMgiDAAAhAEAgDAAAIgwAADIomHApaUAED4uLQUAhMWSPQMAQGQRBgAA\nwgAAQBgAAEQYAABEGAAAZNEwYJ4BAISPeQYAgLBYsmcAAIgswgAAQBgAAAgDAIAIAwCACAMAgAgD\nAIAIAwCACAMAgCwaBixHAQDhYzkKAEBYLNkzAABEFmEAACAMAACEAQBAhAEAQIQBAECEAQBAhAEA\nQIQBAECEAQBAhAEAQIQBAECEAQBAFg0DlrAGgPCxhDUAICyW7BkAACKLMAAAEAYAAMIAACDCAAAg\nwgAAIMIAACDCAAAgwgAAIMIAACDCAAAgwgAAIMIAACDCAAAgwgAAIMIAACDCAAAgwgAAIMIAAKAY\nDIOLFy9q/fr1KioqMrsUAIgbMRcGc+bM0f79+0352x6Px5S/CwBmi1oYlJeXa/bs2crOzg7Y73a7\nlZmZqYyMDNXW1kbrz08KYQAgXkUtDMrKyuR2uwP2jYyMqLKyUm63W52dnaqvr9e5c+d08OBBvfba\na7p8+XJYx6bRRiziexl58faZmvl+oxYGeXl5mjVrVsC+9vZ2zZ07V2lpaUpMTFRJSYlaWlr0wgsv\n6KOPPtK9996r/v5+bdiwQWfOnAnZc4i3Lwisge9l5MXbZ2rq+/VF0cWLF30LFy70bx8+fNi3fv16\n//bBgwd9lZWVt3xcSTx48ODBYxKPUBJkIJvNFpHj/JUHAIBIMfRqotTUVPX29vq3e3t7ZbfbjSwB\nABCEoWGQk5Mjr9ernp4eDQ8Pq6GhQYWFhUaWAAAIImphUFpaqtzcXHV1dcnhcOjAgQNKSEjQ3r17\n5XQ6lZWVpeLiYi1YsCBaJQAAwmTzcQIeAOJezM1AjiXnz5/Xxo0btXbtWn322WdmlwMAUUPPIAyj\no6MqKSlRY2Oj2aUAQFTEXc/gVpfJ+Oabb1RQUKCSkhKjSwUAw8Rdz+CHH35QcnKyXnzxRZ09e1bS\nX8tkzJ8/X8ePH1dqaqqWLFmi+vr6gMHtZ599Vi0tLWaVDQBRZeiks1iQl5ennp6egH3jl8mQ5F8m\n49dff1Vzc7MGBwe1YsUK44sFAIPEXRgE09fXJ4fD4d+22+1qa2vT8uXLtXz5chMrAwBjxN2YQTCR\nWiYDAKyKMBDLZAAAYSCWyQCAuAsDlskAgIni7tJSAMBEcdczAABMRBgAAAgDAABhAAAQYQAAEGEA\nABBhAAAQYQD49fT0TLjPRSS8/PLLampqivhxgUgiDIAos9lsLIaImEcYAONcv35d69atU1ZWloqK\nijQwMBDw/Pnz57V06VL/dk9PjxYtWiRJevfdd/XII48oOztbFRUVQY+flpam/v5+SVJHR4f/PhlX\nr15VeXm5li5dqsWLF+vIkSPReHtASIQBMM7PP/+szZs3q7OzU3fccYc++eSTgOczMzM1PDzsv0FS\nQ0OD/5aoVVVVam9v19mzZzUwMKBvv/12wvFD9RB27dqlJ554Qm1tbWptbdWbb76pP//8M7JvDrgJ\nwgAYx+FwaNmyZZKkdevW6cSJExN+Z+3atWpoaJAkNTY2qri4WJLU2tqqRx99VIsWLVJra6s6OzvD\n/rvHjh3T+++/r4cfflgrVqzQ0NBQwLLqQLRxpzNgnPH/uft8vqD/yRcXF6uoqEirV6+WzWZTenq6\nBgcHtXnzZp06dUqpqal65513NDg4OOG1CQkJGh0dlaQJzzc3NysjIyPC7wgIDz0DYJxLly7p5MmT\nkqQvv/xSeXl5E37n/vvv15QpU7Rz507/KaKxhv3OO+/UH3/8ocOHDwc9flpamjo6OiQp4Aojp9Op\nuro6//bp06cj84aAMBEGwN9sNpvmz5+vjz/+WFlZWfr999+1cePGoL9bXFysQ4cOae3atZKkmTNn\n6tVXX9XChQv11FNPBQwyj7djxw5VV1dryZIlSkhI8Pc83n77bV27dk2LFi3SwoULtWPHjui8SSAE\n7mcAAKBnAAAgDAAAIgwAACIMAAAiDAAAIgwAAJL+D7UFUrHwI1eXAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x108d93390>"
       ]
      }
     ],
     "prompt_number": 12
    }
   ],
   "metadata": {}
  }
 ]
}