LSSTApplications  18.0.0+106,18.0.0+50,19.0.0,19.0.0+1,19.0.0+10,19.0.0+11,19.0.0+13,19.0.0+17,19.0.0+2,19.0.0-1-g20d9b18+6,19.0.0-1-g425ff20,19.0.0-1-g5549ca4,19.0.0-1-g580fafe+6,19.0.0-1-g6fe20d0+1,19.0.0-1-g7011481+9,19.0.0-1-g8c57eb9+6,19.0.0-1-gb5175dc+11,19.0.0-1-gdc0e4a7+9,19.0.0-1-ge272bc4+6,19.0.0-1-ge3aa853,19.0.0-10-g448f008b,19.0.0-12-g6990b2c,19.0.0-2-g0d9f9cd+11,19.0.0-2-g3d9e4fb2+11,19.0.0-2-g5037de4,19.0.0-2-gb96a1c4+3,19.0.0-2-gd955cfd+15,19.0.0-3-g2d13df8,19.0.0-3-g6f3c7dc,19.0.0-4-g725f80e+11,19.0.0-4-ga671dab3b+1,19.0.0-4-gad373c5+3,19.0.0-5-ga2acb9c+2,19.0.0-5-gfe96e6c+2,w.2020.01
LSSTDataManagementBasePackage
Classes | Functions
lsst::afw::math::details Namespace Reference

Classes

struct  AuxFunc1
 Auxiliary struct 1. More...
 
struct  AuxFunc2
 
class  binder2_1
 
class  binder3_1
 
struct  ConstantReg1
 Helpers for constant regions for int2d, int3d: More...
 
struct  ConstantReg2
 
class  FunctionWrapper
 Wrap an integrand in a call to a 1D integrator: romberg() More...
 
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<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)
 

Function Documentation

◆ Aux1()

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

Auxiliary function 1.

Definition at line 626 of file Integrate.h.

626  {
627  return AuxFunc1<UF>(uf);
628 }

◆ Aux2()

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

Auxiliary function 2.

Definition at line 648 of file Integrate.h.

648  {
649  return AuxFunc2<UF>(uf);
650 }

◆ bind21()

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

Definition at line 687 of file Integrate.h.

687  {
688  typedef typename BF::first_argument_type Arg;
689  return binder2_1<BF>(oper, static_cast<Arg>(x));
690 }
x
double x
Definition: ChebyshevBoundedField.cc:277

◆ bind31()

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

Definition at line 730 of file Integrate.h.

730  {
731  typedef typename TF::firstof3_argument_type Arg;
732  return binder3_1<TF>(oper, static_cast<Arg>(x));
733 }
x
double x
Definition: ChebyshevBoundedField.cc:277

◆ Epsilon()

template<class T >
T lsst::afw::math::details::Epsilon ( )
inline

Definition at line 278 of file Integrate.h.

278  {
279  return std::numeric_limits<T>::epsilon();
280 }
std::numeric_limits::epsilon
T epsilon(T... args)

◆ gkp_n()

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

Definition at line 47 of file IntGKPData10.h.

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

◆ 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.

