cern_adapt.h

/*
  -------------------------------------------------------------------
  
  Copyright (C) 2006, 2007, 2008, Andrew W. Steiner
  
  This file is part of O2scl.
  
  O2scl is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.
  
  O2scl is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with O2scl. If not, see <http://www.gnu.org/licenses/>.

  -------------------------------------------------------------------
*/
#ifndef O2SCL_CERN_ADAPT_H
#define O2SCL_CERN_ADAPT_H

#include <o2scl/inte.h>
#include <o2scl/cern_gauss56.h>
 
#ifndef DOXYGENP
namespace o2scl {
#endif

  /** 
      \brief Adaptive integration (CERNLIB)
    
      Uses a base integration object (default is \ref cern_gauss56) to
      perform adaptive integration by automatically subdividing the
      integration interval. At each step, the interval with the
      largest absolute uncertainty is divided in half. The routine
      stops if the absolute tolerance is less than \ref tolx, the
      relative tolerance is less than \ref tolf, or the number of
      segments exceeds the template parameter \c nsub (in which case
      the error handler is called, since the integration may not have
      been successful). The number of segments used in the last
      integration can be obtained from get_nsegments().

      The template parameter \c nsub, is the maximum number of
      subdivisions. It is automatically set to 100 in the original
      CERNLIB routine, and defaults to 100 here. The default base
      integration object is of type \ref cern_gauss56. This is the
      CERNLIB default, but can be modified by calling set_inte().
      
      This class is based on the CERNLIB routines RADAPT and
      DADAPT which are documented at
      http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d102/top.html
      
      \future Allow user to set the initial segments?
  */
  template<class param_t, class func_t=funct<param_t>, size_t nsub=100> 
    class cern_adapt : public inte<param_t,func_t> {

#ifndef DOXYGEN_INTERNAL

      protected:

      /// Lower end of subdivision
      double xlo[nsub];

      /// High end of subdivision
      double xhi[nsub];

      /// Value of integral for subdivision
      double tval[nsub];

      /// Squared error for subdivision
      double ters[nsub];
      
      /// Previous number of subdivisions
      int nter;

      /// Default integration object
      cern_gauss56<param_t,func_t> cg56;
      
      /// The base integration object
      inte<param_t, func_t> *it;
      
#endif
      
      public:
  
      cern_adapt() {
        nseg=1;
        nter=0;

        it=&cg56;
      }

      /// Set the base integration object to use
      int set_inte(inte<param_t,func_t> &i) {
        it=&i;
        return 0;
      }
      
      /** 
          \brief Number of subdivisions 
        
          The options are
          - 0: Use previous binning and do not subdivide further \n
          - 1: Automatic - adapt until tolerance is attained (default) \n
          - n: (n>1) split first in n equal segments, then adapt
          until tolerance is obtained.
      */
      size_t nseg;

      /// Return the number of segments used in the last integration
      size_t get_nsegments() {
        return nter;
      }

      /// Return the ith segment
      int get_ith_segment(size_t i, double &xlow, double &xhigh, double &value,
                          double &errsq) {
        if (i<nter) {
          xlow=xlo[i];
          xhigh=xhi[i];
          value=tval[i];
          errsq=ters[i];
        }
        return 0;
      }
      
      /// Return all of the segments
      template<class vec_t> 
      int get_segments(vec_t &xlow, vec_t &xhigh, vec_t &value, vec_t &errsq) {
        for(int i=0;i<nter;i++) {
          xlow[i]=xlo[i];
          xhigh[i]=xhi[i];
          value[i]=tval[i];
          errsq[i]=ters[i];
        }
        return 0;
      }

      /** \brief Integrate function \c func from \c a to \c b.
       */
      virtual double integ(func_t &func, double a, double b, param_t &pa) {
        double res;
        integ_err(func,a,b,pa,res,this->interror);
        return res;
      }

