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 119 O2SCL_ERR(
"Order cannot be zero in interpm_idw.",
131 template<
class vec2_t>
void set_scales(
size_t n, vec2_t &v) {
133 O2SCL_ERR(
"Scale vector size cannot be zero in interpm_idw.",
136 for(
size_t i=0;i<n;i++) {
138 O2SCL_ERR(
"Scale must be positive and non-zero in interpm_idw.",
154 template<
class vec_vec_t>
155 void set_data(
size_t n_in,
size_t n_out,
size_t n_points,
156 vec_vec_t &vecs,
bool auto_scale_flag=
true) {
159 O2SCL_ERR2(
"Must provide at least three points in ",
163 O2SCL_ERR2(
"Must provide at least one input column in ",
167 O2SCL_ERR2(
"Must provide at least one output column in ",
173 ptrs.resize(n_in+n_out);
174 for(
size_t i=0;i<n_in+n_out;i++) {
175 std::swap(
ptrs[i],vecs[i]);
179 if (auto_scale_flag) {
188 template<
class vec_vec_t>
189 void get_data(
size_t &n_in,
size_t &n_out,
size_t &n_points,
194 vecs.resize(n_in+n_out);
195 for(
size_t i=0;i<n_in+n_out;i++) {
196 std::swap(
ptrs[i],vecs[i]);
211 for(
size_t i=0;i<
nd_in;i++) {
212 scales[i]=fabs(o2scl::vector_max_value<vec_t,double>(
np,
ptrs[i])-
213 o2scl::vector_min_value<vec_t,double>(
np,
ptrs[i]));
226 template<
class vec_vec_t>
235 template<
class vec2_t>
double operator()(
const vec2_t &x)
const {
241 template<
class vec2_t>
double eval(
const vec2_t &x)
const {
244 O2SCL_ERR(
"Data not set in interpm_idw::eval().",
249 std::vector<double> dists(
np);
250 for(
size_t i=0;i<
np;i++) {
255 std::vector<size_t> index;
256 o2scl::vector_smallest_index<std::vector<double>,double,
257 std::vector<size_t> >(dists,
order,index);
260 if (dists[index[0]]<=0.0) {
266 for(
size_t i=0;i<
order;i++) {
267 norm+=1.0/dists[index[i]];
272 for(
size_t i=0;i<
order;i++) {
273 ret+=
ptrs[
nd_in][index[i]]/dists[index[i]];
284 template<
class vec2_t>
void eval_err(
const vec2_t &x,
double &val,
288 O2SCL_ERR(
"Data not set in interpm_idw::eval_err().",
293 std::vector<double> dists(
np);
294 for(
size_t i=0;i<
np;i++) {
299 std::vector<size_t> index;
300 o2scl::vector_smallest_index<std::vector<double>,double,
301 std::vector<size_t> >(dists,
order+1,index);
303 if (dists[index[0]]<=0.0) {
312 std::vector<double> vals(order+1);
314 for(
size_t j=0;j<order+1;j++) {
318 for(
size_t i=0;i<order+1;i++) {
319 if (i!=j) norm+=1.0/dists[index[i]];
324 for(
size_t i=0;i<order+1;i++) {
326 vals[j]+=
ptrs[
nd_in][index[i]]/dists[index[i]];
344 template<
class vec2_t,
class vec3_t>
345 void eval(vec2_t &x, vec3_t &y)
const {
348 O2SCL_ERR(
"Data not set in interpm_idw::eval().",
353 std::cout <<
"interpm_idw: input: ";
354 for(
size_t k=0;k<
nd_in;k++) {
355 std::cout << x[k] <<
" ";
357 std::cout << std::endl;
361 std::vector<double> dists(
np);
362 for(
size_t i=0;i<
np;i++) {
367 std::vector<size_t> index;
368 o2scl::vector_smallest_index<std::vector<double>,double,
369 std::vector<size_t> >(dists,
order,index);
371 for(
size_t i=0;i<
order;i++) {
372 std::cout <<
"interpm_idw: closest point: ";
373 for(
size_t k=0;k<
nd_in;k++) {
374 std::cout <<
ptrs[k][index[i]] <<
" ";
376 std::cout << std::endl;
382 if (dists[index[0]]<=0.0) {
383 for(
size_t i=0;i<
nd_out;i++) {
387 std::cout <<
"interpm_idw: distance zero. " 388 <<
"Returning values at index: " << index[0] << std::endl;
397 for(
size_t i=0;i<
order;i++) {
398 norm+=1.0/dists[index[i]];
401 std::cout <<
"interpm_idw: norm is " << norm << std::endl;
405 for(
size_t j=0;j<
nd_out;j++) {
407 for(
size_t i=0;i<
order;i++) {
408 if (j==0 && verbose>0) {
409 std::cout <<
"interpm_idw: Point: ";
410 for(
size_t k=0;k<
nd_in;k++) {
411 std::cout <<
ptrs[k][index[i]] <<
" ";
413 std::cout << std::endl;
415 y[j]+=
ptrs[
nd_in+j][index[i]]/dists[index[i]];
417 std::cout <<
"interpm_idw: j,order,value,1/dist: " 418 << j <<
" " << i <<
" " 420 << 1.0/dists[index[i]] << std::endl;
425 std::cout <<
"interpm_idw: y[" << j <<
"]: " << y[j]
438 template<
class vec2_t,
class vec3_t,
class vec4_t>
440 std::vector<size_t> &index)
const {
443 O2SCL_ERR(
"Data not set in interpm_idw::eval_err().",
448 std::vector<double> dists(
np);
449 for(
size_t i=0;i<
np;i++) {
454 o2scl::vector_smallest_index<std::vector<double>,double,
455 std::vector<size_t> >(dists,
order+1,index);
457 if (dists[index[0]]<=0.0) {
461 for(
size_t k=0;k<
nd_out;k++) {
469 for(
size_t k=0;k<
nd_out;k++) {
471 std::vector<double> vals(order+1);
473 for(
size_t j=0;j<order+1;j++) {
477 for(
size_t i=0;i<order+1;i++) {
478 if (i!=j) norm+=1.0/dists[index[i]];
483 for(
size_t i=0;i<order+1;i++) {
485 vals[j]+=
ptrs[
nd_in+k][index[i]]/dists[index[i]];
509 template<
class vec2_t,
class vec3_t,
class vec4_t>
510 void eval_err(
const vec2_t &x, vec3_t &val, vec4_t &err)
const {
511 std::vector<size_t> index;
527 template<
class vec3_t>
529 vec3_t &derivs, vec3_t &errs)
const {
532 O2SCL_ERR(
"Data not set in interpm_idw::derivs_err().",
538 for(
size_t i=0;i<
nd_in;i++) {
539 x[i]=
ptrs[i][point_index];
542 double f=
ptrs[nd_in+func_index][point_index];
548 std::vector<double> dists(
np);
549 for(
size_t i=0;i<
np;i++) {
555 std::vector<size_t> index;
556 size_t max_smallest=(nd_in+2)*2;
557 if (max_smallest>np) max_smallest=
np;
558 if (max_smallest<nd_in+1) {
563 std::cout <<
"max_smallest: " << max_smallest << std::endl;
566 o2scl::vector_smallest_index<std::vector<double>,double,
567 std::vector<size_t> >(dists,max_smallest,index);
570 for(
size_t i=0;i<index.size();i++) {
571 std::cout <<
"index[" << i <<
"] = " << index[i] <<
" " 572 << dists[index[i]] << std::endl;
576 std::vector<size_t> index2;
577 for(
size_t i=0;i<max_smallest;i++) {
578 if (dists[index[i]]>0.0) {
579 index2.push_back(index[i]);
580 if (index2.size()==nd_in+1) i=max_smallest;
583 if (index2.size()<nd_in+1) {
584 O2SCL_ERR(
"Couldn't find enough nearby points (2).",
589 for(
size_t i=0;i<index2.size();i++) {
590 std::cout <<
"index2[" << i <<
"] = " << index2[i] <<
" " 591 << dists[index2[i]] << std::endl;
596 std::vector<ubvector> units(nd_in+1);
598 std::vector<double> diff_norms(nd_in+1);
600 std::vector<ubvector> ders(nd_in+1);
602 ubmatrix m(nd_in,nd_in);
606 std::vector<ubvector> ders2(nd_in);
608 for(
size_t i=0;i<nd_in+1;i++) {
611 units[i].resize(nd_in);
612 for(
size_t j=0;j<
nd_in;j++) {
613 units[i][j]=
ptrs[j][index2[i]]-x[j];
617 diff_norms[i]=o2scl::vector_norm<ubvector,double>(units[i]);
618 for(
size_t j=0;j<
nd_in;j++) {
619 units[i][j]/=diff_norms[i];
626 std::cout <<
"Point: ";
627 for(
size_t i=0;i<
nd_in;i++) {
628 std::cout << x[i] <<
" ";
630 std::cout << f << std::endl;
631 for(
size_t j=0;j<nd_in+1;j++) {
632 std::cout <<
"Closest: " << j <<
" " << index2[j] <<
" ";
633 for(
size_t i=0;i<
nd_in;i++) {
634 std::cout <<
ptrs[i][index2[j]] <<
" ";
636 std::cout <<
ptrs[func_index+
nd_in][index2[j]] <<
" " 637 << diff_norms[j] << std::endl;
639 for(
size_t j=0;j<nd_in+1;j++) {
640 std::cout <<
"diff_norm: " << j <<
" " << diff_norms[j]
647 for(
size_t i=0;i<nd_in+1;i++) {
649 ders[i].resize(nd_in);
653 for(
size_t j=0;j<nd_in+1;j++) {
655 for(
size_t k=0;k<
nd_in;k++) {
658 v[jj]=(
ptrs[func_index+
nd_in][index2[j]]-f)/diff_norms[j];
665 std::cout <<
"m:" << std::endl;
667 std::cout <<
"v:" << std::endl;
670 lshh.
solve(nd_in,m,v,ders[i]);
672 std::cout <<
"Derivs: " << i <<
" ";
673 std::cout.setf(std::ios::showpos);
674 for(
size_t j=0;j<
nd_in;j++) {
675 std::cout << ders[i][j] <<
" ";
677 std::cout.unsetf(std::ios::showpos);
678 std::cout << std::endl;
684 for(
size_t i=0;i<
nd_in;i++) {
687 ders2[i].resize(nd_in+1);
688 for(
size_t j=0;j<nd_in+1;j++) {
689 ders2[i][j]=ders[j][i];
700 #ifndef DOXYGEN_INTERNAL 721 template<
class vec2_t>
double dist(
size_t index,
722 const vec2_t &x)
const {
724 size_t nscales=scales.size();
725 for(
size_t i=0;i<
nd_in;i++) {
726 ret+=pow((x[i]-ptrs[i][index])/scales[i%nscales],2.0);
735 #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.
void set_order(size_t n)
Set the number of closest points to use for each interpolation (default 3)
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 with uncertainties.
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.
size_t order
Number of points to include in 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)
void set_scales(size_t n, vec2_t &v)
Set the scales for the distance metric.
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.
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.