Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Random Draws For SEIR Parameters Per Slot #429

Open
wants to merge 10 commits into
base: dev
Choose a base branch
from
132 changes: 132 additions & 0 deletions examples/tutorials/config_sample_2pop_modifiers_test_random.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
name: sample_2pop
setup_name: minimal
start_date: 2020-02-01
end_date: 2020-08-31
nslots: 1

subpop_setup:
geodata: model_input/geodata_sample_2pop.csv
mobility: model_input/mobility_sample_2pop.csv

initial_conditions:
method: SetInitialConditions
initial_conditions_file: model_input/ic_2pop.csv
allow_missing_subpops: TRUE
allow_missing_compartments: TRUE

compartments:
infection_stage: ["S", "E", "I", "R"]

seir:
integration:
method: rk4
dt: 1
parameters:
sigma:
value: 1 / 4
gamma:
value: 1 / 5
Ro:
value:
distribution: truncnorm
mean: 2.5
sd: 0.1
a: 1.1
b: 3
transitions:
- source: ["S"]
destination: ["E"]
rate: ["Ro * gamma"]
proportional_to: [["S"],["I"]]
proportion_exponent: ["1","1"]
- source: ["E"]
destination: ["I"]
rate: ["sigma"]
proportional_to: ["E"]
proportion_exponent: ["1"]
- source: ["I"]
destination: ["R"]
rate: ["gamma"]
proportional_to: ["I"]
proportion_exponent: ["1"]

seir_modifiers:
scenarios:
- Ro_all
modifiers:
Ro_lockdown:
method: SinglePeriodModifier
parameter: Ro
period_start_date: 2020-03-15
period_end_date: 2020-05-01
subpop: "all"
value: 0.4
Ro_relax:
method: SinglePeriodModifier
parameter: Ro
period_start_date: 2020-05-01
period_end_date: 2020-08-31
subpop: "all"
value:
distribution: truncnorm
mean: 0.8
sd: 0.1
a: 0.5
b: 1
Ro_all:
method: StackedModifier
modifiers: ["Ro_lockdown","Ro_relax"]


outcomes:
method: delayframe
outcomes:
incidCase: #incidence of detected cases
source:
incidence:
infection_stage: "I"
probability:
value: 0.5
delay:
value: 5
incidHosp: #incidence of hospitalizations
source:
incidence:
infection_stage: "I"
probability:
value: 0.05
delay:
value: 7
duration:
value: 10
name: currHosp # will track number of current hospitalizations (ie prevalence)
incidDeath: #incidence of deaths
source: incidHosp
probability:
value: 0.2
delay:
value:
distribution: truncnorm
mean: 14
sd: 2
a: 7
b: 21

outcome_modifiers:
scenarios:
- test_limits
modifiers:
# assume that due to limitations in testing, initially the case detection probability was lower
test_limits:
method: SinglePeriodModifier
parameter: incidCase::probability
subpop: "all"
period_start_date: 2020-02-01
period_end_date: 2020-06-01
value:
distribution: truncnorm
mean: 0.5
sd: 0.1
a: 0.1
b: 0.9

Original file line number Diff line number Diff line change
@@ -0,0 +1,134 @@
name: sample_2pop
setup_name: minimal
start_date: 2020-02-01
end_date: 2020-08-31
nslots: 1

subpop_setup:
geodata: model_input/geodata_sample_2pop.csv
mobility: model_input/mobility_sample_2pop.csv

initial_conditions:
method: SetInitialConditions
initial_conditions_file: model_input/ic_2pop.csv
allow_missing_subpops: TRUE
allow_missing_compartments: TRUE

compartments:
infection_stage: ["S", "E", "I", "R"]

seir:
integration:
method: rk4
dt: 1
parameters:
sigma:
value: 1 / 4
gamma:
value: 1 / 5
Ro:
value:
distribution: truncnorm
mean: 2.5
sd: 0.1
a: 1.1
b: 3
transitions:
- source: ["S"]
destination: ["E"]
rate: ["Ro * gamma"]
proportional_to: [["S"],["I"]]
proportion_exponent: ["1","1"]
- source: ["E"]
destination: ["I"]
rate: ["sigma"]
proportional_to: ["E"]
proportion_exponent: ["1"]
- source: ["I"]
destination: ["R"]
rate: ["gamma"]
proportional_to: ["I"]
proportion_exponent: ["1"]

