Versión más rápida de combn

data.table combinations

Question

¿Hay alguna manera de acelerar el comando combn para obtener todas las combinaciones únicas de 2 elementos tomadas de un vector?

Por lo general, esto se configuraría así:

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

Sin embargo, combn es 10 veces más lento (23 segundos frente a 3 segundos en mi computadora) que calcular todas las combinaciones posibles usando data.table.

system.time({
d.2 <- d[, list(neighbor=d$id[-which(d$id==id)]), by=c("id")]
})

Al tratar con vectores muy grandes, estoy buscando una manera de ahorrar memoria calculando solo las combinaciones únicas (como combn ), pero con la velocidad de data.table (consulte el segundo fragmento de código).

Agradezco cualquier ayuda.

Answer 1

Aquí hay una manera de usar la función foverlaps() , ¡que también resulta ser rápida!

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

Tenga en cuenta que foverlaps() no calcula todas las permutaciones. El subconjunto xid != yid es necesario para eliminar las superposiciones . El subconjunto podría manejarse internamente de manera más eficiente implementando el argumento ignoreSelf , similar a IRanges::findOverlaps .

Ahora es solo cuestión de realizar un subconjunto utilizando los identificadores obtenidos:

system.time(ans <- setDT(list(d$id[olaps$xid], d$id[olaps$yid])))
# 0.576 0.047 0.662 

Entonces, totalmente, ~ 1.4 segundos.

La ventaja es que puede hacer lo mismo incluso si su data.table d tiene más de 1 columna para la que debe obtener las combinaciones y usar la misma cantidad de memoria (ya que devolvemos los índices). En ese caso, simplemente harías:

cbind(d[olaps$xid, your_cols, with=FALSE], d[olaps$yid, your_cols, with=FALSE])

Pero se limita a reemplazar solo combn(., 2L) . No más de 2L.

Answer 2

Aquí hay una solución usando 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 

Usar la función Rcpp para obtener los índices y luego formar data.table, funciona mejor.

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 

Answer 3

Podrías usar combnPrim de 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

Answer 4

Una publicación con cualquier variación de la palabra Rápido en el título está incompleta sin puntos de referencia. Antes de publicar cualquier punto de referencia, me gustaría mencionar que desde que se publicó esta pregunta, se han lanzado dos paquetes altamente optimizados, arrangements y RcppAlgos (soy el autor) para generar combinaciones para R Tenga en cuenta que desde la versión 2.3.0 para RcppAlgos podemos aprovechar múltiples hilos para una eficiencia aún mayor.

Para darle una idea de su velocidad sobre combn y gRbase::combnPrim , aquí hay un punto de referencia básico:

## We test generating just over 3 million combinations
choose(25, 10)
[1] 3268760
microbenchmark(arrngmnt = arrangements::combinations(25, 10),
combn = combn(25, 10),
gRBase = gRbase::combnPrim(25, 10),
serAlgos = RcppAlgos::comboGeneral(25, 10),
parAlgos = RcppAlgos::comboGeneral(25, 10, nThreads = 4),
unit = "relative", times = 20)
Unit: relative
expr min lq mean median uq max neval
arrngmnt 2.979378 3.072319 1.898390 3.756307 2.139258 0.4842967 20
combn 226.470755 230.410716 118.157110 232.905393 125.718512 17.7778585 20
gRBase 34.219914 34.209820 18.789954 34.218320 19.934485 3.6455493 20
serAlgos 2.836651 3.078791 2.458645 3.703929 2.231475 1.1652445 20
parAlgos 1.000000 1.000000 1.000000 1.000000 1.000000 1.0000000 20

Ahora, comparamos las otras funciones publicadas para el caso muy específico de producir combinaciones, elija 2 y produzca un objeto data.table .

Las funciones son las siguientes:

