MRD updates merge

- termination rule changes for SCG
- oil flow updates

GPy/examples/dimensionality_reduction.py

@@ -63,7 +63,7 @@ def GPLVM_oil_100(optimize=True):
     m.plot_latent(labels=m.data_labels)
     return m
 
-def BGPLVM_oil(optimize=True, N=100, Q=10, M=15, max_f_eval=300, plot=False):
+def BGPLVM_oil(optimize=True, N=100, Q=10, M=20, max_f_eval=300, plot=False):
     data = GPy.util.datasets.oil()
 
     # create simple GP model
@@ -72,19 +72,19 @@ def BGPLVM_oil(optimize=True, N=100, Q=10, M=15, max_f_eval=300, plot=False):
     m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=kernel, M=M)
     m.data_labels = data['Y'][:N].argmax(axis=1)
 
+    m.constrain('variance', logexp_clipped())
+    m.constrain('length', logexp_clipped())
+    m['lengt'] = 100.
+    m.ensure_default_constraints()
+
     # optimize
     if optimize:
-        m.constrain_fixed('noise', 1. / Y.var())
+        m.unconstrain('noise'); m.constrain_fixed('noise', Y.var() / 100.)
-        m.constrain('variance', logexp_clipped())
+        m.optimize('scg', messages=1, max_f_eval=150)
-        m['lengt'] = 1000
 
-        m.ensure_default_constraints()
-        m.optimize('scg', messages=1, max_f_eval=max(80, max_f_eval))
         m.unconstrain('noise')
-        m.constrain_positive('noise')
+        m.constrain('noise', logexp_clipped())
-        m.optimize('scg', messages=1, max_f_eval=max(0, max_f_eval - 80))
+        m.optimize('scg', messages=1, max_f_eval=max_f_eval)
-    else:
-        m.ensure_default_constraints()
 
     if plot:
         y = m.likelihood.Y[0, :]
@@ -182,7 +182,7 @@ def bgplvm_simulation_matlab_compare():
     Y = sim_data['Y']
     S = sim_data['S']
     mu = sim_data['mu']
-    M, [_, Q] = 20, mu.shape
+    M, [_, Q] = 3, mu.shape
 
     from GPy.models import mrd
     from GPy import kern
@@ -192,7 +192,7 @@ def bgplvm_simulation_matlab_compare():
     m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k,
 #                        X=mu,
 #                        X_variance=S,
-                       _debug=True)
+                       _debug=False)
     m.ensure_default_constraints()
     m.auto_scale_factor = True
     m['noise'] = Y.var() / 100.
@@ -371,7 +371,7 @@ def brendan_faces():
 
     # optimize
     m.ensure_default_constraints()
-    # m.optimize(messages=1, max_f_eval=10000)
+    m.optimize(messages=1, max_f_eval=10000)
 
     ax = m.plot_latent()
     y = m.likelihood.Y[0, :]

GPy/inference/SCG.py

@@ -49,6 +49,7 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
     function_eval = 1
     fnow = fold
     gradnew = gradf(x, *optargs)  # Initial gradient.
+    current_grad = np.dot(gradnew, gradnew)
     gradold = gradnew.copy()
     d = -gradnew  # Initial search direction.
     success = True  # Force calculation of directional derivs.
@@ -110,7 +111,7 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
         iteration += 1
         if display:
             print '\r',
-            print 'i: {0:>5g}  f:{1:> 12e}  b:{2:> 12e} |g|:{3:> 12e}'.format(iteration, fnow, beta, np.dot(gradnew, gradnew)),
+            print 'i: {0:>5g}  f:{1:> 12e}  b:{2:> 12e} |g|:{3:> 12e}'.format(iteration, fnow, beta, current_grad),
             # print 'Iteration:', iteration, ' Objective:', fnow, '  Scale:', beta, '\r',
             sys.stdout.flush()
 
