Needs["Graphics`Animation`"];

numAgents = 21;
roomSize = 10;
idealMin = 3;
idealMax = 4;

insertInitialAgent[agent_, pos_] := initialSociety = 
        ReplacePart[initialSociety, agent, pos]
initialSociety = 
        Table[-1, {3*roomSize}, {3*
              roomSize}];
initialSociety = 
        ReplacePart[initialSociety, 0, 
          Flatten[Table[{i, j}, {i, roomSize + 1, 2*roomSize + 1}, {j, 
                roomSize + 1, 2*roomSize + 1}], 1]];
For[i = 1, 
      i ≤ 
        numAgents, i++, 
While[randPos = {Random[
              Integer, {roomSize + 1, 2*roomSize + 1}], 
            Random[Integer, {roomSize + 1, 2*roomSize + 1}]}; 
        Extract[initialSociety, randPos] =!= 
          0,]; 
insertInitialAgent[
        agent[Random[Integer, {18, 80}], Random[Integer]], 
        randPos];
]

distFromRange[x_, {min_, max_}] := 
      Which[x ≤ min, min - x, x ≥ max, x - max, 
        True, 0]

clusterUtil[agent[myAge_, mySex_], agents_] := 
      Module[{size = Length[agents],  averageAge, 
          averageSex}, 
averageAge = 
          Apply[Plus,  Map[First,  agents]]/
            If[size > 0, size, 1] ; 
averageSex = 
          Apply[Plus,  Map[#[[2]] &,  agents]]/
            If[size > 0, size, 1] ; 
ageDiff = 
          Abs[myAge - 
              averageAge]; 
-2*(ageDiff - (80 - 
                      18))/(80 - 18) - 
          Log[1 + distFromRange[size, {idealMin, idealMax}]]]

(* 
      determining the age difference at which an agent would just as soon be alone,  
      even though the cluster is of ideal size *) 
Solve[-2*(d - (80 - 18))/(80 - 
                18) == 2 - Log[3. ], d]

neighborhood[cell_, radius_] := 
      Map[Function[x, cell + x], 
        Flatten[Table[{i, j}, {i, -radius, radius}, {j, -radius, 
              radius}], 1]]

occupiedCells[society_, cell_, 
        radius_] := 
      Map[#[[2]] &, 
        Cases[Map[{Extract[society, #], #} &, 
            neighborhood[cell, radius]], {_agent, _}]]

emptyCells[society_, cell_, radius_] := 
      Map[#[[2]] &, 
        Cases[Map[{Extract[society, #], #} &, 
            neighborhood[cell, radius]], {0, _}]]

octopus[society_, cells_] := 
      Apply[Union, Map[occupiedCells[society, #, 2] &, cells]]

spreadInk[society_, cell_] := 
      FixedPoint[octopus[society, #] &, {cell}]

agentUtil[me_, society_, cell_] := 
      If[Extract[society, cell] === -1, -Infinity, 
        clusterUtil[me, spreadInk[society, cell]]]

randomize[{}]  :=  {};
\(\(randomize[
          lis_]\  := 
	\(Transpose[
            Sort[Map[{Random[], \ #} &, \ lis]]]\)[\([2]\)];\)\)

(* 
      given a state of the room (society) and the location of an agent,  
      return a new state where that agent has moved to cell of highest utility within a radius of where it' 
        s currently at *) 
update[society_, 
        cell_, radius_] := 
      Module[{me = Extract[society, cell], s = society, possibleBests, 
          bestcell}, 
If[! MatchQ[
              Extract[society, cell], _agent], 
          Return[society]]; 
s = 
          ReplacePart[s, 0, 
            cell]; 
possibleBests = Transpose[Split[
                  Sort[Map[{agentUtil[me, s, #], #} &, 
                      emptyCells[s, cell, 
                        radius]], #1[[1]] > #2[[1]] &], #1[[1]] == #2[[1]] &][[1]]][[2]]; 
bestcell = 
          possibleBests[[Random[
                Integer, {1, 
                  Length[possibleBests]}]]]; 
ReplacePart[
          s, Extract[society, cell], 
          bestcell]
]

updateEach[society_, 
        radius_] := 
      FoldList[update[#1, #2, radius] &, society, 
        randomize[Position[society, _agent]]]

\(imageList[societies_] := 
      Map[Graphics[
            RasterArray[#\  /. 
							{0 -> 
                    RGBColor[0, 0, 1], 
								\(-1\) -> 
                    RGBColor[0, 1, 0], 
								agent[age_, 
                      sex_] -> 
                    Apply[RGBColor, 
                      Table[\((age - 18)\)/\((80 - 
                              18)\), \ {3}]]}]] &, 
			societies]\)

showFrames[societies_] := 
      ShowAnimation[imageList[societies]]

run[n_, radius_] := 
      Nest[Join[#, updateEach[Last[#], radius]] &, {initialSociety}, n] // 
        showFrames

runSparse[n_, radius_] := 
      NestList[Last[updateEach[#, radius]] &, initialSociety, n] // 
        showFrames

runSparse[30, 1]