From 350497c72606f188f83b68588140f0058190559b Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 16 Apr 2013 11:25:51 +0100 Subject: [PATCH] adjusted plotting behaviour in X1d --- GPy/examples/dimensionality_reduction.py | 62 +++++++++---- GPy/models/mrd.py | 112 +++++++++++++++++------ 2 files changed, 128 insertions(+), 46 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 283bdc76..2c7d6bea 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -118,13 +118,13 @@ def mrd_simulation(plot_sim=False): # Y2 -= Y2.mean(0) # make_params = lambda ard: np.hstack([[1], ard, [1, .3]]) - D1, D2, D3, N, M, Q = 6, 7, 8, 150, 18, 5 - x = np.linspace(0, 2 * np.pi, N)[:, None] + D1, D2, D3, N, M, Q = 50, 100, 8, 200, 2, 5 + x = np.linspace(0, 8 * 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.exp(-np.cos(2 * x))) - sS = np.vectorize(lambda x: np.sin(2 * x)) + sS = np.vectorize(lambda x: x * np.sin(2 * x)) s1 = s1(x) s2 = s2(x) @@ -144,20 +144,30 @@ def mrd_simulation(plot_sim=False): S2 = np.hstack([s2, sS]) S3 = np.hstack([s3, sS]) + from GPy.models import mrd + from GPy import kern + reload(mrd); reload(kern) + +# k = kern.rbf(2, ARD=True) + kern.bias(2) + kern.white(2) +# Y1 = np.random.multivariate_normal(np.zeros(N), k.K(S1), D1).T +# 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 + 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 += .1 * np.random.randn(*Y1.shape) - Y2 += .1 * np.random.randn(*Y2.shape) - Y3 += .1 * np.random.randn(*Y3.shape) + Y1 += .5 * np.random.randn(*Y1.shape) + Y2 += .5 * np.random.randn(*Y2.shape) + Y3 += .5 * 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) +# 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, sS] Ylist = [Y1, Y2] @@ -173,21 +183,33 @@ def mrd_simulation(plot_sim=False): ax.plot(x, S, label=lab) ax.legend() for i, Y in enumerate(Ylist): - ax = fig.add_subplot(2, len(Ylist), len(Slist) + i) + ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i) ax.imshow(Y) ax.set_title("Y{}".format(i + 1)) pylab.draw() pylab.tight_layout() - from GPy.models import mrd - from GPy import kern - 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="single", _debug=False) + # k = kern.rbf(Q, ARD=True) + kern.bias(Q) + kern.white(Q) + k = kern.linear(Q, ARD=True) + kern.bias(Q) + kern.white(Q) + m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, initx="concat", _debug=False) m.ensure_default_constraints() - # cstr = "noise|white|variance" - # m.unconstrain(cstr); m.constrain_bounded(cstr, 1e-10, 1.) + for i, Y in enumerate(Ylist): + m.set('{}_noise'.format(i + 1), Y.var() / 100.) + +# import ipdb;ipdb.set_trace() + cstr = "variance" + m.unconstrain(cstr); m.constrain_bounded(cstr, 1e-15, 1.) + +# print "initializing beta" +# cstr = "noise" +# m.unconstrain(cstr); m.constrain_fixed(cstr) +# m.optimize('scg', messages=1, max_f_eval=200) +# +# print "releasing beta" +# cstr = "noise" +# m.unconstrain(cstr); m.constrain_positive(cstr) + m.auto_scale_factor = True diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index 81e46774..943db420 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -5,11 +5,12 @@ Created on 10 Apr 2013 ''' from GPy.core import model from GPy.models.Bayesian_GPLVM import Bayesian_GPLVM -import numpy from GPy.models.sparse_GP import sparse_GP +from GPy.util.linalg import PCA +from scipy import linalg +import numpy import itertools import pylab -from GPy.util.linalg import PCA class MRD(model): """ @@ -23,8 +24,10 @@ class MRD(model): :type names: [str] :param Q: latent dimensionality (will raise :type Q: int - :param init: initialisation method for the latent space - :type init: 'PCA'|'random' + :param initx: initialisation method for the latent space + :type initx: 'PCA'|'random' + :param initz: initialisation method for inducing inputs + :type initz: 'permute'|'random' :param X: Initial latent space :param X_variance: @@ -38,6 +41,7 @@ class MRD(model): :param kernel: kernel to use """ + def __init__(self, *Ylist, **kwargs): if kwargs.has_key("_debug"): self._debug = kwargs['_debug'] @@ -55,24 +59,30 @@ class MRD(model): del kwargs['kernel'] else: k = lambda: None - if kwargs.has_key('init'): - init = kwargs['init'] - del kwargs['init'] + if kwargs.has_key('initx'): + initx = kwargs['initx'] + del kwargs['initx'] else: - init = "PCA" + initx = "PCA" + if kwargs.has_key('initz'): + initz = kwargs['initz'] + del kwargs['initz'] + else: + initz = "permute" try: self.Q = kwargs["Q"] except KeyError: raise ValueError("Need Q for MRD") try: self.M = kwargs["M"] + del kwargs["M"] except KeyError: self.M = 10 self._init = True - 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] + X = self._init_X(initx, Ylist) + Z = self._init_Z(initz, X) + self.bgplvms = [Bayesian_GPLVM(Y, kernel=k(), X=X, Z=Z, M=self.M, **kwargs) for Y in Ylist] del self._init self.gref = self.bgplvms[0] @@ -96,6 +106,26 @@ class MRD(model): if not self._init: raise AttributeError("bgplvm list not initialized") @property + def Z(self): + return self.gref.Z + @Z.setter + def Z(self, Z): + try: + self.propagate_param(Z=Z) + except AttributeError: + if not self._init: + raise AttributeError("bgplvm list not initialized") + @property + def X_variance(self): + return self.gref.X_variance + @X_variance.setter + def X_variance(self, X_var): + try: + self.propagate_param(X_variance=X_var) + 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 @@ -120,6 +150,11 @@ class MRD(model): for g in self.bgplvms: g.__setattr__(key, val) + def randomize(self, initx='concat', initz='permute', *args, **kw): + super(MRD, self).randomize(*args, **kw) + self._init_X(initx, self.Ylist) + self._init_Z(initz, self.X) + 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)], []) @@ -154,14 +189,14 @@ class MRD(model): params = numpy.hstack([X, X_var, Z, numpy.hstack(thetas)]) return params - def _set_var_params(self, g, X, X_var, Z): - g.X = X - g.X_variance = X_var - g.Z = Z - - def _set_kern_params(self, g, p): - g.kern._set_params(p[:g.kern.Nparam]) - g.likelihood._set_params(p[g.kern.Nparam:]) +# def _set_var_params(self, g, X, X_var, Z): +# g.X = X.reshape(self.N, self.Q) +# g.X_variance = X_var.reshape(self.N, self.Q) +# g.Z = Z.reshape(self.M, self.Q) +# +# def _set_kern_params(self, g, p): +# g.kern._set_params(p[:g.kern.Nparam]) +# g.likelihood._set_params(p[g.kern.Nparam:]) def _set_params(self, x): start = 0; end = self.NQ @@ -225,25 +260,50 @@ class MRD(model): self.X = X return X + + def _init_Z(self, init="permute", X=None): + if X is None: + X = self.X + if init in "permute": + Z = numpy.random.permutation(X.copy())[:self.M] + elif init in "random": + Z = numpy.random.randn(self.M, self.Q) * X.var() + self.Z = Z + return Z + def plot_X_1d(self, colors=None): - if colors is None: - colors = pylab.gca()._get_lines.color_cycle - fig = pylab.figure(num="MRD X 1d", figsize=(4 * len(self.bgplvms), (2 * self.X.shape[1]))) + fig = pylab.figure(num="MRD X 1d", figsize=(min(8, (3 * len(self.bgplvms))), min(12, (2 * self.X.shape[1])))) fig.clf() ax1 = fig.add_subplot(self.X.shape[1], 1, 1) + if colors is None: + colors = ax1._get_lines.color_cycle ax1.plot(self.X, c='k', alpha=.3) plots = ax1.plot(self.X.T[0], c=colors.next()) + ax1.fill_between(numpy.arange(self.X.shape[0]), + self.X.T[0] - 2 * numpy.sqrt(self.gref.X_variance.T[0]), + self.X.T[0] + 2 * numpy.sqrt(self.gref.X_variance.T[0]), + facecolor=plots[-1].get_color(), + alpha=.3) + ax1.text(1, 1, r"$\mathbf{{X_{}}}".format(1), + horizontalalignment='right', + verticalalignment='top', + transform=ax1.transAxes) for i in range(self.X.shape[1] - 1): ax = fig.add_subplot(self.X.shape[1], 1, i + 2) ax.plot(self.X, c='k', alpha=.3) plots.extend(ax.plot(self.X.T[i + 1], c=colors.next())) + ax.fill_between(numpy.arange(self.X.shape[0]), + self.X.T[i + 1] - 2 * numpy.sqrt(self.gref.X_variance.T[i + 1]), + self.X.T[i + 1] + 2 * numpy.sqrt(self.gref.X_variance.T[i + 1]), + facecolor=plots[-1].get_color(), + alpha=.3) if i < self.X.shape[1] - 2: ax.set_xticklabels('') ax1.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.) +# 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