ガウゼの法則のシミュレーション 2

発展的なガウゼの法則を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.02  # エージェントBの繁殖確率（Aよりわずかに低い）
p_die = 0.02  # 死亡確率
num_time_steps = 300  # シミュレーションのステップ数

# 資源に関するパラメータ
RESOURCE_MAX = 10  # 各セルの最大資源量
RESOURCE_REGEN = 0.1  # 資源の再生率
RESOURCE_CONSUME_A = 2  # エージェントAが1回のアクションで消費する資源量
RESOURCE_CONSUME_B = 2.5  # エージェントBが1回のアクションで消費する資源量

# グリッドの初期化
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の作成
frames = []
with imageio.get_writer('extended_gause_simulation.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()

        # エージェントの位置をリストアップ
        agent_positions = list(zip(*np.where(grid_agents > 0)))
        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}")
        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.png')
    plt.close()

print("シミュレーションとGIFの作成が完了しました。")