Integrate.h

Go to the documentation of this file.

// -*- LSST-C++ -*-

/* 
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 * 
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program 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.
 * 
 * This program 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 LSST License Statement and 
 * the GNU General Public License along with this program.  If not, 
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */
 
#if !defined(LSST_AFW_MATH_INTEGRATE_H)
#define LSST_AFW_MATH_INTEGRATE_H 1

#include <functional>
#include <vector>
#include <queue>
#include <map>
#include <cmath>
#include <algorithm>
#include <assert.h>
#include <limits>
#include <ostream>
#include <sstream>
#include <complex>
#include <stdexcept>

#include "lsst/pex/exceptions.h"

#include "lsst/afw/math/IntGKPData10.h"

// == The following is Mike Jarvis original comment ==
//
// Basic Usage:
//
// First, define a function object, which should derive from 
// std::unary_function<double,double>.  For example, to integrate a
// Gaussian, use something along the lines of this:
//
// class Gauss :
//   public std::unary_function<double,double>
// {   
//   public :
// 
//     Gauss(double _mu, double _sig) :   
//       mu(_mu), sig(_sig), sigsq(_sig*_sig) {}
// 
//     double operator()(double x) const
//     {
//       const double SQRTTWOPI = 2.50662827463;
//       return exp(-pow(x-mu,2)/2./sigsq)/SQRTTWOPI/sig;
//     }
// 
//   private :
//     double mu,sig,sigsq;
// };
//
// Next, make an IntRegion object with the bounds of the integration region.
// You need to give it the type to use for the bounds and the value of the
// functions (which need to be the same currently - some day I'll allow 
// complex functions...).  
//
// For example, to integrate something from -1 to 1, use:
//
// integ::IntRegion<double> reg1(-1.,1.);
//
// (The integ:: is because everything here is in the integ namespace to
// help prevent name conflicts with other header files.)
//
// If a value is > 1.e10 or < -1.e10, then these values are taken to be
// infinity, rather than the actual value.  
// So to integrate from 0 to infinity:
//
// integ::IntRegion<double> reg2(0.,1.e100);
//
// Or, you can use the variable integ::MOCK_INF which might be clearer.
//
//
// Finally, to perform the integral, the line would be:
//
// double integ1 = int1d(Gauss(0.,1.),reg1,1.e-10,1.e-4);
// double integ2 = int1d(Gauss(0.,2.),reg2,1.e-10,1.e-4);
//
// which should yield 0.68 and 0.5 in our case.
//
// Those last two numbers indicate the precision required.  
// 1.e-10 is the required absolute error, and 
// 1.e-4 is the required relative error.
//
// If you want, you can omit these and 1.e-15,1.e-6 will be used as the 
// default precision (which are generally fine for most purposes).
//
// The absolute error only comes into play for results which are close to 
// 0 to prevent requiring an error of 0 for integrals which evaluate to 0 
// or very close to it.
//
//
//
// Advanced Usage:
//
// When an integration fails to converge with the usual GKP algorithm,
// it splits the region into 2 (or more) and tries again with each sub-region.
// The default is to just bisect the region (or something similarly smart for
// infinite regions), but if you know of a good place to split the region,
// you can tell it using:
//
// reg.AddSplit(10.)
//
// For example, if you know that you have a singularity somewhere in your
// region, it would help the program a lot to split there, so you 
// should add that as a split point.  Zeros can also be good choices.
//
// In addition to the integral being returned from int1d, int2d, or int3d as
// the return value, the value is also stored in the region itself. 
// You can access it using:
//
// reg.Area();
//
// There is also an estimate of the error in the value:
//
// reg.Err();
//
// (It is intended to be an overestimate of the actual error, 
// but it doesn't always get it completely right.)
//
//
// 
// Two- and Three-Dimensional Integrals:
//
// These are slightly more complicated.  The easiest case is when the
// bounds of the integral are a rectangle or 3d box.  In this case,
// you can still use the regular IntRegion.  The only new thing then
// is the definition of the function.  For example, to integrate 
// int(3x^2 + xy + y , x=0..1, y=0..1):
//
// struct Integrand :
//   public std::binary_function<double,double,double>
// {
//   double operator()(double x, double y) const
//   { return x*(3.*x + y) + y; }
// };
//
// integ::IntRegion<double> reg3(0.,1.);
// double integ3 = int2d(Integrand(),reg3,reg3);
//
// (Which should give 1.75 as the result.)
// 
//
//

