% Network of Problem 7.3
% Parameters of the network
K = 4;
N = 3;
mu = 1./[ 2 5 5 ];
m = [ 1 1 1 ];
e = [ 1 0.5 0.2 ];

% === INITIALIZATION ===
temp = ( mu .* m ) ./ e;
lambda = [ min( temp ) ];
bott = min( find( temp == lambda ) );
lambdas = [];
rho = [];
meanK = [];
g = [];
epsilon = 0.02;
err = [ 2 * epsilon ];
corr = [ ];
maxIterations = 10;
iterations = 0;
P_12 = []

fprintf( '==========================================\n' );
fprintf( '=== Results for network of Problem 7.3 ===\n' );
fprintf( '==========================================\n' );
% STEP 2
printf( 'Intermediate results for the bott node %d:\n', bott );
while err(end) > epsilon && iterations < maxIterations

  % STEP 2.1
  lambdas = [ lambdas; lambda(end) .* e ];
  rho = [ rho; lambdas(end,:) ./ ( m .* mu ) ];
  meanK = [ meanK; rho(end,:) ./ ( 1 - (K-1)/K .* rho(end,:) ) ];

  printf( '  rho_%d = %f, meanK_%d = %f\n', bott, rho(end,bott), bott, meanK(end,bott) );

  % STEP 2.2
  g = [ g; sum( meanK(end,:) ) ]; 
  corr = [ corr; K/g(end) ];
  meanK(end,bott) = meanK(end,bott) * corr(end);
  rho(end,bott) = meanK(end,bott) / ( 1 + (K-1)/K * meanK(end,bott) );



  lambda = [ lambda; ( rho(end,bott) * m(bott) * mu(bott) ) / e(bott) ];

  if rho(end,bott) > 1
    fprintf( 'The value of rho_bott is greater than 1: Terminating iteration.\n' );
    rho(end,bott) = 1;
    break;
  end
  
  iterations = iterations + 1;

  % STEP 2.3
  err = [ err; abs( K - g(end) ) ];

end
% STEP 3
lambdas = [ lambdas; lambda(end) * e ];
meanK = [ meanK; meanK(end,:) .* K/g(end) ];
meanT = meanK ./ lambdas;
g = [ g; sum( meanK(end,:) ) ];

% OUTPUT

fprintf( '= Final Results ==========================\n' );
fprintf( ' ITERATIONS: %d\n', iterations );
fprintf( ' BOTT Node: %d\n', bott );
fprintf( '=== Mean response times: meanT_i ===\n' );
fprintf( ' rows = iteration / columns = nodes (1 - 3)\n ' );
meanT
fprintf( '=== Mean number of jobs: meanK_i ===\n' );
fprintf( ' rows = iteration / columns = nodes (1 - 3), column 4: g(lambda) \n ' );
[ meanK, g(1:end,:) ]
fprintf( '=== Utilizations: rho_i ===\n' );
fprintf( ' rows = iteration / columns = nodes (1 - 3)\n ' );
rho
fprintf( '=== Throughputs: lambda and lambda_i ===\n' );
fprintf( ' columns = column one: overall throughput / columns 2-4: nodes 1-3\n ' );
[ lambda, lambdas ]
fprintf( '=== Correction factor: corr ===\n' );
fprintf( ' rows = iterations\n ' );
corr
fprintf( '=== Error ===\n' );
fprintf( ' rows = iterations\n ' );
err

% Network of Problem 7.4
% Parameters of the network
K = 3;
N = 3;
mu = [ 0.4 0.6 0.3 ];
m = [ 2 1 1 ];
e = [ 1 0.514 0.357 ];

% === INITIALIZATION ===
temp = ( mu .* m ) ./ e;
lambda = [ min( temp ) ];
bott = min( find( temp == lambda ) );
lambdas = [];
rho = [];
meanK = [];
g = [];
epsilon = 0.02;
err = [ 2 * epsilon ];
corr = [ ];
maxIterations = 10;
iterations = 0;

fprintf( '==========================================\n' );
fprintf( '=== Results for network of Problem 7.4 ===\n' );
fprintf( '==========================================\n' );
% STEP 2
printf( 'Intermediate results for the bott node %d:\n', bott );
while err(end) > epsilon && iterations < maxIterations

  % STEP 2.1
  lambdas = [ lambdas; lambda(end) .* e ];
  rho = [ rho; lambdas(end,:) ./ ( m .* mu ) ];
  meanK = [ meanK; rho(end,:) ./ ( 1 - (K-1)/K .* rho(end,:) ) ];

  % Calculation of meanK_1 (M/M/2-FCFS node)
  P_12 = [ P_12; 2 * rho(end,1)^2/(1+rho(end,1)) ];
  meanK(end,1) = m(1) * rho(end,1) + rho(end,1)/(1 - ( K - m(1) - 1 )/(K - m(1))) * P_12(end);

  printf( '  rho_%d = %f, meanK_%d = %f\n', bott, rho(end,bott), bott, meanK(end,bott) );

  % STEP 2.2
  g = [ g; sum( meanK(end,:) ) ]; 
  corr = [ corr; K/g(end) ];
  meanK(end,bott) = meanK(end,bott) * corr(end);

  % We use the fact that bott = 1 !
  rho(end,bott) = meanK(end,bott)/(m(bott) + P_12(end));

  lambda = [ lambda; ( rho(end,bott) * m(bott) * mu(bott) ) / e(bott) ];

  if rho(end,bott) > 1
    fprintf( 'The value of rho_bott is greater than 1: Terminating iteration.\n' );
    rho(end,bott) = 1;
    break;
  end
  
  iterations = iterations + 1;

  % STEP 2.3
  err = [ err; abs( K - g(end) ) ];

end
% STEP 3
lambdas = [ lambdas; lambda(end) * e ];
meanK = [ meanK; meanK(end,:) .* K/g(end) ];
meanT = meanK ./ lambdas;
g = [ g; sum( meanK(end,:) ) ];

% OUTPUT

fprintf( '= Final Results ==========================\n' );
fprintf( ' ITERATIONS: %d\n', iterations );
fprintf( ' BOTT Node: %d\n', bott );
fprintf( '=== Mean response times: meanT_i ===\n' );
fprintf( ' rows = iteration / columns = nodes (1 - 3)\n ' );
meanT
fprintf( '=== Mean number of jobs: meanK_i ===\n' );
fprintf( ' rows = iteration / columns = nodes (1 - 3), column 4: g(lambda) \n ' );
[ meanK, g(1:end,:) ]
fprintf( '=== Utilizations: rho_i ===\n' );
fprintf( ' rows = iteration / columns = nodes (1 - 3)\n ' );
rho
fprintf( '=== Throughputs: lambda and lambda_i ===\n' );
fprintf( ' columns = column one: overall throughput / columns 2-4: nodes 1-3\n ' );
[ lambda, lambdas ]
fprintf( '=== Correction factor: corr ===\n' );
fprintf( ' rows = iterations\n ' );
corr
fprintf( '=== Error ===\n' );
fprintf( ' rows = iterations\n ' );
err
