mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-11 13:02:38 +02:00
Merge kern conflicts in examples
This commit is contained in:
commit
bfd99c3607
21 changed files with 156 additions and 157 deletions
|
|
@ -12,7 +12,7 @@ default_seed = np.random.seed(123344)
|
||||||
|
|
||||||
def BGPLVM(seed=default_seed):
|
def BGPLVM(seed=default_seed):
|
||||||
N = 10
|
N = 10
|
||||||
M = 3
|
num_inducing = 3
|
||||||
Q = 2
|
Q = 2
|
||||||
D = 4
|
D = 4
|
||||||
# generate GPLVM-like data
|
# generate GPLVM-like data
|
||||||
|
|
@ -26,7 +26,7 @@ def BGPLVM(seed=default_seed):
|
||||||
# k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
# 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(Q, ARD = False) + GPy.kern.white(Q, 0.00001)
|
||||||
|
|
||||||
m = GPy.models.BayesianGPLVM(Y, Q, kernel=k, M=M)
|
m = GPy.models.BayesianGPLVM(Y, Q, kernel=k, num_inducing=num_inducing)
|
||||||
m.constrain_positive('(rbf|bias|noise|white|S)')
|
m.constrain_positive('(rbf|bias|noise|white|S)')
|
||||||
# m.constrain_fixed('S', 1)
|
# m.constrain_fixed('S', 1)
|
||||||
|
|
||||||
|
|
@ -62,7 +62,7 @@ def GPLVM_oil_100(optimize=True):
|
||||||
m.plot_latent(labels=m.data_labels)
|
m.plot_latent(labels=m.data_labels)
|
||||||
return m
|
return m
|
||||||
|
|
||||||
def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
|
def swiss_roll(optimize=True, N=1000, num_inducing=15, Q=4, sigma=.2, plot=False):
|
||||||
from GPy.util.datasets import swiss_roll_generated
|
from GPy.util.datasets import swiss_roll_generated
|
||||||
from GPy.core.transformations import logexp_clipped
|
from GPy.core.transformations import logexp_clipped
|
||||||
|
|
||||||
|
|
@ -100,11 +100,11 @@ def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
|
||||||
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, Q) * var ** 2,
|
||||||
- (1 - var),
|
- (1 - var),
|
||||||
(1 - var))) + .001
|
(1 - var))) + .001
|
||||||
Z = np.random.permutation(X)[:M]
|
Z = np.random.permutation(X)[:num_inducing]
|
||||||
|
|
||||||
kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
|
kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
|
||||||
|
|
||||||
m = BayesianGPLVM(Y, Q, X=X, X_variance=S, M=M, Z=Z, kernel=kernel)
|
m = BayesianGPLVM(Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel)
|
||||||
m.data_colors = c
|
m.data_colors = c
|
||||||
m.data_t = t
|
m.data_t = t
|
||||||
|
|
||||||
|
|
@ -117,7 +117,7 @@ def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
|
||||||
m.optimize('scg', messages=1)
|
m.optimize('scg', messages=1)
|
||||||
return m
|
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, Q=5, num_inducing=25, max_f_eval=4e3, plot=False, **k):
|
||||||
np.random.seed(0)
|
np.random.seed(0)
|
||||||
data = GPy.util.datasets.oil()
|
data = GPy.util.datasets.oil()
|
||||||
from GPy.core.transformations import logexp_clipped
|
from GPy.core.transformations import logexp_clipped
|
||||||
|
|
@ -128,7 +128,7 @@ def BGPLVM_oil(optimize=True, N=100, Q=5, M=25, max_f_eval=4e3, plot=False, **k)
|
||||||
Yn = Y - Y.mean(0)
|
Yn = Y - Y.mean(0)
|
||||||
Yn /= Yn.std(0)
|
Yn /= Yn.std(0)
|
||||||
|
|
||||||
m = GPy.models.BayesianGPLVM(Yn, Q, kernel=kernel, M=M, **k)
|
m = GPy.models.BayesianGPLVM(Yn, Q, kernel=kernel, num_inducing=num_inducing, **k)
|
||||||
m.data_labels = data['Y'][:N].argmax(axis=1)
|
m.data_labels = data['Y'][:N].argmax(axis=1)
|
||||||
|
|
||||||
# m.constrain('variance|leng', logexp_clipped())
|
# m.constrain('variance|leng', logexp_clipped())
|
||||||
|
|
@ -167,7 +167,7 @@ def oil_100():
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False):
|
def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False):
|
||||||
x = np.linspace(0, 4 * np.pi, N)[:, None]
|
x = np.linspace(0, 4 * np.pi, N)[:, None]
|
||||||
s1 = np.vectorize(lambda x: np.sin(x))
|
s1 = np.vectorize(lambda x: np.sin(x))
|
||||||
s2 = np.vectorize(lambda x: np.cos(x))
|
s2 = np.vectorize(lambda x: np.cos(x))
|
||||||
|
|
@ -227,13 +227,13 @@ def bgplvm_simulation_matlab_compare():
|
||||||
Y = sim_data['Y']
|
Y = sim_data['Y']
|
||||||
S = sim_data['S']
|
S = sim_data['S']
|
||||||
mu = sim_data['mu']
|
mu = sim_data['mu']
|
||||||
M, [_, Q] = 3, mu.shape
|
num_inducing, [_, Q] = 3, mu.shape
|
||||||
|
|
||||||
from GPy.models import mrd
|
from GPy.models import mrd
|
||||||
from GPy import kern
|
from GPy import kern
|
||||||
reload(mrd); reload(kern)
|
reload(mrd); reload(kern)
|
||||||
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
|
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
|
||||||
m = BayesianGPLVM(Y, Q, init="PCA", M=M, kernel=k,
|
m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k,
|
||||||
# X=mu,
|
# X=mu,
|
||||||
# X_variance=S,
|
# X_variance=S,
|
||||||
_debug=False)
|
_debug=False)
|
||||||
|
|
@ -247,8 +247,8 @@ def bgplvm_simulation(optimize='scg',
|
||||||
plot=True,
|
plot=True,
|
||||||
max_f_eval=2e4):
|
max_f_eval=2e4):
|
||||||
# from GPy.core.transformations import logexp_clipped
|
# from GPy.core.transformations import logexp_clipped
|
||||||
D1, D2, D3, N, M, Q = 15, 8, 8, 100, 3, 5
|
D1, D2, D3, N, num_inducing, Q = 15, 8, 8, 100, 3, 5
|
||||||
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot)
|
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot)
|
||||||
|
|
||||||
from GPy.models import mrd
|
from GPy.models import mrd
|
||||||
from GPy import kern
|
from GPy import kern
|
||||||
|
|
@ -258,7 +258,7 @@ def bgplvm_simulation(optimize='scg',
|
||||||
Y = Ylist[0]
|
Y = Ylist[0]
|
||||||
|
|
||||||
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q)
|
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q)
|
||||||
m = BayesianGPLVM(Y, Q, init="PCA", M=M, kernel=k, _debug=True)
|
m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k, _debug=True)
|
||||||
# m.constrain('variance|noise', logexp_clipped())
|
# m.constrain('variance|noise', logexp_clipped())
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
m['noise'] = Y.var() / 100.
|
m['noise'] = Y.var() / 100.
|
||||||
|
|
@ -275,8 +275,8 @@ def bgplvm_simulation(optimize='scg',
|
||||||
return m
|
return m
|
||||||
|
|
||||||
def mrd_simulation(optimize=True, plot=True, plot_sim=True, **kw):
|
def mrd_simulation(optimize=True, plot=True, plot_sim=True, **kw):
|
||||||
D1, D2, D3, N, M, Q = 150, 200, 400, 500, 3, 7
|
D1, D2, D3, N, num_inducing, Q = 150, 200, 400, 500, 3, 7
|
||||||
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
|
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim)
|
||||||
|
|
||||||
from GPy.models import mrd
|
from GPy.models import mrd
|
||||||
from GPy import kern
|
from GPy import kern
|
||||||
|
|
@ -284,7 +284,7 @@ def mrd_simulation(optimize=True, plot=True, plot_sim=True, **kw):
|
||||||
reload(mrd); reload(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))
|
k = kern.linear(Q, [.05] * Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
|
||||||
m = mrd.MRD(Ylist, input_dim=Q, M=M, kernels=k, initx="", initz='permute', **kw)
|
m = mrd.MRD(Ylist, input_dim=Q, num_inducing=num_inducing, kernels=k, initx="", initz='permute', **kw)
|
||||||
|
|
||||||
for i, Y in enumerate(Ylist):
|
for i, Y in enumerate(Ylist):
|
||||||
m['{}_noise'.format(i + 1)] = Y.var() / 100.
|
m['{}_noise'.format(i + 1)] = Y.var() / 100.
|
||||||
|
|
@ -312,7 +312,7 @@ def brendan_faces():
|
||||||
Yn /= Yn.std()
|
Yn /= Yn.std()
|
||||||
|
|
||||||
m = GPy.models.GPLVM(Yn, Q)
|
m = GPy.models.GPLVM(Yn, Q)
|
||||||
# m = GPy.models.BayesianGPLVM(Yn, Q, M=100)
|
# m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=100)
|
||||||
|
|
||||||
# optimize
|
# optimize
|
||||||
m.constrain('rbf|noise|white', GPy.core.transformations.logexp_clipped())
|
m.constrain('rbf|noise|white', GPy.core.transformations.logexp_clipped())
|
||||||
|
|
@ -376,16 +376,16 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True):
|
||||||
# X /= X.std(axis=0)
|
# X /= X.std(axis=0)
|
||||||
#
|
#
|
||||||
# Q = 10
|
# Q = 10
|
||||||
# M = 30
|
# num_inducing = 30
|
||||||
#
|
#
|
||||||
# kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q) + GPy.kern.white(Q)
|
# kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q) + GPy.kern.white(Q)
|
||||||
# m = GPy.models.BayesianGPLVM(X, Q, kernel=kernel, M=M)
|
# m = GPy.models.BayesianGPLVM(X, Q, kernel=kernel, num_inducing=num_inducing)
|
||||||
# # m.scale_factor = 100.0
|
# # m.scale_factor = 100.0
|
||||||
# m.constrain_positive('(white|noise|bias|X_variance|rbf_variance|rbf_length)')
|
# m.constrain_positive('(white|noise|bias|X_variance|rbf_variance|rbf_length)')
|
||||||
# from sklearn import cluster
|
# from sklearn import cluster
|
||||||
# km = cluster.KMeans(M, verbose=10)
|
# km = cluster.KMeans(num_inducing, verbose=10)
|
||||||
# Z = km.fit(m.X).cluster_centers_
|
# Z = km.fit(m.X).cluster_centers_
|
||||||
# # Z = GPy.util.misc.kmm_init(m.X, M)
|
# # Z = GPy.util.misc.kmm_init(m.X, num_inducing)
|
||||||
# m.set('iip', Z)
|
# m.set('iip', Z)
|
||||||
# m.set('bias', 1e-4)
|
# m.set('bias', 1e-4)
|
||||||
# # optimize
|
# # optimize
|
||||||
|
|
|
||||||
|
|
@ -151,8 +151,8 @@ def coregionalisation_sparse(optim_iters=100):
|
||||||
Y2 = -np.sin(X2) + np.random.randn(*X2.shape)*0.05
|
Y2 = -np.sin(X2) + np.random.randn(*X2.shape)*0.05
|
||||||
Y = np.vstack((Y1,Y2))
|
Y = np.vstack((Y1,Y2))
|
||||||
|
|
||||||
M = 40
|
num_inducing = 40
|
||||||
Z = np.hstack((np.random.rand(M,1)*8,np.random.randint(0,2,M)[:,None]))
|
Z = np.hstack((np.random.rand(num_inducing,1)*8,np.random.randint(0,2,num_inducing)[:,None]))
|
||||||
|
|
||||||
k1 = GPy.kern.rbf(1)
|
k1 = GPy.kern.rbf(1)
|
||||||
k2 = GPy.kern.Coregionalise(2,2)
|
k2 = GPy.kern.Coregionalise(2,2)
|
||||||
|
|
@ -261,7 +261,7 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.rbf):
|
||||||
|
|
||||||
return np.array(lls)
|
return np.array(lls)
|
||||||
|
|
||||||
def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100):
|
def sparse_GP_regression_1D(N = 400, num_inducing = 5, optim_iters=100):
|
||||||
"""Run a 1D example of a sparse GP regression."""
|
"""Run a 1D example of a sparse GP regression."""
|
||||||
# sample inputs and outputs
|
# sample inputs and outputs
|
||||||
X = np.random.uniform(-3.,3.,(N,1))
|
X = np.random.uniform(-3.,3.,(N,1))
|
||||||
|
|
@ -271,7 +271,7 @@ def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100):
|
||||||
noise = GPy.kern.white(1)
|
noise = GPy.kern.white(1)
|
||||||
kernel = rbf + noise
|
kernel = rbf + noise
|
||||||
# create simple GP Model
|
# create simple GP Model
|
||||||
m = GPy.models.SparseGPRegression(X, Y, kernel, M=M)
|
m = GPy.models.SparseGPRegression(X, Y, kernel, num_inducing=num_inducing)
|
||||||
|
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
|
|
||||||
|
|
@ -280,7 +280,7 @@ def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100):
|
||||||
m.plot()
|
m.plot()
|
||||||
return m
|
return m
|
||||||
|
|
||||||
def sparse_GP_regression_2D(N = 400, M = 50, optim_iters=100):
|
def sparse_GP_regression_2D(N = 400, num_inducing = 50, optim_iters=100):
|
||||||
"""Run a 2D example of a sparse GP regression."""
|
"""Run a 2D example of a sparse GP regression."""
|
||||||
X = np.random.uniform(-3.,3.,(N,2))
|
X = np.random.uniform(-3.,3.,(N,2))
|
||||||
Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(N,1)*0.05
|
Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(N,1)*0.05
|
||||||
|
|
@ -291,7 +291,7 @@ def sparse_GP_regression_2D(N = 400, M = 50, optim_iters=100):
|
||||||
kernel = rbf + noise
|
kernel = rbf + noise
|
||||||
|
|
||||||
# create simple GP Model
|
# create simple GP Model
|
||||||
m = GPy.models.SparseGPRegression(X,Y,kernel, M = M)
|
m = GPy.models.SparseGPRegression(X,Y,kernel, num_inducing = num_inducing)
|
||||||
|
|
||||||
# contrain all parameters to be positive (but not inducing inputs)
|
# contrain all parameters to be positive (but not inducing inputs)
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
|
|
|
||||||
|
|
@ -69,14 +69,14 @@ class Coregionalise(Kernpart):
|
||||||
else:
|
else:
|
||||||
index2 = np.asarray(index2,dtype=np.int)
|
index2 = np.asarray(index2,dtype=np.int)
|
||||||
code="""
|
code="""
|
||||||
for(int i=0;i<M; i++){
|
for(int i=0;i<num_inducing; i++){
|
||||||
for(int j=0; j<N; j++){
|
for(int j=0; j<N; j++){
|
||||||
target[i+j*M] += B[Nout*index[j]+index2[i]];
|
target[i+j*num_inducing] += B[Nout*index[j]+index2[i]];
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
"""
|
"""
|
||||||
N,M,B,Nout = index.size,index2.size, self.B, self.Nout
|
N,num_inducing,B,Nout = index.size,index2.size, self.B, self.Nout
|
||||||
weave.inline(code,['target','index','index2','N','M','B','Nout'])
|
weave.inline(code,['target','index','index2','N','num_inducing','B','Nout'])
|
||||||
|
|
||||||
|
|
||||||
def Kdiag(self,index,target):
|
def Kdiag(self,index,target):
|
||||||
|
|
@ -91,14 +91,14 @@ class Coregionalise(Kernpart):
|
||||||
index2 = np.asarray(index2,dtype=np.int)
|
index2 = np.asarray(index2,dtype=np.int)
|
||||||
|
|
||||||
code="""
|
code="""
|
||||||
for(int i=0; i<M; i++){
|
for(int i=0; i<num_inducing; i++){
|
||||||
for(int j=0; j<N; j++){
|
for(int j=0; j<N; j++){
|
||||||
dL_dK_small[index[j] + Nout*index2[i]] += dL_dK[i+j*M];
|
dL_dK_small[index[j] + Nout*index2[i]] += dL_dK[i+j*num_inducing];
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
"""
|
"""
|
||||||
N, M, Nout = index.size, index2.size, self.Nout
|
N, num_inducing, Nout = index.size, index2.size, self.Nout
|
||||||
weave.inline(code, ['N','M','Nout','dL_dK','dL_dK_small','index','index2'])
|
weave.inline(code, ['N','num_inducing','Nout','dL_dK','dL_dK_small','index','index2'])
|
||||||
|
|
||||||
dkappa = np.diag(dL_dK_small)
|
dkappa = np.diag(dL_dK_small)
|
||||||
dL_dK_small += dL_dK_small.T
|
dL_dK_small += dL_dK_small.T
|
||||||
|
|
|
||||||
|
|
@ -222,11 +222,11 @@ class kern(Parameterised):
|
||||||
def dK_dtheta(self, dL_dK, X, X2=None):
|
def dK_dtheta(self, dL_dK, X, X2=None):
|
||||||
"""
|
"""
|
||||||
:param dL_dK: An array of dL_dK derivaties, dL_dK
|
:param dL_dK: An array of dL_dK derivaties, dL_dK
|
||||||
:type dL_dK: Np.ndarray (N x M)
|
:type dL_dK: Np.ndarray (N x num_inducing)
|
||||||
:param X: Observed data inputs
|
:param X: Observed data inputs
|
||||||
:type X: np.ndarray (N x input_dim)
|
:type X: np.ndarray (N x input_dim)
|
||||||
:param X2: Observed dara inputs (optional, defaults to X)
|
:param X2: Observed dara inputs (optional, defaults to X)
|
||||||
:type X2: np.ndarray (M x input_dim)
|
:type X2: np.ndarray (num_inducing x input_dim)
|
||||||
"""
|
"""
|
||||||
assert X.shape[1] == self.input_dim
|
assert X.shape[1] == self.input_dim
|
||||||
target = np.zeros(self.num_params)
|
target = np.zeros(self.num_params)
|
||||||
|
|
@ -299,16 +299,16 @@ class kern(Parameterised):
|
||||||
return target
|
return target
|
||||||
|
|
||||||
def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S):
|
def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S):
|
||||||
"""return shapes are N,M,input_dim"""
|
"""return shapes are N,num_inducing,input_dim"""
|
||||||
target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1]))
|
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)]
|
[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
|
return target_mu, target_S
|
||||||
|
|
||||||
def psi2(self, Z, mu, S):
|
def psi2(self, Z, mu, S):
|
||||||
"""
|
"""
|
||||||
:param Z: np.ndarray of inducing inputs (M x input_dim)
|
:param Z: np.ndarray of inducing inputs (num_inducing x input_dim)
|
||||||
:param mu, S: np.ndarrays of means and variances (each N x input_dim)
|
:param mu, S: np.ndarrays of means and variances (each N x input_dim)
|
||||||
:returns psi2: np.ndarray (N,M,M)
|
:returns psi2: np.ndarray (N,num_inducing,num_inducing)
|
||||||
"""
|
"""
|
||||||
target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
|
target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
|
||||||
[p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
|
[p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
|
||||||
|
|
|
||||||
|
|
@ -138,7 +138,7 @@ class linear(Kernpart):
|
||||||
|
|
||||||
def psi2(self, Z, mu, S, target):
|
def psi2(self, Z, mu, S, target):
|
||||||
"""
|
"""
|
||||||
returns N,M,M matrix
|
returns N,num_inducing,num_inducing matrix
|
||||||
"""
|
"""
|
||||||
self._psi_computations(Z, mu, S)
|
self._psi_computations(Z, mu, S)
|
||||||
# psi2_old = self.ZZ * np.square(self.variances) * self.mu2_S[:, None, None, :]
|
# psi2_old = self.ZZ * np.square(self.variances) * self.mu2_S[:, None, None, :]
|
||||||
|
|
@ -168,7 +168,7 @@ class linear(Kernpart):
|
||||||
target += tmp.sum()
|
target += tmp.sum()
|
||||||
|
|
||||||
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
|
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
|
||||||
"""Think N,M,M,input_dim """
|
"""Think N,num_inducing,num_inducing,input_dim """
|
||||||
self._psi_computations(Z, mu, S)
|
self._psi_computations(Z, mu, S)
|
||||||
AZZA = self.ZA.T[:, None, :, None] * self.ZA[None, :, None, :]
|
AZZA = self.ZA.T[:, None, :, None] * self.ZA[None, :, None, :]
|
||||||
AZZA = AZZA + AZZA.swapaxes(1, 2)
|
AZZA = AZZA + AZZA.swapaxes(1, 2)
|
||||||
|
|
@ -184,7 +184,7 @@ class linear(Kernpart):
|
||||||
double factor,tmp;
|
double factor,tmp;
|
||||||
#pragma omp parallel for private(m,mm,q,qq,factor,tmp)
|
#pragma omp parallel for private(m,mm,q,qq,factor,tmp)
|
||||||
for(n=0;n<N;n++){
|
for(n=0;n<N;n++){
|
||||||
for(m=0;m<M;m++){
|
for(m=0;m<num_inducing;m++){
|
||||||
for(mm=0;mm<=m;mm++){
|
for(mm=0;mm<=m;mm++){
|
||||||
//add in a factor of 2 for the off-diagonal terms (and then count them only once)
|
//add in a factor of 2 for the off-diagonal terms (and then count them only once)
|
||||||
if(m==mm)
|
if(m==mm)
|
||||||
|
|
@ -215,9 +215,9 @@ class linear(Kernpart):
|
||||||
'extra_compile_args': ['-fopenmp -O3'], #-march=native'],
|
'extra_compile_args': ['-fopenmp -O3'], #-march=native'],
|
||||||
'extra_link_args' : ['-lgomp']}
|
'extra_link_args' : ['-lgomp']}
|
||||||
|
|
||||||
N,M,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
|
N,num_inducing,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
|
||||||
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
||||||
arg_names=['N','M','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
|
arg_names=['N','num_inducing','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
|
||||||
type_converters=weave.converters.blitz,**weave_options)
|
type_converters=weave.converters.blitz,**weave_options)
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -231,9 +231,9 @@ class linear(Kernpart):
|
||||||
code="""
|
code="""
|
||||||
int n,m,mm,q;
|
int n,m,mm,q;
|
||||||
#pragma omp parallel for private(n,mm,q)
|
#pragma omp parallel for private(n,mm,q)
|
||||||
for(m=0;m<M;m++){
|
for(m=0;m<num_inducing;m++){
|
||||||
for(q=0;q<input_dim;q++){
|
for(q=0;q<input_dim;q++){
|
||||||
for(mm=0;mm<M;mm++){
|
for(mm=0;mm<num_inducing;mm++){
|
||||||
for(n=0;n<N;n++){
|
for(n=0;n<N;n++){
|
||||||
target(m,q) += dL_dpsi2(n,m,mm)*AZA(n,mm,q);
|
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_compile_args': ['-fopenmp -O3'], #-march=native'],
|
||||||
'extra_link_args' : ['-lgomp']}
|
'extra_link_args' : ['-lgomp']}
|
||||||
|
|
||||||
N,M,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
|
N,num_inducing,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
|
||||||
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
||||||
arg_names=['N','M','input_dim','AZA','target','dL_dpsi2'],
|
arg_names=['N','num_inducing','input_dim','AZA','target','dL_dpsi2'],
|
||||||
type_converters=weave.converters.blitz,**weave_options)
|
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))
|
muS_changed = not (np.array_equal(mu, self._mu) and np.array_equal(S, self._S))
|
||||||
if Zv_changed:
|
if Zv_changed:
|
||||||
# Z has changed, compute Z specific stuff
|
# Z has changed, compute Z specific stuff
|
||||||
# self.ZZ = Z[:,None,:]*Z[None,:,:] # M,M,input_dim
|
# self.ZZ = Z[:,None,:]*Z[None,:,:] # num_inducing,num_inducing,input_dim
|
||||||
# self.ZZ = np.empty((Z.shape[0], Z.shape[0], Z.shape[1]), order='F')
|
# 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])]
|
# [tdot(Z[:, i:i + 1], self.ZZ[:, :, i].T) for i in xrange(Z.shape[1])]
|
||||||
self.ZA = Z * self.variances
|
self.ZA = Z * self.variances
|
||||||
|
|
@ -291,5 +291,5 @@ class linear(Kernpart):
|
||||||
self.inner[:, diag_indices[0], diag_indices[1]] += S
|
self.inner[:, diag_indices[0], diag_indices[1]] += S
|
||||||
self._mu, self._S = mu.copy(), S.copy()
|
self._mu, self._S = mu.copy(), S.copy()
|
||||||
if Zv_changed or muS_changed:
|
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 input_dim]!
|
self.ZAinner = np.dot(self.ZA, self.inner).swapaxes(0, 1) # NOTE: self.ZAinner \in [num_inducing x N x input_dim]!
|
||||||
self._psi2 = np.dot(self.ZAinner, self.ZA.T)
|
self._psi2 = np.dot(self.ZAinner, self.ZA.T)
|
||||||
|
|
|
||||||
|
|
@ -113,7 +113,7 @@ class periodic_Matern32(Kernpart):
|
||||||
|
|
||||||
@silence_errors
|
@silence_errors
|
||||||
def dK_dtheta(self,dL_dK,X,X2,target):
|
def dK_dtheta(self,dL_dK,X,X2,target):
|
||||||
"""derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)"""
|
"""derivative of the covariance matrix with respect to the parameters (shape is Nxnum_inducingxNparam)"""
|
||||||
if X2 is None: X2 = X
|
if X2 is None: X2 = X
|
||||||
FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
|
FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
|
||||||
FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
|
FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
|
||||||
|
|
|
||||||
|
|
@ -115,7 +115,7 @@ class periodic_Matern52(Kernpart):
|
||||||
|
|
||||||
@silence_errors
|
@silence_errors
|
||||||
def dK_dtheta(self,dL_dK,X,X2,target):
|
def dK_dtheta(self,dL_dK,X,X2,target):
|
||||||
"""derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)"""
|
"""derivative of the covariance matrix with respect to the parameters (shape is Nxnum_inducingxNparam)"""
|
||||||
if X2 is None: X2 = X
|
if X2 is None: X2 = X
|
||||||
FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
|
FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
|
||||||
FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
|
FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
|
||||||
|
|
|
||||||
|
|
@ -111,7 +111,7 @@ class periodic_exponential(Kernpart):
|
||||||
|
|
||||||
@silence_errors
|
@silence_errors
|
||||||
def dK_dtheta(self,dL_dK,X,X2,target):
|
def dK_dtheta(self,dL_dK,X,X2,target):
|
||||||
"""derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)"""
|
"""derivative of the covariance matrix with respect to the parameters (shape is Nxnum_inducingxNparam)"""
|
||||||
if X2 is None: X2 = X
|
if X2 is None: X2 = X
|
||||||
FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
|
FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
|
||||||
FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
|
FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
|
||||||
|
|
|
||||||
|
|
@ -110,7 +110,7 @@ class rbf(Kernpart):
|
||||||
target(q+1) += var_len3(q)*tmp;
|
target(q+1) += var_len3(q)*tmp;
|
||||||
}
|
}
|
||||||
"""
|
"""
|
||||||
N, M, input_dim = X.shape[0], X.shape[0], self.input_dim
|
N, num_inducing, input_dim = X.shape[0], X.shape[0], self.input_dim
|
||||||
else:
|
else:
|
||||||
code = """
|
code = """
|
||||||
int q,i,j;
|
int q,i,j;
|
||||||
|
|
@ -118,16 +118,16 @@ class rbf(Kernpart):
|
||||||
for(q=0; q<input_dim; q++){
|
for(q=0; q<input_dim; q++){
|
||||||
tmp = 0;
|
tmp = 0;
|
||||||
for(i=0; i<N; i++){
|
for(i=0; i<N; i++){
|
||||||
for(j=0; j<M; j++){
|
for(j=0; j<num_inducing; j++){
|
||||||
tmp += (X(i,q)-X2(j,q))*(X(i,q)-X2(j,q))*dvardLdK(i,j);
|
tmp += (X(i,q)-X2(j,q))*(X(i,q)-X2(j,q))*dvardLdK(i,j);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
target(q+1) += var_len3(q)*tmp;
|
target(q+1) += var_len3(q)*tmp;
|
||||||
}
|
}
|
||||||
"""
|
"""
|
||||||
N, M, input_dim = X.shape[0], X2.shape[0], self.input_dim
|
N, num_inducing, input_dim = X.shape[0], X2.shape[0], self.input_dim
|
||||||
# [np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.input_dim)]
|
# [np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.input_dim)]
|
||||||
weave.inline(code, arg_names=['N', 'M', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3'],
|
weave.inline(code, arg_names=['N','num_inducing','input_dim','X','X2','target','dvardLdK','var_len3'],
|
||||||
type_converters=weave.converters.blitz, **self.weave_options)
|
type_converters=weave.converters.blitz, **self.weave_options)
|
||||||
else:
|
else:
|
||||||
target[1] += (self.variance / self.lengthscale) * np.sum(self._K_dvar * self._K_dist2 * dL_dK)
|
target[1] += (self.variance / self.lengthscale) * np.sum(self._K_dvar * self._K_dist2 * dL_dK)
|
||||||
|
|
@ -191,7 +191,7 @@ class rbf(Kernpart):
|
||||||
target += self._psi2
|
target += self._psi2
|
||||||
|
|
||||||
def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, target):
|
def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, target):
|
||||||
"""Shape N,M,M,Ntheta"""
|
"""Shape N,num_inducing,num_inducing,Ntheta"""
|
||||||
self._psi_computations(Z, mu, S)
|
self._psi_computations(Z, mu, S)
|
||||||
d_var = 2.*self._psi2 / self.variance
|
d_var = 2.*self._psi2 / self.variance
|
||||||
d_length = 2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] / self.lengthscale2) / (self.lengthscale * self._psi2_denom)
|
d_length = 2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] / self.lengthscale2) / (self.lengthscale * self._psi2_denom)
|
||||||
|
|
@ -205,19 +205,18 @@ class rbf(Kernpart):
|
||||||
|
|
||||||
def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, target):
|
def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, target):
|
||||||
self._psi_computations(Z, mu, S)
|
self._psi_computations(Z, mu, S)
|
||||||
term1 = self._psi2_Zdist / self.lengthscale2 # M, M, input_dim
|
term1 = self._psi2_Zdist / self.lengthscale2 # num_inducing, num_inducing, input_dim
|
||||||
term2 = self._psi2_mudist / self._psi2_denom / self.lengthscale2 # N, M, M, input_dim
|
term2 = self._psi2_mudist / self._psi2_denom / self.lengthscale2 # N, num_inducing, num_inducing, input_dim
|
||||||
dZ = self._psi2[:, :, :, None] * (term1[None] + term2)
|
dZ = self._psi2[:, :, :, None] * (term1[None] + term2)
|
||||||
target += (dL_dpsi2[:, :, :, None] * dZ).sum(0).sum(0)
|
target += (dL_dpsi2[:, :, :, None] * dZ).sum(0).sum(0)
|
||||||
|
|
||||||
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
|
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
|
||||||
"""Think N,M,M,input_dim """
|
"""Think N,num_inducing,num_inducing,input_dim """
|
||||||
self._psi_computations(Z, mu, S)
|
self._psi_computations(Z, mu, S)
|
||||||
tmp = self._psi2[:, :, :, None] / self.lengthscale2 / self._psi2_denom
|
tmp = self._psi2[:, :, :, None] / self.lengthscale2 / self._psi2_denom
|
||||||
target_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * self._psi2_mudist).sum(1).sum(1)
|
target_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * self._psi2_mudist).sum(1).sum(1)
|
||||||
target_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*self._psi2_mudist_sq - 1)).sum(1).sum(1)
|
target_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*self._psi2_mudist_sq - 1)).sum(1).sum(1)
|
||||||
|
|
||||||
|
|
||||||
#---------------------------------------#
|
#---------------------------------------#
|
||||||
# Precomputations #
|
# Precomputations #
|
||||||
#---------------------------------------#
|
#---------------------------------------#
|
||||||
|
|
@ -241,41 +240,41 @@ class rbf(Kernpart):
|
||||||
def _psi_computations(self, Z, mu, S):
|
def _psi_computations(self, Z, mu, S):
|
||||||
# here are the "statistics" for psi1 and psi2
|
# here are the "statistics" for psi1 and psi2
|
||||||
if not np.array_equal(Z, self._Z):
|
if not np.array_equal(Z, self._Z):
|
||||||
# Z has changed, compute Z specific stuff
|
#Z has changed, compute Z specific stuff
|
||||||
self._psi2_Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,input_dim
|
self._psi2_Zhat = 0.5*(Z[:,None,:] +Z[None,:,:]) # num_inducing,num_inducing,input_dim
|
||||||
self._psi2_Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,input_dim
|
self._psi2_Zdist = 0.5*(Z[:,None,:]-Z[None,:,:]) # num_inducing,num_inducing,input_dim
|
||||||
self._psi2_Zdist_sq = np.square(self._psi2_Zdist / self.lengthscale) # M,M,input_dim
|
self._psi2_Zdist_sq = np.square(self._psi2_Zdist/self.lengthscale) # num_inducing,num_inducing,input_dim
|
||||||
self._Z = Z
|
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)):
|
if not (np.array_equal(Z, self._Z) and np.array_equal(mu, self._mu) and np.array_equal(S, self._S)):
|
||||||
# something's changed. recompute EVERYTHING
|
#something's changed. recompute EVERYTHING
|
||||||
|
|
||||||
# psi1
|
#psi1
|
||||||
self._psi1_denom = S[:, None, :] / self.lengthscale2 + 1.
|
self._psi1_denom = S[:,None,:]/self.lengthscale2 + 1.
|
||||||
self._psi1_dist = Z[None, :, :] - mu[:, None, :]
|
self._psi1_dist = Z[None,:,:]-mu[:,None,:]
|
||||||
self._psi1_dist_sq = np.square(self._psi1_dist) / self.lengthscale2 / self._psi1_denom
|
self._psi1_dist_sq = np.square(self._psi1_dist)/self.lengthscale2/self._psi1_denom
|
||||||
self._psi1_exponent = -0.5 * np.sum(self._psi1_dist_sq + np.log(self._psi1_denom), -1)
|
self._psi1_exponent = -0.5*np.sum(self._psi1_dist_sq+np.log(self._psi1_denom),-1)
|
||||||
self._psi1 = self.variance * np.exp(self._psi1_exponent)
|
self._psi1 = self.variance*np.exp(self._psi1_exponent)
|
||||||
|
|
||||||
# psi2
|
#psi2
|
||||||
self._psi2_denom = 2.*S[:, None, None, :] / self.lengthscale2 + 1. # N,M,M,input_dim
|
self._psi2_denom = 2.*S[:,None,None,:]/self.lengthscale2+1. # N,num_inducing,num_inducing,input_dim
|
||||||
self._psi2_mudist, self._psi2_mudist_sq, self._psi2_exponent, _ = self.weave_psi2(mu, self._psi2_Zhat)
|
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,input_dim
|
#self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,num_inducing,num_inducing,input_dim
|
||||||
# self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscale2*self._psi2_denom)
|
#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_exponent = np.sum(-self._psi2_Zdist_sq -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,num_inducing,num_inducing
|
||||||
self._psi2 = np.square(self.variance) * np.exp(self._psi2_exponent) # N,M,M
|
self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,num_inducing,num_inducing
|
||||||
|
|
||||||
# store matrices for caching
|
#store matrices for caching
|
||||||
self._Z, self._mu, self._S = Z, mu, S
|
self._Z, self._mu, self._S = Z, mu,S
|
||||||
|
|
||||||
def weave_psi2(self, mu, Zhat):
|
def weave_psi2(self,mu,Zhat):
|
||||||
N, input_dim = mu.shape
|
N,input_dim = mu.shape
|
||||||
M = Zhat.shape[0]
|
num_inducing = Zhat.shape[0]
|
||||||
|
|
||||||
mudist = np.empty((N, M, M, input_dim))
|
mudist = np.empty((N,num_inducing,num_inducing,input_dim))
|
||||||
mudist_sq = np.empty((N, M, M, input_dim))
|
mudist_sq = np.empty((N,num_inducing,num_inducing,input_dim))
|
||||||
psi2_exponent = np.zeros((N, M, M))
|
psi2_exponent = np.zeros((N,num_inducing,num_inducing))
|
||||||
psi2 = np.empty((N, M, M))
|
psi2 = np.empty((N,num_inducing,num_inducing))
|
||||||
|
|
||||||
psi2_Zdist_sq = self._psi2_Zdist_sq
|
psi2_Zdist_sq = self._psi2_Zdist_sq
|
||||||
_psi2_denom = self._psi2_denom.squeeze().reshape(N, self.input_dim)
|
_psi2_denom = self._psi2_denom.squeeze().reshape(N, self.input_dim)
|
||||||
|
|
@ -290,7 +289,7 @@ class rbf(Kernpart):
|
||||||
|
|
||||||
#pragma omp parallel for private(tmp)
|
#pragma omp parallel for private(tmp)
|
||||||
for (int n=0; n<N; n++){
|
for (int n=0; n<N; n++){
|
||||||
for (int m=0; m<M; m++){
|
for (int m=0; m<num_inducing; m++){
|
||||||
for (int mm=0; mm<(m+1); mm++){
|
for (int mm=0; mm<(m+1); mm++){
|
||||||
for (int q=0; q<input_dim; q++){
|
for (int q=0; q<input_dim; q++){
|
||||||
//compute mudist
|
//compute mudist
|
||||||
|
|
@ -325,7 +324,7 @@ class rbf(Kernpart):
|
||||||
#include <math.h>
|
#include <math.h>
|
||||||
"""
|
"""
|
||||||
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
||||||
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'],
|
arg_names=['N','num_inducing','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)
|
type_converters=weave.converters.blitz, **self.weave_options)
|
||||||
|
|
||||||
return mudist, mudist_sq, psi2_exponent, psi2
|
return mudist, mudist_sq, psi2_exponent, psi2
|
||||||
|
|
|
||||||
|
|
@ -122,12 +122,12 @@ class spkern(Kernpart):
|
||||||
int i;
|
int i;
|
||||||
int j;
|
int j;
|
||||||
int N = target_array->dimensions[0];
|
int N = target_array->dimensions[0];
|
||||||
int M = target_array->dimensions[1];
|
int num_inducing = target_array->dimensions[1];
|
||||||
int input_dim = X_array->dimensions[1];
|
int input_dim = X_array->dimensions[1];
|
||||||
//#pragma omp parallel for private(j)
|
//#pragma omp parallel for private(j)
|
||||||
for (i=0;i<N;i++){
|
for (i=0;i<N;i++){
|
||||||
for (j=0;j<M;j++){
|
for (j=0;j<num_inducing;j++){
|
||||||
target[i*M+j] = k(%s);
|
target[i*num_inducing+j] = k(%s);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
%s
|
%s
|
||||||
|
|
@ -149,17 +149,17 @@ class spkern(Kernpart):
|
||||||
"""%(diag_arglist,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed
|
"""%(diag_arglist,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed
|
||||||
|
|
||||||
#here's some code to compute gradients
|
#here's some code to compute gradients
|
||||||
funclist = '\n'.join([' '*16 + 'target[%i] += partial[i*M+j]*dk_d%s(%s);'%(i,theta.name,arglist) for i,theta in enumerate(self._sp_theta)])
|
funclist = '\n'.join([' '*16 + 'target[%i] += partial[i*num_inducing+j]*dk_d%s(%s);'%(i,theta.name,arglist) for i,theta in enumerate(self._sp_theta)])
|
||||||
self._dK_dtheta_code =\
|
self._dK_dtheta_code =\
|
||||||
"""
|
"""
|
||||||
int i;
|
int i;
|
||||||
int j;
|
int j;
|
||||||
int N = partial_array->dimensions[0];
|
int N = partial_array->dimensions[0];
|
||||||
int M = partial_array->dimensions[1];
|
int num_inducing = partial_array->dimensions[1];
|
||||||
int input_dim = X_array->dimensions[1];
|
int input_dim = X_array->dimensions[1];
|
||||||
//#pragma omp parallel for private(j)
|
//#pragma omp parallel for private(j)
|
||||||
for (i=0;i<N;i++){
|
for (i=0;i<N;i++){
|
||||||
for (j=0;j<M;j++){
|
for (j=0;j<num_inducing;j++){
|
||||||
%s
|
%s
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
@ -169,7 +169,7 @@ class spkern(Kernpart):
|
||||||
#here's some code to compute gradients for Kdiag TODO: thius is yucky.
|
#here's some code to compute gradients for Kdiag TODO: thius is yucky.
|
||||||
diag_funclist = re.sub('Z','X',funclist,count=0)
|
diag_funclist = re.sub('Z','X',funclist,count=0)
|
||||||
diag_funclist = re.sub('j','i',diag_funclist)
|
diag_funclist = re.sub('j','i',diag_funclist)
|
||||||
diag_funclist = re.sub('partial\[i\*M\+i\]','partial[i]',diag_funclist)
|
diag_funclist = re.sub('partial\[i\*num_inducing\+i\]','partial[i]',diag_funclist)
|
||||||
self._dKdiag_dtheta_code =\
|
self._dKdiag_dtheta_code =\
|
||||||
"""
|
"""
|
||||||
int i;
|
int i;
|
||||||
|
|
@ -182,17 +182,17 @@ class spkern(Kernpart):
|
||||||
"""%(diag_funclist,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed
|
"""%(diag_funclist,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed
|
||||||
|
|
||||||
#Here's some code to do gradients wrt x
|
#Here's some code to do gradients wrt x
|
||||||
gradient_funcs = "\n".join(["target[i*input_dim+%i] += partial[i*M+j]*dk_dx%i(%s);"%(q,q,arglist) for q in range(self.input_dim)])
|
gradient_funcs = "\n".join(["target[i*input_dim+%i] += partial[i*num_inducing+j]*dk_dx%i(%s);"%(q,q,arglist) for q in range(self.input_dim)])
|
||||||
self._dK_dX_code = \
|
self._dK_dX_code = \
|
||||||
"""
|
"""
|
||||||
int i;
|
int i;
|
||||||
int j;
|
int j;
|
||||||
int N = partial_array->dimensions[0];
|
int N = partial_array->dimensions[0];
|
||||||
int M = partial_array->dimensions[1];
|
int num_inducing = partial_array->dimensions[1];
|
||||||
int input_dim = X_array->dimensions[1];
|
int input_dim = X_array->dimensions[1];
|
||||||
//#pragma omp parallel for private(j)
|
//#pragma omp parallel for private(j)
|
||||||
for (i=0;i<N; i++){
|
for (i=0;i<N; i++){
|
||||||
for (j=0; j<M; j++){
|
for (j=0; j<num_inducing; j++){
|
||||||
%s
|
%s
|
||||||
//if(isnan(target[i*input_dim+2])){printf("%%f\\n",dk_dx2(X[i*input_dim+0], X[i*input_dim+1], X[i*input_dim+2], Z[j*input_dim+0], Z[j*input_dim+1], Z[j*input_dim+2], param[0], param[1], param[2], param[3], param[4], param[5]));}
|
//if(isnan(target[i*input_dim+2])){printf("%%f\\n",dk_dx2(X[i*input_dim+0], X[i*input_dim+1], X[i*input_dim+2], Z[j*input_dim+0], Z[j*input_dim+1], Z[j*input_dim+2], param[0], param[1], param[2], param[3], param[4], param[5]));}
|
||||||
//if(isnan(target[i*input_dim+2])){printf("%%f,%%f,%%i,%%i\\n", X[i*input_dim+2], Z[j*input_dim+2],i,j);}
|
//if(isnan(target[i*input_dim+2])){printf("%%f,%%f,%%i,%%i\\n", X[i*input_dim+2], Z[j*input_dim+2],i,j);}
|
||||||
|
|
@ -208,7 +208,7 @@ class spkern(Kernpart):
|
||||||
int i;
|
int i;
|
||||||
int j;
|
int j;
|
||||||
int N = partial_array->dimensions[0];
|
int N = partial_array->dimensions[0];
|
||||||
int M = 0;
|
int num_inducing = 0;
|
||||||
int input_dim = X_array->dimensions[1];
|
int input_dim = X_array->dimensions[1];
|
||||||
for (i=0;i<N; i++){
|
for (i=0;i<N; i++){
|
||||||
j = i;
|
j = i;
|
||||||
|
|
|
||||||
|
|
@ -139,7 +139,7 @@ class EP(likelihood):
|
||||||
The expectation-propagation algorithm with sparse pseudo-input.
|
The expectation-propagation algorithm with sparse pseudo-input.
|
||||||
For nomenclature see ... 2013.
|
For nomenclature see ... 2013.
|
||||||
"""
|
"""
|
||||||
M = Kmm.shape[0]
|
num_inducing = Kmm.shape[0]
|
||||||
|
|
||||||
#TODO: this doesn't work with uncertain inputs!
|
#TODO: this doesn't work with uncertain inputs!
|
||||||
|
|
||||||
|
|
@ -235,7 +235,7 @@ class EP(likelihood):
|
||||||
The expectation-propagation algorithm with sparse pseudo-input.
|
The expectation-propagation algorithm with sparse pseudo-input.
|
||||||
For nomenclature see Naish-Guzman and Holden, 2008.
|
For nomenclature see Naish-Guzman and Holden, 2008.
|
||||||
"""
|
"""
|
||||||
M = Kmm.shape[0]
|
num_inducing = Kmm.shape[0]
|
||||||
|
|
||||||
"""
|
"""
|
||||||
Prior approximation parameters:
|
Prior approximation parameters:
|
||||||
|
|
@ -258,7 +258,7 @@ class EP(likelihood):
|
||||||
mu = w + P*Gamma
|
mu = w + P*Gamma
|
||||||
"""
|
"""
|
||||||
self.w = np.zeros(self.N)
|
self.w = np.zeros(self.N)
|
||||||
self.Gamma = np.zeros(M)
|
self.Gamma = np.zeros(num_inducing)
|
||||||
mu = np.zeros(self.N)
|
mu = np.zeros(self.N)
|
||||||
P = P0.copy()
|
P = P0.copy()
|
||||||
R = R0.copy()
|
R = R0.copy()
|
||||||
|
|
@ -305,10 +305,10 @@ class EP(likelihood):
|
||||||
dtd1 = Delta_tau*Diag[i] + 1.
|
dtd1 = Delta_tau*Diag[i] + 1.
|
||||||
dii = Diag[i]
|
dii = Diag[i]
|
||||||
Diag[i] = dii - (Delta_tau * dii**2.)/dtd1
|
Diag[i] = dii - (Delta_tau * dii**2.)/dtd1
|
||||||
pi_ = P[i,:].reshape(1,M)
|
pi_ = P[i,:].reshape(1,num_inducing)
|
||||||
P[i,:] = pi_ - (Delta_tau*dii)/dtd1 * pi_
|
P[i,:] = pi_ - (Delta_tau*dii)/dtd1 * pi_
|
||||||
Rp_i = np.dot(R,pi_.T)
|
Rp_i = np.dot(R,pi_.T)
|
||||||
RTR = np.dot(R.T,np.dot(np.eye(M) - Delta_tau/(1.+Delta_tau*Sigma_diag[i]) * np.dot(Rp_i,Rp_i.T),R))
|
RTR = np.dot(R.T,np.dot(np.eye(num_inducing) - Delta_tau/(1.+Delta_tau*Sigma_diag[i]) * np.dot(Rp_i,Rp_i.T),R))
|
||||||
R = jitchol(RTR).T
|
R = jitchol(RTR).T
|
||||||
self.w[i] += (Delta_v - Delta_tau*self.w[i])*dii/dtd1
|
self.w[i] += (Delta_v - Delta_tau*self.w[i])*dii/dtd1
|
||||||
self.Gamma += (Delta_v - Delta_tau*mu[i])*np.dot(RTR,P[i,:].T)
|
self.Gamma += (Delta_v - Delta_tau*mu[i])*np.dot(RTR,P[i,:].T)
|
||||||
|
|
@ -321,7 +321,7 @@ class EP(likelihood):
|
||||||
Diag = Diag0 * Iplus_Dprod_i
|
Diag = Diag0 * Iplus_Dprod_i
|
||||||
P = Iplus_Dprod_i[:,None] * P0
|
P = Iplus_Dprod_i[:,None] * P0
|
||||||
safe_diag = np.where(Diag0 < self.tau_tilde, self.tau_tilde/(1.+Diag0*self.tau_tilde), (1. - Iplus_Dprod_i)/Diag0)
|
safe_diag = np.where(Diag0 < self.tau_tilde, self.tau_tilde/(1.+Diag0*self.tau_tilde), (1. - Iplus_Dprod_i)/Diag0)
|
||||||
L = jitchol(np.eye(M) + np.dot(RPT0,safe_diag[:,None]*RPT0.T))
|
L = jitchol(np.eye(num_inducing) + np.dot(RPT0,safe_diag[:,None]*RPT0.T))
|
||||||
R,info = linalg.lapack.flapack.dtrtrs(L,R0,lower=1)
|
R,info = linalg.lapack.flapack.dtrtrs(L,R0,lower=1)
|
||||||
RPT = np.dot(R,P.T)
|
RPT = np.dot(R,P.T)
|
||||||
Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1)
|
Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1)
|
||||||
|
|
|
||||||
|
|
@ -23,7 +23,7 @@ class BayesianGPLVM(SparseGP, GPLVM):
|
||||||
:type init: 'PCA'|'random'
|
:type init: 'PCA'|'random'
|
||||||
|
|
||||||
"""
|
"""
|
||||||
def __init__(self, likelihood_or_Y, input_dim, X=None, X_variance=None, init='PCA', M=10,
|
def __init__(self, likelihood_or_Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10,
|
||||||
Z=None, kernel=None, oldpsave=10, _debug=False,
|
Z=None, kernel=None, oldpsave=10, _debug=False,
|
||||||
**kwargs):
|
**kwargs):
|
||||||
if type(likelihood_or_Y) is np.ndarray:
|
if type(likelihood_or_Y) is np.ndarray:
|
||||||
|
|
@ -39,7 +39,7 @@ class BayesianGPLVM(SparseGP, GPLVM):
|
||||||
X_variance = np.clip((np.ones_like(X) * 0.5) + .01 * np.random.randn(*X.shape), 0.001, 1)
|
X_variance = np.clip((np.ones_like(X) * 0.5) + .01 * np.random.randn(*X.shape), 0.001, 1)
|
||||||
|
|
||||||
if Z is None:
|
if Z is None:
|
||||||
Z = np.random.permutation(X.copy())[:M]
|
Z = np.random.permutation(X.copy())[:num_inducing]
|
||||||
assert Z.shape[1] == X.shape[1]
|
assert Z.shape[1] == X.shape[1]
|
||||||
|
|
||||||
if kernel is None:
|
if kernel is None:
|
||||||
|
|
|
||||||
|
|
@ -44,7 +44,7 @@ class MRD(Model):
|
||||||
:param kernels: list of kernels or kernel shared for all BGPLVMS
|
:param kernels: list of kernels or kernel shared for all BGPLVMS
|
||||||
:type kernels: [GPy.kern.kern] | GPy.kern.kern | None (default)
|
:type kernels: [GPy.kern.kern] | GPy.kern.kern | None (default)
|
||||||
"""
|
"""
|
||||||
def __init__(self, likelihood_or_Y_list, input_dim, M=10, names=None,
|
def __init__(self, likelihood_or_Y_list, input_dim, num_inducing=10, names=None,
|
||||||
kernels=None, initx='PCA',
|
kernels=None, initx='PCA',
|
||||||
initz='permute', _debug=False, **kw):
|
initz='permute', _debug=False, **kw):
|
||||||
if names is None:
|
if names is None:
|
||||||
|
|
@ -61,13 +61,13 @@ class MRD(Model):
|
||||||
assert not ('kernel' in kw), "pass kernels through `kernels` argument"
|
assert not ('kernel' in kw), "pass kernels through `kernels` argument"
|
||||||
|
|
||||||
self.input_dim = input_dim
|
self.input_dim = input_dim
|
||||||
self.num_inducing = M
|
self.num_inducing = num_inducing
|
||||||
self._debug = _debug
|
self._debug = _debug
|
||||||
|
|
||||||
self._init = True
|
self._init = True
|
||||||
X = self._init_X(initx, likelihood_or_Y_list)
|
X = self._init_X(initx, likelihood_or_Y_list)
|
||||||
Z = self._init_Z(initz, X)
|
Z = self._init_Z(initz, X)
|
||||||
self.bgplvms = [BayesianGPLVM(l, input_dim=input_dim, kernel=k, X=X, Z=Z, M=self.num_inducing, **kw) for l, k in zip(likelihood_or_Y_list, kernels)]
|
self.bgplvms = [BayesianGPLVM(l, input_dim=input_dim, kernel=k, X=X, Z=Z, num_inducing=self.num_inducing, **kw) for l, k in zip(likelihood_or_Y_list, kernels)]
|
||||||
del self._init
|
del self._init
|
||||||
|
|
||||||
self.gref = self.bgplvms[0]
|
self.gref = self.bgplvms[0]
|
||||||
|
|
|
||||||
|
|
@ -26,7 +26,7 @@ class SparseGPClassification(SparseGP):
|
||||||
|
|
||||||
"""
|
"""
|
||||||
|
|
||||||
def __init__(self, X, Y=None, likelihood=None, kernel=None, normalize_X=False, normalize_Y=False, Z=None, M=10):
|
def __init__(self, X, Y=None, likelihood=None, kernel=None, normalize_X=False, normalize_Y=False, Z=None, num_inducing=10):
|
||||||
if kernel is None:
|
if kernel is None:
|
||||||
kernel = kern.rbf(X.shape[1]) + kern.white(X.shape[1],1e-3)
|
kernel = kern.rbf(X.shape[1]) + kern.white(X.shape[1],1e-3)
|
||||||
|
|
||||||
|
|
@ -38,7 +38,7 @@ class SparseGPClassification(SparseGP):
|
||||||
raise Warning, 'likelihood.data and Y are different.'
|
raise Warning, 'likelihood.data and Y are different.'
|
||||||
|
|
||||||
if Z is None:
|
if Z is None:
|
||||||
i = np.random.permutation(X.shape[0])[:M]
|
i = np.random.permutation(X.shape[0])[:num_inducing]
|
||||||
Z = X[i].copy()
|
Z = X[i].copy()
|
||||||
else:
|
else:
|
||||||
assert Z.shape[1]==X.shape[1]
|
assert Z.shape[1]==X.shape[1]
|
||||||
|
|
|
||||||
|
|
@ -26,14 +26,14 @@ class SparseGPRegression(SparseGP):
|
||||||
|
|
||||||
"""
|
"""
|
||||||
|
|
||||||
def __init__(self, X, Y, kernel=None, normalize_X=False, normalize_Y=False, Z=None, M=10, X_variance=None):
|
def __init__(self, X, Y, kernel=None, normalize_X=False, normalize_Y=False, Z=None, num_inducing=10, X_variance=None):
|
||||||
# kern defaults to rbf (plus white for stability)
|
# kern defaults to rbf (plus white for stability)
|
||||||
if kernel is None:
|
if kernel is None:
|
||||||
kernel = kern.rbf(X.shape[1]) + kern.white(X.shape[1], 1e-3)
|
kernel = kern.rbf(X.shape[1]) + kern.white(X.shape[1], 1e-3)
|
||||||
|
|
||||||
# Z defaults to a subset of the data
|
# Z defaults to a subset of the data
|
||||||
if Z is None:
|
if Z is None:
|
||||||
i = np.random.permutation(X.shape[0])[:M]
|
i = np.random.permutation(X.shape[0])[:num_inducing]
|
||||||
Z = X[i].copy()
|
Z = X[i].copy()
|
||||||
else:
|
else:
|
||||||
assert Z.shape[1] == X.shape[1]
|
assert Z.shape[1] == X.shape[1]
|
||||||
|
|
|
||||||
|
|
@ -23,9 +23,9 @@ class SparseGPLVM(SparseGPRegression, GPLVM):
|
||||||
:type init: 'PCA'|'random'
|
:type init: 'PCA'|'random'
|
||||||
|
|
||||||
"""
|
"""
|
||||||
def __init__(self, Y, input_dim, kernel=None, init='PCA', M=10):
|
def __init__(self, Y, input_dim, kernel=None, init='PCA', num_inducing=10):
|
||||||
X = self.initialise_latent(init, input_dim, Y)
|
X = self.initialise_latent(init, input_dim, Y)
|
||||||
SparseGPRegression.__init__(self, X, Y, kernel=kernel, M=M)
|
SparseGPRegression.__init__(self, X, Y, kernel=kernel, num_inducing=num_inducing)
|
||||||
|
|
||||||
def _get_param_names(self):
|
def _get_param_names(self):
|
||||||
return (sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] 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)], [])
|
||||||
|
|
|
||||||
|
|
@ -8,67 +8,67 @@ from GPy.models.bayesian_gplvm import BayesianGPLVM
|
||||||
|
|
||||||
class BGPLVMTests(unittest.TestCase):
|
class BGPLVMTests(unittest.TestCase):
|
||||||
def test_bias_kern(self):
|
def test_bias_kern(self):
|
||||||
N, M, input_dim, D = 10, 3, 2, 4
|
N, num_inducing, input_dim, D = 10, 3, 2, 4
|
||||||
X = np.random.rand(N, input_dim)
|
X = np.random.rand(N, input_dim)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
||||||
Y -= Y.mean(axis=0)
|
Y -= Y.mean(axis=0)
|
||||||
k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
m = BayesianGPLVM(Y, input_dim, kernel=k, M=M)
|
m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
m.randomize()
|
m.randomize()
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
||||||
def test_linear_kern(self):
|
def test_linear_kern(self):
|
||||||
N, M, input_dim, D = 10, 3, 2, 4
|
N, num_inducing, input_dim, D = 10, 3, 2, 4
|
||||||
X = np.random.rand(N, input_dim)
|
X = np.random.rand(N, input_dim)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
||||||
Y -= Y.mean(axis=0)
|
Y -= Y.mean(axis=0)
|
||||||
k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
m = BayesianGPLVM(Y, input_dim, kernel=k, M=M)
|
m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
m.randomize()
|
m.randomize()
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
||||||
def test_rbf_kern(self):
|
def test_rbf_kern(self):
|
||||||
N, M, input_dim, D = 10, 3, 2, 4
|
N, num_inducing, input_dim, D = 10, 3, 2, 4
|
||||||
X = np.random.rand(N, input_dim)
|
X = np.random.rand(N, input_dim)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
||||||
Y -= Y.mean(axis=0)
|
Y -= Y.mean(axis=0)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
m = BayesianGPLVM(Y, input_dim, kernel=k, M=M)
|
m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
m.randomize()
|
m.randomize()
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
||||||
def test_rbf_bias_kern(self):
|
def test_rbf_bias_kern(self):
|
||||||
N, M, input_dim, D = 10, 3, 2, 4
|
N, num_inducing, input_dim, D = 10, 3, 2, 4
|
||||||
X = np.random.rand(N, input_dim)
|
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 = GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
||||||
Y -= Y.mean(axis=0)
|
Y -= Y.mean(axis=0)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
m = BayesianGPLVM(Y, input_dim, kernel=k, M=M)
|
m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
m.randomize()
|
m.randomize()
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
||||||
#@unittest.skip('psi2 cross terms are NotImplemented for this combination')
|
#@unittest.skip('psi2 cross terms are NotImplemented for this combination')
|
||||||
def test_linear_bias_kern(self):
|
def test_linear_bias_kern(self):
|
||||||
N, M, input_dim, D = 30, 5, 4, 30
|
N, num_inducing, input_dim, D = 30, 5, 4, 30
|
||||||
X = np.random.rand(N, input_dim)
|
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 = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
||||||
Y -= Y.mean(axis=0)
|
Y -= Y.mean(axis=0)
|
||||||
k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
m = BayesianGPLVM(Y, input_dim, kernel=k, M=M)
|
m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
m.randomize()
|
m.randomize()
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
|
||||||
|
|
@ -7,7 +7,7 @@ import GPy
|
||||||
|
|
||||||
class GPLVMTests(unittest.TestCase):
|
class GPLVMTests(unittest.TestCase):
|
||||||
def test_bias_kern(self):
|
def test_bias_kern(self):
|
||||||
N, M, input_dim, D = 10, 3, 2, 4
|
N, num_inducing, input_dim, D = 10, 3, 2, 4
|
||||||
X = np.random.rand(N, input_dim)
|
X = np.random.rand(N, input_dim)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
|
|
@ -19,7 +19,7 @@ class GPLVMTests(unittest.TestCase):
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
||||||
def test_linear_kern(self):
|
def test_linear_kern(self):
|
||||||
N, M, input_dim, D = 10, 3, 2, 4
|
N, num_inducing, input_dim, D = 10, 3, 2, 4
|
||||||
X = np.random.rand(N, input_dim)
|
X = np.random.rand(N, input_dim)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
|
|
@ -31,7 +31,7 @@ class GPLVMTests(unittest.TestCase):
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
||||||
def test_rbf_kern(self):
|
def test_rbf_kern(self):
|
||||||
N, M, input_dim, D = 10, 3, 2, 4
|
N, num_inducing, input_dim, D = 10, 3, 2, 4
|
||||||
X = np.random.rand(N, input_dim)
|
X = np.random.rand(N, input_dim)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
|
|
|
||||||
|
|
@ -14,7 +14,7 @@ class MRDTests(unittest.TestCase):
|
||||||
|
|
||||||
def test_gradients(self):
|
def test_gradients(self):
|
||||||
num_m = 3
|
num_m = 3
|
||||||
N, M, input_dim, D = 20, 8, 6, 20
|
N, num_inducing, input_dim, D = 20, 8, 6, 20
|
||||||
X = np.random.rand(N, input_dim)
|
X = np.random.rand(N, input_dim)
|
||||||
|
|
||||||
k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim)
|
k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim)
|
||||||
|
|
@ -23,7 +23,7 @@ class MRDTests(unittest.TestCase):
|
||||||
Ylist = [np.random.multivariate_normal(np.zeros(N), K, input_dim).T for _ in range(num_m)]
|
Ylist = [np.random.multivariate_normal(np.zeros(N), K, input_dim).T for _ in range(num_m)]
|
||||||
likelihood_list = [GPy.likelihoods.Gaussian(Y) for Y in Ylist]
|
likelihood_list = [GPy.likelihoods.Gaussian(Y) for Y in Ylist]
|
||||||
|
|
||||||
m = GPy.models.MRD(likelihood_list, input_dim=input_dim, kernels=k, M=M)
|
m = GPy.models.MRD(likelihood_list, input_dim=input_dim, kernels=k, num_inducing=num_inducing)
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
|
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
|
||||||
|
|
@ -11,7 +11,7 @@ import itertools
|
||||||
from GPy.core import Model
|
from GPy.core import Model
|
||||||
|
|
||||||
class PsiStatModel(Model):
|
class PsiStatModel(Model):
|
||||||
def __init__(self, which, X, X_variance, Z, M, kernel):
|
def __init__(self, which, X, X_variance, Z, num_inducing, kernel):
|
||||||
self.which = which
|
self.which = which
|
||||||
self.X = X
|
self.X = X
|
||||||
self.X_variance = X_variance
|
self.X_variance = X_variance
|
||||||
|
|
@ -64,8 +64,8 @@ class DPsiStatTest(unittest.TestCase):
|
||||||
|
|
||||||
def testPsi0(self):
|
def testPsi0(self):
|
||||||
for k in self.kernels:
|
for k in self.kernels:
|
||||||
m = PsiStatModel('psi0', X=self.X, X_variance=self.X_var, Z=self.Z,
|
m = PsiStatModel('psi1', X=self.X, X_variance=self.X_var, Z=self.Z,
|
||||||
M=self.num_inducing, kernel=k)
|
num_inducing=self.num_inducing, kernel=k)
|
||||||
try:
|
try:
|
||||||
assert m.checkgrad(), "{} x psi0".format("+".join(map(lambda x: x.name, k.parts)))
|
assert m.checkgrad(), "{} x psi0".format("+".join(map(lambda x: x.name, k.parts)))
|
||||||
except:
|
except:
|
||||||
|
|
@ -74,33 +74,33 @@ class DPsiStatTest(unittest.TestCase):
|
||||||
# def testPsi1(self):
|
# def testPsi1(self):
|
||||||
# for k in self.kernels:
|
# for k in self.kernels:
|
||||||
# m = PsiStatModel('psi1', X=self.X, X_variance=self.X_var, Z=self.Z,
|
# m = PsiStatModel('psi1', X=self.X, X_variance=self.X_var, Z=self.Z,
|
||||||
# M=self.M, kernel=k)
|
# num_inducing=self.num_inducing, kernel=k)
|
||||||
# assert m.checkgrad(), "{} x psi1".format("+".join(map(lambda x: x.name, k.parts)))
|
# assert m.checkgrad(), "{} x psi1".format("+".join(map(lambda x: x.name, k.parts)))
|
||||||
|
|
||||||
def testPsi2_lin(self):
|
def testPsi2_lin(self):
|
||||||
k = self.kernels[0]
|
k = self.kernels[0]
|
||||||
m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
|
m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
|
||||||
M=self.num_inducing, kernel=k)
|
num_inducing=self.num_inducing, kernel=k)
|
||||||
assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
|
assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
|
||||||
def testPsi2_lin_bia(self):
|
def testPsi2_lin_bia(self):
|
||||||
k = self.kernels[3]
|
k = self.kernels[3]
|
||||||
m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
|
m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
|
||||||
M=self.num_inducing, kernel=k)
|
num_inducing=self.num_inducing, kernel=k)
|
||||||
assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
|
assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
|
||||||
def testPsi2_rbf(self):
|
def testPsi2_rbf(self):
|
||||||
k = self.kernels[1]
|
k = self.kernels[1]
|
||||||
m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
|
m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
|
||||||
M=self.num_inducing, kernel=k)
|
num_inducing=self.num_inducing, kernel=k)
|
||||||
assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
|
assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
|
||||||
def testPsi2_rbf_bia(self):
|
def testPsi2_rbf_bia(self):
|
||||||
k = self.kernels[-1]
|
k = self.kernels[-1]
|
||||||
m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
|
m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
|
||||||
M=self.num_inducing, kernel=k)
|
num_inducing=self.num_inducing, kernel=k)
|
||||||
assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
|
assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
|
||||||
def testPsi2_bia(self):
|
def testPsi2_bia(self):
|
||||||
k = self.kernels[2]
|
k = self.kernels[2]
|
||||||
m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
|
m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
|
||||||
M=self.num_inducing, kernel=k)
|
num_inducing=self.num_inducing, kernel=k)
|
||||||
assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
|
assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -122,11 +122,11 @@ if __name__ == "__main__":
|
||||||
numpy.random.seed(0)
|
numpy.random.seed(0)
|
||||||
input_dim = 5
|
input_dim = 5
|
||||||
N = 50
|
N = 50
|
||||||
M = 10
|
num_inducing = 10
|
||||||
D = 15
|
D = 15
|
||||||
X = numpy.random.randn(N, input_dim)
|
X = numpy.random.randn(N, input_dim)
|
||||||
X_var = .5 * numpy.ones_like(X) + .1 * numpy.clip(numpy.random.randn(*X.shape), 0, 1)
|
X_var = .5 * numpy.ones_like(X) + .1 * numpy.clip(numpy.random.randn(*X.shape), 0, 1)
|
||||||
Z = numpy.random.permutation(X)[:M]
|
Z = numpy.random.permutation(X)[:num_inducing]
|
||||||
Y = X.dot(numpy.random.randn(input_dim, D))
|
Y = X.dot(numpy.random.randn(input_dim, D))
|
||||||
# kernel = GPy.kern.bias(input_dim)
|
# kernel = GPy.kern.bias(input_dim)
|
||||||
#
|
#
|
||||||
|
|
@ -148,7 +148,7 @@ if __name__ == "__main__":
|
||||||
# m2 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
|
# m2 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
|
||||||
# num_inducing=num_inducing, kernel=GPy.kern.rbf(input_dim))
|
# num_inducing=num_inducing, kernel=GPy.kern.rbf(input_dim))
|
||||||
m3 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
|
m3 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
|
||||||
M=M, kernel=GPy.kern.linear(input_dim, ARD=True, variances=numpy.random.rand(input_dim)))
|
num_inducing=num_inducing, kernel=GPy.kern.linear(input_dim, ARD=True, variances=numpy.random.rand(input_dim)))
|
||||||
m3.ensure_default_constraints()
|
m3.ensure_default_constraints()
|
||||||
# + GPy.kern.bias(input_dim))
|
# + GPy.kern.bias(input_dim))
|
||||||
# m4 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
|
# m4 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
|
||||||
|
|
|
||||||
|
|
@ -8,38 +8,38 @@ from GPy.models.sparse_gplvm import SparseGPLVM
|
||||||
|
|
||||||
class sparse_GPLVMTests(unittest.TestCase):
|
class sparse_GPLVMTests(unittest.TestCase):
|
||||||
def test_bias_kern(self):
|
def test_bias_kern(self):
|
||||||
N, M, input_dim, D = 10, 3, 2, 4
|
N, num_inducing, input_dim, D = 10, 3, 2, 4
|
||||||
X = np.random.rand(N, input_dim)
|
X = np.random.rand(N, input_dim)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
||||||
k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
m = SparseGPLVM(Y, input_dim, kernel=k, M=M)
|
m = SparseGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
m.randomize()
|
m.randomize()
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
||||||
@unittest.skip('linear kernels do not have dKdiag_dX')
|
@unittest.skip('linear kernels do not have dKdiag_dX')
|
||||||
def test_linear_kern(self):
|
def test_linear_kern(self):
|
||||||
N, M, input_dim, D = 10, 3, 2, 4
|
N, num_inducing, input_dim, D = 10, 3, 2, 4
|
||||||
X = np.random.rand(N, input_dim)
|
X = np.random.rand(N, input_dim)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
||||||
k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
m = SparseGPLVM(Y, input_dim, kernel=k, M=M)
|
m = SparseGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
m.randomize()
|
m.randomize()
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
||||||
def test_rbf_kern(self):
|
def test_rbf_kern(self):
|
||||||
N, M, input_dim, D = 10, 3, 2, 4
|
N, num_inducing, input_dim, D = 10, 3, 2, 4
|
||||||
X = np.random.rand(N, input_dim)
|
X = np.random.rand(N, input_dim)
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
K = k.K(X)
|
K = k.K(X)
|
||||||
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
|
||||||
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
|
||||||
m = SparseGPLVM(Y, input_dim, kernel=k, M=M)
|
m = SparseGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
|
||||||
m.ensure_default_constraints()
|
m.ensure_default_constraints()
|
||||||
m.randomize()
|
m.randomize()
|
||||||
self.assertTrue(m.checkgrad())
|
self.assertTrue(m.checkgrad())
|
||||||
|
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue