Support domain inclusion in chained TrueMeasure transformations

mkdocs.yml

@@ -75,6 +75,7 @@ nav:
       - Importance Sampling with True Measures:
         - Statistics for True Measures: demos/statistics_for_TrueMeasure.ipynb
         - Some True Measures: demos/some_true_measures.ipynb
+        - TrueMeasure Domain Inclusion: demos/true_measure_domain_inclusion.ipynb
         - SciPyWrapper dependence and Custom distributions: demos/scipywrapper_dependence_custom/scipywrapper_demo.ipynb
         - ProductMeasure: demos/product_measure.ipynb
         - Acceptance-Rejection Sampling: demos/acceptance_rejection.ipynb

qmcpy/discrete_distribution/dummy_sampler.py

@@ -30,6 +30,16 @@ class DummySampler(AbstractLDDiscreteDistribution):
     """
 
     def __init__(self, dimension=1, replications=None, seed=None, warn=True):
+        r"""
+        Args:
+            dimension (Union[int, list, tuple, np.ndarray]): Dimension of the placeholder sampler. A list, tuple, or array specifies unique coordinate indices. Defaults to `1`.
+            replications (Union[None, int]): Replication metadata preserved when spawning placeholders. `None` records no explicit replication axis. Defaults to `None`.
+            seed (Union[None, int, np.random.SeedSequence]): Seed used to initialize the sampler state and spawn child samplers. Defaults to `None`.
+            warn (bool): Compatibility argument matching other discrete-distribution constructors. It is ignored because `DummySampler` cannot generate samples. Defaults to `True`.
+
+        Raises:
+            ParameterError: If an array-like `dimension` is not one-dimensional with unique entries, if it exceeds the dimension limit, or if `replications` is negative.
+        """
         # Keep the same constructor as other discrete distributions.
         del warn
 

qmcpy/true_measure/abstract_true_measure.py

@@ -6,19 +6,25 @@
 from scipy import sparse
 
 
+def _clip_unit_interval(u):
+    """Clip unit-interval values away from endpoints for stable quantiles."""
+    eps = np.finfo(float).eps
+    return np.clip(u, eps, 1.0 - eps)
+
+
 class AbstractTrueMeasure(object):
 
     def __init__(self):
         prefix = "A concrete implementation of TrueMeasure must have "
         if not hasattr(self, "domain"):
             raise ParameterError(
                 prefix
-                + "self.domain, 2xd ndarray of domain lower bounds (first col) and upper bounds (second col)"
+                + "self.domain, (d, 2) ndarray of domain lower bounds (first col) and upper bounds (second col)"
             )
         if not hasattr(self, "range"):
             raise ParameterError(
                 prefix
-                + "self.range, 2xd ndarray of range lower bounds (first col) and upper bounds (second col)"
+                + "self.range, (d, 2) ndarray of range lower bounds (first col) and upper bounds (second col)"
             )
         if not hasattr(self, "parameters"):
             self.parameters = []
@@ -30,6 +36,52 @@ def _read_only_array(value):
         array.setflags(write=False)
         return array
 
+    @staticmethod
+    def _range_in_domain(transform_range, domain):
+        """Return whether a transform range is contained within a domain."""
+        try:
+            transform_range = np.asarray(transform_range)
+            domain = np.asarray(domain)
+        except (TypeError, ValueError):
+            return False
+
+        if (
+            transform_range.ndim != 2
+            or domain.ndim != 2
+            or transform_range.shape[1] != 2
+            or domain.shape[1] != 2
+            or transform_range.shape[0] == 0
+            or domain.shape[0] == 0
+        ):
+            return False
+
+        if not (
+            np.issubdtype(transform_range.dtype, np.number)
+            and np.issubdtype(domain.dtype, np.number)
+            and np.isrealobj(transform_range)
+            and np.isrealobj(domain)
+            and transform_range.dtype != np.bool_
+            and domain.dtype != np.bool_
+        ):
+            return False
+
+        if np.isnan(transform_range).any() or np.isnan(domain).any():
+            return False
+
+        if np.any(transform_range[:, 0] > transform_range[:, 1]) or np.any(
+            domain[:, 0] > domain[:, 1]
+        ):
+            return False
+
+        try:
+            transform_range, domain = np.broadcast_arrays(transform_range, domain)
+        except ValueError:
+            return False
+
+        lower_bounds_valid = np.all(domain[:, 0] <= transform_range[:, 0])
+        upper_bounds_valid = np.all(transform_range[:, 1] <= domain[:, 1])
+        return bool(lower_bounds_valid and upper_bounds_valid)
+
     def _set_moments(self, mean, variance, standard_deviation, covariance):
         self._mean = self._read_only_array(mean)
         self._variance = self._read_only_array(variance)
@@ -100,11 +152,11 @@ def _parse_sampler(self, sampler):
                 sampler.d
             )  # take the dimension from the sub-sampler (composed transform)
             self.discrete_distrib = self.transform.discrete_distrib
-            if (self.domain != self.transform.range).any():
+            if not self._range_in_domain(self.transform.range, self.domain):
                 self.sub_compatibility_error = True
             if self.transform.sub_compatibility_error:
                 raise ParameterError(
-                    "The sub-transform domain must match the sub-sub-transform range."
+                    "The nested sub-transform range must be contained within its transform domain."
                 )
         else:
             raise ParameterError(
@@ -147,7 +199,7 @@ def _jacobian_transform_r(self, x, return_weights):
         jac = None
         if self.sub_compatibility_error:
             raise ParameterError(
-                "The transform domain must match the sub-transform range."
+                "The sub-transform range must be contained within the transform domain."
             )
         if self.transform == self:  # is \Psi_0
             if return_weights:

qmcpy/true_measure/bernoulli_cont.py

@@ -75,6 +75,7 @@ def _transform(self, x):
         return tf
 
     def _weight(self, x):
+        in_support = np.all((0 <= x) & (x <= 1), axis=-1)
         w = np.zeros(x.shape, dtype=float)
         for j in range(self.d):
             C = (
@@ -83,7 +84,7 @@ def _weight(self, x):
                 else 2 * np.arctanh(1 - 2 * self.l[j]) / (1 - 2 * self.l[j])
             )
             w[..., j] = C * self.l[j] ** x[..., j] * (1 - self.l[j]) ** (1 - x[..., j])
-        return np.prod(w, -1)
+        return np.where(in_support, np.prod(w, -1), 0.0)
 
     def _spawn(self, sampler, dimension):
         if dimension == self.d:  # don't do anything if the dimension doesn't change

qmcpy/true_measure/brownian_motion.py

@@ -1,4 +1,5 @@
 from .gaussian import Gaussian
+from .abstract_true_measure import _clip_unit_interval
 from ..discrete_distribution import DigitalNetB2
 from ..util import ParameterError, ParameterWarning
 import warnings
@@ -252,7 +253,7 @@ def _spawn(self, sampler, dimension):
 
     def _transform(self, x):
         if self.decomp_type == "BROWNIANBRIDGE":
-            z = norm.ppf(x)
+            z = norm.ppf(_clip_unit_interval(x))
             w = self._bridge_transform(z)
             paths = self.drift_time_vec_plus_init + np.sqrt(self.diffusion) * w
             return paths[..., self._output_order]

qmcpy/true_measure/copula.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from .abstract_true_measure import AbstractTrueMeasure
+from .abstract_true_measure import AbstractTrueMeasure, _clip_unit_interval
 from ..util import DimensionError, MethodImplementationError, ParameterError
 
 
@@ -97,11 +97,6 @@ def _unit_weight_with_warning(self, x):
         return np.ones(x.shape[:-1], dtype=float)
 
 
-def _clip_unit_interval(u):
-    eps = np.finfo(float).eps
-    return np.clip(u, eps, 1.0 - eps)
-
-
 def _validate_marginals(marginals):
     try:
         parsed = list(marginals)

qmcpy/true_measure/gaussian.py

@@ -1,4 +1,4 @@
-from .abstract_true_measure import AbstractTrueMeasure
+from .abstract_true_measure import AbstractTrueMeasure, _clip_unit_interval
 from ..util import DimensionError, ParameterError
 from ..discrete_distribution import DigitalNetB2
 import numpy as np
@@ -158,6 +158,7 @@ def mvn_scipy(self, value):
         self._mvn_scipy_cache = value
 
     def _transform(self, x):
+        x = _clip_unit_interval(x)
         return self.mu + np.einsum("...ij,kj->...ik", norm.ppf(x), self.a)
 
     def _weight(self, t):

qmcpy/true_measure/johnsons_su.py

@@ -1,4 +1,4 @@
-from .abstract_true_measure import AbstractTrueMeasure
+from .abstract_true_measure import AbstractTrueMeasure, _clip_unit_interval
 from ..util import DimensionError, ParameterError
 from ..discrete_distribution import DigitalNetB2
 import numpy as np
@@ -92,6 +92,7 @@ def __init__(self, sampler, gamma=1, xi=1, delta=2, lam=2):
         )
 
     def _transform(self, x):
+        x = _clip_unit_interval(x)
         return self._lam * np.sinh((norm.ppf(x) - self._gamma) / self._delta) + self._xi
 
     def _weight(self, x):

qmcpy/true_measure/kumaraswamy.py

@@ -158,13 +158,17 @@ def _transform(self, x):
         return (1 - (1 - x) ** (1 / self.beta)) ** (1 / self.alpha)
 
     def _weight(self, x):
-        return np.prod(
-            self.alpha
-            * self.beta
-            * x ** (self.alpha - 1)
-            * (1 - x**self.alpha) ** (self.beta - 1),
-            -1,
-        )
+        in_support = np.all((0 <= x) & (x <= 1), axis=-1)
+        x_in_support = np.clip(x, 0, 1)
+        with np.errstate(divide="ignore", invalid="ignore"):
+            weight = np.prod(
+                self.alpha
+                * self.beta
+                * x_in_support ** (self.alpha - 1)
+                * (1 - x_in_support**self.alpha) ** (self.beta - 1),
+                -1,
+            )
+        return np.where(in_support, weight, 0.0)
 
     def _spawn(self, sampler, dimension):
         if dimension == self.d:  # don't do anything if the dimension doesn't change

qmcpy/true_measure/product_measure.py

@@ -105,25 +105,22 @@ def __init__(self, sampler, marginals):
         """
         Initialize a product measure from one sampler and several marginals.
 
