add safety check to gsl_ran_flat.

cpp/functions/add_mutation.cc

@@ -26,6 +26,7 @@
 #include <fwdpy11/rng.hpp>
 #include <fwdpy11/types/DiploidPopulation.hpp>
 #include <fwdpy11/policies/mutation.hpp>
+#include <core/internal/gsl_ran_flat.hpp>
 #include <fwdpp/ts/marginal_tree_functions/nodes.hpp>
 #include <fwdpp/ts/table_collection_functions.hpp>
 #include <fwdpp/ts/marginal_tree_functions/samples.hpp>
@@ -243,8 +244,8 @@ namespace
                 return static_cast<std::int64_t>(std::ceil(mutation_node_time));
             }
         // choose a random time
-        auto candidate_time = gsl_ran_flat(rng.get(), mutation_node_time,
-                                           node_table[mutation_node_parent].time);
+        auto candidate_time = fwdpy11_core::internal::gsl_ran_flat(
+            rng, mutation_node_time, node_table[mutation_node_parent].time);
         if (candidate_time < 0.0)
             {
                 return static_cast<std::int64_t>(std::floor(mutation_node_time));
@@ -386,7 +387,8 @@ add_mutation(const fwdpy11::GSLrng_t& rng, const double left, const double right
         empty, pop.mutations, pop.mut_lookup, false, mutation_node_time,
         // The lambdas will simply transfer input parameters to the new object.
         [&rng, &candidate_data]() {
-            return gsl_ran_flat(rng.get(), candidate_data.left, candidate_data.right);
+            return fwdpy11_core::internal::gsl_ran_flat(rng, candidate_data.left,
+                                                        candidate_data.right);
         },
         [&data]() { return data.effect_size; },
         [&data](double) { return data.dominance; }, [&data]() { return data.esizes; },

cpp/gsl/gsl_random.cc

@@ -20,39 +20,43 @@
 #include <pybind11/pybind11.h>
 #include <fwdpy11/rng.hpp>
 #include <gsl/gsl_randist.h>
+#include <core/internal/gsl_ran_flat.hpp>
 
 namespace py = pybind11;
 
 void
 init_gsl_random(py::module& m)
 {
-    m.def("gsl_ran_gaussian_ziggurat",
-          [](const fwdpy11::GSLrng_t& rng, const double sd) {
-              return gsl_ran_gaussian_ziggurat(rng.get(), sd);
-          },
-          "Gaussian deviate with mean zero");
+    m.def(
+        "gsl_ran_gaussian_ziggurat",
+        [](const fwdpy11::GSLrng_t& rng, const double sd) {
+            return gsl_ran_gaussian_ziggurat(rng.get(), sd);
+        },
+        "Gaussian deviate with mean zero");
 
-    m.def("gsl_rng_uniform",
-          [](const fwdpy11::GSLrng_t& rng) {
-              return gsl_rng_uniform(rng.get());
-          },
-          "Uniform deviate on interval [0,1)");
+    m.def(
+        "gsl_rng_uniform",
+        [](const fwdpy11::GSLrng_t& rng) { return gsl_rng_uniform(rng.get()); },
+        "Uniform deviate on interval [0,1)");
 
-    m.def("gsl_ran_flat",
-          [](const fwdpy11::GSLrng_t& rng, double a, double b) {
-              return gsl_ran_flat(rng.get(), a, b);
-          },
-          "Unform deviate on interval [a,b)");
+    m.def(
+        "gsl_ran_flat",
+        [](const fwdpy11::GSLrng_t& rng, double a, double b) {
+            return fwdpy11_core::internal::gsl_ran_flat(rng, a, b);
+        },
+        "Unform deviate on interval [a,b)");
 
-    m.def("gsl_ran_poisson",
-          [](const fwdpy11::GSLrng_t& rng, double mean) {
-              return gsl_ran_poisson(rng.get(), mean);
-          },
-          "Poisson deviate parameterized by mean.");
+    m.def(
+        "gsl_ran_poisson",
+        [](const fwdpy11::GSLrng_t& rng, double mean) {
+            return gsl_ran_poisson(rng.get(), mean);
+        },
+        "Poisson deviate parameterized by mean.");
 
-    m.def("gsl_ran_geometric",
-          [](const fwdpy11::GSLrng_t& rng, const double p) {
-              return gsl_ran_geometric(rng.get(), p);
-          },
-          "Geometric distribution parameterized by success probability.");
+    m.def(
+        "gsl_ran_geometric",
+        [](const fwdpy11::GSLrng_t& rng, const double p) {
+            return gsl_ran_geometric(rng.get(), p);
+        },
+        "Geometric distribution parameterized by success probability.");
 }

cpp/ts/infinite_sites.cc

@@ -1,5 +1,6 @@
 #include <cmath>
 #include <algorithm>
+#include <core/internal/gsl_ran_flat.hpp>
 #include <fwdpp/ts/recycling.hpp>
 #include <fwdpp/ts/mutate_tables.hpp>
 #include <fwdpp/ts/count_mutations.hpp>
@@ -20,7 +21,7 @@ namespace
                               const fwdpy11::mutation_origin_time generation)
     {
         const auto uniform
-            = [left, right, &rng]() { return gsl_ran_flat(rng.get(), left, right); };
+            = [left, right, &rng]() { return fwdpy11_core::internal::gsl_ran_flat(rng, left, right); };
         const auto return_zero = []() { return 0.0; };
         const auto return_zero_h = [](const double) { return 0.0; };
         auto key = fwdpy11::infsites_Mutation(recycling_bin, pop.mutations,

fwdpy11/headers/fwdpy11/regions/Region.hpp

@@ -44,7 +44,12 @@ namespace fwdpy11
         inline double
         operator()(const GSLrng_t& rng) const
         {
-            return gsl_ran_flat(rng.get(), beg, end);
+            auto rv = gsl_ran_flat(rng.get(), beg, end);
+            while (rv == end)
+                {
+                    rv = gsl_ran_flat(rng.get(), beg, end);
+                }
+            return rv;
         }
 
         inline bool

fwdpy11/headers/fwdpy11/regions/UniformS.hpp

@@ -47,10 +47,15 @@ namespace fwdpy11
             return infsites_Mutation(
                 recycling_bin, mutations, lookup_table, false, generation,
                 [this, &rng]() { return region(rng); },
-                [this, &rng]() { return gsl_ran_flat(rng.get(), lo, hi) / scaling; },
-                [this, &rng](const double esize) {
-                    return dominance(rng, esize);
+                [this, &rng]() {
+                    auto rv = gsl_ran_flat(rng.get(), lo, hi);
+                    while (rv == hi)
+                        {
+                            rv = gsl_ran_flat(rng.get(), lo, hi);
+                        }
+                    return rv / scaling;
                 },
+                [this, &rng](const double esize) { return dominance(rng, esize); },
                 this->label());
         }
 

lib/core/internal/gsl_ran_flat.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <fwdpy11/rng.hpp>
+#include <gsl/gsl_randist.h>
+
+namespace fwdpy11_core
+{
+    namespace internal
+    {
+        inline auto
+        gsl_ran_flat(const fwdpy11::GSLrng_t& rng, double lo, double hi)
+        {
+            auto rv = gsl_ran_flat(rng.get(), lo, hi);
+            while (rv == hi)
+                {
+                    rv = gsl_ran_flat(rng.get(), lo, hi);
+                }
+            return rv;
+        }
+    }
+}

lib/mutation_dominance/dependency_injection.cc

@@ -1,6 +1,7 @@
 #include <cmath>
 #include <stdexcept>
 #include <functional>
+#include <core/internal/gsl_ran_flat.hpp>
 #include <gsl/gsl_randist.h> // FIXME: hide in static lib
 #include <fwdpp/util/validators.hpp>
 #include <fwdpy11/mutation_dominance/MutationDominance.hpp>
@@ -46,7 +47,7 @@ namespace fwdpy11
                 throw std::invalid_argument("hi must be > lo");
             }
         return MutationDominance([lo, hi](const GSLrng_t& rng, const double) {
-            return gsl_ran_flat(rng.get(), lo, hi);
+            return fwdpy11_core::internal::gsl_ran_flat(rng, lo, hi);
         });
     }