# Mimic a procedural, recursive clustering algorithm for site percolation using functional programming

Sorry in advance for my logorrhea: I just want to make sure all of the information is here.

Context and Question

I am investigating site percolation on a square lattice. I have a working, depth-first, recursive clustering algorithm. The algorithm as coded is recursive but procedural, so I want to find a way to mimic the recursive nature of the algorithm using Mathematica's functional-programming style. In addition, the Kernel dies if the number $l^2$ of grid points is too large and if the filling fraction is large enough that the percolating cluster is very large; I presume this is due to the large recursion depth of the search. (If coded in C, I can do much larger lattice sizes and filling fractions.) Therefore, I would like a way of getting around this problem of the Kernel crashing. In summary, I would like to:

• Mimic the recursive algorithm with functional-programming.

• Get around the problem of the Kernel crashing due to the large recursion depth of the algorithm (presumably this is solved by fulfilling the previous bullet point).

I have tried to come up with a way of implementing this using functions like Nest or Fold, but I can't see how to mimic the tree-nature of the recursive search, since they seem not to be able to deal with trees.

One last point: this algorithm is fast (EDIT: apparently, no, it's not: see Update below). It may be that implementing a functional programming version in Mathematica will be slower. In that case, I am still interested in the answer for the purpose of honing my functional programming skills.

Algorithm Description

The input is an $l\times l$ array configuration consisting of $-1$'s at filled sites and $0$'s elsewhere. The output consists of a list of {clusterNumber,clusterSize}s and an $l\times l$ array with the $-1$'s replaced by the corresponding cluster number and the $0$'s left alone. Sites belong to the same cluster if they are nearest neighbors (diagonals not allowed). As an example in a $4\times4$ lattice,

inputConfiguration={{-1,-1,-1,0},{-1,0,0,-1},{-1,-1,0,-1},{0,-1,0,0}};
output={{1,7} , {2,2} , {{1,1,1,0},{1,0,0,2},{1,1,0,2},{0,1,0,0}}};


Below, you will see that the configurations are Flattened. This is merely a programming convenience. To go from one version to another, I either Flatten[configuration] or Partition[configuration,l].

Before talking about the algorithm specifically, note that

neighbors[l_]= DeleteCases[{If[-l + # >= 0, -l + #, -1],
If[Quotient[-1 + #, l] == Quotient[#, l], -1 + #, -1],
If[Quotient[1 + #, l] == Quotient[#, l], 1 + #, -1],
If[l + # <= l^2 - 1, l + #, -1]} + 1, 0] & /@ Range[0, l^2 - 1]


is a function that returns a list of the indices of the neighboring sites of each site, where the indices correspond to the flattened configuration. For instance,

neighbors[] = {1,3,6}


because site number 2 is bordered by site 1 on the left, 3 on the right, and 6 below. (This can certainly be implemented in a nicer way, so feel free to fix this, but this is not the focus of this question.)

The heart of the algorithm is the recursive function cluster that accepts the current site index i and cluster number clusterNumber as inputs, sets the current grid value to clusterNumber, then checks the neighbors of site i. If a neighbor site j is a filled site but has not yet been assigned a cluster number, cluster is called with inputs j and clusterNumber. A function call to cluster will terminate once the input site j has no neighbors not yet part of a cluster, and once all recursive calls to cluster have terminated, all sites in the current cluster have been labeled with clusterNumber.

At this point, I Sow both the current cluster number and the current cluster size and then move on. At the end, I Sow the clustered configuration.

Code

In practice, I run the code in a Block in order to set the $RecursionLimit larger, which is necessary for larger lattice sizes and larger filling fractions. For instance, for the$4\times4$configuration above, we would call the function as Block[{$RecursionLimit=10000}, Reap[clustering[inputConfiguration, neighbors,4]]]


Here is the code:

clustering[configuration_, boundaries_, l_] := Module[{
grid = Flatten[configuration]
,cluster
,clusterNumber = 1
,clusterSize
}
,cluster[i_, num_] := Module[{}
,grid[[i]] = num
;clusterSize = clusterSize + 1
;If[grid[[#]] < 0, cluster[#, num]] & /@ boundaries[[i]]
]
;Do[
If[grid[[m]] < 0
,clusterSize = 0
;cluster[m, clusterNumber]
;Sow[{clusterNumber++, clusterSize}]
]
, {m, 1, l^2}]
;Sow[Partition[grid,l]]
]
]


Further Goals

It is also useful to return a list of {clusterNumber,clusterIndices}s, where clusterIndices is a list of the site indices of the cluster labeled clusterNumber. Presumably, given code that mimics the code above with functional programming, I would be able to modify it in order to return the clusterIndices rather than the clusterSize, so this is not a top priority, but perhaps it's important to keep in mind.

Update: Benchmarking the cluster algorithms

I have stripped out the clusterSize and clusterIndex finding from the original code and Virgil's two versions (the image-processing version and the recursive version) in order to benchmark just the clustering algorithms.

Specifically, I have a random configuration function that randomly chooses site indices to fill and returns the grid in the correct format:

randomConfiguration[l_, nf] := Block[
{grid = ConstantArray[0, l^2], config = RandomSample[Range[l^2], nf]}
, grid[[config]] = -1; Partition[grid,l]]


I have then defined the three clustering codes. First, the original, essentially:

march[config_] := Module[{grid = Flatten[config], l = Length@config, clusterSearch, num = 1, boundaries}
, boundaries = DeleteCases[{If[-l + # >= 0, -l + #, -1],
If[Quotient[-1 + #, l] == Quotient[#, l], -1 + #, -1],
If[Quotient[1 + #, l] == Quotient[#, l], 1 + #, -1],
If[l + # <= l^2 - 1, l + #, -1]} + 1, 0] & /@ Range[0, l^2 - 1]
; clusterSearch[i_, num_] := Module[{}
, grid[[i]] = num
; If[grid[[#]] < 0, clusterSearch[#, num]] & /@ boundaries[[i]]
]
; Do[If[grid[[m]] < 0 , clusterSearch[m, num] ; num++ ] , {m, 1, l^2}]
; Partition[grid, l]
]


Then, Virgil's recursive code:

virgilRec[config_] := Module[{output = config, cnum = 0, length = Length@config, csearch}
, csearch[{i_, j_}, cnum_] := If[output[[i, j]] == -1
, output[[i, j]] = cnum
; csearch[#, cnum] & /@ Select[{i, j} + # & /@ {{0, 1}, {0, -1}, {1, 0}, {-1, 0}}, (1 <= First@# <= length && 1 <= Last@# <= length) &]
]
; Scan[
If[output[[Sequence @@ #]] == -1
, cnum++
; csearch[#, cnum]
] &
, SparseArray[config]["NonzeroPositions"]
];
output];


Finally, Virgil's image-processing code:

virgilImage[config_] := Module[{output, csizes, cindices}
, output = MorphologicalComponents[Image@Abs@config, CornerNeighbors -> False]
; output];


For each, I chose 50 random configurations at a filling fraction of $\approx0.593$ for various choices of $l$ and averaged the times taken to complete the clustering process, using

benchmark[f_, l, frac_] := With[{config = randomConfiguration[l,frac]}
, {l, Mean@Table[First@AbsoluteTiming@f@config, {50}]}
]


The results are interesting! In both recursive codes, the clustering time grows as a power law with the grid dimension from the beginning, whereas in the image processing version, the time is flat for a good while. Eventually, the clustering time in virgilImage grows as a power law as well. The slopes of the linear regions of the plots are all basically 2, which indicates that the clustering time grows linearly with the number $l^2$ of grid points for large enough $l$. virgilImage is a clear winner from the get-go.

To do

• Per a comment, see whether virgilImage is faster if we don't convert the grid to an image first.
• See if we can implement SparseArray as a convenient data structure for the grid.
• Take advantage of the fact that in choosing the random configuration, we have the indices of the filled sites from the beginning, which means we can combine the grid "data structure" with the list of filled site indices in order not to have to call the NonzeroPositions Method.
• Wait for more answers (hopefully)! I am still interested in seeing the experts' take on Mathematica-style recursive algorithms.
• "Sorry in advance for my logorrhea…" - in this case, nothing to be sorry for; most askers tend to give much less info that will help us help them. Jun 4 '15 at 20:43
• @J. M. All right then! No more apologies. :) Jun 4 '15 at 20:44

Built-in option

This sidesteps most of your code, so it might not be what you are looking for, but I believe your goal can be achieved with Mathematica's built-in image processing capability, specifically: MorphologicalComponents!

Define a new clustering function

clustering1[config_] := Module[{output, csizes, cindices},
output =  MorphologicalComponents[Image@Abs@config, CornerNeighbors -> False];
csizes = Rest@Sort@Tally@Flatten@output;
cindices = Module[
{sa = SparseArray[output], xx, yy, sa1, sa2},
sa1 = sa["NonzeroValues"];
xx = GatherBy[Range@Length@sa1, sa1[[#]] &];
sa2 = sa["NonzeroPositions"];
yy = sa2[[#]] & /@ xx;
Transpose[{sa1[[xx[[All, 1]]]], yy}]];
];
{csizes, cindices, output}
];


and apply:

inputConfig = {{-1, -1, -1, 0}, {-1, 0, 0, -1}, {-1, -1, 0, -1}, {0, -1, 0, 0}};
clustering1@inputConfig

{{{1, 7}, {2, 2}},
{{1, {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {3, 1}, {3, 2}, {4, 2}}}, {2, {{2, 4}, {3, 4}}}},
{{1, 1, 1, 0}, {1, 0, 0, 2}, {1, 1, 0, 2}, {0, 1, 0, 0}}}


The first item in the output is the list of {clusterNumber, clusterSize}, the second the list of {clusterNumber, clusterIndices}, and the third is the input array with cluster numbers replacing filled sites.

EDIT

It appears that finding the cluster indices with Position is extremely slow with large datasets, although the rest of the code is relatively fast. With thanks to ciao (see below), I've replaced it with a much faster construction.

Mimicking the recursive algorithm

Here is my take on what the OP was actually asking for: a Mathematic approach to a recursive algorithm. The actual algorithm is not much changed - the use of Sow and Reap in the OP's code is a good idea, and I agree that the recursion needed is not amenable to Nest or Fold - all I can do is clean it up a little and exchange the Do for a Scan over the populated sites:

clustering2[config_] := Module[
{output = config,
cnum = 0, length = Length@config,
csearch, clusters, c},

csearch[{i_, j_}, cnum_] := If[
output[[i, j]] == -1,
output[[i, j]] = cnum;
Sow[{i, j}];
csearch[#, cnum] & /@ Select[
{{i, j + 1}, {i, j - 1}, {i + 1, j}, {i - 1, j}},
(1 <= First@# <= length && 1 <= Last@# <= length) &]
];

clusters = Reap[
Scan[
If[output[[Sequence @@ #]] == -1,
cnum++;
c = Reap[csearch[#, cnum]][[2, 1]];
Sow[{cnum, Length@c, c}]] &,
SparseArray[config]["NonzeroPositions"]]
][[2, 1]];

{clusters, output}
];


Applying this to inputConfig gives

{{{1, 7, {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {3, 1}, {3, 2}, {4, 2}}},
{2, 2, {{2, 4}, {3, 4}}}},
{{1, 1, 1, 0}, {1, 0, 0, 2}, {1, 1, 0, 2}, {0, 1, 0, 0}}}


The first item in the output is the list of {clusterNumber, clusterSize, clusterIndices} and the second is the input array with cluster numbers replacing filled sites.

I've run this on a 1000x1000 array with filling factor ~0.5, and although it is about 3 times slower than the first method, it does work.

I'd be interested to know if anybody can think of a more clever way of implementing the algorithm. Mathematica has the SparseArray object, which is a good fit for the input and output configurations we are dealing with here, so one might consider using those from the beginning for the input. Thanks to ciao, who pointed out that one can extract the positions of non-zero elements easily from a SparseArray with ["NonzeroElements"], we can restrict our attention to those only.

• Thank you! This is very beautiful, and I will certainly take the time to pick this apart to learn more about Mathematica's image processing capability, but yes, it's not what I'm looking for. I'm mainly interested in mimicking the recursive, depth-first search algorithm using functional programming (partly as a learning tool, partly for speed improvements, if possible). Also, you're right: larger grids and larger filling fractions are certainly slower. Still, +1! Jun 4 '15 at 21:36
• Actually: during some testing (which perhaps I will post as an update), your method runs significantly faster than mine, as long as I remove the part of the code that finds the cluster indices (the clear rate-limiting step). For a 400 by 400 lattice with filling fraction of 0.593, yours is already 100 times faster (0.01 s versus 1 s). Mine dies already for a 500 by 500 lattice, and for a 10000 by 10000 lattice, yours takes 10 s. I wish I had asked two questions, one for speed and one for functional programming, so that I could accept your answer. Jun 5 '15 at 3:37
• @Virgil: cindices = {First@#, SparseArray[Unitize@Subtract[output, First@#], Automatic, 1][ "NonzeroPositions"]} & /@ csizes should give a nice speed bump for getting indices...
– ciao
Jun 5 '15 at 20:02
• @ciao: Thanks! I don't have too much experience with SparseArrays, so it is nice to know that one can get the non-zero elements. It turns out that this is a little faster: cindices = Union@Flatten[#, 1] & /@ GatherBy[{output[[Sequence @@ #]], #} & /@ SparseArray[output]["NonzeroPositions"], First] Jun 5 '15 at 21:10
• @Virgil: cindices = {#[[1, 1]], #[[All, 2]]} & /@ Rest[SortBy[ GatherBy[ Transpose[{Flatten@output, Tuples@Range@Dimensions[output]}], First], First]] about twice as fast...
– ciao
Jun 5 '15 at 21:29

March - this is not a complete answer, but instead addresses the gathering of indices since that seems to be important. Here's a comparison of some methods, tested only on my cigar-lounge netbook so caveat lector.

OP - Virgil's original (position based), Virgil - Virgil's adaptation of my first comment, CMT. 2/3 - my second and third comments (gathered tuples and pure sparse array), post - the code below inserted into Virgil's first snippet.

I stopped timing on the OP method at 200x200 array - it was nearing a minute on my poor netbook (and "post" was ~600X faster at that data point)...

Fill was ~0.5 for these tests. Times are for complete results from Virgil's code, just with index gathering changed.

"post" code:

cindices = Module[{sa = SparseArray[output], xx, yy, sa1, sa2},
sa1 = sa["NonzeroValues"];
xx = GatherBy[Range@Length@sa1, sa1[[#]] &];
sa2 = sa["NonzeroPositions"];
yy = sa2[[#]] & /@ xx;
Transpose[{sa1[[xx[[All, 1]]]], yy}]];


Benchies: I hope you find use in this.

• Nice! I was in the midst of this. So this is @Virgil's image-processing version with various ways of finding the cluster indices, is that right? 1000 x 1000 in 1 s with sizes, indices, and clustered grid as output is great! Jun 6 '15 at 6:42
• @march: Yes - I clarified post (Meant OP was Virgil's OP, and that times are for complete results, just changing the index gathering part). I'll ponder this further, the "post" idea came to me while cigar smoking, I think it can be improved. Also, no need to convert to image, MC works on arrays, might be worth timing (though I'll not be surprised if converting to image is faster even with that overhead - images are treated as atomic - but worth testing)
– ciao
Jun 6 '15 at 6:52
• This gathering of the site indices can be very useful, so thanks for giving it some thought! Jun 6 '15 at 6:55
• @march: No worries - interesting (and well presented) question - also see update to my earlier comment re: conversion to image...
– ciao
Jun 6 '15 at 6:56
• Nice! I was going to to this, too. It is an interesting question in itself. +1! Jun 6 '15 at 13:14

All done by Morphological Components

clusteringb[config_] := Module[{output, cm, cindices, csizes},
output = MorphologicalComponents[Image@Abs@config, CornerNeighbors -> False];
cm = ComponentMeasurements[ output, {"Label", "Mask", "Count"}][[All, 2]];
{cindices, csizes} = Transpose[{{#1, #2["NonzeroPositions"]}, {#1, #3}} & @@@ cm];
{csizes, cindices, output}];

inputConfig = {{-1, -1, -1, 0}, {-1, 0, 0, -1}, {-1, -1,  0, -1}, {0, -1, 0, 0}};
clusteringb@inputConfig

(* {{{1, 7}, {2, 2}},
{{1, {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {3, 1}, {3, 2}, {4, 2}}},
{2, {{2, 4}, {3, 4}}}},
{{1, 1, 1, 0}, {1, 0, 0, 2}, {1, 1, 0, 2}, {0, 1, 0, 0}}}

• Thank you for your answer! Is the significant difference between this and @Virgil/ciao's answer that you have used ComponentMeasurements? Time for more benchmarking! (I am still interested in seeing people's take on the recursive depth-first search, Mathematica-style; perhaps the question is ill-posed, or perhaps MorphologicalComponents is Mathematica's built-in depth-first clustering algorithm?) Jun 10 '15 at 16:48
• @march The main problem when programming in Mma is finding a way without programming. Mma internal algorithms are usually faster and easier to use than any user-programmed module. You may try the already done depth-first search but I believe the setup is too cumbersome Jun 10 '15 at 17:02
• @beliarius Sure: agreed. I've just benefited greatly from seeing the experts come up with clever (even if not fast!) ways of solving problems, even if there is a built-in option. Jun 10 '15 at 17:19