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

Merge branch 'params' of github.com:SheffieldML/GPy into params

Browse source
This commit is contained in:
James Hensman 2014-03-21 13:57:59 +00:00
commit 3a2e15c508
25 changed files with 198 additions and 151 deletions
Download patch file Download diff file Expand all files Collapse all files

22
GPy/core/parameterization/parameter_core.py
View file

@ -541,12 +541,12 @@ class Constrainable(Nameable, Indexable):
print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name)
which.add(what, self._raveled_index())
def _remove_from_index_operations(self, which, what):
def _remove_from_index_operations(self, which, transforms):
"""
Helper preventing copy code.
Remove given what (transform prior etc) from which param index ops.
"""
if len(what) == 0:
if len(transforms) == 0:
transforms = which.properties()
removed = np.empty((0,), dtype=int)
for t in transforms:
@ -567,10 +567,10 @@ class OptimizationHandlable(Constrainable):
super(OptimizationHandlable, self).__init__(name, default_constraint=default_constraint, *a, **kw)
def transform(self):
[np.put(self._param_array_, ind, c.finv(self._param_array_[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
[np.put(self._param_array_, ind, c.finv(self._param_array_.flat[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
def untransform(self):
[np.put(self._param_array_, ind, c.f(self._param_array_[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
[np.put(self._param_array_, ind, c.f(self._param_array_.flat[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
def _get_params_transformed(self):
# transformed parameters (apply transformation rules)
@ -590,9 +590,9 @@ class OptimizationHandlable(Constrainable):
if self.has_parent() and self.constraints[__fixed__].size != 0:
fixes = np.ones(self.size).astype(bool)
fixes[self.constraints[__fixed__]] = FIXED
self._param_array_[fixes] = p
elif self._has_fixes(): self._param_array_[self._fixes_] = p
else: self._param_array_[:] = p
self._param_array_.flat[fixes] = p
elif self._has_fixes(): self._param_array_.flat[self._fixes_] = p
else: self._param_array_.flat = p
self.untransform()
self._trigger_params_changed()
@ -670,8 +670,8 @@ class OptimizationHandlable(Constrainable):
for pi in self._parameters_:
pislice = slice(pi_old_size, pi_old_size+pi.size)
self._param_array_[pislice] = pi._param_array_.ravel()#, requirements=['C', 'W']).flat
self._gradient_array_[pislice] = pi._gradient_array_.ravel()#, requirements=['C', 'W']).flat
self._param_array_[pislice] = pi._param_array_.flat#, requirements=['C', 'W']).flat
self._gradient_array_[pislice] = pi._gradient_array_.flat#, requirements=['C', 'W']).flat
pi._param_array_.data = parray[pislice].data
pi._gradient_array_.data = garray[pislice].data
@ -878,8 +878,8 @@ class Parameterizable(OptimizationHandlable):
# first connect all children
p._propagate_param_grad(self._param_array_[pslice], self._gradient_array_[pslice])
# then connect children to self
self._param_array_[pslice] = p._param_array_.ravel()#, requirements=['C', 'W']).ravel(order='C')
self._gradient_array_[pslice] = p._gradient_array_.ravel()#, requirements=['C', 'W']).ravel(order='C')
self._param_array_[pslice] = p._param_array_.flat#, requirements=['C', 'W']).ravel(order='C')
self._gradient_array_[pslice] = p._gradient_array_.flat#, requirements=['C', 'W']).ravel(order='C')
if not p._param_array_.flags['C_CONTIGUOUS']:
import ipdb;ipdb.set_trace()

8
GPy/core/sparse_gp.py
View file

@ -64,8 +64,8 @@ class SparseGP(GP):
self.kern.gradient += target
#gradients wrt Z
self.Z.gradient[:,self.kern.active_dims] = self.kern.gradients_X(dL_dKmm, self.Z)
self.Z.gradient[:,self.kern.active_dims] += self.kern.gradients_Z_expectations(
self.Z.gradient = self.kern.gradients_X(dL_dKmm, self.Z)
self.Z.gradient += self.kern.gradients_Z_expectations(
self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X)
else:
#gradients wrt kernel
@ -76,8 +76,8 @@ class SparseGP(GP):
self.kern.update_gradients_full(self.grad_dict['dL_dKmm'], self.Z, None)
self.kern.gradient += target
#gradients wrt Z
self.Z.gradient[:,self.kern.active_dims] = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z)
self.Z.gradient[:,self.kern.active_dims] += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X)
self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z)
self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X)
def _raw_predict(self, Xnew, full_cov=False):
"""

11
GPy/examples/coreg_example.py
View file

@ -23,13 +23,10 @@ K = Bias.prod(Coreg,name='X')
#K.coregion.W = 0
#print K.coregion.W
#print Bias.K(_X,_X)
#print K.K(X,X)
#pb.matshow(K.K(X,X))
Mlist = [GPy.kern.Matern32(1,lengthscale=20.,name="Mat")]
kern = GPy.util.multioutput.LCM(input_dim=1,num_outputs=2,kernels_list=Mlist,name='H')
kern.B.W = 0
@ -37,16 +34,22 @@ kern.B.kappa = 1.
#kern.B.W.fix()
#kern.B.kappa.fix()
#m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2], kernel=kern)
m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1], Y_list=[Y1], kernel=kern)
Z1 = np.array([1.5,2.5])[:,None]
m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1], Y_list=[Y1], Z_list = [Z1], kernel=kern)
#m.optimize()
m.checkgrad(verbose=1)
"""
fig = pb.figure()
ax0 = fig.add_subplot(211)
ax1 = fig.add_subplot(212)
slices = GPy.util.multioutput.get_slices([Y1,Y2])
m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],ax=ax0)
#m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],ax=ax1)
"""

5
GPy/examples/dimensionality_reduction.py
View file

@ -160,6 +160,7 @@ 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 matplotlib import pyplot as plt
from ..util.misc import param_to_array
_np.random.seed(0)
data = GPy.util.datasets.oil()
@ -173,11 +174,11 @@ def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40,
m.optimize('scg', messages=verbose, max_iters=max_iters, gtol=.05)
if plot:
y = m.Y[0, :]
y = m.Y
fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
m.plot_latent(ax=latent_axes)
data_show = GPy.plotting.matplot_dep.visualize.vector_show(y)
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X[0, :], # @UnusedVariable
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(param_to_array(m.X.mean), # @UnusedVariable
m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
raw_input('Press enter to finish')
plt.close(fig)

63
GPy/examples/regression.py
View file

@ -25,80 +25,51 @@ def olympic_marathon_men(optimize=True, plot=True):
return m
def coregionalization_toy2(optimize=True, plot=True):
def coregionalization_toy(optimize=True, plot=True):
"""
A simple demonstration of coregionalization on two sinusoidal functions.
"""
#build a design matrix with a column of integers indicating the output
X1 = np.random.rand(50, 1) * 8
X2 = np.random.rand(30, 1) * 5
index = np.vstack((np.zeros_like(X1), np.ones_like(X2)))
X = np.hstack((np.vstack((X1, X2)), index))
#build a suitable set of observed variables
Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.
Y = np.vstack((Y1, Y2))
#build the kernel
k1 = GPy.kern.RBF(1) + GPy.kern.Bias(1)
k2 = GPy.kern.Coregionalize(2,1)
k = k1**k2
m = GPy.models.GPRegression(X, Y, kernel=k)
m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2])
if optimize:
m.optimize('bfgs', max_iters=100)
if plot:
m.plot(fixed_inputs=[(1,0)])
m.plot(fixed_inputs=[(1,1)], ax=pb.gca())
slices = GPy.util.multioutput.get_slices([X1,X2])
m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0})
m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca())
return m
#FIXME: Needs recovering once likelihoods are consolidated
#def coregionalization_toy(optimize=True, plot=True):
# """
# A simple demonstration of coregionalization on two sinusoidal functions.
# """
# X1 = np.random.rand(50, 1) * 8
# X2 = np.random.rand(30, 1) * 5
# X = np.vstack((X1, X2))
# Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
# Y2 = -np.sin(X2) + np.random.randn(*X2.shape) * 0.05
# Y = np.vstack((Y1, Y2))
#
# k1 = GPy.kern.RBF(1)
# m = GPy.models.GPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1])
# m.constrain_fixed('.*rbf_var', 1.)
# m.optimize(max_iters=100)
#
# fig, axes = pb.subplots(2,1)
# m.plot(fixed_inputs=[(1,0)],ax=axes[0])
# m.plot(fixed_inputs=[(1,1)],ax=axes[1])
# axes[0].set_title('Output 0')
# axes[1].set_title('Output 1')
# return m
def coregionalization_sparse(optimize=True, plot=True):
"""
A simple demonstration of coregionalization on two sinusoidal functions using sparse approximations.
"""
#fetch the data from the non sparse examples
m = coregionalization_toy2(optimize=False, plot=False)
X, Y = m.X, m.Y
#build a design matrix with a column of integers indicating the output
X1 = np.random.rand(50, 1) * 8
X2 = np.random.rand(30, 1) * 5
k = GPy.kern.RBF(1)**GPy.kern.Coregionalize(2)
#build a suitable set of observed variables
Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.
#construct a model
m = GPy.models.SparseGPRegression(X,Y, num_inducing=25, kernel=k)
m.Z[:,1].fix() # don't optimize the inducing input indexes
m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2])
if optimize:
m.optimize('bfgs', max_iters=100, messages=1)
m.optimize('bfgs', max_iters=100)
if plot:
m.plot(fixed_inputs=[(1,0)])
m.plot(fixed_inputs=[(1,1)], ax=pb.gca())
slices = GPy.util.multioutput.get_slices([X1,X2])
m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0})
m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca())
pb.ylim(-3,)
return m

4
GPy/inference/latent_function_inference/expectation_propagation.py
View file

@ -11,9 +11,9 @@ class EP(object):
:param epsilon: Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float)
:type epsilon: float
:param eta: Power EP thing TODO: Ricardo: what, exactly?
:param eta: parameter for fractional EP updates.
:type eta: float64
:param delta: Power EP thing TODO: Ricardo: what, exactly?
:param delta: damping EP updates factor.
:type delta: float64
"""
self.epsilon, self.eta, self.delta = epsilon, eta, delta

3
GPy/inference/latent_function_inference/var_dtc.py
View file

@ -176,7 +176,6 @@ class VarDTC(object):
#construct a posterior object
post = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector, K=Kmm, mean=None, cov=None, K_chol=Lm)
return post, log_marginal, grad_dict
class VarDTCMissingData(object):
@ -365,7 +364,7 @@ class VarDTCMissingData(object):
return post, log_marginal, grad_dict
def _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, VVT_factor, Cpsi1Vf, DBi_plus_BiPBi, psi1, het_noise, uncertain_inputs):
dL_dpsi0 = -0.5 * output_dim * (beta * np.ones([num_data, 1])).flatten()
dL_dpsi0 = -0.5 * output_dim * (beta[:,None] * np.ones([num_data, 1])).flatten()
dL_dpsi1 = np.dot(VVT_factor, Cpsi1Vf.T)
dL_dpsi2_beta = 0.5 * backsub_both_sides(Lm, output_dim * np.eye(num_inducing) - DBi_plus_BiPBi)
if het_noise:

10
GPy/kern/_src/ODE_UY.py
View file

@ -139,11 +139,11 @@ class ODE_UY(Kern):
dVdly = 0.5/np.sqrt(ly)*np.sqrt(2*Vy)
dVdVy = 0.5/np.sqrt(Vy)*np.sqrt(2*ly)
rd=rdist.shape[0]
dktheta1 = np.zeros([rd,rd])
dktheta2 = np.zeros([rd,rd])
dkUdvar = np.zeros([rd,rd])
dkYdvar = np.zeros([rd,rd])
rd=rdist.shape
dktheta1 = np.zeros(rd)
dktheta2 = np.zeros(rd)
dkUdvar = np.zeros(rd)
dkYdvar = np.zeros(rd)
# dk dtheta for UU
UUdtheta1 = lambda dist: np.exp(-lu* dist)*dist + (-dist)*np.exp(-lu* dist)*(1+lu*dist)

6
GPy/kern/_src/add.py
View file

@ -23,7 +23,7 @@ class Add(CombinationKernel):
If a list of parts (of this kernel!) `which_parts` is given, only
the parts of the list are taken to compute the covariance.
"""
assert X.shape[1] == self.input_dim
assert X.shape[1] > max(np.r_[self.active_dims])
if which_parts is None:
which_parts = self.parts
elif not isinstance(which_parts, (list, tuple)):
@ -33,7 +33,7 @@ class Add(CombinationKernel):
@Cache_this(limit=2, force_kwargs=['which_parts'])
def Kdiag(self, X, which_parts=None):
assert X.shape[1] == self.input_dim
assert X.shape[1] > max(np.r_[self.active_dims])
if which_parts is None:
which_parts = self.parts
elif not isinstance(which_parts, (list, tuple)):
@ -172,7 +172,7 @@ class Add(CombinationKernel):
def add(self, other, name='sum'):
if isinstance(other, Add):
other_params = other._parameters_.copy()
other_params = other._parameters_[:]
for p in other_params:
other.remove_parameter(p)
self.add_parameters(*other_params)

3
GPy/kern/_src/linear.py
View file

@ -312,5 +312,4 @@ class Linear(Kern):
return np.dot(ZA, inner).swapaxes(0, 1) # NOTE: self.ZAinner \in [num_inducing x num_data x input_dim]!
def input_sensitivity(self):
if self.ARD: return self.variances
else: return self.variances.repeat(self.input_dim)
return np.ones(self.input_dim) * self.variances

2
GPy/likelihoods/gaussian.py
View file

@ -92,7 +92,7 @@ class Gaussian(Likelihood):
def predictive_variance(self, mu, sigma, predictive_mean=None):
return self.variance + sigma**2
def predictive_quantiles(self, mu, var, quantiles, Y_metadata):
def predictive_quantiles(self, mu, var, quantiles, Y_metadata=None):
return [stats.norm.ppf(q/100.)*np.sqrt(var) + mu for q in quantiles]
def pdf_link(self, link_f, y, Y_metadata=None):

3
GPy/likelihoods/likelihood.py
View file

@ -411,10 +411,7 @@ class Likelihood(Parameterized):
#compute the quantiles by sampling!!!
N_samp = 1000
s = np.random.randn(mu.shape[0], N_samp)*np.sqrt(var) + mu
#ss_f = s.flatten()
#ss_y = self.samples(ss_f, Y_metadata)
ss_y = self.samples(s, Y_metadata)
#ss_y = ss_y.reshape(mu.shape[0], N_samp)
return [np.percentile(ss_y ,q, axis=1)[:,None] for q in quantiles]

38
GPy/likelihoods/mixed_noise.py
View file

@ -11,7 +11,7 @@ import itertools
class MixedNoise(Likelihood):
def __init__(self, likelihoods_list, name='mixed_noise'):
#NOTE at the moment this likelihood only works for using a list of gaussians
super(Likelihood, self).__init__(name=name)
self.add_parameters(*likelihoods_list)
@ -24,11 +24,11 @@ class MixedNoise(Likelihood):
variance = np.zeros(ind.size)
for lik, j in zip(self.likelihoods_list, range(len(self.likelihoods_list))):
variance[ind==j] = lik.variance
return variance[:,None]
return variance
def betaY(self,Y,Y_metadata):
#TODO not here.
return Y/self.gaussian_variance(Y_metadata=Y_metadata)
return Y/self.gaussian_variance(Y_metadata=Y_metadata)[:,None]
def update_gradients(self, gradients):
self.gradient = gradients
@ -39,24 +39,27 @@ class MixedNoise(Likelihood):
return np.array([dL_dKdiag[ind==i].sum() for i in range(len(self.likelihoods_list))])
def predictive_values(self, mu, var, full_cov=False, Y_metadata=None):
if all([isinstance(l, Gaussian) for l in self.likelihoods_list]):
ind = Y_metadata['output_index'].flatten()
_variance = np.array([self.likelihoods_list[j].variance for j in ind ])
if full_cov:
var += np.eye(var.shape[0])*_variance
else:
var += _variance
return mu, var
ind = Y_metadata['output_index'].flatten()
_variance = np.array([self.likelihoods_list[j].variance for j in ind ])
if full_cov:
var += np.eye(var.shape[0])*_variance
else:
raise NotImplementedError
var += _variance
return mu, var
def predictive_variance(self, mu, sigma, **other_shit):
if isinstance(noise_index,int):
_variance = self.variance[noise_index]
else:
_variance = np.array([ self.variance[j] for j in noise_index ])[:,None]
def predictive_variance(self, mu, sigma, Y_metadata):
_variance = self.gaussian_variance(Y_metadata)
return _variance + sigma**2
def predictive_quantiles(self, mu, var, quantiles, Y_metadata):
ind = Y_metadata['output_index'].flatten()
outputs = np.unique(ind)
Q = np.zeros( (mu.size,len(quantiles)) )
for j in outputs:
q = self.likelihoods_list[j].predictive_quantiles(mu[ind==j,:],
var[ind==j,:],quantiles,Y_metadata=None)
Q[ind==j,:] = np.hstack(q)
return [q[:,None] for q in Q.T]
def samples(self, gp, Y_metadata):
"""
@ -75,4 +78,3 @@ class MixedNoise(Likelihood):
_ysim = np.array([np.random.normal(lik.gp_link.transf(gpj), scale=np.sqrt(lik.variance), size=1) for gpj in gp_filtered.flatten()])
Ysim[flt,:] = _ysim.reshape(n1,N2)
return Ysim

14
GPy/models/bayesian_gplvm.py
View file

@ -72,15 +72,19 @@ class BayesianGPLVM(SparseGP):
self.variational_prior.update_gradients_KL(self.X)
def plot_latent(self, plot_inducing=True, *args, **kwargs):
"""
See GPy.plotting.matplot_dep.dim_reduction_plots.plot_latent
"""
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, **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, plot_inducing=plot_inducing, *args, **kwargs)
return dim_reduction_plots.plot_latent(self, labels, which_indices,
resolution, ax, marker, s,
fignum, plot_inducing, legend,
plot_limits, aspect, updates, **kwargs)
def do_test_latents(self, Y):
"""

2
GPy/models/gp_coregionalized_regression.py
View file

@ -36,7 +36,7 @@ class GPCoregionalizedRegression(GP):
#Kernel
if kernel is None:
kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=GPy.kern.rbf(X.shape[1]-1), W_rank=1,name=kernel_name)
kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=kern.RBF(X.shape[1]-1), W_rank=1,name=kernel_name)
#Likelihood
likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list)

16
GPy/models/gplvm.py
View file

@ -67,12 +67,22 @@ class GPLVM(GP):
assert self.likelihood.Y.shape[1] == 2
pb.scatter(self.likelihood.Y[:, 0], self.likelihood.Y[:, 1], 40, self.X[:, 0].copy(), linewidth=0, cmap=pb.cm.jet) # @UndefinedVariable
Xnew = np.linspace(self.X.min(), self.X.max(), 200)[:, None]
mu, var, upper, lower = self.predict(Xnew)
mu, _ = self.predict(Xnew)
pb.plot(mu[:, 0], mu[:, 1], 'k', linewidth=1.5)
def plot_latent(self, *args, **kwargs):
def plot_latent(self, labels=None, which_indices=None,
resolution=50, ax=None, marker='o', s=40,
fignum=None, legend=True,
plot_limits=None,
aspect='auto', updates=False, **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, *args, **kwargs)
return dim_reduction_plots.plot_latent(self, labels, which_indices,
resolution, ax, marker, s,
fignum, False, legend,
plot_limits, aspect, updates, **kwargs)
def plot_magnification(self, *args, **kwargs):
return util.plot_latent.plot_magnification(self, *args, **kwargs)

4
GPy/models/sparse_gp_coregionalized_regression.py
View file

@ -43,14 +43,14 @@ class SparseGPCoregionalizedRegression(SparseGP):
#Kernel
if kernel is None:
kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=GPy.kern.rbf(X.shape[1]-1), W_rank=1,name=kernel_name)
kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=kern.RBF(X.shape[1]-1), W_rank=1,name=kernel_name)
#Likelihood
likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list)
#Inducing inputs list
if len(Z_list):
assert len(Z_list) == self.output_dim, 'Number of outputs do not match length of inducing inputs list.'
assert len(Z_list) == Ny, 'Number of outputs do not match length of inducing inputs list.'
else:
if isinstance(num_inducing,np.int):
num_inducing = [num_inducing] * Ny

