{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Double Digest Problem\n",
    "\n",
    "In diesem Beispiel betrachten wir das Double Digest Problem. \n",
    "Hierfür werden Restriktionsenzyme verwendet um die DNA in Fragmente unterschiedlicher Länge zu zerschneiden.  \n",
    "Mit Hilfe der Fragmente kann man dann eine physikalische Karte erstellen, welche die Reihenfolge der Fragmente darstellt. Für das Double Digest Problem werden zwei Restriktionsenzyme benötigt, Restriktionsenzym A und Restriktionsenzym B.  \n",
    "\n",
    "Das Double Digest Verfahren funktioniert wie folgt:  \n",
    "\n",
    "1. Erzeuge drei Kopien der DNA\n",
    "2. Wende die Restriktionsenzyme auf die DNA Kopien an:  \n",
    "    1. Wende Restriktionsenzym A auf die erste DNA Kopie an\n",
    "    2. Wende Restriktionsenzym B auf die zweite DNA Kopie an\n",
    "    3. Wende beide Restriktionsenzyme A + B auf die dritte Kopie an\n",
    "3. Bestimme für jede der Kopien die Längen der Fragmente und speichere diese in drei Listen, $fragments_a$, $fragments_b$ und $fragments_{ab}$\n",
    "4. Erstelle alle Permutationen für die Listen $fragments_a$ und $fragments_b$\n",
    "5. Berechne für alle möglichen Kombinationen der Permutationen in $fragments_a$ und $fragments_b$\n",
    "    1. die Positionen der Restriktionsstellen in $fragments_a$ und speichere diese in $pos_a$\n",
    "    2. die Positionen der Restriktionsstellen in $fragments_b$ und speichere diese in $pos_b$\n",
    "    3. die Vereinigung der Positionen von $pos_a$ und $pos_b$\n",
    "    4. die Distanzen zwischen allen Positionen in der Vereinigungsmenge\n",
    "    5. Vergleiche die Distanzen mit den Fragmentlängen aus $fragments_{ab}$. Falls die Mengen gleich sind ist eine gültige Lösung gefunden, welche wir in unserer Lösungsmenge speichern\n",
    "    \n",
    "Der Algorithmus findet alle möglichen Lösungen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code Double Digest Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline  \n",
    "import scipy as sp\n",
    "import itertools\n",
    "import pylab as pl\n",
    "import matplotlib.patches as patches\n",
    "\n",
    "'''\n",
    "Finde alle möglichen Lösung mittels Double Digest\n",
    "'''\n",
    "def double_digest(fragments_a, fragments_b, fragments_ab):\n",
    "    fragments_ab = sp.sort(fragments_ab)\n",
    "    #erstelle alle möglichen Permutationen von fragments_a und fragments_b\n",
    "    permsA = itertools.permutations(fragments_a)\n",
    "    permsB = itertools.permutations(fragments_b)\n",
    "    \n",
    "    results = []\n",
    "    for permA, permB in itertools.product(permsA,permsB):\n",
    "        #berechne Positionen für Permutation A\n",
    "        posA = sp.cumsum(permA)\n",
    "        #berechne Positionen für Permutation B\n",
    "        posB = sp.cumsum(permB)\n",
    "        #Bilde Schnittmenge zwischen posA und posB\n",
    "        pos = sp.union1d(posA,posB)\n",
    "        #Berechne inverse CumSum bzw. Distanzen zwischen allen Positionen in \"pos\"\n",
    "        dists = pos.copy()\n",
    "        dists[1:] = sp.diff(pos)\n",
    "\n",
    "        if dists.shape[0]!=fragments_ab.shape[0]:\n",
    "            continue\n",
    "        #Falls die Permutierte Distanzliste gleich der Menge der doppelt verdauten Liste enspricht\n",
    "        #ist eine gültige Lösung gefunden\n",
    "        if sp.all(sp.equal(sp.sort(dists),fragments_ab)):\n",
    "            #if not {\"permA\":permA,\"permB\":permB,\"positions\":pos} in results:\n",
    "            results.append({\"posA\":posA,\"posB\":posB,\"positions\":pos})\n",
    "    return results\n",
    "\n",
    "'''\n",
    "Plotte Marker\n",
    "'''\n",
    "def plot_markers(positions=None,level=1,label=None,color=\"b\",fig_obj=None):\n",
    "    pl.axhline(level,color=color,label=label)\n",
    "    for pos in positions:\n",
    "        fig_obj.add_patch(patches.Rectangle((pos, level-0.25),2,0.5,color=\"k\"))\n",
    "\n",
    "'''\n",
    "Plotte Double Digest Lösung\n",
    "'''\n",
    "def plot_solutions(results=None):\n",
    "    for i,res in enumerate(results):\n",
    "        pl.figure(figsize=(15,5))\n",
    "        fig = pl.subplot(1,1,1)\n",
    "        \n",
    "        plot_markers(positions=res['positions'],level=1,\n",
    "                     label=\"Enzym A+B\",color=\"b\",\n",
    "                     fig_obj=fig)\n",
    "        \n",
    "        plot_markers(positions=res['posB'],level=2,\n",
    "                     label=\"Enzym B\",color=\"r\",\n",
    "                     fig_obj=fig)\n",
    "        \n",
    "        plot_markers(positions=res['posA'],level=3,\n",
    "                     label=\"Enzym A\",color=\"y\",\n",
    "                     fig_obj=fig)\n",
    "        \n",
    "        pl.ylim(0,4)\n",
    "        pl.xlim(0,res['positions'].max())\n",
    "        pl.xlabel(\"Position (kb)\", fontsize=14)\n",
    "        pl.yticks([1,2,3],[\"Enzym A+B\",\"Enzym B\",\"Enzym A\"],fontsize=14)\n",
    "        pl.title(\"Loesung %d\"%(i+1),fontsize=14)\n",
    "        pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beispiel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Positionen Enzym A:\t[ 20 220 620 720]\n",
      "Positionen Enzym B:\t[ 50 300 450 670]\n",
      "Loesung 1:\t\t[ 20  50 220 300 450 620 670 720]\n",
      "\n",
      "Positionen Enzym A:\t[100 500 700 720]\n",
      "Positionen Enzym B:\t[ 50 270 420 670]\n",
      "Loesung 2:\t\t[ 50 100 270 420 500 670 700 720]\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAFgCAYAAABNIolGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+4rXVdJ/z3RxAsbRRDOf5IT06YWpOo5zEdUyGNsPxR\njTOpQwiilnb16Mwzj6P4lORcMFGO5tRMpSWooI0/0kif5MnUYbLUDiBGEooMCCqIooiC/PLz/LHu\n7Vlns8/ZG84+Z33P2a/Xda1rr3Xf3/u+P2ud77r3Xu/zvb+rujsAAAAAI7vTogsAAAAAWI0AAwAA\nABieAAMAAAAYngADAAAAGJ4AAwAAABieAAMAAAAYngADAAAAGJ4AAwDY46rqtKp636LrWJSqen1V\nba2qb1fVpYuuBwD2BgIMAIA9705J3pzkLYsuBAD2FgIMAGA4VfWAqnpPVV033f6squ6/rM3Tquqc\naRTD/66qk6rqgLn1v1BVn6qqG6rqmqr6n1V1yLTuxKq6YNn+jq2qb849PrGqLqiqZ1XV56Y63ltV\nB8+12b+qXldVX5tur6uqP6iqj+zs+XX3r3X37yX5zK69UgCwcQgwAIChVNWdkvx5kkOSHDHd7pvk\nvVVVU5ufTnJGkt9P8iNJnpfkmUlOntZvSvKnmY1yeGiSJyR56x0oZ3OSX0zy80mOTPKIJCfNrf8P\nSY5N8vwkj8nsb6vn3IHjAACr2H/RBQAALPOkJD+W5J9396VJUlXPSXLxtO6DSV6Z5He6+9Rpm89V\n1X9McnpV/d+ZBR53TvKu7r5sarPdiIs12j/Jsd197VTHG5IcN7f+JUlO6e53T+tfmuSoO3AcAGAV\nRmAAAKN5aJIvLoUXSdLdlyT5YpKHTYseleSVVfXNpVuStyW5a5JNSc7PLOi4oKreXVUvqqp73YFa\nLlsKLyZfTHLvJKmqu0/H+sRcnT3/GABYPwIMAGBv0tPPOyX5zSSHzd1+LMmhSa7u7lszu+TjyCSf\nSnJ8ks9W1cOn7b+TpJbt+84rHO/mFY7v7ycAWAC/gAGA0VyY5L5VtXlpQVU9KLPLQj49LTo3yUO6\n++IVbrcks9EQ3f133f2bSf6PzEZP/OK0/dVJDlmaU2Ny2O0pchqZceW076U6a/4xALB+zIEBACzK\nP6uq5aHB1zO79ONTSc6oqpdMy38vs9DiQ9PjVyd5X1VdluQdSW5J8qNJHt3dL6uqxyR5cpKzklyV\n2eSbP5BtAchHktwzyQlV9adJDs9sEtDb6/VJXlZVn5n2/ctJ7pPkSzvbqKp+KMndMgtlDph7HT7d\n3TfdgToAYJ9nBAYAsCiPT3LesttrpnkknpHZKIkPT7crk/zctC7dfVaSn83sG0o+Md1enuTz076v\nTfK4JO9L8tkk/yXJf+ru06ftL0zyoiQvzCws+alM32ByO70ms283OTXJx6Zl70ny7VW2++Pp+f67\nzAKPped/3ztQAwBsCDX9HQAAwDqoqvOS/E13/9qiawGAfYlLSAAA7qCqemCSn07yPzObBPQFmU0m\n+oJF1gUA+yIBBgDAHfedJMck+Z3MLs39dJKndPfWhVYFAPsgl5AAAAAAwzOJJwAAADA8l5AM4OCD\nD+7NmzcvugwAAADY484555yvdPe9VmsnwBjA5s2bs3WrS2UBAADYeKrqsrW0cwkJAAAAMDwBBgAA\nADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAA\nMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAw\nPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8\nAQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwB\nBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEG\nAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYA\nAAAwPAEGAAAAMDwBxgayadOmVFWqKps2bVp0OcDgls4ZzhcAwHrxmYRdIcDYQK666qoV7wOsZOk8\n4XwBAKwXn0nYFQIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAA\nYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAsLMKrqtKrqFW4fW1RN66Gqzq+q\nW6rqwYuuBQAAAPYVix6B8cEk91l2+5mFVrQLqurRSe6d5C1Jjl9wOQAAALDPWHSAcWN3X7nsds3S\nymlExgur6p1V9a2quqSqjp5bv6NRHMdW1TFV9dWqOnD+gFV1RlWdOd0/saouqKrnVtWl0zFOraoD\nqurFVXX5tI/XVtVaXqvjk7wtyalJjqmq/dfnZQIAAICNbW/4gP0bSV6e5BWZBQRvqqqzu/vzSV4y\nrVvyvCSvTLI1yeeSvD7JM5K8I0mq6u5Jfj7Js+e22Ty1eWqS+yV5d2YjQb6U5MgkD5m2/+i0bkVV\nddckz0pyeJJPJrlh2ud7V3uC119/Uc477/DVmu2y171u+8d74pjA3mv+nOF8AQCsB59J2BWLHoFx\nVFV9c9ntlGVt3trdp3f3xUl+PcktSZ6QJN197dLIjSSHZhZ2HN3dF3T3DUnOyCzUWPKcJN9I8v65\nZfslOW7a5qwkH0jyqCS/3N0Xdvd7MgsvjljlufzrJJd393nd3UlOT/L8HTWeRpZsraqtN9988yq7\nBgAAgI1t0SMwzk7ywmXLvr7s8aeW7nT3LVV1dWbzTHxXVW3ObHTEq6fAYckbk5xbVffv7isyCzPe\n3N23zLX5fHdfO/f4qiSf6e6bli3b7pgreH6St849fmuSE6rqvt39xeWNu/sNSd6QJFu2bOlHPOIj\nq+x+1z3ykbWsht1/TGDvNX/OcL4AANaDzySsrFZvksWPwLi+uy9edvvKsjbLhyd05uquqrslOTPJ\nWd198nYNu89Pcm6SY6vqR5NsSfKmNex/p8dcrqoekuRxSU6evoHkliQXZhrdsaPtAAAAgLVZ9AiM\nXTJNrHlGkuuy48s13pjkZUkOTvLR7r5oN5RyfJKPr1DDv0ryvKo6ebqsBAAAALgDFh1gHFhVm5Yt\nu7W7r17j9q9K8tgkT05yUNV3h51cO82BkSRvT/LaJC9K8iu7WO9tVNWdkxyT5KTuvmDZuq9mNi/H\nEUk+tN7HBgAAgI1i0ZeQPDmzb/uYv513O7Z/YpJ7JTl/2T5+calBd1+X2beI3Dj9XG9Pm2q4zTeU\ndPeXMpsAdIeTeQIAAACrq41wZUNV/WWSK7r7BYuuZSVbtmzprVu37vbjzI1QSZJshH974I6bP2c4\nXwAA68FnElZSVed095bV2i36EpLdqqoOSvL4JEcmefiCywEAAADuoH06wMjscpR7Jjlh+fwUAAAA\nwN5jnw4wunvzomsAAAAAdt2iJ/EEAAAAWJUAAwAAABieAAMAAAAYngADAAAAGJ4AAwAAABieAAMA\nAAAYngADAAAAGJ4AAwAAABieAGMDOeSQQ1a8D7CSpfOE8wUAsF58JmFX7L/oAthzrrzyykWXAOxF\nnDMAgPXm7wt2hREYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAA\nwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA\n8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDw\nBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAE\nGAAAAMDwBBgAAADA8AQYAAAAwPAEGAxl06ZNqaps2rRp0aUA7DWWzp3On7Ay7xFgRM5Nt19196Jr\n2PC2bNnSW7duXXQZQ6iq797XNwHWZv7cmTh/wnLeI8CInJu2qapzunvLau2MwAAAAACGJ8AAAAAA\nhifAAAAAAIYnwAAAAACGJ8AAAAAAhifAAAAAAIYnwAAAAACGJ8AAAAAAhifAAAAAAIYnwAAAAACG\nJ8AAAAAAhifAAAAAAIa3sACjqk6rql7h9rFF1XRHVdXhy57DDVX1qao6ZtG1AQAAwL5g/wUf/4NJ\nfmnZspsWUcg6+ZEk1yT53iTPTHJaVX2mu/e6UAYAAABGsuhLSG7s7iuX3a5ZWjmNZnhhVb2zqr5V\nVZdU1dFz63c0iuPYqjqmqr5aVQfOH7CqzqiqM6f7J1bVBVX13Kq6dDrGqVV1QFW9uKoun/bx2qpa\ny2v15ek5XNLdv51ZmPGIdXqtAAAAYMNa9AiMtfiNJC9P8ookxyd5U1Wd3d2fT/KSad2S5yV5ZZKt\nST6X5PVJnpHkHUlSVXdP8vNJnj23zeapzVOT3C/Ju5PcJ8mXkhyZ5CHT9h+d1q1qCjuenuQeSf5+\n1Q0uuig5/PC17Hqf9+H5B14TgDX58PIFzp+wHe8RYETOTbffokdgHFVV31x2O2VZm7d29+ndfXGS\nX09yS5InJEl3X7s0ciPJoZmFHUd39wXdfUOSMzILNZY8J8k3krx/btl+SY6btjkryQeSPCrJL3f3\nhd39nszCiyPW8HwurapvJrkxybuSvKy7t67UcBpZsrWqtt58881r2DUAAABsXIsegXF2khcuW/b1\nZY8/tXSnu2+pqquT3Hu+QVVtzmx0xKunwGHJG5OcW1X37+4rMgsz3tzdt8y1+Xx3Xzv3+Kokn+nu\nm5Yt2+6YO3BEkq8lOTDJo5P8XlVd191vXN6wu9+Q5A1JsmXLls5HPrKG3e/7jqj67v32mgCsyfy5\nM3H+hOW8R4AROTfNWfZa7MiiA4zrp5EVO7N8eEJnbuRIVd0tyZlJzuruk7dr2H1+VZ2b5Niqem+S\nLUmOzvZW2v9Ky/Zbpc4k+d/d/ZXp/j9W1Y9nNirkNgEGAAAAsHaLDjB2yTTXxBlJrkvy/B00e2OS\nlyU5OMlHu/uiPVRektya2TeSAAAAALtg0QHGgVW1admyW7v76jVu/6okj03y5CQH1bZhJ9dOc2Ak\nyduTvDbJi5L8yi7Wu5p7V9X+2XYJyS9lNhcGAAAAsAsWHWA8ObNv+5j3hST3X+P2T0xyryTnL1t+\nXJLTkqS7r6uqdyR5ZqZvI9mN/nH6eUuSy5P8UZITd/MxAQAAYJ9X3b3oGna7qvrLJFd09wsWXctK\ntmzZ0lu3rvhlJRtOzU/iuQH6JsB6qOWTgDl/wna8R4AROTdtU1XndPeW1dotegTGblVVByV5fJIj\nkzx8weUAAAAAd9A+HWAkOS/JPZOc0N0XLLoYAAAA4I7ZpwOM7t686BoAAACAXXenRRcAAAAAsBoB\nBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQZDOeSQQ7b7\nCcDq5s+Zzp9wW94jwIicm26//RddAMy78sorF10CwF7HuRN2znsEGJFz0+1nBAYAAAAwPAEGAAAA\nMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAw\nPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8\nAQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwB\nBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AcYGsmnTplRV\nqiqbNm1adDkrWqpx1PpgI/F+3HvsDed3WCTvERiH9+M2Xovbr7p70TVseFu2bOmtW7fu9uNU1XaP\nR/y3n69xxPpgI/F+3HvsDed3WCTvERiH9+M2Xottquqc7t6yWjsjMAAAAIDhCTAAAACA4QkwAAAA\ngOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA\n4QkwAAAAgOGtGmBU1WlV1SvcPrYnCtxdqur8qrqlqh68Dvu6dO51ubWqvlBVf1hV37cetQIAAMBG\nt9YRGB9Mcp9lt5/ZXUXtblX16CT3TvKWJMevof1pVXXiKs1endnr8oAkz83s9fntXasUAAAASNYe\nYNzY3Vcuu12ztHIaefDCqnpnVX2rqi6pqqPn1u9oFMexVXVMVX21qg6cP2BVnVFVZ073T6yqC6rq\nudNoh29V1alVdUBVvbiqLp/28dqqWstzOj7J25KcmuSYqtp/ja/Dzlw3vS5f6O4PJnlHkkeuw34B\nAABgw1uPD+5LfiPJy5O8IrOA4E1VdXZ3fz7JS6Z1S56X5JVJtib5XJLXJ3lGZh/6U1V3T/LzSZ49\nt83mqc1Tk9wvybszG/HwpSRHJnnItP1Hp3Urqqq7JnlWksOTfDLJDdM+33sHn/dKx3hAkp9O8uG1\ntL/oouTww9fr6DuzfTl75pi317Yax6wPNhLvx73H3nB+h0XyHoFxeD9u47W4vdY6AuOoqvrmstsp\ny9q8tbtP7+6Lk/x6kluSPCFJuvvapZEbSQ7NLOw4ursv6O4bkpyRWaix5DlJvpHk/XPL9kty3LTN\nWUk+kORRSX65uy/s7vdkFl4cscpz+ddJLu/u87q7k5ye5PlrfB125qTpdbkhyWVJvpbkhB01nkas\nbK2qrTfffPM6HB4AAAD2XTX7DL+TBlWnZTavwwuXrfp6d39latNJntPdb5/b7rIkr+/u184t25zk\nE0l+t7tPnlv+8CTnJnlgd19RVX+f5EPd/R+n9ScmeXZ3//DcNv8tyWHd/bi5Zf9jek7/ZifP52+S\n/EV3nzI9PjTJhUke0N1fnJadkO3DhwOTdJKb5pY9pbv/19T+0iRvT/InSSrJDyQ5OcmNSY7o7u/s\nqJ4k2bJlS2/dunVnTdZFVW33eLV/+0WYr3HE+mAj8X7ce+wN53dYJO8RGIf34zZei22q6pzu3rJa\nu7VeQnL9NLJiZ5YPI+jMjfCoqrslOTPJWfPhRZJ09/lVdW6SY6vqvUm2JDk621tp/yst229HBVbV\nQ5I8Lsljq+qkuVX7JTkuydKyP8x0OcvklCRfSPJf55Z9Ydnuvzr3Gn22ql6a5O8yGxHy1zuqCQAA\nAFjdes6BsUPTxJpnJLkuO75c441JXpbk4CQf7e6LdkMpxyf5+Ao1/Kskz6uqk3vmmiTzk5Rel+Sa\nNYQ4826dfn7vrhQMAAAArD3AOLCqNi1bdmt3X73G7V+V5LFJnpzkoLmhMtdOc2Aks0swXpvkRUl+\nZY37XbOqunOSY5Kc1N0XLFv31czm5TgiyYfu4CG+b3qNli4h+e0kVyf52ztcNAAAAJBk7ZN4Pjmz\nb/uYv513O47zxCT3SnL+sn384lKD7r4us8s2bsz2l2+sl6dNNdzmG0q6+0uZTQC6K5N5/kZmz+mL\nSd6X5FtJjuzur+7CPgEAAICsYRLPPamq/jLJFd39gkXXsieZxHMbkwbCOLwf9x57w/kdFsl7BMbh\n/biN12Kb9Z7Ec7eqqoOSPD7JkUkevuByAAAAgMEMEWBkdjnKPZOcsHx+CgAAAIAhAozu3rzoGgAA\nAIBxrXUSTwAAAICFEWAAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AA\nAAAAwxNgbCCHHHLIivdHslTXqPXBRuL9uPfYG87vsEjeIzAO78dtvBa3X3X3omvY8LZs2dJbt25d\ndBkAAACwx1XVOd29ZbV2RmAAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADD\nE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMT\nYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNg\nAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AA\nAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAA\nAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAA\nAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAw6vuXnQNG15VXZfkokXXwXAOTvKV\nRRfBcPQLVqJfsBL9gpXoF6xEv2Ale7JfPLC777Vao/33RCWs6qLu3rLoIhhLVW3VL1hOv2Al+gUr\n0S9YiX7BSvQLVjJiv3AJCQAAADA8AQYAAAAwPAHGGN6w6AIYkn7BSvQLVqJfsBL9gpXoF6xEv2Al\nw/ULk3gCAAAAwzMCAwAAABieAAMAAAAYngBjgarqqKq6qKourqqXL7oe9pyqelNVfbmqLphbds+q\n+quq+uz086BpeVXVf536yaeq6pGLq5zdqap+oKo+XFWfrqp/rKqXTMv1jQ2squ5SVZ+oqvOnfvGb\n0/IfrKqPT//+/6OqDpiWHzg9vnhav3mR9bN7VdV+VXVeVb1veqxfbHBVdWlV/UNVfbKqtk7L/B7Z\n4KrqHlX1rqr6p6q6sKoeq19QVT88nSuWbt+oqpeO3DcEGAtSVfsl+W9JnpLkYUmeXVUPW2xV7EGn\nJTlq2bKXJ/nr7j40yV9Pj5NZHzl0ur0wyR/soRrZ825J8n9198OSPCbJr07nBX1jY7sxyU9298OT\nHJbkqKp6TJJTkryuu38oydeSHD+1Pz7J16blr5vase96SZIL5x7rFyTJEd19WHdvmR77PcLrk3yg\nux+S5OGZnTf0iw2uuy+azhWHJXlUkuuTvCcD9w0BxuI8OsnF3X1Jd9+U5E+TPGPBNbGHdPfZSa5Z\ntvgZSd483X9zkp+bW/6WnvlYkntU1X32TKXsSd39pe4+d7p/XWZ/XNwv+saGNv37fnN6eOfp1kl+\nMsm7puXL+8VSf3lXkidVVe2hctmDqur+SX42yR9Pjyv6BSvze2QDq6q7J3lCkj9Jku6+qbu/Hv2C\n7T0pyee6+7IM3DcEGItzvySXzz2+YlrGxnVId39pun9lkkOm+/rKBjQN735Eko9H39jwpssEPpnk\ny0n+Ksnnkny9u2+Zmsz/23+3X0zrr03y/Xu2YvaQ303ysiTfmR5/f/QLZgHn/1dV51TVC6dlfo9s\nbD+Y5Ookp06XnP1xVd01+gXbe1aSt0/3h+0bAgwYUM++39h3HG9QVXW3JO9O8tLu/sb8On1jY+ru\nW6fhnffPbATfQxZcEgtWVU9N8uXuPmfRtTCcn+juR2Y21PtXq+oJ8yv9HtmQ9k/yyCR/0N2PSPKt\nbLskIIl+sdFN8yU9Pck7l68brW8IMBbnC0l+YO7x/adlbFxXLQ3Bmn5+eVqur2wgVXXnzMKLM7r7\nz6bF+gZJkmnI74eTPDazYZv7T6vm/+2/2y+m9XdP8tU9XCq73+OSPL2qLs3sMtSfzOwad/1ig+vu\nL0w/v5zZteyPjt8jG90VSa7o7o9Pj9+VWaChX7DkKUnO7e6rpsfD9g0BxuL8fZJDp9nCD8hsyM6Z\nC66JxTozyXOn+89N8udzy4+ZZv19TJJr54Z0sQ+Zrkf/kyQXdvdr51bpGxtYVd2rqu4x3f+eJD+V\n2fwoH07yzKnZ8n6x1F+emeRD0/+esA/p7ld09/27e3Nmf0N8qLv/bfSLDa2q7lpV37d0P8mRSS6I\n3yMbWndfmeTyqvrhadGTknw6+gXbPDvbLh9JBu4b5XfX4lTVz2R2/ep+Sd7U3SctuCT2kKp6e5LD\nkxyc5Kokr0ry3iTvSPKAJJcl+Tfdfc30ofb3M/vWkuuTHNfdWxdRN7tXVf1Ekv+V5B+y7Zr2EzKb\nB0Pf2KCq6scym0Brv8z+4+Ed3f3qqnpQZv/zfs8k5yU5urtvrKq7JHlrZnOoXJPkWd19yWKqZ0+o\nqsOT/Ifufqp+sbFN//7vmR7un+Rt3X1SVX1//B7Z0KrqsMwm/D0gySVJjsv0OyX6xYY2hZ2fT/Kg\n7r52WjbsOUOAAQAAAAzPJSQAAADA8AQYAAAAwPAEGAAAAMDwBBgAAADA8AQYAAAAwPAEGADAXqmq\nDq+qrqqDV2n3kar6/T1U0xOr6jNVtd/0+MSqumAn7Y+tqm+uss9fraq/WO9aAWBvI8AAAHaLqjpt\nChi6qm6uqkuq6jXTd86vh79Ncp8kX52Ot6Mw4BeSvGKdjrma30lyUnffuo77/OMkj6qqx6/jPgFg\nryPAAAB2pw9mFjI8KMn/k+TFSV6zHjvu7pu6+8ru7lXaXdPd163HMXemqv5lkockecd67re7b0zy\ntiT/53ruFwD2NgIMAGB3unEKGS7v7rclOSPJzy2trKonVNXHq+rbVXVVVb2uqg5Ytv5jVfXNqrq2\nqj5RVT86rfvuJSRVdXiSU5PcdW7Ux4lTu+0uIamqg6rqzVX1taq6oao+WFU/Mrf+2Ol4T6qqC6rq\nW1X14ar6wVWe63OS/HV337CjBlX1gKr6p+n4+88tf9p06cm3p2M9aNmmZyZ5elV97yo1AMA+S4AB\nAOxJNyS5c5JU1f2S/GWS85I8IsnxSZ6d5D9P6/dP8udJ/ibJw5P8eJLfTbLS5Rl/m+SlSa7PbMTH\nfbLjkR6nTft6RpJHT9t8oKq+Z67NgZlddvK8JI9Nco8kf7jKc3t8kq07WllVD03y0ST/b5Jju/uW\nuWO9KsnQek0WAAADW0lEQVRx07H2S/JnVVVzm29Nsv+0HgA2pP1XbwIAsOuq6tGZRilMi16c5ItJ\nXtzd30lyYVW9PMkfVdWvJ7lLZsHBX3T356Zt/mmlfXf3TVV17exuX7mTGg5N8vQkT+zus6dlv5Tk\n80n+bWbzTSSzv5F+tbsvmtq8Jsmbqqp2csnKA6fns9JxfzzJ+5O8rrtPWrZ6/yQv6e6PztVzSZIn\nZXYJTrr7+un5bd7RcwOAfZ0RGADA7nTUdDnGt5P8XZKzk/zatO6hST42hRdL/ibJAUl+qLuvyWy0\nxFlV9f6q+vdV9YBdrOehSb4z1ZIk6e5rk/xDkofNtbtxKbyYfHGq66Cd7Pt7knx7heX3yyyIOGWF\n8CJTPZ+Yq+ey6XgPW9buhukYALAhCTAAgN3p7CSHJfnhJHfp7l/o7i+vYbtOku4+LrPLPc7ObOTE\nRVX107up1vmRFbfsYN3O/nb6SlYOOL6S5GNJnlVVOwpAdjoR6eSeSa5eQzsA2CcJMACA3en67r64\nuy/r7puXrbswyWOqav7vkZ9IclOSpUtG0t3nd/cp3X14ko8kee4OjnVTZvNH7MyFmf398925JKrq\nnyX5F0k+vfrT2anzcttRE0lyY2bhy9eS/FVV3WPZ+jtlNhfHUj0PSHLfqdalZf88s0tqzt3FGgFg\nryXAAAAW5b9n9kH9v1fVQ6vqZ5P8VpLfn+Z8+MGq+q2q+pdV9cCqOiLJj2XHQcOlSe5SVT81fTPJ\nbb6xo7s/m9nEoH9UVY+vqn+R5PQk38jsq0p3xVmZBTC3MX0zydOSXJvbhhi3JPndqnpsVR2W5M1J\n/jHT/BeTxye5ZKofADYkAQYAsBDd/YUkT8nsG0g+meRNSd6e5ISpyfVJHpzknUk+k9kH+zOSnLKD\n/f1tZt8U8vbMLrV42Q4OfVxmc06cOf383iRH7ezrT9fo9CQPnv9K1mX13ZDkqZmFJfMhxo1JTkry\nliQfz+zvs19YNlnos5O8cRfrA4C9Wu14Im0AAG6PqvqtJPfq7uPXcZ8/mtk3tzx4mnAUADYkIzAA\nANbPyUkuqarV5uK4Pe6b5BjhBQAbnREYAAAAwPCMwAAAAACGJ8AAAAAAhifAAAAAAIYnwAAAAACG\nJ8AAAAAAhifAAAAAAIb3/wNTkjq7WkxMhQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10de4b4d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAFgCAYAAABNIolGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu4bWV9H/rvLyCYmEQxKNtLdMdWRZMG1F2iNSpEQzTx\nkouplxAEURLN0+o5p4cqNpF4HjghsV7atEk0EYwg1ks0RBtpjFAaEzUbEEMkKFIQLyCCIghy0V//\nmGPBZLH2Xgv2Zb57r8/neeaz5hzjHWP85lzvHGvu737HO6u7AwAAADCy71l0AQAAAACrEWAAAAAA\nwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAADtdVZ1cVR9cdB2LUFUHVNVp\nVXV5Vd1YVRdV1TFV5XMZAGzFnosuAABgnXlckquS/GqSLyQ5KMlbM/tcdsIC6wKAoUn6AYDhVNVD\nqur9VXXddPuzqnrwsjbPqqpzqurbVfW/q+r4qtprbv0vVtWnp1EO11TV/6yq/aZ1x1XVBcv2d0RV\nXT/3+LiquqCqnl9Vn5/q+EBV7TvXZs+qemNVfX26vbGq/qCqztrSc+vut3X3v+3us7r7ku5+V5I/\nSPJL2/zCAcBuTIABAAxlupTiz5Psl+SQ6fbAJB+oqpra/EySU5P8fpIfTfLiJM/NNIKhqjYkeVeS\ntyd5VJInJ3nH3ShnY5LnJfmFJIcmeUyS4+fW/7skRyR5SZLHZ/bZ6oV34zg/mOTrd2M7AFg3XEIC\nAIzmqUl+PMk/6+5Lk6SqXpjk4mndR5K8JsnvdfdJ0zafr6p/n+SUqvp/Mws87pHkvd192dTmDiMu\n1mjPJEd097VTHW9JcuTc+lckObG73zetf2WSp9+VA1TVYzMLQX7lbtQHAOuGERgAwGgeleTLS+FF\nknT3JUm+nOTR06LHJXlNVV2/dEvyziT3SrIhyfmZBR0XVNX7quplVXW/u1HLZUvhxeTLSe6fJFV1\n7+lYn5yrs+cfr6aqHpnkQ0netBSCAAArE2AAALuSnn5+T5LfTnLg3O3Hkzw8yVXd/Z3MLvk4NMmn\nkxyV5HNVdcC0/XeT1LJ932OF492ywvG3y+enqto/yVlJ3tXdr9oe+wSA3ZkAAwAYzYVJHlhVG5cW\nVNXDMrss5DPTonOT7N/dF69wuzWZjYbo7r/r7t9O8i8zGz3xvGn7q5LstzSnxuTAu1LkNDLjimnf\nS3XW/OMtqapHZxZevKe7/6+7clwAWK/MgQEALMoPVtXy0OAbmV368ekkp1bVK6bl/zmz0OKj0+PX\nJflgVV2W5N1Jbk3yY0kO6u5jqurxSZ6W5IwkV2Y2+eYP5/YA5Kwk901ybFW9K8nBmU0Cele9Ockx\nVfXZad+/luQBSb6ypQ2q6ken53FmkhOmCUeTJN19xd2oAQDWBSMwAIBFeVKS85bdXj/NI/GczEZJ\nnDndrkjy89O6dPcZSX4us28o+eR0e1WSL0z7vjbJE5N8MMnnkvzHJP9fd58ybX9hkpclOTqzsOSn\nM32DyV30+sy+3eSkJB+flr0/ybe3ss0vZzaPxvMyCzrmbwDAFtT0OQAAgO2gqs5L8jfd/W8WXQsA\n7E5cQgIAcDdV1UOT/EyS/5nZJKAvzWwy0Zcusi4A2B0JMAAA7r7vJjk8ye9ldmnuZ5I8o7s3L7Qq\nANgNuYQEAAAAGJ5JPAEAAIDhuYRkAPvuu29v3Lhx0WUAAADATnfOOed8rbvvt1o7AcYANm7cmM2b\nXSoLAADA+lNVl62lnUtIAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4Qkw\nAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAA\nAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAA\nAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAA\ngOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA\n4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDh\nCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMGBQGzZsSFVlw4YNiy4FANhNLH2+\n8BkDWM2I5wsBBgzqyiuvvMNPAIBtNf+5wmcMYGtGPF8IMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAA\nAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAA\ngOEtLMCoqpOrqle4fXxRNW0PVXV+Vd1aVY9YdC0AAACwu1j0CIyPJHnAstvPLrSibVBVByW5f5I/\nTXLUgssBAACA3caiA4ybuvuKZbdrllZOIzKOrqr3VNW3quqSqjpsbv2WRnEcUVWHV9XVVbX3/AGr\n6tSqOn26f1xVXVBVL6qqS6djnFRVe1XVy6vq8mkfb6iqtbxWRyV5Z5KTkhxeVXtun5cJAAAA1rdd\n4R/Yv5XkVUlenVlA8LaqOru7v5DkFdO6JS9O8pokm5N8PsmbkzwnybuTpKruneQXkrxgbpuNU5tn\nJnlQkvdlNhLkK0kOTbL/tP3HpnUrqqp7JXl+koOTfCrJjdM+P7DaE7zhhoty3nkHr9aMdeaNb7z9\nvv4BAGwP858vEp8xgC0b8Xyx6BEYT6+q65fdTlzW5h3dfUp3X5zkN5PcmuTJSdLd1y6N3Ejy8MzC\njsO6+4LuvjHJqZmFGktemOSbST40t2yPJEdO25yR5MNJHpfk17r7wu5+f2bhxSGrPJdfTnJ5d5/X\n3Z3klCQv2VLjaWTJ5qrafMstt6yyawAAAFjfFj0C4+wkRy9b9o1ljz+9dKe7b62qqzKbZ+I2VbUx\ns9ERr5sChyVvTXJuVT24u7+YWZjx9u6+da7NF7r72rnHVyb5bHffvGzZHY65gpckecfc43ckObaq\nHtjdX17euLvfkuQtSbJp06Z+zGPOWmX3rDePfWzddr/7rMUVAgDsNuY/XyQ+YwBbtnPPF7V6kyx+\nBMYN3X3xstvXlrVZPjyhM1d3VX1/ktOTnNHdJ9yhYff5Sc5NckRV/ViSTUnetob9b/WYy1XV/kme\nmOSE6RtIbk1yYabRHVvaDgAAAFibRY/A2CbTxJqnJrkuW75c461Jjkmyb5KPdfdFO6CUo5J8YoUa\nfinJi6vqhOmyEgAAAOBuWHSAsXdVbVi27DvdfdUat39tkickeVqSfapuG3Zy7TQHRpKcluQNSV6W\n5Ne3sd47qap7JDk8yfHdfcGydVdnNi/HIUk+ur2PDQAAAOvFoi8heVpm3/YxfzvvLmz/lCT3S3L+\nsn08b6lBd1+X2beI3DT93N6eNdVwp28o6e6vZDYB6BYn8wQAAABWV+vhyoaq+sskX+zuly66lpVs\n2rSpN2/evOgyGMzciKKsh/cpALDjzX++SHzGALZsZ54vquqc7t60WrtFX0KyQ1XVPkmelOTQJAcs\nuBwAAADgbtqtA4zMLke5b5Jjl89PAQAAAOw6dusAo7s3LroGAAAAYNstehJPAAAAgFUJMAAAAIDh\nCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTBgUPvtt98dfgIAbKv5\nzxU+YwBbM+L5Ys9FFwCs7Iorrlh0CQDAbsbnC2CtRjxfGIEBAAAADE+AAQAAAAxPgAEAAAAMT4AB\nAAAADE+AAQAAAAxPgAEAAAAMT4ABAAAADE+AAQAAAAxPgAEAAAAMT4ABAAAADE+AAQAAAAxPgAEA\nAAAMT4ABAAAADE+AAQAAAAxPgAEAAAAMT4ABAAAADE+AAQAAAAxPgAEAAAAMT4ABAAAADE+AAQAA\nAAxPgAEAAAAMT4ABAAAADE+AAQAAAAxPgAEAAAAMT4ABAAAADE+AAQAAAAxPgAEAAAAMT4ABAAAA\nDE+AAQAAAAxPgAEAAAAMT4ABAAAADE+AAQAAAAxPgAEAAAAMT4ABAAAADE+AwVA2bNiQqsqGDRsW\nXQrAbZybYOu8RwDuuqVzp/Pn2lV3L7qGdW/Tpk29efPmRZcxhKq67b6+CYzCuQm2znsE4K6bP3cm\n6/v8WVXndPem1doZgQEAAAAMT4ABAAAADE+AAQAAAAxPgAEAAAAMT4ABAAAADE+AAQAAAAxPgAEA\nAAAMT4ABAAAADE+AAQAAAAxPgAEAAAAMT4ABAAAADE+AAQAAAAxvYQFGVZ1cVb3C7eOLqunuqqqD\nlz2HG6vq01V1+KJrAwAAgN3Bngs+/keS/OqyZTcvopDt5EeTXJPk+5I8N8nJVfXZ7t7lQhkAAAAY\nyaIvIbmpu69YdrtmaeU0muHoqnpPVX2rqi6pqsPm1m9pFMcRVXV4VV1dVXvPH7CqTq2q06f7x1XV\nBVX1oqq6dDrGSVW1V1W9vKoun/bxhqpay2v11ek5XNLdv5tZmPGY7fRaAQAAwLq16BEYa/FbSV6V\n5NVJjkrytqo6u7u/kOQV07olL07ymiSbk3w+yZuTPCfJu5Okqu6d5BeSvGBum41Tm2cmeVCS9yV5\nQJKvJDk0yf7T9h+b1q1qCjueneQ+Sf5+1Q0uuig5+OC17Hq3d+b8A68JMAjnJtg67xGAu+7M5Quc\nP1e16BEYT6+q65fdTlzW5h3dfUp3X5zkN5PcmuTJSdLd1y6N3Ejy8MzCjsO6+4LuvjHJqZmFGkte\nmOSbST40t2yPJEdO25yR5MNJHpfk17r7wu5+f2bhxSFreD6XVtX1SW5K8t4kx3T35pUaTiNLNlfV\n5ltuuWUNuwYAAID1a9EjMM5OcvSyZd9Y9vjTS3e6+9aquirJ/ecbVNXGzEZHvG4KHJa8Ncm5VfXg\n7v5iZmHG27v71rk2X+jua+ceX5nks91987JldzjmFhyS5OtJ9k5yUJL/XFXXdfdblzfs7rckeUuS\nbNq0qXPWWWvY/e7vkKrb7rfXBBiEcxNsnfcIwF03f+5M1vn5c9lrsSWLDjBumEZWbM3y4QmduZEj\nVfX9SU5PckZ3n3CHht3nV9W5SY6oqg8k2ZTksNzRSvtfadkeq9SZJP+7u7823f/HqvqJzEaF3CnA\nAAAAANZu0QHGNpnmmjg1yXVJXrKFZm9NckySfZN8rLsv2knlJcl3MvtGEgAAAGAbLDrA2LuqNixb\n9p3uvmqN2782yROSPC3JPnX7sJNrpzkwkuS0JG9I8rIkv76N9a7m/lW1Z26/hORXM5sLAwAAANgG\niw4wnpbZt33M+1KSB69x+6ckuV+S85ctPzLJyUnS3ddV1buTPDfTt5HsQP84/bw1yeVJ/ijJcTv4\nmAAAALDbq+5edA07XFX9ZZIvdvdLF13LSjZt2tSbN6/4ZSXrTs1PArYO+iawa3Bugq3zHgG462r5\nJJ7r+PxZVed096bV2i16BMYOVVX7JHlSkkOTHLDgcgAAAIC7abcOMJKcl+S+SY7t7gsWXQwAAABw\n9+zWAUZ3b1x0DQAAAMC2+55FFwAAAACwGgEGAAAAMDwBBgAAADA8AQYAAAAwPAEGAAAAMDwBBgAA\nADA8AQYAAAAwPAEGAAAAMDwBBkPZb7/97vATYATOTbB13iMAd938OdP5c232XHQBMO+KK65YdAkA\nd+LcBFvnPQJw1zl33nVGYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMT\nYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNg\nAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AA\nAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAAAADDE2AAAAAAwxNgAAAAAMMTYAAA\nAADDE2AAAAAAwxNgAAAAAMMTYDCUDRs2pKqyYcOGRZeycF4LGIf3I2yd98iuY+l35fcFizf6+3HE\n+qq7F13Durdp06bevHnzossYQlXddn+9902vBYzD+xG2zntk1zH/u0r8vmCRRn8/7sz6quqc7t60\nWjsjMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOEJMAAAAIDhCTAAAACA\n4QkwAAAAgOEJMAAAAIDhCTAAAACA4QkwAAAAgOGtGmBU1clV1SvcPr4zCtxRqur8qrq1qh6xHfZ1\n6dzr8p2q+lJV/WFV/cD2qBUAAADWu7WOwPhIkgcsu/3sjipqR6uqg5LcP8mfJjlqDe1PrqrjVmn2\nusxel4ckeVFmr8/vblulAAAAQLL2AOOm7r5i2e2apZXTyIOjq+o9VfWtqrqkqg6bW7+lURxHVNXh\nVXV1Ve09f8CqOrWqTp/uH1dVF1TVi6bRDt+qqpOqaq+qenlVXT7t4w1VtZbndFSSdyY5KcnhVbXn\nGl+Hrbluel2+1N0fSfLuJI/dDvsFAACAdW97/MN9yW8leVWSV2cWELytqs7u7i8kecW0bsmLk7wm\nyeYkn0/y5iTPyewf/amqeyf5hSQvmNtm49TmmUkelOR9mY14+EqSQ5PsP23/sWndiqrqXkmen+Tg\nJJ9KcuO0zw/czee90jEekuRnkpy5lvYXXZQcfPD2Ovqu7vaXzGvitYBxeD/C1nmP7Dru+PHU7wsW\nafT343j1rXUExtOr6vpltxOXtXlHd5/S3Rcn+c0ktyZ5cpJ097VLIzeSPDyzsOOw7r6gu29Mcmpm\nocaSFyb5ZpIPzS3bI8mR0zZnJPlwkscl+bXuvrC7359ZeHHIKs/ll5Nc3t3ndXcnOSXJS9b4OmzN\n8dPrcmOSy5J8PcmxW2o8jVjZXFWbb7nllu1weAAAANh91ezf8FtpUHVyZvM6HL1s1Te6+2tTm07y\nwu4+bW67y5K8ubvfMLdsY5JPJnlTd58wt/yAJOcmeWh3f7Gq/j7JR7v730/rj0vygu5+5Nw2/yXJ\ngd39xLll/216Tv96K8/nb5L8RXefOD1+eJILkzyku788LTs2dwwf9k7SSW6eW/aM7v5fU/tLk5yW\n5E+SVJIfTnJCkpuSHNLd391SPUmyadOm3rx589aarBtVddv91frm7s5rAePwfoSt8x7Zdcz/rhK/\nL1ik0d+PO7O+qjqnuzet1m6tl5DcMI2s2Jrlwwg6cyM8qur7k5ye5Iz58CJJuvv8qjo3yRFV9YEk\nm5Icljtaaf8rLdtjSwVW1f5JnpjkCVV1/NyqPZIcmWRp2R9mupxlcmKSLyX5T3PLvrRs91fPvUaf\nq6pXJvm7zEaE/PWWagIAAABWtz3nwNiiaWLNU5Ncly1frvHWJMck2TfJx7r7oh1QylFJPrFCDb+U\n5MVVdULPXJNkfpLS65Jcs4YQZ953pp/fty0FAwAAAGsPMPauqg3Lln2nu69a4/avTfKEJE9Lss/c\nUJRrpzkwktklGG9I8rIkv77G/a5ZVd0jyeFJju/uC5atuzqzeTkOSfLRu3mIH5heo6VLSH43yVVJ\n/vZuFw0AAAAkWfsknk/L7Ns+5m/n3YXjPCXJ/ZKcv2wfz1tq0N3XZXbZxk254+Ub28uzphru9A0l\n3f2VzCYA3ZbJPH8rs+f05SQfTPKtJId299XbsE8AAAAga5jEc2eqqr9M8sXufumia9mZTOJ5O5OA\n3c5rAePwfoSt8x7ZdYw+aSCsJ6O/H3flSTx3qKraJ8mTkhya5IAFlwMAAAAMZogAI7PLUe6b5Njl\n81MAAAAADBFgdPfGRdcAAAAAjGutk3gCAAAALIwAAwAAABieAAMAAAAYngADAAAAGJ4AAwAAABie\nAAMAAAAYngADAAAAGJ4AAwAAABieAIOh7Lfffnf4uZ55LWAc3o+wdd4ju47535HfFyzW6O/HEeur\n7l50Devepk2bevPmzYsuAwAAAHa6qjqnuzet1s4IDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4\nAgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgC\nDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIM\nAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwA\nAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAA\nAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAA\nYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB4AgwAAABgeAIMAAAAYHgCDAAAAGB41d2LrmHd\nq6rrkly06DoYzr5JvrboIhiOfsFK9AtWol+wEv2ClegXrGRn9ouHdvf9Vmu0586ohFVd1N2bFl0E\nY6mqzfoFy+kXrES/YCX6BSvRL1iJfsFKRuwXLiEBAAAAhifAAAAAAIYnwBjDWxZdAEPSL1iJfsFK\n9AtWol+wEv2ClegXrGS4fmESTwAAAGB4RmAAAAAAwxNgAAAAAMMTYCxQVT29qi6qqour6lWLroed\np6reVlVfraoL5pbdt6r+qqo+N/3cZ1peVfWfpn7y6ap67OIqZ0eqqh+uqjOr6jNV9Y9V9Yppub6x\njlXVPavqk1V1/tQvfnta/iNV9Ynp9//fqmqvafne0+OLp/UbF1k/O1ZV7VFV51XVB6fH+sU6V1WX\nVtU/VNWnqmrztMzfkXWuqu5TVe+tqn+qqgur6gn6BVX1yOlcsXT7ZlW9cuS+IcBYkKraI8l/SfKM\nJI9O8oKqevRiq2InOjnJ05cte1WSv+7uhyf56+lxMusjD59uRyf5g51UIzvfrUn+n+5+dJLHJ/mN\n6bygb6xvNyX5qe4+IMmBSZ5eVY9PcmKSN3b3P0/y9SRHTe2PSvL1afkbp3bsvl6R5MK5x/oFSXJI\ndx/Y3Zumx/6O8OYkH+7u/ZMckNl5Q79Y57r7oulccWCSxyW5Icn7M3DfEGAszkFJLu7uS7r75iTv\nSvKcBdfETtLdZye5Ztni5yR5+3T/7Ul+fm75n/bMx5Pcp6oesHMqZWfq7q9097nT/esy+3DxoOgb\n69r0+71+eniP6dZJfirJe6fly/vFUn95b5KnVlXtpHLZiarqwUl+LskfT48r+gUr83dkHauqeyd5\ncpI/SZLuvrm7vxH9gjt6apLPd/dlGbhvCDAW50FJLp97/MVpGevXft39len+FUn2m+7rK+vQNLz7\nMUk+EX1j3ZsuE/hUkq8m+askn0/yje6+dWoy/7u/rV9M669N8kM7t2J2kjclOSbJd6fHPxT9glnA\n+T+q6pyqOnpa5u/I+vYjSa5KctJ0ydkfV9W9ol9wR89Pctp0f9i+IcCAAfXs+419x/E6VVXfn+R9\nSV7Z3d+cX6dvrE/d/Z1peOeDMxvBt/+CS2LBquqZSb7a3ecsuhaG85Pd/djMhnr/RlU9eX6lvyPr\n0p5JHpvkD7r7MUm+ldsvCUiiX6x303xJz07ynuXrRusbAozF+VKSH557/OBpGevXlUtDsKafX52W\n6yvrSFXdI7Pw4tTu/rNpsb5BkmQa8ntmkidkNmxzz2nV/O/+tn4xrb93kqt3cqnseE9M8uyqujSz\ny1B/KrNr3PWLda67vzT9/Gpm17IfFH9H1rsvJvlid39ievzezAIN/YIlz0hybndfOT0etm8IMBbn\n75M8fJotfK/MhuycvuCaWKzTk7xouv+iJH8+t/zwadbfxye5dm5IF7uR6Xr0P0lyYXe/YW6VvrGO\nVdX9quo+0/3vTfLTmc2PcmaS507NlveLpf7y3CQfnf73hN1Id7+6ux/c3Rsz+wzx0e7+legX61pV\n3auqfmDpfpJDk1wQf0fWte6+IsnlVfXIadFTk3wm+gW3e0Fuv3wkGbhvlL9di1NVP5vZ9at7JHlb\ndx+/4JLYSarqtCQHJ9k3yZVJXpvkA0neneQhSS5L8q+7+5rpH7W/n9m3ltyQ5Mju3ryIutmxquon\nk/yvJP+Q269pPzazeTD0jXWqqn48swm09sjsPx7e3d2vq6qHZfY/7/dNcl6Sw7r7pqq6Z5J3ZDaH\nyjVJnt/dlyymenaGqjo4yb/r7mfqF+vb9Pt///RwzyTv7O7jq+qH4u/IulZVB2Y24e9eSS5JcmSm\nvynRL9a1Kez8QpKHdfe107JhzxkCDAAAAGB4LiEBAAAAhifAAAAAAIYnwAAAAACGJ8AAAAAAhifA\nAAAAAIYnwAAAdklVdXBVdVXtu0q7s6rq93dSTU+pqs9W1R7T4+Oq6oKttD+iqq5fZZ+/UVV/sb1r\nBYBdjQADANghqurkKWDoqrqlqi6pqtdP3zm/PfxtkgckuXo63pbCgF9M8urtdMzV/F6S47v7O9tx\nn3+c5HFV9aTtuE8A2OUIMACAHekjmYUMD0vyH5K8PMnrt8eOu/vm7r6iu3uVdtd093Xb45hbU1X/\nKsn+Sd69Pffb3TcleWeSf7s99wsAuxoBBgCwI900hQyXd/c7k5ya5OeXVlbVk6vqE1X17aq6sqre\nWFV7LVv/8aq6vqqurapPVtWPTetuu4Skqg5OclKSe82N+jhuaneHS0iqap+qentVfb2qbqyqj1TV\nj86tP2I63lOr6oKq+lZVnVlVP7LKc31hkr/u7hu31KCqHlJV/zQdf8+55c+aLj359nSshy3b9PQk\nz66q71ulBgDYbQkwAICd6cYk90iSqnpQkr9Mcl6SxyQ5KskLkvz/0/o9k/x5kr9JckCSn0jypiQr\nXZ7xt0l4fIzDAAADoElEQVRemeSGzEZ8PCBbHulx8rSv5yQ5aNrmw1X1vXNt9s7sspMXJ3lCkvsk\n+cNVntuTkmze0sqqelSSjyX570mO6O5b54712iRHTsfaI8mfVVXNbb45yZ7TegBYl/ZcvQkAwLar\nqoMyjVKYFr08yZeTvLy7v5vkwqp6VZI/qqrfTHLPzIKDv+juz0/b/NNK++7um6vq2tndvmIrNTw8\nybOTPKW7z56W/WqSLyT5lczmm0hmn5F+o7svmtq8Psnbqqq2csnKQ6fns9JxfyLJh5K8sbuPX7Z6\nzySv6O6PzdVzSZKnZnYJTrr7hun5bdzScwOA3Z0RGADAjvT06XKMbyf5uyRnJ/k307pHJfn4FF4s\n+ZskeyX55919TWajJc6oqg9V1f9dVQ/ZxnoeleS7Uy1Jku6+Nsk/JHn0XLublsKLyZenuvbZyr6/\nN8m3V1j+oMyCiBNXCC8y1fPJuXoum4736GXtbpyOAQDrkgADANiRzk5yYJJHJrlnd/9id391Ddt1\nknT3kZld7nF2ZiMnLqqqn9lBtc6PrLh1C+u29tnpa1k54Phako8neX5VbSkA2epEpJP7JrlqDe0A\nYLckwAAAdqQbuvvi7r6su29Ztu7CJI+vqvnPIz+Z5OYkS5eMpLvP7+4Tu/vgJGcledEWjnVzZvNH\nbM2FmX3+uW0uiar6wST/IslnVn86W3Ve7jxqIkluyix8+XqSv6qq+yxb/z2ZzcWxVM9DkjxwqnVp\n2T/L7JKac7exRgDYZQkwAIBF+a+Z/UP9v1bVo6rq55L8TpLfn+Z8+JGq+p2q+ldV9dCqOiTJj2fL\nQcOlSe5ZVT89fTPJnb6xo7s/l9nEoH9UVU+qqn+R5JQk38zsq0q3xRmZBTB3Mn0zybOSXJs7hxi3\nJnlTVT2hqg5M8vYk/5hp/ovJk5JcMtUPAOuSAAMAWIju/lKSZ2T2DSSfSvK2JKclOXZqckOSRyR5\nT5LPZvYP+1OTnLiF/f1tZt8Uclpml1ocs4VDH5nZnBOnTz+/L8nTt/b1p2t0SpJHzH8l67L6bkzy\nzMzCkvkQ46Ykxyf50ySfyOzz2S8umyz0BUneuo31AcAurbY8kTYAAHdFVf1Okvt191HbcZ8/ltk3\ntzximnAUANYlIzAAALafE5JcUlWrzcVxVzwwyeHCCwDWOyMwAAAAgOEZgQEAAAAMT4ABAAAADE+A\nAQAAAAxPgAEAAAAMT4ABAAAADE+AAQAAAAzv/wAUIlv1xqnZjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10df30850>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fragments_a = sp.array([20,100,200,400])\n",
    "fragments_b = sp.array([50,150,220,250])\n",
    "fragments_ab = sp.array([20,30,50,50,80,150,170,170])\n",
    "\n",
    "#Berechne Lösungen\n",
    "results = double_digest(fragments_a,fragments_b,fragments_ab)\n",
    "#Plotte Lösungen\n",
    "plot_solutions(results)\n",
    "\n",
    "for i,result in enumerate(results):\n",
    "    print(\"Positionen Enzym A:\\t%s\"%(result['posA']))\n",
    "    print(\"Positionen Enzym B:\\t%s\"%(result['posB']))\n",
    "    print(\"Loesung %d:\\t\\t%s\"%(i+1,result['positions']))\n",
    "    print"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}