6-exercise-solutions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tskit, msprime\n",
    "import io\n",
    "from IPython.display import SVG\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notebook 1\n",
    "\n",
    "*Exercise*: Can you modify the following Table Collection until the corresponding tree sequence looks like the one in this plot?\n",
    "\n",
    "<img src=\"pics/worksheet1-pic-for-example.png\" width=\"400\" height=\"400\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg baseProfile=\"full\" height=\"200\" version=\"1.1\" width=\"400\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs><style type=\"text/css\"><![CDATA[.background path {fill: #808080; fill-opacity: 0}.background path:nth-child(odd) {fill-opacity: .1}.axes {font-size: 14px}.x-axis .tick .lab {font-weight: bold; dominant-baseline: hanging}.axes, .tree {font-size: 14px; text-anchor: middle}.axes line, .edge {stroke: black; fill: none}.y-axis .grid {stroke: #FAFAFA}.node > .sym {fill: black; stroke: none}.site > .sym {stroke: black}.mut text {fill: red; font-style: italic}.mut.extra text {fill: hotpink}.mut line {fill: none; stroke: none}.mut .sym {fill: none; stroke: red}.mut.extra .sym {stroke: hotpink}.node .mut .sym {stroke-width: 1.5px}.tree text, .tree-sequence text {dominant-baseline: central}.plotbox .lab.lft {text-anchor: end}.plotbox .lab.rgt {text-anchor: start}]]></style></defs><g class=\"tree-sequence\"><g class=\"background\"><path d=\"M20,0 l180,0 l0,138.2 l0,25 l0,5 l-180,0 l0,-5 l0,-25 l0,-138.2z\"/><path d=\"M200,0 l180,0 l0,138.2 l0,25 l0,5 l-180,0 l0,-5 l0,-25 l0,-138.2z\"/></g><g class=\"axes\"><g class=\"x-axis\"><g transform=\"translate(200 200)\"><text class=\"lab\" text-anchor=\"middle\" transform=\"translate(0 -11)\">Genome position</text></g><line x1=\"20\" x2=\"380\" y1=\"163.2\" y2=\"163.2\"/><g class=\"ticks\"><g class=\"tick\" transform=\"translate(20 163.2)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"5\"/><g transform=\"translate(0 6)\"><text class=\"lab\">0.00</text></g></g><g class=\"tick\" transform=\"translate(200 163.2)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"5\"/><g transform=\"translate(0 6)\"><text class=\"lab\">0.50</text></g></g><g class=\"tick\" transform=\"translate(380 163.2)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"5\"/><g transform=\"translate(0 6)\"><text class=\"lab\">1.00</text></g></g></g></g></g><g class=\"plotbox trees\"><g class=\"tree t0\" transform=\"translate(20 0)\"><g class=\"plotbox\"><g class=\"c2 node n9 p0 root\" transform=\"translate(100.5 38.625)\"><g class=\"a9 c2 node n6 p0\" transform=\"translate(31.5 47.3)\"><g class=\"a6 leaf node n2 p0 sample\" transform=\"translate(-14 35.475)\"><path class=\"edge\" d=\"M 0 0 V -35.475 H 14\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">2</text></g><g class=\"a6 leaf node n3 p0 sample\" transform=\"translate(14 35.475)\"><path class=\"edge\" d=\"M 0 0 V -35.475 H -14\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">3</text></g><path class=\"edge\" d=\"M 0 0 V -47.3 H -31.5\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">6</text></g><g class=\"a9 c2 node n7 p0\" transform=\"translate(-31.5 35.475)\"><g class=\"a7 leaf node n4 p0 sample\" transform=\"translate(21 47.3)\"><path class=\"edge\" d=\"M 0 0 V -47.3 H -21\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">4</text></g><g class=\"a7 c2 node n5 p0\" transform=\"translate(-21 23.65)\"><g class=\"a5 leaf node n0 p0 sample\" transform=\"translate(-14 23.65)\"><path class=\"edge\" d=\"M 0 0 V -23.65 H 14\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">0</text></g><g class=\"a5 leaf node n1 p0 sample\" transform=\"translate(14 23.65)\"><path class=\"edge\" d=\"M 0 0 V -23.65 H -14\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">1</text></g><path class=\"edge\" d=\"M 0 0 V -23.65 H 21\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">5</text></g><path class=\"edge\" d=\"M 0 0 V -35.475 H 31.5\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">7</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab\" transform=\"translate(0 -11)\">9</text></g></g></g><g class=\"tree t1\" transform=\"translate(200 0)\"><g class=\"plotbox\"><g class=\"c2 node n10 p0 root\" transform=\"translate(111 26.8)\"><g class=\"a10 leaf node n4 p0 sample\" transform=\"translate(35 94.6)\"><path class=\"edge\" d=\"M 0 0 V -94.6 H -35\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">4</text></g><g class=\"a10 c2 node n8 p0\" transform=\"translate(-35 23.65)\"><g class=\"a8 c2 node n5 p0\" transform=\"translate(-28 47.3)\"><g class=\"a5 leaf node n0 p0 sample\" transform=\"translate(-14 23.65)\"><path class=\"edge\" d=\"M 0 0 V -23.65 H 14\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">0</text></g><g class=\"a5 leaf node n1 p0 sample\" transform=\"translate(14 23.65)\"><path class=\"edge\" d=\"M 0 0 V -23.65 H -14\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">1</text></g><path class=\"edge\" d=\"M 0 0 V -47.3 H 28\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">5</text></g><g class=\"a8 c2 node n6 p0\" transform=\"translate(28 35.475)\"><g class=\"a6 leaf node n2 p0 sample\" transform=\"translate(-14 35.475)\"><path class=\"edge\" d=\"M 0 0 V -35.475 H 14\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">2</text></g><g class=\"a6 leaf node n3 p0 sample\" transform=\"translate(14 35.475)\"><path class=\"edge\" d=\"M 0 0 V -35.475 H -14\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">3</text></g><path class=\"edge\" d=\"M 0 0 V -35.475 H -28\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">6</text></g><path class=\"edge\" d=\"M 0 0 V -23.65 H 35\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">8</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab\" transform=\"translate(0 -11)\">10</text></g></g></g></g></g></svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Modify this code.\n",
    "nodes_ex = io.StringIO(\"\"\"\\\n",
    "id      is_sample   population      time\n",
    "0       1       0               0.0\n",
    "1       1       0               0.0\n",
    "2       1       0               0.0\n",
    "3       1       0               0.0\n",
    "4       1       0               0.0\n",
    "5       0       0               1.0\n",
    "6       0       0               1.5\n",
    "7       0       0               2.0\n",
    "8       0       0               3.0\n",
    "9       0       0               3.5\n",
    "10      0       0               4.0\n",
    "\"\"\")\n",
    "edges_ex = io.StringIO(\"\"\"\\\n",
    "id      left            right           parent  child\n",
    "0       0.00000000      1.00000000      5       0\n",
    "1       0.00000000      1.00000000      5       1\n",
    "2       0.00000000      1.00000000      6       2\n",
    "3       0.00000000      1.00000000      6       3\n",
    "4       0.00000000      0.50000000      7       4\n",
    "5       0.00000000      0.50000000      7       5\n",
    "6       0.50000000      1.00000000      8       5\n",
    "7       0.50000000      1.00000000      8       6\n",
    "8       0.00000000      0.50000000      9       6\n",
    "9       0.00000000      0.50000000      9       7\n",
    "10      0.50000000      1.00000000      10      4\n",
    "11      0.50000000      1.00000000      10      8\n",
    "\"\"\")\n",
    "\n",
    "# GT: I added edges 4 and 10.\n",
    "\n",
    "# Load the tree sequence.\n",
    "ts_ex = tskit.load_text(nodes=nodes_ex, edges=edges_ex, strict=False)\n",
    "\n",
    "# Test by plotting it.\n",
    "SVG(ts_ex.draw_svg())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notebook 2\n",
    "<img src=\"pics/worksheet2-ex.png\" width=\"400\" height=\"400\">\n",
    "\n",
    "*Exercise:* The code below specifies a simulation of two populations in which population 1 and population 0 are created by a population split at time 200.\n",
    "\n",
    "Can you add a third contemporary population, population 2, so that populations 1 and 2 are created by a more recent population split  at time 150? (You may need to create some other ancestral populations too.)\n",
    "\n",
    "*GT: There was a small typo in my 'helper' code that might have made this a bit more frustrating than I meant it to be! Sorry!!*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"margin-left:20px\"><div><style scoped=\"\">\n",
       "            .tskit-table thead tr th:only-of-type {vertical-align: middle;}\n",
       "            .tskit-table thead tr th {text-align: center;vertical-align: top;}\n",
       "            .tskit-table tbody td {text-align: right;padding: 0.5em 0.5em;}\n",
       "            .tskit-table tbody th {padding: 0.5em 0.5em;}\n",
       "        </style><b>Populations (5)</b><table border=\"1\" class=\"tskit-table\"><thead><tr><th>id</th><th>name</th><th>description</th><th>initial_size</th><th>growth_rate</th><th>default_sampling_time</th><th>extra_metadata</th></tr></thead><tbody><tr><td>0</td><td>R</td><td>Plotted in red.</td><td>500.0</td><td>0</td><td>0</td><td>{}</td></tr><tr><td>1</td><td>B</td><td>Plotted in blue.</td><td>500.0</td><td>0</td><td>0</td><td>{}</td></tr><tr><td>2</td><td>G</td><td>Plotted in green.</td><td>500.0</td><td>0</td><td>0</td><td>{}</td></tr><tr><td>3</td><td>AncestralPopulationGB</td><td>Plotted in purple.</td><td>500.0</td><td>0</td><td>1.5e+02</td><td>{}</td></tr><tr><td>4</td><td>AncestralPopulationRGB</td><td>Plotted in purple.</td><td>500.0</td><td>0</td><td>2e+02</td><td>{}</td></tr></tbody></table></div><div><style scoped=\"\">\n",
       "            .tskit-table thead tr th:only-of-type {vertical-align: middle;}\n",
       "            .tskit-table thead tr th {text-align: center;vertical-align: top;}\n",
       "            .tskit-table tbody td {text-align: right;padding: 0.5em 0.5em;}\n",
       "            .tskit-table tbody th {padding: 0.5em 0.5em;}\n",
       "        </style><b>Migration matrix (all zero)</b><table border=\"1\" class=\"tskit-table\"><thead><tr></tr></thead><tbody></tbody></table></div><div><style scoped=\"\">\n",
       "            .tskit-table thead tr th:only-of-type {vertical-align: middle;}\n",
       "            .tskit-table thead tr th {text-align: center;vertical-align: top;}\n",
       "            .tskit-table tbody td {text-align: right;padding: 0.5em 0.5em;}\n",
       "            .tskit-table tbody th {padding: 0.5em 0.5em;}\n",
       "        </style><b>Events (2)</b><table border=\"1\" class=\"tskit-table\"><thead><tr><th>time</th><th>type</th><th>parameters</th><th>effect</th></tr></thead><tbody><tr><td>150</td><td><a href='https://tskit.dev/msprime/docs/latest/api.html#msprime.Demography.add_population_split'>Population Split</a></td><td>derived=[G, B], ancestral=AncestralPopulationGB</td><td>Moves all lineages from derived populations &#x27;G&#x27; and &#x27;B&#x27; to the ancestral &#x27;AncestralPopulationGB&#x27; population. Also set the derived populations to inactive, and all migration rates to and from the derived populations to zero.</td></tr><tr><td>200</td><td><a href='https://tskit.dev/msprime/docs/latest/api.html#msprime.Demography.add_population_split'>Population Split</a></td><td>derived=[AncestralPopulationGB, R], ancestral=AncestralPopulationRGB</td><td>Moves all lineages from derived populations &#x27;AncestralPopulationGB&#x27; and &#x27;R&#x27; to the ancestral &#x27;AncestralPopulationRGB&#x27; population. Also set the derived populations to inactive, and all migration rates to and from the derived populations to zero.</td></tr></tbody></table></div></div>"
      ],
      "text/plain": [
       "Demography(populations=[Population(initial_size=500, growth_rate=0, name='R', description='Plotted in red.', extra_metadata={}, default_sampling_time=None, initially_active=None, id=0), Population(initial_size=500, growth_rate=0, name='B', description='Plotted in blue.', extra_metadata={}, default_sampling_time=None, initially_active=None, id=1), Population(initial_size=500, growth_rate=0, name='G', description='Plotted in green.', extra_metadata={}, default_sampling_time=None, initially_active=None, id=2), Population(initial_size=500, growth_rate=0, name='AncestralPopulationGB', description='Plotted in purple.', extra_metadata={}, default_sampling_time=150, initially_active=False, id=3), Population(initial_size=500, growth_rate=0, name='AncestralPopulationRGB', description='Plotted in purple.', extra_metadata={}, default_sampling_time=200, initially_active=False, id=4)], events=[PopulationSplit(time=150, derived=['G', 'B'], ancestral='AncestralPopulationGB'), PopulationSplit(time=200, derived=['AncestralPopulationGB', 'R'], ancestral='AncestralPopulationRGB')], migration_matrix=array([[0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0.],\n",
       "       [0., 0., 0., 0., 0.]]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Setup.\n",
    "dem = msprime.Demography()\n",
    "dem.add_population(name=\"R\", description=\"Plotted in red.\", initial_size=500)\n",
    "dem.add_population(name=\"B\", description=\"Plotted in blue.\",initial_size=500)\n",
    "dem.add_population(name=\"G\", description=\"Plotted in green.\",initial_size=500)\n",
    "dem.add_population(name=\"AncestralPopulationGB\",\n",
    "                   description=\"Plotted in purple.\", initial_size=500)\n",
    "dem.add_population(name=\"AncestralPopulationRGB\",\n",
    "                   description=\"Plotted in purple.\", initial_size=500)\n",
    "\n",
    "# Add the population split.\n",
    "dem.add_population_split(time=150, derived=[\"G\",\"B\"], ancestral=\"AncestralPopulationGB\")\n",
    "dem.add_population_split(time=200, derived=[\"AncestralPopulationGB\",\"R\"],\n",
    "                         ancestral=\"AncestralPopulationRGB\")\n",
    "dem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate!\n",
    "ts_ex = msprime.sim_ancestry(\n",
    "  samples={\"R\" : 1, \"B\" : 1, \"G\" : 1},\n",
    "  demography=dem,\n",
    "  sequence_length=1000,\n",
    "  random_seed=63\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tree on interval: Interval(left=0.0, right=1000.0)\n"
     ]
    },
    {
     "data": {
      "image/svg+xml": [
       "<svg baseProfile=\"full\" height=\"400\" version=\"1.1\" width=\"500\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs><style type=\"text/css\"><![CDATA[.background path {fill: #808080; fill-opacity: 0}.background path:nth-child(odd) {fill-opacity: .1}.axes {font-size: 14px}.x-axis .tick .lab {font-weight: bold; dominant-baseline: hanging}.axes, .tree {font-size: 14px; text-anchor: middle}.axes line, .edge {stroke: black; fill: none}.y-axis .grid {stroke: #FAFAFA}.node > .sym {fill: black; stroke: none}.site > .sym {stroke: black}.mut text {fill: red; font-style: italic}.mut.extra text {fill: hotpink}.mut line {fill: none; stroke: none}.mut .sym {fill: none; stroke: red}.mut.extra .sym {stroke: hotpink}.node .mut .sym {stroke-width: 1.5px}.tree text, .tree-sequence text {dominant-baseline: central}.plotbox .lab.lft {text-anchor: end}.plotbox .lab.rgt {text-anchor: start}]]></style></defs><g class=\"tree t0\"><g class=\"plotbox\"><g class=\"c2 node n10 p4 root\" transform=\"translate(273.958 26.8)\"><g class=\"a10 c2 node n8 p4\" transform=\"translate(129.375 286.395)\"><g class=\"a8 i0 leaf node n1 p0 sample\" transform=\"translate(-38.3333 55.0046)\"><path class=\"edge\" d=\"M 0 0 V -55.0046 H 38.3333\"/><rect class=\"sym\" height=\"6\" style=\"fill:red;\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">1</text></g><g class=\"a8 i1 leaf node n2 p1 sample\" transform=\"translate(38.3333 55.0046)\"><path class=\"edge\" d=\"M 0 0 V -55.0046 H -38.3333\"/><rect class=\"sym\" height=\"6\" style=\"fill:blue;\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">2</text></g><path class=\"edge\" d=\"M 0 0 V -286.395 H -129.375\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\" style=\"fill:orange;\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">8</text></g><g class=\"a10 c2 node n9 p4\" transform=\"translate(-129.375 173.469)\"><g class=\"a9 i0 leaf node n0 p0 sample\" transform=\"translate(-86.25 167.931)\"><path class=\"edge\" d=\"M 0 0 V -167.931 H 86.25\"/><rect class=\"sym\" height=\"6\" style=\"fill:red;\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">0</text></g><g class=\"a9 c2 node n7 p4\" transform=\"translate(86.25 116.448)\"><g class=\"a7 i2 leaf node n5 p2 sample\" transform=\"translate(57.5 51.4821)\"><path class=\"edge\" d=\"M 0 0 V -51.4821 H -57.5\"/><rect class=\"sym\" height=\"6\" style=\"fill:green;\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">5</text></g><g class=\"a7 c2 node n6 p3\" transform=\"translate(-57.5 10.933)\"><g class=\"a6 i1 leaf node n3 p1 sample\" transform=\"translate(-38.3333 40.5491)\"><path class=\"edge\" d=\"M 0 0 V -40.5491 H 38.3333\"/><rect class=\"sym\" height=\"6\" style=\"fill:blue;\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">3</text></g><g class=\"a6 i2 leaf node n4 p2 sample\" transform=\"translate(38.3333 40.5491)\"><path class=\"edge\" d=\"M 0 0 V -40.5491 H -38.3333\"/><rect class=\"sym\" height=\"6\" style=\"fill:green;\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">4</text></g><path class=\"edge\" d=\"M 0 0 V -10.933 H 57.5\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\" style=\"fill:purple;\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">6</text></g><path class=\"edge\" d=\"M 0 0 V -116.448 H -86.25\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\" style=\"fill:orange;\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">7</text></g><path class=\"edge\" d=\"M 0 0 V -173.469 H 129.375\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\" style=\"fill:orange;\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">9</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\" style=\"fill:orange;\"/><text class=\"lab\" transform=\"translate(0 -11)\">10</text></g></g></g></svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot.\n",
    "colour_map = {0:\"red\", 1:\"blue\", 2:\"green\", 3: \"purple\", 4: \"orange\"}\n",
    "node_colours = {u.id: colour_map[u.population] for u in ts_ex.nodes()}\n",
    "for tree in ts_ex.trees():\n",
    "    print(\"Tree on interval:\", tree.interval)\n",
    "    # The code below will only work in a Jupyter notebook with SVG output enabled.\n",
    "    display(SVG(tree.draw(node_colours=node_colours, width=500, height=400)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notebook 3\n",
    "\n",
    "*Exercise*. For most of their evolutionary history, mutations in your study organisms were well described by the Jukes-Cantor model. \n",
    "However, they were all were exposed to an environmental mutagen 50 generations ago, and from then on all new mutations were of type 👽. All of the nucleotide bases have an equal probability of mutating to an allele of type 👽. Once a nucleotide mutates to a 👽, it cannot mutate back to a normal nucleotide. Modify the code below to simulate this scenario."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The underlying genealogy\n",
    "ts_ex = msprime.sim_ancestry(\n",
    "    10, sequence_length=1000, random_seed=1987, recombination_rate=0.001,\n",
    "    population_size=1e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Modify code here\n",
    "model = msprime.MatrixMutationModel(\n",
    "    alleles = [\"A\", \"C\", \"G\", \"T\", \"👽\"],\n",
    "    root_distribution = [0.25, 0.25, 0.25, 0.25, 0.0],\n",
    "    transition_matrix = [[0.0, 0.0, 0.0, 0.0, 1.0],\n",
    "                         [0.0, 0.0, 0.0, 0.0, 1.0],\n",
    "                         [0.0, 0.0, 0.0, 0.0, 1.0],\n",
    "                         [0.0, 0.0, 0.0, 0.0, 1.0],\n",
    "                         [0.0, 0.0, 0.0, 0.0, 1.0]]\n",
    ")\n",
    "\n",
    "mts_ex = msprime.sim_mutations(ts_ex, rate=1e-6, random_seed=752, start_time=50)\n",
    "\n",
    "mts_ex = msprime.sim_mutations(mts_ex, rate=1e-6, random_seed=851, start_time=0,\n",
    "                            end_time=50, model=model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "╔══╤════╤════╤═════════════╤═════════════╤══════╤════════╗\n",
      "║id│site│node│time         │derived_state│parent│metadata║\n",
      "╠══╪════╪════╪═════════════╪═════════════╪══════╪════════╣\n",
      "║0 │   0│  12│  26.19673780│            👽│    -1│     b''║\n",
      "║1 │   1│ 993│1080.73366012│            C│    -1│     b''║\n",
      "║2 │   2│   7│  45.06634352│            👽│    -1│     b''║\n",
      "║3 │   3│1840│ 508.53139009│            T│    -1│     b''║\n",
      "║4 │   4│   8│ 224.66157800│            G│    -1│     b''║\n",
      "║5 │   5│ 687│6529.29267485│            T│    -1│     b''║\n",
      "║6 │   6│3622│7099.53462718│            A│    -1│     b''║\n",
      "║7 │   7│  13│ 413.05691858│            G│    -1│     b''║\n",
      "║8 │   8│  12│ 175.08036637│            C│    -1│     b''║\n",
      "║9 │   9│3205│1671.37834344│            T│    -1│     b''║\n",
      "║10│  10│2619│2897.51611664│            T│    -1│     b''║\n",
      "║11│  11│  12│  48.97276474│            👽│    -1│     b''║\n",
      "║12│  12│   9│ 716.86318876│            A│    -1│     b''║\n",
      "║13│  13│1511│1651.37247430│            T│    -1│     b''║\n",
      "║14│  14│1006│5547.00775284│            G│    -1│     b''║\n",
      "╚══╧════╧════╧═════════════╧═════════════╧══════╧════════╝\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Check here\n",
    "print(mts_ex.tables.mutations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "35.0\t('C', '👽')\t[0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0]\n",
      "95.0\t('A', 'C')\t[0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1]\n",
      "101.0\t('C', '👽')\t[0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      "118.0\t('A', 'T')\t[1 0 0 1 0 1 0 1 1 0 0 0 0 0 1 0 1 0 1 1]\n",
      "142.0\t('T', 'G')\t[0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0]\n",
      "205.0\t('G', 'T')\t[0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0]\n",
      "268.0\t('T', 'A')\t[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1]\n",
      "333.0\t('A', 'G')\t[0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0]\n",
      "440.0\t('T', 'C')\t[0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0]\n",
      "456.0\t('A', 'T')\t[0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 1 0 0]\n",
      "662.0\t('A', 'T')\t[1 1 1 1 0 1 1 1 0 0 0 0 0 0 1 0 0 1 0 1]\n",
      "667.0\t('T', '👽')\t[0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0]\n",
      "694.0\t('T', 'A')\t[0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]\n",
      "815.0\t('A', 'T')\t[0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0]\n",
      "905.0\t('A', 'G')\t[0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0]\n"
     ]
    }
   ],
   "source": [
    "# Check here\n",
    "for var in mts_ex.variants():\n",
    "    print(var.site.position, var.alleles, var.genotypes, sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notebook 4\n",
    "\n",
    "*Exercise*. Modify the following code to compute $f_2$ statistics in 5Mb windows along the genome between the samples from populations SMALL and BIG, as well as between populations SMALL and ADMIX."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "demography = msprime.Demography()\n",
    "demography.add_population(name=\"SMALL\", initial_size=2000)\n",
    "demography.add_population(name=\"BIG\", initial_size=5000)\n",
    "demography.add_population(name=\"ADMIX\", initial_size=2000)\n",
    "demography.add_population(name=\"ANC\", initial_size=5000)\n",
    "demography.add_admixture(\n",
    "    time=100, derived=\"ADMIX\", ancestral=[\"SMALL\", \"BIG\"],\n",
    "    proportions=[0.5, 0.5])\n",
    "demography.add_population_split(\n",
    "    time=1000, derived=[\"SMALL\", \"BIG\"], ancestral=\"ANC\")\n",
    "\n",
    "ts = msprime.sim_ancestry(\n",
    "    samples={\"SMALL\": 100, \"BIG\": 100, \"ADMIX\" : 100},\n",
    "    demography=demography,\n",
    "    random_seed=2432,\n",
    "    sequence_length=5e7,\n",
    "    recombination_rate=1e-8\n",
    ")\n",
    "ts = msprime.sim_mutations(ts, rate=1e-8, random_seed=6151)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Windowed F2 values')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEWCAYAAAB2X2wCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABGDklEQVR4nO3dd3hU1dbA4d9KIaGEIoReQu8QioAUaSIdBKRdRbFhwd4ulnv16qeioNhFBMWKiAoCSlNaUAFD7x3pEEAglNTZ3x/nBIaQkDaTM5ms93nmycypaybJmn323mdvMcaglFIq7wtwOgCllFKeoQldKaX8hCZ0pZTyE5rQlVLKT2hCV0opP6EJXSml/IQmdOVRItJORLblYH8jIjU8GVMWz79XRG5w6vwZEZHFInK303Eo36QJXV2ViDwjInNSLduRzrIhxpgoY0zt3I0yd4jIZBFJEJGzbo/BIhIiIpNE5G8RiRWRtSLS3el4Vf6jCV1lZCnQWkQCAUSkHBAMNEm1rIa9rb97wxhTxO0xFQgC9gPtgWLA88B3IhLhYJwqH9KErjLyF1YCj7RftwMWAdtSLdtljDkkIh1E5EDKznYVxpMisl5ETovIVBEJdVv/lIgcFpFDInKn+4lFpJiIfCEiMXbp93kRCbDX/S0izeznt9hVNfXt13eJyAz7eYCIjBKRXSJyQkS+E5Fr3M4xzD7WCRF5LjsfkDHmnDHmRWPMXmOMyxgzG9gDNEu9rV2aPyUiDdyWhYvIBREpLSIlRGS2/Z7/sZ9XTOu8IvKiiHzl9jrC/hyC3D6/Sfbne1BE/s/tS7iGiCyxfyfHRWRqdt678i2a0NVVGWMSgBXA9fai64EoYFmqZVcrnQ8CugFVgUbAcAAR6QY8CXQBagKp667fwyrxVsMq/d4G3GGvWwJ0sJ+3B3a7xdPeXg/wEHCTvaw88A/wgX3+esBHwDB7XUkgzeSZFSJSBqgFbEq9zhgTD/wIDHVbPAhYYow5hvU/+RlQBagMXADez2Yok4EkrKunJsCNQEr9+8vAfKAE1nt+L5vnUD7E0YQuIp+KyDER2eih4yXb9ZdrRWSmJ46pACs5piTLdlgJPSrVsiVp7JfiXWPMIWPMSWAWl0r2g4DPjDEbjTHngBdTdrBLkkOAZ4wxscaYvcCbWMk3Jab2bud/ze21e0K/D3jOGHPATqYvAjfbpdibgdnGmKX2uv8Argw+iyftEvYpETmeeqWIBANfA58bY7amc4xv7PeW4l/2MowxJ4wxPxhjzhtjYoFX3N5XptlfKj2AR+0riGPAOLfzJmJ9aZQ3xsQZY5Zl9RzK9zhdQp+MVXLzlAvGmEj70ceDx83vlgJt7aqKcGPMDuAPrLr1a4AGXL2EfsTt+XmgiP28PFbdc4q/3Z6Xwqrq+TvV+gr28yVAO7v+PhD4Dmhj11sXA9ba21UBpqckYWALkAyUSX1++0vlxFXeB8BYY0xx+1HKfYVdHfQlkAA8eJVjLAIKiUhLO95IYLp9jEIi8rFdDXQG63MtnlJVkgVVsD6/w27v/WOgtL3+aUCAlSKyKXV1l8qbgpw8uTFmaeqGIxGpjnVJHI71z3/PVUo6Knf8iZUk7wF+BzDGnBGRQ/ayQ8aYPdk47mGgktvrym7Pj3OpFLnZbf1B+/w7ReQ8VpXKUjueI8AIYJkxJqWkvR+40xjze+qTi8hhoK7b60JY1S5ZJiICTML6ouhhjElMb1tjTLKIfIdV7XIU6yoh1l79BFAbaGmMOSIikcAarOSb2jmgkNvrsm7P9wPxQCljTFIaMRzB+t0hIm2BX0VkqTFmZ2ber/JNTpfQ0zIBeMgY0wyrfvXDLOwbKiLRIrJcRG7ySnT5kDHmAhANPI5V1ZJimb0su71bvgOGi0g9O5m+4HbOZHv9KyISJiJV7HN95bb/EqyScEr1yuJUrwHG28eoAhcbIPva674HeolIWxEpALxE9v8nPsL6cuhtf14Z+QYYDNxiP08RhlVvfsq++nkhjX1TrAWuF5HKIlIMeCZlhTHmMFYd+ZsiUtRuHK4uIu0BRGSgW2PrP4Ah4+om5eN8KqGLSBGgNTBNRNZiXSKWs9f1F5GNaTzmuR2iijGmOVad5Nt2aV95xhKsy3X3utYoe1m2EroxZg7wNrAQ2Gn/dPcQVil0t33eb4BPU8UU5nb+1K8B3gFmAvNFJBZYDrS0z78JGGkf9zBWYjtAFtlfFvdiVZ0ckUt91G9Jbx9jzAr7vZUH3Pv0vw0UxLpCWQ7MvcoxFgBTgfXAKmB2qk1uAwpgXeH8g/UFVs5edy2wQkTOYn0+jxhjdmfi7SofJk5PcGFXucw2xjQQkaLANmNMuQx2y8xxJ9vH/T6nx1JKqbzAp0roxpgzwB4RGQhWvaSINM7Mvnb/3RD7eSmgDZfqXpVSyu853W1xClaDW20ROSAid2HVKd4lIuuw+vH2vdox3NQFou39FgGjjTGa0JVS+Uamq1zsblPRwEFjTK9U64YDY7B7IADvG2MmejBOpZRSGchKt8VHsPrwFk1n/VRjzNX63iqllPKiTCV0u3tTT6y71h73xIlLlSplIiIiPHEopZTKN1atWnXcGBOe1rrMltDfxrqzLOwq2wwQkeuB7cBjxpj9qTcQkRFYN35QuXJloqOjM3l6pZRSYA1Ml966DBtFRaQXcMwYs+oqm80CIowxjYAFwOdpbWSMmWCMaW6MaR4enuYXjFJKqWzKTC+XNkAfEdkLfAt0ErchO+HigELx9suJpDFsqFJKKe/KMKEbY54xxlQ0xkRgjdS20Bhzq/s29gBJKfpgNZ4qpZTKRdkenEtEXgKijTEzgYdFpA/W2Msnsce7zqrExEQOHDhAXFxcdsNSeUBoaCgVK1YkODjY6VCU8iuO3frfvHlzk7pRdM+ePYSFhVGyZEmsweuUvzHGcOLECWJjY6latarT4SiV54jIKnvMqiv41K3/cXFxmsz9nIhQsmRJvQpTygt8KqEDmszzAf0dK+UdPpfQlVKesTvmLN+vOoDTI6qq3OPojEVKKc9LSHIxYeku3l24k4QkF4UKBNKjYY5HpFZ5gJbQU3n33XepW7cuIkKjRo1o2LAhrVu3Zt26dRe3uXDhAu3btyc5ORmXy8XDDz9MgwYNaNiwIddeey179lizsUVERNCuXbvLjh8ZGUmDBg0uW/boo49SoUIFXK5LE8ZMnjyZBx+8cmiciIgIjh+/Ym7iixYvXkyxYsWIjIykUaNG3HDDDRw7dizNY3711Vc0atSI+vXr07hxY+6++25OnToFwJAhQ9ixY0cmPzXlK9buP0Wf95cxdv52utQtQ+0yYYyes5X4pGSnQ1O5QBN6Kh9++CELFizg999/Z8mSJWzYsIH//Oc/jBgx4uI2n376Kf379ycwMJCpU6dy6NAh1q9fz4YNG5g+fTrFixe/uG1sbCz791ujIGzZcmX3fJfLxfTp06lUqRJLliy5Yn12tGvXjrVr17J+/XquvfZaPvjggyu2mTt3LuPGjWPOnDls2rSJ1atX07p1a44ePQrA/fffzxtvvOGReJT3nYtP4qVZm+n34e+cOp/IhGHN+OCWpjzXsy77Tp7nyz/TvVtc+RGfrXL536xNbD50xqPHrFe+KC/0rp/u+vvuu4/du3fTvXt37rzzTlq3bg1Aq1atOHDg0sxkX3/9Nd98Y00DefjwYcqVK0dAgPXdWLFixcuOOWjQIKZOncqTTz7JlClTGDp0KF9++eXF9YsXL6Z+/foMHjyYKVOm0LFjR4+9X2MMsbGx1KhR44p1r7zyCmPHjqVChQoABAYGcuedlyZ+b9euHcOHDycpKYmgIJ/9M1HA4m3HeG76Rg6eusCwVlV4ulttwkKtPv7X1wqnQ+1w3v1tBwOaVqRE4QIOR6u8SUvobsaPH0/58uVZtGgRjz322MXlkyZNonv37gAkJCSwe/duUkaKHDRoELNmzSIyMpInnniCNWvWXHbMAQMG8OOPPwIwa9Ysevfufdn6lCTfr18/fv75ZxIT050sPtOioqKIjIykcuXK/Prrr5cl6hSbNm2iadOm6R4jICCAGjVqXFbVpHzLyXMJPDZ1LcM/+4vQ4ACm3XcdL9/U4GIyT/Fsj7qcjU/ind+0Cs3f+WzR62ol6dy0aNEiJk2axLJl1tzIx48fv6xKpWLFimzbto2FCxeycOFCOnfuzLRp0+jcuTMAJUuWpESJEnz77bfUrVuXQoUKXdw3ISGBX375hbfeeouwsDBatmzJvHnz6NXrsvlDsqxdu3bMnm3NF/z666/z9NNPM378+HS337BhA8OGDSM2NpZXX32VwYMHA1C6dGkOHTpEs2Y6NI8vMcbw09pDvDR7M7FxiTzcuSYjO1YnJCgwze1rlQljSIvKfLX8b267rgrVwovkcsQqt2gJ/SrWr1/P3XffzU8//UTJkiUBKFiw4BU3xYSEhNC9e3fGjBnDs88+y4wZMy5bP3jwYEaOHMnQoUMvWz5v3jxOnTpFw4YNiYiIYNmyZUyZMsWj76FPnz4sXbr0iuX169dn9erVADRs2JC1a9fSvXt3Lly4cHGbuLg4ChYs6NF4VM7sP3me2z/7i0enrqXyNYWY/VA7Hu9SK91knuKxG2oREhTA6DlbcylS5QRN6OnYt28f/fv358svv6RWrVoXl5coUYLk5OSLSX316tUcOnQIsBo4169fT5UqVS47Vr9+/Xj66afp2rXrZcunTJnCxIkT2bt3L3v37mXPnj0sWLCA8+fPZynWlStXctttt6W5btmyZVSvXv2K5c888wxPPvnkZW0D7skcYPv27Vf0yFHOSHYZJi3bw43jlhK99yQv9q7HD/e3pnbZq01RcEl4WAgPdKzB/M1HWb77hJejVU7x2SoXp7300kucOHGCBx54AICgoKCLE3LceOONLFu27GKXwHvuuYf4eGv04BYtWlzR3TAsLIx///vfly07f/48c+fOvawqpHDhwrRt25ZZs2YBVjdD99L+8uXLAWjUqNHFRthBgwbRunXry0rSKXXoxhiKFSvGxIlXTu/ao0cPYmJi6N69O8nJyRQvXpwGDRpc/NI5evQoBQsWpGzZsln/8JRHbTl8hlE/rGfdgdN0rB3O//VrSIXiWb9yuqttVb5ZsY//+3kzM0e2JSBA79j1Nz41ONeWLVuoW7euI/FkxerVqxk3btxlvVWc9NRTTzFs2DAaNWrksWOOGzeOokWLctddd3nsmO7yyu/aSXGJyby/cCfjl+yiWMFg/tu7Hn0al8/R0Akz1hzk0alreWtQY/o3rZjxDsrnXG1wLi2hZ0PTpk3p2LEjycnJBAZeve4yN4wZM8bjxyxevDjDhg3z+HFV5qzYfYJnftzA7uPnGNC0Is/3rOuRLod9Gpfn09/3MGbeNro3KEfBAs7//SrP0Tr0bLrzzjt9Ipl7yx133KH9zx1wJi6RZ6dvYPCE5SQku/jyrha8Oaixx/qPBwQIz/esx+HTcUyM2u2RYyrfof+xSvmIuRuP8N+fNnL8bDz3tKvKY11qUaiA5/9FW1S9hm71y/LRkl0MblGJ0mGhHj+HcoaW0JVy2NEzcdz35Sru+2oVJYuEMGNkG57rWc8ryTzFqO51SEx2MW7Bdq+dQ+U+LaEr5RBjDN/+tZ9Xf9lCfJKLp7vV5p521QgO9H45K6JUYYa1imDyH3u4vXUEdcoW9fo5lfdpCV0pB+yOOcuQCct55scN1C9flHmPXs8DHWrkSjJP8XDnGoSFBvPKzzqnu7/QhJ5KXh8+N6NjhoeH06RJE2rWrEnXrl35448/Lq4fPnw4hQoVIjY29rLjiMjFcxYpYt02Hh0dTf369UlISABg165dVKtWjTNnzrBhwwaGDx+eYYz5UWKyiw8W7aTbO1FsPnyG0f0bMuWeVlQtVTjXYyleqAAPd65J1I7jLN52LNfPrzxPE3oq/jB87tWOOXjwYNasWcOOHTsYNWoU/fv3vyyuGjVq8NNPP108zsKFCy+OyOiuefPmtG/fnrFjxwIwcuRIXnnlFYoWLUrDhg05cOAA+/bt88j78Rfr9p+i93vLGDNvG53rlOa3x9szpEVlR6fkG9aqChElC/HqL1tISnZlvIPyaZmuQxeRQCAaOGiM6ZVqXQjwBdAMOAEMNsbszVFkc0bBkQ05OsQVyjaE7qPTXe0vw+dm9pgdO3ZkxIgRTJgwgXHjxgHWxBZTp07l1ltvZfHixbRp04Y5c+akuf+rr75KkyZNCAoKIikp6bKxanr37s23337L008/neP3k9edT0jirfnb+fT3PYSHhfDxsGZ0re8bd+AWCApgVPc63PfVar6LPsC/WlZ2OiSVA1kpoT8CpFfZdhfwjzGmBjAOeD2ngTnBX4bPzcoxmzZtytatlwZsqlWrFjExMfzzzz9MmTKFIUOGpLtv8eLFGTVqFM8888wVk2g0b96cqKioHL+XvG7p9hhuHLeUicv2MLRFZRY83t5nknmKrvXL0iLiGt5asI3YuJz//SnnZKqELiIVgZ7AK8DjaWzSF3jRfv498L6IiMnJuAJXKUnnprw2fG5Wj5nWr6h///58++23rFixgo8//viq55szZw5lypRh8+bN1K5d++LylKF386t/ziXw8s+b+XH1QaqVKsx3915Hi6rXOB1WmkSE53vVpc/7vzN+yS6e6lrH6ZBUNmW2yuVt4GkgvaHdKgD7AYwxSSJyGigJXNZ6JyIjgBEAlSv7/qVdyvC5c+bMydTwud27d6dMmTLMmDHjYkKHS8PnTp48+bL93IfPBWvAroIFC+YooWf1mGvWrLliTJXBgwfTrFkzbr/99otVSWmZPXs2p0+fZt68efTr14+uXbte/MLKz0Pv7j1+jgEf/cHpC4k81KkGIzvWIDTYt+8qblSxODdFlmdi1B5uaVmF8tkY/Es5L8MqFxHpBRwzxqzK6cmMMROMMc2NMc3Dw8NzejivyqvD52blmEuWLGHChAncc889ly2vUqUKr7zyysWRJtNy4cIFHn/8cT744AMaNmxI3759eeWVVy6uz89D774+dytxicnMfLAtT9xY2+eTeYqnulkl8zHztjkcicquzNShtwH6iMhe4Fugk4h8lWqbg0AlABEJAophNY7mWe7D50ZGRtK8+aXBzVKGzwU4duwYvXv3pkGDBjRq1IigoKB0h88tUODSeBwpw+f27Nnz4rK0hs+tWLHixUdKw2yjRo0uLnv88cfZt28fBQsWzNQxp06dSmRkJLVq1eLVV1/lhx9+SHPUw3vvvTfNcdRTvPzyy/Tr14969eoB8OKLLzJlyhR27LCmOVu0aNFlceQXq/f9w5yNR7jn+mrUK5+3btapULwgd7eryvQ1B1m3/5TT4ajsMMZk+gF0AGansXwkMN5+PgT4LqNjNWvWzKS2efPmK5b5olWrVplbb73V6TAuevLJJ826deucDuOiuLg407JlS5OYmJjuNnnld50VLpfLDPzoD9Ps5QXmbFz6792XxcYlmmYvzzcDP/rDuFwup8NRaQCiTTp5Ndv90EXkJRHpY7+cBJQUkZ1Yjaajsv0Nkwe4D5/rC8aMGePRsdBzat++fYwePTrfjdb425ZjrNx7kkduqEnhkLz53ouEBPFYl1qs3HuSeZuOOh2OyiKfm+CiTp06jt5oobzPGMPWrVv9aoKLZJeh29tLSXIZ5j92fa7ewu9pSckuur8TRWKyi/mPtadAUN59L/7oahNc+NRvKjQ0lBMnTqTZlU75B2MMJ06cIDTUv4Zs/WHVAXYcO8tTXWvn6WQOEBQYwHM967L3xHm+XP630+GoLPCp68KUhr+YmBinQ1FeFBoaesUdtXnZhYRk3lqwnchKxenewLduGsquDrVL065mKd79bQcDmlageCHPTLChvMunEnpwcDBVq1Z1OgylsuSzP/Zw5Ewc7wyJ9Kvqwud61qXHO1G8t3An/+lVz+lwVCbk7WtDpRz2z7kEPlq8i851StOyWkmnw/GoOmWLMqh5Jb74cy97j59zOhyVCZrQlcqBDxbt5Fx8Ek9388/b5R+/sRbBgQGMnrM1442V4zShK5VN+0+e54s//2ZA04rULpveqBh5W+mwUO5vX525m46wcs9Jp8NRGdCErlQ2vbVgOyJWKdaf3d2uGmWLhvLKz5txubQHmi/Lcwl93f5TDJu0Qof5VI7adOg0M9Ye5I42VSlXzL8HsipYIJCnutZm3YHTzFqff0fQzAvyXEJPNoZlO48zVgcQUg4aPWcrxQoGc3+H9Me78Sf9mlSgQYWivD7HGnhM+aY8l9CbVi7Bba2q8MXyv1n19z9Oh6PyoWU7jhO14zgPdqxBsYLBToeTKwIChOd71uPQ6TgmLdvjdDgqHXkuoYM1zGfZoqE88+N6EpJ0HkSVe1wuw+i5W6hQvCC3tqqS8Q5+pFW1knSpV4aPFu/i+Nl4p8NRaciTCb1ISBAv923A9qNn+XjJLqfDUfnIrPWH2HjwDE/cWCvPjHPuSc90r0NcYjLjFmx3OhSVhjyZ0AFuqFeGng3L8d7CneyKOet0OCofiE9KZuz8bdQtV5SbIis4HY4jqoUX4dZWVZiych/bj8Y6HY5KJc8mdIAX+tQjNDiAZ37coN2plNd9s2If+09eYFT3OgQE+M8t/ln1SOeaFAkJ4tVf0pszXjklTyf00mGhPNujLiv3nGRq9H6nw1F+LDYukfcW7qRNjZJcX7OU0+E4qkThAjzUqSaLt8UQtUMH0vMleTqhAwy+thItq17Dq79s4diZuIx3UCobPl6ym5PnEhjVra5fDcCVXbe1rkLlawrxys9bSNarY5+R5xO6iPBa/4bEJ7l4cdYmp8NRfujomTgmLttN78blaVixmNPh+ISQoEBGda/D1iOxTNOrY5+R5xM6WA01D3eqwS8bjrBgs06bpTzr7V+3k+wyPHVjbadD8SndG5SleZUSvLlgO+fik5wOR+EnCR1gxPXVqV0mjP/M2KjDAiiP2XnsLFP/2s8tLatQuWQhp8PxKSLCcz3rEhMb79Pdh0+fT+SjxbtoM3ohXcct5dfNR/12VjS/SegFggIYPaAhR2PjGKPDAigPeWPuVgoVCOKhTjWcDsUnNalcgt6NyzMhajeHT19wOpzL7Dl+jv/+tJFWr/3G63O3UqVkIRJdLu7+IpohE5azbv8pp0P0OJ+asSinmlQuwe3XRfD5n3vpG1mBZlVKOB2SysOi955k/uajPNGlFiWLhDgdjs96umtt5m06wph523hrUKSjsRhjWL77JJOW7ea3rccIDgigT2R57mxTlXrli5KY7OLbv/bzzq/b6fvB7/RuXJ6nbqztN1dfktGlh4iEAkuBEKwvgO+NMS+k2mY4MAY4aC963xgz8WrHbd68uYmOjs5m2Ok7G5/EjW8toUhoELMfaqczlqtsMcZw8/g/2X/yPIuf6kChAn5V9vG40XO2Mn7JLmY/1JYGFXK/4TghycXs9YeYGLWHzYfPcE3hAtzasjK3XleF0mFXTkh+Nj6Jj5fs4pOo3SS7DLddF8FDnWrkiblTRWSVMaZ5musykdAFKGyMOSsiwcAy4BFjzHK3bYYDzY0xD2Y2KG8ldIDfthzlrs+jebxLLR7uXNMr51D+bd6mI9z75Spe7deQf7Ws7HQ4Pu9MXCIdxiymVpkiTLmnVa517Tx5LoFvVvzNF3/+zbHYeGqWLsJdbatyU5MKmRqa4eiZON6av51pq/ZTJCSIkR1rcHvrCJ8e1uFqCT3D4quxpNxbH2w/fLpFoXPdMvRsVI73F+5k5zEdFkBlTVKyizfmbqVaeGEGNa/odDh5QtHQYB67oSbLd5/k1y3HvH6+ncfO8uz0DVz32m+Mnb+dOuWK8vmdLZj/2PUMaVE50wm5TNFQXr+5EXMeuZ5mVUrw2pytdH5zCTPWHMyTd59nqj5CRAJFZC1wDFhgjFmRxmYDRGS9iHwvIpXSOc4IEYkWkeiYGO/eYfZCb2tYgGd1WACVRdNWHWBXzDme7lqHoECtssusoS0qU6N0EV77ZQuJyZ4fBdUYQ9SOGIZ/tpIb3lrC96sO0K9JBeY/dj1f3NmC9rXCs31lULtsGJ/d0YJv7m5JicLBPDp1LX0+WMYfO497+F14V4ZVLpdtLFIcmA48ZIzZ6La8JHDWGBMvIvcCg40xna52LG9WuaSY+tc+/v3DBr1sVpl2PiGJDmMWU+maQnx/33V6V2gWLdx6lDsnR/Ni73oMb1PVI8eMS0xm5tpDfPr7HrYeiaVUkRBuu64Kt7Ss7JXGapfLMHPdIcbM28bBUxfoUDucZ7rX9Zl5Y3NUh57Gwf4LnDfGjE1nfSBw0hhz1ZaR3EjoxhiGfrKcTYfO8Nvj7Sld9MrGEaXcvb9wB2Pnb+f7+66jecQ1ToeT5xhjuHXSCjYfOsPipzrmaAKQ42fj+Wr533y1/G+On02gTtkw7mpblT6R5QkJ8n4dd1xiMl/8uZf3F+7kbHwSA5tV4rEutShbzNk8kqM6dBEJt0vmiEhBoAuwNdU25dxe9gF8Yhg2a1iARsQnuXhhpg4LoK7u5LkExi/ZTZd6ZTSZZ5OI8FyPepy6kMgHi3Zm6xjbjsTy9PfraD16IW//uoNGFYvz9d0tmfNIOwY2r5QryRwgNDiQEddXZ8lTHbmjTVV+XHOADmMX8eb8bZz10TtjM9MXqxzwuV3yDgC+M8bMFpGXgGhjzEzgYRHpAyQBJ4Hh3go4q6qWKswjnWsyZt425m86wo31yzodkvJR7y3cwfmEJP7dTW/xz4l65YsysFlFJv++l1szeYety2VYsiOGT5ftIWrHcUKDAxjYrCJ3tKlKjdJFciHq9JUoXID/9KrH7ddFMGb+Nt5buJNvVuzj0RtqMqRFZYJ9qJ0ly1UunpIbVS4pEpNd9H5vGafOJ7Lg8esJC80f80CqzNt34jyd31rMgKYVGT2gkdPh5HlHz8TRYcxiOtUpzQe3NE13u7jEZH5cfZBPf9/DzmNnKR0Wwu2tI/hXi8qUKOybfcLX7j/Fq79sYeWek1QrVZh/d6/DjfXK5Fp7S46qXPxBcGAAr/W3hgV4Y64OC6CuNHb+NgIDhMe61HI6FL9Qpmgo97avxs8bDrPq75NXrD92Jo4352/jutd+49npGwgNDmDc4MYs+3cnRnas4bPJHCCyUnGmjmjFxNuaIwL3frmKQR//yep9zk9any9K6Cn+N2sTk//Yy/f3XUezKlpHqiwbDpym9/vLGNmxOk91reN0OH7jfEISHccuplyxgkx/oDUiwqZDp5m0bA+z1h0iyWW4oW4Z7mpblZZVr8mTPYqSkl1Mjd7PuAU7OH42nh4Ny/J01zpElCrstXN6tJeLpziR0M/FJ3HjuKUUKhDI7Ifb5lrjSn6w7UgsJQoF58meRLdOXMGmQ6dZ8nRHimp1nEdNi97PU9+v59721Vi//zR/7j5BoQKBF+vHvZn4ctO5+CQmLN3NhKW7SXK5uKVlFR7uXJNrvHClke+rXFIUDgni/25qwI5jZxm/eLfT4fiFZJfhvd920P2dpdzw1hJmrjvkdEhZsnR7DMt2HuehTjU1mXvBgKYVqVeuKB8v2c3eE+d4pnsd/hzVmf/1beA3yRys3PJYl1oseaoDNzerxBd/7qX9G4v4aPEu4hKTcy2OfFVCT/HQlDXM23iEXx5pS43SvnGzQF4UExvPY1PXsmzncXo1KsfBUxdYs+8U/ZpU4H996/t8gnS5DL3eW8aZuER+e6K9XrF5yf6T59l6JJYOtcN9qkeIN+04GsvoOVv5besxyhcL5Ykba9OvSQWPTC6uJfRU/turHgULBPKMDguQbX/sOk6Pd6P4a+9JRvdvyHtDmzDt3ut49IaazFx3iO5vR7Fyz5WNYb7kp3UH2Xz4DE91ra3J3IsqXVOILvXK5JtkDlCzTBiThl/LlHtaUbJICE9MW0fP95Z5fVLt/PMJuwkPC+G5nnX5a+8/TPlrn9Ph5CnJLsM7v+7g1okrCAsNYsbINgxpURkRISgwgEdvqMW0+64jKFAYMuFPxszbSkKS58f1yKm4xGTGzttOgwpF6d2ovNPhKD91XfWS/DSyDe8MiSQ2LpFhk1Zy26cr2XzojFfOly8TOsDAZhVpXb0ko3/ZytEzcU6Hkycci41j2KQVjPt1O30jKzDrwbbULVf0iu2aVi7Bzw+3Y2CzSnywaBcDPvqDXTG+NerlV8v/5uCpC4zqVtcjl8FKpScgQOgbWYHfnmjP8z3rsm7/KWav905bU76sQ0+x9/g5ur69lI61SzN+WDNHY/F1f+w8zsPfriU2LpGX+zZgYPOKmepmNnfjYUb9uIG4xGSe71mPW1pWdrx72ukLibQfs4iGFYrx5V0tHY1F5T+nzycSGCgUCcnepClah56OiFKFeeSGmszddIR5m444HY5PSnYZxi3Yzi2TVlCsYBAzH2zLoGsrZTopd2tQjnmPXs+1Edfw/IyN3P15NMfPxns56qsbv2QXp84n8u9u2udc5b5ihYKzncwzkq8TOsA97apRt1xR/vvTRs7EJTodjk85FhvHrRNX8M5vO+jXpAIzH2ybrSFEyxQN5fM7WvBC73pE7TxOt7eXsnDrUS9EnLHDpy/w6bI93BRZ3pGp0pTypnyf0IMDAxjdvyExsfG8MXdrxjvkE8t2HKfHO1Gs2f8PY25uxFuDIimcg1JFQIBwR5uqzHqwLaWKhHDn5Gien7GBCwm510cXYNyC7RgDT9yoA3Ap/5PvEzpA40rFGd66Kl8t30f0Xt/uaudtyS7DWwu2M+zTFRQvVICZD7ZlYPM0J6DKltplw/jpwTbc0876vHu9F8XGg6c9dvyr2X40lu9XHWDYdVWodI1/zPKulDtN6LYnbqxFheIFGfXjBuKTcrfU6CuOnYnjlonLefe3HQxoWpGZD7ahVhnP33gVEhTIcz3r8fXdLTkXn8xNH/zOh4t3kuzlewLemLuVwgWCeLBjDa+eRymnaEK3FQ4J4v/6NWDnsbN8tHiX0+HkuqgdMfR4N4p1+08zdmBjxg5sTKEC3mm4SdGmRinmPtqOrvXL8sbcbQz9ZDkH/jnvlXOt3GNNXnxfh+o+PZKfUjmhCd1Nx9ql6dO4PB8u2sXOY7FOh5MrkpJdvDl/G7d9upIShQow88E23Nws92a6L16oAO//qwlvDmzM5kNn6P5OFD+tPejRcxhjeG3OFsoWDeVOD81zqZQv0oSeyn9716NQSCCjfvD/YQGOnonjlokreG/hTm5uWpGfHmxDTS9UsWRERBjQrCJzHmlHrTJhPPLtWh6esobTFzzT62jepiOs2XeKx7rUpGABvcVf+S9N6KmUKhLCcz3qEv33P3yz0n+HBVi6PYYe70Sx/sBp3hzYmDG5UMWSkUrXFGLqiFY80aUWP284TPe3l7J894kcHTMx2cUbc7dRs3QRBjTNvSsPpZygCT0NNzerSJsaJXl9zlaOnPavYQGSkl2MnbeN2z9bSakiIcx6qA0DcrGKJSNBgQE81LkmP9zfmpDgQIZ+spzRc7I/HszUv/az+/g5/t2tDkH5aHAolT/pX3gaRIRXbmpIQrKLF2ZudDocjzlyOo5/fbKC9xftZHDzSswY2cZnhw+OrFSc2Q+1Zci1lRi/ZBf9Pvw9y+0a5+KTePvXHVwbUYLOdUt7KVKlfIcm9HRElCrMozfUYt6mo8zdmPeHBVi87Rg93o1i46HTjBvcmNEDGvl8fXLhkCBe69+ICcOacfh0HD3fXcYXf+4ls+MPTYzaw/Gz8YzqXtfx8WOUyg0ZJnQRCRWRlSKyTkQ2icj/0tgmRESmishOEVkhIhFeiTaX3d2uap4fFiAp2cXrc7cy/LO/KB0WwswH29Kvie9UsWTGjfXLMvfRdrSqVpL//rSJOyf/RUzs1ceDOX42nglLd9GtflmaVSmRS5Eq5azMlNDjgU7GmMZAJNBNRFql2uYu4B9jTA1gHPC6R6N0SHBgAK8PaMjxs/G8PifvDQtw+PQFhn6ynI8W72Joi5QqliJOh5UtpcNCmXzHtfyvT33+2HWCbm8v5dfN6Y8H895vO4hLcvFUN73FX+UfGSZ0Y0kZzDrYfqS+5u0LfG4//x7oLH5yjduoYnHuaFOVr1fs4688NCzAom3H6PFOFJsPneGdIZG81r8RocG+XcWSERHh9tYRzH6oLWWKhnL3F9E8O30D5xOSLttu7/FzfL1iH4OvrUT18Lz5BaZUdmSqDl1EAkVkLXAMWGCMWZFqkwrAfgBjTBJwGiiZxnFGiEi0iETHxHh3KiZPeryLPSzAD+t9fliAxGQXo+ds5Y7P/qJM0VBmPtSWvpEVnA7Lo2qWCWP6yNbc274aU1buo9e7y1h/4NTF9WPmbyM4MIBHO9d0LkilHJCphG6MSTbGRAIVgRYi0iA7JzPGTDDGNDfGNA8PD8/OIRxROCSIV/o1YFfMOT5c5LvDAhw6dYEhE5Yzfsku/tWyMjNGtvHbEmpIUCDPdK/L13e35EJiMv0//IP3F+5g9b5/+Hn9Ye5pV5XSRUOdDlOpXJWlXi7GmFPAIqBbqlUHgUoAIhIEFANydkeIj+lQuzR9I8vz4eKd7Djqe8MCLNx6lB7vRrH1sFXF8mq/hnm+iiUzWlcvxdxHrqd7w3KMnb+dIROWU7JwAUa0r+50aErlusz0cgkXkeL284JAFyB1C+FM4Hb7+c3AQuPU3HZe9J9e9SgcEsSoH31nWIDEZBev/bKFOydHU65YQWb5YRVLRooVCubdIZG8PTiSsJAg/t29jtdmhFHKl2Xmr74c8LmIBGJ9AXxnjJktIi8B0caYmcAk4EsR2QmcBIZ4LWIHlSoSwvM96/HktHV8vXIfw1pV8ejx4xKTOXU+kX/OJ3DqfCKnzifwz/lETl2wXv9zLoFTF6zl1naJnL6QQGKy4ZaWlflPr3r5olSeFhHhpiYV6BtZXvucq3wrX08SnR3GGIZNWsna/af49fH2lC12ZT1tYrKLU3ay/SedRHzpeQKnL1g/4xLTv709JCiAEoUKULxQMMULBbs9L8C1ESXoVKeMN9+2UspHXG2SaE3o2fD3iXN0fXspNUuHEVGq8OXJ+XwisfFJ6e4bFCAXE3HxgvbPQsGUKOT+3E7WBQtQorD1Or+WvJVSl7taQteKxmyoUrIwL/auz9j524mNS6RYoQKULFKAGqWLUKyglYBLFA6++Nw9SRcJCdIqAaWUV2hCz6YhLSozpEVlp8NQSqmLdHAupZTyE5rQlVLKT2hCV0opP6EJXSml/IQmdKWU8hOa0JVSyk9oQldKKT+hCV0ppfyEJnSllPITmtCVUspPaEJXSik/oQldKaX8hCZ0pZTyE5rQlVLKT2hCV0opP6EJXSml/IQmdKWU8hOa0JVSyk9kmNBFpJKILBKRzSKySUQeSWObDiJyWkTW2o//eidcpZRS6cnMnKJJwBPGmNUiEgasEpEFxpjNqbaLMsb08nyISimlMiPDErox5rAxZrX9PBbYAlTwdmBKKaWyJkt16CISATQBVqSx+joRWScic0Skfjr7jxCRaBGJjomJyXq0Siml0pXphC4iRYAfgEeNMWdSrV4NVDHGNAbeA2akdQxjzARjTHNjTPPw8PBshqyUUiotmUroIhKMlcy/Nsb8mHq9MeaMMeas/fwXIFhESnk0UqWUUleVmV4uAkwCthhj3kpnm7L2dohIC/u4JzwZqFJKqavLTC+XNsAwYIOIrLWXPQtUBjDGjAduBu4XkSTgAjDEGGM8H65SSqn0ZJjQjTHLAMlgm/eB9z0VlFJKqazTO0WVUspPaEJXSik/oQldKaX8hCZ0pZTyE5rQlVLKT2hCV0opP6EJXSml/IQmdKWU8hOa0JVSyk9oQldKKT+hCV0ppfyEJnSllPITmtCVUspPaEJXSik/oQldKaX8hCZ0pZTyE5rQlVLKT2hCV0opP6EJXSml/IQmdKWU8hOa0JVSyk9kmNBFpJKILBKRzSKySUQeSWMbEZF3RWSniKwXkabeCVcppVR6gjKxTRLwhDFmtYiEAatEZIExZrPbNt2BmvajJfCR/VMppVQuybCEbow5bIxZbT+PBbYAFVJt1hf4wliWA8VFpJzHo1VKKZWuLNWhi0gE0ARYkWpVBWC/2+sDXJn0EZERIhItItExMTFZDFUppdTVZDqhi0gR4AfgUWPMmeyczBgzwRjT3BjTPDw8PDuHUEoplY5MJXQRCcZK5l8bY35MY5ODQCW31xXtZUoppXJJZnq5CDAJ2GKMeSudzWYCt9m9XVoBp40xhz0Yp1JKqQxkppdLG2AYsEFE1trLngUqAxhjxgO/AD2AncB54A6PR6pUXnFiF8wdBaVqQUQ7qHIdhBZzOiqVD4gxxpETN2/e3ERHRztybqW8Jj4WJnaBU3+DKwmSE0ACoGwjiGirCV7lmIisMsY0T2tdZkroSqnMMAZmPADHt8GtP0LlVnAgGvYug71RsHIC/Pm+leDLNb6U4Cu30gSvPEITulKeEvUmbJkJN/4fVO9oLavaznrwDCRegAN/2Ql+Gaz4GP54zy3Bt3NL8EUdfSsqb9IqF6U8Yft8+GYQNBgAAyaCSMb7JF6A/SsvJfgDf4Er0U7wkalK8JrgleVqVS6a0JXKqRO7YEJHKF4Z7poPBQpl7zgJ591K8FFWdY0rESQQykdenuBDwjz6FlTeoXXoSnlLfCx8+y8ICIAhX2U/mYO1b7X21gPsBO9Wgv/zQ/j9HbcEn1JF01ITvAI0oSuVfcbAjPvh+HarEbREhGePX6AQVOtgPcBK8PtXuCX4D+D3t+0E3yRVCb6IZ2NReYImdKWyK2osbJl1eSOoNxUoZJ0n5VwJ59zq4KOsHjQpCb5CUyvB1+pmJXiVL2gdulLZsX0efDMYGt4M/T/JXCOotyWcu7wEf3CV1Rd+8NdQt5fT0SkP0Tp0pTzpxC744R4o2wB6v+sbyRygQGGo3sl6gFW//0Vfq1oovA6UquFsfMrrdAo6pbLiYiNooFXyzUkjqLeFhMHAzyEgCL4bZpXglV/ThK5UZrlcMP0+qxF04GdQoorTEWWseCW4eRIc2wKzHrEacpXf0oSuVGZFvQlbZ0OXly/1PMkLqneCTs/Bhmmw8hOno1FepAldqczYPg8WvQINB8F1I52OJuvaPgG1usO8Z2Bf6gnHlL/QhK5URo7vhB/uthtB3/GdRtCsCAiAfuOhWCWYdjucPeZ0RMoLNKErdTVxZ+xG0CDfbwTNSMHiMPhLuHAKvr8TkpOcjkh5mCZ0pdLjclld/k7shIGT80YjaEbKNoRe46wbkRa+5HQ0ysM0oSuVnqixViPojS9fGl/FH0QOheZ3WuPCbJnldDTKgzShK5WWbXNh0atWI2irB5yOxvO6jYYKzWD6/VYbgfILmtCVSu34DvjxHqt6Iq82gmYkKAQGfQFBBWDqrXrTkZ/QhK6Uu5RG0MBgGJLHG0EzUqwiDJhkTZk382G96cgPaEJXKkXKnaAndlmNoMUrOx2R91XvCB2fg43fW3Oeqjwtw4QuIp+KyDER2ZjO+g4iclpE1tqP/3o+TKVywdIxsO1nazjcqtc7HU3uafu4fdPRs3rTUR6XmRL6ZKBbBttEGWMi7Yf2hVJ5z7Y5sPhVaDQYWt3vdDS5S2868hsZJnRjzFLgZC7EopQzju+AH0dAucb+2wiakYLFYfBXetNRHuepOvTrRGSdiMwRkfrpbSQiI0QkWkSiY2JiPHRqpXLAvRF08NcQXNDpiJxTtgH0ftu66ei3/zkdjcoGTyT01UAVY0xj4D1gRnobGmMmGGOaG2Oah4eHZ+9sSQlwVr8MlAe4XDD9XrsR9HNrqNn8rvEQaH4X/PEubJ7pdDSekxgHC16w7i1Iinc6Gq/JcUI3xpwxxpy1n/8CBItIqRxHlp6dC+DN2vDVAFg3FeLPeu1Uys8tfQO2/QJdX4Gq7ZyOxnd0ew0qNIcZD1jVUXnd8R0w8QZrvtUlr8OEDnB4ndNReUWOE7qIlBWxKh1FpIV9zBM5PW66SteDNo9AzHaYPgLG1oTv77Lu7EtO9NpplZ/Z+gssfg0aD4WW9zkdjW8JCoFBn1+66SgvF5rWfgMft4czB+Ff38G/psH5k/BJJ1j8ut/ljAwniRaRKUAHoBRwFHgBCAYwxowXkQeB+4Ek4ALwuDHmj4xOnONJol0ua0LcDd/Bpulw4R8oeA3U7weNBkGllvmzcSu3GQPHNsP2uVCoJNTrCwVLOB3V1cVst/6hS1aHO+fm73rzq9m9GL7sZ/1PDZiUt/6f4mPh5ydh/bdQpS0M+ASKlrfWnT8Jc562JvwoFwn9PobSdRwNNyuuNkl0hgndW3Kc0N0lJcCu32D9d1b3s6QL1k0hDQdaj9J1PXMeZXG54MBfsHUWbJkN/+y5tC4gGGp1tT73Wt0gONS5ONMSdxo+6WwVAEYs1nrzjES9Cb+9BN1eh1Z55Erm8DqYdof1d9l+FFz/pDUHbGqbZsDPj1tXIJ2etyYuSWs7H+P/Cd1dfCxs/dlK7rsXgXFBmYbQaCA0uBmKVfD8OfODpATYs9RK4lt/gXPHrORdrT3U6QW1e0DsYetz3/gDnD0CIUWhbh/rs49o5/w/i8sFU2+xZh+67SetN8+MlM9sx3wY/jNUbuV0ROkzxrrbdf7z1tXigIkQ0fbq+5w9BrMetW4oq9QKbvrQunLzYfkrobs7eww2/mhVyxxcBYj1C254c96oGnBa/FmrEXrLbOsfOv4MBBeGml2gbm/rZ2ixK/dzJVvJf8M0q6dEQiyElYMGA6ySe7nGzly+Lx5t1ZvnpdKmL7hwCj7pCIkXYMQSCCvjdERXOn8SfnrQSsw1u8JNH0Hhkpnb1xhYPxV+eRpcidDlJaunT4BvjoySfxO6uxO7YMP3VnI/sRMCC0DNG323asAp505YPT+2zoZdiyA53irt1O4OdXpbkyNn5bNKvGDVr6+fZn0puBKhVC2rnaPhQCgR4a13crmtv8C3Q61G0Js+ylv1wb7gyEarp0iFZtbVTWCQ0xFd8vef1hSBZ49aybjV/dn7/Z4+CDMfhF0Lrb/zPu/7ZJWcJnR3xsChNVbpceMP1h9BStVAw5utMTycrhrIbaf2W9VUW2fD379b1VTFKllVKXV7WZeinvgHPn8SNv9kVcvss9vNK7W0Env9/pkvUWWVNoJ6xrqpVs+y1g9Z4904zZUMy96CRa9ZbWY3fwoVmubsmMbAqs9g3vNWHug2GiL/5VMFAE3o6UmraqBIWatqoNFAqwXch36RHmMMxGy71Kh5eK21PLyOVZVSp5f3q0VO7bOvmKZZvWQCgqB6Z6vkXruH54atdW8EvXeJNWSsyr6fn4C/Jlpjqdfr61wcsUes4Rr2LLH+X3u9DaFFPXf8k3usfvj7/rAGLuv9js9UNWlCz4y0qgZK1rRKj40GwjXVnI4wZ1wuOLTamnJs62yr2gmsG0jq9rKqU0rVcCa2IxutqrAN31v9hYMLWzE1GgRVO2T/6sDlsm7r37nAqibIqIFMZSwpASb3gGNb4Z6FEF4r92PY+Sv8eK81KUePMdDkVu8UPlwuWPER/Po/q4DR8y1o0N/z58kiTehZdf4kbJlpJfe/l1nLKjS3Ekz9/lAkm8MW5LbkRNi7zErgW3+2eqEEBFmJrU4vqNPzUt9cX+ByWVU+G76zqmbiTkPhcLsxdZB1OZ2Vf9xFr8GS0dD9DWh5r/fizm9OH4SPr7faVu5ZCCFFcue8yYmw8GVrLtTS9eDmz3Kn/3jMNmuc/EOrrf//nm9CoWu8f950aELPidMHLlUNHN0IEmhNCtBwoNWoGlrMt+rcE85bffK3zLauOOJOQVBBqNHZqk6p1TVv9O5JireulNZ/Z3UzTI6Ha6rbV0yDMu5atvVnq3SujaDesXsJfHkT1LvJqrv29uf7z17rjvCD0dDsDmt4gtxsC0lOgt/HWXeXFrrGqoKp3T33zu9GE7qnHN18qWrg9P5LywOCrT+uoFCrB0hQwVQ/7cfFbTL6meoYQSFXbuP+D3T+pJX0ts6Gnb9ZN1aFFrd7pvSC6p3y9lRqF05ZVUUbvoM9UYCB8k2txN5gABQpffn22giaO6LeskZl7Dbau2PIb5phTZGHgT7vWneuOuXwephxv1W4i7wVur2adtddL9KE7mkpww4cWGmN4pZ0IZM/46y6+pSfJjn7Mbh/SZw9Zh0rrLxVjVK3F1RpYw0J62/OHLrU/fTIBuuKqVoHK7nX6Wn10Pmkk/UloI2g3mUMfHsL7JjnnZuOEi9YsyhFf2p1l7z509zr5no1SfHWIF/Lxln/c33ft67ac4kmdF+VnJSFLwO3n0nxl38xFA63SuLlm/jszRBecWyrfcU0zeo1E1TQahM49bc2guaWuNPW6IUJ5+HepZ7rCRKzzbp9/9gmaP0wdPqPNViYLzkQbQ+/vBOuvQe6/A8KFPb6aTWhK/9mjHXFtP47q93g+ieh+Z1OR5V/HN1kdQ2t0NS+6SgHV4bGwJqvrMGzggtZA2fVvMFzsXpawnmroXb5h1CiqjWVn5eHR9CErpTyrvXfwY/3wHUPWuPLZ0d8LMx+zLriimgH/T+BouU8G6e37ImCnx6wbtJr/RB0fM5rd59fLaHno+tzpZTXNBpkVTv8+b7ViJlVh9ZYXSE3/gAdn7dK+nklmYM10Nv9f0Cz4dZsTxPaW+8pl2lCV0p5RtdXoeK18NNIq6dRZhgDf34IE7tYbUPDf4b2T/lWV+DMCgmz5mS95YdLdygvei1XJ9HQhK6U8oygAtbcrEGhmZvp6PxJmDIU5j1jjdx53zKo0jp3YvWmmjfAA39aY0MtGW31ujq6OVdOrQldKeU5xSpY3QtP7LBGLkyvjW7v7/BRG+smuG6vw5BvHL370uMKloD+E2DwV1ZX2wntrW6Orhx0Vc4ETehKKc+q1h46/9eaGnL5R5evcyVbd1t+3stqNLxrgTU2vb/eyVu3N4xcYd2h/euL8Gk3ayhvL9GErpTyvDaPWvdGLPiPNV45wJnD8EVfWPyqNXvYvUuhfKSTUeaOwqVg0JdWr53j26wrk9VfeOVUmtCVUp4nYk3nVrwKTLsd1n0L49tYM4f1/dCqjggJczrK3CNi9QR6YLl1w1uwd4bi0H7oSinvOboZJnaGxPNQpoE1QqITQ+76kpScm81qphz1QxeRT0XkmIhsTGe9iMi7IrJTRNaLSA6nDFFK+Y0y9azqhvb/hrt/1WQOViL3UptBZqpcJgPdrrK+O1DTfowAPrrKtkqp/KbmDdDxWR31MhdkmNCNMUuBk1fZpC/whbEsB4qLSB66xUsppfyDJxpFKwBug4NzwF6mlFIqF+VqLxcRGSEi0SISHRMTk5unVkopv+eJhH4QqOT2uqK97ArGmAnGmObGmObh4XlkXk6llMojPJHQZwK32b1dWgGnjTGHPXBcpZRSWRCU0QYiMgXoAJQSkQPAC0AwgDFmPPAL0APYCZwH7vBWsEoppdKXYUI3xgzNYL0BRnosIqWUUtmit/4rpZSfcOzWfxGJAf7O5u6lgOMeDCev08/jcvp5XKKfxeX84fOoYoxJs1eJYwk9J0QkOr2xDPIj/Twup5/HJfpZXM7fPw+tclFKKT+hCV0ppfxEXk3oE5wOwMfo53E5/Twu0c/icn79eeTJOnSllFJXyqsldKWUUqloQldKKT+R5xK6iHQTkW32DEmjnI7HSRnNJpWfiEglEVkkIptFZJOIPOJ0TE4SkVARWSki6+zP439Ox+Q0EQkUkTUiMtvpWLwlTyV0EQkEPsCaJakeMFRE6jkblaMmc/XZpPKTJOAJY0w9oBUwMp//bcQDnYwxjYFIoJs9eF5+9giwxekgvClPJXSgBbDTGLPbGJMAfIs1Y1K+lInZpPINY8xhY8xq+3ks1j9uvp1oxZ5B7Kz9Mth+5NseECJSEegJTHQ6Fm/KawldZ0dSGRKRCKAJsMLhUBxlVzGsBY4BC4wx+fnzeBt4GnA5HIdX5bWErtRViUgR4AfgUWPMGafjcZIxJtkYE4k16UwLEWngcEiOEJFewDFjzCqnY/G2vJbQMz07ksp/RCQYK5l/bYz50el4fIUx5hSwiPzb3tIG6CMie7GqaTuJyFfOhuQdeS2h/wXUFJGqIlIAGII1Y5LK50REgEnAFmPMW07H4zQRCReR4vbzgkAXYKujQTnEGPOMMaaiMSYCK2csNMbc6nBYXpGnEroxJgl4EJiH1ej1nTFmk7NROceeTepPoLaIHBCRu5yOyUFtgGFYpa+19qOH00E5qBywSETWYxWEFhhj/La7nrLorf9KKeUn8lQJXSmlVPo0oSullJ/QhK6UUn5CE7pSSvkJTehKKZULsjKYnoiMc+uttV1ETmXmHJrQVaaISLL9x7VRRKaJSKEs7l9eRL63n0e6dykUkT6+OnKmiEToaJbKQyaTyZu7jDGPGWMi7Tt93wMydaOcJnSVWRfsP7AGQAJwX1Z2NsYcMsbcbL+MBHq4rZtpjBntsUgdICJBTsegfFtag+mJSHURmSsiq0QkSkTqpLHrUGBKZs6hCV1lRxRQQ0SuEZEZIrJeRJaLSCMAEWnvdrm4RkTCUkq69h2+LwGD7fWDRWS4iLxv7xshIgvtY/4mIpXt5ZNF5F0R+UNEdovIzWkFJiKP2+fZKCKPuh1zi4h8Yo8NPt++ezL1vmVEZLo9hvg6EWltrwpMa18RWSwib4tINPCIiHS23+8G+/I6xN5ur4i8Zr/faBFpKiLzRGSXiNxnbyMiMsaOe4OIDLaXlxORpW5XR+3s5UPt7TaKyOtu7+GsfZxNIvKriLSw49wtIn3sbQLtbf6yP+d7c/bnoHJgAvCQMaYZ8CTwoftKEakCVAUWZupoxhh96CPDB3DW/hkE/ATcj3Up+IK9vBOw1n4+C2hjPy9i7xMBbLSXDQfedzv2xdf2vrfbz+8EZtjPJwPTsAoh9bCGUU4dYzNgA1DYPu8mrFEXI7DGS4+0t/sOuDWN/adiDeoFEAgUu9q+wGLgQ/t5KNZIoLXs11+4HWsvcL/9fBywHggDwoGj9vIBwAL7vGWAfVh3ez4BPOcWUxhQ3l4fbn+2C4Gb7G0M0N1+Ph2YjzV0bmO3388I4Hn7eQgQDVR1+m8sPzxS/R8UAS4Aa90eW1Jt/2/gvcweX0voKrMKijUUazRWMpkEtAW+BDDGLARKikhR4HfgLRF5GChurCEbMus64Bv7+Zf2OVLMMMa4jDGbsZJeam2B6caYc8YaC/xHoJ29bo8xZq39fBXWP1ZqnYCP7PeTbIw5nYl9p9o/a9vbbbdffw5c77ZdyphDG4AVxphYY0wMEC/WmCttgSn2eY8CS4BrsW7bv0NEXgQaGmus92uBxcaYGPuz/drtXAnAXLdzLTHGJNrPU+K+EbjN/n2uAEoCNdP4PJR3BQCnjF1Xbj/qptpmCJmsbkk5oFKZccHtj+4hY00wkiZj1YffDRQEfk+nXjA74t2eSw72TcYq2Xpi33NZPIYr1fFcV4vFWPWu12ONKjpZRG7L4DyJxi7auZ/LGON+HsG6zE/5fVY1xszP5PtQHmKs4Z33iMhAuFjt1jhlvf1/UwJrvKZM0YSuciIKuAVARDoAx40xZ0SkujFmgzHmdawSZuqEHotVdZCWP7BKJdjHjspiPDeJSCERKQz0y+L+v2FVJaXUMxfLwr7bgAgRqWG/HoZVys6sKKx2hUARCcdK4ivtOtSjxphPsGbbaQqsBNqLSCmxpmUcmsVzzQPuF2u4YUSklv15KS+StAfTuwW4S0TWYVURus/ANgT41u0LOkPaMq9y4kXgU7FG9DsP3G4vf1REOmKVEDcBc7Dqg1MsAkbZl/yvpTrmQ8BnIvIUEAPckdlgjDGrRWQyVsIDmGiMWSPWDEaZ8Qgwwf5HS8ZK7oczee44EbkDmCZWj5e/gPGZjR2rvvs6YB1WPfjTxpgjInI78JSIJAJngduMMYfF6ua5CKu0/bMx5qcsnGsiVvXLahERrM/5pizsr7LBGDM0nVVpdmU0xryY1XPoaItKKeUntMpFKaX8hCZ0pZTyE5rQlVLKT2hCV0opP6EJXSml/IQmdKWU8hOa0JVSyk/8P29I+RLq8sEGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "breakpoints = [5*i*1e6 for i in range(0, 10 + 1)]\n",
    "f2_small_big = ts.f2(sample_sets=[ts.samples(0),\n",
    "                                  ts.samples(1)],\n",
    "                                windows=breakpoints)\n",
    "f2_small_admix = ts.f2(sample_sets=[ts.samples(0),\n",
    "                                  ts.samples(2)],\n",
    "                                  windows=breakpoints)\n",
    "\n",
    "# Plot.\n",
    "names_to_plot = ['f2(SMALL,BIG)', 'f2(SMALL,ADMIX)']\n",
    "lines = plt.plot(\n",
    "    breakpoints[:-1], np.transpose([f2_small_big, f2_small_admix]))\n",
    "plt.legend(lines, names_to_plot)\n",
    "plt.xticks()\n",
    "plt.xlabel(\"Position on chromosome\")\n",
    "plt.title(\"Windowed F2 values\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}