From db5fd17609346b56c14ce07b32fa1268abcdd007 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 7 Mar 2014 16:59:41 +0000 Subject: [PATCH 01/19] slicing support for kernel input dimension --- GPy/core/gp.py | 5 +- GPy/core/parameterization/array_core.py | 2 +- GPy/core/parameterization/parameter_core.py | 6 +- GPy/core/parameterization/variational.py | 16 ++- GPy/core/sparse_gp.py | 8 +- .../exact_gaussian_inference.py | 3 - GPy/kern/_src/add.py | 54 +++++---- GPy/kern/_src/kern.py | 106 ++++++++++++++++-- GPy/kern/_src/stationary.py | 8 +- GPy/util/caching.py | 35 +++--- 10 files changed, 178 insertions(+), 65 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 1add8268..6441561b 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -48,7 +48,7 @@ class GP(Model): self.Y_metadata = None assert isinstance(kernel, kern.Kern) - assert self.input_dim == kernel.input_dim + #assert self.input_dim == kernel.input_dim self.kern = kernel assert isinstance(likelihood, likelihoods.Likelihood) @@ -68,8 +68,9 @@ class GP(Model): def parameters_changed(self): self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y, Y_metadata=self.Y_metadata) + self.likelihood.update_gradients(np.diag(grad_dict['dL_dK'])) self.kern.update_gradients_full(grad_dict['dL_dK'], self.X) - + def log_likelihood(self): return self._log_marginal_likelihood diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index 27801e23..9ce0e8f6 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -16,7 +16,7 @@ class ObservableArray(np.ndarray, Observable): __array_priority__ = -1 # Never give back ObservableArray def __new__(cls, input_array): if not isinstance(input_array, ObservableArray): - obj = np.atleast_1d(input_array).view(cls) + obj = np.atleast_1d(np.require(input_array, dtype=np.float64, requirements=['W', 'C'])).view(cls) else: obj = input_array cls.__name__ = "ObservableArray\n " return obj diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index a78cf02d..351eacef 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -15,7 +15,6 @@ Observable Pattern for patameterization from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED import numpy as np -import itertools __updated__ = '2013-12-16' @@ -43,6 +42,7 @@ class Observable(object): _updated = True def __init__(self, *args, **kwargs): self._observer_callables_ = [] + def __del__(self, *args, **kwargs): del self._observer_callables_ @@ -551,8 +551,8 @@ class OptimizationHandlable(Constrainable, Observable): return p def _set_params_transformed(self, p): - if p is self._param_array_: - p = p.copy() + #if p is self._param_array_: + p = p.copy() if self._has_fixes(): self._param_array_[self._fixes_] = p else: self._param_array_[:] = p self.untransform() diff --git a/GPy/core/parameterization/variational.py b/GPy/core/parameterization/variational.py index 8bc7ca59..4c929cc8 100644 --- a/GPy/core/parameterization/variational.py +++ b/GPy/core/parameterization/variational.py @@ -66,10 +66,10 @@ class VariationalPosterior(Parameterized): def __init__(self, means=None, variances=None, name=None, **kw): super(VariationalPosterior, self).__init__(name=name, **kw) self.mean = Param("mean", means) - self.ndim = self.mean.ndim - self.shape = self.mean.shape self.variance = Param("variance", variances, Logexp()) self.add_parameters(self.mean, self.variance) + self.ndim = self.mean.ndim + self.shape = self.mean.shape self.num_data, self.input_dim = self.mean.shape if self.has_uncertain_inputs(): assert self.variance.shape == self.mean.shape, "need one variance per sample and dimenion" @@ -77,6 +77,18 @@ class VariationalPosterior(Parameterized): def has_uncertain_inputs(self): return not self.variance is None + def __getitem__(self, s): + import copy + n = self.__new__(self.__class__) + dc = copy.copy(self.__dict__) + dc['mean'] = dc['mean'][s] + dc['variance'] = dc['variance'][s] + dc['shape'] = dc['mean'].shape + dc['ndim'] = dc['ndim'] + dc['num_data'], dc['input_dim'] = self.mean.shape[0], self.mean.shape[1] if dc['ndim'] > 1 else 1 + n.__dict__ = dc + return n + class NormalPosterior(VariationalPosterior): ''' diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index f4f34a5e..16b66676 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -64,8 +64,8 @@ class SparseGP(GP): self.kern.gradient += target #gradients wrt Z - self.Z.gradient = self.kern.gradients_X(dL_dKmm, self.Z) - self.Z.gradient += self.kern.gradients_Z_expectations( + self.Z.gradient[:,self.kern.active_dims] = self.kern.gradients_X(dL_dKmm, self.Z) + self.Z.gradient[:,self.kern.active_dims] += self.kern.gradients_Z_expectations( self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X) else: #gradients wrt kernel @@ -77,8 +77,8 @@ class SparseGP(GP): self.kern.gradient += target #gradients wrt Z - self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) - self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) + self.Z.gradient[:,self.kern.active_dims] = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) + self.Z.gradient[:,self.kern.active_dims] += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) def _raw_predict(self, Xnew, full_cov=False): """ diff --git a/GPy/inference/latent_function_inference/exact_gaussian_inference.py b/GPy/inference/latent_function_inference/exact_gaussian_inference.py index 922b52f4..47f6ea09 100644 --- a/GPy/inference/latent_function_inference/exact_gaussian_inference.py +++ b/GPy/inference/latent_function_inference/exact_gaussian_inference.py @@ -49,9 +49,6 @@ class ExactGaussianInference(object): dL_dK = 0.5 * (tdot(alpha) - Y.shape[1] * Wi) - #TODO: does this really live here? - likelihood.update_gradients(np.diag(dL_dK)) - return Posterior(woodbury_chol=LW, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK} diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py index 77fe057d..87dda365 100644 --- a/GPy/kern/_src/add.py +++ b/GPy/kern/_src/add.py @@ -1,12 +1,10 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) -import sys import numpy as np import itertools -from linear import Linear from ...core.parameterization import Parameterized -from ...core.parameterization.param import Param +from ...util.caching import Cache_this from kern import Kern class Add(Kern): @@ -14,19 +12,24 @@ class Add(Kern): assert all([isinstance(k, Kern) for k in subkerns]) if tensor: input_dim = sum([k.input_dim for k in subkerns]) - self.input_slices = [] + self.self.active_dims = [] n = 0 for k in subkerns: - self.input_slices.append(slice(n, n+k.input_dim)) + self.self.active_dims.append(slice(n, n+k.input_dim)) n += k.input_dim else: - assert all([k.input_dim == subkerns[0].input_dim for k in subkerns]) - input_dim = subkerns[0].input_dim - self.input_slices = [slice(None) for k in subkerns] + #assert all([k.input_dim == subkerns[0].input_dim for k in subkerns]) + #input_dim = subkerns[0].input_dim + #self.input_slices = [slice(None) for k in subkerns] + input_dim = reduce(np.union1d, map(lambda x: np.r_[x.active_dims], subkerns)) super(Add, self).__init__(input_dim, 'add') self.add_parameters(*subkerns) - - + + @property + def parts(self): + return self._parameters_ + + @Cache_this(limit=1, force_kwargs=('which_parts',)) def K(self, X, X2=None): """ Compute the kernel function. @@ -37,13 +40,19 @@ class Add(Kern): handLes this as X2 == X. """ assert X.shape[1] == self.input_dim - if X2 is None: - return sum([p.K(X[:, i_s], None) for p, i_s in zip(self._parameters_, self.input_slices)]) - else: - return sum([p.K(X[:, i_s], X2[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]) + which_parts=None + if which_parts is None: + which_parts = self.parts + elif not isinstance(which_parts, (list, tuple)): + # if only one part is given + which_parts = [which_parts] + return sum([p.K(X, X2) for p in which_parts]) - def update_gradients_full(self, dL_dK, X): - [p.update_gradients_full(dL_dK, X[:,i_s]) for p, i_s in zip(self._parameters_, self.input_slices)] + def update_gradients_full(self, dL_dK, X, X2=None): + [p.update_gradients_full(dL_dK, X, X2) for p in self.parts] + + def update_gradients_diag(self, dL_dK, X): + [p.update_gradients_diag(dL_dK, X) for p in self.parts] def gradients_X(self, dL_dK, X, X2=None): """Compute the gradient of the objective function with respect to X. @@ -55,16 +64,17 @@ class Add(Kern): :param X2: Observed data inputs (optional, defaults to X) :type X2: np.ndarray (num_inducing x input_dim)""" - target = np.zeros_like(X) - if X2 is None: - [np.add(target[:,i_s], p.gradients_X(dL_dK, X[:, i_s], None), target[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)] - else: - [np.add(target[:,i_s], p.gradients_X(dL_dK, X[:, i_s], X2[:,i_s]), target[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)] + target = np.zeros(X.shape) + for p in self.parts: + target[:, p.active_dims] += p.gradients_X(dL_dK, X, X2) return target def Kdiag(self, X): + which_parts=None assert X.shape[1] == self.input_dim - return sum([p.Kdiag(X[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]) + if which_parts is None: + which_parts = self.parts + return sum([p.Kdiag(X) for p in which_parts]) def psi0(self, Z, variational_posterior): diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index 47166156..33b9ff08 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -2,13 +2,22 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) import sys -import numpy as np -import itertools -from ...core.parameterization import Parameterized -from ...core.parameterization.param import Param - +from ...core.parameterization.parameterized import ParametersChangedMeta, Parameterized +from ...util.caching import Cache_this +class KernCallsViaSlicerMeta(ParametersChangedMeta): + def __call__(self, *args, **kw): + instance = super(KernCallsViaSlicerMeta, self).__call__(*args, **kw) + instance.K = instance._slice_wrapper(instance.K) + instance.Kdiag = instance._slice_wrapper(instance.Kdiag, True) + instance.update_gradients_full = instance._slice_wrapper(instance.update_gradients_full, False, True) + instance.update_gradients_diag = instance._slice_wrapper(instance.update_gradients_diag, True, True) + instance.gradients_X = instance._slice_wrapper(instance.gradients_X, False, True) + instance.gradients_X_diag = instance._slice_wrapper(instance.gradients_X_diag, True, True) + return instance + class Kern(Parameterized): + __metaclass__ = KernCallsViaSlicerMeta def __init__(self, input_dim, name, *a, **kw): """ The base class for a kernel: a positive definite function @@ -20,11 +29,83 @@ class Kern(Parameterized): Do not instantiate. """ super(Kern, self).__init__(name=name, *a, **kw) - self.input_dim = input_dim - + if isinstance(input_dim, int): + self.active_dims = slice(0, input_dim) + self.input_dim = input_dim + else: + self.active_dims = input_dim + self.input_dim = len(self.active_dims) + self._sliced_X = False + self._sliced_X2 = False + + @Cache_this(limit=10, ignore_args = (0,)) + def _slice_X(self, X): + return X[:, self.active_dims] + + def _slice_wrapper(self, operation, diag=False, derivative=False): + """ + This method wraps the functions in kernel to make sure all kernels allways see their respective input dimension. + The different switches are: + diag: if X2 exists + derivative: if firest arg is dL_dK + """ + if derivative: + if diag: + def x_slice_wrapper(dL_dK, X, *args, **kw): + X = self._slice_X(X) if not self._sliced_X else X + self._sliced_X = True + try: + ret = operation(dL_dK, X, *args, **kw) + except: raise + finally: + self._sliced_X = False + return ret + else: + def x_slice_wrapper(dL_dK, X, X2=None, *args, **kw): + X, X2 = self._slice_X(X) if not self._sliced_X else X, self._slice_X(X2) if X2 is not None and not self._sliced_X2 else X2 + self._sliced_X = True + self._sliced_X2 = True + try: + ret = operation(dL_dK, X, X2, *args, **kw) + except: raise + finally: + self._sliced_X = False + self._sliced_X2 = False + return ret + else: + if diag: + def x_slice_wrapper(X, *args, **kw): + X = self._slice_X(X) if not self._sliced_X else X + self._sliced_X = True + try: + ret = operation(X, *args, **kw) + except: raise + finally: + self._sliced_X = False + return ret + else: + def x_slice_wrapper(X, X2=None, *args, **kw): + X, X2 = self._slice_X(X) if not self._sliced_X else X, self._slice_X(X2) if X2 is not None and not self._sliced_X2 else X2 + self._sliced_X = True + self._sliced_X2 = True + try: + ret = operation(X, X2, *args, **kw) + except: raise + finally: + self._sliced_X = False + self._sliced_X2 = False + return ret + x_slice_wrapper._operation = operation + x_slice_wrapper.__name__ = ("slicer("+operation.__name__ + +(","+str(bool(diag)) if diag else'') + +(','+str(bool(derivative)) if derivative else '') + +')') + x_slice_wrapper.__doc__ = "**sliced**\n\n" + (operation.__doc__ or "") + return x_slice_wrapper + def K(self, X, X2): raise NotImplementedError - def Kdiag(self, Xa): + def Kdiag(self, X): raise NotImplementedError def psi0(self, Z, variational_posterior): raise NotImplementedError @@ -34,13 +115,16 @@ class Kern(Parameterized): raise NotImplementedError def gradients_X(self, dL_dK, X, X2): raise NotImplementedError - def gradients_X_diag(self, dL_dK, X): + def gradients_X_diag(self, dL_dKdiag, X): raise NotImplementedError - + def update_gradients_full(self, dL_dK, X, X2): """Set the gradients of all parameters when doing full (N) inference.""" raise NotImplementedError - + def update_gradients_diag(self, dL_dKdiag, X): + """Set the gradients for all parameters for the derivative of the diagonal of the covariance w.r.t the kernel parameters.""" + raise NotImplementedError + def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): """ Set the gradients of all parameters when doing inference with diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index 44e17d8a..725f8660 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -57,7 +57,7 @@ class Stationary(Kern): if lengthscale.size != input_dim: lengthscale = np.ones(input_dim)*lengthscale else: - lengthscale = np.ones(self.input_dim) + lengthscale = np.ones(self.input_dim) self.lengthscale = Param('lengthscale', lengthscale, Logexp()) self.variance = Param('variance', variance, Logexp()) assert self.variance.size==1 @@ -85,12 +85,14 @@ class Stationary(Kern): Compute the Euclidean distance between each row of X and X2, or between each pair of rows of X if X2 is None. """ + #X, = self._slice_X(X) if X2 is None: Xsq = np.sum(np.square(X),1) r2 = -2.*tdot(X) + (Xsq[:,None] + Xsq[None,:]) util.diag.view(r2)[:,]= 0. # force diagnoal to be zero: sometime numerically a little negative return np.sqrt(r2) else: + #X2, = self._slice_X(X2) X1sq = np.sum(np.square(X),1) X2sq = np.sum(np.square(X2),1) return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])) @@ -124,7 +126,6 @@ class Stationary(Kern): self.lengthscale.gradient = 0. def update_gradients_full(self, dL_dK, X, X2=None): - self.variance.gradient = np.einsum('ij,ij,i', self.K(X, X2), dL_dK, 1./self.variance) #now the lengthscale gradient(s) @@ -136,7 +137,7 @@ class Stationary(Kern): #self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3 tmp = dL_dr*self._inv_dist(X, X2) if X2 is None: X2 = X - self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(X[:,q:q+1] - X2[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in xrange(self.input_dim)]) + self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(self._slice_X(X)[:,q:q+1] - self._slice_X(X2)[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in xrange(self.input_dim)]) else: r = self._scaled_dist(X, X2) self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale @@ -176,7 +177,6 @@ class Stationary(Kern): ret = np.empty(X.shape, dtype=np.float64) [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)] ret /= self.lengthscale**2 - return ret def gradients_X_diag(self, dL_dKdiag, X): diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 250efe11..0b6f7234 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -9,24 +9,27 @@ class Cacher(object): """ - def __init__(self, operation, limit=5, ignore_args=()): + def __init__(self, operation, limit=5, ignore_args=(), force_kwargs=()): self.limit = int(limit) self.ignore_args = ignore_args + self.force_kwargs = force_kwargs self.operation=operation self.cached_inputs = [] self.cached_outputs = [] self.inputs_changed = [] - def __call__(self, *args): + def __call__(self, *args, **kw): """ A wrapper function for self.operation, """ #ensure that specified arguments are ignored + items = sorted(kw.items(), key=lambda x: x[0]) + oa_all = args + tuple(a for _,a in items) if len(self.ignore_args) != 0: - oa = [a for i,a in enumerate(args) if i not in self.ignore_args] + oa = [a for i,a in itertools.chain(enumerate(args), items) if i not in self.ignore_args and i not in self.force_kwargs] else: - oa = args + oa = oa_all # this makes sure we only add an observer once, and that None can be in args observable_args = [] @@ -37,8 +40,13 @@ class Cacher(object): #make sure that all the found argument really are observable: #otherswise don't cache anything, pass args straight though if not all([isinstance(arg, Observable) for arg in observable_args]): - return self.operation(*args) + return self.operation(*args, **kw) + if len(self.force_kwargs) != 0: + # check if there are force args, which force reloading + for k in self.force_kwargs: + if k in kw and kw[k] is not None: + return self.operation(*args, **kw) # TODO: WARNING !!! Cache OFFSWITCH !!! WARNING # return self.operation(*args) @@ -48,7 +56,7 @@ class Cacher(object): i = state.index(True) if self.inputs_changed[i]: #(elements of) the args have changed since we last computed: update - self.cached_outputs[i] = self.operation(*args) + self.cached_outputs[i] = self.operation(*args, **kw) self.inputs_changed[i] = False return self.cached_outputs[i] else: @@ -62,11 +70,11 @@ class Cacher(object): self.cached_outputs.pop(0) #compute - self.cached_inputs.append(args) - self.cached_outputs.append(self.operation(*args)) + self.cached_inputs.append(oa_all) + self.cached_outputs.append(self.operation(*args, **kw)) self.inputs_changed.append(False) [a.add_observer(self, self.on_cache_changed) for a in observable_args] - return self.cached_outputs[-1]#Max says return. + return self.cached_outputs[-1]#return def on_cache_changed(self, arg): """ @@ -90,15 +98,16 @@ class Cache_this(object): """ A decorator which can be applied to bound methods in order to cache them """ - def __init__(self, limit=5, ignore_args=()): + def __init__(self, limit=5, ignore_args=(), force_kwargs=()): self.limit = limit self.ignore_args = ignore_args + self.force_args = force_kwargs self.c = None def __call__(self, f): - def f_wrap(*args): + def f_wrap(*args, **kw): if self.c is None: - self.c = Cacher(f, self.limit, ignore_args=self.ignore_args) - return self.c(*args) + self.c = Cacher(f, self.limit, ignore_args=self.ignore_args, force_kwargs=self.force_args) + return self.c(*args, **kw) f_wrap._cacher = self f_wrap.__doc__ = "**cached**\n\n" + (f.__doc__ or "") return f_wrap From 7e9078b0f9f58a539a1153622b24874a235e21b6 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 10 Mar 2014 16:01:32 +0000 Subject: [PATCH 02/19] merged params here --- GPy/kern/_src/add.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py index 3514a224..604ed103 100644 --- a/GPy/kern/_src/add.py +++ b/GPy/kern/_src/add.py @@ -24,6 +24,9 @@ class Add(Kern): super(Add, self).__init__(input_dim, 'add') self.add_parameters(*subkerns) + @property + def parts(self): + return self._parameters_ def K(self, X, X2=None): """ @@ -107,8 +110,6 @@ class Add(Kern): def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): from static import White, Bias - mu, S = variational_posterior.mean, variational_posterior.variance - for p1, is1 in zip(self._parameters_, self.input_slices): #compute the effective dL_dpsi1. Extra terms appear becaue of the cross terms in psi2! @@ -129,7 +130,6 @@ class Add(Kern): def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): from static import White, Bias - target = np.zeros(Z.shape) for p1, is1 in zip(self._parameters_, self.input_slices): @@ -151,7 +151,6 @@ class Add(Kern): def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): from static import White, Bias - target_mu = np.zeros(variational_posterior.shape) target_S = np.zeros(variational_posterior.shape) for p1, is1 in zip(self._parameters_, self.input_slices): From 2d8246d33f779823ba4b5bf8060c855c888f5147 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 11 Mar 2014 10:24:15 +0000 Subject: [PATCH 03/19] Combination Kernel for add and prod --- GPy/kern/_src/add.py | 51 +++++++------------------- GPy/kern/_src/kern.py | 83 +++++++++++++++++++++++++++++-------------- GPy/kern/_src/prod.py | 51 +++++++++++++------------- 3 files changed, 94 insertions(+), 91 deletions(-) diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py index 604ed103..433a8921 100644 --- a/GPy/kern/_src/add.py +++ b/GPy/kern/_src/add.py @@ -5,40 +5,19 @@ import numpy as np import itertools from ...core.parameterization import Parameterized from ...util.caching import Cache_this -from kern import Kern +from kern import CombinationKernel -class Add(Kern): - def __init__(self, subkerns, tensor): - assert all([isinstance(k, Kern) for k in subkerns]) - if tensor: - input_dim = sum([k.input_dim for k in subkerns]) - self.input_slices = [] - n = 0 - for k in subkerns: - self.input_slices.append(slice(n, n+k.input_dim)) - n += k.input_dim - else: - assert all([k.input_dim == subkerns[0].input_dim for k in subkerns]) - input_dim = subkerns[0].input_dim - self.input_slices = [slice(None) for k in subkerns] - super(Add, self).__init__(input_dim, 'add') - self.add_parameters(*subkerns) +class Add(CombinationKernel): + """ + Add given list of kernels together. + propagates gradients thorugh. + """ + def __init__(self, subkerns, name='add'): + super(Add, self).__init__(subkerns, name) - @property - def parts(self): - return self._parameters_ - - def K(self, X, X2=None): - """ - Compute the kernel function. - - :param X: the first set of inputs to the kernel - :param X2: (optional) the second set of arguments to the kernel. If X2 - is None, this is passed throgh to the 'part' object, which - handLes this as X2 == X. - """ + @Cache_this(limit=2, force_kwargs=['which_parts']) + def K(self, X, X2=None, which_parts=None): assert X.shape[1] == self.input_dim - which_parts=None if which_parts is None: which_parts = self.parts elif not isinstance(which_parts, (list, tuple)): @@ -46,12 +25,6 @@ class Add(Kern): which_parts = [which_parts] return sum([p.K(X, X2) for p in which_parts]) - def update_gradients_full(self, dL_dK, X, X2=None): - [p.update_gradients_full(dL_dK, X, X2) for p in self.parts] - - def update_gradients_diag(self, dL_dK, X): - [p.update_gradients_diag(dL_dK, X) for p in self.parts] - def gradients_X(self, dL_dK, X, X2=None): """Compute the gradient of the objective function with respect to X. @@ -67,8 +40,8 @@ class Add(Kern): target[:, p.active_dims] += p.gradients_X(dL_dK, X, X2) return target - def Kdiag(self, X): - which_parts=None + @Cache_this(limit=2, force_kwargs=['which_parts']) + def Kdiag(self, X, which_parts=None): assert X.shape[1] == self.input_dim if which_parts is None: which_parts = self.parts diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index 96bab646..a1106241 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -2,6 +2,7 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) import sys +import numpy as np from ...core.parameterization.parameterized import ParametersChangedMeta, Parameterized from ...util.caching import Cache_this @@ -14,8 +15,11 @@ class KernCallsViaSlicerMeta(ParametersChangedMeta): instance.update_gradients_diag = instance._slice_wrapper(instance.update_gradients_diag, True, True) instance.gradients_X = instance._slice_wrapper(instance.gradients_X, False, True) instance.gradients_X_diag = instance._slice_wrapper(instance.gradients_X_diag, True, True) + instance.psi0 = instance._slice_wrapper(instance.psi0, False, False) + instance.psi1 = instance._slice_wrapper(instance.psi1, False, False) + instance.psi2 = instance._slice_wrapper(instance.psi2, False, False) return instance - + class Kern(Parameterized): __metaclass__ = KernCallsViaSlicerMeta def __init__(self, input_dim, name, *a, **kw): @@ -37,11 +41,11 @@ class Kern(Parameterized): self.input_dim = len(self.active_dims) self._sliced_X = False self._sliced_X2 = False - - @Cache_this(limit=10, ignore_args = (0,)) + + @Cache_this(limit=10)#, ignore_args = (0,)) def _slice_X(self, X): return X[:, self.active_dims] - + def _slice_wrapper(self, operation, diag=False, derivative=False): """ This method wraps the functions in kernel to make sure all kernels allways see their respective input dimension. @@ -56,7 +60,8 @@ class Kern(Parameterized): self._sliced_X = True try: ret = operation(dL_dK, X, *args, **kw) - except: raise + except: + raise finally: self._sliced_X = False return ret @@ -67,7 +72,8 @@ class Kern(Parameterized): self._sliced_X2 = True try: ret = operation(dL_dK, X, X2, *args, **kw) - except: raise + except: + raise finally: self._sliced_X = False self._sliced_X2 = False @@ -79,7 +85,8 @@ class Kern(Parameterized): self._sliced_X = True try: ret = operation(X, *args, **kw) - except: raise + except: + raise finally: self._sliced_X = False return ret @@ -100,10 +107,18 @@ class Kern(Parameterized): +(","+str(bool(diag)) if diag else'') +(','+str(bool(derivative)) if derivative else '') +')') - x_slice_wrapper.__doc__ = "**sliced**\n\n" + (operation.__doc__ or "") + x_slice_wrapper.__doc__ = "**sliced**\n" + (operation.__doc__ or "") return x_slice_wrapper def K(self, X, X2): + """ + Compute the kernel function. + + :param X: the first set of inputs to the kernel + :param X2: (optional) the second set of arguments to the kernel. If X2 + is None, this is passed throgh to the 'part' object, which + handLes this as X2 == X. + """ raise NotImplementedError def Kdiag(self, X): raise NotImplementedError @@ -179,17 +194,10 @@ class Kern(Parameterized): """ Overloading of the '+' operator. for more control, see self.add """ return self.add(other) - def add(self, other, tensor=False): + def add(self, other, name='add'): """ Add another kernel to this one. - If Tensor is False, both kernels are defined on the same _space_. then - the created kernel will have the same number of inputs as self and - other (which must be the same). - - If Tensor is True, then the dimensions are stacked 'horizontally', so - that the resulting kernel has self.input_dim + other.input_dim - :param other: the other kernel to be added :type other: GPy.kern @@ -197,23 +205,23 @@ class Kern(Parameterized): assert isinstance(other, Kern), "only kernels can be added to kernels..." from add import Add kernels = [] - if not tensor and isinstance(self, Add): kernels.extend(self._parameters_) + if isinstance(self, Add): kernels.extend(self._parameters_) else: kernels.append(self) - if not tensor and isinstance(other, Add): kernels.extend(other._parameters_) + if isinstance(other, Add): kernels.extend(other._parameters_) else: kernels.append(other) - return Add(kernels, tensor) + return Add(kernels, name=name) def __mul__(self, other): """ Here we overload the '*' operator. See self.prod for more information""" return self.prod(other) - def __pow__(self, other): - """ - Shortcut for tensor `prod`. - """ - return self.prod(other, tensor=True) + #def __pow__(self, other): + # """ + # Shortcut for tensor `prod`. + # """ + # return self.prod(other, tensor=True) - def prod(self, other, tensor=False, name=None): + def prod(self, other, name=None): """ Multiply two kernels (either on the same space, or on the tensor product of the input space). @@ -226,4 +234,27 @@ class Kern(Parameterized): """ assert isinstance(other, Kern), "only kernels can be added to kernels..." from prod import Prod - return Prod(self, other, tensor, name) + kernels = [] + if isinstance(self, Prod): kernels.extend(self._parameters_) + else: kernels.append(self) + if isinstance(other, Prod): kernels.extend(other._parameters_) + else: kernels.append(other) + return Prod(self, other, name) + + +class CombinationKernel(Kern): + def __init__(self, kernels, name): + assert all([isinstance(k, Kern) for k in kernels]) + input_dim = reduce(np.union1d, (np.r_[x.active_dims] for x in kernels)) + super(CombinationKernel, self).__init__(input_dim, name) + self.add_parameters(*kernels) + + @property + def parts(self): + return self._parameters_ + + def update_gradients_full(self, dL_dK, X, X2=None): + [p.update_gradients_full(dL_dK, X, X2) for p in self.parts] + + def update_gradients_diag(self, dL_dK, X): + [p.update_gradients_diag(dL_dK, X) for p in self.parts] diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py index 51490687..77b2ea51 100644 --- a/GPy/kern/_src/prod.py +++ b/GPy/kern/_src/prod.py @@ -1,10 +1,12 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) -from kern import Kern import numpy as np +from kern import CombinationKernel +from ...util.caching import Cache_this +import itertools -class Prod(Kern): +class Prod(CombinationKernel): """ Computes the product of 2 kernels @@ -15,34 +17,31 @@ class Prod(Kern): :rtype: kernel object """ - def __init__(self, k1, k2, tensor=False,name=None): - if tensor: - name = k1.name + '_xx_' + k2.name if name is None else name - super(Prod, self).__init__(k1.input_dim + k2.input_dim, name) - self.slice1 = slice(0,k1.input_dim) - self.slice2 = slice(k1.input_dim,k1.input_dim+k2.input_dim) - else: - assert k1.input_dim == k2.input_dim, "Error: The input spaces of the kernels to multiply don't have the same dimension." - name = k1.name + '_x_' + k2.name if name is None else name - super(Prod, self).__init__(k1.input_dim, name) - self.slice1 = slice(0, self.input_dim) - self.slice2 = slice(0, self.input_dim) - self.k1 = k1 - self.k2 = k2 - self.add_parameters(self.k1, self.k2) + def __init__(self, kernels, name='prod'): + super(Prod, self).__init__(kernels, name) - def K(self, X, X2=None): - if X2 is None: - return self.k1.K(X[:,self.slice1], None) * self.k2.K(X[:,self.slice2], None) - else: - return self.k1.K(X[:,self.slice1], X2[:,self.slice1]) * self.k2.K(X[:,self.slice2], X2[:,self.slice2]) + @Cache_this(limit=2, force_kwargs=['which_parts']) + def K(self, X, X2=None, which_parts=None): + assert X.shape[1] == self.input_dim + if which_parts is None: + which_parts = self.parts + elif not isinstance(which_parts, (list, tuple)): + # if only one part is given + which_parts = [which_parts] + return reduce(np.multiply, (p.K(X, X2) for p in which_parts)) - def Kdiag(self, X): - return self.k1.Kdiag(X[:,self.slice1]) * self.k2.Kdiag(X[:,self.slice2]) + @Cache_this(limit=2, force_kwargs=['which_parts']) + def Kdiag(self, X, which_parts=None): + assert X.shape[1] == self.input_dim + if which_parts is None: + which_parts = self.parts + return reduce(np.multiply, (p.Kdiag(X) for p in which_parts)) def update_gradients_full(self, dL_dK, X): - self.k1.update_gradients_full(dL_dK*self.k2.K(X[:,self.slice2]), X[:,self.slice1]) - self.k2.update_gradients_full(dL_dK*self.k1.K(X[:,self.slice1]), X[:,self.slice2]) + for k1,k2 in itertools.combinations(self.parts, 2): + k1._sliced_X = k1._sliced_X2 = k2._sliced_X = k2._sliced_X2 = True + k1.update_gradients_full(dL_dK*k2.K(X, X) + self.k2.update_gradients_full(dL_dK*self.k1.K(X[:,self.slice1]), X[:,self.slice2]) def gradients_X(self, dL_dK, X, X2=None): target = np.zeros(X.shape) From 81d35686d987d45df7bbc9ccd1f292d5e419689e Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 11 Mar 2014 10:24:30 +0000 Subject: [PATCH 04/19] slicing tests and ipdb delete --- GPy/testing/kernel_tests.py | 30 +++++++++++++++++++++++++----- GPy/testing/likelihood_tests.py | 6 +++--- 2 files changed, 28 insertions(+), 8 deletions(-) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index d373a546..2789d1de 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -6,7 +6,9 @@ import numpy as np import GPy import sys -verbose = True +verbose = 0 + + class Kern_check_model(GPy.core.Model): """ @@ -91,7 +93,7 @@ class Kern_check_dKdiag_dX(Kern_check_dK_dX): -def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False): +def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verbose=False): """ This function runs on kernels to check the correctness of their implementation. It checks that the covariance function is positive definite @@ -210,7 +212,7 @@ def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False): -class KernelTestsContinuous(unittest.TestCase): +class KernelGradientTestsContinuous(unittest.TestCase): def setUp(self): self.X = np.random.randn(100,2) self.X2 = np.random.randn(110,2) @@ -220,16 +222,34 @@ class KernelTestsContinuous(unittest.TestCase): def test_Matern32(self): k = GPy.kern.Matern32(2) - self.assertTrue(kern_test(k, X=self.X, X2=self.X2, verbose=verbose)) + self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) def test_Matern52(self): k = GPy.kern.Matern52(2) - self.assertTrue(kern_test(k, X=self.X, X2=self.X2, verbose=verbose)) + self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) #TODO: turn off grad checkingwrt X for indexed kernels liek coregionalize +class KernelTestsMiscellaneous(unittest.TestCase): + def setUp(self): + N, D = 100, 10 + self.X = np.linspace(-np.pi, +np.pi, N)[:,None] * np.ones(D) + self.rbf = GPy.kern.RBF(range(2)) + self.linear = GPy.kern.Linear((3,5,6)) + self.matern = GPy.kern.Matern32(np.array([2,4,7])) + self.sumkern = self.rbf + self.linear + self.sumkern += self.matern + self.sumkern.randomize() + + def test_active_dims(self): + self.assertListEqual(self.sumkern.active_dims.tolist(), range(8)) + + def test_which_parts(self): + self.assertTrue(np.allclose(self.sumkern.K(self.X, which_parts=[self.linear, self.matern]), self.linear.K(self.X)+self.matern.K(self.X))) + self.assertTrue(np.allclose(self.sumkern.K(self.X, which_parts=[self.linear, self.rbf]), self.linear.K(self.X)+self.rbf.K(self.X))) + self.assertTrue(np.allclose(self.sumkern.K(self.X, which_parts=self.sumkern.parts[0]), self.rbf.K(self.X))) if __name__ == "__main__": print "Running unit tests, please be (very) patient..." diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py index 631f2ec2..c71842d8 100644 --- a/GPy/testing/likelihood_tests.py +++ b/GPy/testing/likelihood_tests.py @@ -651,7 +651,7 @@ class LaplaceTests(unittest.TestCase): m2['.*white'].constrain_fixed(1e-6) m2['.*rbf.variance'].constrain_bounded(1e-4, 10) m2.randomize() - + if debug: print m1 print m2 @@ -663,7 +663,7 @@ class LaplaceTests(unittest.TestCase): if debug: print m1 print m2 - + m2[:] = m1[:] #Predict for training points to get posterior mean and variance @@ -702,7 +702,7 @@ class LaplaceTests(unittest.TestCase): m1.randomize() import ipdb;ipdb.set_trace() m2[:] = m1[:] - + np.testing.assert_almost_equal(m1.log_likelihood(), m2.log_likelihood(), decimal=2) #Check they are checkgradding From 3e91ea497d1214bd1ee04612962e6809e5d4814c Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 11 Mar 2014 10:24:51 +0000 Subject: [PATCH 05/19] caching doc --- GPy/util/caching.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 0b6f7234..5de03059 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -109,5 +109,5 @@ class Cache_this(object): self.c = Cacher(f, self.limit, ignore_args=self.ignore_args, force_kwargs=self.force_args) return self.c(*args, **kw) f_wrap._cacher = self - f_wrap.__doc__ = "**cached**\n\n" + (f.__doc__ or "") + f_wrap.__doc__ = "**cached**" + (f.__doc__ or "") return f_wrap From 10608a45656ad61aa34ecd3197c716c11640cb67 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 11 Mar 2014 10:25:21 +0000 Subject: [PATCH 06/19] empty spaces --- GPy/models/sparse_gp_regression.py | 4 ++-- GPy/plotting/matplot_dep/models_plots.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/GPy/models/sparse_gp_regression.py b/GPy/models/sparse_gp_regression.py index 99176601..7edb93e4 100644 --- a/GPy/models/sparse_gp_regression.py +++ b/GPy/models/sparse_gp_regression.py @@ -45,10 +45,10 @@ class SparseGPRegression(SparseGP): assert Z.shape[1] == input_dim likelihood = likelihoods.Gaussian() - + if not (X_variance is None): X = NormalPosterior(X,X_variance) - + SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method=VarDTC()) def _getstate(self): diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index 4ca4441e..86777527 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -56,7 +56,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', if ax is None: fig = pb.figure(num=fignum) ax = fig.add_subplot(111) - + if hasattr(model, 'has_uncertain_inputs') and model.has_uncertain_inputs(): X = model.X.mean X_variance = param_to_array(model.X.variance) From 85a471e0f6340dadbd4fe9002ee7d82b5dc07ef0 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 11 Mar 2014 16:22:45 +0000 Subject: [PATCH 07/19] oh huge bug in checkgrad global --- GPy/core/model.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index a858a62d..6a6fe1ba 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -253,7 +253,7 @@ class Model(Parameterized): sgd.run() self.optimization_runs.append(sgd) - def _checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3): + def _checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3, _debug=False): """ Check the gradient of the ,odel by comparing to a numerical estimate. If the verbose flag is passed, invividual @@ -271,7 +271,7 @@ class Model(Parameterized): and numerical gradients is within of unity. """ x = self._get_params_transformed().copy() - + if not verbose: # make sure only to test the selected parameters if target_param is None: @@ -298,12 +298,12 @@ class Model(Parameterized): dx = dx[transformed_index] gradient = gradient[transformed_index] - + denominator = (2 * np.dot(dx, gradient)) global_ratio = (f1 - f2) / np.where(denominator==0., 1e-32, denominator) gloabl_diff = (f1 - f2) - denominator - - return (np.abs(1. - global_ratio) < tolerance) or (np.abs(gloabl_diff) < tolerance) + + return (np.abs(1. - global_ratio) < tolerance) or (np.abs(gloabl_diff) == 0) else: # check the gradient of each parameter individually, and do some pretty printing try: @@ -349,6 +349,8 @@ class Model(Parameterized): xx[xind] -= 2.*step f2 = self.objective_function(xx) numerical_gradient = (f1 - f2) / (2 * step) + if _debug: + self.gradient[xind] = numerical_gradient if np.all(gradient[xind]==0): ratio = (f1-f2) == gradient[xind] else: ratio = (f1 - f2) / (2 * step * gradient[xind]) difference = np.abs((f1 - f2) / 2 / step - gradient[xind]) @@ -366,7 +368,7 @@ class Model(Parameterized): ng = '%.6f' % float(numerical_gradient) grad_string = "{0:<{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name, r, d, g, ng, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4]) print grad_string - + self._set_params_transformed(x) return ret From 74999a89ad37bc55821fc12ce786400acc5f722f Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 11 Mar 2014 16:23:29 +0000 Subject: [PATCH 08/19] gradient check --- GPy/core/parameterization/param.py | 4 +-- GPy/core/parameterization/parameter_core.py | 39 +++++++++++---------- 2 files changed, 22 insertions(+), 21 deletions(-) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index a2dc9514..8eb10608 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -446,8 +446,8 @@ class ParamConcatenation(object): def untie(self, *ties): [param.untie(*ties) for param in self.params] - def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3): - return self.params[0]._highest_parent_._checkgrad(self, verbose, step, tolerance) + def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3, _debug=False): + return self.params[0]._highest_parent_._checkgrad(self, verbose, step, tolerance, _debug=_debug) #checkgrad.__doc__ = Gradcheckable.checkgrad.__doc__ __lt__ = lambda self, val: self._vals() < val diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 38fe0526..5727bc17 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__ = '2013-12-16' +__updated__ = '2014-03-11' class HierarchyError(Exception): """ @@ -34,7 +34,7 @@ def adjust_name_for_printing(name): class Observable(object): """ Observable pattern for parameterization. - + This Object allows for observers to register with self and a (bound!) function as an observer. Every time the observable changes, it sends a notification with self as only argument to all its observers. @@ -43,10 +43,10 @@ class Observable(object): def __init__(self, *args, **kwargs): super(Observable, self).__init__(*args, **kwargs) self._observer_callables_ = [] - + def add_observer(self, observer, callble, priority=0): self._insert_sorted(priority, observer, callble) - + def remove_observer(self, observer, callble=None): to_remove = [] for p, obs, clble in self._observer_callables_: @@ -58,15 +58,15 @@ class Observable(object): to_remove.append((p, obs, clble)) for r in to_remove: self._observer_callables_.remove(r) - + def notify_observers(self, which=None, min_priority=None): """ Notifies all observers. Which is the element, which kicked off this notification loop. - + NOTE: notifies only observers with priority p > min_priority! ^^^^^^^^^^^^^^^^ - + :param which: object, which started this notification loop :param min_priority: only notify observers with priority > min_priority if min_priority is None, notify all observers in order @@ -88,11 +88,11 @@ class Observable(object): break ins += 1 self._observer_callables_.insert(ins, (p, o, c)) - + class Pickleable(object): """ Make an object pickleable (See python doc 'pickling'). - + This class allows for pickling support by Memento pattern. _getstate returns a memento of the class, which gets pickled. _setstate() (re-)sets the state of the class to the memento @@ -153,7 +153,7 @@ class Pickleable(object): class Parentable(object): """ Enable an Object to have a parent. - + Additionally this adds the parent_index, which is the index for the parent to look for in its parameter list. """ @@ -161,7 +161,7 @@ class Parentable(object): _parent_index_ = None def __init__(self, *args, **kwargs): super(Parentable, self).__init__(*args, **kwargs) - + def has_parent(self): """ Return whether this parentable object currently has a parent. @@ -205,8 +205,8 @@ class Gradcheckable(Parentable): """ def __init__(self, *a, **kw): super(Gradcheckable, self).__init__(*a, **kw) - - def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3): + + def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3, _debug=False): """ Check the gradient of this parameter with respect to the highest parent's objective function. @@ -214,20 +214,21 @@ class Gradcheckable(Parentable): with a stepsize step. The check passes if either the ratio or the difference between numerical and analytical gradient is smaller then tolerance. - + :param bool verbose: whether each parameter shall be checked individually. :param float step: the stepsize for the numerical three point gradient estimate. :param flaot tolerance: the tolerance for the gradient ratio or difference. """ if self.has_parent(): - return self._highest_parent_._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance) - return self._checkgrad(self[''], verbose=verbose, step=step, tolerance=tolerance) - def _checkgrad(self, param): + return self._highest_parent_._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance, _debug=_debug) + return self._checkgrad(self[''], verbose=verbose, step=step, tolerance=tolerance, _debug=_debug) + + def _checkgrad(self, param, verbose=0, step=1e-6, tolerance=1e-3, _debug=False): """ Perform the checkgrad on the model. TODO: this can be done more efficiently, when doing it inside here """ - raise NotImplementedError, "Need log likelihood to check gradient against" + raise HierarchyError, "This parameter is not in a model with a likelihood, and, therefore, cannot be gradient checked!" class Nameable(Gradcheckable): @@ -258,7 +259,7 @@ class Nameable(Gradcheckable): def hierarchy_name(self, adjust_for_printing=True): """ return the name for this object with the parents names attached by dots. - + :param bool adjust_for_printing: whether to call :func:`~adjust_for_printing()` on the names, recursively """ From e078bb47e10f97519805f363416eae0281ab6c20 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 11 Mar 2014 16:23:51 +0000 Subject: [PATCH 09/19] psi_stat_expectaions now working with new parameterized --- GPy/testing/psi_stat_expectation_tests.py | 57 ++++++----------------- 1 file changed, 13 insertions(+), 44 deletions(-) diff --git a/GPy/testing/psi_stat_expectation_tests.py b/GPy/testing/psi_stat_expectation_tests.py index aec0d36d..075800f6 100644 --- a/GPy/testing/psi_stat_expectation_tests.py +++ b/GPy/testing/psi_stat_expectation_tests.py @@ -12,6 +12,7 @@ import numpy from GPy.kern import RBF from GPy.kern import Linear from copy import deepcopy +from GPy.core.parameterization.variational import NormalPosterior __test__ = lambda: 'deep' in sys.argv # np.random.seed(0) @@ -28,53 +29,20 @@ def ard(p): class Test(unittest.TestCase): input_dim = 9 num_inducing = 13 - N = 300 + N = 1000 Nsamples = 1e6 def setUp(self): - i_s_dim_list = [2,4,3] - indices = numpy.cumsum(i_s_dim_list).tolist() - input_slices = [slice(a,b) for a,b in zip([None]+indices, indices)] - #input_slices[2] = deepcopy(input_slices[1]) - input_slice_kern = GPy.kern.kern(9, - [ - RBF(i_s_dim_list[0], np.random.rand(), np.random.rand(i_s_dim_list[0]), ARD=True), - RBF(i_s_dim_list[1], np.random.rand(), np.random.rand(i_s_dim_list[1]), ARD=True), - Linear(i_s_dim_list[2], np.random.rand(i_s_dim_list[2]), ARD=True) - ], - input_slices = input_slices - ) self.kerns = ( -# input_slice_kern, -# (GPy.kern.rbf(self.input_dim, ARD=True) + -# GPy.kern.linear(self.input_dim, ARD=True) + -# GPy.kern.bias(self.input_dim) + -# GPy.kern.white(self.input_dim)), - (#GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True) - GPy.kern.Linear(self.input_dim, np.random.rand(self.input_dim), ARD=True) - +GPy.kern.RBF(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True) -# +GPy.kern.bias(self.input_dim) -# +GPy.kern.white(self.input_dim)), - ), -# (GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True) + -# GPy.kern.bias(self.input_dim, np.random.rand())), -# (GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True) -# +GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True) -# #+GPy.kern.bias(self.input_dim, np.random.rand()) -# #+GPy.kern.white(self.input_dim, np.random.rand())), -# ), -# GPy.kern.white(self.input_dim, np.random.rand())), -# GPy.kern.rbf(self.input_dim), GPy.kern.rbf(self.input_dim, ARD=True), -# GPy.kern.linear(self.input_dim, ARD=False), GPy.kern.linear(self.input_dim, ARD=True), -# GPy.kern.linear(self.input_dim) + GPy.kern.bias(self.input_dim), -# GPy.kern.rbf(self.input_dim) + GPy.kern.bias(self.input_dim), -# GPy.kern.linear(self.input_dim) + GPy.kern.bias(self.input_dim) + GPy.kern.white(self.input_dim), -# GPy.kern.rbf(self.input_dim) + GPy.kern.bias(self.input_dim) + GPy.kern.white(self.input_dim), -# GPy.kern.bias(self.input_dim), GPy.kern.white(self.input_dim), + GPy.kern.RBF(self.input_dim, ARD=True)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim), + GPy.kern.RBF(self.input_dim)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim), + GPy.kern.Linear(self.input_dim) + GPy.kern.Bias(self.input_dim) + GPy.kern.White(self.input_dim), + GPy.kern.Linear(self.input_dim, ARD=True) + GPy.kern.Bias(self.input_dim) + GPy.kern.White(self.input_dim), ) - self.q_x_mean = np.random.randn(self.input_dim) - self.q_x_variance = np.exp(np.random.randn(self.input_dim)) + self.q_x_mean = np.random.randn(self.input_dim)[None] + self.q_x_variance = np.exp(.5*np.random.randn(self.input_dim))[None] self.q_x_samples = np.random.randn(self.Nsamples, self.input_dim) * np.sqrt(self.q_x_variance) + self.q_x_mean + self.q_x = NormalPosterior(self.q_x_mean, self.q_x_variance) self.Z = np.random.randn(self.num_inducing, self.input_dim) self.q_x_mean.shape = (1, self.input_dim) self.q_x_variance.shape = (1, self.input_dim) @@ -114,8 +82,9 @@ class Test(unittest.TestCase): def test_psi2(self): for kern in self.kerns: + kern.randomize() Nsamples = int(np.floor(self.Nsamples/self.N)) - psi2 = kern.psi2(self.Z, self.q_x_mean, self.q_x_variance) + psi2 = kern.psi2(self.Z, self.q_x) K_ = np.zeros((self.num_inducing, self.num_inducing)) diffs = [] for i, q_x_sample_stripe in enumerate(np.array_split(self.q_x_samples, self.Nsamples / Nsamples)): @@ -130,8 +99,8 @@ class Test(unittest.TestCase): pylab.figure(msg) pylab.plot(diffs, marker='x', mew=.2) # print msg, np.allclose(psi2.squeeze(), K_, rtol=1e-1, atol=.1) - self.assertTrue(np.allclose(psi2.squeeze(), K_), - #rtol=1e-1, atol=.1), + self.assertTrue(np.allclose(psi2.squeeze(), K_, + atol=.1, rtol=1), msg=msg + ": not matching") # sys.stdout.write(".") except: From 01f5d789c5999de7df8818ce659c2c0ea0a633fb Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 11 Mar 2014 16:24:09 +0000 Subject: [PATCH 10/19] automatic slicing --- GPy/kern/_src/add.py | 78 ++++++++++++------------------ GPy/kern/_src/kern.py | 107 +++++++++++------------------------------- GPy/kern/_src/rbf.py | 31 ++++++------ 3 files changed, 72 insertions(+), 144 deletions(-) diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py index 433a8921..8a3cefaf 100644 --- a/GPy/kern/_src/add.py +++ b/GPy/kern/_src/add.py @@ -23,7 +23,7 @@ class Add(CombinationKernel): elif not isinstance(which_parts, (list, tuple)): # if only one part is given which_parts = [which_parts] - return sum([p.K(X, X2) for p in which_parts]) + return reduce(np.add, (p.K(X, X2) for p in which_parts)) def gradients_X(self, dL_dK, X, X2=None): """Compute the gradient of the objective function with respect to X. @@ -49,14 +49,14 @@ class Add(CombinationKernel): def psi0(self, Z, variational_posterior): - return np.sum([p.psi0(Z[:, i_s], variational_posterior[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)],0) + return reduce(np.add, (p.psi0(Z, variational_posterior) for p in self.parts)) def psi1(self, Z, variational_posterior): - return np.sum([p.psi1(Z[:, i_s], variational_posterior[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)], 0) + return reduce(np.add, (p.psi1(Z, variational_posterior) for p in self.parts)) def psi2(self, Z, variational_posterior): - psi2 = np.sum([p.psi2(Z[:, i_s], variational_posterior[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)], 0) - + psi2 = reduce(np.add, (p.psi2(Z, variational_posterior) for p in self.parts)) + return psi2 # compute the "cross" terms from static import White, Bias from rbf import RBF @@ -64,18 +64,23 @@ class Add(CombinationKernel): from linear import Linear #ffrom fixed import Fixed - for (p1, i1), (p2, i2) in itertools.combinations(itertools.izip(self._parameters_, self.input_slices), 2): + for p1, p2 in itertools.combinations(self.parts, 2): + i1, i2 = p1.active_dims, p2.active_dims # white doesn;t combine with anything if isinstance(p1, White) or isinstance(p2, White): pass # rbf X bias #elif isinstance(p1, (Bias, Fixed)) and isinstance(p2, (RBF, RBFInv)): elif isinstance(p1, Bias) and isinstance(p2, (RBF, Linear)): - tmp = p2.psi1(Z[:,i2], variational_posterior[:, i_s]) + # manual override for slicing: + p2._sliced_X = p1._sliced_X = True + tmp = p2.psi1(Z[:,i2], variational_posterior[:, i1]) psi2 += p1.variance * (tmp[:, :, None] + tmp[:, None, :]) #elif isinstance(p2, (Bias, Fixed)) and isinstance(p1, (RBF, RBFInv)): elif isinstance(p2, Bias) and isinstance(p1, (RBF, Linear)): - tmp = p1.psi1(Z[:,i1], variational_posterior[:, i_s]) + # manual override for slicing: + p2._sliced_X = p1._sliced_X = True + tmp = p1.psi1(Z[:,i1], variational_posterior[:, i2]) psi2 += p2.variance * (tmp[:, :, None] + tmp[:, None, :]) else: raise NotImplementedError, "psi2 cannot be computed for this kernel" @@ -83,11 +88,10 @@ class Add(CombinationKernel): def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): from static import White, Bias - for p1, is1 in zip(self._parameters_, self.input_slices): - + for p1 in self.parts: #compute the effective dL_dpsi1. Extra terms appear becaue of the cross terms in psi2! eff_dL_dpsi1 = dL_dpsi1.copy() - for p2, is2 in zip(self._parameters_, self.input_slices): + for p2 in self.parts: if p2 is p1: continue if isinstance(p2, White): @@ -95,42 +99,35 @@ class Add(CombinationKernel): elif isinstance(p2, Bias): eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2. else: - eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z[:,is2], variational_posterior[:, is1]) * 2. - - - p1.update_gradients_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z[:,is1], variational_posterior[:, is1]) - + eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z, variational_posterior) * 2. + p1.update_gradients_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior) def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): from static import White, Bias target = np.zeros(Z.shape) - for p1, is1 in zip(self._parameters_, self.input_slices): - + for p1 in self.parts: #compute the effective dL_dpsi1. extra terms appear becaue of the cross terms in psi2! eff_dL_dpsi1 = dL_dpsi1.copy() - for p2, is2 in zip(self._parameters_, self.input_slices): + for p2 in self.parts: if p2 is p1: continue if isinstance(p2, White): continue elif isinstance(p2, Bias): - eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2. + eff_dL_dpsi1 += 0#dL_dpsi2.sum(1) * p2.variance * 2. else: - eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z[:,is2], variational_posterior[:, is2]) * 2. - - - target += p1.gradients_Z_expectations(eff_dL_dpsi1, dL_dpsi2, Z[:,is1], variational_posterior[:, is1]) + eff_dL_dpsi1 += 0#dL_dpsi2.sum(1) * p2.psi1(Z, variational_posterior) * 2. + target[:, p1.active_dims] += p1.gradients_Z_expectations(eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior) return target def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): from static import White, Bias target_mu = np.zeros(variational_posterior.shape) target_S = np.zeros(variational_posterior.shape) - for p1, is1 in zip(self._parameters_, self.input_slices): - + for p1 in self._parameters_: #compute the effective dL_dpsi1. extra terms appear becaue of the cross terms in psi2! eff_dL_dpsi1 = dL_dpsi1.copy() - for p2, is2 in zip(self._parameters_, self.input_slices): + for p2 in self._parameters_: if p2 is p1: continue if isinstance(p2, White): @@ -138,35 +135,20 @@ class Add(CombinationKernel): elif isinstance(p2, Bias): eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2. else: - eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z[:,is2], variational_posterior[:, is2]) * 2. - - - a, b = p1.gradients_qX_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z[:,is1], variational_posterior[:, is1]) - target_mu += a - target_S += b + eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z, variational_posterior) * 2. + a, b = p1.gradients_qX_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior) + target_mu[:, p1.active_dims] += a + target_S[:, p1.active_dims] += b return target_mu, target_S - def input_sensitivity(self): - in_sen = np.zeros((self.num_params, self.input_dim)) - for i, [p, i_s] in enumerate(zip(self._parameters_, self.input_slices)): - in_sen[i, i_s] = p.input_sensitivity() - return in_sen - def _getstate(self): """ Get the current state of the class, here just all the indices, rest can get recomputed """ - return Parameterized._getstate(self) + [#self._parameters_, - self.input_dim, - self.input_slices, - self._param_slices_ - ] + return super(Add, self)._getstate() def _setstate(self, state): - self._param_slices_ = state.pop() - self.input_slices = state.pop() - self.input_dim = state.pop() - Parameterized._setstate(self, state) + super(Add, self)._setstate(state) diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index a1106241..8feb9a04 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -3,25 +3,18 @@ import sys import numpy as np -from ...core.parameterization.parameterized import ParametersChangedMeta, Parameterized +from ...core.parameterization.parameterized import Parameterized +from kernel_slice_operations import KernCallsViaSlicerMeta from ...util.caching import Cache_this -class KernCallsViaSlicerMeta(ParametersChangedMeta): - def __call__(self, *args, **kw): - instance = super(KernCallsViaSlicerMeta, self).__call__(*args, **kw) - instance.K = instance._slice_wrapper(instance.K) - instance.Kdiag = instance._slice_wrapper(instance.Kdiag, True) - instance.update_gradients_full = instance._slice_wrapper(instance.update_gradients_full, False, True) - instance.update_gradients_diag = instance._slice_wrapper(instance.update_gradients_diag, True, True) - instance.gradients_X = instance._slice_wrapper(instance.gradients_X, False, True) - instance.gradients_X_diag = instance._slice_wrapper(instance.gradients_X_diag, True, True) - instance.psi0 = instance._slice_wrapper(instance.psi0, False, False) - instance.psi1 = instance._slice_wrapper(instance.psi1, False, False) - instance.psi2 = instance._slice_wrapper(instance.psi2, False, False) - return instance + class Kern(Parameterized): + #=========================================================================== + # This adds input slice support. The rather ugly code for slicing can be + # found in kernel_slice_operations __metaclass__ = KernCallsViaSlicerMeta + #=========================================================================== def __init__(self, input_dim, name, *a, **kw): """ The base class for a kernel: a positive definite function @@ -40,76 +33,11 @@ class Kern(Parameterized): self.active_dims = input_dim self.input_dim = len(self.active_dims) self._sliced_X = False - self._sliced_X2 = False @Cache_this(limit=10)#, ignore_args = (0,)) def _slice_X(self, X): return X[:, self.active_dims] - def _slice_wrapper(self, operation, diag=False, derivative=False): - """ - This method wraps the functions in kernel to make sure all kernels allways see their respective input dimension. - The different switches are: - diag: if X2 exists - derivative: if firest arg is dL_dK - """ - if derivative: - if diag: - def x_slice_wrapper(dL_dK, X, *args, **kw): - X = self._slice_X(X) if not self._sliced_X else X - self._sliced_X = True - try: - ret = operation(dL_dK, X, *args, **kw) - except: - raise - finally: - self._sliced_X = False - return ret - else: - def x_slice_wrapper(dL_dK, X, X2=None, *args, **kw): - X, X2 = self._slice_X(X) if not self._sliced_X else X, self._slice_X(X2) if X2 is not None and not self._sliced_X2 else X2 - self._sliced_X = True - self._sliced_X2 = True - try: - ret = operation(dL_dK, X, X2, *args, **kw) - except: - raise - finally: - self._sliced_X = False - self._sliced_X2 = False - return ret - else: - if diag: - def x_slice_wrapper(X, *args, **kw): - X = self._slice_X(X) if not self._sliced_X else X - self._sliced_X = True - try: - ret = operation(X, *args, **kw) - except: - raise - finally: - self._sliced_X = False - return ret - else: - def x_slice_wrapper(X, X2=None, *args, **kw): - X, X2 = self._slice_X(X) if not self._sliced_X else X, self._slice_X(X2) if X2 is not None and not self._sliced_X2 else X2 - self._sliced_X = True - self._sliced_X2 = True - try: - ret = operation(X, X2, *args, **kw) - except: raise - finally: - self._sliced_X = False - self._sliced_X2 = False - return ret - x_slice_wrapper._operation = operation - x_slice_wrapper.__name__ = ("slicer("+operation.__name__ - +(","+str(bool(diag)) if diag else'') - +(','+str(bool(derivative)) if derivative else '') - +')') - x_slice_wrapper.__doc__ = "**sliced**\n" + (operation.__doc__ or "") - return x_slice_wrapper - def K(self, X, X2): """ Compute the kernel function. @@ -241,6 +169,21 @@ class Kern(Parameterized): else: kernels.append(other) return Prod(self, other, name) + def _getstate(self): + """ + Get the current state of the class, + here just all the indices, rest can get recomputed + """ + return super(Kern, self)._getstate() + [ + self.active_dims, + self.input_dim, + self._sliced_X] + + def _setstate(self, state): + self._sliced_X = state.pop() + self.input_dim = state.pop() + self.active_dims = state.pop() + super(Kern, self)._setstate(state) class CombinationKernel(Kern): def __init__(self, kernels, name): @@ -258,3 +201,9 @@ class CombinationKernel(Kern): def update_gradients_diag(self, dL_dK, X): [p.update_gradients_diag(dL_dK, X) for p in self.parts] + + def input_sensitivity(self): + in_sen = np.zeros((self.num_params, self.input_dim)) + for i, p in enumerate(self.parts): + in_sen[i, p.active_dims] = p.input_sensitivity() + return in_sen diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index cd6c41e9..7ba1f35d 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -56,28 +56,28 @@ class RBF(Stationary): if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): _, _dpsi1_dvariance, _, _, _, _, _dpsi1_dlengthscale = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) _, _dpsi2_dvariance, _, _, _, _, _dpsi2_dlengthscale = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) - + #contributions from psi0: self.variance.gradient = np.sum(dL_dpsi0) - + #from psi1 self.variance.gradient += np.sum(dL_dpsi1 * _dpsi1_dvariance) if self.ARD: self.lengthscale.gradient = (dL_dpsi1[:,:,None]*_dpsi1_dlengthscale).reshape(-1,self.input_dim).sum(axis=0) else: self.lengthscale.gradient = (dL_dpsi1[:,:,None]*_dpsi1_dlengthscale).sum() - - + #from psi2 self.variance.gradient += (dL_dpsi2 * _dpsi2_dvariance).sum() if self.ARD: self.lengthscale.gradient += (dL_dpsi2[:,:,:,None] * _dpsi2_dlengthscale).reshape(-1,self.input_dim).sum(axis=0) else: self.lengthscale.gradient += (dL_dpsi2[:,:,:,None] * _dpsi2_dlengthscale).sum() - + elif isinstance(variational_posterior, variational.NormalPosterior): - - l2 = self.lengthscale **2 + l2 = self.lengthscale**2 + if l2.size != self.input_dim: + l2 = l2*np.ones(self.input_dim) #contributions from psi0: self.variance.gradient = np.sum(dL_dpsi0) @@ -92,11 +92,9 @@ class RBF(Stationary): else: self.lengthscale.gradient += dpsi1_dlength.sum() self.variance.gradient += np.sum(dL_dpsi1 * psi1) / self.variance - #from psi2 S = variational_posterior.variance _, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - if not self.ARD: self.lengthscale.gradient += self._weave_psi2_lengthscale_grads(dL_dpsi2, psi2, Zdist_sq, S, mudist_sq, l2).sum() else: @@ -112,17 +110,16 @@ class RBF(Stationary): if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): _, _, _, _, _, _dpsi1_dZ, _ = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) _, _, _, _, _, _dpsi2_dZ, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) - + #psi1 grad = (dL_dpsi1[:, :, None] * _dpsi1_dZ).sum(axis=0) - + #psi2 grad += (dL_dpsi2[:, :, :, None] * _dpsi2_dZ).sum(axis=0).sum(axis=1) - + return grad elif isinstance(variational_posterior, variational.NormalPosterior): - l2 = self.lengthscale **2 #psi1 @@ -145,10 +142,10 @@ class RBF(Stationary): # Spike-and-Slab GPLVM if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): ndata = variational_posterior.mean.shape[0] - + _, _, _dpsi1_dgamma, _dpsi1_dmu, _dpsi1_dS, _, _ = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) _, _, _dpsi2_dgamma, _dpsi2_dmu, _dpsi2_dS, _, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) - + #psi1 grad_mu = (dL_dpsi1[:, :, None] * _dpsi1_dmu).sum(axis=1) grad_S = (dL_dpsi1[:, :, None] * _dpsi1_dS).sum(axis=1) @@ -157,11 +154,11 @@ class RBF(Stationary): grad_mu += (dL_dpsi2[:, :, :, None] * _dpsi2_dmu).reshape(ndata,-1,self.input_dim).sum(axis=1) grad_S += (dL_dpsi2[:, :, :, None] * _dpsi2_dS).reshape(ndata,-1,self.input_dim).sum(axis=1) grad_gamma += (dL_dpsi2[:,:,:, None] * _dpsi2_dgamma).reshape(ndata,-1,self.input_dim).sum(axis=1) - + return grad_mu, grad_S, grad_gamma elif isinstance(variational_posterior, variational.NormalPosterior): - + l2 = self.lengthscale **2 #psi1 denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) From 5f2b383510f31d592acf1601970caec173f38530 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Mar 2014 09:51:26 +0000 Subject: [PATCH 11/19] plotting returns --- GPy/plotting/matplot_dep/models_plots.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index 86777527..8b5a40e7 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -68,7 +68,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #work out what the inputs are for plotting (1D or 2D) fixed_dims = np.array([i for i,v in fixed_inputs]) free_dims = np.setdiff1d(np.arange(model.input_dim),fixed_dims) - + plots = {} #one dimensional plotting if len(free_dims) == 1: @@ -89,20 +89,20 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', m, v, lower, upper = model.predict(Xgrid) Y = Y for d in which_data_ycols: - gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], ax=ax, edgecol=linecol, fillcol=fillcol) - ax.plot(X[which_data_rows,free_dims], Y[which_data_rows, d], 'kx', mew=1.5) + plots['gpplot'] = gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], ax=ax, edgecol=linecol, fillcol=fillcol) + plots['dataplot'] = ax.plot(X[which_data_rows,free_dims], Y[which_data_rows, d], 'kx', mew=1.5) #optionally plot some samples if samples: #NOTE not tested with fixed_inputs Ysim = model.posterior_samples(Xgrid, samples) for yi in Ysim.T: - ax.plot(Xnew, yi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25) + plots['posterior_samples'] = ax.plot(Xnew, yi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25) #ax.plot(Xnew, yi[:,None], marker='x', linestyle='--',color=Tango.colorsHex['darkBlue']) #TODO apply this line for discrete outputs. #add error bars for uncertain (if input uncertainty is being modelled) if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs(): - ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(), + plots['xerrorbar'] = ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(), xerr=2 * np.sqrt(X_variance[which_data_rows, free_dims].flatten()), ecolor='k', fmt=None, elinewidth=.5, alpha=.5) @@ -118,7 +118,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims] Zu = Z[:,free_dims] z_height = ax.get_ylim()[0] - ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12) + plots['inducing_inputs'] = ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12) @@ -143,8 +143,8 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', Y = Y for d in which_data_ycols: m_d = m[:,d].reshape(resolution, resolution).T - ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) - ax.scatter(X[which_data_rows, free_dims[0]], X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.) + plots['contour'] = ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) + plots['dataplot'] = ax.scatter(X[which_data_rows, free_dims[0]], X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.) #set the limits of the plot to some sensible values ax.set_xlim(xmin[0], xmax[0]) @@ -157,11 +157,11 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', if hasattr(model,"Z"): #Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims] Zu = Z[:,free_dims] - ax.plot(Zu[:,free_dims[0]], Zu[:,free_dims[1]], 'wo') + plots['inducing_inputs'] = ax.plot(Zu[:,free_dims[0]], Zu[:,free_dims[1]], 'wo') else: raise NotImplementedError, "Cannot define a frame with more than two input dimensions" - + return plots def plot_fit_f(model, *args, **kwargs): """ From 02bce95c41606da02ea8f9548282425fe125fbdf Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Mar 2014 09:51:46 +0000 Subject: [PATCH 12/19] psi stat testing improvements, gradients not working yet --- GPy/testing/psi_stat_expectation_tests.py | 2 +- GPy/testing/psi_stat_gradient_tests.py | 12 +++++++----- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/GPy/testing/psi_stat_expectation_tests.py b/GPy/testing/psi_stat_expectation_tests.py index 075800f6..04167bef 100644 --- a/GPy/testing/psi_stat_expectation_tests.py +++ b/GPy/testing/psi_stat_expectation_tests.py @@ -34,7 +34,7 @@ class Test(unittest.TestCase): def setUp(self): self.kerns = ( - GPy.kern.RBF(self.input_dim, ARD=True)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim), + GPy.kern.RBF([0,1,2], ARD=True)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim), GPy.kern.RBF(self.input_dim)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim), GPy.kern.Linear(self.input_dim) + GPy.kern.Bias(self.input_dim) + GPy.kern.White(self.input_dim), GPy.kern.Linear(self.input_dim, ARD=True) + GPy.kern.Bias(self.input_dim) + GPy.kern.White(self.input_dim), diff --git a/GPy/testing/psi_stat_gradient_tests.py b/GPy/testing/psi_stat_gradient_tests.py index fc189f93..d51cd913 100644 --- a/GPy/testing/psi_stat_gradient_tests.py +++ b/GPy/testing/psi_stat_gradient_tests.py @@ -11,6 +11,7 @@ import itertools from GPy.core import Model from GPy.core.parameterization.param import Param from GPy.core.parameterization.transformations import Logexp +from GPy.core.parameterization.variational import NormalPosterior class PsiStatModel(Model): def __init__(self, which, X, X_variance, Z, num_inducing, kernel): @@ -18,23 +19,24 @@ class PsiStatModel(Model): self.which = which self.X = Param("X", X) self.X_variance = Param('X_variance', X_variance, Logexp()) + self.q = NormalPosterior(self.X, self.X_variance) self.Z = Param("Z", Z) self.N, self.input_dim = X.shape self.num_inducing, input_dim = Z.shape assert self.input_dim == input_dim, "shape missmatch: Z:{!s} X:{!s}".format(Z.shape, X.shape) self.kern = kernel - self.psi_ = self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance) - self.add_parameters(self.X, self.X_variance, self.Z, self.kern) + self.psi_ = self.kern.__getattribute__(self.which)(self.Z, self.q) + self.add_parameters(self.q, self.Z, self.kern) def log_likelihood(self): return self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance).sum() def parameters_changed(self): - psimu, psiS = self.kern.__getattribute__("d" + self.which + "_dmuS")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance) + psimu, psiS = self.kern.__getattribute__("d" + self.which + "_dmuS")(numpy.ones_like(self.psi_), self.Z, self.q) self.X.gradient = psimu self.X_variance.gradient = psiS #psimu, psiS = numpy.ones(self.N * self.input_dim), numpy.ones(self.N * self.input_dim) - try: psiZ = self.kern.__getattribute__("d" + self.which + "_dZ")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance) + try: psiZ = self.kern.__getattribute__("d" + self.which + "_dZ")(numpy.ones_like(self.psi_), self.Z, self.q) except AttributeError: psiZ = numpy.zeros_like(self.Z) self.Z.gradient = psiZ #psiZ = numpy.ones(self.num_inducing * self.input_dim) @@ -176,6 +178,6 @@ if __name__ == "__main__": +GPy.kern.White(input_dim) ) ) - m2.ensure_default_constraints() + #m2.ensure_default_constraints() else: unittest.main() From dfb63860ca9f6b8b10fc4879a21a25733e0277c1 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Mar 2014 12:03:25 +0000 Subject: [PATCH 13/19] psi stat expectations with slices --- GPy/testing/psi_stat_expectation_tests.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/GPy/testing/psi_stat_expectation_tests.py b/GPy/testing/psi_stat_expectation_tests.py index 04167bef..ffbde37c 100644 --- a/GPy/testing/psi_stat_expectation_tests.py +++ b/GPy/testing/psi_stat_expectation_tests.py @@ -34,10 +34,11 @@ class Test(unittest.TestCase): def setUp(self): self.kerns = ( - GPy.kern.RBF([0,1,2], ARD=True)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim), - GPy.kern.RBF(self.input_dim)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim), - GPy.kern.Linear(self.input_dim) + GPy.kern.Bias(self.input_dim) + GPy.kern.White(self.input_dim), - GPy.kern.Linear(self.input_dim, ARD=True) + GPy.kern.Bias(self.input_dim) + GPy.kern.White(self.input_dim), + #GPy.kern.RBF([0,1,2], ARD=True)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim), + #GPy.kern.RBF(self.input_dim)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim), + #GPy.kern.Linear(self.input_dim) + GPy.kern.Bias(self.input_dim) + GPy.kern.White(self.input_dim), + #GPy.kern.Linear(self.input_dim, ARD=True) + GPy.kern.Bias(self.input_dim) + GPy.kern.White(self.input_dim), + GPy.kern.Linear([1,3,6,7], ARD=True) + GPy.kern.RBF([0,5,8], ARD=True) + GPy.kern.White(self.input_dim), ) self.q_x_mean = np.random.randn(self.input_dim)[None] self.q_x_variance = np.exp(.5*np.random.randn(self.input_dim))[None] From 54239555a1a099d423f28458ba8ebc5bb25a33ad Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Mar 2014 12:03:37 +0000 Subject: [PATCH 14/19] psi_stat slices for kernels --- GPy/kern/_src/add.py | 55 +++++++++++++++++++++++++---------------- GPy/kern/_src/kern.py | 18 +++++--------- GPy/kern/_src/linear.py | 3 +-- GPy/kern/_src/rbf.py | 6 ++++- GPy/kern/_src/static.py | 28 +++++++++++++++++++++ 5 files changed, 74 insertions(+), 36 deletions(-) diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py index 8a3cefaf..fdebdfac 100644 --- a/GPy/kern/_src/add.py +++ b/GPy/kern/_src/add.py @@ -17,6 +17,11 @@ class Add(CombinationKernel): @Cache_this(limit=2, force_kwargs=['which_parts']) def K(self, X, X2=None, which_parts=None): + """ + Add all kernels together. + If a list of parts (of this kernel!) `which_parts` is given, only + the parts of the list are taken to compute the covariance. + """ assert X.shape[1] == self.input_dim if which_parts is None: which_parts = self.parts @@ -25,6 +30,22 @@ class Add(CombinationKernel): which_parts = [which_parts] return reduce(np.add, (p.K(X, X2) for p in which_parts)) + @Cache_this(limit=2, force_kwargs=['which_parts']) + def Kdiag(self, X, which_parts=None): + assert X.shape[1] == self.input_dim + if which_parts is None: + which_parts = self.parts + elif not isinstance(which_parts, (list, tuple)): + # if only one part is given + which_parts = [which_parts] + return reduce(np.add, (p.Kdiag(X) for p in which_parts)) + + def update_gradients_full(self, dL_dK, X, X2=None): + [p.update_gradients_full(dL_dK, X, X2) for p in self.parts] + + def update_gradients_diag(self, dL_dK, X): + [p.update_gradients_diag(dL_dK, X) for p in self.parts] + def gradients_X(self, dL_dK, X, X2=None): """Compute the gradient of the objective function with respect to X. @@ -36,18 +57,9 @@ class Add(CombinationKernel): :type X2: np.ndarray (num_inducing x input_dim)""" target = np.zeros(X.shape) - for p in self.parts: - target[:, p.active_dims] += p.gradients_X(dL_dK, X, X2) + [target.__setitem__([Ellipsis, p.active_dims], target[:, p.active_dims]+p.gradients_X(dL_dK, X, X2)) for p in self.parts] return target - @Cache_this(limit=2, force_kwargs=['which_parts']) - def Kdiag(self, X, which_parts=None): - assert X.shape[1] == self.input_dim - if which_parts is None: - which_parts = self.parts - return sum([p.Kdiag(X) for p in which_parts]) - - def psi0(self, Z, variational_posterior): return reduce(np.add, (p.psi0(Z, variational_posterior) for p in self.parts)) @@ -56,7 +68,7 @@ class Add(CombinationKernel): def psi2(self, Z, variational_posterior): psi2 = reduce(np.add, (p.psi2(Z, variational_posterior) for p in self.parts)) - return psi2 + #return psi2 # compute the "cross" terms from static import White, Bias from rbf import RBF @@ -65,23 +77,24 @@ class Add(CombinationKernel): #ffrom fixed import Fixed for p1, p2 in itertools.combinations(self.parts, 2): - i1, i2 = p1.active_dims, p2.active_dims + # i1, i2 = p1.active_dims, p2.active_dims # white doesn;t combine with anything if isinstance(p1, White) or isinstance(p2, White): pass # rbf X bias #elif isinstance(p1, (Bias, Fixed)) and isinstance(p2, (RBF, RBFInv)): elif isinstance(p1, Bias) and isinstance(p2, (RBF, Linear)): - # manual override for slicing: - p2._sliced_X = p1._sliced_X = True - tmp = p2.psi1(Z[:,i2], variational_posterior[:, i1]) + tmp = p2.psi1(Z, variational_posterior) psi2 += p1.variance * (tmp[:, :, None] + tmp[:, None, :]) #elif isinstance(p2, (Bias, Fixed)) and isinstance(p1, (RBF, RBFInv)): elif isinstance(p2, Bias) and isinstance(p1, (RBF, Linear)): - # manual override for slicing: - p2._sliced_X = p1._sliced_X = True - tmp = p1.psi1(Z[:,i1], variational_posterior[:, i2]) + tmp = p1.psi1(Z, variational_posterior) psi2 += p2.variance * (tmp[:, :, None] + tmp[:, None, :]) + elif isinstance(p2, (RBF, Linear)) and isinstance(p1, (RBF, Linear)): + assert np.intersect1d(p1.active_dims, p2.active_dims).size == 0, "only non overlapping kernel dimensions allowed so far" + tmp1 = p1.psi1(Z, variational_posterior) + tmp2 = p2.psi1(Z, variational_posterior) + psi2 += (tmp1[:, :, None] * tmp2[:, None, :]) + (tmp2[:, :, None] * tmp1[:, None, :]) else: raise NotImplementedError, "psi2 cannot be computed for this kernel" return psi2 @@ -98,7 +111,7 @@ class Add(CombinationKernel): continue elif isinstance(p2, Bias): eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2. - else: + else:# np.setdiff1d(p1.active_dims, ar2, assume_unique): # TODO: Careful, not correct for overlapping active_dims eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z, variational_posterior) * 2. p1.update_gradients_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior) @@ -114,9 +127,9 @@ class Add(CombinationKernel): if isinstance(p2, White): continue elif isinstance(p2, Bias): - eff_dL_dpsi1 += 0#dL_dpsi2.sum(1) * p2.variance * 2. + eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2. else: - eff_dL_dpsi1 += 0#dL_dpsi2.sum(1) * p2.psi1(Z, variational_posterior) * 2. + eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z, variational_posterior) * 2. target[:, p1.active_dims] += p1.gradients_Z_expectations(eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior) return target diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index 8feb9a04..8a24e24a 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -15,6 +15,7 @@ class Kern(Parameterized): # found in kernel_slice_operations __metaclass__ = KernCallsViaSlicerMeta #=========================================================================== + _debug=False def __init__(self, input_dim, name, *a, **kw): """ The base class for a kernel: a positive definite function @@ -27,12 +28,12 @@ class Kern(Parameterized): """ super(Kern, self).__init__(name=name, *a, **kw) if isinstance(input_dim, int): - self.active_dims = slice(0, input_dim) + self.active_dims = np.r_[0:input_dim] self.input_dim = input_dim else: - self.active_dims = input_dim + self.active_dims = np.r_[input_dim] self.input_dim = len(self.active_dims) - self._sliced_X = False + self._sliced_X = 0 @Cache_this(limit=10)#, ignore_args = (0,)) def _slice_X(self, X): @@ -60,14 +61,13 @@ class Kern(Parameterized): raise NotImplementedError def gradients_X_diag(self, dL_dKdiag, X): raise NotImplementedError - + def update_gradients_full(self, dL_dK, X, X2): """Set the gradients of all parameters when doing full (N) inference.""" raise NotImplementedError def update_gradients_diag(self, dL_dKdiag, X): """Set the gradients for all parameters for the derivative of the diagonal of the covariance w.r.t the kernel parameters.""" raise NotImplementedError - def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): """ Set the gradients of all parameters when doing inference with @@ -188,7 +188,7 @@ class Kern(Parameterized): class CombinationKernel(Kern): def __init__(self, kernels, name): assert all([isinstance(k, Kern) for k in kernels]) - input_dim = reduce(np.union1d, (np.r_[x.active_dims] for x in kernels)) + input_dim = reduce(np.union1d, (x.active_dims for x in kernels)) super(CombinationKernel, self).__init__(input_dim, name) self.add_parameters(*kernels) @@ -196,12 +196,6 @@ class CombinationKernel(Kern): def parts(self): return self._parameters_ - def update_gradients_full(self, dL_dK, X, X2=None): - [p.update_gradients_full(dL_dK, X, X2) for p in self.parts] - - def update_gradients_diag(self, dL_dK, X): - [p.update_gradients_diag(dL_dK, X) for p in self.parts] - def input_sensitivity(self): in_sen = np.zeros((self.num_params, self.input_dim)) for i, p in enumerate(self.parts): diff --git a/GPy/kern/_src/linear.py b/GPy/kern/_src/linear.py index 60645d11..f2ac0124 100644 --- a/GPy/kern/_src/linear.py +++ b/GPy/kern/_src/linear.py @@ -147,7 +147,6 @@ class Linear(Kern): mu = variational_posterior.mean S = variational_posterior.variance mu2S = np.square(mu)+S - _dpsi2_dvariance, _, _, _, _ = linear_psi_comp._psi2computations(self.variances, Z, mu, S, gamma) grad = np.einsum('n,nq,nq->q',dL_dpsi0,gamma,mu2S) + np.einsum('nm,nq,mq,nq->q',dL_dpsi1,gamma,Z,mu) +\ np.einsum('nmo,nmoq->q',dL_dpsi2,_dpsi2_dvariance) @@ -175,7 +174,7 @@ class Linear(Kern): mu = variational_posterior.mean S = variational_posterior.variance _, _, _, _, _dpsi2_dZ = linear_psi_comp._psi2computations(self.variances, Z, mu, S, gamma) - + grad = np.einsum('nm,nq,q,nq->mq',dL_dpsi1,gamma, self.variances,mu) +\ np.einsum('nmo,noq->mq',dL_dpsi2,_dpsi2_dZ) diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 7ba1f35d..341d46a7 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -19,7 +19,6 @@ class RBF(Stationary): k(r) = \sigma^2 \exp \\bigg(- \\frac{1}{2} r^2 \\bigg) """ - def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='rbf'): super(RBF, self).__init__(input_dim, variance, lengthscale, ARD, name) self.weave_options = {} @@ -81,6 +80,8 @@ class RBF(Stationary): #contributions from psi0: self.variance.gradient = np.sum(dL_dpsi0) + if self._debug: + num_grad = self.lengthscale.gradient.copy() self.lengthscale.gradient = 0. #from psi1 @@ -100,6 +101,8 @@ class RBF(Stationary): else: self.lengthscale.gradient += self._weave_psi2_lengthscale_grads(dL_dpsi2, psi2, Zdist_sq, S, mudist_sq, l2) + if self._debug: + import ipdb;ipdb.set_trace() self.variance.gradient += 2.*np.sum(dL_dpsi2 * psi2)/self.variance else: @@ -150,6 +153,7 @@ class RBF(Stationary): grad_mu = (dL_dpsi1[:, :, None] * _dpsi1_dmu).sum(axis=1) grad_S = (dL_dpsi1[:, :, None] * _dpsi1_dS).sum(axis=1) grad_gamma = (dL_dpsi1[:,:,None] * _dpsi1_dgamma).sum(axis=1) + #psi2 grad_mu += (dL_dpsi2[:, :, :, None] * _dpsi2_dmu).reshape(ndata,-1,self.input_dim).sum(axis=1) grad_S += (dL_dpsi2[:, :, :, None] * _dpsi2_dS).reshape(ndata,-1,self.input_dim).sum(axis=1) diff --git a/GPy/kern/_src/static.py b/GPy/kern/_src/static.py index f344357c..387c92c6 100644 --- a/GPy/kern/_src/static.py +++ b/GPy/kern/_src/static.py @@ -89,3 +89,31 @@ class Bias(Static): def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): self.variance.gradient = dL_dpsi0.sum() + dL_dpsi1.sum() + 2.*self.variance*dL_dpsi2.sum() +class Fixed(Static): + def __init__(self, input_dim, covariance_matrix, variance=1., name='fixed'): + """ + :param input_dim: the number of input dimensions + :type input_dim: int + :param variance: the variance of the kernel + :type variance: float + """ + super(Bias, self).__init__(input_dim, variance, name) + self.fixed_K = covariance_matrix + def K(self, X, X2): + return self.variance * self.fixed_K + + def Kdiag(self, X): + return self.variance * self.fixed_K.diag() + + def update_gradients_full(self, dL_dK, X, X2=None): + self.variance.gradient = np.einsum('ij,ij', dL_dK, self.fixed_K) + + def update_gradients_diag(self, dL_dKdiag, X): + self.variance.gradient = np.einsum('i,i', dL_dKdiag, self.fixed_K) + + def psi2(self, Z, variational_posterior): + return np.zeros((variational_posterior.shape[0], Z.shape[0], Z.shape[0]), dtype=np.float64) + + def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + self.variance.gradient = dL_dpsi0.sum() + From 5027e8e312caa99946f4817a446c12779b5b0934 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Mar 2014 12:03:47 +0000 Subject: [PATCH 15/19] diagonal add kmm --- GPy/inference/latent_function_inference/var_dtc.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index 52f44cdf..6239e5a4 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -3,6 +3,7 @@ from posterior import Posterior from ...util.linalg import jitchol, backsub_both_sides, tdot, dtrtrs, dtrtri, dpotri, dpotrs, symmetrify +from ...util import diag from ...core.parameterization.variational import VariationalPosterior import numpy as np from ...util.misc import param_to_array @@ -28,7 +29,7 @@ class VarDTC(object): def set_limit(self, limit): self.get_trYYT.limit = limit self.get_YYTfactor.limit = limit - + def _get_trYYT(self, Y): return param_to_array(np.sum(np.square(Y))) @@ -77,10 +78,10 @@ class VarDTC(object): num_inducing = Z.shape[0] num_data = Y.shape[0] # kernel computations, using BGPLVM notation - - Kmm = kern.K(Z) +np.eye(Z.shape[0]) * self.const_jitter - Lm = jitchol(Kmm+np.eye(Z.shape[0])*self.const_jitter) + Kmm = kern.K(Z).copy() + diag.add(Kmm, self.const_jitter) + Lm = jitchol(Kmm) # The rather complex computations of A if uncertain_inputs: @@ -169,7 +170,6 @@ class VarDTC(object): Bi, _ = dpotri(LB, lower=1) symmetrify(Bi) Bi = -dpotri(LB, lower=1)[0] - from ...util import diag diag.add(Bi, 1) woodbury_inv = backsub_both_sides(Lm, Bi) @@ -238,7 +238,8 @@ class VarDTCMissingData(object): dL_dKmm = 0 log_marginal = 0 - Kmm = kern.K(Z) + Kmm = kern.K(Z).copy() + diag.add(Kmm, self.const_jitter) #factor Kmm Lm = jitchol(Kmm) if uncertain_inputs: LmInv = dtrtri(Lm) @@ -324,7 +325,6 @@ class VarDTCMissingData(object): Bi, _ = dpotri(LB, lower=1) symmetrify(Bi) Bi = -dpotri(LB, lower=1)[0] - from ...util import diag diag.add(Bi, 1) woodbury_inv_all[:, :, ind] = backsub_both_sides(Lm, Bi)[:,:,None] From 2200c5c30b8c98e653b7a0433ac02dce20835075 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Mar 2014 12:05:54 +0000 Subject: [PATCH 16/19] uncertain_inputs_example plot changed --- GPy/examples/regression.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index cc23410a..7cd1e964 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -468,7 +468,7 @@ def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, opt def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True): """Run a 1D example of a sparse GP regression with uncertain inputs.""" - fig, axes = pb.subplots(1, 2, figsize=(12, 5)) + fig, axes = pb.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True) # sample inputs and outputs S = np.ones((20, 1)) From 53e071b892a5b2bf0ab8efee05b358efed9f4ec8 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Mar 2014 12:06:21 +0000 Subject: [PATCH 17/19] gradient check and debug options --- GPy/core/gp.py | 6 +++--- GPy/core/model.py | 7 ++++++- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 6441561b..2cff4341 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -70,7 +70,7 @@ class GP(Model): self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y, Y_metadata=self.Y_metadata) self.likelihood.update_gradients(np.diag(grad_dict['dL_dK'])) self.kern.update_gradients_full(grad_dict['dL_dK'], self.X) - + def log_likelihood(self): return self._log_marginal_likelihood @@ -186,7 +186,7 @@ class GP(Model): """ assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ..plotting.matplot_dep import models_plots - models_plots.plot_fit_f(self,*args,**kwargs) + return models_plots.plot_fit_f(self,*args,**kwargs) def plot(self, *args, **kwargs): """ @@ -207,7 +207,7 @@ class GP(Model): """ assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ..plotting.matplot_dep import models_plots - models_plots.plot_fit(self,*args,**kwargs) + return models_plots.plot_fit(self,*args,**kwargs) def _getstate(self): """ diff --git a/GPy/core/model.py b/GPy/core/model.py index 6a6fe1ba..710c1b22 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -339,7 +339,7 @@ class Model(Parameterized): print "No free parameters to check" return - gradient = self.objective_function_gradients(x) + gradient = self.objective_function_gradients(x).copy() np.where(gradient == 0, 1e-312, gradient) ret = True for nind, xind in itertools.izip(param_index, transformed_index): @@ -350,7 +350,12 @@ class Model(Parameterized): f2 = self.objective_function(xx) numerical_gradient = (f1 - f2) / (2 * step) if _debug: + for p in self.kern.flattened_parameters: + p._parent_._debug=True self.gradient[xind] = numerical_gradient + self._set_params_transformed(x) + for p in self.kern.flattened_parameters: + p._parent_._debug=False if np.all(gradient[xind]==0): ratio = (f1-f2) == gradient[xind] else: ratio = (f1 - f2) / (2 * step * gradient[xind]) difference = np.abs((f1 - f2) / 2 / step - gradient[xind]) From abc7545e0993dd927b20ee4a5854fd72e8cb0694 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Mar 2014 13:00:02 +0000 Subject: [PATCH 18/19] kernel slicer --- GPy/kern/_src/kernel_slice_operations.py | 108 +++++++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 GPy/kern/_src/kernel_slice_operations.py diff --git a/GPy/kern/_src/kernel_slice_operations.py b/GPy/kern/_src/kernel_slice_operations.py new file mode 100644 index 00000000..c1774a35 --- /dev/null +++ b/GPy/kern/_src/kernel_slice_operations.py @@ -0,0 +1,108 @@ +''' +Created on 11 Mar 2014 + +@author: maxz +''' +from ...core.parameterization.parameterized import ParametersChangedMeta + +class KernCallsViaSlicerMeta(ParametersChangedMeta): + def __call__(self, *args, **kw): + instance = super(ParametersChangedMeta, self).__call__(*args, **kw) + instance.K = _slice_wrapper(instance, instance.K) + instance.Kdiag = _slice_wrapper(instance, instance.Kdiag, True) + instance.update_gradients_full = _slice_wrapper(instance, instance.update_gradients_full, False, True) + instance.update_gradients_diag = _slice_wrapper(instance, instance.update_gradients_diag, True, True) + instance.gradients_X = _slice_wrapper(instance, instance.gradients_X, False, True) + instance.gradients_X_diag = _slice_wrapper(instance, instance.gradients_X_diag, True, True) + instance.psi0 = _slice_wrapper(instance, instance.psi0, False, False) + instance.psi1 = _slice_wrapper(instance, instance.psi1, False, False) + instance.psi2 = _slice_wrapper(instance, instance.psi2, False, False) + instance.update_gradients_expectations = _slice_wrapper(instance, instance.update_gradients_expectations, psi_stat=True) + instance.gradients_Z_expectations = _slice_wrapper(instance, instance.gradients_Z_expectations, psi_stat_Z=True) + instance.gradients_qX_expectations = _slice_wrapper(instance, instance.gradients_qX_expectations, psi_stat=True) + instance.parameters_changed() + return instance + +def _slice_wrapper(kern, operation, diag=False, derivative=False, psi_stat=False, psi_stat_Z=False): + """ + This method wraps the functions in kernel to make sure all kernels allways see their respective input dimension. + The different switches are: + diag: if X2 exists + derivative: if first arg is dL_dK + psi_stat: if first 3 args are dL_dpsi0..2 + psi_stat_Z: if first 2 args are dL_dpsi1..2 + """ + if derivative: + if diag: + def x_slice_wrapper(dL_dK, X): + X = kern._slice_X(X) if not kern._sliced_X else X + kern._sliced_X += 1 + try: + ret = operation(dL_dK, X) + except: + raise + finally: + kern._sliced_X -= 1 + return ret + else: + def x_slice_wrapper(dL_dK, X, X2=None): + X, X2 = kern._slice_X(X) if not kern._sliced_X else X, kern._slice_X(X2) if X2 is not None and not kern._sliced_X else X2 + kern._sliced_X += 1 + try: + ret = operation(dL_dK, X, X2) + except: + raise + finally: + kern._sliced_X -= 1 + return ret + elif psi_stat: + def x_slice_wrapper(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + Z, variational_posterior = kern._slice_X(Z) if not kern._sliced_X else Z, kern._slice_X(variational_posterior) if not kern._sliced_X else variational_posterior + kern._sliced_X += 1 + try: + ret = operation(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior) + except: + raise + finally: + kern._sliced_X -= 1 + return ret + elif psi_stat_Z: + def x_slice_wrapper(dL_dpsi1, dL_dpsi2, Z, variational_posterior): + Z, variational_posterior = kern._slice_X(Z) if not kern._sliced_X else Z, kern._slice_X(variational_posterior) if not kern._sliced_X else variational_posterior + kern._sliced_X += 1 + try: + ret = operation(dL_dpsi1, dL_dpsi2, Z, variational_posterior) + except: + raise + finally: + kern._sliced_X -= 1 + return ret + else: + if diag: + def x_slice_wrapper(X, *args, **kw): + X = kern._slice_X(X) if not kern._sliced_X else X + kern._sliced_X += 1 + try: + ret = operation(X, *args, **kw) + except: + raise + finally: + kern._sliced_X -= 1 + return ret + else: + def x_slice_wrapper(X, X2=None, *args, **kw): + X, X2 = kern._slice_X(X) if not kern._sliced_X else X, kern._slice_X(X2) if X2 is not None and not kern._sliced_X else X2 + kern._sliced_X += 1 + try: + ret = operation(X, X2, *args, **kw) + except: raise + finally: + kern._sliced_X -= 1 + return ret + x_slice_wrapper._operation = operation + x_slice_wrapper.__name__ = ("slicer("+operation.__name__ + +(","+str(bool(diag)) if diag else'') + +(','+str(bool(derivative)) if derivative else '') + +')') + x_slice_wrapper.__doc__ = "**sliced**\n" + (operation.__doc__ or "") + return x_slice_wrapper \ No newline at end of file From b975a45cd23dfdbe4689e8fdc6983816080462fe Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Mar 2014 13:02:52 +0000 Subject: [PATCH 19/19] copy --- GPy/core/parameterization/parameter_core.py | 28 ++++++++++----------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 5727bc17..d1122f79 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-11' +__updated__ = '2014-03-12' class HierarchyError(Exception): """ @@ -796,27 +796,27 @@ class Parameterizable(OptimizationHandlable): """ if not param in self._parameters_: raise RuntimeError, "Parameter {} does not belong to this object, remove parameters directly from their respective parents".format(param._short()) - + start = sum([p.size for p in self._parameters_[:param._parent_index_]]) self._remove_parameter_name(param) self.size -= param.size del self._parameters_[param._parent_index_] - + param._disconnect_parent() 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 = parent._parent_ - + 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 @@ -829,29 +829,26 @@ class Parameterizable(OptimizationHandlable): old_size = 0 self._param_array_ = np.empty(self.size, dtype=np.float64) self._gradient_array_ = np.empty(self.size, dtype=np.float64) - + self._param_slices_ = [] - for i, p in enumerate(self._parameters_): p._parent_ = self p._parent_index_ = i - + pslice = slice(old_size, old_size+p.size) - # first connect all children p._propagate_param_grad(self._param_array_[pslice], self._gradient_array_[pslice]) - # then connect children to self self._param_array_[pslice] = p._param_array_.ravel()#, requirements=['C', 'W']).ravel(order='C') self._gradient_array_[pslice] = p._gradient_array_.ravel()#, requirements=['C', 'W']).ravel(order='C') - + if not p._param_array_.flags['C_CONTIGUOUS']: import ipdb;ipdb.set_trace() p._param_array_.data = self._param_array_[pslice].data p._gradient_array_.data = self._gradient_array_[pslice].data - + self._param_slices_.append(pslice) - + self._add_parameter_name(p, ignore_added_names=ignore_added_names) old_size += p.size @@ -862,12 +859,13 @@ class Parameterizable(OptimizationHandlable): self.parameters_changed() def _pass_through_notify_observers(self, which): self.notify_observers(which) - + #=========================================================================== # TODO: not working yet #=========================================================================== def copy(self): """Returns a (deep) copy of the current model""" + raise NotImplementedError, "Copy is not yet implemented, TODO: Observable hierarchy" import copy from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView from .lists_and_dicts import ArrayList