采用Rcpparamdillo架起R与Armadillo的桥梁

引言

我们都知道R语言的强大,而R最大的瓶颈是运算能力,所以众多R包都采用了c++来优化底层的代码,armadillo是C++的一个强大的线性代数运算库,本博文讲的就是如何把这几个对象链接起来,让我们的R代码可以使用高效的运算库来提高效率。

导入/引用 Armadillo

常规创建RcppArmadillo骨骼包的方式

install
install.packages('RcppArmadillo')
下面这步可以用于观看使用说明

??RcppArmadillo 

Create Rcpp Skeleton

list <- YOUR_VAR
RcppArmadillo.package.skeleton(name = "YOUR_PACKAGE_NAME", list, environment = .GlobalEnv, path = ".", force = FALSE, code_files = character(), example_code = TRUE)

在已有Rcpp项目导入RcppArmadillo

经过以上两步我们可以得到一个具有RcppArmadillo架构的R包,通过观察发现Demo,该过程引入了

#include <RcppArmadillo.h>

那么我们可以从这方面入手
添加 RcppArmadillo.h
在src中的某个全局头文件加入

#include <RcppArmadillo.h>

注意,此时我们应该去掉

#include <Rcpp.h>

因为这两个文件是冲突的,这在后面的编译过程中会出错提醒,如

error: "The file 'Rcpp.h' should not be included. Please correct to include only 'RcppArmadillo.h'."

我们通过观察Demo包发现,在Descirption文件中,LinkingTo字段后面有RcppArmadillo,其实在编译过程中也会报错提醒, 如

The following packages are referenced using Rcpp::depends attributes however are not listed in the Depends, Imports or LinkingTo fields of the package DESCRIPTION file: RcppArmadillo

添加 LinkingTo 描述
所以我们在包文件的根目录下找到DESCRIPTION文件,修改如下
找到

LinkingTo: Rcpp

加上RcppArmadillo

LinkingTo: Rcpp, RcppArmadillo

到此,导入完成,编译成功。

不再重复添加
若重复include了Armadillo.h会导致冲突致使奔溃,注意此问题可以通过编译,但在运行时会报错。

联系R function参数与C++ function参数

方式
在主函数的参数中用使用 arma:: + Armadillo中的类来声明。
比如vec, colvec, rowvec, matrix类等,这里以matrix类为例:

function(arma::mat& X, ...){...}

注意,&符合的位置,这里&表示绑定而不是X的地址,所以是接在声明类后面的。
最后加上与Rcpp一样的注释。

  // [[Rcpp::export]]

官方demo

// [[Rcpp::export]]
List fastLm(const arma::mat& X, const arma::colvec& y) {
    int n = X.n_rows, k = X.n_cols;

    arma::colvec coef = arma::solve(X, y);    // fit model y ~ X
    arma::colvec res  = y - X*coef;           // residuals

    // std.errors of coefficients
    double s2 = std::inner_product(res.begin(), res.end(), res.begin(), 0.0)/(n - k);

    arma::colvec std_err = arma::sqrt(s2 * arma::diagvec(arma::pinv(arma::trans(X)*X)));

    return List::create(Named("coefficients") = coef,
                        Named("stderr")       = std_err,
                        Named("df.residual")  = n - k);
}

使用Armadillo类

引入'arma'命名空间
在头文件加上

using namespace arma;

这样一来就可以直接声明mat, rowvec等Armadillo类了而不用在前部加上'arma::'。

rowvec a(10);
mat b(10, 10);

但是,与上节的主函数参数声明无关,参数仍然需要声明为’arma::‘ + ...
使用openblas作为Armadillo的Blas库
在R project下的src文件夹里新建编译参数文件Makevars。
它的作用与~/.R/Makevars一致,区别也很明显,一个是全局的,即对所有的R文件生效,另一种是对该项目生效,为了让RcppArmadillo使用openblas的Blas库,我们在该文件作如下定义:

OPENBLAS_LIB_PATH = /usr/local/opt/openblas/lib

BLAS_LIBS := -L$(OPENBLAS_LIB_PATH) -lopenblas

PKG_LIBS := $(BLAS_LIBS)

PKG_CXXFLAGS = -O3 -mavx -std=c++11

注意带了参数运算符 $ 需要用 := 来赋值

总结

到此为止,我们就可以导入RcppArmadillo接受R传过来的参数并使用Armadillo类,而对于Armadillo的用法,待我完成GSoC后来填坑。

2017-06-21 18:55 101
Comments
Write a Comment