LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
Integrate.h
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2
3/*
4 * LSST Data Management System
5 * Copyright 2008, 2009, 2010 LSST Corporation.
6 *
7 * This product includes software developed by the
8 * LSST Project (http://www.lsst.org/).
9 *
10 * This program is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 3 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the LSST License Statement and
21 * the GNU General Public License along with this program. If not,
22 * see <http://www.lsstcorp.org/LegalNotices/>.
23 */
24
25#if !defined(LSST_AFW_MATH_INTEGRATE_H)
26#define LSST_AFW_MATH_INTEGRATE_H 1
27/*
28 * Compute 1d
29 */
30
31#include <algorithm>
32#include <cassert>
33#include <cmath>
34#include <complex>
35#include <functional>
36#include <limits>
37#include <map>
38#include <ostream>
39#include <queue>
40#include <sstream>
41#include <stdexcept>
42#include <vector>
43
44#include "lsst/pex/exceptions.h"
45
47
48// == The following is based on Mike Jarvis original comment ==
49//
50// Basic Usage:
51//
52// First, define a function object.
53// For example, to integrate a
54// Gaussian, use something along the lines of this:
55//
56// class Gauss {
57// public :
58// Gauss(double _mu, double _sig) : mu(_mu), sig(_sig) {}
59// double operator()(double x) const {
60// constexpr double inv_sqrt_2pi = 0.3989422804014327;
61// double a = (x - mu) / sig;
62// return inv_sqrt_2pi / sig * std::exp(-0.5 * a * a);
63// }
64// private :
65// double mu,sig;
66// };
67//
68// Next, make an IntRegion object with the bounds of the integration region.
69// You need to give it the type to use for the bounds and the value of the
70// functions (which need to be the same currently - some day I'll allow
71// complex functions...).
72//
73// For example, to integrate something from -1 to 1, use:
74//
75// lsst::afw::math::IntRegion<double> reg1(-1.,1.);
76// If a value is > 1.e10 or < -1.e10, then these values are taken to be
77// infinity, rather than the actual value.
78// So to integrate from 0 to infinity:
79//
80// lsst::afw::math::IntRegion<double> reg2(0.,20);
81//
82// Or, you can use the variable lsst::afw::math:MOCK_INF which might be clearer.
83// The integral might diverge depending on the limits chosen
84//
85// Finally, to perform the integral, the line would be:
86//
87// double integ1 = int1d(Gauss(0.,1.),reg1,1.e-10,1.e-4);
88// double integ2 = int1d(Gauss(0.,2.),reg2,1.e-10,1.e-4);
89//
90// which should yield 0.68 and 0.5 in our case.
91//
92// Those last two numbers indicate the precision required.
93// 1.e-10 is the required absolute error, and
94// 1.e-4 is the required relative error.
95//
96// If you want, you can omit these and 1.e-15,1.e-6 will be used as the
97// default precision (which are generally fine for most purposes).
98//
99// The absolute error only comes into play for results which are close to
100// 0 to prevent requiring an error of 0 for integrals which evaluate to 0
101// or very close to it.
102//
103//
104//
105// Advanced Usage:
106//
107// When an integration fails to converge with the usual GKP algorithm,
108// it splits the region into 2 (or more) and tries again with each sub-region.
109// The default is to just bisect the region (or something similarly smart for
110// infinite regions), but if you know of a good place to split the region,
111// you can tell it using:
112//
113// reg.AddSplit(10.)
114//
115// For example, if you know that you have a singularity somewhere in your
116// region, it would help the program a lot to split there, so you
117// should add that as a split point. Zeros can also be good choices.
118//
119// In addition to the integral being returned from int1d, int2d, or int3d as
120// the return value, the value is also stored in the region itself.
121// You can access it using:
122//
123// reg.Area();
124//
125// There is also an estimate of the error in the value:
126//
127// reg.Err();
128//
129// (It is intended to be an overestimate of the actual error,
130// but it doesn't always get it completely right.)
131
132
133namespace lsst {
134namespace afw {
135namespace math {
136
137double const MOCK_INF = 1.e10;
138
139#ifdef NDEBUG
140#define integ_dbg1 \
141 if (false) (*_dbgout)
142#define integ_dbg2 \
143 if (false) (*(reg.getDbgout()))
144#define integ_dbg3 \
145 if (false) (*(tempreg.getDbgout()))
146#else
147#define integ_dbg1 \
148 if (_dbgout) (*_dbgout)
149#define integ_dbg2 \
150 if (reg.getDbgout()) (*(reg.getDbgout()))
151#define integ_dbg3 \
152 if (tempreg.getDbgout()) (*(tempreg.getDbgout()))
153#endif
154
155//#define COUNTFEVAL
156// If defined, then count the number of function evaluations
157
158namespace details {
159template <class T>
160inline T norm(const T &x) {
161 return x * x;
162}
163using std::norm;
164template <class T>
165inline T real(const T &x) {
166 return x;
167}
168using std::real;
169#ifdef COUNTFEVAL
170int nfeval = 0;
171#endif
172} // namespace details
173
174template <class T>
175struct IntRegion final {
176public:
177 IntRegion(T const a, T const b, std::ostream *dbgout = nullptr)
178 : _a(a), _b(b), _error(0.0), _area(0), _dbgout(dbgout) {}
179
180 IntRegion(IntRegion const &) = default;
181 IntRegion(IntRegion &&) = default;
182 IntRegion &operator=(IntRegion const &) = default;
184 ~IntRegion() = default;
185
186 bool operator<(IntRegion<T> const &r2) const { return _error < r2._error; }
187 bool operator>(IntRegion<T> const &r2) const { return _error > r2._error; }
188
190 assert(children->size() == 0);
191 if (_splitpoints.size() == 0) {
192 Bisect();
193 }
194 if (_splitpoints.size() > 1) {
195 std::sort(_splitpoints.begin(), _splitpoints.end());
196 }
197
198#if 0
199 if (_a > _splitpoints[0] || _b < _splitpoints.back()) {
200 std::cerr << "a, b = " << _a << ', ' << _b << std::endl;
201 std::cerr << "_splitpoints = ";
202 for (size_t i = 0; i<_splitpoints.size(); i++) {
203 std::cerr << _splitpoints[i] << " ";
204 }
206 }
207#endif
208 assert(_splitpoints[0] >= _a);
209 assert(_splitpoints.back() <= _b);
210 children->push_back(IntRegion<T>(_a, _splitpoints[0], _dbgout));
211 for (size_t i = 1; i < _splitpoints.size(); i++) {
212 children->push_back(IntRegion<T>(_splitpoints[i - 1], _splitpoints[i], _dbgout));
213 }
214 children->push_back(IntRegion<T>(_splitpoints.back(), _b, _dbgout));
215 }
216
217 void Bisect() { _splitpoints.push_back((_a + _b) / 2.0); }
218 void AddSplit(const T x) { _splitpoints.push_back(x); }
219 size_t NSplit() const { return _splitpoints.size(); }
220
221 T const &Left() const { return _a; }
222 T const &Right() const { return _b; }
223 T const &Err() const { return _error; }
224 T const &Area() const { return _area; }
225 void SetArea(const T &a, const T &e) {
226 _area = a;
227 _error = e;
228 }
229
230 std::ostream *getDbgout() { return _dbgout; }
231
232private:
233 T _a, _b, _error, _area;
234 std::vector<T> _splitpoints;
235 std::ostream *_dbgout;
236};
237
238double const DEFABSERR = 1.e-15;
239double const DEFRELERR = 1.e-6;
240
241namespace details {
242
243template <class T>
244inline T Epsilon() {
246}
247template <class T>
248inline T MinRep() {
250}
251
252#ifdef EXTRA_PREC_H
253template <>
254inline Quad Epsilon<Quad>() {
255 return 3.08148791094e-33;
256}
257template <>
258inline Quad MinRep<Quad>() {
259 return 2.2250738585072014e-308;
260}
261#endif
262
263template <class T>
264inline T rescale_error(T err, T const &resabs, T const &resasc) {
265 if (resasc != 0.0 && err != 0.0) {
266 T const scale = (200.0 * err / resasc);
267 if (scale < 1.0) {
268 err = resasc * scale * sqrt(scale);
269 } else {
270 err = resasc;
271 }
272 }
273 if (resabs > MinRep<T>() / (50.0 * Epsilon<T>())) {
274 T const min_err = 50.0 * Epsilon<T>() * resabs;
275 if (min_err > err) {
276 err = min_err;
277 }
278 }
279 return err;
280}
281
297template <typename UnaryFunctionT, typename Arg>
298inline bool intGKPNA(UnaryFunctionT func, IntRegion<Arg> &reg,
299 Arg const epsabs, Arg const epsrel,
300 std::map<Arg, Arg> *fxmap = nullptr) {
301 Arg const a = reg.Left();
302 Arg const b = reg.Right();
303
304 Arg const halfLength = 0.5 * (b - a);
305 Arg const absHalfLength = fabs(halfLength);
306 Arg const center = 0.5 * (b + a);
307 Arg const fCenter = func(center);
308#ifdef COUNTFEVAL
309 nfeval++;
310#endif
311
312 assert(gkp_wb<Arg>(0).size() == gkp_x<Arg>(0).size() + 1);
313 Arg area1 = gkp_wb<Arg>(0).back() * fCenter;
314 std::vector<Arg> fv1, fv2;
315 fv1.reserve(2 * gkp_x<Arg>(0).size() + 1);
316 fv2.reserve(2 * gkp_x<Arg>(0).size() + 1);
317 for (size_t k = 0; k < gkp_x<Arg>(0).size(); k++) {
318 Arg const abscissa = halfLength * gkp_x<Arg>(0)[k];
319 Arg const fval1 = func(center - abscissa);
320 Arg const fval2 = func(center + abscissa);
321 area1 += gkp_wb<Arg>(0)[k] * (fval1 + fval2);
322 fv1.push_back(fval1);
323 fv2.push_back(fval2);
324 if (fxmap) {
325 (*fxmap)[center - abscissa] = fval1;
326 (*fxmap)[center + abscissa] = fval2;
327 }
328 }
329#ifdef COUNTFEVAL
330 nfeval += gkp_x<Arg>(0).size() * 2;
331#endif
332
333 integ_dbg2 << "level 0 rule: area = " << area1 << std::endl;
334
335 Arg err = 0;
336 bool calcabsasc = true;
337 Arg resabs = 0.0, resasc = 0.0;
338 for (int level = 1; level < NGKPLEVELS; level++) {
339 assert(gkp_wa<Arg>(level).size() == fv1.size());
340 assert(gkp_wa<Arg>(level).size() == fv2.size());
341 assert(gkp_wb<Arg>(level).size() == gkp_x<Arg>(level).size() + 1);
342 Arg area2 = gkp_wb<Arg>(level).back() * fCenter;
343 // resabs = approximation to integral of abs(f)
344 if (calcabsasc) {
345 resabs = fabs(area2);
346 }
347 for (size_t k = 0; k < fv1.size(); k++) {
348 area2 += gkp_wa<Arg>(level)[k] * (fv1[k] + fv2[k]);
349 if (calcabsasc) {
350 resabs += gkp_wa<Arg>(level)[k] * (fabs(fv1[k]) + fabs(fv2[k]));
351 }
352 }
353 for (size_t k = 0; k < gkp_x<Arg>(level).size(); k++) {
354 Arg const abscissa = halfLength * gkp_x<Arg>(level)[k];
355 Arg const fval1 = func(center - abscissa);
356 Arg const fval2 = func(center + abscissa);
357 Arg const fval = fval1 + fval2;
358 area2 += gkp_wb<Arg>(level)[k] * fval;
359 if (calcabsasc) {
360 resabs += gkp_wb<Arg>(level)[k] * (fabs(fval1) + fabs(fval2));
361 }
362 fv1.push_back(fval1);
363 fv2.push_back(fval2);
364 if (fxmap) {
365 (*fxmap)[center - abscissa] = fval1;
366 (*fxmap)[center + abscissa] = fval2;
367 }
368 }
369#ifdef COUNTFEVAL
370 nfeval += gkp_x<Arg>(level).size() * 2;
371#endif
372 if (calcabsasc) {
373 Arg const mean = area1 * Arg(0.5);
374 // resasc = approximation to the integral of abs(f-mean)
375 resasc = gkp_wb<Arg>(level).back() * fabs(fCenter - mean);
376 for (size_t k = 0; k < gkp_wa<Arg>(level).size(); k++) {
377 resasc += gkp_wa<Arg>(level)[k] * (fabs(fv1[k] - mean) + fabs(fv2[k] - mean));
378 }
379 for (size_t k = 0; k < gkp_x<Arg>(level).size(); k++) {
380 resasc += gkp_wb<Arg>(level)[k] * (fabs(fv1[k] - mean) + fabs(fv2[k] - mean));
381 }
382 resasc *= absHalfLength;
383 resabs *= absHalfLength;
384 }
385 area2 *= halfLength;
386 err = rescale_error(fabs(area2 - area1), resabs, resasc);
387 if (err < resasc) {
388 calcabsasc = false;
389 }
390
391 integ_dbg2 << "at level " << level << " area2 = " << area2;
392 integ_dbg2 << " +- " << err << std::endl;
393
394 // test for convergence.
395 if (err < epsabs || err < epsrel * fabs(area2)) {
396 reg.SetArea(area2, err);
397 return true;
398 }
399 area1 = area2;
400 }
401
402 // failed to converge
403 reg.SetArea(area1, err);
404
405 integ_dbg2 << "Failed to reach tolerance with highest-order GKP rule";
406
407 return false;
408}
409
425template <typename UnaryFunctionT, typename Arg>
426inline void intGKP(UnaryFunctionT func, IntRegion<Arg> &reg,
427 Arg const epsabs, Arg const epsrel,
428 std::map<Arg, Arg> *fxmap = nullptr) {
429 integ_dbg2 << "Start intGKP\n";
430
431 assert(epsabs >= 0.0);
432 assert(epsrel > 0.0);
433
434 // perform the first integration
435 bool done = intGKPNA(func, reg, epsabs, epsrel, fxmap);
436 if (done) return;
437
438 integ_dbg2 << "In adaptive GKP, failed first pass... subdividing\n";
439 integ_dbg2 << "Intial range = " << reg.Left() << ".." << reg.Right() << std::endl;
440
441 int roundoffType1 = 0, errorType = 0;
442 Arg roundoffType2 = 0;
443 size_t iteration = 1;
444
446 allregions.push(reg);
447 Arg finalarea = reg.Area();
448 Arg finalerr = reg.Err();
449 Arg tolerance = std::max(epsabs, epsrel * fabs(finalarea));
450 assert(finalerr > tolerance);
451
452 while (!errorType && finalerr > tolerance) {
453 // Bisect the subinterval with the largest error estimate
454 integ_dbg2 << "Current answer = " << finalarea << " +- " << finalerr;
455 integ_dbg2 << " (tol = " << tolerance << ")\n";
456 IntRegion<Arg> parent = allregions.top();
457 allregions.pop();
458 integ_dbg2 << "Subdividing largest error region ";
459 integ_dbg2 << parent.Left() << ".." << parent.Right() << std::endl;
460 integ_dbg2 << "parent area = " << parent.Area();
461 integ_dbg2 << " +- " << parent.Err() << std::endl;
463 parent.SubDivide(&children);
464 // For "GKP", there are only two, but for GKPOSC, there is one
465 // for each oscillation in region
466
467 // Try to do at least 3x better with the children
468 Arg factor = 3 * children.size() * finalerr / tolerance;
469 Arg newepsabs = fabs(parent.Err() / factor);
470 Arg newepsrel = newepsabs / fabs(parent.Area());
471 integ_dbg2 << "New epsabs, rel = " << newepsabs << ", " << newepsrel;
472 integ_dbg2 << " (" << children.size() << " children)\n";
473
474 Arg newarea = Arg(0.0);
475 Arg newerror = 0.0;
476 for (size_t i = 0; i < children.size(); i++) {
477 IntRegion<Arg> &child = children[i];
478 integ_dbg2 << "Integrating child " << child.Left();
479 integ_dbg2 << ".." << child.Right() << std::endl;
480 bool hasConverged;
481 hasConverged = intGKPNA(func, child, newepsabs, newepsrel);
482 integ_dbg2 << "child (" << i + 1 << '/' << children.size() << ") ";
483 if (hasConverged) {
484 integ_dbg2 << " converged.";
485 } else {
486 integ_dbg2 << " failed.";
487 }
488 integ_dbg2 << " Area = " << child.Area() << " +- " << child.Err() << std::endl;
489
490 newarea += child.Area();
491 newerror += child.Err();
492 }
493 integ_dbg2 << "Compare: newerr = " << newerror;
494 integ_dbg2 << " to parent err = " << parent.Err() << std::endl;
495
496 finalerr += (newerror - parent.Err());
497 finalarea += newarea - parent.Area();
498
499 Arg delta = parent.Area() - newarea;
500 if (newerror <= parent.Err() && fabs(delta) <= parent.Err() && newerror >= 0.99 * parent.Err()) {
501 integ_dbg2 << "roundoff type 1: delta/newarea = ";
502 integ_dbg2 << fabs(delta) / fabs(newarea);
503 integ_dbg2 << ", newerror/error = " << newerror / parent.Err() << std::endl;
504 roundoffType1++;
505 }
506 if (iteration >= 10 && newerror > parent.Err() && fabs(delta) <= newerror - parent.Err()) {
507 integ_dbg2 << "roundoff type 2: newerror/error = ";
508 integ_dbg2 << newerror / parent.Err() << std::endl;
509 roundoffType2 += std::min(newerror / parent.Err() - 1.0, Arg(1.0));
510 }
511
512 tolerance = std::max(epsabs, epsrel * fabs(finalarea));
513 if (finalerr > tolerance) {
514 if (roundoffType1 >= 200) {
515 errorType = 1; // round off error
516 integ_dbg2 << "GKP: Round off error 1\n";
517 }
518 if (roundoffType2 >= 200.0) {
519 errorType = 2; // round off error
520 integ_dbg2 << "GKP: Round off error 2\n";
521 }
522 if (fabs((parent.Right() - parent.Left()) / (reg.Right() - reg.Left())) < Epsilon<double>()) {
523 errorType = 3; // found singularity
524 integ_dbg2 << "GKP: Probable singularity\n";
525 }
526 }
527 for (size_t i = 0; i < children.size(); i++) {
528 allregions.push(children[i]);
529 }
530 iteration++;
531 }
532
533 // Recalculate finalarea in case there are any slight rounding errors
534 finalarea = 0.0;
535 finalerr = 0.0;
536 while (!allregions.empty()) {
537 IntRegion<Arg> const &r = allregions.top();
538 finalarea += r.Area();
539 finalerr += r.Err();
540 allregions.pop();
541 }
542 reg.SetArea(finalarea, finalerr);
543
544 if (errorType == 1) {
546 s << "Type 1 roundoff's = " << roundoffType1;
547 s << ", Type 2 = " << roundoffType2 << std::endl;
548 s << "Roundoff error 1 prevents tolerance from being achieved ";
549 s << "in intGKP\n";
551 } else if (errorType == 2) {
553 s << "Type 1 roundoff's = " << roundoffType1;
554 s << ", Type 2 = " << roundoffType2 << std::endl;
555 s << "Roundoff error 2 prevents tolerance from being achieved ";
556 s << "in intGKP\n";
558 } else if (errorType == 3) {
560 s << "Bad integrand behavior found in the integration interval ";
561 s << "in intGKP\n";
563 }
564}
565
570template <typename UnaryFunctionT>
571struct AuxFunc1 { // f(1/x-1) for int(a..infinity)
572public:
573 AuxFunc1(UnaryFunctionT const &f) : _f(f) {}
574 template <typename Arg>
575 auto operator()(Arg x) const {
576 return _f(1.0 / x - 1.0) / (x * x);
577 }
578private:
579 UnaryFunctionT const &_f;
580};
581
586template <class UF>
587AuxFunc1<UF> inline Aux1(UF uf) {
588 return AuxFunc1<UF>(uf);
589}
590
591template <typename UnaryFunctionT>
592struct AuxFunc2 { // f(1/x+1) for int(-infinity..b)
593public:
594 AuxFunc2(UnaryFunctionT const &f) : _f(f) {}
595 template <typename Arg>
596 auto operator()(Arg x) const {
597 return _f(1.0 / x + 1.0) / (x * x);
598 }
599private:
600 UnaryFunctionT const &_f;
601};
602
607template <class UF>
608AuxFunc2<UF> inline Aux2(UF uf) {
609 return AuxFunc2<UF>(uf);
610}
611
612
613// pulled from MoreFunctional.h. Needed in class Int2DAuxType and Int3DAuxType
614template <class BF>
616public:
617 binder2_1(const BF &oper, typename BF::first_argument_type val) : _oper(oper), _value(val) {}
618 typename BF::result_type operator()(const typename BF::second_argument_type &x) const {
619 return _oper(_value, x);
620 }
621
622protected:
624 typename BF::first_argument_type _value;
625};
626
627template <class BF, class Tp>
628inline binder2_1<BF> bind21(const BF &oper, const Tp &x) {
629 using Arg = typename BF::first_argument_type;
630 return binder2_1<BF>(oper, static_cast<Arg>(x));
631}
632
633template <class BF, class YREG>
635public:
636 Int2DAuxType(BF const &func, YREG const &yreg, typename BF::result_type const &abserr,
637 typename BF::result_type const &relerr)
638 : _func(func), _yreg(yreg), _abserr(abserr), _relerr(relerr) {}
639
640 typename BF::result_type operator()(typename BF::first_argument_type x) const {
641 typename YREG::result_type tempreg = _yreg(x);
642 typename BF::result_type result = int1d(bind21(_func, x), tempreg, _abserr, _relerr);
643 integ_dbg3 << "Evaluated int2dAux at x = " << x;
644 integ_dbg3 << ": f = " << result << " +- " << tempreg.Err() << std::endl;
645 return result;
646 }
647
648private:
649 BF const &_func;
650 YREG const &_yreg;
651 typename BF::result_type _abserr, _relerr;
652};
653
654// pulled from MoreFunctional.h. Needed in class Int3DAuxtype
655template <class TF>
657public:
658 binder3_1(const TF &oper, typename TF::firstof3_argument_type val) : _oper(oper), _value(val) {}
659 typename TF::result_type operator()(typename TF::secondof3_argument_type const &x1,
660 typename TF::thirdof3_argument_type const &x2) const {
661 return _oper(_value, x1, x2);
662 }
663
664protected:
666 typename TF::firstof3_argument_type _value;
667};
668
669template <class TF, class Tp>
670inline binder3_1<TF> bind31(const TF &oper, const Tp &x) {
671 using Arg = typename TF::firstof3_argument_type;
672 return binder3_1<TF>(oper, static_cast<Arg>(x));
673}
674
675template <class TF, class YREG, class ZREG>
677public:
678 Int3DAuxType(const TF &func, const YREG &yreg, const ZREG &zreg, const typename TF::result_type &abserr,
679 const typename TF::result_type &relerr)
680 : _func(func), _yreg(yreg), _zreg(zreg), _abserr(abserr), _relerr(relerr) {}
681
682 typename TF::result_type operator()(typename TF::firstof3_argument_type x) const {
683 typename YREG::result_type tempreg = _yreg(x);
684 typename TF::result_type result =
685 int2d(bind31(_func, x), tempreg, bind21(_zreg, x), _abserr, _relerr);
686 integ_dbg3 << "Evaluated int3dAux at x = " << x;
687 integ_dbg3 << ": f = " << result << " +- " << tempreg.Err() << std::endl;
688 return result;
689 }
690
691private:
692 const TF &_func;
693 const YREG &_yreg;
694 const ZREG &_zreg;
695 typename TF::result_type _abserr, _relerr;
696};
697
698} // end namespace details
699
703template <typename UnaryFunctionT, typename Arg>
704inline Arg int1d(UnaryFunctionT func, IntRegion<Arg> &reg,
705 Arg const &abserr = DEFABSERR,
706 Arg const &relerr = DEFRELERR) {
707 using namespace details;
708
709 integ_dbg2 << "start int1d: " << reg.Left() << ".." << reg.Right() << std::endl;
710
711 if ((reg.Left() <= -MOCK_INF && reg.Right() > 0) || (reg.Right() >= MOCK_INF && reg.Left() < 0)) {
712 reg.AddSplit(0);
713 }
714
715 if (reg.NSplit() > 0) {
717 reg.SubDivide(&children);
718 integ_dbg2 << "Subdivided into " << children.size() << " children\n";
719 Arg answer = Arg();
720 Arg err = 0;
721 for (size_t i = 0; i < children.size(); i++) {
722 IntRegion<Arg> &child = children[i];
723 integ_dbg2 << "i = " << i;
724 integ_dbg2 << ": bounds = " << child.Left() << ", " << child.Right() << std::endl;
725 answer += int1d(func, child, abserr, relerr);
726 err += child.Err();
727 integ_dbg2 << "subint = " << child.Area() << " +- " << child.Err() << std::endl;
728 }
729 reg.SetArea(answer, err);
730 return answer;
731
732 } else {
733 if (reg.Left() <= -MOCK_INF) {
734 integ_dbg2 << "left = -infinity, right = " << reg.Right() << std::endl;
735 assert(reg.Right() <= 0.0);
736 IntRegion<Arg> modreg(1.0 / (reg.Right() - 1.0), 0.0, reg.getDbgout());
737 intGKP(Aux2<UnaryFunctionT>(func), modreg, abserr, relerr);
738 reg.SetArea(modreg.Area(), modreg.Err());
739 } else if (reg.Right() >= MOCK_INF) {
740 integ_dbg2 << "left = " << reg.Left() << ", right = infinity\n";
741 assert(reg.Left() >= 0.0);
742 IntRegion<Arg> modreg(0.0, 1.0 / (reg.Left() + 1.0), reg.getDbgout());
743 intGKP(Aux1<UnaryFunctionT>(func), modreg, abserr, relerr);
744 reg.SetArea(modreg.Area(), modreg.Err());
745 } else {
746 integ_dbg2 << "left = " << reg.Left();
747 integ_dbg2 << ", right = " << reg.Right() << std::endl;
748 intGKP(func, reg, abserr, relerr);
749 }
750 integ_dbg2 << "done int1d answer = " << reg.Area();
751 integ_dbg2 << " +- " << reg.Err() << std::endl;
752 return reg.Area();
753 }
754}
755
756// =============================================================
764template <typename UnaryFunctionT, typename Arg>
765auto integrate(UnaryFunctionT func,
766 Arg const a,
767 Arg const b,
768 double eps = 1.0e-6) {
769 IntRegion<Arg> region(a, b);
770 return int1d(func, region, DEFABSERR, eps);
771}
772
773// =============================================================
778template <typename BinaryFunctionT, typename X, typename Y>
779auto integrate2d(BinaryFunctionT func, X x1, X x2, Y y1, Y y2, double eps = 1.0e-6) {
780 auto outer = [func, x1, x2, eps](auto y) {
781 auto inner = [func, y](auto x) { return func(x, y); };
782 return integrate(inner, x1, x2, eps);
783 };
784 return integrate(outer, y1, y2, eps);
785}
786} // namespace math
787} // namespace afw
788} // namespace lsst
789
790#endif
py::object result
Definition _schema.cc:429
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
#define integ_dbg2
Definition Integrate.h:149
#define integ_dbg3
Definition Integrate.h:151
int y
Definition SpanSet.cc:48
table::Key< int > b
table::Key< int > a
Int2DAuxType(BF const &func, YREG const &yreg, typename BF::result_type const &abserr, typename BF::result_type const &relerr)
Definition Integrate.h:636
BF::result_type operator()(typename BF::first_argument_type x) const
Definition Integrate.h:640
Int3DAuxType(const TF &func, const YREG &yreg, const ZREG &zreg, const typename TF::result_type &abserr, const typename TF::result_type &relerr)
Definition Integrate.h:678
TF::result_type operator()(typename TF::firstof3_argument_type x) const
Definition Integrate.h:682
BF::first_argument_type _value
Definition Integrate.h:624
BF::result_type operator()(const typename BF::second_argument_type &x) const
Definition Integrate.h:618
binder2_1(const BF &oper, typename BF::first_argument_type val)
Definition Integrate.h:617
TF::firstof3_argument_type _value
Definition Integrate.h:666
binder3_1(const TF &oper, typename TF::firstof3_argument_type val)
Definition Integrate.h:658
TF::result_type operator()(typename TF::secondof3_argument_type const &x1, typename TF::thirdof3_argument_type const &x2) const
Definition Integrate.h:659
Reports errors that are due to events beyond the control of the program.
Definition Runtime.h:104
T empty(T... args)
T endl(T... args)
T epsilon(T... args)
T max(T... args)
T min(T... args)
bool intGKPNA(UnaryFunctionT func, IntRegion< Arg > &reg, Arg const epsabs, Arg const epsrel, std::map< Arg, Arg > *fxmap=nullptr)
Non-adaptive integration of the function f over the region 'reg'.
Definition Integrate.h:298
binder3_1< TF > bind31(const TF &oper, const Tp &x)
Definition Integrate.h:670
void intGKP(UnaryFunctionT func, IntRegion< Arg > &reg, Arg const epsabs, Arg const epsrel, std::map< Arg, Arg > *fxmap=nullptr)
An adaptive integration algorithm which computes the integral of f over the region reg.
Definition Integrate.h:426
AuxFunc2< UF > Aux2(UF uf)
Auxiliary function 2.
Definition Integrate.h:608
T rescale_error(T err, T const &resabs, T const &resasc)
Definition Integrate.h:264
AuxFunc1< UF > Aux1(UF uf)
Auxiliary function 1.
Definition Integrate.h:587
binder2_1< BF > bind21(const BF &oper, const Tp &x)
Definition Integrate.h:628
double const DEFABSERR
Definition Integrate.h:238
double const DEFRELERR
Definition Integrate.h:239
double const MOCK_INF
Definition Integrate.h:137
Arg int1d(UnaryFunctionT func, IntRegion< Arg > &reg, Arg const &abserr=DEFABSERR, Arg const &relerr=DEFRELERR)
Front end for the 1d integrator.
Definition Integrate.h:704
auto integrate(UnaryFunctionT func, Arg const a, Arg const b, double eps=1.0e-6)
The 1D integrator.
Definition Integrate.h:765
auto integrate2d(BinaryFunctionT func, X x1, X x2, Y y1, Y y2, double eps=1.0e-6)
The 2D integrator.
Definition Integrate.h:779
T pop(T... args)
T push_back(T... args)
T push(T... args)
T reserve(T... args)
T size(T... args)
T sort(T... args)
ImageT val
Definition CR.cc:146
std::ostream * getDbgout()
Definition Integrate.h:230
void SetArea(const T &a, const T &e)
Definition Integrate.h:225
IntRegion(IntRegion &&)=default
T const & Err() const
Definition Integrate.h:223
bool operator>(IntRegion< T > const &r2) const
Definition Integrate.h:187
T const & Area() const
Definition Integrate.h:224
void AddSplit(const T x)
Definition Integrate.h:218
IntRegion & operator=(IntRegion &&)=default
IntRegion(T const a, T const b, std::ostream *dbgout=nullptr)
Definition Integrate.h:177
IntRegion(IntRegion const &)=default
void SubDivide(std::vector< IntRegion< T > > *children)
Definition Integrate.h:189
T const & Right() const
Definition Integrate.h:222
bool operator<(IntRegion< T > const &r2) const
Definition Integrate.h:186
T const & Left() const
Definition Integrate.h:221
IntRegion & operator=(IntRegion const &)=default
AuxFunc1(UnaryFunctionT const &f)
Definition Integrate.h:573
AuxFunc2(UnaryFunctionT const &f)
Definition Integrate.h:594
T top(T... args)