From f67b86537916fd5b8f8ca1b36f11b74eb290178a Mon Sep 17 00:00:00 2001 From: kat Date: Thu, 28 Sep 2023 14:27:02 +0200 Subject: [PATCH] Add tests for integrands with Gaussian measure (#585) * rename var -> std in uniform_to_gaussian_quadprob * use pytest-cases --- src/probnum/problems/zoo/quad/__init__.py | 20 +- .../zoo/quad/_quadproblems_gaussian.py | 74 ++++-- .../zoo/quad/_quadproblems_uniform.py | 2 +- .../test_quad/test_quadproblems_gaussian.py | 239 ++++++++++++++++++ .../test_quadproblems_gaussian_cases.py | 140 ++++++++++ .../test_quad/test_quadproblems_uniform.py | 151 +++-------- .../test_quadproblems_uniform_cases.py | 111 ++++++++ 7 files changed, 599 insertions(+), 138 deletions(-) create mode 100644 tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian.py create mode 100644 tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian_cases.py create mode 100644 tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform_cases.py diff --git a/src/probnum/problems/zoo/quad/__init__.py b/src/probnum/problems/zoo/quad/__init__.py index b5e8e213b..16320cb2c 100644 --- a/src/probnum/problems/zoo/quad/__init__.py +++ b/src/probnum/problems/zoo/quad/__init__.py @@ -5,4 +5,22 @@ from ._quadproblems_uniform import * # Public classes and functions. Order is reflected in documentation. -__all__ = ["circulargaussian2d", "hennig1d", "hennig2d", "sombrero2d"] +__all__ = [ + "circulargaussian2d", + "hennig1d", + "hennig2d", + "sombrero2d", + "uniform_to_gaussian_quadprob", + "sum_polynomials", + "genz_continuous", + "genz_cornerpeak", + "genz_discontinuous", + "genz_gaussian", + "genz_oscillatory", + "genz_productpeak", + "bratley1992", + "roos_arnold", + "gfunction", + "morokoff_caflisch_1", + "morokoff_caflisch_2", +] diff --git a/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py b/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py index 0d4a764d3..aa300a9b3 100644 --- a/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py +++ b/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py @@ -1,6 +1,6 @@ """Test problems for integration against a Gaussian measure.""" -from typing import Callable +from typing import Callable, Union import numpy as np from scipy import special @@ -16,8 +16,35 @@ ] +# Construct transformation of the integrand +def uniform_to_gaussian_integrand( + fun: Callable[[np.ndarray], np.ndarray], + mean: Union[float, np.floating, np.ndarray] = 0.0, + std: Union[float, np.floating, np.ndarray] = 1.0, +) -> Callable[[np.ndarray], np.ndarray]: + # mean and var should be either one-dimensional + if isinstance(mean, np.ndarray): + if len(mean.shape) != 1: + raise TypeError( + "The mean parameter should be a float or a d-dimensional array." + ) + + if isinstance(std, np.ndarray): + if len(std.shape) != 1: + raise TypeError( + "The std parameter should be a float or a d-dimensional array." + ) + + def new_func(x): + return fun(norm.cdf(x, loc=mean, scale=std)) + + return new_func + + def uniform_to_gaussian_quadprob( - quadprob: QuadratureProblem, mean: FloatLike = 0.0, var: FloatLike = 1.0 + quadprob: QuadratureProblem, + mean: Union[float, np.floating, np.ndarray] = 0.0, + std: Union[float, np.floating, np.ndarray] = 1.0, ) -> QuadratureProblem: r"""Creates a new QuadratureProblem for integration against a Gaussian on :math:`\mathbb{R}^d` by using an existing QuadratureProblem whose integrand is @@ -39,9 +66,12 @@ def uniform_to_gaussian_quadprob( quadprob A QuadratureProblem instance which includes an integrand defined on [0,1]^d mean - Mean of the Gaussian distribution. - var - Diagonal element for the covariance matrix of the Gaussian distribution. + Mean of the Gaussian distribution. If `float`, mean is set to the same value + across all dimensions. Else, specifies the mean as a d-dimensional array. + std + Diagonal element for the covariance matrix of the Gaussian distribution. If + `float`, the covariance matrix has the same diagonal value for all dimensions. + Else, specifies the covariance matrix via a d-dimensional array. Returns ------- @@ -78,27 +108,14 @@ def uniform_to_gaussian_quadprob( if np.ndim(quadprob.solution) != 0: raise ValueError("The solution of quadprob is not a scalar.") - # Construct transformation of the integrand - def uniform_to_gaussian_integrand( - fun: Callable[[np.ndarray], np.ndarray], - mean: FloatLike = 0.0, - var: FloatLike = 1.0, - ) -> Callable[[np.ndarray], np.ndarray]: - # mean and var should be either one-dimensional, or an array of dimension d - if isinstance(mean, float) is False: - raise TypeError("The mean parameter should be a float.") - - if isinstance(var, float) is False or var <= 0.0: - raise TypeError("The variance should be a positive float.") - - def newfunc(x): - return fun(norm.cdf((x - mean) / var)) + dim = lower_bd.shape[0] - return newfunc - - gaussian_measure = GaussianMeasure(mean=mean, cov=var) + cov = std**2 + if isinstance(std, np.ndarray): + cov = np.eye(dim) * cov + gaussian_measure = GaussianMeasure(mean=mean, cov=cov, input_dim=dim) return QuadratureProblem( - fun=uniform_to_gaussian_integrand(fun=quadprob.fun, mean=mean, var=var), + fun=uniform_to_gaussian_integrand(fun=quadprob.fun, mean=mean, std=std), measure=gaussian_measure, solution=quadprob.solution, ) @@ -159,7 +176,7 @@ def sum_polynomials( if len(b.shape) != 2: raise ValueError( - f"Invalid shape {a.shape} for parameter `a`. " + f"Invalid shape {b.shape} for parameter `b`. " f"Expected parameters of shape (p+1)xdim" ) @@ -196,9 +213,12 @@ def integrand(x: np.ndarray) -> np.ndarray: delta = (np.remainder(b, 2) - 1) ** 2 doublefact = special.factorial2(b - 1) solution = np.sum(np.prod(a * delta * (var**b) * doublefact, axis=1)) - mean = np.zeros_like(dim) + if isinstance(var, float): + mean = 0.0 + else: + mean = np.zeros(dim) - gaussian_measure = GaussianMeasure(mean=mean, cov=var) + gaussian_measure = GaussianMeasure(mean=mean, cov=var, input_dim=dim) return QuadratureProblem( fun=integrand, measure=gaussian_measure, diff --git a/src/probnum/problems/zoo/quad/_quadproblems_uniform.py b/src/probnum/problems/zoo/quad/_quadproblems_uniform.py index 250452910..7883e49d5 100644 --- a/src/probnum/problems/zoo/quad/_quadproblems_uniform.py +++ b/src/probnum/problems/zoo/quad/_quadproblems_uniform.py @@ -264,7 +264,7 @@ def fun(x: np.ndarray) -> np.ndarray: return f.reshape(n, 1) if dim == 1: - solution = (np.exp(a * u) - 1.0) / a + solution = ((np.exp(a * u) - 1.0) / a).item() if dim > 1: solution = np.prod((np.exp(a * np.minimum(u, 1.0)) - 1.0) / a) diff --git a/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian.py b/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian.py new file mode 100644 index 000000000..df932036b --- /dev/null +++ b/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian.py @@ -0,0 +1,239 @@ +from itertools import product +from typing import Callable, Optional, Union + +import numpy as np +import pytest +from pytest_cases import parametrize_with_cases + +from probnum.problems import QuadratureProblem +from probnum.problems.zoo.quad import ( + bratley1992, + genz_continuous, + genz_cornerpeak, + genz_discontinuous, + genz_gaussian, + genz_oscillatory, + genz_productpeak, + gfunction, + morokoff_caflisch_1, + morokoff_caflisch_2, + roos_arnold, + sum_polynomials, + uniform_to_gaussian_quadprob, +) +from probnum.quad.integration_measures import LebesgueMeasure +from probnum.quad.typing import DomainLike +from tests.test_problems.test_zoo.test_quad.test_quadproblems_gaussian_cases import ( + GenzStandardNormalCases, + GenzVariedNormalCases, + OtherIntegrandsGaussianCases, +) + +dim_values = [1, 2] +quad_prob_constructor_values = [ + bratley1992, + genz_continuous, + genz_cornerpeak, + genz_discontinuous, + genz_gaussian, + genz_oscillatory, + genz_productpeak, + gfunction, + morokoff_caflisch_1, + morokoff_caflisch_2, + roos_arnold, +] +param_combinations = list(product(dim_values, quad_prob_constructor_values)) + + +@pytest.mark.parametrize("quad_prob_constructor", quad_prob_constructor_values) +@pytest.mark.parametrize("dim", dim_values) +def test_wrapping_all_test_functions_works( + dim: int, + quad_prob_constructor: Callable[[int], QuadratureProblem], +): + """Integration test that wrapping all problems works.""" + quadprob = quad_prob_constructor(dim) + gaussian_quadprob = uniform_to_gaussian_quadprob(quadprob) + assert isinstance(gaussian_quadprob, QuadratureProblem) + + +@pytest.mark.parametrize( + "domain, dim", + [ + ((0.5, 1.0), 1), + ((0.5, 1.0), 2), + ((0.0, 1.5), 1), + ((0.0, 1.5), 2), + ], +) +def test_wrong_measure_bounds_scalar_raises(domain: DomainLike, dim: int): + uniform_measure = LebesgueMeasure(domain=domain, input_dim=dim) + quadprob = genz_continuous(dim) + quadprob.measure = uniform_measure + with pytest.raises(ValueError) as exc_info: + uniform_to_gaussian_quadprob(quadprob) + + assert "[0,1]" in str(exc_info.value) + + +@pytest.mark.parametrize( + "domain", + [ + (np.array([0.0, 0.5]), np.array([1.0, 1.0])), + (np.array([0.0, 0.0]), np.array([1.0, 1.5])), + ], +) +def test_wrong_measure_bounds_array_raises(domain: DomainLike): + lower_bd, _ = domain + dim = len(lower_bd) + uniform_measure = LebesgueMeasure(domain=domain, input_dim=dim) + quadprob = genz_continuous(dim) + quadprob.measure = uniform_measure + with pytest.raises(ValueError) as exc_info: + uniform_to_gaussian_quadprob(quadprob) + + assert "[0,1]" in str(exc_info.value) + + +def test_wrong_data_shapes_for_mean_raises(): + with pytest.raises(TypeError) as exc_info: + uniform_to_gaussian_quadprob( + genz_continuous(2), mean=np.array([[0.0, 0.5]]), std=np.array([0.0, 0.5]) + ) + + assert "mean parameter" in str(exc_info.value) + + +def test_wrong_data_shapes_for_std_raises(): + with pytest.raises(TypeError) as exc_info: + uniform_to_gaussian_quadprob( + genz_continuous(2), mean=np.array([0.0, 0.5]), std=np.array([[0.0, 0.5]]) + ) + + assert "std parameter" in str(exc_info.value) + + +def test_wrong_shaped_a_in_sum_polynomials_raises(): + a = np.array([1.0, 1.0]) + dim = 2 + with pytest.raises(ValueError) as exc_info: + sum_polynomials(dim=dim, a=a) + + assert "Invalid shape" in str(exc_info.value) + assert "parameter `a`" in str(exc_info.value) + + +def test_wrong_shaped_b_in_sum_polynomials_raises(): + b = np.array([1.0, 1.0]) + dim = 2 + with pytest.raises(ValueError) as exc_info: + sum_polynomials(dim=dim, b=b) + + assert "Invalid shape" in str(exc_info.value) + assert "parameter `b`" in str(exc_info.value) + + +def test_mismatch_dim_a_in_sum_polynomials_raises(): + a = np.array([[1.0, 1.0]]) + dim = 1 + with pytest.raises(ValueError) as exc_info: + sum_polynomials(dim=dim, a=a) + + assert "parameter `a`" in str(exc_info.value) + assert f"Expected {dim} columns" in str(exc_info.value) + + +def test_mismatch_dim_b_in_sum_polynomials_raises(): + b = np.array([[1.0, 1.0]]) + dim = 1 + with pytest.raises(ValueError) as exc_info: + sum_polynomials(dim=dim, b=b) + + assert "parameter `b`" in str(exc_info.value) + assert f"Expected {dim} columns" in str(exc_info.value) + + +def test_negative_values_for_b_in_sum_polynomials_raises(): + b = np.array([[-0.5, 1.0]]) + dim = 2 + with pytest.raises(ValueError) as exc_info: + sum_polynomials(dim=dim, b=b) + + assert "negative" in str(exc_info.value) + assert "parameters `b`" in str(exc_info.value) + + +@parametrize_with_cases( + "quadprob, rtol", + cases=( + GenzStandardNormalCases, + GenzVariedNormalCases, + OtherIntegrandsGaussianCases, + ), +) +def test_genz_gaussian(quadprob: QuadratureProblem, rtol: float): + """Compare a Monte Carlo estimator against Gaussian measure with a very large number + of samples to the true value of the integral. + + The former should be an approximation of the latter. + """ + # Number of Monte Carlo samples for the test + n = 10000000 + + # generate some normally distributed points from N(0, 1) + rng = np.random.default_rng(0) + x_gaussian = quadprob.measure.sample(n, rng=rng) # Monte Carlo samples x_1,...,x_n + + # Test that all Monte Carlo estimators approximate the true value of the integral + # (integration against Gaussian) + np.testing.assert_allclose( + np.sum(quadprob.fun(x_gaussian)) / n, + quadprob.solution, + rtol=rtol, + ) + + +@pytest.mark.parametrize( + "dim, a, b, var", + [ + (1, None, None, 1.0), + (2, None, None, 1.0), + (1, np.array([[0.5]]), np.array([[1]]), 0.5), + (1, np.array([[0.5]]), np.array([[1]]), np.array([[0.5]])), + ( + 2, + np.array([[0.5, 1.0], [2.0, 2.0]]), + np.array([[1, 2], [2, 3]]), + np.array([[1.0, 0], [0, 0.5]]), + ), + ], +) +def test_sum_polynomials( + dim: int, + a: Optional[np.ndarray], + b: Optional[np.ndarray], + var: Union[float, np.ndarray], +): + """Compare a Monte Carlo estimator against Gaussian measure with a very large number + of samples to the true value of the integral. + + The former should be an approximation of the latter. + """ + + # Number of Monte Carlo samples for the test + n = 10000000 + + quadprob = sum_polynomials(dim, a, b, var) + + # generate some normally distributed points from N(0, 1) + rng = np.random.default_rng(0) + x_gaussian = quadprob.measure.sample(n, rng=rng) # Monte Carlo samples x_1,...,x_n + + # Test that all Monte Carlo estimators approximate the true value of the integral + # (integration against uniform) + np.testing.assert_allclose( + np.sum(quadprob.fun(x_gaussian)) / n, + quadprob.solution, + atol=1e-02, + ) diff --git a/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian_cases.py b/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian_cases.py new file mode 100644 index 000000000..4061684d4 --- /dev/null +++ b/tests/test_problems/test_zoo/test_quad/test_quadproblems_gaussian_cases.py @@ -0,0 +1,140 @@ +import numpy as np +import pytest + +from probnum.problems.zoo.quad import ( + bratley1992, + genz_continuous, + genz_cornerpeak, + genz_discontinuous, + genz_gaussian, + genz_oscillatory, + genz_productpeak, + gfunction, + morokoff_caflisch_1, + roos_arnold, + uniform_to_gaussian_quadprob, +) + +genz_params = [ + (None, None, 3), + (5, 0.5, 1), + (5.0, 1.0, 1), + (2.0, 0.8, 1), + (5, 0.5, 4), + (5.0, 1.0, 4), + (2.0, 0.8, 4), +] + + +class GenzStandardNormalCases: + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_continuous(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = uniform_to_gaussian_quadprob( + genz_continuous(dim=dim, a=a_vec, u=u_vec) + ) + rtol = 1e-2 + return quadprob, rtol + + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_cornerpeak(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = uniform_to_gaussian_quadprob( + genz_cornerpeak(dim=dim, a=a_vec, u=u_vec) + ) + rtol = 1e-2 + return quadprob, rtol + + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_discontinuous(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = uniform_to_gaussian_quadprob( + genz_discontinuous(dim=dim, a=a_vec, u=u_vec) + ) + rtol = 1e-2 + return quadprob, rtol + + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_gaussian(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = uniform_to_gaussian_quadprob( + genz_gaussian(dim=dim, a=a_vec, u=u_vec) + ) + rtol = 1e-2 + return quadprob, rtol + + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_oscillatory(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = uniform_to_gaussian_quadprob( + genz_oscillatory(dim=dim, a=a_vec, u=u_vec) + ) + rtol = 9e-2 + return quadprob, rtol + + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_productpeak(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = uniform_to_gaussian_quadprob( + genz_productpeak(dim=dim, a=a_vec, u=u_vec) + ) + rtol = 1e-2 + return quadprob, rtol + + +normal_params = [ + (3.0, 2.0, 2), + (0.5, 2.0, 1), + (np.array([0.0, 1.0]), np.array([0.5, 2.0]), 2), +] + + +class GenzVariedNormalCases: + @pytest.mark.parametrize("mean, std, dim", normal_params) + def case_genz_continuous(self, mean, std, dim): + quadprob = uniform_to_gaussian_quadprob( + genz_continuous(dim=dim), mean=mean, std=std + ) + rtol = 1e-2 + return quadprob, rtol + + +dim_params = [1, 2] + + +class OtherIntegrandsGaussianCases: + @pytest.mark.parametrize("dim", dim_params) + def case_bratley1992(self, dim): + quadprob = uniform_to_gaussian_quadprob(bratley1992(dim=dim)) + rtol = 1e-3 + return quadprob, rtol + + @pytest.mark.parametrize("dim", dim_params) + def case_roos_arnold(self, dim): + quadprob = uniform_to_gaussian_quadprob(roos_arnold(dim=dim)) + rtol = 1e-3 + return quadprob, rtol + + @pytest.mark.parametrize("dim", dim_params) + def case_gfunction(self, dim): + quadprob = uniform_to_gaussian_quadprob(gfunction(dim=dim)) + rtol = 1e-3 + return quadprob, rtol + + @pytest.mark.parametrize("dim", dim_params) + def case_morokoff_caflisch_1(self, dim): + quadprob = uniform_to_gaussian_quadprob(morokoff_caflisch_1(dim=dim)) + rtol = 1e-3 + return quadprob, rtol + + @pytest.mark.parametrize("dim", dim_params) + def case_morokoff_caflisch_2(self, dim): + quadprob = uniform_to_gaussian_quadprob(morokoff_caflisch_1(dim=dim)) + rtol = 1e-3 + return quadprob, rtol diff --git a/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform.py b/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform.py index 1d1176634..424959796 100644 --- a/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform.py +++ b/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform.py @@ -1,7 +1,9 @@ +from itertools import product from typing import Callable import numpy as np import pytest +from pytest_cases import parametrize_with_cases from probnum.problems import QuadratureProblem from probnum.problems.zoo.quad import ( @@ -17,6 +19,10 @@ morokoff_caflisch_2, roos_arnold, ) +from tests.test_problems.test_zoo.test_quad.test_quadproblems_uniform_cases import ( + GenzUniformCases, + OtherIntegrandsUniformCases, +) @pytest.mark.parametrize( @@ -72,83 +78,39 @@ def test_integrand_eval_checks( quad_problem.fun(np.full((4, 2), -0.1)) -@pytest.mark.parametrize( - "a, u, dim", - [ - (None, None, 3), - (5, 0.5, 1), - (5.0, 1.0, 1), - (2.0, 0.8, 1), - (5, 0.5, 4), - (5.0, 1.0, 4), - (2.0, 0.8, 4), - ], -) -def test_genz_uniform(a, u, dim): - """Compare a Monte Carlo estimator with a very large number of samples to the true - value of the integral. - - The former should be an approximation of the latter. - """ - - # Number of Monte Carlo samples for the test - n = 10000000 - - # Define parameters a and u - a_vec = np.repeat(a, dim) if a is not None else None - u_vec = np.repeat(u, dim) if u is not None else None - - quadprob_genz_continuous = genz_continuous(dim=dim, a=a_vec, u=u_vec) - quadprob_genz_cornerpeak = genz_cornerpeak(dim=dim, a=a_vec, u=u_vec) - quadprob_genz_discontinuous = genz_discontinuous(dim=dim, a=a_vec, u=u_vec) - quadprob_genz_gaussian = genz_gaussian(dim=dim, a=a_vec, u=u_vec) - quadprob_genz_oscillatory = genz_oscillatory(dim=dim, a=a_vec, u=u_vec) - quadprob_genz_productpeak = genz_productpeak(dim=dim, a=a_vec, u=u_vec) +dim_values = [1, 2] +quad_prob_constructor_values = [ + bratley1992, + genz_continuous, + genz_cornerpeak, + genz_discontinuous, + genz_gaussian, + genz_oscillatory, + genz_productpeak, + gfunction, + morokoff_caflisch_1, + morokoff_caflisch_2, + roos_arnold, +] +param_combinations = list(product(dim_values, quad_prob_constructor_values)) - # generate some uniformly distributed points on [0,1]^d - np.random.seed(0) - x_unif = np.random.uniform( - low=0.0, high=1.0, size=(n, dim) - ) # Monte Carlo samples x_1,...,x_n - # Test that all Monte Carlo estimators approximate the true value of the integral - # (integration against uniform) - np.testing.assert_allclose( - np.sum(quadprob_genz_continuous.fun(x_unif)) / n, - quadprob_genz_continuous.solution, - rtol=1e-03, - ) - np.testing.assert_allclose( - np.sum(quadprob_genz_cornerpeak.fun(x_unif)) / n, - quadprob_genz_cornerpeak.solution, - rtol=3e-03, - ) - np.testing.assert_allclose( - np.sum(quadprob_genz_discontinuous.fun(x_unif)) / n, - quadprob_genz_discontinuous.solution, - rtol=1e-03, - ) - np.testing.assert_allclose( - np.sum(quadprob_genz_gaussian.fun(x_unif)) / n, - quadprob_genz_gaussian.solution, - rtol=2e-03, - ) - np.testing.assert_allclose( - np.sum(quadprob_genz_oscillatory.fun(x_unif)) / n, - quadprob_genz_oscillatory.solution, - rtol=9e-02, - ) - np.testing.assert_allclose( - np.sum(quadprob_genz_productpeak.fun(x_unif)) / n, - quadprob_genz_productpeak.solution, - rtol=1e-03, - ) +@pytest.mark.parametrize("dim, quad_prob_constructor", param_combinations) +def test_integrand_solution_float( + dim: int, quad_prob_constructor: Callable[..., QuadratureProblem] +): + quadprob = quad_prob_constructor(dim) + if np.ndim(quadprob.solution) != 0: + raise ValueError(f"The solution of {quadprob} is not a scalar.") -@pytest.mark.parametrize("dim", [(1), (4), (10)]) -def test_integrands_uniform(dim): - """Compare a Monte Carlo estimator with a very large number of samples to the true - value of the integral. +@parametrize_with_cases( + "quadprob, rtol", + cases=(GenzUniformCases, OtherIntegrandsUniformCases), +) +def test_quadprob_uniform_with_mc(quadprob, rtol): + """Compare a Monte Carlo estimator against uniform measure with a very large number + of samples to the true value of the integral. The former should be an approximation of the latter. """ @@ -156,43 +118,14 @@ def test_integrands_uniform(dim): # Number of Monte Carlo samples for the test n = 10000000 - # generate some uniformly distributed points on [0,1]^d - np.random.seed(0) - x_unif = np.random.uniform( - low=0.0, high=1.0, size=(n, dim) - ) # Monte Carlo samples x_1,...,x_n - - # Set the integrands to test - quadprob_bratley1992 = bratley1992(dim=dim) - quadprob_roos_arnold = roos_arnold(dim=dim) - quadprob_gfunction = gfunction(dim=dim) - quadprob_morokoff_caflisch_1 = morokoff_caflisch_1(dim=dim) - quadprob_morokoff_caflisch_2 = morokoff_caflisch_2(dim=dim) + # generate some uniformly distributed points from [0, 1] + rng = np.random.default_rng(0) + x_uniform = quadprob.measure.sample(n, rng=rng) # Monte Carlo samples x_1,...,x_n # Test that all Monte Carlo estimators approximate the true value of the integral - # (integration against uniform) - np.testing.assert_allclose( - np.sum(quadprob_bratley1992.fun(x_unif)) / n, - quadprob_bratley1992.solution, - rtol=1e-03, - ) - np.testing.assert_allclose( - np.sum(quadprob_roos_arnold.fun(x_unif)) / n, - quadprob_roos_arnold.solution, - rtol=1e-03, - ) - np.testing.assert_allclose( - np.sum(quadprob_gfunction.fun(x_unif)) / n, - quadprob_gfunction.solution, - rtol=1e-03, - ) - np.testing.assert_allclose( - np.sum(quadprob_morokoff_caflisch_1.fun(x_unif)) / n, - quadprob_morokoff_caflisch_1.solution, - rtol=1e-03, - ) + # (integration against [0, 1]) np.testing.assert_allclose( - np.sum(quadprob_morokoff_caflisch_2.fun(x_unif)) / n, - quadprob_morokoff_caflisch_2.solution, - rtol=1e-03, + np.sum(quadprob.fun(x_uniform)) / n, + quadprob.solution, + rtol=rtol, ) diff --git a/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform_cases.py b/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform_cases.py new file mode 100644 index 000000000..ec40e094f --- /dev/null +++ b/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform_cases.py @@ -0,0 +1,111 @@ +import numpy as np +import pytest + +from probnum.problems.zoo.quad import ( + bratley1992, + genz_continuous, + genz_cornerpeak, + genz_discontinuous, + genz_gaussian, + genz_oscillatory, + genz_productpeak, + gfunction, + morokoff_caflisch_1, + morokoff_caflisch_2, + roos_arnold, +) + +genz_params = [ + (None, None, 3), + (5, 0.5, 1), + (5.0, 1.0, 1), + (2.0, 0.8, 1), + (5, 0.5, 4), + (5.0, 1.0, 4), + (2.0, 0.8, 4), +] + + +class GenzUniformCases: + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_continuous(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = genz_continuous(dim=dim, a=a_vec, u=u_vec) + rtol = 1e-3 + return quadprob, rtol + + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_cornerpeak(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = genz_cornerpeak(dim=dim, a=a_vec, u=u_vec) + rtol = 1e-2 + return quadprob, rtol + + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_discontinuous(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = genz_discontinuous(dim=dim, a=a_vec, u=u_vec) + rtol = 2e-3 + return quadprob, rtol + + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_gaussian(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = genz_gaussian(dim=dim, a=a_vec, u=u_vec) + rtol = 1e-2 + return quadprob, rtol + + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_oscillatory(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = genz_oscillatory(dim=dim, a=a_vec, u=u_vec) + rtol = 9e-2 + return quadprob, rtol + + @pytest.mark.parametrize("a, u, dim", genz_params) + def case_genz_productpeak(self, dim, a, u): + a_vec = np.repeat(a, dim) if a is not None else None + u_vec = np.repeat(u, dim) if u is not None else None + quadprob = genz_productpeak(dim=dim, a=a_vec, u=u_vec) + rtol = 1e-3 + return quadprob, rtol + + +dim_params = [1, 2, 10] + + +class OtherIntegrandsUniformCases: + @pytest.mark.parametrize("dim", dim_params) + def case_bratley1992(self, dim): + quadprob = bratley1992(dim=dim) + rtol = 1e-3 + return quadprob, rtol + + @pytest.mark.parametrize("dim", dim_params) + def case_roos_arnold(self, dim): + quadprob = roos_arnold(dim=dim) + rtol = 1e-3 + return quadprob, rtol + + @pytest.mark.parametrize("dim", dim_params) + def case_gfunction(self, dim): + quadprob = gfunction(dim=dim) + rtol = 1e-3 + return quadprob, rtol + + @pytest.mark.parametrize("dim", dim_params) + def case_morokoff_caflisch_1(self, dim): + quadprob = morokoff_caflisch_1(dim=dim) + rtol = 1e-3 + return quadprob, rtol + + @pytest.mark.parametrize("dim", dim_params) + def case_morokoff_caflisch_2(self, dim): + quadprob = morokoff_caflisch_2(dim=dim) + rtol = 1e-3 + return quadprob, rtol