apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-02 08:12:39 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

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

Browse source
This commit is contained in:
Alan Saul 2014-11-05 14:51:23 +00:00
commit 9480f71fe4
16 changed files with 207 additions and 651 deletions
Download patch file Download diff file Expand all files Collapse all files

11
GPy/core/model.py
View file

@ -263,7 +263,7 @@ class Model(Parameterized):
sgd.run() sgd.run()
self.optimization_runs.append(sgd) self.optimization_runs.append(sgd)
def _checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3): def _checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3, df_tolerance=1e-12):
""" """
Check the gradient of the ,odel by comparing to a numerical Check the gradient of the ,odel by comparing to a numerical
estimate. If the verbose flag is passed, individual estimate. If the verbose flag is passed, individual
@ -322,7 +322,7 @@ class Model(Parameterized):
except NotImplementedError: except NotImplementedError:
names = ['Variable %i' % i for i in range(len(x))] names = ['Variable %i' % i for i in range(len(x))]
# Prepare for pretty-printing # Prepare for pretty-printing
header = ['Name', 'Ratio', 'Difference', 'Analytical', 'Numerical'] header = ['Name', 'Ratio', 'Difference', 'Analytical', 'Numerical', 'dF_ratio']
max_names = max([len(names[i]) for i in range(len(names))] + [len(header[0])]) max_names = max([len(names[i]) for i in range(len(names))] + [len(header[0])])
float_len = 10 float_len = 10
cols = [max_names] cols = [max_names]
@ -359,6 +359,8 @@ class Model(Parameterized):
f1 = self._objective(xx) f1 = self._objective(xx)
xx[xind] -= 2.*step xx[xind] -= 2.*step
f2 = self._objective(xx) f2 = self._objective(xx)
df_ratio = np.abs((f1-f2)/min(f1,f2))
df_unstable = df_ratio<df_tolerance
numerical_gradient = (f1 - f2) / (2 * step) numerical_gradient = (f1 - f2) / (2 * step)
if np.all(gradient[xind] == 0): ratio = (f1 - f2) == gradient[xind] if np.all(gradient[xind] == 0): ratio = (f1 - f2) == gradient[xind]
else: ratio = (f1 - f2) / (2 * step * gradient[xind]) else: ratio = (f1 - f2) / (2 * step * gradient[xind])
@ -370,12 +372,15 @@ class Model(Parameterized):
else: else:
formatted_name = "\033[91m {0} \033[0m".format(names[nind]) formatted_name = "\033[91m {0} \033[0m".format(names[nind])
ret &= False ret &= False
if df_unstable:
formatted_name = "\033[94m {0} \033[0m".format(names[nind])
r = '%.6f' % float(ratio) r = '%.6f' % float(ratio)
d = '%.6f' % float(difference) d = '%.6f' % float(difference)
g = '%.6f' % gradient[xind] g = '%.6f' % gradient[xind]
ng = '%.6f' % float(numerical_gradient) ng = '%.6f' % float(numerical_gradient)
grad_string = "{0:<{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name, r, d, g, ng, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4]) df = '%1.e' % float(df_ratio)
grad_string = "{0:<{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}|{5:^{c5}}".format(formatted_name, r, d, g, ng, df, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4], c5=cols[5])
print grad_string print grad_string
self.optimizer_array = x self.optimizer_array = x

7
GPy/core/sparse_gp_mpi.py
View file

@ -44,16 +44,15 @@ class SparseGP_MPI(SparseGP):
super(SparseGP_MPI, self).__init__(X, Y, Z, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata, normalizer=normalizer) super(SparseGP_MPI, self).__init__(X, Y, Z, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata, normalizer=normalizer)
self.update_model(False) self.update_model(False)
self.link_parameter(self.X, index=0)
if variational_prior is not None: if variational_prior is not None:
self.link_parameter(variational_prior) self.link_parameter(variational_prior)
# self.X.fix()
self.mpi_comm = mpi_comm self.mpi_comm = mpi_comm
# Manage the data (Y) division # Manage the data (Y) division
if mpi_comm != None: if mpi_comm != None:
from ..util.mpi import divide_data from ..util.parallel import divide_data
N_start, N_end, N_list = divide_data(Y.shape[0], mpi_comm) N_start, N_end, N_list = divide_data(Y.shape[0], mpi_comm.rank, mpi_comm.size)
self.N_range = (N_start, N_end) self.N_range = (N_start, N_end)
self.N_list = np.array(N_list) self.N_list = np.array(N_list)
self.Y_local = self.Y[N_start:N_end] self.Y_local = self.Y[N_start:N_end]

2
GPy/examples/dimensionality_reduction.py
View file

@ -348,7 +348,7 @@ def ssgplvm_simulation(optimize=True, verbose=1,
D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9 D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9
_, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim)
Y = Ylist[0] Y = Ylist[0]
k = kern.Linear(Q, ARD=True, useGPU=useGPU)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
#k = kern.RBF(Q, ARD=True, lengthscale=10.) #k = kern.RBF(Q, ARD=True, lengthscale=10.)
m = SSGPLVM(Y, Q, init="pca", num_inducing=num_inducing, kernel=k) m = SSGPLVM(Y, Q, init="pca", num_inducing=num_inducing, kernel=k)
m.X.variance[:] = _np.random.uniform(0,.01,m.X.shape) m.X.variance[:] = _np.random.uniform(0,.01,m.X.shape)

32
GPy/examples/regression.py
View file

