From 1d2cbfe44a9f78a1206fdd2366a36e7cad562bd1 Mon Sep 17 00:00:00 2001 From: Daniel Beck Date: Fri, 6 Feb 2015 19:39:46 +1100 Subject: [PATCH 1/5] first attempt --- GPy/kern/_src/prod.py | 15 +++++++++++++-- GPy/testing/kernel_tests.py | 15 +++++++++++++++ 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py index dd9a5fe4..e3776838 100644 --- a/GPy/kern/_src/prod.py +++ b/GPy/kern/_src/prod.py @@ -5,6 +5,7 @@ import numpy as np from kern import CombinationKernel from ...util.caching import Cache_this import itertools +import operator class Prod(CombinationKernel): """ @@ -42,9 +43,19 @@ class Prod(CombinationKernel): return reduce(np.multiply, (p.Kdiag(X) for p in which_parts)) def update_gradients_full(self, dL_dK, X, X2=None): + np.seterr(invalid='raise') k = self.K(X,X2)*dL_dK - for p in self.parts: - p.update_gradients_full(k/p.K(X,X2),X,X2) + try: + for p in self.parts: + p.update_gradients_full(k/p.K(X,X2),X,X2) + except FloatingPointError: + np.seterr(invalid='warn') + print "Gradient warning: falling back to slow version due to zero-valued kernel" + for combination in itertools.combinations(self.parts, len(self.parts) - 1): + prod = reduce(operator.mul, [p.K(X, X2) for p in combination]) + to_update = list(set(self.parts) - set(combination))[0] + to_update.update_gradients_full(dL_dK * prod, X, X2) + def update_gradients_diag(self, dL_dKdiag, X): k = self.Kdiag(X)*dL_dKdiag diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index c1bb9265..387047b6 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -400,12 +400,27 @@ class Coregionalize_weave_test(unittest.TestCase): #reset the weave state for any other tests GPy.util.config.config.set('weave', 'working', 'False') +class KernelTestsProductWithZeroValues(unittest.TestCase): + + def test_zero_valued_kernel(self): + X = np.array([[0,1],[1,0]]) + Y = np.array([[1],[10]]) + lin = GPy.kern.Linear(2) + bias = GPy.kern.Bias(2) + k = lin * bias + #k = lin + m = GPy.models.GPRegression(X, Y, kernel=k) + #m['mul.bias.variance'].constrain_fixed(0) + m.optimize(messages=False) if __name__ == "__main__": print "Running unit tests, please be (very) patient..." unittest.main() + #suite = unittest.TestLoader().loadTestsFromTestCase(KernelTestsProductWithZeroValues) + #unittest.TextTestRunner().run(suite) + # np.random.seed(0) # N0 = 3 # N1 = 9 From 8b4274339ad034aeede9b926ac47bad89ae2f397 Mon Sep 17 00:00:00 2001 From: Daniel Beck Date: Mon, 9 Feb 2015 09:28:53 +1100 Subject: [PATCH 2/5] added decorator that changes numpy invalid op warning to exception --- GPy/kern/_src/prod.py | 20 +++++++++++++++++--- GPy/testing/kernel_tests.py | 7 ++++--- 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py index e3776838..4f9f5ea6 100644 --- a/GPy/kern/_src/prod.py +++ b/GPy/kern/_src/prod.py @@ -7,6 +7,19 @@ from ...util.caching import Cache_this import itertools import operator + +def numpy_invalid_op_as_exception(func): + """ + A decorator that allows catching numpy invalid operations + as exceptions (the default behaviour is raising warnings). + """ + def func_wrapper(*args, **kwargs): + np.seterr(invalid='raise') + func(*args, **kwargs) + np.seterr(invalid='warn') + return func_wrapper + + class Prod(CombinationKernel): """ Computes the product of 2 kernels @@ -42,15 +55,14 @@ class Prod(CombinationKernel): which_parts = self.parts return reduce(np.multiply, (p.Kdiag(X) for p in which_parts)) + @numpy_invalid_op_as_exception def update_gradients_full(self, dL_dK, X, X2=None): - np.seterr(invalid='raise') k = self.K(X,X2)*dL_dK try: for p in self.parts: p.update_gradients_full(k/p.K(X,X2),X,X2) except FloatingPointError: - np.seterr(invalid='warn') - print "Gradient warning: falling back to slow version due to zero-valued kernel" + #print "WARNING: gradient calculation falling back to slow version due to zero-valued kernel" for combination in itertools.combinations(self.parts, len(self.parts) - 1): prod = reduce(operator.mul, [p.K(X, X2) for p in combination]) to_update = list(set(self.parts) - set(combination))[0] @@ -75,3 +87,5 @@ class Prod(CombinationKernel): for p in self.parts: target += p.gradients_X_diag(k/p.Kdiag(X),X) return target + + diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 387047b6..ac6d7ab4 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -408,10 +408,11 @@ class KernelTestsProductWithZeroValues(unittest.TestCase): lin = GPy.kern.Linear(2) bias = GPy.kern.Bias(2) k = lin * bias - #k = lin m = GPy.models.GPRegression(X, Y, kernel=k) - #m['mul.bias.variance'].constrain_fixed(0) - m.optimize(messages=False) + try: + m.optimize() + except np.linalg.LinAlgError: + self.fail("Zero-valued kernel raised exception!") From d6a56a6f0bf234f4c8c9f9f1b595ff5c9305bed0 Mon Sep 17 00:00:00 2001 From: Daniel Beck Date: Mon, 9 Feb 2015 09:35:32 +1100 Subject: [PATCH 3/5] changed operator.mul to np.multiply for consistency --- GPy/kern/_src/prod.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py index 4f9f5ea6..241c2448 100644 --- a/GPy/kern/_src/prod.py +++ b/GPy/kern/_src/prod.py @@ -5,7 +5,6 @@ import numpy as np from kern import CombinationKernel from ...util.caching import Cache_this import itertools -import operator def numpy_invalid_op_as_exception(func): @@ -64,7 +63,7 @@ class Prod(CombinationKernel): except FloatingPointError: #print "WARNING: gradient calculation falling back to slow version due to zero-valued kernel" for combination in itertools.combinations(self.parts, len(self.parts) - 1): - prod = reduce(operator.mul, [p.K(X, X2) for p in combination]) + prod = reduce(np.multiply, [p.K(X, X2) for p in combination]) to_update = list(set(self.parts) - set(combination))[0] to_update.update_gradients_full(dL_dK * prod, X, X2) From fc8705104b05cadc772307536f06f6de803c72bc Mon Sep 17 00:00:00 2001 From: Daniel Beck Date: Mon, 9 Feb 2015 09:41:21 +1100 Subject: [PATCH 4/5] a cleaner test --- GPy/kern/_src/prod.py | 1 - GPy/testing/kernel_tests.py | 8 +++----- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py index 241c2448..5e4c0d29 100644 --- a/GPy/kern/_src/prod.py +++ b/GPy/kern/_src/prod.py @@ -61,7 +61,6 @@ class Prod(CombinationKernel): for p in self.parts: p.update_gradients_full(k/p.K(X,X2),X,X2) except FloatingPointError: - #print "WARNING: gradient calculation falling back to slow version due to zero-valued kernel" for combination in itertools.combinations(self.parts, len(self.parts) - 1): prod = reduce(np.multiply, [p.K(X, X2) for p in combination]) to_update = list(set(self.parts) - set(combination))[0] diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index ac6d7ab4..f9d90607 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -408,11 +408,9 @@ class KernelTestsProductWithZeroValues(unittest.TestCase): lin = GPy.kern.Linear(2) bias = GPy.kern.Bias(2) k = lin * bias - m = GPy.models.GPRegression(X, Y, kernel=k) - try: - m.optimize() - except np.linalg.LinAlgError: - self.fail("Zero-valued kernel raised exception!") + k.update_gradients_full(1, X) + self.assertFalse(np.isnan(k['linear.variances'].gradient), + "Gradient resulted in NaN") From 98c743d157f2954226b0ef5b6d3d1817f28e67f6 Mon Sep 17 00:00:00 2001 From: Daniel Beck Date: Mon, 9 Feb 2015 10:02:26 +1100 Subject: [PATCH 5/5] test + code change in gradients_X --- GPy/kern/_src/prod.py | 15 +++++++++++---- GPy/testing/kernel_tests.py | 22 ++++++++++++---------- 2 files changed, 23 insertions(+), 14 deletions(-) diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py index 5e4c0d29..a3b49973 100644 --- a/GPy/kern/_src/prod.py +++ b/GPy/kern/_src/prod.py @@ -14,8 +14,9 @@ def numpy_invalid_op_as_exception(func): """ def func_wrapper(*args, **kwargs): np.seterr(invalid='raise') - func(*args, **kwargs) + result = func(*args, **kwargs) np.seterr(invalid='warn') + return result return func_wrapper @@ -66,17 +67,23 @@ class Prod(CombinationKernel): to_update = list(set(self.parts) - set(combination))[0] to_update.update_gradients_full(dL_dK * prod, X, X2) - def update_gradients_diag(self, dL_dKdiag, X): k = self.Kdiag(X)*dL_dKdiag for p in self.parts: p.update_gradients_diag(k/p.Kdiag(X),X) + @numpy_invalid_op_as_exception def gradients_X(self, dL_dK, X, X2=None): target = np.zeros(X.shape) k = self.K(X,X2)*dL_dK - for p in self.parts: - target += p.gradients_X(k/p.K(X,X2),X,X2) + try: + for p in self.parts: + target += p.gradients_X(k/p.K(X,X2),X,X2) + except FloatingPointError: + for combination in itertools.combinations(self.parts, len(self.parts) - 1): + prod = reduce(np.multiply, [p.K(X, X2) for p in combination]) + to_update = list(set(self.parts) - set(combination))[0] + target += to_update.gradients_X(dL_dK * prod, X, X2) return target def gradients_X_diag(self, dL_dKdiag, X): diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index f9d90607..415cc7eb 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -400,25 +400,27 @@ class Coregionalize_weave_test(unittest.TestCase): #reset the weave state for any other tests GPy.util.config.config.set('weave', 'working', 'False') + class KernelTestsProductWithZeroValues(unittest.TestCase): - def test_zero_valued_kernel(self): - X = np.array([[0,1],[1,0]]) - Y = np.array([[1],[10]]) - lin = GPy.kern.Linear(2) - bias = GPy.kern.Bias(2) - k = lin * bias - k.update_gradients_full(1, X) - self.assertFalse(np.isnan(k['linear.variances'].gradient), + def setUp(self): + self.X = np.array([[0,1],[1,0]]) + self.k = GPy.kern.Linear(2) * GPy.kern.Bias(2) + + def test_zero_valued_kernel_full(self): + self.k.update_gradients_full(1, self.X) + self.assertFalse(np.isnan(self.k['linear.variances'].gradient), "Gradient resulted in NaN") + def test_zero_valued_kernel_gradients_X(self): + target = self.k.gradients_X(1, self.X) + self.assertFalse(np.any(np.isnan(target)), + "Gradient resulted in NaN") if __name__ == "__main__": print "Running unit tests, please be (very) patient..." unittest.main() - #suite = unittest.TestLoader().loadTestsFromTestCase(KernelTestsProductWithZeroValues) - #unittest.TextTestRunner().run(suite) # np.random.seed(0) # N0 = 3