GKL_FitnessCalculation.m

% clear memory so that we can compare timespans
clear

tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CONSTANTS DECLARATION
% the size of our array of cells
latticeSize=149;
% neighborhood radius
radius = 3;
% rule size
ruleSize= radius*2 +1;
% number of time steps
max=320;
% number of initial configurations
ICsNo = 100;
%CONSTANTS DECLARATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%RULESET AND TRANSFORMATIONS
% the ruleset data object
ruleSet = CARuleset;
% patterns 
ruleSet.patterns = flipud(combn([0 1],ruleSize));
                
% declare transformation rules array - 128 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;
%RULESET AND TRANSFORMATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% generate initial configurations
ICs = generateBinaryInitialConfigurations(ICsNo, latticeSize);

% declare correct final state convergence counter
correctlySynced = 0;

% outer loop - tests the rule over all the initial configurations
icIterator=1;
while(icIterator <= ICsNo)
    % declare main operationl vector (a) and new vector buffer (newa)
    a=zeros(1,latticeSize);
    newa=zeros(1,latticeSize);

    % assign current initial configuration
    a = ICs(icIterator, :);

    % calculate if more 0 or 1 to have a check for final state
    sumInitialState = sum(a);
    disp(['iteration: ' num2str(icIterator)])
    if(sumInitialState > latticeSize/2)
        disp(['should converge to black - sum = ' num2str(sumInitialState)])
    else
        disp(['should converge to white - sum = ' num2str(sumInitialState)])
    end

    % 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;

    % this is the main loop (where the CA processing happens)
    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;

      % check if correctly synced or stable configuration and break if so
      sumCurrentState = sum(newa);
      if(sumCurrentState == 149 && sumInitialState > latticeSize/2)
          % increment correctly synced counter
          correctlySynced = correctlySynced +1;
          
          % some output
          disp(['correctly synced to all black (1) in ' num2str(g) ' time steps'])
          break;
      elseif (sumCurrentState == 0 && sumInitialState < latticeSize/2)
          % increment correctly synced counter
          correctlySynced = correctlySynced +1;
          
          % some output 
          disp(['correctly synced to all white (0) in ' num2str(g) ' time steps'])
          break;
      elseif(GRID(g,:) == GRID(g-1,:))
          % wrong stable config reached do nothing in this case - just break
          % some output first)
          disp(['wrong stable config reached in ' num2str(g) ' time steps'])
          break;
      end

    end

    disp(['iteration complete: ' num2str(icIterator)])
    icIterator = icIterator +1;

    %in case we want to play back
    %movie(M,0)
end

ratio = correctlySynced / ICsNo;

disp(['fitness: ' num2str(ratio)])
toc