From 94c84a23a3d5c9e1ed93b40e5848a83cf791b83f Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 16 May 2014 11:21:08 +0100 Subject: [PATCH] [mrd] plotting, init, inference etc. --- GPy/core/gp.py | 82 ++++++++--- GPy/core/parameterization/parameterized.py | 19 ++- GPy/core/sparse_gp.py | 6 +- GPy/examples/dimensionality_reduction.py | 43 +++++- .../latent_function_inference/__init__.py | 19 +++ GPy/models/bayesian_gplvm.py | 4 +- GPy/models/mrd.py | 132 +++++++++++++----- GPy/testing/pickle_tests.py | 3 +- 8 files changed, 236 insertions(+), 72 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 29d9032a..ca2583e0 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -180,40 +180,80 @@ class GP(Model): return Ysim - def plot_f(self, *args, **kwargs): + def plot_f(self, plot_limits=None, which_data_rows='all', + which_data_ycols='all', fixed_inputs=[], + levels=20, samples=0, fignum=None, ax=None, resolution=None, + plot_raw=True, + linecol=None,fillcol=None, Y_metadata=None, data_symbol='kx'): """ - - Plot the GP's view of the world, where the data is normalized and - before applying a likelihood. - - This is a convenience function: arguments are passed to - GPy.plotting.matplot_dep.models_plots.plot_f_fit - + Plot the GP's view of the world, where the data is normalized and before applying a likelihood. + This is a call to plot with plot_raw=True. + Data will not be plotted in this, as the GP's view of the world + may live in another space, or units then the data. """ assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ..plotting.matplot_dep import models_plots - return models_plots.plot_fit_f(self,*args,**kwargs) + kw = {} + if linecol is not None: + kw['linecol'] = linecol + if fillcol is not None: + kw['fillcol'] = fillcol + return models_plots.plot_fit(self, plot_limits, which_data_rows, + which_data_ycols, fixed_inputs, + levels, samples, fignum, ax, resolution, + plot_raw=plot_raw, Y_metadata=Y_metadata, + data_symbol=data_symbol, **kw) - def plot(self, *args, **kwargs): + def plot(self, plot_limits=None, which_data_rows='all', + which_data_ycols='all', fixed_inputs=[], + levels=20, samples=0, fignum=None, ax=None, resolution=None, + plot_raw=False, + linecol=None,fillcol=None, Y_metadata=None, data_symbol='kx'): """ Plot the posterior of the GP. - - In one dimension, the function is plotted with a shaded region - identifying two standard deviations. - - In two dimsensions, a contour-plot shows the mean predicted - function - - In higher dimensions, use fixed_inputs to plot the GP with some of - the inputs fixed. + - In one dimension, the function is plotted with a shaded region identifying two standard deviations. + - In two dimsensions, a contour-plot shows the mean predicted function + - In higher dimensions, use fixed_inputs to plot the GP with some of the inputs fixed. Can plot only part of the data and part of the posterior functions - using which_data_rows which_data_ycols and which_parts - - This is a convenience function: arguments are passed to - GPy.plotting.matplot_dep.models_plots.plot_fit + using which_data_rowsm which_data_ycols. + :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits + :type plot_limits: np.array + :param which_data_rows: which of the training data to plot (default all) + :type which_data_rows: 'all' or a slice object to slice model.X, model.Y + :param which_data_ycols: when the data has several columns (independant outputs), only plot these + :type which_data_rows: 'all' or a list of integers + :param fixed_inputs: a list of tuple [(i,v), (i,v)...], specifying that input index i should be set to value v. + :type fixed_inputs: a list of tuples + :param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D + :type resolution: int + :param levels: number of levels to plot in a contour plot. + :type levels: int + :param samples: the number of a posteriori samples to plot + :type samples: int + :param fignum: figure to plot on. + :type fignum: figure number + :param ax: axes to plot on. + :type ax: axes handle + :type output: integer (first output is 0) + :param linecol: color of line to plot [Tango.colorsHex['darkBlue']] + :type linecol: + :param fillcol: color of fill [Tango.colorsHex['lightBlue']] + :param levels: for 2D plotting, the number of contour levels to use is ax is None, create a new figure """ assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ..plotting.matplot_dep import models_plots - return models_plots.plot_fit(self,*args,**kwargs) + kw = {} + if linecol is not None: + kw['linecol'] = linecol + if fillcol is not None: + kw['fillcol'] = fillcol + return models_plots.plot_fit(self, plot_limits, which_data_rows, + which_data_ycols, fixed_inputs, + levels, samples, fignum, ax, resolution, + plot_raw=plot_raw, Y_metadata=Y_metadata, + data_symbol=data_symbol, **kw) def input_sensitivity(self): """ diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index cfb2171d..48eb2ddc 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -272,8 +272,11 @@ class Parameterized(Parameterizable): def __setattr__(self, name, val): # override the default behaviour, if setting a param, so broadcasting can by used if hasattr(self, "parameters"): - pnames = self.parameter_names(False, adjust_for_printing=True, recursive=False) - if name in pnames: self.parameters[pnames.index(name)][:] = val; return + try: + pnames = self.parameter_names(False, adjust_for_printing=True, recursive=False) + if name in pnames: self.parameters[pnames.index(name)][:] = val; return + except AttributeError: + pass object.__setattr__(self, name, val); #=========================================================================== @@ -281,11 +284,15 @@ class Parameterized(Parameterizable): #=========================================================================== def __setstate__(self, state): super(Parameterized, self).__setstate__(state) - self._connect_parameters() - self._connect_fixes() - self._notify_parent_change() + try: + self._connect_parameters() + self._connect_fixes() + self._notify_parent_change() + + self.parameters_changed() + except Exception as e: + print "WARNING: caught exception {!s}, trying to continue".format(e) - self.parameters_changed() def copy(self): c = super(Parameterized, self).copy() c._connect_parameters() diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 388b682a..ac752c36 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -66,7 +66,11 @@ class SparseGP(GP): #gradients wrt Z self.Z.gradient = self.kern.gradients_X(dL_dKmm, self.Z) self.Z.gradient += self.kern.gradients_Z_expectations( - self.grad_dict['dL_dpsi0'], self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X) + self.grad_dict['dL_dpsi0'], + self.grad_dict['dL_dpsi1'], + self.grad_dict['dL_dpsi2'], + Z=self.Z, + variational_posterior=self.X) else: #gradients wrt kernel self.kern.update_gradients_diag(self.grad_dict['dL_dKdiag'], self.X) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 8acf5489..2a319e8a 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -303,9 +303,11 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1, k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) inan = _np.random.binomial(1, .6, size=Y.shape).astype(bool) - m = BayesianGPLVM(Y.copy(), Q, init="random", num_inducing=num_inducing, kernel=k) - m.inference_method = VarDTCMissingData() - m.Y[inan] = _np.nan + Y[inan] = _np.nan + + m = BayesianGPLVM(Y.copy(), Q, init="random", num_inducing=num_inducing, + inference_method=VarDTCMissingData(inan=inan), kernel=k) + m.X.variance[:] = _np.random.uniform(0,.01,m.X.shape) m.likelihood.variance = .01 m.parameters_changed() @@ -338,7 +340,40 @@ def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw): print "Optimizing Model:" m.optimize(messages=verbose, max_iters=8e3, gtol=.1) if plot: - m.plot_X_1d("MRD Latent Space 1D") + m.X.plot("MRD Latent Space 1D") + m.plot_scales("MRD Scales") + return m + +def mrd_simulation_missing_data(optimize=True, verbose=True, plot=True, plot_sim=True, **kw): + from GPy import kern + from GPy.models import MRD + from GPy.inference.latent_function_inference.var_dtc import VarDTCMissingData + + D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5 + _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) + + #Ylist = [Ylist[0]] + k = kern.Linear(Q, ARD=True) + inanlist = [] + + for Y in Ylist: + inan = _np.random.binomial(1, .6, size=Y.shape).astype(bool) + inanlist.append(inan) + Y[inan] = _np.nan + + imlist = [] + for inan in inanlist: + imlist.append(VarDTCMissingData(limit=1, inan=inan)) + + m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing, + kernel=k, inference_method=imlist, + initx="random", initz='permute', **kw) + + if optimize: + print "Optimizing Model:" + m.optimize('bfgs', messages=verbose, max_iters=8e3, gtol=.1) + if plot: + m.X.plot("MRD Latent Space 1D") m.plot_scales("MRD Scales") return m diff --git a/GPy/inference/latent_function_inference/__init__.py b/GPy/inference/latent_function_inference/__init__.py index 878c1e4f..8589a60e 100644 --- a/GPy/inference/latent_function_inference/__init__.py +++ b/GPy/inference/latent_function_inference/__init__.py @@ -38,6 +38,25 @@ class LatentFunctionInference(object): """ pass +class InferenceMethodList(LatentFunctionInference, list): + + def on_optimization_start(self): + for inf in self: + inf.on_optimization_start() + + def on_optimization_end(self): + for inf in self: + inf.on_optimization_end() + + def __getstate__(self): + state = [] + for inf in self: + state.append(inf) + return state + + def __setstate__(self, state): + for inf in state: + self.append(inf) from exact_gaussian_inference import ExactGaussianInference from laplace import Laplace diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 2bcbe0b2..26f4cd53 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -83,7 +83,7 @@ class BayesianGPLVM(SparseGP): resolution=50, ax=None, marker='o', s=40, fignum=None, plot_inducing=True, legend=True, plot_limits=None, - aspect='auto', updates=False, **kwargs): + aspect='auto', updates=False, predict_kwargs={}, imshow_kwargs={}): import sys assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ..plotting.matplot_dep import dim_reduction_plots @@ -91,7 +91,7 @@ class BayesianGPLVM(SparseGP): return dim_reduction_plots.plot_latent(self, labels, which_indices, resolution, ax, marker, s, fignum, plot_inducing, legend, - plot_limits, aspect, updates, **kwargs) + plot_limits, aspect, updates, predict_kwargs, imshow_kwargs) def do_test_latents(self, Y): """ diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index 458a70a1..73b267ba 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -11,9 +11,11 @@ from ..core.parameterization.variational import NormalPosterior, NormalPrior from ..core.parameterization import Param, Parameterized from ..inference.latent_function_inference.var_dtc import VarDTCMissingData, VarDTC from ..likelihoods import Gaussian -from GPy.util.initialization import initialize_latent +from ..util.initialization import initialize_latent +from ..core.sparse_gp import SparseGP, GP +from ..inference.latent_function_inference import InferenceMethodList -class MRD(Model): +class MRD(SparseGP): """ Apply MRD to all given datasets Y in Ylist. @@ -39,17 +41,18 @@ class MRD(Model): :param num_inducing: number of inducing inputs to use :param Z: initial inducing inputs :param kernel: list of kernels or kernel to copy for each output - :type kernel: [GPy.kern.kern] | GPy.kern.kern | None (default) - :param :class:`~GPy.inference.latent_function_inference inference_method: the inference method to use - :param :class:`~GPy.likelihoods.likelihood.Likelihood` likelihood: the likelihood to use + :type kernel: [GPy.kernels.kernels] | GPy.kernels.kernels | None (default) + :param :class:`~GPy.inference.latent_function_inference inference_method: + InferenceMethodList of inferences, or one inference method for all + :param :class:`~GPy.likelihoodss.likelihoods.likelihoods` likelihoods: the likelihoods to use :param str name: the name of this model :param [str] Ynames: the names for the datasets given, must be of equal length as Ylist or None """ def __init__(self, Ylist, input_dim, X=None, X_variance=None, initx = 'PCA', initz = 'permute', num_inducing=10, Z=None, kernel=None, - inference_method=None, likelihood=None, name='mrd', Ynames=None): - super(MRD, self).__init__(name) + inference_method=None, likelihoods=None, name='mrd', Ynames=None): + super(GP, self).__init__(name) self.input_dim = input_dim self.num_inducing = num_inducing @@ -63,16 +66,16 @@ class MRD(Model): # sort out the kernels if kernel is None: from ..kern import RBF - self.kern = [RBF(input_dim, ARD=1, lengthscale=fracs[i], name='rbf'.format(i)) for i in range(len(Ylist))] + self.kernels = [RBF(input_dim, ARD=1, lengthscale=fracs[i], name='rbf'.format(i)) for i in range(len(Ylist))] elif isinstance(kernel, Kern): - self.kern = [] + self.kernels = [] for i in range(len(Ylist)): k = kernel.copy() - self.kern.append(k) + self.kernels.append(k) else: assert len(kernel) == len(Ylist), "need one kernel per output" assert all([isinstance(k, Kern) for k in kernel]), "invalid kernel object detected!" - self.kern = kernel + self.kernels = kernel if X_variance is None: X_variance = np.random.uniform(0.1, 0.2, X.shape) @@ -80,32 +83,44 @@ class MRD(Model): self.variational_prior = NormalPrior() self.X = NormalPosterior(X, X_variance) - if likelihood is None: - self.likelihood = [Gaussian(name='Gaussian_noise'.format(i)) for i in range(len(Ylist))] - else: self.likelihood = likelihood + if likelihoods is None: + self.likelihoods = [Gaussian(name='Gaussian_noise'.format(i)) for i in range(len(Ylist))] + else: self.likelihoods = likelihoods if inference_method is None: - self.inference_method= [] + self.inference_method= InferenceMethodList() for y in Ylist: - if np.any(np.isnan(y)): - self.inference_method.append(VarDTCMissingData(limit=1)) + inan = np.isnan(y) + if np.any(inan): + self.inference_method.append(VarDTCMissingData(limit=1, inan=inan)) else: self.inference_method.append(VarDTC(limit=1)) else: + if not isinstance(inference_method, InferenceMethodList): + inference_method = InferenceMethodList(inference_method) self.inference_method = inference_method - self.inference_method.set_limit(len(Ylist)) self.add_parameters(self.X, self.Z) if Ynames is None: Ynames = ['Y{}'.format(i) for i in range(len(Ylist))] + self.names = Ynames + + self.bgplvms = [] + self.num_data = Ylist[0].shape[0] + + for i, n, k, l, Y in itertools.izip(itertools.count(), Ynames, self.kernels, self.likelihoods, self.Ylist): + assert Y.shape[0] == self.num_data, "All datasets need to share the number of datapoints, and those have to correspond to one another" - for i, n, k, l in itertools.izip(itertools.count(), Ynames, self.kern, self.likelihood): p = Parameterized(name=n) p.add_parameter(k) + p.kern = k p.add_parameter(l) - setattr(self, 'Y{}'.format(i), p) + p.likelihood = l self.add_parameter(p) + self.bgplvms.append(p) + + self.posterior = None self._in_init_ = False def parameters_changed(self): @@ -114,13 +129,13 @@ class MRD(Model): self.Z.gradient[:] = 0. self.X.gradient[:] = 0. - for y, k, l, i in itertools.izip(self.Ylist, self.kern, self.likelihood, self.inference_method): + for y, k, l, i in itertools.izip(self.Ylist, self.kernels, self.likelihoods, self.inference_method): posterior, lml, grad_dict = i.inference(k, self.X, self.Z, l, y) self.posteriors.append(posterior) self._log_marginal_likelihood += lml - # likelihood gradients + # likelihoods gradients l.update_gradients(grad_dict.pop('dL_dthetaL')) #gradients wrt kernel @@ -133,13 +148,20 @@ class MRD(Model): #gradients wrt Z self.Z.gradient += k.gradients_X(dL_dKmm, self.Z) self.Z.gradient += k.gradients_Z_expectations( - grad_dict['dL_dpsi1'], grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X) + grad_dict['dL_dpsi0'], + grad_dict['dL_dpsi1'], + grad_dict['dL_dpsi2'], + Z=self.Z, variational_posterior=self.X) dL_dmean, dL_dS = k.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, **grad_dict) self.X.mean.gradient += dL_dmean self.X.variance.gradient += dL_dS # update for the KL divergence + self.posterior = self.posteriors[0] + self.kern = self.kernels[0] + self.likelihood = self.likelihoods[0] + self.variational_prior.update_gradients_KL(self.X) self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X) @@ -207,16 +229,25 @@ class MRD(Model): else: return pylab.gcf() - def plot_X(self, fignum=None, ax=None): - fig = self._handle_plotting(fignum, ax, lambda i, g, ax: ax.imshow(g.X)) - return fig + def predict(self, Xnew, full_cov=False, Y_metadata=None, kern=None, Yindex=0): + """ + Prediction for data set Yindex[default=0]. + This predicts the output mean and variance for the dataset given in Ylist[Yindex] + """ + self.posterior = self.posteriors[Yindex] + self.kern = self.kernels[Yindex] + self.likelihood = self.likelihoods[Yindex] + return super(MRD, self).predict(Xnew, full_cov, Y_metadata, kern) - def plot_predict(self, fignum=None, ax=None, sharex=False, sharey=False, **kwargs): - fig = self._handle_plotting(fignum, - ax, - lambda i, g, ax: ax.imshow(g. predict(g.X)[0], **kwargs), - sharex=sharex, sharey=sharey) - return fig + #=============================================================================== + # TODO: Predict! Maybe even change to several bgplvms, which share an X? + #=============================================================================== + # def plot_predict(self, fignum=None, ax=None, sharex=False, sharey=False, **kwargs): + # fig = self._handle_plotting(fignum, + # ax, + # lambda i, g, ax: ax.imshow(g.predict(g.X)[0], **kwargs), + # sharex=sharex, sharey=sharey) + # return fig def plot_scales(self, fignum=None, ax=None, titles=None, sharex=False, sharey=True, *args, **kwargs): """ @@ -228,16 +259,28 @@ class MRD(Model): """ if titles is None: titles = [r'${}$'.format(name) for name in self.names] - ymax = reduce(max, [np.ceil(max(g.input_sensitivity())) for g in self.bgplvms]) + ymax = reduce(max, [np.ceil(max(g.kernels.input_sensitivity())) for g in self.bgplvms]) def plotf(i, g, ax): ax.set_ylim([0,ymax]) - g.kern.plot_ARD(ax=ax, title=titles[i], *args, **kwargs) + g.kernels.plot_ARD(ax=ax, title=titles[i], *args, **kwargs) fig = self._handle_plotting(fignum, ax, plotf, sharex=sharex, sharey=sharey) return fig - def plot_latent(self, fignum=None, ax=None, *args, **kwargs): - fig = self.gref.plot_latent(fignum=fignum, ax=ax, *args, **kwargs) # self._handle_plotting(fignum, ax, lambda i, g, ax: g.plot_latent(ax=ax, *args, **kwargs)) - return fig + def plot_latent(self, labels=None, which_indices=None, + resolution=50, ax=None, marker='o', s=40, + fignum=None, plot_inducing=True, legend=True, + plot_limits=None, + aspect='auto', updates=False, predict_kwargs=dict(Yindex=0), imshow_kwargs={}): + """ + """ + import sys + assert "matplotlib" in sys.modules, "matplotlib package has not been imported." + from ..plotting.matplot_dep import dim_reduction_plots + + return dim_reduction_plots.plot_latent(self, labels, which_indices, + resolution, ax, marker, s, + fignum, plot_inducing, legend, + plot_limits, aspect, updates, predict_kwargs, imshow_kwargs) def _debug_plot(self): self.plot_X_1d() @@ -252,4 +295,19 @@ class MRD(Model): pylab.draw() fig.tight_layout() + def __getstate__(self): + # TODO: + import copy + state = copy.copy(self.__dict__) + del state['kernels'] + del state['kern'] + del state['likelihood'] + return state + def __setstate__(self, state): + # TODO: + super(MRD, self).__setstate__(state) + self.kernels = [p.kern for p in self.bgplvms] + self.kern = self.kernels[0] + self.likelihood = self.likelihoods[0] + self.parameters_changed() \ No newline at end of file diff --git a/GPy/testing/pickle_tests.py b/GPy/testing/pickle_tests.py index 94504645..13a67744 100644 --- a/GPy/testing/pickle_tests.py +++ b/GPy/testing/pickle_tests.py @@ -4,7 +4,8 @@ Created on 13 Mar 2014 @author: maxz ''' import unittest, itertools -import cPickle as pickle +#import cPickle as pickle +import pickle import numpy as np from GPy.core.parameterization.index_operations import ParameterIndexOperations,\ ParameterIndexOperationsView