From 865e9df255d4e641a91ef433b0c979183a0ba9ce Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 17 Apr 2013 15:45:20 +0100 Subject: [PATCH] BGPLVM still failing, doesn't seem to be numerical : ( --- GPy/examples/dimensionality_reduction.py | 75 +++++++++++++------ GPy/models/Bayesian_GPLVM.py | 34 ++++++--- GPy/models/mrd.py | 92 ++++++++++-------------- 3 files changed, 117 insertions(+), 84 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 04687a35..1ee19e62 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/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)) - s3 = np.vectorize(lambda x: np.sin(2 * x)) - sS = np.vectorize(lambda x:-np.exp(-np.cos(2 * x))) + s2 = np.vectorize(lambda x: np.cos(x)) + s3 = 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): - D1, D2, D3, N, M, Q = 50, 34, 8, 100, 2, 6 +def bgplvm_simulation(burnin='scg', plot_sim=False, max_f_eval=12): + 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.unconstrain(cstr), m.constrain_bounded(cstr, 1e-20, 1.) + m.ensure_default_constraints() + + 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' - m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-12, 1.) - - cstr = 'linear_variance' - m.unconstrain(cstr), m.constrain_positive(cstr) +# cstr = 'variance' +# 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" diff --git a/GPy/models/Bayesian_GPLVM.py b/GPy/models/Bayesian_GPLVM.py index a99f7667..211d21c6 100644 --- a/GPy/models/Bayesian_GPLVM.py +++ b/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))) - - 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):]) + x = np.hstack((self.X.flatten(), self.X_variance.flatten(), sparse_GP._get_params(self))) + return x + 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 diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index 31548d9a..096c9cb9 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -271,14 +271,31 @@ class MRD(model): self.Z = Z return Z - def plot_X_1d(self, colors=None): - fig = pylab.figure(num="MRD X 1d", figsize=(min(8, (3 * len(self.bgplvms))), min(12, (2 * self.X.shape[1])))) + def _handle_plotting(self, fig_num, axes, plotf): + 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): - fig = pylab.figure("MRD X", figsize=(4 * len(self.bgplvms), 3)) - 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() + def plot_X(self, fig_num="MRD Predictions", axes=None): + fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: ax.imshow(g.X)) return fig - def plot_predict(self): - fig = pylab.figure("MRD Predictions", figsize=(4 * len(self.bgplvms), 3)) - 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() + def plot_predict(self, fig_num="MRD Predictions", axes=None): + fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: ax.imshow(g.predict(g.X)[0])) return fig - def plot_scales(self, *args, **kwargs): - fig = pylab.figure("MRD Scales", figsize=(4 * len(self.bgplvms), 3)) - 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() + def plot_scales(self, fig_num="MRD Scales", axes=None, *args, **kwargs): + fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: g.kern.plot_ARD(ax=ax, *args, **kwargs)) return fig - def plot_latent(self, *args, **kwargs): - fig = pylab.figure("MRD Latent Spaces", figsize=(4 * len(self.bgplvms), 3)) - 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() + def plot_latent(self, fig_num="MRD Latent Spaces", axes=None, *args, **kwargs): + fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: g.plot_latent(ax=ax, *args, **kwargs)) return fig def _debug_plot(self): - self.plot_X() self.plot_X_1d() - self.plot_latent() - self.plot_scales() + fig = pylab.figure("MRD DEBUG PLOT", figsize=(4 * len(self.bgplvms), 9)) + 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")