interp.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 /* 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 cache=0;
240  size_t index=this->svx.find_const(x0,cache);
241 
242  double x_lo=(*this->px)[index];
243  double x_hi=(*this->px)[index+1];
244  double y_lo=(*this->py)[index];
245  double y_hi=(*this->py)[index+1];
246  double dx=x_hi-x_lo;
247 
248  return y_lo+(x0-x_lo)/dx*(y_hi-y_lo);
249  }
250 
251  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
252  virtual double deriv(double x0) const {
253 
254  size_t cache=0;
255  size_t index=this->svx.find_const(x0,cache);
256 
257  double x_lo=(*this->px)[index];
258  double x_hi=(*this->px)[index+1];
259  double y_lo=(*this->py)[index];
260  double y_hi=(*this->py)[index+1];
261  double dx=x_hi-x_lo;
262  double dy=y_hi-y_lo;
263 
264  return dy/dx;
265  }
266 
267  /** \brief Give the value of the second derivative
268  \f$ y^{\prime \prime}(x=x_0) \f$ (always zero)
269  */
270  virtual double deriv2(double x0) const {
271  return 0.0;
272  }
273 
274  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
275  virtual double integ(double a, double b) const {
276 
277  size_t cache=0;
278  size_t i, index_a, index_b;
279 
280  bool flip=false;
281  if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
282  ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
283  double tmp=a;
284  a=b;
285  b=tmp;
286  flip=true;
287  }
288 
289  index_a=this->svx.find_const(a,cache);
290  index_b=this->svx.find_const(b,cache);
291 
292  double result=0.0;
293  for(i=index_a; i<=index_b; i++) {
294 
295  double x_lo=(*this->px)[i];
296  double x_hi=(*this->px)[i+1];
297  double y_lo=(*this->py)[i];
298  double y_hi=(*this->py)[i+1];
299  double dx=x_hi-x_lo;
300 
301  if(dx != 0.0) {
302 
303  if (i == index_a || i == index_b) {
304  double x1=(i == index_a) ? a : x_lo;
305  double x2=(i == index_b) ? b : x_hi;
306  double D=(y_hi-y_lo)/dx;
307  result += (x2-x1)*(y_lo+0.5*D*((x2-x_lo)+(x1-x_lo)));
308  } else {
309  result += 0.5*dx*(y_lo+y_hi);
310  }
311 
312  } else {
313  std::string str=((std::string)"Interval of length zero ")+
314  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
315  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
316  " in interp_linear::integ().";
317  O2SCL_ERR(str.c_str(),exc_einval);
318  }
319  }
320 
321  if (flip) result=-result;
322  return result;
323  }
324 
325  /// Return the type, \c "interp_linear".
326  virtual const char *type() const { return "interp_linear"; }
327 
328 #ifndef DOXYGEN_INTERNAL
329 
330  private:
331 
334 
335 #endif
336 
337  };
338 
339  /** \brief Cubic spline interpolation (GSL)
340 
341  See also the \ref intp_section section of the \o2 User's guide.
342 
343  By default, this class uses natural boundary conditions, where
344  the second derivative vanishes at each end point. Extrapolation
345  effectively assumes that the second derivative is linear outside
346  of the endpoints.
347  */
348  template<class vec_t, class vec2_t=vec_t> class interp_cspline :
349  public interp_base<vec_t,vec2_t> {
350 
351 #ifdef O2SCL_NEVER_DEFINED
352  }{
353 #endif
354 
355  public:
356 
358  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
359  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
360  typedef boost::numeric::ublas::slice slice;
361  typedef boost::numeric::ublas::range range;
362 
363 #ifndef DOXYGEN_INTERNAL
364 
365  protected:
366 
367  /// \name Storage for cubic spline interpolation
368  //@{
369  ubvector c;
370  ubvector g;
371  ubvector diag;
372  ubvector offdiag;
373  //@}
374 
375  /// Memory for the tridiagonalization
377 
378  /// Compute coefficients for cubic spline interpolation
379  void coeff_calc(const ubvector &c_array, double dy, double dx,
380  size_t index, double &b, double &c2, double &d) const {
381 
382  double c_i=c_array[index];
383  double c_ip1=c_array[index+1];
384  b=(dy/dx)-dx*(c_ip1+2.0*c_i)/3.0;
385  c2=c_i;
386  d=(c_ip1-c_i)/(3.0*dx);
387 
388  return;
389  }
390 
391 #endif
392 
393  public:
394 
395  /** \brief Create a base interpolation object with natural or
396  periodic boundary conditions
397  */
399  this->min_size=3;
400  }
401 
402  virtual ~interp_cspline() {
403  }
404 
405  /** \brief Initialize interpolation routine
406  */
407  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
408 
409  if (size<this->min_size) {
410  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
411  " than "+szttos(this->min_size)+" in interp_cspline::"+
412  "set().").c_str(),exc_einval);
413  }
414 
415  if (size!=this->sz) {
416  c.resize(size);
417  g.resize(size);
418  diag.resize(size);
419  offdiag.resize(size);
420  p4m.resize(size);
421  }
422 
423  this->px=&xa;
424  this->py=&ya;
425  this->sz=size;
426 
427  this->svx.set_vec(size,xa);
428 
429  /// Natural boundary conditions
430 
431  size_t i;
432  size_t num_points=size;
433  size_t max_index=num_points-1;
434  size_t sys_size=max_index-1;
435 
436  c[0]=0.0;
437  c[max_index]=0.0;
438 
439  for (i=0; i < sys_size; i++) {
440  double h_i=xa[i+1]-xa[i];
441  double h_ip1=xa[i+2]-xa[i+1];
442  double ydiff_i=ya[i+1]-ya[i];
443  double ydiff_ip1=ya[i+2]-ya[i+1];
444  double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
445  double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
446  offdiag[i]=h_ip1;
447  diag[i]=2.0*(h_ip1+h_i);
448  g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
449  }
450 
451  if (sys_size == 1) {
452 
453  c[1]=g[0]/diag[0];
454 
455  return;
456  }
457 
458  ubvector_range cp1(c,range(1,c.size()));
459  o2scl_linalg::solve_tridiag_sym<ubvector,ubvector,ubvector,
460  ubvector_range,o2scl_linalg::ubvector_4_mem,ubvector>
461  (diag,offdiag,g,cp1,sys_size,p4m);
462 
463  return;
464  }
465 
466  /// Give the value of the function \f$ y(x=x_0) \f$ .
467  virtual double eval(double x0) const {
468 
469  size_t cache=0;
470  size_t index=this->svx.find_const(x0,cache);
471 
472  double x_lo=(*this->px)[index];
473  double x_hi=(*this->px)[index+1];
474  double dx=x_hi-x_lo;
475 
476  double y_lo=(*this->py)[index];
477  double y_hi=(*this->py)[index+1];
478  double dy=y_hi-y_lo;
479  double delx=x0-x_lo;
480  double b_i, c_i, d_i;
481 
482  coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
483 
484  return y_lo+delx*(b_i+delx*(c_i+delx*d_i));
485  }
486 
487  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
488  virtual double deriv(double x0) const {
489 
490  size_t cache=0;
491  size_t index=this->svx.find_const(x0,cache);
492 
493  double x_lo=(*this->px)[index];
494  double x_hi=(*this->px)[index+1];
495  double dx=x_hi-x_lo;
496 
497  double y_lo=(*this->py)[index];
498  double y_hi=(*this->py)[index+1];
499  double dy=y_hi-y_lo;
500  double delx=x0-x_lo;
501  double b_i, c_i, d_i;
502 
503  coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
504 
505  return b_i+delx*(2.0*c_i+3.0*d_i*delx);
506  }
507 
508  /** \brief Give the value of the second derivative
509  \f$ y^{\prime \prime}(x=x_0) \f$ .
510  */
511  virtual double deriv2(double x0) const {
512 
513  size_t cache=0;
514  size_t index=this->svx.find_const(x0,cache);
515 
516  double x_lo=(*this->px)[index];
517  double x_hi=(*this->px)[index+1];
518  double dx=x_hi-x_lo;
519 
520  double y_lo=(*this->py)[index];
521  double y_hi=(*this->py)[index+1];
522  double dy=y_hi-y_lo;
523  double delx=x0-x_lo;
524  double b_i, c_i, d_i;
525 
526  coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
527 
528  return 2.0*c_i+6.0*d_i*delx;
529  }
530 
531  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
532  virtual double integ(double a, double b) const {
533 
534  size_t i, index_a, index_b;
535 
536  bool flip=false;
537  if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
538  ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
539  double tmp=a;
540  a=b;
541  b=tmp;
542  flip=true;
543  }
544 
545  size_t cache=0;
546  index_a=this->svx.find_const(a,cache);
547  index_b=this->svx.find_const(b,cache);
548 
549  double result=0.0;
550 
551  for(i=index_a; i<=index_b; i++) {
552 
553  double x_lo=(*this->px)[i];
554  double x_hi=(*this->px)[i+1];
555  double y_lo=(*this->py)[i];
556  double y_hi=(*this->py)[i+1];
557  double dx=x_hi-x_lo;
558  double dy=y_hi-y_lo;
559 
560  if(dx != 0.0) {
561  double b_i, c_i, d_i;
562  coeff_calc(c,dy,dx,i,b_i,c_i,d_i);
563  if (i == index_a || i == index_b) {
564  double x1=(i == index_a) ? a : x_lo;
565  double x2=(i == index_b) ? b : x_hi;
566  result += this->integ_eval(y_lo,b_i,c_i,d_i,x_lo,x1,x2);
567  } else {
568  result += dx*(y_lo+dx*(0.5*b_i+
569  dx*(c_i/3.0+0.25*d_i*dx)));
570  }
571  } else {
572  std::string str=((std::string)"Interval of length zero ")+
573  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
574  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
575  " in interp_cspline::integ().";
576  O2SCL_ERR(str.c_str(),exc_einval);
577  }
578 
579  }
580 
581  if (flip) result*=-1.0;
582 
583  return result;
584  }
585 
586  /// Return the type, \c "interp_cspline".
587  virtual const char *type() const { return "interp_cspline"; }
588 
589 #ifndef DOXYGEN_INTERNAL
590 
591  private:
592 
595  (const interp_cspline<vec_t,vec2_t>&);
596 
597 #endif
598 
599  };
600 
601  /** \brief Cubic spline interpolation with periodic
602  boundary conditions (GSL)
603 
604  See also the \ref intp_section section of the \o2 User's guide.
605  */
606  template<class vec_t, class vec2_t=vec_t> class interp_cspline_peri :
607  public interp_cspline<vec_t,vec2_t> {
608 
609 #ifdef O2SCL_NEVER_DEFINED
610  }{
611 #endif
612 
613  protected:
614 
615  /// Memory for the tridiagonalization
617 
618  public:
619 
621  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
622  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
623  typedef boost::numeric::ublas::slice slice;
624  typedef boost::numeric::ublas::range range;
625 
626  interp_cspline_peri() : interp_cspline<vec_t,vec2_t>() {
627  this->min_size=2;
628  }
629 
630  virtual ~interp_cspline_peri() {
631  }
632 
633  /// Return the type, \c "interp_cspline_peri".
634  virtual const char *type() const { return "interp_cspline_peri"; }
635 
636  /** \brief Initialize interpolation routine
637  */
638  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
639 
640  if (size<this->min_size) {
641  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
642  " than "+szttos(this->min_size)+" in interp_cspline"+
643  "_peri::set().").c_str(),exc_einval);
644  }
645 
646  if (size!=this->sz) {
647  this->c.resize(size);
648  this->g.resize(size);
649  this->diag.resize(size);
650  this->offdiag.resize(size);
651  p5m.resize(size);
652  }
653 
654  this->px=&xa;
655  this->py=&ya;
656  this->sz=size;
657 
658  this->svx.set_vec(size,xa);
659 
660  /// Periodic boundary conditions
661 
662  size_t i;
663  size_t num_points=size;
664  // Engeln-Mullges+Uhlig "n"
665  size_t max_index=num_points-1;
666  // linear system is sys_size x sys_size
667  size_t sys_size=max_index;
668 
669  if (sys_size == 2) {
670 
671  // solve 2x2 system
672 
673  double h0=xa[1]-xa[0];
674  double h1=xa[2]-xa[1];
675  double h2=xa[3]-xa[2];
676  double A=2.0*(h0+h1);
677  double B=h0+h1;
678  double gx[2];
679  double det;
680 
681  gx[0]=3.0*((ya[2]-ya[1])/h1-(ya[1]-ya[0])/h0);
682  gx[1]=3.0*((ya[1]-ya[2])/h2-(ya[2]-ya[1])/h1);
683 
684  det=3.0*(h0+h1)*(h0+h1);
685  this->c[1]=( A*gx[0]-B*gx[1])/det;
686  this->c[2]=(-B*gx[0]+A*gx[1])/det;
687  this->c[0]=this->c[2];
688 
689  return;
690 
691  } else {
692 
693  for (i=0; i < sys_size-1; i++) {
694  double h_i=xa[i+1]-xa[i];
695  double h_ip1=xa[i+2]-xa[i+1];
696  double ydiff_i=ya[i+1]-ya[i];
697  double ydiff_ip1=ya[i+2]-ya[i+1];
698  double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
699  double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
700  this->offdiag[i]=h_ip1;
701  this->diag[i]=2.0*(h_ip1+h_i);
702  this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
703  }
704 
705  i=sys_size-1;
706  {
707  double h_i=xa[i+1]-xa[i];
708  double h_ip1=xa[1]-xa[0];
709  double ydiff_i=ya[i+1]-ya[i];
710  double ydiff_ip1=ya[1]-ya[0];
711  double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
712  double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
713  this->offdiag[i]=h_ip1;
714  this->diag[i]=2.0*(h_ip1+h_i);
715  this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
716  }
717 
718  ubvector_range cp1(this->c,range(1,this->c.size()));
719  o2scl_linalg::solve_cyc_tridiag_sym<ubvector,ubvector,ubvector,
720  ubvector_range,o2scl_linalg::ubvector_5_mem,ubvector>
721  (this->diag,this->offdiag,this->g,cp1,sys_size,p5m);
722  this->c[0]=this->c[max_index];
723 
724  return;
725  }
726 
727  }
728 
729 #ifndef DOXYGEN_INTERNAL
730 
731  private:
732 
736  (const interp_cspline_peri<vec_t,vec2_t>&);
737 
738 #endif
739 
740  };
741 
742  /** \brief Akima spline interpolation (GSL)
743 
744  See also the \ref intp_section section of the \o2 User's guide.
745 
746  This class uses natural boundary conditions, where the second
747  derivative vanishes at each end point. Extrapolation effectively
748  assumes that the second derivative is linear outside of the
749  endpoints. Use \ref interp_akima_peri for perodic boundary
750  conditions.
751  */
752  template<class vec_t, class vec2_t=vec_t> class interp_akima :
753  public interp_base<vec_t,vec2_t> {
754 
755  public:
756 
758  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
759  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
760  typedef boost::numeric::ublas::slice slice;
761  typedef boost::numeric::ublas::range range;
762 
763 #ifndef DOXYGEN_INTERNAL
764 
765  protected:
766 
767  /// \name Storage for Akima spline interpolation
768  //@{
769  ubvector b;
770  ubvector c;
771  ubvector d;
772  ubvector um;
773  //@}
774 
775  /// For initializing the interpolation
776  void akima_calc(const vec_t &x_array, size_t size,
777  ubvector &umx) {
778 
779  for(size_t i=0;i<this->sz-1;i++) {
780 
781  double NE=fabs(umx[3+i]-umx[2+i])+fabs(umx[1+i]-umx[i]);
782 
783  if (NE == 0.0) {
784  b[i]=umx[2+i];
785  c[i]=0.0;
786  d[i]=0.0;
787  } else {
788  double h_i=(*this->px)[i+1]-(*this->px)[i];
789  double NE_next=fabs(umx[4+i]-umx[3+i])+
790  fabs(umx[2+i]-umx[1+i]);
791  double alpha_i=fabs(umx[1+i]-umx[i])/NE;
792  double alpha_ip1;
793  double tL_ip1;
794  if (NE_next == 0.0) {
795  tL_ip1=umx[2+i];
796  } else {
797  alpha_ip1=fabs(umx[2+i]-umx[1+i])/NE_next;
798  tL_ip1=(1.0-alpha_ip1)*umx[2+i]+alpha_ip1*umx[3+i];
799  }
800  b[i]=(1.0-alpha_i)*umx[1+i]+alpha_i*umx[2+i];
801  c[i]=(3.0*umx[2+i]-2.0*b[i]-tL_ip1)/h_i;
802  d[i]=(b[i]+tL_ip1-2.0*umx[2+i])/(h_i*h_i);
803  }
804  }
805  }
806 
807 #endif
808 
809  public:
810 
811  /** \brief Create a base interpolation object with or without
812  periodic boundary conditions
813  */
815  this->min_size=5;
816  }
817 
818  virtual ~interp_akima() {
819  }
820 
821  /** \brief Initialize interpolation routine
822  */
823  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
824 
825  if (size<this->min_size) {
826  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
827  " than "+szttos(this->min_size)+" in interp_akima::"+
828  "set().").c_str(),exc_einval);
829  }
830 
831  if (size!=this->sz) {
832  b.resize(size);
833  c.resize(size);
834  d.resize(size);
835  um.resize(size+4);
836  }
837 
838  this->px=&xa;
839  this->py=&ya;
840  this->sz=size;
841 
842  this->svx.set_vec(size,xa);
843 
844  // Non-periodic boundary conditions
845 
846  ubvector_range m(um,range(2,um.size()));
847  size_t i;
848  for (i=0;i<=size-2;i++) {
849  m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
850  }
851 
852  um[0]=3.0*m[0]-2.0*m[1];
853  um[1]=2.0*m[0]-m[1];
854  m[this->sz-1]=2.0*m[size-2]-m[size-3];
855  m[size]=3.0*m[size-2]-2.0*m[size-3];
856 
857  akima_calc(xa,size,um);
858 
859  return;
860  }
861 
862  /// Give the value of the function \f$ y(x=x_0) \f$ .
863  virtual double eval(double x0) const {
864 
865  size_t cache=0;
866  size_t index=this->svx.find_const(x0,cache);
867 
868  double x_lo=(*this->px)[index];
869  double delx=x0-x_lo;
870  double bb=b[index];
871  double cc=c[index];
872  double dd=d[index];
873 
874  return (*this->py)[index]+delx*(bb+delx*(cc+dd*delx));
875  }
876 
877  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
878  virtual double deriv(double x0) const {
879 
880  size_t cache=0;
881  size_t index=this->svx.find_const(x0,cache);
882 
883  double x_lo=(*this->px)[index];
884  double delx=x0-x_lo;
885  double bb=b[index];
886  double cc=c[index];
887  double dd=d[index];
888 
889  return bb+delx*(2.0*cc+3.0*dd*delx);
890  }
891 
892  /** \brief Give the value of the second derivative
893  \f$ y^{\prime \prime}(x=x_0) \f$ .
894  */
895  virtual double deriv2(double x0) const {
896 
897  size_t cache=0;
898  size_t index=this->svx.find_const(x0,cache);
899 
900  double x_lo=(*this->px)[index];
901  double delx=x0-x_lo;
902  double cc=c[index];
903  double dd=d[index];
904 
905  return 2.0*cc+6.0*dd*delx;
906  }
907 
908  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
909  virtual double integ(double aa, double bb) const {
910 
911  size_t i, index_a, index_b;
912 
913  bool flip=false;
914  if (((*this->px)[0]<(*this->px)[this->sz-1] && aa>bb) ||
915  ((*this->px)[0]>(*this->px)[this->sz-1] && aa<bb)) {
916  double tmp=aa;
917  aa=bb;
918  bb=tmp;
919  flip=true;
920  }
921 
922  size_t cache=0;
923  index_a=this->svx.find_const(aa,cache);
924  index_b=this->svx.find_const(bb,cache);
925 
926  double result=0.0;
927 
928  for(i=index_a; i<=index_b; i++) {
929 
930  double x_lo=(*this->px)[i];
931  double x_hi=(*this->px)[i+1];
932  double y_lo=(*this->py)[i];
933  double dx=x_hi-x_lo;
934 
935  if (dx != 0.0) {
936 
937  if (i==index_a || i==index_b) {
938  double x1=(i==index_a) ? aa : x_lo;
939  double x2=(i==index_b) ? bb : x_hi;
940  result += this->integ_eval(y_lo,b[i],c[i],d[i],x_lo,x1,x2);
941  } else {
942  result+=dx*(y_lo+dx*(0.5*b[i]+dx*(c[i]/3.0+0.25*d[i]*dx)));
943  }
944 
945  } else {
946  double y_hi=(*this->py)[i+1];
947  std::string str=((std::string)"Interval of length zero ")+
948  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
949  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
950  " in interp_akima::integ().";
951  O2SCL_ERR(str.c_str(),exc_einval);
952  }
953  }
954 
955  if (flip) result*=-1.0;
956 
957  return result;
958  }
959 
960  /// Return the type, \c "interp_akima".
961  virtual const char *type() const { return "interp_akima"; }
962 
963 #ifndef DOXYGEN_INTERNAL
964 
965  private:
966 
969 
970 #endif
971 
972  };
973 
974  /** \brief Akima spline interpolation with periodic
975  boundary conditions (GSL)
976 
977  See also the \ref intp_section section of the \o2 User's guide.
978  */
979  template<class vec_t, class vec2_t=vec_t> class interp_akima_peri :
980  public interp_akima<vec_t,vec2_t> {
981 
982  public:
983 
985  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
986  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
987  typedef boost::numeric::ublas::slice slice;
988  typedef boost::numeric::ublas::range range;
989 
990  public:
991 
993  }
994 
995  virtual ~interp_akima_peri() {
996  }
997 
998  /// Return the type, \c "interp_akima_peri".
999  virtual const char *type() const { return "interp_akima_peri"; }
1000 
1001  /// Initialize interpolation routine
1002  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
1003 
1004  if (size<this->min_size) {
1005  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
1006  " than "+szttos(this->min_size)+" in interp_akima"+
1007  "_peri::set().").c_str(),exc_einval);
1008  }
1009 
1010  if (size!=this->sz) {
1011  this->b.resize(size);
1012  this->c.resize(size);
1013  this->d.resize(size);
1014  this->um.resize(size+4);
1015  }
1016 
1017  this->px=&xa;
1018  this->py=&ya;
1019  this->sz=size;
1020 
1021  this->svx.set_vec(size,xa);
1022 
1023  // Periodic boundary conditions
1024 
1025  ubvector_range m(this->um,range(2,this->um.size()));
1026 
1027  // Form the required set of divided differences
1028  for (size_t i=0;i<=size-2;i++) {
1029  m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
1030  }
1031 
1032  this->um[0]=m[this->sz-3];
1033  this->um[1]=m[this->sz-2];
1034  m[this->sz-1]=m[0];
1035  m[this->sz]=m[1];
1036 
1037  this->akima_calc(xa,size,this->um);
1038 
1039  return;
1040  }
1041 
1042 #ifndef DOXYGEN_INTERNAL
1043 
1044  private:
1045 
1048  (const interp_akima_peri<vec_t,vec2_t>&);
1049 
1050 #endif
1051 
1052  };
1053 
1054  /** \brief Steffen's monotonicity-preserving interpolation
1055 
1056  Adapted from the GSL version by J.-F. Caron which
1057  was based on \ref Steffen90 .
1058  */
1059  template<class vec_t, class vec2_t=vec_t> class interp_steffen :
1060  public interp_base<vec_t,vec2_t> {
1061 
1062 #ifdef O2SCL_NEVER_DEFINED
1063  }{
1064 #endif
1065 
1066  public:
1067 
1069  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
1070  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
1071  typedef boost::numeric::ublas::slice slice;
1072  typedef boost::numeric::ublas::range range;
1073 
1074 #ifndef DOXYGEN_INTERNAL
1075 
1076  protected:
1077 
1078  /// \name Storage for cubic spline interpolation
1079  //@{
1080  ubvector a;
1081  ubvector b;
1082  ubvector c;
1083  ubvector d;
1084  ubvector y_prime;
1085  //@}
1086 
1087  /** \brief Flip the sign of x if x and y have different signs
1088  */
1089  double copysign(const double x, const double y) {
1090  if ((x < 0 && y > 0) || (x > 0 && y < 0)) {
1091  return -x;
1092  }
1093  return x;
1094  }
1095 
1096 #endif
1097 
1098  public:
1099 
1100  /** \brief Create a base interpolation object */
1102  this->min_size=3;
1103  }
1104 
1105  virtual ~interp_steffen() {
1106  }
1107 
1108  /** \brief Initialize interpolation routine
1109  */
1110  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
1111 
1112  if (size<this->min_size) {
1113  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
1114  " than "+szttos(this->min_size)+" in interp_steffen::"+
1115  "set().").c_str(),exc_einval);
1116  }
1117 
1118  if (size!=this->sz) {
1119  a.resize(size);
1120  b.resize(size);
1121  c.resize(size);
1122  d.resize(size);
1123  y_prime.resize(size);
1124  }
1125 
1126  this->px=&xa;
1127  this->py=&ya;
1128  this->sz=size;
1129 
1130  this->svx.set_vec(size,xa);
1131 
1132  /*
1133  * first assign the interval and slopes for the left boundary.
1134  * We use the "simplest possibility" method described in the paper
1135  * in section 2.2
1136  */
1137  double h0=(xa[1]-xa[0]);
1138  double s0=(ya[1]-ya[0]) / h0;
1139 
1140  y_prime[0]=s0;
1141 
1142  /* Now we calculate all the necessary s, h, p, and y' variables
1143  from 1 to N-2 (0 to size-2 inclusive) */
1144  for (size_t i=1; i < (size-1); i++) {
1145 
1146  double pi;
1147 
1148  /* equation 6 in the paper */
1149  double hi=(xa[i+1]-xa[i]);
1150  double him1=(xa[i]-xa[i-1]);
1151 
1152  /* equation 7 in the paper */
1153  double si=(ya[i+1]-ya[i]) / hi;
1154  double sim1=(ya[i]-ya[i-1]) / him1;
1155 
1156  /* equation 8 in the paper */
1157  pi=(sim1*hi + si*him1) / (him1 + hi);
1158 
1159  /* This is a C equivalent of the FORTRAN statement below eqn 11 */
1160  y_prime[i]=(copysign(1.0,sim1)+copysign(1.0,si))*
1161  std::min(fabs(sim1),std::min(fabs(si),0.5*fabs(pi)));
1162 
1163  }
1164 
1165  /*
1166  * we also need y' for the rightmost boundary; we use the
1167  * "simplest possibility" method described in the paper in
1168  * section 2.2
1169  *
1170  * y'=s_{n-1}
1171  */
1172  y_prime[size-1]=(ya[size-1]-ya[size-2])/
1173  (xa[size-1]-xa[size-2]);
1174 
1175  /* Now we can calculate all the coefficients for the whole range. */
1176  for (size_t i=0; i < (size-1); i++) {
1177  double hi=(xa[i+1]-xa[i]);
1178  double si=(ya[i+1]-ya[i]) / hi;
1179 
1180  /* These are from equations 2-5 in the paper. */
1181  a[i]=(y_prime[i] + y_prime[i+1]-2*si) / hi / hi;
1182  b[i]=(3*si-2*y_prime[i]-y_prime[i+1]) / hi;
1183  c[i]=y_prime[i];
1184  d[i]=ya[i];
1185  }
1186 
1187  return;
1188  }
1189 
1190  /// Give the value of the function \f$ y(x=x_0) \f$ .
1191  virtual double eval(double x0) const {
1192 
1193  size_t cache=0;
1194  size_t index=this->svx.find_const(x0,cache);
1195  double x_lo=(*this->px)[index];
1196  double delx=x0-x_lo;
1197 
1198  /* Use Horner's scheme for efficient evaluation of polynomials. */
1199  double y = d[index]+delx*(c[index]+delx*(b[index]+delx*a[index]));
1200 
1201  return y;
1202  }
1203 
1204  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1205  virtual double deriv(double x0) const {
1206 
1207  size_t cache=0;
1208  size_t index=this->svx.find_const(x0,cache);
1209  double x_lo=(*this->px)[index];
1210  double delx=x0-x_lo;
1211 
1212  return c[index]+delx*(2.0*b[index]+delx*3.0*a[index]);
1213  }
1214 
1215  /** \brief Give the value of the second derivative
1216  \f$ y^{\prime \prime}(x=x_0) \f$ .
1217  */
1218  virtual double deriv2(double x0) const {
1219 
1220  size_t cache=0;
1221  size_t index=this->svx.find_const(x0,cache);
1222  double x_lo=(*this->px)[index];
1223  double delx=x0-x_lo;
1224 
1225  return 2.0*b[index]+delx*6.0*a[index];
1226  }
1227 
1228  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1229  virtual double integ(double al, double bl) const {
1230 
1231  size_t i, index_a, index_b;
1232 
1233  bool flip=false;
1234  if (((*this->px)[0]<(*this->px)[this->sz-1] && al>bl) ||
1235  ((*this->px)[0]>(*this->px)[this->sz-1] && al<bl)) {
1236  double tmp=al;
1237  al=bl;
1238  bl=tmp;
1239  flip=true;
1240  }
1241 
1242  size_t cache=0;
1243  index_a=this->svx.find_const(al,cache);
1244  index_b=this->svx.find_const(bl,cache);
1245 
1246  double result=0.0;
1247 
1248  for(i=index_a; i<=index_b; i++) {
1249 
1250  double x_lo=(*this->px)[i];
1251  double x_hi=(*this->px)[i+1];
1252  double y_lo=(*this->py)[i];
1253  double y_hi=(*this->py)[i+1];
1254  double dx=x_hi-x_lo;
1255 
1256  if(dx != 0.0) {
1257 
1258  double x1=(i == index_a) ? al-x_lo : 0.0;
1259  double x2=(i == index_b) ? bl-x_lo : x_hi-x_lo;
1260  result += (1.0/4.0)*a[i]*(x2*x2*x2*x2-x1*x1*x1*x1)+
1261  (1.0/3.0)*b[i]*(x2*x2*x2-x1*x1*x1)+
1262  (1.0/2.0)*c[i]*(x2*x2-x1*x1)+d[i]*(x2-x1);
1263 
1264  } else {
1265  std::string str=((std::string)"Interval of length zero ")+
1266  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
1267  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
1268  " in interp_steffen::integ().";
1269  O2SCL_ERR(str.c_str(),exc_einval);
1270  }
1271 
1272  }
1273 
1274  if (flip) result*=-1.0;
1275 
1276  return result;
1277  }
1278 
1279  /// Return the type, \c "interp_steffen".
1280  virtual const char *type() const { return "interp_steffen"; }
1281 
1282 #ifndef DOXYGEN_INTERNAL
1283 
1284  private:
1285 
1287  interp_steffen<vec_t,vec2_t>& operator=
1288  (const interp_steffen<vec_t,vec2_t>&);
1289 
1290 #endif
1291 
1292  };
1293 
1294  /** \brief Monotonicity-preserving interpolation
1295 
1296  \warning This class is experimental. Integrals don't work yet.
1297 
1298  This class uses a method based on cubic Hermite interpolation,
1299  modifying the slopes to guarantee monotonicity. In the
1300  notation of \ref Fritsch80, if
1301  \f[
1302  \alpha_i^2+\beta_i^2 \geq 9
1303  \f]
1304  then \f$ \alpha \f$ and \f$ \beta \f$ are decreased by
1305  the factor \f$ \tau \f$ as described at the end of
1306  section 4.
1307 
1308  \note The results of the interpolation will only be monotonic in
1309  the regions where the original data set is also monotonic. Also,
1310  this interpolation routine does not enforce strict monotonicity,
1311  and the results of the interpolation will be flat where the data
1312  is also flat.
1313 
1314  Based on \ref Fritsch80 .
1315 
1316  \future Convert into fewer loops over the data
1317  */
1318  template<class vec_t, class vec2_t=vec_t> class interp_monotonic :
1319  public interp_base<vec_t,vec2_t> {
1320 
1321 #ifdef O2SCL_NEVER_DEFINED
1322  }{
1323 #endif
1324 
1325  public:
1326 
1328 
1329  protected:
1330 
1331  /// Slopes
1332  ubvector m;
1333  /// Finite differences
1334  ubvector Delta;
1335  /// Ratio
1336  ubvector alpha;
1337  /// Staggered ratio
1338  ubvector beta;
1339 
1340  public:
1341 
1342  interp_monotonic() {
1343  this->min_size=2;
1344  }
1345 
1346  virtual ~interp_monotonic() {
1347  }
1348 
1349  /// Initialize interpolation routine
1350  virtual void set(size_t size, const vec_t &x, const vec2_t &y) {
1351 
1352  // Verify size
1353  if (size<this->min_size) {
1354  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
1355  " than "+szttos(this->min_size)+" in interp_monotonic::"+
1356  "set().").c_str(),exc_einval);
1357  }
1358 
1359  // Setup search_vec object
1360  this->svx.set_vec(size,x);
1361 
1362  // Resize internal vectors
1363  if (this->sz!=size) {
1364  m.resize(size);
1365  Delta.resize(size-1);
1366  alpha.resize(size-1);
1367  beta.resize(size-1);
1368  }
1369 
1370  // Copy pointers
1371  this->px=&x;
1372  this->py=&y;
1373  this->sz=size;
1374 
1375  // Create the interpolation arrays in this sub-interval
1376 
1377  // Compute Delta and m
1378  for(size_t i=0;i<size-1;i++) {
1379  Delta[i]=(y[i+1]-y[i])/(x[i+1]-x[i]);
1380  if (i>0) {
1381  m[i]=(Delta[i]+Delta[i-1])/2.0;
1382  }
1383  }
1384  m[0]=Delta[0];
1385  m[size-1]=Delta[size-2];
1386 
1387  // Check to see if the data is flat anywhere
1388  for(size_t i=0;i<size-1;i++) {
1389  if (y[i]==y[i+1]) {
1390  m[i]=0.0;
1391  m[i+1]=0.0;
1392  }
1393  }
1394 
1395  // Compute alpha and beta
1396  for(size_t i=0;i<size-1;i++) {
1397  alpha[i]=m[i]/Delta[i];
1398  beta[i]=m[i+1]/Delta[i];
1399  }
1400 
1401  // Constrain m to ensure monotonicity
1402  for(size_t i=0;i<size-1;i++) {
1403  double norm2=alpha[i]*alpha[i]+beta[i]*beta[i];
1404  if (norm2>9.0) {
1405  double tau=3.0/sqrt(norm2);
1406  m[i]=tau*alpha[i]*Delta[i];
1407  m[i+1]=tau*beta[i]*Delta[i];
1408  }
1409  }
1410 
1411  return;
1412  }
1413 
1414  /// Give the value of the function \f$ y(x=x_0) \f$ .
1415  virtual double eval(double x0) const {
1416 
1417  size_t cache=0;
1418  size_t index=this->svx.find_const(x0,cache);
1419 
1420  double x_lo=(*this->px)[index];
1421  double x_hi=(*this->px)[index+1];
1422  double y_lo=(*this->py)[index];
1423  double y_hi=(*this->py)[index+1];
1424  double h=x_hi-x_lo;
1425  double t=(x0-x_lo)/h;
1426  double t2=t*t, t3=t2*t;
1427 
1428  double h00=2.0*t3-3.0*t2+1.0;
1429  double h10=t3-2.0*t2+t;
1430  double h01=-2.0*t3+3.0*t2;
1431  double h11=t3-t2;
1432  double interp=y_lo*h00+h*m[index]*h10+y_hi*h01+h*m[index+1]*h11;
1433 
1434  return interp;
1435  }
1436 
1437  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1438  virtual double deriv(double x0) const {
1439 
1440  size_t cache=0;
1441  size_t index=this->svx.find_const(x0,cache);
1442 
1443  double x_lo=(*this->px)[index];
1444  double x_hi=(*this->px)[index+1];
1445  double y_lo=(*this->py)[index];
1446  double y_hi=(*this->py)[index+1];
1447  double h=x_hi-x_lo;
1448  double t=(x0-x_lo)/h;
1449  double t2=t*t;
1450 
1451  double dh00=6.0*t2-6.0*t;
1452  double dh10=3.0*t2-4.0*t+1.0;
1453  double dh01=-6.0*t2+6.0*t;
1454  double dh11=3.0*t2-2.0*t;
1455  double deriv=(y_lo*dh00+h*m[index]*dh10+y_hi*dh01+
1456  h*m[index+1]*dh11)/h;
1457 
1458  return deriv;
1459  }
1460 
1461  /** \brief Give the value of the second derivative
1462  \f$ y^{\prime \prime}(x=x_0) \f$
1463  */
1464  virtual double deriv2(double x0) const {
1465 
1466  size_t cache=0;
1467  size_t index=this->svx.find_const(x0,cache);
1468 
1469  double x_lo=(*this->px)[index];
1470  double x_hi=(*this->px)[index+1];
1471  double y_lo=(*this->py)[index];
1472  double y_hi=(*this->py)[index+1];
1473  double h=x_hi-x_lo;
1474  double t=(x0-x_lo)/h;
1475 
1476  double ddh00=12.0*t-6.0;
1477  double ddh10=6.0*t-4.0;
1478  double ddh01=-12.0*t+6.0;
1479  double ddh11=6.0*t-2.0;
1480  double deriv2=(y_lo*ddh00+h*m[index]*ddh10+y_hi*ddh01+
1481  h*m[index+1]*ddh11)/h/h;
1482 
1483  return deriv2;
1484  }
1485 
1486  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1487  virtual double integ(double a, double b) const {
1488 
1489  size_t i, index_a, index_b;
1490 
1491  bool flip=false;
1492  if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
1493  ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
1494  double tmp=a;
1495  a=b;
1496  b=tmp;
1497  flip=true;
1498  }
1499 
1500  size_t cache=0;
1501  index_a=this->svx.find_const(a,cache);
1502  index_b=this->svx.find_const(b,cache);
1503 
1504  double result=0.0;
1505 
1506  for(i=index_a; i<=index_b; i++) {
1507 
1508  double x_hi=(*this->px)[i+1];
1509  double x_lo=(*this->px)[i];
1510  double y_lo=(*this->py)[i];
1511  double y_hi=(*this->py)[i+1];
1512  double h=x_hi-x_lo;
1513 
1514  if (h != 0.0) {
1515 
1516  if (i == index_a) {
1517  x_lo=a;
1518  }
1519  if (i == index_b) {
1520  x_hi=b;
1521  }
1522 
1523  double t=(x_hi-x_lo)/h;
1524  double t2=t*t, t3=t2*t, t4=t3*t;
1525 
1526  double ih00=t4/2.0-t3+t;
1527  double ih10=t4/4.0-2.0*t3/3.0+t2/2.0;
1528  double ih01=-t4/2.0+t3;
1529  double ih11=t4/4.0-t3/3.0;
1530  double intres=h*(y_lo*ih00+h*m[i]*ih10+y_hi*ih01+
1531  h*m[i+1]*ih11);
1532  result+=intres;
1533 
1534  } else {
1535  std::string str=((std::string)"Interval of length zero ")+
1536  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
1537  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
1538  " in interp_monotonic::integ().";
1539  O2SCL_ERR(str.c_str(),exc_einval);
1540  }
1541 
1542  }
1543 
1544  if (flip) result*=-1.0;
1545 
1546  return result;
1547  }
1548 
1549  /// Return the type, \c "interp_monotonic".
1550  virtual const char *type() const { return "interp_monotonic"; }
1551 
1552 #ifndef DOXYGEN_INTERNAL
1553 
1554  private:
1555 
1558  (const interp_monotonic<vec_t,vec2_t>&);
1559 
1560 #endif
1561 
1562  };
1563 
1564  /** \brief Interpolation class for general vectors
1565 
1566  See also the \ref intp_section section of the \o2 User's guide.
1567 
1568  Interpolation of ublas vector like objects is performed with the
1569  default template parameters, and \ref interp_array is provided
1570  for simple interpolation on C-style arrays.
1571 
1572  The type of interpolation to be performed can be specified using
1573  the set_type() function or in the constructor. The default is
1574  cubic splines with natural boundary conditions.
1575  */
1576  template<class vec_t=boost::numeric::ublas::vector<double>,
1577  class vec2_t=vec_t> class interp {
1578 
1579 #ifndef DOXYGEN_INTERNAL
1580 
1581  protected:
1582 
1583  /// Pointer to base interpolation object
1585 
1586 #endif
1587 
1588  public:
1589 
1590  /// Create with base interpolation object \c it
1591  interp(size_t interp_type=itp_cspline) {
1592  if (interp_type==itp_linear) {
1594  } else if (interp_type==itp_cspline) {
1596  } else if (interp_type==itp_cspline_peri) {
1598  } else if (interp_type==itp_akima) {
1600  } else if (interp_type==itp_akima_peri) {
1602  } else if (interp_type==itp_monotonic) {
1604  } else if (interp_type==itp_steffen) {
1606  } else {
1607  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1608  o2scl::szttos(interp_type)+", in "+
1609  "interp::interp().").c_str(),exc_einval);
1610  }
1611  }
1612 
1613  virtual ~interp() {
1614  delete itp;
1615  }
1616 
1617  /// Give the value of the function \f$ y(x=x_0) \f$ .
1618  virtual double eval(const double x0, size_t n, const vec_t &x,
1619  const vec2_t &y) {
1620  itp->set(n,x,y);
1621  return itp->eval(x0);
1622  }
1623 
1624  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1625  virtual double deriv(const double x0, size_t n, const vec_t &x,
1626  const vec2_t &y) {
1627  itp->set(n,x,y);
1628  return itp->deriv(x0);
1629  }
1630 
1631  /** \brief Give the value of the second derivative
1632  \f$ y^{\prime \prime}(x=x_0) \f$ .
1633  */
1634  virtual double deriv2(const double x0, size_t n, const vec_t &x,
1635  const vec2_t &y) {
1636  itp->set(n,x,y);
1637  return itp->deriv2(x0);
1638  }
1639 
1640  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1641  virtual double integ(const double x1, const double x2, size_t n,
1642  const vec_t &x, const vec2_t &y) {
1643  itp->set(n,x,y);
1644  return itp->integ(x1,x2);
1645  }
1646 
1647  /// Set base interpolation type
1648  void set_type(size_t interp_type) {
1649  delete itp;
1650  if (interp_type==itp_linear) {
1652  } else if (interp_type==itp_cspline) {
1654  } else if (interp_type==itp_cspline_peri) {
1656  } else if (interp_type==itp_akima) {
1658  } else if (interp_type==itp_akima_peri) {
1660  } else if (interp_type==itp_monotonic) {
1662  } else if (interp_type==itp_steffen) {
1664  } else {
1665  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1666  o2scl::szttos(interp_type)+", in "+
1667  "interp::set().").c_str(),exc_einval);
1668  }
1669  return;
1670  }
1671 
1672 #ifndef DOXYGEN_INTERNAL
1673 
1674  private:
1675 
1677  interp<vec_t,vec2_t>& operator=(const interp<vec_t,vec2_t>&);
1678 
1679 #endif
1680 
1681  };
1682 
1683  /** \brief Interpolation class for pre-specified vector
1684 
1685  See also the \ref intp_section section of the \o2 User's guide.
1686 
1687  This interpolation class is intended to be sufficiently general
1688  to handle any vector type. Interpolation of ublas vector like
1689  objects is performed with the default template parameters.
1690 
1691  This class does not double check the vector to ensure that
1692  all of the intervals for the abcissa are all increasing or
1693  all decreasing.
1694 
1695  The type of interpolation to be performed can be specified
1696  using the set_type() function. The default is cubic splines
1697  with natural boundary conditions.
1698 
1699  \future Make version which copies vectors
1700  rather than storing pointers might be better and then
1701  has copy constructors.
1702  */
1703  template<class vec_t=boost::numeric::ublas::vector<double>,
1704  class vec2_t=vec_t> class interp_vec :
1705  public interp_base<vec_t,vec2_t> {
1706 
1707 #ifndef DOXYGEN_INTERNAL
1708 
1709  protected:
1710 
1711  /// Base interpolation object
1713 
1714  /// Interpolation type
1715  size_t itype;
1716 
1717 #endif
1718 
1719  public:
1720 
1721  /// Blank interpolator
1723  itp=0;
1724  itype=itp_cspline;
1725  }
1726 
1727  /// Create with base interpolation object \c it
1728  interp_vec(size_t n, const vec_t &x,
1729  const vec2_t &y, size_t interp_type=itp_cspline) {
1730 
1731  if (x[0]==x[n-1]) {
1732  O2SCL_ERR((((std::string)"Vector endpoints equal (value=")+
1733  o2scl::dtos(x[0])+") in interp_vec()::"+
1734  "interp_vec().").c_str(),exc_einval);
1735  }
1736 
1737  if (interp_type==itp_linear) {
1739  } else if (interp_type==itp_cspline) {
1741  } else if (interp_type==itp_cspline_peri) {
1743  } else if (interp_type==itp_akima) {
1745  } else if (interp_type==itp_akima_peri) {
1747  } else if (interp_type==itp_monotonic) {
1749  } else if (interp_type==itp_steffen) {
1751  } else {
1752  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1753  o2scl::szttos(interp_type)+", in "+
1754  "interp_vec::interp_vec().").c_str(),exc_einval);
1755  }
1756  itype=interp_type;
1757 
1758  itp->set(n,x,y);
1759  }
1760 
1761  virtual ~interp_vec() {
1762  if (itp!=0) {
1763  delete itp;
1764  }
1765  }
1766 
1767  /** \brief Set a new vector to interpolate
1768  */
1769  void set(size_t n, const vec_t &x, const vec2_t &y) {
1770  set(n,x,y,itype);
1771  return;
1772  }
1773 
1774  /** \brief Set a new vector to interpolate
1775  */
1776  void set(size_t n, const vec_t &x,
1777  const vec2_t &y, size_t interp_type) {
1778 
1779  if (x[0]==x[n-1]) {
1780  O2SCL_ERR((((std::string)"Vector endpoints equal (value=")+
1781  o2scl::dtos(x[0])+") in interp_vec()::"+
1782  "interp_vec().").c_str(),exc_einval);
1783  }
1784 
1785  delete itp;
1786  if (interp_type==itp_linear) {
1788  } else if (interp_type==itp_cspline) {
1790  } else if (interp_type==itp_cspline_peri) {
1792  } else if (interp_type==itp_akima) {
1794  } else if (interp_type==itp_akima_peri) {
1796  } else if (interp_type==itp_monotonic) {
1798  } else if (interp_type==itp_steffen) {
1800  } else {
1801  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1802  o2scl::szttos(interp_type)+", in "+
1803  "interp_vec::set().").c_str(),exc_einval);
1804  }
1805  itype=interp_type;
1806 
1807  itp->set(n,x,y);
1808  }
1809 
1810  /** \brief Manually clear the pointer to the user-specified vector
1811  */
1812  void clear() {
1813  if (itp!=0) {
1814  delete itp;
1815  itp=0;
1816  }
1817  return;
1818  }
1819 
1820  /// Give the value of the function \f$ y(x=x_0) \f$ .
1821  virtual double eval(const double x0) const {
1822  if (itp==0) {
1823  O2SCL_ERR("No vector set in interp_vec::eval().",
1824  exc_einval);
1825  }
1826  return itp->eval(x0);
1827  }
1828 
1829  /// Give the value of the function \f$ y(x=x_0) \f$ .
1830  virtual double operator()(double x0) const {
1831  if (itp==0) {
1832  O2SCL_ERR("No vector set in interp_vec::operator().",
1833  exc_einval);
1834  }
1835  return itp->eval(x0);
1836  }
1837 
1838  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1839  virtual double deriv(const double x0) const {
1840  if (itp==0) {
1841  O2SCL_ERR("No vector set in interp_vec::deriv().",
1842  exc_einval);
1843  }
1844  return itp->deriv(x0);
1845  }
1846 
1847  /** \brief Give the value of the second derivative
1848  \f$ y^{\prime \prime}(x=x_0) \f$ .
1849  */
1850  virtual double deriv2(const double x0) const {
1851  if (itp==0) {
1852  O2SCL_ERR("No vector set in interp_vec::deriv2().",
1853  exc_einval);
1854  }
1855  return itp->deriv2(x0);
1856  }
1857 
1858  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1859  virtual double integ(const double x1, const double x2) const {
1860  if (itp==0) {
1861  O2SCL_ERR("No vector set in interp_vec::integ().",
1862  exc_einval);
1863  }
1864  return itp->integ(x1,x2);
1865  }
1866 
1867  /// Return the type, "interp_vec"
1868  virtual const char *type() const {
1869  return "interp_vec";
1870  }
1871 
1872 #ifndef DOXYGEN_INTERNAL
1873 
1874  private:
1875 
1877  interp_vec<vec_t,vec2_t> &operator=
1878  (const interp_vec<vec_t,vec2_t> &it);
1879 
1880 #endif
1881 
1882  };
1883 
1884  /** \brief A specialization of interp for C-style double arrays
1885 
1886  See also the \ref intp_section section of the \o2 User's guide.
1887  */
1888  template<size_t n> class interp_array :
1889  public interp<double[n]> {
1890 
1891  public:
1892 
1893  /// Create with base interpolation objects \c it and \c rit
1894  interp_array(size_t interp_type)
1895  : interp<double[n]>(interp_type) {}
1896 
1897  /** \brief Create an interpolator using the default base
1898  interpolation objects
1899  */
1900  interp_array() : interp<double[n]>() {}
1901 
1902  };
1903 
1904  /** \brief A specialization of \ref o2scl::interp_vec for C-style arrays
1905 
1906  See also the \ref intp_section section of the \o2 User's guide.
1907  */
1908  template<class arr_t> class interp_array_vec :
1909  public interp_vec<arr_t> {
1910 
1911  public:
1912 
1913  /// Create with base interpolation object \c it
1914  interp_array_vec(size_t nv, const arr_t &x, const arr_t &y,
1915  size_t interp_type) :
1916  interp_vec<arr_t>(nv,x,y,interp_type) {}
1917  };
1918 
1919  /** \brief Count level crossings
1920 
1921  This function counts the number of times the function \f$ y(x) =
1922  \mathrm{level} \f$ where the function is defined by the vectors
1923  \c x and \c y.
1924 
1925  The value returned is exactly the same as the size of the
1926  \c locs vector computed by \ref vector_find_level().
1927 
1928  This function will call the error handler if \c n is
1929  less than two.
1930 
1931  If one of the entries in \c y is either larger or smaller than
1932  its neighbors (i.e. if the function \f$ y(x) \f$ has an
1933  extremum), and if the value of \c level is exactly equal to the
1934  extremum, then this is counted as 1 level crossing. The same
1935  applies if either the first or last entry in \c y is exactly
1936  equal to \c level . Finite-precision errors may affect
1937  the result of this function when \c level is nearly
1938  equal to one of the value in the vector \c y.
1939  */
1940  template<class vec_t, class vec2_t> size_t vector_level_count
1941  (double level, size_t n, vec_t &x, vec2_t &y) {
1942 
1943  if (n<=1) {
1944  O2SCL_ERR2("Need at least two data points in ",
1945  "vector_find_count().",exc_einval);
1946  }
1947 
1948  size_t count=0;
1949 
1950  // Handle left boundary
1951  if (y[0]==level) count++;
1952 
1953  // Find intersections
1954  for(size_t i=0;i<n-1;i++) {
1955 
1956  if ((y[i]<level && y[i+1]>=level) ||
1957  (y[i]>level && y[i+1]<=level)) {
1958  count++;
1959  }
1960  }
1961 
1962  return count;
1963  }
1964 
1965  /** \brief Perform inverse linear interpolation
1966 
1967  This function performs inverse linear interpolation of the data
1968  defined by \c x and \c y, finding all points in \c x which have
1969  the property \f$ \mathrm{level} = y(x) \f$. All points for which
1970  this relation holds are put into the vector \c locs. The
1971  previous information contained in vector \c locs before the
1972  function call is destroyed.
1973 
1974  This is the 1-dimensional analog of finding contour levels as
1975  the \ref contour class does for 2 dimensions.
1976 
1977  This function will call the error handler if \c n is
1978  less than two.
1979 
1980  This function is inclusive of the endpoints, so that if either
1981  <tt>y[0]</tt> and/or <tt>y[n-1]</tt> is equal to level, then
1982  <tt>x[0]</tt> and/or <tt>x[n-1]</tt> are added to \c locs,
1983  respectively.
1984  */
1985  template<class vec_t, class vec2_t> void vector_find_level
1986  (double level, size_t n, vec_t &x, vec2_t &y, std::vector<double> &locs) {
1987 
1988  if (n<=1) {
1989  O2SCL_ERR2("Need at least two data points in ",
1990  "vector_find_level().",exc_einval);
1991  }
1992 
1993  // Ensure that the location vector is empty
1994  locs.clear();
1995 
1996  // Handle left boundary
1997  if (y[0]==level) {
1998  locs.push_back(x[0]);
1999  }
2000 
2001  // Find intersections
2002  for(size_t i=0;i<n-1;i++) {
2003 
2004  if ((y[i]<level && y[i+1]>level) ||
2005  (y[i]>level && y[i+1]<level)) {
2006  // For each intersection, add the location using linear
2007  // interpolation
2008  double x0=x[i]+(x[i+1]-x[i])*(level-y[i])/(y[i+1]-y[i]);
2009  locs.push_back(x0);
2010  } else if (y[i+1]==level) {
2011  locs.push_back(x[i+1]);
2012  }
2013  }
2014 
2015  return;
2016  }
2017 
2018  /** \brief Compute derivative at all points from an
2019  interpolation object
2020  */
2021  template<class ovec_t, class vec2_t>
2022  void vector_deriv_interp(size_t n, ovec_t &v, vec2_t &dv,
2023  size_t interp_type=itp_linear) {
2024  ovec_t grid(n);
2025  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2026  interp_vec<ovec_t,ovec_t> oi(n,grid,v,interp_type);
2027  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(((double)i));
2028  return;
2029  }
2030 
2031  /** \brief Compute second derivative at all points from an
2032  interpolation object
2033  */
2034  template<class ovec_t, class vec2_t>
2035  void vector_deriv2_interp(size_t n, ovec_t &v, vec2_t &dv,
2036  size_t interp_type=itp_linear) {
2037  ovec_t grid(n);
2038  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2039  interp_vec<ovec_t,ovec_t> oi(n,grid,v,interp_type);
2040  for(size_t i=0;i<n;i++) dv[i]=oi.deriv2(((double)i));
2041  return;
2042  }
2043 
2044  /** \brief Compute derivative at all points from an
2045  interpolation object
2046  */
2047  template<class vec_t, class vec2_t, class vec3_t>
2048  void vector_deriv_xy_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv,
2049  size_t interp_type=itp_linear) {
2050  interp_vec<vec_t,vec2_t> oi(n,vx,vy,interp_type);
2051  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(vx[i]);
2052  return;
2053  }
2054 
2055  /** \brief Compute second derivative at all points from an
2056  interpolation object
2057  */
2058  template<class vec_t, class vec2_t, class vec3_t>
2059  void vector_deriv2_xy_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv,
2060  size_t interp_type=itp_linear) {
2061  interp_vec<vec_t,vec2_t> oi(n,vx,vy,interp_type);
2062  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(vx[i]);
2063  return;
2064  }
2065 
2066  /** \brief Integral of a vector from interpolation object
2067  */
2068  template<class ovec_t>
2069  double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type) {
2070  ovec_t grid(n);
2071  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2072  interp_vec<ovec_t> oi(n,grid,v,interp_type);
2073  return oi.integ(0.0,((double)(n-1)));
2074  }
2075 
2076  /** \brief Compute the integral over <tt>y(x)</tt> using
2077  interpolation
2078 
2079  Assuming vectors \c y and \c x define a function \f$ y(x) \f$
2080  then this computes
2081  \f[
2082  \int_{x_0}^{x^{n-1}} y(x)~dx
2083  \f]
2084  using interpolation to approximate the result.
2085 
2086  See also \ref vector_deriv_interp() and
2087  \ref vector_deriv2_interp() in \ref vector_derint.h .
2088  */
2089  template<class vec_t, class vec2_t>
2090  double vector_integ_xy_interp(size_t n, const vec_t &x, const vec2_t &y,
2091  size_t interp_type=itp_linear) {
2092 
2093  // Interpolation object
2094  interp<vec_t,vec2_t> si(interp_type);
2095 
2096  // Compute full integral
2097  double total=si.integ(x[0],x[n-1],n,x,y);
2098 
2099  return total;
2100  }
2101 
2102  /** \brief Compute integral over <tt>y(x)</tt> and store result
2103  in a vector using interpolation
2104  */
2105  template<class vec_t, class vec2_t, class vec3_t>
2106  void vector_integ_xy_interp(size_t n, const vec_t &x, const vec2_t &y,
2107  vec3_t &iy, size_t interp_type=itp_linear) {
2108 
2109  // Interpolation object
2110  interp<vec_t,vec2_t> si(interp_type);
2111 
2112  // Compute full integral
2113  iy[0]=0.0;
2114  for(size_t i=1;i<n;i++) {
2115  iy[i]=si.integ(x[0],x[i],n,x,y);
2116  }
2117 
2118  return;
2119  }
2120 
2121  /** \brief Compute the integral of a vector using
2122  interpolation up to a specified upper limit
2123  */
2124  template<class ovec_t>
2125  double vector_integ_ul_interp(size_t n, double x2,
2126  ovec_t &v, size_t interp_type) {
2127  ovec_t grid(n);
2128  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2129  interp_vec<ovec_t> oi(n,grid,v,interp_type);
2130  return oi.integ(0.0,x2);
2131  }
2132 
2133  /** \brief Compute the integral over <tt>y(x)</tt> using
2134  interpolation up to a specified upper limit
2135  */
2136  template<class vec_t, class vec2_t>
2137  double vector_integ_ul_xy_interp(size_t n, const vec_t &x, const vec2_t &y,
2138  double x2, size_t interp_type=itp_linear) {
2139 
2140  // Interpolation object
2141  interp<vec_t,vec2_t> si(interp_type);
2142 
2143  // Compute full integral
2144  double total=si.integ(x[0],x2,n,x,y);
2145 
2146  return total;
2147  }
2148 
2149  /** \brief Compute the endpoints which enclose the regions whose
2150  integral is equal to \c sum
2151 
2152  Defining a new function, \f$ g(y_0) \f$ which takes as input any
2153  y-value, \f$ y_0 \f$ from the function \f$ y(x) \f$ (specified
2154  with the parameters \c x and \c y) and outputs the integral of
2155  the function \f$ y(x) \f$ over all regions where \f$ y(x) > y_0
2156  \f$. This function inverts \f$ g(y) \f$, taking the value
2157  of an integral as input, and returns the corresponding y-value
2158  in the variable \c lev.
2159 
2160  In order to make sure that the intepretation of the integral is
2161  unambiguous, this function requires that the first and last values
2162  of \c y are equal, i.e. <tt>y[0]==y[n-1]</tt>. This function
2163  also requires at least two data points, i.e. \c n cannot be
2164  0 or 1.
2165 
2166  This function is particularly useful, for example, in computing
2167  the region which defines 68\% around a peak of data, thus
2168  providing \f$ 1~\sigma \f$ confidence limits.
2169 
2170  Linear interpolation is used to describe the function \f$ g \f$,
2171  and the precision of this function is limited by this assumption.
2172  This function may also sometimes fail if \c sum is very close to
2173  the minimum or maximum value of the function \f$ g \f$.
2174 
2175  \todo This function may also require that all of the
2176  y-vector values have the same sign or are all positive.
2177  Check this.
2178 
2179  \comment
2180  Note that the two vector types for x and y must be the
2181  same in order to use o2scl_interp.
2182  \endcomment
2183  */
2184  template<class vec_t, class vec2_t> int vector_invert_enclosed_sum
2185  (double sum, size_t n, vec_t &x, vec2_t &y, double &lev,
2186  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2187 
2188  if (n<=1) {
2189  O2SCL_ERR2("Need at least two data points in ",
2190  "vector_invert_enclosed_sum().",exc_einval);
2191  }
2192 
2194 
2195  // Handle boundaries
2196  std::vector<double> x2, y2;
2197  size_t n2;
2198  if (boundaries==1) {
2199  x2.resize(n+1);
2200  y2.resize(n+1);
2201  x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2202  y2[0]=0.0;
2203  for(size_t i=0;i<n;i++) {
2204  x2[i+1]=x[i];
2205  y2[i+1]=y[i];
2206  }
2207  n2=n+1;
2208  } else if (boundaries==2) {
2209  x2.resize(n+1);
2210  y2.resize(n+1);
2211  for(size_t i=0;i<n;i++) {
2212  x2[i]=x[i];
2213  y2[i]=y[i];
2214  }
2215  x2[n]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2216  y2[n]=0.0;
2217  n2=n+1;
2218  } else if (boundaries==3) {
2219  x2.resize(n+2);
2220  y2.resize(n+2);
2221  x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2222  y2[0]=0.0;
2223  for(size_t i=0;i<n;i++) {
2224  x2[i+1]=x[i];
2225  y2[i+1]=y[i];
2226  }
2227  x2[n+1]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2228  y2[n+1]=0.0;
2229  n2=n+2;
2230  } else {
2231  x2.resize(n);
2232  y2.resize(n);
2233  for(size_t i=0;i<n;i++) {
2234  x2[i]=x[i];
2235  y2[i]=y[i];
2236  }
2237  n2=n;
2238  }
2239 
2240  // Construct a sorted list of function values
2241  ubvector ysort(n2);
2242  vector_copy(n2,y2,ysort);
2243  vector_sort<ubvector,double>(ysort.size(),ysort);
2244 
2245  // Create list of y-values to perform y-value and integral
2246  // interpolation. We choose values in between the grid points to
2247  // ensure that we don't accidentally land precisely on a peak or
2248  // valley.
2249  std::vector<double> ylist;
2250  for(size_t i=0;i<ysort.size()-1;i++) {
2251  if (ysort[i]!=ysort[i+1]) {
2252  ylist.push_back((ysort[i+1]+3.0*ysort[i])/4.0);
2253  ylist.push_back((ysort[i+1]*3.0+ysort[i])/4.0);
2254  }
2255  }
2256 
2257  // Interpolation objects. We need two, one for the
2258  // user-specified vector type, and one for std::vector<double>
2260  interp<std::vector<double>,std::vector<double> > itp2(itp_linear);
2261 
2262  // Construct vectors for interpolation
2263  std::vector<double> xi, yi;
2264 
2265  size_t nfail=0;
2266 
2267  for(size_t k=0;k<ylist.size();k++) {
2268  double lev_tmp=ylist[k];
2269  std::vector<double> locs;
2270  vector_find_level(lev_tmp,n2,x2,y2,locs);
2271  if ((locs.size()%2)!=0) {
2272  nfail++;
2273  if (verbose>0) {
2274  std::cout << lev_tmp << " " << 0.0 << " "
2275  << locs.size() << " (fail)" << std::endl;
2276  }
2277  } else {
2278  double sum_temp=0.0;
2279  for(size_t i=0;i<locs.size()/2;i++) {
2280  double x0=locs[2*i];
2281  double x1=locs[2*i+1];
2282  sum_temp+=itp.integ(x0,x1,n2,x2,y2);
2283  }
2284  xi.push_back(sum_temp);
2285  yi.push_back(lev_tmp);
2286  if (verbose>0) {
2287  std::cout << lev_tmp << " " << sum_temp << " "
2288  << locs.size() << " ";
2289  for(size_t i=0;i<locs.size();i++) {
2290  std::cout << locs[i] << " ";
2291  }
2292  std::cout << std::endl;
2293  }
2294  }
2295  }
2296  if (nfail>10) {
2297  if (err_on_fail) {
2298  O2SCL_ERR2("More than 10 failures in ",
2299  "vector_invert_enclosed_sum().",o2scl::exc_efailed);
2300  }
2301  return o2scl::exc_efailed;
2302  }
2303 
2304  lev=itp2.eval(sum,xi.size(),xi,yi);
2305 
2306  return 0;
2307  }
2308 
2309  /** \brief Find the region enclosing an integral
2310  */
2311  template<class vec_t, class vec2_t> int vector_region_int
2312  (size_t n, vec_t &x, vec2_t &y, double intl, std::vector<double> &locs,
2313  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2314 
2315  // Total integral
2316  double total=vector_integ_xy_interp(n,x,y,itp_linear);
2317  if (verbose>0) {
2318  std::cout << "Total integral: " << total << std::endl;
2319  }
2320  // Specified fractional integral
2321  if (verbose>0) {
2322  std::cout << "Target integral: " << intl << std::endl;
2323  }
2324  // Find correct level
2325  double lev;
2326  int ret=vector_invert_enclosed_sum(intl,n,x,y,lev,
2327  boundaries,verbose,err_on_fail);
2328  if (ret!=0) {
2329  if (err_on_fail) {
2330  O2SCL_ERR2("Failed to find a level which enclosed the ",
2331  "specified integral in vector_region_int().",
2333  }
2334  return o2scl::exc_efailed;
2335  }
2336  if (verbose>0) {
2337  std::cout << "Level from vector_invert: " << lev << std::endl;
2338  }
2339 
2340  // Inverse interpolate to find locations corresponding to level
2341  vector_find_level(lev,n,x,y,locs);
2342  if (verbose>0) {
2343  std::cout << "Locations from vector_find_level: ";
2344  for(size_t i=0;i<locs.size();i++) {
2345  std::cout << locs[i];
2346  if (i!=locs.size()-1) {
2347  std::cout << " ";
2348  }
2349  }
2350  std::cout << std::endl;
2351  }
2352  return 0;
2353  }
2354 
2355  /** \brief Find the region enclosing a partial integral
2356  */
2357  template<class vec_t, class vec2_t> int vector_region_fracint
2358  (size_t n, vec_t &x, vec2_t &y, double frac, std::vector<double> &locs,
2359  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2360  double total=vector_integ_xy_interp(n,x,y,itp_linear);
2361  return vector_region_int(n,x,y,frac*total,locs,boundaries,
2362  verbose,err_on_fail);
2363  }
2364 
2365  /** \brief Find the boundaries of the region enclosing a integral
2366 
2367  This function finds the boundaries of the region which
2368  has integral equal to <tt>frac</tt> times the full
2369  integral from the lower x limit to the upper x limit.
2370  */
2371  template<class vec_t, class vec2_t> int vector_bound_fracint
2372  (size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high,
2373  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2374 
2375  std::vector<double> locs;
2376  int ret=vector_region_fracint(n,x,y,frac,locs,boundaries,
2377  verbose,err_on_fail);
2378  if (locs.size()==0 || ret!=0) {
2379  if (err_on_fail) {
2380  O2SCL_ERR2("Zero level crossings or vector_region_fracint() ",
2381  "failed in vector_bound_sigma().",exc_efailed);
2382  }
2383  return o2scl::exc_efailed;
2384  }
2385  // Return min and max location
2386  size_t ix;
2387  vector_min(locs.size(),locs,ix,low);
2388  vector_max(locs.size(),locs,ix,high);
2389  return 0;
2390  }
2391 
2392  /** \brief Find the boundaries of the region enclosing a integral
2393 
2394  This function finds the boundaries of the region which
2395  has integral equal to <tt>frac</tt> times the full
2396  integral from the lower x limit to the upper x limit.
2397  */
2398  template<class vec_t, class vec2_t> int vector_bound_int
2399  (size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high,
2400  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2401 
2402  std::vector<double> locs;
2403  int ret=vector_region_int(n,x,y,frac,locs,boundaries,
2404  verbose,err_on_fail);
2405  if (locs.size()==0 || ret!=0) {
2406  if (err_on_fail) {
2407  O2SCL_ERR2("Zero level crossings or vector_region_int() ",
2408  "failed in vector_bound_sigma().",exc_efailed);
2409  }
2410  return o2scl::exc_efailed;
2411  }
2412  // Return min and max location
2413  size_t ix;
2414  vector_min(locs.size(),locs,ix,low);
2415  vector_max(locs.size(),locs,ix,high);
2416  return 0;
2417  }
2418 
2419 #ifndef DOXYGEN_NO_O2NS
2420 }
2421 #endif
2422 
2423 #endif
virtual const char * type() const
Return the type, "interp_steffen".
Definition: interp.h:1280
virtual const char * type() const
Return the type, "interp_cspline_peri".
Definition: interp.h:634
Interpolation class for pre-specified vector.
Definition: interp.h:1704
o2scl_linalg::ubvector_5_mem p5m
Memory for the tridiagonalization.
Definition: interp.h:616
void vector_deriv2_xy_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv, size_t interp_type=itp_linear)
Compute second derivative at all points from an interpolation object.
Definition: interp.h:2059
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:1618
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:1728
virtual const char * type() const
Return the type, "interp_linear".
Definition: interp.h:326
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:1591
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 vector_deriv_interp(size_t n, ovec_t &v, vec2_t &dv, size_t interp_type=itp_linear)
Compute derivative at all points from an interpolation object.
Definition: interp.h:2022
void akima_calc(const vec_t &x_array, size_t size, ubvector &umx)
For initializing the interpolation.
Definition: interp.h:776
int vector_invert_enclosed_sum(double sum, size_t n, vec_t &x, vec2_t &y, double &lev, int boundaries=0, int verbose=0, bool err_on_fail=true)
Compute the endpoints which enclose the regions whose integral is equal to sum.
Definition: interp.h:2185
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:1888
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:1438
const double pi
Definition: constants.h:62
interp_base< vec_t, vec2_t > * itp
Pointer to base interpolation object.
Definition: interp.h:1584
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:1850
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:379
Akima spline interpolation with periodic boundary conditions (GSL)
Definition: interp.h:979
virtual double integ(double al, double bl) const
Give the value of the integral .
Definition: interp.h:1229
void clear()
Manually clear the pointer to the user-specified vector.
Definition: interp.h:1812
virtual const char * type() const
Return the type, "interp_cspline".
Definition: interp.h:587
virtual const char * type() const
Return the type, "interp_akima".
Definition: interp.h:961
double vector_integ_ul_xy_interp(size_t n, const vec_t &x, const vec2_t &y, double x2, size_t interp_type=itp_linear)
Compute the integral over y(x) using interpolation up to a specified upper limit. ...
Definition: interp.h:2137
int vector_bound_fracint(size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high, int boundaries=0, int verbose=0, bool err_on_fail=true)
Find the boundaries of the region enclosing a integral.
Definition: interp.h:2372
virtual const char * type() const
Return the type, "interp_monotonic".
Definition: interp.h:1550
double copysign(const double x, const double y)
Flip the sign of x if x and y have different signs.
Definition: interp.h:1089
size_t sz
Vector size.
Definition: interp.h:128
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:1415
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:1191
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:878
Monotonicity-preserving interpolation.
Definition: interp.h:1318
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:1914
virtual const char * type() const
Return the type, "interp_akima_peri".
Definition: interp.h:999
double vector_integ_xy_interp(size_t n, const vec_t &x, const vec2_t &y, size_t interp_type=itp_linear)
Compute the integral over y(x) using interpolation.
Definition: interp.h:2090
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:1205
virtual double deriv(const double x0) const
Give the value of the derivative .
Definition: interp.h:1839
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:532
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:1464
interp_vec()
Blank interpolator.
Definition: interp.h:1722
void vector_deriv2_interp(size_t n, ovec_t &v, vec2_t &dv, size_t interp_type=itp_linear)
Compute second derivative at all points from an interpolation object.
Definition: interp.h:2035
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:909
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:1218
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:895
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:1487
Steffen&#39;s monotonicity-preserving interpolation.
Definition: interp.h:1059
interp_akima()
Create a base interpolation object with or without periodic boundary conditions.
Definition: interp.h:814
Cubic spline interpolation with periodic boundary conditions (GSL)
Definition: interp.h:606
virtual double operator()(double x0) const
Give the value of the function .
Definition: interp.h:1830
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:488
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:1641
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:275
Cubic spline for natural boundary conditions.
Definition: interp.h:72
double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type)
Integral of a vector from interpolation object.
Definition: interp.h:2069
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:252
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
ubvector beta
Staggered ratio.
Definition: interp.h:1338
Allocation object for 5 arrays of equal size.
Definition: tridiag_base.h:86
virtual const char * type() const
Return the type, "interp_vec".
Definition: interp.h:1868
interp_array()
Create an interpolator using the default base interpolation objects.
Definition: interp.h:1900
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
int vector_region_fracint(size_t n, vec_t &x, vec2_t &y, double frac, std::vector< double > &locs, int boundaries=0, int verbose=0, bool err_on_fail=true)
Find the region enclosing a partial integral.
Definition: interp.h:2358
size_t itype
Interpolation type.
Definition: interp.h:1715
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:1859
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:863
const vec2_t * py
Dependent vector.
Definition: interp.h:125
o2scl_linalg::ubvector_4_mem p4m
Memory for the tridiagonalization.
Definition: interp.h:376
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
void vector_deriv_xy_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv, size_t interp_type=itp_linear)
Compute derivative at all points from an interpolation object.
Definition: interp.h:2048
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:1648
size_t vector_level_count(double level, size_t n, vec_t &x, vec2_t &y)
Count level crossings.
Definition: interp.h:1941
Akima spline for natural boundary conditions.
Definition: interp.h:76
int vector_bound_int(size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high, int boundaries=0, int verbose=0, bool err_on_fail=true)
Find the boundaries of the region enclosing a integral.
Definition: interp.h:2399
Akima spline interpolation (GSL)
Definition: interp.h:752
interp_steffen()
Create a base interpolation object.
Definition: interp.h:1101
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:511
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
int vector_region_int(size_t n, vec_t &x, vec2_t &y, double intl, std::vector< double > &locs, int boundaries=0, int verbose=0, bool err_on_fail=true)
Find the region enclosing an integral.
Definition: interp.h:2312
interp_cspline()
Create a base interpolation object with natural or periodic boundary conditions.
Definition: interp.h:398
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:1634
ubvector Delta
Finite differences.
Definition: interp.h:1334
Cubic spline interpolation (GSL)
Definition: interp.h:348
interp_array(size_t interp_type)
Create with base interpolation objects it and rit.
Definition: interp.h:1894
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:107
double vector_integ_ul_interp(size_t n, double x2, ovec_t &v, size_t interp_type)
Compute the integral of a vector using interpolation up to a specified upper limit.
Definition: interp.h:2125
search_vec< const vec_t > svx
To perform binary searches.
Definition: interp.h:119
ubvector alpha
Ratio.
Definition: interp.h:1336
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:1821
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:270
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:1986
Linear.
Definition: interp.h:70
ubvector m
Slopes.
Definition: interp.h:1332
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:467
A specialization of o2scl::interp_vec for C-style arrays.
Definition: interp.h:1908
interp_base< vec_t, vec2_t > * itp
Base interpolation object.
Definition: interp.h:1712
Interpolation class for general vectors.
Definition: interp.h:1577
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:1625

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