I wrote an R package that compiles user-written c++ code with the help of Rcpp and RcppEigen. Every time the user wants to create a new model, he/she must write three files:
svol_leverage_bswc.h),svol_leverage_bswc.cpp ),svol_leverage_bswc_export.cpp).Currently, all three files are saved in the src/ directory, and compilation is triggered by calling devtools::load_all(".") After compilation, the functions callable from R are exported to NAMESPACE. Everything compiles just fine in this package framework.
Unfortunately, though, I think I designed this package incorrectly. I expect the users to change these files frequently, and different users to create different models than other users. So, I would like to save these files outside of the package directory (easy), and to compile user-written code with cxxfunplus::cxxfunctionplus() (I stole this idea from RStan's rstan::stan_model()). I also wrote my inline plugin based off of their inline plugin.
So how does my package compilation map to cxxfunction compilation?
My current attempt is below, but I'm pretty sure this is way off. The error I'm getting is
Error in compileCode(f, code, language = language, verbose = verbose) :
file25a041c179ba.cpp:12:10: fatal error: pf/include/pf/bootstrap_filter_with_covariates.h: No such file or directory 12 | #include "pf/include/pf/bootstrap_filter_with_covariates.h" // the boostrap particle filter | ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~compilation terminated.make: *** [/usr/lib/R/etc/Makeconf:177: file25a041c179ba.o] Error 1
Stan has a line like this in their plugin that finds all their header files (e.g. system.file("include", "stan", "math", "prim", package = "StanHeaders", mustWork = TRUE)). However, my package's header files don't appear to be installed anywhere in that same directory. Why don't my headers get installed?
src <- r"(// replace every instance of <TODO> with your own code!
#include "svol_leverage_bswc.h"
Svol_leverageBSWC::Svol_leverageBSWC(double phi, double mu, double sigma, double rho)
{
m_params << phi, mu, sigma, rho;
}
auto Svol_leverageBSWC::q1Samp(const ovec &y1, const cvec& z1) -> svec
{
// phi, mu, sigma, rho
svec x1samp;
x1samp(0) = m_example_stdNormSampler.sample() * m_params(2) / std::sqrt(1.0 - m_params(0) * m_params(0));
return x1samp;
}
auto Svol_leverageBSWC::fSamp(const svec &xtm1, const cvec& zt) -> svec
{
// phi, mu, sigma, rho
svec xt;
double mean = m_params(1) + m_params(0) * (xtm1(0) - m_params(1)) +
zt(0) * m_params(3) * m_params(2) * std::exp(-.5 * xtm1(0));
xt(0) = mean + m_example_stdNormSampler.sample() * m_params(2) * std::sqrt(1.0 - m_params(3) * m_params(3));
return xt;
}
double Svol_leverageBSWC::logGEv(const ovec &yt, const svec &xt, const cvec& zt)
{
return rveval::evalUnivNorm<double>(
yt(0),
0.0,
std::exp(.5 * xt(0)),
true);
}
double Svol_leverageBSWC::logMuEv(const svec &x1, const cvec& z1)
{
// parameter order phi, mu, sigma, rho
return rveval::evalUnivNorm<double>(
x1(0),
0.0,
m_params(2) / std::sqrt(1.0 - m_params(0) * m_params(0)),
true);
}
double Svol_leverageBSWC::logQ1Ev(const svec &x1, const ovec &y1, const cvec& z1)
{
// parameter order phi, mu, sigma, rho
return rveval::evalUnivNorm<double>(
x1(0),
0.0,
m_params(2) / std::sqrt(1.0 - m_params(0) * m_params(0)),
true);
}
)"
header <- r"(// replace every instance of <TODO> with your own code!
#ifndef SVOL_LEVERAGE_BSWC_H
#define SVOL_LEVERAGE_BSWC_H
#include "pf/include/pf/bootstrap_filter_with_covariates.h" // the boostrap particle filter
#include "pf/include/pf/rv_samp.h" // for sampling random numbers
#include "pf/include/pf/rv_eval.h" // for evaluating densities and pmfs
#include "pf/include/pf/resamplers.h" // for resampling classes
using namespace pf;
using namespace pf::filters;
using namespace pf::resamplers;
#define nparts_SVOL_LEVERAGE_BSWC 500 // number of particles
#define dimstate_SVOL_LEVERAGE_BSWC 1 // dimension of state vectors
#define dimobs_SVOL_LEVERAGE_BSWC 1 // dimension of observation vectors
#define dimcov_SVOL_LEVERAGE_BSWC 1 // dimension of covariate vectors
#define dimparam_SVOL_LEVERAGE_BSWC 4 // dimension of parameters
// helper type aliases
using resampT = mn_resamp_fast1<nparts_SVOL_LEVERAGE_BSWC,dimstate_SVOL_LEVERAGE_BSWC,double>;
using svec = Eigen::Matrix <double, dimstate_SVOL_LEVERAGE_BSWC,1>;
using ovec = Eigen::Matrix <double, dimobs_SVOL_LEVERAGE_BSWC,1>;
using cvec = Eigen::Matrix <double, dimcov_SVOL_LEVERAGE_BSWC,1>;
using param_vec = Eigen::Matrix <double, dimparam_SVOL_LEVERAGE_BSWC,1>;
using DynMat = Eigen::Matrix <double, Eigen::Dynamic, Eigen::Dynamic>;
using func = std::function <const DynMat(const svec&, const cvec&)>;
using BasePF = BSFilterWC <nparts_SVOL_LEVERAGE_BSWC, dimstate_SVOL_LEVERAGE_BSWC, dimobs_SVOL_LEVERAGE_BSWC, dimcov_SVOL_LEVERAGE_BSWC, resampT, double, true>;
/**
* @brief a particle filter class template for a TODO model
*
*/
class Svol_leverageBSWC : public BasePF
{
private:
// example parameter
param_vec m_params;
// use this for sampling
rvsamp::UnivNormSampler<double> m_example_stdNormSampler;
// required by algorithm and required to define your model
double logQ1Ev (const svec &x1, const ovec &y1, const cvec &z1 );
double logMuEv (const svec &x1, const cvec &z1 );
double logGEv (const ovec &yt, const svec &xt, const cvec& zt );
auto q1Samp (const ovec &y1, const cvec& z1 ) -> svec;
auto fSamp (const svec &xtm1, const cvec &zt ) -> svec;
public:
// constructor
Svol_leverageBSWC(double phi, double mu, double sigma, double rho);
};
#endif // SVOL_LEVERAGE_BSWC_H)"
# //' @useDynLib pfr, .registration = TRUE
# //' @import RcppEigen
# //' @importFrom Rcpp evalCpp
# //' @export
# // [[Rcpp::export]]
# double svol_leverage_bswc_approx_LL(const Rcpp::NumericVector& obsTS, const Rcpp::NumericVector& params){
#
body <- r"(
// instantiate model with arg params
Svol_leverageBSWC mod(params[0],params[1],params[2],params[3]);
// calculate likelihood by iterating through time series
double ans(0.0);
ovec yt, ytm1;
for (int i = 1; i < obsTS.length(); i++){
yt << obsTS[i]; ytm1 << obsTS[i-1];
mod.filter(yt, ytm1);
ans += mod.getLogCondLike();
}
return ans;
)"
#inline::registerPlugin("pfr", pfrplugin)
thisMod <- cxxfunplus::cxxfunctionplus(signature(), body = body, plugin = "pfr", includes = header, verbose = T)
Following the Rcpp modules vignette and jamming it all in includes= works:
inc <- r"(
#ifndef SVOL_LEVERAGE_BSWC_H
#define SVOL_LEVERAGE_BSWC_H
#include "pf/bootstrap_filter_with_covariates.h" // the boostrap particle filter
#include "pf/rv_samp.h" // for sampling random numbers
#include "pf/rv_eval.h" // for evaluating densities and pmfs
#include "pf/resamplers.h" // for resampling classes
using namespace pf;
using namespace pf::filters;
using namespace pf::resamplers;
#define nparts_SVOL_LEVERAGE_BSWC 500 // number of particles
#define dimstate_SVOL_LEVERAGE_BSWC 1 // dimension of state vectors
#define dimobs_SVOL_LEVERAGE_BSWC 1 // dimension of observation vectors
#define dimcov_SVOL_LEVERAGE_BSWC 1 // dimension of covariate vectors
#define dimparam_SVOL_LEVERAGE_BSWC 4 // dimension of parameters
// helper type aliases
using resampT = mn_resamp_fast1<nparts_SVOL_LEVERAGE_BSWC,dimstate_SVOL_LEVERAGE_BSWC,double>;
using svec = Eigen::Matrix <double, dimstate_SVOL_LEVERAGE_BSWC,1>;
using ovec = Eigen::Matrix <double, dimobs_SVOL_LEVERAGE_BSWC,1>;
using cvec = Eigen::Matrix <double, dimcov_SVOL_LEVERAGE_BSWC,1>;
using param_vec = Eigen::Matrix <double, dimparam_SVOL_LEVERAGE_BSWC,1>;
using DynMat = Eigen::Matrix <double, Eigen::Dynamic, Eigen::Dynamic>;
using func = std::function <const DynMat(const svec&, const cvec&)>;
using BasePF = BSFilterWC <nparts_SVOL_LEVERAGE_BSWC, dimstate_SVOL_LEVERAGE_BSWC, dimobs_SVOL_LEVERAGE_BSWC, dimcov_SVOL_LEVERAGE_BSWC, resampT, double, true>;
/**
* @brief a particle filter class template for a TODO model
*
*/
class Svol_leverageBSWC : public BasePF
{
private:
// example parameter
param_vec m_params;
// use this for sampling
rvsamp::UnivNormSampler<double> m_example_stdNormSampler;
// required by algorithm and required to define your model
double logQ1Ev (const svec &x1, const ovec &y1, const cvec &z1 );
double logMuEv (const svec &x1, const cvec &z1 );
double logGEv (const ovec &yt, const svec &xt, const cvec& zt );
auto q1Samp (const ovec &y1, const cvec& z1 ) -> svec;
auto fSamp (const svec &xtm1, const cvec &zt ) -> svec;
public:
// constructor
Svol_leverageBSWC(double phi, double mu, double sigma, double rho);
};
#endif // SVOL_LEVERAGE_BSWC_H
double svol_leverage_bswc_approx_LL(const Rcpp::NumericVector& obsTS, const Rcpp::NumericVector& params){
// instantiate model with arg params
Svol_leverageBSWC mod(params[0],params[1],params[2],params[3]);
// calculate likelihood by iterating through time series
double ans(0.0);
ovec yt, ytm1;
for (int i = 1; i < obsTS.length(); i++){
yt << obsTS[i]; ytm1 << obsTS[i-1];
mod.filter(yt, ytm1);
ans += mod.getLogCondLike();
}
return ans;
}
Svol_leverageBSWC::Svol_leverageBSWC(double phi, double mu, double sigma, double rho)
{
m_params << phi, mu, sigma, rho;
}
auto Svol_leverageBSWC::q1Samp(const ovec &y1, const cvec& z1) -> svec
{
// phi, mu, sigma, rho
svec x1samp;
x1samp(0) = m_example_stdNormSampler.sample() * m_params(2) / std::sqrt(1.0 - m_params(0) * m_params(0));
return x1samp;
}
auto Svol_leverageBSWC::fSamp(const svec &xtm1, const cvec& zt) -> svec
{
// phi, mu, sigma, rho
svec xt;
double mean = m_params(1) + m_params(0) * (xtm1(0) - m_params(1)) +
zt(0) * m_params(3) * m_params(2) * std::exp(-.5 * xtm1(0));
xt(0) = mean + m_example_stdNormSampler.sample() * m_params(2) * std::sqrt(1.0 - m_params(3) * m_params(3));
return xt;
}
double Svol_leverageBSWC::logGEv(const ovec &yt, const svec &xt, const cvec& zt)
{
return rveval::evalUnivNorm<double>(
yt(0),
0.0,
std::exp(.5 * xt(0)),
true);
}
double Svol_leverageBSWC::logMuEv(const svec &x1, const cvec& z1)
{
// parameter order phi, mu, sigma, rho
return rveval::evalUnivNorm<double>(
x1(0),
0.0,
m_params(2) / std::sqrt(1.0 - m_params(0) * m_params(0)),
true);
}
double Svol_leverageBSWC::logQ1Ev(const svec &x1, const ovec &y1, const cvec& z1)
{
// parameter order phi, mu, sigma, rho
return rveval::evalUnivNorm<double>(
x1(0),
0.0,
m_params(2) / std::sqrt(1.0 - m_params(0) * m_params(0)),
true);
}
)"
inline::registerPlugin("pfr", pfrplugin)
yay <- cxxfunplus::cxxfunctionplus(signature(), plugin = "pfr", includes = inc)