谢林隔离模型Python代码问题：不满主体与空位置交换异常

问题背景 #

实现谢林隔离模型时，矩阵中0代表空房屋，1、2代表不同族群。当主体的相似邻居占比低于阈值self.par时，该主体成为不满主体，需交换至空房屋位置。预期多步迭代后，frac_mean（相似邻居占比的均值）指标应持续下降，但实际未呈现稳定下降趋势，推测核心问题为不满主体与空房屋的交换逻辑错误。

已尝试方案 #

未直接修改原self.array，而是复制数组完成交换后赋值给原数组，但问题仍存在。

原代码 #

import numpy as np
from scipy.signal import correlate2d

class Schelling():
    kernel = [[1,1,1],[1,0,1],[1,1,1]]
    #par = 0.3
    
    def __init__(self, n, par=0.3):
        self.par=par
        probs = [0.1, 0.45, 0.45]
        choices = [0, 1, 2]
        self.array = np.random.choice(choices, (n, n), p=probs)

    def count_neighbours(self):
        a = self.array
        
        empty = a == 0
        red = a == 1
        blue = a == 2

        num_red = correlate2d(red, self.kernel, mode='same', boundary='wrap')
        num_blue = correlate2d(blue, self.kernel, mode='same', boundary='wrap')
        num_neighbours = num_red + num_blue

        frac_red = num_red / num_neighbours
        frac_blue = num_blue / num_neighbours

        frac_red[num_neighbours == 0] = 0
        frac_blue[num_neighbours == 0] = 0

        # 向量式if-else应用
        frac_same = np.where(red, frac_red, frac_blue)

        # 修正空房屋位置的frac_same值
        frac_same[empty] = np.nan

        return empty, frac_red, frac_blue, frac_same, a
    
    def step(self):
        empty, frac_red, frac_blue, frac_same, count_neighbours_list = self.count_neighbours()
        metric=np.nanmean(frac_same)
        unhappy_address = list(zip(*np.array(np.nonzero(frac_same < self.par))))
        np.random.shuffle(unhappy_address)
        empty_address = list(zip(*np.array(np.nonzero(empty))))

        
        # 执行交换直到无法继续
        unhappy_copy=unhappy_address.copy()
        empty_copy=empty_address.copy()
        
        ind=len(unhappy_copy)
        
        #ind=min(len(unhappy_address), len(empty_address))
        for i in range(ind):
            # 添加索引越界检查
            if i == len(empty_address):
                
                break
                
            else:
            
                unhappy_tup_req=unhappy_copy[i]
                emp_tup_req=empty_copy[i]
                
                #count_neighbours_list[emp_tup_req]=count_neighbours_list[unhappy_tup_req]
                #count_neighbours_list[unhappy_tup_req]==0
                
                count_neighbours_list[emp_tup_req], count_neighbours_list[unhappy_tup_req] = count_neighbours_list[unhappy_tup_req], count_neighbours_list[emp_tup_req]

            
        self.array= count_neighbours_list
                

        return unhappy_address, empty_address, count_neighbours_list, metric

问题分析 #

1. 数组引用错误：count_neighbours()返回的a是原self.array的引用，直接对其进行交换操作会修改原数组，导致后续邻居计算基于已修改的数组，逻辑混乱。
2. 配对逻辑不合理：按索引顺序一一配对不满主体和空房屋，未考虑空房屋数量可能少于不满主体数量，且不符合随机分配空房屋的模型逻辑。
3. 循环终止条件错误：if i == len(empty_address): break的判断时机滞后，当i等于空房屋长度时，已经超出空房屋列表的索引范围，会触发索引越界问题。

修复方案与代码 #

修复点 #

• 对原数组进行深拷贝，避免修改原数组的同时影响后续计算
• 随机打乱不满主体和空房屋列表，取两者长度的最小值进行随机配对
• 简化循环逻辑，直接遍历配对后的坐标执行交换

修复后代码 #

import numpy as np
from scipy.signal import correlate2d

class Schelling():
    kernel = [[1,1,1],[1,0,1],[1,1,1]]
    
    def __init__(self, n, par=0.3):
        self.par = par
        probs = [0.1, 0.45, 0.45]
        choices = [0, 1, 2]
        self.array = np.random.choice(choices, (n, n), p=probs)

    def count_neighbours(self):
        a = self.array
        
        empty = a == 0
        red = a == 1
        blue = a == 2

        num_red = correlate2d(red, self.kernel, mode='same', boundary='wrap')
        num_blue = correlate2d(blue, self.kernel, mode='same', boundary='wrap')
        num_neighbours = num_red + num_blue

        frac_red = num_red / num_neighbours
        frac_blue = num_blue / num_neighbours

        frac_red[num_neighbours == 0] = 0
        frac_blue[num_neighbours == 0] = 0

        frac_same = np.where(red, frac_red, frac_blue)
        frac_same[empty] = np.nan

        return empty, frac_red, frac_blue, frac_same, a
    
    def step(self):
        empty, frac_red, frac_blue, frac_same, _ = self.count_neighbours()
        metric = np.nanmean(frac_same)
        
        # 获取不满主体与空房屋的坐标列表
        unhappy_address = list(zip(*np.nonzero(frac_same < self.par)))
        empty_address = list(zip(*np.nonzero(empty)))
        
        # 深拷贝原数组，避免直接修改原数组
        new_array = self.array.copy()
        
        # 确定可交换的最大数量
        swap_count = min(len(unhappy_address), len(empty_address))
        
        # 随机打乱列表，实现随机配对
        np.random.shuffle(unhappy_address)
        np.random.shuffle(empty_address)
        
        # 执行交换：将不满主体移至空房屋，原位置设为0
        for unhappy_tup, emp_tup in zip(unhappy_address[:swap_count], empty_address[:swap_count]):
            new_array[emp_tup] = new_array[unhappy_tup]
            new_array[unhappy_tup] = 0
        
        # 更新原数组
        self.array = new_array

        return unhappy_address, empty_address, new_array, metric

验证说明 #

修复后，每次迭代会将随机选择的不满主体移至随机空房屋位置，frac_mean指标会随着迭代逐步下降，最终趋于稳定，符合谢林隔离模型的预期行为。

内容的提问来源于stack exchange，提问作者Sl30202