Added first draft of functionality for multiple output sympy kernels.

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.

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)]
+        X = sp.symbols('x_:' + str(input_dim))
-        Z = [sp.var('z%i' % i) for i in range(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)]
+        X = sp.symbols('x_:' + str(input_dim))
-        Z = [sp.var('z%i' % i) for i in range(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):

GPy/kern/parts/sympy_helpers.cpp

@@ -1,4 +1,7 @@
 #include <math.h>
+#include <float.h>
+#include <stdlib.h>
+
 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_MAX<xmax ? DBL_MAX : xmax;
+  // Find values where erfcx can be evaluated
+  double t = 3.97886080735226 / (abs(x) + 3.97886080735226);
+  double u = t-0.5;
+  double y = (((((((((u * 0.00127109764952614092 + 1.19314022838340944e-4) * u 
+	      - 0.003963850973605135)   * u - 8.70779635317295828e-4) * u 
+	    + 0.00773672528313526668) * u + 0.00383335126264887303) * u 
+	  - 0.0127223813782122755)  * u - 0.0133823644533460069)  * u 
+	+ 0.0161315329733252248)  * u + 0.0390976845588484035)  * u + 0.00249367200053503304;
+  if (x<xneg)
+    return -INFINITY;
+  else if (x<0)
+    return 2*exp(x*x)-y;
+  else if (x>xmax)
+    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;
+}

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);

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 multiple inputs, call them x_1, z_1, etc
-     - to handle multpile correlated outputs, you'll need to define each covariance function and 'cross' variance function. TODO
+     - 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_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]=='z'],key=lambda z:int(z.name[1:]))
+        self._sp_z= sorted([e for e in sp_vars if e.name[0:2]=='z_'],key=lambda z:int(z.name[2:]))
-        assert all([x.name=='x%i'%i for i,x in enumerate(self._sp_x)])
+        # Check that variable names make sense.
-        assert all([z.name=='z%i'%i for i,z in enumerate(self._sp_z)])
+        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.output_dim = output_dim
-        self.num_params = len(self._sp_theta)
+        # 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)]
+        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)]
+        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!
@@ -101,12 +139,18 @@ class spkern(Kernpart):
 
     def _gen_code(self):
         #generate c functions from sympy objects        
-        (foo_c,self._function_code),(foo_h,self._function_header) = \
+        argument_sequence = self._sp_x+self._sp_z+self._sp_theta
-                codegen([('k',self._sp_k)] \
+        code_list = [('k',self._sp_k)]
-                + [('dk_d%s'%x.name,dx) for x,dx in zip(self._sp_x,self._sp_dk_dx)]\
+        # gradients with respect to covariance input
-                #+ [('dk_d%s'%z.name,dz) for z,dz in zip(self._sp_z,self._sp_dk_dz)]\
+        code_list += [('dk_d%s'%x.name,dx) for x,dx in zip(self._sp_x,self._sp_dk_dx)]
-                + [('dk_d%s'%theta.name,dtheta) for theta,dtheta in zip(self._sp_theta,self._sp_dk_dtheta)]\
+        # gradient with respect to parameters
-                ,"C",'foobar',argument_sequence=self._sp_x+self._sp_z+self._sp_theta)
+        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)
 
+        # 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
-        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)])
-
-        
         self._K_code =\
         """
         int i;
@@ -131,19 +191,19 @@ class spkern(Kernpart):
         //#pragma omp parallel for private(j)
         for (i=0;i<N;i++){
             for (j=0;j<num_inducing;j++){
+%s
                 target[i*num_inducing+j] = k(%s);
             }
         }
         %s
-        """%(arglist,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed
+        """%(precompute_string,arg_string,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed
-
-        # Similar code when only X is provided. 
-        self._K_code_X = self._K_code.replace('Z[', 'X[')
 
         
         # Code to compute diagonal of covariance.
-        diag_arglist = re.sub('Z','X',arglist)
+        diag_arg_string = re.sub('Z','X',arg_string)
-        diag_arglist = re.sub('j','i',diag_arglist)
+        diag_arg_string = re.sub('j','i',diag_arg_string)
+        diag_precompute_string = re.sub('Z','X',precompute_string)
+        diag_precompute_string = re.sub('j','i',diag_precompute_string)
         # Code to do the looping for Kdiag
         self._Kdiag_code =\
         """
@@ -152,13 +212,19 @@ class spkern(Kernpart):
         int input_dim = X_array->dimensions[1];
         //#pragma omp parallel for
         for (i=0;i<N;i++){
+                %s
                 target[i] = k(%s);
         }
         %s
-        """%(diag_arglist,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed
+        """%(diag_precompute_string,diag_arg_string,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed
 
         # Code to compute gradients
-        funclist = '\n'.join([' '*16 + 'target[%i] += partial[i*num_inducing+j]*dk_d%s(%s);'%(i,theta.name,arglist) for i,theta in  enumerate(self._sp_theta)])
+        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)])
+        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_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_func_string = re.sub('Z','X',func_string,count=0)
-        diag_funclist = re.sub('j','i',diag_funclist)
+        diag_func_string = re.sub('j','i',diag_func_string)
-        diag_funclist = re.sub('partial\[i\*num_inducing\+i\]','partial[i]',diag_funclist)
+        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)])
+        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)])
-        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);}"""
 
         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
+
+        # 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):        
-        param = self._param
         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
+        self._weave_inline(self._Kdiag_code, X, target)
-        weave.inline(self._Kdiag_code,arg_names=['target','X','param'],**self.weave_kwargs)
 
     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
+        self._weave_inline(self._dKdiag_dtheta_code, X, target, Z=None, partial=partial)
-        weave.inline(self._dKdiag_dtheta_code,arg_names=['target','X','param','partial'],**self.weave_kwargs)
                
     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
+        self._weave.inline(self._dKdiag_dX_code, X, target, Z, partial)
-        weave.inline(self._dKdiag_dX_code,arg_names=['target','X','param','partial'],**self.weave_kwargs)
 
     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]

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:
+        if x.is_Number:
-            return S.Zero
+            if x is S.NaN:
-        else:
+                return S.NaN
-            return (x*cos(x) - sin(x))/(x*x)
+            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):
+    def eval(cls, arg):
-        if x is S.Zero:
+        if arg.is_Number:
-            return S.One
+            if arg is S.NaN:
-        else:
+                return S.NaN
-            return sin(x)/x
+            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
+