Percolation/Site percolation: Difference between revisions
m (correct D code to match requirements) |
m (correct C code to match requirements) |
||
Line 20: | Line 20: | ||
=={{header|C}}== |
=={{header|C}}== |
||
{{incomplete|C|"p going from 0.0 to 1.0 in 0.1 increments" (not 0.02 increments)}} |
|||
{{trans|D}} |
{{trans|D}} |
||
<lang c>#include <stdio.h> |
<lang c>#include <stdio.h> |
||
Line 27: | Line 26: | ||
#include <string.h> |
#include <string.h> |
||
#include <stdbool.h> |
#include <stdbool.h> |
||
#define N_COLS 15 |
#define N_COLS 15 |
||
#define N_ROWS 15 |
#define N_ROWS 15 |
||
// Probability granularity |
// Probability granularity 0.0, 0.1, ... 1.0 |
||
#define N_STEPS |
#define N_STEPS 11 |
||
// Simulation tries |
// Simulation tries |
||
#define N_TRIES |
#define N_TRIES 100 |
||
typedef unsigned char Cell; |
typedef unsigned char Cell; |
||
enum { EMPTY_CELL = ' ', |
enum { EMPTY_CELL = ' ', |
||
Line 42: | Line 41: | ||
VISITED_CELL = '.' }; |
VISITED_CELL = '.' }; |
||
typedef Cell Grid[N_ROWS][N_COLS]; |
typedef Cell Grid[N_ROWS][N_COLS]; |
||
void initialize(Grid grid, const double probability) { |
void initialize(Grid grid, const double probability) { |
||
for (size_t r = 0; r < N_ROWS; r++) |
for (size_t r = 0; r < N_ROWS; r++) |
||
Line 50: | Line 49: | ||
} |
} |
||
} |
} |
||
void show(Grid grid) { |
void show(Grid grid) { |
||
char line[N_COLS + 3]; |
char line[N_COLS + 3]; |
||
Line 57: | Line 56: | ||
line[N_COLS + 1] = '+'; |
line[N_COLS + 1] = '+'; |
||
line[N_COLS + 2] = '\0'; |
line[N_COLS + 2] = '\0'; |
||
printf("%s\n", line); |
printf("%s\n", line); |
||
for (size_t r = 0; r < N_ROWS; r++) { |
for (size_t r = 0; r < N_ROWS; r++) { |
||
Line 67: | Line 66: | ||
printf("%s\n", line); |
printf("%s\n", line); |
||
} |
} |
||
bool walk(Grid grid, const size_t r, const size_t c) { |
bool walk(Grid grid, const size_t r, const size_t c) { |
||
const size_t bottom = N_ROWS - 1; |
const size_t bottom = N_ROWS - 1; |
||
grid[r][c] = VISITED_CELL; |
grid[r][c] = VISITED_CELL; |
||
if (r < bottom && grid[r + 1][c] == EMPTY_CELL) { // Down. |
if (r < bottom && grid[r + 1][c] == EMPTY_CELL) { // Down. |
||
if (walk(grid, r + 1, c)) |
if (walk(grid, r + 1, c)) |
||
Line 77: | Line 76: | ||
} else if (r == bottom) |
} else if (r == bottom) |
||
return true; |
return true; |
||
if (c && grid[r][c - 1] == EMPTY_CELL) // Left. |
if (c && grid[r][c - 1] == EMPTY_CELL) // Left. |
||
if (walk(grid, r, c - 1)) |
if (walk(grid, r, c - 1)) |
||
return true; |
return true; |
||
if (c < N_COLS - 1 && grid[r][c + 1] == EMPTY_CELL) // Right. |
if (c < N_COLS - 1 && grid[r][c + 1] == EMPTY_CELL) // Right. |
||
if (walk(grid, r, c + 1)) |
if (walk(grid, r, c + 1)) |
||
return true; |
return true; |
||
if (r && grid[r - 1][c] == EMPTY_CELL) // Up. |
if (r && grid[r - 1][c] == EMPTY_CELL) // Up. |
||
if (walk(grid, r - 1, c)) |
if (walk(grid, r - 1, c)) |
||
return true; |
return true; |
||
return false; |
return false; |
||
} |
} |
||
bool percolate(Grid grid) { |
bool percolate(Grid grid) { |
||
const size_t startR = 0; |
const size_t startR = 0; |
||
Line 101: | Line 100: | ||
return false; |
return false; |
||
} |
} |
||
typedef struct { |
typedef struct { |
||
double prob; |
double prob; |
||
size_t count; |
size_t count; |
||
} Counter; |
} Counter; |
||
int main() { |
int main() { |
||
const double probability_step = 1.0 / N_STEPS; |
const double probability_step = 1.0 / (N_STEPS - 1); |
||
Counter counters[N_STEPS]; |
Counter counters[N_STEPS]; |
||
for (size_t i = 0; i < N_STEPS; i++) |
for (size_t i = 0; i < N_STEPS; i++) |
||
counters[i] = (Counter){ i * probability_step, 0 }; |
counters[i] = (Counter){ i * probability_step, 0 }; |
||
bool sample_shown = false; |
bool sample_shown = false; |
||
static Grid grid; |
static Grid grid; |
||
srand(time(NULL)); |
srand(time(NULL)); |
||
for (size_t i = 0; i < N_STEPS; i++) { |
for (size_t i = 0; i < N_STEPS; i++) { |
||
for (size_t t = 0; t < N_TRIES; t++) { |
for (size_t t = 0; t < N_TRIES; t++) { |
||
Line 132: | Line 131: | ||
} |
} |
||
} |
} |
||
printf("\nFraction of %d tries that percolate through:\n", N_TRIES); |
printf("\nFraction of %d tries that percolate through:\n", N_TRIES); |
||
for (size_t i = 0; i < N_STEPS; i++) |
for (size_t i = 0; i < N_STEPS; i++) |
||
printf("%1. |
printf("%1.1f %1.3f\n", counters[i].prob, |
||
counters[i].count / (double)N_TRIES); |
counters[i].count / (double)N_TRIES); |
||
return 0; |
return 0; |
||
} |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre>Percolating sample (15x15, probability = 0. |
<pre>Percolating sample (15x15, probability = 0.40): |
||
+---------------+ |
+---------------+ |
||
| |
|###. # # # #| |
||
|### ## |
|###.. # ##### | |
||
| |
| #. ###### #| |
||
| |
|###.... ######| |
||
|## |
|######. ### # | |
||
| ## |
| #####.###### | |
||
| |
|#......... ## | |
||
|# |
|...#...##.# ## | |
||
| |
|##.#...##.### #| |
||
| |
| ###..# #. # | |
||
|##### |
|# #######. # ##| |
||
|# |
| # ##...#### | |
||
| # |
| ## # .##### | |
||
| |
|#######.## ###| |
||
| |
|# ## .## # # | |
||
+---------------+ |
+---------------+ |
||
Fraction of |
Fraction of 100 tries that percolate through: |
||
0. |
0.0 0.000 |
||
0. |
0.1 0.000 |
||
0. |
0.2 0.000 |
||
0. |
0.3 0.000 |
||
0. |
0.4 0.010 |
||
0. |
0.5 0.070 |
||
0. |
0.6 0.630 |
||
0. |
0.7 0.970 |
||
0. |
0.8 1.000 |
||
0. |
0.9 1.000 |
||
0 |
1.0 1.000</pre> |
||
0.220 0.000 |
|||
0.240 0.000 |
|||
0.260 0.000 |
|||
0.280 0.000 |
|||
0.300 0.000 |
|||
0.320 0.000 |
|||
0.340 0.000 |
|||
0.360 0.000 |
|||
0.380 0.001 |
|||
0.400 0.004 |
|||
0.420 0.007 |
|||
0.440 0.015 |
|||
0.460 0.029 |
|||
0.480 0.056 |
|||
0.500 0.093 |
|||
0.520 0.155 |
|||
0.540 0.226 |
|||
0.560 0.332 |
|||
0.580 0.437 |
|||
0.600 0.564 |
|||
0.620 0.680 |
|||
0.640 0.783 |
|||
0.660 0.862 |
|||
0.680 0.920 |
|||
0.700 0.959 |
|||
0.720 0.981 |
|||
0.740 0.991 |
|||
0.760 0.997 |
|||
0.780 0.999 |
|||
0.800 1.000 |
|||
0.820 1.000 |
|||
0.840 1.000 |
|||
0.860 1.000 |
|||
0.880 1.000 |
|||
0.900 1.000 |
|||
0.920 1.000 |
|||
0.940 1.000 |
|||
0.960 1.000 |
|||
0.980 1.000 |
|||
</pre> |
|||
=={{header|D}}== |
=={{header|D}}== |
Revision as of 21:56, 28 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
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 successfull percolation path through a cell grid graphically.
Show all output on this page.
C
<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: 0.0, 0.1 ... 1.0 nTries = 100; // 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 / cast(double)rng.max; rng.popFront; cell = (r < probability) ? Cell.empty : Cell.filled; }
}
void show(in ref Grid grid) {
immutable static line = '+' ~ "-".replicate(nCols) ~ "+"; line.writeln; foreach (const ref row; grid) writeln('|', cast(BaseType[nCols])row, '|'); line.writeln;
}
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.1f %1.3f", co.prob, co.count / cast(double)nTries); writefln("\nSimulations and grid printing performed" ~ " in %2.3f seconds.", sw.peek.msecs / 1000.0);
}</lang>
- Output:
Percolating sample (15x15, probability = 0.50): +---------------+ |.....#. #### | |...#.#.# # # ##| |####.#.....# #| |##...#.#.#. #| |..##..#..#. # | |#..#.##.##..# #| |#....# ## #. # | |##.#####.... | |....###..##.## | |#.#### .# #### | |#.##.##. ##### | |#....# .## ##| | ######.## ## #| | # ## .# # | | ## # .### ###| +---------------+ 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.000 0.5 0.090 0.6 0.570 0.7 0.950 0.8 1.000 0.9 1.000 1.0 1.000 Simulations and grid printing performed in 0.046 seconds.
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.