#! /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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 ( )