diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py index cb968520..548e2924 100644 --- a/GPy/core/gp_base.py +++ b/GPy/core/gp_base.py @@ -14,8 +14,11 @@ class GPBase(Model): Here we define some functions that are use """ def __init__(self, X, likelihood, kernel, normalize_X=False): + if len(X.shape)==1: + X = X.reshape(-1,1) + warning.warn("One dimension output (N,) being reshaped to (N,1)") self.X = X - assert len(self.X.shape) == 2 + assert len(self.X.shape) == 2, "too many dimensions for X input" self.num_data, self.input_dim = self.X.shape assert isinstance(kernel, kern.kern) self.kern = kernel diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index b60c7479..44f4ae3f 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -464,7 +464,7 @@ def symmetric(k): This constructor builds a covariance function of this form from the initial kernel """ k_ = k.copy() - k_.parts = [symmetric.Symmetric(p) for p in k.parts] + k_.parts = [parts.symmetric.Symmetric(p) for p in k.parts] return k_ def coregionalize(output_dim,rank=1, W=None, kappa=None): diff --git a/GPy/kern/parts/symmetric.py b/GPy/kern/parts/symmetric.py index bbdd5ac0..d836763d 100644 --- a/GPy/kern/parts/symmetric.py +++ b/GPy/kern/parts/symmetric.py @@ -56,7 +56,7 @@ class Symmetric(Kernpart): AX = np.dot(X,self.transform) if X2 is None: X2 = X - ZX2 = AX + AX2 = AX else: AX2 = np.dot(X2, self.transform) self.k.dK_dtheta(dL_dK,X,X2,target) 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/models/__init__.py b/GPy/models/__init__.py index 10ce577b..a8be5890 100644 --- a/GPy/models/__init__.py +++ b/GPy/models/__init__.py @@ -1,18 +1,19 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) -from gp_regression import GPRegression -from gp_classification import GPClassification -from sparse_gp_regression import SparseGPRegression -from svigp_regression import SVIGPRegression -from sparse_gp_classification import SparseGPClassification -from fitc_classification import FITCClassification -from gplvm import GPLVM -from bcgplvm import BCGPLVM -from sparse_gplvm import SparseGPLVM -from warped_gp import WarpedGP -from bayesian_gplvm import BayesianGPLVM -from mrd import MRD -from gradient_checker import GradientChecker -from gp_multioutput_regression import GPMultioutputRegression -from sparse_gp_multioutput_regression import SparseGPMultioutputRegression +from gp_regression import GPRegression; _gp_regression = gp_regression ; del gp_regression +from gp_classification import GPClassification; _gp_classification = gp_classification ; del gp_classification +from sparse_gp_regression import SparseGPRegression; _sparse_gp_regression = sparse_gp_regression ; del sparse_gp_regression +from svigp_regression import SVIGPRegression; _svigp_regression = svigp_regression ; del svigp_regression +from sparse_gp_classification import SparseGPClassification; _sparse_gp_classification = sparse_gp_classification ; del sparse_gp_classification +from fitc_classification import FITCClassification; _fitc_classification = fitc_classification ; del fitc_classification +from gplvm import GPLVM; _gplvm = gplvm ; del gplvm +from bcgplvm import BCGPLVM; _bcgplvm = bcgplvm; del bcgplvm +from sparse_gplvm import SparseGPLVM; _sparse_gplvm = sparse_gplvm ; del sparse_gplvm +from warped_gp import WarpedGP; _warped_gp = warped_gp ; del warped_gp +from bayesian_gplvm import BayesianGPLVM; _bayesian_gplvm = bayesian_gplvm ; del bayesian_gplvm +from mrd import MRD; _mrd = mrd ; del mrd +from gradient_checker import GradientChecker; _gradient_checker = gradient_checker ; del gradient_checker +from gp_multioutput_regression import GPMultioutputRegression; _gp_multioutput_regression = gp_multioutput_regression ; del gp_multioutput_regression +from sparse_gp_multioutput_regression import SparseGPMultioutputRegression; _sparse_gp_multioutput_regression = sparse_gp_multioutput_regression ; del sparse_gp_multioutput_regression + diff --git a/GPy/util/config.py b/GPy/util/config.py index d2ed7543..02796e0b 100644 --- a/GPy/util/config.py +++ b/GPy/util/config.py @@ -5,13 +5,18 @@ import ConfigParser import os config = ConfigParser.ConfigParser() -user_file = os.path.join(os.getenv('HOME'),'.gpy_config.cfg') -default_file = os.path.join('..','gpy_config.cfg') +home = os.getenv('HOME') or os.getenv('USERPROFILE') +user_file = os.path.join(home,'.gpy_config.cfg') +default_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', 'gpy_config.cfg')) +print user_file, os.path.isfile(user_file) +print default_file, os.path.isfile(default_file) # 1. check if the user has a ~/.gpy_config.cfg if os.path.isfile(user_file): config.read(user_file) -else: +elif os.path.isfile(default_file): # 2. if not, use the default one - path = os.path.dirname(__file__) - config.read(os.path.join(path,default_file)) + config.read(default_file) +else: + #3. panic + raise ValueError, "no configuration file found" 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 diff --git a/MANIFEST.in b/MANIFEST.in index c89284cd..8d5b2304 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -2,3 +2,5 @@ include *.txt recursive-include doc *.txt include *.md recursive-include doc *.md +include *.cfg +recursive-include doc *.cfg