Skip to content

Commit

Permalink
🎨 WIP formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
Christopher Lorton committed Feb 3, 2025
1 parent 0a8e6b2 commit 6faff0a
Show file tree
Hide file tree
Showing 18 changed files with 5,229 additions and 1,714 deletions.
1 change: 1 addition & 0 deletions notebooks/.gitignore
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
SIR_output.csv
163 changes: 104 additions & 59 deletions notebooks/2patch_SIR_wbirths_correlation.ipynb
View file

Large diffs are not rendered by default.

15 changes: 6 additions & 9 deletions notebooks/SEIR_wavefrontspeed.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,6 +142,7 @@
"import matplotlib.pyplot as plt\n",
"from matplotlib.figure import Figure\n",
"\n",
"\n",
"class Incubation:\n",
" def __init__(self, model, verbose: bool = False) -> None:\n",
" return\n",
Expand Down Expand Up @@ -183,23 +184,19 @@
"metadata": {},
"outputs": [],
"source": [
"from laser_generic import Model\n",
"from laser_generic import Susceptibility\n",
"from laser_generic import Infection\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"from laser_generic import Infection\n",
"from laser_generic import Model\n",
"from laser_generic import Susceptibility\n",
"\n",
"# Fox density from 1/km^2 to 4+/km^2 (higher in suburban and urban areas)\n",
"# Each patch will represent 5km^2\n",
"scenario = pd.DataFrame(data=np.full(201, 20, dtype=np.int32), columns=[\"population\"])\n",
"\n",
"model = Model(scenario)\n",
"model.components = [\n",
" Susceptibility,\n",
" Incubation,\n",
" Infection,\n",
" Transmission\n",
"]"
"model.components = [Susceptibility, Incubation, Infection, Transmission]"
]
},
{
Expand Down
223 changes: 106 additions & 117 deletions notebooks/SIR_CCS.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,22 +26,20 @@
}
],
"source": [
"from pathlib import Path\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"from laser_core.propertyset import PropertySet\n",
"import matplotlib.pyplot as plt\n",
"import os\n",
"from scipy.optimize import fsolve\n",
"\n",
"from laser_generic import Model\n",
"from laser_generic import Births_ConstantPop\n",
"from laser_generic import Infection\n",
"from laser_generic import Model\n",
"from laser_generic import Susceptibility\n",
"from laser_generic import Transmission\n",
"from laser_generic import Births_ConstantPop\n",
"from laser_generic.importation import Infect_Agents_In_Patch\n",
"\n",
"from laser_generic.utils import set_initial_susceptibility_in_patch\n",
"from laser_generic.utils import seed_infections_in_patch\n",
"\n",
"%load_ext line_profiler\n",
"\n",
Expand Down Expand Up @@ -71,56 +69,56 @@
"source": [
"%%capture\n",
"\n",
"nticks = 50*365\n",
"nticks = 50 * 365\n",
"npatches = 61\n",
"pops = np.logspace(3, 6, npatches)\n",
"scenario = pd.DataFrame({\n",
" \"name\": [str(i) for i in range(npatches)],\n",
" \"population\": pops\n",
"})\n",
"scenario = pd.DataFrame({\"name\": [str(i) for i in range(npatches)], \"population\": pops})\n",
"\n",
"#np.random.seed(5) # Ensure reproducibility\n",
"# np.random.seed(5) # Ensure reproducibility\n",
"nsims = 200\n",
"R0_samples = np.random.uniform(3, 16, nsims)\n",
"infmean_samples = 5+np.random.gamma(2, 10, nsims)\n",
"cbr_samples = 10+np.random.gamma(2, 20, nsims)\n",
"infmean_samples = 5 + np.random.gamma(2, 10, nsims)\n",
"cbr_samples = 10 + np.random.gamma(2, 20, nsims)\n",
"i = 0\n",
"outputs = np.zeros((nsims, nticks, npatches))\n",
"# Create a folder to store the outputs\n",
"output_folder = os.path.abspath(os.path.join(os.getcwd(), '..', '..', 'laser-generic-outputs', 'CCSSIRoutputs2'))\n",
"if not os.path.exists(output_folder):\n",
" os.makedirs(output_folder)\n",
"output_folder = (Path.cwd().parent.parent / \"laser-generic-outputs\" / \"CCSSIRoutputs2\").absolute()\n",
"if not output_folder.exists():\n",
" Path.mkdir(output_folder, parents=True)\n",
"for R0, infmean, cbr in zip(R0_samples, infmean_samples, cbr_samples):\n",
" parameters = PropertySet({\"seed\": np.random.randint(0, 1000000),\n",
" \"nticks\": nticks,\n",
" \"verbose\": True, \n",
" \"beta\": R0/infmean,\n",
" \"inf_mean\": infmean,\n",
" \"cbr\": cbr, \n",
" \"importation_period\": 180, \n",
" \"importation_end\": 20*365\n",
" })\n",
"\n",
" mu = ((1+parameters.cbr/1000)**(1/365)-1) \n",
" parameters = PropertySet(\n",
" {\n",
" \"seed\": np.random.randint(0, 1000000),\n",
" \"nticks\": nticks,\n",
" \"verbose\": True,\n",
" \"beta\": R0 / infmean,\n",
" \"inf_mean\": infmean,\n",
" \"cbr\": cbr,\n",
" \"importation_period\": 180,\n",
" \"importation_end\": 20 * 365,\n",
" }\n",
" )\n",
"\n",
" mu = (1 + parameters.cbr / 1000) ** (1 / 365) - 1\n",
"\n",
" model = Model(scenario, parameters)\n",
" model.components = [Births_ConstantPop,\n",
" model.components = [\n",
" Births_ConstantPop,\n",
" Susceptibility,\n",
" Transmission,\n",
" Infection,\n",
" Infect_Agents_In_Patch,\n",
" ] \n",
" ]\n",
"\n",
" #Start them slightly asynchronously - different initial susceptibilities, infection only in 1 patch\n",
" #Want to see how connectivity drives correlation over time.\n",
" # Start them slightly asynchronously - different initial susceptibilities, infection only in 1 patch\n",
" # Want to see how connectivity drives correlation over time.\n",
" for j in range(npatches):\n",
" set_initial_susceptibility_in_patch(model, j, 1/R0+.1/R0*np.random.normal())\n",
" set_initial_susceptibility_in_patch(model, j, 1 / R0 + 0.1 / R0 * np.random.normal())\n",
"\n",
" model.run()\n",
" outputs[i, :, :] = model.patches.cases\n",
" np.save(f\"{output_folder}/CCSSIRoutputs_{i}.npy\", outputs[i, :, :])\n",
" i+=1\n",
"\n"
" np.save(output_folder / f\"CCSSIRoutputs_{i}.npy\", outputs[i, :, :])\n",
" i += 1"
]
},
{
Expand All @@ -139,7 +137,7 @@
" 'cbr': cbr_samples\n",
"})\n",
"\n",
"params_df.to_csv(os.path.join(output_folder, 'params.csv'), index=False)"
"params_df.to_csv(output_folder / 'params.csv', index = False)"
]
},
{
Expand All @@ -159,12 +157,12 @@
}
],
"source": [
"plt.imshow(outputs[26, 7300:, :].T/pops[:, np.newaxis], aspect='auto', origin='lower')\n",
"plt.colorbar(label='Cases')\n",
"plt.xlabel('Time (days)')\n",
"plt.ylabel('Patch')\n",
"plt.imshow(outputs[26, 7300:, :].T / pops[:, np.newaxis], aspect=\"auto\", origin=\"lower\")\n",
"plt.colorbar(label=\"Cases\")\n",
"plt.xlabel(\"Time (days)\")\n",
"plt.ylabel(\"Patch\")\n",
"plt.yticks(range(0, npatches, 10), np.log10(pops[::10]))\n",
"plt.title('Infection spread over time across patches')\n",
"plt.title(\"Infection spread over time across patches\")\n",
"plt.show()"
]
},
Expand Down Expand Up @@ -261,13 +259,9 @@
}
],
"source": [
"params_df = pd.DataFrame({\n",
" 'R0': R0_samples,\n",
" 'infmean': infmean_samples,\n",
" 'cbr': cbr_samples\n",
"})\n",
"params_df = pd.DataFrame({\"R0\": R0_samples, \"infmean\": infmean_samples, \"cbr\": cbr_samples})\n",
"\n",
"params_df.to_csv('CCSSIRoutputs/params.csv', index=False)"
"params_df.to_csv(\"CCSSIRoutputs/params.csv\", index=False)"
]
},
{
Expand Down Expand Up @@ -319,21 +313,18 @@
" end_output = outputs[sim, -1, :]\n",
" zero_pops = pops[end_output == 0]\n",
" nonzero_pops = pops[end_output != 0]\n",
" \n",
"\n",
" if len(zero_pops) > 0:\n",
" CCS2.append(np.max(zero_pops))\n",
" else:\n",
" CCS2.append(None)\n",
" \n",
"\n",
" if len(nonzero_pops) > 0:\n",
" CCS1.append(np.min(nonzero_pops))\n",
" else:\n",
" CCS1.append(None)\n",
"\n",
"results_df = pd.DataFrame({\n",
" 'largest_zero_pop': CCS2,\n",
" 'smallest_nonzero_pop': CCS1\n",
"})\n",
"results_df = pd.DataFrame({\"largest_zero_pop\": CCS2, \"smallest_nonzero_pop\": CCS1})\n",
"\n",
"print(results_df)"
]
Expand All @@ -358,42 +349,42 @@
"fig, axs = plt.subplots(3, 2, figsize=(15, 15))\n",
"\n",
"# Plot largest_zero_pop against R0, infmean, and cbr\n",
"axs[0, 0].scatter(params_df['R0'], results_df['largest_zero_pop'])\n",
"axs[0, 0].set_xlabel('R0')\n",
"axs[0, 0].set_ylabel('Largest Zero Pop')\n",
"axs[0, 0].set_title('Largest Zero Pop vs R0')\n",
"axs[0, 0].set_yscale('log')\n",
"\n",
"axs[1, 0].scatter(params_df['infmean'], results_df['largest_zero_pop'])\n",
"axs[1, 0].set_xlabel('Infectious Mean Period')\n",
"axs[1, 0].set_ylabel('Largest Zero Pop')\n",
"axs[1, 0].set_title('Largest Zero Pop vs Infectious Mean Period')\n",
"axs[1, 0].set_yscale('log')\n",
"\n",
"axs[2, 0].scatter(params_df['cbr'], results_df['largest_zero_pop'])\n",
"axs[2, 0].set_xlabel('Contact Birth Rate')\n",
"axs[2, 0].set_ylabel('Largest Zero Pop')\n",
"axs[2, 0].set_title('Largest Zero Pop vs Contact Birth Rate')\n",
"axs[2, 0].set_yscale('log')\n",
"axs[0, 0].scatter(params_df[\"R0\"], results_df[\"largest_zero_pop\"])\n",
"axs[0, 0].set_xlabel(\"R0\")\n",
"axs[0, 0].set_ylabel(\"Largest Zero Pop\")\n",
"axs[0, 0].set_title(\"Largest Zero Pop vs R0\")\n",
"axs[0, 0].set_yscale(\"log\")\n",
"\n",
"axs[1, 0].scatter(params_df[\"infmean\"], results_df[\"largest_zero_pop\"])\n",
"axs[1, 0].set_xlabel(\"Infectious Mean Period\")\n",
"axs[1, 0].set_ylabel(\"Largest Zero Pop\")\n",
"axs[1, 0].set_title(\"Largest Zero Pop vs Infectious Mean Period\")\n",
"axs[1, 0].set_yscale(\"log\")\n",
"\n",
"axs[2, 0].scatter(params_df[\"cbr\"], results_df[\"largest_zero_pop\"])\n",
"axs[2, 0].set_xlabel(\"Contact Birth Rate\")\n",
"axs[2, 0].set_ylabel(\"Largest Zero Pop\")\n",
"axs[2, 0].set_title(\"Largest Zero Pop vs Contact Birth Rate\")\n",
"axs[2, 0].set_yscale(\"log\")\n",
"\n",
"# Plot smallest_nonzero_pop against R0, infmean, and cbr\n",
"axs[0, 1].scatter(params_df['R0'], results_df['smallest_nonzero_pop'])\n",
"axs[0, 1].set_xlabel('R0')\n",
"axs[0, 1].set_ylabel('Smallest Nonzero Pop')\n",
"axs[0, 1].set_title('Smallest Nonzero Pop vs R0')\n",
"axs[0, 1].set_yscale('log')\n",
"\n",
"axs[1, 1].scatter(params_df['infmean'], results_df['smallest_nonzero_pop'])\n",
"axs[1, 1].set_xlabel('Infectious Mean Period')\n",
"axs[1, 1].set_ylabel('Smallest Nonzero Pop')\n",
"axs[1, 1].set_title('Smallest Nonzero Pop vs Infectious Mean Period')\n",
"axs[1, 1].set_yscale('log')\n",
"\n",
"axs[2, 1].scatter(params_df['cbr'], results_df['smallest_nonzero_pop'])\n",
"axs[2, 1].set_xlabel('Contact Birth Rate')\n",
"axs[2, 1].set_ylabel('Smallest Nonzero Pop')\n",
"axs[2, 1].set_title('Smallest Nonzero Pop vs Contact Birth Rate')\n",
"axs[2, 1].set_yscale('log')\n",
"axs[0, 1].scatter(params_df[\"R0\"], results_df[\"smallest_nonzero_pop\"])\n",
"axs[0, 1].set_xlabel(\"R0\")\n",
"axs[0, 1].set_ylabel(\"Smallest Nonzero Pop\")\n",
"axs[0, 1].set_title(\"Smallest Nonzero Pop vs R0\")\n",
"axs[0, 1].set_yscale(\"log\")\n",
"\n",
"axs[1, 1].scatter(params_df[\"infmean\"], results_df[\"smallest_nonzero_pop\"])\n",
"axs[1, 1].set_xlabel(\"Infectious Mean Period\")\n",
"axs[1, 1].set_ylabel(\"Smallest Nonzero Pop\")\n",
"axs[1, 1].set_title(\"Smallest Nonzero Pop vs Infectious Mean Period\")\n",
"axs[1, 1].set_yscale(\"log\")\n",
"\n",
"axs[2, 1].scatter(params_df[\"cbr\"], results_df[\"smallest_nonzero_pop\"])\n",
"axs[2, 1].set_xlabel(\"Contact Birth Rate\")\n",
"axs[2, 1].set_ylabel(\"Smallest Nonzero Pop\")\n",
"axs[2, 1].set_title(\"Smallest Nonzero Pop vs Contact Birth Rate\")\n",
"axs[2, 1].set_yscale(\"log\")\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
Expand Down Expand Up @@ -441,28 +432,26 @@
}
],
"source": [
"from mpl_toolkits.mplot3d import Axes3D\n",
"\n",
"fig = plt.figure(figsize=(14, 7))\n",
"\n",
"# Surface plot for largest_zero_pop\n",
"ax1 = fig.add_subplot(121, projection='3d')\n",
"ax1.plot_trisurf(params_df['infmean'], params_df['cbr'], results_df['largest_zero_pop'], cmap='viridis')\n",
"ax1.set_xlabel('Infectious Mean Period')\n",
"ax1.set_ylabel('Contact Birth Rate')\n",
"ax1.set_zlabel('Largest Zero Pop')\n",
"ax1.set_zscale('log')\n",
"ax1 = fig.add_subplot(121, projection=\"3d\")\n",
"ax1.plot_trisurf(params_df[\"infmean\"], params_df[\"cbr\"], results_df[\"largest_zero_pop\"], cmap=\"viridis\")\n",
"ax1.set_xlabel(\"Infectious Mean Period\")\n",
"ax1.set_ylabel(\"Contact Birth Rate\")\n",
"ax1.set_zlabel(\"Largest Zero Pop\")\n",
"ax1.set_zscale(\"log\")\n",
"\n",
"ax1.set_title('Largest Zero Pop vs Infectious Mean Period and Contact Birth Rate')\n",
"ax1.set_title(\"Largest Zero Pop vs Infectious Mean Period and Contact Birth Rate\")\n",
"\n",
"# Surface plot for smallest_nonzero_pop\n",
"ax2 = fig.add_subplot(122, projection='3d')\n",
"ax2.plot_trisurf(params_df['infmean'], params_df['cbr'], results_df['smallest_nonzero_pop'], cmap='viridis')\n",
"ax2.set_xlabel('Infectious Mean Period')\n",
"ax2.set_ylabel('Contact Birth Rate')\n",
"ax2.set_zlabel('Smallest Nonzero Pop')\n",
"ax2.set_title('Smallest Nonzero Pop vs Infectious Mean Period and Contact Birth Rate')\n",
"ax2.set_zscale('log')\n",
"ax2 = fig.add_subplot(122, projection=\"3d\")\n",
"ax2.plot_trisurf(params_df[\"infmean\"], params_df[\"cbr\"], results_df[\"smallest_nonzero_pop\"], cmap=\"viridis\")\n",
"ax2.set_xlabel(\"Infectious Mean Period\")\n",
"ax2.set_ylabel(\"Contact Birth Rate\")\n",
"ax2.set_zlabel(\"Smallest Nonzero Pop\")\n",
"ax2.set_title(\"Smallest Nonzero Pop vs Infectious Mean Period and Contact Birth Rate\")\n",
"ax2.set_zscale(\"log\")\n",
"plt.tight_layout()\n",
"plt.show()"
]
Expand Down Expand Up @@ -502,19 +491,19 @@
"plt.figure(figsize=(12, 6))\n",
"\n",
"plt.subplot(1, 2, 1)\n",
"plt.scatter(1/alpha, results_df['largest_zero_pop'])\n",
"plt.xlabel('Alpha (inf_mean * cbr)')\n",
"plt.ylabel('Largest Zero Pop')\n",
"plt.title('Largest Zero Pop vs Alpha')\n",
"plt.yscale('log')\n",
"plt.scatter(1 / alpha, results_df[\"largest_zero_pop\"])\n",
"plt.xlabel(\"Alpha (inf_mean * cbr)\")\n",
"plt.ylabel(\"Largest Zero Pop\")\n",
"plt.title(\"Largest Zero Pop vs Alpha\")\n",
"plt.yscale(\"log\")\n",
"\n",
"# Plot smallest_nonzero_pop against alpha\n",
"plt.subplot(1, 2, 2)\n",
"plt.scatter(1/alpha, results_df['smallest_nonzero_pop'])\n",
"plt.xlabel('Alpha (inf_mean * cbr)')\n",
"plt.ylabel('Smallest Nonzero Pop')\n",
"plt.title('Smallest Nonzero Pop vs Alpha')\n",
"#plt.yscale('line')\n",
"plt.scatter(1 / alpha, results_df[\"smallest_nonzero_pop\"])\n",
"plt.xlabel(\"Alpha (inf_mean * cbr)\")\n",
"plt.ylabel(\"Smallest Nonzero Pop\")\n",
"plt.title(\"Smallest Nonzero Pop vs Alpha\")\n",
"# plt.yscale('line')\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
Expand Down Expand Up @@ -564,7 +553,7 @@
}
],
"source": [
"plt.hist(10+cbr_samples, bins=50)"
"plt.hist(10 + cbr_samples, bins=50)"
]
},
{
Expand All @@ -573,7 +562,7 @@
"metadata": {},
"outputs": [],
"source": [
"output_folder = \"..\\..\\laser-generic-outputs\\CCSSIRoutputs\""
"output_folder = Path(r\"..\\..\\laser-generic-outputs\\CCSSIRoutputs\")"
]
},
{
Expand Down
Loading

0 comments on commit 6faff0a

Please sign in to comment.