In very, very, very rare cases msprime can generate a zero waiting time until the next coalescent event, resulting in an error being thrown. This script for example, does it:
import msprime as msp
def ancient_sample_test(
num_modern=1000, anc_pop = 0, anc_num = 1, anc_time = 200, split_time_anc = 400,
Ne0 = 10000, Ne1 = 10000, length = 1000):
samples = [msp.Sample(population = 0, time = 0)]*num_modern
samples.extend([msp.Sample(population = anc_pop, time = anc_time)]*(2*anc_num))
pop_config = [msp.PopulationConfiguration(initial_size = Ne0), msp.PopulationConfiguration(initial_size = Ne1)]
divergence = [msp.MassMigration(time = split_time_anc, source = 1, destination = 0, proportion = 1.0)]
seed = 94320219
sims = msp.simulate(
samples=samples,Ne=Ne0,population_configurations=pop_config,
demographic_events=divergence, length=length,
random_seed=seed)
if __name__ == "__main__":
num_ind = 10
ancient_sample_test(
num_modern=100,anc_pop=0,anc_num=num_ind,Ne0=3000,Ne1=3000,anc_time=5419,split_time_anc=5919,length=500)
we get
File "ancient_genotypes_simulation.py", line 19, in <module>
num_modern=100,anc_pop=0,anc_num=num_ind,Ne0=3000,Ne1=3000,anc_time=5419,split_time_anc=5919,length=500)
File "ancient_genotypes_simulation.py", line 14, in ancient_sample_test
random_seed=seed)
File "/home/jk/work/github/msprime/msprime/simulations.py", line 485, in simulate
sim, mutation_generator, 1, provenance_dict, end_time))
File "/home/jk/work/github/msprime/msprime/simulations.py", line 163, in _replicate_generator
sim.run(end_time)
File "/home/jk/work/github/msprime/msprime/simulations.py", line 704, in run
self.ll_sim.run(end_time)
_msprime.LibraryError: The simulation model supplied resulted in a parent node having a time value <= to its child. This can occur either as a result of multiple bottlenecks happening at the same time or because of numerical imprecision with very small population sizes.
What happens is, this line returns exactly zero and so the minimum time until the next event is 0 and naturally enough this gets chosen as the next event. But, we end up with a zero branch length, which tskit won't allow.
We could put in some cludge where we have 1e-200 or something as the time if t_wait
is exactly zero, but I don't think this is a good idea. We can easily get to a case where branch lengths become too short for double precision when growth rates are specified, and it's much better that the simulation ends with an error than the user getting back a bunch of completed trees with miniscule branch lengths.
So; I don't really know how to handle this!
Thanks to David Lawrie for the script and reporting the problem. Pinging @molpopgen for thoughts.