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> 42 #include <o2scl/columnify.h> 43 #include <o2scl/vector.h> 44 #include <o2scl/vec_stats.h> 45 #include <o2scl/mmin_bfgs2.h> 46 #include <o2scl/constants.h> 48 #ifndef DOXYGEN_NO_O2NS 63 template<
class vec_t,
class vec2_t=vec_t,
64 class covar_func_t=std::function<double(double,double)>,
65 class covar_integ_t=std::function<double(double,double,double)> >
68 #ifdef O2SCL_NEVER_DEFINED 76 typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
109 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
110 O2SCL_ERR2(
"Function set(size_t,vec_t,vec_t) unimplemented ",
119 const vec_t &y, covar_func_t &fcovar,
120 covar_func_t &fderiv,
121 covar_func_t &fderiv2,
122 covar_func_t &finteg,
129 covar_func_t &fcovar,
double noise_var) {
133 " than "+
szttos(this->min_size)+
" in interp_krige::"+
141 ubmatrix KXX(n_dim,n_dim);
142 for(
size_t irow=0;irow<n_dim;irow++) {
143 for(
size_t icol=0;icol<n_dim;icol++) {
145 KXX(irow,icol)=KXX(icol,irow);
146 }
else if (irow==icol) {
147 KXX(irow,icol)=fcovar(x[irow],x[icol])+noise_var;
149 KXX(irow,icol)=fcovar(x[irow],x[icol]);
155 ubmatrix inv_KXX(n_dim,n_dim);
160 O2SCL_ERR(
"KXX matrix is singular in interp_krige::set().",
163 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
164 (n_dim,KXX,p,inv_KXX);
168 boost::numeric::ublas::axpy_prod(inv_KXX,y,Kinvf,
true);
178 virtual void set_covar(
size_t n_dim,
const vec_t &x,
const vec_t &y,
179 covar_func_t &fcovar) {
185 virtual double eval(
double x0)
const {
189 for(
size_t i=0;i<this->
sz;i++) {
190 ret+=(*f)(x0,(*this->
px)[i])*Kinvf[i];
197 virtual double deriv(
double x0)
const {
209 virtual double integ(
double a,
double b)
const {
214 virtual const char *
type()
const {
return "interp_krige"; }
216 #ifndef DOXYGEN_INTERNAL 223 (
const interp_krige<vec_t,vec2_t,covar_func_t,covar_integ_t>&);
237 template<
class vec_t,
class vec2_t=vec_t>
244 typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
259 return var*exp(-(x1-x2)*(x1-x2)/len/len);
264 return -2.0*var*exp(-(x1-x2)*(x1-x2)/len/len)/len/len*(x1-
x2);
269 return (4.0*(x1-x2)*(x1-x2)-2.0*len*len)*
270 var*exp(-(x1-x2)*(x1-x2)/len/len)/len/len/len/len;
276 (gsl_sf_erf((x2-x)/len)+gsl_sf_erf((x-x1)/len));
280 std::function<double(double,double)>
ff;
292 size_t size=this->
sz;
295 for(
size_t k=0;k<size;k++) {
303 ubmatrix KXX(size-1,size-1);
304 for(
size_t irow=0;irow<size-1;irow++) {
305 for(
size_t icol=0;icol<size-1;icol++) {
307 KXX(irow,icol)=KXX(icol,irow);
309 KXX(irow,icol)=exp(-var*pow(((*this->
px)[irow]-
310 (*this->
px)[icol])/len,2.0));
316 ubmatrix inv_KXX(size-1,size-1);
321 O2SCL_ERR(
"KXX matrix is singular in interp_krige::set().",
324 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
325 (size-1,KXX,p,inv_KXX);
328 this->
Kinvf.resize(size-1);
329 boost::numeric::ublas::axpy_prod(inv_KXX,y2,this->
Kinvf,
true);
332 double yact=(*this->
py)[k];
333 for(
size_t i=0;i<size-1;i++) {
334 ypred+=exp(-var*pow(((*this->
px)[k]-x2[i])/len,2.0))*this->
Kinvf[i];
337 qual+=pow(yact-ypred,2.0);
359 ff=std::bind(std::mem_fn<
double(
double,
double)>
361 std::placeholders::_1,std::placeholders::_2);
369 virtual void set_noise(
size_t size,
const vec_t &x,
const vec2_t &y,
384 (std::mem_fn<
double(
size_t,
const ubvector &)>
386 std::placeholders::_1,std::placeholders::_2);
388 mp->
mmin(2,p,qual,mf);
394 std::vector<double> diff(size-1);
395 for(
size_t i=0;i<size-1;i++) {
396 diff[i]=fabs(x[i+1]-x[i]);
398 double var_ratio=1.0e2;
402 <std::vector<double>,
double>(size-1,diff)/3.0;
403 double len_max=fabs(x[size-1]-x[0])*3.0;
404 double len_ratio=len_max/len_min;
407 double min_qual=0.0, var_opt=0.0, len_opt=0.0;
410 for(
size_t i=0;i<nvar;i++) {
411 var=var_min*pow(var_ratio,((
double)i)/((
double)nvar-1));
412 for(
size_t j=0;j<nlen;j++) {
413 len=len_min*pow(len_ratio,((
double)j)/((
double)nlen-1));
416 for(
size_t k=0;k<size;k++) {
424 ubmatrix KXX(size-1,size-1);
425 for(
size_t irow=0;irow<size-1;irow++) {
426 for(
size_t icol=0;icol<size-1;icol++) {
428 KXX(irow,icol)=KXX(icol,irow);
430 KXX(irow,icol)=exp(-var*pow((x[irow]-x[icol])/len,2.0));
436 ubmatrix inv_KXX(size-1,size-1);
441 O2SCL_ERR(
"KXX matrix is singular in interp_krige::set().",
444 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
445 (size-1,KXX,p,inv_KXX);
448 this->
Kinvf.resize(size-1);
449 boost::numeric::ublas::axpy_prod(inv_KXX,y2,this->
Kinvf,
true);
453 for(
size_t i=0;i<size-1;i++) {
454 ypred+=exp(-var*pow((x[k]-x2[i])/len,2.0))*this->
Kinvf[i];
457 qual+=pow(yact-ypred,2.0);
461 if ((i==0 && j==0) || qual<min_qual) {
482 virtual void set(
size_t size,
const vec_t &x,
const vec2_t &y) {
483 set_noise(size,x,y,0.0);
490 #ifndef DOXYGEN_NO_O2NS std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct
Multi-dimensional function typedef.
data_t vector_min_value(size_t n, const vec_t &data)
Compute the minimum of the first n elements of a vector.
virtual void 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.
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 .
size_t nvar
Number of variance points to try.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
double qual_fun(size_t nv, const ubvector &x)
Function to optimize the covariance.
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.
covar_func_t * f
Pointer to user-specified covariance function.
mmin_bfgs2 def_mmin
Default minimizer.
std::function< double(double, double)> ff
Functor for the covariance function covar()
invalid argument supplied by user
virtual int mmin(size_t nvar, vec_t &x, double &fmin, func_t &func)=0
Calculate the minimum min of func w.r.t. the array x of size nvar.
const vec_t * px
Independent vector.
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.
virtual void set_covar(size_t n_dim, const vec_t &x, const vec_t &y, covar_func_t &fcovar)
Initialize interpolation routine.
double vector_variance(size_t n, const vec_t &data, double mean)
Compute the variance with specified mean.
double len
The covariance function length scale.
requested feature not (yet) implemented
One-dimensional interpolation using an optimized covariance function.
double var
The covariance function coefficient.
virtual void set_noise(size_t size, const vec_t &x, const vec2_t &y, double noise_var)
Initialize interpolation routine.
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.
covar_integ_t * intp
Pointer to user-specified second derivative.
double deriv(double x1, double x2)
The derivative of the covariance function.
virtual void set_covar_noise(size_t n_dim, const vec_t &x, const vec_t &y, covar_func_t &fcovar, double noise_var)
Initialize interpolation routine.
Multidimensional minimization [abstract base].
mmin_base * mp
Pointer to the user-specified minimizer.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
const vec2_t * py
Dependent vector.
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.
double covar(double x1, double x2)
The covariance function.
virtual double deriv2(double x0) const
Give the value of the second derivative (always zero)
double qual
The quality factor of the optimization.
Multidimensional minimization by the BFGS algorithm (GSL)
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.
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.