Versión más rápida de combn

¿Hay alguna manera de acelerar el comando combn para obtener todas las combinaciones únicas de 2 elementos tomados 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")] }) 

Tratando 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 (vea el segundo fragmento de código).

Agradezco cualquier ayuda.

Puedes 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 

Esta es una forma de utilizar 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 auto superposiciones . El subconjunto podría manejarse internamente de manera más eficiente mediante la implementación del argumento ignoreSelf , similar a IRanges::findOverlaps .

Ahora solo se trata 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 cual debe obtener las combinaciones y utiliza 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 está limitado a reemplazar simplemente combn(., 2L) . No más de 2L.

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 el archivo 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 

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

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

 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 

Ahora, comparamos las otras funciones publicadas para el caso muy específico de producción de 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 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))} 

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

Vemos que la función proporcionada por @AnirbanMukherjee es la más rápida para esta tarea, seguida de RcppAlgos / arrangements (tiempos muy cercanos).

Todos dan el mismo resultado:

 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 la molestia de crear nuevos data.tables y luego organizarlos:

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

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

  • comb2.int usa rep y otras funciones de generación de secuencia para generar el resultado deseado.

  • 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){ # eg 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 <- io[x] }else{ # eg 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 <- io[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 los 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:

tiempo ( combn ) = ~ 5x tiempo ( comb2.mat ) = ~ 80x tiempo ( 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