Working eq_ode1 in sympy now.

GPy/kern/parts/__init__.py

@@ -26,4 +26,5 @@ import rbf
 import rbf_inv
 import spline
 import symmetric
+import sympy_helpers
 import white

GPy/kern/parts/sympy_helpers.cpp

@@ -1,3 +1,4 @@
+#include "Python.h"
 #include <math.h>
 #include <float.h>
 #include <stdlib.h>
@@ -29,24 +30,33 @@ double sinc_grad(double x){
   else 
     return (x*cos(x) - sin(x))/(x*x);
 }
-
 double erfcx(double x){
+  // Based on code by Soren Hauberg 2010 for Octave.
   // compute the scaled complex error function.
+  //return erfc(x)*exp(x*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 t = 3.97886080735226 / (fabs(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;
+		     - 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;
+  y = ((((((((((((y * u - 0.0838864557023001992) * u -		       
+		 0.119463959964325415) * u + 0.0166207924969367356) * u + 
+	       0.357524274449531043) * u + 0.805276408752910567)  * u + 
+	     1.18902982909273333)  * u + 1.37040217682338167)   * u +	
+	   1.31314653831023098)  * u + 1.07925515155856677)   * u +	
+	 0.774368199119538609) * u + 0.490165080585318424)  * u +	
+       0.275374741597376782) * t;
+
   if (x<xneg)
     return -INFINITY;
   else if (x<0)
-    return 2*exp(x*x)-y;
+    return 2.0*exp(x*x)-y;
   else if (x>xmax)
     return 0.0;
   else 
@@ -55,16 +65,19 @@ double erfcx(double x){
 
 double ln_diff_erf(double x0, double x1){
   // stably compute the log of difference between two erfs.
-  if (x1>x0)
-    throw std::runtime_error("Error: second argument must be smaller than first in ln_diff_err");
-  return log(erf(x0) - erf(x1));
-  if (x0==x1)
+  if (x1>x0){
+    PyErr_SetString(PyExc_RuntimeError,"second argument must be smaller than or equal to first in ln_diff_erf");
+    throw 1;
+  }
+  if (x0==x1){
+    PyErr_WarnEx(PyExc_RuntimeWarning,"divide by zero encountered in log", 1);
     return -INFINITY;
-  else if(x0<0 && x1>0 || x0>0 && x1<0)
+  }
+  else if(x0<0 && x1>0 || x0>0 && x1<0) //x0 and x1 have opposite signs
     return log(erf(x0)-erf(x1));
-  else if(x1>0)
-    return log(erfcx(x1)-erfcx(x0)*exp(x1*x1- x0*x0))-x1*x1;
-  else 
+  else if(x0>0) //x0 positive, x1 non-negative
+    return log(erfcx(x1)-erfcx(x0)*exp(x1*x1- x0*x0))-x1*x1; 
+  else //x0 and x1 non-positive
     return log(erfcx(-x0)-erfcx(-x1)*exp(x0*x0 - x1*x1))-x0*x0;
 }
 
@@ -80,26 +93,19 @@ double h(double t, double tprime, double d_i, double d_j, double l){
     sign_val = 0.0;
   else if (t/l < 0)
     sign_val = -1.0;
-  double ln_part_2 = ln_diff_erf(half_l_di, arg_2);
-  
-  return sign_val*exp(half_l_di*half_l_di - d_i*(t-tprime) + ln_part_1 - log(d_i + d_j)) - sign_val*exp(half_l_di*half_l_di - d_i*t - d_j*tprime + ln_part_2 - log(d_i + d_j));
-}
-
-double dh_dl(double t, double tprime, double d_i, double d_j, double l){
-  // compute gradient of h function with respect to lengthscale for sim covariance
-  // TODO a lot of energy wasted recomputing things here, need to do this in a shared way somehow ... perhaps needs rewrite of sympykern.
-  double half_l_di = 0.5*l*d_i;
-  double arg_1 = half_l_di + tprime/l;
-  double arg_2 = half_l_di - (t-tprime)/l;
-  double ln_part_1 = ln_diff_erf(arg_1, arg_2);
   arg_2 = half_l_di - t/l;
   double ln_part_2 = ln_diff_erf(half_l_di, arg_2);
-  double diff_t = t - tprime;
-  double l2 = l*l;
-  double hv = h(t, tprime, d_i, d_j, l);
-  return 0.5*d_i*d_i*l*hv + 2/(sqrt(M_PI)*(d_i+d_j))*((-diff_t/l2-d_i/2)*exp(-diff_t*diff_t/l2)+(-tprime/l2+d_i/2)*exp(-tprime*tprime/l2-d_i*t)-(-t/l2-d_i/2)*exp(-t*t/l2-d_j*tprime)-d_i/2*exp(-(d_i*t+d_j*tprime)));
+  // if either ln_part_1 or ln_part_2 are -inf, don't bother computing rest of that term.
+  double part_1 = 0.0;
+  if(isfinite(ln_part_1))
+    part_1 = sign_val*exp(half_l_di*half_l_di - d_i*(t-tprime) + ln_part_1 - log(d_i + d_j));
+  double part_2 = 0.0;
+  if(isfinite(ln_part_2))
+    part_2 = sign_val*exp(half_l_di*half_l_di - d_i*t - d_j*tprime + ln_part_2 - log(d_i + d_j));
+  return part_1 - part_2;
 }
 
+
 double dh_dd_i(double t, double tprime, double d_i, double d_j, double l){
   double diff_t = (t-tprime);
   double l2 = l*l;
@@ -116,41 +122,52 @@ double dh_dd_i(double t, double tprime, double d_i, double d_j, double l){
   else if (t/l < 0)
     sign_val = -1.0;
   double ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l);
-
-  double base = ((0.5*d_i*l2*(d_i+d_j)-1)*hv 
-		 + (-diff_t*sign_val*exp(half_l_di*half_l_di
-					 -d_i*diff_t
-					 +ln_part_1)
-		    +t*sign_val*exp(half_l_di*half_l_di
-				    -d_i*t-d_j*tprime
-				    +ln_part_2))
-		 + l/sqrt(M_PI)*(-exp(-diff_t*diff_t/l2)
-			       +exp(-tprime*tprime/l2-d_i*t)
-			       +exp(-t*t/l2-d_j*tprime)
-			       -exp(-(d_i*t + d_j*tprime))));
+  double base = (0.5*d_i*l2*(d_i+d_j)-1)*hv;
+  if(isfinite(ln_part_1))
+    base -= diff_t*sign_val*exp(half_l_di*half_l_di
+				-d_i*diff_t
+				+ln_part_1);
+  if(isfinite(ln_part_2))
+    base += t*sign_val*exp(half_l_di*half_l_di
+			   -d_i*t-d_j*tprime
+			   +ln_part_2);
+  base += l/sqrt(M_PI)*(-exp(-diff_t*diff_t/l2)
+			+exp(-tprime*tprime/l2-d_i*t)
+			+exp(-t*t/l2-d_j*tprime)
+			-exp(-(d_i*t + d_j*tprime)));
   return base/(d_i+d_j);
+
 }
 
 double dh_dd_j(double t, double tprime, double d_i, double d_j, double l){
-  double diff_t = (t-tprime);
-  double l2 = l*l;
   double half_l_di = 0.5*l*d_i;
   double hv = h(t, tprime, d_i, d_j, l);
-  double arg_1 = half_l_di + tprime/l;
-  double arg_2 = half_l_di - (t-tprime)/l;
-  double ln_part_1 = ln_diff_erf(arg_1, arg_2);
-  arg_1 = half_l_di;
-  arg_2 = half_l_di - t/l;
   double sign_val = 1.0;
   if(t/l==0)
     sign_val = 0.0;
   else if (t/l < 0)
     sign_val = -1.0;
   double ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l);
-  double base = tprime*sign_val*exp(half_l_di*half_l_di-(d_i*t+d_j*tprime)+ln_part_2)-hv;
+  double base = -hv;
+  if(isfinite(ln_part_2))
+    base += tprime*sign_val*exp(half_l_di*half_l_di-(d_i*t+d_j*tprime)+ln_part_2);
   return base/(d_i+d_j);
 }
 
+double dh_dl(double t, double tprime, double d_i, double d_j, double l){
+  // compute gradient of h function with respect to lengthscale for sim covariance
+  // TODO a lot of energy wasted recomputing things here, need to do this in a shared way somehow ... perhaps needs rewrite of sympykern.
+  double half_l_di = 0.5*l*d_i;
+  double arg_1 = half_l_di + tprime/l;
+  double arg_2 = half_l_di - (t-tprime)/l;
+  double ln_part_1 = ln_diff_erf(arg_1, arg_2);
+  arg_2 = half_l_di - t/l;
+  double ln_part_2 = ln_diff_erf(half_l_di, arg_2);
+  double diff_t = t - tprime;
+  double l2 = l*l;
+  double hv = h(t, tprime, d_i, d_j, l);
+  return 0.5*d_i*d_i*l*hv + 2/(sqrt(M_PI)*(d_i+d_j))*((-diff_t/l2-d_i/2)*exp(-diff_t*diff_t/l2)+(-tprime/l2+d_i/2)*exp(-tprime*tprime/l2-d_i*t)-(-t/l2-d_i/2)*exp(-t*t/l2-d_j*tprime)-d_i/2*exp(-(d_i*t+d_j*tprime)));
+}
 
 double dh_dt(double t, double tprime, double d_i, double d_j, double l){
   return 0.0;

GPy/kern/parts/sympy_helpers.py

@@ -0,0 +1,71 @@
+# Code for testing functions written in sympy_helpers.cpp
+from scipy import weave
+import tempfile
+import os
+import numpy as np
+current_dir = os.path.dirname(os.path.abspath(os.path.dirname(__file__)))
+extra_compile_args = []
+
+weave_kwargs = {
+    'support_code': "",
+    'include_dirs':[tempfile.gettempdir(), current_dir],
+    'headers':['"parts/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 erfcx(x):
+    code = """
+        // Code for computing scaled complementary erf
+        int i;
+        int dim;
+        int elements = Ntarget[0];
+        for (dim=1; dim<Dtarget; dim++)
+          elements *= Ntarget[dim];
+        for (i=0;i<elements;i++) 
+            target[i] = erfcx(x[i]);
+        """
+    x = np.asarray(x)
+    arg_names = ['target','x']
+    target = np.zeros_like(x)
+    weave.inline(code=code, arg_names=arg_names,**weave_kwargs)
+    return target
+
+def ln_diff_erf(x, y):
+    code = """
+        // Code for computing scaled complementary erf
+        int i;
+        int dim;
+        int elements = Ntarget[0];
+        for (dim=1; dim<Dtarget; dim++)
+          elements *= Ntarget[dim];
+        for (i=0;i<elements;i++) 
+          target[i] = ln_diff_erf(x[i], y[i]);
+        """
+    x = np.asarray(x)
+    y = np.asarray(y)
+    assert(x.shape==y.shape)
+    target = np.zeros_like(x)
+    arg_names = ['target','x', 'y']
+    weave.inline(code=code, arg_names=arg_names,**weave_kwargs)
+    return target
+
+def h(t, tprime, d_i, d_j, l):
+    code = """
+        // Code for computing the 1st order ODE h helper function.
+        int i;
+        int dim;
+        int elements = Ntarget[0];
+        for (dim=1; dim<Dtarget; dim++)
+          elements *= Ntarget[dim];
+        for (i=0;i<elements;i++) 
+          target[i] = h(t[i], tprime[i], d_i, d_j, l);
+        """
+    t = np.asarray(t)
+    tprime = np.asarray(tprime)
+    assert(tprime.shape==t.shape)
+    target = np.zeros_like(t)
+    arg_names = ['target','t', 'tprime', 'd_i', 'd_j', 'l']
+    weave.inline(code=code, arg_names=arg_names,**weave_kwargs)
+    return target

GPy/util/symbolic.py

@@ -10,7 +10,7 @@ class ln_diff_erf(Function):
             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)))
+            return 2.*exp(-x0**2)/(sqrt(pi)*(erf(x0)-erf(x1)))
         else:
             raise ArgumentIndexError(self, argindex)