Percolation/Mean cluster density: Difference between revisions
(+ D entry) |
(A little more optimized D entry) |
||
Line 37: | Line 37: | ||
std.range, std.ascii; |
std.range, std.ascii; |
||
alias Cell = |
alias Cell = ubyte; |
||
alias Grid = Cell[][]; |
alias Grid = Cell[][]; |
||
enum Cell notClustered = 1; // Filled cell, but not in in a cluster. |
enum Cell notClustered = 1; // Filled cell, but not in in a cluster. |
||
FP rand(FP = double, UniformRandomNumberGenerator) // |
|||
⚫ | |||
(ref UniformRandomNumberGenerator urng) { |
|||
immutable FP result = urng.front / cast(FP)urng.max; |
|||
urng.popFront; |
|||
return result; |
|||
} |
|||
⚫ | |||
/*nothrow*/ { |
|||
foreach (row; grid) |
foreach (row; grid) |
||
foreach (ref cell; row) |
foreach (ref cell; row) |
||
cell = |
cell = cast(Cell)(rng.rand < prob); |
||
return grid; |
return grid; |
||
} |
} |
||
Line 59: | Line 67: | ||
} |
} |
||
size_t countClusters(Grid grid) pure nothrow { |
size_t countClusters(bool justCount=false)(Grid grid) pure nothrow { |
||
immutable side = grid.length; |
immutable side = grid.length; |
||
static if (justCount) |
|||
⚫ | |||
enum Cell clusterID = 2; |
|||
else |
|||
⚫ | |||
size_t nClusters = 0; |
|||
void walk(in size_t r, in size_t c) nothrow { |
void walk(in size_t r, in size_t c) nothrow { |
||
grid[r][c] = clusterID; |
grid[r][c] = clusterID; // Fill grid. |
||
if (r < side - 1 && grid[r + 1][c] == notClustered) // Down. |
if (r < side - 1 && grid[r + 1][c] == notClustered) // Down. |
||
Line 79: | Line 91: | ||
foreach (immutable c; 0 .. side) |
foreach (immutable c; 0 .. side) |
||
if (grid[r][c] == notClustered) { |
if (grid[r][c] == notClustered) { |
||
static if (!justCount) |
|||
clusterID++; |
|||
nClusters++; |
|||
walk(r, c); |
walk(r, c); |
||
} |
} |
||
return |
return nClusters; |
||
} |
} |
||
double clusterDensity(Grid grid, in double prob) { |
double clusterDensity(Grid grid, in double prob, ref Xorshift rng) { |
||
return grid.initialize(prob).countClusters / |
return grid.initialize(prob, rng).countClusters!true / |
||
cast(double)(grid.length ^^ 2); |
cast(double)(grid.length ^^ 2); |
||
} |
} |
||
void showDemo(in size_t side, in double prob) { |
void showDemo(in size_t side, in double prob, ref Xorshift rng) { |
||
auto grid = new Grid(side, side); |
auto grid = new Grid(side, side); |
||
grid.initialize(prob); |
grid.initialize(prob, rng); |
||
writefln("Found %d clusters in this %d by %d grid:\n", |
writefln("Found %d clusters in this %d by %d grid:\n", |
||
grid.countClusters, side, side); |
grid.countClusters, side, side); |
||
Line 101: | Line 115: | ||
immutable prob = 0.5; |
immutable prob = 0.5; |
||
immutable nIters = 5; |
immutable nIters = 5; |
||
auto rng = Xorshift(unpredictableSeed); |
|||
showDemo(15, prob); |
showDemo(15, prob, rng); |
||
writeln; |
writeln; |
||
foreach (immutable i; iota(4, 14, 2)) { |
foreach (immutable i; iota(4, 14, 2)) { |
||
Line 109: | Line 124: | ||
immutable density = nIters |
immutable density = nIters |
||
.iota |
.iota |
||
.map!(_ => grid.clusterDensity(prob)) |
.map!(_ => grid.clusterDensity(prob, rng)) |
||
//.sum / nIters; |
//.sum / nIters; |
||
.reduce!q{a + b} / nIters; |
.reduce!q{a + b} / nIters; |
||
Line 142: | Line 157: | ||
n_iters= 5, p=0.50, n= 1024, sim=0.06609497 |
n_iters= 5, p=0.50, n= 1024, sim=0.06609497 |
||
n_iters= 5, p=0.50, n= 4096, sim=0.06580237</pre> |
n_iters= 5, p=0.50, n= 4096, sim=0.06580237</pre> |
||
Increasing the maximum index i to 16 |
Increasing the maximum index i to 16: |
||
<pre>n_iters= 5, p=0.50, n= 16, sim=0.07343750 |
<pre>n_iters= 5, p=0.50, n= 16, sim=0.07343750 |
||
n_iters= 5, p=0.50, n= 64, sim=0.07290039 |
n_iters= 5, p=0.50, n= 64, sim=0.07290039 |
Revision as of 15:40, 14 August 2013
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
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.
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 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; size_t nClusters = 0;
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); }
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 maximum index i to 16:
n_iters= 5, p=0.50, n= 16, sim=0.07343750 n_iters= 5, p=0.50, n= 64, sim=0.07290039 n_iters= 5, p=0.50, n= 256, sim=0.06653748 n_iters= 5, p=0.50, n= 1024, sim=0.06601048 n_iters= 5, p=0.50, n= 4096, sim=0.06586764 n_iters= 5, p=0.50, n=16384, sim=0.06579554
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...