我希望能够在不进行完全交叉连接的情况下，使用范围标准来交叉连接两个数据表。实质上，我想要过滤/范围表达式的 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
结论
- Arun提供的Bioconductor树/IRanges解决方案（bioTreeRangeRes）比其他解决方案快两个数量级。但安装似乎更新了其他CRAN库（是我的问题，当安装程序询问时我接受了它），其中一些在加载时无法找到，例如
gtools和gplots。 - BrodieG提供的最快纯数据表选项（dtJoin2）可能不如我需要的那么高效，但至少在内存消耗方面是合理的（我将让它在我的实际场景中运行过夜，大约有100万行）。
- 我尝试改变数据表键（使用日期而不是ID），但没有任何影响。
- 正如预期的那样，在R中显式编写循环（manualIter）非常慢。
dt2为D可能落入的每个区间制作标签：[5,7]是“21”；[7,9]是“21,22”等，具有适当的边缘情况条件。之后，只需将这些标签应用于dt1即可。 - Frank