finite `num_loci` argument produces breakpoints not lying between recombination map positions

[petrelharp]

If I wanted to simulate a region of five bases and equal recombination fractions, I could do this:

r = msprime.RecombinationMap(positions=[0,1,2,3,4],rates=[1.0, 1.0, 1.0, 1.0, 0.0],num_loci=4)
ts = msprime.simulate(100, recombination_map=r)
print(list(ts.breakpoints()))
# [0, 1.0, 2.0, 3.0, 4.0]

and I'd get recombination breakpoints occurring at the positions passed to the RecombinationMap, as expected. However, if recombination rates are not equal in the four regions, unexpected things happen, for instance:

r = msprime.RecombinationMap(positions=[0,1,2,3,4],rates=[100.0, 0.0, 1.0, 1.0, 0.0],num_loci=4)
ts = msprime.simulate(100, recombination_map=r)
print(list(ts.breakpoints()))
# [0, 0.255, 0.51, 0.765, 4.0]

This is strange because (a) all the breakpoints actually occur in the first interval, and (b) there are three non-endpoint breakpoints, just as above, even though one of the interval's rates was 0.0, meaning it should have acted like it had only 3 loci.

I suppose that num_loci works by dividing up the chromosome into segments with equal recombination rates. This means that you can't use this feature to simulate discrete recombination with unequal rates, I think.

[petrelharp]

I will clarify this in the docs in #548.

[petrelharp]

So maybe its NOTABUG - it's not clear what else we would do, because what if num_loci was 7 in the above example, what would we want to happen? - but I'm also not convinced it's all that useful, as-is, either.

[jeromekelleher]

I suppose that num_loci works by dividing up the chromosome into segments with equal recombination rates. This means that you can't use this feature to simulate discrete recombination with unequal rates, I think.

Yes, we calculate the overall recombination rate from the map, and then use this to compute the recombination rate between num_loci discrete chunks. We'd normally imagine num_loci as being very large (~2**32). For the recombination algorithm to work we must have a uniform recombination rate among k loci, otherwise the way of efficiently finding recombination breakpoints doesn't work. We work around this limitation then by rescaling between genetic and physical distances at the end.

It would take a pretty fundamental rethink to do anything else I'd say...

[petrelharp]

It would take a pretty fundamental rethink to do anything else I'd say...

If you think my addition to the documentation has clarified this, then go ahead and close.

[jeromekelleher]

Closed in #548.