116  {
117  // w21a, weights of the 21-point formula for abscissae x10
118  static const T aw21a[5] = {0.032558162307964727478818972459390, 0.075039674810919952767043140916190,
119  0.109387158802297641899210590325805, 0.134709217311473325928054001771707,
120  0.147739104901338491374841515972068};
121  static const std::vector<T> vw21a(aw21a, aw21a + 5);
122 
123  // w43a, weights of the 43-point formula for abscissae x10,x21
124  static const T aw43a[10] = {0.016296734289666564924281974617663, 0.037522876120869501461613795898115,
125  0.054694902058255442147212685465005, 0.067355414609478086075553166302174,
126  0.073870199632393953432140695251367, 0.005768556059769796184184327908655,
127  0.027371890593248842081276069289151, 0.046560826910428830743339154433824,
128  0.061744995201442564496240336030883, 0.071387267268693397768559114425516};
129  static const std::vector<T> vw43a(aw43a, aw43a + 10);
130 
131  // w87a, weights of the 87-point formula for abscissae x10..x43
132  static const T aw87a[21] = {0.008148377384149172900002878448190, 0.018761438201562822243935059003794,
133  0.027347451050052286161582829741283, 0.033677707311637930046581056957588,
134  0.036935099820427907614589586742499, 0.002884872430211530501334156248695,
135  0.013685946022712701888950035273128, 0.023280413502888311123409291030404,
136  0.030872497611713358675466394126442, 0.035693633639418770719351355457044,
137  0.000915283345202241360843392549948, 0.005399280219300471367738743391053,
138  0.010947679601118931134327826856808, 0.016298731696787335262665703223280,
139  0.021081568889203835112433060188190, 0.025370969769253827243467999831710,
140  0.029189697756475752501446154084920, 0.032373202467202789685788194889595,
141  0.034783098950365142750781997949596, 0.036412220731351787562801163687577,
142  0.037253875503047708539592001191226};
143  static const std::vector<T> vw87a(aw87a, aw87a + 21);
144 
145  // w175a, weights of the 175-point formula for abscissae x10..x87
146  static const T aw175a[43] = {
147  0.004074188692082600103872341, 0.009380719100781411819741177, 0.01367372552502614309257226,
148  0.01683885365581896502443513, 0.01846754991021395380770806, 0.001442436294030254510518689,
149  0.006842973011356375224430623, 0.01164020675144415562952859, 0.01543624880585667934076834,
150  0.01784681681970938536027839, 0.0004576754584154192941064591, 0.002699640110136963534043642,
151  0.005473839800559775445647306, 0.008149365848393670964180954, 0.01054078444460191775302900,
152  0.01268548488462691364854167, 0.01459484887823787625642303, 0.01618660123360139484467359,
153  0.01739154947518257137619173, 0.01820611036567589378188459, 0.01862693775152385427017140,
154  0.0001363171844264931156688042, 0.0009035606060432074452289352, 0.002048434635928091782250346,
155  0.003379145025868061092978424, 0.004774978836099393880602724, 0.006164723826122346701274159,
156  0.007505223673194467726814647, 0.008774483993121594091281642, 0.009969018893220443741650500,
157  0.01109746798050614328494604, 0.01216957356300040269315573, 0.01318725270741960360319771,
158  0.01414345539438560032188640, 0.01502629056404634765715107, 0.01582337568571996469999659,
159  0.01652520670998925164398162, 0.01712754985211303089259342, 0.01763120633007834051620255,
160  0.01803849481144435059221422, 0.01834930224922804724856518, 0.01856027463491628805666919,
161  0.01866711437596752016025132};
162  static const std::vector<T> vw175a(aw175a, aw175a + 43);
163 
164  static const std::vector<T>* wa[NGKPLEVELS] = {0, &vw21a, &vw43a, &vw87a, &vw175a};
165 
166  assert(level >= 1 && level < NGKPLEVELS);
167  return *wa[level];
168 }
std::vector
STL class.

◆ 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.

171  {
172  // w10, weights of the 10-point formula
173  static const T aw10b[6] = {0.066671344308688137593568809893332, 0.149451349150580593145776339657697,
174  0.219086362515982043995534934228163, 0.269266719309996355091226921569469,
175  0.295524224714752870173892994651338, 0.0};
176  static const std::vector<T> vw10b(aw10b, aw10b + 6);
177 
178  // w21b, weights of the 21-point formula for abscissae x21
179  static const T aw21b[6] = {0.011694638867371874278064396062192, 0.054755896574351996031381300244580,
180  0.093125454583697605535065465083366, 0.123491976262065851077958109831074,
181  0.142775938577060080797094273138717, 0.149445554002916905664936468389821};
182  static const std::vector<T> vw21b(aw21b, aw21b + 6);
183 
184  // w43b, weights of the 43-point formula for abscissae x43
185  static const T aw43b[12] = {0.001844477640212414100389106552965, 0.010798689585891651740465406741293,
186  0.021895363867795428102523123075149, 0.032597463975345689443882222526137,
187  0.042163137935191811847627924327955, 0.050741939600184577780189020092084,
188  0.058379395542619248375475369330206, 0.064746404951445885544689259517511,
189  0.069566197912356484528633315038405, 0.072824441471833208150939535192842,
190  0.074507751014175118273571813842889, 0.074722147517403005594425168280423};
191  static const std::vector<T> vw43b(aw43b, aw43b + 12);
192 
193  // w87b, weights of the 87-point formula for abscissae x87
194  static const T aw87b[23] = {0.000274145563762072350016527092881, 0.001807124155057942948341311753254,
195  0.004096869282759164864458070683480, 0.006758290051847378699816577897424,
196  0.009549957672201646536053581325377, 0.012329447652244853694626639963780,
197  0.015010447346388952376697286041943, 0.017548967986243191099665352925900,
198  0.019938037786440888202278192730714, 0.022194935961012286796332102959499,
199  0.024339147126000805470360647041454, 0.026374505414839207241503786552615,
200  0.028286910788771200659968002987960, 0.030052581128092695322521110347341,
201  0.031646751371439929404586051078883, 0.033050413419978503290785944862689,
202  0.034255099704226061787082821046821, 0.035262412660156681033782717998428,
203  0.036076989622888701185500318003895, 0.036698604498456094498018047441094,
204  0.037120549269832576114119958413599, 0.037334228751935040321235449094698,
205  0.037361073762679023410321241766599};
206  static const std::vector<T> vw87b(aw87b, aw87b + 23);
207 
208  // w175b, weights of the 175-point formula for abscissae x175
209  static const T aw175b[45] = {
210  0.00003901440664395060055484226, 0.0002787312407993348139961464, 0.0006672934146291138823512275,
211  0.001163051749910003055308177, 0.001738533916888715565420637, 0.002369555040550195269486165,
212  0.003036735036032612976379865, 0.003725394125884235586678819, 0.004424387541776355951933433,
213  0.005125062882078508655875582, 0.005820601038952832388647054, 0.006505667785357412495794539,
214  0.007176258976795165213169074, 0.007829642423806550517967042, 0.008464316525043036555640657,
215  0.009079917879101951192512711, 0.009677029304308397008134699, 0.01025687419430508685393322,
216  0.01082092935368957685380896, 0.01137052956179654620079332, 0.01190655150082017425918844,
217  0.01242924137407215963971061, 0.01293819980598456566872951, 0.01343248644726851711632966,
218  0.01391078107864287587682982, 0.01437154592198526807352501, 0.01481316256915705351874381,
219  0.01523404482122031244940129, 0.01563274056894636856297717, 0.01600802989372176552081826,
220  0.01635901148521841710605397, 0.01668515675646806540972517, 0.01698630892029276765404134,
221  0.01726261506199017604299177, 0.01751439915887340284018269, 0.01774200489411242260549908,
222  0.01794564958245060930824294, 0.01812532816637054898168673, 0.01828078919412312732671357,
223  0.01841158014800801587767163, 0.01851713825962505903167768, 0.01859689376712028821908313,
224  0.01865035760791320287947920, 0.01867717982648033250824110, 0.01868053688133951170552407};
225  static const std::vector<T> vw175b(aw175b, aw175b + 45);
226 
227  static const std::vector<T>* wb[NGKPLEVELS] = {&vw10b, &vw21b, &vw43b, &vw87b, &vw175b};
228 
229  assert(level >= 0 && level < NGKPLEVELS);
230  return *wb[level];
231 }
std::vector
STL class.

◆ 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.

54  {
55  // x10, abscissae common to the 10-, 21-, 43- and 87-point rule
56  static const T ax10[5] = {0.973906528517171720077964012084452, 0.865063366688984510732096688423493,
57  0.679409568299024406234327365114874, 0.433395394129247190799265943165784,
58  0.148874338981631210884826001129720};
59  static const std::vector<T> vx10(ax10, ax10 + 5);
60 
61  // x21, abscissae common to the 21-, 43- and 87-point rule
62  static const T ax21[5] = {0.995657163025808080735527280689003, 0.930157491355708226001207180059508,
63  0.780817726586416897063717578345042, 0.562757134668604683339000099272694,
64  0.294392862701460198131126603103866};
65  static const std::vector<T> vx21(ax21, ax21 + 5);
66 
67  // x43, abscissae common to the 43- and 87-point rule
68  static const T ax43[11] = {0.999333360901932081394099323919911, 0.987433402908088869795961478381209,
69  0.954807934814266299257919200290473, 0.900148695748328293625099494069092,
70  0.825198314983114150847066732588520, 0.732148388989304982612354848755461,
71  0.622847970537725238641159120344323, 0.499479574071056499952214885499755,
72  0.364901661346580768043989548502644, 0.222254919776601296498260928066212,
73  0.074650617461383322043914435796506};
74  static const std::vector<T> vx43(ax43, ax43 + 11);
75 
76  // x87, abscissae of the 87-point rule
77  static const T ax87[22] = {0.999902977262729234490529830591582, 0.997989895986678745427496322365960,
78  0.992175497860687222808523352251425, 0.981358163572712773571916941623894,
79  0.965057623858384619128284110607926, 0.943167613133670596816416634507426,
80  0.915806414685507209591826430720050, 0.883221657771316501372117548744163,
81  0.845710748462415666605902011504855, 0.803557658035230982788739474980964,
82  0.757005730685495558328942793432020, 0.706273209787321819824094274740840,
83  0.651589466501177922534422205016736, 0.593223374057961088875273770349144,
84  0.531493605970831932285268948562671, 0.466763623042022844871966781659270,
85  0.399424847859218804732101665817923, 0.329874877106188288265053371824597,
86  0.258503559202161551802280975429025, 0.185695396568346652015917141167606,
87  0.111842213179907468172398359241362, 0.037352123394619870814998165437704};
88  static const std::vector<T> vx87(ax87, ax87 + 22);
89 
90  // x175, new abscissae of the 175-point rule
91  static const T ax175[44] = {
92  0.9999863601049729677997719, 0.9996988005421428623760900, 0.9987732473138838303043008,
93  0.9969583851613681945576954, 0.9940679446543307664316419, 0.9899673391741082769960501,
94  0.9845657252243253759726856, 0.9778061574856414566055911, 0.9696573134362640567828241,
95  0.9601075259508299668299148, 0.9491605191933640559839528, 0.9368321320434676982049391,
96  0.9231475266471484946315312, 0.9081385954543932335481458, 0.8918414554201179766275301,
97  0.8742940658668577971868106, 0.8555341249920432726468181, 0.8355974653331171075829366,
98  0.8145171513255796294567278, 0.7923233739697216224661771, 0.7690440708082055924983192,
99  0.7447060413364141221086362, 0.7193362474393392286563226, 0.6929630148748741413274435,
100  0.6656169582351898081906710, 0.6373315753624730450361940, 0.6081435407529742253347887,
101  0.5780927509986774074337714, 0.5472221576559162873479056, 0.5155774001694817879753989,
102  0.4832062516695473663839323, 0.4501579202236101928554315, 0.4164822936334372150371142,
103  0.3822292525735955063074986, 0.3474481818229070516859048, 0.3121877718977455143765043,
104  0.2764961311348905899785770, 0.2404211453027010909070665, 0.2040109589145005126324680,
105  0.1673144327616352753945014, 0.1303814601095591675431696, 0.09326308333660176874846400,
106  0.05601141598806823374241435, 0.01867941779948308845140053};
107  static const std::vector<T> vx175(ax175, ax175 + 44);
108 
109  static const std::vector<T>* x[NGKPLEVELS] = {&vx10, &vx21, &vx43, &vx87, &vx175};
110 
111  assert(level >= 0 && level < NGKPLEVELS);
112  return *x[level];
113 }
x
double x
Definition: ChebyshevBoundedField.cc:277
std::vector
STL class.

◆ intGKP()

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

465  {
466  typedef typename UF::result_type UfResult;
467  integ_dbg2 << "Start intGKP\n";
468 
469  assert(epsabs >= 0.0);
470  assert(epsrel > 0.0);
471 
472  // perform the first integration
473  bool done = intGKPNA(func, reg, epsabs, epsrel, fxmap);
474  if (done) return;
475 
476  integ_dbg2 << "In adaptive GKP, failed first pass... subdividing\n";
477  integ_dbg2 << "Intial range = " << reg.Left() << ".." << reg.Right() << std::endl;
478 
479  int roundoffType1 = 0, errorType = 0;
480  UfResult roundoffType2 = 0;
481  size_t iteration = 1;
482 
483  std::priority_queue<IntRegion<UfResult>, std::vector<IntRegion<UfResult> > > allregions;
484  allregions.push(reg);
485  UfResult finalarea = reg.Area();
486  UfResult finalerr = reg.Err();
487  UfResult tolerance = std::max(epsabs, epsrel * fabs(finalarea));
488  assert(finalerr > tolerance);
489 
490  while (!errorType && finalerr > tolerance) {
491  // Bisect the subinterval with the largest error estimate
492  integ_dbg2 << "Current answer = " << finalarea << " +- " << finalerr;
493  integ_dbg2 << " (tol = " << tolerance << ")\n";
494  IntRegion<UfResult> parent = allregions.top();
495  allregions.pop();
496  integ_dbg2 << "Subdividing largest error region ";
497  integ_dbg2 << parent.Left() << ".." << parent.Right() << std::endl;
498  integ_dbg2 << "parent area = " << parent.Area();
499  integ_dbg2 << " +- " << parent.Err() << std::endl;
500  std::vector<IntRegion<UfResult> > children;
501  parent.SubDivide(&children);
502  // For "GKP", there are only two, but for GKPOSC, there is one
503  // for each oscillation in region
504 
505  // Try to do at least 3x better with the children
506  UfResult factor = 3 * children.size() * finalerr / tolerance;
507  UfResult newepsabs = fabs(parent.Err() / factor);
508  UfResult newepsrel = newepsabs / fabs(parent.Area());
509  integ_dbg2 << "New epsabs, rel = " << newepsabs << ", " << newepsrel;
510  integ_dbg2 << " (" << children.size() << " children)\n";
511 
512  UfResult newarea = UfResult(0.0);
513  UfResult newerror = 0.0;
514  for (size_t i = 0; i < children.size(); i++) {
515  IntRegion<UfResult> &child = children[i];
516  integ_dbg2 << "Integrating child " << child.Left();
517  integ_dbg2 << ".." << child.Right() << std::endl;
518  bool hasConverged;
519  hasConverged = intGKPNA(func, child, newepsabs, newepsrel);
520  integ_dbg2 << "child (" << i + 1 << '/' << children.size() << ") ";
521  if (hasConverged) {
522  integ_dbg2 << " converged.";
523  } else {
524  integ_dbg2 << " failed.";
525  }
526  integ_dbg2 << " Area = " << child.Area() << " +- " << child.Err() << std::endl;
527 
528  newarea += child.Area();
529  newerror += child.Err();
530  }
531  integ_dbg2 << "Compare: newerr = " << newerror;
532  integ_dbg2 << " to parent err = " << parent.Err() << std::endl;
533 
534  finalerr += (newerror - parent.Err());
535  finalarea += newarea - parent.Area();
536 
537  UfResult delta = parent.Area() - newarea;
538  if (newerror <= parent.Err() && fabs(delta) <= parent.Err() && newerror >= 0.99 * parent.Err()) {
539  integ_dbg2 << "roundoff type 1: delta/newarea = ";
540  integ_dbg2 << fabs(delta) / fabs(newarea);
541  integ_dbg2 << ", newerror/error = " << newerror / parent.Err() << std::endl;
542  roundoffType1++;
543  }
544  if (iteration >= 10 && newerror > parent.Err() && fabs(delta) <= newerror - parent.Err()) {
545  integ_dbg2 << "roundoff type 2: newerror/error = ";
546  integ_dbg2 << newerror / parent.Err() << std::endl;
547  roundoffType2 += std::min(newerror / parent.Err() - 1.0, UfResult(1.0));
548  }
549 
550  tolerance = std::max(epsabs, epsrel * fabs(finalarea));
551  if (finalerr > tolerance) {
552  if (roundoffType1 >= 200) {
553  errorType = 1; // round off error
554  integ_dbg2 << "GKP: Round off error 1\n";
555  }
556  if (roundoffType2 >= 200.0) {
557  errorType = 2; // round off error
558  integ_dbg2 << "GKP: Round off error 2\n";
559  }
560  if (fabs((parent.Right() - parent.Left()) / (reg.Right() - reg.Left())) < Epsilon<double>()) {
561  errorType = 3; // found singularity
562  integ_dbg2 << "GKP: Probable singularity\n";
563  }
564  }
565  for (size_t i = 0; i < children.size(); i++) {
566  allregions.push(children[i]);
567  }
568  iteration++;
569  }
570 
571  // Recalculate finalarea in case there are any slight rounding errors
572  finalarea = 0.0;
573  finalerr = 0.0;
574  while (!allregions.empty()) {
575  IntRegion<UfResult> const &r = allregions.top();
576  finalarea += r.Area();
577  finalerr += r.Err();
578  allregions.pop();
579  }
580  reg.SetArea(finalarea, finalerr);
581 
582  if (errorType == 1) {
583  std::ostringstream s;
584  s << "Type 1 roundoff's = " << roundoffType1;
585  s << ", Type 2 = " << roundoffType2 << std::endl;
586  s << "Roundoff error 1 prevents tolerance from being achieved ";
587  s << "in intGKP\n";
588  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, s.str());
589  } else if (errorType == 2) {
590  std::ostringstream s;
591  s << "Type 1 roundoff's = " << roundoffType1;
592  s << ", Type 2 = " << roundoffType2 << std::endl;
593  s << "Roundoff error 2 prevents tolerance from being achieved ";
594  s << "in intGKP\n";
595  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, s.str());
596  } else if (errorType == 3) {
597  std::ostringstream s;
598  s << "Bad integrand behavior found in the integration interval ";
599  s << "in intGKP\n";
600  throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, s.str());
601  }
602 }
std::priority_queue::empty
T empty(T... args)
std::endl
T endl(T... args)
std::priority_queue::pop
T pop(T... args)
std::priority_queue::top
T top(T... args)
std::priority_queue::push
T push(T... args)
std::min
T min(T... args)
std::ostringstream::str
T str(T... args)
std::fabs
T fabs(T... args)
std::max
T max(T... args)
std::vector::size
T size(T... args)
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::vector
STL class.
integ_dbg2
#define integ_dbg2
Definition: Integrate.h:183
lsst::afw::math::details::intGKPNA
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 &#39;reg&#39;.
Definition: Integrate.h:333
lsst::pex::exceptions::RuntimeError
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104

