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

GPy/examples/dimensionality_reduction.py

@@ -131,7 +131,7 @@ def BGPLVM_oil(optimize=True, N=100, Q=5, M=25, max_f_eval=4e3, plot=False, **k)
     m = GPy.models.Bayesian_GPLVM(Yn, Q, kernel=kernel, M=M, **k)
     m.data_labels = data['Y'][:N].argmax(axis=1)
 
-    m.constrain('variance|leng', logexp_clipped())
+    # m.constrain('variance|leng', logexp_clipped())
     m['lengt'] = m.X.var(0).max() / m.X.var(0)
     m['noise'] = Yn.var() / 100.
 
@@ -246,7 +246,7 @@ def bgplvm_simulation_matlab_compare():
 def bgplvm_simulation(optimize='scg',
                       plot=True,
                       max_f_eval=2e4):
-    from GPy.core.transformations import logexp_clipped
+#     from GPy.core.transformations import logexp_clipped
     D1, D2, D3, N, M, Q = 15, 8, 8, 100, 3, 5
     slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot)
 
@@ -259,8 +259,8 @@ def bgplvm_simulation(optimize='scg',
 
     k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q)
     m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k, _debug=True)
-    m.constrain('variance|noise', logexp_clipped())
+    # m.constrain('variance|noise', logexp_clipped())
-#     m.ensure_default_constraints()
+    m.ensure_default_constraints()
     m['noise'] = Y.var() / 100.
     m['linear_variance'] = .01
 
@@ -273,8 +273,8 @@ def bgplvm_simulation(optimize='scg',
         pylab.figure(); pylab.axis(); m.kern.plot_ARD()
     return m
 
-def mrd_simulation(optimize=True, plot_sim=False):
+def mrd_simulation(optimize=True, plot_sim=False, **kw):
-    D1, D2, D3, N, M, Q = 150, 250, 30, 300, 3, 7
+    D1, D2, D3, N, M, Q = 150, 250, 30, 200, 3, 7
     slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
 
     from GPy.models import mrd
@@ -284,12 +284,12 @@ def mrd_simulation(optimize=True, plot_sim=False):
     reload(mrd); reload(kern)
 
     k = kern.linear(Q, [0.01] * Q, True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
-    m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, initx="concat", initz='permute')
+    m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, initx="concat", initz='permute', **kw)
 
     for i, Y in enumerate(Ylist):
         m['{}_noise'.format(i + 1)] = Y.var() / 100.
 
-    m.constrain('variance|noise', logexp_clipped())
+    # m.constrain('variance|noise', logexp_clipped())
     m.ensure_default_constraints()
 
     # DEBUG
@@ -305,12 +305,13 @@ def brendan_faces():
     from GPy import kern
     data = GPy.util.datasets.brendan_faces()
     Q = 2
-    # Y = data['Y'][0:-1:2, :]
+    Y = data['Y'][0:-1:10, :]
-    Y = data['Y']
+    # Y = data['Y']
     Yn = Y - Y.mean()
     Yn /= Yn.std()
 
-    m = GPy.models.GPLVM(Yn, Q)#, M=Y.shape[0]/4)
+    m = GPy.models.GPLVM(Yn, Q)
+    # m = GPy.models.Bayesian_GPLVM(Yn, Q, M=100)
 
     # optimize
     m.constrain('rbf|noise|white', GPy.core.transformations.logexp_clipped())
@@ -318,7 +319,7 @@ def brendan_faces():
     m.ensure_default_constraints()
     m.optimize('scg', messages=1, max_f_eval=10000)
 
-    ax = m.plot_latent()
+    ax = m.plot_latent(which_indices=(0,1))
     y = m.likelihood.Y[0, :]
     data_show = GPy.util.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, invert=False, scale=False)
     lvm_visualizer = GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)

GPy/inference/SCG.py