namespace lsst {
namespace afw {
namespace math {

double const MOCK_INF = 1.e10;

#ifdef NDEBUG
#define integ_dbg1 if (false) (*_dbgout)
#define integ_dbg2 if (false) (*(reg.getDbgout()))
#define integ_dbg3 if (false) (*(tempreg.getDbgout()))
#else
#define integ_dbg1 if (_dbgout) (*_dbgout)
#define integ_dbg2 if (reg.getDbgout()) (*(reg.getDbgout()))
#define integ_dbg3 if (tempreg.getDbgout()) (*(tempreg.getDbgout()))
#endif

//#define COUNTFEVAL 
// If defined, then count the number of function evaluations
            
namespace details {            
template <class T> inline T norm(const T& x) { return x*x; }
using std::norm;
template <class T> inline T real(const T& x) { return x; }
using std::real;
#ifdef COUNTFEVAL
    int nfeval = 0;
#endif
}
            
            
template <class T>
struct IntRegion {

public:
    IntRegion(T const a, T const b, std::ostream* dbgout = 0) :
        _a(a), _b(b), _error(0.0), _area(0), _dbgout(dbgout) {}
    
    bool operator<(IntRegion<T> const &r2) const { return _error < r2._error; }
    bool operator>(IntRegion<T> const &r2) const { return _error > r2._error; }
    
    void SubDivide(std::vector<IntRegion<T> >* children) {
        assert(children->size() == 0);
        if (_splitpoints.size() == 0) { Bisect(); }
        if (_splitpoints.size() > 1) {
            std::sort(_splitpoints.begin(), _splitpoints.end());
        }
        
#if 0
        if (_a > _splitpoints[0] || _b < _splitpoints.back()) {
            std::cerr << "a, b = " << _a << ', ' << _b << std::endl;
            std::cerr << "_splitpoints = ";
            for (size_t i = 0; i<_splitpoints.size(); i++)  {
                std::cerr << _splitpoints[i] << "  ";
            }
            std::cerr << std::endl;
        }
#endif
        assert(_splitpoints[0] >= _a);
        assert(_splitpoints.back() <= _b);
        children->push_back(IntRegion<T>(_a, _splitpoints[0], _dbgout));
        for (size_t i = 1; i<_splitpoints.size(); i++) {
            children->push_back(IntRegion<T>(_splitpoints[i-1], _splitpoints[i], _dbgout));
        }
        children->push_back(IntRegion<T>(_splitpoints.back(), _b, _dbgout));
    }
    
    void Bisect() { _splitpoints.push_back((_a + _b)/2.0); }
    void AddSplit(const T x) { _splitpoints.push_back(x); }
    size_t NSplit() const { return _splitpoints.size(); }
    
    T const &Left() const { return _a; }
    T const &Right() const  {return _b; }
    T const &Err() const { return _error; }
    T const &Area() const { return _area; }
    void SetArea(const T& a, const T& e) {
        _area = a;
        _error = e;
    }

