diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index 44f4ae3f..083960b4 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -600,3 +600,20 @@ def ODE_1(input_dim=1, varianceU=1., varianceY=1., lengthscaleU=None, lengthsc """ part = parts.ODE_1.ODE_1(input_dim, varianceU, varianceY, lengthscaleU, lengthscaleY) return kern(input_dim, [part]) + +def ODE_UY(input_dim=2, varianceU=1., varianceY=1., lengthscaleU=None, lengthscaleY=None): + """ + kernel resultiong from a first order ODE with OU driving GP + :param input_dim: the number of input dimension, has to be equal to one + :type input_dim: int + :param input_lengthU: the number of input U length + :param varianceU: variance of the driving GP + :type varianceU: float + :param varianceY: 'variance' of the transfer function + :type varianceY: float + :param lengthscaleY: 'lengthscale' of the transfer function + :type lengthscaleY: float + :rtype: kernel object + """ + part = parts.ODE_UY.ODE_UY(input_dim, varianceU, varianceY, lengthscaleU, lengthscaleY) + return kern(input_dim, [part]) diff --git a/GPy/kern/parts/__init__.py b/GPy/kern/parts/__init__.py index 0a758f1e..f278941a 100644 --- a/GPy/kern/parts/__init__.py +++ b/GPy/kern/parts/__init__.py @@ -14,6 +14,7 @@ import Matern32 import Matern52 import mlp import ODE_1 +#import ODE_UY import periodic_exponential import periodic_Matern32 import periodic_Matern52 @@ -26,4 +27,5 @@ import rbf import rbf_inv import spline import symmetric +import sympy_helpers import white diff --git a/GPy/kern/parts/sympy_helpers.cpp b/GPy/kern/parts/sympy_helpers.cpp index d21d2683..9f30eea9 100644 --- a/GPy/kern/parts/sympy_helpers.cpp +++ b/GPy/kern/parts/sympy_helpers.cpp @@ -1,3 +1,4 @@ +#include "Python.h" #include #include #include @@ -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_MAXxmax) 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; diff --git a/GPy/kern/parts/sympy_helpers.py b/GPy/kern/parts/sympy_helpers.py new file mode 100644 index 00000000..125dac58 --- /dev/null +++ b/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