Add vignettes on varying fitness effects across demes

doc/_toc.yml

@@ -32,6 +32,13 @@ parts:
   - file: short_vignettes/demography_vignettes_intro
   - file: short_vignettes/demes_vignette
   - file: short_vignettes/demography_debugger_vignette
+- caption: Short vignettes - varying fitness conditions across demes
+  chapters:
+  - file: short_vignettes/varying_fitness_conditions_across_demes_intro
+  - file: short_vignettes/multivariate_gaussian_effect_sizes_across_demes
+  - file: short_vignettes/multivariate_lognormal_effect_sizes
+  - file: short_vignettes/custom_multivarate_distributions
+  - file: short_vignettes/special_considerations_demes_models
 - caption: Short vignettes - conditional models
   chapters:
   - file: short_vignettes/trackmutationfates

doc/short_vignettes/custom_multivarate_distributions.md

@@ -0,0 +1,131 @@
+---
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+kernelspec:
+  display_name: Python 3
+  language: python
+  name: python3
+---
+
+(custom_multivariate_distributions_vignette)=
+
+# "Custom" multivariate distributions
+
+```{code-cell} python
+---
+tags: ['hide-input']
+---
+import demes
+import demesdraw
+import fwdpy11
+import numpy as np
+```
+
+The previous two sections cover cases where the methods for generating
+deviates from a multivariate distribution are straightforward and agreed
+upon.
+
+In order to simulate multivariate distributions of effect sizes based on
+{class}`fwdpy11.Sregion` types, we follow a fairly intuitive approach
+described in {cite}`Xue-Kun_Song2000-qn`.  Briefly, the multivariate Gaussian kernel is
+used to produce deviates.  Then, the quantiles from the cummulative distribution
+of each marginal Gaussian are used to generate a deviate from the desired output distribution of interest.
+
+For a simulation with `n` populations we need:
+
+* A {class}`list` of `n` {class}`fwdpy11.Sregion` objects
+* An array of `n` means for the multivariate Gaussian
+* An `n-by-n` covariance matrix for the multivariate
+  Gaussian
+
+We will use the same demographic model as in the previous vignettes.
+The details can be viewed by expanding the "click to show" button below".
+
+```{code-cell}python
+---
+tags: ['hide-input']
+---
+yaml = """
+description: Island model forever
+time_units: generations
+demes:
+  - name: A
+    epochs:
+     - start_size: 100
+  - name: B
+    epochs:
+     - start_size: 100
+migrations:
+  - demes: [A, B]
+    rate: 0.10
+"""
+g = demes.loads(yaml)
+model = fwdpy11.discrete_demography.from_demes(g, burnin=1)
+initial_sizes = [v for v in model.metadata["initial_sizes"].values()]
+pdict = {
+    "nregions": [],
+    "recregions": [],
+    "sregions": None, # Will get filled in below
+    "rates": (0, 5e-3, None),
+    "demography": model,
+    "simlen": model.metadata["total_simulation_length"],
+    "gvalue": fwdpy11.Multiplicative(ndemes=2, scaling=2),
+}
+rng = fwdpy11.GSLrng(123512)
+```
+
+The following generates exponentially distributed effect sizes in each deme
+with a high correlation across demes:
+
+```{code-cell} python
+mvdes = fwdpy11.mvDES(
+    [fwdpy11.ExpS(0, 1, 1, -0.5)] * 2,
+    np.zeros(2),
+    np.matrix([1, 0.9, 0.9, 1]).reshape((2, 2)),
+)
+pdict["sregions"] = [mvdes]
+params = fwdpy11.ModelParams(**pdict)
+pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
+fwdpy11.evolvets(rng, pop, params, 10)
+for i in pop.tables.mutations:
+    print(pop.mutations[i.key].esizes)
+```
+
+We can mix and match our distributions.  Here, the distribution of effect
+sizes in deme 0 is exponential and the distribution in deme 1 is gamma.  The
+two distributions have means with opposite signs and the magnitudes of the
+marginal deviates negatively covary:
+
+```{code-cell} python
+mvdes = fwdpy11.mvDES(
+    [fwdpy11.ExpS(0, 1, 1, -0.5), fwdpy11.GammaS(0, 1, 1, mean=0.1, shape_parameter=1)],
+    np.zeros(2),
+    np.matrix([1, -0.9, -0.9, 1]).reshape((2, 2)),
+)
+pdict["sregions"] = [mvdes]
+params = fwdpy11.ModelParams(**pdict)
+pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
+fwdpy11.evolvets(rng, pop, params, 10)
+for i in pop.tables.mutations:
+    print(pop.mutations[i.key].esizes)
+```
+
+The type {class}`fwdpy11.ConstantS` has intuitive behavior:
+
+```{code-cell} python
+mvdes = fwdpy11.mvDES(
+    [fwdpy11.ExpS(0, 1, 1, -0.5), fwdpy11.ConstantS(0, 1, 1, -0.1)],
+    np.zeros(2),
+    np.matrix([1, -0.9, -0.9, 1]).reshape((2, 2)),
+)
+pdict["sregions"] = [mvdes]
+params = fwdpy11.ModelParams(**pdict)
+pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
+rng = fwdpy11.GSLrng(1010)
+fwdpy11.evolvets(rng, pop, params, 10)
+for i in pop.tables.mutations:
+    print(pop.mutations[i.key].esizes)
+```

doc/short_vignettes/genetic_values_multi_deme_models.md

