遠程/過濾交叉聯接有R data.table

Question

我想交叉聯接兩個數據表沒有評價全交叉連接，在使用過程中的測距標準。實質上，我希望CJ具有過濾/範圍表達。

有人能提出一個高績效的方法避免了全交叉連接？

看下面的測試例子，用邪惡的完成交叉連接。

library(data.table) 

# Test data. 
dt1 <- data.table(id1=1:10, D=2*(1:10), key="id1") 
dt2 <- data.table(id2=21:23, D1=c(5, 7, 10), D2=c(9, 12, 16), key="id2") 

# Desired filtered cross-join data table by hand: D1 <= D & D <= D2. 
dtfDesired <- data.table(
    id1=c(3, 4, 4, 5, 6, 5, 6, 7, 8) 
    , id2=c(rep(21, 2), rep(22, 3), rep(23, 4)) 
    , D1=c(rep(5, 2), rep(7, 3), rep(10, 4)) 
    , D=c(6, 8, 8, 10, 12, 10, 12, 14, 16) 
    , D2=c(rep(9, 2), rep(12, 3), rep(16, 4)) 
) 
setkey(dtfDesired, id1, id2) 

# My "inefficient" programmatic attempt with full cross join. 
fullCJ <- function(dt1, dt2) { 
    # Full cross-product: NOT acceptable with real data! 
    dtCrossAll <- CJ(dt1$id1, dt2$id2) 
    setnames(dtCrossAll, c("id1", "id2")) 
    # Merge all columns. 
    dtf <- merge(dtCrossAll, dt1, by="id1") 
    dtf <- merge(dtf, dt2, by="id2") 
    setkey(dtf, id1, id2) 
    # Reorder columns for convenience. 
    setcolorder(dtf, c("id1", "id2", "D1", "D", "D2")) 
    # Finally, filter the cases I want. 
    dtf[D1 <= D & D <= D2, ] 
} 

dtf <- fullCJ(dt1, dt2) 

# Print results. 
print(dt1) 
print(dt2) 
print(dtfDesired) 
all.equal(dtf, dtfDesired) 

測試數據輸出

> # Print results. 
> print(dt1) 
    id1 D 
1: 1 2 
2: 2 4 
3: 3 6 
4: 4 8 
5: 5 10 
6: 6 12 
7: 7 14 
8: 8 16 
9: 9 18 
10: 10 20 
> print(dt2) 
    id2 D1 D2 
1: 21 5 9 
2: 22 7 12 
3: 23 10 16 
> print(dtfDesired) 
    id1 id2 D1 D D2 
1: 3 21 5 6 9 
2: 4 21 5 8 9 
3: 4 22 7 8 12 
4: 5 22 7 10 12 
5: 5 23 10 10 16 
6: 6 22 7 12 12 
7: 6 23 10 12 16 
8: 7 23 10 14 16 
9: 8 23 10 16 16 
> all.equal(dtf, dtfDesired) 
[1] TRUE 

所以，現在的挑戰是編寫過濾交叉在可擴展到數百萬行的方式加入！

下面是替代實現包括那些答案和意見提出的集合。

# My "inefficient" programmatic attempt looping manually. 
manualIter <- function(dt1, dt2) { 
    id1Match <- NULL; id2Match <- NULL; dtf <- NULL; 
    for (i1 in seq_len(nrow(dt1))) { 
    # Find matches in dt2 of this dt1 row. 
    row1 <- dt1[i1, ] 
    id1 <- row1$id1 
    D <- row1$D 
    dt2Match <- dt2[D1 <= D & D <= D2, ] 
    nMatches <- nrow(dt2Match) 
    if (0 < nMatches) { 
     id1Match <- c(id1Match, rep(id1, nMatches)) 
     id2Match <- c(id2Match, dt2Match$id2) 
    } 
    } 
    # Build the return data.table for the matching ids. 
    dtf <- data.table(id1=id1Match, id2=id2Match) 
    dtf <- merge(dtf, dt1, by="id1") 
    dtf <- merge(dtf, dt2, by="id2") 
    setkey(dtf, id1, id2) 
    # Reorder columns for convenience & consistency. 
    setcolorder(dtf, c("id1", "id2", "D1", "D", "D2")) 
    return(dtf) 
} 

dtJoin1 <- function(dt1, dt2) { 
    dtf <- dt1[, dt2[D1 <= D & D <= D2, list(id2=id2)], by=id1] 
    dtf <- merge(dtf, dt1, by="id1") 
    dtf <- merge(dtf, dt2, by="id2") 
    setkey(dtf, id1, id2) 
    setcolorder(dtf, c("id1", "id2", "D1", "D", "D2")) # Reorder columns for convenience & consistency. 
    return(dtf) 
} 

dtJoin2 <- function(dt1, dt2) { 
    dtf <- dt2[, dt1[D1 <= D & D <= D2, list(id1=id1, D1=D1, D=D, D2=D2)], by=id2] 
    setkey(dtf, id1, id2) 
    setcolorder(dtf, c("id1", "id2", "D1", "D", "D2")) # Reorder columns for convenience & consistency. 
    return(dtf) 
} 

# Install Bioconductor IRanges (see bioTreeRange below). 
source("http://bioconductor.org/biocLite.R") 
biocLite("IRanges") 

# Solution using Bioconductor IRanges. 
bioTreeRange <- function(dt1, dt2) { 
    require(IRanges) 
    ir1 <- IRanges(dt1$D, width=1L) 
    ir2 <- IRanges(dt2$D1, dt2$D2) 
    olaps <- findOverlaps(ir1, ir2, type="within") 
    dtf <- cbind(dt1[queryHits(olaps)], dt2[subjectHits(olaps)]) 
    setkey(dtf, id1, id2) 
    setcolorder(dtf, c("id1", "id2", "D1", "D", "D2")) # Reorder columns for convenience. 
    return(dtf) 
} 

