Integrator.h

/**
 * \file Integrator.h
 * This file contains routines for numerical integration using Gauss-Kronrod quadrature.
 * \sa R Piessens, E de Doncker-Kapenger, C Ueberhuber, D Kahaner, QUADPACK, A Subroutine Package
 *     for Automatic Integration, Springer Verlag, 1983.
 */

#ifndef EIGEN_INTEGRATOR_H
#define EIGEN_INTEGRATOR_H

namespace Eigen
{
    /**
     * \ingroup Quadrature_Module
     *
     * \brief This class provides numerical integration functionality.
     *
     * Memory management and additional information (e.g. number of function calls, estimated error),
     * are provided by the class in a way that can support the future porting of additional quadrature
     * functions from the QUADPACK library.
     *
     * \todo Ensure only appropriates types are used for Scalar, e.g. prohibit integers.
     */
    template <typename Scalar>

    class Integrator
    {
    public:

        /**
         * \brief The local Gauss-Kronrod quadrature rule to use.
         */
        enum QuadratureRule
        {
            GaussKronrod15  =  1,  /**<  Use   7,  15 points. */
            GaussKronrod21  =  2,  /**<  Use  10,  21 points. */
            GaussKronrod31  =  3,  /**<  Use  15,  31 points. */
            GaussKronrod41  =  4,  /**<  Use  20,  41 points. */
            GaussKronrod51  =  5,  /**<  Use  25,  51 points. */
            GaussKronrod61  =  6,  /**<  Use  30,  61 points. */
            GaussKronrod71  =  7,  /**<  Use  35,  71 points. */
            GaussKronrod81  =  8,  /**<  Use  40,  81 points. */
            GaussKronrod91  =  9,  /**<  Use  45,  91 points. */
            GaussKronrod101 = 10,  /**<  Use  50, 101 points. */
            GaussKronrod121 = 11,  /**<  Use  60, 121 points. */
            GaussKronrod201 = 12   /**<  Use 100, 201 points. */
        };

        /**
         * \brief Prepares an Integrator for a call to a quadrature function.
         *
         * \param[in] maxSubintervals The maximum number of subintervals allowed in the subdivision process
         *        of quadrature functions. This corresponds to the amount of memory allocated for said
         *        functions.
         */
        Integrator(const int maxSubintervals)
            : m_maxSubintervals(maxSubintervals)
        {
            assert(maxSubintervals >= 1); // \todo use Eigen assert.

            m_errorListIndices.resize(maxSubintervals, 1);
            m_lowerList.resize(maxSubintervals, 1);
            m_upperList.resize(maxSubintervals, 1);
            m_integralList.resize(maxSubintervals, 1);
            m_errorList.resize(maxSubintervals, 1);
        }

