From 1d98d0a718159d37cddd58c106eb802251ffc8f4 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 11 Mar 2013 11:06:08 +0000 Subject: [PATCH 01/18] FIxed a transpose bug in sparse_GPLVM --- GPy/models/sparse_GPLVM.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/models/sparse_GPLVM.py b/GPy/models/sparse_GPLVM.py index fe7c1c43..542fbe0e 100644 --- a/GPy/models/sparse_GPLVM.py +++ b/GPy/models/sparse_GPLVM.py @@ -43,7 +43,7 @@ class sparse_GPLVM(sparse_GP_regression, GPLVM): def dL_dX(self): dL_dX = self.kern.dKdiag_dX(self.dL_dpsi0,self.X) - dL_dX += self.kern.dK_dX(self.dL_dpsi1,self.X,self.Z) + dL_dX += self.kern.dK_dX(self.dL_dpsi1.T,self.X,self.Z) return dL_dX From cd7539a4292aa2abe3c3faff6933bc1ac8c8e500 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 11 Mar 2013 11:08:53 +0000 Subject: [PATCH 02/18] added simple gplvm_tests --- GPy/testing/gplvm_tests.py | 47 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 GPy/testing/gplvm_tests.py diff --git a/GPy/testing/gplvm_tests.py b/GPy/testing/gplvm_tests.py new file mode 100644 index 00000000..51828768 --- /dev/null +++ b/GPy/testing/gplvm_tests.py @@ -0,0 +1,47 @@ +# Copyright (c) 2012, Nicolo Fusi +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import unittest +import numpy as np +import GPy + +class GPLVMTests(unittest.TestCase): + def test_bias_kern(self): + N, M, Q, D = 10, 3, 2, 4 + X = np.random.rand(N, Q) + k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(N),K,D).T + k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) + m = GPy.models.GPLVM(Y, Q, kernel = k) + m.ensure_default_constraints() + m.randomize() + self.assertTrue(m.checkgrad()) + + def test_linear_kern(self): + N, M, Q, D = 10, 3, 2, 4 + X = np.random.rand(N, Q) + k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(N),K,D).T + k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001) + m = GPy.models.GPLVM(Y, Q, kernel = k) + m.ensure_default_constraints() + m.randomize() + self.assertTrue(m.checkgrad()) + + def test_rbf_kern(self): + N, M, Q, D = 10, 3, 2, 4 + X = np.random.rand(N, Q) + k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(N),K,D).T + k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) + m = GPy.models.GPLVM(Y, Q, kernel = k) + m.ensure_default_constraints() + m.randomize() + self.assertTrue(m.checkgrad()) + +if __name__ == "__main__": + print "Running unit tests, please be (very) patient..." + unittest.main() From 7d4e568d7b10e36207982ae5a78c35777519c3c9 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 11 Mar 2013 11:11:42 +0000 Subject: [PATCH 03/18] added sparse_gplvm_tests -- they fail --- GPy/testing/sparse_gplvm_tests.py | 47 +++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 GPy/testing/sparse_gplvm_tests.py diff --git a/GPy/testing/sparse_gplvm_tests.py b/GPy/testing/sparse_gplvm_tests.py new file mode 100644 index 00000000..72bb5bf8 --- /dev/null +++ b/GPy/testing/sparse_gplvm_tests.py @@ -0,0 +1,47 @@ +# Copyright (c) 2012, Nicolo Fusi, James Hensman +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import unittest +import numpy as np +import GPy + +class sparse_GPLVMTests(unittest.TestCase): + def test_bias_kern(self): + N, M, Q, D = 10, 3, 2, 4 + X = np.random.rand(N, Q) + k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(N),K,D).T + k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) + m = GPy.models.sparse_GPLVM(Y, Q, kernel = k, M=M) + m.ensure_default_constraints() + m.randomize() + self.assertTrue(m.checkgrad()) + + def test_linear_kern(self): + N, M, Q, D = 10, 3, 2, 4 + X = np.random.rand(N, Q) + k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(N),K,D).T + k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001) + m = GPy.models.sparse_GPLVM(Y, Q, kernel = k, M=M) + m.ensure_default_constraints() + m.randomize() + self.assertTrue(m.checkgrad()) + + def test_rbf_kern(self): + N, M, Q, D = 10, 3, 2, 4 + X = np.random.rand(N, Q) + k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(N),K,D).T + k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) + m = GPy.models.sparse_GPLVM(Y, Q, kernel = k, M=M) + m.ensure_default_constraints() + m.randomize() + self.assertTrue(m.checkgrad()) + +if __name__ == "__main__": + print "Running unit tests, please be (very) patient..." + unittest.main() From c561867b3cbd9833ca81725d02af2e79e3474382 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 11 Mar 2013 11:19:20 +0000 Subject: [PATCH 04/18] skipping a test known to fail (linear sparse) --- GPy/testing/sparse_gplvm_tests.py | 1 + 1 file changed, 1 insertion(+) diff --git a/GPy/testing/sparse_gplvm_tests.py b/GPy/testing/sparse_gplvm_tests.py index 72bb5bf8..35fa4fcf 100644 --- a/GPy/testing/sparse_gplvm_tests.py +++ b/GPy/testing/sparse_gplvm_tests.py @@ -18,6 +18,7 @@ class sparse_GPLVMTests(unittest.TestCase): m.randomize() self.assertTrue(m.checkgrad()) + @unittest.skip('linear kernels do not have dKdiag_dX') def test_linear_kern(self): N, M, Q, D = 10, 3, 2, 4 X = np.random.rand(N, Q) From 79ad72c46aedfea11bd60b82bfc667b8776c872e Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 11 Mar 2013 11:51:21 +0000 Subject: [PATCH 05/18] added the outline of a tutorial on 'interacting with models' --- doc/tuto_interacting_with_models.rst | 60 ++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 doc/tuto_interacting_with_models.rst diff --git a/doc/tuto_interacting_with_models.rst b/doc/tuto_interacting_with_models.rst new file mode 100644 index 00000000..370ffd95 --- /dev/null +++ b/doc/tuto_interacting_with_models.rst @@ -0,0 +1,60 @@ +************************************* +Interacting with models +************************************* + +The GPy model class has a set of features which are designed to make it simple to explore the parameter space of the model. By default, the scipy optimisers are used to fit GPy models (via model.optimize()), for which we provide mechanisms for 'free' optimisation: GPy can ensure that naturally positive parameters (such as variances) remain positive. But these mechanisms are much more powerful than simple reparameterisation, as we shall see. + +All of the examples included in GPy return an instance of a model class. We'll use GPy.examples.?? as an example:: + + import pylab as pb + pb.ion() + import GPy + m = GPy.examples.?? + +Examining the model using print +=============================== +To see the current state of the model parameters, and the model's (marginal) likelihood just print the model:: + print m + +?? output + +Getting the model's likelihood and gradients +=========================================== +foobar + +Setting and fetching parameters by name +======================================= +foobar + +Constraining and optimising the model +===================================== +A simple task in GPy is to ensure that the models' variances remain positive during optimisation. the models class has a function called constrain_positive(), which accepts a regex string as above. To constrain the models' variance to be positive:: + m.constrain_positive('variance') + print m + +Now we see that the variance of the model is constrained to be postive. GPy handles the effective change of gradients: see how m.objective_gradients has changed approriately + + +For convenience, we also provide a catch all function which ensures that anything which appears to require positivity is constrianed appropriately:: + m.ensure_default_constraints() + + +Fixing parameters +================= + + +Tying Parameters +================ + +Bounding parameters +=================== + + +Further Reading +=============== +All of the mechansiams for dealing with parameters are baked right into GPy.core.model, from which all of the classes in GPy.models inherrit. To learn how to construct your own model, you might want to read ??link?? creating_new_models. + +By deafult, GPy uses the tnc optimizer (from scipy.optimize.tnc). To use other optimisers, and to control the setting of those optimisers, as well as other funky features like automated restarts and diagnostics, you can read the optimization tutorial ??link??. + + + From 12d6f5056bf5f5e590fcf1c293fabc3c87f24ebf Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Mon, 11 Mar 2013 12:15:59 +0000 Subject: [PATCH 06/18] removed keyname partial --- GPy/kern/Matern32.py | 18 ++++----- GPy/kern/Matern52.py | 18 ++++----- GPy/kern/bias.py | 38 +++++++++--------- GPy/kern/coregionalise.py | 26 ++++++------ GPy/kern/exponential.py | 18 ++++----- GPy/kern/kern.py | 68 ++++++++++++++++---------------- GPy/kern/kernpart.py | 18 ++++----- GPy/kern/linear.py | 50 +++++++++++------------ GPy/kern/periodic_Matern32.py | 16 ++++---- GPy/kern/periodic_Matern52.py | 18 ++++----- GPy/kern/periodic_exponential.py | 18 ++++----- GPy/kern/product.py | 24 +++++------ GPy/kern/product_orthogonal.py | 24 +++++------ GPy/kern/rbf.py | 60 ++++++++++++++-------------- GPy/kern/symmetric.py | 24 +++++------ GPy/kern/white.py | 30 +++++++------- GPy/models/GP.py | 2 +- 17 files changed, 235 insertions(+), 235 deletions(-) diff --git a/GPy/kern/Matern32.py b/GPy/kern/Matern32.py index c175009d..9503361d 100644 --- a/GPy/kern/Matern32.py +++ b/GPy/kern/Matern32.py @@ -76,7 +76,7 @@ class Matern32(kernpart): """Compute the diagonal of the covariance matrix associated to X.""" np.add(target,self.variance,target) - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters.""" if X2 is None: X2 = X dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1)) @@ -84,29 +84,29 @@ class Matern32(kernpart): invdist = 1./np.where(dist!=0.,dist,np.inf) dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3 #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] - target[0] += np.sum(dvar*partial) + target[0] += np.sum(dvar*dL_dK) if self.ARD == True: dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] #dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None] - target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + target[1:] += (dl*dL_dK[:,:,None]).sum(0).sum(0) else: dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist)) * dist2M.sum(-1)*invdist #dl = self.variance*dvar*dist2M.sum(-1)*invdist - target[1] += np.sum(dl*partial) + target[1] += np.sum(dl*dL_dK) - def dKdiag_dtheta(self,partial,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): """derivative of the diagonal of the covariance matrix with respect to the parameters.""" - target[0] += np.sum(partial) + target[0] += np.sum(dL_dKdiag) - def dK_dX(self,partial,X,X2,target): + def dK_dX(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to X.""" if X2 is None: X2 = X dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None] ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf) dK_dX = - np.transpose(3*self.variance*dist*np.exp(-np.sqrt(3)*dist)*ddist_dX,(1,0,2)) - target += np.sum(dK_dX*partial.T[:,:,None],0) + target += np.sum(dK_dX*dL_dK.T[:,:,None],0) - def dKdiag_dX(self,partial,X,target): + def dKdiag_dX(self,dL_dKdiag,X,target): pass def Gram_matrix(self,F,F1,F2,lower,upper): diff --git a/GPy/kern/Matern52.py b/GPy/kern/Matern52.py index 26caad1c..377526d5 100644 --- a/GPy/kern/Matern52.py +++ b/GPy/kern/Matern52.py @@ -74,7 +74,7 @@ class Matern52(kernpart): """Compute the diagonal of the covariance matrix associated to X.""" np.add(target,self.variance,target) - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters.""" if X2 is None: X2 = X dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1)) @@ -82,29 +82,29 @@ class Matern52(kernpart): dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3 dvar = (1+np.sqrt(5.)*dist+5./3*dist**2)*np.exp(-np.sqrt(5.)*dist) dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] - target[0] += np.sum(dvar*partial) + target[0] += np.sum(dvar*dL_dK) if self.ARD: dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] - target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + target[1:] += (dl*dL_dK[:,:,None]).sum(0).sum(0) else: dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist)) * dist2M.sum(-1)*invdist #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist)) * dist2M.sum(-1)*invdist - target[1] += np.sum(dl*partial) + target[1] += np.sum(dl*dL_dKdiag) - def dKdiag_dtheta(self,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): """derivative of the diagonal of the covariance matrix with respect to the parameters.""" - target[0] += np.sum(partial) + target[0] += np.sum(dL_dKdiag) - def dK_dX(self,partial,X,X2,target): + def dK_dX(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to X.""" if X2 is None: X2 = X dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None] ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf) dK_dX = - np.transpose(self.variance*5./3*dist*(1+np.sqrt(5)*dist)*np.exp(-np.sqrt(5)*dist)*ddist_dX,(1,0,2)) - target += np.sum(dK_dX*partial.T[:,:,None],0) + target += np.sum(dK_dX*dL_dK.T[:,:,None],0) - def dKdiag_dX(self,partial,X,target): + def dKdiag_dX(self,dL_dKdiag,X,target): pass def Gram_matrix(self,F,F1,F2,F3,lower,upper): diff --git a/GPy/kern/bias.py b/GPy/kern/bias.py index 91594e4c..07679abd 100644 --- a/GPy/kern/bias.py +++ b/GPy/kern/bias.py @@ -35,16 +35,17 @@ class bias(kernpart): def Kdiag(self,X,target): target += self.variance - def dK_dtheta(self,partial,X,X2,target): - target += partial.sum() + def dK_dtheta(self,dL_dKdiag,X,X2,target): + target += dL_dKdiag.sum() - def dKdiag_dtheta(self,partial,X,target): - target += partial.sum() - def dK_dX(self, partial,X, X2, target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): + target += dL_dKdiag.sum() + + def dK_dX(self, dL_dK,X, X2, target): pass - def dKdiag_dX(self,partial,X,target): + def dKdiag_dX(self,dL_dKdiag,X,target): pass #---------------------------------------# @@ -60,30 +61,29 @@ class bias(kernpart): def psi2(self, Z, mu, S, target): target += self.variance**2 - def dpsi0_dtheta(self, partial, Z, mu, S, target): - target += partial.sum() + def dpsi0_dtheta(self, dL_dpsi0, Z, mu, S, target): + target += dL_dpsi0.sum() - def dpsi1_dtheta(self, partial, Z, mu, S, target): - target += partial.sum() + def dpsi1_dtheta(self, dL_dpsi1, Z, mu, S, target): + target += dL_dpsi1.sum() - def dpsi2_dtheta(self, partial, Z, mu, S, target): - target += 2.*self.variance*partial.sum() + def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, target): + target += 2.*self.variance*dL_dpsi2.sum() - - def dpsi0_dZ(self, partial, Z, mu, S, target): + def dpsi0_dZ(self, dL_dpsi0, Z, mu, S, target): pass - def dpsi0_dmuS(self, partial, Z, mu, S, target_mu, target_S): + def dpsi0_dmuS(self, dL_dpsi0, Z, mu, S, target_mu, target_S): pass - def dpsi1_dZ(self, partial, Z, mu, S, target): + def dpsi1_dZ(self, dL_dpsi1, Z, mu, S, target): pass - def dpsi1_dmuS(self, partial, Z, mu, S, target_mu, target_S): + def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S, target_mu, target_S): pass - def dpsi2_dZ(self, partial, Z, mu, S, target): + def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, target): pass - def dpsi2_dmuS(self, partial, Z, mu, S, target_mu, target_S): + def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S): pass diff --git a/GPy/kern/coregionalise.py b/GPy/kern/coregionalise.py index 2a9177d5..a76bb31e 100644 --- a/GPy/kern/coregionalise.py +++ b/GPy/kern/coregionalise.py @@ -53,7 +53,7 @@ class coregionalise(kernpart): def Kdiag(self,index,target): target += np.diag(self.B)[np.asarray(index,dtype=np.int).flatten()] - def dK_dtheta(self,partial,index,index2,target): + def dK_dtheta(self,dL_dK,index,index2,target): index = np.asarray(index,dtype=np.int) if index2 is None: index2 = index @@ -62,28 +62,28 @@ class coregionalise(kernpart): ii,jj = np.meshgrid(index,index2) ii,jj = ii.T, jj.T - partial_small = np.zeros_like(self.B) + dL_dK_small = np.zeros_like(self.B) for i in range(self.Nout): for j in range(self.Nout): - tmp = np.sum(partial[(ii==i)*(jj==j)]) - partial_small[i,j] = tmp + tmp = np.sum(dL_dK[(ii==i)*(jj==j)]) + dL_dK_small[i,j] = tmp - dkappa = np.diag(partial_small) - partial_small += partial_small.T - dW = (self.W[:,None,:]*partial_small[:,:,None]).sum(0) + dkappa = np.diag(dL_dK_small) + dL_dK_small += dL_dK_small.T + dW = (self.W[:,None,:]*dL_dK_small[:,:,None]).sum(0) target += np.hstack([dW.flatten(),dkappa]) - def dKdiag_dtheta(self,partial,index,target): + def dKdiag_dtheta(self,dL_dKdiag,index,target): index = np.asarray(index,dtype=np.int).flatten() - partial_small = np.zeros(self.Nout) + dL_dKdiag_small = np.zeros(self.Nout) for i in range(self.Nout): - partial_small[i] += np.sum(partial[index==i]) - dW = 2.*self.W*partial_small[:,None] - dkappa = partial_small + dL_dKdiag_small[i] += np.sum(dL_dKdiag[index==i]) + dW = 2.*self.W*dL_dKdiag_small[:,None] + dkappa = dL_dKdiag_small target += np.hstack([dW.flatten(),dkappa]) - def dK_dX(self,partial,X,X2,target): + def dK_dX(self,dL_dK,X,X2,target): pass diff --git a/GPy/kern/exponential.py b/GPy/kern/exponential.py index 366ddf3b..9e50712b 100644 --- a/GPy/kern/exponential.py +++ b/GPy/kern/exponential.py @@ -74,35 +74,35 @@ class exponential(kernpart): """Compute the diagonal of the covariance matrix associated to X.""" np.add(target,self.variance,target) - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters.""" if X2 is None: X2 = X dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1)) invdist = 1./np.where(dist!=0.,dist,np.inf) dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3 dvar = np.exp(-dist) - target[0] += np.sum(dvar*partial) + target[0] += np.sum(dvar*dL_dK) if self.ARD == True: dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None] - target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + target[1:] += (dl*dL_dK[:,:,None]).sum(0).sum(0) else: dl = self.variance*dvar*dist2M.sum(-1)*invdist - target[1] += np.sum(dl*partial) + target[1] += np.sum(dl*dL_dK) - def dKdiag_dtheta(self,partial,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): """derivative of the diagonal of the covariance matrix with respect to the parameters.""" #NB: derivative of diagonal elements wrt lengthscale is 0 - target[0] += np.sum(partial) + target[0] += np.sum(dL_dKdiag) - def dK_dX(self,partial,X,X2,target): + def dK_dX(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to X.""" if X2 is None: X2 = X dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None] ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf) dK_dX = - np.transpose(self.variance*np.exp(-dist)*ddist_dX,(1,0,2)) - target += np.sum(dK_dX*partial.T[:,:,None],0) + target += np.sum(dK_dX*dL_dK.T[:,:,None],0) - def dKdiag_dX(self,partial,X,target): + def dKdiag_dX(self,dL_dKdiag,X,target): pass def Gram_matrix(self,F,F1,lower,upper): diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py index 99ad46ea..c1f5eca9 100644 --- a/GPy/kern/kern.py +++ b/GPy/kern/kern.py @@ -271,10 +271,10 @@ class kern(parameterised): [p.K(X[s1,i_s],X2[s2,i_s],target=target[s1,s2]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)] return target - def dK_dtheta(self,partial,X,X2=None,slices1=None,slices2=None): + def dK_dtheta(self,dL_dK,X,X2=None,slices1=None,slices2=None): """ - :param partial: An array of partial derivaties, dL_dK - :type partial: Np.ndarray (N x M) + :param dL_dK: An array of dL_dK derivaties, dL_dK + :type dL_dK: Np.ndarray (N x M) :param X: Observed data inputs :type X: np.ndarray (N x D) :param X2: Observed dara inputs (optional, defaults to X) @@ -288,16 +288,16 @@ class kern(parameterised): if X2 is None: X2 = X target = np.zeros(self.Nparam) - [p.dK_dtheta(partial[s1,s2],X[s1,i_s],X2[s2,i_s],target[ps]) for p,i_s,ps,s1,s2 in zip(self.parts, self.input_slices, self.param_slices, slices1, slices2)] + [p.dK_dtheta(dL_dK[s1,s2],X[s1,i_s],X2[s2,i_s],target[ps]) for p,i_s,ps,s1,s2 in zip(self.parts, self.input_slices, self.param_slices, slices1, slices2)] return self._transform_gradients(target) - def dK_dX(self,partial,X,X2=None,slices1=None,slices2=None): + def dK_dX(self,dL_dK,X,X2=None,slices1=None,slices2=None): if X2 is None: X2 = X slices1, slices2 = self._process_slices(slices1,slices2) target = np.zeros_like(X) - [p.dK_dX(partial[s1,s2],X[s1,i_s],X2[s2,i_s],target[s1,i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)] + [p.dK_dX(dL_dK[s1,s2],X[s1,i_s],X2[s2,i_s],target[s1,i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)] return target def Kdiag(self,X,slices=None): @@ -307,20 +307,20 @@ class kern(parameterised): [p.Kdiag(X[s,i_s],target=target[s]) for p,i_s,s in zip(self.parts,self.input_slices,slices)] return target - def dKdiag_dtheta(self,partial,X,slices=None): + def dKdiag_dtheta(self,dL_dKdiag,X,slices=None): assert X.shape[1]==self.D - assert len(partial.shape)==1 - assert partial.size==X.shape[0] + assert len(dL_dKdiag.shape)==1 + assert dL_dKdiag.size==X.shape[0] slices = self._process_slices(slices,False) target = np.zeros(self.Nparam) - [p.dKdiag_dtheta(partial[s],X[s,i_s],target[ps]) for p,i_s,s,ps in zip(self.parts,self.input_slices,slices,self.param_slices)] + [p.dKdiag_dtheta(dL_dKdiag[s],X[s,i_s],target[ps]) for p,i_s,s,ps in zip(self.parts,self.input_slices,slices,self.param_slices)] return self._transform_gradients(target) - def dKdiag_dX(self, partial, X, slices=None): + def dKdiag_dX(self, dL_dKdiag, X, slices=None): assert X.shape[1]==self.D slices = self._process_slices(slices,False) target = np.zeros_like(X) - [p.dKdiag_dX(partial[s],X[s,i_s],target[s,i_s]) for p,i_s,s in zip(self.parts,self.input_slices,slices)] + [p.dKdiag_dX(dL_dKdiag[s],X[s,i_s],target[s,i_s]) for p,i_s,s in zip(self.parts,self.input_slices,slices)] return target def psi0(self,Z,mu,S,slices=None): @@ -329,16 +329,16 @@ class kern(parameterised): [p.psi0(Z,mu[s],S[s],target[s]) for p,s in zip(self.parts,slices)] return target - def dpsi0_dtheta(self,partial,Z,mu,S,slices=None): + def dpsi0_dtheta(self,dL_dpsi0,Z,mu,S,slices=None): slices = self._process_slices(slices,False) target = np.zeros(self.Nparam) - [p.dpsi0_dtheta(partial[s],Z,mu[s],S[s],target[ps]) for p,ps,s in zip(self.parts, self.param_slices,slices)] + [p.dpsi0_dtheta(dL_dpsi0[s],Z,mu[s],S[s],target[ps]) for p,ps,s in zip(self.parts, self.param_slices,slices)] return self._transform_gradients(target) - def dpsi0_dmuS(self,partial,Z,mu,S,slices=None): + def dpsi0_dmuS(self,dL_dpsi0,Z,mu,S,slices=None): slices = self._process_slices(slices,False) target_mu,target_S = np.zeros_like(mu),np.zeros_like(S) - [p.dpsi0_dmuS(partial,Z,mu[s],S[s],target_mu[s],target_S[s]) for p,s in zip(self.parts,slices)] + [p.dpsi0_dmuS(dL_dpsi0,Z,mu[s],S[s],target_mu[s],target_S[s]) for p,s in zip(self.parts,slices)] return target_mu,target_S def psi1(self,Z,mu,S,slices1=None,slices2=None): @@ -348,25 +348,25 @@ class kern(parameterised): [p.psi1(Z[s2],mu[s1],S[s1],target[s1,s2]) for p,s1,s2 in zip(self.parts,slices1,slices2)] return target - def dpsi1_dtheta(self,partial,Z,mu,S,slices1=None,slices2=None): + def dpsi1_dtheta(self,dL_dpsi1,Z,mu,S,slices1=None,slices2=None): """N,M,(Ntheta)""" slices1, slices2 = self._process_slices(slices1,slices2) target = np.zeros((self.Nparam)) - [p.dpsi1_dtheta(partial[s2,s1],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[ps]) for p,ps,s1,s2,i_s in zip(self.parts, self.param_slices,slices1,slices2,self.input_slices)] + [p.dpsi1_dtheta(dL_dpsi1[s2,s1],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[ps]) for p,ps,s1,s2,i_s in zip(self.parts, self.param_slices,slices1,slices2,self.input_slices)] return self._transform_gradients(target) - def dpsi1_dZ(self,partial,Z,mu,S,slices1=None,slices2=None): + def dpsi1_dZ(self,dL_dpsi1,Z,mu,S,slices1=None,slices2=None): """N,M,Q""" slices1, slices2 = self._process_slices(slices1,slices2) target = np.zeros_like(Z) - [p.dpsi1_dZ(partial[s2,s1],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[s2,i_s]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)] + [p.dpsi1_dZ(dL_dpsi1[s2,s1],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[s2,i_s]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)] return target - def dpsi1_dmuS(self,partial,Z,mu,S,slices1=None,slices2=None): + def dpsi1_dmuS(self,dL_dpsi1,Z,mu,S,slices1=None,slices2=None): """return shapes are N,M,Q""" slices1, slices2 = self._process_slices(slices1,slices2) target_mu, target_S = np.zeros((2,mu.shape[0],mu.shape[1])) - [p.dpsi1_dmuS(partial[s2,s1],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target_mu[s1,i_s],target_S[s1,i_s]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)] + [p.dpsi1_dmuS(dL_dpsi1[s2,s1],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target_mu[s1,i_s],target_S[s1,i_s]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)] return target_mu, target_S def psi2(self,Z,mu,S,slices1=None,slices2=None): @@ -416,11 +416,11 @@ class kern(parameterised): return target + crossterms - def dpsi2_dtheta(self,partial,partial1,Z,mu,S,slices1=None,slices2=None): + def dpsi2_dtheta(self,dL_dpsi2,partial1,Z,mu,S,slices1=None,slices2=None): """Returns shape (N,M,M,Ntheta)""" slices1, slices2 = self._process_slices(slices1,slices2) target = np.zeros(self.Nparam) - [p.dpsi2_dtheta(partial[s1,s2,s2],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[ps]) for p,i_s,s1,s2,ps in zip(self.parts,self.input_slices,slices1,slices2,self.param_slices)] + [p.dpsi2_dtheta(dL_dpsi2[s1,s2,s2],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[ps]) for p,i_s,s1,s2,ps in zip(self.parts,self.input_slices,slices1,slices2,self.param_slices)] #compute the "cross" terms #TODO: better looping @@ -434,11 +434,11 @@ class kern(parameterised): pass #rbf X bias elif p1.name=='bias' and p2.name=='rbf': - p2.dpsi1_dtheta(partial.sum(1)*p1.variance,Z,mu,S,target[ps2]) - p1.dpsi1_dtheta(partial.sum(1)*p2._psi1,Z,mu,S,target[ps1]) + p2.dpsi1_dtheta(dL_dpsi2.sum(1)*p1.variance,Z,mu,S,target[ps2]) + p1.dpsi1_dtheta(dL_dpsi2.sum(1)*p2._psi1,Z,mu,S,target[ps1]) elif p2.name=='bias' and p1.name=='rbf': - p1.dpsi1_dtheta(partial.sum(1)*p2.variance,Z,mu,S,target[ps1]) - p2.dpsi1_dtheta(partial.sum(1)*p1._psi1,Z,mu,S,target[ps2]) + p1.dpsi1_dtheta(dL_dpsi2.sum(1)*p2.variance,Z,mu,S,target[ps1]) + p2.dpsi1_dtheta(dL_dpsi2.sum(1)*p1._psi1,Z,mu,S,target[ps2]) #rbf X linear elif p1.name=='linear' and p2.name=='rbf': raise NotImplementedError #TODO @@ -469,10 +469,10 @@ class kern(parameterised): # target += (partial.sum(0)[:,:,None] * (tmp[:, None] + tmp[:,:,None]).sum(0)).sum(0).sum(0) return self._transform_gradients(target) - def dpsi2_dZ(self,partial,Z,mu,S,slices1=None,slices2=None): + def dpsi2_dZ(self,dL_dpsi2,Z,mu,S,slices1=None,slices2=None): slices1, slices2 = self._process_slices(slices1,slices2) target = np.zeros_like(Z) - [p.dpsi2_dZ(partial[s1,s2,s2],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[s2,i_s]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)] + [p.dpsi2_dZ(dL_dpsi2[s1,s2,s2],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target[s2,i_s]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)] #compute the "cross" terms #TODO: slices (need to iterate around the input slices also...) @@ -482,9 +482,9 @@ class kern(parameterised): pass #rbf X bias elif p1.name=='bias' and p2.name=='rbf': - target += p2.dpsi1_dX(partial.sum(1)*p1.variance,Z,mu,S) + target += p2.dpsi1_dX(dL_dpsi2.sum(1)*p1.variance,Z,mu,S) elif p2.name=='bias' and p1.name=='rbf': - target += p1.dpsi1_dZ(partial.sum(2)*p2.variance,Z,mu,S) + target += p1.dpsi1_dZ(dL_dpsi2.sum(2)*p2.variance,Z,mu,S) #rbf X linear elif p1.name=='linear' and p2.name=='rbf': raise NotImplementedError #TODO @@ -496,11 +496,11 @@ class kern(parameterised): return target - def dpsi2_dmuS(self,partial,Z,mu,S,slices1=None,slices2=None): + def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,slices1=None,slices2=None): """return shapes are N,M,M,Q""" slices1, slices2 = self._process_slices(slices1,slices2) target_mu, target_S = np.zeros((2,mu.shape[0],mu.shape[1])) - [p.dpsi2_dmuS(partial[s1,s2,s2],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target_mu[s1,i_s],target_S[s1,i_s]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)] + [p.dpsi2_dmuS(dL_dpsi2[s1,s2,s2],Z[s2,i_s],mu[s1,i_s],S[s1,i_s],target_mu[s1,i_s],target_S[s1,i_s]) for p,i_s,s1,s2 in zip(self.parts,self.input_slices,slices1,slices2)] #TODO: there are some extra terms to compute here! return target_mu, target_S diff --git a/GPy/kern/kernpart.py b/GPy/kern/kernpart.py index 3a5486de..30a1cc3d 100644 --- a/GPy/kern/kernpart.py +++ b/GPy/kern/kernpart.py @@ -26,31 +26,31 @@ class kernpart(object): raise NotImplementedError def Kdiag(self,X,target): raise NotImplementedError - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): raise NotImplementedError - def dKdiag_dtheta(self,partial,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): raise NotImplementedError def psi0(self,Z,mu,S,target): raise NotImplementedError - def dpsi0_dtheta(self,partial,Z,mu,S,target): + def dpsi0_dtheta(self,dL_dpsi0,Z,mu,S,target): raise NotImplementedError - def dpsi0_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi0_dmuS(self,dL_dpsi0,Z,mu,S,target_mu,target_S): raise NotImplementedError def psi1(self,Z,mu,S,target): raise NotImplementedError def dpsi1_dtheta(self,Z,mu,S,target): raise NotImplementedError - def dpsi1_dZ(self,partial,Z,mu,S,target): + def dpsi1_dZ(self,dL_dpsi1,Z,mu,S,target): raise NotImplementedError - def dpsi1_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi1_dmuS(self,dL_dpsi1,Z,mu,S,target_mu,target_S): raise NotImplementedError def psi2(self,Z,mu,S,target): raise NotImplementedError - def dpsi2_dZ(self,partial,Z,mu,S,target): + def dpsi2_dZ(self,dL_dpsi2,Z,mu,S,target): raise NotImplementedError - def dpsi2_dtheta(self,partial,Z,mu,S,target): + def dpsi2_dtheta(self,dL_dpsi2,Z,mu,S,target): raise NotImplementedError - def dpsi2_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,target_mu,target_S): raise NotImplementedError def dK_dX(self,X,X2,target): raise NotImplementedError diff --git a/GPy/kern/linear.py b/GPy/kern/linear.py index df2fed46..7d817f62 100644 --- a/GPy/kern/linear.py +++ b/GPy/kern/linear.py @@ -73,16 +73,16 @@ class linear(kernpart): def Kdiag(self,X,target): np.add(target,np.sum(self.variances*np.square(X),-1),target) - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): if self.ARD: product = X[:,None,:]*X2[None,:,:] - target += (partial[:,:,None]*product).sum(0).sum(0) + target += (dL_dK[:,:,None]*product).sum(0).sum(0) else: self._K_computations(X, X2) - target += np.sum(self._dot_product*partial) + target += np.sum(self._dot_product*dL_dK) - def dK_dX(self,partial,X,X2,target): - target += (((X2[:, None, :] * self.variances)) * partial[:,:, None]).sum(0) + def dK_dX(self,dL_dK,X,X2,target): + target += (((X2[:, None, :] * self.variances)) * dL_dK[:,:, None]).sum(0) #---------------------------------------# # PSI statistics # @@ -92,40 +92,40 @@ class linear(kernpart): self._psi_computations(Z,mu,S) target += np.sum(self.variances*self.mu2_S,1) - def dKdiag_dtheta(self,partial, X, target): - tmp = partial[:,None]*X**2 + def dKdiag_dtheta(self,dL_dKdiag, X, target): + tmp = dL_dKdiag[:,None]*X**2 if self.ARD: target += tmp.sum(0) else: target += tmp.sum() - def dpsi0_dtheta(self,partial,Z,mu,S,target): + def dpsi0_dtheta(self,dL_dpsi0,Z,mu,S,target): self._psi_computations(Z,mu,S) - tmp = partial[:, None] * self.mu2_S + tmp = dL_dpsi0[:, None] * self.mu2_S if self.ARD: target += tmp.sum(0) else: target += tmp.sum() - def dpsi0_dmuS(self,partial, Z,mu,S,target_mu,target_S): - target_mu += partial[:, None] * (2.0*mu*self.variances) - target_S += partial[:, None] * self.variances + def dpsi0_dmuS(self,dL_dpsi0, Z,mu,S,target_mu,target_S): + target_mu += dL_dpsi0[:, None] * (2.0*mu*self.variances) + target_S += dL_dpsi0[:, None] * self.variances def psi1(self,Z,mu,S,target): """the variance, it does nothing""" self.K(mu,Z,target) - def dpsi1_dtheta(self,partial,Z,mu,S,target): + def dpsi1_dtheta(self,dL_dpsi1,Z,mu,S,target): """the variance, it does nothing""" - self.dK_dtheta(partial,mu,Z,target) + self.dK_dtheta(dL_dpsi1,mu,Z,target) - def dpsi1_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi1_dmuS(self,dL_dpsi1,Z,mu,S,target_mu,target_S): """Do nothing for S, it does not affect psi1""" self._psi_computations(Z,mu,S) - target_mu += (partial.T[:,:, None]*(Z*self.variances)).sum(1) + target_mu += (dL_dpsi1.T[:,:, None]*(Z*self.variances)).sum(1) - def dpsi1_dZ(self,partial,Z,mu,S,target): - self.dK_dX(partial.T,Z,mu,target) + def dpsi1_dZ(self,dL_dpsi1,Z,mu,S,target): + self.dK_dX(dL_dpsi1.T,Z,mu,target) def psi2(self,Z,mu,S,target): """ @@ -135,25 +135,25 @@ class linear(kernpart): psi2 = self.ZZ*np.square(self.variances)*self.mu2_S[:, None, None, :] target += psi2.sum(-1) - def dpsi2_dtheta(self,partial,Z,mu,S,target): + def dpsi2_dtheta(self,dL_dpsi2,Z,mu,S,target): self._psi_computations(Z,mu,S) - tmp = (partial[:,:,:,None]*(2.*self.ZZ*self.mu2_S[:,None,None,:]*self.variances)) + tmp = (dL_dpsi2[:,:,:,None]*(2.*self.ZZ*self.mu2_S[:,None,None,:]*self.variances)) if self.ARD: target += tmp.sum(0).sum(0).sum(0) else: target += tmp.sum() - def dpsi2_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,target_mu,target_S): """Think N,M,M,Q """ self._psi_computations(Z,mu,S) tmp = self.ZZ*np.square(self.variances) # M,M,Q - target_mu += (partial[:,:,:,None]*tmp*2.*mu[:,None,None,:]).sum(1).sum(1) - target_S += (partial[:,:,:,None]*tmp).sum(1).sum(1) + target_mu += (dL_dpsi2[:,:,:,None]*tmp*2.*mu[:,None,None,:]).sum(1).sum(1) + target_S += (dL_dpsi2[:,:,:,None]*tmp).sum(1).sum(1) - def dpsi2_dZ(self,partial,Z,mu,S,target): + def dpsi2_dZ(self,dL_dpsi2,Z,mu,S,target): self._psi_computations(Z,mu,S) mu2_S = np.sum(self.mu2_S,0)# Q, - target += (partial[:,:,:,None] * (self.mu2_S[:,None,None,:]*(Z*np.square(self.variances)[None,:])[None,None,:,:])).sum(0).sum(1) + target += (dL_dpsi2[:,:,:,None] * (self.mu2_S[:,None,None,:]*(Z*np.square(self.variances)[None,:])[None,None,:,:])).sum(0).sum(1) #---------------------------------------# # Precomputations # diff --git a/GPy/kern/periodic_Matern32.py b/GPy/kern/periodic_Matern32.py index be1148c4..898dff7b 100644 --- a/GPy/kern/periodic_Matern32.py +++ b/GPy/kern/periodic_Matern32.py @@ -101,7 +101,7 @@ class periodic_Matern32(kernpart): FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target) - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" if X2 is None: X2 = X FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) @@ -166,13 +166,13 @@ class periodic_Matern32(kernpart): dK_dper = mdot(dFX_dper,self.Gi,FX2.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX2.T) + mdot(FX,self.Gi,dFX2_dper.T) # np.add(target[:,:,0],dK_dvar, target[:,:,0]) - target[0] += np.sum(dK_dvar*partial) + target[0] += np.sum(dK_dvar*dL_dK) #np.add(target[:,:,1],dK_dlen, target[:,:,1]) - target[1] += np.sum(dK_dlen*partial) + target[1] += np.sum(dK_dlen*dL_dK) #np.add(target[:,:,2],dK_dper, target[:,:,2]) - target[2] += np.sum(dK_dper*partial) + target[2] += np.sum(dK_dper*dL_dK) - def dKdiag_dtheta(self,partial,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): """derivative of the diagonal covariance matrix with respect to the parameters""" FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) @@ -231,6 +231,6 @@ class periodic_Matern32(kernpart): dK_dper = 2* mdot(dFX_dper,self.Gi,FX.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX.T) - target[0] += np.sum(np.diag(dK_dvar)*partial) - target[1] += np.sum(np.diag(dK_dlen)*partial) - target[2] += np.sum(np.diag(dK_dper)*partial) + target[0] += np.sum(np.diag(dK_dvar)*dL_dKdiag) + target[1] += np.sum(np.diag(dK_dlen)*dL_dKdiag) + target[2] += np.sum(np.diag(dK_dper)*dL_dKdiag) diff --git a/GPy/kern/periodic_Matern52.py b/GPy/kern/periodic_Matern52.py index 8d1da8b1..c533961f 100644 --- a/GPy/kern/periodic_Matern52.py +++ b/GPy/kern/periodic_Matern52.py @@ -46,7 +46,7 @@ class periodic_Matern52(kernpart): r = np.sqrt(r1**2 + r2**2) psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2)) return r,omega[:,0:1], psi - + def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2): Gint1 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + 1./(omega1-omega2.T)*( np.sin((omega1-omega2.T)*self.upper+phi1-phi2.T) - np.sin((omega1-omega2.T)*self.lower+phi1-phi2.T) ) Gint2 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + np.cos(phi1-phi2.T)*(self.upper-self.lower) @@ -105,7 +105,7 @@ class periodic_Matern52(kernpart): FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target) - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" if X2 is None: X2 = X FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) @@ -178,13 +178,13 @@ class periodic_Matern52(kernpart): dK_dper = mdot(dFX_dper,self.Gi,FX2.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX2.T) + mdot(FX,self.Gi,dFX2_dper.T) # np.add(target[:,:,0],dK_dvar, target[:,:,0]) - target[0] += np.sum(dK_dvar*partial) + target[0] += np.sum(dK_dvar*dL_dK) #np.add(target[:,:,1],dK_dlen, target[:,:,1]) - target[1] += np.sum(dK_dlen*partial) + target[1] += np.sum(dK_dlen*dL_dK) #np.add(target[:,:,2],dK_dper, target[:,:,2]) - target[2] += np.sum(dK_dper*partial) + target[2] += np.sum(dK_dper*dL_dK) - def dKdiag_dtheta(self,partial,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): """derivative of the diagonal of the covariance matrix with respect to the parameters""" FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) @@ -251,6 +251,6 @@ class periodic_Matern52(kernpart): dG_dper = 1./self.variance*(3*self.lengthscale**5/(400*np.sqrt(5))*dGint_dper + 0.5*dlower_terms_dper) dK_dper = 2*mdot(dFX_dper,self.Gi,FX.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX.T) - target[0] += np.sum(np.diag(dK_dvar)*partial) - target[1] += np.sum(np.diag(dK_dlen)*partial) - target[2] += np.sum(np.diag(dK_dper)*partial) + target[0] += np.sum(np.diag(dK_dvar)*dL_dKdiag) + target[1] += np.sum(np.diag(dK_dlen)*dL_dKdiag) + target[2] += np.sum(np.diag(dK_dper)*dL_dKdiag) diff --git a/GPy/kern/periodic_exponential.py b/GPy/kern/periodic_exponential.py index 7f566f25..b966bbef 100644 --- a/GPy/kern/periodic_exponential.py +++ b/GPy/kern/periodic_exponential.py @@ -101,7 +101,7 @@ class periodic_exponential(kernpart): FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target) - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" if X2 is None: X2 = X FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) @@ -162,11 +162,11 @@ class periodic_exponential(kernpart): dK_dper = mdot(dFX_dper,self.Gi,FX2.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX2.T) + mdot(FX,self.Gi,dFX2_dper.T) - target[0] += np.sum(dK_dvar*partial) - target[1] += np.sum(dK_dlen*partial) - target[2] += np.sum(dK_dper*partial) + target[0] += np.sum(dK_dvar*dL_dK) + target[1] += np.sum(dK_dlen*dL_dK) + target[2] += np.sum(dK_dper*dL_dK) - def dKdiag_dtheta(self,partial,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): """derivative of the diagonal of the covariance matrix with respect to the parameters""" FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) @@ -222,7 +222,7 @@ class periodic_exponential(kernpart): dK_dper = 2*mdot(dFX_dper,self.Gi,FX.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX.T) - target[0] += np.sum(np.diag(dK_dvar)*partial) - target[1] += np.sum(np.diag(dK_dlen)*partial) - target[2] += np.sum(np.diag(dK_dper)*partial) - + target[0] += np.sum(np.diag(dK_dvar)*dL_dKdiag) + target[1] += np.sum(np.diag(dK_dlen)*dL_dKdiag) + target[2] += np.sum(np.diag(dK_dper)*dL_dKdiag) + diff --git a/GPy/kern/product.py b/GPy/kern/product.py index 92522418..3bad51c1 100644 --- a/GPy/kern/product.py +++ b/GPy/kern/product.py @@ -55,7 +55,7 @@ class product(kernpart): self.k2.Kdiag(X,target2) target += target1 * target2 - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters.""" if X2 is None: X2 = X K1 = np.zeros((X.shape[0],X2.shape[0])) @@ -65,13 +65,13 @@ class product(kernpart): k1_target = np.zeros(self.k1.Nparam) k2_target = np.zeros(self.k2.Nparam) - self.k1.dK_dtheta(partial*K2, X, X2, k1_target) - self.k2.dK_dtheta(partial*K1, X, X2, k2_target) + self.k1.dK_dtheta(dL_dK*K2, X, X2, k1_target) + self.k2.dK_dtheta(dL_dK*K1, X, X2, k2_target) target[:self.k1.Nparam] += k1_target target[self.k1.Nparam:] += k2_target - def dK_dX(self,partial,X,X2,target): + def dK_dX(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to X.""" if X2 is None: X2 = X K1 = np.zeros((X.shape[0],X2.shape[0])) @@ -79,19 +79,19 @@ class product(kernpart): self.k1.K(X,X2,K1) self.k2.K(X,X2,K2) - self.k1.dK_dX(partial*K2, X, X2, target) - self.k2.dK_dX(partial*K1, X, X2, target) + self.k1.dK_dX(dL_dK*K2, X, X2, target) + self.k2.dK_dX(dL_dK*K1, X, X2, target) - def dKdiag_dX(self,partial,X,target): + def dKdiag_dX(self,dL_dKdiag,X,target): target1 = np.zeros((X.shape[0],)) target2 = np.zeros((X.shape[0],)) self.k1.Kdiag(X,target1) self.k2.Kdiag(X,target2) - self.k1.dKdiag_dX(partial*target2, X, target) - self.k2.dKdiag_dX(partial*target1, X, target) + self.k1.dKdiag_dX(dL_dKdiag*target2, X, target) + self.k2.dKdiag_dX(dL_dKdiag*target1, X, target) - def dKdiag_dtheta(self,partial,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): """Compute the diagonal of the covariance matrix associated to X.""" target1 = np.zeros((X.shape[0],)) target2 = np.zeros((X.shape[0],)) @@ -100,8 +100,8 @@ class product(kernpart): k1_target = np.zeros(self.k1.Nparam) k2_target = np.zeros(self.k2.Nparam) - self.k1.dKdiag_dtheta(partial*target2, X, k1_target) - self.k2.dKdiag_dtheta(partial*target1, X, k2_target) + self.k1.dKdiag_dtheta(dL_dKdiag*target2, X, k1_target) + self.k2.dKdiag_dtheta(dL_dKdiag*target1, X, k2_target) target[:self.k1.Nparam] += k1_target target[self.k1.Nparam:] += k2_target diff --git a/GPy/kern/product_orthogonal.py b/GPy/kern/product_orthogonal.py index a231cf8b..b0112199 100644 --- a/GPy/kern/product_orthogonal.py +++ b/GPy/kern/product_orthogonal.py @@ -46,7 +46,7 @@ class product_orthogonal(kernpart): self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],target2) target += target1 * target2 - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters.""" if X2 is None: X2 = X K1 = np.zeros((X.shape[0],X2.shape[0])) @@ -54,8 +54,8 @@ class product_orthogonal(kernpart): self.k1.K(X[:,:self.k1.D],X2[:,:self.k1.D],K1) self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],K2) - self.k1.dK_dtheta(partial*K2, X[:,:self.k1.D], X2[:,:self.k1.D], target[:self.k1.Nparam]) - self.k2.dK_dtheta(partial*K1, X[:,self.k1.D:], X2[:,self.k1.D:], target[self.k1.Nparam:]) + self.k1.dK_dtheta(dL_dK*K2, X[:,:self.k1.D], X2[:,:self.k1.D], target[:self.k1.Nparam]) + self.k2.dK_dtheta(dL_dK*K1, X[:,self.k1.D:], X2[:,self.k1.D:], target[self.k1.Nparam:]) def Kdiag(self,X,target): """Compute the diagonal of the covariance matrix associated to X.""" @@ -65,15 +65,15 @@ class product_orthogonal(kernpart): self.k2.Kdiag(X[:,self.k1.D:],target2) target += target1 * target2 - def dKdiag_dtheta(self,partial,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): K1 = np.zeros(X.shape[0]) K2 = np.zeros(X.shape[0]) self.k1.Kdiag(X[:,:self.k1.D],K1) self.k2.Kdiag(X[:,self.k1.D:],K2) - self.k1.dKdiag_dtheta(partial*K2,X[:,:self.k1.D],target[:self.k1.Nparam]) - self.k2.dKdiag_dtheta(partial*K1,X[:,self.k1.D:],target[self.k1.Nparam:]) + self.k1.dKdiag_dtheta(dL_dKdiag*K2,X[:,:self.k1.D],target[:self.k1.Nparam]) + self.k2.dKdiag_dtheta(dL_dKdiag*K1,X[:,self.k1.D:],target[self.k1.Nparam:]) - def dK_dX(self,partial,X,X2,target): + def dK_dX(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to X.""" if X2 is None: X2 = X K1 = np.zeros((X.shape[0],X2.shape[0])) @@ -81,15 +81,15 @@ class product_orthogonal(kernpart): self.k1.K(X[:,0:self.k1.D],X2[:,0:self.k1.D],K1) self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],K2) - self.k1.dK_dX(partial*K2, X[:,:self.k1.D], X2[:,:self.k1.D], target) - self.k2.dK_dX(partial*K1, X[:,self.k1.D:], X2[:,self.k1.D:], target) + self.k1.dK_dX(dL_dK*K2, X[:,:self.k1.D], X2[:,:self.k1.D], target) + self.k2.dK_dX(dL_dK*K1, X[:,self.k1.D:], X2[:,self.k1.D:], target) - def dKdiag_dX(self, partial, X, target): + def dKdiag_dX(self, dL_dKdiag, X, target): K1 = np.zeros(X.shape[0]) K2 = np.zeros(X.shape[0]) self.k1.Kdiag(X[:,0:self.k1.D],K1) self.k2.Kdiag(X[:,self.k1.D:],K2) - self.k1.dK_dX(partial*K2, X[:,:self.k1.D], target) - self.k2.dK_dX(partial*K1, X[:,self.k1.D:], target) + self.k1.dK_dX(dL_dKdiag*K2, X[:,:self.k1.D], target) + self.k2.dK_dX(dL_dKdiag*K1, X[:,self.k1.D:], target) diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py index 16eda459..3c3d59e6 100644 --- a/GPy/kern/rbf.py +++ b/GPy/kern/rbf.py @@ -82,27 +82,27 @@ class rbf(kernpart): def Kdiag(self,X,target): np.add(target,self.variance,target) - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): self._K_computations(X,X2) - target[0] += np.sum(self._K_dvar*partial) + target[0] += np.sum(self._K_dvar*dL_dK) if self.ARD == True: dl = self._K_dvar[:,:,None]*self.variance*self._K_dist2/self.lengthscale - target[1:] += (dl*partial[:,:,None]).sum(0).sum(0) + target[1:] += (dl*dL_dK[:,:,None]).sum(0).sum(0) else: - target[1] += np.sum(self._K_dvar*self.variance*(self._K_dist2.sum(-1))/self.lengthscale*partial) - #np.sum(self._K_dvar*self.variance*self._K_dist2/self.lengthscale*partial) + target[1] += np.sum(self._K_dvar*self.variance*(self._K_dist2.sum(-1))/self.lengthscale*dL_dK) + #np.sum(self._K_dvar*self.variance*self._K_dist2/self.lengthscale*dL_dK) - def dKdiag_dtheta(self,partial,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): #NB: derivative of diagonal elements wrt lengthscale is 0 - target[0] += np.sum(partial) + target[0] += np.sum(dL_dKdiag) - def dK_dX(self,partial,X,X2,target): + def dK_dX(self,dL_dK,X,X2,target): self._K_computations(X,X2) _K_dist = X[:,None,:]-X2[None,:,:] dK_dX = np.transpose(-self.variance*self._K_dvar[:,:,np.newaxis]*_K_dist/self.lengthscale2,(1,0,2)) - target += np.sum(dK_dX*partial.T[:,:,None],0) + target += np.sum(dK_dX*dL_dK.T[:,:,None],0) - def dKdiag_dX(self,partial,X,target): + def dKdiag_dX(self,dL_dKdiag,X,target): pass @@ -113,69 +113,69 @@ class rbf(kernpart): def psi0(self,Z,mu,S,target): target += self.variance - def dpsi0_dtheta(self,partial,Z,mu,S,target): - target[0] += np.sum(partial) + def dpsi0_dtheta(self,dL_dpsi0,Z,mu,S,target): + target[0] += np.sum(dL_dpsi0) - def dpsi0_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi0_dmuS(self,dL_dpsi0,Z,mu,S,target_mu,target_S): pass def psi1(self,Z,mu,S,target): self._psi_computations(Z,mu,S) target += self._psi1 - def dpsi1_dtheta(self,partial,Z,mu,S,target): + def dpsi1_dtheta(self,dL_dpsi1,Z,mu,S,target): self._psi_computations(Z,mu,S) denom_deriv = S[:,None,:]/(self.lengthscale**3+self.lengthscale*S[:,None,:]) d_length = self._psi1[:,:,None]*(self.lengthscale*np.square(self._psi1_dist/(self.lengthscale2+S[:,None,:])) + denom_deriv) - target[0] += np.sum(partial*self._psi1/self.variance) - dpsi1_dlength = d_length*partial[:,:,None] + target[0] += np.sum(dL_dpsi1*self._psi1/self.variance) + dpsi1_dlength = d_length*dL_dpsi1[:,:,None] if not self.ARD: target[1] += dpsi1_dlength.sum() else: target[1:] += dpsi1_dlength.sum(0).sum(0) - def dpsi1_dZ(self,partial,Z,mu,S,target): + def dpsi1_dZ(self,dL_dpsi1,Z,mu,S,target): self._psi_computations(Z,mu,S) denominator = (self.lengthscale2*(self._psi1_denom)) dpsi1_dZ = - self._psi1[:,:,None] * ((self._psi1_dist/denominator)) - target += np.sum(partial.T[:,:,None] * dpsi1_dZ, 0) + target += np.sum(dL_dpsi1.T[:,:,None] * dpsi1_dZ, 0) - def dpsi1_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi1_dmuS(self,dL_dpsi1,Z,mu,S,target_mu,target_S): self._psi_computations(Z,mu,S) tmp = self._psi1[:,:,None]/self.lengthscale2/self._psi1_denom - target_mu += np.sum(partial.T[:, :, None]*tmp*self._psi1_dist,1) - target_S += np.sum(partial.T[:, :, None]*0.5*tmp*(self._psi1_dist_sq-1),1) + target_mu += np.sum(dL_dpsi1.T[:, :, None]*tmp*self._psi1_dist,1) + target_S += np.sum(dL_dpsi1.T[:, :, None]*0.5*tmp*(self._psi1_dist_sq-1),1) def psi2(self,Z,mu,S,target): self._psi_computations(Z,mu,S) target += self._psi2 - def dpsi2_dtheta(self,partial,Z,mu,S,target): + def dpsi2_dtheta(self,dL_dpsi2,Z,mu,S,target): """Shape N,M,M,Ntheta""" self._psi_computations(Z,mu,S) d_var = 2.*self._psi2/self.variance d_length = self._psi2[:,:,:,None]*(0.5*self._psi2_Zdist_sq*self._psi2_denom + 2.*self._psi2_mudist_sq + 2.*S[:,None,None,:]/self.lengthscale2)/(self.lengthscale*self._psi2_denom) - target[0] += np.sum(partial*d_var) - dpsi2_dlength = d_length*partial[:,:,:,None] + target[0] += np.sum(dL_dpsi2*d_var) + dpsi2_dlength = d_length*dL_dpsi2[:,:,:,None] if not self.ARD: target[1] += dpsi2_dlength.sum() else: target[1:] += dpsi2_dlength.sum(0).sum(0).sum(0) - - def dpsi2_dZ(self,partial,Z,mu,S,target): + + def dpsi2_dZ(self,dL_dpsi2,Z,mu,S,target): self._psi_computations(Z,mu,S) term1 = 0.5*self._psi2_Zdist/self.lengthscale2 # M, M, Q term2 = self._psi2_mudist/self._psi2_denom/self.lengthscale2 # N, M, M, Q dZ = self._psi2[:,:,:,None] * (term1[None] + term2) - target += (partial[:,:,:,None]*dZ).sum(0).sum(0) + target += (dL_dpsi2[:,:,:,None]*dZ).sum(0).sum(0) - def dpsi2_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,target_mu,target_S): """Think N,M,M,Q """ self._psi_computations(Z,mu,S) tmp = self._psi2[:,:,:,None]/self.lengthscale2/self._psi2_denom - target_mu += (partial[:,:,:,None]*-tmp*2.*self._psi2_mudist).sum(1).sum(1) - target_S += (partial[:,:,:,None]*tmp*(2.*self._psi2_mudist_sq-1)).sum(1).sum(1) + target_mu += (dL_dpsi2[:,:,:,None]*-tmp*2.*self._psi2_mudist).sum(1).sum(1) + target_S += (dL_dpsi2[:,:,:,None]*tmp*(2.*self._psi2_mudist_sq-1)).sum(1).sum(1) #---------------------------------------# diff --git a/GPy/kern/symmetric.py b/GPy/kern/symmetric.py index d493bfb1..c3b046c7 100644 --- a/GPy/kern/symmetric.py +++ b/GPy/kern/symmetric.py @@ -51,7 +51,7 @@ class symmetric(kernpart): self.k.K(X,AX2,target) self.k.K(AX,AX2,target) - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters.""" AX = np.dot(X,self.transform) if X2 is None: @@ -59,13 +59,13 @@ class symmetric(kernpart): ZX2 = AX else: AX2 = np.dot(X2, self.transform) - self.k.dK_dtheta(partial,X,X2,target) - self.k.dK_dtheta(partial,AX,X2,target) - self.k.dK_dtheta(partial,X,AX2,target) - self.k.dK_dtheta(partial,AX,AX2,target) + self.k.dK_dtheta(dL_dK,X,X2,target) + self.k.dK_dtheta(dL_dK,AX,X2,target) + self.k.dK_dtheta(dL_dK,X,AX2,target) + self.k.dK_dtheta(dL_dK,AX,AX2,target) - def dK_dX(self,partial,X,X2,target): + def dK_dX(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to X.""" AX = np.dot(X,self.transform) if X2 is None: @@ -73,10 +73,10 @@ class symmetric(kernpart): ZX2 = AX else: AX2 = np.dot(X2, self.transform) - self.k.dK_dX(partial, X, X2, target) - self.k.dK_dX(partial, AX, X2, target) - self.k.dK_dX(partial, X, AX2, target) - self.k.dK_dX(partial, AX ,AX2, target) + self.k.dK_dX(dL_dK, X, X2, target) + self.k.dK_dX(dL_dK, AX, X2, target) + self.k.dK_dX(dL_dK, X, AX2, target) + self.k.dK_dX(dL_dK, AX ,AX2, target) def Kdiag(self,X,target): """Compute the diagonal of the covariance matrix associated to X.""" @@ -84,9 +84,9 @@ class symmetric(kernpart): self.K(X,X,foo) target += np.diag(foo) - def dKdiag_dX(self,partial,X,target): + def dKdiag_dX(self,dL_dKdiag,X,target): raise NotImplementedError - def dKdiag_dtheta(self,partial,X,target): + def dKdiag_dtheta(self,dL_dKdiag,X,target): """Compute the diagonal of the covariance matrix associated to X.""" raise NotImplementedError diff --git a/GPy/kern/white.py b/GPy/kern/white.py index b3b00c48..f5d6894a 100644 --- a/GPy/kern/white.py +++ b/GPy/kern/white.py @@ -37,50 +37,50 @@ class white(kernpart): def Kdiag(self,X,target): target += self.variance - def dK_dtheta(self,partial,X,X2,target): + def dK_dtheta(self,dL_dK,X,X2,target): if X.shape==X2.shape: if np.all(X==X2): - target += np.trace(partial) + target += np.trace(dL_dK) - def dKdiag_dtheta(self,partial,X,target): - target += np.sum(partial) + def dKdiag_dtheta(self,dL_dKdiag,X,target): + target += np.sum(dL_dKdiag) - def dK_dX(self,partial,X,X2,target): + def dK_dX(self,dL_dK,X,X2,target): pass - def dKdiag_dX(self,partial,X,target): + def dKdiag_dX(self,dL_dKdiag,X,target): pass def psi0(self,Z,mu,S,target): target += self.variance - def dpsi0_dtheta(self,partial,Z,mu,S,target): - target += partial.sum() + def dpsi0_dtheta(self,dL_dpsi0,Z,mu,S,target): + target += dL_dpsi0.sum() - def dpsi0_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi0_dmuS(self,dL_dpsi0,Z,mu,S,target_mu,target_S): pass def psi1(self,Z,mu,S,target): pass - def dpsi1_dtheta(self,partial,Z,mu,S,target): + def dpsi1_dtheta(self,dL_dpsi1,Z,mu,S,target): pass - def dpsi1_dZ(self,partial,Z,mu,S,target): + def dpsi1_dZ(self,dL_dpsi1,Z,mu,S,target): pass - def dpsi1_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi1_dmuS(self,dL_dpsi1,Z,mu,S,target_mu,target_S): pass def psi2(self,Z,mu,S,target): pass - def dpsi2_dZ(self,partial,Z,mu,S,target): + def dpsi2_dZ(self,dL_dpsi2,Z,mu,S,target): pass - def dpsi2_dtheta(self,partial,Z,mu,S,target): + def dpsi2_dtheta(self,dL_dpsi2,Z,mu,S,target): pass - def dpsi2_dmuS(self,partial,Z,mu,S,target_mu,target_S): + def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,target_mu,target_S): pass diff --git a/GPy/models/GP.py b/GPy/models/GP.py index 08ac1bb1..1d985c33 100644 --- a/GPy/models/GP.py +++ b/GPy/models/GP.py @@ -129,7 +129,7 @@ class GP(model): For the likelihood parameters, pass in alpha = K^-1 y """ - return np.hstack((self.kern.dK_dtheta(partial=self.dL_dK,X=self.X,slices1=self.Xslices,slices2=self.Xslices), self.likelihood._gradients(partial=np.diag(self.dL_dK)))) + return np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK,X=self.X,slices1=self.Xslices,slices2=self.Xslices), self.likelihood._gradients(partial=np.diag(self.dL_dK)))) def _raw_predict(self,_Xnew,slices=None, full_cov=False): """ From f562d6cd46779f2bb1adebb50eed7ef4b0de0c69 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Mon, 11 Mar 2013 12:21:12 +0000 Subject: [PATCH 07/18] deprecated flapack, namespace changed to lapack.flapack --- GPy/kern/Matern52.py | 2 +- GPy/likelihoods/EP.py | 10 +++++----- GPy/util/linalg.py | 6 +++--- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/GPy/kern/Matern52.py b/GPy/kern/Matern52.py index 377526d5..9338db15 100644 --- a/GPy/kern/Matern52.py +++ b/GPy/kern/Matern52.py @@ -90,7 +90,7 @@ class Matern52(kernpart): else: dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist)) * dist2M.sum(-1)*invdist #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist)) * dist2M.sum(-1)*invdist - target[1] += np.sum(dl*dL_dKdiag) + target[1] += np.sum(dl*dL_dK) def dKdiag_dtheta(self,dL_dKdiag,X,target): """derivative of the diagonal of the covariance matrix with respect to the parameters.""" diff --git a/GPy/likelihoods/EP.py b/GPy/likelihoods/EP.py index efd887ae..30b21d9b 100644 --- a/GPy/likelihoods/EP.py +++ b/GPy/likelihoods/EP.py @@ -114,7 +114,7 @@ class EP(likelihood): Sroot_tilde_K = np.sqrt(self.tau_tilde)[:,None]*K B = np.eye(self.N) + np.sqrt(self.tau_tilde)[None,:]*Sroot_tilde_K L = jitchol(B) - V,info = linalg.flapack.dtrtrs(L,Sroot_tilde_K,lower=1) + V,info = linalg.lapack.flapack.dtrtrs(L,Sroot_tilde_K,lower=1) Sigma = K - np.dot(V.T,V) mu = np.dot(Sigma,self.v_tilde) epsilon_np1 = sum((self.tau_tilde-self.np1[-1])**2)/self.N @@ -190,7 +190,7 @@ class EP(likelihood): #Posterior distribution parameters update LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau L = jitchol(LLT) - V,info = linalg.flapack.dtrtrs(L,Kmn,lower=1) + V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1) Sigma_diag = np.sum(V*V,-2) si = np.sum(V.T*V[:,i],-1) mu = mu + (Delta_v-Delta_tau*mu[i])*si @@ -198,8 +198,8 @@ class EP(likelihood): #Sigma recomputation with Cholesky decompositon LLT0 = LLT0 + np.dot(Kmn*tau_tilde[None,:],Kmn.T) L = jitchol(LLT) - V,info = linalg.flapack.dtrtrs(L,Kmn,lower=1) - V2,info = linalg.flapack.dtrtrs(L.T,V,lower=0) + V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1) + V2,info = linalg.lapack.flapack.dtrtrs(L.T,V,lower=0) Sigma_diag = np.sum(V*V,-2) Knmv_tilde = np.dot(Kmn,v_tilde) mu = np.dot(V2.T,Knmv_tilde) @@ -297,7 +297,7 @@ class EP(likelihood): P = (Diag / Diag0)[:,None] * P0 RPT0 = np.dot(R0,P0.T) L = jitchol(np.eye(self.M) + np.dot(RPT0,(1./Diag0 - Diag/(Diag0**2))[:,None]*RPT0.T)) - R,info = linalg.flapack.dtrtrs(L,R0,lower=1) + R,info = linalg.lapack.flapack.dtrtrs(L,R0,lower=1) RPT = np.dot(R,P.T) Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1) self.w = Diag * self.v_tilde diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index 7414eb29..26105789 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -11,7 +11,7 @@ import re import pdb import cPickle import types -import scipy.lib.lapack.flapack +#import scipy.lib.lapack.flapack import scipy as sp def mdot(*args): @@ -101,7 +101,7 @@ def chol_inv(L): """ - return linalg.flapack.dtrtri(L, lower = True)[0] + return linalg.lapack.flapack.dtrtri(L, lower = True)[0] def multiple_pdinv(A): @@ -118,7 +118,7 @@ def multiple_pdinv(A): N = A.shape[-1] chols = [jitchol(A[:,:,i]) for i in range(N)] halflogdets = [np.sum(np.log(np.diag(L[0]))) for L in chols] - invs = [linalg.flapack.dpotri(L[0],True)[0] for L in chols] + invs = [linalg.lapack.flapack.dpotri(L[0],True)[0] for L in chols] invs = [np.triu(I)+np.triu(I,1).T for I in invs] return np.dstack(invs),np.array(halflogdets) From e32afa11e5b437ef5db6bfd015f4f87936723bd0 Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Mon, 11 Mar 2013 12:33:03 +0000 Subject: [PATCH 08/18] added GPy.tests(), removed some useless tests --- GPy/__init__.py | 5 +++++ GPy/testing/bgplvm_tests.py | 5 ++++- GPy/testing/kernel_tests.py | 1 - GPy/testing/unit_tests.py | 11 ----------- setup.py | 6 ++---- 5 files changed, 11 insertions(+), 17 deletions(-) diff --git a/GPy/__init__.py b/GPy/__init__.py index c0772c27..6c43e471 100644 --- a/GPy/__init__.py +++ b/GPy/__init__.py @@ -9,3 +9,8 @@ import util import examples from core import priors import likelihoods +import testing +from numpy.testing import Tester + +def tests(): + Tester(testing).test(verbose=10) diff --git a/GPy/testing/bgplvm_tests.py b/GPy/testing/bgplvm_tests.py index c49bdfda..e3bd2b36 100644 --- a/GPy/testing/bgplvm_tests.py +++ b/GPy/testing/bgplvm_tests.py @@ -12,6 +12,7 @@ class BGPLVMTests(unittest.TestCase): k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) K = k.K(X) Y = np.random.multivariate_normal(np.zeros(N),K,D).T + Y -= Y.mean(axis=0) k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M) m.constrain_positive('(rbf|bias|noise|white|S)') @@ -24,6 +25,7 @@ class BGPLVMTests(unittest.TestCase): k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) K = k.K(X) Y = np.random.multivariate_normal(np.zeros(N),K,D).T + Y -= Y.mean(axis=0) k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001) m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M) m.constrain_positive('(linear|bias|noise|white|S)') @@ -36,13 +38,14 @@ class BGPLVMTests(unittest.TestCase): k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) K = k.K(X) Y = np.random.multivariate_normal(np.zeros(N),K,D).T + Y -= Y.mean(axis=0) k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M) m.constrain_positive('(rbf|bias|noise|white|S)') m.randomize() self.assertTrue(m.checkgrad()) - + if __name__ == "__main__": print "Running unit tests, please be (very) patient..." unittest.main() diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 3d738106..bb809ea6 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -13,7 +13,6 @@ class KernelTests(unittest.TestCase): X = np.random.rand(5,5) Y = np.ones((5,1)) m = GPy.models.GP_regression(X,Y,K) - print m self.assertTrue(m.checkgrad()) def test_coregionalisation(self): diff --git a/GPy/testing/unit_tests.py b/GPy/testing/unit_tests.py index 61fb15bb..55963805 100644 --- a/GPy/testing/unit_tests.py +++ b/GPy/testing/unit_tests.py @@ -177,17 +177,6 @@ class GradientTests(unittest.TestCase): m.approximate_likelihood() self.assertTrue(m.checkgrad()) - def test_warped_GP(self): - xmin, xmax = 1, 2.5*np.pi - b, C, SNR = 1, 0, 0.1 - X = np.linspace(xmin, xmax, 500) - y = b*X + C + 1*np.sin(X) - y += 0.05*np.random.randn(len(X)) - X, y = X[:, None], y[:, None] - m = GPy.models.warpedGP(X, y, warping_terms = 3) - m.constrain_positive('(tanh_a|tanh_b|rbf|white|bias)') - self.assertTrue(m.checkgrad()) - if __name__ == "__main__": print "Running unit tests, please be (very) patient..." diff --git a/setup.py b/setup.py index d24171e2..b701b74d 100644 --- a/setup.py +++ b/setup.py @@ -3,8 +3,6 @@ import os from setuptools import setup -#from numpy.distutils.core import Extension, setup -#from sphinx.setup_command import BuildDoc # Version number version = '0.1.3' @@ -14,12 +12,12 @@ def read(fname): setup(name = 'GPy', version = version, - author = 'James Hensman, Nicolo Fusi, Ricardo Andrade, Nicolas Durrande, Alan Saul, Neil D. Lawrence', + author = read('AUTHORS.txt'), author_email = "james.hensman@gmail.com", description = ("The Gaussian Process Toolbox"), license = "BSD 3-clause", keywords = "machine-learning gaussian-processes kernels", - url = "http://ml.sheffield.ac.uk/GPy/", + url = "http://sheffieldml.github.com/GPy/", packages = ['GPy', 'GPy.core', 'GPy.kern', 'GPy.util', 'GPy.models', 'GPy.inference', 'GPy.examples', 'GPy.likelihoods'], package_dir={'GPy': 'GPy'}, package_data = {'GPy': ['GPy/examples']}, From a86676016247b222664c177bfc4e9dc834500c00 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Mon, 11 Mar 2013 12:39:44 +0000 Subject: [PATCH 09/18] Removed unused partial1 --- GPy/kern/kern.py | 2 +- GPy/models/sparse_GP.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py index f1a5bd45..87e67f33 100644 --- a/GPy/kern/kern.py +++ b/GPy/kern/kern.py @@ -399,7 +399,7 @@ class kern(parameterised): return target - def dpsi2_dtheta(self,dL_dpsi2,partial1,Z,mu,S,slices1=None,slices2=None): + def dpsi2_dtheta(self,dL_dpsi2,Z,mu,S,slices1=None,slices2=None): """Returns shape (N,M,M,Ntheta)""" slices1, slices2 = self._process_slices(slices1,slices2) target = np.zeros(self.Nparam) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index ff00faea..e2019d99 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -208,7 +208,7 @@ class sparse_GP(GP): if self.has_uncertain_inputs: dL_dtheta += self.kern.dpsi0_dtheta(self.dL_dpsi0, self.Z,self.X,self.X_uncertainty) dL_dtheta += self.kern.dpsi1_dtheta(self.dL_dpsi1.T,self.Z,self.X, self.X_uncertainty) - dL_dtheta += self.kern.dpsi2_dtheta(self.dL_dpsi2,self.dL_dpsi1.T, self.Z,self.X, self.X_uncertainty) + dL_dtheta += self.kern.dpsi2_dtheta(self.dL_dpsi2, self.Z,self.X, self.X_uncertainty) else: dL_dtheta += self.kern.dK_dtheta(self.dL_dpsi1,self.Z,self.X) dL_dtheta += self.kern.dKdiag_dtheta(self.dL_dpsi0, self.X) From f98e52ffe8cd5678b178eb69fa771b23af55125f Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Mon, 11 Mar 2013 12:40:14 +0000 Subject: [PATCH 10/18] now running nosetest doesn't run unittests twice --- GPy/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/GPy/__init__.py b/GPy/__init__.py index 6c43e471..fa69dac3 100644 --- a/GPy/__init__.py +++ b/GPy/__init__.py @@ -11,6 +11,8 @@ from core import priors import likelihoods import testing from numpy.testing import Tester +from nose.tools import nottest +@nottest def tests(): Tester(testing).test(verbose=10) From e511bb69cf25cfa3c6c8aea946a85fcdf3f437d2 Mon Sep 17 00:00:00 2001 From: Nicolas Date: Mon, 11 Mar 2013 12:40:29 +0000 Subject: [PATCH 11/18] added in documentation the current status of kernel implementation --- doc/Figures/tick.png | Bin 0 -> 175 bytes doc/GPy.examples.rst | 16 +++++++++++++ doc/GPy.kern.rst | 28 ++++++++++++++++++----- doc/kernel_implementation.rst | 41 +++++++++++++++++++++++++--------- doc/tuto_kernel_overview.rst | 2 +- 5 files changed, 70 insertions(+), 17 deletions(-) create mode 100644 doc/Figures/tick.png diff --git a/doc/Figures/tick.png b/doc/Figures/tick.png new file mode 100644 index 0000000000000000000000000000000000000000..1175c8021717199329a79061bbf1e02c49f677fc GIT binary patch literal 175 zcmeAS@N?(olHy`uVBq!ia0vp@K+MO&3?$hCyB`86tpJ}8*Z=?j1DW&Y%?k(!NJ>g_ zaBwIoDcQDd8xIeUl9Cb&3rkm5*O@bCo;`ae#qL%PRO;aA;uyklJvo7aQH|B$Xadv1 z2w{a1iIyb8h6|zsT?;v#4E-4`_ . The following figure gives a summary of most of them: +Many kernels are already implemented in GPy. A comprehensive list can be found `here `_ and the following figure gives a summary of most of them: .. figure:: Figures/tuto_kern_overview_allkern.png :align: center From da219a66e811767303f858d59d37fbde96c7eb76 Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Mon, 11 Mar 2013 12:40:53 +0000 Subject: [PATCH 12/18] added init --- GPy/testing/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 GPy/testing/__init__.py diff --git a/GPy/testing/__init__.py b/GPy/testing/__init__.py new file mode 100644 index 00000000..e69de29b From 25d73b13e99ecf6471063573d3ab4f2d02b6f587 Mon Sep 17 00:00:00 2001 From: Nicolas Date: Mon, 11 Mar 2013 12:59:39 +0000 Subject: [PATCH 13/18] update in the documentation on kernel implementation --- doc/kernel_implementation.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/doc/kernel_implementation.rst b/doc/kernel_implementation.rst index 888d1ee5..99ee006b 100644 --- a/doc/kernel_implementation.rst +++ b/doc/kernel_implementation.rst @@ -35,4 +35,13 @@ spline |tick| |tick| |tick| |tick| |tick| white |tick| |tick| |tick| |tick| |tick| |tick| |tick| |tick| |tick| |tick| ==================== =========== ====== ======= =========== =============== ======= =========== ====== ====== ======= +Depending on the use, all functions may not be required + + * ``get/set, K, Kdiag``: compulsory + * ``dK_dtheta``: necessary to optimize the model + * ``dKdiag_dtheta``: sparse models, BGPLVM, GPs with uncertain inputs + * ``dK_dX``: sparse models, GPLVM, BGPLVM, GPs with uncertain inputs + * ``dKdiag_dX``: sparse models, BGPLVM, GPs with uncertain inputs + * ``psi0, psi1, psi2``: BGPLVM, GPs with uncertain inputs + .. |tick| image:: Figures/tick.png From 16a23758c6f64b568c2ce985ef98454932ed3350 Mon Sep 17 00:00:00 2001 From: Nicolas Date: Mon, 11 Mar 2013 13:13:18 +0000 Subject: [PATCH 14/18] example files for tutorials are now in Neil's format --- GPy/examples/__init__.py | 1 + GPy/examples/tuto_GP_regression.py | 56 -------- GPy/examples/tuto_kernel_overview.py | 139 ------------------ GPy/examples/tutorials.py | 201 +++++++++++++++++++++++++++ doc/tuto_GP_regression.rst | 2 +- doc/tuto_kernel_overview.rst | 2 +- 6 files changed, 204 insertions(+), 197 deletions(-) delete mode 100644 GPy/examples/tuto_GP_regression.py delete mode 100644 GPy/examples/tuto_kernel_overview.py create mode 100644 GPy/examples/tutorials.py diff --git a/GPy/examples/__init__.py b/GPy/examples/__init__.py index 2f3cf0f4..ce4618ac 100644 --- a/GPy/examples/__init__.py +++ b/GPy/examples/__init__.py @@ -6,3 +6,4 @@ import classification import regression import unsupervised +import tutorials diff --git a/GPy/examples/tuto_GP_regression.py b/GPy/examples/tuto_GP_regression.py deleted file mode 100644 index b3953de0..00000000 --- a/GPy/examples/tuto_GP_regression.py +++ /dev/null @@ -1,56 +0,0 @@ -# The detailed explanations of the commands used in this file can be found in the tutorial section - -import pylab as pb -pb.ion() -import numpy as np -import GPy - -X = np.random.uniform(-3.,3.,(20,1)) -Y = np.sin(X) + np.random.randn(20,1)*0.05 - -kernel = GPy.kern.rbf(D=1, variance=1., lengthscale=1.) - -m = GPy.models.GP_regression(X,Y,kernel) - -print m -m.plot() - -m.constrain_positive('') - -m.unconstrain('') # Required to remove the previous constrains -m.constrain_positive('rbf_variance') -m.constrain_bounded('lengthscale',1.,10. ) -m.constrain_fixed('noise',0.0025) - -m.optimize() - -m.optimize_restarts(Nrestarts = 10) - -########################### -# 2-dimensional example # -########################### - -import pylab as pb -pb.ion() -import numpy as np -import GPy - -# sample inputs and outputs -X = np.random.uniform(-3.,3.,(50,2)) -Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(50,1)*0.05 - -# define kernel -ker = GPy.kern.Matern52(2,ARD=True) + GPy.kern.white(2) - -# create simple GP model -m = GPy.models.GP_regression(X,Y,ker) - -# contrain all parameters to be positive -m.constrain_positive('') - -# optimize and plot -pb.figure() -m.optimize('tnc', max_f_eval = 1000) - -m.plot() -print(m) diff --git a/GPy/examples/tuto_kernel_overview.py b/GPy/examples/tuto_kernel_overview.py deleted file mode 100644 index ebd19d76..00000000 --- a/GPy/examples/tuto_kernel_overview.py +++ /dev/null @@ -1,139 +0,0 @@ -# The detailed explanations of the commands used in this file can be found in the tutorial section - -import pylab as pb -import numpy as np -import GPy -pb.ion() - -ker1 = GPy.kern.rbf(1) # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.) -ker2 = GPy.kern.rbf(D=1, variance = .75, lengthscale=2.) -ker3 = GPy.kern.rbf(1, .5, .5) - -print ker2 -ker1.plot() -ker2.plot() -ker3.plot() - -k1 = GPy.kern.rbf(1,1.,2.) -k2 = GPy.kern.Matern32(1, 0.5, 0.2) - -# Product of kernels -k_prod = k1.prod(k2) -k_prodorth = k1.prod_orthogonal(k2) - -# Sum of kernels -k_add = k1.add(k2) -k_addorth = k1.add_orthogonal(k2) - -pb.figure(figsize=(8,8)) -pb.subplot(2,2,1) -k_prod.plot() -pb.title('prod') -pb.subplot(2,2,2) -k_prodorth.plot() -pb.title('prod_orthogonal') -pb.subplot(2,2,3) -k_add.plot() -pb.title('add') -pb.subplot(2,2,4) -k_addorth.plot() -pb.title('add_orthogonal') -pb.subplots_adjust(wspace=0.3, hspace=0.3) - -k1 = GPy.kern.rbf(1,1.,2) -k2 = GPy.kern.periodic_Matern52(1,variance=1e3, lengthscale=1, period = 1.5, lower=-5., upper = 5) - -k = k1 * k2 # equivalent to k = k1.prod(k2) -print k - -# Simulate sample paths -X = np.linspace(-5,5,501)[:,None] -Y = np.random.multivariate_normal(np.zeros(501),k.K(X),1) - -# plot -pb.figure(figsize=(10,4)) -pb.subplot(1,2,1) -k.plot() -pb.subplot(1,2,2) -pb.plot(X,Y.T) -pb.ylabel("Sample path") -pb.subplots_adjust(wspace=0.3) - -k = (k1+k2)*(k1+k2) -print k.parts[0].name, '\n', k.parts[1].name, '\n', k.parts[2].name, '\n', k.parts[3].name - -k1 = GPy.kern.rbf(1) -k2 = GPy.kern.Matern32(1) -k3 = GPy.kern.white(1) - -k = k1 + k2 + k3 -print k - -k.constrain_positive('var') -k.constrain_fixed(np.array([1]),1.75) -k.tie_param('len') -k.unconstrain('white') -k.constrain_bounded('white',lower=1e-5,upper=.5) -print k - -k_cst = GPy.kern.bias(1,variance=1.) -k_mat = GPy.kern.Matern52(1,variance=1., lengthscale=3) -Kanova = (k_cst + k_mat).prod_orthogonal(k_cst + k_mat) -print Kanova - -# sample inputs and outputs -X = np.random.uniform(-3.,3.,(40,2)) -Y = 0.5*X[:,:1] + 0.5*X[:,1:] + 2*np.sin(X[:,:1]) * np.sin(X[:,1:]) - -# Create GP regression model -m = GPy.models.GP_regression(X,Y,Kanova) -pb.figure(figsize=(5,5)) -m.plot() - -pb.figure(figsize=(20,3)) -pb.subplots_adjust(wspace=0.5) -pb.subplot(1,5,1) -m.plot() -pb.subplot(1,5,2) -pb.ylabel("= ",rotation='horizontal',fontsize='30') -pb.subplot(1,5,3) -m.plot(which_functions=[False,True,False,False]) -pb.ylabel("cst +",rotation='horizontal',fontsize='30') -pb.subplot(1,5,4) -m.plot(which_functions=[False,False,True,False]) -pb.ylabel("+ ",rotation='horizontal',fontsize='30') -pb.subplot(1,5,5) -pb.ylabel("+ ",rotation='horizontal',fontsize='30') -m.plot(which_functions=[False,False,False,True]) - -import pylab as pb -import numpy as np -import GPy -pb.ion() - -ker1 = GPy.kern.rbf(D=1) # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.) -ker2 = GPy.kern.rbf(D=1, variance = .75, lengthscale=3.) -ker3 = GPy.kern.rbf(1, .5, .25) - -ker1.plot() -ker2.plot() -ker3.plot() -#pb.savefig("Figures/tuto_kern_overview_basicdef.png") - -kernels = [GPy.kern.rbf(1), GPy.kern.exponential(1), GPy.kern.Matern32(1), GPy.kern.Matern52(1), GPy.kern.Brownian(1), GPy.kern.bias(1), GPy.kern.linear(1), GPy.kern.spline(1), GPy.kern.periodic_exponential(1), GPy.kern.periodic_Matern32(1), GPy.kern.periodic_Matern52(1), GPy.kern.white(1)] -kernel_names = ["GPy.kern.rbf", "GPy.kern.exponential", "GPy.kern.Matern32", "GPy.kern.Matern52", "GPy.kern.Brownian", "GPy.kern.bias", "GPy.kern.linear", "GPy.kern.spline", "GPy.kern.periodic_exponential", "GPy.kern.periodic_Matern32", "GPy.kern.periodic_Matern52", "GPy.kern.white"] - -pb.figure(figsize=(16,12)) -pb.subplots_adjust(wspace=.5, hspace=.5) -for i, kern in enumerate(kernels): - pb.subplot(3,4,i+1) - kern.plot(x=7.5,plot_limits=[0.00001,15.]) - pb.title(kernel_names[i]+ '\n') - -# actual plot for the noise -i = 11 -X = np.linspace(0.,15.,201) -WN = 0*X -WN[100] = 1. -pb.subplot(3,4,i+1) -pb.plot(X,WN,'b') diff --git a/GPy/examples/tutorials.py b/GPy/examples/tutorials.py new file mode 100644 index 00000000..be550e01 --- /dev/null +++ b/GPy/examples/tutorials.py @@ -0,0 +1,201 @@ +# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + + +""" +Code of Tutorials +""" + +def tuto_GP_regression(): + """The detailed explanations of the commands used in this file can be found in the tutorial section""" + + import pylab as pb + pb.ion() + import numpy as np + import GPy + + X = np.random.uniform(-3.,3.,(20,1)) + Y = np.sin(X) + np.random.randn(20,1)*0.05 + + kernel = GPy.kern.rbf(D=1, variance=1., lengthscale=1.) + + m = GPy.models.GP_regression(X,Y,kernel) + + print m + m.plot() + + m.constrain_positive('') + + m.unconstrain('') # Required to remove the previous constrains + m.constrain_positive('rbf_variance') + m.constrain_bounded('lengthscale',1.,10. ) + m.constrain_fixed('noise',0.0025) + + m.optimize() + + m.optimize_restarts(Nrestarts = 10) + + ########################### + # 2-dimensional example # + ########################### + + import pylab as pb + pb.ion() + import numpy as np + import GPy + + # sample inputs and outputs + X = np.random.uniform(-3.,3.,(50,2)) + Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(50,1)*0.05 + + # define kernel + ker = GPy.kern.Matern52(2,ARD=True) + GPy.kern.white(2) + + # create simple GP model + m = GPy.models.GP_regression(X,Y,ker) + + # contrain all parameters to be positive + m.constrain_positive('') + + # optimize and plot + pb.figure() + m.optimize('tnc', max_f_eval = 1000) + + m.plot() + print(m) + + +def tuto_kernel_overview(): + """The detailed explanations of the commands used in this file can be found in the tutorial section""" + import pylab as pb + import numpy as np + import GPy + pb.ion() + + ker1 = GPy.kern.rbf(1) # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.) + ker2 = GPy.kern.rbf(D=1, variance = .75, lengthscale=2.) + ker3 = GPy.kern.rbf(1, .5, .5) + + print ker2 + ker1.plot() + ker2.plot() + ker3.plot() + + k1 = GPy.kern.rbf(1,1.,2.) + k2 = GPy.kern.Matern32(1, 0.5, 0.2) + + # Product of kernels + k_prod = k1.prod(k2) + k_prodorth = k1.prod_orthogonal(k2) + + # Sum of kernels + k_add = k1.add(k2) + k_addorth = k1.add_orthogonal(k2) + + pb.figure(figsize=(8,8)) + pb.subplot(2,2,1) + k_prod.plot() + pb.title('prod') + pb.subplot(2,2,2) + k_prodorth.plot() + pb.title('prod_orthogonal') + pb.subplot(2,2,3) + k_add.plot() + pb.title('add') + pb.subplot(2,2,4) + k_addorth.plot() + pb.title('add_orthogonal') + pb.subplots_adjust(wspace=0.3, hspace=0.3) + + k1 = GPy.kern.rbf(1,1.,2) + k2 = GPy.kern.periodic_Matern52(1,variance=1e3, lengthscale=1, period = 1.5, lower=-5., upper = 5) + + k = k1 * k2 # equivalent to k = k1.prod(k2) + print k + + # Simulate sample paths + X = np.linspace(-5,5,501)[:,None] + Y = np.random.multivariate_normal(np.zeros(501),k.K(X),1) + + # plot + pb.figure(figsize=(10,4)) + pb.subplot(1,2,1) + k.plot() + pb.subplot(1,2,2) + pb.plot(X,Y.T) + pb.ylabel("Sample path") + pb.subplots_adjust(wspace=0.3) + + k = (k1+k2)*(k1+k2) + print k.parts[0].name, '\n', k.parts[1].name, '\n', k.parts[2].name, '\n', k.parts[3].name + + k1 = GPy.kern.rbf(1) + k2 = GPy.kern.Matern32(1) + k3 = GPy.kern.white(1) + + k = k1 + k2 + k3 + print k + + k.constrain_positive('var') + k.constrain_fixed(np.array([1]),1.75) + k.tie_param('len') + k.unconstrain('white') + k.constrain_bounded('white',lower=1e-5,upper=.5) + print k + + k_cst = GPy.kern.bias(1,variance=1.) + k_mat = GPy.kern.Matern52(1,variance=1., lengthscale=3) + Kanova = (k_cst + k_mat).prod_orthogonal(k_cst + k_mat) + print Kanova + + # sample inputs and outputs + X = np.random.uniform(-3.,3.,(40,2)) + Y = 0.5*X[:,:1] + 0.5*X[:,1:] + 2*np.sin(X[:,:1]) * np.sin(X[:,1:]) + + # Create GP regression model + m = GPy.models.GP_regression(X,Y,Kanova) + pb.figure(figsize=(5,5)) + m.plot() + + pb.figure(figsize=(20,3)) + pb.subplots_adjust(wspace=0.5) + pb.subplot(1,5,1) + m.plot() + pb.subplot(1,5,2) + pb.ylabel("= ",rotation='horizontal',fontsize='30') + pb.subplot(1,5,3) + m.plot(which_functions=[False,True,False,False]) + pb.ylabel("cst +",rotation='horizontal',fontsize='30') + pb.subplot(1,5,4) + m.plot(which_functions=[False,False,True,False]) + pb.ylabel("+ ",rotation='horizontal',fontsize='30') + pb.subplot(1,5,5) + pb.ylabel("+ ",rotation='horizontal',fontsize='30') + m.plot(which_functions=[False,False,False,True]) + + ker1 = GPy.kern.rbf(D=1) # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.) + ker2 = GPy.kern.rbf(D=1, variance = .75, lengthscale=3.) + ker3 = GPy.kern.rbf(1, .5, .25) + + ker1.plot() + ker2.plot() + ker3.plot() + #pb.savefig("Figures/tuto_kern_overview_basicdef.png") + + kernels = [GPy.kern.rbf(1), GPy.kern.exponential(1), GPy.kern.Matern32(1), GPy.kern.Matern52(1), GPy.kern.Brownian(1), GPy.kern.bias(1), GPy.kern.linear(1), GPy.kern.spline(1), GPy.kern.periodic_exponential(1), GPy.kern.periodic_Matern32(1), GPy.kern.periodic_Matern52(1), GPy.kern.white(1)] + kernel_names = ["GPy.kern.rbf", "GPy.kern.exponential", "GPy.kern.Matern32", "GPy.kern.Matern52", "GPy.kern.Brownian", "GPy.kern.bias", "GPy.kern.linear", "GPy.kern.spline", "GPy.kern.periodic_exponential", "GPy.kern.periodic_Matern32", "GPy.kern.periodic_Matern52", "GPy.kern.white"] + + pb.figure(figsize=(16,12)) + pb.subplots_adjust(wspace=.5, hspace=.5) + for i, kern in enumerate(kernels): + pb.subplot(3,4,i+1) + kern.plot(x=7.5,plot_limits=[0.00001,15.]) + pb.title(kernel_names[i]+ '\n') + + # actual plot for the noise + i = 11 + X = np.linspace(0.,15.,201) + WN = 0*X + WN[100] = 1. + pb.subplot(3,4,i+1) + pb.plot(X,WN,'b') diff --git a/doc/tuto_GP_regression.rst b/doc/tuto_GP_regression.rst index 9de79a8c..24e10528 100644 --- a/doc/tuto_GP_regression.rst +++ b/doc/tuto_GP_regression.rst @@ -2,7 +2,7 @@ Gaussian process regression tutorial ************************************* -We will see in this tutorial the basics for building a 1 dimensional and a 2 dimensional Gaussian process regression model, also known as a kriging model. The code shown in this tutorial can be found without the comments at GPy/examples/tuto_GP_regression.py. +We will see in this tutorial the basics for building a 1 dimensional and a 2 dimensional Gaussian process regression model, also known as a kriging model. The code shown in this tutorial can be obtained at GPy/examples/tutorials.py, or by running ``GPy.examples.tutorials.tuto_GP_regression()``. We first import the libraries we will need: :: diff --git a/doc/tuto_kernel_overview.rst b/doc/tuto_kernel_overview.rst index e410696a..dfb7fb3f 100644 --- a/doc/tuto_kernel_overview.rst +++ b/doc/tuto_kernel_overview.rst @@ -2,7 +2,7 @@ **************************** tutorial : A kernel overview **************************** -The aim of this tutorial is to give a better understanding of the kernel objects in GPy and to list the ones that are already implemented. The code shown in this tutorial can be found without the comments at GPy/examples/tuto_kernel_overview.py. +The aim of this tutorial is to give a better understanding of the kernel objects in GPy and to list the ones that are already implemented. The code shown in this tutorial can be obtained at GPy/examples/tutorials.py or by running ``GPy.examples.tutorials.tuto_kernel_overview()``. First we import the libraries we will need :: From 6a330db25336dad8c2668fff0cb7fc8430be41e0 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 11 Mar 2013 13:26:39 +0000 Subject: [PATCH 15/18] americanized spellings --- GPy/likelihoods/Gaussian.py | 4 +- GPy/models/GP.py | 15 ++-- GPy/models/sparse_GP.py | 4 +- GPy/models/uncollapsed_sparse_GP.py | 2 +- GPy/util/Tango.py | 132 ++++++++++++++-------------- GPy/util/datasets.py | 2 +- GPy/util/plot.py | 2 +- 7 files changed, 80 insertions(+), 81 deletions(-) diff --git a/GPy/likelihoods/Gaussian.py b/GPy/likelihoods/Gaussian.py index 4a32f066..a5084cc0 100644 --- a/GPy/likelihoods/Gaussian.py +++ b/GPy/likelihoods/Gaussian.py @@ -8,7 +8,7 @@ class Gaussian(likelihood): self.Z = 0. # a correction factor which accounts for the approximation made N, self.D = data.shape - #normalisation + #normaliztion if normalize: self._mean = data.mean(0)[None,:] self._std = data.std(0)[None,:] @@ -45,7 +45,7 @@ class Gaussian(likelihood): def predictive_values(self,mu,var): """ - Un-normalise the prediction and add the likelihood variance, then return the 5%, 95% interval + Un-normalize the prediction and add the likelihood variance, then return the 5%, 95% interval """ mean = mu*self._std + self._mean true_var = (var + self._variance)*self._std**2 diff --git a/GPy/models/GP.py b/GPy/models/GP.py index 5879a2bf..796ab7d6 100644 --- a/GPy/models/GP.py +++ b/GPy/models/GP.py @@ -30,7 +30,6 @@ class GP(model): .. Note:: Multiple independent outputs are allowed using columns of Y """ - #FIXME normalize vs normalise def __init__(self, X, likelihood, kernel, normalize_X=False, Xslices=None): # parse arguments @@ -41,7 +40,7 @@ class GP(model): assert isinstance(kernel, kern.kern) self.kern = kernel - #here's some simple normalisation for the inputs + #here's some simple normalization for the inputs if normalize_X: self._Xmean = X.mean(0)[None,:] self._Xstd = X.std(0)[None,:] @@ -134,7 +133,7 @@ class GP(model): def _raw_predict(self,_Xnew,slices=None, full_cov=False): """ Internal helper function for making predictions, does not account - for normalisation or likelihood + for normalization or likelihood """ Kx = self.kern.K(self.X,_Xnew, slices1=self.Xslices,slices2=slices) mu = np.dot(np.dot(Kx.T,self.Ki),self.likelihood.Y) @@ -172,10 +171,10 @@ class GP(model): - If a list of booleans, specifying which kernel parts are active If full_cov and self.D > 1, the return shape of var is Nnew x Nnew x self.D. If self.D == 1, the return shape is Nnew x Nnew. - This is to allow for different normalisations of the output dimensions. + This is to allow for different normalizations of the output dimensions. """ - #normalise X values + #normalize X values Xnew = (Xnew.copy() - self._Xmean) / self._Xstd mu, var = self._raw_predict(Xnew, slices, full_cov) @@ -187,7 +186,7 @@ class GP(model): def plot_f(self, samples=0, plot_limits=None, which_data='all', which_functions='all', resolution=None, full_cov=False): """ - Plot the GP's view of the world, where the data is normalised and the likelihood is Gaussian + Plot the GP's view of the world, where the data is normalized and the likelihood is Gaussian :param samples: the number of a posteriori samples to plot :param which_data: which if the training data to plot (default all) @@ -203,7 +202,7 @@ class GP(model): - In higher dimensions, we've no implemented this yet !TODO! Can plot only part of the data and part of the posterior functions using which_data and which_functions - Plot the data's view of the world, with non-normalised values and GP predictions passed through the likelihood + Plot the data's view of the world, with non-normalized values and GP predictions passed through the likelihood """ if which_functions=='all': which_functions = [True]*self.kern.Nparts @@ -221,7 +220,7 @@ class GP(model): Ysim = np.random.multivariate_normal(m.flatten(),v,samples) gpplot(Xnew,m,m-2*np.sqrt(np.diag(v)[:,None]),m+2*np.sqrt(np.diag(v))[:,None]) for i in range(samples): - pb.plot(Xnew,Ysim[i,:],Tango.coloursHex['darkBlue'],linewidth=0.25) + pb.plot(Xnew,Ysim[i,:],Tango.colorsHex['darkBlue'],linewidth=0.25) pb.plot(self.X[which_data],self.likelihood.Y[which_data],'kx',mew=1.5) pb.xlim(xmin,xmax) ymin,ymax = min(np.append(self.likelihood.Y,m-2*np.sqrt(np.diag(v)[:,None]))), max(np.append(self.likelihood.Y,m+2*np.sqrt(np.diag(v)[:,None]))) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index ff00faea..acf4f6c0 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -54,7 +54,7 @@ class sparse_GP(GP): GP.__init__(self, X, likelihood, kernel=kernel, normalize_X=normalize_X, Xslices=Xslices) - #normalise X uncertainty also + #normalize X uncertainty also if self.has_uncertain_inputs: self.X_uncertainty /= np.square(self._Xstd) @@ -228,7 +228,7 @@ class sparse_GP(GP): return dL_dZ def _raw_predict(self, Xnew, slices, full_cov=False): - """Internal helper function for making predictions, does not account for normalisation""" + """Internal helper function for making predictions, does not account for normalization""" Kx = self.kern.K(self.Z, Xnew) mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V) diff --git a/GPy/models/uncollapsed_sparse_GP.py b/GPy/models/uncollapsed_sparse_GP.py index 43624e72..d2638784 100644 --- a/GPy/models/uncollapsed_sparse_GP.py +++ b/GPy/models/uncollapsed_sparse_GP.py @@ -93,7 +93,7 @@ class uncollapsed_sparse_GP(sparse_GP): return A+B+C+D+E def _raw_predict(self, Xnew, slices,full_cov=False): - """Internal helper function for making predictions, does not account for normalisation""" + """Internal helper function for making predictions, does not account for normalization""" Kx = self.kern.K(Xnew,self.Z) mu = mdot(Kx,self.Kmmi,self.q_u_expectation[0]) diff --git a/GPy/util/Tango.py b/GPy/util/Tango.py index 8035ffe6..eeb2e075 100644 --- a/GPy/util/Tango.py +++ b/GPy/util/Tango.py @@ -25,7 +25,7 @@ def fewerXticks(ax=None,divideby=2): ax.set_xticks(ax.get_xticks()[::divideby]) -coloursHex = {\ +colorsHex = {\ "Aluminium6":"#2e3436",\ "Aluminium5":"#555753",\ "Aluminium4":"#888a85",\ @@ -54,9 +54,9 @@ coloursHex = {\ "mediumButter":"#edd400",\ "darkButter":"#c4a000"} -darkList = [coloursHex['darkBlue'],coloursHex['darkRed'],coloursHex['darkGreen'], coloursHex['darkOrange'], coloursHex['darkButter'], coloursHex['darkPurple'], coloursHex['darkChocolate'], coloursHex['Aluminium6']] -mediumList = [coloursHex['mediumBlue'], coloursHex['mediumRed'],coloursHex['mediumGreen'], coloursHex['mediumOrange'], coloursHex['mediumButter'], coloursHex['mediumPurple'], coloursHex['mediumChocolate'], coloursHex['Aluminium5']] -lightList = [coloursHex['lightBlue'], coloursHex['lightRed'],coloursHex['lightGreen'], coloursHex['lightOrange'], coloursHex['lightButter'], coloursHex['lightPurple'], coloursHex['lightChocolate'], coloursHex['Aluminium4']] +darkList = [colorsHex['darkBlue'],colorsHex['darkRed'],colorsHex['darkGreen'], colorsHex['darkOrange'], colorsHex['darkButter'], colorsHex['darkPurple'], colorsHex['darkChocolate'], colorsHex['Aluminium6']] +mediumList = [colorsHex['mediumBlue'], colorsHex['mediumRed'],colorsHex['mediumGreen'], colorsHex['mediumOrange'], colorsHex['mediumButter'], colorsHex['mediumPurple'], colorsHex['mediumChocolate'], colorsHex['Aluminium5']] +lightList = [colorsHex['lightBlue'], colorsHex['lightRed'],colorsHex['lightGreen'], colorsHex['lightOrange'], colorsHex['lightButter'], colorsHex['lightPurple'], colorsHex['lightChocolate'], colorsHex['Aluminium4']] def currentDark(): return darkList[-1] @@ -76,85 +76,85 @@ def nextLight(): return lightList[-1] def reset(): - while not darkList[0]==coloursHex['darkBlue']: + while not darkList[0]==colorsHex['darkBlue']: darkList.append(darkList.pop(0)) - while not mediumList[0]==coloursHex['mediumBlue']: + while not mediumList[0]==colorsHex['mediumBlue']: mediumList.append(mediumList.pop(0)) - while not lightList[0]==coloursHex['lightBlue']: + while not lightList[0]==colorsHex['lightBlue']: lightList.append(lightList.pop(0)) def setLightFigures(): - mpl.rcParams['axes.edgecolor']=coloursHex['Aluminium6'] - mpl.rcParams['axes.facecolor']=coloursHex['Aluminium2'] - mpl.rcParams['axes.labelcolor']=coloursHex['Aluminium6'] - mpl.rcParams['figure.edgecolor']=coloursHex['Aluminium6'] - mpl.rcParams['figure.facecolor']=coloursHex['Aluminium2'] - mpl.rcParams['grid.color']=coloursHex['Aluminium6'] - mpl.rcParams['savefig.edgecolor']=coloursHex['Aluminium2'] - mpl.rcParams['savefig.facecolor']=coloursHex['Aluminium2'] - mpl.rcParams['text.color']=coloursHex['Aluminium6'] - mpl.rcParams['xtick.color']=coloursHex['Aluminium6'] - mpl.rcParams['ytick.color']=coloursHex['Aluminium6'] + mpl.rcParams['axes.edgecolor']=colorsHex['Aluminium6'] + mpl.rcParams['axes.facecolor']=colorsHex['Aluminium2'] + mpl.rcParams['axes.labelcolor']=colorsHex['Aluminium6'] + mpl.rcParams['figure.edgecolor']=colorsHex['Aluminium6'] + mpl.rcParams['figure.facecolor']=colorsHex['Aluminium2'] + mpl.rcParams['grid.color']=colorsHex['Aluminium6'] + mpl.rcParams['savefig.edgecolor']=colorsHex['Aluminium2'] + mpl.rcParams['savefig.facecolor']=colorsHex['Aluminium2'] + mpl.rcParams['text.color']=colorsHex['Aluminium6'] + mpl.rcParams['xtick.color']=colorsHex['Aluminium6'] + mpl.rcParams['ytick.color']=colorsHex['Aluminium6'] def setDarkFigures(): - mpl.rcParams['axes.edgecolor']=coloursHex['Aluminium2'] - mpl.rcParams['axes.facecolor']=coloursHex['Aluminium6'] - mpl.rcParams['axes.labelcolor']=coloursHex['Aluminium2'] - mpl.rcParams['figure.edgecolor']=coloursHex['Aluminium2'] - mpl.rcParams['figure.facecolor']=coloursHex['Aluminium6'] - mpl.rcParams['grid.color']=coloursHex['Aluminium2'] - mpl.rcParams['savefig.edgecolor']=coloursHex['Aluminium6'] - mpl.rcParams['savefig.facecolor']=coloursHex['Aluminium6'] - mpl.rcParams['text.color']=coloursHex['Aluminium2'] - mpl.rcParams['xtick.color']=coloursHex['Aluminium2'] - mpl.rcParams['ytick.color']=coloursHex['Aluminium2'] + mpl.rcParams['axes.edgecolor']=colorsHex['Aluminium2'] + mpl.rcParams['axes.facecolor']=colorsHex['Aluminium6'] + mpl.rcParams['axes.labelcolor']=colorsHex['Aluminium2'] + mpl.rcParams['figure.edgecolor']=colorsHex['Aluminium2'] + mpl.rcParams['figure.facecolor']=colorsHex['Aluminium6'] + mpl.rcParams['grid.color']=colorsHex['Aluminium2'] + mpl.rcParams['savefig.edgecolor']=colorsHex['Aluminium6'] + mpl.rcParams['savefig.facecolor']=colorsHex['Aluminium6'] + mpl.rcParams['text.color']=colorsHex['Aluminium2'] + mpl.rcParams['xtick.color']=colorsHex['Aluminium2'] + mpl.rcParams['ytick.color']=colorsHex['Aluminium2'] def hex2rgb(hexcolor): hexcolor = [hexcolor[1+2*i:1+2*(i+1)] for i in range(3)] r,g,b = [int(n,16) for n in hexcolor] return (r,g,b) -coloursRGB = dict([(k,hex2rgb(i)) for k,i in coloursHex.items()]) +colorsRGB = dict([(k,hex2rgb(i)) for k,i in colorsHex.items()]) -cdict_RB = {'red' :((0.,coloursRGB['mediumRed'][0]/256.,coloursRGB['mediumRed'][0]/256.), - (.5,coloursRGB['mediumPurple'][0]/256.,coloursRGB['mediumPurple'][0]/256.), - (1.,coloursRGB['mediumBlue'][0]/256.,coloursRGB['mediumBlue'][0]/256.)), - 'green':((0.,coloursRGB['mediumRed'][1]/256.,coloursRGB['mediumRed'][1]/256.), - (.5,coloursRGB['mediumPurple'][1]/256.,coloursRGB['mediumPurple'][1]/256.), - (1.,coloursRGB['mediumBlue'][1]/256.,coloursRGB['mediumBlue'][1]/256.)), - 'blue':((0.,coloursRGB['mediumRed'][2]/256.,coloursRGB['mediumRed'][2]/256.), - (.5,coloursRGB['mediumPurple'][2]/256.,coloursRGB['mediumPurple'][2]/256.), - (1.,coloursRGB['mediumBlue'][2]/256.,coloursRGB['mediumBlue'][2]/256.))} +cdict_RB = {'red' :((0.,colorsRGB['mediumRed'][0]/256.,colorsRGB['mediumRed'][0]/256.), + (.5,colorsRGB['mediumPurple'][0]/256.,colorsRGB['mediumPurple'][0]/256.), + (1.,colorsRGB['mediumBlue'][0]/256.,colorsRGB['mediumBlue'][0]/256.)), + 'green':((0.,colorsRGB['mediumRed'][1]/256.,colorsRGB['mediumRed'][1]/256.), + (.5,colorsRGB['mediumPurple'][1]/256.,colorsRGB['mediumPurple'][1]/256.), + (1.,colorsRGB['mediumBlue'][1]/256.,colorsRGB['mediumBlue'][1]/256.)), + 'blue':((0.,colorsRGB['mediumRed'][2]/256.,colorsRGB['mediumRed'][2]/256.), + (.5,colorsRGB['mediumPurple'][2]/256.,colorsRGB['mediumPurple'][2]/256.), + (1.,colorsRGB['mediumBlue'][2]/256.,colorsRGB['mediumBlue'][2]/256.))} -cdict_BGR = {'red' :((0.,coloursRGB['mediumBlue'][0]/256.,coloursRGB['mediumBlue'][0]/256.), - (.5,coloursRGB['mediumGreen'][0]/256.,coloursRGB['mediumGreen'][0]/256.), - (1.,coloursRGB['mediumRed'][0]/256.,coloursRGB['mediumRed'][0]/256.)), - 'green':((0.,coloursRGB['mediumBlue'][1]/256.,coloursRGB['mediumBlue'][1]/256.), - (.5,coloursRGB['mediumGreen'][1]/256.,coloursRGB['mediumGreen'][1]/256.), - (1.,coloursRGB['mediumRed'][1]/256.,coloursRGB['mediumRed'][1]/256.)), - 'blue':((0.,coloursRGB['mediumBlue'][2]/256.,coloursRGB['mediumBlue'][2]/256.), - (.5,coloursRGB['mediumGreen'][2]/256.,coloursRGB['mediumGreen'][2]/256.), - (1.,coloursRGB['mediumRed'][2]/256.,coloursRGB['mediumRed'][2]/256.))} +cdict_BGR = {'red' :((0.,colorsRGB['mediumBlue'][0]/256.,colorsRGB['mediumBlue'][0]/256.), + (.5,colorsRGB['mediumGreen'][0]/256.,colorsRGB['mediumGreen'][0]/256.), + (1.,colorsRGB['mediumRed'][0]/256.,colorsRGB['mediumRed'][0]/256.)), + 'green':((0.,colorsRGB['mediumBlue'][1]/256.,colorsRGB['mediumBlue'][1]/256.), + (.5,colorsRGB['mediumGreen'][1]/256.,colorsRGB['mediumGreen'][1]/256.), + (1.,colorsRGB['mediumRed'][1]/256.,colorsRGB['mediumRed'][1]/256.)), + 'blue':((0.,colorsRGB['mediumBlue'][2]/256.,colorsRGB['mediumBlue'][2]/256.), + (.5,colorsRGB['mediumGreen'][2]/256.,colorsRGB['mediumGreen'][2]/256.), + (1.,colorsRGB['mediumRed'][2]/256.,colorsRGB['mediumRed'][2]/256.))} -cdict_Alu = {'red' :((0./5,coloursRGB['Aluminium1'][0]/256.,coloursRGB['Aluminium1'][0]/256.), - (1./5,coloursRGB['Aluminium2'][0]/256.,coloursRGB['Aluminium2'][0]/256.), - (2./5,coloursRGB['Aluminium3'][0]/256.,coloursRGB['Aluminium3'][0]/256.), - (3./5,coloursRGB['Aluminium4'][0]/256.,coloursRGB['Aluminium4'][0]/256.), - (4./5,coloursRGB['Aluminium5'][0]/256.,coloursRGB['Aluminium5'][0]/256.), - (5./5,coloursRGB['Aluminium6'][0]/256.,coloursRGB['Aluminium6'][0]/256.)), - 'green' :((0./5,coloursRGB['Aluminium1'][1]/256.,coloursRGB['Aluminium1'][1]/256.), - (1./5,coloursRGB['Aluminium2'][1]/256.,coloursRGB['Aluminium2'][1]/256.), - (2./5,coloursRGB['Aluminium3'][1]/256.,coloursRGB['Aluminium3'][1]/256.), - (3./5,coloursRGB['Aluminium4'][1]/256.,coloursRGB['Aluminium4'][1]/256.), - (4./5,coloursRGB['Aluminium5'][1]/256.,coloursRGB['Aluminium5'][1]/256.), - (5./5,coloursRGB['Aluminium6'][1]/256.,coloursRGB['Aluminium6'][1]/256.)), - 'blue' :((0./5,coloursRGB['Aluminium1'][2]/256.,coloursRGB['Aluminium1'][2]/256.), - (1./5,coloursRGB['Aluminium2'][2]/256.,coloursRGB['Aluminium2'][2]/256.), - (2./5,coloursRGB['Aluminium3'][2]/256.,coloursRGB['Aluminium3'][2]/256.), - (3./5,coloursRGB['Aluminium4'][2]/256.,coloursRGB['Aluminium4'][2]/256.), - (4./5,coloursRGB['Aluminium5'][2]/256.,coloursRGB['Aluminium5'][2]/256.), - (5./5,coloursRGB['Aluminium6'][2]/256.,coloursRGB['Aluminium6'][2]/256.))} +cdict_Alu = {'red' :((0./5,colorsRGB['Aluminium1'][0]/256.,colorsRGB['Aluminium1'][0]/256.), + (1./5,colorsRGB['Aluminium2'][0]/256.,colorsRGB['Aluminium2'][0]/256.), + (2./5,colorsRGB['Aluminium3'][0]/256.,colorsRGB['Aluminium3'][0]/256.), + (3./5,colorsRGB['Aluminium4'][0]/256.,colorsRGB['Aluminium4'][0]/256.), + (4./5,colorsRGB['Aluminium5'][0]/256.,colorsRGB['Aluminium5'][0]/256.), + (5./5,colorsRGB['Aluminium6'][0]/256.,colorsRGB['Aluminium6'][0]/256.)), + 'green' :((0./5,colorsRGB['Aluminium1'][1]/256.,colorsRGB['Aluminium1'][1]/256.), + (1./5,colorsRGB['Aluminium2'][1]/256.,colorsRGB['Aluminium2'][1]/256.), + (2./5,colorsRGB['Aluminium3'][1]/256.,colorsRGB['Aluminium3'][1]/256.), + (3./5,colorsRGB['Aluminium4'][1]/256.,colorsRGB['Aluminium4'][1]/256.), + (4./5,colorsRGB['Aluminium5'][1]/256.,colorsRGB['Aluminium5'][1]/256.), + (5./5,colorsRGB['Aluminium6'][1]/256.,colorsRGB['Aluminium6'][1]/256.)), + 'blue' :((0./5,colorsRGB['Aluminium1'][2]/256.,colorsRGB['Aluminium1'][2]/256.), + (1./5,colorsRGB['Aluminium2'][2]/256.,colorsRGB['Aluminium2'][2]/256.), + (2./5,colorsRGB['Aluminium3'][2]/256.,colorsRGB['Aluminium3'][2]/256.), + (3./5,colorsRGB['Aluminium4'][2]/256.,colorsRGB['Aluminium4'][2]/256.), + (4./5,colorsRGB['Aluminium5'][2]/256.,colorsRGB['Aluminium5'][2]/256.), + (5./5,colorsRGB['Aluminium6'][2]/256.,colorsRGB['Aluminium6'][2]/256.))} # cmap_Alu = mpl.colors.LinearSegmentedColormap('TangoAluminium',cdict_Alu,256) # cmap_BGR = mpl.colors.LinearSegmentedColormap('TangoRedBlue',cdict_BGR,256) # cmap_RB = mpl.colors.LinearSegmentedColormap('TangoRedBlue',cdict_RB,256) diff --git a/GPy/util/datasets.py b/GPy/util/datasets.py index 5506fbef..a8ec2539 100644 --- a/GPy/util/datasets.py +++ b/GPy/util/datasets.py @@ -46,7 +46,7 @@ def oil_100(seed=default_seed): return {'X': X, 'Y': Y, 'info': "Subsample of the oil data extracting 100 values randomly without replacement."} def pumadyn(seed=default_seed): - # Data is variance 1, no need to normalise. + # Data is variance 1, no need to normalize. data = np.loadtxt(os.path.join(data_path, 'pumadyn-32nm/Dataset.data.gz')) indices = np.random.permutation(data.shape[0]) indicesTrain = indices[0:7168] diff --git a/GPy/util/plot.py b/GPy/util/plot.py index 8e71764d..295047b1 100644 --- a/GPy/util/plot.py +++ b/GPy/util/plot.py @@ -6,7 +6,7 @@ import Tango import pylab as pb import numpy as np -def gpplot(x,mu,lower,upper,edgecol=Tango.coloursHex['darkBlue'],fillcol=Tango.coloursHex['lightBlue'],axes=None,**kwargs): +def gpplot(x,mu,lower,upper,edgecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue'],axes=None,**kwargs): if axes is None: axes = pb.gca() mu = mu.flatten() From 417dac0080a8f8699d3d6e15662d94c84a783d59 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Mon, 11 Mar 2013 13:54:47 +0000 Subject: [PATCH 16/18] Adding testing file for examples --- GPy/testing/examples_tests.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 GPy/testing/examples_tests.py diff --git a/GPy/testing/examples_tests.py b/GPy/testing/examples_tests.py new file mode 100644 index 00000000..dd85ea34 --- /dev/null +++ b/GPy/testing/examples_tests.py @@ -0,0 +1,26 @@ +# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import unittest +import numpy as np +import GPy + +class ExamplesTests(unittest.TestCase): + def test_check_model_returned(self): + pass + + def test_model_checkgrads(self): + pass + + def test_all_examples(self): + #Load models + + #Loop through models + for model in models: + + self.assertTrue(m.checkgrad()) + + +if __name__ == "__main__": + print "Running unit tests, please be (very) patient..." + unittest.main() From 84119a19b3e82d84d201399404758fb9c8396839 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Mon, 11 Mar 2013 13:55:43 +0000 Subject: [PATCH 17/18] Skipping tests --- GPy/testing/examples_tests.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/GPy/testing/examples_tests.py b/GPy/testing/examples_tests.py index dd85ea34..25cfad04 100644 --- a/GPy/testing/examples_tests.py +++ b/GPy/testing/examples_tests.py @@ -13,12 +13,12 @@ class ExamplesTests(unittest.TestCase): pass def test_all_examples(self): + pass #Load models #Loop through models - for model in models: - - self.assertTrue(m.checkgrad()) + #for model in models: + #self.assertTrue(m.checkgrad()) if __name__ == "__main__": From 05ca5cfe6d3788f065da2e8f420fe5ba301a021a Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Mon, 11 Mar 2013 14:03:23 +0000 Subject: [PATCH 18/18] working on psi cross terms --- GPy/kern/kern.py | 8 ++++---- GPy/testing/bgplvm_tests.py | 26 ++++++++++++++++++++++++++ setup.py | 2 -- 3 files changed, 30 insertions(+), 6 deletions(-) diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py index 87e67f33..8cadf662 100644 --- a/GPy/kern/kern.py +++ b/GPy/kern/kern.py @@ -444,9 +444,9 @@ class kern(parameterised): pass #rbf X bias elif p1.name=='bias' and p2.name=='rbf': - target += p2.dpsi1_dX(dL_dpsi2.sum(1)*p1.variance,Z,mu,S,target) + p2.dpsi1_dX(dL_dpsi2.sum(1)*p1.variance,Z,mu,S,target) elif p2.name=='bias' and p1.name=='rbf': - target += p1.dpsi1_dZ(dL_dpsi2.sum(2)*p2.variance,Z,mu,S,target) + p1.dpsi1_dZ(dL_dpsi2.sum(2)*p2.variance,Z,mu,S,target) #rbf X linear elif p1.name=='linear' and p2.name=='rbf': raise NotImplementedError #TODO @@ -471,9 +471,9 @@ class kern(parameterised): pass #rbf X bias elif p1.name=='bias' and p2.name=='rbf': - target += p2.dpsi1_dmuS(partial.sum(1)*p1.variance,Z,mu,S,target_mu,target_S) + p2.dpsi1_dmuS(partial.sum(1)*p1.variance,Z,mu,S,target_mu,target_S) elif p2.name=='bias' and p1.name=='rbf': - target += p1.dpsi1_dmuS(partial.sum(2)*p2.variance,Z,mu,S,target_mu,target_S) + p1.dpsi1_dmuS(partial.sum(2)*p2.variance,Z,mu,S,target_mu,target_S) #rbf X linear elif p1.name=='linear' and p2.name=='rbf': raise NotImplementedError #TODO diff --git a/GPy/testing/bgplvm_tests.py b/GPy/testing/bgplvm_tests.py index e3bd2b36..80e6fecd 100644 --- a/GPy/testing/bgplvm_tests.py +++ b/GPy/testing/bgplvm_tests.py @@ -45,6 +45,32 @@ class BGPLVMTests(unittest.TestCase): m.randomize() self.assertTrue(m.checkgrad()) + def test_rbf_bias_kern(self): + N, M, Q, D = 10, 3, 2, 4 + X = np.random.rand(N, Q) + k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(N),K,D).T + Y -= Y.mean(axis=0) + k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) + m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M) + m.constrain_positive('(rbf|bias|noise|white|S)') + m.randomize() + self.assertTrue(m.checkgrad()) + + def test_linear_bias_kern(self): + N, M, Q, D = 10, 3, 2, 4 + X = np.random.rand(N, Q) + k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) + K = k.K(X) + Y = np.random.multivariate_normal(np.zeros(N),K,D).T + Y -= Y.mean(axis=0) + k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) + m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M=M) + m.constrain_positive('(linear|bias|noise|white|S)') + m.randomize() + self.assertTrue(m.checkgrad()) + if __name__ == "__main__": print "Running unit tests, please be (very) patient..." diff --git a/setup.py b/setup.py index b701b74d..ca193fbc 100644 --- a/setup.py +++ b/setup.py @@ -32,7 +32,5 @@ setup(name = 'GPy', #setup_requires=['sphinx'], #cmdclass = {'build_sphinx': BuildDoc}, classifiers=[ - "Development Status :: 1 - Alpha", - "Topic :: Machine Learning", "License :: OSI Approved :: BSD License"], )