linemodel_save.m

function linemodel_save(model, filename)
% LINEMODEL_SAVE Save a LINE model to JSON.
%
%   LINEMODEL_SAVE(MODEL, FILENAME) saves the model to the specified JSON
%   file, conforming to the line-model.schema.json specification.
%
% Parameters:
%   model    - Network, LayeredNetwork, Workflow, or Environment object
%   filename - output file path (should end in .json)
%
% Example:
%   model = Network('M/M/1');
%   source = Source(model, 'Source');
%   queue  = Queue(model, 'Queue', SchedStrategy.FCFS);
%   sink   = Sink(model, 'Sink');
%   oclass = OpenClass(model, 'Class1');
%   source.setArrival(oclass, Exp(1.0));
%   queue.setService(oclass, Exp(2.0));
%   P = model.initRoutingMatrix();
%   P{1}(1,2) = 1; P{1}(2,3) = 1;
%   model.link(P);
%   linemodel_save(model, 'mm1.json');
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
if isa(model, 'LayeredNetwork')
    modelMap = layered2json(model);
elseif isa(model, 'Workflow')
    modelMap = workflow2json(model);
elseif isa(model, 'Environment')
    modelMap = environment2json(model);
else
    modelMap = network2json(model);
end
 
% Build the full document
sb = {};
sb{end+1} = '{';
sb{end+1} = '  "format": "line-model",';
sb{end+1} = '  "version": "1.0",';
sb{end+1} = ['  "model": ', encode_value(modelMap, 2)];
sb{end+1} = '}';
jsonStr = strjoin(sb, newline);
 
fid = fopen(filename, 'w');
if fid == -1
    error('linemodel_save:fileOpen', 'Cannot open file: %s', filename);
end
cleanupObj = onCleanup(@() fclose(fid));
fprintf(fid, '%s\n', jsonStr);
end
 
 
% =========================================================================
%  Network serialization
% =========================================================================
 
function result = network2json(model)
% Convert a Network to a containers.Map (preserves key ordering/commas)
result = containers.Map();
result('type') = 'Network';
result('name') = model.getName();
 
nodes = model.getNodes();
classes = model.getClasses();
K = length(classes);
M = length(nodes);
 
