徒に JuliaでPlot した話

はじめに

タイトルは五七五です。
そんなことより，Juliaを学んでいる時に何をベースにして図を描こうか迷いませんか？私は迷いました。 Plots.jl が一般的のような印象を受けますが，他にも様々な可視化パッケージがあるようです。Julia 熟練者からすれば「自分でパッケージを作れば良いじゃないか」「自分でプルリク送れば万事解決」と思われるかもしれませんが，初学者兼非ガチ勢の私にはハードルが高めです。加えて，本質的な進捗に繋がらないので，データ可視化とそのスキル習得には極力時間をかけたくない…という気持ちもあります。

また，Juliaの全ての分野に精通していそう(?)な @bicycle1885 さんが次のように仰っていたので，勝手ながら引用させていただきます。

Juliaでプロットするには何を使えば良いんですかと聞かれると舌を噛み切って死ぬ。

— (｢・ω・)｢ｶﾞｵｰ (@bicycle1885) 2018年7月15日

Q: Juliaでプロットするには何を使ったら良いですか?
A: プロットするにはmatplotlibをPythonで使ってます。

— (｢・ω・)｢ｶﾞｵｰ (@bicycle1885) 2018年9月20日

私はたまに Python を使う程度で知識もないため，Python で一本化する気にもなれませんでした。これらのことから，いくつかのパッケージで最低限のことだけ習得してから考えよう，という発想に至りました。
ということで本記事では，Juliaの可視化パッケージを複数用い，各々で同一の数値計算結果を表示するアニメーションを作成します。そしてパッケージごとの特徴と個人的な感想を共有することを目的とします。

使用したパッケージ

まずはこの記事で扱われる可視化パッケージ3つについて超簡単に書きます。

• Plots.jl
  Juliaの可視化ツールの中で標準っぽいパッケージ。何かしらの既存パッケージを呼び出し(バックエンド)それを Plots.jl が定めるルールで記述する(フロントエンドコード)ことによって図が生成されます。

• PyPlot.jl
  matplotlib.pyplot。もはやこれ以上の説明は不要。

• GMT.jl
  GMT (Generic Mapping Tools)を基礎としています。GMTは地球科学系の地震火山，気象水象，海岸海洋などの分野で多く使われている印象があります。主な特徴としては，出力形式がPostScriptであること，地理的情報を含むデータを描くのが得意であることが挙げられます。

これら以外にも Plotly.jl, Gadfly.jl などなどがあると思いますが，3パターンで描けて一旦満足してしまったので手を止めました。
plotly については，Plots.jl に plotlyjs バックエンドがあるのでそれを使用すると似たような見た目の図ができると思います（適当）。

出力結果を見てみよう

コードがどうこう論じる前にざっくり図を見た方が良いと思うので，先に gif を貼ります。上記の3つのパッケージを使用して，数値計算結果の可視化を行い，アニメーションを作成しました。 数値計算については，@Bimaterial さんの 2次元波動方程式の解の導出と性質にあるコードを使用しました。この記事のおかげでただ図を作成するだけでも楽しくできました。感謝いたします。
下記のアニメーションは，リンク先の記事のコードでファイル出力されているvの平面分布を各出力時刻について描いたものです。
colorbarのラベルの桁が違うのは(おそらく単位を空間スケールと揃えているからだと思いますが)定数倍すればなんとかなるので目をつぶってください。
元記事にあるアニメーションのほうが美しいのでは？とも思うかもしれませんがそれも見なかったことにして以下の3つを見てください。

• Plots.jl
   
• PyPlot.jl

- GMT.jl   

個人的な使用感

とりあえず図は作成できるようになったので，いろいろいじくってみた個人的な感想をそれぞれ以下に述べます。

• Plots.jl
  一般的なバックエンドであるGRを使うと処理が速いです。Plots 専用の記述ルールなのでその習得コストがかかりますが，引数の順番と構造は matplotlib や MATLAB と似ています。図の体裁を整えるためのキーワード引数をある程度覚える必要があって，覚えるまではソースを読むか公式ドキュメントの Attributes ページをひたすら見返します。
  平面分布図を描く時にややハマったことは，plot(x,y,z)的な関数で x,y が Array{T,1} またはStepRange等の 1-dimension である必要があることです。matplotlib や MATLAB のようにArray{T,2}ではありません。x,yがArray{T,1}しか受け付けないなら曲線座標はどうやって描けばいいんだ？という疑問が湧きますが，まだ Julia でその必要性が生じたことがなく未検証です。
  また，公式ドキュメントや Issues で論じられているように，この記事執筆時点では colorbar (着色のスケール)表示の融通が利かないという難を抱えています。これは論文用の図作成には致命的です。しかし，概形・概観を把握するには上の gif にある通り十分であるように思えます。
  良い点としては，多種多様な図を作成できることです。当然バックエンドの範疇を超える図は描けないので，バックエンドによっては対応していない関数があったり，3次元図が得意でなかったりしますが，コードを編集せずともバックエンドを切り替えるだけで何とかなってしまったり，雰囲気の違う図ができて面白いです。

