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 (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的强项，所以任何人都可以通过矢量化这个操作来协助我吗？

在此先感谢。

UPDATE

下面是来自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") 

UPDATE

谢谢大家对你的帮助。我仍然不熟悉矢量方式的做法，我无法及时实施一些更复杂的解决方案。我在Subs的建议下面有一些时间。现在快速完成我所需要的工作。这些数字都是以秒为超过P.

Z的一次迭代

我的方式：59.77

替补：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) 

希望这有助于:)

UPDATE（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中的列具有相同长度的整个新向量。为了整洁，我经常发现创建这些新的向量作为新的数据框列是很有用的，但显然这不是必须的。

第二个观察：使用“paste'n'match”策略将所需的Est.Date列从Z31正确传输到pred。这种任务有其他方法（例如使用合并），但是我采用这种方式，因为它完全是故障安全的并且保证在pred中保留顺序（这很关键）。实质上，粘贴操作只是让您一次匹配多个字段，因为如果粘贴的字符串匹配，则所有组成部分都匹配。我使用〜作为分隔符（前提是我知道该符号不会出现在任何字段中），以避免粘贴操作造成任何歧义。如果您使用空格分隔符，然后粘贴诸如（“AB”，“C”，“D”）之类的东西会得到与粘贴（“A”，“BC”，“D”）相同的结果 - 并且我们希望避免任何头痛！第三个观察：很容易对逻辑运算进行矢量化，例如查看一个矢量是否超过另一个矢量（请参阅pred $ ttValid），或者根据矢量的值选择一个或另一个值（请参阅pred $ ttFiltered）。在目前的情况下，这些可以合并为一条线，但是我更多地将事情分解成一个示范。

第四个观察：创建pred $ cum_tt_Z31的最后一行本质上只是在与$ pred $ LookupNum的每个单独值相对应的行上执行cumsum操作，使用tapply，它允许您在不同的行组上应用相同的函数（在这里，我们正在按$ predupNum进行分组）。 pred $ LookupNum的定义在这里有很大的帮助 - 它是一个数字索引，其中一个块为1，后面跟着一个2的块，依此类推。这意味着从tapply出来的cumsum载体的最终列表可以简单地不列出并放入向量中，并自动按照正确的顺序排列。如果你按照这样的顺序进行排序和分组，那么通常需要一些额外的代码行来正确地重新匹配事务（虽然它并不复杂）。最后的观察：如果最后一点太可怕了，值得强调的是，如果循环中的工作被很好地矢量化，那么在少数情况下（比如说48）的快速循环并不总是灾难性的。 UPDATE部分末尾的“替代方法”表明，cumsum-on-group步骤也可以通过预先准备答案栏（最初是所有的NA），然后逐个浏览48个子集并填写每个块都有相应的cumsum。正如文中所指出的那样，这一步骤的速度大约是巧妙的方法的一半，如果需要更多的子集，这将会相当慢。

如果有人对此类任务有任何后续问题，请不要犹豫，给我留言。

Answer 2

一个快速的解决方案是定义一个向量：使用数字日期14.62

替补循环作为外：

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起步价1的j指数时要小心。在你的例子中，它不会导致这个条件问题。

编辑：

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

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

2）根据日期 - 矢量排序的Z31 data_frame：Z31=Z31[with(Z31, order(-Date)), ]（注以降序因为要比较PRED [，日期索引]> = Z31 [，日期索引]

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