    std::ostream* getDbgout() { return _dbgout; }

private:
    T _a, _b, _error, _area;
    std::vector<T> _splitpoints;
    std::ostream* _dbgout;
    
};
            
            
double const DEFABSERR = 1.e-15;
double const DEFRELERR = 1.e-6;
 
namespace details {            
            
template <class T>
inline T Epsilon() { return std::numeric_limits<T>::epsilon(); }
template <class T>
inline T MinRep()  { return std::numeric_limits<T>::min(); }


#ifdef EXTRA_PREC_H
template <>
inline Quad Epsilon<Quad>() { return 3.08148791094e-33; }
template <>
inline Quad MinRep<Quad>()  { return 2.2250738585072014e-308; }
#endif

    
template <class T>
inline T rescale_error (T err, T const &resabs, T const &resasc) {
    if (resasc != 0.0 && err != 0.0) {
        T const scale = (200.0 * err / resasc);
        if (scale < 1.0) {
            err = resasc * scale * sqrt(scale);
        }
        else {
            err = resasc;
        }
    }
    if (resabs > MinRep<T>() / (50.0 * Epsilon<T>())) {
        T const min_err = 50.0 * Epsilon<T>() * resabs;
        if (min_err > err) {
            err = min_err;
        }
    }
    return err;
}

template <class UF>
inline bool intGKPNA(
                     UF const &func, IntRegion<typename UF::result_type>& reg,
                     typename UF::result_type const epsabs, 
                     typename UF::result_type const epsrel,
                     std::map<typename UF::result_type, typename UF::result_type>* fxmap = 0) {
    
    typedef typename UF::result_type UfResult;
    UfResult const a = reg.Left();
    UfResult const b = reg.Right();
                   
    UfResult const halfLength =  0.5 * (b - a);
    UfResult const absHalfLength = fabs (halfLength);
    UfResult const center = 0.5 * (b + a);
    UfResult const fCenter = func(center);
#ifdef COUNTFEVAL
    nfeval++;
#endif
    
    assert(gkp_wb<UfResult>(0).size() == gkp_x<UfResult>(0).size() + 1);
    UfResult area1 = gkp_wb<UfResult>(0).back() * fCenter;
    std::vector<UfResult> fv1, fv2;
    fv1.reserve(2*gkp_x<UfResult>(0).size() + 1);
    fv2.reserve(2*gkp_x<UfResult>(0).size() + 1);
    for (size_t k = 0; k<gkp_x<UfResult>(0).size(); k++) {
        UfResult const abscissa = halfLength * gkp_x<UfResult>(0)[k];
        UfResult const fval1 = func(center - abscissa);
        UfResult const fval2 = func(center + abscissa);
        area1 += gkp_wb<UfResult>(0)[k] * (fval1 + fval2);
        fv1.push_back(fval1);
        fv2.push_back(fval2);
        if (fxmap) {
            (*fxmap)[center - abscissa] = fval1;
            (*fxmap)[center + abscissa] = fval2;
        }
    }
#ifdef COUNTFEVAL
    nfeval += gkp_x<UfResult>(0).size()*2;
#endif
    
    integ_dbg2 << "level 0 rule: area = " << area1 << std::endl;
    
    UfResult err = 0; 
    bool calcabsasc = true;
    UfResult resabs = 0.0, resasc = 0.0;
    for (int level = 1; level < NGKPLEVELS; level++) {
        assert(gkp_wa<UfResult>(level).size() == fv1.size());
        assert(gkp_wa<UfResult>(level).size() == fv2.size());
        assert(gkp_wb<UfResult>(level).size() == gkp_x<UfResult>(level).size() + 1);
        UfResult area2 = gkp_wb<UfResult>(level).back() * fCenter;
        // resabs = approximation to integral of abs(f)
        if (calcabsasc) {
            resabs = fabs(area2);
        }
        for (size_t k = 0; k < fv1.size(); k++) {
            area2 += gkp_wa<UfResult>(level)[k] * (fv1[k] + fv2[k]);
            if (calcabsasc) {
                resabs += gkp_wa<UfResult>(level)[k] * (fabs(fv1[k]) + fabs(fv2[k]));
            }
        }
        for (size_t k = 0; k < gkp_x<UfResult>(level).size(); k++) {
            UfResult const abscissa = halfLength * gkp_x<UfResult>(level)[k];
            UfResult const fval1 = func(center - abscissa);
            UfResult const fval2 = func(center + abscissa);
            UfResult const fval = fval1 + fval2;
            area2 += gkp_wb<UfResult>(level)[k] * fval;
            if (calcabsasc) {
                resabs += gkp_wb<UfResult>(level)[k] * (fabs(fval1) + fabs(fval2));
            }
            fv1.push_back(fval1);
            fv2.push_back(fval2);
            if (fxmap) {
                (*fxmap)[center - abscissa] = fval1;
                (*fxmap)[center + abscissa] = fval2;
            }
        }
#ifdef COUNTFEVAL
        nfeval += gkp_x<UfResult>(level).size()*2;
#endif
        if (calcabsasc) {
            UfResult const mean = area1*UfResult(0.5);
            // resasc = approximation to the integral of abs(f-mean) 
            resasc = gkp_wb<UfResult>(level).back() * fabs(fCenter - mean);
            for (size_t k = 0; k<gkp_wa<UfResult>(level).size(); k++) {
                resasc += gkp_wa<UfResult>(level)[k] * (fabs(fv1[k] - mean) + fabs(fv2[k] - mean));
            }
            for (size_t k = 0; k<gkp_x<UfResult>(level).size(); k++) {
                resasc += gkp_wb<UfResult>(level)[k] * (fabs(fv1[k] - mean) + fabs(fv2[k] - mean));
            }
            resasc *= absHalfLength;
            resabs *= absHalfLength;
        }
        area2 *= halfLength;
        err = rescale_error (fabs(area2 - area1), resabs, resasc) ;
        if (err < resasc) {
            calcabsasc = false;
        }
        
        integ_dbg2 << "at level " << level << " area2 = " << area2;
        integ_dbg2 << " +- " << err << std::endl;
        
        //   test for convergence.
        if (err < epsabs || err < epsrel * fabs (area2)) {
            reg.SetArea(area2, err);
            return true;
        }
        area1 = area2;
    }
    
    // failed to converge 
    reg.SetArea(area1, err);
    
    integ_dbg2 << "Failed to reach tolerance with highest-order GKP rule";
    
    return false;
}


template <class UF>
inline void intGKP (
                    UF const &func, IntRegion<typename UF::result_type> &reg,
                    typename UF::result_type const epsabs,
                    typename UF::result_type const epsrel,
                    std::map<typename UF::result_type, typename UF::result_type>* fxmap = 0) {

    typedef typename UF::result_type UfResult;
    integ_dbg2 << "Start intGKP\n";
    
    assert(epsabs >= 0.0);
    assert(epsrel > 0.0);
    
    // perform the first integration 
    bool done = intGKPNA(func, reg, epsabs, epsrel, fxmap);
    if (done) return;
    
    integ_dbg2 << "In adaptive GKP, failed first pass... subdividing\n";
    integ_dbg2 << "Intial range = " << reg.Left() << ".." << reg.Right() << std::endl;
    
    int roundoffType1 = 0, errorType = 0;
    UfResult roundoffType2 = 0;
    size_t iteration = 1;
    
    std::priority_queue<IntRegion<UfResult>, std::vector<IntRegion<UfResult> > > allregions;
    allregions.push(reg);
    UfResult finalarea = reg.Area();
    UfResult finalerr = reg.Err();
    UfResult tolerance = std::max(epsabs, epsrel * fabs(finalarea));
    assert(finalerr > tolerance);
    
    while (!errorType && finalerr > tolerance) {
        // Bisect the subinterval with the largest error estimate 
        integ_dbg2 << "Current answer = " << finalarea << " +- " << finalerr;
        integ_dbg2 << "  (tol = " << tolerance << ")\n";
        IntRegion<UfResult> parent = allregions.top(); 
        allregions.pop();
        integ_dbg2 << "Subdividing largest error region ";
        integ_dbg2 << parent.Left() << ".." << parent.Right() << std::endl;
        integ_dbg2 << "parent area = " << parent.Area();
        integ_dbg2 << " +- " << parent.Err() << std::endl;
        std::vector<IntRegion<UfResult> > children;
        parent.SubDivide(&children);
        // For "GKP", there are only two, but for GKPOSC, there is one 
        // for each oscillation in region
        
        // Try to do at least 3x better with the children
        UfResult factor = 3*children.size()*finalerr/tolerance;
        UfResult newepsabs = fabs(parent.Err()/factor);
        UfResult newepsrel = newepsabs/fabs(parent.Area());
        integ_dbg2 << "New epsabs, rel = " << newepsabs << ", " << newepsrel;
        integ_dbg2 << "  (" << children.size() << " children)\n";
        
        UfResult newarea = UfResult(0.0);
        UfResult newerror = 0.0;
        for (size_t i = 0; i<children.size(); i++) {
            IntRegion<UfResult>& child = children[i];
            integ_dbg2 << "Integrating child " << child.Left();
            integ_dbg2 << ".." << child.Right() << std::endl;
            bool hasConverged;
            hasConverged = intGKPNA(func, child, newepsabs, newepsrel);
            integ_dbg2 << "child (" << i + 1 << '/' << children.size() << ") ";
            if (hasConverged) {
                integ_dbg2 << " converged.";
            } else {
                integ_dbg2 << " failed.";
            }
            integ_dbg2 << "  Area = " << child.Area() << " +- " << child.Err() << std::endl;
            
            newarea += child.Area();
            newerror += child.Err();
        }
        integ_dbg2 << "Compare: newerr = " << newerror;
        integ_dbg2 << " to parent err = " << parent.Err() << std::endl;
        
        finalerr += (newerror - parent.Err());
        finalarea += newarea - parent.Area();
        
        UfResult delta = parent.Area() - newarea;
        if (newerror <= parent.Err() && fabs (delta) <=  parent.Err()
            && newerror >= 0.99 * parent.Err()) {
            integ_dbg2 << "roundoff type 1: delta/newarea = ";
            integ_dbg2 << fabs(delta)/fabs(newarea);
            integ_dbg2 << ", newerror/error = " << newerror/parent.Err() << std::endl;
            roundoffType1++;
        }
        if (iteration >= 10 && newerror > parent.Err() && 
            fabs(delta) <= newerror - parent.Err()) {
            integ_dbg2 << "roundoff type 2: newerror/error = ";
            integ_dbg2 << newerror/parent.Err() << std::endl;
            roundoffType2 += std::min(newerror/parent.Err() - 1.0, UfResult(1.0));
        }
        
        tolerance = std::max(epsabs, epsrel * fabs(finalarea));
        if (finalerr > tolerance) {
            if (roundoffType1 >= 200) {
                errorType = 1; // round off error 
                integ_dbg2 << "GKP: Round off error 1\n";
            }
            if (roundoffType2 >= 200.0) {
                errorType = 2; // round off error 
                integ_dbg2 << "GKP: Round off error 2\n";
            }
            if (fabs((parent.Right() - parent.Left())/(reg.Right() - reg.Left())) 
                < Epsilon<double>()) {
                errorType = 3; // found singularity
                integ_dbg2 << "GKP: Probable singularity\n";
            }
        }
        for (size_t i = 0; i<children.size(); i++) {
            allregions.push(children[i]);
        }
        iteration++;
    } 
    
    // Recalculate finalarea in case there are any slight rounding errors
    finalarea = 0.0;
    finalerr = 0.0;
    while (!allregions.empty()) {
        IntRegion<UfResult> const &r = allregions.top();
        finalarea += r.Area();
        finalerr += r.Err();
        allregions.pop();
    }
    reg.SetArea(finalarea, finalerr);
    
    if (errorType == 1) {
        std::ostringstream s;
        s << "Type 1 roundoff's = " << roundoffType1;
        s << ", Type 2 = " << roundoffType2 << std::endl;
        s << "Roundoff error 1 prevents tolerance from being achieved ";
        s << "in intGKP\n";
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, s.str());
    } else if (errorType == 2) {
        std::ostringstream s;
        s << "Type 1 roundoff's = " << roundoffType1;
        s << ", Type 2 = " << roundoffType2 << std::endl;
        s << "Roundoff error 2 prevents tolerance from being achieved ";
        s << "in intGKP\n";
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, s.str());
    } else if (errorType == 3) {
        std::ostringstream s;
        s << "Bad integrand behavior found in the integration interval ";
        s << "in intGKP\n";
        throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, s.str());
    }
}


template <class UF>
struct AuxFunc1 : // f(1/x-1) for int(a..infinity)
        public std::unary_function<typename UF::argument_type, typename UF::result_type> {
public:
    AuxFunc1(const UF& f) : _f(f) {}
    typename UF::result_type operator()(typename UF::argument_type x) const {
        return _f(1.0/x - 1.0)/(x*x);
    }
private:
    UF const &_f;
};

    
template <class UF>
AuxFunc1<UF> inline Aux1(UF uf) { return AuxFunc1<UF>(uf); }

    
template <class UF>
struct AuxFunc2 : // f(1/x+1) for int(-infinity..b)
        public std::unary_function<typename UF::argument_type, typename UF::result_type> {
public:
    AuxFunc2(UF const &f) : _f(f) {}
    typename UF::result_type operator()(typename UF::argument_type x) const {
        return _f(1.0/x + 1.0)/(x*x);
    }
private:
    UF const &_f;
};


template <class UF>
AuxFunc2<UF> inline Aux2(UF uf) { return AuxFunc2<UF>(uf); }



template <class T>
struct ConstantReg1 : public std::unary_function<T, IntRegion<T> > {
    ConstantReg1(T a, T b) : ir(a, b) {}
    ConstantReg1(IntRegion<T> const &r) : ir(r) {}
    IntRegion<T> operator()(T) const { return ir; }
    IntRegion<T> ir;
};
    
template <class T>
struct ConstantReg2 : public std::binary_function<T, T, IntRegion<T> > {
    ConstantReg2(T a, T b) : ir(a, b) {}
    ConstantReg2(IntRegion<T> const &r) : ir(r) {}
    IntRegion<T> operator()(T x, T y) const { return ir; }
    IntRegion<T> ir;
};


// pulled from MoreFunctional.h.  Needed in class Int2DAuxType and Int3DAuxType
template <class BF>
class binder2_1 
    : public std::unary_function<typename BF::second_argument_type,
                                 typename BF::result_type> {
public:
    binder2_1(const BF& oper,
              typename BF::first_argument_type val)
        : _oper(oper), _value(val) {}
    typename BF::result_type 
    operator()(const typename BF::second_argument_type& x) const {
        return _oper(_value, x); 
    }
protected:
    BF _oper;
    typename BF::first_argument_type _value;
};
    
template <class BF, class Tp>
inline binder2_1<BF> bind21(const BF& oper, const Tp& x) {
    typedef typename BF::first_argument_type Arg;
    return binder2_1<BF>(oper, static_cast<Arg>(x));
}

    
template <class BF, class YREG>
class Int2DAuxType : public std::unary_function<typename BF::first_argument_type, typename BF::result_type> {
public:
    Int2DAuxType(BF const &func, YREG const &yreg,
                 typename BF::result_type  const &abserr,
                 typename BF::result_type  const &relerr) :
        _func(func), _yreg(yreg), _abserr(abserr), _relerr(relerr) {}
    
