From db7485b906cb2d1d0fb740002b206c68825c64c5 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 8 Jul 2013 13:06:02 +0100 Subject: [PATCH 1/2] fixed a bug in constructor of periodic_matern52 --- GPy/core/gp_base.py | 24 ++++++++++++++++++------ GPy/examples/regression.py | 6 +++--- GPy/kern/constructors.py | 6 +++--- 3 files changed, 24 insertions(+), 12 deletions(-) diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py index b82f3298..63568968 100644 --- a/GPy/core/gp_base.py +++ b/GPy/core/gp_base.py @@ -91,12 +91,14 @@ class GPBase(Model): else: raise NotImplementedError, "Cannot define a frame with more than two input dimensions" - def plot(self, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, samples=0, fignum=None, ax=None): + def plot(self, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, samples=0, fignum=None, ax=None, fixed_inputs=[], linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue']): """ TODO: Docstrings! :param levels: for 2D plotting, the number of contour levels to use is ax is None, create a new figure + + fixed_inputs: a list of tuple [(i,v), (i,v)...], specifying that input index i should be set to value v. """ # TODO include samples if which_data == 'all': @@ -106,15 +108,25 @@ class GPBase(Model): fig = pb.figure(num=fignum) ax = fig.add_subplot(111) - if self.X.shape[1] == 1: + plotdims = self.input_dim - len(fixed_inputs) + + if plotdims == 1: Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now - Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits) - m, _, lower, upper = self.predict(Xnew, which_parts=which_parts) + fixed_dims = np.array([i for i,v in fixed_inputs]) + freedim = np.setdiff1d(np.arange(self.input_dim),fixed_dims) + + Xnew, xmin, xmax = x_frame1D(Xu[:,freedim], plot_limits=plot_limits) + Xgrid = np.empty((Xnew.shape[0],self.input_dim)) + Xgrid[:,freedim] = Xnew + for i,v in fixed_inputs: + Xgrid[:,i] = v + + m, _, lower, upper = self.predict(Xgrid, which_parts=which_parts) for d in range(m.shape[1]): - gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax) - ax.plot(Xu[which_data], self.likelihood.data[which_data, d], 'kx', mew=1.5) + gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax, edgecol=linecol, fillcol=fillcol) + ax.plot(Xu[which_data,freedim], self.likelihood.data[which_data, d], 'kx', mew=1.5) ymin, ymax = min(np.append(self.likelihood.data, lower)), max(np.append(self.likelihood.data, upper)) ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin) ax.set_xlim(xmin, xmax) diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index 21b435e7..452167ce 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -83,7 +83,7 @@ def coregionalisation_toy2(optim_iters=100): Y = np.vstack((Y1,Y2)) k1 = GPy.kern.rbf(1) + GPy.kern.bias(1) - k2 = GPy.kern.Coregionalise(2,1) + k2 = GPy.kern.coregionalise(2,1) k = k1.prod(k2,tensor=True) m = GPy.models.GPRegression(X,Y,kernel=k) m.constrain_fixed('.*rbf_var',1.) @@ -114,7 +114,7 @@ def coregionalisation_toy(optim_iters=100): Y = np.vstack((Y1,Y2)) k1 = GPy.kern.rbf(1) - k2 = GPy.kern.Coregionalise(2,2) + k2 = GPy.kern.coregionalise(2,2) k = k1.prod(k2,tensor=True) m = GPy.models.GPRegression(X,Y,kernel=k) m.constrain_fixed('.*rbf_var',1.) @@ -149,7 +149,7 @@ def coregionalisation_sparse(optim_iters=100): Z = np.hstack((np.random.rand(num_inducing,1)*8,np.random.randint(0,2,num_inducing)[:,None])) k1 = GPy.kern.rbf(1) - k2 = GPy.kern.Coregionalise(2,2) + k2 = GPy.kern.coregionalise(2,2) k = k1.prod(k2,tensor=True) + GPy.kern.white(2,0.001) m = GPy.models.SparseGPRegression(X,Y,kernel=k,Z=Z) diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index 697f3554..d1e2885a 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -227,7 +227,7 @@ def periodic_Matern52(input_dim, variance=1., lengthscale=None, period=2 * np.pi :param n_freq: the number of frequencies considered for the periodic subspace :type n_freq: int """ - part = parts.periodic_Matern52part(input_dim, variance, lengthscale, period, n_freq, lower, upper) + part = parts.periodic_Matern52.PeriodicMatern52(input_dim, variance, lengthscale, period, n_freq, lower, upper) return kern(input_dim, [part]) def prod(k1,k2,tensor=False): @@ -296,5 +296,5 @@ def independent_outputs(k): """ for sl in k.input_slices: assert (sl.start is None) and (sl.stop is None), "cannot adjust input slices! (TODO)" - parts = [independent_outputs.IndependentOutputs(p) for p in k.parts] - return kern(k.input_dim+1,parts) + _parts = [parts.independent_outputs.IndependentOutputs(p) for p in k.parts] + return kern(k.input_dim+1,_parts) From 58a59dcab081c4f07a3c6992f126cc7d0f5d867e Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 8 Jul 2013 14:20:42 +0100 Subject: [PATCH 2/2] created a hierarchical kernel --- GPy/kern/parts/hierarchical.py | 75 ++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 GPy/kern/parts/hierarchical.py diff --git a/GPy/kern/parts/hierarchical.py b/GPy/kern/parts/hierarchical.py new file mode 100644 index 00000000..98f4a183 --- /dev/null +++ b/GPy/kern/parts/hierarchical.py @@ -0,0 +1,75 @@ +# Copyright (c) 2012, James Hesnsman +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +from kernpart import Kernpart +import numpy as np +from independent_outputs import index_to_slices + +class hierarchical(Kernpart): + """ + A kernel part which can reopresent a hierarchy of indepencnce: a gerenalisation of independent_outputs + + """ + def __init__(self,parts): + self.levels = len(parts) + self.input_dim = parts[0].input_dim + 1 + self.num_params = np.sum([k.num_params for k in parts]) + self.name = 'hierarchy' + self.parts = parts + + self.param_starts = np.hstack((0,np.cumsum([k.num_params for k in self.parts[:-1]]))) + self.param_stops = np.cumsum([k.num_params for k in self.parts]) + + def _get_params(self): + return np.hstack([k._get_params() for k in self.parts]) + + def _set_params(self,x): + [k._set_params(x[start:stop]) for start, stop in zip(self.param_starts, self.param_stops)] + + def _get_param_names(self): + return self.k._get_param_names() + + def _sort_slices(self,X,X2): + slices = [index_to_slices(x) for x in X[-self.levels:].T] + X = X[:-self.levels] + if X2 is None: + slices2 = slices + X2 = X + else: + slices2 = [index_to_slices(x) for x in X2[-self.levels:].T] + X2 = X2[:-self.levels] + return X, X2, slices, slices2 + + def K(self,X,X2,target): + X, X2, slices, slices2 = self._sort_slices(X,X2) + + [[[k.K(X[s],X2[s2],target[s,s2]) for s in slices_i] for s2 in slices_j] for k,slices_i,slices_j in zip(self.parts,slices,slices2)] + + def Kdiag(self,X,target): + raise NotImplementedError + #X,slices = X[:,:-1],index_to_slices(X[:,-1]) + #[[self.k.Kdiag(X[s],target[s]) for s in slices_i] for slices_i in slices] + + def dK_dtheta(self,dL_dK,X,X2,target): + [[[k.dK_dtheta(dL_dK[s,s2],X[s],X2[s2],target[p_start:p_stop]) for s in slices_i] for s2 in slices_j] for k,slices_i,slices_j, p_start, p_stop in zip(self.parts, slices, slices2, self.param_starts, self.param_stops)] + + + def dK_dX(self,dL_dK,X,X2,target): + raise NotImplementedError + #X,slices = X[:,:-1],index_to_slices(X[:,-1]) + #if X2 is None: + #X2,slices2 = X,slices + #else: + #X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1]) + #[[[self.k.dK_dX(dL_dK[s,s2],X[s],X2[s2],target[s,:-1]) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)] +# + def dKdiag_dX(self,dL_dKdiag,X,target): + raise NotImplementedError + #X,slices = X[:,:-1],index_to_slices(X[:,-1]) + #[[self.k.dKdiag_dX(dL_dKdiag[s],X[s],target[s,:-1]) for s in slices_i] for slices_i in slices] + + + def dKdiag_dtheta(self,dL_dKdiag,X,target): + raise NotImplementedError + #X,slices = X[:,:-1],index_to_slices(X[:,-1]) + #[[self.k.dKdiag_dX(dL_dKdiag[s],X[s],target) for s in slices_i] for slices_i in slices]