如何使用“糖”方式对Rcpp :: NumericMatrix进行逻辑运算?

问题描述 投票:2回答:2

我必须将一个矩阵入口与数字进行比较,所以我尝试定义一个Cxx函数,如

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

但是有些例外被抛弃了。是什么原因?那我怎么能以一种整洁的方式做到这一点呢?

r rcpp
2个回答
3
投票

我假设“整洁的方式”你的意思是避免使用Rcpp中提供的句法糖循环。因为糖为比较器提供了一个矢量值而不是矩阵值(参见herehere),我认为你现在可以做的最“整洁的方法”就是(仅)循环(或者)在列(或行)上,即没有必须循环遍历列和行:

// [[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

3
投票

@duckmayr的回答非常重要,并且显示了一个重要的细节:我们可能隐藏一个函数背后的实现细节,因为毕竟Rcpp Sugar等人都为我们做了。

但是如果我们首先将矩阵转换为向量,对该向量进行操作然后恢复矩阵,我们可以依赖@zengchao所需的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()); 
}
© www.soinside.com 2019 - 2024. All rights reserved.