restore efficient genetic maps

cpp/fwdpy11_types/GeneticMapUnit.cc

@@ -1,33 +1,37 @@
+#include "fwdpy11/regions/RecombinationRegions.hpp"
 #include <pybind11/pybind11.h>
-#include <fwdpy11/regions/GeneticMapUnit.hpp>
+#include <core/genetic_maps/regions.hpp>
 
 namespace py = pybind11;
 
 void
 init_GeneticMapUnit(py::module& m)
 {
-    py::class_<fwdpy11::GeneticMapUnit>(m, "GeneticMapUnit");
+    py::class_<fwdpy11::PoissonCrossoverGenerator>(m, "PoissonCrossoverGenerator");
+    py::class_<fwdpy11::NonPoissonCrossoverGenerator>(m, "NonPoissonCrossoverGenerator");
 
-    py::class_<fwdpy11::PoissonInterval, fwdpy11::GeneticMapUnit>(m,
-                                                                  "_ll_PoissonInterval")
+    py::class_<fwdpy11_core::PoissonInterval, fwdpy11::PoissonCrossoverGenerator>(
+        m, "_ll_PoissonInterval")
         .def(py::init<double, double, double, bool>(), py::kw_only(), py::arg("beg"),
              py::arg("end"), py::arg("mean"), py::arg("discrete") = true);
 
-    py::class_<fwdpy11::PoissonPoint, fwdpy11::GeneticMapUnit>(m, "_ll_PoissonPoint")
+    py::class_<fwdpy11_core::PoissonPoint, fwdpy11::PoissonCrossoverGenerator>(
+        m, "_ll_PoissonPoint")
         .def(py::init<double, double, bool>(), py::kw_only(), py::arg("position"),
              py::arg("mean"), py::arg("discrete") = true);
 
-    py::class_<fwdpy11::BinomialInterval, fwdpy11::GeneticMapUnit>(
+    py::class_<fwdpy11_core::BinomialInterval, fwdpy11::NonPoissonCrossoverGenerator>(
         m, "_ll_BinomialInterval")
         .def(py::init<double, double, double, bool>(), py::kw_only(), py::arg("beg"),
              py::arg("end"), py::arg("probability"), py::arg("discrete") = true);
 
-    py::class_<fwdpy11::BinomialPoint, fwdpy11::GeneticMapUnit>(m, "_ll_BinomialPoint")
+    py::class_<fwdpy11_core::BinomialPoint, fwdpy11::NonPoissonCrossoverGenerator>(
+        m, "_ll_BinomialPoint")
         .def(py::init<double, double, bool>(), py::kw_only(), py::arg("position"),
              py::arg("probability"), py::arg("discrete") = true);
 
-    py::class_<fwdpy11::FixedCrossovers, fwdpy11::GeneticMapUnit>(m,
-                                                                  "_ll_FixedCrossovers")
+    py::class_<fwdpy11_core::FixedCrossovers, fwdpy11::NonPoissonCrossoverGenerator>(
+        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);
 }

cpp/regions/RecombinationRegions.cc

@@ -16,16 +16,7 @@ init_RecombinationRegions(py::module& m)
         .def_readonly("weights", &fwdpy11::RecombinationRegions::weights);
 
     py::class_<fwdpy11::GeneralizedGeneticMap, fwdpy11::GeneticMap>(
-        m, "GeneralizedGeneticMap")
-        .def(py::init([](py::list l) {
-            std::vector<std::unique_ptr<fwdpp::genetic_map_unit>> callbacks;
-            for (auto& i : l)
-                {
-                    auto& ref = i.cast<fwdpp::genetic_map_unit&>();
-                    callbacks.emplace_back(ref.clone());
-                }
-            return fwdpy11::GeneralizedGeneticMap(std::move(callbacks));
-        }));
+        m, "GeneralizedGeneticMap");
 
     m.def("dispatch_create_GeneticMap",
           [](py::object o, std::vector<fwdpy11::Region>& regions) {
@@ -36,15 +27,27 @@ init_RecombinationRegions(py::module& m)
               return fwdpy11::RecombinationRegions(o.cast<double>(), regions);
           });
 
