runAnalyzer.m

function runtime = runAnalyzer(self, options)
% RUNTIME = RUNANALYZER(OPTIONS)
% Run the solver
 
if nargin<2
    options = self.getOptions;
end
 
% Wall-clock time-budget launch marker (see options.timeout / lineTimeoutExceeded)
options.timeout_tic = tic;
 
if strcmp(options.lang,'python')
    line_debug(options, 'MVA: using lang=python, delegating to native line_solver');
    [QN,UN,RN,TN,AN,WN,runtime] = PYLINE.getAvg(self.name, self.model, options);
    self.setAvgResults(QN,UN,RN,TN,AN,WN,[],[],runtime,options.method,NaN);
    return
end
 
self.runAnalyzerChecks(options);
 
sn = self.getStruct();
% Finite Capacity Region: MVA does not enforce the aggregate per-region job
% limit and would silently return the unconstrained product-form answer.
if isfield(sn,'nregions') && sn.nregions > 0
    line_error(mfilename,'This model uses a Finite Capacity Region (addRegion), which is not supported by SolverMVA. Use SolverJMT, or setCapacity for a single-station limit.');
end
if isfield(sn,'immfeed') && ~isempty(sn.immfeed) && any(sn.immfeed(:))
    line_warning(mfilename,'SolverMVA does not handle immediate feedback (immfeed); the solver will treat self-loops as class-switching with re-queueing.\n');
end
 
Solver.resetRandomGeneratorSeed(options.seed);
 
% Show library attribution if verbose and not yet shown
if options.verbose ~= VerboseLevel.SILENT && ~GlobalConstants.isLibraryAttributionShown()
    sn = self.getStruct();
    libs = SolverMVA.getLibrariesUsed(sn, options);
    if ~isempty(libs)
        line_printf('The solver will leverage %s.\n', strjoin(libs, ', '));
        GlobalConstants.setLibraryAttributionShown(true);
    end
end
 
iter = 0;
%options.lang='java';
 
