列均值3d矩阵(立方体)Rcpp [英] Column means 3d matrix (cube) Rcpp

问题描述

我有一个程序，其中我需要重复计算Rcpp中立方体X(nRow, nCol, nSlice)的每个切片的列均值，所得的均值形成矩阵M(nCol, nSlice).以下代码产生了错误:

I have a program in which I need to calculate repeatedly the column means of each slice of a cube X(nRow, nCol, nSlice) in Rcpp, with the resulting means forming a matrix M(nCol, nSlice). The following code produced an error:

#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 *&)'

'/Rcpp/internal/Exporter.h:31:31: error: no matching function for call to 'arma::Cube::Cube(SEXPREC*&)'

我不太清楚.当函数的输入是矩阵(并返回向量)时，我没有收到错误.但是，我将上述功能作为主程序的一部分，即

I couldn't quite figure it out. I didn't get the error when the input of the function was a matrix (and returned a vector). However, I included the above function as part of my main program i.e.

#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

该程序已成功编译，但运行速度非常慢(几乎与使用colMeans的R版本一样慢).有没有更好的方法来计算多维数据集上的列均值，为什么会出现该编译错误?

The program compiled successfully, but it is painfully slow (almost as slow as the R version of the program using colMeans). Is there a better way to calculate column means on a cube, and why am I getting that compilation error?

我将不胜感激.

此致

推荐答案

在尝试将arma::cube用作Rcpp函数参数时，我也收到此错误. 基于编译器错误，我认为这是因为当前没有定义Rcpp::wrap<arma::cube>(处理传递给该函数的R对象需要使用它).†阅读了几篇相关的文章在网上的示例中，典型的解决方法似乎是将R array读为NumericVector，并且由于它保留了其dims属性，请使用这些设置您的arma::cube尺寸.尽管缺少wrap专业化 †需要花一两个步骤，但我放进来的Armadillo版本似乎比我的R解决方案快很多:

I also received this error when attempting to use an arma::cube as an Rcpp function parameter. Based on the compiler error, I believe this is because there is no Rcpp::wrap<arma::cube> currently defined (which is needed to handle the R object you would pass to the function).† After reading a couple of related examples online, it looks like the typical workaround is to read in your R array as a NumericVector, and since it retains its dims attribute, use these to set your arma::cube dimensions. Despite the fact that there is an extra step or two required to account for the missing wrap specialization†, the Armadillo version I put together seems to be quite a bit faster than my R solution:

#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

*/

我利用了arma::mat的arma::mean函数重载将默认计算列均值的事实(arma::mean(x.slice(i), 1)将为您提供该切片的行均值).

where I am taking advantage of the fact that the arma::mean function overload for arma::mats will calculate column means by default (arma::mean(x.slice(i), 1) would give you the row means of that slice).

†再三考虑，我不确定这是否与Rcpp::wrap有关-但问题似乎与缺少Exporter<>专业化有关arma::cube-Rcpp的Exporter.h的第31行:

† On second thought, I'm not really sure if this has to do with Rcpp::wrap or not - but the issue seems to be related to a missing Exporter<> specialization for arma::cube - line 31 of Rcpp's Exporter.h:

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

private:
  T t ;
} ;

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

Regardless, NumericVector / setting dimensions approach I used seems to be functional solution for now.

根据您在问题中描述的输出维度，我假设您希望所得矩阵的每一列都是对应数组切片的列均值的向量(列1 =切片1的列均值，依此类推... )，即

Based on the output dimensions you described in your question, I assumed you wanted each column of the resulting matrix to be a vector of column means of the corresponding array slice (column 1 = column means of slice 1, etc...), i.e.

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

但是，如果需要的话，对您来说微不足道.

but it would be trivial for you to alter this if needed.

这篇关于列均值3d矩阵(立方体)Rcpp的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！