ldes_cli.cpp

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
 
/**
 * @file
 * `ldes`: the native LDES engine behind the interface of `jline.cli.LdesCLI`.
 *
 * WHAT THIS BINARY IS FOR. `common/ldes` has been a GraalVM image of the JAVA
 * engine; this program is its replacement, built from
 * `cpp/include/line/solvers/ldes/`. Every client of that binary -- MATLAB's
 * `@SolverLDES/solveCli.m`, native Python's `wrappers/solver_ldes`, and the C++
 * client `wrappers/ldes/solver_ldes.h` -- speaks to it through exactly two
 * artefacts: the ARGUMENT VECTOR and the `ldes-result` DOCUMENT. Neither may
 * drift. A client is not recompiled when the engine changes hands, so an
 * argument this program rejects and the Java one accepted, or a key it spells
 * differently, is a silent breakage in three codebases at once.
 *
 * THE INTERFACE IS THEREFORE COPIED, not designed:
 *
 *   ldes solve <model.json> -o <result.json> [flags]
 *
 * with the flag set of `LdesCLI.handleSolveCommand` and the document shape of
 * `LDESResultIO.write`. Flags this engine cannot honour are REFUSED BY NAME
 * rather than ignored, which is the one deliberate difference from the AOT
 * image it replaces: that image ignored flags it predated, and an ignored
 * `--slotlength` silently ran a continuous-time simulation under a slotted
 * model's name. A refusal is visible; a silent default is not.
 *
 * A LAYEREDNETWORK document (`modelType: LayeredNetwork`) is dispatched to the
 * native LN engine, `ldes_ln_engine_solve`, and answered with the layered
 * document of `LDESResultIO.saveLN` (see `ln_result_json`). Flags that engine
 * has no counterpart for (transient, histograms, busy periods, warm starts,
 * slotted time, replications) are refused by name.
 *
 * WHAT IS NOT HERE. `--rest`/server mode belongs to `line_cli` and is refused
 * with the reason rather than accepted and dropped.
 */
 
#include <chrono>
#include <cmath>
#include <cstdio>
#include <ctime>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>
 
#include "line/io/lqn_json_reader.h"
#include "line/io/network_reader.h"
#include "line/solvers/ldes/ldes_engine.h"
#include "line/solvers/ldes/ldes_ln_engine.h"
#include "line/solvers/wrappers/ldes/ldes_options.h"
#include "line/util/error.h"
 
