Rcpp矩阵:逐行循环,逐列处理

21

这是我第一次尝试Rcpp,但这个非常简单的问题让我困扰。我想使用嵌套的for循环逐列操作矩阵中的每个值。我希望最终的脚本看起来像这样:

src <- '
    Rcpp::NumericMatrix Am(A);
    int nrows = Am.nrow();
    int ncolumns = Am.ncol();
    for (int i = 0; i < ncolumns; i++){
        for (int j = 1; j < nrows; j++){
            Am[j,i] = Am[j,i] + Am[j-1,i];
        }
    }
    return Am;
'
fun <- cxxfunction(signature(A = "numeric"), body = src, plugin="Rcpp")
fun(matrix(1,4,4))

期望的输出应该是这样的:

     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    2    2    2    2
[3,]    3    3    3    3
[4,]    4    4    4    4
明显问题出在这行代码,我不知道如何引用矩阵中的单个元素。
Am[j,i] = Am[j,i] + Am[j-1,i];

如果这是个愚蠢的新手问题,我很抱歉。任何提示都将不胜感激!


3
我以前说过,现在再说一遍:rcpp-devel是提出这些问题的更好地方。 - Dirk Eddelbuettel
3
虽然我理解 rcpp-devel 列表可能会更容易接触到有经验的 rcpp 用户,但在我看来,stackoverflow 更加易于访问。 - jbaums
@jbaums:当然,但是在所有关键的rcpp-devel贡献者中,只有一个人在这里看到问题。问题的关注度更低... - Dirk Eddelbuettel
1个回答

31

在单个[ ]表达式中,您不能使用多个索引。这是一种C语言的限制,我所知道的没有任何一个C++矩阵类系统或库可以克服它。因此使用( )代替。

修复这个错误和你没有实际传递srccxxfunction()的错误,我们得到了以下代码:

R> src <- '
+     Rcpp::NumericMatrix Am(A);
+     int nrows = Am.nrow();
+     int ncolumns = Am.ncol();
+     for (int i = 0; i < ncolumns; i++) {
+         for (int j = 1; j < nrows; j++) {
+             Am(j,i) = Am(j,i) + Am(j-1,i);
+         }
+     }
+     return Am;
+ '
R> fun <- cxxfunction(signature(A = "numeric"), body = src, plugin="Rcpp")
R> fun(matrix(1,4,4))
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    2    2    2    2
[3,]    3    3    3    3
[4,]    4    4    4    4
R> 

最后,注意Rcpp sugar有用于一次处理整个行或列的示例,请参阅邮件列表存档和vignette。

编辑: 仅为明确起见,这里使用了只有一个循环和Rcpp sugar的按列索引方法:

R> src <- '
+     Rcpp::NumericMatrix Am(A);
+     int nrows = Am.nrow();
+     for (int j = 1; j < nrows; j++) {
+         Am(j,_) = Am(j,_) + Am(j-1,_);
+     }
+     return Am;
+ '
R> fun <- cxxfunction(signature(A = "numeric"), body = src, plugin="Rcpp")
R> fun(matrix(1,4,4))
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    2    2    2    2
[3,]    3    3    3    3
[4,]    4    4    4    4
R> 

3
太好了!谢谢你的回复,Dirk。如果将来有任何问题,我会把它们直接发到Rcpp-devel列表中。再次感谢! - Vincent
6
你可以使用Am.row(j)代替Am(j,_) - Artem Klevtsov

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接