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")
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 (withdemography.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: