CG_ClassifySimple.json

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import random as ra\n",
      "class ReservoirSampler(object):\n",
      "    ''' Simple reservoir sampler.  Add as many elements as you want,\n",
      "        but it will only keep up to k, sampled uniformly without\n",
      "        replacement. '''\n",
      "    \n",
      "    def __init__(self, k):\n",
      "        self.k = k\n",
      "        self.r = []\n",
      "        self.n = 0\n",
      "    \n",
      "    def add(self, obj):\n",
      "        if self.n < self.k:\n",
      "            self.r.append(obj)\n",
      "        else:\n",
      "            j = ra.randint(0, self.n)\n",
      "            if j < self.k:\n",
      "                self.r[j] = obj\n",
      "        self.n += 1\n",
      "    \n",
      "    def draw(self): return ra.choice(self.r)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Example of how reservoir sampler works\n",
      "ra.seed(725) # so you and I get same results\n",
      "example = ReservoirSampler(10) # sampler that will keep only 10 samples\n",
      "for i in xrange(1000): example.add(i) # add all numbers in [0, 1000)\n",
      "print example.r # print just the ones we sampled"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[250, 656, 388, 924, 930, 478, 476, 516, 525, 28]\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import re\n",
      "nonAcgt = re.compile('[^ACGTacgt]')\n",
      "def fastaKmerParser(fh, k):\n",
      "    ''' Parse k-mers from FASTA filehandle. '''\n",
      "    kmer, off = [], 0\n",
      "    for ln in fh:\n",
      "        if ln[0] == '>':\n",
      "            kmer, off = [], 0\n",
      "        else:\n",
      "            for c in filter(lambda x: x.isalpha(), ln):\n",
      "                if len(kmer) == k:\n",
      "                    kmer.pop(0)\n",
      "                kmer.append(c.upper())\n",
      "                off += 1\n",
      "                if len(kmer) == k:\n",
      "                    kmerstr = ''.join(kmer)\n",
      "                    if not nonAcgt.search(kmerstr):\n",
      "                        yield (kmerstr, off - k)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def fastaKmerParserIslands(fh, k, isles):\n",
      "    ''' Yield k-mers along with flag indicating whether k-mer lies\n",
      "        entirely within an island (\"island\") or not (\"non-island\") '''\n",
      "    curIsland = 0\n",
      "    for kmer, off in fastaKmerParser(fh, k):\n",
      "        while curIsland < len(isles) and off >= isles[curIsland][1]:\n",
      "            curIsland += 1\n",
      "        if curIsland < len(isles) and off >= isles[curIsland][0]:\n",
      "            if off + k <= isles[curIsland][1]:\n",
      "                yield kmer, \"island\"\n",
      "        else:\n",
      "            yield kmer, \"non-island\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import string\n",
      "def parseIslands(fh, chromosome='chr22'):\n",
      "    ''' Parse a file with island annotations, per the output from Hao Wu's\n",
      "        model-based approach: http://rafalab.jhsph.edu/CGI/.  Only take\n",
      "        records from given chromosome name. '''\n",
      "    isles = []\n",
      "    for ln in fh:\n",
      "        ch, st, en, _ = string.split(ln, '\\t', 3)\n",
      "        if ch == chromosome:\n",
      "            # convert from 1-based closed interval to 0-based right-open\n",
      "            isles.append((int(st)-1, int(en)))\n",
      "    return isles"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import urllib\n",
      "islandFh = urllib.urlopen('http://rafalab.jhsph.edu/CGI/model-based-cpg-islands-hg19.txt')\n",
      "islands = parseIslands(islandFh, chromosome='chr1') # this takes a few seconds\n",
      "islandFh.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "islands[1:10]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "[(28704, 29791),\n",
        " (135085, 135805),\n",
        " (136163, 137362),\n",
        " (137664, 138121),\n",
        " (228634, 228770),\n",
        " (326014, 326481),\n",
        " (326943, 327293),\n",
        " (327549, 328356),\n",
        " (436853, 438161)]"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import gzip\n",
      "def hg19ChrKmers(k, chromosome='chr22'):\n",
      "    # Get island annotations\n",
      "    islandFh = urllib.urlopen('http://rafalab.jhsph.edu/CGI/model-based-cpg-islands-hg19.txt')\n",
      "    islands = parseIslands(islandFh, chromosome=chromosome) # this takes a few seconds\n",
      "    islandFh.close()\n",
      "    # Get FASTA\n",
      "    fastaUrl = 'ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/%s.fa.gz' % chromosome\n",
      "    fastaFn, _ = urllib.urlretrieve(fastaUrl)\n",
      "    with gzip.open(fastaFn) as fastaFh:\n",
      "        # Yield all the k-mer tuples\n",
      "        for r in fastaKmerParserIslands(fastaFh, k, islands): yield r"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from collections import defaultdict\n",
      "def hg19Sample(k, n=100000, chromosome='chr22'):\n",
      "    ''' Given given k, and n, sample n k-mers from both inside\n",
      "        and outside CpG islands, then return histograms of number\n",
      "        of times each k-mer occurs inside and outside. '''\n",
      "    # sample 100K inside and outside k-mers\n",
      "    ins, out = ReservoirSampler(n), ReservoirSampler(n)\n",
      "    for kmer, isle in hg19ChrKmers(k, chromosome=chromosome):\n",
      "        if isle == 'island': ins.add(kmer)\n",
      "        else: out.add(kmer)\n",
      "    return ins, out"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def hist(xs):\n",
      "    ''' Convert iterable to histogram '''\n",
      "    hist = defaultdict(int)\n",
      "    for x in xs: hist[x] += 1\n",
      "    return hist"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ra.seed(723444)\n",
      "q = 'CGCGC'\n",
      "n = 100000\n",
      "ins, out = hg19Sample(len(q), n, 'chr22')\n",
      "inHist, outHist = hist(ins.r), hist(out.r)\n",
      "# print info about inside/outside counts and probabilities\n",
      "print \"inside: %d out of %d\" % (inHist[q], n)\n",
      "print \"outside: %d out of %d\" % (outHist[q], n)\n",
      "print \"p(inside): %0.5f\" % (float(inHist[q]) / (inHist[q] + outHist[q]))\n",
      "print \"p(outside): %0.5f\" % (float(outHist[q]) / (inHist[q] + outHist[q]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "inside: 315 out of 100000\n",
        "outside: 12 out of 100000\n",
        "p(inside): 0.96330\n",
        "p(outside): 0.03670\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Now to build inside and outside Markov chains\n",
      "\n",
      "# compile dinucleotide tables\n",
      "inDinucSamp, outDinucSamp = hg19Sample(2, n=100000, chromosome='chr22')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "def dinucsToTransitionTable(dinucs, stopAfter=None):\n",
      "    ''' Given dinucleotide frequencies, make a transition table. '''\n",
      "    tab = np.zeros((4, 4), dtype=np.float64)\n",
      "    for i in xrange(0, 4):\n",
      "        tot = 0\n",
      "        for j in xrange(0, 4):\n",
      "            tot += dinucs.get('ACGT'[i] + 'ACGT'[j], 0)\n",
      "        if tot > 0:\n",
      "            for j in xrange(0, 4):\n",
      "                tab[i, j] = dinucs.get('ACGT'[i] + 'ACGT'[j], 0) / float(tot)\n",
      "    return tab"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dinucsToTransitionTable(hist(inDinucSamp.r))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "array([[ 0.18357139,  0.27529681,  0.4036642 ,  0.13746761],\n",
        "       [ 0.18946172,  0.3611801 ,  0.25080039,  0.19855779],\n",
        "       [ 0.17342127,  0.33099558,  0.35415423,  0.14142891],\n",
        "       [ 0.09715985,  0.34002913,  0.38002185,  0.18278917]])"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "iTab, nTab = dinucsToTransitionTable(hist(inDinucSamp.r)),\\\n",
      "             dinucsToTransitionTable(hist(outDinucSamp.r))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# transition probabilities inside CpG island\n",
      "iTab"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "array([[ 0.18357139,  0.27529681,  0.4036642 ,  0.13746761],\n",
        "       [ 0.18946172,  0.3611801 ,  0.25080039,  0.19855779],\n",
        "       [ 0.17342127,  0.33099558,  0.35415423,  0.14142891],\n",
        "       [ 0.09715985,  0.34002913,  0.38002185,  0.18278917]])"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# confirm that rows add to 1\n",
      "np.sum(iTab, 1), np.sum(nTab, 1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "(array([ 1.,  1.,  1.,  1.]), array([ 1.,  1.,  1.,  1.]))"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# elementwise log2 of above table\n",
      "np.log2(iTab)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "array([[-2.4455869 , -1.86094019, -1.30877247, -2.8628364 ],\n",
        "       [-2.40002174, -1.4692097 , -1.99538847, -2.33236911],\n",
        "       [-2.52764722, -1.59511613, -1.49755031, -2.82185101],\n",
        "       [-3.36349593, -1.55626975, -1.39584573, -2.45174747]])"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# log ratio table\n",
      "np.log2(iTab) - np.log2(nTab)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "array([[-0.68702616,  0.51274306,  0.49186554, -0.70734831],\n",
        "       [-0.80195953,  0.31979356,  1.99457808, -0.67661388],\n",
        "       [-0.58990756,  0.5264245 ,  0.26127348, -0.59540562],\n",
        "       [-0.87418734,  0.54594298,  0.37698699, -0.69723849]])"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def classify(seq, lrTab):\n",
      "    ''' Classify seq using given log-ratio table '''\n",
      "    bits = 0\n",
      "    nucmap = { 'A':0, 'C':1, 'G':2, 'T':3 }\n",
      "    for dinuc in [ seq[i:i+2] for i in xrange(0, len(seq)-1) ]:\n",
      "        i, j = nucmap[dinuc[0]], nucmap[dinuc[1]]\n",
      "        bits += lrTab[i, j]\n",
      "    return bits"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "classify('CGCGCGCGCGCGCGCGCGCGCGCGCG', np.log2(iTab) - np.log2(nTab))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "32.246609048643101"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "classify('ATTCTACTATCATCTATCTATCTTCT', np.log2(iTab) - np.log2(nTab))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "-9.5012097651362559"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "itest, otest = hg19Sample(100, 1000, 'chr18')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "itestClass = [ classify(x, np.log2(iTab) - np.log2(nTab)) for x in itest.r ]\n",
      "otestClass = [ classify(x, np.log2(iTab) - np.log2(nTab)) for x in otest.r ]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy\n",
      "from matplotlib import pyplot\n",
      "%pylab inline\n",
      "bins = numpy.linspace(-60, 60, 100)\n",
      "pyplot.hist(itestClass, bins, alpha=0.5)\n",
      "pyplot.hist(otestClass, bins, alpha=0.5)\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD9CAYAAAChtfywAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF8BJREFUeJzt3W9sU9f9x/GPaZm6aWwUiThVw2RKCYR/IR2o07QKs9Sg\ndSUNhdGuXWa6dg+m3x7QB4jsyRoeQIyqaqLttEkV3fJbpW5Im9KoSyOg4N+WspIuJe0KHWEj6QKN\nXUjq0tKQBOf8HjDMTYgT2/GfHPv9kpDse33vPSfVPrO+x+cclzHGCABgrRm5bgAAYGoIcgCwHEEO\nAJYjyAHAcgQ5AFiOIAcAy92cyIc8Ho++8pWv6KabbtLMmTPV1tam/v5+PfTQQ/rggw/k8Xi0f/9+\nzZ49O9PtBQCMkdA3cpfLpWAwqOPHj6utrU2SFAgE5PP51NnZqcrKSgUCgYw2FAAwvoRLK2PnDTU1\nNcnv90uS/H6/Ghsb09syAEBCEv5Gfu+992rVqlV64YUXJEnhcFhut1uS5Ha7FQ6HM9dKAEBcCdXI\n33jjDd122206f/68fD6fFi9ePOq8y+WSy+W64brxjgEAJpfM6ikJfSO/7bbbJElz587Vxo0b1dbW\nJrfbrVAoJEnq7e1VUVFR3Mbk67+nnnoq522gb/SP/uXfv2RNGuSff/65Pv30U0nSpUuXdODAAS1f\nvlxVVVVqaGiQJDU0NKi6ujrphwMApm7S0ko4HNbGjRslSVeuXNGjjz6qdevWadWqVdqyZYv27dsX\n+/khACD7Jg3y+fPnq6Oj44bjc+bM0aFDhzLSKFt4vd5cNyFj8rlvEv2zXb73L1kuk0pBJtGbu1wp\n1XsAoJAlm51M0QcAyxHkAGA5ghwALJfQhCCkX21drUKRq7/DL55drEAda9UASA1BniOhSEieao8k\nqbuxO6dtAWA3SisAYDmCHAAsR5ADgOUIcgCwHEEOAJYjyAHAcgQ5AFiOIAcAyxHkAGA5ghwALEeQ\nA4DlCHIAsBxBDgCWI8gBwHIEOQBYjvXILcbmFAAkgtxqbE4BQKK0AgDWI8gBwHIEOQBYjiAHAMsR\n5ABgOYIcACxHkAOA5QhyALAcQQ4AliPIAcByBDkAWI4gBwDLEeQAYLmEgjwajaqiokIbNmyQJPX3\n98vn86m0tFTr1q1TJBLJaCMBAPElFOR79+7VkiVL5HK5JEmBQEA+n0+dnZ2qrKxUIMA62ACQK5MG\n+dmzZ9Xc3KwnnnhCxhhJUlNTk/x+vyTJ7/ersbExs60EAMQ16cYSTz75pJ5++mldvHgxdiwcDsvt\ndkuS3G63wuFw3Ovr6upir71er7xeb+qtLXDOHYEkqb2jPbaxRHt7u7Zu2yrpxt2C2EkImN6CwaCC\nwWDK108Y5K+++qqKiopUUVER9yEulytWchmPM8gxNc4dgSSpta019nogOhB3tyB2EgKmt7Ffcnfu\n3JnU9RMG+dGjR9XU1KTm5mZdvnxZFy9eVE1Njdxut0KhkIqLi9Xb26uioqKUGg8AmLoJa+S7d+9W\nT0+Purq69Pvf/17f/va39bvf/U5VVVVqaGiQJDU0NKi6ujorjQUA3Cip35FfK6HU1tbq4MGDKi0t\n1eHDh1VbW5uRxgEAJjfpYOc1a9as0Zo1ayRJc+bM0aFDhzLWKEyNc+BTij8oKjH4CeSDhIMc9nAO\nfErxB0UlBj+BfMAUfQCwHEEOAJajtDINMYEHQDII8mmICTwAkkFpBQAsR5ADgOUIcgCwHDXyNEtl\noHKiCTwAMBmCPM1SGaicaAIPAEyG0goAWI4gBwDLUVoBsqS2do9CoQFJUnHxFxUI7Mhxi5AvCHIg\nS0KhAXk8dZKk7u66nLYF+YXSCgBYjiAHAMsR5ABgOWrkBc45GYmVFgE7EeQFzjkZiZUWATtRWgEA\nyxHkAGA5SitImnNhMKnwauvOiT0Sk3uQewQ5kuZcGEwqvNq6c2KPxOQe5B6lFQCwHEEOAJajtDLN\nOX/nzYYTAMbDN/Jp7trvvD3VHg0MDUx+AYCCQ5ADgOUIcgCwHEEOAJYjyAHAcgQ5AFiOIAcAyxHk\nAGA5JgRhytic4rp0L6iV6v2c17GoV/6bMMgvX76sNWvWaHBwUENDQ3rggQdUX1+v/v5+PfTQQ/rg\ngw/k8Xi0f/9+zZ49O1ttxjTD5hTXpXtBrVTv57yORb3y34SllVtuuUVHjhxRR0eH3n33XR05ckSt\nra0KBALy+Xzq7OxUZWWlAoHC/QYGALk2aY38S1/6kiRpaGhI0WhUt956q5qamuT3+yVJfr9fjY2N\nmW0lACCuSYN8ZGREK1eulNvt1tq1a7V06VKFw2G53W5JktvtVjgcznhDAQDjm3Swc8aMGero6NAn\nn3yi9evX68iRI6POu1wuuVyuuNfX1dXFXnu9Xnm93pQbC2QTA4bIlmAwqGAwmPL1Cf9q5atf/aq+\n+93vqr29XW63W6FQSMXFxert7VVRUVHc65xBDtiEAUNky9gvuTt37kzq+glLKxcuXFAkEpEkDQwM\n6ODBg6qoqFBVVZUaGhokSQ0NDaqurk6y2QCAdJnwG3lvb6/8fr9GRkY0MjKimpoaVVZWqqKiQlu2\nbNG+fftiPz8EAOTGhEG+fPlyvf322zccnzNnjg4dOpSxRmH6qa2rVSgSksRORcB0wxR9JCQUCbFT\nETBNEeQAYDmCHAAsR5ADgOVY/RDjcg5uSgxwtre3a+vWuv++/oc8npw2BxiFIMe4rg1uXtPa1pq7\nxkwDAwM3xSYHtbYybwLTC6UVALAcQQ4AlqO0ghjnTj+p1sQLcbegdNfPnYt1JXO/VK+D/QhyxDh3\n+km1Jl6IuwWlu37uXKwrmfuleh3sR2kFACxHkAOA5QhyALAcNfIsYYINJsJAJaaCIM8SJthgIgxU\nYioorQCA5QhyALAcpRUUNGdturj4iwoEduS4RUDyCHIUNGdturu7LqdtAVJFaQUALEeQA4DlCHIA\nsBw1ciAHnCsmXn1v1yQgBomnF4IcyAHniomSfZOAGCSeXiitAIDlCHIAsBylFeSEcxGxQtlJKB2c\ntfV01KapdecHghw54VxErFB2EkoHZ209HbVpat35gdIKAFiOIAcAyxHkAGA5auTIufb2dm3dtjX2\nPh2DnwziXTd28lGh/z3yEUGOnBuIDozaPSkdg58M4l03dvJRof898hGlFQCwHEEOAJajtIJpx1kz\nT/dkIWft/Oqz7FqsChjPpN/Ie3p6tHbtWi1dulTLli3Ts88+K0nq7++Xz+dTaWmp1q1bp0gkkvHG\nojBcq5l7qj2x2Z/pcq12fu3fwEA0rfcHcmHSIJ85c6Z+8Ytf6MSJE3rzzTf1y1/+Uu+//74CgYB8\nPp86OztVWVmpQIAp1gCQC5MGeXFxsVauXClJ+vKXv6yysjKdO3dOTU1N8vv9kiS/36/GxsbMthQA\nMK6kBju7u7t1/Phx3X333QqHw3K73ZIkt9utcDickQYCACaW8GDnZ599pk2bNmnv3r2aNWvWqHMu\nl0sul2vc6+rq6mKvvV6vvF5vSg0FMFqquww5r2Owd3oIBoMKBoMpX59QkA8PD2vTpk2qqalRdfXV\nnUzcbrdCoZCKi4vV29uroqKica91BjmA9El1lyHndbbtTJSvxn7J3blzZ1LXT1paMcbo8ccf15Il\nS7Rt27bY8aqqKjU0NEiSGhoaYgEPAMiuSb+Rv/HGG3rppZe0YsUKVVRUSJLq6+tVW1urLVu2aN++\nffJ4PNq/f3/GGwsAuNGkQf6tb31LIyMj4547dOhQ2hsE5IrtO9unWzoWHmPxsuxgZifwX7bvbJ9u\n6Vh4jMXLsoO1VgDAcgQ5AFiO0koCnDu+S+z6jsTly2+2GT+Y3gjyBDh3fJfY9R2Jy5ffbDN+ML1R\nWgEAyxHkAGA5ghwALEeNfIrGDoS2d7SPqqfjOuffKtG/k3O3IGn0QLPzfgxA545zIDQTk36YVDQ5\ngnyKxg6Etra15q4x05zzb5Xo3+nabkHXOAeanfdjADp3nAOhmZj0w6SiyVFaAQDLEeQAYDlKKylw\n1m2piWdXe/t71yfY/Pu92N/+z39u0Z2tK2Of+yT8uf7n8brsNxDIAYI8Bc66LTXx7BoYuHJ9gs17\n1/eJHTRDKrn/+iSVj158MdtNA3KG0goAWI4gBwDLEeQAYDlq5LBKX9+HagxulSSd7T2jxsagJGlw\ncCjuNfmyAiEQD0EOq0RvuqLZXs/V16ejmj3bK0kaGXkr7jX5sgIhEA+lFQCwHEEOAJajtIKc6OuL\nxOrbfX2RhD4nTVwLx9RkehegsfdnAaz0IciRE9HoSKy+fSb6bkKfkyauhWNqMr0L0Nj7swBW+lBa\nAQDLEeQAYDmCHAAsV9A18rG7+8TbfSbVFQ5ZJRGFJpkB00zvLFRICjrIx+7uE2/3mVRXOGSVRBSa\nZAZMM72zUCGhtAIAliPIAcByBV1aQWY5xwgkxgkQH/XyqSHIkTHOMQKJcQLER718aiitAIDlCHIA\nsBxBDgCWo0YOwEq1tXsUCg3E3v/zn+9o8eJySYU3YDrpN/If/ehHcrvdWr58eexYf3+/fD6fSktL\ntW7dOkUi8ZchBYBMCIUG5PHUxf5duGBir50BXwgmDfLHHntMLS0to44FAgH5fD51dnaqsrJSgUAg\nYw0EAExs0iC/5557dOutt4461tTUJL/fL0ny+/1qbGzMTOsAAJNKqUYeDofldrslSW63W+FwOO5n\n6+rqYq+9Xq+8Xm8qj8wKFrnKjcHBwVG7AJ09G4q9Z0cgODknDqV7ByNpdN09m3X2YDCoYDCY8vVT\nHux0uVxyuVxxzzuDfLpjkavcGBnRqF2AotG3Yu/ZEQhOzolD6d7BSLped5eyOzFp7JfcnTt3JnV9\nSj8/dLvdCoWuLvHa29uroqKiVG4DAEiDlIK8qqpKDQ0NkqSGhgZVV6f//xkBAImZNMi///3v65vf\n/KZOnTqlefPm6Te/+Y1qa2t18OBBlZaW6vDhw6qtrc1GWwEA45i0Rv7yyy+Pe/zQoUNpbwzyW19f\nhEFMIAOY2YmsiUZHGMQEMoC1VgDAcgQ5AFiO0grSylkHd07skaiLI3ucE4ek/F9EiyBHWjnr4M6J\nPRJ1cWSPc+KQlP+7DlFaAQDLEeQAYDmCHAAsR40ceeny5Ygag1slSWcj/xd73XfpRO4ahYSMHahM\nxyqHzns6Bz7H7jKUiRUVs4EgR14yM6Oa7fVIkqKnh2Kvz5w5nLtGISFjByrTscqh857OgU/naofp\nelYuUFoBAMsR5ABgOUorQBx9l05QW89Dmd5lKBcIciCO6MxBaut5KNO7DOUCpRUAsBxBDgCWo7SC\nKWPDCCC3CHJMGRtGALlFaQUALEeQA4DlCHIAsFxB1Mhr62oVioQkScWzixWoC+S4RdOTc9Cyry8S\n95zNO/84F9OSRk/0cU4AkqTBK6P/BvE+x2Sh/DTRLkPOxbbG7j400blMKYggD0VC8lR7JEndjd05\nbct05hy0PBN9N+45m3f+cS6mJY2e6OOcACRJI6ej495j7OeYLJSfJtplyLnY1tjdhyY6lymUVgDA\ncgQ5AFgub0oridbB29vbtXXb1quvO9pjJZdC4ax1X3uPyTlr62Nr5/E2sRj7nlo6MiVvgjzROvhA\ndCD2uda21sw3bJpx1rqlG2vhGJ+ztj62dh5vE4ux76mlI1MorQCA5QhyALAcQQ4Alpt2NXLnoKU0\neuCSiT2pSWV1wsHBQWsn/SRjokFMwBbTLsidg5bS6IFLJvakJpXVCUdGZO2kn2RMNIgJ2ILSCgBY\nLuPfyM+fPy9JmjVrlm655ZZMPw4ACk7Gg7x2b62GB4d139336eHNDyd9fbwJPM7jY88VonxZ1KpQ\nTLTwlvPc2AlGE00qcl430WJgTExKjXMRrfb2f8jjSe6aTC6uNaUgb2lp0bZt2xSNRvXEE09ox44b\nGzLv3nkK/zusy0OXU3pGvAk8zuNjz2VLd0d31p8ZT7oXtTKD+V0vjnR3a3Yi/0vMkIkW3nKeGzvB\naKJJRc7rOk82J/QsW3V3B7P+TOciWq2t1Ulfk8nFtVKukUejUf30pz9VS0uLTp48qZdfflnvv//+\nlBpjm+kU5OlmhkZy3YSMinR357oJGRUdHMx1EzIqF0E+naUc5G1tbbrzzjvl8Xg0c+ZMPfzww3rl\nlVfS2TYAQAJSLq2cO3dO8+bNi70vKSnRsWPHbvhc8H+DGr48rLllc1N9FABgAi5jjEnlwj/+8Y9q\naWnRCy+8IEl66aWXdOzYMT333HPXb+5ypaeVAFBgkonmlL+R33777erp6Ym97+npUUlJScoNAQCk\nJuUa+apVq3T69Gl1d3draGhIf/jDH1RVVZXOtgEAEpDyN/Kbb75Zzz//vNavX69oNKrHH39cZWVl\n6WwbACABU5qi/53vfEenTp3Sv/71L/3sZz+LHX/uuedUVlamZcuWjfpteX19vRYuXKjFixfrwIED\nU3l0zj3zzDOaMWOG+vv7Y8fyoX/bt29XWVmZysvL9eCDD+qTTz6JncuH/klX5z8sXrxYCxcu1J49\ne3LdnCnp6enR2rVrtXTpUi1btkzPPvusJKm/v18+n0+lpaVat26dIhG7FwSLRqOqqKjQhg0bJOVX\n/yKRiDZv3qyysjItWbJEx44dS75/Js0OHz5s7r33XjM0NGSMMeajjz4yxhhz4sQJU15eboaGhkxX\nV5dZsGCBiUaj6X58VvznP/8x69evNx6Px/T19Rlj8qd/Bw4ciLV7x44dZseOHcaY/OnflStXzIIF\nC0xXV5cZGhoy5eXl5uTJk7luVsp6e3vN8ePHjTHGfPrpp6a0tNScPHnSbN++3ezZs8cYY0wgEIj9\nd7TVM888Yx555BGzYcMGY4zJq/798Ic/NPv27TPGGDM8PGwikUjS/Ut7kH/ve98zr7/++g3Hd+/e\nbQKBQOz9+vXrzd/+9rd0Pz4rNm/ebN55551RQZ5P/bvmT3/6k3n00UeNMfnTv6NHj5r169fH3tfX\n15v6+voctii9HnjgAXPw4EGzaNEiEwqFjDFXw37RokU5blnqenp6TGVlpTl8+LC5//77jTEmb/oX\niUTM/PnzbziebP/Svvrh6dOn9Ze//EXf+MY35PV69fe//12S9OGHH476VUtJSYnOnTuX7sdn3Cuv\nvKKSkhKtWLFi1PF86Z/Tiy++qPvuu09S/vRvvPkPNvZjPN3d3Tp+/LjuvvtuhcNhud1uSZLb7VY4\nHM5x61L35JNP6umnn9aMGdfjKl/619XVpblz5+qxxx7TXXfdpR//+Me6dOlS0v1LabDT5/MpFArd\ncHzXrl26cuWKPv74Y7355pt66623tGXLFp05c2bc+0zX35lP1L/6+vpR9WEzwU8sbevf7t27YzXI\nXbt26Qtf+IIeeeSRuPeZrv2biI1tTsRnn32mTZs2ae/evZo1a9aocy6Xy9p+v/rqqyoqKlJFRYWC\nweC4n7G5f1euXNHbb7+t559/XqtXr9a2bdsUCIzeMCeR/qUU5AcPHox77le/+pUefPBBSdLq1as1\nY8YMXbhw4YbfnZ89e1a33357Ko/PuHj9e++999TV1aXy8nJJV/vw9a9/XceOHcuL/l3z29/+Vs3N\nzXr99ddjx2zq30QSmf9gm+HhYW3atEk1NTWqrr66mJPb7VYoFFJxcbF6e3tVVFSU41am5ujRo2pq\nalJzc7MuX76sixcvqqamJm/6V1JSopKSEq1evVqStHnzZtXX16u4uDi5/qW75vPrX//a/PznPzfG\nGHPq1Ckzb948Y8z1wbLBwUFz5swZc8cdd5iRkZF0Pz6rxhvstL1/r732mlmyZIk5f/78qOP50r/h\n4WFzxx13mK6uLjM4OGj9YOfIyIipqakx27ZtG3V8+/btsTGN+vp6qwcDrwkGg7EaeT7175577jGn\nTp0yxhjz1FNPme3btyfdv7QH+dDQkPnBD35gli1bZu666y5z5MiR2Lldu3aZBQsWmEWLFpmWlpZ0\nPzrr5s+fHwtyY/Kjf3feeaf52te+ZlauXGlWrlxpfvKTn8TO5UP/jDGmubnZlJaWmgULFpjdu3fn\nujlT8te//tW4XC5TXl4e+2/22muvmb6+PlNZWWkWLlxofD6f+fjjj3Pd1CkLBoOxX63kU/86OjrM\nqlWrzIoVK8zGjRtNJBJJun8pr7UCAJge2LMTACxHkAOA5QhyALAcQQ4AliPIAcByBDkAWO7/AYfn\nY0h9QHSPAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x102da3bd0>"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 25
    }
   ],
   "metadata": {}
  }
 ]
}