new functions mrd init_X update

GPy/examples/dimensionality_reduction.py

@@ -6,7 +6,6 @@ import pylab as pb
 from matplotlib import pyplot as plt, pyplot
 
 import GPy
-from GPy.models.mrd import MRD
 
 default_seed = np.random.seed(123344)
 
@@ -102,10 +101,10 @@ def oil_100():
 
 def mrd_simulation(plot_sim=False):
     # num = 2
-    ard1 = np.array([1., 1, 0, 0], dtype=float)
+#     ard1 = np.array([1., 1, 0, 0], dtype=float)
-    ard2 = np.array([0., 1, 1, 0], dtype=float)
+#     ard2 = np.array([0., 1, 1, 0], dtype=float)
-    ard1[ard1 == 0] = 1E-10
+#     ard1[ard1 == 0] = 1E-10
-    ard2[ard2 == 0] = 1E-10
+#     ard2[ard2 == 0] = 1E-10
 
 #     ard1i = 1. / ard1
 #     ard2i = 1. / ard2
@@ -119,46 +118,74 @@ def mrd_simulation(plot_sim=False):
 #     Y2 -= Y2.mean(0)
 #     make_params = lambda ard: np.hstack([[1], ard, [1, .3]])
 
-    D1, D2, N, M, Q = 5, 10, 150, 15, 3
+    D1, D2, D3, N, M, Q = 5, 5, 5, 150, 18, 5
-    x = np.linspace(0, 4 * np.pi, N)[:, None]
+    x = np.linspace(0, 2 * np.pi, N)[:, None]
 
     s1 = np.vectorize(lambda x: np.sin(x))
     s2 = np.vectorize(lambda x: np.cos(x))
+    s3 = np.vectorize(lambda x: np.cos(4 * x))
     sS = np.vectorize(lambda x: np.sin(2 * x))
 
-    S1 = np.hstack([s1(x), sS(x)]) + .1 * np.random.randn(N, 2)
+    s1 = s1(x)
-    S2 = np.hstack([s2(x), sS(x)]) + .1 * np.random.randn(N, 2)
+    s2 = s2(x)
+    s3 = s3(x)
+    sS = sS(x)
+
+    s1 /= np.abs(s1).max()
+    s2 /= np.abs(s2).max()
+    s3 /= np.abs(s3).max()
+    sS /= np.abs(sS).max()
+
+    S1 = np.hstack([s1, sS])
+    S2 = np.hstack([s2, sS])
+    S3 = np.hstack([s3, sS])
 
     Y1 = S1.dot(np.random.randn(S1.shape[1], D1))
     Y2 = S2.dot(np.random.randn(S2.shape[1], D2))
+    Y3 = S3.dot(np.random.randn(S3.shape[1], D3))
+
+    Y1 += .041 * np.random.randn(*Y1.shape)
+    Y2 += .041 * np.random.randn(*Y2.shape)
+    Y3 += .041 * np.random.randn(*Y3.shape)
 
     Y1 -= Y1.mean(0)
     Y2 -= Y2.mean(0)
+    Y3 -= Y3.mean(0)
+    Y1 /= Y1.std(0)
+    Y2 /= Y2.std(0)
+    Y3 /= Y3.std(0)
+
+    Slist = [s1, s2, s3, sS]
+    Ylist = [Y1, Y2, Y3]
 
     if plot_sim:
         import pylab
