R に Intel MKL を導入する

R

Rを高速化する方法として、Rblas.dllを入れ替えるという手法は以前からあった。
ATLASとかGotoBLASとか。
BLASの高速化 - RjpWiki

今回、Microsoft R Openで使われている Intel MKL をRblas.dllとして利用する方法が紹介されていたので ここにメモをしておく。

教えてもらったツイートはこちら。

ネタ元はこちらのTom Wenseleersさんのコメント。
MRO 3.6 coming?

さらにその元ネタはこちら。
Linking Intel's Math Kernel Library (MKL) to R on Windows - Stack Overflow


1. Intel MKL (Math Kernel Library) をインストールする。配布元はこちら。Customizable Package でOK。
https://software.intel.com/en-us/mkl/choose-download

2.MKLのインストール先の以下の中身を

C:\Program Files (x86)\IntelSWTools\compilers_and_libraries\windows\redist\intel64\mkl
C:\Program Files (x86)\IntelSWTools\compilers_and_libraries\windows\redist\intel64\compiler

Rのインストール先の以下のフォルダにコピーする。

C:\Program Files\R\R-3.x.x\bin\x64

3. 元々入っている R の Rblas.dll と Rlapack.dll を適当にリネームする。
4. C:\Program Files\R\R-3.x.x\bin\x64 内に mkl_rt.dll のコピーを二つ作り、それぞれ Rblas.dll と Rlapack.dll にリネームする。

broomパッケージを使ってみる

R

broomパッケージの上手な使い方を考えてる。

library(stargazer)
library(broom)
library(dplyr)
library(ggplot2)

res1 <- lm(yield ~ block + N + P + K, npk)
res2 <- lm(yield ~ block + N + P * K, npk)

# res <- res1
res <- res2

anova(res)
res %>% 
  aov() %>% 
  tidy()
summary(res)
tidy(res)

# 多重比較
res %>% 
  aov() %>% 
  TukeyHSD() %>% 
  tidy() %>%
  mutate(pair = paste(term, comparison, sep="_")) %>%
  ggplot(aes(colour=cut(adj.p.value, c(0, 0.01, 0.05, 1),
                        label=c("p<0.01","p<0.05","Non-Sig")))) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(pair, ymin=conf.low, ymax=conf.high), width=0.2) +
  geom_point(aes(pair, estimate)) +
  labs(colour="significance", title = "Tukey Honest Significant Differences") +
  coord_flip()

glance(res)
augment(res)

