Repository: popsim-consortium/stdpopsim (language: python)

Issue:
add `uniform` DFE
In #1470 @klohmueller says

> I would also suggest adding a “discrete DFE” as well. We could use one from Table 2 of Kim et al (2017). I will give the parameters here that we inferred from the EUR 1000 Genomes data with mu=1.5e-8. Note, ranges of |s| correspond to: [0, 10e−5), [10e−5, 10e−4), [10e−4, 10e−3), and [10e−3, 10e−2), respectively. The mass in the [10e−2, 1] range is the complement of the total mass of the four aforementioned categories. We then inferred the proportions of mutations in each bin. We inferred the proportions to be: 0.286, 0.099, 0.222, 0.305.

To do this, we'd just need to add a single "uniform" distribution type to the [list of distribution types](https://github.com/popsim-consortium/stdpopsim/blob/main/stdpopsim/dfe.py#L21) (and then the implementation would use four of these). However, @RyanGutenkunst asks

> For this parameterization, 0.088 weight is for the [1e-2, 1] range. Is that uniform on the linear scale or the log scale? In other words, is half the mass assumed to reside between [1e-2, 0.505] or between [1e-2, 1e-1]?

This suggests *also* adding a "log uniform" type? (which would, I think, be a truncated exponential) And I guess also a "negative log uniform"?
add `uniform` DFE
In #1470 @klohmueller says

> I would also suggest adding a “discrete DFE” as well. We could use one from Table 2 of Kim et al (2017). I will give the parameters here that we inferred from the EUR 1000 Genomes data with mu=1.5e-8. Note, ranges of |s| correspond to: [0, 10e−5), [10e−5, 10e−4), [10e−4, 10e−3), and [10e−3, 10e−2), respectively. The mass in the [10e−2, 1] range is the complement of the total mass of the four aforementioned categories. We then inferred the proportions of mutations in each bin. We inferred the proportions to be: 0.286, 0.099, 0.222, 0.305.

To do this, we'd just need to add a single "uniform" distribution type to the [list of distribution types](https://github.com/popsim-consortium/stdpopsim/blob/main/stdpopsim/dfe.py#L21) (and then the implementation would use four of these). However, @RyanGutenkunst asks

> For this parameterization, 0.088 weight is for the [1e-2, 1] range. Is that uniform on the linear scale or the log scale? In other words, is half the mass assumed to reside between [1e-2, 0.505] or between [1e-2, 1e-1]?

This suggests *also* adding a "log uniform" type? (which would, I think, be a truncated exponential) And I guess also a "negative log uniform"?


Discussion / hints:

Fix patch:
diff --git a/stdpopsim/dfe.py b/stdpopsim/dfe.py
index 78ac8e7c1..d2858671e 100644
--- a/stdpopsim/dfe.py
+++ b/stdpopsim/dfe.py
@@ -23,6 +23,7 @@ class MutationType(object):
     - ``g``: gamma, two parameters (mean, shape)
     - ``n``: normal, two parameters (mean, SD)
     - ``w``: Weibull, two parameters (scale, shape)
+    - ``u``: Uniform, two parameters (min, max)
     - ``lp``: positive logNormal, two parameters (mean and sd on log scale; see rlnorm)
     - ``ln``: negative logNormal, two parameters (mean and sd on log scale; see rlnorm)
 
@@ -143,6 +144,19 @@ def __attrs_post_init__(self):
                 f"return {sign}rlnorm(1, {logmean} + log(Q), {logsd});"
             ]
             self.distribution_type = "s"
+        elif self.distribution_type == "u":
+            # Uniform
+            if (
+                len(self.distribution_args) != 2
+                or self.distribution_args[0] > self.distribution_args[1]
+            ):
+                raise ValueError(
+                    "Uniformly-distributed sel. coefs. (distribution_type='u') "
+                    "use a (min, max) parameterisation, with min <= max."
+                )
+            umin, umax = self.distribution_args
+            self.distribution_args = [f"return runif(1, Q * {umin}, Q * {umax});"]
+            self.distribution_type = "s"
         else:
             raise ValueError(
                 f"{self.distribution_type} is not a supported distribution type."


Test patch:
diff --git a/tests/test_dfes.py b/tests/test_dfes.py
index 54b16b305..2c62f6fc3 100644
--- a/tests/test_dfes.py
+++ b/tests/test_dfes.py
@@ -32,9 +32,9 @@ def test_Q_scaled_index(self):
             "n": ([0.5, 1], [0, 1]),
             "ln": ([0.5, 1], []),
             "lp": ([0.5, 1], []),
+            "u": ([-0.5, 1], []),
         }
         for t in mut_params:
-            print(f"{t}\t{mut_params[t]}")
             if t == "f":
                 mt = dfe.MutationType(
                     distribution_type=t,
@@ -157,6 +157,14 @@ def test_create_bad_mutation_type_message(self):
                 distribution_args=[1, -2],
             )
 
+        # Uniformly-distributed selection coefficients
+        for bad_args in ([0], [1, 2, 3], [3, -2]):
+            with pytest.raises(ValueError, match="use a .min, max"):
+                dfe.MutationType(
+                    distribution_type="u",
+                    distribution_args=bad_args,
+                )
+
         # Lognormally-distributed selection coefficients
         for dt in ["lp", "ln"]:
             with pytest.raises(
@@ -204,13 +212,14 @@ def test_mutation_types(self):
             "e": ([0.1], [10], [5000], [0]),
             "n": ([-0.1, 0.2], [0.1, 0.1], [50, 50]),
             "w": ([0.1, 0.2], [0.1, 0.1], [50, 50]),
+            "u": ([0.1, 0.2], [0.1, 0.1], [-5, 50]),
             "lp": ([-0.1, 0.2], [0.1, 0.1], [50, 50]),
             "ln": ([-0.1, 0.2], [0.1, 0.1], [50, 50]),
         }
         for t in mut_params:
             for p in mut_params[t]:
                 mt = dfe.MutationType(distribution_type=t, distribution_args=p)
-                if t in ("lp", "ln"):
+                if t in ("lp", "ln", "u"):
                     assert mt.distribution_type == "s"
                 else:
                     assert mt.distribution_type == t
@@ -225,6 +234,7 @@ def test_bad_mutation_types(self):
             "e": ([], [0, 1], [0.1, 0.4, 0.5], [np.inf]),
             "n": ([], [0.1, -1], [0.1, 0.4, 0.5], [0.1], [0.3, np.inf]),
             "w": ([], [-0.1, 1], [0.1, -1], [0.1, 0.4, 0.5], [0.1], [np.inf, 2.3]),
+            "u": ([], [-0.1, -0.5], [0.1, 0.4, 0.5], [0.1], [2.3, np.inf]),
             "lp": ([], [0.1, -1], [0.1, 0.4, 0.5], [0.1], [0.1, np.inf]),
             "ln": ([], [0.1, -1], [0.1, 0.4, 0.5], [0.1], [0.1, np.inf]),
         }