Странное поведение при инкрементальной выборке с использованием RcppArmadillo :: sample

Question

Я пытаюсь реализовать некоторые отрисовки, используя схему polya urn, используя Rcpp. По сути, у меня есть матрица, из которой я рисую, и вторая матрица с весами, пропорциональными вероятностям. После каждого розыгрыша мне нужно увеличивать вес любой ячейки, которую я нарисовал.

Я столкнулся с некоторыми ошибками индексации, которые побудили меня рассмотреть выборку более широко, и я обнаружил, что моя матрица весов модифицировалась с помощью RcppArmadillo :: sample. Два вопроса (1) это поведение, которое я должен был ожидать, или это ошибка, о которой я должен сообщить где-нибудь? (2) Есть идеи о текущем обходном пути? Вот воспроизводимый пример:

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

using namespace Rcpp ;

// [[Rcpp::export]]
void sampler(int N, int inc, NumericMatrix& weight_matrix, int reps) {

  IntegerVector wm_tmp = seq_along(weight_matrix);
  Rcout << "Initial weight_matrix:\n" << weight_matrix << "\n";

  int x_ind;
  for(int i = 0; i < reps; ++i) {
    x_ind = RcppArmadillo::sample(wm_tmp, 1, true, weight_matrix)(0) - 1;
    Rcout << "Weight matrix after sample: (rep = " << i << ")\n" << weight_matrix << "\n";
    Rcout << "x_ind: " << x_ind  << "\n"; 
    // get indices
    weight_matrix[x_ind] = weight_matrix[x_ind] +  inc;
    Rcout << "Add increment of " << inc << " to weight_matrix:\n" << weight_matrix << "\n";
  }
}
// 
// // [[Rcpp::export]]
// IntegerVector seq_cpp(IntegerMatrix x) {
//   IntegerVector tmp = seq_along(x);
//   IntegerVector ret = RcppArmadillo::sample(tmp, 2, true);
//   return ret;
// }

/*** R
weight_matrix <- matrix(1, 5, 2)
sampler(5, 1, weight_matrix, 3)

weight_matrix <- matrix(1, 5, 2)
sampler(5, 0, weight_matrix, 3)
*/

Спасибо!

Answer 1

Это известное и задокументированное поведение.

Вы могли бы сделать

i) Используйте Rcpp::clone(), чтобы создать отдельную копию вашего SEXP (т.е. NumericMatrix).

ii) Используйте вместо этого матрицу Armadillo и передайте как const arma::mat & m.

Существуют архитектурные причины, связанные с тем, как R организует свою структуру данных, что означает, что мы не можем предоставить вам быстрый доступ (без копий!), А также защитить от записи.