% --- Nodes ---
nodesJson = {};
for i = 1:M
    node = nodes{i};
 
    % Skip implicit ClassSwitch nodes (auto-created by link())
    if isa(node, 'ClassSwitch') && isprop(node, 'autoAdded') && node.autoAdded
        continue;
    end
 
    nj = containers.Map();
    nj('name') = node.name;
    nj('type') = node_type_str(node);
 
    % Scheduling. A queueing Place carries an embedded-queue scheduling strategy
    % just as a Queue does, but Place extends Station rather than Queue, so an
    % isa(node,'Queue') guard drops it (and everything below) on the floor. The
    % JAR (LineModelIO, "Place" writer branch) is the reference here.
    if isa(node, 'Delay')
        nj('scheduling') = 'INF';
    elseif isa(node, 'Place')
        if node.isQueueing() && ~isempty(node.schedStrategy)
            nj('scheduling') = sched_id_to_str(node.schedStrategy);
        end
    elseif isa(node, 'Queue')
        sched = node.schedStrategy;
        if ~isempty(sched)
            nj('scheduling') = sched_id_to_str(sched);
        end
    end
 
    % Servers. An infinite server count is deliberately NOT emitted: in the LINE
    % object model it is not independent state that could be lost. A station has
    % numberOfServers == Inf if and only if it is INF-scheduled -- Queue's ctor
    % sets Inf for SchedStrategy.INF, Queue.setNumberOfServers ignores the call on
    % an INF queue, and Place.installQueueServer resets a non-INF place to 1 -- so
    % 'scheduling' already carries it losslessly and every reader (this one, and
    % the JAR via parseSchedStrategy + the Queue/Place ctors) rebuilds the
    % infinite-server station from it. Emitting "Infinity" here would be a
    % redundant key that current JAR readers reject outright
    % (LineModelIO getAsInt -> NumberFormatException), breaking every INF-Queue
    % and INF-Place model sent to the LDES engine.
    if isa(node, 'Station') && ~isa(node, 'Delay') && ...
            (isa(node, 'Queue') || (isa(node, 'Place') && node.isQueueing()))
        ns = node.numberOfServers;
        if isfinite(ns) && ns > 1
            nj('servers') = ns;
        end
    end
 
    % Buffer. Station, not Queue: Place extends Station directly, so a queueing
    % Place's capacity was dropped by the old isa(node,'Queue') guard.
    if isa(node, 'Station')
        c = node.cap;
        if ~isempty(c) && isfinite(c) && c > 0
            nj('buffer') = c;
        end
    end
 
    % Per-class buffer capacity
    if isa(node, 'Station') && ~isempty(node.classCap)
        ccMap = containers.Map();
        for r = 1:K
            jc = classes{r};
            if r <= length(node.classCap) && isfinite(node.classCap(r))
                ccMap(jc.name) = node.classCap(r);
            end
        end
        if ccMap.Count > 0
            nj('classCap') = ccMap;
        end
    end
 
    % Drop rules
    if isa(node, 'Station') && ~isempty(node.dropRule)
        drMap = containers.Map();
        for r = 1:K
            jc = classes{r};
            if r <= length(node.dropRule)
                dr = node.dropRule(r);
                drStr = droprule_to_str(dr);
                if ~isempty(drStr)
                    drMap(jc.name) = drStr;
                end
            end
        end
        if drMap.Count > 0
            nj('dropRule') = drMap;
        end
    end
 
    % Load-dependent scaling
    if isa(node, 'Station') && ~isempty(node.lldScaling)
        ldMap = containers.Map();
        ldMap('type') = 'loadDependent';
        ldMap('scaling') = node.lldScaling(:)';
        nj('loadDependence') = ldMap;
    end
 
    % Class-dependent scaling beta_{i,r}(n). The handle cannot cross the JSON
    % boundary, so it is materialized over the per-class box lattice exactly as
    % the OI/PAS rate table is (see oi_rate_table); the reader rebuilds a handle
    % that clamps to the cutoffs. Without this the LDES engine, which is a JSON
    % subprocess, would receive no class dependence and silently simulate the
    % unscaled network.
    if isa(node, 'Station') && ~isempty(node.lcdScaling)
        maxc = zeros(1, K);
        for r = 1:K
            if isa(classes{r}, 'ClosedClass') && isfinite(classes{r}.population)
                maxc(r) = round(classes{r}.population);
            else
                maxc(r) = 10;   % open-class saturation cutoff (beta clamped beyond)
            end
        end
        cdMap = containers.Map();
        cdMap('type') = 'classDependent';
        % num2cell keeps a single-class (1x1) vector from collapsing to a scalar.
        cdMap('cutoffs') = num2cell(double(maxc(:)'));
        cdMap('scaling') = cd_scaling_table(node.lcdScaling, maxc, K);
        % Declared peak rate scaling per class (Util = T*S/peak). Broadcast a
        % scalar to K entries so the reader restores a per-class vector.
        pk = node.lcdScalingPeak;
        if isscalar(pk)
            pk = repmat(pk, 1, K);
        end
        cdMap('peak') = num2cell(double(pk(:)'));
        nj('classDependence') = cdMap;
    end
 
    % Service / arrival distributions. A queueing Place carries per-class service
    % processes in its embedded queue, together with the departure discipline of
    % its depository; both are emitted here (JAR LineModelIO "Place" branch is the
    % reference). An ordinary Place has neither.
    svc = containers.Map();
    depDisc = containers.Map();
    for r = 1:K
        jc = classes{r};
        dist = [];
        if isa(node, 'Source')
            try
                dist = node.getArrivalProcess(jc);
            catch
                dist = [];
            end
        elseif isa(node, 'Place')
            if node.isQueueing() && numel(node.serviceProcess) >= jc.index ...
                    && ~isempty(node.serviceProcess{jc.index})
                dist = node.serviceProcess{jc.index};
            end
        elseif isa(node, 'Queue') || isa(node, 'Delay')
            if isa(node, 'Queue') && ~isempty(node.svcRateFun)
                % An OI/PAS queue is parameterized by the total rate function
                % mu(c) alone, carried below as oiServiceRate. getService would
                % return the internally materialized per-class representative,
                % which is derived state, is not settable (Queue.setService
                % rejects a distribution on a PAS/OI queue), and would make the
                % file unloadable.
                dist = [];
            else
                try
                    dist = node.getService(jc);
                catch
                    dist = [];
                end
            end
        end
        if ~isempty(dist) && ~isa(dist, 'Disabled')
            dj = dist2json(dist);
            if ~isempty(dj)
                svc(jc.name) = dj;
                if isa(node, 'Place') && numel(node.departureDiscipline) >= jc.index
                    depDisc(jc.name) = depdisc_to_str(node.departureDiscipline(jc.index));
                end
            end
        end
    end
    if svc.Count > 0
        nj('service') = svc;
    end
    if depDisc.Count > 0
        nj('departureDiscipline') = depDisc;
    end
 
    % Order-independent / pass-and-swap (OI/PAS) service: serialize the total
    % rate function mu(c) as a macrostate table keyed by the per-class counts
    % (mu is permutation-invariant, hence a function of the class counts only).
    % The JAR reconstructs mu(c) from this table; open classes are clamped to
    % the per-class cutoff. PAS additionally carries its swap graph.
    if isa(node, 'Queue') && ~isempty(node.svcRateFun)
        maxc = zeros(1, K);
        for r = 1:K
            if isa(classes{r}, 'ClosedClass') && isfinite(classes{r}.population)
                maxc(r) = round(classes{r}.population);
            else
                maxc(r) = 10;   % open-class saturation cutoff (mu constant beyond)
            end
        end
        rateMap = oi_rate_table(node.svcRateFun, maxc, K);
        nj('oiServiceRate') = rateMap;
        % num2cell keeps jsonencode from collapsing a single-class (1x1) vector
        % to a scalar, which the JAR reader parses as an array.
        nj('oiCutoffs') = num2cell(double(maxc(:)'));
        if node.schedStrategy == SchedStrategy.PAS && ~isempty(node.swapGraph)
            sgRows = cell(1, size(node.swapGraph, 1));
            for sgi = 1:size(node.swapGraph, 1)
                sgRows{sgi} = num2cell(double(node.swapGraph(sgi, :)));
            end
            nj('swapGraph') = sgRows;
        end
    end
 
    % Batch arrivals: the batch-size law released at each arrival epoch, per
    % class. Separate from 'service' above, which only spaces the epochs.
    if isa(node, 'Source') && ~isempty(node.arrivalBatch)
        batchMap = containers.Map();
        for r = 1:K
            if numel(node.arrivalBatch) >= r && ~isempty(node.arrivalBatch{r})
                bj = dist2json(node.arrivalBatch{r});
                if ~isempty(bj)
                    batchMap(classes{r}.name) = bj;
                end
            end
        end
        if batchMap.Count > 0
            nj('arrivalBatch') = batchMap;
        end
    end
 
    % Marked (MMAP) arrival binding: class names ordered by mark
    if isa(node, 'Source') && ~isempty(node.markedClasses)
        markedNames = cell(1, numel(node.markedClasses));
        for km = 1:numel(node.markedClasses)
            markedNames{km} = classes{node.markedClasses(km)}.name;
        end
        nj('markedClasses') = markedNames;
    end
 
    % ClassSwitch matrix
    if isa(node, 'ClassSwitch')
        csm = node.server.csMatrix;
        if ~isempty(csm)
            csDict = containers.Map();
            for ri = 1:K
                row = containers.Map();
                for ci = 1:K
                    if ri <= size(csm,1) && ci <= size(csm,2) && csm(ri,ci) ~= 0
                        row(classes{ci}.name) = csm(ri,ci);
                    end
                end
                if row.Count > 0
                    csDict(classes{ri}.name) = row;
                end
            end
            if csDict.Count > 0
                nj('classSwitchMatrix') = csDict;
            end
        end
    end
 
    % Cache config
    if isa(node, 'Cache')
        cc = containers.Map();
        cc('items') = node.items.nitems;
        ilc = node.itemLevelCap;
        % The replacement policy is emitted verbatim: readers decode every policy
        % natively, and CLIMB is rewritten into its FIFO unit-capacity-list form at
        % solve time by refreshLocalVars, which remaps itemcap and accost together.
        % Remapping here instead would emit the rewritten capacity next to the
        % original accessProb, leaving the two geometries inconsistent.
        if isscalar(ilc)
            cc('capacity') = ilc;
        else
            cc('capacity') = ilc(:)';
        end
        cc('replacement') = repl_to_str(node.replacestrategy);
        if isprop(node,'admissionProb') && ~isempty(node.admissionProb)
            cc('admissionProb') = node.admissionProb;
        end
 
        % Hit/miss class mappings
        hc = full(node.server.hitClass);
        mc = full(node.server.missClass);
        if ~isempty(hc) && any(hc > 0)
            hitMap = containers.Map();
            for hi = 1:length(hc)
                if hc(hi) > 0 && hi <= K && hc(hi) <= K
                    hitMap(classes{hi}.name) = classes{hc(hi)}.name;
                end
            end
            if hitMap.Count > 0
                cc('hitClass') = hitMap;
            end
        end
        if ~isempty(mc) && any(mc > 0)
            missMap = containers.Map();
            for mi = 1:length(mc)
                if mc(mi) > 0 && mi <= K && mc(mi) <= K
                    missMap(classes{mi}.name) = classes{mc(mi)}.name;
                end
            end
            if missMap.Count > 0
                cc('missClass') = missMap;
            end
        end
 
        % Read popularity distributions (setRead)
        if ~isempty(node.popularity)
            popMap = containers.Map();
            for pi = 1:size(node.popularity, 1)
                for pj = 1:size(node.popularity, 2)
                    if pi <= size(node.popularity, 1) && pj <= size(node.popularity, 2) ...
                            && ~isempty(node.popularity{pi, pj})
                        popDist = node.popularity{pi, pj};
                        dj = dist2json(popDist);
                        if ~isempty(dj) && pj <= K
                            popMap(classes{pj}.name) = dj;
                        end
                    end
                end
            end
            if popMap.Count > 0
                cc('popularity') = popMap;
            end
        end
 
        % Access-cost (list-move) structure: per-item graph shared by all
        % classes, or full per-class accessProb matrices. The default
        % super-diagonal is rebuilt by sanitize on load and is not saved.
        if ~isempty(node.graph)
            gArr = cell(1, numel(node.graph));
            for gi = 1:numel(node.graph)
                gArr{gi} = full(node.graph{gi});
            end
            cc('accessGraph') = gArr;
        elseif ~isempty(node.accessProb)
            [Kap, Nap] = size(node.accessProb);
            apArr = cell(1, Kap);
            for k1 = 1:Kap
                rowArr = cell(1, Nap);
                for k2 = 1:Nap
                    if ~isempty(node.accessProb{k1, k2})
                        rowArr{k2} = full(node.accessProb{k1, k2});
                    else
                        rowArr{k2} = [];
                    end
                end
                apArr{k1} = rowArr;
            end
            cc('accessProb') = apArr;
        end
 
        % Initial cache state [class counts | contents | retrieval bitmap]
        cacheState = node.getState;
        if ~isempty(cacheState)
            cc('initialState') = num2cell(double(full(cacheState(1, :))));
        end
 
        nj('cache') = cc;
 
        % Also emit flat cache fields at the node level for the Java LineModelIO
        % reader used by the LDES engine CLI (Python save_model emits the same
        % flat form; MATLAB/Python linemodel_load read the nested 'cache' object
        % above). itemLevelCap is forced to a JSON array via a cell wrapper.
        % The CLIMB -> FIFO-over-unit-lists remap applied to the nested
        % 'replacement' key above is applied here too: emitting CLIMB flat while
        % the nested key says FIFO would have the two readers rebuild different
        % caches from the same file. Both keys are taken from cc.
        nj('numItems') = double(node.items.nitems);
        nj('itemLevelCap') = num2cell(double(ilc(:)'));
        nj('replacementStrategy') = cc('replacement');
        flatKeys = {'hitClass', 'missClass', 'popularity', 'accessGraph', ...
                    'accessProb', 'initialState', 'admissionProb'};
        for fk = 1:numel(flatKeys)
            if isKey(cc, flatKeys{fk})
                nj(flatKeys{fk}) = cc(flatKeys{fk});
            end
        end
 
        % Retrieval system (delayed-hit cache): flat block for the Java
        % LineModelIO reader, mirroring Python save_model. Present only when a
        % retrieval system was configured (setRetrievalSystem). Carries, per
        % arrival class, the retrieval queue node names and the per-item
        % retrieval class into which a miss switches to be fetched.
        if ~isempty(node.retrievalSystemCapacity) && node.retrievalSystemCapacity > 0
            byClass = containers.Map();
            nItemsR = node.items.nitems;
            rc = node.server.retrievalClasses;
            qKeys = node.retrievalSystemQueueIndices.keys();
            for kk = 1:numel(qKeys)
                key0 = qKeys{kk};                 % jobinClass.index - 1 (0-based)
                inIdx = double(key0) + 1;
                if inIdx < 1 || inIdx > K
                    continue;
                end
                entry = containers.Map();
                qidxs = node.retrievalSystemQueueIndices(key0);
                qnames = cell(1, numel(qidxs));
                for qi = 1:numel(qidxs)
                    qnames{qi} = nodes{qidxs(qi)}.name;
                end
                entry('queues') = qnames;
                itemsMap = containers.Map();
                for it = 1:nItemsR
                    if size(rc, 1) >= it && size(rc, 2) >= inIdx
                        rClassIdx = rc(it, inIdx);
                        if rClassIdx > 0 && rClassIdx <= K
                            itemsMap(num2str(it - 1)) = classes{rClassIdx}.name;
                        end
                    end
                end
                if itemsMap.Count > 0
                    entry('items') = itemsMap;
                end
                byClass(classes{inIdx}.name) = entry;
            end
            if byClass.Count > 0
                rsRoot = containers.Map();
                rsRoot('capacity') = double(node.retrievalSystemCapacity);
                rsRoot('byClass') = byClass;
                nj('retrievalSystem') = rsRoot;
            end
        end
    end
 
    % Fork tasksPerLink
    if isa(node, 'Fork')
        if ~isempty(node.output) && isprop(node.output, 'tasksPerLink') && node.output.tasksPerLink > 1
            nj('tasksPerLink') = node.output.tasksPerLink;
        end
    end
 
    % Join paired fork and join strategy
    if isa(node, 'Join')
        if ~isempty(node.joinOf)
            nj('forkNode') = node.joinOf.name;
        end
        % Serialize per-class join strategy if non-default
        if ~isempty(node.input) && isprop(node.input, 'joinStrategy') && ~isempty(node.input.joinStrategy)
            for r = 1:K
                jc = classes{r};
                if r <= length(node.input.joinStrategy) && ~isempty(node.input.joinStrategy{r})
                    js = node.input.joinStrategy{r};
                    if js ~= JoinStrategy.STD
                        if js == JoinStrategy.PARTIAL
                            nj('joinStrategy') = 'PARTIAL';
                        end
                    end
                end
            end
        end
        if ~isempty(node.input) && isprop(node.input, 'joinRequired') && ~isempty(node.input.joinRequired)
            for r = 1:K
                jc = classes{r};
                if r <= length(node.input.joinRequired) && ~isempty(node.input.joinRequired{r})
                    jq = node.input.joinRequired{r};
                    if jq > 0
                        nj('joinQuorum') = jq;
                    end
                end
            end
        end
    end
 
    % DPS/GPS scheduling parameters (per-class weights). Include the priority
    % variants DPSPRIO/GPSPRIO: they carry the same per-class weights via
    % schedStrategyPar, and omitting them silently resets the weights to 1 on
    % reload (e.g. prio_identical's GPSPRIO weights).
    if isa(node, 'Queue') && ~isa(node, 'Delay')
        sched = node.schedStrategy;
        if ~isempty(sched) && (sched == SchedStrategy.DPS || sched == SchedStrategy.GPS || ...
                sched == SchedStrategy.DPSPRIO || sched == SchedStrategy.GPSPRIO)
            sp = containers.Map();
            for r = 1:K
                jc = classes{r};
                try
                    w = node.schedStrategyPar(r);
                    if ~isempty(w) && isfinite(w) && w > 0
                        sp(jc.name) = w;
                    end
                catch
                end
            end
            if sp.Count > 0
                nj('schedParams') = sp;
            end
        end
    end
 
    % Transition modes
    if isa(node, 'Transition')
        modesJson = {};
        nModes = node.getNumberOfModes();
        allNodes = model.getNodes();
        for mi = 1:nModes
            mj = containers.Map();
            if mi <= length(node.modeNames) && ~isempty(node.modeNames{mi})
                mj('name') = node.modeNames{mi};
            else
                mj('name') = sprintf('Mode%d', mi);
            end
            % Distribution. An immediate mode fires with no delay and has no
            % firing distribution: addMode leaves an Exp(1) placeholder in
            % distributions{mi}, which the timing strategy overrides. Emitting
            % that placeholder would present the mode to a reader as a timed one
            % firing at rate 1, so it is omitted here, as the Java and Python
            % serializers do.
            isImmediateMode = mi <= length(node.timingStrategies) && ...
                node.timingStrategies(mi) == TimingStrategy.IMMEDIATE;
            if ~isImmediateMode && mi <= length(node.distributions) && ...
                    ~isempty(node.distributions{mi})
                dj = dist2json(node.distributions{mi});
                if ~isempty(dj)
                    mj('distribution') = dj;
                end
            end
            % Timing strategy
            if mi <= length(node.timingStrategies)
                if node.timingStrategies(mi) == TimingStrategy.TIMED
                    mj('timingStrategy') = 'TIMED';
                else
                    mj('timingStrategy') = 'IMMEDIATE';
                end
            end
            % Number of servers
            if mi <= length(node.numberOfServers) && node.numberOfServers(mi) > 1
                mj('numServers') = node.numberOfServers(mi);
            end
            % Firing priority
            if mi <= length(node.firingPriorities) && node.firingPriorities(mi) > 0
                mj('firingPriority') = node.firingPriorities(mi);
            end
            % Firing weight
            if mi <= length(node.firingWeights) && node.firingWeights(mi) ~= 1.0
                mj('firingWeight') = node.firingWeights(mi);
            end
            % Enabling conditions
            if mi <= length(node.enablingConditions)
                ecMat = node.enablingConditions{mi};
                ecList = {};
                for ni = 1:size(ecMat, 1)
                    for ci = 1:size(ecMat, 2)
                        if ecMat(ni, ci) > 0
                            ec = containers.Map();
                            ec('node') = allNodes{ni}.name;
                            ec('class') = classes{ci}.name;
                            ec('count') = ecMat(ni, ci);
                            ecList{end+1} = ec; %#ok<AGROW>
                        end
                    end
                end
                if ~isempty(ecList)
                    mj('enablingConditions') = ecList;
                end
            end
            % Inhibiting conditions
            if mi <= length(node.inhibitingConditions)
                icMat = node.inhibitingConditions{mi};
                icList = {};
                for ni = 1:size(icMat, 1)
                    for ci = 1:size(icMat, 2)
                        if isfinite(icMat(ni, ci))
                            ic = containers.Map();
                            ic('node') = allNodes{ni}.name;
                            ic('class') = classes{ci}.name;
                            ic('count') = icMat(ni, ci);
                            icList{end+1} = ic; %#ok<AGROW>
                        end
                    end
                end
                if ~isempty(icList)
                    mj('inhibitingConditions') = icList;
                end
            end
            % Firing outcomes
            if mi <= length(node.firingOutcomes)
                foMat = node.firingOutcomes{mi};
                foList = {};
                for ni = 1:size(foMat, 1)
                    for ci = 1:size(foMat, 2)
                        if foMat(ni, ci) ~= 0
                            fo = containers.Map();
                            fo('node') = allNodes{ni}.name;
                            fo('class') = classes{ci}.name;
                            fo('count') = foMat(ni, ci);
                            foList{end+1} = fo; %#ok<AGROW>
                        end
                    end
                end
                if ~isempty(foList)
                    mj('firingOutcomes') = foList;
                end
            end
            modesJson{end+1} = mj; %#ok<AGROW>
        end
        if ~isempty(modesJson)
            nj('modes') = modesJson;
        end
    end
 
    % Initial state for Place nodes (token counts). Wrapped in a cell so that a
    % single-class marking still encodes as a JSON array: the readers index it
    % positionally, and a bare scalar would not survive the round-trip.
    if isa(node, 'Place') && ~isempty(node.state)
        nj('initialState') = num2cell(double(node.state(:)'));
    end
 
    % Prior over the node's initial states (setStatePrior). Emitted here, in the
    % node loop, rather than in a later pass: the routing section below calls
    % model.getStruct(), whose initialization overwrites every statePrior with
    % the default, so a later pass would serialize that default instead of the
    % user's prior.
    % The prior indexes the rows of the node's state space, so it is meaningless
    % without it: the reader validates the two against each other and rejects a
    % mismatch. Emit the pair or neither, as the JAR writer does.
    if isa(node, 'StatefulNode') && ~isempty(node.statePrior)
        prior = double(node.statePrior(:));
        % The trivial prior [1] over a single state is exactly what
        % initDefault/initFromMarginal rebuild from "initialState", so it is not
        % emitted.
        trivialPrior = numel(prior) == 1 && abs(prior(1) - 1.0) < 1e-12;
        if ~trivialPrior
            space = double(full(node.space));
            if isempty(space) || size(space, 1) ~= numel(prior)
                line_warning(mfilename, sprintf(['Node %s carries a state prior over %d states but a ' ...
                    'state space of %d rows; the prior is not saved.'], ...
                    node.getName(), numel(prior), size(space, 1)));
            else
                spaceRows = cell(1, size(space, 1));
                for si = 1:size(space, 1)
                    spaceRows{si} = num2cell(space(si, :));
                end
                nj('stateSpace') = spaceRows;
                % num2cell keeps a single-state (1x1) prior from collapsing to a
                % bare scalar: the reader takes statePrior as an array
                % (getAsJsonArray), so a scalar makes it throw.
                nj('statePrior') = num2cell(prior');
            end
        end
    end
 
    nodesJson{end+1} = nj; %#ok<AGROW>
end
result('nodes') = nodesJson;
 
% --- Classes ---
classesJson = {};
for r = 1:K
    jc = classes{r};
    cj = containers.Map();
    cj('name') = jc.name;
    if isa(jc, 'OpenSignal')
        cj('type') = 'Signal';
        cj('openOrClosed') = 'Open';
        cj('signalType') = SignalType.toText(jc.signalType);
        if ~isempty(jc.targetJobClass)
            cj('targetClass') = jc.targetJobClass.name;
        end
        if ~isempty(jc.removalDistribution)
            cj('removalDistribution') = dist2json(jc.removalDistribution);
        end
        if ~isempty(jc.removalPolicy) && jc.removalPolicy ~= RemovalPolicy.RANDOM
            cj('removalPolicy') = RemovalPolicy.toText(jc.removalPolicy);
        end
    elseif isa(jc, 'ClosedSignal')
        cj('type') = 'Signal';
        cj('openOrClosed') = 'Closed';
        cj('signalType') = SignalType.toText(jc.signalType);
        if ~isempty(jc.refstat) && isprop(jc.refstat, 'name')
            cj('refNode') = jc.refstat.name;
        end
        if ~isempty(jc.targetJobClass)
            cj('targetClass') = jc.targetJobClass.name;
        end
        if ~isempty(jc.removalDistribution)
            cj('removalDistribution') = dist2json(jc.removalDistribution);
        end
        if ~isempty(jc.removalPolicy) && jc.removalPolicy ~= RemovalPolicy.RANDOM
            cj('removalPolicy') = RemovalPolicy.toText(jc.removalPolicy);
        end
    elseif isa(jc, 'Signal')
        % A bare Signal subclasses JobClass directly (OpenSignal/ClosedSignal
        % subclass OpenClass/ClosedClass instead), so it matched neither of the
        % branches above nor the OpenClass branch below and was emitted as a
        % plain "Open" class, losing its signal semantics. It is neither open nor
        % closed, so 'openOrClosed' is omitted; that absence is what tells the
        % reader to rebuild a Signal rather than an OpenSignal/ClosedSignal.
        cj('type') = 'Signal';
        cj('signalType') = SignalType.toText(jc.signalType);
        if ~isempty(jc.targetJobClass)
            cj('targetClass') = jc.targetJobClass.name;
        end
        if ~isempty(jc.removalDistribution)
            cj('removalDistribution') = dist2json(jc.removalDistribution);
        end
        if ~isempty(jc.removalPolicy) && jc.removalPolicy ~= RemovalPolicy.RANDOM
            cj('removalPolicy') = RemovalPolicy.toText(jc.removalPolicy);
        end
    elseif isa(jc, 'SelfLoopingClass')
        % SelfLoopingClass subclasses ClosedClass, so it must be tested BEFORE
        % Closed: otherwise the Closed branch shadows it and the class reloads as
        % an ordinary closed class that no longer self-loops.
        cj('type') = 'SelfLooping';
        cj('population') = jc.population;
        if ~isempty(jc.refstat) && isprop(jc.refstat, 'name')
            cj('refNode') = jc.refstat.name;
        end
    elseif isa(jc, 'ClosedClass')
        cj('type') = 'Closed';
        cj('population') = jc.population;
        if ~isempty(jc.refstat) && isprop(jc.refstat, 'name')
            cj('refNode') = jc.refstat.name;
        end
    elseif isa(jc, 'OpenClass')
        cj('type') = 'Open';
        % An open class's reference station is normally the Source and is
        % re-derived on load, but setReferenceStation may have overridden it.
        % That override is real state, so carry it.
        if ~isempty(jc.refstat) && isprop(jc.refstat, 'name') && ~isa(jc.refstat, 'Source')
            cj('refNode') = jc.refstat.name;
        end
    else
        cj('type') = 'Open';
    end
    if jc.priority ~= 0
        cj('priority') = jc.priority;
    end
    if isprop(jc, 'deadline') && isfinite(jc.deadline)
        cj('deadline') = jc.deadline;
    end
    if jc.isReferenceClass()
        cj('isReferenceClass') = true;
    end
    % Reply signal binding (sn.syncreply). Without it a REPLY signal class is
    % inert after a round-trip: nothing unblocks the servers waiting on it.
    if isprop(jc, 'replySignalClass') && ~isempty(jc.replySignalClass)
        cj('replySignalClass') = jc.replySignalClass.name;
    end
    % Spawn-on-completion binding (sn.classspawn): the class injected at the
    % same station whenever a job of this class completes service.
    if isprop(jc, 'spawnClass') && ~isempty(jc.spawnClass)
        cj('spawnClass') = jc.spawnClass.name;
    end
    % Class-level (global) patience, distinct from the node-scoped 'patience'
    % emitted per Queue. A node-scoped entry overrides this one on load.
    if isprop(jc, 'patience') && ~isempty(jc.patience) && ~isa(jc.patience, 'Disabled')
        cj('patience') = dist2json(jc.patience);
        if ~isempty(jc.impatienceType)
            cj('impatienceType') = ImpatienceType.toText(jc.impatienceType);
        end
    end
    classesJson{end+1} = cj; %#ok<AGROW>
end
result('classes') = classesJson;
 
% --- Routing ---
routingMap = containers.Map();
try
    sn = model.getStruct();
    % Prefer rtorig (original P matrix before ClassSwitch expansion)
    if ~isempty(sn) && isfield(sn, 'rtorig') && iscell(sn.rtorig) && ~isempty(sn.rtorig) && ~isempty(sn.rtorig{1,1})
        P_orig = sn.rtorig;
        M_orig = size(P_orig{1,1}, 1);
        % Identify explicit ClassSwitch nodes (not auto-added) whose switch is
        % carried by a NON-IDENTITY classSwitchMatrix. Only those need same-class
        % collapsing: the matrix is reapplied on load, so emitting cross-class
        % routing too would double-switch. An explicit ClassSwitch with an
        % identity matrix expresses the switch through the routing itself, which
        % must be emitted verbatim to preserve it (mirrors jline.io.LineModelIO).
        nodes = model.getNodes();
        explicit_cs = false(1, M_orig);
        for ii = 1:min(M_orig, length(nodes))
            if isa(nodes{ii}, 'ClassSwitch') && ~nodes{ii}.autoAdded ...
                    && ~isIdentityClassSwitch(nodes{ii}.server.csMatrix, K)
                explicit_cs(ii) = true;
            end
        end
        % For explicit CS sources, compute same-class routing:
        % P_same(s,ii,jj) = sum_r P_orig{r,s}(ii,jj)
        % This avoids saving cross-class entries that would cause
        % double-switching on load.
        cs_same = zeros(K, M_orig, M_orig);
        for ii = 1:M_orig
            if explicit_cs(ii)
                for s = 1:K
                    for jj = 1:M_orig
                        total = 0;
                        for r = 1:K
                            Prs = P_orig{r,s};
                            if issparse(Prs); Prs = full(Prs); end
                            total = total + Prs(ii, jj);
                        end
                        cs_same(s, ii, jj) = total;
                    end
                end
            end
        end
        for r = 1:K
            for s = 1:K
                fromTo = containers.Map();
                Prs = P_orig{r,s};
                if issparse(Prs)
                    Prs = full(Prs);
                end
                for ii = 1:M_orig
                    if explicit_cs(ii)
                        % For explicit CS, use same-class routing only
                        if r == s
                            for jj = 1:M_orig
                                val = cs_same(s, ii, jj);
                                if val > 1e-14
                                    ni = sn.nodenames{ii};
                                    njn = sn.nodenames{jj};
                                    if ~fromTo.isKey(ni)
                                        fromTo(ni) = containers.Map();
                                    end
                                    dest = fromTo(ni);
                                    dest(njn) = val;
                                    fromTo(ni) = dest;
                                end
                            end
                        end
                        % Skip cross-class entries from explicit CS
                        continue;
                    end
                    for jj = 1:M_orig
                        val = Prs(ii, jj);
                        if val > 1e-14
                            ni = sn.nodenames{ii};
                            njn = sn.nodenames{jj};
                            if ~fromTo.isKey(ni)
                                fromTo(ni) = containers.Map();
                            end
                            dest = fromTo(ni);
                            dest(njn) = val;
                            fromTo(ni) = dest;
                        end
                    end
                end
                if fromTo.Count > 0
                    key = char(sprintf('%s,%s', classes{r}.name, classes{s}.name));
                    routingMap(key) = fromTo;
                end
            end
        end
    elseif ~isempty(sn) && isfield(sn, 'rtnodes') && ~isempty(sn.rtnodes)
        % Fallback to rtnodes if rtorig not available
        rt = sn.rtnodes;
        N = sn.nnodes;
        % Identify explicit ClassSwitch node indices (not auto-added). rtnodes
        % folds the CS switching into cross-class entries; emitting those
        % together with the node's classSwitchMatrix double-encodes the switch
        % and double-switches classes on load. Mirror the rtorig branch: for
        % explicit CS sources emit only same-class topology entries.
        nodes = model.getNodes();
        explicit_cs = false(1, N);
        for ii = 1:min(N, length(nodes))
            if isa(nodes{ii}, 'ClassSwitch') && ~nodes{ii}.autoAdded ...
                    && ~isIdentityClassSwitch(nodes{ii}.server.csMatrix, K)
                explicit_cs(ii) = true;
            end
        end
        % cs_same(s,ii,jj) = sum_r rt((ii,r),(jj,s)): destination probability
        % conditioned on the class s the job leaves the CS in (same formula as
        % the rtorig branch).
        cs_same = zeros(K, N, N);
        for ii = 1:N
            if explicit_cs(ii)
                for s = 1:K
                    for jj = 1:N
                        total = 0;
                        for r = 1:K
                            total = total + rt((ii-1)*K+r, (jj-1)*K+s);
                        end
                        cs_same(s, ii, jj) = total;
                    end
                end
            end
        end
        for r = 1:K
            for s = 1:K
                fromTo = containers.Map();
                for ii = 1:N
                    if explicit_cs(ii)
                        % For explicit CS, use same-class routing only
                        if r == s
                            for jj = 1:N
                                val = cs_same(s, ii, jj);
                                if val > 1e-14
                                    ni = sn.nodenames{ii};
                                    njn = sn.nodenames{jj};
                                    if ~fromTo.isKey(ni)
                                        fromTo(ni) = containers.Map();
                                    end
                                    dest = fromTo(ni);
                                    dest(njn) = val;
                                    fromTo(ni) = dest;
                                end
                            end
                        end
                        % Skip cross-class entries from explicit CS
                        continue;
                    end
                    for jj = 1:N
                        val = rt((ii-1)*K+r, (jj-1)*K+s);
                        if val > 1e-14
                            ni = sn.nodenames{ii};
                            njn = sn.nodenames{jj};
                            if ~fromTo.isKey(ni)
                                fromTo(ni) = containers.Map();
                            end
                            dest = fromTo(ni);
                            dest(njn) = val;
                            fromTo(ni) = dest;
                        end
                    end
                end
                if fromTo.Count > 0
                    key = char(sprintf('%s,%s', classes{r}.name, classes{s}.name));
                    routingMap(key) = fromTo;
                end
            end
        end
    end
catch
    % If struct not available, routing stays empty
end
 
routing = containers.Map();
routing('type') = 'matrix';
routing('matrix') = routingMap;
result('routing') = routing;
 
% --- Routing Strategies ---
try
    sn2 = model.getStruct();
    if ~isempty(sn2) && isfield(sn2, 'routing') && ~isempty(sn2.routing)
        routingStrategies = containers.Map();
        stratNames = containers.Map('KeyType','int32','ValueType','char');
        stratNames(int32(RoutingStrategy.RAND)) = 'RAND';
        stratNames(int32(RoutingStrategy.RROBIN)) = 'RROBIN';
        stratNames(int32(RoutingStrategy.WRROBIN)) = 'WRROBIN';
        stratNames(int32(RoutingStrategy.JSQ)) = 'JSQ';
        stratNames(int32(RoutingStrategy.KCHOICES)) = 'KCHOICES';
        stratNames(int32(RoutingStrategy.FIRING)) = 'FIRING';
        stratNames(int32(RoutingStrategy.RL)) = 'RL';
        stratNames(int32(RoutingStrategy.DISABLED)) = 'DISABLED';
        for i = 1:sn2.nnodes
            nodeStrats = containers.Map();
            for r = 1:K
                routVal = int32(sn2.routing(i, r));
                if routVal ~= int32(RoutingStrategy.PROB) && routVal ~= int32(RoutingStrategy.RAND) && stratNames.isKey(routVal)
                    nodeStrats(classes{r}.name) = stratNames(routVal);
                end
            end
            if nodeStrats.Count > 0
                routingStrategies(sn2.nodenames{i}) = nodeStrats;
            end
        end
        if routingStrategies.Count > 0
            result('routingStrategies') = routingStrategies;
        end
 
        % Save WRROBIN weights. Index the node list by the NODE index i:
        % indexing it by the station index silently reads a different node
        % (and skips non-station WRROBIN nodes such as Router), losing the
        % weights on save.
        routingWeights = containers.Map();
        for i = 1:sn2.nnodes
            nodeObj2 = nodes{i};
            nodeClassWeights = containers.Map();
            for r = 1:K
                if int32(sn2.routing(i, r)) == int32(RoutingStrategy.WRROBIN)
                    os = nodeObj2.output.outputStrategy;
                    if size(os,2) >= r
                        osEntry = os{1, r};
                        % osEntry = {className, stratName, forwardLinks}
                        if length(osEntry) >= 3
                            fwdLinks = osEntry{3};
                            destWeights = containers.Map();
                            for fi = 1:length(fwdLinks)
                                link = fwdLinks{fi};
                                % link = {destNode, weight}
                                if iscell(link) && length(link) >= 2 && isa(link{1}, 'Node')
                                    destWeights(link{1}.name) = link{2};
                                end
                            end
                            if destWeights.Count > 0
                                nodeClassWeights(classes{r}.name) = destWeights;
                            end
                        end
                    end
                end
            end
            if nodeClassWeights.Count > 0
                routingWeights(sn2.nodenames{i}) = nodeClassWeights;
            end
        end
        if routingWeights.Count > 0
            result('routingWeights') = routingWeights;
        end
    end
catch
end
 
% --- Setup / Delay-Off, Polling Type and Switchover Times ---
nodesCellTmp = result('nodes');
for i = 1:M
    nodeObj = nodes{i};
    if ~isa(nodeObj, 'Queue') || isa(nodeObj, 'Delay')
        continue;
    end
    njIdx = 0;
    for nj_idx = 1:length(nodesCellTmp)
        if strcmp(nodesCellTmp{nj_idx}('name'), nodeObj.name)
            njIdx = nj_idx;
            break;
        end
    end
    if njIdx == 0
        continue;
    end
    nj = nodesCellTmp{njIdx};
 
    % Setup and delay-off are emitted as a pair: setDelayOff requires both
    % on reload, so a class with only one of the two is not representable.
    setupMap = containers.Map();
    delayOffMap = containers.Map();
    for r = 1:K
        if r <= length(nodeObj.setupTime) && r <= length(nodeObj.delayoffTime)
            suDist = nodeObj.setupTime{1,r};
            doffDist = nodeObj.delayoffTime{1,r};
            if ~isempty(suDist) && ~isempty(doffDist) && ...
                    ~isa(suDist, 'Disabled') && ~isa(doffDist, 'Disabled')
                setupMap(classes{r}.name) = dist2json(suDist);
                delayOffMap(classes{r}.name) = dist2json(doffDist);
            end
        end
    end
    if setupMap.Count > 0
        nj('setupTime') = setupMap;
        nj('delayOffTime') = delayOffMap;
    end
 
    isPolling = SchedStrategy.toId(nodeObj.schedStrategy) == SchedStrategy.POLLING;
 
    % Polling type is written by name: the ids agree with Java but Python
    % assigns them via auto().
    if isPolling && ~isempty(nodeObj.pollingType)
        ptId = PollingType.toId(nodeObj.pollingType{1,1});
        nj('pollingType') = PollingType.toName(ptId);
        if ptId == PollingType.KLIMITED && ~isempty(nodeObj.pollingPar)
            nj('pollingPar') = nodeObj.pollingPar;
        end
    end
 
    % Under POLLING the switchover is indexed by the departing class alone
    % (a 1xK cell) and is written without a "to" field; otherwise it is a
    % KxK cell of (from,to) pairs.
    soTimes = {};
    if ~isempty(nodeObj.switchoverTime)
        [soRows, soCols] = size(nodeObj.switchoverTime);
        if isPolling
            for r = 1:min(K, soCols)
                dist = nodeObj.switchoverTime{1,r};
                if ~isempty(dist) && ~isa(dist, 'Disabled')
                    so = containers.Map();
                    so('from') = classes{r}.name;
                    so('distribution') = dist2json(dist);
                    soTimes{end+1} = so;
                end
            end
        else
            for r = 1:min(K, soRows)
                for s = 1:min(K, soCols)
                    dist = nodeObj.switchoverTime{r,s};
                    if ~isempty(dist) && ~isa(dist, 'Disabled')
                        so = containers.Map();
                        so('from') = classes{r}.name;
                        so('to') = classes{s}.name;
                        so('distribution') = dist2json(dist);
                        soTimes{end+1} = so;
                    end
                end
            end
        end
    end
    if ~isempty(soTimes)
        nj('switchoverTimes') = soTimes;
    end
    nodesCellTmp{njIdx} = nj;
end
result('nodes') = nodesCellTmp;
 
% --- Heterogeneous Server Types ---
try
    nodesCellTmp = result('nodes');
    for i = 1:M
        nodeObj = nodes{i};
        if isa(nodeObj, 'Queue') && nodeObj.isHeterogeneous()
            stArr = {};
            for ti = 1:length(nodeObj.serverTypes)
                st = nodeObj.serverTypes{ti};
                stj = containers.Map();
                stj('name') = st.name;
                stj('count') = st.numOfServers;
                % Compatible classes
                ccNames = {};
                for cci = 1:length(st.compatibleClasses)
                    ccNames{end+1} = st.compatibleClasses{cci}.name; %#ok<AGROW>
                end
                if ~isempty(ccNames)
                    stj('compatibleClasses') = ccNames;
                end
                % Per-class service distributions
                svcMap = containers.Map();
                for r = 1:K
                    jc = classes{r};
                    dist = nodeObj.getHeteroService(jc, st);
                    if ~isempty(dist) && ~isa(dist, 'Disabled')
                        svcMap(jc.name) = dist2json(dist);
                    end
                end
                if svcMap.Count > 0
                    stj('service') = svcMap;
                end
                stArr{end+1} = stj; %#ok<AGROW>
            end
            if ~isempty(stArr)
                for nj_idx = 1:length(nodesCellTmp)
                    nj = nodesCellTmp{nj_idx};
                    if strcmp(nj('name'), nodeObj.name)
                        nj('serverTypes') = stArr;
                        % Scheduling policy
                        policy = nodeObj.getHeteroSchedPolicy();
                        if ~isempty(policy) && policy ~= HeteroSchedPolicy.ORDER
                            nj('heteroSchedPolicy') = HeteroSchedPolicy.toText(policy);
                        end
                        nodesCellTmp{nj_idx} = nj;
                        break;
                    end
                end
            end
        end
    end
    result('nodes') = nodesCellTmp;
catch
end
 
% --- Balking, Retrial, Patience, Orbit Impatience, Immediate Feedback ---
% No try/catch here: a bare catch silently dropped this entire block (every
% balking threshold, retrial delay and patience distribution in the model) on any
% error, including the node_map lookup below, which never existed as a variable
% in this function at all.
nodesCellTmp = result('nodes');
nodeByName = containers.Map();
for i = 1:M
    nodeByName(nodes{i}.name) = nodes{i};
end
for nj_idx = 1:length(nodesCellTmp)
    nj = nodesCellTmp{nj_idx};
    nodeName = nj('name');
    nodeObj = nodeByName(nodeName);
    % Immediate feedback is a per-class node property on any Station.
    if isa(nodeObj, 'Queue')
        ifMap = containers.Map();
        for r = 1:K
            if nodeObj.hasImmediateFeedback(classes{r})
                ifMap(classes{r}.name) = true;
            end
        end
        if ifMap.Count > 0
            nj('immediateFeedback') = ifMap;
        end
    end
    if isa(nodeObj, 'Queue')
        % Balking
        balkJson = containers.Map();
        for r = 1:K
            jc = classes{r};
            if nodeObj.hasBalking(jc)
                [strategy, thresholds] = nodeObj.getBalking(jc);
                bjc = containers.Map();
                switch strategy
                    case BalkingStrategy.QUEUE_LENGTH, bjc('strategy') = 'QUEUE_LENGTH';
                    case BalkingStrategy.EXPECTED_WAIT, bjc('strategy') = 'EXPECTED_WAIT';
                    case BalkingStrategy.COMBINED, bjc('strategy') = 'COMBINED';
                end
                thArr = {};
                for ti = 1:length(thresholds)
                    th = thresholds{ti};
                    tjson = containers.Map();
                    tjson('minJobs') = th{1};
                    if isinf(th{2})
                        tjson('maxJobs') = -1;
                    else
                        tjson('maxJobs') = th{2};
                    end
                    tjson('probability') = th{3};
                    thArr{end+1} = tjson;
                end
                bjc('thresholds') = thArr;
                balkJson(jc.name) = bjc;
            end
        end
        if balkJson.Count > 0
            nj('balking') = balkJson;
        end
        % Retrial
        retrialJson = containers.Map();
        for r = 1:K
            jc = classes{r};
            if nodeObj.hasRetrial(jc)
                [delayDist, maxAttempts] = nodeObj.getRetrial(jc);
                rjc = containers.Map();
                rjc('delay') = dist2json(delayDist);
                rjc('maxAttempts') = maxAttempts;
                retrialJson(jc.name) = rjc;
            end
        end
        if retrialJson.Count > 0
            nj('retrial') = retrialJson;
        end
        % Patience
        patienceJson = containers.Map();
        for r = 1:K
            jc = classes{r};
            patDist = nodeObj.getPatience(jc);
            if ~isempty(patDist) && ~isa(patDist, 'Disabled')
                pjc = containers.Map();
                pjc('distribution') = dist2json(patDist);
                impType = nodeObj.getImpatienceType(jc);
                if ~isempty(impType)
                    pjc('impatienceType') = ImpatienceType.toText(impType);
                end
                patienceJson(jc.name) = pjc;
            end
        end
        if patienceJson.Count > 0
            nj('patience') = patienceJson;
        end
        % Orbit impatience (abandonment from the retrial orbit), distinct from
        % the queue patience above.
        orbitJson = containers.Map();
        for r = 1:K
            jc = classes{r};
            orbDist = nodeObj.getOrbitImpatience(jc);
            if ~isempty(orbDist) && ~isa(orbDist, 'Disabled')
                orbitJson(jc.name) = dist2json(orbDist);
            end
        end
        if orbitJson.Count > 0
            nj('orbitImpatience') = orbitJson;
        end
        % Batch rejection probability (retrial queues), per class
        brpJson = containers.Map();
        for r = 1:K
            jc = classes{r};
            brp = nodeObj.getBatchRejectProbability(jc);
            if ~isempty(brp) && brp > 0
                brpJson(jc.name) = brp;
            end
        end
        if brpJson.Count > 0
            nj('batchRejectProb') = brpJson;
        end
    end
    nodesCellTmp{nj_idx} = nj;
end
result('nodes') = nodesCellTmp;
 
% --- Finite Capacity Regions ---
try
    regions = model.regions;
    if ~isempty(regions)
        fcrArray = {};
        for ri = 1:length(regions)
            reg = regions{ri};
            rj = containers.Map();
            rj('name') = reg.name;
            % Stations with per-class details
            stationsJson = {};
            for ni = 1:length(reg.nodes)
                sj = containers.Map();
                sj('node') = reg.nodes{ni}.name;
                % Per-class classCap
                if isprop(reg, 'classMaxJobs') && ~isempty(reg.classMaxJobs)
                    ccMap = containers.Map();
                    for r = 1:K
                        jc = classes{r};
                        if r <= length(reg.classMaxJobs) && isfinite(reg.classMaxJobs(r))
                            ccMap(jc.name) = reg.classMaxJobs(r);
                        end
                    end
                    if ccMap.Count > 0
                        sj('classCap') = ccMap;
                    end
                end
                % Per-class classWeight
                if isprop(reg, 'classWeight') && ~isempty(reg.classWeight)
                    cwMap = containers.Map();
                    for r = 1:K
                        jc = classes{r};
                        if r <= length(reg.classWeight) && reg.classWeight(r) ~= 1
                            cwMap(jc.name) = reg.classWeight(r);
                        end
                    end
                    if cwMap.Count > 0
                        sj('classWeight') = cwMap;
                    end
                end
                % Per-class classSize
                if isprop(reg, 'classSize') && ~isempty(reg.classSize)
                    csMap = containers.Map();
                    for r = 1:K
                        jc = classes{r};
                        if r <= length(reg.classSize) && reg.classSize(r) ~= 1
                            csMap(jc.name) = reg.classSize(r);
                        end
                    end
                    if csMap.Count > 0
                        sj('classSize') = csMap;
                    end
                end
                stationsJson{end+1} = sj; %#ok<AGROW>
            end
            rj('stations') = stationsJson;
            if isprop(reg, 'globalMaxJobs') && isfinite(reg.globalMaxJobs)
                rj('globalMaxJobs') = reg.globalMaxJobs;
            end
            if isprop(reg, 'globalMaxMemory') && isfinite(reg.globalMaxMemory)
                rj('globalMaxMemory') = reg.globalMaxMemory;
            end
            % Per-class classMaxJobs at region level
            if isprop(reg, 'classMaxJobs') && ~isempty(reg.classMaxJobs)
                cmjMap = containers.Map();
                for r = 1:K
                    jc = classes{r};
                    if r <= length(reg.classMaxJobs) && isfinite(reg.classMaxJobs(r))
                        cmjMap(jc.name) = reg.classMaxJobs(r);
                    end
                end
                if cmjMap.Count > 0
                    rj('classMaxJobs') = cmjMap;
                end
            end
            % Per-class classMaxMemory at region level. The memory budget pairs with
            % the per-station classSize footprint and folds into the equivalent job
            % cap floor(maxMem_r/size_r); dropping it silently leaves the class
            % unconstrained on the reader side.
            if isprop(reg, 'classMaxMemory') && ~isempty(reg.classMaxMemory)
                cmmMap = containers.Map();
                for r = 1:K
                    jc = classes{r};
                    if r <= length(reg.classMaxMemory) && isfinite(reg.classMaxMemory(r)) ...
                            && reg.classMaxMemory(r) >= 0
                        cmmMap(jc.name) = reg.classMaxMemory(r);
                    end
                end
                if cmmMap.Count > 0
                    rj('classMaxMemory') = cmmMap;
                end
            end
            % Drop rule
            if isprop(reg, 'dropRule') && ~isempty(reg.dropRule)
                drMap = containers.Map();
                for r = 1:K
                    jc = classes{r};
                    if r <= length(reg.dropRule)
                        drStr = droprule_to_str(reg.dropRule(r));
                        if ~isempty(drStr)
                            drMap(jc.name) = drStr;
                        end
                    end
                end
                if drMap.Count > 0
                    rj('dropRule') = drMap;
                end
            end
            % Linear constraints (A * x <= b) if present
            if ismethod(reg, 'hasLinearConstraints') && reg.hasLinearConstraints()
                [A, b] = reg.getLinearConstraints();
                if ~isempty(A) && ~isempty(b)
                    % Serialize A row-by-row as cell of arrays for JSON compat
                    Acell = cell(1, size(A,1));
                    for ri = 1:size(A,1)
                        Acell{ri} = A(ri,:);
                    end
                    rj('constraintA') = Acell;
                    rj('constraintB') = b(:)';
                end
            end
            fcrArray{end+1} = rj; %#ok<AGROW>
        end
        if ~isempty(fcrArray)
            result('finiteCapacityRegions') = fcrArray;
        end
    end
catch
end
 
% --- Rewards ---
rewardsJson = rewards2json(model);
if ~isempty(rewardsJson)
    result('rewards') = rewardsJson;
end
end
 
 
% =========================================================================
%  Reward serialization
% =========================================================================
 
function rewardsJson = rewards2json(model)
% Serialize the model's reward definitions in the declarative form
%   {name, type, node, class}
% Only rewards created through a Reward.* template carry the structural
% metadata needed to reproduce them. A reward defined from a bare function
% handle (or via Reward.custom) is not reproducible from JSON: warn and omit
% it rather than emit a reward that would be wrong on reload.
rewardsJson = {};
sn = model.sn;
if isempty(sn) || ~isfield(sn, 'reward') || isempty(sn.reward)
    return;
end
% Emit in name order, so that the array is identical to the one written by the
% Python and JAR writers (the JAR stores its rewards in a HashMap, whose iteration
% order is not insertion order).
rewardNames = cell(1, length(sn.reward));
for i = 1:length(sn.reward)
    rewardNames{i} = sn.reward{i}.name;
end
[~, order] = sort(rewardNames);
for oi = 1:length(order)
    rw = sn.reward{order(oi)};
    descriptor = [];
    if isfield(rw, 'descriptor')
        descriptor = rw.descriptor;
    end
    if isempty(descriptor) || ~isa(descriptor, 'RewardDescriptor')
        line_warning(mfilename, sprintf(['Reward "%s" is defined by a bare function handle and cannot be ' ...
            'serialized to JSON; it is omitted from the saved model. Use a Reward.* template ' ...
            '(Reward.queueLength/utilization/blocking) for a serializable reward.'], rw.name));
        continue;
    end
    if strcmp(descriptor.kind, 'Custom')
        line_warning(mfilename, sprintf(['Reward "%s" is a custom reward wrapping an arbitrary function and ' ...
            'cannot be serialized to JSON; it is omitted from the saved model.'], rw.name));
        continue;
    end
    rj = containers.Map();
    rj('name') = rw.name;
    rj('type') = descriptor.kind;
    if isempty(descriptor.node)
        line_warning(mfilename, sprintf(['Reward "%s" of type %s has no associated node and cannot be ' ...
            'serialized to JSON; it is omitted from the saved model.'], rw.name, descriptor.kind));
        continue;
    end
    rj('node') = descriptor.node.name;
    if ~isempty(descriptor.jobclass)
        rj('class') = descriptor.jobclass.name;
    end
    rewardsJson{end+1} = rj; %#ok<AGROW>
end
end
 
 
% =========================================================================
%  LayeredNetwork serialization
% =========================================================================
 
function result = layered2json(model)
result = containers.Map();
result('type') = 'LayeredNetwork';
result('name') = model.getName();
 
% --- Processors ---
procsJson = {};
hosts = model.hosts;
for i = 1:length(hosts)
    h = hosts{i};
    pj = containers.Map();
    pj('name') = h.name;
    mult = h.multiplicity;
    if ~isfinite(mult)
        pj('multiplicity') = inf_multiplicity();
    elseif mult > 1
        pj('multiplicity') = mult;
    end
    schedStr = h.scheduling;
    if ~isempty(schedStr) && ~strcmpi(schedStr, 'inf')
        pj('scheduling') = upper(schedStr);
    end
    q = h.quantum;
    if q > 0 && q ~= 0.001
        pj('quantum') = q;
    end
    sf = h.speedFactor;
    if sf ~= 1.0
        pj('speedFactor') = sf;
    end
    repl = h.replication;
    if repl > 1
        pj('replication') = repl;
    end
    procsJson{end+1} = pj; %#ok<AGROW>
end
result('hosts') = procsJson;
 
% --- Tasks ---
tasksJson = {};
tasksList = model.tasks;
for i = 1:length(tasksList)
    t = tasksList{i};
    tj = containers.Map();
    tj('name') = t.name;
    if ~isempty(t.parent)
        tj('host') = t.parent.name;
    end
    mult = t.multiplicity;
    if ~isfinite(mult)
        tj('multiplicity') = inf_multiplicity();
    elseif mult > 1
        tj('multiplicity') = mult;
    end
    schedStr = t.scheduling;
    if ~isempty(schedStr)
        tj('scheduling') = upper(schedStr);
    end
    % Think time
    ttMean = t.thinkTimeMean;
    if ~isempty(ttMean) && ttMean > GlobalConstants.FineTol
        if ~isempty(t.thinkTime) && isa(t.thinkTime, 'Distribution')
            tj('thinkTime') = dist2json(t.thinkTime);
        else
            params = containers.Map();
            params('lambda') = 1.0 / ttMean;
            dj = containers.Map();
            dj('type') = 'Exp';
            dj('params') = params;
            tj('thinkTime') = dj;
        end
    end
    % Fan in
    if ~isempty(t.fanInSource) && ischar(t.fanInSource) && ~isempty(t.fanInSource)
        fi = containers.Map();
        fi(t.fanInSource) = t.fanInValue;
        tj('fanIn') = fi;
    end
    % Fan out
    if ~isempty(t.fanOutDest)
        fo = containers.Map();
        for fi_idx = 1:length(t.fanOutDest)
            fo(t.fanOutDest{fi_idx}) = t.fanOutValue(fi_idx);
        end
        tj('fanOut') = fo;
    end
    repl = t.replication;
    if repl > 1
        tj('replication') = repl;
    end
    % FunctionTask detection
    if isa(t, 'FunctionTask')
        tj('taskType') = 'FunctionTask';
    end
    % Setup time / delay-off time (on any Task)
    if ~isempty(t.setupTime) && isa(t.setupTime, 'Distribution')
        stMean = t.setupTimeMean;
        if stMean > GlobalConstants.FineTol
            tj('setupTime') = dist2json(t.setupTime);
        end
    end
    if ~isempty(t.delayOffTime) && isa(t.delayOffTime, 'Distribution')
        dotMean = t.delayOffTimeMean;
        if dotMean > GlobalConstants.FineTol
            tj('delayOffTime') = dist2json(t.delayOffTime);
        end
    end
    % CacheTask detection
    if isa(t, 'CacheTask')
        tj('taskType') = 'CacheTask';
        tj('totalItems') = t.items;
        tj('cacheCapacity') = t.itemLevelCap;
        rs = t.replacestrategy;
        rsNameMap = containers.Map({ReplacementStrategy.RR, ReplacementStrategy.FIFO, ...
            ReplacementStrategy.SFIFO, ReplacementStrategy.LRU}, ...
            {'RR', 'FIFO', 'SFIFO', 'LRU'});
        if rsNameMap.isKey(rs)
            tj('replacementStrategy') = rsNameMap(rs);
        else
            tj('replacementStrategy') = 'FIFO';
        end
    end
    tasksJson{end+1} = tj; %#ok<AGROW>
end
result('tasks') = tasksJson;
 
% --- Entries ---
entriesJson = {};
entriesList = model.entries;
for i = 1:length(entriesList)
    e = entriesList{i};
    ej = containers.Map();
    ej('name') = e.name;
    if ~isempty(e.parent)
        ej('task') = e.parent.name;
    end
    % Entry arrival distribution
    if ~isempty(e.arrival) && isa(e.arrival, 'Distribution')
        ej('arrival') = dist2json(e.arrival);
    end
    % ItemEntry detection
    if isa(e, 'ItemEntry')
        ej('entryType') = 'ItemEntry';
        ej('totalItems') = e.cardinality;
        if ~isempty(e.popularity)
            if isa(e.popularity, 'Distribution')
                ej('accessProb') = dist2json(e.popularity);
            end
        end
    end
    entriesJson{end+1} = ej; %#ok<AGROW>
end
result('entries') = entriesJson;
 
% --- Build reply map: activityName -> entryName ---
replyMap = containers.Map();
for i = 1:length(entriesList)
    e = entriesList{i};
    if ~isempty(e.replyActivity)
        for j = 1:length(e.replyActivity)
            replyMap(e.replyActivity{j}) = e.name;
        end
    end
end
 
% --- Activities ---
actsJson = {};
actsList = model.activities;
for i = 1:length(actsList)
    a = actsList{i};
    aj = containers.Map();
    aj('name') = a.name;
    if ~isempty(a.parent)
        if isa(a.parent, 'Task') || isa(a.parent, 'Entry')
            aj('task') = a.parent.name;
        elseif ischar(a.parent) || isstring(a.parent)
            aj('task') = char(a.parent);
        elseif ischar(a.parentName) && ~isempty(a.parentName)
            aj('task') = a.parentName;
        end
    elseif ~isempty(a.parentName) && ischar(a.parentName)
        aj('task') = a.parentName;
    end
    % Host demand
    if ~isempty(a.hostDemand) && isa(a.hostDemand, 'Distribution')
        if ~isa(a.hostDemand, 'Immediate')
            aj('hostDemand') = dist2json(a.hostDemand);
        end
    elseif ~isempty(a.hostDemandMean) && a.hostDemandMean > GlobalConstants.FineTol
        params = containers.Map();
        params('lambda') = 1.0 / a.hostDemandMean;
        dj = containers.Map();
        dj('type') = 'Exp';
        dj('params') = params;
        aj('hostDemand') = dj;
    end
    % Bound to entry
    if ~isempty(a.boundToEntry)
        aj('boundToEntry') = a.boundToEntry;
    end
    % Replies to entry
    if replyMap.isKey(a.name)
        aj('repliesTo') = replyMap(a.name);
    end
    % Synch calls
    if ~isempty(a.syncCallDests)
        synchCalls = {};
        for j = 1:length(a.syncCallDests)
            sc = containers.Map();
            sc('dest') = a.syncCallDests{j};
            if j <= length(a.syncCallMeans) && a.syncCallMeans(j) ~= 1.0
                sc('mean') = a.syncCallMeans(j);
            end
            synchCalls{end+1} = sc; %#ok<AGROW>
        end
        aj('synchCalls') = synchCalls;
    end
    % Asynch calls
    if ~isempty(a.asyncCallDests)
        asynchCalls = {};
        for j = 1:length(a.asyncCallDests)
            ac = containers.Map();
            ac('dest') = a.asyncCallDests{j};
            if j <= length(a.asyncCallMeans) && a.asyncCallMeans(j) ~= 1.0
                ac('mean') = a.asyncCallMeans(j);
            end
            asynchCalls{end+1} = ac; %#ok<AGROW>
        end
        aj('asynchCalls') = asynchCalls;
    end
    actsJson{end+1} = aj; %#ok<AGROW>
end
result('activities') = actsJson;
 
% --- Precedences ---
precsJson = {};
for i = 1:length(tasksList)
    t = tasksList{i};
    precs = t.precedences;
    if isempty(precs), continue; end
    for j = 1:length(precs)
        p = precs(j);
        pj = containers.Map();
        pj('task') = t.name;
 
        preType  = p.preType;
        postType = p.postType;
 
        % Determine JSON precedence type and collect activity names
        if preType == ActivityPrecedenceType.PRE_SEQ && postType == ActivityPrecedenceType.POST_SEQ
            pj('type') = 'Serial';
            pj('activities') = [p.preActs, p.postActs];
        elseif preType == ActivityPrecedenceType.PRE_SEQ && postType == ActivityPrecedenceType.POST_AND
            pj('type') = 'AndFork';
            pj('activities') = [p.preActs, p.postActs];
        elseif preType == ActivityPrecedenceType.PRE_AND && postType == ActivityPrecedenceType.POST_SEQ
            pj('type') = 'AndJoin';
            pj('activities') = [p.preActs, p.postActs];
        elseif preType == ActivityPrecedenceType.PRE_SEQ && postType == ActivityPrecedenceType.POST_OR
            pj('type') = 'OrFork';
            pj('activities') = [p.preActs, p.postActs];
            if ~isempty(p.postParams)
                pj('probabilities') = p.postParams(:)';
            end
        elseif preType == ActivityPrecedenceType.PRE_OR && postType == ActivityPrecedenceType.POST_SEQ
            pj('type') = 'OrJoin';
            pj('activities') = [p.preActs, p.postActs];
        elseif postType == ActivityPrecedenceType.POST_LOOP
            pj('type') = 'Loop';
            % For Loop, preActs is the trigger, postActs is the loop body
            pj('activities') = p.postActs;
            if ~isempty(p.preActs)
                pj('preActivity') = p.preActs{1};
            end
            if ~isempty(p.postParams)
                pj('loopCount') = p.postParams(1);
            end
        elseif preType == ActivityPrecedenceType.PRE_SEQ && postType == ActivityPrecedenceType.POST_CACHE
            pj('type') = 'CacheAccess';
            pj('activities') = [p.preActs, p.postActs];
        else
            continue;
        end
        precsJson{end+1} = pj; %#ok<AGROW>
    end
end
if ~isempty(precsJson)
    result('precedences') = precsJson;
end
end
 
 
% =========================================================================
%  Workflow serialization
% =========================================================================
 
function result = workflow2json(model)
% Convert a Workflow to a containers.Map for JSON output.
result = containers.Map();
result('type') = 'Workflow';
result('name') = model.getName();
 
% --- Activities ---
actsJson = {};
acts = model.activities;
for i = 1:length(acts)
    act = acts{i};
    aj = containers.Map();
    aj('name') = act.name;
    if ~isempty(act.hostDemand) && isa(act.hostDemand, 'Distribution')
        dj = dist2json(act.hostDemand);
        if ~isempty(dj)
            aj('hostDemand') = dj;
        end
    end
    actsJson{end+1} = aj; %#ok<AGROW>
end
result('activities') = actsJson;
 
% --- Precedences ---
precsJson = {};
precs = model.precedences;
for i = 1:length(precs)
    p = precs(i);
    pj = containers.Map();
 
    % preActs
    preActsJson = {};
    for a = 1:length(p.preActs)
        preActsJson{end+1} = p.preActs{a}; %#ok<AGROW>
    end
    pj('preActs') = preActsJson;
 
    % postActs
    postActsJson = {};
    for a = 1:length(p.postActs)
        postActsJson{end+1} = p.postActs{a}; %#ok<AGROW>
    end
    pj('postActs') = postActsJson;
 
    % preType / postType - convert numeric IDs to JAR-compatible strings
    pj('preType') = prectype_to_str(p.preType);
    pj('postType') = prectype_to_str(p.postType);
 
    % preParams
    if ~isempty(p.preParams)
        pj('preParams') = p.preParams(:)';
    end
 
    % postParams
    if ~isempty(p.postParams)
        pj('postParams') = p.postParams(:)';
    end
 
    precsJson{end+1} = pj; %#ok<AGROW>
end
result('precedences') = precsJson;
end
 
 
% =========================================================================
%  Environment serialization
% =========================================================================
 
function result = environment2json(model)
% Convert an Environment to a containers.Map for JSON output.
result = containers.Map();
result('type') = 'Environment';
result('name') = model.getName();
 
E = height(model.envGraph.Nodes);
result('numStages') = E;
 
% --- Stages ---
stagesJson = {};
for e = 1:E
    sj = containers.Map();
    sj('name') = model.envGraph.Nodes.Name{e};
    % Serialize the stage's Network model
    if e <= length(model.ensemble) && ~isempty(model.ensemble{e})
        sj('model') = network2json(model.ensemble{e});
    end
    stagesJson{end+1} = sj; %#ok<AGROW>
end
result('stages') = stagesJson;
 
% --- Transitions ---
transJson = {};
for e = 1:E
    for h = 1:E
        if ~isempty(model.env) && e <= size(model.env, 1) && h <= size(model.env, 2) ...
                && ~isempty(model.env{e,h}) && ~isa(model.env{e,h}, 'Disabled')
            tj = containers.Map();
            tj('from') = e - 1;  % Convert to 0-indexed for JAR compatibility
            tj('to') = h - 1;    % Convert to 0-indexed for JAR compatibility
            dj = dist2json(model.env{e,h});
            if ~isempty(dj)
                tj('distribution') = dj;
                transJson{end+1} = tj; %#ok<AGROW>
            end
        end
    end
end
result('transitions') = transJson;
 
% --- Node failures ---
% Declarative record of the breakdown/repair macros applied through
% addNodeBreakdown/addNodeRepair. The stages and transitions above already carry
% the full structure losslessly; this record additionally carries the queue-length
% reset policies, which are function handles and are otherwise unrecoverable.
nfJson = {};
for i = 1:length(model.nodeFailures)
    nf = model.nodeFailures{i};
    nj = containers.Map();
    nj('node') = nf.node;
    bj = dist2json(nf.breakdown);
    if isempty(bj)
        line_warning(mfilename, sprintf(['Node failure on "%s" has a breakdown distribution that cannot be ' ...
            'serialized; the nodeFailures entry is omitted.'], nf.node));
        continue;
    end
    nj('breakdownRate') = bj;
    if ~isempty(nf.repair)
        rj = dist2json(nf.repair);
        if isempty(rj)
            line_warning(mfilename, sprintf(['Node failure on "%s" has a repair distribution that cannot be ' ...
                'serialized; the nodeFailures entry is omitted.'], nf.node));
            continue;
        end
        nj('repairRate') = rj;
    end
    dj = dist2json(nf.downService);
    if isempty(dj)
        line_warning(mfilename, sprintf(['Node failure on "%s" has a down-service distribution that cannot ' ...
            'be serialized; the nodeFailures entry is omitted.'], nf.node));
        continue;
    end
    nj('downService') = dj;
    if strcmp(nf.breakdownResetPolicy, 'custom')
        line_warning(mfilename, sprintf(['Node failure on "%s" uses a custom breakdown reset function, which ' ...
            'cannot be serialized to JSON; the saved model falls back to the ''keep'' policy on reload.'], nf.node));
    else
        nj('breakdownResetPolicy') = nf.breakdownResetPolicy;
    end
    if ~isempty(nf.repairResetPolicy)
        if strcmp(nf.repairResetPolicy, 'custom')
            line_warning(mfilename, sprintf(['Node failure on "%s" uses a custom repair reset function, which ' ...
                'cannot be serialized to JSON; the saved model falls back to the ''keep'' policy on reload.'], nf.node));
        else
            nj('repairResetPolicy') = nf.repairResetPolicy;
        end
    end
    nfJson{end+1} = nj; %#ok<AGROW>
end
if ~isempty(nfJson)
    result('nodeFailures') = nfJson;
end
end
 
 
% =========================================================================
%  Multiplicity serialization
% =========================================================================
 
function v = inf_multiplicity()
% Wire sentinel for infinite host/task multiplicity: Java's Integer.MAX_VALUE,
% which the JAR uses as its infinite-multiplicity marker. Written as a literal
% rather than taken from GlobalConstants.MaxInt, which is settable at runtime:
% a sentinel that varies per session would not survive a round-trip between two
% differently configured readers.
v = 2147483647;
end
 
% =========================================================================
%  Distribution serialization
% =========================================================================
 
function d = dist2json(dist)
% Convert a Distribution to a containers.Map for JSON output.
if isempty(dist)
    d = [];
    return;
end
d = containers.Map();
cn = builtin('class', dist);
switch cn
    case 'Disabled'
        d('type') = 'Disabled';
    case 'Immediate'
        d('type') = 'Immediate';
    case 'Exp'
        d('type') = 'Exp';
        params = containers.Map();
        params('lambda') = dist.getParam(1).paramValue;
        d('params') = params;
    case 'Det'
        d('type') = 'Det';
        params = containers.Map();
        params('value') = dist.getParam(1).paramValue;
        d('params') = params;
    case 'Erlang'
        d('type') = 'Erlang';
        params = containers.Map();
        params('lambda') = dist.getParam(1).paramValue;
        params('k') = dist.getParam(2).paramValue;
        d('params') = params;
    case 'HyperExp'
        % HyperExp.m stores the FULL rate vector in both param2 and param3 for
        % the n-phase form, so concatenating them emits a lambda of length 2n
        % that no reader can consume. Emit p and lambda both of length n. The
        % 2-phase form keeps its historical scalar-p storage (param2/param3 are
        % then the two distinct rates).
        d('type') = 'HyperExp';
        params = containers.Map();
        p = dist.getParam(1).paramValue;
        l1 = dist.getParam(2).paramValue;
        l2 = dist.getParam(3).paramValue;
        if isscalar(p)
            params('p') = [p, 1-p];
            params('lambda') = [l1, l2];
        else
            params('p') = p(:)';
            if isscalar(l1) && isscalar(l2)
                params('lambda') = [l1, l2];
            else
                % n-phase: param2 already holds all n rates (param3 duplicates it)
                params('lambda') = l1(:)';
            end
        end
        d('params') = params;
    case 'Gamma'
        d('type') = 'Gamma';
        params = containers.Map();
        params('alpha') = dist.getParam(1).paramValue;
        params('beta') = dist.getParam(2).paramValue;
        d('params') = params;
    case 'Lognormal'
        d('type') = 'Lognormal';
        params = containers.Map();
        params('mu') = dist.getParam(1).paramValue;
        params('sigma') = dist.getParam(2).paramValue;
        d('params') = params;
    case 'Uniform'
        d('type') = 'Uniform';
        params = containers.Map();
        params('a') = dist.getParam(1).paramValue;
        params('b') = dist.getParam(2).paramValue;
        d('params') = params;
    case 'Zipf'
        d('type') = 'Zipf';
        params = containers.Map();
        params('s') = dist.getParam(3).paramValue;
        params('n') = dist.getParam(4).paramValue;
        d('params') = params;
    case 'Pareto'
        d('type') = 'Pareto';
        params = containers.Map();
        params('alpha') = dist.getParam(1).paramValue;
        params('scale') = dist.getParam(2).paramValue;
        d('params') = params;
    case 'Weibull'
        d('type') = 'Weibull';
        params = containers.Map();
        params('alpha') = dist.getParam(1).paramValue;
        params('beta') = dist.getParam(2).paramValue;
        d('params') = params;
    case 'Normal'
        d('type') = 'Normal';
        params = containers.Map();
        params('mu') = dist.getParam(1).paramValue;
        params('sigma') = dist.getParam(2).paramValue;
        d('params') = params;
    case 'Geometric'
        d('type') = 'Geometric';
        params = containers.Map();
        params('p') = dist.getParam(1).paramValue;
        d('params') = params;
    case 'Binomial'
        d('type') = 'Binomial';
        params = containers.Map();
        params('n') = dist.getParam(1).paramValue;
        params('p') = dist.getParam(2).paramValue;
        d('params') = params;
    case 'Poisson'
        d('type') = 'Poisson';
        params = containers.Map();
        params('lambda') = dist.getParam(1).paramValue;
        d('params') = params;
    case 'Bernoulli'
        d('type') = 'Bernoulli';
        params = containers.Map();
        params('p') = dist.getParam(1).paramValue;
        d('params') = params;
    case 'DiscreteUniform'
        d('type') = 'DiscreteUniform';
        params = containers.Map();
        params('min') = dist.getParam(1).paramValue;
        params('max') = dist.getParam(2).paramValue;
        d('params') = params;
    case {'Coxian', 'Cox2'}
        d('type') = 'Coxian';
        params = containers.Map();
        params('mu') = dist.getMu()';
        params('phi') = dist.getPhi()';
        d('params') = params;
    case {'PH', 'APH'}
        % Keep the concrete class: an APH written back as a generic PH is a
        % lossy downgrade, since solver feature sets admit APH but not PH.
        d('type') = cn;
        ph = containers.Map();
        alpha = dist.getInitProb();
        T = dist.getSubgenerator();
        if isvector(alpha)
            ph('alpha') = alpha(:)';
        else
            ph('alpha') = alpha;
        end
        ph('T') = T;
        d('ph') = ph;
    case 'MAP'
        d('type') = 'MAP';
        mapSpec = containers.Map();
        mapSpec('D0') = dist.getParam(1).paramValue;
        mapSpec('D1') = dist.getParam(2).paramValue;
        d('map') = mapSpec;
    case 'DMAP'
        % Discrete-time MAP. Distinct from MAP on the wire: D0+D1 is stochastic,
        % not an infinitesimal generator, so a reader must not rebuild it as MAP.
        d('type') = 'DMAP';
        params = containers.Map();
        params('D0') = dist.getParam(1).paramValue;
        params('D1') = dist.getParam(2).paramValue;
        d('params') = params;
    case {'ME', 'CME'}
        % A CME goes on the wire as its (alpha, A) ME representation: the pair
        % determines the distribution completely, and every reader that accepts
        % ME accepts it. The subclass tag is not preserved by the round trip.
        d('type') = 'ME';
        params = containers.Map();
        alphaME = dist.getParam(1).paramValue;
        params('alpha') = alphaME(:)';
        params('A') = dist.getParam(2).paramValue;
        d('params') = params;
    case 'RAP'
        d('type') = 'RAP';
        params = containers.Map();
        params('H0') = dist.getParam(1).paramValue;
        params('H1') = dist.getParam(2).paramValue;
        d('params') = params;
    case 'BMAP'
        % Batch MAP: D = {D0, D1, ..., Dk}, Dk driving batches of size k. BMAP
        % subclasses MarkedMAP, so it must never be emitted through the MMAP
        % branch: the mark index is a batch size, not a class binding, and the
        % MarkedMAP ctor form would reinterpret the D1k as per-class arrivals.
        d('type') = 'BMAP';
        params = containers.Map();
        Kb = dist.getNumberOfTypes;
        dArr = cell(1, Kb + 1);
        dArr{1} = dist.getParam(1).paramValue;      % D0
        for kb = 1:Kb
            dArr{1+kb} = dist.getParam(2+kb).paramValue;   % Dk, batch size k
        end
        params('D') = dArr;
        d('params') = params;
    case 'MMDP2'
        d('type') = 'MMDP2';
        params = containers.Map();
        params('r0') = dist.getParam(1).paramValue;
        params('r1') = dist.getParam(2).paramValue;
        params('sigma0') = dist.getParam(3).paramValue;
        params('sigma1') = dist.getParam(4).paramValue;
        d('params') = params;
    case 'MarkedMMPP'
        % M3PP: D = {D0, D11, ..., D1K}; the aggregate D1 is rebuilt by the ctor
        % from the K == length(D)-1 form. MarkedMMPP extends MarkovModulated, not
        % MarkedMAP, so it never matched the MMAP branch and was lost entirely.
        d('type') = 'MarkedMMPP';
        params = containers.Map();
        Km = dist.getNumberOfTypes;
        dArr = cell(1, Km + 1);
        dArr{1} = dist.getParam(1).paramValue;      % D0
        for km = 1:Km
            dArr{1+km} = dist.getParam(2+km).paramValue;   % D1k
        end
        params('D') = dArr;
        params('K') = Km;
        d('params') = params;
    case 'EmpiricalCDF'
        % data is [F, x] rows (cdf value, support point), as assembled by the
        % two-argument ctor; emit the two columns separately.
        d('type') = 'EmpiricalCDF';
        params = containers.Map();
        ecdf = dist.data;
        params('F') = ecdf(:, 1)';
        params('x') = ecdf(:, 2)';
        d('params') = params;
    case 'Expolynomial'
        % Nested object, mirroring the Python writer (the density is a Sirio
        % expression string, not a numeric parameter).
        d('type') = 'Expolynomial';
        ep = containers.Map();
        ep('density') = dist.getParam(1).paramValue;
        ep('eft') = dist.getParam(2).paramValue;
        lft = dist.getParam(3).paramValue;
        if isfinite(lft)
            ep('lft') = lft;
        else
            ep('lft') = 'Inf';
        end
        d('expolynomial') = ep;
    case 'MarkedMAP'
        % Marked MAP: {D0, per-mark D1k}; the aggregate D1 is rebuilt on load
        d('type') = 'MMAP';
        mmapSpec = containers.Map();
        mmapSpec('D0') = dist.getParam(1).paramValue;
        Kmarks = dist.getNumberOfTypes;
        d1k = cell(1, Kmarks);
        for km = 1:Kmarks
            d1k{km} = dist.getParam(2+km).paramValue;
        end
        mmapSpec('D1k') = d1k;
        d('mmap') = mmapSpec;
    case 'MMPP2'
        d('type') = 'MMPP2';
        params = containers.Map();
        params('lambda0') = dist.getParam(1).paramValue;
        params('lambda1') = dist.getParam(2).paramValue;
        params('sigma0') = dist.getParam(3).paramValue;
        params('sigma1') = dist.getParam(4).paramValue;
        d('params') = params;
    case 'NHPP'
        d('type') = 'NHPP';
        params = containers.Map();
        % num2cell keeps a single-segment (1x1) rate vector from collapsing to
        % a JSON scalar, which the Gson reader would reject.
        nhppBp = double(dist.getBreakpoints());
        nhppRt = double(dist.getRates());
        params('breakpoints') = num2cell(nhppBp(:)');
        params('rates') = num2cell(nhppRt(:)');
        params('cyclic') = dist.isCyclic();
        d('params') = params;
    case 'DiscreteSampler'
        d('type') = 'DiscreteSampler';
        params = containers.Map();
        params('p') = dist.getParam(1).paramValue(:)';
        params('x') = dist.getParam(2).paramValue(:)';
        d('params') = params;
    case 'Replayer'
        d('type') = 'Replayer';
        params = containers.Map();
        params('fileName') = dist.getParam(1).paramValue;
        try
            params('mean') = dist.getMean();
        catch
        end
        d('params') = params;
        % Save APH fit as fallback
        try
            aphDist = dist.fitAPH();
            if ~isempty(aphDist) && isa(aphDist, 'Distribution')
                ph = containers.Map();
                alpha = aphDist.getParam(1).paramValue;
                T = aphDist.getParam(2).paramValue;
                if isvector(alpha)
                    ph('alpha') = alpha(:)';
                else
                    ph('alpha') = alpha;
                end
                ph('T') = T;
                d('ph') = ph;
            end
        catch
        end
    case 'Prior'
        d('type') = 'Prior';
        alts = {};
        for ai = 1:dist.getNumAlternatives()
            altDist = dist.getAlternative(ai);
            altJson = dist2json(altDist);
            if ~isempty(altJson)
                alts{end+1} = altJson; %#ok<AGROW>
            end
        end
        d('distributions') = alts;
        d('probabilities') = dist.probabilities(:)';
    otherwise
        % No branch matches. Warn and emit the real type name with the first two
        % moments, which the readers reconstruct via APH.fitMeanAndSCV. Emitting
        % 'Exp' + fitMean was a silent degradation on two counts: it discarded the
        % SCV, and it lied about the type, so a reader could not even tell that
        % information had been lost.
        line_warning(mfilename, sprintf(['Distribution "%s" has no JSON representation; ' ...
            'saving its mean and SCV only. The reloaded model will use an APH fitted ' ...
            'to those two moments.\n'], cn));
        m = dist.getMean();
        s = dist.getSCV();
        d('type') = cn;
        params = containers.Map();
        params('mean') = m;
        params('scv') = s;
        d('params') = params;
end
end
 
 
% =========================================================================
%  JSON encoding
% =========================================================================
 
function s = encode_value(val, indent)
% Recursively encode a MATLAB value to JSON string.
if nargin < 2, indent = 0; end
pad  = repmat(' ', 1, indent);
pad2 = repmat(' ', 1, indent + 2);
 
if isa(val, 'containers.Map')
    ks = val.keys();
    if isempty(ks)
        s = '{}';
    else
        parts = cell(1, length(ks));
        for i = 1:length(ks)
            k = ks{i};
            v = val(k);
            parts{i} = sprintf('%s"%s": %s', pad2, json_escape(k), encode_value(v, indent + 2));
        end
        s = sprintf('{\n%s\n%s}', strjoin(parts, sprintf(',\n')), pad);
    end
elseif ischar(val) || isstring(val)
    s = sprintf('"%s"', json_escape(char(val)));
elseif islogical(val) && isscalar(val)
    if val, s = 'true'; else, s = 'false'; end
elseif isnumeric(val) && isscalar(val)
    if isnan(val)
        s = 'null';
    elseif isinf(val)
        if val > 0, s = '"Infinity"'; else, s = '"-Infinity"'; end
    elseif val == floor(val) && abs(val) < 1e15
        s = sprintf('%d', val);
    else
        s = sprintf('%.15g', val);
    end
elseif isnumeric(val) && isvector(val) && ~isscalar(val)
    parts = cell(1, length(val));
    for i = 1:length(val)
        parts{i} = encode_value(val(i), 0);
    end
    s = ['[', strjoin(parts, ', '), ']'];
elseif isnumeric(val) && ismatrix(val) && ~isvector(val)
    rows = cell(1, size(val, 1));
    for i = 1:size(val, 1)
        rows{i} = encode_value(val(i,:), 0);
    end
    s = ['[', strjoin(rows, ', '), ']'];
elseif iscell(val)
    if isempty(val)
        s = '[]';
    else
        parts = cell(1, length(val));
        for i = 1:length(val)
            parts{i} = sprintf('%s%s', pad2, encode_value(val{i}, indent + 2));
        end
        s = sprintf('[\n%s\n%s]', strjoin(parts, sprintf(',\n')), pad);
    end
elseif isstruct(val) && isscalar(val)
    fnames = fieldnames(val);
    if isempty(fnames)
        s = '{}';
    else
        parts = cell(1, length(fnames));
        for i = 1:length(fnames)
            fn = fnames{i};
            fv = val.(fn);
            parts{i} = sprintf('%s"%s": %s', pad2, json_escape(fn), encode_value(fv, indent + 2));
        end
        s = sprintf('{\n%s\n%s}', strjoin(parts, sprintf(',\n')), pad);
    end
else
    s = 'null';
end
end
 
function s = json_escape(str)
% Escape special characters for JSON strings.
s = strrep(str, '\', '\\');
s = strrep(s, '"', '\"');
s = strrep(s, sprintf('\n'), '\n');
s = strrep(s, sprintf('\r'), '\r');
s = strrep(s, sprintf('\t'), '\t');
end
 
 
% =========================================================================
%  Helper functions
% =========================================================================
 
function s = node_type_str(node)
% Get the JSON node type string for a node object.
if isa(node, 'Source'),      s = 'Source';
elseif isa(node, 'Sink'),   s = 'Sink';
elseif isa(node, 'Delay'),  s = 'Delay';
elseif isa(node, 'Cache'),  s = 'Cache';
elseif isa(node, 'Place'),  s = 'Place';
elseif isa(node, 'Transition'), s = 'Transition';
elseif isa(node, 'Queue'),  s = 'Queue';
elseif isa(node, 'Fork'),   s = 'Fork';
elseif isa(node, 'Join'),   s = 'Join';
elseif isa(node, 'Router'), s = 'Router';
elseif isa(node, 'ClassSwitch'), s = 'ClassSwitch';
else,                        s = 'Queue';
end
end
 
function s = sched_id_to_str(id)
% Map a SchedStrategy numeric ID to the wire enum name.
%
% The wire carries enum NAMES, uppercased, matching the JAR enum constants
% (jline.lang.constant.SchedStrategy) one-for-one for all 40 strategies. Do not
% reintroduce a hand-rolled whitelist here: the previous one covered 23 of 40 and
% silently degraded SRPT/FSP/EDD/EDF/FB/LCFSPI/PSJF/LRPT/SETF/LPS/... to FCFS.
% SchedStrategy.toText errors on an unknown id rather than inventing a default.
s = upper(SchedStrategy.toText(SchedStrategy.toId(id)));
end
 
function s = repl_to_str(id)
% Map ReplacementStrategy numeric ID to the wire enum name.
if id == ReplacementStrategy.LRU,   s = 'LRU';
elseif id == ReplacementStrategy.FIFO, s = 'FIFO';
elseif id == ReplacementStrategy.RR, s = 'RR';
elseif id == ReplacementStrategy.SFIFO, s = 'SFIFO';
elseif id == ReplacementStrategy.HLRU, s = 'HLRU';
elseif id == ReplacementStrategy.CLIMB, s = 'CLIMB';
elseif id == ReplacementStrategy.QLRU, s = 'QLRU';
else
    line_error(mfilename, sprintf('Unrecognized replacement strategy id %d.', id));
end
end
 
function s = depdisc_to_str(id)
% Map a DepartureDiscipline numeric ID to the wire name. The names are the JAR
% enum constants (jline.lang.constant.DepartureDiscipline), which is what the
% JAR writer emits and what its reader matches case-insensitively.
if id == DepartureDiscipline.NORMAL,    s = 'Normal';
elseif id == DepartureDiscipline.FIFO,  s = 'FIFO';
else
    line_error(mfilename, sprintf('Unrecognized departure discipline id %d.', id));
end
end
 
function s = droprule_to_str(id)
% Map DropStrategy numeric ID to schema-compatible string.
if id == DropStrategy.DROP,        s = 'drop';
elseif id == DropStrategy.WAITQ,   s = 'waitingQueue';
elseif id == DropStrategy.BAS,     s = 'blockingAfterService';
elseif id == DropStrategy.RETRIAL, s = 'retrial';
elseif id == DropStrategy.RETRIAL_WITH_LIMIT, s = 'retrialWithLimit';
else,                              s = '';
end
end
 
function rateMap = oi_rate_table(muFun, maxc, K)
% Build the OI/PAS macrostate rate table: for every per-class count vector cnt
% on the box lattice 0 <= cnt(r) <= maxc(r), evaluate mu on a canonical ordered
% microstate holding cnt(r) copies of class r (any ordering is valid since mu is
% order-independent). Keyed by the comma-joined 0-based class counts, matching
% the JAR reader (LineModelIO.oiServiceRate). The empty state is omitted.
rateMap = containers.Map('KeyType', 'char', 'ValueType', 'double');
shp = maxc + 1;
total = prod(shp);
for i = 1:total
    li = i - 1; cnt = zeros(1, K);
    for d = 1:K, cnt(d) = mod(li, shp(d)); li = floor(li / shp(d)); end
    if sum(cnt) == 0, continue, end
    micro = repelem(1:K, cnt);
    rate = muFun(micro);
    if ~isfinite(rate), rate = 0; end
    key = strjoin(arrayfun(@(x) sprintf('%d', x), cnt, 'UniformOutput', false), ',');
    rateMap(key) = rate;
end
end
 
function tbl = cd_scaling_table(beta, maxc, K)
% Materialize the class-dependence handle beta(n) over the box lattice
% 0 <= n(r) <= maxc(r). Keyed by the comma-joined 0-based per-class counts,
% matching the JAR reader (LineModelIO, "classDependence"). The handle may
% return a scalar (one scaling shared by every class) or a length-K vector of
% per-class scalings; the scalar form is broadcast to K entries here so that the
% reader is uniform and need not re-derive which form was used.
tbl = containers.Map('KeyType', 'char', 'ValueType', 'any');
shp = maxc + 1;
total = prod(shp);
for i = 1:total
    li = i - 1; n = zeros(1, K);
    for d = 1:K, n(d) = mod(li, shp(d)); li = floor(li / shp(d)); end
    v = beta(n);
    if isscalar(v)
        v = repmat(double(v), 1, K);
    else
        v = double(v(:)');
    end
    v(~isfinite(v)) = 0;
    key = strjoin(arrayfun(@(x) sprintf('%d', x), n, 'UniformOutput', false), ',');
    tbl(key) = num2cell(v);
end
end
 
function s = prectype_to_str(id)
% Map ActivityPrecedenceType numeric ID to JAR-compatible string.
if id == ActivityPrecedenceType.PRE_SEQ,        s = 'pre';
elseif id == ActivityPrecedenceType.PRE_AND,    s = 'pre-AND';
elseif id == ActivityPrecedenceType.PRE_OR,     s = 'pre-OR';
elseif id == ActivityPrecedenceType.POST_SEQ,   s = 'post';
elseif id == ActivityPrecedenceType.POST_AND,   s = 'post-AND';
elseif id == ActivityPrecedenceType.POST_OR,    s = 'post-OR';
elseif id == ActivityPrecedenceType.POST_LOOP,  s = 'post-LOOP';
elseif id == ActivityPrecedenceType.POST_CACHE, s = 'post-CACHE';
else,                                           s = 'pre';
end
end
 
function tf = isIdentityClassSwitch(csm, K)
% True if the classSwitchMatrix CSM is the K x K identity (or empty/unset,
% which defaults to identity). A non-identity matrix carries the class switch
% and must be collapsed to same-class routing on export; an identity matrix
% means the switch is expressed through the routing and must be kept verbatim.
tf = true;
if isempty(csm)
    return;
end
for rr = 1:min(K, size(csm,1))
    for ss = 1:min(K, size(csm,2))
        if abs(csm(rr,ss) - double(rr==ss)) > 1e-12
            tf = false;
            return;
        end
    end
end
end