From e7a9b5350601166177f462a992ef00b9e962bdbe Mon Sep 17 00:00:00 2001 From: Mike Croucher Date: Mon, 7 Sep 2015 13:22:51 +0100 Subject: [PATCH 1/4] Automatic fallback to Numpy if Cython modules not available --- GPy/kern/_src/coregionalize.py | 7 ++++-- GPy/testing/cython_tests.py | 18 ++++++++------ GPy/testing/kernel_tests.py | 44 ++++++++++++++++++++-------------- GPy/util/choleskies.py | 7 ++++-- GPy/util/linalg.py | 8 +++++-- 5 files changed, 53 insertions(+), 31 deletions(-) diff --git a/GPy/kern/_src/coregionalize.py b/GPy/kern/_src/coregionalize.py index b8e1e139..9003f2b9 100644 --- a/GPy/kern/_src/coregionalize.py +++ b/GPy/kern/_src/coregionalize.py @@ -6,7 +6,11 @@ import numpy as np from ...core.parameterization import Param from ...core.parameterization.transformations import Logexp 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): """ @@ -126,4 +130,3 @@ class Coregionalize(Kern): def gradients_X_diag(self, dL_dKdiag, X): return np.zeros(X.shape) - diff --git a/GPy/testing/cython_tests.py b/GPy/testing/cython_tests.py index fd257e5b..30e27fbb 100644 --- a/GPy/testing/cython_tests.py +++ b/GPy/testing/cython_tests.py @@ -2,11 +2,20 @@ import numpy as np import scipy as sp from GPy.util import choleskies 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 """ +@unittest.skipIf(not config.getboolean('cython', 'working'),"Cython modules have not been built on this machine") class CythonTestChols(np.testing.TestCase): def setUp(self): self.flat = np.random.randn(45,5) @@ -20,6 +29,7 @@ class CythonTestChols(np.testing.TestCase): A2 = choleskies._triang_to_flat_cython(self.triang) 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): def setUp(self): 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) 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): def setUp(self): 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) np.testing.assert_allclose(r1, r2) np.testing.assert_allclose(r1, r3) - - - - - - - diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 2f33f99b..ec005b6c 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -6,9 +6,16 @@ import numpy as np import GPy import sys from GPy.core.parameterization.param import Param +from ..util.config import config 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): """ @@ -312,12 +319,12 @@ class KernelGradientTestsContinuous(unittest.TestCase): k = GPy.kern.LinearFull(self.D, self.D-1) k.randomize() self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) - + def test_standard_periodic(self): k = GPy.kern.StdPeriodic(self.D, self.D-1) k.randomize() self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) - + class KernelTestsMiscellaneous(unittest.TestCase): def setUp(self): N, D = 100, 10 @@ -371,6 +378,7 @@ class KernelTestsNonContinuous(unittest.TestCase): X2 = self.X2[self.X2[:,-1]!=2] 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): """ Make sure that the coregionalize kernel work with and without cython enabled @@ -438,28 +446,28 @@ class KernelTestsProductWithZeroValues(unittest.TestCase): "Gradient resulted in NaN") class Kernel_Psi_statistics_GradientTests(unittest.TestCase): - + def setUp(self): from GPy.core.parameterization.variational import NormalPosterior N,M,Q = 100,20,3 - + X = np.random.randn(N,Q) X_var = np.random.rand(N,Q)+0.01 self.Z = np.random.randn(M,Q) self.qX = NormalPosterior(X, X_var) - + self.w1 = np.random.randn(N) 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.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) - + def test_kernels(self): from GPy.kern import RBF,Linear Q = self.Z.shape[1] kernels = [RBF(Q,ARD=True), Linear(Q,ARD=True)] - + for k in kernels: k.randomize() self._test_kernel_param(k) @@ -476,12 +484,12 @@ class Kernel_Psi_statistics_GradientTests(unittest.TestCase): psi0 = kernel.psi0(self.Z, self.qX) psi1 = kernel.psi1(self.Z, self.qX) 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() else: psi2 = kernel.psi2n(self.Z, self.qX) return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3n*psi2).sum() - + def df(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) @@ -492,39 +500,39 @@ class Kernel_Psi_statistics_GradientTests(unittest.TestCase): self.assertTrue(m.checkgrad()) def _test_Z(self, kernel, psi2n=False): - + def f(p): psi0 = kernel.psi0(p, self.qX) psi1 = kernel.psi1(p, self.qX) psi2 = kernel.psi2(p, self.qX) 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() else: psi2 = kernel.psi2n(p, self.qX) return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3n*psi2).sum() - + def df(p): 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 m = GradientChecker(f, df, self.Z.copy()) self.assertTrue(m.checkgrad()) - + def _test_qX(self, kernel, psi2n=False): - + def f(p): self.qX.param_array[:] = p self.qX._trigger_params_changed() psi0 = kernel.psi0(self.Z, self.qX) psi1 = kernel.psi1(self.Z, self.qX) 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() else: psi2 = kernel.psi2n(self.Z, self.qX) return (self.w1*psi0).sum() + (self.w2*psi1).sum() + (self.w3n*psi2).sum() - + def df(p): self.qX.param_array[:] = p self.qX._trigger_params_changed() diff --git a/GPy/util/choleskies.py b/GPy/util/choleskies.py index f079cabd..ca055e08 100644 --- a/GPy/util/choleskies.py +++ b/GPy/util/choleskies.py @@ -4,8 +4,11 @@ import numpy as np from . import linalg from .config import config - -from . import choleskies_cython +try: + from . import choleskies_cython + config.set('cython', 'working', 'True') +except ImportError: + config.set('cython', 'working', 'False') def safe_root(N): i = np.sqrt(N) diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index acceeb16..c2f481f0 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -8,10 +8,14 @@ import numpy as np from scipy import linalg from scipy.linalg import lapack, blas - from .config import config 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): From d904045663a70e52dc4d61c1f284c2050fb80dfc Mon Sep 17 00:00:00 2001 From: Mike Croucher Date: Mon, 7 Sep 2015 13:38:11 +0100 Subject: [PATCH 2/4] Fixes Cython compilation on Mac OS X --- GPy/kern/_src/stationary_utils.h | 2 ++ setup.py | 21 +++++++++++++++++---- 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/GPy/kern/_src/stationary_utils.h b/GPy/kern/_src/stationary_utils.h index 9c9745fb..948362ff 100644 --- a/GPy/kern/_src/stationary_utils.h +++ b/GPy/kern/_src/stationary_utils.h @@ -1,3 +1,5 @@ +#ifndef __APPLE__ #include +#endif 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); diff --git a/setup.py b/setup.py index b806a559..ef51cd3e 100644 --- a/setup.py +++ b/setup.py @@ -12,18 +12,31 @@ version = '0.6.1' def read(fname): return open(os.path.join(os.path.dirname(__file__), fname)).read() -#compile_flags = ["-march=native", '-fopenmp', '-O3', ] -compile_flags = [ '-fopenmp', '-O3', ] +#Mac OS X Clang doesn't support OpenMP th the current time. +#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', sources=['GPy/kern/_src/stationary_cython.c','GPy/kern/_src/stationary_utils.c'], include_dirs=[np.get_include()], extra_compile_args=compile_flags, - extra_link_args = ['-lgomp']), + extra_link_args = link_args), Extension(name='GPy.util.choleskies_cython', sources=['GPy/util/choleskies_cython.c'], include_dirs=[np.get_include()], - extra_link_args = ['-lgomp'], + extra_link_args = link_args, extra_compile_args=compile_flags), Extension(name='GPy.util.linalg_cython', sources=['GPy/util/linalg_cython.c'], From e6b1482d2132a687496b9a2b669c9065f75154b0 Mon Sep 17 00:00:00 2001 From: Zhenwen Dai Date: Mon, 7 Sep 2015 14:04:57 +0100 Subject: [PATCH 3/4] speed tuning for mlp kernel and gauss qudrature for psi-statistics --- GPy/kern/_src/kern.py | 2 ++ GPy/kern/_src/kernel_slice_operations.py | 1 + GPy/kern/_src/mlp.py | 7 +++++++ GPy/kern/_src/psi_comp/gaussherm.py | 5 +++-- 4 files changed, 13 insertions(+), 2 deletions(-) diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index ff5d49d3..531137ae 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -88,6 +88,8 @@ class Kern(Parameterized): return self.psicomp.psicomputations(self, Z, variational_posterior, return_psi2_n=True)[2] def gradients_X(self, dL_dK, X, X2): 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): 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): diff --git a/GPy/kern/_src/kernel_slice_operations.py b/GPy/kern/_src/kernel_slice_operations.py index 8c06d0c0..f274eef3 100644 --- a/GPy/kern/_src/kernel_slice_operations.py +++ b/GPy/kern/_src/kernel_slice_operations.py @@ -19,6 +19,7 @@ class KernCallsViaSlicerMeta(ParametersChangedMeta): put_clean(dct, 'update_gradients_full', _slice_update_gradients_full) put_clean(dct, 'update_gradients_diag', _slice_update_gradients_diag) 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_diag', _slice_gradients_X_diag) put_clean(dct, 'gradients_X_diag', _slice_gradients_X_diag) diff --git a/GPy/kern/_src/mlp.py b/GPy/kern/_src/mlp.py index b65fb2e0..c495b77b 100644 --- a/GPy/kern/_src/mlp.py +++ b/GPy/kern/_src/mlp.py @@ -5,6 +5,7 @@ from .kern import Kern from ...core.parameterization import Param from ...core.parameterization.transformations import Logexp import numpy as np +from ...util.linalg import tdot from ...util.caching import Cache_this four_over_tau = 2./np.pi @@ -40,6 +41,7 @@ class MLP(Kern): self.link_parameters(self.variance, self.weight_variance, self.bias_variance) + @Cache_this(limit=20, ignore_args=()) def K(self, X, X2=None): if X2 is None: 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,:] return self.variance*four_over_tau*np.arcsin(XTX) + @Cache_this(limit=20, ignore_args=()) def Kdiag(self, X): """Compute the diagonal of the covariance matrix for X.""" X_prod = self._comp_prod(X) @@ -73,6 +76,10 @@ class MLP(Kern): """Derivative of the covariance matrix with respect to X""" 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): """Gradient of diagonal of covariance with respect to X""" return self._comp_grads_diag(dL_dKdiag, X)[3] diff --git a/GPy/kern/_src/psi_comp/gaussherm.py b/GPy/kern/_src/psi_comp/gaussherm.py index afbca545..8e54e6a0 100644 --- a/GPy/kern/_src/psi_comp/gaussherm.py +++ b/GPy/kern/_src/psi_comp/gaussherm.py @@ -80,8 +80,9 @@ class PSICOMP_GH(PSICOMP): dL_dkfu = (dL_dpsi1+ 2.*Kfu.dot(dL_dpsi2))*self.weights[i] kern.update_gradients_full(dL_dkfu, X, Z) dtheta += kern.gradient - dX += kern.gradients_X(dL_dkfu, X, Z) - dZ += kern.gradients_X(dL_dkfu.T, Z, X) + dX_i, dZ_i = kern.gradients_X_X2(dL_dkfu, X, Z) + dX += dX_i + dZ += dZ_i dmu += dX dS += dX*self.locs[i]/(2.*S_sq) kern.gradient[:] = dtheta_old From c47c63b4088254755176e81259c071d56f190fdd Mon Sep 17 00:00:00 2001 From: Mike Croucher Date: Mon, 7 Sep 2015 14:07:09 +0100 Subject: [PATCH 4/4] Python 3 fixes --- GPy/core/sparse_gp.py | 4 ++-- GPy/inference/latent_function_inference/laplace.py | 2 +- GPy/likelihoods/gaussian.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index e227625d..9ce5c391 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -49,7 +49,7 @@ class SparseGP(GP): else: #inference_method = ?? 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.num_inducing = Z.shape[0] @@ -159,7 +159,7 @@ class SparseGP(GP): mu = np.dot(psi1_star, la) # TODO: dimensions? 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])) di = np.diag_indices(la.shape[1]) else: diff --git a/GPy/inference/latent_function_inference/laplace.py b/GPy/inference/latent_function_inference/laplace.py index 00a2c2b0..2f089141 100644 --- a/GPy/inference/latent_function_inference/laplace.py +++ b/GPy/inference/latent_function_inference/laplace.py @@ -171,7 +171,7 @@ class Laplace(LatentFunctionInference): #define the objective function (to be maximised) 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)) - print ll + print(ll) if np.isnan(ll): import ipdb; ipdb.set_trace() # XXX BREAKPOINT return -np.inf diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py index 424a7f5a..e1299f73 100644 --- a/GPy/likelihoods/gaussian.py +++ b/GPy/likelihoods/gaussian.py @@ -48,7 +48,7 @@ class Gaussian(Likelihood): def betaY(self,Y,Y_metadata=None): #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) def gaussian_variance(self, Y_metadata=None):