Hatena::ブログ(Diary)

My Life as a Mock Quant このページをアンテナに追加 RSSフィード Twitter

2014-09-23

foreachパッケージで、並列化する/しないを綺麗に書き分ける

| 08:14 | foreachパッケージで、並列化する/しないを綺麗に書き分けるを含むブックマーク

で、「これからRで並列化するならforeachパッケージ使うだろ」という話を書いたが、その続き的な?もの。

並列化する/しないを、どう美しく書き分けるにはどうしたらよいかというお話。

たとえば、自作の関数かなんかで、並列化する/しないを書き分けようとすると、引数かなんかにparallel引数を用意しておいて、以下のように書いてしまうんじゃかろうか。

parallel=TRUE
if(parallel){
  foreach(i=1:3) %dopar% {
    rnorm(i)
  }
}else{
  foreach(i=1:3) %do% {
    rnorm(i)
  }  
}

これはdo/dopar部分だけが違うだけであって、冗長なコードだ。Rだと、最近流行りのパイプ一族系ポスト

を見てもわかるように、%で挟まれた奴らは実はこれ関数なのです。なので、以下のように書くことで、上のコードは

parallel=TRUE
doop <- ifelse(parallel, `%do%`, `%dopar%`)
doop(foreach(i=1:3), {rnorm(i)})

とまで簡略化される。こっちのほうがスマートに(私には)見えるので、これでいく。

2014-09-20

foreachパッケージで並列化する時、現在の"環境"にない変数・関数は、明示的に.export引数にて指定しなければならない

| 15:39 | foreachパッケージで並列化する時、現在の"環境"にない変数・関数は、明示的に.export引数にて指定しなければならないを含むブックマーク

foreachパッケージについて

遅い遅いと巷で噂のR君も並列化すればそれなりにパフォーマンスが期待できるわけです。Rにおいてどうやって並列化をするのかというと、foreachパッケージ&関数で並列化するのが手っ取り早くて、これはざっくりでいうと

foreach(i in 1:3) %dopar%{
    #ここに処理
}

と書くだけで処理の並列化をすることが可能です。ここで%dopar%と書いている箇所を%do%に直せば並列化させないとして処理されるのも良いところだ。この裏側の仕組みとしてはforeachパッケージが各種並列化パッケージであるsnow/multicore/parallel/Rmpiへのフロントエンドとして機能していて、裏側ではそれらのパッケージに処理をブン投げる形になっているそうだ。詳しくは

を読むとよい。上述の書籍を読むに、このforeachパッケージとparallelパッケージを組み合わせた並列化がこれからのデファクトスタンダードになりそうなんで、ここではparallel+foreachパッケージでやってみる。というわけでパッケージのロードをしておこう。

library(foreach)
library(doParallel)

さらに並列化の準備として

registerDoParallel(detectCores())

も実行しておく。これはdetectCores関数で現在のマシンのコア数を取得し、それらを並列環境としてregisterDoParallel関数で登録しますよということだ。

俺がやりたかったことと、ハマったこと

適当な乱数生成をする関数を以下のように2つ用意する。f1のほうが並列化した版で、f2の方が並列化していない版に相当する。

f1 <- function()
{
  foreach(i=1:3) %dopar% {rnorm(i)}
}
f2 <- function()
{
  foreach(i=1:3) %do% {g(i)}
}

これを実行すると、正しく求めたい結果が返ってきて

> f1()
[[1]]
[1] 0.4182321

[[2]]
[1] 0.8150335 1.2575948

[[3]]
[1]  0.5326094 -1.3211135  1.7029272

> f2()
[[1]]
[1] -0.7583267

[[2]]
[1]  0.5325022 -0.4120185

[[3]]
[1]  1.58705166 -0.58542707  0.03252658

となる。一方、実際に並列化させたい処理なんかは、自作の関数を用いることも多々あるわけで、ここではその自作関数の例として関数gを定義して、さきほどと同じことをやろうとする・・・と・・・

g <- function(i){rnorm(i)}
f3 <- function()
{
  foreach(i=1:3) %dopar% {g(i)}
}

↓実行すると・・・

