lsst::afw::math::details Namespace Reference

Classes
struct	AuxFunc1
	Auxiliary struct 1. More...
struct	AuxFunc2
class	binder2_1
class	binder3_1
class	Int2DAuxType
class	Int3DAuxType

Functions
template<class T >
T	norm (const T &x)
template<class T >
T	real (const T &x)
template<class T >
T	Epsilon ()
template<class T >
T	MinRep ()
template<class T >
T	rescale_error (T err, T const &resabs, T const &resasc)
template<typename UnaryFunctionT , typename Arg >
bool	intGKPNA (UnaryFunctionT func, IntRegion< Arg > &reg, Arg const epsabs, Arg const epsrel, std::map< Arg, Arg > *fxmap=nullptr)
	Non-adaptive integration of the function f over the region 'reg'.
template<typename UnaryFunctionT , typename Arg >
void	intGKP (UnaryFunctionT func, IntRegion< Arg > &reg, Arg const epsabs, Arg const epsrel, std::map< Arg, Arg > *fxmap=nullptr)
	An adaptive integration algorithm which computes the integral of f over the region reg.
template<class UF >
AuxFunc1< UF >	Aux1 (UF uf)
	Auxiliary function 1.
template<class UF >
AuxFunc2< UF >	Aux2 (UF uf)
	Auxiliary function 2.
template<class BF , class Tp >
binder2_1< BF >	bind21 (const BF &oper, const Tp &x)
template<class TF , class Tp >
binder3_1< TF >	bind31 (const TF &oper, const Tp &x)
int	gkp_n (int level)
template<class T >
const std::vector< T > &	gkp_x (int level)
template<class T >
const std::vector< T > &	gkp_wa (int level)
template<class T >
const std::vector< T > &	gkp_wb (int level)

Function Documentation

Aux1()

template<class UF >

AuxFunc1< UF > lsst::afw::math::details::Aux1 ( UF uf )

inline

Auxiliary function 1.

Definition at line 587 of file Integrate.h.

                                {
    return AuxFunc1<UF>(uf);
}

Aux2()

template<class UF >

AuxFunc2< UF > lsst::afw::math::details::Aux2 ( UF uf )

inline

Auxiliary function 2.

Definition at line 608 of file Integrate.h.

                                {
    return AuxFunc2<UF>(uf);
}

bind21()

template<class BF , class Tp >

binder2_1< BF > lsst::afw::math::details::bind21 ( const BF & oper, const Tp & x )

inline

Definition at line 628 of file Integrate.h.

                                                         {
    using Arg = typename BF::first_argument_type;
    return binder2_1<BF>(oper, static_cast<Arg>(x));
}

bind31()

template<class TF , class Tp >

binder3_1< TF > lsst::afw::math::details::bind31 ( const TF & oper, const Tp & x )

inline

Definition at line 670 of file Integrate.h.

                                                         {
    using Arg = typename TF::firstof3_argument_type;
    return binder3_1<TF>(oper, static_cast<Arg>(x));
}

Epsilon()

template<class T >

T lsst::afw::math::details::Epsilon ( )

inline

Definition at line 244 of file Integrate.h.

                   {
    return std::numeric_limits<T>::epsilon();
}

gkp_n()

int lsst::afw::math::details::gkp_n ( int level )

inline

Definition at line 47 of file IntGKPData10.h.

                            {
    assert(level >= 0 && level < NGKPLEVELS);
    static const int ngkp[NGKPLEVELS] = {5, 5, 11, 22, 44};
    return ngkp[level];
}

gkp_wa()

template<class T >

const std::vector< T > & lsst::afw::math::details::gkp_wa ( int level )

inline

