Merge pull request #127 from ukhsa-collaboration/dev

Merge dev to master for new version 0.1.9

docs/_toc.yml

@@ -8,9 +8,13 @@ parts:
     chapters:
     - file: md/installation
     - file: md/building_doc
+  - caption: Compartmental models
+    chapters:
+    - file: notebooks/definition
+    - file: notebooks/solution
   - caption: PyGOM workflow
     chapters:
-    - file: notebooks/model_spec
+    - file: notebooks/model_spec_2
     - file: notebooks/insights
       sections:
         - file: notebooks/extract_info

docs/md/intro.md

@@ -6,16 +6,35 @@ PyGOM (Python Generic ODE Model) is a Python package which provides a simple int
 This is backed by a comprehensive and easy to use tool–box implementing functions to easily perform common operations such as parameter estimation and solving for deterministic or stochastic time evolution.
 The package source is freely available (hosted on [GitHub](https://github.com/ukhsa-collaboration/pygom)) and organized in a way that permits easy extension. With both the algebraic and numeric calculations performed automatically (but still accessible), the end user is freed to focus on model development.
 
-## What is new in this release?
-
-The main objective of the current release (0.1.8) is to provide more comprehensive documentation on how to use PyGOM.
-
-The code underlying PyGOM's functionality is largely unchanged since the previous release, barring a few minor bug fixes.
-The only significant changes which previous users should be aware of are:
-- A move away from the {class}`DeterministicOde` class for deterministic simulations and instead employing {class}`SimulateOde` as our do-all class for deterministic or stochastic simulations as well as parameter fitting.
-- Running simulations with random parameters does not require a special simulation function. Instead, PyGOM now recognises the parameter types handed to it (fixed or random) and acts accordingly. This means that stochastic simulations can now be performed with random parameters.
-
-Both these changes are outlined in more detail in the {doc}`Producing forecasts <solving>` section.
+## Release notes:
+
+### [0.1.9] - 2024-11-30 (Latest release)
+
+- Added
+    - Method {func}`solve_stochast` of {class}`SimulateOde` has 1 additional output (for a total of 3): the number of times each event occurs in between each time step. This is useful if you are interested in infection incidence rather than prevalence, for example.
+    - The adaptive tau leap can be bypassed and, instead, a constant step size for the tau leap algorithm can be specified with `SimulateOde._pre_tau`.
+    - When an output of a stochastic simulation needs to be mapped to the user defined time-points, this is performed via linear interpolation.
+- Changed
+    - The {class}`Event` class has replaced the {class}`Transition` class as the fundamental building block. This allows for more flexibility when dealing with correlated transitions, for example.
+    - Internal methods {func}`add_func` and {func}`add_compiled_sympy_object` make the compilation of sympy objects more modular.
+- Deprecated
+    - Models can still be defined via {class}`Transition` objects as before, but users are advised to switch to the {class}`Event` based approach. 
+    - The birth rate state can still be defined as an `origin`, but can now also (and more accurately) be described as a `destination`. In the next version, this may remove the need to specify the transition type since between state transitions will uniquely have an origin and destination, births will uniquely have just a destination and deaths just an origin.
+- Removed
+- Fixed
+    - Minor bug fixes
+
+### [0.1.8] - 2024-08-06
+
+- Added
+    - Comprehensive documentation of how to use PyGOM
+- Changed
+    - Running simulations with random parameters does not require a special simulation function. Instead, PyGOM now recognises the parameter types handed to it (fixed or random) and acts accordingly. This means that stochastic simulations can now be performed with random parameters.
+- Deprecated
+    - {class}`DeterministicOde` is deprecated with {class}`SimulateOde` now performing both deterministic and stochastic simulations.
+- Removed
+- Fixed
+    - Minor bug fixes
 
 ## Using this documentation
 This documentation is built using [JupyterBook](https://jupyterbook.org/en/stable/intro.html).

docs/notebooks/Digraph.gv

@@ -0,0 +1,8 @@
+digraph {
+rankdir=LR	S
+	I
+	R
+	size="5,5"
+	S -> I [label="βSI/N"]
+	I -> R [label="γI"]
+}

docs/notebooks/common_models/Lotka_Volterra.ipynb

@@ -30,22 +30,30 @@
     "import matplotlib.pyplot as plt\n",
     "import math\n",
     "\n",
-    "# population density of predators and prey per square m\n",
-    "x0 = [1, 0.5]\n",
+    "# Initial populations\n",
+    "x0=[40, 20]\n",
     "\n",
-    "# total area we wish to consider\n",
-    "area=200\n",
+    "# Parameters from a snowshoe hare / Canadian lynx system\n",
+    "# https://mc-stan.org/users/documentation/case-studies/lotka-volterra-predator-prey.html\n",
+    "alpha=0.545\n",
+    "beta=0.028\n",
+    "gamma=0.803\n",
+    "delta=0.024\n",
     "\n",
-    "# total animal populations\n",
-    "x0 = [x * area for x in x0]\n",
+    "# scale up the population (this will need scaling in the predation parameters)\n",
+    "scale=10\n",
     "\n",
-    "ode = common_models.Lotka_Volterra({'alpha':0.1,\n",
-    "                                    'beta':0.2,\n",
-    "                                    'gamma':0.3,\n",
-    "                                    'delta':0.25,\n",
-    "                                    'A':area})\n",
+    "x0 = [x * scale for x in x0]\n",
+    "beta=beta/scale\n",
+    "delta=delta/scale\n",
     "\n",
-    "tmax=200                              # maximum time over which to run solver\n",
+    "\n",
+    "ode = common_models.Lotka_Volterra({'alpha':alpha,\n",
+    "                                    'beta':beta,\n",
+    "                                    'gamma':gamma,\n",
+    "                                    'delta':delta})\n",
+    "\n",
+    "tmax=50                               # maximum time over which to run solver\n",
     "dt=0.1                                # timestep\n",
     "n_timestep=math.ceil(tmax/dt)         # number of iterations\n",
     "t = np.linspace(0, tmax, n_timestep)  # times at which solution will be evaluated\n",
@@ -80,7 +88,7 @@
     "ax.set_ylabel(\"Population number\")\n",
     "ax.plot(t, solution[:,0], label=\"prey\")\n",
     "ax.plot(t, solution[:,1], label=\"predator\")\n",
-    "ax.legend(loc=\"upper right\")\n",
+    "ax.legend()\n",
     "plt.show()"
    ]
   },
@@ -102,7 +110,7 @@
     "np.random.seed(1)\n",
     "\n",
     "n_sim=1\n",
-    "solution, simT = ode.solve_stochast(t, n_sim, full_output=True)\n",
+    "solution, jump, simT = ode.solve_stochast(t, n_sim, full_output=True)\n",
     "\n",
     "f, ax = plt.subplots(figsize=(10, 2))\n",
     "\n",
@@ -112,7 +120,7 @@
     "ax.set_ylabel(\"Population number\")\n",
     "ax.plot(t, y[:,0], label=\"prey\")\n",
     "ax.plot(t, y[:,1], label=\"predator\")\n",
-    "ax.legend(loc=\"upper right\")\n",
+    "ax.legend()\n",
     "plt.show()"
    ]
   },
