Percolation/Mean cluster density: Difference between revisions

Content added Content deleted

Inline

Revision as of 16:02, 17 March 2014

Let $c$ be a 2D boolean square matrix of $n\times n$ values of either 1 or 0 where the probability of any value being 1 is $p$ , (and of 0 is therefore $1-p$ ).

Define a cluster of 1's as being a group of 1's connected vertically or horizontally (i.e., using the Von Neumann neighborhood rule) and bounded by either $0$ or by the limits of the matrix. Let the number of such clusters in such a randomly constructed matrix be $C_{n}$ .

Percolation theory states that $K(p)=C_{n}/n^{2}$ as $n$ tends to infinity is a constant.

$K(p)$ is called the mean cluster density and for $p=0.5$ is found numerically to approximate $0.065770$ ...

Task

Any calculation of $C_{n}$ for finite $n$ is subject to randomnes so should be computed as the average of $t$ runs, where $t$ ≥ $5$ .
Show the effect of varying $n$ on the accuracy of simulated $K(p)$ for $p=0.5$ and for values of $n$ up to at least $1000$ .

For extra credit, graphically show clusters in a $15\times 15$ , $p=0.5$ grid.

Show your output here.

See also

s-Cluster on Wolfram mathworld.

C

<lang c>#include <stdio.h>

include <stdlib.h>

int *map, w, ww;

void make_map(double p) { int i, thresh = RAND_MAX * p; i = ww = w * w;

map = realloc(map, i * sizeof(int)); while (i--) map[i] = -(rand() < thresh); }

char alpha[] = "+.ABCDEFGHIJKLMNOPQRSTUVWXYZ" "abcdefghijklmnopqrstuvwxyz";

define ALEN ((int)(sizeof(alpha) - 3))

void show_cluster(void) { int i, j, *s = map;

for (i = 0; i < w; i++) { for (j = 0; j < w; j++, s++) printf(" %c", *s < ALEN ? alpha[1 + *s] : '?'); putchar('\n'); } }

void recur(int x, int v) { if (x >= 0 && x < ww && map[x] == -1) { map[x] = v; recur(x - w, v); recur(x - 1, v); recur(x + 1, v); recur(x + w, v); } }

int count_clusters(void) { int i, cls;

for (cls = i = 0; i < ww; i++) { if (-1 != map[i]) continue; recur(i, ++cls); }

return cls; }

double tests(int n, double p) { int i; double k;

for (k = i = 0; i < n; i++) { make_map(p); k += (double)count_clusters() / ww; } return k / n; }

int main(void) { w = 15; make_map(.5); printf("width=15, p=0.5, %d clusters:\n", count_clusters()); show_cluster();

printf("\np=0.5, iter=5:\n"); for (w = 1<<2; w <= 1<<14; w<<=2) printf("%5d %9.6f\n", w, tests(5, .5));

free(map); return 0; }</lang>

Output:

width=15, p=0.5, 23 clusters:
 A . . . B . C C C C . D . E .
 A . . B B . . . . . . . . . .
 A . . . . . F . . G . H . . I
 . . J J J . . K K . L . M . I
 . J J . . . K K K K . M M . .
 . . . . K K . K . K . M . N .
 O O . K K K . K . . . . N N N
 . O O . K K K K K . P . N . .
 Q . . K K . . . K . P . . . .
 . R . K K . . K K . P . . S .
 . . K K . . . . K . P . . . K
 K K K K K . . K K . . T . . K
 K . K . . . U . K . . T . . .
 K . K K K K . K K K . T . . .
 . . K . K . V . K K . . . W .

p=0.5, iter=5:
    4  0.125000
   16  0.083594
   64  0.064453
  256  0.066864
 1024  0.065922
 4096  0.065836
16384  0.065774

D

Translation of: python

<lang d>import std.stdio, std.algorithm, std.random, std.math, std.array,

      std.range, std.ascii;

alias Cell = ubyte; alias Grid = Cell[][]; enum Cell notClustered = 1; // Filled cell, but not in a cluster.

FP rand(FP = double, UniformRandomNumberGenerator) //

           (ref UniformRandomNumberGenerator urng) {
   immutable FP result = urng.front / FP(urng.max);
   urng.popFront;
   return result;

}

