From aed6bfb9f76925e0e5270f3cd6a8b3a7db40b57a Mon Sep 17 00:00:00 2001
From: Joseph Gibson <joseph.gibson@ukhsa.gov.uk>
Date: Mon, 2 Sep 2024 15:50:49 +0100
Subject: [PATCH] add sympy bug notebook

---
 notebooks/sympy_bug_hunt.ipynb | 214 +++++++++++++++++++++++++++++++++
 1 file changed, 214 insertions(+)
 create mode 100644 notebooks/sympy_bug_hunt.ipynb

diff --git a/notebooks/sympy_bug_hunt.ipynb b/notebooks/sympy_bug_hunt.ipynb
new file mode 100644
index 0000000..568e08c
--- /dev/null
+++ b/notebooks/sympy_bug_hunt.ipynb
@@ -0,0 +1,214 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Debug adaptive tau leap\n",
+    "\n",
+    "First, we set up SIR model with time dependent infectivity"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pygom import SimulateOde, Transition, TransitionType\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import random\n",
+    "import math\n",
+    "import sympy\n",
+    "\n",
+    "###################\n",
+    "# ODE specification\n",
+    "###################\n",
+    "\n",
+    "# Define SIR model\n",
+    "stateList = ['S', 'I', 'R']\n",
+    "paramList = ['beta', 'gamma', 'N']\n",
+    "transitionList = [Transition(origin='S', destination='I', equation='beta*(1+0.99*cos(2*3.14*t/(10)))*S*I/N', transition_type=TransitionType.T),\n",
+    "                  Transition(origin='I', destination='R', equation='gamma*I', transition_type=TransitionType.T)]\n",
+    "\n",
+    "n_pop=1e4  # Total population is fixed\n",
+    "\n",
+    "beta=0.4\n",
+    "gamma=1/4\n",
+    "\n",
+    "params=[('beta', beta), ('gamma', gamma), ('N', n_pop)]\n",
+    "\n",
+    "# Initial conditions\n",
+    "i0=10\n",
+    "x0 = [n_pop-i0, i0, 0]\n",
+    "\n",
+    "\n",
+    "# Set up pygom object (D_F suffix implies Deterministic_Fixed)\n",
+    "model = SimulateOde(stateList, paramList, transition=transitionList)\n",
+    "model.initial_values = (x0, np.float64(0))  # (initial state conditions, initial timepoint)\n",
+    "model.parameters=params"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As part of the class initialisation, the string equations get turned into sympy expressions. For example:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 59,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "I*S*beta*(0.99*cos(0.628*t) + 1)/N\n",
+      "<class 'sympy.core.mul.Mul'>\n"
+     ]
+    }
+   ],
+   "source": [
+    "eqn=model._transitionVector[0]\n",
+    "print(eqn)\n",
+    "print(type(eqn))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can see which variables/states are present in this equation:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 60,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "{N, t, I, S, beta}\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(eqn.free_symbols)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This promisingly includes time, t.\n",
+    "If we want to perform operations with the variables, such as derivatives, we can't define the variable externally"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 67,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "beta\n",
+      "False\n",
+      "0\n"
+     ]
+    }
+   ],
+   "source": [
+    "state=sympy.symbols(\"beta\")\n",
+    "print(state)\n",
+    "print(state in eqn.free_symbols)\n",
+    "print(sympy.diff(eqn, state))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Instead, we have to take them from where they were defined in the class:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 68,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "beta\n",
+      "True\n",
+      "I*S*(0.99*cos(0.628*t) + 1)/N\n"
+     ]
+    }
+   ],
+   "source": [
+    "state=model._paramDict['beta']\n",
+    "print(state)\n",
+    "print(state in eqn.free_symbols)\n",
+    "print(sympy.diff(eqn, state))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Trouble is, when we pull time from the same dictionary, it is not recognised like the other variables"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 69,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "t\n",
+      "False\n",
+      "0\n"
+     ]
+    }
+   ],
+   "source": [
+    "state=model._paramDict['t']\n",
+    "print(state)\n",
+    "print(state in eqn.free_symbols)\n",
+    "print(sympy.diff(eqn, state))"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "pygom_19",
+   "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.11.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}