% Initialization
k=[ 0 0; 1 0; 2 0; 0 1; 1 1; 2 1; 0 2; 1 2; 2 2 ]; % the state space for a single node
mu1=[ 0.4 0.2 ]; % the service rates for job class 1
mu2=[ 0.3 0.4 ]; % the service rates for job class 2
mu=[ mu1; mu2 ]; % the service rates in a single matrix
e1=[ 1 1 ];      % the visit ratios for job class 1 (by hand calculation)
e2=[ 0.5 1 ];    % the visit ratios for job class 2
e=[ e1; e2 ];    % the visit ratios in a single matrix

% calculation of the function F for node 1, which is a Type-3 node (Eq. 7.82)

% "frac" has as many lines as the state space k and contains
% in its first column the quotient of the visit ratio
% and the service rate for job class 1, column 2 for job
% class 2, respectively 
frac=ones( size( k ) ) ;          
frac( :, 1 ) = e1(1) / mu1( 1 );
frac( :, 2 ) = e1(2) / mu1( 2 ); 
% calculation of the function F for a Type-3 node for the whole state space k
F1=prod( ( 1./factorial(k) .* ( frac ).^k )' )';

% "frac" for node 2 respectively (see above)
frac( :, 1 ) = e2(1)/mu2(1);
frac( :, 2 ) = e2(2)/mu2(2);

% calculation of F for the Type 2 node 2
F2=factorial( sum(k') )' .* prod( ( ( 1./factorial(k)) .* ( frac.^k ) )' )';

% F contains the values of function F: in column 1 the values for node 1
% for each possible state in k and in column 2 respectively
F=[ F1 F2 ];

% K is helper matrix which maps its (i,j)-th entry to the corresponding line
% in the state space k for a single job, when there are i jobs of class 1 
% and j jobs of class 2 in the node
K=[];
p=1;
for l = k';
  K( l(1)+1, l(2)+1 ) = p;
  p = p + 1;
end;

% calculation of the normalization constant: again the lines in G correspond to the
% states in the state space k for a single node and column 1 to the values for n=1, 
% columns 2 respectively

% for n=1 we get G1(.) = F1(.)
G=F( :, 1 );
% Initializing the values for n=2 to 0
G(:,2)=0;

% The following three loops iterate all possible values for k^(n-1), S_n and k^(n) as
% shown in Eq. 8.21, where m corresponds to k^(n-1), n to S_n and l to k^(n)
for l = k';
  for m = k';
    for n = k';
      if ( all( m + n == l ) );
	% Note that k maps the state to the corresponding number of the state
        G( K( l(1)+1, l(2)+1 ), 2 ) = G( K( l(1)+1, l(2)+1 ), 2 ) + G( K( m(1)+1, m(2)+1 ), 1 ) * F( K( n(1)+1, n(2)+1 ), 2 );
      end;
    end;
  end;
end;

fprintf( '=== Normalization constant ===\n' )
G(end,end)


% In order to calculate the marginal probabilities the values for G^(i)_N need to be calculated
% be Eq. 8.25. The values are stored in the matrix G2
G2=zeros( size( k ) );

% The values of G^(i)_N(k) are calculated by subtracting from the corresponding value of G(k)
% (see above) the sum as shown on Eq. 8.25. It make, therefore, sense to initialize the values
% to G(k) and then subtract, step by step, the values for which delta(j) (Eq. 8.25) is 1.
% Note that G_N(k) is G_2(k) in our case:
G2=[ G(:,2), G(:,2) ];
% iterate the whole state space
% l is used for k
for l = k';
  % m is used for j
  for m = k';
    % all(m <= l) prevents l from going beyond l and any( m ) checks whether m is not zero (Eq. 8.26)
    if ( all( m <= l ) && any( m ) );
      % n is used for k-j
      n=l - m;
      % subtract the corresponding value
      G2( K(l(1)+1,l(2)+1), : ) = G2( K(l(1)+1,l(2)+1), : ) - F(K(m(1)+1,m(2)+1), : ) .* G2( K(n(1)+1,n(2)+1), : );
    end;
  end;
end;

% Now the marginal probabilities can be calculated according to Eq. 8.23.
% preparation of the values: F_i(k) / G(K)
pi=F ./ G(end,end);
for l = k';
  % m stands for K-k
  m=[ 2; 2 ] - l;
  pi( K( l(1)+1, l(2)+1 ), : ) =  pi( K( l(1)+1, l(2)+1 ), : ) .* G2( K( m(1)+1, m(2)+1 ), : );
end

fprintf( '=== Marginal probabilities ===\n' )
fprintf( ' c1: number of class 1 job \n c2: number of class 2 jobs \n c3: marg. prob. node 1 \n c4: marg. prob. node 2\n' )
[ k, pi ]


% By the marginal probabilities the mean number of jobs in each system 
% can easily be derived
Kmean=[ 0 0; 0 0 ];
% Iterate the state space of a single node
for l = k'
  % Node 1
  Kmean( 1, : ) = Kmean( 1, : ) + pi( K( l(1)+1, l(2)+1 ), 1 ) .* l';
  % Node 2
  Kmean( 2, : ) = Kmean( 2, : ) + pi( K( l(1)+1, l(2)+1 ), 2 ) .* l';
end

fprintf( '=== Mean number of jobs ===\n' )
fprintf( ' c1: job class 1 \n c2: job class 2 \n l1: node 1 \n l2: node 2\n ' )
Kmean

% The throughput according to Eq. 8.28
lambda=[ 
  e(1,1) * G(8,end) / G(end,end), e(1,2) * G(6,end) / G(end,end);
  e(2,1) * G(8,end) / G(end,end), e(2,2) * G(6,end) / G(end,end) ];

fprintf( '=== Throughputs ===\n' )
fprintf( ' c1: job class 1 \n c2: job class 2 \n l1: node 1 \n l2: node 2\n ' )
lambda

% The mean response time by Little's Law
T= Kmean ./ lambda;

fprintf( '=== Mean response time ===\n' )
fprintf( ' c1: job class 1 \n c2: job class 2 \n l1: node 1 \n l2: node 2\n ' )
T

% The utilizations
rho= lambda ./ mu;

fprintf( '=== Utilizations ===\n' )
fprintf( ' c1: job class 1 \n c2: job class 2 \n l1: node 1 \n l2: node 2\n ' )
rho

