Rcpp - извлечение строк из списка матриц/dataframes

Question

В качестве продолжения до this question я решил спуститься по маршруту Rcpp по сравнению со сложным синтаксисом в R. Я думаю, что это обеспечит лучшую читаемость (и, возможно, также будет быстрее).

Предположим, у меня есть список data.frame s (который я могу легко преобразовать в матрицы через as). Учитывая предшествующие answe-r-s, это кажется лучшим подходом.

# input data 
my_list <- vector("list", length= 10) 
set.seed(65L) 
for (i in 1:10) { 
    my_list[[i]] <- data.frame(matrix(rnorm(10000),ncol=10)) 
    # alternatively 
    # my_list[[i]] <- matrix(rnorm(10000),ncol=10) 
} 

Каким образом можно извлечь строки из матриц? Цель состоит в том, чтобы создать список с каждым элементом списка, содержащим список из nr-й строки каждого из данных dataframes исходного списка. Я попробовал несколько различных синтаксисом и продолжаю получать ошибки:

#include <Rcpp.h> 
using namespace Rcpp; 
using namespace std: 

List foo(const List& my_list, const int& n_geo) { 
    int n_list = my_list.size(); 
    std::vector<std::vector<double> > list2(n_geo); 

    // needed code.... 

    return wrap(list2); 
} 

варианты

for (int i = 0; i < n_list; i++) { 
    for (int nr = 0; nr < n_geo; nr++) { 
    list2[nr][i] = my_list[i].row(nr); 
    // or list2[nr].push_back(my_list[i].row(nr)); 
    // or list2[nr].push_back(as<double>(my_list[i].row(nr))); 
    // or list2[nr].push_back(as<double>(my_list[i](nr, _))); 
    } 
} 

// or: 
NumericMatrix a = my_list[1] 
... 
NumericMatrix j = my_list[10] 

for (int nr = 0; nr < n_geo; nr++) { 
    list2[nr][1] = // as above 
} 

Ни один из них не работает для меня. Что я делаю не так? Вот ошибки, которые я получаю из моих предыдущих вариантов синтаксиса.

error: no matching function for call to 'as(Rcpp::Matrix<14>::Row)'

или

error: cannot convert 'Rcpp::Matrix<14>::Row {aka Rcpp::MatrixRow<14>}' to 'double' in assignment

Answer 1

Вот один из способов сделать это:

#include <Rcpp.h> 

// x[[nx]][ny,] -> y[[ny]][[nx]] 

// [[Rcpp::export]] 
Rcpp::List Transform(Rcpp::List x) { 
    R_xlen_t nx = x.size(), ny = Rcpp::as<Rcpp::NumericMatrix>(x[0]).nrow(); 
    Rcpp::List y(ny); 

    for (R_xlen_t iy = 0; iy < ny; iy++) { 
     Rcpp::List tmp(nx); 
     for (R_xlen_t ix = 0; ix < nx; ix++) { 
      Rcpp::NumericMatrix mtmp = Rcpp::as<Rcpp::NumericMatrix>(x[ix]); 
      tmp[ix] = mtmp.row(iy); 
     } 
     y[iy] = tmp; 
    } 

    return y; 
} 

/*** R 

L1 <- lapply(1:10, function(x) { 
    matrix(rnorm(20), ncol = 5) 
}) 

L2 <- lapply(1:nrow(L1[[1]]), function(x) { 
    lapply(L1, function(y) unlist(y[x,])) 
}) 

all.equal(L2, Transform(L1)) 
#[1] TRUE 

microbenchmark::microbenchmark(
    "R" = lapply(1:nrow(L1[[1]]), function(x) { 
     lapply(L1, function(y) unlist(y[x,])) 
    }), 
    "Cpp" = Transform(L1), 
    times = 200L) 

#Unit: microseconds 
#expr min  lq  mean median  uq  max neval 
# R 254.660 316.627 383.92739 347.547 392.7705 1909.097 200 
#Cpp 18.314 26.007 71.58795 30.230 38.8650 945.167 200 

*/ 

---

Я не знаю, как это будет масштабироваться; Я думаю, что это просто неэффективная трансформация. Согласно моему комментарию в верхней части источника, кажется, что вы просто делаете своего рода своп координат - ny-я строка nx-го элемента входного списка становится nx-м элементом ny-го элемента выходной список:

x[[nx]][ny,] -> y[[ny]][[nx]] 

---

чтобы устранить ошибки, которые вы получали, Rcpp::List является родовым объектом - технически Rcpp::Vector<VECSXP> - поэтому, когда вы пытаетесь сделать, например,

my_list[i].row(nr) 

компилятор не знает, что my_list[i] является NumericMatrix. Таким образом, вы должны сделать явное приведение с Rcpp::as<>,

Rcpp::NumericMatrix mtmp = Rcpp::as<Rcpp::NumericMatrix>(x[ix]); 
tmp[ix] = mtmp.row(iy); 

---

Я просто использовал matrix элементы в иллюстративных данных, чтобы упростить вещи. На практике вам, вероятно, лучше принудить data.frame s к matrix объектам непосредственно в R, чем пытаться сделать это на C++; это будет намного проще, и, скорее всего, принуждение просто вызывает базовый код C, поэтому на самом деле нет ничего, что можно было бы сделать, пытаясь сделать это иначе.

---

Я хотел бы также отметить, что если вы используете Rcpp::List однородных типов, вы можете выжать немного больше производительности с Rcpp::ListOf<type>.Это позволит вам пропустить приведенные выше преобразования Rcpp::as<type>:

typedef Rcpp::ListOf<Rcpp::NumericMatrix> MatList; 

// [[Rcpp::export]] 
Rcpp::List Transform2(MatList x) { 
    R_xlen_t nx = x.size(), ny = x[0].nrow(); 
    Rcpp::List y(ny); 

    for (R_xlen_t iy = 0; iy < ny; iy++) { 
     Rcpp::List tmp(nx); 
     for (R_xlen_t ix = 0; ix < nx; ix++) { 
      tmp[ix] = x[ix].row(iy); 
     } 
     y[iy] = tmp; 
    } 

    return y; 
} 

/*** R 

L1 <- lapply(1:10, function(x) { 
    matrix(rnorm(20000), ncol = 100) 
}) 

L2 <- lapply(1:nrow(L1[[1]]), function(x) { 
    lapply(L1, function(y) unlist(y[x,])) 
}) 

microbenchmark::microbenchmark(
    "R" = lapply(1:nrow(L1[[1]]), function(x) { 
     lapply(L1, function(y) unlist(y[x,])) 
    }), 
    "Transform" = Transform(L1), 
    "Transform2" = Transform2(L1), 
    times = 200L) 

#Unit: microseconds 
#  expr  min  lq  mean median  uq  max neval 
#   R 6049.594 6318.822 7604.871 6707.242 8592.510 64005.190 200 
# Transform 928.468 1041.936 3130.959 1166.819 1659.745 71552.284 200 
#Transform2 850.912 957.918 1694.329 1061.183 2856.724 4502.065 200 

*/