3

I have a $w \times h$ array with either $-1$ or $+1$.

I wish to monitor the size of the largest connected clusters with specific patterns that occur during a Monte Carlo simulation, i.e.,

  • a collection of $+1$'s,
  • a checkerboard pattern of $-1$ and $+1$'s, and,
  • horizontal and vertical stripes of $-1$ and $+1$'s.

Moreover, periodic boundary conditions apply such that opposite boundaries are connected.


Example arrays (checkerboard, stripes) can be imported by,

Uncompress[Import["https://pastebin.com/raw/0GaUtiFy"]]
Uncompress[Import["https://pastebin.com/raw/uF5UagqU"]] 

Ideally I am looking for a fast solution, such that monitoring every iteration of the Monte Carlo simulation will not slow it down to a crawl.

The only idea I have had so far was to use ComponentMeasurements but was not able to specify the patterns. Any help would be greatly appreciated.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
user19218
  • 841
  • 4
  • 15

1 Answers1

2

This might not be the fasted method but it will get you started:

Preparation: The working horse is getComponent while periodicBoundaryCorrection only accounts for the fact that MorphologicalComponents does not allow for periodic boundary conditions. See the also here for a more detailed description.

checkermask = Array[Mod[Plus[##], 2] &, {3, 3}];
hormask = Array[Mod[#1, 2] &, {3, 3}];
vermask = Array[Mod[#2, 2] &, {3, 3}];

ClearAll[periodicBoundaryCorrection];
periodicBoundaryCorrection[
 A_?MatrixQ, 
 OptionsPattern[{CornerNeighbors -> True}]
 ] := Module[{maxA, a, b, c, d, pos, r1, r2, r3, r4, r5, r6, edges, α, β, colorcomp, cols},
   maxA = Max[A];

   a = A[[1]];
   b = A[[-1]];
   pos = DeleteCases[Range[Length[a]] Unitize[a b], 0];
   r1 = Sort /@ Transpose[{a[[pos]], b[[pos]]}];

   c = A[[All, 1]];
   d = A[[All, -1]];
   pos = DeleteCases[Range[Length[c]] Unitize[c d], 0];
   r2 = Sort /@ Transpose[{c[[pos]], d[[pos]]}];

   edges = Union[r1, r2];

   If[OptionValue[CornerNeighbors],

    pos = DeleteCases[Range[Length[a] - 1] Unitize[Rest[a] Most[b]], 0];
    r1 = Sort /@ Transpose[{a[[pos + 1]], b[[pos]]}];
    pos = DeleteCases[Range[Length[b] - 1] Unitize[Rest[b] Most[a]], 0];
    r2 = Sort /@ Transpose[{b[[pos + 1]], a[[pos]]}];

    pos = DeleteCases[Range[Length[c] - 1] Unitize[Rest[c] Most[d]], 0];
    r3 = Sort /@ Transpose[{c[[pos + 1]], d[[pos]]}];
    pos = DeleteCases[Range[Length[d] - 1] Unitize[Rest[d] Most[c]], 0];
    r4 = Sort /@ Transpose[{d[[pos + 1]], c[[pos]]}];

    α = c[[1]];
    β = d[[-1]];
    r5 = If[α β != 0, {Sort[{α, β}]}, {}];
    α = c[[-1]];
    β = d[[1]];
    r6 = If[α β != 0, {Sort[{α, β}]}, {}];

    edges = Union[edges, r1, r2, r3, r4, r5, r6]
   ];
   edges++;
   If[Length[edges] == 0,
    A,
    colorcomp = SparseArray`StronglyConnectedComponents[
      SparseArray[
       Join[edges, Transpose[Transpose[edges][[{2, 1}]]]] -> 
        1, {maxA + 1, maxA + 1}, 0]
      ];
    cols = Compile[{{idx, _Integer, 1}, {acc, _Integer, 1}},
       Block[{colors, j, threshold},
        colors = Table[0, {i, 1, acc[[-1]]}];
        j = 0;
        threshold = Compile`GetElement[acc, j + 1];
        colors[[idx]] = Table[
          If[i > threshold,
           j++;
           threshold = Compile`GetElement[acc, j + 1];
           ];
          j, {i, 1, Length[idx]}];
        colors
        ]
       ][
      Join @@ colorcomp,
      Accumulate[Length /@ colorcomp]
      ];
    Compile[{{a, _Integer, 1}, {cols, _Integer, 1}},
      cols[[a + 1]],
      RuntimeAttributes -> {Listable},
      Parallelization -> True
      ][A, cols]
    ]
   ];

getComponent[A_?MatrixQ, mask_] := periodicBoundaryCorrection[
  MorphologicalComponents[
   With[{img = Binarize[Image[A], 0]},
    With[{
      a = ImageCorrelate[img, mask/Total[mask, 2], Padding -> "Periodic"],
      b = ImageCorrelate[img, (1 - mask)/Total[1 - mask, 2], 
        Padding -> "Periodic"]
      },
     With[{ab = ImageMultiply[a, b]},
      Binarize[ImageAdd[ImageSubtract[a, ab], ImageSubtract[b, ab]], 
       0.99]
      ]
     ]
    ]
   ]
  ]

Application:

A = Uncompress[Import["https://pastebin.com/raw/0GaUtiFy"]];
B = Uncompress[Import["https://pastebin.com/raw/uF5UagqU"]];

GraphicsGrid[{
  {
   Image[A],
   Colorize[getComponent[A, checkermask]],
   Colorize[getComponent[A, hormask]],
   Colorize[getComponent[A, vermask]]
   },
  {
   Image[B],
   Colorize[getComponent[B, checkermask]],
   Colorize[getComponent[B, hormask]],
   Colorize[getComponent[B, vermask]]
   }
  }, ImageSize -> Large]

enter image description here

Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309
  • Thank you, this definitely got me started. Reverse engineering your solution I found that getComponent[A_, mask_] := With[{B = -A ListCorrelate[mask, A, {2, 2}]}, ImageData[Binarize[Image[B], 3]]] almost does the trick. The only thing left is to implement the boundary conditions – user19218 Apr 11 '18 at 19:30
  • Ah sorry, I missed that. You have to set Padding - > "Periodic" in ListCorrelate. At least that should make the input for MorphologicalComponents correct. MorphologicalComponents won't match opposing components correctly. It has a Padding options but that does not allow for "Periodic"... – Henrik Schumacher Apr 11 '18 at 19:32
  • No sorry I will clarify, the components should found with periodic boundary conditions, i.e. in the checkerboard example there should only be one cluster found, instead of four – user19218 Apr 11 '18 at 19:40
  • In any case this question already has a answer – user19218 Apr 11 '18 at 19:52
  • Thanks for the hint. It was too late; already added a fix. – Henrik Schumacher Apr 11 '18 at 20:19