From da2a88826d670f4284d466dd291d539b9428cf47 Mon Sep 17 00:00:00 2001 From: Neil Lawrence Date: Mon, 14 Oct 2013 22:09:41 +0100 Subject: [PATCH] Basic sim code functional. --- GPy/core/model.py | 2 +- GPy/kern/constructors.py | 4 +-- GPy/kern/parts/sympykern.py | 67 ++++++++++++++++++++++++++----------- GPy/util/symbolic.py | 12 ++++++- 4 files changed, 62 insertions(+), 23 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 7aff8f4d..c1ab7b6a 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -259,7 +259,7 @@ class Model(Parameterized): these terms are present in the name the parameter is constrained positive. """ - positive_strings = ['variance', 'lengthscale', 'precision', 'kappa'] + positive_strings = ['variance', 'lengthscale', 'precision', 'decay', 'kappa'] # param_names = self._get_param_names() currently_constrained = self.all_constrained_indices() to_make_positive = [] diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index c6a6672f..392f43ba 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -330,11 +330,11 @@ if sympy_available: dist = parse_expr(dist_string) f = variance*sp.exp(-dist/2.) else: - lengthscale = sp.var('lengthscale_i lengthscale_j',positive=True) + lengthscales = sp.var('lengthscale_i lengthscale_j',positive=True) shared_lengthscale = sp.var('shared_lengthscale',positive=True) dist_string = ' + '.join(['(x_%i-z_%i)**2' % (i, i) for i in range(real_input_dim)]) dist = parse_expr(dist_string) - f = scale_i*scale_j*sp.exp(-dist/(2*(shared_lengthscale**2 + lengthscale_i*lengthscale_j))) + f = scale_i*scale_j*sp.exp(-dist/(2*(lengthscale_i**2 + lengthscale_j**2 + shared_lengthscale**2))) 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/parts/sympykern.py b/GPy/kern/parts/sympykern.py index ea603eab..88c179aa 100644 --- a/GPy/kern/parts/sympykern.py +++ b/GPy/kern/parts/sympykern.py @@ -117,6 +117,9 @@ class spkern(Kernpart): return spkern(self._sp_k+other._sp_k) def _gen_code(self): + """Generates the C functions necessary for computing the covariance function using the sympy objects as input.""" + #TODO: maybe generate one C function only to save compile time? Also easier to take that as a basis and hand craft other covariances?? + #generate c functions from sympy objects argument_sequence = self._sp_x+self._sp_z+self._sp_theta code_list = [('k',self._sp_k)] @@ -138,15 +141,20 @@ class spkern(Kernpart): # Substitute any known derivatives which sympy doesn't compute self._function_code = re.sub('DiracDelta\(.+?,.+?\)','0.0',self._function_code) - # This is the basic argument construction for the C code. - #arg_list = (["X[i*input_dim+%s]"%x.name[2:] for x in self._sp_x] - # + ["Z[j*input_dim+%s]"%z.name[2:] for z in self._sp_z]) + + ############################################################ + # This is the basic argument construction for the C code. # + ############################################################ + arg_list = (["X2(i, %s)"%x.name[2:] for x in self._sp_x] + ["Z2(j, %s)"%z.name[2:] for z in self._sp_z]) + + # for multiple outputs need to also provide these arguments reversed. if self.output_dim>1: reverse_arg_list = list(arg_list) reverse_arg_list.reverse() + # Add in any 'shared' parameters to the list. param_arg_list = [shared_params.name for shared_params in self._sp_theta] arg_list += param_arg_list @@ -163,6 +171,15 @@ class spkern(Kernpart): reverse_arg_string = ", ".join(reverse_arg_list) arg_string = ", ".join(arg_list) precompute_string = "\n".join(precompute_list) + + # Code to compute argments string needed when only X is provided. + X_arg_string = re.sub('Z','X',arg_string) + # Code to compute argument string when only diagonal is required. + diag_arg_string = re.sub('int jj','//int jj',X_arg_string) + diag_arg_string = re.sub('j','i',diag_arg_string) + diag_precompute_string = precompute_list[0] + + # Here's the code to do the looping for K self._K_code =\ """ @@ -184,14 +201,28 @@ class spkern(Kernpart): %s """%(precompute_string,arg_string,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed - - # Code to compute diagonal of covariance. - diag_arg_string = re.sub('Z','X',arg_string) - diag_arg_string = re.sub('int jj','//int jj',diag_arg_string) - diag_arg_string = re.sub('j','i',diag_arg_string) - diag_precompute_string = re.sub('int jj','//int jj',precompute_string) - diag_precompute_string = re.sub('Z','X',diag_precompute_string) - diag_precompute_string = re.sub('j','i',diag_precompute_string) + self._K_code_X = """ + // _K_code_X + // Code for computing the covariance function. + int i; + int j; + int N = target_array->dimensions[0]; + int num_inducing = target_array->dimensions[1]; + int input_dim = X_array->dimensions[1]; + //#pragma omp parallel for private(j) + for (i=0;i1: grad_func_list += c_define_output_indices - grad_func_list += [' '*16 + 'TARGET1(%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)] - grad_func_list += [' '*16 + 'TARGET1(%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)] - grad_func_list += ([' '*16 + 'TARGET1(%i) += partial[i*num_inducing+j]*dk_d%s(%s);'%(i,theta.name,arg_string) for i,theta in enumerate(self._sp_theta)]) + grad_func_list += [' '*16 + 'TARGET1(%i+ii) += PARTIAL2(i, 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)] + grad_func_list += [' '*16 + 'TARGET1(%i+jj) += PARTIAL2(i, 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)] + grad_func_list += ([' '*16 + 'TARGET1(%i) += PARTIAL2(i, j)*dk_d%s(%s);'%(i,theta.name,arg_string) for i,theta in enumerate(self._sp_theta)]) grad_func_string = '\n'.join(grad_func_list) self._dK_dtheta_code =\ @@ -241,7 +272,7 @@ class spkern(Kernpart): diag_grad_func_string = re.sub('Z','X',grad_func_string,count=0) diag_grad_func_string = re.sub('int jj','//int jj',diag_grad_func_string) diag_grad_func_string = re.sub('j','i',diag_grad_func_string) - diag_grad_func_string = re.sub('partial\[i\*num_inducing\+i\]','partial[i]',diag_grad_func_string) + diag_grad_func_string = re.sub('PARTIAL2\(i, i\)','PARTIAL1(i)',diag_grad_func_string) self._dKdiag_dtheta_code =\ """ // _dKdiag_dtheta_code @@ -259,7 +290,7 @@ class spkern(Kernpart): gradX_func_list = [] if self.output_dim>1: gradX_func_list += c_define_output_indices - gradX_func_list += ["TARGET2(i, %i) += partial[i*num_inducing+j]*dk_dx_%i(%s);"%(q,q,arg_string) for q in range(self._real_input_dim)] + gradX_func_list += ["TARGET2(i, %i) += PARTIAL2(i, j)*dk_dx_%i(%s);"%(q,q,arg_string) for q in range(self._real_input_dim)] gradX_func_string = "\n".join(gradX_func_list) self._dK_dX_code = \ @@ -284,7 +315,7 @@ class spkern(Kernpart): diag_gradX_func_string = re.sub('Z','X',gradX_func_string,count=0) diag_gradX_func_string = re.sub('int jj','//int jj',diag_gradX_func_string) diag_gradX_func_string = re.sub('j','i',diag_gradX_func_string) - diag_gradX_func_string = re.sub('partial\[i\*num_inducing\+i\]','2*partial[i]',diag_gradX_func_string) + diag_gradX_func_string = re.sub('PARTIAL2\(i, i\)','2*PARTIAL1(i)',diag_gradX_func_string) # Code for gradients of Kdiag wrt X self._dKdiag_dX_code= \ @@ -304,10 +335,8 @@ class spkern(Kernpart): #self._dKdiag_dX_code = self._dKdiag_dX_code.replace('Z[j', 'X[i') # Code to use when only X is provided. - self._K_code_X = self._K_code.replace('Z[', 'X[') self._dK_dtheta_code_X = self._dK_dtheta_code.replace('Z[', 'X[') self._dK_dX_code_X = self._dK_dX_code.replace('Z[', 'X[').replace('+= partial[', '+= 2*partial[') - self._K_code_X = self._K_code.replace('Z2(', 'X2(') self._dK_dtheta_code_X = self._dK_dtheta_code.replace('Z2(', 'X2(') self._dK_dX_code_X = self._dK_dX_code.replace('Z2(', 'X2(') diff --git a/GPy/util/symbolic.py b/GPy/util/symbolic.py index 8b368a77..10c59a5e 100644 --- a/GPy/util/symbolic.py +++ b/GPy/util/symbolic.py @@ -22,9 +22,19 @@ class ln_diff_erf(Function): class sim_h(Function): nargs = 5 + def fdiff(self, argindex=1): + pass + @classmethod def eval(cls, t, tprime, d_i, d_j, l): - return exp((d_j/2*l)**2)/(d_i+d_j)*(exp(-d_j*(tprime - t))*(erf((tprime-t)/l - d_j/2*l) + erf(t/l + d_j/2*l)) - exp(-(d_j*tprime + d_i))*(erf(tprime/l - d_j/2*l) + erf(d_j/2*l))) + # putting in the is_Number stuff forces it to look for a fdiff method for derivative. + return (exp((d_j/2*l)**2)/(d_i+d_j) + *(exp(-d_j*(tprime - t)) + *(erf((tprime-t)/l - d_j/2*l) + + erf(t/l + d_j/2*l)) + - exp(-(d_j*tprime + d_i)) + *(erf(tprime/l - d_j/2*l) + + erf(d_j/2*l)))) class erfc(Function): nargs = 1