mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-10 12:32:40 +02:00
Multioutput models added
This commit is contained in:
Ricardo
parent
1c2a4c5c64
commit
4c7ebb6601
9 changed files with 251 additions and 62 deletions
|
|
@ -173,7 +173,7 @@ class GP(GPBase):
|
|||
This is to allow for different normalizations of the output dimensions.
|
||||
|
||||
"""
|
||||
assert isinstance(self.likelihood,EP_Mixed_Noise)
|
||||
assert hasattr(self,'multioutput')
|
||||
index = np.ones_like(Xnew)*output
|
||||
Xnew = np.hstack((Xnew,index))
|
||||
|
||||
|
|
@ -182,7 +182,10 @@ class GP(GPBase):
|
|||
mu, var = self._raw_predict(Xnew, full_cov=full_cov, which_parts=which_parts)
|
||||
|
||||
# now push through likelihood
|
||||
mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, noise_model = output)
|
||||
if isinstance(self.likelihood,EP_Mixed_Noise):
|
||||
mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, noise_model = output)
|
||||
else:
|
||||
mean, var, _025pm, _975pm = self.likelihood_list[output].predictive_values(mu, var, full_cov)
|
||||
return mean, var, _025pm, _975pm
|
||||
|
||||
def _raw_predict_single_output(self, _Xnew, output=0, which_parts='all', full_cov=False,stop=False):
|
||||
|
|
|
|||
|
|
@ -106,13 +106,16 @@ class GPBase(Model):
|
|||
gpplot(Xnew, m, m - 2 * np.sqrt(np.diag(v)[:, None]), m + 2 * np.sqrt(np.diag(v))[:, None, ], axes=ax)
|
||||
for i in range(samples):
|
||||
ax.plot(Xnew, Ysim[i, :], Tango.colorsHex['darkBlue'], linewidth=0.25)
|
||||
#ax.plot(self.X[which_data], self.likelihood.Y[which_data], 'kx', mew=1.5)
|
||||
#ax.plot(Xu[which_data], self.likelihood.Y[self.likelihood.index==output][:,None], 'kx', mew=1.5)
|
||||
ax.set_xlim(xmin, xmax)
|
||||
ymin, ymax = min(np.append(self.likelihood.Y, m - 2 * np.sqrt(np.diag(v)[:, None]))), max(np.append(self.likelihood.Y, m + 2 * np.sqrt(np.diag(v)[:, None])))
|
||||
ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin)
|
||||
ax.set_ylim(ymin, ymax)
|
||||
|
||||
if hasattr(self,'Z'):
|
||||
Zu = self.Z[self.Z[:,-1]==output,:]
|
||||
Zu = self.Z * self._Xscale + self._Xoffset
|
||||
Zu = self.Z[self.Z[:,-1]==output ,0:1] #??
|
||||
ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12)
|
||||
|
||||
else:
|
||||
raise NotImplementedError, "Cannot define a frame with more than two input dimensions"
|
||||
|
|
@ -132,7 +135,7 @@ class GPBase(Model):
|
|||
fig = pb.figure(num=fignum)
|
||||
ax = fig.add_subplot(111)
|
||||
|
||||
if self.X.shape[1] == 1 and not isinstance(self.likelihood,EP_Mixed_Noise):
|
||||
if self.X.shape[1] == 1 and not hasattr(self,'multioutput'):
|
||||
|
||||
Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now
|
||||
|
||||
|
|
@ -146,7 +149,7 @@ class GPBase(Model):
|
|||
ax.set_xlim(xmin, xmax)
|
||||
ax.set_ylim(ymin, ymax)
|
||||
|
||||
elif self.X.shape[1] == 2 and not isinstance(self.likelihood,EP_Mixed_Noise): # FIXME
|
||||
elif self.X.shape[1] == 2 and not hasattr(self,'multioutput'):
|
||||
resolution = resolution or 50
|
||||
Xnew, _, _, xmin, xmax = x_frame2D(self.X, plot_limits, resolution)
|
||||
x, y = np.linspace(xmin[0], xmax[0], resolution), np.linspace(xmin[1], xmax[1], resolution)
|
||||
|
|
@ -158,17 +161,23 @@ class GPBase(Model):
|
|||
ax.set_xlim(xmin[0], xmax[0])
|
||||
ax.set_ylim(xmin[1], xmax[1])
|
||||
|
||||
elif self.X.shape[1] == 2 and isinstance(self.likelihood,EP_Mixed_Noise):
|
||||
Xu = self.X[self.X[:,-1]==output,:]
|
||||
elif self.X.shape[1] == 2 and hasattr(self,'multioutput'):
|
||||
Xu = self.X[self.X[:,-1]==output,:] #keep the output of interest
|
||||
Xu = self.X * self._Xscale + self._Xoffset
|
||||
Xu = self.X[self.X[:,-1]==output ,0:1]
|
||||
Xu = self.X[self.X[:,-1]==output ,0:1] #get rid of the index column
|
||||
|
||||
Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits)
|
||||
|
||||
m, _, lower, upper = self.predict_single_output(Xnew, which_parts=which_parts,output=output)
|
||||
#if not isinstance(self.likelihood,EP_Mixed_Noise):
|
||||
# m, _, lower, upper = self.predict(np.hstack([Xnew,np.repeat(output,Xnew.size)[:,None]]), which_parts=which_parts)
|
||||
#else:
|
||||
# m, _, lower, upper = self.predict_single_output(Xnew, which_parts=which_parts,output=output)
|
||||
|
||||
for d in range(m.shape[1]):
|
||||
gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax)
|
||||
#ax.plot(Xu[which_data], self.likelihood.data[which_data, d], 'kx', mew=1.5)
|
||||
ax.plot(Xu[which_data], self.likelihood.data[self.likelihood.index==output][:,None], 'kx', mew=1.5)
|
||||
#ax.plot(Xu[which_data], self.likelihood.data[self.likelihood.index==output][:,None], 'kx', mew=1.5)
|
||||
ax.plot(Xu[which_data], self.likelihood_list[output].data, 'kx', mew=1.5)
|
||||
ymin, ymax = min(np.append(self.likelihood.data, lower)), max(np.append(self.likelihood.data, upper))
|
||||
ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin)
|
||||
ax.set_xlim(xmin, xmax)
|
||||
|
|
|
|||
|
|
@ -293,7 +293,7 @@ class SparseGP(GPBase):
|
|||
|
||||
return mean, var, _025pm, _975pm
|
||||
|
||||
def plot(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, fignum=None, ax=None):
|
||||
def plot(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, fignum=None, ax=None, output=None):
|
||||
if ax is None:
|
||||
fig = pb.figure(num=fignum)
|
||||
ax = fig.add_subplot(111)
|
||||
|
|
@ -301,8 +301,8 @@ class SparseGP(GPBase):
|
|||
if which_data is 'all':
|
||||
which_data = slice(None)
|
||||
|
||||
GPBase.plot(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, ax=ax)
|
||||
if self.X.shape[1] == 1:
|
||||
GPBase.plot(self, samples=0, plot_limits=plot_limits, which_data='all', which_parts='all', resolution=None, levels=20, ax=ax, output=output)
|
||||
if self.X.shape[1] == 1 and not hasattr(self,'multioutput'):
|
||||
if self.has_uncertain_inputs:
|
||||
Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now
|
||||
ax.errorbar(Xu[which_data, 0], self.likelihood.data[which_data, 0],
|
||||
|
|
@ -311,10 +311,31 @@ class SparseGP(GPBase):
|
|||
Zu = self.Z * self._Xscale + self._Xoffset
|
||||
ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12)
|
||||
|
||||
elif self.X.shape[1] == 2:
|
||||
elif self.X.shape[1] == 2 and not hasattr(self,'multioutput'):
|
||||
Zu = self.Z * self._Xscale + self._Xoffset
|
||||
ax.plot(Zu[:, 0], Zu[:, 1], 'wo')
|
||||
|
||||
elif self.X.shape[1] == 2 and hasattr(self,'multioutput'):
|
||||
Xu = self.X[self.X[:,-1]==output,:]
|
||||
if self.has_uncertain_inputs:
|
||||
Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now
|
||||
|
||||
Xu = self.X[self.X[:,-1]==output ,0:1] #??
|
||||
|
||||
ax.errorbar(Xu[which_data, 0], self.likelihood.data[which_data, 0],
|
||||
xerr=2 * np.sqrt(self.X_variance[which_data, 0]),
|
||||
ecolor='k', fmt=None, elinewidth=.5, alpha=.5)
|
||||
|
||||
Zu = self.Z[self.Z[:,-1]==output,:]
|
||||
Zu = self.Z * self._Xscale + self._Xoffset
|
||||
Zu = self.Z[self.Z[:,-1]==output ,0:1] #??
|
||||
ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12)
|
||||
#ax.set_ylim(ax.get_ylim()[0],)
|
||||
|
||||
else:
|
||||
raise NotImplementedError, "Cannot define a frame with more than two input dimensions"
|
||||
|
||||
|
||||
|
||||
def predict_single_output(self, Xnew, output=0, which_parts='all', full_cov=False):
|
||||
"""
|
||||
|
|
@ -336,7 +357,7 @@ class SparseGP(GPBase):
|
|||
This is to allow for different normalizations of the output dimensions.
|
||||
|
||||
"""
|
||||
assert isinstance(self.likelihood,EP_Mixed_Noise)
|
||||
assert hasattr(self,'multioutput')
|
||||
index = np.ones_like(Xnew)*output
|
||||
Xnew = np.hstack((Xnew,index))
|
||||
|
||||
|
|
@ -345,6 +366,51 @@ class SparseGP(GPBase):
|
|||
mu, var = self._raw_predict(Xnew, full_cov=full_cov, which_parts=which_parts)
|
||||
|
||||
# now push through likelihood
|
||||
mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, noise_model = output)
|
||||
if isinstance(self.likelihood,EP_Mixed_Noise):
|
||||
mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, noise_model = output)
|
||||
else:
|
||||
mean, var, _025pm, _975pm = self.likelihood_list[output].predictive_values(mu, var, full_cov)
|
||||
return mean, var, _025pm, _975pm
|
||||
|
||||
|
||||
|
||||
def _raw_predict_single_output(self, _Xnew, output=0, X_variance_new=None, which_parts='all', full_cov=False,stop=False):
|
||||
"""
|
||||
Internal helper function for making predictions, does not account
|
||||
for normalization or likelihood
|
||||
"""
|
||||
Bi, _ = dpotri(self.LB, lower=0) # WTH? this lower switch should be 1, but that doesn't work!
|
||||
symmetrify(Bi)
|
||||
Kmmi_LmiBLmi = backsub_both_sides(self.Lm, np.eye(self.num_inducing) - Bi)
|
||||
|
||||
if self.Cpsi1V is None:
|
||||
psi1V = np.dot(self.psi1.T,self.likelihood.V)
|
||||
tmp, _ = dtrtrs(self.Lm, np.asfortranarray(psi1V), lower=1, trans=0)
|
||||
tmp, _ = dpotrs(self.LB, tmp, lower=1)
|
||||
self.Cpsi1V, _ = dtrtrs(self.Lm, tmp, lower=1, trans=1)
|
||||
|
||||
assert hasattr(self,'multioutput')
|
||||
index = np.ones_like(_Xnew)*output
|
||||
_Xnew = np.hstack((_Xnew,index))
|
||||
|
||||
if X_variance_new is None:
|
||||
Kx = self.kern.K(self.Z, _Xnew, which_parts=which_parts)
|
||||
mu = np.dot(Kx.T, self.Cpsi1V)
|
||||
if full_cov:
|
||||
Kxx = self.kern.K(_Xnew, which_parts=which_parts)
|
||||
var = Kxx - mdot(Kx.T, Kmmi_LmiBLmi, Kx) # NOTE this won't work for plotting
|
||||
else:
|
||||
Kxx = self.kern.Kdiag(_Xnew, which_parts=which_parts)
|
||||
var = Kxx - np.sum(Kx * np.dot(Kmmi_LmiBLmi, Kx), 0)
|
||||
else:
|
||||
# assert which_p.Tarts=='all', "swithching out parts of variational kernels is not implemented"
|
||||
Kx = self.kern.psi1(self.Z, _Xnew, X_variance_new) # , which_parts=which_parts) TODO: which_parts
|
||||
mu = np.dot(Kx, self.Cpsi1V)
|
||||
if full_cov:
|
||||
raise NotImplementedError, "TODO"
|
||||
else:
|
||||
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]
|
||||
|
|
|
|||
|
|
@ -18,7 +18,7 @@ class Prod(Kernpart):
|
|||
"""
|
||||
def __init__(self,k1,k2,tensor=False):
|
||||
self.num_params = k1.num_params + k2.num_params
|
||||
self.name = k1.name + '<times>' + k2.name
|
||||
self.name = '['+k1.name + '(x)' + k2.name +']'
|
||||
self.k1 = k1
|
||||
self.k2 = k2
|
||||
if tensor:
|
||||
|
|
|
|||
|
|
@ -249,7 +249,7 @@ class EP_Mixed_Noise(likelihood):
|
|||
self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i]
|
||||
self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i]
|
||||
#Marginal moments
|
||||
self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self._transf_data[i],self.tau_[i],self.v_[i])
|
||||
self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model_list[self.index[i]].moments_match(self._transf_data[i],self.tau_[i],self.v_[i])
|
||||
#Site parameters update
|
||||
Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
|
||||
Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i])
|
||||
|
|
@ -344,7 +344,7 @@ class EP_Mixed_Noise(likelihood):
|
|||
self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i]
|
||||
self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i]
|
||||
#Marginal moments
|
||||
self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self._transf_data[i],self.tau_[i],self.v_[i])
|
||||
self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model_list[self.index[i]].moments_match(self._transf_data[i],self.tau_[i],self.v_[i])
|
||||
#Site parameters update
|
||||
Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
|
||||
Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i])
|
||||
|
|
|
|||
|
|
@ -6,6 +6,7 @@ import numpy as np
|
|||
from ..core import GP
|
||||
from .. import likelihoods
|
||||
from .. import kern
|
||||
from ..util import multioutput
|
||||
|
||||
|
||||
import pylab as pb
|
||||
|
|
@ -29,7 +30,7 @@ class GPMultioutput(GP):
|
|||
|
||||
"""
|
||||
|
||||
def __init__(self,X_list,Y_list,noise_list=[],kernel_list=None,normalize_X=False,normalize_Y=False,W=1): #TODO W
|
||||
def __init__(self,X_list,Y_list,kernel_list=None,normalize_X=False,normalize_Y=False,W=1,mixed_noise_list=[]): #TODO W
|
||||
|
||||
assert len(X_list) == len(Y_list)
|
||||
index = []
|
||||
|
|
@ -40,53 +41,30 @@ class GPMultioutput(GP):
|
|||
i += 1
|
||||
index = np.vstack(index)
|
||||
|
||||
if noise_list == []:
|
||||
likelihood_list = []
|
||||
self.likelihood_list = []
|
||||
if mixed_noise_list == []:
|
||||
for Y in Y_list:
|
||||
likelihood_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y))
|
||||
self.likelihood_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y))
|
||||
|
||||
Y = np.vstack([l_.Y for l_ in self.likelihood_list])
|
||||
likelihood = likelihoods.Gaussian(Y,normalize=False)
|
||||
likelihood.index = index
|
||||
|
||||
else:
|
||||
assert len(Y_list) == len(mixed_noise_list)
|
||||
for noise,Y in zip(mixed_noise_list,Y_list):
|
||||
self.likelihood_list.append(likelihoods.EP(Y,noise))
|
||||
likelihood = likelihoods.EP_Mixed_Noise(Y_list, mixed_noise_list)
|
||||
|
||||
Y = np.vstack([l_.Y for l_ in likelihood_list])
|
||||
likelihood = likelihoods.Gaussian(Y,normalize=False)
|
||||
likelihood.index = index
|
||||
|
||||
X = np.hstack([np.vstack(X_list),index])
|
||||
original_dim = X.shape[1] - 1
|
||||
|
||||
if kernel_list is None:
|
||||
original_dim = X.shape[1]-1
|
||||
kernel_list = [kern.rbf(original_dim) + kern.white(original_dim)]
|
||||
kernel_list = [[kern.rbf(original_dim)],[kern.white(original_dim+1)]]
|
||||
|
||||
mkernel = multioutput.build_cor_kernel(input_dim=original_dim, Nout=len(X_list), CK = kernel_list[0], NC = kernel_list[1], W=1)
|
||||
|
||||
mkernel = kernel_list[0].prod(kern.coregionalise(len(X_list),W),tensor=True)
|
||||
for k in kernel_list[1:]:
|
||||
mkernel += k.prod(kern.coregionalise(len(X_list),W),tensor=True)
|
||||
self.multioutput = True
|
||||
GP.__init__(self, X, likelihood, mkernel, normalize_X=normalize_X)
|
||||
self.ensure_default_constraints()
|
||||
|
||||
|
||||
"""
|
||||
if likelihood is None:
|
||||
noise_model_list = []
|
||||
for Y in Y_list:
|
||||
noise_model_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y))
|
||||
#noise_model_list = [likelihoods.gaussian(variance=1.) for Y in Y_list]
|
||||
#likelihood = likelihoods.EP_Mixed_Noise(Y_list, noise_model_list)
|
||||
|
||||
elif Y_list is not None:
|
||||
if not all(np.vstack(Y_list).flatten() == likelihood.data.flatten()):
|
||||
raise Warning, 'likelihood.data and Y_list values are different.'
|
||||
|
||||
X = np.hstack([np.vstack(X_list),likelihood.index])
|
||||
|
||||
if kernel_list is None:
|
||||
original_dim = X.shape[1]-1
|
||||
kernel_list = [kern.rbf(original_dim) + kern.white(original_dim)]
|
||||
|
||||
mkernel = kernel_list[0].prod(kern.coregionalise(len(X_list),W),tensor=True) #TODO W
|
||||
for k in kernel_list[1:]:
|
||||
mkernel += k.prod(kern.coregionalise(len(X_list),W),tensor=True) #TODO W
|
||||
|
||||
#kern1 = kern.rbf(1) + kern.white(1)
|
||||
#kern2 = kern.coregionalise(2,1)
|
||||
#kern3 = kern1.prod(kern2,tensor=True)
|
||||
"""
|
||||
|
||||
|
|
|
|||
97
GPy/models/sparse_gp_multioutput.py
Normal file
97
GPy/models/sparse_gp_multioutput.py
Normal file
|
|
@ -0,0 +1,97 @@
|
|||
# Copyright (c) 2013, Ricardo Andrade
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
|
||||
import numpy as np
|
||||
from ..core import SparseGP
|
||||
from .. import likelihoods
|
||||
from .. import kern
|
||||
from ..util import multioutput
|
||||
|
||||
|
||||
import pylab as pb
|
||||
|
||||
class SparseGPMultioutput(SparseGP):
|
||||
"""
|
||||
Multiple output Gaussian process
|
||||
|
||||
This is a thin wrapper around the models.GP class, with a set of sensible defaults
|
||||
|
||||
:param X_list: input observations
|
||||
:param Y_list: observed values
|
||||
:param L_list: a GPy likelihood, defaults to Binomial with probit link_function
|
||||
:param kernel_list: a GPy kernel, defaults to rbf
|
||||
:param normalize_X: whether to normalize the input data before computing (predictions will be in original scales)
|
||||
:type normalize_X: False|True
|
||||
:param normalize_Y: whether to normalize the input data before computing (predictions will be in original scales)
|
||||
:type normalize_Y: False|True
|
||||
|
||||
.. Note:: Multiple independent outputs are allowed using columns of Y
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self,X_list,Y_list,kernel_list=None,normalize_X=False,normalize_Y=False,Z_list=None,num_inducing_list=10,X_variance=None,W=1,mixed_noise_list=[]): #TODO W
|
||||
|
||||
assert len(X_list) == len(Y_list)
|
||||
index = []
|
||||
for x,y,j in zip(X_list,Y_list,range(len(X_list))):
|
||||
assert x.shape[0] == y.shape[0]
|
||||
index.append(np.repeat(j,y.size)[:,None])
|
||||
index = np.vstack(index)
|
||||
|
||||
|
||||
self.likelihood_list = []
|
||||
if mixed_noise_list == []:
|
||||
for Y in Y_list:
|
||||
self.likelihood_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y))
|
||||
|
||||
Y = np.vstack([l_.Y for l_ in self.likelihood_list])
|
||||
likelihood = likelihoods.Gaussian(Y,normalize=False)
|
||||
likelihood.index = index
|
||||
|
||||
else:
|
||||
assert len(Y_list) == len(mixed_noise_list)
|
||||
for noise,Y in zip(mixed_noise_list,Y_list):
|
||||
self.likelihood_list.append(likelihoods.EP(Y,noise))
|
||||
likelihood = likelihoods.EP_Mixed_Noise(Y_list, mixed_noise_list)
|
||||
|
||||
"""
|
||||
if noise_list == []:
|
||||
self.likelihood_list = []
|
||||
for Y in Y_list:
|
||||
self.likelihood_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y))
|
||||
|
||||
Y = np.vstack([l_.Y for l_ in self.likelihood_list])
|
||||
likelihood = likelihoods.Gaussian(Y,normalize=False)
|
||||
likelihood.index = index
|
||||
"""
|
||||
X = np.hstack([np.vstack(X_list),index])
|
||||
original_dim = X.shape[1] - 1
|
||||
|
||||
if kernel_list is None:
|
||||
kernel_list = [[kern.rbf(original_dim)],[kern.white(original_dim+1)]]
|
||||
|
||||
mkernel = multioutput.build_cor_kernel(input_dim=original_dim, Nout=len(X_list), CK = kernel_list[0], NC = kernel_list[1], W=1)
|
||||
|
||||
z_index = []
|
||||
if Z_list is None:
|
||||
if isinstance(num_inducing_list,int):
|
||||
num_inducing_list = [num_inducing_list for Xj in X_list]
|
||||
Z_list = []
|
||||
for Xj,nj,j in zip(X_list,num_inducing_list,range(len(X_list))):
|
||||
i = np.random.permutation(Xj.shape[0])[:nj]
|
||||
z_index.append(np.repeat(j,nj)[:,None])
|
||||
Z_list.append(Xj[i].copy())
|
||||
else:
|
||||
assert len(Z_list) == len(X_list)
|
||||
for Zj,Xj,j in zip(Z_list,X_list,range(len(Z_list))):
|
||||
assert Zj.shape[1] == Xj.shape[1]
|
||||
z_index.append(np.repeat(j,Zj.shape[0])[:,None])
|
||||
|
||||
Z = np.hstack([np.vstack(Z_list),np.vstack(z_index)])
|
||||
|
||||
|
||||
self.multioutput = True
|
||||
SparseGP.__init__(self, X, likelihood, mkernel, Z=Z, normalize_X=normalize_X, X_variance=X_variance)
|
||||
self.constrain_fixed('.*iip_\d+_1')
|
||||
self.ensure_default_constraints()
|
||||
|
|
@ -14,3 +14,4 @@ import mocap
|
|||
import visualize
|
||||
import decorators
|
||||
import classification
|
||||
import multioutput
|
||||
|
|
|
|||
35
GPy/util/multioutput.py
Normal file
35
GPy/util/multioutput.py
Normal file
|
|
@ -0,0 +1,35 @@
|
|||
import numpy as np
|
||||
import warnings
|
||||
from .. import kern
|
||||
|
||||
def build_cor_kernel(input_dim, Nout, CK = [], NC = [], W=1):
|
||||
"""
|
||||
Builds an appropiate coregionalized kernel
|
||||
|
||||
:input_dim: Input dimensionality
|
||||
:Nout: Number of outputs
|
||||
:param CK: List of coregionalized kernels (i.e., this will be multiplied by a coregionalise kernel).
|
||||
:param K: List of kernels that won't be multiplied by a coregionalise kernel
|
||||
:W:
|
||||
"""
|
||||
|
||||
for k in CK:
|
||||
if k.input_dim <> input_dim:
|
||||
k.input_dim = input_dim
|
||||
#raise Warning("kernel's input dimension overwritten to fit input_dim parameter.")
|
||||
warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
|
||||
|
||||
for k in NC:
|
||||
if k.input_dim <> input_dim + 1:
|
||||
k.input_dim = input_dim + 1
|
||||
#raise Warning("kernel's input dimension overwritten to fit input_dim parameter.")
|
||||
warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
|
||||
|
||||
kernel = CK[0].prod(kern.coregionalise(Nout,W),tensor=True)
|
||||
for k in CK[1:]:
|
||||
kernel += k.prod(kern.coregionalise(Nout,W),tensor=True)
|
||||
|
||||
for k in NC:
|
||||
kernel += k
|
||||
|
||||
return kernel
|
||||
Loading…
Add table
Add a link
Reference in a new issue