エージェントベースモデルの改善

前提

エージェントベースモデルを動かしていて終了時の時間を取得したいのですが、時間を取得するためにどうした上手くいくのかが分からず詰まってしまっています。

実現したいこと

シミュレーション開始から、エージェントが壁に到達した時の時間をカウントして保存し、最後に出力する様にしたいです。

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

該当のソースコード

python
1import numpy as np
2import matplotlib.pyplot as plt
3import matplotlib.animation as animation
4from IPython.display import HTML
5
6a1 = -0.034298
7b1 = 3.348394
8c1 = 4.252840
9d1 = -0.003423
10
11class Agent:
12    def __init__(self, name = None, pos = None, vel = None, SIZE = 10, VEL = 1,
13            RANGE_SEPARATE = 0.5,color = None):
14        self.name = name
15        self.pos = pos if pos else np.random.uniform(-SIZE, SIZE, 2)
16        self.vel = vel if vel else np.random.uniform(-VEL, VEL, 2)
17        self.SIZE = SIZE
18        self.VEL = VEL
19        self.RANGE_SEPARATE = RANGE_SEPARATE#分離範囲
20        self.vel_tmp = np.zeros(2)#移動
21        self.color = color if color else 'k'
22
23    
24    def distance(self, other):#ユークリッド距離
25        return np.linalg.norm(self.pos - other.pos)
26    
27    def ruleSeparate(self, others, ratio = 1):
28        others = [other for other in others if self.distance(other) < self.RANGE_SEPARATE]#範囲内に入ったエージェントの選別
29        v = np.zeros(2)
30        if not others:
31            return 0
32        for other in others:
33            d = self.pos - other.pos#dが-になる時ある？
34            v += d / self.distance(other)#進行方向の単位ベクトル
35        self.vel_tmp += v / len(others) * ratio#平均を出してそれをratio倍する
36        return v / len(others) * ratio
37    
38    
39    
40    def calculate(self, others):
41        self.ruleSeparate(others, 0.5)
42        
43    def update(self) -> bool:#返り値の型名
44        self.vel += self.vel_tmp#進行速度ベクトル更新
45        v = np.linalg.norm(self.vel)
46       #壁に当たったら反転
47        if (abs(self.pos + self.vel)[0] > self.SIZE):
48            self.vel[0] = -self.vel[0]
49        if (abs(self.pos + self.vel)[1] > self.SIZE):
50            self.vel[1] = -self.vel[1]
51        
52        self.pos += self.vel#位置更新
53        self.vel_tmp = np.zeros(2)##進行速度ベクトル初期化
54        return True
55
56class AgentEdited(Agent):
57    def __init__(self, *args, **kwargs):
58        super(AgentEdited, self).__init__(*args, **kwargs)
59        #クラス(子クラス)で別のクラス(親クラス)を継承できます。継承することで、親クラスのメソッドを子クラスから呼び出せる
60        #super().親クラスのメソッド # python3系での標準の書き方
61        #super(親クラスのオブジェクト, self).親クラスのメソッド # python2系での書き方
62        #不特定多数の引数を受け取ることができる仕組みを作るのが「*args」「**kwargs」
63        self.color = kwargs.get('color', 'r')
64    
65    def ruleSeparate(self, others, ratio = 1):
66        others = [other for other in others if self.distance(other) < self.RANGE_SEPARATE]
67        v  = np.zeros(2)
68        if not others:
69            return 0
70        for other in others:
71            d = self.pos - other.pos
72            vel = self.vel - other.vel
73            TTNP = d / (vel + 1e-5) # 0除算防止のため微小値を挿入
74            DINP = self.distance(other)
75            v += 1/(1+np.exp(-(c1+d1*TTNP))) * 1/(1+np.exp(-(b1+a1*DINP)))
76        self.vel_tmp += v / len(others) * ratio
77        return v / len(others) * ratio
78    
79    def calculate(self, others):
80        self.ruleSeparate(others) # only separate
81
82def calc_rad(pos2, pos1):#角計算
83    return np.arctan2(pos2[1] - pos1[1], pos2[0] - pos1[0])
84
85#np.arctan(x)は引数が一つでarctan(x)をラジアンで返す。返り値は-pi / 2からpi / 2（-90度から90度）の間になる。
86#np.arctan2(y, x)は引数が二つで、arctan(y / x)をラジアンで返す。この角度は、極座標平面において原点から座標(x, y)へのベクトルがx軸の正の方向となす角度（偏角）であり、返り値は-piからpi（-180度から180度）の間になる
87#print(np.degrees(np.arctan2(-1, 1)))
88# -45.0
89#print(np.degrees(np.arctan2(1, -1)))
90# 135.0
91
92def rotate_vec(vec, rad):#ベクトル更新
93    return np.dot(vec, np.array([[np.cos(rad), -np.sin(rad)], [np.sin(rad), np.cos(rad)]]).T)
94
95class AgentGoForWall(AgentEdited):
96    def __init__(self, *args, **kwargs):
97        super(AgentGoForWall, self).__init__(*args, **kwargs)
98        self.color = kwargs.get('color', 'b')
99        max_distance = 0
100        self.goal = None # 角が一番遠い場所
101        for px in [-self.SIZE, self.SIZE]:
102            for py in [-self.SIZE, self.SIZE]:
103                dist = np.linalg.norm(self.pos - np.array([px, py]))
104                if max_distance < dist:
105                    max_distance = dist#最大値の更新
106                    self.goal = np.array([px, py])
107        self.vel = rotate_vec(np.array([self.VEL, 0]), calc_rad(self.goal, self.pos))
108
109    def update(self) -> bool:#"返り値の型名"
110        self.vel += self.vel_tmp
111        self.vel += rotate_vec(np.array([self.VEL, 0]), calc_rad(self.goal, self.pos))
112        
113        v = np.linalg.norm(self.vel)
114       
115        if abs(self.pos + self.vel)[0] > self.SIZE or abs(self.pos + self.vel)[1] > self.SIZE:
116            return False # deactivate this Agent
117        
118        self.pos += self.vel
119        self.vel_tmp = np.zeros(2)
120        return True
121
122class Boids:
123    def __init__(self, AGENT = 20):#Agent数
124        self.agents = [Agent() for _ in range(AGENT)]
125    
126    def calculate(self):
127        for agent in self.agents:
128            agent.calculate([other for other in self.agents if agent != other])
129    
130    def update(self):
131        self.agents = [agent for agent in self.agents if agent.update()]
132    
133    def simulation(self):
134        self.calculate()
135        self.update()
136    
137    def positions(self):#showimageはpos_arrayを渡してた
138      x, y, c = list(), list(), list()
139      for agent in self.agents:
140        x.append(agent.pos[0])
141        y.append(agent.pos[1])
142        c.append(agent.color)
143      return (np.array(x), np.array(y)), c
144
145np.random.seed( 0 )
146B = Boids(AGENT = 2)
147B.agents.append(AgentEdited())
148#B.agents.append(Agent("Crazy", RANGE_SEPARATE = 1.5))
149num=50
150for i in range(num):
151  B.agents.append(AgentGoForWall())
152
153fig, ax = plt.subplots(figsize = (6, 6))
154ax.set_xlim(-10, 10)
155ax.set_ylim(-10, 10)
156ax.grid(True)
157ims = []
158for t in range(1000):
159    B.simulation()
160    plot_data,plot_color = B.positions()
161    im=ax.scatter(*plot_data, marker="o", s=20, c=plot_color)
162    ims.append([im])
163    #print(plot_data)
164    #print(plot_color)
165
166ani = animation.ArtistAnimation(fig, ims, interval=400, repeat=False)
167HTML(ani.to_jshtml())

