runTransientJson.m

function res = runTransientJson(self, numEvents)
% RES = RUNTRANSIENTJSON(NUMEVENTS)
% Run a fully JSON-mediated transient LDES simulation over the horizon
% [0, samples] (or [0, numEvents] when numEvents>0) with trajectory export,
% and return the parsed trajectory. Mirrors the Python-native _run_transient.
%
% RES has fields:
%   .t   - (nTimePoints x 1) time vector, or [] if none
%   .QNt - 1 x Nouter cell; RES.QNt{i} is a 1 x R cell whose {r} entry is the
%          (nTimePoints x 2) [value, time] trajectory of class r at the i-th
%          stateful node, or [] if absent. Nouter matches the engine's outer
%          transient dimension (stateful-node major, as consumed by sample()).
 
options = self.getOptions;
% The horizon is passed explicitly via --timespan; options.samples (the
% steady-state event budget) is left untouched since the engine ignores it
% in transient mode.
S = options.samples;
if nargin >= 2 && ~isempty(numEvents) && numEvents > 0
    S = numEvents;
end
extraFlags = {'--timespan', sprintf('0,%.10g', S), '--trajectory'};
 
% Ensemble transient: a single realization is not the transient mean E[N](t)
% (no time-ergodicity at fixed t). When replications > 1 is requested, drive
% the engine's parallel analyzer via --replications so the per-bucket QNt/UNt/
% TNt are averaged across independent replications. Read from options.replications
% or options.config.replications; default (absent/<=1) keeps the single-path run.
reps = 0;
if isfield(options, 'replications') && ~isempty(options.replications)
    reps = options.replications;
elseif isfield(options, 'config') && isfield(options.config, 'replications') ...
        && ~isempty(options.config.replications)
    reps = options.config.replications;
end
if reps > 1
    extraFlags{end+1} = '--replications';
    extraFlags{end+1} = sprintf('%d', round(reps));
    if isfield(options, 'config') && isfield(options.config, 'numthreads') ...
            && ~isempty(options.config.numthreads) && options.config.numthreads > 0
        extraFlags{end+1} = '--numthreads';
        extraFlags{end+1} = sprintf('%d', round(options.config.numthreads));
    end
end
 
data = self.solveCli(options, extraFlags);
 
res = struct('t', [], 'QNt', {{}});
if ~isstruct(data) || ~isfield(data, 'transient') || isempty(data.transient)
    return;
end
tran = data.transient;
if isfield(tran, 't') && ~isempty(tran.t)
    tvec = ldesJson2mat(tran.t, [], []);
    res.t = tvec(:);
end
if isfield(tran, 'QNt') && ~isempty(tran.QNt)
    sn = self.model.getStruct;
    Nouter = sn.nstateful;
    R = sn.nclasses;
    % Normalize into an Nouter x R cell of [nRows x 2] matrices, then repackage
    % as a 1 x Nouter cell of 1 x R class-cells (the layout sample() consumes).
    grid = ldesTrajCell(tran.QNt, Nouter, R);
    res.QNt = cell(1, Nouter);
    for i = 1:Nouter
        classCells = cell(1, R);
        for r = 1:R
            classCells{r} = grid{i, r};
        end
        res.QNt{i} = classCells;
    end
end
end