◆ intGKPNA()

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

335  {
336  typedef typename UF::result_type UfResult;
337  UfResult const a = reg.Left();
338  UfResult const b = reg.Right();
339 
340  UfResult const halfLength = 0.5 * (b - a);
341  UfResult const absHalfLength = fabs(halfLength);
342  UfResult const center = 0.5 * (b + a);
343  UfResult const fCenter = func(center);
344 #ifdef COUNTFEVAL
345  nfeval++;
346 #endif
347 
348  assert(gkp_wb<UfResult>(0).size() == gkp_x<UfResult>(0).size() + 1);
349  UfResult area1 = gkp_wb<UfResult>(0).back() * fCenter;
350  std::vector<UfResult> fv1, fv2;
351  fv1.reserve(2 * gkp_x<UfResult>(0).size() + 1);
352  fv2.reserve(2 * gkp_x<UfResult>(0).size() + 1);
353  for (size_t k = 0; k < gkp_x<UfResult>(0).size(); k++) {
354  UfResult const abscissa = halfLength * gkp_x<UfResult>(0)[k];
355  UfResult const fval1 = func(center - abscissa);
356  UfResult const fval2 = func(center + abscissa);
357  area1 += gkp_wb<UfResult>(0)[k] * (fval1 + fval2);
358  fv1.push_back(fval1);
359  fv2.push_back(fval2);
360  if (fxmap) {
361  (*fxmap)[center - abscissa] = fval1;
362  (*fxmap)[center + abscissa] = fval2;
363  }
364  }
365 #ifdef COUNTFEVAL
366  nfeval += gkp_x<UfResult>(0).size() * 2;
367 #endif
368 
369  integ_dbg2 << "level 0 rule: area = " << area1 << std::endl;
370 
371  UfResult err = 0;
372  bool calcabsasc = true;
373  UfResult resabs = 0.0, resasc = 0.0;
374  for (int level = 1; level < NGKPLEVELS; level++) {
375  assert(gkp_wa<UfResult>(level).size() == fv1.size());
376  assert(gkp_wa<UfResult>(level).size() == fv2.size());
377  assert(gkp_wb<UfResult>(level).size() == gkp_x<UfResult>(level).size() + 1);
378  UfResult area2 = gkp_wb<UfResult>(level).back() * fCenter;
379  // resabs = approximation to integral of abs(f)
380  if (calcabsasc) {
381  resabs = fabs(area2);
382  }
383  for (size_t k = 0; k < fv1.size(); k++) {
384  area2 += gkp_wa<UfResult>(level)[k] * (fv1[k] + fv2[k]);
385  if (calcabsasc) {
386  resabs += gkp_wa<UfResult>(level)[k] * (fabs(fv1[k]) + fabs(fv2[k]));
387  }
388  }
389  for (size_t k = 0; k < gkp_x<UfResult>(level).size(); k++) {
390  UfResult const abscissa = halfLength * gkp_x<UfResult>(level)[k];
391  UfResult const fval1 = func(center - abscissa);
392  UfResult const fval2 = func(center + abscissa);
393  UfResult const fval = fval1 + fval2;
394  area2 += gkp_wb<UfResult>(level)[k] * fval;
395  if (calcabsasc) {
396  resabs += gkp_wb<UfResult>(level)[k] * (fabs(fval1) + fabs(fval2));
397  }
398  fv1.push_back(fval1);
399  fv2.push_back(fval2);
400  if (fxmap) {
401  (*fxmap)[center - abscissa] = fval1;
402  (*fxmap)[center + abscissa] = fval2;
403  }
404  }
405 #ifdef COUNTFEVAL
406  nfeval += gkp_x<UfResult>(level).size() * 2;
407 #endif
408  if (calcabsasc) {
409  UfResult const mean = area1 * UfResult(0.5);
410  // resasc = approximation to the integral of abs(f-mean)
411  resasc = gkp_wb<UfResult>(level).back() * fabs(fCenter - mean);
412  for (size_t k = 0; k < gkp_wa<UfResult>(level).size(); k++) {
413  resasc += gkp_wa<UfResult>(level)[k] * (fabs(fv1[k] - mean) + fabs(fv2[k] - mean));
414  }
415  for (size_t k = 0; k < gkp_x<UfResult>(level).size(); k++) {
416  resasc += gkp_wb<UfResult>(level)[k] * (fabs(fv1[k] - mean) + fabs(fv2[k] - mean));
417  }
418  resasc *= absHalfLength;
419  resabs *= absHalfLength;
420  }
421  area2 *= halfLength;
422  err = rescale_error(fabs(area2 - area1), resabs, resasc);
423  if (err < resasc) {
424  calcabsasc = false;
425  }
426 
427  integ_dbg2 << "at level " << level << " area2 = " << area2;
428  integ_dbg2 << " +- " << err << std::endl;
429 
430  // test for convergence.
431  if (err < epsabs || err < epsrel * fabs(area2)) {
432  reg.SetArea(area2, err);
433  return true;
434  }
435  area1 = area2;
436  }
437 
438  // failed to converge
439  reg.SetArea(area1, err);
440 
441  integ_dbg2 << "Failed to reach tolerance with highest-order GKP rule";
442 
443  return false;
444 }
b
table::Key< int > b
Definition: TransmissionCurve.cc:467
std::endl
T endl(T... args)
a
table::Key< int > a
Definition: TransmissionCurve.cc:466
std::vector::push_back
T push_back(T... args)
std::fabs
T fabs(T... args)
std::vector::size
T size(T... args)
std::vector
STL class.
integ_dbg2
#define integ_dbg2
Definition: Integrate.h:183
lsst::afw::math::details::rescale_error
T rescale_error(T err, T const &resabs, T const &resasc)
Definition: Integrate.h:298
std::vector::reserve
T reserve(T... args)

