mvaDispatch.m

function out = mvaDispatch(self, sn, options)
% OUT = MVADISPATCH(SN, OPTIONS)
%
% One inner solve of the MVA analyzer dispatch: pick the analyzer that fits
% SN and return its metrics. This is the solve callback that
% @NetworkSolver/fjFixedPoint.m drives on each pass of the fork-join (MMT/HT)
% fixed point; on a model without forks it runs exactly once. The body is the
% dispatch that used to sit inline in runAnalyzer.
%
% OUT has fields QN, UN, RN, TN, CN, XN, lG, runtime, lastiter, method and
% actualmethod ('' when the analyzer did not report one).
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
            line_debug(options, 'Product-form check: hasProductForm=%d (exact method requested)', self.model.hasProductFormSolution);
            % OI/PAS are exactly solvable via solver_mva_oi_analyzer; exempt from the product-form guard (see _kb/06-solver-catalog.md, MVA OI dispatch)
            hasOIorPAS = any(sn.sched == SchedStrategy.OI | sn.sched == SchedStrategy.PAS);
            if strcmp(options.method,'exact')  && ~self.model.hasProductFormSolution && ~hasOIorPAS
                line_error(mfilename,'The exact method requires the model to have a product-form solution. This model does not have one.\nYou can use Network.hasProductFormSolution() to check before running the solver.\n Run the ''mva'' method to obtain an approximation based on the exact MVA algorithm.\n');
            end
            if strcmp(options.method,'mva') && ~self.model.hasProductFormSolution
                line_warning(mfilename,'The exact method requires the model to have a product-form solution. This model does not have one.\nYou can use Network.hasProductFormSolution() to check before running the solver.\nSolverMVA will return an approximation generated by an exact MVA algorithm.');
            end
 
            method = options.method;
 
            % Check for size-based policies (SRPT, PSJF, FB, LRPT, SETF) in multiclass open systems
            queueIdx = find(sn.nodetype == NodeType.Queue);
            isSizeBasedPolicy = false;
            if ~isempty(queueIdx)
                schedType = sn.sched(sn.nodeToStation(queueIdx(1)));
                isSizeBasedPolicy = ismember(schedType, [SchedStrategy.SRPT, SchedStrategy.PSJF, SchedStrategy.FB, SchedStrategy.LRPT, SchedStrategy.SETF]);
            end
 
            % Check for order-independent (OI) queues. An OI station is a PAS/OI
            % queue with an all-zero swap graph and a service-rate function;
            % detect it the same way the NC-oi path does (nc_is_oi_model),
            % not via a nodeparam flag (which is never populated).
            noi_idx = [];
            if ~any(isinf(sn.njobs))
                for ist = 1:sn.nstations
                    if sn.sched(ist) ~= SchedStrategy.PAS && sn.sched(ist) ~= SchedStrategy.OI
                        continue;
                    end
                    ind = sn.stationToNode(ist);
                    if ind < 1 || ind > numel(sn.nodeparam) || ~isstruct(sn.nodeparam{ind}) ...
                            || ~isfield(sn.nodeparam{ind}, 'swapGraph') ...
                            || ~isfield(sn.nodeparam{ind}, 'svcRateFun') || isempty(sn.nodeparam{ind}.svcRateFun)
                        continue;
                    end
                    sg = sn.nodeparam{ind}.swapGraph;
                    if isempty(sg) || any(sg(:) ~= 0)
                        continue;
                    end
                    noi_idx = ist;
                    break;
                end
            end
 
            ci_cache = find(sn.nodetype == NodeType.Cache, 1);
            hasRetrieval = ~isempty(ci_cache) && isfield(sn.nodeparam{ci_cache}, 'retrievalSystemCapacity') ...
                && sn.nodeparam{ci_cache}.retrievalSystemCapacity > 0;
            hasOIStation = any(sn.sched == SchedStrategy.OI | sn.sched == SchedStrategy.PAS);
            if ~isempty(noi_idx) && nc_is_oi_model(sn) && any(strcmpi(options.method,{'default','exact'}))
                % Order-independent queueing network
                line_debug(options, 'MVA: order-independent closed network, routing to solver_mva_oi_analyzer');
                [QN,UN,RN,TN,CN,XN,lG,runtime,lastiter,actualmethod] = solver_mva_oi_analyzer(sn, options);
            elseif hasOIStation
                % MVA supports OI/PAS stations only via the exact path above;
                % see _kb/06-solver-catalog.md (MVA section) for why AMVA cannot.
                line_error(mfilename, sprintf(['SolverMVA supports order-independent (OI) and pass-and-swap (PAS) stations only\n' ...
                    'through its exact order-independent analyzer, which requires method ''default'' or ''exact''\n' ...
                    '(got ''%s''), an empty/zero swap graph at every OI/PAS station, a closed model, and every\n' ...
                    'other station to be product-form (INF, PS, LCFS-PR, SIRO, or class-independent-rate FCFS).\n' ...
                    'Use SolverCTMC or SolverLDES for this model.'], options.method));
            elseif hasRetrieval % delayed-hit (retrieval-system) cache
                if any(sn.nodetype == NodeType.Source)
                    line_debug(options, 'Open delayed-hit retrieval cache, routing to mva_retrieval_analyzer');
                    [QN,UN,RN,TN,CN,XN,lG,hitprob,missprob,delayedprob,hitproblist,itemprob,latency,runtime,actualmethod] = solver_mva_retrieval_analyzer(sn, options);
                else
                    line_debug(options, 'Closed integrated delayed-hit retrieval cache, routing to mva_cacheqn_retrieval_analyzer');
                    [QN,UN,RN,TN,CN,XN,lG,hitprob,missprob,delayedprob,hitproblist,itemprob,latency,runtime,actualmethod] = solver_mva_cacheqn_retrieval_analyzer(sn, options);
                end
                lastiter = NaN;
                for ind = 1:sn.nnodes
                    if sn.nodetype(ind) == NodeType.Cache
                        self.model.nodes{ind}.setResultHitProb(hitprob);
                        self.model.nodes{ind}.setResultMissProb(missprob);
                        self.model.nodes{ind}.setResultDelayedHitProb(delayedprob);
                        self.model.nodes{ind}.setResultHitProbList(hitproblist);
                        self.model.nodes{ind}.setResultItemProb(itemprob);
                        self.model.nodes{ind}.setResultResidT(latency);
                    end
                end
                self.model.refreshStruct(true);
            elseif sn.nclosedjobs == 0 && length(sn.nodetype)==3 && all(sort(sn.nodetype)' == sort([NodeType.Source,NodeType.Queue,NodeType.Sink])) && isSizeBasedPolicy
                % Multiclass open system with size-based scheduling (SRPT, PSJF, FB, LRPT, SETF)
                line_debug(options, 'Size-based scheduling detected (%s), routing to qsys_sizebased_analyzer', char(schedType));
                [QN,UN,RN,TN,CN,XN,lG,runtime,lastiter,actualmethod] = solver_mva_qsys_sizebased_analyzer(sn, options, schedType);
            elseif sn.nclasses==1 && sn.nclosedjobs == 0 && length(sn.nodetype)==3 && all(sort(sn.nodetype)' == sort([NodeType.Source,NodeType.Queue,NodeType.Sink])) % is an open queueing system
                line_debug(options, 'Single-class open queueing system (Source-Queue-Sink), routing to qsys_analyzer');
                [QN,UN,RN,TN,CN,XN,lG,runtime,lastiter,actualmethod] = solver_mva_qsys_analyzer(sn, options);
            elseif sn.nclasses>1 && sn.nclosedjobs == 0 && length(sn.nodetype)==3 && all(sort(sn.nodetype)' == sort([NodeType.Source,NodeType.Queue,NodeType.Sink])) && sn.sched(find(sn.nodetype==NodeType.Queue)) == SchedStrategy.POLLING % is an open polling system
                line_debug(options, 'Multiclass open polling system, routing to polling_analyzer');
                [QN,UN,RN,TN,CN,XN,lG,runtime,lastiter,actualmethod] = solver_mva_polling_analyzer(sn, options);
            elseif sn.nclasses>1 && sn.nclosedjobs == 0 && length(sn.nodetype)==3 && all(sort(sn.nodetype)' == sort([NodeType.Source,NodeType.Queue,NodeType.Sink])) && sn.sched(find(sn.nodetype==NodeType.Queue)) == SchedStrategy.HOL && sn.nservers(sn.nodeToStation(sn.nodetype==NodeType.Queue)) == 1 && all(abs(sn.scv(sn.nodeToStation(sn.nodetype==NodeType.Source),:) - 1) < 1e-6 | ~isfinite(sn.scv(sn.nodeToStation(sn.nodetype==NodeType.Source),:)))
                % Multiclass open HOL (non-preemptive priority) M/G/1: exact Cobham
                % formula; the AMVA path applies the preemptive shadow-server formula.
                line_debug(options, 'Multiclass open HOL priority queue, routing to qsys_prio_analyzer');
                [QN,UN,RN,TN,CN,XN,lG,runtime,lastiter,actualmethod] = solver_mva_qsys_prio_analyzer(sn, options);
            elseif sn.nclasses>1 && sn.nclasses<=3 && sn.nclosedjobs == 0 && length(sn.nodetype)==3 && all(sort(sn.nodetype)' == sort([NodeType.Source,NodeType.Queue,NodeType.Sink])) && sn.sched(find(sn.nodetype==NodeType.Queue)) == SchedStrategy.DPS && sn.nservers(sn.nodeToStation(sn.nodetype==NodeType.Queue)) == 1 && all(abs(sn.scv(sn.nodeToStation(sn.nodetype==NodeType.Source),:) - 1) < 1e-6 | ~isfinite(sn.scv(sn.nodeToStation(sn.nodetype==NodeType.Source),:))) && all(abs(sn.scv(sn.nodeToStation(sn.nodetype==NodeType.Queue),:) - 1) < 1e-6 | ~isfinite(sn.scv(sn.nodeToStation(sn.nodetype==NodeType.Queue),:)))
                % Single open M/M/1-DPS queue: numerically exact DPS via the
                % truncated multiclass CTMC (qsys_mm1_dps). The AMVA-DPS cross-term
                % correction violates the equal-rate conservation law.
                line_debug(options, 'Multiclass open M/M/1-DPS queue, routing to exact qsys_mm1_dps');
                T0dps = tic;
                src_dps = sn.nodeToStation(sn.nodetype==NodeType.Source);
                q_dps = sn.nodeToStation(sn.nodetype==NodeType.Queue);
                lam_dps = sn.rates(src_dps,:);
                mu_dps = sn.rates(q_dps,:);
                w_dps = sn.schedparam(q_dps,:); w_dps(~(w_dps>0)) = 1;
                Tdps = qsys_mm1_dps(lam_dps, mu_dps, w_dps);
                QN = zeros(sn.nstations, sn.nclasses); UN = QN; RN = QN; TN = QN; CN = QN;
                XN = zeros(1, sn.nclasses);
                for rdps = 1:sn.nclasses
                    RN(q_dps,rdps) = Tdps(rdps); CN(q_dps,rdps) = Tdps(rdps);
                    XN(rdps) = lam_dps(rdps);
                    UN(q_dps,rdps) = lam_dps(rdps)/mu_dps(rdps);
                    TN(q_dps,rdps) = lam_dps(rdps); TN(src_dps,rdps) = lam_dps(rdps);
                    QN(q_dps,rdps) = lam_dps(rdps)*Tdps(rdps);
                end
                lG = 0; lastiter = 0; runtime = toc(T0dps);
                actualmethod = 'mm1.dps';
            elseif sn.nclosedjobs == 0 && length(sn.nodetype)==3 && all(sort(sn.nodetype)' == sort([NodeType.Source,NodeType.Cache,NodeType.Sink])) % is a non-rentrant cache
                line_debug(options, 'Non-reentrant cache (Source-Cache-Sink), routing to cache_analyzer');
                % random initialization
                for ind = 1:sn.nnodes
                    if sn.nodetype(ind) == NodeType.Cache
                        prob = self.model.nodes{ind}.server.hitClass;
                        prob(prob>0) = 0.5;
                        self.model.nodes{ind}.setResultHitProb(prob);
                        self.model.nodes{ind}.setResultMissProb(1-prob);
                    end
                end
                self.model.refreshChains();
                % start iteration
                [QN,UN,RN,TN,CN,XN,lG,runtime,lastiter,actualmethod,hitproblist,itemprob] = solver_mva_cache_analyzer(sn, options);
 
                for ind = 1:sn.nnodes
                    if sn.nodetype(ind) == NodeType.Cache
                        self.model.nodes{ind}.setResultHitProbList(hitproblist);
                        self.model.nodes{ind}.setResultItemProb(itemprob);
                        hitClass = self.model.nodes{ind}.getHitClass;
                        missClass = self.model.nodes{ind}.getMissClass;
                        hitprob = zeros(1,length(hitClass));
                        for k=1:length(self.model.nodes{ind}.getHitClass)
                            chain_k = sn.chains(:,k)>0;
                            inchain = sn.chains(chain_k,:)>0;
                            h = hitClass(k);
                            m = missClass(k);
                            if h>0 && m>0
                                hitprob(k) = XN(h) / sum(XN(inchain),"omitnan"); %#ok<NANSUM>
                            end
                        end
                        self.model.nodes{ind}.setResultHitProb(hitprob);
                        self.model.nodes{ind}.setResultMissProb(1-hitprob);
                    end
                end
                %self.model.refreshChains();
                self.model.refreshStruct(true);
            else % queueing network
                if any(sn.nodetype == NodeType.Cache) % if integrated caching-queueing
                    line_debug(options, 'Integrated caching-queueing network, routing to cacheqn_analyzer');
                    [QN,UN,RN,TN,CN,XN,lG,hitprob,missprob,runtime,lastiter] = solver_mva_cacheqn_analyzer(self, options);
                    for ind = 1:sn.nnodes
                        if sn.nodetype(ind) == NodeType.Cache
                            self.model.nodes{ind}.setResultHitProb(hitprob(ind,:));
                            self.model.nodes{ind}.setResultMissProb(missprob(ind,:));
                        end
                    end
                    %self.model.refreshChains();
                    self.model.refreshStruct(true);
                else % ordinary queueing network
                    switch method
                        case {'aba.upper', 'aba.lower', 'bjb.upper', 'bjb.lower', 'pb.upper', 'pb.lower', 'gb.upper', 'gb.lower', 'sb.upper', 'sb.lower', 'mwba.upper', 'mwba.lower'}
                            line_error(mfilename, ['Bound methods have moved to SolverBA. ' ...
                                'Use SolverBA(model,''method'',''%s'') (or LINE with that method) instead of SolverMVA.'], method);
                        case {'marie', 'amva.marie'}
                            line_debug(options, 'Using Marie aggregation-decomposition method');
                            [QN,UN,RN,TN,CN,XN,lG,runtime,lastiter,actualmethod] = solver_mva_marie_analyzer(sn, options);
                        otherwise
                            if ~isempty(sn.lldscaling) || ~isempty(sn.cdscaling) || ~isempty(sn.jdscaling)
                                line_debug(options, 'Load-dependent scaling detected (lldscaling=%d, cdscaling=%d), routing to mvald_analyzer', ~isempty(sn.lldscaling), ~isempty(sn.cdscaling));
                                [QN,UN,RN,TN,CN,XN,lG,runtime,lastiter,actualmethod] = solver_mvald_analyzer(sn, options);
                            else
                                subopts = options;
                                if self.model.hasFork && strcmp(subopts.method,'default')
                                    % force AMVA for a forked default model; see
                                    % _kb/06-solver-catalog.md (MVA section) for why
                                    subopts.method = 'amva';
                                end
                                line_debug(options, 'Standard queueing network, routing to mva_analyzer (method=%s)', subopts.method);
                                [QN,UN,RN,TN,CN,XN,lG,runtime,lastiter,actualmethod] = solver_mva_analyzer(sn, subopts);
                            end
                    end
                end
            end
 
out = struct('QN', QN, 'UN', UN, 'RN', RN, 'TN', TN, 'CN', CN, 'XN', XN, ...
    'lG', lG, 'runtime', runtime, 'lastiter', lastiter, 'method', method, ...
    'actualmethod', '');
if exist('actualmethod','var') && ~isempty(actualmethod)
    out.actualmethod = actualmethod;
end
end