Contact-specific, non-stationary birth/death simulation with rejection sampling and a Gillespie algorithm

[glwagner]

This PR adds a script that attempt to implement independent birth-death simulations of inception and termination for each contact independently. A loop over all contacts is multithreaded via numba.

The script currently runs, but the results are not correct. The incorrectness of the results may also influence the time it takes to run the script (which would otherwise be impressive). EDIT: the results are "more correct" now, though probably still not completely correct due to the time-dependence of the rates. See below. Any input is appreciated. The issue could be a bug, or a conceptual error. I'll keep trying to debug this tomorrow. It seems this is what we need to implement node-specific interventions, via the discussion on #65 .

Second edit: the winning algorithm seems to be rejection sampling, at the moment.

[glwagner]

Here is some analysis of a 1,000,000 contact simulation of individual contacts with a naive Gillespie method. The top panel shows the measured inception rate on a single contact. The bottom panel shows the contact duration averaged over all contacts.

The temporal variation in the mean contact duration is markedly different from the "bulk" Gillespie simulation.

But maybe it's ok, since the measured inception rate is not that far off from what we'd like, in fact.

It's notable that this method is far faster than simulating all the contacts at once. Simulating 1,000,000 edges for 4 days takes just 1.6 seconds on my laptop. This method is easily adapted to permit different inception rates on every edge.

[glwagner]

I have added a new script that performs rejection sampling for each contact. Rejection sampling takes short, constant time steps, sampling along the way to determine whether transition occurs.

The algorithm takes longer than a single Gillespie simulation for a contact due to the need for more steps, but is much shorter than the "bulk" Gillespie simulation for all contacts at once because it avoids expensive array operations (presumably).

Since each contact is simulated, the loop over them can be parallelized with multithreading, and we can support independent rates for each contact.

On my laptop, simulating 1,000,000 contacts with a 10 second microstep for a week takes about 7 minutes of wall time. It produces this plot:

I have also attempted to accurately simulate the non-stationary, single-contact process by modifying the "naive" Gillespie algorithm. In general, I think this requires solving an integral equation for the waiting time (though I admit I am not 100% positive of this conclusion). For this I used an iterative method with Gaussian quadrature. However, my attempts were not successful. It is possible that the iteration is not stable.

Rejection sampling seems like the simplest approach.

I found this reference helpful:

https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1004579

Feedback welcome, including other ideas. Rejection sampling seems simple enough, however. I will work on integrating this script into a new ContactSimulator permitting contact-dependent rates, unless there are objections.

@tapios @lubo93 @odunbar

[tapios]

Thanks, @glwagner! This is looking good and I think suffices for what we need. Solving ODEs for the probabilities is another fast alternative, but if this works well enough, and it seems it does, it's preferable (because of the randomness it leaves in the process).

[odunbar]

For this sampling you need a timestep - is 10 seconds the maximum you've found for decent results? This clearly works, as Tapio said we could do an ODE solve too. Though I think this is not a bottleneck any longer. Great job!

[glwagner]

Here’s a plot generated from a direct simulation using a rejection sampling method over 10,000 contacts. The interval time is 1 hour and the simulation is run for 7 days.

The top panel shows the measured number of inceptions, divided by the interval length (thus, an inception rate) for four contacts. The discrete nature of this panel reflects that, at most, 3 inceptions are observed in an interval for the chosen set of edges.

The middle panel visualizes the variability in mean contact duration over an interval (here, 1 hour) that can occur across edges. There is a burst of activity midday. The variability in the number of inceptions in an interval, plus variability in the lifetime of each event, creates variability in the contact duration over an interval for each edge.

The bottom panel plots the “ensemble-averaged” mean contact duration averaged over all edges. This closely follows lambda(t) / [mu + lambda(t)] , provided t is evaluated at the mid-point of each interval. I suppose this is some kind of quasi-stationary solution that’d we’d obtain if we assumed the birth-death process rapidly relaxed to a quasi-equilibrium over a short timescale associated with mu that changes slowly due to variations in lambda.

[lubo93]

Here’s a plot generated from a direct simulation using a rejection sampling method over 10,000 contacts. The interval time is 1 hour and the simulation is run for 7 days.

The top panel shows the measured number of inceptions, divided by the interval length (thus, an inception rate) for four contacts. The discrete nature of this panel reflects that, at most, 3 inceptions are observed in an interval for the chosen set of edges.

The middle panel visualizes the variability in mean contact duration over an interval (here, 1 hour) that can occur across edges. There is a burst of activity midday. The variability in the number of inceptions in an interval, plus variability in the lifetime of each event, creates variability in the contact duration over an interval for each edge.

The bottom panel plots the “ensemble-averaged” mean contact duration averaged over all edges. This closely follows lambda(t) / [mu + lambda(t)] , provided t is evaluated at the mid-point of each interval. I suppose this is some kind of quasi-stationary solution that’d we’d obtain if we assumed the birth-death process rapidly relaxed to a quasi-equilibrium over a short timescale associated with mu that changes slowly due to variations in lambda.

Great that this approach works! All figures are based on the "integration method", right? What is the performance (simulation time) for 10^6 contacts and 1 day?

What I mean by integration method is this piece: integral += w[i] * diurnal_inception_rate(λmin, λmax, t + x[i] * Δt). Is this what you used for the quoted figures?