GKL_CASimulation.m

% the size of our array of cells
latticeSize=149;
% neighborhood radius
radius = 3;
% rule size
ruleSize= radius*2 +1;
% number of time steps
max=200;

a=zeros(1,latticeSize);
newa=zeros(1,latticeSize);

% the ruleset data object
ruleSet = CARuleset;

% patterns 
ruleSet.patterns = flipud(combn([0 1],ruleSize));
                
% declare transformation rules array - 148 bits (2^7) initialized to zero
gkl = zeros(1,2^ruleSize);
% GKL rule: 
% If  ci(t) = 0, then ci(t + 1) = majority [ci(t), ci-1(t) ci-3(t)];
% If  ci(t) = 1, then ci(t + 1) = majority [ci(t), ci+1(t) ci+3(t)];
for j=1:2^ruleSize,
    pattern = ruleSet.patterns(j,:);
    if pattern(4) == 0,
        if (pattern(4-1) && pattern(4-3)),
            gkl(j) = 1;
        end
    else
        if (pattern(4+1) || pattern(4+3)),
            gkl(j) = 1;
        end
    end
end

% assign resulting transformations
ruleSet.transformations = gkl;

%randomize initial configuration
a = randint(1, latticeSize, [0 1]);
%a(50) = 1;

%assign initial configuration to the grid
GRID(1,:)=a;

% initialize the grid except 1st row (initial configuration)
for i=2:max,
    GRID(i,:)=zeros(1,latticeSize);
end

%spit out first frame
spy(GRID, 'k')
M(1) = getframe;

%brace yourself - this is the main loop
g=1;
while (g<max), 
    
  %run the chosen rule foreach cell for a given time step g
  %boundary cells are being ignored in this example
  for i=1:latticeSize,  
    % retrieve pattern in the local 'hood 
    localPattern = circularSubarray(a, i-radius, i+radius);
    % find the transformation rule for the given local pattern
    for c=1:2^ruleSize,
        if(isequal(ruleSet.patterns(c, :), localPattern))
            newa(i) = ruleSet.transformations(c);
            break
        end
    end
  end
  
  % assign new states
  g=g+1;
  a=newa;
  GRID(g,:)=a;
  
  %push a snapshot to a frame
  spy(GRID, 'k')
  M(g) = getframe;
end

%last frame - final state of the matrix
spy(GRID, 'k')
M(max) = getframe;

%playback
movie(M,0)