Use scope guard to promote GSL errors to exceptions

fwdpy11/headers/fwdpy11/gsl/gsl_error_handler_wrapper.hpp

@@ -1,11 +1,31 @@
 #ifndef FWDPY11_GSL_GSL_ERROR_HANDLER_WRAPPER_HPP
 #define FWDPY11_GSL_GSL_ERROR_HANDLER_WRAPPER_HPP
 
+#include <exception>
+#include <sstream>
+#include <string>
 #include <gsl/gsl_errno.h>
 
 namespace fwdpy11
 {
-    struct gsl_error_handler_wrapper
+
+    class __attribute__((visibility("default"))) GSLError : public std::exception
+    {
+      private:
+        std::string message_;
+
+      public:
+        explicit GSLError(std::string message) : message_(std::move(message))
+        {
+        }
+        virtual const char*
+        what() const noexcept
+        {
+            return message_.c_str();
+        }
+    };
+
+    class gsl_scoped_convert_error_to_exception
     /*!
      * Manages turning off and resetting the GSL error
      * handler.
@@ -14,10 +34,60 @@ namespace fwdpy11
      * so this class is like a "smart pointer" for
      * turning off the handler.
      */
+    {
+      private:
+        gsl_error_handler_t *e, *custom_handler;
+
+        static void
+        gsl_error_to_exception(const char* reason, const char* file, int line,
+                               int gsl_errno)
+        {
+            std::ostringstream o;
+            o << "GSL error raised: " << reason << ", " << file << ", " << line << ", "
+              << gsl_errno;
+            throw fwdpy11::GSLError(o.str());
+        }
+
+      public:
+        gsl_scoped_convert_error_to_exception()
+            : e{gsl_set_error_handler_off()},
+              custom_handler{gsl_set_error_handler(
+                  &gsl_scoped_convert_error_to_exception::gsl_error_to_exception)}
+        {
+        }
+
+        ~gsl_scoped_convert_error_to_exception()
+        {
+            gsl_set_error_handler(e);
+        }
+
+        gsl_scoped_convert_error_to_exception(
+            const gsl_scoped_convert_error_to_exception&)
+            = delete;
+        gsl_scoped_convert_error_to_exception&
+        operator=(const gsl_scoped_convert_error_to_exception&)
+            = delete;
+    };
+
+    struct gsl_scoped_disable_error_handler_wrapper
     {
         gsl_error_handler_t* e;
-        gsl_error_handler_wrapper() : e{ gsl_set_error_handler_off() } {}
-        ~gsl_error_handler_wrapper() { gsl_set_error_handler(e); }
+
+        gsl_scoped_disable_error_handler_wrapper() : e{gsl_set_error_handler_off()}
+        {
+        }
+
+        ~gsl_scoped_disable_error_handler_wrapper()
+        {
+            gsl_set_error_handler(e);
+        }
+
+        gsl_scoped_disable_error_handler_wrapper(
+            const gsl_scoped_disable_error_handler_wrapper&)
+            = delete;
+        gsl_scoped_disable_error_handler_wrapper&
+        operator=(const gsl_scoped_disable_error_handler_wrapper&)
+            = delete;
     };
 } // namespace fwdpy11
 

fwdpy11/headers/fwdpy11/regions/MultivariateGaussianEffects.hpp

@@ -12,6 +12,7 @@
 #include <functional>
 #include <memory>
 #include <fwdpy11/policies/mutation.hpp>
+#include <fwdpy11/gsl/gsl_error_handler_wrapper.hpp>
 #include "Sregion.hpp"
 
 namespace fwdpy11
@@ -70,32 +71,24 @@ namespace fwdpy11
                         "input matrix contains non-finite values");
                 }
 
+            gsl_scoped_disable_error_handler_wrapper gsl_error_scope_guard;
             // Assign the matrix and do the Cholesky decomposition
