用python模拟和绘制SAW路径（自回避随机行走问题）

最新推荐文章于 2024-05-12 20:41:29 发布

*1_*

最新推荐文章于 2024-05-12 20:41:29 发布

阅读量1.1k

点赞数 2

文章标签： python 开发语言

本文链接：https://blog.csdn.net/weixin_58627772/article/details/125286091

版权

自回避随机行走问题（SAW），即在格点上不与历史轨迹相撞的随机行走，本文简要说明一种利用python模拟和绘制SAW路径的方法。（完整代码在文章末尾，文字部分是按思考的顺序，代码是最终版本，中间部分代码会出现和文字有出入的内容，建议先看思路）

该路径生成碰到的主要问题是锁死，原因在于路径闭合后随机走向了空间有限的一侧，导致在达到设定步数前出现无路可走的情况。该方法的核心就是在路径闭合时选择闭合圈的外侧，避免进入圈内锁死。

为简单起见所用的是标准的方形网格，可能的前进方向只有4种，排除上一步的落点对应方向，除第一步外（第一步四向随机）可能前进的方向至多三种。设定一个每步更新的可行方向集合，接下来的限制条件均是从中去除某些项。

（Tips:方向可设为ev=[(1,0),(0,1),(-1,0),(0,-1)]，行进用坐标相加，某方向的取反%4，向右（-1）%4，向左（+1）%4，后面会用上）

ev=[(1,0),(0,1),(-1,0),(0,-1)]

为了避免判断历史轨迹需要遍历，将网格赋值为零矩阵，每一步的落点赋值为1，仅需判断周围格点的值，不为0可把对应方向排除。边界设为-1，避免和路径混淆。

order=np.zeros((2*s+3,2*s+3),dtype=int)       #顺序矩阵
direc=[]                            #方向矩阵
for i in np.arange(2*s+3):          #设定边界为-1
    order[0][i]=-1
    order[-1][i]=-1
    order[i][0]=-1
    order[i][-1]=-1       
nucx=[s+1]                           #路径位置记录
nucy=[s+1]
new=(s+1,s+1)                        #暂时位置
order[s+1][s+1]=1                    #设定起点

经简单推导可知，仅路径前端的前、左、右，斜前方5格内出现1才会出现锁死的可能，故就此进行分类讨论。（示意图如下，红色为路径，黄色为讨论区块）

垂直方向（五角星处）优先判断，出现一个1记一个“wall”（每步清零），且在可行方向集合中除去该方向。

wall=0
avi=[0,1,2,3]
ju=0
if k!=0:                         #第一步完全随机
    avi.remove((direc[-1]+2)%4)
    for i in np.arange(3):
        ju=order[new[0]+ev[(direc[-1]-1+i)%4][0]][new[1]+ev[(direc[-1]-1+i)%4][1]]
        if ju!=0:
            wall+=1
            avi.remove((direc[-1]-1+i)%4)

wall==3即为锁死，可输出锁死时的步数和位置，输出图进行比对，检查出现锁死的原因。（在前期用于修正和添加限制条件，认图找端点困难时步数也一般很大，出现锁死的几率不大，算法接近完整后可换成归零条件，即重新计算至能跑到规定步数为止）

 if wall==3:
     #print(k)
     #print(new)
     restart=True
     break

wall==2即只有一个方向可走，下一步方向唯一，无需再做排除。

wall==1即有一个方向被堵，此时分为两种情况：前方被堵和侧面被堵。其中前方被堵必定出现了路径闭合，只能选择一个方向。经推导易知在闭合圈不包围原点时，应排除与"wall"所在格前一步或后一步相同的方向（在可行方向内）；包围原点时，应排除相反的方向。

此时发现原本的数值矩阵不能用于该判断，再引入一个方向集合和顺序矩阵，顺序矩阵可以代替数值矩阵。方向集合记录每一步的方向（ev集合的序号），顺序矩阵初始也是零矩阵，每步的落点赋值该步的顺序。判断某点前后的方向只需以该点顺序作为索引调取方向集合的值。

判断闭合圈是否包围原点则从发生闭合的端点开始（此处是从"wall"所在格开始），沿路径检查其经过象限。新增一个判断集合，若新元素和集合内已有的倒数第二个元素相同（如出现1 4 1），则该元素和倒数第一个元素都可弃置（对是否围住原点没有影响）。因需和倒数第二个数做比较，判断集合初始为[0,0]。若闭合圈包围原点，则需经过四个象限，即判断集合长度为6。函数如下：

