solver_mam_basic_fj.m

function [QN,UN,RN,TN,CN,XN,totiter] = solver_mam_basic_fj(sn, options)
% [QN,UN,RN,TN,CN,XN,TOTITER] = SOLVER_MAM_BASIC_FJ(SN, OPTIONS)
%
% MAM decomposition solver with fork-join synchronization via mmap_max.
% Follows the solver_mam.m (dec.mmap) iteration pattern with FJ-aware
% traffic analysis using solver_mam_traffic_fj.
%
% Supports open networks with general FJ topologies and heterogeneous
% service distributions at parallel queues.
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
config = options.config;
if ~isfield(config, 'fj_sync_q_len')
    config.fj_sync_q_len = 2;
end
 
PH = sn.proc;
I = sn.nnodes;
M = sn.nstations;
K = sn.nclasses;
C = sn.nchains;
N = sn.njobs';
V = cellsum(sn.visits);
S = 1./sn.rates;
 
QN = zeros(M,K);
UN = zeros(M,K);
RN = zeros(M,K);
TN = zeros(M,K);
CN = zeros(1,K);
XN = zeros(1,K);
 
% Build FJ synchronization map
fjSyncMap = sn_build_fj_sync_map(sn);
 
% Determine per-chain arrival rates
lambda = zeros(1,K);
for c=1:C
    inchain = sn.inchain{c};
    lambdas_inchain = sn.rates(sn.refstat(inchain(1)),inchain);
    lambdas_inchain = lambdas_inchain(isfinite(lambdas_inchain));
    lambda(inchain) = sum(lambdas_inchain);
end
 
% Prepare PH service distributions
pie = {};
D0 = {};
for ist=1:M
    switch sn.sched(ist)
        case SchedStrategy.EXT
            TN(ist,:) = sn.rates(ist,:);
            TN(ist,isnan(TN(ist,:))) = 0;
        case {SchedStrategy.FCFS, SchedStrategy.HOL, SchedStrategy.FCFSPRPRIO}
            for k=1:K
                % divide service time by number of servers and put
                % later a surrogate delay server in tandem to compensate
                PH{ist}{k} = map_scale(PH{ist}{k}, map_mean(PH{ist}{k})/sn.nservers(ist));
                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
        case SchedStrategy.INF
            for k=1:K
                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
        case SchedStrategy.PS
            for k=1:K
                PH{ist}{k} = map_scale(PH{ist}{k}, map_mean(PH{ist}{k})/sn.nservers(ist));
                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
 