-    m.def("dispatch_create_GeneticMap", [](py::object, py::list l) {
-        std::vector<std::unique_ptr<fwdpp::genetic_map_unit>> callbacks;
-        for (auto& i : l)
-            {
-                auto& ref = i.cast<fwdpy11::GeneticMapUnit&>();
-                callbacks.emplace_back(ref.ll_clone());
-            }
-        std::unique_ptr<fwdpy11::GeneralizedGeneticMap> rv(
-            new fwdpy11::GeneralizedGeneticMap(std::move(callbacks)));
-        return rv;
-    });
+    m.def("dispatch_create_GeneticMap_non_Region",
+          [](py::list poisson, py::list non_poisson) {
+              std::vector<std::unique_ptr<fwdpy11::PoissonCrossoverGenerator>>
+                  poisson_callbacks;
+              std::vector<std::unique_ptr<fwdpy11::NonPoissonCrossoverGenerator>>
+                  non_poisson_callbacks;
+              double x = 0.0;
+              for (auto& i : poisson)
+                  {
+                      auto& ref = i.cast<fwdpy11::PoissonCrossoverGenerator&>();
+                      x += ref.mean_number_xovers();
+                      poisson_callbacks.emplace_back(ref.ll_clone());
+                  }
+              for (auto& i : non_poisson)
+                  {
+                      auto& ref = i.cast<fwdpy11::NonPoissonCrossoverGenerator&>();
+                      non_poisson_callbacks.emplace_back(ref.ll_clone());
+                  }
+              std::unique_ptr<fwdpy11::GeneralizedGeneticMap> rv(
+                  new fwdpy11::GeneralizedGeneticMap(std::move(poisson_callbacks),
+                                                     std::move(non_poisson_callbacks)));
+              return rv;
+          });
 }

cpp_neutral_benchmark/cpp_neutral_benchmark.cc

@@ -13,6 +13,7 @@
 #include <fwdpy11/evolvets/recorders.hpp>
 #include <fwdpy11/rng.hpp>
 #include <core/evolve_discrete_demes/evolvets.hpp>
+#include <core/genetic_maps/regions.hpp>
 #include <boost/program_options.hpp>
 
 namespace po = boost::program_options;
