From e437e67d484279eb86ece8232d7ee28fb56c7c37 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 10:29:13 -0300 Subject: [PATCH 01/41] WIP: quadrature refactoring --- gpflow/gaussian_quadrature/__init__.py | 2 + gpflow/gaussian_quadrature/base.py | 31 ++++++++++ gpflow/gaussian_quadrature/gauss_hermite.py | 68 +++++++++++++++++++++ gpflow/quadrature.py | 38 +++++++++++- test.py | 13 ++++ 5 files changed, 151 insertions(+), 1 deletion(-) create mode 100644 gpflow/gaussian_quadrature/__init__.py create mode 100644 gpflow/gaussian_quadrature/base.py create mode 100644 gpflow/gaussian_quadrature/gauss_hermite.py create mode 100644 test.py diff --git a/gpflow/gaussian_quadrature/__init__.py b/gpflow/gaussian_quadrature/__init__.py new file mode 100644 index 000000000..22560966d --- /dev/null +++ b/gpflow/gaussian_quadrature/__init__.py @@ -0,0 +1,2 @@ +from .base import GaussianQuadrature +from .gauss_hermite import NDDiagGHQuadrature diff --git a/gpflow/gaussian_quadrature/base.py b/gpflow/gaussian_quadrature/base.py new file mode 100644 index 000000000..a5ecfd4a8 --- /dev/null +++ b/gpflow/gaussian_quadrature/base.py @@ -0,0 +1,31 @@ +from collections.abc import Iterable +import abc + +import tensorflow as tf + +class GaussianQuadrature(): + + @abc.abstractmethod + def _build_X_W(self, mean, var): + raise NotImplementedError + + def logspace(self, fun, mean, var, *args, **kwargs): + X, W = self._build_X_W(mean, var) + logW = tf.math.log(W) + if isinstance(fun, Iterable): + return [ + tf.reduce_logsumexp(f(X, *args, **kwargs) + logW, axis=0) + for f in fun + ] + return tf.reduce_logsumexp(fun(X, *args, **kwargs) + logW, axis=0) + + def __call__(self, fun, mean, var, *args, **kwargs): + X, W = self._build_X_W(mean, var) + if isinstance(fun, Iterable): + # Maybe this can be better implemented by concating [f1(X), f2(X), ...] + # and sum-reducing all at once + return [ + tf.reduce_sum(f(X, *args, **kwargs) * W, axis=0) + for f in fun + ] + return tf.reduce_sum(fun(X, *args, **kwargs) * W, axis=0) diff --git a/gpflow/gaussian_quadrature/gauss_hermite.py b/gpflow/gaussian_quadrature/gauss_hermite.py new file mode 100644 index 000000000..328c2f668 --- /dev/null +++ b/gpflow/gaussian_quadrature/gauss_hermite.py @@ -0,0 +1,68 @@ +import numpy as np +import tensorflow as tf + +from .base import GaussianQuadrature +from ..config import default_float + + +def gh_points_and_weights(n_gh: int): + r""" + Given the number of Gauss-Hermite points n_gh, + returns the points z and the weights dz to perform the following + uni-dimensional gaussian quadrature: + + q(x) = N(mu, sigma²) + + E_{X~q(x)}[f(X)] = ∫ f(x)q(x)dx = \sum f(mu + sigma*z_k)*dz_k + + :param n_gh: number of Gauss-Hermite points, integer + :returns: points z and weights dz to compute uni-dimensional gaussian expectation, tuple + """ + z, dz = np.polynomial.hermite.hermgauss(n_gh) + z = z*np.sqrt(2) + dz = dz/np.sqrt(np.pi) + z, dz = z.astype(default_float()), dz.astype(default_float()) + return z, dz + +def reshape_Z_dZ(Z: np.array, dZ: np.array, dim: int): + Z = np.swapaxes(np.array(np.meshgrid(*Z)), 0, -1).reshape(-1, dim) + dZ = np.swapaxes(np.array(np.meshgrid(*dZ)), 0, -1).reshape(-1, dim).\ + prod(axis=-1, keepdims=True) + return Z, dZ + +def ndgh_points_and_weights(dim: int, n_gh: int): + r""" + :param n_gh: number of Gauss-Hermite points, integer + :param dim: dimension of the multivariate normal, integer + :returns: points Z, with shape [n_gh**dim, dim], and weights dZ, with shape [n_gh**dim, 1] + """ + z, dz = gh_points_and_weights(n_gh) + Z = [z]*dim + dZ = [dz]*dim + return reshape_Z_dZ(Z, dZ, dim) + +class NDDiagGHQuadrature(GaussianQuadrature): + + def __init__(self, dim, n_gh): + Z, dZ = ndgh_points_and_weights(dim, n_gh) + self.n_gh_total = n_gh ** dim + self.Z = tf.convert_to_tensor(Z) + self.dZ = tf.convert_to_tensor(dZ) + + def _build_X_W(self, mean, var): + new_shape = (self.n_gh_total,) + tuple([1] * (len(mean.shape) - 1)) + (-1,) + Z = tf.reshape(self.Z, new_shape) + dZ = tf.reshape(self.dZ, new_shape) + # Z : [n_gh_total] + [1]*len(batch) + [dims], typically [n_gh_total, 1, dims] + # dZ: [n_gh_total] + [1]*len(batch) + [1], typically [n_gh_total, 1, 1] + + mean = tf.expand_dims(mean, 0) + stddev = tf.expand_dims(tf.sqrt(var), 0) + # mean, stddev: [1] + batch + [dims], typically [1, N, dims] + + X = mean + stddev*Z + W = dZ + # X: [n_gh_total] + batch + [dims], typically [n_gh_total, N, dims] + # W: [n_gh_total] + [1]*len(batch) + [1], typically [n_gh_total, 1, 1] + + return X, W diff --git a/gpflow/quadrature.py b/gpflow/quadrature.py index bc9bdbb38..e4ca6dfb3 100644 --- a/gpflow/quadrature.py +++ b/gpflow/quadrature.py @@ -21,6 +21,7 @@ from .config import default_float from .utilities import to_default_float +from .gaussian_quadrature import NDDiagGHQuadrature def hermgauss(n: int): x, w = np.polynomial.hermite.hermgauss(n) @@ -97,7 +98,7 @@ def mvnquad(func, means, covs, H: int, Din: int = None, Dout=None): return tf.reduce_sum(fX * wr, 0) -def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): +def ndiagquad_old(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): """ Computes N Gaussian expectation integrals of one or more functions using Gauss-Hermite quadrature. The Gaussians must be independent. @@ -164,6 +165,41 @@ def eval_func(f): return eval_func(funcs) +def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): + """ + 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. + + :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 + + Fmu, Fvar, Ys should all have same shape, with overall size `N` + :return: shape is the same as that of the first Fmu + """ + print('Using new quadrature') + n_gh = H + if isinstance(Fmu, (tuple, list)): + dim = len(Fmu) + Fmu = tf.stack(Fmu, axis=-1) + Fvar = tf.stack(Fvar, axis=-1) + else: + dim = 1 + + quadrature = NDDiagGHQuadrature(dim, n_gh) + if logspace: + return quadrature.logspace(funcs, Fmu, Fvar, **Ys) + return quadrature(funcs, Fmu, Fvar, **Ys) + + def ndiag_mc(funcs, S: int, Fmu, Fvar, logspace: bool = False, epsilon=None, **Ys): """ Computes N Gaussian expectation integrals of one or more functions diff --git a/test.py b/test.py new file mode 100644 index 000000000..0ef2814f8 --- /dev/null +++ b/test.py @@ -0,0 +1,13 @@ +import numpy as np +import tensorflow as tf +import gpflow as gpf + +N = 100 +D = 1 + +Fmu = np.zeros((N, D)) * 1. +Fvar= np.ones((N, D)) * 1. +Y = np.ones((N, D)) * 1. + +lik = gpf.likelihoods.Bernoulli(invlink=tf.sigmoid) +print(lik.variational_expectations(Fmu, Fvar, Y).numpy()) From 138b41b270af954ebd6f211ac69d8c1cc19e6726 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 10:38:22 -0300 Subject: [PATCH 02/41] Removing old ndiagquad code --- gpflow/quadrature.py | 68 -------------------------------------------- 1 file changed, 68 deletions(-) diff --git a/gpflow/quadrature.py b/gpflow/quadrature.py index e4ca6dfb3..323659926 100644 --- a/gpflow/quadrature.py +++ b/gpflow/quadrature.py @@ -98,73 +98,6 @@ def mvnquad(func, means, covs, H: int, Din: int = None, Dout=None): return tf.reduce_sum(fX * wr, 0) -def ndiagquad_old(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): - """ - 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. - - :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 - - Fmu, Fvar, Ys should all have same shape, with overall size `N` - :return: shape is the same as that of the first Fmu - """ - if isinstance(Fmu, (tuple, list)): - Din = len(Fmu) - - def unify(f_list): - """Stack a list of means/vars into a full block.""" - return tf.reshape( - tensor=tf.concat([tf.reshape(f, shape=(-1, 1)) for f in f_list], axis=1), - shape=(-1, 1, Din), - ) - - shape = tf.shape(Fmu[0]) - Fmu, Fvar = map(unify, [Fmu, Fvar]) # both [N, 1, Din] - else: - Din = 1 - shape = tf.shape(Fmu) - Fmu, Fvar = [tf.reshape(f, (-1, 1, 1)) for f in [Fmu, Fvar]] - - xn, wn = mvhermgauss(H, Din) - # xn: H**Din x Din, wn: H**Din - - gh_x = xn.reshape(1, -1, Din) # [1, H]**Din x Din - Xall = gh_x * tf.sqrt(2.0 * Fvar) + Fmu # [N, H]**Din x Din - Xs = [Xall[:, :, i] for i in range(Din)] # [N, H]**Din each - - gh_w = wn * np.pi ** (-0.5 * Din) # H**Din x 1 - - for name, Y in Ys.items(): - Y = tf.reshape(Y, (-1, 1)) - Y = tf.tile(Y, [1, H ** Din]) # broadcast Y to match X - # without the tiling, some calls such as tf.where() (in bernoulli) fail - Ys[name] = Y # now [N, H]**Din - - def eval_func(f): - feval = f(*Xs, **Ys) # f should be elementwise: return shape [N, H]**Din - if logspace: - log_gh_w = np.log(gh_w.reshape(1, -1)) - result = tf.reduce_logsumexp(feval + log_gh_w, axis=1) - else: - result = tf.linalg.matmul(feval, gh_w.reshape(-1, 1)) - return tf.reshape(result, shape) - - if isinstance(funcs, Iterable): - return [eval_func(f) for f in funcs] - - return eval_func(funcs) - - def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): """ Computes N Gaussian expectation integrals of one or more functions @@ -185,7 +118,6 @@ def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): Fmu, Fvar, Ys should all have same shape, with overall size `N` :return: shape is the same as that of the first Fmu """ - print('Using new quadrature') n_gh = H if isinstance(Fmu, (tuple, list)): dim = len(Fmu) From 1d9f5b97aeff5ac84478a14d171cb04409363eea Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 10:39:04 -0300 Subject: [PATCH 03/41] deleted test code --- test.py | 13 ------------- 1 file changed, 13 deletions(-) delete mode 100644 test.py diff --git a/test.py b/test.py deleted file mode 100644 index 0ef2814f8..000000000 --- a/test.py +++ /dev/null @@ -1,13 +0,0 @@ -import numpy as np -import tensorflow as tf -import gpflow as gpf - -N = 100 -D = 1 - -Fmu = np.zeros((N, D)) * 1. -Fvar= np.ones((N, D)) * 1. -Y = np.ones((N, D)) * 1. - -lik = gpf.likelihoods.Bernoulli(invlink=tf.sigmoid) -print(lik.variational_expectations(Fmu, Fvar, Y).numpy()) From 4e2233798b4a730936f3210e8cac69285d580e1c Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 10:42:02 -0300 Subject: [PATCH 04/41] formatting and type-hint --- gpflow/gaussian_quadrature/gauss_hermite.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/gpflow/gaussian_quadrature/gauss_hermite.py b/gpflow/gaussian_quadrature/gauss_hermite.py index 328c2f668..5038943d5 100644 --- a/gpflow/gaussian_quadrature/gauss_hermite.py +++ b/gpflow/gaussian_quadrature/gauss_hermite.py @@ -10,11 +10,11 @@ def gh_points_and_weights(n_gh: int): Given the number of Gauss-Hermite points n_gh, returns the points z and the weights dz to perform the following uni-dimensional gaussian quadrature: - + q(x) = N(mu, sigma²) - + E_{X~q(x)}[f(X)] = ∫ f(x)q(x)dx = \sum f(mu + sigma*z_k)*dz_k - + :param n_gh: number of Gauss-Hermite points, integer :returns: points z and weights dz to compute uni-dimensional gaussian expectation, tuple """ @@ -43,7 +43,7 @@ def ndgh_points_and_weights(dim: int, n_gh: int): class NDDiagGHQuadrature(GaussianQuadrature): - def __init__(self, dim, n_gh): + def __init__(self, dim: int, n_gh: int): Z, dZ = ndgh_points_and_weights(dim, n_gh) self.n_gh_total = n_gh ** dim self.Z = tf.convert_to_tensor(Z) From 941548d9b11a4dc857166b21a23b7d0ba600cb8b Mon Sep 17 00:00:00 2001 From: ST John Date: Mon, 8 Jun 2020 15:35:57 +0100 Subject: [PATCH 05/41] merge modules --- gpflow/{gaussian_quadrature => quadrature}/__init__.py | 1 + gpflow/{gaussian_quadrature => quadrature}/base.py | 0 gpflow/{quadrature.py => quadrature/deprecated.py} | 0 gpflow/{gaussian_quadrature => quadrature}/gauss_hermite.py | 0 4 files changed, 1 insertion(+) rename gpflow/{gaussian_quadrature => quadrature}/__init__.py (51%) rename gpflow/{gaussian_quadrature => quadrature}/base.py (100%) rename gpflow/{quadrature.py => quadrature/deprecated.py} (100%) rename gpflow/{gaussian_quadrature => quadrature}/gauss_hermite.py (100%) diff --git a/gpflow/gaussian_quadrature/__init__.py b/gpflow/quadrature/__init__.py similarity index 51% rename from gpflow/gaussian_quadrature/__init__.py rename to gpflow/quadrature/__init__.py index 22560966d..91ed2f0f3 100644 --- a/gpflow/gaussian_quadrature/__init__.py +++ b/gpflow/quadrature/__init__.py @@ -1,2 +1,3 @@ from .base import GaussianQuadrature from .gauss_hermite import NDDiagGHQuadrature +from .deprecated import hermgauss, mvhermgauss, mvnquad, ndiagquad, ndiag_mc diff --git a/gpflow/gaussian_quadrature/base.py b/gpflow/quadrature/base.py similarity index 100% rename from gpflow/gaussian_quadrature/base.py rename to gpflow/quadrature/base.py diff --git a/gpflow/quadrature.py b/gpflow/quadrature/deprecated.py similarity index 100% rename from gpflow/quadrature.py rename to gpflow/quadrature/deprecated.py diff --git a/gpflow/gaussian_quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py similarity index 100% rename from gpflow/gaussian_quadrature/gauss_hermite.py rename to gpflow/quadrature/gauss_hermite.py From f5bc0a6be0042c655d29100c969968fd5cfa4bdf Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 12:06:16 -0300 Subject: [PATCH 06/41] black formatting --- gpflow/gaussian_quadrature/base.py | 12 +++--------- gpflow/gaussian_quadrature/gauss_hermite.py | 21 ++++++++++++--------- gpflow/quadrature.py | 1 + 3 files changed, 16 insertions(+), 18 deletions(-) diff --git a/gpflow/gaussian_quadrature/base.py b/gpflow/gaussian_quadrature/base.py index a5ecfd4a8..4e1aa5ff2 100644 --- a/gpflow/gaussian_quadrature/base.py +++ b/gpflow/gaussian_quadrature/base.py @@ -3,8 +3,8 @@ import tensorflow as tf -class GaussianQuadrature(): +class GaussianQuadrature: @abc.abstractmethod def _build_X_W(self, mean, var): raise NotImplementedError @@ -13,10 +13,7 @@ def logspace(self, fun, mean, var, *args, **kwargs): X, W = self._build_X_W(mean, var) logW = tf.math.log(W) if isinstance(fun, Iterable): - return [ - tf.reduce_logsumexp(f(X, *args, **kwargs) + logW, axis=0) - for f in fun - ] + return [tf.reduce_logsumexp(f(X, *args, **kwargs) + logW, axis=0) for f in fun] return tf.reduce_logsumexp(fun(X, *args, **kwargs) + logW, axis=0) def __call__(self, fun, mean, var, *args, **kwargs): @@ -24,8 +21,5 @@ def __call__(self, fun, mean, var, *args, **kwargs): if isinstance(fun, Iterable): # Maybe this can be better implemented by concating [f1(X), f2(X), ...] # and sum-reducing all at once - return [ - tf.reduce_sum(f(X, *args, **kwargs) * W, axis=0) - for f in fun - ] + return [tf.reduce_sum(f(X, *args, **kwargs) * W, axis=0) for f in fun] return tf.reduce_sum(fun(X, *args, **kwargs) * W, axis=0) diff --git a/gpflow/gaussian_quadrature/gauss_hermite.py b/gpflow/gaussian_quadrature/gauss_hermite.py index 5038943d5..2c29f6f29 100644 --- a/gpflow/gaussian_quadrature/gauss_hermite.py +++ b/gpflow/gaussian_quadrature/gauss_hermite.py @@ -19,17 +19,20 @@ def gh_points_and_weights(n_gh: int): :returns: points z and weights dz to compute uni-dimensional gaussian expectation, tuple """ z, dz = np.polynomial.hermite.hermgauss(n_gh) - z = z*np.sqrt(2) - dz = dz/np.sqrt(np.pi) + z = z * np.sqrt(2) + dz = dz / np.sqrt(np.pi) z, dz = z.astype(default_float()), dz.astype(default_float()) return z, dz + def reshape_Z_dZ(Z: np.array, dZ: np.array, dim: int): Z = np.swapaxes(np.array(np.meshgrid(*Z)), 0, -1).reshape(-1, dim) - dZ = np.swapaxes(np.array(np.meshgrid(*dZ)), 0, -1).reshape(-1, dim).\ - prod(axis=-1, keepdims=True) + dZ = ( + np.swapaxes(np.array(np.meshgrid(*dZ)), 0, -1).reshape(-1, dim).prod(axis=-1, keepdims=True) + ) return Z, dZ + def ndgh_points_and_weights(dim: int, n_gh: int): r""" :param n_gh: number of Gauss-Hermite points, integer @@ -37,12 +40,12 @@ def ndgh_points_and_weights(dim: int, n_gh: int): :returns: points Z, with shape [n_gh**dim, dim], and weights dZ, with shape [n_gh**dim, 1] """ z, dz = gh_points_and_weights(n_gh) - Z = [z]*dim - dZ = [dz]*dim + Z = [z] * dim + dZ = [dz] * dim return reshape_Z_dZ(Z, dZ, dim) -class NDDiagGHQuadrature(GaussianQuadrature): +class NDDiagGHQuadrature(GaussianQuadrature): def __init__(self, dim: int, n_gh: int): Z, dZ = ndgh_points_and_weights(dim, n_gh) self.n_gh_total = n_gh ** dim @@ -50,7 +53,7 @@ def __init__(self, dim: int, n_gh: int): self.dZ = tf.convert_to_tensor(dZ) def _build_X_W(self, mean, var): - new_shape = (self.n_gh_total,) + tuple([1] * (len(mean.shape) - 1)) + (-1,) + new_shape = (self.n_gh_total,) + tuple([1] * (len(mean.shape) - 1)) + (-1,) Z = tf.reshape(self.Z, new_shape) dZ = tf.reshape(self.dZ, new_shape) # Z : [n_gh_total] + [1]*len(batch) + [dims], typically [n_gh_total, 1, dims] @@ -60,7 +63,7 @@ def _build_X_W(self, mean, var): stddev = tf.expand_dims(tf.sqrt(var), 0) # mean, stddev: [1] + batch + [dims], typically [1, N, dims] - X = mean + stddev*Z + X = mean + stddev * Z W = dZ # X: [n_gh_total] + batch + [dims], typically [n_gh_total, N, dims] # W: [n_gh_total] + [1]*len(batch) + [1], typically [n_gh_total, 1, 1] diff --git a/gpflow/quadrature.py b/gpflow/quadrature.py index 323659926..a31b3cedd 100644 --- a/gpflow/quadrature.py +++ b/gpflow/quadrature.py @@ -23,6 +23,7 @@ from .gaussian_quadrature import NDDiagGHQuadrature + def hermgauss(n: int): x, w = np.polynomial.hermite.hermgauss(n) x, w = x.astype(default_float()), w.astype(default_float()) From 24123091c0bc35c9fc82442240acc2536851ad88 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 12:10:10 -0300 Subject: [PATCH 07/41] formatting --- gpflow/quadrature/base.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gpflow/quadrature/base.py b/gpflow/quadrature/base.py index 4e1aa5ff2..486aae566 100644 --- a/gpflow/quadrature/base.py +++ b/gpflow/quadrature/base.py @@ -4,7 +4,8 @@ import tensorflow as tf -class GaussianQuadrature: +class GaussianQuadrature(): + @abc.abstractmethod def _build_X_W(self, mean, var): raise NotImplementedError From 0a64680e6d1c387c1b112898d5f1abed35c19528 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 13:44:24 -0300 Subject: [PATCH 08/41] solving failing tests --- gpflow/quadrature/__init__.py | 4 +- gpflow/quadrature/base.py | 7 +- gpflow/quadrature/deprecated.py | 175 ++++++------------------------ gpflow/quadrature/quadrature.py | 187 ++++++++++++++++++++++++++++++++ 4 files changed, 227 insertions(+), 146 deletions(-) create mode 100644 gpflow/quadrature/quadrature.py diff --git a/gpflow/quadrature/__init__.py b/gpflow/quadrature/__init__.py index 91ed2f0f3..663f98c7e 100644 --- a/gpflow/quadrature/__init__.py +++ b/gpflow/quadrature/__init__.py @@ -1,3 +1,5 @@ from .base import GaussianQuadrature from .gauss_hermite import NDDiagGHQuadrature -from .deprecated import hermgauss, mvhermgauss, mvnquad, ndiagquad, ndiag_mc +from .quadrature import hermgauss, mvhermgauss, mvnquad, ndiagquad, ndiag_mc + +from .deprecated import ndiagquad as ndiagquad_old diff --git a/gpflow/quadrature/base.py b/gpflow/quadrature/base.py index 486aae566..1b6453461 100644 --- a/gpflow/quadrature/base.py +++ b/gpflow/quadrature/base.py @@ -4,8 +4,7 @@ import tensorflow as tf -class GaussianQuadrature(): - +class GaussianQuadrature: @abc.abstractmethod def _build_X_W(self, mean, var): raise NotImplementedError @@ -20,7 +19,9 @@ def logspace(self, fun, mean, var, *args, **kwargs): def __call__(self, fun, mean, var, *args, **kwargs): X, W = self._build_X_W(mean, var) if isinstance(fun, Iterable): - # Maybe this can be better implemented by concating [f1(X), f2(X), ...] + # Maybe this can be better implemented by stacking [f1(X), f2(X), ...] # and sum-reducing all at once + # The problem: there is no garantee that f1(X), f2(X), ... + # have comaptible shapes return [tf.reduce_sum(f(X, *args, **kwargs) * W, axis=0) for f in fun] return tf.reduce_sum(fun(X, *args, **kwargs) * W, axis=0) diff --git a/gpflow/quadrature/deprecated.py b/gpflow/quadrature/deprecated.py index a31b3cedd..1368b89ca 100644 --- a/gpflow/quadrature/deprecated.py +++ b/gpflow/quadrature/deprecated.py @@ -1,112 +1,17 @@ -# Copyright 2017-2018 the GPflow authors. -# -# 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 -from collections.abc import Iterable - import numpy as np import tensorflow as tf -from .config import default_float -from .utilities import to_default_float - -from .gaussian_quadrature import NDDiagGHQuadrature - - -def hermgauss(n: int): - x, w = np.polynomial.hermite.hermgauss(n) - x, w = x.astype(default_float()), w.astype(default_float()) - return x, w - - -def mvhermgauss(H: int, D: int): - """ - 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 = 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 - - -def mvnquad(func, means, covs, H: int, Din: int = None, Dout=None): - """ - 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 means: NxD - :param covs: NxDxD - :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 (N,*Dout) - """ - # 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:]) +from collections.abc import Iterable - 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) +from .quadrature import mvhermgauss def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): """ 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. - :param funcs: the integrand(s): Callable or Iterable of Callables that operates elementwise :param H: number of Gauss-Hermite quadrature points @@ -115,65 +20,51 @@ def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): :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 - Fmu, Fvar, Ys should all have same shape, with overall size `N` :return: shape is the same as that of the first Fmu """ - n_gh = H if isinstance(Fmu, (tuple, list)): - dim = len(Fmu) - Fmu = tf.stack(Fmu, axis=-1) - Fvar = tf.stack(Fvar, axis=-1) - else: - dim = 1 - - quadrature = NDDiagGHQuadrature(dim, n_gh) - if logspace: - return quadrature.logspace(funcs, Fmu, Fvar, **Ys) - return quadrature(funcs, Fmu, Fvar, **Ys) + Din = len(Fmu) + def unify(f_list): + """Stack a list of means/vars into a full block.""" + return tf.reshape( + tensor=tf.concat([tf.reshape(f, shape=(-1, 1)) for f in f_list], axis=1), + shape=(-1, 1, Din), + ) -def ndiag_mc(funcs, S: int, Fmu, Fvar, logspace: bool = False, epsilon=None, **Ys): - """ - Computes N Gaussian expectation integrals of one or more functions - using Monte Carlo samples. The Gaussians must be independent. - - :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 + shape = tf.shape(Fmu[0]) + Fmu, Fvar = map(unify, [Fmu, Fvar]) # both [N, 1, Din] + else: + Din = 1 + shape = tf.shape(Fmu) + Fmu, Fvar = [tf.reshape(f, (-1, 1, 1)) for f in [Fmu, Fvar]] - Fmu, Fvar, Ys should all have same shape, with overall size `N` - :return: shape is the same as that of the first Fmu - """ - N, D = Fmu.shape[0], Fvar.shape[1] + xn, wn = mvhermgauss(H, Din) + # xn: H**Din x Din, wn: H**Din - if epsilon is None: - epsilon = tf.random.normal((S, N, D), dtype=default_float()) + gh_x = xn.reshape(1, -1, Din) # [1, H]**Din x Din + Xall = gh_x * tf.sqrt(2.0 * Fvar) + Fmu # [N, H]**Din x Din + Xs = [Xall[:, :, i] for i in range(Din)] # [N, H]**Din each - mc_x = Fmu[None, :, :] + tf.sqrt(Fvar[None, :, :]) * epsilon - mc_Xr = tf.reshape(mc_x, (S * N, D)) + gh_w = wn * np.pi ** (-0.5 * Din) # H**Din x 1 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 + Y = tf.reshape(Y, (-1, 1)) + Y = tf.tile(Y, [1, H ** Din]) # broadcast Y to match X + # without the tiling, some calls such as tf.where() (in bernoulli) fail + Ys[name] = Y # now [N, H]**Din - def eval_func(func): - feval = func(mc_Xr, **Ys) - feval = tf.reshape(feval, (S, N, -1)) + def eval_func(f): + feval = f(*Xs, **Ys) # f should be elementwise: return shape [N, H]**Din if logspace: - log_S = tf.math.log(to_default_float(S)) - return tf.reduce_logsumexp(feval, axis=0) - log_S # [N, D] + log_gh_w = np.log(gh_w.reshape(1, -1)) + result = tf.reduce_logsumexp(feval + log_gh_w, axis=1) else: - return tf.reduce_mean(feval, axis=0) + result = tf.linalg.matmul(feval, gh_w.reshape(-1, 1)) + return tf.reshape(result, shape) if isinstance(funcs, Iterable): return [eval_func(f) for f in funcs] - else: - return eval_func(funcs) + + return eval_func(funcs) diff --git a/gpflow/quadrature/quadrature.py b/gpflow/quadrature/quadrature.py new file mode 100644 index 000000000..34085eaa9 --- /dev/null +++ b/gpflow/quadrature/quadrature.py @@ -0,0 +1,187 @@ +# Copyright 2017-2018 the GPflow authors. +# +# 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 +from collections.abc import Iterable + +import numpy as np +import tensorflow as tf + +from ..config import default_float +from ..utilities import to_default_float + +from .gauss_hermite import NDDiagGHQuadrature + + +def hermgauss(n: int): + x, w = np.polynomial.hermite.hermgauss(n) + x, w = x.astype(default_float()), w.astype(default_float()) + return x, w + + +def mvhermgauss(H: int, D: int): + """ + 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 = 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 + + +def mvnquad(func, means, covs, H: int, Din: int = None, Dout=None): + """ + 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 means: NxD + :param covs: NxDxD + :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 (N,*Dout) + """ + # 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) + + +def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): + """ + 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. + + :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 + + Fmu, Fvar, Ys should all have same shape, with overall size `N` + :return: shape is the same as that of the first Fmu + """ + 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) + + quadrature = NDDiagGHQuadrature(dim, n_gh) + + if logspace: + result = quadrature.logspace(funcs, Fmu, Fvar, **Ys) + else: + result = quadrature(funcs, Fmu, Fvar, **Ys) + + if isinstance(result, (tuple, list)): + return [tf.reshape(r, shape) for r in result] + return tf.reshape(result, shape) + + +def ndiag_mc(funcs, S: int, Fmu, Fvar, logspace: bool = False, epsilon=None, **Ys): + """ + Computes N Gaussian expectation integrals of one or more functions + using Monte Carlo samples. The Gaussians must be independent. + + :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 + + Fmu, Fvar, Ys should all have same shape, with overall size `N` + :return: shape is the same as that of the first Fmu + """ + N, D = Fmu.shape[0], Fvar.shape[1] + + if epsilon is None: + epsilon = tf.random.normal((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): + 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) From 4592ed8a04c53920cc9b7ccd4360bf2fd963c2ea Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 15:19:08 -0300 Subject: [PATCH 09/41] fixing failing tests --- gpflow/quadrature/deprecated.py | 3 +++ gpflow/quadrature/quadrature.py | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/gpflow/quadrature/deprecated.py b/gpflow/quadrature/deprecated.py index 1368b89ca..45c06a0e0 100644 --- a/gpflow/quadrature/deprecated.py +++ b/gpflow/quadrature/deprecated.py @@ -35,6 +35,9 @@ def unify(f_list): shape = tf.shape(Fmu[0]) Fmu, Fvar = map(unify, [Fmu, Fvar]) # both [N, 1, Din] + + print(Fmu) + print(Fvar) else: Din = 1 shape = tf.shape(Fmu) diff --git a/gpflow/quadrature/quadrature.py b/gpflow/quadrature/quadrature.py index 34085eaa9..e741410fb 100644 --- a/gpflow/quadrature/quadrature.py +++ b/gpflow/quadrature/quadrature.py @@ -126,8 +126,8 @@ def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): Fmu = tf.stack(Fmu, axis=-1) Fvar = tf.stack(Fvar, axis=-1) else: - dim = 1 - shape = tf.shape(Fmu) + dim = tf.shape(Fmu)[-1] + shape = tf.shape(Fmu)[:-1] quadrature = NDDiagGHQuadrature(dim, n_gh) From d1afb320faf3cf72883ddb3e91fe6d4fa93daf38 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 16:20:48 -0300 Subject: [PATCH 10/41] fixes --- gpflow/quadrature/base.py | 2 +- gpflow/quadrature/quadrature.py | 8 ++------ 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/gpflow/quadrature/base.py b/gpflow/quadrature/base.py index 1b6453461..f6ffadd4c 100644 --- a/gpflow/quadrature/base.py +++ b/gpflow/quadrature/base.py @@ -19,7 +19,7 @@ def logspace(self, fun, mean, var, *args, **kwargs): def __call__(self, fun, mean, var, *args, **kwargs): X, W = self._build_X_W(mean, var) if isinstance(fun, Iterable): - # Maybe this can be better implemented by stacking [f1(X), f2(X), ...] + # Maybe this can be better implemented by stacking [f1(X), f2(X), ...] # and sum-reducing all at once # The problem: there is no garantee that f1(X), f2(X), ... # have comaptible shapes diff --git a/gpflow/quadrature/quadrature.py b/gpflow/quadrature/quadrature.py index e741410fb..eaf95fa59 100644 --- a/gpflow/quadrature/quadrature.py +++ b/gpflow/quadrature/quadrature.py @@ -122,12 +122,10 @@ def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): 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 = tf.shape(Fmu)[-1] - shape = tf.shape(Fmu)[:-1] + dim = int(tf.shape(Fmu)[-1]) quadrature = NDDiagGHQuadrature(dim, n_gh) @@ -136,9 +134,7 @@ def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): else: result = quadrature(funcs, Fmu, Fvar, **Ys) - if isinstance(result, (tuple, list)): - return [tf.reshape(r, shape) for r in result] - return tf.reshape(result, shape) + return result def ndiag_mc(funcs, S: int, Fmu, Fvar, logspace: bool = False, epsilon=None, **Ys): From 8ab2f4a83304dcc1be20961aee78f6632cadb74e Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 16:21:04 -0300 Subject: [PATCH 11/41] adapting tests for new syntax, keeping numerical behavior --- tests/gpflow/test_quadrature.py | 58 ++++++++++++++++++++++++++------- 1 file changed, 46 insertions(+), 12 deletions(-) diff --git a/tests/gpflow/test_quadrature.py b/tests/gpflow/test_quadrature.py index 55e36754a..d5b396f07 100644 --- a/tests/gpflow/test_quadrature.py +++ b/tests/gpflow/test_quadrature.py @@ -24,7 +24,13 @@ @pytest.mark.parametrize("var", [np.array([3.0, 3.5])]) def test_diagquad_1d(mu, var): num_gauss_hermite_points = 25 - quad = quadrature.ndiagquad([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, [mu], [var]) + + mu = mu[:, None] + var = var[:, None] + + # quad = quadrature.ndiagquad([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, [mu], [var]) + quad = quadrature.ndiagquad([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, mu, var) + expected = np.exp(mu + var / 2) assert_allclose(quad[0], expected) @@ -37,10 +43,19 @@ def test_diagquad_2d(mu1, var1, mu2, var2): alpha = 2.5 # using logspace=True we can reduce this, see test_diagquad_logspace num_gauss_hermite_points = 35 - quad = quadrature.ndiagquad( - lambda *X: tf.exp(X[0] + alpha * X[1]), num_gauss_hermite_points, [mu1, mu2], [var1, var2], - ) - expected = np.exp(mu1 + var1 / 2 + alpha * mu2 + alpha ** 2 * var2 / 2) + + mu = np.hstack([x[:, None] for x in (mu1, mu2)]) + var = np.hstack([x[:, None] for x in (var1, var2)]) + + # fun = lambda *X: tf.exp(X[0] + alpha * X[1]) + fun = lambda X: tf.exp(X[..., 0:1] + alpha * X[..., 1:2]) + + # quad = quadrature.ndiagquad(fun, num_gauss_hermite_points, [mu1, mu2], [var1, var2],) + quad = quadrature.ndiagquad(fun, num_gauss_hermite_points, mu, var,) + + # expected = np.exp(mu1 + var1 / 2 + alpha * mu2 + alpha ** 2 * var2 / 2) + expected = np.exp(mu[..., 0:1] + var[..., 0:1] / 2 + alpha * mu[..., 1:2] + alpha ** 2 * var[..., 1:2] / 2) + assert_allclose(quad, expected) @@ -51,14 +66,23 @@ def test_diagquad_2d(mu1, var1, mu2, var2): def test_diagquad_logspace(mu1, var1, mu2, var2): alpha = 2.5 num_gauss_hermite_points = 25 + + mu = np.hstack([x[:, None] for x in (mu1, mu2)]) + var = np.hstack([x[:, None] for x in (var1, var2)]) + + # fun = lambda *X: (X[0] + alpha * X[1]) + fun = lambda X: (X[..., 0:1] + alpha * X[..., 1:2]) + + # quad = quadrature.ndiagquad( + # fun, num_gauss_hermite_points, [mu1, mu2], [var1, var2], logspace=True, + # ) quad = quadrature.ndiagquad( - lambda *X: (X[0] + alpha * X[1]), - num_gauss_hermite_points, - [mu1, mu2], - [var1, var2], - logspace=True, + fun, num_gauss_hermite_points, mu, var, logspace=True, ) - expected = mu1 + var1 / 2 + alpha * mu2 + alpha ** 2 * var2 / 2 + + # expected = mu1 + var1 / 2 + alpha * mu2 + alpha ** 2 * var2 / 2 + expected = mu[..., 0:1] + var[..., 0:1] / 2 + alpha * mu[..., 1:2] + alpha ** 2 * var[..., 1:2] / 2 + assert_allclose(quad, expected) @@ -67,8 +91,18 @@ def test_diagquad_logspace(mu1, var1, mu2, var2): def test_diagquad_with_kwarg(mu1, var1): alpha = np.array([2.5, -1.3]) num_gauss_hermite_points = 25 + + mu1 = mu1[:, None] + var1 = var1[:, None] + alpha = alpha[:, None] + + fun = lambda X, Y: tf.exp(X * Y) quad = quadrature.ndiagquad( - lambda X, Y: tf.exp(X * Y), num_gauss_hermite_points, mu1, var1, Y=alpha + fun, num_gauss_hermite_points, mu1, var1, Y=alpha ) + expected = np.exp(alpha * mu1 + alpha ** 2 * var1 / 2) + + print(expected) + assert_allclose(quad, expected) From 98428db5c78288681cec70c1bddc36376ccafa87 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 16:22:01 -0300 Subject: [PATCH 12/41] black formatting --- tests/gpflow/test_quadrature.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/gpflow/test_quadrature.py b/tests/gpflow/test_quadrature.py index d5b396f07..447b25d67 100644 --- a/tests/gpflow/test_quadrature.py +++ b/tests/gpflow/test_quadrature.py @@ -54,7 +54,9 @@ def test_diagquad_2d(mu1, var1, mu2, var2): quad = quadrature.ndiagquad(fun, num_gauss_hermite_points, mu, var,) # expected = np.exp(mu1 + var1 / 2 + alpha * mu2 + alpha ** 2 * var2 / 2) - expected = np.exp(mu[..., 0:1] + var[..., 0:1] / 2 + alpha * mu[..., 1:2] + alpha ** 2 * var[..., 1:2] / 2) + expected = np.exp( + mu[..., 0:1] + var[..., 0:1] / 2 + alpha * mu[..., 1:2] + alpha ** 2 * var[..., 1:2] / 2 + ) assert_allclose(quad, expected) @@ -76,12 +78,12 @@ def test_diagquad_logspace(mu1, var1, mu2, var2): # quad = quadrature.ndiagquad( # fun, num_gauss_hermite_points, [mu1, mu2], [var1, var2], logspace=True, # ) - quad = quadrature.ndiagquad( - fun, num_gauss_hermite_points, mu, var, logspace=True, - ) + quad = quadrature.ndiagquad(fun, num_gauss_hermite_points, mu, var, logspace=True,) # expected = mu1 + var1 / 2 + alpha * mu2 + alpha ** 2 * var2 / 2 - expected = mu[..., 0:1] + var[..., 0:1] / 2 + alpha * mu[..., 1:2] + alpha ** 2 * var[..., 1:2] / 2 + expected = ( + mu[..., 0:1] + var[..., 0:1] / 2 + alpha * mu[..., 1:2] + alpha ** 2 * var[..., 1:2] / 2 + ) assert_allclose(quad, expected) @@ -97,9 +99,7 @@ def test_diagquad_with_kwarg(mu1, var1): alpha = alpha[:, None] fun = lambda X, Y: tf.exp(X * Y) - quad = quadrature.ndiagquad( - fun, num_gauss_hermite_points, mu1, var1, Y=alpha - ) + quad = quadrature.ndiagquad(fun, num_gauss_hermite_points, mu1, var1, Y=alpha) expected = np.exp(alpha * mu1 + alpha ** 2 * var1 / 2) From 695d87fe888bbcc6edc0694ab8be7e027e4bb6e2 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 20:03:08 -0300 Subject: [PATCH 13/41] remove printf --- tests/gpflow/test_quadrature.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/tests/gpflow/test_quadrature.py b/tests/gpflow/test_quadrature.py index 447b25d67..c895aa55c 100644 --- a/tests/gpflow/test_quadrature.py +++ b/tests/gpflow/test_quadrature.py @@ -102,7 +102,4 @@ def test_diagquad_with_kwarg(mu1, var1): quad = quadrature.ndiagquad(fun, num_gauss_hermite_points, mu1, var1, Y=alpha) expected = np.exp(alpha * mu1 + alpha ** 2 * var1 / 2) - - print(expected) - assert_allclose(quad, expected) From 62b79f85339fbc730db87c8030da424d50269086 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 20:44:29 -0300 Subject: [PATCH 14/41] changed code for compiled tf compatibility --- .../notebooks/advanced/ordinal_regression.pct.py | 12 ++++++++++-- gpflow/quadrature/gauss_hermite.py | 16 +++++++--------- gpflow/quadrature/quadrature.py | 2 +- 3 files changed, 18 insertions(+), 12 deletions(-) diff --git a/doc/source/notebooks/advanced/ordinal_regression.pct.py b/doc/source/notebooks/advanced/ordinal_regression.pct.py index 2f85f204f..b762cde06 100644 --- a/doc/source/notebooks/advanced/ordinal_regression.pct.py +++ b/doc/source/notebooks/advanced/ordinal_regression.pct.py @@ -6,7 +6,7 @@ # extension: .py # format_name: percent # format_version: '1.3' -# jupytext_version: 1.3.3 +# jupytext_version: 1.4.2 # kernelspec: # display_name: Python 3 # language: python @@ -74,8 +74,16 @@ def generate_data(num_data): likelihood = gpflow.likelihoods.Ordinal(bin_edges) # build a model with this likelihood -m = gpflow.models.VGP(data=(X, Y), kernel=gpflow.kernels.Matern32(), likelihood=likelihood) +m = gpflow.models.VGP( + data=(tf.convert_to_tensor(X), tf.convert_to_tensor(Y)), + kernel=gpflow.kernels.Matern32(), + likelihood=likelihood +) + +# %% +tf.cast(tf.shape(tf.ones((10,3)))[-1], tf.int32) +# %% # fit the model opt = gpflow.optimizers.Scipy() opt.minimize(m.training_loss, m.trainable_variables, options=dict(maxiter=100)) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index 2c29f6f29..44e085fbd 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -22,14 +22,12 @@ def gh_points_and_weights(n_gh: int): z = z * np.sqrt(2) dz = dz / np.sqrt(np.pi) z, dz = z.astype(default_float()), dz.astype(default_float()) - return z, dz + return tf.convert_to_tensor(z), tf.convert_to_tensor(dz) -def reshape_Z_dZ(Z: np.array, dZ: np.array, dim: int): - Z = np.swapaxes(np.array(np.meshgrid(*Z)), 0, -1).reshape(-1, dim) - dZ = ( - np.swapaxes(np.array(np.meshgrid(*dZ)), 0, -1).reshape(-1, dim).prod(axis=-1, keepdims=True) - ) +def reshape_Z_dZ(zs, dzs, dim: int): + Z = tf.reshape(tf.stack(tf.meshgrid(*zs), axis=-1), (-1, dim)) + dZ = tf.reduce_prod(tf.reshape(tf.stack(tf.meshgrid(*dzs), axis=-1), (-1, dim)), axis=-1, keepdims=True) return Z, dZ @@ -40,9 +38,9 @@ def ndgh_points_and_weights(dim: int, n_gh: int): :returns: points Z, with shape [n_gh**dim, dim], and weights dZ, with shape [n_gh**dim, 1] """ z, dz = gh_points_and_weights(n_gh) - Z = [z] * dim - dZ = [dz] * dim - return reshape_Z_dZ(Z, dZ, dim) + zs = tf.unstack(tf.repeat(tf.expand_dims(z, axis=0), dim, axis=0), axis=0) + dzs = tf.unstack(tf.repeat(tf.expand_dims(dz, axis=0), dim, axis=0), axis=0) + return reshape_Z_dZ(zs, dzs, dim) class NDDiagGHQuadrature(GaussianQuadrature): diff --git a/gpflow/quadrature/quadrature.py b/gpflow/quadrature/quadrature.py index eaf95fa59..016d843f5 100644 --- a/gpflow/quadrature/quadrature.py +++ b/gpflow/quadrature/quadrature.py @@ -125,7 +125,7 @@ def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): Fmu = tf.stack(Fmu, axis=-1) Fvar = tf.stack(Fvar, axis=-1) else: - dim = int(tf.shape(Fmu)[-1]) + dim = tf.shape(Fmu)[-1] quadrature = NDDiagGHQuadrature(dim, n_gh) From 20ad0bc8d82f60e7929ae8eb73e733439b002aa1 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 20:49:17 -0300 Subject: [PATCH 15/41] black --- gpflow/quadrature/gauss_hermite.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index 44e085fbd..aa349060b 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -25,9 +25,13 @@ def gh_points_and_weights(n_gh: int): return tf.convert_to_tensor(z), tf.convert_to_tensor(dz) -def reshape_Z_dZ(zs, dzs, dim: int): - Z = tf.reshape(tf.stack(tf.meshgrid(*zs), axis=-1), (-1, dim)) - dZ = tf.reduce_prod(tf.reshape(tf.stack(tf.meshgrid(*dzs), axis=-1), (-1, dim)), axis=-1, keepdims=True) +def list_to_flat_grid(xs): + return tf.reshape(tf.stack(tf.meshgrid(*xs), axis=-1), (-1, len(xs))) + + +def reshape_Z_dZ(zs, dzs): + Z = list_to_flat_grid(zs) + dZ = tf.reduce_prod(list_to_flat_grid(dzs), axis=-1, keepdims=True) return Z, dZ @@ -40,7 +44,7 @@ def ndgh_points_and_weights(dim: int, n_gh: int): z, dz = gh_points_and_weights(n_gh) zs = tf.unstack(tf.repeat(tf.expand_dims(z, axis=0), dim, axis=0), axis=0) dzs = tf.unstack(tf.repeat(tf.expand_dims(dz, axis=0), dim, axis=0), axis=0) - return reshape_Z_dZ(zs, dzs, dim) + return reshape_Z_dZ(zs, dzs) class NDDiagGHQuadrature(GaussianQuadrature): From d4161d8d6ab02ab91c2f6854dfdfc218cbfc57b5 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 8 Jun 2020 20:56:30 -0300 Subject: [PATCH 16/41] restored to original version --- .../notebooks/advanced/ordinal_regression.pct.py | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/doc/source/notebooks/advanced/ordinal_regression.pct.py b/doc/source/notebooks/advanced/ordinal_regression.pct.py index b762cde06..2f85f204f 100644 --- a/doc/source/notebooks/advanced/ordinal_regression.pct.py +++ b/doc/source/notebooks/advanced/ordinal_regression.pct.py @@ -6,7 +6,7 @@ # extension: .py # format_name: percent # format_version: '1.3' -# jupytext_version: 1.4.2 +# jupytext_version: 1.3.3 # kernelspec: # display_name: Python 3 # language: python @@ -74,16 +74,8 @@ def generate_data(num_data): likelihood = gpflow.likelihoods.Ordinal(bin_edges) # build a model with this likelihood -m = gpflow.models.VGP( - data=(tf.convert_to_tensor(X), tf.convert_to_tensor(Y)), - kernel=gpflow.kernels.Matern32(), - likelihood=likelihood -) - -# %% -tf.cast(tf.shape(tf.ones((10,3)))[-1], tf.int32) +m = gpflow.models.VGP(data=(X, Y), kernel=gpflow.kernels.Matern32(), likelihood=likelihood) -# %% # fit the model opt = gpflow.optimizers.Scipy() opt.minimize(m.training_loss, m.trainable_variables, options=dict(maxiter=100)) From 029ec7a0d1c3f22f259ab5f23250f88219c914b5 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 9 Jun 2020 13:54:52 -0300 Subject: [PATCH 17/41] undoing changes --- tests/gpflow/test_quadrature.py | 59 ++++++++------------------------- 1 file changed, 14 insertions(+), 45 deletions(-) diff --git a/tests/gpflow/test_quadrature.py b/tests/gpflow/test_quadrature.py index c895aa55c..55e36754a 100644 --- a/tests/gpflow/test_quadrature.py +++ b/tests/gpflow/test_quadrature.py @@ -24,13 +24,7 @@ @pytest.mark.parametrize("var", [np.array([3.0, 3.5])]) def test_diagquad_1d(mu, var): num_gauss_hermite_points = 25 - - mu = mu[:, None] - var = var[:, None] - - # quad = quadrature.ndiagquad([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, [mu], [var]) - quad = quadrature.ndiagquad([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, mu, var) - + quad = quadrature.ndiagquad([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, [mu], [var]) expected = np.exp(mu + var / 2) assert_allclose(quad[0], expected) @@ -43,21 +37,10 @@ def test_diagquad_2d(mu1, var1, mu2, var2): alpha = 2.5 # using logspace=True we can reduce this, see test_diagquad_logspace num_gauss_hermite_points = 35 - - mu = np.hstack([x[:, None] for x in (mu1, mu2)]) - var = np.hstack([x[:, None] for x in (var1, var2)]) - - # fun = lambda *X: tf.exp(X[0] + alpha * X[1]) - fun = lambda X: tf.exp(X[..., 0:1] + alpha * X[..., 1:2]) - - # quad = quadrature.ndiagquad(fun, num_gauss_hermite_points, [mu1, mu2], [var1, var2],) - quad = quadrature.ndiagquad(fun, num_gauss_hermite_points, mu, var,) - - # expected = np.exp(mu1 + var1 / 2 + alpha * mu2 + alpha ** 2 * var2 / 2) - expected = np.exp( - mu[..., 0:1] + var[..., 0:1] / 2 + alpha * mu[..., 1:2] + alpha ** 2 * var[..., 1:2] / 2 + quad = quadrature.ndiagquad( + lambda *X: tf.exp(X[0] + alpha * X[1]), num_gauss_hermite_points, [mu1, mu2], [var1, var2], ) - + expected = np.exp(mu1 + var1 / 2 + alpha * mu2 + alpha ** 2 * var2 / 2) assert_allclose(quad, expected) @@ -68,23 +51,14 @@ def test_diagquad_2d(mu1, var1, mu2, var2): def test_diagquad_logspace(mu1, var1, mu2, var2): alpha = 2.5 num_gauss_hermite_points = 25 - - mu = np.hstack([x[:, None] for x in (mu1, mu2)]) - var = np.hstack([x[:, None] for x in (var1, var2)]) - - # fun = lambda *X: (X[0] + alpha * X[1]) - fun = lambda X: (X[..., 0:1] + alpha * X[..., 1:2]) - - # quad = quadrature.ndiagquad( - # fun, num_gauss_hermite_points, [mu1, mu2], [var1, var2], logspace=True, - # ) - quad = quadrature.ndiagquad(fun, num_gauss_hermite_points, mu, var, logspace=True,) - - # expected = mu1 + var1 / 2 + alpha * mu2 + alpha ** 2 * var2 / 2 - expected = ( - mu[..., 0:1] + var[..., 0:1] / 2 + alpha * mu[..., 1:2] + alpha ** 2 * var[..., 1:2] / 2 + quad = quadrature.ndiagquad( + lambda *X: (X[0] + alpha * X[1]), + num_gauss_hermite_points, + [mu1, mu2], + [var1, var2], + logspace=True, ) - + expected = mu1 + var1 / 2 + alpha * mu2 + alpha ** 2 * var2 / 2 assert_allclose(quad, expected) @@ -93,13 +67,8 @@ def test_diagquad_logspace(mu1, var1, mu2, var2): def test_diagquad_with_kwarg(mu1, var1): alpha = np.array([2.5, -1.3]) num_gauss_hermite_points = 25 - - mu1 = mu1[:, None] - var1 = var1[:, None] - alpha = alpha[:, None] - - fun = lambda X, Y: tf.exp(X * Y) - quad = quadrature.ndiagquad(fun, num_gauss_hermite_points, mu1, var1, Y=alpha) - + quad = quadrature.ndiagquad( + lambda X, Y: tf.exp(X * Y), num_gauss_hermite_points, mu1, var1, Y=alpha + ) expected = np.exp(alpha * mu1 + alpha ** 2 * var1 / 2) assert_allclose(quad, expected) From 48516eafec7f74a90bb0d8b1861bc9d1a939b225 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 9 Jun 2020 13:55:07 -0300 Subject: [PATCH 18/41] renaming --- gpflow/quadrature/__init__.py | 6 +- gpflow/quadrature/deprecated.py | 201 ++++++++++++++++++++++++++------ gpflow/quadrature/old.py | 73 ++++++++++++ 3 files changed, 243 insertions(+), 37 deletions(-) create mode 100644 gpflow/quadrature/old.py diff --git a/gpflow/quadrature/__init__.py b/gpflow/quadrature/__init__.py index 663f98c7e..1f0c68057 100644 --- a/gpflow/quadrature/__init__.py +++ b/gpflow/quadrature/__init__.py @@ -1,5 +1,5 @@ from .base import GaussianQuadrature -from .gauss_hermite import NDDiagGHQuadrature -from .quadrature import hermgauss, mvhermgauss, mvnquad, ndiagquad, ndiag_mc +from .gauss_hermite import NDiagGHQuadrature +from .deprecated import hermgauss, mvhermgauss, mvnquad, ndiagquad, ndiag_mc -from .deprecated import ndiagquad as ndiagquad_old +from .old import ndiagquad as ndiagquad_old diff --git a/gpflow/quadrature/deprecated.py b/gpflow/quadrature/deprecated.py index 45c06a0e0..ae043954d 100644 --- a/gpflow/quadrature/deprecated.py +++ b/gpflow/quadrature/deprecated.py @@ -1,17 +1,112 @@ +# Copyright 2017-2018 the GPflow authors. +# +# 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 +from collections.abc import Iterable + import numpy as np import tensorflow as tf -from collections.abc import Iterable +from ..config import default_float +from ..utilities import to_default_float + +from .gauss_hermite import NDiagGHQuadrature + + +def hermgauss(n: int): + x, w = np.polynomial.hermite.hermgauss(n) + x, w = x.astype(default_float()), w.astype(default_float()) + return x, w + + +def mvhermgauss(H: int, D: int): + """ + 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 = 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 -from .quadrature import mvhermgauss + +def mvnquad(func, means, covs, H: int, Din: int = None, Dout=None): + """ + 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 means: NxD + :param covs: NxDxD + :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 (N,*Dout) + """ + # 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) def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): """ 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. + :param funcs: the integrand(s): Callable or Iterable of Callables that operates elementwise :param H: number of Gauss-Hermite quadrature points @@ -20,54 +115,92 @@ def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): :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 + Fmu, Fvar, Ys should all have same shape, with overall size `N` :return: shape is the same as that of the first Fmu """ + n_gh = H if isinstance(Fmu, (tuple, list)): - Din = len(Fmu) + 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) + if tf.rank(Fmu) < 2: + Fmu = tf.expand_dims(Fmu, axis=-1) + Fvar = tf.expand_dims(Fvar, axis=-1) - def unify(f_list): - """Stack a list of means/vars into a full block.""" - return tf.reshape( - tensor=tf.concat([tf.reshape(f, shape=(-1, 1)) for f in f_list], axis=1), - shape=(-1, 1, Din), - ) + def wrapper(old_fun): + def new_fun(X, **Ys): + Xs = tf.unstack(X, axis=-1) + fun_eval = old_fun(*Xs, **Ys) + if tf.rank(fun_eval) < tf.rank(X): + fun_eval = tf.expand_dims(fun_eval, axis=-1) + return fun_eval + return new_fun - shape = tf.shape(Fmu[0]) - Fmu, Fvar = map(unify, [Fmu, Fvar]) # both [N, 1, Din] + if isinstance(funcs, Iterable): + funcs = [wrapper(f) for f in funcs] + else: + funcs = wrapper(funcs) - print(Fmu) - print(Fvar) + quadrature = NDiagGHQuadrature(dim, n_gh) + if logspace: + result = quadrature.logspace(funcs, Fmu, Fvar, **Ys) else: - Din = 1 - shape = tf.shape(Fmu) - Fmu, Fvar = [tf.reshape(f, (-1, 1, 1)) for f in [Fmu, Fvar]] + result = quadrature(funcs, Fmu, Fvar, **Ys) - xn, wn = mvhermgauss(H, Din) - # xn: H**Din x Din, wn: H**Din + if isinstance(result, list): + result = [tf.reshape(r, shape) for r in result] + else: + result = tf.reshape(result, shape) + + return result - gh_x = xn.reshape(1, -1, Din) # [1, H]**Din x Din - Xall = gh_x * tf.sqrt(2.0 * Fvar) + Fmu # [N, H]**Din x Din - Xs = [Xall[:, :, i] for i in range(Din)] # [N, H]**Din each - gh_w = wn * np.pi ** (-0.5 * Din) # H**Din x 1 +def ndiag_mc(funcs, S: int, Fmu, Fvar, logspace: bool = False, epsilon=None, **Ys): + """ + Computes N Gaussian expectation integrals of one or more functions + using Monte Carlo samples. The Gaussians must be independent. + + :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 + + Fmu, Fvar, Ys should all have same shape, with overall size `N` + :return: shape is the same as that of the first Fmu + """ + N, D = Fmu.shape[0], Fvar.shape[1] + + if epsilon is None: + epsilon = tf.random.normal((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(): - Y = tf.reshape(Y, (-1, 1)) - Y = tf.tile(Y, [1, H ** Din]) # broadcast Y to match X - # without the tiling, some calls such as tf.where() (in bernoulli) fail - Ys[name] = Y # now [N, H]**Din + 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(f): - feval = f(*Xs, **Ys) # f should be elementwise: return shape [N, H]**Din + def eval_func(func): + feval = func(mc_Xr, **Ys) + feval = tf.reshape(feval, (S, N, -1)) if logspace: - log_gh_w = np.log(gh_w.reshape(1, -1)) - result = tf.reduce_logsumexp(feval + log_gh_w, axis=1) + log_S = tf.math.log(to_default_float(S)) + return tf.reduce_logsumexp(feval, axis=0) - log_S # [N, D] else: - result = tf.linalg.matmul(feval, gh_w.reshape(-1, 1)) - return tf.reshape(result, shape) + return tf.reduce_mean(feval, axis=0) if isinstance(funcs, Iterable): return [eval_func(f) for f in funcs] - - return eval_func(funcs) + else: + return eval_func(funcs) diff --git a/gpflow/quadrature/old.py b/gpflow/quadrature/old.py new file mode 100644 index 000000000..9a8cec164 --- /dev/null +++ b/gpflow/quadrature/old.py @@ -0,0 +1,73 @@ +import numpy as np +import tensorflow as tf + +from collections.abc import Iterable + +from .deprecated import mvhermgauss + + +def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): + """ + 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. + :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 + Fmu, Fvar, Ys should all have same shape, with overall size `N` + :return: shape is the same as that of the first Fmu + """ + if isinstance(Fmu, (tuple, list)): + Din = len(Fmu) + + def unify(f_list): + """Stack a list of means/vars into a full block.""" + return tf.reshape( + tensor=tf.concat([tf.reshape(f, shape=(-1, 1)) for f in f_list], axis=1), + shape=(-1, 1, Din), + ) + + shape = tf.shape(Fmu[0]) + Fmu, Fvar = map(unify, [Fmu, Fvar]) # both [N, 1, Din] + + print(Fmu) + print(Fvar) + else: + Din = 1 + shape = tf.shape(Fmu) + Fmu, Fvar = [tf.reshape(f, (-1, 1, 1)) for f in [Fmu, Fvar]] + + xn, wn = mvhermgauss(H, Din) + # xn: H**Din x Din, wn: H**Din + + gh_x = xn.reshape(1, -1, Din) # [1, H]**Din x Din + Xall = gh_x * tf.sqrt(2.0 * Fvar) + Fmu # [N, H]**Din x Din + Xs = [Xall[:, :, i] for i in range(Din)] # [N, H]**Din each + + gh_w = wn * np.pi ** (-0.5 * Din) # H**Din x 1 + + for name, Y in Ys.items(): + Y = tf.reshape(Y, (-1, 1)) + Y = tf.tile(Y, [1, H ** Din]) # broadcast Y to match X + # without the tiling, some calls such as tf.where() (in bernoulli) fail + Ys[name] = Y # now [N, H]**Din + + def eval_func(f): + feval = f(*Xs, **Ys) # f should be elementwise: return shape [N, H]**Din + if logspace: + log_gh_w = np.log(gh_w.reshape(1, -1)) + result = tf.reduce_logsumexp(feval + log_gh_w, axis=1) + else: + result = tf.linalg.matmul(feval, gh_w.reshape(-1, 1)) + return tf.reshape(result, shape) + + if isinstance(funcs, Iterable): + return [eval_func(f) for f in funcs] + + return eval_func(funcs) From 6ab2ede1fe840b1f88a8113dd62d4b861bcd0411 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 9 Jun 2020 13:55:45 -0300 Subject: [PATCH 19/41] renaming --- gpflow/quadrature/gauss_hermite.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index aa349060b..78cb08302 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -47,7 +47,7 @@ def ndgh_points_and_weights(dim: int, n_gh: int): return reshape_Z_dZ(zs, dzs) -class NDDiagGHQuadrature(GaussianQuadrature): +class NDiagGHQuadrature(GaussianQuadrature): def __init__(self, dim: int, n_gh: int): Z, dZ = ndgh_points_and_weights(dim, n_gh) self.n_gh_total = n_gh ** dim From 592bfdcdef20ce119db16f57d1dea3cb2de83f87 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 9 Jun 2020 13:55:51 -0300 Subject: [PATCH 20/41] renaming --- gpflow/quadrature/quadrature.py | 183 -------------------------------- 1 file changed, 183 deletions(-) delete mode 100644 gpflow/quadrature/quadrature.py diff --git a/gpflow/quadrature/quadrature.py b/gpflow/quadrature/quadrature.py deleted file mode 100644 index 016d843f5..000000000 --- a/gpflow/quadrature/quadrature.py +++ /dev/null @@ -1,183 +0,0 @@ -# Copyright 2017-2018 the GPflow authors. -# -# 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 -from collections.abc import Iterable - -import numpy as np -import tensorflow as tf - -from ..config import default_float -from ..utilities import to_default_float - -from .gauss_hermite import NDDiagGHQuadrature - - -def hermgauss(n: int): - x, w = np.polynomial.hermite.hermgauss(n) - x, w = x.astype(default_float()), w.astype(default_float()) - return x, w - - -def mvhermgauss(H: int, D: int): - """ - 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 = 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 - - -def mvnquad(func, means, covs, H: int, Din: int = None, Dout=None): - """ - 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 means: NxD - :param covs: NxDxD - :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 (N,*Dout) - """ - # 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) - - -def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): - """ - 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. - - :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 - - Fmu, Fvar, Ys should all have same shape, with overall size `N` - :return: shape is the same as that of the first Fmu - """ - n_gh = H - if isinstance(Fmu, (tuple, list)): - dim = len(Fmu) - Fmu = tf.stack(Fmu, axis=-1) - Fvar = tf.stack(Fvar, axis=-1) - else: - dim = tf.shape(Fmu)[-1] - - quadrature = NDDiagGHQuadrature(dim, n_gh) - - if logspace: - result = quadrature.logspace(funcs, Fmu, Fvar, **Ys) - else: - result = quadrature(funcs, Fmu, Fvar, **Ys) - - return result - - -def ndiag_mc(funcs, S: int, Fmu, Fvar, logspace: bool = False, epsilon=None, **Ys): - """ - Computes N Gaussian expectation integrals of one or more functions - using Monte Carlo samples. The Gaussians must be independent. - - :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 - - Fmu, Fvar, Ys should all have same shape, with overall size `N` - :return: shape is the same as that of the first Fmu - """ - N, D = Fmu.shape[0], Fvar.shape[1] - - if epsilon is None: - epsilon = tf.random.normal((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): - 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) From 89756c215a2ee2109bd0b3bafacf47e5893f7d20 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 9 Jun 2020 15:10:35 -0300 Subject: [PATCH 21/41] reshape kwargs --- gpflow/quadrature/deprecated.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/gpflow/quadrature/deprecated.py b/gpflow/quadrature/deprecated.py index ae043954d..dfeb7f854 100644 --- a/gpflow/quadrature/deprecated.py +++ b/gpflow/quadrature/deprecated.py @@ -128,9 +128,11 @@ def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): else: dim = 1 shape = tf.shape(Fmu) - if tf.rank(Fmu) < 2: - Fmu = tf.expand_dims(Fmu, axis=-1) - Fvar = tf.expand_dims(Fvar, axis=-1) + + 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): def new_fun(X, **Ys): From 0b76b7d8b879b79a9fa4db2648b9a2b7a69ecfb1 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 9 Jun 2020 15:11:44 -0300 Subject: [PATCH 22/41] quadrature along axis=-2, simplified broadcasting --- gpflow/quadrature/base.py | 8 ++++---- gpflow/quadrature/gauss_hermite.py | 25 +++++++++++-------------- 2 files changed, 15 insertions(+), 18 deletions(-) diff --git a/gpflow/quadrature/base.py b/gpflow/quadrature/base.py index f6ffadd4c..51c538c8a 100644 --- a/gpflow/quadrature/base.py +++ b/gpflow/quadrature/base.py @@ -13,8 +13,8 @@ def logspace(self, fun, mean, var, *args, **kwargs): X, W = self._build_X_W(mean, var) logW = tf.math.log(W) if isinstance(fun, Iterable): - return [tf.reduce_logsumexp(f(X, *args, **kwargs) + logW, axis=0) for f in fun] - return tf.reduce_logsumexp(fun(X, *args, **kwargs) + logW, axis=0) + return [tf.reduce_logsumexp(f(X, *args, **kwargs) + logW, axis=-2) for f in fun] + return tf.reduce_logsumexp(fun(X, *args, **kwargs) + logW, axis=-2) def __call__(self, fun, mean, var, *args, **kwargs): X, W = self._build_X_W(mean, var) @@ -23,5 +23,5 @@ def __call__(self, fun, mean, var, *args, **kwargs): # and sum-reducing all at once # The problem: there is no garantee that f1(X), f2(X), ... # have comaptible shapes - return [tf.reduce_sum(f(X, *args, **kwargs) * W, axis=0) for f in fun] - return tf.reduce_sum(fun(X, *args, **kwargs) * W, axis=0) + return [tf.reduce_sum(f(X, *args, **kwargs) * W, axis=-2) for f in fun] + return tf.reduce_sum(fun(X, *args, **kwargs) * W, axis=-2) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index 78cb08302..5f6a667ad 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -53,21 +53,18 @@ def __init__(self, dim: int, n_gh: int): self.n_gh_total = n_gh ** dim self.Z = tf.convert_to_tensor(Z) self.dZ = tf.convert_to_tensor(dZ) + # Z: [n_gh_total, dims] + # dZ: [n_gh_total, 1] def _build_X_W(self, mean, var): - new_shape = (self.n_gh_total,) + tuple([1] * (len(mean.shape) - 1)) + (-1,) - Z = tf.reshape(self.Z, new_shape) - dZ = tf.reshape(self.dZ, new_shape) - # Z : [n_gh_total] + [1]*len(batch) + [dims], typically [n_gh_total, 1, dims] - # dZ: [n_gh_total] + [1]*len(batch) + [1], typically [n_gh_total, 1, 1] - - mean = tf.expand_dims(mean, 0) - stddev = tf.expand_dims(tf.sqrt(var), 0) - # mean, stddev: [1] + batch + [dims], typically [1, N, dims] - - X = mean + stddev * Z - W = dZ - # X: [n_gh_total] + batch + [dims], typically [n_gh_total, N, dims] - # W: [n_gh_total] + [1]*len(batch) + [1], typically [n_gh_total, 1, 1] + # mean, stddev: batch + [dims], typically [N, dims] + mean = tf.expand_dims(mean, -2) + stddev = tf.expand_dims(tf.sqrt(var), -2) + # mean, stddev: batch + [1, dims], typically [N, 1, dims] + + X = mean + stddev * self.Z + W = self.dZ + # X: batch + [n_gh_total, dims], typically [N, n_gh_total, dims] + # W: batch + [n_gh_total, 1], typically [N, n_gh_total, 1] return X, W From da8b7d281afc99cc5b04fd2e020f670e8ef1f11a Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 9 Jun 2020 15:12:26 -0300 Subject: [PATCH 23/41] black --- gpflow/quadrature/deprecated.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gpflow/quadrature/deprecated.py b/gpflow/quadrature/deprecated.py index dfeb7f854..c3f26b9ee 100644 --- a/gpflow/quadrature/deprecated.py +++ b/gpflow/quadrature/deprecated.py @@ -132,7 +132,7 @@ def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): 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()} + Ys = {Yname: tf.reshape(Y, (-1, 1)) for Yname, Y in Ys.items()} def wrapper(old_fun): def new_fun(X, **Ys): @@ -141,6 +141,7 @@ def new_fun(X, **Ys): if tf.rank(fun_eval) < tf.rank(X): fun_eval = tf.expand_dims(fun_eval, axis=-1) return fun_eval + return new_fun if isinstance(funcs, Iterable): From 41eac9b2a11f2d1b91ecf2e8c30444984747abec Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 9 Jun 2020 15:26:06 -0300 Subject: [PATCH 24/41] docs --- gpflow/quadrature/base.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/gpflow/quadrature/base.py b/gpflow/quadrature/base.py index 51c538c8a..8bda7850b 100644 --- a/gpflow/quadrature/base.py +++ b/gpflow/quadrature/base.py @@ -10,6 +10,14 @@ def _build_X_W(self, mean, var): raise NotImplementedError def logspace(self, fun, mean, var, *args, **kwargs): + r""" + Compute the Gaussian log-Expectation of a function f + q(x) = N(mean, var) + E_{X~q(x)}[f(X)] = log∫f(x)q(x)dx + Using the approximation + log \sum exp(f(x_i) + log w_i) + Where x_i, w_i must be provided by some inheriting class + """ X, W = self._build_X_W(mean, var) logW = tf.math.log(W) if isinstance(fun, Iterable): @@ -17,6 +25,14 @@ def logspace(self, fun, mean, var, *args, **kwargs): return tf.reduce_logsumexp(fun(X, *args, **kwargs) + logW, axis=-2) def __call__(self, fun, mean, var, *args, **kwargs): + r""" + Compute the Gaussian Expectation of a function f + q(x) = N(mean, var) + E_{X~q(x)}[f(X)] = ∫f(x)q(x)dx + Using the approximation + \sum f(x_i)*w_i + Where x_i, w_i must be provided by some inheriting class + """ X, W = self._build_X_W(mean, var) if isinstance(fun, Iterable): # Maybe this can be better implemented by stacking [f1(X), f2(X), ...] From 77fcfb170e1ab8e601f9db24410982740e0f6a01 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 9 Jun 2020 15:29:32 -0300 Subject: [PATCH 25/41] docs --- gpflow/quadrature/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gpflow/quadrature/base.py b/gpflow/quadrature/base.py index 8bda7850b..a29140530 100644 --- a/gpflow/quadrature/base.py +++ b/gpflow/quadrature/base.py @@ -13,7 +13,7 @@ def logspace(self, fun, mean, var, *args, **kwargs): r""" Compute the Gaussian log-Expectation of a function f q(x) = N(mean, var) - E_{X~q(x)}[f(X)] = log∫f(x)q(x)dx + E_{X~q(x)}[f(X)] = log ∫f(x)q(x)dx Using the approximation log \sum exp(f(x_i) + log w_i) Where x_i, w_i must be provided by some inheriting class From 1b5c22bcdce8091d9bc302ec522f8213e7189e79 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 9 Jun 2020 17:52:10 -0300 Subject: [PATCH 26/41] helper function --- gpflow/quadrature/gauss_hermite.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index 5f6a667ad..198481006 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -35,6 +35,10 @@ def reshape_Z_dZ(zs, dzs): return Z, dZ +def repeat_as_list(x, n): + return tf.unstack(tf.repeat(tf.expand_dims(x, axis=0), n, axis=0), axis=0) + + def ndgh_points_and_weights(dim: int, n_gh: int): r""" :param n_gh: number of Gauss-Hermite points, integer @@ -42,8 +46,8 @@ def ndgh_points_and_weights(dim: int, n_gh: int): :returns: points Z, with shape [n_gh**dim, dim], and weights dZ, with shape [n_gh**dim, 1] """ z, dz = gh_points_and_weights(n_gh) - zs = tf.unstack(tf.repeat(tf.expand_dims(z, axis=0), dim, axis=0), axis=0) - dzs = tf.unstack(tf.repeat(tf.expand_dims(dz, axis=0), dim, axis=0), axis=0) + zs = repeat_as_list(z, dim) + dzs = repeat_as_list(dz, dim) return reshape_Z_dZ(zs, dzs) From ed49cf04cde9d313afe8844a44297286b1814845 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Thu, 11 Jun 2020 09:24:54 -0300 Subject: [PATCH 27/41] docstrings and typing --- gpflow/quadrature/base.py | 98 ++++++++++++++++++++++-------- gpflow/quadrature/gauss_hermite.py | 70 +++++++++++++++------ gpflow/quadrature/old.py | 4 +- 3 files changed, 129 insertions(+), 43 deletions(-) diff --git a/gpflow/quadrature/base.py b/gpflow/quadrature/base.py index a29140530..5abf65482 100644 --- a/gpflow/quadrature/base.py +++ b/gpflow/quadrature/base.py @@ -1,38 +1,52 @@ -from collections.abc import Iterable -import abc +from typing import Callable, Union, Iterable +import abc import tensorflow as tf +from ..base import TensorType + class GaussianQuadrature: + """ + Abstract class implementing quadrature methods to compute Gaussian Expectations. + Inhering classes must provide the method _build_X_W to create points and weights + to be used for quadrature. + """ + @abc.abstractmethod - def _build_X_W(self, mean, var): + def _build_X_W(self, mean: TensorType, var: TensorType): raise NotImplementedError - def logspace(self, fun, mean, var, *args, **kwargs): - r""" - Compute the Gaussian log-Expectation of a function f - q(x) = N(mean, var) - E_{X~q(x)}[f(X)] = log ∫f(x)q(x)dx - Using the approximation - log \sum exp(f(x_i) + log w_i) - Where x_i, w_i must be provided by some inheriting class - """ - X, W = self._build_X_W(mean, var) - logW = tf.math.log(W) - if isinstance(fun, Iterable): - return [tf.reduce_logsumexp(f(X, *args, **kwargs) + logW, axis=-2) for f in fun] - return tf.reduce_logsumexp(fun(X, *args, **kwargs) + logW, axis=-2) - def __call__(self, fun, mean, var, *args, **kwargs): r""" - Compute the Gaussian Expectation of a function f - q(x) = N(mean, var) - E_{X~q(x)}[f(X)] = ∫f(x)q(x)dx - Using the approximation - \sum f(x_i)*w_i - Where x_i, w_i must be provided by some inheriting class + Compute the Gaussian Expectation of a function f: + + X ~ N(mean, var) + E[f(X)] = ∫f(x, *args, **kwargs)p(x)dx + + Using the formula: + E[f(X)] = sum_{i=1}^{N_quad_points} f(x_i) * w_i + + where x_i, w_i must be provided by the inheriting class through self._build_X_W. + The computations broadcast along batch-dimensions, represented by [b1, b2, ..., bX]. + + :param fun: Callable or Iterable of Callables that operates elementwise, with + signature f(X, *args, **kwargs). Moreover, if must satisfy the shape-mapping: + X shape: [b1, b2, ..., bX, N_quad_points, d], + usually [N, N_quad_points, d] + f(X) shape: [b1, b2, ...., bf, N_quad_points, d'], + usually [N, N_quad_points, 1] or [N, N_quad_points, d] + In most cases, f should only operate over the last dimension of X + :param mean: Array/Tensor with shape [b1, b2, ..., bX, d], usually [N, d], + representing the mean of a d-Variate Gaussian distribution + :param var: Array/Tensor with shape b1, b2, ..., bX, d], usually [N, d], + representing the variance of a d-Variate Gaussian distribution + :param *args: Passed to fun + :param **kargs: Passed to fun + :return: Array/Tensor with shape [b1, b2, ...., bf, N_quad_points, d'], + usually [N, d] or [N, 1] """ + X, W = self._build_X_W(mean, var) if isinstance(fun, Iterable): # Maybe this can be better implemented by stacking [f1(X), f2(X), ...] @@ -41,3 +55,39 @@ def __call__(self, fun, mean, var, *args, **kwargs): # have comaptible shapes return [tf.reduce_sum(f(X, *args, **kwargs) * W, axis=-2) for f in fun] return tf.reduce_sum(fun(X, *args, **kwargs) * W, axis=-2) + + def logspace(self, fun: Union[Callable, Iterable[Callable]], mean, var, *args, **kwargs): + r""" + Compute the Gaussian log-Expectation of a the exponential of a function f: + + X ~ N(mean, var) + log E[exp[f(X)]] = log ∫exp[f(x, *args, **kwargs)]p(x)dx + + Using the formula: + log E[exp[f(X)]] = log sum_{i=1}^{N_quad_points} exp[f(x_i) + log w_i] + + where x_i, w_i must be provided by the inheriting class through self._build_X_W. + The computations broadcast along batch-dimensions, represented by [b1, b2, ..., bX]. + + :param fun: Callable or Iterable of Callables that operates elementwise, with + signature f(X, *args, **kwargs). Moreover, if must satisfy the shape-mapping: + X shape: [b1, b2, ..., bX, N_quad_points, d], + usually [N, N_quad_points, d] + f(X) shape: [b1, b2, ...., bf, N_quad_points, d'], + usually [N, N_quad_points, 1] or [N, N_quad_points, d] + In most cases, f should only operate over the last dimension of X + :param mean: Array/Tensor with shape [b1, b2, ..., bX, d], usually [N, d], + representing the mean of a d-Variate Gaussian distribution + :param var: Array/Tensor with shape b1, b2, ..., bX, d], usually [N, d], + representing the variance of a d-Variate Gaussian distribution + :param *args: Passed to fun + :param **kargs: Passed to fun + :return: Array/Tensor with shape [b1, b2, ...., bf, N_quad_points, d'], + usually [N, d] or [N, 1] + """ + + X, W = self._build_X_W(mean, var) + logW = tf.math.log(W) + if isinstance(fun, Iterable): + return [tf.reduce_logsumexp(f(X, *args, **kwargs) + logW, axis=-2) for f in fun] + return tf.reduce_logsumexp(fun(X, *args, **kwargs) + logW, axis=-2) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index 198481006..e17f3b2d9 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -1,9 +1,13 @@ +from typing import List + import numpy as np import tensorflow as tf from .base import GaussianQuadrature from ..config import default_float +from ..base import TensorType + def gh_points_and_weights(n_gh: int): r""" @@ -11,12 +15,12 @@ def gh_points_and_weights(n_gh: int): returns the points z and the weights dz to perform the following uni-dimensional gaussian quadrature: - q(x) = N(mu, sigma²) - - E_{X~q(x)}[f(X)] = ∫ f(x)q(x)dx = \sum f(mu + sigma*z_k)*dz_k + X ~ N(mean, stddev²) + E[f(X)] = ∫f(x)p(x)dx = sum_{i=1}^{n_gh} f(mean + stddev*z_i)*dz_i - :param n_gh: number of Gauss-Hermite points, integer - :returns: points z and weights dz to compute uni-dimensional gaussian expectation, tuple + :param n_gh: Number of Gauss-Hermite points, integer + :returns: Points z and weights dz, both tensors with shape [n_gh], + to compute uni-dimensional gaussian expectation """ z, dz = np.polynomial.hermite.hermgauss(n_gh) z = z * np.sqrt(2) @@ -25,17 +29,33 @@ def gh_points_and_weights(n_gh: int): return tf.convert_to_tensor(z), tf.convert_to_tensor(dz) -def list_to_flat_grid(xs): +def list_to_flat_grid(xs: List[TensorType]): + """ + :param xs: List with d rank-1 Tensors, with shapes N1, N2, ..., Nd + :return: Tensor with shape [N1*N2*...*Nd, dim] representing the flattened + D-dimensional grid built from the input tensors xs + """ return tf.reshape(tf.stack(tf.meshgrid(*xs), axis=-1), (-1, len(xs))) -def reshape_Z_dZ(zs, dzs): +def reshape_Z_dZ(zs: List[TensorType], dzs: List[TensorType]): + """ + :param zs: List with d rank-1 Tensors, with shapes N1, N2, ..., Nd + :param dzs: List with d rank-1 Tensors, with shapes N1, N2, ..., Nd + :returns: points Z, Tensor with shape [n_gh**dim, dim], + and weights dZ, Tensor with shape [n_gh**dim, 1] + """ Z = list_to_flat_grid(zs) dZ = tf.reduce_prod(list_to_flat_grid(dzs), axis=-1, keepdims=True) return Z, dZ -def repeat_as_list(x, n): +def repeat_as_list(x: TensorType, n: int): + """ + :param x: Array/Tensor to be repeated + :param n: Integer with the number of repetitions + :return: List of n repetitions of Tensor x + """ return tf.unstack(tf.repeat(tf.expand_dims(x, axis=0), n, axis=0), axis=0) @@ -43,7 +63,8 @@ def ndgh_points_and_weights(dim: int, n_gh: int): r""" :param n_gh: number of Gauss-Hermite points, integer :param dim: dimension of the multivariate normal, integer - :returns: points Z, with shape [n_gh**dim, dim], and weights dZ, with shape [n_gh**dim, 1] + :returns: points Z, Tensor with shape [n_gh**dim, dim], + and weights dZ, Tensor with shape [n_gh**dim, 1] """ z, dz = gh_points_and_weights(n_gh) zs = repeat_as_list(z, dim) @@ -53,22 +74,37 @@ def ndgh_points_and_weights(dim: int, n_gh: int): class NDiagGHQuadrature(GaussianQuadrature): def __init__(self, dim: int, n_gh: int): + """ + :param n_gh: number of Gauss-Hermite points, integer + :param dim: dimension of the multivariate normal, integer + """ Z, dZ = ndgh_points_and_weights(dim, n_gh) self.n_gh_total = n_gh ** dim self.Z = tf.convert_to_tensor(Z) self.dZ = tf.convert_to_tensor(dZ) - # Z: [n_gh_total, dims] - # dZ: [n_gh_total, 1] - - def _build_X_W(self, mean, var): - # mean, stddev: batch + [dims], typically [N, dims] + # Z: [n_gh_total, dim] + # dZ: [n_gh_total, 1] + + def _build_X_W(self, mean: TensorType, var: TensorType): + """ + :param mean: Array/Tensor with shape [b1, b2, ..., bX, dim], usually [N, dim], + representing the mean of a dim-Variate Gaussian distribution + :param var: Array/Tensor with shape b1, b2, ..., bX, dim], usually [N, dim], + representing the variance of a dim-Variate Gaussian distribution + :return: points X, Tensor with shape [b1, b2, ..., bX, n_gh_total, dim], + usually [N, n_gh_total, dim], + and weights W, a Tensor with shape [b1, b2, ..., bX, n_gh_total, 1], + usually [N, n_gh_total, 1] + """ + + # mean, stddev: [b1, b2, ..., bX, dim], usually [N, dim] mean = tf.expand_dims(mean, -2) stddev = tf.expand_dims(tf.sqrt(var), -2) - # mean, stddev: batch + [1, dims], typically [N, 1, dims] + # mean, stddev: [b1, b2, ..., bX, 1, dim], usually [N, 1, dim] X = mean + stddev * self.Z W = self.dZ - # X: batch + [n_gh_total, dims], typically [N, n_gh_total, dims] - # W: batch + [n_gh_total, 1], typically [N, n_gh_total, 1] + # X: [b1, b2, ..., bX, n_gh_total, dim], usually [N, n_gh_total, dim] + # W: [b1, b2, ..., bX, n_gh_total, 1], usually [N, n_gh_total, 1] return X, W diff --git a/gpflow/quadrature/old.py b/gpflow/quadrature/old.py index 9a8cec164..34b1c375a 100644 --- a/gpflow/quadrature/old.py +++ b/gpflow/quadrature/old.py @@ -1,8 +1,8 @@ +from collections.abc import Iterable + import numpy as np import tensorflow as tf -from collections.abc import Iterable - from .deprecated import mvhermgauss From 650991dc2a07ef8b7ec27d5ca26ca537dddbd13d Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Sat, 13 Jun 2020 12:44:58 -0300 Subject: [PATCH 28/41] added new and old quadrature equivalence tests --- gpflow/quadrature/__init__.py | 2 - .../gpflow/quadrature/ndiagquad_old.py | 2 +- .../{ => quadrature}/test_quadrature.py | 0 .../quadrature/test_quadrature_equivalence.py | 85 +++++++++++++++++++ 4 files changed, 86 insertions(+), 3 deletions(-) rename gpflow/quadrature/old.py => tests/gpflow/quadrature/ndiagquad_old.py (98%) rename tests/gpflow/{ => quadrature}/test_quadrature.py (100%) create mode 100644 tests/gpflow/quadrature/test_quadrature_equivalence.py diff --git a/gpflow/quadrature/__init__.py b/gpflow/quadrature/__init__.py index 1f0c68057..81ef818e4 100644 --- a/gpflow/quadrature/__init__.py +++ b/gpflow/quadrature/__init__.py @@ -1,5 +1,3 @@ from .base import GaussianQuadrature from .gauss_hermite import NDiagGHQuadrature from .deprecated import hermgauss, mvhermgauss, mvnquad, ndiagquad, ndiag_mc - -from .old import ndiagquad as ndiagquad_old diff --git a/gpflow/quadrature/old.py b/tests/gpflow/quadrature/ndiagquad_old.py similarity index 98% rename from gpflow/quadrature/old.py rename to tests/gpflow/quadrature/ndiagquad_old.py index 34b1c375a..17e532c71 100644 --- a/gpflow/quadrature/old.py +++ b/tests/gpflow/quadrature/ndiagquad_old.py @@ -3,7 +3,7 @@ import numpy as np import tensorflow as tf -from .deprecated import mvhermgauss +from gpflow.quadrature.deprecated import mvhermgauss def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): diff --git a/tests/gpflow/test_quadrature.py b/tests/gpflow/quadrature/test_quadrature.py similarity index 100% rename from tests/gpflow/test_quadrature.py rename to tests/gpflow/quadrature/test_quadrature.py diff --git a/tests/gpflow/quadrature/test_quadrature_equivalence.py b/tests/gpflow/quadrature/test_quadrature_equivalence.py new file mode 100644 index 000000000..9cfc34d8a --- /dev/null +++ b/tests/gpflow/quadrature/test_quadrature_equivalence.py @@ -0,0 +1,85 @@ +# Copyright 2018 the GPflow authors. +# +# 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 numpy as np +import pytest +import tensorflow as tf +from numpy.testing import assert_allclose + +from gpflow.quadrature import ndiagquad +from ndiagquad_old import ndiagquad as ndiagquad_old + + +@pytest.mark.parametrize("mu", [np.array([1.0, 1.3])]) +@pytest.mark.parametrize("var", [np.array([3.0, 3.5])]) +def test_diagquad_1d(mu, var): + num_gauss_hermite_points = 25 + quad = ndiagquad([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, [mu], [var]) + quad_old = ndiagquad([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, [mu], [var]) + assert_allclose(quad[0], quad_old[0]) + + +@pytest.mark.parametrize("mu1", [np.array([1.0, 1.3])]) +@pytest.mark.parametrize("var1", [np.array([3.0, 3.5])]) +@pytest.mark.parametrize("mu2", [np.array([-2.0, 0.3])]) +@pytest.mark.parametrize("var2", [np.array([4.0, 4.2])]) +def test_diagquad_2d(mu1, var1, mu2, var2): + alpha = 2.5 + # using logspace=True we can reduce this, see test_diagquad_logspace + num_gauss_hermite_points = 35 + quad = ndiagquad( + lambda *X: tf.exp(X[0] + alpha * X[1]), num_gauss_hermite_points, [mu1, mu2], [var1, var2], + ) + quad_old = ndiagquad_old( + lambda *X: tf.exp(X[0] + alpha * X[1]), num_gauss_hermite_points, [mu1, mu2], [var1, var2], + ) + assert_allclose(quad, quad_old) + + +@pytest.mark.parametrize("mu1", [np.array([1.0, 1.3])]) +@pytest.mark.parametrize("var1", [np.array([3.0, 3.5])]) +@pytest.mark.parametrize("mu2", [np.array([-2.0, 0.3])]) +@pytest.mark.parametrize("var2", [np.array([4.0, 4.2])]) +def test_diagquad_logspace(mu1, var1, mu2, var2): + alpha = 2.5 + num_gauss_hermite_points = 25 + quad = ndiagquad( + lambda *X: (X[0] + alpha * X[1]), + num_gauss_hermite_points, + [mu1, mu2], + [var1, var2], + logspace=True, + ) + quad_old = ndiagquad_old( + lambda *X: (X[0] + alpha * X[1]), + num_gauss_hermite_points, + [mu1, mu2], + [var1, var2], + logspace=True, + ) + assert_allclose(quad, quad_old) + + +@pytest.mark.parametrize("mu1", [np.array([1.0, 1.3])]) +@pytest.mark.parametrize("var1", [np.array([3.0, 3.5])]) +def test_diagquad_with_kwarg(mu1, var1): + alpha = np.array([2.5, -1.3]) + num_gauss_hermite_points = 25 + quad = ndiagquad( + lambda X, Y: tf.exp(X * Y), num_gauss_hermite_points, mu1, var1, Y=alpha + ) + quad_old = ndiagquad_old( + lambda X, Y: tf.exp(X * Y), num_gauss_hermite_points, mu1, var1, Y=alpha + ) + assert_allclose(quad, quad_old) From d537f67c07fe66d6c3a6039b9467b588abc12007 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Sat, 13 Jun 2020 12:48:04 -0300 Subject: [PATCH 29/41] black --- tests/gpflow/quadrature/test_quadrature_equivalence.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/gpflow/quadrature/test_quadrature_equivalence.py b/tests/gpflow/quadrature/test_quadrature_equivalence.py index 9cfc34d8a..ac3a143e1 100644 --- a/tests/gpflow/quadrature/test_quadrature_equivalence.py +++ b/tests/gpflow/quadrature/test_quadrature_equivalence.py @@ -76,9 +76,7 @@ def test_diagquad_logspace(mu1, var1, mu2, var2): def test_diagquad_with_kwarg(mu1, var1): alpha = np.array([2.5, -1.3]) num_gauss_hermite_points = 25 - quad = ndiagquad( - lambda X, Y: tf.exp(X * Y), num_gauss_hermite_points, mu1, var1, Y=alpha - ) + quad = ndiagquad(lambda X, Y: tf.exp(X * Y), num_gauss_hermite_points, mu1, var1, Y=alpha) quad_old = ndiagquad_old( lambda X, Y: tf.exp(X * Y), num_gauss_hermite_points, mu1, var1, Y=alpha ) From bc98286757178159760b6a8001da74040663a87d Mon Sep 17 00:00:00 2001 From: gustavocmv <47801305+gustavocmv@users.noreply.github.com> Date: Mon, 20 Jul 2020 15:35:12 -0300 Subject: [PATCH 30/41] Removing comments Co-authored-by: Vincent Dutordoir --- gpflow/quadrature/base.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/gpflow/quadrature/base.py b/gpflow/quadrature/base.py index 5abf65482..f61f5c438 100644 --- a/gpflow/quadrature/base.py +++ b/gpflow/quadrature/base.py @@ -49,10 +49,6 @@ def __call__(self, fun, mean, var, *args, **kwargs): X, W = self._build_X_W(mean, var) if isinstance(fun, Iterable): - # Maybe this can be better implemented by stacking [f1(X), f2(X), ...] - # and sum-reducing all at once - # The problem: there is no garantee that f1(X), f2(X), ... - # have comaptible shapes return [tf.reduce_sum(f(X, *args, **kwargs) * W, axis=-2) for f in fun] return tf.reduce_sum(fun(X, *args, **kwargs) * W, axis=-2) From 5621a9a2d22af3d0b77e24a916cb71c073ed325c Mon Sep 17 00:00:00 2001 From: gustavocmv <47801305+gustavocmv@users.noreply.github.com> Date: Mon, 20 Jul 2020 15:38:10 -0300 Subject: [PATCH 31/41] Typo Co-authored-by: Vincent Dutordoir --- gpflow/quadrature/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gpflow/quadrature/base.py b/gpflow/quadrature/base.py index f61f5c438..60846746f 100644 --- a/gpflow/quadrature/base.py +++ b/gpflow/quadrature/base.py @@ -9,7 +9,7 @@ class GaussianQuadrature: """ Abstract class implementing quadrature methods to compute Gaussian Expectations. - Inhering classes must provide the method _build_X_W to create points and weights + Inheriting classes must provide the method _build_X_W to create points and weights to be used for quadrature. """ From 898f4c04970229b3c9ca85717c516a579997bbc8 Mon Sep 17 00:00:00 2001 From: gustavocmv <47801305+gustavocmv@users.noreply.github.com> Date: Mon, 20 Jul 2020 15:41:12 -0300 Subject: [PATCH 32/41] notation Co-authored-by: Vincent Dutordoir --- gpflow/quadrature/gauss_hermite.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index e17f3b2d9..b93d7ae51 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -32,8 +32,8 @@ def gh_points_and_weights(n_gh: int): def list_to_flat_grid(xs: List[TensorType]): """ :param xs: List with d rank-1 Tensors, with shapes N1, N2, ..., Nd - :return: Tensor with shape [N1*N2*...*Nd, dim] representing the flattened - D-dimensional grid built from the input tensors xs + :return: Tensor with shape [N1*N2*...*Nd, d] representing the flattened + d-dimensional grid built from the input tensors xs """ return tf.reshape(tf.stack(tf.meshgrid(*xs), axis=-1), (-1, len(xs))) From a69acdadba64ce47c713e0b65a01c2ed5c0b2138 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Mon, 20 Jul 2020 16:02:49 -0300 Subject: [PATCH 33/41] reshape_Z_dZ return docstring fix --- gpflow/quadrature/gauss_hermite.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index b93d7ae51..a5ea8d0f3 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -42,8 +42,8 @@ def reshape_Z_dZ(zs: List[TensorType], dzs: List[TensorType]): """ :param zs: List with d rank-1 Tensors, with shapes N1, N2, ..., Nd :param dzs: List with d rank-1 Tensors, with shapes N1, N2, ..., Nd - :returns: points Z, Tensor with shape [n_gh**dim, dim], - and weights dZ, Tensor with shape [n_gh**dim, 1] + :returns: points Z, Tensor with shape [N1*N2*...*Nd, d], + and weights dZ, Tensor with shape [N1*N2*...*Nd, 1] """ Z = list_to_flat_grid(zs) dZ = tf.reduce_prod(list_to_flat_grid(dzs), axis=-1, keepdims=True) From 187bf354c4c9754272fdc45d5e04ddca2677011c Mon Sep 17 00:00:00 2001 From: gustavocmv <47801305+gustavocmv@users.noreply.github.com> Date: Tue, 21 Jul 2020 11:31:37 -0300 Subject: [PATCH 34/41] FIX: quad_old computed with the ndiagquad_old Co-authored-by: Vincent Dutordoir --- tests/gpflow/quadrature/test_quadrature_equivalence.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/gpflow/quadrature/test_quadrature_equivalence.py b/tests/gpflow/quadrature/test_quadrature_equivalence.py index ac3a143e1..1d8697f65 100644 --- a/tests/gpflow/quadrature/test_quadrature_equivalence.py +++ b/tests/gpflow/quadrature/test_quadrature_equivalence.py @@ -26,7 +26,7 @@ def test_diagquad_1d(mu, var): num_gauss_hermite_points = 25 quad = ndiagquad([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, [mu], [var]) - quad_old = ndiagquad([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, [mu], [var]) + quad_old = ndiagquad_old([lambda *X: tf.exp(X[0])], num_gauss_hermite_points, [mu], [var]) assert_allclose(quad[0], quad_old[0]) From 8e138379d737ec6e5d2b9a0f337eb93440d988ad Mon Sep 17 00:00:00 2001 From: gustavocmv <47801305+gustavocmv@users.noreply.github.com> Date: Tue, 21 Jul 2020 12:43:48 -0300 Subject: [PATCH 35/41] more readable implementation Co-authored-by: Vincent Dutordoir --- gpflow/quadrature/gauss_hermite.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index a5ea8d0f3..231c5261a 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -56,7 +56,7 @@ def repeat_as_list(x: TensorType, n: int): :param n: Integer with the number of repetitions :return: List of n repetitions of Tensor x """ - return tf.unstack(tf.repeat(tf.expand_dims(x, axis=0), n, axis=0), axis=0) + return [x for _ in range(n)] def ndgh_points_and_weights(dim: int, n_gh: int): From 25245618fc6c9869370c96d39717a79db9c5dfa0 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 21 Jul 2020 13:31:53 -0300 Subject: [PATCH 36/41] tf.ensure_shape added --- gpflow/quadrature/gauss_hermite.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index 231c5261a..bb4c29465 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -78,12 +78,12 @@ def __init__(self, dim: int, n_gh: int): :param n_gh: number of Gauss-Hermite points, integer :param dim: dimension of the multivariate normal, integer """ - Z, dZ = ndgh_points_and_weights(dim, n_gh) + self.dim = dim + self.n_gh = n_gh self.n_gh_total = n_gh ** dim - self.Z = tf.convert_to_tensor(Z) - self.dZ = tf.convert_to_tensor(dZ) - # Z: [n_gh_total, dim] - # dZ: [n_gh_total, 1] + Z, dZ = ndgh_points_and_weights(self.dim, self.n_gh) + self.Z = tf.ensure_shape(Z, (self.n_gh_total, self.dim)) + self.dZ = tf.ensure_shape(dZ, (self.n_gh_total, self.dim)) def _build_X_W(self, mean: TensorType, var: TensorType): """ From 7bb0e9f1e0f2b0e225a2b8a5b3092c4c2f24ba91 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 21 Jul 2020 13:35:01 -0300 Subject: [PATCH 37/41] removed ndiagquad --- gpflow/quadrature/deprecated.py | 64 --------------------------------- 1 file changed, 64 deletions(-) diff --git a/gpflow/quadrature/deprecated.py b/gpflow/quadrature/deprecated.py index c3f26b9ee..ffc0a3b25 100644 --- a/gpflow/quadrature/deprecated.py +++ b/gpflow/quadrature/deprecated.py @@ -99,70 +99,6 @@ def mvnquad(func, means, covs, H: int, Din: int = None, Dout=None): return tf.reduce_sum(fX * wr, 0) -def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): - """ - 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. - - :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 - - Fmu, Fvar, Ys should all have same shape, with overall size `N` - :return: shape is the same as that of the first Fmu - """ - 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): - def new_fun(X, **Ys): - Xs = tf.unstack(X, axis=-1) - fun_eval = old_fun(*Xs, **Ys) - if tf.rank(fun_eval) < tf.rank(X): - fun_eval = tf.expand_dims(fun_eval, axis=-1) - return 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 - - def ndiag_mc(funcs, S: int, Fmu, Fvar, logspace: bool = False, epsilon=None, **Ys): """ Computes N Gaussian expectation integrals of one or more functions From 8e235241a697696e361158c30ff9aa9b4cc69f8a Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 21 Jul 2020 13:35:20 -0300 Subject: [PATCH 38/41] removed ndiagquad --- gpflow/quadrature/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gpflow/quadrature/__init__.py b/gpflow/quadrature/__init__.py index 81ef818e4..42b623f8c 100644 --- a/gpflow/quadrature/__init__.py +++ b/gpflow/quadrature/__init__.py @@ -1,3 +1,3 @@ from .base import GaussianQuadrature from .gauss_hermite import NDiagGHQuadrature -from .deprecated import hermgauss, mvhermgauss, mvnquad, ndiagquad, ndiag_mc +from .deprecated import hermgauss, mvhermgauss, mvnquad, ndiag_mc From d6e7e63d69a85896ea6b1077f7c7d3c4c30b5fc4 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 21 Jul 2020 13:48:04 -0300 Subject: [PATCH 39/41] Revert "removed ndiagquad" This reverts commit 7bb0e9f1e0f2b0e225a2b8a5b3092c4c2f24ba91. --- gpflow/quadrature/deprecated.py | 64 +++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/gpflow/quadrature/deprecated.py b/gpflow/quadrature/deprecated.py index ffc0a3b25..c3f26b9ee 100644 --- a/gpflow/quadrature/deprecated.py +++ b/gpflow/quadrature/deprecated.py @@ -99,6 +99,70 @@ def mvnquad(func, means, covs, H: int, Din: int = None, Dout=None): return tf.reduce_sum(fX * wr, 0) +def ndiagquad(funcs, H: int, Fmu, Fvar, logspace: bool = False, **Ys): + """ + 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. + + :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 + + Fmu, Fvar, Ys should all have same shape, with overall size `N` + :return: shape is the same as that of the first Fmu + """ + 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): + def new_fun(X, **Ys): + Xs = tf.unstack(X, axis=-1) + fun_eval = old_fun(*Xs, **Ys) + if tf.rank(fun_eval) < tf.rank(X): + fun_eval = tf.expand_dims(fun_eval, axis=-1) + return 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 + + def ndiag_mc(funcs, S: int, Fmu, Fvar, logspace: bool = False, epsilon=None, **Ys): """ Computes N Gaussian expectation integrals of one or more functions From 8f7e949a8e67680b1769eb685f5520c4d2c9fb88 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 21 Jul 2020 13:53:43 -0300 Subject: [PATCH 40/41] FIX: shape checking of dZ --- gpflow/quadrature/gauss_hermite.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gpflow/quadrature/gauss_hermite.py b/gpflow/quadrature/gauss_hermite.py index bb4c29465..86aab372b 100644 --- a/gpflow/quadrature/gauss_hermite.py +++ b/gpflow/quadrature/gauss_hermite.py @@ -83,7 +83,7 @@ def __init__(self, dim: int, n_gh: int): self.n_gh_total = n_gh ** dim Z, dZ = ndgh_points_and_weights(self.dim, self.n_gh) self.Z = tf.ensure_shape(Z, (self.n_gh_total, self.dim)) - self.dZ = tf.ensure_shape(dZ, (self.n_gh_total, self.dim)) + self.dZ = tf.ensure_shape(dZ, (self.n_gh_total, 1)) def _build_X_W(self, mean: TensorType, var: TensorType): """ From d4f313c6c6caa0cc80cd29ee1ebd017ae62e69f2 Mon Sep 17 00:00:00 2001 From: Gustavo Carvalho Date: Tue, 21 Jul 2020 14:08:50 -0300 Subject: [PATCH 41/41] Revert "removed ndiagquad" This reverts commit 8e235241a697696e361158c30ff9aa9b4cc69f8a. --- gpflow/quadrature/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gpflow/quadrature/__init__.py b/gpflow/quadrature/__init__.py index 42b623f8c..81ef818e4 100644 --- a/gpflow/quadrature/__init__.py +++ b/gpflow/quadrature/__init__.py @@ -1,3 +1,3 @@ from .base import GaussianQuadrature from .gauss_hermite import NDiagGHQuadrature -from .deprecated import hermgauss, mvhermgauss, mvnquad, ndiag_mc +from .deprecated import hermgauss, mvhermgauss, mvnquad, ndiagquad, ndiag_mc