romberg.hpp

Go to the documentation of this file.

// Copyright (C) 2013-2015 Kasper Kristensen
// License: GPL-2

namespace romberg {
  template<class Type, class F>
  Type integrate(F f, Type a, Type b, int n = 7, int p = 2){
    Type e;
    return CppAD::RombergOne(f, a, b, n, p, e);
  }

  /* Construct 1D slice of vector function */
  template<class Type, class F>
  struct slicefun{
    F* f;
    vector<Type> x;
    int slice;
    slicefun(F* f_, vector<Type> x_, int slice_) {
      f=f_; x=x_; slice=slice_;
    }
    Type operator()(Type y){
      x(slice) = y;
      return f->operator()(x);
    }
  };

  /* integrate along i'th slice, f vector function */
  template<class Type, class F>
  Type integrate_slice(F &f, vector<Type> x, int i, Type a, Type b, int n=7, int p=2){
    return integrate(slicefun<Type, F>(&f, x, i), a, b, n, p);
  }

  template<class Type>
  struct multivariate_integrand{
    virtual Type operator()(vector<Type>) = 0;
  };

  template<class Type>
  struct integratelast : multivariate_integrand<Type>{
    slicefun<Type,multivariate_integrand<Type> >* fslice;
    int n, p;
    Type a, b;
    integratelast(multivariate_integrand<Type> &f, int dim, Type a_, Type b_, int n_ = 7, int p_ = 2){
      int slice = dim - 1;
      vector<Type> x(dim);
      fslice = new slicefun<Type,multivariate_integrand<Type> >(&f, x, slice);
      n=n_; p=p_;
      a=a_; b=b_;
    }
    ~integratelast(){delete fslice;}
    Type operator()(vector<Type> y){
      for(int i=0; i < y.size(); i++) fslice->x[i] = y[i];
      return integrate(*fslice, a, b, n, p);
    }
  };

  template<class Type>
  Type integrate_multi(multivariate_integrand<Type> &f, vector<Type> a, vector<Type> b, int n = 7, int p = 2){
    int dim = a.size();
    integratelast<Type> ffirst(f, dim, a(dim - 1), b(dim - 1), n, p);
    if(dim == 1){
      vector<Type> y(0);
      return ffirst(y);
    } else {
      vector<Type> afirst = a.segment(0, dim - 1);
      vector<Type> bfirst = b.segment(0, dim - 1);
      return integrate_multi(ffirst, afirst, bfirst, n, p);
    }
  }

  template<class Type, class F>
  Type integrate(F f, vector<Type> a, vector<Type> b, int n=7, int p=2){
    // Wrap f in container
    struct integrand : multivariate_integrand<Type>{
      F f;
      integrand(F f_) : f(f_) {}
      Type operator()(vector<Type> x){
        return f(x);
      }
    };
    integrand f_cpy(f);
    return integrate_multi(f_cpy, a, b, n, p);
  }
}