From c63beddcf0e41663ba2d30a6864251dd4dab44a3 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 15 Apr 2013 10:57:27 +0100 Subject: [PATCH] new functions mrd init_X update --- GPy/examples/dimensionality_reduction.py | 81 ++++++++++++++++-------- GPy/inference/SCG.py | 2 +- GPy/models/__init__.py | 4 +- GPy/models/mrd.py | 75 +++++++++++++++++----- 4 files changed, 115 insertions(+), 47 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 2548a040..d578a0c4 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/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) - ard2 = np.array([0., 1, 1, 0], dtype=float) - ard1[ard1 == 0] = 1E-10 - ard2[ard2 == 0] = 1E-10 +# ard1 = np.array([1., 1, 0, 0], dtype=float) +# ard2 = np.array([0., 1, 1, 0], dtype=float) +# ard1[ard1 == 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 - x = np.linspace(0, 4 * np.pi, N)[:, None] + D1, D2, D3, N, M, Q = 5, 5, 5, 150, 18, 5 + 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) - S2 = np.hstack([s2(x), sS(x)]) + .1 * np.random.randn(N, 2) + s1 = s1(x) + 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") - ax = fig.add_subplot(2, 2, 1) - ax.imshow(S1) - ax.set_title("S1") - ax = fig.add_subplot(2, 2, 2) - ax.imshow(S2) - ax.set_title("S2") - ax = fig.add_subplot(2, 2, 3) - ax.imshow(Y1) - ax.set_title("Y1") - ax = fig.add_subplot(2, 2, 4) - ax.imshow(Y2) - ax.set_title("Y2") + import itertools + fig = pylab.figure("MRD Simulation", figsize=(12, 12)) + fig.clf() + ax = fig.add_subplot(2, 1, 1) + labls = sorted(filter(lambda x: x.startswith("s"), locals())) + for S, lab in itertools.izip(Slist, labls): + 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.imshow(Y) + ax.set_title("Y{}".format(i + 1)) pylab.draw() pylab.tight_layout() - k = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q) + GPy.kern.white(Q) - - m = MRD(Y1, Y2, Q=Q, M=M, kernel=k, init="concat", _debug=False) + 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="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) diff --git a/GPy/inference/SCG.py b/GPy/inference/SCG.py index 6531056b..ca42acfe 100644 --- a/GPy/inference/SCG.py +++ b/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() diff --git a/GPy/models/__init__.py b/GPy/models/__init__.py index 91cc60e3..d63adaf1 100644 --- a/GPy/models/__init__.py +++ b/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 -MRD = mrd.MRD -del mrd +from mrd import MRD from generalized_FITC import generalized_FITC diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index fbebd59f..53825a50 100644 --- a/GPy/models/mrd.py +++ b/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 - - X = self._init_X(Ylist, init) + 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] + 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) - self._set_kern_params(g, thetas[s:e].copy()) - g._compute_kernel_matrices() - g._computations() + g._set_params(numpy.hstack([X, X_var, Z, thetas[s:e]])) +# self._set_var_params(g, X, X_var, Z) +# self._set_kern_params(g, thetas[s:e].copy()) +# 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'): - if init in "PCA_concat": - X = PCA(numpy.hstack(Ylist), self.Q)[0] - elif init in "PCA_single": + def _init_X(self, init='PCA', Ylist=None): + if Ylist is None: + Ylist = self.Ylist + 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: