Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question

2 Answers 2

up vote 5 down vote accepted

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 *)
share|improve this answer

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

mat1

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

mat2

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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.