interp2_seq.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 /** \file interp2_seq.h
24  \brief File defining \ref o2scl::interp2_seq
25 */
26 #ifndef O2SCL_INTERP2_SEQ_H
27 #define O2SCL_INTERP2_SEQ_H
28 
29 #include <iostream>
30 #include <string>
31 #include <vector>
32 
33 #include <boost/numeric/ublas/vector.hpp>
34 #include <boost/numeric/ublas/matrix.hpp>
35 #include <boost/numeric/ublas/matrix_proxy.hpp>
36 
37 #include <o2scl/interp.h>
38 #include <o2scl/interp2.h>
39 
40 #ifndef DOXYGEN_NO_O2NS
41 namespace o2scl {
42 #endif
43 
44  /** \brief Two-dimensional interpolation class by successive
45  one-dimensional interpolation
46 
47  This class implements two-dimensional interpolation by iterating
48  the \o2 one-dimensional interpolation routines. Derivatives and
49  integrals along both x- and y-directions can be computed. This
50  class is likely a bit slower than \ref interp2_direct but more
51  flexible.
52 
53  The convention used by this class is that the first (row) index
54  of the matrix enumerates the x coordinate and that the second
55  (column) index enumerates the y coordinate. See the discussion
56  in the User's guide in the section called \ref rowcol_subsect.
57 
58  The function set_data() does not copy the data, it stores
59  pointers to the data. If the data is modified, then the function
60  reset_interp() must be called to reset the interpolation
61  information with the original pointer information. The storage
62  for the data, including the arrays \c x_grid and \c y_grid are
63  all managed by the user.
64 
65  By default, cubic spline interpolation with natural boundary
66  conditions is used. This can be changed by calling set_interp()
67  again with the same data and the new interpolation type.
68 
69  There is an example for the usage of this class given
70  in <tt>examples/ex_interp2_seq.cpp</tt>.
71 
72  Because of the way this class creates pointers to the
73  data, copy construction is not currently allowed.
74 
75  \future Implement an improved caching system in case, for example
76  \c xfirst is true and the last interpolation used the same
77  value of \c x.
78  */
79  template<class vec_t=boost::numeric::ublas::vector<double>,
80  class mat_t=boost::numeric::ublas::matrix<double>,
81  class mat_row_t=boost::numeric::ublas::matrix_row<mat_t> >
82  class interp2_seq : public interp2_base<vec_t,mat_t> {
83 
84 #ifdef O2SCL_NEVER_DEFINED
85  }{
86 #endif
87 
88  public:
89 
91 
92  interp2_seq() {
93  data_set=false;
95  }
96 
97  virtual ~interp2_seq() {
98  for(size_t i=0;i<itps.size();i++) {
99  delete itps[i];
100  delete vecs[i];
101  }
102  itps.clear();
103  }
104 
105  /** \brief Initialize the data for the 2-dimensional interpolation
106 
107  If \c x_first is true, then set_data() creates interpolation
108  objects for each of the rows. Calls to interp() then uses
109  these to create a column at the specified value of \c x. An
110  interpolation object is created at this column to find the
111  value of the function at the specified value \c y. If \c
112  x_first is false, the opposite strategy is employed. These two
113  options may give slightly different results.
114  */
115  void set_data(size_t n_x, size_t n_y, vec_t &x_grid,
116  vec_t &y_grid, mat_t &data,
117  size_t interp_type=itp_cspline) {
118 
119  // Set new data
120  itype=interp_type;
121  nx=n_x;
122  ny=n_y;
123  xfun=&x_grid;
124  yfun=&y_grid;
125  datap=&data;
126 
127  // Set interpolation objects
128 
129  for(size_t i=0;i<itps.size();i++) {
130  delete itps[i];
131  delete vecs[i];
132  }
133  itps.clear();
134 
135  // If we interpolate along the x-axis first, then we want to fix the
136  // first index, to get nx rows of size ny
137  vecs.resize(nx);
138  itps.resize(nx);
139  for(size_t i=0;i<nx;i++) {
140  vecs[i]=new mat_row_t
141  (o2scl::matrix_row<mat_t,mat_row_t>(*datap,i));
143  }
144  data_set=true;
145 
146  return;
147  }
148 
149  /** \brief Reset the stored interpolation since the data has changed
150 
151  This will throw an exception if the set_data() has not been
152  called.
153  */
154  void reset_interp() {
155  if (data_set) {
156 
157  for(size_t i=0;i<itps.size();i++) {
158  delete itps[i];
159  delete vecs[i];
160  }
161  itps.clear();
162 
163  // Set interpolation objects
164  vecs.resize(nx);
165  itps.resize(nx);
166  for(size_t i=0;i<nx;i++) {
167  vecs[i]=new mat_row_t
168  (o2scl::matrix_row<mat_t,mat_row_t>(*datap,i));
170  }
171 
172  } else {
173  O2SCL_ERR("Data not set in interp2_seq::reset_interp().",exc_einval);
174  }
175 
176  return;
177  }
178 
179  /** \brief Perform the 2-d interpolation
180  */
181  double eval(double x, double y) const {
182  if (data_set==false) {
183  O2SCL_ERR("Data not set in interp2_seq::eval().",exc_efailed);
184  }
185  double result;
186  ubvector icol(nx);
187  for(size_t i=0;i<nx;i++) {
188  icol[i]=itps[i]->eval(y);
189  }
190  interp_vec<vec_t,ubvector> six(nx,*xfun,icol,itype);
191  result=six.eval(x);
192  return result;
193  }
194 
195  /** \brief Perform the 2-d interpolation
196  */
197  double operator()(double x, double y) const {
198  return eval(x,y);
199  }
200 
201  /** \brief Compute the partial derivative in the x-direction
202  */
203  double deriv_x(double x, double y) const {
204  if (data_set==false) {
205  O2SCL_ERR("Data not set in interp2_seq::eval().",exc_efailed);
206  return 0.0;
207  }
208  double result;
209  ubvector icol(nx);
210  for(size_t i=0;i<nx;i++) {
211  icol[i]=itps[i]->eval(y);
212  }
213  interp_vec<vec_t,ubvector> six(nx,*xfun,icol,itype);
214  result=six.deriv(x);
215  return result;
216  }
217 
218  /** \brief Compute the partial second derivative in the x-direction
219  */
220  double deriv_xx(double x, double y) const {
221  if (data_set==false) {
222  O2SCL_ERR("Data not set in interp2_seq::deriv_xx().",exc_efailed);
223  return 0.0;
224  }
225  double result;
226  ubvector icol(nx);
227  for(size_t i=0;i<nx;i++) {
228  icol[i]=itps[i]->eval(y);
229  }
230  interp_vec<vec_t,ubvector> six(nx,*xfun,icol,itype);
231  result=six.deriv2(x);
232  return result;
233  }
234 
235  /** \brief Compute the integral in the x-direction between x=x0
236  and x=x1
237  */
238  double integ_x(double x0, double x1, double y) const {
239  if (data_set==false) {
240  O2SCL_ERR("Data not set in interp2_seq::integ_x().",exc_efailed);
241  return 0.0;
242  }
243  double result;
244  ubvector icol(nx);
245  for(size_t i=0;i<nx;i++) {
246  icol[i]=itps[i]->eval(y);
247  }
248  interp_vec<vec_t,ubvector> six(nx,*xfun,icol,itype);
249  result=six.integ(x0,x1);
250  return result;
251 
252  }
253 
254  /** \brief Compute the partial derivative in the y-direction
255  */
256  double deriv_y(double x, double y) const {
257  if (data_set==false) {
258  O2SCL_ERR("Data not set in interp2_seq::deriv_y().",exc_efailed);
259  return 0.0;
260  }
261  double result;
262  ubvector icol(nx);
263  for(size_t i=0;i<nx;i++) {
264  icol[i]=itps[i]->deriv(y);
265  }
266  interp_vec<vec_t,ubvector> six(nx,*xfun,icol,itype);
267  result=six.eval(x);
268  return result;
269  }
270 
271  /** \brief Compute the partial second derivative in the y-direction
272  */
273  double deriv_yy(double x, double y) const {
274  if (data_set==false) {
275  O2SCL_ERR("Data not set in interp2_seq::deriv_yy().",exc_efailed);
276  return 0.0;
277  }
278  double result;
279  ubvector icol(nx);
280  for(size_t i=0;i<nx;i++) {
281  icol[i]=itps[i]->deriv2(y);
282  }
283  interp_vec<vec_t,ubvector> six(nx,*xfun,icol,itype);
284  result=six.eval(x);
285  return result;
286  }
287 
288  /** \brief Compute the integral in the y-direction between y=y0
289  and y=y1
290  */
291  double integ_y(double x, double y0, double y1) const {
292  if (data_set==false) {
293  O2SCL_ERR("Data not set in interp2_seq::integ_y().",exc_efailed);
294  return 0.0;
295  }
296  double result;
297  ubvector icol(nx);
298  for(size_t i=0;i<nx;i++) {
299  icol[i]=itps[i]->integ(y0,y1);
300  }
301  interp_vec<vec_t,ubvector> six(nx,*xfun,icol,itype);
302  result=six.eval(x);
303  return result;
304  }
305 
306  /** \brief Compute the mixed partial derivative
307  \f$ \frac{\partial^2 f}{\partial x \partial y} \f$
308  */
309  double deriv_xy(double x, double y) const {
310  if (data_set==false) {
311  O2SCL_ERR("Data not set in interp2_seq::eval().",exc_efailed);
312  return 0.0;
313  }
314  double result;
315  ubvector icol(nx);
316  for(size_t i=0;i<nx;i++) {
317  icol[i]=itps[i]->deriv(y);
318  }
319  interp_vec<vec_t,ubvector> six(nx,*xfun,icol,itype);
320  result=six.deriv(x);
321  return result;
322  }
323 
324  /** \brief Compute a general interpolation result
325 
326  This computes
327  \f[
328  \frac{\partial^m}{\partial x^m}
329  \frac{\partial^n}{\partial y^n} f(x,y)
330  \f]
331  for \f$ m \in (-1,0,1,2) \f$ and \f$ n \in (-1,0,1,2) \f$ with
332  the notation
333  \f{eqnarray*}
334  \frac{\partial^{-1}}{\partial x^{-1}}
335  &\equiv & \int_{x_0}^{x_1} f~dx \nonumber \\
336  \frac{\partial^0}{\partial x^0} &\equiv &
337  \left.f\right|_{x=x_0} \nonumber \\
338  \frac{\partial^1}{\partial x^1} &\equiv &
339  \left(\frac{\partial f}{\partial x}\right)_{x=x_0} \nonumber \\
340  \frac{\partial^2}{\partial x^2} &\equiv &
341  \left(\frac{\partial^2 f}{\partial x^2}\right)_{x=x_0}
342  \f}
343  and the value of \f$ x_1 \f$ is ignored when \f$ m \geq 0 \f$
344  and the value of \f$ y_1 \f$ is ignored when \f$ n \geq 0 \f$.
345 
346  */
347  double eval_gen(int ix, int iy, double x0, double x1,
348  double y0, double y1) const {
349  if (data_set==false) {
350  O2SCL_ERR("Data not set in interp2_seq::interp_gen().",exc_efailed);
351  return 0.0;
352  }
353  // Initialize to prevent un'inited var. warnings
354  double result=0.0;
355  ubvector icol(nx);
356  for(size_t i=0;i<nx;i++) {
357  if (iy==-1) {
358  icol[i]=itps[i]->integ(y0,y1);
359  } else if (iy==0) {
360  icol[i]=itps[i]->eval(y0);
361  } else if (iy==1) {
362  icol[i]=itps[i]->deriv(y0);
363  } else if (iy==2) {
364  icol[i]=itps[i]->deriv2(y0);
365  } else {
366  O2SCL_ERR2("Invalid value of 'iy' for interp2_seq::",
367  "interp_gen(). (xfirst=true)",exc_einval);
368  }
369  }
370  interp_vec<vec_t,ubvector> six(nx,*xfun,icol,itype);
371  if (ix==-1) {
372  result=six.integ(x0,x1);
373  } else if (ix==0) {
374  result=six.eval(x0);
375  } else if (ix==1) {
376  result=six.deriv(x0);
377  } else if (ix==2) {
378  result=six.deriv2(x0);
379  } else {
380  O2SCL_ERR2("Invalid value of 'ix' for interp2_seq::",
381  "interp_gen(). (xfirst=true)",exc_einval);
382  }
383  return result;
384  }
385 
386 #ifndef DOXYGEN_INTERNAL
387 
388  protected:
389 
390  /// The array of interpolation objects
391  std::vector<interp_vec<vec_t,mat_row_t> *> itps;
392 
393  /// An array of rows
394  std::vector<mat_row_t *> vecs;
395 
396  /// The number of x grid points
397  size_t nx;
398 
399  /// The number of y grid points
400  size_t ny;
401 
402  /// True if the data has been specified by the user
403  bool data_set;
404 
405  /// The x grid
406  vec_t *xfun;
407 
408  /// The y grid
409  vec_t *yfun;
410 
411  /// The data
412  mat_t *datap;
413 
414  /// Interpolation type
415  size_t itype;
416 
417  private:
418 
422  (const interp2_seq<vec_t,mat_t,mat_row_t>&);
423 
424 #endif
425 
426  };
427 
428 #ifndef DOXYGEN_NO_O2NS
429 }
430 #endif
431 
432 #endif
433 
434 
435 
size_t ny
The number of y grid points.
Definition: interp2_seq.h:400
Interpolation class for pre-specified vector.
Definition: interp.h:1795
double deriv_yy(double x, double y) const
Compute the partial second derivative in the y-direction.
Definition: interp2_seq.h:273
double integ_x(double x0, double x1, double y) const
Compute the integral in the x-direction between x=x0 and x=x1.
Definition: interp2_seq.h:238
Two-dimensional interpolation class by successive one-dimensional interpolation.
Definition: interp2_seq.h:82
double deriv_x(double x, double y) const
Compute the partial derivative in the x-direction.
Definition: interp2_seq.h:203
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
size_t itype
Interpolation type.
Definition: interp2_seq.h:415
std::vector< mat_row_t * > vecs
An array of rows.
Definition: interp2_seq.h:394
virtual double deriv2(const double x0) const
Give the value of the second derivative .
Definition: interp.h:1949
invalid argument supplied by user
Definition: err_hnd.h:59
std::vector< interp_vec< vec_t, mat_row_t > * > itps
The array of interpolation objects.
Definition: interp2_seq.h:391
vec_t * yfun
The y grid.
Definition: interp2_seq.h:409
vec_t * xfun
The x grid.
Definition: interp2_seq.h:406
mat_t * datap
The data.
Definition: interp2_seq.h:412
virtual double deriv(const double x0) const
Give the value of the derivative .
Definition: interp.h:1938
generic failure
Definition: err_hnd.h:61
Cubic spline for natural boundary conditions.
Definition: interp.h:73
double deriv_y(double x, double y) const
Compute the partial derivative in the y-direction.
Definition: interp2_seq.h:256
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
virtual double integ(const double x1, const double x2) const
Give the value of the integral .
Definition: interp.h:1958
bool data_set
True if the data has been specified by the user.
Definition: interp2_seq.h:403
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
double operator()(double x, double y) const
Perform the 2-d interpolation.
Definition: interp2_seq.h:197
double eval(double x, double y) const
Perform the 2-d interpolation.
Definition: interp2_seq.h:181
Two-dimensional interpolation base class [abstract].
Definition: interp2.h:40
double eval_gen(int ix, int iy, double x0, double x1, double y0, double y1) const
Compute a general interpolation result.
Definition: interp2_seq.h:347
double deriv_xy(double x, double y) const
Compute the mixed partial derivative .
Definition: interp2_seq.h:309
size_t nx
The number of x grid points.
Definition: interp2_seq.h:397
void reset_interp()
Reset the stored interpolation since the data has changed.
Definition: interp2_seq.h:154
double integ_y(double x, double y0, double y1) const
Compute the integral in the y-direction between y=y0 and y=y1.
Definition: interp2_seq.h:291
double deriv_xx(double x, double y) const
Compute the partial second derivative in the x-direction.
Definition: interp2_seq.h:220
virtual double eval(const double x0) const
Give the value of the function .
Definition: interp.h:1920
static const double x1[5]
Definition: inte_qng_gsl.h:48
void set_data(size_t n_x, size_t n_y, vec_t &x_grid, vec_t &y_grid, mat_t &data, size_t interp_type=itp_cspline)
Initialize the data for the 2-dimensional interpolation.
Definition: interp2_seq.h:115

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