R：for循环中的foreach循环
R: for loop within a foreach loop
编辑：减小了数据集的大小
示例数据：
df <- data.frame(loc.id = rep(1:10, each = 80*36),
year = rep(rep(1980:2015, each = 80), times = 10),
day = rep(rep(1:80, times = 36),times = 10),
rain = runif(10*36*80, min = 0 , max = 5),
swc = runif(10*36*80,min = 0, max = 50),
SW.max = rep(runif(10, min = 100, max = 200), each = 80*36),
SW.ini = runif(10*36*80),
PETc = runif(10*36*80, min = 0 , max = 1.3),
SW = NA,
PAW = NA,
aetc = NA)
df 包含 1980-2015 年 10 个地点的每日数据（80 天）。
对于每个位置 X 年组合，我想做以下计算
list.result <- list() # create a list to store all results
ptm <- proc.time()
n <- 0
for(i in seq_along(unique(df$loc.id))){
location <- unique(df$loc.id)[i]
print(location)
for(j in seq_along(unique(df$year))){
yr <- unique(df$year)[j]
print(yr)
df_year <- df[df$loc.id == location & df$year == yr,] # subset data for location i and year y
# for the first row of data frame, i need to calculate some values
SW.ini <- df_year$SW.ini[1]
SW.max <- df_year$SW.max[1]
df_year$PAW[1] <- SW.ini + df_year$rain[1]
df_year$aetc[1] <- ifelse(df_year$PAW[1] >= df_year$swc[1],
df_year$PETc[1],(df_year$PAW[1]/df_year$swc[1])*df_year$PETc[1])
df_year$aetc[1] <- ifelse(df_year$aetc[1] > df_year$PAW[1], df_year$PAW[1], df_year$aetc[1])
df_year$SW[1] <- SW.ini + df_year$rain[1] - df_year$aetc[1]
df_year$SW[1] <- ifelse(df_year$SW[1] > SW.max, SW.max, ifelse(df_year$SW[1] < 0, 0,df_year$SW[1]))
# for row 2 till row n of df_year, I need to do this:
for (day in 2:nrow(df_year)){
df_year$PAW[day] <- df_year$SW[day - 1] + df_year$rain[day]
df_year$aetc[day] <- ifelse(df_year$PAW[day] >= df_year$swc[day], df_year$PETc[day], (df_year$PAW[day]/df_year$swc[day]) * df_year$PETc[day])
df_year$aetc[day] <- ifelse(df_year$aetc[day] > df_year$PAW[day], df_year$PAW[day],df_year$aetc[day])
df_year$SW[day] <- df_year$SW[day - 1] + df_year$rain[day] - df_year$aetc[day]
df_year$SW[day] <- ifelse(df_year$SW[day] > SW.max,SW.max, ifelse(df_year$SW[day] < 0, 0,df_year$SW[day]))
}
n <- n + 1
list.result[[n]] <- df_year
}}
proc.time() - ptm
user system elapsed
8.64 0.00 8.75
final.dat <- rbindlist(list.result)
这个循环是连续的，我认为它是 R 中 foreach 的一个很好的候选者。我还没有真正使用过
foreach 所以做了一些在线研究让我想到了这个：
library(doParallel)
cl <- makeCluster(4) # if I understood this correctly, it assings number of cores to be used
registerDoParallel(cl)
foreach(i = seq_along(unique(df$loc.id)) %dopar% {
list.result <- list()
for(j in seq_along(1980:2015)){
df_year <- df[df$loc.id == unique(df$loc.id)[i] & df$year == unique(df$year)[j],] # subset data for location i and year y
# for the first row of data frame, i need to calculate some values
SW.ini <- df_year$SW.ini[1]
SW.max <- df_year$SW.max[1]
df_year$PAW[1] <- SW.ini + df_year$rain[1]
df_year$aetc[1] <- ifelse(df_year$PAW[1] >= df_year$swc[1], df_year$PETc[1],(df_year$PAW[1]/df_year$swc[1])*df_year$PETc[1])
df_year$aetc[1] <- ifelse(df_year$aetc[1] > df_year$PAW[1], df_year$PAW[1], df_year$aetc[1])
df_year$SW[1] <- SW.ini + df_year$rain[1] - df_year$aetc[1]
df_year$SW[1] <- ifelse(df_year$SW[1] > SW.max, SW.max, ifelse(df_year$SW[1] < 0, 0,df_year$SW[1]))
# for row 2 till row n of df_year, I need to do this:
for (day in 2:nrow(df_year)){
df_year$PAW[day] <- df_year$SW[day - 1] + df_year$rain[day]
df_year$aetc[day] <- ifelse(df_year$PAW[day] >= df_year$swc[day], df_year$PETc[day], (df_year$PAW[day]/df_year$swc[day]) * df_year$PETc[day])
df_year$aetc[day] <- ifelse(df_year$aetc[day] > df_year$PAW[day], df_year$PAW[day],df_year$aetc[day])
df_year$SW[day] <- df_year$SW[day - 1] + df_year$rain[day] - df_year$aetc[day]
df_year$SW[day] <- ifelse(df_year$SW[day] > SW.max,SW.max, ifelse(df_year$SW[day] < 0, 0,df_year$SW[day]))
}
list.result[[j]] <- df_year
}
dat <- rbindlist(list.result)
fwrite(dat,paste0(i,"dat.csv"))
}
我的问题是：
1) 上面的数据是否适合foreach
2) foreach中有一个for循环。这有意义吗？
3) 如何让上面的foreach 运行 和return 的所有结果都变成
解决你的三个问题：
- 我不这么认为。 （计算效率更高的方法可以完全消除增加更多处理能力的需要。）
- 并行处理中的 for 循环本身没有什么坏处。 （事实上，每个块上需要完成的计算越多，并行处理的可能性就越大方法可以提高性能。）
- (以下方法不适用)
改用 Rcpp 和 data.table
使用 C++ 编译逻辑并使用 data.table 分组操作按组应用它可以从您的基线得到约 2,000 倍 speed-up，远远大于您希望通过并行化获得的结果。
在你原来的例子中，有 39,420,000 行，这在我的机器上执行 1.883 秒；在 28,800 行 的修订版上，执行时间为 0.004 秒
library(data.table)
library(Rcpp)
定义并编译一个 C++ 函数，CalcSW() 在 R 脚本中内联：
请注意：C/C++ 的计数从 0 开始，不像 R，它从 1- 开始- 这就是这里索引不同的原因
Rcpp::cppFunction('
List CalcSW(NumericVector SW_ini,
NumericVector SW_max,
NumericVector rain,
NumericVector swc,
NumericVector PETc) {
int n = SW_ini.length();
NumericVector SW(n);
NumericVector PAW(n);
NumericVector aetc(n);
double SW_ini_glob = SW_ini[0];
double SW_max_glob = SW_max[0];
SW[0] = SW_ini_glob;
PAW[0] = SW[0] + rain[0];
if (PAW[0] > swc[0]){
aetc[0] = PETc[0];
} else {
aetc[0] = PAW[0]/swc[0]*PETc[0];
}
if (aetc[0] > PAW[0]){
aetc[0] = PAW[0];
}
SW[0] = SW[0] + rain[0] - aetc[0];
if(SW[0] > SW_max_glob){
SW[0] = SW_max_glob;
}
if(SW[0] < 0){
SW[0] = 0;
}
for (int i = 1; i < n; i++) {
PAW[i] = SW[i-1] + rain[i];
if (PAW[i] > swc[i]){
aetc[i] = PETc[i];
} else {
aetc[i] = PAW[i]/swc[i]*PETc[i];
}
if (aetc[i] > PAW[i]){
aetc[i] = PAW[i];
}
SW[i] = SW[i-1] + rain[i] - aetc[i];
if(SW[i] > SW_max_glob){
SW[i] = SW_max_glob;
}
if(SW[i] < 0){
SW[i] = 0;
}
}
return Rcpp::List::create(Rcpp::Named("SW") = SW,
Rcpp::Named("PAW") = PAW,
Rcpp::Named("aetc") = aetc);
}')
创建data.table
df <- data.table(loc.id = rep(1:10, each = 80*36),
year = rep(rep(1980:2015, each = 80), times = 10),
day = rep(rep(1:80, times = 36),times = 10),
rain = runif(10*36*80, min = 0 , max = 5),
swc = runif(10*36*80,min = 0, max = 50),
SW_max = rep(runif(10, min = 100, max = 200), each = 80*36),
SW_ini = runif(10*36*80),
PETc = runif(10*36*80, min = 0 , max = 1.3),
SW = as.numeric(NA),
PAW = as.numeric(NA),
aetc = as.numeric(NA))
setkey(df, loc.id, year, day)
对loc.id和year的每个组合在df上执行函数CalcSW()，同时将返回值赋给三列：
system.time({
df[, c("SW","PAW","aetc") := CalcSW(SW_ini,
SW_max,
rain,
swc,
PETc), keyby = .(loc.id, year)]
})
...
user system elapsed
0.004 0.000 0.004
结果：
head(df)
...
loc.id year day rain swc SW_max SW_ini PETc SW PAW aetc
1: 1 1980 1 0.35813251 28.360715 177.3943 0.69116310 0.2870478 1.038675 1.049296 0.01062025
2: 1 1980 2 1.10331116 37.013022 177.3943 0.02742273 0.4412420 2.125335 1.396808 0.01665171
3: 1 1980 3 1.76680011 32.509970 177.3943 0.66273062 1.1071233 3.807561 2.483467 0.08457420
4: 1 1980 4 3.20966558 8.252797 177.3943 0.12220454 0.3496968 6.840713 4.165693 0.17651342
5: 1 1980 5 1.32498191 14.784203 177.3943 0.66381497 1.2168838 7.573160 7.198845 0.59253503
6: 1 1980 6 0.02547458 47.903637 177.3943 0.21871598 1.0864713 7.418750 7.931292 0.17988449
我不是 100% 肯定我完美地实现了你的逻辑，但逻辑应该非常简单，可以调整我可能遗漏的地方，我以与你布局的方式非常相似的方式实现它。
另一个注意事项：使用 auto-indenting 和代码突出显示 （无论您使用的是 RStudio 还是 Emacs） 编写 C++ 更容易如果您创建一个单独的文件，命名为 TestCode.cpp，格式如下。
然后，您可以使用 Rcpp::sourceCpp("TestCode.cpp") 在您的 R 脚本中编译您的函数，或者您可以将除前三行以外的所有内容作为字符串复制并粘贴到 [=34 的参数中=] 就像我上面做的那样。
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
List CalcSW(NumericVector SW_ini,
NumericVector SW_max,
NumericVector rain,
NumericVector swc,
NumericVector PETc) {
int n = SW_ini.length();
NumericVector SW(n);
NumericVector PAW(n);
NumericVector aetc(n);
double SW_ini_glob = SW_ini[0];
double SW_max_glob = SW_max[0];
SW[0] = SW_ini_glob;
PAW[0] = SW[0] + rain[0];
if (PAW[0] > swc[0]){
aetc[0] = PETc[0];
} else {
aetc[0] = PAW[0]/swc[0]*PETc[0];
}
if (aetc[0] > PAW[0]){
aetc[0] = PAW[0];
}
SW[0] = SW[0] + rain[0] - aetc[0];
if(SW[0] > SW_max_glob){
SW[0] = SW_max_glob;
}
if(SW[0] < 0){
SW[0] = 0;
}
for (int i = 1; i < n; i++) {
PAW[i] = SW[i-1] + rain[i];
if (PAW[i] > swc[i]){
aetc[i] = PETc[i];
} else {
aetc[i] = PAW[i]/swc[i]*PETc[i];
}
if (aetc[i] > PAW[i]){
aetc[i] = PAW[i];
}
SW[i] = SW[i-1] + rain[i] - aetc[i];
if(SW[i] > SW_max_glob){
SW[i] = SW_max_glob;
}
if(SW[i] < 0){
SW[i] = 0;
}
}
return Rcpp::List::create(Rcpp::Named("SW") = SW,
Rcpp::Named("PAW") = PAW,
Rcpp::Named("aetc") = aetc);
}
此代码替换内部循环
clamp <- function(x, low, high)
min(high, max(low, x))
fill1 <- function(df) {
rain <- df$rain
swc <- df$swc
PETc <- df$PETc
SW0 <- df$SW.ini[1]
SW.max <- df$SW.max[1]
SW <- PAW <- aetc <- numeric(nrow(df))
for (day in seq_along(rain)) {
PAW[day] <- SW0 + rain[day]
if (PAW[day] >= swc[day]) {
aetc0 <- PETc[day]
} else {
aetc0 <- (PAW[day] / swc[day]) * PETc[day]
}
aetc[day] <- min(PAW[day], aetc0)
SW0 <- SW[day] <- clamp(PAW[day] - aetc[day], 0, SW.max)
}
list(SW = SW, PAW = PAW, aetc = aetc)
}
并且比原始问题中的实现快约 60 倍。请注意，这是 C++ 中采用的方法，即分配和更新新向量，而不是 data.frame 的现有部分；这是性能差异的很大一部分，并且可以在没有 Rcpp 的情况下获得好处。
这是对 location.year x 天矩阵
进行迭代的泛化（非常简单的测试！）
pclamp <- function(x, low, high)
pmin(high, pmax(low, x))
fill2 <- function(rain, swc, PETc, SW0, SW.max) {
SW <- PAW <- aetc <- matrix(0, nrow = nrow(rain), ncol = ncol(rain))
for (day in seq_len(ncol(rain))) {
PAW[, day] <- SW0 + rain[, day]
aetc0 <- PETc[, day]
idx <- PAW[, day] < swc[, day]
aetc0[idx] <- (PAW[idx, day] / swc[idx, day]) * PETc[idx, day]
aetc[, day] <- pmin(PAW[, day], aetc0)
SW0 <- SW[, day] <- pclamp(PAW[, day] - aetc[, day], 0, SW.max)
}
list(SW = SW, PAW = PAW, aetc = aetc)
}
使用原始输入，假设输入按年份、位置和日期排序
days <- 80
rain <- matrix(df$rain, ncol=days, byrow=TRUE)
swc <- matrix(df$swc, ncol=days, byrow=TRUE)
PETc <- matrix(df$PETc, ncol=days, byrow=TRUE)
SW.ini <- df$SW.ini[df$day == 1]
SW.max <- df$SW.max[df$day == 1]
result <- fill2(rain, swc, PETc, SW.ini, SW.max)
在 per-location.date 的基础上，对于问题中的数据子集，它比 fill1() 快大约 15 倍。对示例数据的操作大约需要 10 毫秒，对完整数据大约需要 10 秒——比 Matt 的 C++ 解决方案慢 5 倍，但仍然比原始的和采用基本 R 技术的改进许多不同领域的代码有很大的改进。
编辑：减小了数据集的大小
示例数据：
df <- data.frame(loc.id = rep(1:10, each = 80*36),
year = rep(rep(1980:2015, each = 80), times = 10),
day = rep(rep(1:80, times = 36),times = 10),
rain = runif(10*36*80, min = 0 , max = 5),
swc = runif(10*36*80,min = 0, max = 50),
SW.max = rep(runif(10, min = 100, max = 200), each = 80*36),
SW.ini = runif(10*36*80),
PETc = runif(10*36*80, min = 0 , max = 1.3),
SW = NA,
PAW = NA,
aetc = NA)
df 包含 1980-2015 年 10 个地点的每日数据（80 天）。
对于每个位置 X 年组合，我想做以下计算
list.result <- list() # create a list to store all results
ptm <- proc.time()
n <- 0
for(i in seq_along(unique(df$loc.id))){
location <- unique(df$loc.id)[i]
print(location)
for(j in seq_along(unique(df$year))){
yr <- unique(df$year)[j]
print(yr)
df_year <- df[df$loc.id == location & df$year == yr,] # subset data for location i and year y
# for the first row of data frame, i need to calculate some values
SW.ini <- df_year$SW.ini[1]
SW.max <- df_year$SW.max[1]
df_year$PAW[1] <- SW.ini + df_year$rain[1]
df_year$aetc[1] <- ifelse(df_year$PAW[1] >= df_year$swc[1],
df_year$PETc[1],(df_year$PAW[1]/df_year$swc[1])*df_year$PETc[1])
df_year$aetc[1] <- ifelse(df_year$aetc[1] > df_year$PAW[1], df_year$PAW[1], df_year$aetc[1])
df_year$SW[1] <- SW.ini + df_year$rain[1] - df_year$aetc[1]
df_year$SW[1] <- ifelse(df_year$SW[1] > SW.max, SW.max, ifelse(df_year$SW[1] < 0, 0,df_year$SW[1]))
# for row 2 till row n of df_year, I need to do this:
for (day in 2:nrow(df_year)){
df_year$PAW[day] <- df_year$SW[day - 1] + df_year$rain[day]
df_year$aetc[day] <- ifelse(df_year$PAW[day] >= df_year$swc[day], df_year$PETc[day], (df_year$PAW[day]/df_year$swc[day]) * df_year$PETc[day])
df_year$aetc[day] <- ifelse(df_year$aetc[day] > df_year$PAW[day], df_year$PAW[day],df_year$aetc[day])
df_year$SW[day] <- df_year$SW[day - 1] + df_year$rain[day] - df_year$aetc[day]
df_year$SW[day] <- ifelse(df_year$SW[day] > SW.max,SW.max, ifelse(df_year$SW[day] < 0, 0,df_year$SW[day]))
}
n <- n + 1
list.result[[n]] <- df_year
}}
proc.time() - ptm
user system elapsed
8.64 0.00 8.75
final.dat <- rbindlist(list.result)
这个循环是连续的，我认为它是 R 中 foreach 的一个很好的候选者。我还没有真正使用过 foreach 所以做了一些在线研究让我想到了这个：
library(doParallel)
cl <- makeCluster(4) # if I understood this correctly, it assings number of cores to be used
registerDoParallel(cl)
foreach(i = seq_along(unique(df$loc.id)) %dopar% {
list.result <- list()
for(j in seq_along(1980:2015)){
df_year <- df[df$loc.id == unique(df$loc.id)[i] & df$year == unique(df$year)[j],] # subset data for location i and year y
# for the first row of data frame, i need to calculate some values
SW.ini <- df_year$SW.ini[1]
SW.max <- df_year$SW.max[1]
df_year$PAW[1] <- SW.ini + df_year$rain[1]
df_year$aetc[1] <- ifelse(df_year$PAW[1] >= df_year$swc[1], df_year$PETc[1],(df_year$PAW[1]/df_year$swc[1])*df_year$PETc[1])
df_year$aetc[1] <- ifelse(df_year$aetc[1] > df_year$PAW[1], df_year$PAW[1], df_year$aetc[1])
df_year$SW[1] <- SW.ini + df_year$rain[1] - df_year$aetc[1]
df_year$SW[1] <- ifelse(df_year$SW[1] > SW.max, SW.max, ifelse(df_year$SW[1] < 0, 0,df_year$SW[1]))
# for row 2 till row n of df_year, I need to do this:
for (day in 2:nrow(df_year)){
df_year$PAW[day] <- df_year$SW[day - 1] + df_year$rain[day]
df_year$aetc[day] <- ifelse(df_year$PAW[day] >= df_year$swc[day], df_year$PETc[day], (df_year$PAW[day]/df_year$swc[day]) * df_year$PETc[day])
df_year$aetc[day] <- ifelse(df_year$aetc[day] > df_year$PAW[day], df_year$PAW[day],df_year$aetc[day])
df_year$SW[day] <- df_year$SW[day - 1] + df_year$rain[day] - df_year$aetc[day]
df_year$SW[day] <- ifelse(df_year$SW[day] > SW.max,SW.max, ifelse(df_year$SW[day] < 0, 0,df_year$SW[day]))
}
list.result[[j]] <- df_year
}
dat <- rbindlist(list.result)
fwrite(dat,paste0(i,"dat.csv"))
}
我的问题是：
1) 上面的数据是否适合foreach
2) foreach中有一个for循环。这有意义吗？
3) 如何让上面的foreach 运行 和return 的所有结果都变成
解决你的三个问题：
- 我不这么认为。 （计算效率更高的方法可以完全消除增加更多处理能力的需要。）
- 并行处理中的 for 循环本身没有什么坏处。 （事实上，每个块上需要完成的计算越多，并行处理的可能性就越大方法可以提高性能。）
- (以下方法不适用)
改用 Rcpp 和 data.table
使用 C++ 编译逻辑并使用 data.table 分组操作按组应用它可以从您的基线得到约 2,000 倍 speed-up，远远大于您希望通过并行化获得的结果。
在你原来的例子中，有 39,420,000 行，这在我的机器上执行 1.883 秒；在 28,800 行 的修订版上，执行时间为 0.004 秒
library(data.table)
library(Rcpp)
定义并编译一个 C++ 函数，CalcSW() 在 R 脚本中内联：
请注意：C/C++ 的计数从 0 开始，不像 R，它从 1- 开始- 这就是这里索引不同的原因
Rcpp::cppFunction('
List CalcSW(NumericVector SW_ini,
NumericVector SW_max,
NumericVector rain,
NumericVector swc,
NumericVector PETc) {
int n = SW_ini.length();
NumericVector SW(n);
NumericVector PAW(n);
NumericVector aetc(n);
double SW_ini_glob = SW_ini[0];
double SW_max_glob = SW_max[0];
SW[0] = SW_ini_glob;
PAW[0] = SW[0] + rain[0];
if (PAW[0] > swc[0]){
aetc[0] = PETc[0];
} else {
aetc[0] = PAW[0]/swc[0]*PETc[0];
}
if (aetc[0] > PAW[0]){
aetc[0] = PAW[0];
}
SW[0] = SW[0] + rain[0] - aetc[0];
if(SW[0] > SW_max_glob){
SW[0] = SW_max_glob;
}
if(SW[0] < 0){
SW[0] = 0;
}
for (int i = 1; i < n; i++) {
PAW[i] = SW[i-1] + rain[i];
if (PAW[i] > swc[i]){
aetc[i] = PETc[i];
} else {
aetc[i] = PAW[i]/swc[i]*PETc[i];
}
if (aetc[i] > PAW[i]){
aetc[i] = PAW[i];
}
SW[i] = SW[i-1] + rain[i] - aetc[i];
if(SW[i] > SW_max_glob){
SW[i] = SW_max_glob;
}
if(SW[i] < 0){
SW[i] = 0;
}
}
return Rcpp::List::create(Rcpp::Named("SW") = SW,
Rcpp::Named("PAW") = PAW,
Rcpp::Named("aetc") = aetc);
}')
创建data.table
df <- data.table(loc.id = rep(1:10, each = 80*36),
year = rep(rep(1980:2015, each = 80), times = 10),
day = rep(rep(1:80, times = 36),times = 10),
rain = runif(10*36*80, min = 0 , max = 5),
swc = runif(10*36*80,min = 0, max = 50),
SW_max = rep(runif(10, min = 100, max = 200), each = 80*36),
SW_ini = runif(10*36*80),
PETc = runif(10*36*80, min = 0 , max = 1.3),
SW = as.numeric(NA),
PAW = as.numeric(NA),
aetc = as.numeric(NA))
setkey(df, loc.id, year, day)
对loc.id和year的每个组合在df上执行函数CalcSW()，同时将返回值赋给三列：
system.time({
df[, c("SW","PAW","aetc") := CalcSW(SW_ini,
SW_max,
rain,
swc,
PETc), keyby = .(loc.id, year)]
})
...
user system elapsed
0.004 0.000 0.004
结果：
head(df)
...
loc.id year day rain swc SW_max SW_ini PETc SW PAW aetc
1: 1 1980 1 0.35813251 28.360715 177.3943 0.69116310 0.2870478 1.038675 1.049296 0.01062025
2: 1 1980 2 1.10331116 37.013022 177.3943 0.02742273 0.4412420 2.125335 1.396808 0.01665171
3: 1 1980 3 1.76680011 32.509970 177.3943 0.66273062 1.1071233 3.807561 2.483467 0.08457420
4: 1 1980 4 3.20966558 8.252797 177.3943 0.12220454 0.3496968 6.840713 4.165693 0.17651342
5: 1 1980 5 1.32498191 14.784203 177.3943 0.66381497 1.2168838 7.573160 7.198845 0.59253503
6: 1 1980 6 0.02547458 47.903637 177.3943 0.21871598 1.0864713 7.418750 7.931292 0.17988449
我不是 100% 肯定我完美地实现了你的逻辑，但逻辑应该非常简单，可以调整我可能遗漏的地方，我以与你布局的方式非常相似的方式实现它。
另一个注意事项：使用 auto-indenting 和代码突出显示 （无论您使用的是 RStudio 还是 Emacs） 编写 C++ 更容易如果您创建一个单独的文件，命名为 TestCode.cpp，格式如下。
然后，您可以使用 Rcpp::sourceCpp("TestCode.cpp") 在您的 R 脚本中编译您的函数，或者您可以将除前三行以外的所有内容作为字符串复制并粘贴到 [=34 的参数中=] 就像我上面做的那样。
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
List CalcSW(NumericVector SW_ini,
NumericVector SW_max,
NumericVector rain,
NumericVector swc,
NumericVector PETc) {
int n = SW_ini.length();
NumericVector SW(n);
NumericVector PAW(n);
NumericVector aetc(n);
double SW_ini_glob = SW_ini[0];
double SW_max_glob = SW_max[0];
SW[0] = SW_ini_glob;
PAW[0] = SW[0] + rain[0];
if (PAW[0] > swc[0]){
aetc[0] = PETc[0];
} else {
aetc[0] = PAW[0]/swc[0]*PETc[0];
}
if (aetc[0] > PAW[0]){
aetc[0] = PAW[0];
}
SW[0] = SW[0] + rain[0] - aetc[0];
if(SW[0] > SW_max_glob){
SW[0] = SW_max_glob;
}
if(SW[0] < 0){
SW[0] = 0;
}
for (int i = 1; i < n; i++) {
PAW[i] = SW[i-1] + rain[i];
if (PAW[i] > swc[i]){
aetc[i] = PETc[i];
} else {
aetc[i] = PAW[i]/swc[i]*PETc[i];
}
if (aetc[i] > PAW[i]){
aetc[i] = PAW[i];
}
SW[i] = SW[i-1] + rain[i] - aetc[i];
if(SW[i] > SW_max_glob){
SW[i] = SW_max_glob;
}
if(SW[i] < 0){
SW[i] = 0;
}
}
return Rcpp::List::create(Rcpp::Named("SW") = SW,
Rcpp::Named("PAW") = PAW,
Rcpp::Named("aetc") = aetc);
}
此代码替换内部循环
clamp <- function(x, low, high)
min(high, max(low, x))
fill1 <- function(df) {
rain <- df$rain
swc <- df$swc
PETc <- df$PETc
SW0 <- df$SW.ini[1]
SW.max <- df$SW.max[1]
SW <- PAW <- aetc <- numeric(nrow(df))
for (day in seq_along(rain)) {
PAW[day] <- SW0 + rain[day]
if (PAW[day] >= swc[day]) {
aetc0 <- PETc[day]
} else {
aetc0 <- (PAW[day] / swc[day]) * PETc[day]
}
aetc[day] <- min(PAW[day], aetc0)
SW0 <- SW[day] <- clamp(PAW[day] - aetc[day], 0, SW.max)
}
list(SW = SW, PAW = PAW, aetc = aetc)
}
并且比原始问题中的实现快约 60 倍。请注意，这是 C++ 中采用的方法，即分配和更新新向量，而不是 data.frame 的现有部分；这是性能差异的很大一部分，并且可以在没有 Rcpp 的情况下获得好处。
这是对 location.year x 天矩阵
进行迭代的泛化（非常简单的测试！）pclamp <- function(x, low, high)
pmin(high, pmax(low, x))
fill2 <- function(rain, swc, PETc, SW0, SW.max) {
SW <- PAW <- aetc <- matrix(0, nrow = nrow(rain), ncol = ncol(rain))
for (day in seq_len(ncol(rain))) {
PAW[, day] <- SW0 + rain[, day]
aetc0 <- PETc[, day]
idx <- PAW[, day] < swc[, day]
aetc0[idx] <- (PAW[idx, day] / swc[idx, day]) * PETc[idx, day]
aetc[, day] <- pmin(PAW[, day], aetc0)
SW0 <- SW[, day] <- pclamp(PAW[, day] - aetc[, day], 0, SW.max)
}
list(SW = SW, PAW = PAW, aetc = aetc)
}
使用原始输入，假设输入按年份、位置和日期排序
days <- 80
rain <- matrix(df$rain, ncol=days, byrow=TRUE)
swc <- matrix(df$swc, ncol=days, byrow=TRUE)
PETc <- matrix(df$PETc, ncol=days, byrow=TRUE)
SW.ini <- df$SW.ini[df$day == 1]
SW.max <- df$SW.max[df$day == 1]
result <- fill2(rain, swc, PETc, SW.ini, SW.max)
在 per-location.date 的基础上，对于问题中的数据子集，它比 fill1() 快大约 15 倍。对示例数据的操作大约需要 10 毫秒，对完整数据大约需要 10 秒——比 Matt 的 C++ 解决方案慢 5 倍，但仍然比原始的和采用基本 R 技术的改进许多不同领域的代码有很大的改进。