Cellular automata

Start with the game of life, a cellular automata on a 2D board where live cells remain alive if they have 2 or 3 neighbours, dead cells become alive if they have 3 neighbours (make it easy to change these parameters since interesting behaviour is obtained for different rules). The basic algorithm is:

  1. Start with a set of live cells and calculate the set of neighbouring cells.
  2. Remove the set of live cells from the neighboring cells set, this results in the set of dead cells that can be influenced by the live cells.
  3. For each cell in the live cell set determine the size of the intersection of the set containing its neighbours and the set of live cells. If this number is within the set of numbers for a live cell to remain alive then add the cell to the new set of live cells.
  4. As 3 but for the dead cells so the rule is different.
  5. Now have the set of live cells in the next generation.

%%% Code: life.erl

-author('Matt McDonnell').

% liferule returns true if the cell is alive in the next generation,
% false otherwise.  State is the current cell state,  Neighbours is
% number of live neighbours and LiveNums is a list of number of
% neighbours required for a live cell to live (similarly for DeadNums)
liferule(State, Neighbours, {LiveNums, DeadNums}) -> 
  case State of
    alive_cell -> lists:member(Neighbours,LiveNums);
    dead_cell  -> lists:member(Neighbours,DeadNums)

liferule(State,NumLive) -> liferule(State,NumLive,{[2,3],[3]}).

% neighbours returns a list of coordinates of neighbouring cells
neighbours({X,Y},NeighbourDirs) ->
  lists:map(fun ({X0,Y0}) -> {X0+X, Y0+Y} end, NeighbourDirs).

neighbours({X,Y}) -> 
  neighbours({X,Y}, [{-1,-1},{-1,0},{-1,1},{0,1},{1,1},{1,0},{1,-1},{0,-1}]).

% new_live takes a set of cells all of a given state and returns the
% set of cells that live in the next generation.
new_live(State, Coords, LiveCoords, Acc0) ->
    sets:fold(fun(E,Acc) -> 
       NumLive = num_live(E, LiveCoords),
       Alive = liferule(State, NumLive), 
         Alive == true -> sets:add_element(E,Acc);
         true -> Acc
     end, Acc0, Coords). 

% num_live = number of living cells surrounding the given cell.
num_live(Coord, LiveCoords) ->
  sets:size( sets:intersection(sets:from_list(neighbours(Coord)),

server() ->
    quit -> 
      io:format("Life server exiting normally.~n",[]),
    {PID, {update, Cells}} ->
      NewCells = update(Cells),
      PID ! {life, NewCells},

update(LiveCoords) -> 
  % Cells is a set of {x_coord, y_coord} coordinates of live cells
  NeighbourCoords = sets:from_list(
    sets:fold(fun(E, Acc) -> neighbours(E) ++ Acc end, [], LiveCoords)),
  DeadCoords = sets:subtract(NeighbourCoords, LiveCoords),
  NextLive = new_live(alive_cell, LiveCoords, LiveCoords, sets:new()),
  new_live(dead_cell, DeadCoords, LiveCoords, NextLive).

%%% End Code: life.erl

The above calculates the next generation of live cells from the exisiting generation. This is the heart of the program, but input and output are also required, ie some gui that displays the state of the automaton and allows input of initial patterns, speed of simulation, etc.

Create another module for the gui, that calls the life module to calculate the next generation when needed. Making the program modular like this also leaves open the possiblity of making it suitable for distributed computing ie spread the calculation of the next state over a group of computers. Initially the gui module explicity called the update function in the life module, however I realised that this was not the way to go about things, it is better to run the life module as a server and send messages to it to update the state, and receive the new state as a message back.

GUI is simple: canvas within a window. Cells are represented by (square) rectangles of a given size, the canvas size divided by the cell size gives the range of cells that are plotted. When the state of the system is plotted the set containing the live cells is filtered so that only those cells within the plot window are plotted. Make the window listen for keyboard events eg n for next generation, ! to fill the canvas with a random distribution of cells, +/- to zoom in or out, cursor keys (or jkl;) to move the view window, g to go into a loop. Create the gui using the gs:create_tree function to create the window, packer frame and canvas in a single step. Each element of the tree has a unique name and the event handler listens needs to handle events from the appropriate elements eg mouse-clicks on the canvas to turn cells on and off.

%%% Code: life_gui.erl

-author('Matt McDonnell').

loop(LiveCells) ->
  DrawCells = sets:to_list(LiveCells),
  DrawnCells = lists:map(fun({X,Y})->
      gs:create(rectangle, can1, 
        [{coords,[{4*X,4*Y},{4*X+4,4*Y+4}]},{fill, black}]
      ) end,DrawCells),
  io:format("Number of live cells ~w~n",[sets:size(LiveCells)]),
    {gs,update_b,click,_,_} ->
      lists:foreach(fun (X) -> gs:destroy(X) end, DrawnCells),
      life ! {self(), {update,LiveCells}},
    {gs,quit_b,click,_,_} ->
      io:format("~w exiting normally: button~n",[self()]),
      life ! quit,
    {gs,window,keypress,_,[n,_KeyCode,_Shift,_Control | _ ]} ->
      lists:foreach(fun (X) -> gs:destroy(X) end, DrawnCells),
      life ! {self(), {update,LiveCells}},
    {gs,window,keypress,_,[c,_KeyCode,_Shift,_Control | _ ]} ->
      lists:foreach(fun (X) -> gs:destroy(X) end, DrawnCells),
    {gs,window,keypress,_,[q,_KeyCode,_Shift,_Control | _ ]} ->
      io:format("~w exiting normally: keyboard~n",[self()]),
      life ! quit,
    {gs,can1,buttonpress,_,[_,X,Y | _]} ->
      lists:foreach(fun (C) -> gs:destroy(C) end, DrawnCells),
      NewCell = {trunc(X/4),trunc(Y/4)},
      Exists = sets:is_element(NewCell, LiveCells),
        Exists == true -> 
        true -> 
    _Other -> 
      lists:foreach(fun (C) -> gs:destroy(C) end, DrawnCells),

update_loop() ->
    {life, Cells} ->
    _Other -> update_loop()

start() ->
  register(life, spawn(life, server, [])),

start_win() ->
  WH = [{width,384},{height,384}],
  R = [{window,window,[{map,true},{keypress,true} | WH],
          {packer_y,[{stretch,1},{stretch,4}]}], %Start frame children
          [{button,quit_b, [{label,{text, "Click to kill"}},
              {click, true},{pack_xy,{1,1}}] },
           {button,update_b, [{label,{text, "Update"}},
              {click, true},{pack_xy,{2,1}}] },
           {canvas,can1,  [{buttonpress,true}, {bg, red},
           {pack_xy,{{1,2},2}}] % End of frame children
         } ] } ] } ],

