interp2_planar.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_INTERP2_PLANAR_H
24 #define O2SCL_INTERP2_PLANAR_H
25 
26 /** \file interp2_planar.h
27  \brief File defining \ref o2scl::interp2_planar
28 */
29 
30 #include <iostream>
31 #include <string>
32 #include <cmath>
33 
34 #include <gsl/gsl_combination.h>
35 
36 #include <o2scl/err_hnd.h>
37 #include <o2scl/vector.h>
38 
39 #ifndef DOXYGEN_NO_O2NS
40 namespace o2scl {
41 #endif
42 
43  /** \brief Interpolate among two independent variables with planes
44 
45  \note This class is experimental.
46 
47  This class performs planar interpolation when the data points
48  are not arranged in a specified order (i.e. not on a grid). For
49  a set of data \f$ {x_i,y_i,f_i} \f$, the value of \f$ f \f$ is
50  predicted given a new value of x and y. This interpolation is
51  performed by finding the plane that goes through three closest
52  points in the data set. Distances are determined with
53  \f[
54  d_{ij} = \sqrt{\left(\frac{x_i-x_j}{\Delta x}\right)^2 +
55  \left(\frac{y_i-y_j}{\Delta y}\right)^2}
56  \f]
57  The values \f$ \Delta_x \f$ and \f$ \Delta_y \f$ are specified
58  in \ref x_scale and \ref y_scale, respectively. If these values
59  are negative (the default) then they are computed with \f$
60  \Delta x = x_{\mathrm{max}}-x_{\mathrm{min}} \f$ and \f$ \Delta
61  y = y_{\mathrm{max}}-y_{\mathrm{min}} \f$ .
62 
63  If the x- and y-values of the entire data set lie on a line,
64  then the interpolation will fail and the error handler will be
65  called. Colinearity is defined by a threshold, \ref thresh
66  which defaults to \f$ 10^{-12} \f$. If the denominator,
67  \f[
68  \sum_{k=1}^{3}\varepsilon_{ijk} x_i y_j < \mathrm{thresh}
69  \f]
70  where \f$ \varepsilon_{ijk} \f$ is an anti-symmetric Levi-Cevita
71  tensor, then the points are colinear. The value of \ref thresh
72  can be zero, but if it is negative then it will be reset
73  to the default value for the next interpolation.
74 
75  This class stores pointers to the data, not a copy. The
76  data can be changed between interpolations without an
77  additional call to \ref set_data(), but the scales may
78  need to be recomputed with \ref compute_scale().
79 
80  The vector type can be any type with a suitably defined \c
81  operator[].
82 
83  The interpolation requires at least three points and
84  \ref set_data() will call the error handler if the
85  first argument is less than three.
86 
87  \note This class operates by performing a \f$ {\cal O}(N) \f$
88  brute-force search to find the three closest points. If the
89  three closest points are colinear, then the data are sorted
90  by distance [ \f$ {\cal O}(N \log N) \f$ ], and the closest
91  triplets are enumerated until a non-colinear triplet is found.
92 
93  \note I believe this interpolation is a bit unstable because it
94  doesn't ensure that the user-specified objective point is inside
95  the region determined by the three closest data points, and this
96  can lead to some strong extrapolation.
97 
98  \future Make a parent class for this and \ref o2scl::interp2_neigh.
99  */
100  template<class vec_t> class interp2_planar {
101 
102  protected:
103 
104  /// The scale in the x direction
105  double dx;
106 
107  /// The scale in the y direction
108  double dy;
109 
110  public:
111 
114 
115  interp2_planar() {
116  data_set=false;
117  thresh=1.0e-12;
118  x_scale=-1.0;
119  y_scale=-1.0;
120  dx=0.0;
121  dy=0.0;
122  }
123 
124  /// Threshold for colinearity (default \f$ 10^{-12} \f$)
125  double thresh;
126 
127  /// The user-specified x scale (default -1)
128  double x_scale;
129 
130  /// The user-specified y scale (default -1)
131  double y_scale;
132 
133  /** \brief Initialize the data for the planar interpolation and
134  compute the scaling factors
135  */
136  void set_data(size_t n_points, vec_t &x, vec_t &y, vec_t &f) {
137  if (n_points<3) {
138  O2SCL_ERR2("Must provide at least three points in ",
139  "interp2_planar::set_data()",exc_efailed);
140  }
141  np=n_points;
142  ux=&x;
143  uy=&y;
144  uf=&f;
145  data_set=true;
146 
147  compute_scale();
148 
149  return;
150  }
151 
152  /// Find scaling
153  void compute_scale() {
154  if (x_scale<0.0) {
155  double minx=(*ux)[0], maxx=(*ux)[0];
156  for(size_t i=1;i<np;i++) {
157  if ((*ux)[i]<minx) minx=(*ux)[i];
158  if ((*ux)[i]>maxx) maxx=(*ux)[i];
159  }
160  dx=maxx-minx;
161  } else {
162  dx=x_scale;
163  }
164  if (y_scale<0.0) {
165  double miny=(*uy)[0], maxy=(*uy)[0];
166  for(size_t i=1;i<np;i++) {
167  if ((*uy)[i]<miny) miny=(*uy)[i];
168  if ((*uy)[i]>maxy) maxy=(*uy)[i];
169  }
170  dy=maxy-miny;
171  } else {
172  dy=y_scale;
173  }
174 
175  if (dx<=0.0 || dy<=0.0) {
176  O2SCL_ERR("No scale in interp2_planar::set_data().",exc_einval);
177  }
178 
179  return;
180  }
181 
182  /** \brief Perform the planar interpolation
183  */
184  double eval(double x, double y) const {
185  double x1, x2, x3, y1, y2, y3, f;
186  size_t i1, i2, i3;
187  eval_points(x,y,f,i1,x1,y1,i2,x2,y2,i3,x3,y3);
188  return f;
189  }
190 
191  /** \brief Perform the planar interpolation
192  */
193  double operator()(double x, double y) const {
194  return eval(x,y);
195  }
196 
197  /** \brief Perform the planar interpolation using the first two
198  elements of \c v as input
199  */
200  template<class vec2_t> double operator()(vec2_t &v) const {
201  return eval(v[0],v[1]);
202  }
203 
204  /** \brief Planar interpolation returning the closest points
205 
206  This function interpolates \c x and \c y into the data
207  returning \c f. It also returns the three closest x- and
208  y-values used for computing the plane.
209  */
210  void eval_points(double x, double y, double &f,
211  size_t &i1, double &x1, double &y1,
212  size_t &i2, double &x2, double &y2,
213  size_t &i3, double &x3, double &y3) const {
214 
215  if (data_set==false) {
216  O2SCL_ERR("Data not set in interp2_planar::eval_points().",
217  exc_einval);
218  }
219 
220  // First, we just find the three closest points by
221  // exhaustively searching the data
222 
223  // Put in initial points
224  i1=0; i2=1; i3=2;
225  double c1=sqrt(pow((x-(*ux)[0])/dx,2.0)+pow((y-(*uy)[0])/dy,2.0));
226  double c2=sqrt(pow((x-(*ux)[1])/dx,2.0)+pow((y-(*uy)[1])/dy,2.0));
227  double c3=sqrt(pow((x-(*ux)[2])/dx,2.0)+pow((y-(*uy)[2])/dy,2.0));
228 
229  // Sort initial points
230  if (c2<c1) {
231  if (c3<c2) {
232  // 321
233  swap(i1,c1,i3,c3);
234  } else if (c3<c1) {
235  // 231
236  swap(i1,c1,i2,c2);
237  swap(i2,c2,i3,c3);
238  } else {
239  // 213
240  swap(i1,c1,i2,c2);
241  }
242  } else {
243  if (c3<c1) {
244  // 312
245  swap(i1,c1,i3,c3);
246  swap(i2,c2,i3,c3);
247  } else if (c3<c2) {
248  // 132
249  swap(i3,c3,i2,c2);
250  }
251  // 123
252  }
253 
254  // Go through remaining points and sort accordingly
255  for(size_t j=3;j<np;j++) {
256  size_t i4=j;
257  double c4=sqrt(pow((x-(*ux)[i4])/dx,2.0)+pow((y-(*uy)[i4])/dy,2.0));
258  if (c4<c1) {
259  swap(i4,c4,i3,c3);
260  swap(i3,c3,i2,c2);
261  swap(i2,c2,i1,c1);
262  } else if (c4<c2) {
263  swap(i4,c4,i3,c3);
264  swap(i3,c3,i2,c2);
265  } else if (c4<c3) {
266  swap(i4,c4,i3,c3);
267  }
268  }
269 
270  // Solve for denominator:
271  double a, b, c, den, z1, z2, z3;
272  x1=(*ux)[i1];
273  x2=(*ux)[i2];
274  x3=(*ux)[i3];
275  y1=(*uy)[i1];
276  y2=(*uy)[i2];
277  y3=(*uy)[i3];
278  den=(-x2*y1+x3*y1+x1*y2-x3*y2-x1*y3+x2*y3);
279 
280  // If the points are colinear (if den is too small), then
281  // we do a complete sorting of the distance between the
282  // user-specified point and the data, and then go through
283  // all combinations to find the closest triplet that is
284  // not colinear.
285 
286  if (fabs(den)<thresh) {
287 
288  ubvector dist(np);
289  ubvector_size_t order(np);
290  for(size_t i=0;i<np;i++) {
291  dist[i]=sqrt(pow((x-(*ux)[i])/dx,2.0)+pow((y-(*uy)[i])/dy,2.0));
292  }
293  o2scl::vector_sort_index(np,dist,order);
294 
295  {
296  gsl_combination *cc=gsl_combination_alloc(np,3);
297 
298  gsl_combination_init_first(cc);
299 
300  bool done=false;
301  while (done==false) {
302  int status=gsl_combination_next(cc);
303 
304  if (status!=gsl_failure) {
305  i1=order[gsl_combination_get(cc,0)];
306  i2=order[gsl_combination_get(cc,1)];
307  i3=order[gsl_combination_get(cc,2)];
308  x1=(*ux)[i1];
309  x2=(*ux)[i2];
310  x3=(*ux)[i3];
311  y1=(*uy)[i1];
312  y2=(*uy)[i2];
313  y3=(*uy)[i3];
314  den=(-x2*y1+x3*y1+x1*y2-x3*y2-x1*y3+x2*y3);
315 
316  if (fabs(den)>thresh) {
317  done=true;
318  }
319 
320  } else {
321  /*
322  If we've gone through the entire data set and the
323  whole thing is colinear, then throw an exception.
324  */
325  O2SCL_ERR2("Interpolation failed in.",
326  "interp2_planar::eval_points().",exc_efailed);
327  }
328  }
329 
330  gsl_combination_free(cc);
331  }
332 
333  }
334 
335  // Calculate the function value with the three closest
336  // non-colinear points
337 
338  z1=(*uf)[i1];
339  z2=(*uf)[i2];
340  z3=(*uf)[i3];
341  a=-(-y2*z1+y3*z1+y1*z2-y3*z2-y1*z3+y2*z3)/den;
342  b=-(x2*z1-x3*z1-x1*z2+x3*z2+x1*z3-x2*z3)/den;
343  c=-(x3*y2*z1-x2*y3*z1-x3*y1*z2+x1*y3*z2+x2*y1*z3-x1*y2*z3)/den;
344 
345  f=a*x+b*y+c;
346 
347  return;
348  }
349 
350 #ifndef DOXYGEN_INTERNAL
351 
352  protected:
353 
354  /// The number of points
355  size_t np;
356  /// The x-values
357  vec_t *ux;
358  /// The y-values
359  vec_t *uy;
360  /// The f-values
361  vec_t *uf;
362  /// True if the data has been specified
363  bool data_set;
364 
365  /// Swap points 1 and 2.
366  int swap(size_t &index_1, double &dist_1, size_t &index_2,
367  double &dist_2) const {
368 
369  size_t index_temp;
370  double dist_temp;
371 
372  index_temp=index_1; dist_temp=dist_1;
373  index_1=index_2; dist_1=dist_2;
374  index_2=index_temp; dist_2=dist_temp;
375 
376  return 0;
377  }
378 
379  private:
380 
382  interp2_planar<vec_t>& operator=(const interp2_planar<vec_t>&);
383 
384 #endif
385 
386  };
387 
388 #ifndef DOXYGEN_NO_O2NS
389 }
390 #endif
391 
392 #endif
393 
394 
395 
o2scl::interp2_planar::x_scale
double x_scale
The user-specified x scale (default -1)
Definition: interp2_planar.h:128
o2scl::interp2_planar::np
size_t np
The number of points.
Definition: interp2_planar.h:355
boost::numeric::ublas::vector< double >
o2scl::interp2_planar::uf
vec_t * uf
The f-values.
Definition: interp2_planar.h:361
o2scl::interp2_planar
Interpolate among two independent variables with planes.
Definition: interp2_planar.h:100
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
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::interp2_planar::y_scale
double y_scale
The user-specified y scale (default -1)
Definition: interp2_planar.h:131
o2scl::interp2_planar::eval_points
void eval_points(double x, double y, double &f, size_t &i1, double &x1, double &y1, size_t &i2, double &x2, double &y2, size_t &i3, double &x3, double &y3) const
Planar interpolation returning the closest points.
Definition: interp2_planar.h:210
o2scl::interp2_planar::operator()
double operator()(double x, double y) const
Perform the planar interpolation.
Definition: interp2_planar.h:193
o2scl_inte_qng_coeffs::x2
static const double x2[5]
Definition: inte_qng_gsl.h:66
o2scl::interp2_planar::dx
double dx
The scale in the x direction.
Definition: interp2_planar.h:105
o2scl::vector_sort_index
void vector_sort_index(size_t n, const vec_t &data, vec_size_t &order)
Create a permutation which sorts the first n elements of a vector (in increasing order)
Definition: vector.h:800
o2scl_inte_qng_coeffs::x1
static const double x1[5]
Definition: inte_qng_gsl.h:48
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::interp2_planar::operator()
double operator()(vec2_t &v) const
Perform the planar interpolation using the first two elements of v as input.
Definition: interp2_planar.h:200
o2scl::interp2_planar::eval
double eval(double x, double y) const
Perform the planar interpolation.
Definition: interp2_planar.h:184
o2scl::interp2_planar::set_data
void set_data(size_t n_points, vec_t &x, vec_t &y, vec_t &f)
Initialize the data for the planar interpolation and compute the scaling factors.
Definition: interp2_planar.h:136
o2scl::interp2_planar::thresh
double thresh
Threshold for colinearity (default )
Definition: interp2_planar.h:125
o2scl::interp2_planar::compute_scale
void compute_scale()
Find scaling.
Definition: interp2_planar.h:153
o2scl::interp2_planar::data_set
bool data_set
True if the data has been specified.
Definition: interp2_planar.h:363
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::interp2_planar::dy
double dy
The scale in the y direction.
Definition: interp2_planar.h:108
o2scl::gsl_failure
@ gsl_failure
Failure.
Definition: err_hnd.h:49
o2scl::interp2_planar::swap
int swap(size_t &index_1, double &dist_1, size_t &index_2, double &dist_2) const
Swap points 1 and 2.
Definition: interp2_planar.h:366
o2scl::interp2_planar::uy
vec_t * uy
The y-values.
Definition: interp2_planar.h:359
o2scl::interp2_planar::ux
vec_t * ux
The x-values.
Definition: interp2_planar.h:357
o2scl_inte_qng_coeffs::x3
static const double x3[11]
Definition: inte_qng_gsl.h:94

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