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

.gitignore

@@ -39,3 +39,11 @@ nosetests.xml
 
 #bfgs optimiser leaves this lying around
 iterate.dat
+
+# Nosetests #
+#############
+*.noseids
+
+# git merge files #
+###################
+*.orig

GPy/core/parameterised.py

@@ -251,7 +251,18 @@ class parameterised(object):
 
     def _set_params_transformed(self, x):
         """ takes the vector x, which is then modified (by untying, reparameterising or inserting fixed values), and then call self._set_params"""
+        self._set_params(self._untransform_params(x))
 
+    def _untransform_params(self,x):
+        """
+        The transformation required for _set_params_transformed.
+
+        This moves the vector x seen by the optimiser (unconstrained) to the
+        valid parameter vector seen by the model
+
+        Note:
+          - This function is separate from _set_params_transformed for downstream flexibility
+        """
         # work out how many places are fixed, and where they are. tricky logic!
         fix_places = self.fixed_indices + [t[1:] for t in self.tied_indices]
         if len(fix_places):
@@ -272,7 +283,8 @@ class parameterised(object):
         [np.put(xx,i,t.f(xx[i])) for i,t in zip(self.constrained_indices, self.constraints)]
         if hasattr(self,'debug'):
             stop
-        self._set_params(xx)
+
+        return xx
 
     def _get_param_names_transformed(self):
         """

GPy/core/transformations.py

@@ -39,23 +39,29 @@ class logexp(transformation):
         return '(+ve)'
 
 class logexp_clipped(transformation):
-    def __init__(self):
+    max_bound = 1e300
+    min_bound = 1e-10
+    log_max_bound = np.log(max_bound)
+    log_min_bound = np.log(min_bound)
+    def __init__(self, lower=1e-6):
         self.domain = 'positive'
+        self.lower = lower
     def f(self, x):
-        f = np.log(1. + np.exp(x))
+        exp = np.exp(np.clip(x, self.log_min_bound, self.log_max_bound))
+        f = np.log(1. + exp)
         return f
     def finv(self, f):
-        return np.log(np.exp(f) - 1.)
+        return np.log(np.exp(np.clip(f, self.min_bound, self.max_bound)) - 1.)
     def gradfactor(self, f):
         ef = np.exp(f)
         gf = (ef - 1.) / ef
-        return np.where(f < 1e-6, 0, gf)
+        return np.where(f < self.lower, 0, gf)
-    def initialize(self,f):
+    def initialize(self, f):
-        if np.any(f<0.):
+        if np.any(f < 0.):
             print "Warning: changing parameters to satisfy constraints"
         return np.abs(f)
     def __str__(self):
-        return '(+ve)'
+        return '(+ve_c)'
 
 class exponent(transformation):
     def __init__(self):
@@ -67,7 +73,7 @@ class exponent(transformation):
     def gradfactor(self, f):
         return f
     def initialize(self, f):
-        if np.any(f<0.):
+        if np.any(f < 0.):
             print "Warning: changing parameters to satisfy constraints"
         return np.abs(f)
     def __str__(self):
@@ -83,7 +89,7 @@ class negative_exponent(transformation):
     def gradfactor(self, f):
         return f
     def initialize(self, f):
-        if np.any(f>0.):
+        if np.any(f > 0.):
             print "Warning: changing parameters to satisfy constraints"
         return -np.abs(f)
     def __str__(self):
@@ -114,11 +120,11 @@ class logistic(transformation):
     def finv(self, f):
         return np.log(np.clip(f - self.lower, 1e-10, np.inf) / np.clip(self.upper - f, 1e-10, np.inf))
     def gradfactor(self, f):
-        return (f-self.lower)*(self.upper-f)/self.difference
+        return (f - self.lower) * (self.upper - f) / self.difference
-    def initialize(self,f):
+    def initialize(self, f):
-        if np.any(np.logical_or(f<self.lower,f>self.upper)):
+        if np.any(np.logical_or(f < self.lower, f > self.upper)):
             print "Warning: changing parameters to satisfy constraints"
-        return np.where(np.logical_or(f<self.lower,f>self.upper),self.f(f*0.),f)
+        return np.where(np.logical_or(f < self.lower, f > self.upper), self.f(f * 0.), f)
     def __str__(self):
         return '({},{})'.format(self.lower, self.upper)
 

GPy/examples/dimensionality_reduction.py

@@ -2,13 +2,11 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
-import pylab as pb
+from matplotlib import pyplot as plt
-from matplotlib import pyplot as plt, pyplot
 
 import GPy
 from GPy.models.Bayesian_GPLVM import Bayesian_GPLVM
-from GPy.util.datasets import simulation_BGPLVM
+from GPy.util.datasets import swiss_roll_generated
-from GPy.core.transformations import square, logexp_clipped
 
 default_seed = np.random.seed(123344)
 
@@ -47,10 +45,11 @@ def BGPLVM(seed=default_seed):
 
 def GPLVM_oil_100(optimize=True):
     data = GPy.util.datasets.oil_100()
+    Y = data['X']
 
     # create simple GP model
     kernel = GPy.kern.rbf(6, ARD=True) + GPy.kern.bias(6)
-    m = GPy.models.GPLVM(data['X'], 6, kernel=kernel)
+    m = GPy.models.GPLVM(Y, 6, kernel=kernel)
     m.data_labels = data['Y'].argmax(axis=1)
 
     # optimize
@@ -63,27 +62,83 @@ 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=20, max_f_eval=300, plot=False):
+def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
+    from GPy.util.datasets import swiss_roll
+    from GPy.core.transformations import logexp_clipped
+
+    data = swiss_roll_generated(N=N, sigma=sigma)
+    Y = data['Y']
+    Y -= Y.mean()
+    Y /= Y.std()
+
+    t = data['t']
+    c = data['colors']
+
+    try:
+        from sklearn.manifold.isomap import Isomap
+        iso = Isomap().fit(Y)
+        X = iso.embedding_
+        if Q > 2:
+            X = np.hstack((X, np.random.randn(N, Q - 2)))
+    except ImportError:
+        X = np.random.randn(N, Q)
+
+    if plot:
+        from mpl_toolkits import mplot3d
+        import pylab
+        fig = pylab.figure("Swiss Roll Data")
+        ax = fig.add_subplot(121, projection='3d')
+        ax.scatter(*Y.T, c=c)
+        ax.set_title("Swiss Roll")
+
+        ax = fig.add_subplot(122)
+        ax.scatter(*X.T[:2], c=c)
+        ax.set_title("Initialization")
+
+
+    var = .5
+    S = (var * np.ones_like(X) + np.clip(np.random.randn(N, Q) * var ** 2,
+                                         - (1 - var),
+                                         (1 - var))) + .001
+    Z = np.random.permutation(X)[:M]
+
+    kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
+
+    m = Bayesian_GPLVM(Y, Q, X=X, X_variance=S, M=M, Z=Z, kernel=kernel)
+    m.data_colors = c
+    m.data_t = t
+
+    m.constrain('variance|length', logexp_clipped())
+    m['lengthscale'] = 1. # X.var(0).max() / X.var(0)
+    m['noise'] = Y.var() / 100.
+    m.ensure_default_constraints()
+
+    if optimize:
+        m.optimize('scg', messages=1)
+    return m
+
+def BGPLVM_oil(optimize=True, N=100, Q=5, M=25, max_f_eval=4e3, plot=False, **k):
+    np.random.seed(0)
     data = GPy.util.datasets.oil()
+    from GPy.core.transformations import logexp_clipped
 
     # create simple GP model
     kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
     Y = data['X'][:N]
-    m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=kernel, M=M)
+    Yn = Y - Y.mean(0)
+    Yn /= Yn.std(0)
+
+    m = GPy.models.Bayesian_GPLVM(Yn, Q, kernel=kernel, M=M, **k)
     m.data_labels = data['Y'][:N].argmax(axis=1)
 
-    m.constrain('variance', logexp_clipped())
+    m.constrain('variance|leng', logexp_clipped())
-    m.constrain('length', logexp_clipped())
+    m['lengt'] = m.X.var(0).max() / m.X.var(0)
-    m['lengt'] = 100.
+    m['noise'] = Yn.var() / 100.
+
     m.ensure_default_constraints()
 
     # optimize
     if optimize:
-        m.unconstrain('noise'); m.constrain_fixed('noise', Y.var() / 100.)
-        m.optimize('scg', messages=1, max_f_eval=150)
-
-        m.unconstrain('noise')
-        m.constrain('noise', logexp_clipped())
         m.optimize('scg', messages=1, max_f_eval=max_f_eval)
 
     if plot:
@@ -95,11 +150,6 @@ def BGPLVM_oil(optimize=True, N=100, Q=10, M=20, max_f_eval=300, plot=False):
         lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X[0, :].copy(), m, data_show, latent_axes=latent_axes) # , sense_axes=sense_axes)
         raw_input('Press enter to finish')
         plt.close('all')
-    # # plot
-    # print(m)
-    # 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():
@@ -115,6 +165,8 @@ def oil_100():
     # m.plot_latent(labels=data['Y'].argmax(axis=1))
     return m
 
+
+
 def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False):
     x = np.linspace(0, 4 * np.pi, N)[:, None]
     s1 = np.vectorize(lambda x: np.sin(x))
@@ -127,15 +179,6 @@ def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False):
     s3 = s3(x)
     sS = sS(x)
 
-#     s1 -= s1.mean()
-#     s2 -= s2.mean()
-#     s3 -= s3.mean()
-#     sS -= sS.mean()
-#     s1 /= .5 * (np.abs(s1).max() - np.abs(s1).min())
-#     s2 /= .5 * (np.abs(s2).max() - np.abs(s2).min())
-#     s3 /= .5 * (np.abs(s3).max() - np.abs(s3).min())
-#     sS /= .5 * (np.abs(sS).max() - np.abs(sS).min())
-
     S1 = np.hstack([s1, sS])
     S2 = np.hstack([s2, sS])
     S3 = np.hstack([s3, sS])
@@ -155,16 +198,17 @@ def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False):
     Y2 /= Y2.std(0)
     Y3 /= Y3.std(0)
 
-    slist = [s1, s2, s3, sS]
+    slist = [sS, s1, s2, s3]
+    slist_names = ["sS", "s1", "s2", "s3"]
     Ylist = [Y1, Y2, Y3]
 
     if plot_sim:
         import pylab
         import itertools
-        fig = pylab.figure("MRD Simulation", figsize=(8, 6))
+        fig = pylab.figure("MRD Simulation Data", figsize=(8, 6))
         fig.clf()
         ax = fig.add_subplot(2, 1, 1)
-        labls = sorted(filter(lambda x: x.startswith("s"), locals()))
+        labls = slist_names
         for S, lab in itertools.izip(slist, labls):
             ax.plot(S, label=lab)
         ax.legend()
@@ -178,6 +222,7 @@ def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False):
     return slist, [S1, S2, S3], Ylist
 
 def bgplvm_simulation_matlab_compare():
+    from GPy.util.datasets import simulation_BGPLVM
     sim_data = simulation_BGPLVM()
     Y = sim_data['Y']
     S = sim_data['S']
@@ -187,7 +232,6 @@ def bgplvm_simulation_matlab_compare():
     from GPy.models import mrd
     from GPy import kern
     reload(mrd); reload(kern)
-    # k = kern.rbf(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
     k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
     m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k,
 #                        X=mu,
@@ -197,24 +241,14 @@ def bgplvm_simulation_matlab_compare():
     m.auto_scale_factor = True
     m['noise'] = Y.var() / 100.
     m['linear_variance'] = .01
-
-#     lscstr = 'X_variance'
-#     m[lscstr] = .01
-#     m.unconstrain(lscstr); m.constrain_fixed(lscstr, .1)
-
-#     cstr = 'white'
-#     m.unconstrain(cstr); m.constrain_bounded(cstr, .01, 1.)
-
-#     cstr = 'noise'
-#     m.unconstrain(cstr); m.constrain_bounded(cstr, .01, 1.)
     return m
 
-def bgplvm_simulation(burnin='scg', plot_sim=False,
+def bgplvm_simulation(optimize='scg',
-                      max_burnin=100, true_X=False,
+                      plot=True,
-                      do_opt=True,
+                      max_f_eval=2e4):
-                      max_f_eval=1000):
+    from GPy.core.transformations import logexp_clipped
-    D1, D2, D3, N, M, Q = 15, 8, 8, 350, 3, 6
+    D1, D2, D3, N, M, Q = 15, 8, 8, 100, 3, 5
-    slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
+    slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot)
 
     from GPy.models import mrd
     from GPy import kern
@@ -224,154 +258,60 @@ def bgplvm_simulation(burnin='scg', plot_sim=False,
     Y = Ylist[0]
 
     k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q)
-#     k = kern.white(Q, .00001) + kern.bias(Q)
     m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k, _debug=True)
-    # m.set('noise',)
+    m.constrain('variance|noise', logexp_clipped())
-    m.constrain('variance', logexp_clipped())
+#     m.ensure_default_constraints()
-
-    m.ensure_default_constraints()
     m['noise'] = Y.var() / 100.
-    m['linear_variance'] = .001
+    m['linear_variance'] = .01
-#     m.auto_scale_factor = True
-#     m.scale_factor = 1.
 
-
+    if optimize:
-    if burnin:
+        print "Optimizing model:"
-        print "initializing beta"
+        m.optimize('scg', max_iters=max_f_eval, max_f_eval=max_f_eval, messages=True)
-        cstr = "noise"
+    if plot:
-        m.unconstrain(cstr); m.constrain_fixed(cstr, Y.var() / 70.)
+        import pylab
-        m.optimize(burnin, messages=1, max_f_eval=max_burnin)
+        m.plot_X_1d()
-
+        pylab.figure(); pylab.axis(); m.kern.plot_ARD()
-        print "releasing beta"
-        cstr = "noise"
-        m.unconstrain(cstr);  m.constrain_positive(cstr)
-
-    if true_X:
-        true_X = np.hstack((slist[0], slist[3], 0. * np.ones((N, Q - 2))))
-        m.set('X_\d', true_X)
-        m.constrain_fixed("X_\d")
-
-        cstr = 'X_variance'
-#         m.unconstrain(cstr), m.constrain_fixed(cstr, .0001)
-        m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-7, .1)
-
-#     cstr = 'X_variance'
-#     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-3, 1.)
-
-    # m['X_var'] = np.ones(N * Q) * .5 + np.random.randn(N * Q) * .01
-
-#     cstr = "iip"
-#     m.unconstrain(cstr); m.constrain_fixed(cstr)
-
-#     cstr = 'variance'
-#     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-10, 1.)
-#     cstr = 'X_\d'
-#     m.unconstrain(cstr), m.constrain_bounded(cstr, -10., 10.)
-#
-#     cstr = 'noise'
-#     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-5, 1.)
-#
-#     cstr = 'white'
-#     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-6, 1.)
-#
-#     cstr = 'linear_variance'
-#     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-10, 10.)
-
-#     cstr = 'variance'
-#     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-10, 10.)
-
-#     np.seterr(all='call')
-#     def ipdbonerr(errtype, flags):
-#         import ipdb; ipdb.set_trace()
-#     np.seterrcall(ipdbonerr)
-
-    if do_opt and burnin:
-        try:
-            m.optimize(burnin, messages=1, max_f_eval=max_f_eval)
-        except:
-            pass
-        finally:
-            return m
     return m
 
 def mrd_simulation(plot_sim=False):
-    # num = 2
-#     ard1 = np.array([1., 1, 0, 0], dtype=float)
-#     ard2 = np.array([0., 1, 1, 0], dtype=float)
-#     ard1[ard1 == 0] = 1E-10
-#     ard2[ard2 == 0] = 1E-10
-
-#     ard1i = 1. / ard1
-#     ard2i = 1. / ard2
-
-#     k = GPy.kern.rbf(Q, ARD=True, lengthscale=ard1i) + GPy.kern.bias(Q, 0) + GPy.kern.white(Q, 0.0001)
-#     Y1 = np.random.multivariate_normal(np.zeros(N), k.K(X), D1).T
-#     Y1 -= Y1.mean(0)
-#
-#     k = GPy.kern.rbf(Q, ARD=True, lengthscale=ard2i) + GPy.kern.bias(Q, 0) + GPy.kern.white(Q, 0.0001)
-#     Y2 = np.random.multivariate_normal(np.zeros(N), k.K(X), D2).T
-#     Y2 -= Y2.mean(0)
-#     make_params = lambda ard: np.hstack([[1], ard, [1, .3]])
     D1, D2, D3, N, M, Q = 150, 250, 300, 700, 3, 7
     slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
 
     from GPy.models import mrd
     from GPy import kern
+    from GPy.core.transformations import logexp_clipped
+
     reload(mrd); reload(kern)
 
-#    k = kern.rbf(2, ARD=True) + kern.bias(2) + kern.white(2)
-#     Y1 = np.random.multivariate_normal(np.zeros(N), k.K(S1), D1).T
-#     Y2 = np.random.multivariate_normal(np.zeros(N), k.K(S2), D2).T
-#     Y3 = np.random.multivariate_normal(np.zeros(N), k.K(S3), D3).T
-
-#     Ylist = Ylist[0:2]
-
-    # k = kern.rbf(Q, ARD=True) + kern.bias(Q) + kern.white(Q)
-
     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', _debug=False)
+    m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, initx="concat", initz='permute')
 
     for i, Y in enumerate(Ylist):
         m['{}_noise'.format(i + 1)] = Y.var() / 100.
 
-    m.constrain('variance', logexp_clipped())
+    m.constrain('variance|noise', logexp_clipped())
     m.ensure_default_constraints()
-#     m.auto_scale_factor = True
 
-#     cstr = 'variance'
+    return m
-#     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-12, 1.)
-#
-#     cstr = 'linear_variance'
-#     m.unconstrain(cstr), m.constrain_positive(cstr)
-
-    print "initializing beta"
-    cstr = "noise"
-    m.unconstrain(cstr); m.constrain_fixed(cstr)
-    m.optimize('scg', messages=1, max_f_eval=2e3, gtol=100)
-
-    print "releasing beta"
-    cstr = "noise"
-    m.unconstrain(cstr);  m.constrain(cstr, logexp_clipped())
-
-#     np.seterr(all='call')
-#     def ipdbonerr(errtype, flags):
-#         import ipdb; ipdb.set_trace()
-#     np.seterrcall(ipdbonerr)
-
-    return m # , mtest
-
-def mrd_silhouette():
-
-    pass
 
 def brendan_faces():
+    from GPy import kern
     data = GPy.util.datasets.brendan_faces()
-    Y = data['Y'][0:-1:10, :]
+    Q = 2
-    m = GPy.models.GPLVM(data['Y'], 2)
+    # Y = data['Y'][0:-1:2, :]
+    Y = data['Y']
+    Yn = Y - Y.mean()
+    Yn /= Yn.std()
+
+    m = GPy.models.GPLVM(Yn, Q)#, M=Y.shape[0]/4)
 
     # optimize
+    # m.constrain_fixed('white', 1e-2)
+    # m.constrain_bounded('noise', 1e-6, 10)
+    m.constrain('rbf', GPy.core.transformations.logexp_clipped())
+
     m.ensure_default_constraints()
-    m.optimize(messages=1, max_f_eval=10000)
+    m.optimize('scg', messages=1, max_f_eval=10000)
 
     ax = m.plot_latent()
     y = m.likelihood.Y[0, :]

GPy/inference/SCG.py

@@ -36,7 +36,8 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
     Returns
     x the optimal value for x
     flog : a list of all the objective values
-
+    function_eval number of fn evaluations
+    status: string describing convergence status
     """
     if xtol is None:
         xtol = 1e-6
