r - R - 需要帮助加快 for 循环

Question

我有两个数据框；一个是 48 行长，如下所示：

名称 = Z31

     Est.Date   Site    Cultivar   Planting
1   24/07/2011  Birchip Axe           1
2   08/08/2011  Birchip Bolac         1
3   24/07/2011  Birchip Derrimut      1
4   12/08/2011  Birchip Eaglehawk     1
5   29/07/2011  Birchip Gregory       1
6   29/07/2011  Birchip Lincoln       1
7   23/07/2011  Birchip Mace          1
8   29/07/2011  Birchip Scout         1
9   17/09/2011  Birchip Axe           2
10  19/09/2011  Birchip Bolac         2

另一个是 > 23000 行，包含来自模拟器的输出。它看起来像这样：

名称 = pred

    Date        maxt    mint    Cultivar    Site    Planting    tt  cum_tt
1   5/05/2011   18       6.5    Axe        Birchip  1        12.25  0
2   6/05/2011   17.5     2.5    Axe        Birchip  1        10     0
3   7/05/2011   18       2.5    Axe        Birchip  1        10.25  0
4   8/05/2011   19.5       2    Axe        Birchip  1        10.75  0
5   9/05/2011   17       4.5    Axe        Birchip  1        10.75  0
6   10/05/2011  15.5    -0.5    Axe        Birchip  1        7.5    0
7   11/05/2011  14       5.5    Axe        Birchip  1        9.75   0
8   12/05/2011  19         8    Axe        Birchip  1        13.5   0
9   13/05/2011  18.5     7.5    Axe        Birchip  1        13     0
10  14/05/2011  16       3.5    Axe        Birchip  1        9.75   0

我想要做的是让 cum_tt 列开始将当前行的 tt 列添加到前一行的 cum_tt （累积加法），前提是 pred DF 中的日期等于或大于 Z31 est.Date . 我编写了以下 for 循环：

for (i in 1:nrow(Z31)){
  for (j in 1:nrow(pred)){
    if (Z31[i,]$Site == pred[j,]$Site & Z31[i,]$Cultivar == pred[j,]$Cultivar & Z31[i,]$Planting == pred[j,]$Planting &
        pred[j,]$Date >= Z31[i,]$Est.Date)
    {
      pred[j,]$cum_tt <- pred[j,]$tt + pred[j-1,]$cum_tt
    }
  }
}

这行得通，但它太慢了，运行整个集合大约需要一个小时。我知道循环不是 R 的强项，所以任何人都可以帮助我对这个操作进行矢量化吗？

提前致谢。

更新

这是 dput(Z31) 的输出：

structure(list(Est.Date = structure(c(15179, 15194, 15179, 15198, 15184, 15184, 15178, 15184, 15234, 15236, 15230, 15238, 15229, 15236, 15229, 15231, 15155, 15170, 15160, 15168, 15165, 15159, 15170, 15170, 15191, 15205, 15198, 15203, 15202, 15195, 15203, 15206, 15193, 15193, 15195, 15200, 15193, 15205, 15200, 15205, 15226, 15245, 15231, 15259, 15241, 15241, 15241, 15241), class = "Date"), Site = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("Birchip", "Gatton", "Tarlee"), class = "factor"), Cultivar = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L), .Label = c("Axe", "Bolac", "Derrimut", "Eaglehawk", "Gregory", "Lincoln", "Mace", "Scout"), class = "factor"), Planting = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L)), .Names = c("Est.Date", "Site", "Cultivar", "Planting"), row.names = c(NA, -48L), class = "data.frame")

这里是预告。请注意，这里有一些额外的列。为了便于阅读，我刚刚包括了上面的相关内容。

structure(list(Date = structure(c(15099, 15100, 15101, 15102, 
15103, 15104, 15105, 15106, 15107, 15108, 15109, 15110, 15111, 
15112, 15113, 15114, 15115, 15116, 15117, 15118), class = "Date"), 
    flowering_das = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), Zadok = c(9, 9, 
    9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 11, 11.032, 11.085, 
    11.157), stagename = structure(c(8L, 8L, 8L, 8L, 8L, 8L, 
    8L, 8L, 9L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 1L, 3L, 3L, 3L), .Label = c("emergence", 
    "end_grain_fill", "end_of_juvenil", "floral_initiat", "flowering", 
    "germination", "maturity", "out", "sowing", "start_grain_fi"
    ), class = "factor"), node_no = c(0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 2, 2.032, 2.085, 2.157), maxt = c(18, 
    17.5, 18, 19.5, 17, 15.5, 14, 19, 18.5, 16, 16, 15, 16.5, 
    16.5, 20.5, 23, 25.5, 16.5, 16.5, 15), mint = c(6.5, 2.5, 
    2.5, 2, 4.5, -0.5, 5.5, 8, 7.5, 3.5, 6, 1, 5.5, 2, 7, 7, 
    9, 13.5, 11.5, 8.5), Cultivar = c("Axe", "Axe", "Axe", "Axe", 
    "Axe", "Axe", "Axe", "Axe", "Axe", "Axe", "Axe", "Axe", "Axe", 
    "Axe", "Axe", "Axe", "Axe", "Axe", "Axe", "Axe"), Site = c("Birchip", 
    "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", 
    "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", 
    "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", "Birchip", 
    "Birchip"), Planting = c("1", "1", "1", "1", "1", "1", "1", 
    "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", 
    "1"), `NA` = c("Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", 
    "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out", "Birchip TOS1 Axe.out"
    ), tt = c(12.25, 10, 10.25, 10.75, 10.75, 7.5, 9.75, 13.5, 
    13, 9.75, 11, 8, 11, 9.25, 13.75, 15, 17.25, 15, 14, 11.75
    ), cum_tt = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0)), .Names = c("Date", "flowering_das", "Zadok", 
"stagename", "node_no", "maxt", "mint", "Cultivar", "Site", "Planting", 
NA, "tt", "cum_tt"), row.names = c(NA, 20L), class = "data.frame")

