
e=[ 0.4 0.2 1 1; 0.4 0.2 1 1 ]; % visit ratios: one line per node, c1/2: open network; c3/4: closed network
e_open=e( :, 1:2 );             % visit ratios of the open network
e_closed=e( :, 3:4 );           % visit ratios of the closed network
mu=[ 1.5 1 3 2; 2 0.5 2 1 ];    % service rates: one line per node, c1/2: open network; c3/4: closed network
mu_open=mu( :, 1:2 );           % service rates of the closed network
mu_closed=mu( :, 3:4 );         % service rates of the open network
lambda_open=[ 0.2 0.1 ];        % arrival rates of the open network
K_closed=[ 2 2 ];               % number of two job classes in the closed network

% state space for a single node of the closed network
% note that this sequence is in topological order of the graph shown in Figure 8.7 for the case that
% there are only 2 jobs of job class 4 in the network
k=[ 0 0; 1 0; 0 1; 2 0; 1 1; 0 2; 2 1; 1 2; 2 2 ];


% K is a helper structure which maps its (i,j)-th entry to the line number of the corresponding
% entry in the variable k, when there i jobs of job class 1 and j jobs of class 2 in the node
K=zeros( 3, 3 );
p=1;
for l = k'
  K( l(1) + 1, l(2) + 1 ) = p;
  p = p + 1;
end
  
% STEP 1: Initialization

% calculation of the utiliziation of the open network by Eq. 8.55
utilizations_open=[ [ lambda_open; lambda_open ] .* e_open ./ mu_open ];
% sum of the utilizations of the open network, per node
utilsum_open=sum( utilizations_open' );


% Job class 3 (closed network) -- this matrices also will contain intermediate results: 
meanK3_closed=[ ]; % mean number of jobs for each node
meanT3_closed=[ ]; % mean response time for each node
lambda3= [ ];      % throughput

% the same for job class 4:
meanK4_closed=[ ];
meanT4_closed=[ ];
lambda4= [ ];

% STEP 2
% main MVA loop for the closed network
for l = k'
  % idx of the state l
  idx=K( l(1)+1, l(2)+1 );

  % Job class 3: only calculate if there is a job of class 3 in the node
  if l(1) > 0
    % ojl_idx stands for the "index of the state if there is one job less of class 3 in the network"
    % (note that Matlabs always starts with index 1, which explains why we use +1 since we are addressing with 
    % the entries of the state and there might be no job of a job class in a node)
    ojl_idx=K( l(1), l(2)+1 );

    % STEP 2.1

    % append the new values for the mean response time to the matrix (Eq. 8.53)
    meanT3_closed = [ meanT3_closed; ( 1 ./ mu_closed( :, 1 )' ) .* ( 1 + meanK3_closed( ojl_idx, : ) + meanK4_closed( ojl_idx, : ) ) ./ ( 1 - utilsum_open ) ];

    % STEP 2.2

    % calculate the troughput (Eq. 8.46)
    lambda3=[ lambda3; l(1) ./ ( sum( e_closed( :, 1 )' .* meanT3_closed( end, : ) ) ) ];

    % STEP 2.3

    % calculate the mean number of jos (Eq. 8.47)
    meanK3_closed = [ meanK3_closed; lambda3(end) .* e_closed( :, 1 )' .* meanT3_closed(end,:) ]; 

  else
    % if there is _no_ job of class 3 in the node all performance measures are 0
    meanT3_closed = [ meanT3_closed; [ 0, 0 ] ];
    meanK3_closed = [ meanK3_closed; [ 0, 0 ] ];
    lambda3 = [ lambda3; 0 ];
  end

  % ...and respectively for job class 4
  if l(2) > 0
    ojl_idx=K( l(1)+1, l(2) );

    % STEP 2.1

    meanT4_closed = [ meanT4_closed; ( 1 ./ mu_closed( :, 2 )' ) .* ( 1 + meanK3_closed( ojl_idx, : ) + meanK4_closed( ojl_idx, : ) ) ./ ( 1 - utilsum_open ) ];

    % STEP 2.2

    lambda4=[ lambda4; l(2) ./ ( sum( e_closed( :, 2 )' .* meanT4_closed( end, : ) ) ) ];

    % STEP 2.3

    meanK4_closed = [ meanK4_closed; lambda4(end) .* e_closed( :, 2 )' .* meanT4_closed(end,:) ];

  else
    meanT4_closed = [ meanT4_closed; [ 0, 0 ] ];
    meanK4_closed = [ meanK4_closed; [ 0, 0 ] ];
    lambda4 = [ lambda4; 0 ];
  end
    
end

% STEP 3

% compute the throughputs of the closed network
lambdas_closed = [ lambda3(end,:), lambda4(end,:); lambda3(end,:), lambda4(end,:) ] .* e_closed;

% Compute the mean number of jobs for each node in the open network (Eq. 8.51)
meanK_open= [ ( utilizations_open( :, 1 ) .* ( 1 + meanK3_closed(end,:)' + meanK4_closed( end,:)' ) ./ ( 1 - utilsum_open' ) )';
              ( utilizations_open( :, 2 ) .* ( 1 + meanK3_closed(end,:)' + meanK4_closed( end,:)' ) ./ ( 1 - utilsum_open' ) )' ];

lambdas_open = e_open .* [ lambda_open; lambda_open ];

% Little's Law
meanT_open = meanK_open ./ lambdas_open;

% OUTPUT

fprintf( '=== Mean response times ===\n' );
fprintf( ' one column per job class / one line per node\n' )
[ meanT_open', meanT3_closed(end,:)', meanT4_closed(end,:)' ]

fprintf( '=== Mean number of jobs ===\n' );
fprintf( ' one column per job class / one line per node\n' )
[ meanK_open', meanK3_closed(end,:)', meanK4_closed(end,:)' ]

fprintf( '=== Throughputs ===\n' );
fprintf( ' one column per job class / one line per node\n' )
[ lambdas_open, lambdas_closed ]