@ -14,7 +14,9 @@ import GPy
def olympic_marathon_men(optimize=True, plot=True): def olympic_marathon_men(optimize=True, plot=True):
"""Run a standard Gaussian process regression on the Olympic marathon data.""" """Run a standard Gaussian process regression on the Olympic marathon data."""
try:import pods try:import pods
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets' except ImportError:
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
return
data = pods.datasets.olympic_marathon_men() data = pods.datasets.olympic_marathon_men()
# create simple GP Model # create simple GP Model
@ -85,7 +87,9 @@ def epomeo_gpx(max_iters=200, optimize=True, plot=True):
to load in the data. to load in the data.
""" """
try:import pods try:import pods
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets' except ImportError:
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
return
data = pods.datasets.epomeo_gpx() data = pods.datasets.epomeo_gpx()
num_data_list = [] num_data_list = []
for Xpart in data['X']: for Xpart in data['X']:
@ -130,7 +134,9 @@ def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=1000
log_SNRs = np.linspace(-3., 4., resolution) log_SNRs = np.linspace(-3., 4., resolution)
try:import pods try:import pods
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets' except ImportError:
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
return
data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number) data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number)
# data['Y'] = data['Y'][0::2, :] # data['Y'] = data['Y'][0::2, :]
# data['X'] = data['X'][0::2, :] # data['X'] = data['X'][0::2, :]
@ -212,7 +218,9 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.RBF):
def olympic_100m_men(optimize=True, plot=True): def olympic_100m_men(optimize=True, plot=True):
"""Run a standard Gaussian process regression on the Rogers and Girolami olympics data.""" """Run a standard Gaussian process regression on the Rogers and Girolami olympics data."""
try:import pods try:import pods
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets' except ImportError:
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
return
data = pods.datasets.olympic_100m_men() data = pods.datasets.olympic_100m_men()
# create simple GP Model # create simple GP Model
@ -231,7 +239,9 @@ def olympic_100m_men(optimize=True, plot=True):
def toy_rbf_1d(optimize=True, plot=True): def toy_rbf_1d(optimize=True, plot=True):
"""Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance.""" """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
try:import pods try:import pods
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets' except ImportError:
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
return
data = pods.datasets.toy_rbf_1d() data = pods.datasets.toy_rbf_1d()
# create simple GP Model # create simple GP Model
@ -247,7 +257,9 @@ def toy_rbf_1d(optimize=True, plot=True):
def toy_rbf_1d_50(optimize=True, plot=True): def toy_rbf_1d_50(optimize=True, plot=True):
"""Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance.""" """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
try:import pods try:import pods
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets' except ImportError:
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
return
data = pods.datasets.toy_rbf_1d_50() data = pods.datasets.toy_rbf_1d_50()
# create simple GP Model # create simple GP Model
@ -364,7 +376,9 @@ def toy_ARD_sparse(max_iters=1000, kernel_type='linear', num_samples=300, D=4, o
def robot_wireless(max_iters=100, kernel=None, optimize=True, plot=True): def robot_wireless(max_iters=100, kernel=None, optimize=True, plot=True):
"""Predict the location of a robot given wirelss signal strength readings.""" """Predict the location of a robot given wirelss signal strength readings."""
try:import pods try:import pods
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets' except ImportError:
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
return
data = pods.datasets.robot_wireless() data = pods.datasets.robot_wireless()
# create simple GP Model # create simple GP Model
@ -390,7 +404,9 @@ def robot_wireless(max_iters=100, kernel=None, optimize=True, plot=True):
def silhouette(max_iters=100, optimize=True, plot=True): def silhouette(max_iters=100, optimize=True, plot=True):
"""Predict the pose of a figure given a silhouette. This is a task from Agarwal and Triggs 2004 ICML paper.""" """Predict the pose of a figure given a silhouette. This is a task from Agarwal and Triggs 2004 ICML paper."""
try:import pods try:import pods
except ImportError:print 'pods unavailable, see https://github.com/sods/ods for example datasets' except ImportError:
print 'pods unavailable, see https://github.com/sods/ods for example datasets'
return
data = pods.datasets.silhouette() data = pods.datasets.silhouette()
# create simple GP Model # create simple GP Model

1
GPy/inference/latent_function_inference/__init__.py
View file

@ -66,7 +66,6 @@ from expectation_propagation_dtc import EPDTC
from dtc import DTC from dtc import DTC
from fitc import FITC from fitc import FITC
from var_dtc_parallel import VarDTC_minibatch from var_dtc_parallel import VarDTC_minibatch
from var_dtc_gpu import VarDTC_GPU
# class FullLatentFunctionData(object): # class FullLatentFunctionData(object):
# #

482
GPy/inference/latent_function_inference/var_dtc_gpu.py
View file

@ -1,482 +0,0 @@
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt)
from posterior import Posterior
from ...util.linalg import jitchol, backsub_both_sides, tdot, dtrtrs
from ...util import diag
from ...core.parameterization.variational import VariationalPosterior
import numpy as np
from . import LatentFunctionInference
log_2_pi = np.log(2*np.pi)
from ...util import gpu_init
try:
import scikits.cuda.linalg as culinalg
import pycuda.gpuarray as gpuarray
from scikits.cuda import cublas
from ...util.linalg_gpu import logDiagSum, strideSum, mul_bcast, sum_axis, outer_prod, mul_bcast_first, join_prod, traceDot
except:
pass
class VarDTC_GPU(LatentFunctionInference):
"""
An object for inference when the likelihood is Gaussian, but we want to do sparse inference.
The function self.inference returns a Posterior object, which summarizes
the posterior.
For efficiency, we sometimes work with the cholesky of Y*Y.T. To save repeatedly recomputing this, we cache it.
"""
const_jitter = np.float64(1e-6)
def __init__(self, batchsize=None, gpu_memory=4., limit=1):
self.batchsize = batchsize
self.gpu_memory = gpu_memory
self.midRes = {}
self.batch_pos = 0 # the starting position of the current mini-batch
self.cublas_handle = gpu_init.cublas_handle
# Initialize GPU caches
self.gpuCache = None
def _initGPUCache(self, kern, num_inducing, input_dim, output_dim, Y):
ndata = Y.shape[0]
if self.batchsize==None:
self.batchsize = self._estimateBatchSize(kern, ndata, num_inducing, input_dim, output_dim)
if self.gpuCache == None:
self.gpuCache = {# inference_likelihood
'Kmm_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
'Lm_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
'ones_gpu' :gpuarray.empty(num_inducing, np.float64,order='F'),
'LL_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
'b_gpu' :gpuarray.empty((num_inducing,output_dim),np.float64,order='F'),
'v_gpu' :gpuarray.empty((num_inducing,output_dim),np.float64,order='F'),
'vvt_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
'KmmInvPsi2LLInvT_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
'KmmInvPsi2P_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
'dL_dpsi2R_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
'dL_dKmm_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
'psi1Y_gpu' :gpuarray.empty((num_inducing,output_dim),np.float64,order='F'),
'psi2_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
'beta_gpu' :gpuarray.empty((ndata,),np.float64,order='F'),
'YT_gpu' :gpuarray.to_gpu(np.asfortranarray(Y.T)), # DxN
'betaYT_gpu' :gpuarray.empty(Y.T.shape,np.float64,order='F'), # DxN
# inference_minibatch
'dL_dpsi0_gpu' :gpuarray.empty((self.batchsize,),np.float64,order='F'),
'dL_dpsi1_gpu' :gpuarray.empty((self.batchsize,num_inducing),np.float64,order='F'),
'dL_dpsi2_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
'psi0p_gpu' :gpuarray.empty((self.batchsize,),np.float64,order='F'),
'psi1p_gpu' :gpuarray.empty((self.batchsize,num_inducing),np.float64,order='F'),
'psi2p_gpu' :gpuarray.empty((num_inducing,num_inducing),np.float64,order='F'),
}
self.gpuCache['ones_gpu'].fill(1.0)
YT_gpu = self.gpuCache['YT_gpu']
self._trYYT = cublas.cublasDdot(self.cublas_handle, YT_gpu.size, YT_gpu.gpudata, 1, YT_gpu.gpudata, 1)
def _estimateMemoryOccupation(self, N, M, D):
"""
Estimate the best batch size.
N - the number of total datapoints
M - the number of inducing points
D - the number of observed (output) dimensions
return: the constant memory size, the memory occupation of batchsize=1
unit: GB
"""
return (M+9.*M*M+3*M*D+N+2.*N*D)*8./1024./1024./1024., (4.+3.*M+D+3.*M*M)*8./1024./1024./1024.
def _estimateBatchSize(self, kern, N, M, Q, D):
"""
Estimate the best batch size.
N - the number of total datapoints
M - the number of inducing points
D - the number of observed (output) dimensions
return: the constant memory size, the memory occupation of batchsize=1
unit: GB
"""
if kern.useGPU:
x0,x1 = kern.psicomp.estimateMemoryOccupation(N,M,Q)
else:
x0, x1 = 0.,0.
y0, y1 = self._estimateMemoryOccupation(N, M, D)
opt_batchsize = min(int((self.gpu_memory-y0-x0)/(x1+y1)), N)
return opt_batchsize
def _get_YYTfactor(self, Y):
"""
find a matrix L which satisfies LLT = YYT.
Note that L may have fewer columns than Y.
"""
N, D = Y.shape
if (N>=D):
return Y.view(np.ndarray)
else:
return jitchol(tdot(Y))
def gatherPsiStat(self, kern, X, Z, Y, beta, uncertain_inputs, het_noise):
num_inducing, input_dim = Z.shape[0], Z.shape[1]
num_data, output_dim = Y.shape
trYYT = self._trYYT
psi1Y_gpu = self.gpuCache['psi1Y_gpu']
psi2_gpu = self.gpuCache['psi2_gpu']
beta_gpu = self.gpuCache['beta_gpu']
YT_gpu = self.gpuCache['YT_gpu']
betaYT_gpu = self.gpuCache['betaYT_gpu']
beta_gpu.fill(beta)
betaYT_gpu.fill(0.)
cublas.cublasDaxpy(self.cublas_handle, betaYT_gpu.size, beta, YT_gpu.gpudata, 1, betaYT_gpu.gpudata, 1)
YRY_full = trYYT*beta
if kern.useGPU:
psi1Y_gpu.fill(0.)
psi2_gpu.fill(0.)
psi0_full = 0
for n_start in xrange(0,num_data,self.batchsize):
n_end = min(self.batchsize+n_start, num_data)
ndata = n_end - n_start
X_slice = X[n_start:n_end]
betaYT_gpu_slice = betaYT_gpu[:,n_start:n_end]
if uncertain_inputs:
psi0 = kern.psi0(Z, X_slice)
psi1p_gpu = kern.psi1(Z, X_slice)
psi2p_gpu = kern.psi2(Z, X_slice)
else:
psi0 = kern.Kdiag(X_slice)
psi1p_gpu = kern.K(X_slice, Z)
cublas.cublasDgemm(self.cublas_handle, 'T', 'T', num_inducing, output_dim, ndata, 1.0, psi1p_gpu.gpudata, ndata, betaYT_gpu_slice.gpudata, output_dim, 1.0, psi1Y_gpu.gpudata, num_inducing)
psi0_full += psi0.sum()
if uncertain_inputs:
sum_axis(psi2_gpu,psi2p_gpu,1,1)
else:
cublas.cublasDgemm(self.cublas_handle, 'T', 'N', num_inducing, num_inducing, ndata, beta, psi1p_gpu.gpudata, ndata, psi1p_gpu.gpudata, ndata, 1.0, psi2_gpu.gpudata, num_inducing)
psi0_full *= beta
if uncertain_inputs:
cublas.cublasDscal(self.cublas_handle, psi2_gpu.size, beta, psi2_gpu.gpudata, 1)
else:
psi2_full = np.zeros((num_inducing,num_inducing))
psi1Y_full = np.zeros((output_dim,num_inducing)) # DxM
psi0_full = 0.
YRY_full = 0.
for n_start in xrange(0,num_data,self.batchsize):
n_end = min(self.batchsize+n_start, num_data)
Y_slice = Y[n_start:n_end]
X_slice = X[n_start:n_end]
if het_noise:
b = beta[n_start]
YRY_full += np.inner(Y_slice, Y_slice)*b
else:
b = beta
if uncertain_inputs:
psi0 = kern.psi0(Z, X_slice)
psi1 = kern.psi1(Z, X_slice)
psi2_full += kern.psi2(Z, X_slice)*b
else:
psi0 = kern.Kdiag(X_slice)
psi1 = kern.K(X_slice, Z)
psi2_full += np.dot(psi1.T,psi1)*b
psi0_full += psi0.sum()*b
psi1Y_full += np.dot(Y_slice.T,psi1)*b # DxM
if not het_noise:
YRY_full = trYYT*beta
psi1Y_gpu.set(psi1Y_full)
psi2_gpu.set(psi2_full)
return psi0_full, YRY_full
def inference_likelihood(self, kern, X, Z, likelihood, Y):
"""
The first phase of inference:
Compute: log-likelihood, dL_dKmm
Cached intermediate results: Kmm, KmmInv,
"""
num_inducing, input_dim = Z.shape[0], Z.shape[1]
num_data, output_dim = Y.shape
#see whether we've got a different noise variance for each datum
beta = 1./np.fmax(likelihood.variance, 1e-6)
het_noise = beta.size > 1
if het_noise:
self.batchsize=0
self._initGPUCache(kern, num_inducing, input_dim, output_dim, Y)
if isinstance(X, VariationalPosterior):
uncertain_inputs = True
else:
uncertain_inputs = False
psi1Y_gpu = self.gpuCache['psi1Y_gpu']
psi2_gpu = self.gpuCache['psi2_gpu']
psi0_full, YRY_full = self.gatherPsiStat(kern, X, Z, Y, beta, uncertain_inputs, het_noise)
#======================================================================
# Compute Common Components
#======================================================================
Kmm = kern.K(Z).copy()
Kmm_gpu = self.gpuCache['Kmm_gpu']
Kmm_gpu.set(np.asfortranarray(Kmm))
diag.add(Kmm, self.const_jitter)
ones_gpu = self.gpuCache['ones_gpu']
cublas.cublasDaxpy(self.cublas_handle, num_inducing, self.const_jitter, ones_gpu.gpudata, 1, Kmm_gpu.gpudata, num_inducing+1)
# assert np.allclose(Kmm, Kmm_gpu.get())
# Lm = jitchol(Kmm)
#
Lm_gpu = self.gpuCache['Lm_gpu']
cublas.cublasDcopy(self.cublas_handle, Kmm_gpu.size, Kmm_gpu.gpudata, 1, Lm_gpu.gpudata, 1)
culinalg.cho_factor(Lm_gpu,'L')
# print np.abs(np.tril(Lm)-np.tril(Lm_gpu.get())).max()
# Lambda = Kmm+psi2_full
# LL = jitchol(Lambda)
#
Lambda_gpu = self.gpuCache['LL_gpu']
cublas.cublasDcopy(self.cublas_handle, Kmm_gpu.size, Kmm_gpu.gpudata, 1, Lambda_gpu.gpudata, 1)
cublas.cublasDaxpy(self.cublas_handle, psi2_gpu.size, np.float64(1.0), psi2_gpu.gpudata, 1, Lambda_gpu.gpudata, 1)
LL_gpu = Lambda_gpu
culinalg.cho_factor(LL_gpu,'L')
# print np.abs(np.tril(LL)-np.tril(LL_gpu.get())).max()
# b,_ = dtrtrs(LL, psi1Y_full)
# bbt_cpu = np.square(b).sum()
#
b_gpu = self.gpuCache['b_gpu']
cublas.cublasDcopy(self.cublas_handle, b_gpu.size, psi1Y_gpu.gpudata, 1, b_gpu.gpudata, 1)
cublas.cublasDtrsm(self.cublas_handle , 'L', 'L', 'N', 'N', num_inducing, output_dim, np.float64(1.0), LL_gpu.gpudata, num_inducing, b_gpu.gpudata, num_inducing)
bbt = cublas.cublasDdot(self.cublas_handle, b_gpu.size, b_gpu.gpudata, 1, b_gpu.gpudata, 1)
# print np.abs(bbt-bbt_cpu)
# v,_ = dtrtrs(LL.T,b,lower=False)
# vvt = np.einsum('md,od->mo',v,v)
# LmInvPsi2LmInvT = backsub_both_sides(Lm,psi2_full,transpose='right')
#
v_gpu = self.gpuCache['v_gpu']
cublas.cublasDcopy(self.cublas_handle, v_gpu.size, b_gpu.gpudata, 1, v_gpu.gpudata, 1)
cublas.cublasDtrsm(self.cublas_handle , 'L', 'L', 'T', 'N', num_inducing, output_dim, np.float64(1.0), LL_gpu.gpudata, num_inducing, v_gpu.gpudata, num_inducing)
vvt_gpu = self.gpuCache['vvt_gpu']
cublas.cublasDgemm(self.cublas_handle, 'N', 'T', num_inducing, num_inducing, output_dim, np.float64(1.0), v_gpu.gpudata, num_inducing, v_gpu.gpudata, num_inducing, np.float64(0.), vvt_gpu.gpudata, num_inducing)
LmInvPsi2LmInvT_gpu = self.gpuCache['KmmInvPsi2LLInvT_gpu']
cublas.cublasDcopy(self.cublas_handle, psi2_gpu.size, psi2_gpu.gpudata, 1, LmInvPsi2LmInvT_gpu.gpudata, 1)
cublas.cublasDtrsm(self.cublas_handle , 'L', 'L', 'N', 'N', num_inducing, num_inducing, np.float64(1.0), Lm_gpu.gpudata, num_inducing, LmInvPsi2LmInvT_gpu.gpudata, num_inducing)
cublas.cublasDtrsm(self.cublas_handle , 'r', 'L', 'T', 'N', num_inducing, num_inducing, np.float64(1.0), Lm_gpu.gpudata, num_inducing, LmInvPsi2LmInvT_gpu.gpudata, num_inducing)
#tr_LmInvPsi2LmInvT = cublas.cublasDasum(self.cublas_handle, num_inducing, LmInvPsi2LmInvT_gpu.gpudata, num_inducing+1)
tr_LmInvPsi2LmInvT = float(strideSum(LmInvPsi2LmInvT_gpu, num_inducing+1).get())
# print np.abs(vvt-vvt_gpu.get()).max()
# print np.abs(np.trace(LmInvPsi2LmInvT)-tr_LmInvPsi2LmInvT)
# Psi2LLInvT = dtrtrs(LL,psi2_full)[0].T
# LmInvPsi2LLInvT= dtrtrs(Lm,Psi2LLInvT)[0]
# KmmInvPsi2LLInvT = dtrtrs(Lm,LmInvPsi2LLInvT,trans=True)[0]
# KmmInvPsi2P = dtrtrs(LL,KmmInvPsi2LLInvT.T, trans=True)[0].T
#
KmmInvPsi2LLInvT_gpu = LmInvPsi2LmInvT_gpu # Reuse GPU memory (size:MxM)
cublas.cublasDcopy(self.cublas_handle, psi2_gpu.size, psi2_gpu.gpudata, 1, KmmInvPsi2LLInvT_gpu.gpudata, 1)
cublas.cublasDtrsm(self.cublas_handle , 'L', 'L', 'N', 'N', num_inducing, num_inducing, np.float64(1.0), Lm_gpu.gpudata, num_inducing, KmmInvPsi2LLInvT_gpu.gpudata, num_inducing)
cublas.cublasDtrsm(self.cublas_handle , 'r', 'L', 'T', 'N', num_inducing, num_inducing, np.float64(1.0), LL_gpu.gpudata, num_inducing, KmmInvPsi2LLInvT_gpu.gpudata, num_inducing)
cublas.cublasDtrsm(self.cublas_handle , 'L', 'L', 'T', 'N', num_inducing, num_inducing, np.float64(1.0), Lm_gpu.gpudata, num_inducing, KmmInvPsi2LLInvT_gpu.gpudata, num_inducing)
KmmInvPsi2P_gpu = self.gpuCache['KmmInvPsi2P_gpu']
cublas.cublasDcopy(self.cublas_handle, KmmInvPsi2LLInvT_gpu.size, KmmInvPsi2LLInvT_gpu.gpudata, 1, KmmInvPsi2P_gpu.gpudata, 1)
cublas.cublasDtrsm(self.cublas_handle , 'r', 'L', 'N', 'N', num_inducing, num_inducing, np.float64(1.0), LL_gpu.gpudata, num_inducing, KmmInvPsi2P_gpu.gpudata, num_inducing)
# print np.abs(KmmInvPsi2P-KmmInvPsi2P_gpu.get()).max()
# dL_dpsi2R = (output_dim*KmmInvPsi2P - vvt)/2. # dL_dpsi2 with R inside psi2
#
dL_dpsi2R_gpu = self.gpuCache['dL_dpsi2R_gpu']
cublas.cublasDcopy(self.cublas_handle, vvt_gpu.size, vvt_gpu.gpudata, 1, dL_dpsi2R_gpu.gpudata, 1)
cublas.cublasDaxpy(self.cublas_handle, KmmInvPsi2P_gpu.size, np.float64(-output_dim), KmmInvPsi2P_gpu.gpudata, 1, dL_dpsi2R_gpu.gpudata, 1)
cublas.cublasDscal(self.cublas_handle, dL_dpsi2R_gpu.size, np.float64(-0.5), dL_dpsi2R_gpu.gpudata, 1)
# print np.abs(dL_dpsi2R_gpu.get()-dL_dpsi2R).max()
#======================================================================
# Compute log-likelihood
#======================================================================
if het_noise:
logL_R = -np.log(beta).sum()
else:
logL_R = -num_data*np.log(beta)
# logL_old = -(output_dim*(num_data*log_2_pi+logL_R+psi0_full-np.trace(LmInvPsi2LmInvT))+YRY_full-bbt)/2.-output_dim*(-np.log(np.diag(Lm)).sum()+np.log(np.diag(LL)).sum())
logdetKmm = float(logDiagSum(Lm_gpu,num_inducing+1).get())
logdetLambda = float(logDiagSum(LL_gpu,num_inducing+1).get())
logL = -(output_dim*(num_data*log_2_pi+logL_R+psi0_full-tr_LmInvPsi2LmInvT)+YRY_full-bbt)/2.+output_dim*(logdetKmm-logdetLambda)
# print np.abs(logL_old - logL)
#======================================================================
# Compute dL_dKmm
#======================================================================
# dL_dKmm = -(output_dim*np.einsum('md,od->mo',KmmInvPsi2LLInvT,KmmInvPsi2LLInvT) + vvt)/2.
#
dL_dKmm_gpu = self.gpuCache['dL_dKmm_gpu']
cublas.cublasDgemm(self.cublas_handle, 'N', 'T', num_inducing, num_inducing, num_inducing, np.float64(1.0), KmmInvPsi2LLInvT_gpu.gpudata, num_inducing, KmmInvPsi2LLInvT_gpu.gpudata, num_inducing, np.float64(0.), dL_dKmm_gpu.gpudata, num_inducing)
cublas.cublasDaxpy(self.cublas_handle, dL_dKmm_gpu.size, np.float64(1./output_dim), vvt_gpu.gpudata, 1, dL_dKmm_gpu.gpudata, 1)
cublas.cublasDscal(self.cublas_handle, dL_dKmm_gpu.size, np.float64(-output_dim/2.), dL_dKmm_gpu.gpudata, 1)
# print np.abs(dL_dKmm - dL_dKmm_gpu.get()).max()
#======================================================================
# Compute the Posterior distribution of inducing points p(u|Y)
#======================================================================
post = Posterior(woodbury_inv=KmmInvPsi2P_gpu.get(), woodbury_vector=v_gpu.get(), K=Kmm_gpu.get(), mean=None, cov=None, K_chol=Lm_gpu.get())
#======================================================================
# Compute dL_dthetaL for uncertian input and non-heter noise
#======================================================================
if not het_noise:
dL_dthetaL = (YRY_full + output_dim*psi0_full - num_data*output_dim)/-2.
dL_dthetaL += cublas.cublasDdot(self.cublas_handle,dL_dpsi2R_gpu.size, dL_dpsi2R_gpu.gpudata,1,psi2_gpu.gpudata,1)
dL_dthetaL += cublas.cublasDdot(self.cublas_handle,v_gpu.size, v_gpu.gpudata,1,psi1Y_gpu.gpudata,1)
self.midRes['dL_dthetaL'] = -beta*dL_dthetaL
return logL, dL_dKmm_gpu.get(), post
def inference_minibatch(self, kern, X, Z, likelihood, Y):
"""
The second phase of inference: Computing the derivatives over a minibatch of Y
Compute: dL_dpsi0, dL_dpsi1, dL_dpsi2, dL_dthetaL
return a flag showing whether it reached the end of Y (isEnd)
"""
num_data, output_dim = Y.shape
num_inducing = Z.shape[0]
if isinstance(X, VariationalPosterior):
uncertain_inputs = True
else:
uncertain_inputs = False
beta = 1./np.fmax(likelihood.variance, 1e-6)
het_noise = beta.size > 1
n_start = self.batch_pos
n_end = min(self.batchsize+n_start, num_data)
if n_end==num_data:
isEnd = True
self.batch_pos = 0
else:
isEnd = False
self.batch_pos = n_end
nSlice = n_end-n_start
X_slice = X[n_start:n_end]
if het_noise:
beta = beta[n_start] # nSlice==1
if kern.useGPU:
if not uncertain_inputs:
psi0p_gpu = kern.Kdiag(X_slice)
psi1p_gpu = kern.K(X_slice, Z)
psi2p_gpu = self.gpuCache['psi2p_gpu']
elif het_noise:
psi0p_gpu = kern.psi0(Z, X_slice)
psi1p_gpu = kern.psi1(Z, X_slice)
psi2p_gpu = kern.psi2(Z, X_slice)
elif not uncertain_inputs or het_noise:
if not uncertain_inputs:
psi0 = kern.Kdiag(X_slice)
psi1 = kern.K(X_slice, Z)
elif het_noise:
psi0 = kern.psi0(Z, X_slice)
psi1 = kern.psi1(Z, X_slice)
psi2 = kern.psi2(Z, X_slice)
psi0p_gpu = self.gpuCache['psi0p_gpu']
psi1p_gpu = self.gpuCache['psi1p_gpu']
psi2p_gpu = self.gpuCache['psi2p_gpu']
if psi0p_gpu.shape[0] > nSlice:
psi0p_gpu = psi0p_gpu[:nSlice]
psi1p_gpu = psi1p_gpu.ravel()[:nSlice*num_inducing].reshape(nSlice,num_inducing)
psi0p_gpu.set(np.asfortranarray(psi0))
psi1p_gpu.set(np.asfortranarray(psi1))
if uncertain_inputs:
psi2p_gpu.set(np.asfortranarray(psi2))
#======================================================================
# Compute dL_dpsi
#======================================================================
dL_dpsi2R_gpu = self.gpuCache['dL_dpsi2R_gpu']
v_gpu = self.gpuCache['v_gpu']
dL_dpsi0_gpu = self.gpuCache['dL_dpsi0_gpu']
dL_dpsi1_gpu = self.gpuCache['dL_dpsi1_gpu']
dL_dpsi2_gpu = self.gpuCache['dL_dpsi2_gpu']
betaYT_gpu = self.gpuCache['betaYT_gpu']
betaYT_gpu_slice = betaYT_gpu[:,n_start:n_end]
# Adjust to the batch size
if dL_dpsi0_gpu.shape[0] > nSlice:
dL_dpsi0_gpu = dL_dpsi0_gpu.ravel()[:nSlice]
dL_dpsi1_gpu = dL_dpsi1_gpu.ravel()[:nSlice*num_inducing].reshape(nSlice,num_inducing)
dL_dpsi0_gpu.fill(-output_dim *beta/2.)
cublas.cublasDgemm(self.cublas_handle, 'T', 'T', nSlice, num_inducing, output_dim, 1.0, betaYT_gpu_slice.gpudata, output_dim, v_gpu.gpudata, num_inducing, 0., dL_dpsi1_gpu.gpudata, nSlice)
if uncertain_inputs:
cublas.cublasDcopy(self.cublas_handle, dL_dpsi2R_gpu.size, dL_dpsi2R_gpu.gpudata, 1, dL_dpsi2_gpu.gpudata, 1)
cublas.cublasDscal(self.cublas_handle, dL_dpsi2_gpu.size, beta, dL_dpsi2_gpu.gpudata, 1)
else:
cublas.cublasDgemm(self.cublas_handle, 'N', 'N', nSlice, num_inducing, output_dim, beta, psi1p_gpu.gpudata, nSlice, dL_dpsi2R_gpu.gpudata, num_inducing, 1.0, dL_dpsi1_gpu.gpudata, nSlice)
#======================================================================
# Compute dL_dthetaL
#======================================================================
if het_noise:
betaY = betaYT_gpu_slice.get()
dL_dthetaL = ((np.square(betaY)).sum(axis=-1) + np.square(beta)*(output_dim*psi0p_gpu.get())-output_dim*beta)/2.
dL_dthetaL += -beta*beta*cublas.cublasDdot(self.cublas_handle,dL_dpsi2R_gpu.size, dL_dpsi2R_gpu.gpudata,1,psi2p_gpu.gpudata,1)
dL_dthetaL += -beta*(betaY*np.dot(psi1p_gpu.get(),v_gpu.get())).sum(axis=-1)
if kern.useGPU:
dL_dpsi0 = dL_dpsi0_gpu
dL_dpsi1 = dL_dpsi1_gpu
else:
dL_dpsi0 = dL_dpsi0_gpu.get()
dL_dpsi1 = dL_dpsi1_gpu.get()
if uncertain_inputs:
if kern.useGPU:
dL_dpsi2 = dL_dpsi2_gpu
else:
dL_dpsi2 = dL_dpsi2_gpu.get()
if not het_noise:
if isEnd:
dL_dthetaL = self.midRes['dL_dthetaL']
else:
dL_dthetaL = 0.
if uncertain_inputs:
grad_dict = {'dL_dpsi0':dL_dpsi0,
'dL_dpsi1':dL_dpsi1,
'dL_dpsi2':dL_dpsi2,
'dL_dthetaL':dL_dthetaL}
else:
grad_dict = {'dL_dKdiag':dL_dpsi0,
'dL_dKnm':dL_dpsi1,
'dL_dthetaL':dL_dthetaL}
return isEnd, (n_start,n_end), grad_dict

