{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Case Study: Predicting salmon returns\n",
    "\n",
    "This case study is based on a ModSim student project by Josh Deng and Erika Lu.\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Can we predict salmon populations?\n",
    "\n",
    "Each year the [U.S. Atlantic Salmon Assessment Committee](https://www.nefsc.noaa.gov/USASAC/Reports/USASAC2018-Report-30-2017-Activities.pdf) reports estimates of salmon populations in oceans and rivers in the northeastern United States.  The reports are useful for monitoring changes in these populations, but they generally do not include predictions.\n",
    "\n",
    "The goal of this case study is to model year-to-year changes in population, evaluate how predictable these changes are, and estimate the probability that a particular population will increase or decrease in the next 10 years.\n",
    "\n",
    "As an example, I'll use data from page 18 of the 2017 report, which provides population estimates for the Narraguagus and Sheepscot Rivers in Maine.\n",
    "\n",
    "![USASAC_Report_2017_Page18](data/USASAC_Report_2017_Page18.png)\n",
    "\n",
    "At the end of this notebook, I make some suggestions for extracting data from a PDF document automatically, but for this example I will keep it simple and type it in.\n",
    "\n",
    "Here are the population estimates for the Narraguagus River:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "pops = [2749, 2845, 4247, 1843, 2562, 1774, 1201, 1284, 1287, 2339, 1177, 962, 1176, 2149, 1404, 969, 1237, 1615, 1201];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get this data into a Pandas Series, I'll also make a range of years to use as an index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "years = range(1997, 2016)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop_series = TimeSeries(pops, index=years, dtype=np.float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what it looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_population(series):\n",
    "    plot(series, label='Estimated population')\n",
    "    decorate(xlabel='Year', \n",
    "             ylabel='Population estimate', \n",
    "             title='Narraguacus River',\n",
    "             ylim=[0, 5000])\n",
    "    \n",
    "plot_population(pop_series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modeling changes\n",
    "\n",
    "To see how the population changes from year-to-year, I'll use `ediff1d` to compute the absolute difference between each year and the next.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "abs_diffs = np.ediff1d(pop_series, to_end=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute relative differences by dividing by the original series elementwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "rel_diffs = abs_diffs / pop_series"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can use the `modsim` function `compute_rel_diff`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "rel_diffs = compute_rel_diff(pop_series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These relative differences are observed annual net growth rates.  So let's drop the `0` and save them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "rates = rel_diffs.drop(2015)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple way to model this system is to draw a random value from this series of observed rates each year.  We can use the NumPy function `choice` to make a random choice from a series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.choice(rates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation\n",
    "\n",
    "Now we can simulate the system by drawing random growth rates from the series of observed rates.\n",
    "\n",
    "I'll start the simulation in 2015."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_0 = 2015\n",
    "p_0 = pop_series[t_0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a `System` object with variables `t_0`, `p_0`, `rates`, and `duration=10` years. \n",
    "\n",
    "The series of observed rates is one big parameter of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = System(t_0=t_0,\n",
    "                p_0=p_0,\n",
    "                duration=10,\n",
    "                rates=rates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write an update functon that takes as parameters `pop`, `t`, and `system`.\n",
    "It should choose a random growth rate, compute the change in population, and return the new population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your update function and run it a few times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "update_func1(p_0, t_0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a version of `run_simulation` that stores the results in a `TimeSeries` and returns it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Simulate a queueing system.\n",
    "    \n",
    "    system: System object\n",
    "    update_func: function object\n",
    "    \"\"\"\n",
    "    t_0 = system.t_0\n",
    "    t_end = t_0 + system.duration\n",
    "    \n",
    "    results = TimeSeries()\n",
    "    results[t_0] = system.p_0\n",
    "    \n",
    "    for t in linrange(t_0, t_end):\n",
    "        results[t+1] = update_func(results[t], t, system)\n",
    "\n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `run_simulation` to run generate a prediction for the next 10 years.\n",
    "\n",
    "The plot your prediction along with the original data.  Your prediction should pick up where the data leave off."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get a sense of how much the results vary, we can run the model several times and plot all of the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_many_simulations(system, update_func, iters):\n",
    "    \"\"\"Runs simulations and plots the results.\n",
    "    \n",
    "    system: System object\n",
    "    update_func: function object\n",
    "    iters: number of simulations to run\n",
    "    \"\"\"\n",
    "    for i in range(iters):\n",
    "        results = run_simulation(system, update_func)\n",
    "        plot(results, color='gray', linewidth=5, alpha=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot option `alpha=0.1` makes the lines semi-transparent, so they are darker where they overlap.\n",
    "\n",
    "Run `plot_many_simulations` with your update function and `iters=30`.  Also plot the original data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are highly variable: according to this model, the population might continue to decline over the next 10 years, or it might recover and grow rapidly!\n",
    "\n",
    "It's hard to say how seriously we should take this model.  There are many factors that influence salmon populations that are not included in the model.  For example, if the population starts to grow quickly, it might be limited by resource limits, predators, or fishing.  If the population starts to fall, humans might restrict fishing and stock the river with farmed fish.\n",
    "\n",
    "So these results should probably not be considered useful predictions.  However, there might be something useful we can do, which is to estimate the probability that the population will increase or decrease in the next 10 years.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distribution of net changes\n",
    "\n",
    "To describe the distribution of net changes, write a function called `run_many_simulations` that runs many simulations, saves the final populations in a `ModSimSeries`, and returns the `ModSimSeries`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_many_simulations(system, update_func, iters):\n",
    "    \"\"\"Runs simulations and report final populations.\n",
    "    \n",
    "    system: System object\n",
    "    update_func: function object\n",
    "    iters: number of simulations to run\n",
    "    \n",
    "    returns: series of final populations\n",
    "    \"\"\"\n",
    "    # FILL THIS IN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your function by running it with `iters=5`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_many_simulations(system, update_func1, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run 1000 simulations and describe the distribution of the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "last_pops = run_many_simulations(system, update_func1, 1000)\n",
    "last_pops.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we substract off the initial population, we get the distribution of changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "net_changes = last_pops - p_0\n",
    "net_changes.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The median is negative, which indicates that the population decreases more often than it increases.\n",
    "\n",
    "We can be more specific by counting the number of runs where `net_changes` is positive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(net_changes > 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can use `mean` to compute the fraction of runs where `net_changes` is positive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(net_changes > 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the fraction where it's negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(net_changes < 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, based on observed past changes, this model predicts that the population is more likely to decrease than increase over the next 10 years, by about 2:1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A refined model\n",
    "\n",
    "There are a few ways we could improve the model.\n",
    "\n",
    "1.  It looks like there might be cyclic behavior in the past data, with a period of 4-5 years.  We could extend the model to include this effect.\n",
    "\n",
    "2.  Older data might not be as relevant for prediction as newer data, so we could give more weight to newer data.\n",
    "\n",
    "The second option is easier to implement, so let's try it.\n",
    "\n",
    "I'll use `linspace` to create an array of \"weights\" for the observed rates.  The probability that I choose each rate will be proportional to these weights.\n",
    "\n",
    "The weights have to add up to 1, so I divide through by the total."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = linspace(0, 1, len(rates))\n",
    "weights /= sum(weights)\n",
    "plot(weights)\n",
    "decorate(xlabel='Index into the rates array',\n",
    "         ylabel='Weight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll add the weights to the `System` object, since they are parameters of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "system.weights = weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can pass these weights as a parameter to `np.random.choice` (see the [documentation](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.choice.html))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.choice(system.rates, p=system.weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write an update function that takes the weights into account."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `plot_many_simulations` to plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `run_many_simulations` to collect the results and `describe` to summarize the distribution of net changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does the refined model have much effect on the probability of population decline?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting data from a PDF document\n",
    "\n",
    "The following section uses `tabula-py` to get data from a PDF document.\n",
    "\n",
    "If you don't already have it installed, and you are using Anaconda, you can install it by running the following command in a Terminal or Git Bash:\n",
    "\n",
    "```\n",
    "conda install -c conda-forge tabula-py\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tabula import read_pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = read_pdf('data/USASAC2018-Report-30-2017-Activities-Page11.pdf')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}