%Student: Santiago Pagola
%LiU-ID & Personal Number: 930302-T092, sanpa993
 
 
%clear all variables from workspace, close all figures and clear everything
%in command line
 
clear all
close all
clc
 
%IMPORTANT NOTE: AN EXTERNAL FUNCTION WAS CREATED IN ORDER TO WRAP THE SUMS
%OF THE VECTORS: a = sumk(vector,k). It will be described in the report,
%but basically sums the first k elements of an array called 'vector'.
 
m = 100; %100 nodes on no-buffering assumption
t = 1000; %Time: 1000 slots of time
backlog_array = zeros(size(1:t)); %Array of the values of the system backlog
backlog_estimate = zeros(size(1:t)); %Array of the estimated values of the system backlog
slots = 1:t; %Slots array
n = 0; %Initially, all nodes are unbacklogged
%lambda = 1/exp(1); %Overal arrival rate will be set to 1/exp(1) (max. stable throughput)
n_estimated = 0; %Estimated backlog initialization. This will be updated during runtime
q_r = 0; %Value of q_r is now stabilized by this algorithm
packets_leaving = 1:t; %Array which will count #packets leaving the system
packets_arriving = 1:t; %Array which will count #packets entering the system
state_probs = zeros(size(1,m)); %State probabilitites array
W = zeros(size(1,7));%Array definition for the P.B. delay
D = zeros(size(1,7));%Array definition for LIttle Theorem delay
v = 1; %Auxiliary variable
%main loop
for lambda2 = 0.05:0.05:0.35
    qa = 1-exp(-lambda2/m); %qa is the probability that each unbacklogged node has to transmit in the next slot
    for i = 1:t
    a = rand(1); %realization of Pr_x
    b = rand(1); 
    c = rand(1);
 
    %Estimation of q_r starting from the estimated backlog
    if n_estimated >= 0 && n_estimated < 1
        q_r = 1;
    else
        q_r = 1/n_estimated;
    end
 
    %Now we are going to create an array with the Poisson pmf's taking into
    %account the number of transmitting nodes
    Pr = zeros(size(1:11));
    for j = 0:10
        Pr(j+1) = ((lambda2)^j/factorial(j))*exp(-lambda2);
    end
 
    if 0 <= a && a <= Pr(1)
        n;
    elseif sumk(Pr,1) < a && a <= sumk(Pr,2)
        n = n+1;
    elseif sumk(Pr,2) < a && a <= sumk(Pr,3)
        n = n+2;
    elseif sumk(Pr,3) < a && a <= sumk(Pr,4)
        n = n+3;
    elseif sumk(Pr,4) < a && a <= sumk(Pr,5)
        n = n+4;
    elseif sumk(Pr,5) < a && a <= sumk(Pr,6)
        n = n+5;
    elseif sumk(Pr,6) < a && a <= sumk(Pr,7)
        n = n+6;
    elseif sumk(Pr,7) < a && a <= sumk(Pr,8)
        n = n+7;
    elseif sumk(Pr,8) < a && a <= sumk(Pr,9)
        n = n+8;
    elseif sumk(Pr,9) < a && a <= sumk(Pr,10)
        n = n+9;
    elseif sumk(Pr,10) < a && a <= 1
        n = n+10;
    end
 
    %Now we are going to create an array with the Qr's (%Probabilities that up to 10 backlogged nodes
    %retransmit (up to 10 new arrivals))
    Qr = zeros(size(1:11));
    for j = 0:10
        Qr(j+1) = binopdf(j,n,q_r);
    end
 
    if 0 <= n && n < m %Check if the backlog is less than m = 100
 
        if 0 <= b && b <= sumk(Qr,1)   %CASE 1: IDLE SLOT, NO NEW ARRIVALS, FEEDBACK 0            
                n_estimated=max(lambda2, n_estimated + lambda2 - 1);
                packets_arriving(i) = 0;
                packets_leaving(i) = 0;
 
        elseif sumk(Qr,1) < b && b <= sumk(Qr,2) %CASE 2: NO RETRASMISSION + 1 UNBACKLOGGED NODE TRANSMISION -> SUCCESS! , FEEDBACK 1               
                n = n-1;
                n_estimated = max(lambda2, n_estimated + lambda2 - 1);               
                packets_arriving(i) = 1;
                packets_leaving(i) = 1;
 
        elseif sumk(Qr,2) < b && c <= 1    % CASE 3: COLLISION: RETRANSMISSION OF MORE THAN ...
 
            if sumk(Qr,2) < b && b <= sumk(Qr,3) %... 2 PACKETS
                x = 2;
            elseif sumk(Qr,3) < b && b <= sumk(Qr,4) %...3 PACKETS
                x = 3;
            elseif sumk(Qr,4) < b && b <= sumk(Qr,5) %...AND SO ON
                x = 4;
            elseif sumk(Qr,5) < b && b <= sumk(Qr,6)
                x = 5;
            elseif sumk(Qr,6) < b && b <= sumk(Qr,7)
                x = 6;
            elseif sumk(Qr,7) < b && b <= sumk(Qr,8)
                x = 7;
            elseif sumk(Qr,8) < b && b <= sumk(Qr,9)
                x = 8;
            elseif sumk(Qr,9) < b && b <= sumk(Qr,10)
                x = 9;
            elseif sumk(Qr,10) < b && b <= 1
                x = 10;
            end
 
            %estimation of n for FEEDBACK e
            n_estimated = n_estimated + lambda2 + (exp(1)-2)^-1;
            packets_arriving(i) = x;
            packets_leaving(i) = 0;
 
        end
    end
    backlog_array(i) = n;
    backlog_estimate(i) = n_estimated; %Vector that saves the estimated backlog after each slot
 
    end
    N = mean(backlog_array);
    nelements = hist(backlog_array,m); %In each bin, it counts how many times element i was seen ("numerical plot of the hist function")
    D(v) = N/lambda2;
    W(v) = ((exp(1)-0.5)/(1-lambda2*exp(1)))-(((exp(1)-1)*((exp(1)^lambda2)-1))/(lambda2*(1-((exp(1)-1)*((exp(1)^lambda2)-1)))));
    v = v + 1;
end
figure(1) %Setting up the plotting environment for the backlog of the system
xlabel('Slot number, n')
ylabel('Backlogged packets')
plot(slots, backlog_array)
hold on
plot(slots,backlog_estimate,'g')
title('Real backlog (blue) vs. Pseudo Bayesian Estimated backlog (green)')
 
figure(2) %Setting up the plotting environment for the packets of the system
plot(slots,packets_arriving)
hold on
plot(slots,packets_leaving,'k')
grid on
xlabel('Slot number')
ylabel('Number of packets')
title('Number of packets entering (blue) or leaving (black) the system vs. Slot number')
 
%Question 8
max = max(backlog_array);
 
lambda2array = [0.05 0.1 0.15 0.2 0.25 0.3 0.35];
 
figure(3)
plot(lambda2array,D)
xlabel('Variation of lambda')
ylabel('Average delay')
title('Average delay according to Little Theorem and P-B Stabilization')
hold on
plot(lambda2array,W,'r')
legend('Little Theorem','P-B Stabilization')