Bird and mosquito movement stages

epymorph/data/mm/local_birds.movement

@@ -0,0 +1,19 @@
+
+[move-steps: per-day=2; duration=[1/3,2/3]]
+
+[predef:  function=
+def birds_movement():
+   shape = 3
+   scale = 3
+   centroid = geo['centroid']
+   distance = pairwise_haversine(centroid['longitude'], centroid['latitude'])
+   weibull_distribution = weibull_distribution_prob(distance, shape, scale)
+   result = row_normalize(weibull_distribution)
+   return {'result' : result}
+]
+[mtype:days=all; leave=1; duration=0d; return = 2; function=
+def bird_commuters_source(t,src):
+    n_commuters = np.floor(geo['bird_population'][src]*0.1).astype(SimDType)
+    return np.multinomial(n_commuters,predef['result'][src])
+
+]

epymorph/data/mm/long_local_birds_model.movement

@@ -0,0 +1,29 @@
+[move-steps: per-day=2; duration=[1/3,2/3]]
+
+[predef:  function=
+def long_range_birds_movement():
+   centroid = geo['centroid']
+   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, params['shape'], params['scale'])
+   result_weibull = row_normalize(weibull_distribution)
+   return {'result_powerlaw' :result_powerlaw,
+   'result_weibull' : result_weibull
+   }
+]
+
+
+[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'][src]*0.1).astype(SimDType)
+    return np.multinomial(n_commuters,predef['result_powerlaw'][src])  
+
+]
+
+[mtype:days=all; leave=1; duration=0d; return = 2; function=
+def bird_commuters_source(t,src):
+    n_commuters = np.floor(geo['bird_population'][src]*0.1).astype(SimDType)
+    return np.multinomial(n_commuters,predef['result_weibull'][src])
+
+]

epymorph/data/mm/long_range_dispersal_birds.movement

@@ -0,0 +1,17 @@
+[move-steps : per-day=2; duration =[1/3,2/3]]
+
+
+[predef:  function=
+def long_range_birds_movement():
+   centroid = geo['centroid']
+   distance = pairwise_haversine(centroid['longitude'], centroid['latitude'])
+   powerlaw_distribution = powerlaw_distribution_probability(distance, params['alpha'])
+   result = row_normalize(1/powerlaw_distribution)
+   return {'result' : result}
+]
+
+[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'][src]*0.1).astype(SimDType)
+    return np.multinomial(n_commuters,predef['result'][src])  
+]

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])  
+]

epymorph/data/mm/mosquito_mobility.movement

@@ -0,0 +1,16 @@
+[move-steps: per-day=2; duration=[1/3,2/3]]
+
+[predef:  function=
+def mosquito_movement():
+   centroid = geo['centroid']
+   distance = pairwise_haversine(centroid['longitude'], centroid['latitude'])
+   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)
+    movers = np.poisson(n_commuters*0.4)
+    return np.multinomial(movers,predef['result'][src])
+]

epymorph/data/mm/my_experiment.movement

@@ -0,0 +1,18 @@
+[move-steps:per-day=2; duration=[2/3,1/3]]
+
+[mtype: days = [M,W,F,Su];leave=1;duration=0d;return=2; function=
+def icecube_clause_1(t,src):
+   commuters_clause_1 = np.zeros_like(geo['population'])
+   if (src+1) < commuters_clause_1.size:
+     commuters_clause_1[src+1] = 5000
+   return commuters_clause_1
+]  
+
+[mtype:days=[T,Th,Sa]; leave=1;duration=0d;return=2; function=
+def icecube_clause_2(t,dst):
+   commuters_clause_2 = np.zeros_like(geo['population'])
+   if (dst) >= 1:
+     commuters_clause_2[dst-1] = 50
+   return commuters_clause_2
+
+]

epymorph/initializer.py

