Как да сумирам 10 стъпкови реда от матрица в Rcpp?

Искам да получа следните резултати с помощта на Rcpp. Когато има големи данни, R е бавен. Затова се опитах да кодирам в Rcpp.

x <- matrix(1:150, ncol = 5)
z <- matrix(nrow = nrow(x) / 10, ncol = 5)
for (i in 1:5) {
    for (j in 1:(nrow(x) / 10)) {
    k = (j - 1) * 10 + 1;
    z[j, i] <- sum(x[k:(k+9), i])
    }
}
x
       [,1] [,2] [,3] [,4] [,5]
 [1,]    1   31   61   91  121
 [2,]    2   32   62   92  122
 [3,]    3   33   63   93  123
 [4,]    4   34   64   94  124
 [5,]    5   35   65   95  125
 [6,]    6   36   66   96  126
 [7,]    7   37   67   97  127
 [8,]    8   38   68   98  128
 [9,]    9   39   69   99  129
[10,]   10   40   70  100  130
[11,]   11   41   71  101  131
[12,]   12   42   72  102  132
[13,]   13   43   73  103  133
[14,]   14   44   74  104  134
[15,]   15   45   75  105  135
[16,]   16   46   76  106  136
[17,]   17   47   77  107  137
[18,]   18   48   78  108  138
[19,]   19   49   79  109  139
[20,]   20   50   80  110  140
[21,]   21   51   81  111  141
[22,]   22   52   82  112  142
[23,]   23   53   83  113  143
[24,]   24   54   84  114  144
[25,]   25   55   85  115  145
[26,]   26   56   86  116  146
[27,]   27   57   87  117  147
[28,]   28   58   88  118  148
[29,]   29   59   89  119  149
[30,]   30   60   90  120  150

z
      [,1] [,2] [,3] [,4] [,5]
 [1,]   55  355  655  955 1255
 [2,]  155  455  755 1055 1355
 [3,]  255  555  855 1155 1455

Rcpp на кода, който опитах, е както следва.

#include <Rcpp.h> 
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector mySum(NumericMatrix x) {

    int ncol = x.ncol();
    int nrow = x.nrow();
    int outRow = nrow / 10;
    int i;
    int j;
    int k;
    Rcpp::NumericMatrix z(outRow, ncol);

    for (i = 0; i < ncol; i++) {
        for (j = 0; j < outRow; j++) {
        k = j * 10;
        Rcpp::SubMatrix<REALSXP> sm = x(Range(k, k + 9), i);
        Rcpp::NumericMatrix m(sm);
        double s = Rcpp::sum(m);
        z(j, i) = s;
        }
    }
  return z;
}

Въпреки това не се движи поради грешка. Моля, кажете ми решението.

test.cpp: In function 'Rcpp::NumericVector mySum(Rcpp::NumericMatrix)':
test.cpp:18:59: error: no match for call to '(Rcpp::NumericMatrix {aka Rcpp::Matrix<14>}) (Rcpp::Range, int&)'

person zunda    schedule 01.04.2015    source източник
comment
Какъв е размерът на вашите данни?   -  person David Arenburg    schedule 01.04.2015
comment
Здравейте. Благодаря за отговора. Наскоро трябваше да изчисля матрицата от 62 400 реда * 4 100 колони.   -  person zunda    schedule 01.04.2015


Отговори (2)


Предпочитам да използвам RcppArmadillo, когато правя нещо с матрици, една от причините е, че документацията е толкова добра (http://arma.sourceforge.net/docs.html#accu). Пренаписах леко вашия код и изглежда работи добре:

library(RcppArmadillo)
library(Rcpp)

cppFunction("
NumericMatrix mySum(arma::mat x) {

    int ncol = x.n_cols;
    int nrow = x.n_rows;
    int outRow = nrow / 10;
    int i, j, k;
    NumericMatrix z(outRow, ncol);

    for (i = 0; i < ncol; i++) {
        for (j = 0; j < outRow; j++) {
            k = j * 10;
            arma::mat sm = x(arma::span(k, k+9), i);
            z(j, i) = arma::accu(sm);
        }
    }
    return z;
}
", depends = "RcppArmadillo")

x <- matrix(1:150, ncol = 5)
mySum(x)
     [,1] [,2] [,3] [,4] [,5]
[1,]   55  355  655  955 1255
[2,]  155  455  755 1055 1355
[3,]  255  555  855 1155 1455
person Jeff    schedule 01.04.2015
comment
Благодаря за отговора. Вашият код ще работи перфектно. Искам да науча за RcppArmadillo. - person zunda; 01.04.2015

Всъщност има напълно векторизирана функция в основата R, наречена rowsum, която може да сумира по групи много ефективно (Като странична бележка, R не винаги е бавен, зависи главно от това как го използвате).

x <- matrix(1:150, ncol = 5)
rowsum.default(x, cumsum(seq_len(nrow(x)) %% 10L == 1L), reorder = FALSE)
#   [,1] [,2] [,3] [,4] [,5]
# 1   55  355  655  955 1255
# 2  155  455  755 1055 1355
# 3  255  555  855 1155 1455

Със сигурност ще бъде по-бавно от версията Rcpp, но на моята система матрица от 20 мм редове с 5 колони се изпълнява за по-малко от 3 секунди

x <- matrix(seq_len(1e8), ncol = 5)
dim(x)
## [1] 20000000        5
system.time(mySum(x))
# user  system elapsed 
# 0.72    0.24    0.96 
system.time(rowsum.default(x, cumsum(seq_len(nrow(x)) %% 10L == 1L), reorder = FALSE))
# user  system elapsed 
# 2.77    0.15    2.93 

Редактиране: според вашия коментар, докато тествате върху вашия реален набор от данни rowsum работи дори по-бързо от версията Rcpp

x <- matrix(seq_len(62400*4100), ncol = 4100)
dim(x)
## [1] 62400  4100
system.time(mySum(x))
# user  system elapsed 
# 1.53    1.03    2.57 
system.time(rowsum.default(x, cumsum(seq_len(nrow(x)) %% 10L == 1L), reorder = FALSE))
# user  system elapsed 
# 1.48    0.00    1.50 
person David Arenburg    schedule 01.04.2015
comment
Благодаря за отговора. Вашият код е интересен и полезен. Както казвате, R не е бавен в зависимост от кодирането. - person zunda; 01.04.2015
comment
какво ще кажете за system.time(data.table(x)[, lapply(.SD, sum), by=rep(1:(nrow(x)/10), each=10)]) - person Khashaa; 01.04.2015
comment
@Khashaa ще бъде много по-бавно за много колони, но предполагам, че е бързо, когато няма толкова много колони. - person David Arenburg; 01.04.2015
comment
Редактирах отговора си, така че z е NumericVector. Мисля, че причината, поради която работи по-бавно с големи набори от данни, се дължи на скритото преобразуване между NumericVector и arma::mat. Сега трябва да е по-бързо, опитах се да стартирам вашия пример, но старият ми лаптоп не го оцени. - person Jeff; 01.04.2015
comment
@jeff Не можах да компилирам новата ти версия. - person David Arenburg; 01.04.2015
comment
@DavidArenburg съжалявам, че реших да напиша Vector вместо Matrix по някаква причина.... сега трябва да работи. - person Jeff; 01.04.2015
comment
@jeff Добре, просто го стартирах отново и резултатите останаха абсолютно същите. - person David Arenburg; 01.04.2015
comment
@Джеф, това е добре. Вашето решение все още е страхотно (гласувано) и изненадващо много подценено. - person David Arenburg; 01.04.2015