Percolation/Site percolation: Difference between revisions
m →{{header|FORTRAN}}: Changed language name to "Fortran" |
Updated D entry |
||
Line 289:
foreach (ref row; grid)
foreach (ref cell; row) {
immutable r = rng.front /
rng.popFront;
cell = (r < probability) ? Cell.empty : Cell.filled;
Line 371:
writefln("\nFraction of %d tries that percolate through:", nTries);
foreach (const co; counters)
writefln("%1.3f %1.3f", co.prob, co.count /
writefln("\nSimulations and grid printing performed" ~
|
Revision as of 21:06, 16 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
Given an rectangular array of cells numbered assume is horizontal and is downwards.
Assume that the probability of any cell being filled is a constant where
- The task
Simulate creating the array of cells with probability and then testing if there is a route through adjacent filled cells from any on row to any on row , i.e. testing for site percolation.
Given repeat the percolation times to estimate the proportion of times that the fluid can percolate to the bottom for any given .
Show how the probability of percolating through the random grid changes with going from to in increments and with the number of repetitions to estimate the fraction at any given as .
Use an grid of cells for all cases.
Optionally depict a percolation through a cell grid graphically.
Show all output on this page.
C
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <string.h>
char *cell, *start, *end; int m, n;
void make_grid(int x, int y, double p) { int i, j, thresh = p * RAND_MAX;
m = x, n = y; end = start = realloc(start, (x+1) * (y+1) + 1);
memset(start, 0, m + 1);
cell = end = start + m + 1; for (i = 0; i < n; i++) { for (j = 0; j < m; j++) *end++ = rand() < thresh ? '+' : '.'; *end++ = '\n'; }
end[-1] = 0; end -= ++m; // end is the first cell of bottom row }
int ff(char *p) // flood fill { if (*p != '+') return 0;
*p = '#'; return p >= end || ff(p+m) || ff(p+1) || ff(p-1) || ff(p-m); }
int percolate(void) { int i; for (i = 0; i < m && !ff(cell + i); i++); return i < m; }
int main(void) { make_grid(15, 15, .5); percolate();
puts("15x15 grid:"); puts(cell);
puts("\nrunning 10,000 tests for each case:");
double p; int ip, i, cnt; for (ip = 0; ip <= 10; ip++) { p = ip / 10.; for (cnt = i = 0; i < 10000; i++) { make_grid(15, 15, p); cnt += percolate(); } printf("p=%.1f: %.4f\n", p, cnt / 10000.); }
return 0; }</lang>
- Output:
15x15 grid: .#...##.#.#.#.. ...+.###.####.# ...+..#.+...#.# +..+..##..##### +...+.#....##.. .+..+.##..##.+. ....+.#...##..+ ..+.+.#####.++. +++....#.###.++ .+.+.#.#.##.... ..++.####...++. +.+.+.##..+++.. +..+.+..+.....+ ..........++..+ .+.+.++++.+...+ running 10,000 tests for each case: p=0.0: 0.0000 p=0.1: 0.0000 p=0.2: 0.0000 p=0.3: 0.0000 p=0.4: 0.0032 p=0.5: 0.0902 p=0.6: 0.5771 p=0.7: 0.9587 p=0.8: 0.9996 p=0.9: 1.0000 p=1.0: 1.0000
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <time.h>
- include <string.h>
- include <stdbool.h>
- define N_COLS 15
- define N_ROWS 15
// Probability granularity 0.0, 0.1, ... 1.0
- define N_STEPS 11
// Simulation tries
- define N_TRIES 100
typedef unsigned char Cell; enum { EMPTY_CELL = ' ',
FILLED_CELL = '#', VISITED_CELL = '.' };
typedef Cell Grid[N_ROWS][N_COLS];
void initialize(Grid grid, const double probability) {
for (size_t r = 0; r < N_ROWS; r++) for (size_t c = 0; c < N_COLS; c++) { const double rnd = rand() / (double)RAND_MAX; grid[r][c] = (rnd < probability) ? EMPTY_CELL : FILLED_CELL; }
}
void show(Grid grid) {
char line[N_COLS + 3]; memset(&line[0], '-', N_COLS + 2); line[0] = '+'; line[N_COLS + 1] = '+'; line[N_COLS + 2] = '\0'; printf("%s\n", line); for (size_t r = 0; r < N_ROWS; r++) { putchar('|'); for (size_t c = 0; c < N_COLS; c++) putchar(grid[r][c]); puts("|"); } printf("%s\n", line);
}
bool walk(Grid grid, const size_t r, const size_t c) {
const size_t bottom = N_ROWS - 1; grid[r][c] = VISITED_CELL; if (r < bottom && grid[r + 1][c] == EMPTY_CELL) { // Down. if (walk(grid, r + 1, c)) return true; } else if (r == bottom) return true; if (c && grid[r][c - 1] == EMPTY_CELL) // Left. if (walk(grid, r, c - 1)) return true; if (c < N_COLS - 1 && grid[r][c + 1] == EMPTY_CELL) // Right. if (walk(grid, r, c + 1)) return true; if (r && grid[r - 1][c] == EMPTY_CELL) // Up. if (walk(grid, r - 1, c)) return true; return false;
}
bool percolate(Grid grid) {
const size_t startR = 0; for (size_t c = 0; c < N_COLS; c++) if (grid[startR][c] == EMPTY_CELL) if (walk(grid, startR, c)) return true; return false;
}
typedef struct {
double prob; size_t count;
} Counter;
int main() {
const double probability_step = 1.0 / (N_STEPS - 1); Counter counters[N_STEPS]; for (size_t i = 0; i < N_STEPS; i++) counters[i] = (Counter){ i * probability_step, 0 }; bool sample_shown = false; static Grid grid; srand(time(NULL)); for (size_t i = 0; i < N_STEPS; i++) { for (size_t t = 0; t < N_TRIES; t++) { initialize(grid, counters[i].prob); if (percolate(grid)) { counters[i].count++; if (!sample_shown) { printf("Percolating sample (%dx%d," " probability =%5.2f):\n", N_COLS, N_ROWS, counters[i].prob); show(grid); sample_shown = true; } } } } printf("\nFraction of %d tries that percolate through:\n", N_TRIES); for (size_t i = 0; i < N_STEPS; i++) printf("%1.1f %1.3f\n", counters[i].prob, counters[i].count / (double)N_TRIES); return 0;
} </lang>
- Output:
Percolating sample (15x15, probability = 0.40): +---------------+ |###. # # # #| |###.. # ##### | | #. ###### #| |###.... ######| |######. ### # | | #####.###### | |#......... ## | |...#...##.# ## | |##.#...##.### #| | ###..# #. # | |# #######. # ##| | # ##...#### | | ## # .##### | |#######.## ###| |# ## .## # # | +---------------+ Fraction of 100 tries that percolate through: 0.0 0.000 0.1 0.000 0.2 0.000 0.3 0.000 0.4 0.010 0.5 0.070 0.6 0.630 0.7 0.970 0.8 1.000 0.9 1.000 1.0 1.000
D
<lang d>import std.stdio, std.random, std.array, std.datetime;
enum size_t nCols = 15,
nRows = 15, nSteps = 11, // Probability granularity. nTries = 20_000; // Simulation tries.
alias BaseType = char; enum Cell : BaseType { empty = ' ',
filled = '#', visited = '.' }
alias Grid = Cell[nCols][nRows];
void initialize(ref Grid grid, in double probability,
ref Xorshift rng) { foreach (ref row; grid) foreach (ref cell; row) { immutable r = rng.front / double(rng.max); rng.popFront; cell = (r < probability) ? Cell.empty : Cell.filled; }
}
void show(in ref Grid grid) {
foreach (const ref row; grid) writeln('|', cast(BaseType[nCols])row, '|');
}
bool percolate(ref Grid grid) pure nothrow {
bool walk(in size_t r, in size_t c) nothrow { enum bottom = nRows - 1; grid[r][c] = Cell.visited;
if (r < bottom && grid[r + 1][c] == Cell.empty) { // Down. if (walk(r + 1, c)) return true; } else if (r == bottom) return true;
if (c && grid[r][c - 1] == Cell.empty) // Left. if (walk(r, c - 1)) return true;
if (c < nCols - 1 && grid[r][c + 1] == Cell.empty) // Right. if (walk(r, c + 1)) return true;
if (r && grid[r - 1][c] == Cell.empty) // Up. if (walk(r - 1, c)) return true;
return false; }
enum startR = 0; foreach (immutable c; 0 .. nCols) if (grid[startR][c] == Cell.empty) if (walk(startR, c)) return true; return false;
}
void main() {
static struct Counter { double prob; size_t count; }
StopWatch sw; sw.start;
enum probabilityStep = 1.0 / (nSteps - 1); Counter[nSteps] counters; foreach (immutable i, ref co; counters) co.prob = i * probabilityStep;
Grid grid; bool sampleShown = false; auto rng = Xorshift(unpredictableSeed);
foreach (ref co; counters) { foreach (immutable _; 0 .. nTries) { grid.initialize(co.prob, rng); if (grid.percolate) { co.count++; if (!sampleShown) { writefln("Percolating sample (%dx%d," ~ " probability =%5.2f):", nCols, nRows, co.prob); grid.show; sampleShown = true; } } } } sw.stop;
writefln("\nFraction of %d tries that percolate through:", nTries); foreach (const co; counters) writefln("%1.3f %1.3f", co.prob, co.count / double(nTries));
writefln("\nSimulations and grid printing performed" ~ " in %3.2f seconds.", sw.peek.msecs / 1000.0);
}</lang>
- Output:
Percolating sample (15x15, probability = 0.30): |#..##.#...# ##| |..#### ##... #| |### #####.#. #| | ### # ###. ##| |## ## ###.## | |# ### ####..##| |###### # ###.. | | ############.#| |### ##### #.. | |## # #####.##| | ######### #..#| |# ###### ###. | | # # # ###..#| |#### #### #..##| |######### #.# #| Fraction of 20000 tries that percolate through: 0.000 0.000 0.100 0.000 0.200 0.000 0.300 0.000 0.400 0.003 0.500 0.094 0.600 0.570 0.700 0.956 0.800 1.000 0.900 1.000 1.000 1.000 Simulations and grid printing performed in 1.13 seconds.
Fortran
Please see sample compilation and program execution in comments at top of program. Thank you. This example demonstrates recursion and integer constants of a specific kind. <lang fortran> ! loosely translated from python. ! compilation: gfortran -Wall -std=f2008 thisfile.f08
!$ a=site && gfortran -o $a -g -O0 -Wall -std=f2008 $a.f08 && $a !100 trials per !Fill Fraction goal(%) simulated through paths(%) ! 0 0 ! 10 0 ! 20 0 ! 30 0 ! 40 0 ! 50 6 ! ! ! b b b b h j m m m ! b b b b b h h m m m m m ! b b b h h h m ! b h h h h h h h ! b b h h h h h h h h h ! b b b h h h h h h h h h h ! b b @ h h h h h h h ! @ @ h h h h h h h h ! @ @ @ @ h h h h ! @ @ @ @ h h h h h h ! @ @ @ h h h h h h h ! @ @ @ h h h h h h ! @ h h h h h h ! @ h h h h h h h ! @ @ h h h h h h h h h h ! 60 59 ! 70 97 ! 80 100 ! 90 100 ! 100 100
program percolation_site
implicit none integer, parameter :: m=15,n=15,t=100 !integer, parameter :: m=2,n=2,t=8 integer(kind=1), dimension(m, n) :: grid real :: p integer :: i, ip, trial, successes logical :: success, unseen, q data unseen/.true./ write(6,'(i3,a11)') t,' trials per' write(6,'(a21,a30)') 'Fill Fraction goal(%)','simulated through paths(%)' do ip=0, 10 p = ip/10.0 successes = 0 do trial = 1, t call newgrid(grid, p) success = .false. do i=1, m q = walk(grid, i) ! deliberately compute all paths success = success .or. q end do if ((ip == 6) .and. unseen) then call display(grid) unseen = .false. end if successes = successes + merge(1, 0, success) end do write(6,'(9x,i3,24x,i3)')ip*10,nint(100*real(successes)/real(t)) end do
contains
logical function walk(grid, start) integer(kind=1), dimension(m,n), intent(inout) :: grid integer, intent(in) :: start walk = rwalk(grid, 1, start, int(start+1,1)) end function walk
recursive function rwalk(grid, i, j, k) result(through) logical :: through integer(kind=1), dimension(m,n), intent(inout) :: grid integer, intent(in) :: i, j integer(kind=1), intent(in) :: k logical, dimension(4) :: q !out of bounds through = .false. if (i < 1) return if (m < i) return if (j < 1) return if (n < j) return !visited or non-pore if (1_1 /= grid(i, j)) return !update grid and recurse with neighbors. deny 'shortcircuit' evaluation grid(i, j) = k q(1) = rwalk(grid,i+0,j+1,k) q(2) = rwalk(grid,i+0,j-1,k) q(3) = rwalk(grid,i+1,j+0,k) q(4) = rwalk(grid,i-1,j+0,k) !newly discovered outlet through = (i == m) .or. any(q) end function rwalk
subroutine newgrid(grid, probability) implicit none real :: probability integer(kind=1), dimension(m,n), intent(out) :: grid real, dimension(m,n) :: harvest call random_number(harvest) grid = merge(1_1, 0_1, harvest < probability) end subroutine newgrid
subroutine display(grid) integer(kind=1), dimension(m,n), intent(in) :: grid integer :: i, j, k, L character(len=n*2) :: lineout write(6,'(/)') lineout = ' ' do i=1,m do j=1,n k = j+j L = grid(i,j)+1 lineout(k:k) = ' @abcdefghijklmnopqrstuvwxyz'(L:L) end do write(6,*) lineout end do end subroutine display
end program percolation_site </lang>
J
<lang J> any =: +./ all =: *./
quickCheck =: [: all [: (any"1) 2 *./\ ] NB. a complete path requires connections between all row pairs
percolate =: 15 15&$: : (dyad define) NB. returns 0 iff blocked Use: (N, M) percolate P
NB. make a binary grid GRID =: y (> ?@($&0)) x
NB. compute the return value if. -. quickCheck GRID do. 0 return. end. STARTING_SITES =. 0 ,. ({. GRID) # i. {: x NB. indexes of 1 in head row of GRID any STARTING_SITES check GRID
)
NB. use local copy of GRID. Too slow.
check =: dyad define"1 2 NB. return 1 iff through path found use: START check GRID
GRID =. y LOCATION =. x if. 0 (= #) LOCATION do. 0 return. end. NB. no starting point? 0 if. LOCATION any@:((>: , 0 > [) $) GRID do. 0 return. end. NB. off grid? 0 INDEX =. <LOCATION if. 1 ~: INDEX { GRID do. 0 return. end. NB. fail. either already looked here or non-path if. (>: {. LOCATION) = (# GRID) do. 1 return. end. NB. Success! (display GRID here) G =: GRID =. INDEX (>:@:{)`[`]}GRID any GRID check~ LOCATION +"1 (, -)0 1,:1 0
)
NB. use global GRID. check =: dyad define"1 2 NB. return 1 iff through path found use: START check GRID
LOCATION =. x if. 0 (= #) LOCATION do. 0 return. end. NB. no starting point? 0 if. LOCATION any@:((>: , 0 > [) $) GRID do. 0 return. end. NB. off grid? 0 INDEX =. <LOCATION if. 1 ~: INDEX { GRID do. 0 return. end. NB. fail. either already looked here or non-path if. (>: {. LOCATION) = (# GRID) do. 1 return. end. NB. Success! (display GRID here) GRID =: INDEX (>:@:{)`[`]}GRID any GRID check~ LOCATION +"1 (, -)0 1,:1 0
)
simulate =: 100&$: : ([ %~ [: +/ [: percolate"0 #) NB. return fraction of connected cases. Use: T simulate P </lang>
,. ' P THRU' ; (, 100x&simulate)"0 (i. % <:)11x +-----------+ | P THRU | +-----------+ | 0 0| |1r10 0| | 1r5 0| |3r10 0| | 2r5 1r100| | 1r2 1r20| | 3r5 31r50| |7r10 97r100| | 4r5 1| |9r10 1| | 1 1| +-----------+ NB. example simulate 0.6 0.51 GRID NB. the final grid of the 100 simulated cases. 2 2 2 2 0 2 2 2 2 0 2 2 0 0 2 2 0 0 2 0 0 2 0 2 0 2 0 0 1 0 2 0 1 0 2 2 0 2 2 0 2 2 0 1 0 2 2 0 0 0 2 2 0 2 0 2 0 0 0 1 0 2 2 2 2 2 2 0 2 0 2 0 0 1 1 0 2 0 2 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 1 0 0 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 0 0 0 1 1 0 1 1 0 0 1 1 1 0 1 1 0 0 1 1 1 0 1 0 1 0 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 0 0 1 1 1 0 1 1 1 0 1 1 0 0 0 1 1 1 1 1 0 0 0 0 1 0 1 1 1 1 1 0 0 1 0 1 1 1 1 0 1 1 1 1 1 0 1 0 0 0 1 0 1 1 1 1 1 0 0 1 1 1 1 1 1 (0 ,. 0 6 10 14) check GRID NB. show possible starting points all fail 0 0 0 0 1j1#"1 GRID { '#',~u: 32 16bb7 NB. sample paths with unicode pepper. # # # # # # # # # # # # # # # # · # · # # # # # # · # # # # # # · # # # # # # # # · · # # · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·
Python
<lang python>from random import random import string from pprint import pprint as pp
M, N, t = 15, 15, 100
cell2char = ' #' + string.ascii_letters NOT_VISITED = 1 # filled cell not walked
class PercolatedException(Exception): pass
def newgrid(p):
return [[int(random() < p) for m in range(M)] for n in range(N)] # cell
def pgrid(cell, percolated=None):
for n in range(N): print( '%i) ' % (n % 10) + ' '.join(cell2char[cell[n][m]] for m in range(M))) if percolated: where = percolated.args[0][0] print('!) ' + ' ' * where + cell2char[cell[n][where]])
def check_from_top(cell):
n, walk_index = 0, 1 try: for m in range(M): if cell[n][m] == NOT_VISITED: walk_index += 1 walk_maze(m, n, cell, walk_index) except PercolatedException as ex: return ex return None
def walk_maze(m, n, cell, indx):
# fill cell cell[n][m] = indx # down if n < N - 1 and cell[n+1][m] == NOT_VISITED: walk_maze(m, n+1, cell, indx) # THE bottom elif n == N - 1: raise PercolatedException((m, indx)) # left if m and cell[n][m - 1] == NOT_VISITED: walk_maze(m-1, n, cell, indx) # right if m < M - 1 and cell[n][m + 1] == NOT_VISITED: walk_maze(m+1, n, cell, indx) # up if n and cell[n-1][m] == NOT_VISITED: walk_maze(m, n-1, cell, indx)
if __name__ == '__main__':
sample_printed = False pcount = {} for p10 in range(11): p = p10 / 10.0 pcount[p] = 0 for tries in range(t): cell = newgrid(p) percolated = check_from_top(cell) if percolated: pcount[p] += 1 if not sample_printed: print('\nSample percolating %i x %i, p = %5.2f grid\n' % (M, N, p)) pgrid(cell, percolated) sample_printed = True print('\n p: Fraction of %i tries that percolate through\n' % t ) pp({p:c/float(t) for p, c in pcount.items()})</lang>
- Output:
The Ascii art grid of cells has blanks for cells that were not filled. Filled cells start off as the '#', hash character and are changed to a succession of printable characters by successive tries to navigate from the top, (top - left actually), filled cell to the bottom.
The '!)' row shows where the percolation finished and you can follow the letter backwards from that row, (letter 'c' in this case), to get the route. The program stops after finding its first route through.
Sample percolating 15 x 15, p = 0.40 grid 0) a a a b c # 1) a a # c c # # 2) # # # # c # # # 3) # # # # # c 4) # # c c c c c c 5) # # # # # # c c c 6) # # # c c c 7) # # # # # # # c 8) # # # # # c c c 9) # # # c 0) # # # # # # c c # # 1) # # # # c 2) # # # # # # c c c c 3) # # # # c c c c 4) # # c # !) c p: Fraction of 100 tries that percolate through {0.0: 0.0, 0.1: 0.0, 0.2: 0.0, 0.3: 0.0, 0.4: 0.01, 0.5: 0.11, 0.6: 0.59, 0.7: 0.94, 0.8: 1.0, 0.9: 1.0, 1.0: 1.0}
Note the abrupt change in percolation at around p = 0.6. These abrupt changes are expected.
Racket
<lang racket>#lang racket (require racket/require (only-in racket/fixnum for*/fxvector)) (require (filtered-in (lambda (name) (regexp-replace #rx"unsafe-" name ""))
racket/unsafe/ops))
(define cell-empty 0) (define cell-filled 1) (define cell-wall 2) (define cell-visited 3) (define cell-exit 4)
(define ((percol->generator p)) (if (< (random) p) cell-filled cell-empty))
(define t (make-parameter 1000))
(define ((make-percol-grid M N) p)
(define p->10 (percol->generator p)) (define M+1 (fx+ 1 M)) (define M+2 (fx+ 2 M)) (for*/fxvector #:length (fx* N M+2) ((n (in-range N)) (m (in-range M+2))) (cond [(fx= 0 m) cell-wall] [(fx= m M+1) cell-wall] [else (p->10)])))
(define (cell->str c) (substring " #|+*" c (fx+ 1 c)))
(define ((draw-percol-grid M N) g)
(define M+2 (fx+ M 2)) (for ((row N)) (for ((col (in-range M+2))) (define idx (fx+ (fx* M+2 row) col)) (printf "~a" (cell->str (fxvector-ref g idx)))) (newline)))
(define ((percolate-percol-grid?! M N) g)
(define M+2 (fx+ M 2)) (define N-1 (fx- N 1)) (define max-idx (fx* N M+2)) (define (inner-percolate g idx) (define row (fxquotient idx M+2)) (cond ((fx< idx 0) #f) ((fx>= idx max-idx) #f) ((fx= N-1 row) (fxvector-set! g idx cell-exit) #t) ((fx= cell-filled (fxvector-ref g idx)) (fxvector-set! g idx cell-visited) (or ; gravity first (thanks Mr Newton) (inner-percolate g (fx+ idx M+2)) ; stick-to-the-left (inner-percolate g (fx- idx 1)) (inner-percolate g (fx+ idx 1)) ; go uphill only if we have to! (inner-percolate g (fx- idx M+2)))) (else #f))) (for/first ((m (in-range 1 M+2)) #:when (inner-percolate g m)) g))
(define make-15x15-grid (make-percol-grid 15 15)) (define draw-15x15-grid (draw-percol-grid 15 15)) (define perc-15x15-grid?! (percolate-percol-grid?! 15 15))
(define (display-sample-percolation p)
(printf "Percolation sample: p=~a~%" p) (for*/first ((i (in-naturals)) (g (in-value (make-15x15-grid 0.6))) #:when (perc-15x15-grid?! g)) (draw-15x15-grid g)) (newline))
(display-sample-percolation 0.4)
(for ((p (sequence-map (curry * 1/10) (in-range 0 (add1 10)))))
(define n-percolated-grids (for/sum ((i (in-range (t))) #:when (perc-15x15-grid?! (make-15x15-grid p))) 1)) (define proportion-percolated (/ n-percolated-grids (t))) (printf "p=~a\t->\t~a~%" p (real->decimal-string proportion-percolated 4)))</lang>
- Output:
Percolation sample: p=0.4 |+++ ++++ + +++| | +++ ++ # +| | + ++ ##++| | ## + ###+ | | ###### + #+++#| | ##### + + #| |## # # +++++## | |### # ++ +++# | |## ## +++++# | |# ### ++ + #| | ## ## +++ ##| |## ## +++ # #| |### # +### | |#### ####+ # | |# ## # *# #| p=0 -> 0.0000 p=1/10 -> 0.0000 p=1/5 -> 0.0000 p=3/10 -> 0.0000 p=2/5 -> 0.0030 p=1/2 -> 0.1110 p=3/5 -> 0.5830 p=7/10 -> 0.9530 p=4/5 -> 1.0000 p=9/10 -> 1.0000 p=1 -> 1.0000
Tcl
<lang tcl>package require Tcl 8.6
oo::class create SitePercolation {
variable cells w h constructor {width height probability} {
set w $width set h $height for {set cells {}} {[llength $cells] < $h} {lappend cells $row} { for {set row {}} {[llength $row] < $w} {lappend row $cell} { set cell [expr {rand() < $probability}] } }
} method print {out} {
array set map {0 "#" 1 " " -1 .} puts "+[string repeat . $w]+" foreach row $cells { set s "|" foreach cell $row { append s $map($cell) } puts [append s "|"] } set outline [lrepeat $w "-"] foreach index $out { lset outline $index "." } puts "+[join $outline {}]+"
} method percolate {} {
for {set work {}; set i 0} {$i < $w} {incr i} { if {[lindex $cells 0 $i]} {lappend work 0 $i} } try { my Fill $work return {} } trap PERCOLATED x { return [list $x] }
} method Fill {queue} {
while {[llength $queue]} { set queue [lassign $queue y x] if {$y >= $h} {throw PERCOLATED $x} if {$y < 0 || $x < 0 || $x >= $w} continue if {[lindex $cells $y $x]<1} continue lset cells $y $x -1 lappend queue [expr {$y+1}] $x [expr {$y-1}] $x lappend queue $y [expr {$x-1}] $y [expr {$x+1}] }
}
}
- Demonstrate one run
puts "Sample percolation, 15x15 p=0.6" SitePercolation create bp 15 15 0.6 bp print [bp percolate] bp destroy puts ""
- Collect statistics
apply {{} {
puts "Percentage of tries that percolate, varying p" set tries 100 for {set pint 0} {$pint <= 10} {incr pint} {
set p [expr {$pint * 0.1}] set tot 0 for {set i 0} {$i < $tries} {incr i} { set bp [SitePercolation new 15 15 $p] if {[$bp percolate] ne ""} { incr tot } $bp destroy } puts [format "p=%.2f: %2.1f%%" $p [expr {$tot*100./$tries}]]
}
}}</lang>
- Output:
Sample percolation, 15x15 p=0.6 +...............+ |.##...###.##...| |.#.#####.####..| |............##.| |....###.###.#..| |.#.##..#....#..| |#.........#..#.| |..#...##.##....| |#.#.#....##...#| |###.....#.#...#| |.....##........| |.#.#..## ......| | #..## # .##.#| | # #.# ####...| |# # # # ##...| | ### ## # . | +-------------.-+ Percentage of tries that percolate, varying p p=0.00: 0.0% p=0.10: 0.0% p=0.20: 0.0% p=0.30: 0.0% p=0.40: 0.0% p=0.50: 6.0% p=0.60: 54.0% p=0.70: 98.0% p=0.80: 100.0% p=0.90: 100.0% p=1.00: 100.0%