Add BinomialIntervalMap (#1142)

Closes #1141

cpp/fwdpy11_types/GeneticMapUnit.cc

@@ -1,5 +1,6 @@
 #include "fwdpy11/regions/RecombinationRegions.hpp"
 #include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
 #include <core/genetic_maps/regions.hpp>
 
 namespace py = pybind11;
@@ -34,4 +35,14 @@ init_GeneticMapUnit(py::module& m)
         m, "_ll_FixedCrossovers")
         .def(py::init<double, double, int, bool>(), py::kw_only(), py::arg("beg"),
              py::arg("end"), py::arg("num_xovers"), py::arg("discrete") = true);
+
+    py::class_<fwdpy11_core::BinomialIntervalMap, fwdpy11::NonPoissonCrossoverGenerator>(
+        m, "_ll_BinomialIntervalMap")
+        .def(py::init([](double probability, const std::vector<fwdpy11::Region>& regions,
+                         bool discrete) {
+                 return fwdpy11_core::BinomialIntervalMap(probability, discrete,
+                                                          regions);
+             }),
+             py::kw_only(), py::arg("probability"), py::arg("regions"),
+             py::arg("discrete") = true);
 }

cpptests/test_core_genetic_map_regions.cc

@@ -1,11 +1,13 @@
 #include "fwdpy11/regions/RecombinationRegions.hpp"
+#include "fwdpy11/regions/Region.hpp"
 #include <boost/test/tools/old/interface.hpp>
 #include <boost/test/unit_test.hpp>
 
 #include <boost/test/unit_test_suite.hpp>
 #include <core/genetic_maps/regions.hpp>
 #include <fwdpy11/rng.hpp>
 #include <limits>
+#include <memory>
 
 BOOST_AUTO_TEST_SUITE(test_core_genetic_map_regions)
 
@@ -144,4 +146,23 @@ BOOST_AUTO_TEST_CASE(test_poisson_point)
         bp.size() - 1);
 }
 
+BOOST_AUTO_TEST_CASE(test_binomial_interval_map)
+{
+    std::vector<fwdpy11::Region> regions;
+    regions.push_back(fwdpy11::Region(0, 10, 1e-3, false, 0));
+    regions.push_back(fwdpy11::Region(10, 20, 5e-3, true, 0));
+    auto b = fwdpy11_core::BinomialIntervalMap(0.5, true, regions);
+    BOOST_CHECK_EQUAL(b.left(), 0.0);
+    BOOST_CHECK_EQUAL(b.right(), 20.0);
+    std::vector<std::unique_ptr<fwdpy11::NonPoissonCrossoverGenerator>> callbacks;
+    callbacks.push_back(
+        std::make_unique<fwdpy11_core::BinomialIntervalMap>(std::move(b)));
+    auto map = fwdpy11::GeneralizedGeneticMap({}, std::move(callbacks));
+    auto rng = fwdpy11::GSLrng_t(42);
+    for (int i = 0; i < 100; ++i)
+        {
+            auto bp = map(rng);
+        }
+}
+
 BOOST_AUTO_TEST_SUITE_END()

doc/pages/regiontypes.md

@@ -74,6 +74,11 @@
 
 ```
 
+```{eval-rst}
+.. autoclass:: fwdpy11.BinomialIntervalMap
+
+```
+
 ```{eval-rst}
 .. autoclass:: fwdpy11.MutationDominance
     

fwdpy11/__init__.py

@@ -35,6 +35,7 @@
     PoissonPoint,
     PoissonInterval,
     BinomialInterval,
+    BinomialIntervalMap,
     BinomialPoint,
     FixedCrossovers,
 )

fwdpy11/genetic_map_unit.py

@@ -313,3 +313,45 @@ def __attrs_post_init__(self):
             num_xovers=self.num_xovers,
             discrete=self.discrete,
         )
+
+
+@attr_add_asblack
+@attr_class_pickle_with_super
+@attr_class_to_from_dict
+@attr.s(auto_attribs=True, frozen=True, repr_ns="fwdpy11")
+class BinomialIntervalMap(fwdpy11._fwdpy11._ll_BinomialIntervalMap):
+    """
+    Generate exactly one crossover with a given probability
+    and a non-uniform map to generate positions.
+
+    This class has the following attributes, whose names
+    are also `kwargs` for intitialization.  The attribute names
+    also determine the order of positional arguments:
+
+    :param probability: The probability of a breakpoint
+    :type probability: float
+    :param regions: the map describing intervals and weights
+    :type regions: list[fwdpy11.Region]
+    :param discrete: If `False`, positions are continuous
+                     and uniform from `[beg, end)`.
+                     If `True`, positions take integer values
+                     uniformly from `[beg, end)`.
+    :type discrete: bool
+
+    Unlike :class:`fwdpy11.BinomialInterval`, this class allows
+    a non-uniform map to generate breakpoint positions.
+    The ``regions`` parameter is used to generate a lookup
+    table allowing efficient generation of positions.
+
+    .. versionadded:: 0.20.0
+    """
+    probability: float
+    regions: typing.List[fwdpy11.Region]
+    discrete: bool = attr.ib(kw_only=True, default=False)
+
+    def __attrs_post_init__(self):
+        super(BinomialIntervalMap, self).__init__(
+            probability=self.probability,
+            regions=self.regions,
+            discrete=self.discrete
+        )

