plotting conflict fixed

This commit is contained in:
James Hensman 2014-02-24 13:55:11 +00:00
commit d3eaef5c99
13 changed files with 196 additions and 158 deletions

View file

@ -30,7 +30,10 @@ class GP(Model):
super(GP, self).__init__(name)
assert X.ndim == 2
self.X = ObservableArray(X)
if isinstance(X, ObservableArray):
self.X = self.X = X
else: self.X = ObservableArray(X)
self.num_data, self.input_dim = self.X.shape
assert Y.ndim == 2

View file

@ -28,7 +28,9 @@ class ObservableArray(np.ndarray, Observable):
"""
__array_priority__ = -1 # Never give back ObservableArray
def __new__(cls, input_array):
obj = np.atleast_1d(input_array).view(cls)
if not isinstance(input_array, ObservableArray):
obj = np.atleast_1d(input_array).view(cls)
else: obj = input_array
cls.__name__ = "ObservableArray\n "
return obj

View file

@ -3,21 +3,54 @@ Created on 6 Nov 2013
@author: maxz
'''
import numpy as np
from parameterized import Parameterized
from param import Param
from transformations import Logexp
class Normal(Parameterized):
class VariationalPrior(object):
def KL_divergence(self, variational_posterior):
raise NotImplementedError, "override this for variational inference of latent space"
def update_gradients_KL(self, variational_posterior):
"""
updates the gradients for mean and variance **in place**
"""
raise NotImplementedError, "override this for variational inference of latent space"
class NormalPrior(VariationalPrior):
def KL_divergence(self, variational_posterior):
var_mean = np.square(variational_posterior.mean).sum()
var_S = (variational_posterior.variance - np.log(variational_posterior.variance)).sum()
return 0.5 * (var_mean + var_S) - 0.5 * variational_posterior.input_dim * variational_posterior.num_data
def update_gradients_KL(self, variational_posterior):
# dL:
variational_posterior.mean.gradient -= variational_posterior.mean
variational_posterior.variance.gradient -= (1. - (1. / (variational_posterior.variance))) * 0.5
class VariationalPosterior(Parameterized):
def __init__(self, means=None, variances=None, name=None, **kw):
super(VariationalPosterior, self).__init__(name=name, **kw)
self.mean = Param("mean", means)
self.variance = Param("variance", variances, Logexp())
self.add_parameters(self.mean, self.variance)
self.num_data, self.input_dim = self.mean.shape
if self.has_uncertain_inputs():
assert self.variance.shape == self.mean.shape, "need one variance per sample and dimenion"
def has_uncertain_inputs(self):
return not self.variance is None
class NormalPosterior(VariationalPosterior):
'''
Normal distribution for variational approximations.
NormalPosterior distribution for variational approximations.
holds the means and variances for a factorizing multivariate normal distribution
'''
def __init__(self, means, variances, name='latent space'):
Parameterized.__init__(self, name=name)
self.mean = Param("mean", means)
self.variance = Param('variance', variances, Logexp())
self.add_parameters(self.mean, self.variance)
def plot(self, *args):
"""
@ -30,8 +63,7 @@ class Normal(Parameterized):
from ...plotting.matplot_dep import variational_plots
return variational_plots.plot(self,*args)
class SpikeAndSlab(Parameterized):
class SpikeAndSlabPosterior(VariationalPosterior):
'''
The SpikeAndSlab distribution for variational approximations.
'''
@ -39,11 +71,9 @@ class SpikeAndSlab(Parameterized):
"""
binary_prob : the probability of the distribution on the slab part.
"""
Parameterized.__init__(self, name=name)
self.mean = Param("mean", means)
self.variance = Param('variance', variances, Logexp())
super(SpikeAndSlabPosterior, self).__init__(means, variances, name)
self.gamma = Param("binary_prob",binary_prob,)
self.add_parameters(self.mean, self.variance, self.gamma)
self.add_parameter(self.gamma)
def plot(self, *args):
"""

View file

@ -5,8 +5,9 @@ import numpy as np
from ..util.linalg import mdot
from gp import GP
from parameterization.param import Param
from GPy.inference.latent_function_inference import var_dtc
from ..inference.latent_function_inference import var_dtc
from .. import likelihoods
from parameterization.variational import NormalPosterior
class SparseGP(GP):
"""
@ -45,16 +46,14 @@ class SparseGP(GP):
self.Z = Param('inducing inputs', Z)
self.num_inducing = Z.shape[0]
self.X_variance = X_variance
if self.has_uncertain_inputs():
assert X_variance.shape == X.shape
self.q = NormalPosterior(X, X_variance)
GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name)
GP.__init__(self, self.q.mean, Y, kernel, likelihood, inference_method=inference_method, name=name)
self.add_parameter(self.Z, index=0)
self.parameters_changed()
def has_uncertain_inputs(self):
return not (self.X_variance is None)
return self.q.has_uncertain_inputs()
def parameters_changed(self):
if self.has_uncertain_inputs():
@ -78,10 +77,11 @@ class SparseGP(GP):
mu = np.dot(Kx.T, self.posterior.woodbury_vector)
if full_cov:
Kxx = self.kern.K(Xnew)
var = Kxx - mdot(Kx.T, self.posterior.woodbury_inv, Kx)
#var = Kxx - mdot(Kx.T, self.posterior.woodbury_inv, Kx)
var = Kxx - np.tensordot(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx).T, Kx, [1,0]).swapaxes(1,2)
else:
Kxx = self.kern.Kdiag(Xnew)
var = Kxx - np.sum(Kx * np.dot(self.posterior.woodbury_inv, Kx), 0)
var = (Kxx - np.sum(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx) * Kx[None,:,:], 1)).T
else:
Kx = self.kern.psi1(self.Z, Xnew, X_variance_new)
mu = np.dot(Kx, self.Cpsi1V)
@ -91,7 +91,7 @@ class SparseGP(GP):
Kxx = self.kern.psi0(self.Z, Xnew, X_variance_new)
psi2 = self.kern.psi2(self.Z, Xnew, X_variance_new)
var = Kxx - np.sum(np.sum(psi2 * Kmmi_LmiBLmi[None, :, :], 1), 1)
return mu, var[:,None]
return mu, var
def _getstate(self):

