/ / Zwiększenie wydajności zagnieżdżonej dla pętli w R - r, wydajność, pętla for

Zwiększ wydajność Nested For Loop w R - r, wydajność, for-loop

Poniżej znajduje się podzbiór dużych ramek danych posiadających obserwacje o wielkości 158K o nazwie "sh_data".

Patient_ID Age_in_days DEMAdmNo
396076 28542 0
396076 28570 0
396076 28598 0
396076 28626 0
396076 28654 0
396076 28682 0
396076 28710 0
396076 28738 0
396076 28766 0
396076 28794 0
396076 28822 0
396076 28850 0
396076 28878 0
396076 28906 0
396076 28934 0
396076 28962 0
396076 28990 0
396076 29018 0
396076 29046 0
396076 29074 0
396076 29102 1
396076 29165 0
396076 29200 0
396076 29228 0
396076 29263 0
396076 29200 0
396076 29228 0
396076 29263 0

Próbuję obliczyć liczbę wystąpieńdla rekordu z ostatnich sześciu miesięcy, w którym 3 kolumna wynosi 1 (oznaczona jako LACE_E). Pierwszy rekord, gdzie wiek jest minimalny, wyniesie zero. I dla drugiego rekordu, jeśli różnica wieku w dniach wynosi <= 183 dni, a kolumna 3 dla pierwszego rekordu wynosi zero, to będzie to jeden i tak dalej.

Przeszukuję następujące zapytanie w R:

LACE_E <- numeric(0)

for(i in 1:length(sh_data[,1]))
{
LACE_E[i] = 0
for(j in 1:length(sh_data[,1]))
{
if(sh_data$Patient_ID[i] == sh_data$Patient_ID[j] & sh_data$Age_in_days[i] > sh_data$Age_in_days[j] & (sh_data$Age_in_days[i]- sh_data$Age_in_days[j])<= 183 & sh_data$DEMAdmNo[j] == 1)
{LACE_E[i] = LACE_E[i] + 1}
}
}

Przetwarzanie tego zapytania trwa długo. 1 godzina na przetworzenie 100 wierszy w moim systemie. Proszę pomóż!!

Odpowiedzi:

5 dla odpowiedzi № 1

Nie musisz iść Rcpp lub data.table aby uzyskać bardzo istotne usprawnienia.

Biorąc oryginalne dane i replikując je, aby uzyskać więcej użytecznych synchronizacji:

