Mathematica Programming Challenge

The Results!

Judged by John M. Novak

The Challenge: Aggregation!

Simulate a simple model of aggregation on a 2D grid: start with a single black cell, and at each step fill in a random white neighbor of a randomly selected existing black cell.

The Details

This is a simple model of growth sometimes called the Eden model (specifically, a type B Eden model). The model has been used since the 1960s as a model for such things as tumor growth and growth of cities. Here's how it works. Have a 2D grid of black and white cells (1s and 0s). Start with a single black cell. At each step, select at random an existing black cell that has at least one white neighbor.  Then select a neighboring white cell at random, and make it black. (By neighbor we mean a cell above, below, left or right—not diagonal.) Your function should be Aggregate[state, t]; state should be an array of 0s and 1s, t should be the number of steps, and the result should be another array of 0s and 1s, with exactly t new black cells. Your function should allow the cluster of black cells to get as big as it wants (i.e., the array may grow). Unfilled cells should always be assumed white. See if you can generate a really big cluster. Make a picture of it, using ListDensityPlot. Is it a circle?

If you are feeling adventurous, you could try changing the criterion for growth, and seeing what shapes you can get.  For example, allow growth only when the neighbors around a cell have some particular form. What shapes can you make?

The Judging

The submissions were first examined for elegance, efficiency, and meeting the problem statement. In many cases, a program may demonstrate a rightness, or artistic elegance, in solving a particular task; this is what we ideally strive in the Mathematica programming challenges. Efficiency is not to be sneered at, however; for a problem like this, we would really like to see a linear-time solution.

As it happened, none of the submissions showed a clear advantage over any of the others; some were slightly better in one way, some in another, but none extremely so. As a result, the focus of the judging was broadened to the less definable “Spirit of Mathematica.” This asks that the code take special advantage of Mathematica's environment and philosophy.

Applying these parameters resulted in several submissions that merited special mention in having one characteristic or another, plus a winning entry chosen on the basis of making the best use of Mathematica.

The Solutions: Honorable Mention

Richard Mercer, for the cleanest, most modular code.

Richard's code, pulled out of a very well-documented notebook, is nice, clean, and modular. It does some redundant operations, and isn't the fastest possible approach, but gains kudos for being the easiest for a less-experienced user to read, great for pedagogic purposes. The full notebook included detailed descriptions of the various subfunctions, which enhanced this aspect considerably.

AggregateOne[mat_] := 
	Module[{base,newmat=mat},
		If[NeedRoom[mat],newmat = GrowMat[mat]]; 
		base = NewBase[newmat];
    ReplacePart[newmat,1,base + 
		NewDirection[newmat,base]]
	];

Aggregate[mat_, t_] := Nest[AggregateOne, mat, t];

NeedRoom[mat_]:= 
	Union[Flatten[{
		Take[mat,1],
		Take[mat,-1],
		Take[Transpose[mat],1],
		Take[Transpose[mat],-1]}
	]] =!= {0};

