solver_ldes.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_WRAPPERS_LDES_SOLVER_LDES_H
#define LINE_SOLVERS_WRAPPERS_LDES_SOLVER_LDES_H
 
/**
 * @file
 * @ingroup line_solvers
 * Port of `SolverLDES`, the discrete-event simulator, as its C++ client.
 *
 * WHAT IS PORTED, AND WHY IT IS A CLIENT AND NOT AN ENGINE. LDES has ONE
 * implementation, `jar/src/main/java/jline/solvers/wrappers/ldes/` (~25 kLOC on SSJ),
 * and every codebase that is not the JAR reaches it the same way: MATLAB
 * (`@@SolverLDES/solveCli.m`) and native Python
 * (`wrappers/solver_ldes/solver_ldes.py`) serialize the model to `model.json`,
 * run `solve model.json -o result.json` as a subprocess and parse the result
 * document back. Nothing is marshalled in process. This file is the third such
 * client, and it is the port of SolverLDES for the same reason those two are:
 * a second simulator would answer with different numbers under the same name,
 * and "LDES agrees across the codebases" -- which the parity harness checks and
 * `CLAUDE.md` states as an invariant -- is a property of there being ONE engine.
 *
 * THE PRIMARY ENTRY IS A DOCUMENT, not a model object, and that is a deliberate
 * difference from the two other clients. They hold a live model and serialize
 * it; `solver_ldes_text` / `solver_ldes_file` are given the `model.json`
 * DOCUMENT and forward its bytes UNCHANGED. There is a writer
 * (`io/network_writer.h`), so a round trip through `qn::NetworkStruct<T>` is
 * possible -- `solver_ldes` below is exactly that, for a caller who built the
 * model in C++ -- but the reader that would have to parse the document first is
 * deliberately scoped to the subset the analytical solvers read, so making the
 * trip MANDATORY would degrade exactly the models LDES exists for (a cache with
 * retrieval, an SPN, a polling server, a G-network signal). Forwarding the bytes
 * keeps the client lossless on every model the engine accepts, including those
 * this port cannot itself represent.
 *
 * WHICH BINARY RUNS. `ldes_runners` reproduces `getLdesRunners` and
 * `_build_cli_runners`: the AOT GraalVM binary `common/ldes` first, for its
 * startup time, then `<java> -jar common/ldes.jar`, and each is tried in turn
 * because the AOT image lacks reflective features the JVM has (the fork-join
 * MMT transform deep-copies the model by Java serialization). THE ORDER IS
 * REVERSED, not merely extended, when the run uses something the prebuilt image
 * predates -- `--respt-samples`, `--initsol`, a cache cost cap, a marked class,
 * a load/class-dependent scaling -- because on that image the flag is not
 * refused but IGNORED, and an ignored cost cap silently simulates the uncapped
 * cache. The two other clients inspect the live model object for those; this
 * one inspects the document, which is the same information at the point where
 * it actually crosses.
 *
 * STATE PROBABILITIES ARE READ FROM THE HISTOGRAM, not from the trajectory.
 * `ldes_prob_aggr` / `ldes_prob_sys_aggr` take the target state as an ARGUMENT
 * and weigh it against `stateHistogram`, the exact residence time of every
 * integer joint state the run visited. The no-argument form the reference also
 * offers -- "the probability of the model's CURRENT state" -- has no counterpart
 * here because `qn::NetworkStruct` carries no current state row (only a Cache's
 * `initstate`, a Place's `initmarking` and the `stateprior`/`statespace` pair).
 * The DEFAULT INITIAL state is a different thing and is derivable, by the same
 * `ctmc::analyzer_detail::default_init_state` the `-a prob` arms of CTMC, MVA
 * and NC already answer over; the LDES arm of the CLI still refuses `-a prob`
 * on the older premise that no state exists at all, which is a reachability gap
 * in `cpp/src/cli/line_cli.cpp` rather than a gap in this client.
 * The three other clients answer BOTH forms off the transient `QNt` series,
 * which is a trajectory of interval MEANS, and so report a near-zero
 * probability for a state the chain mostly occupies; see BUG-96 and the note on
 * `LdesResult::QNt`.
 *
 * WHAT THIS CLIENT DOES NOT PORT, and why it refuses instead of approximating:
 *   LayeredNetwork         MATLAB simulates an LQN through the Java ensemble
 *                          backend (`self.obj`), an in-process path with no JSON
 *                          interface. Both `common/ldes` (since 2026-10-31) and
 *                          `LdesCLI` now solve a LayeredNetwork document, so this
 *                          is a client gap, not an engine one: `line-cli -s ldes`
 *                          reaches the C++ LN engine in process instead.
 *   server breakdowns      `linemodel_save` does not serialize `setBreakdown` at
 *                          all, so a breakdown cannot reach this client: the
 *                          refusal the Python client performs is on the OBJECT,
 *                          before serialization, and has no wire counterpart.
 *
 * WHERE THE ENGINE IS, and whether there is one, lives in `ldes_probe.h`:
 * SolverAUTO needs that answer -- LDES leads several of the reference's
 * rankings, and a chooser that named it on a host without the engine would be
 * proposing something that cannot run -- and must not take on this file's HTTP
 * transport and result parser to get it.
 */
 
#include "line/util/line_console.h"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <fstream>
#include <limits>
#include <sstream>
#include <string>
#include <utility>
#include <vector>
 
#include <sys/stat.h>
#include <sys/utsname.h>
#include <unistd.h>
 
#include "json.hpp"
#include "line/io/network_writer.h"
#include "line/solvers/fj_metrics.h"
#include "line/solvers/wrappers/ldes/ldes_options.h"
#include "line/solvers/wrappers/ldes/ldes_probe.h"
#include "line/util/error.h"
#include "line/util/http.h"
#include "line/util/matrix.h"
#include "line/util/subprocess.h"
#include "line/util/tempdir.h"
 
