Use of Cellular Automaton or Discrete Convolution to calculate Potts Energy of a 2D Matrix

Question

In a project of Monte Carlo Modelling of Grain Growth (Link) one of the most costly subroutines is the calculation of the Potts-Energy of an Monte Carlo Grid cell by comparing the center cell value with it's 9 neighbours and calculating a resulting Potential energy which is calculated as

pottsenergy[osub_]:=Total[Map[1-KroneckerDelta[#,osub[[2,2]]]&,osub,{2}],2];

with osub as a 3x3 matrix of integers. The overall algorithm convolves this calculation across the whole system matrix oin which can consist of big 2D matrices in the following way:

energy[o_, indx_, nmax_] := pottsenergy[o[[Sequence @@ energypart[indx, nmax]]]];

with energypart as a routine which returns the boundaries of the center cell position denoted by indx and its 9 neighbours in part notation:

energypart[indx_, {nx_, ny_}] :=Module[{ii, jj, it, ib, jl,jr},(*get part specification of nearest neighbours of location (ii,jj). Reflecting boundary conditions are assumed.*)
ii = indx[[1]]; jj = indx[[2]];
(*rows of top and bottom neighbours*)
it = If[ii == 1, 1, ii - 1];
ib = If[ii == ny, ny, ii + 1];
(*rows of left and right neighbours*)
(*rows of top and bottom neighbours*)
jl = If[jj == 1, 1, jj - 1];
jr = If[jj == nx, nx, jj + 1];
{it ;; ib, jl ;; jr}
];

The actual convolution is then done as

allenergy[oin_, nmax_] := MapIndexed[energy[oin, #2, nmax] &, oin, {2}]

which takes a big 2D matrix and returns a big matrix with the Potts energy values of the input matrix oin:

oin = RandomInteger[31, {100, 100}];
allenergy[oin, Dimensions[oin]] // MatrixForm

Now the main problem is that this implementation is rather slow and cannot be easily sped up by use of a CUDA or parallelize implementation. Looking into the Mathematica documentation I found some promising algorithms like Cellular Automaton or Discrete Convolution which would be definitely faster I suppose. Now the use of these routines are unfortunately non-trivial for above way of folding.

Any ideas how to speed up above algorithm by using the power of Mathematica?

Answer 1

I'm not sure you are handling the top and left edges in the way you really want; they work with the second rather than first elements being treated as "middle" twice, with first elements not treated that way at all.

Here is code that does not do that, hence gives different results than yours on top and left edges. It is around two orders of magnitude faster. I will assume that the minimal value is positive but this can be adjusted if that's not the case. Here are the steps.

(1) Pad on all sides by zero.

(2) Partition with overlaps of length 2 in both dimensions. This gives neighborhood submatrices.

(3) Subtract "center" elements from each submatrix.

(4) Take Sign[Abs[]] of these differences. But we want there to be no contribution from edges that came from padding, so we also multiply by the submatrix to zero those terms (our assumption enforces that the only elements that are zero in the neighborhood matrices are those that came from the padded edge boundaries).

Here is the code, in its entirety.

energy2[mat_] := Module[
  {n = Length[mat], 
   pmat = Partition[ArrayPad[mat, 1, 0], {3, 3}, 1]},
  Map[Total[Flatten[#]] &, Abs[Sign[pmat*(pmat - mat)]], {2}]
  ]

Example:

SeedRandom[11111];
n = 6;
mat = RandomInteger[{1, 5}, {n, n}]

(* Out[486]= {{3, 1, 1, 3, 1, 2}, {4, 4, 3, 3, 3, 4}, {1, 5, 4, 2, 5, 
  3}, {3, 3, 4, 1, 3, 1}, {4, 3, 3, 3, 5, 5}, {4, 1, 1, 3, 4, 4}} *)

e1 = allenergy[mat, Dimensions[mat]]

(* Out[489]= {{2, 4, 3, 2, 3, 3}, {4, 6, 6, 5, 5, 5}, {5, 8, 6, 8, 8, 
  3}, {3, 5, 7, 8, 6, 5}, {3, 5, 4, 5, 7, 4}, {3, 4, 4, 3, 4, 2}} *)

e2 = energy2[mat]

(* Out[490]= {{3, 4, 4, 2, 5, 3}, {4, 6, 6, 5, 5, 5}, {5, 8, 6, 8, 8, 
  3}, {3, 5, 7, 8, 6, 5}, {4, 5, 4, 5, 7, 4}, {2, 4, 4, 3, 4, 2}} *)

Notice that, as claimed, these only differ on top and left edges.

e1 - e2

(* Out[491]= {{-1, 0, -1, 0, -2, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
  0}, {0, 0, 0, 0, 0, 0}, {-1, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0}} *)

To indicate relative speeds, here is a bigger example.

SeedRandom[11111];
n = 1000;
mat = RandomInteger[{1, 5}, {n, n}];

Timing[e1 = allenergy[mat, Dimensions[mat]];]

(* {59.760000, Null} *)

Timing[e2 = energy2[mat];]

(* {0.520000, Null} *)

Check for agreement of top and left edges:

Max[Abs[(e1 - e2)[[2 ;; -1, 2 ;; -1]]]]

(* Out[497]= 0 *)

Answer 2

Before I start, I want to acknowledge that Daniel's answer is faster than mine because mine does not take advantage of the specific form of the energy function. However, the solution with Cellular Automaton seems very cool and works for any energy function so I decided to put it up anyway.

CellularAutomaton can take a function as an argument to evaluate on the group of neighbors, so let's optimize that.

Notice that

pottsenergy[osub_] := Total[Map[1 - KroneckerDelta[#, osub[[2, 2]]] &, osub, {2}], 2]

Is equivalent to

alternativeEnergy[osub_] := Count[osub, Except[osub[[2, 2]]], {2}]

However, the latter is faster:

test = Array[ RandomInteger[{0, 5}, {3, 3}] &, 100000];
pottsenergy /@ test; // AbsoluteTiming
alternativeEnergy /@ test; // AbsoluteTiming
(* {4.044231, Null} *)
(* {0.609035, Null} *)

You did not specify what to do with the border elements, since your construct did not handle them properly, as Daniel pointed out.

I will assume that what you intend is for elements outside the border not to contribute to the Potts-Energy.

To do this, we have to assign a secretvalue to outside-the-border elements.

As such, we have to edit our alternativeEnergy:

secretvalue = -51.3;
alternativeEnergy[osub_] := Count[osub, Except[osub[[2, 2]] | secretvalue], {2}]

Note that secretvalue remains an unevaluated symbol for now. Performance is more or less the same, but can be sped up slightly by letting the secretvalue be something that doesn't show up in the actual matrix (i.e., you can use 0 like Daniel or -51.3).

Then, it's just a matter of knowing how to use CellularAutomaton:

SeedRandom[11111];
n = 6;
(mat = RandomInteger[{1, 5}, {n, n}]) // MatrixForm

Note that CellularAutomaton has an annoying habit of displaying the padded borders, so I clear them out.

ca[mat_] := 
    ArrayPad[
        CellularAutomaton[{alternativeEnergy[#1]&,{},{1,1}},{mat,secretvalue}]//First,
        -1]
ca[mat]//MatrixForm

SeedRandom[11111];
n = 1000;
mat = RandomInteger[{1, 5}, {n, n}];
AbsoluteTiming[ca[mat];]
(* {7.090406, Null} *)