Pythonでパーコレーション理論におけるクラスター形成のアルゴリズムを組む

投稿日：2021/09/18

はじめに

小田垣先生著作の『つながりの物理学ーパーコレーション理論と複雑ネットワーク』の付録ページのクラスター形成のアルゴリズムを読んで実際に作ってみました。Numpyで基本的に組むため、本とは少し定義が異なります。Pythonで組んだコードを探してみた所、Numpyなどの基本ライブラリで構成された文献が見当たらなかったので参考になれば幸いです。

クラスターとは

$N\times N$の正方格子を考えます。この正方格子上の格子点を$(i,j)\quad (i,j=0,\dots,N-1)$で表し、占有確率$p$のとき、$p\times N^2$個の格子点がランダムに占有されているとします。このとき、ある占有されている格子点$(k,l)$があり、$(k\pm 1,l)$、$(k,l\pm 1)$のいずれかが占有されているとき、その占有されている格子と$(k,l)$は同じクラスター (cluster)に所属するといいます。物理的に興味がある量はこのクラスターの平均値である***平均クラスターサイズ (average cluster size)***です。これはこの本に紹介されているように、非磁性原子中に磁性原子を添加した磁性体の磁化率とこの平均クラスターサイズが関係づけられる事があるようです。他にも、有機半導体中のドナー、アクセプターの添加量と平均クラスターサイズを結び付けた報告を近年よく見かけます。この量を数値的に決めるためにはある格子点がクラスターに所属していると判定するプログラムが必要になります。

クラスター形成のアルゴリズム

for文で行、列の順番に格子点にアクセスしていく事を考えます。$(0,j)\quad (j=0,\dots, N-1)$を左から順に調べ、最初の要素にクラスターの番号$1$をつけます。次の要素が占有されていれば同じ番号を付け、占有されてなければ次に進みます。次の格子点が占有されていればクラスターの番号$2$を付け、次に進みます。これを繰り返し、0行目について調べ終わったら1行目に移動します。$(i,0)\quad (i\neq 0)$は$(i-1,0)$を見て、占有されていれば同じクラスター番号を振り、占有されていなければ新たにクラスター番号をふります。$(i,j)\quad (i\neq 0,j\neq 0)$では$(i,j)$が占有さているときは、$(i-1,j)$と$(i,j-1)$の格子点が占有状況によって対処法が変わります。

• 両方占有されているときは$(i,j)$は両点と同じクラスターに所属している事がいえるので$(i-1,j)$と$(i,j-1)$のクラスター番号のうちどちらかを$(i,j)$にふり、$(i-1,j)$と$(i,j-1)$の番号が同じであることを保持する
• $(i-1,j)$か$(i,j-1)$のみが占有されているときは占有されている方のクラスター番号を$(i,j)$にふる
• どちらも占有されていないときは新たに$(i,j)$にクラスター番号をふる

プログラムを組む

まずは正方格子上を占有率$p$でランダムに占有させ、占有されている格子点を$1$、占有されていない格子点を$0$とおき、$0$と$1$のみを要素として持つ格子点の情報を持つ行列を出力するsqGrid関数を作成します。この格子からクラスターを判別するjdgCluster関数に入れ、drawLattice関数でそのクラスターを描画します。占有確率$p$で占有されるという事は$pN^2$個必ず占有される事を意味しないのですが、今回は$pN^2$個の格子点が必ず占有されるように作成しています。また、このために$pN^2$をintで無理やり変換しているので必要であれば書き直して下さい。

import numpy as np
import matplotlib.pyplot as plt

def sqGrid(N,p):
    """
    サイズgridX*gridYの正方格子上にp*gridX*gridY個の粒子を生成
    ------------------------------------------------------------
    gridX：array
        行数
    gridY：array
        列数
    p：float
        [0,1]の実数。占有確率。
    """
    n = int(p*N**2) # 粒子数
    occPoint = np.zeros((N,N)) 
    # p==1は全て占有するため、以降の処理をスキップ
    if p == 1:
        return occPoint + 1
    if p == 0:
        return occPoint
    X = np.random.randint(0,N,n) # 占有される格子のX軸情報
    Y = np.random.randint(0,N,n) # 占有される格子のY軸情報
    for i,j in zip(X,Y):
        occPoint[i,j] = 1
    
    # 重複した座標用の処理
    while int(occPoint.sum()) != n:
        c = n - int(occPoint.sum())
        x = np.random.randint(0,N,c)
        y = np.random.randint(0,N,c)
        for i,j in zip(x,y):
            occPoint[i,j] = 1
    return occPoint

def jdgCluster(lattice):
    cluster = np.zeros_like(lattice)
    clusterNum = 0
    clusterMerge = []
    xnum,ynum = cluster.shape
    for i in range(xnum):
        for j in range(ynum):
            # 0行目の処理
            if i==0:
                # (0,0) の処理
                if j==0:
                    if lattice[i,j]==1:
                        clusterNum += 1
                        cluster[i,j] = clusterNum
                # (0,j) (j≠0)の処理
                else:
                    if lattice[i,j]==1:
                        if lattice[i,j-1]==1:
                            cluster[i,j] = cluster[i,j-1]
                        else:
                            clusterNum += 1
                            cluster[i,j] = clusterNum
            # 1行目以降の処理
            else:
                # (i,0) (i≠0)の処理
                if j==0:
                    if lattice[i,j]==1:
                        if lattice[i-1,j]==1:
                            cluster[i,j] = cluster[i-1,j]
                        else:
                            clusterNum += 1
                            cluster[i,j] = clusterNum
                # (i,j) (i,j≠0)の処理
                else:
                    if lattice[i,j]==1:
                        # (i-1,j)と(i,j-1)がどちらも占有されている
                        if lattice[i-1,j]==1 and lattice[i,j-1]==1:
                            cluster[i,j] = cluster[i,j-1]
                            clusterMerge.append([cluster[i-1,j],cluster[i,j]])
                        # (i-1,j)が占有されている
                        elif lattice[i-1,j]==1:
                            cluster[i,j] = cluster[i-1,j]
                        # (i,j-1)が占有されている
                        elif lattice[i,j-1]==1:
                            cluster[i,j] = cluster[i,j-1]
                        # どちらも占有さてない
                        else:
                            clusterNum += 1
                            cluster[i,j] = clusterNum

    # clusterMergeに保存された同一とみなせるクラスターの番号の振り直し
    if len(clusterMerge)!=0:
        for i in clusterMerge:
            cluster = np.where(cluster==i[0],i[1],cluster)
    # 数字を振りなおす
    idxUni = np.unique(cluster) # 現在のクラスターのインデックス
    idxRe = np.arange(0,len(idxUni))
    for uni,re in zip(idxUni,idxRe):
        cluster = np.where(cluster==uni,re,cluster)
    return cluster

def drawLattice(lattice,N,p,Name='test'):
    img = lattice.copy()
    img = img*(255)/np.amax(img)
    plt.title(f"N={N},p={p}")
    plt.axis("off")
    plt.imshow(img)
    plt.savefig(f'{Name}.png',transparent=True,dpi=300)

# main
N = 100
p = 0.5
sg = sqGrid(N,p)
clu = jdgCluster(sg)
drawLattice(clu,N,p)