interp2_seq.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 /** \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  double result;
354  ubvector icol(nx);
355  for(size_t i=0;i<nx;i++) {
356  if (iy==-1) {
357  icol[i]=itps[i]->integ(y0,y1);
358  } else if (iy==0) {
359  icol[i]=itps[i]->eval(y0);
360  } else if (iy==1) {
361  icol[i]=itps[i]->deriv(y0);
362  } else if (iy==2) {
363  icol[i]=itps[i]->deriv2(y0);
364  } else {
365  O2SCL_ERR2("Invalid value of 'iy' for interp2_seq::",
366  "interp_gen(). (xfirst=true)",exc_einval);
367  }
368  }
369  interp_vec<vec_t,ubvector> six(nx,*xfun,icol,itype);
370  if (ix==-1) {
371  result=six.integ(x0,x1);
372  } else if (ix==0) {
373  result=six.eval(x0);
374  } else if (ix==1) {
375  result=six.deriv(x0);
376  } else if (ix==2) {
377  result=six.deriv2(x0);
378  } else {
379  O2SCL_ERR2("Invalid value of 'ix' for interp2_seq::",
380  "interp_gen(). (xfirst=true)",exc_einval);
381  }
382  return result;
383  }
384 
385 #ifndef DOXYGEN_INTERNAL
386 
387  protected:
388 
389  /// The array of interpolation objects
390  std::vector<interp_vec<vec_t,mat_row_t> *> itps;
391 
392  /// An array of rows
393  std::vector<mat_row_t *> vecs;
394 
395  /// The number of x grid points
396  size_t nx;
397 
398  /// The number of y grid points
399  size_t ny;
400 
401  /// True if the data has been specified by the user
402  bool data_set;
403 
404  /// The x grid
405  vec_t *xfun;
406 
407  /// The y grid
408  vec_t *yfun;
409 
410  /// The data
411  mat_t *datap;
412 
413  /// Interpolation type
414  size_t itype;
415 
416  private:
417 
421  (const interp2_seq<vec_t,mat_t,mat_row_t>&);
422 
423 #endif
424 
425  };
426 
427 #ifndef DOXYGEN_NO_O2NS
428 }
429 #endif
430 
431 #endif
432 
433 
434 
size_t ny
The number of y grid points.
Definition: interp2_seq.h:399
Interpolation class for pre-specified vector.
Definition: interp.h:1685
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:414
std::vector< mat_row_t * > vecs
An array of rows.
Definition: interp2_seq.h:393
virtual double deriv2(const double x0) const
Give the value of the second derivative .
Definition: interp.h:1831
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:390
vec_t * yfun
The y grid.
Definition: interp2_seq.h:408
vec_t * xfun
The x grid.
Definition: interp2_seq.h:405
mat_t * datap
The data.
Definition: interp2_seq.h:411
virtual double deriv(const double x0) const
Give the value of the derivative .
Definition: interp.h:1820
generic failure
Definition: err_hnd.h:61
Cubic spline for natural boundary conditions.
Definition: interp.h:72
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:1840
bool data_set
True if the data has been specified by the user.
Definition: interp2_seq.h:402
#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:396
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:1802
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).