フラクタル 第3回

今度は前回のマンデルブロ集合の変則的な形として、試してみたものです。


今回は

z(n+1) = z(n)^3 + c

という漸化式による収束性について複素平面の図を作ってみました。

できたのは下のような左右対称な図です。

各ｃ（複素数）の点で、ｚの発散しなかった領域が黒い部分です。緑部分は発散領域です。
この場合もフラクタルのような図形になりました。
図の中央が原点で、左右対称な図形となりました。
範囲は実数軸、虚数軸ともに-1.5から1.5までです。

データ点は約14M 個生成して作成しました。
漸化式の繰り返し回数は10回です。

今回はpythonからデータを生成して、MySQLを使ってデータを管理しましたが、30分程度かかりました。
今回の計算結果はWEBを検索しても簡単に見つからなかったので、ちょっとした課題としていいかもしれません。



作った関数はこんなかんじです。

def makeMandelbrot3( itrn, n=100, xrange=None, yrange=None ):
if xrange == None:
xrange = (-2.5, 2.5)
if yrange == None:
yrange = (-2.5, 2.5)
c = np.random.uniform( xrange[0], xrange[1], n) + 1j*np.random.uniform( yrange[0], yrange[1],n)

z = 0
for i in range(itrn):
z = z**3 + c
return c, z


itrnは繰り返し回数で、データ数はnです。
戻り値はcとzなので、
この関数を使用するときは

c, z = makeMandelbrot3( 10, 1000 )

のようにします。
これで繰り返し10回のデータcとzが1000点分生成されます。
メモリ確保の問題さえなければ、nの値を大きくすればその分データが生成されます。

bool

ブール型の演算について

どの言語もそうですが、整数型や浮動小数点型の演算はだいたい同じですが、ブール型の演算は少しずつやり方が異なっているような気がします。

pythonでの扱い方について記録しておきます。

a = True
b = False
a and b

とすると
False

となります。
a+b
とすると
 1
となり、
a+a
とすると
 2
になります。

つまり、ブール型はTrueが1、Falseが0の値になっていると考えて構わないと思います。

それで、今度はnumpy.arrayでの扱いです。
aa = numpy.array( [ a, a, b] )
とすると
array([ True,  True, False], dtype=bool)

となり、bool型として扱われていることがわかります。
で、今度は
aa and aa
とするとエラーになります。

実行した結果のエラーメッセージは下のようになっていました。
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: The truth value of an array with more than one element is ambiguous.
 Use a.any() or a.all()

numpy.arrayは他の型ではそれぞれの要素に関して個別に演算するので、ブール型もそうかと思うとそうでもないようです。
で、今度は

aa = numpy.array( [True, True, False], dtype=bool)
bb = numpy.array( [False, True, False], dtype=bool)
aa + bb
aa * bb
とすると、それぞれ下のようになりました。

>>> aa+bb
array([ True,  True, False], dtype=bool)
>>> aa*bb
array([False,  True, False], dtype=bool)

つまり、要素個別に演算する場合は+と*でandとorの代用する必要があるということです。
andとor,notは単独のTrue, Falseといったブール値にしか使えません。


notで実現していた、反転をさせる場合はどうするかというと、
~aa
とします。
結果は
array([False, False,  True], dtype=bool)
となります。

ブールのarrayは結構有用です。
たとえば、

x = numpy.random.randint( 0, 10, 12 )
c = x>5
print x[c]

というように、条件にあったものだけ抽出する際に使います。
cは
array([ True, False,  True, False,  True, False,  True,  True,  True,
       False, False, False], dtype=bool)
といった風にbool型のarrayになっています。

フラクタル 第２回

以前、マンデルブロ集合の図を作成しました。
今回pythonはあまり関係ないです。

マンデルブロ集合の画像を線分データのベクター画像(svg)として変換しました。
PDF


そのベクターデータを元に彫刻してみました。

工作関連は別なブログに載せてありますので、詳しくはそちらを見て下さい。

フラクタル

pythonで複素数の練習を兼ねて
フラクタル図形を作ってみました。
マンデルブロ集合というものです。


フラクタルが何かについては専門書または他のサイトなどで調べてみてください。
マンデルブロ集合の定義は簡単な複素数漸化式で記述されます。
つまりマンデルブロ図形というのはある複素数cについて、

z(n+1) = z(n)*z(n) + c
z(0) = 0

とした場合、n→∞で発散するか否かを色分けして複素平面に表示したものです。

