interp.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 /* interpolation/linear.c
24  * interpolation/cspline.c
25  * interpolation/akima.c
26  * interpolation/steffen.c
27  *
28  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
29  * Copyright (C) 2014 Jean-Fran├žois Caron
30  *
31  * This program is free software; you can redistribute it and/or modify
32  * it under the terms of the GNU General Public License as published by
33  * the Free Software Foundation; either version 3 of the License, or (at
34  * your option) any later version.
35  *
36  * This program is distributed in the hope that it will be useful, but
37  * WITHOUT ANY WARRANTY; without even the implied warranty of
38  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
39  * General Public License for more details.
40  *
41  * You should have received a copy of the GNU General Public License
42  * along with this program; if not, write to the Free Software
43  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
44  * 02110-1301, USA.
45  */
46 #ifndef O2SCL_INTERP_H
47 #define O2SCL_INTERP_H
48 
49 /** \file interp.h
50  \brief One-dimensional interpolation classes and interpolation types
51 */
52 
53 #include <iostream>
54 #include <string>
55 
56 #include <boost/numeric/ublas/vector.hpp>
57 #include <boost/numeric/ublas/vector_proxy.hpp>
58 
59 #include <o2scl/search_vec.h>
60 #include <o2scl/tridiag.h>
61 
62 #ifndef DOXYGEN_NO_O2NS
63 namespace o2scl {
64 #endif
65 
66  /** \brief Interpolation types
67  */
68  enum {
69  /// Linear
71  /// Cubic spline for natural boundary conditions
73  /// Cubic spline for periodic boundary conditions
75  /// Akima spline for natural boundary conditions
77  /// Akima spline for periodic boundary conditions
79  /// Monotonicity-preserving interpolation (unfinished)
81  /// Steffen's monotonic method
83  };
84 
85  /** \brief Base low-level interpolation class [abstract base]
86 
87  See also the \ref intp_section section of the \o2 User's guide.
88 
89  To interpolate a set vectors \c x and \c y, call set() and then
90  the interpolation functions eval(), deriv(), deriv2() and
91  integ(). If the \c x and \c y vectors do not change, then you
92  may call the interpolation functions multiple times in
93  succession. These classes do not copy the user-specified vectors
94  but store pointers to them. Thus, if the vector is changed
95  without a new call to \ref interp_base::set(), the behavior of
96  the interpolation functions is undefined.
97 
98  \comment
99  AWS, 12/27/13: Copy constructors might be ill-advised for
100  this class since we store pointers. For now, we don't
101  allow the user to use them.
102  \endcomment
103  */
104  template<class vec_t, class vec2_t=vec_t> class interp_base {
105 
106 #ifdef O2SCL_NEVER_DEFINED
107  }{
108 #endif
109 
110 #ifndef DOXYGEN_INTERNAL
111 
112  protected:
113 
114  /** \brief To perform binary searches
115 
116  This pointer is set to zero in the constructor and should be
117  non-zero only if it has been allocated with \c new.
118  */
120 
121  /// Independent vector
122  const vec_t *px;
123 
124  /// Dependent vector
125  const vec2_t *py;
126 
127  /// Vector size
128  size_t sz;
129 
130  /** \brief An internal function to assist in computing the
131  integral for both the cspline and Akima types
132  */
133  double integ_eval(double ai, double bi, double ci, double di, double xi,
134  double a, double b) const {
135 
136  double r1=a-xi;
137  double r2=b-xi;
138  double r12=r1+r2;
139  double bterm=0.5*bi*r12;
140  double cterm=ci*(r1*r1+r2*r2+r1*r2)/3.0;
141  double dterm=0.25*di*r12*(r1*r1+r2*r2);
142 
143  return (b-a)*(ai+bterm+cterm+dterm);
144  }
145 
146 #endif
147 
148  public:
149 
150  interp_base() {
151  sz=0;
152  }
153 
154  virtual ~interp_base() {
155  }
156 
157  /** \brief The minimum size of the vectors to interpolate between
158 
159  This variable must be set in the constructor of the children
160  for access by the class user.
161  */
162  size_t min_size;
163 
164  /// Initialize interpolation routine
165  virtual void set(size_t size, const vec_t &x, const vec2_t &y)=0;
166 
167  /// Give the value of the function \f$ y(x=x_0) \f$ .
168  virtual double eval(double x0) const=0;
169 
170  /// Give the value of the function \f$ y(x=x_0) \f$ .
171  virtual double operator()(double x0) const {
172  return eval(x0);
173  }
174 
175  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
176  virtual double deriv(double x0) const=0;
177 
178  /** \brief Give the value of the second derivative
179  \f$ y^{\prime \prime}(x=x_0) \f$ .
180  */
181  virtual double deriv2(double x0) const=0;
182 
183  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
184  virtual double integ(double a, double b) const=0;
185 
186  /// Return the type
187  virtual const char *type() const=0;
188 
189 #ifndef DOXYGEN_INTERNAL
190 
191  private:
192 
195 
196 #endif
197 
198  };
199 
200  /** \brief Linear interpolation (GSL)
201 
202  See also the \ref intp_section section of the \o2 User's guide.
203 
204  Linear interpolation requires no calls to allocate() or free()
205  as there is no internal storage required.
206  */
207  template<class vec_t, class vec2_t=vec_t> class interp_linear :
208  public interp_base<vec_t,vec2_t> {
209 
210 #ifdef O2SCL_NEVER_DEFINED
211  }{
212 #endif
213 
214  public:
215 
216  interp_linear() {
217  this->min_size=2;
218  }
219 
220  virtual ~interp_linear() {}
221 
222  /// Initialize interpolation routine
223  virtual void set(size_t size, const vec_t &x, const vec2_t &y) {
224  if (size<this->min_size) {
225  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
226  " than "+szttos(this->min_size)+" in interp_linear::"+
227  "set().").c_str(),exc_einval);
228  }
229  this->svx.set_vec(size,x);
230  this->px=&x;
231  this->py=&y;
232  this->sz=size;
233  return;
234  }
235 
236  /// Give the value of the function \f$ y(x=x_0) \f$ .
237  virtual double eval(double x0) const {
238 
239  size_t index=this->svx.find(x0);
240 
241  double x_lo=(*this->px)[index];
242  double x_hi=(*this->px)[index+1];
243  double y_lo=(*this->py)[index];
244  double y_hi=(*this->py)[index+1];
245  double dx=x_hi-x_lo;
246 
247  return y_lo+(x0-x_lo)/dx*(y_hi-y_lo);
248  }
249 
250  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
251  virtual double deriv(double x0) const {
252 
253  size_t index=this->svx.find(x0);
254 
255  double x_lo=(*this->px)[index];
256  double x_hi=(*this->px)[index+1];
257  double y_lo=(*this->py)[index];
258  double y_hi=(*this->py)[index+1];
259  double dx=x_hi-x_lo;
260  double dy=y_hi-y_lo;
261 
262  return dy/dx;
263  }
264 
265  /** \brief Give the value of the second derivative
266  \f$ y^{\prime \prime}(x=x_0) \f$ (always zero)
267  */
268  virtual double deriv2(double x0) const {
269  return 0.0;
270  }
271 
272  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
273  virtual double integ(double a, double b) const {
274 
275  size_t i, index_a, index_b;
276 
277  bool flip=false;
278  if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
279  ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
280  double tmp=a;
281  a=b;
282  b=tmp;
283  flip=true;
284  }
285 
286  index_a=this->svx.find(a);
287  index_b=this->svx.find(b);
288 
289  double result=0.0;
290  for(i=index_a; i<=index_b; i++) {
291 
292  double x_lo=(*this->px)[i];
293  double x_hi=(*this->px)[i+1];
294  double y_lo=(*this->py)[i];
295  double y_hi=(*this->py)[i+1];
296  double dx=x_hi-x_lo;
297 
298  if(dx != 0.0) {
299 
300  if (i == index_a || i == index_b) {
301  double x1=(i == index_a) ? a : x_lo;
302  double x2=(i == index_b) ? b : x_hi;
303  double D=(y_hi-y_lo)/dx;
304  result += (x2-x1)*(y_lo+0.5*D*((x2-x_lo)+(x1-x_lo)));
305  } else {
306  result += 0.5*dx*(y_lo+y_hi);
307  }
308 
309  } else {
310  std::string str=((std::string)"Interval of length zero ")+
311  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
312  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
313  " in interp_linear::integ().";
314  O2SCL_ERR(str.c_str(),exc_einval);
315  }
316  }
317 
318  if (flip) result=-result;
319  return result;
320  }
321 
322  /// Return the type, \c "interp_linear".
323  virtual const char *type() const { return "interp_linear"; }
324 
325 #ifndef DOXYGEN_INTERNAL
326 
327  private:
328 
331 
332 #endif
333 
334  };
335 
336  /** \brief Cubic spline interpolation (GSL)
337 
338  See also the \ref intp_section section of the \o2 User's guide.
339 
340  By default, this class uses natural boundary conditions, where
341  the second derivative vanishes at each end point. Extrapolation
342  effectively assumes that the second derivative is linear outside
343  of the endpoints.
344  */
345  template<class vec_t, class vec2_t=vec_t> class interp_cspline :
346  public interp_base<vec_t,vec2_t> {
347 
348 #ifdef O2SCL_NEVER_DEFINED
349  }{
350 #endif
351 
352  public:
353 
355  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
356  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
357  typedef boost::numeric::ublas::slice slice;
358  typedef boost::numeric::ublas::range range;
359 
360 #ifndef DOXYGEN_INTERNAL
361 
362  protected:
363 
364  /// \name Storage for cubic spline interpolation
365  //@{
366  ubvector c;
367  ubvector g;
368  ubvector diag;
369  ubvector offdiag;
370  //@}
371 
372  /// Memory for the tridiagonalization
374 
375  /// Compute coefficients for cubic spline interpolation
376  void coeff_calc(const ubvector &c_array, double dy, double dx,
377  size_t index, double &b, double &c2, double &d) const {
378 
379  double c_i=c_array[index];
380  double c_ip1=c_array[index+1];
381  b=(dy/dx)-dx*(c_ip1+2.0*c_i)/3.0;
382  c2=c_i;
383  d=(c_ip1-c_i)/(3.0*dx);
384 
385  return;
386  }
387 
388 #endif
389 
390  public:
391 
392  /** \brief Create a base interpolation object with natural or
393  periodic boundary conditions
394  */
396  this->min_size=3;
397  }
398 
399  virtual ~interp_cspline() {
400  }
401 
402  /** \brief Initialize interpolation routine
403  */
404  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
405 
406  if (size<this->min_size) {
407  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
408  " than "+szttos(this->min_size)+" in interp_cspline::"+
409  "set().").c_str(),exc_einval);
410  }
411 
412  if (size!=this->sz) {
413  c.resize(size);
414  g.resize(size);
415  diag.resize(size);
416  offdiag.resize(size);
417  p4m.resize(size);
418  }
419 
420  this->px=&xa;
421  this->py=&ya;
422  this->sz=size;
423 
424  this->svx.set_vec(size,xa);
425 
426  /// Natural boundary conditions
427 
428  size_t i;
429  size_t num_points=size;
430  size_t max_index=num_points-1;
431  size_t sys_size=max_index-1;
432 
433  c[0]=0.0;
434  c[max_index]=0.0;
435 
436  for (i=0; i < sys_size; i++) {
437  double h_i=xa[i+1]-xa[i];
438  double h_ip1=xa[i+2]-xa[i+1];
439  double ydiff_i=ya[i+1]-ya[i];
440  double ydiff_ip1=ya[i+2]-ya[i+1];
441  double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
442  double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
443  offdiag[i]=h_ip1;
444  diag[i]=2.0*(h_ip1+h_i);
445  g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
446  }
447 
448  if (sys_size == 1) {
449 
450  c[1]=g[0]/diag[0];
451 
452  return;
453  }
454 
455  ubvector_range cp1(c,range(1,c.size()));
456  o2scl_linalg::solve_tridiag_sym<ubvector,ubvector,ubvector,
457  ubvector_range,o2scl_linalg::ubvector_4_mem,ubvector>
458  (diag,offdiag,g,cp1,sys_size,p4m);
459 
460  return;
461  }
462 
463  /// Give the value of the function \f$ y(x=x_0) \f$ .
464  virtual double eval(double x0) const {
465 
466  size_t index=this->svx.find(x0);
467 
468  double x_lo=(*this->px)[index];
469  double x_hi=(*this->px)[index+1];
470  double dx=x_hi-x_lo;
471 
472  double y_lo=(*this->py)[index];
473  double y_hi=(*this->py)[index+1];
474  double dy=y_hi-y_lo;
475  double delx=x0-x_lo;
476  double b_i, c_i, d_i;
477 
478  coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
479 
480  return y_lo+delx*(b_i+delx*(c_i+delx*d_i));
481  }
482 
483  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
484  virtual double deriv(double x0) const {
485 
486  size_t index=this->svx.find(x0);
487 
488  double x_lo=(*this->px)[index];
489  double x_hi=(*this->px)[index+1];
490  double dx=x_hi-x_lo;
491 
492  double y_lo=(*this->py)[index];
493  double y_hi=(*this->py)[index+1];
494  double dy=y_hi-y_lo;
495  double delx=x0-x_lo;
496  double b_i, c_i, d_i;
497 
498  coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
499 
500  return b_i+delx*(2.0*c_i+3.0*d_i*delx);
501  }
502 
503  /** \brief Give the value of the second derivative
504  \f$ y^{\prime \prime}(x=x_0) \f$ .
505  */
506  virtual double deriv2(double x0) const {
507 
508  size_t index=this->svx.find(x0);
509 
510  double x_lo=(*this->px)[index];
511  double x_hi=(*this->px)[index+1];
512  double dx=x_hi-x_lo;
513 
514  double y_lo=(*this->py)[index];
515  double y_hi=(*this->py)[index+1];
516  double dy=y_hi-y_lo;
517  double delx=x0-x_lo;
518  double b_i, c_i, d_i;
519 
520  coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
521 
522  return 2.0*c_i+6.0*d_i*delx;
523  }
524 
525  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
526  virtual double integ(double a, double b) const {
527 
528  size_t i, index_a, index_b;
529 
530  bool flip=false;
531  if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
532  ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
533  double tmp=a;
534  a=b;
535  b=tmp;
536  flip=true;
537  }
538 
539  index_a=this->svx.find(a);
540  index_b=this->svx.find(b);
541 
542  double result=0.0;
543 
544  for(i=index_a; i<=index_b; i++) {
545 
546  double x_lo=(*this->px)[i];
547  double x_hi=(*this->px)[i+1];
548  double y_lo=(*this->py)[i];
549  double y_hi=(*this->py)[i+1];
550  double dx=x_hi-x_lo;
551  double dy=y_hi-y_lo;
552 
553  if(dx != 0.0) {
554  double b_i, c_i, d_i;
555  coeff_calc(c,dy,dx,i,b_i,c_i,d_i);
556  if (i == index_a || i == index_b) {
557  double x1=(i == index_a) ? a : x_lo;
558  double x2=(i == index_b) ? b : x_hi;
559  result += this->integ_eval(y_lo,b_i,c_i,d_i,x_lo,x1,x2);
560  } else {
561  result += dx*(y_lo+dx*(0.5*b_i+
562  dx*(c_i/3.0+0.25*d_i*dx)));
563  }
564  } else {
565  std::string str=((std::string)"Interval of length zero ")+
566  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
567  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
568  " in interp_cspline::integ().";
569  O2SCL_ERR(str.c_str(),exc_einval);
570  }
571 
572  }
573 
574  if (flip) result*=-1.0;
575 
576  return result;
577  }
578 
579  /// Return the type, \c "interp_cspline".
580  virtual const char *type() const { return "interp_cspline"; }
581 
582 #ifndef DOXYGEN_INTERNAL
583 
584  private:
585 
588  (const interp_cspline<vec_t,vec2_t>&);
589 
590 #endif
591 
592  };
593 
594  /** \brief Cubic spline interpolation with periodic
595  boundary conditions (GSL)
596 
597  See also the \ref intp_section section of the \o2 User's guide.
598  */
599  template<class vec_t, class vec2_t=vec_t> class interp_cspline_peri :
600  public interp_cspline<vec_t,vec2_t> {
601 
602 #ifdef O2SCL_NEVER_DEFINED
603  }{
604 #endif
605 
606  protected:
607 
608  /// Memory for the tridiagonalization
610 
611  public:
612 
614  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
615  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
616  typedef boost::numeric::ublas::slice slice;
617  typedef boost::numeric::ublas::range range;
618 
619  interp_cspline_peri() : interp_cspline<vec_t,vec2_t>() {
620  this->min_size=2;
621  }
622 
623  virtual ~interp_cspline_peri() {
624  }
625 
626  /// Return the type, \c "interp_cspline_peri".
627  virtual const char *type() const { return "interp_cspline_peri"; }
628 
629  /** \brief Initialize interpolation routine
630  */
631  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
632 
633  if (size<this->min_size) {
634  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
635  " than "+szttos(this->min_size)+" in interp_cspline"+
636  "_peri::set().").c_str(),exc_einval);
637  }
638 
639  if (size!=this->sz) {
640  this->c.resize(size);
641  this->g.resize(size);
642  this->diag.resize(size);
643  this->offdiag.resize(size);
644  p5m.resize(size);
645  }
646 
647  this->px=&xa;
648  this->py=&ya;
649  this->sz=size;
650 
651  this->svx.set_vec(size,xa);
652 
653  /// Periodic boundary conditions
654 
655  size_t i;
656  size_t num_points=size;
657  // Engeln-Mullges+Uhlig "n"
658  size_t max_index=num_points-1;
659  // linear system is sys_size x sys_size
660  size_t sys_size=max_index;
661 
662  if (sys_size == 2) {
663 
664  // solve 2x2 system
665 
666  double h0=xa[1]-xa[0];
667  double h1=xa[2]-xa[1];
668  double h2=xa[3]-xa[2];
669  double A=2.0*(h0+h1);
670  double B=h0+h1;
671  double gx[2];
672  double det;
673 
674  gx[0]=3.0*((ya[2]-ya[1])/h1-(ya[1]-ya[0])/h0);
675  gx[1]=3.0*((ya[1]-ya[2])/h2-(ya[2]-ya[1])/h1);
676 
677  det=3.0*(h0+h1)*(h0+h1);
678  this->c[1]=( A*gx[0]-B*gx[1])/det;
679  this->c[2]=(-B*gx[0]+A*gx[1])/det;
680  this->c[0]=this->c[2];
681 
682  return;
683 
684  } else {
685 
686  for (i=0; i < sys_size-1; i++) {
687  double h_i=xa[i+1]-xa[i];
688  double h_ip1=xa[i+2]-xa[i+1];
689  double ydiff_i=ya[i+1]-ya[i];
690  double ydiff_ip1=ya[i+2]-ya[i+1];
691  double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
692  double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
693  this->offdiag[i]=h_ip1;
694  this->diag[i]=2.0*(h_ip1+h_i);
695  this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
696  }
697 
698  i=sys_size-1;
699  {
700  double h_i=xa[i+1]-xa[i];
701  double h_ip1=xa[1]-xa[0];
702  double ydiff_i=ya[i+1]-ya[i];
703  double ydiff_ip1=ya[1]-ya[0];
704  double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
705  double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
706  this->offdiag[i]=h_ip1;
707  this->diag[i]=2.0*(h_ip1+h_i);
708  this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
709  }
710 
711  ubvector_range cp1(this->c,range(1,this->c.size()));
712  o2scl_linalg::solve_cyc_tridiag_sym<ubvector,ubvector,ubvector,
713  ubvector_range,o2scl_linalg::ubvector_5_mem,ubvector>
714  (this->diag,this->offdiag,this->g,cp1,sys_size,p5m);
715  this->c[0]=this->c[max_index];
716 
717  return;
718  }
719 
720  }
721 
722 #ifndef DOXYGEN_INTERNAL
723 
724  private:
725 
729  (const interp_cspline_peri<vec_t,vec2_t>&);
730 
731 #endif
732 
733  };
734 
735  /** \brief Akima spline interpolation (GSL)
736 
737  See also the \ref intp_section section of the \o2 User's guide.
738 
739  This class uses natural boundary conditions, where the second
740  derivative vanishes at each end point. Extrapolation effectively
741  assumes that the second derivative is linear outside of the
742  endpoints. Use \ref interp_akima_peri for perodic boundary
743  conditions.
744  */
745  template<class vec_t, class vec2_t=vec_t> class interp_akima :
746  public interp_base<vec_t,vec2_t> {
747 
748  public:
749 
751  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
752  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
753  typedef boost::numeric::ublas::slice slice;
754  typedef boost::numeric::ublas::range range;
755 
756 #ifndef DOXYGEN_INTERNAL
757 
758  protected:
759 
760  /// \name Storage for Akima spline interpolation
761  //@{
762  ubvector b;
763  ubvector c;
764  ubvector d;
765  ubvector um;
766  //@}
767 
768  /// For initializing the interpolation
769  void akima_calc(const vec_t &x_array, size_t size,
770  ubvector &umx) {
771 
772  for(size_t i=0;i<this->sz-1;i++) {
773 
774  double NE=fabs(umx[3+i]-umx[2+i])+fabs(umx[1+i]-umx[i]);
775 
776  if (NE == 0.0) {
777  b[i]=umx[2+i];
778  c[i]=0.0;
779  d[i]=0.0;
780  } else {
781  double h_i=(*this->px)[i+1]-(*this->px)[i];
782  double NE_next=fabs(umx[4+i]-umx[3+i])+
783  fabs(umx[2+i]-umx[1+i]);
784  double alpha_i=fabs(umx[1+i]-umx[i])/NE;
785  double alpha_ip1;
786  double tL_ip1;
787  if (NE_next == 0.0) {
788  tL_ip1=umx[2+i];
789  } else {
790  alpha_ip1=fabs(umx[2+i]-umx[1+i])/NE_next;
791  tL_ip1=(1.0-alpha_ip1)*umx[2+i]+alpha_ip1*umx[3+i];
792  }
793  b[i]=(1.0-alpha_i)*umx[1+i]+alpha_i*umx[2+i];
794  c[i]=(3.0*umx[2+i]-2.0*b[i]-tL_ip1)/h_i;
795  d[i]=(b[i]+tL_ip1-2.0*umx[2+i])/(h_i*h_i);
796  }
797  }
798  }
799 
800 #endif
801 
802  public:
803 
804  /** \brief Create a base interpolation object with or without
805  periodic boundary conditions
806  */
808  this->min_size=5;
809  }
810 
811  virtual ~interp_akima() {
812  }
813 
814  /** \brief Initialize interpolation routine
815  */
816  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
817 
818  if (size<this->min_size) {
819  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
820  " than "+szttos(this->min_size)+" in interp_akima::"+
821  "set().").c_str(),exc_einval);
822  }
823 
824  if (size!=this->sz) {
825  b.resize(size);
826  c.resize(size);
827  d.resize(size);
828  um.resize(size+4);
829  }
830 
831  this->px=&xa;
832  this->py=&ya;
833  this->sz=size;
834 
835  this->svx.set_vec(size,xa);
836 
837  // Non-periodic boundary conditions
838 
839  ubvector_range m(um,range(2,um.size()));
840  size_t i;
841  for (i=0;i<=size-2;i++) {
842  m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
843  }
844 
845  um[0]=3.0*m[0]-2.0*m[1];
846  um[1]=2.0*m[0]-m[1];
847  m[this->sz-1]=2.0*m[size-2]-m[size-3];
848  m[size]=3.0*m[size-2]-2.0*m[size-3];
849 
850  akima_calc(xa,size,um);
851 
852  return;
853  }
854 
855  /// Give the value of the function \f$ y(x=x_0) \f$ .
856  virtual double eval(double x0) const {
857 
858  size_t index=this->svx.find(x0);
859 
860  double x_lo=(*this->px)[index];
861  double delx=x0-x_lo;
862  double bb=b[index];
863  double cc=c[index];
864  double dd=d[index];
865 
866  return (*this->py)[index]+delx*(bb+delx*(cc+dd*delx));
867  }
868 
869  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
870  virtual double deriv(double x0) const {
871 
872  size_t index=this->svx.find(x0);
873 
874  double x_lo=(*this->px)[index];
875  double delx=x0-x_lo;
876  double bb=b[index];
877  double cc=c[index];
878  double dd=d[index];
879 
880  return bb+delx*(2.0*cc+3.0*dd*delx);
881  }
882 
883  /** \brief Give the value of the second derivative
884  \f$ y^{\prime \prime}(x=x_0) \f$ .
885  */
886  virtual double deriv2(double x0) const {
887 
888  size_t index=this->svx.find(x0);
889 
890  double x_lo=(*this->px)[index];
891  double delx=x0-x_lo;
892  double cc=c[index];
893  double dd=d[index];
894 
895  return 2.0*cc+6.0*dd*delx;
896  }
897 
898  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
899  virtual double integ(double aa, double bb) const {
900 
901  size_t i, index_a, index_b;
902 
903  bool flip=false;
904  if (((*this->px)[0]<(*this->px)[this->sz-1] && aa>bb) ||
905  ((*this->px)[0]>(*this->px)[this->sz-1] && aa<bb)) {
906  double tmp=aa;
907  aa=bb;
908  bb=tmp;
909  flip=true;
910  }
911 
912  index_a=this->svx.find(aa);
913  index_b=this->svx.find(bb);
914 
915  double result=0.0;
916 
917  for(i=index_a; i<=index_b; i++) {
918 
919  double x_lo=(*this->px)[i];
920  double x_hi=(*this->px)[i+1];
921  double y_lo=(*this->py)[i];
922  double dx=x_hi-x_lo;
923 
924  if (dx != 0.0) {
925 
926  if (i==index_a || i==index_b) {
927  double x1=(i==index_a) ? aa : x_lo;
928  double x2=(i==index_b) ? bb : x_hi;
929  result += this->integ_eval(y_lo,b[i],c[i],d[i],x_lo,x1,x2);
930  } else {
931  result+=dx*(y_lo+dx*(0.5*b[i]+dx*(c[i]/3.0+0.25*d[i]*dx)));
932  }
933 
934  } else {
935  double y_hi=(*this->py)[i+1];
936  std::string str=((std::string)"Interval of length zero ")+
937  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
938  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
939  " in interp_akima::integ().";
940  O2SCL_ERR(str.c_str(),exc_einval);
941  }
942  }
943 
944  if (flip) result*=-1.0;
945 
946  return result;
947  }
948 
949  /// Return the type, \c "interp_akima".
950  virtual const char *type() const { return "interp_akima"; }
951 
952 #ifndef DOXYGEN_INTERNAL
953 
954  private:
955 
958 
959 #endif
960 
961  };
962 
963  /** \brief Akima spline interpolation with periodic
964  boundary conditions (GSL)
965 
966  See also the \ref intp_section section of the \o2 User's guide.
967  */
968  template<class vec_t, class vec2_t=vec_t> class interp_akima_peri :
969  public interp_akima<vec_t,vec2_t> {
970 
971  public:
972 
974  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
975  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
976  typedef boost::numeric::ublas::slice slice;
977  typedef boost::numeric::ublas::range range;
978 
979  public:
980 
982  }
983 
984  virtual ~interp_akima_peri() {
985  }
986 
987  /// Return the type, \c "interp_akima_peri".
988  virtual const char *type() const { return "interp_akima_peri"; }
989 
990  /// Initialize interpolation routine
991  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
992 
993  if (size<this->min_size) {
994  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
995  " than "+szttos(this->min_size)+" in interp_akima"+
996  "_peri::set().").c_str(),exc_einval);
997  }
998 
999  if (size!=this->sz) {
1000  this->b.resize(size);
1001  this->c.resize(size);
1002  this->d.resize(size);
1003  this->um.resize(size+4);
1004  }
1005 
1006  this->px=&xa;
1007  this->py=&ya;
1008  this->sz=size;
1009 
1010  this->svx.set_vec(size,xa);
1011 
1012  // Periodic boundary conditions
1013 
1014  ubvector_range m(this->um,range(2,this->um.size()));
1015 
1016  // Form the required set of divided differences
1017  for (size_t i=0;i<=size-2;i++) {
1018  m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
1019  }
1020 
1021  this->um[0]=m[this->sz-3];
1022  this->um[1]=m[this->sz-2];
1023  m[this->sz-1]=m[0];
1024  m[this->sz]=m[1];
1025 
1026  this->akima_calc(xa,size,this->um);
1027 
1028  return;
1029  }
1030 
1031 #ifndef DOXYGEN_INTERNAL
1032 
1033  private:
1034 
1037  (const interp_akima_peri<vec_t,vec2_t>&);
1038 
1039 #endif
1040 
1041  };
1042 
1043  /** \brief Steffen's monotonicity-preserving interpolation
1044 
1045  Adapted from the GSL version by J.-F. Caron which
1046  was based on \ref Steffen90 .
1047  */
1048  template<class vec_t, class vec2_t=vec_t> class interp_steffen :
1049  public interp_base<vec_t,vec2_t> {
1050 
1051 #ifdef O2SCL_NEVER_DEFINED
1052  }{
1053 #endif
1054 
1055  public:
1056 
1058  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
1059  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
1060  typedef boost::numeric::ublas::slice slice;
1061  typedef boost::numeric::ublas::range range;
1062 
1063 #ifndef DOXYGEN_INTERNAL
1064 
1065  protected:
1066 
1067  /// \name Storage for cubic spline interpolation
1068  //@{
1069  ubvector a;
1070  ubvector b;
1071  ubvector c;
1072  ubvector d;
1073  ubvector y_prime;
1074  //@}
1075 
1076  /** \brief Flip the sign of x if x and y have different signs
1077  */
1078  double copysign(const double x, const double y) {
1079  if ((x < 0 && y > 0) || (x > 0 && y < 0)) {
1080  return -x;
1081  }
1082  return x;
1083  }
1084 
1085 #endif
1086 
1087  public:
1088 
1089  /** \brief Create a base interpolation object */
1091  this->min_size=3;
1092  }
1093 
1094  virtual ~interp_steffen() {
1095  }
1096 
1097  /** \brief Initialize interpolation routine
1098  */
1099  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
1100 
1101  if (size<this->min_size) {
1102  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
1103  " than "+szttos(this->min_size)+" in interp_steffen::"+
1104  "set().").c_str(),exc_einval);
1105  }
1106 
1107  if (size!=this->sz) {
1108  a.resize(size);
1109  b.resize(size);
1110  c.resize(size);
1111  d.resize(size);
1112  y_prime.resize(size);
1113  }
1114 
1115  this->px=&xa;
1116  this->py=&ya;
1117  this->sz=size;
1118 
1119  this->svx.set_vec(size,xa);
1120 
1121  /*
1122  * first assign the interval and slopes for the left boundary.
1123  * We use the "simplest possibility" method described in the paper
1124  * in section 2.2
1125  */
1126  double h0=(xa[1]-xa[0]);
1127  double s0=(ya[1]-ya[0]) / h0;
1128 
1129  y_prime[0]=s0;
1130 
1131  /* Now we calculate all the necessary s, h, p, and y' variables
1132  from 1 to N-2 (0 to size-2 inclusive) */
1133  for (size_t i=1; i < (size-1); i++) {
1134 
1135  double pi;
1136 
1137  /* equation 6 in the paper */
1138  double hi=(xa[i+1]-xa[i]);
1139  double him1=(xa[i]-xa[i-1]);
1140 
1141  /* equation 7 in the paper */
1142  double si=(ya[i+1]-ya[i]) / hi;
1143  double sim1=(ya[i]-ya[i-1]) / him1;
1144 
1145  /* equation 8 in the paper */
1146  pi=(sim1*hi + si*him1) / (him1 + hi);
1147 
1148  /* This is a C equivalent of the FORTRAN statement below eqn 11 */
1149  y_prime[i]=(copysign(1.0,sim1)+copysign(1.0,si))*
1150  std::min(fabs(sim1),std::min(fabs(si),0.5*fabs(pi)));
1151 
1152  }
1153 
1154  /*
1155  * we also need y' for the rightmost boundary; we use the
1156  * "simplest possibility" method described in the paper in
1157  * section 2.2
1158  *
1159  * y'=s_{n-1}
1160  */
1161  y_prime[size-1]=(ya[size-1]-ya[size-2])/
1162  (xa[size-1]-xa[size-2]);
1163 
1164  /* Now we can calculate all the coefficients for the whole range. */
1165  for (size_t i=0; i < (size-1); i++) {
1166  double hi=(xa[i+1]-xa[i]);
1167  double si=(ya[i+1]-ya[i]) / hi;
1168 
1169  /* These are from equations 2-5 in the paper. */
1170  a[i]=(y_prime[i] + y_prime[i+1]-2*si) / hi / hi;
1171  b[i]=(3*si-2*y_prime[i]-y_prime[i+1]) / hi;
1172  c[i]=y_prime[i];
1173  d[i]=ya[i];
1174  }
1175 
1176  return;
1177  }
1178 
1179  /// Give the value of the function \f$ y(x=x_0) \f$ .
1180  virtual double eval(double x0) const {
1181 
1182  size_t index=this->svx.find(x0);
1183  double x_lo=(*this->px)[index];
1184  double delx=x0-x_lo;
1185 
1186  /* Use Horner's scheme for efficient evaluation of polynomials. */
1187  double y = d[index]+delx*(c[index]+delx*(b[index]+delx*a[index]));
1188 
1189  return y;
1190  }
1191 
1192  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1193  virtual double deriv(double x0) const {
1194 
1195  size_t index=this->svx.find(x0);
1196  double x_lo=(*this->px)[index];
1197  double delx=x0-x_lo;
1198 
1199  return c[index]+delx*(2.0*b[index]+delx*3.0*a[index]);
1200  }
1201 
1202  /** \brief Give the value of the second derivative
1203  \f$ y^{\prime \prime}(x=x_0) \f$ .
1204  */
1205  virtual double deriv2(double x0) const {
1206 
1207  size_t index=this->svx.find(x0);
1208  double x_lo=(*this->px)[index];
1209  double delx=x0-x_lo;
1210 
1211  return 2.0*b[index]+delx*6.0*a[index];
1212  }
1213 
1214  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1215  virtual double integ(double al, double bl) const {
1216 
1217  size_t i, index_a, index_b;
1218 
1219  bool flip=false;
1220  if (((*this->px)[0]<(*this->px)[this->sz-1] && al>bl) ||
1221  ((*this->px)[0]>(*this->px)[this->sz-1] && al<bl)) {
1222  double tmp=al;
1223  al=bl;
1224  bl=tmp;
1225  flip=true;
1226  }
1227 
1228  index_a=this->svx.find(al);
1229  index_b=this->svx.find(bl);
1230 
1231  double result=0.0;
1232 
1233  for(i=index_a; i<=index_b; i++) {
1234 
1235  double x_lo=(*this->px)[i];
1236  double x_hi=(*this->px)[i+1];
1237  double y_lo=(*this->py)[i];
1238  double y_hi=(*this->py)[i+1];
1239  double dx=x_hi-x_lo;
1240 
1241  if(dx != 0.0) {
1242 
1243  double x1=(i == index_a) ? al-x_lo : 0.0;
1244  double x2=(i == index_b) ? bl-x_lo : x_hi-x_lo;
1245  result += (1.0/4.0)*a[i]*(x2*x2*x2*x2-x1*x1*x1*x1)+
1246  (1.0/3.0)*b[i]*(x2*x2*x2-x1*x1*x1)+
1247  (1.0/2.0)*c[i]*(x2*x2-x1*x1)+d[i]*(x2-x1);
1248 
1249  } else {
1250  std::string str=((std::string)"Interval of length zero ")+
1251  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
1252  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
1253  " in interp_steffen::integ().";
1254  O2SCL_ERR(str.c_str(),exc_einval);
1255  }
1256 
1257  }
1258 
1259  if (flip) result*=-1.0;
1260 
1261  return result;
1262  }
1263 
1264  /// Return the type, \c "interp_steffen".
1265  virtual const char *type() const { return "interp_steffen"; }
1266 
1267 #ifndef DOXYGEN_INTERNAL
1268 
1269  private:
1270 
1272  interp_steffen<vec_t,vec2_t>& operator=
1273  (const interp_steffen<vec_t,vec2_t>&);
1274 
1275 #endif
1276 
1277  };
1278 
1279  /** \brief Monotonicity-preserving interpolation
1280 
1281  \warning This class is experimental. Integrals don't work yet.
1282 
1283  This class uses a method based on cubic Hermite interpolation,
1284  modifying the slopes to guarantee monotonicity. In the
1285  notation of \ref Fritsch80, if
1286  \f[
1287  \alpha_i^2+\beta_i^2 \geq 9
1288  \f]
1289  then \f$ \alpha \f$ and \f$ \beta \f$ are decreased by
1290  the factor \f$ \tau \f$ as described at the end of
1291  section 4.
1292 
1293  \note The results of the interpolation will only be monotonic in
1294  the regions where the original data set is also monotonic. Also,
1295  this interpolation routine does not enforce strict monotonicity,
1296  and the results of the interpolation will be flat where the data
1297  is also flat.
1298 
1299  Based on \ref Fritsch80 .
1300 
1301  \future Convert into fewer loops over the data
1302  */
1303  template<class vec_t, class vec2_t=vec_t> class interp_monotonic :
1304  public interp_base<vec_t,vec2_t> {
1305 
1306 #ifdef O2SCL_NEVER_DEFINED
1307  }{
1308 #endif
1309 
1310  public:
1311 
1313 
1314  protected:
1315 
1316  /// Slopes
1317  ubvector m;
1318  /// Finite differences
1319  ubvector Delta;
1320  /// Ratio
1321  ubvector alpha;
1322  /// Staggered ratio
1323  ubvector beta;
1324 
1325  public:
1326 
1327  interp_monotonic() {
1328  this->min_size=2;
1329  }
1330 
1331  virtual ~interp_monotonic() {
1332  }
1333 
1334  /// Initialize interpolation routine
1335  virtual void set(size_t size, const vec_t &x, const vec2_t &y) {
1336 
1337  // Verify size
1338  if (size<this->min_size) {
1339  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
1340  " than "+szttos(this->min_size)+" in interp_monotonic::"+
1341  "set().").c_str(),exc_einval);
1342  }
1343 
1344  // Setup search_vec object
1345  this->svx.set_vec(size,x);
1346 
1347  // Resize internal vectors
1348  if (this->sz!=size) {
1349  m.resize(size);
1350  Delta.resize(size-1);
1351  alpha.resize(size-1);
1352  beta.resize(size-1);
1353  }
1354 
1355  // Copy pointers
1356  this->px=&x;
1357  this->py=&y;
1358  this->sz=size;
1359 
1360  // Create the interpolation arrays in this sub-interval
1361 
1362  // Compute Delta and m
1363  for(size_t i=0;i<size-1;i++) {
1364  Delta[i]=(y[i+1]-y[i])/(x[i+1]-x[i]);
1365  if (i>0) {
1366  m[i]=(Delta[i]+Delta[i-1])/2.0;
1367  }
1368  }
1369  m[0]=Delta[0];
1370  m[size-1]=Delta[size-2];
1371 
1372  // Check to see if the data is flat anywhere
1373  for(size_t i=0;i<size-1;i++) {
1374  if (y[i]==y[i+1]) {
1375  m[i]=0.0;
1376  m[i+1]=0.0;
1377  }
1378  }
1379 
1380  // Compute alpha and beta
1381  for(size_t i=0;i<size-1;i++) {
1382  alpha[i]=m[i]/Delta[i];
1383  beta[i]=m[i+1]/Delta[i];
1384  }
1385 
1386  // Constrain m to ensure monotonicity
1387  for(size_t i=0;i<size-1;i++) {
1388  double norm2=alpha[i]*alpha[i]+beta[i]*beta[i];
1389  if (norm2>9.0) {
1390  double tau=3.0/sqrt(norm2);
1391  m[i]=tau*alpha[i]*Delta[i];
1392  m[i+1]=tau*beta[i]*Delta[i];
1393  }
1394  }
1395 
1396  return;
1397  }
1398 
1399  /// Give the value of the function \f$ y(x=x_0) \f$ .
1400  virtual double eval(double x0) const {
1401 
1402  size_t index=this->svx.find(x0);
1403 
1404  double x_lo=(*this->px)[index];
1405  double x_hi=(*this->px)[index+1];
1406  double y_lo=(*this->py)[index];
1407  double y_hi=(*this->py)[index+1];
1408  double h=x_hi-x_lo;
1409  double t=(x0-x_lo)/h;
1410  double t2=t*t, t3=t2*t;
1411 
1412  double h00=2.0*t3-3.0*t2+1.0;
1413  double h10=t3-2.0*t2+t;
1414  double h01=-2.0*t3+3.0*t2;
1415  double h11=t3-t2;
1416  double interp=y_lo*h00+h*m[index]*h10+y_hi*h01+h*m[index+1]*h11;
1417 
1418  return interp;
1419  }
1420 
1421  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1422  virtual double deriv(double x0) const {
1423 
1424  size_t index=this->svx.find(x0);
1425 
1426  double x_lo=(*this->px)[index];
1427  double x_hi=(*this->px)[index+1];
1428  double y_lo=(*this->py)[index];
1429  double y_hi=(*this->py)[index+1];
1430  double h=x_hi-x_lo;
1431  double t=(x0-x_lo)/h;
1432  double t2=t*t;
1433 
1434  double dh00=6.0*t2-6.0*t;
1435  double dh10=3.0*t2-4.0*t+1.0;
1436  double dh01=-6.0*t2+6.0*t;
1437  double dh11=3.0*t2-2.0*t;
1438  double deriv=(y_lo*dh00+h*m[index]*dh10+y_hi*dh01+
1439  h*m[index+1]*dh11)/h;
1440 
1441  return deriv;
1442  }
1443 
1444  /** \brief Give the value of the second derivative
1445  \f$ y^{\prime \prime}(x=x_0) \f$
1446  */
1447  virtual double deriv2(double x0) const {
1448 
1449  size_t index=this->svx.find(x0);
1450 
1451  double x_lo=(*this->px)[index];
1452  double x_hi=(*this->px)[index+1];
1453  double y_lo=(*this->py)[index];
1454  double y_hi=(*this->py)[index+1];
1455  double h=x_hi-x_lo;
1456  double t=(x0-x_lo)/h;
1457 
1458  double ddh00=12.0*t-6.0;
1459  double ddh10=6.0*t-4.0;
1460  double ddh01=-12.0*t+6.0;
1461  double ddh11=6.0*t-2.0;
1462  double deriv2=(y_lo*ddh00+h*m[index]*ddh10+y_hi*ddh01+
1463  h*m[index+1]*ddh11)/h/h;
1464 
1465  return deriv2;
1466  }
1467 
1468  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1469  virtual double integ(double a, double b) const {
1470 
1471  size_t i, index_a, index_b;
1472 
1473  bool flip=false;
1474  if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
1475  ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
1476  double tmp=a;
1477  a=b;
1478  b=tmp;
1479  flip=true;
1480  }
1481 
1482  index_a=this->svx.find(a);
1483  index_b=this->svx.find(b);
1484 
1485  double result=0.0;
1486 
1487  for(i=index_a; i<=index_b; i++) {
1488 
1489  double x_hi=(*this->px)[i+1];
1490  double x_lo=(*this->px)[i];
1491  double y_lo=(*this->py)[i];
1492  double y_hi=(*this->py)[i+1];
1493  double h=x_hi-x_lo;
1494 
1495  if (h != 0.0) {
1496 
1497  if (i == index_a) {
1498  x_lo=a;
1499  }
1500  if (i == index_b) {
1501  x_hi=b;
1502  }
1503 
1504  double t=(x_hi-x_lo)/h;
1505  double t2=t*t, t3=t2*t, t4=t3*t;
1506 
1507  double ih00=t4/2.0-t3+t;
1508  double ih10=t4/4.0-2.0*t3/3.0+t2/2.0;
1509  double ih01=-t4/2.0+t3;
1510  double ih11=t4/4.0-t3/3.0;
1511  double intres=h*(y_lo*ih00+h*m[i]*ih10+y_hi*ih01+
1512  h*m[i+1]*ih11);
1513  result+=intres;
1514 
1515  } else {
1516  std::string str=((std::string)"Interval of length zero ")+
1517  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
1518  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
1519  " in interp_monotonic::integ().";
1520  O2SCL_ERR(str.c_str(),exc_einval);
1521  }
1522 
1523  }
1524 
1525  if (flip) result*=-1.0;
1526 
1527  return result;
1528  }
1529 
1530  /// Return the type, \c "interp_monotonic".
1531  virtual const char *type() const { return "interp_monotonic"; }
1532 
1533 #ifndef DOXYGEN_INTERNAL
1534 
1535  private:
1536 
1539  (const interp_monotonic<vec_t,vec2_t>&);
1540 
1541 #endif
1542 
1543  };
1544 
1545  /** \brief Interpolation class for general vectors
1546 
1547  See also the \ref intp_section section of the \o2 User's guide.
1548 
1549  Interpolation of ublas vector like objects is performed with the
1550  default template parameters, and \ref interp_array is provided
1551  for simple interpolation on C-style arrays.
1552 
1553  The type of interpolation to be performed can be specified using
1554  the set_type() function or in the constructor. The default is
1555  cubic splines with natural boundary conditions.
1556  */
1557  template<class vec_t=boost::numeric::ublas::vector<double>,
1558  class vec2_t=vec_t> class interp {
1559 
1560 #ifndef DOXYGEN_INTERNAL
1561 
1562  protected:
1563 
1564  /// Pointer to base interpolation object
1566 
1567 #endif
1568 
1569  public:
1570 
1571  /// Create with base interpolation object \c it
1572  interp(size_t interp_type=itp_cspline) {
1573  if (interp_type==itp_linear) {
1575  } else if (interp_type==itp_cspline) {
1577  } else if (interp_type==itp_cspline_peri) {
1579  } else if (interp_type==itp_akima) {
1581  } else if (interp_type==itp_akima_peri) {
1583  } else if (interp_type==itp_monotonic) {
1585  } else if (interp_type==itp_steffen) {
1587  } else {
1588  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1589  o2scl::szttos(interp_type)+", in "+
1590  "interp::interp().").c_str(),exc_einval);
1591  }
1592  }
1593 
1594  virtual ~interp() {
1595  delete itp;
1596  }
1597 
1598  /// Give the value of the function \f$ y(x=x_0) \f$ .
1599  virtual double eval(const double x0, size_t n, const vec_t &x,
1600  const vec2_t &y) {
1601  itp->set(n,x,y);
1602  return itp->eval(x0);
1603  }
1604 
1605  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1606  virtual double deriv(const double x0, size_t n, const vec_t &x,
1607  const vec2_t &y) {
1608  itp->set(n,x,y);
1609  return itp->deriv(x0);
1610  }
1611 
1612  /** \brief Give the value of the second derivative
1613  \f$ y^{\prime \prime}(x=x_0) \f$ .
1614  */
1615  virtual double deriv2(const double x0, size_t n, const vec_t &x,
1616  const vec2_t &y) {
1617  itp->set(n,x,y);
1618  return itp->deriv2(x0);
1619  }
1620 
1621  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1622  virtual double integ(const double x1, const double x2, size_t n,
1623  const vec_t &x, const vec2_t &y) {
1624  itp->set(n,x,y);
1625  return itp->integ(x1,x2);
1626  }
1627 
1628  /// Set base interpolation type
1629  void set_type(size_t interp_type) {
1630  delete itp;
1631  if (interp_type==itp_linear) {
1633  } else if (interp_type==itp_cspline) {
1635  } else if (interp_type==itp_cspline_peri) {
1637  } else if (interp_type==itp_akima) {
1639  } else if (interp_type==itp_akima_peri) {
1641  } else if (interp_type==itp_monotonic) {
1643  } else if (interp_type==itp_steffen) {
1645  } else {
1646  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1647  o2scl::szttos(interp_type)+", in "+
1648  "interp::set().").c_str(),exc_einval);
1649  }
1650  return;
1651  }
1652 
1653 #ifndef DOXYGEN_INTERNAL
1654 
1655  private:
1656 
1658  interp<vec_t,vec2_t>& operator=(const interp<vec_t,vec2_t>&);
1659 
1660 #endif
1661 
1662  };
1663 
1664  /** \brief Interpolation class for pre-specified vector
1665 
1666  See also the \ref intp_section section of the \o2 User's guide.
1667 
1668  This interpolation class is intended to be sufficiently general
1669  to handle any vector type. Interpolation of ublas vector like
1670  objects is performed with the default template parameters.
1671 
1672  This class does not double check the vector to ensure that
1673  all of the intervals for the abcissa are all increasing or
1674  all decreasing.
1675 
1676  The type of interpolation to be performed can be specified
1677  using the set_type() function. The default is cubic splines
1678  with natural boundary conditions.
1679 
1680  \future Make version which copies vectors
1681  rather than storing pointers might be better and then
1682  has copy constructors.
1683  */
1684  template<class vec_t=boost::numeric::ublas::vector<double>,
1685  class vec2_t=vec_t> class interp_vec :
1686  public interp_base<vec_t,vec2_t> {
1687 
1688 #ifndef DOXYGEN_INTERNAL
1689 
1690  protected:
1691 
1692  /// Base interpolation object
1694 
1695  /// Interpolation type
1696  size_t itype;
1697 
1698 #endif
1699 
1700  public:
1701 
1702  /// Blank interpolator
1704  itp=0;
1705  itype=itp_cspline;
1706  }
1707 
1708  /// Create with base interpolation object \c it
1709  interp_vec(size_t n, const vec_t &x,
1710  const vec2_t &y, size_t interp_type=itp_cspline) {
1711 
1712  if (x[0]==x[n-1]) {
1713  O2SCL_ERR((((std::string)"Vector endpoints equal (value=")+
1714  o2scl::dtos(x[0])+") in interp_vec()::"+
1715  "interp_vec().").c_str(),exc_einval);
1716  }
1717 
1718  if (interp_type==itp_linear) {
1720  } else if (interp_type==itp_cspline) {
1722  } else if (interp_type==itp_cspline_peri) {
1724  } else if (interp_type==itp_akima) {
1726  } else if (interp_type==itp_akima_peri) {
1728  } else if (interp_type==itp_monotonic) {
1730  } else if (interp_type==itp_steffen) {
1732  } else {
1733  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1734  o2scl::szttos(interp_type)+", in "+
1735  "interp_vec::interp_vec().").c_str(),exc_einval);
1736  }
1737  itype=interp_type;
1738 
1739  itp->set(n,x,y);
1740  }
1741 
1742  virtual ~interp_vec() {
1743  if (itp!=0) {
1744  delete itp;
1745  }
1746  }
1747 
1748  /** \brief Set a new vector to interpolate
1749  */
1750  void set(size_t n, const vec_t &x, const vec2_t &y) {
1751  set(n,x,y,itype);
1752  return;
1753  }
1754 
1755  /** \brief Set a new vector to interpolate
1756  */
1757  void set(size_t n, const vec_t &x,
1758  const vec2_t &y, size_t interp_type) {
1759 
1760  if (x[0]==x[n-1]) {
1761  O2SCL_ERR((((std::string)"Vector endpoints equal (value=")+
1762  o2scl::dtos(x[0])+") in interp_vec()::"+
1763  "interp_vec().").c_str(),exc_einval);
1764  }
1765 
1766  delete itp;
1767  if (interp_type==itp_linear) {
1769  } else if (interp_type==itp_cspline) {
1771  } else if (interp_type==itp_cspline_peri) {
1773  } else if (interp_type==itp_akima) {
1775  } else if (interp_type==itp_akima_peri) {
1777  } else if (interp_type==itp_monotonic) {
1779  } else if (interp_type==itp_steffen) {
1781  } else {
1782  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1783  o2scl::szttos(interp_type)+", in "+
1784  "interp_vec::set().").c_str(),exc_einval);
1785  }
1786  itype=interp_type;
1787 
1788  itp->set(n,x,y);
1789  }
1790 
1791  /** \brief Manually clear the pointer to the user-specified vector
1792  */
1793  void clear() {
1794  if (itp!=0) {
1795  delete itp;
1796  itp=0;
1797  }
1798  return;
1799  }
1800 
1801  /// Give the value of the function \f$ y(x=x_0) \f$ .
1802  virtual double eval(const double x0) const {
1803  if (itp==0) {
1804  O2SCL_ERR("No vector set in interp_vec::eval().",
1805  exc_einval);
1806  }
1807  return itp->eval(x0);
1808  }
1809 
1810  /// Give the value of the function \f$ y(x=x_0) \f$ .
1811  virtual double operator()(double x0) const {
1812  if (itp==0) {
1813  O2SCL_ERR("No vector set in interp_vec::operator().",
1814  exc_einval);
1815  }
1816  return itp->eval(x0);
1817  }
1818 
1819  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1820  virtual double deriv(const double x0) const {
1821  if (itp==0) {
1822  O2SCL_ERR("No vector set in interp_vec::deriv().",
1823  exc_einval);
1824  }
1825  return itp->deriv(x0);
1826  }
1827 
1828  /** \brief Give the value of the second derivative
1829  \f$ y^{\prime \prime}(x=x_0) \f$ .
1830  */
1831  virtual double deriv2(const double x0) const {
1832  if (itp==0) {
1833  O2SCL_ERR("No vector set in interp_vec::deriv2().",
1834  exc_einval);
1835  }
1836  return itp->deriv2(x0);
1837  }
1838 
1839  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1840  virtual double integ(const double x1, const double x2) const {
1841  if (itp==0) {
1842  O2SCL_ERR("No vector set in interp_vec::integ().",
1843  exc_einval);
1844  }
1845  return itp->integ(x1,x2);
1846  }
1847 
1848  /// Return the type, "interp_vec"
1849  virtual const char *type() const {
1850  return "interp_vec";
1851  }
1852 
1853 #ifndef DOXYGEN_INTERNAL
1854 
1855  private:
1856 
1858  interp_vec<vec_t,vec2_t> &operator=
1859  (const interp_vec<vec_t,vec2_t> &it);
1860 
1861 #endif
1862 
1863  };
1864 
1865  /** \brief A specialization of interp for C-style double arrays
1866 
1867  See also the \ref intp_section section of the \o2 User's guide.
1868  */
1869  template<size_t n> class interp_array :
1870  public interp<double[n]> {
1871 
1872  public:
1873 
1874  /// Create with base interpolation objects \c it and \c rit
1875  interp_array(size_t interp_type)
1876  : interp<double[n]>(interp_type) {}
1877 
1878  /** \brief Create an interpolator using the default base
1879  interpolation objects
1880  */
1881  interp_array() : interp<double[n]>() {}
1882 
1883  };
1884 
1885  /** \brief A specialization of \ref o2scl::interp_vec for C-style arrays
1886 
1887  See also the \ref intp_section section of the \o2 User's guide.
1888  */
1889  template<class arr_t> class interp_array_vec :
1890  public interp_vec<arr_t> {
1891 
1892  public:
1893 
1894  /// Create with base interpolation object \c it
1895  interp_array_vec(size_t nv, const arr_t &x, const arr_t &y,
1896  size_t interp_type) :
1897  interp_vec<arr_t>(nv,x,y,interp_type) {}
1898  };
1899 
1900  /** \brief Count level crossings
1901 
1902  This function counts the number of times the function \f$ y(x) =
1903  \mathrm{level} \f$ where the function is defined by the vectors
1904  \c x and \c y.
1905 
1906  The value returned is exactly the same as the size of the
1907  \c locs vector computed by \ref vector_find_level().
1908 
1909  This function will call the error handler if \c n is
1910  less than two.
1911 
1912  If one of the entries in \c y is either larger or smaller than
1913  it's neighbors (i.e. if the function \f$ y(x) \f$ has an
1914  extremum), and if the value of \c level is exactly equal to the
1915  extremum, then this is counted as 1 level crossing. The same
1916  applies if either the first or last entry in \c y is exactly
1917  equal to \c level .
1918  */
1919  template<class vec_t, class vec2_t> size_t vector_level_count
1920  (double level, size_t n, vec_t &x, vec2_t &y) {
1921 
1922  if (n<=1) {
1923  O2SCL_ERR2("Need at least two data points in ",
1924  "vector_find_count().",exc_einval);
1925  }
1926 
1927  size_t count=0;
1928 
1929  // Handle left boundary
1930  if (y[0]==level) count++;
1931 
1932  // Find intersections
1933  for(size_t i=0;i<n-1;i++) {
1934 
1935  if ((y[i]<level && y[i+1]>=level) ||
1936  (y[i]>level && y[i+1]<=level)) {
1937  count++;
1938  }
1939  }
1940 
1941  return count;
1942  }
1943 
1944  /** \brief Perform inverse linear interpolation
1945 
1946  This function performs inverse linear interpolation of the data
1947  defined by \c x and \c y, finding all points in \c x which have
1948  the property \f$ \mathrm{level} = y(x) \f$. All points for which
1949  this relation holds are put into the vector \c locs. The
1950  previous information contained in vector \c locs before the
1951  function call is destroyed.
1952 
1953  This is the 1-dimensional analog of finding contour levels as
1954  the \ref contour class does for 2 dimensions.
1955 
1956  This function will call the error handler if \c n is
1957  less than two.
1958 
1959  This function is inclusive of the endpoints, so that if either
1960  <tt>y[0]</tt> and/or <tt>y[n-1]</tt> is equal to level, then
1961  <tt>x[0]</tt> and/or <tt>x[n-1]</tt> are added to \c locs,
1962  respectively.
1963  */
1964  template<class vec_t, class vec2_t> void vector_find_level
1965  (double level, size_t n, vec_t &x, vec2_t &y, std::vector<double> &locs) {
1966 
1967  if (n<=1) {
1968  O2SCL_ERR2("Need at least two data points in ",
1969  "vector_find_level().",exc_einval);
1970  }
1971 
1972  // Ensure that the location vector is empty
1973  locs.clear();
1974 
1975  // Handle left boundary
1976  if (y[0]==level) {
1977  locs.push_back(x[0]);
1978  }
1979 
1980  // Find intersections
1981  for(size_t i=0;i<n-1;i++) {
1982 
1983  if ((y[i]<level && y[i+1]>level) ||
1984  (y[i]>level && y[i+1]<level)) {
1985  // For each intersection, add the location using linear
1986  // interpolation
1987  double x0=x[i]+(x[i+1]-x[i])*(level-y[i])/(y[i+1]-y[i]);
1988  locs.push_back(x0);
1989  } else if (y[i+1]==level) {
1990  locs.push_back(x[i+1]);
1991  }
1992  }
1993 
1994  return;
1995  }
1996 
1997  /** \brief Compute the integral over <tt>y(x)</tt> using linear
1998  interpolation
1999 
2000  Assuming vectors \c y and \c x define a function \f$ y(x) \f$
2001  then this computes
2002  \f[
2003  \int_{x_0}^{x^{n-1}} y(x)~dx
2004  \f]
2005  using the trapezoidal rule implied by linear interpolation.
2006 
2007  \future It might be faster to compute the sum directly rather
2008  than going through an \ref o2scl::interp object .
2009  */
2010  template<class vec_t, class vec2_t>
2011  double vector_integ_linear(size_t n, vec_t &x, vec2_t &y) {
2012 
2013  // Interpolation object
2015 
2016  // Compute full integral
2017  double total=si.integ(x[0],x[n-1],n,x,y);
2018 
2019  return total;
2020  }
2021 
2022  /** \brief Compute the endpoints which enclose the regions whose
2023  integral is equal to \c sum
2024 
2025  Defining a new function, \f$ g(y_0) \f$ which takes as input any
2026  y-value, \f$ y_0 \f$ from the function \f$ y(x) \f$ (specified
2027  with the parameters \c x and \c y) and outputs the integral of
2028  the function \f$ y(x) \f$ over all regions where \f$ y(x) > y_0
2029  \f$. This function inverts \f$ g(y) \f$, taking the value
2030  of an integral as input, and returns the corresponding y-value
2031  in the variable \c lev.
2032 
2033  In order to make sure that the intepretation of the integral is
2034  unambiguous, this function requires that the first and last values
2035  of \c y are equal, i.e. <tt>y[0]==y[n-1]</tt>. This function
2036  also requires at least two data points, i.e. \c n cannot be
2037  0 or 1.
2038 
2039  This function is particularly useful, for example, in computing
2040  the region which defines 68\% around a peak of data, thus
2041  providing \f$ 1~\sigma \f$ confidence limits.
2042 
2043  Linear interpolation is used to describe the function \f$ g \f$,
2044  and the precision of this function is limited by this assumption.
2045  This function may also sometimes fail if \c sum is very close to
2046  the minimum or maximum value of the function \f$ g \f$.
2047 
2048  \todo This function may also require that all of the
2049  y-vector values have the same sign or are all positive.
2050  Check this.
2051 
2052  \comment
2053  Note that the two vector types for x and y must be the
2054  same in order to use o2scl_interp.
2055  \endcomment
2056  */
2057  template<class vec_t, class vec2_t> void vector_invert_enclosed_sum
2058  (double sum, size_t n, vec_t &x, vec2_t &y, double &lev,
2059  int verbose=0) {
2060 
2061  if (n<=1) {
2062  O2SCL_ERR2("Need at least two data points in ",
2063  "vector_invert_enclosed_sum().",exc_einval);
2064  }
2065 
2066  if (y[0]!=y[n-1]) {
2067  O2SCL_ERR2("The first and last y-values must be equal in ",
2068  "vector_invert_enclosed_sum().",exc_einval);
2069  }
2070 
2071  // Construct a sorted list of function values
2073  ubvector ysort(n);
2074  vector_copy(n,y,ysort);
2075  vector_sort<ubvector,double>(ysort.size(),ysort);
2076 
2077  // Create list of y-values to perform y-value and integral
2078  // interpolation. We choose values in between the grid points to
2079  // ensure that we don't accidentally land precisely on a peak or
2080  // valley.
2081  std::vector<double> ylist;
2082  for(size_t i=0;i<ysort.size()-1;i++) {
2083  if (ysort[i]!=ysort[i+1]) {
2084  ylist.push_back((ysort[i+1]+3.0*ysort[i])/4.0);
2085  ylist.push_back((ysort[i+1]*3.0+ysort[i])/4.0);
2086  }
2087  }
2088 
2089  // Interpolation objects. We need two, one for the
2090  // user-specified vector type, and one for std::vector<double>
2092  interp<std::vector<double>,std::vector<double> > itp2(itp_linear);
2093 
2094  // Construct vectors for interpolation
2095  std::vector<double> xi, yi;
2096 
2097  size_t nfail=0;
2098 
2099  for(size_t k=0;k<ylist.size();k++) {
2100  double lev_tmp=ylist[k];
2101  std::vector<double> locs;
2102  vector_find_level(lev_tmp,n,x,y,locs);
2103  if ((locs.size()%2)!=0) {
2104  nfail++;
2105  } else {
2106  double sum_temp=0.0;
2107  for(size_t i=0;i<locs.size()/2;i++) {
2108  double x0=locs[2*i];
2109  double x1=locs[2*i+1];
2110  sum_temp+=itp.integ(x0,x1,n,x,y);
2111  }
2112  xi.push_back(sum_temp);
2113  yi.push_back(lev_tmp);
2114  }
2115  }
2116  if (nfail>10) {
2117  O2SCL_ERR2("More than 10 failures in ",
2118  "vector_invert_enclosed_sum().",exc_einval);
2119  }
2120 
2121  if (verbose>1) {
2122  std::cout << "i, xi, yi: " << std::endl;
2123  for(size_t i=0;i<xi.size();i++) {
2124  std::cout << i << " " << xi[i] << " " << yi[i] << std::endl;
2125  }
2126  }
2127 
2128  lev=itp2.eval(sum,xi.size(),xi,yi);
2129 
2130  return;
2131  }
2132 
2133  /** \brief Find the region enclosing a partial integral
2134  */
2135  template<class vec_t, class vec2_t> void vector_region_parint
2136  (size_t n, vec_t &x, vec2_t &y, double frac, std::vector<double> &locs,
2137  int verbose=0) {
2138 
2139  if (frac<0.0 || frac>1.0) {
2140  O2SCL_ERR2("Fraction must be between 0 and 1 (exclusive) in ",
2141  "vector_region_parint().",exc_efailed);
2142  }
2143 
2144  // Total integral
2145  double total=vector_integ_linear(n,x,y);
2146  if (verbose>0) {
2147  std::cout << "Total integral: " << total << std::endl;
2148  }
2149  // Specified partial integral
2150  double partial=frac*total;
2151  if (verbose>0) {
2152  std::cout << "Partial integral: " << partial << std::endl;
2153  }
2154  // Find correct level
2155  double lev;
2156  vector_invert_enclosed_sum(partial,n,x,y,lev,verbose);
2157  if (verbose>0) {
2158  std::cout << "Level from vector_invert: " << lev << std::endl;
2159  }
2160  // Inverse interpolate to find locations corresponding to level
2161  vector_find_level(lev,n,x,y,locs);
2162  if (verbose>0) {
2163  std::cout << "Locations from vector_find_level: ";
2164  for(size_t i=0;i<locs.size();i++) {
2165  std::cout << locs[i];
2166  if (i!=locs.size()-1) {
2167  std::cout << " ";
2168  }
2169  }
2170  std::cout << std::endl;
2171  }
2172  return;
2173  }
2174 
2175  /** \brief Find the boundaries of the region enclosing a partial integral
2176  */
2177  template<class vec_t, class vec2_t> void vector_bound_parint
2178  (size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high) {
2179 
2180  std::vector<double> locs;
2181  vector_region_parint(n,x,y,frac,locs);
2182  if (locs.size()==0) {
2183  O2SCL_ERR("Zero level crossings in vector_bound_sigma().",
2184  exc_efailed);
2185  }
2186  // Return min and max location
2187  size_t ix;
2188  vector_min(locs.size(),locs,ix,low);
2189  vector_max(locs.size(),locs,ix,high);
2190  return;
2191  }
2192 
2193 #ifndef DOXYGEN_NO_O2NS
2194 }
2195 #endif
2196 
2197 #endif
void vector_region_parint(size_t n, vec_t &x, vec2_t &y, double frac, std::vector< double > &locs, int verbose=0)
Find the region enclosing a partial integral.
Definition: interp.h:2136
virtual const char * type() const
Return the type, "interp_steffen".
Definition: interp.h:1265
virtual const char * type() const
Return the type, "interp_cspline_peri".
Definition: interp.h:627
Interpolation class for pre-specified vector.
Definition: interp.h:1685
o2scl_linalg::ubvector_5_mem p5m
Memory for the tridiagonalization.
Definition: interp.h:609
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:237
virtual void set(size_t size, const vec_t &x, const vec2_t &y)=0
Initialize interpolation routine.
virtual double eval(const double x0, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the function .
Definition: interp.h:1599
interp_vec(size_t n, const vec_t &x, const vec2_t &y, size_t interp_type=itp_cspline)
Create with base interpolation object it.
Definition: interp.h:1709
virtual const char * type() const
Return the type, "interp_linear".
Definition: interp.h:323
double integ_eval(double ai, double bi, double ci, double di, double xi, double a, double b) const
An internal function to assist in computing the integral for both the cspline and Akima types...
Definition: interp.h:133
virtual double operator()(double x0) const
Give the value of the function .
Definition: interp.h:171
virtual const char * type() const =0
Return the type.
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
interp(size_t interp_type=itp_cspline)
Create with base interpolation object it.
Definition: interp.h:1572
void solve_tridiag_sym(const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N, mem_t &m)
Solve a symmetric tridiagonal linear system with user-specified memory.
Definition: tridiag_base.h:118
void akima_calc(const vec_t &x_array, size_t size, ubvector &umx)
For initializing the interpolation.
Definition: interp.h:769
size_t min_size
The minimum size of the vectors to interpolate between.
Definition: interp.h:162
A specialization of interp for C-style double arrays.
Definition: interp.h:1869
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:1422
const double pi
Definition: constants.h:62
interp_base< vec_t, vec2_t > * itp
Pointer to base interpolation object.
Definition: interp.h:1565
virtual double integ(double a, double b) const =0
Give the value of the integral .
virtual double deriv2(const double x0) const
Give the value of the second derivative .
Definition: interp.h:1831
void coeff_calc(const ubvector &c_array, double dy, double dx, size_t index, double &b, double &c2, double &d) const
Compute coefficients for cubic spline interpolation.
Definition: interp.h:376
Akima spline interpolation with periodic boundary conditions (GSL)
Definition: interp.h:968
virtual double integ(double al, double bl) const
Give the value of the integral .
Definition: interp.h:1215
void clear()
Manually clear the pointer to the user-specified vector.
Definition: interp.h:1793
virtual const char * type() const
Return the type, "interp_cspline".
Definition: interp.h:580
virtual const char * type() const
Return the type, "interp_akima".
Definition: interp.h:950
virtual const char * type() const
Return the type, "interp_monotonic".
Definition: interp.h:1531
double copysign(const double x, const double y)
Flip the sign of x if x and y have different signs.
Definition: interp.h:1078
size_t sz
Vector size.
Definition: interp.h:128
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:1400
invalid argument supplied by user
Definition: err_hnd.h:59
const vec_t * px
Independent vector.
Definition: interp.h:122
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:1180
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:870
Monotonicity-preserving interpolation.
Definition: interp.h:1303
Base low-level interpolation class [abstract base].
Definition: interp.h:104
interp_array_vec(size_t nv, const arr_t &x, const arr_t &y, size_t interp_type)
Create with base interpolation object it.
Definition: interp.h:1895
virtual const char * type() const
Return the type, "interp_akima_peri".
Definition: interp.h:988
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:1193
virtual double deriv(const double x0) const
Give the value of the derivative .
Definition: interp.h:1820
Akima spline for periodic boundary conditions.
Definition: interp.h:78
void vector_min(size_t n, const vec_t &data, size_t &index, data_t &val)
Compute the minimum of the first n elements of a vector.
Definition: vector.h:1208
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:526
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:1447
interp_vec()
Blank interpolator.
Definition: interp.h:1703
virtual double deriv(double x0) const =0
Give the value of the derivative .
generic failure
Definition: err_hnd.h:61
virtual double integ(double aa, double bb) const
Give the value of the integral .
Definition: interp.h:899
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:1205
Steffen&#39;s monotonic method.
Definition: interp.h:82
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:886
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:1469
Steffen&#39;s monotonicity-preserving interpolation.
Definition: interp.h:1048
interp_akima()
Create a base interpolation object with or without periodic boundary conditions.
Definition: interp.h:807
Cubic spline interpolation with periodic boundary conditions (GSL)
Definition: interp.h:599
virtual double operator()(double x0) const
Give the value of the function .
Definition: interp.h:1811
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:484
virtual double integ(const double x1, const double x2, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the integral .
Definition: interp.h:1622
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:273
Cubic spline for natural boundary conditions.
Definition: interp.h:72
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:251
void vector_bound_parint(size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high)
Find the boundaries of the region enclosing a partial integral.
Definition: interp.h:2178
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
ubvector beta
Staggered ratio.
Definition: interp.h:1323
Allocation object for 5 arrays of equal size.
Definition: tridiag_base.h:86
size_t find(const double x0) const
Search an increasing or decreasing vector for the interval containing x0
Definition: search_vec.h:137
virtual const char * type() const
Return the type, "interp_vec".
Definition: interp.h:1849
interp_array()
Create an interpolator using the default base interpolation objects.
Definition: interp.h:1881
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
size_t itype
Interpolation type.
Definition: interp.h:1696
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
virtual double integ(const double x1, const double x2) const
Give the value of the integral .
Definition: interp.h:1840
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:856
const vec2_t * py
Dependent vector.
Definition: interp.h:125
o2scl_linalg::ubvector_4_mem p4m
Memory for the tridiagonalization.
Definition: interp.h:373
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
Allocation object for 4 arrays of equal size.
Definition: tridiag_base.h:70
void set_type(size_t interp_type)
Set base interpolation type.
Definition: interp.h:1629
size_t vector_level_count(double level, size_t n, vec_t &x, vec2_t &y)
Count level crossings.
Definition: interp.h:1920
Akima spline for natural boundary conditions.
Definition: interp.h:76
double vector_integ_linear(size_t n, vec_t &x, vec2_t &y)
Compute the integral over y(x) using linear interpolation.
Definition: interp.h:2011
Akima spline interpolation (GSL)
Definition: interp.h:745
void vector_invert_enclosed_sum(double sum, size_t n, vec_t &x, vec2_t &y, double &lev, int verbose=0)
Compute the endpoints which enclose the regions whose integral is equal to sum.
Definition: interp.h:2058
interp_steffen()
Create a base interpolation object.
Definition: interp.h:1090
virtual double eval(double x0) const =0
Give the value of the function .
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:506
void vector_max(size_t n, const vec_t &data, size_t &index, data_t &val)
Compute the maximum of the first n elements of a vector.
Definition: vector.h:1133
interp_cspline()
Create a base interpolation object with natural or periodic boundary conditions.
Definition: interp.h:395
Cubic spline for periodic boundary conditions.
Definition: interp.h:74
Monotonicity-preserving interpolation (unfinished)
Definition: interp.h:80
virtual double deriv2(const double x0, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the second derivative .
Definition: interp.h:1615
ubvector Delta
Finite differences.
Definition: interp.h:1319
Cubic spline interpolation (GSL)
Definition: interp.h:345
interp_array(size_t interp_type)
Create with base interpolation objects it and rit.
Definition: interp.h:1875
virtual double deriv2(double x0) const =0
Give the value of the second derivative .
void solve_cyc_tridiag_sym(const vec_t &diag, const vec2_t &offdiag, const vec3_t &b, vec4_t &x, size_t N, mem_t &m)
Solve a symmetric cyclic tridiagonal linear system with user specified memory.
Definition: tridiag_base.h:241
void set_vec(size_t nn, const vec_t &x)
Set the vector to be searched.
Definition: search_vec.h:120
search_vec< const vec_t > svx
To perform binary searches.
Definition: interp.h:119
ubvector alpha
Ratio.
Definition: interp.h:1321
static const double x2[5]
Definition: inte_qng_gsl.h:66
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
Linear interpolation (GSL)
Definition: interp.h:207
std::string szttos(size_t x)
Convert a size_t to a string.
virtual double deriv2(double x0) const
Give the value of the second derivative (always zero)
Definition: interp.h:268
void vector_find_level(double level, size_t n, vec_t &x, vec2_t &y, std::vector< double > &locs)
Perform inverse linear interpolation.
Definition: interp.h:1965
Linear.
Definition: interp.h:70
ubvector m
Slopes.
Definition: interp.h:1317
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:464
A specialization of o2scl::interp_vec for C-style arrays.
Definition: interp.h:1889
interp_base< vec_t, vec2_t > * itp
Base interpolation object.
Definition: interp.h:1693
Interpolation class for general vectors.
Definition: interp.h:1558
virtual double deriv(const double x0, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the derivative .
Definition: interp.h:1606

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