View file

@ -89,7 +89,7 @@ def sparse_gplvm_oil(optimize=True, verbose=0, plot=True, N=100, Q=6, num_induci
Y = Y - Y.mean(0)
Y /= Y.std(0)
# Create the model
kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.bias(Q)
kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.Bias(Q)
m = GPy.models.SparseGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing)
m.data_labels = data['Y'][:N].argmax(axis=1)
@ -139,7 +139,7 @@ def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=15, Q=4
(1 - var))) + .001
Z = _np.random.permutation(X)[:num_inducing]
kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.bias(Q, _np.exp(-2)) + GPy.kern.white(Q, _np.exp(-2))
kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.Bias(Q, _np.exp(-2)) + GPy.kern.White(Q, _np.exp(-2))
m = BayesianGPLVM(Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel)
m.data_colors = c
@ -159,28 +159,26 @@ def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=15, Q=4
def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k):
import GPy
from GPy.likelihoods import Gaussian
from matplotlib import pyplot as plt
_np.random.seed(0)
data = GPy.util.datasets.oil()
kernel = GPy.kern.RBF_inv(Q, 1., [.1] * Q, ARD=True) + GPy.kern.bias(Q, _np.exp(-2))
kernel = GPy.kern.RBF(Q, 1., [.1] * Q, ARD=True)# + GPy.kern.Bias(Q, _np.exp(-2))
Y = data['X'][:N]
Yn = Gaussian(Y, normalize=True)
m = GPy.models.BayesianGPLVM(Yn, Q, kernel=kernel, num_inducing=num_inducing, **k)
m = GPy.models.BayesianGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing, **k)
m.data_labels = data['Y'][:N].argmax(axis=1)
m['noise'] = Yn.Y.var() / 100.
m['.*noise.var'] = Y.var() / 100.
if optimize:
m.optimize('scg', messages=verbose, max_iters=max_iters, gtol=.05)
if plot:
y = m.likelihood.Y[0, :]
y = m.Y[0, :]
fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
m.plot_latent(ax=latent_axes)
data_show = GPy.util.visualize.vector_show(y)
lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X[0, :], # @UnusedVariable
data_show = GPy.plotting.matplot_dep.visualize.vector_show(y)
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X[0, :], # @UnusedVariable
m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
raw_input('Press enter to finish')
plt.close(fig)
@ -190,8 +188,8 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False):
_np.random.seed(1234)
x = _np.linspace(0, 4 * _np.pi, N)[:, None]
s1 = _np.vectorize(lambda x: -_np.sin(x))
s2 = _np.vectorize(lambda x: _np.cos(x))
s1 = _np.vectorize(lambda x: -_np.sin(_np.exp(x)))
s2 = _np.vectorize(lambda x: _np.cos(x)**2)
s3 = _np.vectorize(lambda x:-_np.exp(-_np.cos(2 * x)))
sS = _np.vectorize(lambda x: x*_np.sin(x))
@ -328,7 +326,7 @@ def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw):
_, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim)
likelihood_list = [Gaussian(x, normalize=True) for x in Ylist]
k = kern.linear(Q, ARD=True) + kern.bias(Q, _np.exp(-2)) + kern.white(Q, _np.exp(-2))
k = kern.Linear(Q, ARD=True) + kern.Bias(Q, _np.exp(-2)) + kern.White(Q, _np.exp(-2))
m = MRD(likelihood_list, input_dim=Q, num_inducing=num_inducing, kernels=k, initx="", initz='permute', **kw)
m.ensure_default_constraints()
@ -355,15 +353,15 @@ def brendan_faces(optimize=True, verbose=True, plot=True):
m = GPy.models.GPLVM(Yn, Q)
# optimize
m.constrain('rbf|noise|white', GPy.core.transformations.logexp_clipped())
m.constrain('rbf|noise|white', GPy.transformations.LogexpClipped())
if optimize: m.optimize('scg', messages=verbose, max_iters=1000)
if plot:
ax = m.plot_latent(which_indices=(0, 1))
y = m.likelihood.Y[0, :]
data_show = GPy.util.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, order='F', invert=False, scale=False)
GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, order='F', invert=False, scale=False)
GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
raw_input('Press enter to finish')
return m
@ -382,8 +380,8 @@ def olivetti_faces(optimize=True, verbose=True, plot=True):
if plot:
ax = m.plot_latent(which_indices=(0, 1))
y = m.likelihood.Y[0, :]
data_show = GPy.util.visualize.image_show(y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False)
GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False)
GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
raw_input('Press enter to finish')
return m
@ -398,8 +396,8 @@ def stick_play(range=None, frame_rate=15, optimize=False, verbose=True, plot=Tru
Y = data['Y'][range[0]:range[1], :].copy()
if plot:
y = Y[0, :]
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.util.visualize.data_play(Y, data_show, frame_rate)
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.plotting.matplot_dep.visualize.data_play(Y, data_show, frame_rate)
return Y
def stick(kernel=None, optimize=True, verbose=True, plot=True):
@ -410,12 +408,12 @@ def stick(kernel=None, optimize=True, verbose=True, plot=True):
# optimize
m = GPy.models.GPLVM(data['Y'], 2, kernel=kernel)
if optimize: m.optimize(messages=verbose, max_f_eval=10000)
if plot and GPy.util.visualize.visual_available:
if plot and GPy.plotting.matplot_dep.visualize.visual_available:
plt.clf
ax = m.plot_latent()
y = m.likelihood.Y[0, :]
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
raw_input('Press enter to finish')
return m
@ -429,12 +427,12 @@ def bcgplvm_linear_stick(kernel=None, optimize=True, verbose=True, plot=True):
mapping = GPy.mappings.Linear(data['Y'].shape[1], 2)
m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping)
if optimize: m.optimize(messages=verbose, max_f_eval=10000)
if plot and GPy.util.visualize.visual_available:
if plot and GPy.plotting.matplot_dep.visualize.visual_available:
plt.clf
ax = m.plot_latent()
y = m.likelihood.Y[0, :]
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
raw_input('Press enter to finish')
return m
@ -449,12 +447,12 @@ def bcgplvm_stick(kernel=None, optimize=True, verbose=True, plot=True):
mapping = GPy.mappings.Kernel(X=data['Y'], output_dim=2, kernel=back_kernel)
m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping)
if optimize: m.optimize(messages=verbose, max_f_eval=10000)
if plot and GPy.util.visualize.visual_available:
if plot and GPy.plotting.matplot_dep.visualize.visual_available:
plt.clf
ax = m.plot_latent()
y = m.likelihood.Y[0, :]
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
raw_input('Press enter to finish')
return m
@ -480,7 +478,7 @@ def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True):
data = GPy.util.datasets.osu_run1()
Q = 6
kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.bias(Q, _np.exp(-2)) + GPy.kern.white(Q, _np.exp(-2))
kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.Bias(Q, _np.exp(-2)) + GPy.kern.White(Q, _np.exp(-2))
m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20, kernel=kernel)
# optimize
m.ensure_default_constraints()
@ -491,8 +489,8 @@ def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True):
plt.sca(latent_axes)
m.plot_latent()
y = m.likelihood.Y[0, :].copy()
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.util.visualize.lvm_dimselect(m.X[0, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect'])
GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X[0, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
raw_input('Press enter to finish')
return m
@ -511,8 +509,8 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose
if plot:
ax = m.plot_latent()
y = m.likelihood.Y[0, :]
data_show = GPy.util.visualize.skeleton_show(y[None, :], data['skel'])
lvm_visualizer = GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
data_show = GPy.plotting.matplot_dep.visualize.skeleton_show(y[None, :], data['skel'])
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
raw_input('Press enter to finish')
lvm_visualizer.close()

View file

@ -2,7 +2,7 @@
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import numpy as np
from ...util.linalg import pdinv, dpotrs, tdot, dtrtrs, dpotri, symmetrify, jitchol, dtrtri
from ...util.linalg import pdinv, dpotrs, dpotri, symmetrify, jitchol
class Posterior(object):
"""
@ -83,14 +83,15 @@ class Posterior(object):
#LiK, _ = dtrtrs(self.woodbury_chol, self._K, lower=1)
self._covariance = np.tensordot(np.dot(np.atleast_3d(self.woodbury_inv).T, self._K), self._K, [1,0]).T
#self._covariance = self._K - self._K.dot(self.woodbury_inv).dot(self._K)
return self._covariance
return self._covariance.squeeze()
@property
def precision(self):
if self._precision is None:
self._precision = np.zeros(np.atleast_3d(self.covariance).shape) # if one covariance per dimension
for p in xrange(self.covariance.shape[-1]):
self._precision[:,:,p] = pdinv(self.covariance[:,:,p])[0]
cov = np.atleast_3d(self.covariance)
self._precision = np.zeros(cov.shape) # if one covariance per dimension
for p in xrange(cov.shape[-1]):
self._precision[:,:,p] = pdinv(cov[:,:,p])[0]
return self._precision
@property
@ -98,7 +99,10 @@ class Posterior(object):
if self._woodbury_chol is None:
#compute woodbury chol from
if self._woodbury_inv is not None:
_, _, self._woodbury_chol, _ = pdinv(self._woodbury_inv)
winv = np.atleast_3d(self._woodbury_inv)
self._woodbury_chol = np.zeros(winv.shape)
for p in xrange(winv.shape[-1]):
self._woodbury_chol[:,:,p] = pdinv(winv[:,:,p])[2]
#Li = jitchol(self._woodbury_inv)
#self._woodbury_chol, _ = dtrtri(Li)
#W, _, _, _, = pdinv(self._woodbury_inv)
@ -132,7 +136,7 @@ class Posterior(object):
@property
def K_chol(self):
if self._K_chol is None:
self._K_chol = dportf(self._K)
self._K_chol = jitchol(self._K)
return self._K_chol

View file

@ -127,11 +127,12 @@ from GPy.core.model import Model
class Kern_check_model(Model):
"""This is a dummy model class used as a base class for checking that the gradients of a given kernel are implemented correctly. It enables checkgradient() to be called independently on a kernel."""
def __init__(self, kernel=None, dL_dK=None, X=None, X2=None):
from GPy.kern import RBF
Model.__init__(self, 'kernel_test_model')
num_samples = 20
num_samples2 = 10
if kernel==None:
kernel = GPy.kern.rbf(1)
kernel = RBF(1)
if X==None:
X = np.random.randn(num_samples, kernel.input_dim)
if dL_dK==None:

View file

@ -106,51 +106,52 @@ class Linear(Kern):
# variational #
#---------------------------------------#
def psi0(self, Z, mu, S):
return np.sum(self.variances * self._mu2S(mu, S), 1)
def psi0(self, Z, posterior_variational):
return np.sum(self.variances * self._mu2S(posterior_variational), 1)
def psi1(self, Z, mu, S):
return self.K(mu, Z) #the variance, it does nothing
def psi1(self, Z, posterior_variational):
return self.K(posterior_variational.mean, Z) #the variance, it does nothing
def psi2(self, Z, mu, S):
def psi2(self, Z, posterior_variational):
ZA = Z * self.variances
ZAinner = self._ZAinner(mu, S, Z)
ZAinner = self._ZAinner(posterior_variational, Z)
return np.dot(ZAinner, ZA.T)
def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z):
def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, posterior_variational, Z):
mu, S = posterior_variational.mean, posterior_variational.variance
# psi0:
tmp = dL_dpsi0[:, None] * self._mu2S(mu, S)
tmp = dL_dpsi0[:, None] * self._mu2S(posterior_variational)
if self.ARD: grad = tmp.sum(0)
else: grad = np.atleast_1d(tmp.sum())
#psi1
self.update_gradients_full(dL_dpsi1, mu, Z)
grad += self.variances.gradient
#psi2
tmp = dL_dpsi2[:, :, :, None] * (self._ZAinner(mu, S, Z)[:, :, None, :] * (2. * Z)[None, None, :, :])
tmp = dL_dpsi2[:, :, :, None] * (self._ZAinner(posterior_variational, Z)[:, :, None, :] * (2. * Z)[None, None, :, :])
if self.ARD: grad += tmp.sum(0).sum(0).sum(0)
else: grad += tmp.sum()
#from Kmm
self.update_gradients_full(dL_dKmm, Z, None)
self.variances.gradient += grad
def gradients_Z_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z):
def gradients_Z_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, posterior_variational, Z):
# Kmm
grad = self.gradients_X(dL_dKmm, Z, None)
#psi1
grad += self.gradients_X(dL_dpsi1.T, Z, mu)
grad += self.gradients_X(dL_dpsi1.T, Z, posterior_variational.mean)
#psi2
self._weave_dpsi2_dZ(dL_dpsi2, Z, mu, S, grad)
self._weave_dpsi2_dZ(dL_dpsi2, Z, posterior_variational, grad)
return grad
def gradients_muS_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z):
grad_mu, grad_S = np.zeros(mu.shape), np.zeros(mu.shape)
def gradients_q_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, posterior_variational, Z):
grad_mu, grad_S = np.zeros(posterior_variational.mean.shape), np.zeros(posterior_variational.mean.shape)
# psi0
grad_mu += dL_dpsi0[:, None] * (2.0 * mu * self.variances)
grad_mu += dL_dpsi0[:, None] * (2.0 * posterior_variational.mean * self.variances)
grad_S += dL_dpsi0[:, None] * self.variances
# psi1
grad_mu += (dL_dpsi1[:, :, None] * (Z * self.variances)).sum(1)
# psi2
self._weave_dpsi2_dmuS(dL_dpsi2, Z, mu, S, grad_mu, grad_S)
self._weave_dpsi2_dmuS(dL_dpsi2, Z, posterior_variational, grad_mu, grad_S)
return grad_mu, grad_S
@ -159,7 +160,7 @@ class Linear(Kern):
#--------------------------------------------------#
def _weave_dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
def _weave_dpsi2_dmuS(self, dL_dpsi2, Z, pv, target_mu, target_S):
# Think N,num_inducing,num_inducing,input_dim
ZA = Z * self.variances
AZZA = ZA.T[:, None, :, None] * ZA[None, :, None, :]
@ -203,14 +204,15 @@ class Linear(Kern):
'extra_compile_args': ['-fopenmp -O3'], #-march=native'],
'extra_link_args' : ['-lgomp']}
mu = pv.mean
N,num_inducing,input_dim,mu = mu.shape[0],Z.shape[0],mu.shape[1],param_to_array(mu)
weave.inline(code, support_code=support_code, libraries=['gomp'],
arg_names=['N','num_inducing','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
type_converters=weave.converters.blitz,**weave_options)
def _weave_dpsi2_dZ(self, dL_dpsi2, Z, mu, S, target):
AZA = self.variances*self._ZAinner(mu, S, Z)
def _weave_dpsi2_dZ(self, dL_dpsi2, Z, pv, target):
AZA = self.variances*self._ZAinner(pv, Z)
code="""
int n,m,mm,q;
#pragma omp parallel for private(n,mm,q)
@ -232,21 +234,21 @@ class Linear(Kern):
'extra_compile_args': ['-fopenmp -O3'], #-march=native'],
'extra_link_args' : ['-lgomp']}
N,num_inducing,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
mu = param_to_array(mu)
N,num_inducing,input_dim = pv.mean.shape[0],Z.shape[0],pv.mean.shape[1]
mu = param_to_array(pv.mean)
weave.inline(code, support_code=support_code, libraries=['gomp'],
arg_names=['N','num_inducing','input_dim','AZA','target','dL_dpsi2'],
type_converters=weave.converters.blitz,**weave_options)
def _mu2S(self, mu, S):
return np.square(mu) + S
def _mu2S(self, pv):
return np.square(pv.mean) + pv.variance
def _ZAinner(self, mu, S, Z):
def _ZAinner(self, pv, Z):
ZA = Z*self.variances
inner = (mu[:, None, :] * mu[:, :, None])
diag_indices = np.diag_indices(mu.shape[1], 2)
inner[:, diag_indices[0], diag_indices[1]] += S
inner = (pv.mean[:, None, :] * pv.mean[:, :, None])
diag_indices = np.diag_indices(pv.mean.shape[1], 2)
inner[:, diag_indices[0], diag_indices[1]] += pv.variance
return np.dot(ZA, inner).swapaxes(0, 1) # NOTE: self.ZAinner \in [num_inducing x N x input_dim]!

View file

@ -182,7 +182,7 @@ class RBF(Kern):
return grad
def gradients_q_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, posterior_variational):
def update_gradients_q_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, posterior_variational):
mu = posterior_variational.mean
S = posterior_variational.variance
self._psi_computations(Z, mu, S)
@ -195,7 +195,8 @@ class RBF(Kern):
grad_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * self._psi2_mudist).sum(1).sum(1)
grad_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*self._psi2_mudist_sq - 1)).sum(1).sum(1)
return grad_mu, grad_S
posterior_variational.mean.gradient = grad_mu
posterior_variational.variance.gradient = grad_S
def gradients_X(self, dL_dK, X, X2=None):
#if self._X is None or X.base is not self._X.base or X2 is not None:

View file

@ -8,7 +8,7 @@ from ..core import SparseGP
from ..likelihoods import Gaussian
from ..inference.optimization import SCG
from ..util import linalg
from ..core.parameterization.variational import Normal
from ..core.parameterization.variational import NormalPosterior, NormalPrior
class BayesianGPLVM(SparseGP, GPLVM):
"""
@ -29,7 +29,7 @@ class BayesianGPLVM(SparseGP, GPLVM):
self.init = init
if X_variance is None:
X_variance = np.clip((np.ones_like(X) * 0.5) + .01 * np.random.randn(*X.shape), 0.001, 1)
X_variance = np.random.uniform(0,.1,X.shape)
if Z is None:
Z = np.random.permutation(X.copy())[:num_inducing]
@ -40,7 +40,9 @@ class BayesianGPLVM(SparseGP, GPLVM):
if likelihood is None:
likelihood = Gaussian()
self.q = Normal(X, X_variance)
self.q = NormalPosterior(X, X_variance)
self.variational_prior = NormalPrior()
SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method, X_variance, name, **kwargs)
self.add_parameter(self.q, index=0)
#self.ensure_default_constraints()
@ -57,24 +59,15 @@ class BayesianGPLVM(SparseGP, GPLVM):
self.init = state.pop()
SparseGP._setstate(self, state)
def KL_divergence(self):
var_mean = np.square(self.X).sum()
var_S = np.sum(self.X_variance - np.log(self.X_variance))
return 0.5 * (var_mean + var_S) - 0.5 * self.input_dim * self.num_data
def parameters_changed(self):
super(BayesianGPLVM, self).parameters_changed()
self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.q)
self._log_marginal_likelihood -= self.KL_divergence()
dL_dmu, dL_dS = self.kern.gradients_q_variational(posterior_variational=self.q, Z=self.Z, **self.grad_dict)
self.kern.update_gradients_q_variational(posterior_variational=self.q, Z=self.Z, **self.grad_dict)
# dL:
self.q.mean.gradient = dL_dmu
self.q.variance.gradient = dL_dS
# update for the KL divergence
self.variational_prior.update_gradients_KL(self.q)
# dKL:
self.q.mean.gradient -= self.X
self.q.variance.gradient -= (1. - (1. / (self.X_variance))) * 0.5
def plot_latent(self, plot_inducing=True, *args, **kwargs):
"""
@ -147,6 +140,7 @@ class BayesianGPLVM(SparseGP, GPLVM):
"""
See GPy.plotting.matplot_dep.dim_reduction_plots.plot_steepest_gradient_map
"""
import sys
assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
from ..plotting.matplot_dep import dim_reduction_plots

