Merge branch 'mrd' 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/transformations.py

@@ -39,23 +39,26 @@ class logexp(transformation):
         return '(+ve)'
 
 class logexp_clipped(transformation):
-    def __init__(self):
+    def __init__(self, lower=1e-8, upper=1e200):
         self.domain = 'positive'
+        self.lower = lower
+        self.upper = upper
     def f(self, x):
-        f = np.log(1. + np.exp(x))
+        exp = np.exp(x)
+        f = np.log(1. + np.where(exp > self.upper, self.upper, exp))
         return f
     def finv(self, f):
         return np.log(np.exp(f) - 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):
         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):

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,88 @@ 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(0)
+    Y /= Y.std(0)
+
+    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, 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'] = 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):
     data = GPy.util.datasets.oil()
+    from GPy.core.transformations import logexp_clipped
+    np.random.seed(0)
 
     # 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', logexp_clipped())
-    m.constrain('length', logexp_clipped())
+#     m.constrain('length', logexp_clipped())
-    m['lengt'] = 100.
+    m['lengt'] = m.X.var(0).max() / m.X.var(0)
+    m['noise'] = Yn.var() / 100.
+
     m.ensure_default_constraints()
 
     # optimize
     if optimize:
-        m.unconstrain('noise'); m.constrain_fixed('noise', Y.var() / 100.)
+#         m.unconstrain('noise'); m.constrain_fixed('noise')
-        m.optimize('scg', messages=1, max_f_eval=150)
+#         m.optimize('scg', messages=1, max_f_eval=200)
-
+#         m.unconstrain('noise')
-        m.unconstrain('noise')
+#         m.constrain('noise', logexp_clipped())
-        m.constrain('noise', logexp_clipped())
         m.optimize('scg', messages=1, max_f_eval=max_f_eval)
 
     if plot:
@@ -115,6 +175,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))
@@ -178,6 +240,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']
@@ -213,6 +276,8 @@ def bgplvm_simulation(burnin='scg', plot_sim=False,
                       max_burnin=100, true_X=False,
                       do_opt=True,
                       max_f_eval=1000):
+    from GPy.core.transformations import logexp_clipped
+
     D1, D2, D3, N, M, Q = 15, 8, 8, 350, 3, 6
     slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
 
@@ -317,6 +382,8 @@ def mrd_simulation(plot_sim=False):
 
     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)

GPy/inference/SCG.py

@@ -111,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, 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()
 

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,7 +61,7 @@ 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):
@@ -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:
@@ -196,6 +196,7 @@ class kern(parameterised):
         return np.hstack([p._get_params() for p in self.parts])
 
     def _set_params(self, x):
+        x = np.clip(x, -1e300, 1e300)
         [p._set_params(x[s]) for p, s in zip(self.parts, self.param_slices)]
 
     def _get_param_names(self):

GPy/likelihoods/Gaussian.py

@@ -54,8 +54,8 @@ class Gaussian(likelihood):
         x = float(x)
         if self._variance != x:
             self._variance = x
-            self.covariance_matrix = np.eye(self.N) * self._variance
             self.precision = 1. / self._variance
+            self.covariance_matrix = np.eye(self.N) * self._variance
             self.V = (self.precision) * self.Y
 
     def predictive_values(self, mu, var, full_cov):

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):
 
@@ -263,7 +267,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
 
         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

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)
@@ -274,7 +288,7 @@ def cmu_mocap(subject, train_motions, test_motions=[], sample_every=4):
         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 = []
 
@@ -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 + ', '