solver_mam_basic.m

function [QN,UN,RN,TN,CN,XN,totiter] = solver_mam_basic(sn, options)
% [Q,U,R,T,C,X] = SOLVER_MAM_BASIC(QN, PH, OPTIONS)
 
% This solver uses MAM to solve queues in isolation, but simplifies the
% traffic equations by using visits to rescale the flows into the queue
% inputs.
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
config = options.config;
tol = options.tol;
 
PH = sn.proc;
%% generate local state spaces
I = sn.nnodes;
M = sn.nstations;
K = sn.nclasses;
C = sn.nchains;
N = sn.njobs';
V = cellsum(sn.visits);
S = 1./sn.rates;
Strue = S; % service times as declared, used to report utilization below
Lchain = sn_get_demands_chain(sn);
 
% Interference of self-looping classes on the classes they share a station
% with. A self-looping class holds njobs customers at its reference station in
% perpetuity, so a job of any other class arriving there always finds those
% njobs customers resident and, under any work-conserving discipline, must
% share the server with them: its residence time is multiplied by
% (1 + njobs_slc). Modelling this as an inflated service time keeps the SLC out
% of the decomposition (it has no surrogate arrival stream; see below) while
% still charging the other classes for the capacity it holds. On the sanity SLC
% models this reproduces the CTMC solution exactly. Delay stations are exempt:
% with infinite servers there is no contention to share.
slcjobs = zeros(M,1);
for k=1:K
    if sn.isslc(k)
        ist_k = sn.refstat(k);
        if isfinite(sn.nservers(ist_k))
            slcjobs(ist_k) = slcjobs(ist_k) + sn.njobs(k);
        end
    end
end
for ist=1:M
    if slcjobs(ist) > 0
        for k=1:K
            if ~sn.isslc(k)
                S(ist,k) = S(ist,k)*(1+slcjobs(ist));
            end
        end
        % The chain demands drive the throughput fixed point and must be
        % consistent with the inflated service times.
        Lchain(ist,:) = Lchain(ist,:)*(1+slcjobs(ist));
    end
end
 
QN = zeros(M,K);
UN = zeros(M,K);
RN = zeros(M,K);
TN = zeros(M,K);
CN = zeros(1,K);
XN = zeros(1,K);
 
% Track stations using exact MAP/D/c solver (skip post-processing for these)
mapdcStations = false(M, 1);
 
pie = {};
D0 = {};
 
lambda = zeros(1,C);
chainSysArrivals = cell(1,C);
TN_1 = TN+Inf;
 
it = 0;
 
% open queueing system (one node is the external world)
% first build the joint arrival process
for ist=1:M
    switch sn.sched(ist)
        case {SchedStrategy.FCFS, SchedStrategy.HOL, SchedStrategy.FCFSPRPRIO, SchedStrategy.PS}
            for k=1:K
                % Skip Det processes - they will be handled by qsys_mapdc
                if sn.procid(ist, k) == ProcessType.DET
                    % For Det, we don't need MAP representation since we use exact solver
                    pie{ist}{k} = [];
                    D0{ist,k} = [];
                    continue;
                end
                % divide service time by number of servers and put
                % later a surrogate delay server in tandem to compensate.
                % map_scale normalizes through map_normalize, which clips
                % negative entries to zero; on a matrix-exponential or
                % rational representation that is not a repair of numerical
                % noise but a replacement of the process by a different,
                % Markovian one, so those are rescaled by the rate alone.
                if sn.procid(ist, k) == ProcessType.ME || sn.procid(ist, k) == ProcessType.RAP
                    ratio = map_mean(PH{ist}{k}) / (S(ist,k)/sn.nservers(ist));
                    PH{ist}{k} = {PH{ist}{k}{1}*ratio, PH{ist}{k}{2}*ratio};
                else
                    PH{ist}{k} = map_scale(PH{ist}{k}, S(ist,k)/sn.nservers(ist));
                end
                pie{ist}{k} = map_pie(PH{ist}{k});
                D0{ist,k} = PH{ist}{k}{1};
                if any(isnan(D0{ist,k}))
                    D0{ist,k} = -GlobalConstants.Immediate;
                    pie{ist}{k} = 1;
                    PH{ist}{k} = map_exponential(GlobalConstants.Immediate);
                end
            end
    end
end % i
 
isOpen = false;
isClosed = false;
if any(isinf(sn.njobs))
    isOpen = true;
end
if any(isfinite(sn.njobs))
    isClosed = true;
end
isMixed = isOpen & isClosed;
 
% Self-looping classes never leave their reference station, so a chain made
% up of them admits no surrogate arrival stream: the flow-equivalent rate
% that would reproduce a permanently resident customer does not exist, and
% the regula falsi below would drive it upwards without bound. Such chains
% are excluded from the fixed point and pinned by the SLC clamp at the end
% of this function, mirroring solver_mna_closed and
% solver_mam_basic_mmap_closed.
isslcchain = false(1,C);
for c=1:C
    isslcchain(c) = all(sn.isslc(sn.inchain{c}));
end
 
