solver_ctmc_analyzer.m

function [QN,UN,RN,TN,CN,XN,InfGen,StateSpace,StateSpaceAggr,EventFiltration,runtime,fname,sncopy] = solver_ctmc_analyzer(sn, options)
% [QN,UN,RN,TN,CN,XN,INFGEN,STATESPACE,STATESPACEAGGR,EVENTFILTRATION,RUNTIME,FNAME,sn] = SOLVER_CTMC_ANALYZER(sn, OPTIONS)
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
%if options.remote
%    sn.rtfun = {};
%    sn.lst = {};
%    qn_json = jsonencode(sn);
%    sn = NetworkStruct.fromJSON(qn_json)
%return
%end
 
M = sn.nstations;    %number of stations
K = sn.nclasses;    %number of classes
S = sn.nservers;
NK = sn.njobs';  % initial population per class
sched = sn.sched;
 
Tstart = tic;
PH = sn.proc;
 
line_debug('CTMC analyzer starting: nstations=%d, nclasses=%d, njobs=%s', M, K, mat2str(NK));
 
% Note: hide_immediate now selectively preserves Cache immediate transitions
% in solver_ctmc.m, so we no longer need to disable it entirely for Cache nodes
 
line_debug('Building state space and infinitesimal generator via solver_ctmc');
[InfGen,StateSpace,StateSpaceAggr,EventFiltration,arvRates,depRates,sn] = solver_ctmc(sn, options); % sn is updated with the state space
 
% if the initial state does not reflect the final size of the state
% vectors, attempt to correct it
for isf=1:sn.nstateful
    if size(sn.state{isf},2) < size(sn.space{isf},2)
        sn.state{isf} = [zeros(1,size(sn.space{isf},2)-size(sn.state{isf},2)),sn.state{isf}];
    end
end
sncopy = sn;
 
if options.keep
    line_debug('Saving CTMC data to file (options.keep=true)');
    fname = lineTempName;
    save([fname,'.mat'],'InfGen','StateSpace','StateSpaceAggr','EventFiltration')
    line_printf('CTMC infinitesimal generator and state space saved in: ');
    line_printf(strrep(sprintf('%s.mat\n',fname),'\','\\'))
else
    fname = '';
end
 
wset = 1:length(InfGen);
 
line_debug('State space built: %d states, solving CTMC', length(InfGen));
 
use_ctmc_solve_stable = true;
if use_ctmc_solve_stable
    % stable version
    % Note: solver_ctmc now selectively preserves Cache immediate transitions
    % to enable hit/miss rate computation while hiding other immediate transitions
    [probSysState, ~, nConnComp, connComp] = ctmc_solve(InfGen, options);
 
    if nConnComp > 1
        line_debug('CTMC is reducible: %d connected components', nConnComp);
        % the matrix was reducible
        initState = matchrow(StateSpace, cell2mat(sn.state'));
        if initState <= 0
            % Initial state may have been removed by stochcomp (e.g., SPN with
            % immediate ENABLE states). Use the largest connected component.
            compSizes = accumarray(connComp(:), 1);
            [~, largestComp] = max(compSizes);
            wset = find(connComp == largestComp);
        else
            % determine the weakly connected component associated to the initial state
            wset = find(connComp == connComp(initState));
        end
        if initState > 0
            line_debug('Using component %d with %d states (from initial state)', connComp(initState), length(wset));
        else
            line_debug('Using largest component with %d states (initial state removed by stochcomp)', length(wset));
        end
        probSysState = ctmc_solve(InfGen(wset, wset), options);
        InfGen = InfGen(wset, wset);
        % reduce all per-state arrays to the retained component and remap wset
        % to local indices so matrix-form and loop-form estimators stay aligned
        StateSpace = StateSpace(wset,:);
        StateSpaceAggr = StateSpaceAggr(wset,:);
        arvRates = arvRates(wset,:,:);
        depRates = depRates(wset,:,:);
        wset = 1:numel(wset);
    else
        line_debug('CTMC is irreducible, using full state space');
    end
else
    % development version
 
    % we now find the initial state and then solver the CTMC allowing for the
    % case where it is reducible
    initState = matchrow(StateSpace, cell2mat(sn.state'));
    pi0 = zeros(1,length(InfGen)); pi0(initState) = 1.0;
    [pi,pis,~,scc,~] = ctmc_solve_reducible(InfGen, pi0, options);
 
    if size(pis,1)==1
        probSysState = pi;
    else
        wset = scc == scc(initState);
        InfGen = InfGen(wset, wset);
        StateSpace = StateSpace(wset,:);
        probSysState = pis(scc(initState),scc == scc(initState));
    end
end
probSysState(probSysState<GlobalConstants.Zero)=0;
probSysState = probSysState/sum(probSysState);
 
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
 
% Stations that participate in a finite capacity region with a loss (DROP)
% rule can shed offered load exactly like a station-level capacity limit, so
% for such stations only the departure-based (carried) utilization estimator
% is valid. Precompute per-station membership in any DROP region.
inDropRegion = false(1,M);
if isfield(sn,'nregions') && sn.nregions > 0
    for f=1:sn.nregions
        if sn.regionrule(f) == DropStrategy.DROP
            memb = any(sn.region{f} ~= -1, 2); % stations constrained by region f
            memb = memb(:)';
            inDropRegion(1:min(M,numel(memb))) = inDropRegion(1:min(M,numel(memb))) | memb(1:min(M,numel(memb)));
        end
    end
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
        % For a capacity-constrained OPEN class arrivals can be dropped, so the
        % arrival-based estimator UNarv_ik counts the OFFERED load (including
        % dropped jobs) and overestimates utilization; only the departure-based
        % (carried load) estimator is valid for such a class. Closed classes are
        % never dropped (their classcap is finite by population, not by a buffer).
        % A station inside a DROP-rule finite capacity region can also shed load
        % even when its own station capacity is infinite (the limit is on the
        % region), so treat its open classes as droppable too.
        canDropClass = isinf(sn.njobs(:)') & (isfinite(sn.cap(ist)) | isfinite(sn.classcap(ist,:)) | inDropRegion(ist));
        % G-network signals annihilate jobs in place: a removed job leaves
        % without ever completing service, so UNarv_ik again counts offered
        % rather than carried load (M/M/1 with lambda+=0.5, lambda-=0.4, mu=1
        % gives UNarv=0.5 against the exact rho=lambda+/(mu+lambda-)=0.35714).
        % Only the departure-based estimator is valid for a class exposed to a
        % signal at this station.
        signalLoss = ctmc_signal_lossy(sn, arvRates, probSysState, wset, isf);
        canDropClass = canDropClass | signalLoss;
        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.
                            UNdep_ik = TN(ist,k)*map_mean(PH{ist}{k})/S(ist); % this is valid because CS in LINE is in a separate node
                            if canDropClass(k)
                                UN(ist,k) = UNdep_ik;
                            else
                                UNarv_ik = probSysState*arvRates(wset,isf,k)*map_mean(PH{ist}{k})/S(ist);
                                UN(ist,k) = max(UNarv_ik,UNdep_ik);
                            end
                        end
                    end
                else % lld/cd/ljd cases
                    % Busy-server FRACTION under load-dependent scaling (NC
                    % convention): weight the per-class capacity share by the
                    % current scaling and normalize by the peak (effective
                    % capacity); the unweighted share is a P(busy)-style value.
                    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.
                            UNdep_ik = TN(ist,k)*map_mean(PH{ist}{k})/S(ist); % this is valid because CS in LINE is in a separate node
                            if canDropClass(k)
                                UN(ist,k) = UNdep_ik;
                            else
                                UNarv_ik = probSysState*arvRates(wset,isf,k)*map_mean(PH{ist}{k})/S(ist);
                                UN(ist,k) = max(UNarv_ik,UNdep_ik);
                            end
                        end
                    end
                else % lld/cd/ljd cases
                    % Busy-server FRACTION under load-dependent scaling (NC
                    % convention): weight the in-service indicator by the current
                    % scaling and normalize by the peak (effective capacity);
                    % plain sir/S collapses to P(busy).
                    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,~,sir] = 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
                            sirtot = sum(sir);
                            for k=1:K
                                if sirtot > 0
                                    UN(ist,k) = UN(ist,k) + probSysState(st)*(sir(k)/sirtot)*lldnow/ceff;
                                end
                            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
        % (M/Er2/1 with lambda+=0.5, lambda-=0.4 gives 0.34941 against a true
        % 0.37696). Read the busy-server occupancy off the state space instead;
        % it agrees with T*E[S]/c in the exponential case. The lld/cd branches
        % already accumulate in-service occupancy, so only the plain ones need
        % it, and INF stations report the queue length by definition.
        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), matching the T*S/c convention of ordinary multiserver
% stations; the in-service accumulation above would instead report P(n>=1).
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
 
% True BAS queue-length attribution. A job that completed service at a BAS station
% and is held there (blocked marker = 1) awaiting room downstream is ATTRIBUTED TO
% THE DESTINATION, not to the station it physically occupies. That is the JMT
% convention, and it is what LDES already implements (Solver_ssj.hasBASBlocking
% declares the policy on the destination queue), so CTMC was the odd one out inside
% LINE; aligning removes an internal inconsistency rather than deferring to JMT.
%
% This is a REPORTING change only: the chain, the state space and the become-blocked
% edge are untouched. The shift is exactly P(blocked) per station.
%
% Metric conventions, each verified against JMT and LDES on cqn_bas_blocking
% (N=2, Q1 Exp(1), Q2 Exp(2), Q2 cap=1; exact X=6/7, P(blocked)=1/7):
%   QLen  : blocked job counted at the DESTINATION. Q1 11/7->10/7=1.428571 (JMT
%           1.42850, LDES 1.42710), Q2 3/7->4/7=0.571429 (JMT 0.57150, LDES
%           0.57286).
%   RespT : follows from QN by Little below (RN=QN/TN), giving R1=10/6=1.666667
%           (JMT 1.66795, LDES 1.66610) and R2=4/6=0.666667 (JMT 0.66687, LDES
%           0.66873). Both are self-consistent under Little, which therefore cannot
%           arbitrate the convention -- it is a placement choice, not physics.
%   Util  : UNCHANGED. The blocked job holds the server but is not being served, so
%           it does not add to utilization: all engines already satisfy U=X*E[S]
%           (U1=6/7=0.857143 vs JMT 0.85735/LDES 0.85692; U2=3/7=0.428571 vs JMT
%           0.42942/LDES 0.42978). Util is computed from the service process above,
%           not from QN, so it is untouched by construction.
%
% Attribution needs the destination of the blocked job. The marker records only THAT
% a job is blocked, not for which destination, so this is applied only when the
% station has exactly ONE downstream station, which is the case for every BAS model
% in the tree (tandem/cyclic). With several possible destinations the state cannot
% say which one is meant, so the job is left where it sits rather than guessed.
if ~isempty(sn.isbasblocking)
    for ist=1:M
        ind = sn.stationToNode(ist);
        % Identify the blocking station by the dedicated sn.isbasblocking field
        % (set by refreshLocalVars for exactly the blocking/upstream stations,
        % under BOTH the upstream and destination declaration forms), not by
        % re-testing its own BAS drop rule -- the latter fails for a
        % destination-declared BAS and left its QLen in the where-it-sits
        % convention. The field also disambiguates the shared marker column from
        % the polling controller. See BUG-83.
        if numel(sn.isbasblocking) < ind || sn.isbasblocking(ind) ~= 1
            continue % no blocked marker at this station
        end
        dests = find(sn.connmatrix(ind,:) == 1 & sn.isstation(:)' == 1);
        if numel(dests) ~= 1
            continue % ambiguous destination: leave the job where it sits
        end
        jst = sn.nodeToStation(dests(1));
        if isnan(jst) || jst < 1
            continue
        end
        cols = (istSpaceShift(ist)+1):(istSpaceShift(ist)+size(sn.space{ist},2));
        blocked = StateSpace(wset, cols(end)) == 1; % marker is the trailing column
        for k=1:K
            % Only the held job itself moves, not the whole queue at ist: a blocked
            % state holds exactly one completed job, so cap the per-state count at 1.
            shift = probSysState(blocked) * min(StateSpaceAggr(wset(blocked),(ist-1)*K+k), 1);
            if shift > 0
                QN(ist,k) = QN(ist,k) - shift;
                QN(jst,k) = QN(jst,k) + shift;
            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;
 
runtime = toc(Tstart);
 
% now update the routing probabilities in nodes with state-dependent routing
TNcache = zeros(sn.nstateful,K);
XNcache = zeros(sn.nstateful,K);
for k=1:K
    for isf=1:sn.nstateful
        ind = sncopy.statefulToNode(isf);
        if sncopy.nodetype(ind) == NodeType.Cache
            TNcache(isf,k) = probSysState*depRates(wset,isf,k);
            XNcache(isf,k) = probSysState*arvRates(wset,isf,k);
        end
    end
end
 
% updates cache actual hit and miss data + retrieval-system expected latency.
% Hit / miss for class k are derived from departure rates of the configured
% hitClass / missClass at the cache; for retrieval-aware caches the
% miss-class departure rate equals the true miss rate because
% afterEventCache fires the miss event on a retrieval-complete READ.
retrievalLatencyWarned = false;
for k=1:K
    for isf=1:sncopy.nstateful
        ind = sncopy.statefulToNode(isf);
        if sncopy.nodetype(ind) == NodeType.Cache
            np = sncopy.nodeparam{ind};
            if length(np.hitclass)>=k
                h = np.hitclass(k);
                m = np.missclass(k);
                if h>0 && m>0
                    sncopy.nodeparam{ind}.actualhitprob(k) = TNcache(isf,h)/sum(TNcache(isf,[h,m]));
                    sncopy.nodeparam{ind}.actualmissprob(k) = TNcache(isf,m)/sum(TNcache(isf,[h,m]));
                else
                    sncopy.nodeparam{ind}.actualhitprob(k) = NaN;
                    sncopy.nodeparam{ind}.actualmissprob(k) = NaN;
                end
 
                % The Eq. 8 retrieval-system expected latency is not currently
                % implemented; report NaN whenever a retrieval system is
                % configured for this class.
                expectedLatency = NaN;
                if isfield(np, 'retrievalSystemQueueIndices') ...
                        && isKey(np.retrievalSystemQueueIndices, int32(k-1)) ...
                        && ~isempty(np.retrievalSystemQueueIndices(int32(k-1)))
                    if ~retrievalLatencyWarned
                        line_warning(mfilename, 'Retrieval-system expected latency is not currently implemented; reporting NaN.');
                        retrievalLatencyWarned = true;
                    end
                end
                sncopy.nodeparam{ind}.actualresidt(k) = expectedLatency;
            end
        end
    end
end
end