%=================== CDA6530, Fall 2012 =============================
%========= discrete event simulation of the three callers =========== 
%========problem in homework 2 ======================================
%====================================================================

lambda = 0.125;  mu = 0.2;
% theoretical results
Q = [-3*lambda  3*lambda      0       0       0       0;
     mu -mu-2*lambda  2*lambda        0       0       0;
     0 2*mu -2*mu-lambda lambda 0 0;
     0 0 0 -2*mu 2*mu 0;
     0 0 2*lambda lambda -mu-3*lambda mu;
     0 2*lambda 0 0 2*lambda -4*lambda    ];
Pi_theory = zeros(1, 6);
Q_m = [Q(:, 1:5)  ones(6,1)];
B = [0 0 0 0 0 1];
Pi_theory = B * inv(Q_m);

%-------------------------------------------
%--- discrete event simulation -------------
stateN = 6; N = 1000;  % N is the number of transitions.
initState = 1;         % initial system's state
Tran = zeros(N,1); State = zeros(N,1);
currentState = initState; currentTime = 0; 
for i=1:N,    % simulate N transitions
       % first, simulation currentState during time (next event time)
outRate =  - Q(currentState, currentState);
Tran(i) = currentTime - log(rand)/outRate; % exp. distr. with rate of outRate
% next, determine which state transits to?
U = rand; 
vector = Q(currentState,:); vector(currentState) = 0;
for j=1:stateN,
	if U <= sum(vector(1:j))/sum(vector),
		nextState = j; break;
     end
end
State(i) = nextState;
currentState = nextState; currentTime = Tran(i);
end

Time = zeros(6,1); Time(initState) = Tran(1); 
for k=1:6,
    for i=1:N-1,
        if State(i) == k,
            Time(k) = Time(k) + Tran(i+1) - Tran(i);
        end
    end
end

Pi_simul = zeros(6,1);
Pi_simul = Time./Tran(N);

figure; hold on; 
plot(Pi_theory, '-ob'); plot(Pi_simul, '--^k');
legend('Theoretical', 'Simulation');

