r - R - for ループの高速化に助けが必要

Question

2 つのデータフレームがあります。1 つは 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

もう 1 つは > 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

私がしたいのは、pred DFの日付がZ31 est.Date以上の場合にのみ、現在の行のtt列を前の行のcum_ttに追加するためにcum_tt列を開始することです(累積加算) . 次の 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
    }
  }
}

これは機能しますが、非常に遅いため、セット全体を実行するには約 1 時間かかります。ループが 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が提案した方法について、以下にタイミングを示します。私が必要とすることをするのに十分な速さです。これらの数値は、P に対する Z の 1 回の反復の秒数です。

私のやり方: 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 日):

新しい仕様に取り組むためのベクトル化されたコード、つまり、一連の条件 (サイト/品種/植栽) ごとに個別に累積計算を行う必要がある - を以下に示します。

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列はベクトル化されていないコードからの回答と正確に一致します:)

完全を期すために、上記の最後の複雑なタップ行は、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 日):

プロセスのステップをベクトル化するプロセスにより、実行時間が 1 時間近くから合計で 0.16 秒に短縮されました。さまざまなマシンの速度を考慮しても、これは少なくとも 10,000 倍の速度アップであり、ループ構造を維持しながら微調整を行うことで得られる 2 ～ 5 倍の速度を小さくします。

最初の重要な観察: ソリューションでは、すべての行が、ループなしで、Z31 または pred の列と同じ長さのまったく新しいベクトルを作成します。きれいにするために、これらの新しいベクトルを新しいデータ フレーム列として作成すると便利なことがよくありますが、これは厳密には必要ではありません。

2 番目の観察: 必要な Est.Date 列は、「貼り付けと一致」戦略を使用して、Z31 から pred に正しく転送されます。この種のタスクには別のアプローチがあります (例: マージの使用) が、完全にフェイルセーフであり、pred での順序の保持が保証されているため (これは重要です)、このルートを使用します。基本的に、貼り付け操作では一度に複数のフィールドを一致させることができます。貼り付けられた文字列が一致すると、それらのすべての構成部分が一致するためです。~ をセパレーターとして使用し (シンボルがどのフィールドにも表示されないことがわかっている場合)、貼り付け操作があいまいになるのを回避します。スペース区切りを使用する場合、("A B", "C", "D") のように貼り付けると、("A", "B C", "D") を貼り付けた場合と同じ結果になります。

3 番目の観察: あるベクトルが別のベクトルを超えているかどうかを調べたり (pred$ttValid を参照)、ベクトルの値に基づいて値を選択したり (pred$ttFiltered を参照) するなど、論理演算をベクトル化するのは簡単です。現状では、これらを 1 行にまとめることができますが、デモンストレーションとして、もう少し分解してみました。

4 番目の観察: pred$cum_tt_Z31 を作成する最後の行は、本質的に、異なる行グループに同じ関数を適用できるようにする、tapply を使用して、pred$LookupNum のそれぞれの個別の値に対応する行全体で cumsum 操作を行うだけです (ここでは、 pred$LookupNum でグループ化しています)。pred$LookupNum の定義は、ここで非常に役立ちます。これは、1 のブロックの後に 2 のブロックが続く数値インデックスです。これは、tapply から出力される cumsum ベクトルの最終的なリストを単純に非表示にしてベクトルに入れることができ、自動的に正しい順序になることを意味します。タップして、このように順序付けされていないグループで分割する場合、通常、物事を正しく元に戻すために数行のコードを追加する必要があります (トリッキーではありませんが)。

最終的な観察: 最後のタップが怖すぎる場合は、ループ内の作業が適切にベクトル化されていれば、少数のケース (48 と言う) のクイック ループが必ずしも悲惨ではないことを強調する価値があります。UPDATE セクションの最後にある「代替方法」は、回答列 (最初はすべての NA) を事前に準備し、次に 48 のサブセットを 1 つずつ調べて入力することによっても、グループに対する累積ステップを実行できることを示しています。適切な cumsum を持つ各ブロック。ただし、テキストで述べたように、この 1 つのステップは巧妙なタップ アプローチの約半分の速度であり、より多くのサブセットが必要な場合、これはかなり遅くなります。

この種のタスクについてフォローアップの質問がある場合は、遠慮なく私に声をかけてください。

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  

もう 1 つは、から始まるj-1インデックスで使用する場合は注意が必要です。あなたの例では、その条件付きの問題にはつながりません。j1

編集：

データ形式を見て、これらの提案があります。

1) クラスに変換せずDate、代わりに値のベクトルとして保持します。

2) Z31Date-vector に従って data_frame をソートします: Z31=Z31[with(Z31, order(-Date)), ](pred[,Date index]>=Z31[,Date index] を比較するため、降順で注意してください。

3) として 1 番目のループを使用しますpred。最初に pred の Date を取得 ->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))