        /**
         * \brief This function calculates an approximation I' to a given definite integral I, the
         *        integral of f from lowerLimit to upperLimit, hopefully satisfying
         *        abs(I - I') <= max(desiredAbsoluteError, desiredRelativeError * abs(I)).
         *
         * This function is best suited for integrands without singularities or discontinuities, which
         * are too difficult for non-adaptive quadrature, and, in particular, for integrands with
         * oscillating behavior of a non-specific type.
         *
         * \param[in,out] functionType The function type defining the integrand function.
         * \param[in] lowerLimit The lower limit of integration.
         * \param[in] upperLimit The upper limit of integration.
         * \param[in] desiredAbsoluteError The absolute accuracy requested.
         * \param[in] desiredRelativeError The relative accuracy requested.
         *            If desiredAbsoluteError <= 0 and desiredRelativeError < 50 * machinePrecision,
         *            the routine will end with errorCode = 6.
         * \param[in] quadratureRule The local Gauss-Kronrod quadrature rule to use.
         *
         * \returns The approximation to the integral.
         */
        template <typename FunctionType>
        Scalar quadratureAdaptive(const FunctionType& functionType,
                                  const Scalar lowerLimit,
                                  const Scalar upperLimit,
                                  const Scalar desiredAbsoluteError = Scalar(0.),
                                  const Scalar desiredRelativeError = Scalar(0.),
                                  const QuadratureRule quadratureRule = Eigen::Integrator<Scalar>::QuadratureRule(1))
        {
            using std::abs;
            using std::max;

            if ((desiredAbsoluteError <= Scalar(0.) &&
                desiredRelativeError < NumTraits<Scalar>::epsilon()) ||
                m_maxSubintervals < 1)
            {
                m_errorCode = 6;
                return Scalar(0.);
            }

            m_errorCode = 0;
            m_numEvaluations = 0;
            m_lowerList[0] = lowerLimit;
            m_upperList[0] = upperLimit;
            m_integralList[0] = Scalar(0.);
            m_errorList[0] = Scalar(0.);
            m_errorListIndices[0] = 0;
            m_errorListIndices[1] = 1;

            Scalar absDiff = 0.;
            Scalar absResult = 0.;

            // First approximation to the integral
            Scalar integral = quadratureKronrod(functionType, lowerLimit, upperLimit, m_estimatedError, absDiff, absResult, quadratureRule);

            m_numSubintervals = 1;
            m_integralList[0] = integral;
            m_errorList[0] = m_estimatedError;

            // Test on accuracy.
            Scalar errorBound = max(desiredAbsoluteError, desiredRelativeError * abs(integral));

            if (m_maxSubintervals == 1)
            {
                m_errorCode = 1;
            }
            else if (m_estimatedError <= NumTraits<Scalar>::epsilon() * Scalar(50.) * absDiff
                     && m_estimatedError > errorBound)
            {
                m_errorCode = 2;
            }

            if (m_errorCode != 0 ||
                m_estimatedError == Scalar(0.) ||
                (m_estimatedError <= errorBound &&
                 m_estimatedError != absResult))
            {
                if (quadratureRule == GaussKronrod15)
                {
                    m_numEvaluations = m_numEvaluations * 30 + 15;
                }
                else
                {
                    m_numEvaluations = (quadratureRule * 10 + 1) * (m_numEvaluations * 2 + 1);
                }

                return integral;
            }

            // The sum of the integrals over the subintervals.
            Scalar area = integral;

            Scalar errorSum = m_estimatedError;

            // The maximum interval error.
            Scalar errorMax = m_estimatedError;

            // An index into m_errorList at the interval with largest error estimate.
            Index maxErrorIndex = 0;
            Index nrMax = 0;
            
            int roundOff1 = 0;
            int roundOff2 = 0;

            Scalar error1 = 0.;
            Scalar error2 = 0.;
            Scalar absDiff1 = 0.;
            Scalar absDiff2 = 0.;

            // Main loop for the integration
            for (m_numSubintervals = 2; m_numSubintervals <= m_maxSubintervals; ++m_numSubintervals)
            {
                const Index numSubintervalsIndex = m_numSubintervals - 1;

                // Bisect the subinterval with the largest error estimate.
                const Scalar lower1 = m_lowerList[maxErrorIndex];
                const Scalar upper2 = m_upperList[maxErrorIndex];

                const Scalar upper1 = (lower1 + upper2) * Scalar(.5);
                const Scalar lower2 = upper1;

                const Scalar area1 = quadratureKronrod(functionType, lower1, upper1, error1, absResult, absDiff1, quadratureRule);
                const Scalar area2 = quadratureKronrod(functionType, lower2, upper2, error2, absResult, absDiff2, quadratureRule);

                // Improve previous approximations to integral and error and test for accuracy.
                ++(m_numEvaluations);
                
                const Scalar area12 = area1 + area2;
                const Scalar error12 = error1 + error2;
                
                errorSum += error12 - errorMax;
                area += area12 - m_integralList[maxErrorIndex];

                if (absDiff1 != error1 &&
                    absDiff2 != error2)
                {
                    if (abs(m_integralList[maxErrorIndex] - area12) <= abs(area12) * Scalar(1.e-5) &&
                        error12 >= errorMax * Scalar(.99))
                    {
                        ++roundOff1;
                    }
                    
                    if (m_numSubintervals > 10 &&
                        error12 > errorMax)
                    {
                        ++roundOff2;
                    }
                }

                m_integralList[maxErrorIndex] = area1;
                m_integralList[numSubintervalsIndex] = area2;

                errorBound = max(desiredAbsoluteError, desiredRelativeError * abs(area));

                if (errorSum > errorBound)
                {
                    // Set error flag in the case that the number of subintervals has reached the max allowable.
                    if (m_numSubintervals == m_maxSubintervals)
                    {
                        m_errorCode = 1;
                    }
                    // Test for roundoff error and set error flag.
                    else if (roundOff1 >= 6 || roundOff2 >= 20)
                    {
                        m_errorCode = 2;
                    }
                    // Set m_error_code in the case of poor integrand behaviour within the integration range.
                    else if (max(abs(lower1), abs(upper2)) <=
                             (NumTraits<Scalar>::epsilon() * Scalar(100.) + Scalar(1.))
                             * (abs(lower2) + (std::numeric_limits<Scalar>::min)() * Scalar(1.e3)))
                    {
                        m_errorCode = 3;
                    }
                }

                // Append the newly-created intervals to the list.
                if (error2 > error1)
                {
                    m_lowerList[numSubintervalsIndex] = lower1;
                    m_lowerList[maxErrorIndex] = lower2;
                    m_upperList[numSubintervalsIndex] = upper1;
                    m_integralList[maxErrorIndex] = area2;
                    m_integralList[numSubintervalsIndex] = area1;
                    m_errorList[maxErrorIndex] = error2;
                    m_errorList[numSubintervalsIndex] = error1;
                }
                else
                {
                    m_lowerList[numSubintervalsIndex] = lower2;
                    m_upperList[maxErrorIndex] = upper1;
                    m_upperList[numSubintervalsIndex] = upper2;
                    m_errorList[maxErrorIndex] = error1;
                    m_errorList[numSubintervalsIndex] = error2;
                }

                // Maintain the descending ordering in the list of error estimates and select the subinterval
                // with the largest error estimate, (the next subinterval to be bisected).
                quadratureSort(maxErrorIndex, errorMax, nrMax);

                if (m_errorCode != 0 ||
                    errorSum <= errorBound ||
                    m_numSubintervals == m_maxSubintervals)
                {
                    break;
                }
            }

            integral = Scalar(0.);

            for (Index k = 0; k < m_numSubintervals; ++k)
            {
                integral += m_integralList[k];
            }

            m_estimatedError = errorSum;

            if (quadratureRule == GaussKronrod15)
            {
                m_numEvaluations = m_numEvaluations * 30 + 15;
            }
            else
            {
                m_numEvaluations = (quadratureRule * 10 + 1) * (m_numEvaluations * 2 + 1);
            }

            return integral;
        }

