Percolation/Mean cluster density
Percolation Simulation
This is a simulation of aspects of mathematical percolation theory.
Mean run density
2D finite grid simulations
Site percolation | Bond percolation | Mean cluster density
Let be a 2D boolean square matrix of values of either 1 or 0 where the probability of any value being 1 is , (and of 0 is therefore ).
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 or by the limits of the matrix. Let the number of such clusters in such a randomly constructed matrix be .
Percolation theory states that as tends to infinity is a constant.
is called the mean cluster density and for is found numerically to approximate ...
- Task
Any calculation of for finite is subject to randomnes so should be
computed as the average of runs, where ≥ .
Show the effect of varying on the accuracy of simulated for and
for values of up to at least .
For extra credit, graphically show clusters in a , grid.
Show your output here.
- See also
- s-Cluster on Wolfram mathworld.
D
<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...