%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plays Conway's game of life:
%
% 2D 2-state cellular automata using 8 cell neighbourhood and 
% the following update rule (states are 0 ("dead") and 1 ("live")):
%
%    S_{n}   N-live  S_{n+1}
%
%     0        3        1
%     1        2        1
%     1        3        1
%
%     OTHERWISE         0
%
% Here S_{n}, S_{n+1} are the states at the n-th and n+1-st time steps,
% and N-live is the number of live states (1's) in the neighbourhood at 
% step n.
% 
% This version uses 'dead' boundary conditions.  
% 
% Whole-array operations are used, which produces a significant 
% performance enhancement relative to 'life', which uses loops.
%
% Function 'latout' used to generate file output for subsequent 
% visualization with 'xflat2d_rgb'.
% 
% After running this script, use
%
%    xfpp2d_rgb < lifea.dat
%
% to visualize evolution.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Disable paging so that diagnostic (tracing) output appears as script 
% executes
more off;

% Script name (for diagnostic messages)
S = 'lifea';

% Initial data types
id_random = 0;
id_random_frac = 0.5;

id_glider = 1;

% Which initial data type?
id_type = id_random

% Dimensions of arena (lattice)
nx = 100
ny = 100

% Number of generations
nt = 500

% Define arena
c = zeros(nx,ny,nt);

% Define array used to count live cells in each site's neighbourhood
count = zeros(nx-2,ny-2);

% Define array used to temporarily store all lattice values except
% those on the boundaries
cn = zeros(nx-2,ny-2);

% Initialize arena
if id_type == id_random
   % Use random initial conditions in central part of arena
   % (TODO: controlled by id_random_frac), dead elsewhere
   wx = fix(id_random_frac * nx);
   wy = fix(id_random_frac * nx);
   c(2:nx-1,2:ny-1,1) = fix(2 * rand(nx-2,ny-2));
elseif id_type == id_glider
   % Set initial data to a single glider in lower-left hand corner
   % that will then propagate up and to the right
   cx = 5;
   cy = ny - 5;
   c(cx,  cy,  1) = 1;
   c(cx+1,cy,  1) = 1;
   c(cx+2,cy,  1) = 1;
   c(cx+2,cy+1,1) = 1;
   c(cx+1,cy+2,1) = 1;
else
   fprintf('%s: Invalid [id_type] = %d\n', S, id_type);
end

% Define upper and lower limits for updates in x and y directions.
% Sites on boundaries remain dead (0).
x0 = 2;
x1 = nx - 1;
y0 = 2;
y1 = ny - 1;

% Evolve for a total of nt generations
for n = 1 : nt - 1
   fprintf('%s: Step %d of %d\n', S, n, nt);
   % Compute the number of live cells in the neighbourhood ...
   count(:,:) = ...
      c(x0+1:x1+1, y0-1:y1-1, n) + ...
      c(x0+1:x1+1, y0  :y1  , n) + ...
      c(x0+1:x1+1, y0+1:y1+1, n) + ...
      c(x0  :x1  , y0-1:y1-1, n) + ...
      c(x0  :x1  , y0+1:y1+1, n) + ...
      c(x0-1:x1-1, y0-1:y1-1, n) + ...
      c(x0-1:x1-1, y0  :y1  , n) + ...
      c(x0-1:x1-1, y0+1:y1+1, n);
   cn = c(x0:x1, y0:y1, n);
   % ... and compute the new states using the rules of the game.
   cnp1 = (~cn .* (count == 3)) | (cn .* ((count == 2) | (count ==3)));
   c(x0:x1, y0:y1, n+1) = cnp1;
end

% Define colours for latout_rgb ...
rgb = [[0.3 0.3 0.3]; [0.0 1.0 0.0]];

% Output the entire evolution of the arena states ... Note that 
% the 4th argument to latout, which is optional, enables tracing
% messages in the function.  Its value (10 in this case), specifies
% how frequently messages will be output to the terminal in the 
% function's main loop ...
latout('lifea.dat', c, rgb, 10);