45
GPy/plotting/matplot_dep/dim_reduction_plots.py
View file

@ -30,7 +30,8 @@ def most_significant_input_dimensions(model, which_indices):
def plot_latent(model, labels=None, which_indices=None,
resolution=50, ax=None, marker='o', s=40,
fignum=None, plot_inducing=False, legend=True,
aspect='auto', updates=False):
plot_limits=None,
aspect='auto', updates=False, **kwargs):
"""
:param labels: a np.array of size model.num_data containing labels for the points (can be number, strings, etc)
:param resolution: the resolution of the grid on which to evaluate the predictive variance
@ -38,6 +39,8 @@ def plot_latent(model, labels=None, which_indices=None,
if ax is None:
fig = pb.figure(num=fignum)
ax = fig.add_subplot(111)
else:
fig = ax.figure
Tango.reset()
if labels is None:
@ -57,15 +60,28 @@ def plot_latent(model, labels=None, which_indices=None,
def plot_function(x):
Xtest_full = np.zeros((x.shape[0], model.X.shape[1]))
Xtest_full[:, [input_1, input_2]] = x
mu, var, low, up = model.predict(Xtest_full)
_, var = model.predict(Xtest_full)
var = var[:, :1]
return np.log(var)
#Create an IMshow controller that can re-plot the latent space shading at a good resolution
if plot_limits is None:
xmin, ymin = X[:, [input_1, input_2]].min(0)
xmax, ymax = X[:, [input_1, input_2]].max(0)
x_r, y_r = xmax-xmin, ymax-ymin
xmin -= .1*x_r
xmax += .1*x_r
ymin -= .1*y_r
ymax += .1*y_r
else:
try:
xmin, xmax, ymin, ymax = plot_limits
except (TypeError, ValueError) as e:
raise e.__class__, "Wrong plot limits: {} given -> need (xmin, xmax, ymin, ymax)".format(plot_limits)
view = ImshowController(ax, plot_function,
tuple(X[:, [input_1, input_2]].min(0)) + tuple(X[:, [input_1, input_2]].max(0)),
(xmin, ymin, xmax, ymax),
resolution, aspect=aspect, interpolation='bilinear',
cmap=pb.cm.binary)
cmap=pb.cm.binary, **kwargs)
# make sure labels are in order of input:
ulabels = []
@ -99,8 +115,6 @@ def plot_latent(model, labels=None, which_indices=None,
if not np.all(labels == 1.) and legend:
ax.legend(loc=0, numpoints=1)
#ax.set_xlim(xmin[0], xmax[0])
#ax.set_ylim(xmin[1], xmax[1])
ax.grid(b=False) # remove the grid if present, it doesn't look good
ax.set_aspect('auto') # set a nice aspect ratio
@ -108,9 +122,24 @@ def plot_latent(model, labels=None, which_indices=None,
Z = param_to_array(model.Z)
ax.plot(Z[:, input_1], Z[:, input_2], '^w')
ax.set_xlim((xmin, xmax))
ax.set_ylim((ymin, ymax))
try:
fig.canvas.draw()
fig.tight_layout()
fig.canvas.draw()
except Exception as e:
print "Could not invoke tight layout: {}".format(e)
pass
if updates:
ax.figure.canvas.show()
try:
ax.figure.canvas.show()
except Exception as e:
print "Could not invoke show: {}".format(e)
raw_input('Enter to continue')
view.deactivate()
return ax
def plot_magnification(model, labels=None, which_indices=None,
@ -186,7 +215,7 @@ def plot_magnification(model, labels=None, which_indices=None,
ax.plot(model.Z[:, input_1], model.Z[:, input_2], '^w')
if updates:
ax.figure.canvas.show()
fig.canvas.show()
raw_input('Enter to continue')
pb.title('Magnification Factor')

26
GPy/plotting/matplot_dep/latent_space_visualizations/controllers/axis_event_controller.py
View file

@ -33,7 +33,7 @@ class AxisChangedController(AxisEventController):
Constructor
'''
super(AxisChangedController, self).__init__(ax)
self._lim_ratio_threshold = update_lim or .8
self._lim_ratio_threshold = update_lim or .95
self._x_lim = self.ax.get_xlim()
self._y_lim = self.ax.get_ylim()
@ -80,6 +80,10 @@ class AxisChangedController(AxisEventController):
class BufferedAxisChangedController(AxisChangedController):
def __init__(self, ax, plot_function, plot_limits, resolution=50, update_lim=None, **kwargs):
"""
Buffered axis changed controller. Controls the buffer and handles update events for when the axes changed.
Updated plotting will be after first reload (first time will be within plot limits, after that the limits will be buffered)
:param plot_function:
function to use for creating image for plotting (return ndarray-like)
plot_function gets called with (2D!) Xtest grid if replotting required
@ -91,11 +95,13 @@ class BufferedAxisChangedController(AxisChangedController):
"""
super(BufferedAxisChangedController, self).__init__(ax, update_lim=update_lim)
self.plot_function = plot_function
xmin, xmax = self._x_lim # self._compute_buffered(*self._x_lim)
ymin, ymax = self._y_lim # self._compute_buffered(*self._y_lim)
xmin, ymin, xmax, ymax = plot_limits#self._x_lim # self._compute_buffered(*self._x_lim)
# imshow acts on the limits of the plot, this is why we need to override the limits here, to make sure the right plot limits are used:
self._x_lim = xmin, xmax
self._y_lim = ymin, ymax
self.resolution = resolution
self._not_init = False
self.view = self._init_view(self.ax, self.recompute_X(), xmin, xmax, ymin, ymax, **kwargs)
self.view = self._init_view(self.ax, self.recompute_X(buffered=False), xmin, xmax, ymin, ymax, **kwargs)
self._not_init = True
def update(self, ax):
@ -111,14 +117,16 @@ class BufferedAxisChangedController(AxisChangedController):
def update_view(self, view, X, xmin, xmax, ymin, ymax):
raise NotImplementedError('update view given in here')
def get_grid(self):
xmin, xmax = self._compute_buffered(*self._x_lim)
ymin, ymax = self._compute_buffered(*self._y_lim)
def get_grid(self, buffered=True):
if buffered: comp = self._compute_buffered
else: comp = lambda a,b: (a,b)
xmin, xmax = comp(*self._x_lim)
ymin, ymax = comp(*self._y_lim)
x, y = numpy.mgrid[xmin:xmax:1j * self.resolution, ymin:ymax:1j * self.resolution]
return numpy.hstack((x.flatten()[:, None], y.flatten()[:, None]))
def recompute_X(self):
X = self.plot_function(self.get_grid())
def recompute_X(self, buffered=True):
X = self.plot_function(self.get_grid(buffered))
if isinstance(X, (tuple, list)):
for x in X:
x.shape = [self.resolution, self.resolution]

2
GPy/plotting/matplot_dep/latent_space_visualizations/controllers/imshow_controller.py
View file

@ -9,7 +9,7 @@ import numpy
class ImshowController(BufferedAxisChangedController):
def __init__(self, ax, plot_function, plot_limits, resolution=50, update_lim=.5, **kwargs):
def __init__(self, ax, plot_function, plot_limits, resolution=50, update_lim=.8, **kwargs):
"""
:param plot_function:
function to use for creating image for plotting (return ndarray-like)

