various work on BGPLVM oil demo

GPy/core/model.py

@@ -423,6 +423,28 @@ class model(parameterised):
                 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
 
+    def input_sensitivity(self):
+        """
+        return an array describing the sesitivity of the model to each input
+
+        NB. Right now, we're basing this on the lengthscales (or variances) of the kernel.
+        TODO: proper sensitivity analysis
+        """
+
+        if not hasattr(self,'kern'):
+            raise ValueError, "this model has no kernel"
+
+        k = [p for p in self.kern.parts if p.name in ['rbf','linear']]
+        if (not len(k)==1) or (not k[0].ARD):
+            raise ValueError, "cannot determine sensitivity for this kernel"
+        k = k[0]
+
+        if k.name=='rbf':
+            return k.lengthscale
+        elif k.name=='linear':
+            return 1./k.variances
+
+
     def pseudo_EM(self,epsilon=.1,**kwargs):
         """
         TODO: Should this not bein the GP class?
@@ -469,3 +491,4 @@ class model(parameterised):
             iteration += 1
             if stop:
                 print "%s iterations." %iteration
+

GPy/examples/dimensionality_reduction.py

@@ -47,7 +47,7 @@ def GPLVM_oil_100(optimize=True,M=15):
 
     # create simple GP model
     kernel = GPy.kern.rbf(6, ARD = True) + GPy.kern.bias(6)
-    m = GPy.models.Bayesian_GPLVM(data['X'], 6, kernel=kernel, M=M)
+    m = GPy.models.GPLVM(data['X'], 6, kernel=kernel, M=M)
     m.data_labels = data['Y'].argmax(axis=1)
 
     # optimize
@@ -60,22 +60,32 @@ def GPLVM_oil_100(optimize=True,M=15):
     m.plot_latent(labels=m.data_labels)
     return m
 
-def BGPLVM_oil_100(optimize=True):
+def BGPLVM_oil(optimize=True,N=100,Q=10,M=15):
-    data = GPy.util.datasets.oil_100()
+    data = GPy.util.datasets.oil()
 
     # create simple GP model
-    kernel = GPy.kern.rbf(6, ARD = True) + GPy.kern.bias(6)
+    kernel = GPy.kern.rbf(Q, ARD = True) + GPy.kern.bias(Q) + GPy.kern.white(Q,0.001)
-    m = GPy.models.GPLVM(data['X'], 6, kernel = kernel)
+    m = GPy.models.Bayesian_GPLVM(data['X'][:N], Q, kernel = kernel,M=M)
-    m.data_labels = data['Y'].argmax(axis=1)
+    m.data_labels = data['Y'][:N].argmax(axis=1)
+
+    #initial conditions
 
     # optimize
-    m.ensure_default_constraints()
     if optimize:
-        m.optimize(messages=1)
+        m.constrain_fixed('noise',0.05)
+        m.ensure_default_constraints()
+        m.optimize('scg',messages=1)
+        m.unconstrain('noise')
+        m.constrain_positive('noise')
+        m.optimize('scg',messages=1)
+    else:
+        m.ensure_default_constraints()
 
     # plot
     print(m)
-    m.plot_latent(labels=data['Y'].argmax(axis=1))
+    m.plot_latent(labels=m.data_labels)
+    pb.figure()
+    pb.bar(np.arange(m.kern.D),1./m.input_sensitivity())
     return m
 
 def oil_100():

GPy/models/Bayesian_GPLVM.py

@@ -22,12 +22,12 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
     :type init: 'PCA'|'random'
 
     """
-    def __init__(self, Y, Q, X = None, S = None, init='PCA', M=10, Z=None, kernel=None, **kwargs):
+    def __init__(self, Y, Q, X = None, X_uncertainty = None, init='PCA', M=10, Z=None, kernel=None, **kwargs):
         if X == None:
             X = self.initialise_latent(init, Q, Y)
 
-        if S is None:
+        if X_uncertainty is None:
-            S = np.ones_like(X) * 1e-2#
+            X_uncertainty = np.ones_like(X) * 0.5
 
         if Z is None:
             Z = np.random.permutation(X.copy())[:M]
@@ -37,7 +37,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
             kernel = kern.rbf(Q) + kern.white(Q)
 
 
-        sparse_GP.__init__(self, X, Gaussian(Y), kernel, Z=Z, X_uncertainty=S, **kwargs)
+        sparse_GP.__init__(self, X, Gaussian(Y), kernel, Z=Z, X_uncertainty=X_uncertainty, **kwargs)
 
     def _get_param_names(self):
         X_names = sum([['X_%i_%i'%(n,q) for q in range(self.Q)] for n in range(self.N)],[])
@@ -84,6 +84,14 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
     def _log_likelihood_gradients(self):
         return np.hstack((self.dL_dmuS().flatten(), sparse_GP._log_likelihood_gradients(self)))
 
-    def plot_latent(self, *args, **kwargs):
+    def plot_latent(self, which_indices=None,*args, **kwargs):
-        input_1, input_2 = GPLVM.plot_latent(self, *args, **kwargs)
+
+        if which_indices is None:
+            try:
+                input_1, input_2 = np.argsort(self.input_sensitivity())[:2]
+            except:
+                raise ValueError, "cannot Atomatically determine which dimensions to plot, please pass 'which_indices'"
+        else:
+            input_1, input_2 = which_indices
+        GPLVM.plot_latent(self, which_indices=[input_1, input_2],*args, **kwargs)
         pb.plot(self.Z[:, input_1], self.Z[:, input_2], '^w')

GPy/models/GPLVM.py

@@ -67,7 +67,7 @@ class GPLVM(GP):
         """
 
         util.plot.Tango.reset()
-        
+
         if labels is None:
             labels = np.ones(self.N)
         if which_indices is None:
@@ -77,23 +77,19 @@ class GPLVM(GP):
             if self.Q==2:
                 input_1, input_2 = 0,1
             else:
-                #try to find a linear of RBF kern in the kernel
+                try:
-                k = [p for p in self.kern.parts if p.name in ['rbf','linear']]
+                    input_1, input_2 = np.argsort(self.input_sensitivity())[:2]
-                if (not len(k)==1) or (not k[0].ARD):
+                except:
                     raise ValueError, "cannot Atomatically determine which dimensions to plot, please pass 'which_indices'"
-                input_1, input_2 = self.lengthscale_order()
+        else:
-                k = k[0]
+            input_1, input_2 = which_indices
-                if k.name=='rbf':
-                    input_1, input_2 = np.argsort(k.lengthscale)[:2]
-                elif k.name=='linear':
-                    input_1, input_2 = np.argsort(k.variances)[::-1][:2]
 
         #first, plot the output variance as a function of the latent space
         Xtest, xx,yy,xmin,xmax = util.plot.x_frame2D(self.X[:,[input_1, input_2]],resolution=resolution)
-	Xtest_full = np.zeros((Xtest.shape[0], self.X.shape[1]))
+        Xtest_full = np.zeros((Xtest.shape[0], self.X.shape[1]))
-	Xtest_full[:, :2] = Xtest
+        Xtest_full[:, :2] = Xtest
-	mu, var, low, up = self.predict(Xtest_full)
+        mu, var, low, up = self.predict(Xtest_full)
-	var = var.mean(axis=1) # this was var[:, :2] edit by Neil
+        var = var.mean(axis=1) # this was var[:, :2] edit by Neil
         pb.imshow(var.reshape(resolution,resolution).T[::-1,:],extent=[xmin[0],xmax[0],xmin[1],xmax[1]],cmap=pb.cm.binary,interpolation='bilinear')