Извлечь каждый k-й столбец матрицы arma :: mat в rcpp - PullRequest
Новая
       50

Извлечь каждый k-й столбец матрицы arma :: mat в rcpp

0 голосов
/ 17 октября 2018

Я боролся с поднабором столбцов матрицы класса arma :: mat.

Допустим, дано arma::mat X, и я попытался создать вектор индексов IDX, чтобы сделать X.cols(IDX).В частности, индексный вектор имеет каждое k-ое целое число от 1 до p (размерность X).Например, каждый может быть заинтересован во всех четных столбцах IDX=[2,4,6,8, ...].

На основании этой документации , смежных индексов, таких как [0, 1, 2, ..., m-1] может быть легко извлечен с помощью X.cols(0, m - 1), если m <= p.Однако я не смог найти хороший способ для подстановки матрицы с вектором индекса <code>IDX, описанным выше.

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

Мой файл "subset_armamat.cpp" выглядит как 

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
mat subset_armamat(mat X, int k){
  uvec IDX = "every k-th integer from 0 to X.ncols";
  return X.cols(IDX);
}

, а R-код для выполнения определенной функции - 

library("Rcpp")
sourceCpp("subset_armamat.cpp")
subset_armamat(matrix(1:10, 2, 5, byrow = T), 2)

Ожидается, что это приведет к 2-Матрица на-3, так как следующий код R даст 

> matrix(1:10, 2, 5, byrow = T)[,seq(1, 5, by = 2)]
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    6    8   10

Было бы очень полезно, если бы вы дали какой-либо вклад.

ps Я пытался

  • генерировать вектор последовательности seq(1,m) * 2 вручную, но это не работает с X.cols().
  • или найдите индекс, используя find(seq(1,p) % 2 == 0), но оператор % не работает должным образом между seq(1,p) и 2.

Ответы [ 2 ]

0 голосов
/ 17 октября 2018

F.Ответ Приве показал, что вы можете использовать uvec для подмножества матрицы, используя .cols(), даже если она не является непрерывным диапазоном, используя базовую функцию R seq() для генерации последовательности.Далее я покажу, что вы можете сгенерировать последовательность, используя функцию Armadillo;Вы можете использовать arma::regspace() - он «генерирует [s] вектор с регулярно расположенными элементами» ( Источник документации Armadillo ): 

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
mat subset_armamat(mat X, int k) {
    uvec IDX = regspace<uvec>(0,  k,  X.n_cols-1);
    return X.cols(IDX);
}

Для сравнения с вызовом R seq() (где subset_armamatR() - функция из ответа Ф. Приве): 

library("Rcpp")
sourceCpp("subset_armamat.cpp")
mat <- matrix(1:10, 2, 5, byrow = TRUE)
subset_armamat(mat, 2)
#>      [,1] [,2] [,3]
#> [1,]    1    3    5
#> [2,]    6    8   10
subset_armamatR(mat, 2)
#>      [,1] [,2] [,3]
#> [1,]    1    3    5
#> [2,]    6    8   10
library(microbenchmark)
microbenchmark(Rseq = subset_armamatR(mat, 2),
               regspace = subset_armamat(mat, 2))
#> Unit: microseconds
#>      expr     min       lq     mean   median       uq       max neval cld
#>      Rseq 235.535 239.1615 291.1954 241.9850 248.6005  4704.467   100   a
#>  regspace  14.221  15.0225 520.9235  15.8165  16.6740 50408.375   100   a

Обновление: передача по ссылке

Комментарий от hbrerkere требует некоторых кратких дополнительныхобсуждение.Если вы вызываете эту функцию из C ++, вы наберете скорость, изменив mat subset_armamat(mat X, int k) на mat subset_armamat(const mat& X, int k).Подобная передача по ссылке позволяет избежать ненужной копии, и когда вы не собираетесь изменять объект, переданный по ссылке, вы должны использовать const.Однако, если вы вызываете эту функцию из R, вы не можете избежать копирования, так как arma::mat не является собственным типом R (см., Например, этот ответ от Dirk Eddelbuettel (сопровождающий Rcpp и RcppArmadillo). Рассмотрим следующий пример: 

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

// [[Rcpp::export]]
void reference_example(arma::mat& X) {
    X(0, 0) = 42;
}

// [[Rcpp::export]]
void print_reference_example(arma::mat X) {
    reference_example(X);
    Rcpp::Rcout << X << "\n";
}

Затем вызов из R: 

library("Rcpp")
sourceCpp("reference_example.cpp")
mat <- matrix(1:4, 2, 2)
mat
#>      [,1] [,2]
#> [1,]    1    3
#> [2,]    2    4
reference_example(mat)
mat
#>      [,1] [,2]
#> [1,]    1    3
#> [2,]    2    4
print_reference_example(mat)
#>    42.0000    3.0000
#>     2.0000    4.0000
mat
#>      [,1] [,2]
#> [1,]    1    3
#> [2,]    2    4
0 голосов
/ 17 октября 2018

Это работает: 

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
mat subset_armamat(mat X, int k) {

  // Obtain environment containing function
  Rcpp::Environment base("package:base"); 

  // Make function callable from C++
  Rcpp::Function seq = base["seq"];    

  uvec IDX = as<uvec>(seq(0, X.n_cols, k));
  return X.cols(IDX);
}

Я просто вызываю функцию R base::seq() из Rcpp.

...