Autômato Celular Elementar

% Simulates any elementary automaton rule.
%Rev1 speed is ~20k cells/s, rev2 ~200k cells/s, rev3 ~1000k cells/s, rev4 ~2M cells/s on AMD Phenom II X4 955 at 3.2GHz.
%There are multiple versions here, please select one (I recommend rev5).
% It runs a cellular automata with the specified rule (int 1~255), for a specified amount of steps (int > 1) and a initial pattern (line vector).

---------------------------------------------------------------------------------

%Rev 1.0 (reconstruction, for comparison):

function a = cellsim(rule,steps,initial)

cells=[0,0,initial(1,:),0,0]; % Bidimensional array for cells with zeros border.
cells(2,:)=zeros();
rule=de2bi(rule,8,'right-msb'); % Transform decimal rule into binary matrix.
state=0;

for cstep = 2:steps+1
    for k = 2:size(cells,2)-1
        state=bi2de(cells(cstep-1,k-1:k+1),'left-msb')+1;
        % Cell/state check: (checks state of neighbors cells and applies rule)
        if rule(1,state) == 1
            cells(cstep,k)=1;
        else
            cells(cstep,k)=0;
        end
    end
    % Boundaries checks: (resize matrix so it won't overflow)
    if cells(cstep,2) == 1
        cells(:,2:end+1)=cells(:,:); % Move matrix right and fill cells with zeros.
    end
    if cells(cstep,end-1) == 1
        cells(:,end+1)=zeros(); % Add right column filled with zeros.
    end
end

cells=cells(:,3:end-2);

a=cells;

---------------------------------------------------------------------------------

%Rev 2.0: Changed function "bi2de" to a matrix multiplication, made it ~10x faster (from 19k points/s to >200k points/s).

function a = cellsim(rule,steps,initial)

cells=[0,0,initial(1,:),0,0]; % Bidimensional array for cells with zeros border.
cells(2,:)=zeros();
rule=de2bi(rule,8,'left-msb'); % Transform decimal rule into binary matrix.

for cstep = 2:steps+1
    for k = 2:size(cells,2)-1
        % Cell/state check: (checks state of neighbors cells and applies rule)
        if rule(1,8-cells(cstep-1,k-1:k+1)*[4;2;1]) == 1
            cells(cstep,k)=1;
        else
            cells(cstep,k)=0;
        end
    end
    % Boundaries checks: (resize matrix so it won't overflow)
    if cells(cstep,2) == 1
        cells(:,2:end+1)=cells(:,:); % Move matrix right and fill cells with zeros.
    end
    if cells(cstep,end-1) == 1
        cells(:,end+1)=zeros(); % Add right column filled with zeros.
    end
end

cells=cells(:,3:end-2);

a=cells;

---------------------------------------------------------------------------------

%Rev 3.0: Changed loop for determining the value of past cells to a matrix manipulation, on a single turn. About 4x faster. (~1M points/s).

function a = cellsim(rule,steps,initial)

cells=[0,initial(1,:),0]; % Bidimensional array for cells with zeros border.
cells(2,:)=zeros();
rule=de2bi(rule,8,'right-msb'); % Transform decimal rule into binary matrix.

for cstep = 2:steps+1
    past=[cells(cstep-1,:),0,0;0,cells(cstep-1,:),0;0,0,cells(cstep-1,:)]'*[4;2;1]+1;
    for k = 2:size(past,1)-1
        % Cell/state check: (checks state of neighbors cells and applies rule)
        if rule(1,past(k,1)) == 1
            cells(cstep,k-1)=1;
        else
            cells(cstep,k-1)=0;
        end
    end
    % Boundaries checks: (resize matrix so it won't overflow)
    if cells(cstep,1) == 1
        cells(:,2:end+1)=cells(:,:); % Move matrix right and fill cells with zeros.
    end
    if cells(cstep,end) == 1
        cells(:,end+1)=zeros(); % Add right column filled with zeros.
    end
end

cells=cells(:,3:end-2);

a=cells;

---------------------------------------------------------------------------------

%Rev 4.0: Matrix preallocation on both axes.

function a = cellsim(rule,steps,initial)

rule=de2bi(rule,8,'right-msb'); % Transform decimal rule into binary matrix.

cells(steps,size(initial,2)+2*steps)=zeros();
cells(1,steps+1:steps+size(initial,2))=initial(1,:); % Preallocating matrix
start=steps; % Virtual limits
final=size(initial,2)+steps;

for cstep = 2:steps+1
    past=[cells(cstep-1,2:end),0;cells(cstep-1,:);0,cells(cstep-1,1:end-1)]'*[4;2;1]+1;
    for k = start:final
        % Cell/state check: (checks state of neighbors cells and applies rule)
        if rule(1,past(k,1)) == 1
            cells(cstep,k)=1;
        else
            cells(cstep,k)=0;
        end
    end
    if cells(cstep,start) == 1
        start=start-1;
    end
    if cells(cstep,final) == 1
        final=final+1;
    end
end

a=cells(:,start+1:final-1);

---------------------------------------------------------------------------------

%Rev 4.1: Order selection (elementary or second order).

function a = cellsim(rule,order,steps,initial)

rule=de2bi(rule,8,'right-msb'); % Transform decimal rule into binary matrix.

cells(steps,size(initial,2)+2*steps)=zeros();
cells(1,steps+1:steps+size(initial,2))=initial(1,:); % Preallocating matrix
start=steps; % Virtual limits.
final=size(initial,2)+steps;

for cstep = 2:steps
    past=[cells(cstep-1,2:end),0;cells(cstep-1,:);0,cells(cstep-1,1:end-1)]'*[4;2;1]+1;
    for k = start:final+1
        % Cell/state check: (checks state of neighbors cells and applies rule)
        if xor(rule(1,past(k,1)) == 1, (order-1)*cells(max(cstep-2,1),k))
            cells(cstep,k)=1;
        else
            cells(cstep,k)=0;
        end
    end
    if cells(cstep,start) == 1 % Modifying "virtual" limits if necessary.
        start=start-1;
    end
    if cells(cstep,final+1) == 1
        final=final+1;
    end
end

a=cells(:,start+1:final-1); % Cropping matrix.

---------------------------------------------------------------------------------

%Rev 5: Parallel method, very simple & fast. Completely "manual".

function a = cellsim(rule,steps,initial)

rule=de2bi(rule,8,'right-msb'); % Transform decimal rule into binary matrix.

cells=zeros(steps,size(initial,2));
cells(1,:)=initial;

for cstep = 2:steps+1
    past=([cells(cstep-1,2:end),0;cells(cstep-1,1:end);0,cells(cstep-1,1:end-1)]'*[4;2;1]+1)';
    for k = 1:8
        past(past==k)=rule(1,k);
        cells(cstep,:)=past;
    end
end

a=cells;

----------------------------------

%Rev 6:

function a = cellsim(rule,rule2,order,steps,initial)

rule=de2bi(rule,8,'right-msb'); % Transform decimal rule into binary matrix.
rule(2,1:8)=de2bi(rule2,8,'right-msb');

cells=zeros(steps,size(initial,2));
cells(1,:)=initial;

for cstep = 2:steps+1
    past=([cells(cstep-1,2:end),0;cells(cstep-1,1:end);0,cells(cstep-1,1:end-1)]'*[4;2;1]+1)';
    for k = 1:8
        past(past==k)=rule(1,k);
    end
    if order==1 || cstep==2
        cells(cstep,:)=past;
    elseif order==2
        cells(cstep,:)=bitxor(past,cells(cstep-2,:));
    elseif order==3
        tpast=([cells(cstep-2,2:end),0;cells(cstep-2,1:end);0,cells(cstep-2,1:end-1)]'*[4;2;1]+1)';
        for k = 1:8
            tpast(tpast==k)=rule(2,k);
        end
        cells(cstep,:)=bitxor(past,tpast);
    end
end

a=cells;