def roll(circle,k,nucx,nucy):         #判断闭合是否包围原点
    judge=[0,0]
    for i in range(circle,k-1):
        flag=0
        if (nucx[i-1]==s+1 and nucx[i]>s+1 and nucy[i]>s+1) or (nucy[i-1]==s+1 and nucy[i]>s+1 and nucx[i]>s+1):
            judge.append(1)
            flag=1
        elif (nucx[i-1]==s+1 and nucx[i]<s+1 and nucy[i]>s+1) or (nucy[i-1]==s+1 and nucy[i]>s+1 and nucx[i]<s+1):
            judge.append(2)
            flag=1
        elif (nucx[i-1]==s+1 and nucx[i]<s+1 and nucy[i]<s+1) or (nucy[i-1]==s+1 and nucy[i]<s+1 and nucx[i]<s+1):
            judge.append(3)
            flag=1
        elif (nucx[i-1]==s+1 and nucx[i]>s+1 and nucy[i]<s+1) or (nucy[i-1]==s+1 and nucy[i]<s+1 and nucx[i]>s+1):
            judge.append(4)
            flag=1
        if flag==1:
            if judge[-1]==judge[-3]:
                judge.pop(-1)
                judge.pop(-2)
            elif judge[-1]==judge[-2]:
                judge.pop(-1)
    if len(judge)==6:
        return 0
    else:
        return 1

当wall==1是侧面被堵时，则需进一步判断两个可行方向所夹的斜前方格是否为>0。若不>0则两侧互通，都可行；若>0则下一步必将发生路径闭合。判断方法类似，但闭合起点变为该斜方格,当闭合圈不包围原点时排除闭合起点后一步的方向，包围时排除闭合起点前一步的反方向。

if wall==1:
    if order[new[0]+ev[direc[-1]][0]][new[1]+ev[direc[-1]][1]]>1:        #前方被堵
        circle=int(order[new[0]+ev[direc[-1]][0]][new[1]+ev[direc[-1]][1]])
        may1=direc[circle-2]
        may2=direc[circle-1]
        if roll(circle,k,nucx,nucy):
            if avi[0]==may1 or avi[1]==may1:
                may3=may1
            else:
                may3=may2
            avi.remove(may3)
        else:
            if avi[0]==may1 or avi[1]==may1:
                may3=(may1+2)%4
            else:
                may3=(may2+2)%4
            avi.remove(may3)
    elif order[new[0]+ev[(direc[-1]+1)%4][0]][new[1]+ev[(direc[-1]+1)%4][1]]>1:  #左侧被堵
        circle=int(order[new[0]+ev[(direc[-1]-1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]-1)%4][1]+ev[direc[-1]][1]])
        if circle>1:
            if roll(circle,k,nucx,nucy):
                avi.remove(direc[circle-1])
            else:
                avi.remove((direc[circle-2]+2)%4)
    elif order[new[0]+ev[(direc[-1]-1)%4][0]][new[1]+ev[(direc[-1]-1)%4][1]]>1:   #右侧被堵
        circle=int(order[new[0]+ev[(direc[-1]+1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]+1)%4][1]+ev[direc[-1]][1]])
        if circle>1:
            if roll(circle,k,nucx,nucy):
                avi.remove(direc[circle-1])
            else:
                avi.remove((direc[circle-2]+2)%4)

wall==0时，进行类似wall==1侧面被堵时的判断。需进一步判断两个可行方向所夹的斜前方格是否>0，若>0则判断闭合圈是否包围原点，不包围则排除闭合起点的前后两步方向，包围则排除相反方向。