2
GPy/plotting/matplot_dep/models_plots.py
View file

@ -123,6 +123,8 @@ 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]
if isinstance(model,SparseGPCoregionalizedRegression):
Z = Z[Z[:,-1] == Y_metadata['output_index'],:]
Zu = Z[:,free_dims]
z_height = ax.get_ylim()[0]
plots['inducing_inputs'] = ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12)

15
GPy/plotting/matplot_dep/visualize.py
View file

@ -4,6 +4,8 @@ import GPy
import numpy as np
import matplotlib as mpl
import time
from ...util.misc import param_to_array
from GPy.core.parameterization.variational import VariationalPosterior
try:
import visual
visual_available = True
@ -72,12 +74,13 @@ class vector_show(matplotlib_show):
"""
def __init__(self, vals, axes=None):
matplotlib_show.__init__(self, vals, axes)
self.handle = self.axes.plot(np.arange(0, len(vals))[:, None], self.vals.T)[0]
self.handle = self.axes.plot(np.arange(0, len(vals))[:, None], self.vals)
def modify(self, vals):
self.vals = vals.copy()
xdata, ydata = self.handle.get_data()
self.handle.set_data(xdata, self.vals.T)
for handle, vals in zip(self.handle, self.vals.T):
xdata, ydata = handle.get_data()
handle.set_data(xdata, vals)
self.axes.figure.canvas.draw()
@ -91,8 +94,12 @@ class lvm(matplotlib_show):
:param latent_axes: the axes where the latent visualization should be plotted.
"""
if vals == None:
vals = model.X[0]
if isinstance(model.X, VariationalPosterior):
vals = param_to_array(model.X.mean)
else:
vals = param_to_array(model.X)
vals = param_to_array(vals)
matplotlib_show.__init__(self, vals, axes=latent_axes)
if isinstance(latent_axes,mpl.axes.Axes):

