{
 "metadata": {
  "name": "04-diginorm"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Change what's inside the quotes into your MBL STAMPS server username, or something else unique\n",
      "# reminder! hit shift-ENTER to execute cell & move onto next cell\n",
      "username=\"CHANGEME\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# configure paths etc.\n",
      "import sys\n",
      "import os\n",
      "os.environ['PYTHONPATH'] = '/usr/local/share/khmer/python'\n",
      "sys.path.insert(0, '/usr/local/share/khmer/python')\n",
      "\n",
      "import khmer\n",
      "import screed"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cd /mnt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "/mnt\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mkdir -p $username"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cd $username"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "/mnt/CHANGEME"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You should now be in '/mnt/USERNAME', where USERNAME is *not* CHANGEME"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Make two random genomes with a 10:1 abundance"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import random\n",
      "random.seed(1)\n",
      "\n",
      "x = [\"A\"] + [\"G\"] + [\"C\"] + [\"T\"]\n",
      "x = x*1000\n",
      "random.shuffle(x)\n",
      "x = \"\".join(x)\n",
      "\n",
      "y = [\"A\"] + [\"G\"] + [\"C\"] + [\"T\"]\n",
      "y = y*1000\n",
      "random.shuffle(y)\n",
      "y = \"\".join(y)\n",
      "\n",
      "print 'x is', x[:100]\n",
      "print 'y is', y[:100]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "x is CAGACTTGGAAGCTGAGAGTCCGACGTCACTGCCTCAACTCGCGCAAATGTTCCCGCCAAATTGTATCCTAGGGATCTTCCATAAGCTTATATACGGGGG\n",
        "y is ATTAGACCACATGTTGCTAGCTATTCGGAGGAAGGCAGTTTAGCTAAATATGTGAGGGCTGTGGGACGAATATCGGGGGATGACGGTCCCATAGCTAGTT\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "outfp = open('metagenome.fa', 'w')\n",
      "print >>outfp, \">x 1\"\n",
      "print >>outfp, x\n",
      "print >>outfp, \">y 2\"\n",
      "print >>outfp, y\n",
      "outfp.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!python /usr/local/share/2012-paper-diginorm/pipeline/make-biased-reads.py metagenome.fa | head -100000 > reads.fa"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Traceback (most recent call last):\r\n",
        "  File \"/usr/local/share/2012-paper-diginorm/pipeline/make-biased-reads.py\", line 56, in <module>\r\n",
        "    print '>read%d\\n%s' % (i, read)\r\n",
        "IOError: [Errno 32] Broken pipe\r\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(Yes, you should see an error.)\n",
      "\n",
      "## Now, apply digital normalization & error trimming"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!/usr/local/share/khmer/scripts/normalize-by-median.py -k 20 -C 20 -x 1e8 reads.fa --savehash normC20k20.kh"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r\n",
        "PARAMETERS:\r\n",
        " - kmer size =    20 \t\t(-k)\r\n",
        " - n hashes =     4 \t\t(-N)\r\n",
        " - min hashsize = 1e+08 \t(-x)\r\n",
        " - paired = False \t\t(-p)\r\n",
        "\r\n",
        "Estimated memory usage is 4e+08 bytes             (n_hashes x min_hashsize)\r\n",
        "--------\r\n",
        "making hashtable\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "DONE with reads.fa; kept 3787 of 50000 or  7%\r\n",
        "output in reads.fa.keep\r\n",
        "Saving hashfile through reads.fa\r\n",
        "...saving to normC20k20.kh\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "fp rate estimated to be 0.000\r\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!/usr/local/share/khmer/scripts/filter-abund.py normC20k20.kh reads.fa.keep"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "file with ht: normC20k20.kh\r\n",
        "loading hashtable\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "K: 20\r\n",
        "filtering reads.fa.keep\r\n",
        "starting threads\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "starting writer\r\n",
        "loading...\r\n",
        "... filtering 0\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done loading in sequences\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "DONE writing.\r\n",
        "processed 3787 / wrote 3069 / removed 718\r\n",
        "processed 378700 bp / wrote 211885 bp / removed 166815 bp\r\n",
        "discarded 44.0%\r\n",
        "output in reads.fa.keep.abundfilt\r\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!/usr/local/share/khmer/scripts/normalize-by-median.py -k 20 -C 5 -x 1e8 reads.fa.keep.abundfilt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r\n",
        "PARAMETERS:\r\n",
        " - kmer size =    20 \t\t(-k)\r\n",
        " - n hashes =     4 \t\t(-N)\r\n",
        " - min hashsize = 1e+08 \t(-x)\r\n",
        " - paired = False \t\t(-p)\r\n",
        "\r\n",
        "Estimated memory usage is 4e+08 bytes             (n_hashes x min_hashsize)\r\n",
        "--------\r\n",
        "making hashtable\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "DONE with reads.fa.keep.abundfilt; kept 1059 of 3069 or 34%\r\n",
        "output in reads.fa.keep.abundfilt.keep\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "fp rate estimated to be 0.000\r\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ls"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "group0.hist    reads.fa                      reads.group0000.fa\r\n",
        "group1.hist    reads.fa.keep                 reads.group0001.fa\r\n",
        "metagenome.fa  reads.fa.keep.abundfilt       reads.info\r\n",
        "normC20k20.kh  reads.fa.keep.abundfilt.keep\r\n",
        "reads.dist     reads.fa.part\r\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!/usr/local/share/khmer/scripts/abundance-dist-single.py -x 1e8 -k 20 reads.fa.keep reads.fa.keep.hist\n",
      "!/usr/local/share/khmer/scripts/abundance-dist-single.py -x 1e8 -k 20 reads.fa.keep.abundfilt reads.fa.ka.hist\n",
      "!/usr/local/share/khmer/scripts/abundance-dist-single.py -x 1e8 -k 20 reads.fa.keep.abundfilt.keep reads.fa.kak.hist"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r\n",
        "PARAMETERS:\r\n",
        " - kmer size =    20 \t\t(-k)\r\n",
        " - n hashes =     4 \t\t(-N)\r\n",
        " - min hashsize = 1e+08 \t(-x)\r\n",
        "\r\n",
        "Estimated memory usage is 4e+08 bytes (n_hashes x min_hashsize)\r\n",
        "--------\r\n",
        "making hashtable\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "building tracking ht\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "K: 20\r\n",
        "HT sizes: [100000007, 100000037, 100000039, 100000049]\r\n",
        "outputting to reads.fa.keep.hist\r\n",
        "consuming input, round 1 -- reads.fa.keep\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "preparing hist from reads.fa.keep...\r\n",
        "consuming input, round 2 -- reads.fa.keep\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r\n",
        "PARAMETERS:\r\n",
        " - kmer size =    20 \t\t(-k)\r\n",
        " - n hashes =     4 \t\t(-N)\r\n",
        " - min hashsize = 1e+08 \t(-x)\r\n",
        "\r\n",
        "Estimated memory usage is 4e+08 bytes (n_hashes x min_hashsize)\r\n",
        "--------\r\n",
        "making hashtable\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "building tracking ht\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "K: 20\r\n",
        "HT sizes: [100000007, 100000037, 100000039, 100000049]\r\n",
        "outputting to reads.fa.ka.hist\r\n",
        "consuming input, round 1 -- reads.fa.keep.abundfilt\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "preparing hist from reads.fa.keep.abundfilt...\r\n",
        "consuming input, round 2 -- reads.fa.keep.abundfilt\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r\n",
        "PARAMETERS:\r\n",
        " - kmer size =    20 \t\t(-k)\r\n",
        " - n hashes =     4 \t\t(-N)\r\n",
        " - min hashsize = 1e+08 \t(-x)\r\n",
        "\r\n",
        "Estimated memory usage is 4e+08 bytes (n_hashes x min_hashsize)\r\n",
        "--------\r\n",
        "making hashtable\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "building tracking ht\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "K: 20\r\n",
        "HT sizes: [100000007, 100000037, 100000039, 100000049]\r\n",
        "outputting to reads.fa.kak.hist\r\n",
        "consuming input, round 1 -- reads.fa.keep.abundfilt.keep\r\n",
        "preparing hist from reads.fa.keep.abundfilt.keep...\r\n",
        "consuming input, round 2 -- reads.fa.keep.abundfilt.keep\r\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dn1 = numpy.loadtxt('reads.fa.keep.hist')\n",
      "abund = numpy.loadtxt('reads.fa.ka.hist')\n",
      "dn2 = numpy.loadtxt('reads.fa.kak.hist')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(dn1[:,0], dn1[:,1], label='first round')\n",
      "plot(abund[:,0], abund[:,1], label='trim errors')\n",
      "plot(dn2[:,0], dn2[:,1], label='second round')\n",
      "axis(ymax=2500)\n",
      "legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "<matplotlib.legend.Legend at 0x7f60c00f6d10>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD9CAYAAAC7iRw+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXlcVPX6x98z7CgibrggqYkirqQI7khpZimmt9TMFq2b\n2L2aZdtts65t5s203CrNcqnftfKmlmla5FaohWmIipahlguuIDDDzJzfH0dGkJlhgBnmDD7v12te\n45zl+/0enPmc5zzf5/s8OkVRFARBEIRrBr2nByAIgiBULyL8giAI1xgi/IIgCNcYIvyCIAjXGCL8\ngiAI1xgi/IIgCNcYDoX/6NGj9O/fn/bt25OYmMiKFSsAmDZtGhEREcTGxhIbG8u6deus58yZM4eo\nqChiYmLYunWrdXtmZiY33HADrVq14plnnnHT5QiCIAjloXMUx3/ixAlOnDhBly5dyMnJoXv37vzy\nyy+8+eabhISE8Oijj5Y6/tSpU/Tt25cNGzbw+++/M2XKFH7++WcABg8ezL333stNN91EcnIyb731\nFt26dXPv1QmCIAhl8HW0s3HjxjRu3BiABg0a0L59e3bu3AmArftFWloagwYNIjIyksjISBRFIS8v\nj9q1a3PgwAFGjhwJwPDhw0lLSxPhFwRB8ABO+/gPHTpERkYG8fHxALz99tskJCTw+uuvk5ubC8CO\nHTto166d9Zy2bduSlpbGoUOHaNSokXV7TEwMP/74o6uuQRAEQagADi3+YnJzcxk5ciSzZs2iVq1a\npKSk8Pzzz3Px4kUef/xxFi5cyNSpU20+Beh0ujLb7HmXbB0rCIIglE9Fsu+Ua/EXFRUxYsQIxo4d\nS3JyMgCNGjVCp9MRGhrKww8/zKpVqwCIj49n37591nP3799PXFwcrVu35uTJk9bt+/btIyEhwe7g\nS76a/acZr255tcz24ldGhkK7drb3ae31wgsveHwMcn1yfdfatV0L11dRHAq/oiiMHz+eDh068Mgj\nj1i3//XXXwCYTCZWrFjB4MGDAejevTvr168nOzub1NRU9Ho9ISEhAERHR/PJJ5+Qk5PDqlWrrC6j\n8jCajRhMBrv7fXzAZHKqKUEQBIFyXD3btm1j2bJldOrUidjYWABeeeUVPv74Y3bv3o2/vz99+/Yl\nJSUFgPDwcFJSUkhKSsLf35+FCxda25o5cyZ33303Tz/9NKNGjXJ6YtdgNmAw2xd+X18wm51qShAE\nQaAc4e/duzcWi6XM9ltuucXuOZMnT2by5MlltsfExFhDOyuC0Wx0KPzeZPEnJiZ6eghuRa7Pe6nJ\n1wY1//oqisM4/upGp9OV8Vf5vORDSrcU3hn8js1zjh2DhAT1XRAE4VrElnY6wqmoHk9htpixKJZy\nXT3eYvELgqupV68e586d8/QwhGoiLCyMs2fPVrkdTQt/seCXN7krPn7hWuXcuXOViuoQvBNXhbxr\nOkmb0WwEEItfEATBhXiH8Es4pyAIgsvwDuGXcE5B8DqOHDnCHXfcQVhYGG+//TYpKSlMnz7d08Ny\nC4mJiSxatMjTw3Aabfv4Tc75+MXiFwTtsWzZMkJDQzlz5gx6fdVszBYtWrB48WKSkpJcNDrXotPp\nvCrljFj8giC4ha1btxIfH++U6JvKsd4qGq5YXnvXOt4h/A4sfr0eLBaQwAZB0A5JSUls3LiRSZMm\nUadOHbKysrjvvvt47rnnAEhNTSUiIoIFCxbQpk0bxo8fT0FBAQ888AAtWrSgfv369OvXD0VRGDt2\nLNnZ2QwZMoSQkBBmzpxZpj9b7ZlMJpYtW0ZCQgI9evRg+fLl1hvCkiVL6NOnT6k29Ho9v/32GwD3\n3Xcfjz76KCNHjiQ8PJwHHniA48ePW4/NyMjgzjvvpGnTptbCUt4UXaV54ff38Xdo8et0EtIpCFrj\n22+/pU+fPsydO5eLFy8SFRVVxh1y8uRJdu7cyebNm1m4cCGLFy+moKCAPXv2cOrUKV555RV0Oh1L\nly4lMjKStWvXkpuby9SpU232eXV7S5cu5Y033uCdd95h7ty5zJw5k2XLljl9DYsWLWLUqFFkZGRw\n8uRJFixYAKgCn5SURHx8PHv27KGwsJDt27eLq8dVGMwG6gTUodBU6PA4CekUBG1ytRVc8rPZbGba\ntGk0btyYwMBALBYLOTk5HD9+HB8fH3r16lWhvq5u73//+x8PP/ww3bp144YbbiiVSdgZ+vXrx+23\n306DBg0YM2YM33zzDQA7d+4kICCAxx57jAYNGjB9+nSvEn3QuPAbzUZC/EMcunpAJngFwRE6nWte\nlevb/onh4eE0b97c+nn8+PEkJiZy22230bFjxwpHyVzd3vbt2+natav1c9euXdmyZYvT4+7SpYv1\nc+PGja2unrS0NDp37mzdFxQURHR0dIXG6mm0L/wBIQ5dPSATvILgCEVxzcsVlLwR+PqWDioMDg7m\n6aef5vDhwyxevJhHH33UWt/Dx8enXB/61e316tWLXbt2WT/v2rWLvn37AtCsWbNSNULS09PLtGev\nv+La48UUFBSwf/9+h2PTGpoX/joBdcTiFwQvpaR4llc05Msvv+TQoUNYLBZq1aqFv78/gYGBgGqt\n//TTTxXqOzk5mQULFvDTTz+Rnp7OggULGDZsGKBmHj569CgbNmzg6NGjzJgxw+64ryYuLg6DwcCs\nWbM4ffo0zz33nFdN7ILGhd9gUn38YvELgndS0sK/enL3ajdQVlYWAwYMIDQ0lAcffJDp06fTqlUr\nACZMmMDatWupV68eb775Zrl9Adx9991MmTKFiRMnkpKSwiOPPMKYMWMA1T3z7rvvMnXqVG6++WZG\njRrlcKwl29fr9WzcuJFt27bRuXNnAgICKjwf4Wk0nZZ5+Z7lrDm4hv9m/Bfz82a7/sKmTWHXLvVd\nEK4lKhrfLng39v6/K/o90LTFbzQbCfILQq/TY7LY9+VIOKcgCILzaF74/X38CfANkAydgiAILkLT\nwm8wG1Th9wmovnw9CxfChQsuakwQBEF7aFr4jWYjAT4BTln8LnP1vPYa/PCDixoTBEHQHpoX/mq3\n+M+cgcxMFzUmCIKgPTQt/FZXT3VZ/EYj5OaK8AuCUKPRtPA7a/G7bHK3uIixl63CEwRBqAiaF35n\nfPwuC+c8cwbCwsTiFwShRqN54S+2+B1l6HSZxX/mDLRrB0VFkJPjggYFQXBEdnY2ISEhsgitmtG0\n8BtMJXz81TG5e/Ys1K+vir9Y/YJQJVq0aMG3337r8JjIyEhyc3O9Lq2xt6Np4TeajQT4Bqg+/uqY\n3D1z5orwi59fEKpEeWkEPFke0VbCuIqOx5vLO2pe+KvV4i8p/GLxC0KlsVUu8ciRI+j1elauXEmH\nDh0YMGAAf/zxB3q9HovFAkBiYiKvvvoqAwYMoGHDhtx7770UFRUxYcIEmjZtyv33389ff/1lt9/D\nhw/zxBNPcN111/Hggw9a0zoXt/3KK68wcOBA6taty2+//YZer+ejjz4iNjbWmlM/NTWVESNGEBUV\nxYwZMzh37py1DVvHv/baa3Tp0oXQ0FA6depERkaGO/6kLsUrhD/QN7B6Lf7oaBF+QagCjsolrlix\ngtWrV/P111/bfCJYuHAhzz77LGlpaWzZsoXY2Fh69uzJ3r17OXv2LO+//77NPs1mMz179iQmJoZf\nf/2VPn36cPPNN5c6Zu7cuUyePJkzZ84QGRkJwPvvv8+SJUvIyMjg999/5/bbb2f06NGkpqby888/\nM2XKlFJtlDw+IyODJUuW8NVXX3HhwgVWrlxJ/fr1q/rnczuaFn5nUza4dHJXLH5BcCuPPvoorVq1\nIiAgoMw+nU5HcnIy/fr1o1WrVtx0000EBARwzz33UL9+fUaOHMmmTZtstvvtt9/SuXNn7rvvPkJC\nQrjnnnto0KABO3futB4zaNAgbr31Vnx9ffHz8wPgwQcftKZX/t///sctt9zC3/72N5o1a8b06dNZ\nu3at9Ynk6uPNZjOFhYVkZWVhsVho27YtjRs3dvFfzPX4ln+I5/BIOGf9+tCyJZw6Bfn5EBzsgoYF\nwXPoXnTNxKnygmsib+Lj4x3uL1nWMDw8nE6dOlk/N2rUyFoC8Wo2btzIli1bCAsLs24zmUxs3ryZ\nuLg4dDqdzb5Lbtu+fTsJCQnWz61bt8ZkMpGRkUHHjh3LHN+pUydefvllnnrqKY4cOcKDDz7IU089\nRbDGdUPzwl/tFn+9euqdpHVrOHAAYmNd0LAgeA5XCXZFsVcu8eoSia4iKSmJPXv2sG7dOrvH2Oq7\n5LZevXqRlpZm/ZyVlYWPjw/t27e328aYMWMYM2YMf/zxByNHjiQ8PJyHH364KpfidjTt6nE2LbPL\nwzlB/PyCUEUqUy4RypZrdJabbrqJvXv38tFHH3Hu3DkKCwtJTU0t9YRQXnvJycmsX7+ezz//nOPH\nj/PCCy8wZMgQ9HrbUrlr1y7S0tIoKioiKCgIX19fQkJCnB6zp9C08Fvj+J2w+F3q6gHx8wtCFbFV\nLtFWvL69EofF/3a0vyQ+Pj6kpqZy4MABunbtSmRkJP/5z39KiX15bbVs2ZKVK1eydOlS+vXrR6dO\nnUqVerz6+IsXL/L3v/+devXq0b9/f7p3787dd99tc3xaQtOlFzvN78Sy4cv4Yv8XGMwGpidNt3ne\n2LEwYADcc08VOlcUCAiAixchMBA++QQ++wxWrqxCo4LgXqT04rVFtZRePHr0KP3796d9+/YkJiay\nYsUKAHJzc0lOTiYyMpJhw4aRl5dnPWfOnDlERUURExPD1q1brdszMzO54YYbaNWqFc8884xTg6tI\nBa4qW/x5eeDnp4o+iKtHEIQai0Ph9/PzY9asWWRkZPDpp5/y7LPPkpuby/z584mMjCQrK4uIiAgW\nLFgAwKlTp5g3bx6bNm1i/vz5TJo0ydrWY489xpNPPsnOnTv5/vvv2bVrV7mDq9bJ3ZJuHoC2beHw\nYanpKAhCjcOh8Ddu3JguXboA0KBBA9q3b8/OnTvZsWMH48ePJyAggHHjxllnwdPS0hg0aBCRkZH0\n69cPRVGsTwMHDhxg5MiR1K9fn+HDh5eaObeHs/n4XRLOebXwBwVBkybw++9VbFgQBEFbOD25e+jQ\nITIyMujevTs7d+60LleOjo5mx44dgCr87dq1s57Ttm1b0tLSOHToEI0aNbJuj4mJ4ccffyy3T2sc\nvycsfpAJXkEQaiROBdTm5uYycuRIZs2aRe3atSs0iWBrBt7R+dOmTbP+O//PfKvF7ygts0vCOc+e\nVWP4S1Ls5x86tIqNC4IguI7U1FRSU1MrfX65wl9UVMSIESMYO3YsycnJAMTFxZGZmUlsbCyZmZnE\nxcUB6oq2jRs3Ws/dv38/cXFxhISEcPLkSev2ffv2lVodV5KSwv/6y69f8fG7e3LXnsVfYoJaEARB\nCyQmJpKYmGj9/OKLL1bofIeuHkVRGD9+PB06dOCRRx6xbo+Pj2fx4sUUFBSwePFiq4h3796d9evX\nk52dTWpqKnq93rqYITo6mk8++YScnBxWrVpV7rJtRVGczsfvVlePpGcWBKGG4VD4t23bxrJly/j2\n22+JjY0lNjaWr7/+mpSUFLKzs2nbti3Hjx9nwoQJgJpXIyUlhaSkJCZOnMjs2bOtbc2cOZMZM2YQ\nFxdHnz596Natm8OBmRUzep0eH71Pudk53TK5C1dcPRInLQhCDcKhq6d3796lstKV5IsvvrC5ffLk\nyUyePLnM9piYGH7++WenB1YcyglU3+Ru166lt9Wvry7q+usvaNq0ih0IgqAVpk2bxuHDh1m6dKmn\nh1KGxMRExo4dy/jx493Wh2ZTNhS7eQDPhHMWI5E9glDj0HKpR1tpKlyNZoXfIxa/PeEXP78gXLN4\nc4lFe2ha+AN81UIN1ZKds2RmzpJI6gZBqBSLFi2iR48ehIaGEh0dXarw+rp16xg6dCht27Zl1qxZ\npdK+HDhwgMcff5yIiAgiIyP58MMPAbh06RLz58+nY8eODBw4kDVr1ljPWbJkCb179+all16iWbNm\nDBo0iB9++MG6//Tp0zzxxBM0btyYESNGlOrvalJTU60ZCdq0acP48eMxmUwsW7aMhIQEevTowfLl\ny603hCVLltCnT59Sbej1en777TcA7rvvPh599FFryuYHHnigVMbQjIwM7rzzTpo2bWpNZ+Pu/Eua\nFv6KWPwucfVcHccP4uoRhEqQk5PDtGnT+Oijj7hw4QIbNmygRYsWAKxevZqnnnqKf/3rX6SmpvLj\njz/y6quvAqp13atXLxo2bMjevXvZvXu3NXvAjBkzWLlyJZ999hlPPfUUkyZNKhXLXlxp69dffyUh\nIYEnnnjCui8lJYW//vqL3bt3M3ToUObPn+/QnXLy5El27tzJ5s2bWbhwIUuXLuWNN97gnXfeYe7c\nucycOZNly5Y5/fdYtGgRo0aNIiMjg5MnT1rT3CiKQlJSEvHx8ezZs4fCwkK2b9/ufleUoiFKDifj\nVIYS/U60oiiKcvzicaXxzMZ2z5s1S1EmT65Cx0VFiuLjoygmU9l9f/yhKE2aVKFxQXAfGvsJW8nJ\nyVHq16+vrF27VjEajaX23XXXXcry5cutn9PT05WYmBhFURTlq6++Ujp37myzzc6dOyvr16+3fn7m\nmWeUSZMmKYqiKB988IESFhammM1mRVEU5c8//1T8/PyUvLw8paioSAkNDVUOHz5sPbdPnz7K2LFj\nbfbz3XffKTqdTsnOzrZuGzp0qLJw4ULr5/fee08ZOnSote/evXuXakOn01n7u/fee5UhQ4ZY9338\n8cdKfHy8oiiKkpaWpjRv3ty6Lz8/XwkICFAWLVpkc2z2/r8r+j3QtMUf4HPZ1VOOxV/lyd1z5yA0\nVG3oaiIi1FTNFy5UoQNB8CA6nWteFaB+/fosXbqUWbNm0aRJEx555BFOnz4NqCUSU1JSCAsLIyws\njP79+3PkyBFOnjzJd999R8+ePcu0l5uby549e+haIvKua9eubNmyxfq5ffv21oIpTZo0wWQycfLk\nSTIzM7FYLLRq1cp67A033ODQnRIeHk7z5s2tn7dv3+6wb0fodDrrUwuoOdCKXT1paWmlSk0GBQVZ\n0+G4E00Lv7NRPVWe3LU3sQug16uZOmWCV/BWFMU1rwpyyy23sHHjRvbt28fvv//OjBkzALVE4nvv\nvce5c+esr0uXLhEeHk5SUhLbtm0r01ZISAidOnUqldV3165d9O3bt9xxREdHo9frOXz4sHXbTz/9\n5NCdcnV5xV69etntu1mzZqUyE6Snp5dpz95Npnv37vzyyy/WzwUFBeyvBq3xDuF3t8XvSPhB/PyC\nUEEOHjzIt99+i8FgwN/fn4CAAOsq/rFjxzJjxgy2bt2K2Wzm9OnTrF69GlDLJ/7555/MnDmTs2fP\ncubMGaswJicn88Ybb3Dw4EFSU1P5+OOPGTZsWLlj8fPz46abbuLFF1/kxIkTLFu2jN27d1foepKT\nk1mwYAE//fQT6enpLFiwwNp37969OXr0KBs2bODo0aPWG1wxjp4s4uLiMBgMzJo1i9OnT/Pcc89V\nS2EdzQp/yTh+X70vFsWCyWLbrHerxQ8i/IJQQQwGA08//TQNGzakW7du1K1blylTpgDqk8BLL73E\nO++8Q8OGDenRo4c1w6+vry9btmzh+PHjtG/fntjYWPbs2QPAE088wbBhwxg+fDgvv/wyb775Jv36\n9QPKL9E4b948GjVqRJcuXVi1ahUpKSkOx391W3fffTdTpkxh4sSJpKSk8MgjjzBmzBhAdc+8++67\nTJ06lZtvvplRo0Y5XT5Sr9ezceNGtm3bRufOnQkICKBXr17O/ZGrgGZLL67LWsecHXNYN2YdAEEv\nB5HzeA61/GuVOe/DD2HTJvjoo0p2vGQJfPed2pAtPvtM3XfZKhEErSClF68tqqX0oicp6eoBHGbo\nrHI4Z3kWf6tW8McfVehAEARBO3iP8DvI0Ol2V0+TJmq+HkEQhBqAZoW/uOxiMY4sfpdM7tpavFVM\nw4ZqyGdRURU6EQRB0AaaFf6ScfyAmprZUxa/jw80aACnTlWhE0EQBG3gVOlFT2DT1eNOi9+R8MMV\nd0+zZlXoSBBcS1hYmKYzTQquJSwszCXteI/wO4jld7vFD+LnFzTJ2bNnPT0EwQvRrKunZBw/lG/x\ni/ALgiA4h2aFv2RaZijf4q+Sq8deSuaSiPALglBD0K7wW4z46/3h4EEoLHRo8VfJ1ZOfr74HBzs+\nToRfEIQagnaFv9jH/49/wLffOrT4qzS564ybB0T4BUGoMWhW+K0+/kuX4NIl91n85cXwF9O4MZw4\nUclOBEEQtINmhd/q48/Ph/x8sfgFQRBchKaF39/HHwoKoKDAvRa/M8LfuDGcPFmpvOSCIAhaQrPC\nb03ZUFBgtfgLTYU2j61SOKezwh8YCLVqqccLgiB4MZoV/lIWf35+uUnaKu3qcSaUsxhx9wiCUAPQ\ntPAH+Fzl4/ekqwdE+AVBqBFoWvj99X5XfPzuys4pwi8IwjWGZoXfYDIQYNGDxeKUq6faLH4J6RQE\nwcvRrPAbzUaCisU8P19Ny+wui9+ZOH5QI3vE4hcEwcvRtPAHGi3qh3Li+MXHLwiC4DyaFv6Aossx\n88WuHndk5xThFwThGkOzwm8wGwgouuy/KZ7cdXU4p8UCFy6As8UNRPgFQagBaFb4jWYjAUbnLP5K\nu3rOn4eQELUBZxDhFwShBqBp4fc3miE01H25eiri5gH1JqEokJtbic4EQRC0gbaF31CkCrO7LP6K\nCr9OJ1a/IAhej2aF32Ay4GcwqcJcjo+/2ix+kFh+QRC8HofCP27cOMLDw+nYsaN127Rp04iIiCA2\nNpbY2FjWrVtn3TdnzhyioqKIiYlh69at1u2ZmZnccMMNtGrVimeeecapgRnNRvyMZqcsfv3lq7BY\nnGr6ChWJ4S9GLH5BELwch8J///338/XXX5faptPpePTRR0lPTyc9PZ1bbrkFgFOnTjFv3jw2bdrE\n/PnzmTRpkvWcxx57jCeffJKdO3fy/fffs2vXrnIHZjQb8Svp6nFg8UMlQzorY/HLIi5BELwch8Lf\np08fwmyEOio2ctKnpaUxaNAgIiMj6devH4qikJeXB8CBAwcYOXIk9evXZ/jw4aSlpTkclKIoGM1G\nfAoN6uSuohCg6O2mZYZKhnRWJDNnMWLxC4Lg5VTKx//222+TkJDA66+/Tu7lCJcdO3bQrl076zFt\n27YlLS2NQ4cO0ahRI+v2mJgYfvzxR4ftmywm9Do9+kIDBAVBUBCBRYpdVw9UcoK3sj5+EX5BELwY\nJwPYr5CSksLzzz/PxYsXefzxx1m4cCFTp061+RSg0+nKbLN1XEmmTZtGkbkI/XY9qZH7SGzeHIKD\nCTRaynX1VNjiF+EXBMELSU1NJTU1tdLnV1j4i6330NBQHn74YSZOnMjUqVOJj49n48aN1uP2799P\nXFwcISEhnDx50rp93759JCQk2G1/2rRpnCs4x7w580g831i1+IODCTCYxeIXBEEAEhMTSUxMtH5+\n8cUXK3R+hV09f10WPZPJxIoVKxg8eDAA3bt3Z/369WRnZ5OamoperyckJASA6OhoPvnkE3Jycli1\nahXx8fEO+7BW38rPtwq/pix+CecUBMGLcWjxjx49mu+//56cnByaN2/Oiy++SGpqKrt378bf35++\nffuSkpICQHh4OCkpKSQlJeHv78/ChQut7cycOZO7776bp59+mlGjRtGtWzeHgypVb/ey8PsbTNqw\n+OvXV1fuGo3g71/BDgVBEDyPQ+H/+OOPy2wbN26c3eMnT57M5MmTy2yPiYnh559/dnpQ1rKLBQUQ\nHAxBQfgaijBbzFgUC3pd2QeVSgt/ReP49Xpo1Ei1+iMjK9ihIAiC59Hkyt1ShdYvW/y6ggL8ffxd\nt3q3sBCKiqB27YoPUGL5BUHwYrQt/CV8/C7P11Mcw28j8qhcZIJXEAQvRpPCbzCV9fG7PENnZfz7\nxYjwC4LgxWhS+I1mIwG+pX38FBS41uIX4RcE4RpFs8IvFr8gCIJ70KTwW8M53e3jr2hETzESyy8I\nghejSeGvjMVfYeG/cAHq1q3cAMXiFwTBi9Gs8F8dx1/s47eXobPCrp6LF6FOncoNUIRfEAQvRrPC\nb9fid5Wr58KFygt/o0Zw+nQly34JgiB4Fk0Kv8FkwF/vpy6yCgws7eN31eTuxYtqrv/K4O+vnpuT\nU7nzBUEQPIgmhd9oNlLb4qsKrF7vHou/Kq4eEHePIAhei2aFv5b5suCDeyz+qrh6QIRfEASvRbPC\nH2zSqf59sE7uBvoGutbir6yrB0T4BUHwWjQp/AazobTwuyOc0xWuHonlFwTBC9Gk8BvNRmoV2RB+\nBwu4xNUjCILgHJoV/iATZX38rrb4xdUjCMI1iHaFv4gyPn6XWfwmExgMV24slUGEXxAEL0WTwm8w\nGwgsKfyuDue8eBFCQiqXi78YKcYiCIKXoknhN5qNBBmVCk3uVsjir6qbB65Y/IpStXYEQRCqGc0K\nf2CR5YrwBwRAURGBOj/XWfxVmdgFqFUL/PzUSWJBEAQvQrPCH1BkueKD1+kgOJhaJr1rJnddIfwg\nIZ2CIHglmhR+g8mAv7GExQ8QFESwCddk56xqKGcxMsErCIIXoknhN5qNBBjNpYU/OJjgIlzn6qmq\njx9U4f/zz6q3IwiCUI1oVvj9jSbbwu+qyV1XWPzNmonwC4LgdWhW+P0M5tJx9sHBBBkV11j8rnL1\nNG0qwi8IgtehSeE3mA34GYvK+PgDixTthHOCCL8gCF6JJoXfaDbiV1hUxtUTaLRoJ5wTVFfP8eNV\nb0cQBKEa0azw+xrsCL8rwjnF1SMIwjWMJoXfYDLgYygq4+P3N5pdk6vH1a4eWb0rCIIXoUnhN5qN\n+BQay1j8/gazthZwFReCP3u26m0JgiBUExoWfkOZyV1/Y5FrLH5XuXpA3D2CIHgdmhV+/dXCHxyM\nX6HJdRa/K1w9IBO8giB4HZoUfoPZoAr/VT5+P4PRdT5+sfgFQbhG0aTwqxZ/YRmL37ewSFtRPaAK\nv1j8giBvqQgdAAAgAElEQVR4EZoTfkVRMJqNUFBYxsfvU2jf4nda+I1GKLoqVLQqSNoGQRC8DM0J\nv8liwlfviy4/v4zF71NYiNFsRLERPum0qyc3V/XvV6X6VknE4hcEwctwKPzjxo0jPDycjh07Wrfl\n5uaSnJxMZGQkw4YNIy8vz7pvzpw5REVFERMTw9atW63bMzMzueGGG2jVqhXPPPOMwwEZzAb8ffyh\noKCMj1+XX4Cf3nYxFqctflf690EsfkEQvA6Hwn///ffz9ddfl9o2f/58IiMjycrKIiIiggULFgBw\n6tQp5s2bx6ZNm5g/fz6TJk2ynvPYY4/x5JNPsnPnTr7//nt27dplt0+j2UiQzl9VcX//KzuKyy/6\n2i6/6LTF70r/PsjkriAIXodD4e/Tpw9hYWGltu3YsYPx48cTEBDAuHHjSEtLAyAtLY1BgwYRGRlJ\nv379UBTF+jRw4MABRo4cSf369Rk+fLj1HFsYzUZCFH8IDCztjgkOhoICuwXXK2TxuyqUEyA8HHJy\nKjCzLAiC4Fkq7OPfuXMn0dHRAERHR7Njxw5AFf527dpZj2vbti1paWkcOnSIRo0aWbfHxMTw448/\n2m3faDYSavYtO/kaFOTQ4vf1ddLid7Wrx9cXGjSAkydd16YgCIIb8a3oCbYmVu2hszGBWt75M1+Z\nybkdl5hWpJCYmkpiYqK6o9jVY8fi9/Fx0uh2tasHrkzwNmvm2nYFQRBskJqaSmpqaqXPr7Dwx8XF\nkZmZSWxsLJmZmcTFxQEQHx/Pxo0brcft37+fuLg4QkJCOFnCGt63bx8JCQl223/osYc4PPcrpv3p\nA8WiD1bhD/QNtGvxe8TVAzLBKwhCtZKYmHjFKAZefPHFCp1fYVdPfHw8ixcvpqCggMWLF1tFvHv3\n7qxfv57s7GxSU1PR6/WEhIQAqkvok08+IScnh1WrVhEfH2+3faPZSG1brp5iH7+vfYvfI64ekAle\nQRC8CofCP3r0aHr27MnBgwdp3rw5H3zwASkpKWRnZ9O2bVuOHz/OhAkTAAgPDyclJYWkpCQmTpzI\n7Nmzre3MnDmTGTNmEBcXR58+fejWrZvdPo1mIyEmvX0fv499H7/HXT2CIAhegENXz8cff2xz+xdf\nfGFz++TJk5k8eXKZ7TExMfz8889ODchgNlDb7GPb4i+e3K2qxd+0qVNjcZpmzWDLFte2KQiC4CY0\nt3LXaDZSy+xTevEWqOGdBgOBOv+q+/jF4q80aw+uZeDSgZzJP+PpoQiCUEm0KfwmXVmLX6+HwEBq\nK/ZX7npkARdcM5O7uYZcUr5MITQwlL5L+nL84rVxsxOEmob3CD9AUBB1TD52V+561OK/BoT/hdQX\nuLHljay8YyX3dLqH3h/0JutMlqeHJQhCBalwOKe7MZgMBBfZEf7gYELMPtpauQtQrx7k56v5hVyV\n9VNjpP+VzvK9y/k15VcAnuz9JPWC6tFvST++GvMVXRp38fAIBUFwFk1a/EFmXVkfP0BwMLVNem3l\n6gE1tUQNtvrNFjMTvpzAK0mv0LBWQ+v2B7s+yOxBsxm4dCBbs7c6aEEQBC2hTeEvwq7FX9ukp9BU\nWGaXRyd3oUZP8L7707v46f24P/b+MvvuaH8Hy4Yv4/b/u53dJ3Z7YHSCIFQUzQm/wWwgqEixK/y1\nTPqqh3O62tUDNXaC90TeCZ5PfZ4Fty1Ar7P9dRl4/UCmJEzh3Z/erebRCYJQGTQn/EazkUB7wh8U\nRLBJV/lwToMBLBYICHDNYEvixa4ei2LhouGizX2PbXiM8bHj6dCog8M27up4Fyv3raTIXOSOIQqC\n4EI0N7lrFX47Pv5aRXC8suGcxW4eV1XfKokXu3rGfTGO5XuX06Z+G3pE9KBHRA96Nu/J0YtH2X50\nu3VC1xEt6ragbf22bDi8gVvb3FoNoxYEobJoUvgDjBa7rp4go1L57JzucvOA6upxcnWylliXtY7N\nf2wm5/Ecfjv3G9uPbue7I9/x8paXyb6QzRejvqCWfy2n2rqr412s+HWFCL8gaBzNCb/BbCDQaLYv\n/Ka8yrt63DWxC15p8V80XGTClxNYNHQRoYGhxDaJJbZJLA93fxiAgqICgvycD0+9I+YO/rXpX+QZ\n86jtX9tdwxYEoYpo0sfvb8/iDwoiqMi+xV+uq8cdoZzFeOHk7lMbn2JAqwHc1Oomm/srIvoADWs1\npFdkL1YfWO2K4QmC4CY0Kvxmuz7+QINSNYvfXa6eJk1U4a9AoRpPsvmPzaw+sJqZA2e6tN27OtzF\n8r3LXdqmIAiuRZPC72c02XX1BBaZK2/xu9PVU7s2+PnB+fPuad+FFBQV8MDqB5g7eC51A+u6tO3k\n6GS2ZW/j9KXTLm1XEATXoTnhN5gM+BnsC7+/0VJ5i9+drh7wGnfPC6kvcEOTG0iOTnZ527X9azM4\najCf7vvU5W0LguAaNCf8RrMRP0ORfeE3mCqfndOdFj94xQTvrj938eEvHzLnljlu6+OujuLuEQQt\n413CHxSkCn9ls3O608cPml/EZVEsjF89njcHvkmjWo3c1s/N19/MgTMHOHL+iNv6EASh8mhS+H0M\nRXYnd/0NRQ59/A7nVq9xV8/O4zspMhdxV8e73NqPn48ff4v5Gx/vtV3BTRAEz6I54TeYDfgUGu26\nenwLi2xa/DqdWqvFYnHQ+DXu6llzcA1D2w5F546Vy1cxpuMYVvy6wu39CIJQcTQn/EazEZ9Cg33h\nt2PxgxMTvO529Wjc4l9zcA1D2gyplr56Nu/JRcNF9pzcUy39CYLgPJoTfrOhEB06NTTyaoKC8C00\n2kzLDE6EdLrb1aNhiz/7QjbHLx4nISKhWvrT6/SM7jCaFXvF6hcEraE54dcVFmIOspM9MzgYn0KD\nTVcPOGnxu1v4NWrxrz24lsFRg/HR+1Rbn3d3uptle5bZ/f8SBMEzaE/4CwpRAgNt7wwOxqfA4NDV\n49Did7erp0kTOHXKycIA1cvag2u5rc1t1dpnh0Yd6NK4C+/9/F619isIgmM0J/w+BQaUIPvCry8s\ntGtBlhvS6W5Xj58fhIWp4q8hLhkvsSV7Czdff3O19/1S/5d4Zcsr5BflV3vfgiDYRnPCrzcYHFr8\nuoLCyln8iuJ+Vw9ocoJ3428b6d6sO6GBbnzascMNTW6gZ/OezNs5r9r7FgTBNpoTfp9Co+0YfoCg\nIHT5BZWz+A0GNd7THdW3SqLBCd7qjOaxxUv9X+KN7W+Qa8it8LmffAJvvumGQQnCNYz2hN9gJ4Yf\n1O0FBRhNBhQbK7UcTu5Wh7UPmpvgtSgWvsz6str9+yWJaRjDwOsH8taPbzl9jqLAG2/Ak0/CrFmw\nbp0bBygI1xiaE37fwiJ09ix+Hx90/v4EW3wospSt7eownNPd/v1imjXTlMX/058/UTewLq3rtfbo\nOF7o9wKz02ZztuBsucdaLDBlCnz0EWzbBsuWwbhxmrqfCoJXoz3hNxShD7Ij/ADBwYRZAiqeofMa\ntfjXZlV/NI8tWtdrze3tbuc/P/zH4XEGA4weDenpsGULRERAv34wYQKMHavJgClB8Do0J/x+hiJ0\ntRyU7QsKIlTxr3iGTneHchajscndNQc8698vyXN9n2PBrgWcumQ76unCBbjlFvX/cP16qFuiVMCz\nz6rbX3utmgYrCDUY7Qm/0YzenqsHIDiYuma/imforC5Xj4Ymd49dPMYfF/6gZ/Oenh4KAJGhkYzp\nOIbXt71eZp/BAP37Q0wM/N//wdWBXT4+sHw5vP02bN1aTQMWhBqK5oQ/xKRHV46rp47Fr3IWf3UI\nf/PmkJ2tiRKMXx78kkGtB+Gr9/X0UKw83ftpluxewvGLpW+Oc+eq98y331ZF3hbNmsGiRTBmDJwt\nf6pAEAQ7aE74a5t97Uf1gCr8Zt+KW/zV5eqpX18tw3jkiPv7Koc1B9dwW5Tn/fslaRLShHs638P8\nXfOt23Jy4NVXYeZMNcuqI269FUaMUCd7NXBvFQSvRIPCr7cfxw8QFESI2ceuxe9xVw/ADTfAzz9X\nT192yC/K5/s/vmdQ60EeHYct7oi5gzUH11g/v/gijBoF0dHOnf/aa3D4MHz1lZsGKAg1HM0Jfy2z\nvlyLP8TkazNDp8Nwzupy9YAmhH/Tb5vo2qQrYUFhHh2HLeKbxfNn7p/8cf4P9u9XF2lNm+b8+f7+\nMHkyvCcpgAShUlRa+Fu0aEGnTp2IjY2le/fuAOTm5pKcnExkZCTDhg0jLy/PevycOXOIiooiJiaG\nrQ5m52qbfMoXfrO+cuGc1eHqAejaFX76qXr6soOnV+s6wkfvw+CowXyZ9SWPPw5PP616yCrCqFHw\n/feaCqASBK+h0sKv0+lITU0lPT2dHTt2ADB//nwiIyPJysoiIiKCBQsWAHDq1CnmzZvHpk2bmD9/\nPpMmTbLbbrBJV67w1zLpKz656wlXj4ec0GaLmS8OfMGw6GEe6d8ZhrQZwpLta8jMhIcfrvj5tWvD\nnXfCBx+4fmyCUNOpkqvn6rQJO3bsYPz48QQEBDBu3DjS0tIASEtLY9CgQURGRtKvXz8URSE313be\nlmCTzrGPv1j4KzO5W13C37SpOkvpobDO7Ue307h2Y66vd71H+neGG1sM5KfT23jx1bxKp0968EE1\nysdhuU1BEMpQJYs/KSmJYcOGsXr1agB27txJ9OUZuujoaOuTQFpaGu3atbOe27ZtW+u+qynX4g8K\nolaRTrvhnKCKvgfdPav2r2J49HCP9O0sn62oQ52L8QR32FjpNrp2Vb13mza5cGCCcA1Q6QDvbdu2\n0aRJEzIzMxkyZAjdu3e3mTjNHvYKfn98vJBtq1ZBRgaJiYkkJiaWPiA4mFomKreAq7p8/HDF3ZOc\nXH19oj6FfZ75OatHr67WfitCbi48/zyMm30bX2at5fZ2lXNJ6XTw97/Du+/CgAEuHqQgaJjU1FRS\nU1MrfX6lhb9JkyYAtGvXjqFDh7JmzRri4uLIzMwkNjaWzMxM4uLiAIiPj2fjxiuW3f79+637rubh\nsAC6jh8PN95ou+PgYIJMOnIrGs5ZnRY/qMLvAQd0+ol0fPW+dGzUsdr7dpZZs+Cmm2DijbfR+4PX\nsCgW9LrKPXzedRf8619q7ZtGjVw8UEHQKFcbxS+++GKFzq/Ury0/P9/qoz99+jTr169n0KBBxMfH\ns3jxYgoKCli8eDEJCWph7+7du7N+/Xqys7NJTU1Fr9cTEhJis+2gIqXcyd1go2I3qkcTrh5Qhd8D\nrp5V+1cxvN1wu09UnqaoCBYsUNMtX1/vesICw/jpz8r/nUJDYdgw+PBDFw5SEGo4lRL+kydP0qdP\nH7p06cKoUaN47LHHaN68OSkpKWRnZ9O2bVuOHz/OhAkTAAgPDyclJYWkpCQmTpzI7Nmz7bYdWKSU\nu4ArsEix6eO36+qprupbJbnuOjUBzV9/VV+fwOeZn3N79O3V2mdFWLsWrr8e2rdXPw9pO6TUYq7K\n8Pe/qzH9spJXEJyjUq6eli1bsnv37jLbQ0JC+OKLL2yeM3nyZCZPnlxu2wFOWPyBRkvFLP6CAnWn\nv3+5/bsMnU61+tPT1SLs1cCBnAOcKzhHfER8tfRXGRYuhIceuvJ5SJshTFo3iZf6v1TpNhMS1P/a\n77+Hq6eEBEEoi+ZW7gYYzc4Jf0Us/uq29oupZnfPqv2ruL3d7ZX2l7ub335T/xx/+9uVbQkRCWRf\nyObYxWOVblenU0M7ZSWvIDiH5hQiwGgpV/gDjOaKhXNW56rdknTtWq2pG7Tu5nnvPbWYSsmUy756\nX26JuoW1B9dWqe2xY+HLLyVrpyA4g+aE399oLncBl7/BXLFwzupctVuSaszZc/TCUQ6fO0y/6/pV\nS38VxWhUg5xKunmKuS3qtioLf716cNttsHRplZpxKRYL7N2rppweNUotMqORUg3CNY7mhN/PaCp3\nAZe/wVSx7JyecvW0agXnz6t5h93M//b/jyFthuDn4+f2virDF19Au3bQtm3ZfTe3vpnNf2wmvyi/\nSn1MnKhm7ty+vUrNVImsLHjzTRg6FBo0gOHD1Xv/oEHQowf07Am//uq58QkCaFD4Lb4+oHcwrOBg\n/A1FNrNzOhR+T7h69HqIja0Wq//z/dp281w9qVuSuoF16da0G5t+q9oS3J491RQOt9+uLuqqDsxm\ntSD8k0+qN7Z+/WD/ftX1lJGh3ggWLYL77lMXrb3yCiQlwXffVc/4BMEWmhP+Iv9yAo2Cg/E1mGy6\neho3hmO25gg95eqBavHzn750mp//+pmB1w90az+VJStLdXnc7uC+NKRN1cM6AQYPVkszzp6thnka\nyn5NXMLevfDPf6oBWxMnqlFFH32kfv/efRfuuMN2MNeYMWoa6pEj1VKSguAJtCf8geW4KoKD8Ss0\n2nT1dO4Mv/xi4xxPuXqgWvz8aw6uYUCrAQT5OXCReZB334V778VhMrbb2qh+fpPF3tJr54mKgh9/\nhDNn1Dq+rkrdnJ8PS5aoLptbblHnFdLS1O/cv/8NcXGOH1aLSUqCb79V01G/+qqsPxCqH80Jvymg\nHOEPCsKn0GjT4u/SBXbvtvFD8rTwuzmk8/PMzxneTptJ2QwGdVXt3//u+Lio+lHENonl0fWPuqTf\nkBBYuVIt1di9O7zzDhw8WDmRPXECHnlELaf86adqiogjR9TKYS1bVm58HTrADz+oheX791ff3fV0\nIghXoznhNweUs8gqKAjfQiMGGz7+4kfrMotlqztBW0natFETyZw755bmLxousvmPzdwadatb2q8q\nn3+uPom1bl3+scuHL2fD4Q2895NrAvL1enjmGVi2TH3oSkpShfrBB9WbQnmhn7m58MIL6ipjHx91\nLd7atTBkiDqfVFWaNVOfGFJS1DmQyEh44gnVNSYI7sT7hN/PD8VHj9lQVvh1OjvuHk9a/D4+6qBs\nrHR2BSszVpLYIpHQQA/d2MphwQL7k7pXUzewLqtHr+aZb59hyx9bXDaGxERYvBiOHlXr9HbooLps\nWrRQnwaefVZd9Ws0qscXFakhmG3aXFl09p//qMLsagICVH//t9/Cli3qE0mvXupTwOzZchMQ3IP2\nhD+o/KoclqBAdAUFNvcVu3tK4UnhB7e6e977+T0evOFBt7RdVXbtUoWrIpmp29Rvw7Lhy7jz0zv5\n4/wfLh2PTgcxMWq93i+/hNOnYcYMNd5+6lRo2FBdC9C+vRp++tVX6rqAFi1cOgy7tGkDb7yh3qAm\nTVInkBMT1TmLyZPh66/FHSS4Bs0JvyXQSeHPty38No1rT7p6wG2RPXtO7uHYxWPcEnWLy9uuKoqi\nitf06eBXwaUFA68fyJO9nmToJ0PJM+aVf0IlCQhQhfWVV2DnTjh8WA3DXLAANmxQI3E9QUCAGgH1\n/vtqlNCnn6oRay++qIasnj/vmXEJNQfNCb/ihPArgYHoC8q6ekC1+DXl6gG3Rfa89/N7jIsdh6/e\nBQ5nF7N8ueoyue++yp0/OX4y3Zp2455V92BRqqe2YoMGqtslKalaunOKYvfl00+rC9N691YjiuxU\nLhUEp9Ce8AcFln9QcLBdi79tW8jOhkuXSmz0tPC3a6c+v7vw11pQVMCKvSsYHzveZW26itxcdUHT\nnDnOhTfaQqfTMW/wPE5eOsnLm1927QC9FJ0O3noLOnZUVwbnV22hs3ANo0HhLz8W3a92KHkXT9tc\n4u/np+rs3r0lNnpq5W4xvr7qjKILJ3g/3fcp3Zt157q617msTVfxyitqAbUeParWToBvAJ/e8Slz\nd87lx2M/umZwXo5OB/PnQ0SEmg5CfP5CZdCc8DvM03MZfa1atK/V0m7lpjLuHk+u3C2ma1fVkewi\ntDqpm5WlZuF87TXXtNckpAlzB89l7KqxXDJeKv+EawAfHzXhXe3aqmuqqMjTIxK8Da8UfoKC6Fwn\nyq4VWGqCNzcXCgs9L/xDh6oB5S5Yppl5OpODZw4ypM0QFwzMtTz6qBqL3rSp69ocETOCns17MvWb\nqa5r1Mvx9YUVK9TcVGPHOqg1LQg20Jzw6xylZC4mOJj2tVryw7EfbO4uFdL5/fdqiSZXrLipCgMH\nqi6ntLQqN/V++vvcH3u/5jJxrlunJihzotBahZkzaA7rstaxLmud6xv3Uvz91Yifc+fUHEBi+QvO\noj3hD3JO+NsGR/LDsR9QbFjQnTurqW/NZuCbb2DAANcPtKLo9eoSzblzq9SMwWRg6S9LeSD2ARcN\nzDUYjWpag7fecpyTp7KEBoayZNgSHljzADn57k9z7S0EBqprDgoKYMQI9eFWEMpDe8Jfq1b5BwUH\n04BgLIqFoxePltkdGqqG5h0+jHaEH+D++2HNGnXlUCX53/7/0TG8I9fXu96FA6s6b7+tpmW41Y2Z\nIxJbJDK6w2geWvuQzRv+tUpgoGr5BwZKtI/gHJoTfn2wE8IfFISuoIAeET3s+vm7dIEDm47ByZNq\nHL0WqFdPDcVYtKjSTbz787uam9Q9eVLNMvnmm+7va3rSdA6eOcjSPRoqtaUB/P1Vn3/jxmrRF4nz\nFxyhOeH3Ca5d/kHBwZCfT0JEgsMJXsNXG9XVOD4+Lh5lFXj4YXVpqM3iwI45fPYwe0/u1VzBlWee\nURdq2aqu5WoCfQNZdvsyHtvwGEfOH3F/h16Er6+agygmBm66yW15AYUaQI0V/i5doP7PGnLzFNO1\nK4SHq4lgKsj76e8ztvNYAnzd4ESvJD/9pGasfO656uuzc+POPNHzCcauGovZUvEbaE1Gr1fj/Hv1\nUtdSSPF5wRbaE/5aIeUfFBwMBQV0a9qNX07+Yjs3fycLHU5s1J7wg2r1V3CS90z+Gd7/+X0mdJ3g\npkFVHEVRI3imT6/+9XGP9XwMfx9/XtvqogUDNQidTs0mmpgIN9+sLmMRhJJ4r/Dn51PbvzZR9aLY\nfaLsitjrLu4llxBO165kpQx3cuedau6eCuTcnb5lOnfE3EFU/Sg3Dqxi/N//qdEk999f/X3rdXo+\nHPYhc3bMYcfxHdU/AI1TLP7x8eLzF8qiOeH3q+2E8AcFWUMX7Ll7dJs2sid8gO1SjJ4mMFBVy/nz\nnTr88NnDLP1lKS/0e8HNA3Oe/Hx1odbs2Z6bQomoE8E7t7zDmM/HuDWLp7ei06n5kjp1UmsRX5KF\nz8JltCf8tZxYYXvZ4gfUyJ7jNvz833xDTqxGhR9gwgS1OrcTsXdPb3qaKQlTCK8dXg0Dc44ZM1Q/\ncu/enh3HHe3voE9kHx75+hHPDkSjFPv8W7dWK4dJqKcAGhR+39pOCv/lQiw2Lf7CQti2jcBbktxV\n+KrqtGyprij++GOHh/1w9Ad+OPYDU3pMqaaBlU92tlrD9vXXPT0SldmDZpN6JJXP9n3m6aFoEr1e\nze3ftCkMGyZuH0GDwu9syoZi0yWqfhQXCi9wIu/Elf3bt0P79sT0rKtd4Ycrk7x2FiMpisLUb6by\n7/7/JtjPib9LNfHEE/DPf7qnFGFlCAkIYfnw5Uz8aiLHLh7z9HA0iY/PlXKTEREwahSsXn2l3KRw\nbaE54Xc2SVux8Ot1euIj4kk7ViIHzuXVujExcOiQhpex33yz+sv7739t7l61fxWXjJcY22lsNQ/M\nPl98AT/8AI8/7umRlCY+Ip5/dv8noz4dRa5BTFpb+PrCu++qdYT794eZM9WngIcegm3bPD06oTrx\nTuEvYfEDJDRLKO3nvyz8gYGqb3PfPjeM0xUUP4NPngynTpXaZTQbeXLjk7wx4A189NpYgPb11/Dg\ng7BypfpfoDWe7v00HRp14MaPbpR8Pg6oX18V+82b1eCy669XM3z+7W9qvSCh5qM94Q90ogJX06Zw\n4gRs2gRc5ec/c0YNk0xIAOzU4NUSCQlwzz3wj3+U2rxw10Ja12vNgOu1sQ5h40Z1mP/7H3Tv7unR\n2MZH78P8W+dzY6sb6ftBX3H7OEFkpOq6y8hQi8zHxqpPApLps2ajPeHX6co/pn59+Pxz1VG5dSvx\nEfHs+nMXJotJvRn06aMmL8FODV6t8dJLasmwlSsBOFtwlulbpjPjphkeHpjK99/DXXfBZ5+pxb61\njE6n49UbX+X+LvfTe3FvDp456OkheQVBQWox9x9+UG/ysbHqE4FQM9Ge8DtL375qVqrhw6m7N4vm\ndZrz66lfy2TjLJWbX6sEBsIHH2D+5z+Y/tkjtH2nLeNjx9MxvKOnR8a2bXDHHepirT59PD0a53m8\n1+M82/dZEpckkv5XuqeH4zVERal1FaZNU3P8jxjh0sJxglZQNESlhrN6taI0aqQ8PytZmb9jnqJc\nd52iZGRYd586pSh16ihKbq7rxulqdv+1W7ln1T3KnL4Byk99opRDZw55ekiKoijKDz8oSsOGirJ+\nvadHUnlWZqxUGs5oqLy+9XUl/a90xWwxe3pIXkNenqK89ZaiREYqSmKioqxbpygWi6dHJdiiotqp\nu3ySJtDpdJXLs75yJfkpD/Du37vxyIf74dixUi6jhx5SC1+tXq2dEMSjF46y5uAaVu5bycEzB/ln\n93/yUMw9hPXoDy+/rM60eYiTJ9WauR99BEuXqqs+vZm0Y2ks3bOUDYc3cMFwgQGtBjCg1QD6t+xP\nRJ0I9DrvffCtDoqK4JNP1EV7ej1MmaI+GYSEqHV/i19BQc55agXXU1HtrFbh37x5Mw899BAmk4lJ\nkybxz3/+s/RgKiv8QPY7r9Bk8jP4jb1XDVgugaKoueL/8x/VT92jR2WvoPIoisJ7n73H8QbHWX1g\nNUcvHGVw1GCGth3K0LZD8fdR5yTYvl19vt67V60mU42cOaP+uN97D+6+G55+Gpo0cf781NRUEhMT\n3TY+V/D7ud/55rdv+Oa3b9jyxxbOF56naUhTIupE0Dy0Oc3rNCeqXhSxTWJp37B9qUyoxddXUFTA\n7hO72fnnTgqKChjebrimcihVBmf+7xRFdQO9+64ahJaXp75yc9X3sDC1+PuYMWoSWi3dBLzhu1kV\nNLQZJGwAAAlvSURBVC38sbGxzJ49m+uuu46bb76ZrVu30qCEuFVF+M0WM/ffXZu3H15LaK8bbR7z\n1Vdq3vj//EcNX3MlBUUFnC88z/nC85zOP82hs4fIOptF1pksss5mcfjsYQK3BnL/lPsZ2mYoPZr3\nwFdvpw7w1KnqDNsDD6g3gYqobyU4dkz9Mc+dq/rzn3kGmjeveDvTpk1j2rRpLh+fOyk0FXLs4jGO\nXjiqvl88yv6c/aSfSOfw2cO0qd+G2CaxtGvQjv/O/S+Wfhb25+wnukE03Zt1R6/T83nm5zQNacqd\n7e/kjpg7NFcdzRlc8X+3f7+6EH35cvXJ4K67YPRoaNPG8zcBb/xuVoSKame1VSC/cDk3bN++fQEY\nOHAgaWlp3OqiWn0+eh+O3dyDjwOyuOnsdYQFhhEaGFpKXAcPhg2bjAy7M4/t+3KZMPkiFw3nOVd4\njvOF5zlXcI6Lhos0rNWQlnVb0qJuC5rXuQ6dORCTCYp0uWSc+YXdJ9JJP5HO7hO7+TP3T84XnkdB\noW5AGCF+dQn1q0+L0OuJqh/FrS3/RkxCFG0btuYt4ywe7zmNvDz443fVUsrPVx+T69ZVLabatUH3\n+uvqXWrlSjXRfefOakbPESPUXP6oqXZ//1195eSoxb0aNICGDdX3+vVtJ08zGiE9XY3eKH4VFEBy\nMuzapWaSuJYI9A2kdb3WtK7Xusy+gqIC9p7aS/pf6WTmZNKoViOev/V5Ood3JsjvynqTt295my3Z\nW/hvxn/pubgnEXUi6BzemSYhTWhcqzFNQprQpHYTGtZqSLBfMIG+gdaXr94Xi2IhvyifPGMeecY8\ncg25FJgKCPEPISwojLqBdanlVwudp9WzHKKj1cigadPU79Ly5epCsbw89XvVsqW6crhlS7VSmNGo\nfvcKC6+8fH3V38DVbiQ/P3VfyZdOp7qhTKYrL7NZPT4sTP1NhYRU/aZjMl0ZmzPR5t5AtVn8Gzdu\nZNGiRXx8OTfNggULOH78OP/+97+vDKYKFj/Akt1LeHvH26VEPNgvmJCAEAqKCsgz5qGgUNsvhIIL\nIZguheBbFIbeWBedse7l9zqY/E9hrPU75pAjUOcoFDRAZwpCqfUXnG6P7+lY/M/GEnwhFsv5SPLP\n1sV4KYhatdQvXWCg+oUs/jIXFKiPyRbLNIKCppX6UgcFqVkTz59XXwUFam77OnXUdkIDCulbuIEB\n5/9LjzNrmdhqPWtOxVNUdOXH1KCBWm0pJ0d9nT6ttuVr47ZusUC7dmpYZo8e6vv117vGIqvpVpUz\n12eymNiWvY2ss1n8lfsXJy6d4ETeCf7K/YvT+acpNBVSUFSgvpsK0KHDolgI8gsixD+E2v61qe1f\nm0DfQPKMeep3ufAcRrORuoF1SemWwkv9X/LItVWWc+fgyJErhsrvv6uuooAA9TseFKS+BwSowl3s\nOip+z8srLe7FL4tFvSGUvCno9erx586V/j2ZTNMIDp5Wpg2druxNRacDg+HKb9diUcc4cyakpLjl\nT1RlNOvqcVb4BUEQhIqjSVdPXFwcj5dI8JKRkcGgQYNKHaOhACNBEIQaS7XFsYVers23efNmjhw5\nwjfffEN8fHx1dS8IgiBcptosfoC33nqLhx56iKKiIiZNmlQqokcQBEGoHqp15Uq/fv3IzMzk0KFD\nTJo0ybp98+bNtGvXjqioKN5+++3qHJJbGDduHOHh4XTseCXlQm5uLsnJyURGRjJs2DDy8ry3VODR\no0fp378/7du3JzExkRUrVgA14xoLCwuJj4+nS5cuJCQkMGvWLKBmXFtJzGYzsbGxDBkyBKhZ19ei\nRQs6depEbGws3S9nFKxJ13fp0iXuvfde2rRpQ0xMDGlpaRW+Pk0sWZw8eTILFy5k48aNzJ07l5wc\n706pe//99/P111+X2jZ//nwiIyPJysoiIiKCBQsWeGh0VcfPz49Zs2aRkZHBp59+yrPPPktubm6N\nuMbAwEC+++47du/ezffff8+iRYvIysqqEddWktmzZxMTE2MNqKhJ16fT6UhNTSU9PZ0dO3YANev6\nXnjhBSIjI9mzZw979uwhOjq6wtfnceEvGd9/3XXXWeP7vZk+ffoQFhZWatuOHTsYP348AQEBjBs3\nzquvsXHjxnTp0gWABg0a0L59e3bu3FljrjH4crGBvLw8TCYTAQEBNebaAI4dO8ZXX33FAw88YA2o\nqEnXB2UDRWrS9W3cuJF//etfBAYG4uvrS2hoaIWvz+PCv3PnTqKjo62fY2Ji+PFHG8XTvZyS1xkd\nHW21RLydQ4cOkZGRQffu3WvMNVosFjp37kx4eDj/+Mc/iIyMrDHXBjBlyhTeeOMN9PorP/+adH06\nnY6kpCSGDRvG6tWrgZpzfceOHaOwsJCUlBTi4+N5/fXXKSgoqPD1eVz4rxVqYqhqbm4uI0eOZNas\nWdSuXbvGXKNer+eXX37h0KFDzJs3j/T09BpzbWvXrqVRo0bExsaWuqaacn0A27Zt45dffuHVV1/l\n0Ucf5cSJEzXm+goLCzl48CAjRowgNTWVjIwM/vvf/1b4+jwu/HFxcezfv9/6OSMjg4TL1bNqEnFx\ncWRmZgKQmZlJXFych0dUNYqKihgxYgRjx44lOTkZqHnX2KJFCwYPHkxaWlqNubbt27ezevVqWrZs\nyejRo/n2228ZO3Zsjbk+gCaXc1u1a9eOoUOHsmbNmhpzfa1bt6Zt27YMGTKEoKAgRo8ezddff13h\n6/O48F8r8f3x8fEsXryYgoICFi9e7NU3N0VRGD9+PB06dOCRRx6xbq8J15iTk8P58+cBOHPmDBs2\nbCA5OblGXBvAK6+8wtGjR/n999/55JNPSEpKYunSpTXm+vLz88nNzQXg9OnTrF+/nkGDBtWY6wOI\niooiLS0Ni8XCl19+yU033VTx66ts4n9XkpqaqkRHRyvXX3+9Mnv2bE8Pp8qMGjVKadKkieLv769E\nREQoixcvVi5evKgMHTpUad68uZKcnKzkarkyTDls2bJF0el0SufOnZUuXbooXbp0UdatW1cjrnHP\nnj1KbGys0qlTJ2XgwIHKhx9+qCiKUiOu7WpSU1OVIUOGKIpSc67vt99+Uzp37qx07txZSUpKUhYt\nWqQoSs25PkVRlAMHDijx8fFK586dlccee0zJy8ur8PVpqhCLIAiC4H487uoRBEEQqhcRfkEQhGsM\nEX5BEIRrDBF+QRCEawwRfkEQhGsMEX5BEIRrjP8Hfz+VGv1q1TsAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3456210>"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    }
   ],
   "metadata": {}
  }
 ]
}