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

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