solver_ctmc_transient_analyzer.m

function [t,pit,QNt,UNt,RNt,TNt,CNt,XNt,InfGen,StateSpace,StateSpaceAggr,EventFiltration,runtime,fname] = solver_ctmc_transient_analyzer(sn, options)
% [T,PIT,QNT,UNT,RNT,TNT,CNT,XNT,INFGEN,STATESPACE,STATESPACEAGGR,EVENTFILTRATION,RUNTIME,FNAME] = SOLVER_CTMC_TRANSIENT_ANALYZER(QN, OPTIONS)
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
RNt=[]; CNt=[];  XNt=[];
 
M = sn.nstations;    %number of stations
K = sn.nclasses;    %number of classes
fname = '';
Tstart = tic;
S = sn.nservers;
sched = sn.sched;
PH = sn.proc;
 
line_debug('CTMC transient analyzer starting: nstations=%d, nclasses=%d', M, K);
 
[InfGen,StateSpace,StateSpaceAggr,EventFiltration,~,depRates,sn] = solver_ctmc(sn, options); % sn is updated with the state space
 
if options.keep
    fname = lineTempName;
    save([fname,'.mat'],'InfGen','StateSpace','StateSpaceAggr','EventFiltration')
    line_printf('\nCTMC infinitesimal generator and state space saved in: ');
    line_printf([fname, '.mat'])
end
 
state = [];
for ist=1:sn.nnodes
    if sn.isstateful(ist)
        isf = sn.nodeToStateful(ist);
        state = [state,zeros(1,size(sn.space{isf},2)-length(sn.state{isf})),sn.state{isf}];
    end
end
pi0 = zeros(1,length(InfGen));
 
state0 = matchrow(StateSpace, state);
if state0 == -1
    state0 = matchrow(StateSpace, round(state));
    state = round(state);
    if state0 == -1
        line_error(mfilename,'Initial state not contained in the state space.');
    end
end
pi0(state0) = 1; % find initial state and set it to probability 1
 
% Time-varying generator support. When options.config.rate_sched is supplied
% (per-(station,class) rate trajectories, e.g. injected by the coupled LN
% transient), the CTMC becomes time-inhomogeneous: the transitions attributable
% to a scaled (station,class) service are modulated by m(t)=rate(t)/nominal.
% The generator is linear in sn.rates(ist,r), so each scaled component
% Qhat_sc = (Q(scaled) - Q_base)/(probe-1) is extracted with one extra
% solver_ctmc build, and Q(t) = Q_base + sum_sc (m_sc(t)-1) Qhat_sc is
% integrated by non-homogeneous uniformization over piecewise-constant segments.
rate_sched = [];
if isfield(options,'config') && isfield(options.config,'rate_sched') && ~isempty(options.config.rate_sched)
    rate_sched = options.config.rate_sched;
end
if isempty(rate_sched)
    [pit,t] = ctmc_transient(InfGen,pi0,options.timespan(1),options.timespan(2),options.stiff,[],options.timestep);
    mscale = [];
else
    [pit,t,mscale] = local_ctmc_timevarying(sn, options, InfGen, StateSpace, pi0, rate_sched, M, K);
end
pit(pit<GlobalConstants.Zero)=0;
 
QNt = cell(M,K);
UNt = cell(M,K);
%XNt = cell(1,K);
TNt = cell(M,K);
 
if t(1) == 0
    t(1) = GlobalConstants.Zero;
