inte_adapt_cern.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \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/inte.h>
30 #include <o2scl/inte_gauss56_cern.h>
31 #include <o2scl/string_conv.h>
32 
33 #ifndef DOXYGEN_NO_O2NS
34 namespace o2scl {
35 #endif
36 
37  /** \brief Adaptive integration (CERNLIB)
38 
39  Uses a base integration object (default is \ref
40  inte_gauss56_cern) to perform adaptive integration by
41  automatically subdividing the integration interval. At each
42  step, the interval with the largest absolute uncertainty is
43  divided in half. The routine succeeds if the absolute tolerance
44  is less than \ref tol_abs or if the relative tolerance is less
45  than \ref tol_rel, i.e.
46  \f[
47  \mathrm{err}\leq\mathrm{tol\_abs}~\mathrm{or}~
48  \mathrm{err}\leq\mathrm{tol\_rel}\cdot|I|
49  \f]
50  where \f$ I \f$ is the current estimate for the integral and \f$
51  \mathrm{err} \f$ is the current estimate for the uncertainty. If
52  the number of subdivisions exceeds the template parameter \c
53  nsub, the error handler is called, since the integration may not
54  have been successful. The number of subdivisions used in the
55  last integration can be obtained from get_nsubdivisions().
56 
57  The template parameter \c nsub, is the maximum number of
58  subdivisions. It is automatically set to 100 in the original
59  CERNLIB routine, and defaults to 100 here. The default base
60  integration object is of type \ref inte_gauss56_cern. This is the
61  CERNLIB default, but can be modified by calling set_inte().
62 
63  This class is based on the CERNLIB routines RADAPT and
64  DADAPT which are documented at
65  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d102/top.html
66 
67  \future
68  - Allow user to set the initial subdivisions?
69  - It might be interesting to directly compare the performance
70  of this class to \ref o2scl::inte_qag_gsl .
71  - There is a fixme entry in the code which could be resolved.
72  - Output the point where most subdividing was required?
73  */
74  template<class func_t=funct11, size_t nsub=100>
75  class inte_adapt_cern : public inte<func_t> {
76 
77 #ifndef DOXYGEN_INTERNAL
78 
79  protected:
80 
81  /// Lower end of subdivision
82  double xlo[nsub];
83 
84  /// High end of subdivision
85  double xhi[nsub];
86 
87  /// Value of integral for subdivision
88  double tval[nsub];
89 
90  /// Squared error for subdivision
91  double ters[nsub];
92 
93  /// Previous number of subdivisions
95 
96  /// The base integration object
98 
99 #endif
100 
101  public:
102 
103  inte_adapt_cern() {
104  nsubdiv=1;
105  prev_subdiv=0;
106 
107  it=&def_inte;
108  }
109 
110  /// \name Basic usage
111  //@{
112  /** \brief Integrate function \c func from \c a to \c b
113  giving result \c res and error \c err
114  */
115  virtual int integ_err(func_t &func, double a, double b,
116  double &res, double &err) {
117 
118  double tvals=0.0, terss, xlob, xhib, yhib=0.0, te, root=0.0;
119  int i, nsubdivd;
120 
121  if (nsubdiv==0) {
122  if (prev_subdiv==0) {
123  // If the previous binning was requested, but
124  // there is no previous binning stored, then
125  // just shift to automatic binning
126  nsubdivd=1;
127  } else {
128  tvals=0.0;
129  terss=0.0;
130  for(i=0;i<prev_subdiv;i++) {
131  it->integ_err(func,xlo[i],xhi[i],tval[i],te);
132  ters[i]=te*te;
133  tvals+=tval[i];
134  terss+=ters[i];
135  }
136  err=sqrt(2.0*terss);
137  res=tvals;
138  return 0;
139  }
140  }
141 
142  // Ensure we're not asking for too many divisions
143  if (nsub<nsubdiv) {
144  nsubdivd=nsub;
145  } else {
146  nsubdivd=nsubdiv;
147  }
148 
149  // Compute the initial set of intervals and integral values
150  xhib=a;
151  double bin=(b-a)/((double)nsubdivd);
152  for(i=0;i<nsubdivd;i++) {
153  xlo[i]=xhib;
154  xlob=xlo[i];
155  xhi[i]=xhib+bin;
156  if (i==nsubdivd-1) xhi[i]=b;
157  xhib=xhi[i];
158  it->integ_err(func,xlob,xhib,tval[i],te);
159  ters[i]=te*te;
160  }
161  prev_subdiv=nsubdivd;
162 
163  for(size_t iter=1;iter<=nsub;iter++) {
164 
165  // Compute the total value of the integrand
166  // and the squared uncertainty
167  tvals=tval[0];
168  terss=ters[0];
169  for(i=1;i<prev_subdiv;i++) {
170  tvals+=tval[i];
171  terss+=ters[i];
172  }
173 
174  // Output iteration information
175  if (this->verbose>0) {
176  std::cout << "inte_adapt_cern Iter: " << iter;
177  std::cout.setf(std::ios::showpos);
178  std::cout << " Res: " << tvals;
179  std::cout.unsetf(std::ios::showpos);
180  std::cout << " Err: " << sqrt(2.0*terss);
181  if (this->tol_abs>this->tol_rel*fabs(tvals)) {
182  std::cout << " Tol: " << this->tol_abs << std::endl;
183  } else {
184  std::cout << " Tol: " << this->tol_rel*fabs(tvals) << std::endl;
185  }
186  if (this->verbose>1) {
187  char ch;
188  std::cout << "Press a key and type enter to continue. " ;
189  std::cin >> ch;
190  }
191  }
192 
193  // See if we're finished
194  root=sqrt(2.0*terss);
195  if (root<=this->tol_abs || root<=this->tol_rel*fabs(tvals)) {
196  res=tvals;
197  err=root;
198  this->last_iter=iter;
199  return 0;
200  }
201 
202  // Test if we've run out of intervals
203  if (prev_subdiv==nsub) {
204  res=tvals;
205  err=root;
206  this->last_iter=iter;
207  std::string s="Reached maximum number ("+itos(nsub)+
208  ") of subdivisions in inte_adapt_cern::integ_err().";
209  O2SCL_CONV_RET(s.c_str(),exc_etable,this->err_nonconv);
210  }
211 
212  // Find the subdivision with the largest error
213  double bige=ters[0];
214  int ibig=0;
215  for(i=1;i<prev_subdiv;i++) {
216  if (ters[i]>bige) {
217  bige=ters[i];
218  ibig=i;
219  }
220  }
221 
222  // Subdivide that subdivision further
223  xhi[prev_subdiv]=xhi[ibig];
224  double xnew=(xlo[ibig]+xhi[ibig])/2.0;
225  xhi[ibig]=xnew;
226  xlo[prev_subdiv]=xnew;
227  it->integ_err(func,xlo[ibig],xhi[ibig],tval[ibig],te);
228  ters[ibig]=te*te;
229  it->integ_err(func,xlo[prev_subdiv],
230  xhi[prev_subdiv],tval[prev_subdiv],te);
231  ters[prev_subdiv]=te*te;
232  prev_subdiv++;
233 
234  }
235 
236  // FIXME: Should we set an error here, or does this
237  // only happen if we happen to need exactly nsub
238  // intervals?
239  res=tvals;
240  err=root;
241  return 0;
242  }
243  //@}
244 
245  /// \name Integration object
246  //@{
247  /// Set the base integration object to use
249  it=&i;
250  return 0;
251  }
252 
253  /// Default integration object
255  //@}
256 
257  /// \name Subdivisions
258  //@{
259  /** \brief Number of subdivisions
260 
261  The options are
262  - 0: Use previous binning and do not subdivide further \n
263  - 1: Automatic - adapt until tolerance is attained (default) \n
264  - n: (n>1) split first in n equal subdivisions, then adapt
265  until tolerance is obtained.
266  */
267  size_t nsubdiv;
268 
269  /// Return the number of subdivisions used in the last integration
270  size_t get_nsubdivisions() {
271  return prev_subdiv;
272  }
273 
274  /// Return the ith subdivision
275  int get_ith_subdivision(size_t i, double &xlow, double &xhigh,
276  double &value, double &errsq) {
277  if (i<prev_subdiv) {
278  xlow=xlo[i];
279  xhigh=xhi[i];
280  value=tval[i];
281  errsq=ters[i];
282  }
283  return 0;
284  }
285 
286  /// Return all of the subdivisions
287  template<class vec_t>
288  int get_subdivisions(vec_t &xlow, vec_t &xhigh, vec_t &value,
289  vec_t &errsq) {
290 
291  for(int i=0;i<prev_subdiv;i++) {
292  xlow[i]=xlo[i];
293  xhigh[i]=xhi[i];
294  value[i]=tval[i];
295  errsq[i]=ters[i];
296  }
297  return 0;
298  }
299  //@}
300 
301  };
302 
303 #ifndef DOXYGEN_NO_O2NS
304 }
305 #endif
306 
307 #endif
inte_gauss56_cern< func_t > def_inte
Default integration object.
int prev_subdiv
Previous number of subdivisions.
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:292
virtual int integ_err(func_t &func, double a, double b, double &res, double &err)=0
Integrate function func from a to b and place the result in res and the error in err.
double xhi[nsub]
High end of subdivision.
virtual int integ_err(func_t &func, double a, double b, double &res, double &err)
Integrate function func from a to b giving result res and error err.
Adaptive integration (CERNLIB)
inte< func_t > * it
The base integration object.
size_t last_iter
The most recent number of iterations taken.
Definition: inte.h:63
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
Definition: inte.h:81
table table limit exceeded
Definition: err_hnd.h:103
size_t get_nsubdivisions()
Return the number of subdivisions used in the last integration.
double tol_abs
The maximum absolute uncertainty in the value of the integral (default )
Definition: inte.h:73
int get_ith_subdivision(size_t i, double &xlow, double &xhigh, double &value, double &errsq)
Return the ith subdivision.
One-dimensional solver [abstract base].
Definition: root.h:47
Base integration class [abstract base].
Definition: inte.h:44
5,6-point Gaussian quadrature (CERNLIB)
int set_inte(inte< func_t > &i)
Set the base integration object to use.
double ters[nsub]
Squared error for subdivision.
double tval[nsub]
Value of integral for subdivision.
size_t nsubdiv
Number of subdivisions.
double xlo[nsub]
Lower end of subdivision.
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:68
std::string itos(int x)
Convert an integer to a string.
int verbose
Verbosity.
Definition: inte.h:60
int get_subdivisions(vec_t &xlow, vec_t &xhigh, vec_t &value, vec_t &errsq)
Return all of the subdivisions.

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