function [Q, U, R, X,resSEFF,RT_CDF,resSEFF_CDF] = CMCQN_CS_analysis_SEFF(myCQN, SEFF, SEFFgraphs, processors, max_iter, delta_max, RT, RTrange, verbose)
% [Q, U, R, X] = CMCQN_CS_ANALYSIS_SEFF(myCQN) computes the stationary performance
% measures of a Closed Multi-Class Queueing Network with Class Switching
% (CMCQNCS) myCQN. 
% More details on this type of queueing networks can be found 
% on the LINE documentation, available at http://code.google.com/p/line
% 
% Parameters: 
% myCQN:        a CMCQNCS model to analyze
% SEFF:         list of SEFF (LQN intermediate layer processors) in the model
% SEFFgraphs:   activity execution graphs for each SEFF
% processors:   list of LQN processors in the model 
% max_iter:     maximum number of iterations (termination condition)
% delta_max:    maximum difference between two succesive iterations (termination condition)
% RT:           1 if response time percentiles must be computed; 0 otherwise
% RTrange:      percentiles of the response time distribution to be  computed, if any
% verbose:      1 for screen output
% 
% Ouput:
% Q:            mean queue-length for each station and job class 
% U:            utilization for each server
% R:            response time for each job class
% X:            throughput for each job class
% resSEFF:      results (mean response time) for the SEFF (LQN soft processors) defined in the LQN model
% RT_CDF:       response time CDF for the main classes in the LQN
%
% Copyright (c) 2012-2015, Imperial College London 
% All rights reserved.

myCQNCSCox = CMCQNCS2CMCQNCSCox(myCQN);
if verbose > 1 
    disp(myCQNCSCox.toString());
end

classMatch = myCQN.classMatch;
y0 = [];
[Q, ymean_emb, ~, ~] = CMCQN_CS_fluid_ps(myCQNCSCox, max_iter, delta_max, y0, verbose);
if RT >= 1
    RT_CDF = CMCQN_CS_fluid_ps_RT(myCQNCSCox, max_iter, delta_max, ymean_emb{end}, RTrange, verbose);
    if RT == 2
        if ~isempty(SEFFgraphs)
            RT_CDF_SEFF = CMCQN_CS_fluid_ps_RT_SEFF(myCQNCSCox, max_iter, delta_max, ymean_emb{end}, SEFFgraphs, processors, RTrange, verbose);
        else
            disp('No SEFFgraphs definition has been found. Response time distributions at SEFF level will not be computed.');
        end
    else
        RT_CDF_SEFF = [];
    end
else
    RT_CDF = [];
    RT_CDF_SEFF = [];
end


%% assumes the existence of a delay node through which all classes pass
delayNodes = zeros(1,myCQN.M); 
for i = 1:myCQN.M
    if strcmp(myCQN.sched{i},'inf')
        delayNodes(i) = 1;
    end
end
delayRefNodes = zeros(1,myCQN.M); % delay nodes that are also reference nodes
for k = 1:myCQN.K
    delayRefNodes(myCQN.refNodes(k)) = 1;
end

%% computing tput - simpler version: no chain analysis - effective throughput Tput
K = myCQN.K;
rates = myCQN.rates;
refNodes = myCQN.refNodes;
S = myCQN.S;
S(refNodes(refNodes>0)) = myCQN.N*ones(sum(refNodes>0),1); 


%tput for all classes in each station
xi = sum(Q,2); %number of jobs in each station
minxi = min(xi,S); %minimum between number of jobs and servers in each station

minxi(xi>0) = minxi(xi>0)./xi(xi>0); %ratio to determine effective rate
Xfull = Q.*rates.*(minxi*ones(1,K));

%% response times + utilizations
origK = size(classMatch,1);
R = zeros(myCQN.M, myCQN.K);
for i = 1:myCQN.M
    R(i, Xfull(i,:)>0) = Q(i,Xfull(i,:)>0)./Xfull(i,Xfull(i,:)>0);
end

newR = zeros(myCQN.M,origK);
X = zeros(1,origK);
newQ = zeros(myCQN.M,origK);
% determine eventual probability of visiting each station in each class (expected number of visits)
% weight response time of each original class with the expented number of
% visits to each station in each associated artificial class
idxNR = reshape(1:myCQN.M*myCQN.K, myCQN.K, myCQN.M); %indices of non-delay (non-reference) nodes
idxNR = reshape(idxNR(:,delayRefNodes==0),1,[]); 
Ptrans = myCQN.P(idxNR,idxNR); %transient transition matrix for non-reference nodes
eventualVisit = inv( eye(size(Ptrans)) - Ptrans );
idxOrigClasses = zeros(origK,1); 
for k = 1:origK
    idxOrigClasses(k) = find(classMatch(k,:),1);
    refNode = refNodes(idxOrigClasses(k));
    
    eventualVisitProb = reshape( myCQN.P((refNode-1)*myCQN.K+idxOrigClasses(k),idxNR)*eventualVisit, myCQN.K , myCQN.M-sum(delayRefNodes>0) )'; %probability of eventual visit 
    eventualVisitProb = eventualVisitProb(:,classMatch(k,:)==1); 
    
    newR(refNode,k) = sum( R(refNode,classMatch(k,:)==1),2 );
    newR(delayRefNodes==0,k) = sum( R(delayRefNodes==0,classMatch(k,:)==1).*eventualVisitProb,  2 );
    
    X(k) = sum( Xfull(refNode,classMatch(k,:)==1) );
    newQ(:,k) = sum( Q(:,classMatch(k,:)==1),2 );
end
Rfull = R;
R = newR;
Q = newQ;

UN = zeros(myCQN.M*myCQN.K,1);
UN(myCQN.rates>0) = Xfull(myCQN.rates>0) ./ myCQN.rates(myCQN.rates>0);
UN = reshape(UN, myCQN.M, myCQN.K);

U = sum(UN,2);
U(delayNodes==0) = U(delayNodes==0)./myCQN.S(delayNodes==0);

%% Compute response times for SEFF
if ~isempty(SEFF)
    [resSEFF,resSEFF_CDF] = CMCQN_CS_respTime_SEFF(processors, SEFF, Rfull, RT_CDF_SEFF, Xfull, verbose);
else
    resSEFF = [];
    resSEFF_CDF = []; 
    %disp('No SEFF definition has been found. Performance metrics at SEFF level will not be computed.');
end