Skip to content

Instantly share code, notes, and snippets.

@BenLangmead
Last active December 28, 2015 00:29
Show Gist options
  • Save BenLangmead/7413873 to your computer and use it in GitHub Desktop.
Save BenLangmead/7413873 to your computer and use it in GitHub Desktop.
{
"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": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment