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

Issue:
mutation generation bug
I'm running into a bug using msprime.sim_mutations(). This pops up very infrequently so it’s not easy to reproduce, at least in the tree sequences I’m working with. But I was able to find a combo of seed (18112) and tree sequence that will generate the error. The tree sequence can be downloaded [here](http://sesame.uoregon.edu/~chriscs/ )(330Mb). 

MRE:
```
import tskit
import msprime
ts = tskit.load("/home/chriscs/theone.trees")
ts = msprime.sim_mutations(ts, rate=1e-8, random_seed=18112)

```

The error message is: 
“Bug detected in lib/mutgen.c at line 1137. Please report this issue on GitHub, ideally with a reproducible example. (https://github.com/tskit-dev/msprime/issues)
Aborted”

@andrewkern helped me add a printf near the bug assert
```
+++ b/lib/mutgen.c
@@ -1134,6 +1134,7 @@ mutgen_place_mutations(mutgen_t *self, bool discrete_sites)
                 } while (avl_node != NULL && !discrete_sites);

                 time = gsl_ran_flat(self->rng, branch_start, branch_end);
+                printf("DEBUG -- site_left: %lf position: %lf site_right: %lf\n", site_left, position, site_right);
                 tsk_bug_assert(site_left <= position && position < site_right);
                 tsk_bug_assert(branch_start <= time && time < branch_end);
                 if (avl_node != NULL) {
```

The following output shows the variable `position` from the gsl random number generator is hitting `site_right` 
DEBUG -- site_left: 66196102.000000 position: 66196121.000000 site_right: 66196127.000000
DEBUG -- site_left: 66195597.000000 position: 66195999.000000 site_right: 66196083.000000
DEBUG -- site_left: 28918378.000000 position: 28918563.000000 site_right: 28918696.000000
DEBUG -- site_left: 28918378.000000 position: 28918556.000000 site_right: 28918696.000000
DEBUG -- site_left: 28918378.000000 position: 28918623.000000 site_right: 28918696.000000
DEBUG -- site_left: 28918378.000000 position: 28918648.000000 site_right: 28918696.000000
DEBUG -- site_left: 72983639.000000 position: 72983787.000000 site_right: 72983828.000000
DEBUG -- site_left: 72983532.000000 position: 72983540.000000 site_right: 72983543.000000
DEBUG -- site_left: 72983639.000000 position: 72983811.000000 site_right: 72983828.000000
DEBUG -- site_left: 82717896.000000 position: 82717915.000000 site_right: 82717930.000000
DEBUG -- site_left: 82717896.000000 position: 82717914.000000 site_right: 82717930.000000
DEBUG -- site_left: 85189995.000000 **position: 85190021.000000 site_right: 85190021.000000**



Discussion / hints:
in debugging this with @chriscrsmith it looks like `gsl_ran_flat()` is _very rarely_ returning the top limit of its range, even though the documentation says this should not happen
Deleted my previous comment when I read your code in more detail. Something's fishy.
Instead of printing out the suspicious value, could you print out this instead:

(a < b) || (a > b)

where a and b are the position and the right boundary.


Well, it is a GSL bug.  This is a bit nuts--I can't find a bug report in this function since 2006.

 Here's the proof that it is in GSL.

Modify the msprime code to write the rng state to file.  When the error is detected, also write some the lo/hi values to file:

```c
    FILE * fp=NULL;
    FILE * fp2=NULL;
    void * rng_state;

    for (j = 0; j < edges.num_rows; j++) {
        left = edges.left[j];
        edge_right = edges.right[j];
        parent = edges.parent[j];
        child = edges.child[j];
        tsk_bug_assert(child >= 0 && child < (tsk_id_t) nodes.num_rows);
        branch_start = GSL_MAX(start_time, nodes.time[child]);
        branch_end = GSL_MIN(end_time, nodes.time[parent]);
        branch_length = branch_end - branch_start;

        map_index = rate_map_get_index(&self->rate_map, left);
        right = 0;
        while (right != edge_right) {
            right = GSL_MIN(edge_right, map_position[map_index + 1]);
            site_left = discrete_sites ? ceil(left) : left;
            site_right = discrete_sites ? ceil(right) : right;
            mu = branch_length * (site_right - site_left) * map_rate[map_index];
            branch_mutations = gsl_ran_poisson(self->rng, mu);
            for (k = 0; k < branch_mutations; k++) {
                /* Rejection sample positions until we get one we haven't seen before,
                 * unless we are doing discrete sites. Note that in principle this
                 * could lead to an infinite loop here, but in practise we'd need to
                 * use up all of the doubles before it could happen and so we'd
                 * certainly run out of memory first. */
                do {
                    fp=fopen("rng_state.bin","wb");
                    fp2=fopen("mutations.txt","w");
                    gsl_rng_fwrite(fp, self->rng);
                    position = gsl_ran_flat(self->rng, site_left, site_right);
                    if (!(position<site_right)) {
                        fprintf(stdout,"here\n");
                        fflush(stdout);
                        fwrite(&site_left,  sizeof(double),1, fp2);
                        fwrite(&site_right, sizeof(double),1, fp2);
                        fclose(fp);
                        fclose(fp2);
                        exit(0);
                    }
                    tsk_bug_assert(position < site_right);
                    if (discrete_sites) {
                        position = floor(position);
                    }
                    search.position = position;
                    avl_node = avl_search(&self->sites, &search);
                    fclose(fp);
                    fclose(fp2);
                } while (avl_node != NULL && !discrete_sites);
```

Then, a separate program to read the rng state, lo/hi, and make the next call:

```c
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>

int main(int argc, char **argv)
{
    // Do what msprime does.
    gsl_rng * r= gsl_rng_alloc(gsl_rng_default);;
    FILE * fp = fopen("rng_state.bin","rb");
    FILE * fp2 = fopen("mutations.txt","r");
    double a,b, pos;
    gsl_rng_fread(fp,r);
    fread(&a, sizeof(double), 1, fp2);
    fread(&b, sizeof(double), 1, fp2);

    fprintf(stdout, "%lf %lf\n", a, b);

    pos = gsl_ran_flat(r, a, b);
    fprintf(stdout, "%lf %lf %lf %d\n", a, b, pos, (b<pos)||(b>pos));
}
```

To compile the new code:

```sh
gcc -o test test.c -lgsl -lgslcblas
```

The output is:

```
./test
85189995.000000 85190021.000000
85189995.000000 85190021.000000 85190021.000000 0
```
Some system info:

```
msprime on  main [!?] via 🐍 v3.8.10 (fwdpy11_dev) 
❯ gsl-config --version
2.5
```
@andrewkern @chriscrsmith  -- I will be out of town for a few days.  Will figure out how to report a bug to GSL next week.  I have the files, etc., on my drive here.
@chriscrsmith, how often did you hit this? And, did it only happen with the one tree sequence, or with more than one?
@petrelharp this must be a very rare bug.  This GSL function has been called loads of times in msprime, slim, and fwdpy11.  Well not all of them may do runtime checks like msprime, if the bug were common they're probably a be a bunch of weird stuff happening.
while we sort this out, here is a patch on the rejection sampling that we could merge https://github.com/andrewkern/msprime/commit/683f598454651187dfe3315d2280c6bc05ecd813
> while we sort this out, here is a patch on the rejection sampling that we could merge https://github.com/andrewkern/msprime/commit/683f598454651187dfe3315d2280c6bc05ecd813

We should let @jeromekelleher weigh in tomorrow. I imagine we'll want audit the whole code base for calls to this function. I know I'm going to have to do the same on my end.
It might also be useful to check more recent versions of GSL. I'm using the one that comes on pop os 21.04, which is not the latest. However I did read the change logs and there's no indication of anything related to this being touched.  

The problem with waiting for a bug fix is that it takes ages in practice to trickle down to LTS Linux releases that are the basis of the docker images for doing things like the wheel building.
> It might also be useful to check more recent versions of GSL. I'm using the one that comes on pop os 21.04, which is not the latest. However I did read the change logs and there's no indication of anything related to this being touched.  
> 
> The problem with waiting for a bug fix is that it takes ages in practice to trickle down to LTS Linux releases that are the basis of the docker images for doing things like the wheel building.

Both msprime and fwdpy11 manually build GSL when building wheels because the version in the doctor image is too old.  But I know that I'm not using the very latest version in fwdpy11.  

It seems like we should be able to check the provenance information for an msprime tree sequence to see what version its wheels are using?
Well, I have to say I'm perversely delighted by this. A nasty bug that could have taken weeks to track down fully and completely described - thanks for making my life so easy here @chriscrsmith @andrewkern @molpopgen!

Given that it would probably take a decade or more for an updated version of GSL to make it onto all the platforms we're using, the answer is very straightforward I think:

- Create a function ``msp_gsl_rng_flag`` in util.c that looks like
```c
/* [Description of what we're doing, why we're doing it and a link to this issue] */
double 
msp_gsl_rng_flat(gsl_rng_t *rng, double lo, double hi)
{
    double x;
    do {
         x = gsl_rng_flat(rng, lo, hi);
    } while (x >= hi);
    return x
}
```
- Do a global s/gsl_rng_flat/msp_gsl_rng_flat/g in the lib directory

Done.

There's no point in waiting around for an updated GSL, I think, and it would create an annoying version dependency. RNG sampling is nowhere near the critical path in terms of performance, so this is an easy decision IMO.

Anyone want to pick this one up to claim the git-cred? (@chriscrsmith?)
I think there's a potential wrinkle in your prototype? I'm working from memory here, but if low equals high, doesn't the function return low?  That would give us an infinite loop.
I can work with @chriscrsmith on getting the patch you describe done @jeromekelleher 
@andrewkern -- be sure to add a test where lo == hi.  (See my comment above.)
Thanks for your help everyone. I'll work on this right away.

@petrelharp this was popping up every ~100,000 calls to the function for me, however only with specific sets of tree sequences—maybe particularly big tree sequences? I'm certain it's happened on more than one tree sequence, although I didn't save the other ones. I got this one by hitting my tree sequences with mutations for days on the cluster.
Great, thanks @chriscrsmith. A handful of tests for boundary cases would be good as well, as @molpopgen says. We can have a function ``test_gsl_flat_workaround`` or something in test_core.c for this.
The [source code](https://github.com/ampl/gsl/blob/dee1063dd57ec2841981b010eb732901b452c30a/randist/flat.c#L39) is very simple - it's just doing `a * (1 - u) + b * u` where `u` is a uniform between 0 and 1 drawn from `gsl_rng_uniform`. For *that* function, the documentation [explicitly states that 0 is included but 1 is excluded](https://www.gnu.org/software/gsl/doc/html/rng.html#c.gsl_rng_uniform), and indeed, if `u=0` then we get `a`, while if `u=1` then we'd get `b`, so `b` *should* be excluded.

 Here we've got `a = 85189995.000000` and `b = 85190021.000000`, and we're getting `b`; so, we must have `u = 1 - eps`. This seems like floating-point error... but, I don't actually see how floating point error could do this?
> The [source code](https://github.com/ampl/gsl/blob/dee1063dd57ec2841981b010eb732901b452c30a/randist/flat.c#L39) is very simple - it's just doing `a * (1 - u) + b * u` where `u` is a uniform between 0 and 1 drawn from `gsl_rng_uniform`. For _that_ function, the documentation [explicitly states that 0 is included but 1 is excluded](https://www.gnu.org/software/gsl/doc/html/rng.html#c.gsl_rng_uniform), and indeed, if `u=0` then we get `a`, while if `u=1` then we'd get `b`, so `b` _should_ be excluded.
> 
> Here we've got `a = 85189995.000000` and `b = 85190021.000000`, and we're getting `b`; so, we must have `u = 1 - eps`. This seems like floating-point error... but, I don't actually see how floating point error could do this?

This is the issue, and I'll try to look into it next week.  But, FP error is the likely culprit--the way the bits are being shuffled around is leading to a precision loss.  Roughly, it implies something like u*b > next_double_toward_minus_infinity(b), which returns b (b/c not double exist by definition other than b).
@petrelharp -- [these](https://en.cppreference.com/w/c/numeric/math/nextafter) functions help understand what I wrote above.  They are part of the C99 standard.
We've just written things on the board, and are suggesting that `a + u * (b - a)` might be better for this purpose; I've not got a good argument, though.
> We've just written things on the board, and are suggesting that `a + u * (b - a)` might be better for this purpose; I've not got a good argument, though.

It's totally possible that there are better ways to do this.  You can easily check what would happen in this case by:

1. Use C99 to get the largest double smaller than 1.0. Hint: it is not 1.0 - DBL_EPSILON, which is the delta to the first value larger than one.  You need the functions mentioned above.
2. Use that value and apply your function with a and b equal to what triggered the bug in this issue.

I don't have a computer with me for a few days, else I would have done it rather than writing this.

But even if it passes, that's not a sufficient proof that it's a generally better method.
> We've just written things on the board, and are suggesting that `a + u * (b - a)` might be better for this purpose; I've not got a good argument, though.

Oh, another issue is backwards compatibility given the same seed. That's kind of important. The rejection sampling approach would preserve backwards compatibility. This proposal changes the order of operations, and in general would change results for the same random number seed.
I'm with @molpopgen here - fixing this by coming up with a different method could easily end up with introducing other (more frequent) bugs. This is a one-in-a-gazillion thing, so let's just rejection sample it away.
A little note to think about:

There is a potential pathological case for the rejection sampling. The largest number of doubles in an interval of size 1 is between 0 and 1. So, if |b-a| is very small, and a,b >> 1, the rejection sampling risks taking a long time.

Thinking about this makes me wonder is such cases are an easy way to trigger the bug, returning b.

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9f1d2fa30..d0162583e 100644
--- a/CHANGELOG.md
+++ b/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**
diff --git a/lib/mutgen.c b/lib/mutgen.c
index 3bab35bbe..e33ebf869 100644
--- a/lib/mutgen.c
+++ b/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,7 +1125,7 @@ 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);
                     }
@@ -1133,7 +1133,7 @@ mutgen_place_mutations(mutgen_t *self, bool discrete_sites)
                     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) {
diff --git a/lib/util.c b/lib/util.c
index b4149ff61..c20f4bb78 100644
--- a/lib/util.c
+++ b/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;
+}
diff --git a/lib/util.h b/lib/util.h
index 0acbb90b6..8ae092f5a 100644
--- a/lib/util.h
+++ b/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 *
  ***********************************/


Test patch:
diff --git a/lib/tests/test_core.c b/lib/tests/test_core.c
index 3cc8ac461..ac0747e87 100644
--- a/lib/tests/test_core.c
+++ b/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,
     };