msp_gsl_ran_flat()

msp_gsl_ran_flat()

incorporated comments from Jerome and Kevin

incorporated additional edits from Jerome

CHANGELOG.md

@@ -1,5 +1,12 @@
 # Changelog
 
+## [1.1.1] - 2022-02-xx
+
+**Bug fixes**:
+- Fix (very) rare assertion trip caused by underlying GSL bug.
+  ({issue}`1997`, {pr}`2000`, {user}`chriscrsmith `, {user}`molpopgen`,
+  {user}`andrewkern`)
+
 ## [1.1.0] - 2021-12-14
 
 **New features**

lib/mutgen.c

@@ -202,7 +202,7 @@ mutation_matrix_choose_root_state(mutation_model_t *self, gsl_rng *rng, site_t *
 {
     int ret = 0;
     mutation_matrix_t params = self->params.mutation_matrix;
-    double u = gsl_ran_flat(rng, 0.0, 1.0);
+    double u = msp_gsl_ran_flat(rng, 0.0, 1.0);
     size_t j = probability_list_select(u, params.num_alleles, params.root_distribution);
     tsk_bug_assert(j < params.num_alleles);
     site->ancestral_state = params.alleles[j];
@@ -218,7 +218,7 @@ mutation_matrix_transition(mutation_model_t *self, gsl_rng *rng,
 {
     int ret = 0;
     mutation_matrix_t params = self->params.mutation_matrix;
-    double u = gsl_ran_flat(rng, 0.0, 1.0);
+    double u = msp_gsl_ran_flat(rng, 0.0, 1.0);
     double *probs;
     tsk_id_t j, pi;
 
@@ -1125,15 +1125,15 @@ mutgen_place_mutations(mutgen_t *self, bool discrete_sites)
                  * use up all of the doubles before it could happen and so we'd
                  * certainly run out of memory first. */
                 do {
-                    position = gsl_ran_flat(self->rng, site_left, site_right);
+                    position = msp_gsl_ran_flat(self->rng, site_left, site_right);
                     if (discrete_sites) {
                         position = floor(position);
                     }
                     search.position = position;
                     avl_node = avl_search(&self->sites, &search);
                 } while (avl_node != NULL && !discrete_sites);
 
-                time = gsl_ran_flat(self->rng, branch_start, branch_end);
+                time = msp_gsl_ran_flat(self->rng, branch_start, branch_end);
                 tsk_bug_assert(site_left <= position && position < site_right);
                 tsk_bug_assert(branch_start <= time && time < branch_end);
                 if (avl_node != NULL) {

lib/tests/test_core.c

@@ -130,6 +130,39 @@ test_tskit_version(void)
     CU_ASSERT_EQUAL(TSK_VERSION_PATCH, 15);
 }
 
+static void
+test_gsl_ran_flat_patch(void)
+{
+    double left, right, x;
+    int i;
+    gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
+
+    right = 85190021.000000;
+    left = nexttoward(right, 0.0);
+    gsl_rng_set(rng, 12345);
+    x = gsl_ran_flat(rng, left, right);
+    CU_ASSERT(!(right > x));
+    gsl_rng_set(rng, 12345);
+    x = msp_gsl_ran_flat(rng, left, right);
+    CU_ASSERT(right > x);
+    CU_ASSERT(left <= x);
+
+    left = 0.0;
+    right = 1.0;
+    for (i = 0; i < 1000; i++) {
+        x = msp_gsl_ran_flat(rng, left, right);
+        CU_ASSERT(left <= x);
+        CU_ASSERT(right > x);
+    }
+
+    left = 2.0;
+    right = 2.0;
+    x = msp_gsl_ran_flat(rng, left, right);
+    CU_ASSERT(left == x && right == x);
+
+    gsl_rng_free(rng);
+}
+
 int
 main(int argc, char **argv)
 {
@@ -138,6 +171,7 @@ main(int argc, char **argv)
         { "test_strerror_tskit", test_strerror_tskit },
         { "test_probability_list_select", test_probability_list_select },
         { "test_tskit_version", test_tskit_version },
+        { "test_gsl_ran_flat_patch", test_gsl_ran_flat_patch },
         CU_TEST_INFO_NULL,
     };
 

lib/util.c

@@ -25,6 +25,7 @@
 #include <errno.h>
 
 #include <gsl/gsl_math.h>
+#include <gsl/gsl_randist.h>
 
 #include <tskit/core.h>
 #include "util.h"
@@ -712,3 +713,28 @@ fast_search_free(fast_search_t *self)
 extern inline size_t sub_idx_1st_strict_upper_bound(
     const double *base, size_t start, size_t stop, double query);
 extern inline size_t fast_search_idx_strict_upper(fast_search_t *self, double query);
+
+/*   The gsl_ran_flat() function is supposed to output lo<=x<hi, but
+ *   sometimes (with probability <1e-9) we find x=hi.
+ *   This function includes a while loop to ensure x<hi.
+ *   https://github.com/tskit-dev/msprime/issues/1997
+ */
+double
+msp_gsl_ran_flat(gsl_rng *rng, double lo, double hi)
+{
+    double x;
+
+    /* The logic here requires that lo <= hi, but */
+    /* gsl_ran_flat produces output in this case. Keep */
+    /* this assert to make sure that msprime's usage */
+    /* has the required properties.  */
+    tsk_bug_assert(lo <= hi);
+    if (lo == hi) {
+        x = lo;
+    } else {
+        do {
+            x = gsl_ran_flat(rng, lo, hi);
+        } while (x >= hi);
+    }
+    return x;
+}

lib/util.h

@@ -24,6 +24,8 @@
 #include <assert.h>
 #include <stdio.h>
 
+#include <gsl/gsl_randist.h>
+
 #include "tskit.h"
 
 #ifdef __GNUC__
@@ -166,6 +168,8 @@ int fast_search_alloc(fast_search_t *self, const double *values, size_t n_values
 int fast_search_free(fast_search_t *self);
 inline size_t fast_search_idx_strict_upper(fast_search_t *self, double query);
 
+double msp_gsl_ran_flat(gsl_rng *rng, double lo, double hi);
+
 /***********************************
  * INLINE FUNCTION IMPLEMENTATIONS *
  ***********************************/