列表示3d矩阵(立方体)Rcpp

Question

我有一个程序,我需要X(nRow, nCol, nSlice)在Rcpp中重复计算每个立方体切片的列平均值,结果意味着形成一个矩阵M(nCol, nSlice).以下代码产生错误:

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp; 
using namespace arma;

// [[Rcpp::export]]

mat cubeMeans(arma::cube X){
   int nSlice = X.n_slices;
   int nCol = X.n_cols;
   int nRow = X.n_rows;
   arma::vec Vtmp(nCol);
   arma::mat Mtmp(nRow, nCol);
   arma::mat Means(nCol, nSlice);
   for (int i = 0; i < nSlice; i++){
      Mtmp = X.slice(i);
      for(int j = 0; j < nCol; j++){
         Vtmp(j) = sum(Mtmp.col(j))/nRow; 
      }
      Means.col(i) = Vtmp;
   }
  return(wrap(Means));
}

'/Rcpp/internal/Exporter.h:31:31:错误:没有匹配函数来调用'arma :: Cube :: Cube(SEXPREC*&)'

我无法弄明白.当函数的输入是一个矩阵(并返回一个向量)时,我没有得到错误.但是,我将上述功能作为主程序的一部分包括在内

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

mat cubeMeans(arma::cube X){
  int nSlice = X.n_slices;
  ...
  return(Means);
}

// [[Rcpp::export]]

main part of program

该程序编译成功,但速度很慢(几乎与使用的R版本一样慢colMeans).有没有更好的方法来计算多维数据集上的列方式,为什么我得到编译错误？

我很感激任何帮助.

问候,

Answer 1

尝试使用arma::cube作为Rcpp函数参数时,我也收到此错误.根据编译器错误,我相信这是因为Rcpp::wrap<arma::cube>当前没有定义(需要处理你传递给函数的R对象).†在线阅读了几个相关示例后,看起来典型的解决方法是将R array作为a 读入NumericVector,并且由于它保留了其dims属性,因此使用这些来设置arma::cube维度.尽管需要额外的一两步来解释缺失的wrap专业化 †,但我放在一起的Armadillo版本似乎比我的R解决方案要快得多:

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

// [[Rcpp::export]]
arma::mat cube_means(Rcpp::NumericVector vx) {

  Rcpp::IntegerVector x_dims = vx.attr("dim");
  arma::cube x(vx.begin(), x_dims[0], x_dims[1], x_dims[2], false);

  arma::mat result(x.n_cols, x.n_slices);
  for (unsigned int i = 0; i < x.n_slices; i++) {
    result.col(i) = arma::conv_to<arma::colvec>::from(arma::mean(x.slice(i)));  
  }

  return result;
}

/*** R

rcube_means <- function(x) t(apply(x, 2, colMeans))

xl <- array(1:10e4, c(100, 100 ,10))
all.equal(rcube_means(xl), cube_means(xl))
#[1] TRUE

R> microbenchmark::microbenchmark(
    "R Cube Means" = rcube_means(xl),
    "Arma Cube Means" = cube_means(xl),
    times = 200L)
Unit: microseconds
            expr      min       lq      mean   median       uq       max neval
    R Cube Means 6856.691 8204.334 9843.7455 8886.408 9859.385 97857.999   200
 Arma Cube Means  325.499  380.540  643.7565  416.863  459.800  3068.367   200

*/

我正在利用s 的arma::mean函数重载arma::mat将默认计算列均值的事实(arma::mean(x.slice(i), 1)将给出该切片的行方式).

编辑: †第二个想法,我不确定这是否与此Rcpp::wrap有关 - 但问题似乎与Rcpp的Exporter.h第31行缺少Exporter<>专业化有关arma::cube:

template <typename T>
class Exporter{
public:
  Exporter( SEXP x ) : t(x){}
  inline T get(){ return t ; }

private:
  T t ;
} ;

无论如何,NumericVector我使用的/设置尺寸方法似乎是现在的功能解决方案.

根据您在问题中描述的输出维度,我假设您希望得到的矩阵的每一列都是相应数组切片的列均值的向量(第1列=第1列的列均值等),即

R> x <- array(1:27, c(3, 3, 3))
R> rcube_means(x)
     [,1] [,2] [,3]
[1,]    2   11   20
[2,]    5   14   23
[3,]    8   17   26
R> cube_means(x)
     [,1] [,2] [,3]
[1,]    2   11   20
[2,]    5   14   23
[3,]    8   17   26

但如果需要,你可以改变这一点.