代码之家 › 专栏 › 技术社区 › zengc

如何使用“sugar”方式对Rcpp::numerimatrix执行逻辑操作?

rcpp r

zengc · 技术社区 · 6 年前

我必须明智地将矩阵条目与数字进行比较,因此我尝试定义一个Cxx函数,如

src <- '
LogicalMatrix f(NumericMatrix fmat, double t){
  LogicalMatrix result = fmat >= t;
  return result;
}
'
cppFunction(src)

但也有一些例外。原因是什么?那么我怎样才能做到整洁呢?

2 回复 | 直到 6 年前

Josh O'Brien 6 年前

@duckmayr的回答非常准确,并且显示了一个重要的细节:我们还可以隐藏实现细节 函数后面 因为毕竟这就是全部 Rcpp糖 不管怎样,等人都为我们做了。

但是,如果我们首先将矩阵转换为向量,对该向量进行操作,然后恢复矩阵,那么我们可以依赖于@增超所期望的Sugar操作。这是因为在内部,矩阵只是一个具有附加维度的向量(二阶;数组泛化为两个以上)。

但事实证明。。。这个版本(稍微)比循环更昂贵(而且比处理列稍微便宜)。有关完整的详细信息,请参见下面的功能 f3() 可以是:

// [[Rcpp::export]]
LogicalMatrix f3(NumericMatrix fmat, double t) {
    IntegerVector dims = fmat.attr("dim");
    NumericVector v(fmat);
    LogicalVector lv = v >= t;
    return LogicalMatrix(dims[0], dims[1], lv.begin()); 
}

但非显而易见的元素 f2() 保持最快速度:

R> microbenchmark(f(mat, 1.0), f2(mat, 1.0), f3(mat, 1.0), times = 5e4)
Unit: nanoseconds
       expr min  lq    mean median     uq     max neval
  f(mat, 1) 873 992 1322.10   1042 1118.0 1493236 50000
 f2(mat, 1) 823 925 1195.49    975 1049.5 2068214 50000
 f3(mat, 1) 864 977 1288.68   1031 1114.0 1909361 50000
R>

道德的 :简单循环解决方案复制临时对象最少,速度最快。总的来说,三者之间的速度差几乎无关紧要。

对于较大的矩阵 不复制临时文件 变得更重要:

R> mat <- matrix(sqrt(1:1000), 1000)

R> microbenchmark(f(mat, 1.0), f2(mat, 1.0), f3(mat, 1.0), times = 1e3)
Unit: microseconds
       expr   min    lq    mean median     uq    max neval
  f(mat, 1) 3.720 3.895 3.99972 3.9555 4.0425 16.758  1000
 f2(mat, 1) 1.999 2.122 2.23261 2.1760 2.2545 17.325  1000
 f3(mat, 1) 3.921 4.156 4.31034 4.2220 4.3270 19.982  1000
R>

下面是完整代码。

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
LogicalMatrix f(NumericMatrix fmat, double t){
    int n = fmat.nrow(), m = fmat.ncol();
    LogicalMatrix result(n, m);
    for ( int j = 0; j < m; ++j ) {
        result(_, j) = fmat(_, j) >= t;
    }
    return result;
}

// [[Rcpp::export]]
LogicalMatrix f2(NumericMatrix fmat, double t){
    int n = fmat.nrow(), m = fmat.ncol();
    LogicalMatrix result(n, m);
    for ( int i = 0; i < n; ++i ) {
        for ( int j = 0; j < m; ++j ) {
            result(i, j) = fmat(i, j) >= t;
        }
    }
    return result;
}

// [[Rcpp::export]]
LogicalMatrix f3(NumericMatrix fmat, double t) {
    int dims[2] = { fmat.nrow(), fmat.ncol() };
    NumericVector v(fmat);
    LogicalVector lv = v >= t;
    return LogicalMatrix(dims[0], dims[1], lv.begin()); 
}

/*** R
mat <- matrix(c(1,2,3,4), 2, 2)
library(microbenchmark)
microbenchmark(f(mat, 1.0), f2(mat, 1.0), f3(mat, 1.0), times = 1e5)

mat <- matrix(sqrt(1:1000), 1000)
microbenchmark(f(mat, 1.0), f2(mat, 1.0), f3(mat, 1.0), times = 1e3)
*/

编辑: 我们可以再删除一行相对于 f3() 但它在运行时上没有什么区别:

// [[Rcpp::export]]
LogicalMatrix f4(NumericMatrix fmat, double t) {
    IntegerVector dims = fmat.attr("dim");
    LogicalVector lv = NumericVector(fmat) >= t;
    return LogicalMatrix(dims[0], dims[1], lv.begin()); 
}

duckmayr 6 年前

我假设“整洁的方式”是指避免循环,而使用Rcpp中提供的语法糖。因为sugar为向量提供了一个具有一个值的比较器,而不是为矩阵提供了一个值(参见 here 和 here ),我认为目前最“整洁的方式”是(仅)在列(或行)上循环,即不必在列和行上循环:

// [[Rcpp::export]]
LogicalMatrix f(NumericMatrix fmat, double t){
    int n = fmat.nrow(), m = fmat.ncol();
    LogicalMatrix result(n, m);
    for ( int j = 0; j < m; ++j ) {
        result(_, j) = fmat(_, j) >= t;
    }
    return result;
}

> f(fmat, 1.0)
      [,1]  [,2]
[1,]  TRUE FALSE
[2,] FALSE  TRUE

> f(fmat, -1.0)
     [,1] [,2]
[1,] TRUE TRUE
[2,] TRUE TRUE

> f(fmat, 2.0)
      [,1]  [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE

然而,我建议避免额外的循环并不能真正为您带来可读性方面的任何好处(事实上,可能会损害代码某些读者的可读性);考虑循环行和列的函数:

// [[Rcpp::export]]
LogicalMatrix f2(NumericMatrix fmat, double t){
    int n = fmat.nrow(), m = fmat.ncol();
    LogicalMatrix result(n, m);
    for ( int i = 0; i < n; ++i ) {
        for ( int j = 0; j < m; ++j ) {
            result(i, j) = fmat(i, j) >= t;
        }
    }
    return result;
}

我真的看不出这怎么会很难输入,它似乎本质上是性能等效的(平均执行时间略低,但中位数略高——请参见下面的基准测试),至少对于一些读者来说,我敢打赌,它会更清楚地显示出你在做什么。

也就是说,如果跳过一个循环对你有帮助,我认为这是目前你能做的最好的事情。

library(microbenchmark)

> microbenchmark(loop = f(fmat, 1.0), nonloop = f2(fmat, 1.0), times = 1e4)
Unit: microseconds
    expr   min    lq     mean median    uq      max neval cld
    loop 6.564 7.402  9.77116  7.612 8.031 9173.952 10000   a
 nonloop 6.425 7.123 10.01659  7.333 7.682 4377.448 10000   a

> microbenchmark(nonloop = f2(fmat, 1.0), loop = f(fmat, 1.0), times = 1e4)
Unit: microseconds
    expr   min    lq      mean median    uq      max neval cld
 nonloop 6.356 7.124 10.179950  7.333 7.544 4822.066 10000   a
    loop 6.775 7.404  9.588326  7.613 7.892 4278.971 10000   a