@@ -142,7 +150,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.9.15 ('sphinx-doc')",
+   "display_name": "pygom",
    "language": "python",
    "name": "python3"
   },
@@ -156,12 +164,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.19"
-  },
-  "vscode": {
-   "interpreter": {
-    "hash": "4dc1e323c80fe09539c74ad5c5a7c7d8d9ff99e04f7b3dbd3680daf878629d6e"
-   }
+   "version": "3.11.10"
   }
  },
  "nbformat": 4,

docs/notebooks/common_models/SIR.ipynb

@@ -8,12 +8,28 @@
     "\n",
     "{func}`.SIR`\n",
     "\n",
-    "The standard Susceptible-Infected-Recovered (SIR) model, which features heavily throughout this documentation, is defined by the following equations:\n",
+    "The standard Susceptible-Infected-Recovered (SIR) model features heavily throughout this documentation and so we outline its key features for anyone unfamiliar.\n",
+    "\n",
+    "```{warning}\n",
+    "There is some ambiguity in the naming of the Recovered class which is also commonly referred to as Removed.\n",
+    "In the latter sense, there is no distinction made between those who can no longer be infected due to infection acquired immunity and infection induced death.\n",
+    "However, for more complex models, the recovered class may become susceptible again due to the effects of immune waning and deaths versus recoveries are typically important to distinguish in real world applications.\n",
+    "Therefore, in this tutorial, the Recovered class will be reserved solely for those who survive infection.\n",
+    "```\n",
+    "\n",
+    "The assumptions of the SIR model can be stated as:\n",
+    "1) An average member of the population makes contact sufficient to transmit or receive infection with $c$ others per unit time. Each of these events carries a probability, $p$, of transmission such that each individual has an average of $cp = \\beta$ infectious contacts per unit time. This fixed contact rate reflects what is referred to as *standard* incidence, as opposed to *mass-action* incidence, where contacts per person are proportional to the total population size, $N$.\n",
+    "2) The population interacts heterogeneously as if a well mixed continuum.\n",
+    "For instance, a susceptible does not have contacts with other individuals, but with the entire population on average.\n",
+    "3) The infective class recovers at a rate, $\\gamma$.\n",
+    "4) No births, deaths (natural or from disease) or migration mean there is no entry into or departure from the population: $S(t)+I(t)+R(t)=N$.\n",
+    "\n",
+    "Under these assumptions, the rates of change of the population in each compartment (**S**usceptible, **I**nfected and **R**ecovered) are given by the following equations:\n",
     "\n",
     "$$\\begin{aligned}\n",