@@ -0,0 +1,23 @@
+---
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+kernelspec:
+  display_name: Python 3
+  language: python
+  name: python3
+---
+
+(genetic_values_multi_deme_models)=
+
+# Setting up genetic value objects.
+
+NOTE: this will get too long -- split!!
+
+## Case 1: constant effect size, different fitness effects.
+
+## Case 2: different effect sizes, same fitness effects.
+
+## Case 3: different effect sizes, different fitness effects.

doc/short_vignettes/multivariate_gaussian_effect_sizes_across_demes.md

@@ -0,0 +1,146 @@
+---
+jupytext:
+  formats: md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+kernelspec:
+  display_name: Python 3
+  language: python
+  name: python3
+---
+
+(multivariate_gaussian_effect_sizes_across_demes_vignette)=
+
+# The multivariate Gaussian distribution
+
+```{code-cell}python
+---
+tags: ['hide-input']
+---
+import demes
+import demesdraw
+import fwdpy11
+import numpy as np
+```
+
+
+We may model Gaussian effect sizes using the existing {class}`fwdpy11.MultivariateGaussianEffects`
+in conjunction with {class}`fwdpy11.mvDES`.
+
+At this time, it is probably best to look at an example. The following code models Gaussian stabilizing
+selection on a quantitative trait.  The effects sizes within each deme are themselves given by Gaussian
+distributions and there is no correlation in the effect size in the two demes.
+
+We will simulate a demographic model of migration happening into the infinite past of two equal-sized demes:
+
+```{code-cell}python
+---
+tags: ['hide-input']
+---
+yaml = """
+description: Island model forever
+time_units: generations
+demes:
+  - name: ancestor
+    epochs:
+     - start_size: 100
+       end_time: 100
+  - name: A
+    ancestors: [ancestor]
+    epochs:
+     - start_size: 100
+  - name: B
+    ancestors: [ancestor]
+    epochs:
+     - start_size: 100
+migrations:
+  - demes: [A, B]
+    rate: 0.10
+"""
+g = demes.loads(yaml)
+model = fwdpy11.discrete_demography.from_demes(g, burnin=1)
+demesdraw.tubes(g);
+```
+
+Let's set up a dictionary to hold the parameters of our model:
+
+```{code-cell} python
+pdict = {
+    "nregions": [],
+    "recregions": [],
+    "sregions": [
+        fwdpy11.mvDES(
+            fwdpy11.MultivariateGaussianEffects(0, 1, 1, np.identity(3)), np.zeros(3)
+        )
+    ],
+    "rates": (0, 2.5e-3, None),
+    "demography": model,
+    "simlen": model.metadata["total_simulation_length"],
+    "gvalue": fwdpy11.Additive(
+        ndemes=3, scaling=2, gvalue_to_fitness=fwdpy11.GSS(optimum=0.0, VS=10.0)
+    ),
+}
+```
+
+Most of the above is standard.  Let's dissect the new bits:
+
+* An instance of {class}`fwdpy11.mvDES` is our only region with selected mutations.
+* This instance holds an instance of {class}`fwdpy11.MultivariateGaussianEffects`
+  that puts mutations on the interval {math}`[0, 1)` with weight 1 and an identity
+  matrix specifies the correlation in effect sizes between demes 0 and 1.  The
+  identity matrix has the value zero for all off-diagonal elements, meaning
+  no covariance in effect sizes across demes.
+* The final constructor argument specifies the mean of each marginal Gaussian
+  distribution. The means are both zero.
+* Our genetic value type accepts an `ndemes` parameter, telling it that it has
+  to look for deme-specific effect sizes.  This value must be set to the maximum
+  number of demes that will exist during a simulation.
+
+Let's evolve the model now:
+
+```{code-cell} python
+params = fwdpy11.ModelParams(**pdict)
+# TODO: update this once we have a function to pull the sizes
+# automatically from demes-derived models:
+initial_sizes = [v for v in model.metadata["initial_sizes"].values()]
+pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
+rng = fwdpy11.GSLrng(1010)
+fwdpy11.evolvets(rng, pop, params, 10)
+```
+
+Let's extract the effect sizes from each deme:
+
+```{code-cell} python
+assert len(pop.tables.mutations) > 0
+for i in pop.tables.mutations:
+    print(pop.mutations[i.key].esizes)
+```
+
+Let's look at another example where effect sizes covary negatively across demes and raise the mutation rate a bit:
+
+```{code-cell} python
+# Effect sizes across demes will
+# have a correlation coefficient of r=1/2
+cor_matrix = np.array([-0.5]*9).reshape(3,3)
+np.fill_diagonal(cor_matrix, np.array([1.0]*3))
+
+# Get our covariance matrix
+sd = np.array([0.1]*3)
+D = np.identity(3)
+np.fill_diagonal(D, sd)
+vcv_matrix = np.matmul(np.matmul(D, cor_matrix), D)
+pdict["sregions"] = [
+    fwdpy11.mvDES(fwdpy11.MultivariateGaussianEffects(0, 1, 1, vcv_matrix), np.zeros(3))
+]
+params = fwdpy11.ModelParams(**pdict)
+# TODO: update this once we have a function to pull the sizes
+# automatically from demes-derived models:
+initial_sizes = [v for v in model.metadata["initial_sizes"].values()]
+pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
+fwdpy11.evolvets(rng, pop, params, 10)
+for i in pop.tables.mutations:
+    print(pop.mutations[i.key].esizes)
+```
+
+Now we see that the effect sizes often differ in sign between the two demes.