Rcpp: Difference between revisions

From 太極
Jump to navigation Jump to search
Line 595: Line 595:
* [https://github.com/cole-trapnell-lab/monocle3 monocle3] - an analysis toolkit for single-cell RNA-Seq experiments.
* [https://github.com/cole-trapnell-lab/monocle3 monocle3] - an analysis toolkit for single-cell RNA-Seq experiments.
* [https://github.com/immunogenomics/harmony Harmony] - Fast, sensitive and accurate integration of single-cell data
* [https://github.com/immunogenomics/harmony Harmony] - Fast, sensitive and accurate integration of single-cell data
* [https://github.com/MarioniLab/DropletUtils/blob/master/src/downsample_run.cpp DropletUtils]
* [https://github.com/MarioniLab/DropletUtils/blob/master/src/downsample_run.cpp DropletUtils] and [https://github.com/LTLA/scuttle/tree/master/src scuttle]


= Rcpp sugar =
= Rcpp sugar =

Revision as of 15:42, 17 June 2021

Misc

sudo apt-get install libblas-dev liblapack-dev
sudo apt-get install gfortran

Speed Comparison

  • A comparison of high-performance computing techniques in R. It compares Rcpp to an R looping operator (like mapply), a parallelized version of a looping operator (like mcmapply), explicit parallelization, via the parallel package or the ParallelR suite.
  • In the following example, C++ avoids the overhead of creating an intermediate object (eg vector of the same length as the original vector). The c++ uses an intermediate scalar. So C++ wins R over memory management in this case. 
    # http://blog.mckuhn.de/2016/03/avoiding-unnecessary-memory-allocations.html 
library(Rcpp)

`%count<%` <- cppFunction('
size_t count_less(NumericVector x, NumericVector y) {
  const size_t nx = x.size();
  const size_t ny = y.size();
  if (nx > 1 & ny > 1) stop("Only one parameter can be a vector!");
  size_t count = 0;
  if (nx == 1) {
    double c = x[0];
    for (int i = 0; i < ny; i++) count += c < y[i];
  } else {
    double c = y[0];
    for (int i = 0; i < nx; i++) count += x[i] < c;
  }
  return count;
}
')

set.seed(42)

N <- 10^7
v <- runif(N, 0, 10000)

# Testing on my ODroid xu4 running ubuntu 15.10
system.time(sum(v < 5000))
#   user  system elapsed
#  1.135   0.305   1.453
system.time(v %count<% 5000)
#   user  system elapsed
#  0.535   0.000   0.540
  • Why R for data science – and not Python?
    library(Rcpp)
 
bmi_R <- function(weight, height) {
  weight / (height * height)
}
bmi_R(80, 1.85) # body mass index of person with 80 kg and 185 cm
## [1] 23.37473
 
cppFunction("
  float bmi_cpp(float weight, float height) {
    return weight / (height * height);
  }
")
bmi_cpp(80, 1.85) # same with cpp function
## [1] 23.37473
  • Boost the speed of R calls from Rcpp
  • Is the pain worth it?: Can Rcpp speed up Passing Bablok Regression? Several headers are included: vector, cmath, limits, algorithm, numeric.

Use Rcpp in RStudio

RStudio makes it easy to use Rcpp package.

Open RStudio, click New File -> C++ File. It will create a C++ template on the RStudio editor 

#include <Rcpp.h>
using namespace Rcpp;

// Below is a simple example of exporting a C++ function to R. You can
// source this function into an R session using the Rcpp::sourceCpp 
// function (or via the Source button on the editor toolbar)

// For more on using Rcpp click the Help button on the editor toolbar

// [[Rcpp::export]]
int timesTwo(int x) {
   return x * 2;
}

Now in R console, type 

library(Rcpp)
sourceCpp("~/Downloads/timesTwo.cpp")
timesTwo(9)
# [1] 18

See more examples on http://adv-r.had.co.nz/Rcpp.html and Calculating a fuzzy kmeans membership matrix

If we wan to test Boost library, we can try it in RStudio. Consider the following example in stackoverflow.com. 

// [[Rcpp::depends(BH)]]
#include <Rcpp.h>
#include <boost/foreach.hpp>
#include <boost/math/special_functions/gamma.hpp>

#define foreach BOOST_FOREACH

using namespace boost::math;

//[[Rcpp::export]]
Rcpp::NumericVector boost_gamma( Rcpp::NumericVector x ) {
  foreach( double& elem, x ) {
    elem = boost::math::tgamma(elem);
  };

  return x;
}

Then the R console 

boost_gamma(0:10 + 1)
#  [1]       1       1       2       6      24     120     720    5040   40320
# [10]  362880 3628800

identical( boost_gamma(0:10 + 1), factorial(0:10) )
# [1] TRUE

From the useR 2019 tutorial

useR 2019 Rcpp tutorial

Number of packages depends on Rcpp

db <- tools::CRAN_package_db() # added in R 3.4.0
## rows: number of pkgs, cols: different attributes
nTot <- nrow(db)
## all direct Rcpp reverse depends, ie packages using Rcpp
nRcpp <- length(tools::dependsOnPkgs("Rcpp", recursive=FALSE, installed=db))
nCompiled <- table(db[, "NeedsCompilation"])[["yes"]]
propRcpp <- nRcpp / nCompiled * 100
data.frame(tot=nTot, totRcpp = nRcpp, totCompiled = nCompiled,
           RcppPctOfCompiled = propRcpp)
#     tot totRcpp totCompiled RcppPctOfCompiled
# 1 14594    1710        3646          46.90071

evalCpp()

evalCpp("2 + 2")

set.seed(42)
evalCpp("Rcpp::rnorm(2)")

library(Rcpp)
evalCpp("std::numeric_limits<double>::max()")

cppFunction()

https://www.rdocumentation.org/packages/Rcpp/versions/1.0.1/topics/cppFunction

Note that the c/cpp function can be used as a regular R function without the need of .Call() function. 

cppFunction("
int exampleCpp11() {
auto x = 10; // guesses type
return x;
}", plugins=c("cpp11"))
// exampleCpp11()

Rcpp::cppFunction("Rcpp::NumericVector f(int n) { 
  srand(1234); 
  Rcpp::NumericVector res(n); 
  for(int i=0; i<n; i++) res[i] = rand(); 
  return(res); 
}")
//  f(); different result on different platform

cppFunction("int f(int a, int b) { return(a + b); }")
// f(21, 21)

// fibonacci
cppFunction('int g(int n) {
if (n < 2) return(n);
return(g(n-1) + g(n-2)); }')
// sapply(0:10, g)

sourceCpp(): the actual workhorse behind evalCpp() and cppFunction()

https://www.rdocumentation.org/packages/Rcpp/versions/1.0.1/topics/sourceCpp

Create the following in a file called testrcpp.cpp and in R run Rcpp::sourceCpp("testrcpp.cpp"). Any exported functions in .cpp will be available in an R session. 

#include <Rcpp.h>
using namespace Rcpp;
// This is a simple example of exporting a C++ function to R. You can
// source this function into an R session using the Rcpp::sourceCpp()
// [[Rcpp::export]]
NumericVector timesTwo(NumericVector x) {
return x * 2;
}
// [[Rcpp::export]]
int g(int n) {
if (n < 2) return(n);
return(g(n-1) + g(n-2)); 
}
// timesTwo(c(2, 5.3))
// sapply(0:10, g)

g++ command line options

  • header location via -I,
  • library location via -L,
  • library via -llibraryname
g++ -I/usr/include -c qnorm_rmath.cpp
g++ -o qnorm_rmath qnorm_rmath.o -L/usr/lib -lRmath

C++ and R types

http://dirk.eddelbuettel.com/papers/useR2019_rcpp_tutorial.pdf#page=56

5 Types of Vectors

In R everything is a vector 

// Vector
// [[Rcpp::export]]
double getMax(NumericVector v) {
    int n = v.size(); // vectors are describing
    double m = v[0]; // initialize
    for (int i=0; i<n; i++) {
        if v[i] > m {
            Rcpp::Rcout << ”Now ” << m << std::endl;
            m = v[i];
        }
    }
    return(m);
}

// Matrix
// [[Rcpp::export]]
Rcpp::NumericVector colSums(Rcpp::NumericMatrix mat) {
  size_t cols = mat.cols();
  Rcpp::NumericVector res(cols);
  for (size_t i=0; i<cols; i++) {
    res[i] = sum(mat.column(i));
  }
  return(res);
}

STL vectors (it copies data)

cppFunction("double getMax2(std::vector<double> v) {
  int n = v.size(); // vectors are describing
  double m = v[0]; // initialize
  for (int i=0; i<n; i++) {
    if (v[i] > m) {
      m = v[i];
    }
  }
  return(m);
}")

getMax2(c(4,5,2))

Sugar functions

max()

Packages

  • Overview
    1. Rcpp.package.skeleton("mypackage"). This will create a cpp source file src/rcpp_hello_world.cpp.
    2. Run Rcpp::compileAttributes() which generates src/RcppExports.cpp and R/RcppExports.R (do not edit by hand)
    3. Edit DESCRIPTION
    4. optional: Makevars and Makevars.win
    5. NAMESPACE
    6. rcpp_hello_world.Rd
    7. mypackage-package.Rd
  • Writing a package that uses Rcpp vignette
  • Thirteen Simple Steps for Creating An R Package with an External C++ Library
  • Read mypackage/Read-and-delete-me
  • The new package can be installed by devtools::load_all() or install.package("mypackage", repos = NULL, type = "source")
  • From Advanced R. Before building the package, you’ll need to run Rcpp::compileAttributes(). This function scans the C++ files for Rcpp::export attributes and generates the code required to make the functions available in R. Re-run compileAttributes() whenever functions are added, removed, or have their signatures changed. This is done automatically by the devtools package and by Rstudio.
  • Examples

Armadillo

Machine Learning

Examples

Example 1. convolution example

First, Rcpp package should be installed (I am working on Linux system). Next we try one example shipped in Rcpp package.

PS. If R was not available in global environment (such as built by ourselves), we need to modify 'Makefile' file by replacing 'R' command with its complete path (4 places). 

cd ~/R/x86_64-pc-linux-gnu-library/3.0/Rcpp/examples/ConvolveBenchmarks/
make
R

Then type the following in an R session to see how it works. Note that we don't need to issue library(Rcpp) in R. 

dyn.load("convolve3_cpp.so")
x <- .Call("convolve3cpp", 1:3, 4:6)
x # 4 13 28 27 18

If we have our own cpp file, we need to use the following way to create dynamic loaded library file. Note that the character (grave accent) ` is not (single quote)'. If you mistakenly use ', it won't work. 

export PKG_CXXFLAGS=`Rscript -e "Rcpp:::CxxFlags()"`
export PKG_LIBS=`Rscript -e "Rcpp:::LdFlags()"`
R CMD SHLIB xxxx.cpp

Example 2. Use together with inline package

library(inline)
x <- as.numeric(1:10)
n <- as.integer(10)
code2 <- "
      integer i
      do 1 i=1, n 
    1 x(i) = x(i)**3
"
cubefn2 <- cfunction(signature(n="integer", x="numeric"), 
                     implicit = "none", 
                     dim = c("", "(*)"), 
                     code2, 
                     convention=".Fortran")
cubefn2(n, x)$x
#  [1]    1    8   27   64  125  216  343  512  729 1000
library(inline)
src <-'
 Rcpp::NumericVector xa(a);
 Rcpp::NumericVector xb(b);
 int n_xa = xa.size(), n_xb = xb.size();

 Rcpp::NumericVector xab(n_xa + n_xb - 1);
 for (int i = 0; i < n_xa; i++)
 for (int j = 0; j < n_xb; j++)
 xab[i + j] += xa[i] * xb[j];
 return xab;
'
fun <- cxxfunction(signature(a = "numeric", b = "numeric"),
 src, plugin = "Rcpp")
fun(1:3, 1:4) 
# [1]  1  4 10 16 17 12
  • Simulate VAR(1) data. See Section 1.3.3 C++ Solution of the book 'Seamless R and C++ Integration with Rcpp'.
a <- matrix(c(0.5, 0.1, 0.1, 0.5), nrow = 2)
u <- matrix(rnorm(10000), ncol=2)
code <- '
  arma::mat coeff = Rcpp::as<arma::mat> (a);
  arma::mat errors = Rcpp::as<arma::mat>(u);
  int m = errors.n_rows;
  int n = errors.n_cols;
  arma::mat simdata(m, n);
  simdata.row(0) = arma::zeros<arma::mat>(1, n);
  for (int row=1; row<m; row++) {
    simdata.row(row) = simdata.row(row-1)*trans(coeff) +
                       errors.row(row);
  }
  return Rcpp::wrap(simdata);
'
rcppSim <- inline::cxxfunction(signature(a = "numeric", u = "numeric"),
                       code, plugin = "RcppArmadillo")
rcppData <- rcppSim(a, u)

rSim <- function(coeff, errors) {
    simdata <- matrix(0, nrow(errors), ncol(errors))
    for(row in 2:nrow(errors)) {
      simdata[row, ] <- coeff %*% simdata[(row-1), ] + errors[row, ]
    }
    return(simdata)
}
rData <- rSim(a, u)

> benchmark(rcppSim(a, u), rSim(a, u),  
  columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"),
  order= "relative")
           test replications elapsed relative user.self sys.self
1 rcppSim(a, u)          100   0.026    1.000     0.022    0.004
2    rSim(a, u)          100   1.557   59.885     1.557    0.000

Example 3. Calling Rmath functions, random number generation

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector mypnorm(NumericVector x) {
  int n = x.size();
  NumericVector y(n);
  for (int i=0; i<n; i++)
    y[i] = R::pnorm(x[i], 0.0, 1.0, 1, 0);
  return y;
}
sourceCpp("code/using-rmath-rng.cpp")
x <- seq(-2, 2)
mypnorm(x)
pnorm(x)
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix rngCpp(const int N) {
  NumericMatrix X(N, 4);
  X(_, 0) = runif(N);
  X(_, 1) = rnorm(N);
  X(_, 2) = rt(N, 5);
  X(_, 3) = rbeta(N, 1, 1);
  return X;
}
set.seed(42)     # setting seed
 M1 <- rngCpp(5)
 M1
set.seed(42)	  # resetting seed
 M2 <- cbind( runif(5), rnorm(5), rt(5, 5), rbeta(5, 1, 1))
 M2
all.equal(M1, M2)

Example 4: Bootstrap

See Rcpp introduction vignette.

R version 

# Function declaration
bootstrap_r <- function(ds, B = 1000) {
  # Preallocate storage for statistics
  boot_stat <- matrix(NA, nrow = B, ncol = 2)
  # Number of observations
  n <- length(ds)
  # Perform bootstrap
  for(i in seq_len(B)) {
    # Sample initial data
    gen_data <- ds[ sample(n, n, replace=TRUE) ]
    # Calculate sample data mean and SD
    boot_stat[i,] <- c(mean(gen_data),
    sd(gen_data))
  }
  # Return bootstrap result
  return(boot_stat)
}

Rcpp version 

#include <Rcpp.h>
// Function declaration with export tag
// [[Rcpp::export]]
Rcpp::NumericMatrix
bootstrap_cpp(Rcpp::NumericVector ds, int B = 1000) {
  // Preallocate storage for statistics
  Rcpp::NumericMatrix boot_stat(B, 2);
  // Number of observations
  int n = ds.size();
  // Perform bootstrap
  for(int i = 0; i < B; i++) {
    // Sample initial data
    Rcpp::NumericVector gen_data =
    ds[ floor(Rcpp::runif(n, 0, n)) ];
    // Calculate sample mean and std dev
    boot_stat(i, 0) = mean(gen_data);
    boot_stat(i, 1) = sd(gen_data);
  }
  // Return bootstrap results
  return boot_stat;
}

Comparison 

# Set seed to generate data
set.seed(512)
# Generate data
initdata <- rnorm(1000, mean = 21, sd = 10)

# Use the same seed in R and C++
set.seed(883)
# Perform bootstrap with C++ function
result_cpp <- bootstrap_cpp(initdata)

# Set a new _different_ seed for bootstrapping
set.seed(883)
# Perform bootstrap
result_r <- bootstrap_r(initdata)

# Compare output
all.equal(result_r, result_cpp)

library(rbenchmark)
benchmark(r = bootstrap_r(initdata), cpp = bootstrap_cpp(initdata))[, 1:4]

Example 5: Calling R functions from C++

https://gallery.rcpp.org/articles/r-function-from-c++/

Example 6: Using Rcout for output synchronised with R

https://gallery.rcpp.org/articles/using-rcout/

Example 7: Passing user-supplied C++ functions

https://gallery.rcpp.org/articles/passing-cpp-function-pointers/

RcppExamples package

There is no need to install it. Download the source and use Rcpp to run the examples.

  • RNGs.cpp
> library(Rcpp)
> sourceCpp("~/Downloads/RcppExamples/src/RNGs.cpp")
> set.seed(42)
> x <- RcppRNGs(5L)
       rnorm            rt rpois
1  1.3709584   -0.05492941     2
2 -0.5646982   -0.03887139     1
3  0.3631284    1.99724535     1
4  0.6328626    0.97864627     1
5  0.4042683 -216.32581066     0
> set.seed(42)
> y <- data.frame(rnorm=rnorm(5), rt =rt(5,1), rpois=rpois(5,1))
> all.equal(x, y)
  • MatrixExample.cpp. This uses STL transform() algorithm. The cpp needs to 'include' <cmath>.
> sourceCpp("~/Downloads/RcppExamples/src/MatrixExample.cpp")
> M <- matrix((1:16)^2, 4)
> MatrixExample(M)
  • NumericVectorExample.cpp. This uses STL transform() algorithm. The cpp needs to 'include' <cmath>.
> sourceCpp("~/Downloads/RcppExamples/src/NumericVectorExample.cpp")
> NumericVectorExample(seq(1,9)^2)
  • StringVectorExample.cpp. This uses STL transform() algorithm. No extra needs to be 'included'. It transforms characters to lower case.
> sourceCpp("~/Downloads/RcppExamples/src/StringVectorExample.cpp")
> StringVectorExample(c("Tick", "Tack", "Tock"))
  • DataFrameExample.cpp
> sourceCpp("~/Downloads/RcppExamples/src/DataFrameExample.cpp")
> D <- data.frame(a=1:3,
                  b=LETTERS[1:3],
                  c=as.Date("2011-01-01")+0:2,
                  stringsAsFactors=FALSE)
> val <- DataFrameExample(D)
  • RcppListExample.cpp
> sourceCpp("~/Downloads/RcppExamples/src/ListExample.cpp")
> params <- list(method='BFGS',
               tolerance=1.0e-5,
               maxIter=100,
               startDate=as.Date('2006-7-15'))

# call the underlying  C++ function
> result <- ListExamples(params)

Some packages that import Rcpp

Rcpp sugar

The motivation of Rcpp sugar is to bring a subset of the high-level R syntax in C++. For example, consider a simple R function 

foo <- function(x, y) {
  ifelse (x < y, x*x, -(y*y) )
}

With Rcpp instead of writing 

RcppExport SEXP foo( SEXP x, SEXP y){
  Rcpp::NumericVector xx(x), yy(y) ;
  int n = xx.size() ;
  Rcpp::NumericVector res( n ) ;
  double x_ = 0.0, y_ = 0.0 ;
  for( int i=0; i<n; i++){
    x_ = xx[i] ;
    y_ = yy[i] ;
    if( x_ < y_ ){
      res[i] = x_ * x_ ;
    } else {
      res[i] = -( y_ * y_) ;
    }
  }
  return res ;
}

, we can use the ifelse function from R 

RcppExport SEXP foo( SEXP xs, SEXP ys){
  Rcpp::NumericVector x(xs) ;
  Rcpp::NumericVector y(ys) ;
  return Rcpp::wrap( ifelse( x < y, x*x, -(y*y) )) ;
}

Operators

  • Binary arithmetic operators,
NumericVector x;
NumericVector y;

NumericVector res = x + y;
NumericVector res = x * y; // element-wise multiplication

NumericVector res = x + 2.0;

NumericVector res = x * y + 7 / 2.0; // two expressions
  • Binary logical operators: <, <=, >=, ==, !=
  • Unary operators
NumericVector res = -x;
LogicalVector res = ! x;

Functions

  • Functions producing a single logical result
IntegerVector x = seq_len( 1000 );
all( x*x < 3 );

// Respect the concept of missing values (NA) in R,
// expression generated by 'any' or 'all' cannot be converted 
// directly to 'bool'. Instead one must use 'is.true',
// 'if_false' or 'is_na'
book res = is_true( any (x < y) );
  • Functions producing sugar expressions
is_na
seq_along
seq_len
pmin and pmax
ifelse
sapply
lapply
mapply
sign
diff
setdiff
union
intersect
clamp
unique
sort_unique
table
duplicated
  • Mathematical functions
NumericVector x, y;
int k;
double z;

abs( x )
exp( x )
floor( x )
ceil( x )
pow(x, z)
log(x);
log10(x);
sqrt(x);
sin(x); cos(x); tan(x); sinh(x); cosh(x); tanh(x); asin(x); acos(x); atan(x);
gamma(x); lgamma(x); digamma(x); trigamma(x);
factorial(x); choose(n, k); beta(n, k); trunc(x); round(x, k); signif(x,k)
mean(x); var(x); sd(x); sum(x); cumsum(x); min(x), max(x); range(x); 
which_min(x); which_max(x); setequal(x, y);
  • The d/q/p/r statistical functions
NumericVector y1, y2, y3;
x1 = dnorm(y1, 0, 1);
x2 = pnorm(y2, 0, 1);
x3 = qnorm(y3, 0, 1);
x4 = rnorm(n, 0, 1);

// Initialize the state of the random number generator
RcppExport SEXP getRGamma() {
  RNGScope scope;
  NumericVector x = rgamma( 10, 1, 1 );
  return x;
}

Implementation

  • The curiously recurring template pattern
  • The VectorBase class
  • Example : sapply

RcppParallel

rTRNG: parallel RNG in R

rTRNG: ADVANCED PARALLEL RNG IN R

RcppAnnoy

https://github.com/eddelbuettel/rcppannoy

Annoy library is a small and lightweight C++ template header library for very fast approximate nearest neighbors - originally developed to drive the famous Spotify music discovery algorithm.

R6

Handling R6 objects in C++

Retrieved from "https://wiki.taichimd.us/index.php?title=Rcpp&oldid=29399"