-        fig = pylab.figure("MRD Simulation")
+        import itertools
-        ax = fig.add_subplot(2, 2, 1)
+        fig = pylab.figure("MRD Simulation", figsize=(12, 12))
-        ax.imshow(S1)
+        fig.clf()
-        ax.set_title("S1")
+        ax = fig.add_subplot(2, 1, 1)
-        ax = fig.add_subplot(2, 2, 2)
+        labls = sorted(filter(lambda x: x.startswith("s"), locals()))
-        ax.imshow(S2)
+        for S, lab in itertools.izip(Slist, labls):
-        ax.set_title("S2")
+            ax.plot(x, S, label=lab)
-        ax = fig.add_subplot(2, 2, 3)
+        ax.legend()
-        ax.imshow(Y1)
+        for i, Y in enumerate(Ylist):
-        ax.set_title("Y1")
+            ax = fig.add_subplot(2, len(Ylist), len(Slist) + i)
-        ax = fig.add_subplot(2, 2, 4)
+            ax.imshow(Y)
-        ax.imshow(Y2)
+            ax.set_title("Y{}".format(i + 1))
-        ax.set_title("Y2")
         pylab.draw()
         pylab.tight_layout()
 
-    k = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q) + GPy.kern.white(Q)
+    from GPy.models import mrd
-
+    from GPy import kern
-    m = MRD(Y1, Y2, Q=Q, M=M, kernel=k, init="concat", _debug=False)
+    reload(mrd); reload(kern)
+    k = kern.rbf(Q, ARD=True) + kern.bias(Q) + kern.white(Q)
+    m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, init="concat", _debug=False)
     m.ensure_default_constraints()
-    cstr = "noise|white"
+
-    m.unconstrain(cstr); m.constrain_bounded(cstr, 1e-3, 1.)
+    # cstr = "noise|white|variance"
+    # m.unconstrain(cstr); m.constrain_bounded(cstr, 1e-6, 1.)
+
+    m.auto_scale_factor = True
 
 #     fig = pyplot.figure("expected", figsize=(8, 3))
 #     ax = fig.add_subplot(121)

GPy/inference/SCG.py

@@ -105,7 +105,7 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
         iteration += 1
         if display:
             print '\r',
-            print 'Iteration: {0:<5g}  Objective:{1:< 12g}  Scale:{2:< 12g}'.format(iteration, fnow, beta),
+            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()
 

GPy/models/__init__.py

@@ -11,7 +11,5 @@ from warped_GP import warpedGP
 from sparse_GPLVM import sparse_GPLVM
 from uncollapsed_sparse_GP import uncollapsed_sparse_GP
 from Bayesian_GPLVM import Bayesian_GPLVM
-import mrd
+from mrd import MRD
-MRD = mrd.MRD
-del mrd
 from generalized_FITC import generalized_FITC

GPy/models/mrd.py

@@ -8,7 +8,6 @@ from GPy.models.Bayesian_GPLVM import Bayesian_GPLVM
 import numpy
 from GPy.models.sparse_GP import sparse_GP
 import itertools
-from matplotlib import pyplot
 import pylab
 from GPy.util.linalg import PCA
 
@@ -59,6 +58,8 @@ class MRD(model):
         if kwargs.has_key('init'):
             init = kwargs['init']
             del kwargs['init']
+        else:
+            init = "PCA"
         try:
             self.Q = kwargs["Q"]
         except KeyError:
@@ -68,11 +69,11 @@ class MRD(model):
         except KeyError:
             self.M = 10
 
-
+        self._init = True
-        X = self._init_X(Ylist, init)
+        X = self._init_X(init, Ylist)
         Z = numpy.random.permutation(X.copy())[:self.M]
-
         self.bgplvms = [Bayesian_GPLVM(Y, kernel=k(), X=X, Z=Z, **kwargs) for Y in Ylist]
+        del self._init
 
         self.gref = self.bgplvms[0]
         nparams = numpy.array([0] + [sparse_GP._get_params(g).size - g.Z.size for g in self.bgplvms])
@@ -84,6 +85,41 @@ class MRD(model):
 
         model.__init__(self)  # @UndefinedVariable
 
