apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-24 14:15:14 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

[mrd] plotting, init, inference etc.

Browse source
This commit is contained in:
Max Zwiessele 2014-05-16 11:21:08 +01:00
parent efc1f4413c
commit 94c84a23a3
8 changed files with 236 additions and 72 deletions
Download patch file Download diff file Expand all files Collapse all files

82
GPy/core/gp.py
View file

@ -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):
"""

19
GPy/core/parameterization/parameterized.py
View file

@ -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()

6
GPy/core/sparse_gp.py
View file

@ -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)

43
GPy/examples/dimensionality_reduction.py
View file

@ -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

19
GPy/inference/latent_function_inference/__init__.py
View file

@ -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

4
GPy/models/bayesian_gplvm.py
View file

@ -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):
"""

132
GPy/models/mrd.py
View file

@ -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()

3
GPy/testing/pickle_tests.py
View file

@ -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