datawhale12月学习——算法的应用：Percolation

1 题目

1.1 问题概述

详见AlgorithmRunning - Percolation

背景问题

Percolation: 渗透问题

• 给定一个由随机分布的绝缘材料和金属材料组成的复合系统：需要多大比例的材料是金属才能使复合系统成为电导体？
• 给定一个表面有水（或下面有油）的多孔景观，在什么条件下水能够流到底部（或油涌到地表）。

抽象问题为

使用 n×n 网格对渗透系统进行建模。每个格子要么是打开的，也叫做开放的格子，要么是被阻止的。如果一个方格，可以通过一系列相邻（左、右、上、下）开放格子连接到顶行中的开放格子。如果底行中有一个完整的格子，我们说系统会渗透。换句话说：从第一行到最后一行存在一条通路，这条通路上所有的格子都是打开的，这个模型我们叫做渗透。

而科学家关心的是是否存在一个概率值$p^{\ast }$，当格子开启的概率$p$小于$p^{\ast }$时，这个系统一定是非渗透的，否则一定是渗透的

这个问题可以通过蒙特卡洛模拟实验来实现，即进行T次实验，每次实验中

对一个初始化的 n×n 网格来说，每次随机选择一个未打开的格子并打开，直到网格系统联通，记录此时的打开格子数m
则该实验的$p^{\ast }=m/(n\times n)$

T次重复实验后可以得到 $p^{\ast }$的期望和95%置信区间。

1.2 实现要求描述

需要实现下面的两个接口：

• 如何判断系统是否是渗透的：Percolation

  • 初始化系统nxn，默认全部是0，表示关闭
  • 随机选择一个没打开的格子并打开它，1表示打开
  • 判断系统是否渗透
  • 输出渗透系统，并给出一个可能的渗透路径

• 设计模拟一个蒙特卡洛实验：PercolationStats

  • 初始模拟实验，对nxn的系统运行T次
  • 进行T次实验并记录每一次实验的实验结果
  • 计算均值、扰动值、置信区间的值
  • 输出T次实验的渗透阈值和置信区间

2 程序实现

2.1 Percolation

关键功能为：

1. 判断系统是否渗透

目标是判断一个系统是不是渗透，即是否存在第一层到最后一层的连通路径。每一个都有n个格子，这里面存在一个多对多的关系。需要依次判断第一层的每一个格子和最后一层的每一个格子是否能够渗透。其实我们并不关注是谁能进行渗透，只要存在即可，因此可以设置两个虚拟的点，start和end点，start和第一层每个打开的格子都是连通的，end和最后一层的每个打开的格子都是连通的，我们只要判断start和end是否渗透，即表示当前系统是不是渗透。

在本次实现中，由于判断在每次新开启一个格子时进行，因此，我们只需要对每一个新开启点，判断其是否同时与<start,end>联通即可。
判断过程如下：

• 对随机打开的(i,j)，判断周围联通的四个格子[(i-1,j), (i+1, j), (i, j-1), (i, j+1)]是否打开，并将其中打开的格子作为新的查找集
• 对查找集进行广度优先搜索，对查找集中的每一格子都进行上一步的判断，形成新的查找集，直到某一个格子与<start,end>其中一个联通
• 与<start,end>都联通，或查找集为空时，停止查找
• 在查找过程中需要去除查找集中已查找过的点

2. 输出渗透系统，给出从上至下的可能渗透路径

当我们的系统是渗透状态的时候，系统是由0-1组成的二维网格，此时存在一条可以从最上层向下渗透的系统化，我们的目的是可视化这条能渗透的路径，为了简单，我们可以将这条可以渗透的路径全部赋值成2，此时的系统是有0-1-2组成的，其中0表示未被打开，1表示被打开，但是不能被渗透到；2表示一个可以渗透路径上的格子。

注：第一层所有打开的格子都是可渗透的。

则只需要从第一层开始往下查找联通格子即可。若想得到上述矩阵，只需要将联通矩阵与打开矩阵相加。

则实现代码如下：

