PyQtでお手軽GUI開発♪―――は可能だったか？第10回　マンデルブロ編

　これまで主にPyQtを使ってGUIを作る方法に関していろいろ書いてきました。しかしどんなにそれが使いやすくてどんなに美麗なUIができようとも、基本性能がお粗末では実用にはなりません。その基本性能の中でも最も重要な要素と言えるのがプログラムの実行速度でしょう。

　そこでこれよりPyQtの速度はどの程度なのか、また遅ければそれをどのくらい改善できる見込みがあるのかということについて、マンデルブロ集合描画をテーマにお送りしましょう。

マンデルブロ集合

　プログラミングに興味のある方なら多分どこかで下のような画像を見たことがあるのではないでしょうか。

　1985年10月、その頃まだ日経のついてないサイエンス誌にA.K.デュードニーという人の“コンピューターレクリエーション”という名コラムがありましたが、そこで「コンピューターの顕微鏡で，数学で最も複雑な図形を拡大観察する」というタイトルで紹介されて一躍有名になったのがこのマンデルブロ集合です。

　この妖しい美しさに満ちた画像に脳をやられてこちらの世界に来てしまった方も結構いるんじゃないかと思いますが、そうなった最大の理由がこの画像を得るための原理が驚くほど簡単だったからなのではないでしょうか。

　そもそもマンデルブロ集合とは何かというと……

a, b を任意の実数として

X₀ = 0

Y₀ = 0

X_n+1 = X_n² - Y_n² + a

Y_n+1 = 2X_nY_n + b

こういう漸化式があった場合、ベクトル（X_n, Y_n）の絶対値が無限大に発散しない初期値（a, b）の集合をマンデルブロ集合と呼ぶ。

　式の内容は高校で習う程度でどこも難しくありません。そんな単純なルールから上のような摩訶不思議な画像が現れてくるのです。

　さてここで問題になるのが集合の定義の絶対値が無限大に発散しないというところです。こんな条件のため、厳密にはある点がマンデルブロ集合に属するかどうかは一般に有限時間内には決定できません。

　絶対値が2を越えたら発散することは証明されているので、ベクトル長が2を越えたらそこで計算を打ち切ってしまえるのですが、まだ越えてない場合にあと何回計算したらいいかは分からないのです。あと10回で越えるかもしれないし、100万回かかるかもしれないし、永久に越えないかもしれません。

　そこで実際上は近似値として、計算回数が一定回数を越えたらもう発散しないだろうとみなして打ち切ってしまうしかありませんが、綺麗な画像を鑑賞するのが目的なら実はそれで十分なのです。

　というのが、上記の画像には美しい模様がありますが、この色づけはベクトルの絶対値が2を越えるまでの回数を元に色づけしたものだからです。すなわち綺麗なのはマンデルブロ集合ではないところで、マンデルブロ集合そのものは中央にある変な形の真っ黒い穴なのです。

　それはともかく今の話は平面上にある一点についてです。従ってこれを画像にしようと思ったら初期値（a, b）を少しずつ変えていって、各地点ごとにそんな計算をしなければなりません。

　例えば上の画像を得るためには、画像サイズが500*500で各ピクセルごとに打ち切り回数4000回の計算をしているので、単純なかけ算なら10億回の計算をしていることになりますが―――色がついているということは打ち切り回数よりは少ない回数で計算が終わっているということなので、おおむね数億回の計算が行われているわけです。

　要するにちょっと凝った画像を描こうと思うととてつもない計算回数が必要になってくるのです。

　このマンデルブロ集合は大昔に、初めて買ったPC9801VX2（80286 8MHz）でQuickBasicを使って描いた記憶がありますが―――多分サイズが250×250くらいで打ち切り回数は数百回クラスの画像一枚を描くのに1時間くらいはかかっていて、本気で80287を買おうかどうしようか悩んだような記憶もそこはかとあったりします。

　あれからPCの性能は多分数千倍以上の規模で進歩していますが、ではこいつをPyQtで書いてみたらどういうことになるでしょうか？

プログラムの基本構造

　さて、まず以下がマンデルブロ集合生成の基本関数です。

　関数mandelは座標（a, b）を左下として幅 iw、高さ ih ピクセルの画像用のデータを作成します。1ピクセルのサイズが pw で与えられているので、画像の幅と高さが iw*pw ，ih*pw になります。

なおpwは1/2ⁿになるように値を与えています。こうすることでいろいろ計算誤差が減るためです。