    typename BF::result_type operator()(typename BF::first_argument_type x) const  {
        typename YREG::result_type tempreg = _yreg(x);
        typename BF::result_type result = 
            int1d(bind21(_func, x), tempreg, _abserr, _relerr);
        integ_dbg3 << "Evaluated int2dAux at x = " << x;
        integ_dbg3 << ": f = " << result << " +- " << tempreg.Err() << std::endl;
        return result;
    }
    
private:
    BF   const &_func;
    YREG const &_yreg;
    typename BF::result_type _abserr, _relerr;
};
    
    
// pulled from MoreFunctional.h.  Needed in class Int3DAuxtype    
template <class TF>
class binder3_1 
    : public std::binary_function<typename TF::secondof3_argument_type,
                                  typename TF::thirdof3_argument_type,
                                  typename TF::result_type> {
public:
    binder3_1(const TF& oper,
              typename TF::firstof3_argument_type val)
        : _oper(oper), _value(val) {}
    typename TF::result_type 
    operator()(typename TF::secondof3_argument_type const &x1, 
               typename TF::thirdof3_argument_type const &x2) const {
        return _oper(_value, x1, x2); 
    }
protected:
    TF _oper;
    typename TF::firstof3_argument_type _value;
};

template <class TF, class Tp>
inline binder3_1<TF>
bind31(const TF& oper, const Tp& x) {
    typedef typename TF::firstof3_argument_type Arg;
    return binder3_1<TF>(oper, static_cast<Arg>(x));
}


template <class TF, class YREG, class ZREG>
class Int3DAuxType :
        public std::unary_function<typename TF::firstof3_argument_type, typename TF::result_type> {
public:
    Int3DAuxType(const TF& func, const YREG& yreg, const ZREG& zreg, 
                 const typename TF::result_type& abserr,
                 const typename TF::result_type& relerr) :
        _func(func), _yreg(yreg), _zreg(zreg), _abserr(abserr), _relerr(relerr) {}
    
    typename TF::result_type operator()(typename TF::firstof3_argument_type x) const {
        typename YREG::result_type tempreg = _yreg(x);
        typename TF::result_type result = 
            int2d(bind31(_func, x), tempreg, bind21(_zreg, x), _abserr, _relerr);
        integ_dbg3 << "Evaluated int3dAux at x = " << x;
        integ_dbg3 << ": f = " << result << " +- " << tempreg.Err() << std::endl;
        return result;
    }
    
private:
    const TF& _func;
    const YREG& _yreg;
    const ZREG& _zreg;
    typename TF::result_type _abserr, _relerr;
};
    
} // end namespace details



