Percolation/Mean cluster density: Difference between revisions

m
m (Fix Perl 6 -> Raku links)
 
(23 intermediate revisions by 9 users not shown)
Line 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}}==
<langsyntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
 
Line 99 ⟶ 187:
free(map);
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>
Line 127 ⟶ 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}}
<langsyntaxhighlight lang="d">import std.stdio, std.algorithm, std.random, std.math, std.array,
std.range, std.ascii;
 
Line 220 ⟶ 443:
nIters, prob, side, density);
}
}</langsyntaxhighlight>
{{out}}
<pre>Found 26 clusters in this 15 by 15 grid:
Line 252 ⟶ 475:
=={{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.
<langsyntaxhighlight lang="scheme">
(define-constant BLACK (rgb 0 0 0.6))
(define-constant WHITE -1)
Line 290 ⟶ 513:
(writeln 'n n 'Cn Cn 'density (// Cn (* n n) 5) )
(vector->pixels C)) ;; to screen
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 301 ⟶ 524:
 
=={{header|Factor}}==
<langsyntaxhighlight lang="factor">USING: combinators formatting generalizations kernel math
math.matrices random sequences ;
IN: rosetta-code.mean-cluster-density
Line 342 ⟶ 565:
] each ;
 
MAIN: main</langsyntaxhighlight>
{{out}}
<pre>
Line 354 ⟶ 577:
=={{header|Go}}==
{{trans|Python}}
<langsyntaxhighlight lang="go">package main
 
import (
Line 451 ⟶ 674:
fmt.Printf("t=%3d p=%4.2f n=%5d sim=%7.5f\n", t, p, n, sim)
}
}</langsyntaxhighlight>
{{out}}
<pre>
Line 478 ⟶ 701:
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}}==
Line 485 ⟶ 801:
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.
 
<langsyntaxhighlight Jlang="j">congeal=: |.@|:@((*@[*>.)/\.)^:4^:_</langsyntaxhighlight>
 
Example:
 
<langsyntaxhighlight Jlang="j"> M=:0.4>?6 6$0
M
1 0 0 0 0 0
Line 510 ⟶ 826:
71 71 0 0 0 0
0 0 0 149 0 113
131 131 0 149 149 0</langsyntaxhighlight>
 
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.
 
<langsyntaxhighlight Jlang="j">idclust=: $ $ [: (~. i.])&.(0&,)@,@congeal ] * 1 + i.@$</langsyntaxhighlight>
 
Example use:
 
<langsyntaxhighlight Jlang="j"> idclust M
1 0 0 0 0 0
0 0 0 2 0 0
Line 531 ⟶ 847:
CC....
...D.E
FF.DD.</langsyntaxhighlight>
 
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:
 
