BGPLVM still failing, doesn't seem to be numerical : (

GPy/examples/dimensionality_reduction.py

@@ -103,9 +103,9 @@ def oil_100():
 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))
-    s2 = np.vectorize(lambda x: x * np.cos(x))
+    s2 = np.vectorize(lambda x: np.cos(x))
-    s3 = np.vectorize(lambda x: np.sin(2 * x))
+    s3 = np.vectorize(lambda x:-np.exp(-np.cos(2 * x)))
-    sS = np.vectorize(lambda x:-np.exp(-np.cos(2 * x)))
+    sS = np.vectorize(lambda x: np.sin(2 * x))
 
     s1 = s1(x)
     s2 = s2(x)
@@ -162,27 +162,57 @@ def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False):
 
     return slist, [S1, S2, S3], Ylist
 
-def bgplvm_simulation(plot_sim=False):
+def bgplvm_simulation(burnin='scg', plot_sim=False, max_f_eval=12):
-    D1, D2, D3, N, M, Q = 50, 34, 8, 100, 2, 6
+    D1, D2, D3, N, M, Q = 2000, 8, 8, 500, 2, 6
     slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
 
     from GPy.models import mrd
     from GPy import kern
     reload(mrd); reload(kern)
 
-    Y = Ylist[0]
+    Y = Ylist[1]
 
-    k = kern.linear(Q, ARD=True) + kern.bias(Q, .01) + kern.white(Q, .1)
+    k = kern.linear(Q, ARD=True) + kern.bias(Q, .0001) + kern.white(Q, .1)
     m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k)
-    m.ensure_default_constraints()
     m.set('noise', Y.var() / 100.)
-    m.auto_scale_factor = True
+#     m.auto_scale_factor = True
+#     m.scale_factor = 1.
 
-    cstr = 'variance'
+    m.ensure_default_constraints()
-    m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-20, 1.)
+
+    if burnin:
+        print "initializing beta"
+        cstr = "noise"
+        m.unconstrain(cstr); m.constrain_fixed(cstr)
+        m.optimize(burnin, messages=1, max_f_eval=max_f_eval)
+
+        print "releasing beta"
+        cstr = "noise"
+        m.unconstrain(cstr);  m.constrain_positive(cstr)
+
+
+# #     cstr = 'variance'
+# #     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-10, 1.)
+#     cstr = 'X_\d'
+#     m.unconstrain(cstr), m.constrain_bounded(cstr, -100., 100.)
+#
+#     cstr = 'noise'
+#     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-3, 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.)  # m.constrain_positive(cstr)
+#
+#     cstr = 'X_variance'
+#     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-10, 1.)  # m.constrain_positive(cstr)
+
+#     np.seterr(all='call')
+#     def ipdbonerr(errtype, flags):
+#         import ipdb; ipdb.set_trace()
+#     np.seterrcall(ipdbonerr)
 
-    cstr = 'linear_variance'
-    m.unconstrain(cstr), m.constrain_positive(cstr)
 
     return m
 
@@ -204,7 +234,7 @@ def mrd_simulation(plot_sim=False):
 #     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 = 50, 34, 8, 100, 2, 6
+    D1, D2, D3, N, M, Q = 2000, 34, 8, 500, 3, 6
     slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
 
     from GPy.models import mrd
@@ -216,24 +246,23 @@ def mrd_simulation(plot_sim=False):
 #     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]]
+    Ylist = Ylist[0:2]
 
     # k = kern.rbf(Q, ARD=True) + kern.bias(Q) + kern.white(Q)
 
-    k = kern.linear(Q, ARD=True) + kern.bias(Q, .01) + kern.white(Q, .1)
+    k = kern.linear(Q, ARD=True) + kern.bias(Q, .01) + kern.white(Q, .001)
     m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, initx="concat", initz='permute', _debug=False)
-    m.ensure_default_constraints()
-    ardvar = 5. / (m.X.max(axis=0) - m.X.min(axis=0))
 
     for i, Y in enumerate(Ylist):
         m.set('{}_noise'.format(i + 1), Y.var() / 100.)
 
+    m.ensure_default_constraints()
 
