base-R seq 的 Rcpp 版本丢弃值

Question

我写了一个 Rcpp 版本的 base-R seq 函数。

library(Rcpp)

cppFunction('NumericVector seqC(double x, double y, double by) {

  // length of result vector
  int nRatio = (y - x) / by;
  NumericVector anOut(nRatio + 1);

  // compute sequence
  int n = 0;
  for (double i = x; i <= y; i = i + by) {
    anOut[n] = i;
    n += 1;
  }

  return anOut;
}')

对于下面的测试，它工作得很好。

seqC(1, 11, 2)
[1]  1  3  5  7  9 11

seqC(1, 10, 2)
[1]  1  3  5  7  9 11

此外，它（有时）在传递带有十进制数字的值而不是整数。

seqC(0.43, 0.45, 0.001)
[1] 0.430 0.431 0.432 0.433 0.434 0.435 0.436 0.437 0.438 0.439 0.440 0.441 0.442 0.443 0.444 0.445 0.446 0.447 0.448 0.449 0.450

但是，自上次以来，该功能有时似乎无法按预期工作序列的条目被删除（或者更确切地说，输出向量 anOut 没有合适的大小），这 - 根据我相当缺乏的 C++ 技能，可能归因于某种舍入错误。

seqC(0.53, 0.59, 0.001)
 [1] 0.530 0.531 0.532 0.533 0.534 0.535 0.536 0.537 0.538 0.539 0.540 0.541 0.542 0.543 0.544 0.545 0.546 0.547 0.548 0.549 0.550 0.551
[23] 0.552 0.553 0.554 0.555 0.556 0.557 0.558 0.559 0.560 0.561 0.562 0.563 0.564 0.565 0.566 0.567 0.568 0.569 0.570 0.571 0.572 0.573
[45] 0.574 0.575 0.576 0.577 0.578 0.579 0.580 0.581 0.582 0.583 0.584 0.585 0.586 0.587 0.588 0.589

例如，在最后一个示例中，缺少最后一个值 (0.590)。做有人知道如何解决这个问题吗？

Answer 1

“<=”可能会给浮点数带来困难。这是著名问题 "Why are these numbers not equal?" 的变体。此外，矢量长度也有类似的问题，在你的最后一个例子中应该是 60，但它实际上被计算为 59。这很可能是由于转换为整数（通过转换，即，截断）像 59.999999 或类似的值。

这些问题似乎很难修复，所以我重写了相当一部分代码，希望现在功能能按要求运行。

下面的代码应该为基本上任何类型的递增级数（即 y > x、by > 0）提供正确的结果。

cppFunction('NumericVector seqC(double x, double y, double by) {
 NumericVector anOut(1);
 // compute sequence
 double min_by = 1.e-8;
 if (by < min_by) min_by = by/100;
 double i = x + by;
 anOut(0) = x;
 while(i/min_by < y/min_by + 1) { 
  anOut.push_back(i);
  i += by;
 }
return anOut;
}')

希望这对您有所帮助。非常感谢@Konrad Rudolph 指出我之前尝试中的错误！

Answer 2

如其他人所述，您遇到的问题基本上是浮点运算错误。一个常见的解决方法是将 doubles 缩放到足够大的整数，执行任务，然后将结果调整为输入的原始比例。我采用了与@RHertel 略有不同的方法，让缩放量 (adjust) 由增量的精度决定，而不是使用固定量，但想法本质上是相同的。

#include <Rcpp.h>

struct add_multiple {
  int incr;
  int count;
  add_multiple(int incr)
    : incr(incr), count(0)
    {}
  inline int operator()(int d) {
    return d + incr * count++;
  }
};

// [[Rcpp::export]]
Rcpp::NumericVector rcpp_seq(double from_, double to_, double by_ = 1.0) {
  int adjust = std::pow(10, std::ceil(std::log10(10 / by_)) - 1);
  int from = adjust * from_;
  int to = adjust * to_;
  int by = adjust * by_;

  std::size_t n = ((to - from) / by) + 1;
  Rcpp::IntegerVector res = Rcpp::rep(from, n);
  add_multiple ftor(by);

  std::transform(res.begin(), res.end(), res.begin(), ftor);
  return Rcpp::NumericVector(res) / adjust;
}

/*** R

all.equal(seq(.53, .59, .001), seqC(.53, .59, .001)) &&
  all.equal(seq(.53, .59, .001), rcpp_seq(.53, .59, .001))
# [1] TRUE

all.equal(seq(.53, .54, .000001), seqC(.53, .54, .000001)) &&
  all.equal(seq(.53, .54, .000001), rcpp_seq(.53, .54, .000001))
# [1] TRUE 

microbenchmark::microbenchmark(
  "seq" = seq(.53, .54, .000001),
  "seqC" = seqC(0.53, 0.54, 0.000001),
  "rcpp_seq" = rcpp_seq(0.53, 0.54, 0.000001),
  times = 100L)
Unit: microseconds
     expr        min          lq        mean     median         uq        max neval
      seq    896.190   1015.7940   1167.4708   1132.466   1221.624   1651.571   100
     seqC 212293.307 219527.6590 226933.4329 223384.592 227860.410 398462.561   100
 rcpp_seq    182.848    194.1665    225.4338    227.396    244.942    320.436   100

*/

其中 seqC 是@RHertel 修改后的实现，产生了正确的结果。 FWIW 我认为此函数的性能较慢主要是因为 push_back 在 NumericVector 类型上的使用，Rcpp 开发人员强烈建议不要这样做。

base-R seq 的 Rcpp 版本丢弃值

Rcpp version of base-R seq drops values

c++

r

rcpp