Source code for probnum.quad.solvers.bayesian_quadrature

"""Probabilistic numerical methods for solving integrals."""

from typing import Callable, Optional, Tuple
import warnings

import numpy as np

from probnum.quad.solvers.policies import Policy, RandomPolicy
from probnum.quad.solvers.stopping_criteria import (
from probnum.randprocs.kernels import ExpQuad, Kernel
from probnum.randvars import Normal
from probnum.typing import FloatLike, IntLike

from .._integration_measures import IntegrationMeasure, LebesgueMeasure
from .._quad_typing import DomainLike
from ..kernel_embeddings import KernelEmbedding
from .belief_updates import BQBeliefUpdate, BQStandardBeliefUpdate
from .bq_state import BQIterInfo, BQState

class BayesianQuadrature:
    r"""The Bayesian quadrature method.

    Bayesian quadrature solves integrals of the form

    .. math:: F = \int_\Omega f(x) d \mu(x).

        The kernel used for the GP model.
        The integration measure.
        The policy choosing nodes at which to evaluate the integrand.
        The inference method.
        The criterion that determines convergence.

    See Also
    bayesquad : Computes the integral using an acquisition policy.
    bayesquad_from_data : Computes the integral :math:`F` using a given dataset of
                          nodes and function evaluations.


    def __init__(
        kernel: Kernel,
        measure: IntegrationMeasure,
        policy: Optional[Policy],
        belief_update: BQBeliefUpdate,
        stopping_criterion: BQStoppingCriterion,
    ) -> None:
        self.kernel = kernel
        self.measure = measure
        self.policy = policy
        self.belief_update = belief_update
        self.stopping_criterion = stopping_criterion

[docs] @classmethod def from_problem( cls, input_dim: int, kernel: Optional[Kernel] = None, measure: Optional[IntegrationMeasure] = None, domain: Optional[DomainLike] = None, policy: Optional[str] = "bmc", max_evals: Optional[IntLike] = None, var_tol: Optional[FloatLike] = None, rel_tol: Optional[FloatLike] = None, batch_size: IntLike = 1, rng: np.random.Generator = None, ) -> "BayesianQuadrature": r"""Creates an instance of this class from a problem description. Parameters ---------- input_dim The input dimension. kernel The kernel used for the GP model. Defaults to the ``ExpQuad`` kernel. measure The integration measure. Defaults to the Lebesgue measure on the ``domain``. domain The integration bounds. Obsolete if ``measure`` is given. policy The policy choosing nodes at which to evaluate the integrand. Choose ``None`` if you want to integrate from a fixed dataset. max_evals Maximum number of evaluations as stopping criterion. var_tol Variance tolerance as stopping criterion. rel_tol Relative tolerance as stopping criterion. batch_size Batch size used in node acquisition. rng The random number generator. Returns ------- BayesianQuadrature An instance of this class. Raises ------ ValueError If neither a ``domain`` nor a ``measure`` are given. ValueError If Bayesian Monte Carlo ('bmc') is selected as ``policy`` and no random number generator (``rng``) is given. NotImplementedError If an unknown ``policy`` is given. """ # Set up integration measure if domain is None and measure is None: raise ValueError( "You need to either specify an integration domain or a measure." ) if measure is None: measure = LebesgueMeasure(domain=domain, input_dim=input_dim) # Select the kernel if kernel is None: kernel = ExpQuad(input_shape=(input_dim,)) # Select policy if policy is None: # If policy is None, this implies that the integration problem is defined # through a fixed set of nodes and function evaluations which will not # require an acquisition loop. The error handling is done in ``integrate``. pass elif policy == "bmc": if rng is None: errormsg = ( "Policy 'bmc' relies on random sampling, " "thus requires a random number generator ('rng')." ) raise ValueError(errormsg) policy = RandomPolicy(measure.sample, batch_size=batch_size, rng=rng) else: raise NotImplementedError( "Policies other than random sampling are not available at the moment." ) # Select the belief updater belief_update = BQStandardBeliefUpdate() # Select stopping criterion: If multiple stopping criteria are given, BQ stops # once any criterion is fulfilled (logical `or`). def _stopcrit_or(sc1, sc2): if sc1 is None: return sc2 return sc1 | sc2 _stopping_criterion = None if max_evals is not None: _stopping_criterion = _stopcrit_or( _stopping_criterion, MaxNevals(max_evals) ) if var_tol is not None: _stopping_criterion = _stopcrit_or( _stopping_criterion, IntegralVarianceTolerance(var_tol) ) if rel_tol is not None: _stopping_criterion = _stopcrit_or( _stopping_criterion, RelativeMeanChange(rel_tol) ) # If no stopping criteria are given, use some default values. if _stopping_criterion is None: _stopping_criterion = IntegralVarianceTolerance(var_tol=1e-6) | MaxNevals( max_nevals=input_dim * 25 # 25 is an arbitrary value ) # If no policy is given, then the iteration must terminate immediately. if policy is None: _stopping_criterion = ImmediateStop() return cls( kernel=kernel, measure=measure, policy=policy, belief_update=belief_update, stopping_criterion=_stopping_criterion, )
[docs] def bq_iterator( self, bq_state: BQState, info: Optional[BQIterInfo], fun: Optional[Callable], ) -> Tuple[Normal, BQState, BQIterInfo]: """Generator that implements the iteration of the BQ method. This function exposes the state of the BQ method one step at a time while running the loop. Parameters ---------- bq_state State of the Bayesian quadrature methods. Contains the information about the problem and the BQ belief. info The state of the iteration. 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``. Yields ------ new_integral_belief : Updated belief about the integral. new_bq_state : The updated state of the Bayesian quadrature belief. new_info : The updated state of the iteration. """ # Setup iteration info if info is None: info = BQIterInfo.from_bq_state(bq_state) while True: _has_converged = self.stopping_criterion(bq_state, info) info = BQIterInfo.from_stopping_decision(info, has_converged=_has_converged) yield bq_state.integral_belief, bq_state, info # Have we already converged? if _has_converged: break # Select new nodes via policy new_nodes = self.policy(bq_state=bq_state) # Evaluate the integrand at new nodes new_fun_evals = fun(new_nodes) # Update the belief about the integrand _, bq_state = self.belief_update( bq_state=bq_state, new_nodes=new_nodes, new_fun_evals=new_fun_evals, ) # Update the state of the iteration info = BQIterInfo.from_iteration(info=info, dnevals=self.policy.batch_size)
[docs] def integrate( self, fun: Optional[Callable], nodes: Optional[np.ndarray], fun_evals: Optional[np.ndarray], ) -> Tuple[Normal, BQState, BQIterInfo]: """Integrates the function ``fun``. ``fun`` may be analytically given, or numerically in terms of ``fun_evals`` at fixed nodes. This function calls the generator ``bq_iterator`` until the first stopping criterion is met. It immediately stops after processing the initial ``nodes`` if ``policy`` is not available. Parameters ---------- 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``. nodes *shape=(n_eval, input_dim)* -- Optional nodes at which function evaluations are available as ``fun_evals`` from start. fun_evals *shape=(n_eval,)* -- Optional function evaluations at ``nodes`` available from the start. Returns ------- integral_belief : Posterior belief about the integral. bq_state : Final state of the Bayesian quadrature method. Raises ------ ValueError If neither the integrand function (``fun``) nor integrand evaluations (``fun_evals``) are given. ValueError If ``nodes`` are not given and no policy is present. ValueError If dimension of ``nodes`` or ``fun_evals`` is incorrect, or if their shapes do not match. """ # no policy given: Integrate on fixed dataset. if self.policy is None: # nodes must be provided if no policy is given. if nodes is None: raise ValueError("No policy available: Please provide nodes.") # Use fun_evals and disregard fun if both are given if fun is not None and fun_evals is not None: warnings.warn( "No policy available: 'fun_eval' are used instead of 'fun'." ) fun = None # override stopping condition as no policy is given. self.stopping_criterion = ImmediateStop() # Check if integrand function is provided if fun is None and fun_evals is None: raise ValueError( "Please provide an integrand function 'fun' or function values " "'fun_evals'." ) # Setup initial design if nodes is not None and fun_evals is None: fun_evals = fun(nodes) # Check if shapes of nodes and function evaluations match if fun_evals is not None and fun_evals.ndim != 1: raise ValueError( f"fun_evals must be one-dimensional " f"({fun_evals.ndim})." ) if nodes is not None and nodes.ndim != 2: raise ValueError(f"nodes must be two-dimensional ({nodes.ndim}).") if nodes is not None and fun_evals is not None: if nodes.shape[0] != fun_evals.shape[0]: raise ValueError( f"nodes ({nodes.shape[0]}) and fun_evals " f"({fun_evals.shape[0]}) need to contain the same number " f"of evaluations." ) # Setup BQ state: This encodes a zero-mean prior. bq_state = BQState( measure=self.measure, kernel=self.kernel, integral_belief=Normal( 0.0, KernelEmbedding(self.kernel, self.measure).kernel_variance() ), ) if nodes is not None: _, bq_state = self.belief_update( bq_state=bq_state, new_nodes=nodes, new_fun_evals=fun_evals, ) info = None for (_, bq_state, info) in self.bq_iterator(bq_state, info, fun): pass return bq_state.integral_belief, bq_state, info