Mass migration events at time zero are skipped

[nspope]

Ran across this bug (?) while writing some tests-- consider a two-population island model, both populations of size $N$, with a mass migration event at time 0 where 50% of lineages from population B move to population A. Given this model, the expected proportion of (B, B) pairs that haven't coalesced by time $t$ is $0.5 + 0.5 \exp( -\frac{1}{2N} t )$. However, the demography debugger and simulations don't follow this expectation unless I set the time of the mass migration to a very small floating point value:

Note that in the first panel (with demography.add_mass_migration(time=0.0, ...) the demography debugger and simulator match neither one another nor the theoretical curve. Whereas, in the second panel (with demography.add_mass_migration(time=1e-15, ...) everything lines up. This is a bit of an edge case, as the mass migration will effectively just randomly reassign the population of sample nodes -- but silently ignoring it is probably not ideal (and I'm not sure what's going on with the debugger).

Code to reproduce:

import msprime
import numpy as np
import matplotlib.pyplot as plt

print(msprime.__version__) # '1.3.3.dev29+ge90c338' 

fig, axs = plt.subplots(
    1, 2, figsize=(8, 4), sharey=True, sharex=True, 
    constrained_layout=True,
)
for i, admix_time in enumerate([0.0, 1e-15]):
    # isolation model with "admixture event" at time 0 or slightly later
    demog = msprime.Demography()
    demog.add_population(name="A", initial_size=1e4)
    demog.add_population(name="B", initial_size=1e4)
    demog.add_population(name="C", initial_size=1e4)
    demog.add_mass_migration(time=admix_time, source="B", dest="A", proportion=0.5)
    demog.add_population_split(time=1e4, ancestral="C", derived=["A", "B"])
    time_grid = np.linspace(0, 1e4, 7)

    # surviving (B, B) pairs from debugger
    dbg = demog.debug()
    _, surv_debug = dbg.coalescence_rate_trajectory(lineages={"B": 2}, steps=time_grid)

    # surviving (B, B) pairs from simulation
    num_reps = 10000
    ts_gen = msprime.sim_ancestry(
        samples={"B":1}, 
        demography=demog, 
        recombination_rate=1e-8, 
        sequence_length=1e6, 
        num_replicates=num_reps,
    )
    E_cdf = np.zeros(time_grid.size)
    for ts in ts_gen:
        E_cdf[1:] += np.cumsum(
            ts.pair_coalescence_counts(
                time_windows=time_grid, 
                pair_normalise=True,
            )
        )
    surv_sim = 1 - E_cdf / num_reps

    # surviving (B, B) pairs from theory
    surv_theory = 0.5 * (1 + np.exp(-1 / 1e4 / 2 * time_grid))

    axs[i].plot(time_grid, surv_debug, "-o", color="red", label="debugger")
    axs[i].plot(time_grid, surv_sim, "-o", color="blue", label="simulation")
    axs[i].plot(time_grid, surv_theory, "-o", color="black", label="theory")
    axs[i].legend()
    axs[i].set_title(f"MassMig at time {admix_time}")
fig.supylabel("Proportion uncoalesced (B, B) pairs")
fig.supxlabel("Time ago")
plt.savefig("msp_skipped_massmig_bug.png")

[jeromekelleher]

Nice bug report, thanks! What's the output of print(demog.debug())?

[molpopgen]

I think the issue is a bit more complicated -- this model would have been rejected if it were given as a demes.Graph. I was just discussing this with @apragsdale and I discovered that the demes spec doesn't match the actual behavior. The spec implies that the model would be accepted but it wouldn't be and there are tests of that.

[nspope]

@jeromekelleher -- here's the debugger:

DemographyDebugger
╠══════════════════════════════════╗
║ Epoch[0]: [0, 1e+04) generations ║
╠══════════════════════════════════╝
╟    Populations (total=3 active=2)
║    ┌───────────────────────────────────────────────┐
║    │   │     start│       end│growth_rate  │ A │ B │
║    ├───────────────────────────────────────────────┤
║    │  A│   10000.0│   10000.0│ 0           │ 0 │ 0 │
║    │  B│   10000.0│   10000.0│ 0           │ 0 │ 0 │
║    └───────────────────────────────────────────────┘
╟    Events @ generation 1e+04
║    ┌───────────────────────────────────────────────────────────────────────────────┐
║    │   time│type        │parameters       │effect                                  │
║    ├───────────────────────────────────────────────────────────────────────────────┤
║    │  1e+04│Population  │derived=[A, B],  │Moves all lineages from derived         │
║    │       │Split       │ancestral=C      │populations 'A' and 'B' to the          │
║    │       │            │                 │ancestral 'C' population. Also set the  │
║    │       │            │                 │derived populations to inactive, and    │
║    │       │            │                 │all migration rates to and from the     │
║    │       │            │                 │derived populations to zero.            │
║    └───────────────────────────────────────────────────────────────────────────────┘
╠════════════════════════════════════╗
║ Epoch[1]: [1e+04, inf) generations ║
╠════════════════════════════════════╝
╟    Populations (total=3 active=1)
║    ┌───────────────────────────────────────┐
║    │   │     start│       end│growth_rate  │
║    ├───────────────────────────────────────┤
║    │  C│   10000.0│   10000.0│ 0           │
║    └───────────────────────────────────────┘

and the demography itself:

Demography                                                                                                                                                                                     
╟  Populations                                                                                                                                                                                 
║  ┌───────────────────────────────────────────────────────────────────────────────────────────────┐                                                                                           
║  │ id │name  │description  │initial_size  │ growth_rate │  default_sampling_time│extra_metadata  │                                                                                           
║  ├───────────────────────────────────────────────────────────────────────────────────────────────┤                                                                                           
║  │ 0  │A     │             │10000.0       │      0      │                      0│{}              │                                                                                           
║  │ 1  │B     │             │10000.0       │      0      │                      0│{}              │                                     
║  │ 2  │C     │             │10000.0       │      0      │                  1e+04│{}              │                                     
║  └───────────────────────────────────────────────────────────────────────────────────────────────┘                                     
╟  Migration Matrix                                                                                                                      
║  ┌───────────────┐                                                                                                                     
║  │   │ A │ B │ C │                                                                                                                     
║  ├───────────────┤                                                                                                                     
║  │  A│ 0 │ 0 │ 0 │                                                                                                                     
║  │  B│ 0 │ 0 │ 0 │                                                                                                                     
║  │  C│ 0 │ 0 │ 0 │                                                                                                                     
║  └───────────────┘                                                                                                                     
╟  Events                                                                                                                                                      
║  ┌─────────────────────────────────────────────────────────────────────────────────────┐                                                                     
║  │   time│type            │parameters         │effect                                  │                                               
║  ├─────────────────────────────────────────────────────────────────────────────────────┤                                               
║  │      0│Mass Migration  │source=B, dest=A,  │Lineages currently in population B      │                                               
║  │       │                │proportion=0.5     │move to A with probability 0.5          │          
║  │       │                │                   │(equivalent to individuals migrating    │          
║  │       │                │                   │from A to B forwards in time)           │          
║  │┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈│          
║  │  1e+04│Population      │derived=[A, B],    │Moves all lineages from derived         │                                               
║  │       │Split           │ancestral=C        │populations 'A' and 'B' to the          │                                               
║  │       │                │                   │ancestral 'C' population. Also set the  │                    
║  │       │                │                   │derived populations to inactive, and    │                                               
║  │       │                │                   │all migration rates to and from the     │                                               
║  │       │                │                   │derived populations to zero.            │                                               
║  └─────────────────────────────────────────────────────────────────────────────────────┘   

[nspope]

@molpopgen you're right, demography.to_demes() errors out with

ValueError: invalid pulse at time=0.0, which is dest=B's end_time

So maybe the thing to do is error out in msprime as well, and update the deme spec?

[molpopgen]

update the deme spec?

I'm prepping something for later this week.

[jeromekelleher]

Hmm, so it looks like those mass migrations at time 0 aren't happening at all in msprime, they're not picked up by the Debugger. I'm not surprised that the coalescence_rate_trajectory gives a different answer here as it has it's own logic for reasoning about the effects of events, but the output of print(demon.debug()) is a direct view on what the simulation engine does.

We could probably fix both these things, but maybe it's a bit of an odd corner case and we should just raise an error?

[nspope]

I think raising an error is a good move, especially if this is consistent with demes.

[molpopgen]

We could probably fix both these things, but maybe it's a bit of an odd corner case and we should just raise an error?

One simplistic way to raise an error would be to make an internal conversion of the model to a demes Graph (and avoid doing so if the input is such a Graph b/c that's presumably okay).

[petrelharp]

Ideally the error would suggest the preferred way to do this. (Is that: create a new admixed population at time zero?)

[nspope]

(Is that: create a new admixed population at time zero?)

I think introduce the mass migration at a infinitesimal time step is the easier option -- otherwise, one would have to make samples haploid, "migrate" them manually with an RNG, pass the new population-sample-sizes to msprime, and then dump tables and recode populations back to the original assignments. (Unless I'm missing some other shortcut?)

[apragsdale]

I don't think there is a way to specify such a model using the demes framework. Pulses into a deme at its end time are generally disallowed, which might be overly restrictive for continuous time simulators, but makes more sense in discrete simulators. In any case, it's what the spec says.

Ideally the error would suggest the preferred way to do this. (Is that: create a new admixed population at time zero?)

This would also raise an error in demes, given that the admixed population would have its start time and end time both equal to zero, and demes need to have a positive total existence time. I'm not sure how specifying a model using an admixture event directly with msprime model specification would behave, and if those events at time zero would be applied correctly.

[apragsdale]

Admixture events at time zero are not picked up either.

import msprime

demog = msprime.Demography()
demog.add_population(name="A", initial_size=1e4)
demog.add_population(name="B", initial_size=1e4)
demog.add_population(name="C", initial_size=1e4)
demog.add_population(name="D", initial_size=1e4)

demog.add_admixture(0, derived="A", ancestral=["B", "C"], proportions=[0.5, 0.5])
demog.add_population_split(time=1e4, ancestral="D", derived=["B", "C"])

dbg = demog.debug()

In [13]: print(dbg)
DemographyDebugger
╠══════════════════════════════════╗
║ Epoch[0]: [0, 1e+04) generations ║
╠══════════════════════════════════╝
╟    Populations (total=4 active=2)
║    ┌───────────────────────────────────────────────┐
║    │   │     start│       end│growth_rate  │ B │ C │
║    ├───────────────────────────────────────────────┤
║    │  B│   10000.0│   10000.0│ 0           │ 0 │ 0 │
║    │  C│   10000.0│   10000.0│ 0           │ 0 │ 0 │
║    └───────────────────────────────────────────────┘
╟    Events @ generation 1e+04
║    ┌───────────────────────────────────────────────────────────────────────────────┐
║    │   time│type        │parameters       │effect                                  │
║    ├───────────────────────────────────────────────────────────────────────────────┤
║    │  1e+04│Population  │derived=[B, C],  │Moves all lineages from derived         │
║    │       │Split       │ancestral=D      │populations 'B' and 'C' to the          │
║    │       │            │                 │ancestral 'D' population. Also set the  │
║    │       │            │                 │derived populations to inactive, and    │
║    │       │            │                 │all migration rates to and from the     │
║    │       │            │                 │derived populations to zero.            │
║    └───────────────────────────────────────────────────────────────────────────────┘
╠════════════════════════════════════╗
║ Epoch[1]: [1e+04, inf) generations ║
╠════════════════════════════════════╝
╟    Populations (total=4 active=1)
║    ┌───────────────────────────────────────┐
║    │   │     start│       end│growth_rate  │
║    ├───────────────────────────────────────┤
║    │  D│   10000.0│   10000.0│ 0           │
║    └───────────────────────────────────────┘

[petrelharp]

It makes sense that a continuous-time event at zero should have no effect (since the half-open interval [x,x), if defined, is empty), but should raise an error. I suspect that's the source of some of this inconsistency, since stuff was coded thinking about continuous-time events?

But an instantaneous event does happen at a single time; however, is the issue that there is ambiguity around which of various instantaneous events happens first? (In this case, "mass migration" and "t=0 sampling" are the two instantaneous events).

This does feel like a bug/feature request - I don't see why you shouldn't be able to do a mass migration at t=0. (Anything that produced a change in node would not work, since you can't have parent and child at the same time.) But, maybe I'm missing something.

[nspope]

This does feel like a bug/feature request - I don't see why you shouldn't be able to do a mass migration at t=0. (Anything that produced a change in node would not work, since you can't have parent and child at the same time.) But, maybe I'm missing something.

Hmm -- if record_full_arg=True, then we'll end up with non-sample nodes at the time of the mass migration. If these coincide with sample nodes, then presumably we'd have issues as branch lengths would be zero? Unless simultaneous events are already collapsed into the same node by the simulation code?

For example, this errors out:

import msprime
import numpy as np

demog = msprime.Demography()
demog.add_population(name="A", initial_size=1e4)
demog.add_population(name="B", initial_size=1e4)
demog.add_population(name="C", initial_size=1e4)
demog.add_mass_migration(1e-15, source="A", dest="B", proportion=0.5)
demog.add_population_split(time=1e4, ancestral="C", derived=["A", "B"])

ts = msprime.sim_ancestry(
    samples={"A":10}, demography=demog, recombination_rate=1e-8, 
    sequence_length=1e6, record_full_arg=True, random_seed=1,
)

assert np.min(ts.nodes_time[ts.nodes_time > 0]) == 1e-15

tab = ts.dump_tables()
new_nodes_time = tab.nodes.time.copy()
new_nodes_time[new_nodes_time == 1e-15] = 0.0
tab.nodes.time = new_nodes_time
tab.sort()  # library error