Msprime: Fix up time scaling in multiple merger models

[jeromekelleher]

1028 removed the concept of a "reference size" for simulation models entirely, as it turns out that it did nothing. That PR removed all the infrastructure for moving back and forth between time scales. However, the multiple merger coalescents are not quite right afterwards, and some of the code will need to be ported back. See https://github.com/tskit-dev/msprime/pull/1028#issuecomment-637019349

Must also check the run_xi_dirac_recombinations test to reinstate the division by Ne in the mutation rate according to @JereKoskela

I do struggle to understand why if the tests we have are passing without doing any of this, though!

cc @JereKoskela, @TPPSellinger

[TPPSellinger]

I will tomorrow investigate a bit deeper all the test we have and check some other values which could indicate an issue.

[TPPSellinger]

I went through all the tests and you are right they all seem to work (which is good). The reason they all work it's because none of them check the for the population size. We forgot to check if distribution were correct given the population size. This might be a tricky part with multiple merger as @JereKoskela explained. At first we should simply assume constant population size and check if we have the correct scaling. We could check if branch length fits expected value (given population size). If yes, we are good for the constant population case. If not this is quite problematic. Once constant population size works we could check for non constant population size (e.g. expansion) as @JereKoskela suggested.

I will build the test for constant population size. I hope to have results tomorrow or thursday.

[jeromekelleher]

This sounds like an excellent plan, thanks @TPPSellinger!

One thing though - we have a large refactor of the verification.py code underway in #1035, and it would be a good idea to wait until that's merged before starting (or else you'll have some horrible conflicts to resolve). Should be merged in the next hour.

[TPPSellinger]

Okay I will wait.

[JereKoskela]

I do struggle to understand why if the tests we have are passing _without_ doing any of this, though!

I'll try to put together an explanation. Let me know if it clarifies things at all @jeromekelleher.

The rates of common ancestor events (n choose 2 for binary mergers and 0 for everything else for Hudson, say) is determined by the infinite-population coalescent process. The relationship between the population size and the waiting times for common ancestor events is not determined by the coalescent limit, but by the pre-limiting sequence of finite-population models who's ancestries converge to the coalescent.

For example, the Hudson model can be obtained as a limit of Wright-Fisher models, in which case waiting times are proportional to population sizes, exactly as implemented. But you could choose to obtain it as a limit of Moran models instead, in which case waiting times would need to be proportional to the square of the population size. No test based upon normalised branch lengths can detect which choice you have made, because both result in a correct coalescent. Checking that densities of recombination and mutation events are proportional to branch lengths will also not work for the same reason. In fact, you could probably plug in pretty much any population size relationship you like for common ancestor waiting times, and the only thing that would need to change is what you mean by "per-generation" when specifying recombination/mutation/etc rates.

It's the same situation for the multiple merger coalescents. To take the Beta-coalescent as an example, the most natural derivation available (in my opinion) is the Schweinsberg model, for which common ancestor waiting times are proportional to pop_size^(alpha - 1). We would still be simulating the correct coalescent even if we use a different scaling, but would end up with a lot of users obtaining output with a recombination rate that is different to what they thought they specified.

[TPPSellinger]

Couldn't we just chose and specify the scaling we are using (e.g. the one in Schweinsberg model ) ? Then check that given the scaling we have chosen, the obtained branches make sense ?

[JereKoskela]

Couldn't we just chose and specify the scaling we are using (e.g. the one in Schweinsberg model ) ? Then check that given the scaling we have chosen, the obtained branches make sense ?

Yes, provided we make precise what it means for branches to "make sense". Results on quantities like unnormalised branch lengths are very sparse for multiple merger models.

In this newest iteration of the code, every model uses the standard Wright-Fisher scaling. I'm intending to go over the code to change it back for multiple merger models asap.

[TPPSellinger]

I see what you mean. Why not first work with theta, at least for constant population size? I mean if the SFS is good and the theta distribution is also good (given model and model parameters), it would be a first step in the right direction.

PS: should I wait your are done with the code before doing anything ?

[JereKoskela]

PS: should I wait your are done with the code before doing anything ?

Yeah, that's probably wise. I hope to get things done this week, though that may end up slipping a bit. There's not much point testing an implementation we know is not what we want though. I'll keep you informed.

[jeromekelleher]

Thanks @JereKoskela, that all sounds good to me. Let's wait on Jere's updates before adding any further verification tests.

[jeromekelleher]

Closed in #1059