Rcpp 移动平均 - 边界错误导致致命错误

Rcpp moving average - boundary error leads to fatal error

我使用滚动加权移动平均函数,其代码如下所示。它是通过 Rcpp 在 C++ 中编码的。 此功能适用于大多数时间序列,没有循环问题或类似问题。我在下面提供了一个有时会触发致命错误的长度为 2 的时间序列。 我找不到错误的原因。

感谢您的帮助! =)

这是 R 代码:

# Install packages
sourceCpp("partialMA.cpp")
spencer_weights=c( -3, -6, -5, 3, 21, 46, 67, 0, 67, 46, 21, 3, -5, -6, -3)
spencer_ma <- function(x) roll_mean(x,spencer_weights)

x=c(11.026420323685528,0.25933761651337001)
spencer_ma(x) # works
for(i in 1:1000) spencer_ma(x) # triggers the fatal error 

我在下面包含了我的 roll_mean 函数的 C++ 代码:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector roll_mean(const NumericVector& x,
                        const NumericVector& w) {

  int n = x.size();
  int w_size = w.size();
  int size = (w_size - 1) / 2;

  NumericVector res(n);
  int i, ind_x, ind_w;

  double w_sum = Rcpp::sum(w), tmp_wsum, tmp_xwsum, tmp_w;

  // beginning
  for (i = 0; i < size; i++) {
    tmp_xwsum = tmp_wsum = 0;
    for (ind_x = i + size, ind_w = w_size - 1; ind_x >= 0; ind_x--, ind_w--) {
      tmp_w = w[ind_w];
      tmp_wsum += tmp_w;
      tmp_xwsum += x[ind_x] * tmp_w;
    }
    res[i] = tmp_xwsum / tmp_wsum;
  }

  // middle
  int lim2 = n - size;
  for (; i < lim2; i++) {
    tmp_xwsum = 0;
    for (ind_x = i - size, ind_w = 0; ind_w < w_size; ind_x++, ind_w++) {
      tmp_xwsum += x[ind_x] * w[ind_w];
    }
    res[i] = tmp_xwsum / w_sum;
  }

  // end
  for (; i < n; i++) {
    tmp_xwsum = tmp_wsum = 0;
    for (ind_x = i - size, ind_w = 0; ind_x < n; ind_x++, ind_w++) {
      tmp_w = w[ind_w];
      tmp_wsum += tmp_w;
      tmp_xwsum += x[ind_x] * tmp_w;
    }
    res[i] = tmp_xwsum / tmp_wsum;
  }

  return res;
}

A Wild Index Out of Bounds Error Appeared!

您可以通过将元素访问器从 [] 切换到 () 来查明问题。后者有一个内置的边界检查,例如index 介于 0n-1 之间。

运行 带有 built-in 检查的代码给出:

 Error in roll_mean(x, spencer_weights) : 
  Index out of bounds: [index=7; extent=2]. 

因此,所使用的索引大大超过了向量的长度。添加跟踪语句表示它的第一个循环是错误的。

#include <Rcpp.h>
// [[Rcpp::export]]
NumericVector roll_mean(const NumericVector& x,
                        const NumericVector& w) {

  int n = x.size();
  int w_size = w.size();
  int size = (w_size - 1) / 2;

  Rcpp::Rcout << n << ", w_size: " << w_size << ", size: " << size << std::endl;

  NumericVector res(n);
  int i, ind_x, ind_w;

  double w_sum = Rcpp::sum(w), tmp_wsum, tmp_xwsum, tmp_w;

  // beginning
  for (i = 0; i < size; i++) {
    tmp_xwsum = tmp_wsum = 0;

    // Fix this line
    for (ind_x = i + size, ind_w = w_size - 1; ind_x >= 0; ind_x--, ind_w--) { 
      tmp_w = w(ind_w);
      Rcpp::Rcout << "Loop at: " << ind_w << std::endl;
      tmp_wsum += tmp_w;
      tmp_xwsum += x(ind_x) * tmp_w;
    }

    res(i) = tmp_xwsum / tmp_wsum;
  }

  Rcpp::Rcout << "success" << std::endl;
  return res;
}

这就是所有人!