LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Classes | Functions
lsst::afw::math::details Namespace Reference

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.

685  {
686  typedef typename BF::first_argument_type Arg;
687  return binder2_1<BF>(oper, static_cast<Arg>(x));
688 }
x
int x
Definition: GaussianCentroid.cc:36
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.

737  {
738  typedef typename TF::firstof3_argument_type Arg;
739  return binder3_1<TF>(oper, static_cast<Arg>(x));
740 }
x
int x
Definition: GaussianCentroid.cc:36
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.

45  {
46  assert(level >= 0 && level < NGKPLEVELS);
47  static const int ngkp[NGKPLEVELS] = {5, 5, 11, 22, 44};
48  return ngkp[level];
49  }
template<class T >
const std::vector<T>& lsst::afw::math::details::gkp_wa ( int  level)
inline

Definition at line 173 of file IntGKPData10.h.

173  {
174 
175  // w21a, weights of the 21-point formula for abscissae x10
176  static const T aw21a[5] = {
177  0.032558162307964727478818972459390,
178  0.075039674810919952767043140916190,
179  0.109387158802297641899210590325805,
180  0.134709217311473325928054001771707,
181  0.147739104901338491374841515972068
182  };
183  static const std::vector<T> vw21a(aw21a, aw21a + 5);
184 
185  // w43a, weights of the 43-point formula for abscissae x10,x21
186  static const T aw43a[10] = {
187  0.016296734289666564924281974617663,
188  0.037522876120869501461613795898115,
189  0.054694902058255442147212685465005,
190  0.067355414609478086075553166302174,
191  0.073870199632393953432140695251367,
192  0.005768556059769796184184327908655,
193  0.027371890593248842081276069289151,
194  0.046560826910428830743339154433824,
195  0.061744995201442564496240336030883,
196  0.071387267268693397768559114425516
197  };
198  static const std::vector<T> vw43a(aw43a, aw43a + 10);
199 
200  // w87a, weights of the 87-point formula for abscissae x10..x43
201  static const T aw87a[21] = {
202  0.008148377384149172900002878448190,
203  0.018761438201562822243935059003794,
204  0.027347451050052286161582829741283,
205  0.033677707311637930046581056957588,
206  0.036935099820427907614589586742499,
207  0.002884872430211530501334156248695,
208  0.013685946022712701888950035273128,
209  0.023280413502888311123409291030404,
210  0.030872497611713358675466394126442,
211  0.035693633639418770719351355457044,
212  0.000915283345202241360843392549948,
213  0.005399280219300471367738743391053,
214  0.010947679601118931134327826856808,
215  0.016298731696787335262665703223280,
216  0.021081568889203835112433060188190,
217  0.025370969769253827243467999831710,
218  0.029189697756475752501446154084920,
219  0.032373202467202789685788194889595,
220  0.034783098950365142750781997949596,
221  0.036412220731351787562801163687577,
222  0.037253875503047708539592001191226
223  };
224  static const std::vector<T> vw87a(aw87a, aw87a + 21);
225 
226  // w175a, weights of the 175-point formula for abscissae x10..x87
227  static const T aw175a[43] = {
228  0.004074188692082600103872341,
229  0.009380719100781411819741177,
230  0.01367372552502614309257226,
231  0.01683885365581896502443513,
232  0.01846754991021395380770806,
233  0.001442436294030254510518689,
234  0.006842973011356375224430623,
235  0.01164020675144415562952859,
236  0.01543624880585667934076834,
237  0.01784681681970938536027839,
238  0.0004576754584154192941064591,
239  0.002699640110136963534043642,
240  0.005473839800559775445647306,
241  0.008149365848393670964180954,
242  0.01054078444460191775302900,
243  0.01268548488462691364854167,
244  0.01459484887823787625642303,
245  0.01618660123360139484467359,
246  0.01739154947518257137619173,
247  0.01820611036567589378188459,
248  0.01862693775152385427017140,
249  0.0001363171844264931156688042,
250  0.0009035606060432074452289352,
251  0.002048434635928091782250346,
252  0.003379145025868061092978424,
253  0.004774978836099393880602724,
254  0.006164723826122346701274159,
255  0.007505223673194467726814647,
256  0.008774483993121594091281642,
257  0.009969018893220443741650500,
258  0.01109746798050614328494604,
259  0.01216957356300040269315573,
260  0.01318725270741960360319771,
261  0.01414345539438560032188640,
262  0.01502629056404634765715107,
263  0.01582337568571996469999659,
264  0.01652520670998925164398162,
265  0.01712754985211303089259342,
266  0.01763120633007834051620255,
267  0.01803849481144435059221422,
268  0.01834930224922804724856518,
269  0.01856027463491628805666919,
270  0.01866711437596752016025132
271  };
272  static const std::vector<T> vw175a(aw175a, aw175a + 43);
273 
274  static const std::vector<T>* wa[NGKPLEVELS] = {0, &vw21a, &vw43a, &vw87a, &vw175a};
275 
276  assert(level >= 1 && level < NGKPLEVELS);
277  return *wa[level];
278  }
template<class T >
const std::vector<T>& lsst::afw::math::details::gkp_wb ( int  level)
inline