namespace line {
namespace ldes {
 
namespace detail {
 
using Json = nlohmann::json;
 
/**
 * The shortest decimal that reads back as `v`, which is what Python's `repr`
 * emits on the same flag. A fixed `%.10g`, as `solveCli.m` uses, would round a
 * tolerance the caller set exactly, and the engine parses whatever it is given.
 */
inline std::string shortest(double v) {
    char buf[64];
    for (int prec = 15; prec <= 17; ++prec) {
        std::snprintf(buf, sizeof(buf), "%.*g", prec, v);
        if (std::strtod(buf, nullptr) == v) return std::string(buf);
    }
    return std::string(buf);
}
 
/** A JSON array of arrays as a matrix, with `null` read back as NaN. */
inline Matrix<double> mat(const Json& v) {
    if (v.is_null() || !v.is_array() || v.empty()) return Matrix<double>();
    if (!v[0].is_array()) {
        Matrix<double> m(1, v.size());
        for (std::size_t j = 0; j < v.size(); ++j)
            m(0, j) = v[j].is_null() ? std::numeric_limits<double>::quiet_NaN()
                                     : v[j].get<double>();
        return m;
    }
    std::size_t cols = 0;
    for (std::size_t i = 0; i < v.size(); ++i)
        if (v[i].is_array() && v[i].size() > cols) cols = v[i].size();
    Matrix<double> m(v.size(), cols, std::numeric_limits<double>::quiet_NaN());
    for (std::size_t i = 0; i < v.size(); ++i) {
        if (!v[i].is_array()) continue;
        for (std::size_t j = 0; j < v[i].size(); ++j)
            if (!v[i][j].is_null()) m(i, j) = v[i][j].get<double>();
    }
    return m;
}
 
/** A JSON array as a vector of doubles, `null` read back as NaN. */
inline std::vector<double> vec(const Json& v) {
    std::vector<double> out;
    if (!v.is_array()) return out;
    out.reserve(v.size());
    for (std::size_t i = 0; i < v.size(); ++i)
        out.push_back(v[i].is_null() ? std::numeric_limits<double>::quiet_NaN()
                                     : v[i].get<double>());
    return out;
}
 
inline const Json& field(const Json& o, const char* key) {
    static const Json null_value;
    return o.is_object() && o.contains(key) ? o[key] : null_value;
}
 
/** Reads a whole file, or throws naming it. */
inline std::string read_file(const std::string& path) {
    std::ifstream in(path.c_str(), std::ios::binary);
    if (!in) throw InputError("SolverLDES: cannot open the model document '" + path + "'");
    std::ostringstream ss;
    ss << in.rdbuf();
    return ss.str();
}
 
inline void write_file(const std::string& path, const std::string& text) {
    std::ofstream out(path.c_str(), std::ios::binary);
    if (!out) throw InputError("SolverLDES: cannot write '" + path + "'");
    out << text;
}
 
/**
 * True when the DOCUMENT uses something the native engine would MIS-SERVE.
 *
 * The scan is on the bytes that actually cross rather than on a live model, and
 * a substring test is deliberate: every key is unique enough that a false
 * positive only reorders the runners, which costs startup time and changes no
 * number, while a false negative would silently simulate a different model.
 *
 * THE BAR IS "SILENTLY IGNORED", NOT "DIFFERS". A key belongs here only while
 * the native engine would accept it and then not honour it, because that is the
 * failure this guard exists to prevent: an ignored cost cap simulates the
 * uncapped cache and says nothing. A feature the native engine REFUSES by name
 * needs no entry -- the runner list falls through to the jar on failure -- and
 * neither does one it honours.
 *
 * `markedClasses` was removed on 2026-09-12, for three reasons together. The
 * native engine has routed marks since 2026-09-11 and released batches since
 * 2026-09-12, so it no longer ignores them. `common/ldes` is no longer the
 * retired GraalVM image of the jar but a C++ binary built from this same tree
 * (`add_executable(ldes src/cli/ldes_cli.cpp)`), so it cannot predate a document
 * this tree emits. And neither other client ever carried the entry -- MATLAB's
 * `getLdesRunners` reverses for nothing and python's `_build_cli_runners` only
 * for `--respt-samples` -- so keeping it meant
 * one marked model ran on TWO different engines depending on which codebase's
 * client was asked, which is the one thing "there is one LDES" is supposed to
 * rule out.
 *
 * `costCaps` was removed on 2026-10-31: the native engine enforces the caps and
 * now also reports `hitList`, `itemProb` and `listCost` for them.
 *
 * Note that digit-for-digit agreement is NOT the bar, and could not be: the two
 * engines part company on any model whose routing consumes extra draws, which
 * `CLAUDE.md` states outright. `--initsol` and `--busyperiod` happened to clear
 * that higher bar when they were pruned; marks do not have to.
 */
inline bool document_postdates_native(const std::string& doc) {
    static const char* keys[] = {"\"classDependence\"",
                                 "\"loadDependence\"", "\"jointDependence\""};
    for (std::size_t i = 0; i < sizeof(keys) / sizeof(keys[0]); ++i)
        if (doc.find(keys[i]) != std::string::npos) return true;
    return false;
}
 
}  // namespace detail
 
/**
 * The engine flags of one run, after `solve <model> -o <result>`.
 *
 * ONE MAPPING for the subprocess and the REST paths, as in both other clients:
 * the server takes the same long-form flags verbatim, and deriving them twice
 * is how the two transports start disagreeing. A knob is emitted only when it
 * differs from the engine default, so a default run yields the minimal command
 * line an older AOT image still parses -- except `-s` and `--seed`, which are
 * always emitted: the CLI's seed default is -1 (random), so a silent seed makes
 * the run irreproducible and unlike the MATLAB and Python clients.
 */
inline std::vector<std::string> ldes_flags(const LdesOptions& o,
                                           const std::vector<std::string>& extra) {
    std::vector<std::string> f;
    const std::size_t budget = o.events ? o.events : o.samples;
    f.push_back("-s");
    f.push_back(std::to_string(budget));
    f.push_back("--seed");
    f.push_back(std::to_string(o.seed));
    if (!o.method.empty() && o.method != "default") {
        f.push_back("--method");
        f.push_back(o.method);
    }
    if (o.cnvgon) {
        f.push_back("--cnvgon");
        if (o.cnvgtol != 0.05) {
            f.push_back("--cnvgtol");
            f.push_back(detail::shortest(o.cnvgtol));
        }
        if (o.cnvgbatch != 20) {
            f.push_back("--cnvgbatch");
            f.push_back(std::to_string(o.cnvgbatch));
        }
        if (o.cnvgchk != 0) {
            f.push_back("--cnvgchk");
            f.push_back(std::to_string(o.cnvgchk));
        }
    }
    if (o.tranfilter != "mser5") {
        f.push_back("--tranfilter");
        f.push_back(o.tranfilter);
    }
    if (o.warmupfrac != 0.2) {
        f.push_back("--warmupfrac");
        f.push_back(detail::shortest(o.warmupfrac));
    }
    if (o.tranobs > 0) {
        f.push_back("--tranobs");
        f.push_back(std::to_string(o.tranobs));
    }
    if (o.mserbatch != 5) {
        f.push_back("--mserbatch");
        f.push_back(std::to_string(o.mserbatch));
    }
    if (o.cimethod != "obm") {
        f.push_back("--cimethod");
        f.push_back(o.cimethod);
    }
    if (o.obmoverlap != 0.5) {
        f.push_back("--obmoverlap");
        f.push_back(detail::shortest(o.obmoverlap));
    }
    if (o.ciminbatch != 10) {
        f.push_back("--ciminbatch");
        f.push_back(std::to_string(o.ciminbatch));
    }
    if (o.ciminobs != 100) {
        f.push_back("--ciminobs");
        f.push_back(std::to_string(o.ciminobs));
    }
    if (o.spectral_low_freq_frac != 0.25) {
        f.push_back("--spectrallowfreqfrac");
        f.push_back(detail::shortest(o.spectral_low_freq_frac));
    }
    // --slotlength implies --slotted on the CLI side; both are emitted when the
    // length is non-default so the command line states the intent explicitly.
    if (o.slotted) {
        f.push_back("--slotted");
        if (o.slot_length != 1.0) {
            f.push_back("--slotlength");
            f.push_back(detail::shortest(o.slot_length));
        }
    }
    // method='parallel' (alias 'para') is a replication count, not a second engine: the caller's
    // count when set, 8 otherwise, as solveCli.m resolves it. The native engine reads only --replications.
    int reps = o.replications;
    if ((o.method == "parallel" || o.method == "para") && reps <= 1) reps = 8;
    if (reps > 1) {
        f.push_back("--replications");
        f.push_back(std::to_string(reps));
        if (o.numthreads > 0) {
            f.push_back("--numthreads");
            f.push_back(std::to_string(o.numthreads));
        }
    }
    if (o.has_timespan) {
        f.push_back("--timespan");
        f.push_back(detail::shortest(o.t0) + "," + detail::shortest(o.t1));
    }
    if (o.timeout > 0.0 && std::isfinite(o.timeout)) {
        f.push_back("--maxtime");
        f.push_back(detail::shortest(o.timeout));
    }
    if (o.busy_period_orders > 0) {
        f.push_back("--busyperiod");
        f.push_back(std::to_string(o.busy_period_orders));
        for (std::size_t i = 0; i < o.busy_period_subnets.size(); ++i) {
            // station indexes cross the wire zero-based, as the engine holds them
            std::string v;
            for (std::size_t t = 0; t < o.busy_period_subnets[i].size(); ++t)
                v += (t ? "," : "") + std::to_string(o.busy_period_subnets[i][t]);
            f.push_back("--busyperiod-subnet");
            f.push_back(v);
        }
    }
    if (!o.init_sol.empty()) {
        std::string v;
        for (std::size_t i = 0; i < o.init_sol.size(); ++i)
            v += (i ? "," : "") + detail::shortest(o.init_sol[i]);
        f.push_back("--initsol");
        f.push_back(v);
    }
    if (o.init_jobs.rows() > 0) {
        std::string v;
        for (std::size_t r = 0; r < o.init_jobs.rows(); ++r)
            for (std::size_t c = 0; c < o.init_jobs.cols(); ++c)
                v += (r && !c ? ";" : (c ? "," : "")) + detail::shortest(o.init_jobs(r, c));
        f.push_back("--initjobs");
        f.push_back(v);
    }
    if (o.track_jobs) f.push_back("--trackjobs");
    for (std::size_t i = 0; i < extra.size(); ++i) f.push_back(extra[i]);
    return f;
}
 
/** One runner: the argv prefix that runs the engine, and the name of its image. */
struct LdesRunner {
    std::vector<std::string> argv;  ///< up to and including "solve"
    std::string engine;             ///< "native" or "jar"
};
 
/**
 * The runners to try, in order (see the file header for why there are two and
 * when the order flips).
 *
 * @param doc   the model document, scanned for features the AOT image predates
 * @param flags the resolved flags, scanned for the same reason
 * @throws UnsupportedError when no engine can run on this host
 */
inline std::vector<LdesRunner> ldes_runners(const std::string& doc,
                                            const std::vector<std::string>& flags) {
    const std::string& dir = ldes_engine_dir();
    std::vector<LdesRunner> out;
    const std::string native = detail::native_ldes_path(dir);
    if (!native.empty()) {
        LdesRunner r;
        r.argv.push_back(native);
        r.argv.push_back("solve");
        r.engine = "native";
        out.push_back(r);
    }
    const std::string java = detail::find_java();
    const std::string jar = dir.empty() ? std::string() : dir + "/ldes.jar";
    if (!java.empty() && detail::is_file(jar)) {
        LdesRunner r;
        r.argv.push_back(java);
        r.argv.push_back("-jar");
        r.argv.push_back(jar);
        r.argv.push_back("solve");
        r.engine = "jar";
        out.push_back(r);
    }
    if (out.empty()) throw ExternalError(ErrorKind{"LDES"}, detail::no_backend_message(dir));
    if (out.size() > 1) {
        // --initsol and --busyperiod used to be here, because the retired
        // GraalVM image predated them; the C++ engine honours both and agrees
        // with the jar to the last digit, so only --respt-samples remains.
        bool prefer_jar = detail::document_postdates_native(doc);
        for (std::size_t i = 0; i < flags.size() && !prefer_jar; ++i)
            prefer_jar = flags[i] == "--respt-samples";
        if (prefer_jar) std::swap(out[0], out[1]);
    }
    return out;
}
 
/**
 * Solves through an LDES REST server and returns its `result` document.
 *
 * The payload is the model text plus the same flag vector the subprocess would
 * have been given, so a fixed seed gives the same numbers on both transports.
 */
inline detail::Json ldes_solve_rest(const std::string& base_url, const std::string& doc,
                                    const std::vector<std::string>& flags, double timeout) {
    std::string url = base_url;
    while (!url.empty() && url[url.size() - 1] == '/') url.erase(url.size() - 1);
    if (url.size() < 6 || url.compare(url.size() - 6, 6, "/solve") != 0) url += "/api/v1/solve";
 
    detail::Json payload = detail::Json::object();
    detail::Json model = detail::Json::object();
    model["content"] = doc;
    model["base64"] = false;
    payload["model"] = model;
    detail::Json fl = detail::Json::array();
    for (std::size_t i = 0; i < flags.size(); ++i) fl.push_back(flags[i]);
    payload["flags"] = fl;
 
    const int millis = (timeout > 0.0 && std::isfinite(timeout))
                           ? static_cast<int>((timeout + 30.0) * 1000.0)
                           : 600000;
    const http::Response resp = http::post_json(url, payload.dump(), millis);
    if (resp.body.empty())
        throw ExternalError(ErrorKind{"LDES"}, "SolverLDES: the REST server at " + url + " returned HTTP " +
                           std::to_string(resp.status) + " with no body");
    detail::Json out = detail::Json::parse(resp.body, nullptr, false);
    if (out.is_discarded() || !out.is_object())
        throw ExternalError(ErrorKind{"LDES"}, "SolverLDES: the REST server at " + url +
                           " returned a non-JSON body: " + resp.body);
    if (!out.contains("status") || out["status"].get<std::string>() != "ok") {
        const std::string msg = out.contains("message") ? out["message"].get<std::string>()
                                                        : std::string("unspecified error");
        const std::string err = out.contains("stderr") ? out["stderr"].get<std::string>()
                                                       : std::string();
        throw ExternalError(ErrorKind{"LDES"}, "SolverLDES: the REST solve failed: " + msg +
                           (err.empty() ? "" : " (engine stderr: " + err + ")"));
    }
    return out["result"];
}
 
/**
 * Parses one `ldes-result` document.
 *
 * `format` is CHECKED, not assumed: the runner falls through on failure, and a
 * runner that wrote some other JSON to the output path would otherwise be read
 * as a result full of absent metrics rather than as the failure it is.
 */
inline LdesResult parse_ldes_result(const detail::Json& d) {
    if (d.is_object() && d.contains("error"))
        throw NumericError(ErrorKind{"InvalidResult"}, "SolverLDES: the engine reported an error: " +
                           d["error"].get<std::string>());
    if (!d.is_object() || !d.contains("format") || d["format"].get<std::string>() != "ldes-result")
        throw NumericError(ErrorKind{"InvalidResult"},
            "SolverLDES: the engine wrote a document that is not an ldes-result; the run did not "
            "produce metrics");
 
    LdesResult r;
    r.method = d.contains("method") && !d["method"].is_null() ? d["method"].get<std::string>()
                                                              : std::string("default");
    if (d.contains("runtime")) r.runtime = d["runtime"].get<double>();
    if (d.contains("converged")) r.converged = d["converged"].get<bool>();
    if (d.contains("stoppingReason") && !d["stoppingReason"].is_null())
        r.stopping_reason = d["stoppingReason"].get<std::string>();
    // The engine's console is not shown, so its warnings are re-raised here.
    if (d.contains("warnings") && d["warnings"].is_array())
        for (const auto& w : d["warnings"]) {
            r.warnings.push_back(w.get<std::string>());
            std::fprintf(stderr, "[SolverLDES] Warning: %s\n", r.warnings.back().c_str());
        }
    if (d.contains("convergenceBatches")) r.convergence_batches = d["convergenceBatches"].get<long>();
    if (d.contains("totalSimulatedEvents"))
        r.total_simulated_events = d["totalSimulatedEvents"].get<long long>();
 
    const detail::Json& dim = detail::field(d, "dimensions");
    if (dim.is_object()) {
        if (dim.contains("nstations")) r.nstations = dim["nstations"].get<std::size_t>();
        if (dim.contains("nclasses")) r.nclasses = dim["nclasses"].get<std::size_t>();
        if (dim.contains("nchains")) r.nchains = dim["nchains"].get<std::size_t>();
        if (dim.contains("stationNames"))
            for (std::size_t i = 0; i < dim["stationNames"].size(); ++i)
                r.station_names.push_back(dim["stationNames"][i].get<std::string>());
        if (dim.contains("classNames"))
            for (std::size_t i = 0; i < dim["classNames"].size(); ++i)
                r.class_names.push_back(dim["classNames"][i].get<std::string>());
    }
 
    const detail::Json& m = detail::field(d, "metrics");
    r.QN = detail::mat(detail::field(m, "QN"));
    r.UN = detail::mat(detail::field(m, "UN"));
    r.RN = detail::mat(detail::field(m, "RN"));
    r.TN = detail::mat(detail::field(m, "TN"));
    r.AN = detail::mat(detail::field(m, "AN"));
    r.WN = detail::mat(detail::field(m, "WN"));
    r.CN = detail::mat(detail::field(m, "CN"));
    r.XN = detail::mat(detail::field(m, "XN"));
    r.DropRateJoin = detail::mat(detail::field(m, "DropRateJoin"));
 
    const detail::Json& ci = detail::field(d, "confidenceIntervals");
    r.QNCI = detail::mat(detail::field(ci, "QNCI"));
    r.UNCI = detail::mat(detail::field(ci, "UNCI"));
    r.RNCI = detail::mat(detail::field(ci, "RNCI"));
    r.TNCI = detail::mat(detail::field(ci, "TNCI"));
    r.ANCI = detail::mat(detail::field(ci, "ANCI"));
    r.WNCI = detail::mat(detail::field(ci, "WNCI"));
 
    const detail::Json& rp = detail::field(d, "relativePrecision");
    r.QNRelPrec = detail::mat(detail::field(rp, "QNRelPrec"));
    r.UNRelPrec = detail::mat(detail::field(rp, "UNRelPrec"));
    r.RNRelPrec = detail::mat(detail::field(rp, "RNRelPrec"));
    r.TNRelPrec = detail::mat(detail::field(rp, "TNRelPrec"));
 
    const detail::Json& sc = detail::field(d, "sampleCounts");
    r.QNSamples = detail::mat(detail::field(sc, "QNSamples"));
    r.UNSamples = detail::mat(detail::field(sc, "UNSamples"));
    r.RNSamples = detail::mat(detail::field(sc, "RNSamples"));
    r.TNSamples = detail::mat(detail::field(sc, "TNSamples"));
 
    const detail::Json& fcr = detail::field(d, "fcr");
    if (fcr.is_object()) {
        if (fcr.contains("nregions")) r.nregions = fcr["nregions"].get<std::size_t>();
        r.QNfcr = detail::mat(detail::field(fcr, "QNfcr"));
        r.UNfcr = detail::mat(detail::field(fcr, "UNfcr"));
        r.RNfcr = detail::mat(detail::field(fcr, "RNfcr"));
        r.TNfcr = detail::mat(detail::field(fcr, "TNfcr"));
        r.ANfcr = detail::mat(detail::field(fcr, "ANfcr"));
        r.WNfcr = detail::mat(detail::field(fcr, "WNfcr"));
        r.WeightNfcr = detail::mat(detail::field(fcr, "WeightNfcr"));
        r.MemOccNfcr = detail::mat(detail::field(fcr, "MemOccNfcr"));
        r.DropRateNfcr = detail::mat(detail::field(fcr, "DropRateNfcr"));
    }
 
    // Fork-join section metrics: rows named by "forks"/"joins", mapped onto the model's pairs by ldes_fj_metrics.
    const detail::Json& fjb = detail::field(d, "forkJoin");
    if (fjb.is_object()) {
        const detail::Json& fk = detail::field(fjb, "forks");
        const detail::Json& jn = detail::field(fjb, "joins");
        for (std::size_t p = 0; fk.is_array() && p < fk.size(); ++p)
            r.fj_forks.push_back(fk[p].is_string() ? fk[p].get<std::string>() : std::string());
        for (std::size_t p = 0; jn.is_array() && p < jn.size(); ++p)
            r.fj_joins.push_back(jn[p].is_string() ? jn[p].get<std::string>() : std::string());
        r.FJQLen = detail::mat(detail::field(fjb, "FJQLen"));
        r.FJRespT = detail::mat(detail::field(fjb, "FJRespT"));
        r.FJTput = detail::mat(detail::field(fjb, "FJTput"));
    }
 
    const detail::Json& imp = detail::field(d, "impatience");
    if (imp.is_object()) {
        r.renegedCustomers = detail::mat(detail::field(imp, "renegedCustomers"));
        r.avgRenegingWaitTime = detail::mat(detail::field(imp, "avgRenegingWaitTime"));
        r.renegingRate = detail::mat(detail::field(imp, "renegingRate"));
        r.balkedCustomers = detail::mat(detail::field(imp, "balkedCustomers"));
        r.balkingProbability = detail::mat(detail::field(imp, "balkingProbability"));
        r.retriedCustomers = detail::mat(detail::field(imp, "retriedCustomers"));
        r.retrialDropped = detail::mat(detail::field(imp, "retrialDropped"));
        r.avgOrbitSize = detail::mat(detail::field(imp, "avgOrbitSize"));
    }
 
    const detail::Json& cm = detail::field(d, "cacheMetrics");
    if (cm.is_object())
        for (detail::Json::const_iterator it = cm.begin(); it != cm.end(); ++it) {
            LdesCacheMetrics c;
            c.hit = detail::mat(detail::field(it.value(), "hit"));
            c.delayed = detail::mat(detail::field(it.value(), "delayed"));
            c.miss = detail::mat(detail::field(it.value(), "miss"));
            c.latency = detail::mat(detail::field(it.value(), "latency"));
            c.hitList = detail::mat(detail::field(it.value(), "hitList"));
            c.itemProb = detail::mat(detail::field(it.value(), "itemProb"));
            c.listCost = detail::mat(detail::field(it.value(), "listCost"));
            r.cache_metrics[it.key()] = c;
        }
 
    const detail::Json& bp = detail::field(d, "busyPeriods");
    if (bp.is_object()) {
        const detail::Json& targets = detail::field(bp, "targets");
        for (std::size_t i = 0; i < targets.size(); ++i) {
            LdesResult::BusyPeriodTarget t;
            const detail::Json& tj = targets[i];
            if (tj.contains("name") && tj["name"].is_string())
                t.name = tj["name"].get<std::string>();
            const detail::Json& st = detail::field(tj, "stations");
            for (std::size_t k = 0; k < st.size(); ++k)
                t.stations.push_back(static_cast<std::size_t>(st[k].get<double>()));
            if (tj.contains("class") && tj["class"].is_number())
                t.job_class = static_cast<int>(tj["class"].get<double>());
            t.mean = detail::vec(detail::field(tj, "mean"));
            t.count = detail::vec(detail::field(tj, "count"));
            r.busy_periods.push_back(t);
        }
    }
 
    const detail::Json& js = detail::field(d, "jobState");
    if (js.is_object()) {
        r.finalJobs = detail::mat(detail::field(js, "finalJobs"));
        r.jobDepartures = detail::mat(detail::field(js, "jobDepartures"));
    }
 
    const detail::Json& hist = detail::field(d, "stateHistogram");
    if (hist.is_object()) {
        r.histogram_space = detail::mat(detail::field(hist, "space"));
        r.histogram_time = detail::mat(detail::field(hist, "time"));
        r.traj_space = detail::mat(detail::field(hist, "trajSpace"));
        r.traj_time = detail::mat(detail::field(hist, "trajTime"));
    }
 
    // `respTimeSamples` appears at TOP LEVEL under --respt-samples and inside
    // `transient` under --trajectory; the two are the same measurement and
    // whichever is present is read.
    const detail::Json& tran = detail::field(d, "transient");
    if (tran.is_object()) {
        r.t = detail::vec(detail::field(tran, "t"));
        const char* keys[3] = {"QNt", "UNt", "TNt"};
        std::vector<std::vector<Matrix<double>>>* dst[3] = {&r.QNt, &r.UNt, &r.TNt};
        for (int q = 0; q < 3; ++q) {
            const detail::Json& blk = detail::field(tran, keys[q]);
            if (!blk.is_array()) continue;
            for (std::size_t i = 0; i < blk.size(); ++i) {
                std::vector<Matrix<double>> row;
                for (std::size_t k = 0; k < blk[i].size(); ++k)
                    row.push_back(detail::mat(blk[i][k]));
                dst[q]->push_back(row);
            }
        }
    }
    const detail::Json& rts = tran.is_object() && tran.contains("respTimeSamples")
                                  ? tran["respTimeSamples"]
                                  : detail::field(d, "respTimeSamples");
    if (rts.is_array())
        for (std::size_t i = 0; i < rts.size(); ++i) {
            std::vector<std::vector<double>> row;
            for (std::size_t k = 0; k < rts[i].size(); ++k) row.push_back(detail::vec(rts[i][k]));
            r.respTimeSamples.push_back(row);
        }
    return r;
}
 
/**
 * Runs one LDES simulation on a `model.json` DOCUMENT and parses its result.
 *
 * The document is written into a private scratch directory, the engine is run
 * there, and the directory is removed on both the success and the failure path
 * (`TempDir`'s destructor), because a wrapper that leaks one directory per
 * failed solve leaks silently.
 *
 * THE RUNNERS ARE TRIED IN ORDER and the LAST failure is reported in full,
 * merged stdout and stderr, because the engine states why it refused there and
 * an exit code alone is not a diagnosis.
 *
 * @param doc         the model.json text, forwarded byte for byte
 * @param o           the run's knobs
 * @param extra_flags flags an analysis adds (--trajectory, --export-histogram,
 *                    --respt-samples)
 * @return the parsed result, carrying which image produced it
 */
inline LdesResult solver_ldes_text(const std::string& doc, const LdesOptions& o,
                                   const std::vector<std::string>& extra_flags =
                                       std::vector<std::string>()) {
    const std::vector<std::string> flags = ldes_flags(o, extra_flags);
    util::TempDir tmp("ldes");
    const std::string model_path = tmp.file("model.json");
    const std::string result_path = tmp.file("result.json");
    line::util::LineConsole::step("serializing the model to JSON for the LDES engine");
    detail::write_file(model_path, doc);
 
    if (!o.rest_url.empty()) {
        LdesResult r = parse_ldes_result(ldes_solve_rest(o.rest_url, doc, flags, o.timeout));
        r.engine = "rest";
        return r;
    }
 
    const std::vector<LdesRunner> runners = ldes_runners(doc, flags);
    // The cooperative --maxtime stops the event loop; this is the HARD bound and
    // leaves the engine 30 s of grace to write its result and exit, else no
    // result file is produced at all. An infinite budget falls back to 600 s.
    const int hard = (o.timeout > 0.0 && std::isfinite(o.timeout))
                         ? static_cast<int>(o.timeout + 30.0)
                         : 600;
    std::string last_err;
    for (std::size_t i = 0; i < runners.size(); ++i) {
        std::vector<std::string> argv = runners[i].argv;
        argv.push_back(model_path);
        argv.push_back("-o");
        argv.push_back(result_path);
        for (std::size_t j = 0; j < flags.size(); ++j) argv.push_back(flags[j]);
        std::remove(result_path.c_str());
        if (o.verbose) {
            std::string line;
            for (std::size_t j = 0; j < argv.size(); ++j) line += (j ? " " : "") + argv[j];
            std::fprintf(stderr, "SolverLDES command: %s\n", line.c_str());
        }
        line::util::LineConsole::step("running the LDES engine (%s) as a subprocess",
                                      runners[i].engine.c_str());
        const util::ProcResult p = util::capture(argv, hard, true);
        if (p.timedOut) {
            LdesResult r;
            r.timed_out = true;
            r.stopping_reason = "max_time";
            r.engine = runners[i].engine;
            return r;
        }
        if (p.exitCode == 0 && detail::is_file(result_path)) {
            line::util::LineConsole::step("parsing the LDES result document");
            const std::string text = detail::read_file(result_path);
            detail::Json parsed = detail::Json::parse(text, nullptr, false);
            if (parsed.is_discarded())
                throw ExternalError(ErrorKind{"LDES"}, "SolverLDES: the engine wrote a result that is not JSON");
            LdesResult r = parse_ldes_result(parsed);
            r.engine = runners[i].engine;
            return r;
        }
        last_err = "exit code " + std::to_string(p.exitCode) + ": " + util::trim(p.out);
    }
    throw ExternalError(ErrorKind{"LDES"}, "SolverLDES: the engine failed on all " +
                       std::to_string(runners.size()) + " runner(s); " + last_err);
}
 
/** The same, reading the document from a file. */
inline LdesResult solver_ldes_file(const std::string& path, const LdesOptions& o,
                                   const std::vector<std::string>& extra_flags =
                                       std::vector<std::string>()) {
    return solver_ldes_text(detail::read_file(path), o, extra_flags);
}
 
/**
 * One constraint of a joint-state query: a station and the per-class job counts
 * it is required to hold.
 *
 * `station` is 1-BASED, matching `sn.nodeToStation`. `counts` is read up to
 * `nclasses` entries; a shorter vector constrains only the classes it names,
 * which is how a caller asks for "two class-1 jobs here, any number of the
 * rest" without enumerating the remainder.
 */
struct LdesStateQuery {
    std::size_t station = 0;
    std::vector<double> counts;
};
 
/**
 * Residence-time probability of an aggregate joint state, from a parsed result.
 *
 * WHY THE HISTOGRAM AND NOT THE TRAJECTORY. `stateHistogram` is the exact
 * residence time of every integer joint state the run visited, keyed on the
 * state itself (`Solver_ssj.updateRewardStats` accumulates it per marking), so
 * P(state) = t(state) / sum(t) is an unbiased estimate that is exact on the
 * sampled path. The transient `QNt` series is a sequence of INTERVAL MEANS of
 * the queue length; comparing one against an integer state matches only where a
 * bucket mean happens to land on an integer, which is why the three clients
 * that estimated a probability that way reported a near-zero number for a state
 * the chain occupies most of the time (BUG-96).
 *
 * The layout is the aggregated station-major, class-minor one
 * `ctmc_state_space_aggr` builds: column `(i-1)*nclasses + k` is the number of
 * class-k jobs at station i. THE RESOLUTION IS PER CLASS, not per phase: the
 * engine records the integer queue lengths, so a phase-resolved query cannot be
 * answered from this document and the caller must aggregate first.
 *
 * @param r        a result produced with `--export-histogram`
 * @param nclasses the model's class count, which fixes the column stride
 * @param query    the stations to constrain, and to what
 * @return the residence-time fraction, 0 when the state was never visited
 * @throws NumericError when the run carried no histogram (the flag was omitted)
 */
inline double ldes_prob_from_histogram(const LdesResult& r, std::size_t nclasses,
                                       const std::vector<LdesStateQuery>& query) {
    if (r.histogram_space.rows() == 0 || r.histogram_time.rows() == 0)
        throw NumericError(ErrorKind{"Degenerate"},
            "SolverLDES: this result carries no state histogram, so no state probability can be "
            "read from it; the run has to pass --export-histogram");
    if (nclasses == 0)
        throw InputError("SolverLDES: a state probability needs the model's class count");
    const std::size_t nstates = r.histogram_space.rows();
    double total = 0.0;
    for (std::size_t s = 0; s < nstates; ++s) total += r.histogram_time(s, 0);
    if (!(total > 0.0)) return 0.0;
 
    const std::size_t ncols = r.histogram_space.cols();
    for (std::size_t q = 0; q < query.size(); ++q) {
        if (query[q].station == 0)
            throw InputError("SolverLDES: a state query names station 0; stations are 1-based");
        const std::size_t last = (query[q].station - 1) * nclasses + nclasses;
        if (last > ncols)
            throw InputError(
                "SolverLDES: station " + std::to_string(query[q].station) +
                " is past the end of the state histogram, which holds " + std::to_string(ncols) +
                " columns for " + std::to_string(nclasses) + " classes");
    }
 
    double matched = 0.0;
    for (std::size_t s = 0; s < nstates; ++s) {
        bool ok = true;
        for (std::size_t q = 0; q < query.size() && ok; ++q) {
            const std::size_t base = (query[q].station - 1) * nclasses;
            const std::size_t n = std::min(nclasses, query[q].counts.size());
            for (std::size_t k = 0; k < n; ++k)
                if (std::fabs(r.histogram_space(s, base + k) - query[q].counts[k]) > 1e-9) {
                    ok = false;
                    break;
                }
        }
        if (ok) matched += r.histogram_time(s, 0);
    }
    return matched / total;
}
 
/**
 * Port of `getCdfRespT`: the EMPIRICAL response time CDF of one (station,
 * class) pair, built from the per-job samples the engine exports under
 * `--respt-samples`.
 *
 * The returned matrix is (n x 2) with columns [F(t), t], NOT [t, F(t)]. That
 * order is the convention every `getCdfRespT` follows across the codebases
 * (MATLAB `@SolverLDES/getCdfRespT.m`, Python `SolverJMT.getCdfRespT`), and it
 * is the REVERSE of the transient `getTranCdfRespT`, which is a different
 * method with a different contract. Swapping them yields a CDF that reads as a
 * time axis and vice versa, with no error anywhere.
 *
 * A pair the run observed nothing at returns an EMPTY matrix rather than a
 * fabricated law: an analytical fallback carrying the right mean says nothing
 * about the tail, which is the whole reason to ask a simulator for a CDF.
 * Repeated observations are collapsed, keeping the largest CDF value at each
 * distinct time, or the ecdf is multivalued and an interpolating consumer
 * reads a quantile off whichever duplicate it happens to hit.
 */
inline Matrix<double> ldes_cdf_respt(const LdesResult& r, std::size_t station,
                                     std::size_t job_class) {
    if (station >= r.respTimeSamples.size() ||
        job_class >= r.respTimeSamples[station].size())
        return Matrix<double>();
    std::vector<double> x = r.respTimeSamples[station][job_class];
    if (x.empty()) return Matrix<double>();
    std::sort(x.begin(), x.end());
    const std::size_t n = x.size();
    std::vector<double> tu, fu;
    for (std::size_t i = 0; i < n; ++i) {
        const double f = static_cast<double>(i + 1) / static_cast<double>(n);
        if (i + 1 < n && x[i + 1] == x[i]) continue;  // keep the LAST of a tie
        tu.push_back(x[i]);
        fu.push_back(f);
    }
    Matrix<double> out(tu.size(), 2, 0.0);
    for (std::size_t i = 0; i < tu.size(); ++i) {
        out(i, 0) = fu[i];
        out(i, 1) = tu[i];
    }
    return out;
}
 
/**
 * Port of `getProbAggr`: the marginal probability of a per-class job count at
 * one station.
 *
 * The run is a plain steady-state solve with `--export-histogram` added, which
 * is what `getAvgReward` already does, so the estimate is over the same
 * post-warmup path the mean metrics come from.
 *
 * @param doc     the model.json document, forwarded byte for byte
 * @param o       the run's knobs
 * @param station 1-based station index
 * @param counts  per-class job counts to test at that station
 * @param nclasses the model's class count
 */
inline double ldes_prob_aggr(const std::string& doc, const LdesOptions& o, std::size_t station,
                             const std::vector<double>& counts, std::size_t nclasses) {
    std::vector<std::string> flags;
    flags.push_back("--export-histogram");
    const LdesResult r = solver_ldes_text(doc, o, flags);
    std::vector<LdesStateQuery> q(1);
    q[0].station = station;
    q[0].counts = counts;
    return ldes_prob_from_histogram(r, nclasses, q);
}
 
/**
 * Port of `getProbSysAggr`: the joint probability of a whole aggregate state.
 *
 * @param target (nstations x nclasses) job counts, EVERY station constrained
 */
inline double ldes_prob_sys_aggr(const std::string& doc, const LdesOptions& o,
                                 const Matrix<double>& target) {
    if (target.rows() == 0 || target.cols() == 0)
        throw InputError("SolverLDES: getProbSysAggr needs an (nstations x nclasses) target state");
    std::vector<std::string> flags;
    flags.push_back("--export-histogram");
    const LdesResult r = solver_ldes_text(doc, o, flags);
    const std::size_t K = target.cols();
    std::vector<LdesStateQuery> q(target.rows());
    for (std::size_t i = 0; i < target.rows(); ++i) {
        q[i].station = i + 1;
        q[i].counts.resize(K);
        for (std::size_t k = 0; k < K; ++k) q[i].counts[k] = target(i, k);
    }
    return ldes_prob_from_histogram(r, K, q);
}
 
/**
 * `ldesParseForkJoin`: the engine's "forkJoin" rows mapped onto the Fork/Join pairs of `sn` by (fork name,
 * join name), in fork-then-join node order. A pair the engine did not report, and a non-finite entry, read
 * zero. EMPTY when the model has no pair or the engine wrote no block.
 */
template <class T>
solvers::FjMetrics<double> ldes_fj_metrics(const qn::NetworkStruct<T>& sn, const LdesResult& r) {
    if (solvers::fj_sorted_pairs(sn).empty() || r.fj_forks.empty()) return solvers::FjMetrics<double>();
    const std::size_t R = sn.nclasses;
    solvers::FjMetrics<double> out = solvers::fj_metrics_zero<double>(sn, R);
    const Matrix<double>* src[3] = {&r.FJQLen, &r.FJRespT, &r.FJTput};
    Matrix<double>* dst[3] = {&out.QLen, &out.RespT, &out.Tput};
    for (std::size_t p = 0; p < out.npairs(); ++p) {
        const std::string& fname = sn.nodes[out.forks[p] - 1].name;
        const std::string& jname = sn.nodes[out.joins[p] - 1].name;
        std::size_t e = r.fj_forks.size();
        for (std::size_t k = 0; k < r.fj_forks.size() && k < r.fj_joins.size(); ++k)
            if (r.fj_forks[k] == fname && r.fj_joins[k] == jname) {
                e = k;
                break;
            }
        if (e == r.fj_forks.size()) continue;
        for (int m = 0; m < 3; ++m)
            for (std::size_t c = 0; c < R; ++c) {
                std::size_t cc = c;  // the engine's column of class c, matched by name when it reports names
                for (std::size_t k = 0; k < r.class_names.size(); ++k)
                    if (r.class_names[k] == sn.classes[c].name) {
                        cc = k;
                        break;
                    }
                const double v =
                    (e < src[m]->rows() && cc < src[m]->cols()) ? (*src[m])(e, cc) : 0.0;
                (*dst[m])(p, c) = std::isfinite(v) ? v : 0.0;
            }
    }
    return out;
}
 
/**
 * The same, for a model built through the C++ API.
 *
 * The struct is serialized with `io::network_json_envelope`, the writer whose
 * output the reference readers consume, and the document then takes the same
 * path as any other. WHAT THE STRUCT CANNOT CARRY DOES NOT CROSS: a reward built
 * from a lambda and any construct outside the writer's schema are dropped by the
 * writer, with its own warning, which is why a caller holding a `model.json`
 * should pass the DOCUMENT rather than a struct parsed from it -- the round trip
 * can only lose.
 *
 * An OI/PAS station with no finite buffer is refused here, by name, before the
 * document is written: its mu(c) crosses as a table and the engine needs a
 * finite buffer to bound the state (`io::require_finite_pas_buffers`, the
 * Python `save_model` refusal). The in-process engine, which holds mu itself,
 * is not gated.
 */
/**
 * Closed-chain population conservation check on a steady-state QN, ported from MATLAB's parseSteady.
 *
 * A closed chain must conserve its population in the reported queue lengths. Engine builds that drop
 * blocked closed jobs at a finite-capacity station (pre-guard ldes.jar) leak population silently, and
 * the reported means then describe a different model. FCR models are excluded (a WAITQ waiter is
 * parked outside the station rows), and so are fork-join models (a forked job is counted once per
 * sibling task, so the rows sum to more than the population by construction). A chain is flagged when
 * |sum(QN) - N| > max(0.05*N, 0.05), NaN entries being omitted from the sum.
 *
 * @return one message per failing chain, empty when none fails
 */
template <class T>
std::vector<std::string> ldes_closed_population_warnings(const qn::NetworkStruct<T>& sn,
                                                         const Matrix<double>& QN) {
    std::vector<std::string> out;
    if (QN.empty() || !sn.regions.empty()) return out;
    for (std::size_t i = 0; i < sn.nodes.size(); ++i)
        if (sn.nodes[i].nodetype == qn::NodeType::Fork) return out;
    const std::vector<double> njobs = sn.njobs();
    for (std::size_t nc = 0; nc < sn.nchains && nc < sn.chains.size(); ++nc) {
        double njobs_chain = 0.0;
        std::vector<std::size_t> inchain;
        for (std::size_t r = 0; r < sn.chains[nc].size() && r < njobs.size(); ++r)
            if (sn.chains[nc][r]) {
                inchain.push_back(r);
                njobs_chain += njobs[r];
            }
        if (!(std::isfinite(njobs_chain) && njobs_chain > 0.0)) continue;
        double qtot = 0.0;
        for (std::size_t i = 0; i < QN.rows(); ++i)
            for (std::size_t k = 0; k < inchain.size(); ++k)
                if (inchain[k] < QN.cols() && !std::isnan(QN(i, inchain[k]))) qtot += QN(i, inchain[k]);
        if (std::fabs(qtot - njobs_chain) > std::max(0.05 * njobs_chain, 0.05)) {
            char buf[512];
            std::snprintf(buf, sizeof(buf),
                          "LDES returned %.3f jobs for closed chain %zu, whose population is %g: the engine "
                          "dropped blocked jobs at a finite-capacity station, so these averages describe a "
                          "different model. Use SolverCTMC/SolverSSA for capped closed networks, or rebuild "
                          "common/ldes.jar.",
                          qtot, nc + 1, njobs_chain);
            out.push_back(buf);
        }
    }
    return out;
}
 
template <class T>
LdesResult solver_ldes(const qn::NetworkStruct<T>& sn, const LdesOptions& o,
                       const std::vector<std::string>& extra_flags =
                           std::vector<std::string>()) {
    io::require_finite_pas_buffers(sn);
    LdesResult r = solver_ldes_text(io::network_json_envelope(sn).dump(), o, extra_flags);
    // Steady state only, as MATLAB's parseSteady runs it; SILENT mutes it as line_warning does.
    if (!o.has_timespan && !r.timed_out &&
        line::util::LineConsole::get_verbose() != line::util::VerboseLevel::SILENT) {
        const std::vector<std::string> w = ldes_closed_population_warnings(sn, r.QN);
        for (std::size_t k = 0; k < w.size(); ++k)
            std::fprintf(stderr, "[SolverLDES] Warning: %s\n", w[k].c_str());
    }
    return r;
}
 
}  // namespace ldes
}  // namespace line
 
#endif  // LINE_SOLVERS_WRAPPERS_LDES_SOLVER_LDES_H