LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
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 and 2d integral
29  */
30 
31 #include <functional>
32 #include <vector>
33 #include <queue>
34 #include <map>
35 #include <cmath>
36 #include <algorithm>
37 #include <assert.h>
38 #include <limits>
39 #include <ostream>
40 #include <sstream>
41 #include <complex>
42 #include <stdexcept>
43 
44 #include "lsst/pex/exceptions.h"
45 
47 
48 // == The following is Mike Jarvis original comment ==
49 //
50 // Basic Usage:
51 //
52 // First, define a function object, which should derive from
53 // std::unary_function<double,double>. For example, to integrate a
54 // Gaussian, use something along the lines of this:
55 //
56 // class Gauss :
57 // public std::unary_function<double,double>
58 // {
59 // public :
60 //
61 // Gauss(double _mu, double _sig) :
62 // mu(_mu), sig(_sig), sigsq(_sig*_sig) {}
63 //
64 // double operator()(double x) const
65 // {
66 // const double SQRTTWOPI = 2.50662827463;
67 // return exp(-pow(x-mu,2)/2./sigsq)/SQRTTWOPI/sig;
68 // }
69 //
70 // private :
71 // double mu,sig,sigsq;
72 // };
73 //
74 // Next, make an IntRegion object with the bounds of the integration region.
75 // You need to give it the type to use for the bounds and the value of the
76 // functions (which need to be the same currently - some day I'll allow
77 // complex functions...).
78 //
79 // For example, to integrate something from -1 to 1, use:
80 //
81 // integ::IntRegion<double> reg1(-1.,1.);
82 //
83 // (The integ:: is because everything here is in the integ namespace to
84 // help prevent name conflicts with other header files.)
85 //
86 // If a value is > 1.e10 or < -1.e10, then these values are taken to be
87 // infinity, rather than the actual value.
88 // So to integrate from 0 to infinity:
89 //
90 // integ::IntRegion<double> reg2(0.,1.e100);
91 //
92 // Or, you can use the variable integ::MOCK_INF which might be clearer.
93 //
94 //
95 // Finally, to perform the integral, the line would be:
96 //
97 // double integ1 = int1d(Gauss(0.,1.),reg1,1.e-10,1.e-4);
98 // double integ2 = int1d(Gauss(0.,2.),reg2,1.e-10,1.e-4);
99 //
100 // which should yield 0.68 and 0.5 in our case.
101 //
102 // Those last two numbers indicate the precision required.
103 // 1.e-10 is the required absolute error, and
104 // 1.e-4 is the required relative error.
105 //
106 // If you want, you can omit these and 1.e-15,1.e-6 will be used as the
107 // default precision (which are generally fine for most purposes).
108 //
109 // The absolute error only comes into play for results which are close to
110 // 0 to prevent requiring an error of 0 for integrals which evaluate to 0
111 // or very close to it.
112 //
113 //
114 //
115 // Advanced Usage:
116 //
117 // When an integration fails to converge with the usual GKP algorithm,
118 // it splits the region into 2 (or more) and tries again with each sub-region.
119 // The default is to just bisect the region (or something similarly smart for
120 // infinite regions), but if you know of a good place to split the region,
121 // you can tell it using:
122 //
123 // reg.AddSplit(10.)
124 //
125 // For example, if you know that you have a singularity somewhere in your
126 // region, it would help the program a lot to split there, so you
127 // should add that as a split point. Zeros can also be good choices.
128 //
129 // In addition to the integral being returned from int1d, int2d, or int3d as
130 // the return value, the value is also stored in the region itself.
131 // You can access it using:
132 //
133 // reg.Area();
134 //
135 // There is also an estimate of the error in the value:
136 //
137 // reg.Err();
138 //
139 // (It is intended to be an overestimate of the actual error,
140 // but it doesn't always get it completely right.)
141 //
142 //
143 //
144 // Two- and Three-Dimensional Integrals:
145 //
146 // These are slightly more complicated. The easiest case is when the
147 // bounds of the integral are a rectangle or 3d box. In this case,
148 // you can still use the regular IntRegion. The only new thing then
149 // is the definition of the function. For example, to integrate
150 // int(3x^2 + xy + y , x=0..1, y=0..1):
151 //
152 // struct Integrand :
153 // public std::binary_function<double,double,double>
154 // {
155 // double operator()(double x, double y) const
156 // { return x*(3.*x + y) + y; }
157 // };
158 //
159 // integ::IntRegion<double> reg3(0.,1.);
160 // double integ3 = int2d(Integrand(),reg3,reg3);
161 //
162 // (Which should give 1.75 as the result.)
163 //
164 //
165 //
166 
167 namespace lsst {
168 namespace afw {
169 namespace math {
170 
171 double const MOCK_INF = 1.e10;
172 
173 #ifdef NDEBUG
174 #define integ_dbg1 \
175  if (false) (*_dbgout)
176 #define integ_dbg2 \
177  if (false) (*(reg.getDbgout()))
178 #define integ_dbg3 \
179  if (false) (*(tempreg.getDbgout()))
180 #else
181 #define integ_dbg1 \
182  if (_dbgout) (*_dbgout)
183 #define integ_dbg2 \
184  if (reg.getDbgout()) (*(reg.getDbgout()))
185 #define integ_dbg3 \
186  if (tempreg.getDbgout()) (*(tempreg.getDbgout()))
187 #endif
188 
189 //#define COUNTFEVAL
190 // If defined, then count the number of function evaluations
191 
192 namespace details {
193 template <class T>
194 inline T norm(const T &x) {
195  return x * x;
196 }
197 using std::norm;
198 template <class T>
199 inline T real(const T &x) {
200  return x;
201 }
202 using std::real;
203 #ifdef COUNTFEVAL
204 int nfeval = 0;
205 #endif
206 } // namespace details
207 
208 template <class T>
209 struct IntRegion final {
210 public:
211  IntRegion(T const a, T const b, std::ostream *dbgout = 0)
212  : _a(a), _b(b), _error(0.0), _area(0), _dbgout(dbgout) {}
213 
214  IntRegion(IntRegion const &) = default;
215  IntRegion(IntRegion &&) = default;
216  IntRegion &operator=(IntRegion const &) = default;
217  IntRegion &operator=(IntRegion &&) = default;
218  ~IntRegion() = default;
219 
220  bool operator<(IntRegion<T> const &r2) const { return _error < r2._error; }
221  bool operator>(IntRegion<T> const &r2) const { return _error > r2._error; }
222 
223  void SubDivide(std::vector<IntRegion<T> > *children) {
224  assert(children->size() == 0);
225  if (_splitpoints.size() == 0) {
226  Bisect();
227  }
228  if (_splitpoints.size() > 1) {
229  std::sort(_splitpoints.begin(), _splitpoints.end());
230  }
231 
232 #if 0
233  if (_a > _splitpoints[0] || _b < _splitpoints.back()) {
234  std::cerr << "a, b = " << _a << ', ' << _b << std::endl;
235  std::cerr << "_splitpoints = ";
236  for (size_t i = 0; i<_splitpoints.size(); i++) {
237  std::cerr << _splitpoints[i] << " ";
238  }
239  std::cerr << std::endl;
240  }
241 #endif
242  assert(_splitpoints[0] >= _a);
243  assert(_splitpoints.back() <= _b);
244  children->push_back(IntRegion<T>(_a, _splitpoints[0], _dbgout));
245  for (size_t i = 1; i < _splitpoints.size(); i++) {
246  children->push_back(IntRegion<T>(_splitpoints[i - 1], _splitpoints[i], _dbgout));
247  }
248  children->push_back(IntRegion<T>(_splitpoints.back(), _b, _dbgout));
249  }
250 
251  void Bisect() { _splitpoints.push_back((_a + _b) / 2.0); }
252  void AddSplit(const T x) { _splitpoints.push_back(x); }
253  size_t NSplit() const { return _splitpoints.size(); }
254 
255  T const &Left() const { return _a; }
256  T const &Right() const { return _b; }
257  T const &Err() const { return _error; }
258  T const &Area() const { return _area; }
259  void SetArea(const T &a, const T &e) {
260  _area = a;
261  _error = e;
262  }
263 
264  std::ostream *getDbgout() { return _dbgout; }
265 
266 private:
267  T _a, _b, _error, _area;
268  std::vector<T> _splitpoints;
269  std::ostream *_dbgout;
270 };
271 
272 double const DEFABSERR = 1.e-15;
273 double const DEFRELERR = 1.e-6;
274 
275 namespace details {
276 
277 template <class T>
278 inline T Epsilon() {
280 }
281 template <class T>
282 inline T MinRep() {
284 }
285 
286 #ifdef EXTRA_PREC_H
287 template <>
288 inline Quad Epsilon<Quad>() {
289  return 3.08148791094e-33;
290 }
291 template <>
292 inline Quad MinRep<Quad>() {
293  return 2.2250738585072014e-308;
294 }
295 #endif
296 
297 template <class T>
298 inline T rescale_error(T err, T const &resabs, T const &resasc) {
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 }
315 
332 template <class UF>
333 inline bool intGKPNA(UF const &func, IntRegion<typename UF::result_type> &reg,
334  typename UF::result_type const epsabs, typename UF::result_type const epsrel,
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 }
445 
462 template <class UF>
463 inline void intGKP(UF const &func, IntRegion<typename UF::result_type> &reg,
464  typename UF::result_type const epsabs, typename UF::result_type const epsrel,
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 
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;
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) {
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";
589  } else if (errorType == 2) {
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";
596  } else if (errorType == 3) {
598  s << "Bad integrand behavior found in the integration interval ";
599  s << "in intGKP\n";
601  }
602 }
603 
608 template <class UF>
609 struct AuxFunc1 : // f(1/x-1) for int(a..infinity)
610  public std::unary_function<typename UF::argument_type, typename UF::result_type> {
611 public:
612  AuxFunc1(const UF &f) : _f(f) {}
613  typename UF::result_type operator()(typename UF::argument_type x) const {
614  return _f(1.0 / x - 1.0) / (x * x);
615  }
616 
617 private:
618  UF const &_f;
619 };
620 
625 template <class UF>
626 AuxFunc1<UF> inline Aux1(UF uf) {
627  return AuxFunc1<UF>(uf);
628 }
629 
630 template <class UF>
631 struct AuxFunc2 : // f(1/x+1) for int(-infinity..b)
632  public std::unary_function<typename UF::argument_type, typename UF::result_type> {
633 public:
634  AuxFunc2(UF const &f) : _f(f) {}
635  typename UF::result_type operator()(typename UF::argument_type x) const {
636  return _f(1.0 / x + 1.0) / (x * x);
637  }
638 
639 private:
640  UF const &_f;
641 };
642 
647 template <class UF>
648 AuxFunc2<UF> inline Aux2(UF uf) {
649  return AuxFunc2<UF>(uf);
650 }
651 
656 template <class T>
657 struct ConstantReg1 : public std::unary_function<T, IntRegion<T> > {
658  ConstantReg1(T a, T b) : ir(a, b) {}
659  ConstantReg1(IntRegion<T> const &r) : ir(r) {}
660  IntRegion<T> operator()(T) const { return ir; }
662 };
663 
664 template <class T>
665 struct ConstantReg2 : public std::binary_function<T, T, IntRegion<T> > {
666  ConstantReg2(T a, T b) : ir(a, b) {}
667  ConstantReg2(IntRegion<T> const &r) : ir(r) {}
668  IntRegion<T> operator()(T x, T y) const { return ir; }
670 };
671 
672 // pulled from MoreFunctional.h. Needed in class Int2DAuxType and Int3DAuxType
673 template <class BF>
674 class binder2_1 : public std::unary_function<typename BF::second_argument_type, typename BF::result_type> {
675 public:
676  binder2_1(const BF &oper, typename BF::first_argument_type val) : _oper(oper), _value(val) {}
677  typename BF::result_type operator()(const typename BF::second_argument_type &x) const {
678  return _oper(_value, x);
679  }
680 
681 protected:
682  BF _oper;
683  typename BF::first_argument_type _value;
684 };
685 
686 template <class BF, class Tp>
687 inline binder2_1<BF> bind21(const BF &oper, const Tp &x) {
688  typedef typename BF::first_argument_type Arg;
689  return binder2_1<BF>(oper, static_cast<Arg>(x));
690 }
691 
692 template <class BF, class YREG>
693 class Int2DAuxType : public std::unary_function<typename BF::first_argument_type, typename BF::result_type> {
694 public:
695  Int2DAuxType(BF const &func, YREG const &yreg, typename BF::result_type const &abserr,
696  typename BF::result_type const &relerr)
697  : _func(func), _yreg(yreg), _abserr(abserr), _relerr(relerr) {}
698 
699  typename BF::result_type operator()(typename BF::first_argument_type x) const {
700  typename YREG::result_type tempreg = _yreg(x);
701  typename BF::result_type result = int1d(bind21(_func, x), tempreg, _abserr, _relerr);
702  integ_dbg3 << "Evaluated int2dAux at x = " << x;
703  integ_dbg3 << ": f = " << result << " +- " << tempreg.Err() << std::endl;
704  return result;
705  }
706 
707 private:
708  BF const &_func;
709  YREG const &_yreg;
710  typename BF::result_type _abserr, _relerr;
711 };
712 
713 // pulled from MoreFunctional.h. Needed in class Int3DAuxtype
714 template <class TF>
715 class binder3_1 : public std::binary_function<typename TF::secondof3_argument_type,
716  typename TF::thirdof3_argument_type, typename TF::result_type> {
717 public:
718  binder3_1(const TF &oper, typename TF::firstof3_argument_type val) : _oper(oper), _value(val) {}
719  typename TF::result_type operator()(typename TF::secondof3_argument_type const &x1,
720  typename TF::thirdof3_argument_type const &x2) const {
721  return _oper(_value, x1, x2);
722  }
723 
724 protected:
725  TF _oper;
726  typename TF::firstof3_argument_type _value;
727 };
728 
729 template <class TF, class Tp>
730 inline binder3_1<TF> bind31(const TF &oper, const Tp &x) {
731  typedef typename TF::firstof3_argument_type Arg;
732  return binder3_1<TF>(oper, static_cast<Arg>(x));
733 }
734 
735 template <class TF, class YREG, class ZREG>
737  : public std::unary_function<typename TF::firstof3_argument_type, typename TF::result_type> {
738 public:
739  Int3DAuxType(const TF &func, const YREG &yreg, const ZREG &zreg, const typename TF::result_type &abserr,
740  const typename TF::result_type &relerr)
741  : _func(func), _yreg(yreg), _zreg(zreg), _abserr(abserr), _relerr(relerr) {}
742 
743  typename TF::result_type operator()(typename TF::firstof3_argument_type x) const {
744  typename YREG::result_type tempreg = _yreg(x);
745  typename TF::result_type result =
746  int2d(bind31(_func, x), tempreg, bind21(_zreg, x), _abserr, _relerr);
747  integ_dbg3 << "Evaluated int3dAux at x = " << x;
748  integ_dbg3 << ": f = " << result << " +- " << tempreg.Err() << std::endl;
749  return result;
750  }
751 
752 private:
753  const TF &_func;
754  const YREG &_yreg;
755  const ZREG &_zreg;
756  typename TF::result_type _abserr, _relerr;
757 };
758 
759 } // end namespace details
760 
764 template <class UF>
765 inline typename UF::result_type int1d(UF const &func, IntRegion<typename UF::result_type> &reg,
766  typename UF::result_type const &abserr = DEFABSERR,
767  typename UF::result_type const &relerr = DEFRELERR) {
768  typedef typename UF::result_type UfResult;
769  using namespace details;
770 
771  integ_dbg2 << "start int1d: " << reg.Left() << ".." << reg.Right() << std::endl;
772 
773  if ((reg.Left() <= -MOCK_INF && reg.Right() > 0) || (reg.Right() >= MOCK_INF && reg.Left() < 0)) {
774  reg.AddSplit(0);
775  }
776 
777  if (reg.NSplit() > 0) {
779  reg.SubDivide(&children);
780  integ_dbg2 << "Subdivided into " << children.size() << " children\n";
781  UfResult answer = UfResult();
782  UfResult err = 0;
783  for (size_t i = 0; i < children.size(); i++) {
784  IntRegion<UfResult> &child = children[i];
785  integ_dbg2 << "i = " << i;
786  integ_dbg2 << ": bounds = " << child.Left() << ", " << child.Right() << std::endl;
787  answer += int1d(func, child, abserr, relerr);
788  err += child.Err();
789  integ_dbg2 << "subint = " << child.Area() << " +- " << child.Err() << std::endl;
790  }
791  reg.SetArea(answer, err);
792  return answer;
793 
794  } else {
795  if (reg.Left() <= -MOCK_INF) {
796  integ_dbg2 << "left = -infinity, right = " << reg.Right() << std::endl;
797  assert(reg.Right() <= 0.0);
798  IntRegion<UfResult> modreg(1.0 / (reg.Right() - 1.0), 0.0, reg.getDbgout());
799  intGKP(Aux2<UF>(func), modreg, abserr, relerr);
800  reg.SetArea(modreg.Area(), modreg.Err());
801  } else if (reg.Right() >= MOCK_INF) {
802  integ_dbg2 << "left = " << reg.Left() << ", right = infinity\n";
803  assert(reg.Left() >= 0.0);
804  IntRegion<UfResult> modreg(0.0, 1.0 / (reg.Left() + 1.0), reg.getDbgout());
805  intGKP(Aux1<UF>(func), modreg, abserr, relerr);
806  reg.SetArea(modreg.Area(), modreg.Err());
807  } else {
808  integ_dbg2 << "left = " << reg.Left();
809  integ_dbg2 << ", right = " << reg.Right() << std::endl;
810  intGKP(func, reg, abserr, relerr);
811  }
812  integ_dbg2 << "done int1d answer = " << reg.Area();
813  integ_dbg2 << " +- " << reg.Err() << std::endl;
814  return reg.Area();
815  }
816 }
817 
821 template <class BF, class YREG>
822 inline typename BF::result_type int2d(BF const &func, IntRegion<typename BF::result_type> &reg,
823  YREG const &yreg, typename BF::result_type const &abserr = DEFABSERR,
824  typename BF::result_type const &relerr = DEFRELERR) {
825  using namespace details;
826  integ_dbg2 << "Starting int2d: range = ";
827  integ_dbg2 << reg.Left() << ".." << reg.Right() << std::endl;
828  Int2DAuxType<BF, YREG> faux(func, yreg, abserr * 1.0e-3, relerr * 1.0e-3);
829  typename BF::result_type answer = int1d(faux, reg, abserr, relerr);
830  integ_dbg2 << "done int2d answer = " << answer << " +- " << reg.Err() << std::endl;
831  return answer;
832 }
833 
837 template <class TF, class YREG, class ZREG>
838 inline typename TF::result_type int3d(TF const &func, IntRegion<typename TF::result_type> &reg,
839  YREG const &yreg, ZREG const &zreg,
840  typename TF::result_type const &abserr = DEFABSERR,
841  typename TF::result_type const &relerr = DEFRELERR) {
842  using namespace details;
843  integ_dbg2 << "Starting int3d: range = ";
844  integ_dbg2 << reg.Left() << ".." << reg.Right() << std::endl;
845  Int3DAuxType<TF, YREG, ZREG> faux(func, yreg, zreg, abserr * 1.e-3, relerr * 1.e-3);
846  typename TF::result_type answer = int1d(faux, reg, abserr, relerr);
847  integ_dbg2 << "done int3d answer = " << answer << " +- " << reg.Err() << std::endl;
848  return answer;
849 }
850 
854 template <class BF>
855 inline typename BF::result_type int2d(BF const &func, IntRegion<typename BF::result_type> &reg,
857  typename BF::result_type const &abserr = DEFABSERR,
858  typename BF::result_type const &relerr = DEFRELERR) {
859  using namespace details;
860  return int2d(func, reg, ConstantReg1<typename BF::result_type>(yreg), abserr, relerr);
861 }
862 
866 template <class TF>
867 inline typename TF::result_type int3d(TF const &func, IntRegion<typename TF::result_type> &reg,
870  typename TF::result_type const &abserr = DEFABSERR,
871  typename TF::result_type const &relerr = DEFRELERR) {
872  using namespace details;
873  return int3d(func, reg, ConstantReg1<typename TF::result_type>(yreg),
874  ConstantReg2<typename TF::result_type>(zreg), abserr, relerr);
875 }
876 
877 // =============================================================
885 template <typename UnaryFunctionT>
886 typename UnaryFunctionT::result_type integrate(UnaryFunctionT func,
887  typename UnaryFunctionT::argument_type const a,
888  typename UnaryFunctionT::argument_type const b,
889  double eps = 1.0e-6) {
890  typedef typename UnaryFunctionT::argument_type Arg;
891  IntRegion<Arg> region(a, b);
892 
893  return int1d(func, region, DEFABSERR, eps);
894 }
895 
896 namespace details {
897 
907 template <typename BinaryFunctionT>
908 class FunctionWrapper : public std::unary_function<typename BinaryFunctionT::second_argument_type,
909  typename BinaryFunctionT::result_type> {
910 public:
911  FunctionWrapper(BinaryFunctionT func, typename BinaryFunctionT::first_argument_type const x1,
912  typename BinaryFunctionT::first_argument_type const x2, double const eps = 1.0e-6)
913  : _func(func), _x1(x1), _x2(x2), _eps(eps) {}
914  typename BinaryFunctionT::result_type operator()(
915  typename BinaryFunctionT::second_argument_type const y) const {
916  return integrate(std::bind2nd(_func, y), _x1, _x2, _eps);
917  }
918 
919 private:
920  BinaryFunctionT _func;
921  typename BinaryFunctionT::first_argument_type _x1, _x2;
922  double _eps;
923 };
924 
925 } // namespace details
926 
927 // =============================================================
934 template <typename BinaryFunctionT>
935 typename BinaryFunctionT::result_type integrate2d(BinaryFunctionT func,
936  typename BinaryFunctionT::first_argument_type const x1,
937  typename BinaryFunctionT::first_argument_type const x2,
938  typename BinaryFunctionT::second_argument_type const y1,
939  typename BinaryFunctionT::second_argument_type const y2,
940  double eps = 1.0e-6) {
941  using namespace details;
942  // note the more stringent eps requirement to ensure the requested limit
943  // can be reached.
944  FunctionWrapper<BinaryFunctionT> fwrap(func, x1, x2, eps);
945  return integrate(fwrap, y1, y2, eps);
946 }
947 } // namespace math
948 } // namespace afw
949 } // namespace lsst
950 
951 #endif
py::object result
Definition: _schema.cc:430
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
#define integ_dbg2
Definition: Integrate.h:183
#define integ_dbg3
Definition: Integrate.h:185
int y
Definition: SpanSet.cc:49
table::Key< int > b
table::Key< int > a
Wrap an integrand in a call to a 1D integrator: romberg()
Definition: Integrate.h:909
FunctionWrapper(BinaryFunctionT func, typename BinaryFunctionT::first_argument_type const x1, typename BinaryFunctionT::first_argument_type const x2, double const eps=1.0e-6)
Definition: Integrate.h:911
BinaryFunctionT::result_type operator()(typename BinaryFunctionT::second_argument_type const y) const
Definition: Integrate.h:914
Int2DAuxType(BF const &func, YREG const &yreg, typename BF::result_type const &abserr, typename BF::result_type const &relerr)
Definition: Integrate.h:695
BF::result_type operator()(typename BF::first_argument_type x) const
Definition: Integrate.h:699
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:739
TF::result_type operator()(typename TF::firstof3_argument_type x) const
Definition: Integrate.h:743
BF::first_argument_type _value
Definition: Integrate.h:683
BF::result_type operator()(const typename BF::second_argument_type &x) const
Definition: Integrate.h:677
binder2_1(const BF &oper, typename BF::first_argument_type val)
Definition: Integrate.h:676
TF::firstof3_argument_type _value
Definition: Integrate.h:726
binder3_1(const TF &oper, typename TF::firstof3_argument_type val)
Definition: Integrate.h:718
TF::result_type operator()(typename TF::secondof3_argument_type const &x1, typename TF::thirdof3_argument_type const &x2) const
Definition: Integrate.h:719
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:648
AuxFunc1< UF > Aux1(UF uf)
Auxiliary function 1.
Definition: Integrate.h:626
T norm(const T &x)
Definition: Integrate.h:194
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.
Definition: Integrate.h:463
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'.
Definition: Integrate.h:333
T real(const T &x)
Definition: Integrate.h:199
T rescale_error(T err, T const &resabs, T const &resasc)
Definition: Integrate.h:298
binder2_1< BF > bind21(const BF &oper, const Tp &x)
Definition: Integrate.h:687
binder3_1< TF > bind31(const TF &oper, const Tp &x)
Definition: Integrate.h:730
double const DEFABSERR
Definition: Integrate.h:272
double const DEFRELERR
Definition: Integrate.h:273
BinaryFunctionT::result_type integrate2d(BinaryFunctionT func, typename BinaryFunctionT::first_argument_type const x1, typename BinaryFunctionT::first_argument_type const x2, typename BinaryFunctionT::second_argument_type const y1, typename BinaryFunctionT::second_argument_type const y2, double eps=1.0e-6)
The 2D integrator.
Definition: Integrate.h:935
double const MOCK_INF
Definition: Integrate.h:171
BF::result_type int2d(BF const &func, IntRegion< typename BF::result_type > &reg, YREG const &yreg, typename BF::result_type const &abserr=DEFABSERR, typename BF::result_type const &relerr=DEFRELERR)
Front end for the 2d integrator.
Definition: Integrate.h:822
TF::result_type int3d(TF const &func, IntRegion< typename TF::result_type > &reg, YREG const &yreg, ZREG const &zreg, typename TF::result_type const &abserr=DEFABSERR, typename TF::result_type const &relerr=DEFRELERR)
Front end for the 3d integrator.
Definition: Integrate.h:838
UF::result_type int1d(UF const &func, IntRegion< typename UF::result_type > &reg, typename UF::result_type const &abserr=DEFABSERR, typename UF::result_type const &relerr=DEFRELERR)
Front end for the 1d integrator.
Definition: Integrate.h:765
UnaryFunctionT::result_type integrate(UnaryFunctionT func, typename UnaryFunctionT::argument_type const a, typename UnaryFunctionT::argument_type const b, double eps=1.0e-6)
The 1D integrator.
Definition: Integrate.h:886
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:259
IntRegion(IntRegion &&)=default
bool operator>(IntRegion< T > const &r2) const
Definition: Integrate.h:221
void AddSplit(const T x)
Definition: Integrate.h:252
IntRegion(T const a, T const b, std::ostream *dbgout=0)
Definition: Integrate.h:211
T const & Area() const
Definition: Integrate.h:258
size_t NSplit() const
Definition: Integrate.h:253
std::ostream * getDbgout()
Definition: Integrate.h:264
T const & Left() const
Definition: Integrate.h:255
T const & Right() const
Definition: Integrate.h:256
IntRegion(IntRegion const &)=default
void SubDivide(std::vector< IntRegion< T > > *children)
Definition: Integrate.h:223
bool operator<(IntRegion< T > const &r2) const
Definition: Integrate.h:220
IntRegion & operator=(IntRegion const &)=default
IntRegion & operator=(IntRegion &&)=default
T const & Err() const
Definition: Integrate.h:257
UF::result_type operator()(typename UF::argument_type x) const
Definition: Integrate.h:613
UF::result_type operator()(typename UF::argument_type x) const
Definition: Integrate.h:635
Helpers for constant regions for int2d, int3d:
Definition: Integrate.h:657
IntRegion< T > operator()(T) const
Definition: Integrate.h:660
ConstantReg1(IntRegion< T > const &r)
Definition: Integrate.h:659
ConstantReg2(IntRegion< T > const &r)
Definition: Integrate.h:667
IntRegion< T > operator()(T x, T y) const
Definition: Integrate.h:668
T top(T... args)