nc is qt

%% this is the best version! no trade off, version 8 only considers departure points
%% remove final empty queue assumption
%% speed performance
rng(42)

N = 100000;
in = exprnd(1, 1, N);
c = exprnd(3, 1, N);
q = zeros(1, N + 1);
q(1) = 10;

for i = 1:N
    q(i + 1) = max(0, q(i) + in(i) - c(i));
end

cumu_in = cumsum(in);
cumu_out = cumu_in - q(2:(N + 1));

delays = NaN(size(cumu_in));
in_len = length(cumu_in);

for i = 1:in_len
    
    depart_epoch = i; 
    arriv_epoch = i;    
    arriv_bits = cumu_in(i);  
    depart_bits = cumu_out(i);
    flag_break = 0;
    
    while (arriv_bits - depart_bits > 0)
        depart_epoch = depart_epoch + 1;
        if (depart_epoch > in_len)
            flag_break = 1;
            break;            
        end
        depart_bits = cumu_out(depart_epoch);
    end
    
    if (flag_break == 1)
        break;  
    end
    
    delays(i) = depart_epoch - arriv_epoch;        
end

mean(delays)


%% QT
t = c; % inter-arrival time
s = in; % service time

wq = zeros(1, N + 1);
w = zeros(1, N);
wq(1) = 10;

for i = 1:N
    wq(i + 1) = max(0, wq(i) + s(i) - t(i));
    w(i) = wq(i) + s(i);
end

arriv_epochs =[0 cumsum(t)];
depart_epochs = arriv_epochs(1:end - 1) + w;
pkt_num = length(arriv_epochs) - 1;
N_on_arriv = zeros(1, pkt_num);

% 求队长
for i = 1:pkt_num
    arriv_epoch = arriv_epochs(i);
    depart_idx = (i - 1);
    while (depart_idx > 0 && depart_epochs(depart_idx) - arriv_epoch > 0)
        N_on_arriv(i) = N_on_arriv(i) + 1;
        depart_idx = depart_idx - 1;            
    end
end

mean(N_on_arriv)