階乗のループがある複雑ネットワーク解析プログラムを高速化したい

Question

実現したいこと

複雑ネットワークでウィルスが感染する様子を調べたいのですがデータの数が多く処理に時間が掛かってしまうため高速したいです。
###プログラムの概要
1step毎に7000～10000の数のノードデータ(ノードの名称id,x座標、y座標)が記載されたものが1800stepまであり、1step毎に読み込み、感染ノードであれば1を、それ以外のとき0をnode_listの最後に追加します。その後、関数を呼び出しノード同士の座標をそれぞれ比較して、座標の距離を求めて範囲内であれば感染リストに格納します。感染リストは、次のstep時に、参照して、新たに感染リストに追加していきます。stepの最後に時間と感染リストの数を出力します。

発生している問題・エラーメッセージ

1step毎にノードの数の階乗回for文でループしているため、1stepあたり約40秒ほど時間がかかってしまってしまい、1800秒間では約20時間掛かっています。

該当のソースコード

import sys
import random
import math
import time
import re

#２点の差が縦横100mの距離の長さ
#sqrt(100^2+100^2)
D_100=math.sqrt((100*100)+(100*100))

Input = "node_data.txt"
Output="output.txt"
def file_write(output,infeced_list,t):
    with open(output, 'a') as f:
        f.write("%d %d\n" % (t,len(infeced_list)))
    f.close()

def SUMOTraceParser():
    #ファイルオープン
    f1 = open(FCDOutput, 'r')
    node_list = []　# ノードデータ保存リスト((車のid,x,y,レーン))
    infected_list=[]　#感染リスト(id)
    row_no = 0　#行数カウント
    while True:
        line = f1.readline()
        line = line.replace("\n", "") # 改行削除
        #timestep
        if "timestep" in line:
            t = re.sub(r"[^\d.]", "", line)
            t = int(float(t))
            row_no+=1
        #ノードデータ読み込み、node_listに追加
        elif "end" not in line and "timestep" not in line:
            row_no += 1
            # 分割
            l = line.split(',')
            # ダブルコーテション削除
            for i in range(len(l)):
                l[i] = l[i].replace("''", "")
                l[i]=float(l[i])
                l[i]=int(l[i]) 
            # 感染リストにidがあるとき
            if l[0] in infected_list:
                node_list.append((l[0],l[1],l[2], 1)) #1を追加
            else:
                node_list.append((l[0],l[1],l[2], 0)) #0を追加
        #1stepの読み込み終了時
        elif "end" in line:
            test_list=Infect(node_list,t,infected_list) #関数呼び出し
            infected_list=test_list
            file_write(Output,infected_list,t) #ファイルに感染したノードを書き込み
            node_list=[]
            row_no+=1 
        if t == 1800:
            # 強制終了
            sys.exit()
#ノード同士を比較して、感染リストに格納する関数
def Infect(node_list, t, infected_list):
    node_list.sort()　#sort
    node_list2=reversed(node_list) #revresed
    node_id_list=set() #idリスト
　　#node_listを階乗回ループする
    for i in range(len(node_list)):
        for j in range(len(node_list2)-i):
            A=None
            #x,y差を止める
            x_1=node_list[i][1] - node_list[j][1]
            y_1=node_list[i][2] - node_list[j][2]
            #2乗する
            x_2=pow(x_1,2)
            y_2=pow(y_1,2)
            D=math.sqrt(x_2+y_2)　#2点間の距離計算
            if D<D_100:
                if node_list[i][0] != node_list[j][0]:　#同一ノードではないとき
                    # どちらかが1（感染状態）のとき
                    if node_list[i][3] == 1 and node_list[j][3] == 0:
                        A=True
                    if node_list[i][3] == 0 and node_list[j][3] == 1:
                        A=True
                    if A==True:
                        A=None
                        # 確率1%の確率で
                        b=random.random()
                        if b < 0.1:　
                            infected_list.append(node_list[j][0])　#感染リストに追加

            node_id_list.add(node_list[i][0])

    # 感染するノードをランダムに選択
    if t==start_s:
        infected_list=random.sample(node_id_list,1)
    # リストの重複を削除
    infected_list = list(set(infected_list))
    return infected_list

if __name__ == '__main__':
    SUMOTraceParser()

Answer 1

scipy.spatial.KDTreeを使うと距離計算が楽になります。
参考：【Python】scipy KD treeの使い方
以下は人は入れ替わらない前提でのコード例です。

import numpy as np
from scipy import spatial

FIELD_SIZE = 25000
INFECT_DIST = 100
INFECT_RATE = 0.1
N_AGENT = 10000

infecteds_cur = np.zeros(N_AGENT, dtype=int) # 感染状態 1=感染

# 最初の感染者
np.random.seed(110)
init_idx = np.random.randint(N_AGENT)
infecteds_cur[init_idx] = 1

for step in range(1800):
    infecteds_next = infecteds_cur.copy()

    # とりあえず人の位置をランダムに更新
    positions = np.random.uniform(0,FIELD_SIZE, [N_AGENT,2])
    tree = spatial.KDTree(positions)

    # 感染者
    infected_idxs = np.where(infecteds_cur == 1)[0]
    for infected_i in infected_idxs:

        # 近辺の人を取得
        neighbours = np.array(tree.query_ball_point(positions[infected_i], INFECT_DIST))
        neighbours = neighbours[neighbours != infected_i] # 自分は除く

        for n_i in neighbours:
            # 未感染→感染
            if infecteds_cur[n_i] == 0 and np.random.uniform(0,1) < INFECT_RATE:
                infecteds_next[n_i] = 1

    infecteds_cur = infecteds_next # 次の状態へ
    infected_cnt = len(infecteds_cur[infecteds_cur == 1])

    print(step, infected_cnt)

    # 人は治らず入れ替わらないので打ち切り
    if infected_cnt >= N_AGENT:
        break

Answer 2

階乗ではなく2乗ですね。
色々と問題は見受けられますがそれはさておき、
ループ回数を減らす手段としては、以下が思いつきます。
・座標を100m×100mのマス目に区切り、まず全ノードをマスに割り振る
・2×2の4マスごとに、それらのマスにあるノードの間で処理をする
・1マスづつずらして繰り返す
座標の範囲が狭ければ単純に配列にすればよく、範囲が広ければ疎行列的な扱いが必要になるでしょう。
連想配列として扱い、隣のマスを求めるのが面倒なので逆にノードごとに隣を含む4マスに割り振っておくのが楽かなと思います。

Comments

[Dobutori]

座標範囲は縦25000m,横20000mのデータなので、この場合は疎行列の方がいいのでしょうか。

[ikadzuchi]

250*200で50kマスですか。マスあたり(ポインタとして)8バイトとしても高々400kBなので密でじゅうぶん可能だと思います。
ただ改めて考えてみると、密か疎かというのはC言語などでメモリ上のデータ配置を考えるような場合についての話で、Pythonのようなメモリアドレスを直接扱わない言語ですとそもそも考える必要は無かったかもしれません。