Add mixed distribution support to TrueMeasure

docs/api/true_measures.md

@@ -24,6 +24,10 @@ jupyter:
 
 ::: qmcpy.true_measure.product_measure.ProductMeasure
 
+## `Mixture`
+
+::: qmcpy.true_measure.mixture.Mixture
+
 ## `StudentT`
 
 ::: qmcpy.true_measure.student_t.StudentT

makefile

@@ -633,6 +633,16 @@ uml:
 #
 # Use `mkdocs serve` to run a local server. The webpages are stored in a temporary folder and will be deleted when the server is stopped.
 ##########################################################
+# mkdocs is a `docs` extra, installed alongside $(PYTHON) in the qmcpy env.
+# Prefer that colocated binary over a bare PATH lookup: an older `pip install
+# --user` shim earlier on PATH (e.g. left over from a Python version bump
+# that removed the interpreter its shebang points at) can shadow the correct
+# one and fail with "bad interpreter" (exit 126) instead of a clean "not
+# found". Mirror-image of check_pydoclint_dependency's PATH-first order,
+# which instead assumes pydoclint may live in a separate, lighter test-only
+# env rather than this one.
+MKDOCS ?= $(shell test -x "$(dir $(PYTHON))mkdocs" && echo "$(dir $(PYTHON))mkdocs" || command -v mkdocs 2>/dev/null || echo mkdocs)
+
 copydocs:  # mkdocs only looks for content in the docs/ folder, so we have to copy it there
 	@rm -rf docs/paper docs/demos
 	@cp README.md docs/README.md
@@ -670,19 +680,18 @@ runmkdocserve:
 		PORT=$$((PORT+1)); \
 	done; \
 	echo "Starting mkdocs on http://127.0.0.1:$$PORT"; \
-	NO_MKDOCS_2_WARNING=1 JUPYTER_PLATFORM_DIRS=1 mkdocs serve -a 127.0.0.1:$$PORT
-	NO_MKDOCS_2_WARNING=1 JUPYTER_PLATFORM_DIRS=1 mkdocs serve -a 127.0.0.1:$$PORT
+	NO_MKDOCS_2_WARNING=1 JUPYTER_PLATFORM_DIRS=1 $(MKDOCS) serve -a 127.0.0.1:$$PORT
 
 doc: uml copydocs runmkdocserve
 
 docnouml: copydocs runmkdocserve
 
 check_links: copydocs  # internal links + anchors only; fast, no network, safe for CI
-	@NO_MKDOCS_2_WARNING=1 mkdocs build -q -d site
+	@NO_MKDOCS_2_WARNING=1 $(MKDOCS) build -q -d site
 	@$(PYTHON) scripts/check_links.py site
 
 check_links_external: copydocs  # also checks http/https links; slow and network-flaky, run locally
-	@NO_MKDOCS_2_WARNING=1 mkdocs build -q -d site
+	@NO_MKDOCS_2_WARNING=1 $(MKDOCS) build -q -d site
 	@$(PYTHON) scripts/check_links.py site --external
 
 # The targets above check links inside the new site; these check the other

mkdocs.yml

@@ -77,6 +77,7 @@ nav:
         - Some True Measures: demos/some_true_measures.ipynb
         - SciPyWrapper dependence and Custom distributions: demos/scipywrapper_dependence_custom/scipywrapper_demo.ipynb
         - ProductMeasure: demos/product_measure.ipynb
+        - Mixed Distributions: demos/mixed_distributions.ipynb
         - Acceptance-Rejection Sampling: demos/acceptance_rejection.ipynb
         - Copula TrueMeasure Examples: demos/copula_examples.ipynb
   - For Developers:

qmcpy/discrete_distribution/digital_net_b2/digital_net_b2.py

@@ -704,7 +704,17 @@ def _gen_samples(self, n_min, n_max, return_binary, warn):
                 if warn:
                     warnings.warn("DigitalNetB2 in graycode order recommends n_min and n_max be 0 or powers of 2 at which the digital net achieves superior uniformity properties")
             elif self.order == "RADICAL INVERSE":
-                raise ParameterError("DigitalNetB2 in radical inverse order requires n_min and n_max be 0 or powers of 2")
+                message = (
+                    "DigitalNetB2 in radical inverse order cannot generate "
+                    f"n={n_max - n_min} samples with n_min={n_min} and n_max={n_max}: "
+                    "both endpoints must be 0 or powers of 2."
+                )
+                if n_min == 0:
+                    message += (
+                        " When starting at index 0, use a sample count "
+                        "of 0, 1, 2, 4, 8, ... ."
+                    )
+                raise ParameterError(message)
             else:
                 raise ValueError("invalid digital net order")
         r = np.uint64(self.replications)