-        Parameters
-        ----------
-        sampler : AbstractDiscreteDistribution
-            The sampler for the whole product measure. Its dimension must
-            equal the sum of the marginal dimensions.
+        Args:
+            sampler (AbstractDiscreteDistribution): Sampler for the whole product measure. Its dimension must equal the sum of the marginal dimensions.
+            marginals (Union[list, tuple]): Nonempty sequence of independent `AbstractTrueMeasure` instances to place side by side. A marginal may itself be multidimensional.
 
-        marginals : list or tuple of AbstractTrueMeasure
-            Independent true measures to place side by side. A marginal may
-            itself be multidimensional.
+        Raises:
+            ParameterError: If `sampler` is not an `AbstractDiscreteDistribution`, or if `marginals` is empty or contains a non-`AbstractTrueMeasure` value.
+            DimensionError: If a marginal is not dimension-preserving or the sampler dimension differs from the sum of the marginal dimensions.
 
-        Why one sampler?
-        ----------------
-        The product measure should be driven by one total-dimensional QMC
-        point set. We do not generate separate QMC samples from each marginal.
-        Instead, one sample u in [0,1]^d is split into blocks:
+        Note:
+            The product measure is driven by one total-dimensional QMC point
+            set. It does not generate separate QMC samples from each marginal.
+            Instead, one sample u in [0,1]^d is split into blocks:
 
