{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulate the Tree\n",
    "This notebook simulates a phylogenetic tree under the dual-birth model [(Moshiri & Mirarab, 2017)](https://doi.org/10.1093/sysbio/syx088) and performs some basic analyses of the tree."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**Rule 2: Document the Process, Not Just the Results.** Here we describe the steps how to produce the dataset.\n",
    "\n",
    "**Rule 7: Build a Pipeline.** This notebook describes the entire workflow, and its modularity makes it easy to change models or model parameters.\n",
    "\n",
    "**Rule 8: Share and Explain Your Data.** To enable reproducibility we provide a `/intermediate_data` directory with files produced by the workflow.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Tree Simulation Parameters\n",
    "In the dual-birth model, all branches are labeled as either *active* or *inactive* and is parameterized by two rates: the *birth rate* $\\left(\\lambda_b\\right)$, which is the Poisson rate at which active branches split, and the *activation rate* $\\left(\\lambda_a\\right)$, which is the Poisson rate at which inactive branches split, where $\\lambda_a\\le\\lambda_b$. We will start by choosing our desired values for these two rates. For our purposes, we will choose the rate estimates for *Alu* elements as found in the original manuscript [(Moshiri & Mirarab, 2017)](https://doi.org/10.1093/sysbio/syx088).\n",
    "\n",
    "Our tree simulations also require an end criterion: either a tree height or a number of leaves. We will specify that we want trees with *n* = 100 leaves.\n",
    "\n",
    "We will also include import statements here to keep the notebook clean and organized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "OUTPUT_TREE_FILE = \"./intermediate_data/dualbirth.tre\"\n",
    "BIRTH_RATE = 122.03\n",
    "ACTIVATION_RATE = 0.73\n",
    "NUM_LEAVES = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from seaborn import distplot\n",
    "from treesap import dualbirth_tree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate the Tree\n",
    "We can now use [TreeSAP](https://github.com/niemasd/TreeSAP) to simulate the tree, which will return a [TreeSwift](https://github.com/niemasd/TreeSwift) `Tree` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = dualbirth_tree(Lb=BIRTH_RATE, La=ACTIVATION_RATE, end_num_leaves=NUM_LEAVES)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A common file format to represent tree structures is the Newick format. In TreeSwift, `Tree` objects have a `newick()` function to output a Newick string representation of the tree. The `str()` cast automatically calls the `newick()` function, so we can use Python's `print()` function to print the tree in the Newick format conveniently. Thus, we could run `print(tree)` to print a Newick string representation of the tree we have simulated under the dual-birth model, but because of the length, we will only print a prefix to give you an idea of how it looks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[&R] (((((((((((((((((((((((((((((((((99:0.0039061301439481944,98:0.0039061301439481944):0.006016644...\n"
     ]
    }
   ],
   "source": [
    "print('%s...' % tree.newick()[:100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore the Tree\n",
    "Now that we've simulated our tree, we can begin to explore its properties and statistics. For example, consider the following questions:\n",
    "* How many lineages exist at any given \"time\" (where \"time\" is measured in unit of \"expected number of per-site mutations\")?\n",
    "* What is the branch length distribution? What about just internal branches? What about just terminal branches?\n",
    "* What is the height? Diameter?\n",
    "* What is the average branch length? Internal branch length? Terminal branch length?\n",
    "* What is the Colless balance index? Sackin balance index? Gamma statistic? Treeness?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xm4XFWZ7/Hvj4RJEBNIwECAMKpBNLGPiHjVCHoFlMEJcUSkxbnhOoEjtEM3tnIRmvuIQQVsB0C8KCrYIBrwqiCBxAh4kTC1gUCOQoSABkje/mOvgk1Zp84+p2rX3lX1+zxPPbVrj2/tqnNWvXvttZYiAjMzs05sUHUAZmbW/1yYmJlZx1yYmJlZx1yYmJlZx1yYmJlZx1yYmJlZx1yYWGkkvVDSTVXH0QuSFkn6x6rjAJB0tqTPdmlfl0g6ohv7ssHmwsQ6Jul2SS9tnh8Rv4iIp1URUzelQnFNejwoKXKv10jaoeoYJ6PI+4qIAyLinKpjtfqbWnUAZnUXEb8ANgeQNAe4DZgWEY821pFUeH+Spua3rUqR92VWlDMTK42kBZJW5F7fLulDkpZJ+ouk8yRtklv+SklLJa2W9CtJz8otO17SLZIekHSjpFfllk2RdLKkP0m6TdL70q/sqWn5UyR9TdJKSXdK+qykKWnZrpKuSPH8SdJ5HbzlHSX9MsV4qaQZ6RhzUjxHSfov4Gdp/sGSbkjvd5GkZ+TeU0jaNff6CZeuJH0kvZ+7JP1j8/rAdEk/TrFcLWmXybyh/OU7SW9L7++UFPOtkvZJ8/8oaVX+kpikjSV9UdJ/SbpH0hmSNp1MHFZ/Lkys1w4D9gd2Ap4FvA1A0nzg68A7ga2ArwAXSdo4bXcL8ELgKcA/A9+UNCstewdwADAPeA5waNMxzwYeBXYF5gP/E2jUb3wGuBSYDswG/r2D9/ZG4Ehga2Aj4ENNy18MPAN4uaTdge8AxwIzgYuBH0raaLyDSNof+ADw0vSeFrRY7XCy8zQdWA58buJvp6XnAcvIPqNvA+cCz01xvBk4XdLmad2TgN3JPpddge2AT3UpDqsZFybWa6dFxF0RcS/wQ7J/NABHA1+JiKsjYl26Tr8W2BsgIr6btlsfEecBNwN7pW0PA06NiBURcR/ZPzEAJG0DHAgcGxEPRsQq4BSyf7YAjwA7AttGxN8i4v918N7Oiog/RMRfgfNz763hxBTDX4HXAz+OiMsi4hHgi8CmwD4FjnNYOtYNEfEQcGKLdS6MiN+kS1bfahHLZN0WEWdFxDrgPGB74NMRsTYiLgUeBnZVdt3vaOB/RcS9EfEA8C88ft5twLgwsV67Ozf9EOmaPdk/9A+myyerJa0m+0e1LYCkt+Yuga0GngnMSNtuC/wxt9/89I7AhsDK3LZfIcseAD4CCPhNuuT09hLeW6u4tgXuaLyIiPVp+XYFjtPu/RaNZbLuyU3/FSAimudtTpZtPQm4Nnfef5Lm2wByBbzVxR+Bz0XE312OkbQjcCawH/DriFgnaSlZIQCwkuwSVcP2TftdC8xoVbEcEXeTXSZD0v8AfirpyohY3oX39HeHy03fBezZeJF+yW8P3JlmPUT2z7jhqUCj/qnd+62LP5EVLHtExJ3jrWz9z5mJdcuGkjbJPSb6Q+VM4F2SnqfMZpJeIenJwGZk/4hHASQdSZaZNJwPHCNpO0nTgOMaCyJiJVmdyMmStpC0gaRdJL047et1khr/mO9Lx1k/8bc/YecDr5C0n6QNgQ+SFXq/SsuXAm9MNxfsT1bfkt/2SEnPkPQk4JM9iHdCUqZ1JnCKpK0B0ufz8mojs7K4MLFuuZjsl2jjceJENo6IxWQZwulk/9SXkyrnI+JG4GTg12SXWfYEfpnb/EyyAmMZsCTF8iiwLi1/K1mF+I1p3xcAjcr75wJXS1oDXAQcExG3TiT2yYiIm8gqrP+d7Ff8QcBBEfFwWuWYNG818Cbg+7ltLwFOA35Odp6uSovWlh33BB1Hik/S/cBPgb5vd2StyYNj2aCRdABwRkTsWHUsvZBuKb4e2NhtRKwqzkys70naVNKBkqZK2g44Abiw6rjKJOlVqR3HdODzwA9dkFiVXJjYIBBZm4r7yC5z/Z7Bb8/wTmAVWfubdcC7qw3Hhp0vc5mZWcecmZiZWcf6up3JjBkzYs6cOVWHYWbWV6699to/RURXG5D2dWEyZ84cFi9eXHUYZmZ9RdId4681Mb7MZWZmHXNhYmZmHXNhYmZmHXNhYmZmHXNhYmZmHXNhYmZmHXNhYmZmHSutMJH0dUmrJF2fm7elpMsk3Zyep6f5CyT9JY2kt1TSoPerZGY2UMrMTM4G9m+adzxweUTsBlyeXjf8IiLmpcenS4zLzMy6rLTCJCKuBO5tmn0IcE6aPgc4tKzjm5lZ7/S6zmSbNIwqwN3ANrllz5f0W0mXSNpjrB1IOlrSYkmLR0dHSw3WzMyKqawCPrK+7xv9318H7BgRzyYbxvT7bbZbGBEjETEyc2ZX+ykzM7NJ6nVhco+kWQDpeRVARNwfEWvS9MXAhpJm9Dg2MzObpF4XJhcBR6TpI4AfAEh6qiSl6b1SXH/ucWxmZjZJpXVBL+k7wAJghqQVZONynwScL+ko4A7gsLT6a4F3S3oU+CtweHgISDOzvlFaYRIRbxhj0X4t1j0dOL2sWMzMrFxuAW9mZh1zYWJmZh1zYWJmZh1zYWJmNkTmz59fyn5dmJiZDZHbbrutlP2WdjeXmZlVb/78+YyOjrLrrrsCsGbNmlKO48zEzGyAjY6OPqEAmTJlSinHcWZiZjaAGhnJmjVr2HzzzVm0aBEACxYs4Iorruj68ZyZmJkNoHxB0otOcZ2ZmJn1ueZ6EeCxgmTFihU9icGZiZlZn2uuFwF6lpE0ODMxMxsA+XqRKjgzMTOzjrkwMTPrU/Pnz2f27NmltR2ZCBcmZmZ9qtd3bLXjOhMzsz7T3IakV3dstePMxMysz9QpI2lwZmJm1gfybUnqlJE0ODMxM+sD+bYkdcpIGpyZmJnVXCMrmTlzZqVtSdqpJDORdIyk6yXdIOnYNG9LSZdJujk9T68iNjOzuhkdHQWoXTaS1/PCRNIzgXcAewHPBl4paVfgeODyiNgNuDy9NjMzsoJkyZIlVYcxpioyk2cAV0fEQxHxKHAF8GrgEOCctM45wKEVxGZmVhuNRomjo6OsXr266nDaqqIwuR54oaStJD0JOBDYHtgmIlamde4Gtmm1saSjJS2WtLiR+pmZDaJGpfvMmTPZbbfdqg6nrZ5XwEfE7yV9HrgUeBBYCqxrWickxRjbLwQWAoyMjLRcx8ys3+Ur3et0C/BYKqmAj4ivRcQ/RMSLgPuAPwD3SJoFkJ5XVRGbmVkd9EOle15Vd3NtnZ53IKsv+TZwEXBEWuUI4AdVxGZmVqV85411r3TPq6qdyfckbQU8Arw3IlZLOgk4X9JRwB3AYRXFZmZWmTp2lVJEJYVJRLywxbw/A/tVEI6ZWeXq2HnjRLg7FTOzGujXjKTB3amYmVWs3+7casWZiZlZxfrtzq1WXJiYmVUon5X0y51brbgwMTOr0CBkJeA6EzOznsoPcgX0XXuSsTgzMTProfwgV1DPga4mw5mJmVkPNLcjqesgV5PlzMTMrAf6vR3JeJyZmJmVqN9bthflzMTMrESDnpE0ODMxMyvJILRsL8qZiZlZSQalDUkRE8pMJG0AbB4R95cUj5lZX8u3IxmUNiRFjJuZSPq2pC0kbUY2fvuNkj5cfmhmZv0n345k0OtJ8opkJnMj4n5JbwIuAY4HrgW+UGpkZmY119yaHRjYdiTjKVJnsqGkDYFDgYsi4hEgyg3LzKz+mluzw3BlI3lFMpOvALcDvwWulLQj4DoTMxsYrTKMIoY1C2ll3MIkIk4DTsvNukPSS8oLycyst1plGEUMaxbSyriFiaRtgH8Bto2IAyTNBZ4PfK3s4MzMyjTo/WX1UpE6k7OB/wS2Ta//ABxbVkBmZr0yLK3Te6FIncmMiDhf0kcBIuJRSesme0BJTwPOy83aGfgUMA14BzCa5n8sIi6e7HHMzIoY5P6yeqlIYfKgpK1Id3BJ2hv4y2QPGBE3AfPSvqYAdwIXAkcCp0TEFye7bzMzq0aRwuQDwEXALpJ+CcwEXtul4+8H3BIRd0jq0i7NzNprriuxzhW5m+s6SS8GngYIuCm1NemGw4Hv5F6/T9JbgcXAByPivuYNJB0NHA2www47dCkMMxsmrivpviLdqbwaOJisMNkdOEjSfpK27uTAkjZK+/1umvVlYBeyS2ArgZNbbRcRCyNiJCJG/CUws8lq1JUMQ79ZvVDkMtdRZLcC/zy9XkDWncpOkj4dEf8xyWMfAFwXEfcANJ4BJJ0J/GiS+zUzsx4rUphMBZ7R+Gef2p18A3gecCUw2cLkDeQucUmaFREr08tXkXUqaWbWFc29+bqupLuKFCbb57MGYFWad6+kSdWdpB6IXwa8Mzf73yTNI7tr7PamZWZmHRnW3nx7pUhhskjSj3i8buM1ad5mwOrJHDQiHgS2apr3lsnsy8xsLK2yEbdyL0eRwuS9ZAXIC9LrbwDfi4gA3EeXmdWWs5HeKXJrcAAXpIeZWa05G6lGkVuD95Z0jaQ1kh6WtE6Su6A3s1pyNlKNIpe5TidrXPhdYAR4K1l7EzOzWnI20ntFeg0mIpYDUyJiXUScBexfblhmZtZPimQmD6XW6ksl/RtZ6/RChZCZmQ2HIoXCW9J67wMeBLYnu7vLzMwMKHY31x2SNgVmRcQ/9yAmM7MJc0/A1SpyN9dBwFLgJ+n1PEkXlR2YmdlEuCfgahWpMzkR2AtYBBARSyXtVGJMZmaFtGpT4lETq1GkzuSRiGgeWTHKCMbMbCLcpqQ+imQmN0h6IzBF0m7APwG/KjcsM7Ni3KakHopkJu8H9gDWknUZfz9wbJlBmZkVsXr1alavnlR/s9ZlRe7megj4eHqYmdXGtGnTqg7BknELE0m7Ax8C5uTXj4h9ywvLzOyJ8pXtDb4NuD6K1Jl8FzgD+CqwrtxwzMxay1e2N7jSvT6KFCaPRsSXS4/EzGwMjaxk5syZrmyvqSIV8D+U9B5JsyRt2XiUHpmZWTI6OgrgLKTGimQmR6TnD+fmBbBz98MxM3tcvouUmTNnsmTJkqpDsjEUuZvLrd3NrBLuIqV/jFmYSNo3In4m6dWtlkfE/53sQSVNI6vQfyZZlvN24CbgPLK7xm4HDouI+yZ7DDPrL+3u1nIXKfXXrs7kxen5oBaPV3Z43FOBn0TE04FnA78Hjgcuj4jdgMvTazMbEr5bq7+NmZlExAnp+cjmZZImPZ6JpKcALwLelvb/MPCwpEOABWm1c8g6ljxusscxs3przkQaWYjv1upPkx0x8ZQOjrkTMAqcJWmJpK9K2gzYJiJWpnXuBrZptbGkoyUtlrS4cYeHmfWf5kzEWUh/K3I3Vyvq8JjPAd4fEVdLOpWmS1oREZJa9kwcEQuBhQAjIyPuvdis5lrVhYAzkUEz2cykk3/iK4AVEXF1en0BWeFyj6RZAOl5VQfHMLOaaFUXAs5EBk27u7l+R+tCQ4xxCaqIiLhb0h8lPS0ibgL2A25MjyOAk9LzDyZ7DDOrXvMwus5ABlu7y1yd3rHVzvuBb0naCLgVOJIsSzpf0lHAHcBhJR7fzErmNiLDpd3dXHeUddCIWAqMtFi0X1nHNLNyjXV3ltuIDIfJ1pmYmT2B784abpO9m8vMhpzbiVjemJmJpMvT8+d7F46Z9QtnIpbXLjOZJWkf4GBJ59LUtiQiris1MjOrFWci1k67wuRTwCeB2cD/bloWgIftNRsizkSsnXZ3c10AXCDpkxHxmR7GZGY1kc9GnIlYO0XGM/mMpIPJOmcEWBQRPyo3LDOrg3w24kzE2hm3MJH0r8BewLfSrGMk7RMRHys1MjOrlMddt4kocmvwK4B5EbEeQNI5wBLAhYnZAPO46zYRRRstTstNP6WMQMysfjzuuhVVJDP5V2CJpJ+T3R78IjwKopmZ5RSpgP+OpEXAc9Os4yLi7lKjMrPKNPf2a1ZEoe5U0giIF5Uci5nVgHv7tclw31xmQ869/Vo3uNdgsyHnlu3WDW0zE0lTgBsi4uk9isfMesQjIVo3tc1MImIdcJOkHXoUj5n1iOtGrJuK1JlMB26Q9BvgwcbMiDi4tKjMrFT51u2uG7FuKFKYfLL0KMysp9y63bqtSDuTKyTtCOwWET+V9CRgSvmhmVm35etJ3Lrdumncu7kkvQO4APhKmrUd8P1ODyxpiqQlkn6UXp8t6TZJS9NjXqfHMLMncj2JlaXIZa73kvUafDVARNwsaesuHPsY4PfAFrl5H07jqJjZBDW3F2nFbUisLEXamayNiIcbLyRNJRtpcdIkzSbrjfirnezHzB7X3F6kFWckVpYimckVkj4GbCrpZcB7gB92eNwvAR8Bntw0/3OSPgVcDhwfEWs7PI7ZUHF7EatKkczkeGAU+B3wTuBi4BOTPaCkVwKrIuLapkUfBZ5O1qHklsBxY2x/tKTFkhY37kgxM7NqFbmba30aEOtqsstbN0VEJ5e5XgAcLOlAYBNgC0nfjIg3p+VrJZ0FfGiMeBYCCwFGRkY6utxmNijy7UbMqlDkbq5XALcApwGnA8slHTDZA0bERyNidkTMAQ4HfhYRb5Y0Kx1PwKHA9ZM9htmwcbsRq1qROpOTgZdExHIASbsAPwYu6XIs35I0k2wArqXAu7q8f7OBlM9K3G7EqlKkMHmgUZAktwIPdOPgEbEIWJSm9+3GPs2GjbMSq4MxCxNJr06TiyVdDJxPVmfyOuCaHsRmZk1atSVxa3arg3aZyUG56XuAF6fpUWDT0iIys5bmz5/PsmXLmDp16hMKE7cdsToYszCJiCN7GYiZtTc6OsrUqVOZO3eu25JY7YxbZyJpJ+D9wJz8+u6C3qy7xusOxZezrM6KVMB/H/gaWav39eWGYza8xusOxZezrM6KFCZ/i4jTSo/EzNwdivWtIoXJqZJOAC4FHusrKyKuKy0qMzPrK0UKkz2BtwD78vhlrkivzawL3B2K9bsihcnrgJ3z3dCbWXe54aH1uyK9Bl8PTCs7ELNh5zu1rJ8VyUymAf9f0jU8sc7EtwabmRlQrDA5ofQozIaY60tsEBQZz+SKXgRiNqxcX2KDoEgL+Ad4fMz3jYANgQcjYosyAzMbJO1at7tluw2CIpnJY+O0p4GrDgH2LjMos0HTrnW7W7bbIChSZ/KYNFzv91MjxuPLCclscDQykjVr1rh1uw20Ipe5Xp17uQEwAvyttIjMBki+IHH2YYOsSGaSH9fkUeB2sktdZlbA5ptvzooVK6oOw6xURepMPK6JmZm11W7Y3k+12S4i4jMlxGM2MNx+xIZJu8zkwRbzNgOOArYCXJiYteH2IzZM2g3be3JjWtKTgWOAI4FzgZPH2m48kjYBrgQ2Tse/ICJOSCM6nktWUF0LvMWdS1q/ymclbj9iw6BtR4+StpT0WWAZ2T/+50TEcRGxqoNjrgX2jYhnA/OA/SXtDXweOCUidgXuI8uAzPqSsxIbNmMWJpK+AFwDPADsGREnRsR9nR4wMo3WWxumR2N8lAvS/HOAQzs9llmvzZ8/n9mzZ7tVuw2ddpnJB4FtgU8Ad0m6Pz0ekHR/JweVNEXSUmAVcBlwC7A6Ih5Nq6wAthtj26MlLZa0uPHrz6wu3K7EhlW7OpMiY51MSkSsA+ZJmgZcCDx9AtsuBBYCjIyMxDirm/VMvp7E7Ups2JRWYBQREauBnwPPB6ZJahRus4E7KwvMbBJcT2LDrOeFiaSZKSNB0qbAy4DfkxUqr02rHQH8oNexmU2G60nMJtjRY5fMAs6RNIWsMDs/In4k6Ubg3HT32BLgaxXEZjZhricxq6AwiYhlwPwW828F9up1PGaTkR+fpFGQuJ7EhlmldSZm/So/PokzErNqLnOZ9aVW2YjHJzHLODMxK8jZiNnYnJmYFZBvQ+JsxOzvOTMxK8BtSMzac2ZiNobmOhK3ITEbmzMTszG4jsSsOGcmZm34ji2zYpyZmJlZx5yZ2FDL14s0a7QlMbPxOTOxoZavF2nmehKz4pyZ2MBql3U0uCW7WXc4M7GB1S7raHD2YdYdzkxs4DQyEmcdZr3jzMQGjscXMes9ZyY2UDwOu1k1nJnYQHEfWmbVcGZifa35ji33oWVWDWcm1rfmz5/PsmXLHstGwHdnmVXFmYn1rdHRUaZOncrcuXN9x5ZZxSopTCR9HXglsCoinpnmnQi8A2j8zPxYRFxcRXxWjSKNDPN8ScusPqq6zHU2sH+L+adExLz0cEEyZIo0MszzJS2z+qgkM4mIKyXNqeLYVq6JZhd5bmRo1r/qVgH/PknLJH1d0vRWK0g6WtJiSYvzFa9WDxPNLvKcaZj1rzpVwH8Z+AwQ6flk4O3NK0XEQmAhwMjISPQyQGsv32DQ2YXZcKlNZhIR90TEuohYD5wJ7FV1TDYxbjBoNrxqk5lImhURK9PLVwHXVxmPtTbeYFK+u8psOFV1a/B3gAXADEkrgBOABZLmkV3muh14ZxWxWXseTMrMWqnqbq43tJj9tZ4HYpPiO67MrFlt6kzMzKx/1abOxOplrLqRRlsQM7M8ZybW0lh1I64XMbNWnJn0qU5amhfh1uhmNhHOTPpUJy3Ni3AGYmYT4cykT7QaBMqZg5nVhTOTPtGciThzMLM6cWZSY/lsxJmImdWZM5May2cjzkTMrM6cmVRsvL6unI2YWT9wZlIx93VlZoPAmUkFXBdiZoPGmUkFXBdiZoPGmUmJxuvfytmImQ0KZyYlcv9WZjYsnJmUxOOhm9kwcWZSEo+HbmbDxJlJF7SqG/F46GY2TJyZdEGruhHXi5jZMBmIzKTssT3G47uzzGzYDURmUvbYHuNxFmJmw652mYmk/YFTgSnAVyPipLHWXbZsGQsWLHBmYGZWsVplJpKmAP8HOACYC7xB0tyx1l+/fj3gzMDMrGp1y0z2ApZHxK0Aks4FDgFubLXypptu6mzEzKwG6laYbAf8Mfd6BfC8/AqSjgaOTi/XSrq+R7EVNQP4U9VBtFDHuBxTMXWMCeoZl2Mq5mnd3mHdCpNxRcRCYCGApMURMVJxSE9Qx5ignnE5pmLqGBPUMy7HVIykxd3eZ63qTIA7ge1zr2eneWZmVmN1K0yuAXaTtJOkjYDDgYsqjsnMzMZRq8tcEfGopPcB/0l2a/DXI+KGNpss7E1kE1LHmKCecTmmYuoYE9QzLsdUTNdjUkR0e59mZjZk6naZy8zM+pALEzMz61itChNJ+0u6SdJySce3WL6xpPPS8qslzckt+2iaf5OklxfdZ1kxSXqZpGsl/S4975vbZlHa59L02LpHMc2R9Nfccc/IbfMPKdblkk6TpB7F9KZcPEslrZc0rxvnqWBcL5J0naRHJb22adkRkm5OjyNy88s+Vy1jkjRP0q8l3SBpmaTX55adLem23Lma14uY0rJ1ueNelJu/U/qsl6fPfqNexCTpJU3fqb9JOjQtK/s8fUDSjenzuVzSjrllpXyfOomrq9+piKjFg6zC/RZgZ2Aj4LfA3KZ13gOckaYPB85L03PT+hsDO6X9TCmyzxJjmg9sm6afCdyZ22YRMFLBeZoDXD/Gfn8D7A0IuAQ4oBcxNa2zJ3BLN87TBOKaAzwL+Abw2tz8LYFb0/P0ND29R+dqrJh2B3ZL09sCK4Fp6fXZ+XV7dZ7SsjVj7Pd84PA0fQbw7l7F1PQ53gs8qUfn6SW5Y72bx//2Svk+dSGurn2n6pSZPNaVSkQ8DDS6Usk7BDgnTV8A7JdK8UOAcyNibUTcBixP+yuyz1JiioglEXFXmn8DsKmkjSdw7K7HNNYOJc0CtoiIqyL7Fn0DOLSCmN6Qtu2WceOKiNsjYhmwvmnblwOXRcS9EXEfcBmwfy/O1VgxRcQfIuLmNH0XsAroRqd0nZynltJnuy/ZZw3ZZ9+T89TktcAlEfHQBI7dSUw/zx3rKrK2clDe96mjuLr5napTYdKqK5XtxlonIh4F/gJs1WbbIvssK6a81wDXRcTa3LyzUur4yQmmtZ3GtJOkJZKukPTC3PorxtlnmTE1vB74TtO8yZ6nonFNdNtenKtxSdqL7FfoLbnZn0uXKk6Z4A+XTmPaRNJiSVc1LieRfbar02c9mX125TyRZcHN36lenaejyDKNdtt2+n3qNK7HdPqdqlNhMpAk7QF8HnhnbvabImJP4IXp8ZYehbMS2CEi5gMfAL4taYseHbstSc8DHoqIfF9rVZ2nWku/Zv8DODIiGr/KPwo8HXgu2aWU43oY0o6RdRfyRuBLknbp4bHHlM7TnmTt1hp6cp4kvRkYAb5Qxv4na6y4uvGdqlNhUqQrlcfWkTQVeArw5zbbdto9SycxIWk2cCHw1oh4rLSPiDvT8wPAt8nS1NJjSpcB/5yOfS3ZL5Dd0/qzc9v39Dwlf/cLssPzVDSuiW7bi3M1plT4/xj4eERc1ZgfESsjsxY4i+5/p8aU+5xuJavnmk/22U5Ln/WE99lpTMlhwIUR8Ugu1tLPk6SXAh8HDs5djSjr+9RpXN37ThWtXCn7QdYa/1ayCvRGJdIeTeu8lydW4p6fpvfgiRXwt5JVSo27zxJjmpbWf3WLfc5I0xuSXVN+V49imglMSdM7py/cltG6EvDAXsSUXm+QYtm5W+epaFy5dc/m7yvgbyOrLJ2epntyrtrEtBFwOXBsi3VnpWcBXwJO6lFM04GN0/QM4GZS5S/wXZ5YAf+eXsSUm38V8JJenieygvQWUqV22d+nLsTVte9U4YB78QAOBP6Q3vTH07xPk5WkAJukL+jy9AHk//l8PG13E7m7IVrtsxcxAZ8AHgSW5h5bA5sB1wLLyCrmTyX9g+9BTK9Jx1wKXAcclNvnCHB92ufppN4RevTZLQCuatpfx+epYFzPJbvG/CDZr+kbctu+PcW7nCz979W5ahkT8Gbgkabv1Ly07GfA71Jc3wQ271FM+6Tj/jY9H5Xb587ps16ePvuNe/jZzSH7gbJB0z7LPk8/Be7JfT4Xlf196iSubn6n3J2KmZl1rE51JmZm1qda2yzSAAABM0lEQVRcmJiZWcdcmJiZWcdcmJiZWcdcmJiZWcdqNdKiWV1J2orsfnyApwLrgNH0+qGI2KeSwMxqwrcGm02QpBPJesr9YtWxmNWFL3OZdUjSmvS8IHWg+QNJt0o6Sdl4Lb9J41XsktabKel7kq5JjxdU+w7MOufCxKy7ng28C3gGWceUu0fEXsBXgfendU4FTomI55L1SvDVKgI16ybXmZh11zURsRJA0i3ApWn+78gGKAJ4KTA316P+FpI2j4g1PY3UrItcmJh1V37MmvW51+t5/O9tA2DviPhbLwMzK5Mvc5n13qU8fsmLiY5DblZHLkzMeu+fgJE0gt2NZHUsZn3NtwabmVnHnJmYmVnHXJiYmVnHXJiYmVnHXJiYmVnHXJiYmVnHXJiYmVnHXJiYmVnH/hvjufv6IXtcmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tree.ltt(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that, as we move forward in time, the number of lineages grows more and more rapidly. This is because, in the dual-birth model, each lineage is a Poisson process, so the rate at which splitting events occurs increases as more lineages are created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "branch_lengths = dict()\n",
    "branch_lengths['all'] = list(tree.branch_lengths())\n",
    "branch_lengths['internal'] = list(tree.branch_lengths(terminal=False))\n",
    "branch_lengths['terminal'] = list(tree.branch_lengths(internal=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VPW5+PHPkz2EhECC7BBUEFEBEVHcCopVW/dad6XVys+61dvbW7W1dalttd5rtdVeL2or2ha3qli1VkHBlV0WcQFU9jUBAmTP5Pn98T0Dk5CZOTOZSTLJ8369zmvOnDnLMwOZZ77L+X5FVTHGGNN5pbV1AMYYY9qWJQJjjOnkLBEYY0wnZ4nAGGM6OUsExhjTyVkiMMaYTs4SgWlTIqIicnBbxxELEblTRP6a4HP+S0QmJehcJ4rIFyHPV4vIxESc2zvfchEZn6jzmbZniaAT8r4YqkRkj4jsEJHXRGRAW8cVTTK+gH1cc7yIrG/hOVREKrzPu0xEZorIRaH7qOoZqjrV57kiJk5VfU9VD2lJzCHXe1JE7mly/sNUdVYizm/aB0sEnddZqtoV6ANsAf4YbkcRSW+1qDqukd7nfQjwJPCwiNyR6IuISEaiz2k6PksEnZyqVgMvAMOD27xfgf8rIq+LSAUwQUS+LSIfi8guEVknIneG7F/i/VKdJCJrRaRURH4e8nq6iPxMRL4Ukd0isrBJCWSiiKwUkZ0i8oiISKzvQ0T6isg/RGSbiHwtIjeFvHaniDwnIk95118uImNCXh/tvbfdIvK8iDwrIveISB7wL6Cv92t+j4j09Q7LCne+SFS1VFWfBn4I3CYiRV4Ms0TkB976wSIyW0TKvc/yWW/7u95plnixXBQssYjILSKyGfhLmFLM0SLyqVcC/IuI5Hjn/J6IvN/ks1QvhsnAZcBPvev903t9b1WTiGSLyIMistFbHhSRbO+1YGz/KSJbRWSTiHw/5Drf8mLaLSIbROQnfj5Dk3iWCDo5EekCXATMafLSpcCvgXzgfaACuBIoBL4N/FBEzm1yzAm4X7ynAL8UkUO97T8GLgG+BRQAVwGVIcedCRwNjAAuBE6L8T2kAf8ElgD9vOvfLCKh5zkbeMaL/xXgYe/YLOAl3K/0HsA04DwAVa0AzgA2qmpXb9kY6XwxmA5kAGObee1XwJtAd6A/XmlNVU/yXh/pxfKs97y3F/sgYHKY612G+1wPAoYCt0cLUFWnAH8Dfudd76xmdvs5cCwwChjpvZ/Qc/cGuuH+Xa4GHhGR7t5rTwD/T1XzgcOBt6PFZJLDEkHn9bKI7ATKgVOB+5u8Pl1VP1DVBlWtVtVZqrrMe74U94X5jSbH3KWqVaq6BPelPNLb/gPgdlX9Qp0lqloWcty9qrpTVdcC7+C+VGJxNNBTVe9W1VpV/Qp4DLg4ZJ/3VfV1VQ0AT4fEdizuC/kPqlqnqi8C83xcM9z5fFHVOqAU9wXeVB3uS72v99m/38w+oRqAO1S1RlWrwuzzsKquU9XtuAR/SSzxRnAZcLeqblXVbcBdwBUhr9d5r9ep6uvAHtyPheBrw0WkQFV3qOqiBMVkYmSJoPM6V1ULgRzgBmC2iPQOeX1d6M4icoyIvONVvZQD1wLFTc65OWS9EujqrQ8AvowQS7jj/BqEq77ZGVyAnwG9Ilwjx6tP7wts0MajLzZ67z5jDp7PFxHJBHoC25t5+aeAAPO8aqeropxum1fFF0noe1qDe9+J0Nc7X7hzl6lqfcjz0H/f7+BKiWu8qrBxCYrJxMgSQSenqgHvV3AAV7Wz96Umu/4dVwUyQFW7AY/ivqz8WIerkkiWdcDXqloYsuSr6rd8HLsJ6NekXSK0/SJZw/OeA9TTTOlDVTer6jWq2hf4f8CfJHJPIT8xhr6ngUCwiqsC6BJ8ocmPAT/n3ohLxM2dOyJVna+q5wAHAC8Dz/k5ziSeJYJOTpxzcPXRn0XYNR/YrqrVIjIW14bg1+PAr0RkiHe9EcFG0jikiUhOyJKN+zLd7TWY5oprnD5cRI72cb6PcEnwBhHJ8D6L0Hr7LUCRiHSLM95GRKSHiFwGPALc16SKLLjPd0Wkv/d0B+7LuCEkngPjuPT1ItJfRHrg6vWD7QtLgMNEZJTXgHxnk+OiXW8acLuI9BSRYuCXQNQuviKSJSKXiUg3r5psF/veo2lllgg6r3+KyB7cH+CvgUmqujzC/tcBd4vIbtwfeyy/3h7w9n/Tu94TQG5cUbu67aqQ5Uuvnv5MXNvC17i698dxjZQRqWotcD6uIXMncDnwKlDjvf457svuK6/aKd4qlSXe570K12byH6r6yzD7Hg3M9fZ/BfiR1+4B7ot6qhfLhTFc/++4z/8rXDXdPQCqugK4G5gBrMR1DAj1BK4ef6eIvNzMee8BFgBLgWXAouC5fbgCWC0iu3BVjZfF8H5MAolNTGNMYyIyF3hUVf/S1rEY0xqsRGA6PRH5hoj09qqGJuG6sb7R1nEZ01rsLkRjXHfG54A8XNXJBaq6qW1DMqb1WNWQMcZ0clY1ZIwxnVxKVA0VFxdrSUlJW4dhjDEpZeHChaWq2jPafimRCEpKSliwYEFbh2GMMSlFRNZE38uqhowxptOzRGCMMZ2cJQJjjOnkUqKNwBjT8dXV1bF+/Xqqq6MNpGqaysnJoX///mRmZsZ1vCUCY0y7sH79evLz8ykpKUFin6Su01JVysrKWL9+PYMHD47rHFY1ZIxpF6qrqykqKrIkECMRoaioqEUlKUsExph2w5JAfFr6uVki6IhU3WKMMT4kLRGIyCEisjhk2SUiN3sTc7wlIiu9x+7Rz2Zict558M1vtnUUxqSkl19+GRHh888/B2D16tUcfvjhAMyaNYszzzyzLcNLiqQlAm+i8lGqOgo4CjdX6UvArcBMVR0CzPSem0R5912YPh1mzGjrSIxJSdOmTeOEE05g2rRpbR1Kq2mtqqFTcDNJrcHN1TrV2z4VOLeVYugcFi7ct15R0XZxGJOC9uzZw/vvv88TTzzBM88809bhtJrW6j56MW66P4BeIWO9bwZ6NXeAiEwGJgMMHDgw6QF2GDt27FvftAkOjjTnuTHt0803w+LFiT3nqFHw4IOR95k+fTqnn346Q4cOpaioiIULF1JUFO/02qkj6SUCEckCzgaeb/qauskQmm3VVNUpqjpGVcf07Bl18DwTtHPnvvVNNreKMbGYNm0aF198MQAXX3xxp6keao0SwRnAIlXd4j3fIiJ9VHWTiPQBtrZCDJ1HaIlg48a2i8OYFoj2yz0Ztm/fzttvv82yZcsQEQKBACLC9ddf3/rBtLLWaCO4hH3VQgCvAJO89UnA9FaIofPYsQMGDXLrViIwxrcXXniBK664gjVr1rB69WrWrVvH4MGDWbduXVuHlnRJTQQikgecCrwYsvle4FQRWQlM9J6bRNm5Ew46CLKzrURgTAymTZvGeeed12jbd77zHX7729+2UUStJyXmLB4zZozaxDQ+HXYYHHoozJ8P48fD1KlRDzGmPfjss8849NBD2zqMlNXc5yciC1V1TLRj7c7ijmbnTigshG7dYPfuto7GGJMCLBF0NDt2QPfuUFAA5eVtHY0xJgVYIuhIamqgqsqVCAoKYNeuto7IGJMCLBF0JMF7CIIlAksExhgfLBF0JMFEYCUCY0wMLBF0JJWV7jEvzxKBMcY3SwQdSVWVe8zJcYmgshLq69s2JmNSSNeuXaPu8+CDD1IZ/NGVRKHDXyebJYKOJJgIcnNdIgDrQmpMgsWTCAKBQJKiSQxLBB1JaCLo1s2tW/WQMTGbNWsW48eP54ILLmDYsGFcdtllqCp/+MMf2LhxIxMmTGDChAkAvPnmm4wbN47Ro0fz3e9+lz179gBQUlLCLbfcwujRo3n++ecZP348t9xyC2PHjmXo0KG89957gPvlf+KJJzJ69GhGjx7Nhx9+2Orvt7WGoTatobkSgSUCk4raahzqEB9//DHLly+nb9++HH/88XzwwQfcdNNNPPDAA7zzzjsUFxdTWlrKPffcw4wZM8jLy+O+++7jgQce4Je//CUARUVFLFq0CIBHH32U+vp65s2bx+uvv85dd93FjBkzOOCAA3jrrbfIyclh5cqVXHLJJbT2SAqWCDoSSwTGJMzYsWPp378/AKNGjWL16tWccMIJjfaZM2cOn376KccffzwAtbW1jBs3bu/rF110UaP9zz//fACOOuooVq9eDUBdXR033HADixcvJj09nRUrViTrLYVliaAjqa52j5YITKpri3Gom8jOzt67np6eTn0zHS9UlVNPPTXsvAV5eXnNnjP0fL///e/p1asXS5YsoaGhgZycnES9Bd+sjaAjsRKBMUmXn5/Pbq8TxrHHHssHH3zAqlWrAKioqIj5F315eTl9+vQhLS2Np59+uk0ali0RdCSh3Ufz8926JQJjEmry5MmcfvrpTJgwgZ49e/Lkk09yySWXMGLECMaNG8fnn38e0/muu+46pk6dysiRI/n888/3K0W0BhuGuiO54w64+25oaHB3GffoAb//vWt4M6ads2GoW8aGoTZOVZUrDYi4u4sBKiraNiZjTLtniaAjqapy7QMAWVmQmQlen2ZjjAnHEkFHEpoIALp2tRKBMSYqSwQdSdNEkJdnJQJjTFSWCDoSKxEYY+KQ1EQgIoUi8oKIfC4in4nIOBHpISJvichK77F7MmPoVKqrXWNxUF6eJQJjTFTJLhE8BLyhqsOAkcBnwK3ATFUdAsz0nptEsKohY+JSVlbGqFGjGDVqFL1796Zfv357n9fW1ibsOuvWrdtv2IlYnHDCCSxO9BhMJHGICRHpBpwEfA9AVWuBWhE5Bxjv7TYVmAXckqw4OpWqqn2jjoKrGtq6te3iMSZFFBUV7f2CvfPOO+natSs/+clPfB8fCARIT0+Put+AAQN49tln444zWZJZIhgMbAP+IiIfi8jjIpIH9FLVTd4+m4FezR0sIpNFZIGILNi2bVsSw+xArERgTMJNnTqVsWPHMmrUKK677joaGhqor6+nsLCQm2++mREjRjBv3jz69+/Pz372M0aOHMnRRx/NokWL+OY3v8lBBx3EY489BsCqVasYNWoUAI8//jgXXHABp512GkOGDOG2227be83JkyczZswYDjvsMO6+++6kv8dkDjqXAYwGblTVuSLyEE2qgVRVRaTZW5tVdQowBdydxUmMs+OwxmLTQdz8xs0s3pzYKpBRvUfx4OmxDWb3ySef8NJLL/Hhhx+SkZHB5MmTeeaZZ7jwwgspLy/npJNO4sGQAfIGDx7MkiVLuPHGG7n66qt5//332bNnDyNHjuSaa67Z7/xLlixh4cKFZGZmMnToUG688Ub69u3LvffeS48ePaivr2fChAlccMEFDB8+vMWfQTjJTATrgfWqOtd7/gIuEWwRkT6quklE+gBWd5EoViIwJqFmzJjB/PnzGTPGjdJQVVXFgAEDAMjKyuK8885rtP/ZZ58NwBFHHEF9fT15eXnk5eWRlpa2d8KaUBMnTqTAGyBy2LBhrF27lr59+zJt2jSeeOIJ6uvr2bhxI59++mlqJgJV3Swi60TkEFX9AjgF+NRbJgH3eo/TkxVDp9O015CVCEyKivWXe7KoKldddRW/+tWvGm2vr68nNzcXEWm0PTjMdFpaWqNhrNPS0podxrq5oa5XrlzJQw89xLx58ygsLOTyyy+nOjjEfJIku9fQjcDfRGQpMAr4DS4BnCoiK4GJ3nOTCDU1EPIfi7w8qK2Furq2i8mYFDZx4kSee+45SktLAde7aO3atUm95q5du8jPz6egoIBNmzbx73//O6nXgyRPTKOqi4HmRr47JZnX7bRqa90YQ0GhA88VFrZNTMaksCOOOII77riDiRMn0tDQQGZmJo8++ih9+/ZN2jVHjx7N8OHDGTZsGIMGDdo7+1ky2TDUHYUqpKfD7be7oagBHnsMJk+G9euhX7+2jc+YKGwY6paxYagNBAIuGTRXIrAGY2NMBJYIOorg3Y+hiaBrV/doDcbGmAiiJgIRGSoiM0XkE+/5CBG5PfmhmZg0lwisRGBSTCpUVbdHLf3c/JQIHgNuA+q8Cy4FLm7RVU3i1dS4RysRmBSVk5NDWVmZJYMYqSplZWXkhHYdj5GfXkNdVHVek/6y+3eINW3LSgQmxfXv35/169djQ8rELicnh/79+8d9vJ9EUCoiBwEKICIXAJsiH2JaXTARhN5HYCUCk0IyMzMZPHhwW4fRKflJBNfjxvwZJiIbgK+By5IalYldpBKBJQJjTAR+EoGq6kRv5NA0Vd0tIpa22xurGjLGxMlPY/E/AFS1QlV3e9teSF5IJi7NJYLcXBCxEoExJqKwJQIRGQYcBnQTkfNDXioA4m+eNsnRXK8hERuB1BgTVaSqoUOAM4FC4KyQ7buB/QfWNm2ruRIB2AikxpiowiYCVZ0OTBeRcar6USvGZOLRXK8hsBKBMSYqP43FH4vI9bhqor1VQqp6VdKiMrGzEoExJk5+GoufBnoDpwGzgf646iHTnoRLBFYiMMZE4ScRHKyqvwAqVHUq8G3gmOSGZWJmJQJjTJz8JILg9FY7ReRwoBtwQPJCMnFprtcQuBKBJQJjTAR+2gimiEh34BfAK0BX4JdJjcrEzqqGjDFxipoIVPVxb3U2cGBywzFxC9dryKqGjDFRRE0EIlIIXAmUhO6vqjclLywTMysRGGPi5Kdq6HVgDrAMaIjl5CKyGtfDKADUq+oYEekBPItLLKuBC1V1RyznNc2I1FhcWQkNDZBmE9IZY/bnJxHkqOqPW3CNCapaGvL8VmCmqt4rIrd6z29pwfkN7EsEmZmNtwcHnqus3DcstTHGhPB1H4GIXCMifUSkR3BpwTXPAaZ661OBc1twLhNUUwMZGfv/6rc5CYwxUfhJBLXA/cBHwEJvWeDz/Aq8KSILRWSyt62XqgYnttkM9GruQBGZLCILRGSBzVjkQ23t/tVCYENRG2Oi8lM19J+4m8pKo+65vxNUdYOIHAC8JSKfh76oqioizU5QqqpTcBPiMGbMGJvENJra2v17DIFNTmOMicpPiWAVUBnPyVV1g/e4FXgJGAtsEZE+AN7j1njObZoIVyKwqiFjTBR+SgQVwGIReQeoCW6M1n20yYxmecA3gbtxN6VNAu71HqfHGbsJZVVDxpg4+UkEL3tLrHoBL4lI8Dp/V9U3RGQ+8JyIXA2sAS6M49ymKSsRGGPi5OfO4qnR9glz3FfAyGa2lwGnxHNOE0FNjZUIjDFxiTRV5XOqeqGILMP1/mlEVUckNTITGysRGGPiFKlE8CPv8czWCMS0ULReQ1YiMMaEEbbXUEhf/+tUdU3oAlzXOuEZ36I1FluJwBgThp/uo6c2s+2MRAdiWihcIkhPh5wcKxEYY8KK1EbwQ9wv/4NEZGnIS/nAB8kOzMSotjb8WEI2OY0xJoJIbQR/B/4F/BY3MFzQblXdntSoTOzC9RoCm5PAGBNRpDaCclVdDdwObPbaBgYDl3tzFJj2JFzVENicBMaYiPy0EfwDCIjIwbixfwbgSgumPQnXawisRGCMichPImhQ1XrgfOCPqvpfQJ/khmViZiUCY0yc/CSCOhG5BDdd5avetswI+5u2ECkRWInAGBOBn0TwfWAc8GtV/VpEBgNPJzcsEzMrERhj4hSp+2iBqu5S1U+BvSONeslgWqtEZ/yzXkPGmDhFKhHMCq6IyMwmr8UzGqlJpmglAksExpgwIiUCCVlvOkexYNqPhgaorw/fayhYNaQ20ZsxZn+REoGGWW/uuWlLdXXuMVLVUCDgSg3GGNNEpDuLDxCRH+N+/QfX8Z73THpkxr/gF3ykqiFwpYJwpQZjTKcVqUTwGG5coa4h68Hnjyc/NONbtERgcxIYYyIIWyJQ1btaMxDTAjXeVNJ+SgTGGNOEn/sITHtnJQJjTAtYIugIgokgUq8hsBKBMaZZUROBiKS3RiCmBaxEYIxpAT8lgpUicr+IDI/nAiKSLiIfi8ir3vPBIjJXRFaJyLMiEubby/jmt9eQJQJjTDP8JIKRwArgcRGZIyKTRaQghmv8CPgs5Pl9wO9V9WBgB3B1DOcyzYml+6gxxjQRNRGo6m5VfUxVjwNuAe4ANonIVG+OgrBEpD/wbbzupiIiwMnAC94uU4FzWxC/gei9hqxqyBgTga82AhE5W0ReAh4E/gc4EPgn8HqUwx8Efgo0eM+LgJ3e/AYA64F+Ya47WUQWiMiCbdu2RX8nnZmVCIwxLRDpzuKglcA7wP2q+mHI9hdE5KRwB4nImcBWVV0oIuNjDUxVp+BmRGPMmDE2pEUk0XoNZWVBRoaVCIwxzfKTCK5U1fdDN4jI8ar6gareFO4g4HjgbBH5FpADFAAPAYUikuGVCvoDG+KM3QRFKxGI2JwExpiw/DQW/6GZbX+MdpCq3qaq/VW1BLgYeFtVL8OVLi7wdpsETPcZqwknWiIAm5PAGBNWpIlpxgHHAT1DBpwD98u+JfcW3AI8IyL3AB8DT7TgXAb8JQIrERhjwohUNZSFG2AuAzfYXNAu9v2i90VVZ+FNdKOqXwFjYzneRBGt1xDY5DTGmLAiDTo3G5gtIk+q6ppWjMnEyqqGjDEtEKlq6EFVvRl4WET267WjqmcnNTLjX7ReQ+BKBGVlrROPMSalRKoaetp7/O/WCMS0gN8Swdq1rROPMSalRKoaWug9zg5uE5HuwABVXdoKsRm/rLHYGNMCfu4sniUiBSLSA1gEPCYiDyQ/NONbba27VyA9QmcuayMwxoTh5z6Cbqq6CzgfeEpVjwEmJjcsE5OaGlcaEAm/j5UIjDFh+EkEGSLSB7gQeDXJ8Zh41NZGrhYCVyKoqYFAoHViMsakDD+J4G7g38AqVZ0vIgfixh8y7UVtbeQeQ2BzEhhjwoo61pCqPg88H/L8K+A7yQzKxMhPiSB0BNKCWKaTMMZ0dFETgYj0BK4BSkL3V9WrkheWiYnfqiGwdgJjzH78jD46HXgPmAFYBXN75CcR5HujhOzenfx4jDEpxU8i6KKqtyQ9EhO/YK+hSAoL3ePOncmPxxiTUvw0Fr/qzSlg2is/JYJu3dxjeXny4zHGpBQ/ieBHuGRQLSK7RGS3iOxKdmAmBn56DVmJwBgThp9eQ/nR9jFtzEoExpgW8DPEhIjI5SLyC+/5ABGx+QTaEz+JoKDA3XlsJQJjTBN+qob+BIwDLvWe7wEeSVpEJnZ+GovT0lwysBKBMaYJP72GjlHV0SLyMYCq7hCRKN86plX5KRGAqx6yEoExpgk/JYI6EUkHFPbeYNaQ1KhMbPwmgsJCSwTGmP34SQR/AF4CDhCRXwPvA79JalQmNn56DYErEVjVkDGmCT+9hv4mIguBUwABzlXVz6IdJyI5wLtAtnedF1T1DhEZDDwDFAELgStUtbYF78HEUiJYty758RhjUoqfXkNHAEcAW4FZfpKApwY4WVVHAqOA00XkWOA+4PeqejCwA7g6rsjNPrG0EViJwBjTRKTJ67vhxhkaACzFlQaOEJG1wDneZDVhqariehgBZHqLAiezrwfSVOBO4H/jfwvGV68hcCWCHTuSH48xJqVEKhH8ClgADFHV81T1XGAoMB/4tZ+Ti0i6iCzGlSbeAr4EdqpqvbfLeqBfmGMni8gCEVmwbds2f++ms/JbIigudo3F9fXR9zXGdBqREsFE4FZV3dtDSFUDwM/wOVWlqgZUdRTQHxgLDPMbmKpOUdUxqjqmZ8+efg/rfFRjSwQA27cnNyZjTEqJlAhqQ3657+Vtq4nlIqq6E3gHd2NaoYgEq6T6AxtiOZdpIhBwycBPr6GiIvdYWprcmIwxKSVSr6EcETkS1zYQSnA9gSLy7jeoU9WdIpILnIprKH4HuADXc2gSrh3CxKvW63AVS4mgrCx58RhjUk6kRLAJeCDMa5t9nLsPMNW7GS0NeE5VXxWRT4FnROQe4GPgiVgCNk3EkwisRGCMCRE2EajqhJacWFWXAkc2s/0rXHuBSYQar5bOEoExJk5+7iw27VksJQJrIzDGNMMSQaqLJRHk5kJeniUCY0wjlghSXTAR+Ok1BK5UYI3FxpgQke4sHh3pQFVdlPhwTMyCiSAz09/+xcVWIjDGNBKp19D/RHgtOFSEaWvBxmK/JQJLBMaYJpLWa8i0kngSwZdfJi8eY0zK8TP6aBcRuV1EpnjPh4jImckPzfhSXe0ec3L87W8lAmNME34ai/8C1ALHec83APckLSITm3hKBOXlUFeXvJiMMSnFTyI4SFV/B9QBqGol+w87YdpKrIkgeC+BDTxnjPH4SQS13lhBwTmLDyLGQedMEsVTIgCrHjLG7BV1qkrgDuANYICI/A04HvheMoMyMQgmgljaCMASgTFmLz9zFr8lIouAY3FVQj9SVfsWaS+CjcVWIjDGxMlPiQAgBze/cAYwXERQ1XeTF5bxLd6qIbu72BjjiZoIROQ+4CJgORCcrUwBSwTtQbyJYLOfkcSNMZ2BnxLBucAhqmoNxO1RrIkgK8slg02bkheTMSal+Ok19BXgcyAb0+pqaiAjA9LT/R/Tp48lAmPMXn5KBJXAYhGZSUi3UVW9KWlRGf+qq/2XBoL69rVEYIzZy08ieMVbTHtUUxN7IujTB5YvT048xpiUEzERePMNf1NVL2uleEys4k0EmzdDQwOk2ZQUxnR2Eb8FVDUADBIRH9NfmTZRU+P/ZrKgvn2hvt7uJTDGAP6qhr4CPhCRV4CK4EZVfSDSQSIyAHgK6IXrbjpFVR8SkR7As0AJsBq4UFV3xBW9ib9EAK6d4IADEh+TMSal+KkX+BJ41ds3P2SJph74T1Udjrsr+XoRGQ7cCsxU1SHATO+5iVe8jcUAGzcmPh5jTMrxM8TEXeDmJfBGHvVFVTcBm7z13SLyGdAPOAcY7+02FZgF3BJT1GaflpYIjDGdnp+JacaJyKfA597zkSLyp1guIiIlwJHAXKCXlyQANuOqjpo7ZrKILBCRBdu2bYvlcp1LPG0ElgiMMSH8VA09CJwGlAGo6hLgJL8XEJGuwD+Am1V1V+hrqqp4w1s3papTVHWMqo7p2bOn38t1PvGUCLKzoUcPqxpG4q1iAAAdBUlEQVQyxgD+EgGquq7JpoCf40QkE5cE/qaqL3qbt4hIH+/1PsBWn7Ga5sSTCMDuLjbG7OUnEawTkeMAFZFMEfkJ8Fm0g0REgCeAz5r0MHoFmOStTwKmxxizCRVvIujb10oExhjAXyK4Frge19C7ARgFXOfjuOOBK4CTRWSxt3wLuBc4VURWAhO95yZe8fQaAisRGGP2CttrSEQGqOo6bxKay5q8diauS2lYqvo+4ec2PiXWQE0Y8TQWA/Tr5xKB3V1sTKcX6RvgLa+3TyMi8n3goWQFZGIUb9XQwIFQV2elAmNMxETwY+BNERkS3CAit3nbv5HswIxPLUkEAGvXJjYeY0zKCZsIVPV14IfAv0TkcBF5EDgLOElV17dWgCaKqirIzY39uEGD3OOaNYmNxxiTcqINOjcT+D7u7t8DgZNtXKB2pL7eVe/EkwiCJQJLBMZ0epEai3fjbvYSIBvXwLvV6xaqqlrQOiGasKqq3GOXLrEfm58P3btbIjDGhE8EqupnYDnTloKJIJ4SAbjqIUsExnR61m8wlVV6YwDGUyIAOPhgWLUqcfEYY1KSJYJUFkwE8ZYIhgyBr75ybQ3GmE7LEkEqa0kbAcDQoS4JrF6dsJCMManHEkEqa2mJYOhQ97hiRWLiMcakJEsEqaylJYIh3r2ClgiM6dQsEaSylpYIiovdnMVLlyYuJmNMyrFEkMpa2mtIBEaPhkWLEheTMSblWCJIZS29jwBcIli+3A1nbYzplCwRpLKWlggAjjzS9Ryy6iFjOi1LBKmspY3FACee6B5nzmx5PMaYlGSJIJUFSwQRJqaprIT58+Ef/4B//Qs++mhf/gCgVy8YNQr+/e/kxmqMabfCjjVkUkBVlUsCTWYYq6mBF16ARx+FDz90k5CFys6G446DK66ASy+F7NNPh//+bygtdT2JjDGdipUIUlllZaOGYlV48UU46CC4/HLYvBl+9jNXGvj4Y1camD4dbrjBvXbVVXDggfAMF7t2gmefbcM3Y4xpK1YiSGWVlXvbB0pL4eqr4ZVXXE3Pn/8MEyc2Px3x2WfD/ffDW2/Bb34Dl9w7kqNyR9D30afIu/76Vn4Txpi2lrQSgYj8WUS2isgnIdt6iMhbIrLSe+yerOt3Ct7sZMuWwdFHu2r+++93bQLf/GbkOelF3D7vvAPTpsFf068k75N5PP3zz1FtvbdgjGl7yawaehI4vcm2W4GZqjoEmOk9N/GqrGR3oAvHHefaBd59F37yE8iIoZwnAhdfDDfPu4wAaaz9zdNcdBFUVCQvbGNM+5K0RKCq7wLbm2w+B5jqrU8Fzk3W9TuDbWsqWL66CwcfDAsWwNix8Z+r+6G9STv9NG7s9jQvvtDA+PGwdWvCQjXGtGOt3VjcS1U3eeubgV6tfP0O45VX4Kslu9H8AmbOhL59W35OmXQlBeXrmH33bJYvh3HjYOXKlp/XGNO+tVmvIVVV3JzIzRKRySKyQEQWbNu2rRUja//efRcuvBAOyC7nqAkF9OiRoBOfcw4UFHD8l0/xzjuwa5frZjp3boLOb4xpl1o7EWwRkT4A3mPYygdVnaKqY1R1TM+ePVstwPZu6VLX62fwYBhUuIusooLEnTw3F777XXjhBY45opIPP4Ru3WDCBFcCMcZ0TK2dCF4BJnnrk4DprXz9lLZxI5xxBnTt6noIpe3Z5b6pE+mSS2DPHnjjDYYMcTekHX44nHee65JqjOl4ktl9dBrwEXCIiKwXkauBe4FTRWQlMNF7bnyoqoJzz4Xycnj9dRjYL+C+sAsSWCIA+MY3oKjI3ZqMm67g7bfdPQlXXw333Yd1LzWmg0naDWWqekmYl05J1jU7KlV3F/CCBfDSSzBiBFC+x72Y6ESQkeF+/j/zjBuaOieHrl3hn/+E730Pbr3V9Sa6//7I9ykYY1KH/SmngN/8xn0v/+Y3rj0XcC25kPhEAK6dYM8eePPNvZuysuCvf4Ubb4QHHoBJk6CuLvGXNsa0PksE7dyLL8Ltt8Nll8Ett4S8kMxEMGEC9OwJDz7YqB4oLQ0eegjuucclhXPOsRvPjOkILBG0Y4sXuxFCjzkGHn/c3QW8VzITQWYm/PKXbvyJG290g9Ht2AG4GH7+c5gyxTVYT5zoxjkyxqQuSwTt1JYtrptojx7w8svNTDmQzEQAcO218IMfwCOPuDEojj7adVvyXHMNPP+8G9X0mGPcbJfGmNRkiaAdqq527bWlpW7Y6N69m9mpvNw9JisRZGTAY4+5hPPmm7BuHdx1V6Ndzj8fZs92g6COG+cSljEm9dgw1O2Mqvsh/tFH7hf36NFhdkxwiaC0spS15WsbLVsqtlBRW0FFXQUVP+lBxfYpBP74HulZ2aRLOulp6WSkZTDk7i4sW5TPeVPzGTkvn2+Oz6cwN5/8rHwKsgvoltONbtnd9nvMTM9MSOzGmJaxRNDO3H47/O1vrkH2ggsi7NiCRFBaWcqCjQuYv2E+8zfOZ8HGBWzas6nRPrkZufTq2ouuWV3Jy8wjr99Aeny9mfSKLALFAwhogAZtoL6hngqtoN8Rmwn03M2Smt0s/WAXmlYfNY7cjNywSSI/K5/sjGyy0rOaXVSVBm1otCiNt6nq3jgjLZlpmeRk5JCdkU12eja5mbl0z+lOUZcieuT2oCi3iOIuxWRnZMf8WRuTCiwRtCNTprguotdc42YWi6isDNLTIT8/6nkbtIH5G+bz2srXeG3layzatAgAQTik+BBOOfAUjux9JIMLBzOw20AGdhtIcZdiJLR1WhUePARKD4A3w4838eKLcM1kZU9VDT/9xW4uu2oXlYFyyqvLKa+J8Oitr9u1jvLqcvbU7qE2UEtNoMbPRxcTQUiTNNIkDRGhLlCHhh/2aq+++X0pKSyhpLCEwYWDGdJjCCN6jeDQnoeSkxF+3mhj2jvRFLhNdMyYMbpgwYK2DiOpXn8dzjoLTj/dtQtEnVPgBz+A116DTZuafTnQEGDW6ln8ddlfeW3Fa2yr3EaapDGu/zjOOPgMjh94PKP7jKYgO4YSxU9/6rqUlpVFTEBbtsB117mkcOihrsvpqaf6v0yo4K/62kAtNfU11AZqqQ3U7v0iD/1Cb/QcQUT2VmE13d70GnUNddTU11ATqKGyrpIdVTsoqyqjrLKMsqoyNu/ZzJryNazeuZrVO1ezrnwdAQ0AkC7pDC0ayoheIxjTdwzj+o/jqL5HWXIwbU5EFqrqmKj7WSJoewsXupEdDjnENb527erjoDPPhA0bXLedEJ9t+4ynljzFX5f9lfW71lOQXcBZQ8/iW0O+xWkHnUZRl6L4A33nHTj5ZJepzj474q6qLk/dfDN8+aVLBHfd5RqVO4L6hnq+3P4lS7csdcvWpSzZvIQ15WsAyEzLZHSf0Rw34DjGl4znpEEnUZhT2MZRm87GEkGK+OILlwRycmDOnDA9hJpz9NFuTKA33qCmvobnP32eh+c9zNwNc0mXdE4/+HSuHHklZw09i9zM3Ojn86O21l3z8svhf//X1yHV1a4H6n33wbZtMH48/OhHrvSTnp6YsNqTLXu28NH6j/ho3Ud8tP4j5m+cT3V9NWmSxlF9juKUwadw8uCTOX7g8XTJ7NLW4ZoOzhJBClixwn0xBgKuJDBsWAwHDxzI+olj+b/LD2XKoilsrdjK0KKhXHvUtVx6xKX06pqkOX/OPdfd6fb1103ucItszx6XOx5+GNauhX794NJL3TJyZEynSik19TXM3TCXmV/N5O3VbzNn/RzqG+rJTMtk3IBxexPD2H5jyUrPautwTQdjiaCd+/xzOOUUN17PO+/AYYf5O05VeW/Nu/zx1gm8NAwaBM4ceiY3jL2BiQdOJE2SfGvI//2fu9nss89izFxOfb2b2+DJJ+Ff/3LPBwyAb30LTjjB3Zx28MEdNzHsqd3DB2s/YObXM3n767dZtGkRitIlswvHDzieCSUTGF8ynjF9x1j3WtNilgjasTlzXBV/ejrMmAFHHBH9mMq6Sv629G88PP9hlm5ZSvcquLpgPD/84RMc2P3A5AcdtHYtDBoEv/2tG4q0BbZtc0nhtdfgrbdcqQHc3dRjx7pRVgcPhgMPdI8DB0J2B+vBuaNqB7PXzObtr99m1upZLNu6DIC8zDxOGHgC40vGM75kPEf2PtK6r5qYWSJop6ZPd3O/9Onjbtg96KDI+3+942semf8IT3z8BDurdzKi1whuGPAdLjvvDrr85a9uNLrWdtxxbrS5JUsSdspAAD791CXJuXPdsmKFa5YIlZcH3bvvW3Jz3dBIGRnuMS3NlTKCS11d4/Xa2saPTddVXYLOyHCPmZmu8b6gwHWUys+H4mL379enj2vT6dfPJari4paXZEorS5m9ejazVs9i1ppZfLL1EwCy0rMY3Wc0x/Q7hmP7H8vYfmMZXDh4vx5QxoSyRNDOBAJw553uRrExY+DVV6FXmGp8VWXGVzP447w/8uqKV0mTNM4/9HxuGHsDJw48EZk503XDefttN1Joa3v4YTcY3Zw5ri4nSRoa3PBGX38NX33lRrnYvt2Nfxdcqqsbf+EHAvuSQkbG/utZWe5508fgIuLOEQi489XWupLK7t1u2bXLlWS2bt1/gp6CAld6Oegg12328MPdMnSoO3c8tlVs47217zF3/VzmbJjD/A3zqaqvAiA/K58RvUYwstdIRvYeyeEHHM7QoqEU5RZZgjCAJYJ2Zd06N7HMjBnw/e/Dn/7UzCBywOY9m3l6ydP8efGf+bz0c3p26cnkoyZz7Zhr6V/Qf9+Of/wj3HSTu4fAdzejBNq92/0EPvJIV6zphF869fUuIWzaBOvXu0T15ZduWbXKPTY0uH0zM13X4MMPd9WAwQRRUhL75D51gTo+2foJCzYuYMmWJW7ZvITdtbv37tM9pztDi4YypGgIJd1KGNBtAP0L+jOgwD0W5hRaougkLBG0Aw0Nbp7f//xP9wvz979394GF/g3W1NfwzxX/5MnFT/LGqjcIaIBx/cdx7ZhrufCwC5u/Kenaa+G559yNXW31Bx0sFfzsZ66YY18sjVRXu67Bn3zSeFm9et8+ubmu5DB8uOssMHy4WwYPjq1rraqyeudqPiv9jBVlKxotG3ZvoEEbGu2fl5lH/4L+9MnvQ3GXYopzi+mZ19OtdymmZ5d968VdihPX/di0OksEbWzWLJcAFi1y9wn8+c+u2gCgoraCN1a9wctfvMyrK15lZ/VO+uX348qRVzJp5CQOKT4k8slPOsnVS7z3XtLfR1jB0fH+/Gd3X8HjjyenJbe+3sdt1gkWCCTtJofdu11byCefuMfly93junX79snJcR2yhg93/2dKStwyaJDrYRXLx1zfUM+m3ZtYt2sd63etZ12597hrHVsqtlBaWUppZSlllWVhh9nIzcilqEsRRblF+x5D15t5LMwpTH4PNhOV30RgYw0lUCDgesDcfz+8/77r5fL3v8OFFyqrdqxgysLZvLbyNd788k2q66vpkduDcw45h4sPv5hTDzyV9DQfXz6q7tsj4oh0rUDEffkPHgy/+IXrTfTaaz5vi/Zh0SJXj7Z0qcukjz/u+pUm07JlLrnNmwfHHuv6uB4SJSnHKD/fNas0bVrZtcv1yA1NDh984KYobQj5QS8Cffu6/1sHHOAmkuvZs/F6jx7uOl27Qn5+Bn27DmBAtwER4wo0BNhRvWNvYiitLGVbxTaXJKrKGg23sXTLUsqqythetX2/0kZQmqTtHbivuEsxvbv2pk/XPvTu2nu/9QPyDrCusm3MSgQtpOrur3ruOTd94/r10H9wFRdev5wB4+bx0cbZvLvmXTbv2QzAgIIBnDvsXM4bdh4nDjqRjLQYc/GcOW6chqlT4cork/CO4jBtmisVTJjgWsGbawCJxezZ7tbjwkI3f/KTT7pf6K+/7lrak2HePDjtNBf7pZfC00+77f/+t2sLaSN1dW4kkdWr3bJmjXtcu9a1UQSXQCDyebp0CU0OrloqO7vxkpUV+XnottxcyMpuoCGrnNr0MmrSyqiSMiq0lD0NZewJlLGrroydtWWUVpayec9mNu/ZzI7qHfvFJsjeZNE0UfTq2oteeb32PhZ1KbKSRgzaddWQiJwOPASkA4+r6r2R9m9PiaCy0v1ymztXeWf+Ft5etJrtgTVI8Sp6j1yK9F7K5roVe38p9cvvxzdKvsE3BrllaNHQljXU3Xyzu0V361bo1i1B7yoBnnrKzWh/1lkuK8abDF580X0RH3iga4ju39/1Iz3tNDdTz4svxj+CXTjvvQff/rb7OT1zpquHWbHCXWfHDjfjzsknJ/aaCdTQADt37ksKO3a4KqjQ3k7B9eBjZaXrEVVT03hpuq1p9914BLvfFhRA18Iacoq2kFG4ibSCzWjXzQRyNlObvZmq9E1UyGZ2NWxmZ/1mahv2H3k2XdLpmdezUXJomii653Sne273vY+d+Y7tdpsIRCQdWAGcCqwH5gOXqOqn4Y5pSSIINASob6iPuOyurGPNxkpKyyso213BjooKdlZUULqrgm07Kyit2M726m2U12+jklLI2woF6yCzutG1BhcOZkSvEXu79B3Z58jE9vWeOxdOPBEuvNAVP9qbP/0Jrr/ejRlx222uLaO4OHrfycpK17L69NNuqNKxY13JoihkgLyNG+GMM1wWvv12OOccV22TnR1fQ3VlpevaM20a/Pd/u8Qzc6a7KSBo/Xp3zS++gP/6LzdtXEOD27e4OPZrpiBVVyppmjCqq6Gqyi2h680937PHVX01twS75DbsV8OkkLMTum6BvC3QdQtp+VvI6rGF9G6bScvfQkPeFgI5W6jN3EJDWvjhyrOkC3npheRndCc/szsFWYXkZeWRl5lHl4wudMnMIy8rjy6ZXeialUfX7C7ueVY22RlZ5GRmkpOVSXZGJtmZWWSlZ5KZnklWWhaZ6ZlkprnnGWnpiMjeYc6D6yKCNqQRqBcCASFQL2Skp1FYmPwOFu05EYwD7lTV07zntwGo6m/DHRNvIvj237/N6ytfjzfUfRoyyKwrposW0z2nJ70Lijm8/0BGDCxhcPeSvWPUd81KUP14OL/7HTz6KCxY4CqC26NXXoEbbmjc+pmV5fpJqrq/+NDH0P9/Iq6f7R/+4Ooymtq507UbhM6JmZbmSh/B84SeM7je3Guh3zyXXup6QXXv3vw1r70Wnn1237a//93dFWgSQtXl5XBJIvQejubWy3cpu6p3sathCzVp26lL3wG5OyAnzGN2OWRWQlYFZFa49YwEFH3iFJo43Ofh/o8GG++XX7ecYcWxD+cC7TsRXACcrqo/8J5fARyjqjc02W8yMNl7egjwRasGur9ioLSNY4gmFWIEizORUiFGsDgTzW+cg1S1Z7Sd2m2vIVWdAkxp6ziCRGSBn8zallIhRrA4EykVYgSLM9ESHWdbNL9vAEL7svX3thljjGkDbZEI5gNDRGSwiGQBFwPhJ8E1xhiTVK1eNaSq9SJyA/BvXPfRP6vq8taOIw7tppoqglSIESzOREqFGMHiTLSExpkSN5QZY4xJHrtFzxhjOjlLBMYY08lZIvCISA8ReUtEVnqPzdxdBCIyydtnpYhMCtk+S0S+EJHF3nJAguM73Tv/KhHZb45IEckWkWe91+eKSEnIa7d5278QkdMSGVei4hSREhGpCvn8Hm3DGE8SkUUiUu/d9xL6WrP//u0wzkDIZ5nUzhg+4vyxiHwqIktFZKaIDAp5rVU+zxbG2J4+y2tFZJkXy/siMjzktfj/zlXVFtdO8jvgVm/9VuC+ZvbpAXzlPXb31rt7r80CxiQptnTgS+BAIAtYAgxvss91wKPe+sXAs976cG//bGCwd570dhhnCfBJK/w7+4mxBBgBPAVc4Offvz3F6b22J9mfZQxxTgC6eOs/DPk3b5XPsyUxtsPPsiBk/WzgDW+9RX/nViLY5xxgqrc+FTi3mX1OA95S1e2qugN4Czi9FWIbC6xS1a9UtRZ4xos3VGj8LwCniLtn/RzgGVWtUdWvgVXe+dpbnK0laoyqulpVlwJNR8BpzX//lsTZmvzE+Y6qVnpP5+DuHYLW+zxbEmNr8hPnrpCnebB3EokW/Z1bItinl6pu8tY3A83NKNwPCBlEh/XetqC/eEW2XyT4yy3adRvto6r1QDlQ5PPY9hAnwGAR+VhEZovIiW0YYzKOjVVLr5UjIgtEZI6INPejJlFijfNq4F9xHhuvlsQI7eyzFJHrReRLXC3GTbEcG067HWIiGURkBtDcJL8/D32iqioisfarvUxVN4hIPvAP4Apckd34swkYqKplInIU8LKIHNbkF5Dxb5D3//FA4G0RWaaqX7ZlQCJyOTAG+EZbxhFJmBjb1Wepqo8Aj4jIpcDtQIvbVjpViUBVJ6rq4c0s04EtItIHwHvc2swpwg6PoarBx93A30ls9YufYTn27iMiGUA3oMznsW0ep1ekLQNQ1YW4Os6hbRRjMo6NVYuuFfL/8Stc+1WyZtfxFaeITMT94DpbVWtiObaNY2x3n2WIZ9hXhd2yz7I1GkFSYQHup3Fj8e+a2acH8DWuYau7t94DV7Iq9vbJxNV9X5vA2DJwDWmD2deIdFiTfa6ncSPsc976YTRuRPqK5DUWtyTOnsG4cI1lG4AebRFjyL5Psn9j8X7//m31WUaIszuQ7a0XAytp0ujYyv/mR+IS+5Am21vl82xhjO3tsxwSsn4WsMBbb9HfecLfTKouuHrqmd4/9Izgf0hcMfHxkP2uwjXErAK+723LAxYCS4HleLOvJTi+b+Em9PkS+Lm37W7crxeAHOB5L655wIEhx/7cO+4L4Iwkf45xxQl8x/vsFgOLgLPaMMajcXWsFbhS1fJI//7tLU7gOGCZ98WwDLi6jeOcAWzx/m0XA6+09ucZb4zt8LN8KOTv5B1CEkVL/s5tiAljjOnkOlUbgTHGmP1ZIjDGmE7OEoExxnRylgiMMaaTs0RgjDGdnCUCk3JCRoNc4o2+eVySr1ciIp/42O9OEflJkuO4NOT590Tk4WRdz3QelghMKqpS1VGqOhK4Dfht0x28u5Y7mhLg0mg7GRMrSwQm1RUAOwBEZLyIvOeNGf+pt+1lEVkoIstFZHLwIBHZIyK/9koVc0Skl7e9l4i85G1fElLaSBeRx7zzvCkiuX4DFJHLRWSeV4r5PxFJjxLDQd7zZSJyj4js8U51L3Cid57/8Lb1FZE3vPH8fxf/x2g6M0sEJhXlel+GnwOPA78KeW008CNVDY5TdJWqHoW7Q/wmEQmOdJoHzPFKFe8C13jb/wDM9raPxt3FCTAEeERVDwN24u6EjkpEDgUuAo5X1VFAALgsSgwPAQ+p6hG4O4eDbgXe80pDv/e2jfLOfwRwkYiEjjdjjC+WCEwqClYNDcONX/9UyLDf89SNxx50k4gswY0xPwD3hQ5QC7zqrS/EVbsAnAz8L4CqBlS13Nv+taoubmb/aE4BjgLmi8hi7/mBUWIYhxuGA9wAhpHMVNVyVa3GlYIGRdnfmP10xHpU04mo6kciUowbtA7cuDuAqyoCJgLjVLVSRGbhxjoCqNN946sEiP63UBOyHgD8Vg0JMFVVb2vmtVhj8BOX/U2bmFmJwKQ0ERmGm+KvrJmXuwE7vCQwDDjWxyln4qYqRETSRaRbC0OcCVwg3hzW4ubGjvarfQ77qp4uDtm+G8hvYTzG7McSgUlFwTaCxcCzwCRVDTSz3xtAhoh8hmtonePj3D8CJojIMlx1zfAo+zd1u4isDy6q+ilu8pA3RWQpbjrGPlHOcTPwY2//g3GzuIEb3TbgNS7/R9ijjYmRjT5qTDsjIl1w7SAqIhcDl6hq07mfjUkYq080pv05CnjYawDfiRuz35iksRKBMcZ0ctZGYIwxnZwlAmOM6eQsERhjTCdnicAYYzo5SwTGGNPJ/X9qtAek5rvb5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "colors = {'all':'blue', 'internal':'red', 'terminal':'green'}\n",
    "for k in sorted(branch_lengths.keys()):\n",
    "    distplot(branch_lengths[k], kde=True, hist=False, color=colors[k], label=k.capitalize());\n",
    "plt.title(\"Branch Length Distributions\");\n",
    "plt.xlabel(\"Branch Length\");\n",
    "plt.ylabel(\"Kernel Density Estimate\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that, as expected under the dual-birth model, the distribution of internal branch lengths is smaller than that of terminal branch lengths. Specifically, we see that the vast majority of internal branches are very close to 0 length (i.e., are extremely short).\n",
    "\n",
    "We can compare these three branch length distributions numerically as well by comparing their averages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average Branch Length: 0.043235\n",
      "Average Internal Branch Length: 0.009187\n",
      "Average Terminal Branch Length: 0.076943\n"
     ]
    }
   ],
   "source": [
    "print(\"Average Branch Length: %f\" % tree.avg_branch_length())\n",
    "print(\"Average Internal Branch Length: %f\" % tree.avg_branch_length(terminal=False))\n",
    "print(\"Average Terminal Branch Length: %f\" % tree.avg_branch_length(internal=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to looking at the branch length distributions of trees, it is common to look at the *height* (i.e., maximum distance from the root to a leaf), *diameter* (i.e., maximum pairwise distance between leaves), and *treeness* (sum of internal branch lengths divided by sum of all branch lengths). These values give you an idea of the distances between elements in your tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Height: 0.224848\n",
      "Diameter: 0.445494\n",
      "Treeness: 0.105712\n"
     ]
    }
   ],
   "source": [
    "print(\"Height: %f\" % tree.height())\n",
    "print(\"Diameter: %f\" % tree.diameter())\n",
    "print(\"Treeness: %f\" % tree.treeness())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to distances on a phylogenetic tree, the topology of the tree is typically of importance to us. In the case of the dual-birth model, the topology of the tree is significantly affected by the ratio $r={\\lambda_a}/{\\lambda_b}$, and the balance of the tree can be used to estimate this ratio. Two common metrics to do so are the Colless Balance Index and the Sackin Index, both of which are described by [Mir, Rossello, & Rotger (2012)](https://arxiv.org/abs/1202.1223)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Colless Balance Index: 0.318697\n",
      "Sackin Balance Index: 19.060000\n"
     ]
    }
   ],
   "source": [
    "print(\"Colless Balance Index: %f\" % tree.colless())\n",
    "print(\"Sackin Balance Index: %f\" % tree.sackin())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The rates of speciation down the tree are also inherently of interest. Given a phylogenetic tree, we may want to test the null hypothesis that speciation occurs at a constant rate over time. [Pybus & Harvey (2000)](https://doi.org/10.1098/rspb.2000.1278) developed a statistical test, the $\\gamma$-statistic, to compare a given tree against the null hypothesis. The $\\gamma$-statistic is essentially a summary of the information contained in the internode intervals of a phylogeny. If $\\gamma>0$, a phylogeny's internal nodes are closer to the leaves than expected under with a constant rate of speciation, and if $\\gamma<0$, the internal nodes are closer to the root than expected. In this example, we saw earlier that the tree we simulated under the dual-birth model had much longer terminal branches than internal branches, which implies that we would expect a very negative $\\gamma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gamma Statistic: -7.520125\n"
     ]
    }
   ],
   "source": [
    "print(\"Gamma Statistic: %f\" % tree.gamma_statistic())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save Dataset\n",
    "Now that we have finished exploring our tree, we can save it in the Newick format. The tool that we wish to use next doesn't support the \"rooted\" prefix of the Newick format (`[&R]`), so we need to specify to TreeSwift that we wish to omit the \"rooted\" prefix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree.write_tree_newick(OUTPUT_TREE_FILE, hide_rooted_prefix=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next step\n",
    "After you saved the dataset here, go back to the [0-Workflow.ipynb](./0-Workflow.ipynb), or go to [2-SimulateSequences.ipynb](./2-SimulateSequences.ipynb) to run the next step of the analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Author:** [Niema Moshiri](https://niema.net/), UC San Diego, October 2, 2018\n",
    "\n",
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}