+    @property
+    def X(self):
+        return self.gref.X
+    @X.setter
+    def X(self, X):
+        try:
+            self.propagate_param(X=X)
+        except AttributeError:
+            if not self._init:
+                raise AttributeError("bgplvm list not initialized")
+    @property
+    def Ylist(self):
+        return [g.likelihood.Y for g in self.bgplvms]
+    @Ylist.setter
+    def Ylist(self, Ylist):
+        for g, Y in itertools.izip(self.bgplvms, Ylist):
+            g.likelihood.Y = Y
+
+    @property
+    def auto_scale_factor(self):
+        """
+        set auto_scale_factor for all gplvms
+        :param b: auto_scale_factor
+        :type b:
+        """
+        return self.gref.auto_scale_factor
+    @auto_scale_factor.setter
+    def auto_scale_factor(self, b):
+        self.propagate_param(auto_scale_factor=b)
+
+    def propagate_param(self, **kwargs):
+        for key, val in kwargs.iteritems():
+            for g in self.bgplvms:
+                g.__setattr__(key, val)
+
     def _get_param_names(self):
         # X_names = sum([['X_%i_%i' % (n, q) for q in range(self.Q)] for n in range(self.N)], [])
         # S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.Q)] for n in range(self.N)], [])
@@ -129,11 +165,11 @@ class MRD(model):
 
     def _set_params(self, x):
         start = 0; end = self.NQ
-        X = x[start:end].reshape(self.N, self.Q)
+        X = x[start:end]
         start = end; end += start
-        X_var = x[start:end].reshape(self.N, self.Q)
+        X_var = x[start:end]
         start = end; end += self.MQ
-        Z = x[start:end].reshape(self.M, self.Q)
+        Z = x[start:end]
         thetas = x[end:]
 
         if self._debug:
@@ -144,10 +180,14 @@ class MRD(model):
 
         # set params for all:
         for g, s, e in itertools.izip(self.bgplvms, self.nparams, self.nparams[1:]):
-            self._set_var_params(g, X, X_var, Z)
+            g._set_params(numpy.hstack([X, X_var, Z, thetas[s:e]]))
-            self._set_kern_params(g, thetas[s:e].copy())
+#             self._set_var_params(g, X, X_var, Z)
-            g._compute_kernel_matrices()
+#             self._set_kern_params(g, thetas[s:e].copy())
-            g._computations()
+#             g._compute_kernel_matrices()
+#             if self.auto_scale_factor:
+#                 g.scale_factor = numpy.sqrt(g.psi2.sum(0).mean() * g.likelihood.precision)
+# #                 self.scale_factor = numpy.sqrt(self.psi2.sum(0).mean() * self.likelihood.precision)
+#             g._computations()
 
 
     def log_likelihood(self):
@@ -171,15 +211,18 @@ class MRD(model):
                                                 partial=g.partial_for_likelihood)]) \
                               for g in self.bgplvms])))
 
-    def _init_X(self, Ylist, init='PCA_concat'):
+    def _init_X(self, init='PCA', Ylist=None):
-        if init in "PCA_concat":
+        if Ylist is None:
-            X = PCA(numpy.hstack(Ylist), self.Q)[0]
+            Ylist = self.Ylist
-        elif init in "PCA_single":
+        if init in "PCA_single":
             X = numpy.zeros((Ylist[0].shape[0], self.Q))
             for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.Q), len(Ylist)), Ylist):
                 X[:, qs] = PCA(Y, len(qs))[0]
+        elif init in "PCA_concat":
+            X = PCA(numpy.hstack(Ylist), self.Q)[0]
         else:  # init == 'random':
             X = numpy.random.randn(Ylist[0].shape[0], self.Q)
+        self.X = X
         return X
 
     def plot_X(self):
@@ -229,7 +272,7 @@ class MRD(model):
 
     def _debug_optimize(self, opt='scg', maxiters=500, itersteps=10):
         iters = 0
-        optstep = lambda: self.optimize(opt, messages=1, max_iters=itersteps)
+        optstep = lambda: self.optimize(opt, messages=1, max_f_eval=itersteps)
         self._debug_plot()
         raw_input("enter to start debug")
         while iters < maxiters: