基于囚徒困境的空间模型python代码

对于基于囚徒困境的空间模型，代码如下：

import numpy as np
import random

class SpatialPD:
    def __init__(self, n_agents, payoff_matrix, topology):
        self.n_agents = n_agents
        self.payoff_matrix = payoff_matrix
        self.topology = topology
        self.game_history = np.zeros((n_agents, n_agents))
        self.strategies = np.zeros((n_agents, 2))
        self.neighbors = [[] for i in range(n_agents)]
        
        self.strategies[:,0] = 1
        self.strategies[:,1] = 0
        
        self.strategies[random.sample(range(n_agents), n_agents//2),0] = 0
        self.strategies[random.sample(range(n_agents), n_agents//2),1] = 1
        
        if topology == "lattice":
            self.lattice_topology()
        elif topology == "random":
            self.random_topology()
    
    def lattice_topology(self):
        for i in range(self.n_agents):
            if i % int(np.sqrt(self.n_agents)) != 0:
                self.neighbors[i].append(i-1)
            if i % int(np.sqrt(self.n_agents)) != int(np.sqrt(self.n_agents))-1:
                self.neighbors[i].append(i+1)
            if i // int(np.sqrt(self.n_agents)) != 0:
                self.neighbors[i].append(i-int(np.sqrt(self.n_agents)))
            if i // int(np.sqrt(self.n_agents)) != int(np.sqrt(self.n_agents))-1:
                self.neighbors[i].append(i+int(np.sqrt(self.n_agents)))
        
    def random_topology(self):
        for i in range(self.n_agents):
            neighbors = random.sample(range(self.n_agents), 4)
            neighbors.remove(i)
            self.neighbors[i] = neighbors
    
    def play_game(self):
        for i in range(self.n_agents):
            for j in self.neighbors[i]:
                if self.strategies[i,0] == 1 and self.strategies[j,0] == 1:
                    self.game_history[i,j] += self.payoff_matrix[0,0]
                    self.game_history[j,i] += self.payoff_matrix[0,0]
                elif self.strategies[i,0] == 1 and self.strategies[j,1] == 1:
                    self.game_history[i,j] += self.payoff_matrix[0,1]
                    self.game_history[j,i] += self.payoff_matrix[1,0]
                elif self.strategies[i,1] == 1 and self.strategies[j,0] == 1:
                    self.game_history[i,j] += self.payoff_matrix[1,0]
                    self.game_history[j,i] += self.payoff_matrix[0,1]
                elif self.strategies[i,1] == 1 and self.strategies[j,1] == 1:
                    self.game_history[i,j] += self.payoff_matrix[1,1]
                    self.game_history[j,i] += self.payoff_matrix[1,1]
    
    def update_strategy(self, beta):
        for i in range(self.n_agents):
            neighbor_payoffs = []
            for j in self.neighbors[i]:
                neighbor_payoffs.append(self.game_history[i,j])
            max_neighbor_payoff = max(neighbor_payoffs)
            if max_neighbor_payoff > self.game_history[i,i]:
                p = 1 / (1 + np.exp(-beta*(max_neighbor_payoff-self.game_history[i,i])))
                if random.uniform(0,1) < p:
                    self.strategies[i,:] = self.strategies[self.neighbors[i][neighbor_payoffs.index(max_neighbor_payoff)],:]
    
    def run(self, n_iterations, beta):
        for i in range(n_iterations):
            self.play_game()
            self.update_strategy(beta)