Python で実現した私の最初のプロジェクトの 1 つは、スティック パーコレーションのモンテカルロ シミュレーションを行います。コードは継続的に成長しました。最初の部分は、スティック浸透の視覚化でした。width*length の領域に、特定の長さの定義された密度 (スティック/領域) の直線スティックが、ランダムな開始座標と方向でプロットされます。私は gnuplot をよく使うので、生成された (x, y) の開始座標と終了座標をテキスト ファイルに書き込んで、後で gnuplot にしました。
次に、scipy.ndimage.measurements を使用して画像データを分析する良い方法を見つけました。イメージは、グレースケールで ndimage.imread を使用して読み取られます。結果のnumpy配列は、異なるスティック間の接続にのみ関心があるため、ブール値にさらに縮小されます。結果のクラスターは、ndimage.measurements で分析されます。これにより、一方から他方につながる経路があるかどうかを調べることができます。最小化された例はこちらです。
import random
import math
from scipy.ndimage import measurements
from scipy.ndimage import imread
import numpy as np
import matplotlib.pyplot as plt
#dimensions of plot
width = 10
length = 8
stick_length = 1
fig = plt.figure(frameon=False)
ax = fig.add_axes([0, 0, 1, 1])
fig.set_figwidth(width)
fig.set_figheight(length)
ax.axis('off')
file = open("coordinates.txt", "w")
for i in range (300):
# randomly create (x,y) start coordinates in channel and direction
xstart = width * random.random() # xstart = 18
ystart = length * random.random() # ystart = 2
# randomly generate direction of stick from start coordinates and convert from GRAD in RAD
dirgrad = 360 * random.random()
dirrad = math.radians(dirgrad)
# calculate (x,y) end coordinates
xend = xstart + (math.cos(dirrad) * stick_length)
yend = ystart + (math.sin(dirrad) * stick_length)
# write start and end coordinates into text file for gnuplot plotting
file.write(str(i) + ":\t" + str(xstart) + "\t" + str(ystart) + "\t" + str(dirgrad) + ":\t" + str(xend) + "\t" + str(yend) + "\n")
file.write(str(i) + ":\t" + str(xend) + "\t" + str(yend) + "\n\n")
# or plot directly with matplotlib
ax.plot([xstart,xend],[ystart,yend],"black", lw=1)
fig.savefig("testimage.png", dpi=100)
# now read just saved image and do analysis with scipy.ndimage
fig1, ax1 = plt.subplots(1,1)
img_input = imread("testimage.png", flatten = True) # read image to np.ndarray in grey scales
img_bw = img_input < 255 # convert grey scales to b/w (boolean)
labeled_array, num_clusters = measurements.label(img_bw) #labeled_array: labeled clusters in array, num_clusters: number of clusters
area = measurements.sum(img_bw, labeled_array, index=np.arange(labeled_array.max() + 1)) # area of each cluster
areaImg = area[labeled_array] # label each cluster with labelnumber=area
cax = ax1.imshow(areaImg, origin='upper', interpolation='nearest', cmap = 'rainbow')
cbar = fig1.colorbar(cax)
fig1.savefig("testimage_analyzed.png")
これは基本的にうまくいきますが、より多くの異なるスティック密度に対して 1000 回の反復を行う細かいモンテカルロ シミュレーションは、8 時間以上実行されることになります。これは、作成された画像と配列が非常に大きく、何千ものスティックがより高い密度でプロットされているという事実によるものです。その理由は、ピクセル化によるエラーを最小限に抑えながら、さまざまなジオメトリ (たとえば、500 ～ 20000 ピクセルの長さ) をシミュレートしたいからです。
画像データを使用せずにベクトル問題として扱うのが最善の方法だと思いますが、アルゴリズムを開始する方法さえわかりません。また、多くの接続により、データ配列が大きくなる可能性もあります。
上記の方法にとどまると、ファイルにデータを書き込んで再読み込みすることはあまり効果的でないことは明らかです。したがって、これを高速化する方法を探しています。最初のステップとして、matplotlib を使用して画像を作成しましたが、少なくとも各スティックを個別のプロット呼び出しでプロットする場合、スティックの数が多いと最大 10 倍遅くなります。配列内にスティック座標のリストを作成し、1 回のプロット呼び出しで完全なリストをプロットすると、速度が向上する可能性がありますが、画像の書き込みと読み取りのボトルネックが残ります。
スティックの白黒画像を表すブール型のnumpy配列を直接生成する効果的な方法を教えてください。たぶん、座標のリストをプロットし、何らかの方法で図を配列に変換しますか? また、線が PIL 画像にプロットされる興味深い議論も見つけました。これは matplotlib よりも高速でしょうか?
