TMB Documentation  v1.9.11
tiny_ad/integrate/integrate.hpp
1 #ifndef TINY_AD_INTEGRATE_H
2 #define TINY_AD_INTEGRATE_H
3 
4 /* Standalone ? */
5 #ifndef R_RCONFIG_H
6 #include <vector>
7 using std::vector;
8 #include <float.h> // INFINITY etc
9 #endif
10 
16 namespace gauss_kronrod {
17 // Fake that Rmath.h is included - and take explicitly what we need
18 // #ifndef RMATH_H
19 // #define RMATH_H
20 // #endif
21 // *** Update: Rmath.h inclusion now disabled in integrate.cpp
22 
23 /*
24  These fmin2 and fmax2 assume that derivatives are not required. This
25  is OK for integrate since these functions are used to calculate
26  error bounds - nothing else.
27 
28  FIXME: Put all helper functions for integrate in its own namespace
29  to ensure that fmin2, fmax2 etc, are never used elsewhere.
30 */
31 template<class T>
32 double value(T x){ return ((double*) &x)[0]; }
33 template<class S, class T> // FIXME: Would imin2(int,int) conflict with Rf_imin2 ?
34 int imin2(S x, T y) {return (x < y) ? x : y;}
35 template<class S, class T>
36 double fmin2(S x, T y) {return (value(x) < value(y)) ? value(x) : value(y);}
37 template<class S, class T>
38 double fmax2(S x, T y) {return (value(x) < value(y)) ? value(y) : value(x);}
39 
40 #include "integrate.cpp"
41 
49 struct control {
50  int subdivisions;
51  double reltol;
52  double abstol;
53  control(int subdivisions_ = 100,
54  double reltol_ = 1e-4,
55  double abstol_ = 1e-4) :
56  subdivisions(subdivisions_),
57  reltol(reltol_),
58  abstol(abstol_) {}
59 };
60 
75 template<class Integrand>
76 struct Integral {
77  typedef typename Integrand::Scalar Type;
78  // R's integrators require a vectorized integrand
79  struct vectorized_integrand {
80  Integrand f;
81  vectorized_integrand(Integrand f_) : f(f_) {}
82  void operator() (Type *x, int n, void *ex) {
83  for(int i=0; i<n; i++) x[i] = f(x[i]);
84  }
85  } fn;
87  Integrand& integrand(){return fn.f;}
88  // Input to R's integrators
89  Type epsabs, epsrel, result, abserr;
90  int neval, ier, limit, lenw, last;
91  vector<int> iwork;
92  vector<Type> work;
93  void setAccuracy(double epsrel_ = 1e-4, double epsabs_ = 1e-4) {
94  epsabs = epsabs_; epsrel = epsrel_; result = 0; abserr = 1e4;
95  neval = 0; ier = 0; last=0;
96  }
97  void setWorkspace(int subdivisions = 100) {
98  limit = subdivisions; lenw = 4 * limit;
99  iwork.resize(limit);
100  work.resize(lenw);
101  }
102  Type a, b, bound;
103  int inf;
104  void setBounds(Type a_, Type b_) {
105  int a_finite = (a_ != -INFINITY) && (a_ != INFINITY);
106  int b_finite = (b_ != -INFINITY) && (b_ != INFINITY);
107  if ( a_finite && b_finite) { inf = 0; a = a_; b = b_; }
108  else if ( a_finite && !b_finite) { inf = 1; bound = a_; }
109  else if (!a_finite && b_finite) { inf = -1; bound = b_; }
110  else { inf = 2; }
111  }
118  Integral(Integrand f_, Type a_, Type b_,
119  control c = control()
120  ) : fn(f_) {
121  setAccuracy(c.reltol, c.abstol);
122  setWorkspace(c.subdivisions);
123  setBounds(a_, b_);
124  }
125  Type operator()() {
126  if (inf)
127  Rdqagi(fn, NULL, &bound, &inf, &epsabs, &epsrel, &result, &abserr,
128  &neval, &ier, &limit, &lenw, &last, &iwork[0], &work[0]);
129  else
130  Rdqags(fn, NULL, &a, &b, &epsabs, &epsrel, &result, &abserr,
131  &neval, &ier, &limit, &lenw, &last, &iwork[0], &work[0]);
132  return result;
133  }
134 };
135 
162 template<class Integrand>
163 typename Integrand::Scalar integrate(Integrand f,
164  typename Integrand::Scalar a = -INFINITY,
165  typename Integrand::Scalar b = INFINITY,
166  control c = control() ) {
167  Integral< Integrand > I(f, a, b, c);
168  return I();
169 }
170 
185 template<class Integrand>
186 struct mvIntegral {
187  typedef typename Integrand::Scalar Scalar;
188  struct evaluator {
189  typedef typename Integrand::Scalar Scalar;
190  Integrand &f;
191  Scalar &x; // Integration variable
192  evaluator(Integrand &f_,
193  Scalar &x_ ) : f(f_), x(x_) {}
194  Scalar operator()(const Scalar &x_) {
195  x = x_;
196  return f();
197  }
198  } ev;
199  control c;
201  mvIntegral(Integrand &f_,
202  Scalar &x_,
203  Scalar a= -INFINITY,
204  Scalar b= INFINITY,
205  control c_ = control()) :
206  ev(f_, x_), c(c_), I(ev, a, b, c_) {}
207  Scalar operator()() { return I(); }
210  Scalar a= -INFINITY,
211  Scalar b= INFINITY) {
212  return mvIntegral<mvIntegral>(*this, x, a, b, c);
213  }
214 };
215 
216 
217 template<class Integrand>
218 struct mvIntegral0 {
219  typedef typename Integrand::Scalar Scalar;
220  Integrand &f;
221  control c;
222  mvIntegral0(Integrand &f_,
223  control c_) : f(f_), c(c_) {}
225  mvIntegral<Integrand> wrt(Scalar &x,
226  Scalar a= -INFINITY,
227  Scalar b= INFINITY) {
228  return mvIntegral<Integrand>(f, x, a, b, c);
229  }
230 };
258 template<class Integrand>
259 mvIntegral0<Integrand> mvIntegrate(Integrand &f,
260  control c = control() ) {
261  return mvIntegral0<Integrand> (f, c);
262 }
263 
264 } // End namespace gauss_kronrod
265 
266 // using gauss_kronrod::Integral;
267 // using gauss_kronrod::integrate;
268 // using gauss_kronrod::mvIntegrate;
269 
270 #endif
Vector class used by TMB.
Definition: vector.hpp:17
mvIntegral< mvIntegral > wrt(Scalar &x, Scalar a=-INFINITY, Scalar b=INFINITY)
With respect to.
Integral(Integrand f_, Type a_, Type b_, control c=control())
Constructor.
User control parameters for R&#39;s integrate.
Integrand::Scalar integrate(Integrand f, typename Integrand::Scalar a=-INFINITY, typename Integrand::Scalar b=INFINITY, control c=control())
Integrate function over finite or infinite interval.
mvIntegral0< Integrand > mvIntegrate(Integrand &f, control c=control())
Multivariate integration.
Interface to R&#39;s adaptive integrate routine.
Namespace with utility functions for adaptive numerical integration.
Multivariate integral class.
Integrand & integrand()
Return reference to integrand so the user can change parameters.
License: GPL v2