108
GPy/inference/latent_function_inference/var_dtc_parallel.py
View file

@ -38,7 +38,7 @@ class VarDTC_minibatch(LatentFunctionInference):
self.midRes = {} self.midRes = {}
self.batch_pos = 0 # the starting position of the current mini-batch self.batch_pos = 0 # the starting position of the current mini-batch
self.Y_speedup = False # Replace Y with the cholesky factor of YY.T, but the posterior inference will be wrong self.Y_speedup = False # Replace Y with the cholesky factor of YY.T, but the computation of posterior object will be skipped.
def __getstate__(self): def __getstate__(self):
# has to be overridden, as Cacher objects cannot be pickled. # has to be overridden, as Cacher objects cannot be pickled.
@ -76,6 +76,8 @@ class VarDTC_minibatch(LatentFunctionInference):
def gatherPsiStat(self, kern, X, Z, Y, beta, uncertain_inputs): def gatherPsiStat(self, kern, X, Z, Y, beta, uncertain_inputs):
het_noise = beta.size > 1 het_noise = beta.size > 1
assert beta.size == 1
trYYT = self.get_trYYT(Y) trYYT = self.get_trYYT(Y)
if self.Y_speedup and not het_noise: if self.Y_speedup and not het_noise:
@ -83,17 +85,16 @@ class VarDTC_minibatch(LatentFunctionInference):
num_inducing = Z.shape[0] num_inducing = Z.shape[0]
num_data, output_dim = Y.shape num_data, output_dim = Y.shape
if self.batchsize == None: batchsize = num_data if self.batchsize is None else self.batchsize
self.batchsize = num_data
psi2_full = np.zeros((num_inducing,num_inducing)) psi2_full = np.zeros((num_inducing,num_inducing)) # MxM
psi1Y_full = np.zeros((output_dim,num_inducing)) # DxM psi1Y_full = np.zeros((output_dim,num_inducing)) # DxM
psi0_full = 0. psi0_full = 0.
YRY_full = 0. YRY_full = 0.
for n_start in xrange(0,num_data,self.batchsize): for n_start in xrange(0,num_data,batchsize):
n_end = min(self.batchsize+n_start, num_data) n_end = min(batchsize+n_start, num_data)
if (n_end-n_start)==num_data: if batchsize==num_data:
Y_slice = Y Y_slice = Y
X_slice = X X_slice = X
else: else:
@ -168,16 +169,18 @@ class VarDTC_minibatch(LatentFunctionInference):
Kmm = kern.K(Z).copy() Kmm = kern.K(Z).copy()
diag.add(Kmm, self.const_jitter) diag.add(Kmm, self.const_jitter)
KmmInv,Lm,LmInv,_ = pdinv(Kmm) Lm = jitchol(Kmm, maxtries=100)
LmInvPsi2LmInvT = LmInv.dot(psi2_full).dot(LmInv.T) LmInvPsi2LmInvT = backsub_both_sides(Lm,psi2_full,transpose='right')
Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
LInv,LL,LLInv,logdet_L = pdinv(Lambda) LL = jitchol(Lambda, maxtries=100)
b = LLInv.dot(LmInv.dot(psi1Y_full.T)) logdet_L = 2.*np.sum(np.log(np.diag(LL)))
b = dtrtrs(LL,dtrtrs(Lm,psi1Y_full.T)[0])[0]
bbt = np.square(b).sum() bbt = np.square(b).sum()
v = LmInv.T.dot(LLInv.T.dot(b)) v = dtrtrs(Lm,dtrtrs(LL,b,trans=1)[0],trans=1)[0]
dL_dpsi2R = LmInv.T.dot(-LLInv.T.dot(tdot(b)+output_dim*np.eye(input_dim)).dot(LLInv)+output_dim*np.eye(input_dim)).dot(LmInv)/2. tmp = -backsub_both_sides(LL, tdot(b)+output_dim*np.eye(input_dim), transpose='left')
dL_dpsi2R = backsub_both_sides(Lm, tmp+output_dim*np.eye(input_dim), transpose='left')/2.
# Cache intermediate results # Cache intermediate results
self.midRes['dL_dpsi2R'] = dL_dpsi2R self.midRes['dL_dpsi2R'] = dL_dpsi2R
@ -196,14 +199,15 @@ class VarDTC_minibatch(LatentFunctionInference):
# Compute dL_dKmm # Compute dL_dKmm
#====================================================================== #======================================================================
dL_dKmm = dL_dpsi2R - output_dim*KmmInv.dot(psi2_full).dot(KmmInv)/2. dL_dKmm = dL_dpsi2R - output_dim*backsub_both_sides(Lm, LmInvPsi2LmInvT, transpose='left')/2.
#====================================================================== #======================================================================
# Compute the Posterior distribution of inducing points p(u|Y) # Compute the Posterior distribution of inducing points p(u|Y)
#====================================================================== #======================================================================
if not self.Y_speedup or het_noise: if not self.Y_speedup or het_noise:
post = Posterior(woodbury_inv=LmInv.T.dot(np.eye(input_dim)-LInv).dot(LmInv), woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=Lm) wd_inv = backsub_both_sides(Lm, np.eye(input_dim)- backsub_both_sides(LL, np.identity(input_dim), transpose='left'), transpose='left')
post = Posterior(woodbury_inv=wd_inv, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=Lm)
else: else:
post = None post = None
@ -242,7 +246,8 @@ class VarDTC_minibatch(LatentFunctionInference):
YYT_factor = Y YYT_factor = Y
n_start = self.batch_pos n_start = self.batch_pos
n_end = min(self.batchsize+n_start, num_data) batchsize = num_data if self.batchsize is None else self.batchsize
n_end = min(batchsize+n_start, num_data)
if n_end==num_data: if n_end==num_data:
isEnd = True isEnd = True
self.batch_pos = 0 self.batch_pos = 0
@ -250,8 +255,12 @@ class VarDTC_minibatch(LatentFunctionInference):
isEnd = False isEnd = False
self.batch_pos = n_end self.batch_pos = n_end
Y_slice = YYT_factor[n_start:n_end] if batchsize==num_data:
X_slice = X[n_start:n_end] Y_slice = YYT_factor
X_slice =X
else:
Y_slice = YYT_factor[n_start:n_end]
X_slice = X[n_start:n_end]
if not uncertain_inputs: if not uncertain_inputs:
psi0 = kern.Kdiag(X_slice) psi0 = kern.Kdiag(X_slice)
@ -405,3 +414,66 @@ def update_gradients(model, mpi_comm=None):
# dL_dthetaL # dL_dthetaL
model.likelihood.update_gradients(dL_dthetaL) model.likelihood.update_gradients(dL_dthetaL)
def update_gradients_sparsegp(model, mpi_comm=None):
if mpi_comm == None:
Y = model.Y
X = model.X
else:
Y = model.Y_local
X = model.X[model.N_range[0]:model.N_range[1]]
model._log_marginal_likelihood, dL_dKmm, model.posterior = model.inference_method.inference_likelihood(model.kern, X, model.Z, model.likelihood, Y)
het_noise = model.likelihood.variance.size > 1
if het_noise:
dL_dthetaL = np.empty((model.Y.shape[0],))
else:
dL_dthetaL = np.float64(0.)
kern_grad = model.kern.gradient.copy()
kern_grad[:] = 0.
model.Z.gradient = 0.
isEnd = False
while not isEnd:
isEnd, n_range, grad_dict = model.inference_method.inference_minibatch(model.kern, X, model.Z, model.likelihood, Y)
if (n_range[1]-n_range[0])==X.shape[0]:
X_slice = X
elif mpi_comm ==None:
X_slice = model.X[n_range[0]:n_range[1]]
else:
X_slice = model.X[model.N_range[0]+n_range[0]:model.N_range[0]+n_range[1]]
model.kern.update_gradients_diag(grad_dict['dL_dKdiag'], X_slice)
kern_grad += model.kern.gradient
model.kern.update_gradients_full(grad_dict['dL_dKnm'], X_slice, model.Z)
kern_grad += model.kern.gradient
model.Z.gradient += model.kern.gradients_X(grad_dict['dL_dKnm'].T, model.Z, X_slice)
if het_noise:
dL_dthetaL[n_range[0]:n_range[1]] = grad_dict['dL_dthetaL']
else:
dL_dthetaL += grad_dict['dL_dthetaL']
# Gather the gradients from multiple MPI nodes
if mpi_comm != None:
if het_noise:
raise "het_noise not implemented!"
kern_grad_all = kern_grad.copy()
Z_grad_all = model.Z.gradient.copy()
mpi_comm.Allreduce([kern_grad, MPI.DOUBLE], [kern_grad_all, MPI.DOUBLE])
mpi_comm.Allreduce([model.Z.gradient, MPI.DOUBLE], [Z_grad_all, MPI.DOUBLE])
kern_grad = kern_grad_all
model.Z.gradient = Z_grad_all
model.kern.update_gradients_full(dL_dKmm, model.Z, None)
model.kern.gradient += kern_grad
model.Z.gradient += model.kern.gradients_X(dL_dKmm, model.Z)
# dL_dthetaL
model.likelihood.update_gradients(dL_dthetaL)