　計算結果、すなわちベクトル長が2を越えるまでのループ回数を data という二次元リストに保管します。こうしておけば後でそれを元にいろんな色づけができるからです。

def mandel(data, a, b, pw, iw, ih, cmax):
    """マンデルブロ集合を生成する
    
    data    結果を保管する二次元リスト
    a,b     計算開始位置の座標(float)
    pw      ピクセルのサイズ(float)
    iw,ih   計算する領域の幅と高さ(pixel)
    cmax    計算を打ち切る回数
    """
    for j in range(ih):
        b0 = b + pw*(ih - j)                        #スクリーン座標を数学の座標に変換
        for i in range(iw):
            a0 = a + pw*i

            x0, y0 = a0, b0                         #第1回目の計算でこうなるので
            for n in range(cmax):
                x1 = x0**2 - y0**2 + a0
                y1 = 2*x0*y0 + b0
                if math.sqrt(x1**2 + y1**2) > 2:    #ベクトルの長さが2を越えたら抜ける
                    break
                x0, y0 = x1, y1

            data[j][i] = n                          #結果は配列に保管しておく

　関数は見たとおりとてもシンプルです。

ちなみに、トップにある画像のパラメータは、a=-1.9990958532821423, b=2.763181328191422e-08, pw=1/2⁴²=2.2737367544323206e-13, iw=ih= 500, cmax=4000です。

◆計算部分の関数化

　ということなのですが、いろいろなアルゴリズムを比較するに当たっての試行錯誤の結果、実際には下記のように漸化式の計算部分を関数化して呼び出すようにしています。

def calc0(a0, b0, cmax):
    """マンデルブロ漸化式計算ルーチン"""

    x0, y0 = a0, b0
    for n in range(cmax):
        x1 = x0**2 - y0**2 + a0
        y1 = 2*x0*y0 + b0
        if math.sqrt(x1**2 + y1**2) > 2:
            break
        x0 , y0 = x1, y1
    return n


def mandel_p(data, a, b, pw, iw, ih, cmax, calc):
    """マンデルブロ集合のデータ作成"""

    for j in range(ih):
        b0 = b + pw*(ih - j)
        for i in range(iw):
            a0 = a + pw*i

            data[j][i] = calc(a0, b0, cmax)

　というのも、やってみると表示が遅いのでプログレスバーを出したりなどいろいろせねばならず、mandel_p関数内が相当複雑になってしまい、それをコピペしていると後からの修正が大変になってしまったためです。

　また計算速度に関してはほぼcalcの速度に依存しているので、そこだけを抜き出して見せた方が分かりやすいということもあります。

　そこで今後はアルゴリズムに関してはcalcだけを抜き出して説明することにします。

◆画像の色づけ方法

　あともう一点、計算結果を画像にするには発散が確定するまでのループ回数を色にしなければなりません。

　そのためにはループ回数に対応した色を割り当てたカラーパレットを作成して、得られたdataの値を見ながらQImageにsetPixcelColorしていけばいいわけです。

　キモはそのパレットの作成方法ですが、今回は色をHSVで表してH（色相）をぐるぐる循環させながらV（明度）を徐々に下げていってcmaxのところで0にするという方法にしました。多分あちこちでわりとよく使われてる方法です。

昔はこのHSVをRGBに変換するのに結構苦労しました。しかしQtのQColorにはsetHsvという関数があるので今では一発です。
なおパラメータhstartやhdiffでHの初期値や増分を調整できるようにしていますが、これをちょっと変えるだけで画像のイメージががらっと変わって、それだけでしばらく遊べるからなのです。

