tiny_ad/integrate/integrate.hpp

#ifndef TINY_AD_INTEGRATE_H
#define TINY_AD_INTEGRATE_H

/* Standalone ? */
#ifndef R_RCONFIG_H
#include <vector>
using std::vector;
#include <float.h>  // INFINITY etc
#endif

namespace gauss_kronrod {
// Fake that Rmath.h is included - and take explicitly what we need
// #ifndef RMATH_H
// #define RMATH_H
// #endif
// *** Update: Rmath.h inclusion now disabled in integrate.cpp

/*
  These fmin2 and fmax2 assume that derivatives are not required. This
  is OK for integrate since these functions are used to calculate
  error bounds - nothing else.

  FIXME: Put all helper functions for integrate in its own namespace
  to ensure that fmin2, fmax2 etc, are never used elsewhere.
*/
template<class T>
double value(T x){ return ((double*) &x)[0]; }
template<class S, class T> // FIXME: Would imin2(int,int) conflict with Rf_imin2 ?
int imin2(S x, T y) {return (x < y) ? x : y;}
template<class S, class T>
double fmin2(S x, T y) {return (value(x) < value(y)) ? value(x) : value(y);}
template<class S, class T>
double fmax2(S x, T y) {return (value(x) < value(y)) ? value(y) : value(x);}

#include "integrate.cpp"

struct control {
  int subdivisions;
  double reltol;
  double abstol;
  control(int subdivisions_ = 100,
          double reltol_    = 1e-4,
          double abstol_    = 1e-4) :
    subdivisions(subdivisions_),
    reltol(reltol_),
    abstol(abstol_) {}
};

template<class Integrand>
struct Integral {
  typedef typename Integrand::Scalar Type;
  // R's integrators require a vectorized integrand
  struct vectorized_integrand {
    Integrand f;
    vectorized_integrand(Integrand f_) : f(f_) {}
    void operator() (Type *x, int n, void *ex) {
      for(int i=0; i<n; i++) x[i] = f(x[i]);
    }
  } fn;
  Integrand& integrand(){return fn.f;}
  // Input to R's integrators
  Type epsabs, epsrel, result, abserr;
  int neval, ier, limit, lenw, last;
  vector<int> iwork;
  vector<Type> work;
  void setAccuracy(double epsrel_ = 1e-4, double epsabs_ = 1e-4) {
    epsabs = epsabs_; epsrel = epsrel_; result = 0; abserr = 1e4;
    neval = 0; ier = 0; last=0;
  }
  void setWorkspace(int subdivisions = 100) {
     limit = subdivisions; lenw = 4 * limit;
     iwork.resize(limit);
     work.resize(lenw);
  }
  Type a, b, bound;
  int inf;
  void setBounds(Type a_, Type b_) {
    int a_finite = (a_ != -INFINITY) && (a_ != INFINITY);
    int b_finite = (b_ != -INFINITY) && (b_ != INFINITY);
    if      ( a_finite &&  b_finite) { inf =  0; a = a_; b = b_; }
    else if ( a_finite && !b_finite) { inf =  1; bound = a_; }
    else if (!a_finite &&  b_finite) { inf = -1; bound = b_; }
    else                             { inf =  2; }
  }
  Integral(Integrand f_, Type a_, Type b_,
           control c = control()
           ) : fn(f_) {
    setAccuracy(c.reltol, c.abstol);
    setWorkspace(c.subdivisions);
    setBounds(a_, b_);
  }
  Type operator()() {
    if (inf)
      Rdqagi(fn, NULL, &bound, &inf, &epsabs, &epsrel, &result, &abserr,
             &neval, &ier, &limit, &lenw, &last, &iwork[0], &work[0]);
    else
      Rdqags(fn, NULL, &a, &b, &epsabs, &epsrel, &result, &abserr,
             &neval, &ier, &limit, &lenw, &last, &iwork[0], &work[0]);
    return result;
  }
};

template<class Integrand>
typename Integrand::Scalar integrate(Integrand f,
                                     typename Integrand::Scalar a = -INFINITY,
                                     typename Integrand::Scalar b =  INFINITY,
                                     control c = control() ) {
  Integral< Integrand > I(f, a, b, c);
  return I();
}

template<class Integrand>
struct mvIntegral {
  typedef typename Integrand::Scalar Scalar;
  struct evaluator {
    typedef typename Integrand::Scalar Scalar;
    Integrand &f;
    Scalar &x; // Integration variable
    evaluator(Integrand &f_,
              Scalar &x_ ) : f(f_), x(x_) {}
    Scalar operator()(const Scalar &x_) {
      x = x_;
      return f();
    }
  } ev;
  control c;
  Integral<evaluator> I;
  mvIntegral(Integrand &f_,
             Scalar &x_,
             Scalar a= -INFINITY,
             Scalar b=  INFINITY,
             control c_ = control()) :
    ev(f_, x_), c(c_), I(ev, a, b, c_) {}
  Scalar operator()() { return I(); }
  mvIntegral<mvIntegral> wrt(Scalar &x,
                             Scalar a= -INFINITY,
                             Scalar b=  INFINITY) {
    return mvIntegral<mvIntegral>(*this, x, a, b, c);
  }
};


template<class Integrand>
struct mvIntegral0 {
  typedef typename Integrand::Scalar Scalar;
  Integrand &f;
  control c;
  mvIntegral0(Integrand &f_,
              control c_) : f(f_), c(c_) {}
  mvIntegral<Integrand> wrt(Scalar &x,
                            Scalar a= -INFINITY,
                            Scalar b=  INFINITY) {
    return mvIntegral<Integrand>(f, x, a, b, c);
  }
};
template<class Integrand>
mvIntegral0<Integrand> mvIntegrate(Integrand &f,
                                   control c = control() ) {
  return mvIntegral0<Integrand> (f, c);
}

} // End namespace gauss_kronrod

// using gauss_kronrod::Integral;
// using gauss_kronrod::integrate;
// using gauss_kronrod::mvIntegrate;

#endif