在 C++ 代码中使用列表作为输入并使用 Rcpp 调用(列表输入非常慢) [英] Using list as an input in C++ code and calling using Rcpp(List inputs are extremely slow)
本文介绍了在 C++ 代码中使用列表作为输入并使用 Rcpp 调用(列表输入非常慢)的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我正在尝试使用列表(R 对象)作为 C++ 函数的输入，然后使用 R 中的 Rcpp 调用它.此列表包含大量矩阵.我提供的代码是一个玩具示例.我已经编写了一个非常复杂的代码，但效率非常低.在下面的代码中，我想知道是否有一种从列表中提取矩阵的有效方法.
以下是我尝试过的代码.它有效，但它也告诉我下标值不是数组、指针或向量.我正在使用 R studio 编写此代码.当我编译代码时，它可以工作，但是当我将鼠标光标放在编辑器中时，我也看到了红色叉号下标值不是数组、指针或向量".
#include 使用命名空间 Rcpp；//[[Rcpp::depends(RcppArmadillo)]]//[[Rcpp::export]]Rcpp::List tt(Rcpp::List&H T，输入 n){列表 A(n);for(int i=0;i<n;++i){arma::mat htt = ht[i];//这是我看到下标值不是数组、指针或向量的地方arma::mat x = htt * htt.t();A[i] = x;//这是我看到下标值不是数组、指针或向量的地方}列表资源(1);res[0] = A;//这是我看到下标值不是数组、指针或向量的地方返回(资源)；} 同样，这是一个可以在 R 中轻松完成的玩具示例.我想了解如何有效地完成此操作.假设，我希望列表的每个矩阵都乘以它的转置.任何帮助，将不胜感激?以下是我的实际问题.
#include 使用命名空间 Rcpp；//[[Rcpp::depends(RcppArmadillo)]]//[[Rcpp::export]]列表 se_4a(Rcpp::List&H T，const int&n,const int&p,const int&个人电脑，数字矩阵&S1byS0_,数字矩阵&S1byS0c_,数字矩阵&za_,数字矩阵&zb_,数字矩阵&wd_,数字矩阵&一世_，NumericVector&S0c_,NumericVector&伽玛_){列表 A(ht.length());arma::mat S1byS0hat(S1byS0_.begin(),S1byS0_.nrow(),S1byS0_.ncol(),false);arma::mat S1byS0hatc(S1byS0c_.begin(),S1byS0c_.nrow(),S1byS0c_.ncol(),false);arma::mat z(za_.begin(),za_.nrow(),za_.ncol(),false);arma::mat zc(zb_.begin(),zb_.nrow(),zb_.ncol(),false);arma::mat wdM(wd_.begin(),wd_.nrow(),wd_.ncol(),false);arma::mat Ic(I_.begin(),I_.nrow(),I_.ncol(),false);arma::vec S0hatc(S0c_.begin(),S0c_.size(),false);arma::vec gammahat(gammah_.begin(),gammah_.size(),false);Rcpp::List q1hat(n);Rcpp::List q2hat(n);for(int i=0; i < n;++i){arma::mat q11hat(p,n);q11hat.zeros();arma::mat q21hat(p,n);q21hat.zeros();//arma::mat q11hat(q11hata.begin(),q11hata.nrow(),q11hata.ncol(),false);//arma::mat q21hat(q21hata.begin(),q21hata.nrow(),q21hata.ncol(),false);for(int u = 0;u (ht[k]);Rcpp::NumericMatrix htt_R = ht[k];arma::mat htt(htt_R.begin(), htt_R.rows(), htt_R.cols(), false, true);//arma::vec y = httt(_,j);arma::colvec y = htt.col(j);arma::rowvec yy = y.t() * Ic * (zc.row(i).t() - S1byS0hatc.col(u));双 zz = yy(0,0);q += (z.row(j).t() - S1byS0hat.col(k)) *zz * wdM(j,k);如果 (u <= k){qq += (z.row(j).t() - S1byS0hat.col(k)) *exp(arma::as_scalar(gammahat.t()*zc.row(j).t()))*wdM(j,k)/(S0hatc(u)/n);}}}}q11hat.cols(u,u) = -1 * q;q21hat.cols(u,u) = -1 * qq;}q1hat[i] = q11hat/n;q2hat[i] = q21hat/n;}return List::create(Named("A")=q1hat,Named("B")=q2hat);} 现在，让我们从 R 中调用上述函数.
#Calling from Rht <- 列表()for(i in 1 : 100){ht[[i]] <-matrix(runif(10*100), 10, 100)}n <- 100p＜- 10pc <- 10S1byS0 <- 矩阵(norm(10*100),10,100)S1byS0c <- S1byS0za <- 矩阵(rnorm(100*10),100,10)zb <- zawd <- 矩阵(rnorm(100*100),100,100)I <- 矩阵(norm(100),10,10)S0c <- c(norm(100))伽马 <- 矩阵(rnorm(10)，1,10)#调用se_4a函数pp=bench::mark(se_4a(ht,n,p,pc, S1byS0,S1byS0c,za,zb,wd,I,S0c,gammah))bench::mark(se_4a(ht,n,p,pc, S1byS0,S1byS0c,za,zb,wd,I,S0c,gammah))
# tibble: 1 x 13表达式最小中位数`itr/sec`<bch:expr><bch><bch:><dbl>1 se_4a(ht, n, p, pc, S1byS0, S1byS0c, za, zb, wd, I, S0c, gammah) 20.9s 20.9s 0.0479# ... 还有 9 个变量:mem_alloc 、`gc/sec` 、n_itr 、n_gc 、# total_time <bch:tm>, result <list>, memory <list>, time <list>, gc <list>警告信息:有些表达式在每次迭代中都有一个 GC；所以过滤被禁用. 使用矩阵作为参数而不是列表.
#include 使用命名空间 Rcpp；//[[Rcpp::depends(RcppArmadillo)]]//[[Rcpp::export]]列表 se_4c(const arma::mat&H T，const int&n,const int&p,const int&个人电脑，const arma::mat&S1byS0hat，const arma::mat&S1byS0hatc，const arma::mat&z,const arma::mat&zc,const arma::mat&wdm,const arma::mat&我知道了，const arma::vec&S0hatc,const arma::vec&伽玛哈特){Rcpp::List q1hat(n);Rcpp::List q2hat(n);arma::mat q11hat(p,n);arma::mat q21hat(p,n);arma::mat q(p,1);arma::mat qq(p,1);std::vector 解决方案 几点意见:
RcppArmadillo 可以将 Armadillo 的高级构造函数直接用于函数参数.我不确定确切的条件，但 const 引用被这样处理.
从 Rcpp::List 中提取 Rcpp::NumericMatrix 似乎代价高昂.因此，一开始就对矩阵进行一次处理是有意义的.
计算的几个部分与内部循环变量无关.在循环之前计算一次它们是有意义的(参见下面的 foo、bar 和 baz).
我也做了一些清理工作.
替代功能:
//[[Rcpp::export]]列表 se_4b(Rcpp::List&H T，const int&n,const int&p,const int&个人电脑，const arma::mat&S1byS0hat，const arma::mat&S1byS0hatc，const arma::mat&z,const arma::mat&zc,const arma::mat&wdm,const arma::mat&我知道了，const arma::vec&S0hatc,const arma::vec&伽玛哈特){Rcpp::List q1hat(n);Rcpp::List q2hat(n);arma::mat q11hat(p,n);arma::mat q21hat(p,n);arma::mat q(p,1);arma::mat qq(p,1);std::vectorhtt_vec(n);for(int i = 0; i < n; ++i) {Rcpp::NumericMatrix htt_R = ht[i];arma::mat htt(htt_R.begin(), htt_R.rows(), htt_R.cols(), false, true);htt_vec[i] = htt;}for(int i=0; i < n;++i){for(int u = 0;u
与您的函数比较的基准结果:
<代码>>基准::标记(se_4a = se_4a(ht,n,p,pc, S1byS0,S1byS0c,za,zb,wd,I,S0c,gammah),+ se_4b = se_4b(ht,n,p,pc, S1by .... [截断]# 小块:2 x 13表达式 min 中值 `itr/sec` mem_alloc `gc/sec` n_itr n_gc total_time 结果<bch:expr><bch:><bch:><dbl><bch:字节><dbl><int><dbl><bch:tm><列表>1 se_4a 21.97s 21.97s 0.0455 1.54MB 7.97 1 175 21.97s <list…2 se_4b 4.84s 4.84s 0.206 1.54MB 0 1 0 4.84s <列表...# ... 还有 3 个变量:内存 <list>、时间 <list>、gc <list>
I am trying to use a list(R object) as an input for C++ function and later call it using Rcpp from R. This list contains large number of matrices. The code I provide is a toy example. I have a very complicated code that I have already written but very inefficient. In the following code, I want to know if there is an efficient way of extracting matrix from the list.
Following is the code that I have tried. It works but it also tells me that the subscripted value is not an array, pointer or vector.I am using R studio to write this code. When I compile the code it works, but when I put mouse cursor in the editor I also see the red cross saying "subscripted value is not an array, pointer or vector".
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List tt(
Rcpp::List& ht,
int n){
List A(n);
for(int i=0;i<n;++i){
arma::mat htt = ht[i];// this is where I see subscripted value is not an array, pointer or vector
arma::mat x = htt * htt.t();
A[i] = x;//this is where I see subscripted value is not an array, pointer or vector
}
List res(1);
res[0] = A;//this is where I see subscripted value is not an array, pointer or vector
return(res);
}
Again, this is a toy example which could easily be done in R. I would like to get some idea on how this could be done efficiently. Suppose, I want every matrix of a list to be multiplied by it's transpose. Any help would be appreciated? Following is my actual problem.
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
List se_4a(
Rcpp::List& ht,
const int& n,
const int& p,
const int& pc,
NumericMatrix& S1byS0_,
NumericMatrix& S1byS0c_,
NumericMatrix& za_,
NumericMatrix& zb_,
NumericMatrix& wd_,
NumericMatrix& I_,
NumericVector& S0c_,
NumericVector& gammah_){
List A(ht.length());
arma::mat S1byS0hat(S1byS0_.begin(),S1byS0_.nrow(),S1byS0_.ncol(),false);
arma::mat S1byS0hatc(S1byS0c_.begin(),S1byS0c_.nrow(),S1byS0c_.ncol(),false);
arma::mat z(za_.begin(),za_.nrow(),za_.ncol(),false);
arma::mat zc(zb_.begin(),zb_.nrow(),zb_.ncol(),false);
arma::mat wdM(wd_.begin(),wd_.nrow(),wd_.ncol(),false);
arma::mat Ic(I_.begin(),I_.nrow(),I_.ncol(),false);
arma::vec S0hatc(S0c_.begin(),S0c_.size(),false);
arma::vec gammahat(gammah_.begin(),gammah_.size(),false);
Rcpp::List q1hat(n);
Rcpp::List q2hat(n);
for(int i=0; i < n;++i){
arma::mat q11hat(p,n);
q11hat.zeros();
arma::mat q21hat(p,n);
q21hat.zeros();
// arma::mat q11hat(q11hata.begin(),q11hata.nrow(),q11hata.ncol(),false);
// arma::mat q21hat(q21hata.begin(),q21hata.nrow(),q21hata.ncol(),false);
for(int u = 0;u < n; ++u){
// arma::mat q(qa.begin(),qa.nrow(),qa.ncol(),false);
// arma::mat qq(qqa.begin(),qqa.nrow(),qqa.ncol(),false);
arma::mat q(p,1);
q.zeros();
arma::mat qq(p,1);
qq.zeros();
for(int j=0;j <n;++j){
if(j < n){
for(int k = j; k <n;++k){
//NumericMatrix httt = as<NumericMatrix>(ht[k]);
Rcpp::NumericMatrix htt_R = ht[k];
arma::mat htt(htt_R.begin(), htt_R.rows(), htt_R.cols(), false, true);
//arma::vec y = httt(_,j);
arma::colvec y = htt.col(j);
arma::rowvec yy = y.t() * Ic * (zc.row(i).t() - S1byS0hatc.col(u));
double zz = yy(0,0);
q += (z.row(j).t() - S1byS0hat.col(k)) *
zz * wdM(j,k);
if (u <= k){
qq += (z.row(j).t() - S1byS0hat.col(k)) *
exp(arma::as_scalar(gammahat.t()*zc.row(j).t()))*wdM(j,k) / (S0hatc(u)/n);
}
}
}
}
q11hat.cols(u,u) = -1 * q;
q21hat.cols(u,u) = -1 * qq;
}
q1hat[i] = q11hat/n;
q2hat[i] = q21hat/n;
}
return List::create(Named("A")=q1hat,
Named("B")=q2hat);
}
Now, let's call the above function from R.
#Calling from R
ht <- list()
for(i in 1 : 100){
ht[[i]] <-matrix(runif(10*100), 10, 100)
}
n <- 100
p <- 10
pc <- 10
S1byS0 <- matrix(rnorm(10*100),10,100)
S1byS0c <- S1byS0
za <- matrix(rnorm(100*10),100,10)
zb <- za
wd <- matrix(rnorm(100*100),100,100)
I <- matrix(rnorm(100),10,10)
S0c <- c(rnorm(100))
gammah <- matrix(rnorm(10),1,10)
#Calling se_4a function
pp=bench::mark(se_4a(ht,n,p,pc, S1byS0,S1byS0c,za,zb,wd,I,S0c,gammah))
bench::mark(se_4a(ht,n,p,pc, S1byS0,S1byS0c,za,zb,wd,I,S0c,gammah))
# A tibble: 1 x 13
expression min median `itr/sec`
<bch:expr> <bch> <bch:> <dbl>
1 se_4a(ht, n, p, pc, S1byS0, S1byS0c, za, zb, wd, I, S0c, gammah) 20.9s 20.9s 0.0479
# … with 9 more variables: mem_alloc <bch:byt>, `gc/sec` <dbl>, n_itr <int>, n_gc <dbl>,
# total_time <bch:tm>, result <list>, memory <list>, time <list>, gc <list>
Warning message:
Some expressions had a GC in every iteration; so filtering is disabled.
Using matrix as an argument instead of a list.
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
List se_4c(
const arma::mat& ht,
const int& n,
const int& p,
const int& pc,
const arma::mat& S1byS0hat,
const arma::mat& S1byS0hatc,
const arma::mat& z,
const arma::mat& zc,
const arma::mat& wdM,
const arma::mat& Ic,
const arma::vec& S0hatc,
const arma::vec& gammahat){
Rcpp::List q1hat(n);
Rcpp::List q2hat(n);
arma::mat q11hat(p,n);
arma::mat q21hat(p,n);
arma::mat q(p,1);
arma::mat qq(p,1);
std::vector<arma::mat> htt_vec(n);
for(int i = 0; i < n; ++i) {
// Rcpp::NumericMatrix htt_R = ht[i];
// arma::mat htt(htt_R.begin(), htt_R.rows(), htt_R.cols(), false, true);
// htt_vec[i] = htt;
htt_vec[i] = ht.rows(i,i+(p-1));
}
for(int i=0; i < n;++i){
for(int u = 0;u < n; ++u){
q.zeros();
qq.zeros();
arma::mat bar = Ic * (zc.row(i).t() - S1byS0hatc.col(u));
for(int j=0;j <n;++j){
if(j < n){
double foo = exp(arma::as_scalar(gammahat.t()*zc.row(j).t())) / (S0hatc(u)/n);
for(int k = j; k <n;++k){
//arma::mat htt_vec = ht.rows(k,k+(p-1));
arma::colvec y = htt_vec[k].col(j);
arma::rowvec yy = y.t() * bar;
double zz = yy(0,0);
arma::mat baz = (z.row(j).t() - S1byS0hat.col(k)) * wdM(j,k);
q += zz * baz;
if (u <= k){
qq += foo * baz;
}
}
}
}
q11hat.col(u) = -q;
q21hat.col(u) = -qq;
}
q1hat[i] = q11hat/n;
q2hat[i] = q21hat/n;
}
return List::create(Named("A")=q1hat,
Named("B")=q2hat);
}
Following is the simple code that takes about 4 secs. Expected it to faster.
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
List test(
List& ht,
const int& n,
const int& p){
Rcpp::List q1hat(n);
Rcpp::List q2hat(n);
for(int i=0; i < n;++i){
arma::mat q11hat(p,n);
q11hat.zeros();
arma::mat q21hat(p,n);
q21hat.zeros();
for(int u = 0;u < n; ++u){
arma::mat q(p,1);
q.zeros();
arma::mat qq(p,1);
qq.zeros();
for(int j=0;j <n;++j){
if(j < ht.length()){
for(int k = j; k <n;++k){
Rcpp::NumericMatrix htt_R = ht[k];
arma::mat htt(htt_R.begin(), htt_R.rows(), htt_R.cols(), false, true);
}
}
}
q11hat.cols(u,u) = -1 * q;
q21hat.cols(u,u) = -1 * qq;
}
q1hat[i] = q11hat/n;
q2hat[i] = q21hat/n;
}
return List::create(Named("A")=q1hat,
Named("B")=q2hat);
}
The formula that I am trying to implement is as follows:
解决方案
Several comments:
RcppArmadillo can make use of Armadillo's advanced constructors directly for function arguments. I am not sure about the precise condition, but
constreferences get treated like this.Extracting a
Rcpp::NumericMatrixfrom aRcpp::Listseems to be costly. It therefore makes sense to process the matrices once in the beginning.Several parts of your calculations where independent of the inner loop variables. It makes sense to calculate them once before the loop (see
foo,barandbazbelow).I have also done some clean-up.
Alternative function:
// [[Rcpp::export]]
List se_4b(
Rcpp::List& ht,
const int& n,
const int& p,
const int& pc,
const arma::mat& S1byS0hat,
const arma::mat& S1byS0hatc,
const arma::mat& z,
const arma::mat& zc,
const arma::mat& wdM,
const arma::mat& Ic,
const arma::vec& S0hatc,
const arma::vec& gammahat){
Rcpp::List q1hat(n);
Rcpp::List q2hat(n);
arma::mat q11hat(p,n);
arma::mat q21hat(p,n);
arma::mat q(p,1);
arma::mat qq(p,1);
std::vector<arma::mat> htt_vec(n);
for(int i = 0; i < n; ++i) {
Rcpp::NumericMatrix htt_R = ht[i];
arma::mat htt(htt_R.begin(), htt_R.rows(), htt_R.cols(), false, true);
htt_vec[i] = htt;
}
for(int i=0; i < n;++i){
for(int u = 0;u < n; ++u){
q.zeros();
qq.zeros();
arma::mat bar = Ic * (zc.row(i).t() - S1byS0hatc.col(u));
for(int j=0;j <n;++j){
if(j < n){
double foo = exp(arma::as_scalar(gammahat.t()*zc.row(j).t())) / (S0hatc(u)/n);
for(int k = j; k <n;++k){
arma::colvec y = htt_vec[k].col(j);
arma::rowvec yy = y.t() * bar;
double zz = yy(0,0);
arma::mat baz = (z.row(j).t() - S1byS0hat.col(k)) * wdM(j,k);
q += zz * baz;
if (u <= k){
qq += foo * baz;
}
}
}
}
q11hat.col(u) = -q;
q21hat.col(u) = -qq;
}
q1hat[i] = q11hat/n;
q2hat[i] = q21hat/n;
}
return List::create(Named("A")=q1hat,
Named("B")=q2hat);
}
Benchmark result comparing with your function:
> bench::mark(se_4a = se_4a(ht,n,p,pc, S1byS0,S1byS0c,za,zb,wd,I,S0c,gammah),
+ se_4b = se_4b(ht,n,p,pc, S1by .... [TRUNCATED]
# A tibble: 2 x 13
expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc total_time result
<bch:expr> <bch:> <bch:> <dbl> <bch:byt> <dbl> <int> <dbl> <bch:tm> <list>
1 se_4a 21.97s 21.97s 0.0455 1.54MB 7.97 1 175 21.97s <list…
2 se_4b 4.84s 4.84s 0.206 1.54MB 0 1 0 4.84s <list…
# … with 3 more variables: memory <list>, time <list>, gc <list>
这篇关于在 C++ 代码中使用列表作为输入并使用 Rcpp 调用(列表输入非常慢)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
