Adding reliability analysis for ABM

Author: JBris

presentations/te_punaha_matatini/2025_10_calisim_workshop.pptx

@@ -0,0 +1,33 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "71fa9406-b967-4303-baec-9e4eaac19622",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.10.17"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

workshop/4_mesa_agent_based_models/1_geosir_network_reliability.ipynb

@@ -0,0 +1,33 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a14b59ed-5ef6-46fc-b61f-e5714c2d90d0",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.10.17"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

workshop/4_mesa_agent_based_models/2_geosir_network_surrogate.ipynb

@@ -7,6 +7,7 @@
 import mesa
 import mesa_geo as mg
 from shapely.geometry import Point
+import numpy as np
 
 class PersonAgent(mg.GeoAgent):
     """Person Agent."""
@@ -239,6 +240,7 @@ def __init__(
         self.space = mg.GeoSpace(warn_crs_conversion=False)
         self.networks = []
         self.states = []
+        self.r_0_ts = []
         self.network_grid = None
         network_grid_types=[
             "person",
@@ -263,6 +265,9 @@ def __init__(
                 "susceptible": get_susceptible_count,
                 "recovered": get_recovered_count,
                 "dead": get_dead_count,
+                "contact_r_e": get_contact_r_e,
+                "contact_r_0t": get_contact_r_0t,
+                "proxy_contact_r0": get_proxy_contact_r0
             }
         )
 
@@ -359,9 +364,49 @@ def run(self, n_steps) -> None:
         for _ in range(n_steps):
             self.step()
 
+    def calculate_r_e(self) -> float:
+        induced_infections = [
+            a.induced_infections_at_t for a in self.agents_by_type[PersonAgent]
+            if a.infected_others_at_t == True
+        ]
+        n_infections = np.array(induced_infections)
+        if len(n_infections) == 0:
+            effective_r = 0
+        else:
+            effective_r = np.average(n_infections)
+
+        return effective_r
+
+    def calculate_r_0t(self, r_e: int) -> float:
+        S = get_susceptible_count(self)
+        if S == 0:
+            S = 1e-8
+
+        N = self.pop_size
+        initial_outbreak_size = 0
+        r_0_t = r_e * (N - initial_outbreak_size - 0) / S
+
+        self.r_0_ts.append(r_0_t)
+        return r_0_t
+
+    def calculate_proxy_contact_r0(self) -> float:
+        return np.array(self.r_0_ts).mean()
+
+    def resistant_susceptible_ratio(self):
+        S = self.counts["SUSCEPTIBLE"]
+        if S == 0:
+            S = 1e-8
+
+        return self.counts["RESISTANT"] / S
+
     def get_state(self) -> dict:
         state = deepcopy(self.counts)
         state["t"] = self.steps
+
+        state["contact_r_e"] = self.calculate_r_e()
+        state["contact_r_0t"] = self.calculate_r_0t(state["contact_r_e"])
+        state["proxy_contact_r0"] = self.calculate_proxy_contact_r0()
+
         return state
 
     def get_states(self):
@@ -419,7 +464,6 @@ def get_edge_df(self):
         df = pd.DataFrame(rows)
         return df
 
-# Functions needed for datacollector
 def get_infected_count(model):
     return model.counts["infected"]
 
@@ -431,17 +475,15 @@ def get_susceptible_count(model):
 def get_recovered_count(model):
     return model.counts["recovered"]
 
-
 def get_dead_count(model):
     return model.counts["dead"]
 
+def get_contact_r_e(model):
+      return model.calculate_r_e()
 
-model = GeoSir()
-
-for t in range(10):
-    model.run(1)
+def get_contact_r_0t(model):
+    r_e = model.calculate_r_e()
+    return model.calculate_r_0t(r_e)
 
-node_df = model.get_node_df()
-edge_df = model.get_edge_df()
-print(node_df)
-print(edge_df)
+def get_proxy_contact_r0(model):
+    return model.calculate_proxy_contact_r0()