class Percolation(object):"""PercolationArgs:n(int): 构造矩阵的大小，构造成 (n, n)；0表示关闭，1表示打开"""def __init__(self, n):self.n = nself.data = np.zeros((n,n))self.max_try = n*n# 判断系统中放个(row, col)是否是打开状态def is_open(self, row, col):return self.data[row,col] == 1# 打开(row, col)这个格子，状态从0变成1  def open(self, row, col):self.data[row,col] = 1self.lastpt = (row,col)# 返回打开的格子数def number_of_open_sites(self):return self.data.sum()# 返回当前系统的状态，是nxn的矩阵def get_current_status(self):return self.data# 可视化系统，将满足从上到下系统的格子，进行显示，# 你需要将满足从第一行就连通的格子的状态，从1变成2def show_percolates(self):show = np.zeros((self.n,self.n))# 第一层所有打开的格子可渗透show[0,:] = self.data[0,:]cu = np.where(show[0,:]==1)[0].tolist()lis = [(0,each) for each in cu]totlis = []totlis = self._sub_show_percolates(lis,totlis)for (i,j) in totlis:show[i,j] = 1return show + self.datadef _sub_show_percolates(self,lis,totlis):if len(lis) == 0:return totlisnlis = []for (i,j) in lis:newlis = [(max(i-1,0),j), (min(i+1,self.n-1), j), (i, max(j-1,0)), (i, min(j+1,self.n-1))]newlis = list(set(newlis))newlis = [each for each in newlis if self.data[each[0],each[1]]==1]nlis.extend(newlis)nlis = list(set(nlis)-set(totlis))totlis.extend(nlis)#print(nlis)totlis = self._sub_show_percolates(nlis,totlis)return totlis# 返回true or false，表示当前系统是否是渗透的# 以新添加点为起点进行检索def percolates(self):start = Falseend = Falsetry:lis = [self.lastpt]except:return Falsetotlis = lis.copy()start,end = self._sub_percolates(lis,start,end,totlis)return start and enddef _sub_percolates(self,lis,start,end,totlis):if len(lis) == 0:return start,endnlis = []for (i,j) in lis:newlis = [(max(i-1,0),j), (min(i+1,self.n-1), j), (i, max(j-1,0)), (i, min(j+1,self.n-1))]newlis = list(set(newlis))newlis = [each for each in newlis if self.data[each[0],each[1]]==1]for each in newlis:if each[0] == 0:start = Truenewlis.remove(each)if each[0] == self.n-1:end = Truenewlis.remove(each)nlis.extend(newlis)if start and end:return start,endnlis = list(set(nlis)-set(totlis))totlis.extend(nlis)start,end = self._sub_percolates(nlis,start,end,totlis)return start,end# 运行模拟实验，每次打开一个格子，直到系统联通; 返回打开的格子个数def run(self):i = 0while i < self.max_try:[row,col] = np.random.randint(self.n,size = 2)if not self.is_open(row,col):self.open(row,col)i += 1if self.percolates():breakreturn i/self.max_try

2.2 PercolationStats

这一部分没有什么特别的，只需要反复实例化Percolation，并执行蒙特卡洛实验（Percolation.run()），然后记录实验结果。最后进行统计计算即可。

代码实现如下：

class PercolationStats(object):def __init__(self, n, t):self.n = n self.t = tself.record = np.zeros(t)# sample mean of percolation threshold# 渗透系统的阈值def mean(self):self.meanV = self.record.mean()return self.meanV# sample standard deviation of percolation threshold# T次实验渗透系统阈值对应的标准差def stddev(self):self.std = self.record.std()return self.std# low endpoint of 95% confidence interval# 95置信区间的下届def confidenceLow(self):self.lowcon = self.meanV - 1.96*self.std/(self.t**0.5)return self.lowcon# high endpoint of 95% confidence interval# 95置信区间的上届def confidenceHigh(self):self.highcon = self.meanV + 1.96*self.std/(self.t**0.5)return self.highcon# 系统会默认调用这个函数进行评测，这个函数必须实现# 进行t次模拟实验，需要返回5元组# <mean(渗透阈值), std(方差), low_conf(置信区间下界), high_conf(置信区间上界), precolation_status(T次实验中随机一个可视化的状态，需要将能够从上到下渗透的格子从1标记成2)>def run(self):for i in range(self.t):rand = Percolation(self.n)self.record[i] = rand.run()return (self.mean(), self.stddev(), self.confidenceLow(), self.confidenceHigh(), rand)

3 测试结果

3.1 本地测试

N = 15，T=20
则阈值、方差、置信下界、置信上界分别为

0.5866666666666667,
0.058626248703435205,
0.5609725803027954,
0.6123607530305379,

某次实验的可视化结果如下：

附可视化代码

fig,ax = plt.subplots(figsize=(15,15),dpi = 100)
sns.heatmap(t, cmap=sns.diverging_palette(20, 220, n=200), linewidths = 3,annot=True,cbar = False)
plt.show()

3.2 线上测试

线上测试地址

提交名字为：solo一下_tong

测试效果如图