R cpp вложенные циклы по спискам векторов разных типов: Эффективное решение для шаблона? - PullRequest
Новая
       10

R cpp вложенные циклы по спискам векторов разных типов: Эффективное решение для шаблона?

1 голос
/ 19 апреля 2020

Другая интересная проблема программирования шаблонов R cpp: я хочу вычислить попарно полные наблюдения всех столбцов в R data.frame или матрице. Для матриц код прост благодаря макрокомандам RCPP_RETURN : 

// [[Rcpp::plugins(cpp11)]]
#include <Rcpp.h>
using namespace Rcpp;

template <int RTYPE>
IntegerMatrix pwNobsmCppImpl(const Matrix<RTYPE>& x) {
  int l = x.nrow(), col = x.ncol();
  auto isnnanT = (RTYPE == REALSXP) ? [](typename Rcpp::traits::storage_type<RTYPE>::type x) { return x == x; } : 
       [](typename Rcpp::traits::storage_type<RTYPE>::type x) { return x != Vector<RTYPE>::get_na(); };
  IntegerMatrix out = no_init_matrix(col, col);
  for(int j = 0; j != col; ++j) {
    ConstMatrixColumn<RTYPE> colj = x( _ , j);
    int nj = std::count_if(colj.begin(), colj.end(), isnnanT);
    out(j, j) = nj;
    for(int k = j+1; k != col; ++k) {
      ConstMatrixColumn<RTYPE> colk = x( _ , k);
      int njk = 0;
      for(int i = l; i--; ) if(isnnanT(colj[i]) && isnnanT(colk[i])) ++njk; // fastest? or save logical vector with colj Non-missing?
      out(j, k) = out(k, j) = njk;
    }
  }
  out.attr("dimnames") = List::create(colnames(x), colnames(x));
  return out;
}

template <>
IntegerMatrix pwNobsmCppImpl(const Matrix<CPLXSXP>& x) {
  stop("Not supported SEXP type!");
}

template <>
IntegerMatrix pwNobsmCppImpl(const Matrix<VECSXP>& x) {
  stop("Not supported SEXP type!");
}

template <>
IntegerMatrix pwNobsmCppImpl(const Matrix<RAWSXP>& x) {
  stop("Not supported SEXP type!");
}

template <>
IntegerMatrix pwNobsmCppImpl(const Matrix<EXPRSXP>& x) {
  stop("Not supported SEXP type!");
}

// [[Rcpp::export]]
IntegerMatrix pwNobsmCpp(SEXP x){
  RCPP_RETURN_MATRIX(pwNobsmCppImpl, x);
}

Теперь все сложно, но для data.frames / lists, чьи столбцы могут быть разных типов, очень сложно. Мне кажется, что нативное решение R cpp - это оператор switch(TYPEOF(x[j]), ...), который обрабатывает столбцы списка разного типа в каждом конкретном случае. Однако в этом случае для этого потребуются 2 вложенных переключателя, что приведет к огромному количеству дублирования кода. Мне интересно, можно ли этого избежать с помощью некоторого умного программирования шаблонов. Основной код c для списков: 

// [[Rcpp::export]]
IntegerMatrix pwNobslCpp(const List& x) {
  int l = x.size();
  IntegerMatrix out = no_init_matrix(l, l);
  for(int j = 0; j != l; ++j) {
    int RTYPEj = TYPEOF(x[j]); // The code fails because the SEXPTYPE of x[j] is unknown beforehand!
    auto isnnanTj = (RTYPEj == REALSXP) ? [](typename Rcpp::traits::storage_type<RTYPEj>::type x) { return x == x; } : 
      [](typename Rcpp::traits::storage_type<RTYPEj>::type x) { return x != Vector<RTYPEj>::get_na(); };
    Vector<RTYPEj> colj = x[j];
    int nj = std::count_if(colj.begin(), colj.end(), isnnanTj);
    int rowj = colj.size();
    out(j, j) = nj;
    for(int k = j+1; k != l; ++k) {
      int RTYPEk = TYPEOF(x[k]); // The code fails because the SEXPTYPE of x[k] is unknown beforehand!
      auto isnnanTk = (RTYPEk == REALSXP) ? [](typename Rcpp::traits::storage_type<RTYPEk>::type x) { return x == x; } : 
        [](typename Rcpp::traits::storage_type<RTYPEk>::type x) { return x != Vector<RTYPEk>::get_na(); };
      Vector<RTYPEk> colk = x[k];
      if(colk.size() != rowj) stop("All columns need to have the same length!");
      int njk = 0;
      for(int i = rowj; i--; ) if(isnnanTj(colj[i]) && isnnanTk(colk[i])) ++njk; // fastest? or save logical vector with 
j Non-missing?
      out(j, k) = out(k, j) = njk;
    }
  }
  out.attr("dimnames") = List::create(names(x), names(x));
  return out;
}

Этот код, конечно, не работает, поскольку содержащиеся в нем шаблоны R cpp не активируются с известными значениями. Но шаблонирование всего кода также не имеет смысла из-за парного характера вычислений. Есть ли способ сделать это без вложенного переключателя, которого боятся?

...