Есть ли способ ускорить команду combn, чтобы получить все уникальные комбинации из 2 элементов, взятых из вектора?
Обычно это можно настроить следующим образом:
# Get latest version of data.table
library(devtools)
install_github("Rdatatable/data.table", build_vignettes = FALSE)
library(data.table)
# Toy data
d <- data.table(id=as.character(paste0("A", 10001:15000)))
# Transform data
system.time({
d.1 <- as.data.table(t(combn(d$id, 2)))
})
Тем не менее, combn в 10 раз медленнее (23 секунды против 3 секунд на моем компьютере), чем вычисление всех возможных комбинаций с использованием data.table.
system.time({
d.2 <- d[, list(neighbor=d$id[-which(d$id==id)]), by=c("id")]
})
Работа с очень большими векторами, я ищу способ сохранить память, только вычисляя уникальные комбинации (например, combn), но со скоростью data.table(см. второй фрагмент кода).
Я ценю любую помощь.
Вы можете использовать combnPrim из gRbase
source("http://bioconductor.org/biocLite.R")
biocLite("gRbase") # will install dependent packages automatically.
system.time({
d.1 <- as.data.table(t(combn(d$id, 2)))
})
# user system elapsed
# 27.322 0.585 27.674
system.time({
d.2 <- as.data.table(t(combnPrim(d$id,2)))
})
# user system elapsed
# 2.317 0.110 2.425
identical(d.1[order(V1, V2),], d.2[order(V1,V2),])
#[1] TRUE
Здесь используется способ data.table function foverlaps(), который также оказывается быстрым!
require(data.table) ## 1.9.4+
d[, `:=`(id1 = 1L, id2 = .I)] ## add interval columns for overlaps
setkey(d, id1, id2)
system.time(olaps <- foverlaps(d, d, type="within", which=TRUE)[xid != yid])
# 0.603 0.062 0.717
Обратите внимание, что foverlaps() не вычисляет все перестановки. Подмножество xid != yid необходимо для удаления совпадений. Подмножество может быть внутренне обработано более эффективно, реализуя аргумент ignoreSelf - аналогично IRanges::findOverlaps.
Теперь это просто вопрос выполнения подмножества с использованием полученных идентификаторов:
system.time(ans <- setDT(list(d$id[olaps$xid], d$id[olaps$yid])))
# 0.576 0.047 0.662
Итак, полностью, ~ 1,4 секунды.
Преимущество состоит в том, что вы можете сделать то же самое, даже если ваша data.table d имеет более одного столбца, на котором вы должны получить комбинации и использовать тот же объем памяти (так как мы возвращаем индексы). В этом случае вы просто выполните:
cbind(d[olaps$xid, your_cols, with=FALSE], d[olaps$yid, your_cols, with=FALSE])
Но он ограничивается заменой только combn(., 2L). Не более 2L.
Вот решение, использующее Rcpp.
library(Rcpp)
library(data.table)
cppFunction('
Rcpp::DataFrame combi2(Rcpp::CharacterVector inputVector){
int len = inputVector.size();
int retLen = len * (len-1) / 2;
Rcpp::CharacterVector outputVector1(retLen);
Rcpp::CharacterVector outputVector2(retLen);
int start = 0;
for (int i = 0; i < len; ++i){
for (int j = i+1; j < len; ++j){
outputVector1(start) = inputVector(i);
outputVector2(start) = inputVector(j);
++start;
}
}
return(Rcpp::DataFrame::create(Rcpp::Named("id") = outputVector1,
Rcpp::Named("neighbor") = outputVector2));
};
')
# Toy data
d <- data.table(id=as.character(paste0("A", 10001:15000)))
system.time({
d.2 <- d[, list(neighbor=d$id[-which(d$id==id)]), by=c("id")]
})
# 1.908 0.397 2.389
system.time({
d[, `:=`(id1 = 1L, id2 = .I)] ## add interval columns for overlaps
setkey(d, id1, id2)
olaps <- foverlaps(d, d, type="within", which=TRUE)[xid != yid]
ans <- setDT(list(d$id[olaps$xid], d$id[olaps$yid]))
})
# 0.653 0.038 0.705
system.time(ans2 <- combi2(d$id))
# 1.377 0.108 1.495
Используя функцию Rcpp, чтобы получить индексы, а затем сформировать таблицу data.table, работает лучше.
cppFunction('
Rcpp::DataFrame combi2inds(const Rcpp::CharacterVector inputVector){
const int len = inputVector.size();
const int retLen = len * (len-1) / 2;
Rcpp::IntegerVector outputVector1(retLen);
Rcpp::IntegerVector outputVector2(retLen);
int indexSkip;
for (int i = 0; i < len; ++i){
indexSkip = len * i - ((i+1) * i)/2;
for (int j = 0; j < len-1-i; ++j){
outputVector1(indexSkip+j) = i+1;
outputVector2(indexSkip+j) = i+j+1+1;
}
}
return(Rcpp::DataFrame::create(Rcpp::Named("xid") = outputVector1,
Rcpp::Named("yid") = outputVector2));
};
')
system.time({
indices <- combi2inds(d$id)
ans2 <- setDT(list(d$id[indices$xid], d$id[indices$yid]))
})
# 0.389 0.027 0.425
Сообщение с любым изменением слова Fast в заголовке является неполным без контрольных показателей. Прежде чем опубликовать какие-либо тесты, я хотел бы упомянуть, что, поскольку этот вопрос был опубликован, для R были выпущены два высоко оптимизированных пакета, arrangements и RcppAlgos (я автор) для создания комбинаций.
Чтобы дать вам представление об их скорости над combn и gRbase::combnPrim, вот базовый тест:
microbenchmark(arrangements::combinations(20, 10),
combn(20, 10),
gRbase::combnPrim(20, 10),
RcppAlgos::comboGeneral(20, 10),
unit = "relative")
Unit: relative
expr min lq mean median uq max neval
arrangements::combinations(20, 10) 1.364092 1.244705 1.198256 1.265019 1.192174 3.658389 100
combn(20, 10) 82.672684 61.589411 52.670841 59.976063 58.584740 67.596315 100
gRbase::combnPrim(20, 10) 6.650843 5.290714 5.024889 5.303483 5.514129 4.540966 100
RcppAlgos::comboGeneral(20, 10) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 100
Теперь мы сравниваем другие функции, опубликованные для конкретного случая создания комбинаций, выбираем 2 и data.table объект data.table.
Функции следующие:
funAkraf <- function(d) {
a <- comb2.int(length(d$id)) ## comb2.int from the answer given by @akraf
data.table(V1 = d$id[a[,1]], V2 = d$id[a[,2]])
}
funAnirban <- function(d) {
indices <- combi2inds(d$id)
ans2 <- setDT(list(d$id[indices$xid], d$id[indices$yid]))
ans2
}
funArrangements <- function(d) {as.data.table(arrangements::combinations(x = d$id, k = 2))}
funArun <- function(d) {
d[, ':='(id1 = 1L, id2 = .I)] ## add interval columns for overlaps
setkey(d, id1, id2)
olaps <- foverlaps(d, d, type="within", which=TRUE)[xid != yid]
ans <- setDT(list(d$id[olaps$xid], d$id[olaps$yid]))
ans
}
funGRbase <- function(d) {as.data.table(t(gRbase::combnPrim(d$id,2)))}
funOPCombn <- function(d) {as.data.table(t(combn(d$id, 2)))}
funRcppAlgos <- function(d) {as.data.table(RcppAlgos::comboGeneral(d$id, 2))}
И вот контрольные показатели на примере, приведенном OP:
d <- data.table(id=as.character(paste0("A", 10001:15000)))
microbenchmark(funAkraf(d),
funAnirban(d),
funArrangements(d),
funArun(d),
funGRbase(d),
funOPCombn(d),
funRcppAlgos(d),
times = 10, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
funAkraf(d) 2.961790 2.869365 2.612028 2.948955 2.215608 2.352351 10
funAnirban(d) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 10
funArrangements(d) 1.384152 1.427382 1.473522 1.854861 1.258471 1.233715 10
funArun(d) 2.785375 2.543434 2.353724 2.793377 1.883702 2.013235 10
funGRbase(d) 4.309175 3.909820 3.359260 3.921906 2.727707 2.465525 10
funOPCombn(d) 22.810793 21.722210 17.989826 21.492045 14.079908 12.933432 10
funRcppAlgos(d) 1.359991 1.551938 1.434623 1.727857 1.318949 1.176934 10
Мы видим, что функция, предоставляемая @AnirbanMukherjee, является самой быстрой для этой задачи, за которой следуют RcppAlgos/arrangements (очень близкие тайминги).
Все они дают одинаковый результат:
identical(funAkraf(d), funOPCombn(d))
#[1] TRUE
identical(funAkraf(d), funArrangements(d))
#[1] TRUE
identical(funRcppAlgos(d), funArrangements(d))
#[1] TRUE
identical(funRcppAlgos(d), funAnirban(d))
#[1] TRUE
identical(funRcppAlgos(d), funArun(d))
#[1] TRUE
## different order... we must sort
identical(funRcppAlgos(d), funGRbase(d))
[1] FALSE
d1 <- funGRbase(d)
d2 <- funRcppAlgos(d)
## now it the same
identical(d1[order(V1, V2),], d2[order(V1,V2),])
#[1] TRUE
Благодаря @Frank для указания того, как сравнить два data.tables не data.tables в результате создания новых data.tables а затем упорядочивая их:
fsetequal(funRcppAlgos(d), funGRbase(d))
[1] TRUE
Вот два решения base-R, если вы не хотите использовать дополнительные зависимости:
comb2.int использует rep и другие функции генерации последовательности для генерации желаемого результата.
comb2.mat создает матрицу, использует upper.tri() для получения верхнего треугольника и which(..., arr.ind = TRUE) для получения индексов столбца и строки = > всех комбинаций.
comb2.intcomb2.int <- function(n, rep = FALSE){
if(!rep){
# e.g. n=3 => (1,1), (1,2), (1,3), (2,2), (2,3), (3,3)
x <- rep(1:n,(n:1)-1)
i <- seq_along(x)+1
o <- c(0,cumsum((n-2):1))
y <- i-o[x]
}else{
# e.g. n=3 => (1,2), (1,3), (2,3)
x <- rep(1:n,n:1)
i <- seq_along(x)
o <- c(0,cumsum(n:2))
y <- i-o[x]+x-1
}
return(cbind(x,y))
}
comb2.matcomb2.mat <- function(n, rep = FALSE){
# Use which(..., arr.ind = TRUE) to get coordinates.
m <- matrix(FALSE, nrow = n, ncol = n)
idxs <- which(upper.tri(m, diag = rep), arr.ind = TRUE)
return(idxs)
}
combn(.):for(i in 2:8){
# --- comb2.int ------------------
stopifnot(comb2.int(i) == t(combn(i,2)))
# => Equal
# --- comb2.mat ------------------
m <- comb2.mat(i)
colnames(m) <- NULL # difference 1: colnames
m <- m[order(m[,1]),] # difference 2: output order
stopifnot(m == t(combn(i,2)))
# => Equal up to above differences
}
Используйте возвращаемые значения как индексы:
v <- LETTERS[1:5]
c <- comb2.int(length(v))
cbind(v[c[,1]], v[c[,2]])
#> [,1] [,2]
#> [1,] "A" "B"
#> [2,] "A" "C"
#> [3,] "A" "D"
#> [4,] "A" "E"
#> [5,] "B" "C"
#> [6,] "B" "D"
#> [7,] "B" "E"
#> [8,] "C" "D"
#> [9,] "C" "E"
#> [10,] "D" "E"
time (combn) = ~ 5x time (comb2.mat) = ~ 80x time (comb2.int):
library(microbenchmark)
n <- 800
microbenchmark({
comb2.int(n)
},{
comb2.mat(n)
},{
t(combn(n, 2))
})
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> { comb2.int(n) } 4.394051 4.731737 6.350406 5.334463 7.22677 14.68808 100
#> { comb2.mat(n) } 20.131455 22.901534 31.648521 24.411782 26.95821 297.70684 100
#> { t(combn(n, 2)) } 363.687284 374.826268 391.038755 380.012274 389.59960 532.30305 100