updateMetricsDefault.m

function updateMetricsDefault(self, it)
ensemble = self.ensemble;
lqn = self.lqn;
 
% obtain the activity service times
self.servt = zeros(lqn.nidx,1);
self.residt = zeros(lqn.nidx,1);
for r=1:size(self.servt_classes_updmap,1)
    idx = self.servt_classes_updmap(r,1);
    aidx = self.servt_classes_updmap(r,2);
    nodeidx = self.servt_classes_updmap(r,3);
    classidx = self.servt_classes_updmap(r,4);
 
    % store the residence times and tput at this layer to become
    % the servt / tputs of aidx in another layer, as needed
    iter_min = min(30,ceil(self.options.iter_max/4));
    wnd_size = (it-self.averagingstart+1);
 
    % Compute residt from QN/TN_ref instead of WN to avoid
    % fork+loop visit distortion (WN uses visits from DTMC solve
    % which are distorted when Fork non-stochastic rows coexist
    % with loop back-edges in the routing matrix)
    layerIdx = self.idxhash(idx);
    layerSn = ensemble{layerIdx}.getStruct();
    c = find(layerSn.chains(:, classidx), 1);
    refclass_c = layerSn.refclass(c);
    refstat_k = layerSn.refstat(classidx);
 
    if ~isempty(self.averagingstart) && it>=iter_min % assume steady-state
        self.servt(aidx) = 0;
        self.residt(aidx) = 0;
        self.tput(aidx) = 0;
        for w=0:(wnd_size-1)
            self.servt(aidx) = self.servt(aidx) + self.results{end-w,layerIdx}.RN(nodeidx,classidx) / wnd_size;
            TN_ref = self.results{end-w,layerIdx}.TN(refstat_k, refclass_c);
            if TN_ref > GlobalConstants.FineTol
                self.residt(aidx) = self.residt(aidx) + self.results{end-w,layerIdx}.QN(nodeidx,classidx) / TN_ref / wnd_size;
            else
                self.residt(aidx) = self.residt(aidx) + self.results{end-w,layerIdx}.WN(nodeidx,classidx) / wnd_size;
            end
            self.tput(aidx) = self.tput(aidx) + self.results{end-w,layerIdx}.TN(nodeidx,classidx) / wnd_size;
        end
    else
        self.servt(aidx) = self.results{end,layerIdx}.RN(nodeidx,classidx);
        TN_ref = self.results{end,layerIdx}.TN(refstat_k, refclass_c);
        QN_val = self.results{end,layerIdx}.QN(nodeidx,classidx);
        if TN_ref > GlobalConstants.FineTol
            self.residt(aidx) = QN_val / TN_ref;
        else
            self.residt(aidx) = self.results{end,layerIdx}.WN(nodeidx,classidx);
        end
        self.tput(aidx) = self.results{end,layerIdx}.TN(nodeidx,classidx);
    end
 
    % Fix for async-only entry targets: use RN (response time per visit) for residt
    % The host layer closed model incorrectly splits residence time (WN) between
    % activities when an entry only receives async calls (no sync callers).
    % For async-only entries, use RN instead of WN since the async arrivals
    % don't share the closed chain's visit ratio - each async arrival gets
    % the full response time per visit.
    if aidx > lqn.ashift && aidx <= lqn.ashift + lqn.nacts
        % This is an activity - find its bound entry
        for eidx = (lqn.eshift+1):(lqn.eshift+lqn.nentries)
            if full(lqn.graph(eidx, aidx)) > 0
                % Found bound entry - check if async-only
                hasSyncCallers = full(any(lqn.issynccaller(:, eidx)));
                hasAsyncCallers = full(any(lqn.isasynccaller(:, eidx)));
                if hasAsyncCallers && ~hasSyncCallers
                    % Async-only target: use RN (response time per visit)
                    % instead of WN (residence time with visit ratio)
                    self.residt(aidx) = self.servt(aidx);  % servt already has RN
                end
                break;
            end
        end
    end
 
    % Recover from Inf/NaN: snap back to previous iteration's value
    if it > 1
        if (isinf(self.servt(aidx)) || isnan(self.servt(aidx))) && ~isnan(self.servt_prev(aidx))
            self.servt(aidx) = self.servt_prev(aidx);
        end
        if (isinf(self.residt(aidx)) || isnan(self.residt(aidx))) && ~isnan(self.residt_prev(aidx))
            self.residt(aidx) = self.residt_prev(aidx);
        end
        if (isinf(self.tput(aidx)) || isnan(self.tput(aidx))) && ~isnan(self.tput_prev(aidx))
            self.tput(aidx) = self.tput_prev(aidx);
        end
    end
 
    % Apply under-relaxation if enabled and not first iteration
    omega = self.relax_omega;
    if omega < 1.0 && it > 1
        if ~isnan(self.servt_prev(aidx))
            self.servt(aidx) = omega * self.servt(aidx) + (1 - omega) * self.servt_prev(aidx);
        end
        if ~isnan(self.residt_prev(aidx))
            self.residt(aidx) = omega * self.residt(aidx) + (1 - omega) * self.residt_prev(aidx);
        end
        if ~isnan(self.tput_prev(aidx))
            self.tput(aidx) = omega * self.tput(aidx) + (1 - omega) * self.tput_prev(aidx);
        end
    end
    % Store current values for next iteration
    self.servt_prev(aidx) = self.servt(aidx);
    self.residt_prev(aidx) = self.residt(aidx);
    self.tput_prev(aidx) = self.tput(aidx);
 
    % Safeguard against MVA numerical instability producing extreme values
    % (matches Python _update_metrics_default max_servt guard)
    max_servt = 1e10;
    if self.servt(aidx) > 0 && self.servt(aidx) <= max_servt
        self.servtproc{aidx} = Exp.fitMean(self.servt(aidx));
    end
    self.tputproc{aidx} = Exp.fitRate(self.tput(aidx));