-            u = (u_marginal_1, u_marginal_2, ..., u_marginal_k).
+                u = (u_marginal_1, u_marginal_2, ..., u_marginal_k).
 
-        This preserves the intended total-dimensional QMC construction.
+            This preserves the intended total-dimensional QMC construction.
         """
         if not isinstance(marginals, (list, tuple)) or len(marginals) == 0:
             raise ParameterError("ProductMeasure requires a nonempty list of marginals.")

qmcpy/true_measure/scipy_wrapper.py

@@ -1,4 +1,4 @@
-from .abstract_true_measure import AbstractTrueMeasure
+from .abstract_true_measure import AbstractTrueMeasure, _clip_unit_interval
 from ..util import DimensionError, ParameterError
 from ..discrete_distribution.abstract_discrete_distribution import (
     AbstractDiscreteDistribution,
@@ -107,8 +107,7 @@ def transform(self, u):
             )
 
         # Clip so we never hit exactly 0 or 1 inside norm.ppf.
-        eps = np.finfo(float).eps
-        u_clip = np.clip(u, eps, 1.0 - eps)
+        u_clip = _clip_unit_interval(u)
 
         # Map to i.i.d. standard normals.
         z = scipy.stats.norm.ppf(u_clip)
@@ -455,6 +454,7 @@ def _transform(self, x):
         if self._is_joint:
             return self._joint.transform(x)
 
+        x = _clip_unit_interval(x)
         t = np.empty_like(x, dtype=float)
         for j in range(self.d):
             t[..., j] = self.sds[j].ppf(x[..., j])

qmcpy/true_measure/student_t.py

@@ -2,6 +2,7 @@
 import scipy.stats as stats
 
 from ..util import ParameterError, DimensionError
+from .abstract_true_measure import _clip_unit_interval
 from .scipy_wrapper import SciPyWrapper
 
 
@@ -39,8 +40,7 @@ def __init__(self, loc, shape, df):
 
     @staticmethod
     def _clip_u(u):
-        eps = np.finfo(float).eps
-        return np.clip(u, eps, 1.0 - eps)
+        return _clip_unit_interval(u)
 
     def transform(self, u):
         u = np.asarray(u, dtype=float)

qmcpy/true_measure/uniform.py

@@ -100,7 +100,8 @@ def _transform(self, x):
         return x * self.delta + self.a
 
     def _weight(self, x):
-        return np.tile(self.inv_delta_prod, x.shape[:-1])
+        in_support = np.all((self.a <= x) & (x <= self.b), axis=-1)
+        return np.where(in_support, self.inv_delta_prod, 0.0)
 
     def _spawn(self, sampler, dimension):
         if dimension == self.d:  # don't do anything if the dimension doesn't change

qmcpy/true_measure/zero_inflated_exp_uniform.py

@@ -187,6 +187,17 @@ class ZeroInflatedExpUniform(SciPyWrapper):
     """
 
     def __init__(self, sampler, p_zero=0.4, lam=1.5, y_split=None):