qmcpy/integrand/abstract_integrand.py

@@ -256,12 +256,13 @@ def f(self, x: np.ndarray, *args: tuple, **kwargs: dict) -> np.ndarray:
             raise AssertionError
         # function evaluation with chain rule
         i = (None,) * d_indv_ndim + (...,)
+        transformed_shape = batch_shape + (self.true_measure.d,)
         if self.true_measure == self.true_measure.transform:
             # jacobian*weight/pdf will cancel so f(x) = g(\Psi(x))
             xtf = self.true_measure._jacobian_transform_r(
                 xp, return_weights=False
             )  # get transformed samples, equivalent to self.true_measure._transform_r(x)
-            if not (xtf.shape == xp.shape):
+            if not (xtf.shape == transformed_shape):
                 raise AssertionError
             y = self._g(xtf, *args, **kwargs)
         else:  # using importance sampling --> need to compute pdf, jacobian(s), and weight explicitly
@@ -271,7 +272,7 @@ def f(self, x: np.ndarray, *args: tuple, **kwargs: dict) -> np.ndarray:
             xtf, jacobians = self.true_measure.transform._jacobian_transform_r(
                 xp, return_weights=True
             )  # compute recursive transform+jacobian
-            if not (xtf.shape == xp.shape):
+            if not (xtf.shape == transformed_shape):
                 raise AssertionError
             if not (jacobians.shape == batch_shape):
                 raise AssertionError

qmcpy/integrand/sensitivity_indices.py

@@ -3,6 +3,7 @@
 from .keister import Keister
 from .box_integral import BoxIntegral
 from ..discrete_distribution import DigitalNetB2
+from ..util import DimensionError
 import numpy as np
 from itertools import combinations
 import scipy.special
