今回は、3月18日にまとめたシミュレーションの続編です。
実は、やり切った感あったのですが、よく考えてみると、以下二点を深堀する必要がありました。
①　SIRモデルで説明できない感染伝播
②　クラスターで発生する場合の伝播
つまり、①は前回の粒子密度依存性の遷移領域の感染数の推移を見ると、感染数のピークがいくつかに分かれており、これはSIRモデルでは説明できません。また、現実の感染は多くの場合、クラスターで発生しており、なんとなく単調に感染が伝播しているようには見えません。
ということで、今回は、以下を実施しました。
・臨界密度領域での感染伝播
・感受性保持者の分布がクラスター的に分布している場合の感染伝播
なお、コードは前回のものとほぼ同一ですが、感受性保持者の分布と初期値に速度を持たせるように拡張しています。コードは、おまけに載せておきます。
※今回は速度を持たせるシミュレーションは割愛します
###・臨界密度領域での感染伝播
前回のように感受性保持者がランダムに一様に分布しているとき、
①感染範囲；Rc
②感染確率；p
③感受性保持者の密度；d
④感受性保持者の移動速度；mr
を変化させると、感染伝播及び感染率を制御できる。
####代表的な感染伝播
右側の図が青；感受性保持者数、赤；感染者、緑；治癒者であり、このグラフはよく見るSIRモデルで馴染みな絵である。
パラメータ；①Rc=20x20, ②p=0.05, ③d=1000/200x200, ④mr=1.0

####臨界的な感染伝播
以下のように感染確率を下げて行くと、やがて感染伝播が起こらなくなる。その一歩手前の感染伝播は以下のパラメータで実現出来て下記のGifアニメのようになる。
#####パラメータ；①Rc=20x20, ②p=0.0075, ③d=1000/200x200, ④mr=1.0

この図から分かることは、右の図のとおり、この伝播はもはやSIRモデルでは表せない領域に入っており、確率的に伝播が起きていることが分かる。
この図では、立ち上がりこそ指数関数的に立ち上がっているが、その後はだらだらとテイルを引いて感染伝播と治癒が発生している。
#####パラメータ；①Rc=13x13, ②p=0.03, ③d=1000/200x200, ④mr=1.0
感染範囲を小さくしても臨界的な感染伝播は実現できて、以下のようになった。

このGifアニメでは、立ち上がりさえ、指数関数的な増加は無くなり、だらだらと増加してよれよれと減少し、また増加している。
このようにこの領域での感染伝播は感受性保持者の分布に依存しており、ほぼ確率的なプロセスであると言える。
###・クラスター的に分布している場合の感染伝播
実際の分布は通常は村や町ごとに人口密度の揺らぎが発生していたり、人の集会や集まり毎に分布に密度揺らぎがあると想定される。
ということで、感受性保持者の分布がクラスター的に分布している場合の感染伝播をやってみた。
####クラスターが密に存在する場合
今回は簡単のために、まず10x10のグリッド上に分布するクラスターで感染伝播を見た。
#####パラメータ；①Rc=20x20, ②p=0.5, ③d=100x10/200x200, ④mr=1.5
d=100x10/100x100の意味は以下の図のように200x200のエリアに100のクラスターをグリッド上の点に10個ずつ感受性保持者を置くことを意味する。

このパラメータだと、感染伝播は一様分布と変わらずSIRモデルで表現できそうな領域となっている。
#####パラメータ；①Rc=20x20, ②p=0.5, ③d=100x10/200x200, ④mr=1.25
感受性保持者のランダムウォークの速度が下がると、感染者との接触確率が下がり、以下のようにはみだし感受性保持者による感染に依存して、クラスター間の伝播が発生して、以下のように不規則な伝播が発生する。

#####パラメータ；①Rc=20x20, ②p=0.5, ③d=100x10/200x200, ④mr=1.09
そして速度が1.09辺りがクリティカルで伝播するかどうかの境界になる。

そして、これより速度が小さくなると、伝播は消える。つまり感染しなくなる。感染する前に最初の感染者や少なくとも隣のクラスターの一個だけで感染伝播が止まってしまう。
####クラスターが疎に存在する場合
もっとクラスターが大きく、それぞれが離れている場合はどうなるだろうか。すなわち、イメージは大きな村が点在する場合である。
#####パラメータ；①Rc=20x20, ②p=0.5, ③d=4x250/200x200, ④mr=10
d=4x250/100x100の意味は、200x200に4個の250人ずつの感受性保持者を含むクラスターが点在していることを意味する。

