Meosis.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:4c6044d0237c3473eb12db91038c824c57bddea969f0999a0f51e5505f37a04e"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import vcf\n",
      "%pylab inline\n",
      "\n",
      "print \"pyccf version:\", vcf.VERSION"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n",
        "pyccf version: 0.6.7\n"
       ]
      }
     ],
     "prompt_number": 383
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Do some Tt variants lose at meosis?\n",
      "\n",
      "`gatk.raw.vcf` conatins a while genome's work of (naive) variant calls, made assuming both the ancestors and the descendants are diploid. I want to find likely heterozosities in the ancestral strain (M0) and see how each allele breaks out in the descendants. In particular, since the odd m47/m51 patterns might be explaned by some hunk of MIC chormosome note mkaing into to MAC nuclei, I want to see if those two lines share the same het. allele more that we might expect. \n",
      "\n",
      "The \"best\" heteozygote sites will be those where the ancestro was called heterozygous, and no descendants where. So let's start with isolating just these sites"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "snps = vcf.Reader(filename=\"gatk.raw.vcf\")\n",
      "hets = [r for r in snps if [s.sample for s in r.get_hets()] == [\"TtM0\"]]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 384
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = len(hets)\n",
      "print \"{0} hets, meaning theta is ~ {1}\".format(n, n/1e8)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "4901 hets, meaning theta is ~ 4.901e-05\n"
       ]
      }
     ],
     "prompt_number": 385
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def descendant_maf(rec):\n",
      "    n= sum([s.called for s in r.samples[1:]])\n",
      "    if n == 0:\n",
      "        return 0\n",
      "    aaf =  len(r.get_hom_alts())/float(n)\n",
      "    if aaf > 0.5:\n",
      "        aaf = aaf-0.5\n",
      "    return aaf"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 386
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figure()\n",
      "hist([descendant_maf(r) for r in hets])\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF7lJREFUeJzt3W1sU+fdx/Hf6W1rHYUg2hGH2WguiWkwDcGBmlW7kYIg\nKaQiSkUVNdsgaek0JS8Goi/WddqAaloyaVsFbJmqNZsyOpWgdkvQvcXK2OJqoyVUWRBTvQ3D0il2\nHjZGs4Y1LQ8594sUrymQmITYptf3I1k6vnwe/ueyfX7nwSexbNu2BQAwzh3pLgAAkB4EAAAYigAA\nAEMRAABgKAIAAAxFAACAoSYNgPfee09r1qzRypUr5ff79bWvfU2SdP78eZWUlGjp0qUqLS3V8PBw\nYpr6+nr5fD7l5+ero6Mj0d7d3a2CggL5fD7t2LFjllYHAJCsSQPgzjvvVGdnp06ePKlTp06ps7NT\nf/jDH9TQ0KCSkhKdPn1a69evV0NDgyQpEomopaVFkUhEoVBIdXV1unqbQW1trZqamhSNRhWNRhUK\nhWZ/7QAANzTlKaA5c+ZIki5evKgrV65owYIFOnLkiKqrqyVJ1dXVam1tlSS1tbWpqqpKTqdTXq9X\neXl56urq0sDAgEZGRhQMBiVJ27ZtS0wDAEiPKQNgbGxMK1eulMvl0rp167R8+XINDQ3J5XJJklwu\nl4aGhiRJ/f398ng8iWk9Ho/i8fg17W63W/F4/FavCwDgJjimGuGOO+7QyZMn9e9//1sPPfSQOjs7\nJ7xuWZYsy5q1AgEAs2PKALhq/vz5evjhh9Xd3S2Xy6XBwUHl5ORoYGBA2dnZksb37Pv6+hLTxGIx\neTweud1uxWKxCe1ut/uaZeTl5ens2bMzWR8AME5ubq7OnDlz09NNegro3LlziV/4jI6O6je/+Y0C\ngYDKy8vV3NwsSWpublZFRYUkqby8XIcOHdLFixfV29uraDSqYDConJwcZWVlqaurS7Zt6+DBg4lp\nPuzs2bOybZuHbWv37t1pryFTHvQFfUFfTP6Y7o7zpEcAAwMDqq6u1tjYmMbGxrR161atX79egUBA\nlZWVampqktfr1eHDhyVJfr9flZWV8vv9cjgcamxsTJweamxsVE1NjUZHR1VWVqaNGzdOq2AAwK0x\naQAUFBToj3/84zXtd999t44ePXrdaZ555hk988wz17SvWrVKf/rTn6ZZJgDgVuNO4AxVXFyc7hIy\nBn3xX/TFf9EXM2fZtp0x/xDGsix94hNz01qDw+HQqVPdWrJkSVrrAIBkWZal6WzKk/4VUKq8/35/\nWpf/iU/874Q/bQEAH1cZFwDSvDQv/3/SvHwASA2uAQCAoQgAADAUAQAAhiIAAMBQBAAAGIoAAABD\nEQAAYCgCAAAMRQAAgKEIAAAwFAEAAIYiAADAUAQAABiKAAAAQxEAAGAoAgAADEUAAIChCAAAMBQB\nAACGIgAAwFAEAAAYigAAAEMRAABgKAIAAAw1aQD09fVp3bp1Wr58ue6//37t379fkrRnzx55PB4F\nAgEFAgG1t7cnpqmvr5fP51N+fr46OjoS7d3d3SooKJDP59OOHTtmaXUAAMlyTPai0+nUc889p5Ur\nV+rChQtatWqVSkpKZFmWdu3apV27dk0YPxKJqKWlRZFIRPF4XBs2bFA0GpVlWaqtrVVTU5OCwaDK\nysoUCoW0cePGWV05AMCNTXoEkJOTo5UrV0qS5s6dq2XLlikej0uSbNu+Zvy2tjZVVVXJ6XTK6/Uq\nLy9PXV1dGhgY0MjIiILBoCRp27Ztam1tvdXrAgC4CUlfA3jrrbfU09Ojz372s5KkAwcOqLCwUNu3\nb9fw8LAkqb+/Xx6PJzGNx+NRPB6/pt3tdieCBACQHkkFwIULF/Too49q3759mjt3rmpra9Xb26uT\nJ09q0aJFeuqpp2a7TgDALTbpNQBJunTpkrZs2aIvfvGLqqiokCRlZ2cnXn/yySe1efNmSeN79n19\nfYnXYrGYPB6P3G63YrHYhHa3232DJe750HDxBw8AwFXhcFjhcHjmM7InMTY2Zm/dutXeuXPnhPb+\n/v7E8Pe//327qqrKtm3bfvPNN+3CwkL7/ffft//2t7/ZS5YsscfGxmzbtu1gMGgfP37cHhsbszdt\n2mS3t7dfszxJtmSn9ZGVFbC7u7sn6xYAyChTbMpvaNIjgGPHjunFF1/UihUrFAgEJEnf/va39dJL\nL+nkyZOyLEv33nuvnn/+eUmS3+9XZWWl/H6/HA6HGhsbZVmWJKmxsVE1NTUaHR1VWVkZvwACgDSz\nPkiPjDAeFuktJyurSJ2dL6ioqCitdQBAsizLuu4vM6fCncAAYCgCAAAMRQAAgKEIAAAwFAEAAIYi\nAADAUAQAABiKAAAAQxEAAGAoAgAADEUAAIChCAAAMBQBAACGIgAAwFAEAAAYigAAAEMRAABgKAIA\nAAxFAACAoQgAADAUAQAAhiIAAMBQBAAAGIoAAABDEQAAYCgCAAAMRQAAgKEIAAAwFAEAAIaaNAD6\n+vq0bt06LV++XPfff7/2798vSTp//rxKSkq0dOlSlZaWanh4ODFNfX29fD6f8vPz1dHRkWjv7u5W\nQUGBfD6fduzYMUurAwBI1qQB4HQ69dxzz+nNN9/U8ePH9cMf/lB//vOf1dDQoJKSEp0+fVrr169X\nQ0ODJCkSiailpUWRSEShUEh1dXWybVuSVFtbq6amJkWjUUWjUYVCodlfOwDADU0aADk5OVq5cqUk\nae7cuVq2bJni8biOHDmi6upqSVJ1dbVaW1slSW1tbaqqqpLT6ZTX61VeXp66uro0MDCgkZERBYNB\nSdK2bdsS0wAA0iPpawBvvfWWenp6tGbNGg0NDcnlckmSXC6XhoaGJEn9/f3yeDyJaTwej+Lx+DXt\nbrdb8Xj8Vq0DAGAaHMmMdOHCBW3ZskX79u3TvHnzJrxmWZYsy7qFJe350HDxBw8AwFXhcFjhcHjG\n85kyAC5duqQtW7Zo69atqqiokDS+1z84OKicnBwNDAwoOztb0viefV9fX2LaWCwmj8cjt9utWCw2\nod3tdt9giXumvzYAYIDi4mIVFxcnnu/du3da85n0FJBt29q+fbv8fr927tyZaC8vL1dzc7Mkqbm5\nOREM5eXlOnTokC5evKje3l5Fo1EFg0Hl5OQoKytLXV1dsm1bBw8eTEwDAEiPSY8Ajh07phdffFEr\nVqxQIBCQNP4zz6efflqVlZVqamqS1+vV4cOHJUl+v1+VlZXy+/1yOBxqbGxMnB5qbGxUTU2NRkdH\nVVZWpo0bN87yqgEAJmPZV3+nmQHGwyK95WRlFamz8wUVFRWltQ4ASJZlWZrOppw7gQHAUAQAABiK\nAAAAQxEAAGAoAgAADEUAAIChCAAAMBQBAACGIgAAwFAEAAAYigAAAEMRAABgKAIAAAxFAACAoQgA\nADAUAQAAhiIAAMBQBAAAGIoAAABDEQAAYCgCAAAMRQAAgKEIAAAwFAEAAIYiAADAUAQAABiKAAAA\nQxEAAGCoKQPgiSeekMvlUkFBQaJtz5498ng8CgQCCgQCam9vT7xWX18vn8+n/Px8dXR0JNq7u7tV\nUFAgn8+nHTt23OLVAADcrCkD4PHHH1coFJrQZlmWdu3apZ6eHvX09GjTpk2SpEgkopaWFkUiEYVC\nIdXV1cm2bUlSbW2tmpqaFI1GFY1Gr5knACC1pgyAtWvXasGCBde0X92wf1hbW5uqqqrkdDrl9XqV\nl5enrq4uDQwMaGRkRMFgUJK0bds2tba23oLyAQDTNe1rAAcOHFBhYaG2b9+u4eFhSVJ/f788Hk9i\nHI/Ho3g8fk272+1WPB6fQdkAgJlyTGei2tpaffOb35QkfeMb39BTTz2lpqamW1TSng8NF3/wAABc\nFQ6HFQ6HZzyfaQVAdnZ2YvjJJ5/U5s2bJY3v2ff19SVei8Vi8ng8crvdisViE9rdbvcN5r5nOiUB\ngDGKi4tVXFyceL53795pzWdap4AGBgYSw7/85S8TvxAqLy/XoUOHdPHiRfX29ioajSoYDConJ0dZ\nWVnq6uqSbds6ePCgKioqplUwAODWmPIIoKqqSq+++qrOnTunxYsXa+/evQqHwzp58qQsy9K9996r\n559/XpLk9/tVWVkpv98vh8OhxsZGWZYlSWpsbFRNTY1GR0dVVlamjRs3zu6aAR9DWVl3a2Tk7XSX\noXnzFuidd86nuwzMkGVf7+c8aTIeFuktJyurSJ2dL6ioqCitdQDXkwnfkXHWdX8JiPSwrOm9H9wJ\nDACGIgAAwFAEAAAYigAAAEMRAABgKAIAAAxFAACAoQgAADAUAQAAhiIAAMBQBAAAGIoAAABDEQAA\nYCgCAAAMRQAAgKEIAAAwFAEAAIYiAADAUAQAABiKAAAAQxEAAGAoAgAADEUAAIChCAAAMBQBAACG\nIgAAwFAEAAAYigAAAENNGQBPPPGEXC6XCgoKEm3nz59XSUmJli5dqtLSUg0PDydeq6+vl8/nU35+\nvjo6OhLt3d3dKigokM/n044dO27xagAAbtaUAfD4448rFApNaGtoaFBJSYlOnz6t9evXq6GhQZIU\niUTU0tKiSCSiUCikuro62bYtSaqtrVVTU5Oi0aii0eg18wQApNaUAbB27VotWLBgQtuRI0dUXV0t\nSaqurlZra6skqa2tTVVVVXI6nfJ6vcrLy1NXV5cGBgY0MjKiYDAoSdq2bVtiGgBAekzrGsDQ0JBc\nLpckyeVyaWhoSJLU398vj8eTGM/j8Sgej1/T7na7FY/HZ1I3AGCGHDOdgWVZsizrVtTygT0fGi7+\n4AEAuCocDiscDs94PtMKAJfLpcHBQeXk5GhgYEDZ2dmSxvfs+/r6EuPFYjF5PB653W7FYrEJ7W63\n+wZz3zOdkgDAGMXFxSouLk4837t377TmM61TQOXl5WpubpYkNTc3q6KiItF+6NAhXbx4Ub29vYpG\nowoGg8rJyVFWVpa6urpk27YOHjyYmAYAkB5THgFUVVXp1Vdf1blz57R48WI9++yzevrpp1VZWamm\npiZ5vV4dPnxYkuT3+1VZWSm/3y+Hw6HGxsbE6aHGxkbV1NRodHRUZWVl2rhx4+yuGQBgUpZ99Xea\nGWA8LNJbTlZWkTo7X1BRUVFa6wCuJxO+I+MsZdCmw3iWNb33gzuBAcBQBAAAGIoAAABDEQAAYCgC\nAAAMRQAAgKEIAAAwFAEAAIYiAADAUAQAABiKAAAAQxEAAGAoAgAADEUAAIChCAAAMBQBAACGIgAA\nwFAEAAAYigAAAEMRAABgKAIAAAxFAACAoQgAADAUAQAAhiIAAMBQBAAAGIoAAABDEQAAYKgZBYDX\n69WKFSsUCAQUDAYlSefPn1dJSYmWLl2q0tJSDQ8PJ8avr6+Xz+dTfn6+Ojo6ZlY5AGBGZhQAlmUp\nHA6rp6dHJ06ckCQ1NDSopKREp0+f1vr169XQ0CBJikQiamlpUSQSUSgUUl1dncbGxma+BgCAaZnx\nKSDbtic8P3LkiKqrqyVJ1dXVam1tlSS1tbWpqqpKTqdTXq9XeXl5idAAAKTejI8ANmzYoNWrV+vH\nP/6xJGloaEgul0uS5HK5NDQ0JEnq7++Xx+NJTOvxeBSPx2eyeADADDhmMvGxY8e0aNEi/fOf/1RJ\nSYny8/MnvG5ZlizLuuH0k70GAJhdMwqARYsWSZIWLlyoRx55RCdOnJDL5dLg4KBycnI0MDCg7Oxs\nSZLb7VZfX19i2lgsJrfbfZ257vnQcPEHDwDAVeFwWOFweMbzseyPnsRP0rvvvqsrV65o3rx5+s9/\n/qPS0lLt3r1bR48e1T333KOvfvWramho0PDwsBoaGhSJRPT5z39eJ06cUDwe14YNG3TmzJkJRwHj\nw9Mq55bJyipSZ+cLKioqSmsdwPVkwndknHXN9T+kj2VN7/2Y9hHA0NCQHnnkEUnS5cuX9YUvfEGl\npaVavXq1Kisr1dTUJK/Xq8OHD0uS/H6/Kisr5ff75XA41NjYyCkgAEijaR8BzIZM2LvhCACZLBO+\nI+M4Asgk0z0C4E5gADAUAQAAhiIAAMBQBAAAGIoAAABDEQAAYCgCAAAMRQAAgKEIAAAwFAEAAIYi\nAADAUAQAABiKAAAAQxEAAGAoAgAADEUAAIChCAAAMBQBAACGIgAAwFAEAAAYigAAAEM50l0AANzO\nsrLu1sjI2+kuY1oIAACYgfGNv53mKqxpTcUpIAAwFAEAAIbiFBCQpNv5XC9wPQQAJpUJG7158xbo\nnXfOp7UG6fY+1wtcDwGASWXCRm9khI0eMBtSeg0gFAopPz9fPp9P3/nOd1K5aADAR1i2badk9+7K\nlSu67777dPToUbndbj3wwAN66aWXtGzZsv8WY1lK995mVlaRLl/+m959999preOTn5yrd98dSWsN\nUma8J9L/SBpLcw1XpbsvMuH9kCRLKdp03FA4HFZxcXFaa5Ay5TsyvfcjZUcAJ06cUF5enrxer5xO\npx577DG1tbWlavE3ZXzjb6f1MTp6YfZX9LYxpnS/H+n/guOjwuFwuku47aXsGkA8HtfixYsTzz0e\nj7q6ulK1+NuQ9cGeBZCJHBnw+bxDe/fuTXMNt7eUBUCyH5asrM2zXMnk3nvvbFqX/1+ZsteZ7i85\nMtNlpf/zmQmnXqTb+TuSsgBwu93q6+tLPO/r65PH45kwTm5urs6e/b9UlTSFTHhTM6EGKTPqyIQa\npMyoIxNqkDKjjkyoQUp3Hbm5udOaLmUXgS9fvqz77rtPv/3tb/XpT39awWDwmovAAIDUSdkRgMPh\n0A9+8AM99NBDunLlirZv387GHwDSKGVHAACAzJKWPwaXzA1hX/nKV+Tz+VRYWKienp4UV5g6U/XF\nX/7yFz344IO688479b3vfS8NFabOVH3x85//XIWFhVqxYoU+97nP6dSpU2moMjWm6ou2tjYVFhYq\nEAho1apV+t3vfpeGKlMj2RtI33jjDTkcDv3iF79IYXWpNVVfhMNhzZ8/X4FAQIFAQN/61rcmn6Gd\nYpcvX7Zzc3Pt3t5e++LFi3ZhYaEdiUQmjPOrX/3K3rRpk23btn38+HF7zZo1qS4zJZLpi3/84x/2\nG2+8YX/961+3v/vd76ap0tmXTF+89tpr9vDwsG3btt3e3m705+LChQuJ4VOnTtm5ubmpLjMlkumL\nq+OtW7fOfvjhh+2XX345DZXOvmT6orOz0968eXPS80z5EUAyN4QdOXJE1dXVkqQ1a9ZoeHhYQ0ND\nqS511iXTFwsXLtTq1avldDrTVGVqJNMXDz74oObPny9p/HMRi8XSUeqsS6Yv7rrrrsTwhQsX9KlP\nfSrVZaZEsjeQHjhwQI8++qgWLlyYhipTI9m+sG/irH7KA+B6N4TF4/Epx/k4ftmT6QtT3GxfNDU1\nqaysLBWlpVyyfdHa2qply5Zp06ZN2r9/fypLTJlktxdtbW2qra2VlPw9R7ebZPrCsiy99tprKiws\nVFlZmSKRyKTzTPlfA032zfloin0c39SP4zpN1830RWdnp37yk5/o2LFjs1hR+iTbFxUVFaqoqNDv\nf/97bd26VX/9619nubLUS6Yvdu7cqYaGBlnW+N/DuZk94NtJMn1RVFSkvr4+zZkzR+3t7aqoqNDp\n06dvOH7KAyCZG8I+Ok4sFpPb7U5ZjamSTF+YItm+OHXqlL70pS8pFAppwYIFqSwxZW72c7F27Vpd\nvnxZ//rXv3TPPfekosSUSaYvuru79dhjj0mSzp07p/b2djmdTpWXl6e01tmWTF/MmzcvMbxp0ybV\n1dXp/Pnzuvvuu68/01t5kSIZly5dspcsWWL39vba77///pQXgV9//fWP7cW+ZPriqt27d3+sLwIn\n0xd///vf7dzcXPv1119PU5WpkUxfnDlzxh4bG7Nt27a7u7vtJUuWpKPUWXcz3xHbtu2amhr7lVde\nSWGFqZNMXwwODiY+F11dXfZnPvOZSeeZ8iOAG90Q9vzzz0uSvvzlL6usrEy//vWvlZeXp7vuuks/\n/elPU11mSiTTF4ODg3rggQf0zjvv6I477tC+ffsUiUQ0d+7cNFd/ayXTF88++6zefvvtxLlep9Op\nEydOpLPsWZFMX7zyyiv62c9+JqfTqblz5+rQoUNprnp2JNMXpkimL15++WX96Ec/ksPh0Jw5c6b8\nXHAjGAAYKi03ggEA0o8AAABDEQAAYCgCAAAMRQAAgKEIAAAwFAEAAIYiAADAUP8PaIsGNSC1naAA\nAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f6b897d2650>"
       ]
      }
     ],
     "prompt_number": 387
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##It seems like many of the alt alleles are present in only one or two desendants...\n",
      "\n",
      "Are any samples over-represented? Starting with \"singletons\"\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from collections import Counter \n",
      "\n",
      "def extract_singles(rec):\n",
      "    if rec.num_hom_alt == 1:\n",
      "        return rec.get_hom_alts()[0].sample\n",
      "    if rec.num_hom_ref == 1:\n",
      "        return rec.get_hom_refs()[0].sample\n",
      "\n",
      "singletons = Counter([extract_singles(r) for r in hets if extract_singles(r)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 388
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figure()\n",
      "bar(np.arange(11), singletons.values())\n",
      "xticks(np.arange(11), singletons.keys(), rotation=60)\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAEbCAYAAAD+uL7AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtUlHX+B/D3AKISIIIyKJpoiQZeULzlzwseRSzTNDqk\ntspqralbWu2mVpJsZWLhekk9telpyW7iWkRWJkqjZRnkDUtRRFEkQGFAQRCB+fz+YJl0xSB4Zpjv\n+H6dwzkwDM/7+31meM/D95kZdCIiICIim+fQ3AMgIqKGYWETESmChU1EpAgWNhGRIljYRESKYGET\nESnCqSFX8vPzg7u7OxwdHdGiRQukpKTAaDTikUcewdmzZ+Hn54f4+Hh4eHhYerxERLetBh1h63Q6\nGAwGHDp0CCkpKQCAmJgYhIaG4uTJkxg9ejRiYmIsOlAiottdg5dE/vf1NYmJiYiMjAQAREZGIiEh\nQduRERHRDRp8hD1mzBgMGDAA77zzDgAgPz8fer0eAKDX65Gfn2+5URIRUcPWsPft24cOHTrg4sWL\nCA0NRc+ePW/4vk6ng06ns8gAiYioRoMKu0OHDgCA9u3bY/LkyUhJSYFer0deXh58fHyQm5sLb2/v\nm37u7rvvRmZmprYjJiKyc3379sXhw4dvurzeJZGysjKUlJQAAK5cuYKdO3eid+/emDhxIuLi4gAA\ncXFxmDRp0k0/m5mZCRGx2sfSpUvtNs+e58Y85tl6nrU/jhw5Umcf13uEnZ+fj8mTJwMAqqqq8Oij\nj2Ls2LEYMGAAIiIisGnTJvPT+oiIyHLqLeyuXbvWeWju6emJXbt2WWRQRER0M7t6pePrr680nwDV\n8sPd3bPOvJCQEKvNzZpZzGMe82yTTkQs9g8MdDodLLj5OvMAS+RZdx5EdHu7VXfa1RE2EZE9Y2ET\nESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGRIljY\nRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgW\nNhGRIljYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmi\nQYVdXV2Nfv36YcKECQAAo9GI0NBQ+Pv7Y+zYsSguLrboIImIqIGFvWbNGgQEBECn0wEAYmJiEBoa\nipMnT2L06NGIiYmx6CCJiKgBhX3+/Hl8+eWXePzxxyEiAIDExERERkYCACIjI5GQkGDZURIRUf2F\n/cwzz+CNN96Ag8NvV83Pz4derwcA6PV65OfnW26EREQEAHD6vW9u374d3t7e6NevHwwGQ53X0el0\n5qWSukRHR5s/DwkJQUhISGPGSURktwwGwy079no6qV3nqMMLL7yAzZs3w8nJCVevXsXly5fx0EMP\nITU1FQaDAT4+PsjNzcWoUaOQnp5+88Z1OvzO5jVX88BhiTzrzoOIbm+36s7fLezr7dmzB7Gxsfj8\n88+xcOFCeHl5YdGiRYiJiUFxcXGdJx5Z2EREf9ytuvMPPQ+7dulj8eLFSEpKgr+/P5KTk7F48WJt\nRklERLfU4CPsRm2cR9hERH+YJkfYRETUfFjYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmC\nhU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGRIljYRESK\nYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGR\nIljYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESnidwv76tWrGDx4MIKCghAQEIDnn38eAGA0\nGhEaGgp/f3+MHTsWxcXFVhksEdHtTCci8ntXKCsrg4uLC6qqqjBs2DDExsYiMTER7dq1w8KFC7Fi\nxQoUFRUhJibm5o3rdKhn85rS6XQALJFn3XnUxd3dEyUlRRbZtptbW1y+bLTItonoj7tVd9a7JOLi\n4gIAuHbtGqqrq9G2bVskJiYiMjISABAZGYmEhASNh0v/q6asxSIflnogICJt1VvYJpMJQUFB0Ov1\nGDVqFAIDA5Gfnw+9Xg8A0Ov1yM/Pt/hAiYhud071XcHBwQGHDx/GpUuXEBYWhm+++eaG7+t0uv8u\nRRARkSXVW9i12rRpg/Hjx+PAgQPQ6/XIy8uDj48PcnNz4e3tfcufi46ONn8eEhKCkJCQpoyXiMju\nGAwGGAyGeq/3uycdCwoK4OTkBA8PD5SXlyMsLAxLly7F119/DS8vLyxatAgxMTEoLi7mSUcLs9zc\nAFuYHxH95lbd+buFffToUURGRsJkMsFkMmH69Ol47rnnYDQaERERgXPnzsHPzw/x8fHw8PBocKil\nsLAbvfVmnx8R/aZRhW2pUEvmsbAbtfVmnx8R/abRT+sjIiLbwMImIlIEC5uISBEsbCIiRbCwiYgU\nwcImIlIEC5uISBEsbCIiRbCwiYgUwcImIlIEC5uISBEsbCIiRbCwiYgUwcImIlIEC5uISBEsbCIi\nRbCwiYgUwcImIlIEC5uISBEsbCIiRbCwiYgUwcImIlIEC5uISBEsbCIiRbCwiYgUwcImIlIEC5uI\nSBEsbCIiRbCwiYgUwcImIlIEC5uISBEsbCIiRbCwiYgUwcImIlIEC5uISBEsbCIiRdRb2NnZ2Rg1\nahQCAwPRq1cvrF27FgBgNBoRGhoKf39/jB07FsXFxRYfLBHR7UwnIvJ7V8jLy0NeXh6CgoJQWlqK\n4OBgJCQk4N1330W7du2wcOFCrFixAkVFRYiJiblx4zod6tm8pnQ6HQBL5Fl3HnWOwGJzA2xhfkT0\nm1t1Z71H2D4+PggKCgIAuLq64p577kFOTg4SExMRGRkJAIiMjERCQoLGQyYiouv9oTXsrKwsHDp0\nCIMHD0Z+fj70ej0AQK/XIz8/3yIDJCKiGg0u7NLSUoSHh2PNmjVwc3O74Xs6ne6/f7ITEZGlODXk\nSpWVlQgPD8f06dMxadIkADVH1Xl5efDx8UFubi68vb3r/Nno6Gjz5yEhIQgJCWnyoImI7InBYIDB\nYKj3evWedBQRREZGwsvLC6tWrTJfvnDhQnh5eWHRokWIiYlBcXExTzpaEE86Et0+btWd9Rb2d999\nhxEjRqBPnz7mZY/ly5dj0KBBiIiIwLlz5+Dn54f4+Hh4eHg0KNRSWNiN3nqzz4+IftPowrZEqCXz\nWNiN2nqzz4+IftPop/UREZFtYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGRIljY\nRESKYGETUZO5u3ua37VTyw93d8/mnppN4UvTG7blZn/pNl+aTrbMnn/3mgNfmk5EpDgWNhGRIljY\nRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgW\nNhGRIljYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmChU1EpAgWNhGRIljYRESKYGETESmC\nhU1EpIh6C3vWrFnQ6/Xo3bu3+TKj0YjQ0FD4+/tj7NixKC4utuggiYioAYU9c+ZM7Nix44bLYmJi\nEBoaipMnT2L06NGIiYmx2ACJiKiGTkSkvitlZWVhwoQJOHr0KACgZ8+e2LNnD/R6PfLy8hASEoL0\n9PSbN67ToQGb14xOpwNgiTzrzqPOEVhsboAtzI/UZs+/e83hVt3ZqDXs/Px86PV6AIBer0d+fn7T\nRkdERPVyauoGdDrdfx9d6xYdHW3+PCQkBCEhIU2NJCKyKwaDAQaDod7rNXpJxGAwwMfHB7m5uRg1\nahSXRCyMSyJky+z5d685aLokMnHiRMTFxQEA4uLiMGnSpKaNjoiI6lXvEfbUqVOxZ88eFBQUQK/X\n4+WXX8aDDz6IiIgInDt3Dn5+foiPj4eHh8fNG+cRtnYj4BE22TB7/t1rDrfqzgYtiWgdask8e73T\nsLDJltnz715z0HRJhIiIrI+FTUSkCBY2EZEiWNhERIpgYRMRKYKFTUSkCBY2EZEiWNhERIpgYRMR\nKYKFTUSkCBY2EZEiWNhERIpgYRMRKYKF3QTu7p7m/7ij5Ye7u2dzT83qLLUvb9f9SfaJb6/asC3X\n/d60Vsyz97dXtff52Tu+vaq2+PaqRESKY2ETESmChU1EpAgWNpEV8AQ1aYEnHRu2ZZ50tDDOr9Fb\nbva5AfY/P2vjSUciIsWxsImIFMHCJiJSBAubiEgRLGwiIkWwsImIFMHCJiJShJOlA2qen6ktN7e2\nuHzZqPl2iYhsmcUL2xJPpi8p0f5BgIjI1nFJhIhIESxsIiJFsLDptsQ3Y1Lb7Xr7WfzNn+z1zZis\nncc3R2rS1q24P5v/vtIcuD81HgXf/ImISG0sbCIiRTSpsHfs2IGePXuie/fuWLFihVZjIhvA/2Ku\nttt1jdfeNXoNu7q6Gj169MCuXbvg6+uLgQMH4qOPPsI999zz28btfF3LntewmaddFvPUz7M2zdew\nU1JScPfdd8PPzw8tWrTAlClT8NlnnzVpkEREdGuNLuycnBx07tzZ/HWnTp2Qk5OjyaCIiOhmjS5s\nS7xHCBER3Vqj30vE19cX2dnZ5q+zs7PRqVOnG65z1113ITPTMsV+6wcMe86z3IMk8yydxTz186yn\nb9++dV7e6JOOVVVV6NGjB3bv3o2OHTti0KBBN510JCIi7TT6CNvJyQnr1q1DWFgYqqur8dhjj7Gs\niYgsyKIvTSciIu3wlY5kU65du9bcQyCyWVYvbGsd0JtMJqtlXb16FYB15mbtp05ac24AsG7dOhQX\nF1slC7D+/CoqKqyad/ToUavmXbhwAXl5eaisrLRKnrX3Z3NzjI6OjrZ0yJo1a7Bv3z7ce++95jOw\nImKRs7G//PILrl69Cg8PD+h0OlRXV8PBwXKPS1988QVWrVqFf/7zn7h06RLy8/Oh1+vRunVrzbM+\n/PBDLF68GBUVFfD19YWbm5vmGdez5twAYPPmzZg3bx5+/PFHjB49Gu7u7hbJqWXt+X311Vd46623\n8Pbbb+PixYsoLS2Fl5cXWrVqZZG8zZs3Y8OGDZgyZYpVnvnwn//8By+++CK2bNkCLy8v9OzZ06J5\n1t6ftsDihV1YWIj77rsPGRkZ2Lp1Kzp06IDu3btbpLhFBMOHD8euXbtQVFSEfv36wdnZGUDNEbfW\nd9ravH/84x8ICAhAXl4ejh8/jhMnTsDHxweenp6azu+nn37CoUOH4O3tjS+//BIiAgcHB6SlpaFr\n164W2ZfWmpuIYP78+YiPj0dZWRkOHz6M4cOHw8HBwS5uOxHBuHHjEBUVhR49eiA+Ph47d+5ERUUF\n/Pz84OrqqknO9Xnz589HdHQ0OnXqhLNnz8JgMKCwsBDXrl2Dp6e27wkiInjooYewdu1a9OzZE9u2\nbUNlZSV++eUXFBQUwM/PT/M8a+5PW2Hxwk5NTUVAQAA++ugjVFVVYenSpdi7dy8CAwPh7e2N9evX\no7y8XJMbdO3atcjIyMDy5cuxa9cuvP/++3ByckJAQAB0Oh2ysrLg5ORkLvGmOnHiBNLT0/H888+j\nV69eGDJkCBwdHfHzzz/j4MGDCA0N1bRounfvjqSkJHTr1g1dunRBSkoKoqKi0LZtW4SEhGiaZe25\nvfTSS6iursb8+fPh7OyMf//73zh//jxGjhxpkaNDa88vKSkJWVlZWLhwIXr06IFhw4YhOTkZ5eXl\nyMjIwMiRIzV9gHj++eeRnJyM1atXAwAefPBBnDp1Cnv37sWJEycwcuRItGjRQpMsoOavlczMTCxa\ntAgBAQGYM2cO3NzckJ2djSNHjiAoKAht2rTRLM/a+9NmiIWZTCY5ffq0mEwmERGpqKiQp59+Wrp3\n7y6RkZHSuXNnKS4u1iTr3LlzsmfPHikvL5esrCx59913JTw8XObMmSNJSUkyePBgycvL0yRLRKSk\npETGjx8vs2fPloyMDPPl2dnZMnLkSElMTNQsq3b//fDDDxIbGysiIq+++qr4+vrKE088IUuXLpWr\nV69qlldaWir33XefVeZWWVkpcXFxUlBQYL7s4MGDMnToUHnzzTelqqpKqqqqNMsTqZnf/fffb5X5\niYhcvHhRJk6cKOvWrZP09HRZtWqVzJ07VwoKCmTUqFFy+vRpTfOOHj0qffv2lT59+sigQYPklVde\nERGR3NxcCQ0NlU8//VTTvLy8PBk/fryMGzdO7r//fpk7d66IiBiNRpk5c6asWLHCfB/WwsWLF+WB\nBx6w2v60FRY9wr58+TJatWqFtm3bQkRQWVkJZ2dnjBs3DjNmzMC0adMQFRWFUaNGNTnr6tWr8PLy\nQpcuXeDk5AQPDw/07NkTvXr1QnV1NWbNmoXQ0FBMmTJFg5nVnOxwcXFBaGgojh49ioMHD8JoNKJ9\n+/bo0KED9u7dC2dnZwwZMkSTPJ1OB5PJhNatW2P16tXw9PREbGwsli9fjuDgYHh7eyMgIECTLABw\ndnbGhAkTzMswRUVFFpubg4MDevXqhTvuuANVVVUAgI4dO6JDhw749NNP0blzZ3Tp0kWTrFrOzs4I\nCwtDWloaDh06ZNHbDgBcXFzQunVrHDlyBG+99RauXLmCJUuWoGPHjkhKSoKLiwt69+6tWZ63tzfm\nzJmDVq1aYf/+/Vi5ciVcXV3h6uqKffv2oUWLFprOz9XVFYMHD0a7du3g4+OD3NxchISEwMPDA8nJ\nyWjZsiVGjhypWZ6Liws8PT3x448/4u2330ZZWZlF96etsNjzsBMTE7FlyxacOHECK1euvOnGSk9P\nx6xZs/D99983OevAgQM4cuQIZs2aBaBmvbp27dPBwQEGgwF/+tOfcOrUKU1OSNTmzZw5EzqdDgcP\nHoTBYMDZs2dx7Ngx6PV6HDp0CKmpqXBxcWlS1uHDh1FaWoqgoCDzutz+/fuxePFi+Pr64oMPPmjy\nfK6Xm5uLX3/9FWlpaYiMjMS5c+ewbds282Vazg2oe37XW7ZsGTZv3oxffvkFjo6OmuX16dMH7u7u\nSEtLQ1JSEnJzc3HkyBHN51dZWXnD0kN2djZcXV3h6OgId3d3GAwGPPXUU0hNTdXkvll7+x06dAh/\n/vOf4eR042vjUlNTMWfOHOzbt0/TvPT0dEyZMgWOjo7Izc01L9WZTCbs3r0bP/zwgyYnc2vzTpw4\ngalTp+LKlSuorKy02P60NRYpbBFBv379sGbNGhw7dgwHDhxAeHg4XF1d4ePjg+7du+PUqVMAgLvv\nvrvJeSNHjsT8+fMRHh4OACgpKbnhGRSbNm1CdXU1Zs+e3eSsuvLkv2tlFy5cwPnz55GZmYng4GB0\n69atyVmTJ0/GwYMH8dxzz2H06NG45557cOXKFXzwwQeYNGkSvL29zQ9MWggLC0OnTp1QVFQEZ2dn\nbNq0CQ4ODsjIyICI4MSJExgwYIAmcwNunN+YMWPQs2dPmEwmVFVVmc815OXlwcfHR9O8v//97wgL\nC4O/vz/Kyspw/PhxODo6IiMjQ7PbDgCioqIwZMgQBAcH3zQHk8mE1157DR4eHnjyySc1yau9/YqL\ni+Hg4IBVq1bB19cXOp0O5eXliIqKQrdu3TBv3jxN84qKiuDk5ITVq1dDr9dj9+7dyMzMRElJCYYO\nHYphw4Zpmmc0GtGyZUu8+eabaN++PYCa/bl8+XK0adNGs/1pcyyxzhIVFSVTpkwRkZo1rDZt2sjU\nqVPlwQcflAULFkh5eblmWevXr5f/+7//M3+9dOlSmTp1qvTr10+Sk5M1y7lV3ksvvSQREREyaNAg\nMRgMmuclJiZKnz595OGHH5YxY8bI1q1b5dlnn5XU1FTNs9atWycTJ04UEZGCggL5y1/+Itu2bdM8\n53r/O79t27bJvHnz5Ntvv7V43ujRo2Xr1q2yYMEC+f777zXPWrduneh0Opk9e7a8/vrr8v3330tJ\nSYnmOdfnXX/7zZ49+6a16qKiIovmffLJJ5ptvyF5Wq/F2zqLFPY333wjv/76q4iIxMTEyPz580VE\n5Ndff5WwsDD57LPPNMkxmUwyZ84cmTBhgiQnJ8v8+fNl6tSpcvjwYdmwYYMEBwfL+fPnNcmqL2/9\n+vUSHBwsOTk5muXVWr16tWRmZkpaWprce++9otPpZPv27VJWVqZZhslkkpkzZ97woPP+++/LhAkT\nzF+vW7dOKisrNcus9b/zc3BwkO3bt8uVK1c0z6orr3Z/lpaWapqzatUqSUxMlP3798u8efPkscce\nk3/961+SkZEhJpNJXnrpJamoqNDkZFxDbr+33npLs9uvofeXa9euWS3vzTff1CzPVmle2LXPVKi9\nE5aWlt5wh5w3b56sXbtWs7ycnBx57733ZNasWdKxY8cbnnEyY8YM+eqrrzTLao48EZGNGzdKRESE\niIgMHTpUpk+fLgMGDJCNGzdqmnPy5Ek5fvy4VFdXi0jNUUxAQIBcvnxZXnvtNXn00Uc1zatlrflZ\nO6+iokIKCwvNX2/ZskVmzJghUVFR8sADD8jo0aM1zavv9ps2bZpV87S+vzTX/dOWaFrYqampsmnT\nJvPXtUVd+6iempoq/fv312RJ5PTp0/Lqq6+an75z7Ngx+frrr83fv3TpkvTp00fOnTvX5Cxr550+\nfVpefvll+fnnn82XLV26VMLDw2XEiBEiUvM0Kq3+vD19+rQsW7ZMMjMzzZfVPo3u1VdflWXLlklw\ncLBkZWVplmft+Vk775VXXpH09PSbvldZWSlRUVHi6uqq6f609u1nz3m2TNPCHjFihGzdutX89aVL\nl8yfX7hwQRYsWCAbNmzQJGv27Nmi0+lk+vTpsmrVKvMSjEjNjfnYY4/J4sWLNcmydl5tVmRkpMTE\nxEhOTo4cO3ZMhgwZInv27NEko668GTNmyIoVK2644+/YsUN0Op0sX75c8zxrz6859ufrr79+U5H8\n7W9/Mz9PWes8a95+9ppnyzR7HvaGDRuQlpaG2NhYAEB0dDTi4uKwbNky9OjRA4GBgejRowdCQ0O1\niMOwYcOQl5eHO++8E87Ozvjss89gNBrh7+8PnU6HM2fOYMmSJZq90smaebVZnTt3houLCxISEtC+\nfXssWbIE/v7+Gsym7rw777wTd9xxBz7//HMUFhaie/fu6NGjBzIzM7Fy5UrN96W159cc+9PFxQWf\nf/65+b7SokULdO3aFTNmzLBInjVvP3vNs2WaFLaI4N1330VVVRV8fHywevVqZGdnY9GiRfD29saK\nFSswefJk+Pr6ajDkGq1atUJFRQXi4uLQqVMnBAYGIiUlBXv37oWLiwumTp2q6Q1ozbzrs3x9fREY\nGIjk5GTs27cPlZWVmpdMXXmpqanYu3cv3NzcsGDBAovty+aan7Xzau8r1dXVmr5g5VZ51r797CnP\nlmlS2DqdDv3790fr1q2RmJiIxMREfPXVV+jSpQsGDhyIffv2wdPTE927d2/ygE+dOoWCggK0a9cO\nvXr1wvDhw2E0GjFu3Dj07dsXZ86cwfnz5zFixIgmZ1k77/eygoKCkJWVhYKCAgwfPlyDmTVsbmfP\nnrXKvrT2/Jpjf2ZlZeHChQt2efvZQ54KmvzCmTNnzuDDDz/EtGnT0LVrVxw/fhzZ2dkYO3YsgJqX\npw8fPhzbt29H586dmzTYkpIS+Pr64o477sC0adNw1113IS0tDWfOnEG7du2wadMmlJSUoEWLFvDw\n8GhSlrXz7HluzGOereeposlH2IsXL0ZsbCyMRiPOnj2LAQMGICgoCABQXV2Np556CkFBQXjooYea\nPNhr166hqKgIJSUlSE9PR3h4OO666y6cPn0aHTp0QO/evdGxY0fNXpJqzTx7nhvzmGfrecpo6llL\no9EokZGRsmTJEnnjjTdk7ty58t5770lZWZlcu3ZNVq5cqem7dInUPB/zr3/9q/Tr1++GF+FondMc\nefY8N+Yxz9bzbF2TC9tkMslHH30k/fv3lxdeeEE2btwoc+fOlWeeeUa++eYbDYZYo/YFKtu2bZP4\n+Hiprq6WxMREGTlypNx3332yY8cOzbKsnWfPc2Me82w9TyWNXsM+deoUTCaT+Qz7sWPHYDAYMH78\neJSUlOCTTz6Bo6MjXnzxxSb/FbB7927Ex8ejsLAQP//8M1JTU+Hm5gYRQVlZGTZu3Ijdu3cjMTGx\nyVnWzrPnuTGPebaep5zGtPzly5fFzc1NfHx85Nlnn5X169fLE088IWPHjpVp06ZJeXm5XLhwQbNX\njvXv31+SkpJk6NChotfrJTY29oY3609NTdX0fTWsmWfPc2Me82w9TzVO9Vd63WbOnIn9+/dj586d\nWL58OXx9fbFlyxb4+PjgwoULuPPOOzV5QHnttdfg6+uLMWPGoLCwEJ06dcLatWsxYsQIvPDCC7h0\n6RK+/fZbvPPOO8rl2fPcmMc8W89TUlPa3tInBPLy8sTb21tmz54tO3fuvOFfSH388ccycOBAGThw\noKSlpSmXZ89zYx7zbD1PVX94DfvSpUto06YNPvnkE1RXVyM8PBxffPEFVq5cCRcXFyxYsABhYWGa\nPJjMmjULLVu2xMCBA/HTTz/Bx8cH9957L4YPH25+Ok9OTo5mr6C0Zp49z415zLP1PFX9ocK25gmB\nrKwsREREICUlBQCQnJyMvXv3wmg0olu3bhgyZIimL/G1Zp49z415zLP1PJX9ocIODg7GihUrsHTp\nUmRmZuK5557Dk08+iZYtWwIAfvrpJwQGBmryv9sAwGg0wtPT0/z15cuX8fXXX+PAgQMoLCzE008/\njcDAQE2yrJ1nz3NjHvNsPU9ZDV07WbZsmfm/O3z88cfy3Xffmf81VkJCgsTFxcnjjz+u/aKN1KyJ\nX78ufvr0afn4448tkmXtPHueG/OYZ+t5qmnQEXZ+fj769OmDSZMm4eGHH0b//v3h5eUFANiyZQtW\nrlwJoOaf3VryX8vLf//ZrbVYM8+e58Y85tl6nioaVNi2dkLAnu889jw35jHP1vNsXb2FzRMCRES2\noUFH2Dzn2SM7AAAAbUlEQVQhQETU/P7Qs0Rqr1r7J8qZM2eQkpKCRx55xDKjIyIis0a9+RPXlYiI\nrM+hMT9UW9aN6HoiImqkRhV2LR5lExFZT5MKm4iIrIeFTUSkCBY2EZEiWNhERIpgYRMRKYKFTUSk\nCBY2EZEi/h8oTXmmh4vCsgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f6b8a49fe50>"
       ]
      }
     ],
     "prompt_number": 389
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "What about when only two make it through, do you have a over-represented pairs?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def extract_doubles(rec):\n",
      "    if rec.num_hom_alt == 2:\n",
      "        return \"_\".join([s.sample for s in rec.get_hom_alts()])\n",
      "    if rec.num_hom_ref == 2:\n",
      "        return \"_\".join([s.sample for s in rec.get_hom_refs()])\n",
      "\n",
      "doubles = Counter([extract_doubles(r) for r in hets if extract_doubles(r)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 390
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figure()\n",
      "N = len(doubles.keys())\n",
      "bar(np.arange(N), doubles.values())\n",
      "xticks(np.arange(N)-0.5, doubles.keys(), rotation=90)\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAE/CAYAAABMwNhyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9YVGXeP/D3IJT5a/0NOkM7KiCOjIgmmRs5haCyl8Rm\nUei1otRu5a5LbanVtt+wWqG0Zyt72Ge7Hn/vs6JbG5BbLG4G1ZqgWFri9UA6GjMgK6KG+QOVz/cP\n13lEZhAOw5kzx/fruuYKzpz33Pe555wPpzO3cwwiIiAiIt0J8HUHiIioe7DAExHpFAs8EZFOscAT\nEekUCzwRkU6xwBMR6VS7BT4jIwPBwcGwWq2uZeXl5YiNjUVMTAwmTZqEXbt2uZ7Lzs5GeHg4IiMj\nUVxc7FpeUVEBq9WK8PBwZGZmdsNmEBFRG9KOTz75RPbs2SNRUVGuZVOnTpWioiIREfnggw/EZrOJ\niMj+/fslOjpampubxW63y6hRo6SlpUVERCZNmiRlZWUiIjJz5kz58MMP22uWiIi8oN0z+Li4OAwY\nMKDVsmHDhuHUqVMAgJMnT8JoNAIACgoKkJaWhqCgIJjNZoSFhaGsrAx1dXVoampCbGwsAGDevHnI\nz8/vjr9VRER0lcDOBnJycnDnnXfi6aefRktLCz7//HMAQG1tLSZPnuxaz2Qywel0IigoCCaTybXc\naDTC6XR6oetERNSeTn/I+vDDD+PNN9/Et99+i9///vfIyMjojn4REVEXdfoMvry8HP/4xz8AAPff\nfz8eeeQRAJfPzGtqalzrORwOmEwmGI1GOByOVsuvXNa5VlhYGA4ePNjZLhER3dCio6Px5Zdftlne\n6TP4sLAwlJaWAgC2b9+OiIgIAEBycjLy8vLQ3NwMu92O6upqxMbGIiQkBP369UNZWRlEBBs3bkRK\nSorb1z548CBEpFOPF1544YbPaL1/estovX/MqN+Wrx979+51W1PbPYNPS0tDaWkpGhoaEBoaihdf\nfBFvv/02fvGLX+D8+fO45ZZb8PbbbwMALBYLUlNTYbFYEBgYiNzcXBgMBgBAbm4u5s+fj7NnzyIp\nKQkzZszo7N8VIiLqpHYL/KZNm9wuLysrc7v8ueeew3PPPddm+cSJE/HVV18p6B4RESnVIysrK8vX\nnbhi2bJlUNIds9l8w2fUbIsZddtiRvvvka95qp0GEdHMDT8MBgM01B0iIr/gqXbyu2iIiHSKBZ6I\nSKdY4HWuX7+BMBgMHXr06zfQ190lIi/iNXiduzxVtaNjyvEn8ke8Bk9EdINhgSci0ikWeCIinWKB\nJyLSKRZ4IiKdYoEnItIpFngiIp1igSci0ikWeCIinWKBJyLSqXYLfEZGBoKDg2G1WlstX7VqFcaM\nGYOoqCgsXbrUtTw7Oxvh4eGIjIxEcXGxa3lFRQWsVivCw8ORmZnp5U0gIiK3pB2ffPKJ7NmzR6Ki\nolzLtm/fLtOmTZPm5mYREfnXv/4lIiL79++X6OhoaW5uFrvdLqNGjZKWlhYREZk0aZKUlZWJiMjM\nmTPlww8/dNvedbpDCgAQQDr44PgT+SNPx267Z/BxcXEYMGBAq2V/+MMf8OyzzyIoKAgAMGTIEABA\nQUEB0tLSEBQUBLPZjLCwMJSVlaGurg5NTU2IjY0FAMybNw/5+fne/jtFRETX6PQ1+OrqanzyySeY\nPHkybDYbdu/eDQCora2FyWRyrWcymeB0OtssNxqNcDqdXug6ERG1p92bbrtz8eJFnDhxAjt37sSu\nXbuQmpqKQ4cOea1DV99X0GazwWazee21iYj0oKSkBCUlJdddr9MF3mQy4b777gMATJo0CQEBAWho\naIDRaERNTY1rPYfDAZPJBKPRCIfD0Wq50Wj0+Poaugc4EZEmXXvyu2zZMrfrdfoSTUpKCrZv3w4A\nqKqqQnNzMwYPHozk5GTk5eWhubkZdrsd1dXViI2NRUhICPr164eysjKICDZu3IiUlBRlW0VERB3W\n7hl8WloaSktLcfz4cYSGhuLFF19ERkYGMjIyYLVacdNNN2HDhg0AAIvFgtTUVFgsFgQGBiI3N/ff\ndxMCcnNzMX/+fJw9exZJSUmYMWNG928ZEdENjrfs0zneso9I/3jLPiKiGwwLPBGRTrHAExHpFAs8\nEZFOscATEekUCzwRkU6xwBMR6RQLPBGRTrHAExHpFAs8EZFOscATEekUCzwRkU6xwBMR6RQLPBGR\nTrHAExHpFAs8EZFOtVvgMzIyEBwcDKvV2ua51157DQEBAWhsbHQty87ORnh4OCIjI1FcXOxaXlFR\nAavVivDwcGRmZnqx+0RE5Em7BX7BggUoKipqs7ympgbbtm3DD3/4Q9eyyspKbN68GZWVlSgqKsLC\nhQtddxh5/PHHsXr1alRXV6O6utrtaxIRkXe1W+Dj4uIwYMCANst//etf49VXX221rKCgAGlpaQgK\nCoLZbEZYWBjKyspQV1eHpqYmxMbGAgDmzZuH/Px8L24CERG50+lr8AUFBTCZTBg3blyr5bW1tTCZ\nTK7fTSYTnE5nm+VGoxFOp7MLXSYioo4I7MzKZ86cwfLly7Ft2zbXMm/fpDkrK8v1s81mg81m8+rr\nExH5u5KSEpSUlFx3vU4V+IMHD+Lw4cOIjo4GADgcDkycOBFlZWUwGo2oqalxretwOGAymWA0GuFw\nOFotNxqNHtu4usATEVFb1578Llu2zO16nbpEY7VaUV9fD7vdDrvdDpPJhD179iA4OBjJycnIy8tD\nc3Mz7HY7qqurERsbi5CQEPTr1w9lZWUQEWzcuBEpKSld2jgiIrq+dgt8WloapkyZgqqqKoSGhmLt\n2rWtnjcYDK6fLRYLUlNTYbFYMHPmTOTm5rqez83NxSOPPILw8HCEhYVhxowZ3bApRER0NYN4+yJ6\nFxgMBq9f07/RXf4j29Ex5fgT+SNPtZP/kpWISKdY4ImIdIoFnohIp1jgiYh0igWeiEinWOCJiHSK\nBZ6ISKdY4ImIdIoFnohIp1jgiYh0igWeiEinWOCJiHSKBZ6ISKdY4ImIdIoFnohIp1jgiYh0qt0C\nn5GRgeDgYFitVteyxYsXY8yYMYiOjsZ9992HU6dOuZ7Lzs5GeHg4IiMjUVxc7FpeUVEBq9WK8PBw\nZGZmdsNmEBHRtdot8AsWLEBRUVGrZYmJidi/fz/27t2LiIgIZGdnAwAqKyuxefNmVFZWoqioCAsX\nLnTdYeTxxx/H6tWrUV1djerq6javSURE3tdugY+Li8OAAQNaLUtISEBAwOXY7bffDofDAQAoKChA\nWloagoKCYDabERYWhrKyMtTV1aGpqQmxsbEAgHnz5iE/P787toWIiK7SpWvwa9asQVJSEgCgtrYW\nJpPJ9ZzJZILT6Wyz3Gg0wul0dqVZIiLqgEClwd/97ne46aabMGfOHG/2B1lZWa6fbTYbbDabV1+f\niMjflZSUoKSk5LrrKSrw69atwwcffICPPvrItcxoNKKmpsb1u8PhgMlkgtFodF3GubLcaDR6fO2r\nCzwREbV17cnvsmXL3K7X6Us0RUVFWLFiBQoKCtCzZ0/X8uTkZOTl5aG5uRl2ux3V1dWIjY1FSEgI\n+vXrh7KyMogINm7ciJSUlM5vERERdUq7Z/BpaWkoLS1FQ0MDQkNDsWzZMmRnZ6O5uRkJCQkAgDvu\nuAO5ubmwWCxITU2FxWJBYGAgcnNzYTAYAAC5ubmYP38+zp49i6SkJMyYMaP7t4yI6AZnkCtzGTXA\nYDBAQ93Rhct/ZDs6phx/In/kqXbyX7ISEekUCzwRkU6xwBMR6RQLPBGRTrHAExHpFAs8EZFOscAT\nEekUCzwRkU6xwBMR6RQLPBGRTrHAExHpFAs8EZFOscATEekUCzwRkU6xwBMR6VS7BT4jIwPBwcGw\nWq2uZY2NjUhISEBERAQSExNx8uRJ13PZ2dkIDw9HZGQkiouLXcsrKipgtVoRHh6OzMzMbtgMIiK6\nVrsFfsGCBSgqKmq1LCcnBwkJCaiqqkJ8fDxycnIAAJWVldi8eTMqKytRVFSEhQsXur6A/vHHH8fq\n1atRXV2N6urqNq9JRETe126Bj4uLw4ABA1otKywsRHp6OgAgPT0d+fn5AICCggKkpaUhKCgIZrMZ\nYWFhKCsrQ11dHZqamhAbGwsAmDdvnitDRETdp9PX4Ovr6xEcHAwACA4ORn19PQCgtrYWJpPJtZ7J\nZILT6Wyz3Gg0wul0drXfRER0HV36kNVgMLhurE1ERNoS2NlAcHAwjh49ipCQENTV1WHo0KEALp+Z\n19TUuNZzOBwwmUwwGo1wOBytlhuNRo+vn5WV5frZZrPBZrN1totERLpWUlKCkpKS668o12G32yUq\nKsr1++LFiyUnJ0dERLKzs2Xp0qUiIrJ//36Jjo6W8+fPy6FDh2TkyJHS0tIiIiKxsbGyc+dOaWlp\nkZkzZ8qHH37otq0OdIc6CYAA0sEHx5/IH3k6dts9g09LS0NpaSkaGhoQGhqKF198Ec888wxSU1Ox\nevVqmM1mbNmyBQBgsViQmpoKi8WCwMBA5Obmui7f5ObmYv78+Th79iySkpIwY8aMLv31IiKi6zP8\nu/prgsFggIa6owuX/8h2dEw5/kT+yFPt5L9kJSLSKRZ4IiKdYoEnItIpFngiIp1igSci0ikWeCIi\nnWKBJyLSKRZ4IiKdYoEnItIpFngiIp1igSci0ikWeCIinWKBJyLSKRZ4IiKdYoEnItIpFngiIp1S\nXOCzs7MxduxYWK1WzJkzB+fPn0djYyMSEhIQERGBxMREnDx5stX64eHhiIyMRHFxsVc6T0REnim6\no9Phw4dxzz334MCBA7j55pvx4IMPIikpCfv378fgwYOxZMkSvPLKKzhx4gRycnJQWVmJOXPmYNeu\nXXA6nZg2bRqqqqoQEND67wvv6OR9vKMTkf559Y5O/fr1Q1BQEM6cOYOLFy/izJkzGD58OAoLC5Ge\nng4ASE9PR35+PgCgoKAAaWlpCAoKgtlsRlhYGMrLy7uwOUREdD2KCvzAgQPx1FNP4dZbb8Xw4cPR\nv39/JCQkoL6+HsHBwQCA4OBg1NfXAwBqa2thMplceZPJBKfT6YXuExGRJ4oK/MGDB/H666/j8OHD\nqK2txenTp/GnP/2p1ToGg+Hflwfca+85IiLqukAlod27d2PKlCkYNGgQAOC+++7D559/jpCQEBw9\nehQhISGoq6vD0KFDAQBGoxE1NTWuvMPhgNFodPvaWVlZrp9tNhtsNpuSLhIR6VZJSQlKSkquu56i\nD1n37t2LuXPnYteuXejZsyfmz5+P2NhYHDlyBIMGDcLSpUuRk5ODkydPtvqQtby83PUh6zfffNPm\nLJ4fsnofP2Ql0j9PtVPRGXx0dDTmzZuH2267DQEBAZgwYQJ+/vOfo6mpCampqVi9ejXMZjO2bNkC\nALBYLEhNTYXFYkFgYCByc3N5iYaIqJspOoPvLjyD9z6ewRPpn1enSRIRkfaxwBMR6RQLPBGRTrHA\nExHpFAs8EZFOscATEekUCzwRkU6xwBMR6RQLPBGRTrHAExHpFAs8EZFOscATEekUCzwRkU6xwBMR\n6RQLPBGRTrHAExHplOICf/LkSdx///0YM2YMLBYLysrK0NjYiISEBERERCAxMREnT550rZ+dnY3w\n8HBERkaiuLjYK50nIiLPFBf4zMxMJCUl4cCBA9i3bx8iIyORk5ODhIQEVFVVIT4+Hjk5OQCAyspK\nbN68GZWVlSgqKsLChQvR0tLitY0gIqK2FBX4U6dO4dNPP0VGRgYAIDAwED/4wQ9QWFiI9PR0AEB6\nejry8/MBAAUFBUhLS0NQUBDMZjPCwsJQXl7upU0gIiJ3FBV4u92OIUOGYMGCBZgwYQJ+9rOf4fvv\nv0d9fT2Cg4MBAMHBwaivrwcA1NbWwmQyufImkwlOp9ML3SciIk8ClYQuXryIPXv24K233sKkSZPw\nxBNPuC7HXGEwGP59w2f3PD2XlZXl+tlms8FmsynpIhGRbpWUlKCkpOS66ykq8CaTCSaTCZMmTQIA\n3H///cjOzkZISAiOHj2KkJAQ1NXVYejQoQAAo9GImpoaV97hcMBoNLp97asLPBERtXXtye+yZcvc\nrqfoEk1ISAhCQ0NRVVUFAPjHP/6BsWPHYtasWVi/fj0AYP369UhJSQEAJCcnIy8vD83NzbDb7aiu\nrkZsbKySpomIqIMUncEDwKpVqzB37lw0Nzdj1KhRWLt2LS5duoTU1FSsXr0aZrMZW7ZsAQBYLBak\npqbCYrEgMDAQubm57V6+ISKirjOIiPi6E1cYDAZoqDu6cPkPaUfHlONP5I881U7+S1YiIp1igSci\n0ikWeCIinWKBJyLSKRZ4IiKdYoEnItIpFngiIp1igSci0ikWeCIinWKBJyLSKRZ4IiKdYoEnItIp\nFngiIp1igSci0ikWeCIinWKBJyLSqS4V+EuXLiEmJgazZs0CADQ2NiIhIQERERFITEzEyZMnXetm\nZ2cjPDwckZGRKC4u7lqviYjourpU4N944w1YLBbX7fdycnKQkJCAqqoqxMfHIycnBwBQWVmJzZs3\no7KyEkVFRVi4cCFaWlq63nsiIvJIcYF3OBz44IMP8Mgjj7huFVVYWIj09HQAQHp6OvLz8wEABQUF\nSEtLQ1BQEMxmM8LCwlBeXu6F7hMRkSeKC/yTTz6JFStWICDg/16ivr4ewcHBAIDg4GDU19cDAGpr\na2EymVzrmUwmOJ1OpU0TEVEHBCoJbd26FUOHDkVMTAxKSkrcrmMwGFyXbjw9705WVpbrZ5vNBpvN\npqSLRES6VVJS4rH2Xk1Rgd+xYwcKCwvxwQcf4Ny5c/juu+/w05/+FMHBwTh69ChCQkJQV1eHoUOH\nAgCMRiNqampceYfDAaPR6Pa1ry7wRETU1rUnv8uWLXO7nqJLNMuXL0dNTQ3sdjvy8vJwzz33YOPG\njUhOTsb69esBAOvXr0dKSgoAIDk5GXl5eWhubobdbkd1dTViY2OVNE1ERB2k6Az+WlcutzzzzDNI\nTU3F6tWrYTabsWXLFgCAxWJBamoqLBYLAgMDkZub2+7lGyIi6jqDXJkCowEGgwEa6o4uXP5D2tEx\n5fgT+SNPtZP/kpWISKdY4ImIdIoFnohIp1jgiYh0igWeiEinWOCJiHSKBZ6ISKdY4ImIdIoFnohI\np1jgiYh0igWeiEinWOCJiHSKBZ6ISKdY4ImIdIoFnohIp1jgiYh0SlGBr6mpwd13342xY8ciKioK\nb775JgCgsbERCQkJiIiIQGJiIk6ePOnKZGdnIzw8HJGRkSguLvZO74mIyCNFd3Q6evQojh49ivHj\nx+P06dOYOHEi8vPzsXbtWgwePBhLlizBK6+8ghMnTiAnJweVlZWYM2cOdu3aBafTiWnTpqGqqgoB\nAa3/vvCOTt7HOzoR6Z9X7+gUEhKC8ePHAwD69OmDMWPGwOl0orCwEOnp6QCA9PR05OfnAwAKCgqQ\nlpaGoKAgmM1mhIWFoby8XOm2EBFRB3T5Gvzhw4fxxRdf4Pbbb0d9fT2Cg4MBAMHBwaivrwcA1NbW\nwmQyuTImkwlOp7OrTRMRUTsCuxI+ffo0Zs+ejTfeeAN9+/Zt9ZzBYPj35QH3PD2XlZXl+tlms8Fm\ns3Wli0REulNSUoKSkpLrrqe4wF+4cAGzZ8/GT3/6U6SkpAC4fNZ+9OhRhISEoK6uDkOHDgUAGI1G\n1NTUuLIOhwNGo9Ht615d4ImIqK1rT36XLVvmdj1Fl2hEBA8//DAsFgueeOIJ1/Lk5GSsX78eALB+\n/XpX4U9OTkZeXh6am5tht9tRXV2N2NhYJU0TEVEHKZpF89lnn+Guu+7CuHHjXJdasrOzERsbi9TU\nVHz77bcwm83YsmUL+vfvDwBYvnw51qxZg8DAQLzxxhuYPn16285wFo3XcRYNkf55qp2KCnx3YYH3\nPhZ4Iv3z6jRJIiLSPhZ4IiKdYoEnItIpFngiIp1igSci0ikWeCIinWKBJyLSKRZ4IiKdYoEnItIp\nFngiIp1igSci0ikWeCIinWKBJyLSKRZ4IiKdYoEnItIpVQt8UVERIiMjER4ejldeeUXNpomIbjiq\nFfhLly7hl7/8JYqKilBZWYlNmzbhwIEDXX7djtx4Vu8Z5bnOZ7Q8Dlofb2b0+R5pmWoFvry8HGFh\nYTCbzQgKCsJDDz2EgoKCLr+ulncyre/MLPDaLx7MaP890jLVCrzT6URoaKjrd5PJBKfTqVbzREQ3\nHNUK/JWbcxP5Qr9+A2EwGNo8li1b1ur3fv0GXjdzbU5J5tqcL8egO/rmzfH29bj5NVHJ559/LtOn\nT3f9vnz5csnJyWm1zqhRowSX7xDNBx988MFHBx/R0dFu665BxM2tuLvBxYsXMXr0aHz00UcYPnw4\nYmNjsWnTJowZM0aN5omIbjiBqjUUGIi33noL06dPx6VLl/Dwww+zuBMRdSPVzuCJiEhd/JesREQ6\npdolmmsVFRUhPz/fNVXSaDQiJSUFM2bMUPR6L774Iv7f//t/XmvLm/2zWq346quvuj3T3hhc6557\n7sH27ds9Pv/kk09i9uzZuPPOOzvVB6DzY9fQ0IDBgwe7ft+4cSPKy8thtVrxs5/9zO0MrOPHj+Ot\nt96C0WhERkYGsrOzsWPHDlgsFjz33HMYMGCA27a2b9+Od999FzU1NejRowdGjx6NRx55BGFhYR3e\nvuuNXXu8vZ8q7ZuScfDHY6K9jLdrkBb55BJNZmYmqqurMW/ePBiNRgCAw+HAxo0bERYWhjfffLPT\nrxkaGoqamhqvtKUk8+6777ZZZjAYICJ49NFH0dDQ4JVMezyNgdVqdb3uFVVVVYiIiIDBYMC+ffva\nZIYMGYIf/vCH+Ne//oWHHnoIaWlpiImJuW4flIxdTEwMvvjiCwDAyy+/jE8//RRz5szB+++/j9DQ\nUPz+979vk5k5cybGjRuH7777DgcOHIDVasUDDzyAbdu2Yd++fW7/Ed0zzzyDo0ePIj4+Hvn5+Rgx\nYgQiIiLwhz/8Ac8++yxSU1O9Mnbt8eZ+qrRvSsZBy8eEkkx31CAt8kmBDw8PR3V1dZvlIoLw8HB8\n8803bnN9+/b1+Jpnz57FxYsXvdKWkkxQUBDmzJmDgICANpl33nkHp0+f9kpGyRgkJyejb9++eP75\n59GrVy+ICOLi4vDZZ59BRGA2m9tkrhTdqqoq5OXlYfPmzbh48SLmzJmDtLQ0REREuO2DkrG7usDH\nxMTg008/RZ8+fXDhwgXExMTg66+/bpOJjo7G3r17ISIwGo2ora1t89y1oqKiXK918eJF3HXXXdix\nYwdOnDiBO++8E/v37/fK2Km1nyrpm9Jx0PIxoSSjtAb5HeUz25WLioqSsrKyNst37twpUVFRHnOh\noaFSV1fn9jmTyeS1tpRkYmJiZN++fZ3qm5KMkjEQEXn33XflzjvvlPz8fBERMZvNHtcVERk/fnyb\nZV9++aUsXbpURo4c6TGnZOxGjx4tFRUVsnv3bhk7dmyr58aNG+exnePHj8uRI0ekb9++cujQIRER\nOXbsmMd2xo0bJw0NDSIicvjwYbn99ttdz1ksFo/b1NmxU2s/VdI3EWXjoOVjQklG6Xj7G58U+N27\nd8ukSZMkMjJSpk2bJtOmTZPIyEiJjY2V3bt3e8w999xzbt8UEZHFixd7rS0lmdLSUjl8+LDb58rL\ny72WUTIGVzQ1NckTTzwhycnJMnz48HbXdVfgO0LJ2E2dOlVsNpvr4XQ6ReRysZ44caLbzJo1a2Tg\nwIEyatQo2bp1q4wcOVLi4+PFaDTK+vXr3Wby8vLk1ltvlfj4eDGZTPL++++LiEh9fb2kpaW1u12d\nGTu19lMlfRNRNg5aPiaUZLoy3v7Ep9Mk6+rqWn3AMWzYME21pWb/1PTll19i586deOyxxzyu09TU\n1O6lhuvxxthdunQJ586dQ+/evd0+39zcjMDAQAQEBLiuxY8cORJDhgzx+JrHjx/HoUOHEBYW5vGD\n2PZ0ZOy6oivj1pm+KR0HvR0Tetuea/XIysrK8lXjffv2xfDhw9G3b18cOXIEt9xyC3r27NnhfFNT\nE/bv39+hnJK2OpO5cOECNm3ahNraWoSFhWH9+vX47//+b9TV1SEmJsbtTBAlma6MwRUhISEoKytD\nbGysx3VuvvnmNstyc3MxadKkDrWhZLwvXLiAHj16uH6/Urh79erldv0ePXrg4sWL6NGjB26++WaY\nTCb07t0bDQ0NHjO9evXC0KFD2/zRaC9ztY6M3bW6ez9V0jel46DFY6Irx1FXa5Dm+eJ/Gx5//HHX\nz59++qmEhoaKzWYTo9EoW7du9WpOrUxGRobMnj1bZs2aJampqXLvvffKhg0bJDU1VZ5++mmvZZT0\nbeXKlW0eAwcOlNdee01ee+01r2WU9m/79u1iNBpl4MCBkpCQ4LqeLuL5UpFaGSXjoNY+p/Q9UjIO\nWj4m1DqO/JFPCvzVO9HUqVOloqJCREQOHjwoEyZM8GpOrcyVD6eam5tlwIABcu7cORERuXDhglit\nVq9llPStd+/ekpqaKllZWZKVlSUvvPCC9O/f3/W7tzJK+zdx4kT5+uuvpaWlRf7yl7/IqFGjZMeO\nHW1ezxcZJeOg1j6n9D1SMg5aPibUOo78kc8L/LU7VHsf7inJqZW5+tvcEhMTWz3naSaIkoySvh05\nckRmz54tixcvlu+//15Erj/bQklGaf+uPQi//vpriYiIkPfee8/nGSXjoNY+p/Q9UjIOWj4m1DqO\n/JFPCnzPnj0lKipKoqKipHfv3tLY2CgiIhcvXmwzTa6rObUy06dPl6ampjbLa2trZdKkSV7LKB07\nEZH33ntP7rjjDtmyZUuHCoGSjJL+TZw4sc20wpqaGhk3bpz07t3bp5krOjMOau1zSvomomwctHxM\nqH0c+ROfFHi73d7qcf78eRG5PC3unXfe8WpOrYwnp0+flqNHj3ot09W+NTU1yVNPPSVxcXEd7k9n\nMkr6V1xcLF988UWb5SdOnJCXXnrJp5mrdXQcfLHPdeY9UjIOWj4mlGS8uT1axm+T7CZNTU2orq7G\nyJEj0b9B6s1qAAAQ/UlEQVR//27LEPkLtY4JHkf/xydfNubuOzSuaO87NJTk1MosXLgQubm5AIDP\nPvsMc+bMwahRo1BdXY0//vGP+PGPf+yVjJbHQOv9Y0bdttQ6JtQ6jvyRTwp8jx49YDAYkJaWhlmz\nZrm+R6M7cmplPv/8c9fPzz//PPLz8zFhwgQcOnQIDzzwgNudTElGy2Og9f4xo25bah0Tah1Hfqk7\nr/+0p7KyUn77299KTEyMzJ07V7Zu3SoXLlzolpwaGbVmTqi1PUozWu8fM+q1peVZbyLKx86f+KzA\nX23Tpk0yaNAgefXVV7s9110ZtWdOdPf2eCOj9f4x071taXnWm5Lt8Uc++5DV4XBg8+bN+Otf/4oB\nAwbgwQcfxE9+8hP06dPH6zk1MocPH271+/Dhw3HTTTehoaEBpaWlmD17tlcyWh4Df+gfM+q1pdYx\noeZx5G98UuDvuusunD59GqmpqbjvvvswaNCgVt8XMXDgQK/l1MqoRetjoOX+MaN+W1qlt+3xxCcF\n/sqNCDx9CZDdbvdaTq2MWjMTtDwGWu8fM+q2peXZRErHzt9wHryXjB8/vt1P5a/sUF3NEPkLtY4J\nHkeeBVx/le4THx/foWXeyHV35ssvv8Sf//xnnD59GnPnzsVvfvMb7N+/HyaTyeMOpiSj1vZ0JaNm\nW8xo9z1S65hQ+zjyK2p/qisicubMGWloaBCr1SrHjx93Pex2u4wePdqrObUy1+qumQlaHwMt948Z\n9du6mpZmE3lje/yBTwr866+/LmazWW666SYxm82uh9VqlVWrVnk1p1ZG5PIXNq1cuVKmTJkiP/7x\nj2XDhg1uvwSpKxmtj4GW+8eM+m2pcUwoySjdHn/jkwJ/5bun33jjjW7PqZWJi4uTmJgYyc7Olv/9\n3/+VhoaGVmcG3spoeQzUbIsZ7b9Hah0Tah1H/sgnH7LGxMTgiy++UCWnVkatmQlaHgM122JG+++R\nlmcTKR07f+OT76LRo2v/sUV3ZYj8hVrHBI8jz3xyBt+/f3/ExcW5fc5gMKCwsNBrObUyV8THx+Oj\njz667jKlGa2PgZb7x4z6bQHdf0woyXRle/yJT87ghwwZgqefftrjP0zwZk6tzNmzZ3HmzBkcO3YM\njY2NruXfffcdnE6n1zJaHgOt948ZddtS65hQ6zjyS2pe8L9C6T0PleTUyqg1M0HLY6BmW8xo/z3S\n8mwiPd13tT0+KfA/+clPOrRecXFxl3NqZdSamaDlMVCzLWa0/x5peTaR0rHzN5r4umBPtH6G0t73\nUHdXm9352no8i2TGd++R1rfP16+tBp9+VQEREXUfTpP0ErvdjlmzZrl9ztOn8koyRP5CrWOCx5Fn\nLPBeouYsCCJ/oPXZRDcCTRf4ESNGqJbraqZPnz6YOnVqp/JKMh3lizHQYlvM+O49UuuY0NpxpCU+\n+z74pqYmFBUVweFwICAgAKNHj0ZiYiICAjr3scBzzz2H5cuXez1TWlqKkJAQjB49Gp999hk+//xz\nWCwWt3doB4D77rsPf/3rX6/b9rZt25CQkKA4U1hYiMTERPTs2fO6OU8OHTqEL774AmPHjkVkZKTX\nMydPnkRRUZFrDrLJZML06dPRv39/jxkl+4Namc7uC1ofA6X962xGrWNCSeZaSo4Jf+CTAr9lyxas\nXLkS48aNw8cff4w77rgDIoJ9+/bhf/7nfzBu3Di3uUWLFrVZtmHDBsybNw8GgwFvvvmmVzKZmZnY\ntWsXLly4gBkzZuCjjz7CzJkzUVpaivHjx2PlypUKtvqyrn4PyC233IJevXohKSkJaWlpmD59Onr0\n6NFuPiUlBfn5+QCAgoICPPHEE7DZbPjnP/+JZ599FgsWLPBKBrg8tsuWLUNCQgJMJhMAoKamBtu2\nbcMLL7yA9PT0Nhkl+4NaGSX7gpbHQGn/lGQ6yhff3aN0//Y7Ppi5I1FRUfL999+LiMixY8ckISFB\nRET27t0rd9xxh8ec0WiUOXPmyLp162TdunWydu1aGTx4sOt3b2XGjBkjly5dktOnT8sPfvADOX36\ntIiINDc3i8Vi6cqme2XqWWNjo/zxj3+Uu+++W4YMGSKPPvqolJSUdCg/efJkOXTokIhcHnur1eq1\njIhIeHi4nDhxos3yxsZGCQsLc5tRsj+olVGyL2h5DJT2T0mmo3w9rbkz+7e/8dk0ySuXGHr37o1j\nx44BAMaNG4dTp055zFRWVmLw4MEoKipCQkIC5s+fjz59+iA9Pd3jGYSSjMFggMFgQI8ePVw/A0BA\nQIAmPrQZMGAAfv7zn2P79u3Yu3cvxowZg6VLlyI0NPS62ebmZtd1xcGDB3fokpiSzLWuN25K9gc1\nMt7cF7QyBkr7562M1nhj/9Yqn3zImpSUhBkzZuCuu+5CUVERHnjgAQDA8ePH283169cPb7zxBioq\nKjB37lwkJSWhpaXF65n4+HjExcWhubkZv/jFL5CQkOD633JP1/B8ZdiwYcjMzERmZqbHb9Xbt28f\n+vbtCwA4d+4c6urqMGzYMJw/f97jWCjJAMBvfvMbTJw4EYmJia3+V764uBi//e1v3WaU7A9qZZTs\nC1oeA6X9U5LRMqX7t7/x2Yesf/vb33DgwAFER0e7DpSWlhY0Nzd36APElpYW5ObmYufOnfjTn/7U\noTY7mhERlJaWYujQobBYLPjkk0+wc+dOREZGIjk5uWMb6EFHPxDylPn4449x9913d6kPV5w8eRKV\nlZWYMmWKVzONjY34+9//jtraWgCA0WjE9OnTMWDAAI8ZJfuDGhml+4K7MUhMTMTAgQN9Pgae+ne9\n90jJNnVEV48Jb2aUHBNa5rMC7y0VFRWYOHFit2c6q6OfyiuZzXCtPXv2YMKECV3uM1FnKDmOOjKD\nTc3ZRHqnua23Wq0en9uzZw/27NmDiooK13/vvfde13JvZb799ls89NBDuPPOO7F8+XJcuHDB9VxK\nSorbzNXLCwoKEB8fj61btyI5ORlr1651m9mwYQMmTpyIkpISnD17FmfPnsX27dsxYcIErF+/vsPb\nM2vWLK9vj5JMV3KetLc/qJFRsj1r1qxx/exwOBAfH4/+/ftjypQpqKqq8lo7SrZHaVtX73MdPY4W\nLVrU5vGf//mfWLRoEX71q1+5zWzZsgXx8fH4+9//jrfeegu7du3Cxo0bER0djX379nkt4+3x1iqf\nXIN/99132ywzGAwQEdTV1XnM3XbbbZg8eTJuvvlm17Ljx4/jqaeeAnD58oU3MhkZGbj//vtx++23\nY/Xq1Zg6dSoKCwsxePBgHDlyxG3frl6ek5OD7du3Y8SIEWhoaMA999zjdtrVyy+/jIqKijZn6ydO\nnEBsbKzbD4HdbU9jY6PXt0dJRmlOyf6gVkbJ9qxatQoZGRkAgCeffBIPPvggiouLUVhYiMcff9zt\nDSjUGjelbSk5jt577z1MnToViYmJAC5f7srLy8Ntt93msW8vvfQSysrK0KtXLzQ0NGDOnDkoLi7G\nvn378Nhjj2HHjh1eySjdv/2OL6buBAYGyrx582T+/PmtHunp6dK7d2+PuXfeeUfi4uLkb3/7m2uZ\n2Wxuty0lmXHjxrX6fePGjTJmzBj55ptvPE7Punr5hAkTWj0XHR3tNuNp6tmJEyc8Tj1Ta3uUZJTm\nlOwPamW6ui9cO+XO076g1rgpbUvJfnfq1Cn51a9+JQ899JA4nc4OZaKiouTSpUsiInLmzJlW/fE0\nLVVJRun+7W98UuBjYmJk3759bp8zmUztZr/77jvJzMyU+++/Xw4fPnzdHUZJxmKxyNmzZ1st27Zt\nm4waNUpCQkLcZgICAqRPnz7Sp08fCQwMlNraWhEROXfunMd5tevWrZORI0fKY489Ji+//LK8/PLL\n8uijj8qIESNkzZo1Pt0eJRmlOSX7g1oZJdszePBgWbRokfzyl7+UYcOGSXNzs+u5sWPHeq0dpceR\n0vdWybEnIrJ7926x2Wzy6quvyq233truukuWLJGEhAR56aWX5Ec/+pH87ne/ExGRhoYGj8VaSUbp\nGPgbnxT40tJSOXz4sNvnysvLO/QaFRUVMnXqVBk8eHCH2+1o5rXXXpOPP/64zfI9e/bItGnTOtye\nyOWz8X/+858enz9+/Lj8+c9/lpUrV8rKlStl06ZNcvz48Q69dnduj9IxUJJTsj+olVGyPWvXrnX9\no7p169a53s/a2lp55plnvNaO0uOoq/u3kmPv0qVLsmrVKpk7d+511926dausWLGi1c02Ll261KYg\ndyXjzWNcy/x6Fo2IoKmpCf369evWjJbpbXvIPyjd79SYwUb/xycfsi5atMj1YdC1PH0/jNKcWplv\nv/0WS5YsgcPhQFJSEhYvXoygoCAArb/34mpr1qxxfRjncDiQnp6OiooKWCwWrFu3DhEREX41Blrv\nHzPqtnVldo2IuLL33nuv6/vZOzu112q14quvvmqzXMmxp3Ts/I1PzuCDgoIQFRWF1NRUDB8+HABc\nA20wGDx+hYCSnFqZadOmtfpUfs+ePa5P5T19MdLVyx944AEkJCTg4YcfRmFhId566y23sy20PAZa\n7x8z6rYVEBDQZubNzp07MXnyZADuZ960NzPo0UcfRUNDQ5vnlRx7SsfO73TDZZ/rOnbsmOTm5orN\nZpP4+Hh5++233c4m8UZOrYxasy20PAZa7x8z6ralZOaNWjOdlI6dv/H5TbdrampkxYoVMmzYMNmw\nYUO35rozo9ZsC7W2p6sZrfePGXXa6uzMG7VmOl1N6dj5A5/e0amiogJ5eXnYtm0bZs6c2eEPX5Tk\nujvz8MMPY+fOnbDZbK5l06ZNw1/+8hcsWbLEbWbFihWu//287bbb0NTUhIEDB6Kurs7jPSbV2p6u\nZLTeP2bUa6tv3754/fXXsWfPHqSnp+P06dPtrv/66697/ODW0/fIKDn2lG6P3/HFX5Xnn39eJkyY\nIHPnzpX333+/1Zmrt3NqZdSi9THQcv+YUb+tq7W0tMipU6c6nesOWj7GvcknH7IGBARgxIgR6NWr\nV5vnDAaDx++PUJJTK6PWzAQtj4HW+8eMum1peTaR0rHzNz65RHPo0CEAcPumtHcDASU5tTL/9V//\n1e6n8t7KaHkMtN4/ZtRtS61jQq3jyC+p8b8JnixZsqRDy7yR6+6MmrMgOts3tTNqtsWMdt8jrc8m\n6uz2+COfFnh3U5iioqK6JadWRkSdWRBaHwMt948Z9dvS6mwipdvjL3xyiWbVqlV4++23cfDgwVbf\nW93U1IQf/ehHXs2plbmiu2cmaH0MtNw/ZtRvC9DmbKKubI9f8cVflfHjx4vdbpcHH3xQDh8+LHa7\nXex2uzQ0NHg9p1ZGrZkJWh4DrfePGXXb0vJsIqVj5298MovG0z8f7o6cWhm1ZiZoeQzUbIsZ7b9H\nWp5NpHTs/I1PLtEcO3YM//Ef/+FxWtOvf/1rr+XUyqg1M0HLY6D1/jGjbltank2kdOz8jU8K/KVL\nl9DU1KRKTq2M2WwGACxduhSvvPJKq+fcLVOa0fIYqNkWM9p/j9Q6JtQ6jvySaheDrqL0llhKcmpl\n2st6c2aC1sdAy/1jRv22PGV9PZtIT7fla49Pv4tGT9SemUCkdVqfTXRD8MVfFaWfVCvJqZVRa2aC\nlsdAzbaY0f57pOXZRHqbLeOJX9+yT0vUnAVB5A+0PpvoRsBLNF6i5iwIIn+g9dlENwIWeC9RcxYE\nkT/Q+myiGwEv0XgJ/9eSqDVeovG9AF93gIiIugfP4L3k+PHjGDRoULdniPyFWscEjyPPWOCJiHSK\nl2iIiHSKBZ6ISKdY4ImIdIoFnohIp1jgiYh06v8DDd4GiRZ5kRoAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f6b8a53c250>"
       ]
      }
     ],
     "prompt_number": 391
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "doubles['TtM47_TtM51']/float(sum(doubles.values()))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 392,
       "text": [
        "0.9769144144144144"
       ]
      }
     ],
     "prompt_number": 392
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## In 98% of the cases where there 'meosis' split 2:[the rest] the 2 where 47 and 51\n",
      "\n",
      "So somthing is indeed up with those samples. Let's dif a little deeper and see if these two are different in any other ways:\n",
      "\n",
      "### The ancestor don't look any less like heterozygotes:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def is_47_51(rec):\n",
      "    if rec.num_hom_alt == 2:\n",
      "        return \"_\".join([s.sample for s in rec.get_hom_alts()]) == 'TtM47_TtM51'\n",
      "    if rec.num_hom_ref == 2:\n",
      "        return \"_\".join([s.sample for s in rec.get_hom_refs()]) == 'TtM47_TtM51'\n",
      "    return False\n",
      "\n",
      "weird_ones = np.array([is_47_51(r) for r in hets])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 393
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def get_het_qual(rec):\n",
      "    return rec.samples[0].data.PL[1]\n",
      "Q = np.array([get_het_qual(r) for r in hets])\n",
      "Q_w = Q[weird_ones]\n",
      "Q_o = Q[np.invert(weird_ones)]\n",
      "print \"\\t\\tm47/51\\t\\tNot\"\n",
      "print \"Mean\\t\\t\",mean(Q_w), mean(Q_o)\n",
      "print \"std dev\\t\\t\", std(Q_w), std(Q_o)\n",
      "print \"prop. >0\\t\", round(mean(Q_w >0),5), round(mean(Q_o > 0))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\t\tm47/51\t\tNot\n",
        "Mean\t\t0.463400576369 0.286481364498\n",
        "std dev\t\t13.7760060787 15.0504013353\n",
        "prop. >0\t0.00231 0.0\n"
       ]
      }
     ],
     "prompt_number": 394
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### They have pretty similar coverage levels"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "coverage = np.array([r.INFO[\"DP\"] for r in hets])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 395
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "scatter(weird_ones, coverage)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 396,
       "text": [
        "<matplotlib.collections.PathCollection at 0x7f6b8af0f7d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4VPW97/H3xAxChCCKJDhDGySDYSAkqXTAttpQCELU\nlIplG6oExe6e2B6gerx1HzX0QnC3u89Gbaq7G1s2rSK9QNinko2ljFcIGIhaghJokGRIIgiRQIJJ\nyO/8seIUBMkkjLMyzuf1PHmYrLl9Vh7WfGet381hjDGIiEhMirM7gIiI2EdFQEQkhqkIiIjEMBUB\nEZEYpiIgIhLDVARERGJYSEXgxIkTTJw4kczMTLxeLw8++CAARUVFuN1usrKyyMrKYv369cHnFBcX\n4/F4SEtLY8OGDcHtFRUVpKen4/F4WLhwYZh3R0REesIR6jiBlpYWEhIS6Ojo4Ctf+Qo/+9nP2Lhx\nI4MGDeLuu+8+7bFVVVXMmTOHbdu2EQgEmDp1KtXV1TgcDnw+H0888QQ+n4/c3FwWLFjA9OnTP5Wd\nExGRcwv5clBCQgIAbW1tnDx5kiFDhgBwthpSWlpKfn4+TqeTlJQUUlNTKS8vp76+nubmZnw+HwBz\n585l7dq14dgPERHphZCLQGdnJ5mZmSQlJTF58mTGjh0LwOOPP05GRgbz58+nqakJgAMHDuB2u4PP\ndbvdBAKBM7a7XC4CgUC49kVERHoo5CIQFxdHZWUldXV1vPTSS/j9fgoLC6mpqaGyspLhw4dzzz33\nfJpZRUQkzOJ7+oTBgwdz/fXX8/rrr5OdnR3cfuedd3LjjTcC1jf82tra4H11dXW43W5cLhd1dXWn\nbXe5XGe8R2pqKnv37u1pNBGRmDZq1Cj27NnTo+eEdCZw6NCh4KWe1tZWXnjhBbKysmhoaAg+Zs2a\nNaSnpwOQl5fHqlWraGtro6amhurqanw+H8nJySQmJlJeXo4xhpUrVzJz5swz3m/v3r0YY6L255FH\nHrE9QyxmV377f5Tf3p/efHkO6Uygvr6egoICOjs76ezs5LbbbmPKlCnMnTuXyspKHA4HI0eO5Kmn\nngLA6/Uye/ZsvF4v8fHxlJSU4HA4ACgpKWHevHm0traSm5urnkEiIjYKqQikp6ezffv2M7b/13/9\n1yc+5wc/+AE/+MEPzth+1VVX8dZbb/UgooiIfFo0YvhTcGpbSbSJ5uyg/HZT/ugT8mCxSHI4HPTB\nWCIifVpvPjt1JiAiEsNUBEREYpiKgIhIDFMREBGJYSoCIiIxTEVARCSGqQiIiMQwFQERkRimIiAi\nEsNUBEREYpiKgIhIDFMREBGJYSoCIiIxTEVARCSGqQiIiMQwFQERkRimIiAiEsNCKgInTpxg4sSJ\nZGZm4vV6efDBBwE4fPgwOTk5jB49mmnTptHU1BR8TnFxMR6Ph7S0NDZs2BDcXlFRQXp6Oh6Ph4UL\nF4Z5d0REpCdCKgL9+/dn06ZNVFZW8uabb7Jp0yZeeeUVli5dSk5ODrt372bKlCksXboUgKqqKp57\n7jmqqqooKyvjrrvuCi55VlhYyPLly6murqa6upqysrJPb+8i7OjRo7z++uvU1tbaHUVEeqC1tZXt\n27ezZ8+emFvaNuTLQQkJCQC0tbVx8uRJhgwZwrp16ygoKACgoKCAtWvXAlBaWkp+fj5Op5OUlBRS\nU1MpLy+nvr6e5uZmfD4fAHPnzg0+J9pt2bKFESNGM2XKtxk9OpNHHvmJ3ZFEJAR79+5l1Kh0srML\nGD/+Gr71rTvp7Oy0O1bEhFwEOjs7yczMJCkpicmTJzN27FgaGxtJSkoCICkpicbGRgAOHDiA2+0O\nPtftdhMIBM7Y7nK5CAQC4doX2xhjyMv7J44e/Q+OHt3BiRNV/Oxn/8GWLVvsjiYi3fjWt75DY2Mh\nzc1v0dq6h3Xr/sYzzzxjd6yICbkIxMXFUVlZSV1dHS+99BKbNm067X6Hw4HD4Qh7wGjQ2trK4cON\nQF7XliQcjq+ya9cuO2OJSAjeeWcXnZ2zu367iOPHr2fnztg5duN7+oTBgwdz/fXXU1FRQVJSEg0N\nDSQnJ1NfX8+wYcMA6xv+qdfF6+rqcLvduFwu6urqTtvucrnO+j5FRUXB29nZ2WRnZ/c0asQMGDCA\nSy5J4uDBdViFoBFjXsTr/a7d0USkG1deOYZt21bT2XkPcJyLLvozY8f+b7tjhcTv9+P3+8/vRUwI\nDh48aI4cOWKMMaalpcVcc8015i9/+Yu59957zdKlS40xxhQXF5v777/fGGPMzp07TUZGhvnwww/N\n3//+d3PFFVeYzs5OY4wxPp/PbNmyxXR2dpoZM2aY9evXn/F+IcbqUzZv3mwSE5NMYmKm6d//EvPw\nwz+2O5KIhGDPnj1m+PBRJjFxnBkwIMnk599hTp48aXesXunNZ6ej64nn9NZbb1FQUEBnZyednZ3c\ndttt3HvvvRw+fJjZs2ezf/9+UlJSWL16NRdffDEAS5Ys4emnnyY+Pp5ly5Zx3XXXAVYX0Xnz5tHa\n2kpubi6PPfbYGe/ncDiisoX+6NGj7N69m+Tk5NPaPkSkb2ttbeXtt99m0KBBjBo1KmovbffmszOk\nIhBp0VoERETs1JvPTo0YFhGJYSoCIiIxTEVARCSGqQiIiMQwFQERkRimIiAiEsNUBEREYpiKgIhI\nDFMREBGJYSoCIiIxTEVARCSGqQiIiMQwFQERkRimIiAiEsNUBEREYpiKgIhIDFMREBGJYSoCIiIx\nTEVARCSGhVQEamtrmTx5MmPHjmXcuHHBxeGLiopwu91kZWWRlZXF+vXrg88pLi7G4/GQlpbGhg0b\ngtsrKipIT0/H4/GwcOHCMO+OiIj0REgLzTc0NNDQ0EBmZibHjh3jqquuYu3ataxevZpBgwZx9913\nn/b4qqoq5syZw7Zt2wgEAkydOpXq6mocDgc+n48nnngCn89Hbm4uCxYsYPr06aeH0kLzIiI99qkt\nNJ+cnExmZiYAAwcOZMyYMQQCAYCzvmFpaSn5+fk4nU5SUlJITU2lvLyc+vp6mpub8fl8AMydO5e1\na9f2KLCIiIRPj9sE9u3bx44dO5g0aRIAjz/+OBkZGcyfP5+mpiYADhw4gNvtDj7H7XYTCATO2O5y\nuYLFREREIi++Jw8+duwYN998M8uWLWPgwIEUFhby8MMPA/DQQw9xzz33sHz58rAEKyoqCt7Ozs4m\nOzs7LK8rIvJZ4ff78fv95/UaIReB9vZ2Zs2axa233srMmTMBGDZsWPD+O++8kxtvvBGwvuHX1tYG\n76urq8PtduNyuairqzttu8vlOuv7nVoERETkTB//grx48eIev0ZIl4OMMcyfPx+v18uiRYuC2+vr\n64O316xZQ3p6OgB5eXmsWrWKtrY2ampqqK6uxufzkZycTGJiIuXl5RhjWLlyZbCgiIhI5IV0JvDq\nq6/y29/+lvHjx5OVlQXAkiVLePbZZ6msrMThcDBy5EieeuopALxeL7Nnz8br9RIfH09JSQkOhwOA\nkpIS5s2bR2trK7m5uWf0DBIRkcgJqYtopKmLqIhIz31qXURFROSzSUVARCSGqQiIiMSwHo0TkHN7\n99132bFjB8OHD8fn8wUbw0VE+ioVgTB5/vnn+eY3C4iPn8TJkzv55jen8/TTv1AhEJE+Tb2DwsAY\nw+DBw2huLgW+BBzjoosm8N///UsmT55sdzwRiRG9+ezUmUAYtLa2cvz4B8BbwJPAcIwZz7vvvmtz\nMhGRc1PDcBgkJCRw0UWXYhWAycAHtLaWkZaWZnMyEZFzUxEIg46ODlpajgB/AW4HnqRfv6tOmz9J\nRKQvUhEIg87OTsAACcFtTucQ2tvbbcskIhIKFYEw6NevHzfcMIsBA/KBl4mL+zecznJycnLsjiYi\nck7qHRQmJ06c4L77Hmbjxldwu4fz+OPFjB492u5YIhJDevPZqSIgIvIZoQnkbPbBBx+wdetW9u/f\nb3cUEZGQqAiEyebNm/nc564kJ6eQK6/8Ag899CO7I4mIdEuXg8LAGMOwYZ/n0KFfADcC75GQ8EU2\nbnyOSZMm2R1PRGKELgfZpLW1lSNHGoEburYMw+H4Krt27bIzlohIt1QEwmDAgAFceulwoLRrSyPG\n+PF6vXbGEhHplopAGDgcDtate47BgwtJTMygf/8x3HdfIRMnTrQ7mojIOYVUBGpra5k8eTJjx45l\n3LhxPPbYYwAcPnyYnJwcRo8ezbRp02hqago+p7i4GI/HQ1paGhs2bAhur6ioID09HY/Hw8KFC8O8\nO/aZOHEitbW72bTp1+zZ8xaPPPKg3ZFERLoVUsNwQ0MDDQ0NZGZmcuzYMa666irWrl3Lr3/9a4YO\nHcp9993Ho48+ypEjR1i6dClVVVXMmTOHbdu2EQgEmDp1KtXV1TgcDnw+H0888QQ+n4/c3FwWLFjA\n9OnTTw8VZQ3DIiJ9wafWMJycnExmZiYAAwcOZMyYMQQCAdatW0dBQQEABQUFrF27FoDS0lLy8/Nx\nOp2kpKSQmppKeXk59fX1NDc34/P5AJg7d27wOSIiEnk9bhPYt28fO3bsYOLEiTQ2NpKUlARAUlIS\njY2NABw4cAC32x18jtvtJhAInLHd5XIRCATOdx9ERKSXerSozLFjx5g1axbLli1j0KBBp93ncDjC\nupRiUVFR8HZ2djbZ2dlhe20Rkc8Cv9+P3+8/r9cIuQi0t7cza9YsbrvtNmbOnAlY3/4bGhpITk6m\nvr6eYcOGAdY3/FPn0q+rq8PtduNyuairqzttu8vlOuv7nVoERETkTB//grx48eIev0ZIl4OMMcyf\nPx+v18uiRYuC2/Py8lixYgUAK1asCBaHvLw8Vq1aRVtbGzU1NVRXV+Pz+UhOTiYxMZHy8nKMMaxc\nuTL4HBERibyQege98sorXHvttYwfPz54yae4uBifz8fs2bPZv38/KSkprF69mosvvhiAJUuW8PTT\nTxMfH8+yZcu47rrrAKuL6Lx582htbSU3NzfY3fS0UOodJCLSY5pKWkQkhmnuIBER6REVgTA6ePAg\nmzZt4u2337Y7iohISFQEwsTv95OSMobc3AfIyLiWRYvutzuSiEi31CYQBsYYBg68jJaWZ4Ec4Ahx\ncePYtOlZrr32WrvjiUiMUJuATY4fP05LSxMwtWvLEDo7J2lKDBHp81QEwuCCCy4ALgRWdW2pBV5m\n8ODB9oUSEQmBLgeFSWpqOnv31gMXA+8RHx/Hzp1bGT16tN3RRCRG6HKQjfz+9Ywb5wFqGDRoAGvW\n/FYFQET6PBWBMNm3bx81NdUkJKTS1tbOK69sszuSiEi3dDkoDIwxDB48jObmk8AXgF3Ex7fz4otr\n+dKXvmR3PBGJEbocZJOWlhaam48BfwL+AlTR0dGfdevW2ZxMROTcVATCpgP4atftwcAEnE6njXlE\nRLqnIhAGF110EcOGuYEVXVv24XS+wk033WRnLBGRbqkIhMlf/rKOoUMfYcAAN/36jednPysiKyvL\n7lgiIuekhuEwam9vp66ujqFDh56x/KaIyKdN6wmIiMQw9Q4SEZEeUREQEYlhKgIiIjEspCJwxx13\nkJSURHp6enBbUVERbrebrKwssrKyWL9+ffC+4uJiPB4PaWlpbNiwIbi9oqKC9PR0PB4PCxcuDONu\niIhIb4RUBG6//XbKyspO2+ZwOLj77rvZsWMHO3bsYMaMGQBUVVXx3HPPUVVVRVlZGXfddVewoaKw\nsJDly5dTXV1NdXX1Ga8ZzQ4fPswNN/wTQ4a4uPLKCbz22mt2RxIR6VZIReCaa65hyJAhZ2w/Wyt0\naWkp+fn5OJ1OUlJSSE1Npby8nPr6epqbm/H5fADMnTv3M7XoyvTps3j++USamjaze/cDTJlyI/v3\n77c7lojIOZ1Xm8Djjz9ORkYG8+fPp6mpCYADBw7gdruDj3G73QQCgTO2u1wuAoHA+bx9n3HixAm2\nbXsFY34JfA64mRMnvsTvf/97u6OJiJxTfG+fWFhYyMMPPwzAQw89xD333MPy5cvDFqyoqCh4Ozs7\nm+zs7LC9drh1dHR03TqAVQQ6gQC1tbX2hRKRzzy/34/f7z+v1+h1ERg2bFjw9p133smNN94IWN/w\nT/3wq6urw+1243K5qKurO227y+X6xNc/tQj0dQkJCcTFOenszAbmAtuIi6vr04VLRKLfx78gL168\nuMev0evLQfX19cHba9asCfYcysvLY9WqVbS1tVFTU0N1dTU+n4/k5GQSExMpLy/HGMPKlSuZOXNm\nb9++T4mLi+OJJ/6dfv2OAv+N0/kOkyZlBAujiEhfFdKZQH5+Pi+++CKHDh1ixIgRLF68GL/fT2Vl\nJQ6Hg5EjR/LUU08B4PV6mT17Nl6vl/j4eEpKSnA4HACUlJQwb948Wltbyc3NZfr06Z/enkXY5MnX\nkpg4kKamv2NMG7fddn/XAvQiIn2X5g4Kk5SUsezfvwBjvgPsJiHhq2ze/D+MHz/e7mgiEiM0d5BN\nWlpaqK2txpgPgK8DP8eYL7F9+3a7o4mInJOKQBgMGDCAuLgE4P8Bc4BEWltf4NJLL7U5mYjIufW6\nd5D8Q0dHBx0dx4Bq4HbAicMxlGPHjtmcTETk3HQmEDYJwI1YheA3GNPA3/72N5sziYicm84EwqC5\nuRloAX4JXAB8A5iiNgER6fN0JhAG8fHxgAP4aJBcJ1DbtV1EpO9SEQiDhIQErJOqrwD/F5gG1JCR\nkWFrLhHpnjGGJ574JRMn5jBjxjd544037I4UUfqqGgbHjx8H2oEPgDKgEXDG3H8mkWj0gx88zL/+\n61N0dh4FDBs3bmDnztfxeDx2R4sInQmEQWdnJ9af8vWunxpAvYNEosHPf/4knZ03AceBvbS3J/Zq\nDp5opSIQBv3798dqB/jom0M8cCUDBw60L5SIhKStrRV4EKtThxuYT2Xlm/aGiiBdDgqDDz74ALgQ\nmAL0x/qz/pVDhzRlhEjfdwGwBfg81pe5V3nvvUZ7I0WQzgTCIDExEat30AngfwFjgFPXGRCRvqsD\nuAu4Gbgaa+6vBHsjRZDOBMLgvffewyoA/wMkYs0ftJ2amkpbc4lI9xITB3L0aAfWgM9hwFvk5y+0\nOVXk6EwgrByn3O7UVNIiUeC11/7KhRe2Y8395eeLX5zAj3/8Y7tjRYyKQBh8/vOfx2oTGId1KcgL\nbOZrX/uarblEpHter5fp06fSv387gwf348EHvx9TX+BUBMLAWjYzDsgAfgP8MxDHpk2b7IwlIiH4\nyleupbR0AydOTOSDD0Zx0023smHDBrtjRYwWlQmDN998k4yMLwDHsHoHAVxPXFwZJ0+etDGZiHTH\n4bgI60w+H3gX2MYVVwxi79499gbrhd58dqphOAxaWlq6bjVjFQEDvN81iExE+rZ44Bngo+Vuv877\n779sY57IUhEIgwEDBgD9gKlYXc3KgSpbM4lIqDqBK0/5fTyXXbbLrjARF1KbwB133EFSUhLp6enB\nbYcPHyYnJ4fRo0czbdo0mpqagvcVFxfj8XhIS0s77dpaRUUF6enpeDweFi787HTBssYDfAgsArYB\nLuCamGpcEolWTmc/4AHgEFAB/AeTJk2yN1QEhVQEbr/9dsrKyk7btnTpUnJycti9ezdTpkxh6dKl\nAFRVVfHcc89RVVVFWVkZd911V/AaVWFhIcuXL6e6uprq6uozXjNaXXbZZViXgTqAPOCrwBYcDsc5\nnyci9vv6168HtgOpwI1ccMGH3HnnfJtTRU5IReCaa65hyJAhp21bt24dBQUFABQUFLB27VoASktL\nyc/Px+l0kpKSQmpqKuXl5dTX19Pc3IzP5wNg7ty5wedEuyNHjgCtWGcCNwEzgBa1CYhEgT17dgIN\nQBtwlJMnj9He3m5zqsjpdRfRxsZGkpKSAEhKSqKx0Zpr48CBA7jd7uDj3G43gUDgjO0ul4tAINDb\nt+9TrH2/EKtd4HpgBOBUERCJApWV72AN9LwYax6hAdx00yx7Q0VQWBqGHQ5H2C99FBUVBW9nZ2eT\nnZ0d1tcPpy1btmD1CHoea+6R48BYrN5CItK3GawP/xuBfcDrNDcftjVRqPx+P36//7xeo9dFICkp\niYaGBpKTk6mvr2fYsGGA9Q2/trY2+Li6ujrcbjcul6trUNU/trtcrk98/VOLQF9n7UcH1hnA77Cm\nlP4CVp9jEenb4oGVwA1YBeHrWFNI9H0f/4Lcm3UQen05KC8vjxUrVgCwYsUKZs6cGdy+atUq2tra\nqKmpobq6Gp/PR3JyMomJiZSXl2OMYeXKlcHnRLu3334b63TyCmAeMAn4s52RRCRkJ4EU4GXgTayR\n/7EjpDOB/Px8XnzxRQ4dOsSIESP44Q9/yAMPPMDs2bNZvnw5KSkprF69GrDm4Zg9ezZer5f4+HhK\nSkqCl4pKSkqYN28era2t5ObmMn369HO9bdSwppI2WD2EpgCVwGGshiYR6ds6gS8DLYAT6wtd9MxY\ncL40bUQYLFy4kMceewp4Desy0AmsSeRqomo/RGKRwzEIuAirXeBo178fROWx25vPTk0gFwYHDx7E\nahPI6trSH6sYiEjfFwfcCvwn8Hsgyd44EaZpI8LAmiQuAfg5cDfwNyB2ZiEUiW7twG+BHcABrK7e\nsUNnAmFgLUV3HFiMNV5gArF0TVEkuvUDfgJsxGoYvsjeOBGmIhAGx48fx/rQn4C1vOTniLUeBiLR\nqx3I6brtBGZwwQVOG/NElopAGDQ3NwMDgUuBTcB9WHORiEjf5wRKsL7IHQJWcvKkpo2QHoiPj8da\nUGY3VhfRR4DxtmYSkVC1YY32Hw6MxBr0GTvUMBwG1hrDA7D6F/87VlfR39gZSURC1g/wYU3zcjHw\nor1xIkxnAmHgdDqxxgZ8C2u4+UeLzotI33cUeBbYi3X8HrA3ToSpCITBO++8g3Vd8U9Y00i3Y10a\nEpG+LxH4N+B1oAZIszdOhGnEcBh4vV527aoG3gM+WnfhK8CrUbUfIrHI4RiA1S20DWu08CXA36Py\n2NWIYZv8Y+6gC0/ZOsCmNCLSM/FYl2+3YC04f9DeOBGmIhAG1ipE/bGWltwELAFetTWTiPTEf2LN\n9zUD+J7NWSJLvYPCYP/+/VgzEVYC12GdBXTYmklEQhWH1Ric2vX7fhuzRJ7OBMLAGiwWB9wOvA+s\nIdbmHxGJXsewFpL5EdYxvMbeOBGmhuEwSExMpLm5Bath6aO6OhMojar9EIlFDkd/rIVlLsDq2RcP\ntEXlsauGYZscO3YM60/5Z6xri+tQF1GRaHEh1ij/E0AT1gqBsUNFIAysyuvAGiz2V+BeoO6czxGR\nvuIkMLfr9iAgH+t4jg0qAmH1ClYXs7eItYUpRKJXHNZAT4BW4I/E0lTwahMIA2sN5Tjg+8BLWBNR\ntQIvRNV+iMQih2Mg1iWhYVgdO+KB+qg8dm1pE0hJSWH8+PFkZWXh8/kAOHz4MDk5OYwePZpp06bR\n1NQUfHxxcTEej4e0tDQ2bPgsrb41AGviuIeBsWicgEi0MEAK8F1gAdYCUbHjvM8ERo4cSUVFBZdc\ncklw23333cfQoUO57777ePTRRzly5AhLly6lqqqKOXPmsG3bNgKBAFOnTmX37t3ExZ1ei6LzTOAC\n4AjWNUWArwGbomo/RGKRw5GINVrY27XlHuDnUXns2tY76ONvum7dOgoKCgAoKChg7dq1AJSWlpKf\nn4/T6SQlJYXU1FS2bt0ajgh9xKl/h9hpWBKJfi2n3G62LYUdzrsIOBwOpk6dyoQJE/jVr34FQGNj\nI0lJVsNoUlISjY2NABw4cAC32x18rtvtJhAInG+EPqI/cD3WQJMHsL5ZiEjf1441XcQlQDLWovOx\n47ynjXj11VcZPnw4Bw8eJCcnh7S006dhdTgcXZdLzu6T7isqKgrezs7OJjs7+3yjfsrasZaUnI81\naOykvXFEJERO4GpgGRDAGj3camuiUPn9fvx+/3m9xnkXgeHDhwNw2WWX8Y1vfIOtW7eSlJREQ0MD\nycnJ1NfXM2zYMABcLhe1tbXB59bV1eFyuc76uqcWgejQCezCOiMYDHwJrTMsEg3aga2AB+viSApw\n2M5AIfv4F+TFixf3+DXO63JQS0tL17w5cPz4cTZs2EB6ejp5eXmsWLECgBUrVjBz5kwA8vLyWLVq\nFW1tbdTU1FBdXR3sURT9DHA5VjezC4GzFzcR6Wv6AblY3/7fwVppLHac15lAY2Mj3/jGNwDo6Ojg\nW9/6FtOmTWPChAnMnj2b5cuXk5KSwurVqwFr8ZXZs2fj9XqJj4+npKTknJeKoksC8G3gIaxvFX+x\nN46IhOgkMBW4D2uVsVuBn9uaKJI0WCwMrEI2AOs/UDPWJaFjROskVCKxxOFIwLqE+31gH9ao/w+i\n8tjtzWen1hMImw+xvlHMwJo2YgjWwtUi0rfFA7/HWhIWrC9ysdNDSEUgbJxYs4j6sPoce8/9cBHp\nIzo5fa6v2GrP0wRyYdMBfLHrdgKQZWMWEQmdA7gTeAMoBZ60N06EqQiETX/gia7bVYDfvigi0kNX\nAHOApcAkm7NElhqGw+Afs4h6gHqs08sRwK6o2g+RWORwXIjVhvdDYA/wC6AlKo9drSxmq496By3D\n6ir6rr1xRCREA4DpWFO9HAM+b2+cCFMRCJvjWCOG7wV+hTVwTET6vjasjhzvY/UM+hqxNAGkikDY\nJAD/DmwANgKH7I0jIiGKA5YDtwEZwG/QymI2i7Y2gX79+tHe3g4MxDq1PIo1fURtVO2HSCyy1hN4\nGasAAHwP+EVUHrtqE7BJ//79gYuwxgj8HLgf69RSRKLDqZd/YutjUWcCYeD1etm1axfwOayGJQcw\nFHgnqvZDJBY5HBcBI4EfYU0bUQQcjcpjV9NG2KS1tZV/nAlMBRqBn9qaSURCFQfMBJ7GWh52OrDa\n1kSRpCIQBidOnMCahtYPfAC8jdU7aLeNqUQkNCexevS5scb4vGNvnAiLrYtfn5IPP/wQa07yP2D1\nDvobVpdREen7DNYsokXAHcRS91DQmUBYXHjhhVg9guqAW7BWJroKa6k6EenbnMAfgfSu3/diDfqM\nDToTCIMZ51enAAAKhklEQVR/TBtR2LXlT2juIJFocuqa4G22pbCDegeFgVUEnMAjwH6sMQK/BfZF\n1X6IxCKHoz/WVNJLgRpgCXA8Ko9djROwlQNrYYoMoBZrTQER6fucWNNG/AnYiTUHWOywpQiUlZWR\nlpaGx+Ph0UcftSNCWFmDxU5iNQrfBfyaWJuESiR6GeAprC9xvwP+mVhqHI54ETh58iTf+973KCsr\no6qqimeffbZroFX0uuKKK7puDe761wFcYlMaEemZC4C/dt3u6LodfZeCeiviRWDr1q2kpqaSkpKC\n0+nklltuobS0NNIxwuo73/kO1pxBs4HXsRaXednWTCISqjZgATANq4fQG/bGibCIdxENBAKMGDEi\n+Lvb7aa8vDzSMcKqvr4ea7qIF4AXsQacxFYPA5HoFY91Fv9y1+0L7I0TYRE/E7B60ny2OJ1OrD9l\nMtZ85O1YYwVEpO87DjRhHcMfAif4/Odjp00v4mcCLpeL2tra4O+1tbW43e4zHldUVBS8nZ2dTXZ2\ndgTS9U5RURE/+tHPsS4J/RuwHVhDQkKCvcFEpFs5OV/nhRfWYq0JcgHQQUlJic2pQuP3+/H7/ef1\nGhEfJ9DR0cGVV17Jxo0bufzyy/H5fDz77LOMGTPmH6GibJwAwOWXX059/RGsRuFO4CTGtNucSkS6\n8/7773PDDf/E1q0vExcXxw9/+CMefPD/2B2rV3rz2WnLYLH169ezaNEiTp48yfz583nwwQdPDxWF\nRQDg9ddfZ/HixeTm5lJYWNj9E0Skz2hpaaFfv37Ex0fvbDpRUwS6E61FQETEThoxLCIiPaIiICIS\nw1QERERimIqAiEgMUxEQEYlhKgIiIjFMRUBEJIapCIiIxDAVARGRGKYiICISw1QERERimIqAiEgM\nUxEQEYlhKgIiIjFMRUBEJIapCIiIxDAVARGRGKYiICISw1QERERiWK+LQFFREW63m6ysLLKysli/\nfn3wvuLiYjweD2lpaWzYsCG4vaKigvT0dDweDwsXLjy/5CIict56XQQcDgd33303O3bsYMeOHcyY\nMQOAqqoqnnvuOaqqqigrK+Ouu+4KLnxcWFjI8uXLqa6uprq6mrKysvDsRR/j9/vtjtBr0ZwdlN9u\nyh99zuty0NlWtS8tLSU/Px+n00lKSgqpqamUl5dTX19Pc3MzPp8PgLlz57J27drzefs+K5r/I0Vz\ndlB+uyl/9DmvIvD444+TkZHB/PnzaWpqAuDAgQO43e7gY9xuN4FA4IztLpeLQCBwPm8vIiLn6ZxF\nICcnh/T09DN+1q1bR2FhITU1NVRWVjJ8+HDuueeeSGUWEZFwMWFQU1Njxo0bZ4wxpri42BQXFwfv\nu+6668yWLVtMfX29SUtLC25/5plnzHe+852zvt6oUaMMoB/96Ec/+unBz6hRo3r8+R1PL9XX1zN8\n+HAA1qxZQ3p6OgB5eXnMmTOHu+++m0AgQHV1NT6fD4fDQWJiIuXl5fh8PlauXMmCBQvO+tp79uzp\nbSwREemBXheB+++/n8rKShwOByNHjuSpp54CwOv1Mnv2bLxeL/Hx8ZSUlOBwOAAoKSlh3rx5tLa2\nkpuby/Tp08OzFyIi0isOY87SxUdERGJCnxgxfPjwYXJychg9ejTTpk0L9jQ6VW1tLZMnT2bs2LGM\nGzeOxx57zIak/1BWVkZaWhoej4dHH330rI9ZsGABHo+HjIwMduzYEeGE59Zd/t/97ndkZGQwfvx4\nvvzlL/Pmm2/akPKThfL3B9i2bRvx8fH86U9/imC67oWS3+/3k5WVxbhx48jOzo5swG50l//QoUNM\nnz6dzMxMxo0bx29+85vIh/wEd9xxB0lJScFL2GfTl4/d7vL3+NjtcSvCp+Dee+81jz76qDHGmKVL\nl5r777//jMfU19ebHTt2GGOMaW5uNqNHjzZVVVURzfmRjo4OM2rUKFNTU2Pa2tpMRkbGGVn+/Oc/\nmxkzZhhjjNmyZYuZOHGiHVHPKpT8r732mmlqajLGGLN+/fqoy//R4yZPnmyuv/5684c//MGGpGcX\nSv4jR44Yr9dramtrjTHGHDx40I6oZxVK/kceecQ88MADxhgr+yWXXGLa29vtiHuGl156yWzfvj3Y\nmeXj+vKxa0z3+Xt67PaJM4F169ZRUFAAQEFBwVkHkSUnJ5OZmQnAwIEDGTNmDAcOHIhozo9s3bqV\n1NRUUlJScDqd3HLLLZSWlp72mFP3aeLEiTQ1NdHY2GhH3DOEkv/qq69m8ODBgJW/rq7OjqhnFUp+\nsMax3HzzzVx22WU2pPxkoeR/5plnmDVrVnBszdChQ+2Ielah5B8+fDhHjx4F4OjRo1x66aXEx/e6\nCTKsrrnmGoYMGfKJ9/flYxe6z9/TY7dPFIHGxkaSkpIASEpK6vYPvm/fPnbs2MHEiRMjEe8MgUCA\nESNGBH//aEBcd4/pKx+koeQ/1fLly8nNzY1EtJCE+vcvLS2lsLAQINg5oS8IJX91dTWHDx9m8uTJ\nTJgwgZUrV0Y65icKJf+3v/1tdu7cyeWXX05GRgbLli2LdMxe68vHbk+FcuxGrDTn5OTQ0NBwxvaf\n/OQnp/3ucDjOecAeO3aMm2++mWXLljFw4MCw5wxFqB8o5mNt7n3lg6gnOTZt2sTTTz/Nq6+++ikm\n6plQ8i9atIilS5ficDgwxpx1ihO7hJK/vb2d7du3s3HjRlpaWrj66quZNGkSHo8nAgnPLZT8S5Ys\nITMzE7/fz969e8nJyeGNN95g0KBBEUh4/vrqsdsToR67ESsCL7zwwifel5SURENDA8nJydTX1zNs\n2LCzPq69vZ1Zs2Zx6623MnPmzE8rardcLhe1tbXB32tra0+bEuNsj6mrq8PlckUs47mEkh/gzTff\n5Nvf/jZlZWXnPP2MtFDyV1RUcMsttwBWI+X69etxOp3k5eVFNOvZhJJ/xIgRDB06lAEDBjBgwACu\nvfZa3njjjT5RBELJ/9prr/Ev//IvAIwaNYqRI0fyzjvvMGHChIhm7Y2+fOyGqkfHblhbLHrp3nvv\nNUuXLjXGWCOOz9Yw3NnZaW677TazaNGiSMc7Q3t7u7niiitMTU2N+fDDD7ttGN68eXOfalwKJf+7\n775rRo0aZTZv3mxTyk8WSv5TzZs3z/zxj3+MYMJzCyX/rl27zJQpU0xHR4c5fvy4GTdunNm5c6dN\niU8XSv7vf//7pqioyBhjTENDg3G5XOb999+3I+5ZnTrLwcf15WP3I+fK39Njt08Ugffff99MmTLF\neDwek5OTY44cOWKMMSYQCJjc3FxjjDEvv/yycTgcJiMjw2RmZprMzEyzfv162zI///zzZvTo0WbU\nqFFmyZIlxhhjnnzySfPkk08GH/Pd737XjBo1yowfP95UVFTYFfWsuss/f/58c8kllwT/1l/84hft\njHuGUP7+H+lrRcCY0PL/9Kc/NV6v14wbN84sW7bMrqhn1V3+gwcPmhtuuMGMHz/ejBs3zvzud7+z\nM+5pbrnlFjN8+HDjdDqN2+02y5cvj6pjt7v8PT12NVhMRCSG9YneQSIiYg8VARGRGKYiICISw1QE\nRERimIqAiEgMUxEQEYlhKgIiIjFMRUBEJIb9f+6cnMC4sCZJAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f6b8b916b10>"
       ]
      }
     ],
     "prompt_number": 396
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cov_47_51 = coverage[weird_ones]\n",
      "cov_rest = coverage[invert(weird_ones)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 397
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figure()\n",
      "boxplot([cov_47_51, cov_rest])\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFx1JREFUeJzt3W9sU9f9x/HPbZPfqikwWLU4mS+aq8Y0GAKJikw7jS0o\nJLB0jVJRZQ1rkxT6JIyVij6Y9mAjPFiTPaLQNlI1hTWiGn80aQnaVIttwmggLbRR0lVzJdwJmO04\n0RgNSlZaKPj3IIuT8MeBJL5Ojt8vyZVz7et7TJ2PT8753nOtRCKREADASA9kugEAgPQh5AHAYIQ8\nABiMkAcAgxHyAGAwQh4ADJYy5D///HOtW7dOpaWl8vl8+tnPfiZJunz5siorK7V8+XJVVVVpeHg4\nuU9ra6u8Xq+Ki4t14sSJ5Pbe3l6VlJTI6/Vq165daXo7AIDJUob8Qw89pJMnT6q/v19///vfdfLk\nSZ0+fVptbW2qrKzUuXPnVFFRoba2NklSKBTS0aNHFQqFFAgEtGPHDo2X4Tc3N6ujo0PhcFjhcFiB\nQCD97w4Asty0wzVf/epXJUnXrl3TjRs3tHTpUh0/flyNjY2SpMbGRnV1dUmSuru7VV9fr9zcXHk8\nHhUVFamnp0fxeFwjIyPy+/2SpIaGhuQ+AID0mTbkb968qdLSUrlcLm3YsEErV67U0NCQXC6XJMnl\ncmloaEiSNDAwINu2k/vatq1YLHbbdrfbrVgsNtfvBQBwi5zpnvDAAw+ov79fV65c0aZNm3Ty5Mkp\nj1uWJcuy0tZAAMDMTRvy4772ta/pqaeeUm9vr1wulwYHB1VQUKB4PK78/HxJYz30SCSS3Ccajcq2\nbbndbkWj0Snb3W73bccoLS3Vhx9+OJv3AwBZZ82aNerv77/jYymHay5dupSsnLl69ar+9Kc/qays\nTDU1Ners7JQkdXZ2qra2VpJUU1OjI0eO6Nq1azp//rzC4bD8fr8KCgq0ePFi9fT0KJFI6NChQ8l9\nJvvwww+VSCS4zdHte9/bk/E2cON2pxufzbm9peocp+zJx+NxNTY26ubNm7p586ZeeOEFVVRUqKys\nTHV1dero6JDH49GxY8ckST6fT3V1dfL5fMrJyVF7e3tyKKe9vV1NTU26evWqqqurtXnz5lSHBgDM\nASuRSMybpYYty9I8as6C19TUonfeacl0M4Db8NmcW6mykzNeDdbUVJ7pJgB3xGfTOfTkAWCBoycP\nAFmKkAcAgxHyAGAwQh4ADEbIA4DBCHkAMBghDwAGI+QBwGCEPAAYjJAHAIMR8gBgMEIeAAxGyAOA\nwQh5ADAYIQ8ABiPkAcBghDwAGIyQBwCDEfIAYDBCHgAMRsgDgMEIeQAwGCEPAAYj5AHAYIQ8ABiM\nkAcAg6UM+Ugkog0bNmjlypVatWqVDhw4IElqaWmRbdsqKytTWVmZ3nvvveQ+ra2t8nq9Ki4u1okT\nJ5Lbe3t7VVJSIq/Xq127dqXp7QAAJrMSiUTibg8ODg5qcHBQpaWlGh0d1eOPP66uri4dO3ZMixYt\n0u7du6c8PxQKaevWrXr//fcVi8W0ceNGhcNhWZYlv9+vN998U36/X9XV1Xr55Ze1efPmqY2xLKVo\nDgDgDlJlZ8qefEFBgUpLSyVJeXl5WrFihWKxmCTd8QW7u7tVX1+v3NxceTweFRUVqaenR/F4XCMj\nI/L7/ZKkhoYGdXV1zepNAQCmd89j8hcuXFBfX5+eeOIJSdIbb7yhNWvWaPv27RoeHpYkDQwMyLbt\n5D62bSsWi9223e12J78sAADpc08hPzo6qmeffVb79+9XXl6empubdf78efX396uwsFCvvvpqutsJ\nAJiBnOmecP36dW3ZskXPP/+8amtrJUn5+fnJx1966SU9/fTTksZ66JFIJPlYNBqVbdtyu92KRqNT\ntrvd7jser6WlJXm/vLxc5eXl9/WGAMB0wWBQwWDwnp6bcuI1kUiosbFRDz/8sPbt25fcHo/HVVhY\nKEnat2+f3n//ff32t79NTryePXs2OfH6ySefyLIsrVu3TgcOHJDf79dTTz3FxCsAzJFU2ZmyJ3/m\nzBm9++67Wr16tcrKyiRJr732mg4fPqz+/n5ZlqVHHnlEb7/9tiTJ5/Oprq5OPp9POTk5am9vl2VZ\nkqT29nY1NTXp6tWrqq6uvi3gAQBzL2VP3mn05AHg/s24hBIAsLAR8gBgMEIeAAxGyAOAwQh5ADAY\nIQ8ABiPkAcBghDwAx93jGfmYA4Q8AMcR8s4h5AHAYNOuQgkAcyEYnOjB7907sb28fOyG9GDtGgCO\na2qS3nkn060wB2vXAJhXLlzIdAuyByEPwHEeT6ZbkD0YkwfgiMlj8p2dE0HPmHx6EfIAHHFrmE+6\n0ifSiOEaADAYIQ/AcQzPOIcSSgBY4CihBIAsRcgDgMEIeQCOY4Ey5xDyABxHyDuHkAcAg3EyFABH\nsAplZhDyABzBGa+ZwXANABiMkAfgOIZnnMMZrwCwwM34jNdIJKINGzZo5cqVWrVqlQ4cOCBJunz5\nsiorK7V8+XJVVVVpeHg4uU9ra6u8Xq+Ki4t14sSJ5Pbe3l6VlJTI6/Vq165dc/G+AADTSBnyubm5\n2rdvn/7xj3/ob3/7m9566y19/PHHamtrU2Vlpc6dO6eKigq1tbVJkkKhkI4ePapQKKRAIKAdO3Yk\nv12am5vV0dGhcDiscDisQCCQ/ncHAFkuZcgXFBSotLRUkpSXl6cVK1YoFovp+PHjamxslCQ1Njaq\nq6tLktTd3a36+nrl5ubK4/GoqKhIPT09isfjGhkZkd/vlyQ1NDQk9wEApM89T7xeuHBBfX19Wrdu\nnYaGhuRyuSRJLpdLQ0NDkqSBgQHZtp3cx7ZtxWKx27a73W7FYrG5eg8AgLu4pzr50dFRbdmyRfv3\n79eiRYumPGZZlizLmrMGtUwqni0vL1c50/AAMEUwGFTwHteGmDbkr1+/ri1btuiFF15QbW2tpLHe\n++DgoAoKChSPx5Wfny9prIceiUSS+0ajUdm2LbfbrWg0OmW72+2+4/FaOEMCAFK6tQO8d/IpxLdI\nOVyTSCS0fft2+Xw+vfLKK8ntNTU16uzslCR1dnYmw7+mpkZHjhzRtWvXdP78eYXDYfn9fhUUFGjx\n4sXq6elRIpHQoUOHkvsAANInZZ386dOn9d3vflerV69ODsm0trbK7/errq5O//rXv+TxeHTs2DEt\nWbJEkvTaa6/p4MGDysnJ0f79+7Vp0yZJYyWUTU1Nunr1qqqrq5PlmFMaQ508ANy3VNnJyVAAsMBx\n+T8AyFKEPAAYjJAHAIMR8gBgMEIegOO4xqtzCHkAjnvnnUy3IHsQ8gAcd+FCpluQPbjGKwBHTL6Q\n96lTE9d45ULe6UVPHgAMxhmvABxXXs7k61zijFcA84rHk+kWZA9CHoDjmpoy3YLsQcgDgMEIeQCO\nYzzeOYQ8AMdRJ+8c6uQBOGJynXxn58TkK3Xy6UUJJQDHUUI5t1JlJz15AI7gjNfMoCcPwHFNTSxS\nNpc4GQrAvMLJUM4h5AE4juEZ5zBcAwALHMM1AJClCHkAMBghDwAGI+QBwGCEPAAYjJAHAINNG/Lb\ntm2Ty+VSSUlJcltLS4ts21ZZWZnKysr03nvvJR9rbW2V1+tVcXGxTpw4kdze29urkpISeb1e7dq1\na47fBgDgTqYN+RdffFGBQGDKNsuytHv3bvX19amvr0/f//73JUmhUEhHjx5VKBRSIBDQjh07krWb\nzc3N6ujoUDgcVjgcvu01AWQPFidzzrQhv379ei1duvS27XcqvO/u7lZ9fb1yc3Pl8XhUVFSknp4e\nxeNxjYyMyO/3S5IaGhrU1dU1B80HsBCxbo1zZjwm/8Ybb2jNmjXavn27hoeHJUkDAwOybTv5HNu2\nFYvFbtvudrsVi8Vm0WwACxkXDXHOjJYabm5u1i9+8QtJ0s9//nO9+uqr6ujomJMGtYyvPyqpvLxc\n5SxyARiBpYbnTjAYVPAex7xmFPL5+fnJ+y+99JKefvppSWM99EgkknwsGo3Ktm253W5Fo9Ep291u\n9x1fe3LIAwBud2sHeO/evXd97oxCPh6Pq7CwUJL0+9//Pll5U1NTo61bt2r37t2KxWIKh8Py+/2y\nLEuLFy9WT0+P/H6/Dh06pJdffnkmhwawQE3usQeDEz15pNe0IV9fX69Tp07p0qVLWrZsmfbu3atg\nMKj+/n5ZlqVHHnlEb7/9tiTJ5/Oprq5OPp9POTk5am9vl2VZkqT29nY1NTXp6tWrqq6u1ubNm9P7\nzgDMW6wn7xyWGgbguNdfl155JdOtMAdLDQOYV/5XkAcHEPIAYLAZTbwCwP2aXEI5uRiEEsr0Ykwe\ngOOamjjrdS4xJg9gXuGMV+cQ8gAc99lnmW5B9mBMHoAjJo/Jv/8+yxo4hTF5AI5bsoQyyrmUKjvp\nyQNwxOuvS+MrjF+5MtF7r63lxKh0IuQBOKK0dKL3furURMiXlmasSVmB4RoAjvu//5OuXct0K8zB\ncE2WYn0QzCeTJ16vX2fi1SmUUBqMKywCIOQNNuk6LQCyFMM1hplcwfDPf1LBAGQ7Jl4N9tBD0uef\nZ7oVwO1ycqQvv8x0K8zB2jVZ5PXXJyayvvhi4v7rr2e2XcDOnWNXhPJ4pBs3Ju7v3JnZdpmO4RrD\nUIuM+erNN8dukmRZLFLmFIZrDObx8IuE+ePW9eT37Bm7Twnl7KXKTkLeYJs3S4FAplsBjPnOd6QP\nPhi7/8UX0le+MnZ/7Vrp9OnMtcsEhHyWCgbpIWF+siyJX/W5w8QrAGQpJl4N1tZGTx7zx+QxeYll\nDZxCyBusvz/TLQAm/O530h/+MPHz+DVeL10i5NOJkDfM5N7S0BC9JSDbEfKG6e+f+ifx+P0lSwh5\nZFZ/vzQ4OPHz+H3+4kwvQt4wnAyF+aq0dGLRvIsXpYKCie1IH6prAMBg09bJb9u2TX/84x+Vn5+v\njz76SJJ0+fJl/fCHP9TFixfl8Xh07NgxLVmyRJLU2tqqgwcP6sEHH9SBAwdUVVUlSert7VVTU5M+\n//xzVVdXa//+/bc3hjr5OVVUJH3ySaZbAYzhZKj0mVWd/IsvvqjALadNtrW1qbKyUufOnVNFRYXa\n2tokSaFQSEePHlUoFFIgENCOHTuSB25ublZHR4fC4bDC4fBtr4m5l8NgHOaR06fHVkUdXxl1/D4B\nn17Thvz69eu1dOnSKduOHz+uxsZGSVJjY6O6/reAeXd3t+rr65WbmyuPx6OioiL19PQoHo9rZGRE\nfr9fktTQ0JDcB+kzPjYPzAeTV6GUWIXSKTPq6w0NDcnlckmSXC6XhoaGJEkDAwN64oknks+zbVux\nWEy5ubmybTu53e12KxaLzabduAcPPZTpFgATTp2aerWy8funTmWmPdli1n/QW5Yly7Lmoi2SpJbx\nwm5J5eXlKqfu777s3DlxwsnFixO9ph/8YGKZVyATioqkSGTs/pUrUl7exHbcn2AwqODkWukUZhTy\nLpdLg4ODKigoUDweV35+vqSxHnpk/P+ipGg0Ktu25Xa7FZ30FR6NRuV2u+/42pNDHvfm7l+y53Xx\n4iOSpLfeGruNY4IbTjt+XLp5c+LnK1cmtuP+3NoB3rt3712fO6MSypqaGnV2dkqSOjs7VVtbm9x+\n5MgRXbt2TefPn1c4HJbf71dBQYEWL16snp4eJRIJHTp0KLkPZi+RSNzx9q1vee76GOC0J58cq6gZ\nr6oZv//kk5ltl/ES03juuecShYWFidzc3IRt24mDBw8m/vOf/yQqKioSXq83UVlZmfj000+Tz//l\nL3+ZePTRRxOPPfZYIhAIJLd/8MEHiVWrViUeffTRxE9+8pM7HusemoP78OMfZ7oFwIQHH0wkxhYY\nnnp78MFMt2zhS5WdrCcPwBGLFkmjo7dvz8uTRkacb49JuGgIgIy4+3zRTaUaLSYH7k+q7OR0GQBp\nc9fepUWQO4W1awxGoRIAQt5gKaqqgIzasyfTLcgejMkbjIslA9mBC3kDQJYi5AHAYIQ8ABiMkDcY\nk1sACHmDUUKJ+YrPpnOorgHgOCq/5hbVNQCQpQh5ADAYIQ8ABiPkDcbkFgBC3mCsXYP5ivJe51Bd\nYzAqGIDsQHUNAGQpQh4ADEbIA4DBCHmDMbkFgJA3GCWUmK/4bDqH6hoAjqPya25RXQMAWYqQBwCD\nEfIAYDBC3mBMbgFg4tVgTG7BKV//uvTpp+k9xtKl0uXL6T3GQpW2iVePx6PVq1errKxMfr9fknT5\n8mVVVlZq+fLlqqqq0vDwcPL5ra2t8nq9Ki4u1okTJ2ZzaADzyKefjnUo0nlL95eIqWYV8pZlKRgM\nqq+vT2fPnpUktbW1qbKyUufOnVNFRYXa2tokSaFQSEePHlUoFFIgENCOHTt08+bN2b8DAMBdzXpM\n/tY/EY4fP67GxkZJUmNjo7q6uiRJ3d3dqq+vV25urjwej4qKipJfDACA9Jh1T37jxo1au3atfv3r\nX0uShoaG5HK5JEkul0tDQ0OSpIGBAdm2ndzXtm3FYrHZHB4AMI2c2ex85swZFRYW6t///rcqKytV\nXFw85XHLsmRZ1l33T/UYpprpxNb9/BMzsQWYZ1YhX1hYKEn6xje+oWeeeUZnz56Vy+XS4OCgCgoK\nFI/HlZ+fL0lyu92KRCLJfaPRqNxu922v2TKp7q+8vFzl5eWzaaIxxie20onvXGBhCAaDCgaD9/Tc\nGZdQfvbZZ7px44YWLVqk//73v6qqqtKePXv05z//WQ8//LB++tOfqq2tTcPDw2pra1MoFNLWrVt1\n9uxZxWIxbdy4UZ988smU3jwllHfnRDkkJZeYKT6fmZUqO2fckx8aGtIzzzwjSfryyy/1ox/9SFVV\nVVq7dq3q6urU0dEhj8ejY8eOSZJ8Pp/q6urk8/mUk5Oj9vZ2hmsAIM04GWqBoKeE+YzPZ2alpScP\nAOMSsqQ0/2GemPRf3DtCHsCsWUo405NP7yGMxAJlAGAwQh4ADEbIA4DBGJNfIJjYAjAThPwCwcQW\ngJlguAYADEZPHsCcSPcJ7EuXpvf1TUXIA5i1+x1K5OxV5zBcAwAGI+QBwGCEPAAYjJAHAIMR8gAc\nt2dPpluQPVhPfoFgvW4Ad5MqO+nJA4DBCHkAMBghDwAG44zXBYTTxgHcL3ryC0Qicf+3+93v8uXM\nvkdkj5aWTLcge1BdYzCqZTBf8dmcW1TXAECWIuQBwGCEPAAYjJA3GKeOAyDkDUYFA+YrOiDOoboG\nABa4eVNdEwgEVFxcLK/Xq1/96ldOHhoAspJjIX/jxg3t3LlTgUBAoVBIhw8f1scff+zU4bNSMBjM\ndBOAO+Kz6RzHQv7s2bMqKiqSx+NRbm6unnvuOXV3dzt1+KzELxLmKz6bznEs5GOxmJYtW5b82bZt\nxWIxpw6flfg9AuDYAmVWulfXymKp/m0ta+8dtzPBDSek+mzu3Xvnz6bE53MuORbybrdbkUgk+XMk\nEpFt21Oes2bNGr4MHMK/M+YzPp/3Z82aNXd9zLESyi+//FKPPfaY/vKXv+ib3/ym/H6/Dh8+rBUr\nVjhxeADISo715HNycvTmm29q06ZNunHjhrZv307AA0CazauToQAAc4tlDQy0bds2uVwulZSUZLop\nwBSRSEQbNmzQypUrtWrVKh04cCDTTTIePXkD/fWvf1VeXp4aGhr00UcfZbo5QNLg4KAGBwdVWlqq\n0dFRPf744+rq6mLoNo3oyRto/fr1WsoFWzEPFRQUqLS0VJKUl5enFStWaGBgIMOtMhshDyAjLly4\noL6+Pq1bty7TTTEaIQ/AcaOjo3r22We1f/9+5eXlZbo5RiPkATjq+vXr2rJli55//nnV1tZmujnG\nI+QBOCaRSGj79u3y+Xx65ZVXMt2crEDIG6i+vl7f/va3de7cOS1btky/+c1vMt0kQJJ05swZvfvu\nuzp58qTKyspUVlamQCCQ6WYZjRJKADAYPXkAMBghDwAGI+QBwGCEPAAYjJAHAIMR8gBgMEIeAAxG\nyAOAwf4fLsH1RnH4liwAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f6b89e90110>"
       ]
      }
     ],
     "prompt_number": 398
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###_As far as we can tell_, they are no closer together that other heterozygotes:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from collections import defaultdict\n",
      "\n",
      "positions= defaultdict(list)\n",
      "for r in hets:\n",
      "    positions[r.CHROM].append(r.POS)\n",
      "\n",
      "def distance(rec, pos_map=positions):\n",
      "    diffs = [abs(rec.POS - p) for p in pos_map[rec.CHROM] if rec.POS - p != 0]\n",
      "    if diffs:\n",
      "        return min(diffs)\n",
      "    return np.Infinity\n",
      "    \n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 399
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "het_dists = array([distance(r) for r in hets])\n",
      "print mean(het_dists[weird_ones] == inf)\n",
      "print mean(het_dists[invert(weird_ones)] == inf)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.0126801152738\n",
        "0.018635502211\n"
       ]
      }
     ],
     "prompt_number": 400
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figure()\n",
      "boxplot( [\n",
      "            [d for d in het_dists[weird_ones] if d < inf], \n",
      "            [d for d in het_dists[invert(weird_ones)] if d < inf]\n",
      "        ]\n",
      ")\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9M1Pmdx/HnbLFpGiujtMI6Q4rKzLKjFMgtSC7dlJUD\nf7SrNuzBcVfELJtcaXqne5uN2z8s2OSqm8te9kcl2Vy4W7SJdGNaNU1BPLPj7u3dYmq1uWiTpSmL\nwzCSFQdOtwoi7/uDZeArCvgDhh+vRzLxy2e+3/E9k5l5zef7/Xw/X5eZGSIiIp97LN4FiIjI7KJg\nEBERBwWDiIg4KBhERMRBwSAiIg4KBhERcZhSMNy+fZucnByeffZZAK5evUpRURF+v5/i4mJ6e3tj\n6+7btw+fz0dGRgYtLS2x9rNnz5KZmYnP52Pnzp2x9v7+fsrKyvD5fOTn59PR0RG7r6GhAb/fj9/v\n5+DBgw/9ZEVEZHJTCoY33niDQCCAy+UCYP/+/RQVFfHxxx9TWFjI/v37Abh48SK/+MUvuHjxIs3N\nzfzgBz9g5DSJ6upq6uvraWtro62tjebmZgDq6+tJSkqira2NF198kd27dwPD4fOTn/yEM2fOcObM\nGfbu3esIIBERmR6TBkNnZye/+c1veOGFF2Jf8sePH6eyshKAyspKjh49CsCxY8coLy9n0aJFpKWl\nkZ6eTmtrK5FIhGvXrpGXlwfA9u3bY9uMfaySkhJOnToFwIkTJyguLsbtduN2uykqKoqFiYiITJ9J\ng+HFF1/kX/7lX3jssdFVu7u7SU5OBiA5OZnu7m4Aurq68Hq9sfW8Xi/hcHhcu8fjIRwOAxAOh0lN\nTQUgISGBxMREenp67vlYIiIyvSYMhl//+tcsX76cnJwc7jVzhsvliu1iEhGRuS9hojv/+7//m+PH\nj/Ob3/yGmzdv8n//939UVFSQnJzM5cuXSUlJIRKJsHz5cmC4JxAKhWLbd3Z24vV68Xg8dHZ2jmsf\n2ebSpUusWLGCwcFB+vr6SEpKwuPxEAwGY9uEQiHWr18/rsbs7Gx+//vfP9SLICKy0GRlZXH+/Pm7\n32lTFAwG7Tvf+Y6Zmb388su2f/9+MzPbt2+f7d6928zMLly4YFlZWdbf329/+tOfbNWqVTY0NGRm\nZnl5efbRRx/Z0NCQbdq0yZqamszM7MCBA/b973/fzMwOHz5sZWVlZmbW09NjK1eutGg0alevXo0t\n3+k+noJMQU1NTbxLELkrvTcfrYm+OyfsMdxpZJfRK6+8QmlpKfX19aSlpfHuu+8CEAgEKC0tJRAI\nkJCQQF1dXWyburo6duzYwY0bN9i8eTMbN24EoKqqioqKCnw+H0lJSTQ2NgKwbNky9uzZQ25uLgA1\nNTW43e77DUUREblPrs+TY85yuVz3PP4h96+2tpba2tp4lyEyjt6bj9ZE350681kcCgoK4l2CyF3p\nvTlz1GMQEVmA1GMQEZEpUzCIyJzw+uvxrmDhUDCIyJzw+Sw6MgMUDCIi4nBf5zGIiMyk118f7Smc\nPg0jA5O2bYNdu+JW1rynUUkiMicUFMCYWXLkIWlUkoiITJmCQUTmhG3b4l3BwqFdSSIiC5B2JYmI\nyJQpGERExEHBICIiDgoGERFxUDCIiIiDgkFERBwUDCIi4qBgEBERhwmD4ebNm6xbt47s7GwCgQA/\n+tGPgOFrr3q9XnJycsjJyaGpqSm2zb59+/D5fGRkZNDS0hJrP3v2LJmZmfh8Pnbu3Blr7+/vp6ys\nDJ/PR35+Ph0dHbH7Ghoa8Pv9+P1+Dh48+MietIiITMAm8dlnn5mZ2a1bt2zdunX2wQcfWG1trb32\n2mvj1r1w4YJlZWXZwMCAtbe32+rVq21oaMjMzHJzc621tdXMzDZt2mRNTU1mZnbgwAGrrq42M7PG\nxkYrKyszM7Oenh5btWqVRaNRi0ajseU7TeEpiIjIHSb67px0V9KXv/xlAAYGBrh9+zZLly4dCZRx\n6x47dozy8nIWLVpEWloa6enptLa2EolEuHbtGnl5eQBs376do5/PpXv8+HEqKysBKCkp4dSpUwCc\nOHGC4uJi3G43breboqIimpubHzoIRURkYpMGw9DQENnZ2SQnJ/PMM8+wZs0aAN566y2ysrKoqqqi\nt7cXgK6uLrxeb2xbr9dLOBwe1+7xeAiHwwCEw2FSU1MBSEhIIDExkZ6enns+loiITK9Jg+Gxxx7j\n/PnzdHZ28v777xMMBqmurqa9vZ3z58/z+OOP89JLL81ErSIiMgOmfAW3xMREvv3tb/Pb3/6WgpHL\nKAEvvPACzz77LDDcEwiFQrH7Ojs78Xq9eDweOjs7x7WPbHPp0iVWrFjB4OAgfX19JCUl4fF4CI65\nKkcoFGL9+vV3ra22tja2XFBQ4KhPREQgGAw6vlMnNNHBiU8//TR2wPfPf/6zPf300/af//mfFolE\nYuv867/+q5WXl5vZ6MHn/v5++9Of/mSrVq2KHXzOy8uzjz76yIaGhsYdfP7+979vZmaHDx92HHxe\nuXKlRaNRu3r1amz5fg6giIjI3U303TlhjyESiVBZWcnQ0BBDQ0NUVFRQWFjI9u3bOX/+PC6Xi5Ur\nV/L2228DEAgEKC0tJRAIkJCQQF1dHS6XC4C6ujp27NjBjRs32Lx5Mxs3bgSgqqqKiooKfD4fSUlJ\nNDY2ArBs2TL27NlDbm4uADU1Nbjd7vuPSRERuS+6UI+IyAKkC/WIiMiUKRhERMRBwSAiIg4KBhER\ncVAwiIiIg4JBREQcFAwiIuKgYBAREQcFg4iIOCgYRETEQcEgIiIOCgYREXFQMIiIiIOCQUREHBQM\nIiLioGAQEREHBYOIiDgoGERExEHBICIiDhMGw82bN1m3bh3Z2dkEAgF+9KMfAXD16lWKiorw+/0U\nFxfT29sb22bfvn34fD4yMjJoaWmJtZ89e5bMzEx8Ph87d+6Mtff391NWVobP5yM/P5+Ojo7YfQ0N\nDfj9fvx+PwcPHnxkT1pERCZgk/jss8/MzOzWrVu2bt06++CDD+zll1+2V1991czM9u/fb7t37zYz\nswsXLlhWVpYNDAxYe3u7rV692oaGhszMLDc311pbW83MbNOmTdbU1GRmZgcOHLDq6mozM2tsbLSy\nsjIzM+vp6bFVq1ZZNBq1aDQaW77TFJ6CiIjcYaLvzkl3JX35y18GYGBggNu3b7N06VKOHz9OZWUl\nAJWVlRw9ehSAY8eOUV5ezqJFi0hLSyM9PZ3W1lYikQjXrl0jLy8PgO3bt8e2GftYJSUlnDp1CoAT\nJ05QXFyM2+3G7XZTVFREc3PzIw1FEREZb9JgGBoaIjs7m+TkZJ555hnWrFlDd3c3ycnJACQnJ9Pd\n3Q1AV1cXXq83tq3X6yUcDo9r93g8hMNhAMLhMKmpqQAkJCSQmJhIT0/PPR9LRESmV8JkKzz22GOc\nP3+evr4+NmzYwHvvvee43+Vy4XK5pq3AqaitrY0tFxQUUFBQELdaRERmo2AwSDAYnNK6kwbDiMTE\nRL797W9z9uxZkpOTuXz5MikpKUQiEZYvXw4M9wRCoVBsm87OTrxeLx6Ph87OznHtI9tcunSJFStW\nMDg4SF9fH0lJSXg8HseTCIVCrF+//q61jQ0GEREZ784fzXv37r3nuhPuSrpy5UpsxNGNGzc4efIk\nOTk5bNmyhYaGBmB45NC2bdsA2LJlC42NjQwMDNDe3k5bWxt5eXmkpKSwZMkSWltbMTMOHTrE1q1b\nY9uMPNaRI0coLCwEoLi4mJaWFnp7e4lGo5w8eZINGzY84EsiIiJTNWGPIRKJUFlZydDQEENDQ1RU\nVFBYWEhOTg6lpaXU19eTlpbGu+++C0AgEKC0tJRAIEBCQgJ1dXWx3Ux1dXXs2LGDGzdusHnzZjZu\n3AhAVVUVFRUV+Hw+kpKSaGxsBGDZsmXs2bOH3NxcAGpqanC73dP2QoiIyDDX58OW5iyXy8Ucfwoi\nIjNuou9OnfksIiIOCgYREXFQMIiIiIOCQUREHBQMIiLioGAQEREHBYOIiDgoGERExEHBICIiDgoG\nERFxUDCIiIiDgkFERBwUDCIyJ0zxGjPyCCgYRGROUDDMHAWDiIg4TPnSniIiMy0YHO0pjL0SZUHB\n8E2mhy7UIyJzwo4d8M478a5i/tCFekRkzvvkk3hXsHAoGERkTkhLi3cFC8ekwRAKhXjmmWdYs2YN\na9eu5c033wSgtrYWr9dLTk4OOTk5NDU1xbbZt28fPp+PjIwMWlpaYu1nz54lMzMTn8/Hzp07Y+39\n/f2UlZXh8/nIz8+no6Mjdl9DQwN+vx+/38/BgwcfyZMWkbkhGITa2uFbQ8PoskYoTTObRCQSsXPn\nzpmZ2bVr18zv99vFixettrbWXnvttXHrX7hwwbKysmxgYMDa29tt9erVNjQ0ZGZmubm51traamZm\nmzZtsqamJjMzO3DggFVXV5uZWWNjo5WVlZmZWU9Pj61atcqi0ahFo9HY8lhTeAoiMg/U1MS7gvll\nou/OSXsMKSkpZGdnA7B48WKefPJJwuHwSKiMW//YsWOUl5ezaNEi0tLSSE9Pp7W1lUgkwrVr18jL\nywNg+/btHD16FIDjx49TWVkJQElJCadOnQLgxIkTFBcX43a7cbvdFBUV0dzc/NBhKCIi93Zfxxg+\n+eQTzp07R35+PgBvvfUWWVlZVFVV0dvbC0BXVxderze2jdfrJRwOj2v3eDyxgAmHw6SmpgKQkJBA\nYmIiPT0993wsEVl4NDx15kz5PIbr16/z3HPP8cYbb7B48WKqq6v58Y9/DMCePXt46aWXqK+vn7ZC\nJ1JbWxtbLigooEDvIJF5Rx/rhxMMBglO8eDMlILh1q1blJSU8L3vfY9t27YBsHz58tj9L7zwAs8+\n+yww3BMIhUKx+zo7O/F6vXg8Hjo7O8e1j2xz6dIlVqxYweDgIH19fSQlJeHxeBxPJBQKsX79+nH1\njQ0GEZmfgkGFw8O480fz3rFnDN5h0l1JZkZVVRWBQIBdu3bF2iORSGz5V7/6FZmZmQBs2bKFxsZG\nBgYGaG9vp62tjby8PFJSUliyZAmtra2YGYcOHWLr1q2xbRoaGgA4cuQIhYWFABQXF9PS0kJvby/R\naJSTJ0+yYcOG+3gpRGS+0MltM2fSHsOHH37Iz3/+c77xjW+Qk5MDwE9/+lMOHz7M+fPncblcrFy5\nkrfffhuAQCBAaWkpgUCAhIQE6urqcLlcANTV1bFjxw5u3LjB5s2b2bhxIwBVVVVUVFTg8/lISkqi\nsbERgGXLlrFnzx5yc3MBqKmpwe12P/pXQURmPZ3gNnM0JYaIzFp3zpVUUzO8rLmSHp6mxBARkSlT\nj0FE5oS0NO1OepTUYxCROU+HF2eOrscgIrPW2GMMv//98DxJoGMM003BICKz1p0BoFOWZoZ2JYmI\niIOCQRw0nbHMVtp1NHMUDOKgYJDZSsEwcxQMIiLioIPPMu7s0hEa+SGyMCkYRCM/RMRBu5JERMRB\nwSAO2nUkIgoGcVAwyGylEXMzR8EgInOCgmHmKBhERMRBo5JEZNbSUOr40PUYRGROKCjQ7qRHSddj\nEBGRKZs0GEKhEM888wxr1qxh7dq1vPnmmwBcvXqVoqIi/H4/xcXF9Pb2xrbZt28fPp+PjIwMWlpa\nYu1nz54lMzMTn8/Hzp07Y+39/f2UlZXh8/nIz8+no6Mjdl9DQwN+vx+/38/BgwcfyZMWkbkhGBw+\n4bK2Fk6fHl1Wz2Ga2SQikYidO3fOzMyuXbtmfr/fLl68aC+//LK9+uqrZma2f/9+2717t5mZXbhw\nwbKysmxgYMDa29tt9erVNjQ0ZGZmubm51traamZmmzZtsqamJjMzO3DggFVXV5uZWWNjo5WVlZmZ\nWU9Pj61atcqi0ahFo9HY8lhTeAoiMg/U1MS7gvllou/OSXsMKSkpZGdnA7B48WKefPJJwuEwx48f\np7KyEoDKykqOHj0KwLFjxygvL2fRokWkpaWRnp5Oa2srkUiEa9eukZeXB8D27dtj24x9rJKSEk6d\nOgXAiRMnKC4uxu1243a7KSoqorm5+ZEGo4iION3XMYZPPvmEc+fOsW7dOrq7u0lOTgYgOTmZ7u5u\nALq6uvB6vbFtvF4v4XB4XLvH4yEcDgMQDodJTU0FICEhgcTERHp6eu75WCKy8GgU0syZ8nDV69ev\nU1JSwhtvvMFXvvIVx30ulwuXy/XIi5uq2jGzvhUUFFCgd5DIvKOP9cMJBoMEp3hwZkrBcOvWLUpK\nSqioqGDbtm3AcC/h8uXLpKSkEIlEWL58OTDcEwiFQrFtOzs78Xq9eDweOjs7x7WPbHPp0iVWrFjB\n4OAgfX19JCUl4fF4HE8kFAqxfv36cfXVajpQEZEJ3fmjee/YE0PuMOmuJDOjqqqKQCDArl27Yu1b\ntmyhoaEBGB45NBIYW7ZsobGxkYGBAdrb22lrayMvL4+UlBSWLFlCa2srZsahQ4fYunXruMc6cuQI\nhYWFABQXF9PS0kJvby/RaJSTJ0+yYcOG+3w5RETkvkx25PqDDz4wl8tlWVlZlp2dbdnZ2dbU1GQ9\nPT1WWFhoPp/PioqKHKOF/vmf/9lWr15tTzzxhDU3N8faf/vb39ratWtt9erV9g//8A+x9ps3b9pf\n//VfW3p6uq1bt87a29tj9/37v/+7paenW3p6ur3zzjv3dWRdRETubqLvTp35LCKyAOnMZxERmTIF\ng4iIOCgYRETEQcEgIiIOCgYREXFQMIiIiIOCQUREHBQM4qB57kVEwSAOCgYRUTCIiIjDlKfdlvkr\nGBztKYydcLGgQFMdiyxECgYZFwCaxVxkYdOuJBGZE15/Pd4VLBwKBnHQriOZrT6/RLzMAAWDOCgY\nRETHGERk1nr99dGewunToz9ctm2DMReUlEdMF+oRkTmhoEDn2TxKulCPiIhMmYJBROaEbdviXcHC\nMWkwPP/88yQnJ5OZmRlrq62txev1kpOTQ05ODk1NTbH79u3bh8/nIyMjg5aWllj72bNnyczMxOfz\nsXPnzlh7f38/ZWVl+Hw+8vPz6ejoiN3X0NCA3+/H7/dz8ODBh36yIjJ3ZWfHu4IFxCbx/vvv2+9+\n9ztbu3ZtrK22ttZee+21ceteuHDBsrKybGBgwNrb22316tU2NDRkZma5ubnW2tpqZmabNm2ypqYm\nMzM7cOCAVVdXm5lZY2OjlZWVmZlZT0+PrVq1yqLRqEWj0djynabwFERkHqipiXcF88tE352T9hie\nfvppli5derdAGdd27NgxysvLWbRoEWlpaaSnp9Pa2kokEuHatWvk5eUBsH37do5+PtTg+PHjVFZW\nAlBSUsKpU6cAOHHiBMXFxbjdbtxuN0VFRTQ3Nz9o/omIyBQ98HDVt956i4MHD/LUU0/x2muv4Xa7\n6erqIj8/P7aO1+slHA6zaNEivF5vrN3j8RAOhwEIh8OkpqYOF5OQQGJiIj09PXR1dTm2GXksEVk4\nNI9XfDxQMFRXV/PjH/8YgD179vDSSy9RX1//SAu7H7VjJvcpKCigQO8YkXlB83g9OsFgkOAUx/s+\nUDAsX748tvzCCy/w7LPPAsM9gVAoFLuvs7MTr9eLx+Ohs7NzXPvINpcuXWLFihUMDg7S19dHUlIS\nHo/H8SRCoRDr16+/az21ereIiEzozh/Ne8d2we7wQMNVI5FIbPlXv/pVbMTSli1baGxsZGBggPb2\ndtra2sjLyyMlJYUlS5bQ2tqKmXHo0CG2bt0a26ahoQGAI0eOUFhYCEBxcTEtLS309vYSjUY5efIk\nGzZseJBy5T7oBCKZrbQjYOZM2mMoLy/n9OnTXLlyhdTUVPbu3UswGOT8+fO4XC5WrlzJ22+/DUAg\nEKC0tJRAIEBCQgJ1dXW4XC4A6urq2LFjBzdu3GDz5s1s3LgRgKqqKioqKvD5fCQlJdHY2AjAsmXL\n2LNnD7m5uQDU1NTgdrun5UWQUcGgPoAyO+l9OXM0JYY41NZqP67IQjDRd6cm0RON/BARBwWDaOSH\niDhoriRx+OSTeFcgIvGmYBAREQcFgzikpcW7AhGJNx1jEB18FhEHBYPo4LOIOGhXkoiIOCgYxEG7\njmS20nQtM0fBICJzgoJh5igYxEEfPpmtdI7NzNHBZ3H46KN4VyAyauyIuYaG0eHUGjE3vTSJnowb\nrlpTM7ysD5/MJmlp6jU8ShN9dyoYxMHtht7eeFchMkw/WqaPZleVCb3+Ohw9Orzc1zf6gdu2DXbt\niltZIhIn6jGIfpXJnPDFL8LAQLyrmD+0K0mmTPtxZTbRj5bpo11JMmUZGfGuQGTU+fPOIdQjy263\ngmE6KRjE4ZVX4l2ByKjs7NHBEKdPj4ZBdnbcSloQJj3B7fnnnyc5OZnMzMxY29WrVykqKsLv91Nc\nXEzvmGEs+/btw+fzkZGRQUtLS6z97NmzZGZm4vP52LlzZ6y9v7+fsrIyfD4f+fn5dHR0xO5raGjA\n7/fj9/s5ePDgQz9ZmZx+hYkINon333/ffve739natWtjbS+//LK9+uqrZma2f/9+2717t5mZXbhw\nwbKysmxgYMDa29tt9erVNjQ0ZGZmubm51traamZmmzZtsqamJjMzO3DggFVXV5uZWWNjo5WVlZmZ\nWU9Pj61atcqi0ahFo9HY8p2m8BREZB740pfiXcH8MtF356Q9hqeffpqlS5c62o4fP05lZSUAlZWV\nHP18rOOxY8coLy9n0aJFpKWlkZ6eTmtrK5FIhGvXrpGXlwfA9u3bY9uMfaySkhJOnToFwIkTJygu\nLsbtduN2uykqKqK5ufmRhKGIzD2JifGuYOF4oGMM3d3dJCcnA5CcnEx3dzcAXV1d5Ofnx9bzer2E\nw2EWLVqE1+uNtXs8HsLhMADhcJjU1NThYhISSExMpKenh66uLsc2I48lIgvH2HNsurt1js1MeeiD\nzy6XC5fL9ShqeWC1Y64sU1BQQIF2lIvMC7t2jQbA4sWa5PFhBINBglN8AR8oGJKTk7l8+TIpKSlE\nIhGWL18ODPcEQqFQbL3Ozk68Xi8ej4fOzs5x7SPbXLp0iRUrVjA4OEhfXx9JSUl4PB7HkwiFQqxf\nv/6u9dTqkmMi897Nm/GuYG6780fz3rHX8b3DA027vWXLFhoaGoDhkUPbtm2LtTc2NjIwMEB7eztt\nbW3k5eWRkpLCkiVLaG1txcw4dOgQW7duHfdYR44cobCwEIDi4mJaWlro7e0lGo1y8uRJNmzY8CDl\nisgc9frroyez3b49uvz66/Gta76btMdQXl7O6dOnuXLlCqmpqfzkJz/hlVdeobS0lPr6etLS0nj3\n3XcBCAQClJaWEggESEhIoK6uLrabqa6ujh07dnDjxg02b97Mxo0bAaiqqqKiogKfz0dSUhKNjY0A\nLFu2jD179pCbmwtATU0Nbrd7Wl4EGRUMasiqzB46jyE+NCWGONTWDt9EZoOxB59Pn4ZvfWt4WQef\nH56mxBCROUk9hvhQMMi4icpGaKIyibcjR+DXvx79+513hv+9ckXvzemkYBBHABw9ql1JMnukp49e\nzrOjY3Q5PT1eFS0MOsYgDpp2W2Yrlwv0UX90dIxBJjR2V1JHx2iPQbuSJN5++EPnrqSRHsN3vgM/\n+1lcSloQFAwiMmtpV1J8aFeSOKSkwOXL8a5CZNjYHkNHB3z968PL6jE8PF3aU6Zs2TK4ejXeVYiM\np2MMj5aOMciExh5jiEZ1jEFmj7HvTdB7c6YoGERk1hobAPv2aSj1TFEwiMisNXZKjIEBXY9hpigY\nRGTWGns9BpdL12OYKQoGEZm1xvYYQD2GmaJgEN54A957b/Tvkbnuf/97HeATWYgUDILHAyOXuujr\nG132eOJXkwhAfT384Q+jf//Xfw3/29OjHsN00nkM4qCx4jKb6HoM02ei784HurSniIjMX+oxiIN6\nDDKbLFs2fNLlnZYu1Rn6D2vaegxpaWl84xvfICcnh7y8PACuXr1KUVERfr+f4uJiekcuvwTs27cP\nn89HRkYGLS0tsfazZ8+SmZmJz+dj586dsfb+/n7Kysrw+Xzk5+fT0dHxMOXKPXzzm/ClLw3fYHT5\nm9+Mb10iv/wl1NQM32B0+Ze/jG9d891DBYPL5SIYDHLu3DnOnDkDwP79+ykqKuLjjz+msLCQ/fv3\nA3Dx4kV+8YtfcPHiRZqbm/nBD34QS6vq6mrq6+tpa2ujra2N5uZmAOrr60lKSqKtrY0XX3yR3bt3\nP0y5cg9f+9rdg+FrX4tvXSLnzzunxRhZPn8+fjUtCPYQ0tLS7MqVK462J554wi5fvmxmZpFIxJ54\n4gkzM/vpT39q+/fvj623YcMG+5//+R/r6uqyjIyMWPvhw4ft7//+72PrfPTRR2ZmduvWLfvqV786\nroaHfApiZu+9Z1ZTM3yD0eX33otnVSJma9eafeELwzcYXV67Nt6VzX0TfXc+dI/hr/7qr3jqqaf4\nt3/7NwC6u7tJTk4GIDk5me7ubgC6urrwer2xbb1eL+FweFy7x+MhHA4DEA6HSU1NBSAhIYHExESu\nasfiI6dfZTJbJSZCQsLwDUaXExPjW9d891DnMXz44Yc8/vjjfPrppxQVFZGRkeG43+Vy4XK5HqpA\nmX5//KPzcp4jy3/8YzyqEZF4e6hgePzxxwH42te+xne/+13OnDlDcnIyly9fJiUlhUgkwvLly4Hh\nnkAoFIpt29nZidfrxePx0NnZOa59ZJtLly6xYsUKBgcH6evrY9myZePqqB0z5WJBQQEFOl33vjz3\nHHz1q8PLe/fCjh3Dy3oZJd7C4eHJ80aMLH++U0HuQzAYJDjFyaYeOBj+/Oc/c/v2bb7yla/w2Wef\n0dLSQk1NDVu2bKGhoYHdu3fT0NDAtm3bANiyZQt/+7d/yz/90z8RDodpa2sjLy8Pl8vFkiVLaG1t\nJS8vj0OHDvGP//iPsW0aGhrIz8/nyJEjFBYW3rWWWs3FKzIvZWePDlft64MlS0bb5f7c+aN57969\n91z3gYOhu7ub7373uwAMDg7yd3/3dxQXF/PUU09RWlpKfX09aWlpvPvuuwAEAgFKS0sJBAIkJCRQ\nV1cX282UX9KnAAAEXElEQVRUV1fHjh07uHHjBps3b2bjxo0AVFVVUVFRgc/nIykpicbGxgctVyZw\n5IjzguvvvDP875Ur6jVIfH3rW6PBcPr0aCCMnAEt00MnuInjwPPevaNjxnWVLIm3lSuHr/UMwyde\njhyy/PrXob09fnXNB7q0p0xIPQaZrRYvhsc+Hzt5+/bo8uLF8atpIVAwCOfPw+XLo3+PLGu4qsTb\n9eswNDT698jy9evxqWehUDBI7MxngP7+0WWd+Szxph5DfGh2VRERcVCPQfj0U7h5c/TvkeVPP41P\nPSIjLl507kq6fXu0XaaPgkHo64PBwdG/R5b7+uJTj8iIe02coAkVppeCQWLz0cDwL7KRZc1HI/E2\n0kOYars8GjrGICKz1r0OMuvg8/RSMIiIiIPOfBYSEu7eNf/CF5zHHkRmmt6b02faLu0p88PIzKpT\nbReZKTrGEB/qMciEIzz00ko86b05fdRjkAnpAJ/MVhquGh/qMQhf+ILzJKIRjz2mLrvEl3oM00c9\nBpnQ3UJhonYRmd8UDBKbmGyq7SIyv+mjL+oxiIiDgkFERBwUDCIi4jDrg6G5uZmMjAx8Ph+vvvpq\nvMsREZn3ZnUw3L59mx/+8Ic0Nzdz8eJFDh8+zB/+8Id4lzXPBeNdgCxwLpcrdoPbgH1+C45Zvu1Y\nz6UTGx6pWR0MZ86cIT09nbS0NBYtWsTf/M3fcOzYsXiXNS84P1AD3P3DN6APnsw4Mxtz+wJmLsxc\nQDC2PNxujps8OrM6GMLhMKmpqbG/vV4v4XA4jhXNH84P1RfHfPgY8+H7oj54Mm2WLRs+gW2qN7i/\n9V2u4f9D7t+svlCPfqU+nGXLIBq9/+3u92VfuhSuXr3//0cWtmj0/s5edrnu/2xnfYU8mFkdDB6P\nh1AoFPs7FArh9Xod62RlZSlAHrm997V2NKoPoDyY+33fuFz39958kP9jocjKyrrnfbN6rqTBwUGe\neOIJTp06xYoVK8jLy+Pw4cM8+eST8S5NRGTemtU9hoSEBH72s5+xYcMGbt++TVVVlUJBRGSazeoe\ng4iIzLxZPSpJZs7zzz9PcnIymZmZ8S5FxCEUCvHMM8+wZs0a1q5dy5tvvhnvkuY99RgEgA8++IDF\nixezfft2/vd//zfe5YjEXL58mcuXL5Odnc3169f5i7/4C44ePardytNIPQYB4Omnn2bp0qXxLkNk\nnJSUFLKzswFYvHgxTz75JF1dXXGuan5TMIjInPHJJ59w7tw51q1bF+9S5jUFg4jMCdevX+e5557j\njTfeYLEuSD6tFAwiMuvdunWLkpISvve977Ft27Z4lzPvKRhEZFYzM6qqqggEAuzatSve5SwICgYB\noLy8nL/8y7/k448/JjU1lf/4j/+Id0kiAHz44Yf8/Oc/57333iMnJ4ecnByam5vjXda8puGqIiLi\noB6DiIg4KBhERMRBwSAiIg4KBhERcVAwiIiIg4JBREQcFAwiIuKgYBAREYf/B76/OS1Do7m9AAAA\nAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f6b8be73650>"
       ]
      }
     ],
     "prompt_number": 401
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pos_within= defaultdict(list)\n",
      "for b,r in zip(weird_ones, hets):\n",
      "    if b:\n",
      "        pos_within[r.CHROM].append(r.POS)\n",
      "    \n",
      "within_m4751 = array([distance(h, pos_within) for b,h in zip(weird_ones, hets) if b])\n",
      "\n",
      "print mean([d for d in within_m4751 if d < inf])\n",
      "print mean([d for d in het_dists[weird_ones] if d < inf])\n",
      "print mean([d for d in het_dists[invert(weird_ones)] if d < inf])\n",
      "print mean(within_m4751 == Inf)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "13592.7309309\n",
        "5017.5931115\n",
        "6043.52977148\n",
        "0.0403458213256\n"
       ]
      }
     ],
     "prompt_number": 402
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Are there just more close ones?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print mean(het_dists[weird_ones] < 10000)\n",
      "print mean(het_dists[invert(weird_ones)] < 10000)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.857060518732\n",
        "0.844598862919\n"
       ]
      }
     ],
     "prompt_number": 403
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That's a weirdly high number, what the hell is goign on"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print mean(het_dists[weird_ones] < 100)\n",
      "print mean(het_dists[invert(weird_ones)] < 100)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.689337175793\n",
        "0.668982943778\n"
       ]
      }
     ],
     "prompt_number": 404
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print mean(het_dists[weird_ones] < 5)\n",
      "print mean(het_dists[invert(weird_ones)] < 5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.61613832853\n",
        "0.55843335439\n"
       ]
      }
     ],
     "prompt_number": 405
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That seems like way too many too close?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "spaced_hets = [r for r,b in zip(hets,het_dists < 100 ) if b]\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 414
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 414
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 414
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 414
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}