Grid initialize(Grid grid, in double prob, ref Xorshift rng) /*nothrow*/ {

   foreach (row; grid)
       foreach (ref cell; row)
           cell = cast(Cell)(rng.rand < prob);
   return grid;

}

void show(in Grid grid) {

   immutable static cell2char = " #" ~ letters;
   writeln('+', "-".replicate(grid.length), '+');
   foreach (row; grid) {
       write('|');
       row.map!(c => c < cell2char.length ? cell2char[c] : '@').write;
       writeln('|');
   }
   writeln('+', "-".replicate(grid.length), '+');

}

size_t countClusters(bool justCount=false)(Grid grid) pure nothrow {

   immutable side = grid.length;
   static if (justCount)
       enum Cell clusterID = 2;
   else
       Cell clusterID = 1;

   void walk(in size_t r, in size_t c) nothrow {
       grid[r][c] = clusterID; // Fill grid.

       if (r < side - 1 && grid[r + 1][c] == notClustered) // Down.
           walk(r + 1, c);
       if (c < side - 1 && grid[r][c + 1] == notClustered) // Right.
           walk(r, c + 1);
       if (c > 0 && grid[r][c - 1] == notClustered) // Left.
           walk(r, c - 1);
       if (r > 0 && grid[r - 1][c] == notClustered) // Up.
           walk(r - 1, c);
   }

   size_t nClusters = 0;

   foreach (immutable r; 0 .. side)
       foreach (immutable c; 0 .. side)
           if (grid[r][c] == notClustered) {
               static if (!justCount)
                   clusterID++;
               nClusters++;
               walk(r, c);
           }
   return nClusters;

}

double clusterDensity(Grid grid, in double prob, ref Xorshift rng) {

   return grid.initialize(prob, rng).countClusters!true /
          double(grid.length ^^ 2);

}

void showDemo(in size_t side, in double prob, ref Xorshift rng) {

   auto grid = new Grid(side, side);
   grid.initialize(prob, rng);
   writefln("Found %d clusters in this %d by %d grid:\n",
            grid.countClusters, side, side);
   grid.show;

}

void main() {

   immutable prob = 0.5;
   immutable nIters = 5;
   auto rng = Xorshift(unpredictableSeed);

   showDemo(15, prob, rng);
   writeln;
   foreach (immutable i; iota(4, 14, 2)) {
       immutable side = 2 ^^ i;
       auto grid = new Grid(side, side);
       immutable density = nIters
                           .iota
                           .map!(_ => grid.clusterDensity(prob, rng))
                           .sum / nIters;
       writefln("n_iters=%3d, p=%4.2f, n=%5d, sim=%7.8f",
                nIters, prob, side, density);
   }

}</lang>

Output:

Found 26 clusters in this 15 by 15 grid:

+---------------+
| AA B    CCCC  |
|AA D E F CC  G |
|  DDD FF  CC  H|
| I D FF  J   K |
|  L  FF JJJJ   |
|L LLL      J  M|
|LLLLLL   JJJ MM|
|L LL L N  J   M|
|LL    O P J   M|
|LLL QQ R JJ  S |
|LL T  RR  J SSS|
| L   U  V JJ  S|
|  WW  XX JJ YY |
|    XXX   JJ YY|
|ZZ   XXX   JJ  |
+---------------+

n_iters=  5, p=0.50, n=   16, sim=0.09765625
n_iters=  5, p=0.50, n=   64, sim=0.07260742
n_iters=  5, p=0.50, n=  256, sim=0.06679993
n_iters=  5, p=0.50, n= 1024, sim=0.06609497
n_iters=  5, p=0.50, n= 4096, sim=0.06580237

Increasing the index i to 15:

n_iters=  5, p=0.50, n=32768, sim=0.06578374

Python

<lang python>from __future__ import division from random import random import string from math import fsum

n_range, p, t = (2**n2 for n2 in range(4, 14, 2)), 0.5, 5 N = M = 15

NOT_CLUSTERED = 1 # filled but not clustered cell cell2char = ' #' + string.ascii_letters

def newgrid(n, p):

   return [[int(random() < p) for x in range(n)] for y in range(n)]

def pgrid(cell):

   for n in range(N):
       print( '%i)  ' % (n % 10) 
              + ' '.join(cell2char[cell[n][m]] for m in range(M)))

def cluster_density(n, p):

   cc = clustercount(newgrid(n, p))
   return cc / n / n

