Random Draws For SEIR Parameters Per Slot

[TimothyWillard]

Describe your changes.

This pull request addresses the bug described in GH-406 where new seir parameters were not drawn across slots.

Does this pull request make any user interface changes? If so please describe.

The user interface changes are:

• Added a new Parameters.reinitialize_distributions method that recreates the randomly drawn seir parameters so that they do not hold onto the prior random state, and
• SEIR.run_parallel_SEIR will now draw random seir parameters across jobs.

These changes have not been documented yet, @alsnhll where do you think is the best place to describe this behavior (I'll make these changes before marking this PR as ready for review)?

What does your pull request address? Tag relevant issues.

This pull request addresses GH-406.

Tag relevant team members.

@alsnhll

[alsnhll]

I did some testing of how random draws of parameters work for each parameter type, using the config config_sample_2pop_modifiers_test_random.yml which is now in the tutorials folder

In the original code in the main branch, the behavior is:

• spar file, seir:parameters:Ro - same value by slot - Not working as intended
• hpar file, outcomes:incidDeath:delay - varying by slot and by location. Only taking on integar values.
• snpi file, seir_modifier:modifiers:Ro_relax - varying by slot and by location.
• hnpi file, outcome_modifiers:modifiers:test_limits - varying by slot and by location.

Thus variation in spar was the only one not working as intended, while variation in hpar, snpi, and hpi seem to be working as intended.

In this branch, the behavior is:

• spar file, seir:parameters:Ro - varying by slot
• hpar file, outcomes:incidDeath:delay - varying by slot and by location. Only taking on integar values.
• snpi file, seir_modifier:modifiers:Ro_relax - varying by slot and by location.
• hnpi file, outcome_modifiers:modifiers:test_limits - varying by slot and by location.

So now variation in spar seems to be working as intended, and nothing has changed about how hpar, snpi, and hnpi work.

However note that users are likely surprised by the variation by location in hpar, snpi, and hnpi in this case, as I don't think this behavior is documented anywhere. It's likely not clear to most users that these parameters have different values by location by default when they randomly vary, even though when their value is fixed they have the same value by location. Is this what we want the default behavior to be? I know during inference values are fit by location by default but does that make sense as the default for forward simulations? It's also confusing that spar parameters don't vary by location while these other ones do (I know that's a long standing flepimop issue but it re-confuses me every time I encounter it).

If I instead run a config where for both seir and outcome modifiers I use subpop_groups: "all" - config_sample_2pop_modifiers_test_random_subpop_groups.yml- (which should force all subpopulations to have the same parameter value, at least for inference), then I get

• snpi file, seir_modifier:modifiers:Ro_relax - varying by slot but same across locations.
• hnpi file, outcome_modifiers:modifiers:test_limits - varying but same across locations.

which is what we expect. I don't know how to make hpar parameters be grouped by location - can we do this?

To summarize, yes this PR seems to fix the main underlying problem, but this issue has led me to notice a behavior that needs to be documented better and discussed - during forward-only simulations, when hpni, snpi, and hpni are set to be drawn randomly from distributions, they very by slot and location, while spar varies only by slot.

Thoughts from @jcblemai, @shauntruelove, @saraloo on the desired behavior here would be helpful.

[pearsonca]

Nice investigation.

Agree we need to clearly document what matches / doesn't when stochastically simulating. Where in the gitbook do you think it's best to capture that?

Regarding behavior, all the stochastic features should match in the same way in my opinion. Weird that spar seems to deviate, so we should figure out why that is. I kind of think that location shouldn't cause a different sample either, but it's most important that we have consistent behavior.

[jcblemai]

I agree Alison on consistency. The current state is not motivated by anything but historical hastly made decisions during COVID (fixed R0 across subpops, but age structure different in each state). We should change that carefully.

Now:

• @alsnhll I don't see this config in the repository (perhaps I am in the wrong branch), does the outcome section has a parameter modifying file (another weird thing).
• @TimothyWillard An opinion: While this addresses the issue, I believe the underlying problem is the seed issue (see my comment here [Bug]: Sampling parameters from distributions not working as intended #406 (comment): python subprocesses inherit the random process state), and just resetting the distribution leaves other possible bugs:
  • One would expect seeds to be different for each slot -- we can talk about handling that deterministically for reproducibility (but right now it is not done deterministically even if it is the same) but at least different simulations should not share seed.
  • So in addition to allowing other random generators to still be in sync (e.g if one does stochastic simulations), it is pretty fuzzy (for me though, perhaps you got it) how it would work without being at the mercy of a scipy distribution implementation change (right now created distributions uses their own seed ? what happens here that makes the number different) ?
  • Even if we keep this update distribution approach, we need to seed each correctly with a different seed (at least like the emcee, ideally deterministically).

[alsnhll]

Configs are here:
https://github.com/HopkinsIDD/flepiMoP/blob/feature/406/random-draws-per-slot-for-SEIR/examples/tutorials/config_sample_2pop_modifiers_test_random.yml
https://github.com/HopkinsIDD/flepiMoP/blob/feature/406/random-draws-per-slot-for-SEIR/examples/tutorials/config_sample_2pop_modifiers_test_random_subpop_groups.yml

No parameter modifying outcome file

[saraloo]

Hmm yes thanks for this investigation @alsnhll , agree with @jcblemai it's just historical decision making responsible for these inconsistencies. I'm not entirely sure what we should define here...

I see a few arguments either way

• though inconsistent, I think defining spar parameters to be the same across locations makes sense given our typical application of subpops and strata. If these are considered more as biological parameters (eg infectious period) it makes sense to me to keep them as consistent across locations.
• Alternatively one could argue it's possible to think of these as not necessarily biological factors but indeed varying across location (eg contact parameters varying across populations).
• following the same logic then, it's not immediately clear to me why the hpar parameters vary by location. How we define these are indeed a bit muddied - can be combination of both biological severity and reporting etc.

Could use some discussion wrt how to harmonise these definitions. In any case, needs to be made clear in docs etc.

I don't know how to make hpar parameters be grouped by location - can we do this?

Similar inconsistencies in the npi parameters too as alison points out - ie we can group seir_modifiers but not outcome_modifiers. Again this was just development that was implemented at the time in response to something we needed to model, but can work on including this feature to all types of modifiers.

[jcblemai]

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.

[pearsonca]

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

[jcblemai]

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)

[TimothyWillard]

I apologize for the confusion @pearsonca and @jcblemai, but this PR has not been ready for code review but now is. I was working through the GH action failure on c828f1e, which was challenging to debug due to differences in behavior between Mac+Windows vs Linux.

The quick version is both gempyor.seir.run_parallel_SEIR and gempyor.outcomes.run_parallel_outcomes are thin wrappers around their non-parallel counterparts that use tqdm.contrib.concurrent.process_map to parallelize them. This function from tqdm creates a concurrent.futures.ProcessPoolExecutor with the default multiprocessing start method which is spawn for Mac+Windows and fork for Linux.

Reinitializing the parameters distributions was sufficient for addressing the original issue under spawn, but did not address the issue under fork. In fact, the unit tests suggest that sampling for hpar has not been working correctly on linux (like the HPCs). I've added the following unit tests in this PR to demonstrate these issues and allow for testing cross platform:

• Direct testing of Parameters.reinitialize_distributions which reset the random seed held onto by the random samplers in the Parameters class under spawn,
• Overall testing of the flepimop simulate CLI by running it and then manually checking the outputted files for the desired behavior, and
• (new since c828f1e) testing of gempyor.seir.run_parallel_SEIR and gempyor.outcomes.run_parallel_outcomes by running a separate process, force setting the multiprocessing start method with multiprocessing.set_start_method, and then checking the outputted files (same as above point).

There is some duplicated code across these tests that I know I need to consolidate/clean up, but beyond that this is now ready for review. I will also note that I think this PR is a temporary fix to this problem, the correct longer term solution would be to explicitly manage the generator that numpy uses to create the samples as numpy recommends (np.random.seed is considered legacy).

@alsnhll @saraloo To clarify:

So now variation in spar seems to be working as intended, and nothing has changed about how hpar, snpi, and hnpi work.

I'm taking this (and other commentary) to mean that we need similar tests for snpi and hnpi (I suspect that we may not be sampling correctly on Linux due to the above)? I think the usage of subpop_groups: "all" address the location variation issues, but to ensure this we probably need a unit test. Are there other values that this configuration key can take besides "all" that y'all know of? Does this configuration key apply to hpar and spar? Is this a correct summary of where we want things:

	Varies By Slot	Varies By Location By Default	Effect Of subpop_groups: 'all'
spar	Yes	No	N/A
hpar	Yes	Yes	N/A
snpi	Yes	Yes	No Variation By Location
hnpi	Yes	Yes	No Variation By Location

Please feel free to correct me if I'm misunderstanding or if I've not got this table setup correctly (i.e. should be determined by some other variable or something).

[jcblemai]

Impressive, and agree for the Generator thing.

In the meantime though, I would advise random seed from OS random source or something else instead of fixing seed on sim_id which has blind spots as sim_id is used differently in inference and local run. not that important though.

I I think the usage of subpop_groups: "all" address the location variation issues, but to ensure this we probably need a unit test. Are there other values that this configuration key can take besides "all" that y'all know of? Does this configuration key apply to hpar and spar? Is this a correct summary of where we want things:

This key can be a list of list of subpop, which are then grouped together. e.g [["California", "Utah"], ["North Carolina", South Carolina"]] means that "California", "Utah" will have one value, and "North Carolina", South Carolina" another, and the rest of the subpopulation will have independent values. This is a key for modifiers

[TimothyWillard]

Work on this has been paused until after the next release due to higher priority items.