Juliaでジュリア集合を描く

Julia とは科学技術計算に特化したプログラミング言語です。

• MATLAB のような行列計算や数値計算
• C や Fortran に匹敵する速度
• Ruby や Python のような動的型付け
• Lisp のようなマクロ機能

上記のような各言語の特徴的な機能を詰め込んだ欲張り言語です。
詳細は公式ページを参照して下さい。特に他言語とのベンチマークの比較は圧巻です。

さて，Julia という名前を聞くと，フラクタル図形でよく知られているジュリア集合を思い浮かべる方も多いと思います。
本ページでは Julia でジュリア集合を計算してグラフに出力する方法を紹介します。

以下は本ページのコードで得られたジュリア集合の図です。

環境構築

紹介するコードを実行するには以下の環境が必要です。

Julia
http://julialang.org

PyPlot
https://github.com/stevengj/PyPlot.jl

matplotlib
http://matplotlib.org

上記の環境をすでに導入されている方はこのセクションを読み飛ばして下さい。

なお，私の実行環境は以下のようになります。

• OS: Ubuntu 12.04 64bit
• Julia 0.3.0-prerelease+396
• matplotlib 1.3.1
• Python 2.7.3

Julia

Julia の導入方法は http://julialang.org/downloads/ を参照して下さい。

私は Ubuntu 12.04 でナイトリー版の Julia を利用しているので

$ sudo add-apt-repository ppa:staticfloat/julianightlies
$ sudo add-apt-repository ppa:staticfloat/julia-deps
$ sudo apt-get update
$ sudo apt-get install julia

でインストールしました。

PyPlot

Julia でグラフを描画する方法はいくつかありますが本ページでは PyPlot を使います。
PyPlot は，Python のグラフ描画モジュール matplotlib を Julia から利用するためのパッケージです。
以下のコマンドで PyPlot をインストールできます。

julia> Pkg.add("PyPlot")

PyPlot を用いる理由は私自身が Python に慣れているのと比較的高速かつインタラクティブにグラフを描画してくれるためです。

matplotlib

PyPlot を使うためには当然 matplotlib もインストールされている必要があります。
matplotlib の最新版は PyPI から入手可能です。
pip を使えば簡単にインストールできます。

$ pip install matplotlib

なお，事前に Python の科学技術計算ライブラリである numpy や scipy も導入しておく必要があると思います。

ジュリア集合を描く

Wikipedia によれば

ジュリア集合 (Julia set) とは

z_{n+1} = z_n^2 + c

において複素数 $c$ を固定した場合に，この漸化式が無限大に発散しないような初期値 $z_0$ を与える集合である。

定義や特性などは 英語版 Wikipedia が詳しいです。

一見難しそうに見えますがジュリア集合をプロットするアルゴリズムは結構単純です。

1. 適当な複素数 $c = a + bi$ を決めて
2. 適当な複素数 $z_0$ に対して
3. $z_n$ が発散する ($|z_n| > 2$) まで $z_{n+1} = z_n^2 + c$ を繰り返し
4. 発散した時の $n$ を集める
5. ステップ 2 から 4 を十分多く繰り返す
6. 集めた $n$ をプロットする

多分コードを見たほうが分かりやすいと思います。

using PyPlot

# ジュリア集合を計算
function julia_set(c::Complex128)
    SIZE = 500
    x = linspace(-1.5, 1.5, SIZE)
    y = linspace(-1.0, 1.0, SIZE)
    z = zeros(Complex128, SIZE, SIZE)
    for j = 1:SIZE, i = 1:SIZE
        z[i, j] = x[j] + y[i] * im
    end
    k = 2.0
    julia = zeros(Float64, SIZE, SIZE)
    for i = 1:length(z), n = 1:500
        z[i] = z[i]^k + c
        if abs(z[i]) > 2.0
            julia[i] = 1.0 / (10.0+n)
            break
        end
    end
    julia
end

# 等高線 (fill) で色分けプロット
function show_contourf(z)
    levels = plt.MaxNLocator(nbins=50)[:tick_values](minimum(z), maximum(z))
    contourf(z, levels=levels, cmap=get_cmap("gray_r"))
end

# 代表的なジュリア集合6つをプロット
function demo_julia_set()
    c = [complex128(1 - (1+sqrt(5))/2),
          0.28500 + 0.0000im,
          0.28500 + 0.0100im,
         -0.70176 - 0.3842im,
         -0.83500 - 0.2321im,
         -0.80000 + 0.1560im,]
    
    clf()
    for i = 1:length(c)
        @time z = julia_set(c[i])
        subplot(2,3,i)
        show_contourf(z)
        title(string("c = ", c[i]))
        xticks([])
        yticks([])
    end
    subplots_adjust(left=0.05, bottom=0.05, right=0.95, top=0.95,
                    wspace=0.1, hspace=0.1)
end

変数 SIZE を大きくすれば解像度が上がり，計算時間が2乗に比例して増加します。
変数 x と y の範囲を変えれば特定の場所をズームインできます。
組み込みマクロ @time でジュリア集合の計算時間を測定します。
関数 show_contourf のカラーマップを変えれば別の色分けでプロットできます。

関数 demo_julia_set でページ冒頭の図をプロットできます。
私の環境 Thinkpad T510 Intel Core i7 CPU M 620 2.67GHz では，1セット当たり最小で約0.05秒，最大で1秒程度で計算出来ました。

まとめ

Julia でジュリア集合を計算し PyPlot を用いてプロットしました。
動的型付けにも関わらず結構高速に計算できたと思います。

Julia は標準で C と Fortran の機能が呼び出せます。
加えて，PyPlot のバックエンドである PyCall パッケージを用いると Python の機能も利用できます。
このように，Julia は新しい言語でライブラリもまだ少なめですが，既存の言語とのインターフェースがかなり充実してるので，数値計算だけでなく様々なニーズに対応できる汎用的な言語だと感じました。

ジュリア集合やマンデルブロ集合は並列計算しやすいらしいです。
Julia は並列計算もサポートしているので，興味がある方はこれらの集合を並列計算で求めてみてはいかがでしょうか。

参考