elif wall==0:
    if order[new[0]+ev[(direc[-1]+1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]+1)%4][1]+ev[direc[-1]][1]]>1:
        circle=int(order[new[0]+ev[(direc[-1]+1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]+1)%4][1]+ev[direc[-1]][1]])
        if roll(circle,k,nucx,nucy):
            if direc[circle-1] in avi:
                avi.remove(direc[circle-1])
            if direc[circle-2] in avi:
                avi.remove(direc[circle-2])
        else:
            if (direc[circle-1]+2)%4 in avi:
                avi.remove((direc[circle-1]+2)%4)
            if (direc[circle-2]+2)%4 in avi:
                avi.remove((direc[circle-2]+2)%4)
    if order[new[0]+ev[(direc[-1]-1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]-1)%4][1]+ev[direc[-1]][1]]>1:
        circle=int(order[new[0]+ev[(direc[-1]-1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]-1)%4][1]+ev[direc[-1]][1]])
        if roll(circle,k,nucx,nucy):
            if direc[circle-1] in avi:
                avi.remove(direc[circle-1])
            if direc[circle-2] in avi:
                avi.remove(direc[circle-2])
        else:
            if (direc[circle-1]+2)%4 in avi:
                avi.remove((direc[circle-1]+2)%4)
            if (direc[circle-2]+2)%4 in avi:
                avi.remove((direc[circle-2]+2)%4)

最后在可行方向集合avi中随机抽取方向决定下一步方向，保证随机性。更新暂存位置集合new、路径位置集合nuc、方向集合direc，落点处顺序矩阵order的赋值。

ra=random.choice(avi)                          #抽取方向
new=(new[0]+ev[ra][0],new[1]+ev[ra][1])        #更新暂时位置
nucx.append(new[0])                            #更新路径位置
nucy.append(new[1])
order[new[0]][new[1]]=k+2                      #更新顺序矩阵
direc.append(ra)                               #更新方向集合

用matplotlib绘制路径（nuc)。

plt.figure(figsize=(60,60))
plt.plot(nucx,nucy)
plt.grid(True)

十万步效果图：

缺陷

其实不难发现原点周边和边界会有特殊的锁死情况，不能用以上方法排除。边界问题可以用扩大边界的方法回避，原点目前还没有想到解决办法。但碰到因原点锁死的概率并不大，且碰上时步数也一般较小，就先用重算的方式摆烂了（懒，改天有心情有空就来完善）。

目前参数下输出十万步不会超过5秒，边界扩大2倍后输出五十万步在20秒左右。大步数下边界（和内存）是主要限制，numpy.zeros默认浮点型，dtype为int型后可以进一步提高步数上限（主要是顺序矩阵）。若能解决边界锁死的问题可以提高空间利用率，大辐度提高上限，目前的想法是利用回退，但还没写，改天吧（但愿）

附完整代码（后期优化了一下但没完全优化，肯定还有不少简化空间，凑合看吧）：

import matplotlib.pyplot as plt # 导入图形库
import numpy as np
import random
%matplotlib inline
def roll(circle,k,nucx,nucy):         #判断闭合是否包围原点
    judge=[0,0]
    for i in range(circle,k-1):
        flag=0
        if (nucx[i-1]==s+1 and nucx[i]>s+1 and nucy[i]>s+1) or (nucy[i-1]==s+1 and nucy[i]>s+1 and nucx[i]>s+1):
            judge.append(1)
            flag=1
        elif (nucx[i-1]==s+1 and nucx[i]<s+1 and nucy[i]>s+1) or (nucy[i-1]==s+1 and nucy[i]>s+1 and nucx[i]<s+1):
            judge.append(2)
            flag=1
        elif (nucx[i-1]==s+1 and nucx[i]<s+1 and nucy[i]<s+1) or (nucy[i-1]==s+1 and nucy[i]<s+1 and nucx[i]<s+1):
            judge.append(3)
            flag=1
        elif (nucx[i-1]==s+1 and nucx[i]>s+1 and nucy[i]<s+1) or (nucy[i-1]==s+1 and nucy[i]<s+1 and nucx[i]>s+1):
            judge.append(4)
            flag=1
        if flag==1:
            if judge[-1]==judge[-3]:
                judge.pop(-1)
                judge.pop(-2)
            elif judge[-1]==judge[-2]:
                judge.pop(-1)
    if len(judge)==6:
        return 0
    else:
        return 1
    
    
    
