Merge branch 'devel' of github.com:SheffieldML/GPy into devel

GPy/core/sparse_gp.py

@@ -80,8 +80,6 @@ class SparseGP(GPBase):
             self.psi2 = None
 
     def _computations(self):
-
-
         # factor Kmm
         self.Lm = jitchol(self.Kmm)
 
@@ -93,6 +91,8 @@ class SparseGP(GPBase):
                 psi2_beta = self.psi2.sum(0) * self.likelihood.precision
             evals, evecs = linalg.eigh(psi2_beta)
             clipped_evals = np.clip(evals, 0., 1e6) # TODO: make clipping configurable
+            if not np.array_equal(evals, clipped_evals):
+                pass#print evals
             tmp = evecs * np.sqrt(clipped_evals)
             tmp = tmp.T
         else:
@@ -114,7 +114,8 @@ class SparseGP(GPBase):
         # back substutue C into psi1Vf
         tmp, info1 = dtrtrs(self.Lm, np.asfortranarray(self.psi1Vf), lower=1, trans=0)
         self._LBi_Lmi_psi1Vf, _ = dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0)
-        tmp, info2 = dpotrs(self.LB, tmp, lower=1)
+        #tmp, info2 = dpotrs(self.LB, tmp, lower=1)
+        tmp, info2 = dtrtrs(self.LB, self._LBi_Lmi_psi1Vf, lower=1, trans=1)
         self.Cpsi1Vf, info3 = dtrtrs(self.Lm, tmp, lower=1, trans=1)
 
         # Compute dL_dKmm

GPy/core/svigp.py

@@ -92,12 +92,56 @@ class SVIGP(GPBase):
         self._vb_steplength_trace = []
 
     def getstate(self):
-        return GPBase.getstate(self)
-
+        steplength_params = [self.hbar_t, self.tau_t, self.gbar_t, self.gbar_t1, self.gbar_t2, self.hbar_tp, self.tau_tp, self.gbar_tp, self.adapt_param_steplength, self.adapt_vb_steplength, self.vb_steplength, self.param_steplength]
+        return GPBase.getstate(self) + \
+            [self.get_vb_param(),
+             self.Z,
+             self.num_inducing,
+             self.has_uncertain_inputs,
+             self.X_variance,
+             self.X_batch,
+             self.X_variance_batch,
+             steplength_params,
+             self.batchcounter,
+             self.batchsize,
+             self.epochs,
+             self.momentum,
+             self.data_prop,
+             self._param_trace,
+             self._param_steplength_trace,
+             self._vb_steplength_trace,
+             self._ll_trace,
+             self._grad_trace,
+             self.Y,
+             self._permutation,
+             self.iterations
+            ]
 
     def setstate(self, state):
-        return GPBase.setstate(self, state)
-
+        self.iterations = state.pop()
+        self._permutation = state.pop()
+        self.Y = state.pop()
+        self._grad_trace = state.pop()
+        self._ll_trace = state.pop()
+        self._vb_steplength_trace = state.pop()
+        self._param_steplength_trace = state.pop()
+        self._param_trace = state.pop()
+        self.data_prop = state.pop()
+        self.momentum = state.pop()
+        self.epochs = state.pop()
+        self.batchsize = state.pop()
+        self.batchcounter = state.pop()
+        steplength_params = state.pop()
+        (self.hbar_t, self.tau_t, self.gbar_t, self.gbar_t1, self.gbar_t2, self.hbar_tp, self.tau_tp, self.gbar_tp, self.adapt_param_steplength, self.adapt_vb_steplength, self.vb_steplength, self.param_steplength) = steplength_params
+        self.X_variance_batch = state.pop()
+        self.X_batch = state.pop()
+        self.X_variance = state.pop()
+        self.has_uncertain_inputs = state.pop()
+        self.num_inducing = state.pop()
+        self.Z = state.pop()
+        vb_param = state.pop()
+        GPBase.setstate(self, state)
+        self.set_vb_param(vb_param)        
 
     def _compute_kernel_matrices(self):
         # kernel computations, using BGPLVM notation

GPy/models/gplvm.py

@@ -63,6 +63,20 @@ class GPLVM(GP):
 
         return np.hstack((dL_dX.flatten(), GP._log_likelihood_gradients(self)))
 
