Search code examples
rrcpp

call Rcpp function in R


below is the R code where functions stl_example_model and stl_example_jacobian in a dll are referenced.

library(RcppSundials)
library(microbenchmark)

stl_example_model = getNativeSymbolInfo(name = "example_model_stl",
                                        PACKAGE = "RcppSundials")$address
stl_example_jacobian = getNativeSymbolInfo(name = "example_jacobian_stl",
                                           PACKAGE = "RcppSundials")$address

simulationstl = wrap_cvodes(times = 1:5, 
                            states_ = rep(1.0,5), 
                            parameters_ = 0.1, 
                            forcings_data_ =list(cbind(1:3600,1:3600)),
            ..., 
                            model_ = example_model, jacobian_ = example_jacobian) 

I have sourced another function (test.cpp) using the Rcpp::sourceCpp. How do I a make a reference to this function (similar to stl_example_model) using Rcpp?

Test.cpp

#define ARMA_DONT_USE_CXX11
#include <RcppArmadillo.h>
#include <Rcpp.h>
#include <array>
#include <vector>
using namespace std;
using namespace Rcpp;

extern "C" {
  array<vector<double>, 2> example_model_stl(const double& t, const vector<double>& states, 
            const vector<double>& parameters, const vector<double>& forcings) {

    vector<double> derivatives(states.size());  
    for(int i = 0; i < states.size(); i++) {
      derivatives[i] = -states[i]*parameters[0];
    }
    vector<double> observed{forcings[0]};
    array<vector<double>, 2> output{derivatives, observed};
    return output;
  }

  arma::mat example_jacobian_stl(const double& t, const vector<double>& states, 
            const vector<double>& parameters, const vector<double>& forcings) {
    arma::mat output = arma::eye(states.size(), states.size());
    output = -parameters[0]*output;
    return output;
  }

  array<vector<double>, 2> example_dae_stl(const double& t, const vector<double>& states, 
            const vector<double>& derivatives, const vector<double>& parameters, 
            const vector<double>& forcings) {

    vector<double> residues(states.size());  
    residues[0] = -0.04*states[0] + 1e4*states[1]*states[2];
    residues[1] = -residues[0] - 3e7*states[1]*states[1] - derivatives[1];
    residues[0] = residues[0] - derivatives[0];
    residues[2] = states[0] + states[1] + states[2] - 1.0;
    vector<double> observed{forcings[0]};
    array<vector<double>, 2> output{residues, observed};
    return output;
  }

};

Solution

  • It does make a lot of sense to use a package as mentioned by Dirk in the comments. However, it is not strictly necessary.

    1. The PACKAGE argument of getNativeSymbols is actually optional. All loaded DLLs are searched if it is omitted.

    2. Compiling your example produces the warning

      No Rcpp::export attributes or RCPP_MODULE declarations found in source

      As a consequence the produced DLL is not loaded into R's memory. You can change this by adding an exported dummy function.

    Example:

    // [[Rcpp::depends("RcppArmadillo")]]
    #define ARMA_DONT_USE_CXX11
    #include <RcppArmadillo.h>
    #include <array>
    #include <vector>
    using namespace std;
    using namespace Rcpp;
    
    extern "C" {
        array<vector<double>, 2> example_model_stl(const double& t, const vector<double>& states, 
                                                   const vector<double>& parameters, const vector<double>& forcings) {
    
            vector<double> derivatives(states.size());  
            for(int i = 0; i < states.size(); i++) {
                derivatives[i] = -states[i]*parameters[0];
            }
            vector<double> observed{forcings[0]};
            array<vector<double>, 2> output{derivatives, observed};
            return output;
        }
    
        arma::mat example_jacobian_stl(const double& t, const vector<double>& states, 
                                       const vector<double>& parameters, const vector<double>& forcings) {
            arma::mat output = arma::eye(states.size(), states.size());
            output = -parameters[0]*output;
            return output;
        }
    
        array<vector<double>, 2> example_dae_stl(const double& t, const vector<double>& states, 
                                                 const vector<double>& derivatives, const vector<double>& parameters, 
                                                 const vector<double>& forcings) {
    
            vector<double> residues(states.size());  
            residues[0] = -0.04*states[0] + 1e4*states[1]*states[2];
            residues[1] = -residues[0] - 3e7*states[1]*states[1] - derivatives[1];
            residues[0] = residues[0] - derivatives[0];
            residues[2] = states[0] + states[1] + states[2] - 1.0;
            vector<double> observed{forcings[0]};
            array<vector<double>, 2> output{residues, observed};
            return output;
        }
    
    };
    
    // [[Rcpp::export]]
    void dummy() {}
    
    /***R
    getNativeSymbolInfo("example_model_stl")$address
    */
    

    Output:

    > Rcpp::sourceCpp('58881676.cpp')
    
    > getNativeSymbolInfo("example_model_stl")$address
    <pointer: 0x7fd438ed5f80>
    attr(,"class")
    [1] "NativeSymbol"
    

    This way it should be possible to use such functions together with RcppSundials.