Percolation/Bond percolation
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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle M \times N} rectangular array of cells numbered Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathrm{cell}[0..M-1, 0..N-1]} assume Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle M} is horizontal and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N} is downwards. Each Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathrm{cell}[m, n]} is bounded by (horizontal) walls Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathrm{hwall}[m, n]} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathrm{hwall}[m+1, n]} ; (vertical) walls Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathrm{vwall}[m, n]} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathrm{vwall}[m, n+1]}
Assume that the probability of any wall being present is a constant Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle p} where
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0.0 \le p \le 1.0}
Except for the outer horizontal walls at Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m = 0} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m = M} which are always present.
- The task
Simulate pouring a fluid onto the top surface (Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle n = 0} ) where the fluid will enter any empty cell it is adjacent to if there is no wall between where it currently is and the cell on the other side of the (missing) wall.
The fluid does not move beyond the horizontal constraints of the grid.
The fluid may move “up” within the confines of the grid of cells. If the fluid reaches a bottom cell that has a missing bottom wall then the fluid can be said to 'drip' out the bottom at that point.
Given Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle p} repeat the percolation Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle t} times to estimate the proportion of times that the fluid can percolate to the bottom for any given Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle p} .
Show how the probability of percolating through the random grid changes with Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle p} going from Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0.0} to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 1.0} in Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 0.1} increments and with the number of repetitions to estimate the fraction at any given Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle p} as Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle t = 100} .
Use an Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle M=10, N=10} grid of cells for all cases.
Optionally depict fluid successfully percolating through a grid graphically.
Show all output on this page.
D
Compared to the Python entry, the initialize
function is a performance optimization usesul when nTries is high.
<lang d>import std.stdio, std.random, std.array, std.algorithm, std.range,
std.typecons;
struct Grid {
immutable int nr, nc; bool[][] hWalls, vWalls, cells; alias MaybeCol = Nullable!(int, int.min); Xorshift rng;
this(in uint nRows, in uint nCols, in uint seed) /*pure*/ nothrow { nr = nRows; nc = nCols; hWalls = new typeof(hWalls)(nr + 1, nc); vWalls = new typeof(vWalls)(nr, nc + 1); cells = new typeof(cells)(nr, nc); rng.seed = seed; }
void initialize(in double probability) { foreach (ref row; hWalls) foreach (ref x; row) x = uniform(0.0, 1.0, rng) < probability; foreach (ref row; vWalls) foreach (immutable c, ref x; row) x = (c == 0 || c == nc) ? true : (uniform(0.0, 1.0, rng) < probability); foreach (ref row; cells) row[] = false; }
void show(in MaybeCol percolatedCol) const { // Horiz, vert, fill chars. static immutable h = ["+ ", "+-"], v = [" ", "|"], f = [" ", "#", "X"];
foreach (immutable r; 0 .. nr) { writefln("%-(%s%)+", nc.iota.map!(c => h[hWalls[r][c]])); writefln("%-(%s%)", iota(nc + 1) .map!(c => v[vWalls[r][c]] ~ f[c<nc ? cells[r][c] : 0])); } writefln("%-(%s%)+", nc.iota.map!(c => h[hWalls[nr][c]])); if (!percolatedCol.isNull) writeln(" ".replicate(percolatedCol), " ",f[2]); }
MaybeCol pourOnTop() pure nothrow { MaybeCol floodFill(in int r, in int c) pure nothrow { cells[r][c] = true; // Fill cell.
// Bottom. if (r < nr - 1 && !hWalls[r + 1][c] && !cells[r + 1][c]) return floodFill(r + 1, c); else if (r == nr - 1 && !hWalls[r + 1][c]) // THE bottom. return MaybeCol(c);
// Left. if (c && !vWalls[r][c] && !cells[r][c - 1]) return floodFill(r, c - 1);
// Right. if (c < nc - 1 && !vWalls[r][c + 1] && !cells[r][c + 1]) return floodFill(r, c + 1);
// Top. if (r && !hWalls[r][c] && !cells[r - 1][c]) return floodFill(r - 1, c);
return MaybeCol(); }
immutable r = 0; // From top. foreach (immutable c; 0 .. nc) if (!hWalls[r][c]) { immutable percolatedCol = floodFill(r, c); if (!percolatedCol.isNull) return percolatedCol; } return MaybeCol(); }
}
void main() {
enum int nr = 10, nc = 10; // N. rows and columns of the grid. enum int nTries = 1_000; // N. simulations for each probability. enum int nStepsProb = 10; // N. steps of probability.
bool sampleShown = false; alias Pair = Tuple!(double,"prob", uint,"co"); Pair[] results;
foreach (immutable i; 0 .. nStepsProb + 1) { // Count down so sample print is interesting. immutable probability = (nStepsProb - i) / cast(double)nStepsProb; auto grid = Grid(nr, nc, unpredictableSeed); uint lastCount = 0;
foreach (immutable _; 0 .. nTries) { grid.initialize(probability); immutable percolatedCol = grid.pourOnTop; if (!percolatedCol.isNull) { lastCount++; if (!sampleShown) { writefln("Percolating sample %dx%d grid," ~ " probability=%1.1f:", nc,nr,probability); grid.show(percolatedCol); sampleShown = true; } } } results ~= Pair(probability, lastCount); }
writefln("\nFraction of %d tries that percolate," ~ " varying probability p:", nTries); foreach (const r; results) writefln("p=%1.2f: %1.3f", r.prob, r.co / cast(double)nTries);
}</lang>
- Output:
Percolating sample 10x10 grid, probability=0.6: + + + + + +-+-+-+-+-+ |#|# # # #| | | | + + + + + + + + +-+-+ |# #|# #|# #| | | | | + + + +-+-+ +-+ +-+ + |# #|#| | |#| | +-+ +-+ +-+ +-+ + +-+ | |#| # # # | | | +-+-+ + +-+-+-+ + + + | | | |# #| | | | | +-+ + +-+ +-+ + +-+-+ | | |# #| | | | | +-+-+-+ +-+ + +-+ + + | | | |# | | | +-+-+ + +-+-+-+-+-+-+ | |# | | | | + +-+-+ + +-+ + + + + | #| | | | | | +-+ + + +-+ + +-+ + + | | #| | | | + + + + +-+-+-+ +-+ + X Fraction of 1000 tries that percolate, varying probability p: p=1.00: 0.000 p=0.90: 0.000 p=0.80: 0.000 p=0.70: 0.000 p=0.60: 0.007 p=0.50: 0.142 p=0.40: 0.573 p=0.30: 0.929 p=0.20: 0.998 p=0.10: 1.000 p=0.00: 1.000
Python
<lang python>from collections import namedtuple from random import random from pprint import pprint as pp
Grid = namedtuple('Grid', 'cell, hwall, vwall')
M, N, t = 10, 10, 100
class PercolatedException(Exception): pass
HVF = [(' .', ' _'), (':', '|'), (' ', '#')] # Horiz, vert, fill chars
def newgrid(p):
hwall = [[int(random() < p) for m in range(M)] for n in range(N+1)] vwall = [[(1 if m in (0, M) else int(random() < p)) for m in range(M+1)] for n in range(N)] cell = [[0 for m in range(M)] for n in range(N)] return Grid(cell, hwall, vwall)
def pgrid(grid, percolated=None):
cell, hwall, vwall = grid h, v, f = HVF for n in range(N): print(' ' + .join(h[hwall[n][m]] for m in range(M))) print('%i) ' % (n % 10) + .join(v[vwall[n][m]] + f[cell[n][m] if m < M else 0] for m in range(M+1))[:-1]) n = N print(' ' + .join(h[hwall[n][m]] for m in range(M))) if percolated: where = percolated.args[0][0] print('!) ' + ' ' * where + ' ' + f[1])
def pour_on_top(grid):
cell, hwall, vwall = grid n = 0 try: for m in range(M): if not hwall[n][m]: flood_fill(m, n, cell, hwall, vwall) except PercolatedException as ex: return ex return None
def flood_fill(m, n, cell, hwall, vwall):
# fill cell cell[n][m] = 1 # bottom if n < N - 1 and not hwall[n + 1][m] and not cell[n+1][m]: flood_fill(m, n+1, cell, hwall, vwall) # THE bottom elif n == N - 1 and not hwall[n + 1][m]: raise PercolatedException((m, n+1)) # left if m and not vwall[n][m] and not cell[n][m - 1]: flood_fill(m-1, n, cell, hwall, vwall) # right if m < M - 1 and not vwall[n][m + 1] and not cell[n][m + 1]: flood_fill(m+1, n, cell, hwall, vwall) # top if n and not hwall[n][m] and not cell[n-1][m]: flood_fill(m, n-1, cell, hwall, vwall)
if __name__ == '__main__':
sample_printed = False pcount = {} for p10 in range(11): p = (10 - p10) / 10.0 # count down so sample print is interesting pcount[p] = 0 for tries in range(t): grid = newgrid(p) percolated = pour_on_top(grid) if percolated: pcount[p] += 1 if not sample_printed: print('\nSample percolating %i x %i grid' % (M, N)) pgrid(grid, percolated) sample_printed = True print('\n p: Fraction of %i tries that percolate through' % t ) pp({p:c/float(t) for p, c in pcount.items()})</lang>
- Output:
In the Ascii art, cells are either a space or a hash and are surrounded by either '_', '|' for intact walls and '.' and ':' for missing (leaky) walls.
The bottom-most line starting '!)' shows where the fluid can drip out from. (The percolation stops when one route through the bottom is found).
Sample percolating 10 x 10 grid _ _ . _ . _ _ . _ _ 0) | |#:#:#|#| | :#| | | _ _ . _ _ _ . . _ _ 1) | | |#:#| | | |#| : | _ _ _ . _ . . . . _ 2) | | |#:#| : | |#: | | _ _ _ _ . . _ . . . 3) | : : | | | : |#: | | _ _ . _ . . _ . _ _ 4) | : : : | | | |#: : | _ _ _ . _ _ _ . . _ 5) | : | | : | | :#| | | _ _ . . _ _ _ . _ . 6) | : | | : | |#:#:#| | _ . _ _ . _ _ _ . . 7) | : | : | : | | |#: | _ _ _ . . _ _ . . _ 8) | | : | | | |#:#:#: | _ _ _ . . . . _ _ . 9) | : : | : : :#: | : | . _ . _ . . . . _ _ !) # p: Fraction of 100 tries that percolate through {0.0: 1.0, 0.1: 1.0, 0.2: 1.0, 0.3: 1.0, 0.4: 0.9, 0.5: 0.47, 0.6: 0.06, 0.7: 0.0, 0.8: 0.0, 0.9: 0.0, 1.0: 0.0}
Note the abrupt cut-off in percolation at around p = 0.5 which is to be expected.
Racket
<lang racket>#lang racket
(define has-left-wall? (lambda (x) (bitwise-bit-set? x 0))) (define has-right-wall? (lambda (x) (bitwise-bit-set? x 1))) (define has-top-wall? (lambda (x) (bitwise-bit-set? x 2))) (define has-bottom-wall? (lambda (x) (bitwise-bit-set? x 3))) (define has-fluid? (lambda (x) (bitwise-bit-set? x 4)))
(define (walls->cell l? r? t? b?)
(+ (if l? 1 0) (if r? 2 0) (if t? 4 0) (if b? 8 0)))
(define (bonded-percol-grid M N p)
(define rv (make-vector (* M N))) (for* ((idx (in-range (* M N)))) (define left-wall? (or (zero? (modulo idx M)) (has-right-wall? (vector-ref rv (sub1 idx))))) (define right-wall? (or (= (modulo idx M) (sub1 M)) (< (random) p))) (define top-wall? (if (< idx M) (< (random) p) (has-bottom-wall? (vector-ref rv (- idx M))))) (define bottom-wall? (< (random) p)) (define cell-value (walls->cell left-wall? right-wall? top-wall? bottom-wall?)) (vector-set! rv idx cell-value)) rv)
(define (display-percol-grid M . vs)
(define N (/ (vector-length (car vs)) M)) (define-syntax-rule (tab-eol m) (when (= m (sub1 M)) (printf "\t"))) (for ((n N)) (for* ((v vs) (m M)) (when (zero? m) (printf "+")) (printf (match (vector-ref v (+ (* n M) m)) ((? has-top-wall?) "-+") ((? has-fluid?) "#+") (else ".+"))) (tab-eol m)) (newline) (for* ((v vs) (m M)) (when (zero? m) (printf "|")) (printf (match (vector-ref v (+ (* n M) m)) ((and (? has-fluid?) (? has-right-wall?)) "#|") ((? has-right-wall?) ".|") ((? has-fluid?) "##") (else ".."))) (tab-eol m)) (newline)) (for* ((v vs) (m M)) (when (zero? m) (printf "+")) (printf (match (vector-ref v (+ (* (sub1 M) M) m)) ((? has-bottom-wall?) "-+") ((? has-fluid?) "#+") (else ".+"))) (tab-eol m)) (newline))
(define (find-bonded-grid-t/b-path M v)
(define N (/ (vector-length v) M)) (define (flood-cell idx) (cond [(= (quotient idx M) N) #t] ; wootiments! [(has-fluid? (vector-ref v idx)) #f] ; been here [else (define cell (vector-ref v idx)) (vector-set! v idx (bitwise-ior cell 16)) (or (and (not (has-bottom-wall? cell)) (flood-cell (+ idx M))) (and (not (has-left-wall? cell)) (flood-cell (- idx 1))) (and (not (has-right-wall? cell)) (flood-cell (+ idx 1))) (and (not (has-top-wall? cell)) (>= idx M) ; not top row (flood-cell (- idx M))))])) (for/first ((m (in-range M)) #:unless (has-top-wall? (vector-ref v m)) #:when (flood-cell m)) #t))
(define t (make-parameter 1000)) (define (experiment p)
(/ (for*/sum ((sample (in-range (t))) (v (in-value (bonded-percol-grid 10 10 p))) #:when (find-bonded-grid-t/b-path 10 v)) 1) (t)))
(define (main)
(for ((tenths (in-range 0 (add1 10)))) (define p (/ tenths 10)) (define e (experiment p)) (printf "proportion of grids that percolate p=~a : ~a (~a)~%" p e (real->decimal-string e 5))))
(module+ test
(define (make/display/flood/display-bonded-grid M N p attempts (atmpt 1)) (define v (bonded-percol-grid M N p)) (define v+ (vector-copy v)) (cond [(or (find-bonded-grid-t/b-path M v+) (= attempts 0)) (define v* (vector-copy v+)) (define (flood-bonded-grid) (when (find-bonded-grid-t/b-path M v*) (flood-bonded-grid))) (flood-bonded-grid) (display-percol-grid M v v+ v*) (printf "After ~a attempt(s)~%~%" atmpt)] [else (make/display/flood/display-bonded-grid M N p (sub1 attempts) (add1 atmpt))])) (make/display/flood/display-bonded-grid 10 10 0 20) (make/display/flood/display-bonded-grid 10 10 .25 20) (make/display/flood/display-bonded-grid 10 10 .50 20) (make/display/flood/display-bonded-grid 10 10 .75 20000))</lang>
- Output:
Welcome to DrRacket, version 5.3.5 [3m]. Language: racket [custom]; memory limit: 1024 MB. +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ |...................| |##.................| |###################| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ |...................| |##.................| |###################| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ |...................| |##.................| |###################| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ |...................| |##.................| |###################| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ |...................| |##.................| |###################| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ |...................| |##.................| |###################| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ |...................| |##.................| |###################| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ |...................| |##.................| |###################| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ |...................| |##.................| |###################| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ |...................| |##.................| |###################| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+#+#+#+#+#+#+#+ After 1 attempt(s) +.+-+-+.+.+.+-+.+.+.+ +#+-+-+.+.+.+-+.+.+.+ +#+-+-+#+#+#+-+#+#+#+ |...................| |##.................| |##..###############| +.+-+.+-+.+.+-+-+-+.+ +#+-+.+-+.+.+-+-+-+.+ +#+-+#+-+#+#+-+-+-+#+ |.................|.| |##...............|.| |##..##..####.....|#| +.+-+.+.+.+.+-+.+.+.+ +#+-+.+.+.+.+-+.+.+.+ +#+-+#+.+#+#+-+.+.+#+ |.............|.....| |##...........|.....| |######..#####|....#| +.+.+.+.+.+.+.+.+.+.+ +#+.+.+.+.+.+.+.+.+.+ +#+#+#+.+#+#+#+.+.+#+ |.....|...|.|.......| |##...|...|.|.......| |#####|..#|#|##....#| +.+.+.+.+.+.+.+-+-+.+ +#+.+.+.+.+.+.+-+-+.+ +#+#+#+#+#+#+#+-+-+#+ |.|.............|...| |#|.............|...| |#|############.|..#| +.+-+-+.+-+.+.+.+.+.+ +#+-+-+.+-+.+.+.+.+.+ +#+-+-+#+-+#+#+.+.+#+ |...................| |##.................| |##....##..####....#| +.+.+-+.+.+.+.+-+-+.+ +#+.+-+.+.+.+.+-+-+.+ +#+.+-+#+.+#+#+-+-+#+ |...|...|...........| |##.|...|...........| |##.|###|..####..###| +.+.+.+-+.+.+.+.+.+.+ +#+#+.+-+.+.+.+.+.+.+ +#+#+#+-+.+#+#+.+#+#+ |...|...|.........|.| |###|...|.........|.| |###|##.|..####..#|#| +-+.+.+-+-+.+.+.+.+-+ +-+#+.+-+-+.+.+.+.+-+ +-+#+#+-+-+#+#+.+#+-+ |.....|.........|...| |..##.|.........|...| |..###|....####.|###| +.+.+.+.+.+.+.+.+.+.+ +.+#+.+.+.+.+.+.+.+.+ +.+#+#+.+.+#+#+#+#+#+ |.........|.......|.| |..##.....|.......|.| |..####...|#######|#| +.+.+.+-+.+.+-+.+-+.+ +.+#+.+-+.+.+-+.+-+.+ +.+#+#+-+.+#+-+#+-+#+ After 1 attempt(s) +.+.+.+.+-+-+.+-+.+.+ +#+#+#+#+-+-+.+-+.+.+ +#+#+#+#+-+-+#+-+#+#+ |.........|.|.|...|.| |########.|.|.|...|.| |########.|.|#|###|#| +.+-+-+.+-+-+-+.+.+-+ +#+-+-+#+-+-+-+.+.+-+ +#+-+-+#+-+-+-+#+#+-+ |...|...|...|.|.|.|.| |###|..#|...|.|.|.|.| |###|..#|...|.|#|#|.| +-+-+.+.+.+.+-+.+-+.+ +-+-+.+#+#+.+-+.+-+.+ +-+-+.+#+#+.+-+#+-+.+ |.|.|.|...|.|.|.|...| |.|.|.|###|.|.|.|...| |.|.|.|###|.|.|#|...| +.+-+.+-+.+.+.+-+.+-+ +.+-+.+-+#+.+.+-+.+-+ +.+-+.+-+#+.+.+-+.+-+ |.|...|...|.|.....|.| |.|...|###|.|.....|.| |.|...|###|.|.....|.| +.+-+.+.+.+-+-+.+.+.+ +.+-+.+#+#+-+-+.+.+.+ +.+-+.+#+#+-+-+.+.+.+ |.|...|.|.....|.....| |.|...|#|####.|.....| |.|...|#|####.|.....| +-+.+-+.+-+.+-+.+-+-+ +-+.+-+#+-+#+-+#+-+-+ +-+.+-+#+-+#+-+#+-+-+ |.|.|.....|.....|...| |.|.|#####|#####|...| |.|.|#####|#####|...| +-+-+.+.+.+.+-+.+-+-+ +-+-+#+#+#+#+-+#+-+-+ +-+-+#+#+#+#+-+#+-+-+ |...|.|.....|.......| |...|#|#####|..##...| |...|#|#####|..##...| +-+-+-+-+-+-+-+.+-+-+ +-+-+-+-+-+-+-+#+-+-+ +-+-+-+-+-+-+-+#+-+-+ |.|...|.|.|.......|.| |.|...|.|.|######.|.| |.|...|.|.|######.|.| +.+-+-+-+.+.+-+.+.+.+ +.+-+-+-+.+#+-+#+.+.+ +.+-+-+-+.+#+-+#+.+.+ |.|...|.......|.|.|.| |.|...|....##.|#|.|.| |.|...|....##.|#|.|.| +.+.+-+.+.+.+-+-+-+-+ +.+.+-+.+.+#+-+-+-+-+ +.+.+-+.+.+#+-+-+-+-+ |.|.........|.....|.| |.|........#|.....|.| |.|........#|.....|.| +-+.+-+-+-+.+.+.+-+.+ +-+.+-+-+-+#+.+.+-+.+ +-+.+-+-+-+#+.+.+-+.+ After 2 attempt(s) +-+-+-+-+-+-+.+-+-+.+ +-+-+-+-+-+-+#+-+-+.+ +-+-+-+-+-+-+#+-+-+#+ |.|.|...|.|.|.|.|...| |.|.|...|.|.|#|.|...| |.|.|...|.|.|#|.|###| +-+-+-+-+-+-+.+-+-+-+ +-+-+-+-+-+-+#+-+-+-+ +-+-+-+-+-+-+#+-+-+-+ |.|.|.|...|.|...|.|.| |.|.|.|...|.|##.|.|.| |.|.|.|...|.|##.|.|.| +.+.+.+.+.+-+.+-+-+.+ +.+.+.+.+.+-+#+-+-+.+ +.+.+.+.+.+-+#+-+-+.+ |.|.|.|.|...|.|...|.| |.|.|.|.|...|#|...|.| |.|.|.|.|...|#|...|.| +.+-+.+-+-+-+.+-+.+.+ +.+-+.+-+-+-+#+-+.+.+ +.+-+.+-+-+-+#+-+.+.+ |...|...|.|.|...|.|.| |...|...|.|.|###|.|.| |...|...|.|.|###|.|.| +-+-+-+-+-+-+-+.+-+-+ +-+-+-+-+-+-+-+#+-+-+ +-+-+-+-+-+-+-+#+-+-+ |.|.......|.....|.|.| |.|.......|#####|.|.| |.|.......|#####|.|.| +.+-+-+-+.+.+-+-+-+-+ +.+-+-+-+.+#+-+-+-+-+ +.+-+-+-+.+#+-+-+-+-+ |.|.|.|.|.|.|.|.....| |.|.|.|.|.|#|.|.....| |.|.|.|.|.|#|.|.....| +-+-+-+-+-+.+.+-+.+-+ +-+-+-+-+-+#+.+-+.+-+ +-+-+-+-+-+#+.+-+.+-+ |...|.|.|.|.|.|.|.|.| |...|.|.|.|#|.|.|.|.| |...|.|.|.|#|.|.|.|.| +.+.+.+-+-+.+.+-+.+-+ +.+.+.+-+-+#+#+-+.+-+ +.+.+.+-+-+#+#+-+.+-+ |.|.|.|.|.|...|.|...| |.|.|.|.|.|###|.|...| |.|.|.|.|.|###|.|...| +-+-+-+-+-+-+.+.+.+-+ +-+-+-+-+-+-+#+.+.+-+ +-+-+-+-+-+-+#+.+.+-+ |.|.|.|.|.|...|...|.| |.|.|.|.|.|###|...|.| |.|.|.|.|.|###|...|.| +-+-+.+-+-+.+-+-+-+.+ +-+-+.+-+-+#+-+-+-+.+ +-+-+.+-+-+#+-+-+-+.+ |.|.|.|...|...|.|...| |.|.|.|...|###|.|...| |.|.|.|...|###|.|...| +-+-+.+-+.+-+.+.+.+-+ +-+-+.+-+.+-+#+.+.+-+ +-+-+.+-+.+-+#+.+.+-+ After 4611 attempt(s) > (main) proportion of grids that percolate p=0 : 1 (1.00000) proportion of grids that percolate p=1/10 : 1 (1.00000) proportion of grids that percolate p=1/5 : 1 (1.00000) proportion of grids that percolate p=3/10 : 199/200 (0.99500) proportion of grids that percolate p=2/5 : 179/200 (0.89500) proportion of grids that percolate p=1/2 : 451/1000 (0.45100) proportion of grids that percolate p=3/5 : 29/500 (0.05800) proportion of grids that percolate p=7/10 : 1/1000 (0.00100) proportion of grids that percolate p=4/5 : 0 (0.00000) proportion of grids that percolate p=9/10 : 0 (0.00000) proportion of grids that percolate p=1 : 0 (0.00000)
Tcl
<lang tcl>package require Tcl 8.6
- Structure the bond percolation system as a class
oo::class create BondPercolation {
variable hwall vwall cells M N constructor {width height probability} {
set M $height set N $width for {set i 0} {$i <= $height} {incr i} { for {set j 0;set walls {}} {$j < $width} {incr j} { lappend walls [expr {rand() < $probability}] } lappend hwall $walls } for {set i 0} {$i <= $height} {incr i} { for {set j 0;set walls {}} {$j <= $width} {incr j} { lappend walls [expr {$j==0 || $j==$width || rand() < $probability}] } lappend vwall $walls } set cells [lrepeat $height [lrepeat $width 0]]
}
method print Template:Percolated "" {
set nw [string length $M] set grid $cells if {$percolated ne ""} { lappend grid [lrepeat $N 0] lset grid end $percolated 1 } foreach hws $hwall vws [lrange $vwall 0 end-1] r $grid { incr row puts -nonewline [string repeat " " [expr {$nw+2}]] foreach w $hws { puts -nonewline [if {$w} {subst "+-"} {subst "+ "}] } puts "+" puts -nonewline [format "%-*s" [expr {$nw+2}] [expr { $row>$M ? $percolated eq "" ? " " : ">" : "$row)" }]] foreach v $vws c $r { puts -nonewline [if {$v==1} {subst "|"} {subst " "}] puts -nonewline [if {$c==1} {subst "#"} {subst " "}] } puts "" }
}
method percolate {} {
try { for {set i 0} {$i < $N} {incr i} { if {![lindex $hwall 0 $i]} { my FloodFill $i 0 } } return "" } trap PERCOLATED n { return $n }
} method FloodFill {x y} {
# fill cell lset cells $y $x 1 # bottom if {![lindex $hwall [expr {$y+1}] $x]} { if {$y == $N-1} { # THE bottom throw PERCOLATED $x } if {$y < $N-1 && ![lindex $cells [expr {$y+1}] $x]} { my FloodFill $x [expr {$y+1}] } } # left if {![lindex $vwall $y $x] && ![lindex $cells $y [expr {$x-1}]]} { my FloodFill [expr {$x-1}] $y } # right if {![lindex $vwall $y [expr {$x+1}]] && ![lindex $cells $y [expr {$x+1}]]} { my FloodFill [expr {$x+1}] $y } # top if {$y>0 && ![lindex $hwall $y $x] && ![lindex $cells [expr {$y-1}] $x]} { my FloodFill $x [expr {$y-1}] }
}
}
- Demonstrate one run
puts "Sample percolation, 10x10 p=0.5" BondPercolation create bp 10 10 0.5 bp print [bp percolate] bp destroy puts ""
- Collect some aggregate 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 [BondPercolation new 10 10 $p] if {[$bp percolate] ne ""} { incr tot } $bp destroy } puts [format "p=%.2f: %2.1f%%" $p [expr {$tot*100./$tries}]]
}
}}</lang>
- Output:
Sample percolation, 10x10 p=0.5 + + +-+-+-+ +-+ +-+ + 1) |# | | | | | + +-+ + + +-+ + + +-+ 2) |#| | | | | + + +-+ +-+ +-+ + +-+ 3) |# # #|# # #| | | | + +-+ + +-+ +-+ +-+ + 4) |#|# # #| |# | | +-+ + + +-+ +-+-+ +-+ 5) |# # # #| |# | | | +-+-+-+-+ + + + +-+-+ 6) | | | |#| | | +-+-+-+-+-+ + +-+-+ + 7) | | | | |# | | + +-+ +-+-+ +-+ +-+ + 8) | | # | | + +-+-+ +-+ + + + + + 9) | # | + + +-+-+ + +-+-+ + + 10) | | | # | | | + + + + + + +-+ +-+ + > # Percentage of tries that percolate, varying p p=0.00: 100.0% p=0.10: 100.0% p=0.20: 100.0% p=0.30: 100.0% p=0.40: 86.0% p=0.50: 50.0% p=0.60: 6.0% p=0.70: 0.0% p=0.80: 0.0% p=0.90: 0.0% p=1.00: 0.0%