# Copyright 2017 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.
from collections import namedtuple
import tensorflow as tf
import numpy as np
from numpy.testing import assert_almost_equal
import pytest
import gpflow
from gpflow import settings
from gpflow.conditionals import conditional
from gpflow.conditionals import uncertain_conditional
from gpflow.quadrature import mvnquad
from gpflow.test_util import session_context
class MomentMatchingSVGP(gpflow.models.SVGP):
@gpflow.params_as_tensors
def uncertain_predict_f_moment_matching(self, Xmu, Xcov):
return uncertain_conditional(
Xmu, Xcov, self.feature, self.kern, self.q_mu, self.q_sqrt,
mean_function=self.mean_function, white=self.whiten,
full_output_cov=self.full_output_cov)
def uncertain_predict_f_monte_carlo(self, Xmu, Xchol, mc_iter=int(1e6)):
rng = np.random.RandomState(0)
D_in = Xchol.shape[0]
X_samples = Xmu + np.reshape(
Xchol[None, :, :] @ rng.randn(mc_iter, D_in)[:, :, None], [mc_iter, D_in])
F_mu, F_var = self.predict_f(X_samples)
F_samples = F_mu + rng.randn(*F_var.shape) * (F_var ** 0.5)
mean = np.mean(F_samples, axis=0)
covar = np.cov(F_samples.T)
return mean, covar
def gen_L(rng, n, *shape):
return np.array([np.tril(rng.randn(*shape)) for _ in range(n)])
def gen_q_sqrt(rng, D_out, *shape):
return np.array([np.tril(rng.randn(*shape)) for _ in range(D_out)])
def mean_function_factory(rng, mean_function_name, D_in, D_out):
if mean_function_name == "Zero":
return gpflow.mean_functions.Zero(output_dim=D_out)
elif mean_function_name == "Constant":
return gpflow.mean_functions.Constant(c=rng.rand(D_out))
elif mean_function_name == "Linear":
return gpflow.mean_functions.Linear(
A=rng.rand(D_in, D_out), b=rng.rand(D_out))
else:
return None
class Data:
N = 7
N_new = 2
D_out = 3
D_in = 1
rng = np.random.RandomState(1)
X = np.linspace(-5, 5, N)[:, None] + rng.randn(N, 1)
Y = np.hstack([np.sin(X), np.cos(X), X**2])
Xnew_mu = rng.randn(N_new, 1)
Xnew_covar = np.zeros((N_new, 1, 1))
class DataMC1(Data):
Y = np.hstack([np.sin(Data.X), np.sin(Data.X) * 2, Data.X ** 2])
class DataMC2(Data):
N = 7
N_new = 5
D_out = 4
D_in = 2
X = Data.rng.randn(N, D_in)
Y = np.hstack([np.sin(X), np.sin(X)])
Xnew_mu = Data.rng.randn(N_new, D_in)
L = gen_L(Data.rng, N_new, D_in, D_in)
Xnew_covar = np.array([l @ l.T for l in L])
class DataQuadrature:
num_data = 10
num_ind = 10
D_in = 2
D_out = 3
H = 150
rng = np.random.RandomState(1)
Xmu = rng.randn(num_data, D_in)
L = gen_L(rng, num_data, D_in, D_in)
Xvar = np.array([l @ l.T for l in L])
Z = rng.randn(num_ind, D_in)
q_mu = rng.randn(num_ind, D_out)
q_sqrt = gen_q_sqrt(rng, D_out, num_ind, num_ind)
@classmethod
def tensors(cls, white, mean_name):
float_type = settings.float_type
Xmu = tf.placeholder(float_type, [cls.num_data, cls.D_in])
Xvar = tf.placeholder(float_type, [cls.num_data, cls.D_in, cls.D_in])
q_mu = tf.placeholder(float_type, [cls.num_ind, cls.D_out])
q_sqrt = tf.placeholder(float_type, [cls.D_out, cls.num_ind, cls.num_ind])
kern = gpflow.kernels.RBF(cls.D_in)
feat = gpflow.features.InducingPoints(cls.Z)
mean_function = mean_function_factory(cls.rng, mean_name, cls.D_in, cls.D_out)
effective_mean = mean_function or (lambda X: 0.0)
feed_dict = {
Xmu: cls.Xmu,
Xvar: cls.Xvar,
q_mu: cls.q_mu,
q_sqrt: cls.q_sqrt
}
def mean_fn(X):
mean, _ = conditional(X, feat, kern, q_mu, q_sqrt=q_sqrt, white=white)
return mean + effective_mean(X)
def var_fn(X):
_, var = conditional(X, feat, kern, q_mu, q_sqrt=q_sqrt, white=white)
return var
def mean_sq_fn(X):
mean, _ = conditional(X, feat, kern, q_mu, q_sqrt=q_sqrt, white=white)
return (mean + effective_mean(X)) ** 2
Collection = namedtuple('QuadratureCollection',
'Xmu,Xvar,q_mu,q_sqrt,'
'kern,feat,mean_function,'
'feed_dict,mean_fn,'
'var_fn,mean_sq_fn')
return Collection(Xmu=Xmu,
Xvar=Xvar,
q_mu=q_mu,
q_sqrt=q_sqrt,
kern=kern,
feat=feat,
mean_function=mean_function,
feed_dict=feed_dict,
mean_fn=mean_fn,
var_fn=var_fn,
mean_sq_fn=mean_sq_fn)
MEANS = ["Constant", "Linear", "Zero", None]
@pytest.mark.parametrize('white', [True, False])
@pytest.mark.parametrize('mean', MEANS)
def test_no_uncertainty(white, mean):
with session_context() as sess:
m = mean_function_factory(Data.rng, mean, Data.D_in, Data.D_out)
k = gpflow.kernels.RBF(1, variance=Data.rng.rand())
model = MomentMatchingSVGP(
Data.X, Data.Y, k, gpflow.likelihoods.Gaussian(),
mean_function=m, Z=Data.X.copy(), whiten=white)
model.full_output_cov = False
gpflow.train.AdamOptimizer().minimize(model, maxiter=50)
mean1, var1 = model.predict_f(Data.Xnew_mu)
pred_mm = model.uncertain_predict_f_moment_matching(
tf.constant(Data.Xnew_mu), tf.constant(Data.Xnew_covar))
mean2, var2 = sess.run(pred_mm)
assert_almost_equal(mean1, mean2)
for n in range(Data.N_new):
assert_almost_equal(var1[n, :], var2[n, ...])
@pytest.mark.parametrize('white', [True, False])
@pytest.mark.parametrize('mean', MEANS)
def test_monte_carlo_1_din(white, mean):
with session_context() as sess:
k = gpflow.kernels.RBF(1, variance=DataMC1.rng.rand())
m = mean_function_factory(DataMC1.rng, mean, DataMC1.D_in, DataMC1.D_out)
model = MomentMatchingSVGP(
DataMC1.X, DataMC1.Y, k, gpflow.likelihoods.Gaussian(),
Z=DataMC1.X.copy(), whiten=white)
model.full_output_cov = True
gpflow.train.AdamOptimizer().minimize(model, maxiter=50)
pred_mm = model.uncertain_predict_f_moment_matching(
tf.constant(DataMC1.Xnew_mu), tf.constant(DataMC1.Xnew_covar))
mean1, var1 = sess.run(pred_mm)
for n in range(DataMC1.N_new):
mean2, var2 = model.uncertain_predict_f_monte_carlo(
DataMC1.Xnew_mu[n, ...],
DataMC1.Xnew_covar[n, ...] ** 0.5)
assert_almost_equal(mean1[n, ...], mean2, decimal=3)
assert_almost_equal(var1[n, ...], var2, decimal=2)
@pytest.mark.parametrize('white', [True, False])
@pytest.mark.parametrize('mean', MEANS)
def test_monte_carlo_2_din(white, mean):
with session_context() as sess:
k = gpflow.kernels.RBF(DataMC2.D_in, variance=DataMC2.rng.rand())
m = mean_function_factory(DataMC2.rng, mean, DataMC2.D_in, DataMC2.D_out)
model = MomentMatchingSVGP(
DataMC2.X, DataMC2.Y, k, gpflow.likelihoods.Gaussian(),
mean_function=m, Z=DataMC2.X.copy(), whiten=white)
model.full_output_cov = True
gpflow.train.AdamOptimizer().minimize(model)
pred_mm = model.uncertain_predict_f_moment_matching(
tf.constant(DataMC2.Xnew_mu), tf.constant(DataMC2.Xnew_covar))
mean1, var1 = sess.run(pred_mm)
for n in range(DataMC2.N_new):
mean2, var2 = model.uncertain_predict_f_monte_carlo(
DataMC2.Xnew_mu[n, ...],
DataMC2.L[n, ...])
assert_almost_equal(mean1[n, ...], mean2, decimal=2)
assert_almost_equal(var1[n, ...], var2, decimal=2)
@pytest.mark.parametrize('mean', MEANS)
@pytest.mark.parametrize('white', [True, False])
def test_quadrature(white, mean):
with session_context() as session:
c = DataQuadrature
d = c.tensors(white, mean)
quad_args = d.Xmu, d.Xvar, c.H, c.D_in, (c.D_out,)
mean_quad = mvnquad(d.mean_fn, *quad_args)
var_quad = mvnquad(d.var_fn, *quad_args)
mean_sq_quad = mvnquad(d.mean_sq_fn, *quad_args)
mean_analytic, var_analytic = uncertain_conditional(
d.Xmu, d.Xvar, d.feat, d.kern,
d.q_mu, d.q_sqrt,
mean_function=d.mean_function,
full_output_cov=False,
white=white)
mean_quad, var_quad, mean_sq_quad = session.run(
[mean_quad, var_quad, mean_sq_quad], feed_dict=d.feed_dict)
var_quad = var_quad + (mean_sq_quad - mean_quad**2)
mean_analytic, var_analytic = session.run(
[mean_analytic, var_analytic], feed_dict=d.feed_dict)
assert_almost_equal(mean_quad, mean_analytic, decimal=6)
assert_almost_equal(var_quad, var_analytic, decimal=6)
if __name__ == "__main__":
tf.test.main()