it_max = options.iter_max;
for it=1:it_max
    % Initialize departure processes (node-indexed: DEP{ind,r})
    if it == 1
        DEP = cell(I,K);
        for ind=1:I
            isForkJoin = (sn.nodetype(ind) == NodeType.Fork || sn.nodetype(ind) == NodeType.Join);
            if sn.isstation(ind) && ~isForkJoin
                ist = sn.nodeToStation(ind);
                for r=1:K
                    if V(ist,r) > 0 && lambda(r) > 0
                        DEP{ind,r} = map_scale(PH{ist}{r}, 1 / (lambda(r) * V(ist,r)));
                    else
                        DEP{ind,r} = PH{ist}{r};
                    end
                end
            else
                % Non-station nodes or Fork/Join: pass-through
                % map_exponential takes MEAN interarrival time = 1/rate
                for r=1:K
                    if lambda(r) > 0
                        DEP{ind,r} = map_exponential(1/lambda(r));
                    else
                        DEP{ind,r} = map_exponential(1/GlobalConstants.Immediate);
                    end
                end
            end
        end
    end
 
    % Compute arrival processes with FJ synchronization
    ARV = solver_mam_traffic_fj(sn, DEP, config, fjSyncMap);
 
    QN_1 = QN;
    for ist=1:M
        ind = sn.stationToNode(ist);
        switch sn.nodetype(ind)
            case NodeType.Join
                % Join: throughput from synced arrival rate, no queueing
                if ~isempty(ARV{ind}) && iscell(ARV{ind})
                    arrRates = mmap_lambda(ARV{ind});
                    for k=1:K
                        TN(ist,k) = arrRates(k);
                        UN(ist,k) = 0;
                        QN(ist,k) = 0;
                        RN(ist,k) = 0;
                    end
                end
            case NodeType.Queue
                if ~isempty(ARV{ind}) && iscell(ARV{ind})
                    % Compress arrival process if too large
                    if length(ARV{ind}{1}) > config.space_max
                        if options.verbose
                            line_printf('\nArrival process at node %d is now at %d states. Compressing.', ind, length(ARV{ind}{1}));
                        end
                        ARV{ind} = mmap_compress(ARV{ind});
                    end
 
                    switch sn.sched(ist)
                        case {SchedStrategy.FCFS, SchedStrategy.HOL, SchedStrategy.FCFSPRPRIO}
                            [Qret{1:K}, ~] = MMAPPH1FCFS({ARV{ind}{[1,3:end]}}, {pie{ist}{:}}, {D0{ist,:}}, 'ncMoms', 1, 'ncDistr', 2);
                            for k=1:K
                                QN(ist,k) = sum(Qret{k});
                            end
                            TN(ist,:) = mmap_lambda(ARV{ind});
                        case SchedStrategy.PS
                            TN(ist,:) = mmap_lambda(ARV{ind});
                            for k=1:K
                                UN(ist,k) = TN(ist,k) * S(ist,k);
                            end
                            Uden = min([1-GlobalConstants.FineTol, sum(UN(ist,:))]);
                            for k=1:K
                                QN(ist,k) = UN(ist,k)/(1-Uden);
                            end
                    end
 
                    for k=1:K
                        UN(ist,k) = TN(ist,k) * map_mean(PH{ist}{k});
                        % add number of jobs at the surrogate delay server
                        QN(ist,k) = QN(ist,k) + TN(ist,k)*(map_mean(PH{ist}{k})*sn.nservers(ist)) * (sn.nservers(ist)-1)/sn.nservers(ist);
                        RN(ist,k) = QN(ist,k) ./ TN(ist,k);
                    end
                end
            otherwise
                switch sn.sched(ist)
                    case SchedStrategy.INF
                        if ~isempty(ARV{ind}) && iscell(ARV{ind})
                            TN(ist,:) = mmap_lambda(ARV{ind});
                        end
                        for k=1:K
                            if TN(ist,k) > 0
                                UN(ist,k) = S(ist,k)*TN(ist,k);
                                QN(ist,k) = TN(ist,k)*S(ist,k);
                                RN(ist,k) = S(ist,k);
                            end
                        end
                    case SchedStrategy.EXT
                        % Source: already set TN above
                end
        end
    end
 
    % Check convergence
    if it >= 3 && max(abs(QN(:)-QN_1(:))./(QN_1(:)+GlobalConstants.FineTol)) < options.iter_tol
        break;
    end
 
    % Update departure processes
    for ist=1:M
        ind = sn.stationToNode(ist);
        switch sn.nodetype(ind)
            case NodeType.Queue
                if ~isempty(ARV{ind}) && iscell(ARV{ind})
                    switch sn.sched(ist)
                        case {SchedStrategy.FCFS, SchedStrategy.HOL, SchedStrategy.FCFSPRPRIO}
                            [Ret{1:2*K}] = MMAPPH1FCFS({ARV{ind}{[1,3:end]}}, {pie{ist}{:}}, {D0{ist,:}}, 'stDistrPH');
                            for r=1:K
                                % define a ph for the arrival process of class r
                                A = mmap_hide(ARV{ind}, setdiff(1:K,r));
                                tA = sum(A{2},2);
                                pieA = map_pie(A);
 
                                % define a ph for the service process of class r
                                Srv = PH{ist}{r};
                                pieS = map_pie(Srv);
                                tS = sum(Srv{2},2);
 
                                rho = sum(UN(ist,:));
                                AQ = sum(QN(ist,:));
                                Afull = AQ/max(rho, GlobalConstants.FineTol);
                                pfullonarrival = 1 - (Afull)/(1+Afull);
 
                                A = mmap_scale(A, map_mean(A)-map_mean(Srv));
 
                                zAS = 0*tA*pieS;
                                zSA = 0*tS*pieA;
                                zA = 0*A{2};
 
                                DEP0ir = [Srv{1}, sum(rho)*tS*pieA;
                                    zAS, A{1}];
 
                                DEP1ir = [(1-sum(rho))*Srv{2}, zSA;
                                    tA*pieS, zA];
 
                                DEP{ind,r} = map_normalize({DEP0ir, DEP1ir});
                                if V(ist,r) > 0 && lambda(r) > 0
                                    DEP{ind,r} = map_scale(DEP{ind,r}, 1 / (lambda(r) * V(ist,r)));
                                end
                            end
                        case SchedStrategy.PS
                            % For PS, use simpler departure process update
                            for r=1:K
                                if V(ist,r) > 0 && lambda(r) > 0
                                    DEP{ind,r} = map_scale(PH{ist}{r}, 1 / (lambda(r) * V(ist,r)));
                                end
                            end
                    end
                end
            case NodeType.Join
                % Join: departure = synchronized arrival (pass-through)
                % map_exponential takes MEAN interarrival time = 1/rate
                for r=1:K
                    if TN(ist,r) > 0
                        DEP{ind,r} = map_exponential(1/TN(ist,r));
                    end
                end
        end
        % Fork nodes keep their initial departure processes (pass-through)
    end
end
 
totiter = it;
if options.verbose
    line_printf('\nMAM FJ parametric decomposition completed in %d iterations.', it);
end
 
% Post-processing
CN = sum(RN,1);
QN(isnan(QN)) = 0;
RN(isnan(RN)) = 0;
UN(isnan(UN)) = 0;
TN(isnan(TN)) = 0;
end