res %>% 
  tidy(conf.int = TRUE) %>% 
  filter(term != "(Intercept)") %>% 
  ggplot(aes(estimate, term, color = term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) +
  geom_vline(xintercept = 0)

stargazer(res1, res2, type="text", column.labels = c("res1","res2"),)

続々・複数のファイルを一度に読みこむ方法

r

過去2回の続きでこれが最後（多分）。
複数のファイルを一度に読みこむ方法 - bob3’s blog
続・複数のファイルを一度に読みこむ方法 - bob3’s blog

前回、data.table::fread() を foreach で並列化したパターンを追加したのだけど、data.table::fread() はそもそも並列化されているということを思い出しました。
大変な無駄をしていたわけです。

このままだと並列化によって速くなるか否かが不明なままですので、readr::read_csv() を並列化して速さを確認してみます。
結論から言うとやっぱり fread & rbindlist が最速。
でも read_csv() も並列化した分は速くなってます。

f:id:bob3:20190709231809p:plain
read_csvの並列化追加

まあ、通常は fread & rbindlist でいいと思います。

# demo data ---------------------------------------------------------------
library(tidyverse)
library(data.table)
library(doParallel)
library(microbenchmark)

DIR <- "data/"

row_n <- 10000
col_n <- 10
files <- 365
file_names <- paste0(DIR, formatC(1:files, width = 3, flag = "0"), ".csv")

for (i in file_names) {
  rnorm(col_n * row_n) %>%
    round(4) %>%
    matrix(row_n, col_n) %>%
    data.frame() %>%
    fwrite(i)
}

# read_csv + map_dfr -----------------------------------------------------------------
read_csv_and_map_dfr <- function() {
  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map_dfr(
      read_csv,
      progress = FALSE,
      .id = "file_name"
    )
}

# fread + rbindlist ---------------------------------------------------------
fread_and_rbindlist <- function() {
  file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
  tmp <- map(file_names, fread, showProgress = FALSE)
  setattr(tmp, "names", file_names)
  rbindlist(tmp, idcol = "file_name")
}

# 並列化 ---------------------------------------------------------------------
cores <- detectCores(logical = FALSE) # コア数確認
cl <- makeCluster(cores)
registerDoParallel(cl)
clusterEvalQ(cl, c(library(readr), library(foreach)))

my_foreach <- function() {
  file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
  tmp1 <- foreach(i = file_names) %dopar% {
    read_csv(i)
  }
  bind_rows(tmp1, .id = "files") %>%
    left_join(tibble(file_names = file_names, files = as.character(1:length(file_names))), by = "files") %>%
    select(-files)
}

# microbenchimark ---------------------------------------------------------

gc()
gc()

compare <- microbenchmark(
  "read_csv & map_dfr" = read_csv_and_map_dfr(),
  "fread & rbindlist" = fread_and_rbindlist(),
  "doParallel & read_csv & bind_rows" = my_foreach(),
  times = 30L
)

compare

autoplot(compare) + labs(title = "R 3.6.1")

続・複数のファイルを一度に読みこむ方法

r

前々回の続編。

前々回の記事を書いた後の反応にこんなのがありました。

なるほど、並列化という手があったか。

ということで、まだ動作が不安定な vroom を外して、

  • Tidyverse流（read_csv & map_dfr）
  • data.table流（fread & rbindlist）
  • data.table流並列化（fread & rbindlist & foreach）

の3パターンで試します。

デモデータの用意

必要なパッケージを読み込みつつ、実験用のデモデータの用意。

library(tidyverse)
library(data.table)
library(doParallel)
library(microbenchmark)

DIR <- "data/"

row_n <- 10000
col_n <- 10
files <- 365 # 前回は不安定なvroomに合わせて10ファイルだったのを増やした。
file_names <- paste0(DIR, formatC(1:files, width = 3, flag = "0"), ".csv")

for (i in file_names) {
  rnorm(col_n * row_n) %>%
    round(4) %>%
    matrix(row_n, col_n) %>%
    data.frame() %>%
    fwrite(i)
}

Tidyverse流

  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map_dfr(
      read_csv,
      progress = FALSE,
      .id = "file_name"
    ) -> dat

data.table流

  file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
  tmp <- map(file_names, fread, showProgress = FALSE)
  setattr(tmp, "names", file_names)
  dat <- rbindlist(tmp, idcol = "file_name")

data.table流並列化

cores <- detectCores(logical = FALSE) # コア数確認
cl <- makeCluster(cores)
registerDoParallel(cl)
clusterEvalQ(cl, c(library(data.table), library(foreach)))
file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
tmp <- foreach(i = file_names) %dopar% {
    fread(i, showProgress = FALSE)
  }
setattr(tmp, "names", file_names)
rbindlist(tmp, idcol = "files")

比較

ファイル数が増えると Tidyverse 流では遅く感じてしまう。
やはり、data.table流は速い。
期待した並列化だが、バラツキが大きくて一概に速いとは言えない状況。コア数多めのCPU積んでる場合はいいのかもしれない。
書きやすさのTidyverse、速さのdata.tableというところは変わらない模様。

f:id:bob3:20190709010624p:plain
複数ファイルの読み込み速度比較（並列化入り）

# demo data ---------------------------------------------------------------
library(tidyverse)
library(data.table)
library(doParallel)
library(microbenchmark)

DIR <- "data/"

row_n <- 10000
col_n <- 10
files <- 365
file_names <- paste0(DIR, formatC(1:files, width = 3, flag = "0"), ".csv")

for (i in file_names) {
  rnorm(col_n * row_n) %>%
    round(4) %>%
    matrix(row_n, col_n) %>%
    data.frame() %>%
    fwrite(i)
}


# read_csv + map_dfr -----------------------------------------------------------------
read_csv_and_map_dfr <- function() {
  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map_dfr(
      read_csv,
      progress = FALSE,
      .id = "file_name"
    )
}

# fread + rbindlist ---------------------------------------------------------
fread_and_rbindlist <- function() {
  file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
  tmp <- map(file_names, fread, showProgress = FALSE)
  setattr(tmp, "names", file_names)
  rbindlist(tmp, idcol = "file_name")
}

# 並列化 ---------------------------------------------------------------------
cores <- detectCores(logical = FALSE) # コア数確認
cl <- makeCluster(cores)
registerDoParallel(cl)
clusterEvalQ(cl, c(library(data.table), library(foreach)))

my_foreach <- function() {
  file_names <- list.files(DIR, pattern = "*.csv", full.names = TRUE)
  tmp <- foreach(i = file_names) %dopar% {
    fread(i, showProgress = FALSE)
  }
  setattr(tmp, "names", file_names)
  rbindlist(tmp, idcol = "files")
}
# stopCluster(cl)

# microbenchimark ---------------------------------------------------------

gc()
gc()

compare <- microbenchmark(
  "read_csv & map_dfr" = read_csv_and_map_dfr(),
  "fread & rbindlist" = fread_and_rbindlist(),
  "doParallel & fread & rbindlist" = my_foreach(),
  times = 10L
)

compare

autoplot(compare) + labs(title = "R 3.6.1")

base と Tidyverse と data.table と

r

「恋しさと せつなさと 心強さと」といえば『ストリートファイターII MOVIE』の挿入歌ですね。

さて、初心者にRを教えるとき、base準拠にすべきか、Tidyverse準拠にすべきか、という議論があります。
これは私も悩むところでありまして、さらに言えば data.table ってのもあります。

んで、まずはこの三つの流派（？）が存在してそれぞれこんな書き方をするんだよ、というのを示せるといいのかなーと思うわけです。

そこで試しにこんな処理を base, Tidyverse, data.tableのそれぞれで書いてみました。

  1. データセット flights の変数をdep_delay、arr_delay、carrier、air_time、distanceに絞る。
  2. いずれかの変数にNAを含む行を削除。
  3. 時速をkm/hで算出。air_timeは分単位なので60で除して時間単位に、distanceはマイル単位なので1.60934を乗じてkm単位に。
  4. carrier毎にdep_delay、arr_delay、Hourly speedの平均と件数を算出。
  5. 件数の降順に並べる。

こんな風に書いてみましたが、それぞれの書き方の違いの雰囲気ぐらいは伝わるかな？

base

library(nycflights13)
flights.df <- na.omit(as.data.frame(flights)[, c("dep_delay", "arr_delay", "carrier", "air_time", "distance")])
flights.df$HourlySpeed <- (flights.df$distance * 1.60934) / (flights.df$air_time / 60)
flights.df2 <- aggregate(
  flights.df[, c("dep_delay", "arr_delay", "HourlySpeed") ],
  list(carrier = flights.df$carrier),
  mean
)
flights.df2$n <- table(flights.df$carrier)
flights.df2[order(flights.df2$n, decreasing = TRUE), ]

Tidyverse

library(nycflights13)
library(tidyverse)
flights %>%
  select(dep_delay, arr_delay, carrier, air_time, distance) %>%
  drop_na() %>%
  mutate(HourlySpeed = (distance * 1.60934) / (air_time / 60)) %>%
  group_by(carrier) %>%
  summarise(
    mean_dep_delay = mean(dep_delay),
    mean_arr_delay = mean(arr_delay),
    mean_HourlySpeed = mean(HourlySpeed),
    n = n()
  ) %>% 
  arrange(desc(n))

data.table

library(nycflights13)
library(data.table)
flights.dt <- na.omit(as.data.table(flights)[, c("dep_delay", "arr_delay", "carrier", "air_time", "distance")])
setorder(flights.dt[
  , HourlySpeed := (distance * 1.60934) / (air_time / 60)
  ][
    j = list(
      mean_dep_delay = mean(dep_delay),
      mean_arr_delay = mean(arr_delay),
      mean_HourlySpeed = mean(HourlySpeed),
      n = .N
      ),
    by = carrier
    ],
  -n) -> flights.dt2
flights.dt2

実行速度は？

data.table圧勝！……かと思いきや、これくらいのデータ量だとTidyverseと大差ないですね。

f:id:bob3:20190707223349p:plain
実行速度比較

複数のファイルを一度に読みこむ方法

r

Rで複数のファイルを一度に読み込みたい場合がたまにある。
例えば、売り上げのデータが日付ごとに分かれたファイルに収められているような場合。
しょっちゅう発生するわけでもなく、いつも忘れてしまって調べるのに時間がかかるのでここにメモしておく。

ここではそうしたファイルが一つのフォルダに収められているものとする。

方法は大きく三つある

  • Tidyverse流（read_csv & map_dfr）
  • data.table流（fread & rbindlist）
  • vroom流（vroom）

デモデータの用意

必要なパッケージを読み込みつつ、実験用のデモデータの用意。

library(tidyverse)
library(data.table)
library(vroom)
library(microbenchmark)

DIR <- "data/"

row_n <- 10000
col_n <- 10
files <- 10
file_names <- paste0(DIR, formatC(1:files, width = 3, flag = "0"), ".csv")

for (i in file_names) {
  rnorm(col_n * row_n) %>%
    round(4) %>%
    matrix(row_n, col_n) %>%
    data.frame() %>%
    fwrite(i)
}

Tidyverse流

  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map_dfr(read_csv,
      .id = "file_name"
    ) -> dat

data.table流

  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map(fread) %>%
    rbindlist(idcol = "file_name") -> dat

vroom流

  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    vroom(id = "file_name") -> dat

比較

Tidyverseに慣れてる人なら read_csv と map_dfr を組み合わせるのが分かりやすいだろう。
data.table流はとにかく速い。爆速。
vroomはファイル数が多くなるとエラーになる。今後に期待。

f:id:bob3:20190707154441p:plain
複数ファイルの読み込み速度比較

# read_csv + map_dfr -----------------------------------------------------------------
read_csv_and_map_dfr <- function() {
  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map_dfr(read_csv,
      .id = "file_name"
    )
}

# fread + rbindlist ---------------------------------------------------------
fread_and_rbindlist <- function() {
  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    set_names() %>%
    map(fread) %>%
    rbindlist(idcol = "file_name")
}

# vroom -------------------------------------------------------------------
# ファイル数が多いとエラーになる場合がある。
my_vroom <- function() {
  DIR %>%
    list.files(pattern = "*.csv", full.names = TRUE) %>%
    vroom(id = "file_name")
}

# microbenchimark ---------------------------------------------------------

gc()
gc()

compare <- microbenchmark(
  "read_csv & map_dfr" = read_csv_and_map_dfr(),
  "fread & rbindlist" = fread_and_rbindlist(),
  "vroom" = my_vroom(),
  times = 30L
)

compare

autoplot(compare) +
  labs(title = "R 3.6.1")

以上。

R本体のアップデート時に入れるパッケージ

r

まず最初に必ず入れるパッケージ。
あとは必要になったときに入れれば大丈夫。

install.packages(
  c(
    "tidyverse",
    "data.table",
    "psych",
    "vegan",
    "plotly",
    "devtools",
    "colourpicker", # RStudio Addin
    "ggThemeAssist", # RStudio Addin
    "styler" # RStudio Addin
  ),
  Ncpus = getOption("Ncpus", 4L)
)