今回はウラムの螺旋についてです。 
ドラマ「危険なビーナス」を見ていて
気になってしまったので
これを作るプログラムを
作成してみようと思います。

解説動画はこちら



まずウラムの螺旋とは
数学者のウラムさんが発見と言うか作った
素数をとある法則に従って並べた図形です。


1から渦巻状に数値を配置する際に
素数だけ残すと、対角線状に素数が残ります。

まずは1から数値を並べるところを
考えていきましょう。

def uram(N=7):
    # 箱を用意
    data = [ [0 for n1 in range(N)] for n2 in range(N)]

    # 入れる値を生成
    nums,tmp,c = [],[],1
    for num in range(1,N*N+1):
        tmp.append(num)
        if num==c*c:
            nums.append(tmp)
            tmp,c=[],c+1
            #c+=1

    # 起点位置決め
    c_y,c_x = N//2,N//2 if N%2==0 else N//2+1

    # 値の代入
    for i,tmp in enumerate(nums,start=1):
        for j,n in enumerate(tmp):
            sub = 1 if i%2==0 else -1
            if j==0:
                c_y,c_x=c_y,c_x + sub*1
            elif j>len(tmp)//2:
                c_y,c_x=c_y,c_x - sub*1
            else:
                c_y,c_x=c_y - sub*1,c_x
            data[c_y][c_x]=n
    return data

data = uram(5)
for row in data:
    for t in row:
        print('{0:}'.format(t),end='\t')
    print()

最初にN*Nの箱を作ります。

次に入れる数値の生成ですね。
こんな感じでN*Nの数値のリストを作ります。

N = 5
nums,tmp,c = [],[],1
for num in range(1,N*N+1):
    tmp.append(num)
    if num==c*c:
        nums.append(tmp)
        tmp,c=[],c+1
print(nums)

[[1],
[2, 3, 4],
[5, 6, 7, 8, 9],
[10, 11, 12, 13, 14, 15, 16],
[17, 18, 19, 20, 21, 22, 23, 24, 25]]

ここでは螺旋状に数値を入れる際に
折り返しがあるので
入れる数値のリストの半分のところで
折り返すようにしています。

Nが偶数の時はその前の数値から
右に一つずれ、そこから上
真ん中のところで左に折り返し

Nが奇数の時はその前の数値から
左に一つずれ、そこから下
真ん中のところで右に折り返す
ように数値を入れていくと
螺旋状に数値が並んだ
リストが生成できます。

次に素数だけを
残すようにしてみましょう。

sympyのisprimeと言うメソッドで
素数判定ができます。

これを使って素数の部分だけを
残すように変えてみます。

from sympy import isprime

def uram2(N=7):
    # 箱を用意
    data = [ [0 for n1 in range(N)] for n2 in range(N)]

    # 入れる値を生成
    nums,tmp,c = [],[],1
    for num in range(1,N*N+1):
        tmp.append(num)
        if num==c*c:
            nums.append(tmp)
            tmp=[]
            c+=1

    # 起点位置決め
    c_y,c_x = N//2,N//2 if N%2==0 else N//2+1

    # 値の代入
    for i,tmp in enumerate(nums,start=1):
        for j,n in enumerate(tmp):
            text = ' | '
            sub = 1 if i%2==0 else -1
            if j==0:
                c_y,c_x=c_y,c_x + sub*1
            elif  j>len(tmp)//2:
                text = ' - '
                c_y,c_x=c_y,c_x - sub*1
            else:
                c_y,c_x=c_y - sub*1,c_x
            # 素数判定
            if isprime(n):
                data[c_y][c_x]=n
            else:
                data[c_y][c_x]=text
    return data

data = uram2(5)
for row in data:
    for t in row:
        print('{0:}'.format(t),end='\t')
    print()

はいこれで素数だけ残りましたね。

次は塗りつぶしですが
数値を文字列に変えるだけでできます。


これだと使いづらいので
画像にするようにしてみましょう。

matplotlibではリストを投入すると
そのまま画像化できます。

素数を0,他は255を入れた
リストを作ってplotするだけです。

せっかくなので1000x1000
100万までの素数で
ウラムの螺旋を作ってみましょう。

from sympy import isprime
import matplotlib.pyplot as plt
%matplotlib inline

def uram4(N=7):
    # 箱を用意
    data = [ [255 for n1 in range(N)] for n2 in range(N)]

    # 入れる値を生成
    nums,tmp,c = [],[],1
    for num in range(1,N*N+1):
        tmp.append(num)
        if num==c*c:
            nums.append(tmp)
            tmp=[]
            c+=1

    # 起点位置決め
    c_y,c_x = N//2,N//2 if N%2==0 else N//2+1

    # 値の代入
    for i,tmp in enumerate(nums,start=1):
        for j,n in enumerate(tmp):
            sub = 1 if i%2==0 else -1
            if j==0:
                c_y,c_x=c_y,c_x + sub*1
            elif  j>len(tmp)//2:
                c_y,c_x=c_y,c_x - sub*1
            else:
                c_y,c_x=c_y - sub*1,c_x
            if isprime(n):
                data[c_y][c_x]=0

    return data

