% MS-E2170 Simulation
% Exercise 1.2: Monte Carlo simulation for a activity network
%
% Fill the '<your_code_here>' with your own input.
%
% Created: 2018-02-19 Heikki Puustinen


%% Initialization

% Sample size
N = <your_code_here>;

% Define the parameters of the exponential distributions, i.e., 
% mu = [mu1, mu2, ... , mun], where n is the number or activities.
mu = <your_code_here>;

% Initialize the variable where project lenghts are saved.
T = zeros(N,1);

%% Simulation

% Loop through samples
for ii = 1:N
    
    
    % Generate durations for each activity. Hint: exprnd can be used to
    % generate multiple samples from the exp-distribution at once by giving
    % it a vector as input.
    durations = <your_code_here>;
    
    % There are at least two alternative approaches to determine the
    % maximum time of the project. Use one below or make up your own. 
    % 1) Calculate durations for all possible paths and determine the
    %    maximum duration.
    % 2) Calculate the durations to reach each node by determining the 
    %    maximum of edges leading to a same node and add the maximums.
    
    
    % 1) Calculate durations for all possible path. Hint: there are 2^3 = 8
    %    possible paths in the network, e.g., [1 3 5].
    
    % Initialize durations for all possible paths
    t = zeros(8,1);
    
    t(1) = <your_code_here>;
    t(2) = <your_code_here>;
    t(3) = <your_code_here>;
    t(4) = <your_code_here>;
    t(5) = <your_code_here>;
    t(6) = <your_code_here>;
    t(7) = <your_code_here>;
    t(8) = <your_code_here>;
    

    % Determine the maximum duration of the project. Hint: maximimum
    % duration is the maximum of the individual durations of the paths.
    T(ii) = <your_code_here>;
    
    % OR
    
    % 2) Determine the maximum of edges leading to a same node and
    %    calculate their sum
%     t(1) = <your_code_here>;
%     t(2) = <your_code_here>;
%     t(3) = <your_code_here>;
%     T(ii) = <your_code_here>;

end


%% RESULTS

fprintf('Exercise 1.2 results:\n')
fprintf('Expected project duration %.5f \n',mean(T))


%% PLOTTING

% Calculate convergence
T_convergence = zeros(N,1);
for ii = 1:N
    % Calculate average of samples 1 to ii.
    T_convergence(ii) = mean(T(1:ii));
end

% Open new window etc.
f1 = figure('Name','MS-E20170, Exercise 1.2, Convergence of project duration');
ax1 = axes('Parent',f1);
hold(ax1,'on')
title(ax1,'Convergence')
xlabel(ax1,'Sample size')
ylabel(ax1,'Average project duration')

% Plot convergence
plot(ax1,1:N,T_convergence); 

% Calculate cumulative empirical distribution
T_distr = zeros(1500,1);
for ii = 1:1500
    % Determine the number of samples where the project duration is less
    % than ii.
    T_distr(ii) = sum(T < (ii/100));
end

% Open new window etc.
f2 = figure('Name','MS-E20170, Exercise 1.2, Cumulative empirical distribution');
% Move second window a bit to the right.
set(f2,'Position',get(f2,'Position') + [50 -50 0 0])
ax2 = axes('Parent',f2);
hold(ax2,'on');
title(ax2,'Cumulative empirical distribution')
xlabel(ax2,'Project duration')


% Plot cumulative distribution.
plot(ax2,(0.01:0.01:15),T_distr / N);
% Plot also the average project duration.
plot(ax2,[mean(T) mean(T)],[0 1],'r');
% Annotation
text(mean(T),0.2,['Average duration (' num2str(mean(T)) ')'],'Parent',ax2,'Color','r')