Percolation/Mean run density: Difference between revisions
New draft task and Python solution. |
+ D entry |
||
Line 29: | Line 29: | ||
* [[Percolation/Bond percolation]] |
* [[Percolation/Bond percolation]] |
||
* [[Percolation/Site percolation]] |
* [[Percolation/Site percolation]] |
||
=={{header|D}}== |
|||
{{trans|python}} |
|||
<lang d>import std.stdio, std.range, std.algorithm, std.random, std.array, |
|||
std.math; |
|||
enum n = 100, p = 0.5, t = 500; |
|||
auto newVect(in size_t len, in double prob) { |
|||
return len.iota.map!(_ => uniform(0.0, 1.0) < prob).array; |
|||
} |
|||
size_t nRuns(R)(R vec) if (isForwardRange!R) { |
|||
return vec.group.filter!q{ a[0] }.walkLength; |
|||
} |
|||
double meanRunDensity(in size_t n, in double prob) { |
|||
return nRuns(newVect(n, prob)) / cast(double)n; |
|||
} |
|||
void main() { |
|||
foreach (immutable p10; iota(1, 10, 2)) { |
|||
immutable p = p10 / 10.0; |
|||
immutable limit = p * (1 - p); |
|||
writeln; |
|||
foreach (immutable n2; iota(10, 16, 2)) { |
|||
immutable n = 2 ^^ n2; |
|||
immutable sim = t.iota.map!(_ => meanRunDensity(n, p)) |
|||
//.sum / t; |
|||
.reduce!q{a + b} / t; |
|||
writefln("t=%3d, p=%4.2f, n=%5d, p(1-p)=%5.3f," ~ |
|||
" sim=%5.3f, delta=%3.1f%%", t, p, n, limit, sim, |
|||
limit ? abs(sim - limit) / limit * 100 : sim*100); |
|||
} |
|||
} |
|||
}</lang> |
|||
{{out}} |
|||
<pre>t=500, p=0.10, n= 1024, p(1-p)=0.090, sim=0.090, delta=0.5% |
|||
t=500, p=0.10, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.2% |
|||
t=500, p=0.10, n=16384, p(1-p)=0.090, sim=0.090, delta=0.1% |
|||
t=500, p=0.30, n= 1024, p(1-p)=0.210, sim=0.210, delta=0.2% |
|||
t=500, p=0.30, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.0% |
|||
t=500, p=0.30, n=16384, p(1-p)=0.210, sim=0.210, delta=0.0% |
|||
t=500, p=0.50, n= 1024, p(1-p)=0.250, sim=0.250, delta=0.1% |
|||
t=500, p=0.50, n= 4096, p(1-p)=0.250, sim=0.250, delta=0.1% |
|||
t=500, p=0.50, n=16384, p(1-p)=0.250, sim=0.250, delta=0.0% |
|||
t=500, p=0.70, n= 1024, p(1-p)=0.210, sim=0.210, delta=0.2% |
|||
t=500, p=0.70, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.0% |
|||
t=500, p=0.70, n=16384, p(1-p)=0.210, sim=0.210, delta=0.1% |
|||
t=500, p=0.90, n= 1024, p(1-p)=0.090, sim=0.090, delta=0.1% |
|||
t=500, p=0.90, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.1% |
|||
t=500, p=0.90, n=16384, p(1-p)=0.090, sim=0.090, delta=0.1%</pre> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
Revision as of 21:01, 11 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
Let v be a vector of n values of either 1 or 0 where the probability of any value being 1 is p, (and 0 is therefore 1-p).
Define a run of 1's as being a group of consecutive 1's in the vector bounded either by the limits of the vector or by a 0. Let the number of runs in a vector of length n be Rn.
The following vector has R10 = 3
[1 1 0 0 0 1 0 1 1 1]
Percolation theory states that
K(p) = Rn / n = p(1 - p) as n tends to infinity.
- Task
Any calculation of Rn / n for finite n is subject to randomnes so should be computed as the average of t runs, t >= 100.
for values of p of 0.1, 0.3, 0.5, 0.7, and 0.9; show the effect of varying n on the accuracy of simulated K(p).
Show your output here
- See
- s-Run on Wolfram mathworld.
- Percolation/Bond percolation
- Percolation/Site percolation
D
<lang d>import std.stdio, std.range, std.algorithm, std.random, std.array,
std.math;
enum n = 100, p = 0.5, t = 500;
auto newVect(in size_t len, in double prob) {
return len.iota.map!(_ => uniform(0.0, 1.0) < prob).array;
}
size_t nRuns(R)(R vec) if (isForwardRange!R) {
return vec.group.filter!q{ a[0] }.walkLength;
}
double meanRunDensity(in size_t n, in double prob) {
return nRuns(newVect(n, prob)) / cast(double)n;
}
void main() {
foreach (immutable p10; iota(1, 10, 2)) { immutable p = p10 / 10.0; immutable limit = p * (1 - p); writeln; foreach (immutable n2; iota(10, 16, 2)) { immutable n = 2 ^^ n2; immutable sim = t.iota.map!(_ => meanRunDensity(n, p)) //.sum / t; .reduce!q{a + b} / t; writefln("t=%3d, p=%4.2f, n=%5d, p(1-p)=%5.3f," ~ " sim=%5.3f, delta=%3.1f%%", t, p, n, limit, sim, limit ? abs(sim - limit) / limit * 100 : sim*100); } }
}</lang>
- Output:
t=500, p=0.10, n= 1024, p(1-p)=0.090, sim=0.090, delta=0.5% t=500, p=0.10, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.2% t=500, p=0.10, n=16384, p(1-p)=0.090, sim=0.090, delta=0.1% t=500, p=0.30, n= 1024, p(1-p)=0.210, sim=0.210, delta=0.2% t=500, p=0.30, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.0% t=500, p=0.30, n=16384, p(1-p)=0.210, sim=0.210, delta=0.0% t=500, p=0.50, n= 1024, p(1-p)=0.250, sim=0.250, delta=0.1% t=500, p=0.50, n= 4096, p(1-p)=0.250, sim=0.250, delta=0.1% t=500, p=0.50, n=16384, p(1-p)=0.250, sim=0.250, delta=0.0% t=500, p=0.70, n= 1024, p(1-p)=0.210, sim=0.210, delta=0.2% t=500, p=0.70, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.0% t=500, p=0.70, n=16384, p(1-p)=0.210, sim=0.210, delta=0.1% t=500, p=0.90, n= 1024, p(1-p)=0.090, sim=0.090, delta=0.1% t=500, p=0.90, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.1% t=500, p=0.90, n=16384, p(1-p)=0.090, sim=0.090, delta=0.1%
Python
<lang python>from __future__ import division from random import random from math import fsum
n, p, t = 100, 0.5, 500
def newv(n, p):
return [int(random() < p) for i in range(n)]
def runs(v):
return sum((a & ~b) for a, b in zip(v, v[1:] + [0]))
def mean_run_density(n, p):
return runs(newv(n, p)) / n
for p10 in range(1, 10, 2):
p = p10 / 10 limit = p * (1 - p) print() for n2 in range(10, 16, 2): n = 2**n2 sim = fsum(mean_run_density(n, p) for i in range(t)) / t print('t=%3i p=%4.2f n=%5i p(1-p)=%5.3f sim=%5.3f delta=%3.1f%%' % (t, p, n, limit, sim, abs(sim - limit) / limit * 100 if limit else sim * 100))</lang>
- Output:
t=500 p=0.10 n= 1024 p(1-p)=0.090 sim=0.090 delta=0.2% t=500 p=0.10 n= 4096 p(1-p)=0.090 sim=0.090 delta=0.0% t=500 p=0.10 n=16384 p(1-p)=0.090 sim=0.090 delta=0.1% t=500 p=0.30 n= 1024 p(1-p)=0.210 sim=0.210 delta=0.0% t=500 p=0.30 n= 4096 p(1-p)=0.210 sim=0.210 delta=0.0% t=500 p=0.30 n=16384 p(1-p)=0.210 sim=0.210 delta=0.0% t=500 p=0.50 n= 1024 p(1-p)=0.250 sim=0.251 delta=0.3% t=500 p=0.50 n= 4096 p(1-p)=0.250 sim=0.250 delta=0.0% t=500 p=0.50 n=16384 p(1-p)=0.250 sim=0.250 delta=0.0% t=500 p=0.70 n= 1024 p(1-p)=0.210 sim=0.210 delta=0.0% t=500 p=0.70 n= 4096 p(1-p)=0.210 sim=0.210 delta=0.1% t=500 p=0.70 n=16384 p(1-p)=0.210 sim=0.210 delta=0.0% t=500 p=0.90 n= 1024 p(1-p)=0.090 sim=0.091 delta=0.6% t=500 p=0.90 n= 4096 p(1-p)=0.090 sim=0.090 delta=0.2% t=500 p=0.90 n=16384 p(1-p)=0.090 sim=0.090 delta=0.0%