Replaced D reporting channel with W. Added add_reporting_property fun…

…ction so per-tick arrays are contiguous. Added importation function upon elimination. Reduced RI fraction.

nnmm/final_reports.py

@@ -9,8 +9,8 @@ def save_seird( model ):
  'S': model.nodes.S,
  'E': model.nodes.E,
  'I': model.nodes.I,
- 'R': model.nodes.R,
- 'D': model.nodes.D
+ 'W': model.nodes.W,
+ 'R': model.nodes.R
  }
 
  # Save the arrays to an HDF5 file

nnmm/measles_mod.py

@@ -29,10 +29,10 @@ class Model:
  "ticks": int(365*10),
  "cbr": 40, # Nigeria 2015 according to (somewhat random internet source): https://fred.stlouisfed.org/series/SPDYNCBRTINNGA
  "output": Path.cwd() / "outputs",
- "eula_age": 5
+ "eula_age": 5,
+ "prevalence": 0.025 # 2.5% prevalence
 })
 # parameter?
-prevalence = 0.025 # 2.5% prevalence
 
 measles_params = PropertySet({
  "exp_mean": np.float32(7.0),
@@ -42,7 +42,7 @@ class Model:
  "r_naught": np.float32(14.0),
  "seasonality_factor": np.float32(0.125),
  "seasonality_phase": np.float32(182),
- "ri_coverage": np.float32(0.75),
+ "ri_coverage": np.float32(0.25),
 })
 
 network_params = PropertySet({
@@ -120,7 +120,6 @@ def save_pops_in_nodes( model, nn_nodes, initial_populations):
 model.nodes.add_scalar_property("ri_coverages", dtype=np.float32)
 model.nodes.ri_coverages = np.random.rand(len(nn_nodes))
 # ri coverages and init prev seem to be the same "kind of thing"?
-model.nodes.initial_infections = np.uint32(np.round(np.random.poisson(prevalence*initial_populations)))
 
 
 # ## Population per Tick
@@ -170,10 +169,11 @@ def propagate_population(model, tick):
 # In[24]:
 
 
+import pdb
 # consider `step_functions` rather than `phases` for the following
 model.phases = [
- propagate_population,
  ages.update_ages, # type: ignore
+ propagate_population,
  fertility.do_births, # type: ignore
  mortality.do_non_disease_deaths, # type: ignore
  intrahost.do_infection_update, # type: ignore
@@ -188,7 +188,6 @@ def propagate_population(model, tick):
 # ## Running the Simulation
 # 
 # We iterate over the specified number of ticks, keeping track, in `metrics`, of the time spent in each phase at each tick.
-
 model.metrics = []
 for tick in tqdm(range(model.params.ticks)):
  metrics = [tick]
@@ -198,6 +197,10 @@ def propagate_population(model, tick):
  tfinish = datetime.now(tz=None) # noqa: DTZ005
  delta = tfinish - tstart
  metrics.append(delta.seconds * 1_000_000 + delta.microseconds) # delta is a timedelta object, let's convert it to microseconds
+ # age_and_report should go last so we can use tick but for some reason it won't fire when put later.
+ if tick>0 and sum(model.nodes.I[tick-1]) == 0 and sum(model.nodes.E[tick-1]) == 0:
+ print( f"WARNING: Elimination achieved at t={tick}. Reimporting." )
+ init_prev.init(model) 
  model.metrics.append(metrics)
 
 

nnmm/mods/ages.py

@@ -28,23 +28,23 @@
  np.ctypeslib.ndpointer(dtype=np.uint8, flags='C_CONTIGUOUS'), # infectious_timer
  np.ctypeslib.ndpointer(dtype=np.uint8, flags='C_CONTIGUOUS'), # incubation_timer
  np.ctypeslib.ndpointer(dtype=np.uint8, flags='C_CONTIGUOUS'), # immunity
- np.ctypeslib.ndpointer(dtype=np.int32, flags='C_CONTIGUOUS'), # age
+ np.ctypeslib.ndpointer(dtype=np.uint8, flags='C_CONTIGUOUS'),  # susceptibility_timer
  np.ctypeslib.ndpointer(dtype=np.int32, flags='C_CONTIGUOUS'), # expected_lifespan
  np.ctypeslib.ndpointer(dtype=np.uint32, flags='C_CONTIGUOUS'), # infectious_count
  np.ctypeslib.ndpointer(dtype=np.uint32, flags='C_CONTIGUOUS'), # incubating_count
  np.ctypeslib.ndpointer(dtype=np.uint32, flags='C_CONTIGUOUS'), # susceptible_count
- np.ctypeslib.ndpointer(dtype=np.uint32, flags='C_CONTIGUOUS'), # recovered_count
- np.ctypeslib.ndpointer(dtype=np.uint32, flags='C_CONTIGUOUS')  # dead_count
+ np.ctypeslib.ndpointer(dtype=np.uint32, flags='C_CONTIGUOUS'), # waning_count 
+ np.ctypeslib.ndpointer(dtype=np.uint32, flags='C_CONTIGUOUS'), # recovered_count 
  ]
 except Exception as ex:
  print( f"Failed to load libages.so." )
 
 def init( model ):
- model.nodes.add_vector_property("S", model.params.ticks, dtype=np.uint32)
- model.nodes.add_vector_property("E", model.params.ticks, dtype=np.uint32)
- model.nodes.add_vector_property("I", model.params.ticks, dtype=np.uint32)
- model.nodes.add_vector_property("R", model.params.ticks, dtype=np.uint32)
- model.nodes.add_vector_property("D", model.params.ticks, dtype=np.uint32)
+ model.nodes.add_report_property("S", model.params.ticks, dtype=np.uint32)
+ model.nodes.add_report_property("E", model.params.ticks, dtype=np.uint32)
+ model.nodes.add_report_property("I", model.params.ticks, dtype=np.uint32)
+ model.nodes.add_report_property("W", model.params.ticks, dtype=np.uint32) 
+ model.nodes.add_report_property("R", model.params.ticks, dtype=np.uint32) 
 
 
 def update_ages( model, tick ):
@@ -54,11 +54,6 @@ def update_ages( model, tick ):
  model.population.age,
  )
  """
- S = np.zeros( len( model.nodes.S), dtype=np.uint32 )
- E = np.zeros( len( model.nodes.S), dtype=np.uint32 )
- I = np.zeros( len( model.nodes.S), dtype=np.uint32 )
- R = np.zeros( len( model.nodes.S), dtype=np.uint32 )
- D = np.zeros( len( model.nodes.S), dtype=np.uint32 )
  lib.update_ages_and_report(
  ctypes.c_int64(model.population.count),
  len(model.nodes.nn_nodes),
@@ -67,17 +62,13 @@ def update_ages( model, tick ):
  model.population.itimer,
  model.population.etimer,
  model.population.susceptibility,
- model.population.age,
+ model.population.susceptibility_timer,
  model.population.dod,
- S,
- E,
- I,
- R,
- D
+ model.nodes.S[tick],
+ model.nodes.E[tick],
+ model.nodes.I[tick],
+ model.nodes.W[tick],
+ model.nodes.R[tick],
  )
- model.nodes.S[:,tick]=S
- model.nodes.E[:,tick]=E
- model.nodes.I[:,tick]=I
- model.nodes.R[:,tick]=R
- model.nodes.D[:,tick]=D
+
 

nnmm/mods/init_prev.py

@@ -16,22 +16,29 @@
 
 # initial_infections = np.random.randint(0, 11, model.nodes.count, dtype=np.uint32)
 
-@nb.njit((nb.uint32, nb.uint32[:], nb.uint16[:], nb.uint8[:], nb.float32, nb.float32), parallel=True)
-def initialize_infections(count, infections, nodeid, itimer, inf_mean, inf_std):
+@nb.njit((nb.uint32, nb.uint32[:], nb.uint16[:], nb.uint8[:], nb.uint8[:], nb.uint8[:], nb.float32, nb.float32), parallel=True)
+def initialize_infections(count, infections, nodeid, itimer, etimer, susceptibility, inf_mean, inf_std):
 
  for i in nb.prange(count):
- if infections[nodeid[i]] > 0:
- infections[nodeid[i]] -= 1
- itimer[i] = np.maximum(np.uint8(1), np.uint8(np.round(np.random.normal(inf_mean, inf_std)))) # must be at least 1 day
+ # Only seed susceptible uninfecteds
+ if itimer[i] == 0 and etimer[i] == 0 and susceptibility[i] == 1:
+ if infections[nodeid[i]] > 0:
+ infections[nodeid[i]] -= 1
+ itimer[i] = np.maximum(np.uint8(1), np.uint8(np.round(np.random.normal(inf_mean, inf_std)))) # must be at least 1 day
 
  return
 
 def init( model ):
+ init_pop = model.nodes.population[:,0] # should probably be 'now', not t=0
+ model.nodes.initial_infections = np.uint32(np.round(np.random.poisson(model.params.prevalence*init_pop)))
+ #print( f"Seeding infections: {model.nodes.initial_infections}" )
  return initialize_infections(
  np.uint32(model.population.count),
  model.nodes.initial_infections,
  model.population.nodeid,
  model.population.itimer,
+ model.population.etimer,
+ model.population.susceptibility,
  model.params.inf_mean,
  model.params.inf_std
  )

nnmm/mods/update_ages.c

@@ -73,18 +73,18 @@ void update_ages_simd(unsigned long int stop_idx, float *ages) {
 void update_ages_and_report(
  unsigned long int count,
  int num_nodes,
- int *ages,
- u_int16_t *node,
- unsigned char *infectious_timer,
- unsigned char *incubation_timer,
- bool *susceptibility,
- int *age,
- int *dod,
+ int32_t *age,
+ uint16_t *node,
+ unsigned char *infectious_timer, // max 255
+ unsigned char *incubation_timer, // max 255
+ bool *susceptibility, // yes or no
+ unsigned char *susceptibility_timer, // max 255
+ int *dod, // sim day
  uint32_t *susceptible_count,
  uint32_t *incubating_count,
  uint32_t *infectious_count,
- uint32_t *recovered_count,
- uint32_t *dead_count
+ uint32_t *waning_count,
+ uint32_t *recovered_count 
 ) {
  //printf( "%s: count=%ld, num_nodes=%d", __FUNCTION__, count, num_nodes );
  #pragma omp parallel
@@ -94,51 +94,57 @@ void update_ages_and_report(
  int *local_incubating_count = (int*) calloc(num_nodes, sizeof(int));
  int *local_recovered_count = (int*) calloc(num_nodes, sizeof(int));
  int *local_susceptible_count = (int*) calloc(num_nodes, sizeof(int));
- int *local_dead_count = (int*) calloc(num_nodes, sizeof(int));
+ int *local_waning_count = (int*) calloc(num_nodes, sizeof(int));
 
  #pragma omp for
  for (size_t i = 0; i <= count; i++) {
  // Update ages
- if (ages[i] >= 0) {
- ages[i]++;
+ if (age[i] >= 0) {
+ age[i]++;
  }
 
- // Collect report
- if (node[i] >= 0 ) {
+ // Collect report 
+ if (dod[i]>0) {
  int node_id = node[i];
- //printf( "node_id=%d.\n", node_id );
- if (dod[i]==0) {
- local_dead_count[node_id]++;
- } else if (incubation_timer[i] > 0) {
+ //printf( "Found live person at node %d: etimer=%d, itimer=%d, sus=%d.\n", node_id, incubation_timer[i], infectious_timer[i], susceptibility[i] );
+ if (incubation_timer[i] > 0) {
+ //printf( "Found E in node %d.\n", node_id );
  local_incubating_count[node_id]++;
- } else if (infectious_timer[i]) {
+ } else if (infectious_timer[i] > 0) {
+ //printf( "Found I in node %d.\n", node_id );
  local_infectious_count[node_id]++;
  } else if (susceptibility[i]==0) {
- local_recovered_count[node_id]++;
+ //printf( "Found R in node %d.\n", node_id );
+ if (susceptibility_timer[i]>0) {
+ local_waning_count[node_id]++;
+ } else {
+ local_recovered_count[node_id]++;
+ }
  } else {
+ //printf( "Found S in node %d.\n", node_id );
  local_susceptible_count[node_id]++;
  }
  }
  }
 
  // Combine local counts into global counts
- #pragma omp critical
+#pragma omp critical
  {
  for (int j = 0; j < num_nodes; ++j) {
- infectious_count[j] += local_infectious_count[j];
+ susceptible_count[j] += local_susceptible_count[j]; 
  incubating_count[j] += local_incubating_count[j];
+ infectious_count[j] += local_infectious_count[j];
+ waning_count[j] += local_waning_count[j];
  recovered_count[j] += local_recovered_count[j];
- susceptible_count[j] += local_susceptible_count[j];
- dead_count[j] += local_dead_count[j];
  }
  }
 
  // Free local buffers
- free(local_infectious_count);
+ free(local_susceptible_count);
  free(local_incubating_count);
+ free(local_infectious_count);
+ free(local_waning_count);
  free(local_recovered_count);
- free(local_susceptible_count);
- free(local_dead_count);
  }
 }
 

src/idmlaser/numpynumba/population.py

@@ -83,6 +83,11 @@ def add_vector_property(self, name, length: int, dtype=np.uint32, default=0) ->
  setattr(self, name, np.full((self._capacity, length), default, dtype=dtype))
  return
 
+ def add_report_property(self, name, length: int, dtype=np.uint32, default=0) -> None:
+ """Add a vector property to the class"""
+ # initialize the property to a NumPy array with of size self._count, dtype, and default value
+ setattr(self, name, np.full((length, self._capacity), default, dtype=dtype))
+ return
  # Add scalar properties to model.population
  def add_properties_from_schema( self, schema ):
  for name, dtype in schema.items():