interpm_idw.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_INTERPM_IDW_H
24 #define O2SCL_INTERPM_IDW_H
25 
26 /** \file interpm_idw.h
27  \brief File defining \ref o2scl::interpm_idw
28 */
29 
30 #include <iostream>
31 #include <string>
32 #include <cmath>
33 
34 #include <boost/numeric/ublas/matrix.hpp>
35 
36 #include <gsl/gsl_combination.h>
37 
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 
44 #ifndef DOXYGEN_NO_O2NS
45 namespace o2scl {
46 #endif
47 
48  /** \brief Multi-dimensional interpolation by inverse distance
49  weighting
50 
51  This class is experimental, particularly the evaluation of
52  derivatives.
53 
54  This class performs interpolation on a multi-dimensional data
55  set specified as a series of scattered points using the inverse
56  distance-weighted average of nearby points. This should be
57  superior to Gaussian process interpolation as the typical
58  distance betweeen samples becomes smaller than the correlation
59  length. (?)
60 
61  The function \ref set_data() takes as input: the number of input
62  dimensions, the number of output functions, the number of points
63  which specify the data, and a "vector of vectors" which contains
64  the data for all the points. The vector of vectors must be of a
65  type which allows std::swap on individual elements (which are of
66  type <tt>vec_t</tt>).
67 
68  The number of nearby points which are averaged defaults to 3 and
69  can be changed by \ref set_points(). To obtain interpolation
70  uncertainties, this class finds one additional nearby point and
71  returns the standard deviation of the interpolated value for all
72  subsets which are missing one nearby point. One-point
73  interpolation corresponds to nearest-neighbor interpolation.
74 
75  This class requires a distance metric to weight the
76  interpolation, and a Euclidean distance is used. By default, the
77  length scales in each direction are automatically determined by
78  extent of the data (absolute value of max minus min in each
79  direction), but the user can specify the length scales manually
80  in \ref set_scales() .
81 
82  First derivatives can be obtained using \ref derivs_err() , but
83  these derivatives are not obtained from the same approximation
84  used in the interpolation. That is, the derivatives returned are
85  not equal to exact derivatives from the interpolated function
86  (as is the case in, e.g., cubic spline interpolation in one
87  dimension). This will typically only be particularly noticable
88  near discontinuities.
89 
90  If the user specifies an array of pointers, the data can be
91  changed between calls to the interpolation, but data points
92  cannot be added (as set data separately stores the total number
93  of data points) without a new call to \ref set_data(). Also, the
94  automatically-determined length scales may need to be recomputed
95  by calling \ref auto_scale().
96 
97  Increasing the value of \c n_extra away from zero allows the
98  interpolation to ignore points in the data set which are
99  degenerate because they are too close to each other. Points with
100  a distance (normalized by the scales) less than \ref min_dist
101  are automatically considered degenerate and only the single
102  point closest to the requested coordinate is considered.
103  Increasing the value of \c n_extra increases the computational
104  time required to compute the nearest points which are
105  nondegenerate.
106 
107  \todo Make verbose output consistent between the various
108  eval() functions.
109 
110  \future Share code between the various functions
111  */
112  template<class vec_t=boost::numeric::ublas::vector<double> >
113  class interpm_idw {
114 
115  protected:
116 
117  public:
118 
122 
123  interpm_idw() {
124  data_set=false;
125  scales.resize(1);
126  scales[0]=1.0;
127  points=3;
128  verbose=0;
129  n_extra=0;
130  min_dist=1.0e-6;
131  }
132 
133  /** \brief Verbosity parameter (default 0)
134  */
135  int verbose;
136 
137  /** \brief The number of extra nearest neighbors
138  to include to avoid degeneracies (default 0)
139  */
140  size_t n_extra;
141 
142  /** \brief The minimum distance to consider points as
143  non-degenerate (default \f$ 10^{-6} \f$ )
144  */
145  double min_dist;
146 
147  /// \name Get and set functions
148  //@{
149  /** \brief Set the number of closest points to use
150  for each interpolation (default 3)
151  */
152  void set_points(size_t n) {
153  if (n==0) {
154  O2SCL_ERR("Points cannot be zero in interpm_idw.",
156  }
157  points=n;
158  return;
159  }
160 
161  /** \brief Set the scales for the distance metric
162 
163  All the scales must be positive and non-zero. The size of the
164  vector \c (specified in \c n) must be larger than zero.
165  */
166  template<class vec2_t> void set_scales(size_t n, vec2_t &v) {
167  if (n==0) {
168  O2SCL_ERR("Scale vector size cannot be zero in interpm_idw.",
170  }
171  for(size_t i=0;i<n;i++) {
172  if (v[i]<=0.0) {
173  O2SCL_ERR("Scale must be positive and non-zero in interpm_idw.",
175  }
176  }
177  scales.resize(n);
179  return;
180  }
181 
182  /** \brief Initialize the data for the interpolation
183 
184  The object \c vecs should be a vector (of size <tt>n_in+n_out</tt>)
185  of vectors (all of size <tt>n_points</tt>). It may have be
186  any time which allows the use of <tt>std::swap</tt> for
187  each vector in the list.
188  */
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) {
192 
193  if (n_points<3) {
194  O2SCL_ERR2("Must provide at least three points in ",
195  "interpm_idw::set_data()",exc_efailed);
196  }
197  if (n_in<1) {
198  O2SCL_ERR2("Must provide at least one input column in ",
199  "interpm_idw::set_data()",exc_efailed);
200  }
201  if (n_out<1) {
202  O2SCL_ERR2("Must provide at least one output column in ",
203  "interpm_idw::set_data()",exc_efailed);
204  }
205  np=n_points;
206  nd_in=n_in;
207  nd_out=n_out;
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]);
211  }
212  data_set=true;
213 
214  if (auto_scale_flag) {
215  auto_scale();
216  }
217 
218  return;
219  }
220 
221  /** \brief Get the data used for interpolation
222  */
223  template<class vec_vec_t>
224  void get_data(size_t &n_in, size_t &n_out, size_t &n_points,
225  vec_vec_t &vecs) {
226  n_points=np;
227  n_in=nd_in;
228  n_out=nd_out;
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]);
232  }
233  data_set=false;
234  n_points=0;
235  n_in=0;
236  n_out=0;
237  ptrs.clear();
238  return;
239  }
240 
241  /** \brief Automatically determine the length scales from the
242  data
243  */
244  void auto_scale() {
245  scales.resize(nd_in);
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]));
249  }
250  return;
251  }
252 
253  /** \brief Initialize the data for the interpolation
254  for only one output function
255 
256  The object \c vecs should be a vector (of size <tt>n_in+1</tt>)
257  of vectors (all of size <tt>n_points</tt>). It may be
258  any type which allows the use of <tt>std::swap</tt> for
259  each vector in the list.
260  */
261  template<class vec_vec_t>
262  void set_data(size_t n_in, size_t n_points,
263  vec_vec_t &vecs, bool auto_scale=true) {
264  set_data(n_in,1,n_points,vecs,auto_scale);
265  return;
266  }
267  //@}
268 
269  /// \name Evaluate interpolation
270  //@{
271  /** \brief Perform the interpolation over the first function
272  */
273  template<class vec2_t> double operator()(const vec2_t &x) const {
274  return eval(x);
275  }
276 
277  /** \brief Perform the interpolation over the first function
278  */
279  template<class vec2_t> double eval(const vec2_t &x) const {
280 
281  if (data_set==false) {
282  O2SCL_ERR("Data not set in interpm_idw::eval().",
283  exc_einval);
284  }
285 
286  // Compute distances
287  std::vector<double> dists(np);
288  for(size_t i=0;i<np;i++) {
289  dists[i]=dist(i,x);
290  }
291 
292  // Find closest points
293  std::vector<size_t> index;
294  o2scl::vector_smallest_index<std::vector<double>,double,
295  std::vector<size_t> >(dists,points+n_extra,index);
296 
297  if (n_extra>0) {
298  // Remove degenerate points to ensure accurate interpolation
299  bool found=true;
300  while (found==true) {
301  found=false;
302  // Find degenerate points and remove them
303  for(size_t j=0;j<points+n_extra;j++) {
304  for(size_t k=j;k<points+n_extra;k++) {
305  double dist_jk=dist(j,k);
306  if (index.size()>points && dist_jk<min_dist) {
307  found=true;
308  index.erase(index.begin()+j);
309  }
310  }
311  }
312  }
313  }
314 
315  // Check if the closest distance is zero
316  if (dists[index[0]]<=0.0) {
317  return ptrs[nd_in][index[0]];
318  }
319 
320  // Compute normalization
321  double norm=0.0;
322  for(size_t i=0;i<points;i++) {
323  norm+=1.0/dists[index[i]];
324  }
325 
326  // Compute the inverse-distance weighted average
327  double ret=0.0;
328  for(size_t i=0;i<points;i++) {
329  ret+=ptrs[nd_in][index[i]]/dists[index[i]];
330  }
331  ret/=norm;
332 
333  // Return the average
334  return ret;
335  }
336 
337  /** \brief Perform the interpolation over the first function
338  with uncertainty
339  */
340  template<class vec2_t> void eval_err(const vec2_t &x, double &val,
341  double &err) const {
342 
343  if (data_set==false) {
344  O2SCL_ERR("Data not set in interpm_idw::eval_err().",
345  exc_einval);
346  }
347 
348  // Compute distances
349  std::vector<double> dists(np);
350  for(size_t i=0;i<np;i++) {
351  dists[i]=dist(i,x);
352  }
353 
354  // Find closest points
355  std::vector<size_t> index;
356  o2scl::vector_smallest_index<std::vector<double>,double,
357  std::vector<size_t> >(dists,points+1+n_extra,index);
358 
359  if (n_extra>0) {
360  // Remove degenerate points to ensure accurate interpolation
361  bool found=true;
362  while (found==true) {
363  found=false;
364  // Find degenerate points and remove them
365  for(size_t j=0;j<points+1+n_extra;j++) {
366  for(size_t k=j;k<points+1+n_extra;k++) {
367  double dist_jk=dist(j,k);
368  if (index.size()>points+1 && dist_jk<min_dist) {
369  found=true;
370  index.erase(index.begin()+j);
371  }
372  }
373  }
374  }
375  }
376 
377  if (dists[index[0]]<=0.0) {
378 
379  // If the closest distance is zero, just set the value
380  val=ptrs[nd_in][index[0]];
381  err=0.0;
382  return;
383 
384  } else {
385 
386  std::vector<double> vals(points+1);
387 
388  for(size_t j=0;j<points+1;j++) {
389 
390  // Compute normalization
391  double norm=0.0;
392  for(size_t i=0;i<points+1;i++) {
393  if (i!=j) norm+=1.0/dists[index[i]];
394  }
395 
396  // Compute the inverse-distance weighted average
397  vals[j]=0.0;
398  for(size_t i=0;i<points+1;i++) {
399  if (i!=j) {
400  vals[j]+=ptrs[nd_in][index[i]]/dists[index[i]];
401  }
402  }
403  vals[j]/=norm;
404 
405  }
406 
407  val=vals[points];
408  err=o2scl::vector_stddev(vals);
409 
410  }
411 
412  return;
413  }
414 
415  /** \brief Perform the interpolation over all the functions,
416  storing the result in \c y
417  */
418  template<class vec2_t, class vec3_t>
419  void eval(vec2_t &x, vec3_t &y) const {
420 
421  if (data_set==false) {
422  O2SCL_ERR("Data not set in interpm_idw::eval().",
423  exc_einval);
424  }
425 
426  if (verbose>0) {
427  std::cout << "interpm_idw: input: ";
428  for(size_t k=0;k<nd_in;k++) {
429  std::cout << x[k] << " ";
430  }
431  std::cout << std::endl;
432  }
433 
434  // Compute distances
435  std::vector<double> dists(np);
436  for(size_t i=0;i<np;i++) {
437  dists[i]=dist(i,x);
438  }
439 
440  // Find closest points
441  std::vector<size_t> index;
442  o2scl::vector_smallest_index<std::vector<double>,double,
443  std::vector<size_t> >(dists,points,index);
444  if (verbose>0) {
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]] << " ";
449  }
450  std::cout << std::endl;
451  }
452  }
453 
454  if (n_extra>0) {
455  // Remove degenerate points to ensure accurate interpolation
456  bool found=true;
457  while (found==true) {
458  found=false;
459  // Find degenerate points and remove them
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) {
464  found=true;
465  index.erase(index.begin()+j);
466  }
467  }
468  }
469  }
470  }
471 
472  // Check if the closest distance is zero, if so, just
473  // return the value
474  if (dists[index[0]]<=0.0) {
475  for(size_t i=0;i<nd_out;i++) {
476  y[i]=ptrs[nd_in+i][index[0]];
477  }
478  if (verbose>0) {
479  std::cout << "interpm_idw: distance zero. "
480  << "Returning values at index: " << index[0] << std::endl;
481  std::cout << "\t";
482  o2scl::vector_out(std::cout,nd_out,y,true);
483  }
484  return;
485  }
486 
487  // Compute normalization
488  double norm=0.0;
489  for(size_t i=0;i<points;i++) {
490  norm+=1.0/dists[index[i]];
491  }
492  if (verbose>0) {
493  std::cout << "interpm_idw: norm is " << norm << std::endl;
494  }
495 
496  // Compute the inverse-distance weighted averages
497  for(size_t j=0;j<nd_out;j++) {
498  y[j]=0.0;
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]] << " ";
504  }
505  std::cout << std::endl;
506  }
507  y[j]+=ptrs[nd_in+j][index[i]]/dists[index[i]];
508  if (verbose>0) {
509  std::cout << "interpm_idw: j,points,value,1/dist: "
510  << j << " " << i << " "
511  << ptrs[nd_in+j][index[i]] << " "
512  << 1.0/dists[index[i]] << std::endl;
513  }
514  }
515  y[j]/=norm;
516  if (verbose>0) {
517  std::cout << "interpm_idw: y[" << j << "]: " << y[j]
518  << std::endl;
519  }
520  }
521 
522 
523 
524  return;
525  }
526 
527  /** \brief Perform the interpolation over all the functions
528  giving uncertainties and the sorted index vector
529 
530  The vector \c index is automatically resized to
531  a size equal to n_points+1+n_extra
532  be larger than
533  */
534  template<class vec2_t, class vec3_t, class vec4_t>
535  void eval_err_index(const vec2_t &x, vec3_t &val, vec4_t &err,
536  std::vector<size_t> &index) const {
537 
538  if (data_set==false) {
539  O2SCL_ERR("Data not set in interpm_idw::eval_err().",
540  exc_einval);
541  }
542 
543  // Compute distances
544  std::vector<double> dists(np);
545  for(size_t i=0;i<np;i++) {
546  dists[i]=dist(i,x);
547  }
548 
549  // Find closest points, note that index is automatically resized
550  // by the vector_smallest_index function
551  o2scl::vector_smallest_index<std::vector<double>,double,
552  std::vector<size_t> >(dists,points+1+n_extra,index);
553 
554  if (n_extra>0) {
555  // Remove degenerate points to ensure accurate interpolation
556  bool found=true;
557  while (found==true) {
558  found=false;
559  // Find degenerate points and remove them
560  for(size_t j=0;j<points+1+n_extra;j++) {
561  for(size_t k=j;k<points+1+n_extra;k++) {
562  double dist_jk=dist(j,k);
563  if (index.size()>points+1 && dist_jk<min_dist) {
564  found=true;
565  index.erase(index.begin()+j);
566  }
567  }
568  }
569  }
570  }
571 
572  if (dists[index[0]]<=0.0) {
573 
574  // If the closest distance is zero, just set the values and
575  // errors
576  for(size_t k=0;k<nd_out;k++) {
577  val[k]=ptrs[nd_in+k][index[0]];
578  err[k]=0.0;
579  }
580  return;
581 
582  } else {
583 
584  for(size_t k=0;k<nd_out;k++) {
585 
586  std::vector<double> vals(points+1);
587 
588  for(size_t j=0;j<points+1;j++) {
589 
590  // Compute normalization
591  double norm=0.0;
592  for(size_t i=0;i<points+1;i++) {
593  if (i!=j) norm+=1.0/dists[index[i]];
594  }
595 
596  // Compute the inverse-distance weighted average
597  vals[j]=0.0;
598  for(size_t i=0;i<points+1;i++) {
599  if (i!=j) {
600  vals[j]+=ptrs[nd_in+k][index[i]]/dists[index[i]];
601  }
602  }
603  vals[j]/=norm;
604 
605  }
606 
607  // Instead of using the average, we report the value as the
608  // last element in the array, which is the interpolated
609  // value from the closest points
610  val[k]=vals[points];
611 
612  err[k]=o2scl::vector_stddev(vals);
613 
614  }
615 
616  }
617 
618  return;
619  }
620 
621  /** \brief Perform the interpolation over all the functions
622  with uncertainties
623  */
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;
627  return eval_err_index(x,val,err,index);
628  }
629  //@}
630 
631  /// \name Evaluate derivatives
632  //@{
633  /** \brief For one of the functions, compute the partial
634  derivatives (and uncertainties) with respect to all of the
635  inputs at one data point
636 
637  \note This function ignores the points chosen by \ref
638  set_points() and always chooses to average derivative
639  calculations determined from \c n_in+1 combinations of \c n_in
640  points .
641 
642  \todo Use the mechanism provided by <tt>n_extra</tt> above
643  to remove degenerate points.
644 
645  \future This function requires an extra copy from
646  "ders" to "ders2" which could be removed.
647  */
648  template<class vec3_t>
649  void derivs_err(size_t func_index, size_t point_index,
650  vec3_t &derivs, vec3_t &errs) const {
651 
652  if (data_set==false) {
653  O2SCL_ERR("Data not set in interpm_idw::derivs_err().",
654  exc_einval);
655  }
656 
657  // Set x equal to the specified point
658  ubvector x(nd_in);
659  for(size_t i=0;i<nd_in;i++) {
660  x[i]=ptrs[i][point_index];
661  }
662  // Set f equal to the value of the function at the specified point
663  double f=ptrs[nd_in+func_index][point_index];
664 
665  // The linear solver
667 
668  // Compute distances
669  std::vector<double> dists(np);
670  for(size_t i=0;i<np;i++) {
671  dists[i]=dist(i,x);
672  }
673 
674  // Find closest (but not identical) points
675 
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) {
680  O2SCL_ERR("Couldn't find enough nearby points.",o2scl::exc_einval);
681  }
682 
683  if (verbose>0) {
684  std::cout << "max_smallest: " << max_smallest << std::endl;
685  }
686 
687  o2scl::vector_smallest_index<std::vector<double>,double,
688  std::vector<size_t> >(dists,max_smallest,index);
689 
690  if (verbose>0) {
691  for(size_t i=0;i<index.size();i++) {
692  std::cout << "index[" << i << "] = " << index[i] << " "
693  << dists[index[i]] << std::endl;
694  }
695  }
696 
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;
702  }
703  }
704  if (index2.size()<nd_in+1) {
705  O2SCL_ERR("Couldn't find enough nearby points (2).",
707  }
708 
709  if (verbose>0) {
710  for(size_t i=0;i<index2.size();i++) {
711  std::cout << "index2[" << i << "] = " << index2[i] << " "
712  << dists[index2[i]] << std::endl;
713  }
714  }
715 
716  // Unit vector storage
717  std::vector<ubvector> units(nd_in+1);
718  // Difference vector norms
719  std::vector<double> diff_norms(nd_in+1);
720  // Storage for the derivative estimates
721  std::vector<ubvector> ders(nd_in+1);
722  // Matrix of unit vectors
723  ubmatrix m(nd_in,nd_in);
724  // Vector of function value differences
725  ubvector v(nd_in);
726  // Rearranged derivative object
727  std::vector<ubvector> ders2(nd_in);
728 
729  for(size_t i=0;i<nd_in+1;i++) {
730 
731  // Assign unit vector elements
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];
735  }
736 
737  // Normalize the unit vectors
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];
741  }
742 
743  }
744 
745  // Verbose output of the closest points and their norms
746  if (verbose>0) {
747  std::cout << "Point: ";
748  for(size_t i=0;i<nd_in;i++) {
749  std::cout << x[i] << " ";
750  }
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]] << " ";
756  }
757  std::cout << ptrs[func_index+nd_in][index2[j]] << " "
758  << diff_norms[j] << std::endl;
759  }
760  for(size_t j=0;j<nd_in+1;j++) {
761  std::cout << "diff_norm: " << j << " " << diff_norms[j]
762  << std::endl;
763  }
764  // End of verbose output
765  }
766 
767  // Go through each set of points
768  for(size_t i=0;i<nd_in+1;i++) {
769 
770  ders[i].resize(nd_in);
771 
772  // Construct the matrix and vector for the solver
773  size_t jj=0;
774  for(size_t j=0;j<nd_in+1;j++) {
775  if (j!=i) {
776  for(size_t k=0;k<nd_in;k++) {
777  m(jj,k)=units[j][k];
778  }
779  v[jj]=(ptrs[func_index+nd_in][index2[j]]-f)/diff_norms[j];
780  jj++;
781  }
782  }
783 
784  // Solve to compute the derivatives
785  if (verbose>0) {
786  std::cout << "m:" << std::endl;
787  o2scl::matrix_out(std::cout,nd_in,nd_in,m);
788  std::cout << "v:" << std::endl;
789  o2scl::vector_out(std::cout,nd_in,v,true);
790  }
791  lshh.solve(nd_in,m,v,ders[i]);
792  if (verbose>0) {
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] << " ";
797  }
798  std::cout.unsetf(std::ios::showpos);
799  std::cout << std::endl;
800  }
801 
802  // Go to next derivative estimate
803  }
804 
805  for(size_t i=0;i<nd_in;i++) {
806 
807  // Rearrange derivatives
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];
811  }
812 
813  // Compute mean and standard deviation
814  derivs[i]=o2scl::vector_mean(ders2[i]);
815  errs[i]=o2scl::vector_stddev(ders2[i]);
816  }
817 
818  return;
819  }
820  //@}
821 
822 #ifndef DOXYGEN_INTERNAL
823 
824  protected:
825 
826  /// The number of points
827  size_t np;
828  /// The number of dimensions of the inputs
829  size_t nd_in;
830  /// The number of dimensions of the outputs
831  size_t nd_out;
832  /// A vector of pointers holding the data
833  std::vector<vec_t> ptrs;
834  /// True if the data has been specified
835  bool data_set;
836  /// Number of points to include in each interpolation (default 3)
837  size_t points;
838 
839  /// \name Distance determination [protected]
840  //@{
841  /// Distance scales for each coordinate
842  ubvector scales;
843 
844  /** \brief Compute the distance between \c x and the point at
845  index \c index
846  */
847  template<class vec2_t> double dist(size_t index,
848  const vec2_t &x) const {
849  double ret=0.0;
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);
853  }
854  return sqrt(ret);
855  }
856 
857  /** \brief Compute the distance between two points in the
858  data set
859  */
860  double dist(size_t j, size_t k) const {
861  double ret=0.0;
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);
865  }
866  return sqrt(ret);
867  }
868  //@}
869 
870 #endif
871 
872  };
873 
874 #ifndef DOXYGEN_NO_O2NS
875 }
876 #endif
877 
878 #endif
879 
880 
881 
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 ...
Definition: interpm_idw.h:649
double vector_mean(size_t n, const vec_t &data)
Compute the mean of the first n elements of a vector.
Definition: vec_stats.h:55
Multi-dimensional interpolation by inverse distance weighting.
Definition: interpm_idw.h:113
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
invalid argument supplied by user
Definition: err_hnd.h:59
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.
Definition: interpm_idw.h:224
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()
Definition: columnify.h:254
generic failure
Definition: err_hnd.h:61
size_t nd_out
The number of dimensions of the outputs.
Definition: interpm_idw.h:831
double operator()(const vec2_t &x) const
Perform the interpolation over the first function.
Definition: interpm_idw.h:273
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...
Definition: interpm_idw.h:535
double vector_stddev(size_t n, const vec_t &data)
Standard deviation with specified mean.
Definition: vec_stats.h:253
void auto_scale()
Automatically determine the length scales from the data.
Definition: interpm_idw.h:244
ubvector scales
Distance scales for each coordinate.
Definition: interpm_idw.h:842
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
std::vector< vec_t > ptrs
A vector of pointers holding the data.
Definition: interpm_idw.h:833
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
double dist(size_t index, const vec2_t &x) const
Compute the distance between x and the point at index index.
Definition: interpm_idw.h:847
double min_dist
The minimum distance to consider points as non-degenerate (default )
Definition: interpm_idw.h:145
void set_points(size_t n)
Set the number of closest points to use for each interpolation (default 3)
Definition: interpm_idw.h:152
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
void eval_err(const vec2_t &x, vec3_t &val, vec4_t &err) const
Perform the interpolation over all the functions with uncertainties.
Definition: interpm_idw.h:625
size_t np
The number of points.
Definition: interpm_idw.h:827
bool data_set
True if the data has been specified.
Definition: interpm_idw.h:835
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.
Definition: interpm_idw.h:262
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.
Definition: vector.h:2673
size_t nd_in
The number of dimensions of the inputs.
Definition: interpm_idw.h:829
int verbose
Verbosity parameter (default 0)
Definition: interpm_idw.h:135
double dist(size_t j, size_t k) const
Compute the distance between two points in the data set.
Definition: interpm_idw.h:860
void set_scales(size_t n, vec2_t &v)
Set the scales for the distance metric.
Definition: interpm_idw.h:166
size_t points
Number of points to include in each interpolation (default 3)
Definition: interpm_idw.h:837
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.
Definition: interpm_idw.h:190
size_t n_extra
The number of extra nearest neighbors to include to avoid degeneracies (default 0) ...
Definition: interpm_idw.h:140
double eval(const vec2_t &x) const
Perform the interpolation over the first function.
Definition: interpm_idw.h:279
void eval_err(const vec2_t &x, double &val, double &err) const
Perform the interpolation over the first function with uncertainty.
Definition: interpm_idw.h:340
void eval(vec2_t &x, vec3_t &y) const
Perform the interpolation over all the functions, storing the result in y.
Definition: interpm_idw.h:419

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).