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 #include <o2scl/vector.h>
62 
63 #ifndef DOXYGEN_NO_O2NS
64 namespace o2scl {
65 #endif
66 
67  /** \brief Interpolation types
68  */
69  enum {
70  /// Linear
72  /// Cubic spline for natural boundary conditions
74  /// Cubic spline for periodic boundary conditions
76  /// Akima spline for natural boundary conditions
78  /// Akima spline for periodic boundary conditions
80  /// Monotonicity-preserving interpolation (unfinished)
82  /// Steffen's monotonic method
84  /// Nearest-neighbor lookup
86  };
87 
88  /** \brief Base low-level interpolation class [abstract base]
89 
90  See also the \ref intp_section section of the \o2 User's guide.
91 
92  To interpolate a set vectors \c x and \c y, call set() and then
93  the interpolation functions eval(), deriv(), deriv2() and
94  integ(). If the \c x and \c y vectors do not change, then you
95  may call the interpolation functions multiple times in
96  succession. These classes do not copy the user-specified vectors
97  but store pointers to them. Thus, if the vector is changed
98  without a new call to \ref interp_base::set(), the behavior of
99  the interpolation functions is undefined.
100 
101  \comment
102  AWS, 12/27/13: Copy constructors might be ill-advised for
103  this class since we store pointers. For now, we don't
104  allow the user to use them.
105  \endcomment
106  */
107  template<class vec_t, class vec2_t=vec_t> class interp_base {
108 
109 #ifdef O2SCL_NEVER_DEFINED
110  }{
111 #endif
112 
113 #ifndef DOXYGEN_INTERNAL
114 
115  protected:
116 
117  /** \brief To perform binary searches
118 
119  This pointer is set to zero in the constructor and should be
120  non-zero only if it has been allocated with \c new.
121  */
123 
124  /// Independent vector
125  const vec_t *px;
126 
127  /// Dependent vector
128  const vec2_t *py;
129 
130  /// Vector size
131  size_t sz;
132 
133  /** \brief An internal function to assist in computing the
134  integral for both the cspline and Akima types
135  */
136  double integ_eval(double ai, double bi, double ci, double di, double xi,
137  double a, double b) const {
138 
139  double r1=a-xi;
140  double r2=b-xi;
141  double r12=r1+r2;
142  double bterm=0.5*bi*r12;
143  double cterm=ci*(r1*r1+r2*r2+r1*r2)/3.0;
144  double dterm=0.25*di*r12*(r1*r1+r2*r2);
145 
146  return (b-a)*(ai+bterm+cterm+dterm);
147  }
148 
149 #endif
150 
151  public:
152 
153  interp_base() {
154  sz=0;
155  }
156 
157  virtual ~interp_base() {
158  }
159 
160  /** \brief The minimum size of the vectors to interpolate between
161 
162  This variable must be set in the constructor of the children
163  for access by the class user.
164  */
165  size_t min_size;
166 
167  /// Initialize interpolation routine
168  virtual void set(size_t size, const vec_t &x, const vec2_t &y)=0;
169 
170  /// Give the value of the function \f$ y(x=x_0) \f$
171  virtual double eval(double x0) const=0;
172 
173  /// Give the value of the function \f$ y(x=x_0) \f$ .
174  virtual double operator()(double x0) const {
175  return eval(x0);
176  }
177 
178  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
179  virtual double deriv(double x0) const=0;
180 
181  /** \brief Give the value of the second derivative
182  \f$ y^{\prime \prime}(x=x_0) \f$ .
183  */
184  virtual double deriv2(double x0) const=0;
185 
186  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
187  virtual double integ(double a, double b) const=0;
188 
189  /// Return the type
190  virtual const char *type() const=0;
191 
192 #ifndef DOXYGEN_INTERNAL
193 
194  private:
195 
198 
199 #endif
200 
201  };
202 
203  /** \brief Linear interpolation (GSL)
204 
205  See also the \ref intp_section section of the \o2 User's guide.
206 
207  Linear interpolation requires no calls to allocate() or free()
208  as there is no internal storage required.
209  */
210  template<class vec_t, class vec2_t=vec_t> class interp_linear :
211  public interp_base<vec_t,vec2_t> {
212 
213 #ifdef O2SCL_NEVER_DEFINED
214  }{
215 #endif
216 
217  public:
218 
219  interp_linear() {
220  this->min_size=2;
221  }
222 
223  virtual ~interp_linear() {}
224 
225  /// Initialize interpolation routine
226  virtual void set(size_t size, const vec_t &x, const vec2_t &y) {
227  if (size<this->min_size) {
228  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
229  " than "+szttos(this->min_size)+" in interp_linear::"+
230  "set().").c_str(),exc_einval);
231  }
232  this->svx.set_vec(size,x);
233  this->px=&x;
234  this->py=&y;
235  this->sz=size;
236  return;
237  }
238 
239  /// Give the value of the function \f$ y(x=x_0) \f$ .
240  virtual double eval(double x0) const {
241 
242  size_t cache=0;
243  size_t index=this->svx.find_const(x0,cache);
244 
245  double x_lo=(*this->px)[index];
246  double x_hi=(*this->px)[index+1];
247  double y_lo=(*this->py)[index];
248  double y_hi=(*this->py)[index+1];
249  double dx=x_hi-x_lo;
250 
251  return y_lo+(x0-x_lo)/dx*(y_hi-y_lo);
252  }
253 
254  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
255  virtual double deriv(double x0) const {
256 
257  size_t cache=0;
258  size_t index=this->svx.find_const(x0,cache);
259 
260  double x_lo=(*this->px)[index];
261  double x_hi=(*this->px)[index+1];
262  double y_lo=(*this->py)[index];
263  double y_hi=(*this->py)[index+1];
264  double dx=x_hi-x_lo;
265  double dy=y_hi-y_lo;
266 
267  return dy/dx;
268  }
269 
270  /** \brief Give the value of the second derivative
271  \f$ y^{\prime \prime}(x=x_0) \f$ (always zero)
272  */
273  virtual double deriv2(double x0) const {
274  return 0.0;
275  }
276 
277  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
278  virtual double integ(double a, double b) const {
279 
280  size_t cache=0;
281  size_t i, index_a, index_b;
282 
283  bool flip=false;
284  if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
285  ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
286  double tmp=a;
287  a=b;
288  b=tmp;
289  flip=true;
290  }
291 
292  index_a=this->svx.find_const(a,cache);
293  index_b=this->svx.find_const(b,cache);
294 
295  double result=0.0;
296  for(i=index_a; i<=index_b; i++) {
297 
298  double x_lo=(*this->px)[i];
299  double x_hi=(*this->px)[i+1];
300  double y_lo=(*this->py)[i];
301  double y_hi=(*this->py)[i+1];
302  double dx=x_hi-x_lo;
303 
304  if(dx != 0.0) {
305 
306  if (i == index_a || i == index_b) {
307  double x1=(i == index_a) ? a : x_lo;
308  double x2=(i == index_b) ? b : x_hi;
309  double D=(y_hi-y_lo)/dx;
310  result += (x2-x1)*(y_lo+0.5*D*((x2-x_lo)+(x1-x_lo)));
311  } else {
312  result += 0.5*dx*(y_lo+y_hi);
313  }
314 
315  } else {
316  std::string str=((std::string)"Interval of length zero ")+
317  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
318  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
319  " in interp_linear::integ().";
320  O2SCL_ERR(str.c_str(),exc_einval);
321  }
322  }
323 
324  if (flip) result=-result;
325  return result;
326  }
327 
328  /// Return the type, \c "interp_linear".
329  virtual const char *type() const { return "interp_linear"; }
330 
331 #ifndef DOXYGEN_INTERNAL
332 
333  private:
334 
337 
338 #endif
339 
340  };
341 
342  /** \brief Nearest-neighbor interpolation
343 
344  Nearest interpolation requires no calls to allocate() or free()
345  as there is no internal storage required.
346  */
347  template<class vec_t, class vec2_t=vec_t> class interp_nearest_neigh :
348  public interp_base<vec_t,vec2_t> {
349 
350 #ifdef O2SCL_NEVER_DEFINED
351  }{
352 #endif
353 
354  public:
355 
357  this->min_size=1;
358  }
359 
360  virtual ~interp_nearest_neigh() {}
361 
362  /// Initialize interpolation routine
363  virtual void set(size_t size, const vec_t &x, const vec2_t &y) {
364  if (size<this->min_size) {
365  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
366  " than "+szttos(this->min_size)+" in interp_nearest_neigh::"+
367  "set().").c_str(),exc_einval);
368  }
369  this->svx.set_vec(size,x);
370  this->px=&x;
371  this->py=&y;
372  this->sz=size;
373  return;
374  }
375 
376  /// Give the value of the function \f$ y(x=x_0) \f$ .
377  virtual double eval(double x0) const {
378 
379  size_t index=this->svx.ordered_lookup(x0);
380  return (*this->py)[index];
381  }
382 
383  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
384  virtual double deriv(double x0) const {
385  return 0.0;
386  }
387 
388  /** \brief Give the value of the second derivative
389  \f$ y^{\prime \prime}(x=x_0) \f$ (always zero)
390  */
391  virtual double deriv2(double x0) const {
392  return 0.0;
393  }
394 
395  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
396  virtual double integ(double a, double b) const {
397  return 0.0;
398  }
399 
400  /// Return the type, \c "interp_nearest_neigh".
401  virtual const char *type() const { return "interp_nearest_neigh"; }
402 
403 #ifndef DOXYGEN_INTERNAL
404 
405  private:
406 
410  (const interp_nearest_neigh<vec_t,vec2_t>&);
411 
412 #endif
413 
414  };
415 
416  /** \brief Cubic spline interpolation (GSL)
417 
418  See also the \ref intp_section section of the \o2 User's guide.
419 
420  By default, this class uses natural boundary conditions, where
421  the second derivative vanishes at each end point. Extrapolation
422  effectively assumes that the second derivative is linear outside
423  of the endpoints.
424  */
425  template<class vec_t, class vec2_t=vec_t> class interp_cspline :
426  public interp_base<vec_t,vec2_t> {
427 
428 #ifdef O2SCL_NEVER_DEFINED
429  }{
430 #endif
431 
432  public:
433 
435  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
436  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
437  typedef boost::numeric::ublas::slice slice;
438  typedef boost::numeric::ublas::range range;
439 
440 #ifndef DOXYGEN_INTERNAL
441 
442  protected:
443 
444  /// \name Storage for cubic spline interpolation
445  //@{
446  ubvector c;
447  ubvector g;
448  ubvector diag;
449  ubvector offdiag;
450  //@}
451 
452  /// Memory for the tridiagonalization
454 
455  /// Compute coefficients for cubic spline interpolation
456  void coeff_calc(const ubvector &c_array, double dy, double dx,
457  size_t index, double &b, double &c2, double &d) const {
458 
459  double c_i=c_array[index];
460  double c_ip1=c_array[index+1];
461  b=(dy/dx)-dx*(c_ip1+2.0*c_i)/3.0;
462  c2=c_i;
463  d=(c_ip1-c_i)/(3.0*dx);
464 
465  return;
466  }
467 
468 #endif
469 
470  public:
471 
472  /** \brief Create a base interpolation object with natural or
473  periodic boundary conditions
474  */
476  this->min_size=3;
477  }
478 
479  virtual ~interp_cspline() {
480  }
481 
482  /** \brief Initialize interpolation routine
483  */
484  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
485 
486  if (size<this->min_size) {
487  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
488  " than "+szttos(this->min_size)+" in interp_cspline::"+
489  "set().").c_str(),exc_einval);
490  }
491 
492  if (size!=this->sz) {
493  c.resize(size);
494  g.resize(size);
495  diag.resize(size);
496  offdiag.resize(size);
497  p4m.resize(size);
498  }
499 
500  this->px=&xa;
501  this->py=&ya;
502  this->sz=size;
503 
504  this->svx.set_vec(size,xa);
505 
506  /// Natural boundary conditions
507 
508  size_t i;
509  size_t num_points=size;
510  size_t max_index=num_points-1;
511  size_t sys_size=max_index-1;
512 
513  c[0]=0.0;
514  c[max_index]=0.0;
515 
516  for (i=0; i < sys_size; i++) {
517  double h_i=xa[i+1]-xa[i];
518  double h_ip1=xa[i+2]-xa[i+1];
519  double ydiff_i=ya[i+1]-ya[i];
520  double ydiff_ip1=ya[i+2]-ya[i+1];
521  double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
522  double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
523  offdiag[i]=h_ip1;
524  diag[i]=2.0*(h_ip1+h_i);
525  g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
526  }
527 
528  if (sys_size == 1) {
529 
530  c[1]=g[0]/diag[0];
531 
532  return;
533  }
534 
535  ubvector_range cp1(c,range(1,c.size()));
536  o2scl_linalg::solve_tridiag_sym<ubvector,ubvector,ubvector,
537  ubvector_range,o2scl_linalg::ubvector_4_mem,ubvector>
538  (diag,offdiag,g,cp1,sys_size,p4m);
539 
540  return;
541  }
542 
543  /// Give the value of the function \f$ y(x=x_0) \f$ .
544  virtual double eval(double x0) const {
545 
546  size_t cache=0;
547  size_t index=this->svx.find_const(x0,cache);
548 
549  double x_lo=(*this->px)[index];
550  double x_hi=(*this->px)[index+1];
551  double dx=x_hi-x_lo;
552 
553  double y_lo=(*this->py)[index];
554  double y_hi=(*this->py)[index+1];
555  double dy=y_hi-y_lo;
556  double delx=x0-x_lo;
557  double b_i, c_i, d_i;
558 
559  coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
560 
561  return y_lo+delx*(b_i+delx*(c_i+delx*d_i));
562  }
563 
564  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
565  virtual double deriv(double x0) const {
566 
567  size_t cache=0;
568  size_t index=this->svx.find_const(x0,cache);
569 
570  double x_lo=(*this->px)[index];
571  double x_hi=(*this->px)[index+1];
572  double dx=x_hi-x_lo;
573 
574  double y_lo=(*this->py)[index];
575  double y_hi=(*this->py)[index+1];
576  double dy=y_hi-y_lo;
577  double delx=x0-x_lo;
578  double b_i, c_i, d_i;
579 
580  coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
581 
582  return b_i+delx*(2.0*c_i+3.0*d_i*delx);
583  }
584 
585  /** \brief Give the value of the second derivative
586  \f$ y^{\prime \prime}(x=x_0) \f$ .
587  */
588  virtual double deriv2(double x0) const {
589 
590  size_t cache=0;
591  size_t index=this->svx.find_const(x0,cache);
592 
593  double x_lo=(*this->px)[index];
594  double x_hi=(*this->px)[index+1];
595  double dx=x_hi-x_lo;
596 
597  double y_lo=(*this->py)[index];
598  double y_hi=(*this->py)[index+1];
599  double dy=y_hi-y_lo;
600  double delx=x0-x_lo;
601  double b_i, c_i, d_i;
602 
603  coeff_calc(c,dy,dx,index,b_i,c_i,d_i);
604 
605  return 2.0*c_i+6.0*d_i*delx;
606  }
607 
608  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
609  virtual double integ(double a, double b) const {
610 
611  size_t i, index_a, index_b;
612 
613  bool flip=false;
614  if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
615  ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
616  double tmp=a;
617  a=b;
618  b=tmp;
619  flip=true;
620  }
621 
622  size_t cache=0;
623  index_a=this->svx.find_const(a,cache);
624  index_b=this->svx.find_const(b,cache);
625 
626  double result=0.0;
627 
628  for(i=index_a; i<=index_b; i++) {
629 
630  double x_lo=(*this->px)[i];
631  double x_hi=(*this->px)[i+1];
632  double y_lo=(*this->py)[i];
633  double y_hi=(*this->py)[i+1];
634  double dx=x_hi-x_lo;
635  double dy=y_hi-y_lo;
636 
637  if(dx != 0.0) {
638  double b_i, c_i, d_i;
639  coeff_calc(c,dy,dx,i,b_i,c_i,d_i);
640  if (i == index_a || i == index_b) {
641  double x1=(i == index_a) ? a : x_lo;
642  double x2=(i == index_b) ? b : x_hi;
643  result += this->integ_eval(y_lo,b_i,c_i,d_i,x_lo,x1,x2);
644  } else {
645  result += dx*(y_lo+dx*(0.5*b_i+
646  dx*(c_i/3.0+0.25*d_i*dx)));
647  }
648  } else {
649  std::string str=((std::string)"Interval of length zero ")+
650  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
651  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
652  " in interp_cspline::integ().";
653  O2SCL_ERR(str.c_str(),exc_einval);
654  }
655 
656  }
657 
658  if (flip) result*=-1.0;
659 
660  return result;
661  }
662 
663  /// Return the type, \c "interp_cspline".
664  virtual const char *type() const { return "interp_cspline"; }
665 
666 #ifndef DOXYGEN_INTERNAL
667 
668  private:
669 
672  (const interp_cspline<vec_t,vec2_t>&);
673 
674 #endif
675 
676  };
677 
678  /** \brief Cubic spline interpolation with periodic
679  boundary conditions (GSL)
680 
681  See also the \ref intp_section section of the \o2 User's guide.
682  */
683  template<class vec_t, class vec2_t=vec_t> class interp_cspline_peri :
684  public interp_cspline<vec_t,vec2_t> {
685 
686 #ifdef O2SCL_NEVER_DEFINED
687  }{
688 #endif
689 
690  protected:
691 
692  /// Memory for the tridiagonalization
694 
695  public:
696 
698  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
699  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
700  typedef boost::numeric::ublas::slice slice;
701  typedef boost::numeric::ublas::range range;
702 
703  interp_cspline_peri() : interp_cspline<vec_t,vec2_t>() {
704  this->min_size=2;
705  }
706 
707  virtual ~interp_cspline_peri() {
708  }
709 
710  /// Return the type, \c "interp_cspline_peri".
711  virtual const char *type() const { return "interp_cspline_peri"; }
712 
713  /** \brief Initialize interpolation routine
714  */
715  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
716 
717  if (size<this->min_size) {
718  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
719  " than "+szttos(this->min_size)+" in interp_cspline"+
720  "_peri::set().").c_str(),exc_einval);
721  }
722 
723  if (size!=this->sz) {
724  this->c.resize(size);
725  this->g.resize(size);
726  this->diag.resize(size);
727  this->offdiag.resize(size);
728  p5m.resize(size);
729  }
730 
731  this->px=&xa;
732  this->py=&ya;
733  this->sz=size;
734 
735  this->svx.set_vec(size,xa);
736 
737  /// Periodic boundary conditions
738 
739  size_t i;
740  size_t num_points=size;
741  // Engeln-Mullges+Uhlig "n"
742  size_t max_index=num_points-1;
743  // linear system is sys_size x sys_size
744  size_t sys_size=max_index;
745 
746  if (sys_size == 2) {
747 
748  // solve 2x2 system
749 
750  double h0=xa[1]-xa[0];
751  double h1=xa[2]-xa[1];
752  double h2=xa[3]-xa[2];
753  double A=2.0*(h0+h1);
754  double B=h0+h1;
755  double gx[2];
756  double det;
757 
758  gx[0]=3.0*((ya[2]-ya[1])/h1-(ya[1]-ya[0])/h0);
759  gx[1]=3.0*((ya[1]-ya[2])/h2-(ya[2]-ya[1])/h1);
760 
761  det=3.0*(h0+h1)*(h0+h1);
762  this->c[1]=( A*gx[0]-B*gx[1])/det;
763  this->c[2]=(-B*gx[0]+A*gx[1])/det;
764  this->c[0]=this->c[2];
765 
766  return;
767 
768  } else {
769 
770  for (i=0; i < sys_size-1; i++) {
771  double h_i=xa[i+1]-xa[i];
772  double h_ip1=xa[i+2]-xa[i+1];
773  double ydiff_i=ya[i+1]-ya[i];
774  double ydiff_ip1=ya[i+2]-ya[i+1];
775  double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
776  double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
777  this->offdiag[i]=h_ip1;
778  this->diag[i]=2.0*(h_ip1+h_i);
779  this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
780  }
781 
782  i=sys_size-1;
783  {
784  double h_i=xa[i+1]-xa[i];
785  double h_ip1=xa[1]-xa[0];
786  double ydiff_i=ya[i+1]-ya[i];
787  double ydiff_ip1=ya[1]-ya[0];
788  double g_i=(h_i != 0.0) ? 1.0/h_i : 0.0;
789  double g_ip1=(h_ip1 != 0.0) ? 1.0/h_ip1 : 0.0;
790  this->offdiag[i]=h_ip1;
791  this->diag[i]=2.0*(h_ip1+h_i);
792  this->g[i]=3.0*(ydiff_ip1*g_ip1-ydiff_i*g_i);
793  }
794 
795  ubvector_range cp1(this->c,range(1,this->c.size()));
796  o2scl_linalg::solve_cyc_tridiag_sym<ubvector,ubvector,ubvector,
797  ubvector_range,o2scl_linalg::ubvector_5_mem,ubvector>
798  (this->diag,this->offdiag,this->g,cp1,sys_size,p5m);
799  this->c[0]=this->c[max_index];
800 
801  return;
802  }
803 
804  }
805 
806 #ifndef DOXYGEN_INTERNAL
807 
808  private:
809 
813  (const interp_cspline_peri<vec_t,vec2_t>&);
814 
815 #endif
816 
817  };
818 
819  /** \brief Akima spline interpolation (GSL)
820 
821  See also the \ref intp_section section of the \o2 User's guide.
822 
823  This class uses natural boundary conditions, where the second
824  derivative vanishes at each end point. Extrapolation effectively
825  assumes that the second derivative is linear outside of the
826  endpoints. Use \ref interp_akima_peri for perodic boundary
827  conditions.
828  */
829  template<class vec_t, class vec2_t=vec_t> class interp_akima :
830  public interp_base<vec_t,vec2_t> {
831 
832  public:
833 
835  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
836  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
837  typedef boost::numeric::ublas::slice slice;
838  typedef boost::numeric::ublas::range range;
839 
840 #ifndef DOXYGEN_INTERNAL
841 
842  protected:
843 
844  /// \name Storage for Akima spline interpolation
845  //@{
846  ubvector b;
847  ubvector c;
848  ubvector d;
849  ubvector um;
850  //@}
851 
852  /// For initializing the interpolation
853  void akima_calc(const vec_t &x_array, size_t size,
854  ubvector &umx) {
855 
856  for(size_t i=0;i<this->sz-1;i++) {
857 
858  double NE=fabs(umx[3+i]-umx[2+i])+fabs(umx[1+i]-umx[i]);
859 
860  if (NE == 0.0) {
861  b[i]=umx[2+i];
862  c[i]=0.0;
863  d[i]=0.0;
864  } else {
865  double h_i=(*this->px)[i+1]-(*this->px)[i];
866  double NE_next=fabs(umx[4+i]-umx[3+i])+
867  fabs(umx[2+i]-umx[1+i]);
868  double alpha_i=fabs(umx[1+i]-umx[i])/NE;
869  double alpha_ip1;
870  double tL_ip1;
871  if (NE_next == 0.0) {
872  tL_ip1=umx[2+i];
873  } else {
874  alpha_ip1=fabs(umx[2+i]-umx[1+i])/NE_next;
875  tL_ip1=(1.0-alpha_ip1)*umx[2+i]+alpha_ip1*umx[3+i];
876  }
877  b[i]=(1.0-alpha_i)*umx[1+i]+alpha_i*umx[2+i];
878  c[i]=(3.0*umx[2+i]-2.0*b[i]-tL_ip1)/h_i;
879  d[i]=(b[i]+tL_ip1-2.0*umx[2+i])/(h_i*h_i);
880  }
881  }
882  }
883 
884 #endif
885 
886  public:
887 
888  /** \brief Create a base interpolation object with or without
889  periodic boundary conditions
890  */
892  this->min_size=5;
893  }
894 
895  virtual ~interp_akima() {
896  }
897 
898  /** \brief Initialize interpolation routine
899  */
900  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
901 
902  if (size<this->min_size) {
903  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
904  " than "+szttos(this->min_size)+" in interp_akima::"+
905  "set().").c_str(),exc_einval);
906  }
907 
908  if (size!=this->sz) {
909  b.resize(size);
910  c.resize(size);
911  d.resize(size);
912  um.resize(size+4);
913  }
914 
915  this->px=&xa;
916  this->py=&ya;
917  this->sz=size;
918 
919  this->svx.set_vec(size,xa);
920 
921  // Non-periodic boundary conditions
922 
923  ubvector_range m(um,range(2,um.size()));
924  size_t i;
925  for (i=0;i<=size-2;i++) {
926  m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
927  }
928 
929  um[0]=3.0*m[0]-2.0*m[1];
930  um[1]=2.0*m[0]-m[1];
931  m[this->sz-1]=2.0*m[size-2]-m[size-3];
932  m[size]=3.0*m[size-2]-2.0*m[size-3];
933 
934  akima_calc(xa,size,um);
935 
936  return;
937  }
938 
939  /// Give the value of the function \f$ y(x=x_0) \f$ .
940  virtual double eval(double x0) const {
941 
942  size_t cache=0;
943  size_t index=this->svx.find_const(x0,cache);
944 
945  double x_lo=(*this->px)[index];
946  double delx=x0-x_lo;
947  double bb=b[index];
948  double cc=c[index];
949  double dd=d[index];
950 
951  return (*this->py)[index]+delx*(bb+delx*(cc+dd*delx));
952  }
953 
954  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
955  virtual double deriv(double x0) const {
956 
957  size_t cache=0;
958  size_t index=this->svx.find_const(x0,cache);
959 
960  double x_lo=(*this->px)[index];
961  double delx=x0-x_lo;
962  double bb=b[index];
963  double cc=c[index];
964  double dd=d[index];
965 
966  return bb+delx*(2.0*cc+3.0*dd*delx);
967  }
968 
969  /** \brief Give the value of the second derivative
970  \f$ y^{\prime \prime}(x=x_0) \f$ .
971  */
972  virtual double deriv2(double x0) const {
973 
974  size_t cache=0;
975  size_t index=this->svx.find_const(x0,cache);
976 
977  double x_lo=(*this->px)[index];
978  double delx=x0-x_lo;
979  double cc=c[index];
980  double dd=d[index];
981 
982  return 2.0*cc+6.0*dd*delx;
983  }
984 
985  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
986  virtual double integ(double aa, double bb) const {
987 
988  size_t i, index_a, index_b;
989 
990  bool flip=false;
991  if (((*this->px)[0]<(*this->px)[this->sz-1] && aa>bb) ||
992  ((*this->px)[0]>(*this->px)[this->sz-1] && aa<bb)) {
993  double tmp=aa;
994  aa=bb;
995  bb=tmp;
996  flip=true;
997  }
998 
999  size_t cache=0;
1000  index_a=this->svx.find_const(aa,cache);
1001  index_b=this->svx.find_const(bb,cache);
1002 
1003  double result=0.0;
1004 
1005  for(i=index_a; i<=index_b; i++) {
1006 
1007  double x_lo=(*this->px)[i];
1008  double x_hi=(*this->px)[i+1];
1009  double y_lo=(*this->py)[i];
1010  double dx=x_hi-x_lo;
1011 
1012  if (dx != 0.0) {
1013 
1014  if (i==index_a || i==index_b) {
1015  double x1=(i==index_a) ? aa : x_lo;
1016  double x2=(i==index_b) ? bb : x_hi;
1017  result += this->integ_eval(y_lo,b[i],c[i],d[i],x_lo,x1,x2);
1018  } else {
1019  result+=dx*(y_lo+dx*(0.5*b[i]+dx*(c[i]/3.0+0.25*d[i]*dx)));
1020  }
1021 
1022  } else {
1023  double y_hi=(*this->py)[i+1];
1024  std::string str=((std::string)"Interval of length zero ")+
1025  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
1026  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
1027  " in interp_akima::integ().";
1028  O2SCL_ERR(str.c_str(),exc_einval);
1029  }
1030  }
1031 
1032  if (flip) result*=-1.0;
1033 
1034  return result;
1035  }
1036 
1037  /// Return the type, \c "interp_akima".
1038  virtual const char *type() const { return "interp_akima"; }
1039 
1040 #ifndef DOXYGEN_INTERNAL
1041 
1042  private:
1043 
1046 
1047 #endif
1048 
1049  };
1050 
1051  /** \brief Akima spline interpolation with periodic
1052  boundary conditions (GSL)
1053 
1054  See also the \ref intp_section section of the \o2 User's guide.
1055  */
1056  template<class vec_t, class vec2_t=vec_t> class interp_akima_peri :
1057  public interp_akima<vec_t,vec2_t> {
1058 
1059  public:
1060 
1062  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
1063  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
1064  typedef boost::numeric::ublas::slice slice;
1065  typedef boost::numeric::ublas::range range;
1066 
1067  public:
1068 
1070  }
1071 
1072  virtual ~interp_akima_peri() {
1073  }
1074 
1075  /// Return the type, \c "interp_akima_peri".
1076  virtual const char *type() const { return "interp_akima_peri"; }
1077 
1078  /// Initialize interpolation routine
1079  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
1080 
1081  if (size<this->min_size) {
1082  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
1083  " than "+szttos(this->min_size)+" in interp_akima"+
1084  "_peri::set().").c_str(),exc_einval);
1085  }
1086 
1087  if (size!=this->sz) {
1088  this->b.resize(size);
1089  this->c.resize(size);
1090  this->d.resize(size);
1091  this->um.resize(size+4);
1092  }
1093 
1094  this->px=&xa;
1095  this->py=&ya;
1096  this->sz=size;
1097 
1098  this->svx.set_vec(size,xa);
1099 
1100  // Periodic boundary conditions
1101 
1102  ubvector_range m(this->um,range(2,this->um.size()));
1103 
1104  // Form the required set of divided differences
1105  for (size_t i=0;i<=size-2;i++) {
1106  m[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]);
1107  }
1108 
1109  this->um[0]=m[this->sz-3];
1110  this->um[1]=m[this->sz-2];
1111  m[this->sz-1]=m[0];
1112  m[this->sz]=m[1];
1113 
1114  this->akima_calc(xa,size,this->um);
1115 
1116  return;
1117  }
1118 
1119 #ifndef DOXYGEN_INTERNAL
1120 
1121  private:
1122 
1125  (const interp_akima_peri<vec_t,vec2_t>&);
1126 
1127 #endif
1128 
1129  };
1130 
1131  /** \brief Steffen's monotonicity-preserving interpolation
1132 
1133  Adapted from the GSL version by J.-F. Caron which
1134  was based on \ref Steffen90 .
1135  */
1136  template<class vec_t, class vec2_t=vec_t> class interp_steffen :
1137  public interp_base<vec_t,vec2_t> {
1138 
1139 #ifdef O2SCL_NEVER_DEFINED
1140  }{
1141 #endif
1142 
1143  public:
1144 
1146  typedef boost::numeric::ublas::vector_slice<ubvector> ubvector_slice;
1147  typedef boost::numeric::ublas::vector_range<ubvector> ubvector_range;
1148  typedef boost::numeric::ublas::slice slice;
1149  typedef boost::numeric::ublas::range range;
1150 
1151 #ifndef DOXYGEN_INTERNAL
1152 
1153  protected:
1154 
1155  /// \name Storage for cubic spline interpolation
1156  //@{
1157  ubvector a;
1158  ubvector b;
1159  ubvector c;
1160  ubvector d;
1161  ubvector y_prime;
1162  //@}
1163 
1164  /** \brief Flip the sign of x if x and y have different signs
1165  */
1166  double copysign(const double x, const double y) {
1167  if ((x < 0 && y > 0) || (x > 0 && y < 0)) {
1168  return -x;
1169  }
1170  return x;
1171  }
1172 
1173 #endif
1174 
1175  public:
1176 
1177  /** \brief Create a base interpolation object */
1179  this->min_size=3;
1180  }
1181 
1182  virtual ~interp_steffen() {
1183  }
1184 
1185  /** \brief Initialize interpolation routine
1186  */
1187  virtual void set(size_t size, const vec_t &xa, const vec2_t &ya) {
1188 
1189  if (size<this->min_size) {
1190  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
1191  " than "+szttos(this->min_size)+" in interp_steffen::"+
1192  "set().").c_str(),exc_einval);
1193  }
1194 
1195  if (size!=this->sz) {
1196  a.resize(size);
1197  b.resize(size);
1198  c.resize(size);
1199  d.resize(size);
1200  y_prime.resize(size);
1201  }
1202 
1203  this->px=&xa;
1204  this->py=&ya;
1205  this->sz=size;
1206 
1207  this->svx.set_vec(size,xa);
1208 
1209  /*
1210  * first assign the interval and slopes for the left boundary.
1211  * We use the "simplest possibility" method described in the paper
1212  * in section 2.2
1213  */
1214  double h0=(xa[1]-xa[0]);
1215  double s0=(ya[1]-ya[0]) / h0;
1216 
1217  y_prime[0]=s0;
1218 
1219  /* Now we calculate all the necessary s, h, p, and y' variables
1220  from 1 to N-2 (0 to size-2 inclusive) */
1221  for (size_t i=1; i < (size-1); i++) {
1222 
1223  double pi;
1224 
1225  /* equation 6 in the paper */
1226  double hi=(xa[i+1]-xa[i]);
1227  double him1=(xa[i]-xa[i-1]);
1228 
1229  /* equation 7 in the paper */
1230  double si=(ya[i+1]-ya[i]) / hi;
1231  double sim1=(ya[i]-ya[i-1]) / him1;
1232 
1233  /* equation 8 in the paper */
1234  pi=(sim1*hi + si*him1) / (him1 + hi);
1235 
1236  /* This is a C equivalent of the FORTRAN statement below eqn 11 */
1237  y_prime[i]=(copysign(1.0,sim1)+copysign(1.0,si))*
1238  std::min(fabs(sim1),std::min(fabs(si),0.5*fabs(pi)));
1239 
1240  }
1241 
1242  /*
1243  * we also need y' for the rightmost boundary; we use the
1244  * "simplest possibility" method described in the paper in
1245  * section 2.2
1246  *
1247  * y'=s_{n-1}
1248  */
1249  y_prime[size-1]=(ya[size-1]-ya[size-2])/
1250  (xa[size-1]-xa[size-2]);
1251 
1252  /* Now we can calculate all the coefficients for the whole range. */
1253  for (size_t i=0; i < (size-1); i++) {
1254  double hi=(xa[i+1]-xa[i]);
1255  double si=(ya[i+1]-ya[i]) / hi;
1256 
1257  /* These are from equations 2-5 in the paper. */
1258  a[i]=(y_prime[i] + y_prime[i+1]-2*si) / hi / hi;
1259  b[i]=(3*si-2*y_prime[i]-y_prime[i+1]) / hi;
1260  c[i]=y_prime[i];
1261  d[i]=ya[i];
1262  }
1263 
1264  return;
1265  }
1266 
1267  /// Give the value of the function \f$ y(x=x_0) \f$ .
1268  virtual double eval(double x0) const {
1269 
1270  size_t cache=0;
1271  size_t index=this->svx.find_const(x0,cache);
1272  double x_lo=(*this->px)[index];
1273  double delx=x0-x_lo;
1274 
1275  /* Use Horner's scheme for efficient evaluation of polynomials. */
1276  double y = d[index]+delx*(c[index]+delx*(b[index]+delx*a[index]));
1277 
1278  return y;
1279  }
1280 
1281  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1282  virtual double deriv(double x0) const {
1283 
1284  size_t cache=0;
1285  size_t index=this->svx.find_const(x0,cache);
1286  double x_lo=(*this->px)[index];
1287  double delx=x0-x_lo;
1288 
1289  return c[index]+delx*(2.0*b[index]+delx*3.0*a[index]);
1290  }
1291 
1292  /** \brief Give the value of the second derivative
1293  \f$ y^{\prime \prime}(x=x_0) \f$ .
1294  */
1295  virtual double deriv2(double x0) const {
1296 
1297  size_t cache=0;
1298  size_t index=this->svx.find_const(x0,cache);
1299  double x_lo=(*this->px)[index];
1300  double delx=x0-x_lo;
1301 
1302  return 2.0*b[index]+delx*6.0*a[index];
1303  }
1304 
1305  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1306  virtual double integ(double al, double bl) const {
1307 
1308  size_t i, index_a, index_b;
1309 
1310  bool flip=false;
1311  if (((*this->px)[0]<(*this->px)[this->sz-1] && al>bl) ||
1312  ((*this->px)[0]>(*this->px)[this->sz-1] && al<bl)) {
1313  double tmp=al;
1314  al=bl;
1315  bl=tmp;
1316  flip=true;
1317  }
1318 
1319  size_t cache=0;
1320  index_a=this->svx.find_const(al,cache);
1321  index_b=this->svx.find_const(bl,cache);
1322 
1323  double result=0.0;
1324 
1325  for(i=index_a; i<=index_b; i++) {
1326 
1327  double x_lo=(*this->px)[i];
1328  double x_hi=(*this->px)[i+1];
1329  double y_lo=(*this->py)[i];
1330  double y_hi=(*this->py)[i+1];
1331  double dx=x_hi-x_lo;
1332 
1333  if(dx != 0.0) {
1334 
1335  double x1=(i == index_a) ? al-x_lo : 0.0;
1336  double x2=(i == index_b) ? bl-x_lo : x_hi-x_lo;
1337  result += (1.0/4.0)*a[i]*(x2*x2*x2*x2-x1*x1*x1*x1)+
1338  (1.0/3.0)*b[i]*(x2*x2*x2-x1*x1*x1)+
1339  (1.0/2.0)*c[i]*(x2*x2-x1*x1)+d[i]*(x2-x1);
1340 
1341  } else {
1342  std::string str=((std::string)"Interval of length zero ")+
1343  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
1344  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
1345  " in interp_steffen::integ().";
1346  O2SCL_ERR(str.c_str(),exc_einval);
1347  }
1348 
1349  }
1350 
1351  if (flip) result*=-1.0;
1352 
1353  return result;
1354  }
1355 
1356  /// Return the type, \c "interp_steffen".
1357  virtual const char *type() const { return "interp_steffen"; }
1358 
1359 #ifndef DOXYGEN_INTERNAL
1360 
1361  private:
1362 
1364  interp_steffen<vec_t,vec2_t>& operator=
1365  (const interp_steffen<vec_t,vec2_t>&);
1366 
1367 #endif
1368 
1369  };
1370 
1371  /** \brief Monotonicity-preserving interpolation
1372 
1373  \warning This class is experimental. Integrals don't work yet.
1374 
1375  This class uses a method based on cubic Hermite interpolation,
1376  modifying the slopes to guarantee monotonicity. In the
1377  notation of \ref Fritsch80, if
1378  \f[
1379  \alpha_i^2+\beta_i^2 \geq 9
1380  \f]
1381  then \f$ \alpha \f$ and \f$ \beta \f$ are decreased by
1382  the factor \f$ \tau \f$ as described at the end of
1383  section 4.
1384 
1385  \note The results of the interpolation will only be monotonic in
1386  the regions where the original data set is also monotonic. Also,
1387  this interpolation routine does not enforce strict monotonicity,
1388  and the results of the interpolation will be flat where the data
1389  is also flat.
1390 
1391  Based on \ref Fritsch80 .
1392 
1393  \future Convert into fewer loops over the data
1394  */
1395  template<class vec_t, class vec2_t=vec_t> class interp_monotonic :
1396  public interp_base<vec_t,vec2_t> {
1397 
1398 #ifdef O2SCL_NEVER_DEFINED
1399  }{
1400 #endif
1401 
1402  public:
1403 
1405 
1406  protected:
1407 
1408  /// Slopes
1409  ubvector m;
1410  /// Finite differences
1411  ubvector Delta;
1412  /// Ratio
1413  ubvector alpha;
1414  /// Staggered ratio
1415  ubvector beta;
1416 
1417  public:
1418 
1419  interp_monotonic() {
1420  this->min_size=2;
1421  }
1422 
1423  virtual ~interp_monotonic() {
1424  }
1425 
1426  /// Initialize interpolation routine
1427  virtual void set(size_t size, const vec_t &x, const vec2_t &y) {
1428 
1429  // Verify size
1430  if (size<this->min_size) {
1431  O2SCL_ERR((((std::string)"Vector size, ")+szttos(size)+", is less"+
1432  " than "+szttos(this->min_size)+" in interp_monotonic::"+
1433  "set().").c_str(),exc_einval);
1434  }
1435 
1436  // Setup search_vec object
1437  this->svx.set_vec(size,x);
1438 
1439  // Resize internal vectors
1440  if (this->sz!=size) {
1441  m.resize(size);
1442  Delta.resize(size-1);
1443  alpha.resize(size-1);
1444  beta.resize(size-1);
1445  }
1446 
1447  // Copy pointers
1448  this->px=&x;
1449  this->py=&y;
1450  this->sz=size;
1451 
1452  // Create the interpolation arrays in this sub-interval
1453 
1454  // Compute Delta and m
1455  for(size_t i=0;i<size-1;i++) {
1456  Delta[i]=(y[i+1]-y[i])/(x[i+1]-x[i]);
1457  if (i>0) {
1458  m[i]=(Delta[i]+Delta[i-1])/2.0;
1459  }
1460  }
1461  m[0]=Delta[0];
1462  m[size-1]=Delta[size-2];
1463 
1464  // Check to see if the data is flat anywhere
1465  for(size_t i=0;i<size-1;i++) {
1466  if (y[i]==y[i+1]) {
1467  m[i]=0.0;
1468  m[i+1]=0.0;
1469  }
1470  }
1471 
1472  // Compute alpha and beta
1473  for(size_t i=0;i<size-1;i++) {
1474  alpha[i]=m[i]/Delta[i];
1475  beta[i]=m[i+1]/Delta[i];
1476  }
1477 
1478  // Constrain m to ensure monotonicity
1479  for(size_t i=0;i<size-1;i++) {
1480  double norm2=alpha[i]*alpha[i]+beta[i]*beta[i];
1481  if (norm2>9.0) {
1482  double tau=3.0/sqrt(norm2);
1483  m[i]=tau*alpha[i]*Delta[i];
1484  m[i+1]=tau*beta[i]*Delta[i];
1485  }
1486  }
1487 
1488  return;
1489  }
1490 
1491  /// Give the value of the function \f$ y(x=x_0) \f$ .
1492  virtual double eval(double x0) const {
1493 
1494  size_t cache=0;
1495  size_t index=this->svx.find_const(x0,cache);
1496 
1497  double x_lo=(*this->px)[index];
1498  double x_hi=(*this->px)[index+1];
1499  double y_lo=(*this->py)[index];
1500  double y_hi=(*this->py)[index+1];
1501  double h=x_hi-x_lo;
1502  double t=(x0-x_lo)/h;
1503  double t2=t*t, t3=t2*t;
1504 
1505  double h00=2.0*t3-3.0*t2+1.0;
1506  double h10=t3-2.0*t2+t;
1507  double h01=-2.0*t3+3.0*t2;
1508  double h11=t3-t2;
1509  double interp=y_lo*h00+h*m[index]*h10+y_hi*h01+h*m[index+1]*h11;
1510 
1511  return interp;
1512  }
1513 
1514  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1515  virtual double deriv(double x0) const {
1516 
1517  size_t cache=0;
1518  size_t index=this->svx.find_const(x0,cache);
1519 
1520  double x_lo=(*this->px)[index];
1521  double x_hi=(*this->px)[index+1];
1522  double y_lo=(*this->py)[index];
1523  double y_hi=(*this->py)[index+1];
1524  double h=x_hi-x_lo;
1525  double t=(x0-x_lo)/h;
1526  double t2=t*t;
1527 
1528  double dh00=6.0*t2-6.0*t;
1529  double dh10=3.0*t2-4.0*t+1.0;
1530  double dh01=-6.0*t2+6.0*t;
1531  double dh11=3.0*t2-2.0*t;
1532  double deriv=(y_lo*dh00+h*m[index]*dh10+y_hi*dh01+
1533  h*m[index+1]*dh11)/h;
1534 
1535  return deriv;
1536  }
1537 
1538  /** \brief Give the value of the second derivative
1539  \f$ y^{\prime \prime}(x=x_0) \f$
1540  */
1541  virtual double deriv2(double x0) const {
1542 
1543  size_t cache=0;
1544  size_t index=this->svx.find_const(x0,cache);
1545 
1546  double x_lo=(*this->px)[index];
1547  double x_hi=(*this->px)[index+1];
1548  double y_lo=(*this->py)[index];
1549  double y_hi=(*this->py)[index+1];
1550  double h=x_hi-x_lo;
1551  double t=(x0-x_lo)/h;
1552 
1553  double ddh00=12.0*t-6.0;
1554  double ddh10=6.0*t-4.0;
1555  double ddh01=-12.0*t+6.0;
1556  double ddh11=6.0*t-2.0;
1557  double deriv2=(y_lo*ddh00+h*m[index]*ddh10+y_hi*ddh01+
1558  h*m[index+1]*ddh11)/h/h;
1559 
1560  return deriv2;
1561  }
1562 
1563  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1564  virtual double integ(double a, double b) const {
1565 
1566  size_t i, index_a, index_b;
1567 
1568  bool flip=false;
1569  if (((*this->px)[0]<(*this->px)[this->sz-1] && a>b) ||
1570  ((*this->px)[0]>(*this->px)[this->sz-1] && a<b)) {
1571  double tmp=a;
1572  a=b;
1573  b=tmp;
1574  flip=true;
1575  }
1576 
1577  size_t cache=0;
1578  index_a=this->svx.find_const(a,cache);
1579  index_b=this->svx.find_const(b,cache);
1580 
1581  double result=0.0;
1582 
1583  for(i=index_a; i<=index_b; i++) {
1584 
1585  double x_hi=(*this->px)[i+1];
1586  double x_lo=(*this->px)[i];
1587  double y_lo=(*this->py)[i];
1588  double y_hi=(*this->py)[i+1];
1589  double h=x_hi-x_lo;
1590 
1591  if (h != 0.0) {
1592 
1593  if (i == index_a) {
1594  x_lo=a;
1595  }
1596  if (i == index_b) {
1597  x_hi=b;
1598  }
1599 
1600  double t=(x_hi-x_lo)/h;
1601  double t2=t*t, t3=t2*t, t4=t3*t;
1602 
1603  double ih00=t4/2.0-t3+t;
1604  double ih10=t4/4.0-2.0*t3/3.0+t2/2.0;
1605  double ih01=-t4/2.0+t3;
1606  double ih11=t4/4.0-t3/3.0;
1607  double intres=h*(y_lo*ih00+h*m[i]*ih10+y_hi*ih01+
1608  h*m[i+1]*ih11);
1609  result+=intres;
1610 
1611  } else {
1612  std::string str=((std::string)"Interval of length zero ")+
1613  "between ("+o2scl::dtos(x_lo)+","+o2scl::dtos(y_lo)+
1614  ") and ("+o2scl::dtos(x_hi)+","+o2scl::dtos(y_hi)+
1615  " in interp_monotonic::integ().";
1616  O2SCL_ERR(str.c_str(),exc_einval);
1617  }
1618 
1619  }
1620 
1621  if (flip) result*=-1.0;
1622 
1623  return result;
1624  }
1625 
1626  /// Return the type, \c "interp_monotonic".
1627  virtual const char *type() const { return "interp_monotonic"; }
1628 
1629 #ifndef DOXYGEN_INTERNAL
1630 
1631  private:
1632 
1635  (const interp_monotonic<vec_t,vec2_t>&);
1636 
1637 #endif
1638 
1639  };
1640 
1641  /** \brief Interpolation class for general vectors
1642 
1643  See also the \ref intp_section section of the \o2 User's guide.
1644 
1645  Interpolation of ublas vector like objects is performed with the
1646  default template parameters, and \ref interp_array is provided
1647  for simple interpolation on C-style arrays.
1648 
1649  The type of interpolation to be performed can be specified using
1650  the set_type() function or in the constructor. The default is
1651  cubic splines with natural boundary conditions.
1652  */
1653  template<class vec_t=boost::numeric::ublas::vector<double>,
1654  class vec2_t=vec_t> class interp {
1655 
1656 #ifndef DOXYGEN_INTERNAL
1657 
1658  protected:
1659 
1660  /// Pointer to base interpolation object
1662 
1663 #endif
1664 
1665  public:
1666 
1667  /// Create with base interpolation object \c it
1668  interp(size_t interp_type=itp_cspline) {
1669  if (interp_type==itp_linear) {
1671  } else if (interp_type==itp_cspline) {
1673  } else if (interp_type==itp_cspline_peri) {
1675  } else if (interp_type==itp_akima) {
1677  } else if (interp_type==itp_akima_peri) {
1679  } else if (interp_type==itp_monotonic) {
1681  } else if (interp_type==itp_steffen) {
1683  } else if (interp_type==itp_nearest_neigh) {
1685  } else {
1686  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1687  o2scl::szttos(interp_type)+", in "+
1688  "interp::interp().").c_str(),exc_einval);
1689  }
1690  }
1691 
1692  virtual ~interp() {
1693  delete itp;
1694  }
1695 
1696  /** \brief Give the value of the function, \f$ y(x=x_0) \f$ , as
1697  specified as the first \c n elements of vectors \c x and \c y
1698  */
1699  virtual double eval(const double x0, size_t n, const vec_t &x,
1700  const vec2_t &y) {
1701  itp->set(n,x,y);
1702  return itp->eval(x0);
1703  }
1704 
1705  /** \brief Give the value of the derivative, \f$ y^{\prime}(x=x_0)
1706  \f$ , where \f$ y(x) \f$ is specified in the first
1707  \c n elements of vectors \c x and
1708  \c y
1709  */
1710  virtual double deriv(const double x0, size_t n, const vec_t &x,
1711  const vec2_t &y) {
1712  itp->set(n,x,y);
1713  return itp->deriv(x0);
1714  }
1715 
1716  /** \brief Give the value of the second derivative, \f$
1717  y^{\prime\prime}(x=x_0) \f$ , where \f$ y(x) \f$ is specified in
1718  the first \c n elements of vectors \c x and \c y
1719  */
1720  virtual double deriv2(const double x0, size_t n, const vec_t &x,
1721  const vec2_t &y) {
1722  itp->set(n,x,y);
1723  return itp->deriv2(x0);
1724  }
1725 
1726  /** \brief Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$
1727  , where \f$ y(x) \f$ is specified in the first \c n elements of
1728  vectors \c x and \c y
1729  */
1730  virtual double integ(const double x1, const double x2, size_t n,
1731  const vec_t &x, const vec2_t &y) {
1732  itp->set(n,x,y);
1733  return itp->integ(x1,x2);
1734  }
1735 
1736  /// Set base interpolation type
1737  void set_type(size_t interp_type) {
1738  delete itp;
1739  if (interp_type==itp_linear) {
1741  } else if (interp_type==itp_cspline) {
1743  } else if (interp_type==itp_cspline_peri) {
1745  } else if (interp_type==itp_akima) {
1747  } else if (interp_type==itp_akima_peri) {
1749  } else if (interp_type==itp_monotonic) {
1751  } else if (interp_type==itp_steffen) {
1753  } else if (interp_type==itp_nearest_neigh) {
1755  } else {
1756  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1757  o2scl::szttos(interp_type)+", in "+
1758  "interp::set().").c_str(),exc_einval);
1759  }
1760  return;
1761  }
1762 
1763 #ifndef DOXYGEN_INTERNAL
1764 
1765  private:
1766 
1768  interp<vec_t,vec2_t>& operator=(const interp<vec_t,vec2_t>&);
1769 
1770 #endif
1771 
1772  };
1773 
1774  /** \brief Interpolation class for pre-specified vector
1775 
1776  See also the \ref intp_section section of the \o2 User's guide.
1777 
1778  This interpolation class is intended to be sufficiently general
1779  to handle any vector type. Interpolation of ublas vector like
1780  objects is performed with the default template parameters.
1781 
1782  This class does not double check the vector to ensure that
1783  all of the intervals for the abcissa are all increasing or
1784  all decreasing.
1785 
1786  The type of interpolation to be performed can be specified
1787  using the set_type() function. The default is cubic splines
1788  with natural boundary conditions.
1789 
1790  \future Make version which copies vectors
1791  rather than storing pointers might be better and then
1792  has copy constructors.
1793  */
1794  template<class vec_t=boost::numeric::ublas::vector<double>,
1795  class vec2_t=vec_t> class interp_vec :
1796  public interp_base<vec_t,vec2_t> {
1797 
1798 #ifndef DOXYGEN_INTERNAL
1799 
1800  protected:
1801 
1802  /// Base interpolation object
1804 
1805  /// Interpolation type
1806  size_t itype;
1807 
1808 #endif
1809 
1810  public:
1811 
1812  /// Blank interpolator
1814  itp=0;
1815  itype=itp_cspline;
1816  }
1817 
1818  /** \brief Create an interpolation object with interpolation type
1819  \c itp_cspline based on the first \c n entries of vectors
1820  \c x and \c y
1821  */
1822  interp_vec(size_t n, const vec_t &x,
1823  const vec2_t &y, size_t interp_type=itp_cspline) {
1824 
1825  if (x[0]==x[n-1]) {
1826  O2SCL_ERR((((std::string)"Vector endpoints equal (value=")+
1827  o2scl::dtos(x[0])+") in interp_vec()::"+
1828  "interp_vec().").c_str(),exc_einval);
1829  }
1830 
1831  if (interp_type==itp_linear) {
1833  } else if (interp_type==itp_cspline) {
1835  } else if (interp_type==itp_cspline_peri) {
1837  } else if (interp_type==itp_akima) {
1839  } else if (interp_type==itp_akima_peri) {
1841  } else if (interp_type==itp_monotonic) {
1843  } else if (interp_type==itp_steffen) {
1845  } else if (interp_type==itp_nearest_neigh) {
1847  } else {
1848  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1849  o2scl::szttos(interp_type)+", in "+
1850  "interp_vec::interp_vec().").c_str(),exc_einval);
1851  }
1852  itype=interp_type;
1853 
1854  itp->set(n,x,y);
1855  }
1856 
1857  virtual ~interp_vec() {
1858  if (itp!=0) {
1859  delete itp;
1860  }
1861  }
1862 
1863  /** \brief Modify the interpolation object to operate on the first
1864  \c n entries of vectors \c x and \c y
1865  */
1866  void set(size_t n, const vec_t &x, const vec2_t &y) {
1867  set(n,x,y,itype);
1868  return;
1869  }
1870 
1871  /** \brief Set a new vector to interpolate
1872  */
1873  void set(size_t n, const vec_t &x,
1874  const vec2_t &y, size_t interp_type) {
1875 
1876  if (x[0]==x[n-1]) {
1877  O2SCL_ERR((((std::string)"Vector endpoints equal (value=")+
1878  o2scl::dtos(x[0])+") in interp_vec()::"+
1879  "interp_vec().").c_str(),exc_einval);
1880  }
1881 
1882  delete itp;
1883  if (interp_type==itp_linear) {
1885  } else if (interp_type==itp_cspline) {
1887  } else if (interp_type==itp_cspline_peri) {
1889  } else if (interp_type==itp_akima) {
1891  } else if (interp_type==itp_akima_peri) {
1893  } else if (interp_type==itp_monotonic) {
1895  } else if (interp_type==itp_steffen) {
1897  } else if (interp_type==itp_nearest_neigh) {
1899  } else {
1900  O2SCL_ERR((((std::string)"Invalid interpolation type, ")+
1901  o2scl::szttos(interp_type)+", in "+
1902  "interp_vec::set().").c_str(),exc_einval);
1903  }
1904  itype=interp_type;
1905 
1906  itp->set(n,x,y);
1907  }
1908 
1909  /** \brief Manually clear the pointer to the user-specified vector
1910  */
1911  void clear() {
1912  if (itp!=0) {
1913  delete itp;
1914  itp=0;
1915  }
1916  return;
1917  }
1918 
1919  /// Give the value of the function \f$ y(x=x_0) \f$ .
1920  virtual double eval(const double x0) const {
1921  if (itp==0) {
1922  O2SCL_ERR("No vector set in interp_vec::eval().",
1923  exc_einval);
1924  }
1925  return itp->eval(x0);
1926  }
1927 
1928  /// Give the value of the function \f$ y(x=x_0) \f$ .
1929  virtual double operator()(double x0) const {
1930  if (itp==0) {
1931  O2SCL_ERR("No vector set in interp_vec::operator().",
1932  exc_einval);
1933  }
1934  return itp->eval(x0);
1935  }
1936 
1937  /// Give the value of the derivative \f$ y^{\prime}(x=x_0) \f$ .
1938  virtual double deriv(const double x0) const {
1939  if (itp==0) {
1940  O2SCL_ERR("No vector set in interp_vec::deriv().",
1941  exc_einval);
1942  }
1943  return itp->deriv(x0);
1944  }
1945 
1946  /** \brief Give the value of the second derivative
1947  \f$ y^{\prime \prime}(x=x_0) \f$ .
1948  */
1949  virtual double deriv2(const double x0) const {
1950  if (itp==0) {
1951  O2SCL_ERR("No vector set in interp_vec::deriv2().",
1952  exc_einval);
1953  }
1954  return itp->deriv2(x0);
1955  }
1956 
1957  /// Give the value of the integral \f$ \int_a^{b}y(x)~dx \f$ .
1958  virtual double integ(const double x1, const double x2) const {
1959  if (itp==0) {
1960  O2SCL_ERR("No vector set in interp_vec::integ().",
1961  exc_einval);
1962  }
1963  return itp->integ(x1,x2);
1964  }
1965 
1966  /// Return the type, "interp_vec"
1967  virtual const char *type() const {
1968  return "interp_vec";
1969  }
1970 
1971 #ifndef DOXYGEN_INTERNAL
1972 
1973  private:
1974 
1976  interp_vec<vec_t,vec2_t> &operator=
1977  (const interp_vec<vec_t,vec2_t> &it);
1978 
1979 #endif
1980 
1981  };
1982 
1983  /** \brief A specialization of interp for C-style double arrays
1984 
1985  See also the \ref intp_section section of the \o2 User's guide.
1986  */
1987  template<size_t n> class interp_array :
1988  public interp<double[n]> {
1989 
1990  public:
1991 
1992  /// Create with base interpolation objects \c it and \c rit
1993  interp_array(size_t interp_type)
1994  : interp<double[n]>(interp_type) {}
1995 
1996  /** \brief Create an interpolator using the default base
1997  interpolation objects
1998  */
1999  interp_array() : interp<double[n]>() {}
2000 
2001  };
2002 
2003  /** \brief A specialization of \ref o2scl::interp_vec for C-style arrays
2004 
2005  See also the \ref intp_section section of the \o2 User's guide.
2006  */
2007  template<class arr_t> class interp_array_vec :
2008  public interp_vec<arr_t> {
2009 
2010  public:
2011 
2012  /// Create with base interpolation object \c it
2013  interp_array_vec(size_t nv, const arr_t &x, const arr_t &y,
2014  size_t interp_type) :
2015  interp_vec<arr_t>(nv,x,y,interp_type) {}
2016  };
2017 
2018  /** \brief Count level crossings
2019 
2020  This function counts the number of times the function \f$ y(x) =
2021  \mathrm{level} \f$ where the function is defined by the vectors
2022  \c x and \c y.
2023 
2024  The value returned is exactly the same as the size of the
2025  \c locs vector computed by \ref vector_find_level().
2026 
2027  This function will call the error handler if \c n is
2028  less than two.
2029 
2030  If one of the entries in \c y is either larger or smaller than
2031  its neighbors (i.e. if the function \f$ y(x) \f$ has an
2032  extremum), and if the value of \c level is exactly equal to the
2033  extremum, then this is counted as 1 level crossing. The same
2034  applies if either the first or last entry in \c y is exactly
2035  equal to \c level . Finite-precision errors may affect
2036  the result of this function when \c level is nearly
2037  equal to one of the value in the vector \c y.
2038  */
2039  template<class vec_t, class vec2_t> size_t vector_level_count
2040  (double level, size_t n, vec_t &x, vec2_t &y) {
2041 
2042  if (n<=1) {
2043  O2SCL_ERR2("Need at least two data points in ",
2044  "vector_find_count().",exc_einval);
2045  }
2046 
2047  size_t count=0;
2048 
2049  // Handle left boundary
2050  if (y[0]==level) count++;
2051 
2052  // Find intersections
2053  for(size_t i=0;i<n-1;i++) {
2054 
2055  if ((y[i]<level && y[i+1]>=level) ||
2056  (y[i]>level && y[i+1]<=level)) {
2057  count++;
2058  }
2059  }
2060 
2061  return count;
2062  }
2063 
2064  /** \brief Perform inverse linear interpolation
2065 
2066  This function performs inverse linear interpolation of the data
2067  defined by \c x and \c y, finding all points in \c x which have
2068  the property \f$ \mathrm{level} = y(x) \f$. All points for which
2069  this relation holds are put into the vector \c locs. The
2070  previous information contained in vector \c locs before the
2071  function call is destroyed.
2072 
2073  This is the 1-dimensional analog of finding contour levels as
2074  the \ref contour class does for 2 dimensions.
2075 
2076  This function will call the error handler if \c n is
2077  less than two.
2078 
2079  This function is inclusive of the endpoints, so that if either
2080  <tt>y[0]</tt> and/or <tt>y[n-1]</tt> is equal to level, then
2081  <tt>x[0]</tt> and/or <tt>x[n-1]</tt> are added to \c locs,
2082  respectively.
2083  */
2084  template<class vec_t, class vec2_t> void vector_find_level
2085  (double level, size_t n, vec_t &x, vec2_t &y, std::vector<double> &locs) {
2086 
2087  if (n<=1) {
2088  O2SCL_ERR2("Need at least two data points in ",
2089  "vector_find_level().",exc_einval);
2090  }
2091 
2092  // Ensure that the location vector is empty
2093  locs.clear();
2094 
2095  // Handle left boundary
2096  if (y[0]==level) {
2097  locs.push_back(x[0]);
2098  }
2099 
2100  // Find intersections
2101  for(size_t i=0;i<n-1;i++) {
2102 
2103  if ((y[i]<level && y[i+1]>level) ||
2104  (y[i]>level && y[i+1]<level)) {
2105  // For each intersection, add the location using linear
2106  // interpolation
2107  double x0=x[i]+(x[i+1]-x[i])*(level-y[i])/(y[i+1]-y[i]);
2108  locs.push_back(x0);
2109  } else if (y[i+1]==level) {
2110  locs.push_back(x[i+1]);
2111  }
2112  }
2113 
2114  return;
2115  }
2116 
2117  /** \brief Compute derivative at all points from an
2118  interpolation object
2119  */
2120  template<class ovec_t, class vec2_t>
2121  void vector_deriv_interp(size_t n, ovec_t &v, vec2_t &dv,
2122  size_t interp_type=itp_linear) {
2123  ovec_t grid(n);
2124  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2125  interp_vec<ovec_t,ovec_t> oi(n,grid,v,interp_type);
2126  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(((double)i));
2127  return;
2128  }
2129 
2130  /** \brief Compute second derivative at all points from an
2131  interpolation object
2132  */
2133  template<class ovec_t, class vec2_t>
2134  void vector_deriv2_interp(size_t n, ovec_t &v, vec2_t &dv,
2135  size_t interp_type=itp_linear) {
2136  ovec_t grid(n);
2137  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2138  interp_vec<ovec_t,ovec_t> oi(n,grid,v,interp_type);
2139  for(size_t i=0;i<n;i++) dv[i]=oi.deriv2(((double)i));
2140  return;
2141  }
2142 
2143  /** \brief Compute derivative at all points from an
2144  interpolation object
2145  */
2146  template<class vec_t, class vec2_t, class vec3_t>
2147  void vector_deriv_xy_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv,
2148  size_t interp_type=itp_linear) {
2149  interp_vec<vec_t,vec2_t> oi(n,vx,vy,interp_type);
2150  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(vx[i]);
2151  return;
2152  }
2153 
2154  /** \brief Compute second derivative at all points from an
2155  interpolation object
2156  */
2157  template<class vec_t, class vec2_t, class vec3_t>
2158  void vector_deriv2_xy_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv,
2159  size_t interp_type=itp_linear) {
2160  interp_vec<vec_t,vec2_t> oi(n,vx,vy,interp_type);
2161  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(vx[i]);
2162  return;
2163  }
2164 
2165  /** \brief Integral of a vector from interpolation object
2166  */
2167  template<class ovec_t>
2168  double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type) {
2169  ovec_t grid(n);
2170  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2171  interp_vec<ovec_t> oi(n,grid,v,interp_type);
2172  return oi.integ(0.0,((double)(n-1)));
2173  }
2174 
2175  /** \brief Compute the integral over <tt>y(x)</tt> using
2176  interpolation
2177 
2178  Assuming vectors \c y and \c x define a function \f$ y(x) \f$
2179  then this computes
2180  \f[
2181  \int_{x_0}^{x^{n-1}} y(x)~dx
2182  \f]
2183  using interpolation to approximate the result.
2184 
2185  See also \ref vector_deriv_interp() and
2186  \ref vector_deriv2_interp() in \ref vector_derint.h .
2187  */
2188  template<class vec_t, class vec2_t>
2189  double vector_integ_xy_interp(size_t n, const vec_t &x, const vec2_t &y,
2190  size_t interp_type=itp_linear) {
2191 
2192  // Interpolation object
2193  interp<vec_t,vec2_t> si(interp_type);
2194 
2195  // Compute full integral
2196  double total=si.integ(x[0],x[n-1],n,x,y);
2197 
2198  return total;
2199  }
2200 
2201  /** \brief Compute integral over <tt>y(x)</tt> and store result
2202  in a vector using interpolation
2203  */
2204  template<class vec_t, class vec2_t, class vec3_t>
2205  void vector_integ_xy_interp(size_t n, const vec_t &x, const vec2_t &y,
2206  vec3_t &iy, size_t interp_type=itp_linear) {
2207 
2208  // Interpolation object
2209  interp<vec_t,vec2_t> si(interp_type);
2210 
2211  // Compute full integral
2212  iy[0]=0.0;
2213  for(size_t i=1;i<n;i++) {
2214  iy[i]=si.integ(x[0],x[i],n,x,y);
2215  }
2216 
2217  return;
2218  }
2219 
2220  /** \brief Compute the integral of a vector using
2221  interpolation up to a specified upper limit
2222  */
2223  template<class ovec_t>
2224  double vector_integ_ul_interp(size_t n, double x2,
2225  ovec_t &v, size_t interp_type) {
2226  ovec_t grid(n);
2227  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2228  interp_vec<ovec_t> oi(n,grid,v,interp_type);
2229  return oi.integ(0.0,x2);
2230  }
2231 
2232  /** \brief Compute the integral over <tt>y(x)</tt> using
2233  interpolation up to a specified upper limit
2234  */
2235  template<class vec_t, class vec2_t>
2236  double vector_integ_ul_xy_interp(size_t n, const vec_t &x,
2237  const vec2_t &y, double x2,
2238  size_t interp_type=itp_linear) {
2239 
2240  // Interpolation object
2241  interp<vec_t,vec2_t> si(interp_type);
2242 
2243  // Compute full integral
2244  double total=si.integ(x[0],x2,n,x,y);
2245 
2246  return total;
2247  }
2248 
2249  /** \brief Compute the endpoints which enclose the regions whose
2250  integral is equal to \c sum
2251 
2252  Defining a new function, \f$ g(y_0) \f$ which takes as input any
2253  y-value, \f$ y_0 \f$ from the function \f$ y(x) \f$ (specified
2254  with the parameters \c x and \c y) and outputs the integral of
2255  the function \f$ y(x) \f$ over all regions where \f$ y(x) > y_0
2256  \f$. This function inverts \f$ g(y) \f$, taking the value
2257  of an integral as input, and returns the corresponding y-value
2258  in the variable \c lev.
2259 
2260  This function is particularly useful, for example, in computing
2261  the region which defines 68\% around a peak of data, thus
2262  providing \f$ 1~\sigma \f$ confidence limits.
2263 
2264  By default, this function does not allow any enclosed regions to
2265  go beyond the x region specified by the data. In some cases, it
2266  is useful to fix the boundaries to zero to ensure the integral
2267  is well-defined. If \c boundaries is set to 1, then the LHS
2268  boundary is set to zero, if \c boundaries is set to 2, then the
2269  RHS boundary is set to zero, and if \c boundaries is set to 3,
2270  then both boundaries are set to zero.
2271 
2272  Even if the boundaries are set to zero, the region enclosing a
2273  particular integral may not be well-defined, and this function
2274  can fail to find a region given a specified value of \c sum.
2275  Linear interpolation is used to describe the function \f$ g \f$,
2276  and the precision of this function is limited by this
2277  assumption. This function may also sometimes fail if \c sum is
2278  very close to the minimum or maximum value of the function \f$ g
2279  \f$.
2280 
2281  \comment
2282  Note that the two vector types for x and y must be the
2283  same in order to use o2scl_interp.
2284  \endcomment
2285  */
2286  template<class vec_t, class vec2_t> int vector_invert_enclosed_sum
2287  (double sum, size_t n, vec_t &x, vec2_t &y, double &lev,
2288  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2289 
2290  if (n<=1) {
2291  O2SCL_ERR2("Need at least two data points in ",
2292  "vector_invert_enclosed_sum().",exc_einval);
2293  }
2294 
2296 
2297  // Handle boundaries
2298  std::vector<double> x2, y2;
2299  size_t n2;
2300  if (boundaries==1) {
2301  if (verbose>0) {
2302  std::cout << "Fix left boundary to zero." << std::endl;
2303  }
2304  x2.resize(n+1);
2305  y2.resize(n+1);
2306  x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2307  y2[0]=0.0;
2308  for(size_t i=0;i<n;i++) {
2309  x2[i+1]=x[i];
2310  y2[i+1]=y[i];
2311  }
2312  n2=n+1;
2313  } else if (boundaries==2) {
2314  if (verbose>0) {
2315  std::cout << "Fix right boundary to zero." << std::endl;
2316  }
2317  x2.resize(n+1);
2318  y2.resize(n+1);
2319  for(size_t i=0;i<n;i++) {
2320  x2[i]=x[i];
2321  y2[i]=y[i];
2322  }
2323  x2[n]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2324  y2[n]=0.0;
2325  n2=n+1;
2326  } else if (boundaries==3) {
2327  if (verbose>0) {
2328  std::cout << "Fix both boundaries to zero." << std::endl;
2329  }
2330  x2.resize(n+2);
2331  y2.resize(n+2);
2332  x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2333  y2[0]=0.0;
2334  for(size_t i=0;i<n;i++) {
2335  x2[i+1]=x[i];
2336  y2[i+1]=y[i];
2337  }
2338  x2[n+1]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2339  y2[n+1]=0.0;
2340  n2=n+2;
2341  } else {
2342  if (verbose>0) {
2343  std::cout << "No boundary extrapolation." << std::endl;
2344  }
2345  x2.resize(n);
2346  y2.resize(n);
2347  for(size_t i=0;i<n;i++) {
2348  x2[i]=x[i];
2349  y2[i]=y[i];
2350  }
2351  n2=n;
2352  }
2353 
2354  // Construct a sorted list of function values
2355  ubvector ysort(n2);
2356  vector_copy(n2,y2,ysort);
2357  vector_sort<ubvector,double>(ysort.size(),ysort);
2358 
2359  // Create list of y-values to perform y-value and integral
2360  // interpolation. We choose values in between the grid points to
2361  // ensure that we don't accidentally land precisely on a peak or
2362  // valley.
2363  std::vector<double> ylist;
2364  for(size_t i=0;i<ysort.size()-1;i++) {
2365  if (ysort[i]!=ysort[i+1]) {
2366  ylist.push_back((ysort[i+1]+3.0*ysort[i])/4.0);
2367  ylist.push_back((ysort[i+1]*3.0+ysort[i])/4.0);
2368  }
2369  }
2370 
2371  // Interpolation objects. We need two, one for the
2372  // user-specified vector type, and one for std::vector<double>
2374  interp<std::vector<double>,std::vector<double> > itp2(itp_linear);
2375 
2376  // Construct vectors for interpolation
2377  std::vector<double> xi, yi;
2378 
2379  size_t nfail=0;
2380 
2381  for(size_t k=0;k<ylist.size();k++) {
2382  double lev_tmp=ylist[k];
2383  std::vector<double> locs;
2384  vector_find_level(lev_tmp,n2,x2,y2,locs);
2385  if ((locs.size()%2)!=0) {
2386  nfail++;
2387  if (verbose>0) {
2388  std::cout << k << " " << lev_tmp << " " << 0.0 << " "
2389  << locs.size() << " (fail)" << std::endl;
2390  }
2391  } else {
2392  double sum_temp=0.0;
2393  for(size_t i=0;i<locs.size()/2;i++) {
2394  double x0=locs[2*i];
2395  double x1=locs[2*i+1];
2396  sum_temp+=itp2.integ(x0,x1,n2,x2,y2);
2397  }
2398  xi.push_back(sum_temp);
2399  yi.push_back(lev_tmp);
2400  if (verbose>0) {
2401  std::cout << k << " " << lev_tmp << " " << sum_temp << " "
2402  << locs.size() << " ";
2403  for(size_t i=0;i<locs.size();i++) {
2404  std::cout << locs[i] << " ";
2405  }
2406  std::cout << std::endl;
2407  }
2408  }
2409  }
2410  if (nfail>10) {
2411  if (err_on_fail) {
2412  O2SCL_ERR2("More than 10 failures in ",
2413  "vector_invert_enclosed_sum().",o2scl::exc_efailed);
2414  }
2415  return o2scl::exc_efailed;
2416  }
2417 
2418  lev=itp2.eval(sum,xi.size(),xi,yi);
2419 
2420  return 0;
2421  }
2422 
2423  /** \brief Find the region enclosing an integral
2424  */
2425  template<class vec_t, class vec2_t> int vector_region_int
2426  (size_t n, vec_t &x, vec2_t &y, double intl, std::vector<double> &locs,
2427  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2428 
2429  // Total integral
2430  double total=vector_integ_xy_interp(n,x,y,itp_linear);
2431  if (verbose>0) {
2432  std::cout << "Total integral: " << total << std::endl;
2433  }
2434  // Specified fractional integral
2435  if (verbose>0) {
2436  std::cout << "Target integral: " << intl << std::endl;
2437  }
2438  // Find correct level
2439  double lev;
2440  int ret=vector_invert_enclosed_sum(intl,n,x,y,lev,
2441  boundaries,verbose,err_on_fail);
2442  if (ret!=0) {
2443  if (err_on_fail) {
2444  O2SCL_ERR2("Failed to find a level which enclosed the ",
2445  "specified integral in vector_region_int().",
2447  }
2448  return o2scl::exc_efailed;
2449  }
2450  if (verbose>0) {
2451  std::cout << "Level from vector_invert: " << lev << std::endl;
2452  }
2453 
2454  // Inverse interpolate to find locations corresponding to level
2455  vector_find_level(lev,n,x,y,locs);
2456  if (verbose>0) {
2457  std::cout << "Locations from vector_find_level: ";
2458  for(size_t i=0;i<locs.size();i++) {
2459  std::cout << locs[i];
2460  if (i!=locs.size()-1) {
2461  std::cout << " ";
2462  }
2463  }
2464  std::cout << std::endl;
2465  }
2466  return 0;
2467  }
2468 
2469  /** \brief Find the region enclosing a partial integral
2470  */
2471  template<class vec_t, class vec2_t> int vector_region_fracint
2472  (size_t n, vec_t &x, vec2_t &y, double frac, std::vector<double> &locs,
2473  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2474  double total=vector_integ_xy_interp(n,x,y,itp_linear);
2475  return vector_region_int(n,x,y,frac*total,locs,boundaries,
2476  verbose,err_on_fail);
2477  }
2478 
2479  /** \brief Find the boundaries of the region enclosing a integral
2480 
2481  This function finds the boundaries of the region which
2482  has integral equal to <tt>frac</tt> times the full
2483  integral from the lower x limit to the upper x limit.
2484  */
2485  template<class vec_t, class vec2_t> int vector_bound_fracint
2486  (size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high,
2487  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2488 
2489  std::vector<double> locs;
2490  int ret=vector_region_fracint(n,x,y,frac,locs,boundaries,
2491  verbose,err_on_fail);
2492  if (locs.size()==0 || ret!=0) {
2493  if (err_on_fail) {
2494  O2SCL_ERR2("Zero level crossings or vector_region_fracint() ",
2495  "failed in vector_bound_sigma().",exc_efailed);
2496  }
2497  return o2scl::exc_efailed;
2498  }
2499  // Return min and max location
2500  size_t ix;
2501  vector_min(locs.size(),locs,ix,low);
2502  vector_max(locs.size(),locs,ix,high);
2503  return 0;
2504  }
2505 
2506  /** \brief Find the boundaries of the region enclosing a integral
2507 
2508  This function finds the boundaries of the region which
2509  has integral equal to <tt>frac</tt> times the full
2510  integral from the lower x limit to the upper x limit.
2511  */
2512  template<class vec_t, class vec2_t> int vector_bound_int
2513  (size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high,
2514  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2515 
2516  std::vector<double> locs;
2517  int ret=vector_region_int(n,x,y,frac,locs,boundaries,
2518  verbose,err_on_fail);
2519  if (locs.size()==0 || ret!=0) {
2520  if (err_on_fail) {
2521  O2SCL_ERR2("Zero level crossings or vector_region_int() ",
2522  "failed in vector_bound_sigma().",exc_efailed);
2523  }
2524  return o2scl::exc_efailed;
2525  }
2526  // Return min and max location
2527  size_t ix;
2528  vector_min(locs.size(),locs,ix,low);
2529  vector_max(locs.size(),locs,ix,high);
2530  return 0;
2531  }
2532 
2533  /** \brief From an (x,y) pair, create a new (x,y) pair using
2534  interpolation where the new x vector is uniformly spaced
2535  */
2536  template<class vec_t, class vec2_t, class vec3_t, class vec4_t>
2537  void rebin_xy(const vec_t &x, const vec2_t &y,
2538  vec3_t &new_x, vec4_t &new_y, size_t n_pts,
2539  size_t interp_type) {
2540 
2541  if (x.size()!=y.size()) {
2542  O2SCL_ERR2("The x and y vectors must have the same size ",
2543  "in rebin_xy().",o2scl::exc_einval);
2544  }
2545  if (n_pts<2) {
2546  O2SCL_ERR2("Number of points must be at least 2 ",
2547  "in rebin_xy().",o2scl::exc_einval);
2548  }
2549 
2550  // Vector sizes
2551  size_t n=x.size();
2552 
2553  // Create the grid and setup new_x
2554  uniform_grid_end<double> ugx(x[0],x[n-1],n_pts-1);
2555  ugx.vector(new_x);
2556 
2557  // Allocate space and interpolate into new_y
2558  new_y.resize(n_pts);
2559  interp_vec<vec3_t,vec4_t> itp(n,x,y,interp_type);
2560  for(size_t i=0;i<n_pts;i++) {
2561  new_y[i]=itp.eval(new_x[i]);
2562  }
2563 
2564  return;
2565  }
2566 
2567  /** \brief From an (x,y) pair, create a new (x,y) pair using
2568  interpolation where the new x vector is uniformly spaced and
2569  test accuracy
2570  */
2571  template<class vec_t, class vec2_t, class vec3_t, class vec4_t>
2572  int rebin_xy(const vec_t &x, const vec2_t &y,
2573  vec3_t &new_x, vec4_t &new_y, size_t n_pts,
2574  size_t interp_type1, size_t interp_type2,
2575  double acc=1.0e-4) {
2576 
2577  if (x.size()!=y.size()) {
2578  O2SCL_ERR2("The x and y vectors must have the same size ",
2579  "in rebin_xy().",o2scl::exc_einval);
2580  }
2581  if (n_pts<2) {
2582  O2SCL_ERR2("Number of points must be at least 2 ",
2583  "in rebin_xy().",o2scl::exc_einval);
2584  }
2585 
2586  // Vector sizes
2587  size_t n=x.size();
2588 
2589  // Create the grid and setup new_x
2590  uniform_grid_end<double> ugx(x[0],x[n-1],n_pts-1);
2591  ugx.vector(new_x);
2592 
2593  // Allocate space and interpolate into new_y
2594  std::vector<double> new_y2(n_pts);
2595  new_y.resize(n_pts);
2596  interp_vec<vec3_t,vec4_t> itp1(n,x,y,interp_type1);
2597  interp_vec<vec3_t,vec4_t> itp2(n,x,y,interp_type2);
2598  for(size_t i=0;i<n_pts;i++) {
2599  new_y[i]=itp1.eval(new_x[i]);
2600  new_y2[i]=itp2.eval(new_x[i]);
2601  }
2602 
2603  double max_y, min_y;
2604  vector_minmax_value(n_pts,new_y,min_y,max_y);
2605  for(size_t i=0;i<n_pts;i++) {
2606  if (fabs(new_y2[i]-new_y[i])/(max_y-min_y)>acc) {
2607  return 1;
2608  }
2609  }
2610 
2611  return 0;
2612  }
2613 
2614  /** \brief Rebin, rescale, sort, and match to \f$ y=x \f$
2615 
2616  Principally used by \ref linear_or_log() .
2617  */
2618  template<class vec_t, class vec2_t>
2619  double linear_or_log_chi2(const vec_t &x, const vec2_t &y) {
2620 
2621  static const size_t n2=11;
2622 
2623  // Rebin into vectors of length 11
2624  std::vector<double> xn, yn;
2625  rebin_xy(x,y,xn,yn,n2,itp_linear);
2626 
2627  // Rescale to [0,1]
2628  double min_y, max_y;
2629  vector_minmax_value(n2,yn,min_y,max_y);
2630  for(size_t i=0;i<n2;i++) {
2631  yn[i]=(yn[i]-min_y)/(max_y-min_y);
2632  }
2633 
2634  // Sort and match to line
2635  vector_sort_double(n2,yn);
2636  double chi2=0.0;
2637  for(size_t i=0;i<n2;i++) {
2638  double ty=((double)i)/10.0;
2639  chi2+=pow(yn[i]-ty,2.0);
2640  }
2641 
2642  return chi2;
2643  }
2644 
2645  /** \brief Attempt to determine if data represented by (x,y)
2646  would be better plotted on a semi-log or log-log plot
2647 
2648  \note Experimental.
2649 
2650  This function attempts to guess whether the data stored in \c x
2651  and \c y might be best plotted on a log scale. This algorithm
2652  will fail for poorly sampled or highly oscillatory data.
2653  */
2654  template<class vec_t, class vec2_t>
2655  void linear_or_log(vec_t &x, vec2_t &y, bool &log_x, bool &log_y) {
2656  if (x.size()!=y.size()) {
2657  O2SCL_ERR2("The x and y vectors must have the same size ",
2658  "in linear_or_log().",o2scl::exc_einval);
2659  }
2660 
2661  // Vector sizes
2662  size_t n=x.size();
2663  if (n<2) {
2664  O2SCL_ERR2("Vector size must be at least 2 ",
2665  "in linear_or_log().",o2scl::exc_einval);
2666  }
2667 
2668  // Initial values of log_x and log_y
2669  log_x=false;
2670  log_y=false;
2671 
2672  // Check if vectors are positive
2673  bool x_positive=true;
2674  bool y_positive=true;
2675  for(size_t i=0;i<n;i++) {
2676  if (x[i]<=0.0) x_positive=false;
2677  if (y[i]<=0.0) y_positive=false;
2678  }
2679 
2680  if (x_positive==false && y_positive==false) return;
2681 
2682  double chi2=linear_or_log_chi2(x,y);
2683 
2684  // Take the log
2685  std::vector<double> lx(n), ly(n);
2686  if (x_positive) {
2687  for(size_t i=0;i<n;i++) {
2688  lx[i]=log(x[i]);
2689  }
2690  }
2691  if (y_positive) {
2692  for(size_t i=0;i<n;i++) {
2693  ly[i]=log(y[i]);
2694  }
2695  }
2696 
2697  // Depending on whether or not they are positive, find the best match
2698  if (x_positive) {
2699  if (y_positive) {
2700  double chi2_x=linear_or_log_chi2(lx,y);
2701  double chi2_y=linear_or_log_chi2(x,ly);
2702  double chi2_xy=linear_or_log_chi2(lx,ly);
2703  if (chi2_xy<chi2_x && chi2_xy<chi2_y && chi2_xy<chi2) {
2704  log_x=true;
2705  log_y=true;
2706  } else if (chi2_x<chi2_xy && chi2_x<chi2_y && chi2_x<chi2) {
2707  log_x=true;
2708  } else if (chi2_y<chi2_xy && chi2_y<chi2_x && chi2_y<chi2) {
2709  log_y=true;
2710  }
2711  } else {
2712  double chi2_x=linear_or_log_chi2(lx,y);
2713  if (chi2_x<chi2) log_x=true;
2714  }
2715  } else {
2716  double chi2_y=linear_or_log_chi2(x,ly);
2717  if (chi2_y<chi2) log_y=true;
2718  }
2719 
2720  return;
2721  }
2722 
2723  /** \brief Attempt to determine if data stored in \c y
2724  would be better plotted on a semi-log or log-log plot
2725 
2726  \note Experimental.
2727 
2728  \future There is a bit of extra calculation because the
2729  rebin function creates a second new x vector with a
2730  uniform grid, so this could be streamlined.
2731  */
2732  template<class vec_t>
2733  void linear_or_log(vec_t &y, bool &log_y) {
2734  std::vector<double> x(y.size());
2735  for(size_t i=0;i<y.size();i++) {
2736  x[i]=((double)i);
2737  }
2738  bool log_x;
2739  linear_or_log(x,y,log_x,log_y);
2740  return;
2741  }
2742 
2743 #ifndef DOXYGEN_NO_O2NS
2744 }
2745 #endif
2746 
2747 #endif
double linear_or_log_chi2(const vec_t &x, const vec2_t &y)
Rebin, rescale, sort, and match to .
Definition: interp.h:2619
virtual const char * type() const
Return the type, "interp_steffen".
Definition: interp.h:1357
virtual const char * type() const
Return the type, "interp_cspline_peri".
Definition: interp.h:711
Interpolation class for pre-specified vector.
Definition: interp.h:1795
o2scl_linalg::ubvector_5_mem p5m
Memory for the tridiagonalization.
Definition: interp.h:693
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:2158
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:240
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, , as specified as the first n elements of vectors x and y...
Definition: interp.h:1699
interp_vec(size_t n, const vec_t &x, const vec2_t &y, size_t interp_type=itp_cspline)
Create an interpolation object with interpolation type itp_cspline based on the first n entries of ve...
Definition: interp.h:1822
virtual const char * type() const
Return the type, "interp_linear".
Definition: interp.h:329
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:136
virtual double operator()(double x0) const
Give the value of the function .
Definition: interp.h:174
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:1668
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:2121
void akima_calc(const vec_t &x_array, size_t size, ubvector &umx)
For initializing the interpolation.
Definition: interp.h:853
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:2287
size_t min_size
The minimum size of the vectors to interpolate between.
Definition: interp.h:165
A specialization of interp for C-style double arrays.
Definition: interp.h:1987
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:1515
const double pi
Definition: constants.h:62
interp_base< vec_t, vec2_t > * itp
Pointer to base interpolation object.
Definition: interp.h:1661
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:1949
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:456
Akima spline interpolation with periodic boundary conditions (GSL)
Definition: interp.h:1056
virtual double integ(double al, double bl) const
Give the value of the integral .
Definition: interp.h:1306
void clear()
Manually clear the pointer to the user-specified vector.
Definition: interp.h:1911
virtual const char * type() const
Return the type, "interp_cspline".
Definition: interp.h:664
virtual const char * type() const
Return the type, "interp_akima".
Definition: interp.h:1038
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:2236
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:2486
virtual const char * type() const
Return the type, "interp_monotonic".
Definition: interp.h:1627
double copysign(const double x, const double y)
Flip the sign of x if x and y have different signs.
Definition: interp.h:1166
size_t sz
Vector size.
Definition: interp.h:131
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:1492
invalid argument supplied by user
Definition: err_hnd.h:59
const vec_t * px
Independent vector.
Definition: interp.h:125
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:1268
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:396
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:955
Monotonicity-preserving interpolation.
Definition: interp.h:1395
Base low-level interpolation class [abstract base].
Definition: interp.h:107
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:2013
virtual const char * type() const
Return the type, "interp_akima_peri".
Definition: interp.h:1076
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:2189
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:1282
virtual double deriv(const double x0) const
Give the value of the derivative .
Definition: interp.h:1938
Akima spline for periodic boundary conditions.
Definition: interp.h:79
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:1254
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:609
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:1541
interp_vec()
Blank interpolator.
Definition: interp.h:1813
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:2134
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:986
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:1295
Steffen&#39;s monotonic method.
Definition: interp.h:83
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:972
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:1564
Steffen&#39;s monotonicity-preserving interpolation.
Definition: interp.h:1136
interp_akima()
Create a base interpolation object with or without periodic boundary conditions.
Definition: interp.h:891
Cubic spline interpolation with periodic boundary conditions (GSL)
Definition: interp.h:683
virtual const char * type() const
Return the type, "interp_nearest_neigh".
Definition: interp.h:401
virtual double operator()(double x0) const
Give the value of the function .
Definition: interp.h:1929
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:565
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 , where is specified in the first n elements of vectors x and y...
Definition: interp.h:1730
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:278
Cubic spline for natural boundary conditions.
Definition: interp.h:73
double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type)
Integral of a vector from interpolation object.
Definition: interp.h:2168
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:255
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
ubvector beta
Staggered ratio.
Definition: interp.h:1415
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:1967
interp_array()
Create an interpolator using the default base interpolation objects.
Definition: interp.h:1999
void linear_or_log(vec_t &x, vec2_t &y, bool &log_x, bool &log_y)
Attempt to determine if data represented by (x,y) would be better plotted on a semi-log or log-log pl...
Definition: interp.h:2655
#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:2472
size_t itype
Interpolation type.
Definition: interp.h:1806
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:1958
void vector_minmax_value(size_t n, vec_t &data, data_t &min, data_t &max)
Compute the minimum and maximum of the first n elements of a vector.
Definition: vector.h:1275
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:940
const vec2_t * py
Dependent vector.
Definition: interp.h:128
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:377
o2scl_linalg::ubvector_4_mem p4m
Memory for the tridiagonalization.
Definition: interp.h:453
#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:2147
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:1737
size_t vector_level_count(double level, size_t n, vec_t &x, vec2_t &y)
Count level crossings.
Definition: interp.h:2040
Akima spline for natural boundary conditions.
Definition: interp.h:77
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:2513
Akima spline interpolation (GSL)
Definition: interp.h:829
size_t ordered_lookup(const double x0) const
Find the index of x0 in the ordered array x.
Definition: search_vec.h:242
interp_steffen()
Create a base interpolation object.
Definition: interp.h:1178
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:588
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:1179
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:2426
interp_cspline()
Create a base interpolation object with natural or periodic boundary conditions.
Definition: interp.h:475
void rebin_xy(const vec_t &x, const vec2_t &y, vec3_t &new_x, vec4_t &new_y, size_t n_pts, size_t interp_type)
From an (x,y) pair, create a new (x,y) pair using interpolation where the new x vector is uniformly s...
Definition: interp.h:2537
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:384
Cubic spline for periodic boundary conditions.
Definition: interp.h:75
Monotonicity-preserving interpolation (unfinished)
Definition: interp.h:81
virtual double deriv2(double x0) const
Give the value of the second derivative (always zero)
Definition: interp.h:391
Nearest-neighbor lookup.
Definition: interp.h:85
void vector(resize_vec_t &v) const
Fill a vector with the specified grid.
Definition: uniform_grid.h:246
virtual double deriv2(const double x0, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the second derivative, , where is specified in the first n elements of vectors x ...
Definition: interp.h:1720
ubvector Delta
Finite differences.
Definition: interp.h:1411
Cubic spline interpolation (GSL)
Definition: interp.h:425
interp_array(size_t interp_type)
Create with base interpolation objects it and rit.
Definition: interp.h:1993
Nearest-neighbor interpolation.
Definition: interp.h:347
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:2224
search_vec< const vec_t > svx
To perform binary searches.
Definition: interp.h:122
ubvector alpha
Ratio.
Definition: interp.h:1413
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:1920
static const double x1[5]
Definition: inte_qng_gsl.h:48
Linear interpolation (GSL)
Definition: interp.h:210
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:273
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:2085
Linear.
Definition: interp.h:71
ubvector m
Slopes.
Definition: interp.h:1409
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:544
A specialization of o2scl::interp_vec for C-style arrays.
Definition: interp.h:2007
interp_base< vec_t, vec2_t > * itp
Base interpolation object.
Definition: interp.h:1803
void vector_sort_double(size_t n, vec_t &data)
Sort a vector of doubles (in increasing order)
Definition: vector.h:898
Interpolation class for general vectors.
Definition: interp.h:1654
virtual double deriv(const double x0, size_t n, const vec_t &x, const vec2_t &y)
Give the value of the derivative, , where is specified in the first n elements of vectors x and y...
Definition: interp.h:1710
Linear grid with fixed number of bins and fixed endpoint.
Definition: uniform_grid.h:333

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