namespace {
 
using line::ldes::LdesOptions;
using line::ldes::LdesResult;
 
void print_help() {
    std::cout
        << "LDES: discrete-event simulation engine for LINE models.\n\n"
        << "Usage: ldes solve <model.json> [-o <result.json>] [options]\n\n"
        << "Options:\n"
        << "  -s, --samples N          service completions to simulate (default 200000)\n"
        << "  -e, --maxevents N        event budget; overrides --samples\n"
        << "      --maxtime SECONDS    wall-clock budget\n"
        << "      --seed N             run seed (default 23000; -1 draws one)\n"
        << "      --method NAME        engine method (default 'default')\n"
        << "      --cnvgon             stop on the convergence check\n"
        << "      --cnvgtol X          relative tolerance of that check\n"
        << "      --cnvgbatch N        batches before the first check\n"
        << "      --cnvgchk N          events between checks\n"
        << "      --tranfilter NAME    mser5 | fixed | none\n"
        << "      --mserbatch N        MSER batch size\n"
        << "      --warmupfrac X       warmup fraction, for --tranfilter fixed\n"
        << "      --tranobs N          transient observation points over --timespan\n"
        << "      --cimethod NAME      obm | bm | spectral | none\n"
        << "      --obmoverlap X       overlap of the batch means\n"
        << "      --ciminbatch N       batches below which no CI is reported\n"
        << "      --ciminobs N         observations below which no CI is reported\n"
        << "      --spectrallowfreqfrac X   low-frequency fraction of the spectral CI\n"
        << "      --slotted            run on the slot lattice\n"
        << "      --slotlength X       slot length (implies --slotted)\n"
        << "      --replications N     independent replications\n"
        << "      --numthreads N       threads across replications\n"
        << "      --timespan T0,T1     transient run over [T0,T1]\n"
        << "      --busyperiod K       busy-period orders 1..K\n"
        << "      --busyperiod-subnet i,j,...   an extra busy-period target\n"
        << "      --initsol V          initial state, station-major, comma separated\n"
        << "      --initjobs J         job-level initial state, one row per job, rows ';'\n"
        << "                           separated: station,class,remaining,elapsed,qage,sage,order[,id]\n"
        << "      --trackjobs          export the final job-level state and tagged departures\n"
        << "      --export-histogram   fill the state histogram block\n"
        << "      --trajectory         fill the state trajectory block\n"
        << "      --respt-samples      record every per-visit response time\n"
        << "  -h, --help               this text\n";
}
 
[[noreturn]] void fail(const std::string& msg) {
    std::cerr << "Error: " << msg << "\n";
    std::exit(1);
}
 
/** A flag the reference accepts and this engine cannot honour: say so. */
[[noreturn]] void refuse(const std::string& flag, const std::string& why) {
    std::cerr << "Error: " << flag << " is accepted by the Java engine but not by this one: "
              << why << ".\nRun the model through common/ldes.jar for it, rather than "
                        "taking a silently different simulation.\n";
    std::exit(1);
}
 
std::string need_value(int argc, char** argv, int& i, const std::string& flag) {
    if (i + 1 >= argc) fail(flag + " requires a value.");
    return std::string(argv[++i]);
}
 
std::vector<double> parse_doubles(const std::string& csv) {
    std::vector<double> out;
    std::stringstream ss(csv);
    std::string tok;
    while (std::getline(ss, tok, ',')) {
        if (!tok.empty()) out.push_back(std::atof(tok.c_str()));
    }
    return out;
}
 
std::vector<std::size_t> parse_indices(const std::string& csv) {
    std::vector<std::size_t> out;
    std::stringstream ss(csv);
    std::string tok;
    while (std::getline(ss, tok, ',')) {
        if (!tok.empty()) out.push_back(static_cast<std::size_t>(std::atol(tok.c_str())));
    }
    return out;
}
 
/** A JSON number the way the reference writes it: full round-trip precision. */
std::string num(double v) {
    if (std::isnan(v)) return "null";
    if (std::isinf(v)) return v > 0 ? "1e999" : "-1e999";
    char buf[40];
    std::snprintf(buf, sizeof(buf), "%.17g", v);
    return std::string(buf);
}
 
/** A (rows x cols) matrix as the nested array `LDESResultIO` emits. */
std::string mat_json(const line::Matrix<double>& m, std::size_t rows, std::size_t cols) {
    std::string s = "[";
    for (std::size_t i = 0; i < rows; ++i) {
        if (i) s += ", ";
        s += "[";
        for (std::size_t j = 0; j < cols; ++j) {
            if (j) s += ", ";
            const bool inside = (i < m.rows() && j < m.cols());
            s += num(inside ? m(i, j) : 0.0);
        }
        s += "]";
    }
    s += "]";
    return s;
}
 
/** The per-visit response times as the nested `[station][class][sample]` array. */
std::string respt_samples_json(const LdesResult& r) {
    if (r.respTimeSamples.empty()) return std::string();
    std::string s = "[";
    for (std::size_t i = 0; i < r.respTimeSamples.size(); ++i) {
        if (i) s += ", ";
        s += "[";
        for (std::size_t k = 0; k < r.respTimeSamples[i].size(); ++k) {
            if (k) s += ", ";
            s += "[";
            for (std::size_t j = 0; j < r.respTimeSamples[i][k].size(); ++j) {
                if (j) s += ", ";
                s += num(r.respTimeSamples[i][k][j]);
            }
            s += "]";
        }
        s += "]";
    }
    s += "]";
    return s;
}
 
std::string str_json(const std::string& s) {
    std::string out = "\"";
    for (char c : s) {
        if (c == '"' || c == '\\') {
            out += '\\';
            out += c;
        } else if (c == '\n') {
            out += "\\n";
        } else {
            out += c;
        }
    }
    out += "\"";
    return out;
}
 
/**
 * The `ldes-result` document, key for key with `LDESResultIO.write`.
 *
 * The order of the keys is the reference's, and the blocks a run did not fill
 * are still emitted with their zero matrices: every client indexes this
 * document by key and by station row, and a missing block reads as a parse
 * failure rather than as an absent measurement.
 */
std::string result_json(const LdesResult& r, const LdesOptions& o, double runtime,
                        std::size_t events) {
    const std::size_t S = r.nstations, C = r.nclasses;
    std::string s;
    s += "{\n";
    s += "  \"format\": \"ldes-result\",\n";
    s += "  \"version\": \"1.0\",\n";
    s += "  \"solver\": \"SolverLDES\",\n";
    s += "  \"method\": " + str_json(o.method) + ",\n";
    s += "  \"runtime\": " + num(runtime) + ",\n";
    s += "  \"converged\": false,\n";
    s += "  \"stoppingReason\": \"max_events\",\n";
    s += "  \"convergenceBatches\": 0,\n";
    s += "  \"totalSimulatedEvents\": " + std::to_string(events) + ",\n";
 
    s += "  \"dimensions\": {\"nstations\": " + std::to_string(S) +
         ", \"nclasses\": " + std::to_string(C) + ", \"nchains\": " + std::to_string(r.nchains) +
         ", \"stationNames\": [";
    for (std::size_t i = 0; i < r.station_names.size(); ++i) {
        if (i) s += ", ";
        s += str_json(r.station_names[i]);
    }
    s += "], \"classNames\": [";
    for (std::size_t i = 0; i < r.class_names.size(); ++i) {
        if (i) s += ", ";
        s += str_json(r.class_names[i]);
    }
    s += "]},\n";
 
    s += "  \"metrics\": {";
    s += "\"QN\": " + mat_json(r.QN, S, C);
    s += ", \"UN\": " + mat_json(r.UN, S, C);
    s += ", \"RN\": " + mat_json(r.RN, S, C);
    s += ", \"TN\": " + mat_json(r.TN, S, C);
    s += ", \"AN\": " + mat_json(r.AN, S, C);
    s += ", \"WN\": " + mat_json(r.WN, S, C);
    s += ", \"CN\": " + mat_json(r.CN, 1, C);
    s += ", \"XN\": " + mat_json(r.XN, 1, C);
    // The sibling-drop rate at a Join, where LossRate = ArvR - Tput does NOT
    // hold: ArvR counts the SIBLINGS offered and Tput the PARENT jobs released.
    // `LDESResultIO.write` puts it inside `metrics`, so it goes here and not in
    // an optional block of its own.
    if (!r.DropRateJoin.empty()) s += ", \"DropRateJoin\": " + mat_json(r.DropRateJoin, S, C);
    s += "},\n";
 
    // THE OPTIONAL BLOCKS, each written ONLY when the run measured it, exactly
    // as `LDESResultIO.write` gates them. A block emitted empty is not the same
    // as an absent one: `ldes_prob_from_histogram` and `getAvgBusyPeriod` both
    // treat presence as the claim that the measurement was taken.
    // GATED ON THE FLAG, not on the measurement. A transient run computes the
    // histogram whether or not anyone asked for it, and `LDESResultIO` writes
    // the block only when the flag was given; emitting it anyway would make
    // this engine's document differ from the reference's on every transient run.
    if (o.export_histogram && r.histogram_space.rows() > 0 && r.histogram_time.rows() > 0) {
        s += "  \"stateHistogram\": {";
        s += "\"space\": " + mat_json(r.histogram_space, r.histogram_space.rows(),
                                      r.histogram_space.cols());
        s += ", \"time\": " + mat_json(r.histogram_time, r.histogram_time.rows(), 1);
        if (r.traj_space.rows() > 0 && r.traj_time.rows() > 0) {
            s += ", \"trajSpace\": " +
                 mat_json(r.traj_space, r.traj_space.rows(), r.traj_space.cols());
            s += ", \"trajTime\": " + mat_json(r.traj_time, r.traj_time.rows(), 1);
        }
        s += "},\n";
    }
 
    if (!r.busy_periods.empty()) {
        s += "  \"busyPeriods\": {\"orders\": " +
             std::to_string(r.busy_periods[0].mean.size()) + ", \"targets\": [";
        for (std::size_t ti = 0; ti < r.busy_periods.size(); ++ti) {
            const LdesResult::BusyPeriodTarget& t = r.busy_periods[ti];
            if (ti) s += ", ";
            s += "{\"name\": " + str_json(t.name) + ", \"stations\": [";
            for (std::size_t k = 0; k < t.stations.size(); ++k) {
                if (k) s += ", ";
                s += std::to_string(t.stations[k]);
            }
            s += "], \"class\": " + std::to_string(t.job_class) + ", \"mean\": [";
            for (std::size_t k = 0; k < t.mean.size(); ++k) {
                if (k) s += ", ";
                s += num(t.mean[k]);
            }
            s += "], \"count\": [";
            for (std::size_t k = 0; k < t.count.size(); ++k) {
                if (k) s += ", ";
                s += num(t.count[k]);
            }
            s += "]}";
        }
        s += "]},\n";
    }
 
    // The job-level end state, in the --initjobs layout, so a run can be continued from it.
    if (o.track_jobs || o.init_jobs.rows() > 0) {
        s += "  \"jobState\": {\"finalJobs\": " +
             mat_json(r.finalJobs, r.finalJobs.rows(), r.finalJobs.cols());
        s += ", \"jobDepartures\": " +
             mat_json(r.jobDepartures, r.jobDepartures.rows(), r.jobDepartures.cols()) + "},\n";
    }
 
    if (!r.cache_metrics.empty()) {
        s += "  \"cacheMetrics\": {";
        bool first = true;
        for (std::map<std::string, line::ldes::LdesCacheMetrics>::const_iterator it =
                 r.cache_metrics.begin();
             it != r.cache_metrics.end(); ++it) {
            if (!first) s += ", ";
            first = false;
            const line::ldes::LdesCacheMetrics& cm = it->second;
            s += str_json(it->first) + ": {";
            s += "\"hit\": " + mat_json(cm.hit, cm.hit.rows(), cm.hit.cols());
            s += ", \"delayed\": " + mat_json(cm.delayed, cm.delayed.rows(), cm.delayed.cols());
            s += ", \"miss\": " + mat_json(cm.miss, cm.miss.rows(), cm.miss.cols());
            s += ", \"latency\": " + mat_json(cm.latency, cm.latency.rows(), cm.latency.cols());
            s += ", \"hitList\": " + mat_json(cm.hitList, cm.hitList.rows(), cm.hitList.cols());
            s += ", \"itemProb\": " + mat_json(cm.itemProb, cm.itemProb.rows(), cm.itemProb.cols());
            s += ", \"listCost\": " + mat_json(cm.listCost, cm.listCost.rows(), cm.listCost.cols());
            s += "}";
        }
        s += "},\n";
    }
 
    if (r.QNfcr.rows() > 0) {
        const std::size_t G = r.QNfcr.rows();
        s += "  \"fcr\": {\"nregions\": " + std::to_string(G);
        s += ", \"QNfcr\": " + mat_json(r.QNfcr, G, C);
        s += ", \"UNfcr\": " + mat_json(r.UNfcr, G, C);
        s += ", \"RNfcr\": " + mat_json(r.RNfcr, G, C);
        s += ", \"TNfcr\": " + mat_json(r.TNfcr, G, C);
        s += ", \"ANfcr\": " + mat_json(r.ANfcr, G, C);
        s += ", \"WNfcr\": " + mat_json(r.WNfcr, G, C);
        s += ", \"WeightNfcr\": " + mat_json(r.WeightNfcr, G, C);
        s += ", \"MemOccNfcr\": " + mat_json(r.MemOccNfcr, G, C);
        s += ", \"DropRateNfcr\": " + mat_json(r.DropRateNfcr, G, C);
        s += "},\n";
    }
 
    // Fork-join section metrics (getAvgFJ), TOP-LEVEL as the Java engine writes them: one row per Fork/Join pair,
    // named by "forks"/"joins", in fork-then-join node order. Steady state only, and only on a model with a pair.
    if (r.t.empty() && !r.fj_forks.empty() && r.FJTput.rows() == r.fj_forks.size()) {
        const std::size_t P = r.fj_forks.size();
        s += "  \"forkJoin\": {\"forks\": [";
        for (std::size_t p = 0; p < P; ++p) s += (p ? ", " : "") + str_json(r.fj_forks[p]);
        s += "], \"joins\": [";
        for (std::size_t p = 0; p < P; ++p) s += (p ? ", " : "") + str_json(r.fj_joins[p]);
        s += "]";
        s += ", \"FJQLen\": " + mat_json(r.FJQLen, P, C);
        s += ", \"FJRespT\": " + mat_json(r.FJRespT, P, C);
        s += ", \"FJTput\": " + mat_json(r.FJTput, P, C);
        s += "},\n";
    }
 
    // The reference gates the whole transient block on `--trajectory`, not on
    // the run being transient.
    // The samples appear at TOP LEVEL under --respt-samples and INSIDE the
    // transient block under --trajectory. Both spellings are the reference's,
    // and both are read: the C++ client takes whichever is present, MATLAB's
    // `sample()` and `sampleSys()` read the nested one and returned nothing
    // while only the top-level spelling existed.
    const std::string respt_json = respt_samples_json(r);
    if (o.export_trajectory && !r.t.empty()) {
        s += "  \"transient\": {\"t\": [";
        for (std::size_t k = 0; k < r.t.size(); ++k) {
            if (k) s += ", ";
            s += num(r.t[k]);
        }
        s += "]";
        const char* keys[3] = {"QNt", "UNt", "TNt"};
        const std::vector<std::vector<line::Matrix<double>>>* src[3] = {&r.QNt, &r.UNt, &r.TNt};
        for (int q = 0; q < 3; ++q) {
            s += std::string(", \"") + keys[q] + "\": [";
            for (std::size_t i = 0; i < src[q]->size(); ++i) {
                if (i) s += ", ";
                s += "[";
                for (std::size_t k = 0; k < (*src[q])[i].size(); ++k) {
                    if (k) s += ", ";
                    const line::Matrix<double>& mm = (*src[q])[i][k];
                    s += mat_json(mm, mm.rows(), mm.cols());
                }
                s += "]";
            }
            s += "]";
        }
        if (!respt_json.empty()) s += ", \"respTimeSamples\": " + respt_json;
        s += "},\n";
    }
 
    if (o.export_respt && !respt_json.empty())
        s += "  \"respTimeSamples\": " + respt_json + ",\n";
    s += "  \"sampleCounts\": {";
    s += "\"QNSamples\": " + mat_json(r.QNSamples, S, C);
    s += ", \"UNSamples\": " + mat_json(r.UNSamples, S, C);
    s += ", \"RNSamples\": " + mat_json(r.RNSamples, S, C);
    s += ", \"TNSamples\": " + mat_json(r.TNSamples, S, C);
    s += "},\n";
 
    s += "  \"confidenceIntervals\": {";
    s += "\"QNCI\": " + mat_json(r.QNCI, S, C);
    s += ", \"UNCI\": " + mat_json(r.UNCI, S, C);
    s += ", \"RNCI\": " + mat_json(r.RNCI, S, C);
    s += ", \"TNCI\": " + mat_json(r.TNCI, S, C);
    s += ", \"ANCI\": " + mat_json(r.ANCI, S, C);
    s += ", \"WNCI\": " + mat_json(r.WNCI, S, C);
    s += "},\n";
 
    s += "  \"relativePrecision\": {";
    s += "\"QNRelPrec\": " + mat_json(r.QNRelPrec, S, C);
    s += ", \"UNRelPrec\": " + mat_json(r.UNRelPrec, S, C);
    s += ", \"RNRelPrec\": " + mat_json(r.RNRelPrec, S, C);
    s += ", \"TNRelPrec\": " + mat_json(r.TNRelPrec, S, C);
    s += "},\n";
 
    s += "  \"impatience\": {";
    s += "\"renegedCustomers\": " + mat_json(r.renegedCustomers, S, C);
    s += ", \"avgRenegingWaitTime\": " + mat_json(r.avgRenegingWaitTime, S, C);
    s += ", \"renegingRate\": " + mat_json(r.renegingRate, S, C);
    s += ", \"balkedCustomers\": " + mat_json(r.balkedCustomers, S, C);
    s += ", \"balkingProbability\": " + mat_json(r.balkingProbability, S, C);
    s += ", \"retriedCustomers\": " + mat_json(r.retriedCustomers, S, C);
    s += ", \"retrialDropped\": " + mat_json(r.retrialDropped, S, C);
    s += ", \"avgOrbitSize\": " + mat_json(r.avgOrbitSize, S, C);
    s += "}\n";
    s += "}\n";
    return s;
}
 
/** A 1-based (nidx+1) x 1 element column as the 1 x nidx row `LDESResultIO.saveLN` emits. */
std::string ln_row_json(const line::Matrix<double>& m, std::size_t nidx) {
    std::string s = "[[";
    for (std::size_t k = 1; k <= nidx; ++k) {
        if (k > 1) s += ", ";
        s += num(k < m.rows() ? m(k, 0) : 0.0);
    }
    return s + "]]";
}
 
/**
 * The LAYERED `ldes-result` document, key for key with `LDESResultIO.saveLN`, in the JAR's
 * 0-based element index space (C++ element k is entry k-1 of every row and of `cacheMetrics`).
 *
 * UCallLN, TCallLN and UEntryClassLN ARE NOT EMITTED: the native LN engine does not measure the
 * per-call and per-entry-stream thread occupancy, and no subprocess client reads them (they back
 * the JAR-only getLNAvgCallUtil/CallTput/EntryClassUtil). The CI block is all null, as the
 * JAR's layered engine leaves it.
 */
std::string ln_result_json(const line::lqn::LqnStruct<double>& lsn,
                           const line::ldes::engine::LnResult& r, const LdesOptions& o,
                           double runtime, bool respt) {
    const std::size_t n = lsn.nidx;
    // THE DOCUMENT CARRIES THE JAR ENGINE'S RAW CONVENTIONS, which differ from the LnResult table:
    // ULN is MEAN BUSY SERVERS, undivided, because every client divides it by hostReachMult itself
    // (Python getLNAvgTable); emitting the per-server figure would divide twice. QLN at a processor
    // and RLN at a task are 0, as Solver_ssj_ln writes them; the native engine measures both, and
    // LnResult keeps them for an in-process caller.
    line::Matrix<double> QLN = r.QLN, ULN = r.ULN, RLN = r.RLN;
    for (std::size_t k = 1; k <= n; ++k) {
        ULN(k, 0) *= lsn.host_reach_mult(k);
        if (k > lsn.hshift && k <= lsn.hshift + lsn.nhosts) QLN(k, 0) = 0.0;
        if (k > lsn.tshift && k <= lsn.tshift + lsn.ntasks) RLN(k, 0) = 0.0;
    }
    std::string s = "{\"format\": \"ldes-result\", \"version\": \"1.0\", "
                    "\"modelType\": \"LayeredNetwork\", \"solver\": \"LDES\", ";
    s += "\"method\": " + str_json(o.method) + ", \"runtime\": " + num(runtime) + ", ";
    s += "\"dimensions\": {\"nidx\": " + std::to_string(n) +
         ", \"nhosts\": " + std::to_string(lsn.nhosts) +
         ", \"ntasks\": " + std::to_string(lsn.ntasks) +
         ", \"nentries\": " + std::to_string(lsn.nentries) +
         ", \"nacts\": " + std::to_string(lsn.nacts) +
         ", \"ncalls\": " + std::to_string(lsn.ncalls) +
         ", \"tshift\": " + std::to_string(lsn.tshift) +
         ", \"eshift\": " + std::to_string(lsn.eshift) +
         ", \"ashift\": " + std::to_string(lsn.ashift) + ", \"names\": [";
    for (std::size_t k = 1; k <= n; ++k) {
        if (k > 1) s += ", ";
        s += str_json(k < lsn.names.size() ? lsn.names[k] : "Element" + std::to_string(k - 1));
    }
    s += "]}, \"metrics\": {";
    s += "\"QLN\": " + ln_row_json(QLN, n) + ", \"ULN\": " + ln_row_json(ULN, n);
    s += ", \"RLN\": " + ln_row_json(RLN, n) + ", \"WLN\": " + ln_row_json(r.WLN, n);
    s += ", \"TLN\": " + ln_row_json(r.TLN, n) + ", \"ALN\": " + ln_row_json(r.ALN, n);
    s += ", \"ZLN\": " + ln_row_json(r.ZLN, n) + "}, ";
    s += "\"confidenceIntervals\": {\"QLNCI\": null, \"ULNCI\": null, \"RLNCI\": null, "
         "\"TLNCI\": null}";
    if (!r.cache_task.empty()) {
        s += ", \"cacheMetrics\": [";
        for (std::size_t i = 0; i < r.cache_task.size(); ++i) {
            if (i) s += ", ";
            s += "{\"cacheTaskIdx\": " + std::to_string(r.cache_task[i] - 1) +
                 ", \"itemEntryIdx\": " + std::to_string(r.cache_entry[i] - 1) +
                 ", \"hitProb\": " + num(r.cache_hit[i]) + ", \"missProb\": " +
                 num(r.cache_miss[i]) + ", \"delayedProb\": " + num(r.cache_delayed[i]) +
                 ", \"readRate\": " + num(r.cache_read_rate[i]) + "}";
        }
        s += "]";
    }
    if (respt) {
        s += ", \"entryRespTimeSamples\": [";
        for (std::size_t e = 0; e < r.entry_resp_samples.size(); ++e) {
            if (e) s += ", ";
            s += "[";
            for (std::size_t j = 0; j < r.entry_resp_samples[e].size(); ++j) {
                if (j) s += ", ";
                s += num(r.entry_resp_samples[e][j]);
            }
            s += "]";
        }
        s += "]";
    }
    return s + "}\n";
}
 
}  // namespace
 
int main(int argc, char** argv) {
    if (argc < 2) {
        print_help();
        return 1;
    }
    const std::string cmd = argv[1];
    if (cmd == "-h" || cmd == "--help") {
        print_help();
        return 0;
    }
    if (cmd != "solve") {
        std::cerr << "Error: unknown command '" << cmd << "'. The only command is 'solve'.\n";
        return 1;
    }
    if (argc < 3) fail("solve requires a model file.");
 
    const std::string model_path = argv[2];
    std::string output_path;
    LdesOptions o;
    bool have_samples = false, have_events = false;
    std::vector<std::string> given;  // every flag passed, for the layered path's refusals
 
    for (int i = 3; i < argc; ++i) {
        const std::string a = argv[i];
        given.push_back(a);
        if (a == "-o") {
            output_path = need_value(argc, argv, i, a);
        } else if (a == "-s" || a == "--samples") {
            o.samples = static_cast<std::size_t>(std::atol(need_value(argc, argv, i, a).c_str()));
            have_samples = true;
        } else if (a == "-e" || a == "--maxevents") {
            o.events = static_cast<std::size_t>(std::atol(need_value(argc, argv, i, a).c_str()));
            have_events = true;
        } else if (a == "--maxtime") {
            o.timeout = std::atof(need_value(argc, argv, i, a).c_str());
        } else if (a == "--seed") {
            o.seed = std::atol(need_value(argc, argv, i, a).c_str());
        } else if (a == "--method") {
            o.method = need_value(argc, argv, i, a);
        } else if (a == "--cnvgon") {
            o.cnvgon = true;
        } else if (a == "--cnvgtol") {
            o.cnvgtol = std::atof(need_value(argc, argv, i, a).c_str());
        } else if (a == "--cnvgbatch") {
            o.cnvgbatch = std::atoi(need_value(argc, argv, i, a).c_str());
        } else if (a == "--cnvgchk") {
            o.cnvgchk = std::atoi(need_value(argc, argv, i, a).c_str());
        } else if (a == "--tranfilter") {
            o.tranfilter = need_value(argc, argv, i, a);
        } else if (a == "--mserbatch") {
            o.mserbatch = std::atoi(need_value(argc, argv, i, a).c_str());
        } else if (a == "--tranobs") {
            o.tranobs = std::atoi(need_value(argc, argv, i, a).c_str());
            if (o.tranobs <= 0) fail("--tranobs must be positive.");
        } else if (a == "--warmupfrac") {
            o.warmupfrac = std::atof(need_value(argc, argv, i, a).c_str());
        } else if (a == "--cimethod") {
            o.cimethod = need_value(argc, argv, i, a);
        } else if (a == "--obmoverlap") {
            o.obmoverlap = std::atof(need_value(argc, argv, i, a).c_str());
        } else if (a == "--ciminbatch") {
            o.ciminbatch = std::atoi(need_value(argc, argv, i, a).c_str());
        } else if (a == "--ciminobs") {
            o.ciminobs = std::atoi(need_value(argc, argv, i, a).c_str());
        } else if (a == "--spectrallowfreqfrac") {
            o.spectral_low_freq_frac = std::atof(need_value(argc, argv, i, a).c_str());
        } else if (a == "--slotted") {
            o.slotted = true;
        } else if (a == "--slotlength") {
            o.slot_length = std::atof(need_value(argc, argv, i, a).c_str());
            o.slotted = true;  // the reference lets --slotlength imply --slotted
        } else if (a == "--replications") {
            o.replications = std::atoi(need_value(argc, argv, i, a).c_str());
        } else if (a == "--numthreads") {
            o.numthreads = std::atoi(need_value(argc, argv, i, a).c_str());
        } else if (a == "--timespan") {
            const std::vector<double> t = parse_doubles(need_value(argc, argv, i, a));
            if (t.size() != 2) fail("--timespan takes T0,T1.");
            o.has_timespan = true;
            o.t0 = t[0];
            o.t1 = t[1];
        } else if (a == "--busyperiod") {
            o.busy_period_orders = std::atoi(need_value(argc, argv, i, a).c_str());
        } else if (a == "--busyperiod-subnet") {
            o.busy_period_subnets.push_back(parse_indices(need_value(argc, argv, i, a)));
        } else if (a == "--initsol") {
            o.init_sol = parse_doubles(need_value(argc, argv, i, a));
        } else if (a == "--initjobs") {
            const std::string v = need_value(argc, argv, i, a);
            std::vector<std::vector<double>> rows;
            std::size_t at = 0;
            while (at <= v.size()) {
                const std::size_t semi = std::min(v.find(';', at), v.size());
                if (semi > at) rows.push_back(parse_doubles(v.substr(at, semi - at)));
                at = semi + 1;
            }
            const std::size_t ncols = rows.empty() ? 0 : rows[0].size();
            for (std::size_t r = 0; r < rows.size(); ++r)
                if (rows[r].size() != ncols) fail("--initjobs rows must have the same length.");
            o.init_jobs = line::Matrix<double>(rows.size(), ncols, 0.0);
            for (std::size_t r = 0; r < rows.size(); ++r)
                for (std::size_t c = 0; c < ncols; ++c) o.init_jobs(r, c) = rows[r][c];
        } else if (a == "--trackjobs") {
            o.track_jobs = true;
        } else if (a == "--export-histogram") {
            o.export_histogram = true;
        } else if (a == "--trajectory") {
            o.export_trajectory = true;
        } else if (a == "--respt-samples") {
            o.export_respt = true;
        } else if (a == "--rest") {
            refuse(a, "server mode belongs to line-cli, not to this binary");
        } else if (a == "-h" || a == "--help") {
            print_help();
            return 0;
        } else {
            fail("unknown option '" + a + "'.");
        }
    }
    if (have_events && !have_samples) o.samples = o.events;
 
    try {
        // A LayeredNetwork takes the native LN engine (the twin of Solver_ssj_ln) and the layered
        // document, as LdesCLI routes it to its own. That engine reads --samples/--maxevents and
        // --seed; the flags that change WHAT is simulated or reported are refused by name, the
        // estimation and budget knobs are ignored exactly as the JAR's layered engine ignores them.
        if (line::io::is_layered_json(model_path)) {
            static const char* const network_only[] = {
                "--timespan", "--tranobs", "--trajectory", "--export-histogram", "--busyperiod",
                "--busyperiod-subnet", "--initsol", "--initjobs", "--trackjobs", "--slotted",
                "--slotlength", "--replications"};
            for (const std::string& g : given)
                for (const char* f : network_only)
                    if (g == f)
                        fail(g + " is a Network-path option; the layered engine reports the steady "
                                 "state of one replication only.");
            if (o.method != "default")
                fail("method '" + o.method + "' is not a layered LDES method; use 'default'.");
            const line::lqn::LqnStruct<double> lsn =
                line::io::read_layered_model<double>(model_path);
            const std::chrono::steady_clock::time_point t0 = std::chrono::steady_clock::now();
            const line::ldes::engine::LnResult r = line::ldes::ldes_ln_engine_solve(lsn, o);
            const double runtime =
                std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
            const std::string doc = ln_result_json(lsn, r, o, runtime, o.export_respt);
            if (output_path.empty()) {
                std::cout << doc;
            } else {
                std::ofstream out(output_path.c_str());
                if (!out) fail("cannot write '" + output_path + "'.");
                out << doc;
            }
            std::cerr << "LDES analysis [method: default; type: approximate, randomized; lang: cpp; "
                         "native LN engine] completed in " << runtime << "s. Iterations: "
                      << r.events << ".\n";
            return 0;
        }
        line::qn::Network<double> model = line::io::read_network_json<double>(model_path);
        const line::qn::NetworkStruct<double>& sn = model.get_struct();
 
        // WALL time, not std::clock(): the replications run on a thread pool, and CPU time summed over the
        // threads would report the serial cost of a run that took a fraction of it.
        const std::chrono::steady_clock::time_point t_start = std::chrono::steady_clock::now();
        const LdesResult r = line::ldes::ldes_engine_solve(sn, o);
        const double runtime =
            std::chrono::duration<double>(std::chrono::steady_clock::now() - t_start).count();
 
        const std::string doc = result_json(r, o, runtime, o.samples);
        if (output_path.empty()) {
            std::cout << doc;
        } else {
            std::ofstream out(output_path.c_str());
            if (!out) fail("cannot write '" + output_path + "'.");
            out << doc;
        }
        std::cerr << "LDES analysis [method: " << o.method
                  << "; type: approximate, randomized; lang: cpp] completed in " << runtime
                  << "s. Iterations: " << o.samples << ".\n";
        return 0;
    } catch (const line::UnsupportedError& e) {
        std::cerr << "Error: [" << e.id() << "] " << e.what() << "\n";
        return 2;
    } catch (const line::Error& e) {
        // Same exit status as before; only the vocabulary id is added to the line.
        std::cerr << "Error: [" << e.id() << "] " << e.what() << "\n";
        return 1;
    } catch (const std::exception& e) {
        std::cerr << "Error: " << e.what() << "\n";
        return 1;
    }
}