@@ -158,11 +159,11 @@ simulate(const command_line_options& options)
 {
     fwdpy11::DiploidPopulation pop(options.N, 1e7);
     fwdpy11::GSLrng_t rng(options.seed);
-    fwdpy11::PoissonInterval recombination_region(0., 1e7, options.xovers, true);
+    fwdpy11_core::PoissonInterval recombination_region(0., 1e7, options.xovers, true);
 
-    std::vector<std::unique_ptr<fwdpp::genetic_map_unit>> callbacks;
+    std::vector<std::unique_ptr<fwdpy11::PoissonCrossoverGenerator>> callbacks;
     callbacks.emplace_back(recombination_region.ll_clone());
-    fwdpy11::GeneralizedGeneticMap genetic_map(std::move(callbacks));
+    fwdpy11::GeneralizedGeneticMap genetic_map(std::move(callbacks), {});
 
     fwdpy11::MutationRegions mmodel({}, {});
 

fwdpy11/_evolvets.py

@@ -184,8 +184,10 @@ def evolvets(
             )
 
     for r in params.recregions:
-        assert isinstance(r, fwdpy11.Region) or \
-            isinstance(r, fwdpy11._fwdpy11.GeneticMapUnit)
+        assert isinstance(r, fwdpy11._fwdpy11.PoissonCrossoverGenerator) or \
+            isinstance(
+                r, fwdpy11._fwdpy11.NonPoissonCrossoverGenerator) or \
+            isinstance(r, fwdpy11.Region), f"{type(r)}"
         if r.beg < 0:
             raise ValueError(f"{r} has begin value < 0.0")
         if r.end > pop.tables.genome_length:
@@ -214,6 +216,7 @@ def evolvets(
     from ._fwdpy11 import (
         MutationRegions,
         dispatch_create_GeneticMap,
+        dispatch_create_GeneticMap_non_Region,
         evolve_with_tree_sequences,
     )
 
@@ -230,8 +233,16 @@ def evolvets(
             params.rates.neutral_mutation_rate + params.rates.selected_mutation_rate
         )
     mm = MutationRegions.create(pneutral, params.nregions, params.sregions)
-    rm = dispatch_create_GeneticMap(
-        params.rates.recombination_rate, params.recregions)
+
+    if all([isinstance(i, fwdpy11.Region) for i in params.recregions]):
+        rm = dispatch_create_GeneticMap(
+            params.rates.recombination_rate, params.recregions)
+    else:
+        poisson = [i for i in params.recregions if isinstance(
+            i, fwdpy11._fwdpy11.PoissonCrossoverGenerator)]
+        non_poisson = [i for i in params.recregions if isinstance(
+            i, fwdpy11._fwdpy11.NonPoissonCrossoverGenerator)]
+        rm = dispatch_create_GeneticMap_non_Region(poisson, non_poisson)
 
     from ._fwdpy11 import SampleRecorder, _evolve_with_tree_sequences_options
 

fwdpy11/_types/model_params.py

@@ -57,9 +57,11 @@ class MutationAndRecombinationRates(object):
     @selected_mutation_rate.validator
     def validate_individual_rate(self, attribute, value):
         if not np.isfinite(value):
-            raise ValueError(f"{attribute} must be finite, but we got {value} instead")
+            raise ValueError(
+                f"{attribute} must be finite, but we got {value} instead")
         if value < 0.0:
-            raise ValueError(f"{attribute} must be >= 0.0, but we got {value} instead")
+            raise ValueError(
+                f"{attribute} must be >= 0.0, but we got {value} instead")
 
     @recombination_rate.validator
     def validate_recombination_rate(self, attribute, value):
@@ -208,7 +210,7 @@ def validate_recregions(self, attribute, value):
         except TypeError:
             try:
                 for i in value:
-                    attr.validators.instance_of(fwdpy11.GeneticMapUnit)(
+                    attr.validators.instance_of(fwdpy11.PoissonCrossoverGenerator)(
                         self, attribute, i
                     )
                 if not all([i.discrete for i in value]) and not all(
@@ -225,7 +227,7 @@ def validate_recregions(self, attribute, value):
     def rates_validator(self, attribute, value):
         if value.recombination_rate is None:
             for i in self.recregions:
-                if not isinstance(i, fwdpy11.GeneticMapUnit):
+                if not isinstance(i, fwdpy11.PoissonCrossoverGenerator) and not isinstance(i, fwdpy11.NonPoissonCrossoverGenerator):
                     raise ValueError(
                         f"recombination rate of {value.recombination_rate}"
                         " must be paired with"
@@ -267,4 +269,5 @@ def validate_demography(self, attribute, value):
     @simlen.validator
     def validate_simlen(self, attribute, value):
         if value <= 0:
-            raise ValueError(f"{attribute} must be >= 0, but we got {value} instead")
+            raise ValueError(
+                f"{attribute} must be >= 0, but we got {value} instead")

fwdpy11/headers/fwdpy11/regions/RecombinationRegions.hpp

@@ -1,6 +1,7 @@
 #ifndef FWDPY11_RECOMBINATIONREGIONS_HPP
 #define FWDPY11_RECOMBINATIONREGIONS_HPP
 
+#include <cstdlib>
 #include <limits>
 #include <vector>
 #include <algorithm>
@@ -10,6 +11,9 @@
 #include <gsl/gsl_randist.h>
 #include <fwdpy11/rng.hpp>
 #include "Region.hpp"
+#include "fwdpp/util/validators.hpp"
+#include "fwdpy11/gsl/gsl_error_handler_wrapper.hpp"
+#include "fwdpy11/regions/GeneticMapUnit.hpp"
 
 namespace fwdpy11
 {
@@ -71,21 +75,66 @@ namespace fwdpy11
         }
     };
 
+    struct PoissonCrossoverGenerator
+    {
+        // Must return exactly 1 value per call
+        virtual void breakpoint(const GSLrng_t& rng, std::vector<double>& breakpoints)
+            = 0;
+        virtual double left() = 0;
+        virtual double right() = 0;
+        virtual double mean_number_xovers() = 0;
+        virtual ~PoissonCrossoverGenerator() = default;
+        virtual std::unique_ptr<PoissonCrossoverGenerator> ll_clone() = 0;
+    };
+
+    struct NonPoissonCrossoverGenerator
+    {
+        virtual void breakpoint(const GSLrng_t& rng, std::vector<double>& breakpoints)
+            = 0;
+        virtual double left() = 0;
+        virtual double right() = 0;
+        virtual ~NonPoissonCrossoverGenerator() = default;
+        virtual std::unique_ptr<NonPoissonCrossoverGenerator> ll_clone() = 0;
+    };
+
     struct GeneralizedGeneticMap : public GeneticMap
     {
-        std::vector<std::unique_ptr<fwdpp::genetic_map_unit>> callbacks;
-        GeneralizedGeneticMap(std::vector<std::unique_ptr<fwdpp::genetic_map_unit>> c)
-            : callbacks(std::move(c))
+        std::vector<std::unique_ptr<PoissonCrossoverGenerator>> poisson_callbacks;
+        std::vector<std::unique_ptr<NonPoissonCrossoverGenerator>> non_poisson_callbacks;
+        fwdpp::gsl_ran_discrete_t_ptr poisson_lookup;
+        double sum_poisson_means;
+        GeneralizedGeneticMap(
+            std::vector<std::unique_ptr<PoissonCrossoverGenerator>> pc,
+            std::vector<std::unique_ptr<NonPoissonCrossoverGenerator>> nc)
+            : poisson_callbacks(std::move(pc)), non_poisson_callbacks(std::move(nc)),
+              poisson_lookup(nullptr), sum_poisson_means(0.0)
         {
+            std::vector<double> means;
+            for (auto& i : poisson_callbacks)
+                {
+                    sum_poisson_means += i->mean_number_xovers();
+                    means.push_back(i->mean_number_xovers());
+                }
+            if (!means.empty())
+                {
+                    poisson_lookup.reset(
+                        gsl_ran_discrete_preproc(means.size(), means.data()));
+                }
         }
 
         std::vector<double>
         operator()(const GSLrng_t& rng) const final
         {
             std::vector<double> rv;
-            for (auto&& c : callbacks)
+            auto nc = gsl_ran_poisson(rng.get(), sum_poisson_means);
+            for (unsigned i = 0; i < nc; ++i)
+                {
+                    auto region = gsl_ran_discrete(rng.get(), poisson_lookup.get());
+                    poisson_callbacks[region]->breakpoint(rng, rv);
+                }
+            for (auto& i : non_poisson_callbacks)
                 {
-                    c->operator()(rng.get(), rv);
+                    i->breakpoint(rng, rv);
                 }
             if (!rv.empty())
                 {

lib/CMakeLists.txt

@@ -18,9 +18,13 @@ set(EVOLVE_DISCRETE_DEMES_SOURCES
     evolve_discrete_demes/runtime_checks.cc
     evolve_discrete_demes/util.cc)
 
+set(GENETIC_MAP_SOURCES
+    genetic_maps/regions.cc)
+
 set(ALL_SOURCES ${DISCRETE_DEMOGRAPHY_SOURCES}
     ${MUTATION_DOMINANCE_SOURCES}
     ${DEMES_SOURCES}
+    ${GENETIC_MAP_SOURCES}
     ${EVOLVE_DISCRETE_DEMES_SOURCES})
 
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14")