Definition at line 116 of file IntGKPData10.h.

                                             {
    // w21a, weights of the 21-point formula for abscissae x10
    static const T aw21a[5] = {0.032558162307964727478818972459390, 0.075039674810919952767043140916190,
                               0.109387158802297641899210590325805, 0.134709217311473325928054001771707,
                               0.147739104901338491374841515972068};
    static const std::vector<T> vw21a(aw21a, aw21a + 5);
 
    // w43a, weights of the 43-point formula for abscissae x10,x21
    static const T aw43a[10] = {0.016296734289666564924281974617663, 0.037522876120869501461613795898115,
                                0.054694902058255442147212685465005, 0.067355414609478086075553166302174,
                                0.073870199632393953432140695251367, 0.005768556059769796184184327908655,
                                0.027371890593248842081276069289151, 0.046560826910428830743339154433824,
                                0.061744995201442564496240336030883, 0.071387267268693397768559114425516};
    static const std::vector<T> vw43a(aw43a, aw43a + 10);
 
    // w87a, weights of the 87-point formula for abscissae x10..x43
    static const T aw87a[21] = {0.008148377384149172900002878448190, 0.018761438201562822243935059003794,
                                0.027347451050052286161582829741283, 0.033677707311637930046581056957588,
                                0.036935099820427907614589586742499, 0.002884872430211530501334156248695,
                                0.013685946022712701888950035273128, 0.023280413502888311123409291030404,
                                0.030872497611713358675466394126442, 0.035693633639418770719351355457044,
                                0.000915283345202241360843392549948, 0.005399280219300471367738743391053,
                                0.010947679601118931134327826856808, 0.016298731696787335262665703223280,
                                0.021081568889203835112433060188190, 0.025370969769253827243467999831710,
                                0.029189697756475752501446154084920, 0.032373202467202789685788194889595,
                                0.034783098950365142750781997949596, 0.036412220731351787562801163687577,
                                0.037253875503047708539592001191226};
    static const std::vector<T> vw87a(aw87a, aw87a + 21);
 
    // w175a, weights of the 175-point formula for abscissae x10..x87
    static const T aw175a[43] = {
            0.004074188692082600103872341,  0.009380719100781411819741177,  0.01367372552502614309257226,
            0.01683885365581896502443513,   0.01846754991021395380770806,   0.001442436294030254510518689,
            0.006842973011356375224430623,  0.01164020675144415562952859,   0.01543624880585667934076834,
            0.01784681681970938536027839,   0.0004576754584154192941064591, 0.002699640110136963534043642,
            0.005473839800559775445647306,  0.008149365848393670964180954,  0.01054078444460191775302900,
            0.01268548488462691364854167,   0.01459484887823787625642303,   0.01618660123360139484467359,
            0.01739154947518257137619173,   0.01820611036567589378188459,   0.01862693775152385427017140,
            0.0001363171844264931156688042, 0.0009035606060432074452289352, 0.002048434635928091782250346,
            0.003379145025868061092978424,  0.004774978836099393880602724,  0.006164723826122346701274159,
            0.007505223673194467726814647,  0.008774483993121594091281642,  0.009969018893220443741650500,
            0.01109746798050614328494604,   0.01216957356300040269315573,   0.01318725270741960360319771,
            0.01414345539438560032188640,   0.01502629056404634765715107,   0.01582337568571996469999659,
            0.01652520670998925164398162,   0.01712754985211303089259342,   0.01763120633007834051620255,
            0.01803849481144435059221422,   0.01834930224922804724856518,   0.01856027463491628805666919,
            0.01866711437596752016025132};
    static const std::vector<T> vw175a(aw175a, aw175a + 43);
 
    static const std::vector<T>* wa[NGKPLEVELS] = {nullptr, &vw21a, &vw43a, &vw87a, &vw175a};
 
    assert(level >= 1 && level < NGKPLEVELS);
    return *wa[level];
}

gkp_wb()

template<class T >

const std::vector< T > & lsst::afw::math::details::gkp_wb ( int level )

inline

