передача Matrix :: sparseMatrix в Rcpp - PullRequest
Новая
       93

передача Matrix :: sparseMatrix в Rcpp

1 голос
/ 24 марта 2020

, поэтому я действительно запутался в рекомендуемом способе передачи разреженных матриц из R в c ++. У меня сложилось впечатление, что sp_mat был правильным типом аргумента для этого, как показано в приведенном ниже коде 

testCode = '                                                                                                                                                                                                                                                           
#include <RcppArmadillo.h>                                                                                                                                                                                                                                             
// [[Rcpp::depends(RcppArmadillo)]]                                                                                                                                                                                                                                    

// [[Rcpp::export]]                                                                                                                                                                                                                                                    
void testFun(arma::sp_mat F){                                                                                                                                                                                                                                          

  Rcpp::Rcout << "F has " << F.n_rows << " rows" << std::endl;                                                                                                                                                                                                         

}'

Rcpp::sourceCpp(code = testCode)

n = 70000
M = Matrix::sparseMatrix(i=c(n), j=c(n), x=c(1))

testFun(M)

Однако при запуске этого кода выдается следующая ошибка: 

error: SpMat::init(): requested size is too large
Error in testFun(M) : SpMat::init(): requested size is too large
Calls: testFun -> .Call
Execution halted

Я видел пример в https://gallery.rcpp.org/articles/armadillo-sparse-matrix/, но я не уверен, говорит ли он, что каждый раз, когда мы передаем разреженную матрицу в c ++, мы должны использовать предоставленную там функцию? Спасибо за разъяснение!

Ответы [ 2 ]

2 голосов
/ 25 марта 2020

Хорошо, я думаю, что найду ответ. В основном броненосец выдает эту ошибку, если общее количество элементов больше, чем размер переменной, в которой хранится количество элементов, как показано здесь: https://gitlab.com/conradsnicta/armadillo-code/-/blob/9.900.x/include/armadillo_bits/SpMat_meat.hpp

Кто-то понял это раньше и предоставил решение здесь: Большие матрицы в RcppArmadillo через ARMA_64BIT_WORD определяют

1 голос
/ 25 марта 2020

Если вы повторно посмотрите пример basi c в R cpp Gallery и настроите один или два разреженных объекта Matrix, станет очевидно, что высокое значение для j приводит к полному расширение в слоте p (проверьте объект, возвращенный из sparseMatrix).

Итак, вот более простой пример с (довольно высоким) значением i, но низким j. Я думаю, что вы должны быть в состоянии взять это отсюда:

Код 

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp ;

// [[Rcpp::export]]
void convertSparse(S4 mat) {

    // obtain dim, i, p. x from S4 object
    IntegerVector dims = mat.slot("Dim");
    arma::urowvec i = Rcpp::as<arma::urowvec>(mat.slot("i"));
    arma::urowvec p = Rcpp::as<arma::urowvec>(mat.slot("p"));
    arma::vec x     = Rcpp::as<arma::vec>(mat.slot("x"));

    int nrow = dims[0], ncol = dims[1];

    // use Armadillo sparse matrix constructor
    arma::sp_mat res(i, p, x, nrow, ncol);
    Rcout << "SpMat res:\n" << res << std::endl;
}

// [[Rcpp::export]]
void convertSparse2(S4 mat) {         // slight improvement with two non-nested loops

    IntegerVector dims = mat.slot("Dim");
    arma::urowvec i = Rcpp::as<arma::urowvec>(mat.slot("i"));
    arma::urowvec p = Rcpp::as<arma::urowvec>(mat.slot("p"));
    arma::vec x     = Rcpp::as<arma::vec>(mat.slot("x"));

    int nrow = dims[0], ncol = dims[1];
    arma::sp_mat res(nrow, ncol);

    // create space for values, and copy
    arma::access::rw(res.values) = arma::memory::acquire_chunked<double>(x.size() + 1);
    arma::arrayops::copy(arma::access::rwp(res.values), x.begin(), x.size() + 1);

    // create space for row_indices, and copy
    arma::access::rw(res.row_indices) = arma::memory::acquire_chunked<arma::uword>(i.size() + 1);
    arma::arrayops::copy(arma::access::rwp(res.row_indices), i.begin(), i.size() + 1);

    // create space for col_ptrs, and copy
    arma::access::rw(res.col_ptrs) = arma::memory::acquire<arma::uword>(p.size() + 2);
    arma::arrayops::copy(arma::access::rwp(res.col_ptrs), p.begin(), p.size() + 1);

    // important: set the sentinel as well
    arma::access::rwp(res.col_ptrs)[p.size()+1] = std::numeric_limits<arma::uword>::max();

    // set the number of non-zero elements
    arma::access::rw(res.n_nonzero) = x.size();

    Rcout << "SpMat res:\n" << res << std::endl;
}


/*** R
suppressMessages({
  library(methods)
  library(Matrix)
})
i <- c(1,3:6)
j <- c(2,9,6:8)
x <- 5 * (1:5)
A <- sparseMatrix(i, j, x = x)
print(A)
convertSparse(A)

i <- 56789
j <- 87
x <- 42
B <- sparseMatrix(i, j, x=x)
#print(B)
convertSparse(B)
convertSparse2(B)
*/

Выход 

R> Rcpp::sourceCpp("~/git/stackoverflow/60838958/answer.cpp")

R> suppressMessages({library(methods); library(Matrix)})

R> i <- c(1,3:6)

R> j <- c(2,9,6:8)

R> x <- 5 * (1:5)

R> A <- sparseMatrix(i, j, x = x)

R> print(A)
6 x 9 sparse Matrix of class "dgCMatrix"

[1,] . 5 . . .  .  .  .  .
[2,] . . . . .  .  .  .  .
[3,] . . . . .  .  .  . 10
[4,] . . . . . 15  .  .  .
[5,] . . . . .  . 20  .  .
[6,] . . . . .  .  . 25  .

R> convertSparse(A)
SpMat res:
[matrix size: 6x9; n_nonzero: 5; density: 9.26%]

     (0, 1)          5.0000
     (3, 5)         15.0000
     (4, 6)         20.0000
     (5, 7)         25.0000
     (2, 8)         10.0000



R> i <- 56789

R> j <- 87

R> x <- 42

R> B <- sparseMatrix(i, j, x=x)

R> #print(B)
R> convertSparse(B)
SpMat res:
[matrix size: 56789x87; n_nonzero: 1; density: 2.02e-05%]

 (56788, 86)        42.0000



R> convertSparse2(B)
SpMat res:
[matrix size: 56789x87; n_nonzero: 1; density: 2.02e-05%]

 (56788, 86)        42.0000


R>

Редактировать Действительно, хорошее напоминание. Если мы добавим 

#define ARMA_64BIT_WORD 1

перед включением заголовка RcppArmadillo.h, то с большими значениями i и j все будет хорошо. Хвост вывода ниже.

Обновлен вывод (только хвост) 

R> i <- 56789

R> j <- 87654

R> x <- 42

R> B <- sparseMatrix(i, j, x=x)

R> #print(B)
R> convertSparse(B)
SpMat res:
[matrix size: 56789x87654; n_nonzero: 1; density: 2.01e-08%]

 (56788, 87653)     42.0000



R> convertSparse2(B)
SpMat res:
[matrix size: 56789x87654; n_nonzero: 1; density: 2.01e-08%]

 (56788, 87653)     42.0000


R>
...