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> 44 #ifndef DOXYGEN_NO_O2NS 112 template<
class vec_t=boost::numeric::ublas::vector<
double> >
154 O2SCL_ERR(
"Points cannot be zero in interpm_idw.",
166 template<
class vec2_t>
void set_scales(
size_t n, vec2_t &v) {
168 O2SCL_ERR(
"Scale vector size cannot be zero in interpm_idw.",
171 for(
size_t i=0;i<n;i++) {
173 O2SCL_ERR(
"Scale must be positive and non-zero in interpm_idw.",
189 template<
class vec_vec_t>
190 void set_data(
size_t n_in,
size_t n_out,
size_t n_points,
191 vec_vec_t &vecs,
bool auto_scale_flag=
true) {
194 O2SCL_ERR2(
"Must provide at least three points in ",
198 O2SCL_ERR2(
"Must provide at least one input column in ",
202 O2SCL_ERR2(
"Must provide at least one output column in ",
208 ptrs.resize(n_in+n_out);
209 for(
size_t i=0;i<n_in+n_out;i++) {
210 std::swap(
ptrs[i],vecs[i]);
214 if (auto_scale_flag) {
223 template<
class vec_vec_t>
224 void get_data(
size_t &n_in,
size_t &n_out,
size_t &n_points,
229 vecs.resize(n_in+n_out);
230 for(
size_t i=0;i<n_in+n_out;i++) {
231 std::swap(
ptrs[i],vecs[i]);
246 for(
size_t i=0;i<
nd_in;i++) {
247 scales[i]=fabs(o2scl::vector_max_value<vec_t,double>(
np,
ptrs[i])-
248 o2scl::vector_min_value<vec_t,double>(
np,
ptrs[i]));
261 template<
class vec_vec_t>
273 template<
class vec2_t>
double operator()(
const vec2_t &x)
const {
279 template<
class vec2_t>
double eval(
const vec2_t &x)
const {
282 O2SCL_ERR(
"Data not set in interpm_idw::eval().",
287 std::vector<double> dists(
np);
288 for(
size_t i=0;i<
np;i++) {
293 std::vector<size_t> index;
294 o2scl::vector_smallest_index<std::vector<double>,double,
300 while (found==
true) {
305 double dist_jk=
dist(j,k);
308 index.erase(index.begin()+j);
316 if (dists[index[0]]<=0.0) {
322 for(
size_t i=0;i<
points;i++) {
323 norm+=1.0/dists[index[i]];
328 for(
size_t i=0;i<
points;i++) {
329 ret+=
ptrs[
nd_in][index[i]]/dists[index[i]];
340 template<
class vec2_t>
void eval_err(
const vec2_t &x,
double &val,
344 O2SCL_ERR(
"Data not set in interpm_idw::eval_err().",
349 std::vector<double> dists(
np);
350 for(
size_t i=0;i<
np;i++) {
355 std::vector<size_t> index;
356 o2scl::vector_smallest_index<std::vector<double>,double,
362 while (found==
true) {
367 double dist_jk=
dist(j,k);
370 index.erase(index.begin()+j);
377 if (dists[index[0]]<=0.0) {
386 std::vector<double> vals(
points+1);
388 for(
size_t j=0;j<
points+1;j++) {
392 for(
size_t i=0;i<points+1;i++) {
393 if (i!=j) norm+=1.0/dists[index[i]];
398 for(
size_t i=0;i<points+1;i++) {
400 vals[j]+=
ptrs[
nd_in][index[i]]/dists[index[i]];
418 template<
class vec2_t,
class vec3_t>
419 void eval(vec2_t &x, vec3_t &y)
const {
422 O2SCL_ERR(
"Data not set in interpm_idw::eval().",
427 std::cout <<
"interpm_idw: input: ";
428 for(
size_t k=0;k<
nd_in;k++) {
429 std::cout << x[k] <<
" ";
431 std::cout << std::endl;
435 std::vector<double> dists(
np);
436 for(
size_t i=0;i<
np;i++) {
441 std::vector<size_t> index;
442 o2scl::vector_smallest_index<std::vector<double>,double,
443 std::vector<size_t> >(dists,
points,index);
445 for(
size_t i=0;i<
points;i++) {
446 std::cout <<
"interpm_idw: closest point: ";
447 for(
size_t k=0;k<
nd_in;k++) {
448 std::cout <<
ptrs[k][index[i]] <<
" ";
450 std::cout << std::endl;
457 while (found==
true) {
460 for(
size_t j=0;j<points+
n_extra;j++) {
461 for(
size_t k=j;k<points+
n_extra;k++) {
462 double dist_jk=
dist(j,k);
463 if (index.size()>points && dist_jk<
min_dist) {
465 index.erase(index.begin()+j);
474 if (dists[index[0]]<=0.0) {
475 for(
size_t i=0;i<
nd_out;i++) {
479 std::cout <<
"interpm_idw: distance zero. " 480 <<
"Returning values at index: " << index[0] << std::endl;
489 for(
size_t i=0;i<
points;i++) {
490 norm+=1.0/dists[index[i]];
493 std::cout <<
"interpm_idw: norm is " << norm << std::endl;
497 for(
size_t j=0;j<
nd_out;j++) {
499 for(
size_t i=0;i<
points;i++) {
500 if (j==0 && verbose>0) {
501 std::cout <<
"interpm_idw: Point: ";
502 for(
size_t k=0;k<
nd_in;k++) {
503 std::cout <<
ptrs[k][index[i]] <<
" ";
505 std::cout << std::endl;
507 y[j]+=
ptrs[
nd_in+j][index[i]]/dists[index[i]];
509 std::cout <<
"interpm_idw: j,points,value,1/dist: " 510 << j <<
" " << i <<
" " 512 << 1.0/dists[index[i]] << std::endl;
517 std::cout <<
"interpm_idw: y[" << j <<
"]: " << y[j]
534 template<
class vec2_t,
class vec3_t,
class vec4_t>
536 std::vector<size_t> &index)
const {
539 O2SCL_ERR(
"Data not set in interpm_idw::eval_err().",
544 std::vector<double> dists(
np);
545 for(
size_t i=0;i<
np;i++) {
551 o2scl::vector_smallest_index<std::vector<double>,double,
557 while (found==
true) {
562 double dist_jk=
dist(j,k);
565 index.erase(index.begin()+j);
572 if (dists[index[0]]<=0.0) {
576 for(
size_t k=0;k<
nd_out;k++) {
584 for(
size_t k=0;k<
nd_out;k++) {
586 std::vector<double> vals(
points+1);
588 for(
size_t j=0;j<
points+1;j++) {
592 for(
size_t i=0;i<points+1;i++) {
593 if (i!=j) norm+=1.0/dists[index[i]];
598 for(
size_t i=0;i<points+1;i++) {
600 vals[j]+=
ptrs[
nd_in+k][index[i]]/dists[index[i]];
624 template<
class vec2_t,
class vec3_t,
class vec4_t>
625 void eval_err(
const vec2_t &x, vec3_t &val, vec4_t &err)
const {
626 std::vector<size_t> index;
648 template<
class vec3_t>
650 vec3_t &derivs, vec3_t &errs)
const {
653 O2SCL_ERR(
"Data not set in interpm_idw::derivs_err().",
659 for(
size_t i=0;i<
nd_in;i++) {
660 x[i]=
ptrs[i][point_index];
663 double f=
ptrs[nd_in+func_index][point_index];
669 std::vector<double> dists(
np);
670 for(
size_t i=0;i<
np;i++) {
676 std::vector<size_t> index;
677 size_t max_smallest=(nd_in+2)*2;
678 if (max_smallest>np) max_smallest=
np;
679 if (max_smallest<nd_in+1) {
684 std::cout <<
"max_smallest: " << max_smallest << std::endl;
687 o2scl::vector_smallest_index<std::vector<double>,double,
688 std::vector<size_t> >(dists,max_smallest,index);
691 for(
size_t i=0;i<index.size();i++) {
692 std::cout <<
"index[" << i <<
"] = " << index[i] <<
" " 693 << dists[index[i]] << std::endl;
697 std::vector<size_t> index2;
698 for(
size_t i=0;i<max_smallest;i++) {
699 if (dists[index[i]]>0.0) {
700 index2.push_back(index[i]);
701 if (index2.size()==nd_in+1) i=max_smallest;
704 if (index2.size()<nd_in+1) {
705 O2SCL_ERR(
"Couldn't find enough nearby points (2).",
710 for(
size_t i=0;i<index2.size();i++) {
711 std::cout <<
"index2[" << i <<
"] = " << index2[i] <<
" " 712 << dists[index2[i]] << std::endl;
717 std::vector<ubvector> units(nd_in+1);
719 std::vector<double> diff_norms(nd_in+1);
721 std::vector<ubvector> ders(nd_in+1);
723 ubmatrix m(nd_in,nd_in);
727 std::vector<ubvector> ders2(nd_in);
729 for(
size_t i=0;i<nd_in+1;i++) {
732 units[i].resize(nd_in);
733 for(
size_t j=0;j<
nd_in;j++) {
734 units[i][j]=
ptrs[j][index2[i]]-x[j];
738 diff_norms[i]=o2scl::vector_norm<ubvector,double>(units[i]);
739 for(
size_t j=0;j<
nd_in;j++) {
740 units[i][j]/=diff_norms[i];
747 std::cout <<
"Point: ";
748 for(
size_t i=0;i<
nd_in;i++) {
749 std::cout << x[i] <<
" ";
751 std::cout << f << std::endl;
752 for(
size_t j=0;j<nd_in+1;j++) {
753 std::cout <<
"Closest: " << j <<
" " << index2[j] <<
" ";
754 for(
size_t i=0;i<
nd_in;i++) {
755 std::cout <<
ptrs[i][index2[j]] <<
" ";
757 std::cout <<
ptrs[func_index+
nd_in][index2[j]] <<
" " 758 << diff_norms[j] << std::endl;
760 for(
size_t j=0;j<nd_in+1;j++) {
761 std::cout <<
"diff_norm: " << j <<
" " << diff_norms[j]
768 for(
size_t i=0;i<nd_in+1;i++) {
770 ders[i].resize(nd_in);
774 for(
size_t j=0;j<nd_in+1;j++) {
776 for(
size_t k=0;k<
nd_in;k++) {
779 v[jj]=(
ptrs[func_index+
nd_in][index2[j]]-f)/diff_norms[j];
786 std::cout <<
"m:" << std::endl;
788 std::cout <<
"v:" << std::endl;
791 lshh.
solve(nd_in,m,v,ders[i]);
793 std::cout <<
"Derivs: " << i <<
" ";
794 std::cout.setf(std::ios::showpos);
795 for(
size_t j=0;j<
nd_in;j++) {
796 std::cout << ders[i][j] <<
" ";
798 std::cout.unsetf(std::ios::showpos);
799 std::cout << std::endl;
805 for(
size_t i=0;i<
nd_in;i++) {
808 ders2[i].resize(nd_in+1);
809 for(
size_t j=0;j<nd_in+1;j++) {
810 ders2[i][j]=ders[j][i];
822 #ifndef DOXYGEN_INTERNAL 847 template<
class vec2_t>
double dist(
size_t index,
848 const vec2_t &x)
const {
850 size_t nscales=scales.size();
851 for(
size_t i=0;i<
nd_in;i++) {
852 ret+=pow((x[i]-ptrs[i][index])/scales[i%nscales],2.0);
860 double dist(
size_t j,
size_t k)
const {
862 size_t nscales=scales.size();
863 for(
size_t i=0;i<
nd_in;i++) {
864 ret+=pow((ptrs[i][j]-ptrs[i][k])/scales[i%nscales],2.0);
874 #ifndef DOXYGEN_NO_O2NS void derivs_err(size_t func_index, size_t point_index, vec3_t &derivs, vec3_t &errs) const
For one of the functions, compute the partial derivatives (and uncertainties) with respect to all of ...
double vector_mean(size_t n, const vec_t &data)
Compute the mean of the first n elements of a vector.
Multi-dimensional interpolation by inverse distance weighting.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
invalid argument supplied by user
void get_data(size_t &n_in, size_t &n_out, size_t &n_points, vec_vec_t &vecs)
Get the data used for interpolation.
virtual void solve(size_t n, mat_t &A, vec_t &b, vec_t &x)
Solve square linear system of size n.
Generic Householder linear solver.
int matrix_out(std::ostream &os, size_t nrows, size_t ncols, mat_t &A)
A operator for simple matrix output using operator()
size_t nd_out
The number of dimensions of the outputs.
double operator()(const vec2_t &x) const
Perform the interpolation over the first function.
void eval_err_index(const vec2_t &x, vec3_t &val, vec4_t &err, std::vector< size_t > &index) const
Perform the interpolation over all the functions giving uncertainties and the sorted index vector...
double vector_stddev(size_t n, const vec_t &data)
Standard deviation with specified mean.
void auto_scale()
Automatically determine the length scales from the data.
ubvector scales
Distance scales for each coordinate.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
std::vector< vec_t > ptrs
A vector of pointers holding the data.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
double dist(size_t index, const vec2_t &x) const
Compute the distance between x and the point at index index.
double min_dist
The minimum distance to consider points as non-degenerate (default )
void set_points(size_t n)
Set the number of closest points to use for each interpolation (default 3)
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
void eval_err(const vec2_t &x, vec3_t &val, vec4_t &err) const
Perform the interpolation over all the functions with uncertainties.
size_t np
The number of points.
bool data_set
True if the data has been specified.
void set_data(size_t n_in, size_t n_points, vec_vec_t &vecs, bool auto_scale=true)
Initialize the data for the interpolation for only one output function.
void vector_out(std::ostream &os, size_t n, const vec_t &v, bool endline=false)
Output the first n elements of a vector to a stream, os.
size_t nd_in
The number of dimensions of the inputs.
int verbose
Verbosity parameter (default 0)
double dist(size_t j, size_t k) const
Compute the distance between two points in the data set.
void set_scales(size_t n, vec2_t &v)
Set the scales for the distance metric.
size_t points
Number of points to include in each interpolation (default 3)
void set_data(size_t n_in, size_t n_out, size_t n_points, vec_vec_t &vecs, bool auto_scale_flag=true)
Initialize the data for the interpolation.
size_t n_extra
The number of extra nearest neighbors to include to avoid degeneracies (default 0) ...
double eval(const vec2_t &x) const
Perform the interpolation over the first function.
void eval_err(const vec2_t &x, double &val, double &err) const
Perform the interpolation over the first function with uncertainty.
void eval(vec2_t &x, vec3_t &y) const
Perform the interpolation over all the functions, storing the result in y.