1

I am learning to use RcppArmadillo(Rcpp) to make the code run faster. In practice, I usually encounter that I want to call some functions from given packages. The below is a little example. I want to calculate the thresholds of lasso (or scad, mcp). In R we can use the thresh_est function in library(HSDiC), which reads

library(HSDiC) # acquire the thresh_est() function
# usage: thresh_est(z, lambda, tau, a = 3, penalty = c("MCP", "SCAD", "lasso"))
# help: ?thresh_est

z = seq(-5,5,length=500)
thresh = thresh_est(z,lambda=1,tau=1,penalty="lasso")
# thresh = thresh_est(z,lambda=1,tau=1,a=3,penalty="MCP")
# thresh = thresh_est(z,lambda=1,tau=1,a=3.7,penalty="SCAD")
plot(z,thresh)

Then I try to realize the above via RcppArmadillo(Rcpp). According to the answers of teuder as well as coatless and coatless, I write the code (which is saved as thresh_arma.cpp) below:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp; //with this, Rcpp::Function can be Function for short, etc.
using namespace arma;  // with this, arma::vec can be vec for short, etc.
// using namespace std; //std::string

// [[Rcpp::export]]
vec thresh_arma(vec z, double lambda, double a=3, double tau=1, std::string penalty) {
    // Obtaining namespace of the package
    Environment pkg = namespace_env("HSDiC");
    // Environment HSDiC("package:HSDiC");

    // Picking up the function from the package
    Function f = pkg["thresh_est"];
    // Function f = HSDiC["thresh_est"];

    vec thresh;
    //  In R: thresh_est(z, lambda, tau, a = 3, penalty = c("MCP", "SCAD", "lasso"))
    if (penalty == "lasso") {
        thresh = f(_["z"]=z,_["lambda"]=lambda,_["a"]=a,_["tau"]=tau,_["penalty"]="lasso");
    } else if (penalty == "scad") {
        thresh = f(_["z"]=z,_["lambda"]=lambda,_["a"]=a,_["tau"]=tau,_["penalty"]="SCAD");
    } else if (penalty == "mcp") {
        thresh = f(_["z"]=z,_["lambda"]=lambda,_["a"]=a,_["tau"]=tau,_["penalty"]="MCP");
}
return thresh;
}

Then I do the compilation as

library(Rcpp)
library(RcppArmadillo)
sourceCpp("thresh_arma.cpp")

# z = seq(-5,5,length=500)
# thresh = thresh_arma(z,lambda=1,tau=1,penalty="lasso")
# # thresh = thresh_arma(z,lambda=1,a=3,tau=1,penalty="mcp")
# # thresh = thresh_arma(z,lambda=1,a=3.7,tau=1,penalty="scad")
# plot(z,thresh)

However, the compilation fails and I have no idea about the reasons.

John Stone
  • 635
  • 4
  • 13
  • 1
    Your function `f()` is an R function which you access from C++. It will still run at the speed of an R function. Also, for us to help you we need to see the error messages. For which it would help to have a _smaller_ example. It is unclear what is failing, and there is a little too much going on here. Maybe you can edit it down? – Dirk Eddelbuettel Jan 28 '20 at 03:17
  • Tks@DirkEddelbuettel! I'm afraid that my little example `thresh_arma.cpp` is small enough one though, just one package (`library(HSDiC)`) and one function (`thresh_est()`). I saied that **the compilation fails** which means when I run `library(Rcpp);library(RcppArmadillo); sourceCpp("thresh_arma.cpp")`, there are lots of error messages. Maybe I will try another solution way (perhaps manually write the thresholds). Tks again! – John Stone Jan 28 '20 at 09:30
  • You can tackle the error messages one at a time: 1) With C++ you cannot have an argument w/o default value after one w/ a default value. 2) Use `Environment::namespace_env`. 3) The result of calling an R function is a SEXP. I don't know why you want to assign it to an `arma::vec`, but you can do so using `Rcpp::as(f(...))`. Alternativly you could use a `Rcpp::NumericVector`. – Ralf Stubner Jan 28 '20 at 15:11
  • Tks@RalfStubner! I will try what you suggested. – John Stone Jan 29 '20 at 02:07

0 Answers0