-    "\\frac{\\mathrm{d}S}{\\mathrm{d}t} &= - \\frac{\\beta SI}{N} \\\\\n",
-    "\\frac{\\mathrm{d}I}{\\mathrm{d}t} &= \\frac{\\beta SI}{N} - \\gamma I \\\\\n",
-    "\\frac{\\mathrm{d}R}{\\mathrm{d}t} &= \\gamma I\n",
+    "\\frac{\\mathrm{d} S}{\\mathrm{d} t} &= -\\frac{\\beta SI}{N} \\\\\n",
+    "\\frac{\\mathrm{d} I}{\\mathrm{d} t} &= \\frac{\\beta SI}{N} - \\gamma I \\\\\n",
+    "\\frac{\\mathrm{d} R}{\\mathrm{d} t} &= \\gamma I.\n",
     "\\end{aligned}$$\n",
     "\n",
     "We solve deterministically for flu-like parameters:"

docs/notebooks/common_models/SIR_Birth_Death.ipynb

@@ -26,7 +26,11 @@
    "cell_type": "code",
    "execution_count": null,
    "id": "1e5a8099",
-   "metadata": {},
+   "metadata": {
+    "tags": [
+     "hide-output"
+    ]
+   },
    "outputs": [],
    "source": [
     "from pygom import common_models\n",
@@ -45,7 +49,7 @@
     "R0=15\n",
     "beta=R0*gamma\n",
     "\n",
-    "ode = common_models.SIR_Birth_Death({'beta':beta, 'gamma':gamma, 'mu':mu, 'N':n_pop})\n",
+    "ode = common_models.SIR_Birth_Death({'beta':beta, 'gamma':gamma, 'mu':mu})\n",
     "\n",
     "# Time range and increments\n",
     "tmax=365*10                           # maximum time over which to run solver\n",
@@ -58,7 +62,7 @@
     "    s0=math.floor((gamma+mu)*n_pop/beta)\n",
     "    i0=math.floor(mu*(n_pop-s0)*n_pop/(beta*s0))\n",
     "    r0=n_pop-(s0+i0)\n",
-    "    return [s0, i0, r0]\n",
+    "    return [s0, i0, r0, n_pop]\n",
     "\n",
     "x0=sir_bd_endemic_eq(mu, beta, gamma, n_pop)\n",
     "ode.initial_values = (x0, t[0])\n",
@@ -69,27 +73,27 @@
     "n_sim=10\n",
     "np.random.seed(1)\n",
     "\n",
-    "solution, simT = ode.solve_stochast(t, n_sim, full_output=True)\n",
+    "solution, jumps, simT = ode.solve_stochast(t, n_sim, full_output=True)\n",
     "y=np.dstack(solution)\n",
     "\n",
     "############################\n",
     "# try larger population size\n",
     "############################\n",
     "n_pop=1e5\n",
-    "ode = common_models.SIR_Birth_Death({'beta':beta, 'gamma':gamma, 'mu':mu, 'N':n_pop})  # update parameter\n",
+    "ode = common_models.SIR_Birth_Death({'beta':beta, 'gamma':gamma, 'mu':mu})  # update parameter\n",
     "x0=sir_bd_endemic_eq(mu, beta, gamma, n_pop)  # recalculate IC's\n",
     "ode.initial_values = (x0, t[0])\n",
-    "solution_2, simT_2 = ode.solve_stochast(t, n_sim, full_output=True)  # simulate\n",
+    "solution_2, jumps_2, simT_2 = ode.solve_stochast(t, n_sim, full_output=True)  # simulate\n",
     "y_2=np.dstack(solution_2)\n",
     "\n",
     "#################################\n",
     "# try even larger population size\n",
     "#################################\n",
     "n_pop=1e6\n",
-    "ode = common_models.SIR_Birth_Death({'beta':beta, 'gamma':gamma, 'mu':mu, 'N':n_pop})  # update parameter\n",
+    "ode = common_models.SIR_Birth_Death({'beta':beta, 'gamma':gamma, 'mu':mu})  # update parameter\n",
     "x0=sir_bd_endemic_eq(mu, beta, gamma, n_pop)  # recalculate IC's\n",
     "ode.initial_values = (x0, t[0])\n",
-    "solution_3, simT_3 = ode.solve_stochast(t, n_sim, full_output=True)  # simulate\n",
+    "solution_3, jumps_3, simT_3 = ode.solve_stochast(t, n_sim, full_output=True)  # simulate\n",
     "y_3=np.dstack(solution_3)"
    ]
   },