この場合は、クラスターの分布の所為で少し崩れるが一応単一ピークとなる。
#####パラメータ；①Rc=20x20, ②p=0.5, ③d=4x250/200x200, ④mr=7.5
このクラスター分布の場合はこのパラメータ辺りがクラスター間の伝播という意味では臨界的な伝播を与える。

####クラスターが疎で初期感染者が下辺に偏在し速度分布に偏りがある場合
文字で書くと複雑な設定に見えるけど現実にはありそうな話です。
#####パラメータ；①Rc=20x20, ②p=0.5, ③d=4x250/200x200, ④mr=7.5, ⑤mx=1, my=0.5
最後の⑤のパラメータで速度をx方向を大きく、y方向の速度を半分にしています。

#####パラメータ；①Rc=20x20, ②p=0.0075, ③d=4x250/200x200, ④mr=7.5, ⑤mx=1, my=0.75
初期の感染者を左下の村に最初から置くと以下のようになります。そして、感染確率も0.0075（一様分布の臨界的な伝播の値）という小さな値にします。
y方向の速度を0.75にしても、上の村に属する感受性保持者はほぼ感染しないことが分かります。

この領域では、感染数の推移はとても構造を持った複雑になるということが分かる。
###議論
大切なことはSIRモデルで見るような綺麗な感染数の増減ではない領域がある。そして、臨界的な伝播の場合、
①感染者や感受性保持者の移動速度を小さければ感染伝播しない
②感染率を下げると感染伝播は発生しない
③感受性保持者の密度が小さければ感染しない
ということで、感染者と感受性保持者を混ぜる（旅行する）のは危険である。
感染率は、（感染者が）マスクや（感受性保持者が）手洗いなどで０に出来る可能性がある。まあ、無症状感染者もいるのでマスクは全員する必要がある。
さらに、何より密度を一定以下に下げるだけで感染は防げるので集まらないなど、３蜜（crowded, closely, closed）を回避する施策の実践。
という巷で言われている施策の実践で感染回避が出来そうである。
最後に、計算で示してはいませんが、いわゆる集団免疫という概念は、こういう計算では単に密度や感染確率を変化させると、変化するのでほぼ意味の無い概念であることが理解できました。
###まとめ
・臨界的な伝播領域ではSIRモデルは成立しないようだ
・クラスター的な感受性保持者分布では構造的な感染者発生になる
・クラスターが村的に分布する場合、村と村の感染伝播は起こりにくい
・臨界的な伝播の場合、速度が小さいと伝播しない
・集団免疫という概念は意味が無い（感染確率、感染範囲、感受性保持者分布などに依存して変化する）

・感染者に速度を持たせれば当然感染増加速度が速くなるが、シミュレーションしきれていない
###おまけ

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import random
import time

PARTICLE_NO = 1000 # 粒子数
ITERATION = 200 # 最大ループ回数　感染者が０になると止まる
MIN_X, MIN_Y = -100.0, -100.0 # 探索開始時の範囲最小値
MAX_X, MAX_Y = 100.0, 100.0 # 探索開始時の範囲最大値
recovery=30 #一定時間経過したら治癒
p=0.5 #0.03 #probability of infecion
rc=100 #121 #169 #225　感染する範囲円の半径^2

start = time.time()

def plot_particle(sk,positions,elt,r,g,b):
    #fig, ax = plt.subplots()
    el_time = time.time()-start
    fig, (ax1, ax2) = plt.subplots(1, 2, sharey=False,figsize=(8*2, 8))
        
    for j in range(0,PARTICLE_NO):
        x=positions[j]["x"]
        y=positions[j]["y"]
        c=positions[j]["c"]
        s = 5**2
        ax1.scatter(x, y, s, c, marker="o")
    ax1.set_xlim([MIN_X, MAX_X])
    ax1.set_ylim([MIN_Y, MAX_Y])
    ax1.set_xlabel("x")
    ax1.set_ylabel("y")
    ax1.set_title("{:.2f}:InfectionRate;{:.2f} %".format(el_time,(PARTICLE_NO-b[-1])/PARTICLE_NO*100))
    
    ind = np.arange(len(elt))  # the x locations for the groups
    width = 0.3       # the width of the bars

    ax2.set_ylim([0, PARTICLE_NO])
    ax2.set_title("{:.2f}:red_{} green_{} blue_{}".format(el_time,r[-1],g[-1],b[-1]))
    rect1 = ax2.bar(ind, b,width, color="b")
    rect2 = ax2.bar(ind+width, g, width, color="g") #, bottom=b)
    rect3 = ax2.bar(ind+2*width, r,width, color="r") #, bottom=b)
    plt.pause(0.1)
    plt.savefig('./fig/fig{}_.png'.format(sk)) 
    plt.close()
    
