SEIR Movement model

NAU-CCL · Feb 2, 2024 · e0b16c9 · e0b16c9
1 parent 4b5a829
commit e0b16c9
Show file tree

Hide file tree

Showing 11 changed files with 371 additions and 21 deletions.
diff --git a/doc/devlog/example_local_birds.ipynb b/doc/devlog/example_local_birds.ipynb
@@ -37,15 +37,32 @@
    "cell_type": "code",
    "execution_count": 3,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "ename": "KeyboardInterrupt",
+     "evalue": "",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
+      "Cell \u001b[1;32mIn[3], line 13\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,\n\u001b[0;32m      3\u001b[0m     ipm\u001b[38;5;241m=\u001b[39mipm_library[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mno\u001b[39m\u001b[38;5;124m'\u001b[39m](),\n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m     10\u001b[0m     initializer\u001b[38;5;241m=\u001b[39mbird_movement_initializer,\n\u001b[0;32m     11\u001b[0m )\n\u001b[1;32m---> 13\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\\sp2665\\Downloads\\Sindhu\\Research EpyMoRPH\\Epymorph\\epymorph\\engine\\standard_sim.py:151\u001b[0m, in \u001b[0;36mStandardSimulation.run\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    149\u001b[0m     ctx \u001b[38;5;241m=\u001b[39m RumeContext\u001b[38;5;241m.\u001b[39mfrom_config(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_config)\n\u001b[0;32m    150\u001b[0m     ipm_exec \u001b[38;5;241m=\u001b[39m StandardIpmExecutor(ctx)\n\u001b[1;32m--> 151\u001b[0m     movement_exec \u001b[38;5;241m=\u001b[39m \u001b[43mStandardMovementExecutor\u001b[49m\u001b[43m(\u001b[49m\u001b[43mctx\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    153\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m error_gate(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124minitializing the simulation\u001b[39m\u001b[38;5;124m\"\u001b[39m, InitException):\n\u001b[0;32m    154\u001b[0m     ini \u001b[38;5;241m=\u001b[39m ctx\u001b[38;5;241m.\u001b[39minitialize()\n",
+      "File \u001b[1;32mc:\\Users\\sp2665\\Downloads\\Sindhu\\Research EpyMoRPH\\Epymorph\\epymorph\\engine\\mm_exec.py:58\u001b[0m, in \u001b[0;36mStandardMovementExecutor.__init__\u001b[1;34m(self, ctx)\u001b[0m\n\u001b[0;32m     56\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_log \u001b[38;5;241m=\u001b[39m getLogger(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmovement\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m     57\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_clause_masks \u001b[38;5;241m=\u001b[39m {c: c\u001b[38;5;241m.\u001b[39mmask(ctx) \u001b[38;5;28;01mfor\u001b[39;00m c \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_model\u001b[38;5;241m.\u001b[39mclauses}\n\u001b[1;32m---> 58\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_check_predef\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[1;32mc:\\Users\\sp2665\\Downloads\\Sindhu\\Research EpyMoRPH\\Epymorph\\epymorph\\engine\\mm_exec.py:65\u001b[0m, in \u001b[0;36mStandardMovementExecutor._check_predef\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m     63\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m curr_hash \u001b[38;5;241m!=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_predef_hash:\n\u001b[0;32m     64\u001b[0m     \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m---> 65\u001b[0m         \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_predef \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_model\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpredef\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\n\u001b[0;32m     66\u001b[0m         \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_predef_hash \u001b[38;5;241m=\u001b[39m curr_hash\n\u001b[0;32m     67\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     68\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     69\u001b[0m         \u001b[38;5;66;03m# until we can properly validate the MM clauses.\u001b[39;00m\n",
+      "File \u001b[1;32m<string>:6\u001b[0m, in \u001b[0;36mbirds_movement\u001b[1;34m(ctx)\u001b[0m\n",
+      "File \u001b[1;32mc:\\Users\\sp2665\\Downloads\\Sindhu\\Research EpyMoRPH\\Epymorph\\epymorph\\util.py:174\u001b[0m, in \u001b[0;36mweibull_distribution_prob\u001b[1;34m(distance, shape, scale)\u001b[0m\n\u001b[0;32m    171\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mweibull_distribution_prob\u001b[39m(distance: NDArray[N], shape: \u001b[38;5;28mfloat\u001b[39m, scale: \u001b[38;5;28mfloat\u001b[39m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m NDArray[np\u001b[38;5;241m.\u001b[39mfloat64]:\n\u001b[0;32m    172\u001b[0m     result \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mzeros_like(distance, dtype\u001b[38;5;241m=\u001b[39mnp\u001b[38;5;241m.\u001b[39mfloat64)\n\u001b[0;32m    173\u001b[0m     result \u001b[38;5;241m=\u001b[39m ((shape \u001b[38;5;241m/\u001b[39m scale) \u001b[38;5;241m*\u001b[39m ((distance \u001b[38;5;241m/\u001b[39m scale) \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m (shape \u001b[38;5;241m-\u001b[39m \u001b[38;5;241m1\u001b[39m))\n\u001b[1;32m--> 174\u001b[0m               \u001b[38;5;241m*\u001b[39m (\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mexp\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdistance\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mscale\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mshape\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m))\n\u001b[0;32m    175\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m result\n",
+      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
+     ]
+    }
+   ],
    "source": [
     "sim = StandardSimulation(\n",
     "    geo=geo,\n",
     "    ipm=ipm_library['no'](),\n",
     "    mm=mm_library['local_birds'](),\n",
     "    params={\n",
     "        'shape': 3,\n",
-    "        'scale':3\n",
+    "        'scale':100\n",
     "    },\n",
     "    time_frame=TimeFrame.of(\"2023-01-01\", 10),\n",
     "    initializer=bird_movement_initializer,\n",

diff --git a/doc/devlog/example_long_birds.ipynb b/doc/devlog/example_long_birds.ipynb
diff --git a/doc/devlog/example_mosquito.ipynb b/doc/devlog/example_mosquito.ipynb
diff --git a/doc/devlog/probability_vs_distance.ipynb b/doc/devlog/probability_vs_distance.ipynb
diff --git a/doc/devlog/von_mises_distribution.ipynb b/doc/devlog/von_mises_distribution.ipynb
diff --git a/epymorph/data/mm/long_local_birds_model.movement b/epymorph/data/mm/long_local_birds_model.movement
@@ -6,7 +6,7 @@ def long_range_birds_movement():
    distance = pairwise_haversine(centroid['longitude'], centroid['latitude'])
    powerlaw_distribution = powerlaw_distribution_probability(distance, param['alpha'])
    result_powerlaw = row_normalize(1/powerlaw_distribution)
-   weibull_distribution = weibull_distribution_prob(distance, param['shape'], param['scale'])
+   weibull_distribution = weibull_distribution_prob(distance, params['shape'], params['scale'])
    result_weibull = row_normalize(weibull_distribution)
    return {'result_powerlaw' :result_powerlaw,
    'result_weibull' : result_weibull

diff --git a/epymorph/data/mm/long_range_dispersal_birds.movement b/epymorph/data/mm/long_range_dispersal_birds.movement
@@ -5,7 +5,7 @@
 def long_range_birds_movement():
    centroid = geo['centroid']
    distance = pairwise_haversine(centroid['longitude'], centroid['latitude'])
-   powerlaw_distribution = powerlaw_distribution_probability(distance, param['alpha'])
+   powerlaw_distribution = powerlaw_distribution_probability(distance, params['alpha'])
    result = row_normalize(1/powerlaw_distribution)
    return {'result' : result}
 ]

diff --git a/epymorph/data/mm/long_range_seasonal.movement b/epymorph/data/mm/long_range_seasonal.movement
@@ -0,0 +1,14 @@
+[move-steps : per-day=2; duration =[1/3,2/3]]
+
+
+[predef:  function=
+def long_range_birds_movement():
+   tau = np.array([np.pi/2, np.pi,3*np.pi/2,2*np.pi])
+   vonmises_distribution = long_range_von_mises_distribution(params['sij'], params['phiij'], params['deltaij'], tau)
+   return {'vonmises_distribution' :vonmises_distribution}
+]
+[mtype:days=all; leave=1; duration=0d; return = 2; function=
+def long_range_bird_commuters_source(t,src):
+    n_commuters = np.floor(geo['bird_population']*0.1).astype(SimDType)
+    return np.multinomial(n_commuters,predef['vonmises_distribution'][src])  
+]
diff --git a/epymorph/data/mm/mosquito_mobility.movement b/epymorph/data/mm/mosquito_mobility.movement
@@ -4,12 +4,13 @@
 def mosquito_movement():
    centroid = geo['centroid']
    distance = pairwise_haversine(centroid['longitude'], centroid['latitude'])
-   vector_movement_probability = mosquito_movement_probability(distance, param['max_distance'])
+   vector_movement_probability = mosquito_movement_probability(distance, params['max_distance'])
    result = row_normalize(vector_movement_probability)
    return {'result' : result}
 ]
 [mtype:days=all; leave=1; duration=0d; return = 2; function=
 def mosquito_commuters_source(t,src):
     n_commuters = np.floor(geo['bird_population'][src]*0.1).astype(SimDType)
-    return np.multinomial(n_commuters,predef['result'][src])
+    movers = np.poisson(n_commuters*0.4)
+    return np.multinomial(movers,predef['result'][src])
 ]
