Repository: tskit-dev/msprime (language: python)

Issue:
Segfault with `BetaCoalescent` and ploidy greater than diploid
The following snippet gives me a segfault:

```python
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),
)
```


Discussion / hints:
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


Confirmed also on Linux, v 1.3.1. Investigating.
Well the problem is a pretty straighforward. We assume there are at most four "pots" to assign lineages into here: https://github.com/tskit-dev/msprime/blob/657fc8dee8f9753b657efad71a06fb5cacfdaf32/lib/msprime.c#L8013

But a few lines later we set ``num_parental_copies`` to 2 * ploidy, which we then use the number of "pots": 
https://github.com/tskit-dev/msprime/blob/657fc8dee8f9753b657efad71a06fb5cacfdaf32/lib/msprime.c#L8025

@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.

Fix patch:
diff --git a/lib/msprime.c b/lib/msprime.c
index 35560d095..b7b2b30ea 100644
--- a/lib/msprime.c
+++ b/lib/msprime.c
@@ -8010,7 +8010,8 @@ msp_dirac_common_ancestor_event(msp_t *self, population_id_t pop_id, label_id_t
 {
     int ret = 0;
     uint32_t j, n, num_participants, num_parental_copies;
-    avl_tree_t *ancestors, Q[4]; /* MSVC won't let us use num_pots here */
+    avl_tree_t *ancestors;
+    avl_tree_t *Q = NULL;
     avl_node_t *x_node, *y_node;
     segment_t *x, *y;
     lineage_t *x_lin, *y_lin;
@@ -8024,6 +8025,11 @@ msp_dirac_common_ancestor_event(msp_t *self, population_id_t pop_id, label_id_t
     } else {
         num_parental_copies = 2 * self->ploidy;
     }
+    Q = tsk_malloc(num_parental_copies * sizeof(*Q));
+    if (Q == NULL) {
+        ret = MSP_ERR_NO_MEMORY;
+        goto out;
+    }
 
     ancestors = &self->populations[pop_id].ancestors[label];
     n = avl_count(ancestors);
@@ -8080,6 +8086,7 @@ msp_dirac_common_ancestor_event(msp_t *self, population_id_t pop_id, label_id_t
         }
     }
 out:
+    tsk_safe_free(Q);
     return ret;
 }
 
@@ -8229,7 +8236,8 @@ msp_beta_common_ancestor_event(msp_t *self, population_id_t pop_id, label_id_t l
 {
     int ret = 0;
     uint32_t j, n, num_participants, num_parental_copies;
-    avl_tree_t *ancestors, Q[4]; /* MSVC won't let us use num_pots here */
+    avl_tree_t *ancestors;
+    avl_tree_t *Q = NULL;
     double alpha = self->model.params.beta_coalescent.alpha;
     double truncation_point = beta_compute_truncation(self);
     double beta_x, u, increment;
@@ -8241,6 +8249,11 @@ msp_beta_common_ancestor_event(msp_t *self, population_id_t pop_id, label_id_t l
     } else {
         num_parental_copies = 2 * self->ploidy;
     }
+    Q = tsk_malloc(num_parental_copies * sizeof(*Q));
+    if (Q == NULL) {
+        ret = MSP_ERR_NO_MEMORY;
+        goto out;
+    }
 
     for (j = 0; j < num_parental_copies; j++) {
         avl_init_tree(&Q[j], cmp_segment_queue, NULL);
@@ -8300,6 +8313,7 @@ msp_beta_common_ancestor_event(msp_t *self, population_id_t pop_id, label_id_t l
     }
 
 out:
+    tsk_safe_free(Q);
     return ret;
 }
 


Test patch:
diff --git a/lib/tests/test_ancestry.c b/lib/tests/test_ancestry.c
index e165d0bc9..8425a520f 100644
--- a/lib/tests/test_ancestry.c
+++ b/lib/tests/test_ancestry.c
@@ -1515,6 +1515,47 @@ test_multiple_mergers_simulation(void)
     gsl_rng_free(rng);
 }
 
+static void
+test_multiple_mergers_ploidy(void)
+{
+    int ret;
+    size_t j, ploidy;
+    uint32_t n = 10;
+    long seed = 10;
+    msp_t msp;
+    gsl_rng *rng = safe_rng_alloc();
+    tsk_table_collection_t tables;
+
+    for (j = 0; j < 2; j++) {
+        gsl_rng_set(rng, seed);
+        for (ploidy = 1; ploidy < 12; ploidy++) {
+            ret = build_sim(&msp, &tables, rng, 10, 1, NULL, n);
+            CU_ASSERT_EQUAL(ret, 0);
+            CU_ASSERT_EQUAL_FATAL(msp_set_recombination_rate(&msp, 0.1), 0);
+            if (j == 0) {
+                ret = msp_set_simulation_model_dirac(&msp, 0.9, 10);
+            } else {
+                ret = msp_set_simulation_model_beta(&msp, 1.8, 1);
+            }
+            CU_ASSERT_EQUAL(ret, 0);
+
+            msp_set_ploidy(&msp, ploidy);
+            ret = msp_initialise(&msp);
+
+            ret = msp_run(&msp, DBL_MAX, ULONG_MAX);
+            CU_ASSERT_EQUAL_FATAL(ret, 0);
+            CU_ASSERT_TRUE(msp_is_completed(&msp));
+            CU_ASSERT_TRUE(msp.time > 0);
+            msp_verify(&msp, 0);
+
+            ret = msp_free(&msp);
+            CU_ASSERT_EQUAL(ret, 0);
+            tsk_table_collection_free(&tables);
+        }
+    }
+    gsl_rng_free(rng);
+}
+
 static void
 test_multiple_mergers_growth_rate(void)
 {
@@ -4305,6 +4346,7 @@ main(int argc, char **argv)
         { "test_gc_rates", test_gc_rates },
 
         { "test_multiple_mergers_simulation", test_multiple_mergers_simulation },
+        { "test_multiple_mergers_ploidy", test_multiple_mergers_ploidy },
         { "test_multiple_mergers_growth_rate", test_multiple_mergers_growth_rate },
         { "test_dirac_coalescent_bad_parameters", test_dirac_coalescent_bad_parameters },
         { "test_beta_coalescent_bad_parameters", test_beta_coalescent_bad_parameters },