python中的mesa包預製了一些類，提供了一些基礎模型，可以大大簡化abm建模的工作量。在python中實現，也有利於和其它演算法相結合。本文是一次作業，按照個人理解把mesa包教程整理，程式碼壓縮成了兩大部分。如果是新手上手，建議檢視下方的官方連結，在jupyter notebook裡一步一步執行程式碼。

mesa包的英文簡介:

模型構造

建模一個簡單的模型：
強假設（只要有錢就把錢給鄰居）下財富的流動，觀察基尼係數變化

主要的類與類圖

一個未明確任何功能的，簡單的agent模型，含有如下類，分別是：

• 例項化的模型，對於每一個具體問題應當是唯一的
• 主體或智慧體（agent），在啟動模型前應對其進行例項化，通常一個模型有多個agent（實現觀察湧現性等功能）
• 排程器（schedule），控制每一個時間步（ticks）裡agent的動作（step()）順序，如同時（即並行），隨機序列執行。隨機序列是最常見也是最簡單的控制方式
• 網格（grids）提供兩種SingleGrid和MultiGrid，前者一個格子空間上只能有一個agent，而後者一個格子上可以由多個agent
• 資料收集類(DataCollector)，收集模型級變數, agent級別變數和其它變數（下圖省略）
• 控制模型執行類(BatchRunner),使用固定的引數多次生成和執行模型，檢視結果的情況。可消除隨機誤差（下圖省略）
  類圖：
  

程式碼

不考慮agent的位置，隨機把錢給其它agent

注意，下面的程式碼使用的是jupyter notebook，可以複製到自己的jupyter裡面按順序執行

from mesa import Agent, Model
from mesa.time import RandomActivation

先定義MoneyModel類，因為下面MoneyAgent類實現時要呼叫，具體參見上方類圖，二者存在一對多的組合關係

class MoneyModel(Model):
    """A model with some number of agents."""
    def __init__(self, N):
        self.num_agents = N
        self.schedule = RandomActivation(self)
        # Create agents
        for i in range(self.num_agents):
            a = MoneyAgent(i, self)
            self.schedule.add(a)

    def step(self):
        '''Advance the model by one step.'''
        self.schedule.step()

定義MoneyAgent類

class MoneyAgent(Agent):
    """ An agent with fixed initial wealth."""
    def __init__(self, unique_id, model):
        super().__init__(unique_id, model)
        self.wealth = 1
    def step(self):
        # The agent's step will go here.
        # For demonstration purposes we will print the agent's unique_id
        if self.wealth == 0:
            return
        other_agent = self.random.choice(self.model.schedule.agents) # 需要隨機從所有agent中選擇，同時可能會選中自己
        other_agent.wealth += 1
        self.wealth -= 1

建立10個agent，並執行10步

model = MoneyModel(10)
for i in range(10):
    model.step()

輸出結果：

# For a jupyter notebook add the following line:
%matplotlib inline

# The below is needed for both notebooks and scripts
import matplotlib.pyplot as plt

agent_wealth = [a.wealth for a in model.schedule.agents]
plt.hist(agent_wealth)

反覆使用100個模型，每個模型執行10步，統計所有結果的分佈

all_wealth = []
# 反覆使用100個模型，每個模型執行10步，統計所有結果的分佈
for j in range(100):
    # Run the model
    model = MoneyModel(10)
    for i in range(10):
        model.step()

    # Store the results
    for agent in model.schedule.agents:
        all_wealth.append(agent.wealth)

plt.hist(all_wealth, bins=range(max(all_wealth)+1)) # histgram的柱含下界但不含上界

考慮agent移動時

每個agent只給錢到鄰居節點。如何快速獲得鄰居節點，收集資料，反覆執行模型

from mesa.space import MultiGrid
from mesa.datacollection import DataCollector
from mesa.batchrunner import BatchRunner

重新定義類，增添了一些方法
此處進行了修改，定義了 def posible_move(self)和另一種同樣功能的寫法。兩種方法功能都是選擇周圍八個鄰居節點，不過posible_move()選擇的包含自身。在下方進行了展示，可以切換註釋進行嘗試，兩個函式功能和結果完全相同。