template <class UF>
inline typename UF::result_type int1d(
                                      UF const &func, 
                                      IntRegion<typename UF::result_type>& reg,
                                      typename UF::result_type const &abserr = DEFABSERR,
                                      typename UF::result_type const &relerr = DEFRELERR) {
    
    typedef typename UF::result_type UfResult;
    using namespace details;
    
    integ_dbg2 << "start int1d: " << reg.Left() << ".." << reg.Right() << std::endl;
    
    if ((reg.Left() <= -MOCK_INF && reg.Right() > 0) || (reg.Right() >= MOCK_INF && reg.Left() < 0)) { 
        reg.AddSplit(0);
    }

    if (reg.NSplit() > 0) {
        std::vector<IntRegion<UfResult> > children;
        reg.SubDivide(&children);
        integ_dbg2 << "Subdivided into " << children.size() << " children\n";
        UfResult answer = UfResult();
        UfResult err = 0; 
        for (size_t i = 0; i<children.size(); i++) {
            IntRegion<UfResult>& child = children[i];
            integ_dbg2 << "i = " << i;
            integ_dbg2 << ": bounds = " << child.Left() << ", " << child.Right() << std::endl;
            answer += int1d(func, child, abserr, relerr);
            err += child.Err();
            integ_dbg2 << "subint = " << child.Area() << " +- " << child.Err() << std::endl;
        }
        reg.SetArea(answer, err);
        return answer;
        
    }  else {
        if (reg.Left() <= -MOCK_INF) {
            integ_dbg2 << "left = -infinity, right = " << reg.Right() << std::endl;
            assert(reg.Right() <= 0.0);
            IntRegion<UfResult> modreg(1.0/(reg.Right() - 1.0), 0.0, reg.getDbgout());
            intGKP(Aux2<UF>(func), modreg, abserr, relerr);
            reg.SetArea(modreg.Area(), modreg.Err());
        } else if (reg.Right() >= MOCK_INF) {
            integ_dbg2 << "left = " << reg.Left() << ", right = infinity\n";
            assert(reg.Left() >= 0.0);
            IntRegion<UfResult> modreg(0.0, 1.0/(reg.Left() + 1.0), reg.getDbgout());
            intGKP(Aux1<UF>(func), modreg, abserr, relerr);
            reg.SetArea(modreg.Area(), modreg.Err());
        } else {
            integ_dbg2 << "left = " << reg.Left();
            integ_dbg2 << ", right = " << reg.Right() << std::endl;
            intGKP(func, reg, abserr, relerr);
        }
        integ_dbg2 << "done int1d  answer = " << reg.Area();
        integ_dbg2 << " +- " << reg.Err() << std::endl;
        return reg.Area();
    }
}

            
template <class BF, class YREG>
inline typename BF::result_type int2d(
                                      BF const &func,
                                      IntRegion<typename BF::result_type> &reg,
                                      YREG const &yreg,
                                      typename BF::result_type const &abserr = DEFABSERR,
                                      typename BF::result_type const &relerr = DEFRELERR) {
    
    using namespace details;
    integ_dbg2 << "Starting int2d: range = ";
    integ_dbg2 << reg.Left() << ".." << reg.Right() << std::endl;
    Int2DAuxType<BF, YREG> faux(func, yreg, abserr*1.0e-3, relerr*1.0e-3);
    typename BF::result_type answer = int1d(faux, reg, abserr, relerr);
    integ_dbg2 << "done int2d  answer = " << answer << " +- " << reg.Err() << std::endl;
    return answer;
}

