%=================== CDA6530, Fall 2010 =============================
%========= Project 4 simulation code ================================
%===================================================================

% First, analytical results based on simulink code simpleModel.mdl:
eta = 2; N = 400; Omega = power(2, 14); I_0 = 1;
sim simpleModel;
simulTime = 400;
It_theory = zeros(1, simulTime); 
for i=1:simulTime, It_theory(i)=I_t(1,1,i); end

% Second, simulate worm propagation, get average results based on 100
% simulation runs
simulateN = 100; p = 0.3;
It = zeros(simulateN, simulTime); 
% each row of It saves the I(t) for one simulation runs
NodeStatus = zeros(1, N);
% we define a node in one of three status:
%   's':   susceptible
%   't':   being infected, not yet scan others (transit status)
%   'i'    infectious, can scan to infect others
NodeInfectiousTime = zeros(1,N); 
% record when a node begins to scan others

for runs=1:simulateN,   %  simulate worm propagation many runs

    %initialization at the beginning of each simulation run
    for i=1:N, NodeStatus(i) = 's'; end
    NodeInfectiousTime = zeros(1,N); 
    % randomly picks I_0 hosts to be infectious
    for i=1:I_0,  
        index = randi(N);  % a uniform distr. integer between 1 and N. 
        NodeStatus(index) = 'i';
        NodeInfectiousTime(index) = 1;
    end
    for t=1:simulTime, 
        % change 'transit' node to be 'infectious' if its infectious time
        % is now
        for i=1:N,
            if NodeStatus(i)=='t' && NodeInfectiousTime(i)==t,
                NodeStatus(i) = 'i';
            end
        end
                
        for i=1:N,    % check every node 
            if NodeStatus(i)=='i', % this node should scan and infect others
                for j=1:eta,
                    scanIP = randi(power(2,14)) - 1;
                    if scanIP < N*2 && mod(scanIP, 2)==0, % hit 
                        targetIndex = scanIP/2 + 1;
                        % update this target node 
                        if NodeStatus(targetIndex)=='s', % only update if it is 's'
                            NodeStatus(targetIndex) = 't';
                            % comment out one of the following section for 
                            % task #1 or task #2
                            
                            % Task #2
                            sum = 0;  % summation 
                            U = rand; x = 0;
                            while U > sum,  % must be true when k=0.
                                x = x+1;sum = sum + power((1-p), x-1)*p;
                            end
                                                        
                            % Task #1
                            x = 1;   
                                                        
                            NodeInfectiousTime(targetIndex) = t+x;
                        end
                    end
                end
            end
        end
        
        % check how many hosts are infectious at this time
        infectiousNum = 0;
        for i=1:N,
            if NodeStatus(i) ~='s',
                infectiousNum = infectiousNum + 1;
            end
        end
        % save the result
        It(runs, t) = infectiousNum;
    end  % end loop for a discrete time.
end   % end loop for simulation runs

% computer simulation average
It_simul = mean(It);

% draw output results
figure; hold on;
plot(It_theory, 'b'); plot(It_simul, 'r--');
legend('Theory', 'Simulation');
xlabel('Time: t'); ylabel('I(t)');
    
        
        
                
                        
                        
                