• PyPlot.jl
  matplotlib を既習しており，特に不満がないならばこれを使っていれば良いような気がします。
  いつからか Python とほぼ同じようにかけるようになりました。以前はax.add_subplot(...) 等を ax[:add_subplot](...) と書く必要がありましたが，今は前者でも大丈夫です。むしろ今では後者で書くと deprecated と出てきます。Python チョットデキル方々にはこれと PyCall で何ら問題ありません。

• GMT.jl
  コードが他のパッケージ使用時より長くなりがちで，やや面倒くさい。ただ GMT はもともと sh 系を使って描くのが一般的であって，それらに比べたら使いやすいと思います。普段 GMT を使っている人が bash で書いて実行していたものを Julia 経由に乗り換える，というのには十分な動機となり得るでしょう。
  加えて，やはり地理的データの描画に関しては他の追随を許さない美しさを誇ります。PostScript 出力という点からも，既に GMT をある程度使いこなしていれば論文用には GMT.jl で図作成しようか，と感じます。

コード例

ということで，上記のアニメーションを生成したコードをそれぞれ載せます。
元の数値計算の記事でファイル出力しているvを，全てVという変数に格納するように編集しました。すなわち，preallocate してから，for文の中で

if i % 10 == 0
    #= 元コード
    open("dat/V$(i).dat","w") do vo
        write(vo,v)
    end
    =#
    V[:,:,floor(Int64,i/10)+1] = v
end

としています。その他の変数はそのまま使用しているので，元のコードの実行後に以下のコードを実行すれば，上で示したような gif ができます。

ちなみに，速さに関しては特にこだわっていないので，「この書き方は遅い…」と思われた場合は適宜読み替えてください。

Plots.jl

この記事作成時のバージョン: v0.23.2

using Printf
using Plots
gr() # plotlyjs(), pyplot() etc.
clibrary(:colorcet)


"""
Make an animation with Plots.jl
"""
function main()
    # the number of time steps
    nt = size(V)[3]
    # draw each step
    anim = @animate for t = 1:nt
    @printf("%4d / %4d\r", t, nt)
    heatmap(x₁ .- x₁⁰, x₂, permutedims(V[:,:,t],[2 1]),
            axis_ratio=:equal,
            c=:bkr, clims=(-10,10),
            xticks=collect(-0.5:0.5:3.5),
            yticks=([-1.0,-0.5,-0.4,0.0],["-1.0","interface","source","0.0"]),
            tickfont=Plots.font("sans-serif",8),
            xguide="horizontal", yguide="depth",
            guidefont=Plots.font("sans-serif",12),
            colorbar_title="velocity",
            size=(850,200),
            )
    end

    # save
    gifname = "out_Plots.gif"
    if isfile(gifname); rm(gifname); end
    gif(anim, gifname, fps=10)
end

# run
main()

PyPlot.jl

記事作成時のバージョン: v2.8.0

# meshgrid
X1 = repeat(collect(x₁), outer=(1,length(x₂)))
X1 = X1 .- x₁⁰
X2 = repeat(collect(x₂), outer=(1,length(x₁)))
X2 = permutedims(X2,[2 1])


using Printf
using PyPlot, PyCall
anim = pyimport("matplotlib.animation")


"""
Arrange axes
"""
function SetAxes!(ax)
    # axis equal
    ax.axis("scaled")
    # horizontal axis
    ax.set_xlabel("horizontal", fontsize=12)
    ax.set_xticks(collect(Float64, -0.5:0.5:3.5))
    # vertical axis
    ax.set_yticks([-1.0,-0.5,-0.4,0.0])
    ax.set_yticklabels(["-1.0","interface","source","0.0"])
    ax.set_ylabel("depth", fontsize=12)
    # return value
    return ax
end


"""
Draw a 2D distribution at a particular time
using PyPlot
"""
function SnapShot(k, ax, X1, X2, V)
    nt = size(V)[3] # the number of time steps
    @printf("%4d / %4d\r",k,nt-1) # display
    # initialize
    ax.clear()

    # draw
    pcm = ax.pcolormesh(X1, X2, V[:,:,k+1], cmap=:twilight)
    pcm.set_clim(-20., 20.)
    SetAxes!(ax)
    # return value
    return pcm
end


