Работа со значениями NA с использованием Rcpp

Question

Я тестирую часть своего кода, который показан ниже:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix testOutMat(const int& ncols, const int& nrows, const NumericVector& col_prob){

  //Store row and column positions
  NumericVector col_pos = no_init(nrows);
  NumericVector row_pos = no_init(nrows);

  int row_val;
  int nz_counter=0;
  for(int j=0; j<ncols; ++j){ 
    for(int i=0; i<nrows; ++i){
      row_val = R::rbinom(1,col_prob[j]);
      Rcout << "i,j: " << i << "," << j << std::endl;
      Rcout << "val: " << row_val << std::endl;
      if(row_val==1){ //if (i,j)th entry is a 1, save location
        row_pos[i] = i;
        col_pos[i] = j;
        nz_counter += 1;
      } else{ //assign as NA
        row_pos[i] = NA_REAL;
        col_pos[i] = NA_REAL;
      }
      Rcout << "row_pos[i]: " << row_pos[i] << std::endl;
      Rcout << "col_pos[i]: " << col_pos[i] << std::endl;
      Rcout << "num non-zeros: " << nz_counter << std::endl;
    }
  }

  NumericMatrix out = no_init(nz_counter,2);

  Rcout << "Printing output matrix" << std::endl;
  for(int i=0; i<nz_counter; ++i){
    if(!Rcpp::NumericVector::is_na(row_pos[i])){ 
      out(i,0) = row_pos[i];
      out(i,1) = col_pos[i];
    }
    Rcout << "row_pos[i]: " << row_pos[i] << std::endl;
    Rcout << "col_pos[i]: " << col_pos[i] << std::endl; 
  }

  return out;
}

/*** R
set.seed(1)
res <- testOutMat(ncols=5,nrows=5,col_prob = runif(20, 0.1, 0.2))
*/

Исходя из вывода, у меня есть записи (i,j)={(0,0),(3,1)}, отличные от нуля, так что res должна быть 2x2 матрица с 0 0 в первой строке и 3 1 во второй. Однако я получаю что-то совсем другое:

     [,1] [,2]
[1,]   64 1024
[2,]    1    4

Я подозреваю, что это связано с тем, как я обращаюсь NA с. Общая цель функции состоит в том, чтобы сгенерировать индексы строк и столбцов для ненулевых элементов (генерируемых при вызове rbinom).

Я пытался отладить это в течение некоторого времени, и, похоже, не могу получить исправление.

Answer 1

Проблема здесь в том, что вы пишете более row_pos и col_pos снова и снова (ncols раз) без какого-либо отслеживания предыдущего результата. Это, в сочетании с вашим no_init() использованием, является причиной конечного результата, который вы видите. Мы можем немного изменить ваш код, чтобы row_pos и col_pos не были перезаписаны:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerMatrix testOutMat(const int ncols, const int nrows,
                         const NumericVector& col_prob) {

    IntegerMatrix binomial_deviates(nrows, ncols);
    IntegerVector row_positions;
    IntegerVector col_positions;
    int nz_counter = 0;

    for ( int j = 0; j < ncols; ++j ) {
        binomial_deviates(_, j) = rbinom(nrows, 1, col_prob[j]);
        for ( int i = 0; i < nrows; ++i ) {
            if ( binomial_deviates(i, j) == 1 ) {
                row_positions.push_back(i);
                col_positions.push_back(j);
                nz_counter += 1;
            }
        }
    }

    IntegerMatrix out(nz_counter, 2);

    for ( int i = 0; i < nz_counter; ++i ) {
        out(i, 0) = row_positions[i];
        out(i, 1) = col_positions[i];
    }

    return out;
}

/*** R
set.seed(1)
res <- testOutMat(ncols=5,nrows=5,col_prob = runif(20, 0.1, 0.2))
*/

Результат:

> set.seed(1)

> res <- testOutMat(ncols=5,nrows=5,col_prob = runif(20, 0.1, 0.2))
> res
     [,1] [,2]
[1,]    0    0
[2,]    3    1