solver_jmt.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_WRAPPERS_JMT_SOLVER_JMT_H
#define LINE_SOLVERS_WRAPPERS_JMT_SOLVER_JMT_H
 
/**
 * @file
 * @ingroup line_solvers
 * Port of `SolverJMT`, the Java Modelling Tools client.
 *
 * WHAT A JMT SOLVE IS. The model is written out as a JMT document -- a
 * `.jsimg` for the discrete-event engine, a `.jmva` for the analytical one --
 * `jmt.commandline.Jmt` is run on it, and the result document JMT leaves beside
 * the model is parsed back into the same `AvgResult` every other solver in this
 * port returns. The translation is `io/jmt_writer.h` and the SHARED
 * `io/jmva_writer.h`, which SolverLQNS writes through as well; this file is the
 * dispatch, the parse and the metric mapping.
 *
 * THREE WAYS TO REACH JMT, in the order `jmtRun.m` tries them:
 *   1. `options.rest_url`, a JMT REST server (the imperialqore/jmt-rest image).
 *      Nothing runs locally.
 *   2. a local JVM plus `common/JMT.jar`, the default.
 *   3. no JVM, but a usable Docker daemon: the same image, run as a container.
 * Every one of them leaves the result at `<model>-result.jsim` or
 * `<model>-result.jmva`, so the parsers do not know which one ran. That
 * contract is what makes the three interchangeable, and it is why the Docker
 * arm copies the container's output back beside the model.
 *
 * CONSENT IS NOT ASSUMED, FOR EITHER FETCH. The Docker arm is reached only when
 * `LINE_JMT_DOCKER` opts in, and the jar is downloaded only when
 * `LINE_JMT_DOWNLOAD` does: this port has no terminal to ask at -- the reference
 * prompts, and refuses in a `-batch` session -- so an absent variable is a
 * refusal, exactly as an empty answer is there. Fetching 31 MiB because a solver
 * was called is not a decision a library may take on its own.
 *
 * WHAT THE FETCH ADDS OVER THE REFERENCE, once consented to: the transfer is
 * https and stays https across redirects, the bytes are checked against a pinned
 * SHA-256 (`JMT_JAR_SHA256`) before anything is installed, and the file is
 * renamed into place only after it verifies. MATLAB, the JAR and python all
 * trust the transport and stop there, so the C++ edition is the only one that
 * authenticates the jar itself; `LINE_JMT_SHA256` takes a different digest
 * deliberately and `LINE_JMT_URL` a mirror.
 *
 * WHAT IS NOT PORTED, and why it is refused rather than approximated:
 *   getProbAggr / getProbSysAggr   Both weigh the simulated trajectory against
 *                                  `sn.state{isf}`, the model's CURRENT state,
 *                                  which `qn::NetworkStruct` does not carry --
 *                                  the same gap that stops `SolverLDES`'s
 *                                  `getProb` in this port.
 *   the transient ensemble         `default` over a finite timespan averages
 *                                  `sampleSysAggr` over `config.replications`
 *                                  seeds, and `sampleSysAggr` reads the JMT
 *                                  arrival/departure LOG files, which requires
 *                                  a Logger on every station of the model; it is
 *                                  composed from `jmt_sample_sys_aggr` in
 *                                  `jmt_logs.h` (`jmt_transient_replications`).
 */
 
#include "line/util/line_console.h"
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <limits>
#include <map>
#include <sstream>
#include <string>
#include <vector>
 
#include <sys/stat.h>
#include <unistd.h>
 
#include "line/api/sn/sn_predicates.h"
#include "line/io/docker_image.h"
#include "line/io/jmva_writer.h"
#include "line/io/jmt_writer.h"
#include "line/lang/qn/network_struct.h"
#include "line/lang/qn/solver_feature_sets.h"
#include "line/solvers/mva/sn_chain.h"
#include "line/solvers/mva/solver_mva_runner.h"
#include "line/util/error.h"
#include "line/util/http.h"
#include "line/util/sha256.h"
#include "line/util/subprocess.h"
#include "line/util/tempdir.h"
#include "line/util/xml.h"
#include "line/solvers/method_registry.h"
 
namespace line {
namespace jmt {
 
/** The default JMT REST/Docker image, MATLAB's `jmtDockerImage` candidate. */
static const char* JMT_DOCKER_IMAGE = "imperialqore/jmt-rest:latest";
 
/** Where a missing jar is fetched from; `LINE_JMT_URL` names a mirror instead. */
static const char* JMT_JAR_URL = "https://line-solver.sourceforge.net/latest/JMT.jar";
 
/**
 * The digest a fetched jar must have, and the reason the fetch can be trusted.
 *
 * The bytes at `JMT_JAR_URL` are a ROLLING BUILD: the path says `latest`, the
 * jar's own `META-INF/maven/jmt/jmt-light/pom.properties` says `version=latest`,
 * and it is republished without a version number changing anywhere. So there is
 * nothing to pin BY NAME, and this is the digest of the build the port was
 * tested against (published 2026-07-30, 32772660 bytes).
 *
 * WHEN UPSTREAM REPUBLISHES, THE DOWNLOAD FAILS LOUDLY and names both digests.
 * That is the intended failure: the alternative is a solver that silently picks
 * up whatever a `latest` URL happens to serve. Set `LINE_JMT_SHA256` to the new
 * digest to take the new build deliberately, and update this constant when the
 * suite has been run against it.
 */
static const char* JMT_JAR_SHA256 =
    "932cf6b34b04f9ba5b7c3ba0b07dc6c9fa80aa816ffcad7a39ac5d20a31bb9cd";
 
/** The options of one JMT solve, `SolverOptions('JMT')` restricted to what is read. */
struct JmtOptions {
    std::string method = "default";  ///< default | jsim | jmva | `jmva.<alg>`
    double samples = 10000.0;        ///< samples per measure; raised to 5000 below
    long seed = 23000;
    bool keep = false;               ///< keep the scratch directory after the solve
    double confint = 0.99;           ///< 0 disables the confidence interval columns
    double max_simulated_time = std::numeric_limits<double>::infinity();
    double max_events = -1.0;        ///< JSIM event budget, -1 unlimited; set by the samplers as MATLAB's maxEvents
    std::string rest_url;            ///< a JMT REST server; empty selects the local JVM
    std::string container;           ///< `options.config.container`, an image override
    int timeout = 3600;              ///< seconds, for the REST and subprocess arms
    bool verbose = false;
    /**
     * `options.config.replications`: the independent JSIM runs the transient ensemble of `default`
     * averages over a finite timespan, seeded `seed`, `seed+1`, ... No effect on steady-state JSIM.
     */
    int replications = 10;
};
 
namespace detail {
 
inline bool is_file(const std::string& p) {
    struct stat st;
    return !p.empty() && ::stat(p.c_str(), &st) == 0 && S_ISREG(st.st_mode);
}
 
inline bool is_exec(const std::string& p) {
    struct stat st;
    return !p.empty() && ::stat(p.c_str(), &st) == 0 && ::access(p.c_str(), X_OK) == 0;
}
 
inline bool is_dir(const std::string& p) {
    struct stat st;
    return !p.empty() && ::stat(p.c_str(), &st) == 0 && S_ISDIR(st.st_mode);
}
 
inline std::string env_or_empty(const char* name) {
    const char* v = std::getenv(name);
    return v != nullptr ? std::string(v) : std::string();
}
 
inline std::string exe_dir() {
    char buf[4096];
    const ssize_t n = ::readlink("/proc/self/exe", buf, sizeof(buf) - 1);
    if (n <= 0) return std::string();
    std::string p(buf, static_cast<std::size_t>(n));
    const std::size_t s = p.find_last_of('/');
    return s == std::string::npos ? std::string() : p.substr(0, s);
}
 
inline std::string cwd() {
    char buf[4096];
    return ::getcwd(buf, sizeof(buf)) != nullptr ? std::string(buf) : std::string();
}
 
/**
 * The `-Djava.io.tmpdir` switch that keeps one run's JMT scratch private.
 *
 * NOT A TIDINESS MEASURE: without it two JMT runs on one host CORRUPT EACH
 * OTHER. `jmt.commandline.Jmt` hands the engine a copy of the model under
 * `java.io.tmpdir`, and `XMLSimulationOutput.writeAllMeasures` then writes the
 * measures to `<that dir>/res_sim_<the sim's name attribute>`, copies that file
 * to `<model>-result.jsim` and DELETES it. Every LINE-written model carries the
 * same name attribute -- `model.jsim`, from `jmt_write_jsim` and its
 * MATLAB/JAR/python twins -- so with a shared /tmp every concurrent run fights
 * over one path, `/tmp/res_sim_model.jsim`. Two outcomes, both seen in the
 * 2026-09-11 suite: the other run deletes the file first and this one dies
 * `FileNotFoundException ... at jmt.commandline.Jmt.copyFile`, or it copies the
 * OTHER run's measures and reports another model's numbers as its own. The
 * model already sits in a `util::TempDir` of its own, so pointing the JVM's
 * temp dir at it makes every scratch file this run creates private to it.
 */
inline std::string private_tmpdir_flag(const std::string& model_path) {
    const std::size_t s = model_path.find_last_of('/');
    const std::string dir = s == std::string::npos ? cwd() : model_path.substr(0, s);
    return "-Djava.io.tmpdir=" + (dir.empty() ? cwd() : dir);
}
 
/**
 * Port of `jmtGetPath`'s SEARCH; the download is `jmt_download_jar` below.
 *
 * The reference fetches JMT.jar from SourceForge the moment it is missing. This
 * function never does: it reports what is already on disk and returns empty
 * otherwise, so a caller that only wants to know where the jar is -- the CLI's
 * environment report, a test deciding whether to skip -- can ask without the
 * question installing anything. `jmt_run` is what escalates from the empty
 * answer to a fetch, and only with consent.
 *
 * `$LINE_JMT_DIR` overrides the search, for a jar kept outside the tree.
 */
inline std::string jmt_jar_path() {
    const std::string forced = env_or_empty("LINE_JMT_DIR");
    if (!forced.empty()) return is_file(forced + "/JMT.jar") ? forced : std::string();
    const std::string roots[2] = {exe_dir(), cwd()};
    for (int r = 0; r < 2; ++r) {
        if (roots[r].empty()) continue;
        std::string dir = roots[r];
        for (int k = 0; k < 7; ++k) {
            if (is_file(dir + "/JMT.jar")) return dir;
            if (is_file(dir + "/common/JMT.jar")) return dir + "/common";
            dir += "/..";
        }
    }
    return std::string();
}
 
/** The first executable `name` on PATH; empty when there is none. */
inline std::string which_on_path(const std::string& name) {
    const std::string path = env_or_empty("PATH");
    std::size_t b = 0;
    while (b <= path.size()) {
        const std::size_t e = path.find(':', b);
        const std::string dir =
            path.substr(b, e == std::string::npos ? std::string::npos : e - b);
        if (!dir.empty() && is_exec(dir + "/" + name)) return dir + "/" + name;
        if (e == std::string::npos) break;
        b = e + 1;
    }
    return std::string();
}
 
/** `$LINE_JAVA`, then `$JAVA_HOME/bin/java`, then PATH; empty when absent. */
inline std::string find_java() {
    const std::string forced = env_or_empty("LINE_JAVA");
    if (!forced.empty()) return is_exec(forced) ? forced : std::string();
    const std::string home = env_or_empty("JAVA_HOME");
    if (!home.empty() && is_exec(home + "/bin/java")) return home + "/bin/java";
    return which_on_path("java");
}
 
/** `LINE_JMT_DOCKER` in the affirmative; anything else, including unset, is no. */
inline bool docker_consented() {
    const std::string v = env_or_empty("LINE_JMT_DOCKER");
    return v == "1" || v == "true" || v == "TRUE" || v == "yes" || v == "y" || v == "Y";
}
 
/**
 * `LINE_JMT_DOWNLOAD` in the affirmative; anything else, including unset, is no.
 *
 * SAME VARIABLE AS THE OTHER THREE EDITIONS, OPPOSITE DEFAULT. There the jar is
 * fetched unless `LINE_JMT_DOWNLOAD=0` refuses it; here it is fetched only when
 * something affirms it. So `=1` and `=0` mean the same thing everywhere and only
 * the unset case differs, which is the policy difference itself: MATLAB, the JAR
 * and python prompt at a session and this port has nobody to ask.
 */
inline bool download_consented() {
    const std::string v = env_or_empty("LINE_JMT_DOWNLOAD");
    return v == "1" || v == "true" || v == "TRUE" || v == "yes" || v == "y" || v == "Y";
}
 
/**
 * Where a downloaded jar is written: `$LINE_JMT_DIR`, else the tree's `common/`,
 * else the working directory.
 *
 * The walk is `jmt_jar_path`'s, looking for the DIRECTORY the jar belongs in
 * rather than for the jar, so the download lands where the next process's search
 * will find it. `$LINE_JMT_DIR` wins because it is where the caller has said the
 * jar lives, and a download that ignored it would be invisible to the search
 * that runs next.
 */
inline std::string jmt_download_dir() {
    const std::string forced = env_or_empty("LINE_JMT_DIR");
    if (!forced.empty()) return forced;
    const std::string roots[2] = {exe_dir(), cwd()};
    for (int r = 0; r < 2; ++r) {
        if (roots[r].empty()) continue;
        std::string dir = roots[r];
        for (int k = 0; k < 7; ++k) {
            if (is_dir(dir + "/common")) return dir + "/common";
            dir += "/..";
        }
    }
    return cwd();
}
 
/**
 * Fetches `url` to `dest` with curl, or wget, whichever is on PATH.
 *
 * NEITHER IS A DEPENDENCY OF THE BUILD: this arm runs only when the caller has
 * opted in, and a host with neither tool simply does not get the download. That
 * is the whole reason the transfer is delegated rather than written here --
 * libcurl would be a link-time obligation on every binary in the tree, and the
 * port's own `util/http.h` has no TLS and refuses `https://` by name, which is
 * the only scheme SourceForge serves the jar on (plain HTTP answers 301).
 *
 * Both are told to stay on https ACROSS REDIRECTS, not merely on the first
 * request: the published URL is a redirect to begin with, and a downgrade to
 * plaintext part-way through is exactly what must not pass unnoticed. Both
 * inherit stdio, so the caller sees the transfer's own progress, as the Docker
 * arm lets `docker pull` report its own.
 */
inline std::string find_fetcher() {
    const std::string curl = which_on_path("curl");
    return curl.empty() ? which_on_path("wget") : curl;
}
 
inline bool fetch_url(const std::string& tool, const std::string& url, const std::string& dest) {
    const bool is_curl = tool.size() >= 4 && tool.compare(tool.size() - 4, 4, "curl") == 0;
    if (is_curl) {
        std::vector<std::string> argv;
        argv.push_back(tool);
        argv.push_back("-fL");
        argv.push_back("--proto");
        argv.push_back("=https");
        argv.push_back("--proto-redir");
        argv.push_back("=https");
        argv.push_back("--tlsv1.2");
        argv.push_back("--connect-timeout");
        argv.push_back("30");
        argv.push_back("--max-time");
        argv.push_back("1800");
        argv.push_back("-o");
        argv.push_back(dest);
        argv.push_back(url);
        return util::run_inherit(argv) == 0;
    }
    if (tool.empty()) return false;
    std::vector<std::string> argv;
    argv.push_back(tool);
    argv.push_back("--https-only");
    argv.push_back("--timeout=30");
    argv.push_back("--tries=2");
    argv.push_back("-O");
    argv.push_back(dest);
    argv.push_back(url);
    return util::run_inherit(argv) == 0;
}
 
/** The digest a fetched jar must match: `$LINE_JMT_SHA256`, else the pin. */
inline std::string expected_jar_sha256() {
    std::string want = util::trim(env_or_empty("LINE_JMT_SHA256"));
    if (want.empty()) want = JMT_JAR_SHA256;
    for (std::size_t i = 0; i < want.size(); ++i)
        want[i] = static_cast<char>(std::tolower(static_cast<unsigned char>(want[i])));
    return want;
}
 
/**
 * The download `jmtGetPath` does, with the consent and the verification the
 * reference has no way to offer.
 *
 * @return the directory the jar now sits in, or EMPTY when nothing was fetched
 *         because nothing consented or no fetcher exists. A fetch that was
 *         ATTEMPTED and went wrong throws instead: a half-installed or
 *         unrecognised jar must not read as "JMT is simply not available here".
 *
 * The bytes land on a `.part-<pid>` beside the destination and are renamed in
 * only once the digest matches, so a failed or truncated transfer never leaves
 * a `JMT.jar` behind for the next run to find and trust, and two processes
 * fetching at once cannot interleave into one file. The digest is the point of
 * the exercise: curl and wget authenticate the SERVER, and only this
 * authenticates the ARTEFACT.
 */
inline std::string jmt_download_jar() {
    if (!download_consented()) return std::string();
 
    std::string url = env_or_empty("LINE_JMT_URL");
    if (url.empty()) url = JMT_JAR_URL;
    if (url.compare(0, 8, "https://") != 0)
        throw InputError("SolverJMT: refusing to fetch JMT.jar over '" + url +
                         "': the download must be https, so that the transfer is authenticated "
                         "before the digest authenticates the file");
 
    const std::string dir = jmt_download_dir();
    if (!is_dir(dir) || ::access(dir.c_str(), W_OK) != 0)
        throw InputError("SolverJMT: cannot write JMT.jar to '" + dir +
                         "': the directory does not exist or is not writable. Set LINE_JMT_DIR to "
                         "a directory that is");
 
    const std::string tool = find_fetcher();
    if (tool.empty())
        throw InputError("SolverJMT: LINE_JMT_DOWNLOAD is set but neither curl nor wget is on "
                         "PATH, and this port has no TLS client of its own. Install one, or place "
                         "JMT.jar in " +
                         dir + " by hand");
 
    const std::string dest = dir + "/JMT.jar";
    const std::string part = dest + ".part-" + std::to_string(static_cast<long>(::getpid()));
    ::unlink(part.c_str());
    std::printf("\nJava Modelling Tools cannot be found. Fetching %s (about 31 MiB) into %s\n",
                url.c_str(), dir.c_str());
    std::fflush(stdout);
    if (!fetch_url(tool, url, part)) {
        ::unlink(part.c_str());
        throw InputError("SolverJMT: the download of " + url + " failed; see the transfer's own "
                         "message above");
    }
 
    const std::string want = expected_jar_sha256();
    const std::string got = util::sha256_file_hex(part);
    if (got != want) {
        struct stat st;
        const long long bytes = ::stat(part.c_str(), &st) == 0 ? static_cast<long long>(st.st_size)
                                                               : -1;
        ::unlink(part.c_str());
        throw InputError(
            "SolverJMT: the jar fetched from " + url + " is not the one this port pins.\n" +
            "  expected sha256 " + want + "\n" + "  received sha256 " +
            (got.empty() ? std::string("<unreadable>") : got) + " (" + std::to_string(bytes) +
            " bytes)\n" +
            "The file has been discarded. That URL serves a rolling 'latest' build, so an upstream "
            "republication looks exactly like this: check the digest against the publisher and, if "
            "it is the build you want, set LINE_JMT_SHA256 to it");
    }
 
    if (::rename(part.c_str(), dest.c_str()) != 0) {
        ::unlink(part.c_str());
        throw InputError("SolverJMT: could not move the verified jar into place at '" + dest + "'");
    }
    std::printf("Verified sha256 %s and installed %s\n", got.c_str(), dest.c_str());
    std::fflush(stdout);
    return dir;
}
 
inline std::string read_file(const std::string& path) {
    std::ifstream f(path.c_str(), std::ios::binary);
    if (!f) throw InputError("SolverJMT: cannot read '" + path + "'");
    std::ostringstream ss;
    ss << f.rdbuf();
    return ss.str();
}
 
inline void write_file(const std::string& path, const std::string& text) {
    std::ofstream f(path.c_str(), std::ios::binary);
    if (!f) throw InputError("SolverJMT: cannot write '" + path + "'");
    f.write(text.data(), static_cast<std::streamsize>(text.size()));
    if (!f) throw InputError("SolverJMT: write failed on '" + path + "'");
}
 
/** JSON string escaping, for the REST request body. */
inline std::string json_escape(const std::string& s) {
    std::string out;
    out.reserve(s.size() + 16);
    for (std::size_t i = 0; i < s.size(); ++i) {
        const unsigned char c = static_cast<unsigned char>(s[i]);
        switch (c) {
            case '"': out += "\\\""; break;
            case '\\': out += "\\\\"; break;
            case '\n': out += "\\n"; break;
            case '\r': out += "\\r"; break;
            case '\t': out += "\\t"; break;
            default:
                if (c < 0x20) {
                    char buf[8];
                    std::snprintf(buf, sizeof(buf), "\\u%04x", c);
                    out += buf;
                } else {
                    out.push_back(s[i]);
                }
        }
    }
    return out;
}
 
/**
 * The value of a top-level JSON string field, unescaped.
 *
 * The REST response is a small, flat envelope -- `status`, `error`,
 * `raw_output.result_xml`, `raw_output.stdout` -- and pulling four strings out
 * of it does not warrant making this header depend on the JSON library that
 * only the LDES client uses. A malformed response yields an empty string, and
 * the caller then reports that the server sent no result document.
 */
inline std::string json_string_field(const std::string& body, const std::string& key) {
    const std::string needle = "\"" + key + "\"";
    std::size_t p = body.find(needle);
    if (p == std::string::npos) return std::string();
    p = body.find(':', p + needle.size());
    if (p == std::string::npos) return std::string();
    ++p;
    while (p < body.size() && (body[p] == ' ' || body[p] == '\t' || body[p] == '\n')) ++p;
    if (p >= body.size() || body[p] != '"') return std::string();
    ++p;
    std::string out;
    while (p < body.size() && body[p] != '"') {
        if (body[p] == '\\' && p + 1 < body.size()) {
            ++p;
            switch (body[p]) {
                case 'n': out.push_back('\n'); break;
                case 'r': out.push_back('\r'); break;
                case 't': out.push_back('\t'); break;
                case 'u': {
                    if (p + 4 < body.size()) {
                        const long cp = std::strtol(body.substr(p + 1, 4).c_str(), nullptr, 16);
                        if (cp < 0x80) out.push_back(static_cast<char>(cp));
                        p += 4;
                    }
                    break;
                }
                default: out.push_back(body[p]);
            }
        } else {
            out.push_back(body[p]);
        }
        ++p;
    }
    return out;
}
 
}  // namespace detail
 
template <class T>
struct JmtResult;
 
inline util::ProcResult jmt_run_docker(const std::string& image, const std::string& mode,
                                       const std::string& model_path, long seed,
                                       const JmtOptions& opt);
 
/** The suffix `jmt.commandline.Jmt` appends to the model path, per mode. */
inline const char* jmt_result_ext(const std::string& mode) {
    if (mode == "sim") return "jsim";
    if (mode == "mva") return "jmva";
    throw InputError("SolverJMT: unknown JMT analysis mode '" + mode + "'");
}
 
/** True when a local JVM and `common/JMT.jar` are both present. */
inline bool jmt_available() {
    return !detail::jmt_jar_path().empty() && !detail::find_java().empty();
}
 
/**
 * Port of `jmtSolveRest`: POST the model document, write the result document
 * back beside the model so the local parsers are unaffected.
 */
inline util::ProcResult jmt_solve_rest(const std::string& rest_url, const std::string& mode,
                                       const std::string& model_path, long seed,
                                       const JmtOptions& opt) {
    std::string url = rest_url;
    while (!url.empty() && url[url.size() - 1] == '/') url.erase(url.size() - 1);
    // A caller may hand over the full route; anything else gets the default one.
    if (url.size() < 8 || url.substr(url.size() - 8) != "/solve/" + mode.substr(0, 1))
        if (url.find("/api/v") == std::string::npos) url += "/api/v1/solve/" + mode;
 
    std::string body = "{\"model\":{\"content\":\"" +
                       detail::json_escape(detail::read_file(model_path)) +
                       "\",\"base64\":false}";
    // JMVA takes its algorithm and tolerance from the model document, so the
    // seed is meaningful only on the simulation route; the server's option
    // allow-list rejects it on /solve/mva.
    if (mode == "sim") body += ",\"options\":{\"seed\":" + std::to_string(seed) + "}";
    body += "}";
 
    const http::Response resp = http::post_json(url, body, opt.timeout * 1000);
    if (resp.status != 200)
        throw InputError("SolverJMT: the JMT REST server at " + url + " answered HTTP " +
                         std::to_string(resp.status));
    const std::string status = detail::json_string_field(resp.body, "status");
    if (status != "completed") {
        const std::string err = detail::json_string_field(resp.body, "error");
        throw InputError("SolverJMT: the JMT REST solve failed: " +
                         (err.empty() ? std::string("unspecified error") : err));
    }
    const std::string xml = detail::json_string_field(resp.body, "result_xml");
    if (xml.empty())
        throw InputError(
            "SolverJMT: the JMT REST response carries no result document. The server was asked "
            "to include the raw output; check that include_raw_output is not disabled");
    detail::write_file(model_path + "-result." + jmt_result_ext(mode), xml);
    util::ProcResult r;
    r.exitCode = 0;
    r.out = detail::json_string_field(resp.body, "stdout");
    return r;
}
 
/**
 * Port of `jmtRun`: one batch analysis, leaving the result where the JMT CLI
 * itself would leave it.
 *
 * @param mode "sim" or "mva"
 * @param model_path the document JMT is to read; the result lands beside it
 */
inline util::ProcResult jmt_run(const std::string& mode, const std::string& model_path, long seed,
                                const JmtOptions& opt) {
    if (!opt.rest_url.empty()) return jmt_solve_rest(opt.rest_url, mode, model_path, seed, opt);
 
    std::string jar_dir = detail::jmt_jar_path();
    const std::string java = detail::find_java();
    // A jar is worth fetching only when there is a JVM to run it with: the
    // Docker arm below needs none, so downloading one for it would be 31 MiB
    // spent on a file nothing reads. Returns empty unless LINE_JMT_DOWNLOAD
    // consents; see the header note.
    if (jar_dir.empty() && !java.empty()) jar_dir = detail::jmt_download_jar();
    if (!jar_dir.empty() && !java.empty()) {
        std::vector<std::string> argv;
        argv.push_back(java);
        argv.push_back(detail::private_tmpdir_flag(model_path));
        argv.push_back("-cp");
        argv.push_back(jar_dir + "/JMT.jar");
        argv.push_back("jmt.commandline.Jmt");
        argv.push_back(mode);
        argv.push_back(model_path);
        argv.push_back("-seed");
        argv.push_back(std::to_string(seed));
        argv.push_back("--illegal-access=permit");
        const util::ProcResult r = util::capture(argv, opt.timeout, true);
        // A KILLED RUN IS NOT A RESULT. Without this the caller falls through to
        // "no result document", which reads as a solver that refused the model
        // rather than one whose wall-clock budget expired -- and a jsim that is
        // merely SLOW on a loaded host is exactly the case that hits it. Same
        // rule the LQNS, QNS and LDES wrappers already apply to util::capture.
        if (r.timedOut)
            throw ExternalError(ErrorKind{"JMT"}, "SolverJMT: the JMT run did not finish within " +
                               std::to_string(opt.timeout) + "s and was killed");
        return r;
    }
 
    // The Docker arm. The image is used only when it is ALREADY PRESENT or the
    // caller has opted into the pull; see the header note on consent.
    std::string image = opt.container.empty() ? detail::env_or_empty("LINE_JMT_IMAGE")
                                              : opt.container;
    if (image.empty()) image = JMT_DOCKER_IMAGE;
    if (io::docker_daemon_available()) {
        bool have = io::docker_has_local_image(image);
        if (!have && detail::docker_consented() && io::docker_has_storage_for(image))
            have = io::docker_pull(image);
        if (have) return jmt_run_docker(image, mode, model_path, seed, opt);
    }
 
    // Say which half was missing: a jar outside the searched roots reads like a missing JVM otherwise.
    const std::string found =
        std::string(java.empty() ? "no Java runtime (LINE_JAVA, JAVA_HOME or PATH)" : "a Java runtime") +
        (jar_dir.empty() ? " and no JMT.jar (LINE_JMT_DIR, else JMT.jar or common/JMT.jar in the "
                           "executable's or working directory or up to 6 levels above)"
                         : " and JMT.jar in " + jar_dir);
    throw ExternalError(ErrorKind{"JMT"},
        "SolverJMT: no way to reach JMT: found " + found + ". Install a Java runtime and put "
        "JMT.jar in common/ (" + std::string(JMT_JAR_URL) +
        ") or in the directory LINE_JMT_DIR names, or set LINE_JMT_DOWNLOAD=1 to let curl or wget "
        "fetch and verify that jar for you, or "
        "point options.rest_url at a JMT REST server, or set LINE_JMT_DOCKER=1 to allow the " +
        image + " image to be pulled and used");
}
 
/**
 * Run the analysis inside the JMT container and copy the result back beside the
 * model, so the caller sees the layout a local JVM would have produced.
 *
 * The model is STAGED under a fresh directory and the container is given that
 * directory: in `mva` mode JMT rewrites the model file itself, and a confined
 * Docker cannot bind-mount the system temp directory the model normally lives
 * in. The container runs as the calling user so the files it writes are not
 * left owned by root.
 */
inline util::ProcResult jmt_run_docker(const std::string& image, const std::string& mode,
                                       const std::string& model_path, long seed,
                                       const JmtOptions& opt) {
    util::TempDir work("jmt-docker", true);
    const std::size_t slash = model_path.find_last_of('/');
    const std::string base =
        slash == std::string::npos ? model_path : model_path.substr(slash + 1);
    const std::string staged = work.path() + "/" + base;
    detail::write_file(staged, detail::read_file(model_path));
 
    const std::string uid = std::to_string(static_cast<long>(::getuid()));
    const std::string gid = std::to_string(static_cast<long>(::getgid()));
    std::vector<std::string> argv;
    argv.push_back("docker");
    argv.push_back("run");
    argv.push_back("--rm");
    argv.push_back("--user");
    argv.push_back(uid + ":" + gid);
    argv.push_back("-v");
    argv.push_back(work.path() + ":" + work.path());
    argv.push_back("-w");
    argv.push_back(work.path());
    argv.push_back(image);
    argv.push_back(mode);
    argv.push_back(base);
    argv.push_back("-seed");
    argv.push_back(std::to_string(seed));
    const util::ProcResult r = util::capture(argv, opt.timeout, true);
    if (r.timedOut)
        throw ExternalError(ErrorKind{"JMT"}, "SolverJMT: the JMT container run did not finish within " +
                           std::to_string(opt.timeout) + "s and was killed");
 
    const std::string produced = staged + "-result." + jmt_result_ext(mode);
    if (detail::is_file(produced))
        detail::write_file(model_path + "-result." + jmt_result_ext(mode),
                           detail::read_file(produced));
    return r;
}
 
/** One `<measure>` of a JMT result document, by its attributes. */
struct JmtMeasure {
    std::string measure_type, station, job_class, node_type;
    double mean = 0.0, lower = 0.0, upper = 0.0;
    double analyzed_samples = 0.0;
    bool successful = false;
    bool has_bounds = false;
};
 
/**
 * Port of `getResultsJSIM`: every `<measure>` of the result document.
 *
 * A MISSING FILE IS A FAILED SIMULATION, not an empty result: JMT writes the
 * document even when every measure is unsuccessful, so its absence means the
 * run itself did not complete, and returning zeros would report a solved model.
 */
inline std::vector<JmtMeasure> jmt_parse_measures(const std::string& result_path,
                                                  const std::string& command_output) {
    if (!detail::is_file(result_path)) {
        std::string msg =
            "SolverJMT: JMT did not output a result file, the simulation has likely failed.";
        if (!command_output.empty()) msg += " JMT output: " + command_output;
        throw ExternalError(ErrorKind{"JMT"}, msg);
    }
    const std::unique_ptr<xml::Element> doc = xml::parse_file(result_path);
    std::vector<JmtMeasure> out;
    const std::vector<const xml::Element*> ms = doc->by_tag("measure");
    for (std::size_t i = 0; i < ms.size(); ++i) {
        const xml::Element* e = ms[i];
        JmtMeasure m;
        m.measure_type = e->attr("measureType");
        m.station = e->attr("station");
        m.job_class = e->attr("class");
        m.node_type = e->attr("nodeType");
        m.mean = std::strtod(e->attr("meanValue").c_str(), nullptr);
        m.analyzed_samples = std::strtod(e->attr("analyzedSamples").c_str(), nullptr);
        m.successful = e->attr("successful") == "true";
        if (e->has_attr("lowerLimit") && e->has_attr("upperLimit")) {
            m.lower = std::strtod(e->attr("lowerLimit").c_str(), nullptr);
            m.upper = std::strtod(e->attr("upperLimit").c_str(), nullptr);
            m.has_bounds = true;
        }
        out.push_back(m);
    }
    return out;
}
 
/**
 * The result of a JMT solve: the shared `AvgResult` plus what only JMT reports.
 *
 * The FCR rows extend the metric matrices past `nstations`, which is the
 * reference's layout (`getResults.m` allocates `nstations + nregions` rows) and
 * is what makes a region's measures reachable through the same table as a
 * station's.
 */
template <class T>
struct JmtResult {
    mva::AvgResult<T> avg;
    Matrix<T> QCI, UCI, RCI, TCI, ACI;   ///< confidence half-widths, empty when disabled
    Matrix<T> Weight, MemOcc;            ///< FCR-only: weighted and memory occupation
    Matrix<T> TNfcr, DropRateNfcr;       ///< (nregions x nclasses) carried and lost rate
    /** getResults.m `result.Avg.Tard` (stations x classes) and `.SysTard` (1 x classes). */
    Matrix<T> Tard, SysTard;
    Matrix<T> TardCI, SysTardCI;         ///< their half-widths, empty when disabled
    /** Per Cache node (1-based node index), the per-class hit probability. */
    std::map<std::size_t, std::vector<T>> cache_hit_prob;
    /**
     * `result.Prob.logNormConstAggr`, the log normalizing constant JMVA reports.
     *
     * NaN on the simulation path and on any JMVA algorithm that does not
     * compute one -- an AMVA approximation has no G to report -- rather than
     * zero, which is a legitimate value of a log constant.
     */
    T log_norm_const = num_traits<T>::from_double(std::numeric_limits<double>::quiet_NaN());
};
 
/**
 * Port of `getResults.m`: the measures mapped onto the metric matrices.
 *
 * THE RECURRENCE TEST IS WHY `analyzedSamples` IS READ. For a CLOSED class JMT
 * may terminate a measure before it has seen more samples than there are jobs
 * in the chain, and such a measure is not an estimate of anything; the
 * reference zeroes it rather than reporting it, and so does this. An OPEN class
 * has no such bound and is taken as reported.
 */
template <class T>
JmtResult<T> jmt_map_measures(const qn::NetworkStruct<T>& sn,
                              const std::vector<JmtMeasure>& measures, bool confint) {
    const std::size_t M = sn.nstations, K = sn.nclasses, F = sn.regions.size();
    const T zero = num_traits<T>::from_int(0);
    const double nan = std::numeric_limits<double>::quiet_NaN();
    JmtResult<T> res;
    res.avg.QN = Matrix<T>(M + F, K, zero);
    res.avg.UN = Matrix<T>(M + F, K, zero);
    res.avg.RN = Matrix<T>(M + F, K, zero);
    res.avg.TN = Matrix<T>(M + F, K, zero);
    res.avg.AN = Matrix<T>(M + F, K, zero);
    res.avg.WN = Matrix<T>(M + F, K, zero);
    res.Weight = Matrix<T>(M + F, K, num_traits<T>::from_double(nan));
    res.MemOcc = Matrix<T>(M + F, K, num_traits<T>::from_double(nan));
    if (confint) {
        res.QCI = Matrix<T>(M, K, zero);
        res.UCI = Matrix<T>(M, K, zero);
        res.RCI = Matrix<T>(M, K, zero);
        res.TCI = Matrix<T>(M, K, zero);
        res.ACI = Matrix<T>(M, K, zero);
        res.TardCI = Matrix<T>(M, K, zero);
        res.SysTardCI = Matrix<T>(1, K, zero);
    }
    res.Tard = Matrix<T>(M, K, zero);
    res.SysTard = Matrix<T>(1, K, zero);
    // JMT reports no utilization and no arrival rate for a region.
    for (std::size_t f = 0; f < F; ++f)
        for (std::size_t r = 0; r < K; ++r) {
            res.avg.UN(M + f, r) = num_traits<T>::from_double(nan);
            res.avg.AN(M + f, r) = num_traits<T>::from_double(nan);
        }
 
    std::map<std::string, std::size_t> station_of, class_of;
    for (std::size_t i = 1; i <= sn.nodes.size(); ++i) station_of[sn.nodes[i - 1].name] = i;
    for (std::size_t r = 1; r <= K; ++r) class_of[sn.classes[r - 1].name] = r;
 
    // The chain population each closed class's recurrence test compares against.
    std::vector<double> chainpop(K, 0.0);
    for (std::size_t c = 0; c < sn.inchain.size(); ++c) {
        double tot = 0.0;
        for (std::size_t r : sn.inchain[c]) {
            const double n = sn.classes[r - 1].population;
            if (std::isfinite(n)) tot += n;
        }
        for (std::size_t r : sn.inchain[c]) chainpop[r - 1] = tot;
    }
 
    // Fork-join section measures, per (fork node, class), see JmtWriter::save_fj_metrics. 1-based keys.
    std::map<std::pair<std::size_t, std::size_t>, double> fj_respt, fj_qlen;
 
    for (std::size_t i = 0; i < measures.size(); ++i) {
        const JmtMeasure& m = measures[i];
        const double half = m.has_bounds ? (m.upper - m.lower) / 2.0 : 0.0;
 
        if (m.measure_type == "Fork Join Response Time" ||
            m.measure_type == "Fork Join Number of Customers") {
            const auto ni = station_of.find(m.station);
            const auto ci = class_of.find(m.job_class);
            if (ni != station_of.end() && ci != class_of.end()) {
                const std::pair<std::size_t, std::size_t> key(ni->second, ci->second);
                if (m.measure_type == "Fork Join Response Time")
                    fj_respt[key] = m.mean;
                else
                    fj_qlen[key] = m.mean;
            }
            continue;
        }
 
        // A region measure is identified by its nodeType, or by the name the
        // writer gave it -- JMT does not always echo the nodeType back.
        const bool is_fcr =
            m.node_type == "region" || m.station.compare(0, 8, "FCRegion") == 0;
        if (is_fcr && m.station.compare(0, 8, "FCRegion") == 0) {
            const long f = std::strtol(m.station.substr(8).c_str(), nullptr, 10);
            if (f < 1 || static_cast<std::size_t>(f) > F) continue;
            const std::size_t row = M + static_cast<std::size_t>(f) - 1;
            // A region measure is AGGREGATE, not per class. The reference
            // spreads the extensive ones (customers, throughput, arrival rate,
            // weight, memory) evenly over the classes and repeats the intensive
            // ones (response and residence time) unchanged, so that summing a
            // column recovers the region total in both cases.
            const double per_class = m.mean / static_cast<double>(K);
            for (std::size_t r = 0; r < K; ++r) {
                if (m.measure_type == "Number of Customers")
                    res.avg.QN(row, r) = num_traits<T>::from_double(per_class);
                else if (m.measure_type == "Response Time")
                    res.avg.RN(row, r) = num_traits<T>::from_double(m.mean);
                else if (m.measure_type == "Residence Time")
                    res.avg.WN(row, r) = num_traits<T>::from_double(m.mean);
                else if (m.measure_type == "Throughput")
                    res.avg.TN(row, r) = num_traits<T>::from_double(per_class);
                else if (m.measure_type == "Arrival Rate")
                    res.avg.AN(row, r) = num_traits<T>::from_double(per_class);
                else if (m.measure_type == "FCR Capacity")
                    res.Weight(row, r) = num_traits<T>::from_double(per_class);
                else if (m.measure_type == "FCR Memory")
                    res.MemOcc(row, r) = num_traits<T>::from_double(per_class);
            }
            continue;
        }
 
        if (m.measure_type == "Cache Hit Rate") {
            const auto ni = station_of.find(m.station);
            const auto ci = class_of.find(m.job_class);
            if (ni == station_of.end() || ci == class_of.end()) continue;
            if (sn.nodes[ni->second - 1].nodetype != lang::NodeType::Cache) continue;
            const auto cp = sn.nodeparam.find(ni->second);
            if (cp == sn.nodeparam.end()) continue;
            std::vector<T>& hit = res.cache_hit_prob[ni->second];
            if (hit.empty()) hit.assign(K, zero);
            // The measure names the HIT class; the probability belongs to the
            // read class that switches into it.
            for (std::size_t r = 0; r < cp->second.hitclass.size(); ++r)
                if (cp->second.hitclass[r] == ci->second)
                    hit[r] = num_traits<T>::from_double(m.mean);
            continue;
        }
 
        // A SYSTEM measure names no station (saveMetric.m blanks referenceNode), so it is
        // matched on the class alone, with the same recurrence test as a station measure.
        if (m.measure_type == "System Tardiness") {
            const auto ci = class_of.find(m.job_class);
            if (ci == class_of.end()) continue;
            const std::size_t r = ci->second;
            const bool open = !std::isfinite(sn.classes[r - 1].population);
            const bool recurrent = open || m.analyzed_samples > chainpop[r - 1];
            res.SysTard(0, r - 1) = num_traits<T>::from_double(recurrent ? m.mean : 0.0);
            if (confint) res.SysTardCI(0, r - 1) = num_traits<T>::from_double(recurrent ? half : 0.0);
            continue;
        }
 
        const auto ni = station_of.find(m.station);
        const auto ci = class_of.find(m.job_class);
        if (ni == station_of.end() || ci == class_of.end()) continue;
        const std::size_t ist = sn.nodes[ni->second - 1].station;
        if (ist == 0) continue;
        const std::size_t r = ci->second;
        const bool open = !std::isfinite(sn.classes[r - 1].population);
        const bool recurrent = open || m.analyzed_samples > chainpop[r - 1];
        const T v = num_traits<T>::from_double(recurrent ? m.mean : 0.0);
        const T ci_v = num_traits<T>::from_double(recurrent ? half : 0.0);
 
        if (m.measure_type == "Number of Customers") {
            res.avg.QN(ist - 1, r - 1) = v;
            if (confint) res.QCI(ist - 1, r - 1) = ci_v;
        } else if (m.measure_type == "Utilization") {
            res.avg.UN(ist - 1, r - 1) = v;
            if (confint) res.UCI(ist - 1, r - 1) = ci_v;
        } else if (m.measure_type == "Response Time") {
            res.avg.RN(ist - 1, r - 1) = v;
            if (confint) res.RCI(ist - 1, r - 1) = ci_v;
        } else if (m.measure_type == "Throughput") {
            res.avg.TN(ist - 1, r - 1) = v;
            if (confint) res.TCI(ist - 1, r - 1) = ci_v;
        } else if (m.measure_type == "Arrival Rate") {
            res.avg.AN(ist - 1, r - 1) = v;
            if (confint) res.ACI(ist - 1, r - 1) = ci_v;
        } else if (m.measure_type == "Tardiness") {
            res.Tard(ist - 1, r - 1) = v;
            if (confint) res.TardCI(ist - 1, r - 1) = ci_v;
        }
        // `Residence Time` is deliberately neither requested nor read: JMT's
        // definition disagrees with LINE's on class-switching models, and the
        // residence time is recomputed below from the response time.
    }
    // JMT keeps no utilization at an ORDINARY Place (a storage section with no server). A Place is an
    // INF station, so U == Q, the value SolverCTMC, SolverSSA and SolverNC report. A queueing place
    // keeps the busy-server utilization JMT measures for its embedded queue.
    for (std::size_t i = 0; i < sn.nstations; ++i) {
        if (sn.stations[i].nodetype != lang::NodeType::Place || sn.is_queueing_place(i)) continue;
        for (std::size_t r = 0; r < K; ++r) {
            res.avg.UN(i, r) = res.avg.QN(i, r);
            if (confint) res.UCI(i, r) = res.QCI(i, r);
        }
    }
    // getResults.m's FJ block, carried as `avg.fj`: FJRespT and FJQLen as measured at the Fork (FJQLen is an
    // independent measurement, so Little holds only up to noise), FJTput the Join's class throughput.
    if (!solvers::fj_sorted_pairs(sn).empty()) {
        solvers::FjMetrics<T>& fjm = res.avg.fj;
        fjm = solvers::fj_metrics_zero<T>(sn, K);
        for (std::size_t p = 0; p < fjm.npairs(); ++p) {
            const std::size_t jst = sn.nodes[fjm.joins[p] - 1].station;
            for (std::size_t r = 1; r <= K; ++r) {
                const std::pair<std::size_t, std::size_t> key(fjm.forks[p], r);
                const auto rt = fj_respt.find(key);
                if (rt == fj_respt.end() || !std::isfinite(rt->second)) continue;
                fjm.RespT(p, r - 1) = num_traits<T>::from_double(rt->second);
                const auto ql = fj_qlen.find(key);
                if (ql != fj_qlen.end() && std::isfinite(ql->second))
                    fjm.QLen(p, r - 1) = num_traits<T>::from_double(ql->second);
                if (jst != 0 && jst <= M) fjm.Tput(p, r - 1) = res.avg.TN(jst - 1, r - 1);
            }
        }
    }
    return res;
}
 
/**
 * The region loss table, port of the `sn.nregions > 0` tail of `getResults.m`.
 *
 * JMT EXPOSES NO REGION DROP MEASURE, and its region throughput is the CARRIED
 * (admitted) rate, so the offered rate is reconstructed by flow balance: the
 * rate at which stations outside the region route jobs across its boundary,
 * which the drop does not affect. The difference, clamped at zero, is the loss;
 * a region whose rule is not DROP loses nothing by construction and is zeroed.
 */
template <class T>
void jmt_region_losses(const qn::NetworkStruct<T>& sn, JmtResult<T>& res) {
    const std::size_t M = sn.nstations, K = sn.nclasses, F = sn.regions.size();
    if (F == 0) return;
    const T zero = num_traits<T>::from_int(0);
    res.TNfcr = Matrix<T>(F, K, zero);
    res.DropRateNfcr = Matrix<T>(F, K, zero);
 
    for (std::size_t f = 0; f < F; ++f) {
        const typename qn::NetworkStruct<T>::Region& rg = sn.regions[f];
        for (std::size_t r = 1; r <= K; ++r) {
            double acc = 0.0;
            // JMT's region Throughput is one value for all classes and its even split gave every class
            // the same carried rate; the per-class carried rate is the class-r flow LEAVING the region
            // (departures to the Sink are routed to the Source in sn.rt).
            double carried = 0.0;
            for (std::size_t ii = 1; ii <= M; ++ii) {
                if (!(ii <= rg.members.size() && rg.members[ii - 1])) continue;
                for (std::size_t j = 1; j <= M; ++j) {
                    if (j <= rg.members.size() && rg.members[j - 1]) continue;
                    for (std::size_t rp = 1; rp <= K; ++rp) {
                        const std::size_t a = (j - 1) * K + (rp - 1);
                        const std::size_t b = (ii - 1) * K + (r - 1);
                        if (sn.rt.rows() <= a || sn.rt.cols() <= b) continue;
                        const double w = num_traits<T>::to_double(sn.rt(a, b));
                        if (w != 0.0)
                            acc += num_traits<T>::to_double(res.avg.TN(j - 1, rp - 1)) * w;
                        const std::size_t ao = (ii - 1) * K + (r - 1), bo = (j - 1) * K + (rp - 1);
                        if (ao < sn.rt.rows() && bo < sn.rt.cols())
                            carried += num_traits<T>::to_double(res.avg.TN(ii - 1, r - 1)) *
                                       num_traits<T>::to_double(sn.rt(ao, bo));
                    }
                }
            }
            res.TNfcr(f, r - 1) = num_traits<T>::from_double(carried);
            const bool drops = r <= rg.rule.size() && rg.rule[r - 1] == lang::DropStrategy::DROP;
            res.DropRateNfcr(f, r - 1) =
                num_traits<T>::from_double(drops ? std::max(0.0, acc - carried) : 0.0);
        }
    }
}
 
/**
 * Port of `getResultsJMVA.m`: the per-CHAIN answer spread back over the classes.
 *
 * JMVA answers per chain, and each measure is converted to a per-class one by
 * the class's share of its chain at that station: `alpha(i,k)` weighted by the
 * ratio of the class service time to the chain's. The conversions are the
 * reference's, measure by measure, and the two that are not a bare share are
 * the reason this cannot be a generic rescale:
 *   Utilization   is divided by the chain's visits at the REFERENCE station and,
 *                 at a multiserver station, scaled by min(N, c)/c.
 *   Residence time is JMVA's per-CHAIN residence and is converted to LINE's
 *                 response time per visit by dividing by the class visits.
 *
 * THE STATION IS RESOLVED BY NAME, not by the position of the `stationresults`
 * block. The reference indexes `sn.nservers`, `ST` and `sn.visits` with the
 * BLOCK index, but the blocks exclude the Source while those tables do not, so
 * on any open model every station's demand is read one row early. Reading the
 * `station` attribute is the same lookup on a closed model and the correct one
 * on an open model.
 */
template <class T>
JmtResult<T> jmt_parse_jmva(const qn::NetworkStruct<T>& sn, const std::string& result_path,
                            const std::string& command_output) {
    if (!detail::is_file(result_path)) {
        std::string msg =
            "SolverJMT: JMT did not output a result file, the analysis has likely failed.";
        if (!command_output.empty()) msg += " JMT output: " + command_output;
        throw ExternalError(ErrorKind{"JMT"}, msg);
    }
    const std::size_t M = sn.nstations, K = sn.nclasses, C = sn.nchains;
    const T zero = num_traits<T>::from_int(0);
    JmtResult<T> res;
    res.avg.QN = Matrix<T>(M, K, zero);
    res.avg.UN = Matrix<T>(M, K, zero);
    res.avg.RN = Matrix<T>(M, K, zero);
    res.avg.TN = Matrix<T>(M, K, zero);
    res.avg.AN = Matrix<T>(M, K, zero);
    res.avg.WN = Matrix<T>(M, K, zero);
 
    const mva::ChainDemands<T> dem = mva::sn_get_demands_chain(sn);
 
    std::map<std::string, std::size_t> station_of;
    for (std::size_t i = 1; i <= M; ++i)
        station_of[sn.nodes[sn.station_to_node[i - 1] - 1].name] = i;
 
    const std::unique_ptr<xml::Element> doc = xml::parse_file(result_path);
    // The normalizing constant, when the algorithm reports one.
    const std::vector<const xml::Element*> nc = doc->by_tag("normconst");
    res.log_norm_const = nc.empty()
                             ? num_traits<T>::from_double(
                                   std::numeric_limits<double>::quiet_NaN())
                             : num_traits<T>::from_double(
                                   std::strtod(nc[0]->attr("logValue").c_str(), nullptr));
 
    // A Source reports its own arrival rate as its throughput and holds no
    // jobs; JMVA does not model it, so the row is filled from the struct.
    if (sn.sourceIdx != 0)
        for (std::size_t r = 0; r < K; ++r)
            if (!sn.disabled[sn.sourceIdx - 1][r])
                res.avg.TN(sn.sourceIdx - 1, r) = sn.rates(sn.sourceIdx - 1, r);
 
    const std::vector<const xml::Element*> blocks = doc->by_tag("stationresults");
    for (std::size_t b = 0; b < blocks.size(); ++b) {
        const auto si = station_of.find(blocks[b]->attr("station"));
        if (si == station_of.end()) continue;
        const std::size_t i = si->second;
        // The divisor is the capacity `write_jmva` exported, max(nservers, max
        // lldscaling): a load-dependent station carries its c in the scaling and
        // leaves nservers at 1, so reading nservers alone reports U = c*E[busy]/c.
        double ns = sn.stations[i - 1].nservers;
        for (const T& s : sn.stations[i - 1].lldscaling)
            ns = std::max(ns, num_traits<T>::to_double(s));
        const std::vector<const xml::Element*> crs = blocks[b]->by_tag("classresults");
        for (std::size_t c = 1; c <= crs.size() && c <= C; ++c) {
            const std::vector<const xml::Element*> ms = crs[c - 1]->by_tag("measure");
            // A multiserver Queue is written as <ldstation servers="1">, and one
            // ldstation switches JMVA to its load-dependent algorithm, whose
            // Utilization is 1-p_i(0) at EVERY station, delay ones included. That
            // is a different random variable from LINE's E[busy servers], so no
            // rescaling recovers it; derive U from the chain throughput instead,
            // which both JMVA algorithms report alike.
            double chain_tput = std::numeric_limits<double>::quiet_NaN();
            for (std::size_t m = 0; m < ms.size(); ++m) {
                if (ms[m]->attr("measureType") == "Throughput") {
                    const std::string traw = ms[m]->attr("meanValue");
                    chain_tput = traw == "NaN" ? std::numeric_limits<double>::quiet_NaN()
                                               : std::strtod(traw.c_str(), nullptr);
                    break;
                }
            }
            for (std::size_t m = 0; m < ms.size(); ++m) {
                const std::string kind = ms[m]->attr("measureType");
                const std::string raw = ms[m]->attr("meanValue");
                const double val = raw == "NaN" ? std::numeric_limits<double>::quiet_NaN()
                                                : std::strtod(raw.c_str(), nullptr);
                const std::size_t refst = dem.refstatchain[c - 1];
                const double vchain_ref =
                    num_traits<T>::to_double(dem.Vchain(refst - 1, c - 1));
                const double stchain = num_traits<T>::to_double(dem.STchain(i - 1, c - 1));
                for (std::size_t k : sn.inchain[c - 1]) {
                    const double stk = num_traits<T>::to_double(dem.ST(i - 1, k - 1));
                    const double al = num_traits<T>::to_double(dem.alpha(i - 1, k - 1));
                    double v = val;
                    if (kind == "Utilization") {
                        if (vchain_ref == 0.0) continue;
                        v = stk * chain_tput / vchain_ref * al;
                        if (std::isfinite(ns)) v /= ns;
                        res.avg.UN(i - 1, k - 1) = num_traits<T>::from_double(v);
                    } else if (kind == "Throughput") {
                        res.avg.TN(i - 1, k - 1) = num_traits<T>::from_double(val * al);
                    } else if (kind == "Number of Customers") {
                        if (stchain == 0.0 || vchain_ref == 0.0) continue;
                        res.avg.QN(i - 1, k - 1) =
                            num_traits<T>::from_double(val * stk / stchain / vchain_ref * al);
                    } else if (kind == "Residence time" || kind == "Residence Time") {
                        if (stchain == 0.0 || vchain_ref == 0.0) continue;
                        double visits = 0.0;
                        if (c - 1 < sn.visits.size()) {
                            const std::size_t sfi = sn.stateful_of_station(i);
                            if (sfi != 0 && sn.visits[c - 1].rows() >= sfi)
                                visits = num_traits<T>::to_double(sn.visits[c - 1](sfi - 1, k - 1));
                        }
                        if (visits == 0.0) continue;
                        v = (val / visits) * stk / stchain / vchain_ref * al;
                        res.avg.RN(i - 1, k - 1) = num_traits<T>::from_double(v);
                    }
                }
            }
        }
    }
    return res;
}
 
/**
 * Port of `SolverJMT.listValidMethods`.
 *
 * THE TRANSIENT ENSEMBLE IS NOT A METHOD. A single sample path is not E[N](t),
 * there being no time-ergodicity at a fixed t, so `default` over a finite
 * timespan averages `sampleSysAggr` over `config.replications` seeds. It is
 * composed from `jmt_sample_sys_aggr` in `jmt_logs.h`
 * (`jmt_transient_replications`) rather than from `solver_jmt_run_analyzer`.
 * The former `replication` method was that ensemble under a method name, and
 * asking for it now raises `jmt_removed_method_refusal`.
 */
inline std::vector<std::string> jmt_list_valid_methods() {
    // the names and their order are this solver's row of the method registry (doc/line-methods.json)
    return ::line::method_registry::methods("JMT");
}
 
/**
 * The migration sentence for a method name SolverJMT no longer has; empty for any other name.
 *
 * An ERROR and not an alias, as for the removed `'console'` in `Solver.m`: `replication` became the
 * `options.config.replications` knob of `default`, and the sentence says so rather than a bare "unknown".
 */
inline std::string jmt_removed_method_refusal(const std::string& method) {
    if (method == "replication")
        return "SolverJMT: the 'replication' method was removed: use a finite timespan with "
               "options.config.replications";
    return std::string();
}
 
/**
 * The structural half of SolverJMT's method gate; empty when admissible.
 *
 * THE RULES A FLAT FEATURE SET CANNOT STATE, and each is named rather than
 * reported as a bare no:
 *   a removed method name   'replication' is refused with its migration
 *                           sentence (`jmt_removed_method_refusal`).
 *   a single-server station the eight closed-form JMVA algorithms are
 *                           single-server only, which is why `write_jmva`
 *                           refuses the model rather than emitting an
 *                           `<ldstation>` the algorithm cannot read. A server
 *                           count has no feature name; the load-dependent half
 *                           of the same restriction DOES, and rides in
 *                           `jmt_feature_set` as an unset LoadDependence.
 *   the load-dependent shape JMT carries a scaling only as a server count, so
 *                           alpha(n) = min(n,c) is written exactly and nothing
 *                           else is.
 *
 * ONE PREDICATE, TWO CALLERS: `solver_jmt_run_analyzer` raises it, so a caller
 * naming the method by hand gets the sentence rather than a JMT stack trace, and
 * `autosolver::auto_family_refusal` returns it, so findSolver never offers the
 * pair. A second copy of any rule is how the gate and the run drift apart.
 */
/**
 * What the JMVA document cannot state; "" otherwise (MATLAB `jmtMethodRefusal`, isJmva arm).
 *
 * The document carries a station type, a mean demand and a visit count per chain and nothing else.
 * A FCFS station whose classes have different mean service times is not BCMP, and the engine answers
 * for the PS station with the same demands (3+2 jobs, rates 2 and 0.5: QLen 1.99/1.76 against the
 * exact 2.42/1.65). Every caller here names a jmt method, so all three rules key on the jmva* name.
 */
template <class T>
std::string jmva_document_refusal(const qn::NetworkStruct<T>& sn, const std::string& method) {
    if (method.compare(0, 4, "jmva") != 0 && method.compare(0, 8, "jmt.jmva") != 0) return std::string();
    if (api::sn_has_multi_class_heter_fcfs(sn))
        return method + " requires a product-form model, but a FCFS station serves its classes with "
               "different mean service times, and the JMVA engine would answer for the "
               "processor-sharing station with the same demands. Use SolverJMT 'jsim', SolverMVA "
               "(whose FCFS correction handles it), SolverCTMC or SolverLDES.";
    if (api::sn_has_fork_join(sn))
        return "The JMVA document has no fork or join element, so the analytical engine would "
               "solve the model without the synchronisation. Use a simulator (SolverJMT 'jsim', "
               "SolverLDES) or SolverMVA, which approximates fork-join.";
    for (std::size_t i = 0; i < sn.nstations; ++i) {
        const qn::SchedStrategy sc = sn.stations[i].sched;
        const char* sname = sc == qn::SchedStrategy::FCFS ? "FCFS"
                          : sc == qn::SchedStrategy::SIRO ? "SIRO"
                          : sc == qn::SchedStrategy::LCFS ? "LCFS" : nullptr;
        if (sname == nullptr) continue;
        for (std::size_t r = 0; r < sn.nclasses; ++r) {
            const double rate = num_traits<T>::to_double(sn.rates(i, r));
            const double scv = num_traits<T>::to_double(sn.scv(i, r));
            if (!(std::isfinite(rate) && rate > 0.0 && std::isfinite(scv))) continue;
            if (std::fabs(scv - 1.0) <= qn::GlobalConstants::FineTol) continue;
            char buf[64];
            std::snprintf(buf, sizeof(buf), "%.4g", scv);
            return "Station " + sn.stations[i].name + " serves class " + sn.classes[r].name + " under " +
                   sname + " with a non-exponential law (SCV " + buf + "). The JMVA document carries "
                   "the mean demand only, so the analytical engine would answer for the exponential "
                   "law. Use a simulator (SolverJMT 'jsim', SolverLDES), SolverMVA, whose AMVA "
                   "corrections read the SCV, or SolverCTMC or SolverSSA.";
        }
    }
    return std::string();
}
 
template <class T>
std::string jmt_method_refusal(const qn::NetworkStruct<T>& sn, const std::string& method,
                               const JmtOptions& opt) {
    (void)opt;
    {
        const std::string removed = jmt_removed_method_refusal(method);
        if (!removed.empty()) return removed;
    }
    {
        const std::string doc = jmva_document_refusal(sn, method);
        if (!doc.empty()) return doc;
    }
    // LINEAR ADMISSION CONSTRAINTS: no JMT region carries A*n <= b, so the writer dropped them.
    for (const auto& rg : sn.regions)
        if (rg.lincon_A.rows() > 0)
            return "SolverJMT does not support linear admission constraints on a finite capacity region (setLinearConstraints): no JMT region carries them, so the region would be simulated without its A*n <= b limit. Use SolverCTMC, SolverNC or SolverLDES.";
    // EXHAUSTIVE AND K-LIMITED POLLING. The export is JMT's own, but the engine
    // reports a queue length that grows without bound next to a sane response
    // time, and a throughput below the arrival rate (QLen 2674, Tput 0.283
    // against the exact 0.883 and 0.3 on two classes at 0.3/s, Exp(1.5)
    // service, Exp(1) switchover). GATED polling is simulated correctly.
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (const lang::PollingType pt : sn.stations[i].polling_type)
            if (pt == lang::PollingType::EXHAUSTIVE || pt == lang::PollingType::KLIMITED)
                return std::string("SolverJMT does not support ") +
                       (pt == lang::PollingType::EXHAUSTIVE ? "EXHAUSTIVE" : "K-LIMITED") +
                       " polling: JMT's polling server for it reports a queue length that grows "
                       "without bound and a throughput below the arrival rate on a stable station. "
                       "Use SolverMVA, SolverLDES or SolverSSA; JMT simulates GATED polling correctly.";
    // IMMEDIATE FEEDBACK, which no JMT engine has semantics for: a self-looping
    // job KEEPS the server it just used, and both the JSIM routing block and the
    // JMVA visit count answer for a job that leaves and re-queues. This lived in
    // `solver_jmt_run_analyzer` alone until 2026-09-13, which is exactly how
    // `--find-solver-all` came to offer `jmt.*` on a model every JMT run throws
    // on.
    if (sn.has_immediate_feedback())
        return "SolverJMT: JMT has no immediate-feedback semantics (a job re-entering service "
               "while HOLDING the server); use SolverCTMC, SolverSSA or SolverLDES";
    if (qn::jmva_is_closed_only(method)) {
        double maxFinite = 0.0;
        for (std::size_t i = 0; i < sn.nstations; ++i) {
            const double c = static_cast<double>(sn.stations[i].nservers);
            if (std::isfinite(c)) maxFinite = std::max(maxFinite, c);
        }
        if (maxFinite > 1.0) return "SolverJMT: " + method + " does not support multi-server stations";
    }
    // JMVA'S MULTICLASS CoMoM ANSWERS A DIFFERENT MODEL EVERY RUN: when its BTF system is singular,
    // SolverMultiClosedCoMoM perturbs every demand with an unseeded MersenneTwister, re-solves, and
    // reports the result as an exact CoMoM solution. Keyed on the chain count, as in MATLAB
    // jmtMethodRefusal, because the singular case cannot be predicted from the demands; a
    // single-chain model runs SolverSingleClosedRECAL and stays exact. An open or mixed model is
    // left to the OpenClass refusal of jmt_feature_set: the closed-only CoMoM never runs on it.
    const std::vector<double> comom_njobs = sn.njobs();
    const bool comom_closed = std::all_of(comom_njobs.begin(), comom_njobs.end(),
                                          [](double n) { return std::isfinite(n); });
    if ((method == "jmva.comom" || method == "jmt.jmva.comom") && sn.nchains > 1 && comom_closed)
        return method + " is refused on a model with " + std::to_string(sn.nchains) +
               " chains: JMT's multiclass CoMoM engine silently re-solves a randomly perturbed model "
               "whenever its exact linear system is singular, and reports that as an exact CoMoM "
               "solution, so the answer is wrong and different on every run. Use 'jmva.recal', which "
               "is exact and reproducible on the same document, 'jmva'/'jmva.mva' for exact MVA, or "
               "SolverNC for a LINE-native normalizing constant.";
    for (std::size_t i = 0; i < sn.nstations; ++i) {
        const std::vector<T>& alpha = sn.stations[i].lldscaling;
        if (alpha.empty()) continue;
        double c = 0.0;
        for (const T& s : alpha) c = std::max(c, num_traits<T>::to_double(s));
        if (c == 1.0) continue;
        bool ok = (c >= 1.0) && (c == std::floor(c));
        for (std::size_t n = 0; ok && n < alpha.size(); ++n)
            ok = std::abs(num_traits<T>::to_double(alpha[n]) -
                          std::min(static_cast<double>(n + 1), c)) <= qn::GlobalConstants::Zero;
        if (!ok)
            return "SolverJMT: station '" + sn.stations[i].name +
                   "' uses a load-dependent scaling that is not the multiserver encoding alpha(n) "
                   "= min(n,c): JMT has no representation for it, since both the JSIM and the JMVA "
                   "writer carry the scaling as a server count, and the model would be solved at "
                   "the nominal service rate. Use SolverCTMC, SolverNC, SolverMVA or SolverSSA, "
                   "which read sn.lldscaling directly";
    }
    // A BINDING FINITE BUFFER, which neither engine can carry -- JSIM because no
    // JMT drop strategy reproduces LINE's blocking, JMVA because its document has
    // no capacity element at all. The verdict is the WRITER's own
    // (`io::JmtWriter::buffer_capacity_refusal`), asked here without letting it
    // raise, so the gate binds exactly where the writer binds.
    //
    // WHICH ENGINE IS ASKED ABOUT: a method that is neither engine's gets no
    // verdict rather than JSIM's rule applied to a run JSIM is not going to make.
    const bool is_jmva = method.compare(0, 4, "jmva") == 0;
    const bool is_jsim = method == "default" || method == "jsim";
    if (is_jmva || is_jsim) {
        const std::string cap_reason = io::jmt_buffer_capacity_refusal(sn, is_jmva);
        if (!cap_reason.empty()) return cap_reason;
    }
    // A CACHE THE JSIM WRITER REFUSES (FIFO/SFIFO on a list of capacity above 1, HLRU/CLIMB/QLRU, a promotion
    // graph JMT cannot carry, a popularity that is neither Zipf nor uniform): the writer's own verdict
    // (`io::JmtWriter::cache_refusal`), so the gate stops offering jmt.jsim where the run would raise
    // (MATLAB jmtCachePromotion, asked by jmtMethodRefusal and SolverJMT.supports).
    // AN EDD OR EDF STATION WHOSE SERVED CLASSES CARRY NO DUE DATE: the JSIM run would abort on it
    // (jmtMethodRefusal.m asks jmtDeadlineRefusal on the jsim engine only; JMVA has no such discipline).
    if (is_jsim) {
        const std::string deadline_reason = io::jmt_deadline_refusal(sn);
        if (!deadline_reason.empty()) return deadline_reason;
    }
    if (is_jsim) {
        const std::string cache_reason = io::jmt_cache_refusal(sn);
        if (!cache_reason.empty()) return cache_reason;
    }
    return std::string();
}
 
/**
 * Port of `@@SolverJMT/runAnalyzer.m`, the `jsim` and `jmva` arms.
 *
 * The sample count is RAISED to 5000 when the caller asks for less, as the
 * reference does: JMT terminates a measure on its own precision target and
 * needs that many samples per measure before it will report one, so a smaller
 * request silently yields unsuccessful measures rather than a faster run.
 */
template <class T>
JmtResult<T> solver_jmt_run_analyzer(const qn::NetworkStruct<T>& sn, const JmtOptions& opt_in) {
    JmtOptions opt = opt_in;
    if (opt.samples < 5000.0) opt.samples = 5000.0;
    {
        const std::string removed = jmt_removed_method_refusal(opt.method);
        if (!removed.empty()) throw UnsupportedError(ErrorKind{"Method"}, removed);
    }
    const std::vector<std::string> valid = jmt_list_valid_methods();
    if (std::find(valid.begin(), valid.end(), opt.method) == valid.end())
        throw UnsupportedError(ErrorKind{"Method"}, "SolverJMT: unknown method '" + opt.method + "'");
    // runAnalyzerChecks, as MATLAB/JAR/Python ask it: without it a closed-only jmva.* algorithm on an open model
    // reached JMT and came back as "[LINE:external:JMT] ... UnsupportedModelException" instead of a refusal.
    qn::feature_gate("SolverJMT", qn::jmt_feature_set(opt.method), sn);
    // The structural half of the gate, which `auto_family_refusal` also asks:
    // the single-server restriction of the closed-form JMVA algorithms and the
    // load-dependent shape JMT can carry. A second copy of either rule is how
    // the report and the run come to disagree.
    {
        // The id follows the sentence, as in the JAR (jmtStructuralRefusalId), MATLAB and Python: a JSIM writer
        // cache refusal keeps the writer's own LINE:unsupported:Cache, any other structural refusal is Structure.
        const std::string refusal = jmt_method_refusal(sn, opt.method, opt);
        if (!refusal.empty())
            throw UnsupportedError(ErrorKind{refusal == io::jmt_cache_refusal(sn) ? "Cache" : "Structure"}, refusal);
    }
 
    const bool is_mva = opt.method.compare(0, 4, "jmva") == 0;
    const std::string mode = is_mva ? "mva" : "sim";
    util::TempDir work(is_mva ? "jmva" : "jsim");
    if (opt.keep) work.keep();
    const std::string model_path = work.path() + "/model." + (is_mva ? "jmva" : "jsim");
 
    if (is_mva) {
        // `io/jmva_writer.h` is the shared port of `writeJMVA`: SolverLQNS reads
        // the same grammar through `qnsolver`, and a second copy of the writer
        // here would be a second place for the chain aggregation to drift.
        io::write_jmva(sn, model_path, opt.method, static_cast<std::size_t>(opt.samples));
    } else {
        io::JmtWriteOptions wopt;
        wopt.file_name = "model";
        wopt.log_path = sn.log_path;
        wopt.seed = opt.seed;
        wopt.max_samples = opt.samples;
        wopt.max_simulated_time = opt.max_simulated_time;
        wopt.max_events = opt.max_events;
        wopt.sim_conf_int = opt.confint > 0.0 ? opt.confint : 0.99;
        line::util::LineConsole::step("writing the JSIM model file");
        detail::write_file(model_path, io::jmt_write_jsim(sn, wopt));
        line::util::LineConsole::substep("model written to %s", model_path.c_str());
    }
 
    line::util::LineConsole::step("running the JMT simulation engine as a subprocess");
    const util::ProcResult run = jmt_run(mode, model_path, opt.seed, opt);
    const std::string result_path = model_path + "-result." + jmt_result_ext(mode);
 
    line::util::LineConsole::step("parsing the JMT result files");
    JmtResult<T> res;
    if (is_mva) {
        res = jmt_parse_jmva(sn, result_path, run.out);
    } else {
        res = jmt_map_measures(sn, jmt_parse_measures(result_path, run.out), opt.confint > 0.0);
        jmt_region_losses(sn, res);
    }
    // ResidT DERIVED on both arms as getAvg does (JSIM disagrees under class switching, JMVA per-chain); else column zero. Region rows keep JMT total.
    {
        const Matrix<T> WNst = mva::sn_get_residt_from_respt(sn, res.avg.RN);
        for (std::size_t i = 0; i < sn.nstations; ++i)
            for (std::size_t r = 0; r < sn.nclasses; ++r) res.avg.WN(i, r) = WNst(i, r);
    }
    // A SOURCE HAS NO ARRIVALS TO ITSELF. JMT reports an "Arrival Rate" measure
    // at every station the model declares, the Source included, where the
    // reference reports 0: `getAvg.m:197-199` builds a zeroMask over the Source
    // stations and applies it to the ArvR column for EVERY solver. This port's
    // analytical arms already leave that entry at zero, so the JMT reader was
    // the one path that carried JMT's own measure through -- 0.79526 against 0
    // on fcr_mm1kdrop, 0.30214 and 0.4979 on fcr_constraints. Region rows are
    // past `nstations` and are left alone, as the residence time above leaves
    // them.
    for (std::size_t i = 0; i < sn.nstations; ++i) {
        if (sn.stations[i].nodetype != qn::NodeType::Source) continue;
        for (std::size_t r = 0; r < sn.nclasses; ++r)
            res.avg.AN(i, r) = num_traits<T>::from_int(0);
    }
    // A CELL WITH NO RESPONSE TIME HOLDS NO JOBS AND BUSIES NO SERVER.
    // `@@NetworkSolver/getAvg` builds `zeroMask = RN < 10*FineTol` and applies
    // it to the queue length and the utilization for EVERY solver; the
    // analytical arms of this port take it through `mva::filter_metric`, and
    // the JMT reader was the one path that did not. JMT keeps sampling a
    // starved class's utilization long after its response time has been
    // discarded as non-recurrent, so the cell arrives as a small positive
    // number where every other codebase reports 0: 0.0010617 against 0 on
    // prio_hol_closed (PSQueue, Class2), on an otherwise IDENTICAL sample path.
    // Region rows past `nstations` are left alone, as the two rules above are.
    {
        for (std::size_t i = 0; i < sn.nstations; ++i)
            for (std::size_t r = 0; r < sn.nclasses; ++r)
                if (num_traits<T>::to_double(res.avg.RN(i, r)) <
                    10.0 * lang::GlobalConstants::FineTol) {
                    res.avg.QN(i, r) = num_traits<T>::from_int(0);
                    res.avg.UN(i, r) = num_traits<T>::from_int(0);
                }
    }
    res.avg.method = opt_in.method;
    res.avg.actualmethod = opt.method;
    res.avg.iter = 1;
    return res;
}
 
}  // namespace jmt
}  // namespace line
 
#endif  // LINE_SOLVERS_WRAPPERS_JMT_SOLVER_JMT_H