tablas - R: forma más rápida de verificar la presencia de cada elemento de un vector en cada una de las columnas de una matriz

tablas en r

Question

Tengo un vector entero a

a=function(l) as.integer(runif(l,1,600))
a(100)
  [1] 414 476   6  58  74  76  45 359 482 340 103 575 494 323  74 347 157 503 385 518 547 192 149 222 152  67 497 588 388 140 457 429 353
 [34] 484  91 310 394 122 302 158 405  43 300 439 173 375 218 357  98 196 260 588 499 230  22 369  36 291 221 358 296 206  96 439 423 281
 [67] 581 127 178 330 403  91 297 341 280 164 442 114 234  36 257 307 320 307 222  53 327 394 467 480 323  97 109 564 258   2 355 253 596
[100] 215

y una matriz B entera

B=function(c) matrix(as.integer(runif(5*c,1,600)),nrow=5)
B(10)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]  250  411  181  345    4  519  167  395  130   388
[2,]  383  377  555  304  119  317  586  351  136   528
[3,]  238  262  513  476  579  145  461  191  262   302
[4,]  428  467  217  590   50  171  450  189  140   158
[5,]  178   14   31  148  285  365  515   64  166   584

y me gustaría hacer una nueva matriz booleana lxc que muestre si cada elemento vectorial en a está presente en cada columna específica de la matriz B .

Intenté esto con

ispresent1 = function (a,B) { 
out = outer(a, B, FUN = "==" )
apply(out,c(1,3),FUN="any") }

o con

ispresent2 = function (a,B) t(sapply(1:length(a), function(i) apply(B,2,function(x) a[[i]] %in% x)))

pero ninguna de estas formas de hacerlo es muy rápida:

a1=a(1000)
B1=B(20000)
system.time(ispresent1(a1,B1))
   user  system elapsed 
  76.63    1.08   77.84 

system.time(ispresent2(a1,B1))
   user  system elapsed 
  218.10    0.00   230.00 

(en mi aplicación, la matriz B tendría alrededor de 500 000 - 2 millones de columnas)

Probablemente esto sea algo trivial, pero ¿cuál es la forma correcta de hacerlo?

EDITAR: la sintaxis adecuada, como se menciona a continuación, es ispresent = function (a,B) apply(B,2,function(x) { a %in% x } ) , ¡pero la solución Rcpp continuación sigue siendo casi 2 veces más rápida! ¡Gracias por esto!

Answer 1

Después de excavar un poco, y por curiosidad sobre la respuesta Rcpp de @Backlin, escribí un punto de referencia de la solución orignal y nuestras dos soluciones:

Tuve que cambiar un poco la función de Backlin ya que no funcionó en línea en mi ventana de Windows (lo siento si me perdí algo con eso, avíseme si hay algo para adaptar)

Código utilizado:

set.seed(123) # Fix the generator
a=function(l) as.integer(runif(l,1,600))
B=function(c) matrix(as.integer(runif(5*c,1,600)),nrow=5)

ispresent1 = function (a,B) { 
out = outer(a, B, FUN = "==" )
apply(out,c(1,3),FUN="any") }

a1=a(1000)
B1=B(20000)

tensibai <- function(v,m) {
  apply(m,2,function(x) { v %in% x })
}

library(Rcpp)
cppFunction("LogicalMatrix backlin(IntegerVector a,IntegerMatrix B) {
    IntegerVector av(a);
    IntegerMatrix Bm(B);
    int i,j,k;
    LogicalMatrix out(av.size(), Bm.ncol());
    for(i = 0; i < av.size(); i++){
        for(j = 0; j < Bm.ncol(); j++){
            for(k = 0; k < Bm.nrow() && av[i] != Bm(k, j); k++);
            if(k < Bm.nrow()) out(i, j) = true;
        }
    }
    return(out);
}")

Validación:

> identical(ispresent1(a1,B1),tensibai(a1,B1))
[1] TRUE
> identical(ispresent1(a1,B1),backlin(a1,B1))
[1] TRUE

Punto de referencia:

> library(microbenchmark)
> microbenchmark(ispresent1(a1,B1),tensibai(a1,B1),backlin(a1,B1),times=3)

Unit: milliseconds
               expr        min         lq       mean     median         uq        max neval
 ispresent1(a1, B1) 36358.4633 36683.0932 37312.0568 37007.7231 37788.8536 38569.9840     3
   tensibai(a1, B1)   603.6323   645.7884   802.0970   687.9445   901.3294  1114.7144     3
    backlin(a1, B1)   471.5052   506.2873   528.3476   541.0694   556.7689   572.4684     3

La solución de Backlin es ligeramente más rápida, demostrando nuevamente que Rcpp es una buena opción si conoces cpp al principio :)

Answer 2

Rcpp es increíble para problemas como este. Es muy posible que haya alguna manera de hacerlo con data.table o con una función existente, pero con el paquete en inline , toma menos tiempo escribirlo usted mismo que averiguarlo.

require(inline)

ispresent.cpp <- cxxfunction(signature(a="integer", B="integer"),
                             plugin="Rcpp", body=''
    IntegerVector av(a);
    IntegerMatrix Bm(B);
    int i,j,k;
    LogicalMatrix out(av.size(), Bm.ncol());
    for(i = 0; i < av.size(); i++){
        for(j = 0; j < Bm.ncol(); j++){
            for(k = 0; k < Bm.nrow() && av[i] != Bm(k, j); k++);
            if(k < Bm.nrow()) out(i, j) = true;
        }
    }
    return(out);
'')

set.seed(123)
a1 <- a(1000)
B1 <- B(20000)
system.time(res.cpp <- ispresent.cpp(a1, B1))

   user  system elapsed 
  0.442   0.005   0.446

res1 <- ispresent1(a1,B1)
identical(res1, res.cpp)

[1] TRUE

Answer 3

a=function(l) as.integer(runif(l,1,600))
B=function(c) matrix(as.integer(runif(5*c,1,600)),nrow=5)

ispresent1 = function (a,B) { 
  out = outer(a, B, FUN = "==" )
  apply(out,c(1,3),FUN="any") }

ispresent2 = function (a,B) t(sapply(1:length(a), function(i) apply(B,2,function(x) a[[i]] %in% x)))

ispresent3<-function(a,B){
  tf<-matrix((B %in% a),nrow=5)
  sapply(1:ncol(tf),function(x) a %in% B[,x][tf[,x]])
}

a1=a(1000)
B1=B(20000)

> system.time(ispresent1(a1,B1))
   user  system elapsed 
  29.91    0.48   30.44 

> system.time(ispresent2(a1,B1))
   user  system elapsed 
  89.65    0.15   89.83 

> system.time(ispresent3(a1,B1))
   user  system elapsed 
   0.83    0.00    0.86 

res1<-ispresent1(a1,B1)
res3<-ispresent3(a1,B1)

> identical(res1,res3)
[1] TRUE