Msprime: Performance benchmarking for GC

[jeromekelleher]

Current GC implementation takes a simple approach to computing some rates. run some benchmark simulations and see how much of the total time this takes up using Linux perf.

[jeromekelleher]

It would be good to do some perf benchmarking here @fbaumdicker, to see what the cost of the cleft hack is. Ideally we want to run some large sims through a profiler like Linux perf so that we can see the amount of time spent in each function. Is this something you'd like to take on?

[fbaumdicker]

I was rather thinking to look into the no recombination but gene conversion >0 issue. But as long as you are looking into #852, I will have a look at the performance for larger datasets. Do you have any suggestions for the scenarios to look at? Is there some setting you already considered we can compare to?

I was thinking about a thousand individuals with long sequences with a rather short gene conversion track length. Let's say:

recombination rate 1
gene conversion rate 5
gene conversion track length 10
sequence length 10000

Any other suggestions?

[jeromekelleher]

I'm not sure @fbaumdicker, whatever you think would be a few different types of very large GC simulations.

[fbaumdicker]

For a large number of gene conversion events, the time grows quadratic.
The number of gene conversion events itself grows linear with sequence length.

geneconversion_perf.pdf

Also, 99% of the time is spent in get_cleft_total. So we should consider speeding up this function. Maybe with another Fenwick tree.

[jeromekelleher]

Thanks for this @fbaumdicker --- looks like we do need another Fenwick tree to keep track of cleft_total then, right?

[fbaumdicker]

Yes, indeed. I will have a look into that.

[jeromekelleher]

I've added this to the feature list for 1.0 @fbaumdicker, as I think we should be able to fix the perf issues easily enough using a Fenwick tree. What do you think?

[fbaumdicker]

Yes, the Fenwick tree will solve the perf issues. The main issue is currently that we can not simply copy the existing Fenwick tree framework, where each node corresponds to a segment. Instead, we either need a tree where each individual (i.e. a chain of segments) corresponds to a value/treenode or we let the tree node corresponding to the first segment of the chain record the cleft value and set all other values to zero. I currently prefer the second approach.

[jeromekelleher]

The tree would be much smaller (and faster) though if just tracked individuals, wouldn't it? We're going to need an 'individual_t' struct soon anyway, so we could make it have a similar setup to segment_t. If you get something working in the Python version, we can see how it looks?

Also, I think we need a better name than 'cleft'. What does this even mean?? (Perhaps this should spawn a separate issue)

[fbaumdicker]

Yes, it would be much smaller, but it would be more often expanded and decreased in size as each recombination creates a new individual but not necessarily a new segment. An individual_t struct would be really handy for this. Although in the case of gene conversion we will need to link the segments and the individuals as changing the segment length will change the cleft value of the individual in a nontrivial way.

The name "cleft" originates from the ms code and means probably something like "conversion from the left". It is the relative ratio of gene conversion events that start left of the individuals' first segment and do not span the complete chain of segments, which would make the event irrelevant. To compute it we need the value left of the left-most segment_t and the value right of the right-most segment_t and the mean track length parameter.

[jeromekelleher]

Well, have a go with it in Python and we can see what works best before digging in to the C.

[jeromekelleher]

Just doing some some perf benchmarking on the latest code merged in #1196 @fbaumdicker and ... you've nailed it!!! We're easily able to to run simulations at the human chromosome scale - it's basically just as efficient as recombination now. I've run various benchmarks, and all the time is being spent in the low-level stuff like msp_merge_ancestors and fenwick_set_value. As far as I can tell, the time required for running GC events now is basically negligible.

This is a massive step forward, and the culmination of a huge amount of work. Pour yourself a beer!

[fbaumdicker]

Awesome! I hope we can reach a similar performance for non-constant gc_maps.