solver_ctmc_avg_from_pi.m

function [QN,UN,RN,TN,CN,XN] = solver_ctmc_avg_from_pi(sn, pivec, StateSpace, StateSpaceAggr, arvRates, depRates, options)
% 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
        % A class annihilated in place by a G-network signal leaves without a
        % service completion, so its arrival-based estimator counts offered
        % rather than carried load; use the departure-based one there.
        signalLoss = ctmc_signal_lossy(sn, arvRates, probSysState, wset, isf);
        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)
                    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) = signalLoss(k)*UNdep_ik + (1-signalLoss(k))*max(UNarv_ik,UNdep_ik);
                        end
                    end
                else % lld/cd/ljd cases
                    % Weight the per-class capacity share by the current scaling
                    % and normalize by the peak (effective capacity), as
                    % SOLVER_CTMC_ANALYZER does: the unweighted share is a
                    % P(busy)-style value, not a busy-server fraction, and would
                    % overstate the utilization of a load-dependent station.
                    ind = sn.stationToNode(ist);
                    ceff = S(ist);
                    if ~isempty(sn.lldscaling) && ist <= size(sn.lldscaling,1)
                        ceff = max(ceff, max(sn.lldscaling(ist,:)));
                    end
                    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
                            lldnow = 1;
                            if ~isempty(sn.lldscaling) && ist <= size(sn.lldscaling,1)
                                lldnow = sn.lldscaling(ist, min(max(sum(ni),1), size(sn.lldscaling,2)));
                            end
                            for k=1:K
                                UN(ist,k) = UN(ist,k) + probSysState(st)*nir(k)*sn.schedparam(ist,k)/(nir*sn.schedparam(ist,:)')*lldnow/ceff;
                            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)
                    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) = signalLoss(k)*UNdep_ik + (1-signalLoss(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
        % A signal that destroys a partially served job leaves busy time behind
        % with no completion, so T*E[S]/c is exact only for exponential service.
        % Read the busy-server occupancy off the state space instead; the two
        % agree in the exponential case. Mirrors SOLVER_CTMC_ANALYZER.
        if any(signalLoss) && sched(ist) ~= SchedStrategy.INF && ...
                isempty(sn.lldscaling) && isempty(sn.cdscaling)
            UNb = ctmc_signal_busy(sn, ind, ist, sched(ist), S(ist), StateSpace, istSpaceShift, wset, probSysState);
            UN(ist,signalLoss) = UNb(signalLoss);
        end
    end
end
 
% Stations with limited class dependence report utilization as T*S/peak,
% where peak is the user-declared peak rate scaling per class
% (sn.cdscalingpeak, from setClassDependence(beta,peakRatePerClass)). This
% matches the T*S/c convention of ordinary multiserver stations; the
% in-service accumulation above would instead report the busy probability
% P(n>=1), which is a different quantity.
if ~isempty(sn.cdscaling) && ~any(isinf(sn.njobs))
    for ist=1:M
        if length(sn.cdscaling) >= ist && ~isempty(sn.cdscaling{ist})
            for k=1:K
                bmax = sn.cdscalingpeak(ist,k);
                if isfinite(sn.rates(ist,k)) && sn.rates(ist,k) > 0 && bmax > 0
                    UN(ist,k) = TN(ist,k) / sn.rates(ist,k) / bmax;
                else
                    UN(ist,k) = 0;
                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