        /**
         * \brief Returns the estimated absolute error from the last integration.
         *
         * \returns The value returned will only be valid after calling quadratureAdaptive at least once.
         */
        inline Scalar estimatedError()
        {
            return m_estimatedError;
        }

        /**
         * \brief Returns the error code.
         *
         * \returns The value returned will only be valid after calling quadratureAdaptive at least once.
         */
        inline Index errorCode()
        {
            return m_errorCode;
        }

    private:

        /**
         * \brief This routine maintains the descending ordering in the list of the local error
         *        estimates resulting from the interval subdivision process.
         *
         * At each call two error estimates are inserted using the sequential search method, top-down
         * for the largest error estimate and bottom-up for the smallest error estimate.
         *
         * \param[in,out] maxErrorIndex The index to the nrMax-th largest error estimate currently in the list.
         * \param[in,out] errorMax The nrMax-th largest error estimate. errorMaxIndex = errorList(maxError).
         * \param[in,out] nrMax The integer value such that maxError = errorListIndices(nrMax).
         */
        void quadratureSort(Index& maxErrorIndex,
                            Scalar& errorMax,
                            Index& nrMax)
        {
            if (m_numSubintervals <= 2)
            {
                m_errorListIndices[0] = 0;
                m_errorListIndices[1] = 1;
                maxErrorIndex = m_errorListIndices[nrMax];
                errorMax = m_errorList[maxErrorIndex];
                return;
            }

            // This part of the routine is only executed if, due to a difficult integrand, subdivision has
            // increased the error estimate. In the normal case the insert procedure should start after the
            // nrMax-th largest error estimate.
            Index i = 0;
            Index succeed = 0;
            const Scalar errorMaximum = m_errorList[maxErrorIndex];

            if (nrMax != 1)
            {
                for (i = 1; i < nrMax; ++i)
                {
                    succeed = m_errorListIndices[nrMax - 1];

                    if (errorMaximum <= m_errorList[succeed])
                    {
                        break;
                    }

                    m_errorListIndices[nrMax] = succeed;
                    --nrMax;
                }
            }

            // Compute the number of elements in the list to be maintained in descending order. This number
            // depends on the number of subdivisions remaining allowed.
            Index topBegin = m_numSubintervals - 1;
            Index bottomEnd = topBegin - 1;
            Index start = nrMax + 1;

            if (m_numSubintervals > m_maxSubintervals / 2 + 2)
            {
                topBegin = m_maxSubintervals + 3 - m_numSubintervals + 1;
            }

            // Insert errorMax by traversing the list top-down, starting comparison from the element
            // errorlist(m_errorListIndices(nrMax+1)).
            if (start <= bottomEnd)
            {
                for (i = start; i <= bottomEnd; ++i)
                {
                    succeed = m_errorListIndices[i];

                    if (errorMaximum >= m_errorList[succeed])
                    {
                        break;
                    }

                    m_errorListIndices[i - 1] = succeed;
                }
            }

            if (start > bottomEnd)
            {
                m_errorListIndices[bottomEnd] = maxErrorIndex;
                m_errorListIndices[topBegin] = m_numSubintervals - 1;
                maxErrorIndex = m_errorListIndices[nrMax];
                errorMax = m_errorList[maxErrorIndex];
                return;
            }

            // Insert errorMin by traversing the list bottom-up.
            m_errorListIndices[i - 1] = maxErrorIndex;

            Index tempIndex = bottomEnd;
            
            for (Index j = i; j <= bottomEnd; ++j)
            {
                succeed = m_errorListIndices[tempIndex];

                if (m_errorList[m_numSubintervals - 1] < m_errorList[succeed])
                {
                    m_errorListIndices[tempIndex + 1] = m_numSubintervals - 1;
                    maxErrorIndex = m_errorListIndices[nrMax];
                    errorMax = m_errorList[maxErrorIndex];
                    return;
                }

                m_errorListIndices[tempIndex + 1] = succeed;
                --tempIndex;
            }

            m_errorListIndices[i] = m_numSubintervals - 1;
            maxErrorIndex = m_errorListIndices[nrMax];
            errorMax = m_errorList[maxErrorIndex];
            return;
        }

