From 9e83d8162ec3ca406ab08c20fc3bade971a28969 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 22 May 2013 17:14:54 +0100 Subject: [PATCH 1/3] changed likelihood and kernel handling --- GPy/models/mrd.py | 121 ++++++++++++---------------------------------- 1 file changed, 31 insertions(+), 90 deletions(-) diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index 74f526d0..87bcffbf 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -11,6 +11,7 @@ from scipy import linalg import numpy import itertools import pylab +from GPy.kern.kern import kern class MRD(model): """ @@ -19,7 +20,7 @@ class MRD(model): N must be shared across datasets though. :param likelihood_list...: likelihoods of observed datasets - :type likelihood_list: [GPy.likelihood] + :type likelihood_list: [GPy.likelihood] | [Y1..Yy] :param names: names for different gplvm models :type names: [str] :param Q: latent dimensionality (will raise @@ -38,85 +39,32 @@ class MRD(model): number of inducing inputs to use :param Z: initial inducing inputs - :param kernel: - kernel to use + :param kernels: list of kernels or kernel shared for all BGPLVMS + :type kernels: [GPy.kern.kern] | GPy.kern.kern | None (default) """ - def __init__(self,likelihood_list,Q,M=10,names=None,kernels=None,initX='PCA',initz='permute',_debug=False, **kwargs): + def __init__(self, likelihood_or_Y_list, Q, M=10, names=None, + kernels=None, initx='PCA', + initz='permute', _debug=False, **kwargs): if names is None: - self.names = ["{}".format(i + 1) for i in range(len(likelihood_list))] + self.names = ["{}".format(i + 1) for i in range(len(likelihood_or_Y_list))] - #sort out the kernels + # sort out the kernels if kernels is None: - kernels = [None]*len(likelihood_list) - elif isinstance(kernels,kern.kern): - kernels = [kernels.copy() for i in range(len(likelihood_list))] + kernels = [None] * len(likelihood_or_Y_list) + elif isinstance(kernels, kern): + kernels = [kernels.copy() for i in range(len(likelihood_or_Y_list))] else: - assert len(kernels)==len(likelihood_list), "need one kernel per output" - assert all([isinstance(k, kern.kern) for k in kernels]), "invalid kernel object detected!" + assert len(kernels) == len(likelihood_or_Y_list), "need one kernel per output" + assert all([isinstance(k, kern) for k in kernels]), "invalid kernel object detected!" self.Q = Q self.M = M - self.N = self.gref.N - self.NQ = self.N * self.Q - self.MQ = self.M * self.Q + self._debug = _debug self._init = True - X = self._init_X(initx, likelihood_list) + X = self._init_X(initx, likelihood_or_Y_list) Z = self._init_Z(initz, X) - self.bgplvms = [Bayesian_GPLVM(l, k, X=X, Z=Z, M=self.M, **kwargs) for l,k in zip(likelihood_list,kernels)] - - 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]) - self.nparams = nparams.cumsum() - - model.__init__(self) # @UndefinedVariable - - - - def __init__(self, *likelihood_list, **kwargs): - if kwargs.has_key("_debug"): - self._debug = kwargs['_debug'] - del kwargs['_debug'] - else: - self._debug = False - if kwargs.has_key("names"): - self.names = kwargs['names'] - del kwargs['names'] - else: - self.names = ["{}".format(i + 1) for i in range(len(likelihood_list))] - if kwargs.has_key('kernel'): - kernel = kwargs['kernel'] - k = lambda: kernel.copy() - del kwargs['kernel'] - else: - k = lambda: None - if kwargs.has_key('initx'): - initx = kwargs['initx'] - del kwargs['initx'] - else: - 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(initx, likelihood_list) - Z = self._init_Z(initz, X) - self.bgplvms = [Bayesian_GPLVM(Y, kernel=k(), X=X, Z=Z, M=self.M, **kwargs) for Y in likelihood_list] - + self.bgplvms = [Bayesian_GPLVM(l, k, X=X, Z=Z, M=self.M, **kwargs) for l, k in zip(likelihood_or_Y_list, kernels)] del self._init self.gref = self.bgplvms[0] @@ -212,13 +160,6 @@ class MRD(model): X = self.gref.X.ravel() X_var = self.gref.X_variance.ravel() Z = self.gref.Z.ravel() - - if self._debug: - for g in self.bgplvms: - assert numpy.allclose(g.X.ravel(), X) - assert numpy.allclose(g.X_variance.ravel(), X_var) - assert numpy.allclose(g.Z.ravel(), Z) - thetas = [sparse_GP._get_params(g)[g.Z.size:] for g in self.bgplvms] params = numpy.hstack([X, X_var, Z, numpy.hstack(thetas)]) return params @@ -241,12 +182,6 @@ class MRD(model): Z = x[start:end] thetas = x[end:] - if self._debug: - for g in self.bgplvms: - assert numpy.allclose(g.X, self.gref.X) - assert numpy.allclose(g.X_variance, self.gref.X_variance) - assert numpy.allclose(g.Z, self.gref.Z) - # set params for all: for g, s, e in itertools.izip(self.bgplvms, self.nparams, self.nparams[1:]): g._set_params(numpy.hstack([X, X_var, Z, thetas[s:e]])) @@ -259,7 +194,7 @@ class MRD(model): # g._computations() - def update_likelihood_approximation(self):#TODO: object oriented vs script base + def update_likelihood_approximation(self): # TODO: object oriented vs script base for bgplvm in self.bgplvms: bgplvm.update_likelihood_approximation() @@ -287,15 +222,21 @@ class MRD(model): def _init_X(self, init='PCA', likelihood_list=None): if likelihood_list is None: likelihood_list = self.likelihood_list + Ylist = [] + for likelihood_or_Y in likelihood_list: + if type(likelihood_or_Y) is numpy.ndarray: + Ylist.append(likelihood_or_Y) + else: + Ylist.append(likelihood_or_Y.Y) + del likelihood_list if init in "PCA_single": - X = numpy.zeros((likelihood_list[0].Y.shape[0], self.Q)) - for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.Q), len(likelihood_list)), likelihood_list): - X[:, qs] = PCA(Y.Y, len(qs))[0] + 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([l.Y for l in likelihood_list]), self.Q)[0] - #X = PCA(numpy.hstack(likelihood_list), self.Q)[0] + X = PCA(numpy.hstack(Ylist), self.Q)[0] else: # init == 'random': - X = numpy.random.randn(likelihood_list[0].Y.shape[0], self.Q) + X = numpy.random.randn(Ylist[0].shape[0], self.Q) self.X = X return X @@ -334,7 +275,7 @@ class MRD(model): return fig def plot_predict(self, fig_num="MRD Predictions", axes=None, **kwargs): - fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: ax.imshow(g.predict(g.X)[0],**kwargs)) + fig = self._handle_plotting(fig_num, axes, lambda i, g, ax: ax.imshow(g.predict(g.X)[0], **kwargs)) return fig def plot_scales(self, fig_num="MRD Scales", axes=None, *args, **kwargs): From a7692678c944ce2772a9ccaf6e02e51c289cee96 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 22 May 2013 17:16:10 +0100 Subject: [PATCH 2/3] logexp_clipped bounding behaviour --- GPy/core/transformations.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/core/transformations.py b/GPy/core/transformations.py index fcbfb548..c2408cb8 100644 --- a/GPy/core/transformations.py +++ b/GPy/core/transformations.py @@ -49,8 +49,8 @@ class logexp_clipped(transformation): def f(self, x): exp = np.exp(np.clip(x, self.log_min_bound, self.log_max_bound)) f = np.log(1. + exp) - if np.isnan(f).any(): - import ipdb;ipdb.set_trace() +# if np.isnan(f).any(): +# import ipdb;ipdb.set_trace() return f def finv(self, f): return np.log(np.exp(np.clip(f, self.min_bound, self.max_bound)) - 1.) From 26d68861ca35300292b68d6497c5c680fcda5fc6 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 22 May 2013 17:17:20 +0100 Subject: [PATCH 3/3] minor changes, dimensionality reduction tests --- GPy/examples/dimensionality_reduction.py | 18 +++++++++++------- GPy/inference/SCG.py | 17 +++++++++++------ GPy/likelihoods/Gaussian.py | 2 +- GPy/models/Bayesian_GPLVM.py | 6 +++++- GPy/models/GPLVM.py | 2 +- 5 files changed, 29 insertions(+), 16 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 1f30a67f..d7064f6b 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -266,15 +266,19 @@ def bgplvm_simulation(optimize='scg', if optimize: print "Optimizing model:" - m.optimize('scg', max_iters=max_f_eval, max_f_eval=max_f_eval, messages=True) + m.optimize('bfgs', max_iters=max_f_eval, + max_f_eval=max_f_eval, + messages=True, gtol=1e-2) if plot: import pylab m.plot_X_1d() - pylab.figure(); pylab.axis(); m.kern.plot_ARD() + pylab.figure('BGPLVM Simulation ARD Parameters'); + pylab.axis(); + m.kern.plot_ARD() return m def mrd_simulation(optimize=True, plot_sim=False, **kw): - D1, D2, D3, N, M, Q = 150, 250, 30, 200, 3, 7 + D1, D2, D3, N, M, Q = 15, 8, 8, 100, 3, 7 slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim) from GPy.models import mrd @@ -284,20 +288,20 @@ def mrd_simulation(optimize=True, plot_sim=False, **kw): reload(mrd); reload(kern) k = kern.linear(Q, [0.01] * Q, True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) - m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, initx="concat", initz='permute', **kw) + m = mrd.MRD(Ylist, Q=Q, M=M, kernel=k, initx="concat", initz='permute', **kw) for i, Y in enumerate(Ylist): m['{}_noise'.format(i + 1)] = Y.var() / 100. - # m.constrain('variance|noise', logexp_clipped()) + m.constrain('variance|noise', logexp_clipped()) m.ensure_default_constraints() # DEBUG - np.seterr("raise") + # np.seterr("raise") if optimize: print "Optimizing Model:" - m.optimize('scg', messages=1, max_iters=3e3) + m.optimize('bfgs', messages=1, max_iters=3e3) return m diff --git a/GPy/inference/SCG.py b/GPy/inference/SCG.py index d0a30f0d..5b210953 100644 --- a/GPy/inference/SCG.py +++ b/GPy/inference/SCG.py @@ -25,6 +25,13 @@ import numpy as np import sys + +def print_out(len_maxiters, display, fnow, current_grad, beta, iteration): + if display: + print '\r', + print '{0:>0{mi}g} {1:> 12e} {2:> 12e} {3:> 12e}'.format(iteration, float(fnow), float(beta), float(current_grad), mi=len_maxiters), # print 'Iteration:', iteration, ' Objective:', fnow, ' Scale:', beta, '\r', + sys.stdout.flush() + def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xtol=None, ftol=None, gtol=None): """ Optimisation through Scaled Conjugate Gradients (SCG) @@ -65,7 +72,8 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto iteration = 0 if display: - print ' {0:{mi}s} {1:11s} {2:11s} {3:11s}'.format("I", "F", "Scale", "|g|", mi=len(str(maxiters))) + len_maxiters = len(str(maxiters)) + print ' {0:{mi}s} {1:11s} {2:11s} {3:11s}'.format("I", "F", "Scale", "|g|", mi=len_maxiters) # Main optimization loop. while iteration < maxiters: @@ -113,11 +121,7 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto flog.append(fnow) # Current function value iteration += 1 - if display: - print '\r', - print '{0:>0{mi}g} {1:> 12e} {2:> 12e} {3:> 12e}'.format(iteration, float(fnow), float(beta), float(current_grad), mi=len(str(maxiters))), - # print 'Iteration:', iteration, ' Objective:', fnow, ' Scale:', beta, '\r', - sys.stdout.flush() + print_out(len_maxiters, display, fnow, current_grad, beta, iteration) if success: # Test for termination @@ -158,5 +162,6 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto status = "maxiter exceeded" if display: + print_out(len_maxiters, display, fnow, current_grad, beta, iteration) print "" return x, flog, function_eval, status diff --git a/GPy/likelihoods/Gaussian.py b/GPy/likelihoods/Gaussian.py index 98227b16..e08fee90 100644 --- a/GPy/likelihoods/Gaussian.py +++ b/GPy/likelihoods/Gaussian.py @@ -53,7 +53,7 @@ class Gaussian(likelihood): def _set_params(self, x): x = float(x) if self._variance != x: - self.precision = 1. / max(x, 1e-6) + self.precision = 1. / x self.covariance_matrix = np.eye(self.N) * x self.V = (self.precision) * self.Y self._variance = x diff --git a/GPy/models/Bayesian_GPLVM.py b/GPy/models/Bayesian_GPLVM.py index e1e99af9..cc8eede4 100644 --- a/GPy/models/Bayesian_GPLVM.py +++ b/GPy/models/Bayesian_GPLVM.py @@ -93,8 +93,12 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): x = np.hstack((self.X.flatten(), self.X_variance.flatten(), sparse_GP._get_params(self))) return x + def _clipped(self, x): + return x # np.clip(x, -1e100, 1e100) + def _set_params(self, x, save_old=True, save_count=0): # try: + x = self._clipped(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() @@ -176,7 +180,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): # ======================== self.dbound_dmuS = np.hstack((d_dmu, d_dS)) self.dbound_dZtheta = sparse_GP._log_likelihood_gradients(self) - return np.hstack((self.dbound_dmuS.flatten(), self.dbound_dZtheta)) + return self._clipped(np.hstack((self.dbound_dmuS.flatten(), self.dbound_dZtheta))) def plot_latent(self, which_indices=None, *args, **kwargs): diff --git a/GPy/models/GPLVM.py b/GPy/models/GPLVM.py index ff2be732..2525ddf9 100644 --- a/GPy/models/GPLVM.py +++ b/GPy/models/GPLVM.py @@ -119,7 +119,7 @@ class GPLVM(GP): else: x = self.X[index,input_1] y = self.X[index,input_2] - ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), mew=1.3, label=this_label) + ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label) ax.set_xlabel('latent dimension %i'%input_1) ax.set_ylabel('latent dimension %i'%input_2)