Add discretized joint distribution of effect size and dominance.

doc/pages/regiontypes.rst

@@ -23,6 +23,8 @@ Setting up genomic intervals
 
     .. automethod:: mv
 
+.. autoclass:: fwdpy11.DiscreteDESD
+
 .. autoclass:: fwdpy11.mvDES
 
 .. autoclass:: fwdpy11.GeneticMapUnit

doc/pages/tutorial.rst

@@ -222,6 +222,7 @@ Built-in distributions of effect sizes
 * :class:`fwdpy11.GaussianS`
 * :class:`fwdpy11.MultivariateGaussianEffects`
 * :class:`fwdpy11.LogNormalS`
+* :class:`fwdpy11.DiscreteDESD`
 
 Labelling mutations from different regions
 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

fwdpy11/CMakeLists.txt

@@ -1,4 +1,4 @@
-set(FWDPP_TYPES_SOURCES src/fwdpp_types/init.cc 
+set(FWDPP_TYPES_SOURCES src/fwdpp_types/init.cc ;
     src/fwdpp_types/MutationBase.cc
     src/fwdpp_types/HaploidGenome.cc
     src/fwdpp_types/DataMatrix.cc
@@ -45,7 +45,8 @@ set(REGION_SOURCES src/regions/init.cc src/regions/Region.cc src/regions/Sregion
     src/regions/ExpS.cc src/regions/UniformS.cc src/regions/GaussianS.cc src/regions/MutationRegions.cc
     src/regions/RecombinationRegions.cc src/regions/MultivariateGaussianEffects.cc
     src/regions/mvDES.cc
-    src/regions/LogNormalS.cc)
+    src/regions/LogNormalS.cc
+    src/regions/DiscreteDESD.cc)
 
 set(GSL_SOURCES src/gsl/init.cc
     src/gsl/gsl_random.cc)

fwdpy11/regions.py

@@ -21,9 +21,10 @@
 import warnings
 
 import attr
-import fwdpy11._fwdpy11
 import numpy as np
 
+import fwdpy11._fwdpy11
+
 from .class_decorators import (attr_add_asblack, attr_class_pickle_with_super,
                                attr_class_to_from_dict,
                                attr_class_to_from_dict_no_recurse)
@@ -684,3 +685,60 @@ class mvDES(fwdpy11._fwdpy11._ll_mvDES):
 
     def __attrs_post_init__(self):
         super(mvDES, self).__init__(self.des, self.means, self.matrix)
+
+
+@attr_add_asblack
+@attr_class_pickle_with_super
+@attr_class_to_from_dict_no_recurse
+@attr.s(eq=False, **_common_attr_attribs)
+class DiscreteDESD(fwdpy11._fwdpy11._ll_DiscreteDESD):
+    """
+    Discretized distribution of effect sizes and dominance.
+
+    This class allows you to specify a discrete joint
+    distrubtion of effect size and dominance.
+
+    The distribution is specified by a list of tuples.
+    Each tuple contains (effect size, dominance, weight).
+    The weights must all be >= 0 and all values must fe finite.
+
+    This class has the following attributes, whose names
+    are also ``kwargs`` for intitialization.  The attribute names
+    also determine the order of positional arguments:
+
+    :param beg: Beginning of the region
+    :type beg: float
+    :param end: End of the region
+    :type end: float
+    :param weight: Weight on the region
+    :type weight: float
+    :param joint_dist: The joint distribution + weights
+    :type joint_dist: list
+    :param coupled: Specify if weight is function of end-beg or not. Defaults to True
+    :type coupled: bool
+    :param label: Fill :attr:`fwdpy11.Mutation.label` with this value.
+    :type label: numpy.uint16
+    :param scaling: The scaling of the DFE
+    :type scaling: float
+
+    .. versionadded:: 0.10.0
+    """
+
+    beg: float
+    end: float
+    weight: float
+    joint_dist: typing.List[typing.Tuple[float, float, float]]
+    coupled: bool = True
+    label: int = 0
+    scaling: float = 1.0
+
+    def __attrs_post_init__(self):
+        super(DiscreteDESD, self).__init__(
+            self.beg,
+            self.end,
+            self.weight,
+            self.joint_dist,
+            self.coupled,
+            self.label,
+            self.scaling,
+        )

fwdpy11/src/regions/DiscreteDESD.cc

@@ -0,0 +1,129 @@
+#include <vector>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <fwdpp/gsl_discrete.hpp>
+#include <fwdpy11/regions/Sregion.hpp>
+#include <fwdpy11/policies/mutation.hpp>
+
+namespace py = pybind11;
+
+class DiscreteDESD : public fwdpy11::Sregion
+{
+  private:
+    std::vector<double> esize, h, weight;
+    fwdpp::gsl_ran_discrete_t_ptr sh_lookup;
+
+    fwdpp::gsl_ran_discrete_t_ptr
+    init_lookup()
+    {
+        for (auto wi : weight)
+            {
+                if (wi < 0)
+                    {
+                        throw std::invalid_argument("all weights must be >= 0");
+                    }
+                if (!std::isfinite(wi))
+                    {
+                        throw std::invalid_argument("all weights must be finite");
+                    }
+            }
+        return fwdpp::gsl_ran_discrete_t_ptr(
+            gsl_ran_discrete_preproc(weight.size(), weight.data()));
+    }
+
+  public:
+    DiscreteDESD(const fwdpy11::Region& r, const double sc, std::vector<double> esize_,
+                 std::vector<double> h_, std::vector<double> w_)
+        : fwdpy11::Sregion(r, sc, 1), esize{std::move(esize_)}, h{std::move(h_)},
+          weight{std::move(w_)}, sh_lookup{init_lookup()}
+    {
+        if (esize.size() != h.size() || esize.size() != weight.size()
+            || h.size() != weight.size())
+            {
+                throw std::invalid_argument("all arrays must be equal-length");
+            }
+    }
+
+    std::unique_ptr<fwdpy11::Sregion>
+    clone() const override
+    {
+        return std::make_unique<DiscreteDESD>(this->region, this->scaling, this->esize,
+                                              this->h, this->weight);
+    }
+
+    std::uint32_t
+    operator()(fwdpp::flagged_mutation_queue& recycling_bin,
+               std::vector<fwdpy11::Mutation>& mutations,
+               std::unordered_multimap<double, std::uint32_t>& lookup_table,
+               const std::uint32_t generation,
+               const fwdpy11::GSLrng_t& rng) const override
+    {
+        auto idx = gsl_ran_discrete(rng.get(), sh_lookup.get());
+        return fwdpy11::infsites_Mutation(
+            recycling_bin, mutations, lookup_table, false, generation,
+            [this, &rng]() { return region(rng); },
+            [this, idx]() { return esize[idx] / scaling; },
+            [this, idx]() { return h[idx]; }, this->label());
+    }
+
+    double
+    from_mvnorm(const double /*deviate*/, const double P) const override
+    {
+        throw std::runtime_error("not implemented yet");
+        return 1.;
+    }
+
+    std::vector<double>
+    get_dominance() const override
+    {
+        return {h};
+    }
+};
+
+void
+init_DiscreteDESD(py::module& m)
+{
+    py::class_<DiscreteDESD, fwdpy11::Sregion>(m, "_ll_DiscreteDESD")
+        .def(py::init([](double beg, double end, double weight,
+                         std::vector<std::tuple<double, double, double>> joint_dist,
+                         bool coupled, std::uint16_t label, double scaling) {
+                 if (joint_dist.empty())
+                     {
+                         throw std::invalid_argument("empty input list");
+                     }
+                 std::vector<double> esizes, dominance, weights;
+                 for (auto&& t : joint_dist)
+                     {
+                         auto e = std::get<0>(t);
+                         auto h = std::get<1>(t);
+                         auto w = std::get<2>(t);
+                         if (!std::isfinite(e))
+                             {
+                                 throw std::invalid_argument(
+                                     "non-finite effect size input");
+                             }
+                         if (!std::isfinite(h))
+                             {
+                                 throw std::invalid_argument(
+                                     "non-finite dominance input");
+                             }
+                         if (!std::isfinite(w))
+                             {
+                                 throw std::invalid_argument("non-finite weight input");
+                             }
+                         if (w < 0.0)
+                             {
+                                 throw std::invalid_argument("weights must be >= 0.0");
+                             }
+                         esizes.push_back(e);
+                         dominance.push_back(h);
+                         weights.push_back(w);
+                     }
+                 return DiscreteDESD(fwdpy11::Region(beg, end, weight, coupled, label),
+                                     scaling, std::move(esizes), std::move(dominance),
+                                     std::move(weights));
+             }),
+             py::arg("beg"), py::arg("end"), py::arg("weight"), py::arg("joint_dist"),
+             py::arg("coupled"), py::arg("label"), py::arg("scaling"));
+}
+

fwdpy11/src/regions/init.cc

@@ -14,6 +14,7 @@ void init_RecombinationRegions(py::module &);
 void init_MultivariateGaussianEffects(py::module &);
 void init_mvDES(py::module &);
 void init_LogNormalS(py::module &);
+void init_DiscreteDESD(py::module&);
 
 void
 initialize_regions(py::module &m)
@@ -30,4 +31,5 @@ initialize_regions(py::module &m)
     init_MultivariateGaussianEffects(m);
     init_mvDES(m);
     init_LogNormalS(m);
+    init_DiscreteDESD(m);
 }

tests/test_regions.py

@@ -408,6 +408,51 @@ def test_call(self):
         self.assertTrue(np.all(m.esizes <= 0.0))
 
 
+class TestDiscreteDESD(unittest.TestCase):
+    def test_init(self):
+        jdist = [(0.1, 1, 1), (0.2, 0.5, 1)]
+        jd = fwdpy11.DiscreteDESD(0, 1, 2, jdist)
+
+        self.assertEqual(jd.beg, 0)
+        self.assertEqual(jd.end, 1)
+        self.assertEqual(jd.weight, 2)
+        self.assertEqual(jd.joint_dist, jdist)
+
+    def test_pickle(self):
+        jdist = [(0.1, 1, 1), (0.2, 0.5, 1)]
+        jd = fwdpy11.DiscreteDESD(0, 1, 2, jdist)
+
+        pp = pickle.dumps(jd, -1)
+        up = pickle.loads(pp)
+
+        self.assertEqual(jd.beg, up.beg)
+        self.assertEqual(jd.end, up.end)
+        self.assertEqual(jd.weight, up.weight)
+        self.assertEqual(jd.coupled, up.coupled)
+        self.assertEqual(jd.scaling, up.scaling)
+        self.assertEqual(jd.joint_dist, up.joint_dist)
+
+    def test_bad_init_1(self):
+        jdist = [(0.1, 1, np.nan), (0.2, 0.5, 1)]
+        with self.assertRaises(ValueError):
+            fwdpy11.DiscreteDESD(0, 1, 2, jdist)
+
+    def test_bad_init_2(self):
+        jdist = [(0.1, 1, -1), (0.2, 0.5, 1)]
+        with self.assertRaises(ValueError):
+            fwdpy11.DiscreteDESD(0, 1, 2, jdist)
+
+    def test_bad_init_3(self):
+        jdist = [(0.1, np.nan, 1), (0.2, 0.5, 1)]
+        with self.assertRaises(ValueError):
+            fwdpy11.DiscreteDESD(0, 1, 2, jdist)
+
+    def test_bad_init_4(self):
+        jdist = [(np.nan, 1, 1), (0.2, 0.5, 1)]
+        with self.assertRaises(ValueError):
+            fwdpy11.DiscreteDESD(0, 1, 2, jdist)
+
+
 class testPoissonInterval(unittest.TestCase):
     @classmethod
     def setUp(self):