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

@@ -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",
+]

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,

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)
 

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,
+    )