Percolation/Mean cluster density: Difference between revisions
m
→{{header|11l}}: Void
(Updated D entry) |
Alextretyak (talk | contribs) m (→{{header|11l}}: Void) |
||
| (48 intermediate revisions by 19 users not shown) | |||
Line 2:
Let <math>c</math> be a 2D boolean square matrix of <math>n \times n</math> values of either <tt>1</tt> or <tt>0</tt> where the
probability of any value being <tt>1</tt> is <math>p</math>, (and of <tt>0</tt> is therefore <math>1-p</math>).
We define a ''cluster'' of <tt>1</tt>'s as being a group of <tt>1</tt>'s connected vertically or
horizontally (i.e., using the [[wp:Von Neumann neighborhood|Von Neumann neighborhood rule]]) and bounded by either <math>0</math> or by the limits of the matrix.
Let the number of such clusters in such a randomly constructed matrix be <math>C_n</math>.
Percolation theory states that <math>K(p)</math> (the mean cluster density) will satisfy <math>K(p) = C_n / n^2</math>
;Task
Show the effect of varying <math>n</math> on the accuracy of simulated <math>K(p)</math> for <math>p = 0.5</math> and
for values of <math>n</math> up to at least <math>1000</math>.
Any calculation of <math>C_n</math> for finite <math>n</math> is subject to randomness, so an approximation should be
computed as the average of <math>t</math> runs, where <math>t</math> ≥ <math>5</math>.
For extra credit, graphically show clusters in a <math>15\times 15</math>, <math>p=0.5</math> grid.
Line 25 ⟶ 21:
;See also
* [http://mathworld.wolfram.com/s-Cluster.html s-Cluster] on Wolfram mathworld.
=={{header|11l}}==
{{trans|Nim}}
<syntaxhighlight lang="11l">UInt32 seed = 0
F nonrandom()
:seed = 1664525 * :seed + 1013904223
R (:seed >> 16) / Float(FF'FF)
V nn = 15
V tt = 5
V pp = 0.5
V NotClustered = 1
V Cell2Char = ‘ #abcdefghijklmnopqrstuvwxyz’
V NRange = [4, 64, 256, 1024, 4096]
F newGrid(n, p)
R (0 .< n).map(i -> (0 .< @n).map(i -> Int(nonrandom() < @@p)))
F walkMaze(&grid, m, n, idx) -> Void
grid[n][m] = idx
I n < grid.len - 1 & grid[n + 1][m] == NotClustered
walkMaze(&grid, m, n + 1, idx)
I m < grid[0].len - 1 & grid[n][m + 1] == NotClustered
walkMaze(&grid, m + 1, n, idx)
I m > 0 & grid[n][m - 1] == NotClustered
walkMaze(&grid, m - 1, n, idx)
I n > 0 & grid[n - 1][m] == NotClustered
walkMaze(&grid, m, n - 1, idx)
F clusterCount(&grid)
V walkIndex = 1
L(n) 0 .< grid.len
L(m) 0 .< grid[0].len
I grid[n][m] == NotClustered
walkIndex++
walkMaze(&grid, m, n, walkIndex)
R walkIndex - 1
F clusterDensity(n, p)
V grid = newGrid(n, p)
R clusterCount(&grid) / Float(n * n)
F print_grid(grid)
L(row) grid
print(L.index % 10, end' ‘) ’)
L(cell) row
print(‘ ’Cell2Char[cell], end' ‘’)
print()
V grid = newGrid(nn, 0.5)
print(‘Found ’clusterCount(&grid)‘ clusters in this ’nn‘ by ’nn" grid\n")
print_grid(grid)
print()
L(n) NRange
V sum = 0.0
L 0 .< tt
sum += clusterDensity(n, pp)
V sim = sum / tt
print(‘t = #. p = #.2 n = #4 sim = #.5’.format(tt, pp, n, sim))</syntaxhighlight>
{{out}}
<pre>
Found 25 clusters in this 15 by 15 grid
0) a a b c d
1) e e d d d d d d
2) e e e e d d d d
3) e e e e e e e e d d d d
4) e e e e e e e e d d d
5) e e e e e f d
6) g e h e i d
7) g j k k d d
8) l m k k k k k
9) n l m o k k k k p
0) n k k k k k q
1) n r r s k t u
2) r k k k u
3) v r r w k k k x
4) v r r w w w y
t = 5 p = 0.50 n = 4 sim = 0.17500
t = 5 p = 0.50 n = 64 sim = 0.07300
t = 5 p = 0.50 n = 256 sim = 0.06823
t = 5 p = 0.50 n = 1024 sim = 0.06618
t = 5 p = 0.50 n = 4096 sim = 0.06590
</pre>
=={{header|C}}==
<
#include <stdlib.h>
Line 103 ⟶ 187:
free(map);
return 0;
}</
{{out}}
<pre>
Line 131 ⟶ 215:
4096 0.065836
16384 0.065774
</pre>
=={{header|C++}}==
<syntaxhighlight lang="c++">
#include <iostream>
#include <random>
#include <string>
#include <vector>
#include <iomanip>
std::random_device random;
std::mt19937 generator(random());
std::uniform_real_distribution<double> distribution(0.0F, 1.0F);
class Grid {
public:
Grid(const int32_t size, const double probability) {
create_grid(size, probability);
count_clusters();
}
int32_t cluster_count() const {
return clusters;
}
double cluster_density() const {
return (double) clusters / ( grid.size() * grid.size() );
}
void display() const {
for ( uint64_t row = 0; row < grid.size(); ++row ) {
for ( uint64_t col = 0; col < grid.size(); ++col ) {
uint64_t value = grid[row][col];
char ch = ( value < GRID_CHARACTERS.length() ) ? GRID_CHARACTERS[value] : '?';
std::cout << " " << ch;
}
std::cout << std::endl;
}
}
private:
void count_clusters() {
clusters = 0;
for ( uint64_t row = 0; row < grid.size(); ++row ) {
for ( uint64_t col = 0; col < grid.size(); ++col ) {
if ( grid[row][col] == CLUSTERED ) {
clusters += 1;
identify_cluster(row, col, clusters);
}
}
}
}
void identify_cluster(const uint64_t row, const uint64_t col, const uint64_t count) {
grid[row][col] = count;
if ( row < grid.size() - 1 && grid[row + 1][col] == CLUSTERED ) {
identify_cluster(row + 1, col, count);
}
if ( col < grid.size() - 1 && grid[row][col + 1] == CLUSTERED ) {
identify_cluster(row, col + 1, count);
}
if ( col > 0 && grid[row][col - 1] == CLUSTERED ) {
identify_cluster(row, col - 1, count);
}
if ( row > 0 && grid[row - 1][col] == CLUSTERED ) {
identify_cluster(row - 1, col, count);
}
}
void create_grid(int32_t grid_size, double probability) {
grid.assign(grid_size, std::vector<int32_t>(grid_size, 0));
for ( int32_t row = 0; row < grid_size; ++row ) {
for ( int32_t col = 0; col < grid_size; ++col ) {
if ( distribution(generator) < probability ) {
grid[row][col] = CLUSTERED;
}
}
}
}
int32_t clusters;
std::vector<std::vector<int32_t>> grid;
inline static const int CLUSTERED = -1;
inline static const std::string GRID_CHARACTERS = ".ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
};
int main() {
const int32_t size = 15;
const double probability = 0.5;
const int32_t test_count = 5;
Grid grid(size, probability);
std::cout << "This " << size << " by " << size << " grid contains "
<< grid.cluster_count() << " clusters:" << std::endl;
grid.display();
std::cout << "\n p = 0.5, iterations = " << test_count << std::endl;
std::vector<int32_t> grid_sizes = { 10, 100, 1'000, 10'000 };
for ( int32_t grid_size : grid_sizes ) {
double sumDensity = 0.0;
for ( int32_t test = 0; test < test_count; test++ ) {
Grid grid(grid_size, probability);
sumDensity += grid.cluster_density();
}
double result = sumDensity / test_count;
std::cout << " n = " << std::setw(5) << grid_size
<< ", simulations K = " << std::fixed << result << std::endl;
}
}
</syntaxhighlight>
{{ out }}
<pre>
This 15 by 15 grid contains 11 clusters:
A A . B B . B B . . . . . . .
. . B B B B . B B . B B . . .
. C . . B B . B B B B . . . .
. C C . B . B B B B B B . D D
C C . . B B B B . . B . D D .
C C . . B B B . B B B . . D D
C C . . . . B . . . B . D D D
. C . E . . . . D D . . D D .
F . . . G . H . D D D D D D D
F F . F . I . F . D D . D D .
F . F F F . . F . . D D . . D
F . F . F . . F . F . D D D D
F F . . F F F F F F . D . . D
F F F F F . . F . . . D . . D
. F F . . J J . . . . D . K .
p = 0.5, iterations = 5
n = 10, simulation K = 0.088000
n = 100, simulation K = 0.067260
n = 1000, simulation K = 0.066215
n = 10000, simulation K = 0.065777
</pre>
=={{header|D}}==
{{trans|python}}
<
std.range, std.ascii;
Line 142 ⟶ 361:
enum Cell notClustered = 1; // Filled cell, but not in a cluster.
Grid initialize(Grid grid, in double prob, ref Xorshift rng) nothrow {
foreach (row; grid)
foreach (ref cell; row)
cell =
return grid;
}
Line 168 ⟶ 379:
}
size_t countClusters(bool justCount=false)(Grid grid)
pure nothrow @safe @nogc {
immutable side = grid.length;
static if (justCount)
Line 175 ⟶ 387:
Cell clusterID = 1;
void walk(in size_t r, in size_t c) nothrow @safe @nogc {
grid[r][c] = clusterID; // Fill grid.
Line 231 ⟶ 443:
nIters, prob, side, density);
}
}</
{{out}}
<pre>Found 26 clusters in this 15 by 15 grid:
Line 260 ⟶ 472:
Increasing the index i to 15:
<pre>n_iters= 5, p=0.50, n=32768, sim=0.06578374</pre>
=={{header|EchoLisp}}==
We use the canvas bit-map as 2D-matrix. For extra-extra credit, a 800x800 nice cluster tapestry image is shown here : http://www.echolalie.org/echolisp/images/rosetta-clusters-800.png.
<syntaxhighlight lang="scheme">
(define-constant BLACK (rgb 0 0 0.6))
(define-constant WHITE -1)
;; sets pixels to clusterize to WHITE
;; returns bit-map vector
(define (init-C n p )
(plot-size n n)
(define C (pixels->int32-vector )) ;; get canvas bit-map
(pixels-map (lambda (x y) (if (< (random) p) WHITE BLACK )) C)
C )
;; random color for new cluster
(define (new-color)
(hsv->rgb (random) 0.9 0.9))
;; make-region predicate
(define (in-cluster C x y)
(= (pixel-ref C x y) WHITE))
;; paint all adjacents to (x0,y0) with new color
(define (make-cluster C x0 y0)
(pixel-set! C x0 y0 (new-color))
(make-region in-cluster C x0 y0))
;; task
(define (make-clusters (n 400) (p 0.5))
(define Cn 0)
(define C null)
(for ((t 5)) ;; 5 iterations
(plot-clear)
(set! C (init-C n p))
(for* ((x0 n) (y0 n))
#:when (= (pixel-ref C x0 y0) WHITE)
(set! Cn (1+ Cn))
(make-cluster C x0 y0)))
(writeln 'n n 'Cn Cn 'density (// Cn (* n n) 5) )
(vector->pixels C)) ;; to screen
</syntaxhighlight>
{{out}}
<pre>
n 100 Cn 3420 density 0.0684
n 400 Cn 53246 density 0.0665575
n 600 Cn 118346 density 0.06574778
n 800 Cn 212081 density 0.0662753125
n 1000 Cn 330732 density 0.0661464
</pre>
=={{header|Factor}}==
<syntaxhighlight lang="factor">USING: combinators formatting generalizations kernel math
math.matrices random sequences ;
IN: rosetta-code.mean-cluster-density
CONSTANT: p 0.5
CONSTANT: iterations 5
: rand-bit-matrix ( n probability -- matrix )
dupd [ random-unit > 1 0 ? ] curry make-matrix ;
: flood-fill ( x y matrix -- )
3dup ?nth ?nth 1 = [
[ [ -1 ] 3dip nth set-nth ] [
{
[ [ 1 + ] 2dip ]
[ [ 1 - ] 2dip ]
[ [ 1 + ] dip ]
[ [ 1 - ] dip ]
} [ flood-fill ] map-compose 3cleave
] 3bi
] [ 3drop ] if ;
: count-clusters ( matrix -- Cn )
0 swap dup dim matrix-coordinates flip concat [
first2 rot 3dup ?nth ?nth 1 = [ flood-fill 1 + ]
[ 3drop ] if
] with each ;
: mean-cluster-density ( matrix -- mcd )
[ count-clusters ] [ dim first sq / ] bi ;
: simulate ( n -- avg-mcd )
iterations swap [ p rand-bit-matrix mean-cluster-density ]
curry replicate sum iterations / ;
: main ( -- )
{ 4 64 256 1024 4096 } [
[ iterations p ] dip dup simulate
"iterations = %d p = %.1f n = %4d sim = %.5f\n" printf
] each ;
MAIN: main</syntaxhighlight>
{{out}}
<pre>
iterations = 5 p = 0.5 n = 4 sim = 0.13750
iterations = 5 p = 0.5 n = 64 sim = 0.07437
iterations = 5 p = 0.5 n = 256 sim = 0.06786
iterations = 5 p = 0.5 n = 1024 sim = 0.06621
iterations = 5 p = 0.5 n = 4096 sim = 0.06589
</pre>
=={{header|Go}}==
{{trans|Python}}
<syntaxhighlight lang="go">package main
import (
"fmt"
"math/rand"
"time"
)
var (
n_range = []int{4, 64, 256, 1024, 4096}
M = 15
N = 15
)
const (
p = .5
t = 5
NOT_CLUSTERED = 1
cell2char = " #abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"
)
func newgrid(n int, p float64) [][]int {
g := make([][]int, n)
for y := range g {
gy := make([]int, n)
for x := range gy {
if rand.Float64() < p {
gy[x] = 1
}
}
g[y] = gy
}
return g
}
func pgrid(cell [][]int) {
for n := 0; n < N; n++ {
fmt.Print(n%10, ") ")
for m := 0; m < M; m++ {
fmt.Printf(" %c", cell2char[cell[n][m]])
}
fmt.Println()
}
}
func cluster_density(n int, p float64) float64 {
cc := clustercount(newgrid(n, p))
return float64(cc) / float64(n) / float64(n)
}
func clustercount(cell [][]int) int {
walk_index := 1
for n := 0; n < N; n++ {
for m := 0; m < M; m++ {
if cell[n][m] == NOT_CLUSTERED {
walk_index++
walk_maze(m, n, cell, walk_index)
}
}
}
return walk_index - 1
}
func walk_maze(m, n int, cell [][]int, indx int) {
cell[n][m] = indx
if n < N-1 && cell[n+1][m] == NOT_CLUSTERED {
walk_maze(m, n+1, cell, indx)
}
if m < M-1 && cell[n][m+1] == NOT_CLUSTERED {
walk_maze(m+1, n, cell, indx)
}
if m > 0 && cell[n][m-1] == NOT_CLUSTERED {
walk_maze(m-1, n, cell, indx)
}
if n > 0 && cell[n-1][m] == NOT_CLUSTERED {
walk_maze(m, n-1, cell, indx)
}
}
func main() {
rand.Seed(time.Now().Unix())
cell := newgrid(N, .5)
fmt.Printf("Found %d clusters in this %d by %d grid\n\n",
clustercount(cell), N, N)
pgrid(cell)
fmt.Println()
for _, n := range n_range {
M = n
N = n
sum := 0.
for i := 0; i < t; i++ {
sum += cluster_density(n, p)
}
sim := sum / float64(t)
fmt.Printf("t=%3d p=%4.2f n=%5d sim=%7.5f\n", t, p, n, sim)
}
}</syntaxhighlight>
{{out}}
<pre>
Found 29 clusters in this 15 by 15 grid
0) a a a b c
1) d e a c c f
2) g e e e
3) h e i j k
4) l h m m n
5) l o p n n n n
6) q o o o r n
7) o o r r s
8) t t o u r s s s s
9) t t v u u r r s s s s
0) t v u r r s
1) w x r r y y
2) z x r r r r r y
3) A z z B r r r r
4) A A A z C r r r r
t= 5 p=0.50 n= 4 sim=0.16250
t= 5 p=0.50 n= 64 sim=0.07334
t= 5 p=0.50 n= 256 sim=0.06710
t= 5 p=0.50 n= 1024 sim=0.06619
t= 5 p=0.50 n= 4096 sim=0.06585
</pre>
=={{header|Haskell}}==
<syntaxhighlight lang="haskell">{-# language FlexibleContexts #-}
import Data.List
import Data.Maybe
import System.Random
import Control.Monad.State
import Text.Printf
import Data.Set (Set)
import qualified Data.Set as S
type Matrix = [[Bool]]
type Cell = (Int, Int)
type Cluster = Set (Int, Int)
clusters :: Matrix -> [Cluster]
clusters m = unfoldr findCuster cells
where
cells = S.fromList [ (i,j) | (r, i) <- zip m [0..]
, (x, j) <- zip r [0..], x]
findCuster s = do
(p, ps) <- S.minView s
return $ runState (expand p) ps
expand p = do
ns <- state $ extract (neigbours p)
xs <- mapM expand $ S.elems ns
return $ S.insert p $ mconcat xs
extract s1 s2 = (s2 `S.intersection` s1, s2 S.\\ s1)
neigbours (i,j) = S.fromList [(i-1,j),(i+1,j),(i,j-1),(i,j+1)]
n = length m
showClusters :: Matrix -> String
showClusters m = unlines [ unwords [ mark (i,j)
| j <- [0..n-1] ]
| i <- [0..n-1] ]
where
cls = clusters m
n = length m
mark c = maybe "." snd $ find (S.member c . fst) $ zip cls syms
syms = sequence [['a'..'z'] ++ ['A'..'Z']]
------------------------------------------------------------
randomMatrices :: Int -> StdGen -> [Matrix]
randomMatrices n = clipBy n . clipBy n . randoms
where
clipBy n = unfoldr (Just . splitAt n)
randomMatrix n = head . randomMatrices n
tests :: Int -> StdGen -> [Int]
tests n = map (length . clusters) . randomMatrices n
task :: Int -> StdGen -> (Int, Double)
task n g = (n, result)
where
result = mean $ take 10 $ map density $ tests n g
density c = fromIntegral c / fromIntegral n**2
mean lst = sum lst / genericLength lst
main = newStdGen >>= mapM_ (uncurry (printf "%d\t%.5f\n")) . res
where
res = mapM task [10,50,100,500]</syntaxhighlight>
<pre>λ> newStdGen >>= putStrLn . showClusters . randomMatrix 15
. . a a . b b b b . . c c . .
d d . . . . . . b b . c . . .
d . . e . . . b b . c c . f f
d d d . g g . b b . c c c . .
. d . d . . b b . h . . c . i
d d . d d d . . h h h h . . i
d d d d . d . . h . . . . i i
. . . d . d . . h . i i i i i
. j . d . . . . . k . i . i .
. . l . . . . k k k k . . i i
m m . m . . . k k . . n . i .
m m m m m . o . k . n n . . .
. m . m . p . k k . . n . . q
. m m . . . r . k . . n n . q
. m . s s . r r . t . . . . q
λ> take 10 $ tests 15 (mkStdGen 42)
[33,18,26,18,29,14,23,21,18,24]
λ> main
10 0.10100
50 0.07072
100 0.06878
500 0.06676</pre>
=={{header|J}}==
The first thing this task seems to need is some mechanism of identifying "clusters", using "percolation". We can achieve this by assigning every "1" in a matrix a unique integer value and then defining an operation which combines two numbers - doing nothing unless the second one (the one on the right) is non-zero. If it is non-zero we pick the larger of the two values. <code>(*@[ * >.)</code>
Once we have this, we can identify clusters by propagating information in a single direction through the matrix using this operation, rotating the matrix 90 degrees, and then repeating this combination of operations four times. And, finally, by keeping at this until there's nothing more to be done.
<syntaxhighlight lang="j">congeal=: |.@|:@((*@[*>.)/\.)^:4^:_</syntaxhighlight>
Example:
<syntaxhighlight lang="j"> M=:0.4>?6 6$0
M
1 0 0 0 0 0
0 0 0 1 0 0
0 0 0 1 0 0
1 1 0 0 0 0
0 0 0 1 0 1
1 1 0 1 1 0
M*p:i.$M
2 0 0 0 0 0
0 0 0 29 0 0
0 0 0 53 0 0
67 71 0 0 0 0
0 0 0 107 0 113
127 131 0 139 149 0
congeal M*p:i.$M
2 0 0 0 0 0
0 0 0 53 0 0
0 0 0 53 0 0
71 71 0 0 0 0
0 0 0 149 0 113
131 131 0 149 149 0</syntaxhighlight>
We did not have to use primes there - any mechanism for assigning distinct positive integers to the 1s would work. And, in fact, it might be nice if - once we found our clusters - we assigned the smallest distinct positive integers to the clusters. This would allow us to use simple indexing to map the array to characters.
<syntaxhighlight lang="j">idclust=: $ $ [: (~. i.])&.(0&,)@,@congeal ] * 1 + i.@$</syntaxhighlight>
Example use:
<syntaxhighlight lang="j"> idclust M
1 0 0 0 0 0
0 0 0 2 0 0
0 0 0 2 0 0
3 3 0 0 0 0
0 0 0 4 0 5
6 6 0 4 4 0
(idclust M) {'.ABCDEFG'
A.....
...B..
...B..
CC....
...D.E
FF.DD.</syntaxhighlight>
Now we just need a measure of cluster density. Formally cluster density seems to be defined as the number of clusters divided by the total number of elements of the matrix. Thus:
<syntaxhighlight lang="j">K=: (%&#~ }.@~.)&,</syntaxhighlight>
Example use:
<syntaxhighlight lang="j"> K idclust M
0.1666667</syntaxhighlight>
So we can create a word that performs a simulation experiment, given a probability getting a 1 and the number of rows (and columns) of our square matrix M.
<syntaxhighlight lang="j">experiment=: K@ idclust@: > 0 ?@$~ ,~</syntaxhighlight>
Example use:
<syntaxhighlight lang="j"> 0.4 experiment 6
0.1666667
0.4 experiment 6
0.1944444</syntaxhighlight>
The task wants us to perform at least five trials for sizes up to 1000 by 1000 with probability of 1 being 0.5:
<syntaxhighlight lang="j">trials=: 0.5&experiment"0@#</syntaxhighlight>
Example use:
<syntaxhighlight lang="j"> 6 trials 3
0.1111111 0.1111111 0.2222222 0.1111111 0.1111111 0.3333333
6 trials 10
0.16 0.12 0.09 0.1 0.1 0.03
6 trials 30
0.05666667 0.1033333 0.08222222 0.07444444 0.08333333 0.07666667
6 trials 100
0.069 0.0678 0.0666 0.0677 0.0653 0.0739
6 trials 300
0.06563333 0.06663333 0.06713333 0.06727778 0.06658889 0.06664444
6 trials 1000
0.066079 0.066492 0.065847 0.065943 0.066318 0.065998</syntaxhighlight>
Now for averages (these are different trials from the above):
<syntaxhighlight lang="j">mean=: +/%#
mean 8 trials 3
0.1805556
mean 8 trials 10
0.0875
mean 8 trials 30
0.07486111
mean 8 trials 100
0.0690625
mean 8 trials 300
0.06749861
mean 8 trials 1000
0.06616738</syntaxhighlight>
Finally, for the extra credit (thru taken from the [[Loops/Downward_for#J|Loops/Downward for]] task):
<syntaxhighlight lang="j">thru=: <./ + i.@(+*)@-~</syntaxhighlight>
<syntaxhighlight lang="j"> (idclust 0.5 > 0 ?@$~ 15 15) {'.', 'A' thru&.(a.&i.) 'Z'
A.......B..C...
AAAA...D..E.F..
A..A.G.D.D.FFF.
AA..H..DDD.FF.I
AAA...J...FFF..
..AAAA.A.K...AA
LL.A...A..A.AAA
.L.A..AAA.AAAAA
..AA.AAA.AAA.A.
AA.AAAAAA....A.
A.AAAA.AAAA.AA.
AAA...AAA.AAAAA
..AA..A.A...AAA
.M.A.AA.AA..AA.
.MM..A.N..O..A.</syntaxhighlight>
'''Collected definitions'''
<syntaxhighlight lang="j">congeal=: |.@|:@((*@[*>.)/\.)^:4^:_
idclust=: $ $ [: (~. i.])&.(0&,)@,@congeal ] * 1 + i.@$
K=: (%&#~ }.@~.)&,
experiment=: K@ idclust@: > 0 ?@$~ ,~
trials=: 0.5&experiment"0@#
mean=:+/ % #
thru=: <./ + i.@(+*)@-~</syntaxhighlight>
'''Extra Credit'''
<syntaxhighlight lang="j"> M=: (* 1+i.@$)?15 15$2
M
0 2 3 4 0 6 0 8 0 10 11 12 0 0 15
0 0 18 19 20 0 22 0 0 0 0 0 28 29 0
31 32 0 34 35 36 37 38 0 0 0 42 0 0 45
0 0 48 49 0 51 0 0 54 55 0 57 58 0 0
61 62 63 64 0 0 67 0 69 0 71 72 0 74 0
0 0 78 79 0 0 82 0 84 85 86 87 88 0 0
0 92 0 94 0 0 0 0 99 100 101 0 103 0 105
106 107 108 0 0 111 0 0 114 115 116 0 0 0 0
0 0 0 124 125 126 127 0 0 0 0 0 133 134 135
0 0 138 0 0 141 0 143 144 145 0 0 0 0 150
0 152 153 154 0 0 0 158 0 160 0 162 163 164 165
0 167 168 169 170 0 172 173 0 175 176 177 0 0 180
181 182 183 0 0 186 0 188 189 190 191 192 0 194 195
196 197 198 0 200 201 202 0 0 205 0 207 0 0 0
211 212 213 0 0 0 217 218 0 220 221 0 0 224 0
congeal M
0 94 94 94 0 6 0 8 0 12 12 12 0 0 15
0 0 94 94 94 0 94 0 0 0 0 0 29 29 0
32 32 0 94 94 94 94 94 0 0 0 116 0 0 45
0 0 94 94 0 94 0 0 116 116 0 116 116 0 0
94 94 94 94 0 0 82 0 116 0 116 116 0 74 0
0 0 94 94 0 0 82 0 116 116 116 116 116 0 0
0 108 0 94 0 0 0 0 116 116 116 0 116 0 105
108 108 108 0 0 141 0 0 116 116 116 0 0 0 0
0 0 0 141 141 141 141 0 0 0 0 0 221 221 221
0 0 213 0 0 141 0 221 221 221 0 0 0 0 221
0 213 213 213 0 0 0 221 0 221 0 221 221 221 221
0 213 213 213 213 0 221 221 0 221 221 221 0 0 221
213 213 213 0 0 218 0 221 221 221 221 221 0 221 221
213 213 213 0 218 218 218 0 0 221 0 221 0 0 0
213 213 213 0 0 0 218 218 0 221 221 0 0 224 0
(~.@, i. ])congeal M
0 1 1 1 0 2 0 3 0 4 4 4 0 0 5
0 0 1 1 1 0 1 0 0 0 0 0 6 6 0
7 7 0 1 1 1 1 1 0 0 0 8 0 0 9
0 0 1 1 0 1 0 0 8 8 0 8 8 0 0
1 1 1 1 0 0 10 0 8 0 8 8 0 11 0
0 0 1 1 0 0 10 0 8 8 8 8 8 0 0
0 12 0 1 0 0 0 0 8 8 8 0 8 0 13
12 12 12 0 0 14 0 0 8 8 8 0 0 0 0
0 0 0 14 14 14 14 0 0 0 0 0 15 15 15
0 0 16 0 0 14 0 15 15 15 0 0 0 0 15
0 16 16 16 0 0 0 15 0 15 0 15 15 15 15
0 16 16 16 16 0 15 15 0 15 15 15 0 0 15
16 16 16 0 0 17 0 15 15 15 15 15 0 15 15
16 16 16 0 17 17 17 0 0 15 0 15 0 0 0
16 16 16 0 0 0 17 17 0 15 15 0 0 18 0</syntaxhighlight>
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.util.List;
import java.util.concurrent.ThreadLocalRandom;
public final class PercolationMeanCluster {
public static void main(String[] aArgs) {
final int size = 15;
final double probability = 0.5;
final int testCount = 5;
Grid grid = new Grid(size, probability);
System.out.println("This " + size + " by " + size + " grid contains " + grid.clusterCount() + " clusters:");
grid.display();
System.out.println(System.lineSeparator() + " p = 0.5, iterations = " + testCount);
List<Integer> gridSizes = List.of( 10, 100, 1_000, 10_000 );
for ( int gridSize : gridSizes ) {
double sumDensity = 0.0;
for ( int test = 0; test < testCount; test++ ) {
grid = new Grid(gridSize, probability);
sumDensity += grid.clusterDensity();
}
double result = sumDensity / testCount;
System.out.println(String.format("%s%5d%s%.6f", " n = ", gridSize, ", simulation K = ", result));
}
}
}
final class Grid {
public Grid(int aSize, double aProbability) {
createGrid(aSize, aProbability);
countClusters();
}
public int clusterCount() {
return clusterCount;
}
public double clusterDensity() {
return (double) clusterCount / ( grid.length * grid.length );
}
public void display() {
for ( int row = 0; row < grid.length; row++ ) {
for ( int col = 0; col < grid.length; col++ ) {
int value = grid[row][col];
char ch = ( value < GRID_CHARACTERS.length() ) ? GRID_CHARACTERS.charAt(value) : '?';
System.out.print(" " + ch);
}
System.out.println();
}
}
private void countClusters() {
clusterCount = 0;
for ( int row = 0; row < grid.length; row++ ) {
for ( int col = 0; col < grid.length; col++ ) {
if ( grid[row][col] == CLUSTERED ) {
clusterCount += 1;
identifyCluster(row, col, clusterCount);
}
}
}
}
private void identifyCluster(int aRow, int aCol, int aCount) {
grid[aRow][aCol] = aCount;
if ( aRow < grid.length - 1 && grid[aRow + 1][aCol] == CLUSTERED ) {
identifyCluster(aRow + 1, aCol, aCount);
}
if ( aCol < grid[0].length - 1 && grid[aRow][aCol + 1] == CLUSTERED ) {
identifyCluster(aRow, aCol + 1, aCount);
}
if ( aCol > 0 && grid[aRow][aCol - 1] == CLUSTERED ) {
identifyCluster(aRow, aCol - 1, aCount);
}
if ( aRow > 0 && grid[aRow - 1][aCol] == CLUSTERED ) {
identifyCluster(aRow - 1, aCol, aCount);
}
}
private void createGrid(int aGridSize, double aProbability) {
grid = new int[aGridSize][aGridSize];
for ( int row = 0; row < aGridSize; row++ ) {
for ( int col = 0; col < aGridSize; col++ ) {
if ( random.nextDouble(1.0) < aProbability ) {
grid[row][col] = CLUSTERED;
}
}
}
}
private int[][] grid;
private int clusterCount;
private static ThreadLocalRandom random = ThreadLocalRandom.current();
private static final int CLUSTERED = -1;
private static final String GRID_CHARACTERS = ".ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
}
</syntaxhighlight>
{{ out }}
<pre>
This 15 by 15 grid contains 21 clusters:
A . . B B B . . C . D D . E .
. . F . B . G . C . . . E E E
. F F . B . . . C C . H . E .
F F . B B B . I . C C . . E .
. . . . . . . I . . . J . . K
. L . . M M . I . . N . . . K
O . . O . M . I I . . . . K K
O O O O O . I I . K K . . . K
. O . O . P . I . K . K K K K
. . Q . . . I I . K K K K K K
Q Q Q . . I I I I . . . . . K
Q . Q Q . . I . . . . R R . .
. . Q . I I I . . . . R . R R
. . Q . . I I . . . . R R R .
S S . T . . I I I . R R R . U
p = 0.5, iterations = 5
n = 10, simulation K = 0.094000
n = 100, simulation K = 0.070420
n = 1000, simulation K = 0.066056
n = 10000, simulation K = 0.065780
</pre>
=={{header|Julia}}==
{{trans|Python}}
<syntaxhighlight lang="julia">using Printf, Distributions
newgrid(p::Float64, r::Int, c::Int=r) = rand(Bernoulli(p), r, c)
function walkmaze!(grid::Matrix{Int}, r::Int, c::Int, indx::Int)
NOT_CLUSTERED = 1 # const
N, M = size(grid)
dirs = [[1, 0], [-1, 0], [0, 1], [0, -1]]
# fill cell
grid[r, c] = indx
# check for each direction
for d in dirs
rr, cc = (r, c) .+ d
if checkbounds(Bool, grid, rr, cc) && grid[rr, cc] == NOT_CLUSTERED
walkmaze!(grid, rr, cc, indx)
end
end
end
function clustercount!(grid::Matrix{Int})
NOT_CLUSTERED = 1 # const
walkind = 1
for r in 1:size(grid, 1), c in 1:size(grid, 2)
if grid[r, c] == NOT_CLUSTERED
walkind += 1
walkmaze!(grid, r, c, walkind)
end
end
return walkind - 1
end
clusterdensity(p::Float64, n::Int) = clustercount!(newgrid(p, n)) / n ^ 2
function printgrid(G::Matrix{Int})
LETTERS = vcat(' ', '#', 'A':'Z', 'a':'z')
for r in 1:size(G, 1)
println(r % 10, ") ", join(LETTERS[G[r, :] .+ 1], ' '))
end
end
G = newgrid(0.5, 15)
@printf("Found %i clusters in this %i×%i grid\n\n", clustercount!(G), size(G, 1), size(G, 2))
printgrid(G)
println()
const nrange = 2 .^ (4:2:12)
const p = 0.5
const nrep = 5
for n in nrange
sim = mean(clusterdensity(p, n) for _ in 1:nrep)
@printf("nrep = %2i p = %.2f dim = %-13s sim = %.5f\n", nrep, p, "$n × $n", sim)
end</syntaxhighlight>
{{out}}
<pre>Found 20 clusters in this 15×15 grid
1) A B C C D D D
2) E F D D
3) G F D D D D D D
4) G G H F I D D D J
5) G G K L D J
6) G G G G M N N
7) G G G G G G O O O N
8) G G O N N N
9) P P P G G G N N
0) P P P P G Q Q Q N
1) P Q Q Q Q N N
2) P N N N N
3) P P P R S N
4) P R R R S S N N
5) R R R S T N
nrep = 5 p = 0.50 dim = 16 × 16 sim = 0.07500
nrep = 5 p = 0.50 dim = 64 × 64 sim = 0.07178
nrep = 5 p = 0.50 dim = 256 × 256 sim = 0.06690
nrep = 5 p = 0.50 dim = 1024 × 1024 sim = 0.06609
nrep = 5 p = 0.50 dim = 4096 × 4096 sim = 0.06588</pre>
=={{header|Kotlin}}==
{{trans|C}}
<syntaxhighlight lang="scala">// version 1.2.10
import java.util.Random
val rand = Random()
const val RAND_MAX = 32767
lateinit var map: IntArray
var w = 0
var ww = 0
const val ALPHA = "+.ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"
const val ALEN = ALPHA.length - 3
fun makeMap(p: Double) {
val thresh = (p * RAND_MAX).toInt()
ww = w * w
var i = ww
map = IntArray(i)
while (i-- != 0) {
val r = rand.nextInt(RAND_MAX + 1)
if (r < thresh) map[i] = -1
}
}
fun showCluster() {
var k = 0
for (i in 0 until w) {
for (j in 0 until w) {
val s = map[k++]
val c = if (s < ALEN) ALPHA[1 + s] else '?'
print(" $c")
}
println()
}
}
fun recur(x: Int, v: Int) {
if ((x in 0 until ww) && map[x] == -1) {
map[x] = v
recur(x - w, v)
recur(x - 1, v)
recur(x + 1, v)
recur(x + w, v)
}
}
fun countClusters(): Int {
var cls = 0
for (i in 0 until ww) {
if (map[i] != -1) continue
recur(i, ++cls)
}
return cls
}
fun tests(n: Int, p: Double): Double {
var k = 0.0
for (i in 0 until n) {
makeMap(p)
k += countClusters().toDouble() / ww
}
return k / n
}
fun main(args: Array<String>) {
w = 15
makeMap(0.5)
val cls = countClusters()
println("width = 15, p = 0.5, $cls clusters:")
showCluster()
println("\np = 0.5, iter = 5:")
w = 1 shl 2
while (w <= 1 shl 13) {
val t = tests(5, 0.5)
println("%5d %9.6f".format(w, t))
w = w shl 1
}
}</syntaxhighlight>
Sample output:
<pre>
width = 15, p = 0.5, 23 clusters:
A . B . C C . . . . D D . D .
. B B B . . D D D . D D . D D
. B . B B B . . D D D . . D .
. . E . . B . B . . D D D D .
F . . B B B B B . G . D . D D
. . B B . B . . H . D D . D D
. B B . . . I . H H . . D D .
B B B . . J . K . . L . D . B
B . . M . . K K K . . D D . B
B B . . . . . K K . . D . B B
B . . N N . . . . . . . O . .
B . . . . P . . O O O O O . .
. . Q . . P . . . O O O O . .
. . . R . . . . S . O . . T T
. U . . . V . . . . . W . . .
p = 0.5, iter = 5:
4 0.112500
8 0.121875
16 0.075000
32 0.068750
64 0.068164
128 0.065625
256 0.067093
512 0.065815
1024 0.065863
2048 0.065815
4096 0.065764
8192 0.065766
</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">(*Calculate C_n / n^2 for n=1000, 2000, ..., 10 000*)
In[1]:= Table[N[Max@MorphologicalComponents[
RandomVariate[BernoulliDistribution[.5], {n, n}],
CornerNeighbors -> False]/n^2], {n, 10^3, 10^4, 10^3}]
(*Find the average*)
In[2]:= % // MeanAround
(*Show a 15x15 matrix with each cluster given an incrementally higher number, Colorize instead of MatrixForm creates an image*)
In[3]:= MorphologicalComponents[RandomChoice[{0, 1}, {15, 15}], CornerNeighbors -> False] // MatrixForm</syntaxhighlight>
{{out}}
<pre>Out[1]= {0.066339, 0.06568, 0.0656282, 0.0658778, 0.0657444, 0.0658455, 0.06578, 0.0657307, 0.0658186, 0.0657963}
Out[2]= 0.06582 +- 0.00006
Out[3]= (1 1 0 2 2 2 2 2 0 0 2 2 2 2 0
0 0 3 0 0 0 2 2 0 0 0 2 2 2 0
3 3 3 3 3 0 0 2 2 2 2 2 2 0 0
0 0 0 0 0 4 4 0 2 0 0 2 0 0 0
0 7 7 0 5 0 0 0 2 0 0 2 0 0 6
7 7 0 0 0 7 7 0 0 8 0 2 2 0 6
0 7 7 0 7 7 7 0 8 8 8 0 0 9 0
10 0 7 7 7 7 7 7 0 8 0 11 0 9 9
0 0 7 0 7 0 0 7 7 0 11 11 0 0 0
0 0 0 7 7 7 7 7 0 0 11 0 12 12 0
0 13 0 7 0 7 7 0 0 14 0 14 0 0 16
15 0 7 7 0 7 7 7 0 14 14 14 0 16 16
0 0 0 7 7 7 7 0 17 0 14 14 14 0 0
0 18 0 7 7 0 0 0 0 0 0 0 0 19 0
0 0 0 0 7 0 0 20 0 0 21 21 0 0 22)</pre>
=={{header|Nim}}==
{{trans|Go}}
<syntaxhighlight lang="nim">import random, sequtils, strformat
const
N = 15
T = 5
P = 0.5
NotClustered = 1
Cell2Char = " #abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"
NRange = [4, 64, 256, 1024, 4096]
type Grid = seq[seq[int]]
proc newGrid(n: Positive; p: float): Grid =
result = newSeqWith(n, newSeq[int](n))
for row in result.mitems:
for cell in row.mitems:
if rand(1.0) < p: cell = 1
func walkMaze(grid: var Grid; m, n, idx: int) =
grid[n][m] = idx
if n < grid.high and grid[n + 1][m] == NotClustered:
grid.walkMaze(m, n + 1, idx)
if m < grid[0].high and grid[n][m + 1] == NotClustered:
grid.walkMaze(m + 1, n, idx)
if m > 0 and grid[n][m - 1] == NotClustered:
grid.walkMaze(m - 1, n, idx)
if n > 0 and grid[n - 1][m] == NotClustered:
grid.walkMaze(m, n - 1, idx)
func clusterCount(grid: var Grid): int =
var walkIndex = 1
for n in 0..grid.high:
for m in 0..grid[0].high:
if grid[n][m] == NotClustered:
inc walkIndex
grid.walkMaze(m, n, walkIndex)
result = walkIndex - 1
proc clusterDensity(n: int; p: float): float =
var grid = newGrid(n, p)
result = grid.clusterCount() / (n * n)
proc print(grid: Grid) =
for n, row in grid:
stdout.write n mod 10, ") "
for cell in row:
stdout.write ' ', Cell2Char[cell]
stdout.write '\n'
when isMainModule:
randomize()
var grid = newGrid(N, 0.5)
echo &"Found {grid.clusterCount()} clusters in this {N} by {N} grid\n"
grid.print()
echo ""
for n in NRange:
var sum = 0.0
for _ in 1..T:
sum += clusterDensity(n, P)
let sim = sum / T
echo &"t = {T} p = {P:4.2f} n = {n:4} sim = {sim:7.5f}"</syntaxhighlight>
{{out}}
<pre>Found 14 clusters in this 15 by 15 grid
0) a a b c c c c d
1) e e c c c c c c c
2) f e c c c c c
3) f f c c c c c c c c c
4) f c c c c c c c c c c c
5) g c c c c c c c
6) h i c c c c c c c c
7) i i c c c
8) j j j j c c k
9) l j c k k
0) l l l l j j k
1) l l l l l l j j j j j
2) l l l m j j j j
3) l l l l j j
4) n l l l l j j j j j j j j
t = 5 p = 0.50 n = 4 sim = 0.11250
t = 5 p = 0.50 n = 64 sim = 0.07085
t = 5 p = 0.50 n = 256 sim = 0.06702
t = 5 p = 0.50 n = 1024 sim = 0.06619
t = 5 p = 0.50 n = 4096 sim = 0.06587</pre>
=={{header|Perl}}==
{{trans|Raku}}
<syntaxhighlight lang="perl">$fill = 'x';
$D{$_} = $i++ for qw<DeadEnd Up Right Down Left>;
sub deq { defined $_[0] && $_[0] eq $_[1] }
sub perctest {
my($grid) = @_;
generate($grid);
my $block = 1;
for my $y (0..$grid-1) {
for my $x (0..$grid-1) {
fill($x, $y, $block++) if $perc[$y][$x] eq $fill
}
}
($block - 1) / $grid**2;
}
sub generate {
my($grid) = @_;
for my $y (0..$grid-1) {
for my $x (0..$grid-1) {
$perc[$y][$x] = rand() < .5 ? '.' : $fill;
}
}
}
sub fill {
my($x, $y, $block) = @_;
$perc[$y][$x] = $block;
my @stack;
while (1) {
if (my $dir = direction( $x, $y )) {
push @stack, [$x, $y];
($x,$y) = move($dir, $x, $y, $block)
} else {
return unless @stack;
($x,$y) = @{pop @stack};
}
}
}
sub direction {
my($x, $y) = @_;
return $D{Down} if deq($perc[$y+1][$x ], $fill);
return $D{Left} if deq($perc[$y ][$x-1], $fill);
return $D{Right} if deq($perc[$y ][$x+1], $fill);
return $D{Up} if deq($perc[$y-1][$x ], $fill);
return $D{DeadEnd};
}
sub move {
my($dir,$x,$y,$block) = @_;
$perc[--$y][ $x] = $block if $dir == $D{Up};
$perc[++$y][ $x] = $block if $dir == $D{Down};
$perc[ $y][ --$x] = $block if $dir == $D{Left};
$perc[ $y][ ++$x] = $block if $dir == $D{Right};
($x, $y)
}
my $K = perctest(15);
for my $row (@perc) {
printf "%3s", $_ for @$row;
print "\n";
}
printf "𝘱 = 0.5, 𝘕 = 15, 𝘒 = %.4f\n\n", $K;
$trials = 5;
for $N (10, 30, 100, 300, 1000) {
my $total = 0;
$total += perctest($N) for 1..$trials;
printf "𝘱 = 0.5, trials = $trials, 𝘕 = %4d, 𝘒 = %.4f\n", $N, $total / $trials;
}</syntaxhighlight>
{{out}}
<pre> 1 1 1 . . . . 2 2 2 . . . . .
. 1 . 1 1 1 . 2 2 . 2 2 2 . 3
. 1 . . 1 . 2 2 2 2 2 2 . . 3
1 1 1 . 1 . 2 2 . . . . 4 4 .
1 1 1 . 1 . . 2 . . . . . . 1
1 1 1 1 1 . . 2 . . 5 . 6 . .
1 1 . . 1 1 . 2 . 7 . . . 1 1
1 . . . 1 1 . 2 2 . . 8 8 . 1
. 9 9 9 . 1 . . 2 2 . . . 1 1
. . 9 9 . . 10 . . . 11 . 12 . .
. 9 9 . 13 13 . 13 . 14 . . 12 . .
15 . . 13 13 13 13 13 . . . 16 . 17 .
15 . . 13 . 13 . 13 13 . . 16 16 . .
. 18 . . 13 13 13 13 . . . . . 19 19
1 . 1 . . 13 . . . . 20 . 19 19 .
𝘱 = 0.5, 𝘕 = 15, 𝘒 = 0.0889
𝘱 = 0.5, trials = 5, 𝘕 = 10, 𝘒 = 0.0980
𝘱 = 0.5, trials = 5, 𝘕 = 30, 𝘒 = 0.0738
𝘱 = 0.5, trials = 5, 𝘕 = 100, 𝘒 = 0.0670
𝘱 = 0.5, trials = 5, 𝘕 = 300, 𝘒 = 0.0660
𝘱 = 0.5, trials = 5, 𝘕 = 1000, 𝘒 = 0.0661</pre>
=={{header|Phix}}==
{{trans|C}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">grid</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">w</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ww</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">make_grid</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">ww</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">w</span><span style="color: #0000FF;">*</span><span style="color: #000000;">w</span>
<span style="color: #000000;">grid</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ww</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">ww</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">grid</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-(</span><span style="color: #7060A8;">rnd</span><span style="color: #0000FF;">()<</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">alpha</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"+.ABCDEFGHIJKLMNOPQRSTUVWXYZ"</span><span style="color: #0000FF;">&</span>
<span style="color: #008000;">"abcdefghijklmnopqrstuvwxyz"</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">show_cluster</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">ww</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">gi</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">grid</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]+</span><span style="color: #000000;">2</span>
<span style="color: #000000;">grid</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gi</span><span style="color: #0000FF;"><=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">alpha</span><span style="color: #0000FF;">)?</span><span style="color: #000000;">alpha</span><span style="color: #0000FF;">[</span><span style="color: #000000;">gi</span><span style="color: #0000FF;">]:</span><span style="color: #008000;">'?'</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">join_by</span><span style="color: #0000FF;">(</span><span style="color: #000000;">grid</span><span style="color: #0000FF;">,</span><span style="color: #000000;">w</span><span style="color: #0000FF;">,</span><span style="color: #000000;">w</span><span style="color: #0000FF;">,</span><span style="color: #008000;">""</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">recur</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">1</span> <span style="color: #008080;">and</span> <span style="color: #000000;">x</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">ww</span> <span style="color: #008080;">and</span> <span style="color: #000000;">grid</span><span style="color: #0000FF;">[</span><span style="color: #000000;">x</span><span style="color: #0000FF;">]==-</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">grid</span><span style="color: #0000FF;">[</span><span style="color: #000000;">x</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span>
<span style="color: #000000;">recur</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">w</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">recur</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">recur</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">recur</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">w</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">count_clusters</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">cls</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">ww</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">grid</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]=-</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">cls</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #000000;">recur</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">cls</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">cls</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">tests</span><span style="color: #0000FF;">(</span><span style="color: #004080;">int</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">k</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">make_grid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">k</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">count_clusters</span><span style="color: #0000FF;">()/</span><span style="color: #000000;">ww</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">k</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">n</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<span style="color: #000000;">w</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">15</span>
<span style="color: #000000;">make_grid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"width=15, p=0.5, %d clusters:\n"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">count_clusters</span><span style="color: #0000FF;">())</span>
<span style="color: #000000;">show_cluster</span><span style="color: #0000FF;">()</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\np=0.5, iter=5:\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">w</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">4</span>
<span style="color: #008080;">while</span> <span style="color: #000000;">w</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">4096</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%5d %9.6f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">w</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">tests</span><span style="color: #0000FF;">(</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">)})</span>
<span style="color: #000000;">w</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">4</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
width=15, p=0.5, 18 clusters:
..EE.FF...OO..P
.......KK.OOO.P
A..B..I....OO.P
.BBB.G.LL....PP
BBB...J....J.P.
..BB..JJJJJJ.P.
.BBBB.JJJ.J.J..
BBB....JJJ.JJ.Q
BB.D.D.J.JJJ.QQ
..DDDD.JJ.JJJ.Q
.DDDD.JJ..JJ...
C.DDD.J.N.JJ...
C.D...J..JJ....
.....H...J.....
.......M..O...R
p=0.5, iter=5:
4 0.137500
16 0.080469
64 0.068164
256 0.066809
1024 0.066018
4096 0.065777
</pre>
=={{header|Python}}==
<
from random import random
import string
Line 324 ⟶ 1,712:
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))</
{{out}}
Line 353 ⟶ 1,741:
=={{header|Racket}}==
<
(require srfi/14) ; character sets
Line 449 ⟶ 1,837:
(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))</
{{out}}
Line 482 ⟶ 1,870:
| . .. .....| | F HH AAAAA|
8 clusters</pre>
=={{header|Raku}}==
(formerly Perl 6)
{{works with|Rakudo|2017.02}}
<syntaxhighlight lang="raku" line>my @perc;
my $fill = 'x';
enum Direction <DeadEnd Up Right Down Left>;
my $𝘒 = perctest(15);
.fmt("%-2s").say for @perc;
say "𝘱 = 0.5, 𝘕 = 15, 𝘒 = $𝘒\n";
my $trials = 5;
for 10, 30, 100, 300, 1000 -> $𝘕 {
my $𝘒 = ( [+] perctest($𝘕) xx $trials ) / $trials;
say "𝘱 = 0.5, trials = $trials, 𝘕 = $𝘕, 𝘒 = $𝘒";
}
sub infix:<deq> ( $a, $b ) { $a.defined && ($a eq $b) }
sub perctest ( $grid ) {
generate $grid;
my $block = 1;
for ^$grid X ^$grid -> ($y, $x) {
fill( [$x, $y], $block++ ) if @perc[$y; $x] eq $fill
}
($block - 1) / $grid²;
}
sub generate ( $grid ) {
@perc = ();
@perc.push: [ ( rand < .5 ?? '.' !! $fill ) xx $grid ] for ^$grid;
}
sub fill ( @cur, $block ) {
@perc[@cur[1]; @cur[0]] = $block;
my @stack;
my $current = @cur;
loop {
if my $dir = direction( $current ) {
@stack.push: $current;
$current = move $dir, $current, $block
}
else {
return unless @stack;
$current = @stack.pop
}
}
sub direction( [$x, $y] ) {
( Down if @perc[$y + 1][$x] deq $fill ) ||
( Left if @perc[$y][$x - 1] deq $fill ) ||
( Right if @perc[$y][$x + 1] deq $fill ) ||
( Up if @perc[$y - 1][$x] deq $fill ) ||
DeadEnd
}
sub move ( $dir, @cur, $block ) {
my ( $x, $y ) = @cur;
given $dir {
when Up { @perc[--$y; $x] = $block }
when Down { @perc[++$y; $x] = $block }
when Left { @perc[$y; --$x] = $block }
when Right { @perc[$y; ++$x] = $block }
}
[$x, $y]
}
}
</syntaxhighlight>
{{out}}
<pre>. . 1 . 2 . . 3 . . . 4 . . .
2 2 . . 2 2 2 . 5 5 . 4 4 4 4
2 2 2 2 2 . 2 . 5 . . 4 . . 4
2 . 2 . 2 2 . . . . 4 4 4 . 4
. . . . . 2 . . . . 4 4 4 . .
6 6 6 6 . . 7 7 . . 4 . 4 4 .
6 . 6 . . . . . 4 4 4 4 4 . .
6 6 6 . . . 8 8 . 4 4 4 . . 4
6 . . . 9 . . . . . . 4 4 . 4
. 10 . 11 . . 12 12 . 4 . . 4 4 4
11 . 11 11 11 11 . 12 . 4 . 4 4 4 .
11 11 11 11 11 11 11 . . 4 4 . . 4 .
11 11 11 11 . 11 . . . 4 4 4 4 4 .
11 11 11 . 11 11 11 . . 4 4 4 4 . 13
. 11 11 . 11 11 . . . . 4 4 . 14 .
𝘱 = 0.5, 𝘕 = 15, 𝘒 = 0.062222
𝘱 = 0.5, trials = 5, 𝘕 = 10, 𝘒 = 0.114
𝘱 = 0.5, trials = 5, 𝘕 = 30, 𝘒 = 0.082444
𝘱 = 0.5, trials = 5, 𝘕 = 100, 𝘒 = 0.06862
𝘱 = 0.5, trials = 5, 𝘕 = 300, 𝘒 = 0.066889
𝘱 = 0.5, trials = 5, 𝘕 = 1000, 𝘒 = 0.0659358
</pre>
=={{header|Tcl}}==
Note that the queue (variables <code>q</code> and <code>k</code>) used to remember where to find cells when flood-filling the cluster is maintained as a list ''segment''; the front of the list is not trimmed for performance reasons. (This would matter with very long queues, in which case the queue could be shortened occasionally; ''frequent'' trimming is still slower though, because Tcl backs its “list” datatype with arrays and not linked lists.)
{{works with|Tcl|8.6}}
<syntaxhighlight lang="tcl">package require Tcl 8.6
proc determineClusters {w h p} {
# Construct the grid
set grid [lrepeat $h [lrepeat $w 0]]
for {set i 0} {$i < $h} {incr i} {
for {set j 0} {$j < $w} {incr j} {
lset grid $i $j [expr {rand() < $p ? -1 : 0}]
}
}
# Find (and count) the clusters
set cl 0
for {set i 0} {$i < $h} {incr i} {
for {set j 0} {$j < $w} {incr j} {
if {[lindex $grid $i $j] == -1} {
incr cl
for {set q [list $i $j];set k 0} {$k<[llength $q]} {incr k} {
set y [lindex $q $k]
set x [lindex $q [incr k]]
if {[lindex $grid $y $x] != -1} continue
lset grid $y $x $cl
foreach dx {1 0 -1 0} dy {0 1 0 -1} {
set nx [expr {$x+$dx}]
set ny [expr {$y+$dy}]
if {
$nx >= 0 && $ny >= 0 && $nx < $w && $ny < $h &&
[lindex $grid $ny $nx] == -1
} then {
lappend q $ny $nx
}
}
}
}
}
}
return [list $cl $grid]
}
# Print a sample 15x15 grid
lassign [determineClusters 15 15 0.5] n g
puts "15x15 grid, p=0.5, with $n clusters"
puts "+[string repeat - 15]+"
foreach r $g {puts |[join [lmap x $r {format %c [expr {$x==0?32:64+$x}]}] ""]|}
puts "+[string repeat - 15]+"
# Determine the densities as the grid size increases
puts "p=0.5, iter=5"
foreach n {5 30 180 1080 6480} {
set tot 0
for {set i 0} {$i < 5} {incr i} {
lassign [determineClusters $n $n 0.5] nC
incr tot $nC
}
puts "n=$n, K(p)=[expr {$tot/5.0/$n**2}]"
}</syntaxhighlight>
{{out}}
<pre>
15x15 grid, p=0.5, with 21 clusters
+---------------+
| A B CCCCC|
| D A BBB C |
|E B F CCCC|
| B B F CC C|
|BBB B BB CCC|
|B BBBBBB CCCCC|
| B B G C C|
|H II G G J |
|HH II G GG K|
|HH II GGG GG K|
| I G GGGG |
|LL GGG GG M N|
| L G G O P |
|LLLL Q R |
|L L S T UUU|
+---------------+
p=0.5, iter=5
n=5, K(p)=0.184
n=30, K(p)=0.07155555555555557
n=180, K(p)=0.06880246913580246
n=1080, K(p)=0.0661267146776406
n=6480, K(p)=0.06582889898643499
</pre>
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "random" for Random
import "./fmt" for Fmt
var rand = Random.new()
var RAND_MAX = 32767
var list = []
var w = 0
var ww = 0
var ALPHA = "+.ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"
var ALEN = ALPHA.count - 3
var makeList = Fn.new { |p|
var thresh = (p * RAND_MAX).truncate
ww = w * w
var i = ww
list = List.filled(i, 0)
while (i != 0) {
i = i - 1
var r = rand.int(RAND_MAX+1)
if (r < thresh) list[i] = -1
}
}
var showCluster = Fn.new {
var k = 0
for (i in 0...w) {
for (j in 0...w) {
var s = list[k]
k = k + 1
var c = (s < ALEN) ? ALPHA[1 + s] : "?"
System.write(" %(c)")
}
System.print()
}
}
var recur // recursive
recur = Fn.new { |x, v|
if (x >= 0 && x < ww && list[x] == -1) {
list[x] = v
recur.call(x - w, v)
recur.call(x - 1, v)
recur.call(x + 1, v)
recur.call(x + w, v)
}
}
var countClusters = Fn.new {
var cls = 0
for (i in 0...ww) {
if (list[i] == -1) {
cls = cls + 1
recur.call(i, cls)
}
}
return cls
}
var tests = Fn.new { |n, p|
var k = 0
for (i in 0...n) {
makeList.call(p)
k = k + countClusters.call() / ww
}
return k / n
}
w = 15
makeList.call(0.5)
var cls = countClusters.call()
System.print("width = 15, p = 0.5, %(cls) clusters:")
showCluster.call()
System.print("\np = 0.5, iter = 5:")
w = 1 << 2
while (w <= (1 << 13)) {
var t = tests.call(5, 0.5)
Fmt.print("$5d $9.6f", w, t)
w = w << 1
}</syntaxhighlight>
{{out}}
Sample run:
<pre>
width = 15, p = 0.5, 28 clusters:
. A A . . B . C . D . . . . A
A A A . E . F . G . . A A A A
. A A A . F F . . . A A . . .
. A . . . . F . . H . A . A .
. A A . . . . . H H H . A A A
A A A . . I I . . H . J . . A
A . . K . . . . H H . . . A A
A . L . M M . N . . O . . . .
A A . P . . . N . Q . R . S .
. A A . T T . . U . R R R . .
A A . . T T T . . . . R . . .
. A . T T T . . . . . R R R R
. A . . . T . . V . W . R R .
A A . X . T T . . . W . . . .
. . Y . T T . Z . a . . . b b
p = 0.5, iter = 5:
4 0.125000
8 0.081250
16 0.094531
32 0.073242
64 0.067920
128 0.068567
256 0.067404
512 0.066401
1024 0.065870
2048 0.065885
4096 0.065812
8192 0.065795
</pre>
=={{header|zkl}}==
{{trans|C}}
<syntaxhighlight lang="zkl">const X=-1; // the sentinal that marks an untouched cell
var C,N,NN,P;
fcn createC(n,p){
N,P=n,p; NN=N*N;
C=NN.pump(List.createLong(NN),0); // vector of ints
foreach n in (NN){ C[n]=X*(Float.random(1)<=P) } // X is the sentinal
}
fcn showCluster{
alpha:="-ABCDEFGHIJKLMNOPQRSTUVWXYZ" "abcdefghijklmnopqrstuvwxyz";
foreach n in ([0..NN,N]){ C[n,N].pump(String,alpha.get).println() }
}
fcn countClusters{
clusters:=0;
foreach n in (NN){
if(X!=C[n]) continue;
fcn(n,v){
if((0<=n<NN) and C[n]==X){
C[n]=v;
self.fcn(n-N,v); self.fcn(n-1,v); self.fcn(n+1,v); self.fcn(n+N,v);
}
}(n,clusters+=1);
}
clusters
}
fcn tests(N,n,p){
k:=0.0;
foreach z in (n){ createC(N,p); k+=countClusters().toFloat()/NN; }
k/n
}</syntaxhighlight>
<syntaxhighlight lang="zkl">createC(15,0.5);
println("width=%d, p=%.1f, %d clusters:".fmt(N,P,countClusters()));
showCluster();
println("p=0.5, 5 iterations:");
w:=4; do(6){ println("%5d %9.6f".fmt(w,tests(w, 5, 0.5))); w*=4; }</syntaxhighlight>
{{out}}
<pre>
width=15, p=0.5, 16 clusters:
-AAA-BB-BBB---C
------BBBB--D--
E---F---BB--DD-
EE----G-BB---DD
--H-I--J--J--DD
-K--I--JJ-J--D-
-K--I--JJJJ-L--
KK-III-------MM
-K-I--I--NN-I--
I-IIIII-NNN-III
I-II--I-N-N-II-
III-III--NNN-II
I-II-II-O---I--
I-I-IIII-PP-III
I-II--I---P--II
p=0.5, 5 iterations:
4 0.062500
16 0.070312
64 0.067627
256 0.067078
1024 0.065834
4096 0.065771
</pre>
| |||