@@ -39,11 +39,6 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
     function_eval number of fn evaluations
     status: string describing convergence status
     """
-    
-    if display:
-        print "  SCG"
-        print ' {0:{mi}s}   {1:11s}    {2:11s}    {3:11s}'.format("I", "F", "Scale", "|g|", mi=len(str(maxiters)))
-
     if xtol is None:
         xtol = 1e-6
     if ftol is None:
@@ -69,6 +64,9 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
 
     iteration = 0
 
+    if display:
+        print ' {0:{mi}s}   {1:11s}    {2:11s}    {3:11s}'.format("I", "F", "Scale", "|g|", mi=len(str(maxiters)))
+
     # Main optimization loop.
     while iteration < maxiters:
 
@@ -129,10 +127,10 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
 
             else:
                 # Update variables for new position
-                fold = fnew
-                gradold = gradnew
                 gradnew = gradf(x, *optargs)
                 current_grad = np.dot(gradnew, gradnew)
+                gradold = gradnew
+                fold = fnew
                 # If the gradient is zero then we are done.
                 if current_grad <= gtol:
                     status = 'converged'

GPy/likelihoods/Gaussian.py

@@ -53,7 +53,7 @@ class Gaussian(likelihood):
     def _set_params(self, x):
         x = float(x)
         if self._variance != x:
-            self.precision = 1. / x
+            self.precision = 1. / max(x, 1e-6)
             self.covariance_matrix = np.eye(self.N) * x
             self.V = (self.precision) * self.Y
             self._variance = x

GPy/models/Bayesian_GPLVM.py

@@ -14,6 +14,7 @@ import itertools
 from matplotlib.colors import colorConverter
 from matplotlib.figure import SubplotParams
 from GPy.inference.optimization import SCG
+from GPy.util import plot_latent
 
 class Bayesian_GPLVM(sparse_GP, GPLVM):
     """
@@ -37,6 +38,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
 
         if X == None:
             X = self.initialise_latent(init, Q, likelihood.Y)
+        self.init = init
 
         if X_variance is None:
             X_variance = np.clip((np.ones_like(X) * 0.5) + .01 * np.random.randn(*X.shape), 0.001, 1)
@@ -177,18 +179,8 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         self.dbound_dZtheta = sparse_GP._log_likelihood_gradients(self)
         return np.hstack((self.dbound_dmuS.flatten(), self.dbound_dZtheta))
 