@@ -111,7 +112,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, current_grad),
+            print 'Iter: {0:>0{mi}g}  Obj:{1:> 12e}  Scale:{2:> 12e}  |g|:{3:> 12e}'.format(iteration, float(fnow), float(beta), float(current_grad), mi=len(str(maxiters))),
             # print 'Iteration:', iteration, ' Objective:', fnow, '  Scale:', beta, '\r',
             sys.stdout.flush()
 
@@ -130,7 +131,8 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
                 # If the gradient is zero then we are done.
                 if current_grad <= gtol:
                     status = 'converged'
-                    return x, flog, function_eval, status
+                    break
+                    # return x, flog, function_eval, status
 
         # Adjust beta according to comparison ratio.
         if Delta < 0.25:
@@ -147,9 +149,11 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
         elif success:
             gamma = np.dot(gradold - gradnew, gradnew) / (mu)
             d = gamma * d - gradnew
+    else:
+        # If we get here, then we haven't terminated in the given number of
+        # iterations.
+        status = "maxiter exceeded"
 
-    # If we get here, then we haven't terminated in the given number of
+    if display:
-    # iterations.
+        print ""
-    status = "maxiter exceeded"
-
     return x, flog, function_eval, status

GPy/inference/natural_gradient_scg.py

@@ -1,146 +0,0 @@
-#Copyright I. Nabney, N.Lawrence and James Hensman (1996 - 2012)
-
-#Scaled Conjuagte Gradients, originally in Matlab as part of the Netlab toolbox by I. Nabney, converted to python N. Lawrence and given a pythonic interface by James Hensman
-
-#      THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
-#      HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
-#      EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
-#      NOT LIMITED TO, THE IMPLIED WARRANTIES OF
-#      MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-#      PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
-#      REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
-#      DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
-#      EXEMPLARY, OR CONSEQUENTIAL DAMAGES
-#      (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
-#      OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
-#      DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
-#      HOWEVER CAUSED AND ON ANY THEORY OF
-#      LIABILITY, WHETHER IN CONTRACT, STRICT
-#      LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
-#      OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-#      OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-#      POSSIBILITY OF SUCH DAMAGE.
-
-
-import numpy as np
-import sys
-
-def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xtol=1e-6, ftol=1e-6):
-    """
-    Optimisation through Scaled Conjugate Gradients (SCG)
-
-    f: the objective function
-    gradf : the gradient function (should return a 1D np.ndarray)
-    x : the initial condition
-
-    Returns
-    x the optimal value for x
-    flog : a list of all the objective values
-
-    """
-
-    sigma0 = 1.0e-4
-    fold = f(x, *optargs)	# Initial function value.
-    function_eval = 1
-    fnow = fold
-    gradnew = gradf(x, *optargs)	# Initial gradient.
-    gradold = gradnew.copy()
-    d = -gradnew				# Initial search direction.
-    success = True				# Force calculation of directional derivs.
-    nsuccess = 0				# nsuccess counts number of successes.
-    beta = 1.0				# Initial scale parameter.
-    betamin = 1.0e-15 			# Lower bound on scale.
-    betamax = 1.0e100			# Upper bound on scale.
-    status = "Not converged"
-
-    flog = [fold]
-
-    iteration = 0
-
-    # Main optimization loop.
-    while iteration < maxiters:
-
-        # Calculate first and second directional derivatives.
-        if success:
-            mu = np.dot(d, gradnew)
-            if mu >= 0:
-                d = -gradnew
-                mu = np.dot(d, gradnew)
-            kappa = np.dot(d, d)
-            sigma = sigma0/np.sqrt(kappa)
-            xplus = x + sigma*d
-            gplus = gradf(xplus, *optargs)
-            theta = np.dot(d, (gplus - gradnew))/sigma
-
-        # Increase effective curvature and evaluate step size alpha.
-        delta = theta + beta*kappa
-        if delta <= 0:
-            delta = beta*kappa
-            beta = beta - theta/kappa
-
-        alpha = - mu/delta
-
-        # Calculate the comparison ratio.
-        xnew = x + alpha*d
-        fnew = f(xnew, *optargs)
-        function_eval += 1
-
-        if function_eval >= max_f_eval:
-            status = "Maximum number of function evaluations exceeded"
-            return x, flog, function_eval, status
-
-        Delta = 2.*(fnew - fold)/(alpha*mu)
-        if Delta  >= 0.:
-            success = True
-            nsuccess += 1
-            x = xnew
-            fnow = fnew
-        else:
-            success = False
-            fnow = fold
-
-        # Store relevant variables
-        flog.append(fnow)		# Current function value
-
-        iteration += 1
-        if display:
-            print '\r',
-            print 'Iteration: {0:>5g}  Objective:{1:> 12e}  Scale:{2:> 12e}'.format(iteration, fnow, beta),
-            # print 'Iteration:', iteration, ' Objective:', fnow, '  Scale:', beta, '\r',
-            sys.stdout.flush()
-
-        if success:
-            # Test for termination
-            if (np.max(np.abs(alpha*d)) < xtol) or (np.abs(fnew-fold) < ftol):
-                status='converged'
-                return x, flog, function_eval, status
-
-            else:
-                # Update variables for new position
-                fold = fnew
-                gradold = gradnew
-                gradnew = gradf(x, *optargs)
-                # If the gradient is zero then we are done.
-                if np.dot(gradnew,gradnew) == 0:
-                    return x, flog, function_eval, status
-
-        # Adjust beta according to comparison ratio.
-        if Delta < 0.25:
-            beta = min(4.0*beta, betamax)
-        if Delta > 0.75:
-            beta = max(0.5*beta, betamin)
-
-        # Update search direction using Polak-Ribiere formula, or re-start
-        # in direction of negative gradient after nparams steps.
-        if nsuccess == x.size:
-            d = -gradnew
-            nsuccess = 0
-        elif success:
-            gamma = np.dot(gradold - gradnew,gradnew)/(mu)
-            d = gamma*d - gradnew
-
-    # If we get here, then we haven't terminated in the given number of
-    # iterations.
-    status = "maxiter exceeded"
-
-    return x, flog, function_eval, status

