

N=6; % number of nodes
R=3; % number of job classes
K=[ 1 1 1 ]; % number of jobs per job class
maxK = max( K );

% NOTE: this script currently only works for K=[1 1 1]

% visit ratios: columns = nodes / rows = job classes
e= ...
  [ 2 0.5 0.5 0 0 2;
    2 0   0   1 0 2;
    2 0   0   0 1 2 ];

% service rates: columns = nodes / rows = job classes
mu = 1./ [ 7.5 28 12 30 15 2.5;
           7.5 28 12 30 15 2.5;
           7.5 28 12 30 15 2.5; ];

% number of servers: columns = nodes / rows = job classes (just for convenience)
m = [ 2 2 1 1 1 8,
      2 2 1 1 1 8;
      2 2 1 1 1 8 ];

% the following loop calculates the sequence of population vectors
% for the MVA iteration
left = 1;
right = 2;
ks=[ 0 0 0];
while left < right
  ks=[ ks; Inf Inf Inf ];
  for k = ks((left:right),:)'
    k=k';
    for r = (1:R)
      if k(r) + 1 <= K(r)
        k_temp = k;
        k_temp(r) = k_temp(r) + 1;
        [ rows, cols ] = size(ks);
        if ( all( sum( ones(rows,1) * k_temp ~= ks, 2 ) ) )
          right = right + 1;
          ks = [ ks; k_temp ];
        end
      end
    end
    left = left + 1;
  end
end
[ rows, cols ] = size(ks);
ks = ks( find( ~sum( ones(rows,1) * [ Inf, Inf, Inf ] == ks, 2) ), : )


% STEP 1
% mean number of jobs and mean response time (3-dim matrices):
%   1st dim: nodes
%   2nd dim: job classes
%   3nd dim: population vector (via unique number: see k2num)
meanK = zeros( N, R, (maxK+1)^R );
meanT = zeros( N, R, (maxK+1)^R );
% overall throughput (2-dim matrix)
%   1st dim: job classes
%   2nd dim: population vector (via unique number: see k2num)
lambdas = zeros( R, (maxK+1)^R );
% probabilities (3-dim matrix): pi_i( j | v )
%   1st dim: nodes --> i
%   2nd dim: number of jobs in node --> j
%   3nd dim: population vector (via unique number: see k2num) --> v
pi = zeros( N, maxK + 1, (maxK+1)^R );
pi( :, 0 +1, 0 +1 ) = 1;

% STEP 2
% Iterate all population vectors
for  k = ks'

  % take the transpose
  k = k';
  % number of the current population vector
  kNum = k2num( k, R, maxK );

  % 0 population?
  if kNum == 0
    % already handled by STEP 1
    continue
  end

  % calculate the population vectors of the current population with one
  % job less for all job classes
  % actually we need the address of these population vectors
  kOneLessNum=[];
  for r = (0:R-1)
    kOneLessNum = [ kOneLessNum, max( kNum - (maxK+1)^r, 0 ) ];
  end

  % STEP 2.1
  % calculate the mean response times
  meanT( :, :, kNum + 1 ) = ( 1./ ( m .* mu ) )' .*  ...
     ( 1 + reshape( sum( meanK( :, :,  kOneLessNum + 1 ), 2 ), N, R ) ...
     + (m-1)' .* reshape( pi( :, 0+1, kOneLessNum+1 ), N, R ) );


  % node 6 is a special case: it is a M/M/Infinity server
  meanT( 6, :, kNum + 1 ) = 1/mu(1,6);

  % STEP 2.2
  % calculate the throughput
  lambdas( :, kNum + 1 ) = k ./ diag( e * meanT(:,:,kNum+1) )';

  % STEP 2.1 (again)
  % calculate the probabilities
  pi( :, (2:end), kNum+1 ) = (1./(ones( 1, N )' * (1:maxK)))...
   .* sum(( (e./mu .* ( lambdas( :, kNum + 1 ) * ones(1,N)))' .* reshape( pi( :, (1:end-1), kOneLessNum+1 ), N, R ) ), 2 );
  pi( :, 1, kNum+1 ) = 1 - (1./m(1,:))' .* ( sum( (e./mu .* ( lambdas(:, kNum + 1 ) * ones(1,N) ))' , 2 ) + sum( pi( :, (2:end), kNum+ 1 ), 2 ) );

  % STEP 2.3
  % calculate the mean number of jobs
  meanK(:,:,kNum+1) = (ones( N, 1 ) * lambdas(:,kNum+1)') .* e' .* meanT(:,:,kNum+1);

  % some output
  fprintf( '=== Results ===\n' );
  fprintf( 'Population vector:\n' );
  k
  fprintf( 'Probabilities conditioned on current population: rows = nodes / columns = number of jobs in node\n' );
  fprintf( 'NOTE: Negative probabilities are _no_ calculation errors; they are calculated unintentional\n' );
  fprintf( '      and are not needed in the current iteration\n' );
  pi( :, :, kNum+1)
  fprintf( 'Mean response times: rows = nodes / columns = job classes\n' );
  meanT(:,:,kNum+1)
  fprintf( 'Throughputs: rows = job classes\n' );
  lambdas(:,kNum+1)
  fprintf( 'Mean number of jobs: rows = nodes / columns = job classes\n' );
  meanK(:,:,kNum+1)
  
end
