Msprime: Investigating if sims can be done in physical coordinates

[jeromekelleher]

@daniel-goldstein and I are investigating if msprime can be updated to run simulations directly in physical space rather than simulating in genetic space as we currently do. The advantage of working directly in physical coordinates is that it would get rid of a lot of annoying corner cases, especially when it comes to supporting gene conversion. The disadvantage is that it might potentially substantially slow down simulations with variable recombination rates, as we'll need to iterate over subsets of the genetic map each time we insert a segment into the Fenwick tree.

The plan is to make the minimal changes we can to try out simulating directly in physical space for the basic Hudson model, and then do a performance analysis. If it turns out that we lose a lot of performance (more than say, 50%), we go back to the drawing board/abandon the idea.

If not, then we hope to push the changes needed through the whole codebase. There may be some small semantic changes incurred, but if there are we can tag the next release 1.0 and clearly document the semantic changes incurred in the major version bump.

Cleaning up the current mess that we have due to simulating in genetic space would be a great thing to do, and an excellent start for development on the 1.0 series, I think.

cc @fbaumdicker, @andrewkern, @petrelharp, @DomNelson, @ivan-krukov

[daniel-goldstein]

I've started experimenting with this here by changing the semantics of segments to represent discrete physical space. The underlying Fenwick Tree accumulates continuous recombination mass from which we select recombination breakpoints assuming a uniform recombination rate. These breakpoints are then floored since segments have integer ends, but can easily be changed if segments move to continuous physical space.

This doesn't really add any more work so we wouldn't expect the performance to change but as a quick sanity check here's a rough benchmark of recombination and the overall runtime.

Simulation was ran with a few trials of the following:
msprime.simulate(10, Ne=10000, random_seed=1, recombination_map=msprime.RecombinationMap.uniform_map(length=10 ** 8, num_loci=10 ** 8, rate=1e-8)

| | msp_recombination_event | msp_run |
-------------------------|---------------------------|-----------
| master | ~7.7s (20.2%) | ~38s |
| recomb-phys-coord | ~8.3s (21.4%) | ~38s |

Note that with this change, the number of loci must be equal to the sequence length. However, under the continuous model loci lose their utility so just working around this for now by setting the number of loci manually.

Next step: incorporate recombination maps and see if the performance takes a serious hit.

[daniel-goldstein]

I've gotten recombination maps working for discrete physical coordinates. Below are the benchmarks for running stdpopsim HomSap -c chr22 -g HapMapII_GRCh37 <n_samples> -s 1 on my PR and master.

The performance hit is acceptable for moving forward but I think we can get this down a decent amount by caching the recombination mass for segment endpoints when they're created rather than recomputing them when merging ancestors. For now just focusing on getting everything working with these minor changes so we can merge.

[petrelharp]

This is great, and ping me when I should have a serious look through the code. Just to check - the proposal is to switch the underlying simulation method for 1.0? (I support this!)

[andrewkern]

and abandon continuous coords? would it be possible to hold on to the legacy version and have discrete coords as an option?

[daniel-goldstein]

@andrewkern Continuous coords is one of the next steps after this. The plan is to be able to specify either continuous or discrete space, just both in physical coordinates.

[jeromekelleher]

Just to check - the proposal is to switch the underlying simulation method for 1.0? (I support this!)

Yes - we're using the major version bump to introduce some slight semantic differences plus potentially breaking some people's code (but not much, I think)

and abandon continuous coords? would it be possible to hold on to the legacy version and have discrete coords as an option?

As @daniel-goldstein says, the next step will be to introduce an option to simulate on continuous coordinates as well. We have to have this, or we break way too much code.

[andrewkern]

sounds good to me! nice work @daniel-goldstein!!

[petrelharp]

... but, we can't do continuous coordinates, since floats are not continuous!

What you say makes sense, though. I'm just saying that how we describe the 'continuous' option could talk about it as being the "least discrete"?

[petrelharp]

Thinking about this more - could you give a quick run-down of what the proposal is? How are coordinates represented under the hood (int? float?), and on what scale are recombination breakpoints chosen? I've had a look at the code but it'll make a lot more sense if I know what the goal is.

[andrewkern]

one more question I'd like to add to the list-- why the performance hit? shouldn't this be the same complexity?

[jeromekelleher]

... but, we can't do continuous coordinates, since floats are not continuous!
What you say makes sense, though. I'm just saying that how we describe the 'continuous' option could talk about it as being the "least discrete"?

I assumed you were joking here @petrelharp --- sure, we all know that IEEE doubles aren't actually continuous, but for all practical purposes that's what a continuous value on a computer is? Wouldn't this be a very confusing way to talk about double coordinates?

[jeromekelleher]

Thinking about this more - could you give a quick run-down of what the proposal is? How are coordinates represented under the hood (int? float?), and on what scale are recombination breakpoints chosen? I've had a look at the code but it'll make a lot more sense if I know what the goal is.

So the proposal is that the left and right coordinates on the segment_t struct will now be doubles, and will always store physical coordinates. Breakpoints are only every generated as a result of recombination (or GC) events. When we generate a recombination breakpoint, we do this by choosing uniformly from the total recombination mass associated with all segments, and by translating this back into a physical coordinate. This will be a double value, and is naturally continuous. Then, if the is_discrete flag in the recombination map is set, we round this value to the closest integer. Since doubles can represent all integers up to 2^53 exactly, this allows us to work with any conceivable discrete chromosome.

There's a bit of trickiness involved when we're computing what the total recombination mass is, in discrete and continuous space, but we'll figure this out.

Make sense?

[jeromekelleher]

A part I forgot to mention is that the num_loci field in the recombination map will be removed. How to best manage the code breakage is TBD, but I think it's a pretty niche feature and only ever used current when trying to make the simulation work in discrete space.

[jeromekelleher]

one more question I'd like to add to the list-- why the performance hit? shouldn't this be the same complexity?

The perf hit comes because we're now computing the recombination mass associated with each segment, as those segments are being manipulated. In the existing code, we're calling genetic_to_phys on coordinates as they are about to be stored in the edge table. Now we're interacting with the recombination map a lot more.

As @daniel-goldstein says though, we're pretty sure we can do better than the current code. Also, we should still be much faster than any released version of msprime on a recombination map, since master has @ivan-krukov's binary search code that hasn't been released.

Actually, it's probably worth doing this experiment: would you mind adding a line to the plot above for the current release of msprime, for comparison @daniel-goldstein?

[daniel-goldstein]

Updated the plot with the current release installed through stdpopsim.

@andrewkern Like @jeromekelleher said, to minimize code changes on the first step we're doing a lot of conversion between recombination mass and physical position. We think we can bring that down to just one search per breakpoint when we sample them. Since we don't have to translate from genetic to physical coordinates anymore hopefully there will be a negligible performance hit in the end.

[jeromekelleher]

Updated the plot with the current release installed through stdpopsim.

Superb, thanks @daniel-goldstein!

[andrewkern]

thanks for the explanation guys. the new figure looks excellent @daniel-goldstein.

[petrelharp]

Makes sense. I like it.

(... I thought of several other complications, but talked myself out of them. =)

[jeromekelleher]

Closed in #862

[jeromekelleher]

Perf discussion moved to #873