r - 使用 RcppParallel 分组求和

Question

我编写了一个函数，可以对组中的值进行求和。它采用两个长度相同的向量:v 和 g，并且应返回一个与 g 中的唯一元素长度相同的向量。组被编码为从零开始的整数。使用 Rcpp::sourceCpp 代码可以编译，但从 R 调用时(例如 sg(runif(6),rep(0:1,each = 3)))会返回数字(0)。

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

struct SumsInGroups: public Worker
{
  const RVector<double> v;
  const RVector<int> g;

  RVector<double> s;

  SumsInGroups(const NumericVector v, const IntegerVector g, NumericVector s): v(v), g(g), s(s) {}
  SumsInGroups(const SumsInGroups& p, Split): v(p.v), g(p.g), s(p.s) {}

  void operator()(std::size_t begin, std::size_t end) {
    for (std::size_t i = begin; i < end; ++i) {
      if (s[g[i]] != s[g[i]]) s[g[i]] = v[i];
      else s[g[i]] += v[i];
    }
  }

  void join(const SumsInGroups& rhs) {
    for(std::size_t i = 0; i < s.length(); i++) {
      s[i] += rhs.s[i];
    }
  }
};

// [[Rcpp::export]]
RVector<double> sg(NumericVector v, IntegerVector g) {
  NumericVector s;
  SumsInGroups p(v, g, s);
  parallelReduce(0, v.length(), p);
  return p.s;
}

我对 RcppParallel 非常陌生，因此欢迎任何意见和建议。

Answer 1

您需要初始化s。我建议用零初始化。这是对我有用的代码。请注意，由于我用零进行初始化，因此不需要您在运算符 () 中进行检查。

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

struct SumsInGroups: public Worker
{
    const RVector<double> v;
    const RVector<int> g;

    RVector<double> s;

    SumsInGroups(const NumericVector v, const IntegerVector g, NumericVector s): v(v), g(g), s(s) {}
    SumsInGroups(const SumsInGroups& p, Split): v(p.v), g(p.g), s(p.s) {}

    void operator()(std::size_t begin, std::size_t end) {
        for (std::size_t i = begin; i < end; ++i) {
            s[g[i]] += v[i];
        }
    }

    void join(const SumsInGroups& rhs) {
        for(std::size_t i = 0; i < s.length(); i++) {
            s[i] += rhs.s[i];
        }
    }
};

// [[Rcpp::export]]
RVector<double> sg(NumericVector v, IntegerVector g) {
    NumericVector s(*std::max_element(g.begin(), g.end()) + 1);
    SumsInGroups p(v, g, s);
    parallelReduce(0, v.length(), p);
    return p.s;
}

/*** R
set.seed(101)
o <- runif(15)
i <-sample(0:3,15, rep = TRUE)
sg(o, i)
tapply(o, i, sum)
*/