diff --git a/src/probnum/quad/__init__.py b/src/probnum/quad/__init__.py index f1c024efd..c404378d6 100644 --- a/src/probnum/quad/__init__.py +++ b/src/probnum/quad/__init__.py @@ -6,7 +6,7 @@ choosing points to evaluate the integrand based on said model. """ -from probnum.quad.solvers.policies import Policy, RandomPolicy +from probnum.quad.solvers.policies import Policy, RandomPolicy, VanDerCorputPolicy from probnum.quad.solvers.stopping_criteria import ( BQStoppingCriterion, ImmediateStop, @@ -41,6 +41,7 @@ "IntegralVarianceTolerance", "MaxNevals", "RandomPolicy", + "VanDerCorputPolicy", "RelativeMeanChange", ] diff --git a/src/probnum/quad/_bayesquad.py b/src/probnum/quad/_bayesquad.py index 5827cc678..09c05f327 100644 --- a/src/probnum/quad/_bayesquad.py +++ b/src/probnum/quad/_bayesquad.py @@ -30,11 +30,13 @@ def bayesquad( domain: Optional[DomainLike] = None, measure: Optional[IntegrationMeasure] = None, policy: Optional[str] = "bmc", + scale_estimation: Optional[str] = "mle", max_evals: Optional[IntLike] = None, var_tol: Optional[FloatLike] = None, rel_tol: Optional[FloatLike] = None, - batch_size: Optional[IntLike] = 1, + batch_size: IntLike = 1, rng: Optional[np.random.Generator] = np.random.default_rng(), + jitter: FloatLike = 1.0e-8, ) -> Tuple[Normal, BQIterInfo]: r"""Infer the solution of the uni- or multivariate integral :math:`\int_\Omega f(x) d \mu(x)` @@ -78,6 +80,18 @@ def bayesquad( Bayesian Monte Carlo [2]_ ``bmc`` ========================== ======= + ========================== ======= + van Der Corput points ``vdc`` + ========================== ======= + + scale_estimation + Estimation method to use to compute the scale parameter. Defaults to 'mle'. + Options are + + ============================== ======= + Maximum likelihood estimation ``mle`` + ============================== ======= + max_evals Maximum number of function evaluations. var_tol @@ -85,10 +99,13 @@ def bayesquad( rel_tol Tolerance on consecutive updates of the integral mean. batch_size - Number of new observations at each update. + Number of new observations at each update. Defaults to 1. rng Random number generator. Used by Bayesian Monte Carlo other random sampling policies. Optional. Default is `np.random.default_rng()`. + jitter + Non-negative jitter to numerically stabilise kernel matrix inversion. + Defaults to 1e-8. Returns ------- @@ -147,11 +164,13 @@ def bayesquad( measure=measure, domain=domain, policy=policy, + scale_estimation=scale_estimation, max_evals=max_evals, var_tol=var_tol, rel_tol=rel_tol, batch_size=batch_size, rng=rng, + jitter=jitter, ) # Integrate @@ -166,6 +185,8 @@ def bayesquad_from_data( kernel: Optional[Kernel] = None, domain: Optional[DomainLike] = None, measure: Optional[IntegrationMeasure] = None, + scale_estimation: Optional[str] = "mle", + jitter: FloatLike = 1.0e-8, ) -> Tuple[Normal, BQIterInfo]: r"""Infer the value of an integral from a given set of nodes and function evaluations. @@ -184,6 +205,11 @@ def bayesquad_from_data( ``np.ndarray``. Obsolete if ``measure`` is given. measure The integration measure. Defaults to the Lebesgue measure. + scale_estimation + Estimation method to use to compute the scale parameter. Defaults to 'mle'. + jitter + Non-negative jitter to numerically stabilise kernel matrix inversion. + Defaults to 1e-8. Returns ------- @@ -231,6 +257,8 @@ def bayesquad_from_data( measure=measure, domain=domain, policy=None, + scale_estimation=scale_estimation, + jitter=jitter, ) # Integrate diff --git a/src/probnum/quad/_utils.py b/src/probnum/quad/_utils.py index b914cfac3..cd3754d20 100644 --- a/src/probnum/quad/_utils.py +++ b/src/probnum/quad/_utils.py @@ -28,7 +28,7 @@ def as_domain( Parameters ---------- domain - The integration domain as supllied. + The integration domain as supplied. input_dim The input dimensionality as supplied. @@ -44,7 +44,7 @@ def as_domain( ValueError If ``input_dim`` is not positive. If domain has too many or too little elements. - If the bounds of the domain have differening sizes. + If the bounds of the domain have differing sizes. If ``input_dim`` is incompatible with domain bounds. If bounds have wrong shape. If integration domain is empty. diff --git a/src/probnum/quad/solvers/bayesian_quadrature.py b/src/probnum/quad/solvers/bayesian_quadrature.py index aaef0fbe7..36f13de6f 100644 --- a/src/probnum/quad/solvers/bayesian_quadrature.py +++ b/src/probnum/quad/solvers/bayesian_quadrature.py @@ -5,7 +5,7 @@ import numpy as np -from probnum.quad.solvers.policies import Policy, RandomPolicy +from probnum.quad.solvers.policies import Policy, RandomPolicy, VanDerCorputPolicy from probnum.quad.solvers.stopping_criteria import ( BQStoppingCriterion, ImmediateStop, @@ -23,6 +23,8 @@ from .belief_updates import BQBeliefUpdate, BQStandardBeliefUpdate from .bq_state import BQIterInfo, BQState +# pylint: disable=too-many-branches, too-complex + class BayesianQuadrature: r"""The Bayesian quadrature method. @@ -75,11 +77,13 @@ def from_problem( measure: Optional[IntegrationMeasure] = None, domain: Optional[DomainLike] = None, policy: Optional[str] = "bmc", + scale_estimation: Optional[str] = "mle", max_evals: Optional[IntLike] = None, var_tol: Optional[FloatLike] = None, rel_tol: Optional[FloatLike] = None, batch_size: IntLike = 1, rng: np.random.Generator = None, + jitter: FloatLike = 1.0e-8, ) -> "BayesianQuadrature": r"""Creates an instance of this class from a problem description. @@ -97,6 +101,8 @@ def from_problem( policy The policy choosing nodes at which to evaluate the integrand. Choose ``None`` if you want to integrate from a fixed dataset. + scale_estimation + Estimation method to use to compute the scale parameter. Defaults to 'mle'. max_evals Maximum number of evaluations as stopping criterion. var_tol @@ -104,9 +110,12 @@ def from_problem( rel_tol Relative tolerance as stopping criterion. batch_size - Batch size used in node acquisition. + Batch size used in node acquisition. Defaults to 1. rng The random number generator. + jitter + Non-negative jitter to numerically stabilise kernel matrix inversion. + Defaults to 1e-8. Returns ------- @@ -150,14 +159,15 @@ def from_problem( ) raise ValueError(errormsg) policy = RandomPolicy(measure.sample, batch_size=batch_size, rng=rng) - + elif policy == "vdc": + policy = VanDerCorputPolicy(measure=measure, batch_size=batch_size) else: - raise NotImplementedError( - "Policies other than random sampling are not available at the moment." - ) + raise NotImplementedError(f"The given policy ({policy}) is unknown.") # Select the belief updater - belief_update = BQStandardBeliefUpdate() + belief_update = BQStandardBeliefUpdate( + jitter=jitter, scale_estimation=scale_estimation + ) # Select stopping criterion: If multiple stopping criteria are given, BQ stops # once any criterion is fulfilled (logical `or`). diff --git a/src/probnum/quad/solvers/belief_updates/_belief_update.py b/src/probnum/quad/solvers/belief_updates/_belief_update.py index c19ef1e2b..84a1d06f6 100644 --- a/src/probnum/quad/solvers/belief_updates/_belief_update.py +++ b/src/probnum/quad/solvers/belief_updates/_belief_update.py @@ -1,29 +1,35 @@ """Belief updates for Bayesian quadrature.""" import abc -from typing import Tuple +from typing import Optional, Tuple import numpy as np from scipy.linalg import cho_factor, cho_solve +from probnum.quad.kernel_embeddings import KernelEmbedding from probnum.quad.solvers.bq_state import BQState +from probnum.randprocs.kernels import Kernel from probnum.randvars import Normal +from probnum.typing import FloatLike # pylint: disable=too-few-public-methods, too-many-locals class BQBeliefUpdate(abc.ABC): - """Abstract class for the inference scheme.""" - - @abc.abstractmethod - def __call__(self, *args, **kwargs): - pass + """Abstract class for the inference scheme. + Parameters + ---------- + jitter + Non-negative jitter to numerically stabilise kernel matrix inversion. + """ -class BQStandardBeliefUpdate(BQBeliefUpdate): - """Updates integral belief and state using standard Bayesian quadrature based on - standard Gaussian process inference.""" + def __init__(self, jitter: FloatLike) -> None: + if jitter < 0: + raise ValueError(f"Jitter ({jitter}) must be non-negative.") + self.jitter = float(jitter) + @abc.abstractmethod def __call__( self, bq_state: BQState, @@ -51,20 +57,77 @@ def __call__( updated_state : Updated version of ``bq_state`` that contains all updated quantities. """ + raise NotImplementedError - # Update nodes and function evaluations - old_nodes = bq_state.nodes + def _compute_gram_cho_factor(self, gram: np.ndarray) -> np.ndarray: + """Compute the Cholesky decomposition of a positive-definite Gram matrix for use + in scipy.linalg.cho_solve + + .. warning:: + Uses scipy.linalg.cho_factor. The returned matrix is only to be used in + scipy.linalg.cho_solve. + Parameters + ---------- + gram : + symmetric pos. def. kernel Gram matrix :math:`K`, shape (nevals, nevals) + + Returns + ------- + gram_cho_factor : + The upper triangular Cholesky decomposition of the Gram matrix. Other + parts of the matrix contain random data. + """ + return cho_factor(gram + self.jitter * np.eye(gram.shape[0])) + + # pylint: disable=no-self-use + def _gram_cho_solve(self, gram_cho_factor: np.ndarray, z: np.ndarray) -> np.ndarray: + """Wrapper for scipy.linalg.cho_solve. Meant to be used for linear systems of + the gram matrix. Requires the solution of scipy.linalg.cho_factor as input.""" + return cho_solve(gram_cho_factor, z) + + +class BQStandardBeliefUpdate(BQBeliefUpdate): + """Updates integral belief and state using standard Bayesian quadrature based on + standard Gaussian process inference. + + Parameters + ---------- + jitter + Non-negative jitter to numerically stabilise kernel matrix inversion. + scale_estimation + Estimation method to use to compute the scale parameter. + """ + + def __init__(self, jitter: FloatLike, scale_estimation: Optional[str]) -> None: + super().__init__(jitter=jitter) + self.scale_estimation = scale_estimation + + def __call__( + self, + bq_state: BQState, + new_nodes: np.ndarray, + new_fun_evals: np.ndarray, + *args, + **kwargs, + ) -> Tuple[Normal, BQState]: + + # Update nodes and function evaluations nodes = np.concatenate((bq_state.nodes, new_nodes), axis=0) fun_evals = np.append(bq_state.fun_evals, new_fun_evals) - # kernel quantities - if old_nodes.size == 0: - gram = bq_state.kernel.matrix(new_nodes) - kernel_means = bq_state.kernel_embedding.kernel_mean(new_nodes) + # Estimate intrinsic kernel parameters + new_kernel, kernel_was_updated = self._estimate_kernel(bq_state.kernel) + new_kernel_embedding = KernelEmbedding(new_kernel, bq_state.measure) + + # Update gram matrix and kernel mean vector. Recompute everything from + # scratch if the kernel was updated or if these are the first nodes. + if kernel_was_updated or bq_state.nodes.size == 0: + gram = new_kernel.matrix(nodes) + kernel_means = new_kernel_embedding.kernel_mean(nodes) else: - gram_new_new = bq_state.kernel.matrix(new_nodes) - gram_old_new = bq_state.kernel.matrix(new_nodes, old_nodes) + gram_new_new = new_kernel.matrix(new_nodes) + gram_old_new = new_kernel.matrix(new_nodes, bq_state.nodes) gram = np.hstack( ( np.vstack((bq_state.gram, gram_old_new)), @@ -74,47 +137,58 @@ def __call__( kernel_means = np.concatenate( ( bq_state.kernel_means, - bq_state.kernel_embedding.kernel_mean(new_nodes), + new_kernel_embedding.kernel_mean(new_nodes), ) ) - initial_integral_variance = bq_state.kernel_embedding.kernel_variance() - weights = self._solve_gram(gram, kernel_means) + # Cholesky factorisation of the Gram matrix + gram_cho_factor = self._compute_gram_cho_factor(gram) - # integral mean and variance + # Estimate scaling parameter + new_scale_sq = self._estimate_scale(fun_evals, gram_cho_factor, bq_state) + + # Integral mean and variance + weights = self._gram_cho_solve(gram_cho_factor, kernel_means) integral_mean = weights @ fun_evals - integral_variance = initial_integral_variance - weights @ kernel_means + initial_integral_variance = new_kernel_embedding.kernel_variance() + integral_variance = new_scale_sq * ( + initial_integral_variance - weights @ kernel_means + ) - updated_belief = Normal(integral_mean, integral_variance) - updated_state = BQState.from_new_data( + new_belief = Normal(integral_mean, integral_variance) + new_state = BQState.from_new_data( + kernel=new_kernel, + scale_sq=new_scale_sq, nodes=nodes, fun_evals=fun_evals, - integral_belief=updated_belief, + integral_belief=new_belief, prev_state=bq_state, gram=gram, kernel_means=kernel_means, ) - return updated_belief, updated_state - - @staticmethod - def _solve_gram(gram: np.ndarray, rhs: np.ndarray) -> np.ndarray: - """Solve the linear system of the form. - - .. math:: Kx=b, - - Parameters - ---------- - gram : - symmetric pos. def. kernel Gram matrix :math:`K`, shape (nevals, nevals) - rhs : - right-hand-side :math:`b`, matrix or vector, shape (nevals, ...) - - Returns - ------- - x : - The solution to the linear system :math:`K x = b` - """ - jitter = 1.0e-6 - chol_gram = cho_factor(gram + jitter * np.eye(gram.shape[0])) - return cho_solve(chol_gram, rhs) + return new_belief, new_state + + # pylint: disable=no-self-use + def _estimate_kernel(self, kernel: Kernel) -> Tuple[Kernel, bool]: + """Estimate the intrinsic kernel parameters. That is, all parameters except the + scale.""" + new_kernel = kernel + kernel_was_updated = False + return new_kernel, kernel_was_updated + + def _estimate_scale( + self, fun_evals: np.ndarray, gram_cho_factor: np.ndarray, bq_state: BQState + ) -> FloatLike: + """Estimate the scale parameter.""" + if self.scale_estimation is None: + new_scale_sq = bq_state.scale_sq + elif self.scale_estimation == "mle": + new_scale_sq = ( + fun_evals + @ self._gram_cho_solve(gram_cho_factor, fun_evals) + / fun_evals.shape[0] + ) + else: + raise ValueError(f"Scale estimation ({self.scale_estimation}) is unknown.") + return new_scale_sq diff --git a/src/probnum/quad/solvers/bq_state.py b/src/probnum/quad/solvers/bq_state.py index 4b6e380ee..854e90178 100644 --- a/src/probnum/quad/solvers/bq_state.py +++ b/src/probnum/quad/solvers/bq_state.py @@ -9,6 +9,7 @@ from probnum.quad.kernel_embeddings import KernelEmbedding from probnum.randprocs.kernels import Kernel from probnum.randvars import Normal +from probnum.typing import FloatLike # pylint: disable=too-few-public-methods,too-many-instance-attributes @@ -22,12 +23,20 @@ class BQState: The integration measure. kernel The kernel used for BQ. + scale_sq + Square of the kernel scaling parameter. integral_belief Normal distribution over the value of the integral. + previous_integral_beliefs + Integral beliefs on computed on previous iterations. nodes All locations at which function evaluations are available. fun_evals Function evaluations at nodes. + gram + The kernel Gram matrix. + kernel_means + All kernel mean evaluations at ``nodes``. See Also -------- @@ -38,6 +47,7 @@ def __init__( self, measure: IntegrationMeasure, kernel: Kernel, + scale_sq: FloatLike = 1.0, integral_belief: Optional[Normal] = None, previous_integral_beliefs: Tuple[Normal] = (), nodes: Optional[np.ndarray] = None, @@ -48,6 +58,7 @@ def __init__( self.measure = measure self.kernel = kernel self.kernel_embedding = KernelEmbedding(kernel, measure) + self.scale_sq = scale_sq self.integral_belief = integral_belief self.previous_integral_beliefs = previous_integral_beliefs self.input_dim = measure.input_dim @@ -65,6 +76,8 @@ def __init__( @classmethod def from_new_data( cls, + kernel: Kernel, + scale_sq: FloatLike, nodes: np.ndarray, fun_evals: np.ndarray, integral_belief: Normal, @@ -76,6 +89,10 @@ def from_new_data( Parameters ---------- + kernel + The kernel used for BQ. + scale_sq + Square of the kernel scaling parameter. nodes All locations at which function evaluations are available. fun_evals @@ -96,7 +113,8 @@ def from_new_data( """ return cls( measure=prev_state.measure, - kernel=prev_state.kernel, + kernel=kernel, + scale_sq=scale_sq, integral_belief=integral_belief, previous_integral_beliefs=prev_state.previous_integral_beliefs + (prev_state.integral_belief,), diff --git a/src/probnum/quad/solvers/policies/__init__.py b/src/probnum/quad/solvers/policies/__init__.py index 8ad97ab97..61ff04af1 100644 --- a/src/probnum/quad/solvers/policies/__init__.py +++ b/src/probnum/quad/solvers/policies/__init__.py @@ -2,3 +2,4 @@ from ._policy import Policy from ._random_policy import RandomPolicy +from ._van_der_corput_policy import VanDerCorputPolicy diff --git a/src/probnum/quad/solvers/policies/_van_der_corput_policy.py b/src/probnum/quad/solvers/policies/_van_der_corput_policy.py new file mode 100644 index 000000000..9b72c007e --- /dev/null +++ b/src/probnum/quad/solvers/policies/_van_der_corput_policy.py @@ -0,0 +1,99 @@ +"""Van der Corput points for integration on 1D intervals.""" + +from typing import Optional + +import numpy as np + +from probnum.quad._integration_measures import IntegrationMeasure +from probnum.quad.solvers.bq_state import BQState + +from ._policy import Policy + + +class VanDerCorputPolicy(Policy): + r"""Pick nodes from the van der Corput sequence. + + The van der Corput sequence [1]_ is + + .. math:: 0.5, 0.25, 0.75, 0.125, 0.625, \ldots + + If the integration domain is not [0, 1], the van der Corput sequence is linearly + mapped to the domain. The domain must be finite. + + Parameters + ---------- + measure + The integration measure with finite domain. + batch_size + Size of batch of nodes when calling the policy once. + + References + -------- + .. [1] https://en.wikipedia.org/wiki/Van_der_Corput_sequence + """ + + def __init__(self, measure: IntegrationMeasure, batch_size: int) -> None: + super().__init__(batch_size=batch_size) + + if int(measure.input_dim) > 1: + raise ValueError("Policy 'vdc' works only when the input dimension is one.") + + domain_a = measure.domain[0] + domain_b = measure.domain[1] + if np.Inf in [abs(domain_a), abs(domain_b)]: + raise ValueError("Policy 'vdc' works only for bounded domains.") + + self.domain_a = domain_a + self.domain_b = domain_b + + def __call__(self, bq_state: BQState) -> np.ndarray: + n_nodes = bq_state.nodes.shape[0] + vdc_seq = VanDerCorputPolicy.van_der_corput_sequence( + n_nodes + 1, n_nodes + 1 + self.batch_size + ) + transformed_vdc_seq = vdc_seq * (self.domain_b - self.domain_a) + self.domain_a + return transformed_vdc_seq.reshape((self.batch_size, 1)) + + @staticmethod + def van_der_corput_sequence( + n_start: int, n_end: Optional[int] = None + ) -> np.ndarray: + r"""Returns elements ``n_start``, ``n_start + 1``, ..., ``n_end - 1`` in the van + der Corput sequence. + + .. math:: 0.5, 0.25, 0.75, 0.125, 0.625, \ldots + + If no ``n_end`` is given, only a single element in the + sequence is returned. + + Parameters + ---------- + n_start + First element of the van der Corput to be included (inclusive). + n_end + Last element of the van der Corput to be included (exclusive). If not given, + only the ``n_start`` element is returned. + + Returns + ------- + vdc_seq + Array containing elements from ``n_start`` to ``n_end - 1`` of the van der + Corput sequence. + """ + + # pylint: disable=invalid-name + if n_end is None: + n_end = n_start + 1 + vdc_seq = np.zeros((n_end - n_start,)) + ind = 0 + for m in range(n_start, n_end): + q = 0.0 + base_inv = 0.5 + n = m + while n != 0: + q = q + (n % 2) * base_inv + n = n // 2 + base_inv = base_inv / 2.0 + vdc_seq[ind] = q + ind += 1 + return vdc_seq diff --git a/src/probnum/quad/solvers/stopping_criteria/_rel_mean_change.py b/src/probnum/quad/solvers/stopping_criteria/_rel_mean_change.py index 339b5d2c9..ecb097445 100644 --- a/src/probnum/quad/solvers/stopping_criteria/_rel_mean_change.py +++ b/src/probnum/quad/solvers/stopping_criteria/_rel_mean_change.py @@ -31,6 +31,9 @@ def __init__(self, rel_tol: FloatLike): def __call__(self, bq_state: BQState, info: BQIterInfo) -> bool: integral_belief = bq_state.integral_belief + if not bq_state.previous_integral_beliefs: + # On the first iteration there is no previous integral value to use + return False return ( np.abs( (integral_belief.mean - bq_state.previous_integral_beliefs[-1].mean) diff --git a/tests/test_quad/conftest.py b/tests/test_quad/conftest.py index 96954e1e9..48cc4c508 100644 --- a/tests/test_quad/conftest.py +++ b/tests/test_quad/conftest.py @@ -108,7 +108,7 @@ def fixture_measure_params( @pytest.fixture(name="measure") def fixture_measure(measure_params) -> measures.IntegrationMeasure: - """Kernel / covariance function.""" + """Measure.""" name = measure_params.pop("name") if name == "gauss": diff --git a/tests/test_quad/test_bayesian_quadrature.py b/tests/test_quad/test_bayesian_quadrature.py index 092050444..60b7373f1 100644 --- a/tests/test_quad/test_bayesian_quadrature.py +++ b/tests/test_quad/test_bayesian_quadrature.py @@ -9,6 +9,7 @@ ImmediateStop, LebesgueMeasure, RandomPolicy, + VanDerCorputPolicy, ) from probnum.randprocs.kernels import ExpQuad @@ -53,6 +54,17 @@ def test_bq_from_problem_wrong_inputs(input_dim): BayesianQuadrature.from_problem(input_dim=input_dim) +@pytest.mark.parametrize( + "policy, policy_type", [("bmc", RandomPolicy), ("vdc", VanDerCorputPolicy)] +) +def test_bq_from_problem_policy_assignment(policy, policy_type): + """Test if correct policy is assigned from string identifier.""" + bq = BayesianQuadrature.from_problem( + input_dim=1, domain=(0, 1), policy=policy, rng=np.random.default_rng() + ) + assert isinstance(bq.policy, policy_type) + + def test_bq_from_problem_defaults(bq_no_policy, bq): # default policy and stopping criterion diff --git a/tests/test_quad/test_bayesquad/test_bq.py b/tests/test_quad/test_bayesquad/test_bq.py index 043812071..1ce644c73 100644 --- a/tests/test_quad/test_bayesquad/test_bq.py +++ b/tests/test_quad/test_bayesquad/test_bq.py @@ -2,8 +2,9 @@ import numpy as np import pytest -from scipy.integrate import quad +from scipy.integrate import quad as scipyquad +import probnum.quad from probnum.quad import bayesquad, bayesquad_from_data from probnum.quad.kernel_embeddings import KernelEmbedding from probnum.randvars import Normal @@ -26,9 +27,32 @@ def test_type_1d(f1d, kernel, measure, input_dim): @pytest.mark.parametrize("input_dim", [1]) -def test_integral_values_1d(f1d, kernel, measure, input_dim): - """Test numerically that BQ computes 1D integrals correctly.""" - +@pytest.mark.parametrize( + "domain", + [ + ( + 0, + 1, + ), + (-0.5, 1), + (-3.5, -2.9), + ], +) +@pytest.mark.parametrize("var_tol", [None, 1e-7]) +@pytest.mark.parametrize("rel_tol", [None, 1e-7]) +@pytest.mark.parametrize("scale_estimation", [None, "mle"]) +@pytest.mark.parametrize("jitter", [1e-6, 1e-7]) +def test_integral_values_1d( + f1d, kernel, domain, input_dim, scale_estimation, var_tol, rel_tol, jitter +): + """Test numerically that BQ computes 1D integrals correctly for a number of + different parameters. + + The test currently uses van der Corput policy and therefore works only for finite + domains. + """ + + measure = probnum.quad.LebesgueMeasure(input_dim=input_dim, domain=domain) # numerical integral # pylint: disable=invalid-name def integrand(x): @@ -36,19 +60,27 @@ def integrand(x): # pylint: disable=invalid-name bq_integral, _ = bayesquad( - fun=f1d, input_dim=input_dim, kernel=kernel, measure=measure, max_evals=250 + fun=f1d, + input_dim=input_dim, + kernel=kernel, + domain=domain, + policy="vdc", + scale_estimation=scale_estimation, + max_evals=250, + var_tol=var_tol, + rel_tol=rel_tol, + jitter=jitter, ) domain = measure.domain - if domain is None: - domain = (-np.infty, np.infty) - num_integral, _ = quad(integrand, domain[0], domain[1]) + num_integral, _ = scipyquad(integrand, domain[0], domain[1]) np.testing.assert_almost_equal(bq_integral.mean, num_integral, decimal=2) @pytest.mark.parametrize("input_dim", [2, 3, 4]) @pytest.mark.parametrize("measure_name", ["gauss"]) @pytest.mark.parametrize("cov_diagonal", [True]) -def test_integral_values_x2_gaussian(kernel, measure, input_dim): +@pytest.mark.parametrize("scale_estimation", [None, "mle"]) +def test_integral_values_x2_gaussian(kernel, measure, input_dim, scale_estimation): """Test numerical integration of x**2 in higher dimensions.""" # pylint: disable=invalid-name c = np.linspace(0.1, 2.2, input_dim) @@ -60,14 +92,22 @@ def test_integral_values_x2_gaussian(kernel, measure, input_dim): ) fun_evals = fun(nodes) bq_integral, _ = bayesquad_from_data( - nodes=nodes, fun_evals=fun_evals, kernel=kernel, measure=measure + nodes=nodes, + fun_evals=fun_evals, + kernel=kernel, + measure=measure, + scale_estimation=scale_estimation, ) np.testing.assert_almost_equal(bq_integral.mean, true_integral, decimal=2) @pytest.mark.parametrize("input_dim", [2, 3, 4]) @pytest.mark.parametrize("measure_name", ["lebesgue"]) -def test_integral_values_sin_lebesgue(kernel, measure, input_dim): +@pytest.mark.parametrize("scale_estimation", [None, "mle"]) +@pytest.mark.parametrize("jitter", [1e-6, 0.5e-5]) +def test_integral_values_sin_lebesgue( + kernel, measure, input_dim, scale_estimation, jitter +): """Test numerical integration of products of sinusoids.""" # pylint: disable=invalid-name c = np.linspace(0.1, 0.5, input_dim) @@ -85,7 +125,12 @@ def test_integral_values_sin_lebesgue(kernel, measure, input_dim): ) fun_evals = fun(nodes) bq_integral, _ = bayesquad_from_data( - nodes=nodes, fun_evals=fun_evals, kernel=kernel, measure=measure + nodes=nodes, + fun_evals=fun_evals, + kernel=kernel, + measure=measure, + scale_estimation=scale_estimation, + jitter=jitter, ) np.testing.assert_almost_equal(bq_integral.mean, true_integral, decimal=2) @@ -146,3 +191,22 @@ def test_domain_ignored_if_lebesgue(input_dim, measure): nodes=nodes, fun_evals=fun_evals, domain=domain, measure=measure ) assert isinstance(bq_integral, Normal) + + +def test_zero_function_gives_zero_variance_with_mle(): + """Test that BQ variance is zero for zero function when MLE is used to set the + scale parameter.""" + input_dim = 1 + domain = (0, 1) + fun = lambda x: np.zeros(x.shape[0]) + nodes = np.linspace(0, 1, 3).reshape((3, 1)) + fun_evals = fun(nodes) + + bq_integral1, _ = bayesquad( + fun=fun, input_dim=input_dim, domain=domain, scale_estimation="mle" + ) + bq_integral2, _ = bayesquad_from_data( + nodes=nodes, fun_evals=fun_evals, domain=domain, scale_estimation="mle" + ) + assert bq_integral1.var == 0.0 + assert bq_integral2.var == 0.0 diff --git a/tests/test_quad/test_belief_update.py b/tests/test_quad/test_belief_update.py new file mode 100644 index 000000000..7a7041589 --- /dev/null +++ b/tests/test_quad/test_belief_update.py @@ -0,0 +1,12 @@ +"""Test cases for the BQ belief updater.""" + +import pytest + +from probnum.quad.solvers.belief_updates import BQStandardBeliefUpdate + + +def test_belief_update_raises(): + # negative jitter is not allowed + wrong_jitter = -1.0 + with pytest.raises(ValueError): + BQStandardBeliefUpdate(jitter=wrong_jitter, scale_estimation="mle") diff --git a/tests/test_quad/test_bq_state.py b/tests/test_quad/test_bq_state.py index a299a8baa..9294f8844 100644 --- a/tests/test_quad/test_bq_state.py +++ b/tests/test_quad/test_bq_state.py @@ -106,6 +106,7 @@ def test_state_defaults_types(state, request): assert isinstance(s.gram, np.ndarray) assert isinstance(s.kernel_means, np.ndarray) assert isinstance(s.previous_integral_beliefs, tuple) + assert isinstance(s.scale_sq, float) assert s.integral_belief is None @@ -150,9 +151,13 @@ def test_state_from_new_data(state, request): integral = Normal(0, 1) gram = np.eye(new_nevals) kernel_means = np.ones(new_nevals) + kernel = ExpQuad(input_shape=(old_state.input_dim,)) + scale_sq = 1.7 # previously no data given s = BQState.from_new_data( + kernel=kernel, + scale_sq=scale_sq, nodes=x, fun_evals=y, integral_belief=integral, @@ -181,3 +186,4 @@ def test_state_from_new_data(state, request): # values assert s.input_dim == s.measure.input_dim + assert s.scale_sq == scale_sq diff --git a/tests/test_quad/test_policy.py b/tests/test_quad/test_policy.py new file mode 100644 index 000000000..c4e9b0181 --- /dev/null +++ b/tests/test_quad/test_policy.py @@ -0,0 +1,53 @@ +"""Basic tests for BQ policies.""" + +import numpy as np +import pytest + +from probnum.quad import GaussianMeasure, LebesgueMeasure, VanDerCorputPolicy + + +def test_van_der_corput_multi_d_error(): + """Check that van der Corput policy fails in dimensions higher than one.""" + wrong_dimension = 2 + measure = GaussianMeasure(input_dim=wrong_dimension, mean=0.0, cov=1.0) + with pytest.raises(ValueError): + VanDerCorputPolicy(measure, batch_size=1) + + +@pytest.mark.parametrize("domain", [(-np.Inf, 0), (1, np.Inf), (-np.Inf, np.Inf)]) +def test_van_der_corput_infinite_error(domain): + """Check that van der Corput policy fails on infinite domains.""" + measure = LebesgueMeasure(input_dim=1, domain=domain) + with pytest.raises(ValueError): + VanDerCorputPolicy(measure, batch_size=1) + + +@pytest.mark.parametrize("n", [4, 8, 16, 32, 64, 128, 256]) +def test_van_der_corput_full(n): + """Test that the full van der Corput sequence is being computed correctly.""" + # Full sequence + vdc_seq = VanDerCorputPolicy.van_der_corput_sequence(1, n) + expected_seq = np.linspace(1.0 / n, 1.0 - 1.0 / n, n - 1) + np.testing.assert_array_equal(np.sort(vdc_seq), expected_seq) + + +def test_van_der_corput_partial(): + """Test that some partial van der Corput sequences are computed correctly.""" + (n_start, n_end) = (3, 8) + vdc_seq = VanDerCorputPolicy.van_der_corput_sequence(n_start, n_end) + expected_seq = np.array([0.75, 0.125, 0.625, 0.375, 0.875]) + np.testing.assert_array_equal(vdc_seq, expected_seq) + + (n_start, n_end) = (4, 11) + vdc_seq = VanDerCorputPolicy.van_der_corput_sequence(n_start, n_end) + expected_seq = np.array([0.125, 0.625, 0.375, 0.875, 0.0625, 0.5625, 0.3125]) + np.testing.assert_array_equal(vdc_seq, expected_seq) + + +def test_van_der_corput_start_value_only(): + """When no end value is given, test if sequence returns the correct value.""" + (n_start, n_end) = (1, 8) + vdc_seq = VanDerCorputPolicy.van_der_corput_sequence(n_start, n_end) + vdc_seq_single_value = VanDerCorputPolicy.van_der_corput_sequence(n_end - 1) + assert vdc_seq_single_value.shape == (1,) + assert vdc_seq[-1] == vdc_seq_single_value[0]