23 #ifndef O2SCL_INTERPM_KRIGE_H 24 #define O2SCL_INTERPM_KRIGE_H 35 #include <boost/numeric/ublas/matrix.hpp> 36 #include <boost/numeric/ublas/operation.hpp> 38 #include <gsl/gsl_combination.h> 40 #include <o2scl/err_hnd.h> 41 #include <o2scl/vector.h> 42 #include <o2scl/vec_stats.h> 43 #include <o2scl/linear_solver.h> 44 #include <o2scl/columnify.h> 46 #ifndef DOXYGEN_NO_O2NS 59 template<
class vec_t=boost::numeric::ublas::vector<
double>,
60 class mat_col_t=boost::numeric::ublas::matrix_column<
61 boost::numeric::ublas::matrix<
double> >,
62 class covar_func_t=std::function<
double(const vec_t &,const vec_t &)> >
90 template<
class vec_vec_t,
class vec_vec2_t>
91 void set_data(
size_t n_in,
size_t n_out,
size_t n_points,
92 vec_vec_t &x, vec_vec2_t &y,
93 std::vector<covar_func_t> &fcovar) {
96 O2SCL_ERR2(
"Must provide at least three points in ",
100 O2SCL_ERR2(
"Must provide at least one input column in ",
104 O2SCL_ERR2(
"Must provide at least one output column in ",
111 for(
size_t i=0;i<n_points;i++) {
112 if (x[i].size()!=n_in) {
116 std::swap(
ptrs_x[i],x[i]);
121 std::cout <<
"interpm_krige::set_data() : Using " << n_points
122 <<
" points with " <<
nd_in <<
" input variables and\n\t" 123 <<
nd_out <<
" output variables." << std::endl;
129 for(
size_t iout=0;iout<n_out;iout++) {
132 ubmatrix KXX(n_points,n_points);
133 for(
size_t irow=0;irow<n_points;irow++) {
134 for(
size_t icol=0;icol<n_points;icol++) {
136 KXX(irow,icol)=KXX(icol,irow);
138 KXX(irow,icol)=fcovar[iout](
ptrs_x[irow],ptrs_x[icol]);
144 ubmatrix inv_KXX(n_points,n_points);
148 std::cout <<
"interpm_krige::set_data() : " 149 <<
"LU decompose and invert " << iout+1 <<
" of " << n_out
155 "interpm_krige::set_data().",
158 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,mat_col_t>
159 (n_points,KXX,p,inv_KXX);
162 Kinvf[iout].resize(n_points);
163 boost::numeric::ublas::axpy_prod(inv_KXX,y[iout],Kinvf[iout],
true);
172 template<
class vec2_t,
class vec3_t>
173 void eval(
const vec2_t &x, vec3_t &y,
174 std::vector<covar_func_t> &fcovar)
const {
177 O2SCL_ERR(
"Data not set in interpm_krige::eval().",
182 for(
size_t iout=0;iout<
nd_out;iout++) {
184 for(
size_t ipoints=0;ipoints<
np;ipoints++) {
185 y[iout]+=fcovar[iout](x,
ptrs_x[ipoints])*Kinvf[iout][ipoints];
193 #ifndef DOXYGEN_INTERNAL 224 template<
class vec_t=boost::numeric::ublas::vector<
double>,
225 class mat_col_t=boost::numeric::ublas::matrix_column<
226 boost::numeric::ublas::matrix<
double> >,
227 class covar_func_t=std::function<
double(const vec_t &,const vec_t &)> >
247 template<
class vec_vec_t,
class vec_vec2_t>
248 void set_data(
size_t n_in,
size_t n_out,
size_t n_points,
249 vec_vec_t &x, vec_vec2_t &y,
250 std::vector<covar_func_t> &fcovar,
size_t order) {
253 O2SCL_ERR2(
"Must provide at least three points in ",
257 O2SCL_ERR2(
"Must provide at least one input column in ",
261 O2SCL_ERR2(
"Must provide at least one output column in ",
269 for(
size_t i=0;i<n_points;i++) {
270 if (x[i].size()!=n_in) {
274 std::swap(
ptrs_x[i],x[i]);
276 ptrs_y.resize(n_out);
277 for(
size_t i=0;i<n_out;i++) {
278 if (y[i].size()!=n_points) {
282 std::swap(ptrs_y[i],y[i]);
287 std::cout <<
"interpm_krige_nn::set_data() : Using " << n_points
288 <<
" points with " <<
nd_in <<
" input variables and\n\t" 289 <<
nd_out <<
" output variables and order " 290 << norder <<
" ." << std::endl;
298 template<
class vec2_t,
class vec3_t>
299 void eval(
const vec2_t &x, vec3_t &y,
300 std::vector<covar_func_t> &fcovar)
const {
303 O2SCL_ERR(
"Data not set in interpm_krige::eval().",
310 for(
size_t iout=0;iout<
nd_out;iout++) {
316 for(
size_t ip=0;ip<
np;ip++) {
317 dists[ip]=-fcovar[iout](x,
ptrs_x[ip]);
322 ubvector_size_t index;
323 o2scl::vector_smallest_index<ubvector,double,ubvector_size_t>
324 (
np,dists,norder,index);
327 ubvector func(norder);
328 for(
size_t io=0;io<norder;io++) {
329 func[io]=ptrs_y[iout][index[io]];
333 ubmatrix KXX(norder,norder);
334 for(
size_t irow=0;irow<norder;irow++) {
335 for(
size_t icol=0;icol<norder;icol++) {
337 KXX(irow,icol)=KXX(icol,irow);
339 KXX(irow,icol)=fcovar[iout](
ptrs_x[index[irow]],
346 ubmatrix inv_KXX(norder,norder);
352 "interpm_krige_nn::eval().",
355 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,mat_col_t>
356 (norder,KXX,p,inv_KXX);
359 ubvector
Kinvf(norder);
360 boost::numeric::ublas::axpy_prod(inv_KXX,func,Kinvf,
true);
364 for(
size_t ipoints=0;ipoints<norder;ipoints++) {
365 y[iout]+=-dists[index[ipoints]]*Kinvf[ipoints];
376 template<
class vec2_t>
378 std::vector<covar_func_t> &fcovar,
379 ubvector_size_t &index,
380 ubvector_size_t &indep)
const {
382 if (x.size()<
nd_in || index.size()<
np || indep.size()<norder
384 std::cout << x.size() <<
" " <<
nd_in << std::endl;
385 std::cout << index.size() <<
" " <<
np << std::endl;
386 std::cout << indep.size() <<
" " << norder << std::endl;
387 std::cout << iout <<
" " <<
nd_out << std::endl;
388 O2SCL_ERR(
"Vectors not of correct size in find_lin_indep().",
393 while (done==
false) {
396 ubvector func(norder);
397 for(
size_t io=0;io<norder;io++) {
398 func[io]=ptrs_y[iout][index[indep[io]]];
402 ubmatrix KXX(norder,norder);
403 for(
size_t irow=0;irow<norder;irow++) {
404 for(
size_t icol=0;icol<norder;icol++) {
406 KXX(irow,icol)=KXX(icol,irow);
408 KXX(irow,icol)=fcovar[iout](
ptrs_x[index[indep[irow]]],
409 ptrs_x[index[indep[icol]]]);
422 std::cout <<
"Finding new independent rows." << std::endl;
423 for(
size_t j=0;j<norder;j++) {
424 std::cout << indep[j] <<
" " << KXX(j,j) << std::endl;
430 std::cout <<
"Max is: " << max << std::endl;
432 for(
size_t j=0;j<norder;j++) {
436 std::cout <<
"Entry " << j <<
" is zero so replacing " 437 <<
"entry with " << max+1 << std::endl;
442 O2SCL_ERR3(
"Failed to find set of independent points ",
443 "in interpm_krige_nn::find_lin_indep",
444 "(const vec2_t &, size_t).",
457 template<
class vec2_t>
458 double eval(
const vec2_t &x,
size_t iout,
459 std::vector<covar_func_t> &fcovar)
const {
462 O2SCL_ERR(
"Data not set in interpm_krige_nn::eval().",
471 for(
size_t ip=0;ip<
np;ip++) {
472 dists[ip]=-fcovar[iout](x,
ptrs_x[ip]);
475 ubvector_size_t index(np);
476 o2scl::vector_sort_index<ubvector,ubvector_size_t>(
np,dists,index);
481 ubvector_size_t indep(norder);
482 for(
size_t io=0;io<norder;io++) {
486 find_lin_indep(x,iout,fcovar,index,indep);
489 ubvector func(norder);
490 for(
size_t io=0;io<norder;io++) {
491 func[io]=ptrs_y[iout][index[indep[io]]];
495 ubmatrix KXX(norder,norder);
496 for(
size_t irow=0;irow<norder;irow++) {
497 for(
size_t icol=0;icol<norder;icol++) {
499 KXX(irow,icol)=KXX(icol,irow);
501 KXX(irow,icol)=fcovar[iout](
ptrs_x[index[indep[irow]]],
502 ptrs_x[index[indep[icol]]]);
508 ubmatrix inv_KXX(norder,norder);
514 "interpm_krige_nn::eval().",
517 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,mat_col_t>
518 (norder,KXX,p,inv_KXX);
521 ubvector
Kinvf(norder);
522 boost::numeric::ublas::axpy_prod(inv_KXX,func,Kinvf,
true);
526 for(
size_t ipoints=0;ipoints<norder;ipoints++) {
527 ret+=-dists[index[indep[ipoints]]]*Kinvf[ipoints];
537 template<
class vec2_t>
539 std::vector<covar_func_t> &fcovar)
const {
542 O2SCL_ERR(
"Data not set in interpm_krige::eval().",
554 for(
size_t ip=0;ip<
np;ip++) {
555 dists[ip]=-fcovar[iout](x,
ptrs_x[ip]);
559 ubvector_size_t index(np);
560 o2scl::vector_sort_index<ubvector,ubvector_size_t>(
np,dists,index);
565 ubvector_size_t indep(norder);
566 for(
size_t io=0;io<norder;io++) {
574 find_lin_indep(x,iout,fcovar,index,indep);
579 for(
size_t jk=0;jk<norder;jk++) {
582 std::cout <<
"Jackknife: " << jk <<
" matching function value " 583 << ptrs_y[iout][index[jk]] << std::endl;
586 ubvector_size_t indep_jk;
590 ubvector func(norder-1);
591 for(
size_t io=0;io<norder-1;io++) {
592 func[io]=ptrs_y[iout][index[indep_jk[io]]];
596 ubmatrix KXX(norder-1,norder-1);
597 for(
size_t irow=0;irow<norder-1;irow++) {
598 for(
size_t icol=0;icol<norder-1;icol++) {
600 KXX(irow,icol)=KXX(icol,irow);
602 KXX(irow,icol)=fcovar[iout](
ptrs_x[index[indep_jk[irow]]],
603 ptrs_x[index[indep_jk[icol]]]);
610 ubmatrix inv_KXX(norder-1,norder-1);
615 std::cout <<
"Second fli run starting." << std::endl;
616 find_lin_indep(x,iout,fcovar,index,indep);
617 std::cout <<
"Second fli run done." << std::endl;
620 O2SCL_ERR3(
"Failed to find set of independent points ",
621 "in interpm_krige_nn::eval_jackknife",
622 "(const vec2_t &, size_t).",
626 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,mat_col_t>
627 (norder-1,KXX,p,inv_KXX);
630 ubvector
Kinvf(norder-1);
631 boost::numeric::ublas::axpy_prod(inv_KXX,func,Kinvf,
true);
635 for(
size_t ipoints=0;ipoints<norder-1;ipoints++) {
636 ytmp+=-dists[index[indep_jk[ipoints]]]*Kinvf[ipoints];
640 qual+=pow(ptrs_y[iout][index[jk]]-ytmp,2.0);
643 std::cout <<
"Original value: " 644 << ptrs_y[iout][index[jk]] <<
" interpolated: " 645 << ytmp << std::endl;
655 #ifndef DOXYGEN_INTERNAL 680 #ifndef DOXYGEN_NO_O2NS bool data_set
True if the data has been specified.
size_t nd_out
The number of dimensions of the outputs.
void eval(const vec2_t &x, vec3_t &y, std::vector< covar_func_t > &fcovar) const
Perform the interpolation.
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.
Multi-dimensional interpolation by kriging with nearest-neighbor.
double eval_jackknife(const vec2_t &x, size_t iout, std::vector< covar_func_t > &fcovar) const
Use jackknife to match interpolated to real function values.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
int verbose
Verbosity parameter (default 0)
void set_data(size_t n_in, size_t n_out, size_t n_points, vec_vec_t &x, vec_vec2_t &y, std::vector< covar_func_t > &fcovar, size_t order)
Initialize the data for the interpolation.
void find_lin_indep(const vec2_t &x, size_t iout, std::vector< covar_func_t > &fcovar, ubvector_size_t &index, ubvector_size_t &indep) const
Find a set of linearly independent points.
size_t nd_out
The number of dimensions of the outputs.
Multi-dimensional interpolation by kriging.
void set_data(size_t n_in, size_t n_out, size_t n_points, vec_vec_t &x, vec_vec2_t &y, std::vector< covar_func_t > &fcovar)
Initialize the data for the interpolation.
invalid argument supplied by user
A class for representing permutations.
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
int verbose
Verbosity parameter (default 0)
size_t nd_in
The number of dimensions of the inputs.
std::vector< vec_t > ptrs_x
A vector of pointers holding the data.
double eval(const vec2_t &x, size_t iout, std::vector< covar_func_t > &fcovar) const
Perform the interpolation.
#define O2SCL_ERR3(d, d2, d3, n)
Set an error, three-string version.
size_t order
Number of points to include in each interpolation (default 3)
void eval(const vec2_t &x, vec3_t &y, std::vector< covar_func_t > &fcovar) const
Perform the interpolation.
data_t vector_max_value(size_t n, const vec_t &data)
Compute the maximum of the first n elements of a vector.
size_t order
Number of points to include in each interpolation (default 3)
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
bool data_set
True if the data has been specified.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
std::vector< vec_t > ptrs_y
A vector of pointers holding the data.
std::vector< ubvector > Kinvf
Inverse covariance matrix times function vector.
size_t np
The number of points.
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.
size_t np
The number of points.
std::vector< vec_t > ptrs_x
A vector of pointers holding the data.
size_t nd_in
The number of dimensions of the inputs.