Classes
struct	AuxFunc1
	Auxiliary struct 1. More...

struct	AuxFunc2

struct	ConstantReg1
	Helpers for constant regions for int2d, int3d: More...

struct	ConstantReg2

class	binder2_1

class	Int2DAuxType

class	binder3_1

class	Int3DAuxType

class	FunctionWrapper
	Wrap an integrand in a call to a 1D integrator: romberg() More...

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<class UF >
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)
	Non-adaptive integration of the function f over the region 'reg'. More...

template<class UF >
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)
	An adaptive integration algorithm which computes the integral of f over the region reg. More...

template<class UF >
AuxFunc1< UF >	Aux1 (UF uf)
	Auxiliary function 1. More...

template<class UF >
AuxFunc2< UF >	Aux2 (UF uf)
	Auxiliary function 2. More...

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)

template<typename A , typename B >
bool	isSameObject (A const &, B const &)

template<typename A >
bool	isSameObject (A const &a, A const &b)

Function Documentation

template<class UF >

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

inline

Auxiliary function 1.

Definition at line 620 of file Integrate.h.

620 { return AuxFunc1<UF>(uf); }

template<class UF >

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

inline

Auxiliary function 2.

Definition at line 641 of file Integrate.h.

641 { return AuxFunc2<UF>(uf); }

template<class BF , class Tp >

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

inline

Definition at line 685 of file Integrate.h.

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

template<class TF , class Tp >

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

inline

Definition at line 737 of file Integrate.h.

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

template<class T >

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

inline

Definition at line 266 of file Integrate.h.

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

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

inline

Definition at line 45 of file IntGKPData10.h.

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

template<class T >

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

inline

Definition at line 173 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] = {0, &vw21a, &vw43a, &vw87a, &vw175a};
 
         assert(level >= 1 && level < NGKPLEVELS);
         return *wa[level];
     }

template<class T >

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

inline

Definition at line 281 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];
   }

template<class T >

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

inline

Definition at line 52 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];
     }

template<class UF >

void lsst::afw::math::details::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`
	)

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 450 of file Integrate.h.

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

bool lsst::afw::math::details::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`
	)

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 316 of file Integrate.h.

                                                                                         {
     
     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<typename A , typename B >

bool lsst::afw::math::details::isSameObject	(	A const &	,
		B const &
	)

Definition at line 635 of file warpExposure.h.

635 { return false; }

template<typename A >

bool lsst::afw::math::details::isSameObject	(	A const &	a,
		A const &	b
	)

Definition at line 638 of file warpExposure.h.

638 { return &a == &b; }

b

afw::table::Key< double > b

Definition: DoubleGaussianPsf.cc:49

template<class T >

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

inline

Definition at line 268 of file Integrate.h.

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

template<class T >

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

inline

Definition at line 191 of file Integrate.h.

191 { return x*x; }

x

double x

Definition: ChebyshevBoundedField.cc:305

template<class T >

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

inline

Definition at line 193 of file Integrate.h.

193 { return x; }

x

double x

Definition: ChebyshevBoundedField.cc:305

template<class T >

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

inline

Definition at line 280 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;
 }

Classes

Functions

Function Documentation