Return to Answer

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.

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.

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.

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.

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} + # & /@ {{0, 1}, {0, -1}, {1, 0}, {-1, 0}},
       (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}
  ];

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}
  ];

clustering1[config_] := Module[{output, csizes, cindices},
   output =  MorphologicalComponents[Image@Abs@config, CornerNeighbors -> False];
   csizes = Rest@Sort@Tally@Flatten@output;
   cindices = Map[
     Union@Flatten[#, 1] &,
     GatherBy[
      {output[[Sequence @@ #]], #} & /@ SparseArray[output]["NonzeroPositions"], First]
    ];
   {csizes, cindices, output}
  ];

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, I've replaced it with a faster construction.

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}
  ];

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.