@@ -270,6 +270,18 @@ def bottom_locations(ctx: InitContext, attribute: str, num_locations: int, seed_
     return indexed_locations(ctx, selection, seed_size)
 
 
+def bird_movement_initializer(ctx: InitContext) -> NDArray[SimDType]:
+    result = ctx.geo['bird_population'].reshape(
+        (ctx.dim.nodes, ctx.dim.compartments))  # type: ignore
+    return result.astype(SimDType)
+
+
+def mosquito_movement_initializer(ctx: InitContext) -> NDArray[SimDType]:
+    result = ctx.geo['bird_population'].reshape(
+        (ctx.dim.nodes, ctx.dim.compartments))  # type: ignore
+    return result.astype(SimDType)
+
+
 DEFAULT_INITIALIZER = single_location
 
 initializer_library: dict[str, Initializer] = {
@@ -281,6 +293,8 @@ def bottom_locations(ctx: InitContext, attribute: str, num_locations: int, seed_
     'random_locations': random_locations,
     'top_locations': top_locations,
     'bottom_locations': bottom_locations,
+    'bird_movement_intializer': bird_movement_initializer,
+    'mosquito_movement_initializer': mosquito_movement_initializer
 }
 """A library for the built-in initializer functions."""
 

epymorph/simulation.py

@@ -13,8 +13,11 @@
 from numpy.random import SeedSequence
 
 from epymorph.code import ImmutableNamespace, base_namespace
-from epymorph.util import (Event, pairwise_haversine, progress, row_normalize,
-                           subscriptions)
+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)
 
 SimDType = np.int64
 """
@@ -218,6 +221,10 @@ 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,
         'row_normalize': row_normalize,
         # numpy namespace
         'np': ImmutableNamespace({

epymorph/test/util_test.py

@@ -83,3 +83,12 @@ def test_check_ndarray_03(self):
             util.check_ndarray(arr, dtype=np.str_)
         with self.assertRaises(util.NumpyTypeError):
             util.check_ndarray(arr, dimensions=3)
+
+    def test_weibull_distribution(self):
+        act = util.weibull_distribution_prob(
+            np.linspace(0, 2.5, 10), shape=1.5, scale=1)
+        # the expected values were generated by scipy:
+        #   weibull_min.pdf(np.linspace(0, 2.5, 10), c=1.5)
+        exp = np.array([0., 0.68290224, 0.73895749, 0.63991403, 0.49013786,
+                        0.34400327, 0.22519705, 0.13897726, 0.08143303, 0.0455367])
+        np.testing.assert_allclose(act, exp)

epymorph/util.py

@@ -1,4 +1,6 @@
 """epymorph general utility functions and classes."""
+from __future__ import annotations
+
 from contextlib import contextmanager
 from typing import (Any, Callable, Generator, Generic, Iterable, Literal,
                     OrderedDict, TypeVar)
@@ -166,6 +168,45 @@ def row_normalize(arr: NDArray[N], row_sums: NDArray[N] | None = None, dtype: DT
 RADIUS_MI = 3959.87433  # radius of earth in mi
 
 
+def weibull_distribution_prob(distance: NDArray[N], shape: float, scale: float) -> NDArray[np.float64]:
+    """
+    Calculate the Weibull distribution probability density function for each value in `x`.
+    `shape` and `scale` are parameters which alter the shape of the distribution,
+    and each must be greater than zero.
+    https://en.wikipedia.org/wiki/Weibull_distribution
+    """
+    if shape <= 0 or scale <= 0:
+        raise ValueError("`shape` and `scale` must be greater than zero.")
+    x = distance.astype(np.float64)
+    return ((shape / scale)
+            * ((x / scale) ** (shape - 1))
+            * (np.exp(-((x / scale) ** shape))))
+
+
+def powerlaw_distribution_probability(distance: NDArray[N], alpha: float) -> NDArray[np.float64]:
+    result = np.zeros_like(distance, dtype=np.float64)
+    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.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