Drawing Julia sets in Julia language

はじめに

2年前の12月13日付けで公開されていた 「Juliaでジュリア集合を描く」 http://qiita.com/anmitsu/items/0f85a6427b4224a2d56a という記事を、最近目にしました。

Juliaの勉強の題材としてよさげに思えたので読み始めましたが、掲載していたプログラムはそのままでは動きませんし、改善の余地がいくつもありそうでした。

そこで、著者に無断で元記事を添削した顛末を書いてみたのが、以下のエッセイです。

実行環境の紹介

執筆に用いた実行環境を紹介します。

OSは Windows 8.1 64bit です。

PyPlot モジュール を使ってグラフを表示するので、Python と matplotlib を予めインストールする必要があります。

ananonda distribution を用いて、Python 64bit版と数値計算パッケージをインストールしました。numpy, matplotlib, jupyter が一度に入るから楽ですね。( 私は、ここの記事を参考にしました → Python数値計算環境Anacondaの導入 )

> python
Python 2.7.10 |Anaconda 2.4.0 (64-bit)|(default, Nov 7 2015, 13:18:40) [MSC v.1500 64 bit (AMD64)] on win32

>>> import matplotlib
>>> matplotlib.__version__
"1.4.3"
>>>

Julia 64bit 版を、本家のバイナリからインストールしました。→ http://julialang.org/downloads/

In [1]:
versioninfo()
Julia Version 0.4.1
Commit cbe1bee* (2015-11-08 10:33 UTC)
Platform Info:
  System: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i5-5257U CPU @ 2.70GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

Juliaの PyPlotパッケージのインストールは コマンドライン (REPL) から

Pkg.add("PyPlot")

です。→ https://github.com/stevengj/PyPlot.jl#installation

元記事を試してみる

「Juliaでジュリア集合を描く」 http://qiita.com/anmitsu/items/0f85a6427b4224a2d56a (以下、元記事といいます) に掲載されたプログラムを試してみましょう。

まず Julia集合を計算する関数です。

In [2]:
# ジュリア集合を計算
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
Out[2]:
julia_set (generic function with 1 method)

次は等高線をプロットする関数です。 /// なぜ、下を Code セルにしないのかという疑問を先に解消したい方には こちらを参照。

using PyPlot

# 等高線 (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

最後に、メインルーチンです。

In [3]:
# 代表的なジュリア集合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])
        @show minimum(z), maximum(z)
        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
Out[3]:
demo_julia_set (generic function with 1 method)

ここまで構文エラーはありません。 では実行。

>>> demo_julia_set()
LoadError: type PyObject has no field MaxNLocator
while loading In[5], in expression starting on line 1

 in getindex at C:\Users\ ___ \.julia\v0.4\PyCall\src\PyCall.jl:239

 0.014893 seconds (4 allocations: 5.722 MB, 30.87% gc time)
(minimum(z),maximum(z)) = (0.0,0.09090909090909091)

あれれ、エラーです。MaxNLocator が見つかりませんって。 これは、matplotlib の関数 matplotlib.ticker.MaxNLocator です。

あとで説明しますが、show_contourf()を、以下のように修正すると、MaxNLocatorを使うことができます。

In [4]:
function show_contourf(z)
    levels = matplotlib[:ticker][:MaxNLocator](nbins=50)[:tick_values](minimum(z), maximum(z))
    plt[:contourf]( z, levels=levels, cmap=get_cmap("gray_r"))
end
Out[4]:
show_contourf (generic function with 1 method)

修正して、いざ実行。

In [5]:
using PyPlot

demo_julia_set()
  0.004566 seconds (4 allocations: 5.722 MB)
(minimum(z),maximum(z)) = (0.0,0.09090909090909091)
  0.009774 seconds (4 allocations: 5.722 MB)
(minimum(z),maximum(z)) = (0.0,0.09090909090909091)
  0.009198 seconds (6 allocations: 5.722 MB, 50.67% gc time)
(minimum(z),maximum(z)) = (0.0,0.09090909090909091)
  0.004528 seconds (4 allocations: 5.722 MB)
(minimum(z),maximum(z)) = (0.0,0.09090909090909091)
  0.004535 seconds (4 allocations: 5.722 MB)
(minimum(z),maximum(z)) = (0.0,0.09090909090909091)
 
 0.004558 seconds (4 allocations: 5.722 MB)
(minimum(z),maximum(z)) = (0.0,0.09090909090909091)

