diff --git a/GPy/kern/parts/sympykern.py b/GPy/kern/parts/sympykern.py index 2d015b27..a5bb7b1d 100644 --- a/GPy/kern/parts/sympykern.py +++ b/GPy/kern/parts/sympykern.py @@ -1,17 +1,31 @@ -import numpy as np -import sympy as sp -from sympy.utilities.codegen import codegen -from sympy.core.cache import clear_cache +try: + import sympy as sp + sympy_available=True +except ImportError: + sympy_available=False + exit() + +from sympy.core.cache import clear_cache +from sympy.utilities.codegen import codegen + +try: + from scipy import weave + weave_available = True +except ImportError: + weave_available = False -from scipy import weave -import re import os -import sys current_dir = os.path.dirname(os.path.abspath(os.path.dirname(__file__))) +import sys +import numpy as np +import re import tempfile import pdb import ast + from kernpart import Kernpart +from ...core.parameterization import Param +from ...core.parameterization.transformations import Logexp class spkern(Kernpart): """ @@ -75,17 +89,20 @@ class spkern(Kernpart): self.num_split_params = len(self._sp_theta_i) self._split_theta_names = ["%s"%theta.name[:-2] for theta in self._sp_theta_i] for theta in self._split_theta_names: - setattr(self, theta, np.ones(self.output_dim)) + setattr(self, theta, Param(theta, np.ones(self.output_dim), None)) + self.add_parameters(getattr(self, theta)) + + #setattr(self, theta, np.ones(self.output_dim)) self.num_shared_params = len(self._sp_theta) - self.num_params = self.num_shared_params+self.num_split_params*self.output_dim + #self.num_params = self.num_shared_params+self.num_split_params*self.output_dim else: self.num_split_params = 0 self._split_theta_names = [] self._sp_theta = thetas self.num_shared_params = len(self._sp_theta) - self.num_params = self.num_shared_params + #self.num_params = self.num_shared_params # Add parameters to the model. for theta in self._sp_theta: @@ -93,9 +110,12 @@ class spkern(Kernpart): if param is not None: if param.has_key(theta): val = param[theta] - setattr(self, theta.name, val) + #setattr(self, theta.name, val) + setattr(self, theta.name, Param(theta.name, val, None)) + self.add_parameters(getattr(self, theta.name)) + self.parameters_changed() # initializes cache #deal with param - self._set_params(self._get_params()) + #self._set_params(self._get_params()) # Differentiate with respect to parameters. self._sp_dk_dtheta = [sp.diff(k,theta).simplify() for theta in self._sp_theta] @@ -112,26 +132,26 @@ class spkern(Kernpart): # generate the code for the covariance functions self._gen_code() - if weave - if False: - extra_compile_args = ['-ftree-vectorize', '-mssse3', '-ftree-vectorizer-verbose=5'] - else: - extra_compile_args = [] + if weave_available: + if False: + extra_compile_args = ['-ftree-vectorize', '-mssse3', '-ftree-vectorizer-verbose=5'] + else: + extra_compile_args = [] - self.weave_kwargs = { - 'support_code':self._function_code, - 'include_dirs':[tempfile.gettempdir(), os.path.join(current_dir,'parts/')], - 'headers':['"sympy_helpers.h"'], - 'sources':[os.path.join(current_dir,"parts/sympy_helpers.cpp")], - 'extra_compile_args':extra_compile_args, - 'extra_link_args':['-lgomp'], - 'verbose':True} + self.weave_kwargs = { + 'support_code':self._function_code, + 'include_dirs':[tempfile.gettempdir(), os.path.join(current_dir,'parts/')], + 'headers':['"sympy_helpers.h"'], + 'sources':[os.path.join(current_dir,"parts/sympy_helpers.cpp")], + 'extra_compile_args':extra_compile_args, + 'extra_link_args':['-lgomp'], + 'verbose':True} def __add__(self,other): return spkern(self._sp_k+other._sp_k) def _gen_code(self): - #generate c functions from sympy objects + argument_sequence = self._sp_x+self._sp_z+self._sp_theta code_list = [('k',self._sp_k)] # gradients with respect to covariance input @@ -142,193 +162,224 @@ class spkern(Kernpart): if self.output_dim > 1: argument_sequence += self._sp_theta_i + self._sp_theta_j code_list += [('dk_d%s'%theta.name,dtheta) for theta,dtheta in zip(self._sp_theta_i,self._sp_dk_dtheta_i)] + # generate c functions from sympy objects + if weave_available: + code_type = "C" + else: + code_type = "PYTHON" (foo_c,self._function_code), (foo_h,self._function_header) = \ - codegen(code_list, "C",'foobar',argument_sequence=argument_sequence) - #put the header file where we can find it - f = file(os.path.join(tempfile.gettempdir(),'foobar.h'),'w') - f.write(self._function_header) - f.close() + codegen(code_list, + code_type, + self.name, + argument_sequence=argument_sequence) + + # Use weave to compute the underlying functions. + if weave_available: + # put the header file where we can find it + f = file(os.path.join(tempfile.gettempdir(), self.name + '.h'),'w') + f.write(self._function_header) + f.close() + + # 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]) - 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]) - if self.output_dim>1: - reverse_arg_list = list(arg_list) - reverse_arg_list.reverse() + if weave_available: + # arg_list will store the arguments required 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]) - param_arg_list = [shared_params.name for shared_params in self._sp_theta] - arg_list += param_arg_list + # for multiple outputs reverse argument list is also required + if self.output_dim>1: + reverse_arg_list = list(arg_list) + reverse_arg_list.reverse() - precompute_list=[] - if self.output_dim > 1: - reverse_arg_list+=list(param_arg_list) - split_param_arg_list = ["%s1(%s)"%(theta.name[:-2].upper(),index) for index in ['ii', 'jj'] for theta in self._sp_theta_i] - split_param_reverse_arg_list = ["%s1(%s)"%(theta.name[:-2].upper(),index) for index in ['jj', 'ii'] for theta in self._sp_theta_i] - arg_list += split_param_arg_list - reverse_arg_list += split_param_reverse_arg_list - # Extract the right output indices from the inputs. - c_define_output_indices = [' '*16 + "int %s=(int)%s(%s, %i);"%(index, var, index2, self.input_dim-1) for index, var, index2 in zip(['ii', 'jj'], ['X2', 'Z2'], ['i', 'j'])] - precompute_list += c_define_output_indices - reverse_arg_string = ", ".join(reverse_arg_list) - arg_string = ", ".join(arg_list) - precompute_string = "\n".join(precompute_list) - # Here's the code to do the looping for K - self._K_code =\ - """ - // _K_code - // 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;i 1: + reverse_arg_list+=list(param_arg_list) + # For multiple outputs, also need the split parameters. + split_param_arg_list = ["%s1(%s)"%(theta.name[:-2].upper(),index) for index in ['ii', 'jj'] for theta in self._sp_theta_i] + split_param_reverse_arg_list = ["%s1(%s)"%(theta.name[:-2].upper(),index) for index in ['jj', 'ii'] for theta in self._sp_theta_i] + arg_list += split_param_arg_list + reverse_arg_list += split_param_reverse_arg_list + # Extract the right output indices from the inputs. + c_define_output_indices = [' '*16 + "int %s=(int)%s(%s, %i);"%(index, var, index2, self.input_dim-1) for index, var, index2 in zip(['ii', 'jj'], ['X2', 'Z2'], ['i', 'j'])] + precompute_list += c_define_output_indices + reverse_arg_string = ", ".join(reverse_arg_list) + arg_string = ", ".join(arg_list) + precompute_string = "\n".join(precompute_list) + + # Now we use the arguments in code that computes the separate parts. + + # Any precomputations will be done here eventually. + self._precompute = \ + """ + // Precompute code would go here. It will be called when parameters are updated. + """ + + # Here's the code to do the looping for K + self._K_code =\ + """ + // _K_code + // 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;idimensions[0]; - int input_dim = X_array->dimensions[1]; - //#pragma omp parallel for - 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_string = '\n'.join(grad_func_list) - - self._dK_dtheta_code =\ - """ - // _dK_dtheta_code - // Code for computing gradient of covariance with respect to parameters. - int i; - int j; - int N = partial_array->dimensions[0]; - int num_inducing = partial_array->dimensions[1]; - int input_dim = X_array->dimensions[1]; - //#pragma omp parallel for private(j) - for (i=0;idimensions[0]; + int input_dim = X_array->dimensions[1]; + //#pragma omp parallel for + 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_string = '\n'.join(grad_func_list) + + self._dK_dtheta_code =\ + """ + // _dK_dtheta_code + // Code for computing gradient of covariance with respect to parameters. + int i; + int j; + int N = partial_array->dimensions[0]; + int num_inducing = partial_array->dimensions[1]; + int input_dim = X_array->dimensions[1]; + //#pragma omp parallel for private(j) + for (i=0;idimensions[0]; - int input_dim = X_array->dimensions[1]; - for (i=0;idimensions[0]; + int input_dim = X_array->dimensions[1]; + for (i=0;i1: + 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_string = "\n".join(gradX_func_list) + + self._dK_dX_code = \ + """ + // _dK_dX_code + // Code for computing gradient of covariance with respect to inputs. + int i; + int j; + int N = partial_array->dimensions[0]; + int num_inducing = partial_array->dimensions[1]; + int input_dim = X_array->dimensions[1]; + //#pragma omp parallel for private(j) + for (i=0;i1: - 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_string = "\n".join(gradX_func_list) - - self._dK_dX_code = \ - """ - // _dK_dX_code - // Code for computing gradient of covariance with respect to inputs. - int i; - int j; - int N = partial_array->dimensions[0]; - int num_inducing = partial_array->dimensions[1]; - int input_dim = X_array->dimensions[1]; - //#pragma omp parallel for private(j) - for (i=0;idimensions[0]; - int input_dim = X_array->dimensions[1]; - for (int i=0;idimensions[0]; + int input_dim = X_array->dimensions[1]; + for (int i=0;i1: - for i, split_params in enumerate(self._split_theta_names): - start = self.num_shared_params + i*self.output_dim - end = self.num_shared_params + (i+1)*self.output_dim - setattr(self, split_params, param[start:end]) + # if self.output_dim>1: + # for i, split_params in enumerate(self._split_theta_names): + # start = self.num_shared_params + i*self.output_dim + # end = self.num_shared_params + (i+1)*self.output_dim + # setattr(self, split_params, param[start:end]) - def _get_params(self): - params = np.zeros(0) - for shared_params in self._sp_theta: - 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())) - return params + # def _get_params(self): + # params = np.zeros(0) + # for shared_params in self._sp_theta: + # 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())) + # return params - def _get_param_names(self): - if self.output_dim>1: - return [x.name for x in self._sp_theta] + [x.name[:-2] + str(i) for x in self._sp_theta_i for i in range(self.output_dim)] - else: - return [x.name for x in self._sp_theta] + # def _get_param_names(self): + # if self.output_dim>1: + # return [x.name for x in self._sp_theta] + [x.name[:-2] + str(i) for x in self._sp_theta_i for i in range(self.output_dim)] + # else: + # return [x.name for x in self._sp_theta]