LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
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 
133 namespace lsst {
134 namespace afw {
135 namespace math {
136 
137 double 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 
158 namespace details {
159 template <class T>
160 inline T norm(const T &x) {
161  return x * x;
162 }
163 using std::norm;
164 template <class T>
165 inline T real(const T &x) {
166  return x;
167 }
168 using std::real;
169 #ifdef COUNTFEVAL
170 int nfeval = 0;
171 #endif
172 } // namespace details
173 
174 template <class T>
175 struct IntRegion final {
176 public:
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;
183  IntRegion &operator=(IntRegion &&) = 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 
189  void SubDivide(std::vector<IntRegion<T> > *children) {
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  }
205  std::cerr << std::endl;
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 
232 private:
233  T _a, _b, _error, _area;
234  std::vector<T> _splitpoints;
235  std::ostream *_dbgout;
236 };
237 
238 double const DEFABSERR = 1.e-15;
239 double const DEFRELERR = 1.e-6;
240 
241 namespace details {
242 
243 template <class T>
244 inline T Epsilon() {
246 }
247 template <class T>
248 inline T MinRep() {
250 }
251 
252 #ifdef EXTRA_PREC_H
253 template <>
254 inline Quad Epsilon<Quad>() {
255  return 3.08148791094e-33;
256 }
257 template <>
258 inline Quad MinRep<Quad>() {
259  return 2.2250738585072014e-308;
260 }
261 #endif
262 
263 template <class T>
264 inline 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 
297 template <typename UnaryFunctionT, typename Arg>
298 inline 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 
425 template <typename UnaryFunctionT, typename Arg>
426 inline 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;
462  std::vector<IntRegion<Arg> > children;
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 
570 template <typename UnaryFunctionT>
571 struct AuxFunc1 { // f(1/x-1) for int(a..infinity)
572 public:
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  }
578 private:
579  UnaryFunctionT const &_f;
580 };
581 
586 template <class UF>
587 AuxFunc1<UF> inline Aux1(UF uf) {
588  return AuxFunc1<UF>(uf);
589 }
590 
591 template <typename UnaryFunctionT>
592 struct AuxFunc2 { // f(1/x+1) for int(-infinity..b)
593 public:
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  }
599 private:
600  UnaryFunctionT const &_f;
601 };
602 
607 template <class UF>
608 AuxFunc2<UF> inline Aux2(UF uf) {
609  return AuxFunc2<UF>(uf);
610 }
611 
612 
613 // pulled from MoreFunctional.h. Needed in class Int2DAuxType and Int3DAuxType
614 template <class BF>
615 class binder2_1 : public std::unary_function<typename BF::second_argument_type, typename BF::result_type> {
616 public:
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 
622 protected:
623  BF _oper;
624  typename BF::first_argument_type _value;
625 };
626 
627 template <class BF, class Tp>
628 inline 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 
633 template <class BF, class YREG>
634 class Int2DAuxType : public std::unary_function<typename BF::first_argument_type, typename BF::result_type> {
635 public:
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 
648 private:
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
655 template <class TF>
656 class binder3_1 : public std::binary_function<typename TF::secondof3_argument_type,
657  typename TF::thirdof3_argument_type, typename TF::result_type> {
658 public:
659  binder3_1(const TF &oper, typename TF::firstof3_argument_type val) : _oper(oper), _value(val) {}
660  typename TF::result_type operator()(typename TF::secondof3_argument_type const &x1,
661  typename TF::thirdof3_argument_type const &x2) const {
662  return _oper(_value, x1, x2);
663  }
664 
665 protected:
666  TF _oper;
667  typename TF::firstof3_argument_type _value;
668 };
669 
670 template <class TF, class Tp>
671 inline binder3_1<TF> bind31(const TF &oper, const Tp &x) {
672  using Arg = typename TF::firstof3_argument_type;
673  return binder3_1<TF>(oper, static_cast<Arg>(x));
674 }
675 
676 template <class TF, class YREG, class ZREG>
678  : public std::unary_function<typename TF::firstof3_argument_type, typename TF::result_type> {
679 public:
680  Int3DAuxType(const TF &func, const YREG &yreg, const ZREG &zreg, const typename TF::result_type &abserr,
681  const typename TF::result_type &relerr)
682  : _func(func), _yreg(yreg), _zreg(zreg), _abserr(abserr), _relerr(relerr) {}
683 
684  typename TF::result_type operator()(typename TF::firstof3_argument_type x) const {
685  typename YREG::result_type tempreg = _yreg(x);
686  typename TF::result_type result =
687  int2d(bind31(_func, x), tempreg, bind21(_zreg, x), _abserr, _relerr);
688  integ_dbg3 << "Evaluated int3dAux at x = " << x;
689  integ_dbg3 << ": f = " << result << " +- " << tempreg.Err() << std::endl;
690  return result;
691  }
692 
693 private:
694  const TF &_func;
695  const YREG &_yreg;
696  const ZREG &_zreg;
697  typename TF::result_type _abserr, _relerr;
698 };
699 
700 } // end namespace details
701 
705 template <typename UnaryFunctionT, typename Arg>
706 inline Arg int1d(UnaryFunctionT func, IntRegion<Arg> &reg,
707  Arg const &abserr = DEFABSERR,
708  Arg const &relerr = DEFRELERR) {
709  using namespace details;
710 
711  integ_dbg2 << "start int1d: " << reg.Left() << ".." << reg.Right() << std::endl;
712 
713  if ((reg.Left() <= -MOCK_INF && reg.Right() > 0) || (reg.Right() >= MOCK_INF && reg.Left() < 0)) {
714  reg.AddSplit(0);
715  }
716 
717  if (reg.NSplit() > 0) {
718  std::vector<IntRegion<Arg> > children;
719  reg.SubDivide(&children);
720  integ_dbg2 << "Subdivided into " << children.size() << " children\n";
721  Arg answer = Arg();
722  Arg err = 0;
723  for (size_t i = 0; i < children.size(); i++) {
724  IntRegion<Arg> &child = children[i];
725  integ_dbg2 << "i = " << i;
726  integ_dbg2 << ": bounds = " << child.Left() << ", " << child.Right() << std::endl;
727  answer += int1d(func, child, abserr, relerr);
728  err += child.Err();
729  integ_dbg2 << "subint = " << child.Area() << " +- " << child.Err() << std::endl;
730  }
731  reg.SetArea(answer, err);
732  return answer;
733 
734  } else {
735  if (reg.Left() <= -MOCK_INF) {
736  integ_dbg2 << "left = -infinity, right = " << reg.Right() << std::endl;
737  assert(reg.Right() <= 0.0);
738  IntRegion<Arg> modreg(1.0 / (reg.Right() - 1.0), 0.0, reg.getDbgout());
739  intGKP(Aux2<UnaryFunctionT>(func), modreg, abserr, relerr);
740  reg.SetArea(modreg.Area(), modreg.Err());
741  } else if (reg.Right() >= MOCK_INF) {
742  integ_dbg2 << "left = " << reg.Left() << ", right = infinity\n";
743  assert(reg.Left() >= 0.0);
744  IntRegion<Arg> modreg(0.0, 1.0 / (reg.Left() + 1.0), reg.getDbgout());
745  intGKP(Aux1<UnaryFunctionT>(func), modreg, abserr, relerr);
746  reg.SetArea(modreg.Area(), modreg.Err());
747  } else {
748  integ_dbg2 << "left = " << reg.Left();
749  integ_dbg2 << ", right = " << reg.Right() << std::endl;
750  intGKP(func, reg, abserr, relerr);
751  }
752  integ_dbg2 << "done int1d answer = " << reg.Area();
753  integ_dbg2 << " +- " << reg.Err() << std::endl;
754  return reg.Area();
755  }
756 }
757 
758 // =============================================================
766 template <typename UnaryFunctionT, typename Arg>
767 auto integrate(UnaryFunctionT func,
768  Arg const a,
769  Arg const b,
770  double eps = 1.0e-6) {
771  IntRegion<Arg> region(a, b);
772  return int1d(func, region, DEFABSERR, eps);
773 }
774 
775 // =============================================================
780 template <typename BinaryFunctionT, typename X, typename Y>
781 auto integrate2d(BinaryFunctionT func, X x1, X x2, Y y1, Y y2, double eps = 1.0e-6) {
782  auto outer = [func, x1, x2, eps](auto y) {
783  auto inner = [func, y](auto x) { return func(x, y); };
784  return integrate(inner, x1, x2, eps);
785  };
786  return integrate(outer, y1, y2, eps);
787 }
788 } // namespace math
789 } // namespace afw
790 } // namespace lsst
791 
792 #endif
py::object result
Definition: _schema.cc:429
double x
#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:680
TF::result_type operator()(typename TF::firstof3_argument_type x) const
Definition: Integrate.h:684
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:667
binder3_1(const TF &oper, typename TF::firstof3_argument_type val)
Definition: Integrate.h:659
TF::result_type operator()(typename TF::secondof3_argument_type const &x1, typename TF::thirdof3_argument_type const &x2) const
Definition: Integrate.h:660
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)
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:108
class[[deprecated("Removed with no replacement (but see lsst::afw::image::TransmissionCurve). Will be " "removed after v22.")]] FilterProperty final
Describe the properties of a Filter (e.g.
Definition: Filter.h:53
AuxFunc2< UF > Aux2(UF uf)
Auxiliary function 2.
Definition: Integrate.h:608
AuxFunc1< UF > Aux1(UF uf)
Auxiliary function 1.
Definition: Integrate.h:587
T norm(const T &x)
Definition: Integrate.h:160
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
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
T real(const T &x)
Definition: Integrate.h:165
T rescale_error(T err, T const &resabs, T const &resasc)
Definition: Integrate.h:264
binder2_1< BF > bind21(const BF &oper, const Tp &x)
Definition: Integrate.h:628
binder3_1< TF > bind31(const TF &oper, const Tp &x)
Definition: Integrate.h:671
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:706
auto integrate(UnaryFunctionT func, Arg const a, Arg const b, double eps=1.0e-6)
The 1D integrator.
Definition: Integrate.h:767
auto integrate2d(BinaryFunctionT func, X x1, X x2, Y y1, Y y2, double eps=1.0e-6)
The 2D integrator.
Definition: Integrate.h:781
A base class for image defects.
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
T str(T... args)
void SetArea(const T &a, const T &e)
Definition: Integrate.h:225
IntRegion(IntRegion &&)=default
bool operator>(IntRegion< T > const &r2) const
Definition: Integrate.h:187
void AddSplit(const T x)
Definition: Integrate.h:218
T const & Area() const
Definition: Integrate.h:224
size_t NSplit() const
Definition: Integrate.h:219
std::ostream * getDbgout()
Definition: Integrate.h:230
T const & Left() const
Definition: Integrate.h:221
IntRegion(T const a, T const b, std::ostream *dbgout=nullptr)
Definition: Integrate.h:177
T const & Right() const
Definition: Integrate.h:222
IntRegion(IntRegion const &)=default
void SubDivide(std::vector< IntRegion< T > > *children)
Definition: Integrate.h:189
bool operator<(IntRegion< T > const &r2) const
Definition: Integrate.h:186
IntRegion & operator=(IntRegion const &)=default
IntRegion & operator=(IntRegion &&)=default
T const & Err() const
Definition: Integrate.h:223
AuxFunc1(UnaryFunctionT const &f)
Definition: Integrate.h:573
AuxFunc2(UnaryFunctionT const &f)
Definition: Integrate.h:594
T top(T... args)