{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reed-Frost Model in R\n",
    "\n",
    "*Author*: Sangeeta Bhatia\n",
    "\n",
    "*Date*: 2018-10-01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `xcurr` is the number of suspectible individuals at step j.\n",
    "- `ycurr` is the number of infected individuals at step j.\n",
    "- `q` is the probability that a suspectible individual \n",
    "escapes infection.\n",
    "- `ynext` is the number of infected persons at step j + 1.\n",
    "\n",
    "Everyone who is infected at step j becomes immune (or dies) at step j + 1. Thus at step j + 1, the number of susceptibles is `xcurr-ynext`. The following function gives the probability that the number of infected persons at \n",
    "$j+1$ is `ynext`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "reed_frost <- function(xcurr, ycurr, q, ynext) {\n",
    "  pescape <- q^ycurr\n",
    "  pynext <- choose(xcurr, ynext) *\n",
    "      (1 - pescape)^ynext *\n",
    "      pescape^(xcurr - ynext)\n",
    "  pynext\n",
    "}    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of infected people at a time step can be no greater \n",
    "than the number of suspectibles at the current time step.\n",
    "So a sanity check for our code would for this number to be 0.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0"
      ],
      "text/latex": [
       "0"
      ],
      "text/markdown": [
       "0"
      ],
      "text/plain": [
       "[1] 0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reed_frost(10, 5, 0.5, 11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the probability of escaping an infection is 1, then the \n",
    "probability that the number of infected at next time is 0 is 1.\n",
    "Similarly at the other extreme, if probability of escaping an\n",
    "infection is 0, then every susceptible becomes infected and \n",
    "the number of infected at next time is `xcurr`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>1</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 1\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 1\n",
       "2. 0\n",
       "3. 0\n",
       "4. 0\n",
       "5. 0\n",
       "6. 0\n",
       "7. 0\n",
       "8. 0\n",
       "9. 0\n",
       "10. 0\n",
       "11. 0\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " [1] 1 0 0 0 0 0 0 0 0 0 0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>0</li>\n",
       "\t<li>1</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 1\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0\n",
       "2. 0\n",
       "3. 0\n",
       "4. 0\n",
       "5. 0\n",
       "6. 0\n",
       "7. 0\n",
       "8. 0\n",
       "9. 0\n",
       "10. 0\n",
       "11. 1\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " [1] 0 0 0 0 0 0 0 0 0 0 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Probability that ynext = 0 should be 1.\n",
    "## 0 everywhere else.\n",
    "reed_frost(10, 5, 1, 0:10)\n",
    "\n",
    "## Probability that ynext = 10 should be 1.\n",
    "## 0 everywhere else.\n",
    "reed_frost(10, 5, 0, 0:10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That works as expected. \n",
    "\n",
    "Now for various values of q and \n",
    "\\(y_curr\\), we can check the behavior of the model.\n",
    "First we will load the needed libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Attaching package: ‘dplyr’\n",
      "\n",
      "The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    filter, lag\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    intersect, setdiff, setequal, union\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(ggplot2)\n",
    "library(dplyr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pescape <- 0.9\n",
    "ycurrs <- c(20, 40, 60, 80, 100)\n",
    "probs <- purrr::map_dfc(ycurrs,\n",
    "                             ~ reed_frost(xcurr = 100,\n",
    "                                          ycurr = .x,\n",
    "                                          q = pescape,\n",
    "                                          ynext = 0:100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAIVBMVEUAAAAAsPYAv30zMzNN\nTU2jpQDna/Pr6+vy8vL4dm3///9NgKTRAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4\nnO2djXqjyK5FM2c83Td5/we+sRE2P1WFVNqABHt/5yTdHXtFYNaABMZfPwzDuPN1dgEMc4VQ\nJIYBhCIxDCAUiWEAoUgMAwhFYhhAKBLDAEKRGAYQv0h/m9n6uSlAVtS6whaWtS6AI6pQJAQM\nyQpaWNa6AI6oQpEQMCQraGFZ6wI4ogpFQsCQrKCFZa0L4IgqFAkBQ7KCFpa1LoAjqlAkBAzJ\nClpY1roAjqhCkRAwJCtoYVnrAjiiCkVCwJCsoIVlrQvgiCoUCQFDsoIWlrUugCOqUCQEDMkK\nWljWugCOqEKREDAkK2hhWesCOKIKRULAkKyghWWtC+CIKhQJAUOyghaWtS6AI6pQJAQMyQpa\nWNa6AI6oQpEQMCQraGFZ6wI4ogpFQsCQrKCFZa0L4IgqFAkBQ7KCFpa1LoAjqlAkBAzJClpY\n1roAjqhCkRAwJCtoYVnrAjiiCkVCwJCsoIVlrQvgiCoUCQFDsoIWlrUugCOqUCQEDMkKWljW\nugCOqEKREDAkK2hhWesCOKIKRULAkKyghWWtC+CIKhQJAUOyghaWtS6AI6pQJAQMyQpaWNa6\nAI6oQpEQMCQraGFZ6wI4ogpFQsCQrKCFZa0L4IgqFAkBQ7KCFpa1LoAjqlAkBAzJClpY1roA\njqjSJdLj8fj8xbecptxguwhbWNa6UKJspUekx+t/Y3zLacoNtouwhWWtC2bKRjpEery/vOJb\nTlNusF2ELSxqXX/+/Gn/qoNCkRAwJCtoYZHq+u+//8avf/5smAQzZSM+kf79DbQc5pL5/v7e\n/Dp+a+ZXnteX32///e9///vvz5A9a1eGIjHorCT5fmbj6+TbT0mqj0KS/z1zFZGe2diJBzoi\nmKGC1hW2MFVdv9v/68vvt8lXc+TJA244hPtv4dDv///5559RpGbdx4QiIWBIVtDCanUNG/tE\nIXDmCn32RU+PRpP+a9Z9TDj+RsCQrKCFrev6uLOPQt/fv9J81z36Neiff/7b8Ci0SDwhu4Qh\nWUELG+p6WdO3+5ke7P2s/mXBeypUcEgO6l4ivfZF/8gEr1X3MeElQggYkhW0sO9F+2PQZubf\n6+vP6l/mv6Gg0DPf/3sf1b27o8QnZBfZWg7EqzjCgKigdQUszLj/WWtjquu5hyko9D336J9/\nKJIjFMkK8z1dq9Df+Z7HU9faoeHg7SnR98QjiuQJRbLCep9oUuj9BF9d5X2R/PB/n/wzMSn1\ntXaLbC2HYf1uhSJZYfanbCtk3f0o6yooVPRonNd9/359tlsUyR6KZIVZHqxQqNn4dNf13+qM\n62wW9+vOKNLf3z2RXF83sbr9qw4KRULAkKyzCtMdwu1Q11ud5eHckP8tj+qkO5pV1/hVB4Ui\nIWBI1vGFNXdFo0K71LVui1bnhVYeFURqmQRwRBWKhIAhWUcW9lGoINJ8L7RHXUuDSudWx4M6\nObCbzOsmVbd+1UGhSAgYknVMYaq90I51rfZFfwvXKEy7o2fm8zpZDIpkD0Wywgr/VlFoaxaH\nras25Z7lvSt6evQaM7w8mr0rlj1STyiSFbb+J8teaJe61vuilkdyUFdQ6LM8zV8HcEQVioSA\nIVn7FVZuhTYPjrB1qfZFr7QP6rR1ARxRhSIhYEjWHoW1OiLlZBtUl2pfVOuOSm/jo0gdoUhW\nGEIhWF3N00WzzLqjZyhS1A02al3gwhAKoerS7YuemXVHQyhS0A02al24wiodUe8lCq66lvui\nBmt5UCfzusnFQba6AI6oQpEQMCTLD2sfzp1Q12pfVGdNL2QY/kXmdUPWUzuK1BGKpAnucA5U\n13pIV2VVPPqYZK4L4IgqFAkBQ7I8sPLhHOSK0866ygOGLZH+zg7sXiLVbrpFkTpCkdrBdUSo\nuioDhiJr3R09Uz+DpKsL4IgqFAkBQ7L6YOV90ZkrrDXsLrFWI+/pZUH1G0FSpI5QpHW+x7fa\nlQ7nTlxhzWF3gVUbeQ86UaSoG2zUuqyFvQ/giodz562w9lVADZE+mY4ZKFLUDTZqXYbCykdz\ns5bonBW2fQWDWaT6x7dQpI5QpElKDi1HC6esMMUVDHPW4M/So8Xgm1O79s9NoUhDikdzpfHc\nGSusvS8S2PQvo0ELjxpvnbDUBXBEFYqEgCFZTdh38aKFv7UzRUevsMlBXfNquilrfUz3d+uy\nIEtdAEdUoUgIGJLVgimO5vYqTLPCNDujATb5c0mkZXPUMIkideTOIpVHC8cVtrHCtgcMpbrK\n52CV8zpVXRSpCAOigtZVK8y6L4IX1l5higFDoa7COdhXKJJ9OU25p0hd+yJ4Yc0V1j5rVIC9\nvhbOwQ7RDb4Vdf2lSGUYEBW0rnVhBYW0Hh2xwv5b3m9YX1dxzPDKxlsnNHV9fn5MKBIChmRN\nYMUx90mFVVbY+2BOeVA31lXpjmbzOk9dk58fE4qEgCFZH1jf4dxOhZVX2MQgpUMDrNYdWfZF\nrbqmPz8mFAkBQ7JeMPe+CF7YeoVNd0XqfdEQSHdUq2v582NCkRAwJOsJ8++L4IWtVpi1L5qm\n2h1Z5nWVulY/PyYUCQEDskD7InhhyxXW7VG1O3qGIvUupyk3EAm1Lxqy0wrrmS+MeRn0UxvX\nWa5pWNdV/vkxoUgIGIQC3RchCxtQHxbioK7g0Wpexx7JsJymXFok7L5oyC4rzOOR9tzREE7t\n9MtpykVF+l5d0o26Zwl8hXkO6l5RdUf2ulo/PyYUCQFzPfutDnJfNAS9wlw7o+FwbuiR1j/U\nvY2vUlfz58eEIiFgvU9ctUXA+2e5CiugfgAeiUmluuxnkD51tX9+TCgSAtb5vKVBT4dCFFYI\n9KCuVJd5XjeEInXkQiKtR3R7fnq4P76d0TN1kUSg2YEdRTIupynXEak+oou5wpweLU/Czupa\nDxkokn05TbmESBuni+KtMPdB3foNfNO6SuM69kjm5TTlCiJtnS4Kt8JwB3WfsfeWSJzamZfT\nlOQibeyLziusFb9HhXNHmyLpQ5E6klukrX3RaYVV4z+oe6Ut0qJH0u+LBEWR7EkrkmpfdEZh\nzcx2Rp11fU7C1utaeGQziSJ1JKtIun3RCYU1Mz+o66trchK2VFepNzKaRJE6klQk5b7o+MLa\nmR/UddVVu7JOWL7Bt6Aokj0JRVoe1IUprJ1pdyQovEiFMQNF6l1OU/KJZNkXHVpYO+9d0WfC\ncIxI7JF6l9OUVCKZ90VHFbad0qTOXFd5zDCty3kGSVcXwBFVKBICVvg3+77ooMI2sj6oE5SR\nVRkzzOpynUBS1gVwRBWKhIDN/7oedgcpTJHq+VdjXfV3wg51Fe4Dad0X6eoCOKIKRULAZn+z\nDLsPLUyR+nUMWJHKx3Q9JlGkjqQQyTTsPrKwrbSvY4CKBJnX6eoCOKIKRULAxj90Dhj2L0yR\njYvqDHW1xgxDKFIhvuU0JbpIvn3RjoUpsnVxqr6u5phhCEUqxLecpgQWCbAv2qcwbbYuTlXX\n1T6oExh7pHV8y2lKXJEQ+yJhHb/C/pt8xFEdBRWp9LktnNp5ltOUsCL1D+nWrMNXmAi08Y4j\nVV3Lt5SX8lTowBcS4IgqFMmZpzk/OI+OX2HvXVH7HUeaulZvKS/kdVA3Z/XtipR1ARxRhSL5\n8nLnB3JQN+TQFVa7jKGE2q5r/ZbydWpjhk14b10AR1ShSN2ZzBf+IvZFQ45cYVuTuhlKL1Ir\na5F6x3XKugCOqEKRejM/nEN5dOQKs3hEkTZCkToDnC/McrxIutsxbNS1fRJWsuqRKJLEt5ym\nRBFpedIIuozHrDBDdySoZl2Kk7DTD0CasdgjDfEtpylBRFqdNMonUuGde1uoVl2ag7rpmyY4\ntSvEt5ymnC9S+QKGdCKZDuoE5RVp2h3xPFIhvuU05XSRKhcwpBLJfFAnKIrUDEUypDZgyCSS\naVQ3RVXqGvzRHtktRHId1LXr+vz8mFAkZcoHdfC69l5hvR7V6hoN2h7XFXok35ihWdfk58eE\nIulSOaiD13WQSPY7EJfrUk0ZxnzOHw0s5+C7Vdf058eEIqnSPmmURKTO7khQHpGW13pTpEJ8\ny2nKOSItLgXat64dV5h95D1DOURaXV5HkQrxLacpp4ikuIIhhUjdB3WCWtVlu5hhZhJ7pEJ8\ny2nK0SK1Bgw71bXPCnMd1AlqWZfqYoZXqiJxajeNbzlNOVik5oBhp7p2WWHdo7opalGXYcxQ\nFwkQitSRY0WqnTUqwFBV/d1lhSE88ohU6pEA+6JyXeufHxOKVIz2oA5e144ieT50r1ukwh1V\nnzBEd1Sua/3zY0KRSlEf1MHrQq8wQHckqEld6jFD7e7ekHnduq7yz48JRSpEvy8SGKioFwsK\n8428Z6hPXZ4xwxCKVIhvOU05QqTts0a71oVdYZCDuiGfulxjhiEUqRDfcppygEjWndEAAxX1\nYqFgsIO6IVCR2CMV4ltOU/YXqcujkCJBRnWTdIlU65E4tSvEt5ym7CtSz0EdvC7UCkN7JHXp\nxwyS8ieJ8TxSIb7lNGXX9d+3MxpgoKJeLKhIiO5oyKsu/ZihppDAMDW962r+/JhQJInDo2gi\ngbujIc+6AAd1w2VBFKkQ33KastP6/5Yb03Uc1MHrAqywycgbvMIAY4ZhYEeRrpiXQONdh88u\nxh/ZFz2/o9EikuahItLq32Xwja7s/HCPNNkV9d7pMdQeaXpQh15hxiO7fc8gTepq/vyY3Fuk\n6aDOcdfhMCItuyNgYepbM0j2vzhoCEXqCFokz4BhBgMV9WI5YKsLgnCFWUbeQ+qzBvZIxfiW\n0xTw+kd5FEWk9cgbVljnTU5K4dSuEt9ymgJk9Z99LSSYSBOYv54hkMH3JxSpEN9ymoJjwXZG\nrwQQqXzu6AyRqoPvyWVBFKkQ33KaAmNhPQogUuXtEpDCjJcFtc8gAesSFEWyB8KCHtQNOV2k\n2gVBiMIslwU9oxl8U6RCfMtpCoIF3hm9EkakFcxfz3hQp69LMfimSIX4ltMUAGsPj04V6T/5\nQPI4IlWmdhRpI77lNMXJ2uGgbsiJIolAlStUjxapOa5jj9SObzlN8bFmO6NAdS1YJth7V1S+\n0ttZ2HTMoKlrY/DNqV0zvuU0xcWaH9TFqWvJUsMU75bwFTYbMyjqqg++sXXNURTJnm7W+qAu\nRl0llhameQ+sq7D5uSOK1M4tRCpMGELUVWQpYar3kkcRaXmZKkUqxLecpnSySpO6CHWVWTaR\n2u8lP1akao+0uuCbIhXiW05TOli1Sd3ZddVZCpj6veTdha2vZnBM7dZvnaBIhfiW0xQ7q3ra\nKLNIleuBSrDOIgpXM2zU1eqNKJIqvuU0xcyqn35NLJLqoE5gfTWULlFt19WcMlAkVXzLaUqv\nSIXTrxcQSQPrq8Es0sa8jj2SJr7lNMXEmnZHBVROkdTdkcD6akCLxKmdJr7lNMXCeu+KytcC\n5RRJ3x0JrLOIwpsmXCItQ5EK8S2nKQbW1jV1KUUydEcCs//+2k1OHD0SpK4qiiLZo2S1D+oE\nlVkkPcz866tv4Ouc2pXvGESRCvEtpyk6luqdEulEMnZHArP+9vpbymt1bd7ppGQSRSrEt5ym\nqFgqj9KJZO2OBGb97WaRNocMRZMoUiG+5TRlk1W7juHsuiysIszcHQnM+tutImnm3hRJF99y\nmrLF0u2MBlROkawwy4PbNzmhSO1cSSSDR5lE6uqOBGZ47MZNTrpEYo9kiG85TdGJpHobeR6R\n+rojgekfunXnuq4eiVM7Q3zLaUqDpRl5z1BZROrsjgSmf2inSBtTuxrM/pQqiiLZU2dtXMdQ\nQGUTqQ+mf6hdpC6FBNb7xAKKItlTZVkO6gSVQSRHdyQww2M37qW6WmF9B3X2ujZCkTpSYn1P\nP7fSgEogkqc7EpjqUSJQ+16qyxXWOWaw1KUKRepIgTU5pruaSK7uSGCaB+lu620Uqf2BYhSp\nEN9ymrJmTY7pbHd8TCSSB6Z4jPKDJihSOxcRyXznVIo0plOkzouD9HUpQ5E6Mmf1HdONqNgi\n/de4D7EFpnhMr0h9l6vq61KGInVkxuo8phtRoUUShbwebRbWvipoXZdENfjm1M4Y33KaMmXZ\nJ95zVGSRAAd1Amv/eOOqoFVdEuO7+Ox1mVAUyZ6SSL0oimT7hOXPCrO+r9xclw1FkewZWa7u\nSFBRRXKfhJ3mHJFaB3WaukyhSB0Rlq87ElRQkfwnYac5RaTmmEFTlykUqSMDy9kdCSqmSICT\nsNPoeiQdStsjtQffqrosoUgdmYvkQ8UWCZR6YbWbBdVR2qkdReqMbzlN+UF0R4KKJxK0OxpS\nLcyyLxKUdiEpUmd8y2nKD6I7ElQ4kbDd0ZBaYabuSFAvlmZaxx6pL77lNAXSHQ0JJxK4OxqC\nFkk39+bUriu+5VQHdlA3JKxIhxTWKRLgDFK7rh4URTLlexIEL5JInwvrjhDJcFnQDEWRmkki\nEtqjSCJ9Lqx7ftu7MMtlQTOURqTtg7pqXb2hSKYAu6MhcUR6D+qG7mjnwjoO6gSl6JEUY4Za\nXd2hSOqAu6MhAUUSVmSR/GeQanV1hyJpMxl5H7j+bbDeJxZOHUUVSdEbUSRXfMu5melB3dVE\nKp06OqhHskYzZaBIrviWczPTg7qLiVQ8dbRjYebLgj7RzevYI3niW85mlt3RRUWas3YrrHNf\n9Ipy8M2pnSO+5WxldUEQRTJmyurujp6BnUEaQpEK8S1nI+uR94VEqt7dJKZIgLeXT0ORCvEt\nZyPrkfd1RKrf3SSoSFtTO+1B3RCKVIhvORu5sEiNd0vE65EGhdp1qccMhbp8oUjNfL/v6D2D\ndbGKuZdIjnndeFDXrEs/+J7XBQhFakUUWl7LcAmRNt6/t0dhiHkdRWonpEi1C4KuINLW+/d2\nKAwyr6NI7VAkBEz/0M3376UUiT1SPJFal6heSKQ6K5hIqh6JU7twIjXvykCRjAH0SLqpnTEU\nqRDfci7Sft9RcpGqJ2FnLGxhnnndHEWRmgkqUgVmYjVzwnZRPwk7Y0ELQ+yLBFWpy3ZMN8I6\nCyqhKNIy22/gSy2S8pZ10MKAlwVV6jJOGUZYZ0UlFEVapNkdCUzL2g5F2sz8QtVyXda59wjr\nq6iIokjzaO7KkFakjZOwMxZFsqEo0jyauzJkFWnrJOyMhSrsKZCrR6JIllAkBKz9482TsDMW\nqLCXQj+ueR17JEMiiDT4o7hNUHKRdCxMYc43TQzh1E6fACKNBm3fbiuhSIbuSFghRFq/+4jn\nkdo5XyTDLevyiWTpjoQVQaTC+2EpUjsUCQGr/cDUHQkL2yN1pXSHhnVdXQd1Aut9YgFFkSR3\nEMnC8hf2uSxoT5H6xgwC63xeCUWR/tbeCFuHqR6lQ11XpMkx3Y4idQ6+XXUVURSp+kbYOkz3\nMBVqd5FUl6gWWN7Cpt1RL0vRI1Gkec4UyXBQJzD9QzdRe4uku0S1wIogkmJqR5HmoUgI2Pqf\nOg7qhBVCpHXYI7Vzlkjbl3qXYPqHbqIuK5KzR6rdv45Tu3ZOEklxqXcJZnjsFuqSIi3fxmdn\nVe+oyvNI7ZwjkuZS7xLM8uAN1I7bReeYQViOwlbnYM2s+j2+KVI754pkXWnWJzRQ+20XvWMG\nYfUXtr6agSJRpCLM+oQGarftovugTlgJRHJ0RwLzPX2GurVI5u5IYOZn1FEUqRhNj+SZ1wnM\n9ew56rYiaS/1LsHsT6midhHJfK13gXVqj6SZ2rnOIAnM8+QF6q4i9e2LBNb5vBJqD5Hs13oX\nWH2FlW+7ZWK1P7iFIrVztEid3ZHA+p5WRO0gUse13gVWV2GV90xYWBsfJUaR2qFICNjwzXlQ\nJ6yewmrvPjKwtj7ckj1SOxQJARu+3UUkTu0KUYr0eDzW34bol9P4pokCrPN5JRRYJNdJ2Bkr\ngUju3Fakx1Oe4dvP+5tEvZzWN00UYL1PLKCw24XvJOyMFb5HAuTOIv3+f9gbjd/GaJfTdVAn\nMM+TFyjodgE5qBOWsTARqHx/BvjUzn1QJzAERFCZRHpQpDbsPJHatzhBrzD/mEFgAMaIyiTS\nM4+SSP/+RgsQkYz1hc+vPz+jSMf/dmmOfJDfXZHugTL49v22i2ZbpGGy8Bgbo26Rfq7q0duk\nE349QqRXc6R6JEWqRzu1m3zpOLRzXBY0g/mePkNhWO9jOtCBnbWwjdvXaVhb47pJXYBTsfq6\ntKhUh3aPz9cekdzN0QgDMEYUWCRYDu+RDCKxR6pHN2z4fO8ZfwPGDALzI96o/CKVL6+bwxQY\ni0ic2lWjnNo9xk6p54TsRUX6nIOFTgC1MM0tiVUslUc8j7SRIy4RuqZIk+boDJFUN/fWsTQe\nUaSNHHKt3RV7pOlBXVqRVAoNAR3UDaFIhbSXAzSvk5WCgAiKImkP6l5BjRmGUKRCmouB2hfJ\nSoGR8osE6JGUY4ZnYIPvIRSpkNZSwLojWSkokFOk5aXeR4ukmNcJrPlTigQMRbJndQ72YJH0\nHyFGkShScaWgQC6R1ueOjhXJ8GF87JEuItIle6TriMSpHS7HTO1QKw2Ioki28DxSO+d/9KUh\n528X5XeUJ+yRDPsiQSV9IQGOqEKRLKlc6n2cSOp5ncBqPzB0R6q6jKFIhfiW05SzRapdonqY\nSPp9kcAq/26Y143dUdYXEuCIKhTJkLNFMnRHAqv8u0UkmddlfSEBjqhCkQy5o0jjGaSsLyTA\nEVUokjKtG9elE8nQI1EkXSiSLs13lKfrkQxTO4qkC0VSpf1O2CNEMs7rBAYoiD2SKhRJldNF\nsu6LBFb4N/MZJE7tNKFIqpwtkrk7Etj6n+xnkASV9IUEOKIKRdpMa8wgMEhJwtpXJNMZpO26\nOkORCvEtpymniKS4cd1VRZpeqJr1hQQ4ogpF2ojmhluJeqSOU7HNuvpCkQrxLacp9xSpa14n\nsPU/2c8g1epyhCIV4ltOU44W6T/lh5TvKlLfvkhghX8zn0Gq1OUJRSrEt5ymHCzSYJDiVqp7\nitTZHQls9jfjkIEiWUKRqhn3Rdu3JM4hknlcxx7JEIpUjeagTmD+ej6svUTqGHxzaqcPRaom\nhkiwHqn7DJKgkr6QAEdUoUjFbJ+EncEgJQlrAnPM6wQ2+TNF2jUUqRTFSdgZDFKSsD4wz75I\nYNO/uDxK+kJSpDIMiGqx9Ad1AvPX82G9Ya7uSGCzv7luvpXyhfxLkcowIOqGIulTuB1kyhfy\nL0Uqw4CoW4nkOoO0rAsQilSIbzlNOWL928YMAoOUJKwdeiRrd0SR7KFI8xjHDAKDlCSs1dTO\nA3t9Nc/rKJI9FGkW60GdwPz1fFgvmFshgb2+2gff7JHMoUizBBHJf1AnsNdX5zUNk7pQoUiF\n+JbTlLuIBBgzCGz45juDJKhcL+Tk58eEIs3T41F4kew3PCmgkr2Qn58fE4r0zmBQh0fxRTKk\n9jlIiV7Ixc+PCUUa07UvEpi3mikL2yNZ90XVT+bL80Iuf35MKJKkrzsSmL+ed9wXqk7zgziD\nJKgsL+Tq58eEIkmCiITaFw2BnEEakuaFXP38mFAkSQyRYN3REIpEkYowIGrCek8Zzu+RzhaJ\nPVJvKNLboG6P4opk7pE4tesNRfIc040wfz1jsD2SfWpXR4V/IWs/PyYUKYxI0HndkAM3WBsM\niKJI9lxaJNkXwVbYc19EkShSEQZEvVg9bz4qwQAFjd0RaoW9uiMTq9YdDYn9QrZ+fkxuLVLX\nm49KMEBBYJGA87ohoV/I5s+PyZ1FAhzUCcyPOF2k+hmkIZFfyPbPjwlFCiISuEeiSO+fHxOK\nFEUkmded1SNRJF/uKhJozDDEWdds4n3a1I49kis3FQk1Zhjiq2t+Dva8FcapnSf3FAl2UDfE\nVdfiqiDACvv0RjyPRJGKMBTo0iJNpgwUiSIVYSjQlUWazusoEkUqwhAQ6JhhSKgeyS5SuzmS\nxHshBUWR7EGwsGOGId11Fa5TPV6kjXGdJNwLOaIokj0A1nhQF6Ku0nsmDu+Rtk4gSUKssBKK\nItlzMZGK7+I7fGpHkSChSIAEE2kCUzyGIkFyN5HGHilCXXiRllfXsUeiSEWY7+nTWzOEqAvd\nI62uU+XUjiIVYa5nz0beZ9dVfV+5o7D1Fd88j0SRijDPk+cnYU+uq36PE4pkRFEkey4jUuOu\nW4eKpDqoE1SYF3KBcoj09TX9+vX6w+fffv///t92KBIgQUQy90i6MYOgwryQCxRGpK8vMYki\ntbK+LOiSIhmndsrBt6BCvJAFFEqkjz7vP4lhKg1uIVLhsqAr9kgFWPOntxfpY87HltkeafaT\ndu4gUula77Prqt4HsrOw8v0ZKBJFKsL6nhZLpI1bqfYVVrnTCXskilSE9T0tlEhbN/fuKqx2\nzyBO7TbG319iAEVSpfDuo7Pq2vy4iUNFMiTEC1lC+UQaZwnzYcMXRVpmelnQDEWRbDn7hayi\nfCdkR08m4+8vmXhTpEmq74S9lkidPZIlFxfpfUL2szeiSJ/U781wrR7JPrUzdEeCurhI3lAk\nQJQs1QcgHVOYZV4nKIrUBrkJvuU0JbdIug/jsxbWvMF3jWU6gySoK4qkvWxBkeuK1L5b0PF1\nKT8e1lhY+1b5FIkiFWGGx27cLegiIm186ARF4tsoijD9Q7duAXlvkdgjwUORADmvR+oViVM7\ndCgSIKap3RbriB6pIxSpnauKtHVT4mPrUikkrCOmdh2hSO1cUaTaZUEz1JF16Q7qhJV0g7XB\ngCiKZI+OpbpB/pF1KccMwlIXpviI2BLL3B0JiiI1cz2RdB/Zkl8kzYctF1j2eZ2gKFIzFAmQ\nM0RSfWz5mtVxBslWlw4GRDlE+r96zBpQJEDO6JEokqAokj2brPZlQTNU9qkdRRIURbJni7Vx\nWdAMdVBdBoWExR7JhrqMSHEiB3VnlzHN66BuJ/avRz1P+/UIXUnWxBLJ9x8MUzZYuu5IUIfU\nZRozCCvpf/ltMCDqMnsk33KaQpGGKM4gCczymzdQFKmZK4lk+YBymL0AACAASURBVKzyzCJp\nuiOBWX7zBooiNXMVkTSXBc1Qx9Rl9khRmGpeJ7DZ3zqvaRAURWrmIiIZ9kWC2r0u1R0aCqzd\nROqd1wnqviJNby90cZEs3ZGg9q7Lvi8S1l4idZ9B0tZlgQFRu4s03p2rLQtFAqTA6uiOhLVX\nj0SR2iJ9f39XRfr9/9eGLRQJkINF6pvaUaSmSN/fM5MmhtxJpHg90q4i6WHTv7BHaoj0/T03\naWHJDUQSgYweZe2R1Psigc3+xqldp0hf7w+tuK5I5l3RiDpoamdntQrTd0cC66mggrqvSF/z\nL5cUyd4cjagd6+pUSFiNwgzzOoH117FCXVukeo8kilCkCmq/unoP6oRFkWyovad2X5/vFx5/\nBxSpe8wgrD1EcnVHgrq6SPNMDBk/8OXaJ2Tj9Ug7itTbI/nmdYK6q0jKzEX6eieBSNbL6+ao\nlCL1Te2cZ5AERZHUIn3NEl2k3n2RoDL2SHbY6ytFOlakmTx6k3zLacqU1d0dCSrb1M64LxLY\n6ytFOlSkpTlak3zLaUp4kVwKCatcmLU7EtjwjT3SkSL1xrecpkQXyXdQJ6xiYeZ5ncDkO6d2\nxVCkZ+L1SM4xg7D2EAkQitTOSiTbyO7nFJFc8zpBUSQjiiI1s2qMrNPvE0Ty7YsElUskX4+E\nyBVFQmZ1HmnyTRffcpoysJzdkaBy9UjdUztAdyQoitTMUqT5d018y2lKYJE679BQYEELQ8zr\nBHVBkfY7tKNIfTDIvkhYyMIgZ5CGUKRLiaS/S/4WCloXpDsS1qqwroO6IRTpvGGD1aNDRTLc\nJX8LlUakvjHDEIp0ypUNPVfaHSoS5KBuSBqROgffAmOPRJEKCSvSjj2STyRO7Y4XqTe+5TQl\nrkiIeZ2wsCLBQpHaSSXSD8wjXF1Pg6DLiOyRKBIvESrCEGMGQYHqGj5JDMN6BTq1o0i8RGiZ\np0LxtgvglEECW2HP7ijeChMYEHWESOab6Me9ROh1UBdvu9hXJMe+aJh8x1thAgOiDhBJ7lps\nuItQ2BOywDHDkAwiBTmDNOTyIv1uXhWRvt4Sqe9rR5Gs2bFHcs3rKNL75zqRXtsXRTInwdSO\nIqlQGJFkAyuJNPlsJINI7JHU+RzSBRSJPdL7506RvnpECjm1+7wfNtR2MWmOdhp/ezzi1O79\nc7dILyHM9/4Odx5pckwXabuYjhn2Oo/k8WiA+Z4+Q11bpFaP1LNH6ohvOTcz7Y4ibRdHiOSH\nAVFJ69KK1JjajZ9+2TH+pkib2Vck975IYAiIoK4u0jwFS2wnZCmSNnv2SP7uSGCQkgbUnUVS\nZGmO3STfcm4naI+059QOMK8TGKwsirSR0J9Gsbx/XZTtYnktQyiRpmeQoqywFQyIokibWZ2D\nDbJdrK4KokhWGBAVUqSe+JazkfXVDDG2i/X1daF6JIq0+Hk9FMmf2CL5r/sWGKSkAUWRmgl8\nQvbOIrnCqd3s5/XsKFKoS4Ru2yPhYEBU0rrOESnMRavLeZ3AuljFhJ3agU7FDomywlYwIMoh\nEjIx30ZRec/E+dtF+U18wLp8p2KXOX+FVWBAFEWqp/buo9O3i8rbYXF1+U7Fru5id/oKq8GA\nqKCHdhSpkdoby6OItHov39krrAoDokKKFOTe3xSpK+t3xZ69wqowICqkSD23LfYtZzlBe6T9\nRQLf8OTsFVaFAVEUqZXy/RlO3y5275HAH+Jy+gqrwYCokCL1xLecq7RucXL+drH71M61wbJH\nKv28niuL1LxZ0JnbRevudZi6hn2Ra4Pl1K7w83r2EWl5KHfKoV37tlsnbhfN+0BC6pLuKOsG\na4MBUeFEml/PoL+4wbeci0QVqX1HVURd47wu6wZrgwFR8UTq/Kgx33IuQpGybrA2GBC1v0iq\nC+eivbEvaI8UXqTyXVUpkl6k35e3JtLsbqsqkXriW85Vgk7tgvdIlfsTUyS1SK8XuCjS1+dr\njvvabd/bm1O7Smp3+qZIWpHkkKMs0lfPvb/Nx3a+5fxEcZf8c7aL7Y9tOf08EkVq/Nwt0suk\nbpGOntppPm/ilO1C8QFIFMkKA6IO2CPJJ7uYbqIvPhkuXPUt5ztRRdJ8lJizrtlVQeyRjKhD\neiT7Hunz/eATsvcVaX6dKqd2RtTuU7tkIkXtkXYXafHOiawbrA0GRO1+HunrPW2w3/v7WJFE\noKBTu717JIrkQ0U9ITt+O65HUuyKRpjqUTpUmKkdRfKhAl4i9DO9HddhImmaoxGmeZAuqu1i\nWyGBuUrx9UjND4qlSOe9jeLwG0TGFUlxUCcwXy2eqV37I5cp0o3ejxRWJM2YQWD+ej4sE2zj\nw8sp0o1ECtsjHSBS4X3lFMmIiirSCff+1np0OZFKdzqhSEZUTJEOvve3WiGBWR68gQrQIxXv\nvcUeyYgKKdLB9/7WH9QJzPDYLVSAqR1CJE7tgoo0/76vSIYxg8D0D91EtVlahQTWVwNEpGYo\nUlMkZChSKeqDOoF1FuHvkTZCkShSEaZ/6CaqxdKPGQTWW4V7arcRisQeqQgzPHYLFUKkEksN\na3ZHAvPVMkNRJINInNq9kkKk9rxOYN5qJiiKZBHpsPNIRoUEZn9KFXVyj1S9xbd2g904gyQw\nW1FNFEVq5qQrG6wHdQIzP6OOOndqV//QCYpkRN1aJPOYQWDWJzRQNZZNIYFZn9D4GCSKZESF\nE6nvPqtXE8l4UCcw6xMAIrFHEhRFiiiSdcwgMOtvR4jEqd2ACidSb7qWM2qPdJBIgB5JFYp0\ndZGiTu2OEsk/tVOFIp10z4bjROrKZXqkFivZCuuDAVEokX7/01YWSe67EO4m+l37IoH1PrGA\nOm1qt/EhsZoNVtEdCUz3MBXq4iK9DrZrItnvtFrP4/XlMf3WJVJfdySwzueVUGtWj0ICMzx2\n62PLFRusZl4nMGVRGtS1RZLxT32PhBLppc5j8Em+dYnUOa8TWN/TiqgVq+ugTmD6hzbmdcLa\nhKnOIJkL20TdVqShOQKJ9Hjb8xi/XUykvjGDwPQPpUgvGBB1xB4Jd2hHkZow/UMp0gsGRKXq\nkR4/ZZH+/Y0K8MnLI+NzjoiItPvveb1gPsTLI0w19456agcR6fF4zRZQIv3E9EhMOuD3uD16\nmYSo5Pb5v3omj+q5iX4lj0Em76GdY/AtMN/TZ6izpnabrKSHUDYYEJXthKy/R/IMvgXmevYc\nNWH1GzTCfE+fs5JusDYYEJXtEiH3+Ns1ZhCY58kL1IflmHuPMN3DNk7FCivpBmuDAVEZRXKd\nkI0qkmdcN8JUj9o6FSuspBusDQZEZROpGu1y3l2kzcH3srBS1BcHCczy4A0URWrmwGvtgvZI\nmUTSXxwkMMNjt1AUqZkjL1oNOrU7qEcCiGQ4FWsoTImiSM3wbRSHTe38PRJFKqBuJZJ7XyQw\nBERQJ2wX7qkdRSqg7iSSvzsSGIAxol4s975IYAjIyGKPZEM5RELmCJEA8zqB+RFv1JPl744E\ntvUA1b5IWJza2VAUqWOl+RFv1A9kXiewjZ/ruiNhJd1gbTAg6kaHdncXSTmv+xQGC0W6lEhR\neySKZEElrau16aYTKerU7qAeCSKSsTsSmP0pVRRFaube55EOmtoBeiTrvE5g5mfUUaFfyNbP\n66FI/oAUGrL/1M58Bklg1ic0UEFfyLuIBDqok5UCI6EO6obsfx6JIlVR9xAJNWaQlYICwcYM\nQyiSFQZE3UIk2OBbVgoKdJxIhoM6YbFHsqEOEGmQpP2Oc4oESL0uy5hBWJza2VD7i/Q13kHo\n66d+D5R7inRUj2QafAsr6QZrgwFRKJF+/0tVFmnYBYlM1bty3bNHOmpqR5EqMCAKJNLr2Llx\naHeySAGndoNCx2wXFKkCA6IwIsk0J65I4da/HNQdVBeiR+rqjgTW+8QCKtoLOaIokj0A1jhm\nOKou/9Sub14nsM7nlVDBXsg3iiLZk1AkO2sB6zyDJDB/PW9U0BV2lx4p2vqnSJ2JusL2n9p9\nPhvpvPF3wPV/WI9kPqgTFkWyoW5xQjbi+j9oamcfMwiLPZINdYtLhPKufxts/U8dg29hcWpn\nQ1Eke5ys6UnYRCI5QpEoUhHmevbssiCKZEUlrau16VKknswvVM3TI3lCkShSEeZ58sEioaZ2\nrlAkilSEeZ58tEi9rAnMMWYQmO/pM1TQFUaROpKoR+pnfWCewbfAvNVMUEFXGEXqSJKpXedB\nnbDeMNepWIF5nrxABXohZyiHSMjcQ6T1u492rKt3zCAsimRDUSR7elmF98PuV1f34FtYFMmG\nokj2dLJKd2hIIBJ7JBWKItlzM5E4tdOgKJI9KUSC9UiAUCSKVIR1Pu/QHgk2tQOEIlGkIqz3\niYdO7ZysoIVlrQvgiCr3EKmA2qUu175IWC+YuzsSGAIiqKQvJMARVa4uUu3+dbvU5euOhPWE\n+ed1AgMwRhRFaubiIlXvqLpHXc55nbB+IGeQZoUhQpHaubZI9Xt8UyQriiI1Q5EAoUhWGBBF\nkewJLhJ7JDUMiKJI9gTvkTi1U8OAKIpkT/SpHYQVtLCsdQEcUeXqIlVRQesKW1jWugCOqHJd\nkdqfgISuC3BQJ6ykG6wNBkRRJHssrI3P5APXhRgzCCvpBmuDAVEUyR4Da+tTYrF1QQbfQ0Bj\nhiEUiSIVYfqHZhUJNfgeQpEoUhGmf2hSkWCnYodQJIpUhBkem7NHokhmFEWy5/pTO4pkRlEk\ne5SstkKCCrpdsEeyoiiSPTrWxkGdoIJuF5zaWVEUyR4Va2vMIChUXc9jOl7ZYIUBURTJnoAi\nvaYMFMkKA6Iokj3xRAKeiZVk3WBtMCCKItkTr0dCijR0R1k3WBsMiKJI9sSb2mGvDXqalHWD\ntcGAKIpkT8DtAtYjjWeQrr7CBhgQRZHs2WRp9kWCija1o0idKIpkzxZL1R0JKtp2QZE6URTJ\nng2Wbl4nqHDbBXukPhRFsieWSJ8hA2YZObXrQlEke0KJNBnX8YSsFQZEUSR7IvVI08E3RbLC\ngCiKZE+kqR1F8sCAKIpkT6TtgiJ5YEAURbKnwdLviwQVqEeavXMi6wZrgwFRFMmeOsvQHQkq\nztRu/l6+rBusDQZEUSR7qizLvE5QYbaLxbvL4xS2QCWtC+CIKhQJEIpkhQFRFMmeKCItr/em\nSFYYEEWR7AnSI63eOcEeyQoDoiiSPTGmduv3IHFqZ4UBURTJnhjbBVykWS64wgowIIoi2VNi\nWXdFI4oiGVFJ6wI4okp2kczN0YiK0yPNk3WDtcGAKIpkz5plH9eNqBBTu8LdILNusDYYEEWR\n7Iki0grW+bzS/YlDFFZCJa0L4IgqFAmQTlbxjvkRCiuiktYFcESV5CId3SOVb71FkawwIIoi\n2XP61K5yEzuKZIUBURTJnrO3i9rtINkjWWFAFEWyZ87q3BcJKoJInNoBUBTJnhmrtzsSVAiR\nCsm6wdpgQBRFsmfK6p7XCSpAj1RM1g3WBgOiKJI9p4uEndoVk3WDtcGAKIpkz/kiVWBIVtDC\nstYFcEQVv0hn5eXR2UV058+fP2eXwCCTdo909NSu9SlI9mWsf3h51v/y22BA1GX2SL7lNOXE\n7aL5eWLmuoqnYoV1kRXWhgFRFMmekeXaFwnKVlf7k/kokhUGRFEke4TlO4MkKIpkRCWtC+CI\nKglFcs7rBHWqSOyRcCiKZM95IoF7pOLFQcJKusHaYEAURbLnRJHAU7t6sm6wNhgQRZHsOa9H\n2oAhWUELy1oXwBFVMop0wtRuA4ZkBS0sa10AR1RJKRICpWe1DuoEZvnN1e5IWPlXmAIGRFEk\ne34Q+yJBqetqjhkEZvjF9XmdsJJusDYYEEWR7PlBdEeC0tbVHnwLTP97G2eQjIVpQpEoUimQ\ned0QimRFJa0L4IgqFGkjFElQSesCOKIKRdoKe6QBlbQugCOqpBLpjB6JUztBJa0L4IgquUQ6\nY2qngSFZQQvLWhfAEVWSiYRDBa0rbGFZ6wI4okoakZ77oqO3i+2DOoH5apmzkm6wNhgQRZFM\nGe7QgGE9o6lLMWYQmLeaKSvpBmuDAVEUyRLgvG6Ioi7N4Ftgut+5MWYQVtIN1gYDoiiSJVcQ\naWvwLaykG6wNBkRRJEsuINLmqVh1YfpQJIq0SP4eiSJNYEAURbIl/dSOIk1gQBRFsif3dsEe\n6QMDoiiSNp/eKPl2wandGwZEUSRlJlOGG2wXYQvLWhfAEVXCizSd1x23/rXdkcB8tcxZSTdY\nGwyIoki6nCKSel4nMG81U1bSDdYGA6Ioki5niKQ/gyQwfz0fVtIN1gYDoiiSMif0SHCRVGMG\nYSXdYG0wIIoiaXP81A4tkm7wLaykG6wNBkRRJHuS9kjKU7HCSrnCrDAgiiLZk3RqR5FWMCCK\nIm1neZlq0u2CIq1gQBRF2szqgu+s2wV7pCUMiKJIW1m/deKI9W87qBPY1gM4tVvAgCiKtJVT\nRDKOGQTmrWbKSrrB2mBAFEXayhkiWQffAvPX82El3WBtMCCKIm3mhB6JIlVRSesCOKJKZJFO\nmNpRpCoqaV0AR1QJLdIKBkQd1CMZxgzCSrbC+mBAFEVqpXyjk3xTO8vgW1hJN1gbDIiiSI1U\nbhmUbrswnYoVVtIN1gYDoihSPbWbb6XbLihSBQZEUaR6zhGp56BOYLUfUKQKDIiiSPWcIlLX\nmEFg1Z+wRyrDgCiK1MgJPVLf4Ftg9R9xaleEAVEUqZXjp3Y7iWRnJd1gbTAgiiLZQ5GsMCAq\naV0AR1ShSGN26ZE6WGlWmAcGRFGkclqfOZFqamfujoSVdIO1wYAoilRM89NbMm0X9nmdsJJu\nsDYYEEWRSml/DlKi7aLjDJKwkm6wNhgQRZFKOUWk/mO6Ebb+J4rUhAFRFKmUM0RyTBlG2Pqf\nKFITBkRRpGKO75E8c+8RVvg39kgtGBBFkco5fGq3k0ic2rVgQBRFsieVSL2s8CsMAQOiKNIy\n259anqhH6mcl3WBtMCCKIi3S7I4EpmVtZ++pXT8r6QZrgwFRFGme9rxOYEqWIlG3i7CFZa0L\n4IgqFAkBm/+1c8wgrKQbrA0GRFGkea4jUu/gW1hJN1gbDIiiSIuc0SO5uyOBTf/SfSpWWEk3\nWBsMiKJIyxw/tfPP6wQ2/QtFUsCAKIpkD3j9A84gCWz6F4qkgAFRFMmeFCKxR1LAgCiK9M72\nMd0I0z1MhdpPJE7ttmFAFEUao5gyjDDVo3So3XokLyvpBmuDAVEUSaKZe48wzYN02W9q52Ul\n3WBtMCCKIknOFAkFQ7KCFpa1LoAjqlAkBAzJClpY1roAjqhyvkjn9Eigg7oh77pcYwZhJd1g\nbTAgiiK9c8LUDjVmGDLW5Rt8CyvpBmuDAVEUyR4YCzb4HiJ1OU/FCivkCstbF8ARVSgSIBTJ\nCgOiKNJfw0GdwCwPboUi2VFJ6wI4osqpIunHDAIzPLYd9khmVNK6AI6ocqZIhsG3wPQP3Qqn\ndlZU0roAjqhyV5GibhdhC8taF8ARVSgSAoZkBS0sa10AR1S5aY8UdrsIW1jWugCOqHK/qd3Q\nHUXdLsIWlrUugCOq3O48kszr0HUBxgzCirbCRlTSugCOqHI3kcYzSOC6EINvYQVbYW9U0roA\njqhCkQD5gZyKFVawFfZGJa0L4IgqZ4lk7I4EZn/KMhSpF5W0LoAjqpwkknVeJzDzM9bZpUei\nSEYYEHVrkcxnkARmfUIpu0zt2CPZYEAURTpJJEFxamdEJa0L4IgqFAkBQ7KCFpa1LoAjqtyp\nR5peqBp1uwhbWNa6AI6oohTp8Xisvw3pW84Tpnazt05E3S7CFpa1Lrwy5WhFmvxv/CbxLacp\nLtb8zXxx6lqyghaWta5drClEJ9Lj8/UxfhvjW05TIor0nDJQJCsMiEom0uOBEqnroE5gvU98\nZh+RXnNvimSFAVG5RHqZtBbp398Yf+FrzGB8DiYvj8BMORMLpjL5si3SQxTCiCSDb2OVmOA9\nokiMRN8jlUR6xrbn7TyDJLC+pxVRGBbw2iBJ1kMoGwyIynVo90ORymGP1AUDonKJ9J42uMff\nDo9Cbhec2vXAgKhUIiFPyJ4xtVvffCvqdhG2sKx14ZUp5xZv7CvcDjJEXUVW0MKy1gVwRJU7\niFS6QXGEusqsoIVlrQvgiCoUCRCKZIUBUTcUydEdCazvaTuJ9Bl7UyQrDIi6n0ieeZ3AOp+3\nS480OYFEkawwIOp2IrnOIAms94k7TO2mp2IpkhUGRFGkjpXmefICRZGMqKR1ARxRhSJ1hiJ5\nYEDU7UQ6sUcqodgjGVFJ6wI4osrVp3a1DxTj1M6KSloXwBFVLn4eqfoRl1G3i7CFZa0L4Igq\n1xap/qHLUbeLsIVlrQvgiCoUCRCKZIUBUbcSyd0dCcz6BIqEQyWtC+CIKoeI5J/XCcz8jH16\npOV7YimSFQZE3UgkwBkkgdmfssfUbvXucopkhQFRFKljpfkRb1Q/a32fBopkhQFRFKljpfkR\nbxRFMqKS1gVwRJXr9ki1gzpBUSQjKmldAEdUuezUrjpmEBR7JCMqaV0AR1S56nmk+uBbUJza\nGVFJ6wI4ogpFAoQiWWFAFEWyhyJZYUBU0roAjqhyVZH27JHWMCQraGFZ6wI4osreIoHGDLJS\nLA/ebWpXgCFZQQvLWhfAEVV2Fgk1+JaVAiP1bhflG+ZTJCsMiLqFSLBTsbJSUKDe7aLy0RMU\nyQoDoihSx0pDgTq3i9qHuFAkKwyIokgdK033sHZ3JCiKZEQlrQvgiCoX7JE25nWCokhGVNK6\nAI6ocr2p3dYZJEGxRzKiktYFcESV651H2lMkTu1AMCCKItlzvkgVGJIVtLCsdQEcUeV6Iu3Y\nI1VhSFbQwrLWBXBElQuKtN/UrgpDsoIWlrUugCOqXFEkFSpoXWELy1oXwBFVKBIChmQFLSxr\nXQBHVKFIypTndQLz1TJnXWWFNWFAFEWyZ5Ol6Y4EZayrcgZJYDZWM1k3WBsMiKJI9myxVPM6\nQdnqql3TIDATq52sG6wNBkRRJHs2WLozSIKiSEZU0roAjqhCkVShSFnrAjiiCkXShT1S0roA\njqhyJZF27JE4tctaF8ARVS4l0o5TuzYMyQpaWNa6AI6oci2RDKigdYUtLGtdAEdUoUgIGJIV\ntLCsdQEcUYUiIWBIVtDCstYFcESVq4ik744Epa+rNWYQmOlXb7CSbrA2GBBFkeypswzzOkGp\n62oOvgVm+c0bybrB2mBAFEWyp8qynEESlLau9qnYjbo6knWDtcGAKIpkD0WywoCopHUBHFGF\nIm2EIgkqaV0AR1S5hkjskSowICppXQBHVLmISJzalWFAVNK6AI6ochWRzKigdYUtLGtdAEdU\noUgIGJIVtLCsdQEcUYUiIWBIVtDCstYFcESV9CJZm6MRFXS7CFtY1roAjqiSXSTzuG5Eaera\nHjMIrKOAKivpBmuDAVEUyZ41y34CaUQp6lIMvgVm//11VtIN1gYDoiiSPceKpDkVW6urP1k3\nWBsMiKJI9lAkKwyISloXwBFVkou0Z49EkWaopHUBHFElu0h7Tu3YI01RSesCOKJKepF6UZza\nGVFJ6wI4ogpFQsCQrKCFZa0L4IgqFAkBQ7KCFpa1LoAjqiQWqbM7ElTQ7SJsYVnrAjiiSl6R\neud1ggq6XYQtLGtdAEdUSStS9xkkQbXr0o4ZBNZbRYmVdIO1wYAoimTPYSKpB98C662ixEq6\nwdpgQBRFsucokfSnYtd1eZN1g7XBgCiKZM9RPRJFKqCS1gVwRJW8Iu04taNIBVTSugCOqJJY\nJB+KPZIRlbQugCOqUKRiOLVboZLWBXBEFYqEgCFZQQvLWhfAEVVSiuTqjgQVdLsIW1jWugCO\nqJJRJN+8TlBBt4uwhWWtC+CIKglFcp5BElTQ7SJsYVnrAjiiCkUChCJZYUAURbKHIllhQFTS\nugCOqJJQJPZIehgQlbQugCOqZBSJUzs1DIhKWhfAEVVSioRABa0rbGFZ6wI4ogpFQsCQrKCF\nZa0L4IgqFAkBQ7KCFpa1LoAjqiQTCdAdCSrodhG2sKx1ARxRJZdIiHmdoIJuF2ELy1oXwBFV\nUokEOYM0JOp2EbawrHUBHFGFIiFgSFbQwrLWBXBEFYqEgCFZQQvLWhfAEVVSicQeyQoDopLW\nBXBEFb9Ih+bXo7NLYJhCcu2RgKigdYUtLGtdAEdUoUgIGJIVtLCsdQEcUYUiIWBIVtDCstYF\ncESVNCI9pww32C7CFpa1LoAjqmQR6TWvu8F2EbawrHUBHFEliUjAM0hDom4XYQvLWhfAEVUo\nEgKGZAUtLGtdAEdUoUgIGJIVtLCsdQEcUSWJSOyR+mBAVNK6AI6okkUkTu26YEBU0roAjqiS\nRqQXDIgKWlfYwrLWBXBEFYqEgCFZQQvLWhfAEVUoEgKGZAUtLGtdAEdUoUgIGJIVtLCsdQEc\nUYUiIWBIVtDCstYFcEQVioSAIVlBC8taF8ARVSgSAoZkBS0sa10AR1ShSAgYkhW0sKx1ARxR\nhSIhYEhW0MKy1gVwRBWKhIAhWUELy1oXwBFVKBIChmQFLSxrXQBHVKFICBiSFbSwrHUBHFGF\nIiFgSFbQwrLWBXBEFYqEgCFZQQvLWhfAEVUoEgKGZAUtLGtdAEdUoUgIGJIVtLCsdQEcUYUi\nIWBIVtDCstYFcEQVioSAIVlBC8taF8ARVSgSAoZkBS0sa10AR1ShSAgYkhW0sKx1ARxRhSIh\nYEhW0MKy1gVwRBWKhIAhWUELy1oXwBFVKBIChmQFLSxrXQBHVKFICBiSFbSwrHUBHFGFIiFg\nSFbQwrLWBXBEFYqEgCFZQQvLWhfAEVUoEgKGZAUtLGtdAEdUoUgIGJIVtLCsdQEcUYUiIWBI\nVtDCstYFcESVvUX68+eP8qviQT+eJ5t+vy0UyQoDou4h0p9nVF81D/rxPNny+60vpvUJLVbS\nDdYGA6JuIdKfpLG+mMbHN1lJN1gbDIi6jEit/Pnz779n5aHjKgAAAlFJREFUS1FOu65dV0o7\n//574i9vhHW1s69IcZZzmah1hS2MdbVDkYIlamGsqx2KFCxRC2Nd7ewtEsPcIhSJYQChSAwD\nCEViGEAoEsMAsrNIj8dj31/Qk6GoxyNacWNFwep6SF3RVthQi9R0emn7ivQYFzdSXttE0MLk\nW8jaghUlAg11nV/eriI93l8iJcaKL+Tx+RqtuHj/5ZGd9+uPEdbZ/UQaX4BgByo/Y0UxV9pP\nvBVGkc6PbLCxNgypKOJKex9FRaqLIp2e9/YQr7YQG8Uqj9UfAoQinZ1H4U9REmKjWIUibeaG\nIkXt6eNsFMs85l+C5EYihRyOTQ7rQh3yT+fyser6/Kcn1Aq70fg7wHmydabnF8+uZZ5HlJOL\ny7xPfYaq604nZBnmJqFIDAMIRWIYQCgSwwBCkRgGEIrEMIBQJIYBhCIxDCAUiWEAoUgMAwhF\nSpAvvkrhw5coQShS/PAlShCKFD98ic7L4Mfz6/D/L/nXz7+/vn59jf/ExA1foPMyE2m05Wvy\nh+ELRcoQvkDnZSrS189n3yRfF39lQocv0Yl5qzI/yPt8fR/tnVUhow1fohNTFGlyJDcKRJHi\nhy/RianukSY/n35n4oYv0Zn5khegLBIP7RKFL9GZWY26J2MHDhtShS/RqZntcibj71GpySSc\nCR2+QqdmLVJhL8VdUoLwFTo1NOQq4Qt5aijSVcIX8sSw97lO+EqeGIp0nfCVZBhAKBLDAEKR\nGAYQisQwgFAkhgGEIjEMIBSJYQChSAwDCEViGEAoEsMAQpEYBpD/B1UoxZWks8nUAAAAAElF\nTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "colnames(probs) <- ycurrs\n",
    "probs %>%\n",
    "    tidyr::gather(ycurr, prob, factor_key = TRUE) %>%\n",
    "    ggplot(aes(rep(0:100, length(ycurrs)),\n",
    "               log(prob),\n",
    "               col = ycurr)) +\n",
    "    geom_point() +\n",
    "    xlab(\"ynext\") +\n",
    "    ylab(\"log(prob)\")"
   ]
  }
 ],
 "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.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}