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> 45 #include <o2scl/cholesky.h> 47 #ifndef DOXYGEN_NO_O2NS 60 template<
class vec_t=boost::numeric::ublas::vector<
double>,
61 class mat_col_t=boost::numeric::ublas::matrix_column<
62 boost::numeric::ublas::matrix<
double> >,
63 class covar_func_t=std::function<
double(const vec_t &,const vec_t &)> >
91 template<
class vec_vec_t,
class vec_vec2_t>
92 void set_data(
size_t n_in,
size_t n_out,
size_t n_points,
93 vec_vec_t &x, vec_vec2_t &y,
94 std::vector<covar_func_t> &fcovar) {
97 O2SCL_ERR2(
"Must provide at least three points in ",
101 O2SCL_ERR2(
"Must provide at least one input column in ",
105 O2SCL_ERR2(
"Must provide at least one output column in ",
112 for(
size_t i=0;i<n_points;i++) {
113 if (x[i].size()!=n_in) {
117 std::swap(
ptrs_x[i],x[i]);
122 std::cout <<
"interpm_krige::set_data() : Using " << n_points
123 <<
" points with " <<
nd_in <<
" input variables and\n\t" 124 <<
nd_out <<
" output variables." << std::endl;
130 for(
size_t iout=0;iout<n_out;iout++) {
133 ubmatrix KXX(n_points,n_points);
134 for(
size_t irow=0;irow<n_points;irow++) {
135 for(
size_t icol=0;icol<n_points;icol++) {
137 KXX(irow,icol)=KXX(icol,irow);
139 KXX(irow,icol)=fcovar[iout](
ptrs_x[irow],ptrs_x[icol]);
145 ubmatrix &inv_KXX=KXX;
146 o2scl_linalg::cholesky_invert<ubmatrix>(n_points,inv_KXX);
149 Kinvf[iout].resize(n_points);
150 boost::numeric::ublas::axpy_prod(inv_KXX,y[iout],Kinvf[iout],
true);
160 template<
class vec2_t,
class vec3_t>
161 void eval(
const vec2_t &x, vec3_t &y,
162 std::vector<covar_func_t> &fcovar)
const {
165 O2SCL_ERR(
"Data not set in interpm_krige::eval().",
170 for(
size_t iout=0;iout<
nd_out;iout++) {
172 for(
size_t ipoints=0;ipoints<
np;ipoints++) {
173 y[iout]+=fcovar[iout](x,
ptrs_x[ipoints])*Kinvf[iout][ipoints];
181 #ifndef DOXYGEN_INTERNAL 210 template<
class vec_t=boost::numeric::ublas::vector<
double>,
211 class mat_col_t=boost::numeric::ublas::matrix_column<
212 boost::numeric::ublas::matrix<
double> >,
213 class covar_func_t=std::function<
double(const vec_t &,const vec_t &)> >
233 template<
class vec_vec_t,
class vec_vec2_t>
234 void set_data(
size_t n_in,
size_t n_out,
size_t n_points,
235 vec_vec_t &x, vec_vec2_t &y,
236 std::vector<covar_func_t> &fcovar,
size_t order) {
239 O2SCL_ERR2(
"Must provide at least three points in ",
243 O2SCL_ERR2(
"Must provide at least one input column in ",
247 O2SCL_ERR2(
"Must provide at least one output column in ",
255 for(
size_t i=0;i<n_points;i++) {
256 if (x[i].size()!=n_in) {
260 std::swap(
ptrs_x[i],x[i]);
262 ptrs_y.resize(n_out);
263 for(
size_t i=0;i<n_out;i++) {
264 if (y[i].size()!=n_points) {
268 std::swap(ptrs_y[i],y[i]);
273 std::cout <<
"interpm_krige_nn::set_data() : Using " << n_points
274 <<
" points with " <<
nd_in <<
" input variables and\n\t" 275 <<
nd_out <<
" output variables and order " 276 << n_order <<
" ." << std::endl;
285 template<
class vec2_t,
class vec3_t>
286 void eval(
const vec2_t &x, vec3_t &y,
287 std::vector<covar_func_t> &fcovar)
const {
290 O2SCL_ERR(
"Data not set in interpm_krige_nn::eval().",
297 for(
size_t iout=0;iout<
nd_out;iout++) {
302 for(
size_t ip=0;ip<
np;ip++) {
303 dists[ip]=-fcovar[iout](x,
ptrs_x[ip]);
308 ubvector_size_t index;
309 o2scl::vector_smallest_index<ubvector,double,ubvector_size_t>
310 (
np,dists,n_order,index);
313 ubvector func(n_order);
314 for(
size_t io=0;io<n_order;io++) {
315 func[io]=ptrs_y[iout][index[io]];
319 ubmatrix KXX(n_order,n_order);
320 for(
size_t irow=0;irow<n_order;irow++) {
321 for(
size_t icol=0;icol<n_order;icol++) {
323 KXX(irow,icol)=KXX(icol,irow);
325 KXX(irow,icol)=fcovar[iout](
ptrs_x[index[irow]],
333 ubmatrix &inv_KXX=KXX;
334 o2scl_linalg::cholesky_invert<ubmatrix>(n_order,inv_KXX);
337 ubvector
Kinvf(n_order);
338 boost::numeric::ublas::axpy_prod(inv_KXX,func,Kinvf,
true);
342 for(
size_t ipoints=0;ipoints<n_order;ipoints++) {
343 y[iout]+=-dists[index[ipoints]]*Kinvf[ipoints];
359 template<
class vec2_t>
361 std::vector<covar_func_t> &fcovar,
362 ubvector_size_t &index,
363 ubvector_size_t &indep)
const {
365 if (x.size()<
nd_in || index.size()<
np || indep.size()<n_order
367 std::cout << x.size() <<
" " <<
nd_in << std::endl;
368 std::cout << index.size() <<
" " <<
np << std::endl;
369 std::cout << indep.size() <<
" " << n_order << std::endl;
370 std::cout << iout <<
" " <<
nd_out << std::endl;
371 O2SCL_ERR(
"Vectors not of correct size in find_lin_indep().",
376 while (done==
false) {
379 ubvector func(n_order);
380 for(
size_t io=0;io<n_order;io++) {
381 func[io]=ptrs_y[iout][index[indep[io]]];
385 ubmatrix KXX(n_order,n_order);
386 for(
size_t irow=0;irow<n_order;irow++) {
387 for(
size_t icol=0;icol<n_order;icol++) {
389 KXX(irow,icol)=KXX(icol,irow);
391 KXX(irow,icol)=fcovar[iout](
ptrs_x[index[indep[irow]]],
392 ptrs_x[index[indep[icol]]]);
403 std::cout <<
"Finding new independent rows." << std::endl;
404 for(
size_t j=0;j<n_order;j++) {
405 std::cout << indep[j] <<
" " << KXX(j,j) << std::endl;
411 std::cout <<
"Max is: " << max << std::endl;
413 for(
size_t j=0;j<n_order;j++) {
417 std::cout <<
"Entry " << j <<
" is zero so replacing " 418 <<
"entry with " << max+1 << std::endl;
423 O2SCL_ERR3(
"Failed to find set of independent points ",
424 "in interpm_krige_nn::find_lin_indep",
425 "(const vec2_t &, size_t).",
440 template<
class vec2_t>
441 double eval(
const vec2_t &x,
size_t iout,
442 std::vector<covar_func_t> &fcovar)
const {
445 O2SCL_ERR(
"Data not set in interpm_krige_nn::eval().",
454 for(
size_t ip=0;ip<
np;ip++) {
455 dists[ip]=-fcovar[iout](x,
ptrs_x[ip]);
458 ubvector_size_t index(np);
459 o2scl::vector_sort_index<ubvector,ubvector_size_t>(
np,dists,index);
464 ubvector_size_t indep(n_order);
465 for(
size_t io=0;io<n_order;io++) {
469 find_lin_indep(x,iout,fcovar,index,indep);
472 ubvector func(n_order);
473 for(
size_t io=0;io<n_order;io++) {
474 func[io]=ptrs_y[iout][index[indep[io]]];
478 ubmatrix KXX(n_order,n_order);
479 for(
size_t irow=0;irow<n_order;irow++) {
480 for(
size_t icol=0;icol<n_order;icol++) {
482 KXX(irow,icol)=KXX(icol,irow);
484 KXX(irow,icol)=fcovar[iout](
ptrs_x[index[indep[irow]]],
485 ptrs_x[index[indep[icol]]]);
492 ubmatrix &inv_KXX=KXX;
493 o2scl_linalg::cholesky_invert<ubmatrix>(n_order,inv_KXX);
496 ubvector
Kinvf(n_order);
497 boost::numeric::ublas::axpy_prod(inv_KXX,func,Kinvf,
true);
501 for(
size_t ipoints=0;ipoints<n_order;ipoints++) {
502 ret+=-dists[index[indep[ipoints]]]*Kinvf[ipoints];
512 template<
class vec2_t>
514 std::vector<covar_func_t> &fcovar)
const {
517 O2SCL_ERR(
"Data not set in interpm_krige_nn::eval_jackknife().",
529 for(
size_t ip=0;ip<
np;ip++) {
530 dists[ip]=-fcovar[iout](x,
ptrs_x[ip]);
534 ubvector_size_t index(np);
535 o2scl::vector_sort_index<ubvector,ubvector_size_t>(
np,dists,index);
540 ubvector_size_t indep(n_order);
541 for(
size_t io=0;io<n_order;io++) {
549 find_lin_indep(x,iout,fcovar,index,indep);
554 for(
size_t jk=0;jk<n_order;jk++) {
557 std::cout <<
"Jackknife: " << jk <<
" matching function value " 558 << ptrs_y[iout][index[jk]] << std::endl;
561 ubvector_size_t indep_jk;
565 ubvector func(n_order-1);
566 for(
size_t io=0;io<n_order-1;io++) {
567 func[io]=ptrs_y[iout][index[indep_jk[io]]];
571 ubmatrix KXX(n_order-1,n_order-1);
572 for(
size_t irow=0;irow<n_order-1;irow++) {
573 for(
size_t icol=0;icol<n_order-1;icol++) {
575 KXX(irow,icol)=KXX(icol,irow);
577 KXX(irow,icol)=fcovar[iout](
ptrs_x[index[indep_jk[irow]]],
578 ptrs_x[index[indep_jk[icol]]]);
585 ubmatrix &inv_KXX=KXX;
586 o2scl_linalg::cholesky_invert<ubmatrix>(n_order-1,inv_KXX);
589 ubvector
Kinvf(n_order-1);
590 boost::numeric::ublas::axpy_prod(inv_KXX,func,Kinvf,
true);
594 for(
size_t ipoints=0;ipoints<n_order-1;ipoints++) {
595 ytmp+=-dists[index[indep_jk[ipoints]]]*Kinvf[ipoints];
599 qual+=pow(ptrs_y[iout][index[jk]]-ytmp,2.0);
602 std::cout <<
"Original value: " 603 << ptrs_y[iout][index[jk]] <<
" interpolated: " 604 << ytmp << std::endl;
614 #ifndef DOXYGEN_INTERNAL 637 #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
Given covariance function fcovar and input vector x store the result of the interpolation in y...
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
Compute a quality factor for interpolation using jackknife resampling.
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
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.
size_t n_order
The order of the interpolation (specified by set_data() )
double eval(const vec2_t &x, size_t iout, std::vector< covar_func_t > &fcovar) const
Given covariance function fcovar and input vector x return the result of the interpolation for functi...
#define O2SCL_ERR3(d, d2, d3, n)
Set an error, three-string version.
void eval(const vec2_t &x, vec3_t &y, std::vector< covar_func_t > &fcovar) const
Given covariance function fcovar and input vector x store the result of the interpolation in y...
data_t vector_max_value(size_t n, const vec_t &data)
Compute the maximum of the first n elements of a vector.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
bool data_set
True if the data has been specified.
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...
#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.
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.