2015-06-21 70 views
2

我有一個程序,我需要反覆計算Rcpp中立方體X(nRow, nCol, nSlice)的每個切片的列平均值,所得到的平均值形成矩陣M(nCol, nSlice)。下面的代碼產生的錯誤:列意味着3d矩陣(立方體)Rcpp

#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: error: no matching function for call to '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 

成功編譯該程序的一部分,但它是非常緩慢(幾乎爲使用colMeans程序中的R版本慢)。有沒有更好的方法來計算多維數據集中的列方法,以及爲什麼會得到該編譯錯誤?

我很感激任何幫助。

問候,

+1

幾乎和colMeans一樣慢? ...就是這樣。 'colMeans'在C語言中進行了優化。您缺少提供用於運行此程序的代碼可能包含回答問題的線索。順便說一下,R​​中沒有「3d矩陣」。陣列將是另一回事。術語通常很重要。 –

回答

5

嘗試使用arma::cube作爲RCPP函數參數的時候我也收到此錯誤。 基於編譯器錯誤,我相信這是因爲當前沒有定義 Rcpp::wrap<arma::cube>(這是處理要傳遞給函數的R對象所需的)。 †在網上閱讀了幾個相關示例之後,看起來典型的解決方法是將R array的內容讀取爲NumericVector,並且由於它保留了其dims屬性,請使用這些屬性設置arma::cube尺寸。儘管有以佔 失蹤 wrap專業化 †需要一個額外的步驟或兩個,犰狳版本我放在一起似乎是相當快一點比我的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 

*/ 

在哪裏我正在利用arma::matarma::mean函數過載將默認計算列的方式這一事實(arma::mean(x.slice(i), 1)會爲您提供該片的行手段)。


編輯:†關於第二個想法,我真的不知道這是否與Rcpp::wrap或沒有做的事 - 但問題似乎與缺少Exporter<>專業化的arma::cube - 行31 RCPP的Exporter.h:

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 

但是如果你需要改變它,這將是微不足道的。