#! /usr/bin/env python3
#
def percolation_simulation ( m, n, p ):
#*****************************************************************************80
#
## percolation_simulation() simulates and analyzes a percolation system.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 15 July 2022
#
# Author:
#
# Original MATLAB version by Ian Cooper
# Python version by John Burkardt
#
# Input:
#
# integer m, n: the number of rows and columns in the grid.
#
# real p: the probability that a given site should be occupied.
#
from numpy.random import default_rng
import matplotlib.pyplot as plt
import numpy as np
#
# The grid variable U is 0 or 1.
#
rng = default_rng()
R = rng.random ( [ m, n ] )
U = np.zeros ( [ m, n ], dtype = np.int )
U [ R < p ] = 1
#
# Count the number of occupied sites.
#
nosites = np.sum ( U )
#
# Observed probability of a site being occupied.
#
posites = nosites / m / n
#
# cls = copy of u with nonzeros replaced by cluster labels.
#
CLS = components_2d ( U )
#
# Number of clusters is maximum label in CLS.
#
cluster_num = np.max ( CLS )
#
# Number of occupied sites in each cluster.
#
cluster_size = np.zeros ( cluster_num, dtype = np.int )
for c in range ( 0, cluster_num ):
cluster_size[c] = np.sum ( U [ CLS == c ] )
#
# Counter from 1 to max number of occupied sites in a cluster.
#
cluster_size_max = np.max ( cluster_size )
#
# ns = cluster size distribution.
#
ns = np.zeros ( cluster_size_max + 1, dtype = np.int )
for c in range ( 0, cluster_size_max + 1 ):
ns[c] = len ( cluster_size [ cluster_size == c ] )
#
# Probability that a site chosen at random is part of
# an s-site cluster.
#
probros = ns / np.sum ( ns )
#
# Check for spanning clusters.
# Does a "c" occur in column 1 and column n?
# Does a "c" occur in row 1 and row m?
#
spanx = np.zeros ( cluster_num, dtype = np.int )
spany = np.zeros ( cluster_num, dtype = np.int )
for c in range ( 1, cluster_num ):
left = False
for i in range ( 0, m ):
if ( CLS[i,0] == c ):
left = True
break
right = False
for i in range ( 0, m ):
if ( CLS[i,-1] == c ):
right = True
break
if ( left and right ):
spanx[c] = 1
top = False
for j in range ( 0, n ):
if ( CLS[0,j] == c ):
top = True
break
bottom = False
for j in range ( 0, n ):
if ( CLS[-1,j] == c ):
bottom = True
break
if ( top and bottom ):
spany[c] = 1
#
# Output
#
print ( '' )
print ( ' Total number of sites = ', m * n )
print ( ' Number of occupied sites = ', nosites )
print ( ' Observed probability of occupation = ', posites )
print ( ' Number of clusters = ', cluster_num )
print ( ' Mean cluster size = ', np.mean ( cluster_size ) )
print ( '' )
print ( ' Cluster size distribution' )
print ( ' Size Number Probability' )
print ( '' )
for c in range ( 1, cluster_size_max + 1 ):
if ( 0 < ns[c] ):
print ( ' %2d %2d %g' % ( c, ns[c], probros[c] ) )
spanx_num = np.sum ( spanx )
spany_num = np.sum ( spany )
print ( '' )
print ( ' Number of horizontal spanning clusters = ', spanx_num )
print ( ' Number of vertical spanning clusters = ', spany_num )
#
# Plot the occupied cells.
#
plt.clf ( )
plt.axis ( 'equal' )
plt.axis ( 'off' )
outline = True
for i in range ( 0, m ):
for j in range ( 0, n ):
if ( U[i,j] == 0 ):
rgb = 'white'
else:
rgb = 'blue'
cell_ij_fill ( m, i, j, rgb, outline )
plt.title ( 'percolation occupation' )
filename = 'occupation.png'
plt.savefig ( filename )
print ( ' Graphics saved as "' + filename + '"' )
plt.show ( block = False )
plt.close ( )
#
# Second plot.
#
color = rng.random ( [ cluster_num + 1, 3 ] )
for i in range ( 0, cluster_num + 1 ):
color[i,:] = color[i,:] / np.linalg.norm ( color[i,:] )
plt.clf ( )
plt.axis ( 'equal' )
plt.axis ( 'off' )
outline = True
for i in range ( 0, m ):
for j in range ( 0, n ):
c = CLS[i,j]
if ( U[i,j] == 0 ):
rgb = 'white'
else:
rgb = color[c,:]
cell_ij_fill ( m, i, j, rgb, outline )
plt.title ( 'percolation clusters' )
filename = 'clusters.png'
plt.savefig ( filename )
print ( ' Graphics saved as "' + filename + '"' )
plt.show ( block = False )
plt.close ( )
return
def percolation_simulation_test ( ):
#*****************************************************************************80
#
## percolation_simulation_test() tests percolation_simulation().
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 15 July 2022
#
# Author:
#
# John Burkardt
#
import platform
print ( '' )
print ( 'percolation_simulation_test():' )
print ( ' Python version: ' + platform.python_version ( ) )
print ( ' percolation_simulation() simulates percolation.' )
print ( ' A rectangular region is formed by a grid of M x N squares.' )
print ( ' Each square may be porous or solid.' )
print ( ' We are interested in paths of porous squares connecting' )
print ( ' the top and bottom, or the left and right boundaries.' )
m = 32
n = 32
p = 0.60
print ( '' )
print ( ' Grid uses', m, 'rows x', n, 'columns.' )
print ( ' Prescribed probability of occupation =', p )
percolation_simulation ( m, n, p )
#
# Terminate.
#
print ( '' )
print ( 'percolation_simulation_test()' )
print ( ' Normal end of execution.' )
return
def components_2d ( A ):
#*****************************************************************************80
#
## components_2d() assigns contiguous nonzero pixels to a common component.
#
# Discussion:
#
# On input, the A array contains values of 0 or 1.
#
# The 0 pixels are to be ignored. The 1 pixels are to be grouped
# into connected components.
#
# The pixel A(I,J) is "connected" to the pixels A(I-1,J), A(I+1,J),
# A(I,J-1) and A(I,J+1), so most pixels have 4 neighbors.
#
# (Another choice would be to assume that a pixel was connected
# to the other 8 pixels in the 3x3 block containing it.)
#
# On output, COMPONENT_NUM reports the number of components of nonzero
# data, and the array C contains the component assignment for
# each nonzero pixel, and is 0 for zero pixels.
#
# Picture:
#
# Input A:
#
# 0 2 0 0 17 0 3
# 0 0 3 0 1 0 4
# 1 0 4 8 8 0 7
# 3 0 6 45 0 0 0
# 3 17 0 5 9 2 5
#
# Output:
#
# C:
#
# 0 1 0 0 2 0 3
# 0 0 2 0 2 0 3
# 4 0 2 2 2 0 3
# 4 0 2 2 0 0 0
# 4 4 0 2 2 2 2
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 16 July 2022
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer A(M,N), the pixel array.
#
# Output:
#
# integer C(M,N), the component array.
#
import numpy as np
m, n = A.shape
#
# Initialization.
#
C = np.zeros ( [ m, n ], dtype = np.int )
component_num = 0
#
# P is simply used to store the component labels. The dimension used
# here is, of course, usually an absurd overestimate.
#
p = np.arange ( 0, m * n + 1 )
#
# "Read" the array one pixel at a time. If a (nonzero) pixel's north or
# west neighbor already has a label, the current pixel inherits it.
# In case the labels disagree, we need to adjust the P array so we can
# later deal with the fact that the two labels need to be merged.
#
for i in range ( 0, m ):
for j in range ( 0, n ):
if ( i == 0 ):
north = 0
else:
north = C[i-1,j]
if ( j == 0 ):
west = 0
else:
west = C[i,j-1]
if ( A[i,j] != 0 ):
if ( north == 0 ):
if ( west == 0 ):
component_num = component_num + 1
C[i,j] = component_num
else:
C[i,j] = west
elif ( north != 0 ):
if ( west == 0 or west == north ):
C[i,j] = north
else:
C[i,j] = min ( north, west )
if ( north < west ):
p[west] = north
else:
p[north] = west
#
# When a component has multiple labels, have the higher labels
# point to the lowest one.
#
for component in range ( component_num -1, -1, -1 ):
b = component
while ( p[b] != b ):
b = p[b]
p[component] = b
#
# Locate the minimum label for each component.
# Assign these mininum labels new consecutive indices.
#
q = np.zeros ( component_num + 1, dtype = np.int )
i = 0
for component in range ( 0, component_num ):
if ( p[component] == component ):
i = i + 1
q[component] = i
component_num = i
#
# Replace the labels by consecutive labels.
# Need to avoid C(i,j) = 0.
#
i = np.where ( C != 0 )
C[i] = q [ p [ C[i] ] ] - 1
return C
def cell_ij_fill ( m, i, j, rgb, outline ):
#*****************************************************************************80
#
## cell_ij_fill() plots a filled (I,J) cell.
#
# Discussion:
#
# We assume the data is represented in a matrix.
#
# In order to convert between the matrix coordinates and picture
# coordinates, the (I,J) cell will be drawn with the following corners:
#
# (j-1,m-i+1), (j,m-i+1), (j,m-i), (j-1,m-1).
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 17 July 2022
#
# Author:
#
# John Burkardt
#
# Input:
#
# integer M, the maximum row index.
#
# integer I, J, the index of the cell.
#
# color RGB, can be any of the 8 abbreviated color terms
# 'r', 'g', 'b', 'c', 'm', 'y', 'w', 'k', or an RGB triple such as
# [1.0,0.4,0.0]. The square is filled with this color.
#
# bool outline: is True if the cell should be outlined in black.
#
import matplotlib.pyplot as plt
a = j - 1
b = j
c = m - ( i - 1 )
d = m - i
plt.fill ( [ a, b, b, a ], [ c, c, d, d ], color = rgb )
if ( outline ):
plt.plot ( [ a, b, b, a, a ], [ c, c, d, d, c ], 'k-' )
return
def timestamp ( ):
#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
# Licensing:
#
# This code is distributed under the MIT license.
#
# Modified:
#
# 21 August 2019
#
# Author:
#
# John Burkardt
#
import time
t = time.time ( )
print ( time.ctime ( t ) )
return
if ( __name__ == '__main__' ):
timestamp ( )
percolation_simulation_test ( )
timestamp ( )