Definition at line 171 of file IntGKPData10.h.

                                             {
    // w10, weights of the 10-point formula
    static const T aw10b[6] = {0.066671344308688137593568809893332, 0.149451349150580593145776339657697,
                               0.219086362515982043995534934228163, 0.269266719309996355091226921569469,
                               0.295524224714752870173892994651338, 0.0};
    static const std::vector<T> vw10b(aw10b, aw10b + 6);
 
    // w21b, weights of the 21-point formula for abscissae x21
    static const T aw21b[6] = {0.011694638867371874278064396062192, 0.054755896574351996031381300244580,
                               0.093125454583697605535065465083366, 0.123491976262065851077958109831074,
                               0.142775938577060080797094273138717, 0.149445554002916905664936468389821};
    static const std::vector<T> vw21b(aw21b, aw21b + 6);
 
    // w43b, weights of the 43-point formula for abscissae x43
    static const T aw43b[12] = {0.001844477640212414100389106552965, 0.010798689585891651740465406741293,
                                0.021895363867795428102523123075149, 0.032597463975345689443882222526137,
                                0.042163137935191811847627924327955, 0.050741939600184577780189020092084,
                                0.058379395542619248375475369330206, 0.064746404951445885544689259517511,
                                0.069566197912356484528633315038405, 0.072824441471833208150939535192842,
                                0.074507751014175118273571813842889, 0.074722147517403005594425168280423};
    static const std::vector<T> vw43b(aw43b, aw43b + 12);
 
    // w87b, weights of the 87-point formula for abscissae x87
    static const T aw87b[23] = {0.000274145563762072350016527092881, 0.001807124155057942948341311753254,
                                0.004096869282759164864458070683480, 0.006758290051847378699816577897424,
                                0.009549957672201646536053581325377, 0.012329447652244853694626639963780,
                                0.015010447346388952376697286041943, 0.017548967986243191099665352925900,
                                0.019938037786440888202278192730714, 0.022194935961012286796332102959499,
                                0.024339147126000805470360647041454, 0.026374505414839207241503786552615,
                                0.028286910788771200659968002987960, 0.030052581128092695322521110347341,
                                0.031646751371439929404586051078883, 0.033050413419978503290785944862689,
                                0.034255099704226061787082821046821, 0.035262412660156681033782717998428,
                                0.036076989622888701185500318003895, 0.036698604498456094498018047441094,
                                0.037120549269832576114119958413599, 0.037334228751935040321235449094698,
                                0.037361073762679023410321241766599};
    static const std::vector<T> vw87b(aw87b, aw87b + 23);
 
    // w175b, weights of the 175-point formula for abscissae x175
    static const T aw175b[45] = {
            0.00003901440664395060055484226, 0.0002787312407993348139961464, 0.0006672934146291138823512275,
            0.001163051749910003055308177,   0.001738533916888715565420637,  0.002369555040550195269486165,
            0.003036735036032612976379865,   0.003725394125884235586678819,  0.004424387541776355951933433,
            0.005125062882078508655875582,   0.005820601038952832388647054,  0.006505667785357412495794539,
            0.007176258976795165213169074,   0.007829642423806550517967042,  0.008464316525043036555640657,
            0.009079917879101951192512711,   0.009677029304308397008134699,  0.01025687419430508685393322,
            0.01082092935368957685380896,    0.01137052956179654620079332,   0.01190655150082017425918844,
            0.01242924137407215963971061,    0.01293819980598456566872951,   0.01343248644726851711632966,
            0.01391078107864287587682982,    0.01437154592198526807352501,   0.01481316256915705351874381,
            0.01523404482122031244940129,    0.01563274056894636856297717,   0.01600802989372176552081826,
            0.01635901148521841710605397,    0.01668515675646806540972517,   0.01698630892029276765404134,
            0.01726261506199017604299177,    0.01751439915887340284018269,   0.01774200489411242260549908,
            0.01794564958245060930824294,    0.01812532816637054898168673,   0.01828078919412312732671357,
            0.01841158014800801587767163,    0.01851713825962505903167768,   0.01859689376712028821908313,
            0.01865035760791320287947920,    0.01867717982648033250824110,   0.01868053688133951170552407};
    static const std::vector<T> vw175b(aw175b, aw175b + 45);
 
    static const std::vector<T>* wb[NGKPLEVELS] = {&vw10b, &vw21b, &vw43b, &vw87b, &vw175b};
 
    assert(level >= 0 && level < NGKPLEVELS);
    return *wb[level];
}

gkp_x()

template<class T >

const std::vector< T > & lsst::afw::math::details::gkp_x ( int level )

inline

Definition at line 54 of file IntGKPData10.h.

                                            {
    // x10, abscissae common to the 10-, 21-, 43- and 87-point rule
    static const T ax10[5] = {0.973906528517171720077964012084452, 0.865063366688984510732096688423493,
                              0.679409568299024406234327365114874, 0.433395394129247190799265943165784,
                              0.148874338981631210884826001129720};
    static const std::vector<T> vx10(ax10, ax10 + 5);
 
    // x21, abscissae common to the 21-, 43- and 87-point rule
    static const T ax21[5] = {0.995657163025808080735527280689003, 0.930157491355708226001207180059508,
                              0.780817726586416897063717578345042, 0.562757134668604683339000099272694,
                              0.294392862701460198131126603103866};
    static const std::vector<T> vx21(ax21, ax21 + 5);
 
    // x43, abscissae common to the 43- and 87-point rule
    static const T ax43[11] = {0.999333360901932081394099323919911, 0.987433402908088869795961478381209,
                               0.954807934814266299257919200290473, 0.900148695748328293625099494069092,
                               0.825198314983114150847066732588520, 0.732148388989304982612354848755461,
                               0.622847970537725238641159120344323, 0.499479574071056499952214885499755,
                               0.364901661346580768043989548502644, 0.222254919776601296498260928066212,
                               0.074650617461383322043914435796506};
    static const std::vector<T> vx43(ax43, ax43 + 11);
 
    // x87, abscissae of the 87-point rule
    static const T ax87[22] = {0.999902977262729234490529830591582, 0.997989895986678745427496322365960,
                               0.992175497860687222808523352251425, 0.981358163572712773571916941623894,
                               0.965057623858384619128284110607926, 0.943167613133670596816416634507426,
                               0.915806414685507209591826430720050, 0.883221657771316501372117548744163,
                               0.845710748462415666605902011504855, 0.803557658035230982788739474980964,
                               0.757005730685495558328942793432020, 0.706273209787321819824094274740840,
                               0.651589466501177922534422205016736, 0.593223374057961088875273770349144,
                               0.531493605970831932285268948562671, 0.466763623042022844871966781659270,
                               0.399424847859218804732101665817923, 0.329874877106188288265053371824597,
                               0.258503559202161551802280975429025, 0.185695396568346652015917141167606,
                               0.111842213179907468172398359241362, 0.037352123394619870814998165437704};
    static const std::vector<T> vx87(ax87, ax87 + 22);
 
    // x175, new abscissae of the 175-point rule
    static const T ax175[44] = {
            0.9999863601049729677997719,  0.9996988005421428623760900, 0.9987732473138838303043008,
            0.9969583851613681945576954,  0.9940679446543307664316419, 0.9899673391741082769960501,
            0.9845657252243253759726856,  0.9778061574856414566055911, 0.9696573134362640567828241,
            0.9601075259508299668299148,  0.9491605191933640559839528, 0.9368321320434676982049391,
            0.9231475266471484946315312,  0.9081385954543932335481458, 0.8918414554201179766275301,
            0.8742940658668577971868106,  0.8555341249920432726468181, 0.8355974653331171075829366,
            0.8145171513255796294567278,  0.7923233739697216224661771, 0.7690440708082055924983192,
            0.7447060413364141221086362,  0.7193362474393392286563226, 0.6929630148748741413274435,
            0.6656169582351898081906710,  0.6373315753624730450361940, 0.6081435407529742253347887,
            0.5780927509986774074337714,  0.5472221576559162873479056, 0.5155774001694817879753989,
            0.4832062516695473663839323,  0.4501579202236101928554315, 0.4164822936334372150371142,
            0.3822292525735955063074986,  0.3474481818229070516859048, 0.3121877718977455143765043,
            0.2764961311348905899785770,  0.2404211453027010909070665, 0.2040109589145005126324680,
            0.1673144327616352753945014,  0.1303814601095591675431696, 0.09326308333660176874846400,
            0.05601141598806823374241435, 0.01867941779948308845140053};
    static const std::vector<T> vx175(ax175, ax175 + 44);
 
    static const std::vector<T>* x[NGKPLEVELS] = {&vx10, &vx21, &vx43, &vx87, &vx175};
 
    assert(level >= 0 && level < NGKPLEVELS);
    return *x[level];
}

intGKP()

template<typename UnaryFunctionT , typename Arg >

void lsst::afw::math::details::intGKP ( UnaryFunctionT func, IntRegion< Arg > & reg, Arg const epsabs, Arg const epsrel, std::map< Arg, Arg > * fxmap = nullptr )

inline

An adaptive integration algorithm which computes the integral of f over the region reg.

Note

First the non-adaptive GKP algorithm is tried. If that is not accurate enough (according to the absolute and relative accuracies, epsabs and epsrel), the region is split in half, and each new region is integrated. The routine continues by successively splitting the subregion which gave the largest absolute error until the integral converges.

The area and estimated error are returned as reg.Area() and reg.Err() If desired, *retx and *retf return std::vectors of x,f(x) respectively They only include the evaluations in the non-adaptive pass, so they do not give an accurate estimate of the number of function evaluations.

Definition at line 426 of file Integrate.h.

                                                      {
    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;
    Arg roundoffType2 = 0;
    size_t iteration = 1;
 
    std::priority_queue<IntRegion<Arg>, std::vector<IntRegion<Arg> > > allregions;
    allregions.push(reg);
    Arg finalarea = reg.Area();
    Arg finalerr = reg.Err();
    Arg 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<Arg> 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<Arg> > 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
        Arg factor = 3 * children.size() * finalerr / tolerance;
        Arg newepsabs = fabs(parent.Err() / factor);
        Arg newepsrel = newepsabs / fabs(parent.Area());
        integ_dbg2 << "New epsabs, rel = " << newepsabs << ", " << newepsrel;
        integ_dbg2 << "  (" << children.size() << " children)\n";
 
        Arg newarea = Arg(0.0);
        Arg newerror = 0.0;
        for (size_t i = 0; i < children.size(); i++) {
            IntRegion<Arg> &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();
 
        Arg 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, Arg(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<Arg> 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());
    }
}

intGKPNA()

template<typename UnaryFunctionT , typename Arg >

bool lsst::afw::math::details::intGKPNA ( UnaryFunctionT func, IntRegion< Arg > & reg, Arg const epsabs, Arg const epsrel, std::map< Arg, Arg > * fxmap = nullptr )

inline

Non-adaptive integration of the function f over the region 'reg'.

Note

The algorithm computes first a Gaussian quadrature value then successive Kronrod/Patterson extensions to this result. The functions terminates when the difference between successive approximations (rescaled according to rescale_error) is less than either epsabs or epsrel * I, where I is the latest estimate of the integral. The order of the Gauss/Kronron/Patterson scheme is determined by which file is included above. Currently schemes starting with order 1 and order 10 are calculated. There seems to be little practical difference in the integration times using the two schemes, so I haven't bothered to calculate any more.

Definition at line 298 of file Integrate.h.

                                                        {
    Arg const a = reg.Left();
    Arg const b = reg.Right();
 
    Arg const halfLength = 0.5 * (b - a);
    Arg const absHalfLength = fabs(halfLength);
    Arg const center = 0.5 * (b + a);
    Arg const fCenter = func(center);
#ifdef COUNTFEVAL
    nfeval++;
#endif
 
    assert(gkp_wb<Arg>(0).size() == gkp_x<Arg>(0).size() + 1);
    Arg area1 = gkp_wb<Arg>(0).back() * fCenter;
    std::vector<Arg> fv1, fv2;
    fv1.reserve(2 * gkp_x<Arg>(0).size() + 1);
    fv2.reserve(2 * gkp_x<Arg>(0).size() + 1);
    for (size_t k = 0; k < gkp_x<Arg>(0).size(); k++) {
        Arg const abscissa = halfLength * gkp_x<Arg>(0)[k];
        Arg const fval1 = func(center - abscissa);
        Arg const fval2 = func(center + abscissa);
        area1 += gkp_wb<Arg>(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<Arg>(0).size() * 2;
#endif
 
    integ_dbg2 << "level 0 rule: area = " << area1 << std::endl;
 
    Arg err = 0;
    bool calcabsasc = true;
    Arg resabs = 0.0, resasc = 0.0;
    for (int level = 1; level < NGKPLEVELS; level++) {
        assert(gkp_wa<Arg>(level).size() == fv1.size());
        assert(gkp_wa<Arg>(level).size() == fv2.size());
        assert(gkp_wb<Arg>(level).size() == gkp_x<Arg>(level).size() + 1);
        Arg area2 = gkp_wb<Arg>(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<Arg>(level)[k] * (fv1[k] + fv2[k]);
            if (calcabsasc) {
                resabs += gkp_wa<Arg>(level)[k] * (fabs(fv1[k]) + fabs(fv2[k]));
            }
        }
        for (size_t k = 0; k < gkp_x<Arg>(level).size(); k++) {
            Arg const abscissa = halfLength * gkp_x<Arg>(level)[k];
            Arg const fval1 = func(center - abscissa);
            Arg const fval2 = func(center + abscissa);
            Arg const fval = fval1 + fval2;
            area2 += gkp_wb<Arg>(level)[k] * fval;
            if (calcabsasc) {
                resabs += gkp_wb<Arg>(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<Arg>(level).size() * 2;
#endif
        if (calcabsasc) {
            Arg const mean = area1 * Arg(0.5);
            // resasc = approximation to the integral of abs(f-mean)
            resasc = gkp_wb<Arg>(level).back() * fabs(fCenter - mean);
            for (size_t k = 0; k < gkp_wa<Arg>(level).size(); k++) {
                resasc += gkp_wa<Arg>(level)[k] * (fabs(fv1[k] - mean) + fabs(fv2[k] - mean));
            }
            for (size_t k = 0; k < gkp_x<Arg>(level).size(); k++) {
                resasc += gkp_wb<Arg>(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;
}

MinRep()

template<class T >

T lsst::afw::math::details::MinRep ( )

inline

Definition at line 248 of file Integrate.h.

                  {
    return std::numeric_limits<T>::min();
}

norm()

template<class T >

T lsst::afw::math::details::norm ( const T & x )

inline

Definition at line 160 of file Integrate.h.

                          {
    return x * x;
}

real()

template<class T >

T lsst::afw::math::details::real ( const T & x )

inline

Definition at line 165 of file Integrate.h.

                          {
    return x;
}

rescale_error()

template<class T >

T lsst::afw::math::details::rescale_error ( T err, T const & resabs, T const & resasc )

inline

Definition at line 264 of file Integrate.h.

                                                                {
    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;
}