From 02b7aa4c808657c0a3bd22d7cd3641f9f6fb19b4 Mon Sep 17 00:00:00 2001 From: Laasya-73 <77721581+Laasya-73@users.noreply.github.com> Date: Wed, 19 Aug 2026 11:53:19 -0500 Subject: [PATCH 1/6] Add summary statistics to ProductMeasure --- qmcpy/true_measure/product_measure.py | 83 +++++++++++++++++++++++ test/test_product_measure.py | 96 +++++++++++++++++++++++++++ 2 files changed, 179 insertions(+) diff --git a/qmcpy/true_measure/product_measure.py b/qmcpy/true_measure/product_measure.py index b644666a3..dc2317f32 100644 --- a/qmcpy/true_measure/product_measure.py +++ b/qmcpy/true_measure/product_measure.py @@ -1,4 +1,5 @@ import numpy as np +from scipy import sparse from .abstract_true_measure import AbstractTrueMeasure from ..discrete_distribution.abstract_discrete_distribution import ( @@ -46,6 +47,9 @@ class ProductMeasure(AbstractTrueMeasure): Notes ----- + For independent marginal blocks, means, variances, and standard deviations + are concatenated in marginal order, while covariance is block diagonal. + Exact product weights are supported for direct marginal true measures. For recursively composed marginal measures, sampling is supported through QMCPy's recursive transform helper, but exact final-space product weights @@ -177,6 +181,85 @@ def __init__(self, sampler, marginals): super(ProductMeasure, self).__init__() + for statistic in ( + "mean", + "variance", + "standard_deviation", + "covariance", + ): + if all(hasattr(marginal, statistic) for marginal in self.marginals): + self.parameters.append(statistic) + + def _marginal_statistic(self, marginal, marginal_index, statistic): + """Return a statistic or identify the marginal that does not provide it.""" + try: + return getattr(marginal, statistic) + except AttributeError as error: + raise AttributeError( + f"ProductMeasure marginal {marginal_index} " + f"({type(marginal).__name__}) does not provide {statistic}." + ) from error + + def _concatenate_marginal_statistic(self, statistic): + """Concatenate a coordinate-wise statistic in marginal order.""" + values = [] + for marginal_index, marginal in enumerate(self.marginals): + value = self._marginal_statistic( + marginal, marginal_index, statistic + ) + value = np.atleast_1d(np.asarray(value)) + if value.shape != (marginal.d,): + raise DimensionError( + f"ProductMeasure marginal {marginal_index} " + f"({type(marginal).__name__}) {statistic} must have shape " + f"({marginal.d},), got {value.shape}." + ) + values.append(value) + + combined = self._read_only_array(np.concatenate(values)) + return self._scalar_if_univariate(combined) + + @property + def mean(self): + return self._concatenate_marginal_statistic("mean") + + @property + def variance(self): + return self._concatenate_marginal_statistic("variance") + + @property + def standard_deviation(self): + return self._concatenate_marginal_statistic("standard_deviation") + + @property + def covariance(self): + blocks = [] + for marginal_index, marginal in enumerate(self.marginals): + block = self._marginal_statistic( + marginal, marginal_index, "covariance" + ) + if sparse.issparse(block): + block = block.toarray() + block = np.atleast_2d(np.asarray(block)) + expected_shape = (marginal.d, marginal.d) + if block.shape != expected_shape: + raise DimensionError( + f"ProductMeasure marginal {marginal_index} " + f"({type(marginal).__name__}) covariance must have shape " + f"{expected_shape}, got {block.shape}." + ) + blocks.append(block) + + covariance = np.zeros( + (self.d, self.d), dtype=np.result_type(*[block.dtype for block in blocks]) + ) + start = 0 + for block in blocks: + stop = start + block.shape[0] + covariance[start:stop, start:stop] = block + start = stop + return self._read_only_array(covariance) + @staticmethod def _expand_bounds(bounds, dimension, name): """ diff --git a/test/test_product_measure.py b/test/test_product_measure.py index c91b313bd..9e217fd6a 100644 --- a/test/test_product_measure.py +++ b/test/test_product_measure.py @@ -48,6 +48,102 @@ def test_product_measure_replication_shape(): assert x.shape == (r, n, 2) +def test_product_measure_statistics_for_multiple_1d_marginals(): + marginals = [ + Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), + Uniform(DummySampler(1), lower_bound=-1.0, upper_bound=5.0), + ] + tm = ProductMeasure(sampler=DigitalNetB2(2, seed=23), marginals=marginals) + + np.testing.assert_allclose(tm.mean, [10.0, 2.0]) + np.testing.assert_allclose(tm.variance, [4.0 / 3.0, 3.0]) + np.testing.assert_allclose( + tm.standard_deviation, [np.sqrt(4.0 / 3.0), np.sqrt(3.0)] + ) + np.testing.assert_allclose(tm.covariance, np.diag([4.0 / 3.0, 3.0])) + + assert tm.mean.shape == (2,) + assert tm.variance.shape == (2,) + assert tm.standard_deviation.shape == (2,) + assert tm.covariance.shape == (2, 2) + for statistic in ("mean", "variance", "standard_deviation", "covariance"): + assert not getattr(tm, statistic).flags.writeable + + +def test_product_measure_normalizes_scalar_1d_statistics(): + marginals = [ + Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), + Gaussian(DummySampler(1), mean=2.0, covariance=9.0), + ] + for marginal in marginals: + assert isinstance(marginal.mean, float) + assert isinstance(marginal.variance, float) + assert isinstance(marginal.standard_deviation, float) + + tm = ProductMeasure(sampler=DigitalNetB2(2, seed=29), marginals=marginals) + + np.testing.assert_allclose(tm.mean, [10.0, 2.0]) + np.testing.assert_allclose(tm.variance, [4.0 / 3.0, 9.0]) + np.testing.assert_allclose( + tm.standard_deviation, [np.sqrt(4.0 / 3.0), 3.0] + ) + np.testing.assert_allclose(tm.covariance, np.diag([4.0 / 3.0, 9.0])) + assert tm.mean.shape == (2,) + assert tm.variance.shape == (2,) + assert tm.standard_deviation.shape == (2,) + assert tm.covariance.shape == (2, 2) + + +def test_product_measure_statistics_preserve_order_and_covariance_blocks(): + marginals = [ + Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), + Gaussian( + DummySampler(2), + mean=[2.0, 5.0], + covariance=[[2.0, 0.5], [0.5, 3.0]], + ), + ] + tm = ProductMeasure(sampler=DigitalNetB2(3, seed=31), marginals=marginals) + expected_covariance = np.array( + [ + [4.0 / 3.0, 0.0, 0.0], + [0.0, 2.0, 0.5], + [0.0, 0.5, 3.0], + ] + ) + + np.testing.assert_allclose(tm.mean, [10.0, 2.0, 5.0]) + np.testing.assert_allclose(tm.variance, [4.0 / 3.0, 2.0, 3.0]) + np.testing.assert_allclose( + tm.standard_deviation, + [np.sqrt(4.0 / 3.0), np.sqrt(2.0), np.sqrt(3.0)], + ) + np.testing.assert_allclose(tm.covariance, expected_covariance) + + assert tm.mean.shape == (3,) + assert tm.variance.shape == (3,) + assert tm.standard_deviation.shape == (3,) + assert tm.covariance.shape == (3, 3) + assert np.array_equal(tm.covariance[:1, 1:], np.zeros((1, 2))) + assert np.array_equal(tm.covariance[1:, :1], np.zeros((2, 1))) + + +def test_product_measure_missing_marginal_statistic_is_identified(): + marginals = [ + ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5), + Uniform(DummySampler(1), lower_bound=2.0, upper_bound=5.0), + ] + tm = ProductMeasure(sampler=DigitalNetB2(2, seed=23), marginals=marginals) + + np.testing.assert_allclose(tm.mean, [0.4, 3.5]) + assert "covariance" not in tm.parameters + with pytest.raises( + AttributeError, + match=r"marginal 0 \(ZeroInflatedExpUniform\) does not provide covariance", + ): + _ = tm.covariance + + def test_product_measure_marginals_with_different_dimensions(): n = 32 marginals = [ From 1514d5d6a47e12641b8b6383759b2360abc74f56 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 21 Aug 2026 21:16:49 +0000 Subject: [PATCH 2/6] Fix ProductMeasure.covariance to preserve sparse blocks Co-authored-by: fjhickernell <817530+fjhickernell@users.noreply.github.com> --- qmcpy/true_measure/product_measure.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/qmcpy/true_measure/product_measure.py b/qmcpy/true_measure/product_measure.py index dc2317f32..bdce45ac3 100644 --- a/qmcpy/true_measure/product_measure.py +++ b/qmcpy/true_measure/product_measure.py @@ -238,9 +238,8 @@ def covariance(self): block = self._marginal_statistic( marginal, marginal_index, "covariance" ) - if sparse.issparse(block): - block = block.toarray() - block = np.atleast_2d(np.asarray(block)) + if not sparse.issparse(block): + block = np.atleast_2d(np.asarray(block)) expected_shape = (marginal.d, marginal.d) if block.shape != expected_shape: raise DimensionError( @@ -250,11 +249,18 @@ def covariance(self): ) blocks.append(block) + if all(sparse.issparse(b) for b in blocks): + covariance = sparse.block_diag(blocks, format="dia") + covariance.data.setflags(write=False) + return covariance + + dense_blocks = [b.toarray() if sparse.issparse(b) else b for b in blocks] covariance = np.zeros( - (self.d, self.d), dtype=np.result_type(*[block.dtype for block in blocks]) + (self.d, self.d), + dtype=np.result_type(*[b.dtype for b in dense_blocks]), ) start = 0 - for block in blocks: + for block in dense_blocks: stop = start + block.shape[0] covariance[start:stop, start:stop] = block start = stop From 44f1fd6702228e688a3a749902a30f5f1aad7b7f Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 21 Aug 2026 21:39:06 +0000 Subject: [PATCH 3/6] Fix test assertions for sparse covariance matrix Co-authored-by: fjhickernell <817530+fjhickernell@users.noreply.github.com> --- test/test_product_measure.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/test/test_product_measure.py b/test/test_product_measure.py index 9e217fd6a..d4d712feb 100644 --- a/test/test_product_measure.py +++ b/test/test_product_measure.py @@ -1,5 +1,6 @@ import numpy as np import pytest +import scipy.sparse as sp import scipy.stats as stats from qmcpy import ( @@ -60,14 +61,18 @@ def test_product_measure_statistics_for_multiple_1d_marginals(): np.testing.assert_allclose( tm.standard_deviation, [np.sqrt(4.0 / 3.0), np.sqrt(3.0)] ) - np.testing.assert_allclose(tm.covariance, np.diag([4.0 / 3.0, 3.0])) + cov = tm.covariance + cov_dense = cov.toarray() if sp.issparse(cov) else cov + np.testing.assert_allclose(cov_dense, np.diag([4.0 / 3.0, 3.0])) assert tm.mean.shape == (2,) assert tm.variance.shape == (2,) assert tm.standard_deviation.shape == (2,) assert tm.covariance.shape == (2, 2) for statistic in ("mean", "variance", "standard_deviation", "covariance"): - assert not getattr(tm, statistic).flags.writeable + value = getattr(tm, statistic) + flags = value.data.flags if sp.issparse(value) else value.flags + assert not flags.writeable def test_product_measure_normalizes_scalar_1d_statistics(): From 536d3cdd15111e3280ecfddcf60f2bc200e72436 Mon Sep 17 00:00:00 2001 From: Laasya-73 <77721581+Laasya-73@users.noreply.github.com> Date: Mon, 24 Aug 2026 12:17:08 -0500 Subject: [PATCH 4/6] Preserve sparse ProductMeasure covariance --- qmcpy/true_measure/product_measure.py | 14 +++--- test/test_product_measure.py | 71 +++++++++++++++++++++++++-- 2 files changed, 74 insertions(+), 11 deletions(-) diff --git a/qmcpy/true_measure/product_measure.py b/qmcpy/true_measure/product_measure.py index bdce45ac3..c1eeda1e6 100644 --- a/qmcpy/true_measure/product_measure.py +++ b/qmcpy/true_measure/product_measure.py @@ -249,22 +249,24 @@ def covariance(self): ) blocks.append(block) - if all(sparse.issparse(b) for b in blocks): + if any(sparse.issparse(block) for block in blocks): covariance = sparse.block_diag(blocks, format="dia") - covariance.data.setflags(write=False) + data = covariance.data + data.setflags(write=False) + covariance.data = self._read_only_view(data) return covariance - dense_blocks = [b.toarray() if sparse.issparse(b) else b for b in blocks] covariance = np.zeros( (self.d, self.d), - dtype=np.result_type(*[b.dtype for b in dense_blocks]), + dtype=np.result_type(*[block.dtype for block in blocks]), ) start = 0 - for block in dense_blocks: + for block in blocks: stop = start + block.shape[0] covariance[start:stop, start:stop] = block start = stop - return self._read_only_array(covariance) + covariance.setflags(write=False) + return self._read_only_view(covariance) @staticmethod def _expand_bounds(bounds, dimension, name): diff --git a/test/test_product_measure.py b/test/test_product_measure.py index d4d712feb..67475a4bf 100644 --- a/test/test_product_measure.py +++ b/test/test_product_measure.py @@ -92,7 +92,9 @@ def test_product_measure_normalizes_scalar_1d_statistics(): np.testing.assert_allclose( tm.standard_deviation, [np.sqrt(4.0 / 3.0), 3.0] ) - np.testing.assert_allclose(tm.covariance, np.diag([4.0 / 3.0, 9.0])) + np.testing.assert_allclose( + tm.covariance.toarray(), np.diag([4.0 / 3.0, 9.0]) + ) assert tm.mean.shape == (2,) assert tm.variance.shape == (2,) assert tm.standard_deviation.shape == (2,) @@ -123,14 +125,73 @@ def test_product_measure_statistics_preserve_order_and_covariance_blocks(): tm.standard_deviation, [np.sqrt(4.0 / 3.0), np.sqrt(2.0), np.sqrt(3.0)], ) - np.testing.assert_allclose(tm.covariance, expected_covariance) + covariance = tm.covariance.tocsr() + np.testing.assert_allclose(covariance.toarray(), expected_covariance) assert tm.mean.shape == (3,) assert tm.variance.shape == (3,) assert tm.standard_deviation.shape == (3,) - assert tm.covariance.shape == (3, 3) - assert np.array_equal(tm.covariance[:1, 1:], np.zeros((1, 2))) - assert np.array_equal(tm.covariance[1:, :1], np.zeros((2, 1))) + assert covariance.shape == (3, 3) + assert covariance[:1, 1:].nnz == 0 + assert covariance[1:, :1].nnz == 0 + + +def test_product_measure_mixed_covariance_blocks_remain_sparse(): + d = 128 + tm = ProductMeasure( + DummySampler(d + 2), + [ + Uniform(DummySampler(d)), + Gaussian( + DummySampler(2), + covariance=np.array([[1.0, 0.5], [0.5, 1.0]]), + ), + ], + ) + + covariance = tm.covariance + expected = sp.block_diag( + [marginal.covariance for marginal in tm.marginals], format="dia" + ) + + assert sp.issparse(covariance) + assert covariance.format == "dia" + assert covariance.shape == (d + 2, d + 2) + difference = (covariance - expected).tocsr() + difference.eliminate_zeros() + assert difference.nnz == 0 + covariance_csr = covariance.tocsr() + np.testing.assert_allclose( + covariance_csr[-2:, -2:].toarray(), [[1.0, 0.5], [0.5, 1.0]] + ) + assert covariance_csr[:d, d:].nnz == 0 + assert covariance_csr[d:, :d].nnz == 0 + assert not covariance.data.flags.writeable + with pytest.raises(ValueError): + covariance.data.setflags(write=True) + + +def test_product_measure_dense_covariance_cannot_be_made_writeable(): + tm = ProductMeasure( + DummySampler(3), + [ + Gaussian(DummySampler(1), covariance=2.0), + Gaussian( + DummySampler(2), covariance=[[3.0, 0.25], [0.25, 4.0]] + ), + ], + ) + + covariance = tm.covariance + + assert isinstance(covariance, np.ndarray) + np.testing.assert_allclose( + covariance, + [[2.0, 0.0, 0.0], [0.0, 3.0, 0.25], [0.0, 0.25, 4.0]], + ) + assert not covariance.flags.writeable + with pytest.raises(ValueError): + covariance.setflags(write=True) def test_product_measure_missing_marginal_statistic_is_identified(): From 4c1dd27c8f71f7740d09c7295c89c930694dec28 Mon Sep 17 00:00:00 2001 From: Laasya-73 <77721581+Laasya-73@users.noreply.github.com> Date: Sun, 30 Aug 2026 10:44:41 -0500 Subject: [PATCH 5/6] Improve ProductMeasure statistics handling --- qmcpy/true_measure/product_measure.py | 66 ++++++++++++++++++++----- test/test_product_measure.py | 69 ++++++++++++++++++++++++++- 2 files changed, 122 insertions(+), 13 deletions(-) diff --git a/qmcpy/true_measure/product_measure.py b/qmcpy/true_measure/product_measure.py index c1eeda1e6..3f860e84b 100644 --- a/qmcpy/true_measure/product_measure.py +++ b/qmcpy/true_measure/product_measure.py @@ -5,7 +5,7 @@ from ..discrete_distribution.abstract_discrete_distribution import ( AbstractDiscreteDistribution, ) -from ..util import DimensionError, ParameterError +from ..util import DimensionError, ParameterError, _univ_repr class ProductMeasure(AbstractTrueMeasure): @@ -181,6 +181,11 @@ def __init__(self, sampler, marginals): super(ProductMeasure, self).__init__() + self._mean_cache = None + self._variance_cache = None + self._standard_deviation_cache = None + self._covariance_cache = None + for statistic in ( "mean", "variance", @@ -221,18 +226,26 @@ def _concatenate_marginal_statistic(self, statistic): @property def mean(self): - return self._concatenate_marginal_statistic("mean") + if self._mean_cache is None: + self._mean_cache = self._concatenate_marginal_statistic("mean") + return self._mean_cache @property def variance(self): - return self._concatenate_marginal_statistic("variance") + if self._variance_cache is None: + self._variance_cache = self._concatenate_marginal_statistic("variance") + return self._variance_cache @property def standard_deviation(self): - return self._concatenate_marginal_statistic("standard_deviation") + if self._standard_deviation_cache is None: + self._standard_deviation_cache = self._concatenate_marginal_statistic( + "standard_deviation" + ) + return self._standard_deviation_cache - @property - def covariance(self): + def _compute_covariance(self): + """Build and protect the block-diagonal marginal covariance.""" blocks = [] for marginal_index, marginal in enumerate(self.marginals): block = self._marginal_statistic( @@ -250,11 +263,14 @@ def covariance(self): blocks.append(block) if any(sparse.issparse(block) for block in blocks): - covariance = sparse.block_diag(blocks, format="dia") - data = covariance.data - data.setflags(write=False) - covariance.data = self._read_only_view(data) - return covariance + covariance = sparse.block_diag(blocks, format="csr") + # Rebuild from a read-only base so writes cannot be re-enabled. + data = self._read_only_array(covariance.data) + return sparse.csr_matrix( + (data, covariance.indices, covariance.indptr), + shape=covariance.shape, + copy=False, + ) covariance = np.zeros( (self.d, self.d), @@ -268,6 +284,34 @@ def covariance(self): covariance.setflags(write=False) return self._read_only_view(covariance) + @property + def covariance(self): + if self._covariance_cache is None: + self._covariance_cache = self._compute_covariance() + return self._covariance_cache + + def __repr__(self): + """Represent ProductMeasure without expanding marginal sparse matrices.""" + lines = [f"{type(self).__name__} (AbstractTrueMeasure)"] + for parameter in dict.fromkeys(self.parameters): + if parameter == "marginals": + marginals = ", ".join( + f"{type(marginal).__name__}(d={marginal.d})" + for marginal in self.marginals + ) + lines.append(f" {parameter:<15} [{marginals}]") + elif parameter == "covariance" and sparse.issparse(self.covariance): + covariance = self.covariance + summary = ( + f"sparse {covariance.format.upper()}, " + f"shape={covariance.shape}, nnz={covariance.nnz}" + ) + lines.append(f" {parameter:<15} {summary}") + else: + formatted = _univ_repr(self, "AbstractTrueMeasure", [parameter]) + lines.extend(formatted.splitlines()[1:]) + return "\n".join(lines) + @staticmethod def _expand_bounds(bounds, dimension, name): """ diff --git a/test/test_product_measure.py b/test/test_product_measure.py index 67475a4bf..d3cbd8b60 100644 --- a/test/test_product_measure.py +++ b/test/test_product_measure.py @@ -151,11 +151,11 @@ def test_product_measure_mixed_covariance_blocks_remain_sparse(): covariance = tm.covariance expected = sp.block_diag( - [marginal.covariance for marginal in tm.marginals], format="dia" + [marginal.covariance for marginal in tm.marginals], format="csr" ) assert sp.issparse(covariance) - assert covariance.format == "dia" + assert covariance.format == "csr" assert covariance.shape == (d + 2, d + 2) difference = (covariance - expected).tocsr() difference.eliminate_zeros() @@ -171,6 +171,71 @@ def test_product_measure_mixed_covariance_blocks_remain_sparse(): covariance.data.setflags(write=True) +def test_product_measure_sparse_covariance_repr_is_compact(): + d = 128 + tm = ProductMeasure( + DummySampler(d + 2), + [ + Uniform(DummySampler(d)), + Gaussian( + DummySampler(2), + covariance=np.array([[1.0, 0.5], [0.5, 1.0]]), + ), + ], + ) + + representation = repr(tm) + + assert "sparse CSR" in representation + assert "shape=(130, 130)" in representation + assert "nnz=132" in representation + assert "Coords" not in representation + assert "(127, 127)" not in representation + + +def test_product_measure_statistics_are_lazily_cached(): + tm = ProductMeasure( + DummySampler(3), + [ + Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), + Gaussian( + DummySampler(2), + mean=[2.0, 5.0], + covariance=[[2.0, 0.5], [0.5, 3.0]], + ), + ], + ) + + cache = vars(tm) + assert cache["_mean_cache"] is None + assert cache["_variance_cache"] is None + assert cache["_standard_deviation_cache"] is None + assert cache["_covariance_cache"] is None + + mean = tm.mean + np.testing.assert_allclose(mean, [10.0, 2.0, 5.0]) + assert tm.mean is mean + assert cache["_variance_cache"] is None + assert cache["_standard_deviation_cache"] is None + assert cache["_covariance_cache"] is None + + variance = tm.variance + standard_deviation = tm.standard_deviation + covariance = tm.covariance + np.testing.assert_allclose(variance, [4.0 / 3.0, 2.0, 3.0]) + np.testing.assert_allclose( + standard_deviation, + [np.sqrt(4.0 / 3.0), np.sqrt(2.0), np.sqrt(3.0)], + ) + np.testing.assert_allclose( + covariance.toarray(), + [[4.0 / 3.0, 0.0, 0.0], [0.0, 2.0, 0.5], [0.0, 0.5, 3.0]], + ) + assert tm.variance is variance + assert tm.standard_deviation is standard_deviation + assert tm.covariance is covariance + + def test_product_measure_dense_covariance_cannot_be_made_writeable(): tm = ProductMeasure( DummySampler(3), From b13a747fa0696bb4f986752e680142a411fa78fe Mon Sep 17 00:00:00 2001 From: Laasya-73 <77721581+Laasya-73@users.noreply.github.com> Date: Mon, 31 Aug 2026 20:17:49 -0500 Subject: [PATCH 6/6] Organize ProductMeasure tests with unittest --- test/test_product_measure.py | 941 ++++++++++++++++++----------------- 1 file changed, 472 insertions(+), 469 deletions(-) diff --git a/test/test_product_measure.py b/test/test_product_measure.py index d3cbd8b60..f49d5ded5 100644 --- a/test/test_product_measure.py +++ b/test/test_product_measure.py @@ -1,5 +1,6 @@ +import unittest + import numpy as np -import pytest import scipy.sparse as sp import scipy.stats as stats @@ -17,488 +18,490 @@ from qmcpy.util import DimensionError, ParameterError -def test_product_measure_zero_inflated_with_scipy_uniform_shape(): - n = 32 - marginals = [ - ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5), - SciPyWrapper(DummySampler(1), stats.uniform(loc=2.0, scale=3.0)), - ] - tm = ProductMeasure(sampler=DigitalNetB2(2, seed=23), marginals=marginals) - - x = tm(n) - - assert x.shape == (n, 2) - assert np.any(x[:, 0] == 0.0) - assert np.all((2.0 <= x[:, 1]) & (x[:, 1] <= 5.0)) - - -def test_product_measure_replication_shape(): - n = 16 - r = 3 - marginals = [ - ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5), - Uniform(DummySampler(1), lower_bound=2.0, upper_bound=5.0), - ] - tm = ProductMeasure( - sampler=DigitalNetB2(2, seed=23, replications=r), - marginals=marginals, - ) - - x = tm(n) - - assert x.shape == (r, n, 2) - - -def test_product_measure_statistics_for_multiple_1d_marginals(): - marginals = [ - Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), - Uniform(DummySampler(1), lower_bound=-1.0, upper_bound=5.0), - ] - tm = ProductMeasure(sampler=DigitalNetB2(2, seed=23), marginals=marginals) - - np.testing.assert_allclose(tm.mean, [10.0, 2.0]) - np.testing.assert_allclose(tm.variance, [4.0 / 3.0, 3.0]) - np.testing.assert_allclose( - tm.standard_deviation, [np.sqrt(4.0 / 3.0), np.sqrt(3.0)] - ) - cov = tm.covariance - cov_dense = cov.toarray() if sp.issparse(cov) else cov - np.testing.assert_allclose(cov_dense, np.diag([4.0 / 3.0, 3.0])) - - assert tm.mean.shape == (2,) - assert tm.variance.shape == (2,) - assert tm.standard_deviation.shape == (2,) - assert tm.covariance.shape == (2, 2) - for statistic in ("mean", "variance", "standard_deviation", "covariance"): - value = getattr(tm, statistic) - flags = value.data.flags if sp.issparse(value) else value.flags - assert not flags.writeable - - -def test_product_measure_normalizes_scalar_1d_statistics(): - marginals = [ - Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), + +class TestProductMeasure(unittest.TestCase): + def test_zero_inflated_with_scipy_uniform_shape(self): + n = 32 + marginals = [ + ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5), + SciPyWrapper(DummySampler(1), stats.uniform(loc=2.0, scale=3.0)), + ] + tm = ProductMeasure(sampler=DigitalNetB2(2, seed=23), marginals=marginals) + + x = tm(n) + + self.assertEqual(x.shape, (n, 2)) + self.assertTrue(np.any(x[:, 0] == 0.0)) + self.assertTrue(np.all((2.0 <= x[:, 1]) & (x[:, 1] <= 5.0))) + + + def test_replication_shape(self): + n = 16 + r = 3 + marginals = [ + ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5), + Uniform(DummySampler(1), lower_bound=2.0, upper_bound=5.0), + ] + tm = ProductMeasure( + sampler=DigitalNetB2(2, seed=23, replications=r), + marginals=marginals, + ) + + x = tm(n) + + self.assertEqual(x.shape, (r, n, 2)) + + + def test_statistics_for_multiple_1d_marginals(self): + marginals = [ + Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), + Uniform(DummySampler(1), lower_bound=-1.0, upper_bound=5.0), + ] + tm = ProductMeasure(sampler=DigitalNetB2(2, seed=23), marginals=marginals) + + np.testing.assert_allclose(tm.mean, [10.0, 2.0]) + np.testing.assert_allclose(tm.variance, [4.0 / 3.0, 3.0]) + np.testing.assert_allclose( + tm.standard_deviation, [np.sqrt(4.0 / 3.0), np.sqrt(3.0)] + ) + cov = tm.covariance + cov_dense = cov.toarray() if sp.issparse(cov) else cov + np.testing.assert_allclose(cov_dense, np.diag([4.0 / 3.0, 3.0])) + + self.assertEqual(tm.mean.shape, (2,)) + self.assertEqual(tm.variance.shape, (2,)) + self.assertEqual(tm.standard_deviation.shape, (2,)) + self.assertEqual(tm.covariance.shape, (2, 2)) + for statistic in ("mean", "variance", "standard_deviation", "covariance"): + value = getattr(tm, statistic) + flags = value.data.flags if sp.issparse(value) else value.flags + self.assertFalse(flags.writeable) + + + def test_normalizes_scalar_1d_statistics(self): + marginals = [ + Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), Gaussian(DummySampler(1), mean=2.0, covariance=9.0), - ] - for marginal in marginals: - assert isinstance(marginal.mean, float) - assert isinstance(marginal.variance, float) - assert isinstance(marginal.standard_deviation, float) - - tm = ProductMeasure(sampler=DigitalNetB2(2, seed=29), marginals=marginals) - - np.testing.assert_allclose(tm.mean, [10.0, 2.0]) - np.testing.assert_allclose(tm.variance, [4.0 / 3.0, 9.0]) - np.testing.assert_allclose( - tm.standard_deviation, [np.sqrt(4.0 / 3.0), 3.0] - ) - np.testing.assert_allclose( - tm.covariance.toarray(), np.diag([4.0 / 3.0, 9.0]) - ) - assert tm.mean.shape == (2,) - assert tm.variance.shape == (2,) - assert tm.standard_deviation.shape == (2,) - assert tm.covariance.shape == (2, 2) - - -def test_product_measure_statistics_preserve_order_and_covariance_blocks(): - marginals = [ - Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), - Gaussian( - DummySampler(2), - mean=[2.0, 5.0], - covariance=[[2.0, 0.5], [0.5, 3.0]], - ), - ] - tm = ProductMeasure(sampler=DigitalNetB2(3, seed=31), marginals=marginals) - expected_covariance = np.array( - [ - [4.0 / 3.0, 0.0, 0.0], - [0.0, 2.0, 0.5], - [0.0, 0.5, 3.0], ] - ) - - np.testing.assert_allclose(tm.mean, [10.0, 2.0, 5.0]) - np.testing.assert_allclose(tm.variance, [4.0 / 3.0, 2.0, 3.0]) - np.testing.assert_allclose( - tm.standard_deviation, - [np.sqrt(4.0 / 3.0), np.sqrt(2.0), np.sqrt(3.0)], - ) - covariance = tm.covariance.tocsr() - np.testing.assert_allclose(covariance.toarray(), expected_covariance) - - assert tm.mean.shape == (3,) - assert tm.variance.shape == (3,) - assert tm.standard_deviation.shape == (3,) - assert covariance.shape == (3, 3) - assert covariance[:1, 1:].nnz == 0 - assert covariance[1:, :1].nnz == 0 - - -def test_product_measure_mixed_covariance_blocks_remain_sparse(): - d = 128 - tm = ProductMeasure( - DummySampler(d + 2), - [ - Uniform(DummySampler(d)), + for marginal in marginals: + self.assertIsInstance(marginal.mean, float) + self.assertIsInstance(marginal.variance, float) + self.assertIsInstance(marginal.standard_deviation, float) + + tm = ProductMeasure(sampler=DigitalNetB2(2, seed=29), marginals=marginals) + + np.testing.assert_allclose(tm.mean, [10.0, 2.0]) + np.testing.assert_allclose(tm.variance, [4.0 / 3.0, 9.0]) + np.testing.assert_allclose( + tm.standard_deviation, [np.sqrt(4.0 / 3.0), 3.0] + ) + np.testing.assert_allclose( + tm.covariance.toarray(), np.diag([4.0 / 3.0, 9.0]) + ) + self.assertEqual(tm.mean.shape, (2,)) + self.assertEqual(tm.variance.shape, (2,)) + self.assertEqual(tm.standard_deviation.shape, (2,)) + self.assertEqual(tm.covariance.shape, (2, 2)) + + + def test_statistics_preserve_order_and_covariance_blocks(self): + marginals = [ + Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), Gaussian( DummySampler(2), - covariance=np.array([[1.0, 0.5], [0.5, 1.0]]), + mean=[2.0, 5.0], + covariance=[[2.0, 0.5], [0.5, 3.0]], ), - ], - ) - - covariance = tm.covariance - expected = sp.block_diag( - [marginal.covariance for marginal in tm.marginals], format="csr" - ) - - assert sp.issparse(covariance) - assert covariance.format == "csr" - assert covariance.shape == (d + 2, d + 2) - difference = (covariance - expected).tocsr() - difference.eliminate_zeros() - assert difference.nnz == 0 - covariance_csr = covariance.tocsr() - np.testing.assert_allclose( - covariance_csr[-2:, -2:].toarray(), [[1.0, 0.5], [0.5, 1.0]] - ) - assert covariance_csr[:d, d:].nnz == 0 - assert covariance_csr[d:, :d].nnz == 0 - assert not covariance.data.flags.writeable - with pytest.raises(ValueError): - covariance.data.setflags(write=True) - - -def test_product_measure_sparse_covariance_repr_is_compact(): - d = 128 - tm = ProductMeasure( - DummySampler(d + 2), - [ - Uniform(DummySampler(d)), + ] + tm = ProductMeasure(sampler=DigitalNetB2(3, seed=31), marginals=marginals) + expected_covariance = np.array( + [ + [4.0 / 3.0, 0.0, 0.0], + [0.0, 2.0, 0.5], + [0.0, 0.5, 3.0], + ] + ) + + np.testing.assert_allclose(tm.mean, [10.0, 2.0, 5.0]) + np.testing.assert_allclose(tm.variance, [4.0 / 3.0, 2.0, 3.0]) + np.testing.assert_allclose( + tm.standard_deviation, + [np.sqrt(4.0 / 3.0), np.sqrt(2.0), np.sqrt(3.0)], + ) + covariance = tm.covariance.tocsr() + np.testing.assert_allclose(covariance.toarray(), expected_covariance) + + self.assertEqual(tm.mean.shape, (3,)) + self.assertEqual(tm.variance.shape, (3,)) + self.assertEqual(tm.standard_deviation.shape, (3,)) + self.assertEqual(covariance.shape, (3, 3)) + self.assertEqual(covariance[:1, 1:].nnz, 0) + self.assertEqual(covariance[1:, :1].nnz, 0) + + + def test_mixed_covariance_blocks_remain_sparse(self): + d = 128 + tm = ProductMeasure( + DummySampler(d + 2), + [ + Uniform(DummySampler(d)), + Gaussian( + DummySampler(2), + covariance=np.array([[1.0, 0.5], [0.5, 1.0]]), + ), + ], + ) + + covariance = tm.covariance + expected = sp.block_diag( + [marginal.covariance for marginal in tm.marginals], format="csr" + ) + + self.assertTrue(sp.issparse(covariance)) + self.assertEqual(covariance.format, "csr") + self.assertEqual(covariance.shape, (d + 2, d + 2)) + difference = (covariance - expected).tocsr() + difference.eliminate_zeros() + self.assertEqual(difference.nnz, 0) + covariance_csr = covariance.tocsr() + np.testing.assert_allclose( + covariance_csr[-2:, -2:].toarray(), [[1.0, 0.5], [0.5, 1.0]] + ) + self.assertEqual(covariance_csr[:d, d:].nnz, 0) + self.assertEqual(covariance_csr[d:, :d].nnz, 0) + self.assertFalse(covariance.data.flags.writeable) + with self.assertRaises(ValueError): + covariance.data.setflags(write=True) + + + def test_sparse_covariance_repr_is_compact(self): + d = 128 + tm = ProductMeasure( + DummySampler(d + 2), + [ + Uniform(DummySampler(d)), + Gaussian( + DummySampler(2), + covariance=np.array([[1.0, 0.5], [0.5, 1.0]]), + ), + ], + ) + + representation = repr(tm) + + self.assertIn("sparse CSR", representation) + self.assertIn("shape=(130, 130)", representation) + self.assertIn("nnz=132", representation) + self.assertNotIn("Coords", representation) + self.assertNotIn("(127, 127)", representation) + + + def test_statistics_are_lazily_cached(self): + tm = ProductMeasure( + DummySampler(3), + [ + Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), + Gaussian( + DummySampler(2), + mean=[2.0, 5.0], + covariance=[[2.0, 0.5], [0.5, 3.0]], + ), + ], + ) + + cache = vars(tm) + self.assertIsNone(cache["_mean_cache"]) + self.assertIsNone(cache["_variance_cache"]) + self.assertIsNone(cache["_standard_deviation_cache"]) + self.assertIsNone(cache["_covariance_cache"]) + + mean = tm.mean + np.testing.assert_allclose(mean, [10.0, 2.0, 5.0]) + self.assertIs(tm.mean, mean) + self.assertIsNone(cache["_variance_cache"]) + self.assertIsNone(cache["_standard_deviation_cache"]) + self.assertIsNone(cache["_covariance_cache"]) + + variance = tm.variance + standard_deviation = tm.standard_deviation + covariance = tm.covariance + np.testing.assert_allclose(variance, [4.0 / 3.0, 2.0, 3.0]) + np.testing.assert_allclose( + standard_deviation, + [np.sqrt(4.0 / 3.0), np.sqrt(2.0), np.sqrt(3.0)], + ) + np.testing.assert_allclose( + covariance.toarray(), + [[4.0 / 3.0, 0.0, 0.0], [0.0, 2.0, 0.5], [0.0, 0.5, 3.0]], + ) + self.assertIs(tm.variance, variance) + self.assertIs(tm.standard_deviation, standard_deviation) + self.assertIs(tm.covariance, covariance) + + + def test_dense_covariance_cannot_be_made_writeable(self): + tm = ProductMeasure( + DummySampler(3), + [ + Gaussian(DummySampler(1), covariance=2.0), + Gaussian( + DummySampler(2), covariance=[[3.0, 0.25], [0.25, 4.0]] + ), + ], + ) + + covariance = tm.covariance + + self.assertIsInstance(covariance, np.ndarray) + np.testing.assert_allclose( + covariance, + [[2.0, 0.0, 0.0], [0.0, 3.0, 0.25], [0.0, 0.25, 4.0]], + ) + self.assertFalse(covariance.flags.writeable) + with self.assertRaises(ValueError): + covariance.setflags(write=True) + + + def test_missing_marginal_statistic_is_identified(self): + marginals = [ + ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5), + Uniform(DummySampler(1), lower_bound=2.0, upper_bound=5.0), + ] + tm = ProductMeasure(sampler=DigitalNetB2(2, seed=23), marginals=marginals) + + np.testing.assert_allclose(tm.mean, [0.4, 3.5]) + self.assertNotIn("covariance", tm.parameters) + with self.assertRaisesRegex( + AttributeError, + r"marginal 0 \(ZeroInflatedExpUniform\) does not provide covariance", + ): + _ = tm.covariance + + + def test_marginals_with_different_dimensions(self): + n = 32 + marginals = [ Gaussian( DummySampler(2), - covariance=np.array([[1.0, 0.5], [0.5, 1.0]]), + mean=[1.0, -1.0], + covariance=[[2.0, 0.25], [0.25, 1.0]], ), - ], - ) + ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5), + ] + tm = ProductMeasure(sampler=DigitalNetB2(3, seed=31), marginals=marginals) - representation = repr(tm) + x = tm(n) - assert "sparse CSR" in representation - assert "shape=(130, 130)" in representation - assert "nnz=132" in representation - assert "Coords" not in representation - assert "(127, 127)" not in representation + self.assertEqual(tm.d, 3) + np.testing.assert_array_equal(tm.marginal_dimensions, np.array([2, 1])) + self.assertEqual(x.shape, (n, 3)) + self.assertTrue(np.all(np.isfinite(x[:, :2]))) + self.assertTrue(np.all(x[:, 2] >= 0.0)) -def test_product_measure_statistics_are_lazily_cached(): - tm = ProductMeasure( - DummySampler(3), - [ - Uniform(DummySampler(1), lower_bound=8.0, upper_bound=12.0), - Gaussian( + def test_block_split_range_and_weight_product(self): + n = 16 + marginals = [ + Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0), + Uniform( DummySampler(2), - mean=[2.0, 5.0], - covariance=[[2.0, 0.5], [0.5, 3.0]], + lower_bound=[20.0, 30.0], + upper_bound=[24.0, 36.0], ), - ], - ) - - cache = vars(tm) - assert cache["_mean_cache"] is None - assert cache["_variance_cache"] is None - assert cache["_standard_deviation_cache"] is None - assert cache["_covariance_cache"] is None - - mean = tm.mean - np.testing.assert_allclose(mean, [10.0, 2.0, 5.0]) - assert tm.mean is mean - assert cache["_variance_cache"] is None - assert cache["_standard_deviation_cache"] is None - assert cache["_covariance_cache"] is None - - variance = tm.variance - standard_deviation = tm.standard_deviation - covariance = tm.covariance - np.testing.assert_allclose(variance, [4.0 / 3.0, 2.0, 3.0]) - np.testing.assert_allclose( - standard_deviation, - [np.sqrt(4.0 / 3.0), np.sqrt(2.0), np.sqrt(3.0)], - ) - np.testing.assert_allclose( - covariance.toarray(), - [[4.0 / 3.0, 0.0, 0.0], [0.0, 2.0, 0.5], [0.0, 0.5, 3.0]], - ) - assert tm.variance is variance - assert tm.standard_deviation is standard_deviation - assert tm.covariance is covariance - - -def test_product_measure_dense_covariance_cannot_be_made_writeable(): - tm = ProductMeasure( - DummySampler(3), - [ - Gaussian(DummySampler(1), covariance=2.0), - Gaussian( - DummySampler(2), covariance=[[3.0, 0.25], [0.25, 4.0]] + ] + tm = ProductMeasure(sampler=DigitalNetB2(3, seed=41), marginals=marginals) + + u = tm.discrete_distrib.gen_samples(n) + x = tm._transform(u) + x_call, jac = tm(n, return_weights=True) + expected = np.concatenate( + [ + marginals[0]._jacobian_transform_r(u[..., :1], return_weights=False), + marginals[1]._jacobian_transform_r(u[..., 1:], return_weights=False), + ], + axis=-1, + ) + expected_range = np.array([[10.0, 12.0], [20.0, 24.0], [30.0, 36.0]]) + + self.assertEqual(x.shape, (n, 3)) + np.testing.assert_allclose(tm.range, expected_range) + np.testing.assert_allclose(x, expected) + self.assertTrue(np.all((10.0 <= x[:, 0]) & (x[:, 0] <= 12.0))) + self.assertTrue(np.all((20.0 <= x[:, 1]) & (x[:, 1] <= 24.0))) + self.assertTrue(np.all((30.0 <= x[:, 2]) & (x[:, 2] <= 36.0))) + np.testing.assert_allclose(tm._weight(x), 1.0 / (2.0 * 4.0 * 6.0)) + self.assertEqual(x_call.shape, (n, 3)) + np.testing.assert_allclose(jac, 2.0 * 4.0 * 6.0) + + + def test_invalid_inputs(self): + with self.assertRaisesRegex(ParameterError, "nonempty list of marginals"): + ProductMeasure(sampler=DigitalNetB2(1, seed=7), marginals=[]) + + with self.assertRaisesRegex(ParameterError, "marginal"): + ProductMeasure(sampler=DigitalNetB2(1, seed=7), marginals=[object()]) + + with self.assertRaisesRegex(ParameterError, "AbstractDiscreteDistribution"): + ProductMeasure(sampler=object(), marginals=[Uniform(DummySampler(1))]) + + marginals = [Uniform(DummySampler(1))] + with self.assertRaisesRegex(DimensionError, "sum of marginal dimensions"): + ProductMeasure(sampler=DigitalNetB2(2, seed=7), marginals=marginals) + + + def test_rejects_non_dimension_preserving_marginal(self): + marginal = AcceptanceRejection( + DigitalNetB2(2, seed=7), + lambda x: np.ones(len(x)), + 1.0, + 1.0, + ) + + with self.assertRaisesRegex(DimensionError, "dimension-preserving"): + ProductMeasure(DigitalNetB2(2, seed=11), [marginal]) + + + def test_spawn_preserves_marginal_blocks_and_replaces_outer_sampler(self): + marginals = [ + Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0), + Uniform( + DummySampler(2), + lower_bound=[20.0, 30.0], + upper_bound=[24.0, 36.0], ), - ], - ) - - covariance = tm.covariance - - assert isinstance(covariance, np.ndarray) - np.testing.assert_allclose( - covariance, - [[2.0, 0.0, 0.0], [0.0, 3.0, 0.25], [0.0, 0.25, 4.0]], - ) - assert not covariance.flags.writeable - with pytest.raises(ValueError): - covariance.setflags(write=True) - - -def test_product_measure_missing_marginal_statistic_is_identified(): - marginals = [ - ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5), - Uniform(DummySampler(1), lower_bound=2.0, upper_bound=5.0), - ] - tm = ProductMeasure(sampler=DigitalNetB2(2, seed=23), marginals=marginals) - - np.testing.assert_allclose(tm.mean, [0.4, 3.5]) - assert "covariance" not in tm.parameters - with pytest.raises( - AttributeError, - match=r"marginal 0 \(ZeroInflatedExpUniform\) does not provide covariance", - ): - _ = tm.covariance - - -def test_product_measure_marginals_with_different_dimensions(): - n = 32 - marginals = [ - Gaussian( - DummySampler(2), - mean=[1.0, -1.0], - covariance=[[2.0, 0.25], [0.25, 1.0]], - ), - ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5), - ] - tm = ProductMeasure(sampler=DigitalNetB2(3, seed=31), marginals=marginals) - - x = tm(n) - - assert tm.d == 3 - assert np.array_equal(tm.marginal_dimensions, np.array([2, 1])) - assert x.shape == (n, 3) - assert np.all(np.isfinite(x[:, :2])) - assert np.all(x[:, 2] >= 0.0) - - -def test_product_measure_block_split_range_and_weight_product(): - n = 16 - marginals = [ - Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0), - Uniform( - DummySampler(2), - lower_bound=[20.0, 30.0], - upper_bound=[24.0, 36.0], - ), - ] - tm = ProductMeasure(sampler=DigitalNetB2(3, seed=41), marginals=marginals) - - u = tm.discrete_distrib.gen_samples(n) - x = tm._transform(u) - x_call, jac = tm(n, return_weights=True) - expected = np.concatenate( - [ - marginals[0]._jacobian_transform_r(u[..., :1], return_weights=False), - marginals[1]._jacobian_transform_r(u[..., 1:], return_weights=False), - ], - axis=-1, - ) - expected_range = np.array([[10.0, 12.0], [20.0, 24.0], [30.0, 36.0]]) - - assert x.shape == (n, 3) - assert np.allclose(tm.range, expected_range) - assert np.allclose(x, expected) - assert np.all((10.0 <= x[:, 0]) & (x[:, 0] <= 12.0)) - assert np.all((20.0 <= x[:, 1]) & (x[:, 1] <= 24.0)) - assert np.all((30.0 <= x[:, 2]) & (x[:, 2] <= 36.0)) - assert np.allclose(tm._weight(x), 1.0 / (2.0 * 4.0 * 6.0)) - assert x_call.shape == (n, 3) - assert np.allclose(jac, 2.0 * 4.0 * 6.0) - - -def test_product_measure_invalid_inputs(): - with pytest.raises(ParameterError, match="nonempty list of marginals"): - ProductMeasure(sampler=DigitalNetB2(1, seed=7), marginals=[]) - - with pytest.raises(ParameterError, match="marginal"): - ProductMeasure(sampler=DigitalNetB2(1, seed=7), marginals=[object()]) - - with pytest.raises(ParameterError, match="AbstractDiscreteDistribution"): - ProductMeasure(sampler=object(), marginals=[Uniform(DummySampler(1))]) - - marginals = [Uniform(DummySampler(1))] - with pytest.raises(DimensionError, match="sum of marginal dimensions"): - ProductMeasure(sampler=DigitalNetB2(2, seed=7), marginals=marginals) - - -def test_product_measure_rejects_non_dimension_preserving_marginal(): - marginal = AcceptanceRejection( - DigitalNetB2(2, seed=7), - lambda x: np.ones(len(x)), - 1.0, - 1.0, - ) - - with pytest.raises(DimensionError, match="dimension-preserving"): - ProductMeasure(DigitalNetB2(2, seed=11), [marginal]) - - -def test_product_measure_spawn_preserves_marginal_blocks_and_replaces_outer_sampler(): - marginals = [ - Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0), - Uniform( - DummySampler(2), - lower_bound=[20.0, 30.0], - upper_bound=[24.0, 36.0], - ), - ] - tm = ProductMeasure(sampler=DigitalNetB2(3, seed=41), marginals=marginals) + ] + tm = ProductMeasure(sampler=DigitalNetB2(3, seed=41), marginals=marginals) + + spawn = tm.spawn(s=1)[0] + + self.assertIsInstance(spawn, ProductMeasure) + self.assertEqual(spawn.d, 3) + self.assertEqual(spawn.marginals, tm.marginals) + self.assertIsNot(spawn.discrete_distrib, tm.discrete_distrib) + np.testing.assert_array_equal(spawn.marginal_dimensions, np.array([1, 2])) + + with self.assertRaises(DimensionError): + tm.spawn(s=1, dimensions=4) + + + def test_does_not_use_marginal_dummy_sampler_values(self): + marginals = [ + Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0), + Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0), + ] + + with self.assertRaisesRegex(ParameterError, "construction placeholder"): + marginals[0].discrete_distrib(4) + + tm = ProductMeasure(sampler=DigitalNetB2(2, seed=19), marginals=marginals) + x = tm(8) + + self.assertEqual(x.shape, (8, 2)) + self.assertTrue(np.all((0.0 <= x[:, 0]) & (x[:, 0] <= 2.0))) + self.assertTrue(np.all((10.0 <= x[:, 1]) & (x[:, 1] <= 12.0))) + + + def test_same_outer_seed_matches_different_outer_seed_changes(self): + marginals = [ + Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0), + Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0), + ] + + first = ProductMeasure(sampler=DigitalNetB2(2, seed=101), marginals=marginals)(16) + same_outer = ProductMeasure(sampler=DigitalNetB2(2, seed=101), marginals=marginals)(16) + different_outer = ProductMeasure(sampler=DigitalNetB2(2, seed=102), marginals=marginals)(16) + + np.testing.assert_array_equal(first, same_outer) + self.assertFalse(np.array_equal(first, different_outer)) - spawn = tm.spawn(s=1)[0] - - assert isinstance(spawn, ProductMeasure) - assert spawn.d == 3 - assert spawn.marginals == tm.marginals - assert spawn.discrete_distrib is not tm.discrete_distrib - assert np.array_equal(spawn.marginal_dimensions, np.array([1, 2])) - with pytest.raises(DimensionError): - tm.spawn(s=1, dimensions=4) - - -def test_product_measure_does_not_use_marginal_dummy_sampler_values(): - marginals = [ - Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0), - Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0), - ] - - with pytest.raises(ParameterError, match="construction placeholder"): - marginals[0].discrete_distrib(4) - - tm = ProductMeasure(sampler=DigitalNetB2(2, seed=19), marginals=marginals) - x = tm(8) - - assert x.shape == (8, 2) - assert np.all((0.0 <= x[:, 0]) & (x[:, 0] <= 2.0)) - assert np.all((10.0 <= x[:, 1]) & (x[:, 1] <= 12.0)) - - -def test_product_measure_same_outer_seed_matches_different_outer_seed_changes(): - marginals = [ - Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0), - Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0), - ] - - first = ProductMeasure(sampler=DigitalNetB2(2, seed=101), marginals=marginals)(16) - same_outer = ProductMeasure(sampler=DigitalNetB2(2, seed=101), marginals=marginals)(16) - different_outer = ProductMeasure(sampler=DigitalNetB2(2, seed=102), marginals=marginals)(16) - - assert np.array_equal(first, same_outer) - assert not np.array_equal(first, different_outer) - - -def test_product_measure_replication_means_close_to_uniform_targets(): - n = 1024 - r = 4 - marginals = [ - Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0), - Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0), - ] - tm = ProductMeasure( - sampler=DigitalNetB2(2, seed=101, replications=r), - marginals=marginals, - ) - - x = tm(n) - replication_means = x.mean(axis=1) - - assert x.shape == (r, n, 2) - assert np.allclose(replication_means[:, 0], 1.0, atol=0.03) - assert np.allclose(replication_means[:, 1], 11.0, atol=0.03) - - -def test_product_measure_with_scipywrapper_beta_marginal(): - n = 64 - marginals = [ - Uniform(DummySampler(1), lower_bound=-1.0, upper_bound=1.0), - SciPyWrapper(DummySampler(1), stats.beta(a=2.0, b=5.0)), - ] - tm = ProductMeasure(sampler=DigitalNetB2(2, seed=71), marginals=marginals) - - x = tm(n) - - assert x.shape == (n, 2) - assert np.all((-1.0 <= x[:, 0]) & (x[:, 0] <= 1.0)) - assert np.all((0.0 <= x[:, 1]) & (x[:, 1] <= 1.0)) - - -def test_product_measure_matches_equivalent_scipywrapper(): - n = 128 - seed = 55 - scipy_marginals = [stats.norm(loc=0.0, scale=1.0), stats.gamma(a=2.0, scale=1.0)] - product_marginals = [ - SciPyWrapper(DummySampler(1), scipy_marginals[0]), - SciPyWrapper(DummySampler(1), scipy_marginals[1]), - ] - - product_samples = ProductMeasure( - sampler=DigitalNetB2(2, seed=seed), - marginals=product_marginals, - )(n) - scipy_samples = SciPyWrapper(DigitalNetB2(2, seed=seed), scipy_marginals)(n) - - assert np.array_equal(product_samples, scipy_samples) - - -def test_product_measure_with_gaussian_copula_marginal(): - n = 64 - copula = GaussianCopula( - DummySampler(2), - marginals=[stats.beta(a=2.0, b=5.0), stats.gamma(a=3.0, scale=1.0)], - correlation=[[1.0, 0.5], [0.5, 1.0]], - ) - marginals = [copula, Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0)] - tm = ProductMeasure(sampler=DigitalNetB2(3, seed=81), marginals=marginals) - - x = tm(n) - - assert x.shape == (n, 3) - assert np.all((0.0 <= x[:, 0]) & (x[:, 0] <= 1.0)) - assert np.all(x[:, 1] >= 0.0) - assert np.all((10.0 <= x[:, 2]) & (x[:, 2] <= 12.0)) - - -def test_product_measure_recursive_transform_sampling_supported_but_weights_restricted(): - recursive_marginal = Uniform( - Uniform(DummySampler(1), lower_bound=0.0, upper_bound=1.0), - lower_bound=2.0, - upper_bound=4.0, - ) - direct_marginal = Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0) - tm = ProductMeasure( - sampler=DigitalNetB2(2, seed=91), - marginals=[recursive_marginal, direct_marginal], - ) - - x = tm(16) - - assert x.shape == (16, 2) - assert np.all((2.0 <= x[:, 0]) & (x[:, 0] <= 4.0)) - assert np.all((10.0 <= x[:, 1]) & (x[:, 1] <= 12.0)) - with pytest.raises(ParameterError, match="direct marginal"): - tm(16, return_weights=True) + def test_replication_means_close_to_uniform_targets(self): + n = 1024 + r = 4 + marginals = [ + Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0), + Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0), + ] + tm = ProductMeasure( + sampler=DigitalNetB2(2, seed=101, replications=r), + marginals=marginals, + ) + + x = tm(n) + replication_means = x.mean(axis=1) + + self.assertEqual(x.shape, (r, n, 2)) + np.testing.assert_allclose(replication_means[:, 0], 1.0, atol=0.03) + np.testing.assert_allclose(replication_means[:, 1], 11.0, atol=0.03) + + + def test_with_scipywrapper_beta_marginal(self): + n = 64 + marginals = [ + Uniform(DummySampler(1), lower_bound=-1.0, upper_bound=1.0), + SciPyWrapper(DummySampler(1), stats.beta(a=2.0, b=5.0)), + ] + tm = ProductMeasure(sampler=DigitalNetB2(2, seed=71), marginals=marginals) + + x = tm(n) + + self.assertEqual(x.shape, (n, 2)) + self.assertTrue(np.all((-1.0 <= x[:, 0]) & (x[:, 0] <= 1.0))) + self.assertTrue(np.all((0.0 <= x[:, 1]) & (x[:, 1] <= 1.0))) + + + def test_matches_equivalent_scipywrapper(self): + n = 128 + seed = 55 + scipy_marginals = [stats.norm(loc=0.0, scale=1.0), stats.gamma(a=2.0, scale=1.0)] + product_marginals = [ + SciPyWrapper(DummySampler(1), scipy_marginals[0]), + SciPyWrapper(DummySampler(1), scipy_marginals[1]), + ] + + product_samples = ProductMeasure( + sampler=DigitalNetB2(2, seed=seed), + marginals=product_marginals, + )(n) + scipy_samples = SciPyWrapper(DigitalNetB2(2, seed=seed), scipy_marginals)(n) + + np.testing.assert_array_equal(product_samples, scipy_samples) + + + def test_with_gaussian_copula_marginal(self): + n = 64 + copula = GaussianCopula( + DummySampler(2), + marginals=[stats.beta(a=2.0, b=5.0), stats.gamma(a=3.0, scale=1.0)], + correlation=[[1.0, 0.5], [0.5, 1.0]], + ) + marginals = [copula, Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0)] + tm = ProductMeasure(sampler=DigitalNetB2(3, seed=81), marginals=marginals) + + x = tm(n) + + self.assertEqual(x.shape, (n, 3)) + self.assertTrue(np.all((0.0 <= x[:, 0]) & (x[:, 0] <= 1.0))) + self.assertTrue(np.all(x[:, 1] >= 0.0)) + self.assertTrue(np.all((10.0 <= x[:, 2]) & (x[:, 2] <= 12.0))) + + + def test_recursive_transform_sampling_supported_but_weights_restricted(self): + recursive_marginal = Uniform( + Uniform(DummySampler(1), lower_bound=0.0, upper_bound=1.0), + lower_bound=2.0, + upper_bound=4.0, + ) + direct_marginal = Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0) + tm = ProductMeasure( + sampler=DigitalNetB2(2, seed=91), + marginals=[recursive_marginal, direct_marginal], + ) + + x = tm(16) + + self.assertEqual(x.shape, (16, 2)) + self.assertTrue(np.all((2.0 <= x[:, 0]) & (x[:, 0] <= 4.0))) + self.assertTrue(np.all((10.0 <= x[:, 1]) & (x[:, 1] <= 12.0))) + with self.assertRaisesRegex(ParameterError, "direct marginal"): + tm(16, return_weights=True)