interp.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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 in src/base/interp.h
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 
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()));
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 
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()));
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 
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 
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=
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
1410  /// Finite differences
1412  /// Ratio
1414  /// Staggered ratio
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 
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;
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  /// \name A function for inverse interpolation in src/base/interp.h
2019  //@{
2020  /** \brief Count level crossings
2021 
2022  This function counts the number of times the function \f$ y(x) =
2023  \mathrm{level} \f$ where the function is defined by the vectors
2024  \c x and \c y.
2025 
2026  The value returned is exactly the same as the size of the
2027  \c locs vector computed by \ref vector_find_level().
2028 
2029  This function will call the error handler if \c n is
2030  less than two.
2031 
2032  If one of the entries in \c y is either larger or smaller than
2033  its neighbors (i.e. if the function \f$ y(x) \f$ has an
2034  extremum), and if the value of \c level is exactly equal to the
2035  extremum, then this is counted as 1 level crossing. The same
2036  applies if either the first or last entry in \c y is exactly
2037  equal to \c level . Finite-precision errors may affect
2038  the result of this function when \c level is nearly
2039  equal to one of the value in the vector \c y.
2040  */
2041  template<class vec_t, class vec2_t> size_t vector_level_count
2042  (double level, size_t n, vec_t &x, vec2_t &y) {
2043 
2044  if (n<=1) {
2045  O2SCL_ERR2("Need at least two data points in ",
2046  "vector_find_count().",exc_einval);
2047  }
2048 
2049  size_t count=0;
2050 
2051  // Handle left boundary
2052  if (y[0]==level) count++;
2053 
2054  // Find intersections
2055  for(size_t i=0;i<n-1;i++) {
2056 
2057  if ((y[i]<level && y[i+1]>=level) ||
2058  (y[i]>level && y[i+1]<=level)) {
2059  count++;
2060  }
2061  }
2062 
2063  return count;
2064  }
2065  //@}
2066 
2067  /// \name Derivatives and integrals of vectors in src/base/interp.h
2068  //@{
2069  /** \brief Compute derivative at all points from an
2070  interpolation object
2071  */
2072  template<class ovec_t, class vec2_t>
2073  void vector_deriv_interp(size_t n, ovec_t &v, vec2_t &dv,
2074  size_t interp_type=itp_linear) {
2075  ovec_t grid(n);
2076  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2077  interp_vec<ovec_t,ovec_t> oi(n,grid,v,interp_type);
2078  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(((double)i));
2079  return;
2080  }
2081 
2082  /** \brief Compute second derivative at all points from an
2083  interpolation object
2084  */
2085  template<class ovec_t, class vec2_t>
2086  void vector_deriv2_interp(size_t n, ovec_t &v, vec2_t &dv,
2087  size_t interp_type=itp_linear) {
2088  ovec_t grid(n);
2089  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2090  interp_vec<ovec_t,ovec_t> oi(n,grid,v,interp_type);
2091  for(size_t i=0;i<n;i++) dv[i]=oi.deriv2(((double)i));
2092  return;
2093  }
2094 
2095  /** \brief Compute derivative at all points from an
2096  interpolation object
2097  */
2098  template<class vec_t, class vec2_t, class vec3_t>
2099  void vector_deriv_xy_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv,
2100  size_t interp_type=itp_linear) {
2101  interp_vec<vec_t,vec2_t> oi(n,vx,vy,interp_type);
2102  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(vx[i]);
2103  return;
2104  }
2105 
2106  /** \brief Compute second derivative at all points from an
2107  interpolation object
2108  */
2109  template<class vec_t, class vec2_t, class vec3_t>
2110  void vector_deriv2_xy_interp(size_t n, vec_t &vx, vec2_t &vy, vec3_t &dv,
2111  size_t interp_type=itp_linear) {
2112  interp_vec<vec_t,vec2_t> oi(n,vx,vy,interp_type);
2113  for(size_t i=0;i<n;i++) dv[i]=oi.deriv(vx[i]);
2114  return;
2115  }
2116 
2117  /** \brief Integral of a vector from interpolation object
2118  */
2119  template<class ovec_t>
2120  double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type) {
2121  ovec_t grid(n);
2122  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2123  interp_vec<ovec_t> oi(n,grid,v,interp_type);
2124  return oi.integ(0.0,((double)(n-1)));
2125  }
2126 
2127  /** \brief Compute the integral over <tt>y(x)</tt> using
2128  interpolation
2129 
2130  Assuming vectors \c y and \c x define a function \f$ y(x) \f$
2131  then this computes
2132  \f[
2133  \int_{x_0}^{x^{n-1}} y(x)~dx
2134  \f]
2135  using interpolation to approximate the result.
2136 
2137  See also \ref vector_deriv_interp() and
2138  \ref vector_deriv2_interp() in \ref vector_derint.h .
2139  */
2140  template<class vec_t, class vec2_t>
2141  double vector_integ_xy_interp(size_t n, const vec_t &x, const vec2_t &y,
2142  size_t interp_type=itp_linear) {
2143 
2144  // Interpolation object
2145  interp<vec_t,vec2_t> si(interp_type);
2146 
2147  // Compute full integral
2148  double total=si.integ(x[0],x[n-1],n,x,y);
2149 
2150  return total;
2151  }
2152 
2153  /** \brief Compute integral over <tt>y(x)</tt> and store result
2154  in a vector using interpolation
2155  */
2156  template<class vec_t, class vec2_t, class vec3_t>
2157  void vector_integ_xy_interp(size_t n, const vec_t &x, const vec2_t &y,
2158  vec3_t &iy, size_t interp_type=itp_linear) {
2159 
2160  // Interpolation object
2161  interp<vec_t,vec2_t> si(interp_type);
2162 
2163  // Compute full integral
2164  iy[0]=0.0;
2165  for(size_t i=1;i<n;i++) {
2166  iy[i]=si.integ(x[0],x[i],n,x,y);
2167  }
2168 
2169  return;
2170  }
2171 
2172  /** \brief Compute the integral of a vector using
2173  interpolation up to a specified upper limit
2174  */
2175  template<class ovec_t>
2176  double vector_integ_ul_interp(size_t n, double x2,
2177  ovec_t &v, size_t interp_type) {
2178  ovec_t grid(n);
2179  for(size_t i=0;i<n;i++) grid[i]=((double)i);
2180  interp_vec<ovec_t> oi(n,grid,v,interp_type);
2181  return oi.integ(0.0,x2);
2182  }
2183 
2184  /** \brief Compute the integral over <tt>y(x)</tt> using
2185  interpolation up to a specified upper limit
2186  */
2187  template<class vec_t, class vec2_t>
2188  double vector_integ_ul_xy_interp(size_t n, const vec_t &x,
2189  const vec2_t &y, double x2,
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],x2,n,x,y);
2197 
2198  return total;
2199  }
2200  //@}
2201 
2202  /// \name Inverse interpolation and related in src/base/interp.h
2203  //@{
2204  /** \brief Perform inverse linear interpolation
2205 
2206  This function performs inverse linear interpolation of the data
2207  defined by \c x and \c y, finding all points in \c x which have
2208  the property \f$ \mathrm{level} = y(x) \f$. All points for which
2209  this relation holds are put into the vector \c locs. The
2210  previous information contained in vector \c locs before the
2211  function call is destroyed.
2212 
2213  This is the 1-dimensional analog of finding contour levels as
2214  the \ref contour class does for 2 dimensions.
2215 
2216  This function will call the error handler if \c n is
2217  less than two.
2218 
2219  This function is inclusive of the endpoints, so that if either
2220  <tt>y[0]</tt> and/or <tt>y[n-1]</tt> is equal to level, then
2221  <tt>x[0]</tt> and/or <tt>x[n-1]</tt> are added to \c locs,
2222  respectively.
2223  */
2224  template<class vec_t, class vec2_t> void vector_find_level
2225  (double level, size_t n, vec_t &x, vec2_t &y, std::vector<double> &locs) {
2226 
2227  if (n<=1) {
2228  O2SCL_ERR2("Need at least two data points in ",
2229  "vector_find_level().",exc_einval);
2230  }
2231 
2232  // Ensure that the location vector is empty
2233  locs.clear();
2234 
2235  // Handle left boundary
2236  if (y[0]==level) {
2237  locs.push_back(x[0]);
2238  }
2239 
2240  // Find intersections
2241  for(size_t i=0;i<n-1;i++) {
2242 
2243  if ((y[i]<level && y[i+1]>level) ||
2244  (y[i]>level && y[i+1]<level)) {
2245  // For each intersection, add the location using linear
2246  // interpolation
2247  double x0=x[i]+(x[i+1]-x[i])*(level-y[i])/(y[i+1]-y[i]);
2248  locs.push_back(x0);
2249  } else if (y[i+1]==level) {
2250  locs.push_back(x[i+1]);
2251  }
2252  }
2253 
2254  return;
2255  }
2256 
2257  /** \brief Compute the endpoints which enclose the regions whose
2258  integral is equal to \c sum
2259 
2260  Defining a new function, \f$ g(y_0) \f$ which takes as input any
2261  y-value, \f$ y_0 \f$ from the function \f$ y(x) \f$ (specified
2262  with the parameters \c x and \c y) and outputs the integral of
2263  the function \f$ y(x) \f$ over all regions where \f$ y(x) > y_0
2264  \f$. This function inverts \f$ g(y) \f$, taking the value
2265  of an integral as input, and returns the corresponding y-value
2266  in the variable \c lev.
2267 
2268  This function is particularly useful, for example, in computing
2269  the region which defines 68\% around a peak of data, thus
2270  providing \f$ 1~\sigma \f$ confidence limits.
2271 
2272  By default, this function does not allow any enclosed regions to
2273  go beyond the x region specified by the data. In some cases, it
2274  is useful to fix the boundaries to zero to ensure the integral
2275  is well-defined. If \c boundaries is set to 1, then the LHS
2276  boundary is set to zero, if \c boundaries is set to 2, then the
2277  RHS boundary is set to zero, and if \c boundaries is set to 3,
2278  then both boundaries are set to zero.
2279 
2280  Even if the boundaries are set to zero, the region enclosing a
2281  particular integral may not be well-defined, and this function
2282  can fail to find a region given a specified value of \c sum.
2283  Linear interpolation is used to describe the function \f$ g \f$,
2284  and the precision of this function is limited by this
2285  assumption. This function may also sometimes fail if \c sum is
2286  very close to the minimum or maximum value of the function \f$ g
2287  \f$.
2288 
2289  \comment
2290  Note that the two vector types for x and y must be the
2291  same in order to use o2scl_interp.
2292  \endcomment
2293  */
2294  template<class vec_t, class vec2_t> int vector_invert_enclosed_sum
2295  (double sum, size_t n, vec_t &x, vec2_t &y, double &lev,
2296  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2297 
2298  if (n<=1) {
2299  O2SCL_ERR2("Need at least two data points in ",
2300  "vector_invert_enclosed_sum().",exc_einval);
2301  }
2302 
2304 
2305  // Handle boundaries
2306  std::vector<double> x2, y2;
2307  size_t n2;
2308  if (boundaries==1) {
2309  if (verbose>0) {
2310  std::cout << "Fix left boundary to zero." << std::endl;
2311  }
2312  x2.resize(n+1);
2313  y2.resize(n+1);
2314  x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2315  y2[0]=0.0;
2316  for(size_t i=0;i<n;i++) {
2317  x2[i+1]=x[i];
2318  y2[i+1]=y[i];
2319  }
2320  n2=n+1;
2321  } else if (boundaries==2) {
2322  if (verbose>0) {
2323  std::cout << "Fix right boundary to zero." << std::endl;
2324  }
2325  x2.resize(n+1);
2326  y2.resize(n+1);
2327  for(size_t i=0;i<n;i++) {
2328  x2[i]=x[i];
2329  y2[i]=y[i];
2330  }
2331  x2[n]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2332  y2[n]=0.0;
2333  n2=n+1;
2334  } else if (boundaries==3) {
2335  if (verbose>0) {
2336  std::cout << "Fix both boundaries to zero." << std::endl;
2337  }
2338  x2.resize(n+2);
2339  y2.resize(n+2);
2340  x2[0]=x[0]-(x[1]-x[0])/1.0e6;
2341  y2[0]=0.0;
2342  for(size_t i=0;i<n;i++) {
2343  x2[i+1]=x[i];
2344  y2[i+1]=y[i];
2345  }
2346  x2[n+1]=x[n-1]+(x[n-1]-x[n-2])/1.0e6;
2347  y2[n+1]=0.0;
2348  n2=n+2;
2349  } else {
2350  if (verbose>0) {
2351  std::cout << "No boundary extrapolation." << std::endl;
2352  }
2353  x2.resize(n);
2354  y2.resize(n);
2355  for(size_t i=0;i<n;i++) {
2356  x2[i]=x[i];
2357  y2[i]=y[i];
2358  }
2359  n2=n;
2360  }
2361 
2362  // Construct a sorted list of function values
2363  ubvector ysort(n2);
2364  vector_copy(n2,y2,ysort);
2365  vector_sort<ubvector,double>(ysort.size(),ysort);
2366 
2367  // Create list of y-values to perform y-value and integral
2368  // interpolation. We choose values in between the grid points to
2369  // ensure that we don't accidentally land precisely on a peak or
2370  // valley.
2371  std::vector<double> ylist;
2372  for(size_t i=0;i<ysort.size()-1;i++) {
2373  if (ysort[i]!=ysort[i+1]) {
2374  ylist.push_back((ysort[i+1]+3.0*ysort[i])/4.0);
2375  ylist.push_back((ysort[i+1]*3.0+ysort[i])/4.0);
2376  }
2377  }
2378 
2379  // Interpolation objects. We need two, one for the
2380  // user-specified vector type, and one for std::vector<double>
2382  interp<std::vector<double>,std::vector<double> > itp2(itp_linear);
2383 
2384  // Construct vectors for interpolation
2385  std::vector<double> xi, yi;
2386 
2387  size_t nfail=0;
2388 
2389  for(size_t k=0;k<ylist.size();k++) {
2390  double lev_tmp=ylist[k];
2391  std::vector<double> locs;
2392  vector_find_level(lev_tmp,n2,x2,y2,locs);
2393  if ((locs.size()%2)!=0) {
2394  nfail++;
2395  if (verbose>0) {
2396  std::cout << k << " " << lev_tmp << " " << 0.0 << " "
2397  << locs.size() << " (fail)" << std::endl;
2398  }
2399  } else {
2400  double sum_temp=0.0;
2401  for(size_t i=0;i<locs.size()/2;i++) {
2402  double x0=locs[2*i];
2403  double x1=locs[2*i+1];
2404  sum_temp+=itp2.integ(x0,x1,n2,x2,y2);
2405  }
2406  xi.push_back(sum_temp);
2407  yi.push_back(lev_tmp);
2408  if (verbose>0) {
2409  std::cout << k << " " << lev_tmp << " " << sum_temp << " "
2410  << locs.size() << " ";
2411  for(size_t i=0;i<locs.size();i++) {
2412  std::cout << locs[i] << " ";
2413  }
2414  std::cout << std::endl;
2415  }
2416  }
2417  }
2418  if (nfail>10) {
2419  if (err_on_fail) {
2420  O2SCL_ERR2("More than 10 failures in ",
2421  "vector_invert_enclosed_sum().",o2scl::exc_efailed);
2422  }
2423  return o2scl::exc_efailed;
2424  }
2425 
2426  lev=itp2.eval(sum,xi.size(),xi,yi);
2427 
2428  return 0;
2429  }
2430 
2431  /** \brief Find the region enclosing an integral
2432  */
2433  template<class vec_t, class vec2_t> int vector_region_int
2434  (size_t n, vec_t &x, vec2_t &y, double intl, std::vector<double> &locs,
2435  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2436 
2437  // Total integral
2438  double total=vector_integ_xy_interp(n,x,y,itp_linear);
2439  if (verbose>0) {
2440  std::cout << "Total integral: " << total << std::endl;
2441  }
2442  // Specified fractional integral
2443  if (verbose>0) {
2444  std::cout << "Target integral: " << intl << std::endl;
2445  }
2446  // Find correct level
2447  double lev;
2448  int ret=vector_invert_enclosed_sum(intl,n,x,y,lev,
2449  boundaries,verbose,err_on_fail);
2450  if (ret!=0) {
2451  if (err_on_fail) {
2452  O2SCL_ERR2("Failed to find a level which enclosed the ",
2453  "specified integral in vector_region_int().",
2455  }
2456  return o2scl::exc_efailed;
2457  }
2458  if (verbose>0) {
2459  std::cout << "Level from vector_invert: " << lev << std::endl;
2460  }
2461 
2462  // Inverse interpolate to find locations corresponding to level
2463  vector_find_level(lev,n,x,y,locs);
2464  if (verbose>0) {
2465  std::cout << "Locations from vector_find_level: ";
2466  for(size_t i=0;i<locs.size();i++) {
2467  std::cout << locs[i];
2468  if (i!=locs.size()-1) {
2469  std::cout << " ";
2470  }
2471  }
2472  std::cout << std::endl;
2473  }
2474  return 0;
2475  }
2476 
2477  /** \brief Find the region enclosing a partial integral
2478  */
2479  template<class vec_t, class vec2_t> int vector_region_fracint
2480  (size_t n, vec_t &x, vec2_t &y, double frac, std::vector<double> &locs,
2481  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2482  double total=vector_integ_xy_interp(n,x,y,itp_linear);
2483  return vector_region_int(n,x,y,frac*total,locs,boundaries,
2484  verbose,err_on_fail);
2485  }
2486 
2487  /** \brief Find the boundaries of the region enclosing a integral
2488 
2489  This function finds the boundaries of the region which
2490  has integral equal to <tt>frac</tt> times the full
2491  integral from the lower x limit to the upper x limit.
2492  */
2493  template<class vec_t, class vec2_t> int vector_bound_fracint
2494  (size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high,
2495  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2496 
2497  std::vector<double> locs;
2498  int ret=vector_region_fracint(n,x,y,frac,locs,boundaries,
2499  verbose,err_on_fail);
2500  if (locs.size()==0 || ret!=0) {
2501  if (err_on_fail) {
2502  O2SCL_ERR2("Zero level crossings or vector_region_fracint() ",
2503  "failed in vector_bound_sigma().",exc_efailed);
2504  }
2505  return o2scl::exc_efailed;
2506  }
2507  // Return min and max location
2508  size_t ix;
2509  vector_min(locs.size(),locs,ix,low);
2510  vector_max(locs.size(),locs,ix,high);
2511  return 0;
2512  }
2513 
2514  /** \brief Find the boundaries of the region enclosing a integral
2515 
2516  This function finds the boundaries of the region which
2517  has integral equal to <tt>frac</tt> times the full
2518  integral from the lower x limit to the upper x limit.
2519  */
2520  template<class vec_t, class vec2_t> int vector_bound_int
2521  (size_t n, vec_t &x, vec2_t &y, double frac, double &low, double &high,
2522  int boundaries=0, int verbose=0, bool err_on_fail=true) {
2523 
2524  std::vector<double> locs;
2525  int ret=vector_region_int(n,x,y,frac,locs,boundaries,
2526  verbose,err_on_fail);
2527  if (locs.size()==0 || ret!=0) {
2528  if (err_on_fail) {
2529  O2SCL_ERR2("Zero level crossings or vector_region_int() ",
2530  "failed in vector_bound_sigma().",exc_efailed);
2531  }
2532  return o2scl::exc_efailed;
2533  }
2534  // Return min and max location
2535  size_t ix;
2536  vector_min(locs.size(),locs,ix,low);
2537  vector_max(locs.size(),locs,ix,high);
2538  return 0;
2539  }
2540  //@}
2541 
2542  /// \name Binning and log vs. linear in src/base/interp.h
2543  //@{
2544  /** \brief From an (x,y) pair, create a new (x,y) pair using
2545  interpolation where the new x vector is uniformly spaced
2546  */
2547  template<class vec_t, class vec2_t, class vec3_t, class vec4_t>
2548  void rebin_xy(const vec_t &x, const vec2_t &y,
2549  vec3_t &new_x, vec4_t &new_y, size_t n_pts,
2550  size_t interp_type) {
2551 
2552  if (x.size()!=y.size()) {
2553  O2SCL_ERR2("The x and y vectors must have the same size ",
2554  "in rebin_xy().",o2scl::exc_einval);
2555  }
2556  if (n_pts<2) {
2557  O2SCL_ERR2("Number of points must be at least 2 ",
2558  "in rebin_xy().",o2scl::exc_einval);
2559  }
2560 
2561  // Vector sizes
2562  size_t n=x.size();
2563 
2564  // Create the grid and setup new_x
2565  uniform_grid_end<double> ugx(x[0],x[n-1],n_pts-1);
2566  ugx.vector(new_x);
2567 
2568  // Allocate space and interpolate into new_y
2569  new_y.resize(n_pts);
2570  interp_vec<vec3_t,vec4_t> itp(n,x,y,interp_type);
2571  for(size_t i=0;i<n_pts;i++) {
2572  new_y[i]=itp.eval(new_x[i]);
2573  }
2574 
2575  return;
2576  }
2577 
2578  /** \brief From an (x,y) pair, create a new (x,y) pair using
2579  interpolation where the new x vector is uniformly spaced and
2580  test accuracy
2581  */
2582  template<class vec_t, class vec2_t, class vec3_t, class vec4_t>
2583  int rebin_xy(const vec_t &x, const vec2_t &y,
2584  vec3_t &new_x, vec4_t &new_y, size_t n_pts,
2585  size_t interp_type1, size_t interp_type2,
2586  double acc=1.0e-4) {
2587 
2588  if (x.size()!=y.size()) {
2589  O2SCL_ERR2("The x and y vectors must have the same size ",
2590  "in rebin_xy().",o2scl::exc_einval);
2591  }
2592  if (n_pts<2) {
2593  O2SCL_ERR2("Number of points must be at least 2 ",
2594  "in rebin_xy().",o2scl::exc_einval);
2595  }
2596 
2597  // Vector sizes
2598  size_t n=x.size();
2599 
2600  // Create the grid and setup new_x
2601  uniform_grid_end<double> ugx(x[0],x[n-1],n_pts-1);
2602  ugx.vector(new_x);
2603 
2604  // Allocate space and interpolate into new_y
2605  std::vector<double> new_y2(n_pts);
2606  new_y.resize(n_pts);
2607  interp_vec<vec3_t,vec4_t> itp1(n,x,y,interp_type1);
2608  interp_vec<vec3_t,vec4_t> itp2(n,x,y,interp_type2);
2609  for(size_t i=0;i<n_pts;i++) {
2610  new_y[i]=itp1.eval(new_x[i]);
2611  new_y2[i]=itp2.eval(new_x[i]);
2612  }
2613 
2614  double max_y, min_y;
2615  vector_minmax_value(n_pts,new_y,min_y,max_y);
2616  for(size_t i=0;i<n_pts;i++) {
2617  if (fabs(new_y2[i]-new_y[i])/(max_y-min_y)>acc) {
2618  return 1;
2619  }
2620  }
2621 
2622  return 0;
2623  }
2624 
2625  /** \brief Rebin, rescale, sort, and match to \f$ y=x \f$
2626 
2627  Principally used by \ref linear_or_log() .
2628  */
2629  template<class vec_t, class vec2_t>
2630  double linear_or_log_chi2(const vec_t &x, const vec2_t &y) {
2631 
2632  static const size_t n2=11;
2633 
2634  // Rebin into vectors of length 11
2635  std::vector<double> xn, yn;
2636  rebin_xy(x,y,xn,yn,n2,itp_linear);
2637 
2638  // Rescale to [0,1]
2639  double min_y, max_y;
2640  vector_minmax_value(n2,yn,min_y,max_y);
2641  for(size_t i=0;i<n2;i++) {
2642  yn[i]=(yn[i]-min_y)/(max_y-min_y);
2643  }
2644 
2645  // Sort and match to line
2646  vector_sort_double(n2,yn);
2647  double chi2=0.0;
2648  for(size_t i=0;i<n2;i++) {
2649  double ty=((double)i)/10.0;
2650  chi2+=pow(yn[i]-ty,2.0);
2651  }
2652 
2653  return chi2;
2654  }
2655 
2656  /** \brief Attempt to determine if data represented by (x,y)
2657  would be better plotted on a semi-log or log-log plot
2658 
2659  \note Experimental.
2660 
2661  This function attempts to guess whether the data stored in \c x
2662  and \c y might be best plotted on a log scale. This algorithm
2663  will fail for poorly sampled or highly oscillatory data.
2664  */
2665  template<class vec_t, class vec2_t>
2666  void linear_or_log(vec_t &x, vec2_t &y, bool &log_x, bool &log_y) {
2667  if (x.size()!=y.size()) {
2668  O2SCL_ERR2("The x and y vectors must have the same size ",
2669  "in linear_or_log().",o2scl::exc_einval);
2670  }
2671 
2672  // Vector sizes
2673  size_t n=x.size();
2674  if (n<2) {
2675  O2SCL_ERR2("Vector size must be at least 2 ",
2676  "in linear_or_log().",o2scl::exc_einval);
2677  }
2678 
2679  // Initial values of log_x and log_y
2680  log_x=false;
2681  log_y=false;
2682 
2683  // Check if vectors are positive
2684  bool x_positive=true;
2685  bool y_positive=true;
2686  for(size_t i=0;i<n;i++) {
2687  if (x[i]<=0.0) x_positive=false;
2688  if (y[i]<=0.0) y_positive=false;
2689  }
2690 
2691  if (x_positive==false && y_positive==false) return;
2692 
2693  double chi2=linear_or_log_chi2(x,y);
2694 
2695  // Take the log
2696  std::vector<double> lx(n), ly(n);
2697  if (x_positive) {
2698  for(size_t i=0;i<n;i++) {
2699  lx[i]=log(x[i]);
2700  }
2701  }
2702  if (y_positive) {
2703  for(size_t i=0;i<n;i++) {
2704  ly[i]=log(y[i]);
2705  }
2706  }
2707 
2708  // Depending on whether or not they are positive, find the best match
2709  if (x_positive) {
2710  if (y_positive) {
2711  double chi2_x=linear_or_log_chi2(lx,y);
2712  double chi2_y=linear_or_log_chi2(x,ly);
2713  double chi2_xy=linear_or_log_chi2(lx,ly);
2714  if (chi2_xy<chi2_x && chi2_xy<chi2_y && chi2_xy<chi2) {
2715  log_x=true;
2716  log_y=true;
2717  } else if (chi2_x<chi2_xy && chi2_x<chi2_y && chi2_x<chi2) {
2718  log_x=true;
2719  } else if (chi2_y<chi2_xy && chi2_y<chi2_x && chi2_y<chi2) {
2720  log_y=true;
2721  }
2722  } else {
2723  double chi2_x=linear_or_log_chi2(lx,y);
2724  if (chi2_x<chi2) log_x=true;
2725  }
2726  } else {
2727  double chi2_y=linear_or_log_chi2(x,ly);
2728  if (chi2_y<chi2) log_y=true;
2729  }
2730 
2731  return;
2732  }
2733 
2734  /** \brief Refine a vector by interpolating with a second
2735  index vector
2736 
2737  \warning Untested.
2738  */
2739  template<class vec_t, class vec2_t, class data_t>
2740  void vector_refine(size_t n, const vec_t &index, vec2_t &data,
2741  size_t factor, size_t interp_type=itp_linear) {
2742  interp_vec<vec_t,vec2_t> iv(n,index,data,interp_type);
2743  vec2_t copy=data;
2744  data.resize((n-1)*factor+1);
2745  for (size_t j=0;j<n-1;j++) {
2746  for(size_t k=0;k<factor;k++) {
2747  data[j*factor+k]=iv.eval(index[j]+
2748  ((data_t)k)/((data_t)factor)*
2749  (index[j+1]-index[j]));
2750  }
2751  }
2752  data[data.size()-1]=copy[copy.size()-1];
2753  return;
2754  }
2755 
2756  /** \brief Attempt to determine if data stored in \c y
2757  would be better plotted on a semi-log or log-log plot
2758 
2759  \note Experimental.
2760 
2761  \future There is a bit of extra calculation because the
2762  rebin function creates a second new x vector with a
2763  uniform grid, so this could be streamlined.
2764  */
2765  template<class vec_t>
2766  void linear_or_log(vec_t &y, bool &log_y) {
2767  std::vector<double> x(y.size());
2768  for(size_t i=0;i<y.size();i++) {
2769  x[i]=((double)i);
2770  }
2771  bool log_x;
2772  linear_or_log(x,y,log_x,log_y);
2773  return;
2774  }
2775  //@}
2776 
2777 #ifndef DOXYGEN_NO_O2NS
2778 }
2779 #endif
2780 
2781 #endif
o2scl::vector_deriv2_xy_interp
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:2110
o2scl::interp::interp
interp(size_t interp_type=itp_cspline)
Create with base interpolation object it.
Definition: interp.h:1668
o2scl::interp_linear::deriv2
virtual double deriv2(double x0) const
Give the value of the second derivative (always zero)
Definition: interp.h:273
o2scl::vector_refine
void vector_refine(size_t n, const vec_t &index, vec2_t &data, size_t factor, size_t interp_type=itp_linear)
Refine a vector by interpolating with a second index vector.
Definition: interp.h:2740
o2scl::interp_akima_peri::set
virtual void set(size_t size, const vec_t &xa, const vec2_t &ya)
Initialize interpolation routine.
Definition: interp.h:1079
o2scl::interp_vec::type
virtual const char * type() const
Return the type, "interp_vec".
Definition: interp.h:1967
o2scl_linalg::solve_tridiag_sym
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
o2scl::interp_vec::interp_vec
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
o2scl::interp_linear::set
virtual void set(size_t size, const vec_t &x, const vec2_t &y)
Initialize interpolation routine.
Definition: interp.h:226
o2scl::interp_base::py
const vec2_t * py
Dependent vector.
Definition: interp.h:128
o2scl::vector_deriv_xy_interp
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:2099
o2scl::interp_monotonic::deriv2
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:1541
o2scl::dtos
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
Definition: string_conv.h:77
o2scl::interp_vec::operator()
virtual double operator()(double x0) const
Give the value of the function .
Definition: interp.h:1929
o2scl::interp_vec::integ
virtual double integ(const double x1, const double x2) const
Give the value of the integral .
Definition: interp.h:1958
o2scl::interp_cspline_peri::p5m
o2scl_linalg::ubvector_5_mem p5m
Memory for the tridiagonalization.
Definition: interp.h:693
boost::numeric::ublas::vector< double >
o2scl::interp_monotonic::Delta
ubvector Delta
Finite differences.
Definition: interp.h:1411
o2scl::interp_akima
Akima spline interpolation (GSL)
Definition: interp.h:829
o2scl::interp_akima::eval
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:940
o2scl::interp_vec::itype
size_t itype
Interpolation type.
Definition: interp.h:1806
o2scl::interp_nearest_neigh::deriv
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:384
o2scl::itp_akima_peri
@ itp_akima_peri
Akima spline for periodic boundary conditions.
Definition: interp.h:79
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl::interp_cspline_peri::set
virtual void set(size_t size, const vec_t &xa, const vec2_t &ya)
Initialize interpolation routine.
Definition: interp.h:715
o2scl::interp_monotonic::deriv
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:1515
o2scl::rebin_xy
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:2548
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::vector_bound_fracint
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:2494
o2scl::interp_vec::eval
virtual double eval(const double x0) const
Give the value of the function .
Definition: interp.h:1920
o2scl::interp::itp
interp_base< vec_t, vec2_t > * itp
Pointer to base interpolation object.
Definition: interp.h:1661
o2scl::vector_region_fracint
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:2480
o2scl::interp_base::svx
search_vec< const vec_t > svx
To perform binary searches.
Definition: interp.h:122
o2scl::itp_steffen
@ itp_steffen
Steffen's monotonic method.
Definition: interp.h:83
o2scl::interp_cspline
Cubic spline interpolation (GSL)
Definition: interp.h:425
o2scl::linear_or_log_chi2
double linear_or_log_chi2(const vec_t &x, const vec2_t &y)
Rebin, rescale, sort, and match to .
Definition: interp.h:2630
o2scl_const::pi
const double pi
Definition: constants.h:65
o2scl::interp_steffen::deriv
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:1282
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::itp_nearest_neigh
@ itp_nearest_neigh
Nearest-neighbor lookup.
Definition: interp.h:85
o2scl::interp_steffen::deriv2
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:1295
o2scl::interp_array_vec
A specialization of o2scl::interp_vec for C-style arrays.
Definition: interp.h:2007
o2scl::interp_nearest_neigh
Nearest-neighbor interpolation.
Definition: interp.h:347
o2scl::interp_monotonic::eval
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:1492
o2scl::vector_copy
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:124
o2scl::vector_deriv_interp
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:2073
o2scl::interp_cspline::deriv
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:565
o2scl::itp_linear
@ itp_linear
Linear.
Definition: interp.h:71
o2scl::linear_or_log
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:2666
o2scl::interp_cspline::integ
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:609
o2scl::interp_vec::set
void set(size_t n, const vec_t &x, const vec2_t &y)
Modify the interpolation object to operate on the first n entries of vectors x and y.
Definition: interp.h:1866
o2scl::vector_integ_xy_interp
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:2141
o2scl::interp_linear
Linear interpolation (GSL)
Definition: interp.h:210
o2scl::interp_monotonic::integ
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:1564
o2scl::vector_minmax_value
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:1272
o2scl::interp_nearest_neigh::eval
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:377
o2scl::interp_linear::eval
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:240
o2scl::interp::integ
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
o2scl::interp_array::interp_array
interp_array(size_t interp_type)
Create with base interpolation objects it and rit.
Definition: interp.h:1993
o2scl_inte_qng_coeffs::x2
static const double x2[5]
Definition: inte_qng_gsl.h:66
o2scl::interp::set_type
void set_type(size_t interp_type)
Set base interpolation type.
Definition: interp.h:1737
o2scl::interp_vec
Interpolation class for pre-specified vector.
Definition: interp.h:1795
o2scl::interp_monotonic::alpha
ubvector alpha
Ratio.
Definition: interp.h:1413
o2scl::interp_akima::set
virtual void set(size_t size, const vec_t &xa, const vec2_t &ya)
Initialize interpolation routine.
Definition: interp.h:900
o2scl::interp_vec::set
void set(size_t n, const vec_t &x, const vec2_t &y, size_t interp_type)
Set a new vector to interpolate.
Definition: interp.h:1873
o2scl::interp_cspline::coeff_calc
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
o2scl::itp_akima
@ itp_akima
Akima spline for natural boundary conditions.
Definition: interp.h:77
o2scl::interp_nearest_neigh::type
virtual const char * type() const
Return the type, "interp_nearest_neigh".
Definition: interp.h:401
o2scl_inte_qng_coeffs::x1
static const double x1[5]
Definition: inte_qng_gsl.h:48
o2scl::interp_linear::deriv
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:255
o2scl::interp_nearest_neigh::set
virtual void set(size_t size, const vec_t &x, const vec2_t &y)
Initialize interpolation routine.
Definition: interp.h:363
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::vector_region_int
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:2434
o2scl::uniform_grid_end
Linear grid with fixed number of bins and fixed endpoint.
Definition: uniform_grid.h:333
o2scl::interp::deriv2
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
o2scl_linalg::ubvector_5_mem
Allocation object for 5 arrays of equal size.
Definition: tridiag_base.h:86
o2scl::interp::eval
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
o2scl::interp_monotonic
Monotonicity-preserving interpolation.
Definition: interp.h:1395
o2scl::interp_akima::akima_calc
void akima_calc(const vec_t &x_array, size_t size, ubvector &umx)
For initializing the interpolation.
Definition: interp.h:853
o2scl::interp_base::integ_eval
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
o2scl::interp_cspline::deriv2
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:588
o2scl::vector_integ_ul_xy_interp
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:2188
o2scl::vector_find_level
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:2225
o2scl::interp_nearest_neigh::integ
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:396
o2scl::interp_base::min_size
size_t min_size
The minimum size of the vectors to interpolate between.
Definition: interp.h:165
o2scl::interp_akima::type
virtual const char * type() const
Return the type, "interp_akima".
Definition: interp.h:1038
o2scl::search_vec< const vec_t >
o2scl::interp_cspline::set
virtual void set(size_t size, const vec_t &xa, const vec2_t &ya)
Initialize interpolation routine.
Definition: interp.h:484
o2scl::interp_vec::itp
interp_base< vec_t, vec2_t > * itp
Base interpolation object.
Definition: interp.h:1803
o2scl::vector_level_count
size_t vector_level_count(double level, size_t n, vec_t &x, vec2_t &y)
Count level crossings.
Definition: interp.h:2042
o2scl::interp_akima_peri::type
virtual const char * type() const
Return the type, "interp_akima_peri".
Definition: interp.h:1076
o2scl::interp_array_vec::interp_array_vec
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
o2scl::vector_integ_interp
double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type)
Integral of a vector from interpolation object.
Definition: interp.h:2120
o2scl::vector_min
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:1251
o2scl::interp_vec::interp_vec
interp_vec()
Blank interpolator.
Definition: interp.h:1813
o2scl::interp_cspline_peri::type
virtual const char * type() const
Return the type, "interp_cspline_peri".
Definition: interp.h:711
o2scl::interp_monotonic::beta
ubvector beta
Staggered ratio.
Definition: interp.h:1415
o2scl::vector_max
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:1176
o2scl::itp_monotonic
@ itp_monotonic
Monotonicity-preserving interpolation (unfinished)
Definition: interp.h:81
o2scl::interp_steffen::interp_steffen
interp_steffen()
Create a base interpolation object.
Definition: interp.h:1178
o2scl::interp::deriv
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
o2scl::interp_vec::clear
void clear()
Manually clear the pointer to the user-specified vector.
Definition: interp.h:1911
o2scl::interp_akima::deriv
virtual double deriv(double x0) const
Give the value of the derivative .
Definition: interp.h:955
o2scl::interp_cspline::p4m
o2scl_linalg::ubvector_4_mem p4m
Memory for the tridiagonalization.
Definition: interp.h:453
o2scl::vector_sort_double
void vector_sort_double(size_t n, vec_t &data)
Sort a vector of doubles (in increasing order)
Definition: vector.h:895
o2scl::interp_monotonic::m
ubvector m
Slopes.
Definition: interp.h:1409
o2scl::interp_base::operator()
virtual double operator()(double x0) const
Give the value of the function .
Definition: interp.h:174
o2scl::interp_steffen::type
virtual const char * type() const
Return the type, "interp_steffen".
Definition: interp.h:1357
o2scl::interp_akima_peri
Akima spline interpolation with periodic boundary conditions (GSL)
Definition: interp.h:1056
o2scl::interp
Interpolation class for general vectors.
Definition: interp.h:1654
o2scl::interp_steffen::integ
virtual double integ(double al, double bl) const
Give the value of the integral .
Definition: interp.h:1306
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::interp_array::interp_array
interp_array()
Create an interpolator using the default base interpolation objects.
Definition: interp.h:1999
o2scl::interp_monotonic::set
virtual void set(size_t size, const vec_t &x, const vec2_t &y)
Initialize interpolation routine.
Definition: interp.h:1427
o2scl::interp_nearest_neigh::deriv2
virtual double deriv2(double x0) const
Give the value of the second derivative (always zero)
Definition: interp.h:391
o2scl::itp_cspline
@ itp_cspline
Cubic spline for natural boundary conditions.
Definition: interp.h:73
o2scl::interp_monotonic::type
virtual const char * type() const
Return the type, "interp_monotonic".
Definition: interp.h:1627
o2scl::uniform_grid< double >::vector
void vector(resize_vec_t &v) const
Fill a vector with the specified grid.
Definition: uniform_grid.h:246
o2scl_linalg::solve_cyc_tridiag_sym
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
o2scl::interp_cspline_peri
Cubic spline interpolation with periodic boundary conditions (GSL)
Definition: interp.h:683
o2scl::interp_akima::deriv2
virtual double deriv2(double x0) const
Give the value of the second derivative .
Definition: interp.h:972
o2scl::interp_vec::deriv2
virtual double deriv2(const double x0) const
Give the value of the second derivative .
Definition: interp.h:1949
o2scl::interp_steffen::eval
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:1268
o2scl::interp_steffen
Steffen's monotonicity-preserving interpolation.
Definition: interp.h:1136
o2scl::ubvector_range
boost::numeric::ublas::vector_range< boost::numeric::ublas::vector< double > > ubvector_range
A ublas::vector_range typedef for o2scl::tensor_grid and related classes in src/base/tensor_grid....
Definition: tensor_grid.h:73
o2scl::interp_akima::interp_akima
interp_akima()
Create a base interpolation object with or without periodic boundary conditions.
Definition: interp.h:891
o2scl::interp_steffen::set
virtual void set(size_t size, const vec_t &xa, const vec2_t &ya)
Initialize interpolation routine.
Definition: interp.h:1187
o2scl::vector_deriv2_interp
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:2086
o2scl::interp_linear::integ
virtual double integ(double a, double b) const
Give the value of the integral .
Definition: interp.h:278
o2scl::vector_integ_ul_interp
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:2176
o2scl::vector_invert_enclosed_sum
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:2295
o2scl::interp_vec::deriv
virtual double deriv(const double x0) const
Give the value of the derivative .
Definition: interp.h:1938
o2scl::interp_cspline::interp_cspline
interp_cspline()
Create a base interpolation object with natural or periodic boundary conditions.
Definition: interp.h:475
o2scl::vector_bound_int
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:2521
o2scl::interp_base::px
const vec_t * px
Independent vector.
Definition: interp.h:125
o2scl::szttos
std::string szttos(size_t x)
Convert a size_t to a string.
o2scl::interp_base::sz
size_t sz
Vector size.
Definition: interp.h:131
o2scl::interp_akima::integ
virtual double integ(double aa, double bb) const
Give the value of the integral .
Definition: interp.h:986
o2scl::interp_cspline::eval
virtual double eval(double x0) const
Give the value of the function .
Definition: interp.h:544
o2scl::interp_array
A specialization of interp for C-style double arrays.
Definition: interp.h:1987
o2scl::interp_linear::type
virtual const char * type() const
Return the type, "interp_linear".
Definition: interp.h:329
o2scl::itp_cspline_peri
@ itp_cspline_peri
Cubic spline for periodic boundary conditions.
Definition: interp.h:75
o2scl::interp_base
Base low-level interpolation class [abstract base].
Definition: interp.h:107
o2scl::search_vec::set_vec
void set_vec(size_t nn, const vec_t &x)
Set the vector to be searched.
Definition: search_vec.h:107
o2scl::interp_cspline::type
virtual const char * type() const
Return the type, "interp_cspline".
Definition: interp.h:664
o2scl::interp_steffen::copysign
double copysign(const double x, const double y)
Flip the sign of x if x and y have different signs.
Definition: interp.h:1166
o2scl_linalg::ubvector_4_mem
Allocation object for 4 arrays of equal size.
Definition: tridiag_base.h:70

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