Reputation: 7989
I coded a function qSelectMbycol
in Rcpp
that returns the k
th largest element of each column in O(n) time. This function works OK. If I try to do the same but work over rows instead of columns (function qSelectMbyrow
) it returns the error "error: Mat::init(): requested size is not compatible with column vector layout"
. Anybody any thoughts what I am doing wrong? I saved this file as "qselect.cpp
" :
// [[Rcpp::depends(RcppArmadillo)]]
#define RCPP_ARMADILLO_RETURN_COLVEC_AS_VECTOR
#include <RcppArmadillo.h>
using namespace arma;
// [[Rcpp::export]]
arma::vec qSelectMbycol(arma::mat& M, const int k) {
// ARGUMENTS
// M: matrix for which we want to find the k-th largest elements of each column
// k: k-th statistic to look up
arma::mat Y(M.memptr(), M.n_rows, M.n_cols);
// we apply over columns
int c = M.n_cols;
arma::vec out(c);
int i;
for (i = 0; i < c; i++) {
arma::vec y = Y.col(i);
std::nth_element(y.begin(), y.begin() + k - 1, y.end());
out(i) = y(k-1); // the k-th largest value of each column
}
return out;
}
// [[Rcpp::export]]
arma::vec qSelectMbyrow(arma::mat& M, const int k) {
// ARGUMENTS
// M: matrix for which we want to find the k-th largest elements of each row
// k: k-th statistic to look up
arma::mat Y(M.memptr(), M.n_rows, M.n_cols);
// we apply over rows
int r = M.n_rows;
arma::vec out(r);
int i;
for (i = 0; i < r; i++) {
arma::vec y = Y.row(i); // this line throws the error "error: Mat::init(): requested size is not compatible with column vector layout"
std::nth_element(y.begin(), y.begin() + k - 1, y.end());
out(i) = y(k-1); // should give k-th largest value of each row
}
return out;
}
Example:
n=500
p=100
set.seed(1)
M=matrix(rnorm(n, mean = 100, sd = 1),n,1)
library(Rcpp)
library(RcppArmadillo)
Rcpp::sourceCpp('qselect.cpp')
qSelectMbycol(M,5) # works OK
qSelectMbyrow(M,5) # throws error "error: Mat::init(): requested size is not compatible with column vector layout"
I also tried inserting
typedef std::vector<double> stdvec;
and replacing the line setting vector y
by
arma::vec y = arma::conv_to<stdvec>::from(Y.row(i));
in my qSelectMbyrow
function and although the function then runs, it runs slowly compared to applying over columns, and also crashes my R session if I run it 100 times.
Upvotes: 2
Views: 2861
Reputation: 16920
The issue is that an arma::vec
is actually an arma::colvec
(see the docs). So, we can solve this issue by changing
arma::vec y = Y.row(i);
(which is incompatible because it thinks you want a matrix with one column but you're trying to give it a matrix with one row) to
arma::rowvec y = Y.row(i);
Upvotes: 4