template <class TF, class YREG, class ZREG> 
inline typename TF::result_type int3d(
                                      TF const &func,
                                      IntRegion<typename TF::result_type>& reg,
                                      YREG const &yreg,
                                      ZREG const &zreg,
                                      typename TF::result_type const &abserr = DEFABSERR,
                                      typename TF::result_type const &relerr = DEFRELERR) {
    
    using namespace details;
    integ_dbg2 << "Starting int3d: range = ";
    integ_dbg2 << reg.Left() << ".." << reg.Right() << std::endl;
    Int3DAuxType<TF, YREG, ZREG> faux(func, yreg, zreg, abserr*1.e-3, relerr*1.e-3);
    typename TF::result_type answer = int1d(faux, reg, abserr, relerr);
    integ_dbg2 << "done int3d  answer = " << answer << " +- " << reg.Err() << std::endl;
    return answer;
}

template <class BF>
inline typename BF::result_type int2d(
                                      BF const &func,
                                      IntRegion<typename BF::result_type> &reg,
                                      IntRegion<typename BF::result_type> &yreg,
                                      typename BF::result_type const &abserr = DEFABSERR,
                                      typename BF::result_type const &relerr = DEFRELERR) {
    using namespace details;
    return int2d(func, reg, ConstantReg1<typename BF::result_type>(yreg), abserr, relerr); 
}