def makepalette(self, cmax, hstart=0, hdiff=2, s=255):
    """カラーパレットを作成する

    cmax    mandelループ最大カウント数
    hstart  hの初期値
    hdiff   1カウントごとのhの増分
    s       彩度
    """
    self.colortable = [None] * cmax
    hc = 360 // hdiff
    for i in range(cmax):
        h = ((i+hstart) % hc) * hdiff   #hはぐるぐる回る
        v = 255 - (255 * i // cmax)     #Vはだんだん暗くなる
        c = Qg.QColor()
        c.setHsv(h, s, v)
        self.colortable[i] = c

　以上見ての通り、プログラムの基本はすごくシンプルです。こんなものからあんな画像が生まれるというのがまさに神秘といっていいかもしれません。

実行時間計測デコレータ

　というわけでこれからマンデルブロ集合の描画速度に関していろいろ検討していくわけですが、そのためにはまず描画速度をきちんと計測する必要があります。

　ここで実はSpyderには便利なことに最初からプロファイラというものがついています。これはその上でプログラムを実行すると、全ての関数の実行時間を自動的に表にしてくれるものです。これを使えば確かに上記の関数の実行速度を調べられるのですが、このプロファイラはアプリ全体が遅い場合のボトルネックを探すのを主目的に作られているようで、単一の関数をちょこちょこ計測するのにはあまり向いていません。

　そこで自前の速度計測をすることになりますが、Pythonにはtimeモジュールというものがあって、その中のperf_counter関数を使えば簡単です。これは高分解能の経過時間をfloat値で返してくれる物で、開始時と終了時に時間を取得して差分を取ればクロックの分解能ぐらいの精度で経過時間が分かります。

　しかし今後作る関数に一々そんなのを仕込んでいるのも面倒です。そんな場合の強い味方がデコレータです。

　他所のPythonソースを見ると、関数の前に@hogehogeといったような修飾子がついていることがあるでしょう。それがデコレータで、ある関数に後付けで機能追加をする機能です。

　下のfunctime関数はある関数を実行したら実行時間をミリ秒単位で表示してくれるデコレータです。

def functime(func):                                 #“func”という関数を引数に取る
    """関数の実行時間を表示するデコレータ"""

    def dfunc(*arg):                                #関数内で定義された関数“dfunc”
        t = time.perf_counter()                     #現在の時間を取得

        r = func(*arg)                              #functime引数の“func”を実行

        t = (time.perf_counter() - t) * 1000        #実行後時間を取得しミリ秒単位に変換
        print('描画時間: {:.1f} ms'.format(t))      #結果を表示する
        return r                                    #dfuncの結果はfuncの結果

    return dfunc                                    #functimeの結果として“dfunc”を返す

　これを使うときには以下のように適用したい関数の前に@functimeと書くだけです。

@functime
def mandel_p(data, a, b, pw, iw, ih, cmax, calc):

　こうしておけば今後mandel_p関数が実行されたら同時に実行時間が表示されます。

　例えば次のようなパラメータを与えて実行してみます。

mandel_p(data, -1.4765625, -0.9765625, 0.0078125, 250, 250, 500, calc0)

　するとこのような画像が表示されて……

　コンソールには以下のようなメッセージが出力されます。

描画時間: 23352.3 ms

　すなわち上記の画像データを計算するのに23秒ちょっとかかったことが分かります。

◆デコレータの簡単な説明

　ここで軽くデコレータについて筆者の分かる限りの解説をしておきましょう。

　デコレータの公式の説明のところにはいろいろややこしいことが書かれていますが、細かいことを抜きにすればデコレータの作り方は難しくありません。

　まず理解しておくべきことは、デコレータとは関数を引数にとって関数を結果として返す関数である、ということです。Pythonの場合、関数も含めて大体何でもオブジェクトなのでそんな関数でも作ることができます。

　デコレータを作るには前述のfunctime関数のような書き方が一種のテンプレになります。functimeの個々の行が何をやっているかは、横のコメントを見てもらえば理解は難しくないでしょうが、以下もう一度まとめてみると……

まずfunctimeという関数があってfuncという関数引数を取っています。もちろん関数引数というのはコールバック関数など、どこででも使われるテクニックです。
functime関数内でdfuncという関数が定義されています。これもPascalなど多くの言語でごく普通にできていることがPythonでもできるというだけです。
サブ関数dfunc内で、functimeの引数に取られた関数funcが実行されて、その前後で時間を計測しています。すなわちこのdfuncが実行されるとfuncも実行されて、その実行時間が表示されるわけです。

　さて一番分からないのが、

このfunctime関数がその結果としてdfunc関数を返している

　ということでしょう。これは何を意味するのでしょうか？

　それは要するにこうすることでこの関数がデコレータとして機能できるようになるということです。すなわちこのような“関数を結果として返す関数hogehoge”があれば、Pythonという言語にはある関数 f の前に@hogehogeと書いてやることで、関数 f を hogehoge( f ) の結果として戻ってくる関数に差し替えてくれるという機能があるのです。

　上の例ならば mandel_p が functime(mandel_p) の結果として返される dfunc に差し替わってしまうわけです。

そのために普通は引数の数などを一致させておかなければなりませんが、Pythonの場合*argなんて書き方もできるのでこうしておけばそのあたりもあまり気にしなくて構いません。

　そしてデコレータとは大体こんなものだということさえ分かれば、あとは例えば実行前にパラメータを前処理したり、結果に修正を加えたり、関数そのものをMySpecialFunctionに差し替えてしまったりといったこともできるのが分かるでしょう。自分で色々作ってみるのももう難しくはないと思います。

もちろん思慮深い人ならこんな説明では納得いかないでしょう。@hoge云々は仕様なのでそういうものだとしても、例えばdfuncって関数の中で定義されてるけどスコープはどうなってるのかとか、同じデコレータをいろんな関数に適用したら変な干渉はないかとか、考え出したら夜も寝られなくなりそうです。そのへんの説明が公式ドキュメントにある“関数は第一級オブジェクト云々”だと思うのですが、私もよく分かってるわけではないし、ともかくこれで動いてるんで気にしなくていいでしょう（おいこら！）
どうしても気になる人は Python デコレータでググればいろいろ見つかります。

アルゴリズムの効率化

　さて前節で計測したマンデルブロ計算時間ですが、結果は23秒強とあまり芳しいものではありませんでした（それでもPC9801VX2の頃に比べれば格段の差ですがｗ）

　しかしこの例は、マンデルブロ集合の定義を一番素朴に書いたものです。多分見た瞬間においおいと突っ込んだ人も多かったことでしょう。

　まずダメな点ですが、2乗の計算に x0**2 とやっています。しかしこの演算子 ** は累乗数が負だったり整数以外の場合にも対応している汎用的な演算子なので、内部処理が複雑で計算が遅くなるのです。単に2乗を得たい場合なら x0*x0 とかけ算で書いた方が速くなるということが大昔から知られています。

　またベクトル（x1,y1）の絶対値計算に、2乗の和を取ったあとに平方根（sqrt）を計算して2より大きいかどうかチェックしていますが、これは2乗の和が4より大きいという条件にしても同値です。

　というわけで以上の点を修正したのが以下のcalc1です。

def calc1(a0, b0, cmax):
    """累乗を積に変更し、平方根も２乗の和に変更"""

    x0, y0 = a0, b0
    for n in range(cmax):
        x1 = x0*x0 - y0*y0 + a0
        y1 = 2*x0*y0 + b0
        if x1*x1 + y1*y1 > 4:
            break
        x0 , y0 = x1, y1
    return n

　このように変えて同じパラメータで計算してみたら……

描画時間: 7472.1 ms

　おーっと！かなり速くなります！いきなり7秒台とか、前に比べて3倍以上の高速化です！この調子ならもっと何とかなるかもしれません！

　そこで上の関数をよく見てみると―――2乗を計算しているところが、x1を計算しているところとベクトル長を計算しているif文の2カ所にあります。しかしx1の計算のときには前回のベクトル長計算結果を流用することが可能です。

　そんな感じで修正したのが以下のcalc2です。

def calc2(a0, b0, cmax):
    """2乗を2回計算しているのを1回にする"""

    x0, y0 = a0, b0
    xx0, yy0 = a0*a0, b0*b0         #最初の1回
    for n in range(cmax):
        x1 = xx0 - yy0 + a0         #保存していた値を使う
        y1 = 2*x0*y0 + b0
        xx0, yy0 = x1*x1, y1*y1     #ループ中は2乗の計算はここだけ
        if xx0 + yy0 > 4:
            break
        x0 , y0 = x1, y1
    return n

実行結果は……

描画時間: 6633.9 ms

　今度は6秒台と、さっきほどではありませんが少し速くなりました。最初に比べれば3.5倍くらいです。

　しかし、式を眺めていてもこれ以上はどうしようもなさそうです。

複素数でやってみると……

　ところでPythonの公式サイトのFAQにプログラムが遅すぎます。どうしたら速くなりますか？などという項目があります。そのなかに“適切なデータ構造を使ってください。組み込み型や collections を調べてください”というヒントがありました。

　実はこのマンデルブロ集合というのは本来は複素力学系という数学のカテゴリに属する話でした。すなわち「複素数の空間上での関数の反復適用によって定義される力学系の研究」から生まれてきた物なのだそうです。

　この複素数とは a, b が任意の実数、i を虚数（2乗したら-1になるという謎の数）としたら a + bi の形で書き表すことのできる数です。

　このうち実数だけでできた a の項を実数部、虚数との積になっている b の項を虚数部といいます。

　この数は実数部 a を横軸、虚数部 bi を縦軸とした平面上の点で表せますが、その平面を複素平面と呼び、複素数の絶対値とは中心から（a, bi）までの距離、すなわちaとbの2乗の和の平方根で定義されます。

　マンデルブロ集合とは元々この複素平面上で以下のように定義されたものです。

任意の複素数Cがあった場合

Z₀ = 0

Z_n+1 = Z_n² + C

こういう漸化式があった場合、Z_nの絶対値が無限大に発散しない初期値Cの集合をマンデルブロ集合と呼ぶ。

　さてこの複素数の計算なのですが、実は i を単なる変数とみなして (a + bx)² のような普通の整式のように計算するだけです。ただし計算中に i² が出てきたら-1で置き換えてやる必要があるというだけです。

　そこで C = a + bi ， Z_n = X_n + Y_ni として上記の漸化式の部分を書き直してみると……

Z_n+1 = Z_n² + C

　　　= (X_n + Y_ni)² + a + bi

　　　= X_n² + a + (2X_nY_n + b)i + Y_n² i²

　　　= X_n² + a + (2X_nY_n + b)i - Y_n²　　　　← i² を -1 に置き換える

　　　= X_n² - Y_n² + a + (2X_nY_n + b)i

　こんな感じでZ_n+1の実数部が（X_n² - Y_n² + a）、虚数部が（2X_nY_n + b）となり、これを通常の座標平面と見立ててやれば今までの漸化式と同じになるわけです。

　―――要するに何が言いたいかというと、Pythonには組み込みの型に複素数型というのがあるんでこれでやってみようということでした。

　そうやって書いてみたのが以下のcalc3ですが、もう滅茶苦茶シンプルです。

def calc3(a0, b0, cmax):
    """複素数で計算してみる"""

    C = complex(a0, b0)
    Z0 = C
    for n in range(cmax):
        Z1 = Z0*Z0 + C
        if Z1.real*Z1.real + Z1.imag*Z1.imag > 4:
            break
        Z0 = Z1
    return n

　その結果は……

描画時間: 5918.5 ms

　うーむ……それでも5秒台は出ていて確かに高速化はしていますが、ややこしい前振りの割にはちょっと残念な結果だったでしょうか？

やっぱりPythonではダメなのか？

　というわけで、標準のPythonでできることはこのくらいまでのようです。

　だとしたら―――そもそも今の速度とはどの程度のものなのでしょうか？もし他の言語を使って書いてもこの程度というのであれば、もちろんPythonを責めるわけにはいきません。

　そこで使い慣れたDelphi6を使って書いてみます。アルゴリズムは以下のようにcalc2で使った物を利用しました。

Pascalではタプルが使えないんでこういうのでも結構長くなります。

function calc(a0, b0: double): integer;
var x0, y0, x1, y1, xx0, yy0: double;
    n: integer;
begin
    x0 := a0;
    y0 := b0;
    xx0 := a0*a0;
    yy0 := b0*b0;
    for n:=1 to CMAX do begin
        x1 := xx0 - yy0 + a0;
        y1 := 2*x0*y0 + b0;
        xx0 := x1*x1;
        yy0 := y1*y1;
        if (xx0 + yy0) > 4 then break;
        x0 := x1;
        y0 := y1;
    end;
    result := n;
end;

さてその結果は……

経過時間 = 350.06 (ms)

　………………

　…………

　……

ケタが違いますうっっｗ

　同じcalc2に対しては19倍、複素数のcalc3に比べてさえ17倍くらい速いです！実行したら一瞬で画面が現れます。Delphi6って出たのは今から15年以上は前だと思いますが、そんなロートルにこんなにぶっちぎられてていいんでしょうか？ Pythonじゃこの程度しかできないとなるといろいろ考え直さなければならないのですが―――でもこのあたりで第1回目から読んでくれた人なら、張ってあった伏線を思い出してくれるかも知れません。

　そうなのです。公式のFAQにもこそっと書かれていますが、Pythonには強い味方がいたのです。次回はその助っ人Cythonと、そしてもう一人Numbaというのがどの程度の実力を持っていたかについてです。

なお、何か間違いなどがありましたらこちらの方にコメントしてください（質問等でも構いませんが、多分答える能力があまりないのではないかと思います）

2017-04-24

PyQtでお手軽GUI開発♪―――は可能だったか？ 第10回 マンデルブロ編