View file

@ -60,6 +60,8 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
X, Y = param_to_array(model.X, model.Y)
if hasattr(model, 'has_uncertain_inputs') and model.has_uncertain_inputs(): X_variance = model.X_variance
if hasattr(model, 'Z'): Z = param_to_array(model.Z)
#work out what the inputs are for plotting (1D or 2D)
fixed_dims = np.array([i for i,v in fixed_inputs])
free_dims = np.setdiff1d(np.arange(model.input_dim),fixed_dims)
@ -96,10 +98,10 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
#add error bars for uncertain (if input uncertainty is being modelled)
if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs():
ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(),
xerr=2 * np.sqrt(X_variance[which_data_rows, free_dims].flatten()),
ecolor='k', fmt=None, elinewidth=.5, alpha=.5)
#if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs():
# ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(),
# xerr=2 * np.sqrt(X_variance[which_data_rows, free_dims].flatten()),
# ecolor='k', fmt=None, elinewidth=.5, alpha=.5)
#set the limits of the plot to some sensible values
@ -111,7 +113,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
#add inducing inputs (if a sparse model is used)
if hasattr(model,"Z"):
#Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims]
Zu = param_to_array(model.Z[:,free_dims])
Zu = Z[:,free_dims]
z_height = ax.get_ylim()[0]
ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12)
@ -135,7 +137,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
Y = Y
else:
m, _, _, _ = model.predict(Xgrid)
Y = model.data
Y = Y
for d in which_data_ycols:
m_d = m[:,d].reshape(resolution, resolution).T
ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)
@ -151,7 +153,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
#add inducing inputs (if a sparse model is used)
if hasattr(model,"Z"):
#Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims]
Zu = model.Z[:,free_dims]
Zu = Z[:,free_dims]
ax.plot(Zu[:,free_dims[0]], Zu[:,free_dims[1]], 'wo')
else:

View file

@ -30,6 +30,11 @@ class Test(unittest.TestCase):
self.assertListEqual(self.param_index[two].tolist(), [0,3])
self.assertListEqual(self.param_index[one].tolist(), [1])
def test_shift_right(self):
self.param_index.shift_right(5, 2)
self.assertListEqual(self.param_index[three].tolist(), [2,4,9])
self.assertListEqual(self.param_index[two].tolist(), [0,7])
self.assertListEqual(self.param_index[one].tolist(), [3])
def test_index_view(self):
#=======================================================================

View file

@ -10,6 +10,7 @@ from functools import partial
#np.random.seed(300)
#np.random.seed(7)
np.seterr(divide='raise')
def dparam_partial(inst_func, *args):
"""
If we have a instance method that needs to be called but that doesn't
@ -149,9 +150,9 @@ class TestNoiseModels(object):
noise_models = {"Student_t_default": {
"model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": {
"names": ["t_noise"],
"names": [".*t_noise"],
"vals": [self.var],
"constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
"constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
#"constraints": [("t_noise", constrain_positive), ("deg_free", partial(constrain_fixed, value=5))]
},
"laplace": True
@ -159,66 +160,66 @@ class TestNoiseModels(object):
"Student_t_1_var": {
"model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": {
"names": ["t_noise"],
"names": [".*t_noise"],
"vals": [1.0],
"constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
"constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
},
"laplace": True
},
"Student_t_small_deg_free": {
"model": GPy.likelihoods.StudentT(deg_free=1.5, sigma2=self.var),
"grad_params": {
"names": ["t_noise"],
"names": [".*t_noise"],
"vals": [self.var],
"constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
"constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
},
"laplace": True
},
"Student_t_small_var": {
"model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": {
"names": ["t_noise"],
"vals": [0.0001],
"constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
"names": [".*t_noise"],
"vals": [0.001],
"constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
},
"laplace": True
},
"Student_t_large_var": {
"model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
"grad_params": {
"names": ["t_noise"],
"names": [".*t_noise"],
"vals": [10.0],
"constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
"constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
},
"laplace": True
},
"Student_t_approx_gauss": {
"model": GPy.likelihoods.StudentT(deg_free=1000, sigma2=self.var),
"grad_params": {
"names": ["t_noise"],
"names": [".*t_noise"],
"vals": [self.var],
"constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
"constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
},
"laplace": True
},
"Student_t_log": {
"model": GPy.likelihoods.StudentT(gp_link=link_functions.Log(), deg_free=5, sigma2=self.var),
"grad_params": {
"names": ["t_noise"],
"names": [".*t_noise"],
"vals": [self.var],
"constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
"constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
},
"laplace": True
},
"Gaussian_default": {
"model": GPy.likelihoods.Gaussian(variance=self.var),
"grad_params": {
"names": ["variance"],
"names": [".*variance"],
"vals": [self.var],
"constraints": [("variance", constrain_positive)]
"constraints": [(".*variance", constrain_positive)]
},
"laplace": True,
"ep": True
"ep": False # FIXME: Should be True when we have it working again
},
#"Gaussian_log": {
#"model": GPy.likelihoods.gaussian(gp_link=link_functions.Log(), variance=self.var, D=self.D, N=self.N),
@ -252,7 +253,7 @@ class TestNoiseModels(object):
"link_f_constraints": [partial(constrain_bounded, lower=0, upper=1)],
"laplace": True,
"Y": self.binary_Y,
"ep": True
"ep": False # FIXME: Should be True when we have it working again
},
#"Exponential_default": {
#"model": GPy.likelihoods.exponential(),
@ -515,10 +516,10 @@ class TestNoiseModels(object):
#Normalize
Y = Y/Y.max()
white_var = 1e-6
kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
kernel = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
laplace_likelihood = GPy.inference.latent_function_inference.Laplace()
m = GPy.core.GP(X.copy(), Y.copy(), kernel, likelihood=model, inference_method=laplace_likelihood)
m['white'].constrain_fixed(white_var)
m['.*white'].constrain_fixed(white_var)
#Set constraints
for constrain_param, constraint in constraints:
@ -541,9 +542,6 @@ class TestNoiseModels(object):
#NOTE this test appears to be stochastic for some likelihoods (student t?)
# appears to all be working in test mode right now...
if not m.checkgrad():
import ipdb; ipdb.set_trace() # XXX BREAKPOINT
assert m.checkgrad(step=step)
###########
@ -555,10 +553,10 @@ class TestNoiseModels(object):
#Normalize
Y = Y/Y.max()
white_var = 1e-6
kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
kernel = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
ep_inf = GPy.inference.latent_function_inference.EP()
m = GPy.core.GP(X.copy(), Y.copy(), kernel=kernel, likelihood=model, inference_method=ep_inf)
m['white'].constrain_fixed(white_var)
m['.*white'].constrain_fixed(white_var)
for param_num in range(len(param_names)):
name = param_names[param_num]
@ -634,26 +632,24 @@ class LaplaceTests(unittest.TestCase):
Y = Y/Y.max()
#Yc = Y.copy()
#Yc[75:80] += 1
kernel1 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
kernel1 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
#FIXME: Make sure you can copy kernels when params is fixed
#kernel2 = kernel1.copy()
kernel2 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
kernel2 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
gauss_distr1 = GPy.likelihoods.Gaussian(variance=initial_var_guess)
exact_inf = GPy.inference.latent_function_inference.ExactGaussianInference()
m1 = GPy.core.GP(X, Y.copy(), kernel=kernel1, likelihood=gauss_distr1, inference_method=exact_inf)
m1['white'].constrain_fixed(1e-6)
m1['variance'] = initial_var_guess
m1['variance'].constrain_bounded(1e-4, 10)
m1['rbf'].constrain_bounded(1e-4, 10)
m1['.*white'].constrain_fixed(1e-6)
m1['.*rbf.variance'] = initial_var_guess
m1['.*rbf.variance'].constrain_bounded(1e-4, 10)
m1.randomize()
gauss_distr2 = GPy.likelihoods.Gaussian(variance=initial_var_guess)
laplace_inf = GPy.inference.latent_function_inference.Laplace()
m2 = GPy.core.GP(X, Y.copy(), kernel=kernel2, likelihood=gauss_distr2, inference_method=laplace_inf)
m2['white'].constrain_fixed(1e-6)
m2['rbf'].constrain_bounded(1e-4, 10)
m2['variance'].constrain_bounded(1e-4, 10)
m2['.*white'].constrain_fixed(1e-6)
m2['.*rbf.variance'].constrain_bounded(1e-4, 10)
m2.randomize()
if debug: