cmathとmatplotlibでMandelbrot集合を描画する。

cmathに思いをはせて

cmathモジュールの活用例を求めて色々と考えた後、まずはMandelbrot集合を書いてみようと思い立った。 Mandelbrot集合の定義はWikipediaを参照してほしい。

マンデルブロ集合 - Wikipedia

簡単に述べると、ある複素数列が収束する初期値の集合である。

最初は定義をよく確認せず実装してしまいうまくいかず、以下のサイトのお世話になった。 以前からMandelbrot集合の存在は知っていたが定義をきちんと理解したのは今回が初めてかもしれない。

aoking.hatenablog.jp

方針として、複素数列が収束するかを判定する函数、float版rangeとなるジェネレータを実装し、Mandelbrot集合を計算、matplotlibで描画、いう流れである。

import cmath
import decimal
import itertools

import matplotlib.pyplot as plt


def mandelbrot(c):
    z = complex(10, 10)
    for _ in range(1000):
        z = z**2 + c
        if abs(z) > 2:
            return False
    else:
        return True

def drange(x, y, jump):
    while x < y:
        yield float(x)
        x += decimal.Decimal(jump)

xs = []
ys = []
for x, y in itertools.product(drange(-2, 2, 0.01), drange(-2, 2, 0.01)):
    if mandelbrot(complex(x, y)):
        xs.append(x)
        ys.append(y)

plt.axis('off')
plt.plot(xs, ys, ".", color='k')
plt.savefig("mandelbrot.png")

実行結果は以下の通りである。

あまりきれいではない。

解析学っぽいことに触れると、自分の中に位相空間の知識が欠落していることにいやでも気づかされる。 しばらく時間ができそうな気配がするので、位相空間という借金を返済しようかと考えている。

cmathモジュールのジレンマ

無何有の郷へ

Pythonには複素数を計算するための函数を収めた標準モジュールcmathが存在する。

9.3. cmath — 複素数のための数学関数 — Python 3.6.1 ドキュメント

ドキュメントにはこのような一節がある。

mathと同じような関数が選ばれていますが、全く同じではないので注意してください。機能を二つのモジュールに分けているのは、複素数に興味がなかったり、もしかすると複素数とは何かすら知らないようなユーザがいるからです。そういった人たちはむしろ、math.sqrt(-1)が複素数を返すよりも例外を送出してほしいと考えます。また、cmathで定義されている関数は、たとえ結果が実数で表現可能な場合 (虚数部がゼロの複素数) でも、常に複素数を返すので注意してください。

複素数に興味がないプログラマ、複素数を全く知らないプログラマはそれぞれ存在するのか気になるところである。 日々の業務において実数向けの数学モジュールmathすら使う機会がない。最近必要に迫られてmath.cosを利用したぐらいである。

一方で、複素数に興味があるプログラマ、複素数体が代数的閉体であることを知るプログラマも存在する。 複素数に関するプログラミングにニーズがあるとすれば、電気回路、信号理論、制御理論などのFourier変換やLaplace変換が必要となるエンジニア、正則函数と戯れたい数学者がもつだろう。 しかし、彼らが行いたいことを満たすにはcmathは力不足である。彼らが必要なのはNumPyやSciPy、SymPyなどの科学計算ライブラリ、Mathematicaなどの専用のソフトウェアである。

帯に短し襷に長し。 興味がないプログラマには不要の長物であり、必要とするプログラマには機能が足りない。

荘子の逍遥遊篇に次の一説がある。荘子が恵子に無用の用を説く場面である。

今、子に五石の瓠あり。何ぞこれ以て大樽を為りて江や湖に浮かぶことを慮えずして、其の瓠落して容るる所なきを憂うるや。すなわち夫子には猶ほ蓬れたる心有り」と。

一見役に立たないcmathをいかに役に立たせるかが第77回Python Mini Hack-a-thonのテーマであった。

複素数の性質

実数と複素数は似ている要素もあるが違う要素もある。 似ている要素としては、体をなす（四則演算ができる）こと、位相を入れることができる（収束や連続という概念が定義できる）、がある。 位相を入れることができるのでcmathには指数函数や三角函数が存在するのである。

一方で、決定的に異なる要素としては全順序を定義できない、つまり自然な大小関係を定義できないことである。 Pythonもこの事実を受けて、複素数同士を比較すると例外が送出される。

a = complex(1, 2)
b = complex(10, 10)
a < b

