{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 17\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": {
    "collapsed": true
   },
   "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": [
    "### Data\n",
    "\n",
    "We have data from Pacini and Bergman (1986), \"MINMOD: a computer program to calculate insulin sensitivity and pancreatic responsivity from the frequently sampled intravenous glucose tolerance test\", *Computer Methods and Programs in Biomedicine*, 23: 113-122.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('data/glucose_insulin.csv', index_col='time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the glucose time series looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(data.glucose, 'bo', label='glucose')\n",
    "decorate(xlabel='Time (min)',\n",
    "         ylabel='Concentration (mg/dL)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the insulin time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(data.insulin, 'go', label='insulin')\n",
    "decorate(xlabel='Time (min)',\n",
    "         ylabel='Concentration ($\\mu$U/mL)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the book, I put them in a single figure, using `subplot`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "subplot(2, 1, 1)\n",
    "plot(data.glucose, 'bo', label='glucose')\n",
    "decorate(ylabel='mg/dL')\n",
    "\n",
    "subplot(2, 1, 2)\n",
    "plot(data.insulin, 'go', label='insulin')\n",
    "decorate(xlabel='Time (min)',\n",
    "         ylabel='$\\mu$U/mL')\n",
    "\n",
    "savefig('figs/chap08-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolation\n",
    "\n",
    "We have measurements of insulin concentration at discrete points in time, but we need to estimate it at intervening points.  We'll use `interpolate`, which takes a `Series` and returns a function:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The return value from `interpolate` is a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "I = interpolate(data.insulin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the result, `I`, to estimate the insulin level at any point in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "I(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`I` can also take an array of time and return an array of estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_0 = get_first_label(data)\n",
    "t_end = get_last_label(data)\n",
    "ts = linrange(t_0, t_end, endpoint=True)\n",
    "I(ts)\n",
    "type(ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the interpolated values look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(data.insulin, 'go', label='insulin data')\n",
    "plot(ts, I(ts), color='green', label='interpolated')\n",
    "\n",
    "decorate(xlabel='Time (min)',\n",
    "         ylabel='Concentration ($\\mu$U/mL)')\n",
    "\n",
    "savefig('figs/chap08-fig02.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** [Read the documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html) of `scipy.interpolate.interp1d`.  Pass a keyword argument to `interpolate` to specify one of the other kinds of interpolation, and run the code again to see what it looks like. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The glucose minimal model\n",
    "\n",
    "I'll cheat by starting with parameters that fit the data roughly; then we'll see how to improve them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "params = Params(G0 = 290,\n",
    "                k1 = 0.03,\n",
    "                k2 = 0.02,\n",
    "                k3 = 1e-05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a version of `make_system` that takes the parameters and data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_system(params, data):\n",
    "    \"\"\"Makes a System object with the given parameters.\n",
    "    \n",
    "    params: sequence of G0, k1, k2, k3\n",
    "    data: DataFrame with `glucose` and `insulin`\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    G0, k1, k2, k3 = params\n",
    "    \n",
    "    Gb = data.glucose[0]\n",
    "    Ib = data.insulin[0]\n",
    "    \n",
    "    t_0 = get_first_label(data)\n",
    "    t_end = get_last_label(data)\n",
    "\n",
    "    init = State(G=G0, X=0)\n",
    "    \n",
    "    return System(G0=G0, k1=k1, k2=k2, k3=k3,\n",
    "                  init=init, Gb=Gb, Ib=Ib,\n",
    "                  t_0=t_0, t_end=t_end, dt=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the update function.  It uses `unpack` to make the system variables accessible without using dot notation, which makes the translation of the differential equations more readable and less error prone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def update_func(state, t, system):\n",
    "    \"\"\"Updates the glucose minimal model.\n",
    "    \n",
    "    state: State object\n",
    "    t: time in min\n",
    "    system: System object\n",
    "    \n",
    "    returns: State object\n",
    "    \"\"\"\n",
    "    G, X = state\n",
    "    unpack(system)\n",
    "        \n",
    "    dGdt = -k1 * (G - Gb) - X*G\n",
    "    dXdt = k3 * (I(t) - Ib) - k2 * X\n",
    "    \n",
    "    G += dGdt * dt\n",
    "    X += dXdt * dt\n",
    "\n",
    "    return State(G=G, X=X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before running the simulation, it is always a good idea to test the update function using the initial conditions.  In this case we can veryify that the results are at least qualitatively correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "update_func(system.init, system.t_0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now `run_simulation` is pretty much the same as it always is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "        \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    unpack(system)\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.row[t_0] = init\n",
    "    ts = linrange(t_0, t_end, dt)\n",
    "    \n",
    "    for t in ts:\n",
    "        frame.row[t+dt] = update_func(frame.row[t], t, system)\n",
    "    \n",
    "    return frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's how we run it.  `%time` is a Jupyter magic command that runs the function and reports its run time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time results = run_simulation(system, update_func);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are in a `TimeFrame object` with one column per state variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following plot shows the results of the simulation along with the actual glucose data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "subplot(2, 1, 1)\n",
    "\n",
    "plot(results.G, 'b-', label='simulation')\n",
    "plot(data.glucose, 'bo', label='glucose data')\n",
    "decorate(ylabel='mg/dL')\n",
    "\n",
    "subplot(2, 1, 2)\n",
    "\n",
    "plot(results.X, 'g-', label='remote insulin')\n",
    "\n",
    "decorate(xlabel='Time (min)', \n",
    "         ylabel='Arbitrary units')\n",
    "\n",
    "savefig('figs/chap08-fig03.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Under the hood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "%psource interpolate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "**Exercise:**  Our solution to the differential equations is only approximate because we used a finite step size, `dt=2` minutes.\n",
    "\n",
    "If we make the step size smaller, we expect the solution to be more accurate.  Run the simulation with `dt=1` and compare the results.  What is the largest relative error between the two solutions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}