end
 
% Phase-2 support: split activity service times by phase
% Note: Overtaking probability is computed later after entry throughput is available
if self.hasPhase2
    % Reset phase-specific arrays
    self.servt_ph1 = zeros(lqn.nidx, 1);
    self.servt_ph2 = zeros(lqn.nidx, 1);
 
    % Split activity service times by phase
    for a = 1:lqn.nacts
        aidx = lqn.ashift + a;
        if lqn.actphase(a) == 1
            self.servt_ph1(aidx) = self.servt(aidx);
        else
            self.servt_ph2(aidx) = self.servt(aidx);
        end
    end
 
    % Aggregate phase service times to entry level
    for e = 1:lqn.nentries
        eidx = lqn.eshift + e;
        acts = lqn.actsof{eidx};
        for aidx = acts
            a = aidx - lqn.ashift;
            if a > 0 && a <= lqn.nacts
                if lqn.actphase(a) == 1
                    self.servt_ph1(eidx) = self.servt_ph1(eidx) + self.servt_ph1(aidx);
                else
                    self.servt_ph2(eidx) = self.servt_ph2(eidx) + self.servt_ph2(aidx);
                end
            end
        end
    end
end
 
% obtain throughput for activities in thinkt_classes_updmap (needed for async calls)
% this ensures tputproc is set for activities that make async calls from client nodes
for r=1:size(self.thinkt_classes_updmap,1)
    idx = self.thinkt_classes_updmap(r,1);
    aidx = self.thinkt_classes_updmap(r,2);
    nodeidx = self.thinkt_classes_updmap(r,3);
    classidx = self.thinkt_classes_updmap(r,4);
 
    % only update if not already set by servt_classes_updmap processing
    if isempty(self.tputproc) || length(self.tputproc) < aidx || isempty(self.tputproc{aidx})
        iter_min = min(30,ceil(self.options.iter_max/4));
        wnd_size = (it-self.averagingstart+1);
        if ~isempty(self.averagingstart) && it>=iter_min % assume steady-state
            self.tput(aidx) = 0;
            for w=0:(wnd_size-1)
                self.tput(aidx) = self.tput(aidx) + self.results{end-w,self.idxhash(idx)}.TN(nodeidx,classidx) / wnd_size;
            end
        else
            self.tput(aidx) = self.results{end,self.idxhash(idx)}.TN(nodeidx,classidx);
        end
        self.tputproc{aidx} = Exp.fitRate(self.tput(aidx));
    end
end
 
% Obtain the join times for AND-Join activities
% For each AND-join activity, compute the synchronization delay as
% the maximum of predecessor branch service times.
self.joint = zeros(lqn.nidx,1);
joinedacts = find(lqn.actpretype == ActivityPrecedenceType.PRE_AND)';
for aidx = joinedacts
    % Find predecessor activities in the activity graph
    preds = find(lqn.graph(:, aidx) > 0)';
    pred_servts = [];
    for pidx = preds
        if pidx > lqn.ashift && pidx <= lqn.ashift + lqn.nacts
            pred_servts(end+1) = self.servt(pidx); %#ok<AGROW>
        end
    end
    if ~isempty(pred_servts)
        % Join time = maximum of predecessor branch service times
        % This represents the synchronization delay at the AND-join point.
        % For exponentially distributed branch times with means t_1,...,t_k,
        % approximate E[max] using the Clark (1961) two-moment method
        % recursively for k branches.
        if length(pred_servts) == 1
            self.joint(aidx) = pred_servts(1);
        else
            % Recursive two-moment approximation for max of exponentials
            % Start with first two branches, then fold in remaining
            mu_max = pred_servts(1);
            for bi = 2:length(pred_servts)
                t1 = mu_max;
                t2 = pred_servts(bi);
                if t1 > 0 && t2 > 0
                    % For exponentials: E[max(X1,X2)] = E[X1] + E[X2] - E[min(X1,X2)]
                    % E[min(X1,X2)] = 1/(1/t1 + 1/t2) = t1*t2/(t1+t2)
                    mu_max = t1 + t2 - t1*t2/(t1+t2);
                else
                    mu_max = max(t1, t2);
                end
            end
            self.joint(aidx) = mu_max;
        end
    end
end
 
% obtain the call residence time
self.callservt = zeros(lqn.ncalls,1);
self.callresidt = zeros(lqn.ncalls,1);
for r=1:size(self.call_classes_updmap,1)
    idx = self.call_classes_updmap(r,1);
    cidx = self.call_classes_updmap(r,2);
    nodeidx = self.call_classes_updmap(r,3);
    classidx = self.call_classes_updmap(r,4);
    if self.call_classes_updmap(r,3) > 1
        if nodeidx == 1
            self.callservt(cidx) = 0;
        else
            self.callservt(cidx) = self.results{end, self.idxhash(idx)}.RN(nodeidx,classidx) * self.lqn.callproc{cidx}.getMean;
            self.callresidt(cidx) = self.results{end, self.idxhash(idx)}.WN(nodeidx,classidx);
        end
        % Growth rate capping removed - it prevents callservt from converging
        % to the correct value when initial values are near-zero (Immediate)
        % Apply under-relaxation to call service times
        omega = self.relax_omega;
        if omega < 1.0 && it > 1 && ~isnan(self.callservt_prev(cidx))
            self.callservt(cidx) = omega * self.callservt(cidx) + (1 - omega) * self.callservt_prev(cidx);
        end
        self.callservt_prev(cidx) = self.callservt(cidx);
        self.callresidt_prev(cidx) = self.callresidt(cidx);
    end
end
 
% then resolve the entry servt summing up these contributions
entry_servt = self.servtmatrix*[self.residt;self.callresidt(:)];
entry_servt(1:lqn.eshift) = 0;
 
% Update forwarding call service times from the target entry's service time.
% Instead of bundling the FWD target's time into the source entry (which
% overestimates task utilization), we set the FWD call's own service time.
% The FWD call class in the task layer models the additional blocking time
% the caller sees while the forwarding target processes the request.
for cidx = 1:lqn.ncalls
    if lqn.calltype(cidx) == CallType.FWD
        target_eidx = lqn.callpair(cidx, 2);
        fwd_prob = lqn.callproc{cidx}.getMean();
 
        % Get target entry's service time
        target_servt = entry_servt(target_eidx);
 
        % If target entry doesn't have computed service time,
        % compute it directly from its activities' host demands
        if target_servt == 0 || isnan(target_servt)
            target_servt = 0;
            for aidx = lqn.actsof{target_eidx}
                if ~isempty(lqn.hostdem{aidx})
                    target_servt = target_servt + lqn.hostdem{aidx}.getMean();
                end
            end
        end
 
        if target_servt > 0
            self.callservt(cidx) = target_servt * fwd_prob;
            self.callresidt(cidx) = self.callservt(cidx);
            self.callservtproc{cidx} = Exp.fitMean(self.callservt(cidx));
        end
    end
end
 
% For phase-2 entries, the task-layer sync call class sees the rest of the
% entry execution after the call returns, which overcounts the call demand.
% Re-anchor the sync call to the target entry delay before layering in any
% forwarding-chain blocking.
for cidx = 1:lqn.ncalls
    if lqn.calltype(cidx) == CallType.SYNC
        src_aidx = lqn.callpair(cidx, 1);
        src_eidx = [];
        src_tidx = lqn.parent(src_aidx);
        for candidate_eidx = lqn.entriesof{src_tidx}
            if any(lqn.actsof{candidate_eidx} == src_aidx)
                src_eidx = candidate_eidx;
                break;
            end
        end
        if isempty(src_eidx)
            continue;
        end
        if self.hasPhase2 && self.servt_ph2(src_eidx) > GlobalConstants.FineTol
            target_eidx = lqn.callpair(cidx, 2);
            callmean = lqn.callproc{cidx}.getMean();
            target_servt = entry_servt(target_eidx);
            if target_servt == 0 || isnan(target_servt)
                target_servt = 0;
                for aidx = lqn.actsof{target_eidx}
                    if ~isempty(lqn.hostdem{aidx})
                        target_servt = target_servt + lqn.hostdem{aidx}.getMean();
                    end
                end
            end
            self.callservt(cidx) = target_servt * callmean;
            self.callresidt(cidx) = self.callservt(cidx);
        end
    end
end
 
% Add forwarding chain delays to the SYNC calls.
% The caller is blocked until the entire forwarding chain completes,
% so the SYNC call's service time must include all FWD delays.
% Uses BFS to follow chains (e0->e1->e2) and accumulate delays.
for sync_cidx = 1:lqn.ncalls
    if lqn.calltype(sync_cidx) == CallType.SYNC
        total_fwd_delay = 0;
        total_fwd_residt = 0;
        entries_to_scan = lqn.callpair(sync_cidx, 2);
        scanned_entries = [];
        while ~isempty(entries_to_scan)
            eidx = entries_to_scan(1);
            entries_to_scan(1) = [];
            if ismember(eidx, scanned_entries), continue; end
            scanned_entries(end+1) = eidx; %#ok<AGROW>
            for fwd_cidx = 1:lqn.ncalls
                if lqn.calltype(fwd_cidx) == CallType.FWD && lqn.callpair(fwd_cidx, 1) == eidx
                    total_fwd_delay = total_fwd_delay + self.callservt(fwd_cidx);
                    total_fwd_residt = total_fwd_residt + self.callresidt(fwd_cidx);
                    fwd_tgt = lqn.callpair(fwd_cidx, 2);
                    if ~ismember(fwd_tgt, scanned_entries) && ~ismember(fwd_tgt, entries_to_scan)
                        entries_to_scan(end+1) = fwd_tgt; %#ok<AGROW>
                    end
                end
            end
        end
        if total_fwd_delay > 0
            self.callservt(sync_cidx) = self.callservt(sync_cidx) + total_fwd_delay;
            self.callresidt(sync_cidx) = self.callresidt(sync_cidx) + total_fwd_residt;
        end
    end
end
% Recompute entry_servt with forwarding-adjusted callresidt
entry_servt = self.servtmatrix*[self.residt;self.callresidt(:)];
entry_servt(1:lqn.eshift) = 0;
 