+        r"""
+        Args:
+            sampler (Union[AbstractDiscreteDistribution, AbstractTrueMeasure]): One-dimensional sampler for the current construction. The deprecated `y_split` construction also accepts a two-dimensional sampler.
+            p_zero (float): Probability mass at zero, strictly between `0` and `1`. Defaults to `0.4`.
+            lam (float): Rate of the exponential component. Must be positive. Defaults to `1.5`.
+            y_split (Union[None, float]): Deprecated split point for the legacy two-dimensional construction. With a two-dimensional sampler, it must lie strictly between `0` and `1`. With a one-dimensional sampler, it is accepted for backward compatibility, emits a `DeprecationWarning`, and is otherwise ignored. Defaults to `None`.
+
+        Raises:
+            DimensionError: If the sampler dimension is incompatible with the selected construction.
+            ParameterError: If `p_zero`, `lam`, or a two-dimensional `y_split` is outside its valid range.
+        """
         if y_split is not None:
             warnings.warn(
                 "`y_split` is deprecated. The 2D zero-inflated "

scripts/colab_notebooks_manifest.json

@@ -44,6 +44,7 @@
   "demos/talk_paper_demos/Sorokin_random_LD_seq_QMC_fast_kernel_methods_2026/Sorokin_random_LD_seq_QMC_fast_kernel_methods_2026.ipynb",
   "demos/talk_paper_demos/pydata_chi_2023.ipynb",
   "demos/talk_paper_demos/why_add_q_to_mc_blog/why_add_q_to_mc_blog.ipynb",
+  "demos/true_measure_domain_inclusion.ipynb",
   "demos/vectorized_qmc.ipynb",
   "demos/vectorized_qmc_bayes.ipynb"
  ],

test/booktests/tb_true_measure_domain_inclusion.py

@@ -0,0 +1,18 @@
+import unittest
+from testbook import testbook
+from __init__ import TB_TIMEOUT, BaseNotebookTest
+
+
+class NotebookTests(BaseNotebookTest):
+
+    @testbook(
+        "../../demos/true_measure_domain_inclusion.ipynb",
+        execute=True,
+        timeout=TB_TIMEOUT,
+    )
+    def test_true_measure_domain_inclusion_notebook(self, tb):
+        pass
+
+
+if __name__ == "__main__":
+    unittest.main()