Segfault with BetaCoalescent and ploidy greater than diploid

[nspope]

The following snippet gives me a segfault:

import msprime
print(msprime.__version__) # 1.3.2
msprime.sim_ancestry(
    samples=10,
    recombination_rate=1e-8,
    sequence_length=1e6,
    population_size=1e4,
    ploidy=3,
    random_seed=1024,
    model=msprime.BetaCoalescent(alpha=1.5),
)

[andrewkern]

fwiw I can recreate a similar error on msprime 1.3.1 on my mac laptop, although it's telling me abort rather than segfault

[jeromekelleher]

Confirmed also on Linux, v 1.3.1. Investigating.

[jeromekelleher]

Well the problem is a pretty straighforward. We assume there are at most four "pots" to assign lineages into here:

msprime/lib/msprime.c

Line 8013 in 657fc8d

avl_tree_t *ancestors, Q[4]; /* MSVC won't let us use num_pots here */

But a few lines later we set num_parental_copies to 2 * ploidy, which we then use the number of "pots":

msprime/lib/msprime.c

Line 8025 in 657fc8d

num_parental_copies = 2 * self->ploidy;

@JereKoskela can you comment here please? Is the correct thing to do here just to malloc Q with num_parental_copies and not just 4 (easy). We would need to do some statistical validation as well I think, as this has clearly never been tested for e.g. triploids.

[jeromekelleher]

#2308 fixes the segfault. I'd need some help in figuring out whether it's actually doing the right thing, though.

[JereKoskela]

Making num_parental_copies the number of buckets is the right thing to do.

Statistical validation is not totally straightforward because this is the only implementation for ploidy > 2 of which I'm aware. We could estimate moments of TMRCAs of n = 2, for which the expectation and variance are easy to compute for various ploidies, say 1, ..., 6.

[jeromekelleher]

Statistical validation is not totally straightforward because this is the only implementation for ploidy > 2 of which I'm aware. We could estimate moments of TMRCAs of n = 2, for which the expectation and variance are easy to compute for various ploidies, say 1, ..., 6.

That sounds great - would you be able to pick this up if I merge #2308?

I'll push out a bugfix release in the next couple of days, as we don't really want segfault bugs lying around.

[JereKoskela]

Sure, I'll put something together. Go ahead and merge.

[jeromekelleher]

Great, thanks a lot. I'll push out the bugfix once we're reasonably sure that the fix works.

[jeromekelleher]

Closed in #2310 and #2308. Great work everyone!