apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-05 01:32:40 +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 2013-11-27 13:21:14 +00:00
commit 811df6966f
6 changed files with 111 additions and 52 deletions
Download patch file Download diff file Expand all files Collapse all files

41
GPy/_models/gradient_checker.py
View file

@ -28,38 +28,37 @@ class GradientChecker(Model):
:param df: Gradient of function to check
:param x0:
Initial guess for inputs x (if it has a shape (a,b) this will be reflected in the parameter names).
Can be a list of arrays, if takes a list of arrays. This list will be passed
Can be a list of arrays, if f takes a list of arrays. This list will be passed
to f and df in the same order as given here.
If only one argument, make sure not to pass a list!!!
If f takes only one argument, make sure not to pass a list for x0!!!
:type x0: [array-like] | array-like | float | int
:param names:
:param list names:
Names to print, when performing gradcheck. If a list was passed to x0
a list of names with the same length is expected.
:param args: Arguments passed as f(x, *args, **kwargs) and df(x, *args, **kwargs)
:param args kwargs: Arguments passed as f(x, *args, **kwargs) and df(x, *args, **kwargs)
Examples:
---------
from GPy.models import GradientChecker
N, M, Q = 10, 5, 3
from GPy.models import GradientChecker
N, M, Q = 10, 5, 3
Sinusoid:
Sinusoid:
X = numpy.random.rand(N, Q)
grad = GradientChecker(numpy.sin,numpy.cos,X,'x')
grad.checkgrad(verbose=1)
X = numpy.random.rand(N, Q)
grad = GradientChecker(numpy.sin,numpy.cos,X,'sin_in')
grad.checkgrad(verbose=1)
Using GPy:
Using GPy:
X, Z = numpy.random.randn(N,Q), numpy.random.randn(M,Q)
kern = GPy.kern.linear(Q, ARD=True) + GPy.kern.rbf(Q, ARD=True)
grad = GradientChecker(kern.K,
lambda x: 2*kern.dK_dX(numpy.ones((1,1)), x),
x0 = X.copy(),
names='X')
grad.checkgrad(verbose=1)
grad.randomize()
grad.checkgrad(verbose=1)
X, Z = numpy.random.randn(N,Q), numpy.random.randn(M,Q)
kern = GPy.kern.linear(Q, ARD=True) + GPy.kern.rbf(Q, ARD=True)
grad = GradientChecker(kern.K,
lambda x: kern.dK_dX(numpy.ones((1,1)), x),
x0 = X.copy(),
names=['X_input'])
grad.checkgrad(verbose=1)
grad.randomize()
grad.checkgrad(verbose=1)
"""
Model.__init__(self)
if isinstance(x0, (list, tuple)) and names is None:

57
GPy/kern/kern.py
View file

@ -487,12 +487,11 @@ class kern(Parameterized):
p1.psi1(Z, mu, S, psi11)
Mu, Sigma = p1._crossterm_mu_S(Z, mu, S)
Mu, Sigma = Mu.reshape(NM,self.input_dim), Sigma.reshape(NM,self.input_dim)
p2.psi1(Z, Mu, Sigma, psi12)
eK2 = psi12.reshape(N, M, M)
crossterms = eK2 * (psi11[:, :, None] + psi11[:, None, :])
target += crossterms
#import ipdb;ipdb.set_trace()
else:
raise NotImplementedError, "psi2 cannot be computed for this kernel"
return target
@ -540,15 +539,15 @@ class kern(Parameterized):
# turn around to have rbf in front
p1, p2 = self.parts[i2], self.parts[i1]
ps1, ps2 = self.param_slices[i2], self.param_slices[i1]
N, M = mu.shape[0], Z.shape[0]; NM=N*M
psi11 = np.zeros((N, M))
p1.psi1(Z, mu, S, psi11)
Mu, Sigma = p1._crossterm_mu_S(Z, mu, S)
Mu, Sigma = Mu.reshape(NM,self.input_dim), Sigma.reshape(NM,self.input_dim)
tmp1 = np.zeros_like(target[ps1])
tmp2 = np.zeros_like(target[ps2])
# for n in range(N):
@ -559,7 +558,7 @@ class kern(Parameterized):
# Mu, Sigma= Mu.reshape(N,M,self.input_dim), Sigma.reshape(N,M,self.input_dim)
# p2.dpsi1_dtheta((dL_dpsi2[n:n+1,m:m+1,m_prime:m_prime+1]*(psi11[n:n+1,m_prime:m_prime+1]))[0], Z[m:m+1], Mu[n:n+1,m], Sigma[n:n+1,m], target[ps2])
# p2.dpsi1_dtheta((dL_dpsi2[n:n+1,m:m+1,m_prime:m_prime+1]*(psi11[n:n+1,m:m+1]))[0], Z[m_prime:m_prime+1], Mu[n:n+1, m_prime], Sigma[n:n+1, m_prime], target[ps2])#Z[m_prime:m_prime+1], Mu[n+m:(n+m)+1], Sigma[n+m:(n+m)+1], target[ps2])
if isinstance(p1, RBF) and isinstance(p2, RBF):
psi12 = np.zeros((N, M))
p2.psi1(Z, mu, S, psi12)
@ -571,11 +570,11 @@ class kern(Parameterized):
if isinstance(p1, RBF) and isinstance(p2, Linear):
#import ipdb;ipdb.set_trace()
pass
p2.dpsi1_dtheta((dL_dpsi2*(psi11[:,:,None] + psi11[:,None,:])).reshape(NM,M), Z, Mu, Sigma, tmp2)
target[ps1] += tmp1
target[ps2] += tmp2
target[ps2] += tmp2
else:
raise NotImplementedError, "psi2 cannot be computed for this kernel"
@ -615,17 +614,17 @@ class kern(Parameterized):
psi11 = np.zeros((N, M))
psi12 = np.zeros((NM, M))
#psi12_t = np.zeros((N,M))
p1.psi1(Z, mu, S, psi11)
Mu, Sigma = p1._crossterm_mu_S(Z, mu, S)
Mu, Sigma = Mu.reshape(NM,self.input_dim), Sigma.reshape(NM,self.input_dim)
p2.psi1(Z, Mu, Sigma, psi12)
tmp1 = np.zeros_like(target)
p1.dpsi1_dZ((dL_dpsi2*psi12.reshape(N,M,M)).sum(1), Z, mu, S, tmp1)
p1.dpsi1_dZ((dL_dpsi2*psi12.reshape(N,M,M)).sum(2), Z, mu, S, tmp1)
target += tmp1
#p2.dpsi1_dtheta((dL_dpsi2*(psi11[:,:,None] + psi11[:,None,:])).reshape(NM,M), Z, Mu, Sigma, target)
p2.dpsi1_dZ((dL_dpsi2*(psi11[:,:,None] + psi11[:,None,:])).reshape(NM,M), Z, Mu, Sigma, target)
else:
@ -666,21 +665,21 @@ class kern(Parameterized):
psi11 = np.zeros((N, M))
psi12 = np.zeros((NM, M))
#psi12_t = np.zeros((N,M))
p1.psi1(Z, mu, S, psi11)
Mu, Sigma = p1._crossterm_mu_S(Z, mu, S)
Mu, Sigma = Mu.reshape(NM,self.input_dim), Sigma.reshape(NM,self.input_dim)
p2.psi1(Z, Mu, Sigma, psi12)
p1.dpsi1_dmuS((dL_dpsi2*psi12.reshape(N,M,M)).sum(1), Z, mu, S, target_mu, target_S)
p1.dpsi1_dmuS((dL_dpsi2*psi12.reshape(N,M,M)).sum(2), Z, mu, S, target_mu, target_S)
#p2.dpsi1_dtheta((dL_dpsi2*(psi11[:,:,None] + psi11[:,None,:])).reshape(NM,M), Z, Mu, Sigma, target)
p2.dpsi1_dmuS((dL_dpsi2*(psi11[:,:,None])).sum(1)*2, Z, Mu.reshape(N,M,self.input_dim).sum(1), Sigma.reshape(N,M,self.input_dim).sum(1), target_mu, target_S)
else:
raise NotImplementedError, "psi2 cannot be computed for this kernel"
return target_mu, target_S
def plot(self, x=None, plot_limits=None, which_parts='all', resolution=None, *args, **kwargs):
if which_parts == 'all':
which_parts = [True] * self.num_parts
@ -737,15 +736,16 @@ class kern(Parameterized):
else:
raise NotImplementedError, "Cannot plot a kernel with more than two input dimensions"
from GPy.core.model import Model
from ..core.model import Model
class Kern_check_model(Model):
"""This is a dummy model class used as a base class for checking that the gradients of a given kernel are implemented correctly. It enables checkgradient() to be called independently on a kernel."""
def __init__(self, kernel=None, dL_dK=None, X=None, X2=None):
num_samples = 20
num_samples2 = 10
if kernel==None:
import GPy
kernel = GPy.kern.rbf(1)
del GPy
if X==None:
X = np.random.normal(size=(num_samples, kernel.input_dim))
if dL_dK==None:
@ -753,14 +753,14 @@ class Kern_check_model(Model):
dL_dK = np.ones((X.shape[0], X.shape[0]))
else:
dL_dK = np.ones((X.shape[0], X2.shape[0]))
self.kernel=kernel
self.X = X
self.X2 = X2
self.dL_dK = dL_dK
#self.constrained_indices=[]
#self.constraints=[]
Model.__init__(self)
super(Kern_check_model, self).__init__()
def is_positive_definite(self):
v = np.linalg.eig(self.kernel.K(self.X))[0]
@ -768,7 +768,7 @@ class Kern_check_model(Model):
return False
else:
return True
def _get_params(self):
return self.kernel._get_params()
@ -783,7 +783,7 @@ class Kern_check_model(Model):
def _log_likelihood_gradients(self):
raise NotImplementedError, "This needs to be implemented to use the kern_check_model class."
class Kern_check_dK_dtheta(Kern_check_model):
"""This class allows gradient checks for the gradient of a kernel with respect to parameters. """
def __init__(self, kernel=None, dL_dK=None, X=None, X2=None):
@ -798,7 +798,7 @@ class Kern_check_dKdiag_dtheta(Kern_check_model):
Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=None)
if dL_dK==None:
self.dL_dK = np.ones((self.X.shape[0]))
def log_likelihood(self):
return (self.dL_dK*self.kernel.Kdiag(self.X)).sum()
@ -815,7 +815,7 @@ class Kern_check_dK_dX(Kern_check_model):
def _get_param_names(self):
return ['X_' +str(i) + ','+str(j) for j in range(self.X.shape[1]) for i in range(self.X.shape[0])]
def _get_params(self):
return self.X.flatten()
@ -837,7 +837,7 @@ class Kern_check_dKdiag_dX(Kern_check_model):
def _get_param_names(self):
return ['X_' +str(i) + ','+str(j) for j in range(self.X.shape[1]) for i in range(self.X.shape[0])]
def _get_params(self):
return self.X.flatten()
@ -861,13 +861,15 @@ def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False, X_positive=
if X_positive:
X = abs(X)
if output_ind is not None:
X[:, output_ind] = np.random.randint(kern.parts[0].output_dim, X.shape[0])
assert(output_ind<kern.input_dim)
X[:, output_ind] = np.random.randint(low=0,high=kern.parts[0].output_dim, size=X.shape[0])
if X2==None:
X2 = np.random.randn(20, kern.input_dim)
if X_positive:
X2 = abs(X2)
if output_ind is not None:
X2[:, output_ind] = np.random.randint(kern.parts[0].output_dim, X2.shape[0])
assert(output_ind<kern.input_dim)
X2[:, output_ind] = np.random.randint(low=0, high=kern.parts[0].output_dim, size=X2.shape[0])
if verbose:
print("Checking covariance function is positive definite.")
@ -961,3 +963,4 @@ def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False, X_positive=
return False
return pass_checks
del Model

24
GPy/kern/parts/sympy_helpers.cpp
View file

@ -80,7 +80,7 @@ double ln_diff_erf(double x0, double x1){
else //x0 and x1 non-positive
return log(erfcx(-x0)-erfcx(-x1)*exp(x0*x0 - x1*x1))-x0*x0;
}
// TODO: For all these computations of h things are very efficient at the moment. Need to recode sympykern to allow the precomputations to take place and all the gradients to be computed in one function. Not sure of best way forward for that yet. Neil
double h(double t, double tprime, double d_i, double d_j, double l){
// Compute the h function for the sim covariance.
double half_l_di = 0.5*l*d_i;
@ -170,9 +170,27 @@ double dh_dl(double t, double tprime, double d_i, double d_j, double l){
}
double dh_dt(double t, double tprime, double d_i, double d_j, double l){
return 0.0;
// compute gradient of h function with respect to t.
double diff_t = t - tprime;
double half_l_di = 0.5*l*d_i;
double arg_1 = half_l_di + tprime/l;
double arg_2 = half_l_di - diff_t/l;
double ln_part_1 = ln_diff_erf(arg_1, arg_2);
arg_2 = half_l_di - t/l;
double ln_part_2 = ln_diff_erf(half_l_di, arg_2);
return (d_i*exp(ln_part_2-d_i*t - d_j*tprime) - d_i*exp(ln_part_1-d_i*diff_t) + 2*exp(-d_i*diff_t - pow(half_l_di - diff_t/l, 2))/(sqrt(M_PI)*l) - 2*exp(-d_i*t - d_j*tprime - pow(half_l_di - t/l,2))/(sqrt(M_PI)*l))*exp(half_l_di*half_l_di)/(d_i + d_j);
}
double dh_dtprime(double t, double tprime, double d_i, double d_j, double l){
return 0.0;
// compute gradient of h function with respect to tprime.
double diff_t = t - tprime;
double half_l_di = 0.5*l*d_i;
double arg_1 = half_l_di + tprime/l;
double arg_2 = half_l_di - diff_t/l;
double ln_part_1 = ln_diff_erf(arg_1, arg_2);
arg_2 = half_l_di - t/l;
double ln_part_2 = ln_diff_erf(half_l_di, arg_2);
return (d_i*exp(ln_part_1-d_i*diff_t) + d_j*exp(ln_part_2-d_i*t - d_j*tprime) + (-2*exp(-pow(half_l_di - diff_t/l,2)) + 2*exp(-pow(half_l_di + tprime/l,2)))*exp(-d_i*diff_t)/(sqrt(M_PI)*l))*exp(half_l_di*half_l_di)/(d_i + d_j);
}

1
GPy/models.py
View file

@ -20,3 +20,4 @@ from _models.mrd import MRD#; _mrd = mrd; del mrd
from _models.gradient_checker import GradientChecker#; _gradient_checker = gradient_checker ; del gradient_checker
from _models.gp_multioutput_regression import GPMultioutputRegression#; _gp_multioutput_regression = gp_multioutput_regression ; del gp_multioutput_regression
from _models.sparse_gp_multioutput_regression import SparseGPMultioutputRegression#; _sparse_gp_multioutput_regression = sparse_gp_multioutput_regression ; del sparse_gp_multioutput_regression
from _models.gradient_checker import GradientChecker

2
GPy/testing/kernel_tests.py
View file

@ -36,7 +36,7 @@ class KernelTests(unittest.TestCase):
def test_eq_sympykernel(self):
if SYMPY_AVAILABLE:
kern = GPy.kern.eq_sympy(5, 3)
self.assertTrue(GPy.kern.kern_test(kern, output_ind=3, verbose=verbose))
self.assertTrue(GPy.kern.kern_test(kern, output_ind=4, verbose=verbose))
def test_ode1_eqkernel(self):
if SYMPY_AVAILABLE:

38
README.md
View file

@ -9,6 +9,44 @@ A Gaussian processes framework in Python.
Continuous integration status: ![CI status](https://travis-ci.org/SheffieldML/GPy.png)
Getting started
===============
Installing with pip
-------------------
The simplest way to install GPy is using pip. ubuntu users can do:
sudo apt-get install python-pip
pip install gpy
If you'd like to install from source, or want to contribute to the project (e.g. by sending pull requests via github), read on.
Ubuntu
------
For the most part, the developers are using ubuntu. To install the required packages:
sudo apt-get install python-numpy python-scipy python-matplotlib
clone this git repository and add it to your path:
git clone git@github.com:SheffieldML/GPy.git ~/SheffieldML
echo 'PYTHONPATH=$PYTHONPATH:~/SheffieldML' >> ~/.bashrc
Windows
-------
On windows, we recommend the ![anaconda python distribution](http://continuum.io/downloads). We've also had luck with ![enthought](http://www.enthought.com). git clone or unzip the source to a suitable directory, and add an approptiate PYTHONPATH environment variable.
On windows 7 (and possibly earlier versions) there's a bug in scipy version 0.13 which tries to write very long filenames. Reverting to scipy 0.12 seems to do the trick:
conda install scipy=0.12
OSX
---
Everything appears to work out-of-the box using ![enthought](http://www.enthought.com) on osx Mavericks. Download/clone GPy, and then add GPy to your PYTHONPATH
git clone git@github.com:SheffieldML/GPy.git ~/SheffieldML
echo 'PYTHONPATH=$PYTHONPATH:~/SheffieldML' >> ~/.profile
Compiling documentation:
========================