        /**
         * \brief This function calculates an approximation I' to a given definite integral I, the
         *        integral of f from lowerLimit to upperLimit and provides an error estimate.
         *
         * \param[in] f The variable representing the function f(x) be integrated.
         * \param[in] lowerLimit The lower limit of integration.
         * \param[in] upperLimit The upper limit of integration.
         * \param[in,out] errorEstimate Estimate of the modulus of the absolute error, not to exceed
         *             abs(I - I').
         * \param[in,out] absIntegral The approximation to the integral of abs(f) from lowerLimit to
         *             upperLimit.
         * \param[in,out] absDiffIntegral The approximation to the integral of
         *             abs(f - I/(upperLimit - lowerLimit)).
         *
         * \returns The approximation I' to the integral I. It is computed by applying the 15, 21, 31,
         *          41, 51, 61, 71, 81, 91, 101, 121, 201-point kronrod rule obtained by optimal addition
         *          of abscissae to the 7, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 100-point Gauss rule.
         *
         * \detail This series of functions represents a priority queue data structure:
         *         - Apply Gauss-Kronrod on the whole initial interval and estimate the error.
         *         - Split the interval symmetrically in two and again apply Gauss-Kronrod on the intervals,
         *           estimate the errors and add a pair (or tuple) "interval, error" to the priority queue.
         *         - If the total error is smaller than the requested tolerance, or if the maximum number
         *           of subdivisions is reached, discontinue the process.
         *         - Otherwise, pop the top element from the queue (highest error), split the interval in two,
         *           and apply Gauss-Kronrod to the new intervals, add those two new elements to the priority
         *           queue and repeat from the previous step.
         */
        template <typename FunctionType>
        Scalar quadratureKronrod(
            const FunctionType& functionType, const Scalar lowerLimit, const Scalar upperLimit,
            Scalar& estimatedError, Scalar& absIntegral, Scalar& absDiffIntegral,
            const QuadratureRule quadratureRule)
        {
            switch (quadratureRule)
            {
            case GaussKronrod15:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod15,
                QuadratureKronrod<Scalar>::weightsGaussKronrod15, QuadratureKronrod<Scalar>::weightsGauss15,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod21:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod21,
                QuadratureKronrod<Scalar>::weightsGaussKronrod21, QuadratureKronrod<Scalar>::weightsGauss21,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod31:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod31,
                QuadratureKronrod<Scalar>::weightsGaussKronrod31, QuadratureKronrod<Scalar>::weightsGauss31,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod41:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod41,
                QuadratureKronrod<Scalar>::weightsGaussKronrod41, QuadratureKronrod<Scalar>::weightsGauss41,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod51:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod51,
                QuadratureKronrod<Scalar>::weightsGaussKronrod51, QuadratureKronrod<Scalar>::weightsGauss51,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod61:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod61,
                QuadratureKronrod<Scalar>::weightsGaussKronrod61, QuadratureKronrod<Scalar>::weightsGauss61,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod71:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod71,
                QuadratureKronrod<Scalar>::weightsGaussKronrod71, QuadratureKronrod<Scalar>::weightsGauss71,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod81:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod81,
                QuadratureKronrod<Scalar>::weightsGaussKronrod81, QuadratureKronrod<Scalar>::weightsGauss81,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod91:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod91,
                QuadratureKronrod<Scalar>::weightsGaussKronrod91, QuadratureKronrod<Scalar>::weightsGauss91,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod101:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod101,
                QuadratureKronrod<Scalar>::weightsGaussKronrod101, QuadratureKronrod<Scalar>::weightsGauss101,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod121:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod121,
                QuadratureKronrod<Scalar>::weightsGaussKronrod121, QuadratureKronrod<Scalar>::weightsGauss121,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            case GaussKronrod201:
              return quadratureKronrodHelper(
                QuadratureKronrod<Scalar>::abscissaeGaussKronrod201,
                QuadratureKronrod<Scalar>::weightsGaussKronrod201, QuadratureKronrod<Scalar>::weightsGauss201,
                functionType, lowerLimit, upperLimit, estimatedError, absIntegral, absDiffIntegral, quadratureRule);

            default:
              return Scalar(0.);
            }
        }