      /** \brief Integrate function \c func from \c a to \c b
          giving result \c res and error \c err
      */
      virtual int integ_err(func_t &func, double a, double b, param_t &pa,
                            double &res, double &err) 
      {
        
        double tvals=0.0, terss, xlob, xhib, yhib=0.0, te, root=0.0;
        int i, nsegd;

        this->last_conv=0;
  
        if (nseg==0) {
          if (nter==0) {
            // If the previous binning was requested, but
            // there is no previous binning stored, then
            // just shift to automatic binning
            nsegd=1;
          } else {
            tvals=0.0;
            terss=0.0;
            for(i=0;i<nter;i++) {
              it->integ_err(func,xlo[i],xhi[i],pa,tval[i],te);
              ters[i]=te*te;
              tvals+=tval[i];
              terss+=ters[i];
            }
            err=sqrt(2.0*terss);
            res=tvals;
            return 0;
          }
        }
  
        // Ensure we're not asking for too many divisions
        if (nsub<nseg) {
          nsegd=nsub;
        } else {
          nsegd=nseg;
        }

        // Compute the initial set of intervals and integral values
        xhib=a;
        double bin=(b-a)/((double)nsegd);
        for(i=0;i<nsegd;i++) {
          xlo[i]=xhib;
          xlob=xlo[i];
          xhi[i]=xhib+bin;
          if (i==nsegd-1) xhi[i]=b;
          xhib=xhi[i];
          it->integ_err(func,xlob,xhib,pa,tval[i],te);
          ters[i]=te*te;
        }
        nter=nsegd;

        for(size_t iter=1;iter<=nsub;iter++) {

          // Compute the total value of the integrand
          // and the squared uncertainty
          tvals=tval[0];
          terss=ters[0];
          for(i=1;i<nter;i++) {
            tvals+=tval[i];
            terss+=ters[i];
          }
          
          // Output iteration information
          if (this->verbose>0) {
            std::cout << "cern_adapt Iter: " << iter;
            std::cout.setf(std::ios::showpos);
            std::cout << " Res: " << tvals;
            std::cout.unsetf(std::ios::showpos);
            std::cout << " Err: " << sqrt(2.0*terss);
            if (this->tolx>this->tolf*fabs(tvals)) {
              std::cout << " Tol: " << this->tolx << std::endl;
            } else {
              std::cout << " Tol: " << this->tolf*fabs(tvals) << std::endl;
            }
            if (this->verbose>1) {
              char ch;
              std::cout << "Press a key and type enter to continue. " ;
              std::cin >> ch;
            }
          }

          // See if we're finished
          root=sqrt(2.0*terss);
          if (root<=this->tolx || root<=this->tolf*fabs(tvals)) {
            res=tvals;
            err=root;
            this->last_iter=iter;
            return 0;
          }

          // Test if we've run out of intervals
          if (nter==nsub) {
            res=tvals;
            err=root;
            this->last_iter=iter;
            std::string s="Reached maximum number of ";
            s+="subdivisions in cern_adapt::integ_err().";
            this->last_conv=gsl_etable;
            O2SCL_CONV_RET(s.c_str(),gsl_etable,this->err_nonconv);
          }

          // Find the segment with the largest error
          double bige=ters[0];
          int ibig=0;
          for(i=1;i<nter;i++) {
            if (ters[i]>bige) {
              bige=ters[i];
              ibig=i;
            }
          }

          // Subdivide that segment further
          xhi[nter]=xhi[ibig];
          double xnew=(xlo[ibig]+xhi[ibig])/2.0;
          xhi[ibig]=xnew;
          xlo[nter]=xnew;
          it->integ_err(func,xlo[ibig],xhi[ibig],pa,tval[ibig],te);
          ters[ibig]=te*te;
          it->integ_err(func,xlo[nter],xhi[nter],pa,tval[nter],te);
          ters[nter]=te*te;
          nter++;

        }

        // FIXME: Should we set an error here, or does this
        // only happen if we happen to need exactly nsub
        // intervals?
        res=tvals;
        err=root;
        return 0;
      }

    };

#ifndef DOXYGENP
}
#endif

#endif

---

Documentation generated with Doxygen and provided under the GNU Free Documentation License. See License Information for details.

Project hosting provided by , O2scl Sourceforge Project Page