Go to the documentation of this file.
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;
105 matrix_mode=matrix_cholesky;
116 static const size_t matrix_cholesky=0;
118 static const size_t matrix_LU=1;
122 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
123 O2SCL_ERR2(
"Function set(size_t,vec_t,vec_t) unimplemented ",
132 const vec_t &y, covar_func_t &fcovar,
133 covar_func_t &fderiv,
134 covar_func_t &fderiv2,
135 covar_func_t &finteg,
136 double noise_var,
bool err_on_fail=
true) {
137 O2SCL_ERR(
"Function set_covar_di_noise not yet implemented.",
144 covar_func_t &fcovar,
double noise_var,
145 bool err_on_fail=
true) {
147 if (n_dim<this->min_size) {
149 " than "+
szttos(this->min_size)+
" in interp_krige::"+
158 for(
size_t irow=0;irow<n_dim;irow++) {
159 for(
size_t icol=0;icol<n_dim;icol++) {
161 KXX(irow,icol)=KXX(icol,irow);
162 }
else if (irow==icol) {
163 KXX(irow,icol)=fcovar(x[irow],x[icol])+noise_var;
165 KXX(irow,icol)=fcovar(x[irow],x[icol]);
170 if (matrix_mode==matrix_LU) {
180 "in interp_krige::set_covar_noise().",
185 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
186 (n_dim,KXX,p,inv_KXX);
189 this->Kinvf.resize(n_dim);
190 boost::numeric::ublas::axpy_prod(inv_KXX,y,this->Kinvf,
true);
198 O2SCL_ERR2(
"Matrix singular (Cholesky method) ",
199 "in interp_krige::set_covar_noise().",
205 o2scl_linalg::cholesky_invert<ubmatrix>(n_dim,inv_KXX);
209 boost::numeric::ublas::axpy_prod(inv_KXX,y,Kinvf,
true);
222 virtual int set_covar(
size_t n_dim,
const vec_t &x,
const vec_t &y,
223 covar_func_t &fcovar,
bool err_on_fail=
true) {
224 return set_covar_noise(n_dim,x,y,fcovar,0.0,err_on_fail);
228 virtual double eval(
double x0)
const {
232 for(
size_t i=0;i<this->sz;i++) {
233 ret+=(*f)(x0,(*this->px)[i])*Kinvf[i];
240 virtual double deriv(
double x0)
const {
252 virtual double integ(
double a,
double b)
const {
257 virtual const char *
type()
const {
return "interp_krige"; }
259 #ifndef DOXYGEN_INTERNAL
280 template<
class vec_t,
class vec2_t=vec_t>
287 typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
292 std::function<double(
double,
double)>
ff;
334 size_t size=this->
sz;
339 for(
size_t k=0;k<size;k++) {
349 for(
size_t irow=0;irow<size-1;irow++) {
350 for(
size_t icol=0;icol<size-1;icol++) {
352 KXX(irow,icol)=KXX(icol,irow);
354 KXX(irow,icol)=exp(-pow((
x2[irow]-
x2[icol])/
len,2.0)/2.0);
355 if (irow==icol) KXX(irow,icol)+=noise_var;
363 ubmatrix inv_KXX(size-1,size-1);
371 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
372 (size-1,KXX,p,inv_KXX);
375 this->
Kinvf.resize(size-1);
376 boost::numeric::ublas::axpy_prod(inv_KXX,y2,this->
Kinvf,
true);
387 o2scl_linalg::cholesky_invert<ubmatrix>(size-1,inv_KXX);
390 this->
Kinvf.resize(size-1);
391 boost::numeric::ublas::axpy_prod(inv_KXX,y2,this->
Kinvf,
true);
396 double yact=(*this->
py)[k];
397 for(
size_t i=0;i<size-1;i++) {
398 ypred+=exp(-pow(((*this->
px)[k]-
x2[i])/
len,2.0)/2.0)*this->
Kinvf[i];
402 qual+=pow(yact-ypred,2.0);
410 for(
size_t irow=0;irow<size;irow++) {
411 for(
size_t icol=0;icol<size;icol++) {
413 KXX(irow,icol)=KXX(icol,irow);
415 KXX(irow,icol)=exp(-pow(((*this->
px)[irow]-
416 (*this->
px)[icol])/
len,2.0)/2.0);
417 if (irow==icol) KXX(irow,icol)+=noise_var;
434 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
435 (size,KXX,p,inv_KXX);
437 double lndet=o2scl_linalg::LU_lndet<ubmatrix>(size,KXX);
440 this->
Kinvf.resize(size);
441 boost::numeric::ublas::axpy_prod(inv_KXX,*this->
py,this->
Kinvf,
true);
445 for(
size_t i=0;i<size;i++) {
489 virtual int set_noise(
size_t size,
const vec_t &x,
const vec2_t &y,
490 double noise_var,
bool err_on_fail=
true) {
502 std::cout <<
"interp_krige_optim: full minimization"
507 double len_opt=x[1]-x[0];
510 (std::mem_fn<
double(
double,
double,
int &)>
512 std::placeholders::_1,noise_var,std::ref(
success));
520 "interp_krige_optim::set_noise().",
528 std::cout <<
"interp_krige_optim: simple minimization"
533 std::vector<double> diff(size-1);
534 for(
size_t i=0;i<size-1;i++) {
535 diff[i]=fabs(x[i+1]-x[i]);
540 <std::vector<double>,
double>(size-1,diff)/3.0;
541 double len_max=fabs(x[size-1]-x[0])*3.0;
542 double len_ratio=len_max/len_min;
545 std::cout <<
"len (min,max,ratio): " << len_min <<
" "
547 << pow(len_ratio,((
double)1)/((
double)
nlen-1))
552 double min_qual=0.0, len_opt=0.0;
555 std::cout <<
"ilen len qual fail min_qual" << std::endl;
560 for(
size_t j=0;j<
nlen;j++) {
561 len=len_min*pow(len_ratio,((
double)j)/((
double)
nlen-1));
566 if (
success==0 && (min_set==
false ||
qual<min_qual)) {
573 std::cout <<
"interp_krige_optim: ";
575 std::cout << j <<
" " <<
len <<
" " <<
qual <<
" "
576 <<
success <<
" " << min_qual <<
" "
577 << len_opt << std::endl;
588 ff=std::bind(std::mem_fn<
double(
double,
double)>
590 std::placeholders::_1,std::placeholders::_2);
598 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
602 for(
size_t j=0;j<size;j++) {
603 mean_abs+=fabs(y[j]);
615 #ifndef DOXYGEN_NO_O2NS
@ exc_esing
apparent singularity detected
int verbose
Verbosity parameter.
static const size_t matrix_LU
Use LU decomposition.
virtual int set_noise(size_t size, const vec_t &x, const vec2_t &y, double noise_var, bool err_on_fail=true)
Initialize interpolation routine.
const vec_t * py
Dependent vector.
A class for representing permutations.
covar_func_t * df2
Pointer to user-specified second derivative.
double integ(double x, double x1, double x2)
The integral of the covariance function.
covar_integ_t * intp
Pointer to user-specified integral.
virtual double deriv2(double x0) const
Give the value of the second derivative .
@ exc_efailed
generic failure
virtual double eval(double x0) const
Give the value of the function .
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
virtual const char * type() const
Return the type, "interp_krige".
min_base * mp
Pointer to the user-specified minimizer.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
double len
The covariance function length scale.
data_t vector_min_value(size_t n, const vec_t &data)
Compute the minimum of the first n elements of a vector.
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
@ exc_eunimpl
requested feature not (yet) implemented
size_t matrix_mode
Method for matrix inversion.
std::function< double(double, double)> ff
Function object for the covariance.
static const double x2[5]
double deriv(double x1, double x2)
The derivative of the covariance function.
virtual void set(size_t size, const vec_t &x, const vec2_t &y)
Initialize interpolation routine.
static const size_t mode_loo_cv
Leave-one-out cross validation.
One-dimensional interpolation using an optimized covariance function.
virtual double integ(double a, double b) const
Give the value of the integral .
static const double x1[5]
@ exc_einval
invalid argument supplied by user
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 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, bool err_on_fail=true)
Initialize interpolation routine, specifying derivatives and integrals [not yet implemented].
double deriv2(double x1, double x2)
The second derivative of the covariance function.
bool full_min
If true, use the full minimizer.
size_t nlen
Number of length scale points to try when full minimizer is not used (default 20)
static const size_t mode_max_lml
Minus Log-marginal-likelihood.
virtual void set(size_t size, const vec_t &x, const vec2_t &y)
Initialize interpolation routine.
virtual double deriv(double x0) const
Give the value of the 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.
double covar(double x1, double x2)
The covariance function.
Interpolation by Kriging with a user-specified covariance function.
min_brent_gsl def_min
Default minimizer.
double qual
The quality factor of the optimization.
One-dimensional minimization [abstract base].
covar_func_t * f
Pointer to user-specified covariance function.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
std::function< double(double)> funct
One-dimensional function typedef in src/base/funct.h.
covar_func_t * df
Pointer to user-specified derivative.
virtual int min(double &x, double &fmin, func_t &func)=0
Calculate the minimum min of func w.r.t 'x'.
virtual int set_covar(size_t n_dim, const vec_t &x, const vec_t &y, covar_func_t &fcovar, bool err_on_fail=true)
Initialize interpolation routine.
double qual_fun(double x, double noise_var, int &success)
Function to optimize the covariance parameters.
const vec_t * px
Independent vector.
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.
std::string szttos(size_t x)
Convert a size_t to a string.
ubvector Kinvf
Inverse covariance matrix times function vector.
Base low-level interpolation class [abstract base].
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.
size_t mode
Function to minimize (default mode_loo_cv)
One-dimensional minimization using Brent's method (GSL)
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).