lambdas_inchain = cell(1,C);
for c=1:C
    inchain = sn.inchain{c};
    lambdas_inchain{c} = sn.rates(sn.refstat(inchain(1)),inchain);
    %lambdas_inchain{c} = lambdas_inchain{c}(isfinite(lambdas_inchain{c}));
    lambda(c) = sum(lambdas_inchain{c}(isfinite(lambdas_inchain{c})));
    if isslcchain(c)
        lambda(c) = 0;
        lambdas_inchain{c} = zeros(1,length(inchain));
    end
    ist = sn.refstat(inchain(1)); % identical for all classes in the chain
    if isinf(sum(sn.njobs(inchain))) % if open chain
        % ist here is the source
        % assemble a MMAP for the arrival process from all classes
        for k=1:K
            if isnan(PH{ist}{k}{1})
                PH{ist}{k} = map_exponential(Inf); % no arrivals from this class
            end
        end
        inchain = sn.inchain{c};
        k = inchain(1);
        chainSysArrivals{c} = {PH{ist}{k}{1},PH{ist}{k}{2},PH{ist}{k}{2}};
        for ki=2:length(inchain)
            k = inchain(ki);
            if isnan(PH{ist}{k}{1})
                PH{ist}{k} = map_exponential(Inf); % no arrivals from this class
            end
            chainSysArrivals{c} = mmap_super_safe({chainSysArrivals{c},{PH{ist}{k}{1},PH{ist}{k}{2},PH{ist}{k}{2}}}, config.space_max, 'default');
        end
        TN(ist,inchain') = lambdas_inchain{c};
    end
end
 
sd = isfinite(sn.nservers);
 
isclosedchain = false(1,C);
isopenchain = false(1,C);
for c=1:C
    isopenchain(c) = isinf(sum(sn.njobs(sn.inchain{c})));
    isclosedchain(c) = ~isopenchain(c) && ~isslcchain(c);
end
% A mixed network is the only case in which the uniform 1/Umax backoff below is
% unsafe, because there it also scales the open chains, whose arrival rate is
% exogenous data fixed by their source rather than an unknown of the
% decomposition. Purely open and purely closed networks keep the original
% iteration untouched: with a single chain kind the uniform backoff is exactly
% the closed-only backoff, and a closed chain is expected to drive its
% bottleneck to saturation as its population grows.
ismixed = any(isclosedchain) && any(isopenchain);
% In a mixed network the closed chains may claim no more than the capacity the
% open traffic leaves free, and the iterate must stay strictly inside the
% stability region: the queue-length formulas below are of M/G/1 type and
% diverge at total utilization 1.
Ulim = 1 - GlobalConstants.CoarseTol;
 
% Stations already flagged as carrying a non-phase-type service process, scoped
% to this solve so the fixed-point iteration reports each station once while
% every new model is still warned about.
meWarned = false(1, M);
 
while max(max(abs(TN-TN_1))) > tol && it <= options.iter_max %#ok<max>
    it = it + 1;
    TN_1 = TN;
    Umax = max(sum(UN(sd,:),2));
    if ismixed || Umax < 1
        for c=1:C
            inchain = sn.inchain{c};
            if isclosedchain(c)
                Nc = sum(sn.njobs(inchain)); % closed population
                QNc = max(tol, sum(sum(QN(:,inchain),2,"omitnan"))); %#ok<NANSUM>
                TNlb = Nc./sum(Lchain(:,c));
                if it == 1
                    lambda(c) = TNlb; % lower bound
                else
                    lambda(c) = lambda(c) * it/options.iter_max + (Nc / QNc)  * lambda(c) * (options.iter_max-it)/options.iter_max; % iteration-averaged regula falsi;
                end
            end
        end
    end
    if ismixed
        % Back the closed chains alone onto the residual capacity of the
        % busiest finite-server station. Scaling every chain by 1/Umax instead
        % drags the open chains down too, so the open classes report a
        % throughput below their own source rate, and it parks the iterate
        % exactly on U==1, where the queue lengths diverge and, since the
        % backoff then becomes a no-op, the iteration stops moving.
        % Summed rather than multiplied out: with no open chain the empty
        % product would collapse to 1x0 instead of a zero column.
        Uchain  = Lchain(sd,:) .* lambda;
        Uopen   = sum(Uchain(:,~isclosedchain), 2);
        Uclosed = sum(Uchain(:, isclosedchain), 2);
        binding = Uclosed > tol;
        if any(binding)
            theta = min((Ulim - Uopen(binding)) ./ Uclosed(binding));
            if theta < 1
                lambda(isclosedchain) = lambda(isclosedchain) * max(0, theta);
            end
        end
    elseif Umax >= 1
        lambda = lambda * 1/Umax;
    end
 
    for c=1:C
        inchain = sn.inchain{c};
        if isslcchain(c)
            % Vanishing-rate surrogate, as in solver_mam_basic_mmap_inner: a
            % self-looping chain must not inject arrivals into the shared-queue
            % decomposition. An exactly zero rate would give a singular MMAP.
            chainSysArrivals{c} = mmap_exponential(GlobalConstants.Zero);
        elseif ~isinf(sum(sn.njobs(inchain)))
            % Closed chain: Poisson surrogate at the current throughput
            % iterate. Open chains keep the source arrival MMAP assembled
            % before the loop (its rate is re-imposed per station by
            % mmap_scale), so the external arrival variability is not
            % silently discarded.
            chainSysArrivals{c} = mmap_exponential(lambda(c));
        end
        for m=1:M
            TN(m,inchain) = V(m,inchain) .* lambda(c);
        end
    end
 
    for ind=1:I
        if sn.isstation(ind)
            ist = sn.nodeToStation(ind);
            switch sn.nodetype(ind)
                case NodeType.Join
                    for c=1:C
                        inchain = sn.inchain{c};
                        for k=inchain
                            fanin = nnz(sn.rtnodes(:, (ind-1)*K+k));
                            TN(ist,k) = lambda(c)*V(ist,k)/fanin;
                            UN(ist,k) = 0;
                            QN(ist,k) = 0;
                            RN(ist,k) = 0;
                        end
                    end
                otherwise
                    switch sn.sched(ist)
                        case SchedStrategy.INF
                            for c=1:C
                                inchain = sn.inchain{c};
                                for k=inchain
                                    if V(ist,k) == 0
                                        % Class k does not visit this station, so
                                        % its service rate is undefined (NaN) and
                                        % all of its metrics here are null. Guard
                                        % explicitly: S(ist,k)*TN(ist,k) would be
                                        % NaN*0 = NaN, which then propagates into
                                        % the Umax saturation test below (NaN>=1 is
                                        % false, so the guard would be bypassed
                                        % silently). Mirrors the V(ist,r)>0 test in
                                        % solver_mam_basic_mmap_inner.
                                        TN(ist,k) = 0;
                                        UN(ist,k) = 0;
                                        QN(ist,k) = 0;
                                        RN(ist,k) = 0;
                                        continue;
                                    end
                                    TN(ist,k) = lambda(c)*V(ist,k);
                                    % Utilization Law: a c-server station holds TN*S/c of
                                    % its capacity. An INF station has nservers = Inf, so
                                    % dividing would annihilate it (S*TN/Inf = 0); an
                                    % infinite server reports U = QLen = TN*S instead, the
                                    % convention every other solver follows.
                                    UN(ist,k) = S(ist,k)*TN(ist,k);
                                    % TN already carries the visits, so Little's law reads
                                    % QN = TN*S; a second V factor undercounts the delay
                                    % jobs whenever the chain switches class (V<1), which
                                    % skews the closed-chain regula falsi target QNc and,
                                    % through the post-loop rescale, deflates the response
                                    % time of the queueing stations in the same chain. The
                                    % PS and FCFS branches already use the single-V form.
                                    QN(ist,k) = TN(ist,k).*S(ist,k);
                                    RN(ist,k) = QN(ist,k)/TN(ist,k);
                                end
                            end
                        case SchedStrategy.PS
                            for c=1:C
                                inchain = sn.inchain{c};
                                for k=inchain
                                    if V(ist,k) == 0
                                        % See the INF branch: a non-visiting class
                                        % has an undefined rate, and NaN*0 = NaN
                                        % would poison the Uden sum below.
                                        TN(ist,k) = 0;
                                        UN(ist,k) = 0;
                                        continue;
                                    end
                                    TN(ist,k) = lambda(c)*V(ist,k);
                                    % Utilization Law: a c-server station holds TN*S/c of its capacity.
                                    UN(ist,k) = S(ist,k)*TN(ist,k)/sn.nservers(ist);
                                end
                                %Nc = sum(sn.njobs(inchain)); % closed population
                                Uden = min([1-GlobalConstants.FineTol,sum(UN(ist,:))]);
                                for k=inchain
                                    if V(ist,k) == 0
                                        QN(ist,k) = 0;
                                        RN(ist,k) = 0;
                                        continue;
                                    end
                                    %QN(ist,k) = (UN(ist,k)-UN(ist,k)^(Nc+1))/(1-Uden); % geometric bound type approximation
                                    QN(ist,k) = UN(ist,k)/(1-Uden);
                                    RN(ist,k) = QN(ist,k)/TN(ist,k);
                                end
                            end
                        case {SchedStrategy.FCFS, SchedStrategy.HOL, SchedStrategy.FCFSPRPRIO}
                            chainArrivalAtNode = cell(1,C);
                            rates = cell(M,C);
                            mapdcUsed = false;  % Flag for exact MAP/D/c solver
                            finiteCapUsed = false;  % Flag for finite-capacity truncate-and-renormalize
                            finiteCapLossPerClass = [];  % per-class loss, exact MMAP/G/1/K branch only
                            for c=1:C %for each chain
                                rates{ist,c} = V(ist,:) .* lambda(c); % visits of classes within the chain
                                inchain = find(sn.chains(c,:))';
                                markProb = rates{ist,c}(inchain) / sum(rates{ist,c}(inchain));
                                markProb(isnan(markProb)) = 0;
                                chainArrivalAtNode{c} = mmap_mark(chainSysArrivals{c}, markProb);
                                % mmap_normalize clips negative entries to zero and
                                % re-derives the diagonal. On a genuine MAP that
                                % repairs roundoff; on a matrix-exponential or
                                % rational arrival it substitutes a different,
                                % Markovian process, which changes the
                                % autocorrelation that the downstream MMAPPH1FCFS
                                % is there to account for. Same reasoning as the
                                % ME/RAP branch of the service rescaling above.
                                if mam_chain_arrival_is_markovian(sn, c)
                                    chainArrivalAtNode{c} = mmap_normalize(chainArrivalAtNode{c});
                                end
                                % Target the per-class mean inter-arrival times of
                                % THIS chain's classes. mmap_scale reads M(1:C)
                                % where C is the number of marks in the MMAP, i.e.
                                % length(inchain) here, so passing the full 1:K
                                % vector silently targets classes 1..length(inchain)
                                % instead of the chain's own classes. The two
                                % coincide only when the chain occupies the leading
                                % class indices; otherwise every chain but the first
                                % is scaled to another chain's visit rate.
                                % A chain that offers no arrivals here (a self-looping
                                % chain, whose lambda is pinned to zero) has an
                                % infinite target mean: mmap_scale would then scale the
                                % MMAP by mean/Inf = 0 and return an all-zero, singular
                                % MAP. It is already a vanishing-rate arrival process,
                                % so leave it untouched.
                                tgtrates = rates{ist,c}(inchain);
                                if any(tgtrates > 0)
                                    chainArrivalAtNode{c} = mmap_scale(chainArrivalAtNode{c}, 1./tgtrates, 0); % non-iterative approximation
                                end
                                chainIsMarkovian = mam_chain_arrival_is_markovian(sn, c);
                                if c == 1
                                    if chainIsMarkovian
                                        aggrArrivalAtNode = mmap_super_safe({chainArrivalAtNode{c}, mmap_exponential(0,1)}, config.space_max, 'default');
                                        aggrArrivalAtNode = {aggrArrivalAtNode{1} aggrArrivalAtNode{2} aggrArrivalAtNode{2}};
                                    else
                                        % Superposing with a zero-rate exponential is the
                                        % identity on (D0,D1) -- the Kronecker sum with an
                                        % order-1 null adds nothing -- and is used here only
                                        % to reach the {D0,D1,D1} marked shape. What it also
                                        % does is run mmap_normalize inside mmap_super, which
                                        % clips the negative off-diagonal entries a
                                        % matrix-exponential or rational arrival legitimately
                                        % carries and substitutes a Markovian process with a
                                        % different autocorrelation. Build the shape directly
                                        % instead, so the arrival reaching MMAPPH1FCFS is the
                                        % process the user declared.
                                        aggrArrivalAtNode = {chainArrivalAtNode{c}{1}, chainArrivalAtNode{c}{2}, chainArrivalAtNode{c}{2}};
                                    end
                                    lc = map_lambda(chainArrivalAtNode{c});
                                    if lc>0
                                        aggrArrivalAtNode = mmap_scale(aggrArrivalAtNode, 1/lc, 0); % non-iterative approximation
                                    end
                                else
                                    if ~chainIsMarkovian
                                        % mmap_super_safe normalizes, so the ME/RAP component
                                        % of a multi-chain superposition is replaced by a
                                        % Markovian one and its correlation is only
                                        % approximated. Say so rather than report the result
                                        % as if it carried the declared process.
                                        line_warning_always(mfilename, ...
                                            'Chain %d has a matrix-exponential or rational arrival process, which the superposition of several chains at station %s cannot represent exactly. Its autocorrelation is approximated by the closest Markovian arrival process.', ...
                                            c, sn.nodenames{sn.stationToNode(ist)});
                                    end
                                    aggrArrivalAtNode = mmap_super_safe({aggrArrivalAtNode, chainArrivalAtNode{c}}, config.space_max, 'default');
                                end
                            end
                            Qret = cell(1,K);
                            if (sn.sched(ist)==SchedStrategy.HOL && any(sn.classprio ~= sn.classprio(1))) % if priorities are not identical; sched holds numeric ids so == must be used (strcmp on numerics is always false)
                                [uK,iK] = unique(sn.classprio);
                                % BUTools convention: D1=lowest priority, DK=highest priority
                                % LINE convention: lower value = higher priority
                                % unique() returns ascending order, so we need to reverse for BUTools
                                iK = flipud(iK(:));
                                if length(uK) == length(sn.classprio) % if all priorities are different
                                    [Qret{iK'}] = MMAPPH1NPPR({aggrArrivalAtNode{[1;2+iK]}}, {pie{ist}{iK}}, {D0{ist,iK}}, 'ncMoms', 1);
                                else
                                    line_error(mfilename,'Solver MAM requires either identical priorities or all distinct priorities');
                                end
                            elseif (sn.sched(ist)==SchedStrategy.FCFSPRPRIO && any(sn.classprio ~= sn.classprio(1))) % FCFS preemptive resume priority
                                [uK,iK] = unique(sn.classprio);
                                % BUTools convention: D1=lowest priority, DK=highest priority
                                % LINE convention: lower value = higher priority
                                % unique() returns ascending order, so we need to reverse for BUTools
                                iK = flipud(iK(:));
                                if length(uK) == length(sn.classprio) % if all priorities are different
                                    [Qret{iK'}] = MMAPPH1PRPR({aggrArrivalAtNode{[1;2+iK]}}, {pie{ist}{iK}}, {D0{ist,iK}}, 'ncMoms', 1);
                                else
                                    line_error(mfilename,'Solver MAM requires either identical priorities or all distinct priorities');
                                end
                            else
                                aggrUtil = sum(mmap_lambda(aggrArrivalAtNode)./(GlobalConstants.FineTol+sn.rates(ist,1:K)*sn.nservers(ist)), 'omitnan');   %% to debug
                                aggrLambda = mmap_lambda(aggrArrivalAtNode);
                                if aggrUtil < 1-GlobalConstants.FineTol
                                    if any(isinf(N))
                                        % Check for MAP/D/c: single-class open model with deterministic service
                                        isMapDc = (K == 1) && (sn.procid(ist, 1) == ProcessType.DET);
                                        % Check for D/M/c: single-class open, exp service at this queue,
                                        % deterministic source elsewhere.
                                        isDMc = false;
                                        dmcSourceIdx = -1;
                                        if K == 1 && ~isMapDc && sn.procid(ist, 1) == ProcessType.EXP
                                            for jst = 1:sn.nstations
                                                if jst ~= ist && sn.procid(jst, 1) == ProcessType.DET
                                                    isDMc = true;
                                                    dmcSourceIdx = jst;
                                                    break;
                                                end
                                            end
                                        end
                                        % Check for PH/M/c: single-class open, exp service at this queue,
                                        % source has PH inter-arrivals (Erlang/Coxian/MAP/PH/...).
                                        % Covers both c=1 (PH/M/1) and c>1 (PH/M/c via matrix-geometric).
                                        isPhM1 = false;
                                        phM1SourceIdx = -1;
                                        if K == 1 && ~isMapDc && ~isDMc && sn.procid(ist, 1) == ProcessType.EXP ...
                                                && isfinite(sn.nservers(ist)) && sn.nservers(ist) >= 1
                                            for jst = 1:sn.nstations
                                                if jst == ist
                                                    continue;
                                                end
                                                srcProc = sn.procid(jst, 1);
                                                if srcProc ~= ProcessType.EXP && srcProc ~= ProcessType.DET ...
                                                        && srcProc ~= ProcessType.IMMEDIATE && srcProc ~= ProcessType.DISABLED
                                                    % qsys_phmc reads the arrival only through
                                                    % its (pie, D0) phase-type pair, i.e. its
                                                    % renewal marginal. That is exact for a
                                                    % renewal arrival and discards the
                                                    % autocorrelation of any other, so the gate
                                                    % is RENEWAL, not merely non-exponential.
                                                    % Being non-exponential was the old test and
                                                    % it let a correlated MAP, RAP or ME arrival
                                                    % onto the closed form: the answer then came
                                                    % out near the Kingman value implied by the
                                                    % arrival SCV alone. An ME arrival is renewal
                                                    % and stays on this fast path.
                                                    if ~mam_srcproc_is_renewal(sn, jst)
                                                        continue;
                                                    end
                                                    isPhM1 = true;
                                                    phM1SourceIdx = jst;
                                                    break;
                                                elseif srcProc == ProcessType.EXP && sn.nservers(ist) > 1 ...
                                                        && sn.nstations == 2 ...
                                                        && sn.nodetype(sn.stationToNode(jst)) == NodeType.Source
                                                    % M/M/c: the single-fast-server surrogate used by the
                                                    % generic path is inexact for c>1; dispatch to the exact
                                                    % PH/M/c matrix-geometric solver (Erlang-C result)
                                                    isPhM1 = true;
                                                    phM1SourceIdx = jst;
                                                    break;
                                                end
                                            end
                                        end
                                        % Check for finite-buffer FCFS. qsys_phmc, qsys_dmc
                                        % and qsys_mapdc are all infinite-buffer closed
                                        % forms, so letting them win here would silently
                                        % discard sn.cap and report a lossless queue. A
                                        % finite buffer therefore takes precedence over
                                        % them rather than being disabled by them.
                                        isFiniteCap = isfinite(sn.cap(ist));
                                        if isFiniteCap
                                            isPhM1 = false;
                                            isDMc = false;
                                            isMapDc = false;
                                        end
                                        if isPhM1
                                            muQ = 1.0 / S(ist, 1);
                                            try
                                                phPair = sn.proc{phM1SourceIdx}{1};
                                                D0_ph = phPair{1};
                                                D1_ph = phPair{2};
                                                pie_src = map_pie({D0_ph, D1_ph});
                                                phm1Result = qsys_phmc(pie_src, D0_ph, muQ, sn.nservers(ist));
                                                Qret{1} = phm1Result.meanQueueLength;
                                                mapdcUsed = true;
                                                mapdcStations(ist) = true;
                                                line_debug(options, 'Using exact PH/M/%d solver: Q=%.4f, W=%.4f', ...
                                                    sn.nservers(ist), phm1Result.meanQueueLength, phm1Result.meanWaitingTime);
                                            catch
                                                isPhM1 = false;
                                            end
                                        end
                                        if isDMc
                                            muQ = 1.0 / S(ist, 1);
                                            lamD = sn.rates(dmcSourceIdx, 1);
                                            try
                                                dmcResult = qsys_dmc(lamD, muQ, sn.nservers(ist));
                                                Qret{1} = dmcResult.meanQueueLength;
                                                mapdcUsed = true;  % skip surrogate delay correction
                                                mapdcStations(ist) = true;
                                                line_debug(options, 'Using exact D/M/%d solver: Q=%.4f, W=%.4f', ...
                                                    sn.nservers(ist), dmcResult.meanQueueLength, dmcResult.meanWaitingTime);
                                            catch
                                                isDMc = false;
                                            end
                                        end
                                        if isPhM1 || isDMc
                                            % handled above; skip the MAP/D/c, finite-cap, and MMAPPH1FCFS branches
                                        elseif isMapDc
                                            % Use exact MAP/D/c solver from Q-MAM
                                            D0_arr = aggrArrivalAtNode{1};
                                            D1_arr = aggrArrivalAtNode{2};
                                            detServiceTime = S(ist, 1);  % 1/rate = service time
                                            numServers = sn.nservers(ist);
 
                                            mapdcResult = qsys_mapdc(D0_arr, D1_arr, detServiceTime, numServers);
                                            Qret{1} = mapdcResult.meanQueueLength;
                                            % Store result for later use to skip surrogate delay adjustment
                                            mapdcUsed = true;
                                            mapdcStations(ist) = true;  % Mark for skipping post-processing
                                            line_debug('Using exact MAP/D/%d solver: Q=%.4f, W=%.4f', ...
                                                numServers, mapdcResult.meanQueueLength, mapdcResult.meanWaitingTime);
                                        elseif isFiniteCap
                                            % Finite-buffer FCFS. Detect exact M/M/c/K
                                            % (Poisson arrivals + same exp service across
                                            % classes); else use MMAP[K]/PH[K]/1/FCFS
                                            % truncate-and-renormalize approximation.
                                            mapdcUsed = false;
                                            capK = sn.cap(ist);
                                            finiteCapLossPerClass = [];
                                            [isMmck, muMmck] = mam_detect_mmck(sn, ist, K, aggrArrivalAtNode);
                                            if isMmck
                                                aggrLambdaTotal = sum(aggrLambda, 'omitnan');
                                                exactRes = qsys_mmck(aggrLambdaTotal, muMmck, sn.nservers(ist), capK);
                                                finiteCapMeanQ = exactRes.meanQueueLength;
                                                finiteCapLossProb = exactRes.lossProbability;
                                                finiteCapP0 = exactRes.queueLengthDist(1);
                                                line_debug('Using exact M/M/%d/%d: Q=%.4f, ploss=%.4f', ...
                                                    sn.nservers(ist), capK, finiteCapMeanQ, finiteCapLossProb);
                                            elseif sn.nservers(ist) == 1
                                                % Exact MMAP[K]/G/1/K. The embedded chain at
                                                % departure epochs resolves the buffer level
                                                % jointly with the arrival phase, so class k,
                                                % which arrives out of phase i at rate
                                                % (D1k*e)_i, gets its own loss ratio. The
                                                % truncate-renorm fallback below can only
                                                % return one blocking probability for every
                                                % class, making loss rate-proportional and
                                                % blind to interarrival variability.
                                                svcMix = mam_svc_mixture({aggrArrivalAtNode{[1,3:end]}}, ...
                                                    {pie{ist}{:}}, {D0{ist,:}});
                                                exRes = qsys_mmapg1k(aggrArrivalAtNode{1}, ...
                                                    aggrArrivalAtNode(3:end), svcMix, capK);
                                                finiteCapMeanQ = exRes.meanQueueLength;
                                                finiteCapLossProb = exRes.lossAggregate;
                                                finiteCapLossPerClass = exRes.lossRatio;
                                                finiteCapP0 = exRes.p0;
                                                line_debug('Using exact MMAP/G/1/%d: Q=%.4f, ploss=%.4f, per-class ploss=%s', ...
                                                    capK, finiteCapMeanQ, finiteCapLossProb, mat2str(finiteCapLossPerClass, 4));
                                            else
                                                [meanQ_fc, lossProb_fc, p_norm_fc] = mam_truncate_renorm( ...
                                                    {aggrArrivalAtNode{[1,3:end]}}, {pie{ist}{:}}, {D0{ist,:}}, capK);
                                                finiteCapMeanQ = meanQ_fc;
                                                finiteCapLossProb = lossProb_fc;
                                                finiteCapP0 = p_norm_fc(1);
                                                line_debug('Using MMAP/PH/%d/%d truncate-renorm: Q=%.4f, ploss=%.4f', ...
                                                    sn.nservers(ist), capK, meanQ_fc, lossProb_fc);
                                            end
                                            finiteCapUsed = true;
                                            for k=1:K
                                                Qret{k} = NaN;  % filled per class in result loop
                                            end
                                            mapdcStations(ist) = true;  % skip surrogate delay 2nd pass
                                        elseif sn.isfunction(ist)
                                            % Open setup/delay-off station. qbd_setupdelayoff
                                            % races the delay-off against the aggregate
                                            % interarrival Exp(lambda), which is exactly the
                                            % idle seen by a Poisson stream, so it is the right
                                            % surrogate here (unlike the closed case below,
                                            % where the idle is per instance). Without this the
                                            % station fell through to the plain MMAPPH1FCFS
                                            % below and the setup was silently ignored: an
                                            % M/M/1 with setup returned the setup-free
                                            % rho/(1-rho).
                                            mapdcUsed = false;
                                            % nodeparam is NODE-indexed (refreshLocalVars), while
                                            % ist is a STATION index; the two coincide only when
                                            % no non-station node precedes the queue. Reading
                                            % nodeparam{ist} directly errored out on any model
                                            % with e.g. a Router or ClassSwitch ahead of the
                                            % setup station.
                                            alpharate = map_lambda(sn.nodeparam{sn.stationToNode(ist)}{end}.setupTime);
                                            alphascv = map_scv(sn.nodeparam{sn.stationToNode(ist)}{end}.setupTime);
                                            betarate = map_lambda(sn.nodeparam{sn.stationToNode(ist)}{end}.delayoffTime);
                                            betascv = map_scv(sn.nodeparam{sn.stationToNode(ist)}{end}.delayoffTime);
                                            mu_k = zeros(1, K);
                                            lambda_k = zeros(1, K);
                                            active_k = false(1, K);
                                            for k=1:K
                                                pie_k = pie{ist}{k};
                                                if ~isempty(pie_k) && ~isnan(pie_k(1))
                                                    mu_k(k) = 1 / (pie_k * inv(-D0{ist,k}) * ones(size(pie_k))');
                                                    c = find(sn.chains(:,k), 1);
                                                    lambda_k(k) = rates{ist,c}(k);
                                                    active_k(k) = true;
                                                end
                                            end
                                            rho_k = lambda_k ./ mu_k;
                                            rho_k(~active_k) = 0;
                                            rho_total = sum(rho_k);
                                            if rho_total > 0
                                                % Aggregate service rate consistent with the
                                                % aggregate load, then split the aggregate queue
                                                % across classes in proportion to their load.
                                                % aggrLambda is per class; the QBD takes the
                                                % scalar total arrival rate seen by the server.
                                                aggrLambdaTotal = sum(aggrLambda, 'omitnan');
                                                aggrRate = aggrLambdaTotal / rho_total;
                                                Q_total = qbd_setupdelayoff(aggrLambdaTotal, aggrRate, alpharate, alphascv, betarate, betascv);
                                                for k=1:K
                                                    if active_k(k)
                                                        Qret{k} = Q_total * rho_k(k) / rho_total;
                                                    else
                                                        Qret{k} = 0;
                                                    end
                                                end
                                            else
                                                for k=1:K
                                                    Qret{k} = 0;
                                                end
                                            end
                                        else
                                            mapdcUsed = false;
                                            % MMAPPH1FCFS treats service as a renewal
                                            % phase-type (marginal), discarding service
                                            % autocorrelation. For the single-class
                                            % single-server M/MAP/1 case with a genuinely
                                            % correlated service, use the exact MAP/MAP/1
                                            % queue (Q_CT_MAP_MAP_1), which carries the
                                            % service phase across departures.
                                            % A matrix-exponential or rational service
                                            % process is not phase-type, so the (pie,D0)
                                            % pair MMAPPH1FCFS reads does not describe it
                                            % and the PH path would silently return the
                                            % wrong numbers. The RAP/RAP/1 QBD applies
                                            % whenever the SERVICE is ME or RAP, whatever
                                            % the arrival type: the aggregate arrival
                                            % stream is handed over as a (D0,D1) pair, and
                                            % a Poisson stream, a MAP and a PH renewal
                                            % stream are all RAPs.
                                            isMEorRAPsvc = any(sn.procid(ist,:) == ProcessType.ME | ...
                                                sn.procid(ist,:) == ProcessType.RAP);
                                            useRapRap1 = false;
                                            if isMEorRAPsvc
                                                if (K == 1) && (sn.nservers(ist) == 1)
                                                    useRapRap1 = true;
                                                else
                                                    % qbd_raprap1 is a single-class
                                                    % single-server QBD. Say so rather than
                                                    % let an ME model be answered by the
                                                    % phase-type approximation without the
                                                    % user knowing.
                                                    % line_warning_always, not
                                                    % line_warning: this reports that
                                                    % the returned numbers are not
                                                    % exact, so it must be raised for
                                                    % every model it applies to rather
                                                    % than be hidden as a repeat of the
                                                    % previous model's warning.
                                                    if ~meWarned(ist)
                                                        meWarned(ist) = true;
                                                        line_warning_always(mfilename, ...
                                                            'Station %s has a matrix-exponential or rational service process, which the RAP/RAP/1 analysis supports only with a single class at a single server (here %d classes, %g servers). Falling back to the phase-type approximation MMAPPH1FCFS, which is not exact for this service process.', ...
                                                            sn.nodenames{sn.stationToNode(ist)}, K, sn.nservers(ist));
                                                    end
                                                end
                                            end
                                            useMapMap1 = ~useRapRap1 && (K == 1) && (sn.nservers(ist) == 1) && ...
                                                (abs(map_acf(PH{ist}{1}, 1)) > GlobalConstants.CoarseTol);
                                            if useRapRap1
                                                % qbd_raprap1 carries no sojourn-time
                                                % information, so only the queue-length
                                                % mean is taken here; the response time is
                                                % recovered downstream by Little's law.
                                                [~, QNrap] = qbd_raprap1({aggrArrivalAtNode{1}, aggrArrivalAtNode{3}}, ...
                                                    {PH{ist}{1}{1}, PH{ist}{1}{2}});
                                                Qret{1} = QNrap;
                                            elseif useMapMap1
                                                ql = Q_CT_MAP_MAP_1(aggrArrivalAtNode{1}, aggrArrivalAtNode{3}, ...
                                                    PH{ist}{1}{1}, PH{ist}{1}{2}, 'MaxNumComp', 100000);
                                                ql = ql(:);
                                                Qret{1} = sum((0:numel(ql)-1)' .* ql);
                                            else
                                                [Qret{1:K}] = MMAPPH1FCFS({aggrArrivalAtNode{[1,3:end]}}, {pie{ist}{:}}, {D0{ist,:}}, 'ncMoms', 1);
                                            end
                                        end
                                    else % all closed classes
                                        maxLevel = sum(N(isfinite(N)))+1;
                                        D = {aggrArrivalAtNode{[1,3:end]}};
                                        pdistr_k = cell(1,K);
                                        if map_lambda(D)< GlobalConstants.FineTol
                                            for k=1:K
                                                pdistr_k = [1-GlobalConstants.FineTol, GlobalConstants.FineTol];
                                                Qret{k} = GlobalConstants.FineTol / sn.rates(ist);
                                            end
                                        else
                                            if sn.isfunction(ist)
                                                alpharate = map_lambda(sn.nodeparam{sn.stationToNode(ist)}{end}.setupTime);
                                                betarate = map_lambda(sn.nodeparam{sn.stationToNode(ist)}{end}.delayoffTime);
                                                betascv = map_scv(sn.nodeparam{sn.stationToNode(ist)}{end}.delayoffTime);
                                                lambda_k = zeros(1, K);
                                                active_k = false(1, K);
                                                for k=1:K
                                                    pie_k = pie{ist}{k};
                                                    if ~isnan(pie_k(1))
                                                        c = find(sn.chains(:,k), 1);
                                                        lambda_k(k) = rates{ist,c}(k);
                                                        active_k(k) = true;
                                                    end
                                                end
                                                if any(active_k)
                                                    % Per-instance cold-start race. In a closed layer the
                                                    % function instances are the chain's own tokens: each
                                                    % token's idle period between two consecutive requests
                                                    % is the client think time ZT of its chain (the
                                                    % exponential surrogate at the infinite-server client
                                                    % delay), so an arriving request pays the setup iff the
                                                    % delay-off timer expires within Exp(1/ZT):
                                                    %   p_cold = P(delayoff < Exp(1/ZT)) = LST_delayoff(1/ZT).
                                                    % The response time then carries the full setup latency,
                                                    %   R_k = p_cold*E[setup] + S_k,
                                                    % and QN follows by Little's law, which keeps the
                                                    % regula falsi above linear in lambda. The open-arrival
                                                    % surrogate qbd_setupdelayoff(aggrLambda,...) is not
                                                    % used here: it races the delay-off against the
                                                    % aggregate interarrival Exp(lambda) instead of the
                                                    % per-instance idle, and its Poisson clumping creates
                                                    % queueing that cannot occur among closed tokens.
                                                    setupMean = 1/alpharate;
                                                    if betascv == 1.0
                                                        % rate taken directly, as in qbd_setupdelayoff
                                                        betaProc = {-betarate, betarate};
                                                    else
                                                        betaProc = APH.fitMeanAndSCV(1/betarate, betascv).getProcess;
                                                    end
                                                    Tb = betaProc{1};
                                                    pie_b = map_pie(betaProc);
                                                    infstat = isinf(sn.nservers);
                                                    for k=1:K
                                                        if active_k(k) && lambda_k(k) > 0 && isfinite(S(ist,k))
                                                            c = find(sn.chains(:,k), 1);
                                                            inchain_k = find(sn.chains(c,:));
                                                            % Lchain and V are per chain cycle; the token
                                                            % visits this station sum(V(ist,inchain)) times
                                                            % per cycle, so the idle between two consecutive
                                                            % requests to an instance is the chain think
                                                            % demand per visit.
                                                            Vtot = sum(V(ist, inchain_k));
                                                            ZT = sum(Lchain(infstat, c), 'omitnan') / max(Vtot, GlobalConstants.FineTol);
                                                            nu = 1 / max(ZT, GlobalConstants.FineTol);
                                                            pcold = pie_b * ((nu*eye(size(Tb,1)) - Tb) \ (-Tb*ones(size(Tb,1),1)));
                                                            % service part uses S/nservers: the surrogate
                                                            % delay correction below adds the remaining
                                                            % S*(nservers-1)/nservers, restoring R_k =
                                                            % p_cold*E[setup] + S_k.
                                                            Qret{k} = lambda_k(k) * (pcold*setupMean + S(ist,k)/sn.nservers(ist));
                                                        elseif active_k(k)
                                                            % A class that offers no load here (zero visits,
                                                            % or no service process: the chain's client-delay
                                                            % token class) holds no jobs; 0*NaN would stamp
                                                            % NaN into QN, void the unguarded chain sum QNc
                                                            % of the second pass, and the NaN-ignoring max()
                                                            % in the response-time fallback would silently
                                                            % report the bare service time for the served
                                                            % classes.
                                                            Qret{k} = 0;
                                                        else
                                                            % An inactive class holds no jobs here. NaN
                                                            % would propagate into the chain sum used by
                                                            % the closed-population renormalization below
                                                            % (QNc = sum(sum(QN(:,inchain)))), turning the
                                                            % whole row NaN and discarding the queue length;
                                                            % the NaN-ignoring max() in the response-time
                                                            % fallback would then silently report the bare
                                                            % service time. The MMAPPH1FCFS path likewise
                                                            % yields 0 for a zero-population class.
                                                            Qret{k} = 0;
                                                        end
                                                    end
                                                else
                                                    for k=1:K
                                                        Qret{k} = NaN;
                                                    end
                                                end
                                            else
                                                % MMAPPH1FCFS returns one queue-length
                                                % distribution per class (see the for k=1:K loop
                                                % around its 'ncDistr' branch), so all K outputs
                                                % must be captured. Requesting a single output
                                                % yields class 1's distribution alone, which was
                                                % then reused for every class: Qret depended on
                                                % the class only through N(k), so two classes
                                                % with equal populations but different visit
                                                % ratios were handed the same queue length, and
                                                % the throughput fixed point resolved that by
                                                % driving every chain to the same lambda.
                                                pdistr_all = cell(1,K);
                                                [pdistr_all{1:K}] = MMAPPH1FCFS(D, {pie{ist}{:}}, {D0{ist,:}}, 'ncDistr', maxLevel);
                                                for k=1:K
                                                    pdistr = pdistr_all{k};
                                                    pdistr_k = abs(pdistr(1:(N(k)+1)));
                                                    % Truncating at the class population folds
                                                    % the tail beyond level N(k) into that level,
                                                    % so the mean below is E[min(Q_k,N(k))]. The
                                                    % complement must therefore be taken over the
                                                    % truncated vector: reading it from the full
                                                    % distribution instead assigns only the mass
                                                    % above N(k), i.e. P(Q>N(k)) rather than
                                                    % P(Q>=N(k)), losing the levels in between.
                                                    pdistr_k(end) = abs(1-sum(pdistr_k(1:end-1)));
                                                    pdistr_k = pdistr_k / sum(pdistr_k(1:(N(k)+1)));
                                                    Qret{k} = max(0,min(N(k),(0:N(k))*pdistr_k(1:(N(k)+1))'));
                                                end
                                            end
                                        end
                                    end
                                else
                                    for k=1:K
                                        Qret{k} = sn.njobs(k);
                                    end
                                end
                            end
                            if finiteCapUsed
                                % Per-class decomposition under FCFS: wait time in queue is
                                % uniform across classes, so per-class response time is
                                % R_k = W_q + S_k. From E[N_total] (truncate-renorm) and
                                % per-class effective rates lambda_k_eff = lambda_k*(1-ploss):
                                %   W_q = E[N_total]/lambda_total_eff - sum(lambda_k_eff*S_k)/lambda_total_eff
                                %   N_k = lambda_k_eff * R_k
                                lambdaInflow = zeros(1, K);
                                for k=1:K
                                    cidx = find(sn.chains(:,k),1);
                                    lambdaInflow(k) = rates{ist,cidx}(k);
                                end
                                lambdaInflow(isnan(lambdaInflow)) = 0;
                                if ~isempty(finiteCapLossPerClass)
                                    % Exact MMAP/G/1/K branch: loss differs by class
                                    TN_eff = lambdaInflow .* (1 - finiteCapLossPerClass(1:K));
                                else
                                    TN_eff = lambdaInflow * (1 - finiteCapLossProb);
                                end
                                sumTN = sum(TN_eff);
                                if sumTN > 0
                                    Savg_eff = sum(TN_eff .* S(ist,1:K), 'omitnan') / sumTN;
                                    Wq = max(0, finiteCapMeanQ / sumTN - Savg_eff);
                                else
                                    Wq = 0;
                                end
                                for k=1:K
                                    TN(ist,k) = TN_eff(k);
                                    UN(ist,k) = TN(ist,k) * S(ist,k) / sn.nservers(ist);
                                    if TN(ist,k) > 0
                                        RN(ist,k) = Wq + S(ist,k);
                                        QN(ist,k) = TN(ist,k) * RN(ist,k);
                                    else
                                        RN(ist,k) = 0;
                                        QN(ist,k) = 0;
                                    end
                                end
                            else
                                QN(ist,:) = cell2mat(Qret);
                                for k=1:K
                                    c = find(sn.chains(:,k),1);
                                    TN(ist,k) = rates{ist,c}(k);
                                    UN(ist,k) = TN(ist,k) * S(ist,k) / sn.nservers(ist);
                                    if sn.isfunction(ist) && ~isfinite(UN(ist,k))
                                        % A class with no service process here (the chain's
                                        % client-delay token class) has S = NaN, and NaN*0
                                        % poisons the Umax saturation gate above, freezing
                                        % the closed-chain regula falsi at the no-contention
                                        % bound TNlb from it=2. At a setup-bearing function
                                        % station the iteration must keep moving so lambda
                                        % reaches QNc = Nc with the cold-start queue length
                                        % included; the class contributes no utilization.
                                        % Elsewhere the NaN freeze is load-bearing (it pins
                                        % LN class-switching host layers) and is preserved.
                                        UN(ist,k) = 0;
                                    end
                                    QN(ist,k) = Qret{k};
                                    if mapdcUsed
                                        % For MAP/D/c, D/M/c, or PH/M/1, use exact results
                                        % (no surrogate delay adjustment).
                                        if isPhM1
                                            RN(ist,k) = phm1Result.meanSojournTime;
                                        elseif isDMc
                                            RN(ist,k) = dmcResult.meanSojournTime;
                                        else
                                            RN(ist,k) = mapdcResult.meanSojournTime;
                                        end
                                    else
                                        % add number of jobs at the surrogate delay server.
                                        % A class with no service process here has S=NaN, and
                                        % the correction is nil for it (TN=0) and for a single
                                        % server ((nservers-1)/nservers=0), but IEEE makes
                                        % NaN*0 = NaN, which would propagate into the chain sum
                                        % QNc below, void the whole row and, via the NaN-ignoring
                                        % max() in the response-time fallback, silently degrade
                                        % the result to the bare service time.
                                        if isfinite(S(ist,k))
                                            QN(ist,k) = QN(ist,k) + TN(ist,k)*S(ist,k) * (sn.nservers(ist)-1)/sn.nservers(ist);
                                        end
                                        RN(ist,k) = QN(ist,k) ./ TN(ist,k);
                                    end
                                end
                            end
                    end
            end
        else % not a station
            switch sn.nodetype(ind)
                case NodeType.Fork
                    %                    line_error(mfilename,'Fork nodes not supported yet by MAM solver.');
            end
        end
    end
    %it
    %max(max(abs(TN-TN_1)))
end
totiter = it + 2;
CN = sum(RN,1);
QN = abs(QN);
for it=1:2 % second pass to rescale again QN based on RN correction
    for c=1:C
        inchain = sn.inchain{c};
        Nc = sum(sn.njobs(inchain));
        if isfinite(Nc)
            QNc = sum(sum(QN(:,inchain)));
            QN(:,inchain) = QN(:,inchain) * (Nc / QNc);
        end
        for ind=1:I
            for k=inchain
                if sn.isstation(ind)
                    ist = sn.nodeToStation(ind);
                    % Skip stations using exact MAP/D/c solver (already have exact values)
                    if mapdcStations(ist)
                        continue;
                    end
                    if V(ist,k)>0
                        if isinf(sn.nservers(ist))
                            RN(ist,k) = S(ist,k);
                        else
                            RN(ist,k) = max([S(ist,k), QN(ist,k) ./ TN(ist,k)]);
                        end
                    else
                        RN(ist,k) = 0;
                    end
                    QN(ist,k) = RN(ist,k) .* TN(ist,k);
                end
            end
        end
        Nc = sum(sn.njobs(inchain)); % closed population
        if Nc == 0 % if closed chain
            QN(:,c)=0;
            UN(:,c)=0;
            RN(:,c)=0;
            TN(:,c)=0;
            CN(c)=0;
            XN(c)=0;
        end
    end
end
 
% SLC clamp: a self-looping class holds its njobs customers at its reference
% station in perpetuity, so its metrics follow in closed form and supersede
% whatever the decomposition produced. Applied last, after the queue-length
% rescaling above, so that it is not undone.
%
% At a queueing station the class is permanently backlogged, hence a
% work-conserving server never idles and the class takes exactly the capacity
% the other classes leave free: U_slc = 1 - sum_j U_j. This is exact whenever
% the other classes' utilizations are, and it reproduces the CTMC solution of
% the sanity SLC models to 5 significant digits. Note that solver_mna_closed
% and solver_mam_basic_mmap_closed instead pin U_slc = S*njobs*rate, which
% assumes each of the njobs customers owns a server and returns U>1 as soon as
% njobs exceeds nservers; that form is deliberately not reproduced here.
slcAll = find(sn.isslc);
if ~isempty(slcAll)
    % Utilization is the fraction of server time a class actually consumes, so
    % it follows from the declared service time. The inflation applied above
    % models the interference the self-looping customers cause in the queueing
    % delay, not extra work performed by the server, and must not be reported as
    % utilization (nor fed to the leftover-capacity identity below).
    for ist=1:M
        if slcjobs(ist) > 0
            for k=1:K
                if ~sn.isslc(k)
                    UN(ist,k) = Strue(ist,k)*TN(ist,k);
                end
            end
        end
    end
    QN(:,slcAll) = 0;
    UN(:,slcAll) = 0;
    RN(:,slcAll) = 0;
    TN(:,slcAll) = 0;
    for ist=1:M
        slck = slcAll(sn.refstat(slcAll) == ist);
        if isempty(slck)
            continue;
        end
        if isinf(sn.nservers(ist))
            % Delay station: no contention, every customer is always in
            % service, so the class completes at its full aggregate rate.
            for k=slck(:)'
                QN(ist,k) = sn.njobs(k);
                TN(ist,k) = sn.njobs(k)*sn.rates(ist,k);
                RN(ist,k) = Strue(ist,k);
                UN(ist,k) = Strue(ist,k)*TN(ist,k);
            end
        else
            nsrv = sn.nservers(ist);
            Uleft = max(0, 1 - sum(UN(ist,~sn.isslc)));
            % Several self-looping classes at one station share the free
            % capacity in proportion to the service rate they offer.
            w = sn.njobs(slck(:)').*sn.rates(ist,slck(:)');
            if sum(w) <= 0
                continue;
            end
            for i=1:numel(slck)
                k = slck(i);
                UN(ist,k) = min(Uleft*w(i)/sum(w), min(sn.njobs(k),nsrv)/nsrv);
                TN(ist,k) = sn.rates(ist,k)*UN(ist,k)*nsrv;
                QN(ist,k) = sn.njobs(k);
                if TN(ist,k) > 0
                    RN(ist,k) = QN(ist,k) ./ TN(ist,k);
                else
                    RN(ist,k) = 0;
                end
            end
        end
    end
    CN = sum(RN,1);
end
end