clear
a = [0.03125 0.02 0.015 0.015 0.05 0.025]; %these are the probabilities of arrival of each of the processes
pa = [a(1) a(2) a(3) a(4) a(5) a(6)]; %these are the probabilities that there are 0 arrivals after returning from vacation
N = 6;
start_time = 0;
task_start_time = [ 0, 0, 0, 0, 0, 0 ];
total_task_time = [ 500, 1000, 2000, 500, 2000, 1000 ];
prev_analyzed_ratio = [ 0, 0, 0, 0, 0, 0 ];

pe = 0.25;
K = 5;

B = [ 1 zeros(1,K-1)]; %this is the start of the Erlang K memory service
S0 = [ zeros(K-1,1); pe ]; %this is the absorbtion vector of the memory service
S = zeros(K,K);

%this loop builds the Erlang K service matrix
for i=1:K
    S(i,i) = 1-pe;
    if (i < K)        
        S(i,i+1) = pe;
    end
end

num_loops = 0; %this is used to keep track of the number of loops required for convergence
converge = 0; % a flag to determine if we have converged
while (converge == 0)
    pa_old = pa;
    for i = 1:N  %do one loop for every process
        %v_st is the starting vector of the vacation process
        v_st = [ (1-pa(mod(i,N)+1))      pa(mod(i,N)+1)*B     0                                   zeros(1,K)                              0                                   zeros(1,K)                                0                                   zeros(1,K)                                  0                                   zeros(1,K)];
        
        %V is the matrix that represents the vacation before it is done
        V = [ 0                     a(mod(i,N)+1)*B         (1-pa(mod(i+1,N)+1))*(1-a(mod(i,N)+1))   pa(mod(i+1,N)+1)*(1-a(mod(i,N)+1))*B 0                                   zeros(1,K)                                0                                   zeros(1,K)                                  0                                   zeros(1,K);       
              (1-a(mod(i,N)+1))*S0  a(mod(i,N)+1)*(S0*B)+S  zeros(K,1)                          zeros(K,K)                              zeros(K,1)                          zeros(K,K)                                  zeros(K,1)                          zeros(K,K)                                  zeros(K,1)                          zeros(K,K);       
              0                     zeros(1,K)              0                                   a(mod(i+1,N)+1)*B                       (1-pa(mod(i+2,N)+1))*(1-a(mod(i+1,N)+1)) pa(mod(i+2,N)+1)*(1-a(mod(i+1,N)+1))*B 0                                   zeros(1,K)                                  0                                   zeros(1,K);      
              zeros(K,1)            zeros(K,K)              (1-a(mod(i+1,N)+1))*S0              a(mod(i+1,N)+1)*(S0*B)+S                zeros(K,1)                          zeros(K,K)                                  zeros(K,1)                          zeros(K,K)                                  zeros(K,1)                          zeros(K,K);       
              0                     zeros(1,K)              0                                   zeros(1,K)                              0                                   a(mod(i+2,N)+1)*B                       (1-pa(mod(i+3,N)+1))*(1-a(mod(i+2,N)+1)) pa(mod(i+3,N)+1)*(1-a(mod(i+2,N)+1))*B     0                                   zeros(1,K);    
              zeros(K,1)            zeros(K,K)              zeros(K,1)                          zeros(K,K)                              (1-a(mod(i+2,N)+1))*S0              a(mod(i+2,N)+1)*(S0*B)+S                zeros(K,1)                          zeros(K,K)                                      zeros(K,1)                          zeros(K,K);              
              0                     zeros(1,K)              0                                   zeros(1,K)                              0                                   zeros(1,K)                              0                                   a(mod(i+3,N)+1)*B                               (1-pa(mod(i+4,N)+1))*(1-a(mod(i+3,N)+1)) pa(mod(i+4,N)+1)*(1-a(mod(i+3,N)+1))*B;    
              zeros(K,1)            zeros(K,K)              zeros(K,1)                          zeros(K,K)                              zeros(K,1)                          zeros(K,K)                              (1-a(mod(i+3,N)+1))*S0              a(mod(i+3,N)+1)*(S0*B)+S                        zeros(K,1)                          zeros(K,K);
              0                     zeros(1,K)              0                                   zeros(1,K)                              0                                   zeros(1,K)                              0                                   zeros(1,K)                                      0                                   a(mod(i+4,N)+1)*B ;    
              zeros(K,1)            zeros(K,K)              zeros(K,1)                          zeros(K,K)                              zeros(K,1)                          zeros(K,K)                              zeros(K,1)                          zeros(K,K)                                      (1-a(mod(i+4,N)+1))*S0              a(mod(i+4,N)+1)*(S0*B)+S ];
        
        % V0 is the vector representing the end of vacation
        % (i.e.absorbtion)
        V0 = [ 0;         
               zeros(K,1);
               0;         
               zeros(K,1);
               0;         
               zeros(K,1);
               0;         
               zeros(K,1);
               1-a(mod(i+4,N)+1);  
               zeros(K,1)];
        
        %V_prime combines v_st V0 and V in a matrix in order to determine
        %the mean time to absorbtion of the vacation process
        V_prime = [ 0 v_st;
                    V0 V ];
                
        %the probability of an arrival occuring during the vacation period
        %is now calculated.  If this changes it will change the vacation
        %period of the other queues, so the process is iterated until these
        %values do not change.
        pa(i) = calc_phi(V_prime,a(i));
    end

    %now check for convergence
    converge_err = pa_old - pa;
    converge = 1;
    for g = 1:N
        if (abs(converge_err(g)) > 10^-12)
            converge = 0;
        end
    end
    num_loops = num_loops + 1;
end

pa
num_loops

%now that the values for pa have been determined, we can now evaluate each
%of the queues to get the distribution, and from that determine the amount
%of time that the process spends waiting.
for i = 1:N
   %v_st is the starting vector of the vacation process
        v_st = [ (1-pa(mod(i,N)+1))      pa(mod(i,N)+1)*B     0                                   zeros(1,K)                              0                                   zeros(1,K)                                0                                   zeros(1,K)                                  0                                   zeros(1,K)];
        
        %V is the matrix that represents the vacation before it is done
        V = [ 0                     a(mod(i,N)+1)*B         (1-pa(mod(i+1,N)+1))*(1-a(mod(i,N)+1))   pa(mod(i+1,N)+1)*(1-a(mod(i,N)+1))*B 0                                   zeros(1,K)                                0                                   zeros(1,K)                                  0                                   zeros(1,K);       
              (1-a(mod(i,N)+1))*S0  a(mod(i,N)+1)*(S0*B)+S  zeros(K,1)                          zeros(K,K)                              zeros(K,1)                          zeros(K,K)                                  zeros(K,1)                          zeros(K,K)                                  zeros(K,1)                          zeros(K,K);       
              0                     zeros(1,K)              0                                   a(mod(i+1,N)+1)*B                       (1-pa(mod(i+2,N)+1))*(1-a(mod(i+1,N)+1)) pa(mod(i+2,N)+1)*(1-a(mod(i+1,N)+1))*B 0                                   zeros(1,K)                                  0                                   zeros(1,K);      
              zeros(K,1)            zeros(K,K)              (1-a(mod(i+1,N)+1))*S0              a(mod(i+1,N)+1)*(S0*B)+S                zeros(K,1)                          zeros(K,K)                                  zeros(K,1)                          zeros(K,K)                                  zeros(K,1)                          zeros(K,K);       
              0                     zeros(1,K)              0                                   zeros(1,K)                              0                                   a(mod(i+2,N)+1)*B                       (1-pa(mod(i+3,N)+1))*(1-a(mod(i+2,N)+1)) pa(mod(i+3,N)+1)*(1-a(mod(i+2,N)+1))*B     0                                   zeros(1,K);    
              zeros(K,1)            zeros(K,K)              zeros(K,1)                          zeros(K,K)                              (1-a(mod(i+2,N)+1))*S0              a(mod(i+2,N)+1)*(S0*B)+S                zeros(K,1)                          zeros(K,K)                                      zeros(K,1)                          zeros(K,K);              
              0                     zeros(1,K)              0                                   zeros(1,K)                              0                                   zeros(1,K)                              0                                   a(mod(i+3,N)+1)*B                               (1-pa(mod(i+4,N)+1))*(1-a(mod(i+3,N)+1)) pa(mod(i+4,N)+1)*(1-a(mod(i+3,N)+1))*B;    
              zeros(K,1)            zeros(K,K)              zeros(K,1)                          zeros(K,K)                              zeros(K,1)                          zeros(K,K)                              (1-a(mod(i+3,N)+1))*S0              a(mod(i+3,N)+1)*(S0*B)+S                        zeros(K,1)                          zeros(K,K);
              0                     zeros(1,K)              0                                   zeros(1,K)                              0                                   zeros(1,K)                              0                                   zeros(1,K)                                      0                                   a(mod(i+4,N)+1)*B ;    
              zeros(K,1)            zeros(K,K)              zeros(K,1)                          zeros(K,K)                              zeros(K,1)                          zeros(K,K)                              zeros(K,1)                          zeros(K,K)                                      (1-a(mod(i+4,N)+1))*S0              a(mod(i+4,N)+1)*(S0*B)+S ];
        
        % V0 is the vector representing the end of vacation
        % (i.e.absorbtion)
        V0 = [ 0;         
               zeros(K,1);
               0;         
               zeros(K,1);
               0;         
               zeros(K,1);
               0;         
               zeros(K,1);
               1-a(mod(i+4,N)+1);  
               zeros(K,1)];

    [ orderV, null ] = size(V);
    
    V_prime = [ 0 v_st;
                V0 V ];
              
    %this is the probability transition matrix for the given process i
    P = [ 0               a(i)*B          (1-a(i))*v_st           zeros(1,orderV);
          (1-a(i))*S0     a(i)*(S0*B)+S   zeros(K,orderV)         zeros(K,orderV);
          (1-a(i))*V0     a(i)*(V0*B)     (1-a(i))*V              a(i)*V;
          zeros(orderV,1) V0*B            zeros(orderV,orderV)    V ];
    
    %now we find the stationary probability vector
    [ orderP, null ] = size(P);
    pi = [ 1 zeros(1,orderP-1) ];

    error = zeros(1,orderP);
    loop = 1; %this will be set to 1 if the error of any term is greater than 10^-12
    while (loop == 1)
        pi_prime = pi*P;
        error = pi_prime-pi;
        pi = pi_prime;  
        loop = 0;
        for r=1:orderP
            if (abs(error(r)) > 10^-12)
                loop = 1;
            end               
        end
    end
  
    figure
    bar(pi) % display the stationary probability
    probability_of_wait = 0;
    waiting_while_others_served = 0;
    probability_of_no_mem_req = pi(1);
    
    for t = 2:K+1        
        probability_of_wait = probability_of_wait + pi(t);        
    end
    
    for t = K+2:K+orderV+1
        probability_of_no_mem_req = probability_of_no_mem_req + pi(t);
    end
    
    for t = 2+K+orderV:orderP
        probability_of_wait = probability_of_wait + pi(t);
        waiting_while_others_served = waiting_while_others_served + pi(t);
    end
    i
    probability_of_wait % this is the probability that the process will be waiting
    waiting_while_others_served % this is the probability that the process will be waiting while another process is being served
    probability_of_no_mem_req % this is the probability that the process is not waiting because there is no memory request
    pw(i) = waiting_while_others_served;
end

for i=1:N
    %this calculation is done in order to calculate the new partition time
    task_time_remaining(i) = (total_task_time(i) - total_task_time(i)*(prev_analyzed_ratio(i)))/(1-pw(i));    
end

t_partition_new = min(task_time_remaining) % this is the new partition time
finish_time = start_time + t_partition_new

for i=1:N
    ratio_of_task_analyzed(i) = t_partition_new*(1-pw(i)) / total_task_time(i);
    task_time_remaining(i) = total_task_time(i) - total_task_time(i)*(ratio_of_task_analyzed(i)+prev_analyzed_ratio(i));
    expected_end_time(i) = finish_time - task_start_time(i) + task_time_remaining(i);
end

ratio_of_task_analyzed
total_ratio_analyzed = ratio_of_task_analyzed+prev_analyzed_ratio
task_time_remaining
expected_end_time
