R：foreach循环中的for循环
Posted
技术标签:
【中文标题】R：foreach循环中的for循环【英文标题】：R: for loop within a foreach loop 【发布时间】：2018-08-14 15:47:19 【问题描述】：编辑：减小数据集的大小
样本数据：
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 运行并返回所有结果
【问题讨论】：
我会为 1 个位置编写一个函数，然后使用lapply 或 purrr::map 循环遍历所有 3000 个位置。那将摆脱 1 个循环
对于第二个循环，看起来您可以使用Reduce。有关示例，请参见以下链接：***.com/questions/40412516/… | ***.com/questions/34624110/…
我们了解到您的数据集很大，上面的代码很慢。您能否将上述样本的大小从 3900 万行减少到大约 100 行。这将允许其他人运行您的代码并提供经过测试的改进建议。
好的。我能做到。给我1分钟
如果数据是 (year.location) x 天矩阵，那么按天迭代可以跨 year.location 向量化，从而导致 3000 x 15 的加速。
【参考方案1】：
解决您的三个问题：
-
我不这么认为。 （计算效率更高的方法可以完全消除增加更多处理能力的需要。）
并行处理中的 for 循环本身并没有什么坏处。 （事实上，需要对每个块进行的计算越多，并行方法就越有可能提高性能。）时间>
（如果您使用以下方法，则不适用）
改用Rcpp 和data.table
使用 C++ 编译逻辑并使用 data.table 分组操作按组应用它可以比您的基线提高约 2,000 倍的速度，远远超过您希望通过并行化获得的速度。
在您的原始示例中，它有 39,420,000 行，它在我的机器上在 1.883 秒内执行；而在具有 28,800 行的修订版上，这将在 0.004 秒
内执行library(data.table)
library(Rcpp)
在R脚本中定义并编译一个C++函数，CalcSW()内联：
注意：C/C++ 的计数从 0 开始，与 R 不同，R 从 1 开始——这就是索引不同的原因 p>
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);
')
创建数据表
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% 肯定我完美地实现了您的逻辑，但是该逻辑应该很容易调整我可能遗漏的地方，我以与您的布局非常相似的方式实现它。
另一个注意事项：使用自动缩进和代码突出显示编写C++ 会更容易（无论您使用的是 RStudio 还是 Emacs）如果您创建一个单独的文件，命名为 something像TestCode.cppformatted 如下所示。
然后，您可以使用Rcpp::sourceCpp("TestCode.cpp") 在您的 R 脚本中编译您的函数，或者您可以将除前三行以外的所有内容作为字符串复制并粘贴到 Rcpp::cppFunction() 的参数中，就像我在上面所做的那样.
#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);
【讨论】：
谢谢马特。这是一个非常详细的答案。由于我不熟悉 Rcpp，因此我将不得不通过它，并且一旦我设法理解它就会接受（赞成）您的答案。请多多包涵。再次感谢您的宝贵时间。 一点也不着急！我自己不经常使用Rcpp，所以这是复习一些基础知识的好机会。我在这里使用它的唯一原因是因为这个问题有一个元素 （依赖于前一行的计算） 使得 for 循环不可避免——这些是编译后的 c++ 可以真正发挥作用的情况。我写的 99.9% 的代码都是简单的 R + data.table，因为它通常足够快，但是 @f-privé 在 this question 上的回答启发我考虑使用它来解决这类问题。
这被证明是我学到的最有用的东西。所以非常感谢。小点：1）PAW[i] = SW[i-1] + rain[0] 应该是PAW[i] = SW[i-1] + rain[i]，如果我正确理解了这段代码。 2）n = SW_ini.length() 是做什么的？
乐于助人！ 1）我的错字，编辑以反映您的评论。 2) 这是 R 的 length() 函数的 C++ 等价物。它将n定义为一个整数，表示输入向量SW_ini的长度
C 和 C++ 是静态类型的编译语言，而 R 是动态类型的解释语言。预先将逻辑编译为机器代码指令确实可以使其本质上更快地用于像这样需要执行数百万次的简单 for 循环。话虽如此，我建议阅读整页 csgillespie.github.io/efficientR/performance.html（特别关注 profvis），您可以在 R 中做许多其他事情（即使用 data.table 而不是基本 R 数据帧） 以获得数量级的加速。【参考方案2】：
这段代码替换了内部循环
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 day 矩阵进行迭代的概括（非常简单的测试！）
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 技术，这将改进许多不同领域的代码。
【讨论】：
谢谢马丁。让我测试这些解决方案，然后再回复您。问候以上是关于R：foreach循环中的for循环的主要内容，如果未能解决你的问题，请参考以下文章