-            auto error_handler = gsl_set_error_handler_off();
 
             int rv = gsl_matrix_memcpy(matrix.get(), &input_matrix);
             if (rv != GSL_SUCCESS)
                 {
-                    // Reset error handler on the way out
-                    gsl_set_error_handler(error_handler);
                     throw std::runtime_error("failure copying input matrix");
                 }
             rv = gsl_matrix_memcpy(input_matrix_copy.get(), &input_matrix);
             if (rv != GSL_SUCCESS)
                 {
-                    // Reset error handler on the way out
-                    gsl_set_error_handler(error_handler);
                     throw std::runtime_error("failure copying input matrix");
                 }
             rv = gsl_linalg_cholesky_decomp1(matrix.get());
             if (rv == GSL_EDOM)
                 {
-                    // Reset error handler on the way out
-                    gsl_set_error_handler(error_handler);
                     throw std::invalid_argument("Cholesky decomposition failed");
                 }
-            // Reset error handler on the way out
-            gsl_set_error_handler(error_handler);
         }
 
         virtual std::unique_ptr<Sregion>

fwdpy11/src/_fwdpy11.cc

@@ -1,22 +1,13 @@
 #include <stdexcept>
-#include <sstream>
 #include <type_traits>
 #include <gsl/gsl_version.h>
 #include <gsl/gsl_errno.h>
+#include <fwdpy11/gsl/gsl_error_handler_wrapper.hpp>
 #include <pybind11/pybind11.h>
 
 static_assert(GSL_MAJOR_VERSION >= 2, "GSL major version >= 2 required");
 static_assert(GSL_MINOR_VERSION >= 3, "GSL minor version >= 3 required");
 
-void
-gsl_error_to_exception(const char *reason, const char *file, int line, int gsl_errno)
-{
-    std::ostringstream o;
-    o << "GSL error raised: " << reason << ", " << file << ", " << line << ", "
-      << gsl_errno;
-    throw std::runtime_error(o.str());
-}
-
 namespace py = pybind11;
 
 void initialize_fwdpp_types(py::module &);
