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