seir_modifiers:
scenarios:
- Ro_all
modifiers:
Ro_lockdown:
method: SinglePeriodModifier
parameter: Ro
period_start_date: 2020-03-15
period_end_date: 2020-05-01
subpop: "all"
value: 0.4
Ro_relax:
method: SinglePeriodModifier
parameter: Ro
period_start_date: 2020-05-01
period_end_date: 2020-08-31
subpop: "all"
subpop_groups: "all"
value:
distribution: truncnorm
mean: 0.8
sd: 0.1
a: 0.5
b: 1
Ro_all:
method: StackedModifier
modifiers: ["Ro_lockdown","Ro_relax"]


outcomes:
method: delayframe
outcomes:
incidCase: #incidence of detected cases
source:
incidence:
infection_stage: "I"
probability:
value: 0.5
delay:
value: 5
incidHosp: #incidence of hospitalizations
source:
incidence:
infection_stage: "I"
probability:
value: 0.05
delay:
value: 7
duration:
value: 10
name: currHosp # will track number of current hospitalizations (ie prevalence)
incidDeath: #incidence of deaths
source: incidHosp
probability:
value: 0.2
delay:
value:
distribution: truncnorm
mean: 14
sd: 2
a: 7
b: 21

outcome_modifiers:
scenarios:
- test_limits
modifiers:
# assume that due to limitations in testing, initially the case detection probability was lower
test_limits:
method: SinglePeriodModifier
parameter: incidCase::probability
subpop: "all"
subpop_groups: "all"
period_start_date: 2020-02-01
period_end_date: 2020-06-01
value:
distribution: truncnorm
mean: 0.5
sd: 0.1
a: 0.1
b: 0.9

2 changes: 2 additions & 0 deletions flepimop/gempyor_pkg/src/gempyor/outcomes.py
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,8 @@ def onerun_delayframe_outcomes(
load_ID: bool = False,
sim_id2load: int = None,
):
np.random.seed(seed=sim_id2write)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This might fix a seed in a way that is not always desirable: with classical inference, every ith iteration will have the same seed across slots. I would rely on /dev/random for now to be sure things are working correctly and keep the reproducibility of runs as a follow-up issue.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

that suggests a fix is needed in how classical inference works - any inference method should be managing set seeding in a sensible way.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I fully agree. I would suggest a plan as follows:

  1. Make thing works correctly as people expect them (so doing the real random thing that the emcee is doing also here)
  2. Streamline the pathway to call gempyor (Descriptions of processes inside gempyor #463 (comment) shows 3 ways gempyor is called, but we should have one)
  3. Build something with seed so that the run with the above pathway is reproducible.

If we do 3. before 2. there might be places that have unexpected behavior which is IMO more critical than being reproducible (even though that's also very important)


with Timer("buildOutcome.structure"):
parameters = read_parameters_from_config(modinf)

Expand Down
12 changes: 12 additions & 0 deletions flepimop/gempyor_pkg/src/gempyor/parameters.py
Original file line number Diff line number Diff line change
Expand Up @@ -169,6 +169,18 @@ def __init__(
logging.debug(f"Index in arrays are: {self.pnames2pindex}")
logging.debug(f"NPI overlap operation is {self.stacked_modifier_method} ")

def reinitialize_distributions(self) -> None:
"""
Reinitialize all random distributions.

This method reinitializes all random distributions for the parameters contained
within this class. The random seed for each distribution is captured on
initialization so this method will change those seeds.
"""
for pn in self.pnames:
if "dist" in self.pdata[pn]:
self.pdata[pn]["dist"] = self.pconfig[pn]["value"].as_random_distribution()

def get_pnames2pindex(self) -> dict:
"""
Read the `pnames2pindex` attribute.
Expand Down
5 changes: 4 additions & 1 deletion flepimop/gempyor_pkg/src/gempyor/seir.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
import itertools
import logging
import random
import time

import numpy as np
Expand Down Expand Up @@ -257,7 +258,9 @@ def onerun_SEIR(
sim_id2load: int = None,
config=None,
):
np.random.seed()
np.random.seed(seed=sim_id2write)
modinf.parameters.reinitialize_distributions()

npi = None
if modinf.npi_config_seir:
npi = build_npi_SEIR(
Expand Down
Loading
Loading