template <class TF>
inline typename TF::result_type int3d(
                                      TF const &func,
                                      IntRegion<typename TF::result_type> &reg,
                                      IntRegion<typename TF::result_type> &yreg,
                                      IntRegion<typename TF::result_type> &zreg,
                                      typename TF::result_type const &abserr = DEFABSERR,
                                      typename TF::result_type const &relerr = DEFRELERR) {
    using namespace details;
    return int3d(func, reg, ConstantReg1<typename TF::result_type>(yreg), 
                 ConstantReg2<typename TF::result_type>(zreg), abserr, relerr);
}
            




// =============================================================
template<typename UnaryFunctionT>
typename UnaryFunctionT::result_type integrate(UnaryFunctionT func,
                                               typename UnaryFunctionT::argument_type const a,
                                               typename UnaryFunctionT::argument_type const b,
                                               double eps = 1.0e-6)  {
    
    typedef typename UnaryFunctionT::argument_type Arg;
    IntRegion<Arg> region(a, b);
    
    return int1d(func, region, DEFABSERR, eps);
}
            


namespace details {            

template<typename BinaryFunctionT>
class FunctionWrapper :
        public std::unary_function<typename BinaryFunctionT::second_argument_type,
                                   typename BinaryFunctionT::result_type> {
public:
    FunctionWrapper(BinaryFunctionT func,
                    typename BinaryFunctionT::first_argument_type const x1,
                    typename BinaryFunctionT::first_argument_type const x2,
                    double const eps = 1.0e-6) :
        _func(func), _x1(x1), _x2(x2), _eps(eps) {}
    typename BinaryFunctionT::result_type operator()(typename
                                                     BinaryFunctionT::second_argument_type const y) const {
        return integrate(std::bind2nd(_func, y), _x1, _x2, _eps);
    }
private:
    BinaryFunctionT _func;
    typename BinaryFunctionT::first_argument_type _x1, _x2;
    double _eps;
};

} // end of namespace afw::math::details



// =============================================================
template<typename BinaryFunctionT>
typename BinaryFunctionT::result_type integrate2d(BinaryFunctionT func,
                                                  typename BinaryFunctionT::first_argument_type const x1,
                                                  typename BinaryFunctionT::first_argument_type const x2,
                                                  typename BinaryFunctionT::second_argument_type const y1,
                                                  typename BinaryFunctionT::second_argument_type const y2,
                                                  double eps = 1.0e-6) {
    using namespace details;
    // note the more stringent eps requirement to ensure the requested limit
    // can be reached.
    FunctionWrapper<BinaryFunctionT> fwrap(func, x1, x2, eps);
    return integrate(fwrap, y1, y2, eps);
}

            
}}} // end namespaces lsst/afw/math

#endif