inte_adapt_cern.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \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=funct, size_t nsub=100, class fp_t=double>
75  class inte_adapt_cern : public inte<func_t,fp_t> {
76 
77 #ifndef DOXYGEN_INTERNAL
78 
79  protected:
80 
81  /// Lower end of subdivision
82  fp_t xlo[nsub];
83 
84  /// High end of subdivision
85  fp_t xhi[nsub];
86 
87  /// Value of integral for subdivision
88  fp_t tval[nsub];
89 
90  /// Squared error for subdivision
91  fp_t 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, fp_t a, fp_t b,
116  fp_t &res, fp_t &err) {
117 
118  fp_t 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  fp_t bin=(b-a)/((fp_t)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*std::abs(tvals)) {
182  std::cout << " Tol: " << this->tol_abs << std::endl;
183  } else {
184  std::cout << " Tol: " << this->tol_rel*std::abs(tvals)
185  << std::endl;
186  }
187  if (this->verbose>1) {
188  char ch;
189  std::cout << "Press a key and type enter to continue. " ;
190  std::cin >> ch;
191  }
192  }
193 
194  // See if we're finished
195  root=sqrt(2.0*terss);
196  if (root<=this->tol_abs || root<=this->tol_rel*std::abs(tvals)) {
197  res=tvals;
198  err=root;
199  this->last_iter=iter;
200  return 0;
201  }
202 
203  // Test if we've run out of intervals
204  if (prev_subdiv==nsub) {
205  res=tvals;
206  err=root;
207  this->last_iter=iter;
208  std::string s="Reached maximum number ("+itos(nsub)+
209  ") of subdivisions in inte_adapt_cern::integ_err().";
210  O2SCL_CONV_RET(s.c_str(),exc_etable,this->err_nonconv);
211  }
212 
213  // Find the subdivision with the largest error
214  fp_t bige=ters[0];
215  int ibig=0;
216  for(i=1;i<prev_subdiv;i++) {
217  if (ters[i]>bige) {
218  bige=ters[i];
219  ibig=i;
220  }
221  }
222 
223  // Subdivide that subdivision further
224  xhi[prev_subdiv]=xhi[ibig];
225  fp_t xnew=(xlo[ibig]+xhi[ibig])/2.0;
226  xhi[ibig]=xnew;
227  xlo[prev_subdiv]=xnew;
228  it->integ_err(func,xlo[ibig],xhi[ibig],tval[ibig],te);
229  ters[ibig]=te*te;
230  it->integ_err(func,xlo[prev_subdiv],
231  xhi[prev_subdiv],tval[prev_subdiv],te);
232  ters[prev_subdiv]=te*te;
233  prev_subdiv++;
234 
235  }
236 
237  // FIXME: Should we set an error here, or does this
238  // only happen if we happen to need exactly nsub
239  // intervals?
240  res=tvals;
241  err=root;
242  return 0;
243  }
244  //@}
245 
246  /// \name Integration object
247  //@{
248  /// Set the base integration object to use
250  it=&i;
251  return 0;
252  }
253 
254  /// Default integration object
256  //@}
257 
258  /// \name Subdivisions
259  //@{
260  /** \brief Number of subdivisions
261 
262  The options are
263  - 0: Use previous binning and do not subdivide further \n
264  - 1: Automatic - adapt until tolerance is attained (default) \n
265  - n: (n>1) split first in n equal subdivisions, then adapt
266  until tolerance is obtained.
267  */
268  size_t nsubdiv;
269 
270  /// Return the number of subdivisions used in the last integration
271  size_t get_nsubdivisions() {
272  return prev_subdiv;
273  }
274 
275  /// Return the ith subdivision
276  int get_ith_subdivision(size_t i, fp_t &xlow, fp_t &xhigh,
277  fp_t &value, fp_t &errsq) {
278  if (i<prev_subdiv) {
279  xlow=xlo[i];
280  xhigh=xhi[i];
281  value=tval[i];
282  errsq=ters[i];
283  }
284  return 0;
285  }
286 
287  /// Return all of the subdivisions
288  template<class vec_t>
289  int get_subdivisions(vec_t &xlow, vec_t &xhigh, vec_t &value,
290  vec_t &errsq) {
291 
292  for(int i=0;i<prev_subdiv;i++) {
293  xlow[i]=xlo[i];
294  xhigh[i]=xhi[i];
295  value[i]=tval[i];
296  errsq[i]=ters[i];
297  }
298  return 0;
299  }
300  //@}
301 
302  };
303 
304 #ifndef DOXYGEN_NO_O2NS
305 }
306 #endif
307 
308 #endif
fp_t tol_abs
The maximum absolute uncertainty in the value of the integral (default )
Definition: inte.h:80
int get_ith_subdivision(size_t i, fp_t &xlow, fp_t &xhigh, fp_t &value, fp_t &errsq)
Return the ith subdivision.
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:297
int prev_subdiv
Previous number of subdivisions.
fp_t tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:75
size_t get_nsubdivisions()
Return the number of subdivisions used in the last integration.
fp_t tval[nsub]
Value of integral for subdivision.
virtual int integ_err(func_t &func, fp_t a, fp_t b, fp_t &res, fp_t &err)=0
Integrate function func from a to b and place the result in res and the error in err.
Adaptive integration (CERNLIB)
table table limit exceeded
Definition: err_hnd.h:103
int set_inte(inte< func_t, fp_t > &i)
Set the base integration object to use.
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
Definition: inte.h:88
inte_gauss56_cern< func_t, fp_t > def_inte
Default integration object.
int get_subdivisions(vec_t &xlow, vec_t &xhigh, vec_t &value, vec_t &errsq)
Return all of the subdivisions.
One-dimensional solver [abstract base].
Definition: root.h:47
size_t last_iter
The most recent number of iterations taken.
Definition: inte.h:70
Base integration class [abstract base].
Definition: inte.h:51
fp_t ters[nsub]
Squared error for subdivision.
fp_t xlo[nsub]
Lower end of subdivision.
inte< func_t, fp_t > * it
The base integration object.
int verbose
Verbosity.
Definition: inte.h:67
fp_t xhi[nsub]
High end of subdivision.
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.
std::string itos(int x)
Convert an integer to a string.
size_t nsubdiv
Number of subdivisions.

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