Msprime: factor of two in `instantaneous_bottleneck`?

[petrelharp]

For unrelated reasons, @andrewkern and I were looking at the code, and were confused about this calculation:

/* Instantaneous bottlenecks. At a time T1 we have a burst of coalescence
 * equivalent to what would happen in time T2.
 */

msp_instantaneous_bottleneck(msp_t *self, demographic_event_t *event)
   /// ... snip ...
    /* Now we implement the Kingman coalescent for these lineages until we have
     * exceeded T2. This is based on the algorithm from Hudson 1990.
     */
    j = n - 1;
    t = 0.0;
    parent = (tsk_id_t) n;
    while (j > 0) {
        rate = j + 1;
        rate = rate * j;
        t += msp_get_common_ancestor_waiting_time_from_rate(
            self, &self->populations[population_id], rate);
        if (t >= T2) {
            break;
        }

Compare this to what happens in the standard coalescent, which is here:

msp_std_get_common_ancestor_waiting_time(
    msp_t *self, population_id_t pop_id, label_id_t label)
{
    population_t *pop = &self->populations[pop_id];
    double n = (double) avl_count(&pop->ancestors[label]);
    double lambda = n * (n - 1.0) / 2.0;

    return msp_get_common_ancestor_waiting_time_from_rate(self, pop, lambda);

Notice that the second bit of code passes in n (n-1)/2 for the rate (lambda); while the first bit of code passes in n (n-1), despite saying that it's just doing Kingman for time T2. This looks suspicious?

[jeromekelleher]

This is quite likely wrong - Konrad and I have been talking about getting this tested and documented for about 4 years now, but never actually gotten to it.

So, suspicion is well founded, I think.

[jeromekelleher]

See also #9 (note the small issue number!)

[JereKoskela]

Dividing the rate by two in msp_std_get_common_ancestor_waiting_time is a change I made some time ago, and it should be there when simulating the standard Kingman coalescent.

The aim is to simulate mergers at rate n (n - 1) / 2, or 1 per pair of lineages, and recombinations at rate ρ / 2 per link (of which there are k, say). That can be done by simulating waiting times with rates n(n - 1) and k ρ, and dividing the resulting minimum by two to correct for doubling the rates. ms did something like that, which is why it was also the way msprime was originally implemented, as far as I understand.

When we got rid of time scaling, it also didn't make sense to use the population-rescaled recombination rate ρ / 2 anymore. Instead, the recombination rate is just r per link per generation. So the factors of 1/2 don't cancel in the rates anymore, which is why it needs to go into the common ancestor rate explicitly.

[petrelharp]

The factor of 2 in n (n-1) / 2 is definitely right. Good work with that, everything is so much easier to think about.