s=5000           #边界
target=100000    #目标步数
restart=True     #是否重新循环
while restart:
    restart=False
    order=np.zeros((2*s+3,2*s+3),dtype=int)       #顺序矩阵
    direc=[]                            #方向矩阵
    for i in np.arange(2*s+3):          #设定边界为-1
        order[0][i]=-1
        order[-1][i]=-1
        order[i][0]=-1
        order[i][-1]=-1
    ev=[(1,0),(0,1),(-1,0),(0,-1)]       
    nucx=[s+1]                           #路径位置记录
    nucy=[s+1]
    new=(s+1,s+1)                        #暂时位置
    order[s+1][s+1]=1                    #设定起点
    for k in np.arange(target):
        wall=0                         
        avi=[0,1,2,3]                    #可行方向集合
        ju=0                             #判断wall
        if k!=0:                         #第一步完全随机
            avi.remove((direc[-1]+2)%4)
            for i in np.arange(3):
                ju=order[new[0]+ev[(direc[-1]-1+i)%4][0]][new[1]+ev[(direc[-1]-1+i)%4][1]]
                if ju!=0:
                    wall+=1
                    avi.remove((direc[-1]-1+i)%4)

            if wall==3:
                #print(k)
                #print(new)
                restart=True            #重新循环
                break
            if wall==1:
                if order[new[0]+ev[direc[-1]][0]][new[1]+ev[direc[-1]][1]]>1:        #前方被堵
                    circle=int(order[new[0]+ev[direc[-1]][0]][new[1]+ev[direc[-1]][1]])
                    may1=direc[circle-2]
                    may2=direc[circle-1]
                    if roll(circle,k,nucx,nucy):
                        if avi[0]==may1 or avi[1]==may1:
                            may3=may1
                        else:
                            may3=may2
                        avi.remove(may3)
                    else:
                        if avi[0]==may1 or avi[1]==may1:
                            may3=(may1+2)%4
                        else:
                            may3=(may2+2)%4
                        avi.remove(may3)
                elif order[new[0]+ev[(direc[-1]+1)%4][0]][new[1]+ev[(direc[-1]+1)%4][1]]>1:  #左侧被堵
                    circle=int(order[new[0]+ev[(direc[-1]-1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]-1)%4][1]+ev[direc[-1]][1]])
                    if circle>1:
                        if roll(circle,k,nucx,nucy):
                            avi.remove(direc[circle-1])
                        else:
                            avi.remove((direc[circle-2]+2)%4)
                elif order[new[0]+ev[(direc[-1]-1)%4][0]][new[1]+ev[(direc[-1]-1)%4][1]]>1:   #右侧被堵
                    circle=int(order[new[0]+ev[(direc[-1]+1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]+1)%4][1]+ev[direc[-1]][1]])
                    if circle>1:
                        if roll(circle,k,nucx,nucy):
                            avi.remove(direc[circle-1])
                        else:
                            avi.remove((direc[circle-2]+2)%4)
            elif wall==0:
                if order[new[0]+ev[(direc[-1]+1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]+1)%4][1]+ev[direc[-1]][1]]>1:
                    circle=int(order[new[0]+ev[(direc[-1]+1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]+1)%4][1]+ev[direc[-1]][1]])
                    if roll(circle,k,nucx,nucy):
                        if direc[circle-1] in avi:
                            avi.remove(direc[circle-1])
                        if direc[circle-2] in avi:
                            avi.remove(direc[circle-2])
                    else:
                        if (direc[circle-1]+2)%4 in avi:
                            avi.remove((direc[circle-1]+2)%4)
                        if (direc[circle-2]+2)%4 in avi:
                            avi.remove((direc[circle-2]+2)%4)
                if order[new[0]+ev[(direc[-1]-1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]-1)%4][1]+ev[direc[-1]][1]]>1:
                    circle=int(order[new[0]+ev[(direc[-1]-1)%4][0]+ev[direc[-1]][0]][new[1]+ev[(direc[-1]-1)%4][1]+ev[direc[-1]][1]])
                    if roll(circle,k,nucx,nucy):
                        if direc[circle-1] in avi:
                            avi.remove(direc[circle-1])
                        if direc[circle-2] in avi:
                            avi.remove(direc[circle-2])
                    else:
                        if (direc[circle-1]+2)%4 in avi:
                            avi.remove((direc[circle-1]+2)%4)
                        if (direc[circle-2]+2)%4 in avi:
                            avi.remove((direc[circle-2]+2)%4)

        ra=random.choice(avi)                          #抽取方向
        new=(new[0]+ev[ra][0],new[1]+ev[ra][1])        #更新暂时位置
        nucx.append(new[0])                            #更新路径位置
        nucy.append(new[1])
        order[new[0]][new[1]]=k+2                      #更新顺序矩阵
        direc.append(ra)                               #更新方向集合

plt.figure(figsize=(60,60))
plt.plot(nucx,nucy)
plt.grid(True)