Как я могу разместить в Rcpp строки или столбцы bigstatsr::FBM и сохранить их в векторе?

Question

У меня есть функция, которая вычисляет основную сводную статистику из строк (или столбцов) данной матрицы, и сейчас я пытаюсь также использовать эту функцию с 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 делает задачу в 5-10 раз медленнее на моем ПК (4 процессора и использование окон, к сожалению). использование большей матрицы X2 и оставление block.size по умолчанию занимает в 10 раз больше времени с 4 ядрами вместо непараллельной версии.

Результат для 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

Answer 1

Если у вас есть

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 (количество рядов) с умом. функция big_apply() очень полезно, если вы хотите применить функцию R(cpp) к блокам FBM,

Изменить: Конечно, параллелизм будет медленнее для небольших матриц, из-за перевеса параллелизма (обычно, 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

Чем больше ваша матрица, тем больше вы выиграете от параллелизма.