mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-18 13:55:14 +02:00
Replaced Q by input_dim
This commit is contained in:
Alan Saul
parent
312cfebcb1
commit
97f3357b6d
22 changed files with 271 additions and 271 deletions
|
|
@ -126,7 +126,7 @@ class GP(GPBase):
|
|||
Arguments
|
||||
---------
|
||||
:param Xnew: The points at which to make a prediction
|
||||
:type Xnew: np.ndarray, Nnew x self.Q
|
||||
:type Xnew: np.ndarray, Nnew x self.input_dim
|
||||
:param which_parts: specifies which outputs kernel(s) to use in prediction
|
||||
:type which_parts: ('all', list of bools)
|
||||
:param full_cov: whether to return the folll covariance matrix, or just the diagonal
|
||||
|
|
|
|||
|
|
@ -13,7 +13,7 @@ class GPBase(model.model):
|
|||
def __init__(self, X, likelihood, kernel, normalize_X=False):
|
||||
self.X = X
|
||||
assert len(self.X.shape) == 2
|
||||
self.N, self.Q = self.X.shape
|
||||
self.N, self.input_dim = self.X.shape
|
||||
assert isinstance(kernel, kern.kern)
|
||||
self.kern = kernel
|
||||
self.likelihood = likelihood
|
||||
|
|
@ -25,8 +25,8 @@ class GPBase(model.model):
|
|||
self._Xstd = X.std(0)[None, :]
|
||||
self.X = (X.copy() - self._Xmean) / self._Xstd
|
||||
else:
|
||||
self._Xmean = np.zeros((1,self.Q))
|
||||
self._Xstd = np.ones((1,self.Q))
|
||||
self._Xmean = np.zeros((1,self.input_dim))
|
||||
self._Xstd = np.ones((1,self.input_dim))
|
||||
|
||||
model.model.__init__(self)
|
||||
|
||||
|
|
|
|||
|
|
@ -13,15 +13,15 @@ class sparse_GP(GPBase):
|
|||
Variational sparse GP model
|
||||
|
||||
:param X: inputs
|
||||
:type X: np.ndarray (N x Q)
|
||||
:type X: np.ndarray (N x input_dim)
|
||||
:param likelihood: a likelihood instance, containing the observed data
|
||||
:type likelihood: GPy.likelihood.(Gaussian | EP | Laplace)
|
||||
:param kernel : the kernel (covariance function). See link kernels
|
||||
:type kernel: a GPy.kern.kern instance
|
||||
:param X_variance: The uncertainty in the measurements of X (Gaussian variance)
|
||||
:type X_variance: np.ndarray (N x Q) | None
|
||||
:type X_variance: np.ndarray (N x input_dim) | None
|
||||
:param Z: inducing inputs (optional, see note)
|
||||
:type Z: np.ndarray (M x Q) | None
|
||||
:type Z: np.ndarray (M x input_dim) | None
|
||||
:param M : Number of inducing points (optional, default 10. Ignored if Z is not None)
|
||||
:type M: int
|
||||
:param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales)
|
||||
|
|
@ -152,7 +152,7 @@ class sparse_GP(GPBase):
|
|||
return A + B + C + D + self.likelihood.Z
|
||||
|
||||
def _set_params(self, p):
|
||||
self.Z = p[:self.M * self.Q].reshape(self.M, self.Q)
|
||||
self.Z = p[:self.M * self.input_dim].reshape(self.M, self.input_dim)
|
||||
self.kern._set_params(p[self.Z.size:self.Z.size + self.kern.Nparam])
|
||||
self.likelihood._set_params(p[self.Z.size + self.kern.Nparam:])
|
||||
self._compute_kernel_matrices()
|
||||
|
|
@ -252,9 +252,9 @@ class sparse_GP(GPBase):
|
|||
Arguments
|
||||
---------
|
||||
:param Xnew: The points at which to make a prediction
|
||||
:type Xnew: np.ndarray, Nnew x self.Q
|
||||
:type Xnew: np.ndarray, Nnew x self.input_dim
|
||||
:param X_variance_new: The uncertainty in the prediction points
|
||||
:type X_variance_new: np.ndarray, Nnew x self.Q
|
||||
:type X_variance_new: np.ndarray, Nnew x self.input_dim
|
||||
:param which_parts: specifies which outputs kernel(s) to use in prediction
|
||||
:type which_parts: ('all', list of bools)
|
||||
:param full_cov: whether to return the folll covariance matrix, or just the diagonal
|
||||
|
|
|
|||
|
|
@ -14,20 +14,20 @@ default_seed = np.random.seed(123344)
|
|||
def BGPLVM(seed=default_seed):
|
||||
N = 10
|
||||
M = 3
|
||||
Q = 2
|
||||
input_dim = 2
|
||||
D = 4
|
||||
# generate GPLVM-like data
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N), K, D).T
|
||||
|
||||
k = GPy.kern.rbf(Q, ARD=True) + GPy.kern.linear(Q, ARD=True) + GPy.kern.rbf(Q, ARD=True) + GPy.kern.white(Q)
|
||||
# k = GPy.kern.rbf(Q) + GPy.kern.rbf(Q) + GPy.kern.white(Q)
|
||||
# k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
# k = GPy.kern.rbf(Q, ARD = False) + GPy.kern.white(Q, 0.00001)
|
||||
k = GPy.kern.rbf(input_dim, ARD=True) + GPy.kern.linear(input_dim, ARD=True) + GPy.kern.rbf(input_dim, ARD=True) + GPy.kern.white(input_dim)
|
||||
# k = GPy.kern.rbf(input_dim) + GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim)
|
||||
# k = GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
# k = GPy.kern.rbf(input_dim, ARD = False) + GPy.kern.white(input_dim, 0.00001)
|
||||
|
||||
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=k, M=M)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel=k, M=M)
|
||||
m.constrain_positive('(rbf|bias|noise|white|S)')
|
||||
# m.constrain_fixed('S', 1)
|
||||
|
||||
|
|
@ -63,7 +63,7 @@ def GPLVM_oil_100(optimize=True):
|
|||
m.plot_latent(labels=m.data_labels)
|
||||
return m
|
||||
|
||||
def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
|
||||
def swiss_roll(optimize=True, N=1000, M=15, input_dim=4, sigma=.2, plot=False):
|
||||
from GPy.util.datasets import swiss_roll
|
||||
from GPy.core.transformations import logexp_clipped
|
||||
|
||||
|
|
@ -79,10 +79,10 @@ def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
|
|||
from sklearn.manifold.isomap import Isomap
|
||||
iso = Isomap().fit(Y)
|
||||
X = iso.embedding_
|
||||
if Q > 2:
|
||||
X = np.hstack((X, np.random.randn(N, Q - 2)))
|
||||
if input_dim > 2:
|
||||
X = np.hstack((X, np.random.randn(N, input_dim - 2)))
|
||||
except ImportError:
|
||||
X = np.random.randn(N, Q)
|
||||
X = np.random.randn(N, input_dim)
|
||||
|
||||
if plot:
|
||||
from mpl_toolkits import mplot3d
|
||||
|
|
@ -98,14 +98,14 @@ def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
|
|||
|
||||
|
||||
var = .5
|
||||
S = (var * np.ones_like(X) + np.clip(np.random.randn(N, Q) * var ** 2,
|
||||
S = (var * np.ones_like(X) + np.clip(np.random.randn(N, input_dim) * var ** 2,
|
||||
- (1 - var),
|
||||
(1 - var))) + .001
|
||||
Z = np.random.permutation(X)[:M]
|
||||
|
||||
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(input_dim, ARD=True) + GPy.kern.bias(input_dim, np.exp(-2)) + GPy.kern.white(input_dim, np.exp(-2))
|
||||
|
||||
m = Bayesian_GPLVM(Y, Q, X=X, X_variance=S, M=M, Z=Z, kernel=kernel)
|
||||
m = Bayesian_GPLVM(Y, input_dim, X=X, X_variance=S, M=M, Z=Z, kernel=kernel)
|
||||
m.data_colors = c
|
||||
m.data_t = t
|
||||
|
||||
|
|
@ -118,18 +118,18 @@ def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
|
|||
m.optimize('scg', messages=1)
|
||||
return m
|
||||
|
||||
def BGPLVM_oil(optimize=True, N=100, Q=5, M=25, max_f_eval=4e3, plot=False, **k):
|
||||
def BGPLVM_oil(optimize=True, N=100, input_dim=5, M=25, max_f_eval=4e3, plot=False, **k):
|
||||
np.random.seed(0)
|
||||
data = GPy.util.datasets.oil()
|
||||
from GPy.core.transformations import logexp_clipped
|
||||
|
||||
# create simple GP model
|
||||
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(input_dim, ARD=True) + GPy.kern.bias(input_dim, np.exp(-2)) + GPy.kern.white(input_dim, np.exp(-2))
|
||||
Y = data['X'][:N]
|
||||
Yn = Y - Y.mean(0)
|
||||
Yn /= Yn.std(0)
|
||||
|
||||
m = GPy.models.Bayesian_GPLVM(Yn, Q, kernel=kernel, M=M, **k)
|
||||
m = GPy.models.Bayesian_GPLVM(Yn, input_dim, kernel=kernel, M=M, **k)
|
||||
m.data_labels = data['Y'][:N].argmax(axis=1)
|
||||
|
||||
# m.constrain('variance|leng', logexp_clipped())
|
||||
|
|
@ -168,7 +168,7 @@ def oil_100():
|
|||
|
||||
|
||||
|
||||
def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False):
|
||||
def _simulate_sincos(D1, D2, D3, N, M, input_dim, plot_sim=False):
|
||||
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))
|
||||
|
|
@ -228,13 +228,13 @@ def bgplvm_simulation_matlab_compare():
|
|||
Y = sim_data['Y']
|
||||
S = sim_data['S']
|
||||
mu = sim_data['mu']
|
||||
M, [_, Q] = 3, mu.shape
|
||||
M, [_, input_dim] = 3, mu.shape
|
||||
|
||||
from GPy.models import mrd
|
||||
from GPy import kern
|
||||
reload(mrd); reload(kern)
|
||||
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
|
||||
m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k,
|
||||
k = kern.linear(input_dim, ARD=True) + kern.bias(input_dim, np.exp(-2)) + kern.white(input_dim, np.exp(-2))
|
||||
m = Bayesian_GPLVM(Y, input_dim, init="PCA", M=M, kernel=k,
|
||||
# X=mu,
|
||||
# X_variance=S,
|
||||
_debug=False)
|
||||
|
|
@ -248,8 +248,8 @@ def bgplvm_simulation(optimize='scg',
|
|||
plot=True,
|
||||
max_f_eval=2e4):
|
||||
# from GPy.core.transformations import logexp_clipped
|
||||
D1, D2, D3, N, M, Q = 15, 8, 8, 100, 3, 5
|
||||
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot)
|
||||
D1, D2, D3, N, M, input_dim = 15, 8, 8, 100, 3, 5
|
||||
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, input_dim, plot)
|
||||
|
||||
from GPy.models import mrd
|
||||
from GPy import kern
|
||||
|
|
@ -258,8 +258,8 @@ def bgplvm_simulation(optimize='scg',
|
|||
|
||||
Y = Ylist[0]
|
||||
|
||||
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q)
|
||||
m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k, _debug=True)
|
||||
k = kern.linear(input_dim, ARD=True) + kern.bias(input_dim, np.exp(-2)) + kern.white(input_dim, np.exp(-2)) # + kern.bias(input_dim)
|
||||
m = Bayesian_GPLVM(Y, input_dim, init="PCA", M=M, kernel=k, _debug=True)
|
||||
# m.constrain('variance|noise', logexp_clipped())
|
||||
m.ensure_default_constraints()
|
||||
m['noise'] = Y.var() / 100.
|
||||
|
|
@ -279,16 +279,16 @@ def bgplvm_simulation(optimize='scg',
|
|||
return m
|
||||
|
||||
def mrd_simulation(optimize=True, plot=True, plot_sim=True, **kw):
|
||||
D1, D2, D3, N, M, Q = 150, 200, 400, 500, 3, 7
|
||||
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
|
||||
D1, D2, D3, N, M, input_dim = 150, 200, 400, 500, 3, 7
|
||||
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, input_dim, plot_sim)
|
||||
|
||||
from GPy.models import mrd
|
||||
from GPy import kern
|
||||
|
||||
reload(mrd); reload(kern)
|
||||
|
||||
k = kern.linear(Q, [.05] * Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
|
||||
m = mrd.MRD(Ylist, Q=Q, M=M, kernels=k, initx="", initz='permute', **kw)
|
||||
k = kern.linear(input_dim, [.05] * input_dim, ARD=True) + kern.bias(input_dim, np.exp(-2)) + kern.white(input_dim, np.exp(-2))
|
||||
m = mrd.MRD(Ylist, input_dim=input_dim, M=M, kernels=k, initx="", initz='permute', **kw)
|
||||
|
||||
for i, Y in enumerate(Ylist):
|
||||
m['{}_noise'.format(i + 1)] = Y.var() / 100.
|
||||
|
|
@ -309,14 +309,14 @@ def mrd_simulation(optimize=True, plot=True, plot_sim=True, **kw):
|
|||
def brendan_faces():
|
||||
from GPy import kern
|
||||
data = GPy.util.datasets.brendan_faces()
|
||||
Q = 2
|
||||
input_dim = 2
|
||||
Y = data['Y'][0:-1:10, :]
|
||||
# Y = data['Y']
|
||||
Yn = Y - Y.mean()
|
||||
Yn /= Yn.std()
|
||||
|
||||
m = GPy.models.GPLVM(Yn, Q)
|
||||
# m = GPy.models.Bayesian_GPLVM(Yn, Q, M=100)
|
||||
m = GPy.models.GPLVM(Yn, input_dim)
|
||||
# m = GPy.models.Bayesian_GPLVM(Yn, input_dim, M=100)
|
||||
|
||||
# optimize
|
||||
m.constrain('rbf|noise|white', GPy.core.transformations.logexp_clipped())
|
||||
|
|
@ -379,11 +379,11 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True):
|
|||
# X -= X.mean(axis=0)
|
||||
# X /= X.std(axis=0)
|
||||
#
|
||||
# Q = 10
|
||||
# input_dim = 10
|
||||
# M = 30
|
||||
#
|
||||
# kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q) + GPy.kern.white(Q)
|
||||
# m = GPy.models.Bayesian_GPLVM(X, Q, kernel=kernel, M=M)
|
||||
# kernel = GPy.kern.rbf(input_dim, ARD=True) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim)
|
||||
# m = GPy.models.Bayesian_GPLVM(X, input_dim, kernel=kernel, M=M)
|
||||
# # m.scale_factor = 100.0
|
||||
# m.constrain_positive('(white|noise|bias|X_variance|rbf_variance|rbf_length)')
|
||||
# from sklearn import cluster
|
||||
|
|
|
|||
|
|
@ -95,11 +95,11 @@ class opt_SGD(Optimizer):
|
|||
i = 0
|
||||
|
||||
for s in param_shapes:
|
||||
N, Q = s
|
||||
X = x[i:i+N*Q].reshape(N, Q)
|
||||
N, input_dim = s
|
||||
X = x[i:i+N*input_dim].reshape(N, input_dim)
|
||||
X = X[samples]
|
||||
subset = np.append(subset, X.flatten())
|
||||
i += N*Q
|
||||
i += N*input_dim
|
||||
|
||||
subset = np.append(subset, x[i:])
|
||||
|
||||
|
|
@ -167,17 +167,17 @@ class opt_SGD(Optimizer):
|
|||
# self.model.constrained_positive_indices = p
|
||||
self.model.constrained_indices = c
|
||||
|
||||
def get_param_shapes(self, N = None, Q = None):
|
||||
def get_param_shapes(self, N = None, input_dim = None):
|
||||
model_name = self.model.__class__.__name__
|
||||
if model_name == 'GPLVM':
|
||||
return [(N, Q)]
|
||||
return [(N, input_dim)]
|
||||
if model_name == 'Bayesian_GPLVM':
|
||||
return [(N, Q), (N, Q)]
|
||||
return [(N, input_dim), (N, input_dim)]
|
||||
else:
|
||||
raise NotImplementedError
|
||||
|
||||
def step_with_missing_data(self, f_fp, X, step, shapes):
|
||||
N, Q = X.shape
|
||||
N, input_dim = X.shape
|
||||
|
||||
if not sp.sparse.issparse(self.model.likelihood.Y):
|
||||
Y = self.model.likelihood.Y
|
||||
|
|
@ -269,7 +269,7 @@ class opt_SGD(Optimizer):
|
|||
self.model.likelihood._offset = 0.0
|
||||
self.model.likelihood._scale = 1.0
|
||||
|
||||
N, Q = self.model.X.shape
|
||||
N, input_dim = self.model.X.shape
|
||||
D = self.model.likelihood.Y.shape[1]
|
||||
num_params = self.model._get_params()
|
||||
self.trace = []
|
||||
|
|
@ -302,7 +302,7 @@ class opt_SGD(Optimizer):
|
|||
self.model.likelihood._set_params(sigma)
|
||||
|
||||
if missing_data:
|
||||
shapes = self.get_param_shapes(N, Q)
|
||||
shapes = self.get_param_shapes(N, input_dim)
|
||||
f, step, Nj = self.step_with_missing_data(f_fp, X, step, shapes)
|
||||
else:
|
||||
self.model.likelihood.YYT = np.dot(self.model.likelihood.Y, self.model.likelihood.Y.T)
|
||||
|
|
|
|||
|
|
@ -300,15 +300,15 @@ class kern(parameterised):
|
|||
return target
|
||||
|
||||
def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S):
|
||||
"""return shapes are N,M,Q"""
|
||||
"""return shapes are N,M,input_dim"""
|
||||
target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1]))
|
||||
[p.dpsi1_dmuS(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
return target_mu, target_S
|
||||
|
||||
def psi2(self, Z, mu, S):
|
||||
"""
|
||||
:param Z: np.ndarray of inducing inputs (M x Q)
|
||||
:param mu, S: np.ndarrays of means and variances (each N x Q)
|
||||
:param Z: np.ndarray of inducing inputs (M x input_dim)
|
||||
:param mu, S: np.ndarrays of means and variances (each N x input_dim)
|
||||
:returns psi2: np.ndarray (N,M,M)
|
||||
"""
|
||||
target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
|
||||
|
|
|
|||
|
|
@ -168,7 +168,7 @@ class linear(kernpart):
|
|||
target += tmp.sum()
|
||||
|
||||
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
|
||||
"""Think N,M,M,Q """
|
||||
"""Think N,M,M,input_dim """
|
||||
self._psi_computations(Z, mu, S)
|
||||
AZZA = self.ZA.T[:, None, :, None] * self.ZA[None, :, None, :]
|
||||
AZZA = AZZA + AZZA.swapaxes(1, 2)
|
||||
|
|
@ -192,11 +192,11 @@ class linear(kernpart):
|
|||
else
|
||||
factor = 2.0*dL_dpsi2(n,m,mm);
|
||||
|
||||
for(q=0;q<Q;q++){
|
||||
for(q=0;q<input_dim;q++){
|
||||
|
||||
//take the dot product of mu[n,:] and AZZA[:,m,mm,q] TODO: blas!
|
||||
tmp = 0.0;
|
||||
for(qq=0;qq<Q;qq++){
|
||||
for(qq=0;qq<input_dim;qq++){
|
||||
tmp += mu(n,qq)*AZZA(qq,m,mm,q);
|
||||
}
|
||||
|
||||
|
|
@ -215,9 +215,9 @@ class linear(kernpart):
|
|||
'extra_compile_args': ['-fopenmp -O3'], #-march=native'],
|
||||
'extra_link_args' : ['-lgomp']}
|
||||
|
||||
N,M,Q = mu.shape[0],Z.shape[0],mu.shape[1]
|
||||
N,M,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
|
||||
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
||||
arg_names=['N','M','Q','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
|
||||
arg_names=['N','M','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
|
||||
type_converters=weave.converters.blitz,**weave_options)
|
||||
|
||||
|
||||
|
|
@ -232,7 +232,7 @@ class linear(kernpart):
|
|||
int n,m,mm,q;
|
||||
#pragma omp parallel for private(n,mm,q)
|
||||
for(m=0;m<M;m++){
|
||||
for(q=0;q<Q;q++){
|
||||
for(q=0;q<input_dim;q++){
|
||||
for(mm=0;mm<M;mm++){
|
||||
for(n=0;n<N;n++){
|
||||
target(m,q) += dL_dpsi2(n,m,mm)*AZA(n,mm,q);
|
||||
|
|
@ -249,9 +249,9 @@ class linear(kernpart):
|
|||
'extra_compile_args': ['-fopenmp -O3'], #-march=native'],
|
||||
'extra_link_args' : ['-lgomp']}
|
||||
|
||||
N,M,Q = mu.shape[0],Z.shape[0],mu.shape[1]
|
||||
N,M,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
|
||||
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
||||
arg_names=['N','M','Q','AZA','target','dL_dpsi2'],
|
||||
arg_names=['N','M','input_dim','AZA','target','dL_dpsi2'],
|
||||
type_converters=weave.converters.blitz,**weave_options)
|
||||
|
||||
|
||||
|
|
@ -278,7 +278,7 @@ class linear(kernpart):
|
|||
muS_changed = not (np.array_equal(mu, self._mu) and np.array_equal(S, self._S))
|
||||
if Zv_changed:
|
||||
# Z has changed, compute Z specific stuff
|
||||
# self.ZZ = Z[:,None,:]*Z[None,:,:] # M,M,Q
|
||||
# self.ZZ = Z[:,None,:]*Z[None,:,:] # M,M,input_dim
|
||||
# self.ZZ = np.empty((Z.shape[0], Z.shape[0], Z.shape[1]), order='F')
|
||||
# [tdot(Z[:, i:i + 1], self.ZZ[:, :, i].T) for i in xrange(Z.shape[1])]
|
||||
self.ZA = Z * self.variances
|
||||
|
|
@ -291,5 +291,5 @@ class linear(kernpart):
|
|||
self.inner[:, diag_indices[0], diag_indices[1]] += S
|
||||
self._mu, self._S = mu.copy(), S.copy()
|
||||
if Zv_changed or muS_changed:
|
||||
self.ZAinner = np.dot(self.ZA, self.inner).swapaxes(0, 1) # NOTE: self.ZAinner \in [M x N x Q]!
|
||||
self.ZAinner = np.dot(self.ZA, self.inner).swapaxes(0, 1) # NOTE: self.ZAinner \in [M x N x input_dim]!
|
||||
self._psi2 = np.dot(self.ZAinner, self.ZA.T)
|
||||
|
|
|
|||
|
|
@ -205,13 +205,13 @@ class rbf(kernpart):
|
|||
|
||||
def dpsi2_dZ(self,dL_dpsi2,Z,mu,S,target):
|
||||
self._psi_computations(Z,mu,S)
|
||||
term1 = self._psi2_Zdist/self.lengthscale2 # M, M, Q
|
||||
term2 = self._psi2_mudist/self._psi2_denom/self.lengthscale2 # N, M, M, Q
|
||||
term1 = self._psi2_Zdist/self.lengthscale2 # M, M, input_dim
|
||||
term2 = self._psi2_mudist/self._psi2_denom/self.lengthscale2 # N, M, M, input_dim
|
||||
dZ = self._psi2[:,:,:,None] * (term1[None] + term2)
|
||||
target += (dL_dpsi2[:,:,:,None]*dZ).sum(0).sum(0)
|
||||
|
||||
def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,target_mu,target_S):
|
||||
"""Think N,M,M,Q """
|
||||
"""Think N,M,M,input_dim """
|
||||
self._psi_computations(Z,mu,S)
|
||||
tmp = self._psi2[:,:,:,None]/self.lengthscale2/self._psi2_denom
|
||||
target_mu += -2.*(dL_dpsi2[:,:,:,None]*tmp*self._psi2_mudist).sum(1).sum(1)
|
||||
|
|
@ -242,9 +242,9 @@ class rbf(kernpart):
|
|||
#here are the "statistics" for psi1 and psi2
|
||||
if not np.array_equal(Z, self._Z):
|
||||
#Z has changed, compute Z specific stuff
|
||||
self._psi2_Zhat = 0.5*(Z[:,None,:] +Z[None,:,:]) # M,M,Q
|
||||
self._psi2_Zdist = 0.5*(Z[:,None,:]-Z[None,:,:]) # M,M,Q
|
||||
self._psi2_Zdist_sq = np.square(self._psi2_Zdist/self.lengthscale) # M,M,Q
|
||||
self._psi2_Zhat = 0.5*(Z[:,None,:] +Z[None,:,:]) # M,M,input_dim
|
||||
self._psi2_Zdist = 0.5*(Z[:,None,:]-Z[None,:,:]) # M,M,input_dim
|
||||
self._psi2_Zdist_sq = np.square(self._psi2_Zdist/self.lengthscale) # M,M,input_dim
|
||||
self._Z = Z
|
||||
|
||||
if not (np.array_equal(Z, self._Z) and np.array_equal(mu, self._mu) and np.array_equal(S, self._S)):
|
||||
|
|
@ -258,9 +258,9 @@ class rbf(kernpart):
|
|||
self._psi1 = self.variance*np.exp(self._psi1_exponent)
|
||||
|
||||
#psi2
|
||||
self._psi2_denom = 2.*S[:,None,None,:]/self.lengthscale2+1. # N,M,M,Q
|
||||
self._psi2_denom = 2.*S[:,None,None,:]/self.lengthscale2+1. # N,M,M,input_dim
|
||||
self._psi2_mudist, self._psi2_mudist_sq, self._psi2_exponent, _ = self.weave_psi2(mu,self._psi2_Zhat)
|
||||
#self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,Q
|
||||
#self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,input_dim
|
||||
#self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscale2*self._psi2_denom)
|
||||
#self._psi2_exponent = np.sum(-self._psi2_Zdist_sq -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M
|
||||
self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M
|
||||
|
|
@ -269,11 +269,11 @@ class rbf(kernpart):
|
|||
self._Z, self._mu, self._S = Z, mu,S
|
||||
|
||||
def weave_psi2(self,mu,Zhat):
|
||||
N,Q = mu.shape
|
||||
N,input_dim = mu.shape
|
||||
M = Zhat.shape[0]
|
||||
|
||||
mudist = np.empty((N,M,M,Q))
|
||||
mudist_sq = np.empty((N,M,M,Q))
|
||||
mudist = np.empty((N,M,M,input_dim))
|
||||
mudist_sq = np.empty((N,M,M,input_dim))
|
||||
psi2_exponent = np.zeros((N,M,M))
|
||||
psi2 = np.empty((N,M,M))
|
||||
|
||||
|
|
@ -284,7 +284,7 @@ class rbf(kernpart):
|
|||
if self.ARD:
|
||||
lengthscale2 = self.lengthscale2
|
||||
else:
|
||||
lengthscale2 = np.ones(Q)*self.lengthscale2
|
||||
lengthscale2 = np.ones(input_dim)*self.lengthscale2
|
||||
code = """
|
||||
double tmp;
|
||||
|
||||
|
|
@ -292,7 +292,7 @@ class rbf(kernpart):
|
|||
for (int n=0; n<N; n++){
|
||||
for (int m=0; m<M; m++){
|
||||
for (int mm=0; mm<(m+1); mm++){
|
||||
for (int q=0; q<Q; q++){
|
||||
for (int q=0; q<input_dim; q++){
|
||||
//compute mudist
|
||||
tmp = mu(n,q) - Zhat(m,mm,q);
|
||||
mudist(n,m,mm,q) = tmp;
|
||||
|
|
@ -325,7 +325,7 @@ class rbf(kernpart):
|
|||
#include <math.h>
|
||||
"""
|
||||
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
||||
arg_names=['N','M','Q','mu','Zhat','mudist_sq','mudist','lengthscale2','_psi2_denom','psi2_Zdist_sq','psi2_exponent','half_log_psi2_denom','psi2','variance_sq'],
|
||||
arg_names=['N','M','input_dim','mu','Zhat','mudist_sq','mudist','lengthscale2','_psi2_denom','psi2_Zdist_sq','psi2_exponent','half_log_psi2_denom','psi2','variance_sq'],
|
||||
type_converters=weave.converters.blitz,**self.weave_options)
|
||||
|
||||
return mudist,mudist_sq, psi2_exponent, psi2
|
||||
|
|
|
|||
|
|
@ -22,13 +22,13 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
|
||||
:param Y: observed data (np.ndarray) or GPy.likelihood
|
||||
:type Y: np.ndarray| GPy.likelihood instance
|
||||
:param Q: latent dimensionality
|
||||
:type Q: int
|
||||
:param input_dim: latent dimensionality
|
||||
:type input_dim: int
|
||||
:param init: initialisation method for the latent space
|
||||
:type init: 'PCA'|'random'
|
||||
|
||||
"""
|
||||
def __init__(self, likelihood_or_Y, Q, X=None, X_variance=None, init='PCA', M=10,
|
||||
def __init__(self, likelihood_or_Y, input_dim, X=None, X_variance=None, init='PCA', M=10,
|
||||
Z=None, kernel=None, oldpsave=10, _debug=False,
|
||||
**kwargs):
|
||||
if type(likelihood_or_Y) is np.ndarray:
|
||||
|
|
@ -37,7 +37,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
likelihood = likelihood_or_Y
|
||||
|
||||
if X == None:
|
||||
X = self.initialise_latent(init, Q, likelihood.Y)
|
||||
X = self.initialise_latent(init, input_dim, likelihood.Y)
|
||||
self.init = init
|
||||
|
||||
if X_variance is None:
|
||||
|
|
@ -48,7 +48,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
assert Z.shape[1] == X.shape[1]
|
||||
|
||||
if kernel is None:
|
||||
kernel = kern.rbf(Q) + kern.white(Q)
|
||||
kernel = kern.rbf(input_dim) + kern.white(input_dim)
|
||||
|
||||
self.oldpsave = oldpsave
|
||||
self._oldps = []
|
||||
|
|
@ -78,8 +78,8 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
self._oldps.insert(0, p.copy())
|
||||
|
||||
def _get_param_names(self):
|
||||
X_names = sum([['X_%i_%i' % (n, q) for q in range(self.Q)] for n in range(self.N)], [])
|
||||
S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.Q)] for n in range(self.N)], [])
|
||||
X_names = sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
|
||||
S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
|
||||
return (X_names + S_names + sparse_GP._get_param_names(self))
|
||||
|
||||
def _get_params(self):
|
||||
|
|
@ -101,10 +101,10 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
def _set_params(self, x, save_old=True, save_count=0):
|
||||
# try:
|
||||
x = self._clipped(x)
|
||||
N, Q = self.N, self.Q
|
||||
self.X = x[:self.X.size].reshape(N, Q).copy()
|
||||
self.X_variance = x[(N * Q):(2 * N * Q)].reshape(N, Q).copy()
|
||||
sparse_GP._set_params(self, x[(2 * N * Q):])
|
||||
N, input_dim = self.N, self.input_dim
|
||||
self.X = x[:self.X.size].reshape(N, input_dim).copy()
|
||||
self.X_variance = x[(N * input_dim):(2 * N * input_dim)].reshape(N, input_dim).copy()
|
||||
sparse_GP._set_params(self, x[(2 * N * input_dim):])
|
||||
# self.oldps = x
|
||||
# except (LinAlgError, FloatingPointError, ZeroDivisionError):
|
||||
# print "\rWARNING: Caught LinAlgError, continueing without setting "
|
||||
|
|
@ -131,7 +131,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
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.Q * self.N
|
||||
return 0.5 * (var_mean + var_S) - 0.5 * self.input_dim * self.N
|
||||
|
||||
def log_likelihood(self):
|
||||
ll = sparse_GP.log_likelihood(self)
|
||||
|
|
@ -196,16 +196,16 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
"""
|
||||
assert not self.likelihood.is_heteroscedastic
|
||||
N_test = Y.shape[0]
|
||||
Q = self.Z.shape[1]
|
||||
means = np.zeros((N_test, Q))
|
||||
covars = np.zeros((N_test, Q))
|
||||
input_dim = self.Z.shape[1]
|
||||
means = np.zeros((N_test, input_dim))
|
||||
covars = np.zeros((N_test, input_dim))
|
||||
|
||||
dpsi0 = -0.5 * self.D * self.likelihood.precision
|
||||
dpsi2 = self.dL_dpsi2[0][None, :, :] # TODO: this may change if we ignore het. likelihoods
|
||||
V = self.likelihood.precision * Y
|
||||
dpsi1 = np.dot(self.Cpsi1V, V.T)
|
||||
|
||||
start = np.zeros(self.Q * 2)
|
||||
start = np.zeros(self.input_dim * 2)
|
||||
|
||||
for n, dpsi1_n in enumerate(dpsi1.T[:, :, None]):
|
||||
args = (self.kern, self.Z, dpsi0, dpsi1_n, dpsi2)
|
||||
|
|
@ -222,12 +222,12 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
"""
|
||||
Plot latent space X in 1D:
|
||||
|
||||
-if fig is given, create Q subplots in fig and plot in these
|
||||
-if ax is given plot Q 1D latent space plots of X into each `axis`
|
||||
-if fig is given, create input_dim subplots in fig and plot in these
|
||||
-if ax is given plot input_dim 1D latent space plots of X into each `axis`
|
||||
-if neither fig nor ax is given create a figure with fignum and plot in there
|
||||
|
||||
colors:
|
||||
colors of different latent space dimensions Q
|
||||
colors of different latent space dimensions input_dim
|
||||
"""
|
||||
import pylab
|
||||
if ax is None:
|
||||
|
|
@ -245,7 +245,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
elif isinstance(ax, (tuple, list)):
|
||||
ax = ax[i]
|
||||
else:
|
||||
raise ValueError("Need one ax per latent dimnesion Q")
|
||||
raise ValueError("Need one ax per latent dimnesion input_dim")
|
||||
ax.plot(self.X, c='k', alpha=.3)
|
||||
plots.extend(ax.plot(x, self.X.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
|
||||
ax.fill_between(x,
|
||||
|
|
@ -262,7 +262,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
return fig
|
||||
|
||||
def __getstate__(self):
|
||||
return (self.likelihood, self.Q, self.X, self.X_variance,
|
||||
return (self.likelihood, self.input_dim, self.X, self.X_variance,
|
||||
self.init, self.M, self.Z, self.kern,
|
||||
self.oldpsave, self._debug)
|
||||
|
||||
|
|
@ -271,12 +271,12 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
|
||||
def _debug_filter_params(self, x):
|
||||
start, end = 0, self.X.size,
|
||||
X = x[start:end].reshape(self.N, self.Q)
|
||||
X = x[start:end].reshape(self.N, self.input_dim)
|
||||
start, end = end, end + self.X_variance.size
|
||||
X_v = x[start:end].reshape(self.N, self.Q)
|
||||
start, end = end, end + (self.M * self.Q)
|
||||
Z = x[start:end].reshape(self.M, self.Q)
|
||||
start, end = end, end + self.Q
|
||||
X_v = x[start:end].reshape(self.N, self.input_dim)
|
||||
start, end = end, end + (self.M * self.input_dim)
|
||||
Z = x[start:end].reshape(self.M, self.input_dim)
|
||||
start, end = end, end + self.input_dim
|
||||
theta = x[start:]
|
||||
return X, X_v, Z, theta
|
||||
|
||||
|
|
@ -414,24 +414,24 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
# caxkmmdl = divider.append_axes("right", "5%", pad="1%")
|
||||
# cbarkmmdl = pylab.colorbar(imkmmdl, cax=caxkmmdl)
|
||||
|
||||
# Qleg = ax1.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
|
||||
# loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.15, 1, 1.15),
|
||||
# input_dimleg = ax1.legend(Xlatentplts, [r"$input_dim_{}$".format(i + 1) for i in range(self.input_dim)],
|
||||
# loc=3, ncol=self.input_dim, bbox_to_anchor=(0, 1.15, 1, 1.15),
|
||||
# borderaxespad=0, mode="expand")
|
||||
ax2.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
|
||||
loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
|
||||
ax2.legend(Xlatentplts, [r"$input_dim_{}$".format(i + 1) for i in range(self.input_dim)],
|
||||
loc=3, ncol=self.input_dim, bbox_to_anchor=(0, 1.1, 1, 1.1),
|
||||
borderaxespad=0, mode="expand")
|
||||
ax3.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
|
||||
loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
|
||||
ax3.legend(Xlatentplts, [r"$input_dim_{}$".format(i + 1) for i in range(self.input_dim)],
|
||||
loc=3, ncol=self.input_dim, bbox_to_anchor=(0, 1.1, 1, 1.1),
|
||||
borderaxespad=0, mode="expand")
|
||||
ax4.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
|
||||
loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
|
||||
ax4.legend(Xlatentplts, [r"$input_dim_{}$".format(i + 1) for i in range(self.input_dim)],
|
||||
loc=3, ncol=self.input_dim, bbox_to_anchor=(0, 1.1, 1, 1.1),
|
||||
borderaxespad=0, mode="expand")
|
||||
ax5.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
|
||||
loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
|
||||
ax5.legend(Xlatentplts, [r"$input_dim_{}$".format(i + 1) for i in range(self.input_dim)],
|
||||
loc=3, ncol=self.input_dim, bbox_to_anchor=(0, 1.1, 1, 1.1),
|
||||
borderaxespad=0, mode="expand")
|
||||
Lleg = ax1.legend()
|
||||
Lleg.draggable()
|
||||
# ax1.add_artist(Qleg)
|
||||
# ax1.add_artist(input_dimleg)
|
||||
|
||||
indicatorKL, = ax1.plot(kllls[self.showing, 0], kllls[self.showing, 2], 'o', c=KL.get_color())
|
||||
indicatorLL, = ax1.plot(kllls[self.showing, 0], kllls[self.showing, 1] - kllls[self.showing, 2], 'o', c=LL.get_color())
|
||||
|
|
|
|||
|
|
@ -20,35 +20,35 @@ class GPLVM(GP):
|
|||
|
||||
:param Y: observed data
|
||||
:type Y: np.ndarray
|
||||
:param Q: latent dimensionality
|
||||
:type Q: int
|
||||
:param input_dim: latent dimensionality
|
||||
:type input_dim: int
|
||||
:param init: initialisation method for the latent space
|
||||
:type init: 'PCA'|'random'
|
||||
|
||||
"""
|
||||
def __init__(self, Y, Q, init='PCA', X = None, kernel=None, normalize_Y=False):
|
||||
def __init__(self, Y, input_dim, init='PCA', X = None, kernel=None, normalize_Y=False):
|
||||
if X is None:
|
||||
X = self.initialise_latent(init, Q, Y)
|
||||
X = self.initialise_latent(init, input_dim, Y)
|
||||
if kernel is None:
|
||||
kernel = kern.rbf(Q, ARD=Q>1) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
|
||||
kernel = kern.rbf(input_dim, ARD=input_dim>1) + kern.bias(input_dim, np.exp(-2)) + kern.white(input_dim, np.exp(-2))
|
||||
likelihood = Gaussian(Y, normalize=normalize_Y)
|
||||
GP.__init__(self, X, likelihood, kernel, normalize_X=False)
|
||||
self._set_params(self._get_params())
|
||||
|
||||
def initialise_latent(self, init, Q, Y):
|
||||
def initialise_latent(self, init, input_dim, Y):
|
||||
if init == 'PCA':
|
||||
return PCA(Y, Q)[0]
|
||||
return PCA(Y, input_dim)[0]
|
||||
else:
|
||||
return np.random.randn(Y.shape[0], Q)
|
||||
return np.random.randn(Y.shape[0], input_dim)
|
||||
|
||||
def _get_param_names(self):
|
||||
return sum([['X_%i_%i'%(n,q) for q in range(self.Q)] for n in range(self.N)],[]) + GP._get_param_names(self)
|
||||
return sum([['X_%i_%i'%(n,q) for q in range(self.input_dim)] for n in range(self.N)],[]) + GP._get_param_names(self)
|
||||
|
||||
def _get_params(self):
|
||||
return np.hstack((self.X.flatten(), GP._get_params(self)))
|
||||
|
||||
def _set_params(self,x):
|
||||
self.X = x[:self.N*self.Q].reshape(self.N,self.Q).copy()
|
||||
self.X = x[:self.N*self.input_dim].reshape(self.N,self.input_dim).copy()
|
||||
GP._set_params(self, x[self.X.size:])
|
||||
|
||||
def _log_likelihood_gradients(self):
|
||||
|
|
|
|||
|
|
@ -20,15 +20,15 @@ class generalized_FITC(sparse_GP):
|
|||
Naish-Guzman, A. and Holden, S. (2008) implemantation of EP with FITC.
|
||||
|
||||
:param X: inputs
|
||||
:type X: np.ndarray (N x Q)
|
||||
:type X: np.ndarray (N x input_dim)
|
||||
:param likelihood: a likelihood instance, containing the observed data
|
||||
:type likelihood: GPy.likelihood.(Gaussian | EP)
|
||||
:param kernel : the kernel/covariance function. See link kernels
|
||||
:type kernel: a GPy kernel
|
||||
:param X_variance: The variance in the measurements of X (Gaussian variance)
|
||||
:type X_variance: np.ndarray (N x Q) | None
|
||||
:type X_variance: np.ndarray (N x input_dim) | None
|
||||
:param Z: inducing inputs (optional, see note)
|
||||
:type Z: np.ndarray (M x Q) | None
|
||||
:type Z: np.ndarray (M x input_dim) | None
|
||||
:param M : Number of inducing points (optional, default 10. Ignored if Z is not None)
|
||||
:type M: int
|
||||
:param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales)
|
||||
|
|
@ -44,7 +44,7 @@ class generalized_FITC(sparse_GP):
|
|||
self._set_params(self._get_params())
|
||||
|
||||
def _set_params(self, p):
|
||||
self.Z = p[:self.M*self.Q].reshape(self.M, self.Q)
|
||||
self.Z = p[:self.M*self.input_dim].reshape(self.M, self.input_dim)
|
||||
self.kern._set_params(p[self.Z.size:self.Z.size+self.kern.Nparam])
|
||||
self.likelihood._set_params(p[self.Z.size+self.kern.Nparam:])
|
||||
self._compute_kernel_matrices()
|
||||
|
|
|
|||
|
|
@ -23,8 +23,8 @@ class MRD(model):
|
|||
:type likelihood_list: [GPy.likelihood] | [Y1..Yy]
|
||||
:param names: names for different gplvm models
|
||||
:type names: [str]
|
||||
:param Q: latent dimensionality (will raise
|
||||
:type Q: int
|
||||
:param input_dim: latent dimensionality (will raise
|
||||
:type input_dim: int
|
||||
:param initx: initialisation method for the latent space
|
||||
:type initx: 'PCA'|'random'
|
||||
:param initz: initialisation method for inducing inputs
|
||||
|
|
@ -45,7 +45,7 @@ class MRD(model):
|
|||
:param kernels: list of kernels or kernel shared for all BGPLVMS
|
||||
:type kernels: [GPy.kern.kern] | GPy.kern.kern | None (default)
|
||||
"""
|
||||
def __init__(self, likelihood_or_Y_list, Q, M=10, names=None,
|
||||
def __init__(self, likelihood_or_Y_list, input_dim, M=10, names=None,
|
||||
kernels=None, initx='PCA',
|
||||
initz='permute', _debug=False, **kw):
|
||||
if names is None:
|
||||
|
|
@ -61,14 +61,14 @@ class MRD(model):
|
|||
assert all([isinstance(k, kern) for k in kernels]), "invalid kernel object detected!"
|
||||
assert not ('kernel' in kw), "pass kernels through `kernels` argument"
|
||||
|
||||
self.Q = Q
|
||||
self.input_dim = input_dim
|
||||
self.M = M
|
||||
self._debug = _debug
|
||||
|
||||
self._init = True
|
||||
X = self._init_X(initx, likelihood_or_Y_list)
|
||||
Z = self._init_Z(initz, X)
|
||||
self.bgplvms = [Bayesian_GPLVM(l, Q=Q, kernel=k, X=X, Z=Z, M=self.M, **kw) for l, k in zip(likelihood_or_Y_list, kernels)]
|
||||
self.bgplvms = [Bayesian_GPLVM(l, input_dim=input_dim, kernel=k, X=X, Z=Z, M=self.M, **kw) for l, k in zip(likelihood_or_Y_list, kernels)]
|
||||
del self._init
|
||||
|
||||
self.gref = self.bgplvms[0]
|
||||
|
|
@ -76,8 +76,8 @@ class MRD(model):
|
|||
self.nparams = nparams.cumsum()
|
||||
|
||||
self.N = self.gref.N
|
||||
self.NQ = self.N * self.Q
|
||||
self.MQ = self.M * self.Q
|
||||
self.NQ = self.N * self.input_dim
|
||||
self.MQ = self.M * self.input_dim
|
||||
|
||||
model.__init__(self) # @UndefinedVariable
|
||||
self._set_params(self._get_params())
|
||||
|
|
@ -143,8 +143,8 @@ class MRD(model):
|
|||
self._init_Z(initz, self.X)
|
||||
|
||||
def _get_param_names(self):
|
||||
# X_names = sum([['X_%i_%i' % (n, q) for q in range(self.Q)] for n in range(self.N)], [])
|
||||
# S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.Q)] for n in range(self.N)], [])
|
||||
# X_names = sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
|
||||
# S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
|
||||
n1 = self.gref._get_param_names()
|
||||
n1var = n1[:self.NQ * 2 + self.MQ]
|
||||
map_names = lambda ns, name: map(lambda x: "{1}_{0}".format(*x),
|
||||
|
|
@ -170,9 +170,9 @@ class MRD(model):
|
|||
return params
|
||||
|
||||
# def _set_var_params(self, g, X, X_var, Z):
|
||||
# g.X = X.reshape(self.N, self.Q)
|
||||
# g.X_variance = X_var.reshape(self.N, self.Q)
|
||||
# g.Z = Z.reshape(self.M, self.Q)
|
||||
# g.X = X.reshape(self.N, self.input_dim)
|
||||
# g.X_variance = X_var.reshape(self.N, self.input_dim)
|
||||
# g.Z = Z.reshape(self.M, self.input_dim)
|
||||
#
|
||||
# def _set_kern_params(self, g, p):
|
||||
# g.kern._set_params(p[:g.kern.Nparam])
|
||||
|
|
@ -235,13 +235,13 @@ class MRD(model):
|
|||
Ylist.append(likelihood_or_Y.Y)
|
||||
del likelihood_list
|
||||
if init in "PCA_concat":
|
||||
X = PCA(numpy.hstack(Ylist), self.Q)[0]
|
||||
X = PCA(numpy.hstack(Ylist), self.input_dim)[0]
|
||||
elif init in "PCA_single":
|
||||
X = numpy.zeros((Ylist[0].shape[0], self.Q))
|
||||
for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.Q), len(Ylist)), Ylist):
|
||||
X = numpy.zeros((Ylist[0].shape[0], self.input_dim))
|
||||
for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.input_dim), len(Ylist)), Ylist):
|
||||
X[:, qs] = PCA(Y, len(qs))[0]
|
||||
else: # init == 'random':
|
||||
X = numpy.random.randn(Ylist[0].shape[0], self.Q)
|
||||
X = numpy.random.randn(Ylist[0].shape[0], self.input_dim)
|
||||
self.X = X
|
||||
return X
|
||||
|
||||
|
|
@ -252,7 +252,7 @@ class MRD(model):
|
|||
if init in "permute":
|
||||
Z = numpy.random.permutation(X.copy())[:self.M]
|
||||
elif init in "random":
|
||||
Z = numpy.random.randn(self.M, self.Q) * X.var()
|
||||
Z = numpy.random.randn(self.M, self.input_dim) * X.var()
|
||||
self.Z = Z
|
||||
return Z
|
||||
|
||||
|
|
@ -265,7 +265,7 @@ class MRD(model):
|
|||
elif isinstance(axes, (tuple, list)):
|
||||
axes = axes[i]
|
||||
else:
|
||||
raise ValueError("Need one axes per latent dimension Q")
|
||||
raise ValueError("Need one axes per latent dimension input_dim")
|
||||
plotf(i, g, axes)
|
||||
pylab.draw()
|
||||
if axes is None:
|
||||
|
|
|
|||
|
|
@ -17,25 +17,25 @@ class sparse_GPLVM(sparse_GP_regression, GPLVM):
|
|||
|
||||
:param Y: observed data
|
||||
:type Y: np.ndarray
|
||||
:param Q: latent dimensionality
|
||||
:type Q: int
|
||||
:param input_dim: latent dimensionality
|
||||
:type input_dim: int
|
||||
:param init: initialisation method for the latent space
|
||||
:type init: 'PCA'|'random'
|
||||
|
||||
"""
|
||||
def __init__(self, Y, Q, kernel=None, init='PCA', M=10):
|
||||
X = self.initialise_latent(init, Q, Y)
|
||||
def __init__(self, Y, input_dim, kernel=None, init='PCA', M=10):
|
||||
X = self.initialise_latent(init, input_dim, Y)
|
||||
sparse_GP_regression.__init__(self, X, Y, kernel=kernel,M=M)
|
||||
|
||||
def _get_param_names(self):
|
||||
return (sum([['X_%i_%i'%(n,q) for q in range(self.Q)] for n in range(self.N)],[])
|
||||
return (sum([['X_%i_%i'%(n,q) for q in range(self.input_dim)] for n in range(self.N)],[])
|
||||
+ sparse_GP_regression._get_param_names(self))
|
||||
|
||||
def _get_params(self):
|
||||
return np.hstack((self.X.flatten(), sparse_GP_regression._get_params(self)))
|
||||
|
||||
def _set_params(self,x):
|
||||
self.X = x[:self.X.size].reshape(self.N,self.Q).copy()
|
||||
self.X = x[:self.X.size].reshape(self.N,self.input_dim).copy()
|
||||
sparse_GP_regression._set_params(self, x[self.X.size:])
|
||||
|
||||
def log_likelihood(self):
|
||||
|
|
|
|||
|
|
@ -7,67 +7,67 @@ import GPy
|
|||
|
||||
class BGPLVMTests(unittest.TestCase):
|
||||
def test_bias_kern(self):
|
||||
N, M, Q, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
Y -= Y.mean(axis=0)
|
||||
k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M)
|
||||
k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel = k, M=M)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
||||
def test_linear_kern(self):
|
||||
N, M, Q, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
Y -= Y.mean(axis=0)
|
||||
k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M)
|
||||
k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel = k, M=M)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
||||
def test_rbf_kern(self):
|
||||
N, M, Q, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
Y -= Y.mean(axis=0)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel = k, M=M)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
||||
def test_rbf_bias_kern(self):
|
||||
N, M, Q, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
Y -= Y.mean(axis=0)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel = k, M=M)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
||||
#@unittest.skip('psi2 cross terms are NotImplemented for this combination')
|
||||
def test_linear_bias_kern(self):
|
||||
N, M, Q, D = 30, 5, 4, 30
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 30, 5, 4, 30
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
Y -= Y.mean(axis=0)
|
||||
k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M)
|
||||
k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel = k, M=M)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
|
|
|||
|
|
@ -7,37 +7,37 @@ import GPy
|
|||
|
||||
class GPLVMTests(unittest.TestCase):
|
||||
def test_bias_kern(self):
|
||||
N, M, Q, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.GPLVM(Y, Q, kernel = k)
|
||||
k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.GPLVM(Y, input_dim, kernel = k)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
||||
def test_linear_kern(self):
|
||||
N, M, Q, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.GPLVM(Y, Q, kernel = k)
|
||||
k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.GPLVM(Y, input_dim, kernel = k)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
||||
def test_rbf_kern(self):
|
||||
N, M, Q, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.GPLVM(Y, Q, kernel = k)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.GPLVM(Y, input_dim, kernel = k)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
|
|
|||
|
|
@ -14,16 +14,16 @@ class MRDTests(unittest.TestCase):
|
|||
|
||||
def test_gradients(self):
|
||||
num_m = 3
|
||||
N, M, Q, D = 20, 8, 6, 20
|
||||
X = np.random.rand(N, Q)
|
||||
N, M, input_dim, D = 20, 8, 6, 20
|
||||
X = np.random.rand(N, input_dim)
|
||||
|
||||
k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q)
|
||||
k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim)
|
||||
K = k.K(X)
|
||||
|
||||
Ylist = [np.random.multivariate_normal(np.zeros(N), K, D).T for _ in range(num_m)]
|
||||
likelihood_list = [GPy.likelihoods.Gaussian(Y) for Y in Ylist]
|
||||
|
||||
m = GPy.models.MRD(likelihood_list, Q=Q, kernels=k, M=M)
|
||||
m = GPy.models.MRD(likelihood_list, input_dim=input_dim, kernels=k, M=M)
|
||||
m.ensure_default_constraints()
|
||||
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
|
|
|||
|
|
@ -16,25 +16,25 @@ class PsiStatModel(model):
|
|||
self.X = X
|
||||
self.X_variance = X_variance
|
||||
self.Z = Z
|
||||
self.N, self.Q = X.shape
|
||||
self.M, Q = Z.shape
|
||||
assert self.Q == Q, "shape missmatch: Z:{!s} X:{!s}".format(Z.shape, X.shape)
|
||||
self.N, self.input_dim = X.shape
|
||||
self.M, input_dim = Z.shape
|
||||
assert self.input_dim == input_dim, "shape missmatch: Z:{!s} X:{!s}".format(Z.shape, X.shape)
|
||||
self.kern = kernel
|
||||
super(PsiStatModel, self).__init__()
|
||||
self.psi_ = self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance)
|
||||
def _get_param_names(self):
|
||||
Xnames = ["{}_{}_{}".format(what, i, j) for what, i, j in itertools.product(['X', 'X_variance'], range(self.N), range(self.Q))]
|
||||
Znames = ["Z_{}_{}".format(i, j) for i, j in itertools.product(range(self.M), range(self.Q))]
|
||||
Xnames = ["{}_{}_{}".format(what, i, j) for what, i, j in itertools.product(['X', 'X_variance'], range(self.N), range(self.input_dim))]
|
||||
Znames = ["Z_{}_{}".format(i, j) for i, j in itertools.product(range(self.M), range(self.input_dim))]
|
||||
return Xnames + Znames + self.kern._get_param_names()
|
||||
def _get_params(self):
|
||||
return numpy.hstack([self.X.flatten(), self.X_variance.flatten(), self.Z.flatten(), self.kern._get_params()])
|
||||
def _set_params(self, x, save_old=True, save_count=0):
|
||||
start, end = 0, self.X.size
|
||||
self.X = x[start:end].reshape(self.N, self.Q)
|
||||
self.X = x[start:end].reshape(self.N, self.input_dim)
|
||||
start, end = end, end + self.X_variance.size
|
||||
self.X_variance = x[start: end].reshape(self.N, self.Q)
|
||||
self.X_variance = x[start: end].reshape(self.N, self.input_dim)
|
||||
start, end = end, end + self.Z.size
|
||||
self.Z = x[start: end].reshape(self.M, self.Q)
|
||||
self.Z = x[start: end].reshape(self.M, self.input_dim)
|
||||
self.kern._set_params(x[end:])
|
||||
def log_likelihood(self):
|
||||
return self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance).sum()
|
||||
|
|
@ -43,24 +43,24 @@ class PsiStatModel(model):
|
|||
try:
|
||||
psiZ = self.kern.__getattribute__("d" + self.which + "_dZ")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance)
|
||||
except AttributeError:
|
||||
psiZ = numpy.zeros(self.M * self.Q)
|
||||
psiZ = numpy.zeros(self.M * self.input_dim)
|
||||
thetagrad = self.kern.__getattribute__("d" + self.which + "_dtheta")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance).flatten()
|
||||
return numpy.hstack((psimu.flatten(), psiS.flatten(), psiZ.flatten(), thetagrad))
|
||||
|
||||
class DPsiStatTest(unittest.TestCase):
|
||||
Q = 5
|
||||
input_dim = 5
|
||||
N = 50
|
||||
M = 10
|
||||
D = 20
|
||||
X = numpy.random.randn(N, Q)
|
||||
X = numpy.random.randn(N, input_dim)
|
||||
X_var = .5 * numpy.ones_like(X) + .4 * numpy.clip(numpy.random.randn(*X.shape), 0, 1)
|
||||
Z = numpy.random.permutation(X)[:M]
|
||||
Y = X.dot(numpy.random.randn(Q, D))
|
||||
# kernels = [GPy.kern.linear(Q, ARD=True, variances=numpy.random.rand(Q)), GPy.kern.rbf(Q, ARD=True), GPy.kern.bias(Q)]
|
||||
Y = X.dot(numpy.random.randn(input_dim, D))
|
||||
# kernels = [GPy.kern.linear(input_dim, ARD=True, variances=numpy.random.rand(input_dim)), GPy.kern.rbf(input_dim, ARD=True), GPy.kern.bias(input_dim)]
|
||||
|
||||
kernels = [GPy.kern.linear(Q), GPy.kern.rbf(Q), GPy.kern.bias(Q),
|
||||
GPy.kern.linear(Q) + GPy.kern.bias(Q),
|
||||
GPy.kern.rbf(Q) + GPy.kern.bias(Q)]
|
||||
kernels = [GPy.kern.linear(input_dim), GPy.kern.rbf(input_dim), GPy.kern.bias(input_dim),
|
||||
GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim),
|
||||
GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim)]
|
||||
|
||||
def testPsi0(self):
|
||||
for k in self.kernels:
|
||||
|
|
@ -108,31 +108,31 @@ if __name__ == "__main__":
|
|||
import sys
|
||||
interactive = 'i' in sys.argv
|
||||
if interactive:
|
||||
# N, M, Q, D = 30, 5, 4, 30
|
||||
# X = numpy.random.rand(N, Q)
|
||||
# k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
# N, M, input_dim, D = 30, 5, 4, 30
|
||||
# X = numpy.random.rand(N, input_dim)
|
||||
# k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
# K = k.K(X)
|
||||
# Y = numpy.random.multivariate_normal(numpy.zeros(N), K, D).T
|
||||
# Y -= Y.mean(axis=0)
|
||||
# k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
# m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=k, M=M)
|
||||
# k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
# m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel=k, M=M)
|
||||
# m.ensure_default_constraints()
|
||||
# m.randomize()
|
||||
# # self.assertTrue(m.checkgrad())
|
||||
numpy.random.seed(0)
|
||||
Q = 5
|
||||
input_dim = 5
|
||||
N = 50
|
||||
M = 10
|
||||
D = 15
|
||||
X = numpy.random.randn(N, Q)
|
||||
X = numpy.random.randn(N, input_dim)
|
||||
X_var = .5 * numpy.ones_like(X) + .1 * numpy.clip(numpy.random.randn(*X.shape), 0, 1)
|
||||
Z = numpy.random.permutation(X)[:M]
|
||||
Y = X.dot(numpy.random.randn(Q, D))
|
||||
# kernel = GPy.kern.bias(Q)
|
||||
Y = X.dot(numpy.random.randn(input_dim, D))
|
||||
# kernel = GPy.kern.bias(input_dim)
|
||||
#
|
||||
# kernels = [GPy.kern.linear(Q), GPy.kern.rbf(Q), GPy.kern.bias(Q),
|
||||
# GPy.kern.linear(Q) + GPy.kern.bias(Q),
|
||||
# GPy.kern.rbf(Q) + GPy.kern.bias(Q)]
|
||||
# kernels = [GPy.kern.linear(input_dim), GPy.kern.rbf(input_dim), GPy.kern.bias(input_dim),
|
||||
# GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim),
|
||||
# GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim)]
|
||||
|
||||
# for k in kernels:
|
||||
# m = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
|
||||
|
|
@ -140,18 +140,18 @@ if __name__ == "__main__":
|
|||
# assert m.checkgrad(), "{} x psi1".format("+".join(map(lambda x: x.name, k.parts)))
|
||||
#
|
||||
# m0 = PsiStatModel('psi0', X=X, X_variance=X_var, Z=Z,
|
||||
# M=M, kernel=GPy.kern.linear(Q))
|
||||
# M=M, kernel=GPy.kern.linear(input_dim))
|
||||
# m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
|
||||
# M=M, kernel=kernel)
|
||||
# m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
|
||||
# M=M, kernel=kernel)
|
||||
# m2 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
|
||||
# M=M, kernel=GPy.kern.rbf(Q))
|
||||
# M=M, kernel=GPy.kern.rbf(input_dim))
|
||||
m3 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
|
||||
M=M, kernel=GPy.kern.linear(Q, ARD=True, variances=numpy.random.rand(Q)))
|
||||
M=M, kernel=GPy.kern.linear(input_dim, ARD=True, variances=numpy.random.rand(input_dim)))
|
||||
m3.ensure_default_constraints()
|
||||
# + GPy.kern.bias(Q))
|
||||
# + GPy.kern.bias(input_dim))
|
||||
# m4 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
|
||||
# M=M, kernel=GPy.kern.rbf(Q) + GPy.kern.bias(Q))
|
||||
# M=M, kernel=GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim))
|
||||
else:
|
||||
unittest.main()
|
||||
|
|
|
|||
|
|
@ -7,38 +7,38 @@ import GPy
|
|||
|
||||
class sparse_GPLVMTests(unittest.TestCase):
|
||||
def test_bias_kern(self):
|
||||
N, M, Q, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.sparse_GPLVM(Y, Q, kernel = k, M=M)
|
||||
k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.sparse_GPLVM(Y, input_dim, kernel = k, M=M)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
||||
@unittest.skip('linear kernels do not have dKdiag_dX')
|
||||
def test_linear_kern(self):
|
||||
N, M, Q, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.sparse_GPLVM(Y, Q, kernel = k, M=M)
|
||||
k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.sparse_GPLVM(Y, input_dim, kernel = k, M=M)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
||||
def test_rbf_kern(self):
|
||||
N, M, Q, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
N, M, input_dim, D = 10, 3, 2, 4
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
m = GPy.models.sparse_GPLVM(Y, Q, kernel = k, M=M)
|
||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||
m = GPy.models.sparse_GPLVM(Y, input_dim, kernel = k, M=M)
|
||||
m.ensure_default_constraints()
|
||||
m.randomize()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
|
|
|||
|
|
@ -144,23 +144,23 @@ class GradientTests(unittest.TestCase):
|
|||
|
||||
def test_GPLVM_rbf_bias_white_kern_2D(self):
|
||||
""" Testing GPLVM with rbf + bias and white kernel """
|
||||
N, Q, D = 50, 1, 2
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q, 0.5, 0.9*np.ones((1,))) + GPy.kern.bias(Q, 0.1) + GPy.kern.white(Q, 0.05)
|
||||
N, input_dim, D = 50, 1, 2
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.rbf(input_dim, 0.5, 0.9*np.ones((1,))) + GPy.kern.bias(input_dim, 0.1) + GPy.kern.white(input_dim, 0.05)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
m = GPy.models.GPLVM(Y, Q, kernel = k)
|
||||
m = GPy.models.GPLVM(Y, input_dim, kernel = k)
|
||||
m.ensure_default_constraints()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
||||
def test_GPLVM_rbf_linear_white_kern_2D(self):
|
||||
""" Testing GPLVM with rbf + bias and white kernel """
|
||||
N, Q, D = 50, 1, 2
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.linear(Q) + GPy.kern.bias(Q, 0.1) + GPy.kern.white(Q, 0.05)
|
||||
N, input_dim, D = 50, 1, 2
|
||||
X = np.random.rand(N, input_dim)
|
||||
k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim, 0.1) + GPy.kern.white(input_dim, 0.05)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N),K,D).T
|
||||
m = GPy.models.GPLVM(Y, Q, init = 'PCA', kernel = k)
|
||||
m = GPy.models.GPLVM(Y, input_dim, init = 'PCA', kernel = k)
|
||||
m.ensure_default_constraints()
|
||||
self.assertTrue(m.checkgrad())
|
||||
|
||||
|
|
|
|||
|
|
@ -162,19 +162,19 @@ def multiple_pdinv(A):
|
|||
return np.dstack(invs),np.array(halflogdets)
|
||||
|
||||
|
||||
def PCA(Y, Q):
|
||||
def PCA(Y, input_dim):
|
||||
"""
|
||||
Principal component analysis: maximum likelihood solution by SVD
|
||||
|
||||
Arguments
|
||||
---------
|
||||
:param Y: NxD np.array of data
|
||||
:param Q: int, dimension of projection
|
||||
:param input_dim: int, dimension of projection
|
||||
|
||||
Returns
|
||||
-------
|
||||
:rval X: - NxQ np.array of dimensionality reduced data
|
||||
W - QxD mapping from X to Y
|
||||
:rval X: - Nxinput_dim np.array of dimensionality reduced data
|
||||
W - input_dimxD mapping from X to Y
|
||||
"""
|
||||
if not np.allclose(Y.mean(axis=0), 0.0):
|
||||
print "Y is not zero mean, centering it locally (GPy.util.linalg.PCA)"
|
||||
|
|
@ -182,7 +182,7 @@ def PCA(Y, Q):
|
|||
#Y -= Y.mean(axis=0)
|
||||
|
||||
Z = linalg.svd(Y-Y.mean(axis=0), full_matrices = False)
|
||||
[X, W] = [Z[0][:,0:Q], np.dot(np.diag(Z[1]), Z[2]).T[:,0:Q]]
|
||||
[X, W] = [Z[0][:,0:input_dim], np.dot(np.diag(Z[1]), Z[2]).T[:,0:input_dim]]
|
||||
v = X.std(axis=0)
|
||||
X /= v;
|
||||
W *= v;
|
||||
|
|
|
|||
|
|
@ -14,10 +14,10 @@ def plot_latent(model, labels=None, which_indices=None, resolution=50, ax=None,
|
|||
if labels is None:
|
||||
labels = np.ones(model.N)
|
||||
if which_indices is None:
|
||||
if model.Q==1:
|
||||
if model.input_dim==1:
|
||||
input_1 = 0
|
||||
input_2 = None
|
||||
if model.Q==2:
|
||||
if model.input_dim==2:
|
||||
input_1, input_2 = 0,1
|
||||
else:
|
||||
try:
|
||||
|
|
@ -55,7 +55,7 @@ def plot_latent(model, labels=None, which_indices=None, resolution=50, ax=None,
|
|||
m = marker
|
||||
|
||||
index = np.nonzero(labels==ul)[0]
|
||||
if model.Q==1:
|
||||
if model.input_dim==1:
|
||||
x = model.X[index,input_1]
|
||||
y = np.zeros(index.size)
|
||||
else:
|
||||
|
|
|
|||
|
|
@ -74,7 +74,7 @@ class lvm(data_show):
|
|||
self.called = False
|
||||
self.move_on = False
|
||||
self.latent_index = latent_index
|
||||
self.latent_dim = model.Q
|
||||
self.latent_dim = model.input_dim
|
||||
|
||||
# The red cross which shows current latent point.
|
||||
self.latent_values = vals
|
||||
|
|
@ -113,7 +113,7 @@ class lvm(data_show):
|
|||
# A click in the bar chart axis for selection a dimension.
|
||||
if self.sense_axes != None:
|
||||
self.sense_axes.cla()
|
||||
self.sense_axes.bar(np.arange(self.model.Q),1./self.model.input_sensitivity(),color='b')
|
||||
self.sense_axes.bar(np.arange(self.model.input_dim),1./self.model.input_sensitivity(),color='b')
|
||||
|
||||
if self.latent_index[1] == self.latent_index[0]:
|
||||
self.sense_axes.bar(np.array(self.latent_index[0]),1./self.model.input_sensitivity()[self.latent_index[0]],color='y')
|
||||
|
|
@ -128,11 +128,11 @@ class lvm(data_show):
|
|||
|
||||
class lvm_subplots(lvm):
|
||||
"""
|
||||
latent_axes is a np array of dimension np.ceil(Q/2),
|
||||
latent_axes is a np array of dimension np.ceil(input_dim/2),
|
||||
one for each pair of the latent dimensions.
|
||||
"""
|
||||
def __init__(self, vals, model, data_visualize, latent_axes=None, sense_axes=None):
|
||||
self.nplots = int(np.ceil(model.Q/2.))+1
|
||||
self.nplots = int(np.ceil(model.input_dim/2.))+1
|
||||
assert len(latent_axes)==self.nplots
|
||||
if vals==None:
|
||||
vals = model.X[0, :]
|
||||
|
|
@ -140,7 +140,7 @@ class lvm_subplots(lvm):
|
|||
|
||||
for i, axis in enumerate(latent_axes):
|
||||
if i == self.nplots-1:
|
||||
if self.nplots*2!=model.Q:
|
||||
if self.nplots*2!=model.input_dim:
|
||||
latent_index = [i*2, i*2]
|
||||
lvm.__init__(self, self.latent_vals, model, data_visualize, axis, sense_axes, latent_index=latent_index)
|
||||
else:
|
||||
|
|
@ -174,7 +174,7 @@ class lvm_dimselect(lvm):
|
|||
def on_click(self, event):
|
||||
|
||||
if event.inaxes==self.sense_axes:
|
||||
new_index = max(0,min(int(np.round(event.xdata-0.5)),self.model.Q-1))
|
||||
new_index = max(0,min(int(np.round(event.xdata-0.5)),self.model.input_dim-1))
|
||||
if event.button == 1:
|
||||
# Make it red if and y-axis (red=port=left) if it is a left button click
|
||||
self.latent_index[1] = new_index
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue