solver_ssa.m

function [pi,SSq,arvRates,depRates,tranSysState,tranSync,sn]=solver_ssa(sn, init_state, options, eventCache)
% [PI,SSQ,ARVRATES,DEPRATES,TRANSYSSTATE,QN]=SOLVER_SSA(QN,OPTIONS)
 
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
% by default the jobs are all initialized in the first valid state
 
% Impatience support checks (mirror SolverCTMC): reneging supports only
% exponential (memoryless) patience; balking supports only QUEUE_LENGTH.
if isfield(sn,'impatienceClass') && ~isempty(sn.impatienceClass)
    badRenege = (sn.impatienceClass==ImpatienceType.RENEGING) & (sn.impatienceType~=ProcessType.EXP);
    if any(badRenege(:))
        line_error(mfilename,'SolverSSA supports only exponential (memoryless) patience for reneging. Use SolverLDES or SolverJMT for phase-type patience.');
    end
end
if isfield(sn,'balkingStrategy') && ~isempty(sn.balkingStrategy)
    badBalk = (sn.balkingStrategy~=0) & (sn.balkingStrategy~=BalkingStrategy.QUEUE_LENGTH);
    if any(badBalk(:))
        line_error(mfilename,'SolverSSA supports only QUEUE_LENGTH balking. Use SolverLDES or SolverJMT for wait-time-based balking.');
    end
end
if isfield(sn,'retrialProc') && ~isempty(sn.retrialProc)
    hasRetrial = ~cellfun(@isempty, sn.retrialProc);
    if any(hasRetrial(:))
        if any(hasRetrial(:) & (sn.retrialType(:)~=ProcessType.EXP))
            line_error(mfilename,'SolverSSA supports only exponential (memoryless) retrial delay. Use SolverLDES or SolverMAM for phase-type retrials.');
        end
        if any(hasRetrial(:) & (sn.retrialMaxAttempts(:)>=0))
            line_error(mfilename,'SolverSSA supports only unlimited retrials (maxAttempts=-1). Use SolverLDES for finite max-attempts.');
        end
        for ii = find(any(hasRetrial,2))'
            served = 0;
            for rr = 1:sn.nclasses
                if ~isempty(sn.proc{ii}{rr}) && ~any(any(isnan(sn.proc{ii}{rr}{1})))
                    served = served + 1;
                end
            end
            if served > 1
                line_error(mfilename,'SolverSSA supports retrial only for single-class stations. Use SolverLDES for multi-class retrial.');
            end
        end
    end
end
 
if ~isfield(options,'seed')
    options.seed = 23000;
end
% Handle parallel computing toolbox gracefully - get worker index
if isMATLABReleaseOlderThan("R2022b")
    % Use labindex for older MATLAB versions
    try
        if ~isempty(getCurrentTask())
            lab_idx = labindex();  %#ok<DLABINDEX>
        else
            lab_idx = 1;
        end
    catch
        line_warning(mfilename,'Parallel Computing Toolbox not available or not running in parallel mode. Using labindex = 1.');
        lab_idx = 1;
    end
else
    % Use spmdIndex for R2022b and newer (labindex is deprecated)
    try
        lab_idx = spmdIndex;
        if isempty(lab_idx) || lab_idx == 0
            lab_idx = 1;
        end
    catch
        line_warning(mfilename,'Parallel Computing Toolbox not available or not running in parallel mode. Using labindex = 1.');
        lab_idx = 1;
    end
end
Solver.resetRandomGeneratorSeed(options.seed + lab_idx - 1);
 
%% generate local state spaces
%nstations = sn.nstations;
nstateful = sn.nstateful;
%init_nserver = sn.nservers; % restore Inf at delay nodes
R = sn.nclasses;
N = sn.njobs';
nnodes = sn.nnodes;
sync = sn.sync;
gsync = sn.gsync;
 
line_debug('SSA solver starting: nstateful=%d, nclasses=%d, njobs=%s, samples=%d', nstateful, R, mat2str(N), options.samples);
csmask = sn.csmask;
 
cutoff = options.cutoff;
if isscalar(cutoff)
    cutoff = cutoff * ones(sn.nstations, sn.nclasses);
end
 
%%
Np = N';
capacityc = zeros(sn.nnodes, sn.nclasses);
original_classcap = sn.classcap; % preserve original classcap for class switching scenarios
for ind=1:sn.nnodes
    if sn.isstation(ind) % place jobs across stations
        ist = sn.nodeToStation(ind);
        %isf = sn.nodeToStateful(ind);
        for r=1:sn.nclasses %cut-off open classes to finite capacity
            c = find(sn.chains(:,r));
            % Check if visits is 0, but also preserve capacity for classes that can
            % receive jobs via class switching (indicated by non-zero original classcap)
            if isfield(sn,'fjclassmap') && ~isempty(sn.fjclassmap) && length(sn.fjclassmap) >= r && sn.fjclassmap(r) > 0
                % FJ auxiliary sibling class: the per-tag multiplicity
                % (tasksPerLink) can exceed the chain job population, so the
                % adapter-set classcap is the correct capacity bound
                capacityc(ind,r) = original_classcap(ist,r);
            elseif ~isempty(sn.visits{c}) && sn.visits{c}(ist,r) == 0 && original_classcap(ist,r) == 0
                capacityc(ind,r) = 0;
            elseif ~isempty(sn.proc) && ~isempty(sn.proc{ist}{r}) && any(any(isnan(sn.proc{ist}{r}{1}))) && sn.nodetype(ind) ~= NodeType.Place % disabled (but not Place nodes)
                capacityc(ind,r) = 0;
            else
                if isinf(N(r))
                    capacityc(ind,r) =  min(cutoff(ist,r), sn.classcap(ist,r));
                else
                    % closed classes: enumerate up to the chain population, but never
                    % beyond the class capacity at this station (finite-buffer stations)
                    capacityc(ind,r) =  min(sum(sn.njobs(sn.chains(c,:))), sn.classcap(ist,r));
                end
            end
        end
        % never raise the station capacity above its configured total capacity
        capacity_sum = min(sum(capacityc(ind,:)), sn.cap(ist));
        if sn.sched(ist) == SchedStrategy.PAS
            % Pass-and-swap stations are a single order-independent buffer of
            % total size cap, not a sum of per-class buffers. Preserve the
            % original total capacity, which is also the physical width used to
            % size the ordered-list state (State.fromMarginal Wpas = cap).
            capacity_sum = sn.cap(ist);
        end
        if isinf(sn.nservers(ist))
            sn.nservers(ist) = capacity_sum;
        end
        sn.cap(ist,:) = capacity_sum;
        sn.classcap(ist,:) = capacityc(ind,:);
    end
end
% Signal (G-network negative-customer) classes never occupy a station; cap
% their per-station capacity at 0 (except at the Source/EXT node) so the state
% space does not enumerate unreachable states holding signal jobs.
if isfield(sn,'issignal') && ~isempty(sn.issignal) && any(sn.issignal)
    for ii = 1:sn.nstations
        if sn.sched(ii) ~= SchedStrategy.EXT
            sn.classcap(ii, sn.issignal(:)') = 0;
        end
    end
end
 
% Heterogeneous servers (single-class, ORDER policy) -> load-dependent rate
% mu(n) = sum of the first min(n,c) server rates (mirrors SolverCTMC). Multi-
% class heterogeneous servers are rejected (need per-server-type state).
for ind = 1:sn.nnodes
    if sn.isstation(ind) && isfield(sn,'nodeparam') && numel(sn.nodeparam) >= ind ...
            && ~isempty(sn.nodeparam{ind}) && isstruct(sn.nodeparam{ind}) ...
            && isfield(sn.nodeparam{ind},'nservertypes') && sn.nodeparam{ind}.nservertypes > 0
        ist = sn.nodeToStation(ind);
        np = sn.nodeparam{ind};
        served = [];
        for r = 1:sn.nclasses
            if ~isempty(sn.proc{ist}{r}) && ~any(any(isnan(sn.proc{ist}{r}{1}))) && sn.rates(ist,r) > 0
                served(end+1) = r; %#ok<AGROW>
            end
        end
        if numel(served) > 1
            line_error(mfilename,'SolverSSA supports heterogeneous servers only for single-class stations. Use SolverJMT or SolverLDES for multi-class heterogeneous servers.');
        end
        if numel(served) == 1
            r = served;
            srvrates = [];
            for t = 1:np.nservertypes
                if np.servercompat(t,r) && np.heterorates(t,r) > 0
                    srvrates = [srvrates, repmat(np.heterorates(t,r), 1, np.serverspertype(t))]; %#ok<AGROW>
                end
            end
            c = numel(srvrates);
            mu_base = sn.rates(ist,r);
            if c > 0 && mu_base > 0
                if isempty(sn.lldscaling)
                    sn.lldscaling = ones(sn.nstations, max([c, sum(sn.njobs(isfinite(sn.njobs))), 1]));
                elseif size(sn.lldscaling,2) < c
                    sn.lldscaling(:, (size(sn.lldscaling,2)+1):c) = repmat(sn.lldscaling(:,end), 1, c-size(sn.lldscaling,2));
                end
                for n = 1:size(sn.lldscaling,2)
                    mun = sum(srvrates(1:min(n,c)));
                    sn.lldscaling(ist,n) = mun / (mu_base * min(n,c));
                end
            end
        end
    end
end
 
% Finite Capacity Region: precompute per-region member stations and caps. The
% simulation blocks (holds upstream) any arrival that would push the aggregate
% region population past a cap; see the FCR gate in the event loop below. This
% is the SSA counterpart of the state-space filter in solver_ctmc.m.
fcrOn = isfield(sn,'nregions') && sn.nregions > 0;
if fcrOn
    Kfcr = sn.nclasses;
    fcrMembers   = cell(sn.nregions,1);
    fcrMemberMask= cell(sn.nregions,1);
    fcrClassCap  = cell(sn.nregions,1);
    fcrGlobalCap = inf(sn.nregions,1);
    fcrMemCap    = inf(sn.nregions,1);
    fcrSz        = cell(sn.nregions,1);
    fcrA         = cell(sn.nregions,1);
    fcrb         = cell(sn.nregions,1);
    for f = 1:sn.nregions
        Rmat = sn.region{f};                          % M x (K+1)
        % membership: any job-count cap OR the region memory budget set on the
        % station row (a memory-only region has all job-count entries at -1)
        memvecFCR = -ones(sn.nstations,1);
        if isfield(sn,'regionmaxmem') && numel(sn.regionmaxmem) >= f && ~isempty(sn.regionmaxmem{f})
            memvecFCR = sn.regionmaxmem{f}(:);
        end
        mask = (any(Rmat ~= -1, 2) | memvecFCR ~= -1)';
        fcrMemberMask{f} = mask;
        fcrMembers{f} = find(mask);
        ccap = inf(1,Kfcr);
        for r = 1:Kfcr
            cv = Rmat(fcrMembers{f}, r); cv = cv(cv ~= -1);
            if ~isempty(cv); ccap(r) = min(cv); end
        end
        fcrClassCap{f} = ccap;
        gv = Rmat(fcrMembers{f}, Kfcr+1); gv = gv(gv ~= -1);
        if ~isempty(gv); fcrGlobalCap(f) = min(gv); end
        if isfield(sn,'regionmaxmem') && numel(sn.regionmaxmem) >= f && ~isempty(sn.regionmaxmem{f})
            mv = sn.regionmaxmem{f}(fcrMembers{f}); mv = mv(mv ~= -1);
            if ~isempty(mv); fcrMemCap(f) = min(mv); end
        end
        fcrSz{f} = sn.regionsz(f,:);
        if isfield(sn,'regionlincon') && size(sn.regionlincon,1) >= f && ~isempty(sn.regionlincon{f,1})
            fcrA{f} = sn.regionlincon{f,1};
            fcrb{f} = sn.regionlincon{f,2};
        end
    end
    % WAITQ rule per (region, class): a refused entry parks in a per-region
    % FIFO of (class, destination) tokens (JMT waiting-queue semantics) and
    % re-enters head-of-line as capacity frees; DROP keeps censoring.
    fcrRule = false(sn.nregions, Kfcr);
    if isfield(sn,'regionrule') && ~isempty(sn.regionrule)
        for f = 1:sn.nregions
            for r = 1:Kfcr
                fcrRule(f,r) = sn.regionrule(f,r) ~= DropStrategy.DROP;
            end
        end
    end
    fcrBuf = cell(sn.nregions,1);
    for f = 1:sn.nregions
        fcrBuf{f} = zeros(1,0);
    end
end
 
%%
if any(isinf(Np))
    Np(isinf(Np)) = 0;
end
 
init_state_hashed = ones(1,nstateful); % pick the first state in init_state{i}
 
%%
arvRatesSamples = zeros(options.samples,nstateful,R);
depRatesSamples = zeros(options.samples,nstateful,R);
A = length(sync);
G = length(gsync);
samples_collected = 1;
nir = {};
% fill stateCell with initial states
cur_state = cell(nstateful,1); % cell array with current stateful node states
for ind=1:sn.nnodes
    if sn.isstateful(ind)
        isf = sn.nodeToStateful(ind);
        cur_state{isf} = init_state{isf}(init_state_hashed(isf),:);
        if sn.isstation(ind)
            ist = sn.nodeToStation(ind);
            [~,nir{ist}] = State.toMarginal(sn, ind, init_state{isf}(init_state_hashed(isf),:));
            nir{ist} = nir{ist}(:);
        end
    end
end
cur_state_1 = cur_state;
% generate state vector
state = cell2mat(cur_state');
% create function to determine lengths of stateful node states
statelen = cellfun(@length, cur_state);
% data structures to save transient information - pre-allocate for all samples
nSamples = options.samples;
tranSync = zeros(nSamples,1);
tranState = zeros(1+length(state), nSamples);
tranState(1:(1+length(state)),1) = [0, state]';
SSq = zeros(length(cell2mat(nir')), nSamples);
SSq(:,1) = cell2mat(nir');
local = sn.nnodes+1;
last_node_a = 0; % active in the last occurred synchronization
last_node_p = 0; % passive in the last occurred synchronization
for act=1:A
    node_a{act} = sync{act}.active{1}.node;
    node_p{act} = sync{act}.passive{1}.node;
    class_a{act} = sync{act}.active{1}.class;
    class_p{act} = sync{act}.passive{1}.class;
    event_a{act} = sync{act}.active{1}.event;
    event_p{act} = sync{act}.passive{1}.event;
    outprob_a{act} = [];
    outprob_p{act} = [];
    % Immediate-feedback self-loop detection (loop-invariant): a departure
    % that routes back to the same FCFS-family station and arrives as a class
    % flagged in sn.immfeed must hold the server instead of re-queueing behind
    % the waiting jobs. The departure half is then computed with noPromote so
    % the vacated server is re-entered by the fed-back job (see
    % State.afterEventStation).
    immfeed_selfloop{act} = false;
    if event_a{act}==EventType.DEP && node_p{act}==node_a{act} ...
            && node_a{act}>=1 && node_a{act}<=sn.nnodes && sn.isstation(node_a{act}) ...
            && isfield(sn,'immfeed') && ~isempty(sn.immfeed)
        istA_if = sn.nodeToStation(node_a{act});
        if istA_if>=1 && istA_if<=size(sn.immfeed,1) ...
                && class_p{act}>=1 && class_p{act}<=size(sn.immfeed,2) ...
                && sn.immfeed(istA_if, class_p{act})
            immfeed_selfloop{act} = true;
        end
    end
end
enabled_next_states = cell(1,A);
 
%% Start main simulation loop
isSimulation = true; % allow state vector to grow, e.g. for FCFS buffers
% Precompute the loop-invariant afterEvent context. This MUST follow the
% preamble above that rewrites sn.nservers/cap/classcap (Inf servers and
% open-class cutoffs), so that blocking and rate computations match the
% per-call derivation it replaces.
aectx = State.afterEventInit(sn);
samples_collected = 1;
cur_time = 0;
use_inline = true; % true = stable version, false = dev version
 
try
    while samples_collected < options.samples && cur_time <= options.timespan(2) && ~lineTimeoutExceeded(options)
        %% This section corresponds to solver_ssa_findenabled in Java
        %% Inlined for performance reasons
        if use_inline
            enabled_sync = []; % row is action label, col1=rate, col2=new state
            enabled_rates = [];
            enabled_fcr = zeros(0,4); % [region class dest isSwitch] FCR marker per transition
            ctr = 1;
            A = length(sync);
            G = length(gsync);
            % FCR: current aggregate per-class population of each region, used by
            % the arrival gate below to block entries that would exceed a cap.
            if fcrOn
                xcurFCR = cell(sn.nregions,1);
                for f = 1:sn.nregions
                    xf = zeros(1,sn.nclasses);
                    for i = fcrMembers{f}
                        ind_i = sn.stationToNode(i);
                        isf_i = sn.stationToStateful(i);
                        [~, nir_i] = State.toMarginalAggr(sn, ind_i, cur_state{isf_i});
                        xf = xf + nir_i(:)';
                    end
                    xcurFCR{f} = xf;
                end
            end
            for act=1:A
                isf_a = sn.nodeToStateful(node_a{act});
 
                enabled_next_states{act} = cur_state;
                update_cond_a = true;
                if update_cond_a
                    [enabled_next_states{act}{isf_a}, rate_a{act}, outprob_a{act}, eventCache] =  State.afterEvent(sn, node_a{act}, cur_state{isf_a}, event_a{act}, class_a{act}, isSimulation, eventCache, aectx, immfeed_selfloop{act});
                end
 
                if isempty(enabled_next_states{act}{isf_a}) || isempty(rate_a{act})
                    continue
                end
 
                for ia=1:size(enabled_next_states{act}{isf_a},1) % for all possible new states, check if they are enabled
                    % if the transition cannot occur
                    if isnan(rate_a{act}(ia)) || rate_a{act}(ia) == 0 % handles degenerate rate values
                        % set the transition with a zero rate so that it is
                        % never selected
                        rate_a{act}(ia) = 1e-38; % ~ zero in 32-bit precision
                    end
 
                    if enabled_next_states{act}{isf_a}(ia,:) == -1 % hash not found
                        continue
                    end
                    update_cond_p = true; %samples_collected == 1 || ((node_p{act} == last_node_a || node_p{act} == last_node_p)) || isempty(outprob_a{act}) || isempty(outprob_p{act});
 
                    if rate_a{act}(ia)>0
                        if node_p{act} ~= local
                            if node_p{act} == node_a{act} %self-loop, active and passive are the same
                                isf_p = isf_a;
                                if update_cond_p
                                    [enabled_next_states{act}{isf_p}, ~, outprob_p{act}, eventCache] =  State.afterEvent(sn, node_p{act}, enabled_next_states{act}{isf_p}, event_p{act}, class_p{act}, isSimulation, eventCache, aectx);
                                end
                            else % departure
                                isf_p = sn.nodeToStateful(node_p{act});
                                if update_cond_p
                                    [enabled_next_states{act}{isf_p}, ~, outprob_p{act}, eventCache] =  State.afterEvent(sn, node_p{act}, enabled_next_states{act}{isf_p}, event_p{act}, class_p{act}, isSimulation, eventCache, aectx);
                                end
                            end
                            if ~isempty(enabled_next_states{act}{isf_p})
                                if sn.isstatedep(node_a{act},3)
                                    prob_sync_p{act} = sync{act}.passive{1}.prob(cur_state, enabled_next_states{act}); %state-dependent
                                else
                                    prob_sync_p{act} = sync{act}.passive{1}.prob;
                                end
                            else
                                prob_sync_p{act} = 0;
                            end
                        end
                        if ~isempty(enabled_next_states{act}{isf_a})
                            if node_p{act} == local
                                prob_sync_p{act} = 1;
                            end
                            if ~isnan(rate_a{act})
                                if all(~cellfun(@isempty,enabled_next_states{act}))
                                    % FCR gate: an arrival that would push a region's
                                    % aggregate population past a cap is censored for
                                    % DROP classes and parked in the region FIFO for
                                    % WAITQ classes (active part applied only). A
                                    % class-switching hop between two members of the
                                    % same region is an exit followed by a gated
                                    % re-entry (JMT ClassSwitch-outside-region parity).
                                    blockFCR = false;
                                    fcrMark = [0 0 0 0]; % [region class dest isSwitch]
                                    if fcrOn && node_p{act} ~= local && node_p{act} <= sn.nnodes
                                        jp = sn.nodeToStation(node_p{act});
                                        if jp > 0
                                            ja = sn.nodeToStation(node_a{act});
                                            cc = class_p{act};
                                            for f = 1:sn.nregions
                                                mmask = fcrMemberMask{f};
                                                if mmask(jp) && (ja <= 0 || ja > numel(mmask) || ~mmask(ja))
                                                    xn = xcurFCR{f}; xn(cc) = xn(cc) + 1;
                                                    if xn(cc) > fcrClassCap{f}(cc) || sum(xn) > fcrGlobalCap(f) ...
                                                            || (xn * fcrSz{f}(:) > fcrMemCap(f)) ...
                                                            || (~isempty(fcrA{f}) && any(fcrA{f} * xn(:) > fcrb{f}(:)))
                                                        if fcrRule(f,cc)
                                                            fcrMark = [f cc node_p{act} 0]; % park in FIFO
                                                        else
                                                            fcrMark = [f cc node_p{act} 2]; % DROP: destroyed
                                                        end
                                                        break;
                                                    end
                                                elseif mmask(jp) && ja > 0 && ja <= numel(mmask) && mmask(ja) ...
                                                        && cc ~= class_a{act}
                                                    if fcrRule(f,cc)
                                                        fcrMark = [f cc node_p{act} 1]; % exit + gated re-entry
                                                    else
                                                        fcrMark = [f cc node_p{act} 3]; % exit + gated re-entry, DROP on refusal
                                                    end
                                                    break;
                                                end
                                            end
                                        end
                                    end
                                    if fcrMark(1) > 0
                                        % suppress the passive application: the job
                                        % leaves the upstream node and its entry is
                                        % resolved at application time
                                        enabled_next_states{act}{isf_p} = cur_state{isf_p};
                                    end
                                    if event_a{act} == EventType.DEP && ~blockFCR
                                        node_a_sf{act} = isf_a;
                                        node_p_sf{act} = isf_p;
                                        depRatesSamples(samples_collected,node_a_sf{act},class_a{act}) = depRatesSamples(samples_collected,node_a_sf{act},class_a{act}) + outprob_a{act} * outprob_p{act} * rate_a{act}(ia) * prob_sync_p{act};
                                        arvRatesSamples(samples_collected,node_p_sf{act},class_p{act}) = arvRatesSamples(samples_collected,node_p_sf{act},class_p{act}) + outprob_a{act} * outprob_p{act} * rate_a{act}(ia) * prob_sync_p{act};
                                    end
                                    % simulate also self-loops as we need to log them
                                    %if any(~cellfun(@isequal,new_state{act},cur_state))
                                    if node_p{act} < local && ~sn.csmask(class_a{act}, class_p{act}) && sn.nodetype(node_p{act})~=NodeType.Source && (rate_a{act}(ia) * prob_sync_p{act} >0)
                                        line_error(mfilename,sprintf('Error: state-dependent routing at node %d (%s) violates the class switching mask (node %d -> node %d, class %d -> class %d).', node_a{act}, sn.nodenames{node_a{act}}, node_a{act}, node_p{act}, class_a{act}, class_p{act}));
                                    end
                                    if ~blockFCR
                                        enabled_rates(ctr) = rate_a{act}(ia) * prob_sync_p{act};
                                        enabled_sync(ctr) = act;
                                        enabled_fcr(ctr,:) = fcrMark;
                                        ctr = ctr + 1;
                                    end
                                end
                            end
                        end
                    end
                end
            end
            gctr_start = ctr;
 
            for gact=1:G % event at node ind with global side-effects
                gind = gsync{gact}.active{1}.node; % node index for global event
                [enabled_next_states{A+gact}, outrate, outprob] = State.afterGlobalEvent(sn, gind, cur_state, gsync{gact}, isSimulation);
                for ia=find(outrate .* outprob)
                    enabled_rates(ctr) = outrate(ia) * outprob(ia);
                    enabled_sync(ctr) = A+gact;
                    ctr = ctr + 1;
 
                    % Record departure/arrival rates for FIRE events at Places
                    if gsync{gact}.active{1}.event == EventType.FIRE
                        mode = gsync{gact}.active{1}.mode;
                        % Get enabling/firing conditions to determine affected classes
                        enabling_m = sn.nodeparam{gind}.enabling{mode};
                        firing_m = sn.nodeparam{gind}.firing{mode};
 
                        for j=1:length(gsync{gact}.passive)
                            pev = gsync{gact}.passive{j};
                            % Decode linear index to (node, class) - pev.node is a linear index from find() on enabling/firing matrix
                            [pev_node, pev_class] = ind2sub([sn.nnodes, R], pev.node);
                            if pev.event == EventType.PRE
                                % Departure from input Place (consuming tokens)
                                if pev_node <= length(sn.nodeToStateful) && ~isnan(sn.nodeToStateful(pev_node)) && sn.nodeToStateful(pev_node) > 0
                                    ep_isf = sn.nodeToStateful(pev_node);
                                    % Record departures for the specific class from this PRE event
                                    depRatesSamples(samples_collected, ep_isf, pev_class) = ...
                                        depRatesSamples(samples_collected, ep_isf, pev_class) + outrate(ia) * outprob(ia);
                                end
                            elseif pev.event == EventType.POST
                                % Arrival at output Place (producing tokens)
                                if pev_node <= length(sn.nodeToStateful) && ~isnan(sn.nodeToStateful(pev_node)) && sn.nodeToStateful(pev_node) > 0
                                    fp_isf = sn.nodeToStateful(pev_node);
                                    % Record arrivals for the specific class from this POST event
                                    arvRatesSamples(samples_collected, fp_isf, pev_class) = ...
                                        arvRatesSamples(samples_collected, fp_isf, pev_class) + outrate(ia) * outprob(ia);
                                end
                            end
                        end
                    end
                end
            end
 
            % fork firing synchronizations (native fork-join support): the
            % firing consumes the parent held at the Fork and emits one
            % sibling per branch atomically; indices beyond gctr_start take
            % the global-event path in the bookkeeping below
            FJ = 0;
            if isfield(sn,'fjsync') && ~isempty(sn.fjsync)
                FJ = length(sn.fjsync);
            end
            for fjact=1:FJ
                [fjStates, fjrate, fjprob] = State.afterFJEvent(sn, sn.fjsync{fjact}, cur_state, isSimulation);
                if ~isempty(fjStates)
                    enabled_next_states{A+G+fjact} = fjStates{1};
                    enabled_rates(ctr) = fjrate(1) * fjprob(1);
                    enabled_sync(ctr) = A+G+fjact;
                    ctr = ctr + 1;
                    fjentry = sn.fjsync{fjact};
                    isf_fork = sn.nodeToStateful(fjentry.fork);
                    depRatesSamples(samples_collected, isf_fork, fjentry.class) = ...
                        depRatesSamples(samples_collected, isf_fork, fjentry.class) + fjrate(1) * fjprob(1);
                    for b=1:length(fjentry.branchheads)
                        isf_bh = sn.nodeToStateful(fjentry.branchheads(b));
                        arvRatesSamples(samples_collected, isf_bh, fjentry.auxclasses(b)) = ...
                            arvRatesSamples(samples_collected, isf_bh, fjentry.auxclasses(b)) + fjrate(1) * fjprob(1);
                    end
                end
            end
        else
            [enabled_next_states,enabled_rates,enabled_sync,gctr_start,depRatesSamples,arvRatesSamples,outprob_a,outprob_p,rate_a, eventCache] = solver_ssa_findenabled(sn,node_a,enabled_next_states,cur_state,outprob_a,event_a,class_a,isSimulation,node_p,local,outprob_p,event_p,class_p,sync,gsync,depRatesSamples,samples_collected,arvRatesSamples,last_node_a,last_node_p,eventCache);
        end
        %% Gillespie direct method
        tot_rate = sum(enabled_rates);
        cum_rate = cumsum(enabled_rates) / tot_rate;
        selected_transition = 1 + max([0,find( rand > cum_rate )]); % select action
 
        % Update record of last active/passive pair
        if isempty(enabled_sync)
            line_error(mfilename,'SSA simulation entered a deadlock before collecting all samples, no synchronization is enabled.');
        end
        if selected_transition < gctr_start
            % regular event pair
            last_node_a = node_a{enabled_sync(selected_transition)};
            last_node_p = node_p{enabled_sync(selected_transition)};
        else % global event
            last_node_a = NaN;
            last_node_p = NaN;
        end
 
        %% Update paddings
        % This part is needed to ensure that when the state vector grows the
        % padding of zero is done on the left (e.g., for FCFS buffers)
        for ind=1:sn.nnodes
            if sn.isstation(ind)
                isf = sn.nodeToStateful(ind);
                deltalen = length(cur_state{isf}) - statelen(isf);
                if deltalen>0
                    statelen(isf) = length(cur_state{isf});
                    % here do padding
                    if ind==1
                        shift = 0;
                    else
                        shift = sum(statelen(1:isf-1));
                    end
                    pad = zeros(deltalen, size(tranState,2));
                    tranState = [tranState(1:(shift+1), :); pad ; tranState((shift+1+deltalen):end, :)];
                end
            end
        end
 
        %% Simulate the time increment
        state = cell2mat(cur_state');
        dt = -(log(rand)/tot_rate);
        cur_time = cur_time + dt;
 
        %% Save simulation output data
        tranState(1:(1+length(state)),samples_collected) = [dt, state]';
        tranSync(samples_collected,1) = enabled_sync(selected_transition);
        for ind=1:sn.nnodes
            if sn.isstation(ind)
                isf = sn.nodeToStateful(ind);
                ist = sn.nodeToStation(ind);
                [~,nir{ist}] = State.toMarginal(sn, ind, cur_state{isf});
                nir{ist}=nir{ist}(:);
            end
        end
        SSq(:,samples_collected) = cell2mat(nir');
 
        %% Update current state and sample counter
        cur_state_1 = cur_state;
        cur_state = enabled_next_states{enabled_sync(selected_transition)};
 
        % KCHOICES with memory: record the destination just chosen so the
        % next dispatch consults it. Only applies to Router-typed active
        % nodes that have KCHOICES routing with withMemory=true. This is the
        % MATLAB-side mirror of LDES's kchoicesLastSelected tracking.
        if selected_transition < gctr_start
            actSync = enabled_sync(selected_transition);
            if ~isempty(node_a{actSync}) && ~isempty(node_p{actSync}) ...
                    && node_a{actSync} >= 1 && node_a{actSync} <= sn.nnodes
                aNode = node_a{actSync};
                aClass = class_a{actSync};
                if sn.nodetype(aNode) == NodeType.Router && aClass >= 1 ...
                        && sn.routing(aNode, aClass) == RoutingStrategy.KCHOICES
                    np_a = [];
                    if iscell(sn.nodeparam) && aNode <= numel(sn.nodeparam) ...
                            && iscell(sn.nodeparam{aNode}) ...
                            && aClass <= numel(sn.nodeparam{aNode})
                        np_a = sn.nodeparam{aNode}{aClass};
                    end
                    if ~isempty(np_a) && isfield(np_a, 'withMemory') && np_a.withMemory
                        isf_a = sn.nodeToStateful(aNode);
                        Rcl = sn.nclasses;
                        if ~isempty(cur_state{isf_a}) && length(cur_state{isf_a}) >= Rcl
                            cur_state{isf_a}(end - Rcl + aClass) = node_p{actSync};
                        end
                    end
                end
            end
        end
 
        %% FCR WAITQ bookkeeping: park blocked entries in the region FIFO,
        % release FIFO heads freed by the applied transition (strict FIFO,
        % head-of-line), and resolve pending class-switch re-entries
        if fcrOn && use_inline
            pk = [0 0 0 0];
            if selected_transition <= size(enabled_fcr,1)
                pk = enabled_fcr(selected_transition,:);
            end
            if pk(1) > 0 && pk(4) == 0
                % blocked entry: park the (class, destination) token
                fcrBuf{pk(1)}(end+1) = (pk(3)-1)*R + pk(2);
            end
            % pk(4)==2: DROP, the refused job was destroyed (active part only)
            % release cascade over all regions
            [cur_state, fcrBuf, eventCache] = fcr_release(sn, cur_state, fcrBuf, ...
                fcrMembers, fcrClassCap, fcrGlobalCap, fcrMemCap, fcrSz, fcrA, fcrb, ...
                isSimulation, eventCache, aectx);
            if pk(1) > 0 && (pk(4) == 1 || pk(4) == 3)
                % class-switching hop: gate the re-entry after the cascade
                f_ = pk(1); cls_ = pk(2); dest_ = pk(3);
                xf_ = fcr_regionpop(sn, cur_state, fcrMembers{f_});
                xn_ = xf_; xn_(cls_) = xn_(cls_) + 1;
                admitted = false;
                if ~fcr_violates(xn_, fcrClassCap{f_}, fcrGlobalCap(f_), fcrMemCap(f_), fcrSz{f_}, fcrA{f_}, fcrb{f_})
                    isf_d = sn.nodeToStateful(dest_);
                    [ns_, ~, ~, eventCache] = State.afterEvent(sn, dest_, cur_state{isf_d}, EventType.ARV, cls_, isSimulation, eventCache, aectx);
                    if ~isempty(ns_)
                        cur_state{isf_d} = ns_(1,:);
                        admitted = true;
                    end
                end
                if ~admitted && pk(4) == 1
                    fcrBuf{f_}(end+1) = (dest_-1)*R + cls_;
                end
                % pk(4)==3 refused: DROP, the switching job is destroyed
            end
        end
 
        samples_collected = samples_collected + 1;
 
        %% Print progress
        print_progress(options,samples_collected);
    end
    % Print newline after progress counter
    if options.verbose
        line_printf('\n');
    end
catch ME
    getReport(ME)
end
 
% Trim pre-allocated arrays to actual number of samples collected
samples_collected = samples_collected - 1;  % Adjust for the increment at end of loop
tranState = tranState(:, 1:samples_collected);
tranSync = tranSync(1:samples_collected, :);
SSq = SSq(:, 1:samples_collected);
 
% Warmup discard: drop the first floor(warmupfrac * samples) samples so that
% steady-state averages are not biased by transient behaviour. Controlled by
% options.config.warmupfrac (default 0.0 = no discard).
warmupfrac = 0.0;
if isfield(options, 'config') && isfield(options.config, 'warmupfrac')
    warmupfrac = max(0.0, min(0.99, options.config.warmupfrac));
end
if warmupfrac > 0 && samples_collected > 1
    nDrop = floor(warmupfrac * samples_collected);
    if nDrop > 0 && nDrop < samples_collected
        keep = (nDrop+1):samples_collected;
        tranState = tranState(:, keep);
        tranSync = tranSync(keep, :);
        SSq = SSq(:, keep);
        if exist('arvRatesSamples', 'var') && ~isempty(arvRatesSamples) ...
                && size(arvRatesSamples, 1) >= samples_collected
            arvRatesSamples = arvRatesSamples(keep, :, :);
        end
        if exist('depRatesSamples', 'var') && ~isempty(depRatesSamples) ...
                && size(depRatesSamples, 1) >= samples_collected
            depRatesSamples = depRatesSamples(keep, :, :);
        end
        samples_collected = numel(keep);
    end
end
 
tranState = tranState';
 
 
[u,ui,uj] = unique(tranState(:,2:end),'rows');
statesz = cellfun(@length, cur_state_1)';
tranSysState = cell(1,length(cur_state)+1);
tranSysState{1} = cumsum(tranState(:,1));
for j=1:length(statesz)
    tranSysState{1+j} = tranState(:,1+(1+sum(statesz(1:(j-1)))):(1+sum(statesz(1:j))));
end
arvRates = zeros(size(u,1),sn.nstateful,R);
depRates = zeros(size(u,1),sn.nstateful,R);
 
pi = zeros(1,size(u,1));
for s=1:size(u,1)
    pi(s) = sum(tranState(uj==s,1));
end
SSq = SSq(:,ui)'; % we restrict to unique states in the simulation
 
for ind=1:sn.nnodes
    if sn.isstateful(ind)
        isf = sn.nodeToStateful(ind);
        if sn.isstation(ind)
            ist = sn.nodeToStation(ind);
            %K = sn.phasessz(ist,:);
            %Ks = sn.phaseshift(ist,:);
        end
        for s=1:size(u,1)
            for r=1:R
                arvRates(s,isf,r) = arvRatesSamples(ui(s),isf,r); % for each unique state, one (any) sample of the rate is enough here
                depRates(s,isf,r) = depRatesSamples(ui(s),isf,r); % for each unique state, one (any) sample of the rate is enough here
            end
        end
    end
end
pi = pi/sum(pi);
%sn.nservers = init_nserver; % restore Inf at delay nodes
end
 
function print_progress(options,samples_collected)
if options.verbose && ~batchStartupOptionUsed
    if samples_collected == 1e2
        line_printf(sprintf('\nSSA samples: %6d',samples_collected));
    elseif options.verbose == 2
        if samples_collected == 0
            line_printf(sprintf('\nSSA samples: %6d',samples_collected));
        else
            line_printf(sprintf('\b\b\b\b\b\b%6d',samples_collected));
        end
    elseif mod(samples_collected,1e2)==0 || options.verbose == 2
        line_printf(sprintf('\b\b\b\b\b\b%6d',samples_collected));
    end
end
end
function x = fcr_regionpop(sn, cur_state, members)
% X=FCR_REGIONPOP(SN,CUR_STATE,MEMBERS) per-class population of a finite
% capacity region given the current state cells
x = zeros(1, sn.nclasses);
for i = members
    ind_i = sn.stationToNode(i);
    isf_i = sn.stationToStateful(i);
    [~, nir_i] = State.toMarginalAggr(sn, ind_i, cur_state{isf_i});
    x = x + nir_i(:)';
end
end
 
function tf = fcr_violates(xn, ccap, gcap, memcap, sz, A, b)
% TF=FCR_VIOLATES(...) true if population vector xn breaks any admission
% constraint of the region
tf = any(xn > ccap) || sum(xn) > gcap || (xn * sz(:) > memcap);
if ~tf && ~isempty(A)
    tf = any(A * xn(:) > b(:));
end
end
 
function [cur_state, fcrBuf, eventCache] = fcr_release(sn, cur_state, fcrBuf, ...
    fcrMembers, fcrClassCap, fcrGlobalCap, fcrMemCap, fcrSz, fcrA, fcrb, ...
    isSimulation, eventCache, aectx)
% FCR_RELEASE strict-FIFO head-of-line release of parked region tokens:
% admit heads while the admission constraints permit, applying the arrival
% to the destination station state
K = sn.nclasses;
progress = true;
while progress
    progress = false;
    for f = 1:length(fcrBuf)
        if isempty(fcrBuf{f})
            continue
        end
        x = fcr_regionpop(sn, cur_state, fcrMembers{f});
        tok = fcrBuf{f}(1);
        dest = floor((tok-1)/K) + 1;
        r = mod(tok-1, K) + 1;
        xn = x; xn(r) = xn(r) + 1;
        if fcr_violates(xn, fcrClassCap{f}, fcrGlobalCap(f), fcrMemCap(f), fcrSz{f}, fcrA{f}, fcrb{f})
            continue % head-of-line: this region's FIFO stays blocked
        end
        isf_d = sn.nodeToStateful(dest);
        [ns, ~, ~, eventCache] = State.afterEvent(sn, dest, cur_state{isf_d}, EventType.ARV, r, isSimulation, eventCache, aectx);
        if isempty(ns)
            continue % destination cannot accept (e.g. station capacity)
        end
        cur_state{isf_d} = ns(1,:);
        fcrBuf{f}(1) = [];
        progress = true;
    end
end
end