LoginSignup
20
15

More than 5 years have passed since last update.

@skotaro(斉藤 耕太郎)

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

Last updated at Posted at 2018-02-04

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

きっかけ

最近の日本語インターネットでのJuliaの盛り上がりに一役買っていると思われる黒木さんのこのJupyter nobebookではJuliaで様々なグラフを描く方法が紹介されていますが、PyPlot (matplotlib) でヒストグラムのアニメーションを描くとむちゃくちゃ遅いという話になっています。

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

環境

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

データ準備

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

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

フレーム更新毎にFigureとAxesをゼロから作り直す方法

Pyplotベースの方法なので暗黙的にFigureとAxesの作成が行われてます1。@timeの測定結果はコードの下に示しました。75秒かかります。

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秒ほどかかっています。

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の一部（arrayのindexingやメソッドの表記、文字列のクオーテーションなど）をJulia向けに翻訳しただけです。ヒストグラム化にはPyCall経由でnumpy.histogramを使っています2。また、長方形の更新用のglobal変数が増えています。あらかじめ描画してある長方形の高さを更新しているだけなので一気に効率が良くなり15秒ですみます。

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が頂点を変更する以外の仕事をしているからでしょう。

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の最後で紹介されている爆速のGRバックエンドを使った例では、単純なヒストグラム描画関数で更新する方法でも10秒程度しかかからないので、わざわざこの記事に書いてあるようなことはする必要はないとも言えます。ただ、現時点ではGRはあまりドキュメントが充実しておらずmatplotlibほどノウハウも共有されていないので、この記事のような面倒なことが必要だとしてもPyPlotのほうが早く目的の絵を書けることもあるかもしれません。

  1. もしかしたらFigureは毎回作成されているわけではないかもしれません。matplotlib.pyplotはこの辺がすべて暗黙的に行われるのでデバッグに苦労することがあります。 ↩

  2. StatsBase.jlのHistogramで規格化する方法がわからなかったため。 ↩

20
15
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
20
15