%%% End Code: life_gui.erl

Think about extending this problem to make use of the OTP design principles. The cellular automata itself is a gen_server that receives casts containing the current state of the system and then sends a reply containing the next state of the system. The gui is a "view" that displays the current state of the model and sends commands to the controller that distributes the computation of the next state to the set of gen_servers running the cellular automata update rule.

Some problems are introduced in using a set of processes to calculate the next generation of the system. The most obvious one is how to split up the domain of interest into a set of distinct domains. The simple method of splitting a square domain into four equally sized square regions is an acceptable starting point, however if most of the population is in one of the subdomains not much is accomplished. Ideally the partitioning should result in equal populations in each subdomain so the computational load is split equally amongst the processes. A possible solution is to introduce extra processes to monitor the distribution of population over the domain periodically and divide the domain into subdomains accordingly. Treemap distribution of domains? This problem does not arise in the more general case of simulations covering the entire domain (eg fluid dynamics simulations) where all cells in the simulation have to be considered (although in some problems the simulation could be sped up by ignoring cells with a value below a threshold, ie make the matrix that needs to be multiplied as sparse as possible).

The next problem is how to deal with the edges of the subdomains, as cells on the edge of a domain are influenced by cells in two separate subdomains. Consider the case where a cell is influenced by all cells that are distance D from the cell (in the game of life cellular automata D=1). The next state of all cells at a distance of >D from the edge just depend on the state of the cells in their subdomain. The cells lying within D of the edge form a separate domain that needs to be evolved separately, with the slight complication that the number of neighbours that each cell has within the subdomain must also be included in the count. For example when D=1 the number of neighbours that an edge cell has will be the number of neighbours in the first subdomain plus the number of neighbours in the second subdomain, take the number of neighbours on the same edge (to avoid double counting the edge neighbours). Alternatively, just form a domain consisting of the 2D+1 wide strips bounding the subdomains and only stitch the central strips onto the next generation.

Constrained Generating Procedures (CGPs)

John H. Holland's book "Emergence: from Chaos to Order" discusses the emergence of complexity and hierachical structures from systems consisting of a large number of interacting agents governed by simple rules. Examples are cellular automata, games such as draughts (checkers) or chess, statistical mechanics, etc.

In the case of a cellular automaton such as the game of life the individual agents are the cells on the board, each receiving inputs from its neighbouring cells. The natural way to think of modelling this in Erlang is to have each cell represented by a separate process, rather than the set of cells used in the implementation above. (Ignore for the moment the memory cost, roughly a couple of kB per process compared to a few bytes for the cell coordinates and state in the original implementation).

Consider the game of life again. Each cell is a process that changes state only after receiving messages from all the cells that are neighbouring it. For a live cell the neighbours are the surrounding cells. Dead cells only have the set of neighbouring live cells as their neighbours, to avoid having to simulate the entire infinite 2D grid. The state of the server controlling the simulation is specified by two sets of {coordinates, PID}, one set for the live cells and one for the dead cells. As each new live cell is added to the model (e.g. by clicking on the view window) new processes are spawned for the live cell and any surrounding dead cells not present in the dead cell set. A message sent to all neighbouring cells that are already in the dead cell set, adding the new live cell to their list of neighbours.

When the update command is received by the server it sends an update message to all of the cells in each set. All of the live cells send "I'm alive!" to their list of neighbours, all of the dead cells send "I'm not!" to theirs. Once a cell has received all of its messages it changes state and sends its new state to the server. Once the server has received messages from all of the cells it updates the view and waits for the next command. Additionally, a live cell that becomes dead and has no live neighbours should tell the server to remove it from both lists. A dead cell that becomes alive should send the coordinates of its neighbours and receive the PIDs from the server corresponding to any live neighbours that have come into existence.

Have the life_cgp module run as a server that receives update messages from the gui and replies with the new list of live cells. This is the same interface as the original life module, so from the point of view of the gui nothing has changed except the name of the server module spawned in the start function of the gui. The life_cgp server keeps dictionaries of {{X,Y}, PID} for the live and dead cells, and adds or removes cells from them in response to requests from the gui or messages from the cell processes.

Tiling the CGPs

An advantage of the CGP formulation of cellular automata is that small systems can be combined to form larger units that are still CGPs, giving an intrinsic sense of hierarchy. In the case of the life cellular automata the individual cell can be modelled as a CGP having two possible states and (up to) eight inputs from the states of its neighbouring cells. A larger system can be constructed by forming a 2x2 square of cells and treating this square as the CGP. The CGP now has 2^4=16 possible states and eight inputs from the neighbouring squares (also having 16 possible states each).