@@ -125,8 +126,9 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
                 fold = fnew
                 gradold = gradnew
                 gradnew = gradf(x, *optargs)
+                current_grad = np.dot(gradnew, gradnew)
                 # If the gradient is zero then we are done.
-                if (np.dot(gradnew, gradnew)) <= gtol:
+                if current_grad <= gtol:
                     status = 'converged'
                     return x, flog, function_eval, status
 

GPy/models/Bayesian_GPLVM.py

@@ -136,14 +136,16 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
                 self._savedparams.append([self.f_call, self._get_params()])
                 self._savedgradients.append([self.f_call, self._log_likelihood_gradients()])
                 self._savedpsiKmm.append([self.f_call, [self.Kmm, self.dL_dKmm]])
-                sf2 = self.scale_factor ** 2
+#                 sf2 = self.scale_factor ** 2
                 if self.likelihood.is_heteroscedastic:
                     A = -0.5 * self.N * self.D * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.V * self.likelihood.Y)
-                    B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A) * sf2)
+#                     B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A) * sf2)
+                    B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A))
                 else:
                     A = -0.5 * self.N * self.D * (np.log(2.*np.pi) + np.log(self.likelihood._variance)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
-                    B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A) * sf2)
+#                     B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A) * sf2)
-                C = -self.D * (np.sum(np.log(np.diag(self.LB))) + 0.5 * self.M * np.log(sf2))
+                    B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A))
+                C = -self.D * (np.sum(np.log(np.diag(self.LB)))) # + 0.5 * self.M * np.log(sf2))
                 D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V))
                 self._savedABCD.append([self.f_call, A, B, C, D])
 

GPy/models/mrd.py

@@ -273,7 +273,7 @@ class MRD(model):
 
     def _handle_plotting(self, fig_num, axes, plotf):
         if axes is None:
-            fig = pylab.figure(num=fig_num, figsize=(4 * len(self.bgplvms), 3 * len(self.bgplvms)))
+            fig = pylab.figure(num=fig_num, figsize=(4 * len(self.bgplvms), 3))
         for i, g in enumerate(self.bgplvms):
             if axes is None:
                 ax = fig.add_subplot(1, len(self.bgplvms), i + 1)

GPy/models/sparse_GP.py

@@ -104,7 +104,7 @@ class sparse_GP(GP):
                 tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
                 self.A = tdot(tmp)
             else:
-#                 tmp = self.psi1 * (np.sqrt(self.likelihood.precision) / sf)
+                # tmp = self.psi1 * (np.sqrt(self.likelihood.precision) / sf)
                 tmp = self.psi1 * (np.sqrt(self.likelihood.precision))
                 tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
                 self.A = tdot(tmp)
@@ -163,7 +163,7 @@ class sparse_GP(GP):
         else:
             # likelihood is not heterscedatic
             self.partial_for_likelihood = -0.5 * self.N * self.D * self.likelihood.precision + 0.5 * self.likelihood.trYYT * self.likelihood.precision ** 2
-#             self.partial_for_likelihood += 0.5 * self.D * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision * sf2)
+            #  self.partial_for_likelihood += 0.5 * self.D * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision * sf2)
             self.partial_for_likelihood += 0.5 * self.D * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision)
             self.partial_for_likelihood += self.likelihood.precision * (0.5 * np.sum(self.A * self.DBi_plus_BiPBi) - np.sum(np.square(self._LBi_Lmi_psi1V)))
 
@@ -177,7 +177,7 @@ class sparse_GP(GP):
 #             B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A) * sf2)
             B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A))
         else:
-            A = -0.5 * self.N * self.D * (np.log(2.*np.pi) + np.log(self.likelihood._variance)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
+            A = -0.5 * self.N * self.D * (np.log(2.*np.pi) - np.log(self.likelihood.precision)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
 #             B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A) * sf2)
             B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A))
 #         C = -self.D * (np.sum(np.log(np.diag(self.LB))) + 0.5 * self.M * np.log(sf2))