% this block fixes the problem that ResidT is scaled so that the
% task has Vtask=1, but in call servt the entries need to have Ventry=1
for eidx=(lqn.eshift+1):(lqn.eshift+lqn.nentries)
    tidx = lqn.parent(eidx); % task of entry
    hidx = lqn.parent(tidx); %host of entry
    if ~self.ignore(tidx) && ~self.ignore(hidx)
        % Check if this entry has sync callers (which create closed classes)
        hasSyncCallers = full(any(lqn.issynccaller(:, eidx)));
 
        if hasSyncCallers
            % Original logic for entries with sync callers
            % get class in host layer of task and entry
            tidxclass = ensemble{self.idxhash(hidx)}.attribute.tasks(find(ensemble{self.idxhash(hidx)}.attribute.tasks(:,2) == tidx),1);
            eidxclass = ensemble{self.idxhash(hidx)}.attribute.entries(find(ensemble{self.idxhash(hidx)}.attribute.entries(:,2) == eidx),1);
            task_tput  = sum(self.results{end,self.idxhash(hidx)}.TN(ensemble{self.idxhash(hidx)}.attribute.clientIdx,tidxclass));
            entry_tput = sum(self.results{end,self.idxhash(hidx)}.TN(ensemble{self.idxhash(hidx)}.attribute.clientIdx,eidxclass));
            if entry_tput > GlobalConstants.Zero
                self.servt(eidx) = entry_servt(eidx) * task_tput / entry_tput;
                self.residt(eidx) = entry_servt(eidx) * task_tput / entry_tput;
            else
                self.servt(eidx) = entry_servt(eidx);
                self.residt(eidx) = entry_servt(eidx);
            end
        else
            % For async-only targets, use entry_servt directly
            % No throughput ratio scaling needed since there are no closed classes
            self.servt(eidx) = entry_servt(eidx);
            self.residt(eidx) = entry_servt(eidx);
        end
    end
end
 
% Phase-2 support: compute overtaking probability and apply correction
% This must happen AFTER entry throughput is available (computed above)
if self.hasPhase2
    for e = 1:lqn.nentries
        eidx = lqn.eshift + e;
        tidx = lqn.parent(eidx);
 
        if self.servt_ph2(eidx) > GlobalConstants.FineTol
            if lqn.isref(tidx) || ~full(any(lqn.issynccaller(:, eidx)))
                self.residt(eidx) = self.servt(eidx);
                continue;
            end
            % Get entry throughput (use task throughput as approximation if entry not available)
            if self.tput(eidx) > GlobalConstants.FineTol
                entry_tput = self.tput(eidx);
            elseif self.tput(tidx) > GlobalConstants.FineTol
                entry_tput = self.tput(tidx);
            else
                entry_tput = 0;
            end
 
            % Compute overtaking probability now that throughput is available
            if entry_tput > GlobalConstants.FineTol
                self.prOvertake(e) = self.overtake_prob(eidx);
            else
                self.prOvertake(e) = 0;
            end
 
            % Caller's response time = phase-1 only + P(overtake) * phase-2
            % Phase-2 only delays if overtaking occurs
            overtake_delay = self.prOvertake(e) * self.servt_ph2(eidx);
 
            % The caller sees phase-1 + overtaking correction (not full phase-2)
            % self.servt(eidx) remains unchanged (phase-1 + phase-2) for utilization calculation
            self.residt(eidx) = self.servt_ph1(eidx) + overtake_delay;
        end
    end
end
 
%self.servt(lqn.eshift+1:lqn.eshift+lqn.nentries) = entry_servt(lqn.eshift+1:lqn.eshift+lqn.nentries);
%entry_servt((lqn.ashift+1):end) = 0;
for r=1:size(self.call_classes_updmap,1)
    cidx = self.call_classes_updmap(r,2);
    eidx = lqn.callpair(cidx,2);
    if self.call_classes_updmap(r,3) > 1
        if self.servt(eidx) > 0
            self.servtproc{eidx} = Exp.fitMean(self.servt(eidx));
        end
    end
end
 
% determine call response times processes
for r=1:size(self.call_classes_updmap,1)
    cidx = self.call_classes_updmap(r,2);
    eidx = lqn.callpair(cidx,2);
    if self.call_classes_updmap(r,3) > 1
        if it==1
            % note that respt is per visit, so number of calls is 1
            self.callservt(cidx) = self.servt(eidx);
            self.callservtproc{cidx} = self.servtproc{eidx};
        else
            % note that respt is per visit, so number of calls is 1
            if self.callservt(cidx) > 0
                self.callservtproc{cidx} = Exp.fitMean(self.callservt(cidx));
            end
        end
    end
end
 
self.ensemble = ensemble;
end