switch options.lang
    case 'java'
        line_debug(options, 'MVA: using lang=java, delegating to JLINE');
        sn = getStruct(self); % doesn't need initial state
        jmodel = LINE2JLINE(self.model);
        %M = jmodel.getNumberOfStatefulNodes;
        M = jmodel.getNumberOfStations;
        R = jmodel.getNumberOfClasses;
        jsolver = JLINE.SolverMVA(jmodel, options);
        % The JLINE solver object is rebuilt on every call, so the JAR-side
        % fork-join (MMT) iterate cannot survive an outer iteration on its own.
        % Carry it across here, exactly as the lang='matlab' branch does with
        % self.fjForkLambda, so SolverLN layers warm-start identically in both
        % languages.
        if options.config.fj_warmstart && ~isempty(self.fjForkLambda)
            jsolver.setForkWarmStart(JLINE.from_line_matrix(self.fjForkLambda));
        end
        [QN,UN,RN,WN,AN,TN] = JLINE.arrayListToResults(jsolver.getAvgTable);
        if self.model.hasFork
            self.fjForkLambda = JLINE.from_jline_matrix(jsolver.getForkWarmStart());
        end
        runtime = jsolver.result.runtime;
        CN = [];
        XN = [];
        QN = reshape(QN',R,M)';
        UN = reshape(UN',R,M)';
        RN = reshape(RN',R,M)';
        TN = reshape(TN',R,M)';
        WN = reshape(WN',R,M)';
        AN = reshape(AN',R,M)';
        lG = NaN;
        lastiter = NaN;
        for ind = 1:sn.nnodes
            if sn.nodetype(ind) == NodeType.Cache
                jnode = jmodel.getNodeByIndex(ind-1);
                self.model.nodes{ind}.setResultHitProb(JLINE.from_jline_matrix(jnode.getHitRatio()));
                self.model.nodes{ind}.setResultMissProb(JLINE.from_jline_matrix(jnode.getMissRatio()));
                % Retrieval-cache extras (delayed-hit ratio, per-list hit ratio
                % and expected latency) so getAvgCacheTable matches the native path.
                self.model.nodes{ind}.setResultDelayedHitProb(JLINE.from_jline_matrix(jnode.getDelayedHitRatio()));
                self.model.nodes{ind}.setResultHitProbList(JLINE.from_jline_matrix(jnode.getHitRatioByList()));
                self.model.nodes{ind}.setResultItemProb(JLINE.from_jline_matrix(jnode.getItemProb()));
                self.model.nodes{ind}.setResultResidT(JLINE.from_jline_matrix(jnode.getResidT()));
            end
        end
        %self.model.refreshChains();
        self.model.refreshStruct(true);
        self.setAvgResults(QN,UN,RN,TN,AN,WN,CN,XN,runtime,options.method,lastiter);
        self.result.Prob.logNormConstAggr = lG;
        return
    case 'matlab'
        line_debug(options, 'MVA: using lang=matlab');
        sn = getStruct(self); % doesn't need initial state
        forkLoop = true;
        forkIter = 0;
        % create artificial classes arrival rates. Under an outer iteration
        % (e.g. SolverLN) the MMT fixed point is re-solved once per outer
        % iteration on a model whose parameters move only slightly; restarting
        % from FineTol each time discards the previous converged point and
        % forces the inner loop to re-converge from scratch. Warm-start from
        % the retained iterate when it is conformant.
        forkLambda = GlobalConstants.FineTol * ones(1, 2*sn.nclasses*sum(sn.nodetype==NodeType.Fork));
        if options.config.fj_warmstart && ~isempty(self.fjForkLambda) ...
                && isequal(size(self.fjForkLambda), size(forkLambda))
            forkLambda = self.fjForkLambda;
        end
        QN = GlobalConstants.Immediate * ones(1, sn.nclasses);
        QN_1 = 0*QN;
        UN = 0*QN;
        while (forkLoop && forkIter < options.iter_max)
            if self.model.hasFork
                forkIter = forkIter + 1;
                line_debug(options, 'Fork-join iteration %d', forkIter);
                if forkIter == 1
                    switch options.config.fork_join
                        case {'heidelberger-trivedi', 'ht'}
                            [nonfjmodel, fjclassmap, fjforkmap, fj_auxiliary_delays] = ModelAdapter.ht(self.model);
                            line_debug(options, 'Fork-join method: heidelberger-trivedi');
                        case {'mmt', 'default', 'fjt'}
                            % Reuse the MMT transformation across outer
                            % iterations when nothing structural has changed
                            % since it was built (SolverLN re-solves each layer
                            % once per fixed-point iteration, changing only the
                            % rates). The transformation is a function of the
                            % fork topology, so structVersion is the precise
                            % condition; rebuilding costs a full model.copy()
                            % every time. The cached model still holds the
                            % previous iteration's rates and sync delays, so it
                            % must be re-fed from the base model first.
                            cacheUsable = ~isempty(self.mmtCache) && ...
                                self.mmtCache.structVersion == self.model.structVersion && ...
                                ModelAdapter.refreshServicesFromBase(self.mmtCache.nonfjmodel, self.mmtCache.prov);
                            if cacheUsable
                                nonfjmodel = self.mmtCache.nonfjmodel;
                                fjclassmap = self.mmtCache.fjclassmap;
                                fjforkmap = self.mmtCache.fjforkmap;
                                fanout = self.mmtCache.fanout;
                                outer_forks = self.mmtCache.outer_forks;
                                parent_forks = self.mmtCache.parent_forks;
                                nonfjmodel.refreshRates();
                                line_debug(options, 'Fork-join method: mmt (cached transformation)');
                            else
                                [nonfjmodel, fjclassmap, fjforkmap, fanout, prov] = ModelAdapter.mmt(self.model, forkLambda);
                                line_debug(options, 'Fork-join method: mmt');
                                [outer_forks, parent_forks] = ModelAdapter.sortForks(sn, fjforkmap, fjclassmap, nonfjmodel);
                                self.mmtCache = struct('nonfjmodel', nonfjmodel, 'prov', prov, ...
                                    'fjclassmap', fjclassmap, 'fjforkmap', fjforkmap, 'fanout', fanout, ...
                                    'outer_forks', outer_forks, 'parent_forks', parent_forks, ...
                                    'structVersion', self.model.structVersion);
                            end
                    end
                elseif ~strcmp(options.config.fork_join, 'heidelberger-trivedi') & ~strcmp(options.config.fork_join, 'ht')
                    %line_printf('Fork-join iteration %d\n',forkIter);
                    nonfjSource = nonfjmodel.getSource;
                    for r=1:length(fjclassmap) % r is the auxiliary class
                        s = fjclassmap(r);
                        if s>0
                            if fanout(r)>0
                                if ~nonfjSource.arrivalProcess{r}.isDisabled
                                    nonfjSource.arrivalProcess{r}.setRate((fanout(r)-1)*forkLambda(r));
                                end
                            end
                        end
                    end
                    nonfjmodel.refreshRates();
                end
                sn = nonfjmodel.getStruct(false); % this ensures that we solve nonfjmodel instead of the original model
                line_debug(options, 'Fork-join iter %d: rebuilt nonfjmodel struct (nstations=%d, nclasses=%d)', forkIter, sn.nstations, sn.nclasses);
                % Mixed absolute/relative convergence test on the MMT iterate.
                % A purely relative test, max|1 - QN_1./QN| < CoarseTol, cannot
                % certify convergence to zero. When the auxiliary class
                % throughput at the join is ~0, forkLambda decays geometrically
                % under the mean() damping above, so QN halves on every
                % iteration and the relative change stays pinned at 100% down to
                % denormals. The loop then ran to options.iter_max and returned a
                % non-converged fixed point. An entry below the absolute floor is
                % numerically zero and has converged.
                qn_converged = abs(QN_1 - QN) <= GlobalConstants.Zero + GlobalConstants.CoarseTol*abs(QN);
                if isequal(size(QN_1), size(QN)) && all(qn_converged(:)) && (forkIter > 2)
                    line_debug(options, 'Fork-join iter %d: converged (mixed abs/rel test)', forkIter);
                    forkLoop = false;
                else
                    if self.model.hasOpenClasses
                        sourceIndex = self.model.getSource.index;
                        UNnosource = UN; UNnosource(sourceIndex,:) = 0;
                        if any(find(sum(UNnosource(:,isinf(sn.njobs(1:size(QN,2)))),2)>0.99 * sn.nservers))
                            line_warning(mfilename,'The model may be unstable: the utilization of station %i for open classes exceeds 99 percent.\n',maxpos(sum(UNnosource,2)));
                        end
                    end
                    QN_1 = QN;
                end
            else
                forkLoop = false;
            end
            line_debug(options, 'Product-form check: hasProductForm=%d (exact method requested)', self.model.hasProductFormSolution);
            if strcmp(options.method,'exact')  && ~self.model.hasProductFormSolution
                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
                % An OI/PAS station carries a rank-rate function mu(n) of the whole
                % per-class occupancy. The AMVA iteration only ever sees sn.rates
                % (the single-job rates), so it cannot represent such a station and
                % would silently return a zero queue-length there. MVA therefore
                % supports OI/PAS stations only via the exact path above.
                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)
                                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')
                                    % The fork-join transform (ht/mmt) yields a mixed
                                    % model with auxiliary near-zero-rate open classes.
                                    % The default dispatch would route it to exact mixed
                                    % MVA (pfqn_mvamx), which degenerates to zero here;
                                    % the fork-join approximation requires the AMVA
                                    % linearizer, matching the pre-mixed-exact behaviour.
                                    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
            if self.model.hasFork
                line_debug(options, 'Fork-join post-processing: method=%s, computing sync delays', options.config.fork_join);
                nonfjstruct = sn;
                sn = self.getStruct;
                % Pre-compute linked routing matrix once (loop-invariant)
                switch options.config.fork_join
                    case {'mmt', 'default', 'fjt'}
                        Pcs = cell2mat(nonfjmodel.getLinkedRoutingMatrix);
                end
                for f=find(sn.nodetype == NodeType.Fork)'
                    switch options.config.fork_join
                        case {'mmt', 'default', 'fjt'}
                            TNfork = zeros(1,sn.nclasses);
                            for c=1:sn.nchains
                                inchain = find(sn.chains(c,:));
                                for r=inchain(:)'
                                    TNfork(r) =  (sn.nodevisits{c}(parent_forks(f),r) / sum(sn.visits{c}(sn.stationToStateful(sn.refstat(r)),inchain))) * sum(TN(sn.refstat(r),inchain));
                                end
                            end
                            % find join associated to fork f
                            joinIdx = find(sn.fj(f,:));
                            forkauxclasses = find(fjforkmap==f);
                            for s=forkauxclasses(:)'
                                r = fjclassmap(s); % original class associated to auxiliary class s
                                if isempty(joinIdx)
                                    forkLambda(s) = mean([forkLambda(s); TNfork(r)],1);
                                else
                                    joinStat = sn.nodeToStation(joinIdx);
                                    TN(sn.nodeToStation(joinIdx),r) = TN(sn.nodeToStation(joinIdx),r) + sum(TN(sn.nodeToStation(joinIdx), find(fjclassmap == r))) - TN(sn.nodeToStation(joinIdx), s);
                                    forkLambda(s) = mean([forkLambda(s); TN(sn.nodeToStation(joinIdx),r)],1);
                                end
                                if isempty(joinIdx) || ~outer_forks(f, r)
                                    % No join nodes for this fork, no synchronisation delay
                                    continue;
                                end
                                % Find the parallel paths coming out of the fork
                                ri = ModelAdapter.findPathsCS(sn, Pcs, f, joinIdx, r, [r,s], QN, TN, 0, fjclassmap, fjforkmap, nonfjmodel);
                                if isempty(ri)
                                    % No routing from fork for this class - set sync delay to 0 (Immediate)
                                    % Matches JAR behavior where empty Matrix gives syncDelay = 0
                                    syncDelay = 0;
                                else
                                    lambdai = 1./ri;
                                    d0 = 0;
                                    parallel_branches = length(ri);
                                    for pow=0:(parallel_branches - 1)
                                        current_sum = sum(1./sum(nchoosek(lambdai, pow + 1),2));
                                        d0 = d0 + (-1)^pow * current_sum;
                                    end
                                    syncDelay = d0*sn.nodeparam{f}.fanOut - mean(ri);
                                end
                                % Set the synchronisation delays
                                nonfjmodel.nodes{joinIdx}.setService(nonfjmodel.classes{s}, Exp.fitMean(syncDelay));
                                if outer_forks(f, r)
                                    nonfjmodel.nodes{joinIdx}.setService(nonfjmodel.classes{r}, Exp.fitMean(syncDelay));
                                end
                            end
                        case {'heidelberger-trivedi', 'ht'}
                            joinIdx = find(sn.fj(f,:));
                            for c=1:sn.nchains
                                inchain = find(sn.chains(c,:));
                                for r=inchain(:)'
                                    if sn.nodevisits{c}(f,r) == 0
                                        continue;
                                    end
                                    % Obtain the response times on the parallel branches
                                    ri = RN(:, find(fjclassmap == r));
                                    ri(isnan(ri) | isinf(ri)) = 0;
                                    ri = sum(ri, 1, "omitnan") - RN(nonfjstruct.nodeToStation(fj_auxiliary_delays{joinIdx}), find(fjclassmap == r)) - RN(nonfjstruct.nodeToStation(joinIdx), find(fjclassmap == r));
                                    lambdai = 1./ri;
                                    d0 = 0;
                                    parallel_branches = length(self.model.nodes{f}.output.outputStrategy{r}{3});
                                    for pow=0:(parallel_branches - 1)
                                        current_sum = sum(1./sum(nchoosek(lambdai, pow + 1),2));
                                        d0 = d0 + (-1)^pow * current_sum;
                                    end
                                    di = d0*sn.nodeparam{f}.fanOut - ri;
                                    r0 = sum(RN(:, inchain), 2);
                                    r0(isnan(r0) | isinf(r0)) = 0;
                                    r0 = sum(r0, 1, "omitnan") - RN(nonfjstruct.nodeToStation(joinIdx), r);
                                    % Update the delays at the join node and at the auxiliary delay
                                    nonfjmodel.nodes{joinIdx}.setService(nonfjmodel.classes{r}, Exp.fitMean(d0*sn.nodeparam{f}.fanOut));
                                    idx = 1;
                                    for s=find(fjclassmap == r)
                                        nonfjmodel.nodes{joinIdx}.setService(nonfjmodel.classes{s}, Exp.fitMean(di(idx)));
                                        idx = idx + 1;
                                        nonfjmodel.nodes{fj_auxiliary_delays{joinIdx}}.setService(nonfjmodel.classes{s}, Exp.fitMean(r0));
                                    end
 
                                end
                            end
                    end
                end
                % Batch refreshRates after all sync delay updates (moved out of inner loops for performance)
                switch options.config.fork_join
                    case {'mmt', 'default', 'fjt'}
                        nonfjmodel.refreshRates();
                end
                switch options.config.fork_join
                    case {'heidelberger-trivedi', 'ht'}
                        nonfjmodel.refreshStruct();
                        % Delete the queue lengths, response times, throughputs and utilizations of the original classes at the join nodes
                        QN(nonfjstruct.nodeToStation(find(sn.nodetype == NodeType.Join)), nonzeros(fjclassmap)) = 0;
                        RN(nonfjstruct.nodeToStation(find(sn.nodetype == NodeType.Join)), nonzeros(fjclassmap)) = 0;
                        % Save the throughputs of the original classes at the join node
                        TN_orig = TN(nonfjstruct.nodeToStation(find(sn.nodetype == NodeType.Join)), nonzeros(fjclassmap));
                        TN(nonfjstruct.nodeToStation(find(sn.nodetype == NodeType.Join)), nonzeros(fjclassmap)) = 0;
                        UN(nonfjstruct.nodeToStation(find(sn.nodetype == NodeType.Join)), nonzeros(fjclassmap)) = 0;
 
                        % Remove the times at the auxiliary delay
                        QN(nonfjstruct.nodeToStation(cell2mat(fj_auxiliary_delays)),:) = [];
                        UN(nonfjstruct.nodeToStation(cell2mat(fj_auxiliary_delays)),:) = [];
                        RN(nonfjstruct.nodeToStation(cell2mat(fj_auxiliary_delays)),:) = [];
                        TN(nonfjstruct.nodeToStation(cell2mat(fj_auxiliary_delays)),:) = [];
                        % merge back artificial classes into their original classes
                        for r=1:length(fjclassmap)
                            s = fjclassmap(r);
                            if s>0
                                QN(:,s) = QN(:,s) + QN(:,r);
                                UN(:,s) = UN(:,s) + UN(:,r);
                                % Add all throughputs of the auxiliary classes to facilitate the computation of the response times
                                TN(:,s) = TN(:,s) + TN(:,r);
                                RN(:,s) = QN(:,s) ./ TN(:,s);
                            end
                        end
                        % Re-set the throughputs for the original classes
                        TN(nonfjstruct.nodeToStation(find(sn.nodetype == NodeType.Join)), nonzeros(fjclassmap)) = TN_orig;
                    case {'mmt', 'default', 'fjt'}
                        TN_orig = TN([nonfjstruct.nodeToStation(find(sn.nodetype == NodeType.Join)), nonfjstruct.nodeToStation(find(sn.nodetype == NodeType.Source))], nonzeros(fjclassmap));
                        % merge back artificial classes into their original classes
                        for r=1:length(fjclassmap)
                            s = fjclassmap(r);
                            if s>0
                                QN(:,s) = QN(:,s) + QN(:,r);
                                UN(:,s) = UN(:,s) + UN(:,r);
                                TN(:,s) = TN(:,s) + TN(:,r);
                                %RN(:,s) = RN(:,s) + RN(:,r);
                                % for i=find(snorig.nodetype == NodeType.Delay | snorig.nodetype == NodeType.Queue)'
                                %     TN(snorig.nodeToStation(i),s) = TN(snorig.nodeToStation(i),s) + TN(snorig.nodeToStation(i),r);
                                % end
                                RN(:,s) = QN(:,s) ./ TN(:,s);
                                %CN(:,s) = CN(:,s) + CN(:,r);
                                %XN(:,s) = XN(:,s) + XN(:,r);
                            end
                        end
                        TN([nonfjstruct.nodeToStation(find(sn.nodetype == NodeType.Join)), nonfjstruct.nodeToStation(find(sn.nodetype == NodeType.Source))], nonzeros(fjclassmap)) = TN_orig;
                end
                QN(:,fjclassmap>0) = [];
                UN(:,fjclassmap>0) = [];
                RN(:,fjclassmap>0) = [];
                TN(:,fjclassmap>0) = [];
                CN(:,fjclassmap>0) = [];
                XN(:,fjclassmap>0) = [];
            end
            iter = iter + lastiter;
            % Cap accumulated iterations for stability (Python parity)
            if iter > 10000
                iter = 10000;
                break; % Exit forkLoop early
            end
        end
 
        % The fork-join loop previously exhausted options.iter_max silently:
        % convergence was only ever reported through line_debug, so a
        % non-converged MMT fixed point was returned as a normal result.
        if self.model.hasFork
            if forkLoop && forkIter >= options.iter_max
                line_warning(mfilename,'The fork-join (%s) fixed point did not converge in options.iter_max=%d iterations; returning the interim solution.\n', options.config.fork_join, options.iter_max);
            end
            % Retain the MMT iterate so that a subsequent runAnalyzer call on
            % this solver (an outer LN iteration) resumes the fixed point here.
            self.fjForkLambda = forkLambda;
        end
 
        sn = self.model.getStruct();
 
        % Compute average residence time at steady-state
        AN = sn_get_arvr_from_tput(sn, TN, self.getAvgTputHandles());
        WN = sn_get_residt_from_respt(sn, RN, self.getAvgResidTHandles());
        if strcmp(method,'default') && exist('actualmethod','var')
            self.setAvgResults(QN,UN,RN,TN,AN,WN,CN,XN,runtime,['default/' actualmethod],iter);
        else
            self.setAvgResults(QN,UN,RN,TN,AN,WN,CN,XN,runtime,method,iter);
        end
        self.result.Prob.logNormConstAggr = lG;
        if lineTimeoutExceeded(options)
            self.result.Avg.timedOut = true;
            line_warning(mfilename,'Solver stopped after the wall-clock time budget (options.timeout=%gs) was exceeded; returning the interim solution.\n', options.timeout);
        end
end
end