而且現在下面是一個更大的數據有點基準設置量值比我真正的根本情景小2-3個數量級。真正的場景在完全交叉連接巨大的內存分配上失敗。

set.seed(1) 
n1 <- 10000 
n2 <- 1000 
dtbig1 <- data.table(id1=1:n1, D=1:n1, key="id1") 
dtbig2 <- data.table(id2=1:n2, D1=sort(sample(1:n1, n2)), key="id2") 
dtbig2$D2 <- with(dtbig2, D1 + 100) 

library("microbenchmark") 
mbenchmarkRes <- microbenchmark(
    fullCJRes <- fullCJ(dtbig1, dtbig2) 
    , manualIterRes <- manualIter(dtbig1, dtbig2) 
    , dtJoin1Res <- dtJoin1(dtbig1, dtbig2) 
    , dtJoin2Res <- dtJoin2(dtbig1, dtbig2) 
    , bioTreeRangeRes <- bioTreeRange(dtbig1, dtbig2) 
    , times=3, unit="s", control=list(order="inorder", warmup=1) 
) 
mbenchmarkRes$expr <- c("fullCJ", "manualIter", "dtJoin1", "dtJoin2", "bioTreeRangeRes") # Shorten names for better display. 

# Print microbenchmark 
print(mbenchmarkRes, order="median") 

現在目前的基準測試結果我在我的機器上得到：

> print(mbenchmarkRes, order="median") 
Unit: seconds 
      expr  min   lq  median   uq  max neval 
bioTreeRangeRes 0.05833279 0.05843753 0.05854227 0.06099377 0.06344527  3 
     dtJoin2 1.20519664 1.21583650 1.22647637 1.23606216 1.24564796  3 
      fullCJ 4.00370434 4.03572702 4.06774969 4.17001658 4.27228347  3 
     dtJoin1 8.02416333 8.03504136 8.04591938 8.20015977 8.35440016  3 
     manualIter 8.69061759 8.69716448 8.70371137 8.76859060 8.83346982  3 

結論

• 從阿倫的Bioconductor的樹/ IRanges溶液（bioTreeRangeRes）是兩個數量級快比替代品。但安裝似乎已經更新了其他的CRAN庫（我的錯，我在安裝時接受了這個問題）;當加載它們時，它們中的一些不能再被找到 - 例如，gtools和gplots。
• 從BrodieG（dtJoin2）最快的純data.table選擇可能是效率不高，因爲我需要它，但至少是合理的內存消耗方面（我將讓它在我的真實的情景〜100萬在夜間運行行）。
• 我試着改變數據表鍵（使用日期而不是id）;它沒有任何影響。
• 正如預期的那樣，在R（manualIter）中顯式編寫循環將進行爬網。

Answer 1

近日，重疊聯接在data.table實施。這是dt1'開始和結束點相同的特殊情況。你可以抓住從GitHub的項目頁面的最新版本嘗試了這一點：

require(data.table) ## 1.9.3+ 
dt1[, DD := D] ## duplicate column D to create intervals 
setkey(dt2, D1,D2) ## key needs to be set for 2nd argument 
foverlaps(dt1, dt2, by.x=c("D", "DD"), by.y=key(dt2), nomatch=0L) 

# id2 D1 D2 id1 D DD 
# 1: 21 5 9 3 6 6 
# 2: 21 5 9 4 8 8 
# 3: 22 7 12 4 8 8 
# 4: 22 7 12 5 10 10 
# 5: 23 10 16 5 10 10 
# 6: 22 7 12 6 12 12 
# 7: 23 10 16 6 12 12 
# 8: 23 10 16 7 14 14 
# 9: 23 10 16 8 16 16 

這裏是在相同的數據結果基準您在您的文章已經表明：

# Unit: seconds 
#    expr   min   lq  median   uq   max neval 
#   olaps 0.03600603 0.03971068 0.04341533 0.04857602 0.05373671  3 
# bioTreeRangeRes 0.11356837 0.11673968 0.11991100 0.12499391 0.13007681  3 
#   dtJoin2 2.61679908 2.70327940 2.78975971 2.86864832 2.94753693  3 
#   fullCJ 4.45173294 4.75271285 5.05369275 5.08333291 5.11297307  3 
#   dtJoin1 16.51898878 17.39207632 18.26516387 18.60092303 18.93668220  3 
#  manualIter 29.36023340 30.13354967 30.90686594 33.55910653 36.21134712  3 

其中dt_olaps是：

dt_olaps <- function(dt1, dt2) { 
    dt1[, DD := D] 
    setkey(dt2, D1,D2) 
    foverlaps(dt1, dt2, by.x=c("D","DD"), by.y=key(dt2), nomatch=0L) 
} 

Answer 2

這似乎是一個問題，可以從使用interval trees算法中受益很多。 bioconductor包IRanges提供了一個非常好的實現。

# Installation 
source("http://bioconductor.org/biocLite.R") 
biocLite("IRanges") 

# solution 
require(IRanges) 
ir1 <- IRanges(dt1$D, width=1L) 
ir2 <- IRanges(dt2$D1, dt2$D2) 

olaps <- findOverlaps(ir1, ir2, type="within") 
cbind(dt1[queryHits(olaps)], dt2[subjectHits(olaps)]) 

    id1 D id2 D1 D2 
1: 3 6 21 5 9 
2: 4 8 21 5 9 
3: 4 8 22 7 12 
4: 5 10 22 7 12 
5: 5 10 23 10 16 
6: 6 12 22 7 12 
7: 6 12 23 10 16 
8: 7 14 23 10 16 
9: 8 16 23 10 16