GPy/kern/kern.py

@@ -61,12 +61,12 @@ class kern(parameterised):
 
                 ax.bar(np.arange(len(ard_params)) - 0.4, ard_params)
                 ax.set_xticks(np.arange(len(ard_params)))
-                ax.set_xticklabels([r"${}$".format(i + 1) for i in range(len(ard_params))])
+                ax.set_xticklabels([r"${}$".format(i) for i in range(len(ard_params))])
         return ax
 
     def _transform_gradients(self, g):
         x = self._get_params()
-        [np.put(x,i,x*t.gradfactor(x[i])) for i,t in zip(self.constrained_indices, self.constraints)]
+        [np.put(x, i, x * t.gradfactor(x[i])) for i, t in zip(self.constrained_indices, self.constraints)]
         [np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]]
         if len(self.tied_indices) or len(self.fixed_indices):
             to_remove = np.hstack((self.fixed_indices + [t[1:] for t in self.tied_indices]))
@@ -88,7 +88,7 @@ class kern(parameterised):
         """
         return self.add(other)
 
-    def add(self, other,tensor=False):
+    def add(self, other, tensor=False):
         """
         Add another kernel to this one. Both kernels are defined on the same _space_
         :param other: the other kernel to be added
@@ -103,7 +103,7 @@ class kern(parameterised):
             newkern = kern(D, self.parts + other.parts, self_input_slices + other_input_slices)
 
             # transfer constraints:
