From b518c343b4f88890e7e84357e228b812a1aae5ba Mon Sep 17 00:00:00 2001 From: Izaac Molina <39788504+IzMo2000@users.noreply.github.com> Date: Tue, 30 Apr 2024 12:04:57 -0700 Subject: [PATCH] Error handling - IPM Simulation (#104) Implemented error handling for common simulation errors, namely: 1. Incorrect parameters causing a negative or zero rate 2. Divide by zero errors, mainly occurring when the population of a node is 0 3. Incorrect parameters causing a negative probability in fork definitions 4. Incompatible clause/parameter values for MM simulations For more details, check devlog 2024-04-16.ipynb --- doc/devlog/2024-04-04-draw-demo.ipynb | 2 + doc/devlog/2024-04-16.ipynb | 796 ++++++++++++++++++++++++++ doc/devlog/README.md | 1 + epymorph/__init__.py | 6 + epymorph/engine/ipm_exec.py | 85 ++- epymorph/error.py | 69 +++ epymorph/movement/movement_model.py | 8 +- epymorph/test/simulate_test.py | 179 ++++++ 8 files changed, 1143 insertions(+), 3 deletions(-) create mode 100644 doc/devlog/2024-04-16.ipynb diff --git a/doc/devlog/2024-04-04-draw-demo.ipynb b/doc/devlog/2024-04-04-draw-demo.ipynb index 5d04d666..b827ba8b 100644 --- a/doc/devlog/2024-04-04-draw-demo.ipynb +++ b/doc/devlog/2024-04-04-draw-demo.ipynb @@ -7,6 +7,8 @@ "# 2024/04/04 Draw Module Demo\n", "_Author: Izaac Molina_\n", "\n", + "Provides demos on using the `draw` module in epymorph\n", + "\n", "## Goal and Purpose\n", "The traditional way to design a model in epymorph is to first draw a ipm model, then translate that model into an ipm object via code. This project seeks to go in the other direction: turning an ipm object into a model. This is great both for educational and debugging purposes as it allows users to present their model more clearly as well as check to ensure they have implemented the model correctly in epymorph\n", "\n", diff --git a/doc/devlog/2024-04-16.ipynb b/doc/devlog/2024-04-16.ipynb new file mode 100644 index 00000000..d6f570de --- /dev/null +++ b/doc/devlog/2024-04-16.ipynb @@ -0,0 +1,796 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 2024/04/16 Error Handling\n", + "_Author: Izaac Molina and Tyler Coles_\n", + "\n", + "Details additional exceptions implemented to catch common simulation errors" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Error Handling\n", + "\n", + "There are (at least) four error cases we're interested in trapping --\n", + "1. something causes an IPM expression to go to zero or a negative value and results in a poisson draw with a lambda <= 0, and\n", + "2. something causes the population of a node to go to zero (when a divisor isn't \"protected\" by a construction like Max(1, ...)) causing a divide by zero.\n", + "3. similar to issue 1, a negative parameter can result in a negative probability in fork definitions\n", + "4. an incompatible clause and parameter value are used in the movement model execution, leading to exceptions\n", + "\n", + "## First, we define ipm functions for testing purposes" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "from epymorph import *\n", + "from epymorph.compartment_model import (CompartmentModel, compartment,\n", + " create_model, create_symbols, edge,\n", + " fork, param)\n", + "from epymorph.data_shape import Shapes\n", + "from epymorph.params import ParamValue\n", + "\n", + "\n", + "def load_ipm() -> CompartmentModel:\n", + " \"\"\"Load the 'sirs' IPM.\"\"\"\n", + " symbols = create_symbols(\n", + " compartments=[\n", + " compartment('S'),\n", + " compartment('I'),\n", + " compartment('R'),\n", + " ],\n", + " attributes=[\n", + " param('beta', shape=Shapes.TxN), # infectivity\n", + " param('gamma', shape=Shapes.TxN), # progression from infected to recovered\n", + " param('xi', shape=Shapes.TxN) # progression from recovered to susceptible\n", + " ])\n", + "\n", + " [S, I, R] = symbols.compartment_symbols\n", + " [β, γ, ξ] = symbols.attribute_symbols\n", + "\n", + " # LOOK!\n", + " # N is NOT protected by Max(1, ...) here\n", + " # This isn't necessary for Case 1, but is necessary for Case 2.\n", + " N = S + I + R\n", + "\n", + " return create_model(\n", + " symbols=symbols,\n", + " transitions=[\n", + " edge(S, I, rate=β * S * I / N),\n", + " edge(I, R, rate=γ * I),\n", + " edge(R, S, rate=ξ * R)\n", + " ])" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def load_ipm_fork() -> CompartmentModel:\n", + " \"\"\"Load the 'sirs' IPM.\"\"\"\n", + " symbols = create_symbols(\n", + " compartments=[\n", + " compartment('S'),\n", + " compartment('I'),\n", + " compartment('R'),\n", + " ],\n", + " attributes=[\n", + " param('beta', shape=Shapes.TxN), # infectivity\n", + " param('gamma', shape=Shapes.TxN), # progression from infected to recovered\n", + " param('xi', shape=Shapes.TxN), # progression from recovered to susceptible\n", + " param('prob1', shape=Shapes.TxN),\n", + " param('prob2', shape=Shapes.TxN)\n", + " ])\n", + "\n", + " [S, I, R] = symbols.compartment_symbols\n", + " [β, γ, ξ, prob1, prob2] = symbols.attribute_symbols\n", + "\n", + " # LOOK!\n", + " # N is NOT protected by Max(1, ...) here\n", + " # This isn't necessary for Case 1, but is necessary for Case 2.\n", + " N = S + I + R\n", + "\n", + " return create_model(\n", + " symbols=symbols,\n", + " transitions=[\n", + " fork(\n", + " edge(S, I, rate=(β * S * I / N) * prob1),\n", + " edge(S, R, rate=(β * S * I / N) * prob2),\n", + " ),\n", + " edge(I, R, rate=γ * I),\n", + " edge(R, S, rate=ξ * R)\n", + " ])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Case 1:\n", + "\n", + "If we just provide a negative value to beta in a classic SIRS model, this causes the first error.\n", + "\n", + "To handle this, the exception type `IpmSimLesstThanZeroException` is used in `ipm_exec.py` to catch the less than zero error." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running simulation (StandardSimulation):\n", + "• 2015-01-01 to 2015-05-31 (150 days)\n", + "• 6 geo nodes\n", + "| | 0% \r" + ] + }, + { + "ename": "IpmSimLessThanZeroException", + "evalue": "\nLess than zero rate detected. When providing or defining ipm parameters, ensure that\nthey will not result in a negative rate. Note: this can often happen unintentionally\nif a function is given as a parameter.\n\nShowing current Node : Timestep\n1: 0\n\nShowing current compartment values\nS: 9687648\nI: 5\nR: 0\n\nShowing current ipm params\nbeta: 0.4\ngamma: -0.2\nxi: 0.01\n\n", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mIpmSimLessThanZeroException\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[3], line 23\u001b[0m\n\u001b[0;32m 7\u001b[0m sim \u001b[38;5;241m=\u001b[39m StandardSimulation(\n\u001b[0;32m 8\u001b[0m geo\u001b[38;5;241m=\u001b[39mpei_geo,\n\u001b[0;32m 9\u001b[0m ipm\u001b[38;5;241m=\u001b[39mload_ipm(),\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 19\u001b[0m rng\u001b[38;5;241m=\u001b[39mdefault_rng(\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m 20\u001b[0m )\n\u001b[0;32m 22\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m sim_messaging(sim):\n\u001b[1;32m---> 23\u001b[0m out \u001b[38;5;241m=\u001b[39m \u001b[43msim\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\standard_sim.py:180\u001b[0m, in \u001b[0;36mStandardSimulation.run\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 178\u001b[0m \u001b[38;5;66;03m# Then do IPM\u001b[39;00m\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the IPM\u001b[39m\u001b[38;5;124m\"\u001b[39m, IpmSimException, AttributeException):\n\u001b[1;32m--> 180\u001b[0m tick_events, tick_prevalence \u001b[38;5;241m=\u001b[39m \u001b[43mipm_exec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mapply\u001b[49m\u001b[43m(\u001b[49m\u001b[43mworld\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 181\u001b[0m out\u001b[38;5;241m.\u001b[39mincidence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_events\n\u001b[0;32m 182\u001b[0m out\u001b[38;5;241m.\u001b[39mprevalence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_prevalence\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:142\u001b[0m, in \u001b[0;36mStandardIpmExecutor.apply\u001b[1;34m(self, world, tick)\u001b[0m\n\u001b[0;32m 139\u001b[0m cohorts \u001b[38;5;241m=\u001b[39m world\u001b[38;5;241m.\u001b[39mget_cohort_array(node)\n\u001b[0;32m 140\u001b[0m effective \u001b[38;5;241m=\u001b[39m cohorts\u001b[38;5;241m.\u001b[39msum(axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, dtype\u001b[38;5;241m=\u001b[39mSimDType)\n\u001b[1;32m--> 142\u001b[0m occurrences \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_events\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43meffective\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 143\u001b[0m cohort_deltas \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_distribute(cohorts, occurrences)\n\u001b[0;32m 144\u001b[0m world\u001b[38;5;241m.\u001b[39mapply_cohort_delta(node, cohort_deltas)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:175\u001b[0m, in \u001b[0;36mStandardIpmExecutor._events\u001b[1;34m(self, node, tick, effective_pop)\u001b[0m\n\u001b[0;32m 173\u001b[0m \u001b[38;5;66;03m# check for < 0 rate, throw error in this case\u001b[39;00m\n\u001b[0;32m 174\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m rate \u001b[38;5;241m<\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[1;32m--> 175\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m IpmSimLessThanZeroException(\n\u001b[0;32m 176\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_default_error_args(rate_args, node, tick)\n\u001b[0;32m 177\u001b[0m )\n\u001b[0;32m 178\u001b[0m occur[index] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx\u001b[38;5;241m.\u001b[39mrng\u001b[38;5;241m.\u001b[39mpoisson(rate \u001b[38;5;241m*\u001b[39m tick\u001b[38;5;241m.\u001b[39mtau)\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mcase\u001b[39;00m _ForkedTrx(size, rate_lambda, prob_lambda):\n\u001b[0;32m 180\u001b[0m \u001b[38;5;66;03m# get rate from lambda expression, catch divide by zero error\u001b[39;00m\n", + "\u001b[1;31mIpmSimLessThanZeroException\u001b[0m: \nLess than zero rate detected. When providing or defining ipm parameters, ensure that\nthey will not result in a negative rate. Note: this can often happen unintentionally\nif a function is given as a parameter.\n\nShowing current Node : Timestep\n1: 0\n\nShowing current compartment values\nS: 9687648\nI: 5\nR: 0\n\nShowing current ipm params\nbeta: 0.4\ngamma: -0.2\nxi: 0.01\n\n" + ] + } + ], + "source": [ + "from functools import partial\n", + "\n", + "from epymorph.initializer import single_location\n", + "\n", + "pei_geo = geo_library['pei']()\n", + "\n", + "sim = StandardSimulation(\n", + " geo=pei_geo,\n", + " ipm=load_ipm(),\n", + " mm=mm_library['no'](),\n", + " params={\n", + " 'phi': 40.0,\n", + " 'beta': 0.4,\n", + " 'gamma': -1 / 5, # NEGATIVE NUMBER! WATCH OUT!\n", + " 'xi': 1 / 100,\n", + " },\n", + " time_frame=TimeFrame.of(\"2015-01-01\", 150),\n", + " initializer=partial(single_location, location=1, seed_size=5),\n", + " rng=default_rng(1)\n", + ")\n", + "\n", + "with sim_messaging(sim):\n", + " out = sim.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Test using function as a parameter" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running simulation (StandardSimulation):\n", + "• 2015-01-01 to 2015-05-31 (150 days)\n", + "• 6 geo nodes\n", + "| | 0% \r" + ] + }, + { + "ename": "IpmSimLessThanZeroException", + "evalue": "\nLess than zero rate detected. When providing or defining ipm parameters, ensure that\nthey will not result in a negative rate. Note: this can often happen unintentionally\nif a function is given as a parameter.\n\nShowing current Node : Timestep\n0: 0\n\nShowing current compartment values\nS: 18799002\nI: 9990\nR: 0\n\nShowing current ipm params\nbeta: -0.35\ngamma: 0.16666666666666666\nxi: 0.011111111111111112\n\n", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mIpmSimLessThanZeroException\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[4], line 35\u001b[0m\n\u001b[0;32m 31\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m x\n\u001b[0;32m 34\u001b[0m \u001b[38;5;66;03m# The function is our beta value!\u001b[39;00m\n\u001b[1;32m---> 35\u001b[0m \u001b[43mrun_and_plot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mbeta_fn\u001b[49m\u001b[43m)\u001b[49m\n", + "Cell \u001b[1;32mIn[4], line 20\u001b[0m, in \u001b[0;36mrun_and_plot\u001b[1;34m(beta)\u001b[0m\n\u001b[0;32m 2\u001b[0m sim \u001b[38;5;241m=\u001b[39m StandardSimulation(\n\u001b[0;32m 3\u001b[0m ipm\u001b[38;5;241m=\u001b[39mipm_library[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msirs\u001b[39m\u001b[38;5;124m'\u001b[39m](),\n\u001b[0;32m 4\u001b[0m mm\u001b[38;5;241m=\u001b[39mmm_library[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpei\u001b[39m\u001b[38;5;124m'\u001b[39m](),\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 16\u001b[0m initializer\u001b[38;5;241m=\u001b[39mpartial(single_location, location\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, seed_size\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m10_000\u001b[39m),\n\u001b[0;32m 17\u001b[0m )\n\u001b[0;32m 19\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m sim_messaging(sim):\n\u001b[1;32m---> 20\u001b[0m \u001b[43msim\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\standard_sim.py:180\u001b[0m, in \u001b[0;36mStandardSimulation.run\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 178\u001b[0m \u001b[38;5;66;03m# Then do IPM\u001b[39;00m\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the IPM\u001b[39m\u001b[38;5;124m\"\u001b[39m, IpmSimException, AttributeException):\n\u001b[1;32m--> 180\u001b[0m tick_events, tick_prevalence \u001b[38;5;241m=\u001b[39m \u001b[43mipm_exec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mapply\u001b[49m\u001b[43m(\u001b[49m\u001b[43mworld\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 181\u001b[0m out\u001b[38;5;241m.\u001b[39mincidence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_events\n\u001b[0;32m 182\u001b[0m out\u001b[38;5;241m.\u001b[39mprevalence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_prevalence\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:142\u001b[0m, in \u001b[0;36mStandardIpmExecutor.apply\u001b[1;34m(self, world, tick)\u001b[0m\n\u001b[0;32m 139\u001b[0m cohorts \u001b[38;5;241m=\u001b[39m world\u001b[38;5;241m.\u001b[39mget_cohort_array(node)\n\u001b[0;32m 140\u001b[0m effective \u001b[38;5;241m=\u001b[39m cohorts\u001b[38;5;241m.\u001b[39msum(axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, dtype\u001b[38;5;241m=\u001b[39mSimDType)\n\u001b[1;32m--> 142\u001b[0m occurrences \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_events\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43meffective\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 143\u001b[0m cohort_deltas \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_distribute(cohorts, occurrences)\n\u001b[0;32m 144\u001b[0m world\u001b[38;5;241m.\u001b[39mapply_cohort_delta(node, cohort_deltas)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:175\u001b[0m, in \u001b[0;36mStandardIpmExecutor._events\u001b[1;34m(self, node, tick, effective_pop)\u001b[0m\n\u001b[0;32m 173\u001b[0m \u001b[38;5;66;03m# check for < 0 rate, throw error in this case\u001b[39;00m\n\u001b[0;32m 174\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m rate \u001b[38;5;241m<\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[1;32m--> 175\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m IpmSimLessThanZeroException(\n\u001b[0;32m 176\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_default_error_args(rate_args, node, tick)\n\u001b[0;32m 177\u001b[0m )\n\u001b[0;32m 178\u001b[0m occur[index] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx\u001b[38;5;241m.\u001b[39mrng\u001b[38;5;241m.\u001b[39mpoisson(rate \u001b[38;5;241m*\u001b[39m tick\u001b[38;5;241m.\u001b[39mtau)\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mcase\u001b[39;00m _ForkedTrx(size, rate_lambda, prob_lambda):\n\u001b[0;32m 180\u001b[0m \u001b[38;5;66;03m# get rate from lambda expression, catch divide by zero error\u001b[39;00m\n", + "\u001b[1;31mIpmSimLessThanZeroException\u001b[0m: \nLess than zero rate detected. When providing or defining ipm parameters, ensure that\nthey will not result in a negative rate. Note: this can often happen unintentionally\nif a function is given as a parameter.\n\nShowing current Node : Timestep\n0: 0\n\nShowing current compartment values\nS: 18799002\nI: 9990\nR: 0\n\nShowing current ipm params\nbeta: -0.35\ngamma: 0.16666666666666666\nxi: 0.011111111111111112\n\n" + ] + } + ], + "source": [ + "def run_and_plot(beta: ParamValue) -> None:\n", + " sim = StandardSimulation(\n", + " ipm=ipm_library['sirs'](),\n", + " mm=mm_library['pei'](),\n", + " geo=geo_library['pei'](),\n", + " params={\n", + " # IPM params\n", + " 'beta': beta,\n", + " 'gamma': 1 / 6,\n", + " 'xi': 1 / 90,\n", + " # MM params\n", + " 'theta': 0.1,\n", + " 'move_control': 0.9,\n", + " },\n", + " time_frame=TimeFrame.of(\"2015-01-01\", 150),\n", + " initializer=partial(single_location, location=0, seed_size=10_000),\n", + " )\n", + "\n", + " with sim_messaging(sim):\n", + " sim.run()\n", + "\n", + "\n", + "def beta_fn(t, n):\n", + " # NEGATIVE VALUES PRODUCED BY FUNCTION\n", + " x = -0.35 + -0.05 * np.sin(-2 * np.pi * (t / dim.days))\n", + " cutoff = 50 + (n * 3)\n", + " if t > cutoff:\n", + " pop = geo['population'][n]\n", + " cut = 0.3 if pop < 9_000_000 else 0.25\n", + " x -= cut\n", + " return x\n", + "\n", + "\n", + "# The function is our beta value!\n", + "run_and_plot(beta_fn)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Test with fork" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running simulation (StandardSimulation):\n", + "• 2015-01-01 to 2015-05-31 (150 days)\n", + "• 6 geo nodes\n", + "| | 0% \r" + ] + }, + { + "ename": "IpmSimLessThanZeroException", + "evalue": "\nLess than zero rate detected. When providing or defining ipm parameters, ensure that\nthey will not result in a negative rate. Note: this can often happen unintentionally\nif a function is given as a parameter.\n\nShowing current Node : Timestep\n1: 0\n\nShowing current compartment values\nS: 9687648\nI: 5\nR: 0\n\nShowing current ipm params\nbeta: -0.4\ngamma: 0.2\nxi: 0.01\nprob1: 0.5\nprob2: 0.5\n\n", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mIpmSimLessThanZeroException\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[5], line 26\u001b[0m\n\u001b[0;32m 8\u001b[0m sim \u001b[38;5;241m=\u001b[39m StandardSimulation(\n\u001b[0;32m 9\u001b[0m geo\u001b[38;5;241m=\u001b[39mpei_geo,\n\u001b[0;32m 10\u001b[0m ipm\u001b[38;5;241m=\u001b[39mload_ipm_fork(),\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 22\u001b[0m rng\u001b[38;5;241m=\u001b[39mdefault_rng(\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m 23\u001b[0m )\n\u001b[0;32m 25\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m sim_messaging(sim):\n\u001b[1;32m---> 26\u001b[0m out \u001b[38;5;241m=\u001b[39m \u001b[43msim\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\standard_sim.py:180\u001b[0m, in \u001b[0;36mStandardSimulation.run\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 178\u001b[0m \u001b[38;5;66;03m# Then do IPM\u001b[39;00m\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the IPM\u001b[39m\u001b[38;5;124m\"\u001b[39m, IpmSimException, AttributeException):\n\u001b[1;32m--> 180\u001b[0m tick_events, tick_prevalence \u001b[38;5;241m=\u001b[39m \u001b[43mipm_exec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mapply\u001b[49m\u001b[43m(\u001b[49m\u001b[43mworld\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 181\u001b[0m out\u001b[38;5;241m.\u001b[39mincidence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_events\n\u001b[0;32m 182\u001b[0m out\u001b[38;5;241m.\u001b[39mprevalence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_prevalence\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:142\u001b[0m, in \u001b[0;36mStandardIpmExecutor.apply\u001b[1;34m(self, world, tick)\u001b[0m\n\u001b[0;32m 139\u001b[0m cohorts \u001b[38;5;241m=\u001b[39m world\u001b[38;5;241m.\u001b[39mget_cohort_array(node)\n\u001b[0;32m 140\u001b[0m effective \u001b[38;5;241m=\u001b[39m cohorts\u001b[38;5;241m.\u001b[39msum(axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, dtype\u001b[38;5;241m=\u001b[39mSimDType)\n\u001b[1;32m--> 142\u001b[0m occurrences \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_events\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43meffective\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 143\u001b[0m cohort_deltas \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_distribute(cohorts, occurrences)\n\u001b[0;32m 144\u001b[0m world\u001b[38;5;241m.\u001b[39mapply_cohort_delta(node, cohort_deltas)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:190\u001b[0m, in \u001b[0;36mStandardIpmExecutor._events\u001b[1;34m(self, node, tick, effective_pop)\u001b[0m\n\u001b[0;32m 188\u001b[0m \u001b[38;5;66;03m# check for < 0 base, throw error in this case\u001b[39;00m\n\u001b[0;32m 189\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m rate \u001b[38;5;241m<\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[1;32m--> 190\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m IpmSimLessThanZeroException(\n\u001b[0;32m 191\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_default_error_args(rate_args, node, tick)\n\u001b[0;32m 192\u001b[0m )\n\u001b[0;32m 193\u001b[0m base \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx\u001b[38;5;241m.\u001b[39mrng\u001b[38;5;241m.\u001b[39mpoisson(rate \u001b[38;5;241m*\u001b[39m tick\u001b[38;5;241m.\u001b[39mtau)\n\u001b[0;32m 194\u001b[0m prob \u001b[38;5;241m=\u001b[39m prob_lambda(rate_args)\n", + "\u001b[1;31mIpmSimLessThanZeroException\u001b[0m: \nLess than zero rate detected. When providing or defining ipm parameters, ensure that\nthey will not result in a negative rate. Note: this can often happen unintentionally\nif a function is given as a parameter.\n\nShowing current Node : Timestep\n1: 0\n\nShowing current compartment values\nS: 9687648\nI: 5\nR: 0\n\nShowing current ipm params\nbeta: -0.4\ngamma: 0.2\nxi: 0.01\nprob1: 0.5\nprob2: 0.5\n\n" + ] + } + ], + "source": [ + "from functools import partial\n", + "\n", + "from epymorph import *\n", + "from epymorph.initializer import single_location\n", + "\n", + "pei_geo = geo_library['pei']()\n", + "\n", + "sim = StandardSimulation(\n", + " geo=pei_geo,\n", + " ipm=load_ipm_fork(),\n", + " mm=mm_library['no'](),\n", + " params={\n", + " 'phi': 40.0,\n", + " 'beta': -0.4, # NEGATIVE NUMBER! WATCH OUT!\n", + " 'gamma': 1 / 5,\n", + " 'xi': 1 / 100,\n", + " 'prob1': 1 / 2,\n", + " 'prob2': 1 / 2\n", + " },\n", + " time_frame=TimeFrame.of(\"2015-01-01\", 150),\n", + " initializer=partial(single_location, location=1, seed_size=5),\n", + " rng=default_rng(1)\n", + ")\n", + "\n", + "with sim_messaging(sim):\n", + " out = sim.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Case 2:\n", + "\n", + "If we construct our geo to force a population at a node to be zero, we get the second error. (Note: this can also happen if a movement model moves everyone out of a node, so it's not as easy to fix as checking the geo before-hand.)\n", + "\n", + "To handle this, the exception type `IpmSimNaNException` is used in `ipm_exec.py` to catch the divide by zero error, as it results in a NaN (not a number) rate." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running simulation (StandardSimulation):\n", + "• 2015-01-01 to 2015-05-31 (150 days)\n", + "• 3 geo nodes\n", + "| | 0% \r" + ] + }, + { + "ename": "IpmSimNaNException", + "evalue": "\nNaN (not a number) rate detected. This is often the result of a divide by zero error.\nWhen constructing the IPM, ensure that no edge transitions can result in division by zero\nThis commonly occurs when defining an S->I edge that is (some rate / sum of the compartments)\nTo fix this, change the edge to define the S->I edge as (some rate / Max(1/sum of the the compartments))\nSee examples of this in the provided example ipm definitions in the data/ipms folder.\n\nShowing current Node : Timestep\n0: 0\n\nShowing current compartment values\nS: 0\nI: 0\nR: 0\n\nShowing current ipm params\nbeta: 0.4\ngamma: 0.2\nxi: 0.011111111111111112\n\nShowing current corresponding transition\nS->I: I*S*beta/(I + R + S)\n\n", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mFloatingPointError\u001b[0m Traceback (most recent call last)", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:167\u001b[0m, in \u001b[0;36mStandardIpmExecutor._events\u001b[1;34m(self, node, tick, effective_pop)\u001b[0m\n\u001b[0;32m 166\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m--> 167\u001b[0m rate \u001b[38;5;241m=\u001b[39m \u001b[43mrate_lambda\u001b[49m\u001b[43m(\u001b[49m\u001b[43mrate_args\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 168\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m (\u001b[38;5;167;01mZeroDivisionError\u001b[39;00m, \u001b[38;5;167;01mFloatingPointError\u001b[39;00m):\n", + "File \u001b[1;32m:3\u001b[0m, in \u001b[0;36m_lambdifygenerated\u001b[1;34m(_Dummy_46)\u001b[0m\n\u001b[0;32m 2\u001b[0m [S, I, R, beta, gamma, xi] \u001b[38;5;241m=\u001b[39m _Dummy_46\n\u001b[1;32m----> 3\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mI\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mS\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mbeta\u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mI\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mR\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mS\u001b[49m\u001b[43m)\u001b[49m\n", + "\u001b[1;31mFloatingPointError\u001b[0m: invalid value encountered in scalar divide", + "\nDuring handling of the above exception, another exception occurred:\n", + "\u001b[1;31mIpmSimNaNException\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[6], line 41\u001b[0m\n\u001b[0;32m 25\u001b[0m sim \u001b[38;5;241m=\u001b[39m StandardSimulation(\n\u001b[0;32m 26\u001b[0m geo\u001b[38;5;241m=\u001b[39mmy_geo,\n\u001b[0;32m 27\u001b[0m ipm\u001b[38;5;241m=\u001b[39mload_ipm(),\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 37\u001b[0m rng\u001b[38;5;241m=\u001b[39mdefault_rng(\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m 38\u001b[0m )\n\u001b[0;32m 40\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m sim_messaging(sim):\n\u001b[1;32m---> 41\u001b[0m out \u001b[38;5;241m=\u001b[39m \u001b[43msim\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\standard_sim.py:180\u001b[0m, in \u001b[0;36mStandardSimulation.run\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 178\u001b[0m \u001b[38;5;66;03m# Then do IPM\u001b[39;00m\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the IPM\u001b[39m\u001b[38;5;124m\"\u001b[39m, IpmSimException, AttributeException):\n\u001b[1;32m--> 180\u001b[0m tick_events, tick_prevalence \u001b[38;5;241m=\u001b[39m \u001b[43mipm_exec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mapply\u001b[49m\u001b[43m(\u001b[49m\u001b[43mworld\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 181\u001b[0m out\u001b[38;5;241m.\u001b[39mincidence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_events\n\u001b[0;32m 182\u001b[0m out\u001b[38;5;241m.\u001b[39mprevalence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_prevalence\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:142\u001b[0m, in \u001b[0;36mStandardIpmExecutor.apply\u001b[1;34m(self, world, tick)\u001b[0m\n\u001b[0;32m 139\u001b[0m cohorts \u001b[38;5;241m=\u001b[39m world\u001b[38;5;241m.\u001b[39mget_cohort_array(node)\n\u001b[0;32m 140\u001b[0m effective \u001b[38;5;241m=\u001b[39m cohorts\u001b[38;5;241m.\u001b[39msum(axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, dtype\u001b[38;5;241m=\u001b[39mSimDType)\n\u001b[1;32m--> 142\u001b[0m occurrences \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_events\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43meffective\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 143\u001b[0m cohort_deltas \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_distribute(cohorts, occurrences)\n\u001b[0;32m 144\u001b[0m world\u001b[38;5;241m.\u001b[39mapply_cohort_delta(node, cohort_deltas)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:169\u001b[0m, in \u001b[0;36mStandardIpmExecutor._events\u001b[1;34m(self, node, tick, effective_pop)\u001b[0m\n\u001b[0;32m 167\u001b[0m rate \u001b[38;5;241m=\u001b[39m rate_lambda(rate_args)\n\u001b[0;32m 168\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m (\u001b[38;5;167;01mZeroDivisionError\u001b[39;00m, \u001b[38;5;167;01mFloatingPointError\u001b[39;00m):\n\u001b[1;32m--> 169\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m IpmSimNaNException(\n\u001b[0;32m 170\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_zero_division_args(\n\u001b[0;32m 171\u001b[0m rate_args, node, tick, t)\n\u001b[0;32m 172\u001b[0m )\n\u001b[0;32m 173\u001b[0m \u001b[38;5;66;03m# check for < 0 rate, throw error in this case\u001b[39;00m\n\u001b[0;32m 174\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m rate \u001b[38;5;241m<\u001b[39m \u001b[38;5;241m0\u001b[39m:\n", + "\u001b[1;31mIpmSimNaNException\u001b[0m: \nNaN (not a number) rate detected. This is often the result of a divide by zero error.\nWhen constructing the IPM, ensure that no edge transitions can result in division by zero\nThis commonly occurs when defining an S->I edge that is (some rate / sum of the compartments)\nTo fix this, change the edge to define the S->I edge as (some rate / Max(1/sum of the the compartments))\nSee examples of this in the provided example ipm definitions in the data/ipms folder.\n\nShowing current Node : Timestep\n0: 0\n\nShowing current compartment values\nS: 0\nI: 0\nR: 0\n\nShowing current ipm params\nbeta: 0.4\ngamma: 0.2\nxi: 0.011111111111111112\n\nShowing current corresponding transition\nS->I: I*S*beta/(I + R + S)\n\n" + ] + } + ], + "source": [ + "from functools import partial\n", + "\n", + "import numpy as np\n", + "\n", + "from epymorph import *\n", + "from epymorph.geo.spec import NO_DURATION, AttribDef, StaticGeoSpec\n", + "from epymorph.geo.static import StaticGeo\n", + "from epymorph.initializer import single_location\n", + "\n", + "my_geo = StaticGeo(\n", + " spec=StaticGeoSpec(\n", + " attributes=[\n", + " AttribDef('label', dtype=str, shape=Shapes.N),\n", + " AttribDef('population', dtype=str, shape=Shapes.N),\n", + " ],\n", + " time_period=NO_DURATION,\n", + " ),\n", + " values={\n", + " 'label': np.array(['a', 'b', 'c']),\n", + " 'population': np.array([0, 10, 20], dtype=np.int64),\n", + " },\n", + ")\n", + "\n", + "\n", + "sim = StandardSimulation(\n", + " geo=my_geo,\n", + " ipm=load_ipm(),\n", + " mm=mm_library['no'](),\n", + " params={\n", + " 'phi': 40.0,\n", + " 'beta': 0.4,\n", + " 'gamma': 1 / 5,\n", + " 'xi': 1 / 90,\n", + " },\n", + " time_frame=TimeFrame.of(\"2015-01-01\", 150),\n", + " initializer=partial(single_location, location=1, seed_size=5),\n", + " rng=default_rng(1)\n", + ")\n", + "\n", + "with sim_messaging(sim):\n", + " out = sim.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## test with fork" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running simulation (StandardSimulation):\n", + "• 2015-01-01 to 2015-05-31 (150 days)\n", + "• 3 geo nodes\n", + "| | 0% \r" + ] + }, + { + "ename": "IpmSimNaNException", + "evalue": "\nNaN (not a number) rate detected. This is often the result of a divide by zero error.\nWhen constructing the IPM, ensure that no edge transitions can result in division by zero\nThis commonly occurs when defining an S->I edge that is (some rate / sum of the compartments)\nTo fix this, change the edge to define the S->I edge as (some rate / Max(1/sum of the the compartments))\nSee examples of this in the provided example ipm definitions in the data/ipms folder.\n\nShowing current Node : Timestep\n0: 0\n\nShowing current compartment values\nS: 0\nI: 0\nR: 0\n\nShowing current ipm params\nbeta: 0.4\ngamma: 0.2\nxi: 0.011111111111111112\nprob1: 0.5\nprob2: 0.5\n\nShowing current corresponding fork transition\nS->(I, R): I*S*beta*(prob1 + prob2)/(I + R + S)\n\n", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mFloatingPointError\u001b[0m Traceback (most recent call last)", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:182\u001b[0m, in \u001b[0;36mStandardIpmExecutor._events\u001b[1;34m(self, node, tick, effective_pop)\u001b[0m\n\u001b[0;32m 181\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m--> 182\u001b[0m rate \u001b[38;5;241m=\u001b[39m \u001b[43mrate_lambda\u001b[49m\u001b[43m(\u001b[49m\u001b[43mrate_args\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 183\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m (\u001b[38;5;167;01mZeroDivisionError\u001b[39;00m, \u001b[38;5;167;01mFloatingPointError\u001b[39;00m):\n", + "File \u001b[1;32m:3\u001b[0m, in \u001b[0;36m_lambdifygenerated\u001b[1;34m(_Dummy_49)\u001b[0m\n\u001b[0;32m 2\u001b[0m [S, I, R, beta, gamma, xi, prob1, prob2] \u001b[38;5;241m=\u001b[39m _Dummy_49\n\u001b[1;32m----> 3\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mI\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mS\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mbeta\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mprob1\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mprob2\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mI\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mR\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mS\u001b[49m\u001b[43m)\u001b[49m\n", + "\u001b[1;31mFloatingPointError\u001b[0m: invalid value encountered in scalar divide", + "\nDuring handling of the above exception, another exception occurred:\n", + "\u001b[1;31mIpmSimNaNException\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[7], line 43\u001b[0m\n\u001b[0;32m 25\u001b[0m sim \u001b[38;5;241m=\u001b[39m StandardSimulation(\n\u001b[0;32m 26\u001b[0m geo\u001b[38;5;241m=\u001b[39mmy_geo,\n\u001b[0;32m 27\u001b[0m ipm\u001b[38;5;241m=\u001b[39mload_ipm_fork(),\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 39\u001b[0m rng\u001b[38;5;241m=\u001b[39mdefault_rng(\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m 40\u001b[0m )\n\u001b[0;32m 42\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m sim_messaging(sim):\n\u001b[1;32m---> 43\u001b[0m out \u001b[38;5;241m=\u001b[39m \u001b[43msim\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\standard_sim.py:180\u001b[0m, in \u001b[0;36mStandardSimulation.run\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 178\u001b[0m \u001b[38;5;66;03m# Then do IPM\u001b[39;00m\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the IPM\u001b[39m\u001b[38;5;124m\"\u001b[39m, IpmSimException, AttributeException):\n\u001b[1;32m--> 180\u001b[0m tick_events, tick_prevalence \u001b[38;5;241m=\u001b[39m \u001b[43mipm_exec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mapply\u001b[49m\u001b[43m(\u001b[49m\u001b[43mworld\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 181\u001b[0m out\u001b[38;5;241m.\u001b[39mincidence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_events\n\u001b[0;32m 182\u001b[0m out\u001b[38;5;241m.\u001b[39mprevalence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_prevalence\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:142\u001b[0m, in \u001b[0;36mStandardIpmExecutor.apply\u001b[1;34m(self, world, tick)\u001b[0m\n\u001b[0;32m 139\u001b[0m cohorts \u001b[38;5;241m=\u001b[39m world\u001b[38;5;241m.\u001b[39mget_cohort_array(node)\n\u001b[0;32m 140\u001b[0m effective \u001b[38;5;241m=\u001b[39m cohorts\u001b[38;5;241m.\u001b[39msum(axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, dtype\u001b[38;5;241m=\u001b[39mSimDType)\n\u001b[1;32m--> 142\u001b[0m occurrences \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_events\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43meffective\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 143\u001b[0m cohort_deltas \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_distribute(cohorts, occurrences)\n\u001b[0;32m 144\u001b[0m world\u001b[38;5;241m.\u001b[39mapply_cohort_delta(node, cohort_deltas)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:184\u001b[0m, in \u001b[0;36mStandardIpmExecutor._events\u001b[1;34m(self, node, tick, effective_pop)\u001b[0m\n\u001b[0;32m 182\u001b[0m rate \u001b[38;5;241m=\u001b[39m rate_lambda(rate_args)\n\u001b[0;32m 183\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m (\u001b[38;5;167;01mZeroDivisionError\u001b[39;00m, \u001b[38;5;167;01mFloatingPointError\u001b[39;00m):\n\u001b[1;32m--> 184\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m IpmSimNaNException(\n\u001b[0;32m 185\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_zero_division_args(\n\u001b[0;32m 186\u001b[0m rate_args, node, tick, t)\n\u001b[0;32m 187\u001b[0m )\n\u001b[0;32m 188\u001b[0m \u001b[38;5;66;03m# check for < 0 base, throw error in this case\u001b[39;00m\n\u001b[0;32m 189\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m rate \u001b[38;5;241m<\u001b[39m \u001b[38;5;241m0\u001b[39m:\n", + "\u001b[1;31mIpmSimNaNException\u001b[0m: \nNaN (not a number) rate detected. This is often the result of a divide by zero error.\nWhen constructing the IPM, ensure that no edge transitions can result in division by zero\nThis commonly occurs when defining an S->I edge that is (some rate / sum of the compartments)\nTo fix this, change the edge to define the S->I edge as (some rate / Max(1/sum of the the compartments))\nSee examples of this in the provided example ipm definitions in the data/ipms folder.\n\nShowing current Node : Timestep\n0: 0\n\nShowing current compartment values\nS: 0\nI: 0\nR: 0\n\nShowing current ipm params\nbeta: 0.4\ngamma: 0.2\nxi: 0.011111111111111112\nprob1: 0.5\nprob2: 0.5\n\nShowing current corresponding fork transition\nS->(I, R): I*S*beta*(prob1 + prob2)/(I + R + S)\n\n" + ] + } + ], + "source": [ + "\n", + "from functools import partial\n", + "\n", + "import numpy as np\n", + "\n", + "from epymorph import *\n", + "from epymorph.geo.spec import NO_DURATION, AttribDef, StaticGeoSpec\n", + "from epymorph.geo.static import StaticGeo\n", + "from epymorph.initializer import single_location\n", + "\n", + "my_geo = StaticGeo(\n", + " spec=StaticGeoSpec(\n", + " attributes=[\n", + " AttribDef('label', dtype=str, shape=Shapes.N),\n", + " AttribDef('population', dtype=str, shape=Shapes.N),\n", + " ],\n", + " time_period=NO_DURATION,\n", + " ),\n", + " values={\n", + " 'label': np.array(['a', 'b', 'c']),\n", + " 'population': np.array([0, 10, 20], dtype=np.int64),\n", + " },\n", + ")\n", + "\n", + "\n", + "sim = StandardSimulation(\n", + " geo=my_geo,\n", + " ipm=load_ipm_fork(),\n", + " mm=mm_library['no'](),\n", + " params={\n", + " 'phi': 40.0,\n", + " 'beta': 0.4,\n", + " 'gamma': 1 / 5,\n", + " 'xi': 1 / 90,\n", + " 'prob1': 1 / 2,\n", + " 'prob2': 1 / 2\n", + " },\n", + " time_frame=TimeFrame.of(\"2015-01-01\", 150),\n", + " initializer=partial(single_location, location=1, seed_size=5),\n", + " rng=default_rng(1)\n", + ")\n", + "\n", + "with sim_messaging(sim):\n", + " out = sim.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Case 3:\n", + "\n", + "If we just provide a negative value to `prob1`, which is used in the demo fork ipm, we get a negative probability\n", + "\n", + "To handle this, the exception type `IpmSimInvalidProbsException` is used in `ipm_exec.py` to catch the less than zero error." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running simulation (StandardSimulation):\n", + "• 2015-01-01 to 2015-05-31 (150 days)\n", + "• 6 geo nodes\n", + "| | 0% \r" + ] + }, + { + "ename": "IpmSimInvalidProbsException", + "evalue": "\nInvalid probabilities for fork definition detected. Probabilities for a \ngiven tick should always be nonnegative and sum to 1\n\nShowing current Node : Timestep\n0: 0\n\nShowing current compartment values\nS: 18811310\nI: 0\nR: 0\n\nShowing current ipm params\nbeta: 0.4\ngamma: 0.2\nxi: 0.011111111111111112\nprob1: -0.25\nprob2: 0.2\n\nShowing current corresponding fork transition and probabilities\nS->(I, R): I*S*beta*(prob1 + prob2)/(I + R + S)\nProbabilities: prob1/(prob1 + prob2), prob2/(prob1 + prob2)\n\n", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mIpmSimInvalidProbsException\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[8], line 18\u001b[0m\n\u001b[0;32m 1\u001b[0m sim \u001b[38;5;241m=\u001b[39m StandardSimulation(\n\u001b[0;32m 2\u001b[0m geo\u001b[38;5;241m=\u001b[39mpei_geo,\n\u001b[0;32m 3\u001b[0m ipm\u001b[38;5;241m=\u001b[39mload_ipm_fork(),\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 14\u001b[0m rng\u001b[38;5;241m=\u001b[39mdefault_rng(\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m 15\u001b[0m )\n\u001b[0;32m 17\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m sim_messaging(sim):\n\u001b[1;32m---> 18\u001b[0m out \u001b[38;5;241m=\u001b[39m \u001b[43msim\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\standard_sim.py:180\u001b[0m, in \u001b[0;36mStandardSimulation.run\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 178\u001b[0m \u001b[38;5;66;03m# Then do IPM\u001b[39;00m\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the IPM\u001b[39m\u001b[38;5;124m\"\u001b[39m, IpmSimException, AttributeException):\n\u001b[1;32m--> 180\u001b[0m tick_events, tick_prevalence \u001b[38;5;241m=\u001b[39m \u001b[43mipm_exec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mapply\u001b[49m\u001b[43m(\u001b[49m\u001b[43mworld\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 181\u001b[0m out\u001b[38;5;241m.\u001b[39mincidence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_events\n\u001b[0;32m 182\u001b[0m out\u001b[38;5;241m.\u001b[39mprevalence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_prevalence\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:142\u001b[0m, in \u001b[0;36mStandardIpmExecutor.apply\u001b[1;34m(self, world, tick)\u001b[0m\n\u001b[0;32m 139\u001b[0m cohorts \u001b[38;5;241m=\u001b[39m world\u001b[38;5;241m.\u001b[39mget_cohort_array(node)\n\u001b[0;32m 140\u001b[0m effective \u001b[38;5;241m=\u001b[39m cohorts\u001b[38;5;241m.\u001b[39msum(axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, dtype\u001b[38;5;241m=\u001b[39mSimDType)\n\u001b[1;32m--> 142\u001b[0m occurrences \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_events\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43meffective\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 143\u001b[0m cohort_deltas \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_distribute(cohorts, occurrences)\n\u001b[0;32m 144\u001b[0m world\u001b[38;5;241m.\u001b[39mapply_cohort_delta(node, cohort_deltas)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:197\u001b[0m, in \u001b[0;36mStandardIpmExecutor._events\u001b[1;34m(self, node, tick, effective_pop)\u001b[0m\n\u001b[0;32m 195\u001b[0m \u001b[38;5;66;03m# check for negative probs\u001b[39;00m\n\u001b[0;32m 196\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28many\u001b[39m(n \u001b[38;5;241m<\u001b[39m \u001b[38;5;241m0\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m n \u001b[38;5;129;01min\u001b[39;00m prob):\n\u001b[1;32m--> 197\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m IpmSimInvalidProbsException(\n\u001b[0;32m 198\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_invalid_prob_args(rate_args, node, tick, t)\n\u001b[0;32m 199\u001b[0m )\n\u001b[0;32m 200\u001b[0m stop \u001b[38;5;241m=\u001b[39m index \u001b[38;5;241m+\u001b[39m size\n\u001b[0;32m 201\u001b[0m occur[index:stop] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx\u001b[38;5;241m.\u001b[39mrng\u001b[38;5;241m.\u001b[39mmultinomial(\n\u001b[0;32m 202\u001b[0m base, prob)\n", + "\u001b[1;31mIpmSimInvalidProbsException\u001b[0m: \nInvalid probabilities for fork definition detected. Probabilities for a \ngiven tick should always be nonnegative and sum to 1\n\nShowing current Node : Timestep\n0: 0\n\nShowing current compartment values\nS: 18811310\nI: 0\nR: 0\n\nShowing current ipm params\nbeta: 0.4\ngamma: 0.2\nxi: 0.011111111111111112\nprob1: -0.25\nprob2: 0.2\n\nShowing current corresponding fork transition and probabilities\nS->(I, R): I*S*beta*(prob1 + prob2)/(I + R + S)\nProbabilities: prob1/(prob1 + prob2), prob2/(prob1 + prob2)\n\n" + ] + } + ], + "source": [ + "sim = StandardSimulation(\n", + " geo=pei_geo,\n", + " ipm=load_ipm_fork(),\n", + " mm=mm_library['no'](),\n", + " params={\n", + " 'beta': 0.4,\n", + " 'gamma': 1 / 5,\n", + " 'xi': 1 / 90,\n", + " 'prob1': -1 / 4, # WATCHOUT: negative value\n", + " 'prob2': 1 / 5\n", + " },\n", + " time_frame=TimeFrame.of(\"2015-01-01\", 150),\n", + " initializer=partial(single_location, location=1, seed_size=5),\n", + " rng=default_rng(1)\n", + ")\n", + "\n", + "with sim_messaging(sim):\n", + " out = sim.run()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "ename": "IpmSimInvalidProbsException", + "evalue": "\nInvalid probabilities for fork definition detected. Probabilities for a \ngiven tick should always be nonnegative and sum to 1\n\nShowing current Node : Timestep\n0: 0\n\nShowing current compartment values\nS: 18811310\nI: 0\nR: 0\nH: 0\n\nShowing current ipm params\nbeta: 0.4\ngamma: 0.2\nxi: 0.01\nhospitalization_prob: -0.2\nhospitalization_duration: 15.0\n\nShowing current corresponding fork transition and probabilities\nI->(H, R): I*gamma\nProbabilities: hospitalization_prob, 1 - hospitalization_prob\n\n", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mIpmSimInvalidProbsException\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[9], line 17\u001b[0m\n\u001b[0;32m 1\u001b[0m sim \u001b[38;5;241m=\u001b[39m StandardSimulation(\n\u001b[0;32m 2\u001b[0m geo\u001b[38;5;241m=\u001b[39mgeo_library[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpei\u001b[39m\u001b[38;5;124m'\u001b[39m](),\n\u001b[0;32m 3\u001b[0m ipm\u001b[38;5;241m=\u001b[39mipm_library[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msirh\u001b[39m\u001b[38;5;124m'\u001b[39m](),\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 14\u001b[0m rng\u001b[38;5;241m=\u001b[39mdefault_rng(\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m 15\u001b[0m )\n\u001b[1;32m---> 17\u001b[0m \u001b[43msim\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\standard_sim.py:180\u001b[0m, in \u001b[0;36mStandardSimulation.run\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 178\u001b[0m \u001b[38;5;66;03m# Then do IPM\u001b[39;00m\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the IPM\u001b[39m\u001b[38;5;124m\"\u001b[39m, IpmSimException, AttributeException):\n\u001b[1;32m--> 180\u001b[0m tick_events, tick_prevalence \u001b[38;5;241m=\u001b[39m \u001b[43mipm_exec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mapply\u001b[49m\u001b[43m(\u001b[49m\u001b[43mworld\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 181\u001b[0m out\u001b[38;5;241m.\u001b[39mincidence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_events\n\u001b[0;32m 182\u001b[0m out\u001b[38;5;241m.\u001b[39mprevalence[tick\u001b[38;5;241m.\u001b[39mindex] \u001b[38;5;241m=\u001b[39m tick_prevalence\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:142\u001b[0m, in \u001b[0;36mStandardIpmExecutor.apply\u001b[1;34m(self, world, tick)\u001b[0m\n\u001b[0;32m 139\u001b[0m cohorts \u001b[38;5;241m=\u001b[39m world\u001b[38;5;241m.\u001b[39mget_cohort_array(node)\n\u001b[0;32m 140\u001b[0m effective \u001b[38;5;241m=\u001b[39m cohorts\u001b[38;5;241m.\u001b[39msum(axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, dtype\u001b[38;5;241m=\u001b[39mSimDType)\n\u001b[1;32m--> 142\u001b[0m occurrences \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_events\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43meffective\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 143\u001b[0m cohort_deltas \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_distribute(cohorts, occurrences)\n\u001b[0;32m 144\u001b[0m world\u001b[38;5;241m.\u001b[39mapply_cohort_delta(node, cohort_deltas)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\ipm_exec.py:197\u001b[0m, in \u001b[0;36mStandardIpmExecutor._events\u001b[1;34m(self, node, tick, effective_pop)\u001b[0m\n\u001b[0;32m 195\u001b[0m \u001b[38;5;66;03m# check for negative probs\u001b[39;00m\n\u001b[0;32m 196\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28many\u001b[39m(n \u001b[38;5;241m<\u001b[39m \u001b[38;5;241m0\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m n \u001b[38;5;129;01min\u001b[39;00m prob):\n\u001b[1;32m--> 197\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m IpmSimInvalidProbsException(\n\u001b[0;32m 198\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_invalid_prob_args(rate_args, node, tick, t)\n\u001b[0;32m 199\u001b[0m )\n\u001b[0;32m 200\u001b[0m stop \u001b[38;5;241m=\u001b[39m index \u001b[38;5;241m+\u001b[39m size\n\u001b[0;32m 201\u001b[0m occur[index:stop] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx\u001b[38;5;241m.\u001b[39mrng\u001b[38;5;241m.\u001b[39mmultinomial(\n\u001b[0;32m 202\u001b[0m base, prob)\n", + "\u001b[1;31mIpmSimInvalidProbsException\u001b[0m: \nInvalid probabilities for fork definition detected. Probabilities for a \ngiven tick should always be nonnegative and sum to 1\n\nShowing current Node : Timestep\n0: 0\n\nShowing current compartment values\nS: 18811310\nI: 0\nR: 0\nH: 0\n\nShowing current ipm params\nbeta: 0.4\ngamma: 0.2\nxi: 0.01\nhospitalization_prob: -0.2\nhospitalization_duration: 15.0\n\nShowing current corresponding fork transition and probabilities\nI->(H, R): I*gamma\nProbabilities: hospitalization_prob, 1 - hospitalization_prob\n\n" + ] + } + ], + "source": [ + "\n", + "sim = StandardSimulation(\n", + " geo=geo_library['pei'](),\n", + " ipm=ipm_library['sirh'](),\n", + " mm=mm_library['no'](),\n", + " params={\n", + " 'beta': 0.4,\n", + " 'gamma': 1 / 5,\n", + " 'xi': 1 / 100,\n", + " 'hospitalization_prob': -1 / 5,\n", + " 'hospitalization_duration': 15\n", + " },\n", + " time_frame=TimeFrame.of(\"2015-01-01\", 150),\n", + " initializer=partial(single_location, location=1, seed_size=5),\n", + " rng=default_rng(1)\n", + ")\n", + "\n", + "sim.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Case 4:\n", + "\n", + "Issues can arise when an incompatible movement model parameter and clause combination are in a given simulation, resulting in\n", + "simulation errors. An example of this is in the `pei` movement model, where `move_control` should be between 0 and 1\n", + "\n", + "This handling is a little more general because it's difficult to tell whether the error results from an incorrect parameter, clause, or both" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running simulation (StandardSimulation):\n", + "• 2015-01-01 to 2015-05-31 (150 days)\n", + "• 6 geo nodes\n", + "| | 0% \r" + ] + }, + { + "ename": "MmSimException", + "evalue": "Error from applying clause 'commuters': see exception trace", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\movement\\movement_model.py:113\u001b[0m, in \u001b[0;36mDynamicTravelClause.requested\u001b[1;34m(self, ctx, predef, tick)\u001b[0m\n\u001b[0;32m 112\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m--> 113\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_requested\u001b[49m\u001b[43m(\u001b[49m\u001b[43mctx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpredef\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 114\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m e:\n\u001b[0;32m 115\u001b[0m \u001b[38;5;66;03m# NOTE: catching KeyError here will be necessary (to get nice error messages)\u001b[39;00m\n\u001b[0;32m 116\u001b[0m \u001b[38;5;66;03m# until we can properly validate the MM clauses.\u001b[39;00m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\movement\\compile.py:165\u001b[0m, in \u001b[0;36m_adapt_move_function..fn_arity1\u001b[1;34m(ctx, predef, tick)\u001b[0m\n\u001b[0;32m 163\u001b[0m \u001b[38;5;129m@wraps\u001b[39m(fn)\n\u001b[0;32m 164\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfn_arity1\u001b[39m(ctx: MovementContext, predef: PredefParams, tick: Tick) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m NDArray[SimDType]:\n\u001b[1;32m--> 165\u001b[0m requested \u001b[38;5;241m=\u001b[39m \u001b[43mfn\u001b[49m\u001b[43m(\u001b[49m\u001b[43mctx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpredef\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 166\u001b[0m np\u001b[38;5;241m.\u001b[39mfill_diagonal(requested, \u001b[38;5;241m0\u001b[39m)\n", + "File \u001b[1;32m:3\u001b[0m, in \u001b[0;36mcommuters\u001b[1;34m(ctx, predef, t)\u001b[0m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\movement\\compile.py:77\u001b[0m, in \u001b[0;36m_movement_global_namespace..as_simdtype..wrapped_func\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m 75\u001b[0m \u001b[38;5;129m@wraps\u001b[39m(func)\n\u001b[0;32m 76\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mwrapped_func\u001b[39m(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m---> 77\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 78\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m np\u001b[38;5;241m.\u001b[39misscalar(result):\n", + "File \u001b[1;32mnumpy\\\\random\\\\_generator.pyx:2996\u001b[0m, in \u001b[0;36mnumpy.random._generator.Generator.binomial\u001b[1;34m()\u001b[0m\n", + "File \u001b[1;32m_common.pyx:391\u001b[0m, in \u001b[0;36mnumpy.random._common.check_array_constraint\u001b[1;34m()\u001b[0m\n", + "File \u001b[1;32m_common.pyx:377\u001b[0m, in \u001b[0;36mnumpy.random._common._check_array_cons_bounded_0_1\u001b[1;34m()\u001b[0m\n", + "\u001b[1;31mValueError\u001b[0m: p < 0, p > 1 or p contains NaNs", + "\nThe above exception was the direct cause of the following exception:\n", + "\u001b[1;31mMmSimException\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[10], line 23\u001b[0m\n\u001b[0;32m 6\u001b[0m sim \u001b[38;5;241m=\u001b[39m StandardSimulation(\n\u001b[0;32m 7\u001b[0m geo\u001b[38;5;241m=\u001b[39mgeo_library[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpei\u001b[39m\u001b[38;5;124m'\u001b[39m](),\n\u001b[0;32m 8\u001b[0m ipm\u001b[38;5;241m=\u001b[39mipm_library[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpei\u001b[39m\u001b[38;5;124m'\u001b[39m](),\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 19\u001b[0m rng\u001b[38;5;241m=\u001b[39mdefault_rng(\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m 20\u001b[0m )\n\u001b[0;32m 22\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m sim_messaging(sim):\n\u001b[1;32m---> 23\u001b[0m out \u001b[38;5;241m=\u001b[39m \u001b[43msim\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\standard_sim.py:176\u001b[0m, in \u001b[0;36mStandardSimulation.run\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 173\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m tick \u001b[38;5;129;01min\u001b[39;00m ctx\u001b[38;5;241m.\u001b[39mclock():\n\u001b[0;32m 174\u001b[0m \u001b[38;5;66;03m# First do movement\u001b[39;00m\n\u001b[0;32m 175\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the movement model\u001b[39m\u001b[38;5;124m\"\u001b[39m, MmSimException, AttributeException):\n\u001b[1;32m--> 176\u001b[0m \u001b[43mmovement_exec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mapply\u001b[49m\u001b[43m(\u001b[49m\u001b[43mworld\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 178\u001b[0m \u001b[38;5;66;03m# Then do IPM\u001b[39;00m\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the IPM\u001b[39m\u001b[38;5;124m\"\u001b[39m, IpmSimException, AttributeException):\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\mm_exec.py:92\u001b[0m, in \u001b[0;36mStandardMovementExecutor.apply\u001b[1;34m(self, world, tick)\u001b[0m\n\u001b[0;32m 90\u001b[0m \u001b[38;5;28;01mcontinue\u001b[39;00m\n\u001b[0;32m 91\u001b[0m local_array \u001b[38;5;241m=\u001b[39m world\u001b[38;5;241m.\u001b[39mget_local_array()\n\u001b[1;32m---> 92\u001b[0m travelers \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_travelers\u001b[49m\u001b[43m(\u001b[49m\u001b[43mclause\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlocal_array\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 93\u001b[0m returns \u001b[38;5;241m=\u001b[39m clause\u001b[38;5;241m.\u001b[39mreturns(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx, tick)\n\u001b[0;32m 94\u001b[0m return_tick \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx\u001b[38;5;241m.\u001b[39mresolve_tick(tick, returns)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\mm_exec.py:130\u001b[0m, in \u001b[0;36mStandardMovementExecutor._travelers\u001b[1;34m(self, clause, tick, local_cohorts)\u001b[0m\n\u001b[0;32m 122\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m 123\u001b[0m \u001b[38;5;124;03mCalculate the number of travelers resulting from this movement clause for this tick.\u001b[39;00m\n\u001b[0;32m 124\u001b[0m \u001b[38;5;124;03mThis evaluates the requested number movers, modulates that based on the available movers,\u001b[39;00m\n\u001b[0;32m 125\u001b[0m \u001b[38;5;124;03mthen selects exactly which individuals (by compartment) should move.\u001b[39;00m\n\u001b[0;32m 126\u001b[0m \u001b[38;5;124;03mReturns an (N,N,C) array; from-source-to-destination-by-compartment.\u001b[39;00m\n\u001b[0;32m 127\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m 128\u001b[0m _, N, C, _ \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx\u001b[38;5;241m.\u001b[39mdim\u001b[38;5;241m.\u001b[39mTNCE\n\u001b[1;32m--> 130\u001b[0m clause_movers \u001b[38;5;241m=\u001b[39m \u001b[43mclause\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrequested\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_ctx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_predef\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 131\u001b[0m np\u001b[38;5;241m.\u001b[39mfill_diagonal(clause_movers, \u001b[38;5;241m0\u001b[39m)\n\u001b[0;32m 132\u001b[0m clause_sum \u001b[38;5;241m=\u001b[39m clause_movers\u001b[38;5;241m.\u001b[39msum(axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m, dtype\u001b[38;5;241m=\u001b[39mSimDType)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\movement\\movement_model.py:123\u001b[0m, in \u001b[0;36mDynamicTravelClause.requested\u001b[1;34m(self, ctx, predef, tick)\u001b[0m\n\u001b[0;32m 119\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m e:\n\u001b[0;32m 120\u001b[0m \u001b[38;5;66;03m# NOTE: catching exceptions here is necessary to get nice error messages\u001b[39;00m\n\u001b[0;32m 121\u001b[0m \u001b[38;5;66;03m# for some value error cause by incorrect parameter and/or clause definition\u001b[39;00m\n\u001b[0;32m 122\u001b[0m msg \u001b[38;5;241m=\u001b[39m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mError from applying clause \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mname\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m: see exception trace\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m--> 123\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m MmSimException(msg) \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01me\u001b[39;00m\n", + "\u001b[1;31mMmSimException\u001b[0m: Error from applying clause 'commuters': see exception trace" + ] + } + ], + "source": [ + "from functools import partial\n", + "\n", + "from epymorph import *\n", + "from epymorph.initializer import single_location\n", + "\n", + "sim = StandardSimulation(\n", + " geo=geo_library['pei'](),\n", + " ipm=ipm_library['pei'](),\n", + " mm=mm_library['pei'](),\n", + " params={\n", + " 'infection_duration': 40.0,\n", + " 'immunity_duration': 0.4,\n", + " 'humidity': 20.2,\n", + " 'move_control': 100, # NOTICE: Invalid move_control value\n", + " 'theta': 5\n", + " },\n", + " time_frame=TimeFrame.of(\"2015-01-01\", 150),\n", + " initializer=partial(single_location, location=1, seed_size=5),\n", + " rng=default_rng(1)\n", + ")\n", + "\n", + "with sim_messaging(sim):\n", + " out = sim.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Other example: Negative theta" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running simulation (StandardSimulation):\n", + "• 2015-01-01 to 2015-05-31 (150 days)\n", + "• 6 geo nodes\n", + "| | 0% \r" + ] + }, + { + "ename": "MmSimException", + "evalue": "Error from applying clause 'dispersers': see exception trace", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\movement\\movement_model.py:113\u001b[0m, in \u001b[0;36mDynamicTravelClause.requested\u001b[1;34m(self, ctx, predef, tick)\u001b[0m\n\u001b[0;32m 112\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m--> 113\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_requested\u001b[49m\u001b[43m(\u001b[49m\u001b[43mctx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpredef\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 114\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m e:\n\u001b[0;32m 115\u001b[0m \u001b[38;5;66;03m# NOTE: catching KeyError here will be necessary (to get nice error messages)\u001b[39;00m\n\u001b[0;32m 116\u001b[0m \u001b[38;5;66;03m# until we can properly validate the MM clauses.\u001b[39;00m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\movement\\compile.py:165\u001b[0m, in \u001b[0;36m_adapt_move_function..fn_arity1\u001b[1;34m(ctx, predef, tick)\u001b[0m\n\u001b[0;32m 163\u001b[0m \u001b[38;5;129m@wraps\u001b[39m(fn)\n\u001b[0;32m 164\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfn_arity1\u001b[39m(ctx: MovementContext, predef: PredefParams, tick: Tick) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m NDArray[SimDType]:\n\u001b[1;32m--> 165\u001b[0m requested \u001b[38;5;241m=\u001b[39m \u001b[43mfn\u001b[49m\u001b[43m(\u001b[49m\u001b[43mctx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpredef\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 166\u001b[0m np\u001b[38;5;241m.\u001b[39mfill_diagonal(requested, \u001b[38;5;241m0\u001b[39m)\n", + "File \u001b[1;32m:3\u001b[0m, in \u001b[0;36mdispersers\u001b[1;34m(ctx, predef, t)\u001b[0m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\movement\\compile.py:77\u001b[0m, in \u001b[0;36m_movement_global_namespace..as_simdtype..wrapped_func\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m 75\u001b[0m \u001b[38;5;129m@wraps\u001b[39m(func)\n\u001b[0;32m 76\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mwrapped_func\u001b[39m(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m---> 77\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 78\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m np\u001b[38;5;241m.\u001b[39misscalar(result):\n", + "File \u001b[1;32mnumpy\\\\random\\\\_generator.pyx:3230\u001b[0m, in \u001b[0;36mnumpy.random._generator.Generator.poisson\u001b[1;34m()\u001b[0m\n", + "File \u001b[1;32m_common.pyx:883\u001b[0m, in \u001b[0;36mnumpy.random._common.disc\u001b[1;34m()\u001b[0m\n", + "File \u001b[1;32m_common.pyx:680\u001b[0m, in \u001b[0;36mnumpy.random._common.discrete_broadcast_d\u001b[1;34m()\u001b[0m\n", + "File \u001b[1;32m_common.pyx:408\u001b[0m, in \u001b[0;36mnumpy.random._common.check_array_constraint\u001b[1;34m()\u001b[0m\n", + "\u001b[1;31mValueError\u001b[0m: lam < 0 or lam contains NaNs", + "\nThe above exception was the direct cause of the following exception:\n", + "\u001b[1;31mMmSimException\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[11], line 23\u001b[0m\n\u001b[0;32m 6\u001b[0m sim \u001b[38;5;241m=\u001b[39m StandardSimulation(\n\u001b[0;32m 7\u001b[0m geo\u001b[38;5;241m=\u001b[39mgeo_library[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpei\u001b[39m\u001b[38;5;124m'\u001b[39m](),\n\u001b[0;32m 8\u001b[0m ipm\u001b[38;5;241m=\u001b[39mipm_library[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpei\u001b[39m\u001b[38;5;124m'\u001b[39m](),\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 19\u001b[0m rng\u001b[38;5;241m=\u001b[39mdefault_rng(\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m 20\u001b[0m )\n\u001b[0;32m 22\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m sim_messaging(sim):\n\u001b[1;32m---> 23\u001b[0m out \u001b[38;5;241m=\u001b[39m \u001b[43msim\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\standard_sim.py:176\u001b[0m, in \u001b[0;36mStandardSimulation.run\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 173\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m tick \u001b[38;5;129;01min\u001b[39;00m ctx\u001b[38;5;241m.\u001b[39mclock():\n\u001b[0;32m 174\u001b[0m \u001b[38;5;66;03m# First do movement\u001b[39;00m\n\u001b[0;32m 175\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the movement model\u001b[39m\u001b[38;5;124m\"\u001b[39m, MmSimException, AttributeException):\n\u001b[1;32m--> 176\u001b[0m \u001b[43mmovement_exec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mapply\u001b[49m\u001b[43m(\u001b[49m\u001b[43mworld\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 178\u001b[0m \u001b[38;5;66;03m# Then do IPM\u001b[39;00m\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mexecuting the IPM\u001b[39m\u001b[38;5;124m\"\u001b[39m, IpmSimException, AttributeException):\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\mm_exec.py:92\u001b[0m, in \u001b[0;36mStandardMovementExecutor.apply\u001b[1;34m(self, world, tick)\u001b[0m\n\u001b[0;32m 90\u001b[0m \u001b[38;5;28;01mcontinue\u001b[39;00m\n\u001b[0;32m 91\u001b[0m local_array \u001b[38;5;241m=\u001b[39m world\u001b[38;5;241m.\u001b[39mget_local_array()\n\u001b[1;32m---> 92\u001b[0m travelers \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_travelers\u001b[49m\u001b[43m(\u001b[49m\u001b[43mclause\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlocal_array\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 93\u001b[0m returns \u001b[38;5;241m=\u001b[39m clause\u001b[38;5;241m.\u001b[39mreturns(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx, tick)\n\u001b[0;32m 94\u001b[0m return_tick \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx\u001b[38;5;241m.\u001b[39mresolve_tick(tick, returns)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\engine\\mm_exec.py:130\u001b[0m, in \u001b[0;36mStandardMovementExecutor._travelers\u001b[1;34m(self, clause, tick, local_cohorts)\u001b[0m\n\u001b[0;32m 122\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m 123\u001b[0m \u001b[38;5;124;03mCalculate the number of travelers resulting from this movement clause for this tick.\u001b[39;00m\n\u001b[0;32m 124\u001b[0m \u001b[38;5;124;03mThis evaluates the requested number movers, modulates that based on the available movers,\u001b[39;00m\n\u001b[0;32m 125\u001b[0m \u001b[38;5;124;03mthen selects exactly which individuals (by compartment) should move.\u001b[39;00m\n\u001b[0;32m 126\u001b[0m \u001b[38;5;124;03mReturns an (N,N,C) array; from-source-to-destination-by-compartment.\u001b[39;00m\n\u001b[0;32m 127\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m 128\u001b[0m _, N, C, _ \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ctx\u001b[38;5;241m.\u001b[39mdim\u001b[38;5;241m.\u001b[39mTNCE\n\u001b[1;32m--> 130\u001b[0m clause_movers \u001b[38;5;241m=\u001b[39m \u001b[43mclause\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrequested\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_ctx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_predef\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtick\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 131\u001b[0m np\u001b[38;5;241m.\u001b[39mfill_diagonal(clause_movers, \u001b[38;5;241m0\u001b[39m)\n\u001b[0;32m 132\u001b[0m clause_sum \u001b[38;5;241m=\u001b[39m clause_movers\u001b[38;5;241m.\u001b[39msum(axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m, dtype\u001b[38;5;241m=\u001b[39mSimDType)\n", + "File \u001b[1;32mc:\\Users\\izaac\\Documents\\EpiMoRPH\\Epymorph\\epymorph\\movement\\movement_model.py:123\u001b[0m, in \u001b[0;36mDynamicTravelClause.requested\u001b[1;34m(self, ctx, predef, tick)\u001b[0m\n\u001b[0;32m 119\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m e:\n\u001b[0;32m 120\u001b[0m \u001b[38;5;66;03m# NOTE: catching exceptions here is necessary to get nice error messages\u001b[39;00m\n\u001b[0;32m 121\u001b[0m \u001b[38;5;66;03m# for some value error cause by incorrect parameter and/or clause definition\u001b[39;00m\n\u001b[0;32m 122\u001b[0m msg \u001b[38;5;241m=\u001b[39m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mError from applying clause \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mname\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m: see exception trace\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m--> 123\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m MmSimException(msg) \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01me\u001b[39;00m\n", + "\u001b[1;31mMmSimException\u001b[0m: Error from applying clause 'dispersers': see exception trace" + ] + } + ], + "source": [ + "\n", + "from functools import partial\n", + "\n", + "from epymorph import *\n", + "from epymorph.initializer import single_location\n", + "\n", + "sim = StandardSimulation(\n", + " geo=geo_library['pei'](),\n", + " ipm=ipm_library['pei'](),\n", + " mm=mm_library['pei'](),\n", + " params={\n", + " 'infection_duration': 40.0,\n", + " 'immunity_duration': 0.4,\n", + " 'humidity': 20.2,\n", + " 'move_control': 0.2,\n", + " 'theta': -5\n", + " },\n", + " time_frame=TimeFrame.of(\"2015-01-01\", 150),\n", + " initializer=partial(single_location, location=1, seed_size=5),\n", + " rng=default_rng(1)\n", + ")\n", + "\n", + "with sim_messaging(sim):\n", + " out = sim.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Extending this work\n", + "Overall, this analysis has resulted in error handling for common ipm simulation errors. However, if more are to arise in the future, they can be expanded in the `epymorph/error.py` module (if the error is related to ipms)\n", + "by inhereting from `IpmSimException`, or `IpmSimExceptionWithFields` if simulation values are to be printed along with the error. See the definition of `IpmSimNanException` in `epymorph/error.py` for an example.\n", + "\n", + "For errors not related to IPMs, note all classes defined in `error.py` that inherit from `SimulationException` as potential points of extension for future simulation-based errors." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "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.9" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/doc/devlog/README.md b/doc/devlog/README.md index d1bb438c..b12aca3f 100644 --- a/doc/devlog/README.md +++ b/doc/devlog/README.md @@ -47,6 +47,7 @@ This folder is a handy place to put Jupyter notebooks or other documents which h | 2024-03-01.ipynb | Tyler | | Getting the indices of IPM events and compartments by name with wildcard support. | | 2024-03-13.ipynb | Tyler | | Showing off movement data collection (NEW!) | | 2024-04-04-draw-demo.ipynb | Izaac | | Showing the new draw module for visualising IPMs (NEW!) | +| 2024-04-16.ipynb | Izaac | | Showing error handling for common ipm errors (NEW!)| ## Contributing diff --git a/epymorph/__init__.py b/epymorph/__init__.py index dcd342bf..bc1bb14c 100644 --- a/epymorph/__init__.py +++ b/epymorph/__init__.py @@ -1,5 +1,7 @@ """epymorph's main package and main exports""" +from numpy import seterr + import epymorph.compartment_model as IPM from epymorph.data import geo_library, ipm_library, mm_library from epymorph.data_shape import Shapes @@ -10,6 +12,10 @@ from epymorph.proxy import dim, geo from epymorph.simulation import SimDType, TimeFrame, default_rng +# set numpy errors to raise exceptions instead of warnings, useful for catching +# simulation errrors +seterr(all='raise') + __all__ = [ 'IPM', 'ipm_library', diff --git a/epymorph/engine/ipm_exec.py b/epymorph/engine/ipm_exec.py index 0f0ff09b..2bf6e027 100644 --- a/epymorph/engine/ipm_exec.py +++ b/epymorph/engine/ipm_exec.py @@ -12,6 +12,8 @@ TransitionDef, exogenous_states) from epymorph.engine.context import RumeContext, Tick from epymorph.engine.world import World +from epymorph.error import (IpmSimInvalidProbsException, + IpmSimLessThanZeroException, IpmSimNaNException) from epymorph.simulation import SimDType from epymorph.sympy_shim import SympyLambda, lambdify, lambdify_list from epymorph.util import index_of @@ -160,12 +162,41 @@ def _events(self, node: int, tick: Tick, effective_pop: NDArray[SimDType]) -> ND for t in self._trxs: match t: case _IndependentTrx(rate_lambda): - rate = rate_lambda(rate_args) + # get rate from lambda expression, catch divide by zero error + try: + rate = rate_lambda(rate_args) + except (ZeroDivisionError, FloatingPointError): + raise IpmSimNaNException( + self._get_zero_division_args( + rate_args, node, tick, t) + ) + # check for < 0 rate, throw error in this case + if rate < 0: + raise IpmSimLessThanZeroException( + self._get_default_error_args(rate_args, node, tick) + ) occur[index] = self._ctx.rng.poisson(rate * tick.tau) case _ForkedTrx(size, rate_lambda, prob_lambda): - rate = rate_lambda(rate_args) + # get rate from lambda expression, catch divide by zero error + try: + rate = rate_lambda(rate_args) + except (ZeroDivisionError, FloatingPointError): + raise IpmSimNaNException( + self._get_zero_division_args( + rate_args, node, tick, t) + ) + # check for < 0 base, throw error in this case + if rate < 0: + raise IpmSimLessThanZeroException( + self._get_default_error_args(rate_args, node, tick) + ) base = self._ctx.rng.poisson(rate * tick.tau) prob = prob_lambda(rate_args) + # check for negative probs + if any(n < 0 for n in prob): + raise IpmSimInvalidProbsException( + self._get_invalid_prob_args(rate_args, node, tick, t) + ) stop = index + size occur[index:stop] = self._ctx.rng.multinomial( base, prob) @@ -199,6 +230,56 @@ def _events(self, node: int, tick: Tick, effective_pop: NDArray[SimDType]) -> ND desired, available) return occur + def _get_default_error_args(self, rate_attrs: list, node: int, tick: Tick) -> list[tuple[str, dict]]: + arg_list = [] + arg_list.append(("Node : Timestep", {node: tick.step})) + arg_list.append(("compartment values", { + name: value for (name, value) in zip(self._ctx.ipm.compartment_names, + rate_attrs[:self._ctx.dim.compartments]) + })) + arg_list.append(("ipm params", { + attribute.name: value for (attribute, value) in zip(self._ctx.ipm.attributes, + rate_attrs[self._ctx.dim.compartments:]) + })) + + return arg_list + + def _get_invalid_prob_args(self, rate_attrs: list, node: int, tick: Tick, + transition: _ForkedTrx) -> list[tuple[str, dict]]: + arg_list = self._get_default_error_args(rate_attrs, node, tick) + + transition_index = self._trxs.index(transition) + corr_transition = self._ctx.ipm.transitions[transition_index] + if isinstance(corr_transition, ForkDef): + to_compartments = ", ".join([str(edge.compartment_to) + for edge in corr_transition.edges]) + from_compartment = corr_transition.edges[0].compartment_from + arg_list.append(("corresponding fork transition and probabilities", + { + f"{from_compartment}->({to_compartments})": corr_transition.rate, + f"Probabilities": ', '.join([str(expr) for expr in corr_transition.probs]), + })) + + return arg_list + + def _get_zero_division_args(self, rate_attrs: list, node: int, tick: Tick, + transition: _IndependentTrx | _ForkedTrx) -> list[tuple[str, dict]]: + arg_list = self._get_default_error_args(rate_attrs, node, tick) + + transition_index = self._trxs.index(transition) + corr_transition = self._ctx.ipm.transitions[transition_index] + if isinstance(corr_transition, EdgeDef): + arg_list.append(("corresponding transition", { + f"{corr_transition.compartment_from}->{corr_transition.compartment_to}": corr_transition.rate})) + if isinstance(corr_transition, ForkDef): + to_compartments = ", ".join([str(edge.compartment_to) + for edge in corr_transition.edges]) + from_compartment = corr_transition.edges[0].compartment_from + arg_list.append(("corresponding fork transition", { + f"{from_compartment}->({to_compartments})": corr_transition.rate})) + + return arg_list + def _distribute(self, cohorts: NDArray[SimDType], events: NDArray[SimDType]) -> NDArray[SimDType]: """Distribute all events across a location's cohorts and return the compartment deltas for each.""" x = cohorts.shape[0] diff --git a/epymorph/error.py b/epymorph/error.py index 43d5592c..53d0ce6b 100644 --- a/epymorph/error.py +++ b/epymorph/error.py @@ -2,6 +2,7 @@ A common exception framework for epymorph. """ from contextlib import contextmanager +from textwrap import dedent from typing import Self @@ -99,6 +100,74 @@ class IpmSimException(SimulationException): """Exception during IPM processing.""" +class IpmSimExceptionWithFields(IpmSimException): + """ + Exception during IPM processing where it is appropriate to show specific + fields within the simulation. + To create a new error with fields, create a subclass of this and set the + displayed error message along with the fields to print. + See 'IpmSimNaNException' for an example. + """ + + display_fields: tuple[str, dict] | list[tuple[str, dict]] + + def __init__(self, message: str, display_fields: tuple[str, dict] | list[tuple[str, dict]]): + super().__init__(message) + if isinstance(display_fields, tuple): + display_fields = [display_fields] + self.display_fields = display_fields + + def __str__(self): + msg = super().__str__() + fields = "" + for name, values in self.display_fields: + fields += f"Showing current {name}\n" + for key, value in values.items(): + fields += f"{key}: {value}\n" + fields += "\n" + return f"{msg}\n{fields}" + + +class IpmSimNaNException(IpmSimExceptionWithFields): + """Exception for handling NaN (not a number) rate values""" + + def __init__(self, display_fields: tuple[str, dict] | list[tuple[str, dict]]): + msg = ''' + NaN (not a number) rate detected. This is often the result of a divide by zero error. + When constructing the IPM, ensure that no edge transitions can result in division by zero + This commonly occurs when defining an S->I edge that is (some rate / sum of the compartments) + To fix this, change the edge to define the S->I edge as (some rate / Max(1/sum of the the compartments)) + See examples of this in the provided example ipm definitions in the data/ipms folder. + ''' + msg = dedent(msg) + super().__init__(msg, display_fields) + + +class IpmSimLessThanZeroException(IpmSimExceptionWithFields): + """ Exception for handling less than 0 rate values """ + + def __init__(self, display_fields: tuple[str, dict] | list[tuple[str, dict]]): + msg = ''' + Less than zero rate detected. When providing or defining ipm parameters, ensure that + they will not result in a negative rate. Note: this can often happen unintentionally + if a function is given as a parameter. + ''' + msg = dedent(msg) + super().__init__(msg, display_fields) + + +class IpmSimInvalidProbsException(IpmSimExceptionWithFields): + """ Exception for handling invalid probability values """ + + def __init__(self, display_fields: tuple[str, dict] | list[tuple[str, dict]]): + msg = ''' + Invalid probabilities for fork definition detected. Probabilities for a + given tick should always be nonnegative and sum to 1 + ''' + msg = dedent(msg) + super().__init__(msg, display_fields) + + class MmSimException(SimulationException): """Exception during MM processing.""" diff --git a/epymorph/movement/movement_model.py b/epymorph/movement/movement_model.py index a70fdeff..e4efb989 100644 --- a/epymorph/movement/movement_model.py +++ b/epymorph/movement/movement_model.py @@ -11,7 +11,8 @@ from numpy.typing import NDArray from epymorph.compartment_model import CompartmentModel -from epymorph.error import AttributeException, MmValidationException +from epymorph.error import (AttributeException, MmSimException, + MmValidationException) from epymorph.geo.geo import Geo from epymorph.movement.parser import Attribute as MmAttribute from epymorph.params import ContextParams @@ -115,6 +116,11 @@ def requested(self, ctx: MovementContext, predef: PredefParams, tick: Tick) -> N # until we can properly validate the MM clauses. msg = f"Missing attribute {e} required by movement model clause '{self.name}'." raise AttributeException(msg) from None + except Exception as e: + # NOTE: catching exceptions here is necessary to get nice error messages + # for some value error cause by incorrect parameter and/or clause definition + msg = f"Error from applying clause '{self.name}': see exception trace" + raise MmSimException(msg) from e def returns(self, ctx: MovementContext, tick: Tick) -> TickDelta: return self._returns(ctx, tick) diff --git a/epymorph/test/simulate_test.py b/epymorph/test/simulate_test.py index ede1093c..21434f27 100644 --- a/epymorph/test/simulate_test.py +++ b/epymorph/test/simulate_test.py @@ -5,6 +5,14 @@ import numpy as np from epymorph import * +from epymorph.compartment_model import (CompartmentModel, compartment, + create_model, create_symbols, edge, + param) +from epymorph.error import (IpmSimInvalidProbsException, + IpmSimLessThanZeroException, IpmSimNaNException, + MmSimException) +from epymorph.geo.spec import NO_DURATION, AttribDef, StaticGeoSpec +from epymorph.geo.static import StaticGeo from epymorph.initializer import single_location @@ -70,3 +78,174 @@ def test_pei(self): out2.prevalence, "Running the sim twice with the same RNG should yield the same prevalence." ) + + def test_less_than_zero_err(self): + """Test exception handling for a negative rate value due to a negative parameter""" + sim = StandardSimulation( + geo=geo_library['pei'](), + ipm=ipm_library['pei'](), + mm=mm_library['pei'](), + params={ + 'infection_duration': 1 / 4, + 'immunity_duration': -1 / 100, # notice the negative parameter + 'move_control': 0.9, + 'theta': 0.1, + }, + time_frame=TimeFrame.of("2015-01-01", 10), + initializer=partial(single_location, location=0, seed_size=10_000), + rng=default_rng(42) + ) + + with self.assertRaises(IpmSimLessThanZeroException) as e: + sim.run() + + err_msg = str(e.exception) + + self.assertIn("Less than zero rate detected", err_msg) + self.assertIn("Showing current Node : Timestep", err_msg) + self.assertIn("S: ", err_msg) + self.assertIn("I: ", err_msg) + self.assertIn("R: ", err_msg) + self.assertIn("infection_duration: 0.25", err_msg) + self.assertIn("immunity_duration: -0.01", err_msg) + self.assertIn("humidity: 0.01003", err_msg) + + def test_divide_by_zero_err(self): + """Test exception handling for a divide by zero (NaN) error""" + def load_ipm() -> CompartmentModel: + """Load the 'sirs' IPM.""" + symbols = create_symbols( + compartments=[ + compartment('S'), + compartment('I'), + compartment('R'), + ], + attributes=[ + param('beta', shape=Shapes.TxN), # infectivity + # progression from infected to recovered + param('gamma', shape=Shapes.TxN), + # progression from recovered to susceptible + param('xi', shape=Shapes.TxN) + ]) + + [S, I, R] = symbols.compartment_symbols + [β, γ, ξ] = symbols.attribute_symbols + + # N is NOT protected by Max(1, ...) here + N = S + I + R + + return create_model( + symbols=symbols, + transitions=[ + edge(S, I, rate=β * S * I / N), + edge(I, R, rate=γ * I), + edge(R, S, rate=ξ * R) + ]) + + my_geo = StaticGeo( + spec=StaticGeoSpec( + attributes=[ + AttribDef('label', dtype=str, shape=Shapes.N), + AttribDef('population', dtype=str, shape=Shapes.N), + ], + time_period=NO_DURATION, + ), + values={ + 'label': np.array(['a', 'b', 'c']), + 'population': np.array([0, 10, 20], dtype=np.int64), + }, + ) + sim = StandardSimulation( + geo=my_geo, + ipm=load_ipm(), + mm=mm_library['no'](), + params={ + 'phi': 40.0, + 'beta': 0.4, + 'gamma': 1 / 5, + 'xi': 1 / 100, + }, + time_frame=TimeFrame.of("2015-01-01", 150), + initializer=partial(single_location, location=1, seed_size=5), + rng=default_rng(1) + ) + + with self.assertRaises(IpmSimNaNException) as e: + sim.run() + + err_msg = str(e.exception) + + self.assertIn("NaN (not a number) rate detected", err_msg) + self.assertIn("Showing current Node : Timestep", err_msg) + self.assertIn("S: 0", err_msg) + self.assertIn("I: 0", err_msg) + self.assertIn("R: 0", err_msg) + self.assertIn("beta: 0.4", err_msg) + self.assertIn("gamma: 0.2", err_msg) + self.assertIn("xi: 0.01", err_msg) + self.assertIn("S->I: I*S*beta/(I + R + S)", err_msg) + + def test_negative_probs_error(self): + """Test for handling negative probability error""" + sim = StandardSimulation( + geo=geo_library['pei'](), + ipm=ipm_library['sirh'](), + mm=mm_library['no'](), + params={ + 'beta': 0.4, + 'gamma': 1 / 5, + 'xi': 1 / 100, + 'hospitalization_prob': -1 / 5, + 'hospitalization_duration': 15 + }, + time_frame=TimeFrame.of("2015-01-01", 150), + initializer=partial(single_location, location=1, seed_size=5), + rng=default_rng(1) + ) + + with self.assertRaises(IpmSimInvalidProbsException) as e: + sim.run() + + err_msg = str(e.exception) + + self.assertIn("Invalid probabilities for fork definition detected.", err_msg) + self.assertIn("Showing current Node : Timestep", err_msg) + self.assertIn("S: ", err_msg) + self.assertIn("I: ", err_msg) + self.assertIn("R: ", err_msg) + self.assertIn("beta: 0.4", err_msg) + self.assertIn("gamma: 0.2", err_msg) + self.assertIn("xi: 0.01", err_msg) + self.assertIn("hospitalization_prob: -0.2", err_msg) + self.assertIn("hospitalization_duration: 15", err_msg) + + self.assertIn("I->(H, R): I*gamma", err_msg) + self.assertIn( + "Probabilities: hospitalization_prob, 1 - hospitalization_prob", err_msg) + + def test_mm_clause_error(self): + """Test for handling invalid movement model clause application""" + + sim = StandardSimulation( + geo=geo_library['pei'](), + ipm=ipm_library['pei'](), + mm=mm_library['pei'](), + params={ + 'infection_duration': 40.0, + 'immunity_duration': 0.4, + 'humidity': 20.2, + 'move_control': 0.4, + 'theta': -5 + }, + time_frame=TimeFrame.of("2015-01-01", 150), + initializer=partial(single_location, location=1, seed_size=5), + rng=default_rng(1) + ) + + with self.assertRaises(MmSimException) as e: + sim.run() + + err_msg = str(e.exception) + + self.assertIn( + "Error from applying clause 'dispersers': see exception trace", err_msg)