lsoda

A C++ header library for using the libsoda-cxx C++ library with R

https://github.com/mclements/lsoda

Repository

A C++ header library for using the libsoda-cxx C++ library with R

README.org

* lsoda: C++ header library for ordinary differential equations

** Summary

- This R package provides a C++ header file for using the ~lsoda~ function from ~ODEPACK~. This package fills a gap to allow ~LSODA::ode()~ to be called from C++ in other packages.
- We provide an example R function ~lsoda::ode()~, which is similar to ~deSolve::lsoda()~. The former function implements the update of the state vector ~y~ and allows for calculated values, but does not implement events. This function is a proof-of-principle and there is no reason to prefer this function to the one from ~deSolve~.
- C++ signatures for ~LSODA::ode()~:

#+begin_src Cpp :exports code :eval yes
  typedef void (*LSODA_ODE_SYSTEM_TYPE)(double t, double *y, double *dydt, void *);

  template
  Rcpp::NumericMatrix ode(Vector y,
			  Vector times,
			  LSODA_ODE_SYSTEM_TYPE func,
			  size_t nout = 0, // default value => y.size()
			  void* data = (void*) nullptr,
			  double rtol=1e-6, double atol = 1e-6);
  
  template
  Rcpp::NumericMatrix ode(Vector y,
			  Vector times,
			  Functor functor,
			  double rtol=1e-6, double atol = 1e-6)
#+end_src

- For the first ~ode()~ signature, the template argument ~Vector~ can include ~std::vector~, ~arma::vec~ or ~Eigen::VectorXd~. The variable ~nout~ should be the length of the ~dydt~ from ~func~, which is a vector of derivatives (first) /and possibly any other results/.
- For the second signature, the functor takes the signature ~Vector(double,Vector)~ for template type ~Vector~ (again including ~std::vector~, ~arma::vec~ or ~Eigen::VectorXd~). The functor can be a class with a method ~(Vector)operator()(double,Vector)~, or a function. The functor should return a Vector of derivatives (first) /and possibly any other results/. 
- The ~LSODA::ode()~ functions return a ~Rcpp::NumericMatrix~ with the first column being time, then with columns for the state vectors, then with columns for the other results.
- Note: The C++ version of ~lsoda~ is a different version to the Fortran code  used in ~deSolve~ and ~scipy~. The differences are minor.


** Installation

- The only dependency is ~Rcpp~.
- Install the R package from GitHub. The package has been submitted to CRAN.

#+begin_src R :session *R* :exports code :eval yes
  devtools::install_github("mclements/lsoda")
#+end_src


** Examples

*** Example: set-up and calling from R

- We will use the ~Rcpp~, ~lsoda~ and ~deSolve~ packages.
- For the examples, we will use specific times and initial values.
- We first show how to call the ordinary differential equation solver from R and compare the results with those from ~deSolve::lsoda()~. 

#+begin_src R :session *R* :results output :exports both :eval yes
  library(Rcpp)
  library(lsoda)
  library(deSolve)
  times = c(0,0.4*10^(0:10))
  y = c(1,0,0)
  func = function(t,y,...) {
      ydot = rep(0,3)
      ydot[1] = 1.0E4 * y[2] * y[3] - .04E0 * y[1]
      ydot[3] = 3.0E7 * y[2] * y[2]
      ydot[2] = -1.0 * (ydot[1] + ydot[3])
      list(ydot, sum(y))
  }
  lsoda1 = lsoda::ode(y,times,func, rtol=1e-8, atol=1e-8)
  deSolve1 = deSolve::ode(y,times, func, rtol=1e-8, atol=1e-8)
  range(lsoda1 - deSolve1)
#+end_src

#+RESULTS:
: 
: Attaching package: ‘deSolve’
: 
: The following object is masked from ‘package:lsoda’:
: 
:     ode
: [1] -9.226927e-08  9.227054e-08

- The execution times are similar for ~lsoda::ode()~ and ~deSolve::ode()~:

#+begin_src R :session *R* :results output :exports both :eval no
  library(microbenchmark)
  microbenchmark(lsoda::ode(y,times,func))
  microbenchmark(deSolve::ode(y,times,func))
#+end_src

