notebooks/coalescent/r.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kingman coalescent and the Newick format\n",
    "\n",
    "Authors:\n",
    "- Alex Zarebski @aezarebski\n",
    "- Gerry Tonkin-Hill @gtonkinhill\n",
    "\n",
    "Date: 2018-10-03\n",
    "\n",
    "In this notebook we implement the Kingman coalescent and implement some functions for working with trees inspired by Newick format. Newick format is a widely used way to represent tree data structures. Having the genealogy in Newick format makes it easy to read into `ape` --- a popular package in R for working with genealogies --- and use the visualisation functionality it provides. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(ape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model and implementation\n",
    "\n",
    "Given we have `k` copies of the gene in a population of size `pop_size` the probability of at least one pair coming from the same parent is *approximately* `0.25 * k * (k - 1) / pop_size`. Using discrete generations would suggest a geometric number of generations until the first coalescence where the probability of coalesence in each generation is this value. We can approximate the geometric distribution with an exponential distribution with this rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "coalescent_rate <- function(k, pop_size) {\n",
    "    0.25 * k * (k - 1) / pop_size\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following functions, `leaf_node` and `branch_node` are helpers to work with trees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "leaf_node <- function(name, time) {\n",
    "    list(type = \"leaf\", name = name, time = time)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "branch_node <- function(name, children, time) {\n",
    "    list(type = \"branch\",\n",
    "         name = name,\n",
    "         children = children,\n",
    "         time = time,\n",
    "         lengths = c(time - children[[1]]$time,\n",
    "                     time - children[[2]]$time))\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Taking two nodes and linking them as the children of a parent is achieved with the following function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "binary_parent <- function(child1, child2, time) {\n",
    "    parent_name <- paste(child1$name, child2$name, sep = \"-\")\n",
    "    branch_node(parent_name, list(child1, child2), time)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start by setting up a little sample population in a larger population to work on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "current_time <- 0\n",
    "sampled_population <- list(leaf_node(\"beth\", current_time),\n",
    "                           leaf_node(\"gerry\", current_time),\n",
    "                           leaf_node(\"morty\", current_time),\n",
    "                           leaf_node(\"summer\", current_time))\n",
    "population_size <- 100\n",
    "\n",
    "k <- +Inf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Until the population has multiple individuals who have not coalesed continue to coalese individuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "while (k > 2) {\n",
    "    k <- length(sampled_population)\n",
    "    coalescent_time <- rexp(1, coalescent_rate(k, population_size))\n",
    "    current_time <- current_time + coalescent_time\n",
    "    ixs <- sample.int(k, size = 2)\n",
    "    parent_node <- binary_parent(sampled_population[[ixs[1]]], sampled_population[[ixs[2]]], current_time)\n",
    "    if (k > 2) {\n",
    "        sampled_population <- c(list(parent_node), sampled_population[-ixs])\n",
    "    } else {\n",
    "        sampled_population <- parent_node\n",
    "    }\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function recursively constructs a Newick representation of the tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "newick_helper <- function(node) {\n",
    "    if (node$type == \"leaf\") {\n",
    "        node$name\n",
    "    } else if (node$type == \"branch\") {\n",
    "        sprintf(\"(%s:%f,%s:%f)%s\",\n",
    "                newick_helper(node$children[[1]]),\n",
    "                node$lengths[1],\n",
    "                newick_helper(node$children[[2]]),\n",
    "                node$lengths[2],\n",
    "                node$name)\n",
    "    }\n",
    "}\n",
    "\n",
    "newick <- function(node) {\n",
    "    sprintf(\"%s;\", newick_helper(node))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "demo_tree <- read.tree(text = newick(sampled_population))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAIAAAByhViMAAAABmJLR0QA/wD/AP+gvaeTAAAg\nAElEQVR4nO3deZzVVf348XNnB2aGZdhBUBZxDbfU0dBwC5dfYYF7bomo37Lom6ZlUtZXRc1J\n6+vS15TEJXDJUL+C9CvNBS0zNzREEBBEZd9mYJZ7f3+MjoMPfq4MA2+fz7+8Z84593yuyn35\nuTNjJpfLJQAAtn55rX0AAAA2DWEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDC\nDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELY\nAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7\nAIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEH\nABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewA\nAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0A\nQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMA\nCELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAA\nQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAg\nCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAE\nIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAUDLWvvGYxeeMWLXAb3bFhcU\nlrTbbud9/+MXt9Zk359wy65dMpnMH5fWNF/1+y90yWQyExe/OzhlaO9MJjNt2fLfjBnZt3Np\nQVG7HSq/9sfXVqWU5k65Ydju25UUFnbstcPZVz7UfJNPtyqltGbeY+edfvSAXl1KCoo6dus3\n7PjvTZu9uvmEqYf2yWQyf1229LJRh3UtK97t3Kc22evFZ5EDAFpMzdIp/UoKOu54/EPPzFpZ\nXVe9YuHEy45PKe0y6oGmOfuUF+cX9ajLbrDwwA7F+UXda98bPG+bsrz80qrjDr763ulrauve\n+Ne9PYry23b7xoKpPz7k29fMfGf1+lUL/+uw3imln85c3rTJp1s1738vrSjM73HAqKnPzFxd\nu37uv6Z8tW9ZQXHvu99c0zTn/G3K8wo63HvWV6594LnqunXL1ze0zOvHJyPsAKAFPffzPVNK\n57yyrPngj7fvvf3uh9Xncrlcrnb1Myml8j4/aj6hrnpmJpMp7/OT9wbqexXnp5R++MibTXOq\n+ndIKe31laubRpa8eGxKaa/Lnv8sq2qWTOlelF/ac+SSuvdbbcVr41JKvQ78w3sDDX1KCvKL\nuo3+07xP/pLQgnwUCwAtqKhTUUrp7xMebfbRa/rFzDdmPjs1P6WU0ur5N6aUulQe1nzV2kU3\n5XK57kPfHaxZOnnh+oayXudcfmCPpjlv1WZTSmfdfGbTyPplq1JKhe0LP8uq+08b/VZtw7ET\nr60oeD8S2vf7fml+3tIXbm58uG7ZQ/PX1ReVV1731T6f4jWh5RS09gEA2ELV1tZefvnlgwYN\nystzF+BTKikpOezkm4/85T4PXnp0jz8deNKx3/j6iGP237Fb8znz7nompdT/1H7NBxdNeyKl\ntMPp7w6unDkhpbT92Wc2n3Pf0pri8v2+1bNd08ibDyxMKe12UPdPvSpXv+LsKQvyCjpcse8G\nh0yZgrL8zPJc3bs7zxqfUhpw6o/8k7GlEXYAbNwVV1wxduzY1j7FVm/atGmT//3aHddcddPv\nb7v64nOvvvjcfvsf898Tbhq2XVnjhH9Omp9SGj64U/NVMyfMTSmduFPHxodzbnklpbTvyPdv\nj61f+cjM6roug89qvurvkxeklM7sXfqpV1UvvnNpXUNKKyoKN9Js7Xse9O7O42eklPbZMEbZ\nEgg7ADZu4MCBKaUxY8ZUVla29lm2ViUlJUOHDs3Lzz/pB+NO+sG4d/795PjfXv2Ta+762uDX\n5iz5R6+ivJTSjW+szstve2LXts0X3vLv5flFPb5e0abx4V/+vCilNLpXadOEVa/flFLqe9zg\n5qtuWLC6uLxyt3aFn3pV3doXU0oDTnh01u0HfMh1PTL1zZTS6L7ln+jVYDMQdgBsXOMnsJWV\nlSNHjmztswTRdYf9zr96vy8sHHD4pGevXrD6l/3aN6x7/ZnVtW06H1men2maVrP4rsnL1pX1\nGVXw7ljDjYvWFrcfsut7xZZSmjvhhZTSF4/epmlk/cpHX1xb13mXUe8NfJpVeUVdU0r5bfI/\n9DqyNyxaW1y+z16lhR86jVbgw3EAaCk3nnvMTjufm91wsHZNfUppcGlhSilbvyKllF9Q0ezr\nDVd+47u5XK77lw9tfFyz5L4F6+vL+57RfJMnJi9IKY3u/cG7cX2O2f2zrGrX/YxOhXmL/nx/\n81X1Na8O7b/7eTe90vhw3dLJ89fVl/UZ/fFfBzYbYQcALWXPnqteefnXB/xH1YwFyxty2TVL\n5v/xN2OOmzp/wPBfndy1bUqpsHT3r3Zpu/adCb/5vy/XZuvffvWpHx9/4JphPVJK2x7ft3GT\nFTNvSyn1PX635jtfv3BNcfl+g5vfjbvt+ZTSXiP6fJZV+UW9H7j4kFXzxg0f+/u5S9dk62te\nfvzukXvsM31Z5uCjms5za0qp73F7bLKXiU1H2AFAS9nrgil/+d0lpf+8acgOPYsKinsM3Ofn\nk2ad9+s/zrj3u01zJky/Y3jlgAuO+EJpaZcjzrlm9wsnj1hcm1I6YKcOjRMafwZig89Plz/8\nak1dWd9RzZ/ryckLUkqjtin7LKtSSpUXTX1s/KXVU8YN7lNR1K7i0NN+UTrsh8/OfWpY97bN\nd957hF90siXK5HK51j4DAFuiu+6665hjjpk0aZLvsYOthTt2AABBCDsAgCCEHQBAEMIOACAI\nYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh\n7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCE\nHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISw\nAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2\nAABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIO\nACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgB\nAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsA\ngCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAARR0NoH+Hypra29/PLLBw0alJcnqYEt3fTp01NK\n2Wy2tQ8CfFzCbrO64oorxo4d29qnAPgEZsyY0dpHAD4uYbdZDRw4MKU0ZsyYysrK1j4LwEeY\nPn16VVXVzjvv3NoHAT4uYbdZNX4CW1lZOXLkyNY+C8BHq6qq8q0jsBXxrysAQBDCDgAgCGEH\nABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewA\nAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0A\nQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMA\nCELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAA\nQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAg\nCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAE\nIewAAIIQdgAAQQg7AIAghB0AQBDCDgAgCGEHABCEsAMACELYAQAEIewAAIIQdgAAQQg7AIAg\nhB0AQBDCDgAgCGEHABCEsAMACELYAUA0uWzNpCt/2tqnoBUIOwCI5h+XHDjq17WtfQpagbAD\ngFDq1j4//LJnOw0+pLUPQisQdgDQsqYM7Z3JZKYtW/6bMSP7di4tKGq3Q+XX/vjaqpTS3Ck3\nDNt9u5LCwo69djj7yoc+sHDRU5POHnlY/+6digsLyyt6HfB/Tp/4j3eaT5h6aJ9MJvPXZUsv\nG3VY17Li3c596vIBHYtKd1tU2zD3gYMzmcyBt702rn/HTCbz15XrP7D5BbtUZDKZsc8sbtFr\nZzMTdgDQsv4ye1VefumM/xhZd8B/vrxoxdy/37bq2QdP+tLpCx++aNSDtdc+/OKqZfPO22Xt\nDecf8bNXVzSteub6M/rud9zTbb5059+eX15d/er0STu+8+AJlf1/89Ky93eeuTKvoMOKH59Y\nOvzKectXPXLV3he8tvz3u3ZOKd29uDqXyz160oADh/VMKd26YE3zI7395A/GzVjWY8gvf7ZX\nl831MrA5CDsAaFENd7xTnW1Y89ZZE8YcvW+7woLeux19/jZl1W/fM/zqimm/Pnf7LqVFZT1H\n/3L/lNID985vXLPq9f8Z8u2bKyp/9sytF++9/TZtCwu7b7//tVNuzTasGTvy+vd2zv5hcXUm\nr3jq4b/9zpGD2xQUdyjKSynd/Mbq/MKK/1PRpnHSdqfsnFJ64d433j9Rrva7I27IKyi/6Y/n\nbLZXgc1D2AFAC6pZOnnh+oayXudcfmCPpsG3arMppbNuPrNpZP2yVSmlwvaFjQ9vP/4n67K5\nC+74XvP36eKOX9mhbeHK2VWND9cte2j+uvqi8srrvtqnaU59zay/raxt2/WEosy7I+0HnpZS\neuOeGU1z5kw6ceKitbt8574jKko26bXS+oQdALSglTMnpJS2P/vM5oP3La0pLt/vWz3bNY28\n+cDClNJuB3VPKWXr3rrgn4uLyvY+t2/ZB3bL5lK2YU2ucedZ41NKA079UfP38rVv3ZTL5Tru\nclTTSEnHYQPaFKyeP+HdHeqXnHTm/UWluz1w2QGb6hrZchS09gEA2EJls9mU0vTp01v7IFux\nkpKSjn96JaW078j3b6qtX/nIzOq6LoPPaj7z75MXpJTO7F2aUqpePHFVfbaiz9mZDXfL1i6c\nva6+qPyAxvE542eklPY5tV/zOYumPZ5S2vak5oOZs3qWnjfnz2+sb9imOP+Fq4dPX7X+mDvu\n3qY4f9NdKFsKYQfAxs2aNSulVFVVVVVV1dpn2Yqd2q1dSml0r9KmkVWv35RS6nvc4ObTbliw\nuri8crd2hSmlurUvpZTK+vdJG1q98MaGXK7rLqMbHz4y9c2U0ui+5c3nvHrb3JTSkCFdmw8O\nOap37pqXxr9dfWGXhcPHPl2+7em3H9d/k1wdWxphB8DGnX/++dlsdtCgQXl5vm/nUyopKTpn\n5Iji9kN2bVfYNDh3wgsppS8evU3TyPqVj764tq7zLqMaHxYU90kp1a344C8oeX7c3SmlIy8f\nmlJKKXvDorXF5fvsVVrYfM7EmSsymcyp3do2H9z2m7uma1569Pll+00dMX999qr/vargAzcD\niULYAbBxRUVFF198cWufYutWs+SeBevruww6o/ngE5MXpJRG9/7gPbw+x+ze+LBdz3O2Lbnk\nrZeurc8d3lRgtaue/ub4WeXbnXTd/t1TSuuWTp6/rr7zgNEbPF+u/r6l6wpL9xxQssH7e/v+\np6d057x77h5x+yvbHH7j93fsuKkvlC2F/wgDgJayYuZtKaW+x+/WfPD6hWuKy/cb3Pwe3m3P\np5T2GvHuZ695BRUP/GrkumVTDvjudbMWLqurXvHvJ+8bucdh75TuedeTNxZmGne+NaXU97g9\nmu+cy9XlpZRrWD1vdV3NyrVN48UdDhnUtnD2hPNWZypuu/PkFrpYtgTCDgBaypxbXkkf+NR1\n+cOv1tSV9R3VfNqTkxeklEZt8/7PwO48+o7n7rq6/T/+e6+B3dt06nXYt37R+RtjZ8x//LDu\nbZvvvPeIDb4PL5PX5r6LRpRlZw/sWrH/yY80/9LZvUpz2dxeP7p/SHnRpr1GtiiZXC7X2mf4\nHLnrrruOOeaYSZMmjRw5srXPAsDnRbb2zaHd+v2zzVfmzf9jRYF7OpH5HjsACG7CmQc9vqb4\nxqdvVXXh+RsMADGtq8uuWDSr6jtDT7t1zpg7nz1j+/atfSJanDt2ABBQw/r5vdv3W5HaDNpz\n6G+nzTrj4L6tfSI2B2EHAAHlF/dZsq6+tU/B5uajWACAIIQdAEAQwg4AIAhhBwAQhLADAAhC\n2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEI\nOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhh\nBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHs\nAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQd\nAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLAD\nAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYA\nAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4A\nIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEA\nBCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCA\nIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQ\nhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACC\nEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQ\nwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC\n2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEI\nOwCAIIQdAEAQwg4AIAhhBwAQhLADAAhC2AEABCHsAACCEHYAAEEIOwCAIIQdAEAQwg4AIAhh\nBwAQhLADAAhC2AEABCHsAKDFPXfPr4780q4dS0vadui25yEn/GnmynkPHJ7JZI6atqBpzpp5\nj513+tEDenUpKSjq2K3fsOO/N2326uabTD20TyaT+euypZeNOqxrWfFu5z610cFx/TtmMpm/\nrlz/gTNcsEtFJpMZ+8zizXC9tBZhBwAt6+lffmP3EWNq9j77H3MWr5j//FWnDjz1i0MmX/pc\nSulbe3RunDP/ocu2HTj09tldrpv8xJKaNc9Nvb54+s1H7LzTPYvWNu3zl5kr8wo6rPjxiaXD\nr5y3fNUjV+290cEDh/VMKd26YE3zM7z95A/GzVjWY8gvf7ZXl8135Wx+OTajSZMmpZQmTZrU\n2gcBYDOpXvynNvmZii/8pPngw6dsn1IqKNm2IZfL5XI1S6Z0L8ov7TlySV1D05wVr41LKfU6\n8A/vDTT0KSnIL+o2+k/zmu20kcG3nh6ZUtrjkn+9Pyu7/tge7fIKyh9cUrOpr48tizt2ANCC\nnvrP82oackf/7tvNBwefu0NKqV2P0Y1vw/efNvqt2oZjJ15bUfD++3L7ft8vzc9b+sLNjQ/X\nLXto/rr6ovLK677ap2nORgfbDzwtpfTGPTOaRuZMOnHiorW7fOe+IypKWuAS2YIIOwBoQZdM\nnp/JK/7ZFzo3H8xkUkqp1+FfTinl6lecPWVBXkGHK/bttsHKTEFZfibl6hofrZw1PqU04NQf\nNX/n3uhgScdhA9oUrJ4/ofFhtn7JSWfeX1S62wOXHbAJr4stk7ADgJZSt+afj6xY16biaz2L\nNnjDfefRt1NKO528XUqpevGdS+sasvUrKgrzMhtaVNtQ3OGgxiVzxs9IKe1zar/m+2x0MKXM\nWT1L16/48xvrG1JKL1w9fPqq9cN/e/c2xfktdqFsKQpa+wAAsBWora29/PLLBw0alJf3Ce6J\n5Nbcn1Iq7nDwB8afuHl2Sum0ge1TSnVrX0wpDTjh0Vm3f9gdtUemvplSGt23/CMHU0pDjuqd\nu+al8W9XX9hl4fCxT5dve/rtx/X/+Mdm6yXsAOCjXXHFFWPHjv10a/OLNvhB1Lq1z533yrKC\nkn7DOpWklPKKuqaU8tt8+O207A2L1haX77NXaeFHDaaU0rbf3DVd89Kjzy/bb+qI+euzV/3v\nVQWZT3d2tjLCDgA+2sCBA1NKY8aMqays/Pir8rOzRhx/Uc3SB1M6umnwltNHrKrPduhzVuOt\nv3bdz+hU+PNFf74/pf2b5tTXvHroLsfudeEdV56xY0pp3dLJ89fVdx4wuvnmGx1s1L7/6Snd\nOe+eu0fc/so2h9/4/R07fqKLZesl7ADgozV+AltZWTly5MhPtLDqd9d/b9rvzrj+mKtPP2j1\nvBcnXPWdf+3QP6XZvY5494PX/KLeD1x8yH4/GTd87I6/Ovcbfdrn//upB388atT0Zdv98Ki+\njXNWzLw1pdT3uD2a77zRwUbFHQ4Z1Lbw1Qnn5eV3nnznyZ/8ctla+eEJAGhB33nwmUvP+trD\nF32jU1nZAcf+IO/Q6y7ILE4pDT5lu6Y5lRdNfWz8pdVTxg3uU1HUruLQ035ROuyHz859alj3\nto0T5tzySkpp7xF9mu+80cEmZ/cqzWVze/3o/iHlRS10aWyB3LEDgBaUV9jtwuvvu/D690f+\n64LXM3ltfrxTp+bTvnTKhQ+fcuH/b5P9/ueV3P98rMFG2do3711c3a7HVx+86Iuf9uBsldyx\nA4CWcunxXxnYb3i22cjaN++5ZO7KnkOv3altC95bmXDmQY+vKf7VI7c2/43HfB74+w0ALeWg\nAWtee/1Ph5z/uzdW1DTUrn3hkT98be9T8jrsNemeb7bE062ry65YNKvqO0NPu3XOmDufPWP7\n9i3xLGzJhB0AtJR9f/63e67+z7qpl+3Yo7ywTafDTr2k4ugLn3n98f3aF2/y52pYP793WVHn\n7fa46dmy306bddUIv7ju88j32AFAy8n/+pirvj7mqs3xTMV9lqyr3wxPxJbMHTsAgCCEHQBA\nEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAI\nQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABB\nCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAI\nYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh\n7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCE\nHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISw\nAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2\nAABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAARR0NoH+HzJZrMppenTp7f2QQD4ZBr/\n6G78Yxy2WMJus5o1a1ZKqaqqqqqqqrXPAsAnNmPGjNY+AnwYYbdZnX/++dlsdtCgQXl5PgQH\n2JpMnz69qqpq5513bu2DwIcRdptVUVHRxRdf3NqnAODTqKqq8p/lbOH8AwoAEISwAwAIQtgB\nAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsA\ngCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcA\nEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAA\nghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBA\nEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAI\nQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABB\nCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgK47NLIAAAgoSURB\nVCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAaHHf712el1cwd33DJt/5\nnF5leXmFC2uzm3xntkbCDgBaVi5bM/7ttSUdh21bnP/Zt5p05U/ff9iw5vZ3qks6HdmryBs6\nKQk7AGhpNYvvWF6fLe9/ymff6h+XHDjq17VND6vfuX1VfbZ9/5M++87EIOwAoGUtfe7ulNJ2\n39z5M+5Tt/b54Zc922nwIe/v/K97U0rbnrTTZ9yZMIQdALSsOb+blVLab1j5rZecvUe/boUF\nhZ16Djr5wt/WZHPNp62Z99h5px89oFeXkoKijt36DTv+e9Nmr2766uUDOhaV7raotmHuAwdn\nMpkDb3stpTTn5tdSSl86vPOd/3XOHgN6FBcUlHfdbsS519T4jrvPK2EHAC3r/unvZDL5q749\nYl7f4fc8MbN6xYKrj28/4fLRQ85/omnO/Icu23bg0Ntnd7lu8hNLatY8N/X64uk3H7HzTvcs\nWts44YLXlv9+184ppbsXV+dyuUdPGpBSmvz04kxeYf13D3qy7cH3//P1lcvm/fTLmXt+/b0j\nrn+lVa6UVpfJ5XIfPQsAPt8mTpx43HHHjRkzprKy8uOvKikpOWLY0K7t2i+ry1706KKfH9C9\ncTxb905pm+4NJbuuX/N8Smnd0qnb9TxyTeevz533h4qCd++5rJx9RYcBP+x14B8WPHJs48iX\nO5Y8vra0ev2SokzjLusqStotq8t+87ZXbz1xYOOc1W9cVt7nRxU73rbk5RM3zZWzVSlo7QMA\nwFZg1qxZKaWqqqqqqqpPtPCBSectq8u273dBU9WllPIKu+7UtvCfq19Yl00leen+00a/Vdvw\nrYnXNlVdSql9v++X5l+49IWbUzo2pVRfM+tvK2tLe57wbtWlVL144rK6bNuuJ4x/r+qaFJSW\nfaqrZKsn7ADgo51//vnZbHbQoEF5eZ/gu5hKSkp2KboupbTrRad94Etv1Tbk5ZeV5KVc/Yqz\npyzIK+hwxb7dNpiRKSjLzyzP1TU+WvvWTblcruMuRzV9fdlLE1NKg87+QfMDrX3jqZRS7+F9\nP8G1EYiwA4CPVlRUdPHFF3+KhY8eNyal9LVDezYfrK+esXB9Q9vOh6WUqhffubSuIaUVFYUb\nScb2PQ9q/ItF0x5PKW17Ur+mL825eVZKaf8Ttm0+//VbZ6WU9jx6m09xVAIQdgDQgh584p2U\n0tAOxc0HV79xY0qp0xdOSSnVrX0xpTTghEdn3X7Ah+zz6m1zU0pDhnRtGpn6+DsppdG9N/jU\n9a9T30wpnbmNj2I/p/xULAC0mNz6W96uTinlZzYYfnHcwymlQy75Ykopr6hrSim/zUf8Tykm\nzlyRyWRO7db2vZ3rfvfW2uL2B+zStvk9moYbF60tKtt7z9LCTXYJbFWEHQC0lOrFE5fUNaSU\nfjtjedNgfc2sb/1hdkmnQ6/dt1tKqV33MzoV5i368/3NF9bXvDq0/+7n3fTeby3J1d+3dF1h\n6Z4DSt7NuJol975d21C+7RnNV61bOnn++vqyPme25DWxRRN2ANBSlr98V0pp14vOvPXIY+/7\n++y6bP1br04/99Ahc+rLr5x2R1l+JqWUX9T7gYsPWTVv3PCxv5+7dE22vublx+8eucc+05dl\nDj7q3Z+ByOXq8lLKNayet7quZuXalNKKf9+RUup7/OANnu7fE1JKfUbusZkvky2HsAOAlvLa\nDa+klM4+5+rHrj5w3CkHlhaVDNjn6y92PPKeZ1779h6dm6ZVXjT1sfGXVk8ZN7hPRVG7ikNP\n+0XpsB8+O/epYd3f/eA1k9fmvotGlGVnD+xasf/Jj6SU5oyfmVL64vANfkhizi0vp5T2GNFn\nc10fWxy/oBgAIAh37AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAA\nghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBA\nEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAI\nQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABB\nCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAI\nYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh\n7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCE\nHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISw\nAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2\nAABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIO\nACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgB\nAAQh7AAAghB2AABBCDsAgCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAghB2AABBCDsA\ngCCEHQBAEMIOACAIYQcAEISwAwAIQtgBAAQh7AAAgvh/NP+UVreeXP4AAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(demo_tree)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}