Percolation/Mean cluster density: Difference between revisions
(→Tcl: Added implementation) |
m (→{{header|Tcl}}: Some notes, less obfuscated code) |
||
Line 480: | Line 480: | ||
=={{header|Tcl}}== |
=={{header|Tcl}}== |
||
Note that the queue (variables <code>q</code> and <code>k</code>) used to remember where to find cells when flood-filling the cluster is maintained as a list ''segment''; the front of the list is not trimmed for performance reasons. (This would matter with very long queues, in which case the queue could be shortened occasionally; ''frequent'' trimming is still slower though, because Tcl backs its “list” datatype with arrays and not linked lists.) |
|||
{{works with|Tcl|8.6}} |
{{works with|Tcl|8.6}} |
||
<lang tcl>package require Tcl 8.6 |
<lang tcl>package require Tcl 8.6 |
||
Line 497: | Line 498: | ||
if {[lindex $grid $i $j] == -1} { |
if {[lindex $grid $i $j] == -1} { |
||
incr cl |
incr cl |
||
for {set q [list $i $j];set k |
for {set q [list $i $j];set k 0} {$k<[llength $q]} {incr k} { |
||
set y [lindex $q |
set y [lindex $q $k] |
||
set x [lindex $q [incr k]] |
set x [lindex $q [incr k]] |
||
if {[lindex $grid $y $x] != -1} continue |
if {[lindex $grid $y $x] != -1} continue |
Revision as of 19:08, 29 March 2014
You are encouraged to solve this task according to the task description, using any language you may know.
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 ). We 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 (the mean cluster density) will satisfy as tends to infinity. For , is found numerically to approximate ...
- Task
Show the effect of varying on the accuracy of simulated for and for values of up to at least . Any calculation of for finite is subject to randomness, so an approximation should be computed as the average of runs, where ≥ .
For extra credit, graphically show clusters in a , grid.
Show your output here.
- See also
- s-Cluster on Wolfram mathworld.
C
<lang c>#include <stdio.h>
- include <stdlib.h>
int *map, w, ww;
void make_map(double p) { int i, thresh = RAND_MAX * p; i = ww = w * w;
map = realloc(map, i * sizeof(int)); while (i--) map[i] = -(rand() < thresh); }
char alpha[] = "+.ABCDEFGHIJKLMNOPQRSTUVWXYZ" "abcdefghijklmnopqrstuvwxyz";
- define ALEN ((int)(sizeof(alpha) - 3))
void show_cluster(void) { int i, j, *s = map;
for (i = 0; i < w; i++) { for (j = 0; j < w; j++, s++) printf(" %c", *s < ALEN ? alpha[1 + *s] : '?'); putchar('\n'); } }
void recur(int x, int v) { if (x >= 0 && x < ww && map[x] == -1) { map[x] = v; recur(x - w, v); recur(x - 1, v); recur(x + 1, v); recur(x + w, v); } }
int count_clusters(void) { int i, cls;
for (cls = i = 0; i < ww; i++) { if (-1 != map[i]) continue; recur(i, ++cls); }
return cls; }
double tests(int n, double p) { int i; double k;
for (k = i = 0; i < n; i++) { make_map(p); k += (double)count_clusters() / ww; } return k / n; }
int main(void) { w = 15; make_map(.5); printf("width=15, p=0.5, %d clusters:\n", count_clusters()); show_cluster();
printf("\np=0.5, iter=5:\n"); for (w = 1<<2; w <= 1<<14; w<<=2) printf("%5d %9.6f\n", w, tests(5, .5));
free(map); return 0; }</lang>
- Output:
width=15, p=0.5, 23 clusters: A . . . B . C C C C . D . E . A . . B B . . . . . . . . . . A . . . . . F . . G . H . . I . . J J J . . K K . L . M . I . J J . . . K K K K . M M . . . . . . K K . K . K . M . N . O O . K K K . K . . . . N N N . O O . K K K K K . P . N . . Q . . K K . . . K . P . . . . . R . K K . . K K . P . . S . . . K K . . . . K . P . . . K K K K K K . . K K . . T . . K K . K . . . U . K . . T . . . K . K K K K . K K K . T . . . . . K . K . V . K K . . . W . p=0.5, iter=5: 4 0.125000 16 0.083594 64 0.064453 256 0.066864 1024 0.065922 4096 0.065836 16384 0.065774
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 / 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('+', "-".replicate(grid.length), '+'); foreach (row; grid) { write('|'); row.map!(c => c < cell2char.length ? cell2char[c] : '@').write; writeln('|'); } writeln('+', "-".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 / 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; 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...
Racket
<lang racket>#lang racket (require srfi/14) ; character sets
- much faster than safe fixnum functions
(require
racket/require ; for fancy require clause below (filtered-in (lambda (name) (regexp-replace #rx"unsafe-" name "")) racket/unsafe/ops) ; these aren't in racket/unsafe/ops (only-in racket/fixnum for/fxvector in-fxvector fxvector-copy))
- ...(but less safe). if in doubt use this rather than the one above
- (require racket/fixnum)
(define t (make-parameter 5))
(define (build-random-grid p M N)
(define p-num (numerator p)) (define p-den (denominator p)) (for/fxvector #:length (fx* M N) ((_ (in-range (* M N)))) (if (< (random p-den) p-num) 1 0)))
(define letters
(sort (char-set->list (char-set-intersection char-set:letter ; char-set:ascii )) char<?))
(define n-letters (length letters)) (define cell->char
(match-lambda (0 #\space) (1 #\.) (c (list-ref letters (modulo (- c 2) n-letters)))))
(define (draw-percol-grid M N . gs)
(for ((r N)) (for ((g gs)) (define row-str (list->string (for/list ((idx (in-range (* r M) (* (+ r 1) M)))) (cell->char (fxvector-ref g idx))))) (printf "|~a| " row-str)) (newline)))
(define (count-clusters! M N g)
(define (gather-cluster! k c) (when (fx= 1 (fxvector-ref g k)) (define k-r (fxquotient k M)) (define k-c (fxremainder k M)) (fxvector-set! g k c) (define-syntax-rule (gather-surrounds range? k+) (let ((idx k+)) (when (and range? (fx= 1 (fxvector-ref g idx))) (gather-cluster! idx c)))) (gather-surrounds (fx> k-r 0) (fx- k M)) (gather-surrounds (fx> k-c 0) (fx- k 1)) (gather-surrounds (fx< k-c (fx- M 1)) (fx+ k 1)) (gather-surrounds (fx< k-r (fx- N 1)) (fx+ k M)))) (define-values (rv _c) (for/fold ((rv 0) (c 2)) ((pos (in-range (fx* M N))) #:when (fx= 1 (fxvector-ref g pos))) (gather-cluster! pos c) (values (fx+ rv 1) (fx+ c 1)))) rv)
(define (display-sample-clustering p)
(printf "Percolation cluster sample: p=~a~%" p) (define g (build-random-grid p 15 15)) (define g+ (fxvector-copy g)) (define g-count (count-clusters! 15 15 g+)) (draw-percol-grid 15 15 g g+) (printf "~a clusters~%" g-count))
(define (experiment p n t)
(printf "Experiment: ~a ~a ~a\t" p n t) (flush-output) (define sum-Cn (for/sum ((run (in-range t))) (printf "[~a" run) (flush-output) (define g (build-random-grid p n n)) (printf "*") (flush-output) (define Cn (count-clusters! n n g)) (printf "]") (flush-output) Cn)) (printf "\tmean K(p) = ~a~%" (real->decimal-string (/ sum-Cn t (sqr n)) 6)))
(module+ main
(t 10) (for ((n (in-list '(4000 1000 750 500 400 300 200 100 15)))) (experiment 1/2 n (t))) (display-sample-clustering 1/2))
(module+ test
(define grd (build-random-grid 1/2 1000 1000)) (/ (for/sum ((g (in-fxvector grd)) #:when (zero? g)) 1) (fxvector-length grd)) (display-sample-clustering 1/2))</lang>
- Output:
Run from DrRacket, which runs the test and main modules. From the command line, you'll want two commands: ``racket percolation_m_c_d.rkt`` and ``raco test percolation_m_c_d.rkt`` for the same result.
Experiment: 1/2 4000 10 [0*][1*][2*][3*][4*][5*][6*][7*][8*][9*] mean K(p) = 0.065860 Experiment: 1/2 1000 10 [0*][1*][2*][3*][4*][5*][6*][7*][8*][9*] mean K(p) = 0.066130 Experiment: 1/2 750 10 [0*][1*][2*][3*][4*][5*][6*][7*][8*][9*] mean K(p) = 0.066195 Experiment: 1/2 500 10 [0*][1*][2*][3*][4*][5*][6*][7*][8*][9*] mean K(p) = 0.066522 Experiment: 1/2 400 10 [0*][1*][2*][3*][4*][5*][6*][7*][8*][9*] mean K(p) = 0.066778 Experiment: 1/2 300 10 [0*][1*][2*][3*][4*][5*][6*][7*][8*][9*] mean K(p) = 0.066813 Experiment: 1/2 200 10 [0*][1*][2*][3*][4*][5*][6*][7*][8*][9*] mean K(p) = 0.067908 Experiment: 1/2 100 10 [0*][1*][2*][3*][4*][5*][6*][7*][8*][9*] mean K(p) = 0.069980 Experiment: 1/2 15 10 [0*][1*][2*][3*][4*][5*][6*][7*][8*][9*] mean K(p) = 0.089778 Percolation cluster sample: p=1/2 |. ... . . | |A BBB A A | |... .. .... | |AAA AA AAAA | |. . .... ... | |A A AAAA AAA | |. . . .........| |A A C AAAAAAAAA| | ... .. ....| | AAA AA AAAA| |.. ......... ..| |AA AAAAAAAAA AA| | . ... | | A AAA | |. .. .. | |D AA AA | | .. ... . .. | | AA AAA E AA | |. .. .. . . | |F AA AA A A | |. ........ . ..| |F AAAAAAAA A AA| |.. . .... ... | |FF A AAAA AAA | | . . . .... | | F G A AAAA | |.... .. .. . .| |FFFF HH AA A A| | . .. .....| | F HH AAAAA| 8 clusters
Tcl
Note that the queue (variables q
and k
) used to remember where to find cells when flood-filling the cluster is maintained as a list segment; the front of the list is not trimmed for performance reasons. (This would matter with very long queues, in which case the queue could be shortened occasionally; frequent trimming is still slower though, because Tcl backs its “list” datatype with arrays and not linked lists.)
<lang tcl>package require Tcl 8.6
proc determineClusters {w h p} {
# Construct the grid set grid [lrepeat $h [lrepeat $w 0]] for {set i 0} {$i < $h} {incr i} {
for {set j 0} {$j < $w} {incr j} { lset grid $i $j [expr {rand() < $p ? -1 : 0}] }
} # Find (and count) the clusters set cl 0 for {set i 0} {$i < $h} {incr i} {
for {set j 0} {$j < $w} {incr j} { if {[lindex $grid $i $j] == -1} { incr cl for {set q [list $i $j];set k 0} {$k<[llength $q]} {incr k} { set y [lindex $q $k] set x [lindex $q [incr k]] if {[lindex $grid $y $x] != -1} continue lset grid $y $x $cl foreach dx {1 0 -1 0} dy {0 1 0 -1} { set nx [expr {$x+$dx}] set ny [expr {$y+$dy}] if { $nx >= 0 && $ny >= 0 && $nx < $w && $ny < $h && [lindex $grid $ny $nx] == -1 } then { lappend q $ny $nx } } } } }
} return [list $cl $grid]
}
- Print a sample 15x15 grid
lassign [determineClusters 15 15 0.5] n g puts "15x15 grid, p=0.5, with $n clusters" puts "+[string repeat - 15]+" foreach r $g {puts |[join [lmap x $r {format %c [expr {$x==0?32:64+$x}]}] ""]|} puts "+[string repeat - 15]+"
- Determine the densities as the grid size increases
puts "p=0.5, iter=5" foreach n {5 30 180 1080 6480} {
set tot 0 for {set i 0} {$i < 5} {incr i} {
lassign [determineClusters $n $n 0.5] nC incr tot $nC
} puts "n=$n, K(p)=[expr {$tot/5.0/$n**2}]"
}</lang>
- Output:
15x15 grid, p=0.5, with 21 clusters +---------------+ | A B CCCCC| | D A BBB C | |E B F CCCC| | B B F CC C| |BBB B BB CCC| |B BBBBBB CCCCC| | B B G C C| |H II G G J | |HH II G GG K| |HH II GGG GG K| | I G GGGG | |LL GGG GG M N| | L G G O P | |LLLL Q R | |L L S T UUU| +---------------+ p=0.5, iter=5 n=5, K(p)=0.184 n=30, K(p)=0.07155555555555557 n=180, K(p)=0.06880246913580246 n=1080, K(p)=0.0661267146776406 n=6480, K(p)=0.06582889898643499