Descriptions of processes inside gempyor

[jcblemai]

Here lie a description of how data is used and shared inside gempyor.

• How EMCEE inference works Descriptions of processes inside gempyor #463 (comment)
• How classical inference works Descriptions of processes inside gempyor #463 (comment)
• How the parameters for the model integrations are computed: Descriptions of processes inside gempyor #463 (comment)
• The 3 ways (should be one) to run the full gempyor model: Descriptions of processes inside gempyor #463 (comment)

What remains to be written (ask and I'll):

• How outcomes work
• How modifiers are built

[jcblemai]

SEIR Parameters

Are build in several steps:

1. we find out which parameters are (file: spar -- config: seir: parameter:)
2. we compute the modifiers that applies to these parameters (file: snpi -- config: seir_modifiers)
3. we apply the modifiers to the parameters
4. we parse the formulas from the config to find out the unique parameters (e.g if you have beta and gamma from step 3, but only their products is used, it will get computed) (config only: seir: compartments)
5. we apply the formulas from step 4 to the parameters from step 3.

In the code, it is this workflow (which is repeated all over the place, in notebooks (Cause you want to plot the time series sometimes), in seir.py (function called by one_simulation_legacy and the simulate script) and inference.py:one_simulation (for classical inference) and inference.py:simulation_atomic (for emcee inference)

# draw some parameters
# parameters are drawn > parameters array, same for all subpop and times (except if the parameter is a time series), an array of dimensions n_parameters, n_times, n_subpop (unsure of the order)
p_draw = modinf.parameters.parameters_quick_draw(n_days=modinf.n_days, nsubpops=modinf.nsubpops)
# p_draw > SEE SCREENSHOT 2!!! A  SINGLE LINE

# build a modifier object
npi_seir = seir.build_npi_SEIR(modinf=modinf, load_ID=False, sim_id2load=None, config=config)

# apply the NPI to the parameters
# Modifiers are applied, this changes the array (not the dimensions) so parameters varies in time and subpop according to their modifiers
reduced_parameters = modinf.parameters.parameters_reduce(p_draw, npi_seir)
# reduced_parameters > SEE SCREENSHOT 2!!! A TIMESERIE

# parse the compartments:
unique_strings, transition_array,proportion_array, proportion_info  = self.modinf.compartments.get_transition_array()

# parse the unique_string and compute the parsed parameter objects
# Then, the compartments module computes the way to aggregate parameters for the sampler, e.g if on transition has a rate of `beta * # phi / gamma` we do NOT want to pass beta, phi, and gamma to the simulator because
#  * it's a jit compiled function so parsing this formula there is a big hassle and waste of time
#  * we don't use these values separately anyway
#  So instead compartments compute unique_strings (each one being a formula, a constant, or a a parameter name), and then if you give a parameter array (as above, modified by NPI) to the compartments module it will pre-compute these unique_strings (which are different parameters in some sort) by parsing each formula. This array is called, bad name, parsed_parameter and has size n_unique_strings, n_days, n_subpop
parsed_parameters = modinf.compartments.parse_parameters(reduced_parameters, modinf.parameters.pnames, unique_strings)
# parsed_parameters > SEE SCREENSHOT 3!!! now it's a time series that is computed from a formula

[pearsonca]

What would it take, conceptually, to a build / compile / etc action, after which running the model anytime would look like:

model = read_model_obj("output/compiled.thing") # or likely also an element from a config that specifies the location compiled models get sent
# maybe something like parameters = model.transform(raw_parameters) which applies all the modifiers
timeseries = model.run(parameters)

so then when I want to do simulation troubleshooting at the REPL, say, I just have parameters = dict(beta= ..., gamma = ...) by hand and I run it. In simulate via the pipeline, flepimop does some work to pick out the proper parameters, but then just runs the above, gets the timeseries back and then manages where to save what.

Calibration also looks like flepimop picking out the proper parameters, feeding model.run(...), getting back the timeseries, applying the observation transformation, calculating the likelihood, rinse and repeat, and then storing results. In this version, model doesn't know anything about which parameters are the inference targets or not.

[jcblemai]

currently this look like

# create object and stores an attribute static_sim_arguments that has 
gempyor_inference = GempyorInference(config_filepath=config_filepath) 
# find example value for parameters, here sample from config.
parameter_array = gempyor_inference.inferpar.draw_initial(n_draw=1)  

# Simulate
#.   A. If you want logoss:
logloss = gempyor_inference.get_logloss_as_single_number(parameter_array)
#.   B. If you want Timeserie instead
outcomes = gempyor_inference.simulate_proposal(parameter_array)

As you see it's very ugly, but also quite close to what you are saying. But it still needs the config because some objects are difficult to work with at the moment (you'd want to rewrite Modifiers), and it does not write the compiled thing to disk (unsure everything is pickable here, especially the JIT part/also might not be the best for memory.

In this version, model doesn't know anything about which parameters are the inference targets or not.

But then the parameter_array is huge because potentially all seeding row and initial conditions are inference target. Though I see what you mean.

[jcblemai]

EMCEE inference

entrypoint: calibrate.py

1. Build-up

it works fully through inference.py. First, you build an object inference.GempyorInference, which will hold a ModelInfo object (this nesting is really bad -- ModelInfo was though as a pickable thing but you really want the axes only: time, subpop and compartments) . The object that are created inside are:

• ModelInfo: pointer to all objects in gempyor (seeding, compartments ...)
• inference_parameter.InferenceParameters (only if emcee in config: inference:)
  • reads the config, and take each seir and outcomes modifiers with a perturbation as something to propose new values using inference.
  • Useful attribute: stores a data frame that contains, for each parameter fitted it's name, the subpop it applies to, and its type (currently only seir_modifier or outcome_modifier). This data frame is e.g the key to the HDF file.
  • useful methods: if you give it a proposal, it checks that it obeys the bounds of each of the parameters
  • limitations: does not allow for resume yet.
• logloss.LogLoss (only if emcee in config: inference:)
  • this object stores the logloss, which is a collection of Statistics
  • each statistics: a outcome column and a ground truth column, some llik function, perhaps some aggregation, and sometime a regularizer
  • logloss can return the logloss (=llik + regularizers) by subpop and statistics, or summed over one these, or as a single number (what is used for emcee)
• static_sim_argument (dictionary) (only if emcee in config: inference:)
  • created by inference.get_static_arguments(modinf).
  • This asks each module in ModelInfo for some information and stores in in the following dictionary that is picked.

    static_sim_arguments = {
        "snpi_df_ref": snpi_df_ref,
        "hnpi_df_ref": hnpi_df_ref,
        "p_draw": p_draw,
        "unique_strings": unique_strings,
        "transition_array": transition_array,
        "proportion_array": proportion_array,
        "proportion_info": proportion_info,
        "initial_conditions": initial_conditions,
        "seeding_data": seeding_data,
        "seeding_amounts": seeding_amounts,
        "outcomes_parameters": outcomes_parameters,
    }

2. Usage

The function to simulate a proposal is GempyorInference.simulate_proposal, which takes as input a proposal (numpy array or list of le/size gempyor_inference.InferenceParameter.get_dim(), the number of fitted parameters), then it:

1. uses InferenceParameter.inject_proposal to infect the proposal array into the dataframes used to build modifiers (snpi_df_ref and hnpi_df_ref) (Ideally, we want each object to handle the proposal injection instead of InferenceParameter to know how each object is using the proposal as it is now)
2. then uses the worker that computes a simulation: outcomes_df = simulation_atomic(** static_sim_arguments, modinf=self.modinf, save=self.save) with save as a flag if we need to save to disk in classical inference format all files (snpi, hosp, seir, ...). simulation atomic does a few things (convert seeding from Dict to NumbaDict, run the seir simulation, compute the outcomes) using the static parameters (so you cannot change seeding and all). This returns the outcome dataframes.

So GempyorInference.simulate_proposal gets a proposal and returns an outcome data frame (and, if called with save=True, save everything on disk). There are wrappers around that functions:

• GempyorInference.getLogLoss which runs the proposal then return total logloss ll_total, logloss array logloss, and just the regularizations.
• GempyorInference.get_logloss_as_single_number which just return the total logloss as emcee just need one number to accept and reject.

Hence one call EMCEE using

          sampler = emcee.EnsembleSampler(
              nwalkers,
              gempyor_inference.inferpar.get_dim(),
              gempyor_inference.get_logloss_as_single_number,  # function that takes a proposal and return a logloss
              pool=pool,
              backend=backend,         # backend is the place where it'll save the HDF file.
              moves=moves,
          )

Notes

This is all done in object, without performance penalty (in inference.py, there are commented lines to call get_logloss_as_single_number without using the gempyor_inference object at all anymore (using it only to build the static arguments and the different object). It did not seem faster from benchmark but might be preferred.

[jcblemai]

Classical inference

It works fully through inference.py. First, you build an object inference.GempyorInference, which will hold a ModelInfo object:

• ModelInfo: pointer to all objects in gempyor (seeding, compartments ...)

Step 1. initialize the files

If not resume, R inference calls:

gempyor_inference_runner$one_simulation(
          sim_id2write=0,
          load_ID=Fase)

to get one run from config. This writes index 0 of all hnpi, snpi, spar, hosp, seir, ... files in the outputs.
If resumes, these files are taken from the resume location

Step 2: generate proposal

R inference reads the snpi, seed and hnpi files, and writes them again but this time putting proposed parameters inside, then call

  gempyor_inference_runner$one_simulation(
          sim_id2write=this_index,
          load_ID=TRUE,
          sim_id2load=this_index)

So gempyor does one simulation reading the files that inference has modified, and writes the same exact files, but this time also writing an hosp and seir files containing the outcomes

Inside gempyor, the first time one_simulation is run, it builds the compartments structure only (this can be quite long), while all the rest (spar, snpi, ...) is read from the file at this_index

the rest

Rinference read the hosp file, compute the llik and accept/reject it's proposed parameters. This is in a loop that update this_index

[jcblemai]

The 3 ways to run everything

There are 3 ways to run a full gempyor simulation. A reminder that this means: take a config, or a config + some files, and build all modelInfo object, run a RK4 integration to get all the model states in time, then run the observations model to get an hosp file/outcomes).

Ideally these 3 way are consolidated into one so there are no bugs. The first way does not use inference.py (which could be called the interface file of gempyor), the two other does.

The way simulate.py does it:

# [create a modelInfo]
seir.run_parallel_SEIR(...)
outcomes.run_parallel_outcomes(...)

Simply runs the two "runner functions" from the seir and outcomes module. These functions build modifiers, compartments, and all. In inference.py, inference.py:one_simulation_legacy also uses these two functions only to run a module

The way classical inference does it

It uses inference.py:one_simulation which is better than the above because it can of several simulations, the built-up is done only once, and it saves intermediate steps in memory for debug. This can read or not some input files.

The way emcee does it

It uses inference.py:simulate_proposal (more precisely, inference.py:get_logloss_as_single_number) where most arguments are pickled and just the proposal is inputed.

[pearsonca]

Thanks for this.

I think the keys here are:

• figure out how to pickle an object that works a parameters => outputs mapping transformation and does nothing else
• figure out how to represent inference as working on a reference set of parameters, patching those with sampling draws for the fitted bits, and then applying the transformer => evaluating outputs => minimizing some loss function => repeat (ideally cleverly)
• pulling file handling out of these activities and making that the sole domain of the pipeline

[jcblemai]

• figure out how to pickle an object that works a parameters => outputs mapping transformation and does nothing else

In what sense pickle? There are a few methods that don't depend on their self so do you mean to group them together as one method ? You need as ingredients from the static_sim_arguments p_draw (parameters, can be varied), [s/h]npi_df_ref (Dataframe), unique_strings (list), transition_array (array), proportion_array (array), proportion_info (array). The last 4 stores the transformation of parameters spatial timeseries (created from p_draw + npi_df) into the shape and format used by transition in the model. So you want to pickle all that (like in static_sim_arguments and take the functions highlighted above)

• figure out how to represent inference as working on a reference set of parameters, patching those with sampling draws for the fitted bits, and then applying the transformer => evaluating outputs => minimizing some loss function => repeat (ideally cleverly)

Yes, there is some work to do here. Right now inferenceParameter inject the parameters into the [s/h]npi_df_ref dataframes (because Modifier code is a mess). However, I believe each object should be able to handle an inference payload (wether plug-in or not) and deal with it with a fast function to update it's state.

• pulling file handling out of these activities and making that the sole domain of the pipeline

Agree. Right now emcee writes the hdf files but we can use our own backend-class and i think that would be the way to go