<langsyntaxhighlight Jlang="j">K=: (%&#~ }.@~.)&,</langsyntaxhighlight>
 
Example use:
 
<langsyntaxhighlight Jlang="j"> K idclust M
0.1666667</langsyntaxhighlight>
 
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.
 
<langsyntaxhighlight Jlang="j">experiment=: K@ idclust@: > 0 ?@$~ ,~</langsyntaxhighlight>
 
Example use:
 
<langsyntaxhighlight Jlang="j"> 0.4 experiment 6
0.1666667
0.4 experiment 6
0.1944444</langsyntaxhighlight>
 
The task wants us to perform at least five trials for sizes up to 1000 by 1000 with probability of 1 being 0.5:
 
<langsyntaxhighlight Jlang="j">trials=: 0.5&experiment"0@#</langsyntaxhighlight>
 
Example use:
 
<langsyntaxhighlight Jlang="j"> 6 trials 3
0.1111111 0.1111111 0.2222222 0.1111111 0.1111111 0.3333333
6 trials 10
Line 570 ⟶ 886:
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</langsyntaxhighlight>
 
Now for averages (these are different trials from the above):
 
<langsyntaxhighlight Jlang="j">mean=: +/%#
mean 8 trials 3
0.1805556
Line 586 ⟶ 902:
0.06749861
mean 8 trials 1000
0.06616738</langsyntaxhighlight>
 
Finally, for the extra credit (thru taken from the [[Loops/Downward_for#J|Loops/Downward for]] task):
 
<langsyntaxhighlight Jlang="j">thru=: <./ + i.@(+*)@-~</langsyntaxhighlight>
 
<langsyntaxhighlight Jlang="j"> (idclust 0.5 > 0 ?@$~ 15 15) {'.', 'A' thru&.(a.&i.) 'Z'
A.......B..C...
AAAA...D..E.F..
Line 607 ⟶ 923:
..AA..A.A...AAA
.M.A.AA.AA..AA.
.MM..A.N..O..A.</langsyntaxhighlight>
 
'''Collected definitions'''
 
<langsyntaxhighlight Jlang="j">congeal=: |.@|:@((*@[*>.)/\.)^:4^:_
idclust=: $ $ [: (~. i.])&.(0&,)@,@congeal ] * 1 + i.@$
 
Line 621 ⟶ 937:
mean=:+/ % #
 
thru=: <./ + i.@(+*)@-~</langsyntaxhighlight>
 
'''Extra Credit'''
 
<langsyntaxhighlight Jlang="j"> M=: (* 1+i.@$)?15 15$2
M
0 2 3 4 0 6 0 8 0 10 11 12 0 0 15
Line 673 ⟶ 989:
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</langsyntaxhighlight>
 
=={{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}}
 
<langsyntaxhighlight lang="julia">using Printf, Distributions
 
newgrid(p::Float64, r::Int, c::Int=r) = rand(Bernoulli(p), r, c)
Line 728 ⟶ 1,177:
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</langsyntaxhighlight>
 
{{out}}
Line 757 ⟶ 1,206:
=={{header|Kotlin}}==
{{trans|C}}
<langsyntaxhighlight lang="scala">// version 1.2.10
 
import java.util.Random
Line 836 ⟶ 1,285:
w = w shl 1
}
}</langsyntaxhighlight>
 
Sample output:
Line 871 ⟶ 1,320:
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}}
<langsyntaxhighlight lang="perl">$fill = 'x';
$D{$_} = $i++ for qw<DeadEnd Up Right Down Left>;
 
Line 945 ⟶ 1,525:
$total += perctest($N) for 1..$trials;
printf "𝘱 = 0.5, trials = $trials, 𝘕 = %4d, 𝘒 = %.4f\n", $N, $total / $trials;
}</langsyntaxhighlight>
{{out}}
<pre> 1 1 1 . . . . 2 2 2 . . . . .
Line 972 ⟶ 1,552:
=={{header|Phix}}==
{{trans|C}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<lang Phix>sequence grid
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
integer w, ww
<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>
procedure make_grid(atom p)
ww = w*w
<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>
grid = repeat(0,ww)
<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>
for i=1 to ww do
<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>
grid[i] = -(rnd()<p)
<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>
end for
<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>
end procedure
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
constant alpha = "+.ABCDEFGHIJKLMNOPQRSTUVWXYZ"&
"abcdefghijklmnopqrstuvwxyz"
<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>
procedure show_cluster()
for i=1 to ww do
<span style="color: #008080;">procedure</span> <span style="color: #000000;">show_cluster</span><span style="color: #0000FF;">()</span>
integer gi = grid[i]+2
<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>
grid[i] = iff(gi<=length(alpha)?alpha[gi]:'?')
<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>
end for
<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>
puts(1,join_by(grid,w,w,""))
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end procedure
<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>
procedure recur(integer x, v)
if x>=1 and x<=ww and grid[x]==-1 then
<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>
grid[x] = v
<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>
recur(x-w, v)
<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>
recur(x-1, v)
<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>
recur(x+1, v)
<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>
recur(x+w, v)
<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>
end if
<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>
end procedure
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
function count_clusters()
integer cls = 0
<span style="color: #008080;">function</span> <span style="color: #000000;">count_clusters</span><span style="color: #0000FF;">()</span>
for i=1 to ww do
<span style="color: #004080;">integer</span> <span style="color: #000000;">cls</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
if grid[i]=-1 then
<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>
cls += 1
<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>
recur(i, cls)
<span style="color: #000000;">cls</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
end if
<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>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
return cls
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<span style="color: #008080;">return</span> <span style="color: #000000;">cls</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
function tests(int n, atom p)
atom k = 0
<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>
for i=1 to n do
<span style="color: #004080;">atom</span> <span style="color: #000000;">k</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
make_grid(p)
<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>
k += count_clusters()/ww
<span style="color: #000000;">make_grid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
end for
<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>
return k / n
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end function
<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>
procedure main()
w = 15
<span style="color: #008080;">procedure</span> <span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
make_grid(0.5)
<span style="color: #000000;">w</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">15</span>
printf(1,"width=15, p=0.5, %d clusters:\n", count_clusters())
<span style="color: #000000;">make_grid</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">)</span>
show_cluster()
<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>
printf(1,"\np=0.5, iter=5:\n")
w = 4
<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>
while w<=4096 do
<span style="color: #000000;">w</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">4</span>
printf(1,"%5d %9.6f\n", {w, tests(5,0.5)})
<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>
w *= 4
<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>
end while
<span style="color: #000000;">w</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">4</span>
end procedure
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
main()</lang>
<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>
Line 1,067 ⟶ 1,650:
 
=={{header|Python}}==
<langsyntaxhighlight lang="python">from __future__ import division
from random import random
import string
Line 1,129 ⟶ 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))</langsyntaxhighlight>
 
{{out}}
Line 1,158 ⟶ 1,741:
 
=={{header|Racket}}==
<langsyntaxhighlight lang="racket">#lang racket
(require srfi/14) ; character sets
 
Line 1,254 ⟶ 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))</langsyntaxhighlight>
 
{{out}}
Line 1,292 ⟶ 1,875:
{{works with|Rakudo|2017.02}}
 
<syntaxhighlight lang="raku" perl6line>my @perc;
my $fill = 'x';
 
Line 1,358 ⟶ 1,941:
}
}
</syntaxhighlight>
</lang>
{{out}}
<pre>. . 1 . 2 . . 3 . . . 4 . . .
Line 1,387 ⟶ 1,970:
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}}
<langsyntaxhighlight lang="tcl">package require Tcl 8.6
 
proc determineClusters {w h p} {
Line 1,441 ⟶ 2,024:
}
puts "n=$n, K(p)=[expr {$tot/5.0/$n**2}]"
}</langsyntaxhighlight>
{{out}}
<pre>
Line 1,468 ⟶ 2,051:
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}}
<langsyntaxhighlight lang="zkl">const X=-1; // the sentinal that marks an untouched cell
var C,N,NN,P;
fcn createC(n,p){
Line 1,500 ⟶ 2,204:
foreach z in (n){ createC(N,p); k+=countClusters().toFloat()/NN; }
k/n
}</langsyntaxhighlight>
<langsyntaxhighlight 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; }</langsyntaxhighlight>
{{out}}
<pre>
1,480

edits