Julia - 配列の実装について

  • 2つの配列を結合するときは何故 push! よりも append! が推奨されるのか?
  • sizehint! を何故使うべきなのか?

とかあたりの話。今回は一次元配列 (Vector) のみの話で、多次元配列は扱いません。

最近今更ながらデータ構造とアルゴリズムを勉強しているのですが、その手の本を読むと,

  • 配列はサイズが固定だけれども、各要素には O(1) でアクセスできる
  • リストはサイズを可変にできるけれども、各要素へのアクセスは O(n) かかる*1

と書いてあります。 一方で Julia で配列と呼ばれているものは push!pop! 使って自由に要素数変えることができるのでリストなのかと思いきや、ベンチマークを見る限り各要素へのアクセスは O(1) で行えるようなので配列のようでもあります。 配列とリストの良いとこ取りしているような性質を持っているけど、結局のところ Julia の配列って何なの?と疑問で夜も寝られない生活が続いたので、疑問を払拭すべく調べてみました。

概要

他の言語と仕組みに大差ないと思うので重要なことを手短に言うと

  • push!, pop! の計算量は O(1)
  • Julia の配列の growth factor は 2。
  • Julia の配列の capacitypop! で小さくならない。
  • 配列の大体のサイズがわかっているときは sizehint! を使うべし。
  • 頻繁に push!, pop! を使用する場合は DataStructures.jl から適切なデータ構造を使うべし。

※ 今回は Julia 1.1.0 のソースコードを元にしました。 今後のバージョンアップで変更になる可能性は十分にあります。

環境

  • ArchLinux
  • Julia 1.1.0
julia> versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-4460T CPU @ 1.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, haswell)
Environment:
  JULIA_EDITOR = nvim
  JULIA_SHELL = /bin/bash

動的配列

要素の追加

push!

Julia の場合、配列を作ると( [a, b], rand(10)とか) まずはその分のメモリが確保されデータが保存されます。

例としてここでは v = [a, b] という配列を作ったとします。 ここに push!(v, c) とすると、

  1. 元の配列の2倍の大きさのメモリ領域を確保 (_aligned_malloc, _aligned_realloc, malloc, realloc あたりを使っているっぽい)
  2. 新しく確保した領域に元のデータをコピー
  3. 3番目の要素に c を代入

という流れで要素が追加されます。

f:id:goropikarikun:20190321120035p:plain

ここでさらに push!(v, d) とすると確保した領域の4番目に空きがあるのでそこに d が代入されます。 さらに push!(v, e) とすると今度は空き領域がないので、また 1 ~ 3 の手順に従って要素が追加されます。 以降一連の操作の繰り返しで空き領域がなくなったら倍々で配列が大きくなっていきます。 Julia の場合は2倍ずつ大きくなっていきますが Java だったら1.5倍、Python だったら約1.125倍と拡張のされ方は言語によってまちまちです。この何倍されるかという量は growth factor と呼ばれます。 また動的にサイズが変わる配列は動的配列と呼ばれます。*2

今回は、配列*3全体の大きさを表す変数を capacity、要素数の変数を length と呼ぶことにします。

配列サイズが大きくなったら capacity が更新され、要素が追加されたら length が 1 増やされます。

f:id:goropikarikun:20190321120213p:plain

en.wikipedia.org

www.geeksforgeeks.org

さて、本当に capacity が倍々に増えるのか実験してみましょう。

Julia の関数に capacity を取得する関数はおそらく用意されていないので、C言語経由で読み取ってみます。

header = if ispath("/usr/include/julia/julia.h")
    "/usr/include/julia"
elseif ispath("$(dirname(Base.julia_cmd()[1]))/../include/julia/julia.h")
    "$(dirname(Base.julia_cmd()[1]))/../include/julia"
else
    println("Can't find julia.h. Please input the path to the file.")
    path = readline()
    while (!ispath(path) )
        println("Can't find julia.h. Please input the path to the file.")
        path = readline()
    end
    path
end
C_code = raw"""
#include "julia.h"

int jl_array_maxsize(jl_array_t *a) {
    return a->maxsize; // maxsize が上記の capacity に相当する
}
"""

