From e79794f6ef7f337dc3a500f13fe864b79293e8c5 Mon Sep 17 00:00:00 2001 From: Neil Lawrence Date: Thu, 14 Nov 2013 08:47:16 +0000 Subject: [PATCH] Part implementation of ode_eq functionality. Not yet numerically stable or efficient (some horrible use of cut and paste to get things working ...) --- GPy/kern/parts/sympy_helpers.cpp | 106 +++++++++++++++- GPy/kern/parts/sympy_helpers.h | 7 ++ GPy/util/symbolic.py | 203 ++++++++++++++++++++++++++++--- 3 files changed, 299 insertions(+), 17 deletions(-) diff --git a/GPy/kern/parts/sympy_helpers.cpp b/GPy/kern/parts/sympy_helpers.cpp index e4df4d80..d21d2683 100644 --- a/GPy/kern/parts/sympy_helpers.cpp +++ b/GPy/kern/parts/sympy_helpers.cpp @@ -1,7 +1,8 @@ #include #include #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 @@ -14,6 +15,7 @@ double DiracDelta(double x,int foo){ }; double sinc(double x){ + // compute the sinc function if (x==0) return 1.0; else @@ -21,6 +23,7 @@ double sinc(double x){ } double sinc_grad(double x){ + // compute the gradient of the sinc function. if (x==0) return 0.0; else @@ -28,6 +31,7 @@ double sinc_grad(double x){ } double erfcx(double x){ + // compute the scaled complex error function. double xneg=-sqrt(log(DBL_MAX/2)); double xmax = 1/(sqrt(M_PI)*DBL_MIN); xmax = DBL_MAXx0) + 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) - return INFINITY; + 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; + 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; } + +double h(double t, double tprime, double d_i, double d_j, double l){ + // Compute the h function for the sim covariance. + 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 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, 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))); +} + +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; + double hv = h(t, tprime, d_i, d_j, l); + 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_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 = ((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)))); + 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; + return base/(d_i+d_j); +} + + +double dh_dt(double t, double tprime, double d_i, double d_j, double l){ + return 0.0; +} + +double dh_dtprime(double t, double tprime, double d_i, double d_j, double l){ + return 0.0; +} diff --git a/GPy/kern/parts/sympy_helpers.h b/GPy/kern/parts/sympy_helpers.h index 56220167..5e58d5d2 100644 --- a/GPy/kern/parts/sympy_helpers.h +++ b/GPy/kern/parts/sympy_helpers.h @@ -7,3 +7,10 @@ double sinc_grad(double x); double erfcx(double x); double ln_diff_erf(double x0, double x1); + +double h(double t, double tprime, double d_i, double d_j, double l); +double dh_dl(double t, double tprime, double d_i, double d_j, double l); +double dh_dd_i(double t, double tprime, double d_i, double d_j, double l); +double dh_dd_j(double t, double tprime, double d_i, double d_j, double l); +double dh_dt(double t, double tprime, double d_i, double d_j, double l); +double dh_dtprime(double t, double tprime, double d_i, double d_j, double l); diff --git a/GPy/util/symbolic.py b/GPy/util/symbolic.py index 0b5ca381..d546f940 100644 --- a/GPy/util/symbolic.py +++ b/GPy/util/symbolic.py @@ -1,4 +1,4 @@ -from sympy import Function, S, oo, I, cos, sin, asin, log, erf,pi,exp +from sympy import Function, S, oo, I, cos, sin, asin, log, erf,pi,exp,sqrt,sign class ln_diff_erf(Function): @@ -19,15 +19,84 @@ class ln_diff_erf(Function): if x0.is_Number and x1.is_Number: return log(erf(x0)-erf(x1)) -class sim_h(Function): +class dh_dd_i(Function): + nargs = 5 + @classmethod + def eval(cls, t, tprime, d_i, d_j, l): + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + + diff_t = (t-tprime) + l2 = l*l + h = h(t, tprime, d_i, d_j, l) + half_l_di = 0.5*l*d_i + arg_1 = half_l_di + tprime/l + arg_2 = half_l_di - (t-tprime)/l + ln_part_1 = ln_diff_erf(arg_1, arg_2) + arg_1 = half_l_di + arg_2 = half_l_di - t/l + sign_val = sign(t/l) + ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l) + + base = ((0.5*d_i*l2*(d_i+d_j)-1)*h + + (-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(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) + +class dh_dd_j(Function): + nargs = 5 + @classmethod + def eval(cls, t, tprime, d_i, d_j, l): + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + diff_t = (t-tprime) + l2 = l*l + half_l_di = 0.5*l*d_i + h = h(t, tprime, d_i, d_j, l) + arg_1 = half_l_di + tprime/l + arg_2 = half_l_di - (t-tprime)/l + ln_part_1 = ln_diff_erf(arg_1, arg_2) + arg_1 = half_l_di + arg_2 = half_l_di - t/l + sign_val = sign(t/l) + ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l) + sign_val = sign(t/l) + base = tprime*sign_val*exp(half_l_di*half_l_di-(d_i*t+d_j*tprime)+ln_part_2)-h + return base/(d_i+d_j) + +class dh_dl(Function): + nargs = 5 + @classmethod + def eval(cls, t, tprime, d_i, d_j, l): + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + + diff_t = (t-tprime) + l2 = l*l + h = h(t, tprime, d_i, d_j, l) + return 0.5*d_i*d_i*l*h + 2./(sqrt(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))) + +class dh_dt(Function): nargs = 5 - - def fdiff(self, argindex=1): - pass - @classmethod def eval(cls, t, tprime, d_i, d_j, l): - # putting in the is_Number stuff forces it to look for a fdiff method for derivative. if (t.is_Number and tprime.is_Number and d_i.is_Number @@ -40,13 +109,119 @@ class sim_h(Function): or l is S.NaN): return S.NaN else: - 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)))) + half_l_di = 0.5*l*d_i + arg_1 = half_l_di + tprime/l + arg_2 = half_l_di - (t-tprime)/l + ln_part_1 = ln_diff_erf(arg_1, arg_2) + arg_1 = half_l_di + arg_2 = half_l_di - t/l + sign_val = sign(t/l) + ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l) + + + 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))).diff(t) + +class dh_dtprime(Function): + nargs = 5 + @classmethod + def eval(cls, t, tprime, d_i, d_j, l): + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + if (t is S.NaN + or tprime is S.NaN + or d_i is S.NaN + or d_j is S.NaN + or l is S.NaN): + return S.NaN + else: + half_l_di = 0.5*l*d_i + arg_1 = half_l_di + tprime/l + arg_2 = half_l_di - (t-tprime)/l + ln_part_1 = ln_diff_erf(arg_1, arg_2) + arg_1 = half_l_di + arg_2 = half_l_di - t/l + sign_val = sign(t/l) + ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l) + + + 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))).diff(tprime) + + +class h(Function): + nargs = 5 + def fdiff(self, argindex=5): + t, tprime, d_i, d_j, l = self.args + if argindex == 1: + return dh_dt(t, tprime, d_i, d_j, l) + elif argindex == 2: + return dh_dtprime(t, tprime, d_i, d_j, l) + elif argindex == 3: + return dh_dd_i(t, tprime, d_i, d_j, l) + elif argindex == 4: + return dh_dd_j(t, tprime, d_i, d_j, l) + elif argindex == 5: + return dh_dl(t, tprime, d_i, d_j, l) + + + @classmethod + def eval(cls, t, tprime, d_i, d_j, l): + # putting in the is_Number stuff forces it to look for a fdiff method for derivative. If it's left out, then when asking for self.diff, it just does the diff on the eval symbolic terms directly. We want to avoid that because we are looking to ensure everything is numerically stable. Maybe it's because of the if statement that this happens? + if (t.is_Number + and tprime.is_Number + and d_i.is_Number + and d_j.is_Number + and l.is_Number): + if (t is S.NaN + or tprime is S.NaN + or d_i is S.NaN + or d_j is S.NaN + or l is S.NaN): + return S.NaN + else: + half_l_di = 0.5*l*d_i + arg_1 = half_l_di + tprime/l + arg_2 = half_l_di - (t-tprime)/l + ln_part_1 = ln_diff_erf(arg_1, arg_2) + arg_1 = half_l_di + arg_2 = half_l_di - t/l + sign_val = sign(t/l) + ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l) + + + 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))) + + + # 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