Lb/kinetic model parameter changes

[lubo93]

"simulate_simple_epidemic.py" and "scenarios.py" now contain the new NYC data

Some first test showed that the kinetic model produces meaningful results for the updated parameters.

However, changing the maximum in "DiurnalMeanContactRate" from 22 to 3 didn't change the results. Could you please check if changes in the DiurnalMeanContactRate affect the simulation results?

[glwagner]

@lubo93 I have pushed some changes to fast_contact_simulator.py, though I'm not sure they solve your issue. We also will either have to remove contact_simulator.py, or refactor the code so that we don't have to maintain both.

I created an example called changing_max_contact_rate.py. It produces this plot:

The typical number or active contacts and the mean contact duration does trend as expect. But is the value of the mean contact duration correct? Note that this is the average contact on each edge, averaged over all edges. So I suppose it may be small because many edges are inactive.

I haven't run your example yet, but I'll give that a shot.

[glwagner]

An important note is that the changes I think (I hope) the changes I pushed ensure that the ContactSimulator returns correct results even when rates are very long, such that the time-step is longer than the mean contact interval. In other words, the discrete nature of the Gillespie stepping should eventually become apparent as samples are returned with ever greater frequency.

The new example changing_mean_event_duration.py explores changing mean event durations. If the mean event duration is very long, we see that mean contact durations "relax" following a midday peak:

Equilibrium is not reached at night.

[glwagner]

@lubo93 this is what I get with

mean_contact_rate = DiurnalMeanContactRate(maximum=40, minimum=2),

whereas with

mean_contact_rate = DiurnalMeanContactRate(maximum=35, minimum=2),

and with

mean_contact_rate = DiurnalMeanContactRate(maximum=34, minimum=2),

but with

mean_contact_rate = DiurnalMeanContactRate(maximum=33, minimum=2),

the epidemic dies! (Well, for this choice of random seed.)

Smaller values of maximum produce aborted epidemics. This seems reasonable? We have a small network and the initial number of infected is merely 2.

As a side note, we also have some pretty printing, as inspired by @lubo93:

Epidemic day:  15.247, wall_time:  0.142 s, mean(w_ji):   2 min, statuses:  S  295 | E  220 | I  176 | H    5 | R  302 | D    2 |
Epidemic day:  15.364, wall_time:  0.316 s, mean(w_ji):   4 min, statuses:  S  292 | E  217 | I  180 | H    5 | R  304 | D    2 |
Epidemic day:  15.499, wall_time:  1.569 s, mean(w_ji):  22 min, statuses:  S  284 | E  219 | I  174 | H    6 | R  315 | D    2 |
Epidemic day:  15.621, wall_time:  2.317 s, mean(w_ji):  33 min, statuses:  S  272 | E  225 | I  172 | H    6 | R  323 | D    2 |
Epidemic day:  15.744, wall_time:  2.391 s, mean(w_ji):  34 min, statuses:  S  261 | E  229 | I  173 | H    6 | R  329 | D    2 |
Epidemic day:  15.869, wall_time:  1.606 s, mean(w_ji):  23 min, statuses:  S  254 | E  224 | I  178 | H    6 | R  336 | D    2 |
Epidemic day:  15.979, wall_time:  0.292 s, mean(w_ji):   4 min, statuses:  S  251 | E  223 | I  175 | H    7 | R  342 | D    2 |
Epidemic day:  16.102, wall_time:  0.149 s, mean(w_ji):   2 min, statuses:  S  250 | E  219 | I  177 | H    7 | R  345 | D    2 |
Epidemic day:  16.216, wall_time:  0.136 s, mean(w_ji):   2 min, statuses:  S  250 | E  215 | I  179 | H    7 | R  347 | D    2 |
Epidemic day:  16.367, wall_time:  0.288 s, mean(w_ji):   4 min, statuses:  S  248 | E  210 | I  182 | H    7 | R  351 | D    2 |
Epidemic day:  16.496, wall_time:  1.579 s, mean(w_ji):  22 min, statuses:  S  237 | E  209 | I  186 | H    7 | R  359 | D    2 |
Epidemic day:  16.616, wall_time:  2.351 s, mean(w_ji):  33 min, statuses:  S  230 | E  210 | I  178 | H    7 | R  373 | D    2 |
Epidemic day:  16.749, wall_time:  2.561 s, mean(w_ji):  34 min, statuses:  S  225 | E  207 | I  180 | H    7 | R  379 | D    2 |
Epidemic day:  16.874, wall_time:  1.833 s, mean(w_ji):  22 min, statuses:  S  220 | E  203 | I  184 | H    8 | R  383 | D    2 |
Epidemic day:  16.997, wall_time:  0.326 s, mean(w_ji):   4 min, statuses:  S  218 | E  194 | I  193 | H    8 | R  385 | D    2 |
Epidemic day:  17.125, wall_time:  0.149 s, mean(w_ji):   2 min, statuses:  S  217 | E  188 | I  199 | H    8 | R  386 | D    2 |
Epidemic day:  17.236, wall_time:  0.143 s, mean(w_ji):   2 min, statuses:  S  217 | E  178 | I  204 | H    8 | R  391 | D    2 |
Epidemic day:  17.367, wall_time:  0.297 s, mean(w_ji):   4 min, statuses:  S  216 | E  176 | I  201 | H    8 | R  397 | D    2 |
Epidemic day:  17.498, wall_time:  1.562 s, mean(w_ji):  22 min, statuses:  S  210 | E  173 | I  201 | H    7 | R  406 | D    3 |
Epidemic day:  17.621, wall_time:  2.369 s, mean(w_ji):  34 min, statuses:  S  201 | E  174 | I  200 | H    7 | R  415 | D    3 |
Epidemic day:  17.741, wall_time:  2.530 s, mean(w_ji):  34 min, statuses:  S  194 | E  177 | I  193 | H    6 | R  426 | D    4 |
Epidemic day:  17.849, wall_time:  1.777 s, mean(w_ji):  23 min, statuses:  S  194 | E  170 | I  190 | H    6 | R  436 | D    4 |
Epidemic day:  17.998, wall_time:  0.316 s, mean(w_ji):   4 min, statuses:  S  193 | E  165 | I  189 | H    6 | R  443 | D    4 |
Epidemic day:  18.110, wall_time:  0.148 s, mean(w_ji):   2 min, statuses:  S  193 | E  163 | I  185 | H    6 | R  449 | D    4 |
Epidemic day:  18.232, wall_time:  0.139 s, mean(w_ji):   2 min, statuses:  S  191 | E  158 | I  184 | H    5 | R  458 | D    4 |
Epidemic day:  18.360, wall_time:  0.306 s, mean(w_ji):   4 min, statuses:  S  190 | E  153 | I  186 | H    6 | R  461 | D    4 |
Epidemic day:  18.497, wall_time:  1.716 s, mean(w_ji):  22 min, statuses:  S  187 | E  152 | I  183 | H    7 | R  467 | D    4 |
Epidemic day:  18.615, wall_time:  2.772 s, mean(w_ji):  33 min, statuses:  S  177 | E  158 | I  178 | H    8 | R  475 | D    4 |

Finally, I found that the fast version of the contact simulator produces different results from the "slow" contact simulator for the same setup. We need tests for these things to nail this down. The random number generator may, in fact, be different with numba, which is annoying. But that's for a future PR.

[lubo93]

Ok, this looks good and it's good that smaller values of lambda_max constrain the outbreak. I will do some further tests.

[glwagner]

By the way, the pretty printing above illustrates #57