如何在 Rcpp 中对 bigstatsr::FBM 的行或列进行子集并将它们存储在向量中?
How can I subset rows or columns of a bigstatsr::FBM in Rcpp and store them in a vector?
我有一个函数可以根据给定矩阵的行(或列)计算基本汇总统计信息,我现在正尝试将此函数与 bigstatsr::FBM 一起使用(我知道使用列应该更有效率)。
我想将行/列存储在向量中的原因是我想用 std::nth_element 计算分位数。如果有一种不同的方法可以在不使用矢量的情况下做到这一点,我会同样高兴。
这是我用于正则矩阵的代码。
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Rcpp;
// [[Rcpp::export]]
Eigen::MatrixXd summaryC(Eigen::MatrixXd x,int nrow) {
Eigen::MatrixXd result(nrow, 5);
int indices[6] = {-1, 0, 249, 500, 750, 999};
for (int i = 0; i < nrow; i++) {
Eigen::VectorXd v = x.row(i);
for (int q = 0; q < 5; ++q) {
std::nth_element(v.data() + indices[q] + 1,
v.data() + indices[q+1],
v.data() + v.size());
result(i,q) = v[indices[q+1]];
}
}
return result;
}
/*** R
x <- matrix(as.numeric(1:1000000), ncol = 1000)
summaryC(x = x, nrow = 1000)
***/
但是我很难用 FBM 做到这一点,因为我没有完全掌握 FBM - Pointer 工作原理的复杂性。
我尝试了以下但没有成功:
// [[Rcpp::depends(BH, bigstatsr, RcppEigen)]]
// [[Rcpp::plugins(cpp11)]]
#include <bigstatsr/BMAcc.h>
#include <RcppEigen.h>
// [[Rcpp::export]]
Eigen::MatrixXd summaryCbig(Environment fbm,int nrow, Eigen::VecttorXi ind_col) {
Eigen::MatrixXd result(nrow, 5);
XPtr<FBM> xpMat = fbm["address"];
BMAcc<double> macc(xpMat);
int indices[6] = {-1, 0, 249, 500, 750, 999};
for (int i = 0; i < nrow; i++) {
Eigen::VectorXd v = macc.row(i); // this does not work
Eigen::VectorXd v = macc(i,_); // this does not work
SubBMAcc<double> maccr(XPtr, i, ind_col -1); // This did not work with Eigen::VectorXi, but works with const NumericVector&
Eigen::VectorXd v = maccr // this does not work even for appropriate ind_col
for (int q = 0; q < 5; ++q) {
std::nth_element(v.data() + indices[q] + 1,
v.data() + indices[q+1],
v.data() + v.size());
macc(i,q) = v[indices[q+1]];
}
}
}
/*** R
x <- matrix(as.numeric(1:1000000), ncol = 1000)
summaryCbig(x = x, nrow = 1000, ind_col = 1:1000)
***/
任何帮助将不胜感激,谢谢!
更新 - big_apply - 方法
我用两个不同大小的矩阵 X1 和 X2 实施了两次该方法。 X1 的代码:
X1 <- FBM(1000, 1000, init 1e6)
X2 <- FBM(10000, 10000, init = 9999)
library(bigstatsr)
microbenchmark::microbenchmark(
big_apply(X, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X1[ind, ])
}, a.combine = "rbind", ind = rows_along(X), ncores = nb_cores(), block.size = 500),
big_apply(X, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X1[ind, ])
}, a.combine = "rbind", ind = rows_along(X), ncores = 1, block.size = 500),
times = 5
)
当使用 X1 和 block.size = 500 时,有 4 个内核而不是 1 个内核会使我的 PC 上的任务慢 5-10 倍(不幸的是,4 CPU 和使用 windows ).
使用更大的矩阵 X2 并将 block.size 保留为默认值,使用 4 个内核而不是非并行版本需要 10 倍的时间。
X2 的结果:
min lq mean median uq max neval
16.149055 19.13568 19.369975 20.139363 20.474103 20.951676 5
1.297259 2.67385 2.584647 2.858035 2.867537 3.226552 5
假设你有
library(bigstatsr)
X <- FBM(1000, 1000, init = 1:1e6)
我不会重新发明轮子并使用:
big_apply(X, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X[ind, ])
}, a.combine = "rbind", ind = rows_along(X), ncores = nb_cores(), block.size = 500)
明智地选择block.size
(行数)。
如果您想将 R(cpp) 函数应用于 FBM
.
的块,函数 big_apply()
非常有用
编辑: 当然,对于小矩阵,并行性会变慢,因为并行性的开销(通常为 1-3 秒)。查看 X1 和 X2 的结果:
library(bigstatsr)
X1 <- FBM(1000, 1000, init = 1e6)
microbenchmark::microbenchmark(
PAR = big_apply(X1, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X[ind, ])
}, a.combine = "rbind", ind = rows_along(X1), ncores = nb_cores(), block.size = 500),
SEQ = big_apply(X1, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X[ind, ])
}, a.combine = "rbind", ind = rows_along(X1), ncores = 1, block.size = 500),
times = 5
)
Unit: milliseconds
expr min lq mean median uq max neval cld
PAR 1564.20591 1602.0465 1637.77552 1629.9803 1651.04509 1741.59974 5 b
SEQ 68.92936 69.1002 76.70196 72.9173 85.31751 87.24543 5 a
X2 <- FBM(10000, 10000, init = 9999)
microbenchmark::microbenchmark(
PAR = big_apply(X2, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X[ind, ])
}, a.combine = "rbind", ind = rows_along(X2), ncores = nb_cores(), block.size = 500),
SEQ = big_apply(X2, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X[ind, ])
}, a.combine = "rbind", ind = rows_along(X2), ncores = 1, block.size = 500),
times = 5
)
Unit: seconds
expr min lq mean median uq max neval cld
PAR 4.757409 4.958869 5.071982 5.083381 5.218098 5.342153 5 a
SEQ 10.842828 10.846281 11.177460 11.360162 11.416967 11.421065 5 b
您的矩阵越大,您从并行性中获得的收益就越大。
我有一个函数可以根据给定矩阵的行(或列)计算基本汇总统计信息,我现在正尝试将此函数与 bigstatsr::FBM 一起使用(我知道使用列应该更有效率)。 我想将行/列存储在向量中的原因是我想用 std::nth_element 计算分位数。如果有一种不同的方法可以在不使用矢量的情况下做到这一点,我会同样高兴。
这是我用于正则矩阵的代码。
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Rcpp;
// [[Rcpp::export]]
Eigen::MatrixXd summaryC(Eigen::MatrixXd x,int nrow) {
Eigen::MatrixXd result(nrow, 5);
int indices[6] = {-1, 0, 249, 500, 750, 999};
for (int i = 0; i < nrow; i++) {
Eigen::VectorXd v = x.row(i);
for (int q = 0; q < 5; ++q) {
std::nth_element(v.data() + indices[q] + 1,
v.data() + indices[q+1],
v.data() + v.size());
result(i,q) = v[indices[q+1]];
}
}
return result;
}
/*** R
x <- matrix(as.numeric(1:1000000), ncol = 1000)
summaryC(x = x, nrow = 1000)
***/
但是我很难用 FBM 做到这一点,因为我没有完全掌握 FBM - Pointer 工作原理的复杂性。
我尝试了以下但没有成功:
// [[Rcpp::depends(BH, bigstatsr, RcppEigen)]]
// [[Rcpp::plugins(cpp11)]]
#include <bigstatsr/BMAcc.h>
#include <RcppEigen.h>
// [[Rcpp::export]]
Eigen::MatrixXd summaryCbig(Environment fbm,int nrow, Eigen::VecttorXi ind_col) {
Eigen::MatrixXd result(nrow, 5);
XPtr<FBM> xpMat = fbm["address"];
BMAcc<double> macc(xpMat);
int indices[6] = {-1, 0, 249, 500, 750, 999};
for (int i = 0; i < nrow; i++) {
Eigen::VectorXd v = macc.row(i); // this does not work
Eigen::VectorXd v = macc(i,_); // this does not work
SubBMAcc<double> maccr(XPtr, i, ind_col -1); // This did not work with Eigen::VectorXi, but works with const NumericVector&
Eigen::VectorXd v = maccr // this does not work even for appropriate ind_col
for (int q = 0; q < 5; ++q) {
std::nth_element(v.data() + indices[q] + 1,
v.data() + indices[q+1],
v.data() + v.size());
macc(i,q) = v[indices[q+1]];
}
}
}
/*** R
x <- matrix(as.numeric(1:1000000), ncol = 1000)
summaryCbig(x = x, nrow = 1000, ind_col = 1:1000)
***/
任何帮助将不胜感激,谢谢!
更新 - big_apply - 方法
我用两个不同大小的矩阵 X1 和 X2 实施了两次该方法。 X1 的代码:
X1 <- FBM(1000, 1000, init 1e6)
X2 <- FBM(10000, 10000, init = 9999)
library(bigstatsr)
microbenchmark::microbenchmark(
big_apply(X, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X1[ind, ])
}, a.combine = "rbind", ind = rows_along(X), ncores = nb_cores(), block.size = 500),
big_apply(X, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X1[ind, ])
}, a.combine = "rbind", ind = rows_along(X), ncores = 1, block.size = 500),
times = 5
)
当使用 X1 和 block.size = 500 时,有 4 个内核而不是 1 个内核会使我的 PC 上的任务慢 5-10 倍(不幸的是,4 CPU 和使用 windows ). 使用更大的矩阵 X2 并将 block.size 保留为默认值,使用 4 个内核而不是非并行版本需要 10 倍的时间。
X2 的结果:
min lq mean median uq max neval
16.149055 19.13568 19.369975 20.139363 20.474103 20.951676 5
1.297259 2.67385 2.584647 2.858035 2.867537 3.226552 5
假设你有
library(bigstatsr)
X <- FBM(1000, 1000, init = 1:1e6)
我不会重新发明轮子并使用:
big_apply(X, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X[ind, ])
}, a.combine = "rbind", ind = rows_along(X), ncores = nb_cores(), block.size = 500)
明智地选择block.size
(行数)。
如果您想将 R(cpp) 函数应用于 FBM
.
big_apply()
非常有用
编辑: 当然,对于小矩阵,并行性会变慢,因为并行性的开销(通常为 1-3 秒)。查看 X1 和 X2 的结果:
library(bigstatsr)
X1 <- FBM(1000, 1000, init = 1e6)
microbenchmark::microbenchmark(
PAR = big_apply(X1, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X[ind, ])
}, a.combine = "rbind", ind = rows_along(X1), ncores = nb_cores(), block.size = 500),
SEQ = big_apply(X1, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X[ind, ])
}, a.combine = "rbind", ind = rows_along(X1), ncores = 1, block.size = 500),
times = 5
)
Unit: milliseconds
expr min lq mean median uq max neval cld
PAR 1564.20591 1602.0465 1637.77552 1629.9803 1651.04509 1741.59974 5 b
SEQ 68.92936 69.1002 76.70196 72.9173 85.31751 87.24543 5 a
X2 <- FBM(10000, 10000, init = 9999)
microbenchmark::microbenchmark(
PAR = big_apply(X2, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X[ind, ])
}, a.combine = "rbind", ind = rows_along(X2), ncores = nb_cores(), block.size = 500),
SEQ = big_apply(X2, a.FUN = function(X, ind) {
matrixStats::rowQuantiles(X[ind, ])
}, a.combine = "rbind", ind = rows_along(X2), ncores = 1, block.size = 500),
times = 5
)
Unit: seconds
expr min lq mean median uq max neval cld
PAR 4.757409 4.958869 5.071982 5.083381 5.218098 5.342153 5 a
SEQ 10.842828 10.846281 11.177460 11.360162 11.416967 11.421065 5 b
您的矩阵越大,您从并行性中获得的收益就越大。