#+RESULTS:
: Unit: milliseconds
:                        expr      min       lq     mean   median       uq      max neval
:  lsoda::ode(y, times, func) 3.071473 3.164576 3.284264 3.217925 3.279549 4.419929   100
: Unit: milliseconds
:                          expr      min       lq     mean   median       uq      max neval
:  deSolve::ode(y, times, func) 2.930345 2.965816 3.106746 2.986001 3.075975 4.216918   100

*** Example: historical C interface with ~(void*) data~

- The translation from Fortran to C, and then to C++, uses the following C interface:

#+begin_src R :session *R* :results output :exports both :eval yes
  sourceCpp(code="
  // [[Rcpp::depends(lsoda)]]
  #include \"lsoda.h\"
  void fex(double t, double*  y, double* ydot, void* data) {
    ydot[0] = 1.0E4 * y[1] * y[2] - .04E0 * y[0];
    ydot[2] = 3.0E7 * y[1] * y[1];
    ydot[1] = -1.0 * (ydot[0] + ydot[2]);
    ydot[3] = y[0]+y[1]+y[2];
  }
  // [[Rcpp::export]]
  Rcpp::NumericMatrix test_lsoda_1(std::vector y, std::vector times) {
    return LSODA::ode(y, times, fex, 4, (void*) nullptr, 1e-12, 1e-12);
  }")
  test_lsoda_1(y,times) |> print(digits=12)
#+end_src

#+RESULTS:
#+begin_example
       time                y1                y2              y3 res1
 [1,] 0e+00 1.00000000000e+00 0.00000000000e+00 0.0000000000000    1
 [2,] 4e-01 9.85172113866e-01 3.38639537906e-05 0.0147940221806    1
 [3,] 4e+00 9.05518678607e-01 2.24047568782e-05 0.0944589166360    1
 [4,] 4e+01 7.15827068759e-01 9.18553476609e-06 0.2841637457064    1
 [5,] 4e+02 4.50518668519e-01 3.22290144228e-06 0.5494781085800    1
 [6,] 4e+03 1.83202257815e-01 8.94237125505e-07 0.8167968479480    1
 [7,] 4e+04 3.89833771094e-02 1.62176831695e-07 0.9610164607137    1
 [8,] 4e+05 4.93827453625e-03 1.98499409412e-08 0.9950617056138    1
 [9,] 4e+06 5.16809611427e-04 2.06829453100e-09 0.9994831883203    1
[10,] 4e+07 5.20307254747e-05 2.08133601326e-10 0.9999479690664    1
[11,] 4e+08 5.20770650260e-06 2.08309331905e-11 0.9999947922727    1
[12,] 4e+09 5.20830345924e-07 2.08332245574e-12 0.9999994791676    1
#+end_example

- The second signature attempts to simplify the use of the function. See the following examples. 

*** Example: Using a functor class

#+begin_src R :session *R* :results output :exports both :eval yes
  sourceCpp(code="
  // [[Rcpp::depends(lsoda)]]
  #include \"lsoda.h\"
  class Functor {
  public:
  Functor() {}
    std::vector operator()(double t, std::vector y) {
      std::vector ydot(4);
      ydot[0] = 1.0E4 * y[1] * y[2] - .04E0 * y[0];
      ydot[2] = 3.0E7 * y[1] * y[1];
      ydot[1] = -1.0 * (ydot[0] + ydot[2]);
      ydot[3] = y[0]+y[1]+y[2];
      return ydot;
    }
  };
  // [[Rcpp::export]]
  Rcpp::NumericMatrix test_lsoda_3(std::vector y, std::vector times) {
    Functor functor;
    return LSODA::ode(y, times, functor, 1.0e-10, 1.0e-10);
  }")
  test_lsoda_3(c(1,0,0),times) |> print(digits=12)
#+end_src

#+RESULTS:
#+begin_example
       time                y1                y2              y3 res1
 [1,] 0e+00 1.00000000000e+00 0.00000000000e+00 0.0000000000000    1
 [2,] 4e-01 9.85172113672e-01 3.38639537739e-05 0.0147940223743    1
 [3,] 4e+00 9.05518679079e-01 2.24047552536e-05 0.0944589161661    1
 [4,] 4e+01 7.15827070247e-01 9.18553499650e-06 0.2841637442181    1
 [5,] 4e+02 4.50518669367e-01 3.22290137707e-06 0.5494781077317    1
 [6,] 4e+03 1.83202258840e-01 8.94237129398e-07 0.8167968469226    1
 [7,] 4e+04 3.89833775436e-02 1.62176833676e-07 0.9610164602796    1
 [8,] 4e+05 4.93827482381e-03 1.98499421037e-08 0.9950617053263    1
 [9,] 4e+06 5.16809933531e-04 2.06829582151e-09 0.9994831879982    1
[10,] 4e+07 5.20309786314e-05 2.08134614063e-10 0.9999479688132    1
[11,] 4e+08 5.20788235135e-06 2.08316365935e-11 0.9999947920968    1
[12,] 4e+09 5.20971334176e-07 2.08388640948e-12 0.9999994790266    1
#+end_example


*** Example: Using a lambda function

- We can also use a functor that is a lambda function:

#+begin_src R :session *R* :results output :exports both :eval yes
  sourceCpp(code="
  // [[Rcpp::depends(lsoda)]]
  #include \"lsoda.h\"
  auto lambda = [](double t, std::vector y) {
      std::vector ydot(4);
      ydot[0] = 1E4 * y[1] * y[2] - .04E0 * y[0];
      ydot[2] = 3.0E7 * y[1] * y[1];
      ydot[1] = -1.0 * (ydot[0] + ydot[2]);
      ydot[3] = y[0]+y[1]+y[2];
      return ydot;
    };
  // [[Rcpp::export]]
  Rcpp::NumericMatrix test_lsoda_4(std::vector y,
                                   std::vector times,
                                   double rtol = 1-6, double atol = 1e-6) {
    return LSODA::ode(y, times, lambda, rtol, atol);
  }")
  test_lsoda_4(c(1,0,0),times,rtol=1e-10,atol=1e-10)
#+end_src

#+RESULTS:
#+begin_example
       time           y1           y2         y3 res1
 [1,] 0e+00 1.000000e+00 0.000000e+00 0.00000000    1
 [2,] 4e-01 9.851721e-01 3.386395e-05 0.01479402    1
 [3,] 4e+00 9.055187e-01 2.240476e-05 0.09445892    1
 [4,] 4e+01 7.158271e-01 9.185535e-06 0.28416374    1
 [5,] 4e+02 4.505187e-01 3.222901e-06 0.54947811    1
 [6,] 4e+03 1.832023e-01 8.942371e-07 0.81679685    1
 [7,] 4e+04 3.898338e-02 1.621768e-07 0.96101646    1
 [8,] 4e+05 4.938275e-03 1.984994e-08 0.99506171    1
 [9,] 4e+06 5.168099e-04 2.068296e-09 0.99948319    1
[10,] 4e+07 5.203098e-05 2.081346e-10 0.99994797    1
[11,] 4e+08 5.207882e-06 2.083164e-11 0.99999479    1
[12,] 4e+09 5.209713e-07 2.083886e-12 0.99999948    1
#+end_example

*** Example: Using ~RcppArmadillo~

#+begin_src R :session *R* :results output :exports both :eval yes
  sourceCpp(code="
  // [[Rcpp::depends(lsoda)]]
  // [[Rcpp::depends(RcppArmadillo)]]
  #include \"RcppArmadillo.h\"
  #include \"lsoda.h\"
  auto lambda = [](double t, arma::vec y) {
      arma::vec ydot(4);
      ydot[0] = 1E4 * y[1] * y[2] - .04E0 * y[0];
      ydot[2] = 3.0E7 * y[1] * y[1];
      ydot[1] = -1.0 * (ydot[0] + ydot[2]);
      ydot[3] = arma::sum(y);
      return ydot;
    };
  // [[Rcpp::export]]
  Rcpp::NumericMatrix test_lsoda_5(arma::vec y,
                                   arma::vec times,
                                   double rtol = 1-6, double atol = 1e-6) {
    return LSODA::ode(y, times, lambda, rtol, atol);
  }")
  test_lsoda_5(c(1,0,0),times,rtol=1e-12,atol=1e-12) |> print(digits=12)
#+end_src

#+RESULTS:
#+begin_example
       time                y1                y2              y3 res1
 [1,] 0e+00 1.00000000000e+00 0.00000000000e+00 0.0000000000000    1
 [2,] 4e-01 9.85172113866e-01 3.38639537906e-05 0.0147940221806    1
 [3,] 4e+00 9.05518678607e-01 2.24047568782e-05 0.0944589166360    1
 [4,] 4e+01 7.15827068759e-01 9.18553476609e-06 0.2841637457064    1
 [5,] 4e+02 4.50518668519e-01 3.22290144228e-06 0.5494781085800    1
 [6,] 4e+03 1.83202257815e-01 8.94237125505e-07 0.8167968479480    1
 [7,] 4e+04 3.89833771094e-02 1.62176831695e-07 0.9610164607137    1
 [8,] 4e+05 4.93827453625e-03 1.98499409412e-08 0.9950617056138    1
 [9,] 4e+06 5.16809611427e-04 2.06829453100e-09 0.9994831883203    1
[10,] 4e+07 5.20307254747e-05 2.08133601326e-10 0.9999479690664    1
[11,] 4e+08 5.20770650260e-06 2.08309331905e-11 0.9999947922727    1
[12,] 4e+09 5.20830345924e-07 2.08332245574e-12 0.9999994791676    1
#+end_example


*** Example: Using ~RcppEigen~

#+begin_src R :session *R* :results output :exports both :eval yes
  sourceCpp(code="
  // [[Rcpp::depends(lsoda)]]
  // [[Rcpp::depends(RcppEigen)]]
  #include \"RcppEigen.h\"
  #include \"lsoda.h\"
  auto lambda = [](double t, Eigen::VectorXd y) {
      Eigen::VectorXd ydot(4);
      ydot[0] = 1E4 * y[1] * y[2] - .04E0 * y[0];
      ydot[2] = 3.0E7 * y[1] * y[1];
      ydot[1] = -1.0 * (ydot[0] + ydot[2]);
      ydot[3] = y.sum();
      return ydot;
    };
  // [[Rcpp::export]]
  Rcpp::NumericMatrix test_lsoda_6(Eigen::VectorXd y,
                                   Eigen::VectorXd times,
                                   double rtol = 1-6, double atol = 1e-6) {
    return LSODA::ode(y, times, lambda, rtol, atol);
  }")
  test_lsoda_6(c(1,0,0),times,rtol=1e-12,atol=1e-12) |> print(digits=12)
#+end_src

#+RESULTS:
#+begin_example
Registered S3 methods overwritten by 'RcppEigen':
  method               from         
  predict.fastLm       RcppArmadillo
  print.fastLm         RcppArmadillo
  summary.fastLm       RcppArmadillo
  print.summary.fastLm RcppArmadillo
       time                y1                y2              y3 res1
 [1,] 0e+00 1.00000000000e+00 0.00000000000e+00 0.0000000000000    1
 [2,] 4e-01 9.85172113866e-01 3.38639537906e-05 0.0147940221806    1
 [3,] 4e+00 9.05518678607e-01 2.24047568782e-05 0.0944589166360    1
 [4,] 4e+01 7.15827068759e-01 9.18553476609e-06 0.2841637457064    1
 [5,] 4e+02 4.50518668519e-01 3.22290144228e-06 0.5494781085800    1
 [6,] 4e+03 1.83202257815e-01 8.94237125505e-07 0.8167968479480    1
 [7,] 4e+04 3.89833771094e-02 1.62176831695e-07 0.9610164607137    1
 [8,] 4e+05 4.93827453625e-03 1.98499409412e-08 0.9950617056138    1
 [9,] 4e+06 5.16809611427e-04 2.06829453100e-09 0.9994831883203    1
[10,] 4e+07 5.20307254747e-05 2.08133601326e-10 0.9999479690664    1
[11,] 4e+08 5.20770650260e-06 2.08309331905e-11 0.9999947922727    1
[12,] 4e+09 5.20830345924e-07 2.08332245574e-12 0.9999994791676    1
#+end_example

Owner

GitHub Events

Committers

Last synced: 8 months ago

Issues and Pull Requests

Last synced: 7 months ago

Packages