-    def plot_latent(self, which_indices=None, *args, **kwargs):
+    def plot_latent(self, *args, **kwargs):
-
+        plot_latent.plot_latent_indices(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
-        ax = GPLVM.plot_latent(self, which_indices=[input_1, input_2], *args, **kwargs)
-        ax.plot(self.Z[:, input_1], self.Z[:, input_2], '^w')
-        return ax
 
     def do_test_latents(self, Y):
         """
@@ -262,6 +254,14 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
         return fig
 
+    def __getstate__(self):
+        return (self.likelihood, self.Q, self.X, self.X_variance,
+                self.init, self.M, self.Z, self.kern,
+                self.oldpsave, self._debug)
+
+    def __setstate__(self, state):
+        self.__init__(*state)
+
     def _debug_filter_params(self, x):
         start, end = 0, self.X.size,
         X = x[start:end].reshape(self.N, self.Q)

GPy/models/GPLVM.py

@@ -11,6 +11,8 @@ from ..util.linalg import pdinv, PCA
 from GP import GP
 from ..likelihoods import Gaussian
 from .. import util
+from GPy.util import plot_latent
+
 
 class GPLVM(GP):
     """
@@ -60,65 +62,5 @@ class GPLVM(GP):
         mu, var, upper, lower = self.predict(Xnew)
         pb.plot(mu[:,0], mu[:,1],'k',linewidth=1.5)
 
-    def plot_latent(self, labels=None, which_indices=None, resolution=50, ax=None):
+    def plot_latent(self, *args, **kwargs):
-        """
+        util.plot_latent.plot_latent(self, *args, **kwargs)
-        :param labels: a np.array of size self.N 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:
-            ax = pb.gca()
-        util.plot.Tango.reset()
-
-        if labels is None:
-            labels = np.ones(self.N)
-        if which_indices is None:
-            if self.Q==1:
-                input_1 = 0
-                input_2 = None
-            if self.Q==2:
-                input_1, input_2 = 0,1
-            else:
-                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
-
-        #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[:, :2] = Xtest
-        mu, var, low, up = self.predict(Xtest_full)
-        var = var[:, :1]
-        ax.imshow(var.reshape(resolution, resolution).T,
-                  extent=[xmin[0], xmax[0], xmin[1], xmax[1]], cmap=pb.cm.binary,interpolation='bilinear',origin='lower')
-
-        for i,ul in enumerate(np.unique(labels)):
-            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
-
-            index = np.nonzero(labels==ul)[0]
-            if self.Q==1:
-                x = self.X[index,input_1]
-                y = np.zeros(index.size)
-            else:
-                x = self.X[index,input_1]
-                y = self.X[index,input_2]
-            ax.plot(x,y,marker='o',color=util.plot.Tango.nextMedium(),mew=0,label=this_label,linewidth=0)
-
-        ax.set_xlabel('latent dimension %i'%input_1)
-        ax.set_ylabel('latent dimension %i'%input_2)
-
-        if not np.all(labels==1.):
-            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
-        return ax

GPy/testing/__init__.py

@@ -7,6 +7,6 @@ import unittest
 import sys
 
 def deepTest(reason):
-    if 'deep' in sys.argv:
+    if reason:
         return lambda x:x
     return unittest.skip("Not deep scanning, enable deepscan by adding 'deep' argument")

GPy/testing/psi_stat_expactation_tests.py

@@ -6,10 +6,9 @@ Created on 26 Apr 2013
 import unittest
 import GPy
 import numpy as np
-import sys
+from GPy import testing
-from .. import testing
 
-__test__ = True
+__test__ = False
 np.random.seed(0)
 
 def ard(p):
@@ -20,7 +19,7 @@ def ard(p):
         pass
     return ""
 
-@testing.deepTest
+@testing.deepTest(__test__)
 class Test(unittest.TestCase):
     D = 9
     M = 4
@@ -29,13 +28,22 @@ class Test(unittest.TestCase):
 
     def setUp(self):
         self.kerns = (
-                      GPy.kern.rbf(self.D), GPy.kern.rbf(self.D, ARD=True),
+#                       (GPy.kern.rbf(self.D, ARD=True) +
-                      GPy.kern.linear(self.D, ARD=False), GPy.kern.linear(self.D, ARD=True),
+#                        GPy.kern.linear(self.D, ARD=True) +
-                      GPy.kern.linear(self.D) + GPy.kern.bias(self.D),
+#                        GPy.kern.bias(self.D) +
-                      GPy.kern.rbf(self.D) + GPy.kern.bias(self.D),
+#                        GPy.kern.white(self.D)),
-                      GPy.kern.linear(self.D) + GPy.kern.bias(self.D) + GPy.kern.white(self.D),
+                      (GPy.kern.rbf(self.D, np.random.rand(), np.random.rand(self.D), ARD=True) +
-                      GPy.kern.rbf(self.D) + GPy.kern.bias(self.D) + GPy.kern.white(self.D),
+                       GPy.kern.rbf(self.D, np.random.rand(), np.random.rand(self.D), ARD=True) +
-                      GPy.kern.bias(self.D), GPy.kern.white(self.D),
+                       GPy.kern.linear(self.D, np.random.rand(self.D), ARD=True) +
+                       GPy.kern.bias(self.D) +
+                       GPy.kern.white(self.D)),
+#                       GPy.kern.rbf(self.D), GPy.kern.rbf(self.D, ARD=True),
+#                       GPy.kern.linear(self.D, ARD=False), GPy.kern.linear(self.D, ARD=True),
+#                       GPy.kern.linear(self.D) + GPy.kern.bias(self.D),
+#                       GPy.kern.rbf(self.D) + GPy.kern.bias(self.D),
+#                       GPy.kern.linear(self.D) + GPy.kern.bias(self.D) + GPy.kern.white(self.D),
+#                       GPy.kern.rbf(self.D) + GPy.kern.bias(self.D) + GPy.kern.white(self.D),
+#                       GPy.kern.bias(self.D), GPy.kern.white(self.D),
                       )
         self.q_x_mean = np.random.randn(self.D)
         self.q_x_variance = np.exp(np.random.randn(self.D))
@@ -64,8 +72,9 @@ class Test(unittest.TestCase):
             K_ /= self.Nsamples / Nsamples
             msg = "psi1: " + "+".join([p.name + ard(p) for p in kern.parts])
             try:
-#                 pylab.figure(msg)
+                import pylab
-#                 pylab.plot(diffs)
+                pylab.figure(msg)
+                pylab.plot(diffs)
                 self.assertTrue(np.allclose(psi1.squeeze(), K_,
                                             rtol=1e-1, atol=.1),
                                 msg=msg + ": not matching")
@@ -90,8 +99,9 @@ class Test(unittest.TestCase):
             K_ /= self.Nsamples / Nsamples
             msg = "psi2: {}".format("+".join([p.name + ard(p) for p in kern.parts]))
             try:
-#                 pylab.figure(msg)
+                import pylab
-#                 pylab.plot(diffs)
+                pylab.figure(msg)
+                pylab.plot(diffs)
                 self.assertTrue(np.allclose(psi2.squeeze(), K_,
                                             rtol=1e-1, atol=.1),
                                 msg=msg + ": not matching")
@@ -104,7 +114,9 @@ class Test(unittest.TestCase):
                 pass
 
 if __name__ == "__main__":
-    import sys;sys.argv = ['',
+    import sys
+    __test__ = 'deep' in sys.argv
+    sys.argv = ['',
          'Test.test_psi0',
          'Test.test_psi1',
          'Test.test_psi2',

GPy/util/mocap.py

@@ -1,6 +1,8 @@
 import os
 import numpy as np
 import math
+from GPy.util import datasets as dat
+import urllib2
 
 class vertex:
     def __init__(self, name, id, parents=[], children=[], meta = {}):
@@ -687,3 +689,71 @@ def read_connections(file_name, point_names):
     
   
 skel = acclaim_skeleton()
+
+
+
+    
+def fetch_data(base_url = 'http://mocap.cs.cmu.edu:8080/subjects', skel_store_dir = '.', motion_store_dir = '.', subj_motions = None, store_motions = True, return_motions = True, messages = True):
+    ''' 
+    Download and store the skel. and motions indicated in a tuple (A,B) where A is a list of skeletons and B
+    the corresponding 2-D list of motions, ie B_ij is the j-th motion to download for skeleton A_i
+    The method can optionally store the fetched data and / or return them as arrays.
+    If the data are already stored, they are not fetched but just retrieved.
+
+    e.g.
+    # Download the data, do not return anything
+    GPy.util.mocap.fetch_data(subj_motions = (['35'],[['01','02','03']]), return_motions = False)
+    # Fetch and return the data in a list. Do not store them anywhere
+    GPy.util.mocap.fetch_data(subj_motions = (['35'],[['01','02','03']]), return_motions = True, store_motions = False)
+
+    In both cases above, if the data do exist in the given skel_store_dir and motion_store_dir, they are just loaded from there.
+    '''
+    
+    subjects = subj_motions[0]
+    motions = subj_motions[1]
+    all_skels = []
+    
+    assert len(subjects) == len(motions)
+    
+    if return_motions:
+        all_motions = [list() for _ in range(len(subjects))]
+    else:
+        all_motions = []
+            
+    for i in range(len(subjects)):
+        cur_skel_suffix = '/' + subjects[i] + '/' 
+        cur_skel_dir = skel_store_dir + cur_skel_suffix
+        cur_skel_file = cur_skel_dir + subjects[i] + '.asf'
+        cur_skel_url = base_url + cur_skel_suffix + subjects[i] + '.asf'
+        
+        if os.path.isfile(cur_skel_file):
+            if return_motions:
+                with open(cur_skel_file, 'r') as f:
+                    cur_skel_data = f.read()
+        else:
+            if store_motions:
+                if not os.path.isdir(cur_skel_dir):
+                    os.mkdir(cur_skel_dir)
+                if not os.path.isdir(motion_store_dir + cur_skel_suffix):
+                    os.mkdir(motion_store_dir + cur_skel_suffix)
+            cur_skel_data = dat.fetch_dataset(cur_skel_url, cur_skel_file, store_motions, messages)
+        
+        if return_motions:
+            all_skels.append(cur_skel_data)
+        
+        for j in range(len(motions[i])):
+            cur_motion_url = base_url + cur_skel_suffix + subjects[i] + '_' + motions[i][j] + '.amc'
+            cur_motion_file = motion_store_dir + cur_skel_suffix  + subjects[i] + '_' + motions[i][j] + '.amc'
+            if os.path.isfile(cur_motion_file):
+                with open(cur_motion_file, 'r') as f:
+                    if return_motions:
+                        cur_motion_data = f.read()
+            else:
+                cur_motion_data = dat.fetch_dataset(cur_motion_url, cur_motion_file, store_motions, messages)
+
+            if return_motions:
+                all_motions[i].append(cur_motion_data)
+
+    return all_skels, all_motions
+
+

GPy/util/mocap_fetch.py

@@ -1,13 +0,0 @@
-import GPy
-import urllib2
-
-# TODO...
-class mocap_fetch(base_url = 'http://mocap.cs.cmu.edu:8080/subjects/', skel_store_dir = './', motion_store_dir = './'):
-    def __init__(self):
-        self.base_url = base_url
-        self.store_dir = store_dir
-        self.motion_dict = []
-    
-    def fetch_motions(self, motion_dict = None):
-        response = urllib2.urlopen(...)
-        html = response.read()

GPy/util/plot_latent.py

@@ -0,0 +1,91 @@
+import pylab as pb
+import numpy as np
+from .. import util
+
+def plot_latent(model, labels=None, which_indices=None, resolution=50, ax=None, marker='o', s=40):
+    """
+    :param labels: a np.array of size model.N 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:
+        ax = pb.gca()
+    util.plot.Tango.reset()
+
+    if labels is None:
+        labels = np.ones(model.N)
+    if which_indices is None:
+        if model.Q==1:
+            input_1 = 0
+            input_2 = None
+        if model.Q==2:
+            input_1, input_2 = 0,1
+        else:
+            try:
+                input_1, input_2 = np.argsort(model.input_sensitivity())[:2]
+            except:
+                raise ValueError, "cannot Atomatically determine which dimensions to plot, please pass 'which_indices'"
+    else:
+        input_1, input_2 = 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]))
+    Xtest_full[:, :2] = Xtest
+    mu, var, low, up = model.predict(Xtest_full)
+    var = var[:, :1]
+    ax.imshow(var.reshape(resolution, resolution).T,
+              extent=[xmin[0], xmax[0], xmin[1], xmax[1]], cmap=pb.cm.binary,interpolation='bilinear',origin='lower')
+
+    # make sure labels are in order of input:
+    ulabels = []
+    for lab in labels:
+        if not lab in ulabels:
+            ulabels.append(lab)
+
+    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
+        if len(marker) == len(ulabels):
+            m = marker[i]
+        else:
+            m = marker
+
+        index = np.nonzero(labels==ul)[0]
+        if model.Q==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.):
+        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
+    return ax
+
+
+def plot_latent_indices(model, which_indices=None, *args, **kwargs):
+
+    if which_indices is None:
+        try:
+            input_1, input_2 = np.argsort(model.input_sensitivity())[:2]
+        except:
+            raise ValueError, "cannot Automatically determine which dimensions to plot, please pass 'which_indices'"
+    else:
+        input_1, input_2 = which_indices
+    ax = plot_latent(model, which_indices=[input_1, input_2], *args, **kwargs)
+    # TODO: Here test if there are inducing points...
+    ax.plot(model.Z[:, input_1], model.Z[:, input_2], '^w')
+    return ax