> f3()
Error in { : 
  task 1 failed - " 関数 "g" を見つけることができませんでした "
Called from: foreach(i = 1:3) %dopar% {
    g(i)
}

となって怒られる。これはどうも処理を並列化する際に、並列化するスレッドに対して適切な環境*1が渡っていないからのように見える。

解決策

この問題を解決するには明示的にgを並列計算する側の環境上にもexport(輸出)しますよーというのを、その名の通りの.export引数で明示すればいい。

f4 <- function()
{
  foreach(i=1:3, .export="g") %dopar% {g(i)}
}

こうすることで、ちゃんと結果を返すようになる。

> f4()
[[1]]
[1] -0.07732738

[[2]]
[1] -0.6724723  0.5245355

[[3]]
[1]  0.2095755 -0.2227492 -0.4902599

として実行することができる。もうちょっといろいろいじってみたところ、どうもforeachパッケージを用いると

  • foreachに実際に処理を移譲するコードを描いた環境の情報(ここだと関数f4内部で定義されている変数関数のみ)

しか並列計算させる環境の上には持っていってくれないということだ。同様のことはパッケージについても言えて、自作の並列化処理において、例えばXXX, YYYパッケージを用いているならば、foreach関数引数として

foreach(i=1:10, .packages=c('XXX', 'YYY'))

というようにそれらを「明示的に使いますよ」という記述が必要になる。

もうちょっと楽に

今は関数g一個を並列化させる環境上にexportさせればよかったが、もっと多数の変数関数をexportする場合に、いちいちこれを全部書きあげるってのはとても面倒だ。なので、こういう面倒なことをするのが億劫な人、特に俺は

ls(envir=parent.frame())

を.exportの引数に入れておくことにしたい。こうすると親の環境(メインの環境だと思われるところ)の情報(関数とか変数)を全部持って行ってくれるので、楽。

f4 <- function()
{
  foreach(i=1:3, .export=ls(envir=parent.frame())) %dopar% {g(i)}
}

これでもちゃんと動作することは確認できて

> f4()
[[1]]
[1] -0.04673488

[[2]]
[1] -0.3213969 -1.3645549

[[3]]
[1] -1.9164478  0.3883101  2.2534236

となる。実際にls(envir=parent.frame())がどんな出力を返しているのかは上でいう関数f4の中でprint(ls(envir=parent.frame()))というように出力させてみるといい。出力文字列の中にgという文字がみつかるはずだ。




・・・みたいな話も上述の並列化本に書いてあった・・・ちゃんと読んでから並列化すればよかったorz

*1:自分で作ったRの変数関数の集合というか、定義したもろもろのものと考えればOK

774774 2014/09/20 17:06 「上述の書籍」とは?

teramonagiteramonagi 2014/09/21 07:59 「Rによるハイパフォーマンスコンピューティング」です!

2014-09-14

Rはクラス(class)がイケてないので、毎度クロージャ(Closure, 閉包)でごまかす俺が酷い目にあった件

| 14:03 | Rはクラス(class)がイケてないので、毎度クロージャ(Closure, 閉包)でごまかす俺が酷い目にあった件を含むブックマーク

Rのクラスはあまりイケてないというか、

にも記載があるように、S3/S4/S5/S6という4つの書き方が乱立している状況で、どれを使ったらいいんじゃい状態なので、自分でコードを書かなきゃいけない時は俺はほとんど使っていない。なんで、クラスを使った際と似たような動作、クラス/メンバー変数を持たせた場合のような動作をさせるために、その代替としてクロージャーを使うわけです。例えばクラス変数Xとしてその値を3に固定して適当なシミュレーションを行いたいような状況の場合、

make_simulator <- function(x)
{
  X <- x
  function(){
    #一様乱数を一個もってきてそれをX倍して返却するシミュレーション
    X*runif(1)
  }
}

というように一旦コンストラクタ的な"関数を生成する関数"を定義して、それを用いて、その中の変数Xを固定するようにするわけです。実際には以下のように使用する。

> simulator <- make_simulator(3)
> simulator()
[1] 1.932478
> simulator()
[1] 0.1749089
> simulator()
[1] 0.9162946

自作するコードにおいては、だいたい計算条件を"あらかじめ固定"してシミュレーションさせることが多いので、これで事足りていたんだが、クロージャの中で値(a)の更新する、つまりメンバ変数的な振る舞いをさせるために以下のようなコードを書いて実行させてみると・・・

make_f <- function()
{
  #メンバー変数的にふるまって欲しいmake_fの環境に束縛されている変数a
  a <- 1
  function()
  {
    print(a)
    a <- a + 1
    print(a)
  }  
}
> f()
[1] 1
[1] 2
> f()
[1] 1
[1] 2
> f()
[1] 1
[1] 2

うへぇ、aの値が更新されているようなコードを書いたつもりが、実は更新されていないという罠。こういうときはちゃんと”一個上の環境(ここではmake_f関数)のaに代入する"という意味で、永続代入演算子<<-を使わないとだめだということです。

make_f <- function()
{
  a <- 1
  function()
  {
    print(a)
    a <<- a + 1
    print(a)
  }  
}
> f <- make_f()
> f()
[1] 1
[1] 2
> f()
[1] 2
[1] 3
> f()
[1] 3
[1] 4

ちゃんとfを実行するたびにaの値が(内部でインクリメントされて)更新されているのがわかる。

要するに、

に近しいお話。

2014-09-13

数値シミュレーションの結果をlistで受けてからのrbind_allで捌いて可視化するのが俺のデファクトスタンダード

| 08:36 | 数値シミュレーションの結果をlistで受けてからのrbind_allで捌いて可視化するのが俺のデファクトスタンダードを含むブックマーク

例えば、二次元ランダムウォークシミュレーションの結果を返す関数を作る。

randomwalk2D <- function(){cbind(x=cumsum(rnorm(100)), y=cumsum(rnorm(100)))}

これはまぁ、二次元平面上にプロットすると原点(0,0)を始点とする二次元ランダムウォークになる。

library(dplyr)
randomwalk2D() %>% plot(type="l")

f:id:teramonagi:20140913083634p:image

最近はこの手の数値シミュレーション結果を複数回走らせた結果を捌く場合、

  • 一旦リストで受けておいて、その結果をdata.frameにdplyrパッケージのrbind_all関数を使ってdata.frameに倒してから可視化なりする

ってのが、主に俺の中でのデファクトスタンダードになっている。例えば以下のように書くということだ。

library(dplyr)
res <- list()
for(i in 1:10)
{
  res[[length(res)+1]] <- data.frame(randomwalk2D(), group=LETTERS[i])
}
df <- res %>% rbind_all 

あぁ、こう書かなくても

df <- lapply(1:10, function(i)data.frame(randomwalk2D(), group=LETTERS[i])) %>% rbind_all

の方が楽でいいか。

この結果、各数値シミュレーション(ここでは二次元ランダムウォーク)の結果がアルファベットで区別されて、以下のようなdata.frameとして出力される。

> head(df)
           x         y group
1  0.3718975 0.9598459     A
2 -2.0962722 1.1514429     A
3 -2.5984220 2.0756806     A
4 -1.7751281 1.0325865     A
5 -1.4004595 0.6961468     A
6 -2.3383467 1.0103687     A

これはおもにggplot2で可視化し易くて

ggplot(df, aes(x=x, y=y)) + 
  geom_path( aes(color=group)) +
  geom_point(aes(color=group), size=3)

f:id:teramonagi:20140913083633p:image

となって、ビューチフルだなと。もっといいやり方があっても良さそうだが、今んとこの俺俺デファクトスタンダードである。

〜追記〜

「スキあらば並列計算をするおじさん」が現れてコードをforeachで並列化してくれた。感謝感謝。

2014-09-06

tupleの値をseqenceで取得したい

| 13:10 | tupleの値をseqenceで取得したいを含むブックマーク

動的に型情報などを持ってくるReflectionを使えばいいのか。

open Microsoft.FSharp.Reflection
let x = (1,2,3,4,5)
FSharpValue.GetTupleFields x 
    |> Array.iter (fun z -> printfn "%d" (unbox<int> z))

CのDLLをF#から呼びたい俺が辿った軌跡

| 11:56 | CのDLLをF#から呼びたい俺が辿った軌跡を含むブックマーク

はじめに

Cで作ったDLLをF#から呼び出す方法を探していたら、どうもplatform invokeだのmarshallingだのを気にしないとだめらしい。C#を使わない俺はその辺よくわかってないんだが、F#でやっている例として

がよくまとめっていて良さそうだなと思ってマネしようと思った。んで、そこで使っているF# PowerPackが、

を見ると、F# PowerPack自体はより細かいパッケージへと分裂していったようだ。あ、あれ?どう真似したらいいの、ぼく?

そして、俺の旅がはじまった。

注意

C++だといろいろめんどいので、全部CのコードとしてDLLからexportすることにした。

試行錯誤をまとめたコードとその結果

sandbox/CallCppFromFs at master ? teramonagi/sandbox ? GitHub

にある。俺用。

はじめは

のあたりを見ていたのだが、カーベジコレクションされないようにメモリ確保するのどうやんのかよくわかんなくて、適当にコード書いたら動いたのがあったのでそれでいきたい。Cでやむなく書く部分ってのは、大体落っこちてる巣値計算ライブラリを叩くためだったりするわけなので、

  • 整/小数
  • それらの配列
  • 文字

というネイティブな型が受け渡しできればよいだろう。その辺をまとめておくと以下の表。

C++F#
intint
double*float[]/nativeint
charchar
char*string(入力のみ?)
struct[<Struct; StructLayout(LayoutKind.Sequential)>]属性を付けてtypeする
関数ポインタdelegate(stdcallのみっぽい)

C++の中で確保したメモリは

  • nativeint型(System.IntPtrの型略称)

になっていて、こいつをちゃんとその他の型のポインタにキャストしたい場合は

  • NativePtr.ofNativeInt

を用いる。そして、その要素を取得したい場合は

  • NativePtr.get

を用いる。

また

でnativeptr<T'>にキャストできる模様。

関数ポインタはC側で

  • typedef double(__stdcall *Function)(double, double);

のようにtypedefしたものをF#側で

  • type Func = delegate of float * float -> float

のようにdelegateを使って宣言する。括弧つき

  • type Func = delegate of (float * float) -> float

にするとダメだった、意味が変わってくるのかな?

この時の呼び出し規約はどうもstdcallのみっぽい。

これは、F#3.1の言語仕様のドラフトっぽいもの

をみると、呼び出し規約をいじれる属性であるUnmanagedFunctionPointerAttributeを使うなと書いてある。なんで、関数ポインタとしてはstdcallしか使えないっぽい。

参考

とりあえずこいつの18章に比較的詳しく記載があるのは見た。構造体やクラスの受け渡しもちゃんと載ってる。

monoのページだけど、マーシャリング周りの理解に良い。

ちゃんとやるならこれなのか?サンプルコードがググってもかからなくて、使い方がわからん。。。

データ「C++⇔F#」間相互変換周りはこの辺を読む必要がある。