76
GPy/kern/_src/hierarchical.py
View file

@ -1,76 +0,0 @@
# Copyright (c) 2012, James Hesnsman
# Licensed under the BSD 3-clause license (see LICENSE.txt)
from kernpart import Kernpart
import numpy as np
from independent_outputs import index_to_slices
class Hierarchical(Kernpart):
"""
A kernel part which can reopresent a hierarchy of indepencnce: a generalisation of independent_outputs
"""
def __init__(self,parts,name='hierarchy'):
self.levels = len(parts)
self.input_dim = parts[0].input_dim + 1
self.num_params = np.sum([k.num_params for k in parts])
self.name = name
self.parts = parts
self.param_starts = np.hstack((0,np.cumsum([k.num_params for k in self.parts[:-1]])))
self.param_stops = np.cumsum([k.num_params for k in self.parts])
def _get_params(self):
return np.hstack([k._get_params() for k in self.parts])
def _set_params(self,x):
[k._set_params(x[start:stop]) for k, start, stop in zip(self.parts, self.param_starts, self.param_stops)]
def _get_param_names(self):
return sum([[str(i)+'_'+k.name+'_'+n for n in k._get_param_names()] for i,k in enumerate(self.parts)],[])
def _sort_slices(self,X,X2):
slices = [index_to_slices(x) for x in X[:,-self.levels:].T]
X = X[:,:-self.levels]
if X2 is None:
slices2 = slices
X2 = X
else:
slices2 = [index_to_slices(x) for x in X2[:,-self.levels:].T]
X2 = X2[:,:-self.levels]
return X, X2, slices, slices2
def K(self,X,X2,target):
X, X2, slices, slices2 = self._sort_slices(X,X2)
[[[[k.K(X[s],X2[s2],target[s,s2]) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices_,slices2_)] for k, slices_, slices2_ in zip(self.parts,slices,slices2)]
def Kdiag(self,X,target):
raise NotImplementedError
#X,slices = X[:,:-1],index_to_slices(X[:,-1])
#[[self.k.Kdiag(X[s],target[s]) for s in slices_i] for slices_i in slices]
def _param_grad_helper(self,dL_dK,X,X2,target):
X, X2, slices, slices2 = self._sort_slices(X,X2)
[[[[k._param_grad_helper(dL_dK[s,s2],X[s],X2[s2],target[p_start:p_stop]) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices_, slices2_)] for k, p_start, p_stop, slices_, slices2_ in zip(self.parts, self.param_starts, self.param_stops, slices, slices2)]
def gradients_X(self,dL_dK,X,X2,target):
raise NotImplementedError
#X,slices = X[:,:-1],index_to_slices(X[:,-1])
#if X2 is None:
#X2,slices2 = X,slices
#else:
#X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
#[[[self.k.gradients_X(dL_dK[s,s2],X[s],X2[s2],target[s,:-1]) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)]
#
def dKdiag_dX(self,dL_dKdiag,X,target):
raise NotImplementedError
#X,slices = X[:,:-1],index_to_slices(X[:,-1])
#[[self.k.dKdiag_dX(dL_dKdiag[s],X[s],target[s,:-1]) for s in slices_i] for slices_i in slices]
def dKdiag_dtheta(self,dL_dKdiag,X,target):
raise NotImplementedError
#X,slices = X[:,:-1],index_to_slices(X[:,-1])
#[[self.k.dKdiag_dX(dL_dKdiag[s],X[s],target) for s in slices_i] for slices_i in slices]

13
GPy/kern/_src/kern.py
View file

@ -132,13 +132,20 @@ class Kern(Parameterized):
""" """
raise NotImplementedError raise NotImplementedError
def plot(self, *args, **kwargs): def plot(self, x=None, fignum=None, ax=None, title=None, plot_limits=None, resolution=None, **mpl_kwargs):
""" """
See GPy.plotting.matplot_dep.plot plot this kernel.
:param x: the value to use for the other kernel argument (kernels are a function of two variables!)
:param fignum: figure number of the plot
:param ax: matplotlib axis to plot on
:param title: the matplotlib title
:param plot_limits: the range over which to plot the kernel
:resolution: the resolution of the lines used in plotting
:mpl_kwargs avalid keyword arguments to pass through to matplotlib (e.g. lw=7)
""" """
assert "matplotlib" in sys.modules, "matplotlib package has not been imported." assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
from ...plotting.matplot_dep import kernel_plots from ...plotting.matplot_dep import kernel_plots
kernel_plots.plot(self,*args) kernel_plots.plot(self, x, fignum, ax, title, plot_limits, resolution, **mpl_kwargs)
def plot_ARD(self, *args, **kw): def plot_ARD(self, *args, **kw):
""" """

5
GPy/models/bayesian_gplvm.py
View file

@ -7,7 +7,6 @@ from ..core.sparse_gp_mpi import SparseGP_MPI
from ..likelihoods import Gaussian from ..likelihoods import Gaussian
from ..core.parameterization.variational import NormalPosterior, NormalPrior from ..core.parameterization.variational import NormalPosterior, NormalPrior
from ..inference.latent_function_inference.var_dtc_parallel import VarDTC_minibatch from ..inference.latent_function_inference.var_dtc_parallel import VarDTC_minibatch
from ..inference.latent_function_inference.var_dtc_gpu import VarDTC_GPU
import logging import logging
class BayesianGPLVM(SparseGP_MPI): class BayesianGPLVM(SparseGP_MPI):
@ -68,14 +67,12 @@ class BayesianGPLVM(SparseGP_MPI):
if isinstance(inference_method,VarDTC_minibatch): if isinstance(inference_method,VarDTC_minibatch):
inference_method.mpi_comm = mpi_comm inference_method.mpi_comm = mpi_comm
if kernel.useGPU and isinstance(inference_method, VarDTC_GPU):
kernel.psicomp.GPU_direct = True
super(BayesianGPLVM,self).__init__(X, Y, Z, kernel, likelihood=likelihood, super(BayesianGPLVM,self).__init__(X, Y, Z, kernel, likelihood=likelihood,
name=name, inference_method=inference_method, name=name, inference_method=inference_method,
normalizer=normalizer, mpi_comm=mpi_comm, normalizer=normalizer, mpi_comm=mpi_comm,
variational_prior=self.variational_prior, variational_prior=self.variational_prior,
) )
self.link_parameter(self.X, index=0)
def set_X_gradients(self, X, X_grad): def set_X_gradients(self, X, X_grad):
"""Set the gradients of the posterior distribution of X in its specific form.""" """Set the gradients of the posterior distribution of X in its specific form."""

1
GPy/models/bayesian_gplvm_minibatch.py
View file

@ -6,7 +6,6 @@ from .. import kern
from ..likelihoods import Gaussian from ..likelihoods import Gaussian
from ..core.parameterization.variational import NormalPosterior, NormalPrior from ..core.parameterization.variational import NormalPosterior, NormalPrior
from ..inference.latent_function_inference.var_dtc_parallel import VarDTC_minibatch from ..inference.latent_function_inference.var_dtc_parallel import VarDTC_minibatch
from ..inference.latent_function_inference.var_dtc_gpu import VarDTC_GPU
import logging import logging
from GPy.models.sparse_gp_minibatch import SparseGPMiniBatch from GPy.models.sparse_gp_minibatch import SparseGPMiniBatch

12
GPy/models/sparse_gp_regression.py
View file

@ -4,12 +4,13 @@
import numpy as np import numpy as np
from ..core import SparseGP from ..core import SparseGP
from ..core.sparse_gp_mpi import SparseGP_MPI
from .. import likelihoods from .. import likelihoods
from .. import kern from .. import kern
from ..inference.latent_function_inference import VarDTC from ..inference.latent_function_inference import VarDTC
from ..core.parameterization.variational import NormalPosterior from ..core.parameterization.variational import NormalPosterior
class SparseGPRegression(SparseGP): class SparseGPRegression(SparseGP_MPI):
""" """
Gaussian Process model for regression Gaussian Process model for regression
@ -48,7 +49,14 @@ class SparseGPRegression(SparseGP):
if not (X_variance is None): if not (X_variance is None):
X = NormalPosterior(X,X_variance) X = NormalPosterior(X,X_variance)
SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method=VarDTC(), normalizer=normalizer) SparseGP_MPI.__init__(self, X, Y, Z, kernel, likelihood, inference_method=VarDTC(), normalizer=normalizer)
def parameters_changed(self):
from ..inference.latent_function_inference.var_dtc_parallel import update_gradients_sparsegp,VarDTC_minibatch
if isinstance(self.inference_method,VarDTC_minibatch):
update_gradients_sparsegp(self, mpi_comm=self.mpi_comm)
else:
super(SparseGPRegression, self).parameters_changed()
class SparseGPRegressionUncertainInput(SparseGP): class SparseGPRegressionUncertainInput(SparseGP):
""" """

2
GPy/models/ss_gplvm.py
View file

@ -8,7 +8,6 @@ from .. import kern
from ..likelihoods import Gaussian from ..likelihoods import Gaussian
from ..core.parameterization.variational import SpikeAndSlabPrior, SpikeAndSlabPosterior from ..core.parameterization.variational import SpikeAndSlabPrior, SpikeAndSlabPosterior
from ..inference.latent_function_inference.var_dtc_parallel import update_gradients, VarDTC_minibatch from ..inference.latent_function_inference.var_dtc_parallel import update_gradients, VarDTC_minibatch
from ..inference.latent_function_inference.var_dtc_gpu import VarDTC_GPU
from ..kern._src.psi_comp.ssrbf_psi_gpucomp import PSICOMP_SSRBF_GPU from ..kern._src.psi_comp.ssrbf_psi_gpucomp import PSICOMP_SSRBF_GPU
class SSGPLVM(SparseGP_MPI): class SSGPLVM(SparseGP_MPI):
@ -72,6 +71,7 @@ class SSGPLVM(SparseGP_MPI):
super(SSGPLVM,self).__init__(X, Y, Z, kernel, likelihood, variational_prior=self.variational_prior, inference_method=inference_method, name=name, mpi_comm=mpi_comm, normalizer=normalizer, **kwargs) super(SSGPLVM,self).__init__(X, Y, Z, kernel, likelihood, variational_prior=self.variational_prior, inference_method=inference_method, name=name, mpi_comm=mpi_comm, normalizer=normalizer, **kwargs)
# self.X.unfix() # self.X.unfix()
# self.X.variance.constrain_positive() # self.X.variance.constrain_positive()
self.link_parameter(self.X, index=0)
if self.group_spike: if self.group_spike:
[self.X.gamma[:,i].tie('tieGamma'+str(i)) for i in xrange(self.X.gamma.shape[1])] # Tie columns together [self.X.gamma[:,i].tie('tieGamma'+str(i)) for i in xrange(self.X.gamma.shape[1])] # Tie columns together

41
GPy/plotting/matplot_dep/kernel_plots.py
View file

@ -99,7 +99,26 @@ def plot_ARD(kernel, fignum=None, ax=None, title='', legend=False, filtering=Non
return ax return ax
def plot(kernel, x=None, plot_limits=None, resolution=None, *args, **kwargs):
def plot(kernel,x=None, fignum=None, ax=None, title=None, plot_limits=None, resolution=None, **mpl_kwargs):
"""
plot a kernel.
:param x: the value to use for the other kernel argument (kernels are a function of two variables!)
:param fignum: figure number of the plot
:param ax: matplotlib axis to plot on
:param title: the matplotlib title
:param plot_limits: the range over which to plot the kernel
:resolution: the resolution of the lines used in plotting
:mpl_kwargs avalid keyword arguments to pass through to matplotlib (e.g. lw=7)
"""
fig, ax = ax_default(fignum,ax)
if title is None:
ax.set_title('%s kernel' % kernel.name)
else:
ax.set_title(title)
if kernel.input_dim == 1: if kernel.input_dim == 1:
if x is None: if x is None:
x = np.zeros((1, 1)) x = np.zeros((1, 1))
@ -117,10 +136,10 @@ def plot(kernel, x=None, plot_limits=None, resolution=None, *args, **kwargs):
Xnew = np.linspace(xmin, xmax, resolution or 201)[:, None] Xnew = np.linspace(xmin, xmax, resolution or 201)[:, None]
Kx = kernel.K(Xnew, x) Kx = kernel.K(Xnew, x)
pb.plot(Xnew, Kx, *args, **kwargs) ax.plot(Xnew, Kx, **mpl_kwargs)
pb.xlim(xmin, xmax) ax.set_xlim(xmin, xmax)
pb.xlabel("x") ax.set_xlabel("x")
pb.ylabel("k(x,%0.1f)" % x) ax.set_ylabel("k(x,%0.1f)" % x)
elif kernel.input_dim == 2: elif kernel.input_dim == 2:
if x is None: if x is None:
@ -142,11 +161,11 @@ def plot(kernel, x=None, plot_limits=None, resolution=None, *args, **kwargs):
Xnew = np.vstack((xx.flatten(), yy.flatten())).T Xnew = np.vstack((xx.flatten(), yy.flatten())).T
Kx = kernel.K(Xnew, x) Kx = kernel.K(Xnew, x)
Kx = Kx.reshape(resolution, resolution).T Kx = Kx.reshape(resolution, resolution).T
pb.contour(xx, xx, Kx, vmin=Kx.min(), vmax=Kx.max(), cmap=pb.cm.jet, *args, **kwargs) # @UndefinedVariable ax.contour(xx, yy, Kx, vmin=Kx.min(), vmax=Kx.max(), cmap=pb.cm.jet, **mpl_kwargs) # @UndefinedVariable
pb.xlim(xmin[0], xmax[0]) ax.set_xlim(xmin[0], xmax[0])
pb.ylim(xmin[1], xmax[1]) ax.set_ylim(xmin[1], xmax[1])
pb.xlabel("x1") ax.set_xlabel("x1")
pb.ylabel("x2") ax.set_ylabel("x2")
pb.title("k(x1,x2 ; %0.1f,%0.1f)" % (x[0, 0], x[0, 1])) ax.set_title("k(x1,x2 ; %0.1f,%0.1f)" % (x[0, 0], x[0, 1]))
else: else:
raise NotImplementedError, "Cannot plot a kernel with more than two input dimensions" raise NotImplementedError, "Cannot plot a kernel with more than two input dimensions"

35
GPy/util/mpi.py
View file

@ -1,35 +0,0 @@
"""
The tools for mpi
"""
try:
import numpy as np
from mpi4py import MPI
numpy_to_MPI_typemap = {
np.dtype(np.float64) : MPI.DOUBLE,
np.dtype(np.float32) : MPI.FLOAT,
np.dtype(np.int) : MPI.INT,
np.dtype(np.int8) : MPI.CHAR,
np.dtype(np.uint8) : MPI.UNSIGNED_CHAR,
np.dtype(np.int32) : MPI.INT,
np.dtype(np.uint32) : MPI.UNSIGNED_INT,
}
except:
pass
def divide_data(datanum, comm):
residue = (datanum)%comm.size
datanum_list = np.empty((comm.size),dtype=np.int32)
for i in xrange(comm.size):
if i<residue:
datanum_list[i] = int(datanum/comm.size)+1
else:
datanum_list[i] = int(datanum/comm.size)
if comm.rank<residue:
size = datanum/comm.size+1
offset = size*comm.rank
else:
size = datanum/comm.size
offset = size*comm.rank+residue
return offset, offset+size, datanum_list

30
GPy/util/parallel.py
View file

@ -1,7 +1,7 @@
""" """
The module of tools for parallelization (MPI) The module of tools for parallelization (MPI)
""" """
import numpy as np
try: try:
from mpi4py import MPI from mpi4py import MPI
def get_id_within_node(comm=MPI.COMM_WORLD): def get_id_within_node(comm=MPI.COMM_WORLD):
@ -9,5 +9,33 @@ try:
nodename = MPI.Get_processor_name() nodename = MPI.Get_processor_name()
nodelist = comm.allgather(nodename) nodelist = comm.allgather(nodename)
return len([i for i in nodelist[:rank] if i==nodename]) return len([i for i in nodelist[:rank] if i==nodename])
numpy_to_MPI_typemap = {
np.dtype(np.float64) : MPI.DOUBLE,
np.dtype(np.float32) : MPI.FLOAT,
np.dtype(np.int) : MPI.INT,
np.dtype(np.int8) : MPI.CHAR,
np.dtype(np.uint8) : MPI.UNSIGNED_CHAR,
np.dtype(np.int32) : MPI.INT,
np.dtype(np.uint32) : MPI.UNSIGNED_INT,
}
except: except:
pass pass
def divide_data(datanum, rank, size):
assert rank<size and datanum>0
residue = (datanum)%size
datanum_list = np.empty((size),dtype=np.int32)
for i in xrange(size):
if i<residue:
datanum_list[i] = int(datanum/size)+1
else:
datanum_list[i] = int(datanum/size)
if rank<residue:
size = datanum/size+1
offset = size*rank
else:
size = datanum/size
offset = size*rank+residue
return offset, offset+size, datanum_list