◆ MinRep()

template<class T >
T lsst::afw::math::details::MinRep ( )
inline

Definition at line 282 of file Integrate.h.

282  {
283  return std::numeric_limits<T>::min();
284 }
std::numeric_limits::min
T min(T... args)

◆ norm()

template<class T >
T lsst::afw::math::details::norm ( const T &  x)
inline

Definition at line 194 of file Integrate.h.

194  {
195  return x * x;
196 }
x
double x
Definition: ChebyshevBoundedField.cc:277

◆ real()

template<class T >
T lsst::afw::math::details::real ( const T &  x)
inline

Definition at line 199 of file Integrate.h.

199  {
200  return x;
201 }
x
double x
Definition: ChebyshevBoundedField.cc:277

◆ rescale_error()

template<class T >
T lsst::afw::math::details::rescale_error ( T  err,
T const &  resabs,
T const &  resasc 
)
inline

Definition at line 298 of file Integrate.h.

298  {
299  if (resasc != 0.0 && err != 0.0) {
300  T const scale = (200.0 * err / resasc);
301  if (scale < 1.0) {
302  err = resasc * scale * sqrt(scale);
303  } else {
304  err = resasc;
305  }
306  }
307  if (resabs > MinRep<T>() / (50.0 * Epsilon<T>())) {
308  T const min_err = 50.0 * Epsilon<T>() * resabs;
309  if (min_err > err) {
310  err = min_err;
311  }
312  }
313  return err;
314 }
lsst::afw.display.ds9.scale
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:109
std::sqrt
T sqrt(T... args)
Generated on Wed Jan 8 2020 09:12:18 for LSSTApplications by   doxygen 1.8.13