Definition at line 281 of file IntGKPData10.h.

281  {
282 
283  // w10, weights of the 10-point formula
284  static const T aw10b[6] = {
285  0.066671344308688137593568809893332,
286  0.149451349150580593145776339657697,
287  0.219086362515982043995534934228163,
288  0.269266719309996355091226921569469,
289  0.295524224714752870173892994651338,
290  0.0
291  };
292  static const std::vector<T> vw10b(aw10b, aw10b + 6);
293 
294  // w21b, weights of the 21-point formula for abscissae x21
295  static const T aw21b[6] = {
296  0.011694638867371874278064396062192,
297  0.054755896574351996031381300244580,
298  0.093125454583697605535065465083366,
299  0.123491976262065851077958109831074,
300  0.142775938577060080797094273138717,
301  0.149445554002916905664936468389821
302  };
303  static const std::vector<T> vw21b(aw21b, aw21b + 6);
304 
305  // w43b, weights of the 43-point formula for abscissae x43
306  static const T aw43b[12] = {
307  0.001844477640212414100389106552965,
308  0.010798689585891651740465406741293,
309  0.021895363867795428102523123075149,
310  0.032597463975345689443882222526137,
311  0.042163137935191811847627924327955,
312  0.050741939600184577780189020092084,
313  0.058379395542619248375475369330206,
314  0.064746404951445885544689259517511,
315  0.069566197912356484528633315038405,
316  0.072824441471833208150939535192842,
317  0.074507751014175118273571813842889,
318  0.074722147517403005594425168280423
319  };
320  static const std::vector<T> vw43b(aw43b, aw43b + 12);
321 
322  // w87b, weights of the 87-point formula for abscissae x87
323  static const T aw87b[23] = {
324  0.000274145563762072350016527092881,
325  0.001807124155057942948341311753254,
326  0.004096869282759164864458070683480,
327  0.006758290051847378699816577897424,
328  0.009549957672201646536053581325377,
329  0.012329447652244853694626639963780,
330  0.015010447346388952376697286041943,
331  0.017548967986243191099665352925900,
332  0.019938037786440888202278192730714,
333  0.022194935961012286796332102959499,
334  0.024339147126000805470360647041454,
335  0.026374505414839207241503786552615,
336  0.028286910788771200659968002987960,
337  0.030052581128092695322521110347341,
338  0.031646751371439929404586051078883,
339  0.033050413419978503290785944862689,
340  0.034255099704226061787082821046821,
341  0.035262412660156681033782717998428,
342  0.036076989622888701185500318003895,
343  0.036698604498456094498018047441094,
344  0.037120549269832576114119958413599,
345  0.037334228751935040321235449094698,
346  0.037361073762679023410321241766599
347  };
348  static const std::vector<T> vw87b(aw87b, aw87b + 23);
349 
350  // w175b, weights of the 175-point formula for abscissae x175
351  static const T aw175b[45] = {
352  0.00003901440664395060055484226,
353  0.0002787312407993348139961464,
354  0.0006672934146291138823512275,
355  0.001163051749910003055308177,
356  0.001738533916888715565420637,
357  0.002369555040550195269486165,
358  0.003036735036032612976379865,
359  0.003725394125884235586678819,
360  0.004424387541776355951933433,
361  0.005125062882078508655875582,
362  0.005820601038952832388647054,
363  0.006505667785357412495794539,
364  0.007176258976795165213169074,
365  0.007829642423806550517967042,
366  0.008464316525043036555640657,
367  0.009079917879101951192512711,
368  0.009677029304308397008134699,
369  0.01025687419430508685393322,
370  0.01082092935368957685380896,
371  0.01137052956179654620079332,
372  0.01190655150082017425918844,
373  0.01242924137407215963971061,
374  0.01293819980598456566872951,
375  0.01343248644726851711632966,
376  0.01391078107864287587682982,
377  0.01437154592198526807352501,
378  0.01481316256915705351874381,
379  0.01523404482122031244940129,
380  0.01563274056894636856297717,
381  0.01600802989372176552081826,
382  0.01635901148521841710605397,
383  0.01668515675646806540972517,
384  0.01698630892029276765404134,
385  0.01726261506199017604299177,
386  0.01751439915887340284018269,
387  0.01774200489411242260549908,
388  0.01794564958245060930824294,
389  0.01812532816637054898168673,
390  0.01828078919412312732671357,
391  0.01841158014800801587767163,
392  0.01851713825962505903167768,
393  0.01859689376712028821908313,
394  0.01865035760791320287947920,
395  0.01867717982648033250824110,
396  0.01868053688133951170552407
397  };
398  static const std::vector<T> vw175b(aw175b, aw175b + 45);
399 
400  static const std::vector<T>* wb[NGKPLEVELS] = {&vw10b, &vw21b, &vw43b, &vw87b, &vw175b};
401 
402  assert(level >= 0 && level < NGKPLEVELS);
403  return *wb[level];
404  }
template<class T >
const std::vector<T>& lsst::afw::math::details::gkp_x ( int  level)
inline

