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

This commit is contained in:
Max Zwiessele 2015-09-07 14:11:06 +01:00
commit 2f0d3b5dcd
14 changed files with 89 additions and 41 deletions

View file

@ -49,7 +49,7 @@ class SparseGP(GP):
else: else:
#inference_method = ?? #inference_method = ??
raise NotImplementedError("what to do what to do?") raise NotImplementedError("what to do what to do?")
print("defaulting to ", inference_method, "for latent function inference") print(("defaulting to ", inference_method, "for latent function inference"))
self.Z = Param('inducing inputs', Z) self.Z = Param('inducing inputs', Z)
self.num_inducing = Z.shape[0] self.num_inducing = Z.shape[0]
@ -160,7 +160,7 @@ class SparseGP(GP):
mu = np.dot(psi1_star, la) # TODO: dimensions? mu = np.dot(psi1_star, la) # TODO: dimensions?
if full_cov: if full_cov:
raise NotImplementedError, "Full covariance for Sparse GP predicted with uncertain inputs not implemented yet." raise NotImplementedError("Full covariance for Sparse GP predicted with uncertain inputs not implemented yet.")
var = np.empty((Xnew.shape[0], la.shape[1], la.shape[1])) var = np.empty((Xnew.shape[0], la.shape[1], la.shape[1]))
di = np.diag_indices(la.shape[1]) di = np.diag_indices(la.shape[1])
else: else:

View file

@ -171,7 +171,7 @@ class Laplace(LatentFunctionInference):
#define the objective function (to be maximised) #define the objective function (to be maximised)
def obj(Ki_f, f): def obj(Ki_f, f):
ll = -0.5*np.sum(np.dot(Ki_f.T, f)) + np.sum(likelihood.logpdf(f, Y, Y_metadata=Y_metadata)) ll = -0.5*np.sum(np.dot(Ki_f.T, f)) + np.sum(likelihood.logpdf(f, Y, Y_metadata=Y_metadata))
print ll print(ll)
if np.isnan(ll): if np.isnan(ll):
import ipdb; ipdb.set_trace() # XXX BREAKPOINT import ipdb; ipdb.set_trace() # XXX BREAKPOINT
return -np.inf return -np.inf

View file

@ -6,7 +6,11 @@ import numpy as np
from ...core.parameterization import Param from ...core.parameterization import Param
from ...core.parameterization.transformations import Logexp from ...core.parameterization.transformations import Logexp
from ...util.config import config # for assesing whether to use cython from ...util.config import config # for assesing whether to use cython
from . import coregionalize_cython try:
from . import coregionalize_cython
config.set('cython', 'working', 'True')
except ImportError:
config.set('cython', 'working', 'False')
class Coregionalize(Kern): class Coregionalize(Kern):
""" """
@ -126,4 +130,3 @@ class Coregionalize(Kern):
def gradients_X_diag(self, dL_dKdiag, X): def gradients_X_diag(self, dL_dKdiag, X):
return np.zeros(X.shape) return np.zeros(X.shape)

View file

@ -88,6 +88,8 @@ class Kern(Parameterized):
return self.psicomp.psicomputations(self, Z, variational_posterior, return_psi2_n=True)[2] return self.psicomp.psicomputations(self, Z, variational_posterior, return_psi2_n=True)[2]
def gradients_X(self, dL_dK, X, X2): def gradients_X(self, dL_dK, X, X2):
raise NotImplementedError raise NotImplementedError
def gradients_X_X2(self, dL_dK, X, X2):
return self.gradients_X(dL_dK, X, X2), self.gradients_X(dL_dK.T, X2, X)
def gradients_XX(self, dL_dK, X, X2): def gradients_XX(self, dL_dK, X, X2):
raise(NotImplementedError, "This is the second derivative of K wrt X and X2, and not implemented for this kernel") raise(NotImplementedError, "This is the second derivative of K wrt X and X2, and not implemented for this kernel")
def gradients_XX_diag(self, dL_dKdiag, X): def gradients_XX_diag(self, dL_dKdiag, X):

View file

