23 #ifndef O2SCL_INTERP_KRIGE_H 24 #define O2SCL_INTERP_KRIGE_H 33 #include <gsl/gsl_sf_erf.h> 35 #include <boost/numeric/ublas/vector.hpp> 36 #include <boost/numeric/ublas/matrix.hpp> 37 #include <boost/numeric/ublas/operation.hpp> 38 #include <boost/numeric/ublas/vector_proxy.hpp> 40 #include <o2scl/interp.h> 41 #include <o2scl/cholesky.h> 43 #include <o2scl/columnify.h> 44 #include <o2scl/vector.h> 45 #include <o2scl/vec_stats.h> 46 #include <o2scl/min_brent_gsl.h> 47 #include <o2scl/constants.h> 49 #ifndef DOXYGEN_NO_O2NS 64 template<
class vec_t,
class vec2_t=vec_t,
65 class covar_func_t=std::function<double(double,double)>,
66 class covar_integ_t=std::function<double(double,double,double)> >
69 #ifdef O2SCL_NEVER_DEFINED 77 typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
114 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
115 O2SCL_ERR2(
"Function set(size_t,vec_t,vec_t) unimplemented ",
124 const vec_t &y, covar_func_t &fcovar,
125 covar_func_t &fderiv,
126 covar_func_t &fderiv2,
127 covar_func_t &finteg,
129 O2SCL_ERR(
"Function set_covar_di_noise not yet implemented.",
136 covar_func_t &fcovar,
double noise_var,
137 bool err_on_fail=
true) {
141 " than "+
szttos(this->min_size)+
" in interp_krige::"+
149 ubmatrix KXX(n_dim,n_dim);
150 for(
size_t irow=0;irow<n_dim;irow++) {
151 for(
size_t icol=0;icol<n_dim;icol++) {
153 KXX(irow,icol)=KXX(icol,irow);
154 }
else if (irow==icol) {
155 KXX(irow,icol)=fcovar(x[irow],x[icol])+noise_var;
157 KXX(irow,icol)=fcovar(x[irow],x[icol]);
162 if (matrix_mode!=0) {
165 ubmatrix inv_KXX(n_dim,n_dim);
172 "in interp_krige::set_covar_noise().",
177 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
178 (n_dim,KXX,p,inv_KXX);
181 this->Kinvf.resize(n_dim);
182 boost::numeric::ublas::axpy_prod(inv_KXX,y,this->Kinvf,
true);
190 O2SCL_ERR2(
"Matrix singular (Cholesky method) ",
191 "in interp_krige::set_covar_noise().",
196 ubmatrix inv_KXX=KXX;
197 o2scl_linalg::cholesky_invert<ubmatrix>(n_dim,inv_KXX);
201 boost::numeric::ublas::axpy_prod(inv_KXX,y,Kinvf,
true);
214 virtual int set_covar(
size_t n_dim,
const vec_t &x,
const vec_t &y,
215 covar_func_t &fcovar) {
220 virtual double eval(
double x0)
const {
224 for(
size_t i=0;i<this->
sz;i++) {
225 ret+=(*f)(x0,(*this->
px)[i])*Kinvf[i];
232 virtual double deriv(
double x0)
const {
244 virtual double integ(
double a,
double b)
const {
249 virtual const char *
type()
const {
return "interp_krige"; }
251 #ifndef DOXYGEN_INTERNAL 258 (
const interp_krige<vec_t,vec2_t,covar_func_t,covar_integ_t>&);
272 template<
class vec_t,
class vec2_t=vec_t>
279 typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
283 std::function<double(double,double)> ff;
293 return exp(-(x1-x2)*(x1-x2)/len/len/2.0);
298 return -exp(-(x1-x2)*(x1-x2)/len/len/2.0)/len/len*(x1-
x2);
303 return ((x1-x2)*(x1-x2)-len*len)*
304 exp(-(x1-x2)*(x1-x2)/len/len/2.0)/len/len/len/len;
325 size_t size=this->
sz;
328 static const size_t loo_cv=1;
330 static const size_t lml=2;
336 for(
size_t k=0;k<size;k++) {
345 ubmatrix KXX(size-1,size-1);
346 for(
size_t irow=0;irow<size-1;irow++) {
347 for(
size_t icol=0;icol<size-1;icol++) {
349 KXX(irow,icol)=KXX(icol,irow);
351 KXX(irow,icol)=exp(-pow((x2[irow]-x2[icol])/len,2.0)/2.0);
352 if (irow==icol) KXX(irow,icol)+=noise_var;
360 ubmatrix inv_KXX(size-1,size-1);
368 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
369 (size-1,KXX,p,inv_KXX);
372 this->
Kinvf.resize(size-1);
373 boost::numeric::ublas::axpy_prod(inv_KXX,y2,this->
Kinvf,
true);
383 ubmatrix &inv_KXX=KXX;
384 o2scl_linalg::cholesky_invert<ubmatrix>(size-1,inv_KXX);
387 this->
Kinvf.resize(size-1);
388 boost::numeric::ublas::axpy_prod(inv_KXX,y2,this->
Kinvf,
true);
393 double yact=(*this->
py)[k];
394 for(
size_t i=0;i<size-1;i++) {
395 ypred+=exp(-pow(((*this->
px)[k]-x2[i])/len,2.0)/2.0)*this->
Kinvf[i];
399 qual+=pow(yact-ypred,2.0);
403 }
else if (mode==lml) {
406 ubmatrix KXX(size,size);
407 for(
size_t irow=0;irow<size;irow++) {
408 for(
size_t icol=0;icol<size;icol++) {
410 KXX(irow,icol)=KXX(icol,irow);
412 KXX(irow,icol)=exp(-pow(((*this->
px)[irow]-
413 (*this->
px)[icol])/len,2.0)/2.0);
414 if (irow==icol) KXX(irow,icol)+=noise_var;
423 ubmatrix inv_KXX(size,size);
431 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
432 (size,KXX,p,inv_KXX);
434 double lndet=o2scl_linalg::LU_lndet<ubmatrix>(size,KXX);
437 this->
Kinvf.resize(size);
438 boost::numeric::ublas::axpy_prod(inv_KXX,*this->
py,this->
Kinvf,
true);
442 for(
size_t i=0;i<size;i++) {
443 qual+=-0.5*(*this->
py)[i]*this->
Kinvf[i];
473 virtual int set_noise(
size_t size,
const vec_t &x,
const vec2_t &y,
487 (std::mem_fn<
double(
double,
double,
int &)>
489 std::placeholders::_1,noise_var,std::ref(success));
496 std::cout <<
"interp_krige_optim: simple minimization" 501 std::vector<double> diff(size-1);
502 for(
size_t i=0;i<size-1;i++) {
503 diff[i]=fabs(x[i+1]-x[i]);
506 <std::vector<double>,
double>(size-1,diff)/3.0;
507 double len_max=fabs(x[size-1]-x[0])*3.0;
508 double len_ratio=len_max/len_min;
511 std::cout <<
"len: " << len_min <<
" " 513 << pow(len_ratio,((
double)1)/((
double)nlen-1))
518 double min_qual=0.0, len_opt=0.0;
521 std::cout <<
"ilen len qual fail min_qual" << std::endl;
526 for(
size_t j=0;j<nlen;j++) {
527 len=len_min*pow(len_ratio,((
double)j)/((
double)nlen-1));
530 qual=qual_fun(len,noise_var,success);
532 if (success==0 && (min_set==
false || qual<min_qual)) {
539 std::cout <<
"interp_krige_optim: ";
541 std::cout << j <<
" " << len <<
" " << qual <<
" " 542 << success <<
" " << min_qual <<
" " 543 << len_opt << std::endl;
554 ff=std::bind(std::mem_fn<
double(
double,
double)>
556 std::placeholders::_1,std::placeholders::_2);
564 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
566 for(
size_t j=0;j<size;j++) {
567 mean_abs+=fabs(y[j]);
570 set_noise(size,x,y,mean_abs/1.0e8);
577 #ifndef DOXYGEN_NO_O2NS double qual_fun(double x, double noise_var, int &success)
Function to optimize the covariance parameters.
data_t vector_min_value(size_t n, const vec_t &data)
Compute the minimum of the first n elements of a vector.
covar_func_t * df
Pointer to user-specified derivative.
void vector_copy_jackknife(const vec_t &src, size_t iout, vec2_t &dest)
From a given vector, create a new vector by removing a specified element.
virtual double integ(double a, double b) const
Give the value of the integral .
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
size_t min_size
The minimum size of the vectors to interpolate between.
double integ(double x, double x1, double x2)
The integral of the covariance function.
virtual int set_covar_noise(size_t n_dim, const vec_t &x, const vec_t &y, covar_func_t &fcovar, double noise_var, bool err_on_fail=true)
Initialize interpolation routine.
covar_func_t * f
Pointer to user-specified covariance function.
min_base * mp
Pointer to the user-specified minimizer.
int matrix_mode
Method for matrix inversion.
invalid argument supplied by user
const vec_t * px
Independent vector.
apparent singularity detected
A class for representing permutations.
Base low-level interpolation class [abstract base].
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
Interpolation by Kriging with a user-specified covariance function.
double len
The covariance function length scale.
requested feature not (yet) implemented
One-dimensional interpolation using an optimized covariance function.
virtual double deriv(double x0) const
Give the value of the derivative .
double deriv2(double x1, double x2)
The second derivative of the covariance function.
covar_func_t * df2
Pointer to user-specified second derivative.
One-dimensional minimization using Brent's method (GSL)
covar_integ_t * intp
Pointer to user-specified second derivative.
double deriv(double x1, double x2)
The derivative of the covariance function.
min_brent_gsl def_min
Default minimizer.
virtual int set_covar(size_t n_dim, const vec_t &x, const vec_t &y, covar_func_t &fcovar)
Initialize interpolation routine.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
virtual int set_noise(size_t size, const vec_t &x, const vec2_t &y, double noise_var)
Initialize interpolation routine.
const vec2_t * py
Dependent vector.
int cholesky_decomp(const size_t M, mat_t &A, bool err_on_fail=true)
Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix...
virtual int min(double &x, double &fmin, func_t &func)=0
Calculate the minimum min of func w.r.t 'x'.
virtual const char * type() const
Return the type, "interp_linear".
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
bool full_min
If true, use the full minimizer.
ubvector Kinvf
Inverse covariance matrix times function vector.
size_t nlen
Number of length scale points to try.
One-dimensional minimization [abstract base].
double covar(double x1, double x2)
The covariance function.
virtual double deriv2(double x0) const
Give the value of the second derivative (always zero)
std::function< double(double)> funct
One-dimensional function typedef.
double qual
The quality factor of the optimization.
virtual int set_covar_di_noise(size_t n_dim, const vec_t &x, const vec_t &y, covar_func_t &fcovar, covar_func_t &fderiv, covar_func_t &fderiv2, covar_func_t &finteg, double noise_var)
Initialize interpolation routine, specifying derivatives and integrals [not yet implemented].
int diagonal_has_zero(const size_t N, mat_t &A)
Return 1 if at least one of the elements in the diagonal is zero.
int verbose
Verbosity parameter.
static const double x2[5]
virtual double eval(double x0) const
Give the value of the function .
static const double x1[5]
std::string szttos(size_t x)
Convert a size_t to a string.