-            newkern.constrained_indices = self.constrained_indices + [x+self.Nparam for x in other.constrained_indices]
+            newkern.constrained_indices = self.constrained_indices + [x + self.Nparam for x in other.constrained_indices]
             newkern.constraints = self.constraints + other.constraints
             newkern.fixed_indices = self.fixed_indices + [self.Nparam + x for x in other.fixed_indices]
             newkern.fixed_values = self.fixed_values + other.fixed_values
@@ -113,7 +113,7 @@ class kern(parameterised):
             assert self.D == other.D
             newkern = kern(self.D, self.parts + other.parts, self.input_slices + other.input_slices)
             # transfer constraints:
-            newkern.constrained_indices = self.constrained_indices + [i+self.Nparam  for i in other.constrained_indices]
+            newkern.constrained_indices = self.constrained_indices + [i + self.Nparam  for i in other.constrained_indices]
             newkern.constraints = self.constraints + other.constraints
             newkern.fixed_indices = self.fixed_indices + [self.Nparam + x for x in other.fixed_indices]
             newkern.fixed_values = self.fixed_values + other.fixed_values
@@ -126,7 +126,7 @@ class kern(parameterised):
         """
         return self.prod(other)
 
-    def prod(self, other,tensor=False):
+    def prod(self, other, tensor=False):
         """
         multiply two kernels (either on the same space, or on the tensor product of the input space)
         :param other: the other kernel to be added
@@ -136,12 +136,12 @@ class kern(parameterised):
         K2 = other.copy()
 
         slices = []
-        for sl1, sl2 in itertools.product(K1.input_slices,K2.input_slices):
+        for sl1, sl2 in itertools.product(K1.input_slices, K2.input_slices):
-            s1, s2 = [False]*K1.D, [False]*K2.D
+            s1, s2 = [False] * K1.D, [False] * K2.D
             s1[sl1], s2[sl2] = [True], [True]
-            slices += [s1+s2]
+            slices += [s1 + s2]
 
-        newkernparts = [prod(k1, k2,tensor) for k1, k2 in itertools.product(K1.parts, K2.parts)]
+        newkernparts = [prod(k1, k2, tensor) for k1, k2 in itertools.product(K1.parts, K2.parts)]
 
         if tensor:
             newkern = kern(K1.D + K2.D, newkernparts, slices)
@@ -176,8 +176,8 @@ class kern(parameterised):
         prev_constr_ind = [K1.constrained_indices] + [K1.Nparam + i for i in K2.constrained_indices]
         prev_constr = K1.constraints + K2.constraints
 
