diff --git a/docs/source/api/problems/zoo.quad.rst b/docs/source/api/problems/zoo.quad.rst
new file mode 100644
index 000000000..8b263a767
--- /dev/null
+++ b/docs/source/api/problems/zoo.quad.rst
@@ -0,0 +1,6 @@
+Quadrature
+----------
+
+.. automodapi:: probnum.problems.zoo.quad
+ :no-heading:
+ :headings: "*"
diff --git a/docs/source/api/problems/zoo.rst b/docs/source/api/problems/zoo.rst
index 296dacc14..33047d42b 100644
--- a/docs/source/api/problems/zoo.rst
+++ b/docs/source/api/problems/zoo.rst
@@ -19,3 +19,8 @@ probnum.problems.zoo
:hidden:
zoo.linalg
+
+.. toctree::
+ :hidden:
+
+ zoo.quad
diff --git a/src/probnum/problems/_problems.py b/src/probnum/problems/_problems.py
index 54a2f2398..eb86ef208 100644
--- a/src/probnum/problems/_problems.py
+++ b/src/probnum/problems/_problems.py
@@ -8,7 +8,7 @@
import numpy as np
import scipy.sparse
-from probnum import linops, randvars
+from probnum import linops, quad, randvars
from probnum.typing import FloatLike
@@ -225,51 +225,46 @@ class QuadratureProblem:
Compute the integral
- .. math::
- \int_\Omega f(x) \, \text{d} \mu(x)
+ .. math::
+ \int_\Omega f(x) \, \text{d} \mu(x)
- for a function :math:`f: \Omega \rightarrow \mathbb{R}`.
- For the time being, :math:`\mu` is the Lebesgue measure.
- Solved by quadrature rules in :mod:`probnum.quad`.
+ for a function :math:`f: \Omega \rightarrow \mathbb{R}` w.r.t. the measure
+ :math:`\mu`.
Parameters
----------
- integrand
- Function to be integrated.
- lower_bd
- A number or a vector representing the lower bounds of the integrals.
- upper_bd
- A number or a vector representing the upper bounds of the integrals.
- output_dim
- Output dimension of the integrand.
+ fun
+ Function to be integrated. It needs to accept a shape=(n_eval, input_dim)
+ ``np.ndarray`` and return a shape=(n_eval,) ``np.ndarray``.
+ measure
+ The integration measure.
solution
- Closed form, analytic solution to the problem.
+ Analytic value of the integral or precise numerical solution.
Used for testing and benchmarking.
Examples
--------
>>> import numpy as np
- >>> def integrand(x):
- ... return np.linalg.norm(x)**2
- >>> lower_bd = 0.41
- >>> upper_bd = 4.32
- >>> qp1d = QuadratureProblem(integrand, lower_bd=lower_bd, upper_bd=upper_bd)
- >>> np.round(qp1d.integrand(0.2), 2)
- 0.04
- >>> qp1d.lower_bd
- 0.41
+ >>> from probnum.quad.integration_measures import LebesgueMeasure
+ >>>
+ >>> def fun(x):
+ ... return np.linalg.norm(x, axis=1)**2
>>>
- >>> lower_bd = [0., 0.]
- >>> upper_bd = [1., 1.]
- >>> qp2d = QuadratureProblem(integrand, lower_bd=lower_bd, upper_bd=upper_bd)
- >>> qp2d.upper_bd
- [1.0, 1.0]
+ >>> measure1d = LebesgueMeasure(domain=(0, 1), input_dim=1)
+ >>> qp1d = QuadratureProblem(fun, measure=measure1d, solution=1/3)
+ >>> np.round(qp1d.fun(np.array([[0.2]]))[0], 2)
+ 0.04
+ >>> measure2d = LebesgueMeasure(domain=(0, 1), input_dim=2)
+ >>> qp2d = QuadratureProblem(fun, measure=measure2d, solution=None)
+ >>> np.round(qp2d.fun(np.array([[0.2, 0.2]]))[0], 2)
+ 0.08
"""
- integrand: Callable[[np.ndarray], Union[float, np.ndarray]]
- lower_bd: Union[FloatLike, np.ndarray]
- upper_bd: Union[FloatLike, np.ndarray]
- output_dim: Optional[int] = 1
+ fun: Callable[[np.ndarray], np.ndarray]
+ measure: quad.integration_measures.IntegrationMeasure
# For testing and benchmarking
- solution: Optional[Union[float, np.ndarray, randvars.RandomVariable]] = None
+ solution: Optional[Union[float, np.ndarray, randvars.RandomVariable]]
+
+ def __post_init__(self):
+ self.input_dim = self.measure.input_dim
diff --git a/src/probnum/problems/zoo/quad/__init__.py b/src/probnum/problems/zoo/quad/__init__.py
new file mode 100644
index 000000000..9d0ffd79f
--- /dev/null
+++ b/src/probnum/problems/zoo/quad/__init__.py
@@ -0,0 +1,6 @@
+"""Test problems for numerical integration/ quadrature."""
+
+from ._emukit_problems import circulargaussian2d, hennig1d, hennig2d, sombrero2d
+
+# Public classes and functions. Order is reflected in documentation.
+__all__ = ["circulargaussian2d", "hennig1d", "hennig2d", "sombrero2d"]
diff --git a/src/probnum/problems/zoo/quad/_emukit_problems.py b/src/probnum/problems/zoo/quad/_emukit_problems.py
new file mode 100644
index 000000000..4d32666b4
--- /dev/null
+++ b/src/probnum/problems/zoo/quad/_emukit_problems.py
@@ -0,0 +1,197 @@
+"""Toy integrands from Emukit."""
+
+# The integrands are re-implementations from scratch based on the Emukit docs.
+# There is no guarantee that they are identical to the Emukit implementations.
+
+from typing import Optional
+
+import numpy as np
+
+from probnum.problems import QuadratureProblem
+from probnum.quad.integration_measures import LebesgueMeasure
+from probnum.typing import FloatLike
+
+
+def hennig1d() -> QuadratureProblem:
+ r"""The univariate hennig function integrated wrt the Lebesgue measure. [1]_
+
+ The integrand is
+
+ .. math::
+ f(x) = e^{-x^2 -\sin^2(3x)}
+
+ on the domain :math:`\Omega=[-3, 3]`.
+
+ Returns
+ -------
+ quad_problem:
+ The quadrature problem.
+
+ References
+ ----------
+ .. [1] Emukit docs on `hennig1d `__.
+
+ """ # pylint: disable=line-too-long
+
+ def fun(x):
+ return np.exp(-x[:, 0] ** 2 - np.sin(3.0 * x[:, 0]) ** 2)
+
+ measure = LebesgueMeasure(input_dim=1, domain=(-3, 3))
+ return QuadratureProblem(fun=fun, measure=measure, solution=1.1433287777179366)
+
+
+def hennig2d(c: Optional[np.ndarray] = None) -> QuadratureProblem:
+ r"""The two-dimensional hennig function integrated wrt the Lebesgue measure. [1]_
+
+ The integrand is
+
+ .. math::
+ f(x) = e^{-x^{\intercal}c x -\sin(3\|x\|^2)}
+
+ on the domain :math:`\Omega=[-3, 3]^2`. Above, :math:`c` is the ``c`` parameter.
+
+ Parameters
+ ----------
+ c
+ A positive definite matrix of shape (2, 2). Defaults to [[1, .5], [.5, 1]].
+
+ Returns
+ -------
+ quad_problem:
+ The quadrature problem.
+
+ References
+ ----------
+ .. [1] Emukit docs on `hennig2d `__ .
+
+ """ # pylint: disable=line-too-long
+
+ solution = None
+ if c is None:
+ c = np.array([[1, 0.5], [0.5, 1]])
+ solution = 3.525721820076955
+
+ if c.shape != (2, 2):
+ raise ValueError(f"'c' must be a (2, 2) array. Found shape is {c.shape}.")
+
+ eigvals = np.linalg.eigvals(c)
+ if np.any(eigvals <= 0):
+ raise ValueError("'c' must be positive definite.")
+
+ def fun(x):
+ return np.exp(-np.sum((x @ c) * x, axis=1) - np.sin(3 * np.sum(x**2, axis=1)))
+
+ measure = LebesgueMeasure(input_dim=2, domain=(-3, 3))
+ return QuadratureProblem(fun=fun, measure=measure, solution=solution)
+
+
+def sombrero2d(w: Optional[FloatLike] = None) -> QuadratureProblem:
+ r"""The two-dimensional sombrero function integrated wrt the Lebesgue
+ measure. [1]_
+
+ The integrand is
+
+ .. math::
+ f(x) = \frac{\operatorname{sin}(\pi r w)}{\pi r w}
+
+ on the domain :math:`\Omega=[-3, 3]^2`. Above, :math:`w` is the ``w``
+ parameter and :math:`r=\|x\|` is the norm of the input vector :math:`x`.
+
+ Parameters
+ ----------
+ w
+ The positive frequency parameter. Defaults to 1.0.
+
+ Returns
+ -------
+ quad_problem:
+ The quadrature problem.
+
+ References
+ ----------
+ .. [1] Emukit docs on `sombrero2d `__ .
+
+ """ # pylint: disable=line-too-long
+
+ solution = None
+ if w is None:
+ w = 1.0
+ solution = 0.85225026427372
+
+ if w <= 0:
+ raise ValueError(f"The 'w' parameter must be positive ({w}).")
+
+ w = float(w)
+
+ def fun(x):
+ r_scaled = (np.pi * w) * np.sqrt((x * x).sum(axis=1))
+ f = np.sin(r_scaled) / r_scaled
+ f[np.isnan(f)] = 1.0
+ return f
+
+ measure = LebesgueMeasure(input_dim=2, domain=(-3, 3))
+ return QuadratureProblem(fun=fun, measure=measure, solution=solution)
+
+
+def circulargaussian2d(
+ m: Optional[FloatLike] = None, v: Optional[FloatLike] = None
+) -> QuadratureProblem:
+ r"""The two-dimensional circular Gaussian integrated wrt the Lebesgue
+ measure. [1]_
+
+ The integrand is
+
+ .. math::
+ f(x) = (2\pi v)^{-\frac{1}{2}} r^2 e^{-\frac{(r - m)^2}{2 v}}
+
+ on the domain :math:`\Omega=[-3, 3]^2`. Above, :math:`v` is the ``v``
+ parameter, :math:`m` is the ``m`` parameter and :math:`r = \|x\|` is the
+ norm of the input vector :math:`x`.
+
+ Parameters
+ ----------
+ m
+ The non-negative mean of the circular Gaussian in units of radius.
+ Defaults to 0.0.
+ v
+ The positive variance of the circular Gaussian. Defaults to 1.0.
+
+ Returns
+ -------
+ quad_problem:
+ The quadrature problem.
+
+ References
+ ----------
+ .. [1] Emukit docs on `circulargaussian2d `__ .
+
+ """ # pylint: disable=line-too-long
+
+ _v = 1.0
+ _m = 0.0
+
+ solution = None
+ if m is None and v is None:
+ v, m = _v, _m
+ solution = 4.853275495632483
+
+ if m is None:
+ m = _m
+ if v is None:
+ v = _v
+
+ if m < 0:
+ raise ValueError(f"'m' ({m}) must be non-negative.")
+
+ if v <= 0:
+ raise ValueError(f"'v' ({v}) must be positive.")
+
+ m, v = float(m), float(v)
+
+ def fun(x):
+ r = np.linalg.norm(x, axis=1)
+ rel_square_diff = (r - m) ** 2 / (2.0 * v)
+ return r**2 * np.exp(-rel_square_diff) / np.sqrt(2.0 * np.pi * v)
+
+ measure = LebesgueMeasure(input_dim=2, domain=(-3, 3))
+ return QuadratureProblem(fun=fun, measure=measure, solution=solution)
diff --git a/tests/test_problems/test_zoo/test_quad/test_problems.py b/tests/test_problems/test_zoo/test_quad/test_problems.py
new file mode 100644
index 000000000..59fee2edc
--- /dev/null
+++ b/tests/test_problems/test_zoo/test_quad/test_problems.py
@@ -0,0 +1,94 @@
+"""Tests for quadrature problems."""
+
+import numpy as np
+import pytest
+
+from probnum.problems.zoo.quad import circulargaussian2d, hennig1d, hennig2d, sombrero2d
+from probnum.quad.integration_measures import LebesgueMeasure
+
+
+@pytest.fixture(params=[hennig1d, hennig2d, sombrero2d, circulargaussian2d])
+def problem(request):
+ """quad problems with default values"""
+ return request.param()
+
+
+@pytest.fixture(
+ params=[
+ (hennig2d, dict(c=np.array([[2, 0.1], [0.1, 1]]))),
+ (circulargaussian2d, dict(m=1.0, v=1.5)),
+ (sombrero2d, dict(w=2.0)),
+ ],
+ ids=["hennig2d", "circulargaussian2d", "sombrero2d"],
+)
+def problem_parametrized(request):
+ """quad problems with custom parameters"""
+ return request.param[0](**request.param[1])
+
+
+def test_problem_solution_type(problem):
+ assert isinstance(problem.solution, float)
+
+
+def test_problem_parametrized_solution_type(problem_parametrized):
+ # problems that do not have an analytic solution currently do not provide any.
+ assert problem_parametrized.solution is None
+
+
+@pytest.mark.parametrize("num_dat", [1, 5])
+def test_problem_fun_shapes(problem, num_dat):
+ input_dim = problem.input_dim
+ res = problem.fun(np.ones([num_dat, input_dim]))
+ assert res.shape == (num_dat,)
+
+
+def test_problem_correct_solution_value(problem):
+ s = problem.solution
+ m = problem.measure
+
+ # 1e3 samples for MC estimator.
+ # Test might not work for high dimensions (>5). Testing integrals in high dimensions
+ # is intrinsically hard.
+ if s is not None and isinstance(s, float):
+ x = m.sample(int(1e3), rng=np.random.default_rng(0))
+ f = problem.fun(x)
+ s_test = f.mean()
+ c_test = f.std()
+
+ # scale MC estimator with domain volume for unnormalized Lebesgue measure.
+ if isinstance(m, LebesgueMeasure) and not m.normalized:
+ volume = np.prod(m.domain[1] - m.domain[0])
+ s_test *= volume
+ c_test *= volume
+
+ # Check if integral lies in a 95% confidence interval of the MC estimator.
+ assert s_test - 2 * c_test < s < s_test + 2 * c_test
+
+
+def test_problem_sombrero2d_raises():
+
+ # frequency is negative
+ with pytest.raises(ValueError):
+ sombrero2d(w=-2.0)
+
+
+def test_problem_hennig2d_raises():
+
+ # wrong c shape
+ with pytest.raises(ValueError):
+ hennig2d(c=np.eye(3))
+
+ # c not pos def
+ with pytest.raises(ValueError):
+ hennig2d(c=-np.eye(2))
+
+
+def test_problem_circulargaussian2d_raises():
+
+ # mean m negative
+ with pytest.raises(ValueError):
+ circulargaussian2d(m=-2.0)
+
+ # mean v non-positive
+ with pytest.raises(ValueError):
+ circulargaussian2d(v=0.0)