R - 需要帮助加快 for 循环答案

【问题标题】：R - Need help speeding up a for loopR - 需要帮助加快 for 循环
【发布时间】：2012-06-10 01:34:57
【问题描述】：

我有两个数据框；一个是 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 行，包含来自模拟器的输出。它看起来像这样：

name = 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

我想要做的是只有当 pred DF 中的日期等于或大于 Z31估计日期。我写了以下 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

【问题讨论】：

您能否用dput(Z31) 和dput(pred) 的结果更新您的问题？在阅读...中的文本后，我无法让您的代码按原样运行
仅供参考，这看起来像 data.table() 和 roll = TRUE 会更快。
如上关于 dput 输出。看看 data.table ......这可能是午餐后要探索的东西:)
dput(head(pred,20)) 应该可以解决问题 :)
bonk 当然...给你 :) 现在去吃饭。

标签： r for-loop dataframe vectorization data.table

【解决方案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 中的顺序（这是至关重要的）。本质上，粘贴操作只是让您一次匹配多个字段，因为如果粘贴的字符串匹配，那么它们的所有组成部分都匹配。我使用 ~ 作为分隔符（前提是我知道该符号不会出现在任何字段中）以避免粘贴操作导致任何歧义。如果您使用空格分隔符，则将 ("AB", "C", "D") 之类的内容粘贴在一起会得到与粘贴 ("A", "BC", "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 方法的一半，如果需要更多的子集，这将相当慢。

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

【讨论】：

看起来很有希望。 match 粘贴向量的策略在时间上看起来是 O(n)，而使用单独访问数据帧的循环策略是 O(n^2)。
看起来 Est.Date 已经是日期格式（原始帖子没有说清楚），所以在我上面的解决方案中 Z31$estdate_tidy 可以替换为 Z31$Est.Date （类似于 pred $date_tidy 和 pred$Date，如果 pred$Date 已经是日期而不是字符格式）:)
感谢 DWin - 是的，预计匹配策略在这里是最佳的。总体上应该实现至少 1000 倍的加速。
希望我有足够的分数来为下面的讨论做出贡献……但我认为循环在这里是一个糟糕的策略！很高兴看到一些时间比较......
我确实简要地查看了合并，但不确定如何保证它正确保留顺序，这在这里很重要。我最初的几次尝试似乎破坏了订单，所以我转而采用了一种保证有效的方法。什么合并语法在这里可以正常工作？从道德上讲，这两种方法无论如何都是一条线（Z31_m 和 pred_m 是为了清楚起见而单独定义的）。

【解决方案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 开始的j-1 时必须小心。在您的示例中，它不会导致该条件问题。

编辑：

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

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

2) 根据日期向量对Z31 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 进行，您可以将运行时间与其他解决方案进行比较。

【讨论】：

感谢 subs，这确实提高了执行速度，每次外部迭代我从 60 秒下降到大约 15 秒。我希望能够在没有循环的情况下做到这一点，因为它们'太慢了。不过，这是一个很大的进步。我不担心 j-1 索引，因为在数据集中我使用的第一行永远不会是真的。
循环并不是低效的。这取决于你在里面做什么。另一种方法是根据您的Date 对您的Z31 数据框进行排序，并且由于pred 中的Date 已经排序，您可以检查使用算法从索引开始比较，而不是遍历所有Z 行 * P 行
再次感谢订阅者。我需要继续我正在做的事情，所以我无法实现二进制排序技术。不过，您帮助我将其加快到可接受的水平。上述问题中的一些时间安排。

【解决方案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))

【讨论】：

仍然包含 for (i in 1:nrow(pred)) 的“解决方案”并没有真正的帮助 - 我们已经有了一个完全矢量化的解决方案。
问题是加速循环，对吗？好吧，我认为这可以加快速度，但我无法测试它是否确实如此。但确实应该避免使用循环。
我认为这个问题是......“我知道循环不是 R 的强项，所以任何人都可以帮助我对这个操作进行矢量化吗？” :)
@Dirk 我可以看到你在这里尝试了什么。不过，我怀疑更短更简单的data.table 解决方案是可能的。甚至可能是一个或两个衬垫，它应该比 Tim 使用 paste 的矢量化更快（这通常是不能扩展到 1e6、1e7 或 1e9 行的部分）。有一天我会回到这个。我们可以使用1e7 示例来演示时间安排。