r data.table combinations combn

Versión más rápida de combn



data.table combinations (5)

¿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.


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.


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


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


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



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