更改

添加2,918字节 、 2020年10月16日 (五) 17:40
创建页面,内容为“==简介== 关于渗流模型的简介我们可以看这个页面:渗流模型 ==构建最基本的渗流模型== 调用包并设定参数 <syntaxhighlight…”
==简介==

关于渗流模型的简介我们可以看这个页面:[[渗流模型]]

==构建最基本的渗流模型==

调用包并设定参数

<syntaxhighlight lang="python">
from pylab import *
from scipy.ndimage import measurements
import statsmodels.api as sm

L = 100
r = rand(L,L)
p = 0.4
z = r < p
</syntaxhighlight>

绘图

[[File:Percolation_figure_1.png|500px]]

<syntaxhighlight lang="python">
imshow(z, origin='lower', interpolation='nearest')
colorbar()
title("Matrix")
</syntaxhighlight>

==标记不同的cluster==

[[File:Percolation_figure_2.png|500px]]

<syntaxhighlight lang="python">
lw, num = measurements.label(z)
imshow(lw, origin='lower', interpolation='nearest')
colorbar()
title("Labeled clusters")
</syntaxhighlight>

可以随机洗牌cluster label使得颜色更清楚

[[File:Percolation_figure_3.png|500px]]

<syntaxhighlight lang="python">
b = arange(lw.max() + 1)
shuffle(b)
shuffledLw = b[lw]
imshow(shuffledLw, origin='lower', interpolation='nearest')
colorbar()
title("Labeled clusters")
</syntaxhighlight>

还可以令cluster的颜色与大小正比

[[File:Percolation_figure_4.png|500px]]

<syntaxhighlight lang="python">
area = measurements.sum(z, lw, index=arange(lw.max() + 1))
areaImg = area[lw]
im3 = imshow(areaImg, origin='lower', interpolation='nearest')
colorbar()
title("Clusters by area")
</syntaxhighlight>

==给最大cluster加一个方框==

[[File:Percolation_figure_5.png|500px]]

<syntaxhighlight lang="python">
im3 = imshow(areaImg, origin='lower', interpolation='nearest')
colorbar()
title("Clusters by area")
sliced = measurements.find_objects(areaImg == areaImg.max())
if(len(sliced) > 0):
sliceX = sliced[0][1]
sliceY = sliced[0][0]
plotxlim=im3.axes.get_xlim()
plotylim=im3.axes.get_ylim()
plot([sliceX.start, sliceX.start, sliceX.stop, sliceX.stop, sliceX.start], \
[sliceY.start, sliceY.stop, sliceY.stop, sliceY.start, sliceY.start], \
color="red")
xlim(plotxlim)
ylim(plotylim)
</syntaxhighlight>

==拟合cluster size的分布==

[[File:Percolation_figure_6.png|500px]]

<syntaxhighlight lang="python">
def RankOrderPlot(data):
d=array(data)
d = d[d>0]
t=array(sorted(d,key=lambda x:-x))
r=array(range(1,len(d)+1))
y = log(t)
x = log(r)
X = sm.add_constant(x, prepend=True)
res = sm.OLS(y,X).fit()
C,beta = res.params
plt.plot(r,t,"o",color="b")
plt.plot(r,exp(C)*r**beta,"r-")
plt.xscale('log')
plt.yscale('log')
plt.text(max(r)/2,max(t)/100,"beta = " + str(round(beta,2)))
print beta

RankOrderPlot(area)
</syntaxhighlight>

[[category:python]]
[[category:旧词条迁移]]