"""
Make an animation with PyPlot.jl
"""
function main()
    fig = plt.figure(figsize=(8.0,2.2))
    ax = fig.add_subplot(111)

    # colorbar
    pcm = SnapShot(0, ax, X1, X2, V) # initial step
    cbar = fig.colorbar(pcm, ticks=-20:10:20, shrink=0.75)
    cbar.ax.set_ylabel("velocity", fontsize=10)
    # animation test
    #wanim = anim.FuncAnimation(fig, SnapShot, fargs=(ax, X1, X2, V), interval=50, frames=5)
    # animation
    wanim = anim.FuncAnimation(fig, SnapShot, fargs=(ax, X1, X2, V), interval=100, frames=size(V)[3])
    wanim.save("./out_PyPlot.gif", writer="imagemagick")
end

# run
main()

GMT.jl

記事作成時のバージョン: v0.7.0
関数名が他と競合しがちなのでusing GMT: GMTとしています。
GMT.jl 自体にはアニメーション出力する機能はない(と思う)ので convert と ffmpeg その他を使っています。
convert だけでも gif 作成できるだろ，と思っても怒らないでください。

# meshgrid
X1 = repeat(collect(x₁), outer=(1,length(x₂)))
X1 = X1 .- x₁⁰
X2 = repeat(collect(x₂), outer=(1,length(x₁)))
X2 = permutedims(X2,[2 1])

figdir="fig" # output path
nt = size(V,3) # number of time steps


using Printf
using GMT: GMT

"""
Draw a 2D distribution at a particular time
using GMT
"""
function SnapShot(X1, X2, Vt, cpt::GMT.GMTcpt, dx)
    # options
    proj="X12/3"
    region="-0.5/3.5/-1.0/0.0"
    axes="xa0.5f0.5+lhorizontal ya0.5f0.5+ldepth SW"
    cbD="jBR+w2.9/0.3+o-0.8/0.05"
    cbB="xa5f5+lvelocity"
    # makegrd
    G = GMT.surface([X1[:] X2[:] Vt[:]], R=region, I=dx)
    # draw
    GMT.grdview(G, J=proj, R=region, B=axes, C=cpt, Q="i")
    GMT.colorbar!(J=proj, R=G, B=cbB, D=cbD, C=cpt)
end


"""
Print ps file in all time
"""
function PrintAll()
    # initialize
    if !isdir(figdir)
        mkdir(figdir)
    else
        cd(figdir)
        rm.( filter(x->occursin(r"step\d+\.ps",x), readdir()) )
        cd("../")
    end

    # makecpt
    clim="-10/10"
    # -D option in makecpt() doesn't work at GNT .
    # see https://github.com/GenericMappingTools/GMT.jl/issues/179
    cpt = GMT.gmt("makecpt -Csplit -T$clim -D")
    # fotmat of tick labels
    GMT.gmt("set FONT_ANNOT_PRIMARY 8p")
    GMT.gmt("set FONT_LABEL 10p")

    #k = 50 # test
    # print all
    for k = 1:nt
        @printf("print %4d / %4d\r", k, nt)
        SnapShot(X1, X2, V[:,:,k], cpt, dx)

        pstmp = GMT.fname_out(Dict())[1]
        file = "step"*@sprintf("%04d",k-1)
        mv(pstmp, joinpath(figdir,file*".ps"))
    end
end


"""
convert ps - eps - png - gif
ps2eps, convert (ImageMagick), and ffmpeg are required.
"""
function PS2gif()
    cd(figdir)

    # convert
    for k = 1:nt
        @printf("convert %4d / %4d\r", k, nt)

        file = "step"*@sprintf("%04d",k-1)
        run(`ps2eps $file.ps -f -q`) # ps - eps
        run(`convert -density 300 $file.eps $file.png`) # eps - png
    end

    # make an animation
    run(`ffmpeg -i step%04d.png -vf palettegen palette.png -y`)
    run(`ffmpeg -r 10 -i step%04d.png -i palette.png -filter_complex paletteuse out_GMT.gif -y`)
    run(`rm palette.png`)

    cd("../")
end


PrintAll()
PS2gif()

おわりに

ここまできて言うのもアレですが，PlotsとGMTに関してはバージョンが若く，v1.0.0 に達していません。今後大きく仕様が変更される可能性もあるため，そこは Issues 等をチェックしながら最新を追い求めていく必要があります。1年間以上更新されていない Plot 系のパッケージもあり，何かを習得してもいずれ陳腐なものになってしまうこともあるかもしれません。完全に個人の勝手な予想ですが，本記事で示した3つのパッケージは陳腐化しにくいものだと思い選定しました。他にも良いパッケージがあればこの計算をベンチマークとして同様のことをやってみたいと思います。

##参考
PlotsGallery.ipynb --- @goropikari さん
Plots/GR: グラフ package のおすすめ --- http://www.cas.cmc.osaka-u.ac.jp/~paoon/