RcppArmadillo中稀疏和密集矩陣的模板函數

Question

我試圖定義一個模板函數，它可以使用RcppArmadillo處理稀疏和密集矩陣輸入。我發送一個密集或稀疏矩陣C++和回R的非常簡單的情況下，像這樣的工作：

library(inline); library(Rcpp); library(RcppArmadillo) 

sourceCpp(code = " 
// [[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 
using namespace Rcpp ; 
using namespace arma ; 

template <typename T> T importexport_template(const T X) { 
    T ret = X ; 
    return ret ; 
}; 

//[[Rcpp::export]] 
SEXP importexport(SEXP X) { 
    return wrap(importexport_template(X)) ; 
}") 

library(Matrix) 
X <- diag(3) 
X_sp <- as(X, "dgCMatrix") 

importexport(X) 
##  [,1] [,2] [,3] 
##[1,] 1 0 0 
##[2,] 0 1 0 
##[3,] 0 0 1 
importexport(X_sp) 
##3 x 3 sparse Matrix of class "dgCMatrix" 
##   
##[1,] 1 . . 
##[2,] . 1 . 
##[3,] . . 1 

，我理解這意味着模板基本工作原理（即密集R-矩陣通過對Rcpp::as的隱式調用，以及相應的暗指Rcpp:wrap然後做正確的事情並返回稠密以用於稀疏稀疏，稀疏R矩陣變爲arma::sp_mat-對象。

的實際功能我嘗試寫需求，當然多個參數，而這也正是我失敗 - 做這樣的事情：

sourceCpp(code =" 
// [[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 

using namespace Rcpp ; 
using namespace arma ; 

template <typename T> T scalarmult_template(const T X, double scale) { 
    T ret = X * scale; 
    return ret; 
}; 

//[[Rcpp::export]] 
SEXP scalarmult(SEXP X, double scale) { 
    return wrap(scalarmult_template(X, scale)) ; 
}") 

失敗，因爲編譯器不知道如何在編譯時解決*爲SEXPREC* const。 所以我想我需要的東西，如開關語句in this Rcpp Gallery snippet正確分派到特定的模板功能，但我不知道該怎麼寫了，似乎比INTSXP等

比較複雜，我想我知道如何種訪問類型我需要這樣一個switch語句，例如：

sourceCpp(code =" 
// [[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 

using namespace Rcpp ; 
using namespace arma ; 

//[[Rcpp::export]] 
SEXP printtype(SEXP Xr) { 
    Rcpp::Rcout << TYPEOF(Xr) << std::endl ; 
    return R_NilValue; 
}") 
printtype(X) 
##14 
##NULL 
printtype(X_sp) 
##25 
##NULL 

但我不明白如何從那裏繼續。 scalarmult_template適用於稀疏矩陣和密集矩陣的版本是什麼樣的？

Answer 1

基於@ KevinUshey的評論回答我自己的問題。我做矩陣乘法爲三種情況：密密，疏密，以及「indMatrix」 -dense：

library(inline) 
library(Rcpp) 
library(RcppArmadillo) 
library(Matrix) 
library(rbenchmark) 

sourceCpp(code=" 
// [[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 

using namespace Rcpp ; 
using namespace arma ; 

arma::mat matmult_sp(const arma::sp_mat X, const arma::mat Y){ 
    arma::mat ret = X * Y; 
    return ret; 
}; 
arma::mat matmult_dense(const arma::mat X, const arma::mat Y){ 
    arma::mat ret = X * Y; 
    return ret; 
}; 
arma::mat matmult_ind(const SEXP Xr, const arma::mat Y){ 
    // pre-multplication with index matrix is a permutation of Y's rows: 
    S4 X(Xr); 
    arma::uvec perm = X.slot("perm"); 
    arma::mat ret = Y.rows(perm - 1); 
    return ret; 
}; 

//[[Rcpp::export]] 
arma::mat matmult_cpp(SEXP Xr, const arma::mat Y) { 
    if (Rf_isS4(Xr)) { 
     if(Rf_inherits(Xr, "dgCMatrix")) { 
      return matmult_sp(as<arma::sp_mat>(Xr), Y) ; 
     } ; 
     if(Rf_inherits(Xr, "indMatrix")) { 
      return matmult_ind(Xr, Y) ; 
     } ; 
     stop("unknown class of Xr") ; 
    } else { 
     return matmult_dense(as<arma::mat>(Xr), Y) ; 
    } 
}") 

n <- 10000 
d <- 20 
p <- 30 

X <- matrix(rnorm(n*d), n, d) 
X_sp <- as(diag(n)[,1:d], "dgCMatrix") 
X_ind <- as(sample(1:d, n, rep=TRUE), "indMatrix") 
Y <- matrix(1:(d*p), d, p) 

matmult_cpp(as(X_ind, "ngTMatrix"), Y) 
## Error: unknown class of Xr 

all.equal(X%*%Y, matmult_cpp(X, Y)) 
## [1] TRUE 

all.equal(as.vector(X_sp%*%Y), 
      as.vector(matmult_cpp(X_sp, Y))) 
## [1] TRUE 

all.equal(X_ind%*%Y, matmult_cpp(X_ind, Y)) 
## [1] TRUE 

---

編輯：這已經變成一個Rcpp Gallery post。