@ -19,6 +19,7 @@ class KernCallsViaSlicerMeta(ParametersChangedMeta):
put_clean(dct, 'update_gradients_full', _slice_update_gradients_full) put_clean(dct, 'update_gradients_full', _slice_update_gradients_full)
put_clean(dct, 'update_gradients_diag', _slice_update_gradients_diag) put_clean(dct, 'update_gradients_diag', _slice_update_gradients_diag)
put_clean(dct, 'gradients_X', _slice_gradients_X) put_clean(dct, 'gradients_X', _slice_gradients_X)
put_clean(dct, 'gradients_X_X2', _slice_gradients_X)
put_clean(dct, 'gradients_XX', _slice_gradients_XX) put_clean(dct, 'gradients_XX', _slice_gradients_XX)
put_clean(dct, 'gradients_XX_diag', _slice_gradients_X_diag) put_clean(dct, 'gradients_XX_diag', _slice_gradients_X_diag)
put_clean(dct, 'gradients_X_diag', _slice_gradients_X_diag) put_clean(dct, 'gradients_X_diag', _slice_gradients_X_diag)

View file

@ -5,6 +5,7 @@ from .kern import Kern
from ...core.parameterization import Param from ...core.parameterization import Param
from ...core.parameterization.transformations import Logexp from ...core.parameterization.transformations import Logexp
import numpy as np import numpy as np
from ...util.linalg import tdot
from ...util.caching import Cache_this from ...util.caching import Cache_this
four_over_tau = 2./np.pi four_over_tau = 2./np.pi
@ -40,6 +41,7 @@ class MLP(Kern):
self.link_parameters(self.variance, self.weight_variance, self.bias_variance) self.link_parameters(self.variance, self.weight_variance, self.bias_variance)
@Cache_this(limit=20, ignore_args=())
def K(self, X, X2=None): def K(self, X, X2=None):
if X2 is None: if X2 is None:
X_denom = np.sqrt(self._comp_prod(X)+1.) X_denom = np.sqrt(self._comp_prod(X)+1.)
@ -51,6 +53,7 @@ class MLP(Kern):
XTX = self._comp_prod(X,X2)/X_denom[:,None]/X2_denom[None,:] XTX = self._comp_prod(X,X2)/X_denom[:,None]/X2_denom[None,:]
return self.variance*four_over_tau*np.arcsin(XTX) return self.variance*four_over_tau*np.arcsin(XTX)
@Cache_this(limit=20, ignore_args=())
def Kdiag(self, X): def Kdiag(self, X):
"""Compute the diagonal of the covariance matrix for X.""" """Compute the diagonal of the covariance matrix for X."""
X_prod = self._comp_prod(X) X_prod = self._comp_prod(X)
@ -73,6 +76,10 @@ class MLP(Kern):
"""Derivative of the covariance matrix with respect to X""" """Derivative of the covariance matrix with respect to X"""
return self._comp_grads(dL_dK, X, X2)[3] return self._comp_grads(dL_dK, X, X2)[3]
def gradients_X_X2(self, dL_dK, X, X2):
"""Derivative of the covariance matrix with respect to X"""
return self._comp_grads(dL_dK, X, X2)[3:]
def gradients_X_diag(self, dL_dKdiag, X): def gradients_X_diag(self, dL_dKdiag, X):
"""Gradient of diagonal of covariance with respect to X""" """Gradient of diagonal of covariance with respect to X"""
return self._comp_grads_diag(dL_dKdiag, X)[3] return self._comp_grads_diag(dL_dKdiag, X)[3]

View file

@ -80,8 +80,9 @@ class PSICOMP_GH(PSICOMP):
dL_dkfu = (dL_dpsi1+ 2.*Kfu.dot(dL_dpsi2))*self.weights[i] dL_dkfu = (dL_dpsi1+ 2.*Kfu.dot(dL_dpsi2))*self.weights[i]
kern.update_gradients_full(dL_dkfu, X, Z) kern.update_gradients_full(dL_dkfu, X, Z)
dtheta += kern.gradient dtheta += kern.gradient
dX += kern.gradients_X(dL_dkfu, X, Z) dX_i, dZ_i = kern.gradients_X_X2(dL_dkfu, X, Z)
dZ += kern.gradients_X(dL_dkfu.T, Z, X) dX += dX_i
dZ += dZ_i
dmu += dX dmu += dX
dS += dX*self.locs[i]/(2.*S_sq) dS += dX*self.locs[i]/(2.*S_sq)
kern.gradient[:] = dtheta_old kern.gradient[:] = dtheta_old

View file

@ -1,3 +1,5 @@
#ifndef __APPLE__
#include <omp.h> #include <omp.h>
#endif
void _grad_X(int N, int D, int M, double*X, double* X2, double* tmp, double* grad); void _grad_X(int N, int D, int M, double*X, double* X2, double* tmp, double* grad);
void _lengthscale_grads(int N, int D, int M, double* X, double* X2, double* tmp, double* grad); void _lengthscale_grads(int N, int D, int M, double* X, double* X2, double* tmp, double* grad);

