apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-06-02 14:45:15 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

weird merge

Browse source
This commit is contained in:
James Hensman 2013-11-14 08:58:26 +00:00
commit 2387861e25
3 changed files with 299 additions and 17 deletions
Download patch file Download diff file Expand all files Collapse all files

106
GPy/kern/parts/sympy_helpers.cpp
View file

@ -1,7 +1,8 @@
#include <math.h>
#include <float.h>
#include <stdlib.h>
#include <iostream>
#include <stdexcept>
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_MAX<xmax ? DBL_MAX : xmax;
@ -50,12 +54,108 @@ 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)
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;
}

7
GPy/kern/parts/sympy_helpers.h
View file

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

203
GPy/util/symbolic.py
View file

@ -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