1 #ifndef TINY_AD_INTEGRATE_H 2 #define TINY_AD_INTEGRATE_H 32 double value(T x){
return ((
double*) &x)[0]; }
33 template<
class S,
class T>
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);}
40 #include "integrate.cpp" 53 control(
int subdivisions_ = 100,
54 double reltol_ = 1e-4,
55 double abstol_ = 1e-4) :
56 subdivisions(subdivisions_),
75 template<
class Integrand>
77 typedef typename Integrand::Scalar Type;
79 struct vectorized_integrand {
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]);
89 Type epsabs, epsrel, result, abserr;
90 int neval, ier, limit, lenw, last;
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;
97 void setWorkspace(
int subdivisions = 100) {
98 limit = subdivisions; lenw = 4 * limit;
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_; }
121 setAccuracy(c.reltol, c.abstol);
122 setWorkspace(c.subdivisions);
127 Rdqagi(fn, NULL, &bound, &inf, &epsabs, &epsrel, &result, &abserr,
128 &neval, &ier, &limit, &lenw, &last, &iwork[0], &work[0]);
130 Rdqags(fn, NULL, &a, &b, &epsabs, &epsrel, &result, &abserr,
131 &neval, &ier, &limit, &lenw, &last, &iwork[0], &work[0]);
162 template<
class Integrand>
164 typename Integrand::Scalar a = -INFINITY,
165 typename Integrand::Scalar b = INFINITY,
185 template<
class Integrand>
187 typedef typename Integrand::Scalar Scalar;
189 typedef typename Integrand::Scalar Scalar;
192 evaluator(Integrand &f_,
193 Scalar &x_ ) : f(f_), x(x_) {}
194 Scalar operator()(
const Scalar &x_) {
206 ev(f_, x_), c(c_), I(ev, a, b, c_) {}
207 Scalar operator()() {
return I(); }
211 Scalar b= INFINITY) {
217 template<
class Integrand>
219 typedef typename Integrand::Scalar Scalar;
222 mvIntegral0(Integrand &f_,
227 Scalar b= INFINITY) {
258 template<
class Integrand>
261 return mvIntegral0<Integrand> (f, c);
Vector class used by TMB.
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'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'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.