        template <typename FunctionType, int numKronrodRows, int numGaussRows, int alignment>
        Scalar quadratureKronrodHelper(Array<Scalar, numKronrodRows, 1, alignment, numKronrodRows, 1> abscissaeGaussKronrod,
                                       Array<Scalar, numKronrodRows, 1, alignment, numKronrodRows, 1> weightsGaussKronrod,
                                       Array<Scalar, numGaussRows, 1, alignment, numGaussRows, 1> weightsGauss,
                                       const FunctionType& functionType,
                                       const Scalar lowerLimit,
                                       const Scalar upperLimit,
                                       Scalar& estimatedError,
                                       Scalar& absIntegral,
                                       Scalar& absDiffIntegral,
                                       const QuadratureRule quadratureRule)
        {
            using std::abs;
            using std::min;
            using std::max;
            using std::pow;

            // Half-length of the interval.
            const Scalar halfLength = (upperLimit - lowerLimit) * Scalar(.5);

            // Midpoint of the interval.
            const Scalar center = (lowerLimit + upperLimit) * Scalar(.5);
            const Scalar fCenter = functionType(center);

            Index size1 = weightsGaussKronrod.size() - 1;
            Index size2 = weightsGauss.size() - 1;

            Array<Scalar, numKronrodRows - 1, 1> f1Array;
            Array<Scalar, numKronrodRows - 1, 1> f2Array;

            // The result of the Gauss formula.
            Scalar resultGauss = Scalar(0.);

            if (quadratureRule % 2 != 0)
            {
                resultGauss = weightsGauss[size2] * fCenter;
            }

            // The result of the Kronrod formula.
            Scalar resultKronrod = weightsGaussKronrod[size1] * fCenter;            

            absIntegral = abs(resultKronrod);

            for (Index j = 1; j < weightsGaussKronrod.size() - weightsGauss.size(); ++j)
            {
                const Index jj = j * 2 - 1;
                const Scalar abscissa = halfLength * abscissaeGaussKronrod[jj];

                const Scalar f1 = functionType(center - abscissa);
                const Scalar f2 = functionType(center + abscissa);

                f1Array[jj] = f1;
                f2Array[jj] = f2;
                
                const Scalar funcSum = f1 + f2;

                resultGauss += weightsGauss[j - 1] * funcSum;
                resultKronrod += weightsGaussKronrod[jj] * funcSum;

                absIntegral += weightsGaussKronrod[jj] * (abs(f1) + abs(f2));
            }

            for (Index j = 0; j < weightsGauss.size(); ++j)
            {
                const Index jj = j * 2;
                const Scalar abscissa = halfLength * abscissaeGaussKronrod[jj];

                const Scalar f1 = functionType(center - abscissa);
                const Scalar f2 = functionType(center + abscissa);

                f1Array[jj] = f1;
                f2Array[jj] = f2;

                const Scalar funcSum = f1 + f2;

                resultKronrod += weightsGaussKronrod[jj] * funcSum;

                absIntegral += weightsGaussKronrod[jj] * (abs(f1) + abs(f2));
            }

            // Approximation to the mean value of f over the interval (lowerLimit, upperLimit),
            // i.e. I / (upperLimit - lowerLimit)
            Scalar resultMeanKronrod = resultKronrod * Scalar(.5);

            absDiffIntegral = weightsGaussKronrod[size1] * (abs(fCenter - resultMeanKronrod));

            absDiffIntegral += (((f1Array.head(size1) - resultMeanKronrod).abs()
                                + (f2Array.head(size1) - resultMeanKronrod).abs())
                                * weightsGaussKronrod.head(size1)).sum();

            Scalar result = resultKronrod * halfLength;
            absIntegral *= abs(halfLength);
            absDiffIntegral *= abs(halfLength);
            estimatedError = abs((resultKronrod - resultGauss) * halfLength);

            if (absDiffIntegral != Scalar(0.) &&
                estimatedError != Scalar(0.))
            {
                estimatedError = absDiffIntegral * min(Scalar(1.), pow((estimatedError * Scalar(200.) / absDiffIntegral), Scalar(1.5)));
            }

            if (absIntegral > (std::numeric_limits<Scalar>::min)() / (NumTraits<Scalar>::epsilon() * Scalar(50.)))
            {
                estimatedError = max(NumTraits<Scalar>::epsilon() * Scalar(50.) * absIntegral, estimatedError);
            }

            return result;
        }


