Percolation/Mean cluster density: Difference between revisions
(New draft task and Python solution.) |
(+ D entry) |
||
Line 31: | Line 31: | ||
;See: |
;See: |
||
* [http://mathworld.wolfram.com/s-Cluster.html s-Cluster] on Wolfram mathworld. |
* [http://mathworld.wolfram.com/s-Cluster.html s-Cluster] on Wolfram mathworld. |
||
=={{header|D}}== |
|||
{{trans|python}} |
|||
<lang d>import std.stdio, std.algorithm, std.random, std.math, std.array, |
|||
std.range, std.ascii; |
|||
alias Cell = uint; |
|||
alias Grid = Cell[][]; |
|||
enum Cell notClustered = 1; // Filled cell, but not in in a cluster. |
|||
Grid initialize(Grid grid, in double prob) /*nothrow*/ { |
|||
foreach (row; grid) |
|||
foreach (ref cell; row) |
|||
cell = uniform(cast(Cell)0, cast(Cell)2); |
|||
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(Grid grid) pure nothrow { |
|||
immutable side = grid.length; |
|||
Cell clusterID = 1; |
|||
void walk(in size_t r, in size_t c) nothrow { |
|||
grid[r][c] = clusterID; |
|||
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) { |
|||
clusterID++; |
|||
walk(r, c); |
|||
} |
|||
return clusterID - 1; |
|||
} |
|||
double clusterDensity(Grid grid, in double prob) { |
|||
return grid.initialize(prob).countClusters / |
|||
cast(double)(grid.length ^^ 2); |
|||
} |
|||
void showDemo(in size_t side, in double prob) { |
|||
auto grid = new Grid(side, side); |
|||
grid.initialize(prob); |
|||
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; |
|||
showDemo(15, prob); |
|||
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)) |
|||
//.sum / nIters; |
|||
.reduce!q{a + b} / nIters; |
|||
writefln("n_iters=%3d, p=%4.2f, n=%5d, sim=%7.8f", |
|||
nIters, prob, side, density); |
|||
} |
|||
}</lang> |
|||
{{out}} |
|||
<pre>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</pre> |
|||
Increasing the maximum index i to 16 (the grid takes about 1 GB): |
|||
<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= 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</pre> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
Revision as of 14:22, 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 = uint; alias Grid = Cell[][]; enum Cell notClustered = 1; // Filled cell, but not in in a cluster.
Grid initialize(Grid grid, in double prob) /*nothrow*/ {
foreach (row; grid) foreach (ref cell; row) cell = uniform(cast(Cell)0, cast(Cell)2); 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(Grid grid) pure nothrow {
immutable side = grid.length; Cell clusterID = 1;
void walk(in size_t r, in size_t c) nothrow { grid[r][c] = clusterID;
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) { clusterID++; walk(r, c); } return clusterID - 1;
}
double clusterDensity(Grid grid, in double prob) {
return grid.initialize(prob).countClusters / cast(double)(grid.length ^^ 2);
}
void showDemo(in size_t side, in double prob) {
auto grid = new Grid(side, side); grid.initialize(prob); 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;
showDemo(15, prob); 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)) //.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 (the grid takes about 1 GB):
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...