@@ -35,17 +26,6 @@ void init_array_proxies(py::module &m);
 
 PYBIND11_MODULE(_fwdpy11, m)
 {
-    auto handler = gsl_set_error_handler_off();
-
-    auto custom_handler = gsl_set_error_handler(&gsl_error_to_exception);
-
-    // When the module exits, restore the old handler.
-    // We re-use variables to suppress unused variable warnings.
-    py::module::import("atexit").attr("register")(
-        py::cpp_function{[&handler, &custom_handler]() {
-            custom_handler = gsl_set_error_handler(handler);
-        }});
-
     initialize_fwdpp_types(m);
     initialize_fwdpp_functions(m);
     initialize_fwdpy11_types(m);
@@ -60,6 +40,8 @@ PYBIND11_MODULE(_fwdpy11, m)
     init_discrete_demography(m);
     init_array_proxies(m);
 
+    py::register_exception<fwdpy11::GSLError>(m, "GSLError");
+
     m.def(
         "pybind11_version",
         []() {

fwdpy11/src/evolve_population/with_tree_sequences.cc

@@ -39,6 +39,7 @@
 #include <fwdpy11/regions/MutationRegions.hpp>
 #include <fwdpy11/regions/RecombinationRegions.hpp>
 #include <fwdpy11/discrete_demography/DiscreteDemography.hpp>
+#include <fwdpy11/gsl/gsl_error_handler_wrapper.hpp>
 #include <fwdpy11/samplers.hpp>
 #include <fwdpp/ts/recycling.hpp>
 #include "util.hpp"
@@ -177,6 +178,8 @@ evolve_with_tree_sequences(
     const bool preserve_first_generation,
     const fwdpy11::DiploidPopulation_temporal_sampler &post_simplification_recorder)
 {
+    fwdpy11::gsl_scoped_convert_error_to_exception gsl_error_scope_guard;
+
     if (gvalue_pointers.genetic_values.empty())
         {
             throw std::invalid_argument("empty list of genetic values");

tests/gsl_error.cc

@@ -4,8 +4,8 @@
 
 PYBIND11_MODULE(gsl_error, m)
 {
-    m.def("trigger_error", []() {
-        fwdpy11::gsl_error_handler_wrapper h;
+    m.def("trigger_error_handle_locally", []() {
+        fwdpy11::gsl_scoped_disable_error_handler_wrapper h;
         gsl_matrix* m = gsl_matrix_alloc(2, 3); //non-square
         int rv = gsl_matrix_transpose(m);
         if (rv != GSL_SUCCESS)
@@ -16,10 +16,21 @@ PYBIND11_MODULE(gsl_error, m)
         gsl_matrix_free(m);
     });
 
-    // Causes a memory leak!
-    m.def("trigger_error_not_handled", []() {
+    m.def("trigger_error", []() {
+        fwdpy11::gsl_scoped_convert_error_to_exception gsl_error_scope;
         gsl_matrix* m = gsl_matrix_alloc(2, 3); //non-square
-        gsl_matrix_transpose(m);
-        gsl_matrix_free(m);
+        try
+            {
+                gsl_matrix_transpose(m);
+                gsl_matrix_free(m);
+            }
+        catch (fwdpy11::GSLError& e)
+            {
+                throw e;
+            }
+        catch (...)
+            {
+                throw std::runtime_error("unexcpected exception");
+            }
     });
 }

tests/test_GSLerror.py

@@ -7,15 +7,15 @@
 class testGSLerror(unittest.TestCase):
     def testErrorHandler(self):
         with self.assertRaises(RuntimeError):
-            gsl_error.trigger_error()
+            gsl_error.trigger_error_handle_locally()
         try:
-            gsl_error.trigger_error()
+            gsl_error.trigger_error_handle_locally()
         except RuntimeError as e:
             self.assertEqual(str(e), "matrix not square")
 
     def testUnhandledGSLError(self):
-        with self.assertRaises(RuntimeError):
-            gsl_error.trigger_error_not_handled()
+        with self.assertRaises(fwdpy11.GSLError):
+            gsl_error.trigger_error()
 
 
 if __name__ == "__main__":

tests/test_TableCollection_fs.py

@@ -19,9 +19,9 @@
 
 import unittest
 
-import numpy as np
-
 import fwdpy11
+import msprime
+import numpy as np
 
 # NOTE: these tests can all be improved in the future
 # by having a conversion from msprime/tskit that lifts
@@ -42,8 +42,6 @@ def fs_from_ndarray(gm):
 class TestSingleDemeCase(unittest.TestCase):
     @classmethod
     def setUpClass(self):
-        import msprime
-
         Ne = 1000
         Nr = 100.0
         ts = msprime.simulate(
@@ -128,8 +126,6 @@ def test_separated_windows(self):
 class TestTwoDemeCase(unittest.TestCase):
     @classmethod
     def setUpClass(self):
-        import msprime
-
         Ne = 1000
         Nr = 100.0
         nodes_per_deme = 1000

tests/test_tree_sequences.py

@@ -23,11 +23,11 @@
 import unittest
 from collections import namedtuple
 
-import numpy as np
-import tskit
-
 import fwdpy11
 import fwdpy11.tskit_tools
+import msprime
+import numpy as np
+import tskit
 
 
 class Recorder(object):
@@ -1579,8 +1579,6 @@ def test_recapitation(self):
         Recapitation w/o recombination means we can predict
         exactly how many new nodes there will be.
         """
-        import msprime
-
         coalesced_ts = msprime.simulate(
             Ne=self.pop.N,
             from_ts=self.tskit_ts,
@@ -1602,8 +1600,6 @@ def test_invalid_preservation(self):
 
 class TestInvalidAttemptAtRecapitation(unittest.TestCase):
     def test_after_starting_from_msprime(self):
-        import msprime
-
         ts = msprime.simulate(100, Ne=100)
         pop = fwdpy11.DiploidPopulation.create_from_tskit(ts)
 

tests/test_ts_from_msprime.py

@@ -19,19 +19,14 @@
 
 import unittest
 
-import numpy as np
-
 import fwdpy11
-
-# NOTE: msprime is imported on a per-test basis so as not
-# to confict w/other tests of the GSL error handling machinery
+import msprime
+import numpy as np
 
 
 class TestConversion(unittest.TestCase):
     @classmethod
     def setUp(self):
-        import msprime
-
         self.ts = msprime.simulate(10, recombination_rate=0.025, Ne=1000)
 
     # def testGetTablesDiscretizeTime(self):
@@ -53,8 +48,6 @@ def testCreateDiploidPopulation(self):
 
 class TestConversionFromMultipleDemes(unittest.TestCase):
     def test_deme_field_of_metadata(self):
-        import msprime
-
         nodes_per_deme = 500
         Ne = 3 * nodes_per_deme // 2
         Nr = 50.0