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>
46 #include <o2scl/min_brent_gsl.h>
47 #include <o2scl/cblas.h>
49 #ifndef DOXYGEN_NO_O2NS
70 template<
class vec_t=boost::numeric::ublas::vector<
double>,
71 class mat_t=boost::numeric::ublas::vector<
double>,
72 class mat_row_t=boost::numeric::ublas::matrix_row
73 <boost::numeric::ublas::vector<
double> > >
81 typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
82 typedef boost::numeric::ublas::matrix_row<ubmatrix> ubmatrix_row;
119 template<
class mat2_row_t,
class mat2_t,
class func_vec_t>
121 mat_t &user_x, mat2_t &user_y,
123 const vec_t &noise_var,
bool rescale=
false,
124 bool err_on_fail=
true) {
127 O2SCL_ERR2(
"Must provide at least two points in ",
131 O2SCL_ERR2(
"Must provide at least one input column in ",
135 O2SCL_ERR2(
"Must provide at least one output column in ",
138 if (noise_var.size()<1) {
146 if (user_x.size1()!=n_points || user_x.size2()!=n_in) {
154 if (user_y.size2()!=n_points || user_y.size1()!=n_out) {
164 std::cout <<
"interpm_krige::set_data_noise() : Using " << n_points
165 <<
" points with " <<
nd_in <<
" input variables and\n\t"
166 <<
nd_out <<
" output variables." << std::endl;
172 for(
size_t j=0;j<n_in;j++) {
175 for(
size_t i=1;i<n_points;i++) {
177 if (val>
max[j])
max[j]=val;
178 if (val<
min[j])
min[j]=val;
180 for(
size_t i=0;i<n_points;i++) {
185 std::cout <<
"interpm_krige::set_data_noise() rescaled."
191 size_t n_covar=fcovar.size();
194 for(
size_t iout=0;iout<n_out;iout++) {
196 size_t icovar=iout % n_covar;
197 size_t inoise=iout & noise_var.size();
200 mat2_row_t yiout(user_y,iout);
204 for(
size_t irow=0;irow<n_points;irow++) {
205 mat_row_t xrow(
x,irow);
206 for(
size_t icol=0;icol<n_points;icol++) {
207 mat_row_t xcol(
x,icol);
209 KXX(irow,icol)=KXX(icol,irow);
210 }
else if (irow==icol) {
211 KXX(irow,icol)=fcovar[icovar](xrow,xcol)+
214 KXX(irow,icol)=fcovar[icovar](xrow,xcol);
222 ubmatrix inv_KXX(n_points,n_points);
229 "in interpm_krige::set_data_noise().",
234 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
235 (n_points,KXX,p,inv_KXX);
238 Kinvf[iout].resize(n_points);
240 o2scl_cblas::o2cblas_NoTrans,
241 n_points,n_points,1.0,inv_KXX,
242 yiout,0.0,
Kinvf[iout]);
251 O2SCL_ERR2(
"Matrix singular (Cholesky method) ",
252 "in interpm_krige::set_data_noise().",
258 o2scl_linalg::cholesky_invert<ubmatrix>(n_points,inv_KXX);
261 Kinvf[iout].resize(n_points);
263 o2scl_cblas::o2cblas_NoTrans,
264 n_points,n_points,1.0,inv_KXX,
265 yiout,0.0,
Kinvf[iout]);
270 std::cout <<
"interpm_krige::set_data_noise() finished " << iout+1
271 <<
" of " << n_out <<
"." << std::endl;
277 std::cout <<
"interpm_krige::set_data_noise() done."
290 template<
class mat2_row_t,
class mat2_t,
class func_vec_t>
291 int set_data(
size_t n_in,
size_t n_out,
size_t n_points,
292 mat_t &user_x, mat2_t &user_y,
293 func_vec_t &fcovar,
bool rescale=
false,
294 bool err_on_fail=
true) {
298 return set_data_noise<mat2_row_t,mat2_t,func_vec_t>
299 (n_in,n_out,n_points,user_x,user_y,fcovar,
300 noise_vec,rescale,err_on_fail);
306 template<
class vec2_t,
class vec3_t,
class vec_func_t>
307 void eval(
const vec2_t &x0, vec3_t &y0, vec_func_t &fcovar) {
310 O2SCL_ERR(
"Data not set in interpm_krige::eval().",
313 if (fcovar.size()==0) {
314 O2SCL_ERR(
"No covariance functions in interpm_krige::eval().",
323 for(
size_t iin=0;iin<
nd_in;iin++) {
324 x0p[iin]=((x0[iin]-
min[iin])/(
max[iin]-
min[iin])-0.5)*2.0;
328 for(
size_t iout=0;iout<
nd_out;iout++) {
329 size_t icovar=iout % fcovar.size();
331 for(
size_t ipoints=0;ipoints<
np;ipoints++) {
332 mat_row_t xrow(
x,ipoints);
333 double covar_val=fcovar[icovar](xrow,x0p);
334 y0[iout]+=covar_val*
Kinvf[iout][ipoints];
341 for(
size_t iout=0;iout<
nd_out;iout++) {
342 size_t icovar=iout % fcovar.size();
344 for(
size_t ipoints=0;ipoints<
np;ipoints++) {
345 mat_row_t xrow(
x,ipoints);
346 y0[iout]+=fcovar[icovar](xrow,x0)*
Kinvf[iout][ipoints];
356 #ifndef DOXYGEN_INTERNAL
388 template<
class vec_t=boost::numeric::ublas::vector<
double>,
389 class mat_t=boost::numeric::ublas::vector<
double>,
390 class mat_row_t=boost::numeric::ublas::matrix_row
391 <boost::numeric::ublas::vector<
double> > >
399 typedef boost::numeric::ublas::matrix_column<ubmatrix> ubmatrix_column;
402 std::vector<std::function<double(
const mat_row_t &,
const vec_t &)> >
408 std::vector<std::function<double(
const mat_row_t &,
const mat_row_t &)> >
427 template<
class vec2_t,
class vec3_t>
428 double covar(
const vec2_t &
x1,
const vec3_t &
x2,
size_t sz,
double len2) {
430 for(
size_t i=0;i<sz;i++) {
431 ret+=pow(
x1[i]-
x2[i],2.0);
433 ret=exp(-ret/len2/len2/2.0);
442 template<
class vec3_t>
443 double qual_fun(
double xlen,
double noise_var,
size_t iout,
450 size_t size=this->
x.size1();
454 for(
size_t ell=0;ell<
loo_npts;ell++) {
466 for(
size_t irow=0;irow<size-1;irow++) {
468 for(
size_t icol=0;icol<size-1;icol++) {
471 KXX(irow,icol)=KXX(icol,irow);
473 KXX(irow,icol)=
covar<
476 (xrow,xcol,this->
nd_in,xlen);
477 if (irow==icol) KXX(irow,icol)+=noise_var;
486 std::cout <<
"Performing Cholesky decomposition with size "
487 << size-1 << std::endl;
495 std::cout <<
"Performing matrix inversion with size "
496 << size-1 << std::endl;
499 o2scl_linalg::cholesky_invert<ubmatrix>(size-1,KXX);
502 this->
Kinvf[iout].resize(size-1);
504 o2scl_cblas::o2cblas_NoTrans,
505 size-1,size-1,1.0,KXX,
506 y,0.0,this->
Kinvf[iout]);
515 std::cout <<
"Performing LU decomposition with size "
516 << size-1 << std::endl;
524 std::cout <<
"Performing matrix inversion with size "
525 << size-1 << std::endl;
527 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
528 (size-1,KXX,p,inv_KXX);
531 this->
Kinvf[iout].resize(size-1);
533 o2scl_cblas::o2cblas_NoTrans,
534 size-1,size-1,1.0,inv_KXX,
535 y,0.0,this->
Kinvf[iout]);
541 for(
size_t i=0;i<size-1;i++) {
543 mat_row_t xcol(this->
x,row);
544 ypred+=covar<matrix_row_gen<matrix_view_omit_row<mat_t> >,
545 mat_row_t>(xrow,xcol,this->
nd_in,xlen)*
546 this->
Kinvf[iout][i];
549 std::cout <<
"act,pred: " << yact <<
" " << ypred << std::endl;
552 ret+=pow(yact-ypred,2.0);
557 std::cout <<
"ret: " << ret << std::endl;
562 std::cout <<
"Creating covariance matrix with size "
563 << size << std::endl;
568 for(
size_t irow=0;irow<size;irow++) {
569 mat_row_t xrow(this->
x,irow);
570 for(
size_t icol=0;icol<size;icol++) {
571 mat_row_t xcol(this->
x,icol);
573 KXX(irow,icol)=KXX(icol,irow);
575 KXX(irow,icol)=covar<mat_row_t,mat_row_t>(xrow,xcol,
577 if (irow==icol) KXX(irow,icol)+=noise_var;
590 std::cout <<
"Performing Cholesky decomposition with size "
591 << size << std::endl;
599 std::cout <<
"Performing matrix inversion with size "
600 << size << std::endl;
602 o2scl_linalg::cholesky_invert<ubmatrix>(size,KXX);
604 lndet=o2scl_linalg::cholesky_lndet<ubmatrix>(size,KXX);
607 this->
Kinvf[iout].resize(size);
609 o2scl_cblas::o2cblas_NoTrans,
611 y,0.0,this->
Kinvf[iout]);
620 std::cout <<
"Performing LU decomposition with size "
621 << size << std::endl;
629 std::cout <<
"Performing matrix inversion with size "
630 << size << std::endl;
632 o2scl_linalg::LU_invert<ubmatrix,ubmatrix,ubmatrix_column>
633 (size,KXX,p,inv_KXX);
635 lndet=o2scl_linalg::LU_lndet<ubmatrix>(size,KXX);
638 this->
Kinvf[iout].resize(size);
640 o2scl_cblas::o2cblas_NoTrans,
641 size,size,1.0,inv_KXX,
642 y,0.0,this->
Kinvf[iout]);
649 for(
size_t i=0;i<size;i++) {
650 ret+=0.5*y[i]*this->
Kinvf[iout][i];
710 template<
class mat2_row_t,
class mat2_t,
class vec2_t,
713 mat_t &user_x, mat2_t &user_y,
714 const vec2_t &noise_var,
const vec3_t &len_precompute,
715 bool rescale=
false,
bool err_on_fail=
true) {
718 O2SCL_ERR2(
"Must provide at least two points in ",
719 "interpm_krige_optim::set_data_noise()",
exc_efailed);
722 O2SCL_ERR2(
"Must provide at least one input column in ",
723 "interpm_krige_optim::set_data_noise()",
exc_efailed);
726 O2SCL_ERR2(
"Must provide at least one output column in ",
727 "interpm_krige_optim::set_data_noise()",
exc_efailed);
729 if (noise_var.size()<1) {
738 std::swap(this->
x,user_x);
743 std::cout <<
"interpm_krige_optim::set_data_noise() : Using "
745 <<
" points with\n " << n_in <<
" input variables and "
746 << this->
nd_out <<
" output variables." << std::endl;
750 this->
min.resize(n_in);
751 this->
max.resize(n_in);
752 for(
size_t j=0;j<n_in;j++) {
753 this->
min[j]=(this->
x)(0,j);
754 this->
max[j]=(this->
x)(0,j);
755 for(
size_t i=1;i<n_points;i++) {
756 double val=(this->
x)(i,j);
757 if (val>this->
max[j]) this->
max[j]=val;
758 if (val<this->
min[j]) this->
min[j]=val;
763 for(
size_t j=0;j<n_in;j++) {
764 for(
size_t i=0;i<n_points;i++) {
765 (this->
x)(i,j)=(((this->
x)(i,j)-this->
min[j])/
766 (this->
max[j]-this->
min[j])-0.5)*2.0;
770 std::cout <<
"interpm_krige_optim::set_data_noise() rescaled."
777 this->
Kinvf.resize(n_out);
784 for(
size_t iout=0;iout<n_out;iout++) {
787 mat2_row_t yiout(user_y,iout);
789 if (iout<len_precompute.size()) {
792 std::cout <<
"interp_krige_optim: precomputed length "
793 << len_precompute[iout] << std::endl;
795 len[iout]=len_precompute[iout];
796 size_t mode_temp=
mode;
804 std::cout <<
"interp_krige_optim: full minimization"
811 len_opt=(this->
max[iout]-this->
min[iout])/((
double)this->
np);
813 len_opt=2.0/((double)this->
np);
817 (std::mem_fn<
double(
double,
double,
size_t,mat2_row_t &,
int &)>
819 mat_row_t>::qual_fun<mat2_row_t>),
820 this,std::placeholders::_1,noise_var[iout],iout,yiout,
829 "interp_krige_optim::set_noise().",
837 std::cout <<
"interp_krige_optim::set_data_noise() : "
838 <<
"simple minimization" << std::endl;
845 len_avg=(this->
max[iout]-this->
min[iout])/((
double)this->
np);
847 len_avg=2.0/((double)this->
np);
850 len_min=len_avg/1.0e2;
856 std::cout <<
"iout, len (min,max,step): " << iout
858 << pow(len_ratio,((
double)1)/((
double)
nlen-1))
863 double min_qual=0.0, len_opt=0.0;
866 std::cout <<
"ilen len qual fail min_qual len_opt" << std::endl;
871 for(
size_t j=0;j<
nlen;j++) {
872 double xlen=
len_min*pow(len_ratio,((
double)j)/((
double)
nlen-1));
877 if (
success==0 && (min_set==
false ||
qual[iout]<min_qual)) {
884 std::cout <<
"interp_krige_optim: ";
886 std::cout << j <<
" " << xlen <<
" " <<
qual[iout] <<
" "
887 <<
success <<
" " << min_qual <<
" "
888 << len_opt << std::endl;
894 std::cout <<
"interp_krige_optim: ";
896 std::cout <<
" " << len_opt <<
" " << min_qual << std::endl;
905 ff1[iout]=std::bind(std::mem_fn<
double(
const mat_row_t &,
909 mat_row_t>::
covar<mat_row_t,
911 std::placeholders::_1,std::placeholders::_2,
913 ff2[iout]=std::bind(std::mem_fn<
double(
const mat_row_t &,
917 mat_row_t>::
covar<mat_row_t,
919 std::placeholders::_1,std::placeholders::_2,
932 template<
class mat2_row_t,
class mat2_t,
class vec2_t>
933 int set_data(
size_t n_in,
size_t n_out,
size_t n_points,
934 mat_t &user_x, mat2_t &user_y,
935 const vec2_t &len_precompute,
936 bool rescale=
false,
bool err_on_fail=
true) {
940 return set_data_noise<mat2_row_t,mat2_t,vec_t,vec2_t>
941 (n_in,n_out,n_points,user_x,
942 user_y,noise_vec,len_precompute,rescale,err_on_fail);
957 template<
class vec_t=boost::numeric::ublas::vector<
double>,
958 class mat_t=boost::numeric::ublas::vector<
double> >
978 template<
class mat2_t,
class vec_func_t>
979 void set_data(
size_t n_in,
size_t n_out,
size_t n_points,
980 mat_t &user_x, mat_t &user_y,
981 vec_func_t &fcovar,
size_t order) {
984 O2SCL_ERR2(
"Must provide at least two points in ",
988 O2SCL_ERR2(
"Must provide at least one input column in ",
992 O2SCL_ERR2(
"Must provide at least one output column in ",
1000 for(
size_t i=0;i<n_points;i++) {
1001 if (
x[i].size()!=n_in) {
1005 std::swap(
x[i],
x[i]);
1008 for(
size_t i=0;i<n_out;i++) {
1009 if (
y[i].size()!=n_points) {
1013 std::swap(
y[i],
y[i]);
1018 std::cout <<
"interpm_krige_nn::set_data() : Using " << n_points
1019 <<
" points with " <<
nd_in <<
" input variables and\n\t"
1020 <<
nd_out <<
" output variables and order "
1021 <<
n_order <<
" ." << std::endl;
1030 template<
class vec2_t,
class vec3_t,
class vec_func_t>
1031 void eval(
const vec2_t &x0, vec3_t &y0,
1032 vec_func_t &fcovar)
const {
1035 O2SCL_ERR(
"Data not set in interpm_krige_nn::eval().",
1042 for(
size_t iout=0;iout<
nd_out;iout++) {
1047 for(
size_t ip=0;ip<
np;ip++) {
1048 dists[ip]=-fcovar[iout](x0,
x[ip]);
1054 o2scl::vector_smallest_index<ubvector,double,ubvector_size_t>
1059 for(
size_t io=0;io<
n_order;io++) {
1060 func[io]=
y[iout][index[io]];
1065 for(
size_t irow=0;irow<
n_order;irow++) {
1066 for(
size_t icol=0;icol<
n_order;icol++) {
1068 KXX(irow,icol)=KXX(icol,irow);
1070 KXX(irow,icol)=fcovar[iout](
x[index[irow]],
1079 o2scl_linalg::cholesky_invert<ubmatrix>(
n_order,inv_KXX);
1083 boost::numeric::ublas::axpy_prod(inv_KXX,func,
Kinvf,
true);
1087 for(
size_t ipoints=0;ipoints<
n_order;ipoints++) {
1088 y0[iout]+=-dists[index[ipoints]]*
Kinvf[ipoints];
1104 template<
class vec2_t,
class vec_func_t>
1112 std::cout <<
x2.size() <<
" " <<
nd_in << std::endl;
1113 std::cout << index.size() <<
" " <<
np << std::endl;
1114 std::cout << indep.size() <<
" " <<
n_order << std::endl;
1115 std::cout << iout <<
" " <<
nd_out << std::endl;
1116 O2SCL_ERR(
"Vectors not of correct size in find_lin_indep().",
1121 while (done==
false) {
1125 for(
size_t io=0;io<
n_order;io++) {
1126 func[io]=
y[iout][index[indep[io]]];
1131 for(
size_t irow=0;irow<
n_order;irow++) {
1132 for(
size_t icol=0;icol<
n_order;icol++) {
1134 KXX(irow,icol)=KXX(icol,irow);
1136 KXX(irow,icol)=fcovar[iout](
x2[index[indep[irow]]],
1137 x2[index[indep[icol]]]);
1148 std::cout <<
"Finding new independent rows." << std::endl;
1149 for(
size_t j=0;j<
n_order;j++) {
1150 std::cout << indep[j] <<
" " << KXX(j,j) << std::endl;
1156 std::cout <<
"Max is: " <<
max << std::endl;
1158 for(
size_t j=0;j<
n_order;j++) {
1159 if (KXX(j,j)==0.0) {
1162 std::cout <<
"Entry " << j <<
" is zero so replacing "
1163 <<
"entry with " <<
max+1 << std::endl;
1168 O2SCL_ERR3(
"Failed to find set of independent points ",
1169 "in interpm_krige_nn::find_lin_indep",
1170 "(const vec2_t &, size_t).",
1185 template<
class vec2_t,
class vec_func_t>
1187 vec_func_t &fcovar)
const {
1190 O2SCL_ERR(
"Data not set in interpm_krige_nn::eval().",
1199 for(
size_t ip=0;ip<
np;ip++) {
1200 dists[ip]=-fcovar[iout](
x2,
x2[ip]);
1204 o2scl::vector_sort_index<ubvector,ubvector_size_t>(
np,dists,index);
1210 for(
size_t io=0;io<
n_order;io++) {
1218 for(
size_t io=0;io<
n_order;io++) {
1219 func[io]=
y[iout][index[indep[io]]];
1224 for(
size_t irow=0;irow<
n_order;irow++) {
1225 for(
size_t icol=0;icol<
n_order;icol++) {
1227 KXX(irow,icol)=KXX(icol,irow);
1229 KXX(irow,icol)=fcovar[iout](
x2[index[indep[irow]]],
1230 x2[index[indep[icol]]]);
1238 o2scl_linalg::cholesky_invert<ubmatrix>(
n_order,inv_KXX);
1242 boost::numeric::ublas::axpy_prod(inv_KXX,func,
Kinvf,
true);
1246 for(
size_t ipoints=0;ipoints<
n_order;ipoints++) {
1247 ret+=-dists[index[indep[ipoints]]]*
Kinvf[ipoints];
1257 template<
class vec2_t,
class func_vec_t>
1259 func_vec_t &fcovar)
const {
1262 O2SCL_ERR(
"Data not set in interpm_krige_nn::eval_jackknife().",
1274 for(
size_t ip=0;ip<
np;ip++) {
1275 dists[ip]=-fcovar[iout](
x2,
x2[ip]);
1280 o2scl::vector_sort_index<ubvector,ubvector_size_t>(
np,dists,index);
1286 for(
size_t io=0;io<
n_order;io++) {
1299 for(
size_t jk=0;jk<
n_order;jk++) {
1302 std::cout <<
"Jackknife: " << jk <<
" matching function value "
1303 <<
y[iout][index[jk]] << std::endl;
1311 for(
size_t io=0;io<
n_order-1;io++) {
1312 func[io]=
y[iout][index[indep_jk[io]]];
1317 for(
size_t irow=0;irow<
n_order-1;irow++) {
1318 for(
size_t icol=0;icol<
n_order-1;icol++) {
1320 KXX(irow,icol)=KXX(icol,irow);
1322 KXX(irow,icol)=fcovar[iout](
x2[index[indep_jk[irow]]],
1323 x2[index[indep_jk[icol]]]);
1331 o2scl_linalg::cholesky_invert<ubmatrix>(
n_order-1,inv_KXX);
1335 boost::numeric::ublas::axpy_prod(inv_KXX,func,
Kinvf,
true);
1339 for(
size_t ipoints=0;ipoints<
n_order-1;ipoints++) {
1340 ytmp+=-dists[index[indep_jk[ipoints]]]*
Kinvf[ipoints];
1344 qual+=pow(
y[iout][index[jk]]-ytmp,2.0);
1347 std::cout <<
"Original value: "
1348 <<
y[iout][index[jk]] <<
" interpolated: "
1349 << ytmp << std::endl;
1359 #ifndef DOXYGEN_INTERNAL
1372 std::vector<vec_t>
x;
1374 std::vector<vec_t>
y;
1382 #ifndef DOXYGEN_NO_O2NS