demography gotcha with two independent populations

[JoseBlanca]

Hi,

It is my first day with msprime. I'm trying to prepare some simulations for a population genetics class.
I've read the documentation and I've prepared my first simulation. With one population everything has gone well. Then I was thinking about running two independent populations with different sizes to compare the results. The simulation won't finished in that case. After some thought I've realized that the problem (I believe) is that if the populations are completely independent they trees will never coalesce and the simulation will run forever, msprime will stalled at 100% CPU for ever. If this is really the case maybe a check when the demography is first used or a warning in the documentation could be useful. Maybe I'm wrong, I'm not completely sure.
I'm attaching a code example with the problem.
Thanks for your great work. Best,

José Blanca

import time

import msprime

recomb_rate = 1e-8
pop_size = 100
num_indis_to_sample = 5
ploidy = 2
seq_length_in_bp = 100

do_it_wrong = True

demography = msprime.Demography()
demography.add_population(name="pop1", initial_size=pop_size)

samples_pop1 = msprime.SampleSet(
    num_samples=num_indis_to_sample, population="pop1", ploidy=ploidy
)
samples = [samples_pop1]
time1 = time.time()
tree_seqs = msprime.sim_ancestry(
    samples=samples,
    demography=demography,
    recombination_rate=recomb_rate,
    sequence_length=seq_length_in_bp,
)
time2 = time.time()
time_ellapsed = time2 - time1
print(f"{time_ellapsed=}")

if not do_it_wrong:
    samples_pop0 = msprime.SampleSet(
        num_samples=num_indis_to_sample, population="pop0", ploidy=ploidy
    )
samples_pop2 = msprime.SampleSet(
    num_samples=num_indis_to_sample, population="pop2", ploidy=ploidy
)
samples_pop3 = msprime.SampleSet(
    num_samples=num_indis_to_sample, population="pop3", ploidy=ploidy
)
demography = msprime.Demography()
if not do_it_wrong:
    demography.add_population(name="pop0", initial_size=pop_size)
demography.add_population(name="pop2", initial_size=pop_size)
demography.add_population(name="pop3", initial_size=pop_size)
if not do_it_wrong:
    demography.add_population_split(100, derived=["pop2", "pop3"], ancestral="pop0")
samples = [samples_pop2, samples_pop3]
time3 = time.time()
tree_seqs = msprime.sim_ancestry(
    samples=samples,
    demography=demography,
    recombination_rate=recomb_rate,
    sequence_length=seq_length_in_bp,
)
time4 = time.time()
time_ellapsed = time4 - time3
print(f"{time_ellapsed=}")

[JoseBlanca]

I have just seen that there is an infinite-waintg-time LibraryError, but I did not get it.

https://tskit.dev/msprime/docs/stable/ancestry.html#infinite-waiting-time
Best,

José Blanca

[petrelharp]

Good suggestion. This looks related to #649, since running this simplified code:

import msprime

num_indis_to_sample = 5
pop_size = 10
recomb_rate = 0.0

samples_pop2 = msprime.SampleSet(
    num_samples=num_indis_to_sample, population="pop2",
)
samples_pop3 = msprime.SampleSet(
    num_samples=num_indis_to_sample, population="pop3",
)
demography = msprime.Demography()
demography.add_population(name="pop2", initial_size=pop_size)
demography.add_population(name="pop3", initial_size=pop_size)
samples = [samples_pop2, samples_pop3]

tree_seqs = msprime.sim_ancestry(
    samples=samples,
    demography=demography,
    recombination_rate=recomb_rate,
    sequence_length=10,
)

gets

msprime._msprime.LibraryError: Infinite waiting time until next simulation event.

but setting the recombination rate to something nonzero does not get the error. I think we should be able to catch this problem - it just takes some thinking about how to do it.

We also have a stub for discussion of this in the documentation (issue #1527).

[jeromekelleher]

Hi @JoseBlanca, welcome to msprime 👋

This is a common gotcha, and we have a section in the docs for it. We haven't gotten around to filling out the details though. If you'd like to help out, then getting some more details into that section would be very helpful.

Basically it's tricky 😄