inte_adapt_cern.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 /** \file inte_adapt_cern.h
24  \brief File defining \ref o2scl::inte_adapt_cern
25 */
26 #ifndef O2SCL_CERN_ADAPT_H
27 #define O2SCL_CERN_ADAPT_H
28 
29 #include <o2scl/misc.h>
30 #include <o2scl/inte.h>
31 #include <o2scl/inte_gauss56_cern.h>
32 #include <o2scl/string_conv.h>
33 #include <o2scl/vector_derint.h>
34 
35 #ifndef DOXYGEN_NO_O2NS
36 namespace o2scl {
37 #endif
38 
39 #ifdef O2SCL_NEVER_DEFINED
40 
41  // AWS 10/8/19: I think these are unnecessary now, so I'm
42  // removing them with the new ifdef given above.
43 
44  /** \brief A simple precision-agnostic Newton-Cotes integrator
45 
46  \note Experimental.
47  */
48  template<class func_t=funct, class fp_t=double>
49  class inte_newton_cotes : public inte<func_t,fp_t> {
50 
51  public:
52 
53  /** \brief Integrate function \c func from \c a to \c b
54  giving result \c res and error \c err
55  */
56  virtual int integ_err(func_t &func, fp_t a, fp_t b,
57  fp_t &res, fp_t &err) {
58 
59  std::vector<fp_t> x(8), y(8);
60  fp_t seven=7;
61  fp_t dx=b-a;
62  dx/=seven;
63  for(size_t i=0;i<8;i++) {
64  x[i]=((fp_t)i)*dx+a;
65  y[i]=func(x[i]);
66  }
67  res=o2scl::vector_integ_extended8<std::vector<fp_t>,
68  fp_t>(8,y)*dx;
69  fp_t res2=o2scl::vector_integ_extended4<std::vector<fp_t>,
70  fp_t>(8,y)*dx;
71  err=fabs(res2-res);
72 
73  return 0;
74  }
75 
76  };
77 
78  /** \brief A simple precision-agnostic Newton-Cotes integrator
79 
80  \note Experimental.
81  */
82  template<class func_t=funct, class fp_t=double>
83  class inte_newton_cotes2 : public inte<func_t,fp_t> {
84 
85  public:
86 
87  /** \brief Integrate function \c func from \c a to \c b
88  giving result \c res and error \c err
89  */
90  virtual int integ_err(func_t &func, fp_t a, fp_t b,
91  fp_t &res, fp_t &err) {
92 
93  std::vector<fp_t> x(11), y(11);
94  fp_t ten=10;
95  fp_t dx=b-a;
96  dx/=ten;
97  for(size_t i=0;i<11;i++) {
98  x[i]=((fp_t)i)*dx+a;
99  y[i]=func(x[i]);
100  }
101  fp_t prenum=5;
102  fp_t preden=299376;
103  fp_t c1=16067;
104  fp_t c2=106300;
105  fp_t c3=48525;
106  fp_t c4=272400;
107  fp_t c5=260550;
108  fp_t c6=427368;
109  res=dx*prenum/preden*(c1*(y[0]+y[10])+
110  c2*(y[1]+y[9])-
111  c3*(y[2]+y[8])+
112  c4*(y[3]+y[7])-
113  c5*(y[4]+y[6])+
114  c6*y[5]);
115  fp_t nine=9;
116  fp_t dx2=b-a;
117  dx2/=nine;
118  for(size_t i=0;i<10;i++) {
119  x[i]=((fp_t)i)*dx2+a;
120  y[i]=func(x[i]);
121  }
122  prenum=9;
123  preden=89600;
124  c1=2857;
125  c2=15741;
126  c3=1080;
127  c4=19344;
128  c5=5778;
129  fp_t res2=dx2*prenum/preden*(c1*(y[0]+y[9])+
130  c2*(y[1]+y[8])+
131  c3*(y[2]+y[7])+
132  c4*(y[3]+y[6])+
133  c5*(y[4]+y[5]));
134  err=fabs(res2-res);
135 
136  return 0;
137  }
138 
139  };
140 
141  /** \brief A simple precision-agnostic Newton-Cotes integrator
142 
143  \note Experimental.
144  */
145  template<class func_t=funct, class fp_t=double>
146  class inte_newton_cotes_open : public inte<func_t,fp_t> {
147 
148  public:
149 
150  /** \brief Integrate function \c func from \c a to \c b
151  giving result \c res and error \c err
152  */
153  virtual int integ_err(func_t &func, fp_t a, fp_t b,
154  fp_t &res, fp_t &err) {
155 
156  std::vector<fp_t> x(7), y(7);
157  fp_t eight=8;
158  fp_t dx=b-a;
159  dx/=eight;
160  //std::cout << a << " " << b << std::endl;
161  for(size_t i=0;i<7;i++) {
162  x[i]=((fp_t)i+1)*dx+a;
163  y[i]=func(x[i]);
164  //std::cout << x[i] << " " << y[i] << std::endl;
165  }
166  fp_t prenum=8;
167  fp_t preden=945;
168  fp_t c1=460;
169  fp_t c2=-954;
170  fp_t c3=2196;
171  fp_t c4=-2459;
172  res=dx*prenum/preden*(c1*y[0]+c2*y[1]+c3*y[2]+c4*y[3]+
173  c3*y[4]+c2*y[5]+c1*y[6]);
174  //std::cout << res << std::endl;
175  fp_t seven=7;
176  fp_t dx2=b-a;
177  dx2/=seven;
178  for(size_t i=0;i<6;i++) {
179  x[i]=((fp_t)i+1)*dx2+a;
180  y[i]=func(x[i]);
181  //std::cout << x[i] << " " << y[i] << std::endl;
182  }
183  prenum=7;
184  preden=1440;
185  c1=611;
186  c2=-453;
187  c3=562;
188  fp_t res2=dx2*prenum/preden*(c1*(y[0]+y[5])+
189  c2*(y[1]+y[4])+
190  c3*(y[2]+y[3]));
191  //std::cout << res2 << std::endl;
192  //exit(-1);
193  err=fabs(res2-res);
194 
195  return 0;
196  }
197 
198  };
199 
200 #endif
201 
202  /** \brief Adaptive integration (CERNLIB)
203 
204  Uses a base integration object (default is \ref
205  inte_gauss56_cern) to perform adaptive integration by
206  automatically subdividing the integration interval. At each
207  step, the interval with the largest absolute uncertainty is
208  divided in half. The routine succeeds if the absolute tolerance
209  is less than \ref tol_abs or if the relative tolerance is less
210  than \ref tol_rel, i.e.
211  \f[
212  \mathrm{err}\leq\mathrm{tol\_abs}~\mathrm{or}~
213  \mathrm{err}\leq\mathrm{tol\_rel}\cdot|I|
214  \f]
215  where \f$ I \f$ is the current estimate for the integral and \f$
216  \mathrm{err} \f$ is the current estimate for the uncertainty. If
217  the number of subdivisions exceeds the template parameter \c
218  nsub, the error handler is called, since the integration may not
219  have been successful. The number of subdivisions used in the
220  last integration can be obtained from get_nsubdivisions().
221 
222  The template parameter \c nsub, is the maximum number of
223  subdivisions. It is automatically set to 100 in the original
224  CERNLIB routine, and defaults to 100 here. The default base
225  integration object is of type \ref inte_gauss56_cern. This is the
226  CERNLIB default, but can be modified by calling set_inte().
227 
228  This class is based on the CERNLIB routines RADAPT and
229  DADAPT which are documented at
230  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d102/top.html
231 
232  \future
233  - Allow user to set the initial subdivisions?
234  - It might be interesting to directly compare the performance
235  of this class to \ref o2scl::inte_qag_gsl .
236  - There is a fixme entry in the code which could be resolved.
237  - Output the point where most subdividing was required?
238  */
239  template<class func_t=funct,
240  class def_inte_t=inte_gauss56_cern<funct,double,
241  inte_gauss56_coeffs_double>,
242  size_t nsub=100, class fp_t=double>
243  class inte_adapt_cern : public inte<func_t,fp_t> {
244 
245 #ifndef DOXYGEN_INTERNAL
246 
247  protected:
248 
249  /// Lower end of subdivision
250  fp_t xlo[nsub];
251 
252  /// High end of subdivision
253  fp_t xhi[nsub];
254 
255  /// Value of integral for subdivision
256  fp_t tval[nsub];
257 
258  /// Squared error for subdivision
259  fp_t ters[nsub];
260 
261  /// Previous number of subdivisions
263 
264  /// The base integration object
266 
267 #endif
268 
269  public:
270 
271  inte_adapt_cern() {
272  nsubdiv=1;
273  prev_subdiv=0;
274 
275  it=&def_inte;
276  }
277 
278  /// \name Basic usage
279  //@{
280  /** \brief Integrate function \c func from \c a to \c b
281  giving result \c res and error \c err
282  */
283  virtual int integ_err(func_t &func, fp_t a, fp_t b,
284  fp_t &res, fp_t &err) {
285 
286  fp_t tvals=0.0, terss, xlob, xhib, yhib=0.0, te, root=0.0;
287  int i, nsubdivd;
288  fp_t two=2.0;
289 
290  if (nsubdiv==0) {
291  if (prev_subdiv==0) {
292  // If the previous binning was requested, but
293  // there is no previous binning stored, then
294  // just shift to automatic binning
295  nsubdivd=1;
296  } else {
297  tvals=0.0;
298  terss=0.0;
299  for(i=0;i<prev_subdiv;i++) {
300  it->integ_err(func,xlo[i],xhi[i],tval[i],te);
301  ters[i]=te*te;
302  tvals+=tval[i];
303  terss+=ters[i];
304  }
305  err=sqrt(two*terss);
306  res=tvals;
307  return 0;
308  }
309  }
310 
311  // Ensure we're not asking for too many divisions
312  if (nsub<nsubdiv) {
313  nsubdivd=nsub;
314  } else {
315  nsubdivd=nsubdiv;
316  }
317 
318  // Compute the initial set of intervals and integral values
319  xhib=a;
320  fp_t bin=(b-a)/((fp_t)nsubdivd);
321  for(i=0;i<nsubdivd;i++) {
322  xlo[i]=xhib;
323  xlob=xlo[i];
324  xhi[i]=xhib+bin;
325  if (i==nsubdivd-1) xhi[i]=b;
326  xhib=xhi[i];
327  it->integ_err(func,xlob,xhib,tval[i],te);
328  ters[i]=te*te;
329  }
330  prev_subdiv=nsubdivd;
331 
332  for(size_t iter=1;iter<=nsub;iter++) {
333 
334  // Compute the total value of the integrand
335  // and the squared uncertainty
336  tvals=tval[0];
337  terss=ters[0];
338  for(i=1;i<prev_subdiv;i++) {
339  tvals+=tval[i];
340  terss+=ters[i];
341  }
342 
343  // Output iteration information
344  if (this->verbose>0) {
345  std::cout << "inte_adapt_cern Iter: " << iter;
346  std::cout.setf(std::ios::showpos);
347  std::cout << " Res: " << tvals;
348  std::cout.unsetf(std::ios::showpos);
349  std::cout << " Err: " << sqrt(two*terss);
350  if (this->tol_abs>this->tol_rel*abs(tvals)) {
351  std::cout << " Tol: " << this->tol_abs << std::endl;
352  } else {
353  std::cout << " Tol: " << this->tol_rel*abs(tvals)
354  << std::endl;
355  }
356  if (this->verbose>1) {
357  char ch;
358  std::cout << "Press a key and type enter to continue. " ;
359  std::cin >> ch;
360  }
361  }
362 
363  // See if we're finished
364  root=sqrt(two*terss);
365  if (root<=this->tol_abs || root<=this->tol_rel*abs(tvals)) {
366  res=tvals;
367  err=root;
368  this->last_iter=iter;
369  return 0;
370  }
371 
372  // Test if we've run out of intervals
373  if (prev_subdiv==nsub) {
374  res=tvals;
375  err=root;
376  this->last_iter=iter;
377  std::string s="Reached maximum number ("+itos(nsub)+
378  ") of subdivisions in inte_adapt_cern::integ_err().";
379  O2SCL_CONV_RET(s.c_str(),exc_etable,this->err_nonconv);
380  }
381 
382  // Find the subdivision with the largest error
383  fp_t bige=ters[0];
384  int ibig=0;
385  for(i=1;i<prev_subdiv;i++) {
386  if (ters[i]>bige) {
387  bige=ters[i];
388  ibig=i;
389  }
390  }
391 
392  // Subdivide that subdivision further
393  xhi[prev_subdiv]=xhi[ibig];
394  fp_t xnew=(xlo[ibig]+xhi[ibig])/two;
395  xhi[ibig]=xnew;
396  xlo[prev_subdiv]=xnew;
397  it->integ_err(func,xlo[ibig],xhi[ibig],tval[ibig],te);
398  ters[ibig]=te*te;
399  it->integ_err(func,xlo[prev_subdiv],
401  ters[prev_subdiv]=te*te;
402  prev_subdiv++;
403 
404  }
405 
406  // FIXME: Should we set an error here, or does this
407  // only happen if we happen to need exactly nsub
408  // intervals?
409  res=tvals;
410  err=root;
411  return 0;
412  }
413  //@}
414 
415  /// \name Integration object
416  //@{
417  /// Set the base integration object to use
419  it=&i;
420  return 0;
421  }
422 
423  /// Default integration object
424  def_inte_t def_inte;
425  //@}
426 
427  /// \name Subdivisions
428  //@{
429  /** \brief Number of subdivisions
430 
431  The options are
432  - 0: Use previous binning and do not subdivide further \n
433  - 1: Automatic - adapt until tolerance is attained (default) \n
434  - n: (n>1) split first in n equal subdivisions, then adapt
435  until tolerance is obtained.
436  */
437  size_t nsubdiv;
438 
439  /// Return the number of subdivisions used in the last integration
440  size_t get_nsubdivisions() {
441  return prev_subdiv;
442  }
443 
444  /// Return the ith subdivision
445  int get_ith_subdivision(size_t i, fp_t &xlow, fp_t &xhigh,
446  fp_t &value, fp_t &errsq) {
447  if (i<prev_subdiv) {
448  xlow=xlo[i];
449  xhigh=xhi[i];
450  value=tval[i];
451  errsq=ters[i];
452  }
453  return 0;
454  }
455 
456  /// Return all of the subdivisions
457  template<class vec_t>
458  int get_subdivisions(vec_t &xlow, vec_t &xhigh, vec_t &value,
459  vec_t &errsq) {
460 
461  for(int i=0;i<prev_subdiv;i++) {
462  xlow[i]=xlo[i];
463  xhigh[i]=xhi[i];
464  value[i]=tval[i];
465  errsq[i]=ters[i];
466  }
467  return 0;
468  }
469  //@}
470 
471  };
472 
473 #ifdef O2SCL_NEVER_DEFINED
474 
475  /** \brief An experimental adaptive integrator from 0 to \$ \infty \f$
476  based on \ref o2scl::inte_adapt_cern
477 
478  (This class is now replaced by inte_il in inte.h)
479  */
480  template<class func_t=funct,
481  class def_inte_t=inte_gauss56_cern<funct,double,
482  inte_gauss56_coeffs_double>,
483  size_t nsub=100,
484  class fp_t=double>
485  class inte_qagil_cern : public inte<func_t,fp_t> {
486 
487  protected:
488 
489  /// A pointer to the user-specified function
490  func_t *user_func;
491 
492  /// The upper limit
493  fp_t upper_limit;
494 
495  /// Transform to \f$ t \in (0,1] \f$
496  virtual fp_t transform(fp_t t) {
497  fp_t x=upper_limit-(1-t)/t, y;
498  y=(*user_func)(x);
499  /*
500  if (false) {
501  std::cout << upper_limit << " " << t << " " << y/t/t << std::endl;
502  char ch;
503  std::cin >> ch;
504  }
505  */
506  return y/t/t;
507  }
508 
509  public:
510 
511  inte_qagil_cern() {
512  it=&def_inte;
513  fo=std::bind(std::mem_fn<fp_t(fp_t)>(&inte_qagil_cern::transform),
514  this,std::placeholders::_1);
515  }
516 
517  /** \brief Internal function type based on floating-point type
518 
519  \comment
520  This type must be public so the user can change the
521  base integration object
522  \endcomment
523  */
524  typedef std::function<fp_t(fp_t)> internal_funct;
525 
526  protected:
527 
528  /// The base integration object
529  inte<internal_funct,fp_t> *it;
530 
531  /// Function object
532  internal_funct fo;
533 
534  public:
535 
536  /** \brief Integrate function \c func from \c a to \c b
537  giving result \c res and error \c err
538  */
539  virtual int integ_err(func_t &func, fp_t a, fp_t b,
540  fp_t &res, fp_t &err) {
541  user_func=&func;
542  upper_limit=b;
543  int ret=it->integ_err(fo,0.0,1.0,res,err);
544  return ret;
545  }
546 
547  /// \name Integration object
548  //@{
549  /// Set the base integration object to use
550  int set_inte(inte<internal_funct,fp_t> &i) {
551  it=&i;
552  return 0;
553  }
554 
555  /// Default integration object
556  inte_adapt_cern<internal_funct,def_inte_t,nsub,fp_t> def_inte;
557  //@}
558 
559  };
560 
561 #endif
562 
563 #ifndef DOXYGEN_NO_O2NS
564 }
565 #endif
566 
567 #endif
o2scl::inte_adapt_cern::it
inte< func_t, fp_t > * it
The base integration object.
Definition: inte_adapt_cern.h:265
o2scl::inte_adapt_cern::ters
fp_t ters[nsub]
Squared error for subdivision.
Definition: inte_adapt_cern.h:259
o2scl::inte_adapt_cern::tval
fp_t tval[nsub]
Value of integral for subdivision.
Definition: inte_adapt_cern.h:256
o2scl::inte_adapt_cern::nsubdiv
size_t nsubdiv
Number of subdivisions.
Definition: inte_adapt_cern.h:437
o2scl::inte_adapt_cern::prev_subdiv
int prev_subdiv
Previous number of subdivisions.
Definition: inte_adapt_cern.h:262
o2scl::inte_adapt_cern::xlo
fp_t xlo[nsub]
Lower end of subdivision.
Definition: inte_adapt_cern.h:250
o2scl::inte_adapt_cern::integ_err
virtual int integ_err(func_t &func, fp_t a, fp_t b, fp_t &res, fp_t &err)
Integrate function func from a to b giving result res and error err.
Definition: inte_adapt_cern.h:283
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::exc_etable
@ exc_etable
table table limit exceeded
Definition: err_hnd.h:103
o2scl::root
One-dimensional solver [abstract base].
Definition: root.h:48
o2scl::inte_adapt_cern::set_inte
int set_inte(inte< func_t, fp_t > &i)
Set the base integration object to use.
Definition: inte_adapt_cern.h:418
o2scl::inte_adapt_cern::get_nsubdivisions
size_t get_nsubdivisions()
Return the number of subdivisions used in the last integration.
Definition: inte_adapt_cern.h:440
o2scl::inte< funct, double >::last_iter
size_t last_iter
The most recent number of iterations taken.
Definition: inte.h:70
O2SCL_CONV_RET
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:297
o2scl::inte< funct, double >::tol_abs
double tol_abs
The maximum absolute uncertainty in the value of the integral (default )
Definition: inte.h:80
o2scl::inte_adapt_cern::get_ith_subdivision
int get_ith_subdivision(size_t i, fp_t &xlow, fp_t &xhigh, fp_t &value, fp_t &errsq)
Return the ith subdivision.
Definition: inte_adapt_cern.h:445
o2scl::inte< funct, double >::verbose
int verbose
Verbosity.
Definition: inte.h:67
o2scl::inte_adapt_cern::def_inte
def_inte_t def_inte
Default integration object.
Definition: inte_adapt_cern.h:424
o2scl::itos
std::string itos(int x)
Convert an integer to a string.
o2scl::inte< funct, double >::tol_rel
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:75
o2scl::inte
Base integration class [abstract base].
Definition: inte.h:51
o2scl::inte_adapt_cern
Adaptive integration (CERNLIB)
Definition: inte_adapt_cern.h:243
o2scl::inte_adapt_cern::get_subdivisions
int get_subdivisions(vec_t &xlow, vec_t &xhigh, vec_t &value, vec_t &errsq)
Return all of the subdivisions.
Definition: inte_adapt_cern.h:458
o2scl::funct
std::function< double(double)> funct
One-dimensional function typedef in src/base/funct.h.
Definition: funct.h:48
o2scl::inte_adapt_cern::xhi
fp_t xhi[nsub]
High end of subdivision.
Definition: inte_adapt_cern.h:253

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