JuliaでPyPlot.jl (matplotlib) を使ったヒストグラムのアニメーションを効率的に作成する

Juliaでも`PyPlot.jl`によってmatplotlibが使えるそうです。この記事はJuliaでアニメーションを作る話がきっかけになっていますが、実質は「matplotlibの`plt.hist`あるいは`ax.hist`は高コストなメソッドなのでヒストグラムのアニメーションを作るなら直接長方形を更新していくべき」という話です。

# きっかけ
最近の日本語インターネットでのJuliaの盛り上がりに一役買っていると思われる黒木さんの[このJupyter nobebook](http://nbviewer.jupyter.org/gist/genkuroki/3d6dbf52a3e52eb7c664bc88632c81d3#GIF%E5%8B%95%E7%94%BB%E4%BD%9C%E6%88%90)ではJuliaで様々なグラフを描く方法が紹介されていますが、PyPlot (matplotlib) でヒストグラムのアニメーションを描くとむちゃくちゃ遅いという話になっています。

コードを見ると、`FuncAnimation`の更新用関数で`Figure`や`Axes`オブジェクトを指定せずに`plt.hist`や`plt.plot`によって図を毎回ゼロから作り直していました。`Figure`や`Axes`の作成はmatplotlibによるグラフ作成の中でもコストのかかる部分です。初めはこれが遅い原因かなと思ったのですが、試してみたら実は__`hist`メソッドが毎回大量の長方形を作っている部分がボトルネック__でした。調べて見ると、[公式マニュアル](https://matplotlib.org/2.0.2/users/path_tutorial.html#compound-paths)にもヒストグラムのアニメーションは長方形の高さを直接更新していくのが効率的でよいと書いてあります。Examplesにも[Animated histogram](https://matplotlib.org/gallery/animation/histogram.html)として同じ方法が紹介されています。この記事の結論は「matplotlibでヒストグラムのアニメーションを作るときは[Animated histogram](https://matplotlib.org/gallery/animation/histogram.html)の方法を(Juliaで使うなら翻訳して)丸写しする」となりますが、比較のためもとの黒木さんのJupyter notebookでの方法と`Figure`及び`Axes`を先に作って`ax.hist`でフレームを更新する方法も示しておきます。

# 環境
Julia 0.6.2をJupyter notebookで実行しています。`PyPlot`で呼んでいるmatplotlibや`PyCall`で呼んでいるnumpyはJuliaの環境構築の際に追加されたものをそのまま使っています。

# データ準備
おそらく財産分布に関するデータを作っているのだと思います。最終的にヒストグラムにしたいデータは`Asset`という65536行41列のArrayです。

```julia
using Distributions

function asset_update!(asset, taxfunc, niters)
    n = endof(asset)
    ur = UnitRange{Int32}(1:n)
    for iter in 1:niters
        k = rand(ur)
        l = rand(ur)
        d = taxfunc(asset[k])
        asset[k] -= d
        asset[l] += d
    end
end

function sim!(asset, taxfunc; niters=5*10^5, nframes=100, thin=niters÷nframes)
    n = endof(asset)
    T = eltype(asset)
    niters = thin * nframes
    Asset = zeros(T, n, nframes+1) #Array{T, 2}(n, nframes+1)
    Asset[:,1] .= asset
    for t in 1:nframes
        asset_update!(asset, taxfunc, thin)
        Asset[:,t+1] .= asset
    end
    Asset
end

function Elog(Asset)
    n, L = size(Asset)
    ElogA = Array{eltype(Asset),1}(L)
    for l in 1:L
        ElogA[l] = mean(log(Asset[k,l]) for k in 1:n)
    end
    ElogA
end

mutable struct ConstantRateTax{T<:AbstractFloat}; rate::T; end
ConstantRateTax(; rate = 0.2) = ConstantRateTax(rate)
(f::ConstantRateTax)(x) = f.rate * x

srand(2018)
n = 2^16
asset = rand(Exponential(1.0), n)
taxfunc = ConstantRateTax(rate=0.05)
nframes=40
@time Asset = sim!(asset, taxfunc, niters=4*10^6, nframes=nframes);
```

# 作成されるアニメーション
どの方法でも出力されるアニメーションは以下のようなものです。うじょうじょ動いているヒストグラムの描画が高速化の鍵になります。
![hoge.gif](https://qiita-image-store.s3.amazonaws.com/0/222849/5d653724-154c-a3f0-a569-2641f14ad714.gif)


# フレーム更新毎に`Figure`と`Axes`をゼロから作り直す方法
Pyplotベースの方法なので暗黙的に`Figure`と`Axes`の作成が行われてます[^1]。`@time`の測定結果はコードの下に示しました。75秒かかります。

```julia
import PyPlot

PyPlot.rc("figure", titlesize=10)
PyPlot.rc("lines", linewidth=1)
PyPlot.rc("lines", markersize=5)
PyPlot.rc("axes", grid=true)
PyPlot.rc("axes", labelsize=10)
PyPlot.rc("grid", linestyle=":")
PyPlot.rc("xtick", labelsize=8)
PyPlot.rc("ytick", labelsize=8)
PyPlot.rc("legend", fontsize=8)

using PyCall
anim = pyimport("matplotlib.animation")

using Distributions

function PyPlot_plot_Asset(t, Asset, ElogA; xmax=4.0, ymax=2.0)
    n, L = size(Asset)
    bins = max(41, Int(fld(sqrt(n), 2)))
    
    PyPlot.clf()

    PyPlot.subplot(121)
    PyPlot.plt[:hist](@view(Asset[:,t+1]), normed=true, bins=bins, alpha=0.5, range=(0,xmax), label="asset dist")
    fitdist = fit(Gamma, @view(Asset[:,t+1]))
    x = linspace(0, xmax, 201)
    PyPlot.plot(x, pdf.(fitdist, x), label="Gamma dist")
    PyPlot.ylim(0,ymax)
    PyPlot.legend()
    PyPlot.title("t = $t")
    
    PyPlot.subplot(122)
    PyPlot.plot(0:t, ElogA[1:t+1])
    PyPlot.xlim(0, L-1)
    PyPlot.ylim(minimum(ElogA), 0)
    PyPlot.title("mean of log(asset)")
    
    PyPlot.tight_layout()
    PyPlot.plot()
end

function PyPlot_makehtml(Asset)
    L = size(Asset)[2]-1
    ElogA = Elog(Asset)
    
    # t = 0～L について1フレーム分の画像を作成
    update(t) = PyPlot_plot_Asset(t, Asset, ElogA)
    
    # アニメーション作成開始
    fig = PyPlot.figure(figsize=(7, 2.5))
    
    # t の動かし方
    frames = [0;0;0;0;0;0;0;0;0:L;L;L;L;L;L;L;L;L;L:-1:0]
    
    # アニメ―ションオブジェクトの作成
    myanim = anim[:FuncAnimation](fig, update, frames=frames, interval=100)
    
    # Javascript animation の作成
    myanim_html = myanim[:to_jshtml]()
    
    # 無駄な表示を削除
    PyPlot.clf()
    
    return myanim_html
end

@time myanim_html = PyPlot_makehtml(Asset)
display("text/html", myanim_html)
```
```:output
74.597071 seconds (13.09 M allocations: 289.686 MiB, 0.28% gc time)
```


# 先に`Figure`と`Axes`を作っておき更新毎に`hist`を実行する方法
更新用関数の外でglobalとして`fig`と`axes`を宣言して、更新用関数の中で`axes[1][:hist]`を呼んでヒストグラムを描画しています。少なくとも`Axes`が何度も新規作成されることはなくなりましたが、__`hist`メソッドによってフレーム更新毎に大量の長方形が作られ`Axes`に`add_patch`され表示領域やtickなども全て設定し直されている__ため、あまり時間短縮の効果は見られず65秒ほどかかっています。

```julia
import PyPlot
PyPlot.rc("figure", titlesize=10)
PyPlot.rc("lines", linewidth=1)
PyPlot.rc("lines", markersize=5)
PyPlot.rc("axes", grid=true)
PyPlot.rc("axes", labelsize=10)
PyPlot.rc("grid", linestyle=":")
PyPlot.rc("xtick", labelsize=8)
PyPlot.rc("ytick", labelsize=8)
PyPlot.rc("legend", fontsize=8)

using PyCall
anim = pyimport("matplotlib.animation")

using Distributions

global fig
global axes

# あらかじめFigureとAxesを用意
fig, axes = PyPlot.subplots(nrows=1, ncols=2, figsize=(7, 2.5))

function PyPlot_plot_Asset(t, Asset, ElogA; xmax=4.0, ymax=2.0)
    global axes
    
    n, L = size(Asset)
    bins = max(41, Int(fld(sqrt(n), 2)))
    
    # 前のフレームのヒストグラムと線を消す
    axes[1][:patches] = []
    axes[1][:lines] = []
    axes[2][:lines] = []
    
    axes[1][:hist](@view(Asset[:,t+1]), normed=true, bins=bins, alpha=0.5, range=(0,xmax), label="asset dist", color="C0")
    #  colorをカラーサイクルの最初の色に指定しないと毎回色が変わってしまう
    fitdist = fit(Gamma, @view(Asset[:,t+1]))
    x = linspace(0, xmax, 201)
    axes[1][:plot](x, pdf.(fitdist, x), label="Gamma dist", color="C1")
    axes[1][:set](ylim=(0, ymax), title="t = $t")
    
    if !isa(axes[1][:get_legend](), PyObject)
        axes[1][:legend]()
    end
        
    axes[2][:plot](0:t, ElogA[1:t+1], color="C0")
    axes[2][:set](xlim=(0, L-1), ylim=(minimum(ElogA), 0), title="mean of log(asset)")
    
    PyPlot.tight_layout()
    
end

function PyPlot_makehtml(Asset)
    global fig
    L = size(Asset)[2]-1
    ElogA = Elog(Asset)
    
    # t = 0～L について1フレーム分の画像を作成
    update(t) = PyPlot_plot_Asset(t, Asset, ElogA)
    
    # アニメーション作成開始
    
    # t の動かし方
    frames = [0;0;0;0;0;0;0;0;0:L;L;L;L;L;L;L;L;L;L:-1:0]
    
    # アニメ―ションオブジェクトの作成
    myanim = anim[:FuncAnimation](fig, update, frames=frames, interval=100)
    
    # Javascript animation の作成
    myanim_html = myanim[:to_jshtml]()

    # 無駄な表示を削除
    PyPlot.clf()
    
    return myanim_html
end

@time myanim_html = PyPlot_makehtml(Asset)
display("text/html", myanim_html)
```
```:output
66.233906 seconds (11.02 M allocations: 180.657 MiB, 0.06% gc time)
```

# 先に手作りした雛形の長方形(`Patch`)の頂点を更新する方法
[Animated histogram](https://matplotlib.org/gallery/animation/histogram.html)の一部（arrayのindexingやメソッドの表記、文字列のクオーテーションなど）をJulia向けに翻訳しただけです。ヒストグラム化には`PyCall`経由で`numpy.histogram`を使っています[^2]。また、長方形の更新用のglobal変数が増えています。あらかじめ描画してある長方形の高さを更新しているだけなので一気に効率が良くなり15秒ですみます。

```julia
import PyPlot
PyPlot.rc("figure", titlesize=10)
PyPlot.rc("lines", linewidth=1)
PyPlot.rc("lines", markersize=5)
PyPlot.rc("axes", grid=true)
PyPlot.rc("axes", labelsize=10)
PyPlot.rc("grid", linestyle=":")
PyPlot.rc("xtick", labelsize=8)
PyPlot.rc("ytick", labelsize=8)
PyPlot.rc("legend", fontsize=8)

using PyCall
anim = pyimport("matplotlib.animation")
patches = pyimport("matplotlib.patches")
path = pyimport("matplotlib.path")
np = pyimport("numpy")

using Distributions

global fig
global axes
global patch
global verts
global top
global bottom

# データをヒストグラム化し、FigureにFuncAnimationで更新するための大量の長方形を描画しておく
# binsとxmaxはここではひとまず手動で指定
weights, edges = np[:histogram](Asset[:, 1], bins=128, density=true, range=(0, 4))

left = edges[1:end-1]
right = edges[2:end]
bottom = zeros(size(left))
top = bottom + weights
nrects = size(left)[1]
nverts = nrects * (1 + 3 + 1)
verts = zeros((nverts, 2))

codes = ones(Int8, nverts) * path[:Path][:LINETO]
codes[1:5:end] = path[:Path][:MOVETO]
codes[5:5:end] = path[:Path][:CLOSEPOLY] 

verts[1:5:end, 1] = left
verts[1:5:end, 2] = bottom
verts[2:5:end, 1] = left
verts[2:5:end, 2] = top
verts[3:5:end, 1] = right
verts[3:5:end, 2] = top
verts[4:5:end, 1] = right
verts[4:5:end, 2] = bottom
# 各長方形5行目の頂点(vertex)はCLOSEPOLYコードを指定するために必要だが実際の座標は無視されるので(left, bottom)を代入する必要はない

fig, axes= PyPlot.subplots(nrows=1, ncols=2, figsize=(7, 2.5), subplot_kw=Dict([("title", "dummy for tight_layout")]))
# tight_layoutは実行時に描画されているArtistを全て表示するため、先にsubplotのダミーのタイトルを描画しておく

barpath = path[:Path](verts, codes)
patch = patches[:PathPatch](barpath, edgecolor="none", alpha=0.5, label="asset dist")
# 上で指定したラベルがlegendとして使われる。edgecolorをなくすために文字列で"none"と指定する

axes[1][:add_patch](patch)

PyPlot.tight_layout()
# アニメーション中に描画領域は変わらないので全ての構成要素を揃えた後に一度tight_layoutすればよい

function PyPlot_plot_Asset(t, Asset, ElogA; xmax=4.0, ymax=2.0)
    global axes
    global patch
    global verts
    global top
    global bottom
    
    n, L = size(Asset)
    bins = max(41, Int(fld(sqrt(n), 2)))
    
    axes[1][:lines] = []
    axes[2][:lines] = []
    
    # 各フレームで使うAssetのデータをnumpyでヒストグラム化、長方形の高さを更新する
    weights, edges = np[:histogram](Asset[:, t+1], bins=bins, density=true, range=(0, xmax))
    top = bottom + weights
    verts[2:5:end, 2] = top
    verts[3:5:end, 2] = top
    
    fitdist = fit(Gamma, @view(Asset[:,t+1]))
    x = linspace(0, xmax, 201)
    axes[1][:plot](x, pdf.(fitdist, x), label="Gamma dist", color="C1")
    # "C1"で使用しているmatplotlibスタイルでのカラーサイクルの2番目の色を指定
    # 同じsubplotに何度もplotしているのでcolorを指定しないとカラーサイクルに従い毎回色が変わってしまう
    axes[1][:set](ylim=(0, ymax), title="t = $t")
    
    # legendは最初に一度だけ作る。こうしないと更新のたびに凡例が追加されてしまう。
    if !isa(axes[1][:get_legend](), PyObject)
        axes[1][:legend]()
    end
        
    axes[2][:plot](0:t, ElogA[1:t+1], color="C0")
    axes[2][:set](xlim=(0, L-1), ylim=(minimum(ElogA), 0), title="mean of log(asset)")
    
    return [patch;]
end

function PyPlot_makehtml(Asset)
    global fig
    L = size(Asset)[2]-1
    ElogA = Elog(Asset)
    
    # t = 0～L について1フレーム分の画像を作成
    update(t) = PyPlot_plot_Asset(t, Asset, ElogA)
    
    # アニメーション作成開始
    
    # t の動かし方
    frames = [0;0;0;0;0;0;0;0;0:L;L;L;L;L;L;L;L;L;L:-1:0]
    
    # アニメ―ションオブジェクトの作成
    myanim = anim[:FuncAnimation](fig, update, frames=frames, interval=100)
    
    # Javascript animation の作成
    myanim_html = myanim[:to_jshtml]()

    # 無駄な表示を削除
    PyPlot.clf()
    
    return myanim_html
end

@time myanim_html = PyPlot_makehtml(Asset)
display("text/html", myanim_html)
```

```:output
15.712500 seconds (10.89 M allocations: 224.286 MiB, 0.26% gc time)
```

# `ax.hist`で雛形を作り`patch.set_height`メソッドで高さを更新する方法
@ceptree さんが以前interactiveなbarプロットで試したことがあるという方法をhist向けにアレンジしました。20秒ちょいです。非常にシンプルに書けるわりに効果が大きいです。patchの頂点を一気に変更する方法より若干遅いのは長方形`set_height`がforループになっているせいか、あるいは`set_height`が頂点を変更する以外の仕事をしているからでしょう。

```julia
import PyPlot
PyPlot.rc("figure", titlesize=10)
PyPlot.rc("lines", linewidth=1)
PyPlot.rc("lines", markersize=5)
PyPlot.rc("axes", grid=true)
PyPlot.rc("axes", labelsize=10)
PyPlot.rc("grid", linestyle=":")
PyPlot.rc("xtick", labelsize=8)
PyPlot.rc("ytick", labelsize=8)
PyPlot.rc("legend", fontsize=8)

using PyCall
anim = pyimport("matplotlib.animation")
patches = pyimport("matplotlib.patches")
path = pyimport("matplotlib.path")
np = pyimport("numpy")

using Distributions

global fig
global axes
global patches

fig, axes = PyPlot.subplots(nrows=1, ncols=2, figsize=(7, 2.5), subplot_kw=Dict([("title", "dummy for tight_layout")]))
_, _, patches = axes[1][:hist](Asset[:, 1], bins=128, density=true, range=(0, 4), alpha=0.5, label="asset dist")

PyPlot.tight_layout()
# アニメーション中に描画領域は変わらないので全ての構成要素を揃えた後に一度tight_layoutすればよい

function PyPlot_plot_Asset(t, Asset, ElogA; xmax=4.0, ymax=2.0)
    global axes
    global patches
    
    n, L = size(Asset)
    bins = max(41, Int(fld(sqrt(n), 2)))
    
    axes[1][:lines] = []
    axes[2][:lines] = []
    
    # 各フレームで使うAssetのデータをnumpyでヒストグラム化、長方形の高さを更新する
    weights, edges = np[:histogram](Asset[:, t+1], bins=bins, density=true, range=(0, xmax))

    for (patch, weight) in zip(patches, weights)
        patch[:set_height](weight)
    end
    
    fitdist = fit(Gamma, @view(Asset[:,t+1]))
    x = linspace(0, xmax, 201)
    axes[1][:plot](x, pdf.(fitdist, x), label="Gamma dist", color="C1")
    # "C1"で使用しているmatplotlibスタイルでのカラーサイクルの2番目の色を指定
    # 同じsubplotに何度もplotしているのでcolorを指定しないとカラーサイクルに従い毎回色が変わってしまう
    axes[1][:set](ylim=(0, ymax), title="t = $t")
    
    # legendは最初に一度だけ作る。こうしないと更新のたびに凡例が追加されてしまう。
    if !isa(axes[1][:get_legend](), PyObject)
        axes[1][:legend]()
    end
        
    axes[2][:plot](0:t, ElogA[1:t+1], color="C0")
    axes[2][:set](xlim=(0, L-1), ylim=(minimum(ElogA), 0), title="mean of log(asset)")
    

    return [patch;]
    
end

#############################################


function PyPlot_makehtml(Asset)
    global fig
    L = size(Asset)[2]-1
    ElogA = Elog(Asset)
    
    # t = 0～L について1フレーム分の画像を作成
    update(t) = PyPlot_plot_Asset(t, Asset, ElogA)
    
    # アニメーション作成開始
    
    # t の動かし方
    frames = [0;0;0;0;0;0;0;0;0:L;L;L;L;L;L;L;L;L;L:-1:0]
    
    # アニメ―ションオブジェクトの作成
    myanim = anim[:FuncAnimation](fig, update, frames=frames, interval=100)
#     myanim[:save]("hoge.gif", writer="imagemagick")
    # Javascript animation の作成
    myanim_html = myanim[:to_jshtml]()

    # 無駄な表示を削除
    PyPlot.clf()
    
    return myanim_html
end

@time myanim_html = PyPlot_makehtml(Asset)
display("text/html", myanim_html)
```
```:output
 23.607039 seconds (11.37 M allocations: 236.429 MiB, 0.24% gc time)
```

# まとめ
matplotlibでbinの変更がないヒストグラムのアニメーションを作るには、高コストの`hist`メソッドで更新するのではなく、一度`hist`で描画した長方形の高さを更新していく方法をとるのがよいです。しかし、冒頭の[Jupyter nobebook](http://nbviewer.jupyter.org/gist/genkuroki/3d6dbf52a3e52eb7c664bc88632c81d3#GIF%E5%8B%95%E7%94%BB%E4%BD%9C%E6%88%90)の最後で紹介されている__爆速のGRバックエンドを使った例では、単純なヒストグラム描画関数で更新する方法でも10秒程度しかかからないので、わざわざこの記事に書いてあるようなことはする必要はない__とも言えます。ただ、現時点ではGRはあまりドキュメントが充実しておらずmatplotlibほどノウハウも共有されていないので、この記事のような面倒なことが必要だとしても`PyPlot`のほうが早く目的の絵を書けることもあるかもしれません。


[^1]: もしかしたら`Figure`は毎回作成されているわけではないかもしれません。matplotlib.pyplotはこの辺がすべて暗黙的に行われるのでデバッグに苦労することがあります。
[^2]: `StatsBase.jl`の`Histogram`で規格化する方法がわからなかったため。