更新

感谢大家的帮助。我对矢量处理方式仍然很陌生，我无法及时实施一些更复杂的解决方案。对于 Subs 建议的方式，我在下面有一些时间安排。它现在足够快，可以做我需要的。对于 Z 对 P 的一次迭代，这些数字以秒为单位。

我的方式：59.77

替补：14.62

使用数字日期的潜艇：11.12

Answer 1

当然，我们可以在几秒钟内做到这一点......我在这里的第一个答案，所以温柔点！

## first make sure we have dates in a suitable format for comparison
## by using strptime, creating the columns estdate_tidy and date_tidy
## in Z31 and pred respectively

Z31$estdate_tidy = strptime(as.character(Z31$Est.Date), "%d/%m/%Y")
pred$date_tidy = strptime(as.character(pred$Date), "%d/%m/%Y")

## now map the estdate_tidy column over to pred using the match command -
## Z31_m and pred_m are dummy variables that hopefully make this clear

Z31_m = paste(Z31$Site, Z31$Cultivar, Z31$Planting)
pred_m = paste(pred$Site, pred$Cultivar, pred$Planting)
pred$estdate_tidy = Z31$estdate_tidy[match(pred_m, Z31_m)]

## then define a ttfilter column that copies tt, except for being 0 when
## estdate_tidy is after date_tidy (think this is what you described)

pred$ttfilter = ifelse(pred$date_tidy >= pred$estdate_tidy, pred$tt, 0)

## finally use cumsum function to sum this new column up (looks like you
## wanted the answer in cum_tt so I've put it there)

pred$cum_tt = cumsum(pred$ttfilter)

希望这可以帮助 ：）

更新（6 月 7 日）：

处理新规范的矢量化代码 - 即应针对每组条件（场地/栽培品种/种植）分别完成 cumsum - 如下所示：

Z31$Lookup=with(Z31,paste(Site,Cultivar,Planting,sep="~"))
Z31$LookupNum=match(Z31$Lookup,unique(Z31$Lookup))
pred$Lookup=with(pred,paste(Site,Cultivar,Planting,sep="~"))
pred$LookupNum=match(pred$Lookup,unique(pred$Lookup))

pred$Est.Date = Z31$Est.Date[match(pred$Lookup,Z31$Lookup)]
pred$ttValid = (pred$Date>=pred$Est.Date)
pred$ttFiltered = ifelse(pred$ttValid, pred$tt, 0)

### now fill in cumsum of ttFiltered separately for each LookupNum
pred$cum_tt_Z31 = as.vector(unlist(tapply(pred$ttFiltered,
                                          pred$LookupNum,cumsum)))

在我的机器上运行时间为 0.16 秒，最后pred$cum_tt_Z31一列与非矢量化代码的答案完全匹配:)

为了完整起见，值得注意的是，上面最后一个复杂的 tapply 行可以用以下更简单的方法替换，在 48 种可能的情况下进行短循环：

pred$cum_tt_Z31 = rep(NA, nrow(pred))
for (lookup in unique(pred$Lookup)) {
    subs = which(pred$Lookup==lookup)
    pred$cum_tt_Z31[subs] = cumsum(pred$ttFiltered[subs])
}

运行时间仅略微增加到 0.25 秒左右，因为这里的循环非常小，并且在循环内完成的工作是矢量化的。

认为我们已经破解了它！:)

关于矢量化的一些快速观察（6 月 8 日）：

将过程步骤矢量化的过程使运行时间从接近一个小时减少到总共 0.16 秒。即使考虑到不同的机器速度，这也是一个至少 10,000 倍的加速，这使一个人通过进行细微调整但仍保留循环结构可能获得的 2-5 倍相形见绌。

