Wektoryzacja pętli for zawierającej instrukcję i funkcję

Question

Powtarzalny przykład kodu próbuję wektoryzacji.

W szczególności chciałbym się dowiedzieć, czy istnieje sposób na wektorowanie tej części.

nrow(iris[ which(iris$Sepal.Length > cutOffs[plotPoint] & 
        iris$Sepal.Width > cutOffs[plotPoint]), ]) 

Powiedzmy, że było użyć biblioteki plyr lub jakąś formę zastosowania, nie ma prawdopodobnie znacznie przyspieszyć, co jest naprawdę czego szukam. Zasadniczo próbuję sprawdzić, czy istnieje technika wektoryzacji, którą przeoczyłem lub której nie udało mi się przeszukać.

UPDATE:

Unit: milliseconds 
    expr   min   lq  mean  median   uq   max neval 
    op() 33663.39700 33663.39700 33663.39700 33663.39700 33663.39700 33663.39700  1 
    jr() 3976.53088 3976.53088 3976.53088 3976.53088 3976.53088 3976.53088  1 
    dd() 4253.21050 4253.21050 4253.21050 4253.21050 4253.21050 4253.21050  1 
exp() 5085.45331 5085.45331 5085.45331 5085.45331 5085.45331 5085.45331  1 
nic() 8719.82043 8719.82043 8719.82043 8719.82043 8719.82043 8719.82043  1 
    sg() 16.66177 16.66177 16.66177 16.66177 16.66177 16.66177  1 

Bardziej realistyczne zbliżanie co ja rzeczywiście robi to za

# generate data 
numObs <- 1e5 
iris <- data.frame(Sepal.Length = sample(1:numObs), Sepal.Width = sample(1:numObs)) 

cutOffs <- 1:(numObs*0.01) 

plotOutput <- matrix(nrow=length(cutOffs), ncol=2) 
colnames(plotOutput) <- c("x","y") 
plotOutput[,"y"] <- cutOffs 

a następnie w zależności od tego szczególna metoda kto woli.

Ogólnie rzecz biorąc, będzie używany na zestawy danych z 50 000 - 200 000 punktów.

Był duży skok z użyciem

sum(Sepal.Length > cutOffs[plotPoint] & Sepal.Width > cutOffs[plotPoint]) 

co jest, co mi brakuje jako bardziej optymalnego podejścia w pierwszej kolejności.

Zdecydowanie najlepszą odpowiedzią jest sg (sgibb). Kluczem jest to, że liczy się tylko najniższa z dwóch wartości w każdym rzędzie. Po dokonaniu tego mentalnego przeskoku istnieje tylko jeden wektor do rozwiązania i wektoryzacja jest dość prosta.

# cutOff should be lower than the lowest of Sepal.Length & Sepal.Width 
    m <- pmin(iris$Sepal.Length, iris$Sepal.Width) 

Answer 1

Chciałbym dodać kolejną odpowiedź:

sg <- function() { 
    # cutOff should be lower than the lowest of Sepal.Length & Sepal.Width 
    m <- pmin(iris$Sepal.Length, iris$Sepal.Width) 
    ms <- sort.int(m) 
    # use `findInterval` to find all the indices 
    # (equal to "how many numbers below") lower than the threshold 
    plotOutput[,"x"] <- length(ms)-findInterval(cutOffs, ms) 
    plotOutput 
} 

Takie podejście pozwala uniknąć for lub outer pętli i jest 4x razy szybciej niż na @ nicola podejścia:

microbenchmark(sg(), nic(), dd()) 
#Unit: microseconds 
# expr  min  lq  mean median  uq  max neval 
# sg() 88.726 104.5805 127.3172 123.2895 144.2690 232.441 100 
# nic() 474.315 526.7780 625.0021 602.3685 706.7530 997.412 100 
# dd() 669.841 736.7800 887.4873 847.7730 976.6445 2800.930 100 

identical(sg(), dd()) 
# [1] TRUE 

Answer 2

nie usunąć pętlę for, ale zakładam, że to daje pewne przyspieszenie - nie krępuj się odniesienia i daj nam znać jak to porównać na danych rzeczywistych:

