请求将R中行依赖式for循环代码转换为Rcpp实现

完整Rcpp代码 #

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix process_matrix(NumericMatrix temp, double const1) {
  int n = temp.nrow();
  
  for (int p = 1; p <= n; ++p) {
    if (p == 1) {
      // 对应R代码中p==1的计算逻辑
      temp(p, 4) = std::max(std::min(temp(p, 2), temp(p, 1)), 0.0);
      temp(p, 5) = std::max(temp(p, 3) + (0.0 - const1), 0.0);
      temp(p, 6) = temp(p, 1) - temp(p, 4) - temp(p, 5);
    } else {
      // 对应R代码中p>1的计算逻辑
      temp(p, 4) = std::max(std::min(temp(p, 2), temp(p, 1) + temp(p-1, 6)), 0.0);
      temp(p, 5) = std::max(temp(p, 3) + (temp(p-1, 6) - const1), 0.0);
      temp(p, 6) = temp(p-1, 6) + temp(p, 1) - temp(p, 4) - temp(p, 5);
    }
  }
  
  return temp;
}

关键说明与使用方法 #

• 索引对应：Rcpp的NumericMatrix支持1-based索引，和你的R代码完全一致，直接对应temp[p,4]（R）与temp(p,4)（Rcpp）的写法，无需转换行列号。
• 逻辑复刻：循环分支、计算逻辑完全照搬原R代码，你可以逐行对比两者的对应关系，快速理解Rcpp的写法逻辑。
• R中调用：将上述代码保存为.cpp文件后，用sourceCpp("你的文件名.cpp")加载，然后直接调用process_matrix(temp, const1)即可得到与原循环相同的结果。
• 性能优势：这类依赖前序结果的循环无法并行，但Rcpp避免了R解释型循环的开销，对于数千行的矩阵，速度会比原R循环提升数倍。

结果验证代码（R端） #

# 原输入数据
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

# 原R循环计算结果
temp_r <- temp
for (p in 1:nrow(temp_r)) {
  if (p==1) {
    temp_r[p, 4] <- max(min(temp_r[p, 2], temp_r[p, 1]), 0)
    temp_r[p, 5] <- max(temp_r[p, 3] + (0 - const1), 0)
    temp_r[p, 6] <- temp_r[p, 1] - temp_r[p, 4] - temp_r[p, 5]
  }
  if (p>1) {
    temp_r[p, 4] <- max(min(temp_r[p, 2], temp_r[p, 1] + temp_r[p-1, 6]), 0)
    temp_r[p, 5] <- max(temp_r[p, 3] + (temp_r[p-1, 6] - const1), 0)
    temp_r[p, 6] <- temp_r[p-1, 6] + temp_r[p, 1] - temp_r[p, 4] - temp_r[p, 5]
  }
}

# Rcpp计算结果
temp_rcpp <- process_matrix(temp, const1)

# 验证结果一致
all.equal(temp_r, temp_rcpp)

内容的提问来源于stack exchange，提问作者Greg