Integrated ContactSimulator and KineticModel simulation

[glwagner]

This PR integrates ContactSimulator into examples/simulate_kinetic_model.py and extends the example to include a growth phase and a socially-distanced phase.

It also refactors the ContactSimulator and the KineticModel somewhat, and nukes contact_generator.py since it isn't needed anymore.

The example can use some work. For one, I get a KeyError from inside Gillespie_simple_contagion for some random seeds. What causes this? Something must be wrong about the inputs.

The annoying thing is that the simulation runs smoothly sometimes, too (as it does for the current choice of random seed).

Other than that things look fairly nice.

[odunbar]

This looks clean, thanks Greg!

For one, I get a KeyError from inside Gillespie_simple_contagion for some random seeds.

This is strange... does it occur on the first solve (before we perform swapping etc with the hospital beds)

[lubo93]

Ok, thanks a lot. I will test it and see if I can reproduce the KeyError. If it's a KeyError from Gillespie_simple_contagion, it's probably the initial condition/status list. Depending on how you change the statuses, it could be the case that you're missing some nodes in the IC. This can be checked by comparing the node IDs of the network with the node IDs in IC.

[glwagner]

@lubo93, I'm not sure I completely understand, but should we write a check in Gillespie_simple_contagion that throws an error if the input is invalid?

[odunbar]

if the input is invalid?

But what do you mean by invalid? As far as I see we always produce a legitimate input of random infections...

An invalid input is either of the wrong size, or if we input some invalid value (not a str, or not one of the letters S,E,I,H,R,D,P) or an invalid key ( not an int, or out of range?). And why would this only appear for different random seeds? As the only thing that is done randomly in the initial condition is a choice of which 10 values = 'I'.

I cannot see how this involves the ICs. But rather it something in the Gillespie method. I think it best we investigate in another PR so we can merge this?

[glwagner]

Hmm yes... I agree it is some kind of problem with the Gillespie method, but the issue is that I can't tell if its invalid input or not. It's possible that invalid input still "works" for some epidemic trajectories, but not for others (somehow).

Another possibility is that there is a bug in simulation.py.

[glwagner]

Well, it looks like len(statuses) and len(contact_network) are the same...

[glwagner]

I've added some verbose printing of the initial state. There is a node in the contact network for every entry in statuses, and there is an entry in statuses for every node in contact_network. This seems like it might be a correct initial condition?

[glwagner]

What happens when the epidemic dies? Does Gillespie_simple_contagion automatically quit, or does it throw an error?

[odunbar]

What happens when the epidemic dies? Does Gillespie_simple_contagion automatically quit, or does it throw an error?

This is a possible thing. I can imagine the timestep for Gillespie blows up to infinity.

[lubo93]

I am working on it and try to do some checks. What I couldn't find is the function that keeps track of the evolution of the SEIRHD compartments. I used something like:

"initialization" time_arr = [] states_arr = {} for s in return_statuses: states_arr['%s'%s] = []

... times, states = res.summary()

for s in return_statuses: states_arr['%s'%s].extend(states['%s'%s])

time_arr.extend(times+i*deltat)

It would be nice to have access to time_arr and states_arr to plot something:

plt.plot(time_arr, states_arr['D'])