diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index a826cdf7..00a80c7b 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -52,10 +52,7 @@ class SparseGP(GP): self.parameters_changed() def has_uncertain_inputs(self): - if isinstance(self.X, VariationalPosterior): - return True - else: - return False + return isinstance(self.X, VariationalPosterior) def parameters_changed(self): self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, self.Z, self.likelihood, self.Y) diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index 3b8e391b..a2a83929 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -162,6 +162,8 @@ class Matern52(Stationary): k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) \ \ \ \ \ \\text{ where } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} } """ + def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Mat52'): + super(Matern52, self).__init__(input_dim, variance, lengthscale, ARD, name) def K_of_r(self, r): return self.variance*(1+np.sqrt(5.)*r+5./3*r**2)*np.exp(-np.sqrt(5.)*r) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 40cd66dd..e5985145 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -4,6 +4,7 @@ import unittest import numpy as np import GPy +import sys verbose = True @@ -14,106 +15,223 @@ except ImportError: SYMPY_AVAILABLE=False -class KernelTests(unittest.TestCase): - def test_kerneltie(self): - K = GPy.kern.rbf(5, ARD=True) - K.tie_params('.*[01]') - K.constrain_fixed('2') - X = np.random.rand(5,5) - Y = np.ones((5,1)) - m = GPy.models.GPRegression(X,Y,K) - self.assertTrue(m.checkgrad()) +class Kern_check_model(GPy.core.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 + checkgrad() to be called independently on a kernel. + """ + def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): + GPy.core.Model.__init__(self, 'kernel_test_model') + if kernel==None: + kernel = GPy.kern.RBF(1) + if X is None: + X = np.random.randn(20, kernel.input_dim) + if dL_dK is None: + if X2 is None: + dL_dK = np.ones((X.shape[0], X.shape[0])) + else: + dL_dK = np.ones((X.shape[0], X2.shape[0])) - def test_rbfkernel(self): - kern = GPy.kern.rbf(5) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) + self.kernel = kernel + self.X = GPy.core.parameterization.Param('X',X) + self.X2 = X2 + self.dL_dK = dL_dK - def test_rbf_sympykernel(self): - if SYMPY_AVAILABLE: - kern = GPy.kern.rbf_sympy(5) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) + def is_positive_definite(self): + v = np.linalg.eig(self.kernel.K(self.X))[0] + if any(v<-10*sys.float_info.epsilon): + return False + else: + return True - def test_eq_sympykernel(self): - if SYMPY_AVAILABLE: - kern = GPy.kern.eq_sympy(5, 3) - self.assertTrue(GPy.kern.kern_test(kern, output_ind=4, verbose=verbose)) + def log_likelihood(self): + return np.sum(self.dL_dK*self.kernel.K(self.X, self.X2)) - def test_ode1_eqkernel(self): - if SYMPY_AVAILABLE: - kern = GPy.kern.ode1_eq(3) - self.assertTrue(GPy.kern.kern_test(kern, output_ind=1, verbose=verbose, X_positive=True)) +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): + Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) + self.add_parameter(self.kernel) - def test_rbf_invkernel(self): - kern = GPy.kern.rbf_inv(5) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) + def parameters_changed(self): + return self.kernel.update_gradients_full(self.dL_dK, self.X, self.X2) - def test_Matern32kernel(self): - kern = GPy.kern.Matern32(5) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) - def test_Matern52kernel(self): - kern = GPy.kern.Matern52(5) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) +class Kern_check_dKdiag_dtheta(Kern_check_model): + """ + This class allows gradient checks of the gradient of the diagonal of a + kernel with respect to the parameters. + """ + def __init__(self, kernel=None, dL_dK=None, X=None): + Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=None) + self.add_parameter(self.kernel) - def test_linearkernel(self): - kern = GPy.kern.linear(5) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) + def parameters_changed(self): + self.kernel.update_gradients_diag(self.dL_dK, self.X) - def test_periodic_exponentialkernel(self): - kern = GPy.kern.periodic_exponential(1) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) + def log_likelihood(self): + return (np.diag(self.dL_dK)*self.kernel.Kdiag(self.X)).sum() - def test_periodic_Matern32kernel(self): - kern = GPy.kern.periodic_Matern32(1) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) + def parameters_changed(self): + return self.kernel.update_gradients_diag(np.diag(self.dL_dK), self.X) - def test_periodic_Matern52kernel(self): - kern = GPy.kern.periodic_Matern52(1) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) +class Kern_check_dK_dX(Kern_check_model): + """This class allows gradient checks for the gradient of a kernel with respect to X. """ + def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): + Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) + self.add_parameter(self.X) - def test_rational_quadratickernel(self): - kern = GPy.kern.rational_quadratic(1) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) + def parameters_changed(self): + self.X.gradient = self.kernel.gradients_X(self.dL_dK, self.X, self.X2) - def test_gibbskernel(self): - kern = GPy.kern.gibbs(5, mapping=GPy.mappings.Linear(5, 1)) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) +class Kern_check_dKdiag_dX(Kern_check_dK_dX): + """This class allows gradient checks for the gradient of a kernel diagonal with respect to X. """ + def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): + Kern_check_dK_dX.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=None) - def test_heterokernel(self): - kern = GPy.kern.hetero(5, mapping=GPy.mappings.Linear(5, 1), transform=GPy.core.transformations.logexp()) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) + def log_likelihood(self): + return (np.diag(self.dL_dK)*self.kernel.Kdiag(self.X)).sum() - def test_mlpkernel(self): - kern = GPy.kern.mlp(5) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) + def parameters_changed(self): + self.X.gradient = self.kernel.gradients_X_diag(self.dL_dK, self.X) - def test_polykernel(self): - kern = GPy.kern.poly(5, degree=4) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) +def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False): + """ + This function runs on kernels to check the correctness of their + implementation. It checks that the covariance function is positive definite + for a randomly generated data set. - def test_fixedkernel(self): - """ - Fixed effect kernel test - """ - X = np.random.rand(30, 4) - K = np.dot(X, X.T) - kernel = GPy.kern.fixed(4, K) - kern = GPy.kern.poly(5, degree=4) - self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) + :param kern: the kernel to be tested. + :type kern: GPy.kern.Kernpart + :param X: X input values to test the covariance function. + :type X: ndarray + :param X2: X2 input values to test the covariance function. + :type X2: ndarray - # def test_coregionalization(self): - # X1 = np.random.rand(50,1)*8 - # X2 = np.random.rand(30,1)*5 - # index = np.vstack((np.zeros_like(X1),np.ones_like(X2))) - # X = np.hstack((np.vstack((X1,X2)),index)) - # Y1 = np.sin(X1) + np.random.randn(*X1.shape)*0.05 - # Y2 = np.sin(X2) + np.random.randn(*X2.shape)*0.05 + 2. - # Y = np.vstack((Y1,Y2)) + """ + pass_checks = True + if X==None: + X = np.random.randn(10, kern.input_dim) + if output_ind is not None: + X[:, output_ind] = np.random.randint(kern.output_dim, X.shape[0]) + if X2==None: + X2 = np.random.randn(20, kern.input_dim) + if output_ind is not None: + X2[:, output_ind] = np.random.randint(kern.output_dim, X2.shape[0]) + + if verbose: + print("Checking covariance function is positive definite.") + result = Kern_check_model(kern, X=X).is_positive_definite() + if result and verbose: + print("Check passed.") + if not result: + print("Positive definite check failed for " + kern.name + " covariance function.") + pass_checks = False + return False + + if verbose: + print("Checking gradients of K(X, X) wrt theta.") + result = Kern_check_dK_dtheta(kern, X=X, X2=None).checkgrad(verbose=verbose) + if result and verbose: + print("Check passed.") + if not result: + print("Gradient of K(X, X) wrt theta failed for " + kern.name + " covariance function. Gradient values as follows:") + Kern_check_dK_dtheta(kern, X=X, X2=None).checkgrad(verbose=True) + pass_checks = False + return False + + if verbose: + print("Checking gradients of K(X, X2) wrt theta.") + result = Kern_check_dK_dtheta(kern, X=X, X2=X2).checkgrad(verbose=verbose) + if result and verbose: + print("Check passed.") + if not result: + print("Gradient of K(X, X) wrt theta failed for " + kern.name + " covariance function. Gradient values as follows:") + Kern_check_dK_dtheta(kern, X=X, X2=X2).checkgrad(verbose=True) + pass_checks = False + return False + + if verbose: + print("Checking gradients of Kdiag(X) wrt theta.") + result = Kern_check_dKdiag_dtheta(kern, X=X).checkgrad(verbose=verbose) + if result and verbose: + print("Check passed.") + if not result: + print("Gradient of Kdiag(X) wrt theta failed for " + kern.name + " covariance function. Gradient values as follows:") + Kern_check_dKdiag_dtheta(kern, X=X).checkgrad(verbose=True) + pass_checks = False + return False + + if verbose: + print("Checking gradients of K(X, X) wrt X.") + try: + result = Kern_check_dK_dX(kern, X=X, X2=None).checkgrad(verbose=verbose) + except NotImplementedError: + result=True + if verbose: + print("gradients_X not implemented for " + kern.name) + if result and verbose: + print("Check passed.") + if not result: + print("Gradient of K(X, X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:") + Kern_check_dK_dX(kern, X=X, X2=None).checkgrad(verbose=True) + pass_checks = False + return False + + if verbose: + print("Checking gradients of K(X, X2) wrt X.") + try: + result = Kern_check_dK_dX(kern, X=X, X2=X2).checkgrad(verbose=verbose) + except NotImplementedError: + result=True + if verbose: + print("gradients_X not implemented for " + kern.name) + if result and verbose: + print("Check passed.") + if not result: + print("Gradient of K(X, X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:") + Kern_check_dK_dX(kern, X=X, X2=X2).checkgrad(verbose=True) + pass_checks = False + return False + + if verbose: + print("Checking gradients of Kdiag(X) wrt X.") + try: + result = Kern_check_dKdiag_dX(kern, X=X).checkgrad(verbose=verbose) + except NotImplementedError: + result=True + if verbose: + print("gradients_X not implemented for " + kern.name) + if result and verbose: + print("Check passed.") + if not result: + print("Gradient of Kdiag(X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:") + Kern_check_dKdiag_dX(kern, X=X).checkgrad(verbose=True) + pass_checks = False + return False + + return pass_checks + + +class KernelTestsContinuous(unittest.TestCase): + def setUp(self): + self.X = np.random.randn(100,2) + self.X2 = np.random.randn(110,2) + + def test_Matern32(self): + k = GPy.kern.Matern32(2) + self.assertTrue(kern_test(k, X=self.X, X2=self.X2, verbose=verbose)) + + def test_Matern52(self): + k = GPy.kern.Matern52(2) + self.assertTrue(kern_test(k, X=self.X, X2=self.X2, verbose=verbose)) + + #TODO: turn off grad checkingwrt X for indexed kernels liek coregionalize - # k1 = GPy.kern.rbf(1) + GPy.kern.bias(1) - # k2 = GPy.kern.coregionalize(2,1) - # kern = k1**k2 - # self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) if __name__ == "__main__": diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index 22b4f86c..97fe2446 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -75,7 +75,7 @@ def jitchol(A, maxtries=5): raise linalg.LinAlgError, "not pd: non-positive diagonal elements" jitter = diagA.mean() * 1e-6 - return jitchol(A+np.eye(A.shape[0])*jitter, maxtries-1) + return jitchol(A + np.eye(A.shape[0])*jitter, maxtries-1) #def jitchol(A, maxtries=5): # A = np.ascontiguousarray(A)