end
for k=1:K
    %    XNt(k) = pi*arvRates(:,sn.refstat(k),k);
    for ist=1:M
        %occupancy_t = cumsum(pit.*[0;diff(t)],1)./t;
        occupancy_t = pit;
        TNt{ist,k} = occupancy_t*depRates(:,ist,k);
        % Under a time-varying rate schedule, the throughput of a scaled
        % (station,class) is modulated by its multiplier m(t); depRates is the
        % nominal-rate departure rate, so scale it pointwise in time.
        if ~isempty(mscale)
            TNt{ist,k} = TNt{ist,k} .* squeeze(mscale(ist,k,:));
        end
        % A Source station holds an unbounded (Inf) population, so its
        % StateSpaceAggr column is Inf and QNt/UNt are undefined
        % (pit*Inf yields Inf/NaN). Mirror the steady-state analyzer, which
        % leaves the Source QN/UN at 0 and reports only its throughput.
        if sn.nodetype(sn.stationToNode(ist)) == NodeType.Source
            QNt{ist,k} = zeros(size(pit,1),1);
            UNt{ist,k} = zeros(size(pit,1),1);
            continue
        end
        qlenAt_t = pit*StateSpaceAggr(:,(ist-1)*K+k);
        %QNt{i,k} = cumsum(qlenAt_t.*[0;diff(t)])./t;
        QNt{ist,k} = qlenAt_t;
        switch sched(ist)
            case SchedStrategy.INF
                UNt{ist,k} = QNt{ist,k};
            case {SchedStrategy.FCFS, SchedStrategy.HOL, SchedStrategy.SIRO, SchedStrategy.SEPT, SchedStrategy.LEPT, SchedStrategy.SJF}
                if ~isempty(PH{ist}{k})
                    UNt{ist,k} = occupancy_t*min(StateSpaceAggr(:,(ist-1)*K+k),S(ist))/S(ist);
                end
            case SchedStrategy.PS
                uik = min(StateSpaceAggr(:,(ist-1)*K+k),S(ist)) .* StateSpaceAggr(:,(ist-1)*K+k) ./ sum(StateSpaceAggr(:,((ist-1)*K+1):(ist*K)),2);
                uik(isnan(uik))=0;
                utilAt_t = pit * uik / S(ist);
                %UNt{i,k} = cumsum(utilAt_t.*[0;diff(t)])./t;
                UNt{ist,k} = utilAt_t;
            case SchedStrategy.DPS
                w = sn.schedparam(ist,:);
                nik = S(ist) * w(k) * StateSpaceAggr(:,(ist-1)*K+k) ./ sum(repmat(w,size(StateSpaceAggr,1),1).*StateSpaceAggr(:,((ist-1)*K+1):(ist*K)),2);
                nik(isnan(nik))=0;
                UNt{ist,k} = occupancy_t*nik;
            otherwise
                if ~isempty(PH{ist}{k})
                    ind = sn.stationToNode(ist);
                    line_warning(mfilename,'Transient utilization not support yet for station %s, returning an approximation.\n',sn.nodenames{ind});
                    UNt{ist,k} = occupancy_t*min(StateSpaceAggr(:,(ist-1)*K+k),S(ist))/S(ist);
                end
        end
    end
end
runtime = toc(Tstart);
 
%if options.verbose
%    line_printf('\nCTMC analysis completed. Runtime: %f seconds.\n',runtime);
%end
end
 
function [pit, t, mscale] = local_ctmc_timevarying(sn, options, Qbase, StateSpace, pi0, rate_sched, M, K)
% Integrate the time-inhomogeneous forward equation dpi/dt = pi Q(t) by
% non-homogeneous uniformization over a piecewise-constant grid. The generator
% is Q(t) = Qbase + sum_sc (m_sc(t)-1) Qhat_sc, where Qhat_sc is the linear
% component of Qbase attributable to the scaled (station,class) rate,
% extracted by a single probe rebuild (Q is linear in sn.rates).
ts = options.timespan;
Ngrid = 100;
if isfield(options,'config') && isfield(options.config,'ctmc_tv_ngrid') && ~isempty(options.config.ctmc_tv_ngrid)
    Ngrid = options.config.ctmc_tv_ngrid;
end
t = linspace(ts(1), ts(2), Ngrid)';
nt = numel(t);
nS = size(Qbase,1);
 