試したこと

class AgentGoForWall(AgentEdited):
    def __init__(self, *args, **kwargs):
        super(AgentGoForWall, self).__init__(*args, **kwargs)
        self.color = kwargs.get('color', 'b')
        max_distance = 0
        self.goal = None # 角が一番遠い場所
        for px in [-self.SIZE, self.SIZE]:
            for py in [-self.SIZE, self.SIZE]:
                dist = np.linalg.norm(self.pos - np.array([px, py]))
                if max_distance < dist:
                    max_distance = dist#最大値の更新
                    self.goal = np.array([px, py])
        self.vel = rotate_vec(np.array([self.VEL, 0]), calc_rad(self.goal, self.pos))

    def update(self) -> bool:#"返り値の型名"
        self.vel += self.vel_tmp
        self.vel += rotate_vec(np.array([self.VEL, 0]), calc_rad(self.goal, self.pos))
        
        v = np.linalg.norm(self.vel)
       
        if abs(self.pos + self.vel)[0] > self.SIZE or abs(self.pos + self.vel)[1] > self.SIZE:
            return False # deactivate this Agent

cntという空のリストを作って、ここのdeacticvate時に追加すればいいと思ったのですが、カウントを開始する際にどこで開始すれば上手くいくのかが分からなかったです。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

ps_aux_grep

2022/11/17 09:39