36
GPy/testing/kernel_tests.py
View file

@ -152,7 +152,12 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb
if verbose:
print("Checking gradients of Kdiag(X) wrt theta.")
result = Kern_check_dKdiag_dtheta(kern, X=X).checkgrad(verbose=verbose)
try:
result = Kern_check_dKdiag_dtheta(kern, X=X).checkgrad(verbose=verbose)
except NotImplementedError:
result=True
if verbose:
print("update_gradients_diag not implemented for " + kern.name)
if result and verbose:
print("Check passed.")
if not result:
@ -240,9 +245,22 @@ class KernelGradientTestsContinuous(unittest.TestCase):
def test_Add(self):
k = GPy.kern.Matern32(2, active_dims=[2,3]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D)
k += GPy.kern.Matern32(2, active_dims=[2,3]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D)
k.randomize()
self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
def test_Add_dims(self):
k = GPy.kern.Matern32(2, active_dims=[2,self.D]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D)
k.randomize()
self.assertRaises(AssertionError, k.K, self.X)
k = GPy.kern.Matern32(2, active_dims=[2,self.D-1]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D)
k.randomize()
# assert it runs:
try:
k.K(self.X)
except AssertionError:
raise AssertionError, "k.K(X) should run on self.D-1 dimension"
def test_Matern52(self):
k = GPy.kern.Matern52(self.D)
k.randomize()
@ -329,17 +347,11 @@ class KernelTestsNonContinuous(unittest.TestCase):
kern = GPy.kern.IndependentOutputs(k, -1, name='ind_split')
self.assertTrue(check_kernel_gradient_functions(kern, X=self.X, X2=self.X2, verbose=verbose, fixed_X_dims=-1))
class test_ODE_UY(unittest.TestCase):
def setUp(self):
self.k = GPy.kern.ODE_UY(2)
self.X = np.random.randn(50,2)
self.X[:,1] = np.random.randint(0,2,50)
i = np.argsort(X[:,1])
self.X = self.X[i]
self.Y = np.random.randn(50, 1)
def checkgrad(self):
m = GPy.models.GPRegression(X,Y,kernel=k)
self.assertTrue(m.checkgrad())
def test_ODE_UY(self):
kern = GPy.kern.ODE_UY(2, active_dims=[0, self.D])
X = self.X[self.X[:,-1]!=2]
X2 = self.X2[self.X2[:,-1]!=2]
self.assertTrue(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1))
if __name__ == "__main__":

5
GPy/testing/parameterized_tests.py
View file

@ -121,6 +121,11 @@ class ParameterizedTest(unittest.TestCase):
self.test1.randomize()
self.assertEqual(val, self.white.variance)
def test_randomize(self):
ps = self.test1.param.view(np.ndarray).copy()
self.test1.param.randomize()
self.assertFalse(np.all(ps==self.test1.param))
def test_fixing_randomize_parameter_handling(self):
self.rbf.fix(warning=True)
val = float(self.rbf.variance)

2
GPy/util/multioutput.py
View file

@ -56,8 +56,6 @@ def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='X'):
warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
K = kernel.prod(GPy.kern.Coregionalize(1, num_outputs, active_dims=[input_dim], rank=W_rank,W=W,kappa=kappa,name='B'),name=name)
#K = kernel * GPy.kern.Coregionalize(1, num_outputs, active_dims=[input_dim], rank=W_rank,W=W,kappa=kappa,name='B')
#K = kernel ** GPy.kern.Coregionalize(input_dim, num_outputs,W_rank,W,kappa, name= 'B')
K['.*variance'] = 1.
K['.*variance'].fix()
return K