merged

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)

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)

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)

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]