あれれ、枠しか出力されません。 z の最大最小の値が、c の値によらず一定なのもおかしいです。どんどんメモリを確保していく割には、計算時間が短すぎる気がします。

添削開始

二重ループ内の break

元記事によると、Julia集合を計算するアルゴリズムは、以下の通りです。

  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$ をプロットする

// 数式が正しく表示されることを確かめたくて、わざわざ書いてみました。

上のアルゴリズムを実装しているのは、関数 demo_julia_set() の以下の部分です。

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

上のループを見ると

for i = 1:length(z), n = 1:500
    end

二重ループを 一文で書けるのですね。しびれます。

しかし、ループ内で breakしたら、二重ループ全体から抜けてしまうのではないかな? 短いプログラムを書いて試してみましょう。

In [6]:
# i=j=0 # a1
for i=1:3, j=1:3
    @show i,j     
    if j==2 
        break
    end
end
# @show i, j # a2
(i,j) = (1,1)
(i,j) = (1,2)

一文に二重ループを書くと、後ろに書いた変数が先に増加することが分かります。 心配した通り break すると、二重ループ全体から抜けてしまいました。ですから、今回のプログラムには不適です。

( ちなみに、ループを抜けた時点の変数 i, j の値を知ろうとして、ループ直後 (上プログラム片 #a2の位置)に @show i,j を書いてもエラーが出てしまいます。変数 i, j のスコープ (生存範囲)がループ内だからです。そうしたかったら、ループの始まる前に ( #a1 の位置 )、変数 i, j を作っておかないといけません )

さて、ここの二重ループは一文にまとめず、for文を入れ子にすべきです。

初期値配列 z[i] の寸法を小さくして、試してみます。

In [7]:
zc = Complex( 0.28500 , 0.0100 )
    z = [ Complex(x,y) for x in linspace(-0.5, 0.5, 3 ), y in linspace(-0.5, 0.5, 3 ) ]
    julia = zeros(Float64, size(z))
    nmax = 500
    for i = eachindex(z)
        julia[i] = nmax
        @time for n = 1:500
            if abs2(z[i]) > 4.0
                julia[i] = n
                break
            end
            z[i] = z[i] * z[i] + zc
        end
    end
    julia
 
Out[7]:
3x3 Array{Float64,2}:
 76.0  11.0  20.0
 22.0  20.0  22.0
 20.0  11.0  76.0
 0.012142 seconds (1.53 k allocations: 71.213 KB)
  0.000011 seconds (128 allocations: 3.656 KB)
  0.000010 seconds (116 allocations: 3.313 KB)
  0.000004 seconds (62 allocations: 1.766 KB)
  0.000005 seconds (116 allocations: 3.313 KB)
  0.000003 seconds (62 allocations: 1.766 KB)
  0.000004 seconds (116 allocations: 3.313 KB)
  0.000005 seconds (128 allocations: 3.656 KB)
  0.000016 seconds (452 allocations: 12.938 KB)

i=1 で終わることなく、全ての i について計算されてました。よいですね。

ついでに、いくつか気づいたところを書き換えました。

■ Python の list comprehension と同じように

z = [ Complex(x,y) for x in linspace(-0.5, 0.5, 3 ), y in linspace(-0.5, 0.5, 3 ) ]

とすると、2次元配列を作れます。 複素数は x + y * im と書くより、Complex(x,y)と書いたほうが効率がよいとマニュアルに書いてあります。 → http://docs.julialang.org/en/release-0.4/manual/complex-and-rational-numbers/#complex-numbers

size(z)で、配列 z の次元・寸法がとれます。つまり、 julia = zeros(Float64, size(z)) の文で、配列 z と同じ次元・寸法を持つ Float64 配列を確保できるわけです。

eachindex(z) は、配列 z の全要素について、効率的に添字を巡ります。 配列の次元・寸法が変わっても、同じ構文で書けるのもよいですね。

abs(zz) > 2.0 から abs2(zz) > 4.0 へ。 → http://docs.julialang.org/en/release-0.4/manual/performance-tips/#tweaks

z^2 から z * z へ。べき乗は k = 2 の決め打ちだから、自乗がよいと思います。

初期値配列の確保・初期化

パラメータ cについて Julia集合を計算する関数 julia_set(c) では、初期値配列 z[i] の確保・初期化を関数冒頭で毎回行っています。

初期値配列z[i] は、どのパラメータ c についても共通・不変のはずですから、この関数を呼出す側で一度だけ確保・初期化すればよいですね (参考 Pre-allocating outputs )。 んで、初期値配列を毎回確保・初期化しているのは z[i]を書き換えているからです。Julia条件を満足する繰返し数 n が必要なだけなら、z[i] を書き換えたり、z を保存する必要はありません。

以上の方針で書き換えてみましょう。

In [8]:
zc = Complex( 0.28500 , 0.0100 )
    const zs = [ Complex(x,y) for x in linspace(-0.5, 0.5, 3 ), y in linspace(-0.5, 0.5, 3 ) ]
    nmax = 500
    julia = zeros(Float64, size(zs))
    julia[:] = nmax
    for i = eachindex(zs)
        zz=zs[i]
        @time for n = 1:nmax
            if abs2(zz) > 4.0
                julia[i]=n
                break
            end
            zz = zz*zz + zc
        end
    end
    julia
 
Out[8]:
3x3 Array{Float64,2}:
 76.0  11.0  20.0
 22.0  20.0  22.0
 20.0  11.0  76.0
 0.000025 seconds (303 allocations: 8.281 KB)
  0.000011 seconds (87 allocations: 2.375 KB)
  0.000009 seconds (79 allocations: 2.156 KB)
  0.000005 seconds (43 allocations: 1.172 KB)
  0.000005 seconds (79 allocations: 2.156 KB)
  0.000003 seconds (43 allocations: 1.172 KB)
  0.000006 seconds (79 allocations: 2.156 KB)
  0.000006 seconds (87 allocations: 2.375 KB)
  0.000018 seconds (303 allocations: 8.281 KB)

const は、一度設定した値を書き換えないことを示すキーワードです。 実際 zs の値は書き換えませんから、const zs に反したという警告は出ません。いいですね。

毎回の関数呼び出しで確保されるメモリの寸法は、元の 2/3 くらいに少なくなりました。 漸化式

zz = zz * zz + zc

の計算 (複素数演算) で、メモリ確保が起こるようです。

等高線プロット

matplotllib APIの呼出し

PyPlotのドキュメント 曰く

Only the currently documented matplotlib.pyplot API is exported.

すなわち、matplotlib.pyplot に記述されたAPIのみが輸出されているとのことです。

であれば、matplotlib.pyplot.contourf は、contourf() でそのまま呼び出せるはずなのですが、どうもうまくいきません。 plt[:contourf] のようにすると呼び出せます。

MaxNLocatormatplotlib.ticker クラスのメソッドですから、そのまま呼び出せません。matplotlib[:ticker][:MaxNLocator] のように呼び出します。

プロット範囲を任意にする

初期値集合は、必ずしも正方形でないので、任意の x, y について、等高線を描けるとよいですね。

配列 x, y, z の三つ組で contourf を呼び出す場合、配列 x, y, z の次元・寸法を同じにします。 python では mgrid 関数を用いて、各軸の座標値から2次元配列を作るのが通例でしょう。

ここでは、初期値配列(2次元) zs が計算してありますから、これからその実部・虚部の値が入った2次元配列を作りましょう。どちらも one-line です。

zs = [ Complex(x,y) for x in linspace(-1.5, 1.5, 600 ),  y in linspace(-1.0, 1.0, 400 ) ]
    xs = real(zs)
    ys = imag(zs)

最終版

以上を踏まえて、プログラムを完成しましょう。

まず、Julia集合を計算する関数です。 初期値 配列 zs[i]、結果配列 ju[i]は、どちらも呼出し側で領域確保してから呼び出します。

In [9]:
function compute_julia_set(zc,ju,zs,nmax)
    @assert size(ju)==size(zs) "-- ju: size error"
    @inbounds ju[:] = 1.0 / (10.0+nmax)
    @time for i=eachindex(zs)
        zz=zs[i]
        for n=1:nmax
            if abs2(zz) > 4.0
                ju[i]=1.0 / (10.0+n)
                break
            end
            zz=zz*zz + zc
        end
    end
end
Out[9]:
compute_julia_set (generic function with 1 method)

上を呼び出すメイン関数です。 等高線プロットも含めてしまいます。subplot の使い方にもご注目。

In [10]:
function do_compute_julia_set()
    const zs=[ Complex(x,y) for 
        x in linspace(-1.5, 1.5, 600 ), 
        y in linspace(-1.0, 1.0, 400 ) ]
    const xs=real(zs)
    const ys=imag(zs)
    @time ju=zeros(Float64, size(zs))
    fig = matplotlib[:pyplot][:gcf]()
    fig[:set_size_inches](18.5, 10.5)
    nmax=500
    result=[]
    nc=1
    for zc in [ 
            Complex( 1.0 - (1.0+sqrt(5.0))/2.0, 0.0),
            Complex( 0.28500,  0.0000 ),
            Complex( 0.28500,  0.0100 ),
            Complex(-0.70176, -0.3842 ),
            Complex(-0.83500, -0.2321 ),
            Complex(-0.80000,  0.1560 ) ]
        @time compute_julia_set( zc, ju, zs, nmax )
        # println( minimum(ju) )
        # println( maximum(ju) )
        push!( result, deepcopy(ju) )
        levels = matplotlib[:ticker][:MaxNLocator](nbins=50)[:tick_values](minimum(ju), maximum(ju))
        plt[:subplot](2,3,nc)
        plt[:contourf](xs,ys,ju, levels=levels, cmap=get_cmap("gray_r"))
        plt[:title](string("c = ", zc ))
        plt[:xticks]([])
        plt[:yticks]([])
        nc = nc + 1
    end
    # println( typeof(result) )
    # println( size(result) )
    plt[:subplots_adjust](left=0.05, bottom=0.05, right=0.95, top=0.95, wspace=0.1, hspace=0.1)
    # savefig("result6.jpg")
end
Out[10]:
do_compute_julia_set (generic function with 1 method)

上で、変数 result には、各パラメータ c に対するJulia集合が追記されます。単にグラフ表示するだけなら不要ですが、あとで使いたくなるかもしれないしね。

では実行しましょう。

In [11]:
using PyPlot

do_compute_julia_set()
  0.000473 seconds (4 allocations: 1.831 MB)
  0.180337 seconds
  0.180966 seconds (293 allocations: 24.031 KB)
  0.008609 seconds
  0.009236 seconds (305 allocations: 15.313 KB)
  0.016115 seconds
  0.016696 seconds (295 allocations: 14.656 KB)
  0.008826 seconds
  0.009471 seconds (305 allocations: 15.313 KB)
  0.008703 seconds
  0.009290 seconds (298 allocations: 14.766 KB)
 
 0.036813 seconds
  0.037502 seconds (302 allocations: 15.203 KB)

大成功です。

Juliaのコンソールでは関数を再定義できない

ここからは編集後記に相当するものです。

Ruby, Pythonなどのスクリプト言語では、アイデアをその場で試せるのが売りです。コマンドラインから打ち込んだ文は即座に実行され、結果が得られます。 ユーザ関数を再定義すれば、古い定義は消去され、新しい定義で上書きされます。Julia でも当然そうだと、私は完全に思い込んでいました。

へんな挙動に気づいたのは、元記事を試してみる の節で、関数 show_contourf() の不具合に気づき、書き直そうとしていたときでした。この関数を上書きしても、動作が同じなのです。Jupyter上の不具合なのかと想像して、Juliaのコマンドライン (REPL)で作業してみても、やはり動作は同じです。困りましたねえ。

日本語記事を少し当たってみましたが、この件に関する記述はなかなか見つかりません。しかたがないので、Julia function redefinition といったキーワードでぐぐってみると「関数再定義ができない」と騒いでいる人が何人か見つかり、それに対する回答は「Juliaの仕様です」。(えー)

Julia manualには長いこと undocumentedだったそうですが、Julia 0.4 manual の FAQに載っています。

32.2.1 How do I delete an object in memory?

Julia does not have an analog of MATLAB’s clear function; once a name is defined in a Julia session (technically, in module Main), it is always present.

元記事を試してみる の記述で、不具合のある関数show_contourf() を Codeセルにしないで、地の文章 (markdown)で書いたのは、そういうわけなのでした。

「仕様」であれば仕方ないのですが、この仕様が続く限りは、注意しないといけませんね。

まとめ

この記事のエッセンスをまとめておきましょう。

  • Juliaでジュリア集合を描く、を動作可能としました。
  • 一文で書いた二重ループから break すると、二重ループ全体から抜けてしまいます。
  • 配列確保を少なくするように工夫しましょう。
  • PyPlot で matplotlib APIを呼ぶ方法を確認しました。
  • Juliaのコンソールでは関数を再定義できません。
  • プログラムと実行結果を同時に見せるのに、Jupyter - nbviewer は便利です。

以上の文章は、 『Qiita 決定版?「Juliaでジュリア集合を描く」』http://qiita.com/tenfu2tea/items/6bcbdd7586ea070cc25d の本文です。コメントは Qiitaのコメント欄へどうぞ。