afterGlobalEvent.m

function [outglspace, outrate, outprob] =  afterGlobalEvent(sn, ind, glspace, glevent, isSimulation)
% [OUTGLSPACE, OUTRATE, OUTPROB] =  AFTERGLOBALEVENT(QN, IND, GLSPACE, GLEVENT, ISSIMULATION)
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
outspace = []; % state space of ind node
outrate = []; % rate of synchronization
outprob = []; % probability of synchronization
outglspace = glspace; % new global state after synchronization
 
R = sn.nclasses;
%phasessz = sn.phasessz;
%phaseshift = sn.phaseshift;
if sn.nodetype(ind) == NodeType.Transition % same isa(glevent,'ModeEvent')
    % ModeEvent
    isf = sn.nodeToStateful(ind);
    inspace = glspace{isf};
    V = sum(sn.nvars(ind,:));
    % in this case service is always immediate so sum(K)=1
    space_var = inspace(:,(end-V+1):end); % local state variables
    if sn.nodetype(ind) == NodeType.Transition
        fK = sn.nodeparam{ind}.firingphases;
        fKs = [0,cumsum(fK,2)];
        nmodes = sn.nodeparam{ind}.nmodes;
        mode = glevent.active{1}.mode;
        space_buf = inspace(:,1:nmodes); % idle servers count put in buf
        space_srv = inspace(:,(nmodes+1):(nmodes+sum(fK))); % enabled servers' phases
 
        switch glevent.active{1}.event
            case EventType.ENABLE
                enabling_m = sn.nodeparam{ind}.enabling{mode}; % enabling requirement for mode m
                ep_space = zeros(sn.nnodes,R);
                for j=1:length(glevent.passive)
                    ep_ind = glevent.passive{j}.node; % index of enabling place
                    ep_isf = sn.nodeToStateful(ep_ind);
                    K = ones(1,R);
                    Ks = [0,cumsum(K,2)];
                    ep_space_buf = glspace{ep_isf};
                    ep_space_srv = zeros(1,R);
                    ep_space_var = [];
                    [~,ep_space(ep_ind,1:R)] = State.toMarginalAggr(sn,ep_ind, glspace{ep_isf},K,Ks,ep_space_buf,ep_space_srv,ep_space_var);
                end
                m  = mode;
                if any(ep_space < enabling_m)
                    % disable all servers in this mode
                    space_buf(m) = sn.nodeparam{ind}.nmodeservers(m);
                    space_srv((fKs(m)+1):(fKs(m)+fK(m))) = 0;
                    outspace = [space_buf, space_srv, space_var];
                    outrate = GlobalConstants.Immediate;
                    outprob = 1.0;
                else
                    % find enabling degree, i.e., running servers
                    en_degree_m = 1;
                    while all(ep_space >= en_degree_m * enabling_m)
                        en_degree_m = en_degree_m + 1;
                    end
                    en_degree_m = min(en_degree_m - 1, sn.nodeparam{ind}.nmodeservers(m));
                    running_m = sum(space_srv((fKs(m)+1):(fKs(m)+fK(m))));
                    if running_m == en_degree_m
                        % all running as expected, do nothing
                        outspace = [];
                        outrate = [];
                        outprob = 1.0;
                    elseif running_m < en_degree_m
                        % fewer expected, start en_degree_m - running_m
                        space_buf(m) = space_buf(m) - (en_degree_m - running_m);
                        pentry = sn.nodeparam{ind}.firingpie{mode};
                        nadd = en_degree_m - running_m;
                        combs = sortrows(multichoose(fK(m), nadd),'descend');
 
                        for i = 1:size(combs,1)
                            comb = combs(i,:);
                            space_srv_k = space_srv;
                            for k = 1:fK(m)
                                space_srv_k(fKs(m)+k) = space_srv_k(fKs(m)+k) + comb(k);
                            end
                            outspace(end+1,:) = [space_buf, space_srv_k, space_var];
                            outrate(end+1,:) = GlobalConstants.Immediate;
 
                            % compute multinomial coefficient
                            % multinomial probability: n! / (k1! * k2! * ...) * p1^k1 * p2^k2 * ...
                            logprob = factln(nadd);
                            for k = 1:fK(m)
                                if pentry(k)>0
                                    logprob = logprob + comb(k)*log(pentry(k)) - factln(comb(k));                                
                                elseif (pentry(k)==0 && comb(k)==0)
                                    continue
                                else % (pentry(k)==0 && comb(k)>0)
                                    logprob = -Inf; % set zero outprob
                                end
                            end
                            outprob(end+1,:) = exp(logprob);
                        end
                        %outprob = outprob / sum(outprob);
                    else %running_m > en_degree_m
                        % more than expected, stop running_m - en_degree_m
                        % chosen at random
                        diff = en_degree_m - running_m;
                        srv_vec = space_srv((fKs(m)+1):(fKs(m)+fK(m)));
                        n = numel(srv_vec);
 
                        % Generate all combinations of indices to remove
                        idx_combinations = nchoosek(1:n, diff);
 
                        % Generate all derived vectors removing diff
                        for i = 1:size(idx_combinations, 1)
                            idx = true(1, n);
                            idx(idx_combinations(i, :)) = false;
                            outspace(i,:) = [space_buf, srv_vec(idx), space_var];
                            outrate(i,1) = GlobalConstants.Immediate;
                            outprob(i,1) = 1 / size(idx_combinations, 1);
                        end
                    end
                end
                % update new state space
                outglspace{isf} = outspace;
            case EventType.FIRE
                %% Update transition servers
                fK = sn.nodeparam{ind}.firingphases;
                fKs = [0,cumsum(fK,2)];
                mode = glevent.active{1}.mode;
                % find transition server counts
                [~,nim,~,kim] = State.toMarginal(sn,ind,inspace,fK,fKs,space_buf,space_srv,space_var);
                % state of transition mode
                enabling_m = sn.nodeparam{ind}.enabling{mode}; % enabling requirement for mode m
                firing_m = sn.nodeparam{ind}.firing{mode}; % firing requirement for mode m
                % find enabling degree, i.e., running servers
                ep_space = zeros(sn.nnodes,R);
                for j=1:length(glevent.passive)
                    ep_ind = glevent.passive{j}.node; % index of enabling place
                    ep_isf = sn.nodeToStateful(ep_ind);
                    K = ones(1,R);
                    Ks = [0,cumsum(K,2)];
                    ep_space_buf = outglspace{ep_isf}; % this equals glspace at this point
                    ep_space_srv = zeros(1,R);
                    ep_space_var = [];
                    [~,ep_space(ep_ind,1:R)] = State.toMarginalAggr(sn,ep_ind, glspace{ep_isf},K,Ks,ep_space_buf,ep_space_srv,ep_space_var);
                end
                en_degree_m = 1;
                while all(ep_space >= en_degree_m * enabling_m)
                    en_degree_m = en_degree_m + 1;
                end
                %en_degree_m = min(en_degree_m - 1, sn.nodeparam{ind}.nmodeservers(m));
                en_degree_m = min(en_degree_m - 1, nim(mode)); % the actual enabling degree depends on servers actually in execution, though this should coincide
 
                % Create a new state for all possible phase transitions
                % in the exponential case, this will be a single state
                % We will later do the Cartesian product with the new place
                % states
                for k=1:K(mode) % new phase
                    % we do not support MAP firing yet, so we use the total
                    % departure rate from the current phase k in firingproc
                    rate_kd = sum(sn.nodeparam{ind}.firingproc{mode}{2}(k,:)).*kim(:,mode,k);
                    space_buf_kd = space_buf;
                    space_srv_kd = space_srv;
                    % decrease by one the firing server
                    space_srv_kd(fKs(mode)+k) = space_srv_kd(fKs(mode)+k) - 1;
                    % move the server back to the disabled pool
                    space_buf_kd(mode) = space_buf_kd(mode) + 1;
                    space_var_kd = [];
                    outspace = [outspace; space_buf_kd, space_srv_kd, space_var_kd]; %#ok<AGROW>
                    outrate = [outrate; rate_kd]; %#ok<AGROW>
                    outprob = [outprob; 1.0]; %#ok<AGROW>
                end
                outglspace{isf} = outspace;
 
                %% Process ID_PRE events, i.e., consume from all input places
                for j=1:length(glevent.passive)
                    if glevent.passive{j}.event == EventType.PRE
                        ep_ind = glevent.passive{j}.node; % index of enabling place
                        ep_isf = sn.nodeToStateful(ep_ind);
                        ep_space_buf = outglspace{ep_isf}; % this equals glspace at this point
                        ep_ist = sn.nodeToStation(ep_ind);
                        % Places support only for FCFS, LCFS, SIRO
                        consume_set = en_degree_m * enabling_m; % jobs to be consumed
                        switch sn.sched(ep_ist)
                            case SchedStrategy.FCFS
                                for k = find(consume_set)
                                    idx = find(ep_space_buf == k);
                                    ep_space_buf(idx(end - consume_set(k) + 1:end)) = []; % Remove rightmost matches
                                end
                            case SchedStrategy.LCFS
                                for k = find(consume_set)
                                    idx = find(ep_space_buf == k);
                                    ep_space_buf(idx(1:consume_set(k))) = []; % Remove leftmost matches
                                end
                            case SchedStrategy.SIRO
                                ep_space_buf = ep_space_buf - consume_set; % Decrease counts directly
                            otherwise
                                line_error(mfilename,sprintf('Scheduling strategy %s is unsupported at places.', SchedStrategy.toText(sn.sched(ist))));
                        end
                        % update state of enabling place
                        outglspace{ep_isf} = ep_space_buf;
                    end
                end
 
                %% Process ID_POST events, i.e., produce to all output places
                for j=1:length(glevent.passive)
                    if glevent.passive{j}.event == EventType.POST
                        fp_ind = glevent.passive{j}.node; % index of place firing to
                        fp_isf = sn.nodeToStateful(fp_ind);
                        fp_space_buf = outglspace{fp_isf}; % this might have been already modified by the ID_PRE events
                        % Places support only for FCFS, LCFS, SIRO
                        produce_set = firing_m; % jobs to be fired
                        switch sn.sched(ep_ist)
                            case {SchedStrategy.FCFS, SchedStrategy.LCFS}
                                for k = find(produce_set)
                                    fp_space_buf = [repmat(k, 1, find(produce_set == k)), fp_space_buf];
                                end
                            case SchedStrategy.SIRO
                                fp_space_buf = fp_space_buf + produce_set; % Decrease counts directly
                            otherwise
                                line_error(mfilename,sprintf('Scheduling strategy %s is unsupported at places.', SchedStrategy.toText(sn.sched(ist))));
                        end
                        % update state of enabling place
                        outglspace{fp_isf} = fp_space_buf;
                    end
                end
        end
    end
end
 
% TODO: if there are FIRE events with Immediate rate it
% is unclear if the ENABLE event should get priority.
% The order of choice of what fires could change the system behavior.
% On the other hand, this could cause perpetual loops if the enabling is
% applied with priority on the same state with a FIRE?
if isSimulation
    if size(outspace,1) > 1
        tot_rate = sum(outrate);
        cum_rate = cumsum(outrate) / tot_rate;
        firing_ctr = 1 + max([0,find( rand > cum_rate' )]); % select action
        outspace = outspace(firing_ctr,:);
        outrate = sum(outrate);
        outprob = outprob(firing_ctr,:);
    end
end
end