Source code for gpflow.quadrature.deprecated

# Copyright 2017-2020 The GPflow Contributors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import itertools
import warnings
from functools import wraps
from typing import Callable, Iterable, List, Optional, Tuple, Union

import numpy as np
import tensorflow as tf

from ..base import AnyNDArray, TensorType
from ..config import default_float
from ..experimental.check_shapes import check_shapes
from ..utilities import to_default_float
from .gauss_hermite import NDiagGHQuadrature


[docs]@check_shapes( "return[0]: [n_quad_points]", "return[1]: [n_quad_points]", ) def hermgauss(n: int) -> Tuple[AnyNDArray, AnyNDArray]: # Type-ignore is because for some versions mypy can't find np.polynomial.hermite x, w = np.polynomial.hermite.hermgauss(n) x, w = x.astype(default_float()), w.astype(default_float()) return x, w
[docs]@check_shapes( "return[0]: [n_quad_points, D]", "return[1]: [n_quad_points]", ) def mvhermgauss(H: int, D: int) -> Tuple[AnyNDArray, AnyNDArray]: """ Return the evaluation locations 'xn', and weights 'wn' for a multivariate Gauss-Hermite quadrature. The outputs can be used to approximate the following type of integral: int exp(-x)*f(x) dx ~ sum_i w[i,:]*f(x[i,:]) :param H: Number of Gauss-Hermite evaluation points. :param D: Number of input dimensions. Needs to be known at call-time. :return: eval_locations 'x' (H**DxD), weights 'w' (H**D) """ gh_x, gh_w = hermgauss(H) x: AnyNDArray = np.array(list(itertools.product(*(gh_x,) * D))) # H**DxD w = np.prod(np.array(list(itertools.product(*(gh_w,) * D))), 1) # H**D return x, w
[docs]@check_shapes( "means: [N, Din]", "covs: [N, Din, Din]", "return: [N, Dout...]", ) def mvnquad( func: Callable[[tf.Tensor], tf.Tensor], means: TensorType, covs: TensorType, H: int, Din: Optional[int] = None, Dout: Optional[Tuple[int, ...]] = None, ) -> tf.Tensor: """ Computes N Gaussian expectation integrals of a single function 'f' using Gauss-Hermite quadrature. :param f: integrand function. Takes one input of shape ?xD. :param H: Number of Gauss-Hermite evaluation points. :param Din: Number of input dimensions. Needs to be known at call-time. :param Dout: Number of output dimensions. Defaults to (). Dout is assumed to leave out the item index, i.e. f actually maps (?xD)->(?x*Dout). :return: quadratures """ # Figure out input shape information if Din is None: Din = means.shape[1] if Din is None: raise ValueError( "If `Din` is passed as `None`, `means` must have a known shape. " "Running mvnquad in `autoflow` without specifying `Din` and `Dout` " "is problematic. Consider using your own session." ) # pragma: no cover xn, wn = mvhermgauss(H, Din) N = means.shape[0] # transform points based on Gaussian parameters cholXcov = tf.linalg.cholesky(covs) # NxDxD Xt = tf.linalg.matmul( cholXcov, tf.tile(xn[None, :, :], (N, 1, 1)), transpose_b=True ) # NxDxH**D X = 2.0 ** 0.5 * Xt + tf.expand_dims(means, 2) # NxDxH**D Xr = tf.reshape(tf.transpose(X, [2, 0, 1]), (-1, Din)) # (H**D*N)xD # perform quadrature fevals = func(Xr) if Dout is None: Dout = tuple((d if type(d) is int else d.value) for d in fevals.shape[1:]) if any([d is None for d in Dout]): raise ValueError( "If `Dout` is passed as `None`, the output of `func` must have known " "shape. Running mvnquad in `autoflow` without specifying `Din` and `Dout` " "is problematic. Consider using your own session." ) # pragma: no cover fX = tf.reshape( fevals, ( H ** Din, N, ) + Dout, ) wr = np.reshape(wn * np.pi ** (-Din * 0.5), (-1,) + (1,) * (1 + len(Dout))) return tf.reduce_sum(fX * wr, 0)
[docs]@check_shapes( "Fmu: [broadcast Din, N...]", "Fvar: [broadcast Din, N...]", "Ys.values(): [N...]", "return: [broadcast Dout, N...]", ) def ndiagquad( funcs: Union[Callable[..., tf.Tensor], Iterable[Callable[..., tf.Tensor]]], H: int, Fmu: Union[TensorType, Tuple[TensorType, ...], List[TensorType]], Fvar: Union[TensorType, Tuple[TensorType, ...], List[TensorType]], logspace: bool = False, **Ys: TensorType, ) -> tf.Tensor: """ Computes N Gaussian expectation integrals of one or more functions using Gauss-Hermite quadrature. The Gaussians must be independent. The means and variances of the Gaussians are specified by Fmu and Fvar. The N-integrals are assumed to be taken wrt the last dimensions of Fmu, Fvar. `Fmu`, `Fvar`, `Ys` should all have same shape, with overall size `N`. :param funcs: the integrand(s): Callable or Iterable of Callables that operates elementwise :param H: number of Gauss-Hermite quadrature points :param Fmu: array/tensor or `Din`-tuple/list thereof :param Fvar: array/tensor or `Din`-tuple/list thereof :param logspace: if True, funcs are the log-integrands and this calculates the log-expectation of exp(funcs) :param Ys: arrays/tensors; deterministic arguments to be passed by name :return: shape is the same as that of the first Fmu """ warnings.warn( "Please use gpflow.quadrature.NDiagGHQuadrature instead " "(note the changed convention of how multi-dimensional quadrature is handled)", DeprecationWarning, ) n_gh = H if isinstance(Fmu, (tuple, list)): dim = len(Fmu) shape = tf.shape(Fmu[0]) Fmu = tf.stack(Fmu, axis=-1) Fvar = tf.stack(Fvar, axis=-1) else: dim = 1 shape = tf.shape(Fmu) Fmu = tf.reshape(Fmu, (-1, dim)) Fvar = tf.reshape(Fvar, (-1, dim)) Ys = {Yname: tf.reshape(Y, (-1, 1)) for Yname, Y in Ys.items()} def wrapper(old_fun: Callable[..., tf.Tensor]) -> Callable[..., tf.Tensor]: @wraps(old_fun) def new_fun(X: TensorType, **Ys: TensorType) -> tf.Tensor: Xs = tf.unstack(tf.expand_dims(X, axis=-2), axis=-1) fun_eval = old_fun(*Xs, **Ys) return tf.cond( pred=tf.less(tf.rank(fun_eval), tf.rank(X)), true_fn=lambda: fun_eval[..., tf.newaxis], false_fn=lambda: fun_eval, ) return new_fun if isinstance(funcs, Iterable): funcs = [wrapper(f) for f in funcs] else: funcs = wrapper(funcs) quadrature = NDiagGHQuadrature(dim, n_gh) if logspace: result = quadrature.logspace(funcs, Fmu, Fvar, **Ys) else: result = quadrature(funcs, Fmu, Fvar, **Ys) if isinstance(result, list): result = [tf.reshape(r, shape) for r in result] else: result = tf.reshape(result, shape) return result
[docs]@check_shapes( "Fmu: [N, Din]", "Fvar: [N, Din]", "Ys.values(): [broadcast N, .]", "return: [broadcast n_funs, N, P]", ) def ndiag_mc( funcs: Union[Callable[..., tf.Tensor], Iterable[Callable[..., tf.Tensor]]], S: int, Fmu: TensorType, Fvar: TensorType, logspace: bool = False, epsilon: Optional[TensorType] = None, **Ys: TensorType, ) -> tf.Tensor: """ Computes N Gaussian expectation integrals of one or more functions using Monte Carlo samples. The Gaussians must be independent. `Fmu`, `Fvar`, `Ys` should all have same shape, with overall size `N`. :param funcs: the integrand(s): Callable or Iterable of Callables that operates elementwise :param S: number of Monte Carlo sampling points :param Fmu: array/tensor :param Fvar: array/tensor :param logspace: if True, funcs are the log-integrands and this calculates the log-expectation of exp(funcs) :param Ys: arrays/tensors; deterministic arguments to be passed by name :return: shape is the same as that of the first Fmu """ N, D = tf.shape(Fmu)[0], tf.shape(Fvar)[1] if epsilon is None: epsilon = tf.random.normal(shape=[S, N, D], dtype=default_float()) mc_x = Fmu[None, :, :] + tf.sqrt(Fvar[None, :, :]) * epsilon mc_Xr = tf.reshape(mc_x, (S * N, D)) for name, Y in Ys.items(): D_out = Y.shape[1] # we can't rely on broadcasting and need tiling mc_Yr = tf.tile(Y[None, ...], [S, 1, 1]) # [S, N, D]_out Ys[name] = tf.reshape(mc_Yr, (S * N, D_out)) # S * [N, _]out def eval_func(func: Callable[..., tf.Tensor]) -> tf.Tensor: feval = func(mc_Xr, **Ys) feval = tf.reshape(feval, (S, N, -1)) if logspace: log_S = tf.math.log(to_default_float(S)) return tf.reduce_logsumexp(feval, axis=0) - log_S # [N, D] else: return tf.reduce_mean(feval, axis=0) if isinstance(funcs, Iterable): return [eval_func(f) for f in funcs] else: return eval_func(funcs)