View file

@ -48,7 +48,7 @@ class Gaussian(Likelihood):
def betaY(self,Y,Y_metadata=None): def betaY(self,Y,Y_metadata=None):
#TODO: ~Ricardo this does not live here #TODO: ~Ricardo this does not live here
raise RuntimeError, "Please notify the GPy developers, this should not happen" raise RuntimeError("Please notify the GPy developers, this should not happen")
return Y/self.gaussian_variance(Y_metadata) return Y/self.gaussian_variance(Y_metadata)
def gaussian_variance(self, Y_metadata=None): def gaussian_variance(self, Y_metadata=None):

View file

@ -2,11 +2,20 @@ import numpy as np
import scipy as sp import scipy as sp
from GPy.util import choleskies from GPy.util import choleskies
import GPy import GPy
from ..util.config import config
import unittest
try:
from . import linalg_cython
config.set('cython', 'working', 'True')
except ImportError:
config.set('cython', 'working', 'False')
""" """
These tests make sure that the opure python and cython codes work the same These tests make sure that the opure python and cython codes work the same
""" """
@unittest.skipIf(not config.getboolean('cython', 'working'),"Cython modules have not been built on this machine")
class CythonTestChols(np.testing.TestCase): class CythonTestChols(np.testing.TestCase):
def setUp(self): def setUp(self):
self.flat = np.random.randn(45,5) self.flat = np.random.randn(45,5)
@ -20,6 +29,7 @@ class CythonTestChols(np.testing.TestCase):
A2 = choleskies._triang_to_flat_cython(self.triang) A2 = choleskies._triang_to_flat_cython(self.triang)
np.testing.assert_allclose(A1, A2) np.testing.assert_allclose(A1, A2)
@unittest.skipIf(not config.getboolean('cython', 'working'),"Cython modules have not been built on this machine")
class test_stationary(np.testing.TestCase): class test_stationary(np.testing.TestCase):
def setUp(self): def setUp(self):
self.k = GPy.kern.RBF(10) self.k = GPy.kern.RBF(10)
@ -49,6 +59,7 @@ class test_stationary(np.testing.TestCase):
g2 = self.k._lengthscale_grads_cython(self.dKxz, self.X, self.Z) g2 = self.k._lengthscale_grads_cython(self.dKxz, self.X, self.Z)
np.testing.assert_allclose(g1, g2) np.testing.assert_allclose(g1, g2)
@unittest.skipIf(not config.getboolean('cython', 'working'),"Cython modules have not been built on this machine")
class test_choleskies_backprop(np.testing.TestCase): class test_choleskies_backprop(np.testing.TestCase):
def setUp(self): def setUp(self):
a =np.random.randn(10,12) a =np.random.randn(10,12)
@ -61,10 +72,3 @@ class test_choleskies_backprop(np.testing.TestCase):
r3 = GPy.util.choleskies.choleskies_cython.backprop_gradient_par_c(self.dL, self.L) r3 = GPy.util.choleskies.choleskies_cython.backprop_gradient_par_c(self.dL, self.L)
np.testing.assert_allclose(r1, r2) np.testing.assert_allclose(r1, r2)
np.testing.assert_allclose(r1, r3) np.testing.assert_allclose(r1, r3)

View file

