如何解决关于在RcppArmadillo中使用.shed_row / .shed_col
我现在正在尝试将R代码转换为Rcpp代码。
R代码是
hh <- function(A,B,j){
aa <- A[j,-j] %*% B[j,-j] ## A and B are n*m matrixes
return(aa)
}
> set.seed(123)
> A <- matrix(runif(30),5,6)
> B <- matrix(rnorm(30),6)
> j <- 2
> hh(A,j)
> [,1]
> [1,] 0.9702774
我的Rcpp代码是
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
double hh(arma::mat A,arma::mat B,arma::uword j){
arma::mat Bj = B.shed_col(j); /* error occurs */
arma::mat Ak = A.row(j);
double aa = Ak.shed_col(j) * arma::trans(Bj.row(j)); /* error occurs */
return aa;
}
该错误应该与.shed_row / .shed_col的使用有关。我用Google.shed_row搜索,但是还没有一个想法来解决我在这里遇到的问题。你有什么主意吗?预先谢谢你!
进一步更新:
现在我们考虑在函数的for循环中使用.shed_row/.shed_col
。
具体来说,我的Rcpp代码如下
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat ggcpp(arma::mat A,arma::mat B){
/*we assume here A and B are n*n square matrices.*/
int Ac = A.n_cols;
int Bc = B.n_cols;
arma::mat C(Ac,Bc);
for(arma::uword i = 0; i < Ac; i++){
A.shed_col(i);
for(arma::uword j = 0; j < Bc; j++){
B.shed_col(j);
C(i,j) = arma::as_scalar(A.row(i) * B.row(j).t());
}
}
return C;
}
等效的R代码如下
gg <- function(A,B){
ac <- ncol(A)
bc <- ncol(B)
C <- matrix(NA,ac,bc)
for(i in 1:ac){
for(j in 1:bc){
C[i,j] <- A[i,-i] %*% B[j,-j]
}
}
return(C)
}
我的R代码有效。已经过测试。但是,我在使用Rcpp代码时遇到了麻烦。
我尝试了很多方法,主要遇到两个错误:
- ggcpp(a1,a2)中的错误:as_scalar():不兼容的尺寸
- ggcpp(a1,a2)中的错误:Mat :: shed_col():索引超出范围
在这里,a1
和a2
是两个随机生成的6 * 6矩阵。
你有什么主意吗?感激!!!
解决方法
docs在这里很好。行/列被删除。因此,在更改返回类型和矩阵乘法的类型并记住索引在c ++中从零开始之后,就足够了
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat hhcpp(arma::mat A,arma::mat B,arma::uword j){
A.shed_col(j); B.shed_col(j);
arma::mat aa = A.row(j) * B.row(j).t();
return aa;
}
/***R
set.seed(123)
A <- matrix(runif(30),5,6)
B <- matrix(rnorm(30),6)
j <- 2
hh <- function(A,B,j){
aa <- A[j,-j] %*% B[j,-j] ## A and B are n*m matrixes
return(aa)
}
hh(A,j)
hhcpp(A,j-1)
*/
重新更新;另一种策略是选择要保留的矩阵列,而不是删除列。这是通过使用arma::regspace
沿着列创建一个序列,然后保留不使用find
删除的列来完成的。
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat ggcpp(arma::mat A,arma::mat B){
int Ac = A.n_cols;
int Bc = B.n_cols;
arma::mat C(Ac,Bc);
arma::vec aRange = arma::regspace<arma::vec>(0,Ac-1);
arma::vec bRange = arma::regspace<arma::vec>(0,Bc-1);
for(int i = 0; i < Ac; i++){
arma::mat Atemp = A.submat(i* arma::ones<arma::uvec>(1),find(aRange != i)) ;
for(int j = 0; j < Bc; j++){
arma::mat Btemp = B.submat(j* arma::ones<arma::uvec>(1),find(bRange != j)) ;
C(i,j) = arma::as_scalar(Atemp * Btemp.t());
}
}
return C;
}
/***R
# Square matrices only!!
set.seed(123)
A <- matrix(runif(30),5)
B <- matrix(rnorm(30),5)
j <- 2
gg <- function(A,B){
ac <- ncol(A)
bc <- ncol(B)
C <- matrix(NA,ac,bc)
for(i in 1:ac){
for(j in 1:bc){
C[i,j] <- A[i,-i] %*% B[j,-j]
}
}
return(C)
}
all.equal(gg(A,B),ggcpp(A,B))
*/
,
回答进一步更新:
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::NumericMatrix row_erase (Rcpp::NumericMatrix& x,Rcpp::IntegerVector& rowID) {
rowID = rowID.sort();
Rcpp::NumericMatrix x2(Rcpp::Dimension(x.nrow()- rowID.size(),x.ncol()));
int iter = 0;
int del = 1; // to count deleted elements
for (int i = 0; i < x.nrow(); i++) {
if (i != rowID[del - 1]) {
x2.row(iter) = x.row(i);
iter++;
} else {
del++;
}
}
return x2;
}
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::NumericMatrix col_erase (Rcpp::NumericMatrix& x,Rcpp::IntegerVector& colID) {
colID = colID.sort();
Rcpp::NumericMatrix x2(Rcpp::Dimension(x.nrow(),x.ncol()- colID.size()));
int iter = 0;
int del = 1;
for (int i = 0; i < x.ncol(); i++) {
if (i != colID[del - 1]) {
x2.column(iter) = x.column(i);
iter++;
} else {
del++;
}
}
return x2;
}
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat ggcpp(arma::mat A,arma::mat B){
Rcpp::NumericMatrix AA = Rcpp::wrap(A);
Rcpp::NumericMatrix BB = Rcpp::wrap(B);
Rcpp::NumericMatrix AAi;
Rcpp::NumericMatrix BBj;
Rcpp::IntegerVector AiV;
Rcpp::IntegerVector BjV;
arma::mat Ai;
arma::mat Bj;
unsigned int Ac = A.n_cols;
unsigned int Bc = B.n_cols;
arma::mat C(Ac,Bc);
for(arma::uword i = 0; i < Ac; i++){
AiV = {i};
AAi = col_erase(AA,AiV);
Ai = Rcpp::as<arma::mat>(AAi);
for(arma::uword j = 0; j < Bc; j++){
BjV = {j};
BBj = col_erase(BB,BjV);
Bj = Rcpp::as<arma::mat>(BBj);
C(i,j) = arma::as_scalar(Ai.row(i) * Bj.row(j).t());
}
}
return C;
}
注意:row_erase
和col_erase
是从here借来的。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。