複素数平面の点をランダムに選び、n回漸化式を解いてz(n)を作ります。
その時の、
log |z(n)|
の値を画素の輝度値としました。
forループを2重に使うと時間がかかるので、numpy.arrayで一気に計算します。

m = 10000
c = numpy.random.uniform( -2.5,2.5,m ) + 1j * numpy.random.uniform( -2.5,2.5,m )
z = numpy.zeros( m, dtype=numpy.complex )
for i in range(16):
     z = z**2 + c

と言う感じです。

opencvを使って画像化します。
しかし、案外時間が掛かるので、zの漸化式計算と、画像化を分けて処理しました。

Mandelbrot set

30Mポイントのサンプリングをして、そこから画像化しましたが、埋め切れていません。画像を拡大してみるとポツポツと黒い点が見えます。
実は細くつながっています。が、データ点数（サンプリング）で埋め切れていないのと画像化の解像度が荒いため、つながっているようには見えません。

複素数

python numpyでの複素数の扱い方についてです。


複素数は型としてはnumpy.complexと定義されています。
これを利用してみます。

つまり、

numpy.array( [3, 2, 1], dtype=numpy.complex )

とすると

　array([ 3.+0.j,  2.+0.j,  1.+0.j])

というふうに複素数配列となります。


複素数の演算もできます
a = numpy.array( [3, 2, 1], dtype=numpy.complex )

a ** 2

とすると2乗します。

実行結果は

　array([ 9.+0.j,  4.+0.j,  1.+0.j])

となります。

指数計算はそのまま掛けて計算されます。

つまり、python風な数式で表すと

z=x+jy

とすると

z**2 = (x+jy) * (x+jy) = x**2 - y**2 + j 2*x*y

となります。

複素数には色々と実数にはなかった計算方式が出てくるので注意が必要です。

実数部は

numpy.real(z)

虚数部は

numpy.imag(z)

とします。

これらの関数はzがarrayでもOKです。その場合、出力もarrayになります。

複素共役は

numpy.conj(z)

とします。

複素共役が何かわからない場合は検索してみて下さい。

誰かがわかりやすく説明してくれていると思います。

絶対値はnumpy.abs(z)

で計算できます。

複素数の乱数発生は簡単にはできませんでした。

一様分布の場合

z = numpy.random.uniform( -5, 5, 10 ) + 1j * numpy.random.uniform( -5, 5, 10 )

としました。

結果は

array([-4.48658518-4.19169279j,  2.39542899+4.87030019j,

        4.03412469-1.00727465j,  1.72313031+0.99345025j,

        0.68053779-0.385337j  ,  2.08333535-0.32240792j,

       -0.65638447-0.27364433j,  1.36415791-1.89435328j,

        0.19786534-4.51657661j,  2.44853396+3.04850488j])

となりました。

もっとスマートな方法があれば、是非ご教示下さい。

svgwriteのエラー

pythonでSVGファイルを生成する必要があったので、
svgwriteというモジュールをインストールしてみました。

インストールは簡単にできます。
pip install svgwrite
です。


が
サンプルを実行してみようとすると
エラー…

import下だけなのにエラーが出ました。
エラーメッセージを貼り付けると下のようにxml関連らしいです。

>>> import svgwrite
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "c:\Python26\Lib\site-packages\svgwrite\__init__.py", line 51, in <module
>
    from svgwrite.drawing import Drawing
  File "c:\Python26\Lib\site-packages\svgwrite\drawing.py", line 29, in <module>

    from svgwrite.container import SVG, Defs
  File "c:\Python26\Lib\site-packages\svgwrite\container.py", line 27, in <modul
e>
    from svgwrite.base import BaseElement
  File "c:\Python26\Lib\site-packages\svgwrite\base.py", line 12, in <module>
    from svgwrite.etree import etree
  File "c:\Python26\Lib\site-packages\svgwrite\etree.py", line 33, in <module>
    original_serialize_xml = etree._serialize_xml
AttributeError: 'module' object has no attribute '_serialize_xml'


下のサイトを参考に少し修正しました。
https://bitbucket.org/mozman/svgwrite/pull-request/3/fix-etreepy/diff


pythonインストールフォルダに
Lib/site-packages
というフォルダがあります。
この下にいわゆるライブラリがインストールされます。

インストールされていれば、
この下にsvgwriteというフォルダが有ると思います。
その中にetree.py
というファイルがあるので、

33行目
original_serialize_xml = etree._serialize_xml

を削除して、最終行に移動します。
現状のエラーでは循環呼び出しをしてしまったらしいです。

これで、ひとまずimportは出来ました。

使い方などはまた次の機会にレポートします。