@ -6,9 +6,16 @@ import numpy as np
import GPy import GPy
import sys import sys
from GPy.core.parameterization.param import Param from GPy.core.parameterization.param import Param
from ..util.config import config
verbose = 0 verbose = 0
try:
from . import linalg_cython
config.set('cython', 'working', 'True')
except ImportError:
config.set('cython', 'working', 'False')
class Kern_check_model(GPy.core.Model): class Kern_check_model(GPy.core.Model):
""" """
@ -312,12 +319,12 @@ class KernelGradientTestsContinuous(unittest.TestCase):
k = GPy.kern.LinearFull(self.D, self.D-1) k = GPy.kern.LinearFull(self.D, self.D-1)
k.randomize() k.randomize()
self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
def test_standard_periodic(self): def test_standard_periodic(self):
k = GPy.kern.StdPeriodic(self.D, self.D-1) k = GPy.kern.StdPeriodic(self.D, self.D-1)
k.randomize() k.randomize()
self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
class KernelTestsMiscellaneous(unittest.TestCase): class KernelTestsMiscellaneous(unittest.TestCase):
def setUp(self): def setUp(self):
N, D = 100, 10 N, D = 100, 10
@ -371,6 +378,7 @@ class KernelTestsNonContinuous(unittest.TestCase):
X2 = self.X2[self.X2[:,-1]!=2] X2 = self.X2[self.X2[:,-1]!=2]
self.assertTrue(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1)) self.assertTrue(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1))
@unittest.skipIf(not config.getboolean('cython', 'working'),"Cython modules have not been built on this machine")
class Coregionalize_cython_test(unittest.TestCase): class Coregionalize_cython_test(unittest.TestCase):
""" """
Make sure that the coregionalize kernel work with and without cython enabled Make sure that the coregionalize kernel work with and without cython enabled
@ -438,28 +446,28 @@ class KernelTestsProductWithZeroValues(unittest.TestCase):
"Gradient resulted in NaN") "Gradient resulted in NaN")
class Kernel_Psi_statistics_GradientTests(unittest.TestCase): class Kernel_Psi_statistics_GradientTests(unittest.TestCase):
def setUp(self): def setUp(self):
from GPy.core.parameterization.variational import NormalPosterior from GPy.core.parameterization.variational import NormalPosterior
N,M,Q = 100,20,3 N,M,Q = 100,20,3
X = np.random.randn(N,Q) X = np.random.randn(N,Q)
X_var = np.random.rand(N,Q)+0.01 X_var = np.random.rand(N,Q)+0.01
self.Z = np.random.randn(M,Q) self.Z = np.random.randn(M,Q)
self.qX = NormalPosterior(X, X_var) self.qX = NormalPosterior(X, X_var)
self.w1 = np.random.randn(N) self.w1 = np.random.randn(N)
self.w2 = np.random.randn(N,M) self.w2 = np.random.randn(N,M)
self.w3 = np.random.randn(M,M) self.w3 = np.random.randn(M,M)
self.w3 = self.w3+self.w3.T self.w3 = self.w3+self.w3.T
self.w3n = np.random.randn(N,M,M) self.w3n = np.random.randn(N,M,M)
self.w3n = self.w3n+np.swapaxes(self.w3n, 1,2) self.w3n = self.w3n+np.swapaxes(self.w3n, 1,2)
def test_kernels(self): def test_kernels(self):
from GPy.kern import RBF,Linear from GPy.kern import RBF,Linear
Q = self.Z.shape[1] Q = self.Z.shape[1]
kernels = [RBF(Q,ARD=True), Linear(Q,ARD=True)] kernels = [RBF(Q,ARD=True), Linear(Q,ARD=True)]
for k in kernels: for k in kernels:
k.randomize() k.randomize()
self._test_kernel_param(k) self._test_kernel_param(k)
@ -476,12 +484,12 @@ class Kernel_Psi_statistics_GradientTests(unittest.TestCase):
psi0 = kernel.psi0(self.Z, self.qX) psi0 = kernel.psi0(self.Z, self.qX)
psi1 = kernel.psi1(self.Z, self.qX) psi1 = kernel.psi1(self.Z, self.qX)
if not psi2n: if not psi2n:
psi2 = kernel.psi2(self.Z, self.qX) psi2 = kernel.psi2(self.Z, self.qX)
return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3*psi2).sum() return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3*psi2).sum()
else: else:
psi2 = kernel.psi2n(self.Z, self.qX) psi2 = kernel.psi2n(self.Z, self.qX)
return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3n*psi2).sum() return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3n*psi2).sum()
def df(p): def df(p):
kernel.param_array[:] = p kernel.param_array[:] = p
kernel.update_gradients_expectations(self.w1, self.w2, self.w3 if not psi2n else self.w3n, self.Z, self.qX) kernel.update_gradients_expectations(self.w1, self.w2, self.w3 if not psi2n else self.w3n, self.Z, self.qX)
@ -492,39 +500,39 @@ class Kernel_Psi_statistics_GradientTests(unittest.TestCase):
self.assertTrue(m.checkgrad()) self.assertTrue(m.checkgrad())
def _test_Z(self, kernel, psi2n=False): def _test_Z(self, kernel, psi2n=False):
def f(p): def f(p):
psi0 = kernel.psi0(p, self.qX) psi0 = kernel.psi0(p, self.qX)
psi1 = kernel.psi1(p, self.qX) psi1 = kernel.psi1(p, self.qX)
psi2 = kernel.psi2(p, self.qX) psi2 = kernel.psi2(p, self.qX)
if not psi2n: if not psi2n:
psi2 = kernel.psi2(p, self.qX) psi2 = kernel.psi2(p, self.qX)
return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3*psi2).sum() return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3*psi2).sum()
else: else:
psi2 = kernel.psi2n(p, self.qX) psi2 = kernel.psi2n(p, self.qX)
return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3n*psi2).sum() return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3n*psi2).sum()
def df(p): def df(p):
return kernel.gradients_Z_expectations(self.w1, self.w2, self.w3 if not psi2n else self.w3n, p, self.qX) return kernel.gradients_Z_expectations(self.w1, self.w2, self.w3 if not psi2n else self.w3n, p, self.qX)
from GPy.models import GradientChecker from GPy.models import GradientChecker
m = GradientChecker(f, df, self.Z.copy()) m = GradientChecker(f, df, self.Z.copy())
self.assertTrue(m.checkgrad()) self.assertTrue(m.checkgrad())
def _test_qX(self, kernel, psi2n=False): def _test_qX(self, kernel, psi2n=False):
def f(p): def f(p):
self.qX.param_array[:] = p self.qX.param_array[:] = p
self.qX._trigger_params_changed() self.qX._trigger_params_changed()
psi0 = kernel.psi0(self.Z, self.qX) psi0 = kernel.psi0(self.Z, self.qX)
psi1 = kernel.psi1(self.Z, self.qX) psi1 = kernel.psi1(self.Z, self.qX)
if not psi2n: if not psi2n:
psi2 = kernel.psi2(self.Z, self.qX) psi2 = kernel.psi2(self.Z, self.qX)
return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3*psi2).sum() return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3*psi2).sum()
else: else:
psi2 = kernel.psi2n(self.Z, self.qX) psi2 = kernel.psi2n(self.Z, self.qX)
return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3n*psi2).sum() return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3n*psi2).sum()
def df(p): def df(p):
self.qX.param_array[:] = p self.qX.param_array[:] = p
self.qX._trigger_params_changed() self.qX._trigger_params_changed()