AddZeros[mat_]:= Prepend[Append[#,0],0]& /@ mat; 
GrowMat[mat_]:=
	Transpose[AddZeros[Transpose[AddZeros[mat]]]];

Adjacents[mat_, base_] := 
	Part[mat,##]& @@@ ((base + #)& /@ dirs);

RandomFromList[lst_] := 
  Part[lst,Random[Integer,{1,Length[lst]}]];

NewBase[mat_] := 
	Module[{nospace=True,base},
		While[nospace,
			base = RandomFromList[Position[mat,1,{-1}]];
			nospace = (Adjacents[mat,base]=={1,1,1,1})
		];
		base
	];

dirs = {{-1,0},{0,-1},{0,1},{1,0}};

NewDirection[mat_,base_] := 
	Part[dirs,RandomFromList[ 
		Position[Adjacents[mat,base],0,{-1}]]
	] //First

Seth Chandler for the most concise code.

Seth Chandler created the most compact code of the submitted entries; while it does a lot of data manipulation which costs it in terms of efficiency, it does use clever Mathematica programming to keep the command-count down! Seth also makes use of NestList to generate the steps, allowing us to just map a plotting command across the output to get an animation.

Aggregate[state_, t_] := 
      Block[{onestep, 
          expandrules}, 
expandrules = {{x : {___, 1, ___}, 
                rest__} :> {Table[0, {Length[x]}], x, rest}, {start__, 
                x : {___, 1, ___}} :> {start, x, 
                Table[0, {Length[x]}]}}; 
onestep[astate_] := 
          Module[{newstate, pm, bc}, 
newstate = 
              With[{s = astate //. expandrules}, 
                Transpose[
                  Transpose[s] //. expandrules]]; 
pm = 
              Partition[newstate, {3, 3}, {1, 1}, {{2, 2}, {2, 2}}, 
                0]; 
wc = 
              Map[{#1, Count[Extract[pm, #], 1, {2}]} &, 
                Position[pm, {_, {_, 0, _}, _}, {2}], 
                1]; 
ReplacePart[newstate, 
              1, Extract[wc, 
                  With[{lwc = Last /@ wc}, 
                    With[{r = Random[Integer, {1, Plus @@ lwc}]}, 
                      Position[
                        Partition[FoldList[Plus, 0, lwc], 2, 
                          1], {_?(# < 
                                  r &),  _?(r <= # &)}]]]][[1, 
                  1]]]]; 
NestList[onestep, state, 
          t]]

David DeBrota for the most efficient code.

David DeBrota's code uses clever manipulation of data structures to create a fairly efficient approach, which turned in timings several times better than other competitors. However, the end result of this approach is that the code goes straight out and randomly chooses a fillable boundary cell, instead of first selecting a black cell, then randomly choosing a fillable neighboring cell. This variant is called the 'type A' Eden model, as opposed to the 'type B' model posed in problem. This code also works inside a fixed grid size, and does not grow the grid as the aggregate outgrows its bounds.

Aggregate[state_, t_] := 
      Module[{x, ii, m, n, c, fillable}, 
m = 
          Length[state]; 
        n = Length[state[[1]]]; 
x = 
          Map[Prepend[Append[#, 0], 0] &, 
            Prepend[Append[state, Table[0, {n}]], 
              Table[0, {n}]]]; 
fillable = 
          Select[Union[
              Flatten[Map[{{#[[1]] - 
                          1, #[[2]]}, {#[[1]], #[[2]] - 
                          1}, {#[[1]] + 
                          1, #[[2]]}, {#[[1]], #[[2]] + 1}} &, 
                  Position[x, 1]], 
                1]], (#[[1]] ≠ 
                    1  &&  #[[1]] ≠ 
                    m + 2  &&  #[[2]] ≠ 
                    1  &&  #[[2]] ≠ 
                    n + 2) &]; 
For[ii = 1, 
          ii ≤ t, ii++; 
          If[Length[fillable] > 0, 
            c = fillable[[Floor[Random[]*Length[fillable]] + 
                    1]]; 
x[[c[[1]], 
                  c[[2]]]] = 1; 
fillable = 
              Drop[fillable, 
                Flatten[
                  Position[fillable, {c[[1]], c[[2]]}, 1, 
                    1]]]; 
If[
              c[[1]] != 2  &&  
                x[[c[[1]] - 1, c[[2]]]] == 0, 
              fillable = 
                Append[fillable, {c[[1]] - 1, 
                    c[[2]]}]]; 
If[
              c[[1]] ≠ m + 1  &&  
                x[[c[[1]] + 1, c[[2]]]] == 0, 
              fillable = Append[fillable, {c[[1]] + 1, c[[2]]}]]; 
            If[c[[2]] != 2  &&  
                x[[c[[1]], c[[2]] - 1]] == 0, 
              fillable = Append[fillable, {c[[1]], c[[2]] - 1}]]; 
            If[c[[2]] ≠ n + 1  &&  
                x[[c[[1]], c[[2]] + 1]] == 0, 
              fillable = 
                Append[fillable, {c[[1]], 
                    c[[2]] + 1}]]; 
fillable = 
              Union[fillable]
]]; 
Map[
          Take[#, {2, m + 1}] &, 
          Take[x, {2, n + 1}]]
]

This code was turning in timings of just over a second for a thousand generations, approximately an order of magnitude better than most competitors… and the time seemed to scale reasonably well.

Olivier Gerard for an all-around promising entry.

This entry was unfortunately received after the deadline, and there was not time to analyze it in detail. However, it shows both clever use of Mathematica programming and appears to be competitively efficient.

FixArrayRight[kk_List]  :=  
      Module[{m, M, 
          r}, 
{m, 
            M}  =  (Length /@ kk) // {Min@#, 
                Max@#} &; 
r = 
          Array[0 &, M - m]; 
Map[Take[Join[#, r], M] &,  
          kk]]

FixArrayLeft[kk_List]  :=  
      Module[{m, M, 
          r}, 
{m, 
            M}  =  (Length /@ kk) // {Min@#, 
                Max@#} &; 
r = 
          Array[0 &, M - m]; 
Map[
          Take[Join[r, #], -M] &,  kk]]

Aggregate[state_List]  :=  
      Module[ {base, hvf,   nl,  cl, tl, zp, 
          zi}, 
hvf = { Transpose,  
              Identity}[[Random[
              Integer,  {1, 2}]]]; 
base = 
          hvf[state]; 
nl = 
          Flatten@Position[
              Plus @@@ base, _?Positive]; 
cl = 
          Flatten[Split[
                Join[{0}, 
                  base[[tl = 
                      nl[[Random[
                          Integer, {1, 
                            Length[nl]}]]]]], {0}]] //. {{{0, 
                      d : ___, 0}, r : __List} -> {{0, d, 2}, 
                    r}, {q : __List, {0, d : ___, 0}} -> {q, {2, d, 
                      0}}, {0, d : ___, 0} -> {2, d, 
                    2}, {0} -> {2}}]; 
zp = 
          Flatten[Position[cl, 2]]  /.  {1 -> 
                0, ( Length[base[[tl]]] + 
                    2) -> -1}; 
zi = 
          zp[[Random[
              Integer, {1, Length[zp]}]]]; 
hvf[
          Switch[zi, 0, 
            FixArrayLeft@
              ReplacePart[base, Prepend[base[[tl]], 1], 
                tl], -1, 
            FixArrayRight@
              ReplacePart[base, Append[base[[tl]], 1], tl], _, 
            ReplacePart[base, 
              1, {tl, zi - 1}]]]
]

Aggregate[state_List, t_Integer] := 
      Nest[Aggregate, state, t];

The Winner

Robert Nachbar for the most general and "Mathematica-like" code!

Effective use of Mathematica depends on a willingness to experiment; generalization is a good way to do that. While the general approach is technically an optional part of this competition, Robert is the only entrant who actually tackled this, resulting in pattern-based code that employs a good variety of Mathematica functionality. It does suffer from some efficiency issues, but makes up for that in originality. It also makes use of typesetting functionality, and uses various message, option handling, and error checking constructs, as a good application should.

growState[
        m_] := m /. {f : {___, 1, ___}, 
                  r___} :> {Table[0, {Length[f]}], f, 
                  r} /. {r___, l : {___, 1, ___}} :> {r, l, 
                Table[0, {Length[l]}]} /. 
          a : {___, {1, ___}, ___} :> Map[Prepend[#, 0] &, a] /. 
        a : {___, {___, 1}, ___} :> 
          Map[Append[#, 0] &, a]

neighborhood[
        p_] := p + # & /@ {{-1, 0}, {+1, 0}, {0, -1}, {0, 
            1}}

selectOne[l_] := 
      l[[Random[
          Integer, {1, 
            Length[l]}]]]
selectOne[{}] := $Failed

accrete::nopat = "Could not find cell with desired \
neighborhood."

accrete[state_, pat_] := 
      Module[{map = growState[state], blackCell, whiteCell}, 
        blackCell = selectOne[Position[Partition[map, {3, 3}, {1, 1}], pat]]; 
        If[blackCell === $Failed, Message[accrete::nopat]; Throw[map], 
          blackCell += {1, 1}]; 
        whiteCell = 
          Select[neighborhood[blackCell], 
              Part[map, Sequence @@ #] == 0 &] // selectOne; 
        map[[Sequence @@ whiteCell]] = 
          1; map]

ClearAll[Aggregate]
Options[Aggregate]={GrowthPattern->(({{_, 0, _}, {_, 1, _}, {_, _, _}}))|(({{_, _, _}, {_, 1, _}, {_, 0, _}}))|(({{_, _, _}, {0, 1, _}, {_, _, _}}))|(({{_, _, _}, {_, 1, 0}, {_, _, _}}))};
Aggregate[state_, 
        t_Integer?NonNegative, options___?OptionQ] := 
      Module[{pat}, 
        pat = GrowthPattern /. {options} /. Options[Aggregate]; 
        Catch[Nest[accrete[#, pat] &, state, t]]]

Here's a demonstration of the code using an alternate set of rules for growth; observe that the growth no longer has the roughly circular form of the original growth pattern.

init  =  PadLeft[{{1}},  {10, 10}, 0, {5, 5}];

ListDensityPlot[1-Aggregate[init,100,GrowthPattern->(({{_, 0, _}, {0, 1, 0}, {_, 0, _}}))|(({{_, 0, _}, {0, 1, 1}, {_, 0, _}}))|(({{_, 1, _}, {0, 1, 0}, {_, 0, _}}))|(({{_, 1, _}, {0, 1, 1}, {_, 0, _}}))|(({{_, 1, _}, {0, 1, 1}, {_, 1, _}}))], Mesh -> False]

The Wolfram Research Way

Here's an alternate formulation devised by Wolfram Research prior to the contest. (This was originally presented at the conference as belonging to Oyvind Tafjord; after the presentation, it was determined that the core algorithm was written by Stephen Wolfram. John Novak threw the grid-to-state conversion on top of the vergy elegant core algorithm to match the terms of the problem as stated.):

Aggregate[grid_,  t_]  :=  
      statetogrid[Nest[AMStep,  gridtostate[grid],  t]]

AMStep[
        c_]  :=  (If[! MemberQ[c,  #],  Append[c,  #],  
              AMStep[c]] &)[
        f[c] + f[{{1, 0}, {0, 1}, {-1, 0}, {0, -1}}]]

f[a_]  :=  a[[Random[Integer,  {1,  Length[a]}]]]

gridtostate[g_]  :=  Position[g,  1]

statetogrid[
        s_]  := 
Module[{ minx,  maxx,  miny,  
          maxy}, 
{{miny,  maxy}, {minx,  
              maxx}}  =  {Min[#],  Max[#]} &  /@  
            Transpose[s]; 
ReplacePart[
          Table[0,  {maxy  -  miny  +  1},  {maxx  -  minx +  1}],  
          1,  ((#  -  {miny,  minx}) + 1  &)  /@  
            s]]

This approach certainly has advantages in terms of elegance, though it isn't the most efficient.