% Solution to project M6.2, chapter 6.4 of Dekking. Concocted by David Dobor (TA) 
% for CIS-2033, Spring 2013. Uncomment lines 16 and 18 to replicate 
% the figures in the book almost exactly. 
function single_server_queue

%------------------------------generate data------------------------------%
n = 50;                    % number of customers to be processed 
R = 2 + (5-2) * rand(1,n); % water requirements for the customers
T = -log(rand(1,n))/0.5;   % exp(.5) inter-arrival times of customers
v1 = 2; v2 = 3;            % two pump capacities; can vary these (e.g. 1, 5)
t = 100;                   % time horizon (in minutes). (side note: E[sum(T)], since each of the 50 T has E[T} = 2)

% uncomment the following two lines for R and T to replicate the example
% from the book
% T - inter-arrival times from the book:
% T = [0.24 1.97 1.73 2.82 1.01 1.09];
% R - water requirements (bucket sizes) from the book:
% R = [4.39 4.00 2.33 4.03 4.17 4.24];

% the service times for for the two pump capacities:
S_v1 = R/v1; 
S_v2 = R/v2;
%-------------------------------------------------------------------------%

%---compute and plot average waiting times for the two pump capacities----%
W1 = computeWaitingTimes(T, S_v1);
avgWT1 = computeAverageWTs(W1);

W2 = computeWaitingTimes(T, S_v2);
avgWT2 = computeAverageWTs(W2);

% plot the results
figure(1)
plot((1:length(avgWT1)), avgWT1, (1:length(avgWT2)), avgWT2);
title('\it{Average Waiting Times for Two Pump Capacities}')
ylabel('\it{Average Waiting Time, W_n}'); xlabel('\it{Customer}');
hleg = legend(sprintf('v_1 = %d', v1), sprintf('v_2 = %d ', v2));
set(hleg,'Location','NorthWest','FontAngle','italic','TextColor',[.3 .2 .1]);
%-------------------------------------------------------------------------%

%-------------------compute and plot work-in-system-----------------------%
figure(2)
workInSystem (T, S_v1)
set(gca, 'YLim', [-.2 Inf], 'YLim', [-.2 Inf])
set(gca,'Title',text('String',sprintf('Work-In-System for v_1 = %d ', v1),'FontAngle','italic'))
set(gca,'XLabel',text('String','Time Horizon','FontAngle','italic'))
set(gca,'YLabel',text('String','Work-In-System (min)','FontAngle','italic'))

figure(3)
workInSystem (T, S_v2)
set(gca, 'YLim', [-.2 Inf], 'YLim', [-.2 Inf])
set(gca,'Title',text('String',sprintf('Work-In-System for v_2 = %d ', v2),'FontAngle','italic'))
set(gca,'XLabel',text('String','Time Horizon','FontAngle','italic'))
set(gca,'YLabel',text('String','Work-In-System (min)','FontAngle','italic'))

end
%-------------------------------------------------------------------------%


%--------------------------- helper functions ----------------------------%
%-------------------------------------------------------------------------%
function [W] = computeWaitingTimes(T, S)
W = zeros(1, length(T));
for i = 2:length(W)
    W(i) = max(W(i-1) + S(i-1) - T(i),0);
end
end

function [avgWTs] = computeAverageWTs(W)
avgWTs = zeros(1, length(W));
avgWTs(1) = W(1);
for i = 2:length(W)
    avgWTs(i) = sum(W(1:i))/i;
end
end

function workInSystem (T, service_times)
%Record at every moment how much time it would take to serve everyone
% present at that moment. This is called work-in-system at time t and tells
% how long the wait would be if somebody were to arrive at time t.
% input: T - inter-arrival times, and sevice_times.

% work-in-system will have discontinuties at arrival times:
arrival_times = cumsum(T);
% work-in-system just before the arrival of the next customer:
work_just_before = zeros(1, length(arrival_times));
% work-in-system just after the arrival of the next customer:
work_just_after = zeros(1, length(arrival_times)); 
% first jump:
work_just_after(1) = service_times(1);


for i = 2:length(T)
    work_just_before(i) = max(work_just_after(i-1) - T(i), 0);
    work_just_after(i) = work_just_before(i) + service_times(i);
    
    if work_just_before(i) == 0
        t = linspace(arrival_times(i-1), arrival_times(i-1) + work_just_after(i-1));
        w = linspace(work_just_after(i-1), 0);
        plot(t, w, 'LineWidth', 1.3)
        hold on
        t = linspace(arrival_times(i-1) + work_just_after(i-1), arrival_times(i), 2);
        w = zeros(1,length(t));
        plot(t, w, 'LineWidth', 1.3)
        hold on
    else
        t = linspace(arrival_times(i-1), arrival_times(i));
        w = linspace(work_just_after(i-1), work_just_before(i));
        plot(t, w, 'LineWidth', 1.3)
        hold on
    end
    % also plot the dashed lines at jumps
    t = linspace(arrival_times(i) - 0.01, arrival_times(i) + 0.01);
    w = linspace(work_just_before(i), work_just_after(i));
    plot(t, w, 'b:');
    hold on
end
hold off
end