        /**
         * \brief An Array of dimension m_maxSubintervals for error estimates.
         *
         * The first k elements are indices to the error estimates over the subintervals, such that
         * errorList(errorListIndices(0)), ..., errorList(errorListIndices(k - 1)) forms a decreasing
         * sequence, with k = m_numSubintervals if m_numSubintervals <= (m_maxSubintervals/2 + 2),
         * otherwise k = m_maxSubintervals + 1 - m_numSubintervals.
         */
        Array<Index, Dynamic, 1> m_errorListIndices;

        /**
         * \brief An Array of dimension m_maxSubintervals for subinterval left endpoints.
         *
         * The first m_numSubintervals elements are the lower end points of the subintervals in the
         * partition of the given integration range (lowerLimit, upperLimit).
         */
        Array<Scalar, Dynamic, 1> m_lowerList;

        /**
         * \brief An Array of dimension m_maxSubintervals for subinterval upper endpoints.
         *
         * The first m_numSubintervals elements are the upper end points of the subintervals in the
         * partition of the given integration range (lowerLimit, upperLimit).
         */
        Array<Scalar, Dynamic, 1> m_upperList;

        /**
         * \brief An Array of dimension m_maxSubintervals for integral approximations.
         *
         * The first m_numSubintervals elements are the integral approximations on the subintervals.
         */
        Array<Scalar, Dynamic, 1> m_integralList;