funAkraf <- function(d) {
a <- comb2.int(length(d$id)) ## comb2.int from the answer given by @akraf
setDT(list(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
}
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
}
funArrangements <- function(d) {
a <- arrangements::combinations(x = d$id, k = 2)
setDT(list(a[, 1], a[, 2]))
}
funGRbase <- function(d) {
a <- gRbase::combnPrim(d$id,2)
setDT(list(a[1, ], a[2, ]))
}
funOPCombn <- function(d) {
a <- combn(d$id, 2)
setDT(list(a[1, ], a[2, ]))
}
funRcppAlgos <- function(d) {
a <- RcppAlgos::comboGeneral(d$id, 2, nThreads = 4)
setDT(list(a[, 1], a[, 2]))
}

Benchmark con datos OP

Y aquí están los puntos de referencia sobre el ejemplo dado por el 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) 3.220550 2.971264 2.815023 2.665616 2.344018 3.383673 10
funAnirban(d) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 10
funArrangements(d) 1.464730 1.689231 1.834650 1.960233 1.932361 1.693305 10
funArun(d) 3.256889 2.908075 2.634831 2.729180 2.432277 2.193849 10
funGRbase(d) 3.513847 3.340637 3.327845 3.196399 3.291480 3.129362 10
funOPCombn(d) 30.310469 26.255374 21.656376 22.386270 18.527904 15.626261 10
funRcppAlgos(d) 1.676808 1.956696 1.943773 2.085968 1.949133 1.804180 10

Vemos que la función provista por @AnirbanMukherjee es la más rápida para esta tarea, seguida de RcppAlgos / arrangements . Para esta tarea, nThreads no tiene ningún efecto ya que el vector pasado es un character , que no es seguro para subprocesos. ¿Qué pasa si en su lugar convertimos id a un factor?

Puntos de referencia con factores (es decir, variables categóricas)

dFac <- d
dFac$id <- as.factor(dFac$id)
library(microbenchmark)
microbenchmark(funAkraf(dFac),
funAnirban(dFac),
funArrangements(dFac),
funArun(dFac),
funGRbase(dFac),
funOPCombn(dFac),
funRcppAlgos(dFac),
times = 10, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
funAkraf(dFac) 10.898202 10.949896 7.589814 10.01369 8.050005 5.557014 10
funAnirban(dFac) 3.104212 3.337344 2.317024 3.00254 2.471887 1.530978 10
funArrangements(dFac) 2.054116 2.058768 1.858268 1.94507 2.797956 1.691875 10
funArun(dFac) 10.646680 12.905119 7.703085 11.50311 8.410893 3.802155 10
funGRbase(dFac) 16.523356 21.609917 12.991400 19.73776 13.599870 6.498135 10
funOPCombn(dFac) 108.301876 108.753085 64.338478 95.56197 65.494335 28.183104 10
funRcppAlgos(dFac) 1.000000 1.000000 1.000000 1.00000 1.000000 1.000000 10 

Ahora, vemos que RcppAlgos es aproximadamente 2x RcppAlgos más rápido que cualquier otra solución. En particular, la solución RcppAlgos es aproximadamente 3x que la solución más rápida que ofrecía Anirban. Cabe señalar que este aumento en la eficiencia fue posible porque las variables de factor son realmente integers debajo del capó con algunos attributes adicionales.

Confirmar la igualdad

Todos dan el mismo resultado también. La única advertencia es que la solución gRbase no admite factores. Es decir, si pasa un factor , se convertirá en character . Por lo tanto, todas las soluciones darán el mismo resultado si pasara dFac excepto la solución gRbase :

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''s the same
identical(d1[order(V1, V2),], d2[order(V1,V2),])
#[1] TRUE

Gracias a @Frank por señalar cómo comparar dos data.tables sin pasar por los dolores de crear nuevos data.tables y luego organizarlos:

fsetequal(funRcppAlgos(d), funGRbase(d))
[1] TRUE

Answer 5

Aquí hay dos soluciones de base R si no desea utilizar dependencias adicionales:

• comb2.int usa rep y otras funciones generadoras de secuencia para generar la salida deseada.

• comb2.mat crea una matriz, usa upper.tri() para obtener el triángulo superior y which(..., arr.ind = TRUE) para obtener los índices de columna y fila => todas las combinaciones.

Posibilidad 1: comb2.int

comb2.int <- function(n, rep = FALSE){
if(!rep){
# e.g. n=3 => (1,2), (1,3), (2,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,1), (1,2), (1,3), (2,2), (2,3), (3,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))
}

Posibilidad 2: comb2.mat

comb2.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)
}

Las funciones dan el mismo resultado que 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
}

¡Pero tengo otros elementos en mi vector que no son enteros secuenciales!

Use los valores de retorno como índices:

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"

Punto de referencia:

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