如何使用Rcpp在C++中用另一个矩阵的值替换矩阵元素?

5

我有一个使用Rcpp编写的小型C++函数,用于将一个矩阵中的元素替换为另一个矩阵中的值。对于单个单元格或者一列来说,它都能正常工作,如下所示:

cppFunction('NumericMatrix changeC(NumericMatrix one, NumericMatrix two) {
NumericMatrix a = one;
NumericMatrix b = two;
b(_,1) = a(_,1);
return b;
}')

changeC(g,f)

如果原来的矩阵f如下:
      [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    6    6    6    6    6    6
[2,]    6    6    6    6    6    6
[3,]    6    6    6    6    6    6
[4,]    6    6    6    6    6    6
[5,]    6    6    6    6    6    6
[6,]    6    6    6    6    6    6

并且 g 矩阵如下所示:

        [,1] [,2] [,3] [,4] [,5] [,6]
 [1,]    5    5    5    5    5    5
 [2,]    5    5    5    5    5    5
 [3,]    5    5    5    5    5    5
 [4,]    5    5    5    5    5    5
 [5,]    5    5    5    5    5    5
 [6,]    5    5    5    5    5    5

当我运行 changeC(g,f) 函数时,我得到了预期的结果:
      [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    6    5    6    6    6    6
[2,]    6    5    6    6    6    6
[3,]    6    5    6    6    6    6
[4,]    6    5    6    6    6    6
[5,]    6    5    6    6    6    6
[6,]    6    5    6    6    6    6

但是我真正想做的是用另一个矩阵的子集替换另一个矩阵的子集(例如,将一个矩阵(3*3)的第1到3行、第1到3列替换为另一个矩阵的第3到6行、第3到6列(也是3*3)。我已经尝试过:

cppFunction('NumericMatrix changeC(NumericMatrix one, NumericMatrix two) {
NumericMatrix a = one;
NumericMatrix b = two;
b( Range(0,2), Range(0,2)) = a( Range(3,5), Range(3,5));
return b;
}')

但是这段代码无法编译。尽管如此:

cppFunction('NumericMatrix changeC(NumericMatrix one, NumericMatrix two) {
NumericMatrix a = one;
NumericMatrix b = two;
b = a( Range(3,5), Range(3,5));
return b;
}')

我的代码是这样的:f[1:3, 1:3] = g[4:6, 4:6],但是它运行得比较慢,尤其是对于非常大的矩阵(因此我使用了Rcpp)。请问我做错了什么?在R中,我会这样做。

谢谢您的帮助。

编辑1

经过一番尝试,我已经成功让我的矩阵向东和向西移动了(我认为向北和向南也可能类似,可能需要一个两步的方法来实现东北、西北等方向的移动?)。

func <- 'NumericMatrix eastC(NumericMatrix a) {
int acoln=a.ncol();
NumericMatrix out(a.nrow(),a.ncol()) ;
for (int j = 0;j < acoln;j++) {
if (j > 0) {
out(_,j) = a(_,j-1);
 } else {
  out(_,j) = a(_,0);
 }
 }
 return out ;
 }'
 cppFunction(func)

欢迎进行任何改进。我希望将第一列保留为零而不是0列。有什么想法吗?


1
如果您使用几个嵌套的for循环会发生什么? - IRTFM
1
在从R转向C++时,您必须执行许多R可以本地执行的步骤。例如- R是矢量化的,而C++则不是。正如42所指出的那样,通过两个循环(行,列)来完成这个过程是一种选择。另一个选项是使用RcppEigen中的块操作。详见:块操作教程 - alexwhitworth
1
此外,您应该传递引用并在矩阵上使用“const”(在您的第一个示例中为“one”),以避免其被更改。 - alexwhitworth
艾力克斯,42岁,我一定想探索你所说的内容,但不确定该如何做。你能发表一个答案吗? - AntonyDW
我今天下午可能有时间这样做。 - alexwhitworth
显示剩余2条评论
1个回答

7

我认为Rcpp的subMatrix不支持那种分配方式。

建议您尝试使用RcppArmadillo和Armadillo子矩阵视图

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace arma;

// [[Rcpp::export]]
mat example( mat m1, mat m2) {
  m1.submat( 0,0, 2,2) = m2.submat( 3,3, 5,5 );
  return m1;
}

/*** R
m1 <- matrix(1,6,6)
m2 <- matrix(-1,6,6)
example(m1, m2)
*/

> m1 <- matrix(1,6,6)

> m2 <- matrix(-1,6,6)

> example(m1, m2)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   -1   -1   -1    1    1    1
[2,]   -1   -1   -1    1    1    1
[3,]   -1   -1   -1    1    1    1
[4,]    1    1    1    1    1    1
[5,]    1    1    1    1    1    1
[6,]    1    1    1    1    1    1

我想尝试RcppArmadillo - 但是在哪里编写此代码?在R中?还是在另一个编辑器中编写并进行编译?当我尝试在R中编写时,我得到的都是错误信息 - 可以看出我没有使用Armadillo的经验。我的C ++代码是在R中编写和编译的。 - AntonyDW
1
创建一个 test.cpp 文件,并使用 Rcpp 包中的 'sourceCpp("test.cpp")'。此外,阅读 Rcpp 文档。 :) - Thell
我白费了劲。谢谢你的指点,我试了一下,现在明白了。非常感谢,我整天都在苦苦挣扎如何读取这种类型的代码。 - AntonyDW
使用Armadillo肯定有潜力,现在我可以使用它了 - 我不认为它会比我的上面的代码快多少,但肯定更加可控。 - AntonyDW
2
我应该更快,正如Alex在评论中所说; 使用引用传递(并不要忘记const)。 - Thell
最终,我使用了嵌套的for循环来将矩阵向左或向右移动一步(模拟E或W)- 这非常快 - 以及Armadillo subset.matrix来将矩阵向上或向下移动一行(模拟S或N)- 比for循环方法慢但比它快一些。混合使用两种方法来实现SE,NW等方向。真的很高兴。谢谢。 - AntonyDW

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