From 4757265b240924af09ac0032984c5c8d85f4d805 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Fri, 14 Mar 2014 15:23:32 +0000 Subject: [PATCH] Added a hack fix as suggested by max, zeroing any negative values (should really be numerically negative values on diagonal) --- GPy/kern/_src/independent_outputs.py | 3 ++- GPy/kern/_src/stationary.py | 28 +++++++++++++++------------- GPy/testing/kernel_tests.py | 19 +++++++++++++------ 3 files changed, 30 insertions(+), 20 deletions(-) diff --git a/GPy/kern/_src/independent_outputs.py b/GPy/kern/_src/independent_outputs.py index 1848bf6a..85cb95bc 100644 --- a/GPy/kern/_src/independent_outputs.py +++ b/GPy/kern/_src/independent_outputs.py @@ -8,7 +8,7 @@ import itertools def index_to_slices(index): """ - take a numpy array of integers (index) and return a nested list of slices such that the slices describe the start, stop points for each integer in the index. + take a numpy array of integers (index) and return a nested list of slices such that the slices describe the start, stop points for each integer in the index. e.g. >>> index = np.asarray([0,0,0,1,1,1,2,2,2]) @@ -40,6 +40,7 @@ class IndependentOutputs(CombinationKernel): The index of the functions is given by the last column in the input X the rest of the columns of X are passed to the underlying kernel for computation (in blocks). + Kern is wrapped with a slicer metaclass """ def __init__(self, kern, index_dim=-1, name='independ'): assert isinstance(index_dim, int), "IndependentOutputs kernel is only defined with one input dimension being the indeces" diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index df7ba058..a9e837a9 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -15,21 +15,21 @@ class Stationary(Kern): """ Stationary kernels (covariance functions). - Stationary covariance fucntion depend only on r, where r is defined as + Stationary covariance fucntion depend only on r, where r is defined as r = \sqrt{ \sum_{q=1}^Q (x_q - x'_q)^2 } - The covariance function k(x, x' can then be written k(r). + The covariance function k(x, x' can then be written k(r). In this implementation, r is scaled by the lengthscales parameter(s): - r = \sqrt{ \sum_{q=1}^Q \frac{(x_q - x'_q)^2}{\ell_q^2} }. - + r = \sqrt{ \sum_{q=1}^Q \frac{(x_q - x'_q)^2}{\ell_q^2} }. + By default, there's only one lengthscale: seaprate lengthscales for each - dimension can be enables by setting ARD=True. + dimension can be enables by setting ARD=True. To implement a stationary covariance function using this class, one need - only define the covariance function k(r), and it derivative. + only define the covariance function k(r), and it derivative. ... def K_of_r(self, r): @@ -37,10 +37,10 @@ class Stationary(Kern): def dK_dr(self, r): return bar - The lengthscale(s) and variance parameters are added to the structure automatically. - + The lengthscale(s) and variance parameters are added to the structure automatically. + """ - + def __init__(self, input_dim, variance, lengthscale, ARD, active_dims, name): super(Stationary, self).__init__(input_dim, active_dims, name) self.ARD = ARD @@ -57,7 +57,7 @@ class Stationary(Kern): if lengthscale.size != input_dim: lengthscale = np.ones(input_dim)*lengthscale else: - lengthscale = np.ones(self.input_dim) + lengthscale = np.ones(self.input_dim) self.lengthscale = Param('lengthscale', lengthscale, Logexp()) self.variance = Param('variance', variance, Logexp()) assert self.variance.size==1 @@ -95,7 +95,9 @@ class Stationary(Kern): #X2, = self._slice_X(X2) X1sq = np.sum(np.square(X),1) X2sq = np.sum(np.square(X2),1) - return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])) + r2 = -2.*np.dot(X, X2.T) + X1sq[:,None] + X2sq[None,:] + r2[r2<0] = 0. # A bit hacky + return np.sqrt(r2) @Cache_this(limit=5, ignore_args=()) def _scaled_dist(self, X, X2=None): @@ -133,7 +135,7 @@ class Stationary(Kern): if self.ARD: #rinv = self._inv_dis# this is rather high memory? Should we loop instead?t(X, X2) #d = X[:, None, :] - X2[None, :, :] - #x_xl3 = np.square(d) + #x_xl3 = np.square(d) #self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3 tmp = dL_dr*self._inv_dist(X, X2) if X2 is None: X2 = X @@ -247,7 +249,7 @@ class Matern52(Stationary): .. math:: - k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) + k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) """ def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, active_dims=None, name='Mat52'): super(Matern52, self).__init__(input_dim, variance, lengthscale, ARD, active_dims, name) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index b057f8ef..3b17a3d0 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -120,6 +120,8 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb if verbose: print("Checking covariance function is positive definite.") + #if isinstance(kern, GPy.kern.IndependentOutputs): + #import ipdb; ipdb.set_trace() # XXX BREAKPOINT result = Kern_check_model(kern, X=X).is_positive_semi_definite() if result and verbose: print("Check passed.") @@ -306,17 +308,22 @@ class KernelTestsNonContinuous(unittest.TestCase): D = self.D self.X = np.random.randn(N,D) self.X2 = np.random.randn(N1,D) - self.X_block = np.zeros((N+N1, D+D+1)) + #self.X_block = np.zeros((N+N1, D+D+1)) + #self.X_block[0:N, 0:D] = self.X + #self.X_block[N:N+N1, D:D+D] = self.X2 + #self.X_block[0:N, -1] = 0 + #self.X_block[N:N+N1, -1] = 1 + self.X_block = np.zeros((N+N1, D+1)) self.X_block[0:N, 0:D] = self.X - self.X_block[N:N+N1, D:D+D] = self.X2 - self.X_block[0:N, -1] = 1 - self.X_block[N:N+1, -1] = 2 + self.X_block[N:N+N1, 0:D] = self.X2 + self.X_block[0:N, -1] = 0 + self.X_block[N:N+N1, -1] = 1 self.X_block = self.X_block[self.X_block.argsort(0)[:, -1], :] - + def test_IndependentOutputs(self): k = GPy.kern.RBF(self.D) kern = GPy.kern.IndependentOutputs(k, -1) - self.assertTrue(check_kernel_gradient_functions(kern, X=self.X_block, X2=self.X_block, verbose=verbose)) + self.assertTrue(check_kernel_gradient_functions(kern, X=self.X_block, verbose=verbose)) if __name__ == "__main__": print "Running unit tests, please be (very) patient..."