文档建议反馈控制台
首页
学习
活动
专区
圈层
工具
MCP广场
文章/答案/技术大牛
发布
首页
学习
活动
专区
圈层
工具
MCP广场
返回腾讯云官网
社区首页 >问答首页 >具有初始条件的循环的dplyr解

问具有初始条件的循环的dplyr解
EN

Stack Overflow用户
提问于 2018-03-05 19:51:17
回答 3查看 301关注 0票数 2

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

代码语言：javascript
运行
复制
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的变量，计算如下：

代码语言：javascript
运行
复制
SW.ini <- 2 # setting some initial values 
SW.max <- 5
SW.min <- 0

第一天，

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

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

代码语言：javascript
运行
复制
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`

重复第2天

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，SW(day2) = SW.min`

下面是我优雅的for循环来完成这个任务：

代码语言：javascript
运行
复制
      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]))
      }

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

代码语言：javascript
运行
复制
      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))

所以我想做一件事

代码语言：javascript
运行
复制
                df %>% group_by(year) # and then run the above function for each year.

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

谢谢

r
dplyr
data.table
rcpp
purrr
EN

回答 3

Stack Overflow用户

回答已采纳

发布于 2018-03-06 22:51:47

注意:我最初在您的后续问题R: for loop within a foreach loop上发布了这个答案，但是在看到这个问题之后，这个答案在这里似乎更有意义。(在我的回答中，我没有提到任何与并行化有关的问题，这是你后续行动的主题)。

使用Rcpp和data.table

用C++编译逻辑并使用data.table分组操作按组应用它，比您的基线快了大约2,000倍，远远超过您通过并行化所希望得到的速度。

在最初的示例中，有39,420,000行，它以1.883秒在我的计算机上执行；在有28,800行的修改后的机器上，这在0.004秒中执行。

代码语言：javascript
运行
复制
library(data.table)
library(Rcpp)

定义并编译一个C++函数，即在R脚本中内联的CalcSW()：

注意：C/C++中的计数从0__开始，而R__从1__开始--这就是为什么这里的指数是不同的。

代码语言：javascript
运行
复制
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

代码语言：javascript
运行
复制
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)

对CalcSW()和year的每个组合在df上执行函数loc.id，同时将返回的值赋值给三列：

代码语言：javascript
运行
复制
system.time({
  df[,  c("SW","PAW","aetc") := CalcSW(SW_ini,
                                       SW_max,
                                       rain,
                                       swc,
                                       PETc), keyby = .(loc.id, year)]
})

..。

代码语言：javascript
运行
复制
   user  system elapsed 
  0.004   0.000   0.004

结果：

代码语言：javascript
运行
复制
head(df)

..。

代码语言：javascript
运行
复制
   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

我不是百分之百肯定我完美地实现了你的逻辑，但是逻辑应该是非常简单的来调整我可能遗漏了什么东西，我实现它的方式非常类似于你的布局。

另一个注意事项:使用自动缩进和代码突出显示(无论是使用C++还是Emacs)编写RStudio要容易得多，如果您创建了一个单独的文件，名为TestCode.cpp，格式如下。

然后，您可以使用Rcpp::sourceCpp("TestCode.cpp")在R脚本中编译函数，或者像我前面所做的那样，复制并粘贴除前三行作为字符串的字符串之外的所有内容。

代码语言：javascript
运行
复制
#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);
}
票数 5
EN

Stack Overflow用户

发布于 2018-03-05 21:13:06

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

代码语言：javascript
运行
复制
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)
票数 1
EN

Stack Overflow用户

发布于 2018-05-21 20:27:31

马特的data.table插图非常好地说明了data.table可以有多快，因为它在没有拷贝和数据移动的情况下完成了计算。

然而，要回答关于使用管道的关键问题，您可以使用group_by和do来完成您所追求的目标(尽管比data.table慢得多)

下面我设置了马特做的同样的虚拟数据。然后，我使用您的函数(但是将大小写固定在PETc上)。虽然不快，但很容易跟上。

代码语言：javascript
运行
复制
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
票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：

https://stackoverflow.com/questions/49118364

复制
相关文章

相似问题

领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2025 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 深公网安备号 44030502008569

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号 | 京公网安备号11010802020287

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档