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 行，包含来自模拟器的输出。它看起来像这样:

名称 = 预

    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")
这是pred。注意这里有一些额外的列。为了便于阅读，我刚刚包括了上面的相关内容。

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 应针对每组条件( field /栽培品种/种植)单独完成 - 如下所示:

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 中的顺序(这很关键)。本质上，粘贴操作只是让您一次匹配多个字段，因为如果粘贴的字符串匹配，则它们的所有组成部分都匹配。我使用 ~ 作为分隔符(前提是我知道该符号不会出现在任何字段中)以避免粘贴操作导致任何歧义。如果您使用空格分隔符，那么将诸如 ("AB", "C", "D") 之类的内容粘贴在一起会产生与粘贴 ("A", "BC", "D") 相同的结果 - 我们想要避免任何头痛!

第三个观察:很容易向量化逻辑操作，例如查看一个向量是否超过另一个(参见 pred$ttValid)，或者根据向量的值选择一个或另一个值(参见 pred$ttFiltered)。在目前的情况下，这些可以合并为一行，但作为演示，我将其分解了一点。

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

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

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