        /**
         * \brief An Array of dimension m_maxSubintervals for error estimates.
         *
         * The first m_numSubintervals elements of which are the moduli of the absolute error estimates
         * on the subintervals.
         */
        Array<Scalar, Dynamic, 1> m_errorList;

        /**
         * \brief Gives an upper bound on the number of subintervals. Must be at least 1.
         */
        Index m_maxSubintervals;

        /**
         * \brief The number of integrand evaluations.
         */
        Index m_numEvaluations;

        /**
         * \brief Estimate of the modulus of the absolute error, which should equal or exceed abs(I - I').
         */
        Scalar m_estimatedError;

        /**
         * \brief Error messages generated by the routine.
         *
         * errorCode = 0 Indicates normal and reliable termination of the routine. (It is assumed that
         *               the requested accuracy has been achieved.)
         * errorCode > 0 Any errorCode greater than zero indicates abnormal termination of the routine.
         *               (The estimates for integral and m_estimatedError are less reliable and the
         *               requested accuracy has not been achieved.)
         * errorCode = 1 The maximum number of subdivisions allowed has been achieved. One can allow more
         *               subdivisions by increasing the value of m_maxSubintervals. However, if this
         *               yields no improvement it is advised to analyze the integrand in order to
         *               determine the integration difficulaties. If the position of a local difficulty
         *               can be determined, (i.e. singularity or discontinuity within the interval), one
         *               will probably gain from splitting up the interval at this point and calling the
         *               integrator on the subranges. If possible, an appropriate special-purpose
         *               integrator should be used which is designed for handling the type of difficulty
         *               involved.
         * errorCode = 2 The occurrence of roundoff error is detected, preventing the requested tolerance
         *               from being achieved.
         * errorCode = 3 Extremely bad integrand behaviour occurs at points in the integration interval.
         * errorCode = 4 Roundoff error on extrapolation
         * errorCode = 5 Divergent integral (or very slowly convergent integral)
         * errorCode = 6 The input is invalid, because (desiredAbsoluteError <= 0 and
         *               desiredRealtiveError < 50 * relativeMachineAccuracy, or
         *               m_maxSubintervals < 1.
         * errorCode = 7 Applies to (D)QAWF only - limiting number of cycles has been attained
         *
         * \todo make relativeMachineAccuracy a member variable.
         */
        Index m_errorCode;

        /**
         * \brief The number of subintervals actually produced in the subdivision process.
         */
        Index m_numSubintervals;

    };
}
#endif // EIGEN_INTEGRATOR_H