From de0a5d0e70643ddd4a2d2901c740041af81ca981 Mon Sep 17 00:00:00 2001 From: Neil Lawrence Date: Wed, 9 Oct 2013 12:07:39 +0100 Subject: [PATCH] Some fixes and changes to the sympykern. --- GPy/kern/constructors.py | 17 ++++++++++------- GPy/kern/kern.py | 10 +++++----- GPy/kern/parts/kernpart.py | 2 -- GPy/kern/parts/sympykern.py | 22 ++++++++++++---------- GPy/testing/kernel_tests.py | 2 +- 5 files changed, 28 insertions(+), 25 deletions(-) diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index a1252052..62c29744 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -317,20 +317,23 @@ if sympy_available: """ Exponentiated quadratic with multiple outputs. """ - X = sp.symbols('x_:' + str(input_dim)) - Z = sp.symbols('z_:' + str(input_dim)) + real_input_dim = input_dim + if output_dim>1: + real_input_dim -= 1 + X = sp.symbols('x_:' + str(real_input_dim)) + Z = sp.symbols('z_:' + str(real_input_dim)) variance = sp.var('variance',positive=True) if ARD: - lengthscales = [sp.var('lengthscale%i_i lengthscale%i_j' % i, positive=True) for i in range(input_dim)] - dist_string = ' + '.join(['(x_%i-z_%i)**2/(lengthscale%i_i*lengthscale%i_j)' % (i, i, i) for i in range(input_dim)]) + lengthscales = [sp.var('lengthscale%i_i lengthscale%i_j' % i, positive=True) for i in range(real_input_dim)] + dist_string = ' + '.join(['(x_%i-z_%i)**2/(lengthscale%i_i*lengthscale%i_j)' % (i, i, i) for i in range(real_input_dim)]) dist = parse_expr(dist_string) f = variance*sp.exp(-dist/2.) else: lengthscale = sp.var('lengthscale_i lengthscale_j',positive=True) - dist_string = ' + '.join(['(x_%i-z_%i)**2' % (i, i) for i in range(input_dim)]) + dist_string = ' + '.join(['(x_%i-z_%i)**2' % (i, i) for i in range(real_input_dim)]) dist = parse_expr(dist_string) - f = variance*sp.exp(-dist/(2*lengthscale**2)) - return kern(input_dim, [spkern(input_dim, f, name='eq_sympy')]) + f = variance*sp.exp(-dist/(2*lengthscale_i*lengthscale_j)) + return kern(input_dim, [spkern(input_dim, f, output_dim=output_dim, name='eq_sympy')]) def sinc(input_dim, ARD=False, variance=1., lengthscale=1.): """ diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py index ff7dd1c1..08f36109 100644 --- a/GPy/kern/kern.py +++ b/GPy/kern/kern.py @@ -658,7 +658,7 @@ class Kern_check_dKdiag_dX(Kern_check_model): def _set_params(self, x): self.X=x.reshape(self.X.shape) -def kern_test(kern, X=None, X2=None, verbose=False): +def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False): """This function runs on kernels to check the correctness of their implementation. It checks that the covariance function is positive definite for a randomly generated data set. :param kern: the kernel to be tested. @@ -672,12 +672,12 @@ def kern_test(kern, X=None, X2=None, verbose=False): pass_checks = True if X==None: X = np.random.randn(10, kern.input_dim) - for ind in kern.output_indicator: - X[:, ind] = np.random.randint(kern.output_dim, X.shape[0]) + if output_ind is not None: + X[:, output_ind] = np.random.randint(kern.output_dim, X.shape[0]) if X2==None: X2 = np.random.randn(20, kern.input_dim) - for ind in kern.output_indicator: - X2[:, ind] = np.random.randint(kern.output_dim, X2.shape[0]) + if output_ind is not None: + X2[:, output_ind] = np.random.randint(kern.output_dim, X2.shape[0]) if verbose: print("Checking covariance function is positive definite.") diff --git a/GPy/kern/parts/kernpart.py b/GPy/kern/parts/kernpart.py index 95deeb81..f6777083 100644 --- a/GPy/kern/parts/kernpart.py +++ b/GPy/kern/parts/kernpart.py @@ -12,8 +12,6 @@ class Kernpart(object): Do not instantiate. """ - # stores indices of any inputs that are for indicating outputs - self.output_indicator = [] # the input dimensionality for the covariance self.input_dim = input_dim # the number of optimisable parameters diff --git a/GPy/kern/parts/sympykern.py b/GPy/kern/parts/sympykern.py index a9f73436..09ab9934 100644 --- a/GPy/kern/parts/sympykern.py +++ b/GPy/kern/parts/sympykern.py @@ -44,7 +44,6 @@ class spkern(Kernpart): assert len(self._sp_x)==len(self._sp_z) self.input_dim = len(self._sp_x) if output_dim > 1: - self.output_indicator=[self.input_dim] self.input_dim += 1 assert self.input_dim == input_dim @@ -84,7 +83,7 @@ class spkern(Kernpart): if param is not None: if param.has_key(theta): val = param[theta] - setattr(self, theta, val) + setattr(self, theta.name, val) #deal with param self._set_params(self._get_params()) @@ -146,7 +145,7 @@ class spkern(Kernpart): reverse_arg_list = list(arg_list) reverse_arg_list.reverse() - param_arg_list = ["param[%i]"%i for i in range(self.num_shared_params)] + param_arg_list = [shared_params.name for shared_params in self._sp_theta] arg_list += param_arg_list precompute_list=[] @@ -201,11 +200,12 @@ class spkern(Kernpart): """%(diag_precompute_string,diag_arg_string,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed # Code to compute gradients - func_list = ([' '*16 + 'target[%i] += partial[i*num_inducing+j]*dk_d%s(%s);'%(i,theta.name,arg_string) for i,theta in enumerate(self._sp_theta)]) + func_list = [] if self.output_dim>1: func_list += [' '*16 + "int %s=(int)%s[%s*input_dim+output_dim];"%(index, var, index2) for index, var, index2 in zip(['ii', 'jj'], ['X', 'Z'], ['i', 'j'])] func_list += [' '*16 + 'target[%i+ii] += partial[i*num_inducing+j]*dk_d%s(%s);'%(self.num_shared_params+i*self.output_dim, theta.name, arg_string) for i, theta in enumerate(self._sp_theta_i)] func_list += [' '*16 + 'target[%i+jj] += partial[i*num_inducing+j]*dk_d%s(%s);'%(self.num_shared_params+i*self.output_dim, theta.name, reverse_arg_string) for i, theta in enumerate(self._sp_theta_i)] + func_list += ([' '*16 + 'target[%i] += partial[i*num_inducing+j]*dk_d%s(%s);'%(i,theta.name,arg_string) for i,theta in enumerate(self._sp_theta)]) func_string = '\n'.join(func_list) self._dK_dtheta_code =\ @@ -290,7 +290,9 @@ class spkern(Kernpart): #TODO: insert multiple functions here via string manipulation #TODO: similar functions for psi_stats def _get_arg_names(self, Z=None, partial=None): - arg_names = ['target','X','param'] + arg_names = ['target','X'] + for shared_params in self._sp_theta: + arg_names += [shared_params.name] if Z is not None: arg_names += ['Z'] if partial is not None: @@ -301,7 +303,9 @@ class spkern(Kernpart): return arg_names def _weave_inline(self, code, X, target, Z=None, partial=None): - param, output_dim = self._shared_params, self.output_dim + output_dim = self.output_dim + for shared_params in self._sp_theta: + locals()[shared_params.name] = getattr(self, shared_params.name) # Need to extract parameters first for split_params in self._split_theta_names: @@ -369,9 +373,7 @@ class spkern(Kernpart): def _set_params(self,param): assert param.size == (self.num_params) for i, shared_params in enumerate(self._sp_theta): - start = i - end = i+1 - setattr(self, shared_params, param[start:end]) + setattr(self, shared_params.name, param[i]) if self.output_dim>1: for i, split_params in enumerate(self._split_theta_names): @@ -383,7 +385,7 @@ class spkern(Kernpart): def _get_params(self): params = np.zeros(0) for shared_params in self._sp_theta: - params = np.hstack((params, getattr(self, shared_params))) + params = np.hstack((params, getattr(self, shared_params.name))) if self.output_dim>1: for split_params in self._split_theta_names: params = np.hstack((params, getattr(self, split_params).flatten())) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 5c45ae20..f64dac2b 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -34,7 +34,7 @@ class KernelTests(unittest.TestCase): self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) def test_eq_sympykernel(self): - kern = GPy.kern.eq_sympy(5, 3) + kern = GPy.kern.eq_sympy(5, 3, output_ind=4) self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose)) def test_sinckernel(self):