d <- read.table(head = TRUE, text =
"Patient_ID Age_in_days DEMAdmNo
396076 28542 0
396076 28570 0
396076 28598 0
396076 28626 0
396076 28654 0
396076 28682 0
396076 28710 0
396076 28738 0
396076 28766 0
396076 28794 0
396076 28822 0
396076 28850 0
396076 28878 0
396076 28906 0
396076 28934 0
396076 28962 0
396076 28990 0
396076 29018 0
396076 29046 0
396076 29074 0
396076 29102 1
396076 29165 0
396076 29200 0
396076 29228 0
396076 29263 0
396076 29200 0
396076 29228 0
396076 29263 0 ")

d <- rbind(d, d, d, d, d, d, d, d, d, d)

Twój oryginalny kod jako funkcja i przebieg czasowy:

f0 <- function(sh_data) {
LACE_E <- numeric(0)

for(i in 1:length(sh_data[,1])) {
LACE_E[i] = 0
for(j in 1:length(sh_data[,1])) {
if(sh_data$Patient_ID[i] == sh_data$Patient_ID[j] &
sh_data$Age_in_days[i] > sh_data$Age_in_days[j] &
(sh_data$Age_in_days[i]- sh_data$Age_in_days[j])<= 183 &
sh_data$DEMAdmNo[j] == 1) {
LACE_E[i] = LACE_E[i] + 1
}
}
}
}

system.time(v0 <- f0(d))
##   user  system elapsed
##  4.803   0.007   4.812

Profilowanie pokazuje około 90% czasu spędzonego na wydobywaniu kolumn $ w wewnętrznej pętli:

Rprof()
v0 <- f0(d)
Rprof(NULL)
head(summaryRprof()$by.total)
## "f0"                  4.94    100.00      0.60    12.15
## "$"                   4.24     85.83      0.72    14.57
## "$.data.frame"        3.52     71.26      0.36     7.29
## "[["                  3.16     63.97      0.46     9.31
## "[[.data.frame"       2.70     54.66      0.96    19.43
## "%in%"                0.92     18.62      0.22     4.45

Usunięcie ekstraktów z pętli znacznie poprawia wydajność:

f1 <- function(sh_data) {
LACE_E <- numeric(0)

Patient_ID <- sh_data$Patient_ID
Age_in_days <- sh_data$Age_in_days
DEMAdmNo <- sh_data$DEMAdmNo
for(i in 1:length(sh_data[,1])) {
LACE_E[i] = 0
for(j in 1:length(sh_data[,1])) {
if(Patient_ID[i] == Patient_ID[j] &
Age_in_days[i] > Age_in_days[j] &
(Age_in_days[i]- Age_in_days[j])<= 183 &
DEMAdmNo[j] == 1) {
LACE_E[i] = LACE_E[i] + 1
}
}
}
}

system.time(v1 <- f1(d))
##   user  system elapsed
##  0.163   0.000   0.164

Niemal zawsze jest to zły idead, który zaczyna od pustego wyniku i rozwija go; Wstępna alokacja wyniku jest lepszą praktyką. W tym przypadku algorytm już istnieje O(n^2) więc nie zauważasz zbyt wiele, ale to ma znaczenie, zwłaszcza po dodaniu innych ulepszeń. f2 wstępnie alokuje wynik:

f2 <- function(sh_data) {
n <- nrow(sh_data)
LACE_E <- numeric(n)

Patient_ID <- sh_data$Patient_ID
Age_in_days <- sh_data$Age_in_days
DEMAdmNo <- sh_data$DEMAdmNo
for(i in 1:n) {
LACE_E[i] = 0
for(j in 1:n) {
if(Patient_ID[i] == Patient_ID[j] &
Age_in_days[i] > Age_in_days[j] &
(Age_in_days[i]- Age_in_days[j])<= 183 &
DEMAdmNo[j] == 1) {
LACE_E[i] = LACE_E[i] + 1
}
}
}
}

system.time(v2 <- f2(d))
##   user  system elapsed
##  0.147   0.000   0.148

Korzystanie z właściwego operatora logicznego && zamiast & ulepsza to jeszcze bardziej:

f3 <- function(sh_data) {
n <- nrow(sh_data)
LACE_E <- numeric(n)

Patient_ID <- sh_data$Patient_ID
Age_in_days <- sh_data$Age_in_days
DEMAdmNo <- sh_data$DEMAdmNo
for(i in 1:n) {
LACE_E[i] = 0
for(j in 1:n) {
if(Patient_ID[i] == Patient_ID[j] &&
Age_in_days[i] > Age_in_days[j] &&
(Age_in_days[i] - Age_in_days[j]) <= 183 &&
DEMAdmNo[j] == 1) {
LACE_E[i] = LACE_E[i] + 1
}
}
}
}

system.time(v3 <- f3(d))
##   user  system elapsed
##  0.108   0.002   0.111

To są wszystkie kroki, które należy podjąć, aby przejść do Rcpp, ale nie musisz iść Rcpp aby je wziąć.

Aby uzyskać nieco większą szybkość, możesz skompilować bajt:

f3c <- compiler::cmpfun(f3)
system.time(v3 <- f3c(d))
##   user  system elapsed
## 0.036   0.000   0.036

Te obliczenia zostały wykonane w R 3.1.3. ZA microbenchmark streszczenie:

microbenchmark(f0(d), f1(d), f2(d), f3(d), f3c(d), times = 10)
## Unit: milliseconds
##   expr        min        lq       mean     median         uq        max  neval  cld
##   f0(d) 5909.39756 5924.8493 5963.63608 5947.23469 6011.94567 6048.03571    10    d
##   f1(d)  196.16466  197.3252  200.22471  197.93345  202.49236  210.22011    10   c
##   f2(d)  187.68169  190.5644  194.02454  192.47596  195.63821  204.27415    10   c
##   f3(d)  109.17816  110.6695  112.55218  111.93915  114.43341  116.92342    10  b
##  f3c(d)   37.37348   38.8757   39.34564   39.58563   40.50597   40.58568    10 a

R.version$version.string
## [1] "R version 3.1.3 Patched (2015-03-16 r68072)"

Wydanie R 3.2.0 w kwietniu ma szereg ulepszeń w tłumaczeniu i silniku kodu bajtowego, co jeszcze bardziej poprawia wydajność:

## Unit: milliseconds
##    expr        min         lq       mean     median         uq        max neval  cld
##   f0(d) 4351.33908 4429.71559 4471.32960 4479.13901 4499.39769 4601.05390    10    d
##   f1(d)  183.57765  184.68961  190.10391  187.30951  199.56235  200.57238    10   c
##   f2(d)  177.47063  181.09790  189.78291  185.58951  190.34782  233.90264    10   c
##   f3(d)  105.79767  108.02553  114.48950  110.17056  112.85710  149.42474    10  b
##  f3c(d)   14.41182   14.43227   14.70098   14.49289   14.84504   15.67709    10 a

R.version$version.string
## [1] "R Under development (unstable) (2015-03-24 r68072)"

Tak dobre praktyki programowania R i korzystanie z narzędzi do analizy wydajności może zająć dużo czasu. Jeśli chcesz dalszego ulepszenia, możesz przejść do Rcpp ale to może wystarczyć do twoich celów.


2 dla odpowiedzi nr 2

Myślę, że lepiej to osiągnąć przy użyciu Rcpp i data.table. Naprawdę nie musisz wykonywać pętli for w tym problemie.

Moje sugeruję następujące podejście?

Stwórz nowy source.cpp plik następujący (przykładowy katalog to C: Projekty)

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
List myFunction(NumericVector x,NumericVector y) {
const int n(x.size());
NumericVector res(n);
// x is age_in_days
// y is DEMAAdmNo
for (int i=1; i<n; i++)  {
res[i]=0;
for (int j=1; j<j; j++) {
if ( (x[i]>x[j]) & ((x[i]-x[j])<=183) & (y[j]==1)) {
res[i]=res[i]+1;
}
}
}
return Rcpp::List::create(_["res"] = res);
}

jeśli nie masz Rcpp pakiet zainstalowany, zrób to i wczytaj utworzony powyżej plik cpp:

Rcpp::sourceCpp("C:/Projects/source.cpp")

Następnie w głównym pliku wykonaj następujące czynności:

library(data.table) #If not installed, do install.packages("data.table")
sh_data=fread("C:/Projects/data3.csv") #Please put your correct file path here
sh_data[, LACE_E := myFunction(Age_in_days, DEMAdmNo), by=Patient_ID]

Nie mogłem zweryfikować liczb, ponieważ nie określiłeś żądanego wyniku, więc dostosuj if oświadczenie w cpp plik.

W każdym razie połączenie Rcpp i data.table zaoszczędzi ci dużo czasu. Wysoce rekomendowane.

Mam nadzieję że to pomoże.