diff --git a/docs/source/api/problems/zoo.quad.rst b/docs/source/api/problems/zoo.quad.rst index 8b263a767..7c655e800 100644 --- a/docs/source/api/problems/zoo.quad.rst +++ b/docs/source/api/problems/zoo.quad.rst @@ -1,5 +1,4 @@ Quadrature ----------- .. automodapi:: probnum.problems.zoo.quad :no-heading: diff --git a/src/probnum/problems/zoo/quad/__init__.py b/src/probnum/problems/zoo/quad/__init__.py index 9d0ffd79f..987103016 100644 --- a/src/probnum/problems/zoo/quad/__init__.py +++ b/src/probnum/problems/zoo/quad/__init__.py @@ -2,5 +2,8 @@ from ._emukit_problems import circulargaussian2d, hennig1d, hennig2d, sombrero2d +from ._quadproblems_gaussian import * +from ._quadproblems_uniform import * + # Public classes and functions. Order is reflected in documentation. __all__ = ["circulargaussian2d", "hennig1d", "hennig2d", "sombrero2d"] diff --git a/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py b/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py new file mode 100644 index 000000000..31ab179cf --- /dev/null +++ b/src/probnum/problems/zoo/quad/_quadproblems_gaussian.py @@ -0,0 +1,197 @@ +"""Test problems for integration against a Gaussian measure.""" + +from typing import Callable + +import numpy as np +from scipy import special +from scipy.stats import norm + +from probnum.problems import QuadratureProblem +from probnum.typing import FloatLike + +__all__ = [ + "uniform_to_gaussian_quadprob", + "sum_polynomials", +] + + +def uniform_to_gaussian_quadprob( + quadprob: QuadratureProblem, mean: FloatLike = 0.0, var: FloatLike = 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 + suitable for integration against the Lebesgue measure on :math:`[0,1]^d`. + + The multivariate Gaussian is of the form :math:`\mathcal{N}(mean \cdot (1, \dotsc, + 1)^\top, var^2 \cdot I_d)`. + + Using the change of variable formula, we have that + + .. math:: \int_{[0,1]^d} f(x) dx = \int_{\mathbb{R}^d} h(x) \phi(x) dx + + where :math:`h(x)=f(\Phi((x-mean)/var))`, :math:`\phi(x)` is the Gaussian + probability density function and :math:`\Phi(x)` an elementwise application of the + Gaussian cummulative distribution function. See [1]_. + + Parameters + ---------- + 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. + + Returns + ------- + problem + A new Quadrature Problem instance with a transformed integrand taking inputs in + :math:`\mathbb{R}^d`. + + Raises + ------ + ValueError + If the original quadrature problem is over a domain other than [0, 1]^d or if it + does not have a scalar solution. + + Example + ------- + uniform_to_gaussian_quadprob(genz_continuous(1)) + + References + ---------- + .. [1] Si, S., Oates, C. J., Duncan, A. B., Carin, L. & Briol. F-X. (2021). Scalable + control variates for Monte Carlo methods via stochastic optimization. Proceedings + of the 14th Monte Carlo and Quasi-Monte Carlo Methods (MCQMC) conference 2020. + arXiv:2006.07487. + """ + + dim = len(quadprob.lower_bd) + + # Check that the original quadrature problem was defined on [0,1]^d + if np.any(quadprob.lower_bd != 0.0): + raise ValueError("quadprob is not an integration problem over [0,1]^d") + if np.any(quadprob.upper_bd != 1.0): + raise ValueError("quadprob is not an integration problem over [0,1]^d") + + # Check that the original quadrature problem has a scalar valued solution + 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( + func: 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 func(norm.cdf((x - mean) / var)) + + return newfunc + + return QuadratureProblem( + integrand=uniform_to_gaussian_integrand( + func=quadprob.integrand, mean=mean, var=var + ), + lower_bd=np.broadcast_to(-np.Inf, dim), + upper_bd=np.broadcast_to(np.Inf, dim), + output_dim=None, + solution=quadprob.solution, + ) + + +def sum_polynomials( + dim: int, a: np.ndarray = None, b: np.ndarray = None, var: FloatLike = 1.0 +) -> QuadratureProblem: + r"""Quadrature problem with an integrand taking the form of a sum of polynomials + over :math:`\mathbb{R}^d`. + + .. math:: f(x) = \sum_{j=0}^p \prod_{i=1}^dim a_{ji} x_i^{b_ji} + + The integrand is integrated against a multivariate normal :math:`\mathcal{N}(0,var * + I_d)`. See [1]_. + + Parameters + ---------- + dim + Dimension d of the integration domain + a + 2d-array of size (p+1)xd giving coefficients of the polynomials. + b + 2d-array of size (p+1)xd giving orders of the polynomials. All entries + should be natural numbers. + var + diagonal elements of the covariance function. + + Returns + ------- + f + array of size (n,1) giving integrand evaluations at points in 'x'. + + Raises + ------ + ValueError + If the given parameters have the wrong shape or contain invalid values. + + References + ---------- + .. [1] Si, S., Oates, C. J., Duncan, A. B., Carin, L. & Briol. F-X. (2021). Scalable + control variates for Monte Carlo methods via stochastic optimization. Proceedings + of the 14th Monte Carlo and Quasi-Monte Carlo Methods (MCQMC) conference 2020. + arXiv:2006.07487. + """ + + # Specify default values of parameters a and u + if a is None: + a = np.broadcast_to(1.0, (1, dim)) + if b is None: + b = np.broadcast_to(1, (1, dim)) + + # Check that the parameters have valid values and shape + if a.shape[1] != dim: + raise ValueError( + f"Invalid shape {a.shape} for parameter `a`. Expected {dim} columns." + ) + + if b.shape[1] != dim: + raise ValueError( + f"Invalid shape {b.shape} for parameter `b`. Expected {dim} columns." + ) + + if np.any(b < 0): + raise ValueError("The parameters `b` must be non-negative.") + + def integrand(x: np.ndarray) -> np.ndarray: + n = x.shape[0] + + # Compute function values + f = np.sum( + np.prod( + a[np.newaxis, :, :] * (x[:, np.newaxis, :] ** b[np.newaxis, :, :]), + axis=2, + ), + axis=1, + ) + + # Return function values as a 2d array with one column. + return f.reshape((n, 1)) + + # Return function values as a 2d array with one column. + delta = (np.remainder(b, 2) - 1) ** 2 + doublefact = special.factorial2(b - 1) + solution = np.sum(np.prod(a * delta * (var**b) * doublefact, axis=1)) + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(-np.Inf, dim), + upper_bd=np.broadcast_to(np.Inf, dim), + output_dim=None, + solution=solution, + ) diff --git a/src/probnum/problems/zoo/quad/_quadproblems_uniform.py b/src/probnum/problems/zoo/quad/_quadproblems_uniform.py new file mode 100644 index 000000000..c4f577d3c --- /dev/null +++ b/src/probnum/problems/zoo/quad/_quadproblems_uniform.py @@ -0,0 +1,772 @@ +"""Test problems for integration against the Lebesgue measure.""" + +import itertools + +import numpy as np +from scipy.stats import norm + +from probnum.problems import QuadratureProblem + +__all__ = [ + "genz_continuous", + "genz_cornerpeak", + "genz_discontinuous", + "genz_gaussian", + "genz_oscillatory", + "genz_productpeak", + "bratley1992", + "roos_arnold", + "gfunction", + "morokoff_caflisch_1", + "morokoff_caflisch_2", +] + + +def genz_continuous( + dim: int, a: np.ndarray = None, u: np.ndarray = None +) -> QuadratureProblem: + r"""Genz 'continuous' test function on :math:`[0,1]^d`. + + .. math:: f(x) = \exp(- \sum_{i=1}^d a_i |x_i - u_i|). + + Parameters + ---------- + dim + Dimension of the domain + a + First set of parameters of shape (dim,) affecting the difficulty of the + integration problem. See [1]_. + u + Second set of parameters of shape (dim,) affecting the difficulty of the + integration problem. All entries should be in [0,1]. See [1]_. + + Returns + ------- + problem + The :class:`QuadratureProblem` corresponding to the Genz 'continuous' test + function with the given parameters. + + Raises + ------ + ValueError + If any of the parameters have invalid shapes or values. + + References + ---------- + .. [1] https://www.sfu.ca/~ssurjano/cont.html + """ + + # Specify default values of parameters a and u + if a is None: + a = np.broadcast_to(5.0, dim) + + if u is None: + u = np.broadcast_to(0.5, dim) + + # Check that the parameters have valid values and shape + if a.shape != (dim,): + raise ValueError( + f"Invalid shape {a.shape} for parameter `a`. Expected {(dim,)}." + ) + + if u.shape != (dim,): + raise ValueError( + f"Invalid shape {u.shape} for parameter `u`. Expected {(dim,)}." + ) + + if np.any(u < 0.0) or np.any(u > 1): + raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].") + + def integrand(x: np.ndarray) -> np.ndarray: + nonlocal a, u + n = x.shape[0] + + # Check that the input points have valid values and shape + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Reshape u into an (n,dim) array with identical rows + u_reshaped = np.repeat(u.reshape([1, dim]), n, axis=0) + + # Compute function values + f = np.exp(-np.sum(a * np.abs(x - u_reshaped), axis=1)) + + return f.reshape((n, 1)) + + solution = np.prod((2.0 - np.exp(-1.0 * a * u) - np.exp(a * (u - 1))) / a) + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) + + +def genz_cornerpeak( + dim: int, a: np.ndarray = None, u: np.ndarray = None +) -> QuadratureProblem: + r"""Genz 'corner peak' test function on :math:`[0,1]^d`. + + .. math:: f(x) = (1+\sum_{i=1}^d a_i x_i)^{-d+1} + + Parameters + ---------- + dim + Dimension of the domain + a + First set of parameters of shape (dim,) affecting the difficulty of the + integration problem. See [1]_. + u + Second set of parameters of shape (dim,) affecting the difficulty of the + integration problem. All entries should be in [0,1]. See [1]_. + + References + ---------- + .. [1] https://www.sfu.ca/~ssurjano/copeak.html + """ + + # Specify default values of parameters a and u + if a is None: + a = np.broadcast_to(5.0, dim) + + if u is None: + u = np.broadcast_to(0.5, dim) + + # Check that the parameters have valid values and shape + if a.shape != (dim,): + raise ValueError( + f"Invalid shape {a.shape} for parameter `a`. Expected {(dim,)}." + ) + + if u.shape != (dim,): + raise ValueError( + f"Invalid shape {u.shape} for parameter `u`. Expected {(dim,)}." + ) + + if np.any(u < 0.0) or np.any(u > 1): + raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].") + + def integrand(x: np.ndarray) -> np.ndarray: + nonlocal a, u + n = x.shape[0] + + # Check that the input points have valid values + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Compute function values + f = (1.0 + np.sum(a * x, axis=1)) ** (-dim - 1) + + return f.reshape((n, 1)) + + # Calculate closed-form solution of the integral + solution = 0.0 + for k in range(0, dim + 1): + subsets_k = list(itertools.combinations(range(dim), k)) + for subset in subsets_k: + a_subset = a[np.asarray(subset, dtype=int)] + solution = solution + ((-1.0) ** (k + dim)) * ( + 1.0 + np.sum(a) - np.sum(a_subset) + ) ** (-1) + solution = solution / (np.prod(a) * np.math.factorial(dim)) + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) + + +def genz_discontinuous( + dim: int, a: np.ndarray = None, u: np.ndarray = None +) -> QuadratureProblem: + r"""Genz 'discontinuous' test function on :math:`[0,1]^d`. + + .. math:: + f(x) = + \begin{cases} + 0 & \text{if any } x_i > u_i \\ + \exp(\sum_{i=1}^d a_i x_i) & \text{otherwise} + \end{cases} + + Parameters + ---------- + dim + Dimension of the domain + a + First set of parameters of shape (dim,) affecting the difficulty of the + integration problem. See [1]_. + u + Second set of parameters of shape (dim,) affecting the difficulty of the + integration problem. All entries should be in [0,1]. See [1]_. + + References + ---------- + .. [1] https://www.sfu.ca/~ssurjano/disc.html + """ + + # Specify default values of parameters a and u + if a is None: + a = np.broadcast_to(5.0, dim) + + if u is None: + u = np.broadcast_to(0.5, dim) + + # Check that the parameters have valid values and shape + if a.shape != (dim,): + raise ValueError( + f"Invalid shape {a.shape} for parameter `a`. Expected {(dim,)}." + ) + + if u.shape != (dim,): + raise ValueError( + f"Invalid shape {u.shape} for parameter `u`. Expected {(dim,)}." + ) + + if np.any(u < 0.0) or np.any(u > 1): + raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].") + + def integrand(x: np.ndarray) -> np.ndarray: + nonlocal a, u + n = x.shape[0] + + # Check that the input points have valid values + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Compute function values + f = np.exp(np.sum(a * x, axis=1)) + # Set function to zero whenever x_i > u_i for i =1,..,min(2,d) + f[np.any(x - u > 0, axis=1)] = 0 + + return f.reshape(n, 1) + + if dim == 1: + solution = (np.exp(a * u) - 1.0) / a + if dim > 1: + solution = np.prod((np.exp(a * np.minimum(u, 1.0)) - 1.0) / a) + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) + + +def genz_gaussian( + dim: int, a: np.ndarray = None, u: np.ndarray = None +) -> QuadratureProblem: + r"""Genz 'Gaussian' test function on :math:`[0,1]^d`. + + .. math:: f(x) = \exp(- \sum_{i=1}^d a_i^2 (x_i - u_i)^2). + + Parameters + ---------- + dim + Dimension of the domain + a + First set of parameters of shape (dim,) affecting the difficulty of the + integration problem. See [1]_. + u + Second set of parameters of shape (dim,) affecting the difficulty of the + integration problem. All entries should be in [0,1]. See [1]_. + + References + ---------- + .. [1] https://www.sfu.ca/~ssurjano/gaussian.html + """ + + # Specify default values of parameters a and u + if a is None: + a = np.broadcast_to(5.0, dim) + + if u is None: + u = np.broadcast_to(0.5, dim) + + # Check that the parameters have valid values and shape + if a.shape != (dim,): + raise ValueError( + f"Invalid shape {a.shape} for parameter `a`. Expected {(dim,)}." + ) + + if u.shape != (dim,): + raise ValueError( + f"Invalid shape {u.shape} for parameter `u`. Expected {(dim,)}." + ) + + if np.any(u < 0.0) or np.any(u > 1): + raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].") + + def integrand(x: np.ndarray) -> np.ndarray: + nonlocal a, u + n = x.shape[0] + + # Check that the input points have valid values + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Reshape u into an (n,dim) array with identical rows + u_reshaped = np.repeat(u.reshape([1, dim]), n, axis=0) + + # Compute function values + f = np.exp(-np.sum((a * (x - u_reshaped)) ** 2, axis=1)) + + return f.reshape((n, 1)) + + solution = np.pi ** (dim / 2) * np.prod( + (norm.cdf(np.sqrt(2) * a * (1.0 - u)) - norm.cdf(-np.sqrt(2) * a * u)) / a + ) + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) + + +def genz_oscillatory( + dim: int, a: np.ndarray = None, u: np.ndarray = None +) -> QuadratureProblem: + r"""Genz 'oscillatory' test function on :math:`[0,1]^d`. + + .. math:: f(x) = \cos( 2 \pi u_1 + \sum_{i=1}^d a_i x_i). + + + Parameters + ---------- + dim + Dimension of the domain + a + First set of parameters of shape (dim,) affecting the difficulty of the + integration problem. See [1]_. + u + Second set of parameters of shape (dim,) affecting the difficulty of the + integration problem. All entries should be in [0,1]. See [1]_. + + References + ---------- + .. [1] https://www.sfu.ca/~ssurjano/oscil.html + """ + + # pylint: disable=too-complex + + # Specify default values of parameters a and u + if a is None: + a = np.broadcast_to(5.0, dim) + + if u is None: + u = np.broadcast_to(0.5, dim) + + # Check that the parameters have valid values and shape + if a.shape != (dim,): + raise ValueError( + f"Invalid shape {a.shape} for parameter `a`. Expected {(dim,)}." + ) + + if u.shape != (dim,): + raise ValueError( + f"Invalid shape {u.shape} for parameter `u`. Expected {(dim,)}." + ) + + if np.any(u < 0.0) or np.any(u > 1): + raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].") + + def integrand(x: np.ndarray) -> np.ndarray: + nonlocal a, u + n = x.shape[0] + + # Check that the input points have valid values + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Compute function values + f = np.cos(2.0 * np.pi * u[0] + np.sum(a * x, axis=1)) + + return f.reshape((n, 1)) + + # Calculate closed-form solution to the integral + dim_modulo4 = np.remainder(dim, 4) + + def hfunc(x): + if dim_modulo4 == 1: + return np.sin(x) + if dim_modulo4 == 2: + return -np.cos(x) + if dim_modulo4 == 3: + return -np.sin(x) + + assert dim_modulo4 == 0 + + return np.cos(x) + + solution = 0.0 + for k in range(0, dim + 1): + subsets_k = list(itertools.combinations(range(dim), k)) + for subset in subsets_k: + a_subset = a[np.asarray(subset, dtype=int)] + solution = solution + ((-1.0) ** k) * hfunc( + (2.0 * np.pi * u[0]) + np.sum(a) - np.sum(a_subset) + ) + + solution = solution / np.prod(a) + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) + + +def genz_productpeak( + dim: int, a: np.ndarray = None, u: np.ndarray = None +) -> QuadratureProblem: + r"""Genz 'Product Peak' test function on :math:`[0,1]^d`. + + .. math:: f(x) = \prod_{i=1}^d ( a_i^{-2} + (x_i-u_i)^2)^{-1}. + + + Parameters + ---------- + dim + Dimension of the domain + a + First set of parameters of shape (dim,) affecting the difficulty of the + integration problem. See [1]_. + u + Second set of parameters of shape (dim,) affecting the difficulty of the + integration problem. All entries should be in [0,1]. See [1]_. + + References + ---------- + .. [1] https://www.sfu.ca/~ssurjano/prpeak.html + """ + + # Specify default values of parameters a and u + if a is None: + a = np.broadcast_to(5.0, dim) + + if u is None: + u = np.broadcast_to(0.5, dim) + + # Check that the parameters have valid values and shape + if a.shape != (dim,): + raise ValueError( + f"Invalid shape {a.shape} for parameter `a`. Expected {(dim,)}." + ) + + if u.shape != (dim,): + raise ValueError( + f"Invalid shape {u.shape} for parameter `u`. Expected {(dim,)}." + ) + + if np.any(u < 0.0) or np.any(u > 1.0): + raise ValueError("The parameters `u` must lie in the interval [0.0, 1.0].") + + def integrand(x: np.ndarray) -> np.ndarray: + nonlocal a, u + n = x.shape[0] + + # Check that the input points have valid values + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Reshape u into an (n,dim) array with identical rows + u_reshaped = np.repeat(u.reshape([1, dim]), n, axis=0) + + # Compute function values + f = np.prod(1.0 / (1.0 / a**2 + (x - u_reshaped) ** 2), axis=1) + + return f.reshape((n, 1)) + + solution = np.prod(a * (np.arctan(a * (1.0 - u)) - np.arctan(-1.0 * a * u))) + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) + + +def bratley1992(dim: int) -> QuadratureProblem: + r"""'Bratley 1992' test function on :math:`[0,1]^d`. + + .. math:: f(x) = \sum_{i=1}^d (-1)^i \prod_{j=1}^i x_j. + + See [1]_, [2]_. + + Parameters + ---------- + dim + Dimension of the domain + + References + ---------- + .. [1] Bratley, P., Fox, B. L., & Niederreiter, H. (1992). Implementation and tests + of low-discrepancy sequences. ACM Transactions on Modeling and Computer + Simulation (TOMACS), 2(3), 195-213. + .. [2] https://www.sfu.ca/~ssurjano/bratleyetal92.html + """ + + def integrand(x: np.ndarray) -> np.ndarray: + n = x.shape[0] + + # Check that the input points have valid values + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Compute function values + f = np.sum( + ((-1.0) ** np.arange(1, dim + 1)) * np.cumprod(x, axis=1), + axis=1, + ) + + return f.reshape((n, 1)) + + solution = -(1.0 / 3) * (1.0 - ((-0.5) ** dim)) + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) + + +def roos_arnold(dim: int) -> QuadratureProblem: + r"""'Roos & Arnold 1963' test function on :math:`[0,1]^d`. + + .. math:: f(x) = \prod_{i=1}^d |4 x_i - 2 |. + + See [1]_, [2]_. + + Parameters + ---------- + dim + Dimension of the domain + + References + ---------- + .. [1] Roos, P., & Arnold, L. (1963). Numerische experimente zur mehrdimensionalen + quadratur. Springer. + .. [2] https://www.sfu.ca/~ssurjano/roosarn63.html + """ + + def integrand(x: np.ndarray) -> np.ndarray: + n = x.shape[0] + + # Check that the input points have valid values + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Compute function values + f = np.prod(np.abs(4.0 * x - 2.0), axis=1) + + return f.reshape((n, 1)) + + solution = 1.0 + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) + + +def gfunction(dim: int) -> QuadratureProblem: + r"""'G-function' test function on :math:`[0,1]^d`. + + .. math:: f(x) = \prod_{i=1}^d \frac{|4 x_i - 2 |+a_i}{1+a_i} + + where :math:`a_i = \frac{i-2}{2}` for all :math:`i = 1, \dotsc, d`. See [1]_, [2]_. + + Parameters + ---------- + dim + Dimension of the domain + + References + ---------- + .. [1] Marrel, A., Iooss, B., Laurent, B., & Roustant, O. (2009). Calculations of + sobol indices for the gaussian process metamodel. Reliability Engineering & + System Safety, 94(3), 742-751. + .. [2] https://www.sfu.ca/~ssurjano/gfunc.html + """ + + def integrand(x: np.ndarray) -> np.ndarray: + n = x.shape[0] + + # Check that the input points have valid values + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Compute function values + a = np.atleast_2d(((np.arange(dim) + 1.0) - 2.0) / 2.0) + f = np.prod((np.abs(4.0 * x - 2.0) + a) / (1.0 + a), axis=1) + + return f.reshape((n, 1)) + + solution = 1.0 + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) + + +def morokoff_caflisch_1(dim: int) -> QuadratureProblem: + r"""'Morokoff & Caflisch 1995' test function number 1 on :math:`[0,1]^d`. + + .. math:: f(x) = (1+1/d)^d \prod_{i=1}^d x_i^{1/d} + + See [1]_, [2]_, [3]_. + + Parameters + ---------- + dim + Dimension of the domain + + References + ---------- + .. [1] Morokoff, W. J., & Caflisch, R. E. (1995). Quasi-monte carlo integration. + Journal of computational physics, 122(2), 218-230. + .. [2] Gerstner, T., & Griebel, M. (1998). Numerical integration using sparse grids. + Numerical algorithms, 18(3-4), 209-232. + .. [3] https://www.sfu.ca/~ssurjano/morcaf95a.html + """ + + def integrand(x: np.ndarray) -> np.ndarray: + n = x.shape[0] + + # Check that the input points have valid values + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Compute function values + f = ((1.0 + 1.0 / dim) ** (dim)) * np.prod(x ** (1.0 / dim), axis=1) + + return f.reshape((n, 1)) + + solution = 1.0 + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) + + +def morokoff_caflisch_2(dim: int) -> QuadratureProblem: + r"""'Morokoff & Caflisch 1995' test function number 2 on :math:`[0,1]^d`. + + .. math:: f(x) = \frac{1}{(d-0.5)^d} \prod_{i=1}^d (d-x_i) + + See [1]_, [2]_. + + Parameters + ---------- + dim + Dimension of the domain + + References + ---------- + .. [1] Morokoff, W. J., & Caflisch, R. E. (1995). Quasi-monte carlo integration. + Journal of computational physics, 122(2), 218-230. + .. [2] https://www.sfu.ca/~ssurjano/morcaf95b.html + """ + + def integrand(x: np.ndarray) -> np.ndarray: + n = x.shape[0] + + # Check that the input points have valid values + if x.shape != (n, dim): + raise ValueError( + f"Invalid shape {x.shape} for input points `x`. Expected (n, {dim})." + ) + + if np.any(x < 0.0) or np.any(x > 1): + raise ValueError("The input points `x` must lie in the box [0.0, 1.0]^d.") + + # Compute function values + f = (1.0 / ((dim - 0.5) ** dim)) * np.prod(dim - x, axis=1) + + return f.reshape((n, 1)) + + solution = 1.0 + + return QuadratureProblem( + integrand=integrand, + lower_bd=np.broadcast_to(0.0, dim), + upper_bd=np.broadcast_to(1.0, dim), + output_dim=None, + solution=solution, + ) diff --git a/tests/test_problems/test_zoo/test_quad/__init__.py b/tests/test_problems/test_zoo/test_quad/__init__.py new file mode 100644 index 000000000..e69de29bb 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 new file mode 100644 index 000000000..b062b43a7 --- /dev/null +++ b/tests/test_problems/test_zoo/test_quad/test_quadproblems_uniform.py @@ -0,0 +1,193 @@ +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, +) + + +@pytest.mark.parametrize( + "genz_problem", + [ + genz_continuous, + genz_cornerpeak, + genz_discontinuous, + genz_gaussian, + genz_oscillatory, + genz_productpeak, + ], +) +def test_genz_uniform_param_checks(genz_problem): + with pytest.raises(ValueError): + genz_problem(2, a=np.ones(shape=(1,))) + + with pytest.raises(ValueError): + genz_problem(3, u=np.ones(shape=(2, 1))) + + with pytest.raises(ValueError): + genz_problem(3, u=np.full((3,), 1.1)) + + with pytest.raises(ValueError): + genz_problem(3, u=np.full((3,), -0.1)) + + +@pytest.mark.parametrize( + "quad_problem_constructor", + [ + bratley1992, + genz_continuous, + genz_cornerpeak, + genz_discontinuous, + genz_gaussian, + genz_oscillatory, + genz_productpeak, + gfunction, + morokoff_caflisch_1, + morokoff_caflisch_2, + roos_arnold, + ], +) +def test_integrand_eval_checks(quad_problem_constructor): + quad_problem = quad_problem_constructor(2) + + with pytest.raises(ValueError): + quad_problem.integrand(np.zeros((4, 3))) + + with pytest.raises(ValueError): + quad_problem.integrand(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) + + # 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.integrand(x_unif)) / n, + quadprob_genz_continuous.solution, + rtol=1e-03, + ) + np.testing.assert_allclose( + np.sum(quadprob_genz_cornerpeak.integrand(x_unif)) / n, + quadprob_genz_cornerpeak.solution, + rtol=3e-03, + ) + np.testing.assert_allclose( + np.sum(quadprob_genz_discontinuous.integrand(x_unif)) / n, + quadprob_genz_discontinuous.solution, + rtol=1e-03, + ) + np.testing.assert_allclose( + np.sum(quadprob_genz_gaussian.integrand(x_unif)) / n, + quadprob_genz_gaussian.solution, + rtol=2e-03, + ) + np.testing.assert_allclose( + np.sum(quadprob_genz_oscillatory.integrand(x_unif)) / n, + quadprob_genz_oscillatory.solution, + rtol=9e-02, + ) + np.testing.assert_allclose( + np.sum(quadprob_genz_productpeak.integrand(x_unif)) / n, + quadprob_genz_productpeak.solution, + rtol=1e-03, + ) + + +@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. + + The former should be an approximation of the latter. + """ + + # 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) + + # Test that all Monte Carlo estimators approximate the true value of the integral + # (integration against uniform) + np.testing.assert_allclose( + np.sum(quadprob_bratley1992.integrand(x_unif)) / n, + quadprob_bratley1992.solution, + rtol=1e-03, + ) + np.testing.assert_allclose( + np.sum(quadprob_roos_arnold.integrand(x_unif)) / n, + quadprob_roos_arnold.solution, + rtol=1e-03, + ) + np.testing.assert_allclose( + np.sum(quadprob_gfunction.integrand(x_unif)) / n, + quadprob_gfunction.solution, + rtol=1e-03, + ) + np.testing.assert_allclose( + np.sum(quadprob_morokoff_caflisch_1.integrand(x_unif)) / n, + quadprob_morokoff_caflisch_1.solution, + rtol=1e-03, + ) + np.testing.assert_allclose( + np.sum(quadprob_morokoff_caflisch_2.integrand(x_unif)) / n, + quadprob_morokoff_caflisch_2.solution, + rtol=1e-03, + )