-    cstr = 'variance'
+#     cstr = 'variance'
-    m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-12, 1.)
+#     m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-12, 1.)
-
+#
-    cstr = 'linear_variance'
+#     cstr = 'linear_variance'
-    m.unconstrain(cstr), m.constrain_positive(cstr)
+#     m.unconstrain(cstr), m.constrain_positive(cstr)
 
 #     print "initializing beta"
 #     cstr = "noise"

GPy/models/Bayesian_GPLVM.py

@@ -9,6 +9,7 @@ from sparse_GP import sparse_GP
 from GPy.util.linalg import pdinv
 from ..likelihoods import Gaussian
 from .. import kern
+from numpy.linalg.linalg import LinAlgError
 
 class Bayesian_GPLVM(sparse_GP, GPLVM):
     """
@@ -22,7 +23,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
     :type init: 'PCA'|'random'
 
     """
-    def __init__(self, Y, Q, X=None, X_variance=None, init='PCA', M=10, Z=None, kernel=None, **kwargs):
+    def __init__(self, Y, Q, X=None, X_variance=None, init='PCA', M=10, Z=None, kernel=None, oldpsave=5, **kwargs):
         if X == None:
             X = self.initialise_latent(init, Q, Y)
 
@@ -36,9 +37,21 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         if kernel is None:
             kernel = kern.rbf(Q) + kern.white(Q)
 
+        self.oldpsave = oldpsave
+        self._oldps = []
 
         sparse_GP.__init__(self, X, Gaussian(Y), kernel, Z=Z, X_variance=X_variance, **kwargs)
 
+    @property
+    def oldps(self):
+        return self._oldps
+    @oldps.setter
+    def oldps(self, p):
+        if len(self._oldps) == (self.oldpsave + 1):
+            self._oldps.pop()
+        # if len(self._oldps) == 0 or not np.any([np.any(np.abs(p - op) > 1e-5) for op in self._oldps]):
+        self._oldps.insert(0, p.copy())
+
     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)], [])
