Dev

Data Generation.ipynb

@@ -0,0 +1,191 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "dea6bfd2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt\n",
+    "from IPython.display import clear_output\n",
+    "from scipy.integrate import solve_ivp"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "12e528e1",
+   "metadata": {},
+   "source": [
+    "# SIR Model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "889995b4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "'''\n",
+    "Showcase Samples:\n",
+    "1. (S_0, I_0, R_0) = (9900, 100, 0) + (b, g) = (1.5, 0.75).\n",
+    "2. (S_0, I_0, R_0) = (9900, 100, 0) + (b, g) = (3.5, 1.75).\n",
+    "\n",
+    "Standard Settings:\n",
+    "1. 20x noisy observations per unit time with discretization 1.\n",
+    "2. {0.5, 2, 8} interval lengths - we'll generate datasets to 20.00.\n",
+    "3. {0.05, 0.15} alpha noise level multipliers.\n",
+    "4. 10x randomized seeds.\n",
+    "'''\n",
+    "# generate separate datasets for each possible combination, using maximum precision\n",
+    "for seed in range(10):\n",
+    "    for alpha in [0.05, 0.15]:\n",
+    "        \n",
+    "        # two possible combinations of params (beta, gamma)\n",
+    "        for params in [(1.5, 0.75), (3.5, 1.75)]:\n",
+    "    \n",
+    "            # set a seed for reproducibility\n",
+    "            np.random.seed(seed)\n",
+    "\n",
+    "            # what are our parameters + population size?\n",
+    "            b, g, N = params[0], params[1], 10000\n",
+    "\n",
+    "            # encode ODEs for solve_ivp data-generation processes\n",
+    "            def SIR(t, y):\n",
+    "\n",
+    "                # unpack y\n",
+    "                S, I, R = tuple(y)\n",
+    "\n",
+    "                # dS/dt = -b*S*I/N; dI/dt = b*I*S/N - g*I; dR/dt = g*I\n",
+    "                dSdt = -b*S*I/N\n",
+    "                dIdt = +b*S*I/N - g*I\n",
+    "                dRdt = +g*I\n",
+    "\n",
+    "                # return only the derivatives\n",
+    "                return np.array([dSdt, dIdt, dRdt])\n",
+    "    \n",
+    "            # generate our data\n",
+    "            t_start, t_end = 0.0, 20.0\n",
+    "            t_steps = np.linspace(start=t_start, stop=t_end, num=20001)\n",
+    "            SIR_init = np.array([9900, 100, 0])\n",
+    "            X = solve_ivp(fun=SIR, t_span=(t_start, t_end), y0=SIR_init, \n",
+    "                          t_eval=t_steps, atol=1e-10, rtol=1e-10).y.T\n",
+    "\n",
+    "            # compute appropriate noise levels based on alpha choice\n",
+    "            sigmas = alpha * (X.max(axis=0) - X.min(axis=0))\n",
+    "            \n",
+    "            # create a copy of X + add the noise\n",
+    "            X_noised = X.copy()\n",
+    "            X_noised += np.random.normal(loc=0.0, scale=sigmas, size=X.shape)\n",
+    "            \n",
+    "            # save time, X_true, X_noised as a .csv for our gold-standard dataset\n",
+    "            data = np.hstack([t_steps.reshape(-1, 1), X_noised, X])\n",
+    "            cols = [\"t\", \"S_obs\", \"I_obs\", \"R_obs\", \"S_true\", \"I_true\", \"R_true\"]\n",
+    "            \n",
+    "            # save this dataset to our gold-standard datasets list\n",
+    "            df = pd.DataFrame(data=data, columns=cols)\n",
+    "            df.to_csv(f\"data/SIR_beta={b}_gamma={g}_alpha={alpha}_seed={seed}.csv\", index=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ecf5e40f",
+   "metadata": {},
+   "source": [
+    "# Lorenz Model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "9de8613c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "'''\n",
+    "Showcase Samples:\n",
+    "1. (X_0, Y_0, Z_0) = (5, 5, 5) + (beta, rho, sigma) = (8/3, 28.0, 10.0).\n",
+    "2. (X_0, Y_0, Z_0) = (5, 5, 5) + (beta, rho, sigma) = (8/3, 6.0, 10.0).\n",
+    "\n",
+    "Standard Settings:\n",
+    "1. 20x noisy observations per unit time with discretization 1.\n",
+    "2. {0.5, 2, 8} interval lengths - we'll generate datasets to 20.00.\n",
+    "3. {0.05, 0.15} alpha noise level multipliers.\n",
+    "4. 10x randomized seeds.\n",
+    "'''\n",
+    "# generate separate datasets for each possible combination, using maximum precision\n",
+    "for seed in range(10):\n",
+    "    for alpha in [0.05, 0.15]:\n",
+    "        \n",
+    "        # two possible combinations of params (beta, gamma)\n",
+    "        for params in [(8/3, 28.0, 10.0), (8/3, 6.0, 10.0)]:\n",
+    "            \n",
+    "             # set a seed for reproducibility\n",
+    "            np.random.seed(seed)\n",
+    "\n",
+    "            # what are our parameters?\n",
+    "            beta, rho, sigma = params\n",
+    "\n",
+    "            # encode ODEs for solve_ivp data-generation processes\n",
+    "            def lorenz(t, y):\n",
+    "\n",
+    "                # unpack y\n",
+    "                X, Y, Z = tuple(y)\n",
+    "\n",
+    "                # dXdt = sigma * (Y-X); dYdt = x(rho - z) - y; dZdt = xy - beta*z\n",
+    "                dXdt = sigma * (Y-X)\n",
+    "                dYdt = (X * (rho - Z)) - Y\n",
+    "                dZdt = (X*Y) - (beta*Z)\n",
+    "\n",
+    "                # return only the derivatives\n",
+    "                return np.array([dXdt, dYdt, dZdt])\n",
+    "            \n",
+    "            # generate our data\n",
+    "            t_start, t_end = 0.0, 20.0\n",
+    "            t_steps = np.linspace(start=t_start, stop=t_end, num=20001)\n",
+    "            y_init = np.array([5.0, 5.0, 5.0])\n",
+    "            X = solve_ivp(fun=lorenz, t_span=(t_start, t_end), y0=y_init, \n",
+    "                          t_eval=t_steps, atol=1e-10, rtol=1e-10).y.T\n",
+    "            \n",
+    "            # compute appropriate noise levels based on alpha choice\n",
+    "            sigmas = alpha * (X.max(axis=0) - X.min(axis=0))\n",
+    "            \n",
+    "            # create a copy of X + add the noise\n",
+    "            X_noised = X.copy()\n",
+    "            X_noised += np.random.normal(loc=0.0, scale=sigmas, size=X.shape)\n",
+    "            \n",
+    "            # save time, X_true, X_noised as a .csv for our gold-standard dataset\n",
+    "            data = np.hstack([t_steps.reshape(-1, 1), X_noised, X])\n",
+    "            cols = [\"t\", \"X_obs\", \"Y_obs\", \"Z_obs\", \"X_true\", \"Y_true\", \"Z_true\"]\n",
+    "            \n",
+    "            # save this dataset to our gold-standard datasets list\n",
+    "            df = pd.DataFrame(data=data, columns=cols)\n",
+    "            df.to_csv(f\"data/LORENZ_rho={rho}_alpha={alpha}_seed={seed}.csv\", index=False)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3.10 (Afterburner)",
+   "language": "python",
+   "name": "afterburner"
+  },
+  "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.10.14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}