for(i in seq_along(cutOffs)) { 
    x <- cutOffs[i] 
    plotOutput[i, "x"] <- with(iris, sum(Sepal.Length > x & Sepal.Width > x)) 
} 

---

Oto mały punkt odniesienia przy użyciu przykładowych danych (co jest zapewne bardzo małe, ale może dać pewne wskazówki):

library(microbenchmark) 
microbenchmark(op(), jr(), dd(), exp(), nic()) 
Unit: microseconds 
    expr  min  lq median  uq  max neval 
    op() 6745.428 7079.8185 7378.9330 9188.0175 11936.173 100 
    jr() 1335.931 1405.2030 1466.9180 1728.6595 4692.748 100 
    dd() 684.786 711.6005 758.7395 923.6670 4473.725 100 
exp() 1928.083 2066.0395 2165.6985 2392.7030 5392.475 100 
nic() 383.007 402.5495 439.3835 541.6395 851.488 100 

Funkcje stosowane w teście są zdefiniowane następująco:

op <- function(){ 
    for(plotPoint in 1:length(cutOffs)) 
    { 
    plotOutput[plotPoint, "x"] <- 
     nrow(iris[ which(iris$Sepal.Length > cutOffs[plotPoint] & 
         iris$Sepal.Width > cutOffs[plotPoint]), ]) 
    } 
    plotOutput 
} 

jr <- function() { 
    cbind(x = sapply(cutOffs, counts), y = plotOutput[,"y"]) 
} 

dd <- function() { 
    for(i in seq_along(cutOffs)) { 
    x <- cutOffs[i] 
    plotOutput[i, "x"] <- with(iris, sum(Sepal.Length > x & Sepal.Width > x)) 
    } 
    plotOutput 
} 

exp <- function() { 
    data_frame(y=cutOffs) %>% 
    rowwise() %>% 
    mutate(x = sum(iris$Sepal.Length > y & iris$Sepal.Width > y)) 
} 

nic <- function() { 
    plotOutput[,"x"]<-colSums(outer(1:nrow(iris),1:length(cutOffs),function(x,y) iris$Sepal.Length[x] > cutOffs[y] & iris$Sepal.Width[x] > cutOffs[y])) 
} 

Edytuj notatkę: wliczone podejście przez @nicola który jest obecnie najszybciej

Answer 3

Chyba coś takiego:

counts <- function(x) sum(iris$Sepal.Length > x & iris$Sepal.Width > x) 
cbind(x = sapply(cutOffs, counts), y = plotOutput[,"y"]) 

i po prostu sprawdź:

res <- cbind(x=sapply(cutOffs,counts), y=plotOutput[,"y"]) 
identical(plotOutput,res) 
[1] TRUE 

Answer 4

Możesz użyć dplyr

library(dplyr) 
data_frame(y=cutOffs) %>% 
    rowwise() %>% 
    mutate(x = sum(iris$Sepal.Length > y & iris$Sepal.Width > y)) 

Answer 5

Można użyć outer:

plotOutput[,"x"]<-colSums(outer(1:nrow(iris),1:length(cutOffs),function(x,y) iris$Sepal.Length[x] > cutOffs[y] & iris$Sepal.Width[x] > cutOffs[y])) 

Answer 6

Inną możliwością oparciu o pmin, cut i table

brk <- c(cutOffs, Inf) 
rev(cumsum(rev(table(cut(pmin(iris$Sepal.Length, iris$Sepal.Width), brk))))) 

Mniejsza przykład, który może być łatwiejsze w użyciu, jeśli chcesz pracować przez kod "od wewnątrz":

set.seed(1) 
df <- data.frame(x = sample(1:10, 6), y = sample(1:10, 6)) 
cutOffs <- seq(from = 2, to = 8, by = 2) 
brk <- c(cutOffs, Inf) 

rev(cumsum(rev(table(cut(pmin(df$x, df$y), brk))))) 
# (2,4] (4,6] (6,8] (8,Inf] 
#  4  2  1  0 

Ie, cztery wiersze z obu wartości> 2, dwa wiersze z obu wartości> 4 et.c