R：设置初始条件的for循环的dplyr解决方案

Question

我有一个一年有 40 天的数据和一些数据

set.seed(123)
df <- data.frame(day = 1:40,rain = runif(40,min = 0, max = 3), petc = runif(40, min = 0.3, max = 8),swc = runif(40, min = 27.01, max = 117.43))

我想为每一天计算另一个名为 aetc 的变量，计算如下：

SW.ini <- 2 # setting some initial values 
SW.max <- 5
SW.min <- 0

第 1 天，

1) 确定一个名为PAW(day1) = SW.ini + rain(day1)

的变量

2) 如果PAW(day1) >= SWC(day1), aetc(day1) = petc(day1);

If `PAW(day1) < SWC(day1), aetc(day1) = PAW(day1)/SWC(day1) * petc(day1)`

3) 检查是否 aetc(day1) > PAW(day1). If yes, aetc(day1) = paw(day1)

4) 更新SW(day1) = SW.ini + rain(day1) - aetc(day1)

5) 如果 SW(day1) > SW.max, SW(day1) = SW.max. Similarly ifSW(day1) < SW.min, SW(day1) = SW.min`

第二天重复

1) 确定PAW(day2) = SW(day1) + rain(day2)
2）如果PAW(day2) >= SWC(day2), aetc(day2) = petc(day2)；如果PAW(day2) < SWC(day2), aetc(day2) = PAW(day2)/SWC(day2) * petc(day2)

3) 检查是否 aetc(day2) > PAW(day2)。如果是，aetc(day2) = paw(day2)

4) 更新SW(day2) = SW(day1) + rain(day2) - aetc(day2)

5) 如果 SW(day2) > SW.max, SW(day2) = SW.max. Similarly ifSW(day2) < SW.min, 西南（第 2 天）= SW.min`

这是我的优雅 for 循环：

      df$PAW <- NA
      df$aetc <- NA
      df$SW <- NA

      df$PAW[1] <- SW.ini + df$rain[1]

      df$aetc[1] <- ifelse(df$PAW[1] >= df$swc[1], df$petc[1],(df$PAW[1]/df$swc[1])*df$petc[1])
      df$aetc[1] <- ifelse(df$aetc[1] > df$PAW[1], df$PAW[1], df$aetc[1])
      df$SW[1] <- SW.ini + df$rain[1] -  df$aetc[1]
      df$SW[1] <- ifelse(df$SW[1] > SW.max, SW.max, ifelse(df$SW[1] < 0, 0,df$SW[1]))

      for (day in 2:nrow(df)){

        df$PAW[day] <- df$SW[day - 1] + df$rain[day]
        df$aetc[day] <- ifelse(df$PAW[day] >= df$swc[day], df$petc[day], (df$PAW[day]/df$swc[day]) * df$petc[day])
        df$aetc[day] <- ifelse(df$aetc[day] > df$PAW[day], df$PAW[day],df$aetc[day])
        df$SW[day] <- df$SW[day - 1] + df$rain[day] -  df$aetc[day]
        df$SW[day] <- ifelse(df$SW[day] > SW.max,SW.max, ifelse(df$SW[day] < 0, 0,df$SW[day]))
      }

我的问题是这只是一年的数据，我想要运行多年。

      set.seed(123)
      df <- data.frame(year = 1980:2015, day = rep(1:40, each = 36),rain = 
      runif(40*36,min = 0, max = 3), petc = runif(40*36, min = 0.3, max = 8),swc = runif(40*36, min = 27.01, max = 117.43))

所以我想做类似

的事情

                df %>% group_by(year) # and then run the above function for each year.

是否有 dplyr 或任何其他解决方案？

谢谢

Answer 1

您可以将代码包装在另一个 for 循环中，并将每年的 df 保存在列表中：

library(tidyverse)
lst <- vector("list", length(unique(df$year)))
for (i in seq_along(unique(df$year))) {
    df_year <- df %>% filter(year == unique(df$year)[[i]])

    # rest of code with df_year replacing df

    lst[[i]] <- df_year
}
final_df <- bind_rows(lst)

Answer 2

Note: I originally posted this answer on your follow up question, , but after seeing this one, it seems this answer is far more relevant here. (I don't address anything related to parallelizing in my answer, which was the topic of your follow up).

使用 `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[0];

    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[0];

    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);
}

Answer 3

Matt 的 data.table 插图很好地说明了 data.table 有多快，因为它在没有副本和移动数据的情况下就地进行计算。

但是，要回答有关使用管道的问题的关键，您可以使用 group_by 和 do 来完成您想要的（尽管比 data.table 慢得多）

下面我设置了与 Matt 相同的虚拟数据。然后我使用你的函数（但大小写固定在 PETc 上）。它并不快，但很容易跟上。

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) 
                 )

my_fun <- function(df){
  SW.ini <- 2 # setting some initial values 
  SW.max <- 5
  SW.min <- 0

  df$PAW <- NA
  df$aetc <- NA
  df$SW <- NA

  df$PAW[1] <- SW.ini + df$rain[1]

  df$aetc[1] <- ifelse(df$PAW[1] >= df$swc[1], df$PETc[1],(df$PAW[1]/df$swc[1])*df$PETc[1])
  df$aetc[1] <- ifelse(df$aetc[1] > df$PAW[1], df$PAW[1], df$aetc[1])
  df$SW[1] <- SW.ini + df$rain[1] -  df$aetc[1]
  df$SW[1] <- ifelse(df$SW[1] > SW.max, SW.max, ifelse(df$SW[1] < 0, 0,df$SW[1]))

  for (day in 2:nrow(df)){

    df$PAW[day] <- df$SW[day - 1] + df$rain[day]
    df$aetc[day] <- ifelse(df$PAW[day] >= df$swc[day], df$PETc[day], (df$PAW[day]/df$swc[day]) * df$PETc[day])
    df$aetc[day] <- ifelse(df$aetc[day] > df$PAW[day], df$PAW[day],df$aetc[day])
    df$SW[day] <- df$SW[day - 1] + df$rain[day] -  df$aetc[day]
    df$SW[day] <- ifelse(df$SW[day] > SW.max,SW.max, ifelse(df$SW[day] < 0, 0,df$SW[day]))
  }
  return(df)
}


library(tictoc)
library(tidyverse)


tic()
df  %>% 
  group_by(year) %>%
  do(my_fun(.)) -> 
  out
toc()
#> 5.075 sec elapsed

R：设置初始条件的for循环的dplyr解决方案

R: dplyr solution for for-loop with initial conditions set

r

rcpp

dplyr

data.table

purrr

使用 `Rcpp` 和 `data.table`

R：设置初始条件的for循环的dplyr解决方案

R: dplyr solution for for-loop with initial conditions set

r

rcpp

dplyr

data.table

purrr

使用 Rcpp 和 data.table

使用 `Rcpp` 和 `data.table`