Percolation/Mean cluster density

Mean cluster density

Let c be a 2D boolean square matrix of n by n values of either 1 or 0 where the probability of any value being 1 is p, (and of 0 is therefore 1-p).

Define a cluster of 1's as being a group of 1's connected vertically or horizontally and bounded by either 0 or by the limits of the matrix. Let the number of such clusters in such a randomly constructed matrix be Cn.

Percolation theory states that

   K(p) = Cn / n / n  as n tends to infinity is a constant.

K(p) is called the mean cluster density and for p = 0.5 is found numerically to approximate 0.065770...

Task

Any calculation of Cn for finite n is subject to randomnes so should be computed as the average of t runs, t >= 5.

Show the effect of varying n on the accuracy of simulated K(p) for p = 0.5 and for values of n up to at least 1000

For extra credit graphically show clusters in a 15x15 p=0.5 grid

Show your output here.

See

s-Cluster on Wolfram mathworld.

Python

<lang python>from __future__ import division from random import random import string from math import fsum

n_range, p, t = (2**n2 for n2 in range(4, 14, 2)), 0.5, 5 N = M = 15

NOT_CLUSTERED = 1 # filled but not clustered cell cell2char = ' #' + string.ascii_letters

def newgrid(n, p):

   return [[int(random() < p) for x in range(n)] for y in range(n)]

def pgrid(cell):

   for n in range(N):
       print( '%i)  ' % (n % 10) 
              + ' '.join(cell2char[cell[n][m]] for m in range(M)))

def cluster_density(n, p):

   cc = clustercount(newgrid(n, p))
   return cc / n / n

def clustercount(cell):

   walk_index = 1
   for n in range(N):
       for m in range(M):
           if cell[n][m] == NOT_CLUSTERED:
               walk_index += 1
               walk_maze(m, n, cell, walk_index)
   return walk_index - 1

def walk_maze(m, n, cell, indx):

   # fill cell 
   cell[n][m] = indx
   # down
   if n < N - 1 and cell[n+1][m] == NOT_CLUSTERED:
       walk_maze(m, n+1, cell, indx)
   # right
   if m < M - 1 and cell[n][m + 1] == NOT_CLUSTERED:
       walk_maze(m+1, n, cell, indx)
   # left
   if m and cell[n][m - 1] == NOT_CLUSTERED:
       walk_maze(m-1, n, cell, indx)
   # up
   if n and cell[n-1][m] == NOT_CLUSTERED:
       walk_maze(m, n-1, cell, indx)

if __name__ == '__main__':

   cell = newgrid(n=N, p=0.5)
   print('Found %i clusters in this %i by %i grid\n' 
         % (clustercount(cell), N, N))
   pgrid(cell)
   print()
   
   for n in n_range:
       N = M = n
       sim = fsum(cluster_density(n, p) for i in range(t)) / t
       print('t=%3i p=%4.2f n=%5i sim=%7.5f'
             % (t, p, n, sim))</lang>

Output:

Found 20 clusters in this 15 by 15 grid

0)  a a     b     c       d d d d
1)  a a   e     f   g g   d      
2)        e   f f f       d      
3)  h h   e     f   i i   d d    
4)        e       j     d d d d  
5)    k k   k   k   l         d d
6)  k k k k k   k   l   m   n    
7)  k k k   k   k   l     o   p p
8)      k   k k k   l l l   q    
9)  k     k k k k     l     q q q
0)  k   k k k         l     q q q
1)  k k     k k k       r r      
2)  k     k k       r r r   s s  
3)  k k k k   r r r r r     s s  
4)  k   k   t   r   r r     s s  

t=  5 p=0.50 n=   16 sim=0.08984
t=  5 p=0.50 n=   64 sim=0.07310
t=  5 p=0.50 n=  256 sim=0.06706
t=  5 p=0.50 n= 1024 sim=0.06612
t=  5 p=0.50 n= 4096 sim=0.06587

As n increases, the sim result gets closer to 0.065770...