added swissroll example

GPy/examples/dimensionality_reduction.py

@@ -2,13 +2,10 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
-import pylab as pb
-from matplotlib import pyplot as plt, pyplot
+from matplotlib import pyplot as plt
 
 import GPy
 from GPy.models.Bayesian_GPLVM import Bayesian_GPLVM
-from GPy.util.datasets import simulation_BGPLVM
-from GPy.core.transformations import square, logexp_clipped
 
 default_seed = np.random.seed(123344)
 
@@ -47,10 +44,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,22 +61,66 @@ 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):
+    from GPy.util.datasets import swiss_roll
+    from GPy.core.transformations import logexp_clipped
+
+    data = swiss_roll(N=N)
+    Y = data['Y']
+    Y -= Y.mean(0)
+    Y /= Y.std(0)
+
+    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)
+
+    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.constrain('variance|length', logexp_clipped())
+    m['lengthscale'] = X.var(0) / X.var(0).max()
+    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=10, M=10, max_f_eval=1e3, plot=False, **k):
     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)
+    Y -= Y.mean(0)
+
+    m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=kernel, M=M, **k)
     m.data_labels = data['Y'][:N].argmax(axis=1)
 
     m.constrain('variance', logexp_clipped())
     m.constrain('length', logexp_clipped())
-    m['lengt'] = 100.
+    m['lengt'] = 1.
+    m['noise'] = Y.var() / 100.
+
     m.ensure_default_constraints()
 
     # optimize
     if optimize:
+        m.unconstrain('X'); m.constrain_fixed('X')
+        m.optimize('scg', messages=1, max_f_eval=10)
+        m.unconstrain('X'); m.constrain('X_var', logexp_clipped())
+
         m.unconstrain('noise'); m.constrain_fixed('noise', Y.var() / 100.)
         m.optimize('scg', messages=1, max_f_eval=150)
 
@@ -115,6 +157,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 +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']
@@ -213,6 +258,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 +364,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)