solver_ctmc_avg_from_pi.m

function [QN,UN,RN,TN,CN,XN] = solver_ctmc_avg_from_pi(sn, pivec, StateSpace, StateSpaceAggr, arvRates, depRates)
% SOLVER_CTMC_AVG_FROM_PI  Map a state distribution to mean performance metrics.
%
% [QN,UN,RN,TN,CN,XN] = SOLVER_CTMC_AVG_FROM_PI(SN, PIVEC, STATESPACE,
%   STATESPACEAGGR, ARVRATES, DEPRATES)
%
% Given an arbitrary probability vector PIVEC over the enumerated CTMC state
% space of SN (rows of STATESPACE / STATESPACEAGGR), returns the per-(station,
% class) mean queue length QN, utilization UN, response time RN, throughput TN,
% system response time CN and system throughput XN. The discipline-aware mapping
% is identical to the steady-state reduction performed by solver_ctmc_analyzer;
% it is factored here so that callers holding their own distribution (e.g. the
% SolverENV state-vector analyzer, which time-averages a transient distribution)
% can reuse it without re-solving for the stationary vector.
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
M = sn.nstations;
K = sn.nclasses;
S = sn.nservers;
NK = sn.njobs';
sched = sn.sched;
PH = sn.proc;
 
probSysState = pivec(:)';
probSysState(probSysState<GlobalConstants.Zero) = 0;
if sum(probSysState) > 0
    probSysState = probSysState/sum(probSysState);
end
wset = 1:size(StateSpace,1);
 
XN = NaN*zeros(1,K);
UN = NaN*zeros(M,K);
QN = NaN*zeros(M,K);
RN = NaN*zeros(M,K);
TN = NaN*zeros(M,K);
CN = NaN*zeros(1,K);
 
istSpaceShift = zeros(1,M);
for ist=1:M
    if ist==1
        istSpaceShift(ist) = 0;
    else
        istSpaceShift(ist) = istSpaceShift(ist-1) + size(sn.space{ist-1},2);
    end
end
 
for k=1:K
    refsf = sn.stationToStateful(sn.refstat(k));
    XN(k) = probSysState*arvRates(wset,refsf,k);
end
 
for ist=1:M
    isf = sn.stationToStateful(ist);
    ind = sn.stationToNode(ist);
    for k=1:K
        TN(ist,k) = probSysState*depRates(wset,isf,k);
        QN(ist,k) = probSysState*StateSpaceAggr(wset,(ist-1)*K+k);
    end
    if sn.nodetype(ind) ~= NodeType.Source
        switch sched(ist)
            case SchedStrategy.INF
                for k=1:K
                    UN(ist,k) = QN(ist,k);
                end
            case {SchedStrategy.PS, SchedStrategy.DPS, SchedStrategy.GPS, SchedStrategy.LPS}
                if isempty(sn.lldscaling) && isempty(sn.cdscaling) && ~sn_has_joint_dependence(sn)
                    for k=1:K
                        if ~isempty(PH{ist}{k})
                            % There are cases where due to remove of
                            % immediate transitions or due to cutoff the
                            % utilization estimator based on arrivals or
                            % departures can be under-estimated. E.g.:
                            % UNarv_ik doesn't work well with example_stateDependentRouting_3
                            % UNdep_ik doesn't work well with test_OQN_JMT_6
                            % Therefore, we take the maximum of the two
                            % Note: the two estimators are normally
                            % identical.
                            UNarv_ik = probSysState*arvRates(wset,isf,k)*map_mean(PH{ist}{k})/S(ist);
                            UNdep_ik = TN(ist,k)*map_mean(PH{ist}{k})/S(ist); % this is valid because CS in LINE is in a separate node
                            UN(ist,k) = max(UNarv_ik,UNdep_ik);
                        end
                    end
                else % lld/cd/ljd cases
                    ind = sn.stationToNode(ist);
                    UN(ist,1:K) = 0;
                    for st = wset
                        [ni,nir] = State.toMarginal(sn, ind, StateSpace(st,(istSpaceShift(ist)+1):(istSpaceShift(ist)+size(sn.space{ist},2))));
                        if ni>0
                            for k=1:K
                                UN(ist,k) = UN(ist,k) + probSysState(st)*nir(k)*sn.schedparam(ist,k)/(nir*sn.schedparam(ist,:)');
                            end
                        end
                    end
                end
            case SchedStrategy.PAS
                % Pass-and-swap (order-independent): utilization is the
                % time-average number of in-service jobs per class divided by
                % the number of servers. sir already counts the positions whose
                % marginal rate increment is positive (toMarginal PAS branch),
                % so a single job served by multiple server types still counts
                % as one in-service job (not 1/rate).
                ind = sn.stationToNode(ist);
                UN(ist,1:K) = 0;
                for st = wset
                    [~,~,sir] = State.toMarginal(sn, ind, StateSpace(st,(istSpaceShift(ist)+1):(istSpaceShift(ist)+size(sn.space{ist},2))));
                    for k=1:K
                        UN(ist,k) = UN(ist,k) + probSysState(st)*sir(k)/S(ist);
                    end
                end
            otherwise
                if isempty(sn.lldscaling) && isempty(sn.cdscaling) && ~sn_has_joint_dependence(sn)
                    for k=1:K
                        if ~isempty(PH{ist}{k})
                            % There are cases where due to remove of
                            % immediate transitions or due to cutoff the
                            % utilization estimator based on arrivals or
                            % departures can be under-estimated. E.g.:
                            % UNarv_ik doesn't work well with example_stateDependentRouting_3
                            % UNdep_ik doesn't work well with test_OQN_JMT_6
                            % Therefore, we take the maximum of the two
                            % Note: the two estimators are normally
                            % identical.
                            UNarv_ik = probSysState*arvRates(wset,isf,k)*map_mean(PH{ist}{k})/S(ist);
                            UNdep_ik = TN(ist,k)*map_mean(PH{ist}{k})/S(ist); % this is valid because CS in LINE is in a separate node
                            UN(ist,k) = max(UNarv_ik,UNdep_ik);
                        end
                    end
                else % lld/cd/ljd cases
                    ind = sn.stationToNode(ist);
                    UN(ist,1:K) = 0;
                    for st = wset
                        [ni,~,sir] = State.toMarginal(sn, ind, StateSpace(st,(istSpaceShift(ist)+1):(istSpaceShift(ist)+size(sn.space{ist},2))));
                        if ni>0
                            for k=1:K
                                UN(ist,k) = UN(ist,k) + probSysState(st)*sir(k)/S(ist);
                            end
                        end
                    end
                end
        end
    end
end
 
for k=1:K
    for ist=1:M
        if TN(ist,k)>0
            RN(ist,k) = QN(ist,k)./TN(ist,k);
        else
            RN(ist,k)=0;
        end
    end
    CN(k) = NK(k)./XN(k);
end
 
QN(isnan(QN))=0;
CN(isnan(CN))=0;
RN(isnan(RN))=0;
UN(isnan(UN))=0;
XN(isnan(XN))=0;
TN(isnan(TN))=0;
end