data = uram4(1000)
plt.figure(figsize=(10,10))
plt.imshow(data)
plt.axis("off")
plt.show()

メモリとCPUが強力であれば
もっと大きな螺旋も描けると思います。

これ以上はなかなかに痺れますね。

素数の法則は見つかっていないと
思いますが、なんかロマンがある図形ですね。

お暇な方は試してみてください。
それでは。

タグ ：

#プログラミング

#Python

#可視化

#ウラムの螺旋

#matplotlib

最近見ているドラマでフラクタル図形が
取り上げられていたので
気になっちゃいました。

今回は色々なフラクタル図形を作図します。

解説動画はこちら




さてフラクタル図形と言うものは
どう言う図形でしょうか？

フランスの数学者ブノワ・マンデルブロが
導入した幾何学の概念で

タグ ：

#プログラミング

#Python

#可視化

#フラクタル図形

#危険なビーナス

#matplotlib

前回作った関数を進化させてみました。


解説動画はこちら




前回はこちら
動くおっぱい関数

さて
Python言語ではmatplotlibというライブラリで
作図、可視化を行うことができます。

ということなので
曲線を描くことなどが出来てしまうわけです。

前回は黒い曲線で可視化しましたが
それだと味気ないので
色を塗ってみました。

ただし、これは静止画ではありません。
動くんです！！！

なので毎回色ぬりする範囲を
計算しないといけません。

うまく塗る範囲を考えてみました。

出来上がった関数はこちら

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from ipywidgets import interact, FloatSlider, IntSlider
import warnings
warnings.simplefilter('ignore')
%matplotlib inline

def oppai(y,t):
    x_1 = (1.5 * np.exp((0.12*np.sin(t)-0.5) * (y + 0.16 *np.sin(t)) ** 2)) / (1 + np.exp(-20 * (5 * y + np.sin(t))))
    x_2 = ((1.5 + 0.8 * (y + 0.2*np.sin(t)) ** 3) * (1 + np.exp(20 * (5 * y +np.sin(t)))) ** -1)
    x_3 = (1+np.exp(-(100*(y+1)+16*np.sin(t))))
    x_4 = (0.2 * (np.exp(-(y + 1) ** 2) + 1)) / (1 + np.exp(100 * (y + 1) + 16*np.sin(t)))
    x_5 = (0.1 / np.exp(2 * (10 * y + 1.2*(2+np.sin(t))*np.sin(t)) ** 4))
    x = x_1 + (x_2 / x_3) + x_4 + x_5
    return x
 
t = FloatSlider(min=0.1, max=5.0, step=0.1, value=0)
y = np.arange(-3, 3.01, 0.01)
 
@interact(t=t)
def plot_oppai(t):
    x = oppai(y,t)
    plt.figure(figsize=(10,9))
    plt.axes().set_aspect('equal', 'datalim')
    plt.grid()
    b_chiku = (0.1 / np.exp(2 * (10 * y + 1.2*(2+np.sin(t))*np.sin(t)) ** 4))
    b_index = [i for i ,n in enumerate(b_chiku>3.08361524e-003) if n]
    x_2,y_2 = x[b_index],y[b_index]
    plt.axes().set_aspect('equal', 'datalim')
    plt.plot(x, y, '#F5D1B7')
    plt.fill_between(x, y, facecolor='#F5D1B7', alpha=1)
    plt.plot(x_2, y_2, '#F8ABA6')
    plt.fill_between(x_2, y_2, facecolor='#F8ABA6', alpha=1)
    plt.show()

いやー素晴らしい曲線ですよね。

x,y座標に描く曲線
これは前回と一緒です。

色ぬりの指定は

#F5D1B7

のように16進数のwebカラー表記で
書いて指定します。

デフォルトは黒いようなので
肌色を指定しています。

webで検索すれば好みの色に該当する
数値がわかるはずです。

次のポイントは
色ぬりの範囲指定です。

plt.fill_between

これで色ぬりができるんですが
範囲はxの最小値からxの最大値のようですね。

まず全体を肌色で塗り込みます。

後はB地区の計算です。

B地区の範囲を計算している部分は
def oppaiのx_5を計算している部分です。

これが曲線全体における
一番の盛り上がりの部分を
計算しているところです。

このdef oppaiではyの値を元に
x座標を計算して出力していますので
盛り上がり部分の数値がわかれば
塗る範囲も決まります。

この数値を計算すると盛り上がり部分が
ある一定の数以上の値になる部分でした。

ここを
x,yそれぞれをインデックスで取り
その部分だけピンク色を指定して
塗っています。

ようやく完成です。

tの値を変えると
プルンプルン動きます。


いやー今年一番の感動です。
暇でよかった!!!!

matplotlibであれば
あなたの夢も可視化できるかもしれませんね

今回はこれまでです
それでは。

タグ ：

#Python

#プログラミング

#可視化

#matplotlib