Definition at line 52 of file IntGKPData10.h.

52  {
53 
54  // x10, abscissae common to the 10-, 21-, 43- and 87-point rule
55  static const T ax10[5] = {
56  0.973906528517171720077964012084452,
57  0.865063366688984510732096688423493,
58  0.679409568299024406234327365114874,
59  0.433395394129247190799265943165784,
60  0.148874338981631210884826001129720
61  };
62  static const std::vector<T> vx10(ax10, ax10 + 5);
63 
64  // x21, abscissae common to the 21-, 43- and 87-point rule
65  static const T ax21[5] = {
66  0.995657163025808080735527280689003,
67  0.930157491355708226001207180059508,
68  0.780817726586416897063717578345042,
69  0.562757134668604683339000099272694,
70  0.294392862701460198131126603103866
71  };
72  static const std::vector<T> vx21(ax21, ax21 + 5);
73 
74  // x43, abscissae common to the 43- and 87-point rule
75  static const T ax43[11] = {
76  0.999333360901932081394099323919911,
77  0.987433402908088869795961478381209,
78  0.954807934814266299257919200290473,
79  0.900148695748328293625099494069092,
80  0.825198314983114150847066732588520,
81  0.732148388989304982612354848755461,
82  0.622847970537725238641159120344323,
83  0.499479574071056499952214885499755,
84  0.364901661346580768043989548502644,
85  0.222254919776601296498260928066212,
86  0.074650617461383322043914435796506
87  };
88  static const std::vector<T> vx43(ax43, ax43 + 11);
89 
90  // x87, abscissae of the 87-point rule
91  static const T ax87[22] = {
92  0.999902977262729234490529830591582,
93  0.997989895986678745427496322365960,
94  0.992175497860687222808523352251425,
95  0.981358163572712773571916941623894,
96  0.965057623858384619128284110607926,
97  0.943167613133670596816416634507426,
98  0.915806414685507209591826430720050,
99  0.883221657771316501372117548744163,
100  0.845710748462415666605902011504855,
101  0.803557658035230982788739474980964,
102  0.757005730685495558328942793432020,
103  0.706273209787321819824094274740840,
104  0.651589466501177922534422205016736,
105  0.593223374057961088875273770349144,
106  0.531493605970831932285268948562671,
107  0.466763623042022844871966781659270,
108  0.399424847859218804732101665817923,
109  0.329874877106188288265053371824597,
110  0.258503559202161551802280975429025,
111  0.185695396568346652015917141167606,
112  0.111842213179907468172398359241362,
113  0.037352123394619870814998165437704
114  };
115  static const std::vector<T> vx87(ax87, ax87 + 22);
116 
117  // x175, new abscissae of the 175-point rule
118  static const T ax175[44] = {
119  0.9999863601049729677997719,
120  0.9996988005421428623760900,
121  0.9987732473138838303043008,
122  0.9969583851613681945576954,
123  0.9940679446543307664316419,
124  0.9899673391741082769960501,
125  0.9845657252243253759726856,
126  0.9778061574856414566055911,
127  0.9696573134362640567828241,
128  0.9601075259508299668299148,
129  0.9491605191933640559839528,
130  0.9368321320434676982049391,
131  0.9231475266471484946315312,
132  0.9081385954543932335481458,
133  0.8918414554201179766275301,
134  0.8742940658668577971868106,
135  0.8555341249920432726468181,
136  0.8355974653331171075829366,
137  0.8145171513255796294567278,
138  0.7923233739697216224661771,
139  0.7690440708082055924983192,
140  0.7447060413364141221086362,
141  0.7193362474393392286563226,
142  0.6929630148748741413274435,
143  0.6656169582351898081906710,
144  0.6373315753624730450361940,
145  0.6081435407529742253347887,
146  0.5780927509986774074337714,
147  0.5472221576559162873479056,
148  0.5155774001694817879753989,
149  0.4832062516695473663839323,
150  0.4501579202236101928554315,
151  0.4164822936334372150371142,
152  0.3822292525735955063074986,
153  0.3474481818229070516859048,
154  0.3121877718977455143765043,
155  0.2764961311348905899785770,
156  0.2404211453027010909070665,
157  0.2040109589145005126324680,
158  0.1673144327616352753945014,
159  0.1303814601095591675431696,
160  0.09326308333660176874846400,
161  0.05601141598806823374241435,
162  0.01867941779948308845140053
163  };
164  static const std::vector<T> vx175(ax175, ax175 + 44);
165 
166  static const std::vector<T>* x[NGKPLEVELS] = {&vx10, &vx21, &vx43, &vx87, &vx175};
167 
168  assert(level >= 0 && level < NGKPLEVELS);
169  return *x[level];
170  }
x
int x
Definition: GaussianCentroid.cc:36
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.

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

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

280  {
281  if (resasc != 0.0 && err != 0.0) {
282  T const scale = (200.0 * err / resasc);
283  if (scale < 1.0) {
284  err = resasc * scale * sqrt(scale);
285  }
286  else {
287  err = resasc;
288  }
289  }
290  if (resabs > MinRep<T>() / (50.0 * Epsilon<T>())) {
291  T const min_err = 50.0 * Epsilon<T>() * resabs;
292  if (min_err > err) {
293  err = min_err;
294  }
295  }
296  return err;
297 }
Generated on Wed Sep 16 2015 13:35:45 for LSSTApplications by   doxygen 1.8.5