lib/core/genetic_maps/regions.hpp

@@ -1,5 +1,7 @@
 #pragma once
 
+#include "fwdpp/gsl_discrete.hpp"
+#include "fwdpy11/regions/Region.hpp"
 #include "fwdpy11/rng.hpp"
 #include <algorithm>
 #include <gsl/gsl_randist.h>
@@ -105,6 +107,28 @@ namespace fwdpy11_core
         virtual std::unique_ptr<fwdpy11::NonPoissonCrossoverGenerator> ll_clone();
     };
 
+    struct BinomialIntervalMap : public fwdpy11::NonPoissonCrossoverGenerator
+    {
+        struct Segment
+        {
+            double left;
+            double right;
+        };
+        double probability;
+        bool discrete;
+        fwdpp::gsl_ran_discrete_t_ptr lookup;
+        std::vector<fwdpy11::Region> regions;
+        std::vector<Segment> segments;
+        BinomialIntervalMap(double probability, bool discrete,
+                            const std::vector<fwdpy11::Region>&);
+
+        virtual double left();
+        virtual double right();
+        virtual void breakpoint(const fwdpy11::GSLrng_t& rng,
+                                std::vector<double>& breakpoints);
+        virtual std::unique_ptr<fwdpy11::NonPoissonCrossoverGenerator> ll_clone();
+    };
+
     struct BinomialPoint : public fwdpy11::NonPoissonCrossoverGenerator
     {
         double position;

lib/genetic_maps/regions.cc

@@ -1,5 +1,6 @@
 #include "fwdpy11/regions/RecombinationRegions.hpp"
 #include <cmath>
+#include <gsl/gsl_randist.h>
 #include <limits>
 #include <memory>
 #include <stdexcept>
@@ -130,4 +131,73 @@ namespace fwdpy11_core
         return std::unique_ptr<fwdpy11::NonPoissonCrossoverGenerator>(
             new FixedCrossovers(*this));
     }
+
+    BinomialIntervalMap::BinomialIntervalMap(double probability, const bool discrete,
+                                             const std::vector<fwdpy11::Region>& regions)
+        : probability(probability), discrete{discrete}, lookup{nullptr},
+          regions{regions}, segments{}
+    {
+        validate_parameter(probability, 0.0, 1.0);
+        std::vector<double> weights;
+        weights.reserve(regions.size());
+        segments.reserve(regions.size());
+        for (const auto& r : this->regions)
+            {
+                validate_interval(r.beg, r.end, discrete);
+                auto w = r.weight;
+                if (r.coupled)
+                    {
+                        w += (r.end - r.beg);
+                    }
+                weights.push_back(w);
+                segments.push_back(Segment{r.beg, r.end});
+            }
+        lookup.reset(gsl_ran_discrete_preproc(segments.size(), weights.data()));
+    }
+
+    double
+    BinomialIntervalMap::left()
+    {
+        auto rv = std::numeric_limits<double>::max();
+        for (const auto& s : this->segments)
+            {
+                rv = std::min(s.left, rv);
+            }
+        return rv;
+    }
+
+    double
+    BinomialIntervalMap::right()
+    {
+        auto rv = std::numeric_limits<double>::min();
+        for (const auto& s : this->segments)
+            {
+                rv = std::max(s.right, rv);
+            }
+        return rv;
+    }
+
+    void
+    BinomialIntervalMap::breakpoint(const fwdpy11::GSLrng_t& rng,
+                                    std::vector<double>& breakpoints)
+    {
+        if (gsl_rng_uniform(rng.get()) <= probability)
+            {
+                auto index = gsl_ran_discrete(rng.get(), this->lookup.get());
+                auto bp = gsl_ran_flat(rng.get(), this->segments[index].left,
+                                       this->segments[index].right);
+                if (discrete)
+                    {
+                        bp = std::floor(bp);
+                    }
+                breakpoints.push_back(bp);
+            }
+    }
+
+    std::unique_ptr<fwdpy11::NonPoissonCrossoverGenerator>
+    BinomialIntervalMap::ll_clone()
+    {
+        return std::unique_ptr<fwdpy11::NonPoissonCrossoverGenerator>(
+            new BinomialIntervalMap(this->probability, this->discrete, this->regions));
+    }
 }

tests/test_ModelParams_with_recregions.py

@@ -45,7 +45,9 @@ def _validate_expected(regions, numbers):
     ([fwdpy11.PoissonInterval(0, 0.1, 1e-3),
      fwdpy11.PoissonPoint(0.1, 1e-3),
      fwdpy11.BinomialInterval(0.1, 0.2, 1e-3),
-     fwdpy11.BinomialPoint(0.2, 1e-3)], (2, 2))
+     fwdpy11.BinomialPoint(0.2, 1e-3)], (2, 2)),
+    ([fwdpy11.BinomialIntervalMap(0.5, [fwdpy11.Region(0, 1, 1e-5),
+                                        fwdpy11.Region(1, 2, 2)])], (0, 1))
 ])
 def test_model_params_with_recregions(data):
     regions, expected = data