-        prev_constr_fix = K1.fixed_indices + [arr + K1.Nparam for arr in K2.fixed_indices]
+        # prev_constr_fix = K1.fixed_indices + [arr + K1.Nparam for arr in K2.fixed_indices]
-        prev_constr_fix_values = K1.fixed_values + K2.fixed_values
+        # prev_constr_fix_values = K1.fixed_values + K2.fixed_values
 
         # follow the previous ties
         for arr in prev_ties:
@@ -189,8 +189,8 @@ class kern(parameterised):
             index = np.where(index_param == i)[0]
             if index.size > 1:
                 self.tie_params(index)
-        for i,t in zip(prev_constr_ind,prev_constr):
+        for i, t in zip(prev_constr_ind, prev_constr):
-            self.constrain(np.where(index_param == i)[0],t)
+            self.constrain(np.where(index_param == i)[0], t)
 
     def _get_params(self):
         return np.hstack([p._get_params() for p in self.parts])
@@ -208,15 +208,15 @@ class kern(parameterised):
         return sum([[name + '_' + n for n in k._get_param_names()] for name, k in zip(names, self.parts)], [])
 
     def K(self, X, X2=None, which_parts='all'):
-        if which_parts=='all':
+        if which_parts == 'all':
-            which_parts = [True]*self.Nparts
+            which_parts = [True] * self.Nparts
         assert X.shape[1] == self.D
         if X2 is None:
             target = np.zeros((X.shape[0], X.shape[0]))
             [p.K(X[:, i_s], None, target=target) for p, i_s, part_i_used in zip(self.parts, self.input_slices, which_parts) if part_i_used]
         else:
             target = np.zeros((X.shape[0], X2.shape[0]))
-            [p.K(X[:, i_s], X2[:,i_s], target=target) for p, i_s, part_i_used in zip(self.parts, self.input_slices, which_parts) if part_i_used]
+            [p.K(X[:, i_s], X2[:, i_s], target=target) for p, i_s, part_i_used in zip(self.parts, self.input_slices, which_parts) if part_i_used]
         return target
 
     def dK_dtheta(self, dL_dK, X, X2=None):
@@ -248,8 +248,8 @@ class kern(parameterised):
         return target
 
     def Kdiag(self, X, which_parts='all'):
-        if which_parts=='all':
+        if which_parts == 'all':
-            which_parts = [True]*self.Nparts
+            which_parts = [True] * self.Nparts
         assert X.shape[1] == self.D
         target = np.zeros(X.shape[0])
         [p.Kdiag(X[:, i_s], target=target) for p, i_s, part_on in zip(self.parts, self.input_slices, which_parts) if part_on]
@@ -270,22 +270,22 @@ class kern(parameterised):
 
     def psi0(self, Z, mu, S):
         target = np.zeros(mu.shape[0])
