I am now trying to convert R code to Rcpp code.
The R code is
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),5,6)
> j <- 2
> hh(A,B,j)
> [,1]
> [1,] 0.9702774
My Rcpp code is
#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;
}
The error should be about the use of .shed_row/.shed_col. I have Googled .shed_row, however, did not have an idea yet to address the issue I encountered here. Do you have any idea? Thank you in advance!
Further Update:
Now we consider using .shed_row/.shed_col
in for-loop in the function.
Specifically, my Rcpp code is the following
#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;
}
The equivalent R code is the following
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)
}
My R code works. It has been tested. However, I am having troubles with the Rcpp code.
I tried many ways and mainly encountered two errors :
Here, a1
and a2
are two randomly generated 6*6 matrices.
Do you have any idea? Appreciated!!!
Answer to Further Update:
#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;
}
Note: row_erase
and col_erase
are borrowed from here.