From bd1fb56e6c58eaf348322df4283e2ad8bfafad04 Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Fri, 2 Jan 2015 15:07:19 -0800 Subject: [PATCH] re-implemented warpedGP for new release of GPy --- GPy/kern/__init__.py | 2 +- GPy/models/warped_gp.py | 90 ++++++++++++++++------------------- GPy/util/warping_functions.py | 54 +++++++++++++-------- 3 files changed, 76 insertions(+), 70 deletions(-) diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py index c400277c..7a7c7ad8 100644 --- a/GPy/kern/__init__.py +++ b/GPy/kern/__init__.py @@ -1,7 +1,7 @@ from _src.kern import Kern from _src.rbf import RBF from _src.linear import Linear, LinearFull -from _src.static import Bias, White +from _src.static import Bias, White, Fixed from _src.brownian import Brownian from _src.stationary import Exponential, OU, Matern32, Matern52, ExpQuad, RatQuad, Cosine from _src.mlp import MLP diff --git a/GPy/models/warped_gp.py b/GPy/models/warped_gp.py index 4b982ed2..5bc9a417 100644 --- a/GPy/models/warped_gp.py +++ b/GPy/models/warped_gp.py @@ -1,7 +1,6 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) - import numpy as np from ..util.warping_functions import * from ..core import GP @@ -10,14 +9,16 @@ from GPy.util.warping_functions import TanhWarpingFunction_d from GPy import kern class WarpedGP(GP): - def __init__(self, X, Y, kernel=None, warping_function=None, warping_terms=3, normalize_X=False, normalize_Y=False): + def __init__(self, X, Y, kernel=None, warping_function=None, warping_terms=3): if kernel is None: - kernel = kern.rbf(X.shape[1]) + kernel = kern.RBF(X.shape[1]) if warping_function == None: self.warping_function = TanhWarpingFunction_d(warping_terms) self.warping_params = (np.random.randn(self.warping_function.n_terms * 3 + 1,) * 1) + else: + self.warping_function = warping_function self.scale_data = False if self.scale_data: @@ -25,10 +26,10 @@ class WarpedGP(GP): self.has_uncertain_inputs = False self.Y_untransformed = Y.copy() self.predict_in_warped_space = False - likelihood = likelihoods.Gaussian(self.transform_data(), normalize=normalize_Y) + likelihood = likelihoods.Gaussian() - GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X) - self._set_params(self._get_params()) + GP.__init__(self, X, self.transform_data(), likelihood=likelihood, kernel=kernel) + self.link_parameter(self.warping_function) def _scale_data(self, Y): self._Ymax = Y.max() @@ -38,62 +39,55 @@ class WarpedGP(GP): def _unscale_data(self, Y): return (Y + 0.5) * (self._Ymax - self._Ymin) + self._Ymin - def _set_params(self, x): - self.warping_params = x[:self.warping_function.num_parameters] - Y = self.transform_data() - self.likelihood.set_data(Y) - GP._set_params(self, x[self.warping_function.num_parameters:].copy()) + def parameters_changed(self): + self.Y[:] = self.transform_data() + super(WarpedGP, self).parameters_changed() - def _get_params(self): - return np.hstack((self.warping_params.flatten().copy(), GP._get_params(self).copy())) + Kiy = self.posterior.woodbury_vector.flatten() - def _get_param_names(self): - warping_names = self.warping_function._get_param_names() - param_names = GP._get_param_names(self) - return warping_names + param_names - - def transform_data(self): - Y = self.warping_function.f(self.Y_untransformed.copy(), self.warping_params).copy() - return Y - - def log_likelihood(self): - ll = GP.log_likelihood(self) - jacobian = self.warping_function.fgrad_y(self.Y_untransformed, self.warping_params) - return ll + np.log(jacobian).sum() - - def _log_likelihood_gradients(self): - ll_grads = GP._log_likelihood_gradients(self) - alpha = np.dot(self.Ki, self.likelihood.Y.flatten()) - warping_grads = self.warping_function_gradients(alpha) - - warping_grads = np.append(warping_grads[:, :-1].flatten(), warping_grads[0, -1]) - return np.hstack((warping_grads.flatten(), ll_grads.flatten())) - - def warping_function_gradients(self, Kiy): - grad_y = self.warping_function.fgrad_y(self.Y_untransformed, self.warping_params) - grad_y_psi, grad_psi = self.warping_function.fgrad_y_psi(self.Y_untransformed, self.warping_params, + grad_y = self.warping_function.fgrad_y(self.Y_untransformed) + grad_y_psi, grad_psi = self.warping_function.fgrad_y_psi(self.Y_untransformed, return_covar_chain=True) djac_dpsi = ((1.0 / grad_y[:, :, None, None]) * grad_y_psi).sum(axis=0).sum(axis=0) dquad_dpsi = (Kiy[:, None, None, None] * grad_psi).sum(axis=0).sum(axis=0) - return -dquad_dpsi + djac_dpsi + warping_grads = -dquad_dpsi + djac_dpsi + + self.warping_function.psi.gradient[:] = warping_grads[:, :-1] + self.warping_function.d.gradient[:] = warping_grads[0, -1] + + + def transform_data(self): + Y = self.warping_function.f(self.Y_untransformed.copy()).copy() + return Y + + def log_likelihood(self): + ll = GP.log_likelihood(self) + jacobian = self.warping_function.fgrad_y(self.Y_untransformed) + return ll + np.log(jacobian).sum() def plot_warping(self): - self.warping_function.plot(self.warping_params, self.Y_untransformed.min(), self.Y_untransformed.max()) + self.warping_function.plot(self.Y_untransformed.min(), self.Y_untransformed.max()) - def predict(self, Xnew, which_parts='all', full_cov=False, pred_init=None): + def predict(self, Xnew, which_parts='all', pred_init=None): # normalize X values - Xnew = (Xnew.copy() - self._Xoffset) / self._Xscale - mu, var = GP._raw_predict(self, Xnew, full_cov=full_cov, which_parts=which_parts) + # Xnew = (Xnew.copy() - self._Xoffset) / self._Xscale + mu, var = GP._raw_predict(self, Xnew) # now push through likelihood - mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov) + mean, var = self.likelihood.predictive_values(mu, var) if self.predict_in_warped_space: - mean = self.warping_function.f_inv(mean, self.warping_params, y=pred_init) - var = self.warping_function.f_inv(var, self.warping_params) + mean = self.warping_function.f_inv(mean, y=pred_init) + var = self.warping_function.f_inv(var) if self.scale_data: mean = self._unscale_data(mean) - - return mean, var, _025pm, _975pm + + return mean, var + +if __name__ == '__main__': + X = np.random.randn(100, 1) + Y = np.sin(X) + np.random.randn(100, 1)*0.05 + + m = WarpedGP(X, Y) diff --git a/GPy/util/warping_functions.py b/GPy/util/warping_functions.py index a0a385e0..a7547be6 100644 --- a/GPy/util/warping_functions.py +++ b/GPy/util/warping_functions.py @@ -1,17 +1,18 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) - import numpy as np +from GPy.core.parameterization import Parameterized, Param +from ..core.parameterization.transformations import Logexp -class WarpingFunction(object): +class WarpingFunction(Parameterized): """ abstract function for warping z = f(y) """ - def __init__(self): - raise NotImplementedError + def __init__(self, name): + super(WarpingFunction, self).__init__(name=name) def f(self,y,psi): """function transformation @@ -34,9 +35,10 @@ class WarpingFunction(object): def _get_param_names(self): raise NotImplementedError - def plot(self, psi, xmin, xmax): + def plot(self, xmin, xmax): + psi = self.psi y = np.arange(xmin, xmax, 0.01) - f_y = self.f(y, psi) + f_y = self.f(y) from matplotlib import pyplot as plt plt.figure() plt.plot(y, f_y) @@ -50,6 +52,7 @@ class TanhWarpingFunction(WarpingFunction): """n_terms specifies the number of tanh terms to be used""" self.n_terms = n_terms self.num_parameters = 3 * self.n_terms + super(TanhWarpingFunction, self).__init__(name='warp_tanh') def f(self,y,psi): """ @@ -163,8 +166,18 @@ class TanhWarpingFunction_d(WarpingFunction): """n_terms specifies the number of tanh terms to be used""" self.n_terms = n_terms self.num_parameters = 3 * self.n_terms + 1 + self.psi = np.ones((self.n_terms, 3)) - def f(self,y,psi): + super(TanhWarpingFunction_d, self).__init__(name='warp_tanh') + self.psi = Param('psi', self.psi) + self.psi[:, :2].constrain_positive() + + self.d = Param('%s' % ('d'), 1.0, Logexp()) + self.link_parameter(self.psi) + self.link_parameter(self.d) + + + def f(self,y): """ Transform y with f using parameter vector psi psi = [[a,b,c]] @@ -175,9 +188,9 @@ class TanhWarpingFunction_d(WarpingFunction): #1. check that number of params is consistent # assert psi.shape[0] == self.n_terms, 'inconsistent parameter dimensions' # assert psi.shape[1] == 4, 'inconsistent parameter dimensions' - mpsi = psi.copy() - d = psi[-1] - mpsi = mpsi[:self.num_parameters-1].reshape(self.n_terms, 3) + + d = self.d + mpsi = self.psi #3. transform data z = d*y.copy() @@ -187,7 +200,7 @@ class TanhWarpingFunction_d(WarpingFunction): return z - def f_inv(self, z, psi, max_iterations=1000, y=None): + def f_inv(self, z, max_iterations=1000, y=None): """ calculate the numerical inverse of f @@ -198,12 +211,12 @@ class TanhWarpingFunction_d(WarpingFunction): z = z.copy() if y is None: y = np.ones_like(z) - + it = 0 update = np.inf while it == 0 or (np.abs(update).sum() > 1e-10 and it < max_iterations): - update = (self.f(y, psi) - z)/self.fgrad_y(y, psi) + update = (self.f(y) - z)/self.fgrad_y(y) y -= update it += 1 if it == max_iterations: @@ -212,7 +225,7 @@ class TanhWarpingFunction_d(WarpingFunction): return y - def fgrad_y(self, y, psi, return_precalc = False): + def fgrad_y(self, y,return_precalc = False): """ gradient of f w.r.t to y ([N x 1]) @@ -221,9 +234,8 @@ class TanhWarpingFunction_d(WarpingFunction): """ - mpsi = psi.copy() - d = psi[-1] - mpsi = mpsi[:self.num_parameters-1].reshape(self.n_terms, 3) + d = self.d + mpsi = self.psi # vectorized version @@ -240,7 +252,7 @@ class TanhWarpingFunction_d(WarpingFunction): return GRAD - def fgrad_y_psi(self, y, psi, return_covar_chain = False): + def fgrad_y_psi(self, y, return_covar_chain = False): """ gradient of f w.r.t to y and psi @@ -248,10 +260,10 @@ class TanhWarpingFunction_d(WarpingFunction): """ - mpsi = psi.copy() - mpsi = mpsi[:self.num_parameters-1].reshape(self.n_terms, 3) - w, s, r, d = self.fgrad_y(y, psi, return_precalc = True) + mpsi = self.psi + + w, s, r, d = self.fgrad_y(y, return_precalc = True) gradients = np.zeros((y.shape[0], y.shape[1], len(mpsi), 4)) for i in range(len(mpsi)):