const Clib = tempname()
open(`gcc  -I $header -fPIC -O3  -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code) 
end
arraycap(v) = ccall((:jl_array_maxsize, Clib), Int32, (Any,), v)
julia> v = Vector{Int}(undef, 0)
0-element Array{Int64,1}

julia> for i in 1:10
           push!(v, i)
           println("# of elements = $(lpad(i, 3)): cacacity = $(arraycap(v))")
       end
# of elements =   1: cacacity = 4
# of elements =   2: cacacity = 4
# of elements =   3: cacacity = 4
# of elements =   4: cacacity = 4
# of elements =   5: cacacity = 8
# of elements =   6: cacacity = 8
# of elements =   7: cacacity = 8
# of elements =   8: cacacity = 8
# of elements =   9: cacacity = 16
# of elements =  10: cacacity = 16

期待通り capacity が4, 8, 16 と倍々になることが確認できました。*4

append!

push!append! の docstring を見ると2つの配列を結合する場合は push! ではなく append! を使うことが推奨されています。

append!(v, w) としたとき裏では (w の capacity) <= (v の capacity) のとき 2 * (v の capacity) の大きさの capacity をもつ配列ができ、 (w の capacity) > (v の capacity) のときは (w の capacity) + (v の capacity) の大きさの capacity をもつ配列ができます。 w の要素を push!v に追加していった場合、配列の生成とデータのコピーが複数回実行される可能性がありますが、append! を使った場合は配列の生成とデータのコピーが1回で済ませることができます。

要素の削除

要素の削除に関しては追加と比べると単純でpop!(v) とすると単に length が 1 減ります。 length が減るだけで末端のデータを消すこともなければ新しく小さい領域を確保してくることもありません。 そのためポインタで指定すれば値を読み取ることができます。

f:id:goropikarikun:20190321170012p:plain

julia> v = rand(3)
3-element Array{Float64,1}:
 0.742018634943469 
 0.5140336610167988
 0.7877759863795775

julia> arraycap(v)
3

julia> pop!(v)
0.7877759863795775

julia> v[3] # アクセスはできない
ERROR: BoundsError: attempt to access 2-element Array{Float64,1} at index [3]
Stacktrace:
 [1] getindex(::Array{Float64,1}, ::Int64) at ./array.jl:729
 [2] top-level scope at none:0

julia> unsafe_load(pointer(v), 3) # pointer からいけば読み取れる
0.7877759863795775

julia> arraycap(v)
3 # capacity は変わらない。

先の例のように capacity = 4 に対して無駄な領域が1つ2つあったところで今どきのコンピュータに積まれているメモリ量を考えればどうってことありませんが、大きな配列から pop! で要素を大量に削除した場合無駄な領域が無視できなくなってきます。

このとき、言語によっては length が capacity に対してある一定以下になると領域を縮小する仕組みが備わっているためメモリの無駄を抑えることができます。追加するときは領域を拡大していきましたが、その反対のことをします。

一方で Julia では pop! でいくら length が小さくなっても領域を自動的に縮小する機構は備わっていません。

f:id:goropikarikun:20190321165903p:plain
※ 先に書いたとおり pop! しても要素が自体が消えることはありませんが、図ではわかりやすさのため消しています。

pop! を使ってもメモリの使用量が減らないことを確認してみます。

julia> run(`free -m`)
              total        used        free      shared  buff/cache   available
Mem:           7852        3922        1075         295        2853        3329
Swap:          4095         486        3609
Process(`free -m`, ProcessExited(0))

julia> n = 50_000_000
50000000

julia> v = Vector{Int}(undef, 0)
0-element Array{Int64,1}

julia> for i in 1:n
           push!(v, i*10)
       end

julia> run(`free -m`)
              total        used        free      shared  buff/cache   available
Mem:           7852        4289         714         288        2847        2969
Swap:          4095         486        3609
Process(`free -m`, ProcessExited(0))

julia> varinfo()
name                    size summary                        
–––––––––––––––– ––––––––––– –––––––––––––––––––––––––––––––
Base                         Module                         
C_code             159 bytes String                         
Clib                24 bytes String                         
Core                         Module                         
InteractiveUtils 163.804 KiB Module                         
Main                         Module                         
ans                286 bytes Base.Process                   
arraycap             0 bytes typeof(arraycap)               
header              69 bytes String                         
n                    8 bytes Int64                          
v                381.470 MiB 50000000-element Array{Int64,1}

julia> for i in 1:n
           pop!(v)
       end

julia> println(arraycap(v))
57944740

julia> varinfo()
name                    size summary                 
–––––––––––––––– ––––––––––– ––––––––––––––––––––––––
Base                         Module                  
C_code             159 bytes String                  
Clib                24 bytes String                  
Core                         Module                  
InteractiveUtils 163.804 KiB Module                  
Main                         Module                  
ans                  0 bytes Nothing                 
arraycap             0 bytes typeof(arraycap)        
header              69 bytes String                  
n                    8 bytes Int64                   
v                   40 bytes 0-element Array{Int64,1}

julia> run(`free -m`)
              total        used        free      shared  buff/cache   available
Mem:           7852        4310         694         287        2847        2949
Swap:          4095         486        3609
Process(`free -m`, ProcessExited(0))

varinfo で見ると v のサイズが小さくなっているように見えますが、free の結果から pop! の前後でほとんどメモリ使用量が変わっていないことがわかると思います。

メモリ上からデータが消えたわけではないので、パラメターをいじってやれば元の配列が復活します。

julia> C_code2 = raw"""
       #include "julia.h"

       void jl_set_array(jl_array_t *a, int n) {
           a->nrows = n;
           a->length = n;
           a->maxsize = n;
           a->offset = 0;

           return;
       }
       """;

julia> const Clib2 = tempname()
"/tmp/julianSkWap"

julia> open(`gcc  -I $header -fPIC -O3  -msse3 -xc -shared -o $(Clib2 * "." * Libdl.dlext) -`, "w") do f
           print(f, C_code2) 
       end

julia> setarraypar(v, n) = ccall((:jl_set_array, Clib2), Nothing, (Any,Int32,), v, n)
setarraypar (generic function with 1 method)

julia> setarraypar(v, Int32(n))

julia> v
50000000-element Array{Int64,1}:
        10
        20
        30
        40
        50
        60
        70
        80
        90
       100499999920
 499999930
 499999940
 499999950
 499999960
 499999970
 499999980
 499999990
 500000000

sizehint! で capacity を設定する

連続して push! を使っても自動的に capacity が大きくなっていきますが、最終的な配列に入る要素数が見積もれるのなら sizehint! を使って capacity を設定しましょう。 こうすることでデータのコピー回数を減らすことができます。

sizehint! を使った場合とそうでない場合を比較してみます。

julia> using BenchmarkTools

julia> @benchmark begin
           v = Vector{Int}(undef, 0)
           for i in 1:1_000_000
               push!(v, i)
           end
       end
BenchmarkTools.Trial: 
  memory estimate:  9.00 MiB
  allocs estimate:  20
  --------------
  minimum time:     7.161 ms (0.00% GC)
  median time:      11.246 ms (0.00% GC)
  mean time:        10.909 ms (4.27% GC)
  maximum time:     30.634 ms (62.88% GC)
  --------------
  samples:          440
  evals/sample:     1

julia> @benchmark begin
           v = Vector{Int}(undef, 0)
           sizehint!(v, 1_000_000)
           for i in 1:1_000_000
               push!(v, i)
           end
       end
julia> @benchmark begin
                  v = Vector{Int}(undef, 0)
                  sizehint!(v, 1_000_000)
                  for i in 1:1_000_000
                      push!(v, i)
                  end
              end
BenchmarkTools.Trial: 
  memory estimate:  7.63 MiB
  allocs estimate:  2
  --------------
  minimum time:     8.224 ms (0.00% GC)
  median time:      8.445 ms (0.00% GC)
  mean time:        8.513 ms (0.00% GC)
  maximum time:     10.265 ms (0.00% GC)
  --------------
  samples:          540
  evals/sample:     1

sizehint! を使ったほうが約20%早くなりました。 また sizehint! を使うとメモリの割り当て、GC が少なくなっていることがわかると思います。

sizehint! は capacity を大きくもできますが、小さくもできます*5 。 pop! 使って無駄な領域が大きくなった配列の capacity を小さくしてみます。

julia> n = 50_000_000
50000000

julia> v = Vector{Int}(undef, 0)
0-element Array{Int64,1}

julia> for i in 1:n
           push!(v, i*10)
       end

julia> for i in 1:n
           pop!(v)
       end

julia> run(`free -m`);
              total        used        free      shared  buff/cache   available
Mem:           7852        4217         617         311        3017        2974
Swap:          4095         484        3611

julia> sizehint!(v, 10)
0-element Array{Int64,1}

julia> GC.gc()

julia> run(`free -m`);
              total        used        free      shared  buff/cache   available
Mem:           7852        3839         972         331        3039        3331
Swap:          4095         484        3611

sizehint! をすると新たにメモリ割り当てがされるので無駄に大きくなった配列は GC の対象になり、そのうち開放されます。

DataStructures.jl の Deque

今回は詳しくは扱いませんが普通の Julia の配列で pushfirst!, popfirst! をしたときの計算量は O(n) ですが、 DataStructures.jl の Deque を使うと pushfirst!, popfirst! の計算量を O(1) に抑えることができます。 もちろん push!, pop! の計算量は O(1) のままです。

加えて Julia の配列では pop! を使っても capacity が小さくなることはありませんでしたが Deque では自動的に配列サイズも小さくなっていきます。このとき、普通の配列をsizehint! 使って小さくしたときとは違って新たなメモリの割り当ては行われません。

頻繁に push!pop! を使うことが想定される場合は普通の配列では非効率的なことがあるのでパフォーマンスを求める場合は用途にあったデータ構造を適宜選びましょう。

番外編: memory overcommit

この章は Julia とは直接の関係はありません。

sizehint! を使うと裏では realloc とか使ってメモリを確保しているはずなのですが、明らかにメモリ足りてないだろという数値を設定してもエラーが出ないときがあります。

julia> v = Vector{Int}(undef, 0)
0-element Array{Int64,1}

julia> sizehint!(v, 1024^3)
ERROR: OutOfMemoryError()
Stacktrace:
 [1] sizehint!(::Array{Int64,1}, ::Int64) at ./array.jl:1022
 [2] top-level scope at none:0

julia> sizehint!(v, 50*1024^2)
0-element Array{Int64,1}

julia> w = Vector{Int}(undef, 0)
0-element Array{Int64,1}

julia> sizehint!(w, 50*1024^2)
0-element Array{Int64,1}

本当にメモリを確保しているというのならば、今頃 800MB くらいメモリを消費していてもよさそうなものですが、システムモニターを見る限りメモリ消費量は全然増えていません。 ですが、C言語の本を読むと malloc はメモリ確保できなかったら NULL を返すと書いてあります。 もしかして、capacity 変わっていないじゃないかと疑ってみたのですが調べてみるとちゃんと変わっています。

julia> v = Vector{Int}(undef, 0)
0-element Array{Int64,1}

julia> sizehint!(v, 50*1024^2)
0-element Array{Int64,1}

julia> arraycap(v)
52428800

sizehint! 使って push! したとき GC がほぼ 0 になったことを考えるとちゃんと領域を確保しているようなのですが、システムモニターで見る限りはほぼ変わっておらず。

Julia がうまいこと処理しちゃってるのかと思い直接C言語malloc を使ってみても状況は変わらず。

どういうことかなぁと困っていたのですが、どうやら memory overcommit なるものが関係しているようです。

proger.blog10.fc2.com

上記の記事によるとまずは仮想メモリに割り当てられるとのこと。htop で見ると、なるほど確かに Julia のところが真っ赤になってる。

f:id:goropikarikun:20190322022611p:plain

malloc が失敗しなくても安心するなよという教訓を得ました。

まとめ

  • 用途に応じて最適なデータ構造を選ぶことが大事
  • pop! 使いすぎるとメモリを無駄遣いする可能性がある
  • Juliaと仲良くなるためにはC言語とも仲良くする必要がある

参考

gist216bda6b0d3fb339a238eb99179c5722

*1:より正確には片方向リストなら k 番目にアクセスする場合 O(k), 双方向リストならノード数 n のとき O( min(k, n-k) ) なわけですが、片方向リストの最悪の場合と言うことで。

*2:より正確にはサイズは動的に変わっているけれどもバックエンドに静的配列を使用しているもの。

*3:malloc で確保した領域を配列を読んでいいのかいささか疑問ですが便宜上配列と呼ぶことにします。

*4:初期の配列サイズが小さいときに限り capacity は倍々でなく 4 にセットされる。

*5:length より小さくはできません

Julia - 並列計算用のリダクション演算子を作る

ふと「そういえば、Juliaのリダクション (reduction) は分かりづらかったなぁ」と思い出したので複数の変数をリダクションする方法と自前のリダクション演算子(reduction operator, reducer)の作り方について紹介します。

環境

  • OS: ArchLinux
  • Julia 1.1.0

Julia での for 文の並列化

Julia で手っ取り早く for 文を並列化する方法は for の前に @threads または @distributed をつけることです。(望みの結果を得られるかどうかは別にして)

スレッド並列(スレッド並列では並列数を環境変数で指定する必要がある。)

julia> using Base.Threads

julia> @threads for itr in 1:n
          body
       end

プロセス並列 (Julia v0.6 以前は@parallel を使用したが、 Julia v0.7 から @distributed に変更された)

julia> using Distributed

julia> addprocs() # CPUの論理プロセッサ数分プロセスを立ち上げる

julia> @distributed for itr in 1:n
           body
       end

このとき、スレッド並列では出来ませんが、プロセス並列ではリダクションが使えます。

公式ドキュメントの例

julia> nheads = @distributed (+) for i = 1:200000000
           Int(rand(Bool))
       end
100000133

上の例でのリダクション演算子+ です。今回はこのリダクション演算子に注目していきます。

そもそもリダクションって何?

冒頭から散々リダクションという単語を使っていますが、そもそもリダクションって何でしょう? Juliaの公式ドキュメントには次のようにあります。

Many iterations run independently over several processes, and then their results are combined using some function. The combination process is called a reduction, since it is generally tensor-rank-reducing: a vector of numbers is reduced to a single number, or a matrix is reduced to a single row or column, etc.

Parallel Computing · The Julia Language

 n 人が持っていたものを合体させれば全体としては 1 個になって総数が減る(reduction)でしょ、といった感じでしょうか。

どの変数がリダクションされるの?

私が Julia で初めてリダクションを使った時に思ったことは「どの変数がリダクションされるの?」でした。

たとえば、OpenMPの場合では、「変数 sum+ でリダクションしてください」と明示的に宣言するのでわかりやすいです。

#pragma omp parallel for reduction(+:sum)
for(i=0; i<N; i++) {
    sum += a[i] * b[i];
}

OpenMP* 4.5 による新しいレベルの並列プログラミング p.35 より

一方で上記のコードを Julia で書くと以下のようになります。

sum = @distributed (+) for i in 1:N
    a[i] + b[i]
end

この例では for 文の中に1文しかないのでまだいいですが、もっと for 文の中身が複雑になったらどうすればよいのでしょうか? 例えば、並列化される for 文の中にもう一つ for 文を入れて、その途中にリダクションしたい変数が入っているケースとかはどうでしょう。

sum = @distributed (+) for i in 1:N # 並列化のために分割される for 文
    for j in 1:N # こちらの for 文は分割されない
         ︙
        foo = ... # リダクションしたい変数end
end

なかなか難しそうですが、公式ドキュメントの例の直下に書いてある様に、並列化されている for 文の end の上にリダクションしたい変数をもってくればOKです。すなわち、上記の場合だと以下のようにすれば各プロセスで計算した foo を合算できます。

sum = @distributed (+) for i in 1:N
    for j in 1:N
         ︙
        foo = ... # リダクションしたい変数end
    foo # ここにfooと書くだけ
end

使えるリダクション演算子

OpenMPで使えるリダクション演算子は大抵使えます。 https://computing.llnl.gov/tutorials/openMP/#REDUCTION

logical and と logical or に関してはそのままでは使えませんが、以下のようにすることで同じことができます。

# logical and
julia> @distributed (*)  for i in 1:10
           rand(Bool)
       end

# logical or
julia> @distributed (x,y) -> x||y  for i in 1:10
           rand(Bool)
       end

その他にも使えるものはありますが、リダクション演算子が満たすべき条件の詳細は後述の「自前のリダクション演算子を定義する」で紹介します。

複数の変数をリダクションする

2つの変数をリダクションしたい場合、単純に end の上に a, b と書くとエラーが出ます。

julia> @distributed (+)  for i in 1:10
           a = rand()
           b = randn()
           a, b
       end
ERROR: MethodError: no method matching +(::Tuple{Float64,Float64}, ::Tuple{Float64,Float64})
...

ですが、タプルでなく配列にすると通るようになります。

julia> a, b = rand(10), randn(10);

julia> @distributed (+)  for i in 1:10
           [a[i], b[i]]
       end
2-element Array{Float64,1}:
 3.634179493372805 
 0.7752150063169639

julia> sum(a)
3.634179493372805

julia> sum(b)
0.7752150063169639

この例では配列にするだけで上手く行きましたが、一般にどのような関数ならばリダクション演算子として使えるのか、いまいちピンと来ませんね。

自前のリダクション演算子を定義する

公式ドキュメントには載っていない*1ようなので試行錯誤しながらリダクション演算子が満たすべき条件を探っていきます。 今回は簡単のため + とほぼ同じ働きをするリダクション演算子を作ってみます。

docs.julialang.org

julia> using Distributed

julia> addprocs(2); # ここでは例として2並列とする。

julia> @everywhere function reducer end # メソッドを持たない関数を定義。これを最終的にリダクション演算子にする

julia> @distributed (reducer)  for i in 1:10
           i
       end
ERROR: On worker 2:
MethodError: no method matching reducer(::Int64, ::Int64)
...

他のリダクション演算子からなんとなく想像は付いていましたが、どうやらリダクション演算子は二項演算のようですね。 引数に何が入ってくるのか知りたいので @show 使って見てみましょう。

julia> @everywhere function reducer(x::Int64, y::Int64)
           @show x
           @show y
           
           return nothing
       end

julia> @distributed (reducer)  for i in 1:10
           i
       end
      From worker 2:   x = 1
      From worker 2:   y = 2
      From worker 3:   x = 6
      From worker 3:   y = 7
ERROR: On worker 2:
MethodError: no method matching reducer(::Nothing, ::Int64)
...

reducer の返り値は nothing なので、想像するに reducer(reducer(1,2), 3) という計算をしていそうです。 返り値を Int16 にしてリトライ。

julia> @everywhere function reducer(x::Int64, y::Int64)
           @show x
           @show y
                  
           return Int16(x+y)
       end
WARNING: Method definition reducer(Int64, Int64) in module Main at REPL[5]:2 overwritten at REPL[7]:2.

julia> @distributed (reducer)  for i in 1:10
           i
       end
      From worker 2:   x = 1
      From worker 2:   y = 2
      From worker 3:   x = 6
      From worker 3:   y = 7
ERROR: On worker 2:
MethodError: no method matching reducer(::Int16, ::Int64)
Closest candidates are:
  reducer(::Int64, ::Int64) at REPL[7]:2 (method too new to be called from this world context.)
...

想像があたっている気がしないでもない。 第一引数に Int16 を取る reducer を定義してみます。

julia> @everywhere function reducer(x::Int16, y::Int64)
           @show x
           @show y
                  
           return Int16(x+y)
       end

julia> @distributed (reducer)  for i in 1:10
           i
       end
      From worker 3:   x = 6
      From worker 3:   y = 7
      From worker 2:   x = 1
      From worker 2:   y = 2
      From worker 3:   x = 13
      From worker 3:   y = 8
      From worker 3:   x = 21
      From worker 3:   y = 9
      From worker 3:   x = 30
      From worker 3:   y = 10
      From worker 2:   x = 3
      From worker 2:   y = 3
      From worker 2:   x = 6
      From worker 2:   y = 4
      From worker 2:   x = 10
      From worker 2:   y = 5
ERROR: MethodError: no method matching reducer(::Int16, ::Int16)
Closest candidates are:
  reducer(::Int16, ::Int64) at REPL[9]:2
...

だいぶ進みました。 今までは ERROR: On worker 2: と出ていましたが、今回は worker と出ていないので最後の計算はマスタープロセスで実行されているようです。 reducer にメソッドを追加します。

julia> function reducer(x::Int16, y::Int16) # マスタープロセスのみで定義。
           @show x
           @show y
                  
           return Int16(x+y)
       end
reducer (generic function with 3 methods)

julia> @distributed (reducer)  for i in 1:10
           i
       end
      From worker 3:   x = 6
      From worker 3:   y = 7
      From worker 3:   x = 13
      From worker 2:   x = 1
      From worker 2:   y = 2
      From worker 2:   x = 3
      From worker 2:   y = 3
      From worker 2:   x = 6
      From worker 2:   y = 4
      From worker 2:   x = 10
      From worker 3:   y = 8
      From worker 3:   x = 21
      From worker 3:   y = 9
      From worker 3:   x = 30
      From worker 3:   y = 10
x = 15 # マスタープロセスからの出力
y = 40 # マスタープロセスからの出力
      From worker 2:   y = 5
55

結果として 55 が返ってきたので動作としてはよさげです。

今までのことを図にまとめるとこのような感じでしょうか。

f:id:goropikarikun:20190211175311p:plain

結論

以上の結果より、次のような並列計算をするとき

x = @distributed (reducer)  for i in 1:n
    ...
    f(i) # reduced value
end

f(i) の型を X, reducer(f(i), f(j)) の型を Y とすると、リダクション演算子 reducer は以下の3つのメソッドを定義すればよさそうです。

@everywhere reducer(::X, ::X)
@everywhere reducer(::Y, ::X)
reducer(::Y, ::Y)

大抵の場合においては X = Y だと思うので、reducer(::X, ::X) だけ定義すれば十分だと思いますがトリッキーなことをしたい方は上の3つを定義してください。

ちなみに範囲が 1:1 で、もはや逐次の時は reducer を使わないので f(1) が返ってきます。

最後に配列中の最大値と最小値を自作のリダクション演算子を使って求めてみます。

julia> x = rand(100);

julia> @everywhere function reducer(x::Float64, y::Float64)
           return [max(x,y), min(x,y)]
       end

julia> @everywhere function reducer(v::Vector{Float64}, y::Float64)
           return [max(v[1],y), min(v[2],y)]
       end

julia> function reducer(v1::Vector{Float64}, v2::Vector{Float64})
           return [max(v1[1], v2[1]), min(v1[2], v2[2])]
       end
reducer (generic function with 3 methods)

julia> maxval, minval = @distributed (reducer) for item in x
           item
       end
2-element Array{Float64,1}:
 0.992511684241252  
 0.01586818629402975

julia> maximum(x)
0.992511684241252

julia> minimum(x)
0.01586818629402975

出てきた数値を見る限り望み通りに動いていそうです。

終わりに

一応それっぽいリダクション演算子を作ることができましたが、Julia 本体のソースコードをちゃんと読んだわけではないので、本当に望みのリダクション演算子が作れているのか確証はありません。

make_preduce_body という関数の定義を見る限り、今回紹介した作り方で問題ないと思いますが保証はしません。

github.com

参考

*1:正確にはありますけどわかりづらい

Julia - パッケージを効率よく読み込む ~ Requires.jl のすすめ ~

Julia 1.0 が出てから半年経って Julia に慣れてきた頃なのか、私のパッケージの作り方の記事が地味に見られているようなので、パッケージを作るときに知っておいて損はない Requires.jl について紹介します。

goropikari.hatenablog.com

環境

  • Julia 1.1.0
  • Requires.jl v0.5.2

三種の神器

Julia のパッケージを作るときに役に立つパッケージはいくつかありますが、私にとっての三種の神器は Revise.jl, PkgTemplates.jl, Requires.jl です。

github.com

github.com

github.com

前半2つはパッケージ製作者目線でないと困るもので、Requires.jl はパッケージの利用者目線で使ってくれないと困るものです。

Requires.jl

Requires.jl はパッケージを効率よく読み込むためのパッケージです。

例えば次のようなパッケージを作っていたとします。

module SamplePackage

using Plots

f(x) = ... # Plots に依存しない関数
g(x) = ... # Plots に依存しない関数

function h(x)
    Plots.plot(...)
end

end # module

ここで利用者は関数 f, g はよく使うけれど、可視化のための関数 h は滅多に使わないとしましょう。

ご存知の方も多いと思いますが Plots は Julia のパッケージの中でも重量級なので読み込みだけでも5秒近くかかります。仮にこのパッケージから Plots を抜いた場合の読み込み時間が 0.1秒だったとすると、利用者からするとろくに使いもしない関数のせいで読み込み速度が遅くされていると感じるでしょう。

Requires を使うと Plots を読み込む前は h は利用できないけれど、読み込んだら h が利用できるようになる、という風にすることができます。これにより無駄なパッケージの読み込みを抑えることができます。

使い方

公式の README を見ればわかると思いますが、オプション扱いにしたいパッケージを __init__() に入れていきます。

function __init__()
    @require パッケージ名 = "パッケージのUUID" (左で指定したパッケージが読み込まれたら実行される部分)
end

パッケージの UUID は General を直接見るか、]add でパッケージを追加したあと Project.toml を覗いてください。

github.com

書き方として重量級パッケージに依存する部分と依存しない部分をファイルで分けてしまうのが楽かなと個人的には思います。

# SamplePackage.jl
module SamplePackage

using Requires
include("noplot.jl")

function __init__()
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plot.jl")
end

end # module
# noplot.jl
export hello

hello() = println("Hellow World!")
# plot.jl
using .Plots
export sinplot

function sinplot()
    println("Plot sin fn.")
    plot(sin)
end

ディレクトリ構成

SamplePackage/
├── Manifest.toml
├── Project.toml
└── src
    ├── noplot.jl
    ├── plot.jl
    └── SamplePackage.jl

実行例

julia> using SamplePackage

julia> hello()
Hellow World!

julia> sinplot()
ERROR: UndefVarError: sinplot not defined
Stacktrace:
 [1] top-level scope at none:0

julia> using Plots

julia> sinplot() # Plots を読み込むと使えるようになる。
Plot sin fn.

plot.jl 中の using .Plots. を付け忘れないように!完全に Plots をオプション扱いにしてパッケージの依存関係に入れたくない場合、. を忘れると以下のような警告が出ます。

julia> using Plots
┌ Warning: Package SamplePackage does not have Plots in its dependencies:
│ - If you have SamplePackage checked out for development and have
│   added Plots as a dependency but haven't updated your primary
│   environment's manifest file, try `Pkg.resolve()`.
│ - Otherwise you may need to report an issue with SamplePackage
└ Loading Plots into SamplePackage from project dependency, future warnings for SamplePackage are suppressed.

読み込み時間の比較

上記のパッケージと以下のパッケージの読み込み時間を比較してみます。

module SamplePackage2
export hello, sinplot

using Plots

hello() = println("Hellow World!")

function sinplot()
    println("Plot sin fn.")
    plot(sin)
end

end # module

読み込み時間

julia> @time using SamplePackage
  0.336427 seconds (455.48 k allocations: 24.660 MiB)

julia> @time using SamplePackage2
  5.265099 seconds (6.52 M allocations: 366.485 MiB, 5.65% gc time)

SamplePackage の方では Plots を読み込んでこないので読み込み時間を抑えることができました。

製作しているパッケージが計算が主でグラフに出すことはオプション機能とするならば、Requires.jl を使って Plots 依存の部分を分離すると利用者のストレスを減らすことができると思います。 他にも DifferentialEquations.jl などもヘビー級のパッケージですが、これらの重いパッケージに依存しない部分だけでも十分に利用価値があるならば、 Requires.jl を使って分離した方がユーザー思いのパッケージになると思います。

Julia - PackageCompiler.jl を使って Plots.jl を早くする

PackageCompiler.jl の README によるとついに Plots.jl がコンパイルできるようになったらしいので試してみました。

github.com

環境

各パッケージのバージョンは最後に記載の Manifest.toml, Project.toml をご覧ください。 Unix系のOSならばこれらのファイルを使って今回の私の結果を再現できると思います。

Plots.jl をコンパイルする

Windows の場合は管理者モードで Julia を起動して ImageMagick も入れる。 macOS の場合は事前に gcc を入れておく。

import Pkg
Pkg.add.(["PackageCompiler", "Plots"]) # mac の場合はこの後に build しておくと良さそうです。
using PackageCompiler
compile_package("Plots", force=false)

無事に終了すると次のように出ます。

All done
┌ Info: Not replacing system image.
└ You can start julia with `julia -J /home/goropikari/.julia/packages/PackageCompiler/oT98U/sysimg/sys.so` at a posix shell to load the compiled files.
"/home/goropikari/.julia/packages/PackageCompiler/oT98U/sysimg/sys.so"

指示に従い

julia -J /home/goropikari/.julia/packages/PackageCompiler/oT98U/sysimg/sys.so

で Julia を起動します。こうすると Plots をコンパイルした system image を読み込んできます。 引数取らずに普通に Julia を起動するとデフォルトの system image で起動します。

スピードを比較

以下のコードを3回実行し、その実行時間を計測しました。

julia> @time begin
           using Plots
           plot(rand(10))
       end

加えて、PackageCompiler を使うと Plots の読み込み等は早くなりますが、Julia 本体の起動が遅くなるので Julia の起動に掛かる時間も計測しました(Linux のみ)。

Linux

Machine Spec

julia> versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-4460T CPU @ 1.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, haswell)

デフォルト

  • コードの実行時間
31.221949 seconds
28.103826 seconds
30.471024 seconds
平均 29.932 seconds
  • Julia の起動までに掛かる時間

計測方法は以下のよう。

time julia -e 'exit()'

結果

0m0.178s
0m0.182s
0m0.166s
平均 0.175 seconds

PackageCompiler

  • コードの実行時間
0.003226 seconds
0.003102 seconds
0.003285 seconds
平均 0.003204 seconds
  • Julia の起動までに掛かる時間
0m1.418s
0m1.405s
0m1.411s
平均 1.411 seconds

Windows 10

ArchLinux とデュアルブートしている Win7 で検証したかったのですが、そもそもコンパイルに失敗したので別機で検証しました。

Machine Spec

julia> versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i3-3110M CPU @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, ivybridge)

デフォルト

  • コードの実行時間
44.759472 seconds
42.556326 seconds
37.836190 seconds
平均 41.717329 seconds

PackageCompiler

  • コードの実行時間
2.274401 seconds
2.682438 seconds
2.648822 seconds
平均 2.535220 seconds

macOS Mojave

Machine Spec

julia> versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i7-4771 CPU @ 3.50GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, haswell)

デフォルト

  • コードの実行時間
19.305938 seconds
19.382959 seconds
20.343827 seconds
平均 19.677575 seconds
  • Julia の起動までに掛かる時間
0m0.172s
0m0.185s
0m0.177s
平均 0.178s

PackageCompiler

  • コードの実行時間
1.168776 seconds
1.226039 seconds
1.186678 seconds
平均 1.193831 seconds
  • Julia の起動までに掛かる時間
0m1.098s
0m1.117s
0m1.123s
平均 1.113s

Jupyter Notebook に Kernel を追加する

コンパイル時に force=true にしていないとデフォルトで読み込む system image がコンパイルしたものにならないので、このままでは jupyter を使ったとき Plots の読み込みを早くすることができません。

そこでコンパイルした system image を読み込めるように kernel を追加します。

※注意: ここでは IJulia で Jupyter を入れたとします。 Anaconda 等でも同様にできると思います。

Linux

~/.local/share/jupyter/kernelsjulia-1.1 というフォルダがあるのでこれをコピーします。今回はフォルダ名を julia-1.1-plots としました。

cp -r ~/.local/share/jupyter/kernels/julia-1.1 ~/.local/share/jupyter/kernels/julia-1.1-plots

次に julia-1.1-plots 中の kernel.json-J/home/ユーザー名/.julia/packages/PackageCompiler/oT98U/sysimg/sys.so を以下のように追加します。 ついでに display_name も同じだと分かりづらいので適当に変えておきましょう。

{
  "display_name": "Julia 1.1.0 plot",
  "argv": [
    "/home/ユーザー名/.local/julia/julia-1.1.0/bin/julia",
    "-J",
    "/home/ユーザー名/.julia/packages/PackageCompiler/oT98U/sysimg/sys.so",
    "-i",
    "--startup-file=yes",
    "--color=yes",
    "/home/ユーザー名/.julia/packages/IJulia/4UizY/src/kernel.jl",
    "{connection_file}"
  ],
  "language": "julia",
  "env": {},
  "interrupt_mode": "signal"
}

f:id:goropikarikun:20190220181937p:plain

上手くいくとこのように kernel が追加されます。

Windows

C:\Users\ユーザー名\AppData\Roaming\jupyter\kernelsjulia-1.1 というフォルダがあるのでこれを複製します。今回はフォルダ名を julia-1.1-plots としました。

次に julia-1.1-plots 中の kernel.json-JC:\\Users\\ユーザー名\\.julia\\packages\\PackageCompiler\\oT98U\\sysimg\\sys.dll を以下のように追加します。 ついでに display_name も同じだと分かりづらいので適当に変えておきましょう。

{
  "display_name": "Julia 1.1.0 plot",
  "argv": [
    "C:\\Users\\ユーザー名\\AppData\\Local\\Julia-1.1.0\\bin\\julia.exe",
    "-J",
    "C:\\Users\\ユーザー名\\.julia\\packages\\PackageCompiler\\oT98U\\sysimg\\sys.dll",
    "-i",
    "--startup-file=yes",
    "--color=yes",
    "C:\\Users\\ユーザー名\\.julia\\packages\\IJulia\\4UizY\\src\\kernel.jl",
    "{connection_file}"
  ],
  "language": "julia",
  "env": {},
  "interrupt_mode": "message"
}

f:id:goropikarikun:20190220181937p:plain

上手くいくとこのように kernel が追加されます。

f:id:goropikarikun:20190220181720p:plain 試しにプロットしてみましたが、どうやら Jupyter 通すと実行速度が遅くなるようです。

まとめ

単純にパッケージの読み込み・グラフの描写に掛かる時間のみに着目すると PackageCompiler を使った場合、Linux ではデフォルトよりも9342倍早くなるという結果になりました。 PackageCompiler を使った場合は Julia の起動時間が遅くなりますが、それを考慮しても21倍早くなるという結果になりました。

mac, Windows の場合ではパッケージの読み込み・グラフの描写に掛かる時間のみに着目しても16倍程度しか早くならず Linux と比べるとインパクトが少ないですが、元が元なのでやって損はないです。

Julia の起動が遅いのは地味にストレスなのでプロットするか否かに応じて system image を変える方法で運用していくのが良いかと思います。毎回 Plots.jl を読み込むと言うなら

compile_package("Plots", force=true)

とするのが楽です。

参考

nextjournal.com

jupyter-client.readthedocs.io

Manifest.toml, Project.toml

PackageCompiler_Plots

Julia - IJulia のセルで Bash コマンドを使う

環境

  • Julia 1.1.0
  • OS: ArchLinux

IJulia で ; から始めると Shell mode になるので Bash*1 コマンドが使えるけれども、2行以上書くと実行できなくなる。

f:id:goropikarikun:20190201223916p:plain

f:id:goropikarikun:20190201223928p:plain

run を使えば長々書くことも出来るけれども、それは面倒だなぁというときは以下のマクロを定義する。

# https://github.com/JuliaLang/IJulia.jl/blob/90ed075e40feb97efea3e47c899dc88907963fd1/src/magics.jl#L348
macro bash_str(s) open(`bash`,"w",stdout) do io; print(io, s); end; end

定義すると以下のようにして Bash コマンドを実行できる。

bash"""
cmd1
cmd2
...
"""

f:id:goropikarikun:20190201223955p:plain

こんなマクロわざわざ覚えてられるかという場合は、notebook のセルで %%bash と入力して実行する。 そうすると上記のマクロが出てくるのでコピペして実行すれば良い。

gist3b63ae95ed6df1b16db9197ae876c9a9

*1:Shell に Bash を指定した場合