23 #ifndef O2SCL_INTERPM_IDW_H
24 #define O2SCL_INTERPM_IDW_H
34 #include <boost/numeric/ublas/matrix.hpp>
36 #include <gsl/gsl_combination.h>
38 #include <o2scl/err_hnd.h>
39 #include <o2scl/vector.h>
40 #include <o2scl/vec_stats.h>
41 #include <o2scl/linear_solver.h>
42 #include <o2scl/columnify.h>
43 #include <o2scl/table.h>
45 #ifndef DOXYGEN_NO_O2NS
112 template<
class mat_t=const_matrix_view_table<> >
159 O2SCL_ERR(
"Points cannot be zero in interpm_idw.",
171 template<
class vec2_t>
void set_scales(
size_t n, vec2_t &v) {
173 O2SCL_ERR(
"Scale vector size cannot be zero in interpm_idw.",
176 for(
size_t i=0;i<n;i++) {
178 O2SCL_ERR(
"Scale must be positive and non-zero in interpm_idw.",
195 void set_data(
size_t n_in,
size_t n_out,
size_t n_points,
196 mat_t &dat,
bool auto_scale_flag=
true) {
203 O2SCL_ERR2(
"Must provide at least one input column in ",
207 O2SCL_ERR2(
"Must provide at least one output column in ",
216 if (auto_scale_flag) {
225 void get_data(
size_t &n_in,
size_t &n_out,
size_t &n_points,
243 for(
size_t i=0;i<
nd_in;i++) {
244 double min=
data(i,0), max=min;
245 for(
size_t j=1;j<
np;j++) {
246 double val=
data(i,j);
247 if (val>max) max=val;
248 if (val<min) min=val;
274 template<
class vec2_t>
double operator()(
const vec2_t &x)
const {
280 template<
class vec2_t>
double eval(
const vec2_t &x)
const {
283 O2SCL_ERR(
"Data not set in interpm_idw::eval().",
288 std::vector<double> dists(
np);
289 for(
size_t i=0;i<
np;i++) {
294 std::vector<size_t> index;
295 o2scl::vector_smallest_index<std::vector<double>,double,
301 while (found==
true) {
306 double dist_jk=
dist(j,k);
309 index.erase(index.begin()+j);
317 if (dists[index[0]]<=0.0) {
323 for(
size_t i=0;i<
points;i++) {
324 norm+=1.0/dists[index[i]];
329 for(
size_t i=0;i<
points;i++) {
330 ret+=
data(
nd_in,index[i])/dists[index[i]];
341 template<
class vec2_t>
void eval_err(
const vec2_t &x,
double &val,
345 O2SCL_ERR(
"Data not set in interpm_idw::eval_err().",
350 std::vector<double> dists(
np);
351 for(
size_t i=0;i<
np;i++) {
356 std::vector<size_t> index;
357 o2scl::vector_smallest_index<std::vector<double>,double,
363 while (found==
true) {
368 double dist_jk=
dist(j,k);
371 index.erase(index.begin()+j);
378 if (dists[index[0]]<=0.0) {
387 std::vector<double> vals(
points+1);
389 for(
size_t j=0;j<
points+1;j++) {
393 for(
size_t i=0;i<
points+1;i++) {
394 if (i!=j) norm+=1.0/dists[index[i]];
399 for(
size_t i=0;i<
points+1;i++) {
401 vals[j]+=
data(
nd_in,index[i])/dists[index[i]];
419 template<
class vec2_t,
class vec3_t>
420 void eval(
const vec2_t &x, vec3_t &y)
const {
423 O2SCL_ERR(
"Data not set in interpm_idw::eval().",
428 std::cout <<
"interpm_idw: input: ";
429 for(
size_t k=0;k<
nd_in;k++) {
430 std::cout << x[k] <<
" ";
432 std::cout << std::endl;
436 std::vector<double> dists(
np);
437 for(
size_t i=0;i<
np;i++) {
442 std::vector<size_t> index;
443 o2scl::vector_smallest_index<std::vector<double>,double,
444 std::vector<size_t> >(dists,
points,index);
446 for(
size_t i=0;i<
points;i++) {
447 std::cout <<
"interpm_idw: closest point: ";
448 for(
size_t k=0;k<
nd_in;k++) {
449 std::cout <<
data(k,index[i]) <<
" ";
451 std::cout << std::endl;
458 while (found==
true) {
463 double dist_jk=
dist(j,k);
466 index.erase(index.begin()+j);
475 if (dists[index[0]]<=0.0) {
476 for(
size_t i=0;i<
nd_out;i++) {
480 std::cout <<
"interpm_idw: distance zero. "
481 <<
"Returning values at index: " << index[0] << std::endl;
490 for(
size_t i=0;i<
points;i++) {
491 norm+=1.0/dists[index[i]];
494 std::cout <<
"interpm_idw: norm is " << norm << std::endl;
498 for(
size_t j=0;j<
nd_out;j++) {
500 for(
size_t i=0;i<
points;i++) {
502 std::cout <<
"interpm_idw: Point: ";
503 for(
size_t k=0;k<
nd_in;k++) {
504 std::cout <<
data(k,index[i]) <<
" ";
506 std::cout << std::endl;
508 y[j]+=
data(
nd_in+j,index[i])/dists[index[i]];
510 std::cout <<
"interpm_idw: j,points,value,1/dist: "
511 << j <<
" " << i <<
" "
513 << 1.0/dists[index[i]] << std::endl;
518 std::cout <<
"interpm_idw: y[" << j <<
"]: " << y[j]
533 template<
class vec2_t,
class vec3_t,
class vec4_t>
535 std::vector<size_t> &index)
const {
538 O2SCL_ERR(
"Data not set in interpm_idw::eval_err().",
543 std::vector<double> dists(
np);
544 for(
size_t i=0;i<
np;i++) {
550 o2scl::vector_smallest_index<std::vector<double>,double,
556 while (found==
true) {
561 double dist_jk=
dist(j,k);
564 index.erase(index.begin()+j);
571 if (dists[index[0]]<=0.0) {
575 for(
size_t k=0;k<
nd_out;k++) {
583 for(
size_t k=0;k<
nd_out;k++) {
585 std::vector<double> vals(
points+1);
587 for(
size_t j=0;j<
points+1;j++) {
591 for(
size_t i=0;i<
points+1;i++) {
592 if (i!=j) norm+=1.0/dists[index[i]];
597 for(
size_t i=0;i<
points+1;i++) {
599 vals[j]+=
data(
nd_in+k,index[i])/dists[index[i]];
623 template<
class vec2_t,
class vec3_t,
class vec4_t>
624 void eval_err(
const vec2_t &x, vec3_t &val, vec4_t &err)
const {
625 std::vector<size_t> index;
647 template<
class vec3_t>
649 vec3_t &derivs, vec3_t &errs)
const {
652 O2SCL_ERR(
"Data not set in interpm_idw::derivs_err().",
658 for(
size_t i=0;i<
nd_in;i++) {
659 x[i]=
data(i,point_index);
662 double f=
data(
nd_in+func_index,point_index);
668 std::vector<double> dists(
np);
669 for(
size_t i=0;i<
np;i++) {
675 std::vector<size_t> index;
676 size_t max_smallest=(
nd_in+2)*2;
677 if (max_smallest>
np) max_smallest=
np;
678 if (max_smallest<
nd_in+1) {
683 std::cout <<
"max_smallest: " << max_smallest << std::endl;
686 o2scl::vector_smallest_index<std::vector<double>,double,
687 std::vector<size_t> >(dists,max_smallest,index);
690 for(
size_t i=0;i<index.size();i++) {
691 std::cout <<
"index[" << i <<
"] = " << index[i] <<
" "
692 << dists[index[i]] << std::endl;
696 std::vector<size_t> index2;
697 for(
size_t i=0;i<max_smallest;i++) {
698 if (dists[index[i]]>0.0) {
699 index2.push_back(index[i]);
700 if (index2.size()==
nd_in+1) i=max_smallest;
703 if (index2.size()<
nd_in+1) {
704 O2SCL_ERR(
"Couldn't find enough nearby points (2).",
709 for(
size_t i=0;i<index2.size();i++) {
710 std::cout <<
"index2[" << i <<
"] = " << index2[i] <<
" "
711 << dists[index2[i]] << std::endl;
716 std::vector<ubvector> units(
nd_in+1);
718 std::vector<double> diff_norms(
nd_in+1);
720 std::vector<ubvector> ders(
nd_in+1);
726 std::vector<ubvector> ders2(
nd_in);
728 for(
size_t i=0;i<
nd_in+1;i++) {
731 units[i].resize(
nd_in);
732 for(
size_t j=0;j<
nd_in;j++) {
733 units[i][j]=
data(j,index2[i])-x[j];
737 diff_norms[i]=o2scl::vector_norm<ubvector,double>(units[i]);
738 for(
size_t j=0;j<
nd_in;j++) {
739 units[i][j]/=diff_norms[i];
746 std::cout <<
"Point: ";
747 for(
size_t i=0;i<
nd_in;i++) {
748 std::cout << x[i] <<
" ";
750 std::cout << f << std::endl;
751 for(
size_t j=0;j<
nd_in+1;j++) {
752 std::cout <<
"Closest: " << j <<
" " << index2[j] <<
" ";
753 for(
size_t i=0;i<
nd_in;i++) {
754 std::cout <<
data(i,index2[j]) <<
" ";
756 std::cout <<
data(
nd_in+func_index,index2[j]) <<
" "
757 << diff_norms[j] << std::endl;
759 for(
size_t j=0;j<
nd_in+1;j++) {
760 std::cout <<
"diff_norm: " << j <<
" " << diff_norms[j]
767 for(
size_t i=0;i<
nd_in+1;i++) {
769 ders[i].resize(
nd_in);
773 for(
size_t j=0;j<
nd_in+1;j++) {
775 for(
size_t k=0;k<
nd_in;k++) {
778 v[jj]=(
data(
nd_in+func_index,index2[j])-f)/diff_norms[j];
785 std::cout <<
"m:" << std::endl;
787 std::cout <<
"v:" << std::endl;
792 std::cout <<
"Derivs: " << i <<
" ";
793 std::cout.setf(std::ios::showpos);
794 for(
size_t j=0;j<
nd_in;j++) {
795 std::cout << ders[i][j] <<
" ";
797 std::cout.unsetf(std::ios::showpos);
798 std::cout << std::endl;
804 for(
size_t i=0;i<
nd_in;i++) {
807 ders2[i].resize(
nd_in+1);
808 for(
size_t j=0;j<
nd_in+1;j++) {
809 ders2[i][j]=ders[j][i];
821 #ifndef DOXYGEN_INTERNAL
846 template<
class vec2_t>
double dist(
size_t index,
847 const vec2_t &x)
const {
849 size_t nscales=
scales.size();
850 for(
size_t i=0;i<
nd_in;i++) {
859 double dist(
size_t j,
size_t k)
const {
861 size_t nscales=
scales.size();
862 for(
size_t i=0;i<
nd_in;i++) {
873 #ifndef DOXYGEN_NO_O2NS