EnsekiTT Blog

EnsekiTTが書くブログです。

PythonとMatplotlibを使ってシダを描画するのとプログラムの呼び出し回数の話

つまりなにしたの?

先日のライフゲームと同じ要領で「プログラムでシダを描画する」をやってみた2018。
ついでに確率的に実行されているシダの関数呼び出しがどんな分布を描くのか調査した。
f:id:ensekitt:20180526042808j:plain

シダってなに?

2013年末〜2014年頃にプログラミングする界隈でちょっと流行ってQiitaや個人Blogにいろいろな環境でやってみた記事が上がっていた。
2018年、「そういえば、そんな課題あったな?」と思いだしたので描画ネタとして今回ピックアップした。
目新しいことは特に無いのでQiitaを見に行くといろいろな実装が見られて面白い。
あと元ネタになったYouTubeの動画は消えてしまっているみたい。

qiita.com

ただ、いろいろな実装があるものの、一度その実装を見ちゃうと引っ張られるのでほぼ写経になった。

d.hatena.ne.jp

シダを描く関数

def shida(k, x, y):
    if 0 < k:
        for point in shida(k-1, W1x(x, y), W1y(x, y)):
            yield point
        if random.random() < 0.3:
            for point in shida(k-1, W2x(x, y), W2y(x, y)):
                yield point
        if random.random() < 0.3:
            for point in shida(k-1, W3x(x, y), W3y(x, y)):
                yield point
        if random.random() < 0.3:
            for point in shida(k-1, W4x(x, y), W4y(x, y)):
                yield point
    else:
        yield x * WIDTH*0.9 + WIDTH * 0.5, HEIGHT - y * HEIGHT * 0.9

乱数を使ってW1〜W4が確率的に呼び出されるような形になっている。
yieldを使ってforループを回すことでshidaを呼び出したら都度pointを返す形になっている。
kが0未満になるまで再起していって、0未満になると引数で渡されたx、yに描画用の係数をかけて値を返している。

0.3未満の乱数で実行するので何度もやっていたら、shidaは1回の呼び出しに対して概ね1.9回ずつ実行される。
仮にkの初期値が20だったら、1.9^20→375899回呼び出される気がする。

呼び出し回数があっているか確認した

呼び出すと都度ポイントを返してくれるので、ループ回数を数えれば良いはず。
kの初期値を20にして、シダを1000回生成してそのヒストグラムを出してみた。
f:id:ensekitt:20180526040923p:plain
なんとなく近い感じ。
平均値は367894、中央値は337747だった。

実行例

f:id:ensekitt:20180526041715g:plain

コード

呼び出し回数をカウントするコードはコメントアウトしてあるのでそのまま動くコード

import random
import copy
import matplotlib.pyplot as plt

N = 20
xm = 0
ym = 0.5
h = 0.6

WIDTH = 500
HEIGHT = 500

W1x = lambda x, y: 0.836 * x + 0.044 * y
W1y = lambda x, y: -0.044 * x + 0.836 * y + 0.169
W2x = lambda x, y: -0.141 * x + 0.302 * y
W2y = lambda x, y: 0.302 * x + 0.141 * y + 0.127
W3x = lambda x, y: 0.141 * x - 0.302 * y
W3y = lambda x, y: 0.302 * x + 0.141 * y + 0.169
W4x = lambda x, y: 0
W4y = lambda x, y: 0.175337 * y

def shida(k, x, y):
    if 0 < k:
        for point in shida(k-1, W1x(x, y), W1y(x, y)):
            yield point
        if random.random() < 0.3:
            for point in shida(k-1, W2x(x, y), W2y(x, y)):
                yield point
        if random.random() < 0.3:
            for point in shida(k-1, W3x(x, y), W3y(x, y)):
                yield point
        if random.random() < 0.3:
            for point in shida(k-1, W4x(x, y), W4y(x, y)):
                yield point
    else:
        yield x * WIDTH*0.9 + WIDTH * 0.5, HEIGHT - y * HEIGHT * 0.9

epoch_list = []
for i in range(1):
    env_grid = [[0 for i in range(WIDTH)] for j in range(HEIGHT)]
    show_grid = copy.copy(env_grid)

    fig, ax = plt.subplots(1,1, figsize=(8,8))
    imgs = ax.imshow(env_grid, cmap="GnBu", vmax=1, vmin=0)

    for epoch, p in enumerate(shida(N, 0, 0)):
        # print(epoch, int(p[0]), int(p[1]))
        env_grid[int(p[1])][int(p[0])] = 1
        if epoch % 1000 == 0:
            show_grid = copy.copy(env_grid)
            # 表示する
            imgs.set_data(env_grid)
            plt.pause(.001)
    print(i, epoch)
    epoch_list.append(epoch)
    #plt.clear()
#with open('epoch.txt','w') as f:
#    f.write(str(epoch_list))
#plt.hist(epoch_list)
#plt.show()

やってみてどうだった?

実行回数は結構ばらつくんだなーってヒストグラムを見たときには思ったけど、ツリー状に再起していくんだから序盤の乱数が0.3より小さい数字だらけだと簡単に増えるのは直感的にも普通だった。

クリエイティブ・コモンズ・ライセンス
この 作品 は クリエイティブ・コモンズ 表示 4.0 国際 ライセンスの下に提供されています。