@@ -109,7 +110,7 @@ def __init__(self, integrand: AbstractIntegrand, indices: Union[str, np.ndarray]
 
         Args:
             integrand (AbstractIntegrand): Integrand to find sensitivity
-                indices of.
+                indices of. Its true measure must preserve the sampler dimension.
             indices (Union[str, np.ndarray]): Bool array with shape $(\dots,d)$ where each
                 length $d$ vector item indicates which dimensions are active in
                 the subset.
@@ -119,6 +120,10 @@ def __init__(self, integrand: AbstractIntegrand, indices: Union[str, np.ndarray]
         """
         self.parameters = ["indices"]
         self.integrand = integrand
+        if integrand.discrete_distrib.d != integrand.d:
+            raise DimensionError(
+                "SensitivityIndices requires a dimension-preserving true measure."
+            )
         self.dtilde = self.integrand.d
         if not (self.dtilde > 1):
             raise AssertionError("SensitivityIndices does not make sense for d=1")

qmcpy/stopping_criterion/abstract_cub_bayes_ld_g.py

@@ -390,7 +390,7 @@ def integrate(self, resume: Union[None, Data] = None) -> tuple:
             data.n_min = 0
             data.n_max = self.n_init
             data.solution_indv = np.tile(np.nan, self.integrand.d_indv)
-            data.xfull = np.empty((0, self.integrand.d))
+            data.xfull = np.empty((0, self.discrete_distrib.d))
             data.yfull = np.empty(self.integrand.d_indv + (0,))
             data.bounds_half_width = np.tile(np.inf, self.integrand.d_indv)
             data.muhat = np.tile(np.nan, self.integrand.d_indv)
@@ -497,7 +497,7 @@ def _validate_resume(self, data):
         )
         n_total = int(data.n_total)
         output_shape = self.integrand.d_indv + (n_total,)
-        self._validate_resume_shape("xfull", data.xfull, (n_total, self.integrand.d))
+        self._validate_resume_shape("xfull", data.xfull, (n_total, self.discrete_distrib.d))
         self._validate_resume_shape("yfull", data.yfull, output_shape)
         self._validate_resume_shape("_ytildefull", data._ytildefull, output_shape)
         self._validate_resume_shape("n", data.n, self.integrand.d_indv)

qmcpy/stopping_criterion/abstract_cub_qmc_ld_g.py

@@ -211,7 +211,7 @@ def _validate_resume(self, data):
         )
         n_total = int(data.n_total)
         output_shape = self.integrand.d_indv + (n_total,)
-        self._validate_resume_shape("xfull", data.xfull, (n_total, self.integrand.d))
+        self._validate_resume_shape("xfull", data.xfull, (n_total, self.discrete_distrib.d))
         self._validate_resume_shape("yfull", data.yfull, output_shape)
         self._validate_resume_shape("_ytildefull", data._ytildefull, output_shape)
         self._validate_resume_shape("_kappanumap", data._kappanumap, output_shape)
@@ -280,7 +280,7 @@ def integrate(self, resume: Union[None, Data] = None) -> tuple:
             data.n_min = 0
             data.n_max = self.n_init
             data.solution_indv = np.tile(np.nan, self.integrand.d_indv)
-            data.xfull = np.empty((0, self.integrand.d))
+            data.xfull = np.empty((0, self.discrete_distrib.d))
             data.yfull = np.empty(self.integrand.d_indv + (0,))
             if self.ncv > 0:
                 data.ycvfull = np.empty(self.integrand.d_indv + (self.ncv, 0))

qmcpy/stopping_criterion/cub_mc_clt.py

@@ -230,7 +230,7 @@ def integrate(self, resume: Union[None, Data] = None) -> tuple:
         if resume is not None:
             raise ParameterError("CubMCCLT does not support resume.")
         data = Data(parameters=["solution", "bound_low", "bound_high", "bound_diff", "n_total", "time_integrate"])
-        data.xfull = np.empty((0, self.integrand.d))
+        data.xfull = np.empty((0, self.discrete_distrib.d))
         data.yfull = np.empty(0)
         if self.ncv > 0:
             data.ycvfull = np.empty((self.ncv, 0))

qmcpy/stopping_criterion/cub_mc_clt_vec.py

@@ -302,7 +302,7 @@ def integrate(self, resume: Union[None, Data] = None) -> tuple:
             data.n_min = 0
             data.n_max = self.n_init
             data.solution_indv = np.tile(np.nan, self.integrand.d_indv)
-            data.xfull = np.empty((0, self.integrand.d))
+            data.xfull = np.empty((0, self.discrete_distrib.d))
             data.yfull = np.empty(self.integrand.d_indv + (0,))
         first_resume_iter = resume is not None
         while True:

qmcpy/stopping_criterion/cub_mc_g.py

@@ -347,7 +347,7 @@ def integrate(self, resume: Union[None, Data] = None) -> tuple:
         if resume is not None:
             raise ParameterError("CubMCG does not support resume.")
         data = Data(parameters=["solution", "bound_low", "bound_high", "bound_diff", "n_total", "time_integrate"])
-        data.xfull = np.empty((0, self.integrand.d))
+        data.xfull = np.empty((0, self.discrete_distrib.d))
         data.yfull = np.empty(0)
         if self.ncv > 0:
             data.ycvfull = np.empty((self.ncv, 0))

qmcpy/stopping_criterion/cub_qmc_rep_student_t.py

@@ -334,7 +334,7 @@ def integrate(self, resume: Union[None, Data] = None) -> tuple:
             data.n_min = 0
             data.n_max = self.n_init
             data.solution_indv = np.tile(np.nan, self.integrand.d_indv)
-            data.xfull = np.empty((self.discrete_distrib.replications, 0, self.integrand.d))
+            data.xfull = np.empty((self.discrete_distrib.replications, 0, self.discrete_distrib.d))
             data.yfull = np.empty(self.integrand.d_indv + (self.discrete_distrib.replications, 0))
             data._ysums = np.zeros(self.integrand.d_indv + (self.discrete_distrib.replications,), dtype=float)
         first_resume_iter = resume is not None
@@ -410,10 +410,10 @@ def _validate_resume(self, data):
             raise ParameterError(
                 "resume data must include at least n_init samples per replication."
             )
-        if data.xfull.shape != (replications, n_rep_max, self.integrand.d):
+        if data.xfull.shape != (replications, n_rep_max, self.discrete_distrib.d):
             raise ParameterError(
                 "resume data xfull shape must be (%d, %d, %d); got %s."
-                % (replications, n_rep_max, self.integrand.d, data.xfull.shape)
+                % (replications, n_rep_max, self.discrete_distrib.d, data.xfull.shape)
             )
         expected_y_shape = self.integrand.d_indv + (replications, n_rep_max)
         if np.shape(data.yfull) != expected_y_shape:

qmcpy/stopping_criterion/pf_gp_ci.py

@@ -308,7 +308,7 @@ def __init__(
         self.integrand = integrand
         self.true_measure = self.integrand.true_measure
         self.discrete_distrib = self.integrand.discrete_distrib
-        self.d = self.integrand.d
+        self.d = self.discrete_distrib.d
         self.sampler = self.d
         self.failure_threshold = failure_threshold
         self.failure_above_threshold = failure_above_threshold

qmcpy/true_measure/__init__.py

@@ -14,6 +14,7 @@
 from .johnsons_su import JohnsonsSU
 from .scipy_wrapper import SciPyWrapper
 from .matern_gp import MaternGP
+from .mixture import Mixture
 from .student_t import StudentT
 from .student_t_copula import StudentTCopula
 from .uniform_triangle import UniformTriangle