diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py index 2a0ba7c1..750788e9 100644 --- a/GPy/core/gp_base.py +++ b/GPy/core/gp_base.py @@ -116,12 +116,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': @@ -131,15 +133,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 8ab00e10..d1e2885a 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -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)" - new_parts = [parts.independent_outputs.IndependentOutputs(p) for p in k.parts] - return kern(k.input_dim+1,new_parts) + _parts = [parts.independent_outputs.IndependentOutputs(p) for p in k.parts] + return kern(k.input_dim+1,_parts) 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]