runAnalyzer.m

function runtime = runAnalyzer(self, options)
% RUNTIME = RUN()
% 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, 'NC: 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
 
iter = NaN;
self.runAnalyzerChecks(options);
 
sn = self.getStruct();
if isfield(sn,'immfeed') && ~isempty(sn.immfeed) && any(sn.immfeed(:))
    line_warning(mfilename,'SolverNC 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 = SolverNC.getLibrariesUsed(sn, options);
    if ~isempty(libs)
        line_printf('The solver will leverage %s.\n', strjoin(libs, ', '));
        GlobalConstants.setLibraryAttributionShown(true);
    end
end
 
origmethod = options.method;
 
% Warn (debug level) when the requested method is unrecognized. With solver
% checks disabled (e.g. when SolverNC is used as a SolverLN layer backend,
% where enableChecks=false bypasses runAnalyzerChecks) an unknown method name
% silently falls through to the default normalizing-constant analyzer instead
% of raising an error, which masks typos such as 'adaptive'.
validMethods = [self.listValidMethods(), {'comomld'}]; % comomld is auto-selected internally from 'default'
if ~any(strcmpi(origmethod, validMethods))
    line_debug(options, 'NC: unrecognized method ''%s'', falling back to the default normalizing-constant analyzer (nc_analyzer/comom).', origmethod);
end
 
%options.lang = 'java';
 
switch options.lang
    case 'java'
        line_debug(options, 'NC: using lang=java, delegating to JLINE');
        sn = self.getStruct;
        jmodel = LINE2JLINE(self.model);
        %M = jmodel.getNumberOfStatefulNodes;
        M = jmodel.getNumberOfStations;
        R = jmodel.getNumberOfClasses;
        jsolver = JLINE.SolverNC(jmodel, options);
        [QN,UN,RN,WN,AN,TN] = JLINE.arrayListToResults(jsolver.getAvgTable);
        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, 'NC: using lang=matlab');
        sn = getStruct(self); % doesn't need initial state
 
        % Order-independent (OI) closed network: solved exactly by the
        % balanced-fairness normalizing constant (pfqn_oi_nc) with the OI
        % functional-server (pfqn_oi_fnc) identity for the mean queue lengths.
        % Intercept before the multiserver->lldscaling conversion below, since
        % an OI station is multiserver but is not a plain min(n,c) load-
        % dependent station.
        if nc_is_oi_model(sn) && any(strcmpi(options.method,{'default','exact'}))
            line_debug(options, 'NC: order-independent closed network, routing to solver_nc_oi_analyzer');
            [QN,UN,RN,TN,CN,XN,lG,runtime,iter,actualmethod] = solver_nc_oi_analyzer(sn, options);
            AN = sn_get_arvr_from_tput(sn, TN, self.getAvgTputHandles());
            if strcmp(origmethod,'default') && ~strcmp(actualmethod,'default')
                self.setAvgResults(QN,UN,RN,TN,AN,[],CN,XN,runtime,['default/' actualmethod],iter);
            else
                self.setAvgResults(QN,UN,RN,TN,AN,[],CN,XN,runtime,actualmethod,iter);
            end
            self.result.Prob.logNormConstAggr = real(lG);
            return
        end
 
        switch options.method
            case 'default'
                if sn.nstations == 2 && ~any(sn.nodetype == NodeType.Cache) && any(sn.nodetype == NodeType.Delay) && any(sn.nservers(isfinite(sn.nservers))>1)
                    % 2-station Delay+multiserver (the topology of every SolverLN
                    % layer submodel). When the model is product-form the
                    % multiserver can be solved exactly via load-dependent CoMoM
                    % (convert to mu(n)=min(n,c) and route to the ncld analyzer),
                    % which matches CTMC/MVA exactly. Non-product-form cases (e.g.
                    % LN layers with heterogeneous per-class FCFS rates) have no
                    % exact CoMoM-LD solution and fall back to Seidmann's
                    % approximation (comom).
                    if self.model.hasProductFormSolution && isempty(sn.lldscaling)
                        Nt = sum(sn.njobs);
                        if isfinite(Nt)
                            sn.lldscaling = ones(sn.nstations,Nt);
                            for i=1:sn.nstations
                                if sn.nservers(i) > 1 && isfinite(sn.nservers(i))
                                    % The queueing solve reads mu(n) from lldscaling, so the
                                    % server count must be kept: utilization is the fraction
                                    % of the c servers busy, and c is not recoverable from
                                    % lldscaling once Nt<c (min(1:Nt,c) is then just 1:Nt).
                                    % Zeroing it to 1 made solver_ncld normalize U by
                                    % max(lldscaling)=min(Nt,c) instead of c, overstating U
                                    % by c/Nt whenever the population is below the servers.
                                    sn.lldscaling(i,:) = min(1:Nt,sn.nservers(i));
                                end
                            end
                        end
                        line_debug(options, 'Default method: 2-station Delay+multiserver product-form network, using exact load-dependent comomld');
                    else
                        options.method = 'comom'; % non-product-form (e.g. LN layer submodels): Seidmann approximation
                        line_debug(options, 'Default method: 2-station Delay+multiserver non-product-form network, using comom');
                    end
                end
            case {'exact','is'}
                % 'is' (importance sampling) on an ORDINARY product-form model is
                % a normalizing-constant estimator and needs the same model as
                % 'exact': convert multiserver stations to load-dependent here
                % too, so it routes to solver_ncld_analyzer -> pfqn_ncld ->
                % pfqn_ld_is rather than falling through to Seidmann's
                % approximation, which would faithfully estimate a DIFFERENT
                % (approximated) model.
                %
                % OI / pass-and-swap models are the exception: an OI station is
                % multiserver but is not a plain min(n,c) load-dependent station,
                % and a P&S tandem has only per-communicating-class product form
                % (so hasProductFormSolution is false). Leave them untouched for
                % solver_nc_analyzer to route to pfqn_pas_is / pfqn_oi_is.
                if strcmpi(options.method,'is') && nc_is_pas_model(sn)
                    % no-op: handled by solver_nc_analyzer (pfqn_pas_is)
                elseif ~self.model.hasProductFormSolution
                    line_error(mfilename,'The %s method requires the model to have a product-form solution. This model does not have one. You can use Network.hasProductFormSolution() to check before running the solver.', options.method);
                elseif isempty(sn.lldscaling) && any(sn.nservers(isfinite(sn.nservers)) > 1)
                    % if exact/is is requested and does not override a lldscaling assigment
                    % ONLY when a genuine multiserver station is present. Setting
                    % lldscaling=ones on an all-single-server model is a semantic
                    % no-op but forces the load-dependent path (solver_ncld_analyzer),
                    % whose recovery mishandles a single-station-confined
                    % (self-looping) closed chain -- it returned QLen=[1,0,0] for a
                    % Delay+2Queue model with self-looping classes where the exact
                    % answer (CTMC / solver_nc_analyzer) is [0.71,0.14,0.14]. With
                    % no multiserver there is nothing to convert, so leave
                    % lldscaling empty and route to the standard normalizing-constant
                    % path.
                    Nt = sum(sn.njobs);
                    if isfinite(Nt)
                        % trasform multi-server nodes into lld nodes; the server count is
                        % kept so that utilization stays normalized by c (see the
                        % 'default' branch above)
                        sn.lldscaling = ones(sn.nstations,Nt);
                        for i=1:sn.nstations
                            if sn.nservers(i) > 1 && isfinite(sn.nservers(i))
                                sn.lldscaling(i,:) = min(1:Nt,sn.nservers(i));
                            end
                        end
                    end
                    line_debug(options, 'Exact method: converted multiserver stations to load-dependent (Nt=%d)', Nt);
                end
        end
 
        % The feature gate ran in runAnalyzerChecks above, through
        % supportsModelMethod, and reports which features are unsupported. A
        % coarse supports(self.model) repeat here would only re-reject the same
        % models with a message that names nothing.
 
        Solver.resetRandomGeneratorSeed(options.seed);
 
        ci_cache = find(sn.nodetype == NodeType.Cache, 1);
        hasRetrieval = ~isempty(ci_cache) && isfield(sn.nodeparam{ci_cache}, 'retrievalSystemCapacity') ...
            && sn.nodeparam{ci_cache}.retrievalSystemCapacity > 0;
        if hasRetrieval % delayed-hit (retrieval-system) cache
            if any(sn.nodetype == NodeType.Source)
                line_debug(options, 'Open delayed-hit retrieval cache, routing to nc_retrieval_analyzer');
                [QN,UN,RN,TN,CN,XN,lG,hitprob,missprob,delayedprob,hitproblist,itemprob,latency,runtime,actualmethod] = solver_nc_retrieval_analyzer(sn, options);
            else
                line_debug(options, 'Closed integrated delayed-hit retrieval cache, routing to nc_cacheqn_retrieval_analyzer');
                [QN,UN,RN,TN,CN,XN,lG,hitprob,missprob,delayedprob,hitproblist,itemprob,latency,runtime,actualmethod] = solver_nc_cacheqn_retrieval_analyzer(sn, options);
            end
            iter = 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.Cache,NodeType.Sink])) % is a non-rentrant cache
            line_debug(options, 'Non-reentrant cache (Source-Cache-Sink), routing to nc_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,pij,runtime,actualmethod,hitproblist,itemprob] = solver_nc_cache_analyzer(sn, options);
            self.result.Prob.itemProb = pij;
            for ind = 1:sn.nnodes
                if sn.nodetype(ind) == NodeType.Cache
                    self.model.nodes{ind}.setResultHitProbList(hitproblist);
                    self.model.nodes{ind}.setResultItemProb(itemprob);
                    %prob = self.model.nodes{ind}.server.hitClass;
                    %prob(prob>0) = 0.5;
                    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)
                        %                for k=1:length(self.model.nodes{ind}.server.hitClass)
                        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");
                        end
                    end
                    self.model.nodes{ind}.setResultHitProb(hitprob);
                    self.model.nodes{ind}.setResultMissProb(1-hitprob);
                end
            end
            self.model.refreshChains;
        else % queueing network
            if any(sn.nodetype == NodeType.Cache) % if integrated caching-queueing
                line_debug(options, 'Integrated caching-queueing network, routing to nc_cacheqn_analyzer');
                [QN,UN,RN,TN,CN,XN,lG,hitprob,missprob,runtime,iter,actualmethod] = solver_nc_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.refreshStruct(true);  % Force refresh to get updated actualhitprob/actualmissprob
                sn = self.model.sn;
            else % ordinary queueing network
                % Check for open model with single FCR containing single Delay (loss network)
                if ~sn_has_closed_classes(sn) && sn.nregions == 1
                    regionMatrix = sn.region{1};
                    % Find stations in FCR (those with non-negative constraints)
                    stationsInFCR = find(any(regionMatrix(:,1:end-1) >= 0, 2) | regionMatrix(:,end) >= 0);
                    if length(stationsInFCR) == 1 && isinf(sn.nservers(stationsInFCR(1)))
                        % Single delay node in FCR - check drop rule
                        if all(sn.regionrule(1,:) == DropStrategy.DROP)
                            % Use loss network solver
                            line_debug(options, 'Open model with single FCR + Delay node (DROP), routing to nc_lossn_analyzer');
                            [QN,UN,RN,TN,CN,XN,lG,runtime,iter,actualmethod] = solver_nc_lossn_analyzer(sn, options);
                            AN = sn_get_arvr_from_tput(sn, TN, self.getAvgTputHandles());
                            if strcmp(origmethod,'default') && exist('actualmethod','var') && ~strcmp(actualmethod,'default')
                                self.setAvgResults(QN,UN,RN,TN,AN,[],CN,XN,runtime,['default/' actualmethod],iter);
                            else
                                self.setAvgResults(QN,UN,RN,TN,AN,[],CN,XN,runtime,actualmethod,iter);
                            end
                            self.result.Prob.logNormConstAggr = real(lG);
                            return
                        else
                            % WAITQ (blocking) not supported - error and stop
                            line_error(mfilename, 'SolverNC does not support finite capacity regions with WAITQ (blocking) policy. Use DROP policy instead.');
                        end
                    end
                end
                % Residual FCR (not the single-Delay loss-network case dispatched
                % above): NC does not enforce the aggregate region limit; reject
                % instead of silently returning the unconstrained answer.
                if isfield(sn,'nregions') && sn.nregions > 0
                    line_error(mfilename,'This model uses a Finite Capacity Region (addRegion) on queueing stations, which is not supported by SolverNC. Use SolverJMT, or setCapacity for a single-station limit.');
                end
                if ~isempty(sn.lldscaling) || ~isempty(sn.cdscaling)
                    line_debug(options, 'Load-dependent scaling detected, routing to ncld_analyzer');
                    [QN,UN,RN,TN,CN,XN,lG,runtime,iter,actualmethod] = solver_ncld_analyzer(sn, options);
                else
                    switch options.method
                        case 'exact'
                            if ~sn_has_open_classes(sn)
                                % This branch is reached only with EMPTY lldscaling
                                % (the outer test routes any load-dependent model to
                                % ncld). Genuine multiserver stations were converted
                                % to lldscaling above, so an empty lldscaling here
                                % means an all-single-server closed model, which the
                                % standard normalizing-constant path solves exactly.
                                % Routing it to solver_ncld_analyzer instead
                                % mishandles a single-station-confined (self-looping)
                                % closed chain -- returning QLen=[1,0,0] for a
                                % Delay+2Queue model with self-looping classes where
                                % the exact answer is [0.71,0.14,0.14].
                                line_debug(options, 'NC method=exact, single-server closed model, routing to nc_analyzer');
                                [QN,UN,RN,TN,CN,XN,lG,runtime,iter,actualmethod] = solver_nc_analyzer(sn, options);
                            else%if ~snHasClosedClasses(sn)
                                line_debug(options, 'NC method=exact, open model, routing to nc_analyzer');
                                [QN,UN,RN,TN,CN,XN,lG,runtime,iter,actualmethod] = solver_nc_analyzer(sn, options);
                            end
                        case {'rd','nrp','nrl','comomld'}
                            line_debug(options, 'NC method=%s, routing to ncld_analyzer', options.method);
                            [QN,UN,RN,TN,CN,XN,lG,runtime,iter,actualmethod] = solver_ncld_analyzer(sn, options);
                        otherwise
                            line_debug(options, 'NC method=%s, routing to nc_analyzer', options.method);
                            [QN,UN,RN,TN,CN,XN,lG,runtime,iter,actualmethod] = solver_nc_analyzer(sn, options);
                    end
                end
            end
        end
        % Compute average arrival rate at steady-state
        AN = sn_get_arvr_from_tput(sn, TN, self.getAvgTputHandles());
end
if strcmp(origmethod,'default') && exist('actualmethod','var') && ~strcmp(actualmethod,'default')
    self.setAvgResults(QN,UN,RN,TN,AN,[],CN,XN,runtime,['default/' actualmethod],iter);
else
    self.setAvgResults(QN,UN,RN,TN,AN,[],CN,XN,runtime,actualmethod,iter);
end
 
self.result.Prob.logNormConstAggr = real(lG);
if lineTimeoutExceeded(options)
    self.result.Avg.timedOut = true;
    line_warning(mfilename,'Solver exceeded the wall-clock time budget (options.timeout=%gs).\n', options.timeout);
end
end