+    def jacobian(self,X):
+        target = np.zeros((X.shape[0],X.shape[1],self.output_dim))
+        for i in range(self.output_dim):
+        	target[:,:,i]=self.kern.dK_dX(np.dot(self.Ki,self.likelihood.Y[:,i])[None, :],X,self.X)
+        return target
+   
+    def magnification(self,X):
+        target=np.zeros(X.shape[0])
+        J = np.zeros((X.shape[0],X.shape[1],self.output_dim))
+    	J=self.jacobian(X)
+        for i in range(X.shape[0]):
+		    target[i]=np.sqrt(pb.det(np.dot(J[i,:,:],np.transpose(J[i,:,:]))))
+        return target
+
     def plot(self):
         assert self.likelihood.Y.shape[1] == 2
         pb.scatter(self.likelihood.Y[:, 0], self.likelihood.Y[:, 1], 40, self.X[:, 0].copy(), linewidth=0, cmap=pb.cm.jet)
@@ -72,3 +86,6 @@ class GPLVM(GP):
 
     def plot_latent(self, *args, **kwargs):
         return util.plot_latent.plot_latent(self, *args, **kwargs)
+
+    def plot_magnification(self, *args, **kwargs):
+        return util.plot_latent.plot_magnification(self, *args, **kwargs)

GPy/util/datasets/connections.txt

@@ -0,0 +1,22 @@
+LFHD, RFHD
+RFHD, RBHD
+RBHD, LBHD
+LBHD, LFHD
+LELB, LWRB
+LWRB, LFIN
+LELB, LSHO
+LSHO, RSHO
+RSHO, STRN
+LSHO, STRN
+RSHO, RELB
+RELB, RWRB
+RWRB, RFIN
+LSHO, LFWT
+RSHO, RFWT
+LFWT, RFWT
+LFWT, LKNE
+RFWT, RKNE
+LKNE, LHEE
+RKNE, RHEE
+RMT5, RHEE
+LMT5, LHEE

GPy/util/plot_latent.py

@@ -98,3 +98,79 @@ def plot_latent(model, labels=None, which_indices=None,
         ax.figure.canvas.show()
         raw_input('Enter to continue')
     return ax
+
+def plot_magnification(model, labels=None, which_indices=None, 
+                resolution=60, ax=None, marker='o', s=40, 
+                fignum=None, plot_inducing=False, legend=True,
+                aspect='auto', updates=False):
+    """
+    :param labels: a np.array of size model.num_data containing labels for the points (can be number, strings, etc)
+    :param resolution: the resolution of the grid on which to evaluate the predictive variance
+    """
+    if ax is None:
+        fig = pb.figure(num=fignum)
+        ax = fig.add_subplot(111)
+    util.plot.Tango.reset()
+
+    if labels is None:
+        labels = np.ones(model.num_data)
+
+    input_1, input_2 = most_significant_input_dimensions(model, which_indices)
+
+    # first, plot the output variance as a function of the latent space
+    Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(model.X[:, [input_1, input_2]], resolution=resolution)
+    Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1]))
+    def plot_function(x):
+        Xtest_full[:, [input_1, input_2]] = x
+        mf=model.magnification(Xtest_full)
+        return mf
+    view = ImshowController(ax, plot_function, tuple(xmin) + tuple(xmax),
+                            resolution, aspect=aspect, interpolation='bilinear',
+                            cmap=pb.cm.gray)
+
+    # make sure labels are in order of input:
+    ulabels = []
+    for lab in labels:
+        if not lab in ulabels:
+            ulabels.append(lab)
+
+    marker = itertools.cycle(list(marker))
+
+    for i, ul in enumerate(ulabels):
+        if type(ul) is np.string_:
+            this_label = ul
+        elif type(ul) is np.int64:
+            this_label = 'class %i' % ul
+        else:
+            this_label = 'class %i' % i
+        m = marker.next()
+
+        index = np.nonzero(labels == ul)[0]
+        if model.input_dim == 1:
+            x = model.X[index, input_1]
+            y = np.zeros(index.size)
+        else:
+            x = model.X[index, input_1]
+            y = model.X[index, input_2]
+        ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label)
+
+    ax.set_xlabel('latent dimension %i' % input_1)
+    ax.set_ylabel('latent dimension %i' % input_2)
+
+    if not np.all(labels == 1.) and legend:
+        ax.legend(loc=0, numpoints=1)
+
+    ax.set_xlim(xmin[0], xmax[0])
+    ax.set_ylim(xmin[1], xmax[1])
+    ax.grid(b=False) # remove the grid if present, it doesn't look good
+    ax.set_aspect('auto') # set a nice aspect ratio
+
+    if plot_inducing:
+        ax.plot(model.Z[:, input_1], model.Z[:, input_2], '^w')
+
+    if updates:
+        ax.figure.canvas.show()
+        raw_input('Enter to continue')
+
+    pb.title('Magnification Factor')
+    return ax