gistfile1.txt

{
 "metadata": {
  "name": "",
  "signature": "sha256:97307bd8345790dff116e1c0e235202ead1dd71febdf54759ace1b8bc3e5ccb1"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Write linear model to generate gene expression data from SNPs\n",
      "%pylab inline\n",
      "\n",
      "import numpy as np\n",
      "import GPy\n",
      "import pylab as pb\n",
      "#Columns are samples\n",
      "X = np.genfromtxt('ALLsnpspanama.csv', delimiter=',')[1:,1:]\n",
      "#Y = WX + N  , W is sparse\n",
      "\n",
      "SNPsconsider=np.array(range(10))*1000\n",
      "X=X[SNPsconsider,:]\n",
      "ng=1\n",
      "W=np.zeros((ng, X.shape[0]))\n",
      "\n",
      "\n",
      "W[:,[1]]=20\n",
      "W[:,[5]]=30\n",
      "W[:,[3]]=50\n",
      "\n",
      "print W.shape\n",
      "print X.shape\n",
      "print np.dot(W,X).shape\n",
      "\n",
      "Y=np.dot(W,X) + np.random.normal(size=(ng,182))\n",
      "\n",
      "print np.dot(W,X).shape\n",
      "print np.random.normal(size=(1,182)).shape\n",
      "pb.matshow(Y)\n",
      "pb.colorbar()\n",
      "print Y.shape\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n",
        "(1, 10)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "(10, 182)\n",
        "(1, 182)\n",
        "(1, 182)\n",
        "(1, 182)\n",
        "(1, 182)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "WARNING: pylab import has clobbered these variables: ['step']\n",
        "`%matplotlib` prevents importing * from pylab and numpy\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAyMAAACBCAYAAADT72W9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFCtJREFUeJzt3X9QlNe9x/HPoxiNxhp7i4thSXCCBolU6SXx/hGr1gCa\nFqqTjBebH1zFpKPjvWmwo/FXizYRMLYpxSbXVpvBdAb1TooYR3esmS5jM3MljWDsoNFJ4IqU4BjE\nGFGpy94/tLuch92FLcJKfL9mzsyeZ5/nOYfd5+zul3POcyyv1+sVAAAAAPSzQZGuAAAAAIA7E8EI\nAAAAgIggGAEAAAAQEQQjAAAAACKCYAQAAABARBCMAAAAAJAk3W1ZsgKkr3/9631SnsWtfQEAAABI\nkmVZeiXA9rWS+iJsiLrlZwQAAAAwYN3dj2URjAAAAADwIRgBAAAAEBFf68eyCEYAAAAA+NAzAgAA\nACAiCEYAAAAARATDtAAAAABERH/2jLDoIQAAAACfuwOkQIqLi5WcnKxJkyapuLhYktTS0qK0tDRN\nmDBB6enpam1tDVkWwQgAAAAAn54EI3/961+1bds2ffDBBzp27Jj27dunTz75RIWFhUpLS9OpU6c0\na9YsFRYWhiyLYAQAAACAz9eiuia7kydPaurUqRo2bJgGDx6s6dOn65133tHevXuVk5MjScrJydGe\nPXtClkUwAgAAAMDn7qFdk92kSZN0+PBhtbS0qK2tTfv379fZs2fV3Nwsh8MhSXI4HGpubg5ZFhPY\nAQAAAPjcPUxy//1GCiYxMVErV65Uenq6RowYoSlTpmjw4MHGPpZlybKskGVZXq/XeysqDQAAAGBg\nsyxL3vsDbD8jhQob1qxZI6fTqeLiYrndbsXExKipqUkzZ87UyZMngx7HMC0AAAAAfkMDpADOnTsn\nSTpz5oz+8Ic/6Ac/+IGysrJUWloqSSotLdXcuXNDFkXPCAAAAABJN3tGkgNsP961Z+Tb3/62Pv/8\ncw0ZMkSvv/66Zs6cqZaWFs2fP19nzpxRfHy8du/erXvvvTd4eQQjAAAAAKSbwci/Bdj+v6GHaf2z\nmMAOAAAAwC/IsKy+QDACAAAAwG9Y/xVFMAIAAADAj54RAAAAABExov+KIhgBAAAA4EfPCAAAAICI\nYM4IAAAAgIhgmBYAAACAiOjBMK2PP/5Y2dnZvvynn36qDRs26MKFC9q2bZuio6MlSQUFBZo9e3bQ\n87DoIQAAAABJNxc93BBg+0+CL3rY0dGh2NhYVVVV6Xe/+51GjhypvLy8HpVHzwgAAAAAvzAnsB86\ndEgJCQmKi4uT1+sNa6X2QWFWDQAAAMBX2YgAKYSdO3dqwYIFkm70rJSUlGjy5MnKzc1Va2tryGMZ\npgUAAABA0s1hWtsk90nJ/bF/+/p3Aw/Tam9vV2xsrGpraxUdHa1z58755ousW7dOTU1N2r59e/Dy\nCEYAAAAASDeDkbcDbH82cDBSUVGhN998Uy6Xq8tz9fX1yszM1PHjx4OWxzAtAAAAAH5hDNMqKyvz\nDdGSpKamJt/j8vJyJScnhyyKnhEAAAAAkm72jOwPsP2Jrj0jly9f1gMPPKC6ujqNHDlSkvTcc8+p\npqZGlmVp3Lhx2rp1qxwOR/DyCEYAAAAASDeDkfcCbJ8V/Na+vcGtfQEAAAD4sQI7AAAAgIgY1n9F\nEYwAAAAA8Atz0cPeIBgBAAAA4EcwAgAAACAi7um/oghGAAAAAPhcp2cEAAAAQCRcGxpoXfSOPinr\nlqzA7nK5lJiYqPHjx6uoqOhWnBIYUOLj4/XNb35TKSkpevTRRyVJLS0tSktL04QJE5Senq7W1tYI\n1xLoG4sWLZLD4TBW2Q11/RcUFGj8+PFKTEzUwYMHI1FloM8Eag/5+flyOp1KSUlRSkqKDhw44HuO\n9oDb0ZWhw7ukQFpbW/XUU09p4sSJSkpK0pEjR8L+/dPrYMTj8WjZsmVyuVyqra1VWVmZTpw40dvT\nAgOKZVlyu92qrq5WVVWVJKmwsFBpaWk6deqUZs2apcLCwgjXEugbCxculMvlMrYFu/5ra2u1a9cu\n1dbWyuVyaenSpero6Jv/tgGREKg9WJalvLw8VVdXq7q6WnPmzJFEe8Dt65qGdkmBvPjii3riiSd0\n4sQJffTRR0pMTAz790+vg5GqqiolJCQoPj5eQ4YMUXZ2tioqKnp7WmDAsa9KunfvXuXk5EiScnJy\ntGfPnkhUC+hz06ZN0+jRo41twa7/iooKLViwQEOGDFF8fLwSEhJ8ATzwVRCoPUiBV66mPeB21a67\nuiS7ixcv6vDhw1q0aJEkKSoqSqNGjQr790+vg5HGxkbFxcX58k6nU42Njb09LTCgWJalxx9/XKmp\nqfrtb38rSWpubpbD4ZAkORwONTc3R7KKQL8Kdv3/7W9/k9Pp9O3HdwbuFCUlJZo8ebJyc3N9w1Zo\nD7hdXdNdXZJdXV2doqOjtXDhQn3rW9/S888/r8uXL4f9+6fXwYhlWb09BTDgvf/++6qurtaBAwf0\n61//WocPHzaetyyLtoI7VnfXP20DX3VLlixRXV2dampqNHbsWC1fvjzovrQH3A6uaHiXZHf9+nUd\nPXpUS5cu1dGjRzVixIguQ7J68vun13fTio2NVUNDgy/f0NBgRPnAnWDs2LGSpOjoaM2bN09VVVVy\nOBz67LPPFBMTo6amJo0ZMybCtQT6T7Dr3/6dcfbsWcXGxkaqmkC/6Pz5v3jxYmVmZkqiPeD2dU1D\nddR9SUfdXwbdx+l0yul06pFHHpEkPfXUUyooKFBMTExYv3963TOSmpqq06dPq76+Xu3t7dq1a5ey\nsrJ6e1pgwGhra9OlS5ckSZcvX9bBgweVnJysrKwslZaWSpJKS0s1d+7cSFYT6FfBrv+srCzt3LlT\n7e3tqqur0+nTp313oAO+qpqamnyPy8vLfXfaoj3gdtWuuzRpxr/oufwHfMkuJiZGcXFxOnXqlCTp\n0KFDevjhh5WZmRnW759e94xERUVpy5YtysjIkMfjUW5uriZOnNjb0wIDRnNzs+bNmyfpRpfl008/\nrfT0dKWmpmr+/Pnavn274uPjtXv37gjXFOgbCxYsUGVlpc6fP6+4uDht2LBBL7/8csDrPykpSfPn\nz1dSUpKioqL0xhtvMCwFXyn29rB+/Xq53W7V1NTIsiyNGzdOW7dulUR7wO2rLcCwrEBKSkr09NNP\nq729XQ8++KDeeusteTyesH7/WN5At3cAAAAAcMexLEsV3vQu279vHQx4V7jeYgV2AAAAAD7B1hXp\nCwQjAAAAAHzadHe/ldXtBHaXy6XExESNHz9eRUVF/VEnAAAAABHSrqFdUl8JGYx4PB4tW7ZMLpdL\ntbW1Kisr04kTJ/qsMgAAAAAiqycrsN8qIYdpVVVVKSEhQfHx8ZKk7OxsVVRUGHfLemDG/TpT2RDk\nDAAAAMCd6gF5vfWRrkTYbps5I42NjYqLi/PlnU6njhw5YuxzprJBa72rtGPG7/Wc+xm98t5G4/m1\ns1Yb+VfeMZ9f9eRPjHyM9TMj39Lp8RVvvvHcJmukkS+UuaLpy/tsM/5tf+3yjFeM/OaV68y6Fpnn\nWzdvs3mCZxRSyZOLjfx5a7uRz+/89A9tBw8zs1byz80N2XlG1htr3grQ+nmxuf/O/zKyv8k2938h\n2/Za2da48T5i7l/wkx8Z+dWlr5sH7On02PY6eR8yz7V+0gojn7/QHA64/C3zfbrfMt+n6+bp9YUt\n7/GuMvL2a1S/N7Nr3zKvWY/twtn4vv8atU6ar9vbi2/8ba9KWiPpY1vZP/tLgVmXVPMaa5a5MNCW\nSvO1WTF9vZEvOpJv5FdONfPxlplv8Jqv3ca/mO3NesT/WntXrzWeW/2qeWxB5gYjrx+bWe9+2zW5\n6W0j/4aeNfJL3Obx1pe2a/J7tmGi7pVG1v7abHr+p+b+2171Pcz3mlfJT9/fZO5bbmbXbTbfR7tX\ndpnX1Ip/N+uSYHsfXthp+9vuNbO5GVuM/BidM/Ibl5nv2+ot/vdmsK1FRMlj5PMth5H3vmpbCdos\nSit/mW/kN1n3mTsYn9LSKu9ls641N+o6Y7Hk3iatnGKez/7a1Nmu0Uu6x/e45KT5nuuXZtZaYntd\np7xpq2uzkVvmNW9ducUabe5+6Hkz//ivFMpbetHIL9xj1ueNuWabWKo3Oh271HjurNf2/VNhfv+8\nYjvXWpmfk7+R2X5fmG17bVz55vHediN/j2V+Vr0sW3v6nnm8/XN95Wb/85uetR171cyu/R/zM/d+\nW9kvHLJ9zj5ulvXMB+b5fpxqto+fv2e+Fitmme2zu//4vv4js36rf+m/Ru2foS+lmp8F96rVyOdb\nyTcf3fyWeMb2BWlfX+4/zOzG779knt8yG8Hn9uumS3v9u5Gzt9ehutapruZnhfaE/t1h/41lXTTf\nN2+5bX/bn74uNfT3pfW57RqebbZv73GzDdl/HPwqxf/4xT/bznXezGpuvnnuH5rXzOr/tn2X/sj2\nXVrc+Zo3jx0o+nPOSMhgpKf3uq7MP6zW+ouqzD8sjXRLU2bcgqoBAAAAA4f7rCS5I1yL3uvLOSJ2\nIeeMxMbGqqHBPwSroaFBTqezy37T86fp3vhRmp4/jUAEAAAAd6QZTkma0SkNTP05ZyTkoofXr1/X\nQw89pPfee0/33XefHn30UZWVlRlzRmbMmKHKyso+qyAAAAAwEE2fPl1utzvS1QiLZVn6T++mLttL\nrBUBFz30eDxKTU2V0+nUu+++q/z8fG3btk3R0dGSpIKCAs2ePTtoeSGHaUVFRWnLli3KyMiQx+NR\nbm6uEYhIGnAvMAAAAIDgroXRE1JcXKykpCRdunRJ0o1gJi8vT3l5ed0ceUO3ix7OmTNHc+bM6XGF\nAAAAAAxcPZ0zcvbsWe3fv19r1qzRL37xC0mS1+sN2IMSTLeLHgIAAAC4c/R0zshLL72k1157TYMG\n+UMKy7JUUlKiyZMnKzc3V62trQGP/QeCEQAAAAA+bRqu/3PX68P8/b5kt2/fPo0ZM0YpKSlGT8iS\nJUtUV1enmpoajR07VsuXL+9ybGchJ7ADAAAAuHNYlqU53ne6bD9gPWkEHatXr9bbb7+tqKgoXb16\nVV988YWefPJJ7dixw7dPfX29MjMzdfz48aDl0TMCAAAAwKddQ7sku40bN6qhoUF1dXXauXOnvvOd\n72jHjh1qamry7VNeXq7k5OQux3bW7QR2AAAAAHeONg0Pa3+v1+tbLH3FihU6duyYLMvSuHHjtHXr\n1pDHMkwLAAAAgKQbw7T+1fvnLts/tB4L6y5ZPUXPCAAAAACfcNYZ6S2CEQAAAAA+V8IcptUbBCMA\nAAAAfK71cNHDW4FgBAAAAIBPsEUO+wLBCAAAAACfa+0EIwAAAAAi4MqXzBkBAAAAEAHtV/tvzggr\nsAMAAADwu3pX12Tf5epVTZ06VVOmTFFSUpJWrVolSWppaVFaWpomTJig9PR0tba2hiyKRQ8BAAAA\nSLqx6KGOBQgPJltdFj1sa2vT8OHDdf36dT322GPavHmz9u7dq2984xtasWKFioqKdOHCBRUWFgYt\nj54RAAAAAH5XAqQAhg+/Mbekvb1dHo9Ho0eP1t69e5WTkyNJysnJ0Z49e0IWRTACAAAAwO9agBRA\nR0eHpkyZIofDoZkzZ+rhhx9Wc3OzHA6HJMnhcKi5uTlkUUxgBwAAAOB3VdIxt/SRO+RugwYNUk1N\njS5evKiMjAz96U9/Mp63LOvGsK8QCEYAAAAA+H0p6cEZN9I//H590N1HjRql7373u/rwww/lcDj0\n2WefKSYmRk1NTRozZkzIohimBQAAAMDvaoBkc/78ed+dsq5cuaI//vGPSklJUVZWlkpLSyVJpaWl\nmjt3bsii6BkBAAAA4Bcg+LBrampSTk6OOjo61NHRoWeffVazZs1SSkqK5s+fr+3btys+Pl67d+8O\neR5u7QsAAABA0s1b+xYHCA9e7Hpr31uBnhEAAAAAfkFu5dsXCEYAAAAA+AW5lW9fIBgBAAAA4Hep\n/4oiGAEAAADgR88IAAAAgIhgzggAAACAiKBnBAAAAEBEfNl/RbECOwAAAAC/HqzAvmjRIjkcDiUn\nJ/u25efny+l0KiUlRSkpKXK5XN0WRTACAAAAwK8HwcjChQu7BBuWZSkvL0/V1dWqrq7W7Nmzuy2K\nYVoAAAAA/Hpwa99p06apvr6+y/ZwV2mnZwQAAACA37UAqYdKSko0efJk5ebmqrW1tdv9CUYAAAAA\n+F2R9LlbOpvvTz2wZMkS1dXVqaamRmPHjtXy5cu7PYZhWgAAAAD8rkkaPEO6Z4Z/24X13R42ZswY\n3+PFixcrMzOz22PoGQEAAADg92WA1ANNTU2+x+Xl5cadtoKhZwQAAACAX4C7Z9ktWLBAlZWVOn/+\nvOLi4rR+/Xq53W7V1NTIsiyNGzdOW7du7fY8ljfcKe8AAAAAvpIsy5KsAOGB1wr7Tlk9Qc8IAAAA\nAL9+7KogGAEAAADQyd/7rSSCEQAAAACdXOm3kghGAAAAAHTyRb+VRDACAAAAoBN6RgAAAABEBMEI\nAAAAgIggGAEAAAAQEZf6raRB/VYSAAAAgAHgSoDUlcvlUmJiosaPH6+ioqJ/qiSCEQAAAACdtAVI\nJo/Ho2XLlsnlcqm2tlZlZWU6ceJE2CURjAAAAADo5FKAZKqqqlJCQoLi4+M1ZMgQZWdnq6KiIuyS\nCEYAAAAAdNL9MK3GxkbFxcX58k6nU42NjWGXxAR2AAAAAJ1ckfSJpE+D7mFZ1i0piWAEAAAAQCcr\nu2wZPXq0kY+NjVVDQ4Mv39DQIKfTGXZJDNMCAAAAIEnyer0BU0tLi7FfamqqTp8+rfr6erW3t2vX\nrl3KysoKuzx6RgAAAACEJSoqSlu2bFFGRoY8Ho9yc3M1ceLEsM9jeb1ebx/UDwAAAABCYpgWAAAA\ngIggGAEAAAAQEQQjAAAAACKCYAQAAABARBCMAAAAAIgIghEAAAAAEUEwAgAAACAi/h8rt5vIXkip\nugAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xb98d350>"
       ]
      }
     ],
     "prompt_number": 119
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Do GPy ML based regression\n",
      "#Start with no SNPs, add SNPs by greedy forward selection - add best SNP at each step, is there a maxima?\n",
      "#Y ~ N(0, s^2 XX^T + e^2 I )\n",
      "\n",
      "X=X.transpose()\n",
      "Y=Y.transpose()\n",
      "#print X.shape\n",
      "#Y=Y.reshape((182,1))\n",
      "#print Y.shape\n",
      "chosenSNPS=[]\n",
      "X2=np.zeros((182,0))\n",
      "#print X2.shape\n",
      "#print X2.dtype\n",
      "k = GPy.kern.linear(input_dim=0)\n",
      "m = GPy.models.GPRegression(X2,Y,k)\n",
      "m.optimize(messages=False)\n",
      "#m['.*variance']=1\n",
      "#m['linear_variance']=10000\n",
      "bestll=m.log_likelihood()\n",
      "print m\n",
      "print bestll\n",
      "#provides initial parameters for 0 SNPS\n",
      "oldparams=m['.*variance']\n",
      "\n",
      "#print SNPsconsider\n",
      "#print np.array(set(SNPsconsider).difference(set(chosenSNPS)))\n",
      "for step in range(100):\n",
      "    hasbeenadded=False\n",
      "    #bestll=np.inf\n",
      "    bestSNP=np.nan\n",
      "    X2=np.hstack([X2, np.zeros((182, 1))])\n",
      "    tempbestll=bestll\n",
      "    for SNP in np.array(list(set(range(X.shape[1])).difference(set(chosenSNPS)))):\n",
      "    #for SNP in np.array(list(set(SNPsconsider).difference(set(chosenSNPS)))):\n",
      "        X2[:,-1]=X[:,SNP]\n",
      "        \n",
      "        k = GPy.kern.linear(input_dim=step+1)\n",
      "        tempm = GPy.models.GPRegression(X2,Y,k)\n",
      "        tempm.optimize(messages=False)\n",
      "        #m.X=X2\n",
      "        #m.kern.X=X2\n",
      "        #m.input_dim=step+1\n",
      "        #m.kern.input_dim=step+1\n",
      "        #m['.*variance']=1\n",
      "        #m['linear_variance']=10000\n",
      "        #m.optimize(messages=False)\n",
      "        #print m\n",
      "        #m.update_likelihood_approximation()\n",
      "        #print m.log_likelihood()\n",
      "        if tempm.log_likelihood() > tempbestll:\n",
      "            hasbeenadded=True\n",
      "            tempbestll=tempm.log_likelihood()\n",
      "            bestSNP=SNP\n",
      "            #print SNP\n",
      "    if hasbeenadded==False:\n",
      "        print \"converged: \" + str(len(chosenSNPS)) + \"SNPs\" \n",
      "        X2=X2[:,:-1]\n",
      "        break\n",
      "    X2[:,-1]=X[:,bestSNP]\n",
      "\n",
      "    #print \"X2\" + str(X2.shape)\n",
      "    k = GPy.kern.linear(input_dim=step+1)\n",
      "    tempm = GPy.models.GPRegression(X2,Y,k)\n",
      "    tempm.optimize(messages=False)\n",
      "    print tempm.log_likelihood()\n",
      "    if tempm.log_likelihood()>bestll:\n",
      "        bestll=tempm.log_likelihood()\n",
      "        m=tempm\n",
      "        oldparams=m['.*variance']\n",
      "        chosenSNPS.append(bestSNP)\n",
      "    else:\n",
      "        print \"converged: \" + str(len(chosenSNPS)) + \"SNPs\" \n",
      "        X2=X2[:,:-1]\n",
      "        break\n",
      "    #print X2.shape\n",
      "    \n",
      "pb.matshow(X2[:,:].transpose())\n",
      "pb.colorbar()\n",
      "pb.show()\n",
      "\n",
      "#m.plot()\n",
      "#print m\n",
      "\n",
      "#print m.predict(np.array([1]))\n",
      "#print m.predict(np.array([2]))\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Log-likelihood: -9.939e+02\n",
        "\n",
        "       Name        |    Value    |  Constraints  |  Ties  |  prior  \n",
        "--------------------------------------------------------------------\n",
        "  linear_variance  |   1.0000    |     (+ve)     |        |         \n",
        "  noise_variance   |  3046.3900  |     (+ve)     |        |         \n",
        "\n",
        "-993.939978346\n",
        "-862.503700329"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "-699.313601323"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "-269.822065657"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "converged: 3SNPs"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGWVJREFUeJzt3X1wVOX99/HP0V2Qp1ueI9nd373ExM1mgBjNQxmLoiMG\nmDGF0tHQVigGmqFkkGlnxGntVJx7xHTambak7UQHrQwawtSWYAuLk7RBxcIqgvprUINjdLNoplGQ\nJyVkm/sPYDe7e/bJZNkY3q+ZM3LO+e51HZNzknz32ut7GX19fX0CAAAAgAy6KtMXAAAAAAAkJgAA\nAAAyjsQEAAAAQMaRmAAAAADIOBITAAAAABlHYgIAAAAg40hMAAAAAKTk/vvvV1ZWlmbOnBkzZu3a\ntcrLy1NhYaEOHTqUsE0SEwAAAAApWbFihTweT8zzu3bt0tGjR9Xe3q4nnnhCq1evTtgmiQkAAACA\nlMyZM0cTJkyIeX7nzp1avny5JKmsrEwnTpxQV1dX3DZJTAAAAAAMKr/fL4fDEdy32+3q7OyM+xoS\nEwAAAACDrq+vL2zfMIy48SQmAAAAACRJowxDhsk2bty4lNqx2Wzy+XzB/c7OTtlstrivsXylKwYA\nAAAw7Hwp6f+ZHH/49OmU2qmoqFBdXZ0qKyu1f/9+jR8/XllZWXFfQ2ICAAAAIOj/JBGzdOlS7d27\nV93d3XI4HNqwYYPOnz8vSaqurtbChQu1a9cu5ebmasyYMXr66acTtmn0RX74CwAAAMAVyTAMPWly\nfJWi54wMNkZMAAAAAASNylC/JCYAAAAAgpL5KFc6UJULAAAAQNAok82Mx+NRfn6+8vLyVFtbG3X+\n+PHjWrx4sQoLC1VWVqZ///vfcfslMQEAAAAQlExiEggEVFNTI4/Ho7a2NjU0NOjIkSNhMY899phu\nuukmvfnmm9qyZYseeOCBuP2SmAAAAAAISiYx8Xq9ys3NldPplNVqVWVlpZqamsJijhw5ottvv12S\n5HK51NHRof/85z8x+yUxAQAAABA0zmSL5Pf75XA4gvt2u11+vz8sprCwUH/5y18kXUhkPvzwQ3V2\ndsbsl8QEAAAAQNBoky2SYRgJ23nooYd04sQJFRUVqa6uTkVFRbr66qtjxlOVCwAAAEDQKIv0yn+l\nff2XLYlYwsRms8nn8wX3fT6f7HZ7WMy4ceP01FNPBfenT5+unJycmP2ywCIAAAAASRdGQnqujT4+\n4vPwBRZ7e3vlcrnU0tKi7OxslZaWqqGhQW63Oxjz+eefa9SoURoxYoSefPJJ7du3T3/6059i9s2I\nCQAAAIAg68jEMRaLRXV1dSovL1cgEFBVVZXcbrfq6+slSdXV1Wpra9MPfvADGYahGTNmaPPmzXHb\nZMQEAAAAgKQLIyZ9/2Ny/KPwEZN0YMQEAAAAQEgSIybpQGICAAAAIGRMZrolMQEAAAAQck1mumUd\nEwAAAAAhI002Ex6PR/n5+crLy1NtbW3U+e7ubs2fP1833nijZsyYEbcil8TkdwAAAAAXGYahvoUm\nx3eFT34PBAJyuVxqbm6WzWZTSUlJVLngRx55ROfOndPGjRvV3d0tl8ulrq4uWSzmH9pixAQAAABA\nyDUmWwSv16vc3Fw5nU5ZrVZVVlaqqakpLGbatGk6efKkJOnkyZOaNGlSzKREYo4JAAAAgP6SqMrl\n9/vlcDiC+3a7XQcOHAiLWbVqle644w5lZ2fr1KlT2r59e9w2GTEBAAAAEDLGZItgGEbCZh577DHd\neOONOnbsmA4fPqw1a9bo1KlTMeMZMQEAAAAQco3U2im1+mOH2Gw2+Xy+4L7P55Pdbg+LefXVV/Wz\nn/1MknT99ddr+vTpevfdd1VcXGzaJokJAAAAgJCR0tzrL2yXbHgtPKS4uFjt7e3q6OhQdna2Ghsb\n1dDQEBaTn5+v5uZm3XLLLerq6tK7776rnJycmN2SmAAAAAAISWKOicViUV1dncrLyxUIBFRVVSW3\n2636+npJUnV1tX76059qxYoVKiws1H//+1/98pe/1MSJE2O2SblgAAAAAJIulgt+zOT4T8PLBacD\nIyYAAAAAQpIYMUkHEhMAAAAAIRlKTCgXDAAAACAkiXLBkuTxeJSfn6+8vDzV1tZGnf/Vr36loqIi\nFRUVaebMmbJYLDpx4kTMbpljAgAAAEDSxTkmDSbHl4bPMQkEAnK5XGpubpbNZlNJSYkaGhrkdrtN\n2/3b3/6m3/zmN2pubo7ZNyMmAAAAAEJGmmwRvF6vcnNz5XQ6ZbVaVVlZqaampphNPvfcc1q6dGnc\nbklMAAAAAIRcY7JF8Pv9cjgcwX273S6/33xFxrNnz2rPnj1asmRJ3G6Z/A4AAAAgJMackv4Mw0i6\nuRdeeEHf/OY3NX78+LhxJCYAAAAAQkZKrQel1jdih9hsNvl8vuC+z+eT3W43jd22bVvCj3FJTH4H\nAAAAcJFhGOo7bHL8xvDJ7729vXK5XGppaVF2drZKS0tNJ79//vnnysnJUWdnp0aNGhW370GZY5Ko\nVBgw3DmdTs2aNUtFRUUqLS2VJH322WeaN2+ebrjhBt11111xy+MBX2f333+/srKyNHPmzOCxePf/\nxo0blZeXp/z8fL344ouZuGQgbcyeh0ceeUR2uz1YNnX37t3BczwPGJKSKBdssVhUV1en8vJyFRQU\n6N5775Xb7VZ9fb3q6+uDcTt27FB5eXnCpEQahBGTVEuFAcPR9OnTdfDgQU2cODF47MEHH9TkyZP1\n4IMPqra2VsePH9fjjz+ewasE0uPll1/W2LFjtWzZMr399tuSYt//bW1t+u53v6vXXntNfr9fd955\np9577z1ddRW1WDA8mD0PGzZs0Lhx4/TjH/84LJbnAUORYRjq+8jk+P+Ej5ikw4Dv/FRLhQHDVeTD\nunPnTi1fvlyStHz5cu3YsSMTlwWk3Zw5czRhwoSwY7Hu/6amJi1dulRWq1VOp1O5ubnyer2X/ZqB\ndDF7HiTzP+h4HjBkJVGVKx0GnJikUioMGK4Mw9Cdd96p4uJiPfnkk5Kkrq4uZWVlSZKysrLU1dWV\nyUsELqtY9/+xY8fCJkfyOwNXik2bNqmwsFBVVVXBjzbyPGCoOjcmerscBpyYpFIqDBiu9u3bp0OH\nDmn37t36/e9/r5dffjnsvGEYPCu4YiW6/3k2MNytXr1aH3zwgQ4fPqxp06bpJz/5ScxYngcMBedG\njojazCQzz7y1tVVFRUWaMWOG5s6dG7ffAZcLTqVUGDBcTZs2TZI0ZcoULV68WF6vV1lZWfrkk090\n3XXX6eOPP9bUqVMzfJXA5RPr/o/8ndHZ2SmbzZapywQui/4//1euXKm7775bEs8Dhq6eq80SkZ6w\nvUAgoJqamrB55hUVFWHzzE+cOKE1a9Zoz549stvt6u7ujtvvgEdMiouL1d7ero6ODvX09KixsVEV\nFRUDbRb42jh79qxOnTolSTpz5oxefPFFzZw5UxUVFXrmmWckSc8884wWLVqUycsELqtY939FRYW2\nbdumnp4effDBB2pvbw9WsgOGq48//jj477/+9a/Bil08Dxiqzmlk1BYpmXnmzz33nJYsWRIctJg8\neXLcfgc8YtK/VFggEFBVVRUVuXBF6erq0uLFiyVdqOn9ve99T3fddZeKi4t1zz33aPPmzXI6ndq+\nfXuGrxRIj6VLl2rv3r3q7u6Ww+HQo48+qoceesj0/i8oKNA999yjgoICWSwW/eEPf+CjKxhWIp+H\nDRs2qLW1VYcPH5ZhGJo+fXqwlCrPA4aqL5S4tK/ZPPMDBw6ExbS3t+v8+fO6/fbbderUKT3wwAO6\n7777YrbJAosAAAAAJF2Y5/S/fddHHZ9hvB9WXe7555+Xx+MJFv3ZunWrDhw4oE2bNgVjampq9MYb\nb6ilpUVnz57V7Nmz9fe//115eXmmfQ94xAQAAADA8NGjEXq99Yxebz0bMyaZeeYOh0OTJ0/WqFGj\nNGrUKN1666168803YyYmjJgAAAAAkHRhxOSVvpujjn/TOBg2YtLb2yuXy6WWlhZlZ2ertLQ0apH1\nd955RzU1NdqzZ4/OnTunsrIyNTY2qqCgwLTvhJPfkykDBgAAAGB46NGIqC1S/3nmBQUFuvfee+V2\nu1VfXx+cR5Wfn6/58+dr1qxZKisr06pVq2ImJVKCEZNAICCXyxVWBiwyEwIAAAAwPBiGod19c6OO\nLzBale4PWsUdMUmmDBgAAACA4SOZEZN0iDv5PZkyYIbhlPRhGi4NAAAA+Dr7v+rr68j0RaTsrEZn\npN+4iYlhGHrppZeUlZWlqVOnav369SZRH0r6haRWSXOjzv5CG8L2N+gXcS8oMj6eyLZS7SvVvlNt\nL9X2B9L3QL/Ogx2fzrZS/Tqms/1Y92CrzJ6G1NtLZCD3VKL+B/v7lOh5jXctg9FfKq9NZKA/xwbz\n52CqbQ+0r6/afqsG/kyk+2f0QKV6zw/keRtIX8kYyr8jMvl8Rbb/1Z+fVn2VJyLdz2sqUv0+DDQ+\nnX8/JjKwr/vg/p6+XM4lOULi8Xi0bt06BQIBrVy5MipXaG1t1be+9S3l5ORIkpYsWaKHH344Zntx\nExObzabx48dr27ZtWrZsmWkZsIvdSuq4+F/nxQ0AAAC4cnRIuvD38Ndbj8lK75ECgYBqamrC5qJX\nVFREzUW/7bbbtHPnzqT6jZuYFBcXq7u7W2fOnFFfX58aGxvV0NAQFfcL7b2Y+38oaW/YuVQz5ngu\n97vJ6X7XJZW+Bmqw31lP5Z2Pgb6zN9BrT/VdmsF8pz3T714PpXcxExnofTGQ79tg35OJ4gf7Hd90\njioM3jvtrdr7FUbV+59P9es61EYdLtdrh0J/8b5vqfad7nfmE4l3n13u392X8x5N9wjGYP/MT+V6\nBnqPxLs2p6Tw0bDwv4u/Ls4msfJ7/7nokoJz0SMTk1QmzMed/H6pDNiyZcvU3t4eLANmxpl0l8Dw\n58z0BQBDjjPTFwAMIc5MXwAQV49GRm2RzOai+/3+sBjDMPTqq6+qsLBQCxcuVFtbW9x+E678PmPG\nDGVlZenYsWN69tlnNWbMGK1duzYqzpmoIeAK4sz0BQBDjjPTFwAMIc5MXwAQVzJzTAzDSBhz0003\nyefzafTo0dq9e7cWLVqk9957L2Z8wsTEarXq5z//udavX69//etfuvnmmzVv3rywkZPWfvFO8bgB\nAADgytMhaTjMMflCo/Vha4c+bI1deddms8nn8wX3zeaijxs3LvjvBQsW6Ec/+pE+++wzTZw40bTN\nhInJddddF1yhcezYsXK73Tp27FhYYjI3USMAAADAMOeUNBzmmJzTCF039wZdN/eG4LFXNrwUFlNc\nXKz29nZ1dHQoOzvbdC56V1eXpk6dKsMw5PV61dfXFzMpkZJITJYuXaq9e/fq008/1bRp03T+/Hlt\n3bo1Kq5DF74ZA5kwNNgyXTJ3IBPE48V+lWsZ7Mm4gzlhNVK6J30O5oS+WG11yHzkMN3ft8EsCZru\nvlO5FjOD+bwO9P8l3fGD/fp0tWUm8pnIZKnxVA3292kwi7+k2tblfl5TmWicSLqLSQyk/6/+de2Q\n5Ez7tQ+my122eSiVCx7o76Ovo2Sqcl2ai15eXq5AIKCqqiq53W7V19dLkqqrq/XnP/9Zf/zjH2Wx\nWDR69Ght27YtfpuJOr2U+Zw+fVpz587Vww8/rLFjx0bFdYiPcAGXdIjnAeivQzwTQEiHeCIwlCW7\njsmCBQu0YMGCsGPV1dXBf69Zs0Zr1qxJut+EiYkknT9/XkuWLNH3v/99LVq0KOp8q0KrmPCwAQAA\n4MrUoUszTb7OvhiKK79LF2oPV1VVqaCgQOvWrTONmavQGqZ7SUoAAABwRXIq/A36r+8ck0ww+hKs\nevLKK6/o1ltv1axZs4JlwTZu3Kj58+dLkubOnau9e7+eX3QAAAAgXW677Ta1trZm+jJSYhiGFvVF\nL6i+w1gatViix+PRunXrFAgEtHLlSq1fv960zddee02zZ8/W9u3b9e1vfzt234kSEwAAAABXBsMw\ndFdfU9TxF41vhSUmgUBALpdLzc3NstlsKikpUUNDQ9Ri7IFAQPPmzdPo0aO1YsUKLVmyJGbfcVd+\nBwAAAHBl6dGIqC2S1+tVbm6unE6nrFarKisr1dQUndBs2rRJ3/nOdzRlypSE/ZKYAAAAAAg6p5FR\nWyS/3y+HwxHct9vt8vv9UTFNTU1avXq1pMSrxSdVlQsAAADAlcFshCRSoiRDktatW6fHH39chmGo\nr68vao5KJBITAAAAAEFnNVpnWl/X2dbXY8bYbDb5fL7gvs/nk91uD4s5ePCgKisrJUnd3d3avXu3\nrFarKioqTNtk8jsAAAAASRdGQq7v+9+o4+8bM8JGPHp7e+VyudTS0qLs7GyVlpaaTn6/ZMWKFbr7\n7rvjVuVixAQAAABAkNmckkgWi0V1dXUqLy9XIBBQVVWV3G636uvrJYWvAJ8sRkwAAAAASLowYjIp\n0Bl1/NOr7QnniAwUIyYAAAAAgs59mXjEJB1ITAAAAAAE9XyZuCpXOrCOCQAAAICgntOjozYzHo9H\n+fn5ysvLU21tbdT5pqYmFRYWqqioSDfffLP+8Y9/xO2XOSYAAAAAJF1cn+Rdk/TAZYTNMQkEAnK5\nXGpubpbNZlNJSUlUVa4zZ85ozJgxkqS3335bixcv1tGjR2P2zYgJAAAAgJAvTbYIXq9Xubm5cjqd\nslqtqqysVFNTU1jMpaREkk6fPq3JkyfH7ZbEBAAAAEDIFyZbBL/fL4fDEdy32+3y+/1RcTt27JDb\n7daCBQv0u9/9Lm63JCYAAAAAQs6YbBEMw0iqqUWLFunIkSN64YUXdN9998WNpSoXAAAAgJAvJb3V\nKr3dGjPEZrPJ5/MF930+n+x2e8z4OXPmqLe3V59++qkmTZpkGkNiAgAAACDkS0k3zL2wXfLchrCQ\n4uJitbe3q6OjQ9nZ2WpsbFRDQ0NYzPvvv6+cnBwZhqE33nhDkmImJRKJCQAAAID+TicOsVgsqqur\nU3l5uQKBgKqqquR2u1VfXy9Jqq6u1vPPP68tW7bIarVq7Nix2rZtW9w2KRcMAAAAQNLFuSP1JulB\ndXi54HRgxAQAAABAiEl54MuBxAQAAABAiEl54MuBcsEAAAAAQpIoFyxJHo9H+fn5ysvLU21tbdT5\nZ599VoWFhZo1a5ZuueUWvfXWW3G7ZcQEAAAAQEgSIyaBQEA1NTVqbm6WzWZTSUmJKioq5Ha7gzE5\nOTl66aWXdO2118rj8eiHP/yh9u/fH7NNRkwAAAAAhJwz2SJ4vV7l5ubK6XTKarWqsrJSTU1NYTGz\nZ8/WtddeK0kqKytTZ2dn3G5JTAAAAACEnDbZIvj9fjkcjuC+3W6X3++P2eTmzZu1cOHCuN3yUS4A\nAAAAIUlU5TIMI+nm/vnPf+qpp57Svn374saRmAAAAAAI+VLSf1ql7taYITabTT6fL7jv8/lkt9uj\n4t566y2tWrVKHo9HEyZMiNstCywCAAAAkHRxJGSBSXqwO3yBxd7eXrlcLrW0tCg7O1ulpaVqaGgI\nm/z+0Ucf6Y477tDWrVv1jW98I2HfjJgAAAAACDGZUxLJYrGorq5O5eXlCgQCqqqqktvtVn19vSSp\nurpajz76qI4fP67Vq1dLkqxWq7xeb8w2GTEBAAAAIOniiEmxSXrweviISTowYgIAAAAgxKQ88OVA\nYgIAAAAgJImPcqUDiQkAAACAkCTKBacDCywCAAAACPnSZDPh8XiUn5+vvLw81dbWRp1/5513NHv2\nbF1zzTX69a9/nbBbRkwAAAAAhJxKHBIIBFRTU6Pm5mbZbDaVlJSooqIirFzwpEmTtGnTJu3YsSOp\nbhkxAQAAABDSa7JF8Hq9ys3NldPplNVqVWVlpZqamsJipkyZouLiYlmt1qS6JTEBAAAAkBK/3y+H\nwxHct9vt8vv9A2qTj3IBAAAA6Od8wgjDMAa9VxITAAAAAP2ckvSKpH0xI2w2m3w+X3Df5/PJbrcP\nqFcSEwAAAAD9nJV008Xtkl+GRRQXF6u9vV0dHR3Kzs5WY2OjGhoaTFtLdsV4EhMAAAAA/XyRMMJi\nsaiurk7l5eUKBAKqqqqS2+1WfX29JKm6ulqffPKJSkpKdPLkSV111VX67W9/q7a2No0dO9a0TaMv\n2RQGAAAAwLB2Ye7IQZMzNyc98vFVMWICAAAAoJ/EIybpQGICAAAAoJ/MJCasYwIAAACgny9Mtmge\nj0f5+fnKy8tTbW2taczatWuVl5enwsJCHTp0KG6vJCYAAAAA+jlpsoULBAKqqamRx+NRW1ubGhoa\ndOTIkbCYXbt26ejRo2pvb9cTTzyh1atXx+2VxAQAAABAP4lHTLxer3Jzc+V0OmW1WlVZWammpqaw\nmJ07d2r58uWSpLKyMp04cUJdXV0xeyUxAQAAANBP4sTE7/fL4XAE9+12u/x+f8KYzs7OmL0y+R0A\nAABAP9Ef3Yp0oaxwYpElhuO9jsQEAAAAQD8PRx2JXBTRZrPJ5/MF930+n+x2e9yYzs5O2Wy2mL3y\nUS4AAAAAki6McJhtp06dCosrLi5We3u7Ojo61NPTo8bGRlVUVITFVFRUaMuWLZKk/fv3a/z48crK\nyorZNyMmAAAAAFJisVhUV1en8vJyBQIBVVVVye12q76+XpJUXV2thQsXateuXcrNzdWYMWP09NNP\nx23T6Ev32vIAAAAAkAAf5QIAAACQcSQmAAAAADKOxAQAAABAxpGYAAAAAMg4EhMAAAAAGUdiAgAA\nACDjSEwAAAAAZByJCQAAAICM+/+7cvzcC3B2TwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xc825dd0>"
       ]
      }
     ],
     "prompt_number": 120
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print chosenSNPS\n",
      "print m\n",
      "print m.kern\n",
      "m.plot()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[3, 5, 1]\n",
        "\n",
        "Log-likelihood: -2.698e+02\n",
        "\n",
        "       Name        |   Value    |  Constraints  |  Ties  |  prior  \n",
        "-------------------------------------------------------------------\n",
        "  linear_variance  |  615.4866  |     (+ve)     |        |         \n",
        "  noise_variance   |   0.9353   |     (+ve)     |        |         \n",
        "\n",
        "       Name        |   Value    |  Constraints  |  Ties  \n",
        "---------------------------------------------------------\n",
        "  linear_variance  |  615.4866  |               |        \n",
        "\n"
       ]
      },
      {
       "ename": "NotImplementedError",
       "evalue": "Cannot define a frame with more than two input dimensions",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-121-acb802d768d1>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mm\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mkern\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m \u001b[0mm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[1;32m/is/ei/jmcmurr/virtualenv/local/lib/python2.7/site-packages/GPy-0.4.9-py2.7.egg/GPy/core/gp_base.pyc\u001b[0m in \u001b[0;36mplot\u001b[1;34m(self, plot_limits, which_data_rows, which_data_ycols, which_parts, fixed_inputs, levels, samples, fignum, ax, resolution, plot_raw, linecol, fillcol)\u001b[0m\n\u001b[0;32m    230\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0mcontour\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mscatter\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    231\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 232\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"Cannot define a frame with more than two input dimensions\"\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    233\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    234\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mgetstate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mNotImplementedError\u001b[0m: Cannot define a frame with more than two input dimensions"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEIxJREFUeJzt3W9Ilff/x/HXiXNuRET/HJLnHLA8Bz1iamBZRHHaCF1j\n3qhu2K1mIiJEtFsbdSPtRsvujbxj0B/WShpbYDA7QdFpUDmjoqBamNR2PDCZK3EsKDtdvxv7fXVm\nHY96PNZ7zwcIXpyP1/XeB3nu6vwpl+M4jgAApsya6QEAAOlH3AHAIOIOAAYRdwAwiLgDgEHEHQAM\nGjfu27dvV3Z2tpYtW/bWNTt37lQwGFRJSYlu3bqV1gEBABM3btxramoUiUTe+nhHR4cePnyo7u5u\nHT58WA0NDWkdEAAwcePGfe3atVqwYMFbHz979qy2bdsmSSovL9fAwID6+vrSNyEAYMKm/Jx7PB6X\n3+8fPvb5fOrt7Z3qaQEAU5CWF1Rf/xsMXC5XOk4LAJgk91RP4PV6FYvFho97e3vl9XrHrAsEAurp\n6Znq5QDgPyUvL08PHz6c8M9N+c69qqpK33zzjSSps7NT8+fPV3Z29ph1PT09chyHL8fR3r17Z3yG\nd+WLvWAv2IvkX5O9KR73zn3r1q26fPmy+vv75ff71dTUpKGhIUlSfX29Nm7cqI6ODgUCAc2ZM0fH\njh2b1CAAgPQZN+5tbW3jnqSlpSUtwwAA0oNPqM6AcDg80yO8M9iLEezFCPZi6lyO42TkH+twuVzK\n0KUAwIzJtpM7dwAwiLgDgEHEHQAMIu4AYBBxBwCDiDsAGETcAcAg4g4ABhF3ADCIuAOAQcQdAAwi\n7gBgEHEHAIOIOwAYRNwBwCDiDgAGEXcAMIi4A4BBxB0ADCLuAGAQcQcAg4g7ABhE3AHAIOIOAAYR\ndwAwiLgDgEHEHQAMIu4AYBBxBwCDiDsAGETcAcAg4g4ABhF3ADCIuAOAQePGPRKJqKCgQMFgUM3N\nzWMe7+/vV2VlpUpLS1VUVKTjx49Px5wAgAlwOY7jvO3BRCKh/Px8XbhwQV6vVytWrFBbW5tCodDw\nmsbGRj1//lxfffWV+vv7lZ+fr76+Prnd7tEXcrmU5FIAgDeYbDuT3rl3dXUpEAgoNzdXHo9H1dXV\nam9vH7Vm8eLFGhwclCQNDg5q0aJFY8IOAMispBWOx+Py+/3Dxz6fTz///POoNXV1dfrwww+Vk5Oj\nv/76S9999930TAoASFnSuLtcrnFPsH//fpWWlioajaqnp0cbNmzQ7du3NXfu3DFrGxsbh78Ph8MK\nh8MTHhgALItGo4pGo1M+T9K4e71exWKx4eNYLCafzzdqzdWrV7Vnzx5JUl5enpYsWaIHDx6orKxs\nzPn+HXcAwFiv3/g2NTVN6jxJn3MvKytTd3e3Hj9+rBcvXuj06dOqqqoataagoEAXLlyQJPX19enB\ngwdaunTppIYBAKRH0jt3t9utlpYWVVRUKJFIqLa2VqFQSK2trZKk+vp67d69WzU1NSopKdGrV690\n8OBBLVy4MCPDAwDeLOlbIdN6Id4KCQATNi1vhQQAvJ+IOwAYRNwBwCDiDgAGEXcAMIi4A4BBxB0A\nDCLuAGAQcQcAg4g7ABhE3AHAIOIOAAYRdwAwiLgDgEHEHQAMIu4AYBBxBwCDiDsAGETcAcAg4g4A\nBhF3ADCIuAOAQcQdAAwi7gBgEHEHAIOIOwAYRNwBwCDiDgAGEXcAMIi4A4BBxB0ADCLuAGAQcQcA\ng4g7ABhE3AHAoHHjHolEVFBQoGAwqObm5jeuiUajWr58uYqKihQOh9M9IwBgglyO4zhvezCRSCg/\nP18XLlyQ1+vVihUr1NbWplAoNLxmYGBAa9as0fnz5+Xz+dTf36+srKyxF3K5lORSAIA3mGw7k965\nd3V1KRAIKDc3Vx6PR9XV1Wpvbx+15tSpU9q8ebN8Pp8kvTHsAIDMShr3eDwuv98/fOzz+RSPx0et\n6e7u1pMnT7R+/XqVlZXpxIkT0zMpACBl7mQPulyucU8wNDSkmzdv6uLFi3r27JlWr16tVatWKRgM\npm1IAMDEJI271+tVLBYbPo7FYsNPv/yP3+9XVlaWZs+erdmzZ2vdunW6ffv2G+Pe2Ng4/H04HObF\nVwB4TTQaVTQanfJ5kr6g+vLlS+Xn5+vixYvKycnRypUrx7yg+ssvv2jHjh06f/68nj9/rvLycp0+\nfVqFhYWjL8QLqgAwYZNtZ9I7d7fbrZaWFlVUVCiRSKi2tlahUEitra2SpPr6ehUUFKiyslLFxcWa\nNWuW6urqxoQdAJBZSe/c03oh7twBYMKm5a2QAID3E3EHAIOIOwAYRNwBwCDiDgAGEXcAMIi4A4BB\nxB0ADCLuAGAQcQcAg4g7ABhE3AHAIOIOAAYRdwAwiLgDgEHEHQAMIu4AYBBxBwCDiDsAGETcAcAg\n4g4ABhF3ADCIuAOAQcQdAAwi7gBgEHEHAIOIOwAYRNwBwCDiDgAGEXcAMIi4A4BBxB0ADCLuAGAQ\ncQcAg4g7ABhE3AHAoHHjHolEVFBQoGAwqObm5reuu379utxut86cOZPWAQEAE5c07olEQjt27FAk\nEtG9e/fU1tam+/fvv3HdF198ocrKSjmOM23DAgBSkzTuXV1dCgQCys3NlcfjUXV1tdrb28esO3To\nkLZs2aIPPvhg2gYFAKQuadzj8bj8fv/wsc/nUzweH7Omvb1dDQ0NkiSXyzUNYwIAJiJp3FMJ9a5d\nu3TgwAG5XC45jsPTMgDwDnAne9Dr9SoWiw0fx2Ix+Xy+UWtu3Lih6upqSVJ/f7/OnTsnj8ejqqqq\nMedrbGwc/j4cDiscDk9hdACwJxqNKhqNTvk8LifJrfbLly+Vn5+vixcvKicnRytXrlRbW5tCodAb\n19fU1OjTTz/Vpk2bxl7o/+/sAQCpm2w7k965u91utbS0qKKiQolEQrW1tQqFQmptbZUk1dfXT25a\nAMC0SnrnntYLcecOABM22XbyCVUAMIi4A4BBxB0ADCLuAGAQcQcAg4g7ABhE3AHAIOIOAAYRdwAw\niLgDgEHEHQAMIu4AYBBxBwCDiDsAGETcAcAg4g4ABhF3ADCIuAOAQcQdAAwi7gBgEHEHAIOIOwAY\nRNwBwCDiDgAGEXcAMIi4A4BBxB0ADCLuAGAQcQcAg4g7ABhE3AHAIOIOAAYRdwAwiLgDgEHEHQAM\nIu4AYFBKcY9EIiooKFAwGFRzc/OYx0+ePKmSkhIVFxdrzZo1unPnTtoHBQCkzuU4jpNsQSKRUH5+\nvi5cuCCv16sVK1aora1NoVBoeM21a9dUWFioefPmKRKJqLGxUZ2dnaMv5HJpnEsBAF4z2XaOe+fe\n1dWlQCCg3NxceTweVVdXq729fdSa1atXa968eZKk8vJy9fb2TngQAED6jBv3eDwuv98/fOzz+RSP\nx9+6/siRI9q4cWN6pgMATIp7vAUulyvlk126dElHjx7VlStX3vh4Y2Pj8PfhcFjhcDjlcwPAf0E0\nGlU0Gp3yecaNu9frVSwWGz6OxWLy+Xxj1t25c0d1dXWKRCJasGDBG8/177gDAMZ6/ca3qalpUucZ\n92mZsrIydXd36/Hjx3rx4oVOnz6tqqqqUWt+++03bdq0Sd9++60CgcCkBgEApM+4d+5ut1stLS2q\nqKhQIpFQbW2tQqGQWltbJUn19fXat2+fnj59qoaGBkmSx+NRV1fX9E4OAHircd8KmbYL8VZIAJiw\naXsrJADg/UPcAcAg4g4ABhF3ADCIuAOAQcQdAAwi7gBgEHEHAIOIOwAYRNwBwCDiDgAGEXcAMIi4\nA4BBxB0ADCLuAGAQcQcAg4g7ABhE3AHAIOIOAAYRdwAwiLgDgEHEHQAMIu4AYBBxBwCDiDsAGETc\nAcAg4g4ABhF3ADCIuAOAQcQdAAwi7gBgEHEHAIOIOwAYRNwBwCDiDgAGjRv3SCSigoICBYNBNTc3\nv3HNzp07FQwGVVJSolu3bqV9SADAxCSNeyKR0I4dOxSJRHTv3j21tbXp/v37o9Z0dHTo4cOH6u7u\n1uHDh9XQ0DCtA1sQjUZneoR3Bnsxgr0YwV5MXdK4d3V1KRAIKDc3Vx6PR9XV1Wpvbx+15uzZs9q2\nbZskqby8XAMDA+rr65u+iQ3gF3cEezGCvRjBXkxd0rjH43H5/f7hY5/Pp3g8Pu6a3t7eNI8JAJiI\npHF3uVwpncRxnEn9HABgeriTPej1ehWLxYaPY7GYfD5f0jW9vb3yer1jzpWXl0f0/6WpqWmmR3hn\nsBcj2IsR7MU/8vLyJvVzSeNeVlam7u5uPX78WDk5OTp9+rTa2tpGramqqlJLS4uqq6vV2dmp+fPn\nKzs7e8y5Hj58OKkBAQATlzTubrdbLS0tqqioUCKRUG1trUKhkFpbWyVJ9fX12rhxozo6OhQIBDRn\nzhwdO3YsI4MDAN7O5bz+hDkA4L2X9k+o8qGnEePtxcmTJ1VSUqLi4mKtWbNGd+7cmYEpMyOV3wtJ\nun79utxut86cOZPB6TInlX2IRqNavny5ioqKFA6HMztgBo23F/39/aqsrFRpaamKiop0/PjxzA+Z\nIdu3b1d2draWLVv21jUT7qaTRi9fvnTy8vKcR48eOS9evHBKSkqce/fujVrz448/Oh9//LHjOI7T\n2dnplJeXp3OEd0Yqe3H16lVnYGDAcRzHOXfu3H96L/63bv369c4nn3zifP/99zMw6fRKZR+ePn3q\nFBYWOrFYzHEcx/njjz9mYtRpl8pe7N271/nyyy8dx/lnHxYuXOgMDQ3NxLjT7qeffnJu3rzpFBUV\nvfHxyXQzrXfufOhpRCp7sXr1as2bN0/SP3th9fMBqeyFJB06dEhbtmzRBx98MANTTr9U9uHUqVPa\nvHnz8LvSsrKyZmLUaZfKXixevFiDg4OSpMHBQS1atEhud9KXCd9ba9eu1YIFC976+GS6mda486Gn\nEansxb8dOXJEGzduzMRoGZfq70V7e/vwX19h8W2zqexDd3e3njx5ovXr16usrEwnTpzI9JgZkcpe\n1NXV6e7du8rJyVFJSYm+/vrrTI/5zphMN9P6v0E+9DRiIv9Nly5d0tGjR3XlypVpnGjmpLIXu3bt\n0oEDB+RyueQ4zpjfEQtS2YehoSHdvHlTFy9e1LNnz7R69WqtWrVKwWAwAxNmTip7sX//fpWWlioa\njaqnp0cbNmzQ7du3NXfu3AxM+O6ZaDfTGvd0fujpfZfKXkjSnTt3VFdXp0gkkvSPZe+zVPbixo0b\nqq6ulvTPC2nnzp2Tx+NRVVVVRmedTqnsg9/vV1ZWlmbPnq3Zs2dr3bp1un37trm4p7IXV69e1Z49\neyT980GeJUuW6MGDByorK8vorO+CSXUzba8IOI4zNDTkLF261Hn06JHz/PnzcV9QvXbtmtkXEVPZ\ni19//dXJy8tzrl27NkNTZkYqe/Fvn332mfPDDz9kcMLMSGUf7t+/73z00UfOy5cvnb///tspKipy\n7t69O0MTT59U9uLzzz93GhsbHcdxnN9//93xer3On3/+ORPjZsSjR49SekE11W6m9c6dDz2NSGUv\n9u3bp6dPnw4/z+zxeNTV1TWTY0+LVPbivyCVfSgoKFBlZaWKi4s1a9Ys1dXVqbCwcIYnT79U9mL3\n7t2qqalRSUmJXr16pYMHD2rhwoUzPPn02Lp1qy5fvqz+/n75/X41NTVpaGhI0uS7yYeYAMAg/pk9\nADCIuAOAQcQdAAwi7gBgEHEHAIOIOwAYRNwBwCDiDgAG/R9h/b4XlfACjQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x39bfcd0>"
       ]
      }
     ],
     "prompt_number": 121
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "test=np.zeros((1,66))\n",
      "test[0,1]=1\n",
      "print test.shape\n",
      "print m.input_dim\n",
      "#m.predict(test.transpose())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(1, 66)\n",
        "1\n"
       ]
      }
     ],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## Compute the biased estimator of HSIC.\n",
      "# @param x The data.\n",
      "# @param y The labels.\n",
      "# @param kernelx The kernel on the data, default to linear kernel.\n",
      "# @param kernely The kernel on the labels, default to linear kernel.\n",
      "#\n",
      "\n",
      "def BiasedHSIC(x, y, kernelx, kernely):\n",
      "    #nx = kernelx.input_dim\n",
      "    #ny = kernely.input_dim\n",
      "    #print x.shape\n",
      "    nx=float(x.shape[0])\n",
      "    ny=float(y.shape[0])\n",
      "    \n",
      "    assert nx == ny, \\\n",
      "           \"Argument 1 and 2 have different number of data points\"\n",
      "\n",
      "    kMat=kernelx.K(x)\n",
      "    #if len(nx) > 1:\n",
      "    #    kMat = kernelx.Dot(x, x)\n",
      "    #else:\n",
      "    #    kMat = numpy.outerproduct(x, x)\n",
      "\n",
      "    #print kMat\n",
      "    hlhMat = ComputeHLH(y, kernely)\n",
      "    #print hlhMat\n",
      "    return numpy.sum(numpy.sum(kMat * hlhMat)) / ((nx-1)*(nx-1))\n",
      "        \n",
      "## Compute HLH give the labels.\n",
      "# @param y The labels.\n",
      "# @param kernely The kernel on the labels, default to linear kernel.\n",
      "#\n",
      "def ComputeHLH(y, kernely):\n",
      "    #ny = kernely.input_dim\n",
      "    ny=float(y.shape[0])\n",
      "    #if len(ny) > 1:\n",
      "    #    lMat = kernely.Dot(y, y)\n",
      "    #else:\n",
      "    #    lMat = numpy.outerproduct(y, y)\n",
      "    #print y.shape\n",
      "    lMat=kernely.K(y)\n",
      "    #print lMat\n",
      "    sL = numpy.sum(lMat, axis=1)\n",
      "    ssL = numpy.sum(sL)\n",
      "    # hlhMat\n",
      "    return lMat - numpy.add.outer(sL, sL)/ny + ssL/(ny*ny)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 113
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Do GPy ML based regression\n",
      "#Start with no SNPs, add SNPs by greedy forward selection - add best SNP at each step, is there a maxima?\n",
      "#Y ~ N(0, s^2 XX^T + e^2 I )\n",
      "\n",
      "#X=X.transpose()\n",
      "#Y=Y.transpose()\n",
      "#print X.shape\n",
      "#print Y.shape\n",
      "chosenSNPS2=[]\n",
      "X2=np.zeros((182,0))\n",
      "#print X2.shape\n",
      "#print X2.dtype\n",
      "\n",
      "Y=Y.transpose()\n",
      "#X=X.transpose()\n",
      "print X2.shape\n",
      "print Y.shape\n",
      "k2 = GPy.kern.linear(input_dim=0)\n",
      "m2 = GPy.models.GPRegression(X2,Y,k2)\n",
      "m2.optimize(messages=False)\n",
      "#m['.*variance']=1\n",
      "#m['linear_variance']=10000\n",
      "bestll=m2.log_likelihood()\n",
      "X=X.transpose()\n",
      "oldHSIC=0\n",
      "#print m\n",
      "#print bestll\n",
      "#provides initial parameters for 0 SNPS\n",
      "oldparams=m['.*variance']\n",
      "kernelx2 = GPy.kern.linear(input_dim=ng)\n",
      "kernely2 = GPy.kern.linear(input_dim=X.shape[1])\n",
      "\n",
      "#print X.shape\n",
      "for step in range(500):\n",
      "    hasbeenadded=False\n",
      "    #bestll=np.inf\n",
      "    bestSNP=np.nan\n",
      "    X2=np.hstack([X2, np.zeros((182, 1))])\n",
      "    tempbestHSIC=oldHSIC\n",
      "\n",
      "    for SNP in np.array(list(set(range(X.shape[1])).difference(set(chosenSNPS2)))):\n",
      "    #for SNP in np.array(list(set(SNPsconsider).difference(set(chosenSNPS2)))):\n",
      "\n",
      "        X2[:,-1]=X[:,SNP]\n",
      "        #m2.X=X2\n",
      "        #m2.kern.X=X2\n",
      "        #m2.input_dim=step+1\n",
      "        #m2.kern.input_dim=step+1\n",
      "\n",
      "        #m2.update_likelihood_approximation()\n",
      " \n",
      "        #m2._Xoffset = np.zeros((1, m2.input_dim))\n",
      "        #m2._Xscale = np.ones((1, m2.input_dim))\n",
      "        k2 = GPy.kern.linear(input_dim=step+1)\n",
      "\n",
      "        tempm2 = GPy.models.GPRegression(X2,Y,k2)\n",
      "        tempm2.optimize(messages=False)\n",
      "        residuals=tempm2.predict(X2)[0]\n",
      "\n",
      "        #residuals=m2.predict(X2)[0]\n",
      "\n",
      "        newHSIC=BiasedHSIC(residuals, X, kernelx2, kernely2)\n",
      "        if newHSIC > tempbestHSIC:\n",
      "            hasbeenadded=True\n",
      "            #oldHSIC=newHSIC\n",
      "            tempbestHSIC=newHSIC\n",
      "            bestSNP=SNP\n",
      "            \n",
      "    if hasbeenadded==False:\n",
      "        print \"converged: \" + str(len(chosenSNPS2)) + \"SNPs\" \n",
      "        X2=X2[:,:-1]\n",
      "        break\n",
      "    X2[:,-1]=X[:,bestSNP]\n",
      "\n",
      "    #print \"X2\" + str(X2.shape)\n",
      "\n",
      "    k2 = GPy.kern.linear(input_dim=step+1)\n",
      "\n",
      "    tempm2 = GPy.models.GPRegression(X2,Y,k2)\n",
      "    tempm2.optimize(messages=False)\n",
      "    residuals=tempm2.predict(X2)[0]\n",
      "\n",
      "    newHSIC=BiasedHSIC(residuals, X, kernelx2, kernely2)\n",
      "    if newHSIC > oldHSIC:\n",
      "        #hasbeenadded=True\n",
      "        oldHSIC=newHSIC\n",
      "        m2=tempm2\n",
      "        print newHSIC\n",
      "        print tempm2.log_likelihood()\n",
      "        chosenSNPS2.append(bestSNP)\n",
      "    else:\n",
      "        print \"HSIC not improved\"\n",
      "        print \"converged: \" + str(len(chosenSNPS2)) + \"SNPs\" \n",
      "        X2=X2[:,:-1]\n",
      "        break\n",
      "    #print X2.shape\n",
      "    \n",
      "pb.matshow(X2[:,:].transpose())\n",
      "pb.colorbar()\n",
      "pb.show()\n",
      "\n",
      "#m.plot()\n",
      "#print m\n",
      "\n",
      "#print m.predict(np.array([1]))\n",
      "#print m.predict(np.array([2]))\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(182, 0)\n",
        "(182, 1)\n",
        "371.728645593"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "-862.503700329\n",
        "converged: 1SNPs"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFn9JREFUeJzt3X9sVfX9x/HXwXvB8uPLbyu9924XvNfeS4B6tQUNY6Jx\nFki8g7mYsmUidnjD1q8zW+KcWzJxRq1xixt1SzH4g22WmsxxcWkvpmxFxcBVBPU7mBZj9faqjRUq\nPyrUXvv9A7y3995zf9RSLivPR/IJnHPf/XxOyjkt77zP5/Mx+vv7+wUAAAAABTSq0BcAAAAAACQm\nAAAAAAqOxAQAAABAwZGYAAAAACg4EhMAAAAABUdiAgAAAKDgSEwAAAAADMqtt96q4uJizZ07N2PM\n7bffLrfbrbKyMu3duzdnnyQmAAAAAAZl9erVCoVCGT9vamrSwYMH1dbWpg0bNmjt2rU5+yQxAQAA\nADAoixYt0uTJkzN+vnXrVq1atUqStGDBAnV3d6uzszNrnyQmAAAAAM6oaDQqh8MRP7bb7ero6Mj6\nNSQmAAAAAM64/v7+pGPDMLLGk5gAAAAAkCQVGYYMkzZhwoRB9WOz2RSJROLHHR0dstlsWb/G8pWu\nGAAAAMCIc0LSfSbnf3Xs2KD68fv9qqurU1VVlXbt2qVJkyapuLg469eQmAAAAACI+588YlauXKkd\nO3aoq6tLDodD69at0+effy5JCgQCWrZsmZqamuRyuTRu3Dg98cQTOfs0+lNf/gIAAABwXjIMQ4+Z\nnF+j9DkjZxoVEwAAAABxRQUal8QEAAAAQFw+r3INB1blAgAAABBXZNLMhEIheTweud1u1dbWpn1+\n+PBhrVixQmVlZVqwYIH+/e9/Zx2XxAQAAABAXD6JSSwWU01NjUKhkPbv36+GhgYdOHAgKeb+++/X\n5Zdfrtdff12bNm3ST37yk6zjkpgAAAAAiMsnMQmHw3K5XHI6nbJaraqqqlIwGEyKOXDggK655hpJ\nUmlpqdrb2/Xxxx9nHJfEBAAAAEDcBJOWKhqNyuFwxI/tdrui0WhSTFlZmZ599llJpxKZ9957Tx0d\nHRnHJTEBAAAAEDfWpKUyDCNnP3fddZe6u7vl8/lUV1cnn8+nCy64IGM8q3IBAAAAiCuySC99Ie0c\nuG1JyhYmNptNkUgkfhyJRGS325NiJkyYoMcffzx+PHPmTM2aNSvjuGywCAAAAEDSqUpI78T086M/\nTd5gsa+vT6Wlpdq+fbtKSko0f/58NTQ0yOv1xmM+/fRTFRUVafTo0Xrssce0c+dOPfnkkxnHpmIC\nAAAAIM46JneMxWJRXV2dKisrFYvFVF1dLa/Xq/r6eklSIBDQ/v37dcstt8gwDM2ZM0cbN27M2icV\nEwAAAACSTlVM+r9mcv795IrJcKBiAgAAACAhj4rJcCAxAQAAAJAwrjDDkpgAAAAASLiwMMOyjwkA\nAACAhDEmzUQoFJLH45Hb7VZtbW3a511dXVqyZIkuu+wyzZkzJ+uKXBKT3wEAAACcZhiG+peZnG9K\nnvwei8VUWlqqlpYW2Ww2VVRUpC0XfM899+jkyZN64IEH1NXVpdLSUnV2dspiMX9pi4oJAAAAgIQL\nTVqKcDgsl8slp9Mpq9WqqqoqBYPBpJgZM2boyJEjkqQjR45o6tSpGZMSiTkmAAAAAAbKY1WuaDQq\nh8MRP7bb7dq9e3dSzJo1a3TttdeqpKRER48e1TPPPJO1TyomAAAAABLGmbQUhmHk7Ob+++/XZZdd\npg8++ED79u3Tj3/8Yx09ejRjPBUTAAAAAAkXSq0dUms0c4jNZlMkEokfRyIR2e32pJiXX35Zv/zl\nLyVJl1xyiWbOnKm33npL5eXlpn2SmAAAAABIGCMtvuRU+9K6V5JDysvL1dbWpvb2dpWUlKixsVEN\nDQ1JMR6PRy0tLVq4cKE6Ozv11ltvadasWRmHJTEBAAAAkJDHHBOLxaK6ujpVVlYqFoupurpaXq9X\n9fX1kqRAIKC7775bq1evVllZmb744gs99NBDmjJlSsY+WS4YAAAAgKTTywXfb3L+7uTlgocDFRMA\nAAAACXlUTIYDiQkAAACAhAIlJiwXDAAAACAhj+WCJSkUCsnj8cjtdqu2tjbt84cfflg+n08+n09z\n586VxWJRd3d3xmGZYwIAAABA0uk5Jg0m51cmzzGJxWIqLS1VS0uLbDabKioq1NDQIK/Xa9rvP/7x\nDz3yyCNqaWnJODYVEwAAAAAJY0xainA4LJfLJafTKavVqqqqKgWDwYxdPv3001q5cmXWYUlMAAAA\nACRcaNJSRKNRORyO+LHdblc0ar4jY09Pj7Zt26Ybb7wx67BMfgcAAACQkGFOyUCGYeTd3XPPPadv\nfOMbmjRpUtY4EhMAAAAACWOk1j1S62uZQ2w2myKRSPw4EonIbrebxm7evDnna1wSk98BAAAAnGYY\nhvr3mZy/LHnye19fn0pLS7V9+3aVlJRo/vz5ppPfP/30U82aNUsdHR0qKirKOjYVEwAAAAAJebzK\nZbFYVFdXp8rKSsViMVVXV8vr9aq+vl6SFAgEJElbtmxRZWVlzqREomICAAAA4DTDMNT/vsn5ryVX\nTIYDFRMAAAAACSarcJ0NJCYAAAAA4k7m8SrXcDgj+5jk2o4eGOmcTqfmzZsnn8+n+fPnS5IOHTqk\nb33rW7r00kt1/fXXq7u7u8BXCQyPW2+9VcXFxZo7d278XLb7/4EHHpDb7ZbH49Hzzz9fiEsGho3Z\n83DPPffIbrfL5/PJ5/Opubk5/hnPA85FJ8eMTmtm8skBWltb5fP5NGfOHC1evDjruEOeYzLY7eiB\nkWjmzJnas2ePpkyZEj935513atq0abrzzjtVW1urw4cP68EHHyzgVQLD48UXX9T48eN188036803\n35SU+f7fv3+/vve97+mVV15RNBrVddddp7ffflujRrHfL0YGs+dh3bp1mjBhgn76058mxfI84Fxk\nGIY+7h+fdn66cSxpjkk+OUB3d7cWLlyobdu2yW63q6urS9OmTcs49pDv/MFuRw+MVKk5/tatW7Vq\n1SpJ0qpVq7Rly5ZCXBYw7BYtWqTJkycnnct0/weDQa1cuVJWq1VOp1Mul0vhcPisXzMwXMyeB8l8\n0jDPA85VJzUmraXKJwd4+umndeONN8b3N8mWlEhnIDEZzHb0wEhlGIauu+46lZeX67HHHpMkdXZ2\nqri4WJJUXFyszs7OQl4icFZluv8/+OCDpA24+J2B88X69etVVlam6urq+KuNPA84V32morSWKp8c\noK2tTYcOHdI111yj8vJy/fnPf8467pATk8FsRw+MVDt37tTevXvV3NysRx99VC+++GLS54Zh8Kzg\nvJXr/ufZwEi3du1avfvuu9q3b59mzJihn/3sZxljeR5wLsinYpLPvfr555/rtddeU1NTk7Zt26bf\n/OY3amtryxg/5FW5BrMdPTBSzZgxQ5I0ffp0rVixQuFwWMXFxfroo4908cUX68MPP9RFF11U4KsE\nzp5M93/q74yOjg7ZbLZCXSZwVgz8+f/DH/5QN9xwgySeB5y7ejVar7Ye16utPRlj8skBHA6Hpk2b\npqKiIhUVFemb3/ymXn/9dbndbtM+h1wxKS8vV1tbm9rb29Xb26vGxkb5/f6hdgv81+jp6dHRo0cl\nScePH9fzzz+vuXPnyu/366mnnpIkPfXUU1q+fHkhLxM4qzLd/36/X5s3b1Zvb6/effddtbW1xVey\nA0aqDz/8MP73v//97/EVu3gecK7q0VjNXjxdN9/z9XhLlU8O8O1vf1svvfSSYrGYenp6tHv3bs2e\nPTvjuEOumGTajh44X3R2dmrFihWSpL6+Pn3/+9/X9ddfr/Lyct10003auHGjnE6nnnnmmQJfKTA8\nVq5cqR07dqirq0sOh0P33nuv7rrrLtP7f/bs2brppps0e/ZsWSwW/fGPf+TVFYwoqc/DunXr1Nra\nqn379skwDM2cOVP19fWSeB5w7uqV+fLAA2XKAb68vwOBgDwej5YsWaJ58+Zp1KhRWrNmTdbEZMjL\nBQMAAAAYGQzDUHP/4rTzS41W09XlziR2fgcAAAAQl0/FZDiQmAAAAACI69HYgoybc/J7PlvNAwAA\nABgZTmp0WjOTK09obW3VxIkT5fP55PP5dN9992UdN2vFJBaLqaamJmmreb/fz+R2AAAAYITqNdm3\nJFW+ecLVV1+trVu35jVu1opJPlvNAwAAABg5elSU1lLlmycMZsJ81oqJ2Vbzu3fvTooxDKek9/Ie\nEAAAADg/fF39/e2FvohBy6dikl+eYOjll19WWVmZbDabHn744a++j0l+a2m/J+nXkp6UdEvap7/W\nuqTjdfp11t5S47NJ7WuwYw127MH2N9j+hzL2UL/PZzp+OPsa7PdxOPvPdA8+KbOnYfD95TKUeyrX\n+Gf63ynX85rtWs7EeIP52lyG+nPsTP4cHGzfQx3rq/b/pIb+TAz3z+ihGuw9P5TnbShj5eNc/h1R\nyOcrtf+v/vw8qa/yRAz38zoYg/13GGr8cP7/MZehfd/P7O/psyXTnJKB8skTLr/8ckUiEY0dO1bN\nzc1avny53n777YzxWROTfLaaP6VVUvfpP52nGwAAAHD+aJd06v/D/90+01i919qu91ozvxWVT54w\nYcKE+N+XLl2qH/3oRzp06JCmTJli2mfWxGTgVvMlJSVqbGxUQ0ODSeRinfqnWJytOwAAAGDEckpK\n/v/wjkJcxpCd1GhdvPhSXbz40vi5l9a9kBSTT57Q2dmpiy66SIZhKBwOq7+/P2NSIuWx83tzc7Pu\nuOOO+Fbzv/jFL5I+X7x4sXbs+O/8pgMAAADD5eqrr1Zra2uhL2NQDMPQ//Y/lHZ+vXFn2kR2szyh\nvr5ekhQIBPToo4/qT3/6kywWi8aOHavf/e53uvLKKzOPnSsxAQAAAHB+MAxDt/U/knZ+g3HHoFbY\n+irY+R0AAABA3GcF2vmdxAQAAABAXD6rcg2HrBssAgAAADi/9GpMWjMTCoXk8XjkdrtVW1ubsb9X\nXnlFFotFzz77bNZxqZgAAAAAiOvJ41WuWCymmpoatbS0yGazqaKiQn6/X16vNy3u5z//uZYsWZJz\njgoVEwAAAABxvRqd1lKFw2G5XC45nU5ZrVZVVVUpGAymxa1fv17f/e53NX369JzjkpgAAAAAiDup\nMWktVTQalcPhiB/b7XZFo9G0mGAwqLVr10rKvVs8r3IBAAAAiDOrkKTKlWRI0h133KEHH3xQhmGo\nv78/56tcJCYAAAAA4no0VsdbX1VP66sZY2w2myKRSPw4EonIbrcnxezZs0dVVVWSpK6uLjU3N8tq\ntcrv95v2yQaLAAAAACSdqoRc0v9/aeffMeYkVTz6+vpUWlqq7du3q6SkRPPnz1dDQ0Pa5PcvrV69\nWjfccIO+853vZBybigkAAACAOLM5JaksFovq6upUWVmpWCym6upqeb1e1dfXS5ICgcCgx6ViAgAA\nAEDSqYrJ1FhH2vlPLrDnnCMyVFRMAAAAAMSdPJG7YjIcSEwAAAAAxPWeyL0q13BgHxMAAAAAcb3H\nxqY1M6FQSB6PR263W7W1tWmfB4NBlZWVyefz6YorrtA///nPrOMyxwQAAACApNP7k7xlkh6UGklz\nTGKxmEpLS9XS0iKbzaaKioq0VbmOHz+ucePGSZLefPNNrVixQgcPHsw4NhUTAAAAAAknTFqKcDgs\nl8slp9Mpq9WqqqoqBYPBpJgvkxJJOnbsmKZNm5Z1WBITAAAAAAmfmbQU0WhUDocjfmy32xWNRtPi\ntmzZIq/Xq6VLl+oPf/hD1mFJTAAAAAAkHDdpKQzDyKur5cuX68CBA3ruuef0gx/8IGssq3IBAAAA\nSDgh6Y1W6c3WjCE2m02RSCR+HIlEZLfbM8YvWrRIfX19+uSTTzR16lTTGBITAAAAAAknJF26+FT7\n0tPrkkLKy8vV1tam9vZ2lZSUqLGxUQ0NDUkx77zzjmbNmiXDMPTaa69JUsakRCIxAQAAADDQsdwh\nFotFdXV1qqysVCwWU3V1tbxer+rr6yVJgUBAf/vb37Rp0yZZrVaNHz9emzdvztonywUDAAAAkHR6\n7ki9SXoQSF4ueDhQMQEAAACQYLI88NlAYgIAAAAgwWR54LOB5YIBAAAAJOSxXLAkhUIheTweud1u\n1dbWpn3+17/+VWVlZZo3b54WLlyoN954I+uwVEwAAAAAJORRMYnFYqqpqVFLS4tsNpsqKirk9/vl\n9XrjMbNmzdILL7ygiRMnKhQK6bbbbtOuXbsy9knFBAAAAEDCSZOWIhwOy+Vyyel0ymq1qqqqSsFg\nMCnmqquu0sSJEyVJCxYsUEdHR9ZhSUwAAAAAJBwzaSmi0agcDkf82G63KxqNZuxy48aNWrZsWdZh\neZULAAAAQEIeq3IZhpF3d//617/0+OOPa+fOnVnjSEwAAAAAJJyQ9HGr1NWaMcRmsykSicSPI5GI\n7HZ7Wtwbb7yhNWvWKBQKafLkyVmHZYNFAAAAAJJOV0KWmqQHzckbLPb19am0tFTbt29XSUmJ5s+f\nr4aGhqTJ7++//76uvfZa/eUvf9GVV16Zc2wqJgAAAAASTOaUpLJYLKqrq1NlZaVisZiqq6vl9XpV\nX18vSQoEArr33nt1+PBhrV27VpJktVoVDocz9knFBAAAAICk0xWTcpP04NXkislwoGICAAAAIMFk\neeCzgcQEAAAAQEIer3INBxITAAAAAAl5LBc8HNhgEQAAAEDCCZNmIhQKyePxyO12q7a2Nu3z//zn\nP7rqqqt04YUX6re//W3OYamYAAAAAEg4mjskFouppqZGLS0tstlsqqiokN/vT1oueOrUqVq/fr22\nbNmS17BUTAAAAAAk9Jm0FOFwWC6XS06nU1arVVVVVQoGg0kx06dPV3l5uaxWa17DkpgAAAAAGJRo\nNCqHwxE/ttvtikajQ+qTV7kAAAAADPB5zgjDMM74qCQmAAAAAAY4KuklSTszRthsNkUikfhxJBKR\n3W4f0qgkJgAAAAAG6JF0+en2pYeSIsrLy9XW1qb29naVlJSosbFRDQ0Npr3lu2M8iQkAAACAAT7L\nGWGxWFRXV6fKykrFYjFVV1fL6/Wqvr5ekhQIBPTRRx+poqJCR44c0ahRo/T73/9e+/fv1/jx4037\nNPrzTWEAAAAAjGin5o7sMfnkirwrH18VFRMAAAAAA+SumAwHEhMAAAAAAxQmMWEfEwAAAAADfGbS\n0oVCIXk8HrndbtXW1prG3H777XK73SorK9PevXuzjkpiAgAAAGCAIyYtWSwWU01NjUKhkPbv36+G\nhgYdOHAgKaapqUkHDx5UW1ubNmzYoLVr12YdlcQEAAAAwAC5KybhcFgul0tOp1NWq1VVVVUKBoNJ\nMVu3btWqVaskSQsWLFB3d7c6OzszjkpiAgAAAGCA3IlJNBqVw+GIH9vtdkWj0ZwxHR0dGUdl8jsA\nAACAAdJf3Up1alnh3FKXGM72dSQmAAAAAAb4VdqZ1E0RbTabIpFI/DgSichut2eN6ejokM1myzgq\nr3IBAAAAkHSqwmHWjh49mhRXXl6utrY2tbe3q7e3V42NjfL7/Ukxfr9fmzZtkiTt2rVLkyZNUnFx\nccaxqZgAAAAAGBSLxaK6ujpVVlYqFoupurpaXq9X9fX1kqRAIKBly5apqalJLpdL48aN0xNPPJG1\nT6N/uPeWBwAAAIAceJULAAAAQMGRmAAAAAAoOBITAAAAAAVHYgIAAACg4EhMAAAAABQciQkAAACA\ngiMxAQAAAFBwJCYAAAAACu7/AaDC2kW2NU70AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x8b8b290>"
       ]
      }
     ],
     "prompt_number": 125
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print chosenSNPS2\n",
      "print m2\n",
      "pb.matshow(X[:,chosenSNPS2].transpose())\n",
      "pb.colorbar()\n",
      "pb.show()\n",
      "\n",
      "m2.plot()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[3]\n",
        "\n",
        "Log-likelihood: -8.625e+02\n",
        "\n",
        "       Name        |    Value    |  Constraints  |  Ties  |  prior  \n",
        "--------------------------------------------------------------------\n",
        "  linear_variance  |  1215.0372  |     (+ve)     |        |         \n",
        "  noise_variance   |  717.4651   |     (+ve)     |        |         \n",
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAACICAYAAAASyf8xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFn9JREFUeJzt3X9sVfX9x/HXwXvB8uPLbyu9924XvNfeS4B6tQUNY6Jx\nFki8g7mYsmUidnjD1q8zW+KcWzJxRq1xixt1SzH4g22WmsxxcWkvpmxFxcBVBPU7mBZj9faqjRUq\nPyrUXvv9A7y3995zf9RSLivPR/IJnHPf/XxOyjkt77zP5/Mx+vv7+wUAAAAABTSq0BcAAAAAACQm\nAAAAAAqOxAQAAABAwZGYAAAAACg4EhMAAAAABUdiAgAAAKDgSEwAAAAADMqtt96q4uJizZ07N2PM\n7bffLrfbrbKyMu3duzdnnyQmAAAAAAZl9erVCoVCGT9vamrSwYMH1dbWpg0bNmjt2rU5+yQxAQAA\nADAoixYt0uTJkzN+vnXrVq1atUqStGDBAnV3d6uzszNrnyQmAAAAAM6oaDQqh8MRP7bb7ero6Mj6\nNSQmAAAAAM64/v7+pGPDMLLGk5gAAAAAkCQVGYYMkzZhwoRB9WOz2RSJROLHHR0dstlsWb/G8pWu\nGAAAAMCIc0LSfSbnf3Xs2KD68fv9qqurU1VVlXbt2qVJkyapuLg469eQmAAAAACI+588YlauXKkd\nO3aoq6tLDodD69at0+effy5JCgQCWrZsmZqamuRyuTRu3Dg98cQTOfs0+lNf/gIAAABwXjIMQ4+Z\nnF+j9DkjZxoVEwAAAABxRQUal8QEAAAAQFw+r3INB1blAgAAABBXZNLMhEIheTweud1u1dbWpn1+\n+PBhrVixQmVlZVqwYIH+/e9/Zx2XxAQAAABAXD6JSSwWU01NjUKhkPbv36+GhgYdOHAgKeb+++/X\n5Zdfrtdff12bNm3ST37yk6zjkpgAAAAAiMsnMQmHw3K5XHI6nbJaraqqqlIwGEyKOXDggK655hpJ\nUmlpqdrb2/Xxxx9nHJfEBAAAAEDcBJOWKhqNyuFwxI/tdrui0WhSTFlZmZ599llJpxKZ9957Tx0d\nHRnHJTEBAAAAEDfWpKUyDCNnP3fddZe6u7vl8/lUV1cnn8+nCy64IGM8q3IBAAAAiCuySC99Ie0c\nuG1JyhYmNptNkUgkfhyJRGS325NiJkyYoMcffzx+PHPmTM2aNSvjuGywCAAAAEDSqUpI78T086M/\nTd5gsa+vT6Wlpdq+fbtKSko0f/58NTQ0yOv1xmM+/fRTFRUVafTo0Xrssce0c+dOPfnkkxnHpmIC\nAAAAIM46JneMxWJRXV2dKisrFYvFVF1dLa/Xq/r6eklSIBDQ/v37dcstt8gwDM2ZM0cbN27M2icV\nEwAAAACSTlVM+r9mcv795IrJcKBiAgAAACAhj4rJcCAxAQAAAJAwrjDDkpgAAAAASLiwMMOyjwkA\nAACAhDEmzUQoFJLH45Hb7VZtbW3a511dXVqyZIkuu+wyzZkzJ+uKXBKT3wEAAACcZhiG+peZnG9K\nnvwei8VUWlqqlpYW2Ww2VVRUpC0XfM899+jkyZN64IEH1NXVpdLSUnV2dspiMX9pi4oJAAAAgIQL\nTVqKcDgsl8slp9Mpq9WqqqoqBYPBpJgZM2boyJEjkqQjR45o6tSpGZMSiTkmAAAAAAbKY1WuaDQq\nh8MRP7bb7dq9e3dSzJo1a3TttdeqpKRER48e1TPPPJO1TyomAAAAABLGmbQUhmHk7Ob+++/XZZdd\npg8++ED79u3Tj3/8Yx09ejRjPBUTAAAAAAkXSq0dUms0c4jNZlMkEokfRyIR2e32pJiXX35Zv/zl\nLyVJl1xyiWbOnKm33npL5eXlpn2SmAAAAABIGCMtvuRU+9K6V5JDysvL1dbWpvb2dpWUlKixsVEN\nDQ1JMR6PRy0tLVq4cKE6Ozv11ltvadasWRmHJTEBAAAAkJDHHBOLxaK6ujpVVlYqFoupurpaXq9X\n9fX1kqRAIKC7775bq1evVllZmb744gs99NBDmjJlSsY+WS4YAAAAgKTTywXfb3L+7uTlgocDFRMA\nAAAACXlUTIYDiQkAAACAhAIlJiwXDAAAACAhj+WCJSkUCsnj8cjtdqu2tjbt84cfflg+n08+n09z\n586VxWJRd3d3xmGZYwIAAABA0uk5Jg0m51cmzzGJxWIqLS1VS0uLbDabKioq1NDQIK/Xa9rvP/7x\nDz3yyCNqaWnJODYVEwAAAAAJY0xainA4LJfLJafTKavVqqqqKgWDwYxdPv3001q5cmXWYUlMAAAA\nACRcaNJSRKNRORyO+LHdblc0ar4jY09Pj7Zt26Ybb7wx67BMfgcAAACQkGFOyUCGYeTd3XPPPadv\nfOMbmjRpUtY4EhMAAAAACWOk1j1S62uZQ2w2myKRSPw4EonIbrebxm7evDnna1wSk98BAAAAnGYY\nhvr3mZy/LHnye19fn0pLS7V9+3aVlJRo/vz5ppPfP/30U82aNUsdHR0qKirKOjYVEwAAAAAJebzK\nZbFYVFdXp8rKSsViMVVXV8vr9aq+vl6SFAgEJElbtmxRZWVlzqREomICAAAA4DTDMNT/vsn5ryVX\nTIYDFRMAAAAACSarcJ0NJCYAAAAA4k7m8SrXcDgj+5jk2o4eGOmcTqfmzZsnn8+n+fPnS5IOHTqk\nb33rW7r00kt1/fXXq7u7u8BXCQyPW2+9VcXFxZo7d278XLb7/4EHHpDb7ZbH49Hzzz9fiEsGho3Z\n83DPPffIbrfL5/PJ5/Opubk5/hnPA85FJ8eMTmtm8skBWltb5fP5NGfOHC1evDjruEOeYzLY7eiB\nkWjmzJnas2ePpkyZEj935513atq0abrzzjtVW1urw4cP68EHHyzgVQLD48UXX9T48eN188036803\n35SU+f7fv3+/vve97+mVV15RNBrVddddp7ffflujRrHfL0YGs+dh3bp1mjBhgn76058mxfI84Fxk\nGIY+7h+fdn66cSxpjkk+OUB3d7cWLlyobdu2yW63q6urS9OmTcs49pDv/MFuRw+MVKk5/tatW7Vq\n1SpJ0qpVq7Rly5ZCXBYw7BYtWqTJkycnnct0/weDQa1cuVJWq1VOp1Mul0vhcPisXzMwXMyeB8l8\n0jDPA85VJzUmraXKJwd4+umndeONN8b3N8mWlEhnIDEZzHb0wEhlGIauu+46lZeX67HHHpMkdXZ2\nqri4WJJUXFyszs7OQl4icFZluv8/+OCDpA24+J2B88X69etVVlam6urq+KuNPA84V32morSWKp8c\noK2tTYcOHdI111yj8vJy/fnPf8467pATk8FsRw+MVDt37tTevXvV3NysRx99VC+++GLS54Zh8Kzg\nvJXr/ufZwEi3du1avfvuu9q3b59mzJihn/3sZxljeR5wLsinYpLPvfr555/rtddeU1NTk7Zt26bf\n/OY3amtryxg/5FW5BrMdPTBSzZgxQ5I0ffp0rVixQuFwWMXFxfroo4908cUX68MPP9RFF11U4KsE\nzp5M93/q74yOjg7ZbLZCXSZwVgz8+f/DH/5QN9xwgySeB5y7ejVar7Ye16utPRlj8skBHA6Hpk2b\npqKiIhUVFemb3/ymXn/9dbndbtM+h1wxKS8vV1tbm9rb29Xb26vGxkb5/f6hdgv81+jp6dHRo0cl\nScePH9fzzz+vuXPnyu/366mnnpIkPfXUU1q+fHkhLxM4qzLd/36/X5s3b1Zvb6/effddtbW1xVey\nA0aqDz/8MP73v//97/EVu3gecK7q0VjNXjxdN9/z9XhLlU8O8O1vf1svvfSSYrGYenp6tHv3bs2e\nPTvjuEOumGTajh44X3R2dmrFihWSpL6+Pn3/+9/X9ddfr/Lyct10003auHGjnE6nnnnmmQJfKTA8\nVq5cqR07dqirq0sOh0P33nuv7rrrLtP7f/bs2brppps0e/ZsWSwW/fGPf+TVFYwoqc/DunXr1Nra\nqn379skwDM2cOVP19fWSeB5w7uqV+fLAA2XKAb68vwOBgDwej5YsWaJ58+Zp1KhRWrNmTdbEZMjL\nBQMAAAAYGQzDUHP/4rTzS41W09XlziR2fgcAAAAQl0/FZDiQmAAAAACI69HYgoybc/J7PlvNAwAA\nABgZTmp0WjOTK09obW3VxIkT5fP55PP5dN9992UdN2vFJBaLqaamJmmreb/fz+R2AAAAYITqNdm3\nJFW+ecLVV1+trVu35jVu1opJPlvNAwAAABg5elSU1lLlmycMZsJ81oqJ2Vbzu3fvTooxDKek9/Ie\nEAAAADg/fF39/e2FvohBy6dikl+eYOjll19WWVmZbDabHn744a++j0l+a2m/J+nXkp6UdEvap7/W\nuqTjdfp11t5S47NJ7WuwYw127MH2N9j+hzL2UL/PZzp+OPsa7PdxOPvPdA8+KbOnYfD95TKUeyrX\n+Gf63ynX85rtWs7EeIP52lyG+nPsTP4cHGzfQx3rq/b/pIb+TAz3z+ihGuw9P5TnbShj5eNc/h1R\nyOcrtf+v/vw8qa/yRAz38zoYg/13GGr8cP7/MZehfd/P7O/psyXTnJKB8skTLr/8ckUiEY0dO1bN\nzc1avny53n777YzxWROTfLaaP6VVUvfpP52nGwAAAHD+aJd06v/D/90+01i919qu91ozvxWVT54w\nYcKE+N+XLl2qH/3oRzp06JCmTJli2mfWxGTgVvMlJSVqbGxUQ0ODSeRinfqnWJytOwAAAGDEckpK\n/v/wjkJcxpCd1GhdvPhSXbz40vi5l9a9kBSTT57Q2dmpiy66SIZhKBwOq7+/P2NSIuWx83tzc7Pu\nuOOO+Fbzv/jFL5I+X7x4sXbs+O/8pgMAAADD5eqrr1Zra2uhL2NQDMPQ//Y/lHZ+vXFn2kR2szyh\nvr5ekhQIBPToo4/qT3/6kywWi8aOHavf/e53uvLKKzOPnSsxAQAAAHB+MAxDt/U/knZ+g3HHoFbY\n+irY+R0AAABA3GcF2vmdxAQAAABAXD6rcg2HrBssAgAAADi/9GpMWjMTCoXk8XjkdrtVW1ubsb9X\nXnlFFotFzz77bNZxqZgAAAAAiOvJ41WuWCymmpoatbS0yGazqaKiQn6/X16vNy3u5z//uZYsWZJz\njgoVEwAAAABxvRqd1lKFw2G5XC45nU5ZrVZVVVUpGAymxa1fv17f/e53NX369JzjkpgAAAAAiDup\nMWktVTQalcPhiB/b7XZFo9G0mGAwqLVr10rKvVs8r3IBAAAAiDOrkKTKlWRI0h133KEHH3xQhmGo\nv78/56tcJCYAAAAA4no0VsdbX1VP66sZY2w2myKRSPw4EonIbrcnxezZs0dVVVWSpK6uLjU3N8tq\ntcrv95v2yQaLAAAAACSdqoRc0v9/aeffMeYkVTz6+vpUWlqq7du3q6SkRPPnz1dDQ0Pa5PcvrV69\nWjfccIO+853vZBybigkAAACAOLM5JaksFovq6upUWVmpWCym6upqeb1e1dfXS5ICgcCgx6ViAgAA\nAEDSqYrJ1FhH2vlPLrDnnCMyVFRMAAAAAMSdPJG7YjIcSEwAAAAAxPWeyL0q13BgHxMAAAAAcb3H\nxqY1M6FQSB6PR263W7W1tWmfB4NBlZWVyefz6YorrtA///nPrOMyxwQAAACApNP7k7xlkh6UGklz\nTGKxmEpLS9XS0iKbzaaKioq0VbmOHz+ucePGSZLefPNNrVixQgcPHsw4NhUTAAAAAAknTFqKcDgs\nl8slp9Mpq9WqqqoqBYPBpJgvkxJJOnbsmKZNm5Z1WBITAAAAAAmfmbQU0WhUDocjfmy32xWNRtPi\ntmzZIq/Xq6VLl+oPf/hD1mFJTAAAAAAkHDdpKQzDyKur5cuX68CBA3ruuef0gx/8IGssq3IBAAAA\nSDgh6Y1W6c3WjCE2m02RSCR+HIlEZLfbM8YvWrRIfX19+uSTTzR16lTTGBITAAAAAAknJF26+FT7\n0tPrkkLKy8vV1tam9vZ2lZSUqLGxUQ0NDUkx77zzjmbNmiXDMPTaa69JUsakRCIxAQAAADDQsdwh\nFotFdXV1qqysVCwWU3V1tbxer+rr6yVJgUBAf/vb37Rp0yZZrVaNHz9emzdvztonywUDAAAAkHR6\n7ki9SXoQSF4ueDhQMQEAAACQYLI88NlAYgIAAAAgwWR54LOB5YIBAAAAJOSxXLAkhUIheTweud1u\n1dbWpn3+17/+VWVlZZo3b54WLlyoN954I+uwVEwAAAAAJORRMYnFYqqpqVFLS4tsNpsqKirk9/vl\n9XrjMbNmzdILL7ygiRMnKhQK6bbbbtOuXbsy9knFBAAAAEDCSZOWIhwOy+Vyyel0ymq1qqqqSsFg\nMCnmqquu0sSJEyVJCxYsUEdHR9ZhSUwAAAAAJBwzaSmi0agcDkf82G63KxqNZuxy48aNWrZsWdZh\neZULAAAAQEIeq3IZhpF3d//617/0+OOPa+fOnVnjSEwAAAAAJJyQ9HGr1NWaMcRmsykSicSPI5GI\n7HZ7Wtwbb7yhNWvWKBQKafLkyVmHZYNFAAAAAJJOV0KWmqQHzckbLPb19am0tFTbt29XSUmJ5s+f\nr4aGhqTJ7++//76uvfZa/eUvf9GVV16Zc2wqJgAAAAASTOaUpLJYLKqrq1NlZaVisZiqq6vl9XpV\nX18vSQoEArr33nt1+PBhrV27VpJktVoVDocz9knFBAAAAICk0xWTcpP04NXkislwoGICAAAAIMFk\neeCzgcQEAAAAQEIer3INBxITAAAAAAl5LBc8HNhgEQAAAEDCCZNmIhQKyePxyO12q7a2Nu3z//zn\nP7rqqqt04YUX6re//W3OYamYAAAAAEg4mjskFouppqZGLS0tstlsqqiokN/vT1oueOrUqVq/fr22\nbNmS17BUTAAAAAAk9Jm0FOFwWC6XS06nU1arVVVVVQoGg0kx06dPV3l5uaxWa17DkpgAAAAAGJRo\nNCqHwxE/ttvtikajQ+qTV7kAAAAADPB5zgjDMM74qCQmAAAAAAY4KuklSTszRthsNkUikfhxJBKR\n3W4f0qgkJgAAAAAG6JF0+en2pYeSIsrLy9XW1qb29naVlJSosbFRDQ0Npr3lu2M8iQkAAACAAT7L\nGWGxWFRXV6fKykrFYjFVV1fL6/Wqvr5ekhQIBPTRRx+poqJCR44c0ahRo/T73/9e+/fv1/jx4037\nNPrzTWEAAAAAjGin5o7sMfnkirwrH18VFRMAAAAAA+SumAwHEhMAAAAAAxQmMWEfEwAAAAADfGbS\n0oVCIXk8HrndbtXW1prG3H777XK73SorK9PevXuzjkpiAgAAAGCAIyYtWSwWU01NjUKhkPbv36+G\nhgYdOHAgKaapqUkHDx5UW1ubNmzYoLVr12YdlcQEAAAAwAC5KybhcFgul0tOp1NWq1VVVVUKBoNJ\nMVu3btWqVaskSQsWLFB3d7c6OzszjkpiAgAAAGCA3IlJNBqVw+GIH9vtdkWj0ZwxHR0dGUdl8jsA\nAACAAdJf3Up1alnh3FKXGM72dSQmAAAAAAb4VdqZ1E0RbTabIpFI/DgSichut2eN6ejokM1myzgq\nr3IBAAAAkHSqwmHWjh49mhRXXl6utrY2tbe3q7e3V42NjfL7/Ukxfr9fmzZtkiTt2rVLkyZNUnFx\nccaxqZgAAAAAGBSLxaK6ujpVVlYqFoupurpaXq9X9fX1kqRAIKBly5apqalJLpdL48aN0xNPPJG1\nT6N/uPeWBwAAAIAceJULAAAAQMGRmAAAAAAoOBITAAAAAAVHYgIAAACg4EhMAAAAABQciQkAAACA\ngiMxAQAAAFBwJCYAAAAACu7/AaDC2kW2NU70AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x91e6690>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD7CAYAAACL+TRnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0W/WZN/CvVsu7JS+yLMlLvBsntmlISCHgDCQkLJlQ\naEo6pcywnJbp6ZTSaWHaoYR2IOGdznQYpgw9M+kp7ekU+pbSZKZNSmHqlLXhJbYTYsd2vEqyJC+S\nV21X9/7eP2RfW9hOZNnW+nzO4RxHP1n3ibGeXD33uc9PwhhjIIQQktCk0Q6AEELIxqNkTwghSYCS\nPSGEJAFK9oQQkgQo2RNCSBKgZE8IIUlAHq0DNzY2or29PVqHJ4SQuHTjjTeipaVl1d8XtTP79vZ2\nMMbi9r8nn3wy6jFQ/NGPg+KPv//iOXbGGE6fPh1WzqUyDiGEJAFK9oQQkgQo2Yepubk52iGsCcUf\nXRR/9MRz7GshYYxFZTaORCJBlA5NCCFxK9zcSWf2hBCSBEJK9vfffz+0Wi02b94sPnb48GEYDAY0\nNTWhqakJJ0+eFNeOHDmCyspK1NTU4PXXX1//qAkhhKxKSGWct956CxkZGfj85z+P8+fPAwCeeuop\nZGZm4tFHHw16bkdHBz772c/igw8+gMViwc0334zu7m5IpcH/rlAZhxBCVm9Dyzg7d+6EWq1e8vhy\nBzx+/DgOHToEhUKB0tJSVFRU4MyZM6sOjBBCyPpZU83++eefR0NDAx544AFMTEwAAIaHh2EwGMTn\nGAwGWCyWtUVJCCFkTcIel/Dwww/j29/+NgDgiSeewNe+9jUcO3Zs2edKJJJlHz98+LD4dXNzc9K2\nRBFCyEpaWlrCGo/wcWEn+4KCAvHrBx98EHfccQcAQK/Xw2QyiWtmsxl6vX7Z11ic7AkhhCz18RPh\np556KqzXCbuMY7Vaxa9fe+01sVNn//79ePnll+Hz+dDf34+enh5s27Yt3MMQQghZByGd2R86dAin\nT5/G2NgYjEYjnnrqKbS0tKCtrQ0SiQRlZWX44Q9/CACoq6vDwYMHUVdXB7lcjhdeeGHFMg4hhJDI\noDtoCSEkjtAdtIQQQlZEyZ4QQpJA1HaqIoQQcmV+XkDHwCjMIzPw+PiwX4eSPSGExBDGGGyOWZy7\nZIeHE+DleBTmZaFYn7umZhdK9oQQEmUuD4ezXVY4Z3zw+gWkpihgLNRAJlu/Sjsle0IIiTBBYOgc\nGMWgfRoejgdjQIleg02a7A07JiV7QgiJAJtjBucvjcDF8fByAgrUGWsuzawGJXtCCNkAbq8/UJqZ\n9sLjF6BKkaO4UL2upZnVoGRPCCHrQBAYuobGMWCbhNvHQ2BAcZEaZeqsaIcGgJI9IYSEbWzChdYe\nO9w+P3wcQ646HcaiyJVmVoOSPSGEhMjr86O1y4axaQ+8nACFQo7SouiVZlaDkj0hhKxAEBh6TA70\nWSfg8fHghbnSTE5mtENbNUr2hBCyyPikG609Nri9PDwcD01O7JZmVoOSPSEkqXl9frR22zE+5YZn\nrjRTUqSGPA5KM6tByZ4QklQExnBpyIE+6yTcHA9eYDAWqlFqzIh2aBuKkj0hJOGNT7nR1m3DrNcP\nHydAk5MBQ5Em7kszq0HJnhCScLw+P9p67BibnCvNyGUo0WtQmGClmdWgZE8IiXtLSjM8g1GX+KWZ\n1aBkTwiJS4tLM16fgFx18pVmVoOSPSEkLiwuzXg5AfK5rplkLs2sBiV7QkhMEhhDr9mBvuHArBmO\nZyim0kzYKNkTQmJGUGmGE6DJTkeRTgMplWbWjJI9ISRqPl6akSnkKNGpUSin0sx6o2RPCImY+dJM\n7/AE3F4efgHUNRMhlOwJIRvKMeVG68dKM3pdLpVmIoySPSFkXfk4Hq3dtqDSTHFhDgoVsmiHltQo\n2RNC1kRgDH0WB3otga4Zv8BgSIJZM/GGkj0hZNUcc10zMx4/PJyA3Bzqmol1lOwJIVfk43i09dgw\nOrFwQ5OxMAdaKs3EDUr2hJAlFpdmXF4//AJQrMuh0kwco2RPCAGwqDQzNwZYTTc0JRRK9oQkKbE0\nM+mB18dDJpejWEelmURFyZ6QJMEYQ6/FgV7LBNw+AX6ewajLQamBSjPJgJI9IQlscWnG6xOgzklH\nEd3QlJQo2ROSQDi/gNZuK5VmPubBQ3+OF3/6KuTyQMrz+/344r134T9/fjzKkUVOSNOG7r//fmi1\nWmzevFl8zOFwYPfu3aiqqsKePXswMTEhrh05cgSVlZWoqanB66+/vv5RE0IABEozlywO/O5Pffj1\nWz349duXwGQKlBryUL1Ji4riXCgp0eP9d05ja7UOfr8ffr8fW6t1eP+d03jw0J9HO7yICSnZ/9Vf\n/RVOnToV9NjRo0exe/dudHd346abbsLRo0cBAB0dHXjllVfQ0dGBU6dO4a//+q8hCML6R05IknJM\nufG//68fJ97uwS9bujEw4kKRToPK0gLUlWuRkZYS7RBjimN8DADg5zhsrdYFkj7HBa0lg5DKODt3\n7sTAwEDQYydOnMDp06cBAPfddx+am5tx9OhRHD9+HIcOHYJCoUBpaSkqKipw5swZXHvttesePCHJ\ngEoza3Op56L49XySX24t0YVds7fb7dBqtQAArVYLu90OABgeHg5K7AaDARaLZY1hEpI8GGPoHXai\n1+yEa24McDF1zYTt5n378fr/vLbiWrJYlwu0Eonkspv8rrR2+PBh8evm5mY0NzevRziExB3ntBtt\n3XZMezj4OAE52dQ1s16OfP/FFZP9ke+/GOFoVu/Me2/hg/feXvPrhJ3stVotbDYbCgsLYbVaUVBQ\nAADQ6/UwmUzi88xmM/R6/bKvsTjZE5JMxNKM0w2PP1CaKdGpUUClmXX3hXs/ddm1H//iNxGMZvW2\n7diJbTt2in/+9395NqzXCXvvr/379+Oll14CALz00ks4cOCA+PjLL78Mn8+H/v5+9PT0YNu2beEe\nhpCEIMzd0LS4a4aXKlBanI+aTYWoLM5L+q6ZjTI+NhrWWqIJ6cz+0KFDOH36NMbGxmA0GvGd73wH\njz/+OA4ePIhjx46htLQUv/jFLwAAdXV1OHjwIOrq6iCXy/HCCy9ctsRDSKJyTrvR2m3HjJuDh2bN\nRM3QQJ/4dWZWDgBgempiyVqikzDGWFQOLJEgSocmZEMsLs14/QKkchlKdGo6Y4+yb371izjx6svI\nzMrB2+29AIDrG8oxPTWB/Xfdg2fioG6/WH1xTli5k5I9IWEKjAF2otcyAZfXD04AjIU5yEqnPvdY\n882vfhHf/d4PIJMF/uHleR5P/O2X4i7RA5TsCYkI57QbZ7tsmPX44fULyMlKR2FeBpVmSMSEm+xp\nNg4hl+HjeLT32DDidMPD8ZAp5Cgt0qBQHnZvAyFRQcmekEXmSzOXLE64vbxYmiktphuaSHyjZE+S\nnnPajbPddsx6OHh8AtRZadDTDU0kwVCyJ0nHx/Fo7bFhjEozJA54fX5cuGRFa4cJrZ3msF+Hkj1J\neAJj6DM7cckyATfHg+MZlWZIzBIEhn7zGFo7zWjtNKHjkg2cn1/z61I3DklIjmk3WrvmSjNcoDSj\ny8+kG/xITBpzzqC104y2ThPaLpoxOe0JWi8vzkNTrRGNtQZ8bk89deOQ5DVfmhmdK83I5TKU6nOp\nNEMAxF6fvdvD4aOeYZztMKGt0wyTzRm0nqdOF5N7Y40B2Zmpaz4mJXsSlwTG0Gty4NLwBNw+Hvxc\n10wZlWbIx8zfQdvyxim81Ra4g3ZnYzmmJgMjEyKR8HlBQO/QmFh3v9hng59f2NQpNUWB+qoiNNUa\n0FRnhEGbs+6fQqmMQ+LG+Nzm2bMePzw+AbnqdGhzM6g0Qy7r2quKMTM9BQDIyg7MxplP9BmZWXj/\nwtCGHNc+NiXW3dsvWjDj8oprUokElaX5aKw1oqnWgOoyLRTy0MZq0E1VJOF4fX6c7bJhbMoDL8dD\nqZCjRK9BoYxKMyR0n7zhz/D6b34NYCHJL15bL7NuL851DaOtM3D2PjwyGbSuzc1EU10guTdUG5AR\n4bEalOxJzBAEhq6hcQzYJuH2+cHzQHGRBpuKM6MdGoljFdW1YrJfbi1cPC+ge2AErZ0mtHaY0TVg\nhyAsnHGnqZRoqNGjqdaAxlojdPlZUf0USmUcElUjzlm0XxqB2+uHhxOQr8lAvjqdSjNk3XAch6by\n/GXXWntHoVAoQnodxhiso1Ni3f1clwUuj09cl0olqCnToqkucGG1qqQAsg34FEplHBIXXB4OrV1W\nOGd98HIClAo5iovUkFNphmyQv3vki5dd+94Pjq24Pj3rQftFC1o7A10z9vHpoPWigmw0zdXdN1cX\nIT01dieeUrInG4rnBVwcHMfQyBQ8Ph48A0p0GpRpsqMdGkkS75x+I+Q1zs/jYp9drLv3DI5g8Ul0\nZnoKGqoNaKozoKnWiILc+CkxUrIn6254bBof9Y7A7Rfg5QRoNZkwFuVSaYZExY033YL/eS2wk151\nbT0AoKvzIwDADX92C0xWZ6Du3mnG+W4LPF6/+L1ymRS15YViS+QmYx5k0vj8FErJnqzZrJtDa7cV\nzhkfPByPtFQljFrNhtQrCVmtsvJKAEBlTR3+78k/YnLajYO33QDrUC9a+2bw8FMvBz2/WKdGY23g\nzL2+sgipqtBq+rGOkj1ZNZ4XcKFvFOaxGXg4PwAJivUabKLSDIlB77/7RwBAz8UOfPm7r6DXNAbr\nUODmqnHLRdTW7kFDjWGua8aAPHVi3phH3TjkihhjsIxOo6N/VCzN6PKyoM5KpdIMiTmMMQxYHOJF\n1V99714Ay+caqVSGtr4xSKXx83tM3ThkXc26ffiwy4bJWR+8fgHpqUoYCnPj6k1BkodjchZtnea5\nYWJmOKdcC4s5lcBE97Lft2vPrUnzO01n9gQA4OcFnO8dgXV8Fh6OByQSlBZpkKKk8wESezw+Dhd6\nrIFxBB0mDA47gtY12WniILHNVTrs2VYJnzd4kqQyRYUPLlrE4Wjxgs7syaowxjBkm8TFoXF4fDy8\nfgFFBTko1lPXDIk9gsDQZ14YJNbRa4XfvzBILEUpR33lwiCxYp1a/D32+/1LEj0A+LyepDrhpGSf\nRKZmvTjbZcW0m4PXz5CZpkIRlWZIjBp1zIh197aLZkzNLCRsiQSoKM4X+91rNxVCoVj+DP3eu/au\neIx779qLnx9fuQ8/kVCyT2CcX8C5XjvsDhc8Ph4SqRQlejXy8+l/O4k9Lo8PH3UPi5MizbbgoWX5\nmozAmXutEVtq9MjOCG3G+4VzreLXqtQ0AIDH7VqylujoXZ9AGGPoG3ai1zIBlzew/Z6hkEozJDbx\ngoBLg6Ni3f1inx28EDzjfUu1fq7n3QB9mDPebz/waZx49WWoUtPwpw4TAGB7nREetwu3H/j0uv19\nYh1doI1zzmk3WrtsmPH44eUE5GSloTA/E1JK7iQG2caCB4ktnfFeINbdq8sKIF+ni6extlPVWoR7\ngZaSfZzxcTzaemwYnfTA4/NDJpejRKeGcoV6JSHRNOPy4lyXZa4t0gTr6FTQemFellh331KtR0Za\n7A4SixXUjZOgBMZwyeRAv3USbh8Pv8BgLMxBqSEx7/Ij8c3P8+juHxFnzXT3j0BYlJjS05RoqNbP\nTYo0ojA/K4rRJhdK9jFofNKFth47Zr2B0owmOx1FOg2VZkjMYYxheGRSrLuf67bA7eHEdZlUiqvK\nC8W6e+UGzXgnV0bJPgZ4fH60ddswPuWBlxMgV8hRUqSm7fdITJqa8aD9olnsmhl1zAStG7Q5geRe\nZ8TmqiKkqZRRipQsRsk+CgTG0DUwjkH7FNwcD15gKNapUWqMn9nYJHlwfh6dvTax7n5paDRoxntW\nugoNtQuDxAo09HsciyjZR8iIcxbnLo3A5fPDxwnIVWfCUKShlkgScxhjGLI6xQ08zncPw+tbNONd\nLkVdeaFYd99kzKMb8+IAJftVOHXqFG655RYxQTPG8Lvf/Q579y69Q8/t5XD24sL2eylKOYw62n6P\nxKaJKRfaLprR2hEozzgmZ4PWS4o0aKw14Oo6I66q0EGVkhgz3pMJJfsQnTp1Cvv27cOXv/xlPPfc\ncwCAr3zlK3j++edx8uRJ7N69J3j7PQEo0dP2eyQ2eX1+dFwKDBJru2hGn2ksaD0nK1UcJNZYY0Bu\nTnqUIl0fb7e8getuvCnoRO2d02/i+uaboxxZ5Kw52ZeWliIrKwsymQwKhQJnzpyBw+HAZz7zGQwO\nDqK0tBS/+MUvkJOTsx7xRs0HH3wAAHj++efFx+a//umrpzCtKoM2l7bfI7FJEBgGhsfR1hGou1+4\nZIWP48V1pUKGqyqL0FQTuLBaqk+cEuPbLW/gi5+/G5+59wH8/T98DwDwD3//t3jlp8fw4k9+mTQJ\nf803VZWVleHDDz+ERqMRH/vGN76BvLw8fOMb38Czzz4Lp9OJo0ePBh84zm6qUiqV4Dhu2TW5XIG2\nvtEIR0TI5Y1PzIp197ZOMyam3UHrm4x54qyZuopCKBWJ+UG/cVMe/P7ANYfP3PsAAOCVnx4DAMjl\ncrT1ja34vbEoqjdVffzAJ06cwOnTpwEA9913H5qbm5ck+3jA8wI+6hvF8PgM/Lyw4vMYW3mNkEjx\neDl81DMcqLtfNGFo2Bm0npuTLo4iaKjWIycrLUqRRtiiTyjzSX65tUS35mQvkUhw8803QyaT4Qtf\n+AIeeugh2O12aLVaAIBWq4Xdbl9zoJEgbr83MAa3j4eX41GUnw1jUS5237ofr//Pa8t+30379kc4\nUkICpZle0+jcmbsJHb22oBnvqhQ5Nlfp0Ti3v6px0Yz3ZLJ7336cPPHqimvJYs3J/p133oFOp8Po\n6Ch2796NmpqaoHWJRBLTv2DTLh/Odlkx5eLg9QvISE2BXqtZ0kp2x50HV0z2d9x5MBKhEoIRx3Sg\n373DhPaLFkzNBs94ryxZGCRWs0kLhZxmJt124NMrJvvbkmjq5ZqTvU6nAwDk5+fjzjvvxJkzZ6DV\namGz2VBYWAir1YqCgoJlv/fw4cPi183NzWhubl5rOFfk5wWcu2SHbW7GO6RSlBapkZd3+R9F+9kP\nLru2a/e+9Q6VELjcPpzvtoh3q1rsk0HrBZrMoEFiWRmqKEUau1a61naltVhx5r238MF7b6/5ddZ0\ngdblcoHneWRmZmJ2dhZ79uzBk08+iTfeeAO5ubl47LHHcPToUUxMTETtAi1jDIO2SXQNOeD2BWa8\n6wuykJWhWtUnjs2lGjBh+dq8RCrF+QHHsmuErAbPC+gZHBGTe1ffSNCM9zSVEluqi8S2yKKC7Jj+\n5BwL6ovVAFbKNRJ8NORcYS02ReUCrd1ux5133gkgsM/jX/zFX2DPnj3YunUrDh48iGPHjomtl5E0\nMeNBa5cN0x4OPk5AVmbamgeJyeUKcL7A7O2MzMCkvpnpKXGNkHBZRyfFunv7RQtm3T5xTSqVoHbT\nwiCxqnWc8Z4sFAr5imfwigTtQFrOmv6mZWVlaGtrW/K4RqPBG29Ebl9HH8ejvceGkUkPvD4eUpkM\nxUU5KFjH/5H77rgTJ159GRmZWXjnXD8A4LotZZiZnsK+O+5ct+OQxDc/4/1sR2B/VdtY8Ix3XX42\nrq4zoLHWiC3VRUhPpRnva7H3jk/hv3/1CgDg0H0PAQB+/tJ/iGvJIi7/WRMYQ58lsP3efGnGsMEz\n3ud3tFm828075/rjdrcbEjl+nsfFPrs4SKxnYHTJjPdAx0ygNFOYRzPe19NtBz6N//7VKzh030P4\n5nf+j/j4z1/6j6S6QBs3O1WNT7nR1m3HrJeDlxOgzk6HNjeDZryTmMMYg8U+Idbdz3cPL5nxXluu\nFZN7RUk+ZFKambSREmlcQsJtS+j1+dHabcf4lBseToBCLkNxkQYKOb0pSOyZnHGj/aIlcMdqhxmj\nzuAZ78ZCtThIrL6yCKkqus5DwhP32xIKAkOPyYE+6wQ8Ph5+nqG4SINSI22/R2IPx/Ho7LOJdfde\n08dmvGeoAqWZOiOaag3IU9PvMYmuqCb7Eecs2i+NwO0LbL+Xm5NOg8RITGKMYXDYIdbdP+qxBs14\nV8hlqKuYn/FuQJmBZryT2BLVZP+ni3aU6NS0JyWJSc7JuRnvnYGzd8ekK2i9VK8R6+5XVeqgUlJp\nhsSuqCb7TYbcaB6ekCBenx8XLlnR1mnC2Q4zBizjQevqrDSx372x1gBNdnzPeCfJJWZq9vEgka7o\nk8B1on7zmDgC+MIlKzj/woz3FIUc9VU6NM6VZkpoG0kSxyjZh2h+A4R7Pv8QvvXdQK/u0098Ay//\n5D+SagOEeDfmnBHr7m0XzZic9gStlxfniXX32vLEnfFOkk9UWy8/GpqIxqHDsve6BphNgwCAez4f\nuAvv5Z8E7sIzGEtw6p32qMVGVub2zM14n6u7D1mD56DkqdODtt/LzkyNUqSEXB5jDOMTLjQ36OO7\n9TLWjYyOiF/PJ/nl1kh08YKA3qExse5+sc8WtPFMaooC9VVFuLoukOAN2hwqzZCYNePyYcjqgEwq\nQapShuKCzLBfi5J9iHweT1hrZOPZx6bEunt7lxnTs15xTSqRoLqsQKy7V5fRjHcSuzi/gAHLOPx+\nHiqlDHmZKtyxoxyqlLWnakr2Ibr/4a/gR//+LyuukchxuX0412VB69z+qsMjwTPetXmZYt29odqA\njHQaJEZik8AYhkemMDXjRopCivQUBW5o0EOzAeVESvYharh6W1hrZO14XkD3wIhYd7/Yb4cgLNQs\n01RKNNToxR2adPnZUYyWkMubmHLDYp+EQi5BWoocVUYNSnTGDZ/zRRdoQ3TztfWwDZuXXSssMuCN\n9z+KcESJizEG6+iUWHc/12WByxM8472mTIumubp7VUkB3ZhHYpbby6Hf7IBUwpCqlEOXm4bN5VrI\nw/ydDXfjJzqzD9H42MoXYS+3RkIzPesJDBK7GNhf1T4+HbSu12ajsSZQmtlMM95JDON5AYPDTnh8\nHFIUMmgylNi3vQwZqdG9w5qSfYh4fuFmm5JNlQCAwb6eJWskNJyfR1e/Ha0dgbr7pcHgGe+Z6Slo\nqDaI+6sW5IbfhUDIRmKMwT4+g/GJWaQopEhTynBNTQEKNbE1/I6SfYhuP/BpnHj1ZaSoUnHizfcB\nANtqDfB63Lg9iTZACBdjDGbbhHhR9Xy3BR7vwiAxuUyKq8p1Yt19kzGPZryTmDU968WQ1Qm5LNAS\nuUmXg10N+pgefkfJPkRdHYGavNfjxtHDj4tfL14jwSan3XODxAL7q445Z4PWi3Vq8YYmmvFOYpmX\n86Pf4gR4HiqFDPnqVBy4vgJKRfy08VKyD9HU1EJ73/z+lcutJTMf50dHry0wjqDDhF7TWNB6TmYq\nGmoWBonRjHcSqwTGYLFPYXquJTIzVYGbmwzIzlBFO7SwUbIPUdPW7bBaTCuuJaP5Ge/zG3hc6LHC\nywXPeK+v1M1NijSiVJ8b0x9zSfJijME55cbwyCRUShlSlTLUluSiuMCYMHdYU+tliPx+P66uLITA\n+4Mel8rkONtjg1yeHP9uOiZn5waJBe5YdU4Fz3gvM+QG6u61RtRV6JCiTI6fC4k/Li+HQUugJVKl\nkEOfn4GryvLDbomMFGq93GBvt7yxJNEDgMD78XbLG2i+eW8Uotp4Hh+HCz1Wse4+YHEErWuy0xYG\nidUaoM5Ki1KkhFyef64l0reoJfLW7WVIS5JrRXRmH6I7d38SPV0dy65VVtfhtd+/G+GINoYgMPSZ\nx9DWacbZDhM6eq3w+xcGiaUo5aivLBK7Zop16oT5mEsSi9gS6ZyFSilFaoocjRVa5Kvj+4SEzuw3\nmGt2Rvw6O1sNAJicdC5Zi0ejjhlxFEHbRTOmZhYGu0kkQGVJvlh3r91UCEUcdSCQ5DI164XJNgGF\nFFApZSjX5WBXo37DRxHEA0r2IfrE9k/CYh5CdrYaf2y7BAC4obECk5NOfGL7J6Mc3eq4PD581D2M\n1rlNPMy24E9Y+ZoMse6+pUaP7Aya8U5ik5fzY3DYCcYLSFHKUJCtwoHryuOqJTJSqIyzCt/86hfx\n3e/9ADJZ4BeJ53k88bdfwjPffzHKkV0eLwi4NDgq1t07e+3gheAZ71uq9eL+qnqa8U5ilCAwWOyT\nmHZ5oJpriWysKkROHLdErla4ZRxK9gnKNjYl1t3PdVkw4wqe8V5ZWoCr6wxorDWiuqwAchmdCZHY\nM98SaR2dQopcApVShpqSXBRrs5P2hIRq9klu1u0NDBKba4u0jgbf6FWYlyXOmdlSrUdGGg0SI7HJ\n5eEwOLzQEmnIy8Anb6ikyaZrRMk+Tvl5Ht39I2LdvXtgJGjGe3qaEg3V+rlNPIwozM+KYrSErMzP\nCxgadsLH+aGUS5OuJTJSKNnHCcYYhkcmxbp7e5cFbg8nrsukUlxVUYjGWgOurjOiojifzoRITBIY\nw4hjFg7nDFKUUqSnyHFtXSHyc+K7JTLWUbKPYfMz3ufHEYw4gme8G7Q54gYem6uKkKZSRilSQi5v\nxuXFkHUCirkpkWW6bOxqKKKWyAiiZB9DOD+PzvlBYp1mXBoaweLrMFnpKjTULgwSK9DQjHcSm7yc\nHwMWJ5jAQ6WQIz9HhQPXV0Ahp0+b0ULJPooYYzBZnWLd/aOe4eAZ73Ip6soLxbr7JmMeDRIjMUkQ\nGCwjk5ie9SBFHmiJvOlqI7Jps/eYQcl+Fdajz35iyoW2i2bx7H18InjGe0mRRqy7X1WhgyqFLlKR\n2MMYw8SUG8OjU0hRSJGqlKGuJBfGgqykbYmMdRuW7E+dOoVHHnkEPM/jwQcfxGOPPbZRh4qIb371\nizjx6sv4w+9P4u32PgDA9Q2bMD03y36lhO/j/Oi4ZAvU3S+a0ffxGe9ZqQuDxGoMyM1J39i/CCFh\nWmiJBFKVMujzMrCDWiLjxobcVMXzPKqrq/HGG29Ar9fjmmuuwc9//nPU1tYuHDjObqpqqtCC8wVu\nTMrMygYAMdErlClovWQHEDjj6beMo63DjNaLZlzoGYaPW9ijVqmQ4arKIjTVBAaJleo1dCZEYpKf\nF2CyOuHIx7CDAAARVUlEQVT1+ZGikEKdrkRTtY5aIqMspm6qOnPmDCoqKlBaWgoAuOeee3D8+PGg\nZB9v2KLxAtMf25lKEHi8+X4XWufO3iem3EHrm4x5i2a8F0KpoOoZiT3zLZHjzhmoFFKkpcixrUaL\nAjV92kwEG5J1LBYLjEaj+GeDwYA//elPG3GoiLl53x049d+/WnZNqq7B93/8v+Kfc3PSxRHADdV6\n5NCMdxKj5lsi5VIgNUWGMl0OtUQmqA1J9olWlhAEhvsf+e6Kyd647RA+UV+Cxrn9VY00453EKC/n\nx9DwBASeX5gSSS2RSWFDkr1er4fJtLBfq8lkgsFgWPK8H/zzEfHra3Zcj207dm5EOGEZcUyLHTPt\nnWZ0/OE/V3xuOfsAT37pSxGMjpDQzLdEzri8SJFLkJmqwK4mA7VExpGWlha0tLSs+XU25AKt3+9H\ndXU13nzzTRQVFWHbtm0xf4HW5fHhfPdwoO7eaYbZHhzbwG//bsXvlcrkONc/tuI6IZGy0BI5iRSF\nDGlKGaqpJTKhxNQFWrlcjn/7t3/DLbfcAp7n8cADD8TcxVmeF9AzuDBIrKtvJGjGe5pKiS3VRWJb\n5G2/fxJ+zgcAyMoJ7FQ1NRHYqYpaz0g0ub2BlkgJ5loiczOw44Yq+r0kQTasLWTfvn3Yt2/fRr18\nWGyjUzjbaRIHic26fOKaVCpB7aaFQWJVpQVBb5Zb938KJ159GVk5arzVGtipamdTBaYmnNh3x6ci\n/nchyWu5lsi915QhPZVaIsnKEnrzkhmXF+e6LGjtNKG1wwzb2FTQui4/W9zAY0t1EdJTL1/HjNed\nqkh8Y4xh1DGLsUVTIjeXF1BLZJKinaoQmPHe1T+C1g4TWjvN6BkYgbDor5eRloKGGr1YminMoxnv\nJDZNu7wwzU2JVCkCUyIrizXUEkliq2YfKYwxWOwTYt39fPfwkhnv9RWFYnKvKMmHTEp1TBJ7fByP\nwWGn2BKZn63C/k9uQooyrt+iJIbE3W/S5Ix7bvu9wNn7qGMmaN1YqBbr7vWVRUilW7tJDAq0RE5h\n1uWBQi5BpkqB5kZ9Um2cTSIr5pM9x/Ho7LOJdfde02jwjPcM1dx8dyOaag3IU2dEL1hCVrC4JVKl\nkCE1RYZaowbFWiO1RJKIiLlkzxjDkNUp1t0/6hmG17cw410hl6GuYn7GuwFlBprxTmKT28dhyOIA\nAKgUgSmRn7qhCnJqiSRREBPJ3jnlmrtb1YS2Tgsck8Ez3kv1GrHuflWlDiollWZI7OF5AUM2Jzxe\nP1RyKdQZStxCLZEkRkQ12f/o1XfR2mlGv3k86HF1VppYd2+sMUCdTYPESOxhjGHUOdcSOTcl8poq\nLbQaaokksSeqrZeltwZm46Qo5Kiv0ol195IimvFOYtOMyweT1Qm5XAKVXIaSwixUF+dSKZFETFy2\nXt59SxOaag2oLacZ7yQ2cX4BA8MOCH4eKqUMudkq3EEtkSQOJdRNVYSslcAYhkemMD3jhlIR2Di7\nsVILdWZqtEMjBECcntkTEm2MMUzOeGCxT0IpD2ycXV2sQUkhtUSSxELJniQdj49Dv9kBKQCVUoai\nvHR86oZKaokkCY2SPUl4PC9gyDoBr49DikKGnHQF9m4rQwa1RJIkQsmeJBzGGMacsxhxzkClkCFN\nKcUnqgpQqEmn0gxJWpTsSUKYdfswZHVCJgmUZkq0WbhhcxW1RBIyh5I9iUt+XsCgxQHOzyNFIUVe\npgq3X7sJqhT6lSZkOfTOIHGBMQb7+AwcE7NzG3gocF19EXKzqSWSkFBQsicxa2rWC5NtAsq5DTzK\n9TnY1ainDTwICQMlexIzvJwfQ8MTgQ08FFJo1Wk4cF05lApZtEMjJO5RsidR8/G7VbNSaQMPQjYK\nJXsSMYwxTE57YBmZQIpCFrhb1Uh3qxISCZTsyYby+vwYsDgAxpCaIoMuN5028CAkCijZk3XF8wLM\n9km4PD6kyKXIyVDilmtKkJ6qjHZohCQ1SvZkTRhjcEy6YBufDmzgoZChYVMeivIyqTRDSAyhZE9W\nzeXhMDjsEO9WNeRn4Lq6SsioNENIzKJkT67IzwswWZ3w+vxIUUihTlfi1u1lSFPRIDFC4gUle7JE\nYG9VF8Yc0+LdqttqtChQ096qhMQrSvYEADA9v7eqTIJUhQylhVm4cUsV3a1KSIKgZJ+kOL+AQasT\nPBcozeRlqbCf9lYlJGHROztJCHODxCYmZ6FUSJGuUuD6eh1ys2iQGCHJgJJ9ApueGySmkEmQqpSh\n3JCD8kY9tUQSkoQo2ScQL+fHgMUJCAJUShm0mjQcuL4CCjm1RBKS7CjZxzFBmBsk5vIgRS5BZqoC\nN11tRHZ6SrRDI4TEGEr2cWS5QWI1xbko1tIgMULI5VGyj3Eenx+DFgcABpVChqI8GiRGCFm9sDPG\n4cOHYTAY0NTUhKamJpw8eVJcO3LkCCorK1FTU4PXX399XQJNFjwvYHDYic4+O/qHRuGemcUt15Tg\nUzdU4dYd5WisLKRETwhZtbDP7CUSCR599FE8+uijQY93dHTglVdeQUdHBywWC26++WZ0d3dDKqUE\ntRzGGMYnXbCPTUOlkEKlDAwS0+dnRTs0QkgCWVMZhzG25LHjx4/j0KFDUCgUKC0tRUVFBc6cOYNr\nr712LYdKKC4Ph0GLAzJpYJCYMT8T199Ag8QIIRtnTcn++eefx09+8hNs3boV//RP/4ScnBwMDw8H\nJXaDwQCLxbLmQOOZnxcwNOyEb+5uVXVGCm7bUYbUFBokRgiJjMsm+927d8Nmsy15/Omnn8bDDz+M\nb3/72wCAJ554Al/72tdw7NixZV9npU6RH/zzEfHra3Zcj207doYceCybHyQ2Oj4NlVKKdJUc22tp\nkBghZPVaWlrQ0tKy5teRsOVqMas0MDCAO+64A+fPn8fRo0cBAI8//jgAYO/evXjqqaewffv24ANL\nJPhoaGKth44Z84PEZFIgVSlDWWE2qkpyaZAYIWRdSSSSZUvoVxJ2GcdqtUKn0wEAXnvtNWzevBkA\nsH//fnz2s5/Fo48+CovFgp6eHmzbti3cw8Qszi9gYNgJwe9HilyKvJxUGiRGCIlZYWemxx57DG1t\nbZBIJCgrK8MPf/hDAEBdXR0OHjyIuro6yOVyvPDCCwlxw4/AGGxjM5icCgwSy0iRY+dmGiRGCIkP\n61LGCevAcVDGmZ71YsjqhEImQVqKHOWGHGwqUlNphhASNREv4yQiL+fH4LATAs9DpZBDq07FnTsr\naZAYISTuJXWyFwQGy8gkZlxecZDYnzXRIDFCSOJJqmTPGMPEtAfDiwaJ1RZrUKzNTojrCoQQspKE\nT/Zenx/9lnFIAKjkMhTl0yAxQkjySbhkz/MCzLZJzHq8SFHIkJOuxN5rSpGeqox2aIQQEjVxn+wZ\nY3BMumAbm0aKQoo0pQwN5Xkoysuk0gwhhMyJy2Tv8gYGiUklDCqFHMUFmbiOBokRQsiK4iLZ87yA\nQasTXl/gblVNhhK3bi9DmooGiRFCSChiMtkzxjDqmMGocxYqhRRpKjmuqS5AoSYj2qERQkhciplk\nP+PyYdDqgFwKpCrlKC3Mxo1b9JBKqe5OCCFrFdVkf2lwDJzfjxSFDPlZKvz5J8tpkBghhGyAqM7G\nGZt00SAxQghZhXBn40Q12Ufp0IQQErfCzZ3Uq0gIIUmAkj0hhCQBSvaEEJIEKNkTQkgSoGRPCCFJ\ngJI9IYQkAUr2hBCSBCjZE0JIEqBkTwghSYCSPSGEJAFK9oQQkgQo2RNCSBKgZE8IIUmAkn2YWlpa\noh3CmlD80UXxR088x74WlOzDFO+/MBR/dFH80RPPsa8FJXtCCEkClOwJISQJRG2nqsbGRrS3t0fj\n0IQQErduvPHGsEpRUUv2hBBCIofKOIQQkgQo2RNCSBKIWLJ3OBzYvXs3qqqqsGfPHkxMTCx5jslk\nwq5du3DVVVehvr4e//qv/xqp8FZ06tQp1NTUoLKyEs8+++yyz/mbv/kbVFZWoqGhAa2trRGO8PKu\nFP/PfvYzNDQ0YMuWLbjuuutw7ty5KES5slB+/gDwwQcfQC6X41e/+lUEo7u8UGJvaWlBU1MT6uvr\n0dzcHNkAr+BK8Y+NjWHv3r1obGxEfX09fvzjH0c+yBXcf//90Gq12Lx584rPieX37ZXiD+t9yyLk\n61//Onv22WcZY4wdPXqUPfbYY0ueY7VaWWtrK2OMsenpaVZVVcU6OjoiFeISfr+flZeXs/7+fubz\n+VhDQ8OSeH7zm9+wffv2McYYe//999n27dujEeqyQon/3XffZRMTE4wxxk6ePBl38c8/b9euXey2\n225jv/zlL6MQ6VKhxO50OlldXR0zmUyMMcZGR0ejEeqyQon/ySefZI8//jhjLBC7RqNhHMdFI9wl\n/vjHP7KzZ8+y+vr6Zddj+X3L2JXjD+d9G7Ez+xMnTuC+++4DANx333349a9/veQ5hYWFaGxsBABk\nZGSgtrYWw8PDkQpxiTNnzqCiogKlpaVQKBS45557cPz48aDnLP57bd++HRMTE7Db7dEId4lQ4t+x\nYweys7MBBOI3m83RCHVZocQPAM8//zzuvvtu5OfnRyHK5YUS+3/913/hrrvugsFgAADk5eVFI9Rl\nhRK/TqfD1NQUAGBqagq5ubmQy+XRCHeJnTt3Qq1Wr7gey+9b4Mrxh/O+jViyt9vt0Gq1AACtVnvF\nH+zAwABaW1uxffv2SIS3LIvFAqPRKP7ZYDDAYrFc8TmxkjBDiX+xY8eO4dZbb41EaCEJ9ed//Phx\nPPzwwwAAiUQS0RhXEkrsPT09cDgc2LVrF7Zu3Yqf/vSnkQ5zRaHE/9BDD+HChQsoKipCQ0MDnnvu\nuUiHGbZYft+uVqjv23X9Z3j37t2w2WxLHn/66aeD/iyRSC77ppyZmcHdd9+N5557DhkZGesZ4qqE\nmjjYx7pXYyXhrCaOP/zhD/jRj36Ed955ZwMjWp1Q4n/kkUdw9OhRSCQSMMaW/L+IllBi5zgOZ8+e\nxZtvvgmXy4UdO3bg2muvRWVlZQQivLxQ4n/mmWfQ2NiIlpYW9Pb2Yvfu3Whvb0dmZmYEIly7WH3f\nrsZq3rfrmux///vfr7im1Wphs9lQWFgIq9WKgoKCZZ/HcRzuuusufO5zn8OBAwfWM7xV0+v1MJlM\n4p9NJpP4kXul55jNZuj1+ojFeDmhxA8A586dw0MPPYRTp05d9qNjpIUS/4cffoh77rkHQOCC4cmT\nJ6FQKLB///6IxvpxocRuNBqRl5eH1NRUpKam4oYbbkB7e3tMJPtQ4n/33XfxrW99CwBQXl6OsrIy\ndHV1YevWrRGNNRyx/L4N1arft+t2ReEKvv71r7OjR48yxhg7cuTIshdoBUFg9957L3vkkUciFdZl\ncRzHNm3axPr7+5nX673iBdr33nsvpi70hBL/4OAgKy8vZ++9916UolxZKPEv9pd/+Zfs1VdfjWCE\nKwsl9s7OTnbTTTcxv9/PZmdnWX19Pbtw4UKUIg4WSvxf/epX2eHDhxljjNlsNqbX69n4+Hg0wl1W\nf39/SBdoY+19O+9y8Yfzvo1Ysh8fH2c33XQTq6ysZLt372ZOp5MxxpjFYmG33norY4yxt956i0kk\nEtbQ0MAaGxtZY2MjO3nyZKRCXNZvf/tbVlVVxcrLy9kzzzzDGGPsxRdfZC+++KL4nC996UusvLyc\nbdmyhX344YfRCnVZV4r/gQceYBqNRvx5X3PNNdEMd4lQfv7zYinZMxZa7P/4j//I6urqWH19PXvu\nueeiFeqyrhT/6Ogou/3229mWLVtYfX09+9nPfhbNcIPcc889TKfTMYVCwQwGAzt27FhcvW+vFH84\n71sal0AIIUmA7qAlhJAkQMmeEEKSACV7QghJApTsCSEkCVCyJ4SQJEDJnhBCkgAle0IISQKU7Akh\nJAn8f1DV3LgUNhRHAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xc824510>"
       ]
      }
     ],
     "prompt_number": 126
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}