interpm_idw.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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 #include <o2scl/table.h>
44 
45 #ifndef DOXYGEN_NO_O2NS
46 namespace o2scl {
47 #endif
48 
49  /** \brief Multi-dimensional interpolation by inverse distance
50  weighting
51 
52  This class is experimental, particularly the evaluation of
53  derivatives.
54 
55  This class performs interpolation on a multi-dimensional data
56  set specified as a series of scattered points using the inverse
57  distance-weighted average of nearby points.
58 
59  The function \ref set_data() takes as input: the number of input
60  dimensions, the number of output functions, the number of points
61  which specify the data, and a "vector of vectors" which contains
62  the data for all the points. The vector of vectors must be of a
63  type which allows std::swap on individual elements (which are of
64  type <tt>vec_t</tt>). The vector of vectors should be structured
65  so that the first index selects the input or output function
66  and the second index selects the data point.
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 mat_t=const_matrix_view_table<> >
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  dist_expo=2.0;
132  }
133 
134  /** \brief Exponent in computing distance (default 2.0)
135  */
136  double dist_expo;
137 
138  /** \brief Verbosity parameter (default 0)
139  */
140  int verbose;
141 
142  /** \brief The number of extra nearest neighbors
143  to include to avoid degeneracies (default 0)
144  */
145  size_t n_extra;
146 
147  /** \brief The minimum distance to consider points as
148  non-degenerate (default \f$ 10^{-6} \f$ )
149  */
150  double min_dist;
151 
152  /// \name Get and set functions
153  //@{
154  /** \brief Set the number of closest points to use
155  for each interpolation (default 3)
156  */
157  void set_points(size_t n) {
158  if (n==0) {
159  O2SCL_ERR("Points cannot be zero in interpm_idw.",
161  }
162  points=n;
163  return;
164  }
165 
166  /** \brief Set the scales for the distance metric
167 
168  All the scales must be positive and non-zero. The size of the
169  vector \c (specified in \c n) must be larger than zero.
170  */
171  template<class vec2_t> void set_scales(size_t n, vec2_t &v) {
172  if (n==0) {
173  O2SCL_ERR("Scale vector size cannot be zero in interpm_idw.",
175  }
176  for(size_t i=0;i<n;i++) {
177  if (v[i]<=0.0) {
178  O2SCL_ERR("Scale must be positive and non-zero in interpm_idw.",
180  }
181  }
182  scales.resize(n);
184  return;
185  }
186 
187  /** \brief Initialize the data for the interpolation
188 
189  The object \c vecs should be a matrix with a
190  first index of size <tt>n_in+n_out</tt> and a second
191  index of size <tt>n_points</tt>. It may have be
192  any type which allows the use of <tt>operator(,)</tt>
193  and <tt>std::swap</tt>.
194  */
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) {
197 
198  if (n_points<points) {
199  O2SCL_ERR2("Not enough points provided in ",
200  "interpm_idw::set_data()",exc_efailed);
201  }
202  if (n_in<1) {
203  O2SCL_ERR2("Must provide at least one input column in ",
204  "interpm_idw::set_data()",exc_efailed);
205  }
206  if (n_out<1) {
207  O2SCL_ERR2("Must provide at least one output column in ",
208  "interpm_idw::set_data()",exc_efailed);
209  }
210  np=n_points;
211  nd_in=n_in;
212  nd_out=n_out;
213  std::swap(data,dat);
214  data_set=true;
215 
216  if (auto_scale_flag) {
217  auto_scale();
218  }
219 
220  return;
221  }
222 
223  /** \brief Get the data used for interpolation
224  */
225  void get_data(size_t &n_in, size_t &n_out, size_t &n_points,
226  mat_t &dat) {
227  n_points=np;
228  n_in=nd_in;
229  n_out=nd_out;
230  std::swap(data,dat);
231  data_set=false;
232  n_points=0;
233  n_in=0;
234  n_out=0;
235  return;
236  }
237 
238  /** \brief Automatically determine the length scales from the
239  data
240  */
241  void auto_scale() {
242  scales.resize(nd_in);
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;
249  }
250  scales[i]=max-min;
251  }
252  return;
253  }
254 
255  /** \brief Initialize the data for the interpolation
256  for only one output function
257 
258  The object \c vecs should be a vector (of size <tt>n_in+1</tt>)
259  of vectors (all of size <tt>n_points</tt>). It may be
260  any type which allows the use of <tt>std::swap</tt> for
261  each vector in the list.
262  */
263  void set_data(size_t n_in, size_t n_points,
264  mat_t &dat, bool auto_scale=true) {
265  set_data(n_in,1,n_points,dat,auto_scale);
266  return;
267  }
268  //@}
269 
270  /// \name Evaluate interpolation
271  //@{
272  /** \brief Perform the interpolation over the first function
273  */
274  template<class vec2_t> double operator()(const vec2_t &x) const {
275  return eval(x);
276  }
277 
278  /** \brief Perform the interpolation over the first function
279  */
280  template<class vec2_t> double eval(const vec2_t &x) const {
281 
282  if (data_set==false) {
283  O2SCL_ERR("Data not set in interpm_idw::eval().",
284  exc_einval);
285  }
286 
287  // Compute distances
288  std::vector<double> dists(np);
289  for(size_t i=0;i<np;i++) {
290  dists[i]=dist(i,x);
291  }
292 
293  // Find closest points
294  std::vector<size_t> index;
295  o2scl::vector_smallest_index<std::vector<double>,double,
296  std::vector<size_t> >(dists,points+n_extra,index);
297 
298  if (n_extra>0) {
299  // Remove degenerate points to ensure accurate interpolation
300  bool found=true;
301  while (found==true) {
302  found=false;
303  // Find degenerate points and remove them
304  for(size_t j=0;j<points+n_extra;j++) {
305  for(size_t k=j;k<points+n_extra;k++) {
306  double dist_jk=dist(j,k);
307  if (index.size()>points && dist_jk<min_dist) {
308  found=true;
309  index.erase(index.begin()+j);
310  }
311  }
312  }
313  }
314  }
315 
316  // Check if the closest distance is zero
317  if (dists[index[0]]<=0.0) {
318  return data(nd_in,index[0]);
319  }
320 
321  // Compute normalization
322  double norm=0.0;
323  for(size_t i=0;i<points;i++) {
324  norm+=1.0/dists[index[i]];
325  }
326 
327  // Compute the inverse-distance weighted average
328  double ret=0.0;
329  for(size_t i=0;i<points;i++) {
330  ret+=data(nd_in,index[i])/dists[index[i]];
331  }
332  ret/=norm;
333 
334  // Return the average
335  return ret;
336  }
337 
338  /** \brief Perform the interpolation over the first function
339  with uncertainty
340  */
341  template<class vec2_t> void eval_err(const vec2_t &x, double &val,
342  double &err) const {
343 
344  if (data_set==false) {
345  O2SCL_ERR("Data not set in interpm_idw::eval_err().",
346  exc_einval);
347  }
348 
349  // Compute distances
350  std::vector<double> dists(np);
351  for(size_t i=0;i<np;i++) {
352  dists[i]=dist(i,x);
353  }
354 
355  // Find closest points
356  std::vector<size_t> index;
357  o2scl::vector_smallest_index<std::vector<double>,double,
358  std::vector<size_t> >(dists,points+1+n_extra,index);
359 
360  if (n_extra>0) {
361  // Remove degenerate points to ensure accurate interpolation
362  bool found=true;
363  while (found==true) {
364  found=false;
365  // Find degenerate points and remove them
366  for(size_t j=0;j<points+1+n_extra;j++) {
367  for(size_t k=j;k<points+1+n_extra;k++) {
368  double dist_jk=dist(j,k);
369  if (index.size()>points+1 && dist_jk<min_dist) {
370  found=true;
371  index.erase(index.begin()+j);
372  }
373  }
374  }
375  }
376  }
377 
378  if (dists[index[0]]<=0.0) {
379 
380  // If the closest distance is zero, just set the value
381  val=data(nd_in,index[0]);
382  err=0.0;
383  return;
384 
385  } else {
386 
387  std::vector<double> vals(points+1);
388 
389  for(size_t j=0;j<points+1;j++) {
390 
391  // Compute normalization
392  double norm=0.0;
393  for(size_t i=0;i<points+1;i++) {
394  if (i!=j) norm+=1.0/dists[index[i]];
395  }
396 
397  // Compute the inverse-distance weighted average
398  vals[j]=0.0;
399  for(size_t i=0;i<points+1;i++) {
400  if (i!=j) {
401  vals[j]+=data(nd_in,index[i])/dists[index[i]];
402  }
403  }
404  vals[j]/=norm;
405 
406  }
407 
408  val=vals[points];
409  err=o2scl::vector_stddev(vals);
410 
411  }
412 
413  return;
414  }
415 
416  /** \brief Perform the interpolation over all the functions,
417  at point \c x, storing the result in \c y
418  */
419  template<class vec2_t, class vec3_t>
420  void eval(const vec2_t &x, vec3_t &y) const {
421 
422  if (data_set==false) {
423  O2SCL_ERR("Data not set in interpm_idw::eval().",
424  exc_einval);
425  }
426 
427  if (verbose>0) {
428  std::cout << "interpm_idw: input: ";
429  for(size_t k=0;k<nd_in;k++) {
430  std::cout << x[k] << " ";
431  }
432  std::cout << std::endl;
433  }
434 
435  // Compute distances
436  std::vector<double> dists(np);
437  for(size_t i=0;i<np;i++) {
438  dists[i]=dist(i,x);
439  }
440 
441  // Find closest points
442  std::vector<size_t> index;
443  o2scl::vector_smallest_index<std::vector<double>,double,
444  std::vector<size_t> >(dists,points,index);
445  if (verbose>0) {
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]) << " ";
450  }
451  std::cout << std::endl;
452  }
453  }
454 
455  if (n_extra>0) {
456  // Remove degenerate points to ensure accurate interpolation
457  bool found=true;
458  while (found==true) {
459  found=false;
460  // Find degenerate points and remove them
461  for(size_t j=0;j<points+n_extra;j++) {
462  for(size_t k=j;k<points+n_extra;k++) {
463  double dist_jk=dist(j,k);
464  if (index.size()>points && dist_jk<min_dist) {
465  found=true;
466  index.erase(index.begin()+j);
467  }
468  }
469  }
470  }
471  }
472 
473  // Check if the closest distance is zero, if so, just
474  // return the value
475  if (dists[index[0]]<=0.0) {
476  for(size_t i=0;i<nd_out;i++) {
477  y[i]=data(nd_in+i,index[0]);
478  }
479  if (verbose>0) {
480  std::cout << "interpm_idw: distance zero. "
481  << "Returning values at index: " << index[0] << std::endl;
482  std::cout << "\t";
483  o2scl::vector_out(std::cout,nd_out,y,true);
484  }
485  return;
486  }
487 
488  // Compute normalization
489  double norm=0.0;
490  for(size_t i=0;i<points;i++) {
491  norm+=1.0/dists[index[i]];
492  }
493  if (verbose>0) {
494  std::cout << "interpm_idw: norm is " << norm << std::endl;
495  }
496 
497  // Compute the inverse-distance weighted averages
498  for(size_t j=0;j<nd_out;j++) {
499  y[j]=0.0;
500  for(size_t i=0;i<points;i++) {
501  if (j==0 && verbose>0) {
502  std::cout << "interpm_idw: Point: ";
503  for(size_t k=0;k<nd_in;k++) {
504  std::cout << data(k,index[i]) << " ";
505  }
506  std::cout << std::endl;
507  }
508  y[j]+=data(nd_in+j,index[i])/dists[index[i]];
509  if (verbose>0) {
510  std::cout << "interpm_idw: j,points,value,1/dist: "
511  << j << " " << i << " "
512  << data(nd_in+j,index[i]) << " "
513  << 1.0/dists[index[i]] << std::endl;
514  }
515  }
516  y[j]/=norm;
517  if (verbose>0) {
518  std::cout << "interpm_idw: y[" << j << "]: " << y[j]
519  << std::endl;
520  }
521  }
522 
523  return;
524  }
525 
526  /** \brief Perform the interpolation over all the functions
527  giving uncertainties and the sorted index vector
528 
529  The vector \c index is automatically resized to
530  a size equal to n_points+1+n_extra
531  be larger than
532  */
533  template<class vec2_t, class vec3_t, class vec4_t>
534  void eval_err_index(const vec2_t &x, vec3_t &val, vec4_t &err,
535  std::vector<size_t> &index) const {
536 
537  if (data_set==false) {
538  O2SCL_ERR("Data not set in interpm_idw::eval_err().",
539  exc_einval);
540  }
541 
542  // Compute distances
543  std::vector<double> dists(np);
544  for(size_t i=0;i<np;i++) {
545  dists[i]=dist(i,x);
546  }
547 
548  // Find closest points, note that index is automatically resized
549  // by the vector_smallest_index function
550  o2scl::vector_smallest_index<std::vector<double>,double,
551  std::vector<size_t> >(dists,points+1+n_extra,index);
552 
553  if (n_extra>0) {
554  // Remove degenerate points to ensure accurate interpolation
555  bool found=true;
556  while (found==true) {
557  found=false;
558  // Find degenerate points and remove them
559  for(size_t j=0;j<points+1+n_extra;j++) {
560  for(size_t k=j;k<points+1+n_extra;k++) {
561  double dist_jk=dist(j,k);
562  if (index.size()>points+1 && dist_jk<min_dist) {
563  found=true;
564  index.erase(index.begin()+j);
565  }
566  }
567  }
568  }
569  }
570 
571  if (dists[index[0]]<=0.0) {
572 
573  // If the closest distance is zero, just set the values and
574  // errors
575  for(size_t k=0;k<nd_out;k++) {
576  val[k]=data(nd_in+k,index[0]);
577  err[k]=0.0;
578  }
579  return;
580 
581  } else {
582 
583  for(size_t k=0;k<nd_out;k++) {
584 
585  std::vector<double> vals(points+1);
586 
587  for(size_t j=0;j<points+1;j++) {
588 
589  // Compute normalization
590  double norm=0.0;
591  for(size_t i=0;i<points+1;i++) {
592  if (i!=j) norm+=1.0/dists[index[i]];
593  }
594 
595  // Compute the inverse-distance weighted average
596  vals[j]=0.0;
597  for(size_t i=0;i<points+1;i++) {
598  if (i!=j) {
599  vals[j]+=data(nd_in+k,index[i])/dists[index[i]];
600  }
601  }
602  vals[j]/=norm;
603 
604  }
605 
606  // Instead of using the average, we report the value as the
607  // last element in the array, which is the interpolated
608  // value from the closest points
609  val[k]=vals[points];
610 
611  err[k]=o2scl::vector_stddev(vals);
612 
613  }
614 
615  }
616 
617  return;
618  }
619 
620  /** \brief Perform the interpolation over all the functions
621  with uncertainties
622  */
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;
626  return eval_err_index(x,val,err,index);
627  }
628  //@}
629 
630  /// \name Evaluate derivatives
631  //@{
632  /** \brief For one of the functions, compute the partial
633  derivatives (and uncertainties) with respect to all of the
634  inputs at one data point
635 
636  \note This function ignores the points chosen by \ref
637  set_points() and always chooses to average derivative
638  calculations determined from \c n_in+1 combinations of \c n_in
639  points .
640 
641  \todo Use the mechanism provided by <tt>n_extra</tt> above
642  to remove degenerate points.
643 
644  \future This function requires an extra copy from
645  "ders" to "ders2" which could be removed.
646  */
647  template<class vec3_t>
648  void derivs_err(size_t func_index, size_t point_index,
649  vec3_t &derivs, vec3_t &errs) const {
650 
651  if (data_set==false) {
652  O2SCL_ERR("Data not set in interpm_idw::derivs_err().",
653  exc_einval);
654  }
655 
656  // Set x equal to the specified point
657  ubvector x(nd_in);
658  for(size_t i=0;i<nd_in;i++) {
659  x[i]=data(i,point_index);
660  }
661  // Set f equal to the value of the function at the specified point
662  double f=data(nd_in+func_index,point_index);
663 
664  // The linear solver
666 
667  // Compute distances
668  std::vector<double> dists(np);
669  for(size_t i=0;i<np;i++) {
670  dists[i]=dist(i,x);
671  }
672 
673  // Find closest (but not identical) points
674 
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) {
679  O2SCL_ERR("Couldn't find enough nearby points.",o2scl::exc_einval);
680  }
681 
682  if (verbose>0) {
683  std::cout << "max_smallest: " << max_smallest << std::endl;
684  }
685 
686  o2scl::vector_smallest_index<std::vector<double>,double,
687  std::vector<size_t> >(dists,max_smallest,index);
688 
689  if (verbose>0) {
690  for(size_t i=0;i<index.size();i++) {
691  std::cout << "index[" << i << "] = " << index[i] << " "
692  << dists[index[i]] << std::endl;
693  }
694  }
695 
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;
701  }
702  }
703  if (index2.size()<nd_in+1) {
704  O2SCL_ERR("Couldn't find enough nearby points (2).",
706  }
707 
708  if (verbose>0) {
709  for(size_t i=0;i<index2.size();i++) {
710  std::cout << "index2[" << i << "] = " << index2[i] << " "
711  << dists[index2[i]] << std::endl;
712  }
713  }
714 
715  // Unit vector storage
716  std::vector<ubvector> units(nd_in+1);
717  // Difference vector norms
718  std::vector<double> diff_norms(nd_in+1);
719  // Storage for the derivative estimates
720  std::vector<ubvector> ders(nd_in+1);
721  // Matrix of unit vectors
722  ubmatrix m(nd_in,nd_in);
723  // Vector of function value differences
724  ubvector v(nd_in);
725  // Rearranged derivative object
726  std::vector<ubvector> ders2(nd_in);
727 
728  for(size_t i=0;i<nd_in+1;i++) {
729 
730  // Assign unit vector elements
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];
734  }
735 
736  // Normalize the unit vectors
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];
740  }
741 
742  }
743 
744  // Verbose output of the closest points and their norms
745  if (verbose>0) {
746  std::cout << "Point: ";
747  for(size_t i=0;i<nd_in;i++) {
748  std::cout << x[i] << " ";
749  }
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]) << " ";
755  }
756  std::cout << data(nd_in+func_index,index2[j]) << " "
757  << diff_norms[j] << std::endl;
758  }
759  for(size_t j=0;j<nd_in+1;j++) {
760  std::cout << "diff_norm: " << j << " " << diff_norms[j]
761  << std::endl;
762  }
763  // End of verbose output
764  }
765 
766  // Go through each set of points
767  for(size_t i=0;i<nd_in+1;i++) {
768 
769  ders[i].resize(nd_in);
770 
771  // Construct the matrix and vector for the solver
772  size_t jj=0;
773  for(size_t j=0;j<nd_in+1;j++) {
774  if (j!=i) {
775  for(size_t k=0;k<nd_in;k++) {
776  m(jj,k)=units[j][k];
777  }
778  v[jj]=(data(nd_in+func_index,index2[j])-f)/diff_norms[j];
779  jj++;
780  }
781  }
782 
783  // Solve to compute the derivatives
784  if (verbose>0) {
785  std::cout << "m:" << std::endl;
786  o2scl::matrix_out(std::cout,nd_in,nd_in,m);
787  std::cout << "v:" << std::endl;
788  o2scl::vector_out(std::cout,nd_in,v,true);
789  }
790  lshh.solve(nd_in,m,v,ders[i]);
791  if (verbose>0) {
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] << " ";
796  }
797  std::cout.unsetf(std::ios::showpos);
798  std::cout << std::endl;
799  }
800 
801  // Go to next derivative estimate
802  }
803 
804  for(size_t i=0;i<nd_in;i++) {
805 
806  // Rearrange derivatives
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];
810  }
811 
812  // Compute mean and standard deviation
813  derivs[i]=o2scl::vector_mean(ders2[i]);
814  errs[i]=o2scl::vector_stddev(ders2[i]);
815  }
816 
817  return;
818  }
819  //@}
820 
821 #ifndef DOXYGEN_INTERNAL
822 
823  protected:
824 
825  /// The number of points
826  size_t np;
827  /// The number of dimensions of the inputs
828  size_t nd_in;
829  /// The number of dimensions of the outputs
830  size_t nd_out;
831  /// The copy of the data
832  mat_t data;
833  /// True if the data has been specified
834  bool data_set;
835  /// Number of points to include in each interpolation (default 3)
836  size_t points;
837 
838  /// \name Distance determination [protected]
839  //@{
840  /// Distance scales for each coordinate
842 
843  /** \brief Compute the distance between \c x and the point at
844  index \c index
845  */
846  template<class vec2_t> double dist(size_t index,
847  const vec2_t &x) const {
848  double ret=0.0;
849  size_t nscales=scales.size();
850  for(size_t i=0;i<nd_in;i++) {
851  ret+=pow((x[i]-data(i,index))/scales[i%nscales],dist_expo);
852  }
853  return sqrt(ret);
854  }
855 
856  /** \brief Compute the distance between two points in the
857  data set
858  */
859  double dist(size_t j, size_t k) const {
860  double ret=0.0;
861  size_t nscales=scales.size();
862  for(size_t i=0;i<nd_in;i++) {
863  ret+=pow((data(i,j)-data(i,k))/scales[i%nscales],dist_expo);
864  }
865  return sqrt(ret);
866  }
867  //@}
868 
869 #endif
870 
871  };
872 
873 #ifndef DOXYGEN_NO_O2NS
874 }
875 #endif
876 
877 #endif
878 
boost::numeric::ublas::matrix< double >
o2scl::interpm_idw::data
mat_t data
The copy of the data.
Definition: interpm_idw.h:832
o2scl::interpm_idw::eval
void eval(const vec2_t &x, vec3_t &y) const
Perform the interpolation over all the functions, at point x, storing the result in y.
Definition: interpm_idw.h:420
boost::numeric::ublas::vector< double >
o2scl::interpm_idw::nd_out
size_t nd_out
The number of dimensions of the outputs.
Definition: interpm_idw.h:830
o2scl::matrix_out
int matrix_out(std::ostream &os, size_t nrows, size_t ncols, const mat_t &A)
A operator for simple matrix output using operator()
Definition: columnify.h:254
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl::interpm_idw
Multi-dimensional interpolation by inverse distance weighting.
Definition: interpm_idw.h:113
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::interpm_idw::nd_in
size_t nd_in
The number of dimensions of the inputs.
Definition: interpm_idw.h:828
o2scl::interpm_idw::get_data
void get_data(size_t &n_in, size_t &n_out, size_t &n_points, mat_t &dat)
Get the data used for interpolation.
Definition: interpm_idw.h:225
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::vector_mean
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
o2scl::interpm_idw::eval_err
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:624
o2scl::interpm_idw::data_set
bool data_set
True if the data has been specified.
Definition: interpm_idw.h:834
o2scl::vector_copy
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:124
o2scl::interpm_idw::set_data
void set_data(size_t n_in, size_t n_points, mat_t &dat, bool auto_scale=true)
Initialize the data for the interpolation for only one output function.
Definition: interpm_idw.h:263
o2scl::interpm_idw::points
size_t points
Number of points to include in each interpolation (default 3)
Definition: interpm_idw.h:836
o2scl::interpm_idw::eval_err_index
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:534
o2scl::interpm_idw::n_extra
size_t n_extra
The number of extra nearest neighbors to include to avoid degeneracies (default 0)
Definition: interpm_idw.h:145
o2scl::interpm_idw::set_points
void set_points(size_t n)
Set the number of closest points to use for each interpolation (default 3)
Definition: interpm_idw.h:157
o2scl::vector_stddev
double vector_stddev(size_t n, const vec_t &data)
Standard deviation with specified mean.
Definition: vec_stats.h:253
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::interpm_idw::eval_err
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:341
o2scl::interpm_idw::np
size_t np
The number of points.
Definition: interpm_idw.h:826
o2scl::interpm_idw::dist
double dist(size_t j, size_t k) const
Compute the distance between two points in the data set.
Definition: interpm_idw.h:859
o2scl::interpm_idw::eval
double eval(const vec2_t &x) const
Perform the interpolation over the first function.
Definition: interpm_idw.h:280
o2scl::interpm_idw::dist
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:846
o2scl_linalg::linear_solver_HH
Generic Householder linear solver.
Definition: linear_solver.h:109
o2scl::interpm_idw::set_data
void set_data(size_t n_in, size_t n_out, size_t n_points, mat_t &dat, bool auto_scale_flag=true)
Initialize the data for the interpolation.
Definition: interpm_idw.h:195
o2scl::interpm_idw::auto_scale
void auto_scale()
Automatically determine the length scales from the data.
Definition: interpm_idw.h:241
o2scl::interpm_idw::dist_expo
double dist_expo
Exponent in computing distance (default 2.0)
Definition: interpm_idw.h:136
o2scl::interpm_idw::min_dist
double min_dist
The minimum distance to consider points as non-degenerate (default )
Definition: interpm_idw.h:150
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::interpm_idw::scales
ubvector scales
Distance scales for each coordinate.
Definition: interpm_idw.h:841
o2scl::interpm_idw::verbose
int verbose
Verbosity parameter (default 0)
Definition: interpm_idw.h:140
o2scl_linalg::linear_solver_HH::solve
virtual void solve(size_t n, mat_t &A, vec_t &b, vec_t &x)
Solve square linear system of size n.
Definition: linear_solver.h:115
o2scl::vector_out
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:3145
o2scl::interpm_idw::derivs_err
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:648
o2scl::interpm_idw::set_scales
void set_scales(size_t n, vec2_t &v)
Set the scales for the distance metric.
Definition: interpm_idw.h:171
o2scl::interpm_idw::operator()
double operator()(const vec2_t &x) const
Perform the interpolation over the first function.
Definition: interpm_idw.h:274

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