class MoneyAgent(Agent):
    """A model with some number of agents."""
        
    def __init__(self, unique_id, model):
        super().__init__(unique_id, model)
        self.wealth = 1
            
    def move(self):
        possible_steps = self.model.grid.get_neighborhood(
            self.pos,
            moore=True,   # moore包含周邊八個結構，如果使用Neumann只包含正交的四個位置
            include_center=False) # 不包含自己
        new_position = self.random.choice(possible_steps)
        # 另一種同樣功能寫法
        # new_position = self.posible_move()
        self.model.grid.move_agent(self, new_position)
        
    def posible_move(self):
        neighbors = []
        x, y = self.pos
        for dx in [-1, 0, 1]:
            for dy in [-1, 0, 1]:
                neighbors.append((x+dx, y+dy))
        return neighbors
    
    def give_money(self):
        cellmates = self.model.grid.get_cell_list_contents([self.pos]) # 判斷是否存在鄰居節點
        if len(cellmates) > 1:  # 如果存在則隨機給鄰居錢
            other = self.random.choice(cellmates)
            other.wealth += 1
            self.wealth -= 1
            
    def step(self):
        self.move()      # f
        if self.wealth > 0:
            self.give_money()

定義基尼係數計算函式

def compute_gini(model):  # 計算基尼係數
    agent_wealths = [agent.wealth for agent in model.schedule.agents]
    x = sorted(agent_wealths)
    N = model.num_agents
    B = sum( xi * (N-i) for i,xi in enumerate(x) ) / (N * sum(x))
    return (1 + (1/N) - 2 * B)

class MoneyModel(Model):
    """A model with some number of agents."""
    def __init__(self, N, width, height):
        self.num_agents = N
        self.grid = MultiGrid(width, height, True)  # True意味著邊界是迴圈的，從一端出去，將從二維平面相反邊進入
        self.schedule = RandomActivation(self)   # 序列隨機呼叫
        self.running = True

        # Create agents
        for i in range(self.num_agents):
            a = MoneyAgent(i, self)
            self.schedule.add(a)
            # Add the agent to a random grid cell
            x = self.random.randrange(self.grid.width)
            y = self.random.randrange(self.grid.height)
            self.grid.place_agent(a, (x, y))

        self.datacollector = DataCollector(
            model_reporters={"Gini": compute_gini}, # 計算模型的基尼係數
            agent_reporters={"Wealth": "wealth"}) 

    def step(self):
        self.datacollector.collect(self)
        self.schedule.step()

生成一個10*10的panel，放入50個agents，根據上方MoneyModel類例項化裡面定義的介面，介面是迴圈的。
執行模型20步

model = MoneyModel(50, 10, 10)
for i in range(20):
    model.step()

計算每個位置上點的密度

import numpy as np

agent_counts = np.zeros((model.grid.width, model.grid.height))
for cell in model.grid.coord_iter():
    cell_content, x, y = cell
    agent_count = len(cell_content)
    agent_counts[x][y] = agent_count
plt.imshow(agent_counts, interpolation='nearest') #根據每個位置（實際上是二維陣列計數）上agents
plt.colorbar()

# If running from a text editor or IDE, remember you'll need the following:
# plt.show()

結果：

與上面類似，變動為讓模型執行100步

model = MoneyModel(50, 10, 10)
for i in range(100):
    model.step()

輸出基尼係數影像

gini = model.datacollector.get_model_vars_dataframe() # 上一個程式碼模組，執行一百步過程中基尼係數的變化
gini.plot()

輸出模型未執行時的點的狀態，只看header（預設ID前5個的agent）

agent_wealth = model.datacollector.get_agent_vars_dataframe()
agent_wealth.head()

結果

檢視執行結束後（第99步執行完成）agent財富的分佈影像

end_wealth = agent_wealth.xs(99, level="Step")["Wealth"]
end_wealth.hist(bins=range(agent_wealth.Wealth.max()+1))

結果

one_agent_wealth = agent_wealth.xs(14, level="AgentID") # 使用angentID14觀察這個agent的財富變化
one_agent_wealth.Wealth.plot()

結果

如果想要讓模型多次執行，除了for迴圈例項化model外，一種更便捷的方法由mesa提供

# 如果想要讓模型多次執行


fixed_params = {"width": 10,
               "height": 10}
variable_params = {"N": range(10, 500, 10)}

batch_run = BatchRunner(MoneyModel,
                        variable_params,
                        fixed_params,
                        iterations=5,
                        max_steps=100,
                        model_reporters={"Gini": compute_gini})
batch_run.run_all()

收集每個模型的最終執行結果

# BatchRunner的資料收集
run_data = batch_run.get_model_vars_dataframe()
run_data.head()
plt.scatter(run_data.N, run_data.Gini)