# 粒子の位置更新関数
def update_position(positions,velocity):
    x0 = []
    y0 = []
    for i in range(PARTICLE_NO):
        c=positions[i]["c"]
        t_time = positions[i]["t"]  #初期値０，感染は感染時時間
        k_time = time.time()-start  #経過時間
        s = positions[i]["flag"]    #感染なし０，感染：１
        if s == 1 and c == "red":   #感染済な場合
            if k_time-t_time>recovery:  #一定時間経過したら治癒
                #print("inside",i,s,c,k_time-t_time)
                c = "blue"
                positions[i]["c"] = "green"
                positions[i]["flag"] = 1   #ただし、感染履歴ありのまま
        if c == "red":  #感染redなら位置情報取得
            x0.append(positions[i]["x"])
            y0.append(positions[i]["y"])
            #print("4",i,s,c,t_time)
    #print(x0,y0)   
    position = []
    for j in range(PARTICLE_NO):
        x=positions[j]["x"]
        y=positions[j]["y"]
        c=positions[j]["c"]
        s = positions[j]["flag"]
        t_time = positions[j]["t"]
        for k in range(len(x0)):
            if (x-x0[k])**2+(y-y0[k])**2 < rc and random.uniform(0,1)<p:
                if s ==0:
                    c = "red"
                    t_time = time.time()-start
                    s = 1
                    positions[j]["flag"]=s
                else:
                    continue
        vx = velocity[j]["x"]+1.085*random.uniform(-1, 1) #係数が粒子の運動性の大きさ
        vy = velocity[j]["y"]+1.085*random.uniform(-1, 1)
        new_x = x + vx
        new_y = y + vy
        p_color = c
        s=s
        
        position.append({"x": new_x, "y": new_y, "c": p_color, "t": t_time,"flag":s})
        velocity.append({"x": vx, "y": vy})

    return position, velocity, x0

def count_brg(position):
    r=0
    g=0
    b=0
    for j in range(len(position)):
        if position[j]["c"] == "red":
            r += 1
        elif position[j]["c"] == "green":
            g += 1
        else:
            b += 1
    return r,g,b        

def main():
    # 時間計測開始
    #start = time.time()
    xy_min, xy_max = -32, 32
    # 各粒子の初期位置, 速度, personal best, global best 及びsearch space設定
    position = []
    velocity = []  #速度は使えるように拡張
    
    # 初期位置, 初期速度
    #position.append({"x": random.uniform(MIN_X, MAX_X), "y": random.uniform(MIN_Y, MAX_Y), "c": "red", "t":0, "flag":1})
    position.append({"x": 0, "y": 0, "c": "red", "t":0, "flag":1}) #真ん中（0，0）に初期感染者を1人置く
    velocity.append({"x": 0, "y": 0}) #感染者の初速度0としている
    for i in range(0,10): #x方向に集落を10個並べる
        for j in range(0,10): #xy方向に集落をメッシュ10で並べる
            for k in range(0,10): #１集落辺り10個の感受性保持者を分布
                s=k+j*10+i*100;
                position.append({"x": 10+(-100+i*20)+random.uniform(MIN_X/100, MAX_X/100), "y":10+(-100+j*20)+ random.uniform(MIN_Y/100, MAX_Y/100), "c": "blue", "t": 0, "flag":0})
                velocity.append({"x": 0, "y": 0})
    print(len(position))
    sk = 0
    red=[]
    green=[]
    blue=[]
    elapsed_time = []
    while sk < ITERATION:
        position,velocity, x0 = update_position(position,velocity) ######
        r,g,b = count_brg(position)
        red.append(r)
        green.append(g)
        blue.append(b)
        el_time=time.time()-start
        elapsed_time.append(el_time)
        #print("{:.2f}:red_{} green_{} blue_{}".format(el_time,r,g,b))
        plot_particle(sk,position,elapsed_time,red,green,blue)
        if x0==[]:
            break
        sk += 1

    # 時間計測終了
    process_time = time.time() - start
    print("time:", process_time)
    
if __name__ == '__main__':
    main()