From 9680a139d43ac208063a309afc9bae36b4a46978 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 18 Mar 2014 12:28:46 +0000 Subject: [PATCH 1/2] changed the way the Gaussian likelihood interfaces, to enable mixed_noise things --- GPy/core/gp.py | 6 ++-- .../latent_function_inference/dtc.py | 32 ++++++++----------- .../exact_gaussian_inference.py | 5 ++- .../latent_function_inference/var_dtc.py | 4 +-- GPy/likelihoods/gaussian.py | 4 +-- 5 files changed, 23 insertions(+), 28 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 70b7d695..e052ff35 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -42,10 +42,8 @@ class GP(Model): assert Y.shape[0] == self.num_data _, self.output_dim = self.Y.shape - if Y_metadata is None: - Y_metadata = {} - else: - self.Y_metadata = Y_metadata + #TODO: check the type of this is okay? + self.Y_metadata = Y_metadata assert isinstance(kernel, kern.Kern) #assert self.input_dim == kernel.input_dim diff --git a/GPy/inference/latent_function_inference/dtc.py b/GPy/inference/latent_function_inference/dtc.py index 5ebc5e53..89140ce2 100644 --- a/GPy/inference/latent_function_inference/dtc.py +++ b/GPy/inference/latent_function_inference/dtc.py @@ -91,12 +91,8 @@ class vDTC(object): def __init__(self): self.const_jitter = 1e-6 - def inference(self, kern, X, X_variance, Z, likelihood, Y): - assert X_variance is None, "cannot use X_variance with DTC. Try varDTC." - - #TODO: MAX! fix this! - from ...util.misc import param_to_array - Y = param_to_array(Y) + def inference(self, kern, X, Z, likelihood, Y): + #assert X_variance is None, "cannot use X_variance with DTC. Try varDTC." num_inducing, _ = Z.shape num_data, output_dim = Y.shape @@ -109,15 +105,14 @@ class vDTC(object): Kmm = kern.K(Z) Knn = kern.Kdiag(X) Knm = kern.K(X, Z) - U = Knm - Uy = np.dot(U.T,Y) + KnmY = np.dot(Knm.T,Y) - #factor Kmm + #factor Kmm Kmmi, L, Li, _ = pdinv(Kmm) # Compute A - LiUTbeta = np.dot(Li, U.T)*np.sqrt(beta) - A_ = tdot(LiUTbeta) + LiKmnbeta = np.dot(Li, Knm.T)*np.sqrt(beta) + A_ = tdot(LiKmnbeta) trace_term = -0.5*(np.sum(Knn)*beta - np.trace(A_)) A = A_ + np.eye(num_inducing) @@ -125,7 +120,7 @@ class vDTC(object): LA = jitchol(A) # back substutue to get b, P, v - tmp, _ = dtrtrs(L, Uy, lower=1) + tmp, _ = dtrtrs(L, KnmY, lower=1) b, _ = dtrtrs(LA, tmp*beta, lower=1) tmp, _ = dtrtrs(LA, b, lower=1, trans=1) v, _ = dtrtrs(L, tmp, lower=1, trans=1) @@ -145,19 +140,18 @@ class vDTC(object): LAL = Li.T.dot(A).dot(Li) dL_dK = Kmmi - 0.5*(vvT_P + LAL) - # Compute dL_dU + # Compute dL_dKnm vY = np.dot(v.reshape(-1,1),Y.T) - #dL_dU = vY - np.dot(vvT_P, U.T) - dL_dU = vY - np.dot(vvT_P - Kmmi, U.T) - dL_dU *= beta + dL_dKmn = vY - np.dot(vvT_P - Kmmi, Knm.T) + dL_dKmn *= beta #compute dL_dR - Uv = np.dot(U, v) - dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - 1./beta + np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1) )*beta**2 + Knmv = np.dot(Knm, v) + dL_dR = 0.5*(np.sum(Knm*np.dot(Knm,P), 1) - 1./beta + np.sum(np.square(Y), 1) - 2.*np.sum(Knmv*Y, 1) + np.sum(np.square(Knmv), 1) )*beta**2 dL_dR -=beta*trace_term/num_data dL_dthetaL = likelihood.exact_inference_gradients(dL_dR) - grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':np.zeros_like(Knn) + -0.5*beta, 'dL_dKnm':dL_dU.T, 'dL_dthetaL':dL_dthetaL} + grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':np.zeros_like(Knn) + -0.5*beta, 'dL_dKnm':dL_dKmn.T, 'dL_dthetaL':dL_dthetaL} #construct a posterior object post = Posterior(woodbury_inv=Kmmi-P, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=L) diff --git a/GPy/inference/latent_function_inference/exact_gaussian_inference.py b/GPy/inference/latent_function_inference/exact_gaussian_inference.py index e76575c6..ca1b92c6 100644 --- a/GPy/inference/latent_function_inference/exact_gaussian_inference.py +++ b/GPy/inference/latent_function_inference/exact_gaussian_inference.py @@ -3,6 +3,7 @@ from posterior import Posterior from ...util.linalg import pdinv, dpotrs, tdot +from ...util import diag import numpy as np log_2_pi = np.log(2*np.pi) @@ -41,7 +42,9 @@ class ExactGaussianInference(object): K = kern.K(X) - Wi, LW, LWi, W_logdet = pdinv(K + likelihood.covariance_matrix(Y, Y_metadata)) + Ky = K.copy() + diag.add(Ky, likelihood.gaussian_variance(Y, Y_metadata)) + Wi, LW, LWi, W_logdet = pdinv(Ky) alpha, _ = dpotrs(LW, YYT_factor, lower=1) diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index 9b6e26c0..97d54624 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -48,7 +48,7 @@ class VarDTC(object): def get_VVTfactor(self, Y, prec): return Y * prec # TODO chache this, and make it effective - def inference(self, kern, X, Z, likelihood, Y): + def inference(self, kern, X, Z, likelihood, Y, Y_metadata=None): if isinstance(X, VariationalPosterior): uncertain_inputs = True psi0 = kern.psi0(Z, X) @@ -65,7 +65,7 @@ class VarDTC(object): _, output_dim = Y.shape #see whether we've got a different noise variance for each datum - beta = 1./np.fmax(likelihood.gaussian_variance(Y_metadata), 1e-6) + beta = 1./np.fmax(likelihood.gaussian_variance(Y, Y_metadata), 1e-6) # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency! #self.YYTfactor = self.get_YYTfactor(Y) #VVT_factor = self.get_VVTfactor(self.YYTfactor, beta) diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py index 101aac4b..79d62bb7 100644 --- a/GPy/likelihoods/gaussian.py +++ b/GPy/likelihoods/gaussian.py @@ -50,8 +50,8 @@ class Gaussian(Likelihood): if isinstance(gp_link, link_functions.Identity): self.log_concave = True - def covariance_matrix(self, Y, Y_metadata=None): - return np.eye(Y.shape[0]) * self.variance + def gaussian_variance(self, Y, Y_metadata=None): + return self.variance def update_gradients(self, grad): self.variance.gradient = grad From 24b43c490caa1d22703959e537ada28edb74cae2 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 18 Mar 2014 16:30:46 +0000 Subject: [PATCH 2/2] fixes now hierarchical, maybe need to be restructured as lookup from constraints --- GPy/core/parameterization/param.py | 2 +- GPy/core/parameterization/parameter_core.py | 47 ++++++++++++++------- GPy/testing/parameterized_tests.py | 46 +++++++++++++------- 3 files changed, 62 insertions(+), 33 deletions(-) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 324593f9..b73e7dfa 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -226,7 +226,7 @@ class Param(OptimizationHandlable, ObsAr): # Constrainable #=========================================================================== def _ensure_fixes(self): - self._fixes_ = numpy.ones(self._realsize_, dtype=bool) + if not self._has_fixes(): self._fixes_ = numpy.ones(self._realsize_, dtype=bool) #=========================================================================== # Convenience diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 0aab890c..d4779127 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -16,7 +16,7 @@ Observable Pattern for patameterization from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED import numpy as np -__updated__ = '2014-03-17' +__updated__ = '2014-03-18' class HierarchyError(Exception): """ @@ -377,7 +377,7 @@ class Constrainable(Nameable, Indexable): # Ensure that the fixes array is set: # Parameterized: ones(self.size) # Param: ones(self._realsize_ - self._fixes_ = np.ones(self.size, dtype=bool) + if not self._has_fixes(): self._fixes_ = np.ones(self.size, dtype=bool) def _set_fixed(self, index): self._ensure_fixes() @@ -398,7 +398,7 @@ class Constrainable(Nameable, Indexable): self._fixes_ = None def _has_fixes(self): - return hasattr(self, "_fixes_") and self._fixes_ is not None + return hasattr(self, "_fixes_") and self._fixes_ is not None and self._fixes_.size == self.size #=========================================================================== # Prior Operations @@ -576,14 +576,22 @@ class OptimizationHandlable(Constrainable): # transformed parameters (apply transformation rules) p = self._param_array_.copy() [np.put(p, ind, c.finv(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] - if self._has_fixes(): + if self.has_parent() and self.constraints[__fixed__].size != 0: + fixes = np.ones(self.size).astype(bool) + fixes[self.constraints[__fixed__]] = FIXED + return p[fixes] + elif self._has_fixes(): return p[self._fixes_] return p def _set_params_transformed(self, p): if p is self._param_array_: p = p.copy() - if self._has_fixes(): self._param_array_[self._fixes_] = p + if self.has_parent() and self.constraints[__fixed__].size != 0: + fixes = np.ones(self.size).astype(bool) + fixes[self.constraints[__fixed__]] = FIXED + self._param_array_[fixes] = p + elif self._has_fixes(): self._param_array_[self._fixes_] = p else: self._param_array_[:] = p self.untransform() self._trigger_params_changed() @@ -770,11 +778,11 @@ class Parameterizable(OptimizationHandlable): Add all parameters to this param class, you can insert parameters at any given index using the :func:`list.insert` syntax """ - # if param.has_parent(): - # raise AttributeError, "parameter {} already in another model, create new object (or copy) for adding".format(param._short()) if param in self._parameters_ and index is not None: self.remove_parameter(param) self.add_parameter(param, index) + elif param.has_parent(): + raise HierarchyError, "parameter {} already in another model ({}), create new object (or copy) for adding".format(param._short(), param._highest_parent_._short()) elif param not in self._parameters_: if param.has_parent(): parent = param._parent_ @@ -798,13 +806,19 @@ class Parameterizable(OptimizationHandlable): param.add_observer(self, self._pass_through_notify_observers, -np.inf) - self.size += param.size + parent = self + while parent is not None: + parent.size += param.size + parent = parent._parent_ + + self._connect_parameters() + + self._highest_parent_._connect_parameters(ignore_added_names=_ignore_added_names) + self._highest_parent_._notify_parent_change() + self._highest_parent_._connect_fixes() - self._connect_parameters(ignore_added_names=_ignore_added_names) - self._notify_parent_change() - self._connect_fixes() else: - raise RuntimeError, """Parameter exists already added and no copy made""" + raise HierarchyError, """Parameter exists already and no copy made""" def add_parameters(self, *parameters): @@ -830,17 +844,18 @@ class Parameterizable(OptimizationHandlable): param.remove_observer(self, self._pass_through_notify_observers) self.constraints.shift_left(start, param.size) - self._connect_fixes() self._connect_parameters() self._notify_parent_change() parent = self._parent_ while parent is not None: - parent._connect_fixes() - parent._connect_parameters() - parent._notify_parent_change() + parent.size -= param.size parent = parent._parent_ + self._highest_parent_._connect_parameters() + self._highest_parent_._connect_fixes() + self._highest_parent_._notify_parent_change() + def _connect_parameters(self, ignore_added_names=False): # connect parameterlist to this parameterized object # This just sets up the right connection for the params objects diff --git a/GPy/testing/parameterized_tests.py b/GPy/testing/parameterized_tests.py index 81c2dfdd..cd5127c8 100644 --- a/GPy/testing/parameterized_tests.py +++ b/GPy/testing/parameterized_tests.py @@ -34,9 +34,9 @@ class ParameterizedTest(unittest.TestCase): self.param = Param('param', np.random.rand(25,2), Logistic(0, 1)) self.test1 = GPy.core.Parameterized("test model") - self.test1.add_parameter(self.white) - self.test1.add_parameter(self.rbf, 0) - self.test1.add_parameter(self.param) + self.test1.kern = self.rbf+self.white + self.test1.add_parameter(self.test1.kern) + self.test1.add_parameter(self.param, 0) x = np.linspace(-2,6,4)[:,None] y = np.sin(x) @@ -45,22 +45,24 @@ class ParameterizedTest(unittest.TestCase): def test_add_parameter(self): self.assertEquals(self.rbf._parent_index_, 0) self.assertEquals(self.white._parent_index_, 1) + self.assertEquals(self.param._parent_index_, 0) pass def test_fixes(self): self.white.fix(warning=False) - self.test1.remove_parameter(self.test1.param) + self.test1.remove_parameter(self.param) self.assertTrue(self.test1._has_fixes()) from GPy.core.parameterization.transformations import FIXED, UNFIXED self.assertListEqual(self.test1._fixes_.tolist(),[UNFIXED,UNFIXED,FIXED]) - - self.test1.add_parameter(self.white, 0) + self.test1.kern.add_parameter(self.white, 0) self.assertListEqual(self.test1._fixes_.tolist(),[FIXED,UNFIXED,UNFIXED]) + self.test1.kern.rbf.fix() + self.assertListEqual(self.test1._fixes_.tolist(),[FIXED]*3) def test_remove_parameter(self): from GPy.core.parameterization.transformations import FIXED, UNFIXED, __fixed__, Logexp self.white.fix() - self.test1.remove_parameter(self.white) + self.test1.kern.remove_parameter(self.white) self.assertIs(self.test1._fixes_,None) self.assertListEqual(self.white._fixes_.tolist(), [FIXED]) @@ -81,7 +83,12 @@ class ParameterizedTest(unittest.TestCase): self.assertListEqual(self.white._fixes_.tolist(), [FIXED]) self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops) - self.assertListEqual(self.test1.constraints[Logexp()].tolist(), [0,1]) + self.assertListEqual(self.test1.constraints[Logexp()].tolist(), range(self.param.size, self.param.size+self.rbf.size)) + + def test_remove_parameter_param_array_grad_array(self): + val = self.test1.kern._param_array_.copy() + self.test1.kern.remove_parameter(self.white) + self.assertListEqual(self.test1.kern._param_array_.tolist(), val[:2].tolist()) def test_add_parameter_already_in_hirarchy(self): self.assertRaises(HierarchyError, self.test1.add_parameter, self.white._parameters_[0]) @@ -91,28 +98,35 @@ class ParameterizedTest(unittest.TestCase): self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) self.assertListEqual(self.rbf.constraints.indices()[0].tolist(), range(2)) from GPy.core.parameterization.transformations import Logexp - kern = self.rbf+self.white + kern = self.test1.kern + self.test1.remove_parameter(kern) self.assertListEqual(kern.constraints[Logexp()].tolist(), range(3)) def test_constraints(self): self.rbf.constrain(GPy.transformations.Square(), False) - self.assertListEqual(self.test1.constraints[GPy.transformations.Square()].tolist(), range(2)) - self.assertListEqual(self.test1.constraints[GPy.transformations.Logexp()].tolist(), [2]) + self.assertListEqual(self.test1.constraints[GPy.transformations.Square()].tolist(), range(self.param.size, self.param.size+self.rbf.size)) + self.assertListEqual(self.test1.constraints[GPy.transformations.Logexp()].tolist(), [self.param.size+self.rbf.size]) - self.test1.remove_parameter(self.rbf) + self.test1.kern.remove_parameter(self.rbf) self.assertListEqual(self.test1.constraints[GPy.transformations.Square()].tolist(), []) def test_constraints_views(self): - self.assertEqual(self.white.constraints._offset, 2) - self.assertEqual(self.rbf.constraints._offset, 0) - self.assertEqual(self.param.constraints._offset, 3) + self.assertEqual(self.white.constraints._offset, self.param.size+self.rbf.size) + self.assertEqual(self.rbf.constraints._offset, self.param.size) + self.assertEqual(self.param.constraints._offset, 0) def test_fixing_randomize(self): self.white.fix(warning=True) - val = float(self.test1.white.variance) + val = float(self.white.variance) self.test1.randomize() self.assertEqual(val, self.white.variance) + def test_fixing_randomize_parameter_handling(self): + self.rbf.fix(warning=True) + val = float(self.rbf.variance) + self.test1.kern.randomize() + self.assertEqual(val, self.rbf.variance) + def test_fixing_optimize(self): self.testmodel.kern.lengthscale.fix() val = float(self.testmodel.kern.lengthscale)