Notebook

Drawing Julia sets in Julia language¶

はじめに¶

２年前の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集合を計算するアルゴリズムは、以下の通りです。

適当な複素数 $c=a+bi$ を決めて
適当な複素数 $z_0$ に対して
$z_n$ が発散する ( $|z_n|>2$ ) まで $z_{n+1}=z_n^2+c$ を繰り返し
発散した時の $n$ を集める
ステップ 2 から 4 を十分多く繰り返す
集めた $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] のようにすると呼び出せます。

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

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

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

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

ここでは、初期値配列(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のコメント欄へどうぞ。