ここで言う「時間」を定義しないといけないと思います．何度のupdateが実行されたか，を1単位時間にしたいのでしょうか？それともプログラムの動作時間を利用したいのでしょうか？

lack

2022/11/17 16:37 編集

時間の定義が明瞭でなくて申し訳ありませんでした。 updateが実行された回数を単位時間としてカウントして、各エージェントのゴール到達時間を空のリストに格納したいです。 def update(self) -> bool: cnt_collect=[] ... cnt+=1 if abs(self.pos + self.vel)[0] > self.SIZE or abs(self.pos + self.vel)[1] > self.SIZE: cnt_collect.append(cnt) この様な感じでリストに格納して最後にそれを出力出来る様にしたいです。

行動規範の内容に同意します

回答1件

ベストアンサー

あんまり良い解決手法が思いつかなかったので苦肉の策になります．グローバル変数survivalを辞書型で宣言して扱うようにします．

Python
1import numpy as np
2
3survival = dict() # グローバルで用意
4
5class Agent:
6
7# (略)
8
9class AgentGoForWall(AgentEdited):
10    def __init__(self, *args, **kwargs):
11        super(AgentGoForWall, self).__init__(*args, **kwargs)
12        self.cnt = 0 # カウント初期化
13
14        # (略)
15
16    def update(self) -> bool:
17        self.cnt += 1 # カウントアップ
18        self.vel += self.vel_tmp
19        self.vel += rotate_vec(np.array([self.VEL, 0]), calc_rad(self.goal, self.pos))
20        
21        v = np.linalg.norm(self.vel)
22       
23        if abs(self.pos + self.vel)[0] > self.SIZE or abs(self.pos + self.vel)[1] > self.SIZE:
24            survival[id(self)] = self.cnt # カウント数を記録
25            return False # deactivate this Agent
26
27# (略)
28
29ani = animation.ArtistAnimation(fig, ims, interval=400, repeat=False)
30HTML(ani.to_jshtml())
31
32import _ctypes
33for ptr, life in survival.items():
34    agent = _ctypes.PyObj_FromPtr(ptr) # idから復元
35    print(agent.name + ": " + life)

このように書いておけば，ライブラリ_ctpyesを使ってidからAgentオブジェクトを復元できます．
最後にカウントを挿入するのではなく，ログ的に記録し続けることもできます．

Python
1class AgentGoForWall(AgentEdited):
2    def __init__(self, *args, **kwargs):
3        super(AgentGoForWall, self).__init__(*args, **kwargs)
4        survival[id(self)] = 0 # カウント初期化
5
6        # (略)
7
8    def update(self) -> bool:
9        survival[id(self)] += 1 # カウントアップ
10        self.vel += self.vel_tmp
11        self.vel += rotate_vec(np.array([self.VEL, 0]), calc_rad(self.goal, self.pos))
12        
13        v = np.linalg.norm(self.vel)
14       
15        if abs(self.pos + self.vel)[0] > self.SIZE or abs(self.pos + self.vel)[1] > self.SIZE:
16            return False # deactivate this Agent

こちらの方であれば，壁までの到達時間の記録だけでなく，座標を記録させ続けることだってできます．今回は単純に0で初期化して1ずつ増やしましたが，survival[id(self)] = [self.pos]で初期化してupdate毎にsurvival[id(self)].append(self.pos)で座標の記録が可能です．最終的にlen(survival[id(self)])が記録した座標の個数，すなわち壁に到達するまでにかかった時間になります．

Python
1class AgentGoForWall(AgentEdited):
2    def __init__(self, *args, **kwargs):
3        super(AgentGoForWall, self).__init__(*args, **kwargs)
4        survival[id(self)] = [self.pos] # カウント初期化
5
6        # (略)
7
8    def update(self) -> bool:
9        self.vel += self.vel_tmp
10        self.vel += rotate_vec(np.array([self.VEL, 0]), calc_rad(self.goal, self.pos))
11        
12        v = np.linalg.norm(self.vel)
13       
14        if abs(self.pos + self.vel)[0] > self.SIZE or abs(self.pos + self.vel)[1] > self.SIZE:
15            return False # deactivate this Agent
16        
17        self.pos += self.vel
18        survival[id(self)].append(self.pos) # カウントアップ
19        self.vel_tmp = np.zeros(2)
20        return True
21
22# (略)
23
24ani = animation.ArtistAnimation(fig, ims, interval=400, repeat=False)
25HTML(ani.to_jshtml())
26
27import _ctypes
28for ptr, log in survival.items():
29    agent = _ctypes.PyObj_FromPtr(ptr) # idから復元
30    print(agent.name + ": " + len(log))