View file

@ -4,8 +4,11 @@
import numpy as np import numpy as np
from . import linalg from . import linalg
from .config import config from .config import config
try:
from . import choleskies_cython from . import choleskies_cython
config.set('cython', 'working', 'True')
except ImportError:
config.set('cython', 'working', 'False')
def safe_root(N): def safe_root(N):
i = np.sqrt(N) i = np.sqrt(N)

View file

@ -8,10 +8,14 @@
import numpy as np import numpy as np
from scipy import linalg from scipy import linalg
from scipy.linalg import lapack, blas from scipy.linalg import lapack, blas
from .config import config from .config import config
import logging import logging
from . import linalg_cython
try:
from . import linalg_cython
config.set('cython', 'working', 'True')
except ImportError:
config.set('cython', 'working', 'False')
def force_F_ordered_symmetric(A): def force_F_ordered_symmetric(A):

View file

@ -12,18 +12,31 @@ version = '0.6.1'
def read(fname): def read(fname):
return open(os.path.join(os.path.dirname(__file__), fname)).read() return open(os.path.join(os.path.dirname(__file__), fname)).read()
#compile_flags = ["-march=native", '-fopenmp', '-O3', ] #Mac OS X Clang doesn't support OpenMP th the current time.
compile_flags = [ '-fopenmp', '-O3', ] #This detects if we are building on a Mac
def ismac():
platform = sys.platform
ismac = False
if platform[:6] == 'darwin':
ismac = True
return ismac
if ismac():
compile_flags = [ '-O3', ]
link_args = ['']
else:
compile_flags = [ '-fopenmp', '-O3', ]
link_args = ['-lgomp']
ext_mods = [Extension(name='GPy.kern._src.stationary_cython', ext_mods = [Extension(name='GPy.kern._src.stationary_cython',
sources=['GPy/kern/_src/stationary_cython.c','GPy/kern/_src/stationary_utils.c'], sources=['GPy/kern/_src/stationary_cython.c','GPy/kern/_src/stationary_utils.c'],
include_dirs=[np.get_include()], include_dirs=[np.get_include()],
extra_compile_args=compile_flags, extra_compile_args=compile_flags,
extra_link_args = ['-lgomp']), extra_link_args = link_args),
Extension(name='GPy.util.choleskies_cython', Extension(name='GPy.util.choleskies_cython',
sources=['GPy/util/choleskies_cython.c'], sources=['GPy/util/choleskies_cython.c'],
include_dirs=[np.get_include()], include_dirs=[np.get_include()],
extra_link_args = ['-lgomp'], extra_link_args = link_args,
extra_compile_args=compile_flags), extra_compile_args=compile_flags),
Extension(name='GPy.util.linalg_cython', Extension(name='GPy.util.linalg_cython',
sources=['GPy/util/linalg_cython.c'], sources=['GPy/util/linalg_cython.c'],