投稿者: 
さらに発展的なガウゼの法則をpythonで学ぶ
以前の投稿でガウゼの法則のシミュレーションを行った．
今回はそこにさらに時間差という要素を加えてみた．
コードは下記．
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import imageio
from tqdm import tqdm
# パラメータの設定
GRID_SIZE = 50 # グリッドのサイズ
initial_agents_A = 300 # エージェントAの初期数
initial_agents_B = 300 # エージェントBの初期数
p_reproduce_A = 0.04 # エージェントAの繁殖確率
p_reproduce_B = 0.035 # エージェントBの繁殖確率
p_die = 0.01 # 死亡確率
num_time_steps = 200 # シミュレーションのステップ数
# 資源に関するパラメータ
RESOURCE_MAX = 10 # 各セルの最大資源量
RESOURCE_REGEN = 0.1 # 資源の再生率
RESOURCE_CONSUME_A = 2 # エージェントAが1回のアクションで消費する資源量
RESOURCE_CONSUME_B = 2.5 # エージェントBが1回のアクションで消費する資源量
# 日中と夜間のサイクルに関するパラメータ
DAY_LENGTH = 10 # 日中のステップ数
NIGHT_LENGTH = 10 # 夜間のステップ数
# グリッドの初期化
grid_agents = np.zeros((GRID_SIZE, GRID_SIZE), dtype=int) # 0: Empty, 1: Agent A, 2: Agent B
grid_resources = np.random.uniform(5, RESOURCE_MAX, (GRID_SIZE, GRID_SIZE)) # 資源量の初期化
# エージェントの初期配置
available_positions = [(i, j) for i in range(GRID_SIZE) for j in range(GRID_SIZE)]
initial_positions = np.random.choice(len(available_positions), initial_agents_A + initial_agents_B, replace=False)
positions_A = [available_positions[idx] for idx in initial_positions[:initial_agents_A]]
positions_B = [available_positions[idx] for idx in initial_positions[initial_agents_A:initial_agents_A + initial_agents_B]]
for pos in positions_A:
grid_agents[pos] = 1
for pos in positions_B:
grid_agents[pos] = 2
# 隣接セルを取得する関数
def get_neighbors(i, j, grid_size):
neighbors = []
for x in [-1, 0, 1]:
for y in [-1, 0, 1]:
if x == 0 and y == 0:
continue
ni, nj = i + x, j + y
if 0 <= ni < grid_size and 0 <= nj < grid_size:
neighbors.append((ni, nj))
return neighbors
# カラーマップの設定（0: White, 1: Blue, 2: Red）
cmap_agents = mcolors.ListedColormap(['white', 'blue', 'red'])
cmap_resources = plt.cm.Greens
# エージェント数の履歴を記録
history_A = []
history_B = []
# シミュレーションとGIFの作成
with imageio.get_writer('extended_gause_simulation_day_night.gif', mode='I', fps=10) as writer:
for t in tqdm(range(num_time_steps), desc="Simulating"):
new_grid_agents = grid_agents.copy()
new_grid_resources = grid_resources.copy()
# 日中か夜間かを判定
cycle_length = DAY_LENGTH + NIGHT_LENGTH
cycle_step = t % cycle_length
if cycle_step < DAY_LENGTH:
current_phase = 'Day'
active_agent = 1 # エージェントAが活動
else:
current_phase = 'Night'
active_agent = 2 # エージェントBが活動
# エージェントの位置をリストアップ
if active_agent == 1:
agent_positions = list(zip(*np.where(grid_agents == 1)))
else:
agent_positions = list(zip(*np.where(grid_agents == 2)))
np.random.shuffle(agent_positions) # ランダムな順序で処理
for pos in agent_positions:
i, j = pos
agent_type = grid_agents[i, j]
# 死亡判定
if np.random.rand() < p_die:
new_grid_agents[i, j] = 0
continue
# 資源の消費
if agent_type == 1:
consume = RESOURCE_CONSUME_A
else:
consume = RESOURCE_CONSUME_B
available_resource = new_grid_resources[i, j]
if available_resource >= consume:
new_grid_resources[i, j] -= consume
consumed = consume
elif available_resource > 0:
new_grid_resources[i, j] = 0
consumed = available_resource
else:
consumed = 0 # 資源がない場合
# 繁殖判定（資源を消費した場合のみ）
if consumed > 0 and np.random.rand() < (p_reproduce_A if agent_type == 1 else p_reproduce_B):
neighbors = get_neighbors(i, j, GRID_SIZE)
np.random.shuffle(neighbors)
for ni, nj in neighbors:
if new_grid_agents[ni, nj] == 0 and new_grid_resources[ni, nj] > 0:
new_grid_agents[ni, nj] = agent_type
break
# 移動判定
if np.random.rand() < 0.3: # 移動確率
neighbors = get_neighbors(i, j, GRID_SIZE)
# 資源が多いセルを優先
neighbors_sorted = sorted(neighbors, key=lambda x: new_grid_resources[x], reverse=True)
for ni, nj in neighbors_sorted:
if new_grid_agents[ni, nj] == 0:
new_grid_agents[ni, nj] = agent_type
new_grid_agents[i, j] = 0
break
# 資源の再生
new_grid_resources += RESOURCE_REGEN
new_grid_resources = np.clip(new_grid_resources, 0, RESOURCE_MAX)
grid_agents = new_grid_agents.copy()
grid_resources = new_grid_resources.copy()
# エージェント数の記録
count_A = np.sum(grid_agents == 1)
count_B = np.sum(grid_agents == 2)
history_A.append(count_A)
history_B.append(count_B)
# グリッドのプロット
fig, ax = plt.subplots(figsize=(6, 6))
ax.imshow(grid_resources, cmap=cmap_resources, vmin=0, vmax=RESOURCE_MAX, alpha=0.3)
ax.imshow(grid_agents, cmap=cmap_agents, vmin=0, vmax=2, alpha=0.7)
ax.set_title(f"Time Step: {t+1} - {current_phase}")
ax.axis('off')
plt.tight_layout()
# フレームをGIFに追加
fig.canvas.draw()
image = np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8')
image = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))
writer.append_data(image)
plt.close()
# 最終的なエージェント数のプロット
plt.figure(figsize=(8, 4))
plt.plot(history_A, label='Agent A', color='blue')
plt.plot(history_B, label='Agent B', color='red')
plt.xlabel('Time Steps')
plt.ylabel('Number of Agents')
plt.title('Population Dynamics Over Time')
plt.legend()
plt.tight_layout()
plt.savefig('agent_population_history_day_night.png')
plt.close()
print("シミュレーションとGIFの作成が完了しました。")
