1.Introduction_to_msprime_SOLUTION.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Practical 1. Recap of the neutral coalescent and an intro to coalescence simulations with msprime\n",
    "---\n",
    "\n",
    "At the end of this practical you should:\n",
    "\n",
    "- Understand what a coalescent simulation is and what the relevant data structures are.\n",
    "- Know how to run and analyse coalescent simulations with [`msprime`](https://tskit.dev/msprime/docs/stable/intro.html) and [`tskit`](https://tskit.dev/tskit/docs/stable/python-api.html#trees-and-tree-sequences).\n",
    "- Understand the extent to which the variance of the coalescent depends on the two major axis of sampling: # of loci, # of individuals ([Felsenstein 2004](https://academic.oup.com/mbe/article-lookup/doi/10.1093/molbev/msj079)).\n",
    "- Understand why it is natural (and helpful) to treat mutations separately from ancestry."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 Getting started"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing python tools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> We will generally import python tools in the first code cell of a jupyter practical notebook. \n",
    "    \n",
    "There are several tools we will use throughout the course: <code>msprime</code> is a powerful and well documented coalescent simulator. In fact, this notebook is based largely on the \"Running Simulations\" section in the <a href=\"https://tskit.dev/msprime/docs/stable/mutations.html\"><code>msprime</code> manual</a>.\n",
    "\n",
    "<a href=\"https://numpy.org/\"><code>numpy</code></a> is an open source python library that is widely used for scientific computing.  If you are new to `python` or need a refresher about basic <code>python</code> concepts, complete Python-novice, please read sections 1. - 4.3 of the official <a href=\"https://docs.python.org/3/tutorial/index.html\">Python documentation</a>.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import msprime\n",
    "import tskit\n",
    "from IPython.display import SVG, set_matplotlib_formats\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "from collections import Counter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 Build your own simulator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Population genetic models are algorithmic by nature: forwards in time they specify how a population evolves from one generation to the next, backwards in time models are framed in terms of the probabilities of the next event in the ancestry of a sample. It is unsurprising that a large number of population genetic simulators (both forward and backwards) exist. Simulating under the neutral coalescent, the simplest backwards in time model, is straightforward in principle.\n",
    "\n",
    "In fact, if we're only interested in coalescence times, a coalescent simulation is just a set of random draws of coalescence times from an exponential distribution. We can do this with the `random` method in `numpy`. Let's draw 1,000 pairwise coalescence times for a population (of diploid individuals) with an effective size of $N_e= 10,000$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_replicates = 1_000\n",
    "pair_times_array = np.random.exponential(20_000, num_replicates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the result which we have assigned to a variable `pair_times_array` is a numpy array of floats:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(pair_times_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use python indexing to examine the first 10 simulation replicates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([17309.47701672, 15514.09887427,  3858.06923658, 25025.88505369,\n",
       "       31922.7415045 ,  7958.93119912,   235.13348422,  9898.93960206,\n",
       "        2903.52506035,  9562.77298053])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pair_times_array[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> We can use the <a href=\"https://seaborn.pydata.org/generated/seaborn.displot.html#seaborn.displot/\"><code>displot()</code></a> plot method to plot the distribution of pairwise coalescence times. Unsurprisingly (given that we've drawn from an exponential distribution), the pairwise coalescence times are exponentially distributed:\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFuCAYAAAC/a8I8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVvklEQVR4nO3de5AlZX3G8efZCywKBAgj2ctsuERRyijIBGWxDKLRjRdQgwIliooFxMJLjGWktCw0ZaJijPe4qyJeEFeuApYoAivqysKuLriwoNyKXRbdQcLNGOKyv/zR78hhnJk9c+nz6575fqq65pw+3f3+9tTsU++83f22I0IAgN6blV0AAMxUBDAAJCGAASAJAQwASQhgAEhCAANAksYFsO0zbW+xvX6Kjveo7XVluXgqjgkAU8FNuw7Y9vMkPSzpqxHx9Ck43sMRsfPkKwOAqdW4HnBEXC3pvs51tvezfZnttbZ/ZPupSeUBwJRpXACPYrmkt0bEwZLeJelz49h3nu01tq+x/YpaqgOACZiTXcD22N5Z0hJJ59oeWr1j+exVkj44wm53R8SLy+vFEbHZ9r6SrrT9i4i4re66AWB7Gh/Aqnrp90fEgcM/iIgLJF0w1s4Rsbn8vN32SkkHSSKAAaRr/BBERDwo6Q7br5YkV57Zzb62d7c91FveU9Jhkm6qrVgAGIfGBbDtcyT9VNL+tjfZPlHSayWdaPt6STdKOqrLwz1N0pqy31WSPhwRBDCARmjcZWgAMFM0rgcMADNFo07CLV26NC677LLsMmpzzLKfSpJWnHxociUAeswjrWxUD/jee+/NLgEAeqZRAQwAMwkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSTIsAXti/WLa7Whb2L84uFwAkNWw+4InavGmjjlm2qqttV5y8pOZqAKA706IHDABtRAADQBICGACSEMAAkIQABoAkBDAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJLUHsO3Ztn9u+9K62wKANulFD/jtkjb0oB0AaJVaA9j2IkkvlfTFOtsBgDaquwf8CUnvlrRttA1sn2R7je01g4ODNZcDAM1RWwDbfpmkLRGxdqztImJ5RAxExEBfX19d5QBA49TZAz5M0pG275T0TUlH2P56je0BQKvUFsARcVpELIqIvSUdK+nKiDi+rvYAoG24DhgAkszpRSMRsVLSyl60BQBtQQ8YAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkGTmBfCsObLd1bKwf3F2tQCmsZ48E65Rtm3VMctWdbXpipOX1FwMgJls5vWAAaAhCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkIQABoAkBDAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkIQABoAktQWw7Xm2r7V9ve0bbX+grrYAoI3m1HjsRyQdEREP254r6ce2vxsR19TYJgC0Rm0BHBEh6eHydm5Zoq72AKBtah0Dtj3b9jpJWyRdHhGrR9jmJNtrbK8ZHByssxwAaJRaAzgiHo2IAyUtknSI7aePsM3yiBiIiIG+vr46ywGARunJVRARcb+klZKW9qI9AGiDOq+C6LO9W3m9k6QXSrq5rvYAoG3qvApivqSv2J6tKui/FRGX1tgeALRKnVdB3CDpoLqODwBtx51wAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkIQABoAkBDAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYDHMmuObHe9LOxfnF0xgBap86nI7bdtq45ZtqrrzVecvKTGYgBMN/SAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIElXAWz7im7WAQC6N+atyLbnSXqCpD1t7y7J5aNdJS2ouTYAmNa2NxfEyZLeoSps1+qxAH5Q0mfrKwsApr8xAzgiPinpk7bfGhGf7lFNADAjdDUbWkR82vYSSXt37hMRX62pLgCY9roKYNtfk7SfpHWSHi2rQxIBDAAT1O18wAOSDoiIqLMYAJhJur0OeL2kv6izEACYabrtAe8p6Sbb10p6ZGhlRBxZS1UAMAN0G8Cn11kEAMxE3V4F8cO6CwGAmabbqyAeUnXVgyTtIGmupN9FxK51FQYA0123PeBdOt/bfoWkQ+ooCABmignNhhYRF0k6YmpLAYCZpdshiFd1vJ2l6rpgrgkebtYc2R71472O+3dJkk9ZogWL+nX3xrt6VRmABur2KoiXd7zeKulOSUdNeTVtt22rjlm2atSPb9z8oCTp8GWrtOLkJb2qCkBDdTsG/Ma6CwGAmabbCdkX2b7Q9hbbv7F9vu1FdRcHANNZtyfhvizpYlXzAi+UdElZBwCYoG4DuC8ivhwRW8tylqS+GusCgGmv2wC+1/bxtmeX5XhJv62zMACY7roN4DdJeo2kX0u6R9LRkjgxBwCT0O1laP8q6YSI+G9Jsr2HpI+pCmYAwAR02wN+xlD4SlJE3CfpoHpKAoCZodsAnlUeSy/pjz3gbnvPAIARdBui/yFple3zVN2C/BpJH6qtKgCYAbq9E+6rtteomoDHkl4VETfVWhkATHNdDyOUwCV0AWCKTGg6SgDA5NUWwLb7bV9le4PtG22/va62AKCN6rySYaukf46In9neRdJa25czdgwAldp6wBFxT0T8rLx+SNIGVRP5AADUozFg23urunFj9QifnWR7je01g4ODvSgHABqh9gC2vbOk8yW9IyIeHP55RCyPiIGIGOjrY4I1ADNHrQFse66q8D07Ii6osy0AaJs6r4KwpC9J2hARH6+rHQBoqzp7wIdJep2kI2yvK8tLamwPAFqltsvQIuLHqm5bBgCMgDvhACAJAQwASQhgAEhCAANAEgIYAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkI4Cyz5sh2V8vC/sXZ1QKoQW3PhMN2bNuqY5at6mrTFScvqbkYABnoAQNAEgIYAJIQwACQhAAGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICuA3G8fgiHmEEtAePJGqDcTy+SOIRRkBb0AMGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIQgADQBICGACSEMAAkIQAno7G8QQNnp4B5OGJGNPROJ6gwdMzgDz0gAEgSW0BbPtM21tsr6+rDQBoszp7wGdJWlrj8QGg1WoL4Ii4WtJ9dR0fANoufQzY9km219heMzg4mF0OAPRMegBHxPKIGIiIgb6+vuxyAKBn0gMYAGYqAhgAktR5Gdo5kn4qaX/bm2yfWFdbANBGtd0JFxHH1XVsAJgOGIIAgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQE80/EIeyANj6Wf6XiEPZCGHjAAJCGAASAJAQwASQhgAEhCAANAEgIYAJIQwACQhABGbRb2L+YmD2AM3IiB2mzetJGbPIAx0AMGgCQEMAAkIYABIAkBDABJCGAASEIAA0ASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYABIwmxo6N6sObKdXQUwbRDA6N62rV1PLykxxSSwPQxBAEASAhgAkhDAAJCEAAaAJAQwACQhgAEgCQEMAEkIYLTOwv7Fst3VsrB/cXa5wKi4EQOts3nTxq5vCOFmEDQZAYxm4DZnzEAEMJphHLc506vFdMEYMAAkIYABIAkBjOmtjC1z1QSaiDFgTG9MoYkGowcMAEkIYABIQgADQBICGOg0jpN2nLDDZHESDujUwhtCFvYv1uZNG7vadsGift298a5GHBsEMNB6dc6Nwbwb9WIIAgCS1BrAtpfavsX2rbbfU2dbQM+NY7x4zg7zmnFDyDhvTGmK8UxBOp7vOnscv7YhCNuzJX1W0t9J2iTpOtsXR8RNdbUJ9NQ4x4vHdUPIPz6vngCs88aUccxoN97x4vEOhXS97Ti+5zrGuOscAz5E0q0Rcbsk2f6mpKMkEcDA9rTwZCA1j58jYsoPKkm2j5a0NCLeXN6/TtKzI+LUYdudJOmk8nZ/SbdMoLk9Jd07iXIzUHNvtLFmqZ11U/Po7o2IpcNX1tkDHqlf/ydpHxHLJS2fVEP2mogYmMwxeo2ae6ONNUvtrJuax6/Ok3CbJPV3vF8kaXON7QFAq9QZwNdJerLtfWzvIOlYSRfX2B4AtEptQxARsdX2qZK+J2m2pDMj4saampvUEEYSau6NNtYstbNuah6n2k7CAQDGxp1wAJCEAAaAJK0O4OxbnW33277K9gbbN9p+e1m/h+3Lbf+q/Ny9Y5/TSr232H5xx/qDbf+ifPYpl9tzbO9oe0VZv9r23lNU+2zbP7d9aRtqtr2b7fNs31y+70NbUPM/ld+L9bbPsT2viTXbPtP2FtvrO9b1pE7bJ5Q2fmX7hEnWfEb5/bjB9oW2d2tSzSOKiFYuqk7s3SZpX0k7SLpe0gE9rmG+pGeV17tI+qWkAyR9VNJ7yvr3SPpIeX1AqXNHSfuU+meXz66VdKiq66e/K+nvy/q3SPp8eX2spBVTVPs7JX1D0qXlfaNrlvQVSW8ur3eQtFuTa5a0UNIdknYq778l6Q1NrFnS8yQ9S9L6jnW11ylpD0m3l5+7l9e7T6LmF0maU15/pGk1j/jvmIr/zBlL+dK+1/H+NEmnJdf0bVVzX9wiaX5ZN1/SLSPVqOoKkUPLNjd3rD9O0rLObcrrOaru2vEk61wk6QpJR+ixAG5szZJ2VRVmHra+yTUvlLSx/EedI+nSEhCNrFnS3np8mNVeZ+c25bNlko6baM3DPnulpLObVvPwpc1DEEO/4EM2lXUpyp8oB0laLWmviLhHksrPJ5XNRqt5YXk9fP3j9omIrZIekPTnkyz3E5LeLWlbx7om17yvpEFJXy7DJl+0/cQm1xwRd0v6mKS7JN0j6YGI+H6Tax6mF3XW+X/4Tap6tI2uuc0B3NWtzr1ge2dJ50t6R0Q8ONamI6yLMdaPtc+E2H6ZpC0RsbbbXUZpv2c1q+qBPEvSf0XEQZJ+p+rP4tGk11zGTI9S9SfvAklPtH38WLuM0n4vv+duTGWdtdRv+72Stko6exLt96TmNgdwI251tj1XVfieHREXlNW/sT2/fD5f0payfrSaN5XXw9c/bh/bcyT9maT7JlHyYZKOtH2npG9KOsL21xte8yZJmyJidXl/nqpAbnLNL5R0R0QMRsQfJF0gaUnDa+7Uizqn/P9wOSn2MkmvjTJG0OSa2xzA6bc6lzOmX5K0ISI+3vHRxZKGzo6eoGpseGj9seUM6z6Snizp2vIn3kO2n1OO+fph+wwd62hJV3b8Yo1bRJwWEYsiYm9V39mVEXF8w2v+taSNtvcvq16galrTxtasaujhObafUNp6gaQNDa+5Uy/q/J6kF9nevfzF8KKybkJsL5X0L5KOjIj/GfZvaWTNU3YCKmOR9BJVVx7cJum9Ce0/V9WfHzdIWleWl6gaK7pC0q/Kzz069nlvqfcWlTOuZf2ApPXls8/osbsU50k6V9Ktqs7Y7juF9R+ux07CNbpmSQdKWlO+64tUnYFues0fkHRzae9rqs7CN65mSeeoGqf+g6oe3om9qlPVWO2tZXnjJGu+VdX47LqyfL5JNY+0cCsyACRp8xAEALQaAQwASQhgAEhCAANAEgIYAJIQwBiVqxnI3tLxfoHt83rY/nw/NlvbgO1P9artqWb7DbYXdLz/ou0DpvD4f237rKk6HnqDy9AwqjK/xaUR8fSk9s+Q9OOI+PZ2N24A27Mj4tFRPlsp6V0RsabG9n8g6U0RcVddbWBq0QPGWD4saT/b68pcq3sPzb9aenQX2b7E9h22T7X9zjJZzjW29yjb7Wf7Mttrbf/I9lPL+le7mif3ettXj9L+P0i6rGx/eEdv+HRX88GutH277beNtLPtE23/smz3BdufKev7bJ9v+7qyHLa949o+3va15btYZnt2Wf+w7Q/aXi3pUNvvL8dcb3u5K0eruuD/7LL/TqWNgXKM41zNSbve9kc62nzY9ofKd3SN7b22891douruRrTFVN1VxTL9Fv3pFIV/fK9qbttbVc2D3KdqtqhTymf/qWpiIqm6i+rJ5fWzVd3SKUm/kLSwvN5thLb3kbS24/3heuyuvdMlrVJ1Z9mekn4rae6w/RdIulPVdJBzJf1I0mfKZ9+Q9NzyerGqW8lHPa6kp6kKt7llu89Jen15HZJe09Fu5x1jX5P08vJ6paSBjs9WqgrlBapuW+5TNenQlZJe0XHsof0/Kul9Y313qub5uCT794al+6W2pyJjRrgqIh5SdT/9A6pCSqoC4hmuZolbIunc6lZ7SVW4SdJPJJ1l+1uqJqoZbr6qKShH852IeETSI7a3SNpLj59a8BBJP4yI+yTJ9rmSnlI+e6GkAzpq2tX2LmMc9wWSDpZ0XdlnJz02Oc2jqiZjGvJ82++W9ARV4X9jx/cykr+RtDIiBkudZ6uabPwiSf+nah5hSVqraq5pafTvbouqQEdLEMCYjEc6Xm/reL9N1e/WLEn3R8SBw3eMiFNsP1vSSyWts31gRPy2Y5Pfq7ofv5u2H9Wf/i6PNG3gkFmqJtv+/eN2qMJ1pONa0lci4rQRjvW/UcZ9bc9T1TseiIiNtk/fzr9he3X+IUrXtqOWsb67eaq+N7QEY8AYy0OqhhgmJKq5ke+w/Wqpmj3O9jPL6/0iYnVEvF/V0wb6h+3+S1VDHhN1raS/dTVr1RxV48lDvi/p1KE3tg/czrGukHS07SeV7few/ZcjbDcUtveW3v/RHZ+N9l2uLnXuWcaVj5P0w7GKGeO7e4qqiWXQEgQwRlV6VT8pJ3zOmOBhXivpRNvXq/pz/Kiy/oyhE0+Srlb1zK7Otn8n6TbbfzXB2u+W9G+qAu4HqqavfKB8/DZJA64e3niTpFO2c6ybJL1P0vdt3yDpclVDJMO3u1/SF1QNwVykasrUIWdJ+vzQSbiOfe5R9cicq1R9Bz+L7V/1Mdp393xJ39nOvmgQLkNDY9l+paSDI+J9E9x/54h4uPSAL5R0ZkRcOKVFNoTtHVX1nJ8b1SN00AL0gNFYJSzvnMQhTre9TtWf5Xeo6pVOV4tVPcWY8G0ResAAkIQeMAAkIYABIAkBDABJCGAASEIAA0CS/wfk/BaLj2n3EQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = sns.displot(x = pair_times_array, stat = \"density\", binwidth = 5_000)\n",
    "plt.xlabel(\"times (in generations)\");\n",
    "plt.axvline(x= np.mean(pair_times_array)) \n",
    "plt.ylabel(\"count\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 1:</b> Find the <a href=\"https://numpy.org/doc/stable/\"><code>numpy<code></a> methods for computing the mean and standard deviation of an array and check how well the mean of the 1_000 pairwise coalescence times simulated above matches the analytic expectation.\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean pairwise coalescence time is 19036.529880950442\n",
      "The standard deviation of the pairwise coalescence time is 18513.07885092935\n"
     ]
    }
   ],
   "source": [
    "#SOL\n",
    "print(\"The mean pairwise coalescence time is\", np.mean(pair_times_array))\n",
    "print(\"The standard deviation of the pairwise coalescence time is\", pair_times_array.std())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>PopGen:</b> A basic property of the coalescent is that the rate is proportional to the number of pairs of lineages, and so changes after each coalescent event. When there are  $i$ lineages the rate of coalescence is $\\lambda = \\binom i 2 = \\frac{i (i-1)}{2}$.\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> Writing your own function is easy! The code chunk below defines a coalescence time simulator as a function. <code>coaltimes</code> has three arguments: <code>samplesize</code>,  <code>Ne</code> and  <code>replicates</code>, i.e. for a sample of size $n$ and a given $N_e$, <code>coaltimes</code> uses <a href=\"https://numpy.org/doc/stable/reference/random/generated/numpy.random.exponential.html\"><code>np.random.exponential()</code></a> to draw $n-1$ random coalescence times:\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def coaltimes(samplesize, Ne, replicates):\n",
    "    list_of_coal_intervals = [np.random.exponential(2* Ne/int(i*(i-1)/2), replicates) for i in range(2,samplesize+1)]\n",
    "    return np.asarray(list_of_coal_intervals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 2:</b> Which <code>numpy</code> functions and methods does <code>coaltimes</code> use? The waiting time for most other possible processes that affect the ancestry of a sample (e.g. mutation, migration, recombination events) depend on the number of lineages that are present (rather than the number of possible pairs). Copy <code>coaltimes</code> in the cell below and modify it to define a new function that draws waiting times at a rate proportional to the number of lineages (rather than the number of pairs).\n",
    "    \n",
    "<b>Hint:</b> Read the docs of <a href=\"https://numpy.org/doc/stable/reference/random/generated/numpy.random.exponential.html\"><code>np.random.exponential()</code></a> to understand what input parameters this function takes.\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL\n",
    "def coaltimes(samplesize, Ne, replicates):\n",
    "    list_of_coal_intervals = [np.random.exponential(2* Ne/int(i), replicates) for i in range(2,samplesize+1)]\n",
    "    return np.asarray(list_of_coal_intervals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 3:</b>  Use <code>coaltimes</code> to simulate 1,000 replicate sets of coalescence times for a sample of $n=9$ sequences (don't forget to assign the output to a new, helpful variable name).\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL\n",
    "n9_times_array = coaltimes(9, 10_000, num_replicates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 4:</b> Use the funtion <code>type()</code> and the method <code>shape</code> to check what form the output of <code>coaltimes</code> takes.\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#SOL\n",
    "type(n9_times_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 1000)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#SOL\n",
    "n9_times_array.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 5:</b> Compute the $t_{mrca}$ of each simulation replicate (as the sum of internode intervals) and use <code>sns.displot()</code> to visualize the distribution of $t_{mrca}$ arcoss replicates. What is the mean coalescence time for each internode interval? \n",
    "    \n",
    "<b>Hint: </b> Assign the $t_{mrca}$ to a new variable <code>tmrca</code>\n",
    "       </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10374.01597845,  6374.91800691,  4935.44026639,  4131.43502186,\n",
       "        3276.71178545,  2764.094745  ,  2502.62862818,  2238.286014  ])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "##SOL\n",
    "tmrca = sum(n9_times_array)\n",
    "np.mean(n9_times_array, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWEAAAFuCAYAAABQqakCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZAElEQVR4nO3dfbRldX3f8fcHZhQNWqCMOgwzopaYoI2gUxR0pUSNAWrEWBRYPqJZgNFEm2QlmrjU2pUuja2tihFIRMQSxQckaPApKoKhIAMdkAfRUVzOCAmDVBC11oFv/9j7lsPlztwzzN33d86979dae9398Nt7f8+Zy4d9f2fv30lVIUlqY7fWBUjScmYIS1JDhrAkNWQIS1JDhrAkNWQIS1JDUxnCSc5McmuSaxfoeHcn2dhPFyzEMSVpHJnG+4ST/DpwF3B2VT1xAY53V1XtueuVSdLOmcor4aq6GLh9dF2SxyX5XJIrk1yS5FcalSdJY5vKEN6OM4Dfr6qnAH8M/NVO7LtHkg1JLkvy/EGqk6Q5rGhdwEJIsidwOPDxJDOrH9xvewHwtjl2+0FV/VY/v66qbk7yWODLSb5RVd8Zum5JWhIhTHdF/6OqOnj2hqo6DzhvRztX1c39z+8muQg4BDCEJQ1uSXRHVNWdwE1JXgiQzpPG2TfJ3klmrpr3BZ4OXD9YsZI0YipDOMlHgP8JPD7JliSvAl4MvCrJ1cB1wDFjHu5XgQ39fl8B3l5VhrCkRTGVt6hJ0lIxlVfCkrRUTN0Hc0ceeWR97nOfa12GJGX+JvObuivh2267rXUJkrRgpi6EJWkpMYQlqSFDWJIaMoQlqSFDWJIaMoQlqSFDWJIaMoQlqSFDWJIaMoQlqSFDWJIaMoQlqSFDWJIaMoSnwJq160gy9rRm7brWJUsa02DjCSfZA7iY7luPVwCfqKq3zGoT4N3A0cBPgVdU1VVD1TStbt6ymeNOv3Ts9ueefPiA1UhaSEMO6v5z4JlVdVeSlcDXkny2qi4baXMUcGA/PRV4f/9TkpaFwbojqnNXv7iyn2Z/od0xwNl928uAvZKsHqomSZo0g/YJJ9k9yUbgVuCLVXX5rCZrgM0jy1v6dbOPc1KSDUk2bN26dbB6JWmxDRrCVXV3VR0M7A8cmuSJs5rM9R1N9/v656o6o6rWV9X6VatWDVCpJLWxKHdHVNWPgIuAI2dt2gKsHVneH7h5MWqSpEkwWAgnWZVkr37+IcCzgW/OanYB8LJ0ngbcUVW3DFWTJE2aIe+OWA18KMnudGH/sar6TJJTAKrqNOBCutvTNtHdonbigPVI0sQZLISr6hrgkDnWnzYyX8BrhqpBkiadT8xJUkOGsCQ1ZAhLUkOGsCQ1ZAgvRbutcNQ1aUoMeYuaWrlnm6OuSVPCK+EGdnZ8YElLl1fCDTg+sKQZXglLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsGC3FSQZe1qzdl3riqUlY0XrAjQB7tnGcadfOnbzc08+fMBipOXFK2FJasgQlqSGDGFJasgQlqSGBgvhJGuTfCXJDUmuS/K6OdockeSOJBv76c1D1SNJk2jIuyO2AX9UVVcleRhwZZIvVtX1s9pdUlXPHbAOSZpYg10JV9UtVXVVP/9j4AZgzVDnk6RptCh9wkkOAA4BLp9j82FJrk7y2SRPWIx6JGlSDP6wRpI9gU8Cr6+qO2dtvgp4dFXdleRo4HzgwDmOcRJwEsC6dT6tJWnpGPRKOMlKugA+p6rOm729qu6sqrv6+QuBlUn2naPdGVW1vqrWr1q1asiSJWlRDXl3RIAPADdU1bu20+ZRfTuSHNrX88OhapKkSTNkd8TTgZcC30iysV/3Z8A6gKo6DTgWeHWSbcDPgOOrqgasSZImymAhXFVfAzJPm1OBU4eqQZImnU/MSVJDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNWQIS1JDhrAkNTRYCCdZm+QrSW5Icl2S183RJknek2RTkmuSPHmoeiRpEq0Y8NjbgD+qqquSPAy4MskXq+r6kTZHAQf201OB9/c/JWlZGOxKuKpuqaqr+vkfAzcAa2Y1OwY4uzqXAXslWT1UTZI0aRalTzjJAcAhwOWzNq0BNo8sb+H+QU2Sk5JsSLJh69atg9UpSYtt8BBOsifwSeD1VXXn7M1z7FL3W1F1RlWtr6r1q1atGqJMSWpi0BBOspIugM+pqvPmaLIFWDuyvD9w85A1SdIkGfLuiAAfAG6oqndtp9kFwMv6uySeBtxRVbcMVZMkTZoh7454OvBS4BtJNvbr/gxYB1BVpwEXAkcDm4CfAicOWI8kTZzBQriqvsbcfb6jbQp4zVA1SNKk84m5BbBm7TqSjD1J0owhuyOWjZu3bOa40y8du/25Jx8+YDWSpolXwpLUkCEsSQ0ZwpLUkCEsSQ0ZwpLUkCEsSQ0ZwpLUkCEsSQ0ZwpLUkCEsSQ0ZwpLUkCEsSQ0ZwpLUkCEsSQ0ZwpLUkCEsSQ0ZwpLUkCEsSQ0ZwpLUkCEsSQ2NFcJJvjTOOknSztnhty0n2QN4KLBvkr2Bme9rfziw38C1SdKSN99X3p8MvJ4ucK/k3hC+E3jfcGVJ0vKwwxCuqncD707y+1X13kWqSZKWjfmuhAGoqvcmORw4YHSfqjp7oLokaVkYK4STfBh4HLARuLtfXYAhLEm7YKwQBtYDB1VVDVmMJC03494nfC3wqCELkaTlaNwr4X2B65N8Hfj5zMqqet4gVUnSMjFuCL91yCIkabka9+6Irw5diCQtR+PeHfFjurshAB4ErAR+UlUPH6owSVoOxr0SftjocpLnA4cOUZAkLScPaBS1qjofeObCliJJy8+43REvGFncje6+Ye8ZlqRdNO7dEb89Mr8N+B5wzIJXI0nLzLh9wicOXYgkLUfjDuq+f5JPJbk1yT8n+WSS/YcuTpKWunE/mPsgcAHduMJrgE/367YryZl9aF+7ne1HJLkjycZ+evPOFC5JS8G4Ibyqqj5YVdv66Sxg1Tz7nAUcOU+bS6rq4H5625i1SNKSMW4I35bkJUl276eXAD/c0Q5VdTFw+y5XKElL2Lgh/ErgRcA/AbcAxwIL8WHdYUmuTvLZJE/YXqMkJyXZkGTD1q1bF+C0kjQZxg3h/wS8vKpWVdUj6EL5rbt47quAR1fVk4D3Audvr2FVnVFV66tq/apV8/WCSNL0GDeEf62q/vfMQlXdDhyyKyeuqjur6q5+/kJgZZJ9d+WYkjRtxg3h3fqvvAcgyT6M/6DHnJI8Kkn6+UP7WnbYzyxJS824QfpfgUuTfILuceUXAX+xox2SfAQ4Atg3yRbgLXSjr1FVp9H1K786yTbgZ8Dxfn2SpOVm3Cfmzk6ygW7QngAvqKrr59nnhHm2nwqcOm6hkrQUjd2l0IfuDoNXy8RuK+h7ksay3/5r+cHm7w9YkDS9dqlfV8vUPds47vRLx25+7smHD1iMNN0e0HjCkqSFYQhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhreP2AP+NOa9aua12xtGgcwEfDc8Afabu8EpakhgxhSWrIEJakhgxhSWrIEJakhgxhSWrIEJakhgxhSWrIEJakhgxhSWrIEJakhgxhSWrIEJakhgxhSWrIEJakhgxhSWrIEJakhgxhSWrIEJakhgxhSWrIEJakhgxhSWrIEJakhgxhSWposBBOcmaSW5Ncu53tSfKeJJuSXJPkyUPVIkmTasgr4bOAI3ew/SjgwH46CXj/gLVI0kQaLISr6mLg9h00OQY4uzqXAXslWT1UPZI0iVr2Ca8BNo8sb+nXSdKy0TKEM8e6mrNhclKSDUk2bN26deCyJGnxtAzhLcDakeX9gZvnalhVZ1TV+qpav2rVqkUpTpIWQ8sQvgB4WX+XxNOAO6rqlob1SNKiWzHUgZN8BDgC2DfJFuAtwEqAqjoNuBA4GtgE/BQ4cahaJGlSDRbCVXXCPNsLeM1Q55ekaeATc5LUkCEsSQ0ZwpLUkCEsSQ0ZwnNYs3YdScaetMB2W7FT7/+atetaVyw9YIPdHTHNbt6ymeNOv3Ts9ueefPiA1SxD92zz/dey4ZWwJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4awJDVkCEtSQ4OGcJIjk9yYZFOSN8yx/YgkdyTZ2E9vHrIeSZo0K4Y6cJLdgfcBvwlsAa5IckFVXT+r6SVV9dyh6tAysNsKkozdfL/91/KDzd8fsCBpfIOFMHAosKmqvguQ5KPAMcDsEJZ2zT3bOO70S8dufu7Jhw9YjLRzhuyOWANsHlne0q+b7bAkVyf5bJInzHWgJCcl2ZBkw9atW4eoVctJf+U87rRm7brWFWsJG/JKeK6/D2vW8lXAo6vqriRHA+cDB95vp6ozgDMA1q9fP/sY0s7xylkTZMgr4S3A2pHl/YGbRxtU1Z1VdVc/fyGwMsm+A9YkSRNlyBC+AjgwyWOSPAg4HrhgtEGSR6X/RCXJoX09PxywJkmaKIN1R1TVtiSvBT4P7A6cWVXXJTml334acCzw6iTbgJ8Bx1eV3Q2Slo0h+4RnuhgunLXutJH5U4FTh6xBkiaZT8xJUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1tCxCeM3adTv1TQqStFgGHUVtUty8ZbPfpCBpIi2LK2FJmlSGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsCQ1ZAhLUkOGsLTAdvbhoDVr17UuWQ0ti4c1pMXkw0HaGV4JS1JDhrAkNWQIS/PZbcWwA0Dt5PHtQ15a7BOW5nPPtmH7eIc+viaaV8KS1JAhLEkNGcKS1JAhLEkNGcKS1JAhLEkNGcKS1JAhLEkNGcLSEueobpPNJ+akJc5R3SabISxNm36sCS0NhrA0bRxrYkmxT1iSGho0hJMcmeTGJJuSvGGO7Unynn77NUmePGQ9kjRpBgvhJLsD7wOOAg4CTkhy0KxmRwEH9tNJwPuHqkfSmKZ8fONpuxtkyD7hQ4FNVfVdgCQfBY4Brh9pcwxwdlUVcFmSvZKsrqpbBqxL0o5MeZ/ztN0Nki7/BjhwcixwZFX9br/8UuCpVfXakTafAd5eVV/rl78E/GlVbZh1rJPorpQBHg/cOEjRC29f4LbWRewC62/L+tuar/7bqurIXT3JkFfCc91DMzvxx2lDVZ0BnLEQRS2mJBuqan3rOh4o62/L+ttarPqH/GBuC7B2ZHl/4OYH0EaSlqwhQ/gK4MAkj0nyIOB44IJZbS4AXtbfJfE04A77gyUtJ4N1R1TVtiSvBT4P7A6cWVXXJTml334acCFwNLAJ+Clw4lD1NDJ1XSizWH9b1t/WotQ/2AdzkqT5+cScJDVkCEtSQ4bwPJKsTfKVJDckuS7J6/r1+yT5YpJv9z/3Htnnjf2j2Dcm+a2R9U9J8o1+23vSD4WV5MFJzu3XX57kgAV+Dbsn+V/9fdlTVXt/jr2SfCLJN/t/h8Om6TUk+Q/97861ST6SZI9Jrj/JmUluTXLtyLpFqTfJy/tzfDvJyxew/nf2vz/XJPlUkr0mpv6qctrBBKwGntzPPwz4Ft1j2H8JvKFf/wbgHf38QcDVwIOBxwDfAXbvt30dOIzu/ujPAkf1638POK2fPx44d4Ffwx8Cfwt8pl+emtr7434I+N1+/kHAXtPyGoA1wE3AQ/rljwGvmOT6gV8HngxcO7Ju8HqBfYDv9j/37uf3XqD6nwOs6OffMUn1Nw+5aZuAvwN+k+6pvdX9utXAjf38G4E3jrT/fP8PuRr45sj6E4DTR9v08yvontLJAtW7P/Al4JncG8JTUXt/zIfThVhmrZ+K10AXwpv7/zBXAJ/pA2Gi6wcO4L4hNni9o236bacDJyxE/bO2/Q5wzqTUb3fETuj/7DgEuBx4ZPX3NPc/H9E3m/mPbsaWft2afn72+vvsU1XbgDuAf7lAZf934E+Ae0bWTUvtAI8FtgIf7LtU/ibJL03La6iqHwD/Bfg+cAvdvfBfmJb6RyxGvds71kJ7Jd2V7X1qmXXORavfEB5Tkj2BTwKvr6o7d9R0jnW1g/U72meXJHkucGtVXTnuLtupY9FrH7GC7k/L91fVIcBP6P4c3p6Jeg193+kxdH/q7gf8UpKX7GiX7dTS8t9gRxay3sFfR5I/B7YB5+xCLQtavyE8hiQr6QL4nKo6r1/9z0lW99tXA7f267f3KPaWfn72+vvsk2QF8C+A2xeg9KcDz0vyPeCjwDOT/I8pqX3GFmBLVV3eL3+CLpSn5TU8G7ipqrZW1S+A84DDp6j+GYtR76DDGPQflD0XeHH1/QWTUL8hPI/+E9EPADdU1btGNl0AzHz6+XK6vuKZ9cf3n6A+hm6s5K/3f8L9OMnT+mO+bNY+M8c6FvjyyC/JA1ZVb6yq/avqALoPEL5cVS+ZhtpHXsM/AZuTPL5f9Sy64VCn5TV8H3hakof2530WcMMU1T9jMer9PPCcJHv3f0E8p1+3y5IcCfwp8Lyq+ums19W2/l3pvF8OE/AMuj8prgE29tPRdH1AXwK+3f/cZ2SfP6f7lPVG+k9U+/XrgWv7bady7xOLewAfp3t8++vAYwd4HUdw7wdz01b7wcCG/t/gfLpPnqfmNQD/Efhmf+4P030SP7H1Ax+h67/+Bd3V3asWq166/tpN/XTiAta/ia6/dmM/nTYp9fvYsiQ1ZHeEJDVkCEtSQ4awJDVkCEtSQ4awJDVkCGtO6UYu+72R5f2SfGIRz7869476tj7Jexbr3AstySuS7Dey/DdJDlrA4//rJGct1PG0uLxFTXPqx8n4TFU9sdH53wl8rar+bt7GEyDJ7lV193a2XQT8cVVtGPD8/wC8sqq+P9Q5NAyvhLU9bwcel2RjPxbrATPjs/ZXducn+XSSm5K8Nskf9gPsXJZkn77d45J8LsmVSS5J8iv9+hemG1v36iQXb+f8/x74XN/+iJGr4remGy/2oiTfTfIHc+2c5FVJvtW3++skp/brVyX5ZJIr+unp8x03yUuSfL1/L05Psnu//q4kb0tyOXBYkjf3x7w2yRnpHEt30/85/f4P6c+xvj/GCenGrL02yTtGznlXkr/o36PLkjxynvfu03RPRWraLPTTTU5LY+L+Qxn+/2W68XA30Y2vvIpuFKlT+m3/jW6QI+ierDqwn38q3eOdAN8A1vTze81x7scAV44sH8G9T/u9FbiU7qmzfYEfAitn7b8f8D264SNXApcAp/bb/hZ4Rj+/ju5x9O0eF/hVuoBb2bf7K+Bl/XwBLxo57+hTZB8GfrufvwhYP7LtIrpg3o/useZVdAMVfRl4/sixZ/b/S+BNO3rv6MYJ+XTr3xunnZ8G+7ZlLXlfqaof0z1ffwddUEEXEr+WbtS5w4GPd4/eA13AAfwjcFaSj9ENaDPbarrhK7fn76vq58DPk9wKPJL7Djt4KPDVqrodIMnHgV/utz0bOGikpocnedgOjvss4CnAFf0+D+HewWvuphvYacZvJPkT4KF0/wO4buR9mcu/AS6qqq19nefQDUh+PvB/6cYeBriSbgxr2P57dytdqGvKGMJ6oH4+Mn/PyPI9dL9XuwE/qqqDZ+9YVackeSrw74CNSQ6uqh+ONPkZ3fP545z7bu7/ezzXkIIzdqMbkPtn99mhC9i5jhvgQ1X1xjmO9X+q7wdOsgfdVfL6qtqc5K3zvIb56vxF9Ze4I7Xs6L3bg+5905SxT1jb82O67oYHpLoxl29K8kLoRqNL8qR+/nFVdXlVvZnuWwnWztr9W3TdHw/U14F/m240qxV0/cszvgC8dmYhycHzHOtLwLFJHtG33yfJo+doNxO4t/V/BRw7sm177+XlfZ379v3MJwBf3VExO3jvfplusBlNGUNYc+qvrv6x/xDonQ/wMC8GXpXkaro/zY/p179z5sMo4GK67/gaPfdPgO8k+VcPsPYfAP+ZLuT+gW7oyzv6zX8ArE/3hY/XA6fMc6zrgTcBX0hyDfBFuu6S2e1+BPw1XXfM+cAVI5vPAk6b+WBuZJ9b6L5e5yt078FVNf/dINt7734D+Pt59tUE8hY1TaQkvwM8pare9AD337Oq7uqvhD8FnFlVn1rQIidEkgfTXUE/o7qv29EU8UpYE6kPzO/twiHemmQj3Z/oN9FdnS5V6+i+CdkAnkJeCUtSQ14JS1JDhrAkNWQIS1JDhrAkNWQIS1JD/w/I0ot36IWu8gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "##SOL\n",
    "p = sns.displot(x = tmrca, stat = \"density\", label = 'span of trees', binwidth = 5_000)\n",
    "plt.xlabel(\"times (in generations)\");\n",
    "plt.ylabel(\"count\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 6:</b> What is the mean $t_{mrca}$ arcoss replicates?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean tmrca is 36597.53044624645\n"
     ]
    }
   ],
   "source": [
    "#SOL\n",
    "print(\"The mean tmrca is\", np.mean(tmrca))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 7:</b> How well does this match the theoretical expectation $E[t_{mrca}]$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>PopGen hint:</b>\n",
    "work out the expected times between coalescent intervals first and use the fact that $E[t_{mrca}] = \\sum_{i=2}^{n-1} E[t_i]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35555.555555555555"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#SOL\n",
    "expected_intervals = [2/int(i * (i-1)) for i in range(2, 10)]\n",
    "sum(expected_intervals) * 20_000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 8:</b> What is missing from this simulation?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL\n",
    "#Topologies! A more realistic model of the population."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 A simple `msprime` ancestry simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can set up a single simulation using `sim_ancestry()` which is a method of the coalescent simulator `msprime`. The `samples` argument of `sim_ancestry` defines the sampled individuals whose history will be simulated. In the simplest case, we provide a single integer which defines the number of samples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> Typing the variable assigned to this <code>msprime</code> simulation generates a useful, tabular summary of the simulation output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "            <div>\n",
       "              <style>\n",
       "                .tskit-table thead tr th {text-align: left;padding: 0.5em 0.5em;}\n",
       "                .tskit-table tbody tr td {padding: 0.5em 0.5em;}\n",
       "                .tskit-table tbody tr td:first-of-type {text-align: left;}\n",
       "                .tskit-details-label {vertical-align: top; padding-right:5px;}\n",
       "                .tskit-table-set {display: inline-flex;flex-wrap: wrap;margin: -12px 0 0 -12px;width: calc(100% + 12px);}\n",
       "                .tskit-table-set-table {margin: 12px 0 0 12px;}\n",
       "                details {display: inline-block;}\n",
       "                summary {cursor: pointer; outline: 0; display: list-item;}\n",
       "              </style>\n",
       "              <div class=\"tskit-table-set\">\n",
       "                <div class=\"tskit-table-set-table\">\n",
       "                  <table class=\"tskit-table\">\n",
       "                    <thead>\n",
       "                      <tr>\n",
       "                        <th style=\"padding:0;line-height:21px;\">\n",
       "                          <img style=\"height: 32px;display: inline-block;padding: 3px 5px 3px 0;\" src=\"https://raw.githubusercontent.com/tskit-dev/administrative/main/tskit_logo.svg\"/>\n",
       "                          <a target=\"_blank\" href=\"https://tskit.readthedocs.io/en/latest/python-api.html#the-treesequence-class\"> Tree Sequence </a>\n",
       "                        </th>\n",
       "                      </tr>\n",
       "                    </thead>\n",
       "                    <tbody>\n",
       "                      <tr><td>Trees</td><td>1</td></tr>\n",
       "                      <tr><td>Sequence Length</td><td>1.0</td></tr>\n",
       "                      <tr><td>Sample Nodes</td><td>18</td></tr>\n",
       "                      <tr><td>Total Size</td><td>3.4 KiB</td></tr>\n",
       "                      <tr>\n",
       "                        <td>Metadata</td><td style=\"text-align: left;\">No Metadata</td></tr>\n",
       "                    </tbody>\n",
       "                  </table>\n",
       "                </div>\n",
       "                <div class=\"tskit-table-set-table\">\n",
       "                  <table class=\"tskit-table\">\n",
       "                    <thead>\n",
       "                      <tr>\n",
       "                        <th style=\"line-height:21px;\">Table</th>\n",
       "                        <th>Rows</th>\n",
       "                        <th>Size</th>\n",
       "                        <th>Has Metadata</th>\n",
       "                      </tr>\n",
       "                    </thead>\n",
       "                    <tbody>\n",
       "                    \n",
       "                  <tr>\n",
       "                    <td>Edges</td>\n",
       "                      <td>34</td>\n",
       "                      <td>956 Bytes</td>\n",
       "                      <td style=\"text-align: center;\">\n",
       "                        \n",
       "                      </td>\n",
       "                    </tr>\n",
       "                \n",
       "                  <tr>\n",
       "                    <td>Individuals</td>\n",
       "                      <td>9</td>\n",
       "                      <td>156 Bytes</td>\n",
       "                      <td style=\"text-align: center;\">\n",
       "                        \n",
       "                      </td>\n",
       "                    </tr>\n",
       "                \n",
       "                  <tr>\n",
       "                    <td>Migrations</td>\n",
       "                      <td>0</td>\n",
       "                      <td>4 Bytes</td>\n",
       "                      <td style=\"text-align: center;\">\n",
       "                        \n",
       "                      </td>\n",
       "                    </tr>\n",
       "                \n",
       "                  <tr>\n",
       "                    <td>Mutations</td>\n",
       "                      <td>0</td>\n",
       "                      <td>8 Bytes</td>\n",
       "                      <td style=\"text-align: center;\">\n",
       "                        \n",
       "                      </td>\n",
       "                    </tr>\n",
       "                \n",
       "                  <tr>\n",
       "                    <td>Nodes</td>\n",
       "                      <td>35</td>\n",
       "                      <td>844 Bytes</td>\n",
       "                      <td style=\"text-align: center;\">\n",
       "                        \n",
       "                      </td>\n",
       "                    </tr>\n",
       "                \n",
       "                  <tr>\n",
       "                    <td>Populations</td>\n",
       "                      <td>1</td>\n",
       "                      <td>216 Bytes</td>\n",
       "                      <td style=\"text-align: center;\">\n",
       "                        ✅\n",
       "                      </td>\n",
       "                    </tr>\n",
       "                \n",
       "                  <tr>\n",
       "                    <td>Provenances</td>\n",
       "                      <td>1</td>\n",
       "                      <td>966 Bytes</td>\n",
       "                      <td style=\"text-align: center;\">\n",
       "                        \n",
       "                      </td>\n",
       "                    </tr>\n",
       "                \n",
       "                  <tr>\n",
       "                    <td>Sites</td>\n",
       "                      <td>0</td>\n",
       "                      <td>8 Bytes</td>\n",
       "                      <td style=\"text-align: center;\">\n",
       "                        \n",
       "                      </td>\n",
       "                    </tr>\n",
       "                \n",
       "                    </tbody>\n",
       "                  </table>\n",
       "                </div>\n",
       "              </div>\n",
       "            </div>\n",
       "            "
      ],
      "text/plain": [
       "<tskit.trees.TreeSequence at 0x7fc95c808a00>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simple_ts = msprime.sim_ancestry(9);\n",
    "simple_ts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> Strictly speaking, <code>sim_ancestry()</code> always returns a treesequence (more on these practical 3). However, in this case we have only simulated a single tree. We can use the method <code>.first()</code> to extract this tree:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_tree = simple_ts.first()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> We can plot our simulated tree using the method <code>.draw_svg()</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" baseProfile=\"full\" height=\"500\" version=\"1.1\" width=\"500\"><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 n34 p0 root\" transform=\"translate(280.347 26.8)\"><g class=\"a34 c2 node n31 p0\" transform=\"translate(118.194 420.406)\"><g class=\"a31 c2 node n20 p0\" transform=\"translate(-46.3194 15.9837)\"><g class=\"a20 i1 leaf node n2 p0 sample\" transform=\"translate(-12.7778 5.01018)\"><path class=\"edge\" d=\"M 0 0 V -5.01018 H 12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">2</text></g><g class=\"a20 i1 leaf node n3 p0 sample\" transform=\"translate(12.7778 5.01018)\"><path class=\"edge\" d=\"M 0 0 V -5.01018 H -12.7778\"/><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 -15.9837 H 46.3194\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">20</text></g><g class=\"a31 c2 node n27 p0\" transform=\"translate(46.3194 9.85502)\"><g class=\"a27 i8 leaf node n17 p0 sample\" transform=\"translate(22.3611 11.1388)\"><path class=\"edge\" d=\"M 0 0 V -11.1388 H -22.3611\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">17</text></g><g class=\"a27 c2 node n25 p0\" transform=\"translate(-22.3611 2.14952)\"><g class=\"a25 i6 leaf node n12 p0 sample\" transform=\"translate(19.1667 8.98932)\"><path class=\"edge\" d=\"M 0 0 V -8.98932 H -19.1667\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">12</text></g><g class=\"a25 c2 node n23 p0\" transform=\"translate(-19.1667 2.09815)\"><g class=\"a23 i2 leaf node n4 p0 sample\" transform=\"translate(-12.7778 6.89118)\"><path class=\"edge\" d=\"M 0 0 V -6.89118 H 12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">4</text></g><g class=\"a23 i3 leaf node n7 p0 sample\" transform=\"translate(12.7778 6.89118)\"><path class=\"edge\" d=\"M 0 0 V -6.89118 H -12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">7</text></g><path class=\"edge\" d=\"M 0 0 V -2.09815 H 19.1667\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">23</text></g><path class=\"edge\" d=\"M 0 0 V -2.14952 H 22.3611\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">25</text></g><path class=\"edge\" d=\"M 0 0 V -9.85502 H -46.3194\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">27</text></g><path class=\"edge\" d=\"M 0 0 V -420.406 H -118.194\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">31</text></g><g class=\"a34 c2 node n33 p0\" transform=\"translate(-118.194 371.898)\"><g class=\"a33 c2 node n30 p0\" transform=\"translate(-83.0556 52.3049)\"><g class=\"a30 c2 node n19 p0\" transform=\"translate(-33.5417 13.3297)\"><g class=\"a19 i0 leaf node n0 p0 sample\" transform=\"translate(-12.7778 3.86737)\"><path class=\"edge\" d=\"M 0 0 V -3.86737 H 12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">0</text></g><g class=\"a19 i4 leaf node n9 p0 sample\" transform=\"translate(12.7778 3.86737)\"><path class=\"edge\" d=\"M 0 0 V -3.86737 H -12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">9</text></g><path class=\"edge\" d=\"M 0 0 V -13.3297 H 33.5417\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">19</text></g><g class=\"a30 c2 node n26 p0\" transform=\"translate(33.5417 7.44939)\"><g class=\"a26 i0 leaf node n1 p0 sample\" transform=\"translate(-28.75 9.74772)\"><path class=\"edge\" d=\"M 0 0 V -9.74772 H 28.75\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">1</text></g><g class=\"a26 c2 node n24 p0\" transform=\"translate(28.75 2.58928)\"><g class=\"a24 i6 leaf node n13 p0 sample\" transform=\"translate(19.1667 7.15845)\"><path class=\"edge\" d=\"M 0 0 V -7.15845 H -19.1667\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">13</text></g><g class=\"a24 c2 node n21 p0\" transform=\"translate(-19.1667 1.5807)\"><g class=\"a21 i2 leaf node n5 p0 sample\" transform=\"translate(-12.7778 5.57775)\"><path class=\"edge\" d=\"M 0 0 V -5.57775 H 12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">5</text></g><g class=\"a21 i5 leaf node n11 p0 sample\" transform=\"translate(12.7778 5.57775)\"><path class=\"edge\" d=\"M 0 0 V -5.57775 H -12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">11</text></g><path class=\"edge\" d=\"M 0 0 V -1.5807 H 19.1667\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">21</text></g><path class=\"edge\" d=\"M 0 0 V -2.58928 H -28.75\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">24</text></g><path class=\"edge\" d=\"M 0 0 V -7.44939 H -33.5417\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">26</text></g><path class=\"edge\" d=\"M 0 0 V -52.3049 H 83.0556\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">30</text></g><g class=\"a33 c2 node n32 p0\" transform=\"translate(83.0556 38.6265)\"><g class=\"a32 c2 node n22 p0\" transform=\"translate(-46.3194 24.3525)\"><g class=\"a22 i3 leaf node n6 p0 sample\" transform=\"translate(-12.7778 6.52302)\"><path class=\"edge\" d=\"M 0 0 V -6.52302 H 12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">6</text></g><g class=\"a22 i8 leaf node n16 p0 sample\" transform=\"translate(12.7778 6.52302)\"><path class=\"edge\" d=\"M 0 0 V -6.52302 H -12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">16</text></g><path class=\"edge\" d=\"M 0 0 V -24.3525 H 46.3194\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">22</text></g><g class=\"a32 c2 node n29 p0\" transform=\"translate(46.3194 14.7382)\"><g class=\"a29 i5 leaf node n10 p0 sample\" transform=\"translate(22.3611 16.1373)\"><path class=\"edge\" d=\"M 0 0 V -16.1373 H -22.3611\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">10</text></g><g class=\"a29 c2 node n28 p0\" transform=\"translate(-22.3611 3.25)\"><g class=\"a28 i7 leaf node n15 p0 sample\" transform=\"translate(19.1667 12.8873)\"><path class=\"edge\" d=\"M 0 0 V -12.8873 H -19.1667\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">15</text></g><g class=\"a28 c2 node n18 p0\" transform=\"translate(-19.1667 10.5739)\"><g class=\"a18 i4 leaf node n8 p0 sample\" transform=\"translate(-12.7778 2.31347)\"><path class=\"edge\" d=\"M 0 0 V -2.31347 H 12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">8</text></g><g class=\"a18 i7 leaf node n14 p0 sample\" transform=\"translate(12.7778 2.31347)\"><path class=\"edge\" d=\"M 0 0 V -2.31347 H -12.7778\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">14</text></g><path class=\"edge\" d=\"M 0 0 V -10.5739 H 19.1667\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">18</text></g><path class=\"edge\" d=\"M 0 0 V -3.25 H 22.3611\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">28</text></g><path class=\"edge\" d=\"M 0 0 V -14.7382 H -46.3194\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">29</text></g><path class=\"edge\" d=\"M 0 0 V -38.6265 H -83.0556\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">32</text></g><path class=\"edge\" d=\"M 0 0 V -371.898 H 118.194\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">33</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab\" transform=\"translate(0 -11)\">34</text></g></g></g></svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SVG(simple_tree.draw_svg(size=(500, 500)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 9:</b>  Why does the tree we simulated contain 18 sample nodes and 35 nodes in total? How would you change the arguments of <code>sim_ancestry</code> to match the simulation we have assumed for our DIY coalescence time simulator above?\n",
    "\n",
    "<b>Hint: </b> Check what arguments <a href=\"https://tskit.dev/msprime/docs/stable/ancestry.html\"><code>sim_ancestry</code></a> takes. \n",
    "        </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" baseProfile=\"full\" height=\"500\" version=\"1.1\" width=\"500\"><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 n16 p0 root\" transform=\"translate(218.056 26.8)\"><g class=\"a16 c2 node n14 p0\" transform=\"translate(-108.611 317.195)\"><g class=\"a14 i5 leaf node n5 p0 sample\" transform=\"translate(38.3333 124.205)\"><path class=\"edge\" d=\"M 0 0 V -124.205 H -38.3333\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">5</text></g><g class=\"a14 c2 node n9 p0\" transform=\"translate(-38.3333 117.776)\"><g class=\"a9 i0 leaf node n0 p0 sample\" transform=\"translate(-25.5556 6.42851)\"><path class=\"edge\" d=\"M 0 0 V -6.42851 H 25.5556\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">0</text></g><g class=\"a9 i1 leaf node n1 p0 sample\" transform=\"translate(25.5556 6.42851)\"><path class=\"edge\" d=\"M 0 0 V -6.42851 H -25.5556\"/><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 -117.776 H 38.3333\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">9</text></g><path class=\"edge\" d=\"M 0 0 V -317.195 H 108.611\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">14</text></g><g class=\"a16 c2 node n15 p0\" transform=\"translate(108.611 182.002)\"><g class=\"a15 c2 node n12 p0\" transform=\"translate(-89.4444 216.609)\"><g class=\"a12 i2 leaf node n2 p0 sample\" transform=\"translate(-38.3333 42.7885)\"><path class=\"edge\" d=\"M 0 0 V -42.7885 H 38.3333\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">2</text></g><g class=\"a12 c2 node n11 p0\" transform=\"translate(38.3333 3.06194)\"><g class=\"a11 i3 leaf node n3 p0 sample\" transform=\"translate(-25.5556 39.7265)\"><path class=\"edge\" d=\"M 0 0 V -39.7265 H 25.5556\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">3</text></g><g class=\"a11 i8 leaf node n8 p0 sample\" transform=\"translate(25.5556 39.7265)\"><path class=\"edge\" d=\"M 0 0 V -39.7265 H -25.5556\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">8</text></g><path class=\"edge\" d=\"M 0 0 V -3.06194 H -38.3333\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">11</text></g><path class=\"edge\" d=\"M 0 0 V -216.609 H 89.4444\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">12</text></g><g class=\"a15 c2 node n13 p0\" transform=\"translate(89.4444 188.634)\"><g class=\"a13 i7 leaf node n7 p0 sample\" transform=\"translate(38.3333 70.7644)\"><path class=\"edge\" d=\"M 0 0 V -70.7644 H -38.3333\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">7</text></g><g class=\"a13 c2 node n10 p0\" transform=\"translate(-38.3333 61.7752)\"><g class=\"a10 i4 leaf node n4 p0 sample\" transform=\"translate(-25.5556 8.9892)\"><path class=\"edge\" d=\"M 0 0 V -8.9892 H 25.5556\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">4</text></g><g class=\"a10 i6 leaf node n6 p0 sample\" transform=\"translate(25.5556 8.9892)\"><path class=\"edge\" d=\"M 0 0 V -8.9892 H -25.5556\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">6</text></g><path class=\"edge\" d=\"M 0 0 V -61.7752 H 38.3333\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">10</text></g><path class=\"edge\" d=\"M 0 0 V -188.634 H -89.4444\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">13</text></g><path class=\"edge\" d=\"M 0 0 V -182.002 H -108.611\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">15</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab\" transform=\"translate(0 -11)\">16</text></g></g></g></svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#SOL: define ploidy = 1. The point is to realize that msprime and pop gen models in general assume diploidy.\n",
    "simple_ts = msprime.sim_ancestry(9, ploidy=1);\n",
    "simple_tree = simple_ts.first()\n",
    "SVG(simple_tree.draw_svg(size=(500, 500)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 10:</b> Use the <a href=\"https://tskit.dev/tskit/docs/stable/python-api.html#the-tree-class\"><code>.root</code></a> attribute to find the rootnode of this tree and the <code>time()</code> method to find the $t_{mrca}$.\n",
    "        </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The tmrca is 2.857082985053013\n"
     ]
    }
   ],
   "source": [
    "#SOL\n",
    "root = simple_tree.root\n",
    "#Having defined the root, we can use the `time` method to find the time of the root node, aka the $T_{mrca}$:\n",
    "print(\"The tmrca is\", simple_tree.time(root))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use a `for` loop to loop over all nodes in the tree to find the associated coalescence times. Note the use of the `.nodes()` method in the loop below which returns a `generator` object. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "generator"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(simple_tree.nodes())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.041610295455245354, 0.058185051298223545, 0.25714077270780517, 0.2769600325367253, 0.4580420862953657, 0.8039479305969621, 1.6790245389157845, 2.857082985053013]\n"
     ]
    }
   ],
   "source": [
    "coal_times =[simple_tree.time(node) for node in simple_tree.nodes()]\n",
    "print(sorted(coal_times))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.4  A simple `msprime` mutation simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have so far only simulated genealogies, i.e. genetic ancestry without worrying about mutations. \n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>PopGen:</b> One of the most useful properties of the coalescent is that it allows us to separate the stochatic process of genetic ancestry (i.e. drift backwards) from the stochatic process of mutations. We can sprinkle mutations onto a simulated tree in a second step.\n",
    "        </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> We use the method <code>sim_mutations</code> to sprinkle mutations onto our simlulated ancestry:\n",
    "        </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_with_muts_1 = msprime.sim_mutations(simple_ts, rate=0.2)\n",
    "ts_with_muts_2 = msprime.sim_mutations(simple_ts, rate=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 11:</b> Visualize the trees for both replicates and check how they differ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" baseProfile=\"full\" height=\"500\" version=\"1.1\" width=\"500\"><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=\"M56.8,0 l423.2,0 l0,438.2 l0,25 l0,5 l-423.2,0 l0,-5 l0,-25 l0,-438.2z\"/></g><g class=\"axes\"><g class=\"x-axis\"><g transform=\"translate(268.4,500)\"><text class=\"lab\" text-anchor=\"middle\" transform=\"translate(0 -11)\">Genome position</text></g><line x1=\"56.8\" x2=\"480\" y1=\"463.2\" y2=\"463.2\"/><g class=\"ticks\"><g class=\"tick\" transform=\"translate(56.8 463.2)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"5\"/><g transform=\"translate(0,6)\"><text class=\"lab\">0</text></g></g><g class=\"tick\" transform=\"translate(480 463.2)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"5\"/><g transform=\"translate(0,6)\"><text class=\"lab\">1</text></g></g></g><g class=\"site s0\" transform=\"translate(56.8 463.2)\"><line class=\"sym\" x1=\"0\" x2=\"0\" y1=\"0\" y2=\"-10\"/><g class=\"mut m1\"><polyline class=\"sym\" points=\"2.5,-6.5 0,-1.5 -2.5,-6.5\"/></g><g class=\"mut m0\"><polyline class=\"sym\" points=\"2.5,-10.5 0,-5.5 -2.5,-10.5\"/></g></g></g><g class=\"y-axis\"><g transform=\"translate(0,215.7)\"><text class=\"lab\" text-anchor=\"middle\" transform=\"translate(11) rotate(-90)\">Time</text></g><line x1=\"56.8\" x2=\"56.8\" y1=\"421.4\" y2=\"10\"/><g class=\"ticks\"><g class=\"tick\" transform=\"translate(56.8 421.4)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.00</text></g></g><g class=\"tick\" transform=\"translate(56.8 415.653)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.04</text></g></g><g class=\"tick\" transform=\"translate(56.8 413.364)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.06</text></g></g><g class=\"tick\" transform=\"translate(56.8 385.886)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.26</text></g></g><g class=\"tick\" transform=\"translate(56.8 383.148)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.28</text></g></g><g class=\"tick\" transform=\"translate(56.8 358.138)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.46</text></g></g><g class=\"tick\" transform=\"translate(56.8 310.364)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.80</text></g></g><g class=\"tick\" transform=\"translate(56.8 189.505)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">1.68</text></g></g><g class=\"tick\" transform=\"translate(56.8 26.8)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">2.86</text></g></g></g></g></g><g class=\"plotbox trees\"><g class=\"tree t0\" transform=\"translate(56.8 0)\"><g class=\"plotbox\"><g class=\"c2 node n16 p0 root\" transform=\"translate(184.989 26.8)\"><g class=\"a16 c2 node n14 p0\" transform=\"translate(-90.4778 283.564)\"><g class=\"a14 i5 leaf node n5 p0 sample\" transform=\"translate(31.9333 111.036)\"><path class=\"edge\" d=\"M 0 0 V -111.036 H -31.9333\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">5</text></g><g class=\"a14 c2 node n9 p0\" transform=\"translate(-31.9333 105.289)\"><g class=\"a9 i0 leaf node n0 p0 sample\" transform=\"translate(-21.2889 5.74692)\"><path class=\"edge\" d=\"M 0 0 V -5.74692 H 21.2889\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">0</text></g><g class=\"a9 i1 leaf node n1 p0 sample\" transform=\"translate(21.2889 5.74692)\"><path class=\"edge\" d=\"M 0 0 V -5.74692 H -21.2889\"/><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 -105.289 H 31.9333\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">9</text></g><path class=\"edge\" d=\"M 0 0 V -283.564 H 90.4778\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">14</text></g><g class=\"a16 c2 node n15 p0\" transform=\"translate(90.4778 162.705)\"><g class=\"a15 c2 m0 node n12 p0 s0\" transform=\"translate(-74.5111 193.643)\"><g class=\"a12 i2 leaf node n2 p0 sample\" transform=\"translate(-31.9333 38.2518)\"><path class=\"edge\" d=\"M 0 0 V -38.2518 H 31.9333\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">2</text></g><g class=\"a12 c2 node n11 p0\" transform=\"translate(31.9333 2.7373)\"><g class=\"a11 i3 leaf node n3 p0 sample\" transform=\"translate(-21.2889 35.5145)\"><path class=\"edge\" d=\"M 0 0 V -35.5145 H 21.2889\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">3</text></g><g class=\"a11 i8 leaf node n8 p0 sample\" transform=\"translate(21.2889 35.5145)\"><path class=\"edge\" d=\"M 0 0 V -35.5145 H -21.2889\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">8</text></g><path class=\"edge\" d=\"M 0 0 V -2.7373 H -31.9333\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">11</text></g><path class=\"edge\" d=\"M 0 0 V -193.643 H 74.5111\"/><g class=\"mut m0 s0\" transform=\"translate(0 -157.652)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"157.652\"/><path class=\"sym\" d=\"M -3,-3 l 6,6 M -3,3 l 6,-6\"/><text class=\"lab lft\" transform=\"translate(-5 0)\">0</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">12</text></g><g class=\"a15 c2 m1 node n13 p0 s0\" transform=\"translate(74.5111 168.633)\"><g class=\"a13 i7 leaf node n7 p0 sample\" transform=\"translate(31.9333 63.2615)\"><path class=\"edge\" d=\"M 0 0 V -63.2615 H -31.9333\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">7</text></g><g class=\"a13 c2 node n10 p0\" transform=\"translate(-31.9333 55.2254)\"><g class=\"a10 i4 leaf node n4 p0 sample\" transform=\"translate(-21.2889 8.03611)\"><path class=\"edge\" d=\"M 0 0 V -8.03611 H 21.2889\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">4</text></g><g class=\"a10 i6 leaf node n6 p0 sample\" transform=\"translate(21.2889 8.03611)\"><path class=\"edge\" d=\"M 0 0 V -8.03611 H -21.2889\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">6</text></g><path class=\"edge\" d=\"M 0 0 V -55.2254 H 31.9333\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">10</text></g><path class=\"edge\" d=\"M 0 0 V -168.633 H -74.5111\"/><g class=\"mut m1 s0\" transform=\"translate(0 -129.702)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"129.702\"/><path class=\"sym\" d=\"M -3,-3 l 6,6 M -3,3 l 6,-6\"/><text class=\"lab rgt\" transform=\"translate(5 0)\">1</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">13</text></g><path class=\"edge\" d=\"M 0 0 V -162.705 H -90.4778\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">15</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab\" transform=\"translate(0 -11)\">16</text></g></g></g></g></g></svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" baseProfile=\"full\" height=\"500\" version=\"1.1\" width=\"500\"><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=\"M56.8,0 l423.2,0 l0,438.2 l0,25 l0,5 l-423.2,0 l0,-5 l0,-25 l0,-438.2z\"/></g><g class=\"axes\"><g class=\"x-axis\"><g transform=\"translate(268.4,500)\"><text class=\"lab\" text-anchor=\"middle\" transform=\"translate(0 -11)\">Genome position</text></g><line x1=\"56.8\" x2=\"480\" y1=\"463.2\" y2=\"463.2\"/><g class=\"ticks\"><g class=\"tick\" transform=\"translate(56.8 463.2)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"5\"/><g transform=\"translate(0,6)\"><text class=\"lab\">0</text></g></g><g class=\"tick\" transform=\"translate(480 463.2)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"5\"/><g transform=\"translate(0,6)\"><text class=\"lab\">1</text></g></g></g><g class=\"site s0\" transform=\"translate(56.8 463.2)\"><line class=\"sym\" x1=\"0\" x2=\"0\" y1=\"0\" y2=\"-10\"/><g class=\"mut m3\"><polyline class=\"sym\" points=\"2.5,-6.5 0,-1.5 -2.5,-6.5\"/></g><g class=\"mut m2\"><polyline class=\"sym\" points=\"2.5,-10.5 0,-5.5 -2.5,-10.5\"/></g><g class=\"mut m1\"><polyline class=\"sym\" points=\"2.5,-14.5 0,-9.5 -2.5,-14.5\"/></g><g class=\"mut m0\"><polyline class=\"sym\" points=\"2.5,-18.5 0,-13.5 -2.5,-18.5\"/></g></g></g><g class=\"y-axis\"><g transform=\"translate(0,215.7)\"><text class=\"lab\" text-anchor=\"middle\" transform=\"translate(11) rotate(-90)\">Time</text></g><line x1=\"56.8\" x2=\"56.8\" y1=\"421.4\" y2=\"10\"/><g class=\"ticks\"><g class=\"tick\" transform=\"translate(56.8 421.4)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.00</text></g></g><g class=\"tick\" transform=\"translate(56.8 415.653)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.04</text></g></g><g class=\"tick\" transform=\"translate(56.8 413.364)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.06</text></g></g><g class=\"tick\" transform=\"translate(56.8 385.886)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.26</text></g></g><g class=\"tick\" transform=\"translate(56.8 383.148)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.28</text></g></g><g class=\"tick\" transform=\"translate(56.8 358.138)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.46</text></g></g><g class=\"tick\" transform=\"translate(56.8 310.364)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">0.80</text></g></g><g class=\"tick\" transform=\"translate(56.8 189.505)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">1.68</text></g></g><g class=\"tick\" transform=\"translate(56.8 26.8)\"><line x1=\"0\" x2=\"-5\" y1=\"0\" y2=\"0\"/><g transform=\"translate(-6,0)\"><text class=\"lab\" text-anchor=\"end\">2.86</text></g></g></g></g></g><g class=\"plotbox trees\"><g class=\"tree t0\" transform=\"translate(56.8 0)\"><g class=\"plotbox\"><g class=\"c2 node n16 p0 root\" transform=\"translate(184.989 26.8)\"><g class=\"a16 c2 node n14 p0\" transform=\"translate(-90.4778 283.564)\"><g class=\"a14 i5 leaf node n5 p0 sample\" transform=\"translate(31.9333 111.036)\"><path class=\"edge\" d=\"M 0 0 V -111.036 H -31.9333\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">5</text></g><g class=\"a14 c2 node n9 p0\" transform=\"translate(-31.9333 105.289)\"><g class=\"a9 i0 leaf node n0 p0 sample\" transform=\"translate(-21.2889 5.74692)\"><path class=\"edge\" d=\"M 0 0 V -5.74692 H 21.2889\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">0</text></g><g class=\"a9 i1 leaf node n1 p0 sample\" transform=\"translate(21.2889 5.74692)\"><path class=\"edge\" d=\"M 0 0 V -5.74692 H -21.2889\"/><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 -105.289 H 31.9333\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">9</text></g><path class=\"edge\" d=\"M 0 0 V -283.564 H 90.4778\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">14</text></g><g class=\"a16 c2 m0 node n15 p0 s0\" transform=\"translate(90.4778 162.705)\"><g class=\"a15 c2 m2 node n12 p0 s0\" transform=\"translate(-74.5111 193.643)\"><g class=\"a12 i2 leaf node n2 p0 sample\" transform=\"translate(-31.9333 38.2518)\"><path class=\"edge\" d=\"M 0 0 V -38.2518 H 31.9333\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">2</text></g><g class=\"a12 c2 node n11 p0\" transform=\"translate(31.9333 2.7373)\"><g class=\"a11 i3 leaf m3 node n3 p0 s0 sample\" transform=\"translate(-21.2889 35.5145)\"><path class=\"edge\" d=\"M 0 0 V -35.5145 H 21.2889\"/><g class=\"mut m3 s0\" transform=\"translate(0 -14.9975)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"14.9975\"/><path class=\"sym\" d=\"M -3,-3 l 6,6 M -3,3 l 6,-6\"/><text class=\"lab lft\" transform=\"translate(-5 0)\">3</text></g><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">3</text></g><g class=\"a11 i8 leaf node n8 p0 sample\" transform=\"translate(21.2889 35.5145)\"><path class=\"edge\" d=\"M 0 0 V -35.5145 H -21.2889\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">8</text></g><path class=\"edge\" d=\"M 0 0 V -2.7373 H -31.9333\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">11</text></g><path class=\"edge\" d=\"M 0 0 V -193.643 H 74.5111\"/><g class=\"mut m2 s0\" transform=\"translate(0 -30.684)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"30.684\"/><path class=\"sym\" d=\"M -3,-3 l 6,6 M -3,3 l 6,-6\"/><text class=\"lab lft\" transform=\"translate(-5 0)\">2</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">12</text></g><g class=\"a15 c2 m1 node n13 p0 s0\" transform=\"translate(74.5111 168.633)\"><g class=\"a13 i7 leaf node n7 p0 sample\" transform=\"translate(31.9333 63.2615)\"><path class=\"edge\" d=\"M 0 0 V -63.2615 H -31.9333\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">7</text></g><g class=\"a13 c2 node n10 p0\" transform=\"translate(-31.9333 55.2254)\"><g class=\"a10 i4 leaf node n4 p0 sample\" transform=\"translate(-21.2889 8.03611)\"><path class=\"edge\" d=\"M 0 0 V -8.03611 H 21.2889\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">4</text></g><g class=\"a10 i6 leaf node n6 p0 sample\" transform=\"translate(21.2889 8.03611)\"><path class=\"edge\" d=\"M 0 0 V -8.03611 H -21.2889\"/><rect class=\"sym\" height=\"6\" width=\"6\" x=\"-3\" y=\"-3\"/><text class=\"lab\" transform=\"translate(0 11)\">6</text></g><path class=\"edge\" d=\"M 0 0 V -55.2254 H 31.9333\"/><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab lft\" transform=\"translate(-3 -7.0)\">10</text></g><path class=\"edge\" d=\"M 0 0 V -168.633 H -74.5111\"/><g class=\"mut m1 s0\" transform=\"translate(0 -27.3696)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"27.3696\"/><path class=\"sym\" d=\"M -3,-3 l 6,6 M -3,3 l 6,-6\"/><text class=\"lab rgt\" transform=\"translate(5 0)\">1</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">13</text></g><path class=\"edge\" d=\"M 0 0 V -162.705 H -90.4778\"/><g class=\"mut m0 s0\" transform=\"translate(0 -119.177)\"><line x1=\"0\" x2=\"0\" y1=\"0\" y2=\"119.177\"/><path class=\"sym\" d=\"M -3,-3 l 6,6 M -3,3 l 6,-6\"/><text class=\"lab rgt\" transform=\"translate(5 0)\">0</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab rgt\" transform=\"translate(3 -7.0)\">15</text></g><circle class=\"sym\" cx=\"0\" cy=\"0\" r=\"3\"/><text class=\"lab\" transform=\"translate(0 -11)\">16</text></g></g></g></g></g></svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#SOL\n",
    "display(SVG(ts_with_muts_1.draw_svg(size=(500, 500),y_axis=True)), \n",
    "        SVG(ts_with_muts_2.draw_svg(size=(500, 500),y_axis=True)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default `sim_ancestry` and `sim_mutations` use different random seeds. This means that every ancestry and mutation simulation will be different. It is useful to be able to replicate the same simulation; setting the same random seed guarantees this. Try this by rerunning the code chunk above!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.5 Using `msprime` to find the distribution of genealogies and coalescent times"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have so far only simulated single genealogies. We are usually interested in using coalescent simulations to understand the **distribution** of genealogies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> The <code>num_replicates</code> argument of <code>sim_ancestry()</code> provides a convenient and efficient way to iterate over a number of replicate simulations. Let's use this to examine the distribution of the time to the most recent common ancestor.\n",
    "        </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_replicates = 1_000\n",
    "tmrca_msprime = np.zeros(num_replicates)\n",
    "replicates = msprime.sim_ancestry(9, ploidy=1, num_replicates=num_replicates, random_seed=42)\n",
    "for replicate_index, ts in enumerate(replicates):\n",
    "    tree = ts.first()\n",
    "    tmrca_msprime[replicate_index] = tree.time(tree.root)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the output of `sim_ancestry()` is a generator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "generator"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(replicates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> We will use <code>for</code> loops like the one above for populating numpy arrays with simulation results throughout the course, so it's worth understanding the syntax used: <code>np.zeros()</code> sets up an empty <code>numpy</code> array. It’s important to note that the replicate simulations are generated lazily here on demand — the set of simulations we have specified using <code>sim_ancestry</code> and assigned to the <code>replicates</code> variable is also a <code>generator</code>. This means that simulations are only specified but not actually run when we define <code>replicates</code>. \n",
    "\n",
    "Inside the <code>for</code> loop we use the built-in <code>enumerate()</code> function to easily get the index of each replicate in the list. Each replicate simulation is done sequentially as we go through the loop. We therefore only use the memory that we need to represent the results of one simulation.\n",
    "        </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Note:</b> the type of  <code>sim_ancestry()</code> changes when we use the  <code>num_replicates </code> argument. If num_replicates is not specified or None, we return an instance of <code>tskit.TreeSequence</code>. If it is specified, we return an iterator over a set of <code>tskit.TreeSequence </code> instances.\n",
    "        </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 12:</b>  What is the mean $t_{mrca}$? How does the distribution compare to the one we obtained using our own coalescence time simulator <code>coaltimes</code> at the start of this practical?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean Tmrca is 35990.76569596264\n",
      "The mean Tmrca is 36597.53044624645\n"
     ]
    }
   ],
   "source": [
    "#SOL\n",
    "print(\"The mean Tmrca is\", tmrca_msprime.mean()*20000)\n",
    "print(\"The mean Tmrca is\", tmrca.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFuCAYAAAC/a8I8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWvUlEQVR4nO3df7BndX3f8eeLXYhCZJCyjdxddCkh2m2LP3KD8ccka4ydJTGikSg0ao06SFs0NpOJ2HEsrZOORtvGKmazQURTI1ERiroBE+Kixoi7EFBAcTaAsoLZBSOIpcDiu3+cc/Hr9e7u9+69Zz/fu/f5mDnD+fX9fN/3O3dffO75ns/npKqQJB14h7QuQJKWKwNYkhoxgCWpEQNYkhoxgCWpkZWtC5ivDRs21OWXX966jOZe+sd/C8Cfv/YZjSuRNIbMtXPJ9YDvuuuu1iVI0qJYcgEsSQcLA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGA52nt1BRJFrysnZpq/aNIamzJTcje2jfuvJNav37B7WTLlgW3IWlpswcsSY0YwJLUiAEsSY0YwJLUyKABnGRDkpuTbE9yzhzH1ye5J8l1/fKWIeuRpEky2F0QSVYA5wHPA3YAW5NcVlU3zTr1c1X1/KHqkKRJNWQP+GRge1XdUlUPAhcBpw74fpK0pAwZwKuB20e2d/T7ZntGkuuT/EWSfzFXQ0nOTLItybZdu3YNUaskHXBDBnDm2Feztq8FnlBVTwbeDVw6V0NVtamqpqtqetWqVYtbpSQ1MmQA7wCOG9leA9wxekJV3VtV9/Xrm4FDkxwzYE2SNDGGDOCtwIlJjk9yGHA6cNnoCUkelyT9+sl9PXcPWJMkTYzB7oKoqt1JzgauAFYAF1TVjUnO6o9vBE4D/l2S3cD9wOlVNfsyhSQdlAadjKe/rLB51r6NI+vvAd4zZA2SNKkcCSdJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIAdzIoUCS/V6uuuoqrrrqKpKwdmqq9Y8jaT+sbF3AcvUQUOvX7/frX3rUUQD8+fr1ZMuWxShJ0gFmD1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJamRQQM4yYYkNyfZnuScvZz3c0keTnLakPVI0iQZLICTrADOA04B1gFnJFm3h/PeDlwxVC2SNImG7AGfDGyvqluq6kHgIuDUOc57HXAxsHPAWiRp4gwZwKuB20e2d/T7HpFkNfAiYOOAdUjSRBoygDPHvpq1/YfAG6vq4b02lJyZZFuSbbt27Vqs+iSpqSGfiLEDOG5kew1wx6xzpoGLkgAcA/xKkt1VdenoSVW1CdgEMD09PTvEJWlJGjKAtwInJjke+BZwOvBvRk+oquNn1pNcCHxydvhK0sFqsACuqt1Jzqa7u2EFcEFV3ZjkrP64130lLWuDPpSzqjYDm2ftmzN4q+qVQ9YiSZPGkXCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNLJsAXjs1RZIFL5K0WFa2LuBA+cadd1Lr1y+4nWzZsuA2JAmWUQ9YkiaNASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktTIoAGcZEOSm5NsT3LOHMdPTfLlJNcl2Zbk2UPWI0mTZLChyElWAOcBzwN2AFuTXFZVN42cdiVwWVVVkpOAjwBPGqomSZokQ/aATwa2V9UtVfUgcBFw6ugJVXVfVVW/eQRQSNIyMWQArwZuH9ne0e/7EUlelORrwKeAV83VUJIz+0sU23bt2jVIsZJ0oA0ZwHPN3fhjPdyquqSqngS8EHjrXA1V1aaqmq6q6VWrVi1ulZLUyJABvAM4bmR7DXDHnk6uqs8CJyQ5ZsCaJGliDBnAW4ETkxyf5DDgdOCy0ROS/HT6Wc6TPA04DLh7wJokaWIMdhdEVe1OcjZwBbACuKCqbkxyVn98I/Bi4BVJHgLuB1468qWcJB3UBn0iRlVtBjbP2rdxZP3twNuHrEGSJtVYlyCSXDnOPknS+PbaA07yKOBw4Jgkj+WHdzYcCUwNXJskHdT2dQnitcAb6ML2Gn4YwPfSjXKTJO2nvQZwVb0LeFeS11XVuw9QTZK0LIz1JVxVvTvJM4G1o6+pqg8OVJckHfTGCuAkfwqcAFwHPNzvLsAAngCHAv3t1AvyhGOP5bY79jhWRtIiG/c2tGlgnffoTqaHgFq/fsHtZMuWBbchaXzjjoS7AXjckIVI0nIzbg/4GOCmJF8CHpjZWVUvGKQqSVoGxg3gc4csQpKWo3Hvgrhq6EIkabkZ9y6I7/HDuXwPo/vi/ftVdeRQhUnSwW7cHvBjRreTvJDukUOSpP20X/MBV9WlwC8tbimStLyMewni10c2D6G7L9h7giVpAca9C+LXRtZ3A7cx6wnHkqT5Gfca8G8NXYgkLTfjTsi+JsklSXYm+YckFydZM3RxknQwG/dLuPfTPVBzClgNfKLfJ0naT+MG8Kqqen9V7e6XC4FVA9YlSQe9cQP4riQvS7KiX16Gj4+XpAUZN4BfBbwE+DZwJ3Aa4BdzkrQA496G9lbg31bVPwIkORp4J10wS5L2w7g94JNmwhegqr4DPHWYkiRpeRg3gA/pH0sPPNIDHrf3LEmaw7gh+t+BLyT5GN0Q5JcAvz9YVZK0DIw7Eu6DSbbRTcAT4Ner6qZBK5Okg9zYlxH6wDV0JWmR7Nd0lJKkhTOAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJamRQQM4yYYkNyfZnuScOY7/ZpIv98sXkjx5yHokaZIMFsBJVgDnAacA64AzkqybddqtwC9W1UnAW4FNQ9UjSZNmyB7wycD2qrqlqh4ELgJOHT2hqr5QVf/Yb34RWDNgPZI0UYYM4NXA7SPbO/p9e/Jq4C/mOpDkzCTbkmzbtWvXIpYoSe0MGcCZY1/NeWLyHLoAfuNcx6tqU1VNV9X0qlWrFrFESWpn5YBt7wCOG9leA9wx+6QkJwHnA6dU1d0D1iNJE2XIHvBW4MQkxyc5DDgduGz0hCSPBz4OvLyqvj5gLZI0cQbrAVfV7iRnA1cAK4ALqurGJGf1xzcCbwH+CfDeJAC7q2p6qJokaZIMeQmCqtoMbJ61b+PI+muA1wxZgyRNKkfCSVIjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAecSiQZMHL2qmp1j+KtCQM+kgiLS0PAbV+/YLbyZYtC25DWg7sAUtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSI4MGcJINSW5Osj3JOXMcf1KSv03yQJLfHbIWSZo0gwVwkhXAecApwDrgjCTrZp32HeD1wDuHqkMH3qFAkgUva6emWv8o0qBWDtj2ycD2qroFIMlFwKnATTMnVNVOYGeSXx2wDh1gDwG1fv2C28mWLQtuQ5pkQ16CWA3cPrK9o983b0nOTLItybZdu3YtSnGS1NqQAZw59tX+NFRVm6pquqqmV61atcCyJGkyDBnAO4DjRrbXAHcM+H6StKQMGcBbgROTHJ/kMOB04LIB30+SlpTBvoSrqt1JzgauAFYAF1TVjUnO6o9vTPI4YBtwJPCDJG8A1lXVvUPVJUmTYsi7IKiqzcDmWfs2jqx/m+7ShCQtO46Ek6RGDGBJasQAlqRGDGBJasQAlqRGDGBJasQAlqRGDGBJasQAlqRGDGBJasQAlqRGDGBJasQA1sRarGfL+Xw5TapBZ0OTFmKxni0HPl9Ok8kesCQ1YgBLUiMGsCQ1YgBLUiMGsCQ1YgBLUiMGsCQ1YgBLUiMGsCQ1YgBLUiMGsCQ1YgBLUiMGsCQ1YgBrWVisqS2d1lKLyekotSws1tSWTmupxWQPWJIaMYAlqREDWJIaMYAlqREDWJIaMYAlqREDWJIaMYCleXBAhxaTAzGkeXBAhxaTPWBJasQAlqRGDGBJasQAlqRGDGBJasQAlqRGDGBJasQAlqRGDGBJasQAlhpYrCHNDmte2hyKLDWwWEOawWHNS5k9YGmJm7QJgtZOTU1UPZPMHrC0xC1Wb/qwLVtIsuB2wAmLxjVoACfZALwLWAGcX1Vvm3U8/fFfAf4v8MqqunbImiTNzZneDrzBLkEkWQGcB5wCrAPOSLJu1mmnACf2y5nAHw1VjyRNmiGvAZ8MbK+qW6rqQeAi4NRZ55wKfLA6XwSOSnLsgDVJWoYm9bp0qmpRG3yk4eQ0YENVvabffjnw9Ko6e+ScTwJvq6rP99tXAm+sqm2z2jqTrocM8ETgbuCuQQrfP8cwWfXA5NU0afXA5NU0afWANY1jnHruqqoNs3cOeQ14rqv5s9N+nHOoqk3ApkdelGyrqumFlbd4Jq0emLyaJq0emLyaJq0esKZxLKSeIS9B7ACOG9leA9yxH+dI0kFpyADeCpyY5PgkhwGnA5fNOucy4BXp/DxwT1XdOWBNkjQxBrsEUVW7k5wNXEF3G9oFVXVjkrP64xuBzXS3oG2nuw3tt8ZsftO+TzmgJq0emLyaJq0emLyaJq0esKZx7Hc9g30JJ0naO4ciS1IjBrAkNbKkAjjJhiQ3J9me5JwJqOeCJDuT3NC6FoAkxyX5TJKvJrkxyW9PQE2PSvKlJNf3Nf2X1jVBN1Izyd/196I3l+S2JF9Jcl2Sbft+xfCSHJXkY0m+1v9OPaNhLU/sP5uZ5d4kb2hVz0hd/7H/vb4hyYeTPGper18q14D7oc1fB55Hd/vaVuCMqrqpYU2/ANxHN5rvX7aqY6SeY4Fjq+raJI8BrgFe2PgzCnBEVd2X5FDg88Bv9yMfm0nyO8A0cGRVPb9lLX09twHTVTUxAwySfAD4XFWd39/JdHhVfbdxWTNZ8C26gV3faFjHarrf53VVdX+SjwCbq+rCcdtYSj3gcYY2H1BV9VngOy1rGFVVd85MZlRV3wO+CqxuXFNV1X395qH90vT/+knWAL8KnN+yjkmW5EjgF4D3AVTVg5MQvr3nAn/fMnxHrAQenWQlcDjzHMewlAJ4NXD7yPYOGofLJEuyFngqcHXjUmb+3L8O2An8ZVW1rukPgd8DftC4jlEFfDrJNf3Q+9b+GbALeH9/qeb8JEe0Lqp3OvDh1kVU1beAdwLfBO6kG8fw6fm0sZQCeKxhy4IkPwlcDLyhqu5tXU9VPVxVT6Eb6XhykmaXa5I8H9hZVde0qmEPnlVVT6ObIfA/9Je3WloJPA34o6p6KvB9YBK+dzkMeAHw0Qmo5bF0f4UfD0wBRyR52XzaWEoB7LDlMfTXWS8GPlRVH29dz6j+T9gtwI9NSnIAPQt4QX/N9SLgl5L874b1AFBVd/T/3QlcQnfJraUdwI6Rv1Y+RhfIrZ0CXFtV/9C6EOCXgVuraldVPQR8HHjmfBpYSgE8ztDmZa3/wut9wFer6n+0rgcgyaokR/Xrj6b7pf1aq3qq6k1Vtaaq1tL9Dv11Vc2r17LYkhzRf2lK/2f+vwaa3llTVd8Gbk/yxH7Xc4FmX+aOOIMJuPzQ+ybw80kO7//tPZfue5exLZlHEu1paHPLmpJ8GFgPHJNkB/Cfq+p9DUt6FvBy4Cv9NVeA/1RVm9uVxLHAB/pvrg8BPlJVE3Hr1wT5KeCS7t8wK4E/q6rL25YEwOuAD/UdnlsYf6qAQSQ5nO4uqNe2rGNGVV2d5GPAtcBu4O+Y57DkJXMbmiQdbJbSJQhJOqgYwJLUiAEsSY0YwJLUiAEsSY0YwJqXfoasfz+yPdXfinOg3v/YmRnMkkwn+V8H6r0XW5JXJpka2T4/ybpFbP9fJblwsdrT4vM2NM1LP8fEJ1vN/pbkHcDnq+r/tHj/+Uqyoqoe3sOxLcDvVtVg008m+SvgVVX1zaHeQ/vPHrDm623ACf2crO9IsnZmPuS+R3dpkk8kuTXJ2Ul+p5/M5YtJju7POyHJ5f3EM59L8qR+/2/086pen+Sze3j/FwOX9+evH+kNn5tufuYtSW5J8vq5Xpzk1Um+3p/3J0ne0+9fleTiJFv75Vn7ajfJy9LNdXxdkj/uB5uQ5L4k/zXJ1cAzkrylb/OGJJvSOY1uOswP9a9/dP8e030bZ6SbH/iGJG8fec/7kvx+/xl9MclP7eOz+wTdiD9NoqpycRl7AdYCN8y1DbyS7gGrjwFWAfcAZ/XH/ifd5EAAVwIn9utPpxsODPAVYHW/ftQc7308cM3I9nq63jjAucAXgJ8AjgHuBg6d9fop4DbgaLppMT8HvKc/9mfAs/v1x9MN595ju8A/pwu3Q/vz3gu8ol8v4CUj73v0yPqfAr/Wr2+hmwOY0e2+zm/2n+FK4K/p5nWeaXvm9X8AvHlvnx3d6MhPtP69cZl7WTJDkbVkfKa6uYi/l+QeupCCLiBOSjdT2zOBj/ZDb6ELN4C/AS5MN7H1XBMJHUs3ReKefKqqHgAeSLKTbojvjpHjJwNXVdV3AJJ8FPiZ/tgvA+tGajpyZn6GPbT7XOBnga39ax5NN90mwMN0EyLNeE6S36ObL/Zo4MaRz2UuPwdsqapdfZ0fopub91LgQWBmKPc1dENzYc+f3U66QNcEMoC12B4YWf/ByPYP6H7fDgG+W930lD+iqs5K8nS6ydKvS/KUqrp75JT7gb098mX0vR/mx3+/55rSdMYhwDOq6v4feUEXrnO1G+ADVfWmOdr6f9Vf9033iJr30vV0b09y7j5+hn3V+VD1XduRWvb22T2K7nPTBPIasObre3SXGPZLdfMT35rkN6CbwS3Jk/v1E6rq6qp6C3AXPzr9KHSPpFq7v+8NfAn4xSSPTfcEgxePHPs0cPbMRpKn7KOtK4HTkvzT/vyjkzxhjvNmwvauvvd/2sixPX2WV/d1HtNfVz4DuGpvxezls/sZGs+spj0zgDUvfa/qb/ovfN6xn838JvDqJNfT/Tk+82ipd8x88QR8Frh+1nt/H/j7JD+9n7V/C/hvdAH3V3TTK97TH349MJ3ky0luAs7aR1s3AW+me4rFl4G/pLtEMvu87wJ/QncJ5lK6aVVnXAhsnPkSbuQ1dwJvAj5D9xlcW/u+62NPn91zgE/t47VqxNvQtKQkeRHws1X15v18/U9W94DQlXQTn19QVZcsapETIslP0PWcn11Vu1vXox9nD1hLSh+Wty2giXPTzZV8A3ArXa/0YPV44BzDd3LZA5akRuwBS1IjBrAkNWIAS1IjBrAkNWIAS1Ij/x90LhaxkeRbvgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#SOL\n",
    "sns.displot(x = tmrca_msprime, stat = \"density\", label = 'span of trees', binwidth = 0.5, color = \"r\")\n",
    "plt.axvline(x= np.mean(tmrca_msprime)) \n",
    "plt.xlabel(\"times (in generations)\");\n",
    "plt.ylabel(\"count\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFuCAYAAAC/a8I8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWbklEQVR4nO3df7DddX3n8ecrCRS1MsiSVsyPwlKsm93FX9ewYrdFLTuhtY22tECrrlUnZHejdTqdijuOy67THV13t7WCDSki2toFK0KjptCWimhTNMFFJCBOFhlzG7sJWFFcFwy894/zDXu83CTnXu43n3PvfT5mvsP3x+d8vu9z5ubF537P9/u5qSokSUffktYFSNJiZQBLUiMGsCQ1YgBLUiMGsCQ1sqx1ATO1bt26uuGGG1qXIUkzkel2zrsR8P3339+6BEmaE/MugCVpoTCAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGuk1gJOsS3JPkt1JLj5Em7OT3J5kV5LP9FmPJI2T3p6ES7IUuAw4B5gEdiTZWlV3DbU5AXg/sK6qvp7kR/qqR5LGTZ8j4LXA7qq6t6oeAa4G1k9p86vAx6vq6wBVta/HeiRprPQZwCuAPUPbk92+Yc8GnpHk5iS3JXltj/VI0ljpczKe6SafmPr3j5YBLwReDjwF+Nskt1bVV3+go2QDsAFg9erVPZQqSUdfnyPgSWDV0PZKYO80bW6oqu9W1f3ALcBzp3ZUVVuqaqKqJpYvX95bwZJ0NPUZwDuA05OcmuRY4AJg65Q2fwb8yyTLkjwVOBO4u8eaJGls9BbAVXUA2ATcyCBUP1pVu5JsTLKxa3M3cANwB/AF4IqqurOvmhaCFatWk6TXZcUqL/NIR0Pm25+ln5iYqJ07d7Yuo5kknH/59l7Pcc1FZzHffi6kMbcwJmSXpIXCAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWrEAJakRgxgSWqk1wBOsi7JPUl2J7l4muNnJ3kwye3d8o4+65GkcbKsr46TLAUuA84BJoEdSbZW1V1Tmn62ql7RVx2SNK76HAGvBXZX1b1V9QhwNbC+x/NJ0rzSZwCvAPYMbU92+6Z6cZIvJfnzJP90uo6SbEiyM8nO/fv391GrJB11fQZwptlXU7a/CPxYVT0XeB9w/XQdVdWWqpqoqonly5fPbZWS1EifATwJrBraXgnsHW5QVd+uqoe69W3AMUlO6rEmSRobfQbwDuD0JKcmORa4ANg63CDJM5OkW1/b1fNAjzVJ0tjo7S6IqjqQZBNwI7AUuLKqdiXZ2B3fDJwH/JskB4DvARdU1dTLFJK0IPUWwPD4ZYVtU/ZtHlq/FLi0zxokaVz5JJwkNWIAS1IjBrAkNWIAz6EVq1aTpNdF0sLR65dwi83eyT2cf/n2Xs9xzUVn9dq/pKPHEbAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNWIAS1IjBrAkNdJrACdZl+SeJLuTXHyYdi9K8miS8/qsR5LGSW8BnGQpcBlwLrAGuDDJmkO0ezdwY1+1SNI46nMEvBbYXVX3VtUjwNXA+mnavQm4FtjXYy2aiSXLSNLrsmLV6tbvUmpuWY99rwD2DG1PAmcON0iyAngV8DLgRT3Wopl47ADnX76911Ncc9FZvfYvzQd9joAzzb6asv17wFur6tHDdpRsSLIzyc79+/fPVX2S1FSfI+BJYNXQ9kpg75Q2E8DVSQBOAn42yYGqun64UVVtAbYATExMTA1xSZqX+gzgHcDpSU4F/g64APjV4QZVderB9SRXAZ+cGr6StFD1FsBVdSDJJgZ3NywFrqyqXUk2dsc393VuSZoP+hwBU1XbgG1T9k0bvFX1uj5rkaRx45NwktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktSIASxJjRjAktTISAGc5KZR9kmSRrfscAeTHAc8FTgpyTOAdIeOB57Vc22StKAdNoCBi4C3MAjb2/j/Afxt4LL+ypKkhe+wAVxV7wXem+RNVfW+o1STJC0KRxoBA1BV70tyFnDK8Guq6sM91SVJC95IAZzkj4DTgNuBR7vdBRjAkjRLIwUwMAGsqarqsxhJWkxGvQ/4TuCZfRYiSYvNqCPgk4C7knwBePjgzqr6hV6qkqRFYNQAvqTPIiRpMRr1LojP9F2IJC02o94F8R0Gdz0AHAscA3y3qo7vqzBJWuhGHQE/fXg7ySuBtX0UJEmLxaxmQ6uq64GXzW0pkrS4jHoJ4heHNpcwuC/4iPcEJ1kHvBdYClxRVe+acnw98E7gMeAA8Jaq+txopUvS/DbqXRA/P7R+ALgPWH+4FyRZymDCnnOASWBHkq1VdddQs5uArVVVSc4APgo8Z8SaJGleG/Ua8K/Pou+1wO6quhcgydUMQvvxAK6qh4baP40RRtWStFCMOiH7yiTXJdmX5H8nuTbJyiO8bAWwZ2h7sts3te9XJfkK8Cng9Yc4/4YkO5Ps3L9//yglS9LYG/VLuA8CWxnMC7wC+ES373Ayzb4njHCr6rqqeg7wSgbXg5/4oqotVTVRVRPLly8fsWRJGm+jBvDyqvpgVR3olquAIyXhJLBqaHslsPdQjavqFuC0JCeNWJMkzWujBvD9SV6dZGm3vBp44Aiv2QGcnuTUJMcCFzAYRT8uyY8nSbf+AgYPeRypX0laEEa9C+L1wKXA7zK4jLAdOOwXc1V1IMkm4EYGt6FdWVW7kmzsjm8Gfgl4bZLvA98DznfKS0mLxagB/E7gX1fVPwAkORH4rxziS7ODqmobsG3Kvs1D6+8G3j2TgiVpoRj1EsQZB8MXoKq+CTy/n5IkaXEYNYCXdH+WHnh8BDzq6FmSNI1RQ/S/AduTfIzBNeBfAX6nt6okaREY9Um4DyfZyWACngC/OOWRYknSDI18GaELXENXkubIrKajlCQ9eQawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAEtSIwawJDViAKuNJctI0tuyYtXq1u9QOqJlrQvQIvXYAc6/fHtv3V9z0Vm99S3NFUfAktSIASxJjRjAktSIASxJjRjAktSIASxJjfQawEnWJbknye4kF09z/NeS3NEt25M8t896JGmc9BbASZYClwHnAmuAC5OsmdLsa8BPV9UZwDuBLX3VI0njps8R8Fpgd1XdW1WPAFcD64cbVNX2qvqHbvNWYGWP9UjSWOkzgFcAe4a2J7t9h/IG4M+nO5BkQ5KdSXbu379/DkuUpHb6DOBMs6+mbZi8lEEAv3W641W1paomqmpi+fLlc1iiJLXT51wQk8Cqoe2VwN6pjZKcAVwBnFtVD/RYjySNlT5HwDuA05OcmuRY4AJg63CDJKuBjwOvqaqv9liLJI2d3kbAVXUgySbgRmApcGVV7UqysTu+GXgH8I+A9ycBOFBVE33VJEnjpNfpKKtqG7Btyr7NQ+tvBN7YZw2SNK58Ek6SGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJakRA1iSGjGAJamRRRPAK1atJkmviyTNxLLWBRwteyf3cP7l23s9xzUXndVr/5IWlkUzApakcWMAS1IjBrAkNWIAS1IjBrAkNWIAS1IjvQZwknVJ7kmyO8nF0xx/TpK/TfJwkt/qsxZJGje93QecZClwGXAOMAnsSLK1qu4aavZN4M3AK/uqQ5LGVZ8j4LXA7qq6t6oeAa4G1g83qKp9VbUD+H6PdWgxWrKs9ycfV6xa3fpdap7r80m4FcCeoe1J4MzZdJRkA7ABYPVqf+g1gscO+OSjxl6fI+DpJkeo2XRUVVuqaqKqJpYvX/4ky5Kk8dBnAE8Cq4a2VwJ7ezyfJM0rfQbwDuD0JKcmORa4ANja4/kkaV7p7RpwVR1Isgm4EVgKXFlVu5Js7I5vTvJMYCdwPPBYkrcAa6rq233VJUnjotfpKKtqG7Btyr7NQ+t/z+DShCQtOj4JJ0mNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLEmNGMDSbC1ZRpJelxWrVrd+l+rRstYFSPPWYwc4//LtvZ7imovO6rV/teUIWJIaMYAlqREDWJIaMYAlqREDWJIaMYAlqREDWBpn3mu8oHkfsDTOvNd4QXMELEmNGMCS1IgBLEmNGMCS1IgBLEmNGMCS1IgBLC123mvcjPcBS4ud9xo30+sIOMm6JPck2Z3k4mmOJ8nvd8fvSPKCPuuRpHHSWwAnWQpcBpwLrAEuTLJmSrNzgdO7ZQPwB33VI0njps8R8Fpgd1XdW1WPAFcD66e0WQ98uAZuBU5IcnKPNUlaoFasWj3vrmWnqua0w8c7Ts4D1lXVG7vt1wBnVtWmoTafBN5VVZ/rtm8C3lpVO6f0tYHBCBngJ4B7eim6HycB97cuYg74PsaL72O8HOl93F9V66bu7PNLuEyzb2raj9KGqtoCbJmLoo62JDuraqJ1HU+W72O8+D7Gy2zfR5+XICaBVUPbK4G9s2gjSQtSnwG8Azg9yalJjgUuALZOabMVeG13N8S/AB6sqm/0WJMkjY3eLkFU1YEkm4AbgaXAlVW1K8nG7vhmYBvws8Bu4P8Av95XPQ3Ny0sn0/B9jBffx3iZ1fvo7Us4SdLh+SiyJDViAEtSIwZwT5JcmWRfkjtb1/JkJFmV5NNJ7k6yK8lvtK5pNpIcl+QLSb7UvY//2Lqm2UqyNMn/7O6jn5eS3Jfky0luT7LzyK8YT0lOSPKxJF/p/o28eEav9xpwP5L8FPAQgyf9/lnremarezLx5Kr6YpKnA7cBr6yquxqXNiNJAjytqh5KcgzwOeA3uicw55UkvwlMAMdX1Sta1zMbSe4DJqpqXj+EkeRDwGer6orubq+nVtW3Rn29I+CeVNUtwDdb1/FkVdU3quqL3fp3gLuBFW2rmrnucfeHus1jumXejT6SrAR+DriidS2LXZLjgZ8CPgBQVY/MJHzBANYMJDkFeD7w+calzEr3q/vtwD7gL6tqPr6P3wN+G3iscR1PVgF/keS2bqqB+egfA/uBD3aXhK5I8rSZdGAAayRJfhi4FnhLVX27dT2zUVWPVtXzGDxxuTbJvLo0lOQVwL6quq11LXPgJVX1AgYzIv677pLdfLMMeAHwB1X1fOC7wBOm3T0cA1hH1F0zvRb4SFV9vHU9T1b3a+LNwBMmRxlzLwF+obt+ejXwsiR/3Lak2amqvd1/9wHXMZg9cb6ZBCaHfpP6GINAHpkBrMPqvrz6AHB3Vf331vXMVpLlSU7o1p8C/AzwlaZFzVBVva2qVlbVKQwe7f/rqnp147JmLMnTui906X5l/1fAvLtbqKr+HtiT5Ce6XS8HZvTltH+SqCdJ/gdwNnBSkkngP1TVB9pWNSsvAV4DfLm7fgrw76tqW7uSZuVk4EPdHwpYAny0qubtbVzz3I8C1w3+384y4E+q6oa2Jc3am4CPdHdA3MsMp1PwNjRJasRLEJLUiAEsSY0YwJLUiAEsSY0YwJLUiAGsGelmf/q3Q9vPSvKxo3j+kw/OApZkIsnvH61zz7Ukr0vyrKHtK5KsmcP+/3mSq+aqP809b0PTjHTzQXyy1QxvSd4DfK6q/qzF+WcqydKqevQQx24GfquqepuOMclfAa+vqq/3dQ7NniNgzdS7gNO6eVzfk+SUg3MedyO665N8IsnXkmxK8pvdRCW3Jjmxa3dakhu6iVg+m+Q53f5fTnJnN2fvLYc4/y8BN3Ttzx4aDV/SzcF8c5J7k7x5uhcneUOSr3bt/jDJpd3+5UmuTbKjW15ypH6TvLqbY/j2JJd3D3mQ5KEk/ynJ54EXJ3lH1+edSbZk4DwGU0p+pHv9U7pzTHR9XJjBfLl3Jnn30DkfSvI73Wd0a5IfPcJn9wkGT81pHFWVi8vIC3AKcOd028DrGPyB1acDy4EHgY3dsd9lMJEPwE3A6d36mQweqQX4MrCiWz9hmnOfCtw2tH02g9E4wCXAduCHgJOAB4Bjprz+WcB9wIkMpqP8LHBpd+xPgJ/s1lczePT6kP0C/4RBuB3TtXs/8NpuvYBfGTrviUPrfwT8fLd+M4M5cRne7ur8evcZLgP+msEczAf7Pvj6/wK8/XCfHYMnGT/R+ufGZfrFR5E11z5dg3mDv5PkQQYhBYOAOKObVe0s4E+7R1FhEG4AfwNcleSjwHST/pzMYPq/Q/lUVT0MPJxkH4NHXieHjq8FPlNV3wRI8qfAs7tjPwOsGarp+IPzFRyi35cDLwR2dK95CoNpLgEeZTB50UEvTfLbwFMZhP+uoc9lOi8Cbq6q/V2dH2Ew7+z1wCPAwUeobwPO6dYP9dntYxDoGkMGsObaw0Prjw1tP8bg520J8K0aTAv5A6pqY5IzGUw4fnuS51XVA0NNvgccN+K5H+WJP9/h0JYAL66q7/3ACwbhOl2/AT5UVW+bpq//W9113yTHMRgdT1TVniSXHOE9HKnO71c3tB2q5XCf3XEMPjeNIa8Ba6a+w+ASw6zUYC7hryX5ZRjMtpbkud36aVX1+ap6B3A/sGrKy7/K4JLHbH0B+Okkz0iyjMH15IP+Ath0cCPJ847Q103AeUl+pGt/YpIfm6bdwbC9vxv9nzd07FCf5ee7Ok/qritfCHzmcMUc5rN7NvNwprHFwgDWjHSjqr/pvvB5zyy7+TXgDUm+xODX8fXd/vcc/OIJuAX40pRzfxf4X0l+fJa1/x3wnxkE3F8xmDrwwe7wm4GJJHckuQvYeIS+7gLezuCvOtwB/CWDSyRT230L+EMGl2CuB3YMHb4K2HzwS7ih13wDeBvwaQafwRfryHd9HOqzeynwqSO8Vo14G5rmlSSvAl5YVW+f5et/uAZ/mHMZg4nAr6yq6+a0yDGR5IcYjJx/sqoOtK5HT+QIWPNKF5b3PYkuLslgXuM7ga8xGJUuVKuBiw3f8eUIWJIacQQsSY0YwJLUiAEsSY0YwJLUiAEsSY38P+cU87V0soZhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "##SOL\n",
    "sns.displot(x = tmrca/20_000, stat = \"density\", label = 'span of trees', binwidth = 0.5)\n",
    "plt.xlabel(\"times (in generations)\");\n",
    "plt.ylabel(\"count\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.6 Counting trees..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have **still** only considered coalescence times and ignored topologies. How random are tree topologies under the standard neutral coalescent? The number of possible tree topologies grows very quickly with sample size and is **enormous** even for relatively small samples.\n",
    "\n",
    "We can generate all trees using `tskit.all_trees(i, span=1)`. <b>Note:</b> by default, this count includes all trees, not just bifurcating trees."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> Use <a href=\"https://seaborn.pydata.org/generated/seaborn.lineplot.html\"><code>sns.lineplot()</code></a> for plotting simple functions and datasets.\n",
    "        </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For a samples of 7 lineages the number of unique tree topologies is 39208\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEGCAYAAABlxeIAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiaklEQVR4nO3de2yc1bnv8e9jO/f7xQmuHUgg4ZKEJBTvNGxOeygBEpKqUG2QUqklqiKlQvSI6mydLeg/u/0jUvmjZR+kDVJaegjs7kIObUXUTArsUE5vNMHQJBPnQgwB7GR8ycWJnYsTe57zx6yBsTOxx5fxjGd+H2k07zzvu95ZLxHz+H3Wmlnm7oiIiJTkugMiIpIflBBERARQQhARkUAJQUREACUEEREJynLdgYGaOXOmz507N9fdEBEZUd57770T7l6ebt+ITQhz586lpqYm190QERlRzOyTq+1TyUhERIB+JAQzKzWzv5vZ78Lr6Wb2ppkdCc/TUo590szqzOywma1Kid9uZtGw7xkzsxAfY2avhPguM5s7hNcoIiIZ6M8dwuPAwZTXTwA73X0BsDO8xswWAuuARcBq4FkzKw1tngM2AgvCY3WIbwBOu/t84GngqQFdjYiIDFhGCcHMqoC1wM9Twg8AW8L2FuDBlPjL7t7h7keBOmC5mVUAk939HU/8XsaLPdokz/UqsDJ59yAiIsMj0zuEfwP+BYinxGa7ewwgPM8K8UqgPuW4hhCrDNs9493auHsncAaY0bMTZrbRzGrMrKalpSXDrouISCb6TAhm9jWg2d3fy/Cc6f6y917ivbXpHnDf7O7V7l5dXp521pSIiAxQJtNO7wS+bmZrgLHAZDP7D6DJzCrcPRbKQc3h+AZgTkr7KuB4iFeliae2aTCzMmAKcGqA1yQiIgPQ5x2Cuz/p7lXuPpfEYPFb7v4tYBuwPhy2HngtbG8D1oWZQ/NIDB7vDmWlNjNbEcYHHunRJnmuh8J76He5RURSxOPOpu0H2H/sTFbOP5gvpv0Y2GpmG4BPgYcB3L3WzLYCB4BO4DF37wptHgVeAMYBO8ID4HngJTOrI3FnsG4Q/RIRKUh/r2/lZ386yi0Vk1lcOWXIz9+vhODubwNvh+2TwMqrHLcJ2JQmXgMsThO/SEgoIiKSXiQaY3RpCfcsnJ2V8+ubyiIiI4C7syMa48sLZjJ57KisvIcSgojICLCnvpXjZy6y5taKrL2HEoKIyAgQicYYVWpZKxeBEoKISN5zdyLRRv7b/JlMGZedchEoIYiI5L29DWc41nohq+UiUEIQEcl7O0K56L6F12T1fZQQRETymLuzPRrjzvkzmTI+e+UiUEIQEclr0WNnaDh9gTWLs1suAiUEEZG8tj0ao6zEuG9R9mYXJSkhiIjkqcTsohj/OH8mU8ePzvr7KSGIiOSp2uNnqT91gbW3ZncwOUkJQUQkT22Pxigtyf7soiQlBBGRPPRZueiGGUybkP1yESghiIjkpdrjZ/nk5PmsfxktlRKCiEge2rE/US5atWh4ykWghCAikneSv110x/UzmD5M5SLIICGY2Vgz221me82s1sx+FOI/NLNjZrYnPNaktHnSzOrM7LCZrUqJ325m0bDvmbCUJmG5zVdCfJeZzc3CtYqIjAgHY20cPXGO+4dpdlFSJncIHcDd7r4UWAasNrMVYd/T7r4sPCIAZraQxBKYi4DVwLNmVhqOfw7YSGKd5QVhP8AG4LS7zweeBp4a9JWJiIxQkWiMEmNYy0WQQULwhPbwclR4eC9NHgBedvcOdz8K1AHLzawCmOzu77i7Ay8CD6a02RK2XwVWJu8eRESKSXJ20YrrZzBz4phhfe+MxhDMrNTM9gDNwJvuvivs+p6Z7TOzX5jZtBCrBOpTmjeEWGXY7hnv1sbdO4EzwIw0/dhoZjVmVtPS0pJJ10VERpTDTW18dOLcsM4uSsooIbh7l7svA6pI/LW/mET55wYSZaQY8JNweLq/7L2XeG9tevZjs7tXu3t1eXl5Jl0XERlRIvtyUy6Cfs4ycvdW4G1gtbs3hUQRB34GLA+HNQBzUppVAcdDvCpNvFsbMysDpgCn+tM3EZGRLvlT18vnTad80vCWiyCzWUblZjY1bI8D7gEOhTGBpG8A+8P2NmBdmDk0j8Tg8W53jwFtZrYijA88AryW0mZ92H4IeCuMM4iIFI0Pmtr5sOUca3NQLgIoy+CYCmBLmClUAmx199+Z2UtmtoxEaedj4LsA7l5rZluBA0An8Ji7d4VzPQq8AIwDdoQHwPPAS2ZWR+LOYN3gL01EZGTZHo1hBqsWD3+5CDJICO6+D7gtTfzbvbTZBGxKE68BFqeJXwQe7qsvIiKFbEc0xvK505k1aWxO3l/fVBYRyQNHmto40tzO2iW5KReBEoKISF5IlotW52B2UZISgohIHohEY/zDddOZNTk35SJQQhARybm65jY+aGpnzTD/dlFPSggiIjkWiTYCcH+OppsmKSGIiORYJBqj+rppzM5huQiUEEREcurDlnYONbbl5LeLelJCEBHJoci+GMCwr32QjhKCiEgORfY3cvt106iYMi7XXVFCEBHJlaMnznEwdjYvykWghCAikjORaCgX5ei3i3pSQhARyZHt+2Lcdu1UvjA19+UiUEIQEcmJj0+c40DsbM5+6jodJQQRkRyI7E/OLlJCEBEpapFojGVzplKZJ+UiUEIQERl2n548z/5jZ3P+20U9KSGIiAyz7Z/NLsqfchFktqbyWDPbbWZ7zazWzH4U4tPN7E0zOxKep6W0edLM6szssJmtSonfbmbRsO+ZsLYyYf3lV0J8l5nNzcK1iojkhR37YyytmsKc6eNz3ZVuMrlD6ADudvelwDJgtZmtAJ4Adrr7AmBneI2ZLSSxJvIiYDXwbFiPGeA5YCOwIDxWh/gG4LS7zweeBp4a/KWJiOSf+lPn2ddwJm++jJaqz4TgCe3h5ajwcOABYEuIbwEeDNsPAC+7e4e7HwXqgOVmVgFMdvd33N2BF3u0SZ7rVWBl8u5BRKSQJL+MNiITAoCZlZrZHqAZeNPddwGz3T0GEJ5nhcMrgfqU5g0hVhm2e8a7tXH3TuAMMCNNPzaaWY2Z1bS0tGR0gSIi+SQSjXFrZf6ViyDDhODuXe6+DKgi8df+4l4OT/eXvfcS761Nz35sdvdqd68uLy/vo9ciIvml/tR59uZpuQj6OcvI3VuBt0nU/ptCGYjw3BwOawDmpDSrAo6HeFWaeLc2ZlYGTAFO9advIiL57vf7Eyuj5dO3k1NlMsuo3Mymhu1xwD3AIWAbsD4cth54LWxvA9aFmUPzSAwe7w5lpTYzWxHGBx7p0SZ5roeAt8I4g4hIwdgejbG4cjLXzsi/chFAWQbHVABbwkyhEmCru//OzN4BtprZBuBT4GEAd681s63AAaATeMzdu8K5HgVeAMYBO8ID4HngJTOrI3FnsG4oLk5EJF8ca73AnvpW/teqm3LdlavqMyG4+z7gtjTxk8DKq7TZBGxKE68Brhh/cPeLhIQiIlKIdoTZRflaLgJ9U1lEZFhEojEWVkxm7swJue7KVSkhiIhk2fHWC7z/aStrl+Tv3QEoIYiIZN2OMLsoX1ZGuxolBBGRLItEY9x8zSSuL5+Y6670SglBRCSLYmcu8N4np/N6MDlJCUFEJIuSX0Zbk+fjB6CEICKSVcly0Q15Xi4CJQQRkaxpOnuRmk9O591COFejhCAikiU7ojHcYe2S/J5dlKSEICKSJZFoIzfOnsj8WZNy3ZWMKCGIiGRB89mLvPvJqbz9qet0lBBERLLg97WNiXKREoKISHHbvi/G/FkTWTB7ZJSLQAlBRGTINbddZPfHI6tcBEoIIiJD7vXaphFXLgIlBBGRIRfZF+OG8gncODv/v4yWKpMlNOeY2R/M7KCZ1ZrZ4yH+QzM7ZmZ7wmNNSpsnzazOzA6b2aqU+O1mFg37nglLaRKW23wlxHeZ2dwsXKuISNadaO9g19GTrLm1gvARN2JkcofQCfyzu98CrAAeM7OFYd/T7r4sPCIAYd86YBGwGng2LL8J8BywkcQ6ywvCfoANwGl3nw88DTw1+EsTERl+v9/fSNwZceMHkEFCcPeYu78fttuAg0BlL00eAF529w53PwrUAcvNrAKY7O7vuLsDLwIPprTZErZfBVbaSEutIiIkfrvo+pkTuPmakTO7KKlfYwihlHMbsCuEvmdm+8zsF2Y2LcQqgfqUZg0hVhm2e8a7tXH3TuAMMCPN+280sxozq2lpaelP10VEsu5kewd/+2hklougHwnBzCYCvwa+7+5nSZR/bgCWATHgJ8lD0zT3XuK9tekecN/s7tXuXl1eXp5p10VEhsXrtU0jtlwEGSYEMxtFIhn80t1/A+DuTe7e5e5x4GfA8nB4AzAnpXkVcDzEq9LEu7UxszJgCnBqIBckIpIrkWiMuTPGc0vFyCsXQWazjAx4Hjjo7j9NiaemwG8A+8P2NmBdmDk0j8Tg8W53jwFtZrYinPMR4LWUNuvD9kPAW2GcQURkRDh17hLvjOByEUBZBsfcCXwbiJrZnhD7AfBNM1tGorTzMfBdAHevNbOtwAESM5Qec/eu0O5R4AVgHLAjPCCRcF4yszoSdwbrBnNRIiLD7fXaRrriPmLLRZBBQnD3P5O+xh/ppc0mYFOaeA2wOE38IvBwX30REclXkWiM62aMZ9EXJue6KwOmbyqLiAzS6XOX+OuHI7tcBEoIIiKD9saBUC4aIUtlXo0SgojIIG2PNjJn+jgWV47cchEoIYiIDErr+Uv8te7EiC8XgRKCiMigvHGgic64j7ifuk5HCUFEZBAi0RhV08Zxa+WUXHdl0JQQREQG6Mz5y/ylQMpFoIQgIjJgbxxo5HLXyP4yWiolBBGRAdqxv5HKqeNYWjXyy0WghCAiMiBnLlzmT0daWHPrNQVRLgIlBBGRAfmvA01c7nLuL5ByESghiIgMSCQa4wtTxnLbnKm57sqQUUIQEemnsxcv86cjJ7i/QGYXJSkhiIj0086DTVzqihfM7KIkJQQRkX7avq+RigIrF4ESgohIv7RdvMwfj7SwevE1lJQUTrkIMltCc46Z/cHMDppZrZk9HuLTzexNMzsSnqeltHnSzOrM7LCZrUqJ325m0bDvmbCUJmG5zVdCfJeZzc3CtYqIDNrOg81c6owXxG8X9ZTJHUIn8M/ufguwAnjMzBYCTwA73X0BsDO8JuxbBywCVgPPmllpONdzwEYS6ywvCPsBNgCn3X0+8DTw1BBcm4jIkItEY1wzeSxfvHZa3wePMH0mBHePufv7YbsNOAhUAg8AW8JhW4AHw/YDwMvu3uHuR4E6YLmZVQCT3f0dd3fgxR5tkud6FVhphTR0LyIFob2jk7c/KMxyEfRzDCGUcm4DdgGz3T0GiaQBzAqHVQL1Kc0aQqwybPeMd2vj7p3AGWBGmvffaGY1ZlbT0tLSn66LiAzazoNNXOosvNlFSRknBDObCPwa+L67n+3t0DQx7yXeW5vuAffN7l7t7tXl5eV9dVlEZEhFojFmTRpD9XWFVy6CDBOCmY0ikQx+6e6/CeGmUAYiPDeHeAMwJ6V5FXA8xKvSxLu1MbMyYApwqr8XIyKSLec6Onn7cAv3F2i5CDKbZWTA88BBd/9pyq5twPqwvR54LSW+Lswcmkdi8Hh3KCu1mdmKcM5HerRJnush4K0wziAikhfeOtRMRwGXiwDKMjjmTuDbQNTM9oTYD4AfA1vNbAPwKfAwgLvXmtlW4ACJGUqPuXtXaPco8AIwDtgRHpBIOC+ZWR2JO4N1g7ssEZGhFYnGKJ80huq503PdlazpMyG4+59JX+MHWHmVNpuATWniNcDiNPGLhIQiIpJvzl/q5A+Hm3n49jmUFmi5CPRNZRGRPr11qJmLlwu7XARKCCIifYpEY8ycOIbl8wq3XARKCCIivTp/qZM/HGph9eLZBV0uAiUEEZFevX24hQuXu1izuLDLRaCEICLSq+3RGDMmjC74chEoIYiIXNWFS128dbCZVYuvoay08D8uC/8KRUQG6P990MyFy10F+VPX6SghiIhcxfZoI9MnjOZLRVAuAiUEEZG0Ll7uYufBJlYtml0U5SJQQhARSevtwy2cv9RV8F9GS6WEICKSRiQaY9r4Udxx/RVLsxQsJQQRkR4+LxcVx+yipOK5UhGRDP3xgxbOFVm5CJQQRESuEInGmDp+FHfcUDzlIlBCEBHp5uLlLv7rYDP3LZzNqCIqF4ESgohIN38+coL2js6iKxdBZkto/sLMms1sf0rsh2Z2zMz2hMealH1PmlmdmR02s1Up8dvNLBr2PROW0SQstflKiO8ys7lDfI0iIhmLRGNMGTeKO+fPzHVXhl0mdwgvAKvTxJ9292XhEQEws4Uklr9cFNo8a2al4fjngI0k1lhekHLODcBpd58PPA08NcBrEREZlI7OLt480MS9RVguggwSgrv/kcQ6x5l4AHjZ3Tvc/ShQByw3swpgsru/4+4OvAg8mNJmS9h+FViZvHsQERlOfz5ygraOzqL57aKeBpMCv2dm+0JJaVqIVQL1Kcc0hFhl2O4Z79bG3TuBM0DaoX0z22hmNWZW09LSMoiui4hcaXs0xuSxZUVZLoKBJ4TngBuAZUAM+EmIp/vL3nuJ99bmyqD7Znevdvfq8vLyfnVYRKQ3lzrjoVx0DaPLiq9cBANMCO7e5O5d7h4HfgYsD7sagDkph1YBx0O8Kk28WxszKwOmkHmJSkRkSPyl7gRtFztZu+SaXHclZwaUEMKYQNI3gOQMpG3AujBzaB6JwePd7h4D2sxsRRgfeAR4LaXN+rD9EPBWGGcQERk226MxJo0p3nIRQFlfB5jZr4C7gJlm1gD8K3CXmS0jUdr5GPgugLvXmtlW4ADQCTzm7l3hVI+SmLE0DtgRHgDPAy+ZWR2JO4N1Q3BdIiIZu9QZ543aRu5dOJsxZaV9NyhQfSYEd/9mmvDzvRy/CdiUJl4DLE4Tvwg83Fc/RESy5S8fnuDsxeL8Mlqq4hw5ERFJsSOUi758Y/GWi0AJQUSK3OWuOG8caOKeIi8XgRKCiBS5v354ktbzl7l/cfHOLkpSQhCRohbZF2PimDK+cqO+26SEICJF63JXnNcPNLLyllmMHVXc5SJQQhCRIva3jxLlomKfXZSkhCAiRSsSjTFhdCn/XeUiQAlBRIpUZ1ec12ubuPuW2SoXBUoIIlKU/vbRKU6du8TaWzW7KEkJQUSKUmR/jPGjS7nrplm57kreUEIQkaLT2RXn9f2N3H2zZhelUkIQkaKz++gpTp67pNlFPSghiEjR2R6NMW5UKV9VuagbJQQRKSpdcef12kS5aNxolYtSKSGISFHZffQUJ9pVLkpHCUFEikokGmPsqBK+erO+jNZTnwnBzH5hZs1mtj8lNt3M3jSzI+F5Wsq+J82szswOm9mqlPjtZhYN+54JS2kSltt8JcR3mdncIb5GEREgUS7asb+Rr940i/Gj+1wfrOhkcofwArC6R+wJYKe7LwB2hteY2UISS2AuCm2eNbNkke45YCOJdZYXpJxzA3Da3ecDTwNPDfRiRER68+7HpzjR3qFy0VX0mRDc/Y8k1jpO9QCwJWxvAR5Mib/s7h3ufhSoA5abWQUw2d3fcXcHXuzRJnmuV4GVybsHEZGhtCMaY0xZCXffrNlF6Qx0DGG2u8cAwnPyv24lUJ9yXEOIVYbtnvFubdy9EzgDzBhgv0RE0oqnlIsmjFG5KJ2hHlRO95e99xLvrc2VJzfbaGY1ZlbT0tIywC6KSDGq+eQ0zW0d3K/fLrqqgSaEplAGIjw3h3gDMCfluCrgeIhXpYl3a2NmZcAUrixRAeDum9292t2ry8s1Q0BEMheJxhhdVsLKW2bnuit5a6AJYRuwPmyvB15Lia8LM4fmkRg83h3KSm1mtiKMDzzSo03yXA8Bb4VxBhGRIZEoF8W468ZyJqpcdFV9/pcxs18BdwEzzawB+Ffgx8BWM9sAfAo8DODutWa2FTgAdAKPuXtXONWjJGYsjQN2hAfA88BLZlZH4s5g3ZBcmYhI8P6np2k628HaJZpd1Js+E4K7f/Mqu1Ze5fhNwKY08RpgcZr4RUJCERHJhu2hXKTZRb3TN5VFpKDF486OaCNfWVDOpLGjct2dvKaEICIF7e/1p2k8e5G1SzS7qC9KCCJS0Lbva2R0qWYXZUIJQUQKVnJ20VdunMlklYv6pIQgIgVrT0MrsTMXuX+xZhdlQglBRApWZF+MUaXGPQtVLsqEEoKIFCT3xG8XfXlBOVPGqVyUCSUEESlIexvOcKz1gn7quh+UEESkIEWiiXLRvZpdlDElBBEpOO7O9n0x7pw/kynjVS7KlBKCiBScfSoXDYgSgogUnEg0RlmJcZ9mF/WLEoKIFBR3J7I/US6aOn50rrszoighiEhB2X/sLPWnLrBW5aJ+U0IQkYKyPRqjtMS4V+WiflNCEJGC4e5EojH+8YYZTJugclF/KSGISMGoPX6WT0+dV7logAaVEMzsYzOLmtkeM6sJselm9qaZHQnP01KOf9LM6szssJmtSonfHs5TZ2bPhHWXRUT6JRLKRfct0toHAzEUdwhfdfdl7l4dXj8B7HT3BcDO8BozW0hiveRFwGrgWTMrDW2eAzYCC8Jj9RD0S0SKSLJcdMf1M5iuctGAZKNk9ACwJWxvAR5Mib/s7h3ufhSoA5abWQUw2d3fcXcHXkxpIyKSkQOxs3x88ry+jDYIg00IDrxhZu+Z2cYQm+3uMYDwnFzVuhKoT2nbEGKVYbtn/ApmttHMasyspqWlZZBdF5FCkiwXrVqk2UUDVTbI9ne6+3EzmwW8aWaHejk23biA9xK/Mui+GdgMUF1dnfYYESk+iXJRIyuun86MiWNy3Z0Ra1B3CO5+PDw3A78FlgNNoQxEeG4OhzcAc1KaVwHHQ7wqTVxEJCOHGts4euKcykWDNOCEYGYTzGxSchu4D9gPbAPWh8PWA6+F7W3AOjMbY2bzSAwe7w5lpTYzWxFmFz2S0kZEpE+RaIwSg1WaXTQogykZzQZ+G2aIlgH/6e6/N7N3ga1mtgH4FHgYwN1rzWwrcADoBB5z965wrkeBF4BxwI7wEBHpk7uzPRrjS/NmMFPlokEZcEJw94+ApWniJ4GVV2mzCdiUJl4DLB5oX0SkeB1uauOjlnN85855ue7KiKdvKovIiBaJNlJisFrlokFTQhCRES0SjbF83nTKJ6lcNFhKCCIyYn3Q1EZdc7tmFw0RJQQRGbG274thBqsXq1w0FAb7xTQRkWHV3tFJtOEM+xpaeeXdev5h7nRmTRqb624VBCUEEclbHZ1dHIq1sa+hlT31iSRQ19KOh98pqJo2jsdXLshtJwuIEoKI5IWuuPNRSzt76lvZF+4ADsbauNQVB2DGhNEsnTOVtUsqWFo1lSVVU/QzFUNMCUFEhp27c6z1AnvDX/176lvZf+wM5y4lvqs6YXQpt1ZN4Tt3zmXpnMSHf+XUcWiplOxSQhCRrDvZ3sG+hjPsbWhlb7gDOHnuEgCjS0u4pWIS/3R7FUuqprK0agrXl0+ktEQf/sNNCUFEhlR7Ryf7j5357IN/b0MrDacvAGAG88sn8tWbZ7G0agpLqqZyc8UkxpSV9nFWGQ5KCCIyYJc64xxqPMve+lb2hrr/kebug75Lq6byyB3XsaRqKosrpzBxjD528pX+ZUQkI/G489GJ9s9m++ytv3LQd0nVFNbcqkHfkUoJQUSukBz0Ta377z92lvaOTqD7oO+SqqksnaNB30KghCAinDp3qduA776GVk60dx/0/cZtlSydo0HfQqaEIFJkziUHfRs+r/vXn+o+6HvXTRr0LUZKCCIF7LNB34Yz7KtvZW9DK3XN7cR7DPp+60vXsXSOBn2LXd78y5vZauB/A6XAz939xznuksiI4e50xZ2PT577/MteDWc4ePysBn0lY3mREMysFPh34F6gAXjXzLa5+4Hc9kwKQTzudMadzniczrjT1eVcjsfpijudXWFfVzw8J47rijuXuxIfspfjcbq6Pm+fbNMVj39+TFc4X9jfFY9zOd593+UQ/7x9yr4e50v258r2Pfvx+b5UGvSVgciLhAAsB+rCspyY2cvAAyTWXx5SW9+tZ/OfPhpQW3fv+6De2ues8efNk9fw+evkfu/+usf7pV57n2177Oeq+3s5Z49+cpW2V+sLQDwOnfH4Z+WR4WQGo0pKKC0xykqNshKjrLSEshKjtMQYVRr2fba/5LPt0aNKKRtb9tmxyXapx/Q8X3K7YspYls2ZqkFfGZB8SQiVQH3K6wbgSz0PMrONwEaAa6+9dkBvNG3CaG6aPWlAbROdGHjTwTYf7F93ydbJ03z+2rq95rP9lvb47rEex/Q4yZXvlb5duv70+V499/e4DjNjVGmaD+ASo7S0hFE9PpwTx5Z8/oEbtktL7LMP9yvOFz7Mu+0rKaFEH8YyAuVLQkj3f88Vf9e5+2ZgM0B1dfWA/u67d+Fs7l04eyBNRUQKWr6smNYAzEl5XQUcz1FfRESKUr4khHeBBWY2z8xGA+uAbTnuk4hIUcmLkpG7d5rZ94DXSUw7/YW71+a4WyIiRSUvEgKAu0eASK77ISJSrPKlZCQiIjmmhCAiIoASgoiIBEoIIiICgA325xhyxcxagE8G2HwmcGIIuzMS6JqLg665OAzmmq9z9/J0O0ZsQhgMM6tx9+pc92M46ZqLg665OGTrmlUyEhERQAlBRESCYk0Im3PdgRzQNRcHXXNxyMo1F+UYgoiIXKlY7xBERKQHJQQREQGKLCGY2Rwz+4OZHTSzWjN7PNd9yjYzG2tmu81sb7jmH+W6T8PBzErN7O9m9rtc92W4mNnHZhY1sz1mVpPr/mSbmU01s1fN7FD4f/qOXPcpm8zspvBvm3ycNbPvD+l7FNMYgplVABXu/r6ZTQLeAx509yFfuzlfWGJdyQnu3m5mo4A/A4+7+99y3LWsMrP/CVQDk939a7nuz3Aws4+Bancvii9pmdkW4E/u/vOwjsp4d2/NcbeGhZmVAseAL7n7QL+ge4WiukNw95i7vx+224CDJNZzLlie0B5ejgqPgv4rwMyqgLXAz3PdF8kOM5sMfAV4HsDdLxVLMghWAh8OZTKAIksIqcxsLnAbsCvHXcm6UD7ZAzQDb7p7oV/zvwH/AsRz3I/h5sAbZvaemW3MdWey7HqgBfg/oTT4czObkOtODaN1wK+G+qRFmRDMbCLwa+D77n421/3JNnfvcvdlJNaqXm5mi3Pcpawxs68Bze7+Xq77kgN3uvsXgfuBx8zsK7nuUBaVAV8EnnP324BzwBO57dLwCOWxrwP/d6jPXXQJIdTRfw380t1/k+v+DKdwS/02sDq3PcmqO4Gvh3r6y8DdZvYfue3S8HD34+G5GfgtsDy3PcqqBqAh5W73VRIJohjcD7zv7k1DfeKiSghhgPV54KC7/zTX/RkOZlZuZlPD9jjgHuBQTjuVRe7+pLtXuftcErfVb7n7t3LcrawzswlhogShdHIfsD+3vcoed28E6s3sphBaCRTs5JAevkkWykWQR2sqD5M7gW8D0VBTB/hBWM+5UFUAW8KshBJgq7sXzVTMIjIb+G3ibx7KgP9099/ntktZ9z+AX4YSykfAd3Lcn6wzs/HAvcB3s3L+Ypp2KiIiV1dUJSMREbk6JQQREQGUEEREJFBCEBERQAlBREQCJQQREQGUEEREJFBCEBkiZjY3/C7/z8LaE2+Eb4eLjAhKCCJDawHw7+6+CGgF/im33RHJnBKCyNA66u57wvZ7wNzcdUWkf5QQRIZWR8p2F8X3e2EygikhiIgIoIQgIiKBfu1UREQA3SGIiEighCAiIoASgoiIBEoIIiICKCGIiEighCAiIoASgoiIBP8fZMFFOrYS/nUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_trees = [len(list(enumerate(tskit.all_trees(i, span=1)))) for i in range(2,8)]\n",
    "print(\"For a samples of 7 lineages the number of unique tree topologies is\", num_trees[-1])\n",
    "n_samp = np.arange(2,8)\n",
    "plt.xlabel('n')\n",
    "sns.lineplot(x = n_samp, y =  num_trees)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 13:</b>  Use the method <a href=\"https://tskit.dev/tskit/docs/stable/topological-analysis.html#sec-tree-ranks\"><code>rank</code></a> to find the shape of each of the 1,000 simulated trees. Set up an empty numpy array and use a for loop over <code>enumerate(replicates)</code> to add the shape integer of each simulated tree to this array (analogous to the loop at the start of this section).\n",
    "    \n",
    "<b>Hint: </b> Note that <code>replicates</code> is a generator object, which needs to be re-evaluated if you want to use the same set of simulations again.\n",
    "        </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL rank generates a tuple, the first element is an integer which uniquely records the unlabelled shape of a tree.\n",
    "num_replicates = 1_000\n",
    "rep_shapes = np.zeros(num_replicates)\n",
    "replicates = msprime.sim_ancestry(9, ploidy=1, num_replicates=num_replicates, random_seed=42)\n",
    "for replicate_index, ts in enumerate(replicates):\n",
    "    tree = ts.first()\n",
    "    rep_shapes[replicate_index] = tree.rank()[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> <a href=\"https://realpython.com/python-counter/\"><code>Counter</code></a> is extremely useful. The output is a python dictionary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 14:</b> Use a <code>Counter</code> to find out how many unique tree shapes you have simulated. How does this compare to the number possible bifurcating tree shapes, which for $n=9$ is 46 (Felsenstein, 2004)?\n",
    "        </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of unique tree shapes is 46\n"
     ]
    }
   ],
   "source": [
    "#SOL\n",
    "count_dict = Counter(rep_shapes)\n",
    "print(\"The number of unique tree shapes is\", len(count_dict))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Counter` returns a python dictionary. `most_common()` returns the key and value (in this case the count) of a `Counter`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "753.0\n",
      "91\n"
     ]
    }
   ],
   "source": [
    "most_freq_rank, count = count_dict.most_common()[0]\n",
    "print(most_freq_rank)\n",
    "print(count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">    \n",
    "    <b>Exercise 15:</b> If we focus on the root partition and consider a sample of $n=9$, how likely are asymmetric trees, i.e. with a (8,1) compared to very symmetric trees (5,4)? Check the frequency of all root partitions in the 1,000 replicates.\n",
    "\n",
    "Use the function <code>get_root_bipartition</code> defined below on every replicate and store the result in an empty array <code>rep_partitions</code>!\n",
    "        </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_root_bipartition(tree):\n",
    "    root_children = tree.children(tree.root)\n",
    "    left_clade_n = len(sorted(tree.leaves(root_children[0])))\n",
    "    right_clade_n = len(sorted(tree.leaves(root_children[1])))\n",
    "    return sorted((left_clade_n, right_clade_n))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_root_bipartition(simple_tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SOL\n",
    "num_replicates = 1_000\n",
    "rep_partitions = np.zeros(num_replicates)\n",
    "replicates = msprime.sim_ancestry(9, ploidy=1, num_replicates=num_replicates, random_seed=42)\n",
    "for replicate_index, ts in enumerate(replicates):\n",
    "    tree = ts.first()\n",
    "    rep_partitions[replicate_index] = get_root_bipartition(tree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try to use the code chunk below to count the number of root partitions and print the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 229 genealogies with a root partition 3.0\n",
      "There are 254 genealogies with a root partition 4.0\n",
      "There are 247 genealogies with a root partition 1.0\n",
      "There are 270 genealogies with a root partition 2.0\n"
     ]
    }
   ],
   "source": [
    "counts = Counter(rep_partitions)\n",
    "count_array = np.zeros(4)\n",
    "for index, root_partition in enumerate(counts):\n",
    "    count_array[index] =  counts[root_partition]\n",
    "    print(\"There are\", counts[root_partition], \"genealogies with a root partition\", root_partition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Python recipe:</b> Statistical tests with <code>scipy</code>. We can check whether counts are significantly different from each other with a $\\chi^2$ test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.4640000000000004\n"
     ]
    }
   ],
   "source": [
    "#By hand good for the soul\n",
    "chi_square = ((count_array - 250)**2)/250;\n",
    "print(chi_square.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Power_divergenceResult(statistic=3.4640000000000004, pvalue=0.3254613028395365)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#SOL scipy\n",
    "from scipy.stats import chisquare\n",
    "chisquare(count_array)"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}