def clustercount(cell):

   walk_index = 1
   for n in range(N):
       for m in range(M):
           if cell[n][m] == NOT_CLUSTERED:
               walk_index += 1
               walk_maze(m, n, cell, walk_index)
   return walk_index - 1

def walk_maze(m, n, cell, indx):

   # fill cell 
   cell[n][m] = indx
   # down
   if n < N - 1 and cell[n+1][m] == NOT_CLUSTERED:
       walk_maze(m, n+1, cell, indx)
   # right
   if m < M - 1 and cell[n][m + 1] == NOT_CLUSTERED:
       walk_maze(m+1, n, cell, indx)
   # left
   if m and cell[n][m - 1] == NOT_CLUSTERED:
       walk_maze(m-1, n, cell, indx)
   # up
   if n and cell[n-1][m] == NOT_CLUSTERED:
       walk_maze(m, n-1, cell, indx)

if __name__ == '__main__':

   cell = newgrid(n=N, p=0.5)
   print('Found %i clusters in this %i by %i grid\n' 
         % (clustercount(cell), N, N))
   pgrid(cell)
   print()
   
   for n in n_range:
       N = M = n
       sim = fsum(cluster_density(n, p) for i in range(t)) / t
       print('t=%3i p=%4.2f n=%5i sim=%7.5f'
             % (t, p, n, sim))</lang>

Output:

Found 20 clusters in this 15 by 15 grid

0)  a a     b     c       d d d d
1)  a a   e     f   g g   d      
2)        e   f f f       d      
3)  h h   e     f   i i   d d    
4)        e       j     d d d d  
5)    k k   k   k   l         d d
6)  k k k k k   k   l   m   n    
7)  k k k   k   k   l     o   p p
8)      k   k k k   l l l   q    
9)  k     k k k k     l     q q q
0)  k   k k k         l     q q q
1)  k k     k k k       r r      
2)  k     k k       r r r   s s  
3)  k k k k   r r r r r     s s  
4)  k   k   t   r   r r     s s  

t=  5 p=0.50 n=   16 sim=0.08984
t=  5 p=0.50 n=   64 sim=0.07310
t=  5 p=0.50 n=  256 sim=0.06706
t=  5 p=0.50 n= 1024 sim=0.06612
t=  5 p=0.50 n= 4096 sim=0.06587

As n increases, the sim result gets closer to 0.065770...

Racket

<lang racket>#lang racket (require srfi/14) ; character sets

much faster than safe fixnum functions

