Msprime: API for running simulations in parallel?

[jeromekelleher]

If #1257 gets resolved and we have a deterministic random seed for each replicate simulation, then we have a fairly straightforward path to running simulations in parallel. I'm guessing we could just add a num_threads argument to sim_ancestry, and still return an iterator over the simulated tree sequences. E.g.

num_replicates = 1000
tmrca = np.zeros(num_replicates)

for index, ts in enumerate(msprime.sim_ancestry(10**7, num_relicates=num_relicates, num_threads=8)):
     tree = ts.first()
     tmrca[index] = tree.time(tree.root)

Behind the scenes I think we can do all the work in simulating and creating a tree sequence in worker threads, making them available (in arbitrary order) to the main thread as they are ready.

By using threads, there will also be GIL issues, but I think this is probably a second-order effect, as if the simulations are small enough that the time we spend in Python becomes significant, then this pattern is useless anyway. We could also use processes, but I wonder if the overhead of either pickling the resulting tree sequences, or working around the need to pickle them is actually worth it (but, now we can serialise to an FD, maybe this is simpler than I think). In either case,if we use threads or processes the top-level API would be the same I think.

The question is, is this worth it? For this pattern to actually be useful we'd need to be in a fairly sweet-spot in which the processing time is dominated by the actual simulation, and the processing time for each replicate is small. I've chosen the example above as a case where this would probably be true.

Would we be better off just documenting some example workflows where both the simulation and analysis are done in threads? (Like #59 originally advocated?)

What about Dask integration: perhaps our time would be better spent making this possible, rather than providing our (limited) threading APIs?

Any thoughts @molpopgen, @benjeffery, @grahamgower, @petrelharp ?

[jeromekelleher]

Looks like doing simulations using Dask futures would be basically the same as the ideal recipe for using concurrent.futures. So, as in nearly all cases the user will probably want to do the analysis of the resulting tree sequence in the context of the worker process, I think we're probably better off focusing efforts on a "Running simulations in parallel" tutorial which would show how to use both concurrent.futures and Dask futures, and not bother with the num_threads API at all.

[molpopgen]

I can say that working via processes speeds up a lot of my work flows. I don't think that the pickle step has been a problem. Sometimes I return a tree sequence, but I may also just return the result of an analysis.

While it is appealing to try to support threads, it seems hard!

[jeromekelleher]

While it is appealing to try to support threads, it seems hard!

I think it'd work well - I've done it a lot with tsinfer, and it all works nicely once the low-level C objects are independent. As in, we create one Simulator and TableCollection object per-thread and only access them from that thread. The final tree sequence returned to the calling thread is then an independent copy of the tables.

[molpopgen]

While it is appealing to try to support threads, it seems hard!

I think it'd work well - I've done it a lot with tsinfer, and it all works nicely once the low-level C objects are independent. As in, we create one Simulator and TableCollection object per-thread and only access them from that thread. The final tree sequence returned to the calling thread is then an independent copy of the tables.

So does this work by letting python manage the thread creation/destruction, this avoiding having to worry about threading framework?

[grahamgower]

I think we're probably better off focusing efforts on a "Running simulations in parallel" tutorial which would show how to use both concurrent.futures and Dask futures, and not bother with the num_threads API at all.

I tend to agree with this sentiment. Consider another (common?) use case: simulating replicates with some parameter drawn from a distribution. It would be ludicrously complicated to support such flexibility in the API (each param must now accept a list!), but this is a trivial change when using a wrapper. In this case, it can also be useful to have the replicates returned in the same order as the parameters were drawn, which concurrent.futures' map() will do by default.

[benjeffery]

Agree with the concurrent/dask futures tutorial as a way to go.

If we were to do num_threads or other dev work, as it isn't a breaking change would it be best left to a post-1.0 release? Otherwise I would feel that 1.0 is feature-creeping.

[jeromekelleher]

Closing this in the light of #1287 as it turns out we can't really do entirely independent replicates without hurting performance.