diff --git a/epymorph/simulation.py b/epymorph/simulation.py
@@ -13,8 +13,8 @@
 from numpy.random import SeedSequence
 
 from epymorph.code import ImmutableNamespace, base_namespace
-from epymorph.util import (Event, mosquito_movement_probability,
-                           pairwise_haversine,
+from epymorph.util import (Event, long_range_von_mises_distribution,
+                           mosquito_movement_probability, pairwise_haversine,
                            powerlaw_distribution_probability, progress,
                            row_normalize, subscriptions,
                            weibull_distribution_prob)
@@ -221,6 +221,7 @@ def epymorph_namespace() -> dict[str, Any]:
         'SimDType': SimDType,
         # our utility functions
         'pairwise_haversine': pairwise_haversine,
+        'long_range_von_mises_distribution': long_range_von_mises_distribution,
         'powerlaw_distribution_probability': powerlaw_distribution_probability,
         'weibull_distribution_prob': weibull_distribution_prob,
         'mosquito_movement_probability': mosquito_movement_probability,

diff --git a/epymorph/util.py b/epymorph/util.py
@@ -177,18 +177,28 @@ def weibull_distribution_prob(distance: NDArray[N], shape: float, scale: float)
 
 def powerlaw_distribution_probability(distance: NDArray[N], alpha: float) -> NDArray[np.float64]:
     result = np.zeros_like(distance, dtype=np.float64)
-    result = (1 / (distance ** alpha))  # type: ignore
+    min_distance = 22530 * 0.00062
+    normalize_distance = distance / min_distance
+    result = ((alpha - 1) / min_distance) * \
+        ((normalize_distance / min_distance)**-alpha)  # type: ignore
     return result  # type:ignore
 
 
 def mosquito_movement_probability(distance: NDArray[N], max_distance: float) -> NDArray[np.float64]:
     result = np.zeros_like(distance, dtype=np.float64)
     max_distance_mosquito = max_distance * 0.00062
     result = ((max_distance_mosquito) - (distance)) / (max_distance_mosquito)
-    result = np.clip(result, 0, 1)
+    result = np.maximum(result, 0)
     return result  # type:ignore
 
 
+def long_range_von_mises_distribution(sij: float, phiij: float, deltaij: float, tau: NDArray[np.float64]) -> NDArray[np.float64]:
+    exponent_term = np.exp(np.cos(tau - phiij) - 1)
+    tij = sij * (exponent_term / deltaij**2)
+    tij /= np.max(tij)
+    return tij
+
+
 def pairwise_haversine(longitudes: NDArray[np.float64], latitudes: NDArray[np.float64]) -> NDArray[np.float64]:
     """Compute the distances in miles between all pairs of coordinates."""
     # https://www.themathdoctors.org/distances-on-earth-2-the-haversine-formula