%
% gene3.m  Steve Cox  
%
% run and display the evolution of a gene net from a 
% number of initial states
%
% usage:	gene3(wire, rule, nr, ipr)
%
% where		wire is an n-by-3 matrix
%		rule is an n-by-1 vector
%		nr = number of runs
%		ipr = iterations per run
%
% subfunctions:  rtrans
%

function gene3(wire, rule, nr, ipr)

n = length(rule);

for k=1:n,
    ruletab(k,:) = rtrans(rule(k));	% build the rule table
end

disp(' ')
disp('The Rule Table is')
disp(' ')
disp(ruletab)
disp(' ')

for uni=1:nr,

    s = round(rand(1,n));		% initial state

    for j=1:ipr,

        hot = find(s>0.1);		% genes that are on

        plot(hot+10*(uni-1),-j*ones(size(hot)),'s','markersize',2)
        hold on

        for k=1:n,

            pat = s(wire(k,:));		% input pattern at gene k

            cc = 1 + pat*[4 2 1]';	% associated column in ruletab

            ns(k) = ruletab(k,cc);	% next state of gene k
    
        end

        s = ns;		% copy state

    end  % ipr loop

end  % nr loop

hold off
axis off

return

%
% convert the rule number, rnum, to its binary row, brule
%

function brule = rtrans(rnum)
brule = zeros(1,8);
h = floor(log2(rnum));
brule(h+1) = 1;
rnum = rnum - 2^h;

while rnum > 0

      h = floor(log2(rnum));
      brule(h+1) = 1;
      rnum = rnum - 2^h; 

end