@@ -143,7 +147,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "pygom_development",
+   "display_name": "pygom",
    "language": "python",
    "name": "python3"
   },
@@ -157,7 +161,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.19"
+   "version": "3.11.10"
   }
  },
  "nbformat": 4,

docs/notebooks/common_models/SIS_Periodic.ipynb

@@ -58,7 +58,8 @@
     "\n",
     "# Initial conditions\n",
     "i0=0.1*n_pop\n",
-    "x0 = [i0, t[0]]\n",
+    "s0=n_pop-i0\n",
+    "x0 = [s0, i0]\n",
     "\n",
     "ode.initial_values = (x0, t[0])\n",
     "solution=ode.solve_determ(t[1::])"
@@ -86,15 +87,15 @@
     "f, ax = plt.subplots(figsize=(10, 4))\n",
     "\n",
     "ax.set_xlabel(\"Time\")\n",
-    "ax.set_ylabel(\"I\", rotation=0)\n",
-    "ax.plot(t, solution[:,0])\n",
+    "ax.set_title(\"Infected\")\n",
+    "ax.plot(t, solution[:,1])\n",
     "plt.show()"
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.9.15 ('sphinx-doc')",
+   "display_name": "pygom",
    "language": "python",
    "name": "python3"
   },
@@ -108,12 +109,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.19"
-  },
-  "vscode": {
-   "interpreter": {
-    "hash": "4dc1e323c80fe09539c74ad5c5a7c7d8d9ff99e04f7b3dbd3680daf878629d6e"
-   }
+   "version": "3.11.10"
   }
  },
  "nbformat": 4,