要进行的第一个关键观察：在解决方案中，每一行都创建 - 没有循环 - 一个全新的向量，其长度与 Z31 或 pred 中的列相同。为了简洁起见，我经常发现将这些新向量创建为新的数据框列很有用，但显然这并不是绝对必要的。

第二个观察结果：所需的 Est.Date 列使用“paste'n'match”策略正确地从 Z31 转移到 pred。这类任务有其他方法（例如使用合并），但我采用这条路线，因为它是完全故障安全的，并保证保留 pred 中的顺序（这是至关重要的）。本质上，粘贴操作只是让您一次匹配多个字段，因为如果粘贴的字符串匹配，那么它们的所有组成部分都匹配。我使用 ~ 作为分隔符（前提是我知道该符号不会出现在任何字段中）以避免粘贴操作导致任何歧义。如果您使用空格分隔符，则将 ("A B", "C", "D") 之类的内容粘贴在一起会得到与粘贴 ("A", "B C", "D" 相同的结果

第三个观察：向量化逻辑操作很容易，例如查看一个向量是否超过另一个（参见 pred$ttValid），或根据向量的值选择一个或另一个值（参见 pred$ttFiltered）。在目前的情况下，这些可以组合成一条线，但我把事情分解了一点作为演示。

第四个观察：创建 pred$cum_tt_Z31 的最后一行本质上只是使用 tapply 对对应于 pred$LookupNum 的每个单独值的行执行 cumsum 操作，它允许您对不同的行组应用相同的函数（这里，我们'按 pred$LookupNum 分组）。pred$LookupNum 的定义在这里有很大帮助——它是一个数字索引，包含一个 1 的块，然后是一个 2 的块，依此类推。这意味着来自 tapply 的 cumsum 向量的最终列表可以简单地不列出并放入向量中，并自动按正确的顺序排列。如果您执行 tapply 并按未像这样排序的组进行拆分，您通常需要几行额外的代码才能正确地重新匹配事物（尽管这并不棘手）。

最后的观察：如果最后的 tapply 太可怕了，那么值得强调的是，如果循环内的工作被很好地矢量化，那么在少数情况下（比如 48 个）的快速循环并不总是灾难性的。UPDATE 部分末尾的“替代方法”表明，cumsum-on-groups 步骤也可以通过预先准备一个答案列（最初是所有 NA）然后将 48 个子集一一遍历并填写每个块都有适当的 cumsum。然而，正如文中所指出的，这一步的速度大约是聪明的 tapply 方法的一半，如果需要更多的子集，这将相当慢。

如果有人对此类任务有任何后续问题，请随时给我留言。

Answer 2

一个快速的解决方案是将循环外的向量定义为：

 temp_cumtt=c(rep(0,nrow(pred)))

然后使用这个：

if (Z31[i,2] == pred[j,5] & Z31[i,3] == pred[j,4] & Z31[i,4] == pred[j,6] & pred[j,1] >= Z31[i,1]){
   temp_cumtt[j]=pred[j,7] + pred[j-1,8]
}

而不是直接更新 data.frame 列。

退出循环后，您可以更新该列：

 pred$cum_tt = temp_cumtt  

j-1另一件事是在使用with index of jstarting at时必须小心1。在您的示例中，它不会导致该条件问题。

编辑：

现在查看您的数据格式，我有这些建议。

1）不要转换为Date类，而是将其保留为值向量。

2）Z31根据Date-vector对data_frame进行排序：（Z31=Z31[with(Z31, order(-Date)), ]注意按降序排列，因为要比较pred[,Date index]>=Z31[,Date index]

3）使用第一个循环作为pred。首先取 pred -> 的日期pred[i,1]并尝试二进制排序并找到它满足 Z31 中的哪个索引，然后从该索引开始向下列表。如果Date条件满足，则检查其余条件并temp_cumtt[i]像以前一样填写。

这应该非常快（因为二进制排序仅在 48 行上Z31，您可以将运行时间与其他解决方案进行比较。

Answer 3

让我们使用data.table，它应该会加快速度。

Z31 <- data.table(Z31,key="Site,Cultivar,Planting")
pred <- data.table(pred)

## First, let's create an extra column in `pred` to see the corresponding date from `Z31` 
## Note 1: The JT is necessary since both sets have the same column names
## Note 2: I needed to use as.integer on Planting to make it work

pred[,Z31Est.Date:={JT=J(Site,Cultivar,as.integer(Planting)); Z31[JT,Est.Date][[4]]}]

## Now we can see for each row whether the date in `pred` is higher than or equal to that from `Z31`.

pred[,DateTrue:=Date>=Z31Est.Date]

## Finally, we only have to add up `pred[i,tt]` and `pred[i-1,cum_tt]` for each row where `DateTrue` equals `TRUE`.

for (i in 1:nrow(pred)) set(pred,i,13L,if(pred[i,DateTrue]) pred[i-1,cum_tt]+pred[i,tt] else(0))