(require

 racket/require ; for fancy require clause below
 (filtered-in
         (lambda (name) (regexp-replace #rx"unsafe-" name ""))
         racket/unsafe/ops)
 ; these aren't in racket/unsafe/ops
 (only-in racket/fixnum for/fxvector in-fxvector fxvector-copy))

...(but less safe). if in doubt use this rather than the one above
(require racket/fixnum)

(define t (make-parameter 5))

(define (build-random-grid p M N)

 (define p-num (numerator p))
 (define p-den (denominator p))
 (for/fxvector #:length (fx* M N) ((_ (in-range (* M N))))
               (if (< (random p-den) p-num) 1 0)))

(define letters

 (sort (char-set->list (char-set-intersection
                        char-set:letter
                        ; char-set:ascii
                        )) char<?))

(define n-letters (length letters)) (define cell->char

 (match-lambda
   (0 #\space) (1 #\.)
   (c (list-ref letters (modulo (- c 2) n-letters)))))

(define (draw-percol-grid M N . gs)

 (for ((r N))
   (for ((g gs))
     (define row-str
       (list->string
        (for/list ((idx (in-range (* r M) (* (+ r 1) M))))
          (cell->char (fxvector-ref g idx)))))
     (printf "|~a| " row-str))
   (newline)))

(define (count-clusters! M N g)

 (define (gather-cluster! k c)
   (when (fx= 1 (fxvector-ref g k))
     (define k-r (fxquotient k M))
     (define k-c (fxremainder k M))
     (fxvector-set! g k c)       
     (define-syntax-rule (gather-surrounds range? k+)
       (let ((idx k+))
         (when (and range? (fx= 1 (fxvector-ref g idx)))
           (gather-cluster! idx c))))
     (gather-surrounds (fx> k-r 0) (fx- k M))
     (gather-surrounds (fx> k-c 0) (fx- k 1))
     (gather-surrounds (fx< k-c (fx- M 1)) (fx+ k 1))
     (gather-surrounds (fx< k-r (fx- N 1)) (fx+ k M))))
 
 (define-values (rv _c)
   (for/fold ((rv 0) (c 2))
     ((pos (in-range (fx* M N)))
      #:when (fx= 1 (fxvector-ref g pos)))
     (gather-cluster! pos c)
     (values (fx+ rv 1) (fx+ c 1))))
 rv)

(define (display-sample-clustering p)

 (printf "Percolation cluster sample: p=~a~%" p)
 (define g (build-random-grid p 15 15))
 (define g+ (fxvector-copy g))
 (define g-count (count-clusters! 15 15 g+))
 (draw-percol-grid 15 15 g g+)
 (printf "~a clusters~%" g-count))

(define (experiment p n t)

 (printf "Experiment: ~a ~a ~a\t" p n t) (flush-output)
 (define sum-Cn
   (for/sum ((run (in-range t)))
     (printf "[~a" run) (flush-output)
     (define g (build-random-grid p n n))
     (printf "*") (flush-output)
     (define Cn (count-clusters! n n g))
     (printf "]") (flush-output)
     Cn))
 (printf "\tmean K(p) = ~a~%" (real->decimal-string (/ sum-Cn t (sqr n)) 6)))

(module+ main

 (t 10)
 (for ((n (in-list '(4000 1000 750 500 400 300 200 100 15))))
   (experiment 1/2 n (t)))
 (display-sample-clustering 1/2))

(module+ test

 (define grd (build-random-grid 1/2 1000 1000))
 (/ (for/sum ((g (in-fxvector grd)) #:when (zero? g)) 1) (fxvector-length grd))
 (display-sample-clustering 1/2))</lang>

Output:

Run from DrRacket, which runs the test and main modules. From the command line, you'll want two commands: ``racket percolation_m_c_d.rkt`` and ``raco test percolation_m_c_d.rkt`` for the same result.

Experiment: 1/2 4000 10	[0*][1*][2*][3*][4*][5*][6*][7*][8*][9*]	mean K(p) = 0.065860
Experiment: 1/2 1000 10	[0*][1*][2*][3*][4*][5*][6*][7*][8*][9*]	mean K(p) = 0.066130
Experiment: 1/2 750 10	[0*][1*][2*][3*][4*][5*][6*][7*][8*][9*]	mean K(p) = 0.066195
Experiment: 1/2 500 10	[0*][1*][2*][3*][4*][5*][6*][7*][8*][9*]	mean K(p) = 0.066522
Experiment: 1/2 400 10	[0*][1*][2*][3*][4*][5*][6*][7*][8*][9*]	mean K(p) = 0.066778
Experiment: 1/2 300 10	[0*][1*][2*][3*][4*][5*][6*][7*][8*][9*]	mean K(p) = 0.066813
Experiment: 1/2 200 10	[0*][1*][2*][3*][4*][5*][6*][7*][8*][9*]	mean K(p) = 0.067908
Experiment: 1/2 100 10	[0*][1*][2*][3*][4*][5*][6*][7*][8*][9*]	mean K(p) = 0.069980
Experiment: 1/2 15 10	[0*][1*][2*][3*][4*][5*][6*][7*][8*][9*]	mean K(p) = 0.089778
Percolation cluster sample: p=1/2
|.  ...     . . | |A  BBB     A A | 
|...    .. .... | |AAA    AA AAAA | 
|. .   .... ... | |A A   AAAA AAA | 
|. . . .........| |A A C AAAAAAAAA| 
| ...   ..  ....| | AAA   AA  AAAA| 
|.. ......... ..| |AA AAAAAAAAA AA| 
|     . ...     | |     A AAA     | 
|. ..  ..       | |D AA  AA       | 
|   .. ... . .. | |   AA AAA E AA | 
|.  ..  ..  . . | |F  AA  AA  A A | 
|. ........ . ..| |F AAAAAAAA A AA| 
|.. .  .... ... | |FF A  AAAA AAA | 
| .  .  . ....  | | F  G  A AAAA  | 
|.... .. ..  . .| |FFFF HH AA  A A| 
|  .  ..   .....| |  F  HH   AAAAA| 
8 clusters