From 966fe4934541a43476984efa46b1207215d45d8a Mon Sep 17 00:00:00 2001 From: Neil Lawrence Date: Tue, 8 Oct 2013 08:25:26 +0100 Subject: [PATCH] Added first draft of functionality for multiple output sympy kernels. --- GPy/inference/scg.py | 2 +- GPy/kern/constructors.py | 20 +-- GPy/kern/parts/sympy_helpers.cpp | 36 +++++ GPy/kern/parts/sympy_helpers.h | 3 + GPy/kern/parts/sympykern.py | 226 ++++++++++++++++++++++--------- GPy/util/symbolic.py | 85 ++++++++++-- 6 files changed, 281 insertions(+), 91 deletions(-) diff --git a/GPy/inference/scg.py b/GPy/inference/scg.py index f4c7c9c4..252f348e 100644 --- a/GPy/inference/scg.py +++ b/GPy/inference/scg.py @@ -62,7 +62,7 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=np.inf, display=True, fnow = fold gradnew = gradf(x, *optargs) # Initial gradient. if any(np.isnan(gradnew)): - raise UnexpectedInfOrNan + raise UnexpectedInfOrNan, "Gradient contribution resulted in a NaN value" current_grad = np.dot(gradnew, gradnew) gradold = gradnew.copy() d = -gradnew # Initial search direction. diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index a8ec1d4b..e6952186 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -298,17 +298,17 @@ if sympy_available: """ Radial Basis Function covariance. """ - X = [sp.var('x%i' % i) for i in range(input_dim)] - Z = [sp.var('z%i' % i) for i in range(input_dim)] + X = sp.symbols('x_:' + str(input_dim)) + Z = sp.symbols('z_:' + str(input_dim)) variance = sp.var('variance',positive=True) if ARD: lengthscales = [sp.var('lengthscale_%i' % i, positive=True) for i in range(input_dim)] - dist_string = ' + '.join(['(x%i-z%i)**2/lengthscale_%i**2' % (i, i, i) for i in range(input_dim)]) + dist_string = ' + '.join(['(x_%i-z_%i)**2/lengthscale_%i**2' % (i, i, i) for i in range(input_dim)]) dist = parse_expr(dist_string) f = variance*sp.exp(-dist/2.) else: lengthscale = sp.var('lengthscale',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(input_dim)]) dist = parse_expr(dist_string) f = variance*sp.exp(-dist/(2*lengthscale**2)) return kern(input_dim, [spkern(input_dim, f, name='rbf_sympy')]) @@ -318,23 +318,23 @@ if sympy_available: TODO: Not clear why this isn't working, suggests argument of sinc is not a number. sinc covariance funciton """ - X = [sp.var('x%i' % i) for i in range(input_dim)] - Z = [sp.var('z%i' % i) for i in range(input_dim)] + X = sp.symbols('x_:' + str(input_dim)) + Z = sp.symbols('z_:' + str(input_dim)) variance = sp.var('variance',positive=True) if ARD: lengthscales = [sp.var('lengthscale_%i' % i, positive=True) for i in range(input_dim)] - dist_string = ' + '.join(['(x%i-z%i)**2/lengthscale_%i**2' % (i, i, i) for i in range(input_dim)]) + dist_string = ' + '.join(['(x_%i-z_%i)**2/lengthscale_%i**2' % (i, i, i) for i in range(input_dim)]) dist = parse_expr(dist_string) f = variance*sinc(sp.pi*sp.sqrt(dist)) else: lengthscale = sp.var('lengthscale',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(input_dim)]) dist = parse_expr(dist_string) f = variance*sinc(sp.pi*sp.sqrt(dist)/lengthscale) return kern(input_dim, [spkern(input_dim, f, name='sinc')]) - def sympykern(input_dim, k,name=None): + def sympykern(input_dim, k=None, output_dim=1, name=None, param=None): """ A base kernel object, where all the hard work in done by sympy. @@ -349,7 +349,7 @@ if sympy_available: - to handle multiple inputs, call them x1, z1, etc - to handle multpile correlated outputs, you'll need to define each covariance function and 'cross' variance function. TODO """ - return kern(input_dim, [spkern(input_dim, k,name)]) + return kern(input_dim, [spkern(input_dim, k=k, output_dim=output_dim, name=name, param=param)]) del sympy_available def periodic_exponential(input_dim=1, variance=1., lengthscale=None, period=2 * np.pi, n_freq=10, lower=0., upper=4 * np.pi): diff --git a/GPy/kern/parts/sympy_helpers.cpp b/GPy/kern/parts/sympy_helpers.cpp index 76dba4eb..e4df4d80 100644 --- a/GPy/kern/parts/sympy_helpers.cpp +++ b/GPy/kern/parts/sympy_helpers.cpp @@ -1,4 +1,7 @@ #include +#include +#include + double DiracDelta(double x){ // TODO: this doesn't seem to be a dirac delta ... should return infinity. Neil if((x<0.000001) & (x>-0.000001))//go on, laugh at my c++ skills @@ -23,3 +26,36 @@ double sinc_grad(double x){ else return (x*cos(x) - sin(x))/(x*x); } + +double erfcx(double x){ + double xneg=-sqrt(log(DBL_MAX/2)); + double xmax = 1/(sqrt(M_PI)*DBL_MIN); + xmax = DBL_MAXxmax) + return 0.0; + else + return y; +} + +double ln_diff_erf(double x0, double x1){ + if (x0==x1) + return INFINITY; + else if(x0<0 && x1>0 || x0>0 && x1<0) + return log(erf(x0)-erf(x1)); + else if(x1>0) + return log(erfcx(x1)-erfcx(x0)*exp(x1*x1)- x0*x0)-x1*x1; + else + return log(erfcx(-x0)-erfcx(-x1)*exp(x0*x0 - x1*x1))-x0*x0; +} diff --git a/GPy/kern/parts/sympy_helpers.h b/GPy/kern/parts/sympy_helpers.h index d5b495ca..56220167 100644 --- a/GPy/kern/parts/sympy_helpers.h +++ b/GPy/kern/parts/sympy_helpers.h @@ -4,3 +4,6 @@ double DiracDelta(double x, int foo); double sinc(double x); double sinc_grad(double x); + +double erfcx(double x); +double ln_diff_erf(double x0, double x1); diff --git a/GPy/kern/parts/sympykern.py b/GPy/kern/parts/sympykern.py index 9755e37b..dc6a5390 100644 --- a/GPy/kern/parts/sympykern.py +++ b/GPy/kern/parts/sympykern.py @@ -9,6 +9,7 @@ import sys current_dir = os.path.dirname(os.path.abspath(os.path.dirname(__file__))) import tempfile import pdb +import ast from kernpart import Kernpart class spkern(Kernpart): @@ -16,41 +17,78 @@ class spkern(Kernpart): A kernel object, where all the hard work in done by sympy. :param k: the covariance function - :type k: a positive definite sympy function of x1, z1, x2, z2... + :type k: a positive definite sympy function of x_0, z_0, x_1, z_1, x_2, z_2... To construct a new sympy kernel, you'll need to define: - a kernel function using a sympy object. Ensure that the kernel is of the form k(x,z). - that's it! we'll extract the variables from the function k. Note: - - to handle multiple inputs, call them x1, z1, etc - - to handle multpile correlated outputs, you'll need to define each covariance function and 'cross' variance function. TODO + - to handle multiple inputs, call them x_1, z_1, etc + - to handle multpile correlated outputs, you'll need to add parameters with an index, such as lengthscale_i and lengthscale_j. """ - def __init__(self,input_dim,k,name=None,param=None): + def __init__(self,input_dim, k=None, output_dim=1, name=None, param=None): if name is None: self.name='sympykern' else: self.name = name + if k is None: + raise ValueError, "You must provide an argument for the covariance function." self._sp_k = k sp_vars = [e for e in k.atoms() if e.is_Symbol] - self._sp_x= sorted([e for e in sp_vars if e.name[0]=='x'],key=lambda x:int(x.name[1:])) - self._sp_z= sorted([e for e in sp_vars if e.name[0]=='z'],key=lambda z:int(z.name[1:])) - assert all([x.name=='x%i'%i for i,x in enumerate(self._sp_x)]) - assert all([z.name=='z%i'%i for i,z in enumerate(self._sp_z)]) + self._sp_x= sorted([e for e in sp_vars if e.name[0:2]=='x_'],key=lambda x:int(x.name[2:])) + self._sp_z= sorted([e for e in sp_vars if e.name[0:2]=='z_'],key=lambda z:int(z.name[2:])) + # Check that variable names make sense. + assert all([x.name=='x_%i'%i for i,x in enumerate(self._sp_x)]) + assert all([z.name=='z_%i'%i for i,z in enumerate(self._sp_z)]) assert len(self._sp_x)==len(self._sp_z) self.input_dim = len(self._sp_x) + if output_dim > 1: + self.input_dim += 1 assert self.input_dim == input_dim - self._sp_theta = sorted([e for e in sp_vars if not (e.name[0]=='x' or e.name[0]=='z')],key=lambda e:e.name) - self.num_params = len(self._sp_theta) + self.output_dim = output_dim + # extract parameter names + thetas = sorted([e for e in sp_vars if not (e.name[0:2]=='x_' or e.name[0:2]=='z_')],key=lambda e:e.name) + + + # Look for parameters with index. + if self.output_dim>1: + self._sp_theta_i = sorted([e for e in thetas if (e.name[-2:]=='_i')], key=lambda e:e.name) + self._sp_theta_j = sorted([e for e in thetas if (e.name[-2:]=='_j')], key=lambda e:e.name) + # Make sure parameter appears with both indices! + assert len(self._sp_theta_i)==len(self._sp_theta_j) + assert all([theta_i.name[:-2]==theta_j.name[:-2] for theta_i, theta_j in zip(self._sp_theta_i, self._sp_theta_j)]) + + # Extract names of shared parameters + self._sp_theta = [theta for theta in thetas if theta not in self._sp_theta_i and theta not in self._sp_theta_j] + + self.num_split_params = len(self._sp_theta_i) + self._split_param_names = ["%s"%theta.name[:-2] for theta in self._sp_theta_i] + for params in self._split_param_names: + setattr(self, params, 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 + + else: + self.num_split_params = 0 + self._split_param_names = [] + self._sp_theta = thetas + self.num_shared_params = len(self._sp_theta) + self.num_params = self.num_shared_params #deal with param if param is None: param = np.ones(self.num_params) + assert param.size==self.num_params self._set_params(param) #Differentiate! self._sp_dk_dtheta = [sp.diff(k,theta).simplify() for theta in self._sp_theta] + if self.output_dim > 1: + self._sp_dk_dtheta_i = [sp.diff(k,theta).simplify() for theta in self._sp_theta_i] + self._sp_dk_dx = [sp.diff(k,xi).simplify() for xi in self._sp_x] #self._sp_dk_dz = [sp.diff(k,zi) for zi in self._sp_z] @@ -72,8 +110,8 @@ class spkern(Kernpart): def compute_psi_stats(self): #define some normal distributions - mus = [sp.var('mu%i'%i,real=True) for i in range(self.input_dim)] - Ss = [sp.var('S%i'%i,positive=True) for i in range(self.input_dim)] + mus = [sp.var('mu_%i'%i,real=True) for i in range(self.input_dim)] + Ss = [sp.var('S_%i'%i,positive=True) for i in range(self.input_dim)] normals = [(2*sp.pi*Si)**(-0.5)*sp.exp(-0.5*(xi-mui)**2/Si) for xi, mui, Si in zip(self._sp_x, mus, Ss)] #do some integration! @@ -100,13 +138,19 @@ class spkern(Kernpart): def _gen_code(self): - #generate c functions from sympy objects - (foo_c,self._function_code),(foo_h,self._function_header) = \ - codegen([('k',self._sp_k)] \ - + [('dk_d%s'%x.name,dx) for x,dx in zip(self._sp_x,self._sp_dk_dx)]\ - #+ [('dk_d%s'%z.name,dz) for z,dz in zip(self._sp_z,self._sp_dk_dz)]\ - + [('dk_d%s'%theta.name,dtheta) for theta,dtheta in zip(self._sp_theta,self._sp_dk_dtheta)]\ - ,"C",'foobar',argument_sequence=self._sp_x+self._sp_z+self._sp_theta) + #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 + code_list += [('dk_d%s'%x.name,dx) for x,dx in zip(self._sp_x,self._sp_dk_dx)] + # gradient with respect to parameters + code_list += [('dk_d%s'%theta.name,dtheta) for theta,dtheta in zip(self._sp_theta,self._sp_dk_dtheta)] + # gradient with respect to multiple output parameters + 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)] + (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) @@ -115,12 +159,28 @@ class spkern(Kernpart): # Substitute any known derivatives which sympy doesn't compute self._function_code = re.sub('DiracDelta\(.+?,.+?\)','0.0',self._function_code) - # Here's the code to do the looping for K - arglist = ", ".join(["X[i*input_dim+%s]"%x.name[1:] for x in self._sp_x] - + ["Z[j*input_dim+%s]"%z.name[1:] for z in self._sp_z] - + ["param[%i]"%i for i in range(self.num_params)]) + # 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]) + if self.output_dim>1: + reverse_arg_list = list(arg_list) + reverse_arg_list.reverse() - + param_arg_list = ["param[%i]"%i for i in range(self.num_shared_params)] + arg_list += param_arg_list + + precompute_list=[] + if self.output_dim > 1: + reverse_arg_list+=list(param_arg_list) + split_param_arg_list = ["%s[%s]"%(theta.name[:-2],index) for index in ['ii', 'jj'] for theta in self._sp_theta_i] + split_param_reverse_arg_list = ["%s[%s]"%(theta.name[:-2],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 + precompute_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'])] + 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 =\ """ int i; @@ -131,19 +191,19 @@ class spkern(Kernpart): //#pragma omp parallel for private(j) for (i=0;idimensions[1]; //#pragma omp parallel for for (i=0;i1: + 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_string = '\n'.join(func_list) self._dK_dtheta_code =\ """ @@ -174,15 +240,13 @@ class spkern(Kernpart): } } %s - """%(funclist,"/*"+str(self._sp_k)+"*/") # adding a string representation forces recompile when needed + """%(func_string,"/*"+str(self._sp_k)+"*/") # adding a string representation forces recompile when needed - # Similar code when only X is provided, change argument lists. - self._dK_dtheta_code_X = self._dK_dtheta_code.replace('Z[', 'X[') # Code to compute gradients for Kdiag TODO: needs clean up - diag_funclist = re.sub('Z','X',funclist,count=0) - diag_funclist = re.sub('j','i',diag_funclist) - diag_funclist = re.sub('partial\[i\*num_inducing\+i\]','partial[i]',diag_funclist) + diag_func_string = re.sub('Z','X',func_string,count=0) + diag_func_string = re.sub('j','i',diag_func_string) + diag_func_string = re.sub('partial\[i\*num_inducing\+i\]','partial[i]',diag_func_string) self._dKdiag_dtheta_code =\ """ int i; @@ -192,13 +256,10 @@ class spkern(Kernpart): %s } %s - """%(diag_funclist,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed + """%(diag_func_string,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed # Code for gradients wrt X - gradient_funcs = "\n".join(["target[i*input_dim+%i] += partial[i*num_inducing+j]*dk_dx%i(%s);"%(q,q,arglist) for q in range(self.input_dim)]) - if False: - gradient_funcs += """if(isnan(target[i*input_dim+2])){printf("%%f\\n",dk_dx2(X[i*input_dim+0], X[i*input_dim+1], X[i*input_dim+2], Z[j*input_dim+0], Z[j*input_dim+1], Z[j*input_dim+2], param[0], param[1], param[2], param[3], param[4], param[5]));} - if(isnan(target[i*input_dim+2])){printf("%%f,%%f,%%i,%%i\\n", X[i*input_dim+2], Z[j*input_dim+2],i,j);}""" + gradient_funcs = "\n".join(["target[i*input_dim+%i] += partial[i*num_inducing+j]*dk_dx%i(%s);"%(q,q,arg_string) for q in range(self.input_dim)]) self._dK_dX_code = \ """ @@ -216,8 +277,6 @@ class spkern(Kernpart): %s """%(gradient_funcs,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed - # Create code for call when just X is passed as argument. - self._dK_dX_code_X = self._dK_dX_code.replace('Z[', 'X[').replace('+= partial[', '+= 2*partial[') diag_gradient_funcs = re.sub('Z','X',gradient_funcs,count=0) diag_gradient_funcs = re.sub('j','i',diag_gradient_funcs) @@ -235,52 +294,85 @@ class spkern(Kernpart): """%(diag_gradient_funcs,"/*"+str(self._sp_k)+"*/") #adding a # string representation forces recompile when needed Get rid # of Zs in argument for diagonal. TODO: Why wasn't - # diag_funclist called here? Need to check that. + # diag_func_string called here? Need to check that. #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[') + #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'] + if Z is not None: + arg_names += ['Z'] + if partial is not None: + arg_names += ['partial'] + if self.output_dim>1: + arg_names += self._split_param_names + arg_names += ['output_dim'] + return arg_names + + def _weave_inline(self, code, X, target, Z=None, partial=None): + param, output_dim = self._shared_params, self.output_dim - def K(self,X,Z,target): - param = self._param + # Need to extract parameters first + for split_params in self._split_param_names: + locals()[split_params] = getattr(self, split_params) + arg_names = self._get_arg_names(Z, partial) + weave.inline(code=code, arg_names=arg_names,**self.weave_kwargs) + + def K(self,X,Z,target): if Z is None: - weave.inline(self._K_code_X,arg_names=['target','X','param'],**self.weave_kwargs) + self._weave_inline(self._K_code_X, X, target) else: - weave.inline(self._K_code,arg_names=['target','X','Z','param'],**self.weave_kwargs) + self._weave_inline(self._K_code, X, target, Z) + def Kdiag(self,X,target): - param = self._param - weave.inline(self._Kdiag_code,arg_names=['target','X','param'],**self.weave_kwargs) + self._weave_inline(self._Kdiag_code, X, target) def dK_dtheta(self,partial,X,Z,target): - param = self._param if Z is None: - weave.inline(self._dK_dtheta_code_X, arg_names=['target','X','param','partial'],**self.weave_kwargs) + self._weave_inline(self._dK_dtheta_code_X, X, target, Z, partial) else: - weave.inline(self._dK_dtheta_code, arg_names=['target','X','Z','param','partial'],**self.weave_kwargs) - + self._weave_inline(self._dK_dtheta_code, X, target, Z, partial) + def dKdiag_dtheta(self,partial,X,target): - param = self._param - weave.inline(self._dKdiag_dtheta_code,arg_names=['target','X','param','partial'],**self.weave_kwargs) - + self._weave_inline(self._dKdiag_dtheta_code, X, target, Z=None, partial=partial) + def dK_dX(self,partial,X,Z,target): - param = self._param if Z is None: - weave.inline(self._dK_dX_code_X,arg_names=['target','X','param','partial'],**self.weave_kwargs) + self._weave_inline(self._dK_dX_code_X, X, target, Z, partial) else: - weave.inline(self._dK_dX_code,arg_names=['target','X','Z','param','partial'],**self.weave_kwargs) + self._weave_inline(self._dK_dX_code, X, target, Z, partial) def dKdiag_dX(self,partial,X,target): - param = self._param - weave.inline(self._dKdiag_dX_code,arg_names=['target','X','param','partial'],**self.weave_kwargs) + self._weave.inline(self._dKdiag_dX_code, X, target, Z, partial) def _set_params(self,param): #print param.flags['C_CONTIGUOUS'] - self._param = param.copy() + assert param.size == (self.num_params) + self._shared_params = param[0:self.num_shared_params] + if self.output_dim>1: + for i, split_params in enumerate(self._split_param_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): - return self._param + params = self._shared_params + if self.output_dim>1: + for split_params in self._split_param_names: + params = np.hstack((params, getattr(self, split_params).flatten())) + return params def _get_param_names(self): - return [x.name for x in self._sp_theta] + 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] diff --git a/GPy/util/symbolic.py b/GPy/util/symbolic.py index f4f5fda0..8b368a77 100644 --- a/GPy/util/symbolic.py +++ b/GPy/util/symbolic.py @@ -1,32 +1,91 @@ -from sympy import Function, S, oo, I, cos, sin +from sympy import Function, S, oo, I, cos, sin, asin, log, erf,pi,exp +class ln_diff_erf(Function): + nargs = 2 + + def fdiff(self, argindex=2): + if argindex == 2: + x0, x1 = self.args + return -2*exp(-x1**2)/(sqrt(pi)*(erf(x0)-erf(x1))) + elif argindex == 1: + x0, x1 = self.args + return 2*exp(-x0**2)/(sqrt(pi)*(erf(x0)-erf(x1))) + else: + raise ArgumentIndexError(self, argindex) + + @classmethod + def eval(cls, x0, x1): + if x0.is_Number and x1.is_Number: + return log(erf(x0)-erf(x1)) + +class sim_h(Function): + nargs = 5 + + @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))) + +class erfc(Function): + nargs = 1 + + @classmethod + def eval(cls, arg): + return 1-erf(arg) + +class erfcx(Function): + nargs = 1 + + @classmethod + def eval(cls, arg): + return erfc(arg)*exp(arg*arg) + class sinc_grad(Function): nargs = 1 def fdiff(self, argindex=1): - return ((2-x*x)*sin(self.args[0]) - 2*x*cos(x))/(x*x*x) + if argindex==1: + # Strictly speaking this should be computed separately, as it won't work when x=0. See http://calculus.subwiki.org/wiki/Sinc_function + return ((2-x*x)*sin(self.args[0]) - 2*x*cos(x))/(x*x*x) + else: + raise ArgumentIndexError(self, argindex) + @classmethod def eval(cls, x): - if x is S.Zero: - return S.Zero - else: - return (x*cos(x) - sin(x))/(x*x) + if x.is_Number: + if x is S.NaN: + return S.NaN + elif x is S.Zero: + return S.Zero + else: + return (x*cos(x) - sin(x))/(x*x) class sinc(Function): nargs = 1 def fdiff(self, argindex=1): - return sinc_grad(self.args[0]) + if argindex==1: + return sinc_grad(self.args[0]) + else: + raise ArgumentIndexError(self, argindex) + @classmethod - def eval(cls, x): - if x is S.Zero: - return S.One - else: - return sin(x)/x - + def eval(cls, arg): + if arg.is_Number: + if arg is S.NaN: + return S.NaN + elif arg is S.Zero: + return S.One + else: + return sin(arg)/arg + + if arg.func is asin: + x = arg.args[0] + return x / arg + def _eval_is_real(self): return self.args[0].is_real +