-        [p.psi0(Z[:,i_s], mu[:,i_s], S[:,i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
+        [p.psi0(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
         return target
 
     def dpsi0_dtheta(self, dL_dpsi0, Z, mu, S):
         target = np.zeros(self.Nparam)
-        [p.dpsi0_dtheta(dL_dpsi0, Z[:,i_s], mu[:,i_s], S[:,i_s], target[ps]) for p, ps, i_s in zip(self.parts, self.param_slices, self.input_slices)]
+        [p.dpsi0_dtheta(dL_dpsi0, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, ps, i_s in zip(self.parts, self.param_slices, self.input_slices)]
         return self._transform_gradients(target)
 
     def dpsi0_dmuS(self, dL_dpsi0, Z, mu, S):
         target_mu, target_S = np.zeros_like(mu), np.zeros_like(S)
-        [p.dpsi0_dmuS(dL_dpsi0, Z[:,i_s], mu[:,i_s], S[:,i_s], target_mu[:,i_s], target_S[:,i_s]) for p, i_s in zip(self.parts, self.input_slices)]
+        [p.dpsi0_dmuS(dL_dpsi0, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
         return target_mu, target_S
 
     def psi1(self, Z, mu, S):
         target = np.zeros((mu.shape[0], Z.shape[0]))
-        [p.psi1(Z[:,i_s], mu[:,i_s], S[:,i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
+        [p.psi1(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
         return target
 
     def dpsi1_dtheta(self, dL_dpsi1, Z, mu, S):
@@ -314,7 +314,7 @@ class kern(parameterised):
         [p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
 
         # compute the "cross" terms
-        #TODO: input_slices needed
+        # TODO: input_slices needed
         for p1, p2 in itertools.combinations(self.parts, 2):
             # white doesn;t combine with anything
             if p1.name == 'white' or p2.name == 'white':
@@ -335,9 +335,9 @@ class kern(parameterised):
                 target += p2.variance * (tmp[:, :, None] + tmp[:, None, :])
             # rbf X linear
             elif p1.name == 'linear' and p2.name == 'rbf':
-                raise NotImplementedError  # TODO
+                raise NotImplementedError # TODO
             elif p2.name == 'linear' and p1.name == 'rbf':
-                raise NotImplementedError  # TODO
+                raise NotImplementedError # TODO
             else:
                 raise NotImplementedError, "psi2 cannot be computed for this kernel"
         return target
@@ -365,7 +365,7 @@ class kern(parameterised):
                 p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1._psi1 * 2., Z, mu, S, target[ps2])
             # linear X bias
             elif p1.name == 'bias' and p2.name == 'linear':
-                p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1.variance * 2., Z, mu, S, target[ps2])  # [ps1])
+                p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1.variance * 2., Z, mu, S, target[ps2]) # [ps1])
                 psi1 = np.zeros((mu.shape[0], Z.shape[0]))
                 p2.psi1(Z, mu, S, psi1)
                 p1.dpsi1_dtheta(dL_dpsi2.sum(1) * psi1 * 2., Z, mu, S, target[ps1])
@@ -376,9 +376,9 @@ class kern(parameterised):
                 p2.dpsi1_dtheta(dL_dpsi2.sum(1) * psi1 * 2., Z, mu, S, target[ps2])
             # rbf X linear
             elif p1.name == 'linear' and p2.name == 'rbf':
-                raise NotImplementedError  # TODO
+                raise NotImplementedError # TODO
             elif p2.name == 'linear' and p1.name == 'rbf':
-                raise NotImplementedError  # TODO
+                raise NotImplementedError # TODO
             else:
                 raise NotImplementedError, "psi2 cannot be computed for this kernel"
 
@@ -389,7 +389,7 @@ class kern(parameterised):
         [p.dpsi2_dZ(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
 
         # compute the "cross" terms
-        #TODO: we need input_slices here.
+        # TODO: we need input_slices here.
         for p1, p2 in itertools.combinations(self.parts, 2):
             # white doesn;t combine with anything
             if p1.name == 'white' or p2.name == 'white':
@@ -406,9 +406,9 @@ class kern(parameterised):
                 p1.dpsi1_dZ(dL_dpsi2.sum(1).T * p2.variance, Z, mu, S, target)
             # rbf X linear
             elif p1.name == 'linear' and p2.name == 'rbf':
-                raise NotImplementedError  # TODO
+                raise NotImplementedError # TODO
             elif p2.name == 'linear' and p1.name == 'rbf':
-                raise NotImplementedError  # TODO
+                raise NotImplementedError # TODO
             else:
                 raise NotImplementedError, "psi2 cannot be computed for this kernel"
 
@@ -419,7 +419,7 @@ class kern(parameterised):
         [p.dpsi2_dmuS(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
 
         # compute the "cross" terms
-        #TODO: we need input_slices here.
+        # TODO: we need input_slices here.
         for p1, p2 in itertools.combinations(self.parts, 2):
             # white doesn;t combine with anything
             if p1.name == 'white' or p2.name == 'white':
@@ -436,9 +436,9 @@ class kern(parameterised):
                 p1.dpsi1_dmuS(dL_dpsi2.sum(1).T * p2.variance * 2., Z, mu, S, target_mu, target_S)
             # rbf X linear
             elif p1.name == 'linear' and p2.name == 'rbf':
-                raise NotImplementedError  # TODO
+                raise NotImplementedError # TODO
             elif p2.name == 'linear' and p1.name == 'rbf':
-                raise NotImplementedError  # TODO
+                raise NotImplementedError # TODO
             else:
                 raise NotImplementedError, "psi2 cannot be computed for this kernel"
 

GPy/likelihoods/Gaussian.py

@@ -14,7 +14,7 @@ class Gaussian(likelihood):
     def __init__(self, data, variance=1., normalize=False):
         self.is_heteroscedastic = False
         self.Nparams = 1
-        self.Z = 0.  # a correction factor which accounts for the approximation made
+        self.Z = 0. # a correction factor which accounts for the approximation made
         N, self.D = data.shape
 
         # normalization
@@ -53,10 +53,10 @@ class Gaussian(likelihood):
     def _set_params(self, x):
         x = float(x)
         if self._variance != x:
-            self._variance = x
+            self.precision = 1. / x
-            self.covariance_matrix = np.eye(self.N) * self._variance
+            self.covariance_matrix = np.eye(self.N) * x
-            self.precision = 1. / self._variance
             self.V = (self.precision) * self.Y
+            self._variance = x
 
     def predictive_values(self, mu, var, full_cov):
         """
@@ -69,6 +69,7 @@ class Gaussian(likelihood):
                 # Note. for D>1, we need to re-normalise all the outputs independently.
                 # This will mess up computations of diag(true_var), below.
                 # note that the upper, lower quantiles should be the same shape as mean
+            # Augment the output variance with the likelihood variance and rescale.
             true_var = (var + np.eye(var.shape[0]) * self._variance) * self._scale ** 2
             _5pc = mean - 2.*np.sqrt(np.diag(true_var))
             _95pc = mean + 2.*np.sqrt(np.diag(true_var))

GPy/models/Bayesian_GPLVM.py

@@ -27,7 +27,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
 
     """
     def __init__(self, Y, Q, X=None, X_variance=None, init='PCA', M=10,
-                 Z=None, kernel=None, oldpsave=5, _debug=False,
+                 Z=None, kernel=None, oldpsave=10, _debug=False,
                  **kwargs):
         if X == None:
             X = self.initialise_latent(init, Q, Y)
@@ -87,19 +87,19 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         return x
 
     def _set_params(self, x, save_old=True, save_count=0):
-        try:
+#         try:
             N, Q = self.N, self.Q
             self.X = x[:self.X.size].reshape(N, Q).copy()
             self.X_variance = x[(N * Q):(2 * N * Q)].reshape(N, Q).copy()
             sparse_GP._set_params(self, x[(2 * N * Q):])
-            self.oldps = x
+#             self.oldps = x
-        except (LinAlgError, FloatingPointError, ZeroDivisionError):
+#         except (LinAlgError, FloatingPointError, ZeroDivisionError):
-            print "\rWARNING: Caught LinAlgError, continueing without setting            "
+#             print "\rWARNING: Caught LinAlgError, continueing without setting            "
-            if self._debug:
+#             if self._debug:
-                self._savederrors.append(self.f_call)
+#                 self._savederrors.append(self.f_call)
-            if save_count > 10:
+#             if save_count > 10:
-                raise
+#                 raise
-            self._set_params(self.oldps[-1], save_old=False, save_count=save_count + 1)
+#             self._set_params(self.oldps[-1], save_old=False, save_count=save_count + 1)
 
     def dKL_dmuS(self):
         dKL_dS = (1. - (1. / (self.X_variance))) * 0.5
@@ -167,8 +167,12 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
 #         d_dmu = (dL_dmu).flatten()
 #         d_dS = (dL_dS).flatten()
         # ========================
-        dbound_dmuS = np.hstack((d_dmu, d_dS))
+        self.dbound_dmuS = np.hstack((d_dmu, d_dS))
-        return np.hstack((dbound_dmuS.flatten(), sparse_GP._log_likelihood_gradients(self)))
+        self.dbound_dZtheta = sparse_GP._log_likelihood_gradients(self)
+        return np.hstack((self.dbound_dmuS.flatten(), self.dbound_dZtheta))
+
+    def _log_likelihood_normal_gradients(self):
+        Si, _, _, _ = pdinv(self.X_variance)
 
     def plot_latent(self, which_indices=None, *args, **kwargs):
 
@@ -204,23 +208,25 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         else:
             colors = iter(colors)
         plots = []
+        x = np.arange(self.X.shape[0])
         for i in range(self.X.shape[1]):
             if axes is None:
                 ax = fig.add_subplot(self.X.shape[1], 1, i + 1)
             else:
                 ax = axes[i]
             ax.plot(self.X, c='k', alpha=.3)
-            plots.extend(ax.plot(self.X.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
+            plots.extend(ax.plot(x, self.X.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
-            ax.fill_between(np.arange(self.X.shape[0]),
+            ax.fill_between(x,
                             self.X.T[i] - 2 * np.sqrt(self.X_variance.T[i]),
                             self.X.T[i] + 2 * np.sqrt(self.X_variance.T[i]),
                             facecolor=plots[-1].get_color(),
                             alpha=.3)
             ax.legend(borderaxespad=0.)
+            ax.set_xlim(x.min(), x.max())
             if i < self.X.shape[1] - 1:
                 ax.set_xticklabels('')
         pylab.draw()
-        fig.tight_layout(h_pad=.01)  # , rect=(0, 0, 1, .95))
+        fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
         return fig
 
     def _debug_filter_params(self, x):
@@ -259,11 +265,11 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         kllls = np.array(self._savedklll)
         LL, = ax1.plot(kllls[:, 0], kllls[:, 1] - kllls[:, 2], '-', label=r'$\log p(\mathbf{Y})$', mew=1.5)
         KL, = ax1.plot(kllls[:, 0], kllls[:, 2], '-', label=r'$\mathcal{KL}(p||q)$', mew=1.5)
-        L, = ax1.plot(kllls[:, 0], kllls[:, 1], '-', label=r'$L$', mew=1.5)  # \mathds{E}_{q(\mathbf{X})}[p(\mathbf{Y|X})\frac{p(\mathbf{X})}{q(\mathbf{X})}]
+        L, = ax1.plot(kllls[:, 0], kllls[:, 1], '-', label=r'$L$', mew=1.5) # \mathds{E}_{q(\mathbf{X})}[p(\mathbf{Y|X})\frac{p(\mathbf{X})}{q(\mathbf{X})}]
 
         param_dict = dict(self._savedparams)
         gradient_dict = dict(self._savedgradients)
-        kmm_dict = dict(self._savedpsiKmm)
+#         kmm_dict = dict(self._savedpsiKmm)
         iters = np.array(param_dict.keys())
         ABCD_dict = np.array(self._savedABCD)
         self.showing = 0
@@ -407,7 +413,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
 
         # parameter changes
         # ax2 = pylab.subplot2grid((4, 1), (1, 0), 3, 1, projection='3d')
-        button_options = [0, 0]  # [0]: clicked -- [1]: dragged
+        button_options = [0, 0] # [0]: clicked -- [1]: dragged
 
         def update_plots(event):
             if button_options[0] and not button_options[1]:
@@ -479,4 +485,4 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         cidp = figs[0].canvas.mpl_connect('button_press_event', onclick)
         cidd = figs[0].canvas.mpl_connect('motion_notify_event', motion)
 
-        return ax1, ax2, ax3, ax4, ax5  # , ax6, ax7
+        return ax1, ax2, ax3, ax4, ax5 # , ax6, ax7

GPy/models/GPLVM.py

@@ -28,7 +28,7 @@ class GPLVM(GP):
         if X is None:
             X = self.initialise_latent(init, Q, Y)
         if kernel is None:
-            kernel = kern.rbf(Q) + kern.bias(Q)
+            kernel = kern.rbf(Q, ARD=Q>1) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
         likelihood = Gaussian(Y, normalize=normalize_Y)
         GP.__init__(self, X, likelihood, kernel, **kwargs)
 

GPy/models/mrd.py

@@ -93,7 +93,7 @@ class MRD(model):
         self.NQ = self.N * self.Q
         self.MQ = self.M * self.Q
 
-        model.__init__(self)  # @UndefinedVariable
+        model.__init__(self) # @UndefinedVariable
 
     @property
     def X(self):
@@ -255,7 +255,7 @@ class MRD(model):
                 X[:, qs] = PCA(Y, len(qs))[0]
         elif init in "PCA_concat":
             X = PCA(numpy.hstack(Ylist), self.Q)[0]
-        else:  # init == 'random':
+        else: # init == 'random':
             X = numpy.random.randn(Ylist[0].shape[0], self.Q)
         self.X = X
         return X

GPy/models/sparse_GP.py

@@ -76,7 +76,7 @@ class sparse_GP(GP):
 #                 psi2_beta_scaled = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1) / sf2)).sum(0)
                 psi2_beta_scaled = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1))).sum(0)
                 evals, evecs = linalg.eigh(psi2_beta_scaled)
-                clipped_evals = np.clip(evals, 0., 1e15)  # TODO: make clipping configurable
+                clipped_evals = np.clip(evals, 0., 1e6) # TODO: make clipping configurable
                 if not np.allclose(evals, clipped_evals):
                     print "Warning: clipping posterior eigenvalues"
                 tmp = evecs * np.sqrt(clipped_evals)

GPy/util/datasets.py

@@ -4,6 +4,7 @@ import numpy as np
 import GPy
 import scipy.sparse
 import scipy.io
+import cPickle as pickle
 data_path = os.path.join(os.path.dirname(__file__), 'datasets')
 default_seed = 10000
 
@@ -96,16 +97,29 @@ def stick():
     lbls = 'connect'
     return {'Y': Y, 'connect' : connect, 'info': "Stick man data from Ohio."}
 
+def swiss_roll_generated(N=1000, sigma=0.0):
+    with open(os.path.join(data_path, 'swiss_roll.pickle')) as f:
+        data = pickle.load(f)
+    Na = data['Y'].shape[0]
+    perm = np.random.permutation(np.r_[:Na])[:N]
+    Y = data['Y'][perm, :]
+    t = data['t'][perm]
+    c = data['colors'][perm, :]
+    so = np.argsort(t)
+    Y = Y[so, :]
+    t = t[so]
+    c = c[so, :]
+    return {'Y':Y, 't':t, 'colors':c}
 
 def swiss_roll_1000():
     mat_data = scipy.io.loadmat(os.path.join(data_path, 'swiss_roll_data'))
     Y = mat_data['X_data'][:, 0:1000].transpose()
     return {'Y': Y, 'info': "Subsample of the swiss roll data extracting only the first 1000 values."}
 
-def swiss_roll():
+def swiss_roll(N=3000):
     mat_data = scipy.io.loadmat(os.path.join(data_path, 'swiss_roll_data.mat'))
-    Y = mat_data['X_data'][:, 0:3000].transpose()
+    Y = mat_data['X_data'][:, 0:N].transpose()
-    return {'Y': Y, 'info': "The first 3,000 points from the swiss roll data of Tennenbaum, de Silva and Langford (2001)."}
+    return {'Y': Y, 'X': mat_data['X_data'], 'info': "The first 3,000 points from the swiss roll data of Tennenbaum, de Silva and Langford (2001)."}
 
 def toy_rbf_1d(seed=default_seed):
     np.random.seed(seed=seed)
@@ -270,13 +284,13 @@ def cmu_mocap(subject, train_motions, test_motions=[], sample_every=4):
 
     end_ind = 0
     for i in range(len(temp_Y)):
-        start_ind = end_ind 
+        start_ind = end_ind
         end_ind += temp_Y[i].shape[0]
         Y[start_ind:end_ind, :] = temp_Y[i]
         lbls[start_ind:end_ind, :] = temp_lbls[i]
-    if len(test_motions)>0:
+    if len(test_motions) > 0:
         temp_Ytest = []
-        temp_lblstest = [] 
+        temp_lblstest = []
 
         testexlbls = np.eye(len(test_motions))
         tot_test_length = 0
@@ -292,7 +306,7 @@ def cmu_mocap(subject, train_motions, test_motions=[], sample_every=4):
 
         end_ind = 0
         for i in range(len(temp_Ytest)):
-            start_ind = end_ind 
+            start_ind = end_ind
             end_ind += temp_Ytest[i].shape[0]
             Ytest[start_ind:end_ind, :] = temp_Ytest[i]
             lblstest[start_ind:end_ind, :] = temp_lblstest[i]
@@ -304,7 +318,7 @@ def cmu_mocap(subject, train_motions, test_motions=[], sample_every=4):
     for motion in train_motions:
         info += motion + ', '
     info = info[:-2]
-    if len(test_motions)>0:
+    if len(test_motions) > 0:
         info += '. Test motions: '
         for motion in test_motions:
             info += motion + ', '

GPy/util/visualize.py

@@ -44,7 +44,7 @@ class vector_show(data_show):
 
 
 class lvm(data_show):
-    def __init__(self, vals, model, data_visualize, latent_axes=None, latent_index=[0,1]):
+    def __init__(self, vals, model, data_visualize, latent_axes=None, sense_axes=None, latent_index=[0,1]):
         """Visualize a latent variable model
 
         :param model: the latent variable model to visualize.
@@ -71,7 +71,7 @@ class lvm(data_show):
         self.data_visualize = data_visualize
         self.model = model
         self.latent_axes = latent_axes
-
+        self.sense_axes = sense_axes
         self.called = False
         self.move_on = False
         self.latent_index = latent_index
@@ -81,10 +81,12 @@ class lvm(data_show):
         self.latent_values = vals
         self.latent_handle = self.latent_axes.plot([0],[0],'rx',mew=2)[0]
         self.modify(vals)
+        self.show_sensitivities()
 
     def modify(self, vals):
         """When latent values are modified update the latent representation and ulso update the output visualization."""
         y = self.model.predict(vals)[0]
+        print y
         self.data_visualize.modify(y)
         self.latent_handle.set_data(vals[self.latent_index[0]], vals[self.latent_index[1]])
         self.axes.figure.canvas.draw()
@@ -99,6 +101,7 @@ class lvm(data_show):
         if event.inaxes!=self.latent_axes: return
         self.move_on = not self.move_on
         self.called = True
+
     def on_move(self, event):
         if event.inaxes!=self.latent_axes: return
         if self.called and self.move_on:
@@ -107,22 +110,54 @@ class lvm(data_show):
             self.latent_values[self.latent_index[1]]=event.ydata
             self.modify(self.latent_values)
 
+    def show_sensitivities(self):
+        # A click in the bar chart axis for selection a dimension.
+        if self.sense_axes != None:
+            self.sense_axes.cla()
+            self.sense_axes.bar(np.arange(self.model.Q),1./self.model.input_sensitivity(),color='b')
+
+            if self.latent_index[1] == self.latent_index[0]:
+                self.sense_axes.bar(np.array(self.latent_index[0]),1./self.model.input_sensitivity()[self.latent_index[0]],color='y')
+                self.sense_axes.bar(np.array(self.latent_index[1]),1./self.model.input_sensitivity()[self.latent_index[1]],color='y')
+
+            else:
+                self.sense_axes.bar(np.array(self.latent_index[0]),1./self.model.input_sensitivity()[self.latent_index[0]],color='g')
+                self.sense_axes.bar(np.array(self.latent_index[1]),1./self.model.input_sensitivity()[self.latent_index[1]],color='r')
+
+            self.sense_axes.figure.canvas.draw()
+
+
 class lvm_subplots(lvm):
     """
-    latent_axes is a np array of dimension np.ceil(Q/2) + 1,
+    latent_axes is a np array of dimension np.ceil(Q/2),
-    one for each pair of the axes, and the last one for the sensitiity bar chart
+    one for each pair of the latent dimensions.
     """
-    def __init__(self, vals, model, data_visualize, latent_axes=None, latent_index=[0,1]):
+    def __init__(self, vals, model, data_visualize, latent_axes=None, sense_axes=None):
-        lvm.__init__(self, vals, model,data_visualize,latent_axes,[0,1])
         self.nplots = int(np.ceil(model.Q/2.))+1
-        lvm.__init__(self,model,data_visualize,latent_axes,latent_index)
+        assert len(latent_axes)==self.nplots
-        self.latent_values = np.zeros(2*np.ceil(self.model.Q/2.)) # possibly an extra dimension on this
+        if vals==None:
-        assert latent_axes.size == self.nplots
+            vals = model.X[0, :]
+        self.latent_values = vals 
+
+        for i, axis in enumerate(latent_axes):
+            if i == self.nplots-1:
+                if self.nplots*2!=model.Q:
+                    latent_index = [i*2, i*2]
+                lvm.__init__(self, self.latent_vals, model, data_visualize, axis, sense_axes, latent_index=latent_index)
+            else:
+                latent_index = [i*2, i*2+1]
+                lvm.__init__(self, self.latent_vals, model, data_visualize, axis, latent_index=latent_index)
+
 
 
 class lvm_dimselect(lvm):
     """
     A visualizer for latent variable models which allows selection of the latent dimensions to use by clicking on a bar chart of their length scales.
+
+    For an example of the visualizer's use try:
+    
+    GPy.examples.dimensionality_reduction.BGPVLM_oil()
+
     """
     def __init__(self, vals, model, data_visualize, latent_axes=None, sense_axes=None, latent_index=[0, 1]):
         if latent_axes==None and sense_axes==None:
@@ -133,24 +168,9 @@ class lvm_dimselect(lvm):
         else:
             self.sense_axes = sense_axes
         
-        lvm.__init__(self,vals,model,data_visualize,latent_axes,latent_index)
+        lvm.__init__(self,vals,model,data_visualize,latent_axes,sense_axes,latent_index)
-        self.show_sensitivities()
         print "use left and right mouse butons to select dimensions"
 
-    def show_sensitivities(self):
-        # A click in the bar chart axis for selection a dimension.
-        self.sense_axes.cla()
-        self.sense_axes.bar(np.arange(self.model.Q),1./self.model.input_sensitivity(),color='b')
-
-        if self.latent_index[1] == self.latent_index[0]:
-            self.sense_axes.bar(np.array(self.latent_index[0]),1./self.model.input_sensitivity()[self.latent_index[0]],color='y')
-            self.sense_axes.bar(np.array(self.latent_index[1]),1./self.model.input_sensitivity()[self.latent_index[1]],color='y')
-            
-        else:
-            self.sense_axes.bar(np.array(self.latent_index[0]),1./self.model.input_sensitivity()[self.latent_index[0]],color='g')
-            self.sense_axes.bar(np.array(self.latent_index[1]),1./self.model.input_sensitivity()[self.latent_index[1]],color='r')
-
-        self.sense_axes.figure.canvas.draw()
 
     def on_click(self, event):
 
@@ -177,12 +197,6 @@ class lvm_dimselect(lvm):
         self.called = True
 
 
-    def on_move(self, event):
-        if event.inaxes!=self.latent_axes: return
-        if self.called and self.move_on:
-            self.latent_values[self.latent_index[0]]=event.xdata
-            self.latent_values[self.latent_index[1]]=event.ydata            
-            self.modify(self.latent_values)
 
     def on_leave(self,event):
         latent_values = self.latent_values.copy()
@@ -193,7 +207,7 @@ class lvm_dimselect(lvm):
 
 class image_show(data_show):
     """Show a data vector as an image."""
-    def __init__(self, vals, axes=None, dimensions=(16,16), transpose=False, invert=False, scale=False, palette=[], presetMean = 0., presetSTD = -1., selectImage = 0):
+    def __init__(self, vals, axes=None, dimensions=(16,16), transpose=False, invert=False, scale=False, palette=[], presetMean = 0., presetSTD = -1., selectImage=0):
         data_show.__init__(self, vals, axes)
         self.dimensions = dimensions
         self.transpose = transpose
@@ -214,15 +228,30 @@ class image_show(data_show):
     def modify(self, vals):
         self.set_image(vals)
         self.handle.set_array(self.vals)
-        self.axes.figure.canvas.draw() # Teo - original line: plt.show()
+        self.axes.figure.canvas.draw() 
 
     def set_image(self, vals):
         dim = self.dimensions[0] * self.dimensions[1]
-        self.vals = np.reshape(vals[0,dim*self.selectImage+np.array(range(dim))], self.dimensions, order='F')
+        nImg = np.sqrt(vals[0,].size/dim)
+        if nImg > 1 and nImg.is_integer(): # Show a mosaic of images
+            nImg = np.int(nImg)
+            self.vals = np.zeros((self.dimensions[0]*nImg, self.dimensions[1]*nImg))
+            for iR in range(nImg):
+                for iC in range(nImg):
+                    currImgId = iR*nImg + iC
+                    currImg = np.reshape(vals[0,dim*currImgId+np.array(range(dim))], self.dimensions, order='F')
+                    firstRow = iR*self.dimensions[0]
+                    lastRow = (iR+1)*self.dimensions[0]
+                    firstCol = iC*self.dimensions[1]
+                    lastCol = (iC+1)*self.dimensions[1]
+                    self.vals[firstRow:lastRow, firstCol:lastCol] = currImg
+
+        else: 
+            self.vals = np.reshape(vals[0,dim*self.selectImage+np.array(range(dim))], self.dimensions, order='F')
         if self.transpose:
             self.vals = self.vals.T.copy()
-        if not self.scale:
+        # if not self.scale:
-            self.vals = self.vals
+        #     self.vals = self.vals
         if self.invert:
             self.vals = -self.vals
 
@@ -318,7 +347,7 @@ class stick_show(mocap_data_show):
 
     def process_values(self, vals):
         self.vals = vals.reshape((3, vals.shape[1]/3)).T.copy()
-    
+
 class skeleton_show(mocap_data_show):
     """data_show class for visualizing motion capture data encoded as a skeleton with angles."""
     def __init__(self, vals, skel, padding=0, axes=None):