TypeError: unorderable types: complex() < complex()

異なる性質を持つので、実数と複素数を分けて考える必要がある。

複素数と画像変換

今回は、複素数が実数平面上の点として表せることに着目して、ある種の画像変換に活用することにした。 画像のピクセルを座標とみなして、適当に生成したランダムの複素数を掛け算して画像に何かしらの処理を加える。 画像の枠からはみ出ると面白くないので無理やり枠内に収める処理を入れてある。トーラスっぽくなっている。

import itertools
import random
import time

from PIL import Image

image = Image.open("sample.jpg").convert('RGB')
size = image.size

for _ in range(5):
    image2 = Image.new('RGB', size)
    c = complex(random.gauss(0,1), random.gauss(0, 1))

    for x, y in itertools.product(range(size[0]), range(size[1])):
        pixel_data = image.getpixel((x,y))
        convert_complex = complex(x, y) * c
        convert_x = int(convert_complex.real) % size[0]
        convert_y = int(convert_complex.imag) % size[1]
        try:
            image2.putpixel((convert_x, convert_y), pixel_data)
        except IndexError:
            pass
    else:
        image2.save(str(int(time.time())) + ".png")

複素数と大学入試

実際に興味があったのは複素積分である。 まともにやると手も足も出ない積分が複素積分だと鮮やかに解けるのを学部2年で知った。 それを再体験したいのであるが、どうやってプログラミングで行うのか見当もつかない（線積分となるので）。 そこで、適当な大学入試問題を解くことにした。

2002年一橋大学の入学試験の数学の問2である。 問題文をTeX形式で入力するのは骨が折れるので適当なサイトを参照されたい。

import sympy

r = sympy.Symbol("r", real=True, positive=True)
t = sympy.Symbol("t", real=True)
a = sympy.Symbol("a")
f = sympy.Symbol("f", real=True)

a = r * (sympy.cos(t) + sympy.I * sympy.sin(t))
f = sympy.im(a + (1/a))
sympy.simplify(f)

最後に解くべき不等式を導くところまでできた。 SymPyで不等式を扱う方法がイマイチ理解できず最後は手計算で解いた。 自分の受けた課程では複素数平面は範囲外であったが、習っていれば高校二年生でも解けそうな雰囲気はある。

結論

現状、複素数を使わない場合はcmathはmathをC言語で実装した高速版ではなく複素数向けであることを理解する、複素数を本格的に使いたい場合は科学技術計算ライブラリを検討する、が現実的な解である。まだまだ私は恵子のように固定概念に囚われているようである。

PythonicではなくPythonesqueなクラスを書く その2

It’s…

youtu.be

最近、いかに役に立たないプログラムを書くかということを考えている。 世の中では企業が自社の開発の成果の一部をオープンソースとして切り出してたくさんの人に使われている。 そうすると自分も何か公開するには何か役に立つものを公開しなくてはいけないのではないか、という強迫概念に囚われてしまう。 業務では業務のプログラムを書いているので、ここでは何も役に立たないプログラミングをしてみたい。

今回は、Pythonesqueなクラスとして毎回変化する文字列のクラスを実装する。 文字列が評価されるたびに内部で書き換わるという役に立たない代物である。 現在のPythonではstrを直接継承してクラスを実装することが可能であるが、collections.UserStringを継承するほうがイタズラがしやすいと思われる。

import collections
import random
import string


class StrangeString(collections.UserString):


    def __str__(self):
        self.data = "".join([random.choice(string.ascii_letters) for _ in range(random.randint(1, 10))])
        return str(self.data)


if __name__ == "__main__":
    s = StrangeString("abc")
    for _ in range(10):
        print(s, s.upper(), s.lower())

実行結果は以下のとおりである。 実装からわかるように結果はランダムに変化する。

GCPp UAaiHZoGk UaTvs
ZoBIW olrixlE FUfCS
GTR S Itsq
CImAlqaHY ZBNY JoYBhYe
wdNtMJIAa EvQgkA oSfSRkj
DfwsTVxuaT cRZRtmOl BRKZjSav
GZW jcXG bpw
GpeN ecM EbwyaDsqK
naU rpRCnwhe vP
ZCD tK vqwJrSHkg

なるほど無秩序な結果となったが、無秩序な中にも秩序を求めたい。 s.upper()としたら結果が毎回ランダムになろうとも大文字になって欲しい。s.title()はタイトルであるべきだ。