Percolation/Mean run density
Let be a vector of values of either 1 or 0 where the probability of any value being 1 is , (and 0 is therefore ). 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 be .
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
The following vector has
[1 1 0 0 0 1 0 1 1 1] ^^^ ^ ^^^^^
Percolation theory states that
- Task
Any calculation of for finite is subject to randomness so should be computed as the average of runs, where .
For values of of 0.1, 0.3, 0.5, 0.7, and 0.9, show the effect of varying on the accuracy of simulated .
Show your output here
- See also
- s-Run on Wolfram mathworld.
C++
<lang C++>#include <algorithm>
- include <random>
- include <vector>
- include <iostream>
- include <numeric>
- include <iomanip>
using VecIt = std::vector<int>::const_iterator ;
//creates vector of length n, based on probability p for 1 std::vector<int> createVector( int n, double p ) {
std::vector<int> result( n ) ; std::random_device rd ; std::mt19937 gen( rd( ) ) ; std::uniform_real_distribution<> dis( 0 , 1 ) ; for ( int i = 0 ; i < n ; i++ ) { double number = dis( gen ) ; if ( number <= p )
result[ i ] = 1 ;
else
result[ i ] = 0 ;
} return result ;
}
//find number of 1 runs in the vector int find_Runs( const std::vector<int> & numberVector ) {
int runs = 0 ; VecIt found = numberVector.begin( ) ; while ( ( found = std::find( found , numberVector.end( ) , 1 ) )
!= numberVector.end( ) ) {
runs++ ; while ( found != numberVector.end( ) && ( *found == 1 ) )
std::advance( found , 1 ) ;
if ( found == numberVector.end( ) )
break ;
} return runs ;
}
int main( ) {
std::cout << "t = 100\n" ; std::vector<double> p_values { 0.1 , 0.3 , 0.5 , 0.7 , 0.9 } ; for ( double p : p_values ) { std::cout << "p = " << p << " , K(p) = " << p * ( 1 - p ) << std::endl ; for ( int n = 10 ; n < 100000 ; n *= 10 ) {
std::vector<double> runsFound ; for ( int i = 0 ; i < 100 ; i++ ) { std::vector<int> ones_and_zeroes = createVector( n , p ) ; runsFound.push_back( find_Runs( ones_and_zeroes ) / static_cast<double>( n ) ) ; } double average = std::accumulate( runsFound.begin( ) , runsFound.end( ) , 0.0 ) / runsFound.size( ) ; std::cout << " R(" << std::setw( 6 ) << std::right << n << ", p) = " << average << std::endl ;
} } return 0 ;
}</lang>
- Output:
t = 100 p = 0.1 , K(p) = 0.09 R( 10, p) = 0.088 R( 100, p) = 0.0931 R( 1000, p) = 0.09013 R( 10000, p) = 0.089947 p = 0.3 , K(p) = 0.21 R( 10, p) = 0.225 R( 100, p) = 0.2089 R( 1000, p) = 0.21043 R( 10000, p) = 0.20991 p = 0.5 , K(p) = 0.25 R( 10, p) = 0.271 R( 100, p) = 0.253 R( 1000, p) = 0.25039 R( 10000, p) = 0.250278 p = 0.7 , K(p) = 0.21 R( 10, p) = 0.264 R( 100, p) = 0.2155 R( 1000, p) = 0.20829 R( 10000, p) = 0.209977 p = 0.9 , K(p) = 0.09 R( 10, p) = 0.167 R( 100, p) = 0.0928 R( 1000, p) = 0.09071 R( 10000, p) = 0.090341
D
<lang d>import std.stdio, std.range, std.algorithm, std.random, std.math;
enum n = 100, p = 0.5, t = 500;
double meanRunDensity(in size_t n, in double prob) {
//return n.iota.map!(_ => uniform01 < prob).array // .uniq.sum / cast(double)n; return n.iota.map!(_ => uniform(0.0, 1.0) < prob).array .uniq.count!q{a} / cast(double)n;
}
void main() {
foreach (immutable p; iota(0.1, 1.0, 0.2)) { 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.5f, " ~ "sim=%5.5f, 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.09000, sim=0.08949, delta=0.6% t=500, p=0.10, n= 4096, p(1-p)=0.09000, sim=0.08976, delta=0.3% t=500, p=0.10, n=16384, p(1-p)=0.09000, sim=0.08988, delta=0.1% t=500, p=0.30, n= 1024, p(1-p)=0.21000, sim=0.20979, delta=0.1% t=500, p=0.30, n= 4096, p(1-p)=0.21000, sim=0.21020, delta=0.1% t=500, p=0.30, n=16384, p(1-p)=0.21000, sim=0.21005, delta=0.0% t=500, p=0.50, n= 1024, p(1-p)=0.25000, sim=0.25016, delta=0.1% t=500, p=0.50, n= 4096, p(1-p)=0.25000, sim=0.25026, delta=0.1% t=500, p=0.50, n=16384, p(1-p)=0.25000, sim=0.24990, delta=0.0% t=500, p=0.70, n= 1024, p(1-p)=0.21000, sim=0.21050, delta=0.2% t=500, p=0.70, n= 4096, p(1-p)=0.21000, sim=0.20993, delta=0.0% t=500, p=0.70, n=16384, p(1-p)=0.21000, sim=0.21009, delta=0.0% t=500, p=0.90, n= 1024, p(1-p)=0.09000, sim=0.09019, delta=0.2% t=500, p=0.90, n= 4096, p(1-p)=0.09000, sim=0.09047, delta=0.5% t=500, p=0.90, n=16384, p(1-p)=0.09000, sim=0.09007, delta=0.1%
FORTRAN
<lang> ! loosely translated from python. We do not need to generate and store the entire vector at once. ! compilation: gfortran -Wall -std=f2008 -o thisfile thisfile.f08
program percolation_mean_run_density
implicit none integer :: i, p10, n2, n, t real :: p, limit, sim, delta data n,p,t/100,0.5,500/ write(6,'(a3,a5,4a7)')'t','p','n','p(1-p)','sim','delta'!a little tight! do p10=1,10,2 p = p10/10.0 limit = p*(1-p) write(6,'()') do n2=10,15,2 n = 2**n2 sim = 0 do i=1,t sim = sim + mean_run_density(n,p) end do sim = sim/t if (limit /= 0) then delta = abs(sim-limit)/limit else delta = sim end if delta = delta * 100 write(6,'(i3,f5.2,i7,2f7.3,f5.1)')t,p,n,limit,sim,delta end do end do
contains
integer function runs(n, p) integer, intent(in) :: n real, intent(in) :: p real :: harvest logical :: q integer :: count, i count = 0 q = .false. do i=1,n call random_number(harvest) if (harvest < p) then q = .true. else if (q) count = count+1 q = .false. end if end do runs = count end function runs
real function mean_run_density(n, p) integer, intent(in) :: n real, intent(in) :: p mean_run_density = real(runs(n,p))/real(n) end function mean_run_density
end program percolation_mean_run_density </lang>
$ ./f t p n p(1-p) sim delta 500 0.10 1024 0.090 0.090 0.2 500 0.10 4096 0.090 0.090 0.2 500 0.10 16384 0.090 0.090 0.0 500 0.30 1024 0.210 0.210 0.2 500 0.30 4096 0.210 0.210 0.0 500 0.30 16384 0.210 0.210 0.0 500 0.50 1024 0.250 0.250 0.1 500 0.50 4096 0.250 0.250 0.1 500 0.50 16384 0.250 0.250 0.1 500 0.70 1024 0.210 0.210 0.1 500 0.70 4096 0.210 0.210 0.1 500 0.70 16384 0.210 0.210 0.0 500 0.90 1024 0.090 0.090 0.1 500 0.90 4096 0.090 0.090 0.4 500 0.90 16384 0.090 0.090 0.1
Icon and Unicon
The following works in both languages:
<lang unicon>procedure main(A)
t := integer(A[2]) | 500
write(left("p",8)," ",left("n",8)," ",left("p(1-p)",10)," ",left("SimK(p)",10)) every (p := 0.1 | 0.3 | 0.5 | 0.7 | 0.9, n := 1000 | 2000 | 3000) do { Ka := 0.0 every !t do { every (v := "", !n) do v ||:= |((?0.1 > p,"0")|"1") R := 0 v ? while tab(upto('1')) do R +:= (tab(many('1')), 1) Ka +:= real(R)/n } write(left(p,8)," ",left(n,8)," ",left(p*(1-p),10)," ",left(Ka/t, 10)) }
end</lang>
Output:
->pmrd p n p(1-p) SimK(p) 0.1 1000 0.09000000 0.09021400 0.1 2000 0.09000000 0.08984799 0.1 3000 0.09000000 0.08993666 0.3 1000 0.21 0.21080999 0.3 2000 0.21 0.209953 0.3 3000 0.21 0.210564 0.5 1000 0.25 0.250024 0.5 2000 0.25 0.25007399 0.5 3000 0.25 0.24975266 0.7 1000 0.21 0.21098799 0.7 2000 0.21 0.20987700 0.7 3000 0.21 0.21047333 0.9 1000 0.08999999 0.09016400 0.9 2000 0.08999999 0.09004800 0.9 3000 0.08999999 0.09023200 ->
J
<lang J> NB. translation of python
NB. 'N P T' =: 100 0.5 500 NB. silliness
newv =: (> ?@(#&0))~ NB. generate a random binary vector. Use: N newv P runs =: {: + [: +/ 1 0&E. NB. add the tail to the sum of 1 0 occurrences Use: runs V mean_run_density =: [ %~ [: runs newv NB. perform experiment. Use: N mean_run_density P
main =: 3 : 0 NB. Use: main T
T =. y smoutput' T P N P(1-P) SIM DELTA' for_P. 10 %~ >: +: i. 5 do. LIMIT =. (* -.) P smoutput for_N. 2 ^ 10 + +: i. 3 do. SIM =. T %~ +/ ".@:('N mean_run_density P'"_)^:(<T)0 smoutput 4 5j2 6 6j3 6j3 4j1 ": T, P, N, LIMIT, SIM, SIM (100 * [)`(|@:(- % ]))@.(0 ~: ]) LIMIT end. end. EMPTY
) </lang> Session:
main 9 T P N P(1-P) SIM DELTA 9 0.10 1024 0.090 0.081 0.1 9 0.10 4096 0.090 0.079 0.1 9 0.10 16384 0.090 0.079 0.1 9 0.30 1024 0.210 0.190 0.1 9 0.30 4096 0.210 0.186 0.1 9 0.30 16384 0.210 0.187 0.1 9 0.50 1024 0.250 0.220 0.1 9 0.50 4096 0.250 0.224 0.1 9 0.50 16384 0.250 0.223 0.1 9 0.70 1024 0.210 0.186 0.1 9 0.70 4096 0.210 0.186 0.1 9 0.70 16384 0.210 0.189 0.1 9 0.90 1024 0.090 0.081 0.1 9 0.90 4096 0.090 0.080 0.1 9 0.90 16384 0.090 0.079 0.1 main 500 T P N P(1-P) SIM DELTA 500 0.10 1024 0.090 0.089 0.0 500 0.10 4096 0.090 0.090 0.0 500 0.10 16384 0.090 0.090 0.0 500 0.30 1024 0.210 0.210 0.0 500 0.30 4096 0.210 0.210 0.0 500 0.30 16384 0.210 0.210 0.0 500 0.50 1024 0.250 0.250 0.0 500 0.50 4096 0.250 0.250 0.0 500 0.50 16384 0.250 0.250 0.0 500 0.70 1024 0.210 0.210 0.0 500 0.70 4096 0.210 0.210 0.0 500 0.70 16384 0.210 0.210 0.0 500 0.90 1024 0.090 0.091 0.0 500 0.90 4096 0.090 0.090 0.0 500 0.90 16384 0.090 0.090 0.0
Perl
<lang perl>sub R {
my ($n, $p) = @_; my $r = join , map { rand() < $p ? 1 : 0 } 1 .. $n; 0+ $r =~ s/1+//g;
}
use constant t => 100;
printf "t= %d\n", t; for my $p (qw(.1 .3 .5 .7 .9)) {
printf "p= %f, K(p)= %f\n", $p, $p*(1-$p); for my $n (qw(10 100 1000)) { my $r; $r += R($n, $p) for 1 .. t; $r /= $n; printf " R(n, p)= %f\n", $r / t; }
}</lang>
- Output:
t= 100 p= 0.100000, K(p)= 0.090000 R(n, p)= 0.095000 R(n, p)= 0.088100 R(n, p)= 0.089420 p= 0.300000, K(p)= 0.210000 R(n, p)= 0.225000 R(n, p)= 0.208800 R(n, p)= 0.210020 p= 0.500000, K(p)= 0.250000 R(n, p)= 0.289000 R(n, p)= 0.249900 R(n, p)= 0.248980 p= 0.700000, K(p)= 0.210000 R(n, p)= 0.262000 R(n, p)= 0.213200 R(n, p)= 0.209690 p= 0.900000, K(p)= 0.090000 R(n, p)= 0.177000 R(n, p)= 0.096200 R(n, p)= 0.091730
Perl 6
<lang perl6>sub R($n, $p) { [+] ((rand < $p) xx $n).squish }
say 't= ', constant t = 100;
for .1, .3 ... .9 -> $p {
say "p= $p, K(p)= {$p*(1-$p)}"; for 10, 100, 1000 -> $n {
printf " R(%6d, p)= %f\n", $n, t R/ [+] R($n, $p)/$n xx t
}
}</lang>
- Output:
t= 100 p= 0.1, K(p)= 0.09 R( 10, p)= 0.088000 R( 100, p)= 0.085600 R( 1000, p)= 0.089150 p= 0.3, K(p)= 0.21 R( 10, p)= 0.211000 R( 100, p)= 0.214600 R( 1000, p)= 0.211160 p= 0.5, K(p)= 0.25 R( 10, p)= 0.279000 R( 100, p)= 0.249200 R( 1000, p)= 0.250870 p= 0.7, K(p)= 0.21 R( 10, p)= 0.258000 R( 100, p)= 0.215400 R( 1000, p)= 0.209560 p= 0.9, K(p)= 0.09 R( 10, p)= 0.181000 R( 100, p)= 0.094500 R( 1000, p)= 0.091330
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%
Racket
<lang racket>#lang racket (require racket/fixnum) (define t (make-parameter 100))
(define (Rn v)
(define (inner-Rn rv idx b-1) (define b (fxvector-ref v idx)) (define rv+ (if (and (= b 1) (= b-1 0)) (add1 rv) rv)) (if (zero? idx) rv+ (inner-Rn rv+ (sub1 idx) b))) (inner-Rn 0 (sub1 (fxvector-length v)) 0))
(define ((make-random-bit-vector p) n)
(for/fxvector #:length n ((i n)) (if (<= (random) p) 1 0)))
(define (Rn/n l->p n) (/ (Rn (l->p n)) n))
(for ((p (in-list '(1/10 3/10 1/2 7/10 9/10))))
(define l->p (make-random-bit-vector p)) (define Kp (* p (- 1 p))) (printf "p = ~a\tK(p) =\t~a\t~a~%" p Kp (real->decimal-string Kp 4)) (for ((n (in-list '(10 100 1000 10000)))) (define sum-Rn/n (for/sum ((i (in-range (t)))) (Rn/n l->p n))) (define sum-Rn/n/t (/ sum-Rn/n (t))) (printf "mean(R_~a/~a) =\t~a\t~a~%" n n sum-Rn/n/t (real->decimal-string sum-Rn/n/t 4))) (newline))
(module+ test
(require rackunit) (check-eq? (Rn (fxvector 1 1 0 0 0 1 0 1 1 1)) 3))</lang>
- Output:
p = 1/10 K(p) = 9/100 0.0900 mean(R_10/10) = 3/40 0.0750 mean(R_100/100) = 221/2500 0.0884 mean(R_1000/1000) = 4469/50000 0.0894 mean(R_10000/10000) = 90313/1000000 0.0903 p = 3/10 K(p) = 21/100 0.2100 mean(R_10/10) = 231/1000 0.2310 mean(R_100/100) = 1049/5000 0.2098 mean(R_1000/1000) = 131/625 0.2096 mean(R_10000/10000) = 209873/1000000 0.2099 p = 1/2 K(p) = 1/4 0.2500 mean(R_10/10) = 297/1000 0.2970 mean(R_100/100) = 1263/5000 0.2526 mean(R_1000/1000) = 24893/100000 0.2489 mean(R_10000/10000) = 124963/500000 0.2499 p = 7/10 K(p) = 21/100 0.2100 mean(R_10/10) = 131/500 0.2620 mean(R_100/100) = 2147/10000 0.2147 mean(R_1000/1000) = 1049/5000 0.2098 mean(R_10000/10000) = 210453/1000000 0.2105 p = 9/10 K(p) = 9/100 0.0900 mean(R_10/10) = 169/1000 0.1690 mean(R_100/100) = 119/1250 0.0952 mean(R_1000/1000) = 4503/50000 0.0901 mean(R_10000/10000) = 89939/1000000 0.0899
Tcl
<lang tcl>proc randomString {length probability} {
for {set s ""} {[string length $s] < $length} {} {
append s [expr {rand() < $probability}]
} return $s
}
- By default, [regexp -all] gives the number of times that the RE matches
proc runs {str} {
regexp -all {1+} $str
}
- Compute the mean run density
proc mrd {t p n} {
for {set i 0;set total 0.0} {$i < $t} {incr i} {
set run [randomString $n $p] set total [expr {$total + double([runs $run])/$n}]
} return [expr {$total / $t}]
}
- Parameter sweep with nested [foreach]
set runs 500 foreach p {0.10 0.30 0.50 0.70 0.90} {
foreach n {1024 4096 16384} {
set theory [expr {$p * (1 - $p)}] set sim [mrd $runs $p $n] set diffpc [expr {abs($theory-$sim)*100/$theory}] puts [format "t=%d, p=%.2f, n=%5d, p(1-p)=%.3f, sim=%.3f, delta=%.2f%%" \ $runs $p $n $theory $sim $diffpc]
} puts ""
}</lang>
- Output:
t=500, p=0.10, n= 1024, p(1-p)=0.090, sim=0.090, delta=0.07% t=500, p=0.10, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.06% t=500, p=0.10, n=16384, p(1-p)=0.090, sim=0.090, delta=0.17% t=500, p=0.30, n= 1024, p(1-p)=0.210, sim=0.210, delta=0.23% t=500, p=0.30, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.09% t=500, p=0.30, n=16384, p(1-p)=0.210, sim=0.210, delta=0.01% t=500, p=0.50, n= 1024, p(1-p)=0.250, sim=0.250, delta=0.10% t=500, p=0.50, n= 4096, p(1-p)=0.250, sim=0.250, delta=0.07% t=500, p=0.50, n=16384, p(1-p)=0.250, sim=0.250, delta=0.08% t=500, p=0.70, n= 1024, p(1-p)=0.210, sim=0.211, delta=0.33% t=500, p=0.70, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.00% t=500, p=0.70, n=16384, p(1-p)=0.210, sim=0.210, delta=0.01% t=500, p=0.90, n= 1024, p(1-p)=0.090, sim=0.091, delta=1.61% t=500, p=0.90, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.08% t=500, p=0.90, n=16384, p(1-p)=0.090, sim=0.090, delta=0.09%