@@ -54,14 +67,19 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         ===============================================================
 
         """
-        return np.hstack((self.X.flatten(), self.X_variance.flatten(), sparse_GP._get_params(self)))
+        x = np.hstack((self.X.flatten(), self.X_variance.flatten(), sparse_GP._get_params(self)))
-
+        return x
-    def _set_params(self, x):
-        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):])
 
+    def _set_params(self, x, save_old=True):
+        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
+        except (LinAlgError, FloatingPointError):
+            print "\rWARNING: Caught LinAlgError, reconstructing old state            "
+            self._set_params(self.oldps[-1], save_old=False)
 
     def dKL_dmuS(self):
         dKL_dS = (1. - (1. / self.X_variance)) * 0.5

GPy/models/mrd.py

@@ -271,14 +271,31 @@ class MRD(model):
         self.Z = Z
         return Z
 
-    def plot_X_1d(self, colors=None):
+    def _handle_plotting(self, fig_num, axes, plotf):
-        fig = pylab.figure(num="MRD X 1d", figsize=(min(8, (3 * len(self.bgplvms))), min(12, (2 * self.X.shape[1]))))
+        if axes is None:
+            fig = pylab.figure(num=fig_num, figsize=(4 * len(self.bgplvms), 3 * len(self.bgplvms)))
+        for i, g in enumerate(self.bgplvms):
+            if axes is None:
+                ax = fig.add_subplot(1, len(self.bgplvms), i + 1)
+            else:
+                ax = axes[i]
+            plotf(i, g, ax)
+        pylab.draw()
+        if axes is None:
+            fig.tight_layout()
+            return fig
+        else:
+            return pylab.gcf()
+
+    def plot_X_1d(self, fig_num="MRD X 1d", axes=None, colors=None):
+        fig = pylab.figure(num=fig_num, figsize=(min(8, (3 * len(self.bgplvms))), min(12, (2 * self.X.shape[1]))))
         if colors is None:
             colors = pylab.gca()._get_lines.color_cycle
             pylab.clf()
         plots = []
         for i in range(self.X.shape[1]):
-            ax = fig.add_subplot(self.X.shape[1], 1, i + 1)
+            if axes is None:
+                ax = fig.add_subplot(self.X.shape[1], 1, i + 1)
             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)))
             ax.fill_between(numpy.arange(self.X.shape[0]),
@@ -289,72 +306,41 @@ class MRD(model):
             ax.legend(borderaxespad=0.)
             if i < self.X.shape[1] - 1:
                 ax.set_xticklabels('')
-#         ax1.legend(plots, [r"$\mathbf{{X_{}}}$".format(i + 1) for i in range(self.X.shape[1])],
-#                    bbox_to_anchor=(0., 1 + .01 * self.X.shape[1],
-#                                    1., 1. + .01 * self.X.shape[1]), loc=3,
-#                    ncol=self.X.shape[1], mode="expand", borderaxespad=0.)
         pylab.draw()
         fig.tight_layout(h_pad=.01)  # , rect=(0, 0, 1, .95))
         return fig
 
-    def plot_X(self):
+    def plot_X(self, fig_num="MRD Predictions", axes=None):
-        fig = pylab.figure("MRD X", figsize=(4 * len(self.bgplvms), 3))
+        fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: ax.imshow(g.X))
-        fig.clf()
-        for i, g in enumerate(self.bgplvms):
-            ax = fig.add_subplot(1, len(self.bgplvms), i + 1)
-            ax.imshow(g.X)
-        pylab.draw()
-        fig.tight_layout()
         return fig
 
-    def plot_predict(self):
+    def plot_predict(self, fig_num="MRD Predictions", axes=None):
-        fig = pylab.figure("MRD Predictions", figsize=(4 * len(self.bgplvms), 3))
+        fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: ax.imshow(g.predict(g.X)[0]))
-        fig.clf()
-        for i, g in enumerate(self.bgplvms):
-            ax = fig.add_subplot(1, len(self.bgplvms), i + 1)
-            ax.imshow(g.predict(g.X)[0])
-        pylab.draw()
-        fig.tight_layout()
         return fig
 
-    def plot_scales(self, *args, **kwargs):
+    def plot_scales(self, fig_num="MRD Scales", axes=None, *args, **kwargs):
-        fig = pylab.figure("MRD Scales", figsize=(4 * len(self.bgplvms), 3))
+        fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: g.kern.plot_ARD(ax=ax, *args, **kwargs))
-        fig.clf()
-        for i, g in enumerate(self.bgplvms):
-            ax = fig.add_subplot(1, len(self.bgplvms), i + 1)
-            g.kern.plot_ARD(ax=ax, *args, **kwargs)
-        pylab.draw()
-        fig.tight_layout()
         return fig
 
-    def plot_latent(self, *args, **kwargs):
+    def plot_latent(self, fig_num="MRD Latent Spaces", axes=None, *args, **kwargs):
-        fig = pylab.figure("MRD Latent Spaces", figsize=(4 * len(self.bgplvms), 3))
+        fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: g.plot_latent(ax=ax, *args, **kwargs))
-        fig.clf()
-        for i, g in enumerate(self.bgplvms):
-            ax = fig.add_subplot(1, len(self.bgplvms), i + 1)
-            g.plot_latent(ax=ax, *args, **kwargs)
-        pylab.draw()
-        fig.tight_layout()
         return fig
 
     def _debug_plot(self):
-        self.plot_X()
         self.plot_X_1d()
-        self.plot_latent()
+        fig = pylab.figure("MRD DEBUG PLOT", figsize=(4 * len(self.bgplvms), 9))
-        self.plot_scales()
+        fig.clf()
+        axes = [fig.add_subplot(3, len(self.bgplvms), i + 1) for i in range(len(self.bgplvms))]
+        self.plot_X(axes=axes)
+        axes = [fig.add_subplot(3, len(self.bgplvms), i + len(self.bgplvms) + 1) for i in range(len(self.bgplvms))]
+        self.plot_latent(axes=axes)
+        axes = [fig.add_subplot(3, len(self.bgplvms), i + 2 * len(self.bgplvms) + 1) for i in range(len(self.bgplvms))]
+        self.plot_scales(axes=axes)
+        pylab.draw()
+        fig.tight_layout()
 
-    def _debug_optimize(self, opt='scg', maxiters=500, itersteps=10):
+    def _debug_optimize(self, opt='scg', maxiters=5000, itersteps=10):
         iters = 0
-
-        import multiprocessing
-        class M(multiprocessing.Process):
-            def __init__(self, q, *args, **kw):
-                self.q = q
-                super(M, self).__init__(*args, **kw)
-                pass
-            def run(self):
-                pass
-
         optstep = lambda: self.optimize(opt, messages=1, max_f_eval=itersteps)
         self._debug_plot()
         raw_input("enter to start debug")