Percolation/Mean cluster density

Let $c$ be a 2D boolean square matrix of $n\times 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 (i.e., using the Von Neumann neighborhood rule) 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 $C_{n}$ .

Percolation theory states that $K(p)=C_{n}/n^{2}$ 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 $C_{n}$ for finite $n$ is subject to randomnes so should be computed as the average of $t$ runs, where $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 $15\times 15$ , $p=0.5$ grid.

Show your output here.

See also

s-Cluster on Wolfram mathworld.

D

Translation of: python

<lang d>import std.stdio, std.algorithm, std.random, std.math, std.array,

      std.range, std.ascii;

alias Cell = ubyte; alias Grid = Cell[][]; enum Cell notClustered = 1; // Filled cell, but not in a cluster.

FP rand(FP = double, UniformRandomNumberGenerator) //

           (ref UniformRandomNumberGenerator urng) {
   immutable FP result = urng.front / cast(FP)urng.max;
   urng.popFront;
   return result;

}

Grid initialize(Grid grid, in double prob, ref Xorshift rng) /*nothrow*/ {

   foreach (row; grid)
       foreach (ref cell; row)
           cell = cast(Cell)(rng.rand < prob);
   return grid;

}

void show(in Grid grid) {

   immutable static cell2char = " #" ~ letters;
   writeln('+', std.array.replicate("-", grid.length), '+');
   foreach (row; grid) {
       write('|');
       row.map!(c => c < cell2char.length ? cell2char[c] : '@').write;
       writeln('|');
   }
   writeln('+', std.array.replicate("-", grid.length), '+');

}

size_t countClusters(bool justCount=false)(Grid grid) pure nothrow {

   immutable side = grid.length;
   static if (justCount)
       enum Cell clusterID = 2;
   else
       Cell clusterID = 1;

   void walk(in size_t r, in size_t c) nothrow {
       grid[r][c] = clusterID; // Fill grid.

       if (r < side - 1 && grid[r + 1][c] == notClustered) // Down.
           walk(r + 1, c);
       if (c < side - 1 && grid[r][c + 1] == notClustered) // Right.
           walk(r, c + 1);
       if (c > 0 && grid[r][c - 1] == notClustered) // Left.
           walk(r, c - 1);
       if (r > 0 && grid[r - 1][c] == notClustered) // Up.
           walk(r - 1, c);
   }

   size_t nClusters = 0;

   foreach (immutable r; 0 .. side)
       foreach (immutable c; 0 .. side)
           if (grid[r][c] == notClustered) {
               static if (!justCount)
                   clusterID++;
               nClusters++;
               walk(r, c);
           }
   return nClusters;

}

double clusterDensity(Grid grid, in double prob, ref Xorshift rng) {

   return grid.initialize(prob, rng).countClusters!true /
          cast(double)(grid.length ^^ 2);

}

void showDemo(in size_t side, in double prob, ref Xorshift rng) {

   auto grid = new Grid(side, side);
   grid.initialize(prob, rng);
   writefln("Found %d clusters in this %d by %d grid:\n",
            grid.countClusters, side, side);
   grid.show;

}

void main() {

   immutable prob = 0.5;
   immutable nIters = 5;
   auto rng = Xorshift(unpredictableSeed);

   showDemo(15, prob, rng);
   writeln;
   foreach (immutable i; iota(4, 14, 2)) {
       immutable side = 2 ^^ i;
       auto grid = new Grid(side, side);
       immutable density = nIters
                           .iota
                           .map!(_ => grid.clusterDensity(prob, rng))
                           //.sum / nIters;
                           .reduce!q{a + b} / nIters;
       writefln("n_iters=%3d, p=%4.2f, n=%5d, sim=%7.8f",
                nIters, prob, side, density);
   }

}</lang>

Output:

Found 26 clusters in this 15 by 15 grid:

+---------------+
| AA B    CCCC  |
|AA D E F CC  G |
|  DDD FF  CC  H|
| I D FF  J   K |
|  L  FF JJJJ   |
|L LLL      J  M|
|LLLLLL   JJJ MM|
|L LL L N  J   M|
|LL    O P J   M|
|LLL QQ R JJ  S |
|LL T  RR  J SSS|
| L   U  V JJ  S|
|  WW  XX JJ YY |
|    XXX   JJ YY|
|ZZ   XXX   JJ  |
+---------------+

n_iters=  5, p=0.50, n=   16, sim=0.09765625
n_iters=  5, p=0.50, n=   64, sim=0.07260742
n_iters=  5, p=0.50, n=  256, sim=0.06679993
n_iters=  5, p=0.50, n= 1024, sim=0.06609497
n_iters=  5, p=0.50, n= 4096, sim=0.06580237

Increasing the index i to 15:

n_iters=  5, p=0.50, n=32768, sim=0.06578374

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...