% Build the per-schedule generator component and multiplier trajectory.
probe = 2.0;
nsc = numel(rate_sched);
Qhat = cell(1,nsc);
mtraj = ones(nt, nsc);
scStation = zeros(1,nsc); scClass = zeros(1,nsc);
for s = 1:nsc
    ist = rate_sched(s).station;
    r = rate_sched(s).class;
    scStation(s) = ist; scClass(s) = r;
    % probe rebuild with the (ist,r) service process time-scaled by `probe`;
    % the reachable state space is rate-independent, so ordering matches Qbase.
    % solver_ctmc builds transitions from the process representation (sn.proc
    % MAP D0/D1 and sn.mu), NOT from sn.rates, so all rate-carrying fields must
    % be scaled. Time-scaling a MAP/PH by f multiplies D0 and D1 (and the phase
    % rates mu) by f, scaling every completion/phase rate by f.
    snp = sn;
    snp.rates(ist,r) = snp.rates(ist,r) * probe;
    if iscell(snp.proc) && numel(snp.proc) >= ist && iscell(snp.proc{ist}) && numel(snp.proc{ist}) >= r
        pr = snp.proc{ist}{r};
        if iscell(pr)
            for z = 1:numel(pr)
                if isnumeric(pr{z})
                    pr{z} = pr{z} * probe;
                end
            end
            snp.proc{ist}{r} = pr;
        end
    end
    if iscell(snp.mu) && numel(snp.mu) >= ist && iscell(snp.mu{ist}) && numel(snp.mu{ist}) >= r
        snp.mu{ist}{r} = snp.mu{ist}{r} * probe;
    end
    Qp = solver_ctmc(snp, options);
    if size(Qp,1) ~= nS
        line_error(mfilename, 'rate_sched probe changed the CTMC state-space size; cannot build time-varying generator.');
    end
    Qhat{s} = (Qp - Qbase) / (probe - 1);
    % multiplier m(t) = rate(t)/nominal (nominal defaults to sn.rates(ist,r))
    if isfield(rate_sched(s),'nominal') && ~isempty(rate_sched(s).nominal)
        nominal = rate_sched(s).nominal;
    else
        nominal = sn.rates(ist,r);
    end
    seg_t = rate_sched(s).tgrid(:)';
    seg_r = rate_sched(s).rates(:)';
    mtraj(:,s) = interp1(seg_t, seg_r, min(max(t,seg_t(1)),seg_t(end)), 'linear') / nominal;
end
 
% Propagate one segment per grid interval, with the multiplier frozen at the
% interval midpoint (second-order in dt). Over a segment the generator Qk is
% constant, so the exact forward solution is pi(t+dt) = pi(t) expm(Qk dt).
% We deliberately use the matrix exponential rather than ctmc_uniformization:
% an LN layer can carry a near-instantaneous reply-signal sentinel rate (~1e9),
% making the generator stiff (q dt ~ 1e8); uniformization then splits the step
% into >1e5 sub-segments and leaks all probability mass to zero. expm of a
% valid generator is exactly stochastic, so mass is conserved for any
% stiffness. For a large sentinel rate the corresponding transient state simply
% empties within the segment, which is the intended near-instantaneous limit.
pit = zeros(nt, nS);
pit(1,:) = pi0(:)';
Qfull = full(Qbase);
QhatFull = cell(1,nsc);
for s = 1:nsc
    QhatFull{s} = full(Qhat{s});
end
for k = 1:nt-1
    dt = t(k+1) - t(k);
    Qk = Qfull;
    for s = 1:nsc
        mk = 0.5*(mtraj(k,s) + mtraj(k+1,s));
        Qk = Qk + (mk - 1) * QhatFull{s};
    end
    pit(k+1,:) = pit(k,:) * expm(Qk * dt);
end
 
% Per-(station,class) throughput multiplier over time (1 where not scaled).
mscale = ones(M, K, nt);
for s = 1:nsc
    mscale(scStation(s), scClass(s), :) = reshape(mtraj(:,s), 1, 1, nt);
end
end