将 R 中的 for 循环转换为 Rcpp

Question

我一直在尝试使用更高效的数据结构和并行处理以及其他一些东西。我在将脚本从大约 60 秒内运行到大约 9 秒内运行方面取得了很好的进展。

不过，我一生都无法理解的一件事是在 Rcpp 中编写一个循环。具体来说，一个循环根据前一行结果逐行计算并随时更新数据。

想知道是否有人可以将我的代码转换为 Rcpp，这样我就可以通过一个我非常熟悉的示例进行反向工程并弄清楚它是如何完成的。

这是一个循环，计算每行 3 个变量的结果。第 1 行必须单独计算，然后第 2 行以后根据当前行和先前行的值进行计算。

这个示例代码只有 6 行长，但我的原始代码有数千行：

temp <- matrix(c(0, 0, 0, 2.211, 2.345, 0, 0.8978, 1.0452, 1.1524, 0.4154, 
                 0.7102, 0.8576, 0, 0, 0, 1.7956, 1.6348, 0, 
                 rep(NA, 18)), ncol=6, nrow=6)
const1 <- 0.938

for (p in 1:nrow(temp)) {
  if (p==1) {
    temp[p, 4] <- max(min(temp[p, 2], 
                          temp[p, 1]),
                      0)
    temp[p, 5] <- max(temp[p, 3] + (0 - const1), 
                      0)
    temp[p, 6] <- temp[p, 1] - temp[p, 4] - temp[p, 5]
  }
  if (p>1) {
    temp[p, 4] <- max(min(temp[p, 2], 
                          temp[p, 1] + temp[p-1, 6]),
                      0)
    temp[p, 5] <- max(temp[p, 3] + (temp[p-1, 6] - const1),
                      0)
    temp[p, 6] <- temp[p-1, 6] + temp[p, 1] - temp[p, 4] - temp[p, 5]
  }
}

在此先感谢，希望这需要具有 Rcpp 技能的人只需一两分钟！

编辑：感谢您的帮助。只是想知道如果 x 是 6 个向量的列表，而不是 6 列的矩阵，如何布置它......我在想这样的事情，但不确定如何让它工作：

List getResult(  ???  x, double const1) {
  for (int p=1; p<x.nrow(); p++){
    x[3](p) = std::max(std::min(x[p](p), x[0](p) + x[5](p - 1)), 0.0);
    x[4](p) = std::max(x[2](p) + (a[5](p - 1) - const1), 0.0);
    x[5](p) = x[5](p - 1) + x[0](p) - x[3](p) - x[4](p);
  }
  return x
}

Answer 1

这是一个示例 Rcpp 等效代码：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix getResult(NumericMatrix x, double const1){
  for (int p = 0; p < x.nrow(); p++){
    if (p == 0){
      x(p, 3) = std::max(std::min(x(p, 1), x(p, 0)), 0.0);
      x(p, 4) = std::max(x(p, 2) + (0.0 - const1), 0.0);
      x(p, 5) = x(p, 0) - x(p, 3) - x(p, 4);
    }
    if (p > 0){
      x(p, 3) = std::max(std::min(x(p, 1), x(p, 0) + x(p - 1, 5)), 0.0);
      x(p, 4) = std::max(x(p, 2) + (x(p - 1, 5) - const1), 0.0);
      x(p, 5) = x(p - 1, 5) + x(p, 0) - x(p, 3) - x(p, 4);
    }
  }
  return x;
}

几点注意事项：

• 将其保存在文件中并在会话中执行Rcpp::sourceCpp("myCode.cpp") 以编译它并使其在会话中可用。
• 这里我们使用NumericMatrix来表示矩阵。
• 您会看到我们分别调用了std::max 和std::min。这些函数需要两种常见的数据类型，即如果我们使用max(x, y)，则x 和y 必须属于同一类型。数字矩阵条目是double（我相信），所以你需要提供一个double；因此，从 0（C++ 中的 int）更改为 0.0（双精度）
• 在 C++ 中，索引从 0 而不是 1 开始。因此，您将像 temp[1, 4] 这样的 R 代码转换为 temp(0, 3)
• 查看http://adv-r.had.co.nz/Rcpp.html 了解更多信息以支持您的开发