Versione più veloce di Combn

C’è un modo per accelerare il comando combn per ottenere tutte le combinazioni uniche di 2 elementi presi da un vettore?

Di solito questo sarebbe impostato in questo modo:

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

Tuttavia, combn è 10 volte più lento (23 secondi contro 3 secondi sul mio computer) rispetto al calcolo di tutte le combinazioni possibili usando data.table.

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

Trattando con vettori molto grandi, sto cercando un modo per risparmiare memoria calcolando solo le combinazioni univoche (come combn ), ma con la velocità di data.table (vedi il secondo snippet di codice).

Apprezzo qualsiasi aiuto.

Potresti usare combnPrim da 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 

Ecco un modo utilizzando la funzione foverlaps() , che risulta anche essere veloce!

 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 

Nota che foverlaps() non calcola tutte le permutazioni. Il sottoinsieme xid != yid è necessario per rimuovere l' auto-sovrapposizione . Il sottoinsieme potrebbe essere gestito internamente in modo più efficiente implementando ignoreSelf argomento ignoreSelf - simile a IRanges::findOverlaps .

Ora si tratta solo di eseguire un sottoinsieme usando gli ID ottenuti:

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

Così totalmente, ~ 1,4 secondi.


Il vantaggio è che puoi fare la stessa cosa anche se data.table d ha più di una colonna sulla quale devi ottenere le combinazioni e usare la stessa quantità di memoria (dato che restituiamo gli indici). In tal caso, devi solo fare:

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

Ma si limita a sostituire solo combn(., 2L) . Non più di 2L.

Ecco una soluzione che utilizza 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 

Usando la funzione Rcpp per ottenere gli indici e quindi formare data.table, funziona meglio.

 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 

Un post con qualsiasi variazione della parola Veloce nel titolo è incompleto senza benchmark. Prima di pubblicare qualsiasi benchmark, vorrei solo ricordare che da quando questa domanda è stata postata, due pacchetti altamente ottimizzati, arrangements e RcppAlgos (sono l’autore) per generare combinazioni sono stati rilasciati per R

Per darti un’idea della loro velocità rispetto a combn e gRbase::combnPrim , ecco un punto di riferimento di base:

 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 

Ora, confrontiamo le altre funzioni pubblicate per il caso specifico di combinazioni di produzione, scegliamo 2 e produciamo un object data.table .

Le funzioni sono le seguenti:

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

E qui ci sono i benchmark sull'esempio dato dall'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 

Vediamo che la funzione fornita da @AnirbanMukherjee è la più veloce per questa attività, seguita da RcppAlgos / arrangements (intervalli molto ravvicinati).

Tutti danno lo stesso risultato:

 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 

Grazie a @Frank per aver indicato come confrontare due data.tables senza passare attraverso la data.tables di creare nuovi data.tables e quindi disponendoli:

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

Ecco due soluzioni di base R se non si desidera utilizzare dipendenze aggiuntive:

  • comb2.int utilizza rep e altre funzioni generatrici di sequenze per generare l’output desiderato.

  • comb2.mat crea una matrice, usa upper.tri() per ottenere il triangolo superiore e which(..., arr.ind = TRUE) per ottenere gli indici di colonna e riga => tutte le combinazioni.

Possibilità 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)) } 

Possibilità 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) } 

Le funzioni danno lo stesso risultato di 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 } 

Ma ho altri elementi nel mio vettore rispetto agli interi sequenziali!

Utilizza i valori di ritorno come indici:

 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" 

Indice di riferimento:

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