感染伝播モデル

　待ち時間がたくさんできたので、数学やらネットワーク解析やらの基礎学習が捗る。が、数学疲れるのでまったり感染伝播モデルを実装してみた。

　2次元正方格子の各交点に確率qで木が生えている。中央で火を発生させた時、前後左右の4方向に木があれば燃え移る。燃え移った先の木の前後左右に木があればまた燃え移る。周りに木がなくなれば燃え移りは止まる。
　参考書に忠実に従っているわけではないのだが、とりあえず上記設定で、燃えた木の割合の平均θ(q) が急上昇するような相転移がみられるかどうかを調べる。

------------------------------------------------------------
#-*- coding:utf-8 -*-
#argv[1]: edge
#argv[2]: number of each simulation
#argv[3]: step size of q

import sys
import numpy as np
import matplotlib.pyplot as plt

def zeropad(data):
    x, y=data.shape
    data=np.c_[np.zeros(x), data, np.zeros(x)]
    data=np.r_[np.zeros(y+2).reshape(1, -1), data, np.zeros(y+2).reshape(1, -1)]
    return data

def spread(data, x, y):
    if data[x-1, y]==1:
        data[x-1, y]=2
        data=spread(data, x-1, y)

    if data[x+1, y]==1:
        data[x+1, y]=2
        data=spread(data, x+1, y)

    if data[x, y-1]==1:
        data[x, y-1]=2
        data=spread(data, x, y-1)

    if data[x, y+1]==1:
        data[x, y+1]=2
        data=spread(data, x, y+1)

    return data
               
N=np.int(sys.argv[1])
count=np.int(sys.argv[2])
q_list=np.arange(0, 1, np.float(sys.argv[3]))

if (N%2)!=1:
    print('N must be an odd integer.')
    quit()

theta=[]
for q in q_list:
    results=[]
    for i in range(count):
        data=np.random.choice(2, N**2, p=[1-q, q]).reshape(N, N)
        data=zeropad(data)
        center=np.array(data.shape)/2
        data[center[0], center[1]]=2
        out=spread(data, N/2, N/2)
        out=(out==2)
        results.append(np.mean(out[1:N+1, 1:N+1]))

    theta.append(np.mean(results))

plt.plot(q_list, theta)
plt.show()
-------------------------------------------------------------
Input:  edge:37, count:50, step size: 0.05
Output:

 それっぽくなった。ちなみにedge数37だと1/2の確率で再帰潜り過ぎでボイコットされる。39は不可であった。

お次はこれを立法格子へ応用。やることは一緒。
------------------------------------------------------------
import sys
import numpy as np
import matplotlib.pyplot as plt

def zeropad3d(data):
    x, y, z=np.array(data.shape)
    temp=np.zeros((x+2, y+2, z+2))
    temp[1:x+1, 1:y+1, 1:z+1]=data
    return temp

def spread3d(data, x, y, z):
    if data[x-1, y, z]==1:
        data[x-1, y, z]=2
        data=spread3d(data, x-1, y, z)

    if data[x+1, y, z]==1:
        data[x+1, y, z]=2
        data=spread3d(data, x+1, y, z)

    if data[x, y-1, z]==1:
        data[x, y-1, z]=2
        data=spread3d(data, x, y-1, z)

    if data[x, y+1, z]==1:
        data[x, y+1, z]=2
        data=spread3d(data, x, y+1, z)

    if data[x, y, z-1]==1:
        data[x, y, z-1]=2
        data=spread3d(data, x, y, z-1)

    if data[x, y, z+1]==1:
        data[x, y, z+1]=2
        data=spread3d(data, x, y, z+1)

    return data
              
N=np.int(sys.argv[1])
count=np.int(sys.argv[2])
q_list=np.arange(0, 1, np.float(sys.argv[3]))

if (N%2)!=1:
    print('N must be an odd integer.')
    quit()

theta=[]
for q in q_list:
    results=[]
    for i in range(count):
        data=np.random.choice(2, N**3, p=[1-q, q]).reshape(N, N, N)
        data=zeropad3d(data)
        center=np.array(data.shape)/2
        data[center[0], center[1], center[2]]=2
        out=spread3d(data, center[0], center[1], center[2])
        out=(out==2)
        results.append(np.mean(out[1:N+1, 1:N+1, 1:N+1]))

    theta.append(np.mean(results))

plt.plot(q_list, theta)
plt.show()
------------------------------------------------------------
Input: edge:9, count:50, step size:0.05
Output:

もうedge数9の時点でドキドキが足りない。しかし3*3*3のルービックキューブを27個積んであると思うと思ったよりでかい。

で、特に意味はないが4次元格子でも作ってみた。現実世界の時間に当てはめた場合にはそらもうSFよ。
実装はほんとにちょっと書き足すだけである。
------------------------------------------------------------
import sys
import numpy as np
import matplotlib.pyplot as plt

def zeropad4d(data):
    w, x, y, z=np.array(data.shape)
    temp=np.zeros((w+2, x+2, y+2, z+2))
    temp[1:w+1, 1:x+1, 1:y+1, 1:z+1]=data
    return temp

def spread4d(data, w, x, y, z):
    if data[w-1, x, y, z]==1:
        data[w-1, x, y, z]=2
        data=spread4d(data, w-1, x, y, z)

    if data[w+1, x, y, z]==1:
        data[w+1, x, y, z]=2
        data=spread4d(data, w+1, x, y, z)

    if data[w, x-1, y, z]==1:
        data[w, x-1, y, z]=2
        data=spread4d(data, w, x-1, y, z)

    if data[w, x+1, y, z]==1:
        data[w, x+1, y, z]=2
        data=spread4d(data, w, x+1, y, z)

    if data[w, x, y-1, z]==1:
        data[w, x, y-1, z]=2
        data=spread4d(data, w, x, y-1, z)

    if data[w, x, y+1, z]==1:
        data[w, x, y+1, z]=2
        data=spread4d(data, w, x, y+1, z)

    if data[w, x, y, z-1]==1:
        data[w, x, y, z-1]=2
        data=spread4d(data, w, x, y, z-1)

    if data[w, x, y, z+1]==1:
        data[w, x, y, z+1]=2
        data=spread4d(data, w, x, y, z+1)

    return data
              
N=np.int(sys.argv[1])
count=np.int(sys.argv[2])
q_list=np.arange(0, 1, np.float(sys.argv[3]))

if (N%2)!=1:
    print('N must be an odd integer.')
    quit()

theta=[]
for q in q_list:
    results=[]
    for i in range(count):
        data=np.random.choice(2, N**4, p=[1-q, q]).reshape(N, N, N, N)
        data=zeropad4d(data)
        center=np.array(data.shape)/2
        data[center[0], center[1], center[2], center[3]]=2
        out=spread4d(data, center[0], center[1], center[2], center[3])
        out=(out==2)
        results.append(np.mean(out[1:N+1, 1:N+1, 1:N+1, 1:N+1]))

    theta.append(np.mean(results))

plt.plot(q_list, theta)
plt.show()
------------------------------------------------------------
Input: edge:5,  count:50, step size:0.05
Output:

悲しみのedge数5。
5*5*5のルービックキューブ5個並んでるの想像するとそのしょーもなさが身にしみる。