Ta notatka będzie interesująca dla tych, którzy używają biblioteki do przetwarzania danych tabelarycznych w R — data.table, i być może będą zadowoleni z elastyczności jej zastosowania na różnych przykładach.
Zainspirowany dobrym przykładem , i mając nadzieję, że już przeczytaliście jego artykuł, proponuję głębiej zbadać kwestie optymalizacji kodu i wydajności na podstawie data.table.
Wprowadzenie: skąd pochodzi data.table?
Najlepiej rozpocząć znajomość z biblioteką nieco z daleka, a mianowicie od struktur danych, z których może być utworzony obiekt data.table (dalej, DT).
Tablica
Kod
## arrays ---------
arrmatr <- array(1:20, c(4,5))
class(arrmatr)
typeof(arrmatr)
is.array(arrmatr)
is.matrix(arrmatr)
Jedną z takich struktur jest tablica (?base::array). Podobnie jak w innych językach, tablice są tutaj wielowymiarowe. Ciekawym jest to, że na przykład tablica dwuwymiarowa zaczyna dziedziczyć właściwości od klasy macierzy (?base::matrix), a tablica jednowymiarowa, co również jest ważne, nie dziedziczy od wektora (?base::vector).
Zarazem trzeba rozumieć, że typ danych zawartych w jakimkolwiek obiekcie należy sprawdzać funkcją base::typeof, która zwraca wewnętrzny opis typu zgodnie z R Internals — ogólnym protokołem języka, związanym z pierwotnym C.
Jeszcze jedna komenda, do określenia klasy obiektu, base::class, zwraca w przypadku wektorów typ wektorowy (różni się on nazwą od wewnętrznego, ale pozwala również zrozumieć typ danych).
Lista
Z tablicy dwuwymiarowej, zwanej również macierzą, można przejść do listy (?base::list).
Kod
## lists ------------------
mylist <- as.list(arrmatr)
is.vector(mylist)
is.list(mylist)
Przy tym następuje kilka rzeczy jednocześnie:
- Zwiniecie drugiego wymiaru macierzy, czyli, uzyskujemy jednocześnie i listę, i wektor.
- Lista, w ten sposób, dziedziczy od tych klas. Należy pamiętać, że elementowi listy będzie odpowiadała jedna (skalarna) wartość z komórki macierzy-tablicy.
Dzięki temu, że lista to również wektor, można do niej stosować niektóre funkcje dla wektorów.
Ramka danych
Z listy, macierzy lub wektora można przejść do ramki danych (?base::data.frame).
Kod
## data.frames ------------
df <- as.data.frame(arrmatr)
df2 <- as.data.frame(mylist)
is.list(df)
df$V6 <- df$V1 + df$V2
Co w niej interesujące: ramka danych dziedziczy od listy! Kolumny ramki danych są komórkami listy. To będzie ważne w przyszłości, gdy będziemy korzystać z funkcji stosowanych do list.
data.table
Uzyskać DT (?data.table::data.table) można z ramki danych, listy, wektora lub macierzy. Na przykład, w ten sposób (in place).
Kod
## data.tables -----------------------
library(data.table)
data.table::setDT(df)
is.list(df)
is.data.frame(df)
is.data.table(df)
Fajne jest to, że, podobnie jak ramka danych, DT dziedziczy właściwości listy.
DT i pamięć
W przeciwieństwie do wszystkich innych obiektów w R base, DT są przekazywane przez referencję. Jeśli potrzebujesz skopiować je do nowego obszaru pamięci, potrzebna jest funkcja data.table::copy lub musisz zrobić wybór ze starego obiektu.
Kod
df2 <- df
df[V1 == 1, V2 := 999]
data.table::fsetdiff(df, df2)
df2 <- data.table::copy(df)
df[V1 == 2, V2 := 999]
data.table::fsetdiff(df, df2)
Na tym wprowadzenie się kończy. DT to kontynuacja rozwoju struktur danych w R, która zachodzi głównie poprzez rozszerzenie i przyspieszenie operacji wykonywanych na obiektach klasy data.frame. Przy tym zachowuje dziedziczenie od innych prymitywów.
Kilka przykładów użycia właściwości data.table
Jak lista…
Iterowanie po wierszach data.frame lub DT nie jest najlepszym pomysłem, ponieważ kod pętli w języku R jest znacznie wolniejszy C, a przejście w pętli przez kolumny, których zazwyczaj jest znacznie mniej, jest jak najbardziej możliwe. Przechodząc przez kolumny, pamiętamy, że każda kolumna to element listy, który zazwyczaj zawiera wektor. Operacje na wektorach są dobrze wektoryzowane w podstawowych funkcjach języka. Można również używać operatorów wyboru, charakterystycznych dla list i wektorów: `[[`, `$`.
Kod
## operations on data.tables ------------
#using list properties
df$'V1'[1]
df[['V1']]
df[[1]][1]
sapply(df, class)
sapply(df, function(x) sum(is.na(x)))
Wektoryzacja
Jeśli istnieje potrzeba przejścia przez wiersze dużego DT, najlepszym rozwiązaniem będzie napisanie funkcji z wektoryzacją. Jeśli to nie jest możliwe, należy pamiętać, że pętla wciąż DT jest szybsza od pętli w R, ponieważ wykonuje się na C.
Spróbujemy na większym przykładzie z 100K wierszy. Będziemy wydobywać pierwszą literę ze słów, które wchodzą w wektor-kolumnę w.
Zaktualizowane
Kod
library(magrittr)
library(microbenchmark)
## Większy przykład ----
rown <- 100000
dt %
.[, d := 1 + b + c + rnorm(nrow(.))]
# wektoryzacja
microbenchmark({
dt[
, first_l := unlist(strsplit(w, split = ' ', fixed = T))[1]
, by = 1:nrow(dt)
]
})
# drugi
first_l_f %
do.call(rbind, .) %>%
`[`(,1)
}
dt[, first_l := NULL]
microbenchmark({
dt[
, first_l := .(first_l_f(w))
]
})
# trzeci
first_l_f2 %
unlist %>%
matrix(nrow = 3) %>%
`[`(1,)
}
dt[, first_l := NULL]
microbenchmark({
dt[
, first_l := .(first_l_f2(w))
]
})
Pierwszy bieg z iteracją po wierszach:
Jednostka: milisekundy
wyrażenie min
{ dt[, `:=`(first_l, unlist(strsplit(w, split = " ", fixed = T))[1]), by = 1:nrow(dt)] } 439.6217
lq mean median uq max neval
451.9998 460.1593 456.2505 460.9147 621.4042 100
Drugi bieg, gdzie wektoryzacja odbywa się przez przekształcenie listy w macierz i pobieranie elementów z przekroju o indeksie 1 (to właśnie jest wektoryzacja). Poprawię się: wektoryzacja na poziomie funkcji strsplit, która potrafi przyjmować wektor jako wejście. Okazuje się, że procedura przekształcania listy w macierz jest znacznie cięższa niż sama wektoryzacja, ale w każdym razie znacznie szybsza niż niewektoryzowana wersja.
Jednostka: milisekundy
expr min lq mean median uq max neval
{ dt[, `:=`(first_l, .(first_l_f(w)))] } 93.07916 112.1381 161.9267 149.6863 185.9893 442.5199 100
Przyspieszenie według mediana w 3 razy.
Trzeci bieg, gdzie zmieniono schemat przekształcania w macierz.
Jednostka: milisekundy
expr min lq mean median uq max neval
{ dt[, `:=`(first_l, .(first_l_f2(w)))] } 32.60481 34.13679 40.4544 35.57115 42.11975 222.972 100
Przyspieszenie według mediana w 13 razy.
Z tym trzeba eksperymentować, im więcej — tym lepiej.
Jeszcze jeden przykład z wektoryzacją, gdzie także tekst, ale zbliżony do warunków rzeczywistych: różna długość słów, różna liczba słów. Należy wyciągnąć pierwsze 3 słowa. Oto jak:
Tutaj już poprzednia funkcja nie działa, ponieważ wektory mają różne długości, a my określiliśmy rozmiar macierzy. Przeróbimy to, poszukując w internecie.
Kod
# fourth
rown <- 100000
words <-
sapply(
seq_len(rown)
, function(x){
nwords <- rbinom(1, 10, 0.5)
paste(
sapply(
seq_len(nwords)
, function(x){
paste(sample(letters, rbinom(1, 10, 0.5), replace = T), collapse = '')
}
)
, collapse = ' '
)
}
)
dt <-
data.table(
w = words
, a = sample(letters, rown, replace = T)
, b = runif(rown, -3, 3)
, c = runif(rown, -3, 3)
, e = rnorm(rown)
) %>%
.[, d := 1 + b + c + rnorm(nrow(.))]
first_l_f3 <- function(sd, n)
{
l <- strsplit(sd, split = ' ', fixed = T)
maxl <- max(lengths(l))
sapply(l, "length<-", maxl) %>%
`[`(n,) %>%
as.character
}
microbenchmark({
dt[
, (paste0('w_', 1:3)) := lapply(1:3, function(x) first_l_f3(w, x))
]
})
dt[
, (paste0('w_', 1:3)) := lapply(1:3, function(x) first_l_f3(w, x))
]
Jednostka: milisekundy
expr min lq mean median
{ dt[, `:=`((paste0(«w_», 1:3)), strsplit(w, split = " ", fixed = T))] } 851.7623 916.071 1054.5 1035.199
uq max neval
1178.738 1356.816 100
Skrypt działał ze średnią prędkością 1 sekundy. Niezbyt źle.
Połączone jednym łańcuchem…
Z obiektami DT można pracować, używając chaining. Wygląda to jak przyczepienie składni nawiasów z prawej strony, w zasadzie to jest coś na kształt cukru.
Kod
# chaining
res1 <- dt[a == 'a'][sample(.N, 100)]
res2 <- dt[, .N, a][, N]
res3 <- dt[, coefficients(lm(e ~ d))[1], a][, .(letter = a, coef = V1)]
Leci w rurach…
Podobne operacje można wykonać przez piping, wygląda podobnie, ale jest funkcjonalnie bogatsze, ponieważ można używać wszelkich metod, a nie tylko DT. Wyprowadźmy współczynniki regresji logistycznej dla naszych syntetycznych danych z szeregiem filtrów na DT.
Kod
# piping
samplpe_b <- dt[a %in% head(letters), sample(b, 1)]
res4 <-
dt %>%
.[a %in% head(letters)] %>%
.[,
{
dt0 <- .SD[1:100]
quants <-
dt0[, c] %>%
quantile(seq(0.1, 1, 0.1), na.rm = T)
.(q = quants)
}
, .(cond = b > samplpe_b)
] %>%
glm(
cond ~ q -1
, family = binomial(link = "logit")
, data = .
) %>%
summary %>%
.[[12]]
Statystyka, uczenie maszynowe i inne w DT
Można używać funkcji lambda, ale czasami lepiej stworzyć je osobno, opisać cały pipeline analizy danych i naprzód — działają wewnątrz DT. Przykład wzbogacony o wszystkie wymienione wcześniej funkcje, plus kilka przydatnych rzeczy z arsenału DT (takich jak odwołanie się do samego DT wewnątrz DT za pomocą linku, dołączonych czasami nie sekwencyjnie, ale żeby było).
Kod
# function
rm(lm_preds)
lm_preds <- function(
sd, by, n
)
{
if(
n < 100 |
!by[['a']] %in% head(letters, 4)
)
{
res <-
list(
low = NA
, mean = NA
, high = NA
, coefs = NA
)
} else {
lmm <-
lm(
d ~ c + b
, data = sd
)
preds <-
stats::predict.lm(
lmm
, sd
, interval = "prediction"
)
res <-
list(
low = preds[, 2]
, mean = preds[, 1]
, high = preds[, 3]
, coefs = coefficients(lmm)
)
}
res
}
res5 <-
dt %>%
.[e < 0] %>%
.[.[, .I[b > 0]]] %>%
.[, `:=` (
low = as.numeric(lm_preds(.SD, .BY, .N)[[1]])
, mean = as.numeric(lm_preds(.SD, .BY, .N)[[2]])
, high = as.numeric(lm_preds(.SD, .BY, .N)[[3]])
, coef_c = as.numeric(lm_preds(.SD, .BY, .N)[[4]][1])
, coef_b = as.numeric(lm_preds(.SD, .BY, .N)[[4]][2])
, coef_int = as.numeric(lm_preds(.SD, .BY, .N)[[4]][3])
)
, a
] %>%
.[!is.na(mean), -'e', with = F]
# plot
plo <-
res5 %>%
ggplot +
facet_wrap(~ a) +
geom_ribbon(
aes(
x = c * coef_c + b * coef_b + coef_int
, ymin = low
, ymax = high
, fill = a
)
, size = 0.1
, alpha = 0.1
) +
geom_point(
aes(
x = c * coef_c + b * coef_b + coef_int
, y = mean
, color = a
)
, size = 1
) +
geom_point(
aes(
x = c * coef_c + b * coef_b + coef_int
, y = d
)
, size = 1
, color = 'black'
) +
theme_minimal()
print(plo)
Podsumowanie
Mam nadzieję, że udało mi się stworzyć spójną, ale oczywiście niepełną, wizję obiektu takiego jak data.table, zaczynając od jego właściwości związanych z dziedziczeniem po klasach R, a kończąc na jego własnych cechach i otoczeniu z elementów tidyverse. Mam nadzieję, że to pomoże wam lepiej poznać i zastosować tę bibliotekę do pracy i rozrywki.
Dziękuję!
Pełny kod
Kod
## load libs ----------------
library(data.table)
library(ggplot2)
library(magrittr)
library(microbenchmark)
## arrays ---------
arrmatr <- array(1:20, c(4,5))
class(arrmatr)
typeof(arrmatr)
is.array(arrmatr)
is.matrix(arrmatr)
## lists ------------------
mylist <- as.list(arrmatr)
is.vector(mylist)
is.list(mylist)
## data.frames ------------
df <- as.data.frame(arrmatr)
is.list(df)
df$V6 <- df$V1 + df$V2
## data.tables -----------------------
data.table::setDT(df)
is.list(df)
is.data.frame(df)
is.data.table(df)
df2 <- df
df[V1 == 1, V2 := 999]
data.table::fsetdiff(df, df2)
df2 <- data.table::copy(df)
df[V1 == 2, V2 := 999]
data.table::fsetdiff(df, df2)
## operations on data.tables ------------
#using list properties
df$'V1'[1]
df[['V1']]
df[[1]][1]
sapply(df, class)
sapply(df, function(x) sum(is.na(x)))
## Bigger example ----
rown <- 100000
dt <-
data.table(
w = sapply(seq_len(rown), function(x) paste(sample(letters, 3, replace = T), collapse = ' '))
, a = sample(letters, rown, replace = T)
, b = runif(rown, -3, 3)
, c = runif(rown, -3, 3)
, e = rnorm(rown)
) %>%
.[, d := 1 + b + c + rnorm(nrow(.))]
# vectorization
# zero - for loop
microbenchmark({
for(i in 1:nrow(dt))
{
dt[
i
, first_l := unlist(strsplit(w, split = ' ', fixed = T))[1]
]
}
})
# first
microbenchmark({
dt[
, first_l := unlist(strsplit(w, split = ' ', fixed = T))[1]
, by = 1:nrow(dt)
]
})
# second
first_l_f <- function(sd)
{
strsplit(sd, split = ' ', fixed = T) %>%
do.call(rbind, .) %>%
`[`(,1)
}
dt[, first_l := NULL]
microbenchmark({
dt[
, first_l := .(first_l_f(w))
]
})
# third
first_l_f2 <- function(sd)
{
strsplit(sd, split = ' ', fixed = T) %>%
unlist %>%
matrix(nrow = 3) %>%
`[`(1,)
}
dt[, first_l := NULL]
microbenchmark({
dt[
, first_l := .(first_l_f2(w))
]
})
# fourth
rown <- 100000
words <-
sapply(
seq_len(rown)
, function(x){
nwords <- rbinom(1, 10, 0.5)
paste(
sapply(
seq_len(nwords)
, function(x){
paste(sample(letters, rbinom(1, 10, 0.5), replace = T), collapse = '')
}
)
, collapse = ' '
)
}
)
dt <-
data.table(
w = words
, a = sample(letters, rown, replace = T)
, b = runif(rown, -3, 3)
, c = runif(rown, -3, 3)
, e = rnorm(rown)
) %>%
.[, d := 1 + b + c + rnorm(nrow(.))]
first_l_f3 <- function(sd, n)
{
l <- strsplit(sd, split = ' ', fixed = T)
maxl <- max(lengths(l))
sapply(l, "length<-", maxl) %>%
`[`(n,) %>%
as.character
}
microbenchmark({
dt[
, (paste0('w_', 1:3)) := lapply(1:3, function(x) first_l_f3(w, x))
]
})
dt[
, (paste0('w_', 1:3)) := lapply(1:3, function(x) first_l_f3(w, x))
]
# chaining
res1 <- dt[a == 'a'][sample(.N, 100)]
res2 <- dt[, .N, a][, N]
res3 <- dt[, coefficients(lm(e ~ d))[1], a][, .(letter = a, coef = V1)]
# piping
samplpe_b <- dt[a %in% head(letters), sample(b, 1)]
res4 <-
dt %>%
.[a %in% head(letters)] %>%
.[,
{
dt0 <- .SD[1:100]
quants <-
dt0[, c] %>%
quantile(seq(0.1, 1, 0.1), na.rm = T)
.(q = quants)
}
, .(cond = b > samplpe_b)
] %>%
glm(
cond ~ q -1
, family = binomial(link = "logit")
, data = .
) %>%
summary %>%
.[[12]]
# function
rm(lm_preds)
lm_preds <- function(
sd, by, n
)
{
if(
n < 100 |
!by[['a']] %in% head(letters, 4)
)
{
res <-
list(
low = NA
, mean = NA
, high = NA
, coefs = NA
)
} else {
lmm <-
lm(
d ~ c + b
, data = sd
)
preds <-
stats::predict.lm(
lmm
, sd
, interval = "prediction"
)
res <-
list(
low = preds[, 2]
, mean = preds[, 1]
, high = preds[, 3]
, coefs = coefficients(lmm)
)
}
res
}
res5 <-
dt %>%
.[e < 0] %>%
.[.[, .I[b > 0]]] %>%
.[, `:=` (
low = as.numeric(lm_preds(.SD, .BY, .N)[[1]])
, mean = as.numeric(lm_preds(.SD, .BY, .N)[[2]])
, high = as.numeric(lm_preds(.SD, .BY, .N)[[3]])
, coef_c = as.numeric(lm_preds(.SD, .BY, .N)[[4]][1])
, coef_b = as.numeric(lm_preds(.SD, .BY, .N)[[4]][2])
, coef_int = as.numeric(lm_preds(.SD, .BY, .N)[[4]][3])
)
, a
] %>%
.[!is.na(mean), -'e', with = F]
# plot
plo <-
res5 %>%
ggplot +
facet_wrap(~ a) +
geom_ribbon(
aes(
x = c * coef_c + b * coef_b + coef_int
, ymin = low
, ymax = high
, fill = a
)
, size = 0.1
, alpha = 0.1
) +
geom_point(
aes(
x = c * coef_c + b * coef_b + coef_int
, y = mean
, color = a
)
, size = 1
) +
geom_point(
aes(
x = c * coef_c + b * coef_b + coef_int
, y = d
)
, size = 1
, color = 'black'
) +
theme_minimal()
print(plo)
Źródło: habr.com
