solver_facade.cpp

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
 
/**
 * @file
 * The one translation unit that instantiates the solver stack for the facade.
 *
 * `line/solvers/solver.h` declares a `double`-only, non-template solver API so
 * that this file can carry every body: a caller including that header pays for
 * none of the template instantiation below, and `line_mp_api` pays for it once
 * for the CLI, the tests and the examples together.
 *
 * Every entry point forwards to the runner the CLI calls, so a program written
 * against the facade and `line-cli -s <solver>` on the same model produce the
 * same numbers by construction. Nothing is computed here beyond the reshaping
 * the tables need.
 */
 
#include "line/solvers/solver.h"
 
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <set>
#include <sstream>
#include <utility>
 
#include "line/api/sn/sn_fj_nodevisits_mmt.h"
#include "line/api/sn/sn_zero_capacity.h"
#include "line/solvers/auto/auto_methods.h"
#include "line/solvers/auto/solver_auto.h"
#include "line/lang/qn/nodes.h"
#include "line/solvers/ag/ag_dispatch.h"
#include "line/solvers/ba/solver_ba_runner.h"
#include "line/solvers/ctmc/solver_ctmc_analyzer.h"
#include "line/solvers/ctmc/solver_ctmc_cdf.h"
#include "line/solvers/ctmc/solver_ctmc_getters.h"
#include "line/solvers/ctmc/solver_ctmc_mdd_analyzer.h"
#include "line/solvers/ctmc/solver_ctmc_prob.h"
#include "line/solvers/ctmc/solver_ctmc_symbolic.h"
#include "line/solvers/ctmc/solver_ctmc_waitq.h"
#include "line/solvers/fluid/fluid_runner.h"
#include "line/solvers/fluid/solver_fluid.h"
#include "line/solvers/map_env.h"
#include "line/solvers/map_env_stages.h"
#include "line/solvers/method_registry.h"
#include "line/solvers/mam/solver_mam_runner.h"
#include "line/solvers/mva/solver_mva_runner.h"
#include "line/solvers/nc/solver_nc_runner.h"
#include "line/solvers/solver_chain_tables.h"
#include "line/solvers/ssa/ssa_dispatch.h"
#include "line/solvers/wrappers/jmt/jmt_logs.h"
#include "line/solvers/wrappers/jmt/solver_jmt.h"
#include "line/lang/qn/state.h"
#include "line/solvers/wrappers/ldes/ldes_warm_start.h"
#include "line/solvers/wrappers/ldes/solver_ldes.h"
#include "line/solvers/wrappers/lqns/lqns_qnsolver.h"
 
namespace line {
 
// ---------------------------------------------------------------------------
// AvgTable
// ---------------------------------------------------------------------------
 
double AvgTable::get(const std::string& col, const std::string& station,
                     const std::string& jobclass) const {
    for (std::size_t i = 0; i < Station.size(); ++i) {
        if (Station[i] != station || JobClass[i] != jobclass) continue;
        if (col == "QLen") return QLen[i];
        if (col == "Util") return Util[i];
        if (col == "RespT") return RespT[i];
        if (col == "ResidT") return ResidT[i];
        if (col == "ArvR") return ArvR[i];
        if (col == "Tput") return Tput[i];
        const std::vector<double>& ci = column(col);
        return i < ci.size() ? ci[i] : std::numeric_limits<double>::quiet_NaN();
    }
    return std::numeric_limits<double>::quiet_NaN();
}
 
std::vector<double> AvgTable::column(const std::string& col) const {
    if (col == "QLen") return QLen;
    if (col == "Util") return Util;
    if (col == "RespT") return RespT;
    if (col == "ResidT") return ResidT;
    if (col == "ArvR") return ArvR;
    if (col == "Tput") return Tput;
    if (col == "QLenCI") return QLenCI;
    if (col == "UtilCI") return UtilCI;
    if (col == "RespTCI") return RespTCI;
    if (col == "ResidTCI") return ResidTCI;
    if (col == "ArvRCI") return ArvRCI;
    if (col == "TputCI") return TputCI;
    throw InputError("AvgTable::column: unknown column '" + col + "'");
}
 
namespace {
 
template <class Keep>
AvgTable filter_avg_table(const AvgTable& source, Keep keep) {
    AvgTable out;
    out.solver = source.solver;
    out.method = source.method;
    out.iter = source.iter;
    out.has_lognormconst = source.has_lognormconst;
    out.lognormconst = source.lognormconst;
    out.ListCost = source.ListCost;
    out.warning = source.warning;
    out.SysClass = source.SysClass;
    out.SysRespT = source.SysRespT;
    out.SysTput = source.SysTput;
    for (std::size_t i = 0; i < source.Station.size(); ++i) {
        if (!keep(i)) continue;
        out.Station.push_back(source.Station[i]);
        out.JobClass.push_back(source.JobClass[i]);
        out.QLen.push_back(source.QLen[i]);
        out.Util.push_back(source.Util[i]);
        out.RespT.push_back(source.RespT[i]);
        out.ResidT.push_back(source.ResidT[i]);
        out.ArvR.push_back(source.ArvR[i]);
        out.Tput.push_back(source.Tput[i]);
        if (!source.has_ci()) continue;
        out.QLenCI.push_back(source.QLenCI[i]);
        out.UtilCI.push_back(source.UtilCI[i]);
        out.RespTCI.push_back(source.RespTCI[i]);
        out.ResidTCI.push_back(source.ResidTCI[i]);
        out.ArvRCI.push_back(source.ArvRCI[i]);
        out.TputCI.push_back(source.TputCI[i]);
    }
    return out;
}
 
}  // namespace
 
AvgTable AvgTable::filter_by(const std::string& name) const {
    return filter_avg_table(*this, [&](std::size_t i) {
        return Station[i] == name || JobClass[i] == name;
    });
}
 
AvgTable AvgTable::filter_by(const std::string& station, const std::string& jobclass) const {
    return filter_avg_table(*this, [&](std::size_t i) {
        return Station[i] == station && JobClass[i] == jobclass;
    });
}
 
AvgTable AvgTable::filter_by(const Node& node) const {
    return filter_avg_table(*this, [&](std::size_t i) {
        return Station[i] == node.get_name();
    });
}
 
AvgTable AvgTable::filter_by(const ::line::JobClass& jobclass) const {
    return filter_avg_table(*this, [&](std::size_t i) {
        return JobClass[i] == jobclass.get_name();
    });
}
 
AvgTable AvgTable::filter_by(const Node& node, const ::line::JobClass& jobclass) const {
    return filter_by(node.get_name(), jobclass.get_name());
}
 
AvgTable AvgTable::filter_by(const ::line::JobClass& jobclass, const Node& node) const {
    return filter_by(node, jobclass);
}
 
AvgTable AvgTable::get(const Node& node) const { return filter_by(node); }
AvgTable AvgTable::get(const ::line::JobClass& jobclass) const { return filter_by(jobclass); }
AvgTable AvgTable::get(const Node& node, const ::line::JobClass& jobclass) const {
    return filter_by(node, jobclass);
}
AvgTable AvgTable::get(const ::line::JobClass& jobclass, const Node& node) const {
    return filter_by(node, jobclass);
}
AvgTable AvgTable::tget(const Node& node) const { return filter_by(node); }
AvgTable AvgTable::tget(const ::line::JobClass& jobclass) const { return filter_by(jobclass); }
AvgTable AvgTable::tget(const Node& node, const ::line::JobClass& jobclass) const {
    return filter_by(node, jobclass);
}
AvgTable AvgTable::tget(const ::line::JobClass& jobclass, const Node& node) const {
    return filter_by(node, jobclass);
}
AvgTable AvgTable::operator()(const Node& node) const { return filter_by(node); }
AvgTable AvgTable::operator()(const ::line::JobClass& jobclass) const {
    return filter_by(jobclass);
}
AvgTable AvgTable::operator()(const Node& node, const ::line::JobClass& jobclass) const {
    return filter_by(node, jobclass);
}
AvgTable AvgTable::operator()(const ::line::JobClass& jobclass, const Node& node) const {
    return filter_by(node, jobclass);
}
 
void AvgTable::print(std::ostream& out) const {
    const std::ios::fmtflags flags = out.flags();
    const std::streamsize precision = out.precision();
    if (has_ci()) {
        // MATLAB's getAvgTable under options.confint: every cell reads `mean +- half-width`.
        const auto cell = [](double m, double h) {
            std::ostringstream c;
            c << std::setprecision(6) << m << " +- " << std::setprecision(3) << h;
            return c.str();
        };
        out << std::left << std::setw(16) << "Station" << std::setw(14) << "JobClass" << std::right;
        for (const char* h : {"QLen", "Util", "RespT", "ResidT", "ArvR", "Tput"}) out << std::setw(22) << h;
        out << '\n';
        for (std::size_t i = 0; i < Station.size(); ++i)
            out << std::left << std::setw(16) << Station[i] << std::setw(14) << JobClass[i]
                << std::right << std::setw(22) << cell(QLen[i], QLenCI[i])
                << std::setw(22) << cell(Util[i], UtilCI[i]) << std::setw(22) << cell(RespT[i], RespTCI[i])
                << std::setw(22) << cell(ResidT[i], ResidTCI[i]) << std::setw(22) << cell(ArvR[i], ArvRCI[i])
                << std::setw(22) << cell(Tput[i], TputCI[i]) << '\n';
        out.flags(flags);
        return;
    }
    out << std::left << std::setw(16) << "Station" << std::setw(14) << "JobClass"
        << std::right << std::setw(12) << "QLen" << std::setw(12) << "Util"
        << std::setw(12) << "RespT" << std::setw(12) << "ResidT"
        << std::setw(12) << "ArvR" << std::setw(12) << "Tput" << '\n';
    out << std::setprecision(6);
    for (std::size_t i = 0; i < Station.size(); ++i)
        out << std::left << std::setw(16) << Station[i] << std::setw(14) << JobClass[i]
            << std::right << std::setw(12) << QLen[i] << std::setw(12) << Util[i]
            << std::setw(12) << RespT[i] << std::setw(12) << ResidT[i]
            << std::setw(12) << ArvR[i] << std::setw(12) << Tput[i] << '\n';
    out.flags(flags);
    out.precision(precision);
}
 
void AvgTable::print() const { print(std::cout); }
 
std::ostream& operator<<(std::ostream& out, const AvgTable& table) {
    table.print(out);
    return out;
}
 
void AvgSysTable::print(std::ostream& out) const {
    const std::ios::fmtflags flags = out.flags();
    const std::streamsize precision = out.precision();
    std::size_t w = 12;
    for (const std::string& s : JobClasses) w = std::max(w, s.size() + 2);
    out << std::left << std::setw(10) << "Chain" << std::setw(int(w)) << "JobClasses" << std::right
        << std::setw(12) << "SysRespT" << std::setw(12) << "SysTput" << '\n';
    out << std::setprecision(6);
    for (std::size_t c = 0; c < Chain.size(); ++c)
        out << std::left << std::setw(10) << Chain[c] << std::setw(int(w)) << JobClasses[c] << std::right
            << std::setw(12) << SysRespT[c] << std::setw(12) << SysTput[c] << '\n';
    out.flags(flags);
    out.precision(precision);
}
 
void AvgSysTable::print() const { print(std::cout); }
 
std::ostream& operator<<(std::ostream& out, const AvgSysTable& table) {
    table.print(out);
    return out;
}
 
void AvgChainTable::print(std::ostream& out) const {
    const std::ios::fmtflags flags = out.flags();
    const std::streamsize precision = out.precision();
    std::size_t ws = 10, wj = 12;
    for (std::size_t i = 0; i < Station.size(); ++i) {
        ws = std::max(ws, Station[i].size() + 2);
        wj = std::max(wj, JobClasses[i].size() + 2);
    }
    out << std::left << std::setw(int(ws)) << "Station" << std::setw(10) << "Chain" << std::setw(int(wj))
        << "JobClasses" << std::right;
    for (const char* c : {"QLen", "Util", "RespT", "ResidT", "ArvR", "Tput"}) out << std::setw(12) << c;
    out << '\n' << std::setprecision(6);
    for (std::size_t i = 0; i < Station.size(); ++i)
        out << std::left << std::setw(int(ws)) << Station[i] << std::setw(10) << Chain[i] << std::setw(int(wj))
            << JobClasses[i] << std::right << std::setw(12) << QLen[i] << std::setw(12) << Util[i]
            << std::setw(12) << RespT[i] << std::setw(12) << ResidT[i] << std::setw(12) << ArvR[i]
            << std::setw(12) << Tput[i] << '\n';
    out.flags(flags);
    out.precision(precision);
}
 
void AvgChainTable::print() const { print(std::cout); }
 
std::ostream& operator<<(std::ostream& out, const AvgChainTable& table) {
    table.print(out);
    return out;
}
 
void AvgFJTable::print(std::ostream& out) const {
    const std::ios::fmtflags flags = out.flags();
    const std::streamsize precision = out.precision();
    std::size_t wf = 6, wj = 6, wc = 10;
    for (std::size_t i = 0; i < Fork.size(); ++i) {
        wf = std::max(wf, Fork[i].size() + 2);
        wj = std::max(wj, Join[i].size() + 2);
        wc = std::max(wc, JobClass[i].size() + 2);
    }
    out << std::left << std::setw(int(wf)) << "Fork" << std::setw(int(wj)) << "Join" << std::setw(int(wc))
        << "JobClass" << std::right << std::setw(14) << "FJQLen" << std::setw(14) << "FJRespT" << std::setw(14)
        << "FJTput" << '\n';
    out << std::setprecision(6);
    for (std::size_t i = 0; i < Fork.size(); ++i)
        out << std::left << std::setw(int(wf)) << Fork[i] << std::setw(int(wj)) << Join[i] << std::setw(int(wc))
            << JobClass[i] << std::right << std::setw(14) << FJQLen[i] << std::setw(14) << FJRespT[i]
            << std::setw(14) << FJTput[i] << '\n';
    out.flags(flags);
    out.precision(precision);
}
 
void AvgFJTable::print() const { print(std::cout); }
 
std::ostream& operator<<(std::ostream& out, const AvgFJTable& table) {
    table.print(out);
    return out;
}
 
namespace {
 
/** Fill the label columns and the per-class system metrics of a table. */
void fill_table(AvgTable& t, const qn::NetworkStruct<double>& sn, const mva::AvgResult<double>& r) {
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (std::size_t c = 0; c < sn.nclasses; ++c) {
            const double q = r.QN(i, c), u = r.UN(i, c), rt = r.RN(i, c);
            const double w = r.WN(i, c), a = r.AN(i, c), x = r.TN(i, c);
            if (q == 0.0 && u == 0.0 && rt == 0.0 && w == 0.0 && a == 0.0 && x == 0.0) continue;
            t.Station.push_back(sn.stations[i].name);
            t.JobClass.push_back(sn.classes[c].name);
            t.QLen.push_back(q);
            t.Util.push_back(u);
            t.RespT.push_back(rt);
            t.ResidT.push_back(w);
            t.ArvR.push_back(a);
            t.Tput.push_back(x);
        }
    for (std::size_t c = 0; c < sn.nclasses && c < r.CN.size(); ++c) {
        t.SysClass.push_back(sn.classes[c].name);
        t.SysRespT.push_back(r.CN[c]);
        t.SysTput.push_back(c < r.XN.size() ? r.XN[c] : 0.0);
    }
    t.iter = r.iter;
    t.warning = r.warning;
    t.ListCost = r.listcost;
    if (r.lognormconst.has_value()) {
        t.has_lognormconst = true;
        t.lognormconst = r.lognormconst.value();
    }
}
 
/**
 * The same fill for a solver whose solution type is not `mva::AvgResult`.
 *
 * THE SSA AND FLUID ARMS OWE BOTH DERIVED COLUMNS. Neither analyzer reports a
 * residence time or an arrival rate of its own, and neither is a copy of its
 * neighbour: ResidT is the per-JOB time and RespT the per-VISIT one, and they
 * agree only where every station is visited once per cycle. ArvR is a FLOW and
 * parts from the throughput at any station a job leaves by another route.
 * `sn_get_residt_from_respt` and `sn_get_arvr_from_tput` are the reference's own
 * conversions and are what the CLI's `-s fluid` and `-s ssa` arms apply. A
 * Source's ArvR stays zero: nothing arrives TO it.
 */
template <class Sol>
void fill_table_sim(AvgTable& t, const qn::NetworkStruct<double>& sn, const Sol& r) {
    const std::size_t M = sn.nstations, K = sn.nclasses;
    Matrix<double> RN(M, K), TN(M, K);
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t c = 0; c < K; ++c) {
            RN(i, c) = r.RN(i, c);
            TN(i, c) = r.TN(i, c);
        }
    const Matrix<double> WN = mva::sn_get_residt_from_respt<double>(sn, RN);
    const Matrix<double> AN = mva::sn_get_arvr_from_tput<double>(sn, TN);
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t c = 0; c < K; ++c) {
            const double q = r.QN(i, c), u = r.UN(i, c), rt = r.RN(i, c), x = r.TN(i, c);
            const bool is_source = sn.stations[i].sched == lang::SchedStrategy::EXT;
            const double w = WN(i, c), a = is_source ? 0.0 : AN(i, c);
            // The all-zero row test over ALL SIX columns, as `avg_rows` and the
            // reference make it: a station only ARRIVALS reach still has a row.
            if (q == 0.0 && u == 0.0 && rt == 0.0 && w == 0.0 && a == 0.0 && x == 0.0) continue;
            t.Station.push_back(sn.stations[i].name);
            t.JobClass.push_back(sn.classes[c].name);
            t.QLen.push_back(q);
            t.Util.push_back(u);
            t.RespT.push_back(rt);
            t.ResidT.push_back(w);
            t.ArvR.push_back(a);
            t.Tput.push_back(x);
        }
    for (std::size_t c = 0; c < sn.nclasses && c < r.CN.size(); ++c) {
        t.SysClass.push_back(sn.classes[c].name);
        t.SysRespT.push_back(r.CN[c]);
        t.SysTput.push_back(c < r.XN.size() ? r.XN[c] : 0.0);
    }
}
 
/**
 * The LDES fill. Its solution type carries the six metrics separately and its
 * per-class CN/XN are (1 x nclasses) MATRICES rather than vectors, so neither
 * `fill_table` nor `fill_table_sim` can read it.
 *
 * ResidT IS RECOMPUTED RATHER THAN TAKEN, which is what the CLI's `-a avg` LDES
 * arm does and for the same reason: the engine counts ONE VISIT PER STATION, so
 * the residence time it reports is its response time wherever a visit ratio is
 * not 1.
 */
void fill_table_ldes(AvgTable& t, const qn::NetworkStruct<double>& sn,
                     const ldes::LdesResult& r) {
    const Matrix<double> WNfix = mva::sn_get_residt_from_respt<double>(sn, r.RN);
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (std::size_t c = 0; c < sn.nclasses; ++c) {
            const double q = r.QN(i, c), u = r.UN(i, c), rt = r.RN(i, c);
            const double w = WNfix(i, c), a = r.AN(i, c), x = r.TN(i, c);
            if (q == 0.0 && u == 0.0 && rt == 0.0 && w == 0.0 && a == 0.0 && x == 0.0) continue;
            t.Station.push_back(sn.stations[i].name);
            t.JobClass.push_back(sn.classes[c].name);
            t.QLen.push_back(q);
            t.Util.push_back(u);
            t.RespT.push_back(rt);
            t.ResidT.push_back(w);
            t.ArvR.push_back(a);
            t.Tput.push_back(x);
        }
    for (std::size_t c = 0; c < sn.nclasses && c < r.CN.cols(); ++c) {
        t.SysClass.push_back(sn.classes[c].name);
        t.SysRespT.push_back(r.CN(0, c));
        t.SysTput.push_back(c < r.XN.cols() ? r.XN(0, c) : 0.0);
    }
}
 
/**
 * Append the confidence half-widths of the rows a fill just kept.
 *
 * The rows are re-selected by the SAME all-zero test over the six means, so row k here is row k
 * of the table. The ResidT half-width is the RespT one scaled by W/R: the residence time is the
 * response time times a constant visit ratio, so the interval scales with it.
 */
void fill_table_ci(AvgTable& t, const qn::NetworkStruct<double>& sn, const Matrix<double>& QN,
                   const Matrix<double>& UN, const Matrix<double>& RN, const Matrix<double>& WN,
                   const Matrix<double>& AN, const Matrix<double>& TN, const Matrix<double>& QCI,
                   const Matrix<double>& UCI, const Matrix<double>& RCI, const Matrix<double>& ACI,
                   const Matrix<double>& TCI) {
    const auto at = [](const Matrix<double>& m, std::size_t i, std::size_t c) {
        return i < m.rows() && c < m.cols() ? m(i, c) : 0.0;
    };
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (std::size_t c = 0; c < sn.nclasses; ++c) {
            const double q = QN(i, c), u = UN(i, c), rt = RN(i, c), w = WN(i, c), a = AN(i, c), x = TN(i, c);
            if (q == 0.0 && u == 0.0 && rt == 0.0 && w == 0.0 && a == 0.0 && x == 0.0) continue;
            t.QLenCI.push_back(at(QCI, i, c));
            t.UtilCI.push_back(at(UCI, i, c));
            t.RespTCI.push_back(at(RCI, i, c));
            t.ResidTCI.push_back(rt > 0.0 ? at(RCI, i, c) * w / rt : 0.0);
            t.ArvRCI.push_back(at(ACI, i, c));
            t.TputCI.push_back(at(TCI, i, c));
        }
}
 
/**
 * The caller-facing `map_env` knobs, as the gate wants them.
 *
 * They live on `SolverOptions` and NOT on `MvaOptions`/`NcSolverOptions`/
 * `FluidOptions`, because no `runAnalyzer` in the reference reads them: the
 * fallback is decided one level ABOVE the runner, in `getAvg`, which is what
 * this facade's `avg_table` plays the part of here.
 */
solvers::MapEnvConfig map_env_config(const SolverOptions& o) {
    solvers::MapEnvConfig c;
    if (!o.map_env.empty()) c.mode = o.map_env;
    if (!o.map_env_method.empty()) c.method = o.map_env_method;
    if (o.map_env_maxstages) c.max_stages = o.map_env_maxstages;
    return c;
}
 
ctmc::CtmcOptions ctmc_options(const SolverOptions& o) {
    ctmc::CtmcOptions c;
    if (!o.method.empty()) c.method = o.method;
    if (o.cutoff > 0.0) c.cutoff = o.cutoff;
    if (!o.cutoff_vec.empty()) c.cutoff_vec = o.cutoff_vec;
    if (o.state_max) c.state_max = o.state_max;
    if (o.timestep > 0.0) c.timestep = o.timestep;
    c.force = o.force;
    return c;
}
 
/**
 * The symbolic knobs, which live on their own struct rather than on
 * CtmcOptions: the backend and its timeout configure the computer-algebra
 * request, not the generator build. An empty or non-positive field keeps the
 * engine default, as everywhere else here.
 */
ctmc::CtmcSymbolicOptions ctmc_symbolic_options(const SolverOptions& o) {
    ctmc::CtmcSymbolicOptions s;
    if (!o.symbolic.empty()) s.backend = o.symbolic;
    if (o.symbolic_timeout > 0) s.timeout_s = o.symbolic_timeout;
    return s;
}
 
fluid::FluidOptions fluid_options(const SolverOptions& o) {
    fluid::FluidOptions f;
    if (!o.method.empty()) f.method = o.method;
    if (o.tol >= 0.0) f.tol = o.tol;
    if (o.iter_tol >= 0.0) f.iter_tol = o.iter_tol;
    if (o.iter_max >= 0) f.iter_max = static_cast<std::size_t>(o.iter_max);
    if (o.timespan_end > 0.0) f.timespan_end = o.timespan_end;
    if (!o.init_sol.empty()) f.init_sol = o.init_sol;
    if (o.seed) f.seed = o.seed;
    if (!o.highvar.empty()) f.highvar = o.highvar;
    if (o.timestep > 0.0) f.timestep = o.timestep;
    f.stiff = o.stiff;
    return f;
}
 
/**
 * The JMT client's options.
 *
 * `timespan_end` becomes `max_simulated_time`, which is the finite horizon the
 * replicated arms REQUIRE, and `replications` is `options.config.replications`, the
 * size of their ensemble. `iter_max` means nothing to JMT and is not read.
 */
jmt::JmtOptions jmt_options(const SolverOptions& o) {
    jmt::JmtOptions j;
    if (!o.method.empty()) j.method = o.method;
    if (o.samples) j.samples = static_cast<double>(o.samples);
    if (o.seed) j.seed = static_cast<long>(o.seed);
    if (o.replications > 0) j.replications = o.replications;
    if (o.timespan_end > 0.0) j.max_simulated_time = o.timespan_end;
    j.confint = o.confint;
    j.keep = o.keep;
    j.verbose = o.verbose;
    return j;
}
 
/** The LDES engine's options, shared by the average solve and `cdf_respt`. */
ldes::LdesOptions ldes_options(const SolverOptions& o) {
    ldes::LdesOptions opt;
    if (!o.method.empty()) opt.method = o.method;
    if (o.samples) opt.samples = o.samples;
    if (o.events) opt.events = o.events;
    if (o.seed) opt.seed = static_cast<long>(o.seed);
    if (o.timeout > 0.0) opt.timeout = o.timeout;
    if (o.confint > 0.0) opt.confint = o.confint;
    if (!o.cimethod.empty()) opt.cimethod = o.cimethod;
    if (o.ciminbatch >= 0) opt.ciminbatch = o.ciminbatch;
    if (o.ciminobs >= 0) opt.ciminobs = o.ciminobs;
    if (o.obmoverlap >= 0.0) opt.obmoverlap = o.obmoverlap;
    if (o.spectral_low_freq_frac > 0.0) opt.spectral_low_freq_frac = o.spectral_low_freq_frac;
    if (!o.tranfilter.empty()) opt.tranfilter = o.tranfilter;
    if (o.mserbatch > 0) opt.mserbatch = o.mserbatch;
    if (o.warmupfrac >= 0.0) opt.warmupfrac = o.warmupfrac;
    opt.cnvgon = o.cnvgon;
    if (o.cnvgtol > 0.0) opt.cnvgtol = o.cnvgtol;
    if (o.cnvgbatch > 0) opt.cnvgbatch = o.cnvgbatch;
    if (o.cnvgchk >= 0) opt.cnvgchk = o.cnvgchk;
    if (!o.init_sol.empty()) opt.init_sol = o.init_sol;
    if (o.init_jobs.rows() > 0) opt.init_jobs = o.init_jobs;
    opt.track_jobs = o.track_jobs;
    opt.verbose = o.verbose;
    return opt;
}
 
/**
 * The facade name of an AUTO method family, or a refusal naming what the family needs.
 *
 * Mirrors the CLI's `auto_cli_token_of_family`: a family names its engine outright. ln, env and
 * uq take a LayeredNetwork, an Environment or a sweep, none of which a Network facade holds.
 */
std::string facade_name_of_family(const std::string& fam) {
    if (fam == "mva" || fam == "nc" || fam == "ctmc" || fam == "mam" || fam == "ag" ||
        fam == "ba" || fam == "ssa" || fam == "jmt" || fam == "lqns") {
        std::string up = fam;
        std::transform(up.begin(), up.end(), up.begin(), ::toupper);
        return up;
    }
    if (fam == "fld") return "FLD";
    if (fam == "ldes") {
        if (!ldes::ldes_is_available())
            throw ExternalError(ErrorKind{"LDES"},
                "SolverAUTO: method family 'ldes' names the discrete-event engine, and none was "
                "found (common/ldes or common/ldes.jar, or $LINE_LDES_DIR); 'ssa' is the simulator "
                "this port builds in process");
        return "LDES";
    }
    if (fam == "ln")
        throw UnsupportedError(ErrorKind{"Option"},
            "SolverAUTO: method family 'ln' names the layered solver, which takes a "
            "LayeredNetwork rather than a Network");
    if (fam == "env")
        throw UnsupportedError(ErrorKind{"Option"},
            "SolverAUTO: method family 'env' names the random-environment solver, which takes an "
            "Environment rather than a Network");
    if (fam == "uq")
        throw UnsupportedError(ErrorKind{"Option"},
            "SolverAUTO: method family 'uq' names the uncertainty-quantification solver, which "
            "this facade does not wrap");
    throw InputError("SolverAUTO: no engine stands behind method family '" + fam + "'");
}
 
/** A row matrix as a vector; empty stays empty, meaning not computed. */
std::vector<double> row_of(const Matrix<double>& m) {
    std::vector<double> v;
    if (m.rows() > 0)
        for (std::size_t c = 0; c < m.cols(); ++c) v.push_back(m(0, c));
    return v;
}
 
/** The model's Cache nodes with nothing measured yet: caps, item counts and item sizes only. */
solvers::CacheMetrics<double> empty_cache_metrics(const qn::NetworkStruct<double>& sn) {
    return solvers::cache_metrics_of(sn, std::vector<double>(), std::vector<double>(),
                                     std::vector<double>(), std::vector<double>(), Matrix<double>(),
                                     Matrix<double>(), std::vector<double>());
}
 
/** JMT's per-node hit probabilities; the miss is 1 - hit, as SolverJMT/runAnalyzer.m derives it. */
solvers::CacheMetrics<double> cache_metrics_of_jmt(const qn::NetworkStruct<double>& sn,
                                                   const std::map<std::size_t, std::vector<double> >& hit) {
    solvers::CacheMetrics<double> out = empty_cache_metrics(sn);
    for (std::size_t c = 0; c < out.caches.size(); ++c) {
        solvers::CacheNodeMetrics<double>& m = out.caches[c];
        const std::map<std::size_t, std::vector<double> >::const_iterator it = hit.find(m.node);
        if (it == hit.end()) continue;
        m.hitprob = it->second;
        for (std::size_t r = 0; r < m.hitprob.size(); ++r) m.missprob.push_back(1.0 - m.hitprob[r]);
        solvers::cache_mark_nonreading(sn, m);
    }
    return out;
}
 
/** The LDES `cacheMetrics` block, keyed by Cache NAME (see `CacheNodeMetrics::name`). */
solvers::CacheMetrics<double> cache_metrics_of_ldes(const qn::NetworkStruct<double>& sn,
                                                    const std::map<std::string, ldes::LdesCacheMetrics>& cm) {
    solvers::CacheMetrics<double> out = empty_cache_metrics(sn);
    for (std::size_t c = 0; c < out.caches.size(); ++c) {
        solvers::CacheNodeMetrics<double>& m = out.caches[c];
        const std::map<std::string, ldes::LdesCacheMetrics>::const_iterator it = cm.find(m.name);
        if (it == cm.end()) continue;
        m.hitprob = row_of(it->second.hit);
        m.missprob = row_of(it->second.miss);
        m.delayedprob = row_of(it->second.delayed);
        m.latency = row_of(it->second.latency);
        m.hitproblist = it->second.hitList;
        m.itemprob = it->second.itemProb;
        m.listcost = row_of(it->second.listCost);
        solvers::cache_mark_nonreading(sn, m);
    }
    return out;
}
 
/**
 * One engine's average table, by facade name. `o.method` is that engine's own method name.
 * `cache` receives what the engine observed at the Cache nodes, empty where it analyzes none.
 */
void solve_named(const std::string& name, const SolverOptions& o,
                 const qn::NetworkStruct<double>& sn, AvgTable& t,
                 solvers::CacheMetrics<double>& cache, solvers::FjMetrics<double>* fj = nullptr) {
    cache = solvers::CacheMetrics<double>();
    if (fj) *fj = solvers::FjMetrics<double>();
    t.solver = name;
    if (o.confint < 0.0 || o.confint >= 1.0)
        throw InputError("SolverOptions: confint must be a confidence level in (0,1), or 0 to disable; got " +
                         std::to_string(o.confint));
    if (name == "MVA") {
        mva::MvaOptions opt;
        if (!o.method.empty()) opt.method = o.method;
        if (o.tol >= 0.0) opt.tol = o.tol;
        if (o.iter_tol >= 0.0) opt.iter_tol = o.iter_tol;
        if (o.iter_max >= 0) opt.iter_max = o.iter_max;
        if (!o.multiserver.empty()) opt.multiserver = o.multiserver;
        if (!o.highvar.empty()) opt.highvar = o.highvar;
        if (!o.np_priority.empty()) opt.np_priority = o.np_priority;
        if (!o.fork_join.empty()) opt.fork_join = o.fork_join;
        Matrix<double> init;
        const mva::AvgResult<double> r = solvers::run_avg<double>(
            sn, "SolverMVA", qn::mva_feature_set(mva::resolve_method(sn, opt.method), sn),
            map_env_config(o),
            [&opt, &init](const qn::NetworkStruct<double>& m) {
                return mva::solver_mva_run_analyzer(m, opt, init);
            },
            solvers::mva_stage_fn(opt), opt.method);
        t.method = r.actualmethod;
        fill_table(t, sn, r);
        cache = r.cache;
        if (fj) *fj = r.fj;
    } else if (name == "NC") {
        nc::NcSolverOptions opt;
        if (!o.method.empty()) opt.method = o.method;
        if (o.tol >= 0.0) opt.tol = o.tol;
        if (o.iter_tol >= 0.0) opt.iter_tol = o.iter_tol;
        if (o.iter_max >= 0) opt.iter_max = o.iter_max;
        if (!o.highvar.empty()) opt.highvar = o.highvar;
        if (!o.multiserver.empty()) opt.multiserver = o.multiserver;
        if (!o.fork_join.empty()) opt.fork_join = o.fork_join;
        if (o.samples) opt.samples = o.samples;
        if (o.seed) opt.seed = o.seed;
        const mva::AvgResult<double> r = solvers::run_avg<double>(
            sn, "SolverNC", qn::nc_feature_set(opt.method), map_env_config(o),
            [&opt](const qn::NetworkStruct<double>& m) { return nc::solver_nc_run_analyzer(m, opt); },
            solvers::nc_stage_fn(opt), opt.method);
        t.method = r.actualmethod;
        fill_table(t, sn, r);
        cache = r.cache;
        if (fj) *fj = r.fj;
    } else if (name == "CTMC") {
        const ctmc::CtmcOptions opt = ctmc_options(o);
        // `mdd` never forms the generator, so it has its own analyzer, as the CLI's `-s ctmc` arm routes it.
        const mva::AvgResult<double> r = ctmc::is_stateless_method(opt.method)
                                             ? ctmc::solver_ctmc_mdd_run_analyzer(sn, opt)
                                             : ctmc::solver_ctmc_run_analyzer_any(sn, opt);
        t.method = r.actualmethod;
        fill_table(t, sn, r);
        cache = r.cache;
    } else if (name == "MAM") {
        mam::MamOptions opt;
        if (!o.method.empty()) opt.method = o.method;
        if (o.tol >= 0.0) opt.tol = o.tol;
        if (o.iter_max >= 0) opt.iter_max = o.iter_max;
        const mva::AvgResult<double> r = mam::solver_mam_run_analyzer(sn, opt);
        t.method = r.actualmethod;
        fill_table(t, sn, r);
        cache = r.cache;
    } else if (name == "BA") {
        // DELIBERATELY NOT WRAPPED IN `run_avg`. A bound request must be answered
        // with a bound, and the environment image is an APPROXIMATION of the
        // model, so bounds computed on it do not bracket the original. The
        // reference refuses SolverBA inside `needsMapEnv` itself; here the
        // exclusion is the absence of the wrapper, so do not "fix" the
        // inconsistency with the arms above.
        ba::BaOptions opt;
        if (!o.method.empty()) opt.method = o.method;
        opt.level = o.level;
        const mva::AvgResult<double> r = ba::solver_ba_run_analyzer(sn, opt);
        t.method = r.actualmethod;
        fill_table(t, sn, r);
        cache = r.cache;
    } else if (name == "SSA") {
        ssa::SsaOptions opt;
        if (!o.method.empty()) opt.method = o.method;
        if (o.samples) opt.samples = o.samples;
        if (o.events) opt.samples = o.events;
        if (o.seed) opt.seed = o.seed;
        if (o.cutoff > 0.0) opt.cutoff = o.cutoff;
        if (o.warmupfrac >= 0.0) opt.warmupfrac = o.warmupfrac;
        if (!o.state_space_gen.empty()) opt.state_space_gen = o.state_space_gen;
        opt.verbose = o.verbose;
        // The derived columns are taken on the MEASURED cache split, as the
        // CLI's `-s ssa` arm does: the base struct carries `link()`'s even
        // hit/miss offer, so ArvR and ResidT read off it are not the model's.
        std::vector<ssa::SsaCacheRatio> cacheratio;
        const ssa::SsaSolution r = ssa::solver_ssa(sn, opt, &cacheratio);
        t.method = r.method;
        fill_table_sim(t, ssa::sn_with_ssa_cache_split<double>(sn, cacheratio), r);
        cache = ssa::cache_metrics_of_ssa<double>(sn, cacheratio);
    } else if (name == "FLD" || name == "FLUID") {
        const fluid::FluidOptions opt = fluid_options(o);
        // The fluid arm cannot go through `run_avg`: its runner reports a
        // `FluidSolution` and the driver an `AvgResult`, so the two halves fill
        // the table through different helpers. The DECISION is the same one
        // `run_avg` makes, asked here directly rather than duplicated.
        const solvers::MapEnvConfig mecfg = map_env_config(o);
        const solvers::MapEnvDecision d =
            solvers::needs_map_env(qn::fluid_feature_set(opt.method), sn, mecfg);
        if (d.needed) {
            const mva::AvgResult<double> r = solvers::map_env_approx<double>(
                sn, "SolverFLD", mecfg, solvers::fluid_stage_fn(opt), opt.method);
            t.method = r.actualmethod;
            fill_table(t, sn, r);
            cache = r.cache;
            if (fj) *fj = r.fj;
        } else {
            // The refreshed struct carries the CONVERGED cache split; ArvR read off `sn` would report
            // link()'s offered one instead (cache_replc_routing: 0.5/0.5 for 0.4/0.6), as in `-s fluid`.
            qn::NetworkStruct<double> refreshed;
            const fluid::FluidSolution r = fluid::solver_fluid_run_analyzer(
                sn, opt, static_cast<qn::NetworkStruct<double>*>(nullptr), &refreshed, &cache);
            t.method = r.method;
            t.iter = static_cast<int>(r.iters);
            t.warning = r.warning;
            if (fj) *fj = r.fj;
            const bool has_ref = refreshed.nstations == sn.nstations && refreshed.nclasses == sn.nclasses;
            fill_table_sim(t, has_ref ? refreshed : sn, r);
        }
    } else if (name == "JMT") {
        const jmt::JmtOptions opt = jmt_options(o);
        const jmt::JmtResult<double> r = jmt::solver_jmt_run_analyzer(sn, opt);
        t.method = r.avg.actualmethod;
        fill_table(t, sn, r.avg);
        cache = r.cache_hit_prob.empty() ? r.avg.cache : cache_metrics_of_jmt(sn, r.cache_hit_prob);
        if (fj) *fj = r.avg.fj;
        if (o.confint > 0.0 && !r.QCI.empty())
            fill_table_ci(t, sn, r.avg.QN, r.avg.UN, r.avg.RN, r.avg.WN, r.avg.AN, r.avg.TN, r.QCI,
                          r.UCI, r.RCI, r.ACI, r.TCI);
    } else if (name == "LDES") {
        const ldes::LdesOptions opt = ldes_options(o);
        const ldes::LdesResult r = ldes::solver_ldes(sn, opt);
        t.method = opt.method;
        fill_table_ldes(t, sn, r);
        cache = cache_metrics_of_ldes(sn, r.cache_metrics);
        if (fj) *fj = ldes::ldes_fj_metrics(sn, r);
        if (o.confint > 0.0 && !r.QNCI.empty())
            fill_table_ci(t, sn, r.QN, r.UN, r.RN, mva::sn_get_residt_from_respt<double>(sn, r.RN), r.AN,
                          r.TN, r.QNCI, r.UNCI, r.RNCI, r.ANCI, r.TNCI);
    } else if (name == "AG") {
        ag::AgOptions opt;
        if (!o.method.empty() && o.method != "default") opt.method = o.method;
        if (o.tol >= 0.0) opt.tol = o.tol;
        if (o.iter_max >= 0) opt.iter_max = o.iter_max;
        const mva::AvgResult<double> r = ag::solver_ag_run_analyzer(sn, opt);
        t.method = r.actualmethod;
        fill_table(t, sn, r);
        cache = r.cache;
    } else if (name == "LQNS") {
        // On a flat Network the lqns family's methods are the qns ones, served by `qnsolver`.
        lqns::QnsOptions opt;
        if (!o.method.empty()) opt.method = o.method;
        if (!o.multiserver.empty()) opt.multiserver = o.multiserver;
        opt.keep = o.keep;
        const mva::AvgResult<double> r = lqns::solve_network_run_analyzer(sn, opt);
        t.method = r.actualmethod;
        fill_table(t, sn, r);
        cache = r.cache;
    } else {
        throw UnsupportedError(ErrorKind{"Solver"}, "avg_table: unknown solver '" + name + "'");
    }
}
 
/**
 * SolverAUTO's average table, as `-s auto` computes it (line_cli.cpp, `choose_auto_plan`).
 *
 * `o.method` is resolved first: a method FAMILY ('nc', 'nc.comom', 'bound') names the engine and
 * passes the submethod to it; an INTENT ('default', 'sim', 'exact', 'fast', 'accurate', 'heur')
 * picks the ranking, whose proposed order is a RETRY LIST. Only the first attempt carries the
 * method the ranking gated on. Input and internal errors propagate; the other categories move on,
 * and when every candidate fails the chosen engine's error is reported.
 */
void solve_auto(const SolverOptions& o, const qn::NetworkStruct<double>& sn, AvgTable& t,
                solvers::CacheMetrics<double>& cache, solvers::FjMetrics<double>* fj = nullptr) {
    const autosolver::AutoToken tok = autosolver::auto_resolve_token(o.method, sn);
    SolverOptions oo = o;
    if (!tok.is_intent) {
        oo.method = tok.submethod;
        solve_named(facade_name_of_family(tok.family), oo, sn, t, cache, fj);
        return;
    }
    const autosolver::AutoChoice c = autosolver::auto_choose_solver_mode(sn, "getAvgTable", tok.mode);
    const std::vector<autosolver::AutoSolver> order =
        autosolver::auto_proposed_solvers(sn, "getAvgTable", tok.mode);
    std::string first_error, first_category, first_kind;
    for (std::size_t i = 0; i < order.size(); ++i) {
        std::string name = autosolver::auto_solver_name(order[i]);
        std::transform(name.begin(), name.end(), name.begin(), ::toupper);
        if (name == "FLUID") name = "FLD";
        oo.method = (i == 0 && !c.method.empty()) ? c.method : std::string("default");
        try {
            AvgTable attempt;
            solvers::CacheMetrics<double> attempt_cache;
            solvers::FjMetrics<double> attempt_fj;
            solve_named(name, oo, sn, attempt, attempt_cache, &attempt_fj);
            t = attempt;
            cache = attempt_cache;
            if (fj) *fj = attempt_fj;
            return;
        } catch (const Error& e) {
            const std::string cat = e.category();
            if (cat == "input" || cat == "internal") throw;
            if (first_error.empty()) {
                first_error = name + ": " + e.what();
                first_category = cat;
                first_kind = e.kind();
            }
        } catch (const std::exception& e) {
            if (first_error.empty()) first_error = name + ": " + e.what();
        }
    }
    const std::string all_failed =
        "SolverAUTO: every candidate refused this run. The chosen engine reported -- " + first_error;
    if (first_category == "numerical") throw NumericError(ErrorKind{first_kind}, all_failed);
    if (first_category == "external") throw ExternalError(ErrorKind{first_kind}, all_failed);
    throw UnsupportedError(ErrorKind{first_kind.empty() ? std::string("Solver") : first_kind},
                           all_failed);
}
 
/** Defined below `avg_result_of_table`, which it reads its solve back through. */
void solve_zero_capacity(const std::string& name, const SolverOptions& o,
                         const qn::NetworkStruct<double>& sn, const api::ZeroCapacity& zc,
                         AvgTable& t, solvers::CacheMetrics<double>& cache,
                         solvers::FjMetrics<double>* fj);
 
}  // namespace
 
// ---------------------------------------------------------------------------
// NetworkSolver
// ---------------------------------------------------------------------------
 
const AvgTable& NetworkSolver::avg_table() {
    // options.cache: a stored result is reused only while caching is on, as getAvg does.
    if (solved_ && opts_.cache) return table_;
    AvgTable t;
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    // a pairwise switchover at a non-polling queue is kept in the model but no solver simulates it, as the reference warns
    for (const auto& st : sn.stations)
        for (const auto& row : st.switchover_pair)
            if (std::any_of(row.begin(), row.end(), [](const lang::Distrib<double>& d) {
                    return !d.disabled && d.type != lang::ProcessType::IMMEDIATE;
                })) {
                std::cerr << "[LINE] Warning: Station " << st.name << " declares class-to-class switchover "
                          << "times, which no LINE solver simulates at a non-polling queue: every solver "
                          << "answers as if they were zero. Only POLLING switchovers are solved." << std::endl;
                break;
            }
    solvers::CacheMetrics<double> cache;
    solvers::FjMetrics<double> fj;
    // A closed chain visiting a station with zero servers is trapped there;
    // see sn_zero_capacity. This is `@NetworkSolver/getAvg`'s hook, ahead of the
    // map_env funnel each engine applies, so every Network solver takes it.
    const api::ZeroCapacity zc = api::sn_zero_capacity(sn);
    if (!zc.reason.empty()) throw UnsupportedError(ErrorKind{"Structure"}, zc.reason);
    if (zc.any())
        solve_zero_capacity(name_, opts_, sn, zc, t, cache, &fj);
    else if (name_ == "AUTO")
        solve_auto(opts_, sn, t, cache, &fj);
    else
        solve_named(name_, opts_, sn, t, cache, &fj);
    table_ = t;
    cache_ = cache;
    fj_ = fj;
    solved_ = true;
    return table_;
}
 
std::string NetworkSolver::method_used() { return avg_table().method; }
 
namespace {
 
/**
 * The solved table laid back onto the station x class matrices that getAvgSys and getAvgChain read;
 * a row the table dropped was all-zero, which is what the fresh matrices already hold.
 */
mva::AvgResult<double> avg_result_of_table(const qn::NetworkStruct<double>& sn, const AvgTable& t) {
    const std::size_t M = sn.nstations, K = sn.nclasses;
    mva::AvgResult<double> r;
    r.QN = Matrix<double>(M, K);
    r.UN = Matrix<double>(M, K);
    r.RN = Matrix<double>(M, K);
    r.TN = Matrix<double>(M, K);
    r.AN = Matrix<double>(M, K);
    r.WN = Matrix<double>(M, K);
    for (std::size_t row = 0; row < t.size(); ++row)
        for (std::size_t i = 0; i < M; ++i) {
            if (sn.stations[i].name != t.Station[row]) continue;
            for (std::size_t c = 0; c < K; ++c) {
                if (sn.classes[c].name != t.JobClass[row]) continue;
                r.QN(i, c) = t.QLen[row];
                r.UN(i, c) = t.Util[row];
                r.RN(i, c) = t.RespT[row];
                r.TN(i, c) = t.Tput[row];
                r.AN(i, c) = t.ArvR[row];
                r.WN(i, c) = t.ResidT[row];
            }
        }
    return r;
}
 
/**
 * `getAvgZeroCapacity`: solve the model without the chains trapped at a
 * zero-capacity station, whose jobs hold no other resource, then park each
 * trapped population there.
 *
 * The reduced model is a COPY of the struct, refreshed as `get_struct` does,
 * where the reference edits the model and restores it on cleanup: the user's
 * model is never touched, so there is nothing to restore. The trapped classes
 * have population 0 and a zero-server station one server (no remaining chain
 * visits it); a rate-zero station keeps its servers. When every class is
 * trapped nothing is solved. At the dead station the trapped chain holds its
 * whole population, split by visits, with RespT = ResidT = Inf for the classes
 * that visit it; every other metric of the chain is 0, its system response
 * time Inf and its system throughput 0.
 */
void solve_zero_capacity(const std::string& name, const SolverOptions& o,
                         const qn::NetworkStruct<double>& sn, const api::ZeroCapacity& zc,
                         AvgTable& t, solvers::CacheMetrics<double>& cache,
                         solvers::FjMetrics<double>* fj) {
    const std::size_t M = sn.nstations, K = sn.nclasses;
    const double inf = std::numeric_limits<double>::infinity();
    std::vector<bool> trappedcls(K, false);
    std::size_t ntrapped = 0;
    for (std::size_t c = 0; c < sn.nchains; ++c)
        if (zc.deadstat[c] > 0)
            for (std::size_t r = 0; r < K; ++r)
                if (sn.chains[c][r] && !trappedcls[r]) {
                    trappedcls[r] = true;
                    ++ntrapped;
                }
    mva::AvgResult<double> r;
    AvgTable solved;
    if (ntrapped == K) {
        r.QN = Matrix<double>(M, K, 0.0);
        r.UN = r.QN;
        r.RN = r.QN;
        r.TN = r.QN;
        r.AN = r.QN;
        r.WN = r.QN;
        solved.solver = name;
        solved.method = o.method.empty() ? std::string("default") : o.method;
        cache = solvers::CacheMetrics<double>();
        if (fj) *fj = solvers::FjMetrics<double>();
    } else {
        qn::NetworkStruct<double> reduced = sn;
        for (std::size_t k = 0; k < K; ++k)
            if (trappedcls[k]) reduced.classes[k].population = 0.0;
        for (std::size_t c = 0; c < sn.nchains; ++c) {
            const std::size_t i = zc.deadstat[c];
            if (i > 0 && sn.stations[i - 1].nservers == 0.0) reduced.stations[i - 1].nservers = 1.0;
        }
        reduced.refresh_struct();
        api::sn_fj_nodevisits_mmt(reduced);
        if (name == "AUTO")
            solve_auto(o, reduced, solved, cache, fj);
        else
            solve_named(name, o, reduced, solved, cache, fj);
        r = avg_result_of_table(sn, solved);
    }
    std::vector<double> CN(K, 0.0), XN(K, 0.0);
    for (std::size_t row = 0; row < solved.SysClass.size(); ++row)
        for (std::size_t k = 0; k < K; ++k)
            if (sn.classes[k].name == solved.SysClass[row]) {
                CN[k] = solved.SysRespT[row];
                XN[k] = solved.SysTput[row];
            }
    const std::vector<double> st2sf = sn.station_to_stateful();
    const std::vector<double> njobs = sn.njobs();
    for (std::size_t c = 0; c < sn.nchains; ++c) {
        if (zc.deadstat[c] == 0) continue;
        const std::size_t i = zc.deadstat[c] - 1;
        const std::size_t sf = static_cast<std::size_t>(st2sf[i]) - 1;
        double pop = 0.0, vtot = 0.0;
        for (std::size_t k = 0; k < K; ++k)
            if (sn.chains[c][k]) {
                pop += njobs[k];
                vtot += sn.visits[c](sf, k);
            }
        for (std::size_t k = 0; k < K; ++k) {
            if (!sn.chains[c][k]) continue;
            for (std::size_t j = 0; j < M; ++j) {
                r.QN(j, k) = 0.0;
                r.UN(j, k) = 0.0;
                r.RN(j, k) = 0.0;
                r.TN(j, k) = 0.0;
                r.AN(j, k) = 0.0;
                r.WN(j, k) = 0.0;
            }
            const double v = sn.visits[c](sf, k);
            r.QN(i, k) = pop * v / vtot;
            if (v > 0.0) {
                r.RN(i, k) = inf;
                r.WN(i, k) = inf;
            }
            XN[k] = 0.0;
            CN[k] = inf;
        }
    }
    r.CN = CN;
    r.XN = XN;
    r.iter = solved.iter;
    r.warning = solved.warning;
    r.listcost = solved.ListCost;
    if (solved.has_lognormconst) r.lognormconst = solved.lognormconst;
    t = AvgTable();
    t.solver = solved.solver;
    t.method = solved.method;
    fill_table(t, sn, r);
}
 
}  // namespace
 
AvgSysTable NetworkSolver::avg_sys_table() {
    const AvgTable& t = avg_table();
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const mva::AvgResult<double> r = avg_result_of_table(sn, t);
    const solvers::SysResult<double> sys = solvers::solver_get_avg_sys<double>(sn, r);
    AvgSysTable out;
    out.Chain = solvers::chain_names(sn.nchains);
    out.JobClasses = solvers::chain_class_labels<double>(sn);
    out.SysRespT = sys.CN;
    out.SysTput = sys.XN;
    return out;
}
 
AvgChain NetworkSolver::avg_chain() {
    const AvgTable& t = avg_table();
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const solvers::ChainResult<double> c = solvers::solver_get_avg_chain<double>(sn, avg_result_of_table(sn, t));
    AvgChain out;
    out.QN = c.QN;
    out.UN = c.UN;
    out.RN = c.RN;
    out.WN = c.WN;
    out.AN = c.AN;
    out.TN = c.TN;
    return out;
}
 
AvgChainTable NetworkSolver::avg_chain_table() {
    // getAvgChainTable.m: row-major over (station, chain), every row kept.
    const AvgChain c = avg_chain();
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const std::vector<std::string> chains = solvers::chain_names(sn.nchains);
    const std::vector<std::string> labels = solvers::chain_class_labels<double>(sn);
    AvgChainTable out;
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (std::size_t k = 0; k < sn.nchains; ++k) {
            out.Station.push_back(sn.stations[i].name);
            out.Chain.push_back(chains[k]);
            out.JobClasses.push_back(labels[k]);
            out.QLen.push_back(c.QN(i, k));
            out.Util.push_back(c.UN(i, k));
            out.RespT.push_back(c.RN(i, k));
            out.ResidT.push_back(c.WN(i, k));
            out.ArvR.push_back(c.AN(i, k));
            out.Tput.push_back(c.TN(i, k));
        }
    return out;
}
 
const solvers::CacheMetrics<double>& NetworkSolver::cache_metrics() {
    avg_table();
    return cache_;
}
 
const solvers::FjMetrics<double>& NetworkSolver::avg_fj() {
    avg_table();
    return fj_;
}
 
AvgFJTable NetworkSolver::avg_fj_table() {
    const solvers::FjMetrics<double>& fj = avg_fj();
    AvgFJTable out;
    if (fj.empty()) return out;
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    for (std::size_t p = 0; p < fj.npairs(); ++p)
        for (std::size_t r = 0; r < fj.QLen.cols(); ++r) {
            const double q = fj.QLen(p, r), rt = fj.RespT(p, r), x = fj.Tput(p, r);
            if (!(q > 0.0 || rt > 0.0 || x > 0.0)) continue;  // getAvgFJTable keeps a row when any value is > 0
            out.Fork.push_back(sn.nodes[fj.forks[p] - 1].name);
            out.Join.push_back(sn.nodes[fj.joins[p] - 1].name);
            out.JobClass.push_back(sn.classes[r].name);
            out.FJQLen.push_back(q);
            out.FJRespT.push_back(rt);
            out.FJTput.push_back(x);
        }
    return out;
}
 
namespace {
 
/** The named Cache's per-class hit (or miss) vector, refused when the solve measured none. */
std::vector<double> cache_ratio(const solvers::CacheMetrics<double>& cm, const std::string& solver,
                                const std::string& cache, bool hit) {
    for (std::size_t c = 0; c < cm.caches.size(); ++c) {
        if (cm.caches[c].name != cache) continue;
        const std::vector<double>& v = hit ? cm.caches[c].hitprob : cm.caches[c].missprob;
        if (v.empty())
            throw UnsupportedError(ErrorKind{"Analysis"},
                                   "Solver" + solver + " measured no " + (hit ? "hit" : "miss") +
                                       " ratio at Cache '" + cache + "' with this method");
        return v;
    }
    throw InputError("'" + cache + "' is not a Cache node of this model");
}
 
}  // namespace
 
std::vector<double> NetworkSolver::hit_ratio(const std::string& cache) {
    return cache_ratio(cache_metrics(), name_, cache, true);
}
 
std::vector<double> NetworkSolver::miss_ratio(const std::string& cache) {
    return cache_ratio(cache_metrics(), name_, cache, false);
}
 
std::vector<double> NetworkSolver::hit_ratio(const Node& cache) { return hit_ratio(cache.get_name()); }
 
std::vector<double> NetworkSolver::miss_ratio(const Node& cache) { return miss_ratio(cache.get_name()); }
 
std::vector<std::string> NetworkSolver::list_valid_methods() const {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    if (name_ == "MVA") return mva::list_valid_methods(sn);
    if (name_ == "NC") return nc::list_valid_methods();
    if (name_ == "CTMC") return ctmc::list_valid_methods();
    if (name_ == "MAM") return mam::list_valid_methods();
    // The model-aware overload, not the bare one: the reduction bounds are
    // derived for a single-class closed network of single servers and the three
    // open-network bounds for its mirror image, so the list a caller may act on
    // depends on the model. `SolverBA` gates on exactly this list.
    if (name_ == "BA") return ba::list_valid_methods(sn);
    if (name_ == "FLD") return fluid::fluid_list_valid_methods();
    if (name_ == "SSA") return ssa::list_valid_methods();
    if (name_ == "JMT") return jmt::jmt_list_valid_methods();
    if (name_ == "LDES") return ldes::list_valid_methods();
    // AUTO answers about every family it can delegate to, gated by each one's
    // feature set on this model; see `auto_methods.h`.
    if (name_ == "AUTO") return autosolver::auto_list_valid_methods(sn);
    throw UnsupportedError(ErrorKind{"Solver"}, "list_valid_methods: unknown solver '" + name_ + "'");
}
 
// ---------------------------------------------------------------------------
// SolverCTMC
// ---------------------------------------------------------------------------
 
Matrix<double> SolverCTMC::state_space() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const ctmc::CtmcSolution<double> d = ctmc::solver_ctmc_analyzer(sn, ctmc_options(opts_));
    return ctmc::ctmc_get_state_space(sn, d).flat;
}
 
Matrix<double> SolverCTMC::generator() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    ctmc::CtmcOptions opt = ctmc_options(opts_);
    opt.keep_filtration = true;
    const ctmc::CtmcSolution<double> d = ctmc::solver_ctmc_analyzer(sn, opt);
    return ctmc::ctmc_get_infgen(sn, d).Q;
}
 
double SolverCTMC::prob_aggr(std::size_t node, const std::vector<double>& state) {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const ctmc::CtmcSolution<double> d = ctmc::solver_ctmc_analyzer(sn, ctmc_options(opts_));
    const Matrix<double> A = ctmc::ctmc_get_state_space_aggr(sn, d);
    const std::size_t ist = sn.nodes[node - 1].station;
    if (ist == 0) throw InputError("prob_aggr: the node is not a station");
    const std::size_t K = sn.nclasses;
    if (state.size() != K) throw InputError("prob_aggr: the state must have one entry per class");
    double p = 0.0;
    for (std::size_t s = 0; s < A.rows(); ++s) {
        bool hit = true;
        for (std::size_t k = 0; k < K && hit; ++k) hit = A(s, (ist - 1) * K + k) == state[k];
        if (hit) p += d.pi[s];
    }
    return p;
}
 
std::vector<double> SolverCTMC::marg_aggr(std::size_t node) {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const ctmc::CtmcSolution<double> d = ctmc::solver_ctmc_analyzer(sn, ctmc_options(opts_));
    const Matrix<double> A = ctmc::ctmc_get_state_space_aggr(sn, d);
    const std::size_t ist = sn.nodes[node - 1].station;
    if (ist == 0) throw InputError("marg_aggr: the node is not a station");
    const std::size_t K = sn.nclasses;
    std::size_t nmax = 0;
    for (std::size_t s = 0; s < A.rows(); ++s) {
        double tot = 0.0;
        for (std::size_t k = 0; k < K; ++k) tot += A(s, (ist - 1) * K + k);
        nmax = std::max(nmax, static_cast<std::size_t>(tot));
    }
    std::vector<double> pmf(nmax + 1, 0.0);
    for (std::size_t s = 0; s < A.rows(); ++s) {
        double tot = 0.0;
        for (std::size_t k = 0; k < K; ++k) tot += A(s, (ist - 1) * K + k);
        pmf[static_cast<std::size_t>(tot)] += d.pi[s];
    }
    return pmf;
}
 
std::vector<std::vector<CdfCurve> > SolverCTMC::cdf_respt() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    std::vector<std::vector<CdfCurve> > out(sn.nstations, std::vector<CdfCurve>(sn.nclasses));
    const std::vector<std::vector<ctmc::CdfCurve<double> > > R =
        ctmc::solver_ctmc_cdf_respt(sn, ctmc_options(opts_));
    for (std::size_t i = 0; i < R.size() && i < out.size(); ++i)
        for (std::size_t c = 0; c < R[i].size() && c < out[i].size(); ++c) {
            out[i][c].t = R[i][c].t;
            out[i][c].F = R[i][c].F;
        }
    return out;
}
 
Matrix<double> SolverCTMC::state_space_aggr() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const ctmc::CtmcSolution<double> d = ctmc::solver_ctmc_analyzer(sn, ctmc_options(opts_));
    return ctmc::ctmc_get_state_space_aggr(sn, d);
}
 
std::size_t StateView::column(std::size_t node, std::size_t cls) const {
    if (node == 0 || node > sn_->nodes.size())
        throw InputError("StateView: node index " + std::to_string(node) + " does not exist");
    const std::size_t ist = sn_->nodes[node - 1].station;
    if (ist == 0)
        throw InputError("StateView: node '" + sn_->nodes[node - 1].name +
                         "' is not a station, so the aggregate state holds no count for it");
    if (cls == 0 || cls > sn_->nclasses)
        throw InputError("StateView: class index " + std::to_string(cls) + " does not exist");
    return (ist - 1) * sn_->nclasses + (cls - 1);
}
 
double StateView::at(std::size_t node, std::size_t cls) const {
    return (*row_)[column(node, cls)];
}
 
double StateView::total(std::size_t node) const {
    double tot = 0.0;
    for (std::size_t k = 1; k <= sn_->nclasses; ++k) tot += (*row_)[column(node, k)];
    return tot;
}
 
namespace {
 
/**
 * The 1-based rows on which `fn` holds, as the n x 1 index matrix the CTMC
 * getters resolve; empty for an empty predicate or when no state satisfies it,
 * as the reference's `ctmc_predicate_states` returns `[]`.
 */
Matrix<double> predicate_states(const StatePredicate& fn, const qn::NetworkStruct<double>& sn,
                                 const Matrix<double>& aggr, std::size_t n, const char* name) {
    if (!fn) return Matrix<double>();
    if (aggr.rows() != n)
        throw InputError(std::string("The aggregated state space has ") +
                         std::to_string(aggr.rows()) + " rows where the generator has " +
                         std::to_string(n) + ", so the predicate for set " + name +
                         " cannot be resolved state by state.");
    std::vector<double> hit;
    std::vector<double> row(aggr.cols());
    for (std::size_t s = 0; s < n; ++s) {
        for (std::size_t c = 0; c < aggr.cols(); ++c) row[c] = aggr(s, c);
        if (fn(StateView(row, sn))) hit.push_back(static_cast<double>(s + 1));
    }
    Matrix<double> out(hit.size(), hit.empty() ? 0 : 1);
    for (std::size_t i = 0; i < hit.size(); ++i) out(i, 0) = hit[i];
    return out;
}
 
FirstPassage to_first_passage(const ctmc::CtmcFirstPassage& r) {
    FirstPassage out;
    out.t = r.t;
    out.F = r.F;
    out.f = r.f;
    out.reachprob = r.reachprob;
    out.tabooprob = r.tabooprob;
    out.source = r.source;
    out.target = r.target;
    out.taboo = r.taboo;
    return out;
}
 
FirstPassageMoments to_first_passage_moments(const ctmc::CtmcFirstPassageMoments<double>& r) {
    FirstPassageMoments out;
    out.m = r.m;
    out.mall = r.mall;
    out.p = r.p;
    out.source = r.source;
    out.target = r.target;
    out.taboo = r.taboo;
    return out;
}
 
}  // namespace
 
FirstPassage SolverCTMC::cdf_first_passt(const Matrix<double>& A, const Matrix<double>& B,
                                         const Matrix<double>& C,
                                         const std::vector<double>& tset) {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const ctmc::CtmcSolution<double> d = ctmc::solver_ctmc_analyzer(sn, ctmc_options(opts_));
    return to_first_passage(ctmc::ctmc_cdf_firstpasst(sn, d, A, B, C, tset));
}
 
FirstPassage SolverCTMC::cdf_first_passt(const StatePredicate& A, const StatePredicate& B,
                                         const StatePredicate& C,
                                         const std::vector<double>& tset) {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const ctmc::CtmcSolution<double> d = ctmc::solver_ctmc_analyzer(sn, ctmc_options(opts_));
    const qn::NetworkStruct<double>& ssn = d.solved_sn ? *d.solved_sn : sn;
    const Matrix<double> aggr = ctmc::ctmc_get_state_space_aggr(sn, d);
    const std::size_t n = d.chain.Q.rows();
    return to_first_passage(ctmc::ctmc_cdf_firstpasst(
        sn, d, predicate_states(A, ssn, aggr, n, "A"), predicate_states(B, ssn, aggr, n, "B"),
        predicate_states(C, ssn, aggr, n, "C"), tset));
}
 
FirstPassageMoments SolverCTMC::first_passt_moments(const Matrix<double>& A,
                                                    const Matrix<double>& B, std::size_t nmax,
                                                    const Matrix<double>& C) {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const ctmc::CtmcSolution<double> d = ctmc::solver_ctmc_analyzer(sn, ctmc_options(opts_));
    return to_first_passage_moments(ctmc::ctmc_firstpasst_moments(sn, d, A, B, nmax, C));
}
 
FirstPassageMoments SolverCTMC::first_passt_moments(const StatePredicate& A,
                                                    const StatePredicate& B, std::size_t nmax,
                                                    const StatePredicate& C) {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const ctmc::CtmcSolution<double> d = ctmc::solver_ctmc_analyzer(sn, ctmc_options(opts_));
    const qn::NetworkStruct<double>& ssn = d.solved_sn ? *d.solved_sn : sn;
    const Matrix<double> aggr = ctmc::ctmc_get_state_space_aggr(sn, d);
    const std::size_t n = d.chain.Q.rows();
    return to_first_passage_moments(ctmc::ctmc_firstpasst_moments(
        sn, d, predicate_states(A, ssn, aggr, n, "A"), predicate_states(B, ssn, aggr, n, "B"),
        nmax, predicate_states(C, ssn, aggr, n, "C")));
}
 
SymbolicSolution SolverCTMC::symbolic_solution() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const ctmc::CtmcSymbolicSolution<double> r =
        ctmc::ctmc_symbolic_solution(sn, ctmc_options(opts_), ctmc_symbolic_options(opts_));
    SymbolicSolution out;
    out.pi = r.pi;
    out.num = r.num;
    out.den = r.den;
    out.engine = r.engine;
    out.symbols = r.symbols;
    out.rate0 = r.rate0;
    return out;
}
 
void SolverCTMC::print_inf_gen(const Matrix<double>& Q, const Matrix<double>& space) {
    std::printf("%8s  %-28s %s\n", "State", "Marking", "Rates (to: rate)");
    for (std::size_t i = 0; i < Q.rows(); ++i) {
        std::string mark;
        if (i < space.rows()) {
            std::ostringstream os;
            os << "[";
            for (std::size_t c = 0; c < space.cols(); ++c) os << (c ? " " : "") << space(i, c);
            os << "]";
            mark = os.str();
        }
        std::printf("%8zu  %-28s", i + 1, mark.c_str());
        for (std::size_t j = 0; j < Q.cols(); ++j)
            if (i != j && Q(i, j) != 0.0) std::printf("  %zu: %.6g", j + 1, Q(i, j));
        std::printf("\n");
    }
}
 
// ---------------------------------------------------------------------------
// SolverFLD
// ---------------------------------------------------------------------------
 
TranAvg SolverFLD::tran_avg() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const std::vector<fluid::FluidTranPoint> pts =
        fluid::solver_fluid_tran_avg(sn, fluid_options(opts_), 100);
    TranAvg out;
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (std::size_t c = 0; c < sn.nclasses; ++c)
            out.label.push_back(sn.stations[i].name + "/" + sn.classes[c].name);
    for (std::size_t s = 0; s < pts.size(); ++s) {
        out.t.push_back(pts[s].t);
        std::vector<double> row;
        for (std::size_t i = 0; i < sn.nstations; ++i)
            for (std::size_t c = 0; c < sn.nclasses; ++c) row.push_back(pts[s].QN(i, c));
        out.QNt.push_back(row);
    }
    return out;
}
 
std::vector<std::vector<CdfCurve> > SolverFLD::cdf_respt() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    std::vector<std::vector<CdfCurve> > out(sn.nstations, std::vector<CdfCurve>(sn.nclasses));
    const std::vector<std::vector<fluid::FluidPassage> > R =
        fluid::solver_fluid_cdf_respt(sn, fluid_options(opts_));
    for (std::size_t i = 0; i < R.size() && i < out.size(); ++i)
        for (std::size_t c = 0; c < R[i].size() && c < out[i].size(); ++c) {
            out[i][c].t = R[i][c].t;
            out[i][c].F = R[i][c].cdf;
        }
    return out;
}
 
// ---------------------------------------------------------------------------
// SolverLDES
// ---------------------------------------------------------------------------
 
std::vector<std::vector<CdfCurve> > SolverLDES::cdf_respt() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const ldes::LdesResult r =
        ldes::solver_ldes(sn, ldes_options(opts_), std::vector<std::string>(1, "--respt-samples"));
    if (r.respTimeSamples.empty())
        throw ExternalError(ErrorKind{"LDES"},
                            "SolverLDES::cdf_respt needs the per-job response times the engine records under "
                            "--respt-samples and the run returned none; raise the samples so completions are observed");
    std::vector<std::vector<CdfCurve> > out(sn.nstations, std::vector<CdfCurve>(sn.nclasses));
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (std::size_t c = 0; c < sn.nclasses; ++c) {
            const Matrix<double> m = ldes::ldes_cdf_respt(r, i, c);
            for (std::size_t j = 0; j < m.rows(); ++j) {
                out[i][c].F.push_back(m(j, 0));
                out[i][c].t.push_back(m(j, 1));
            }
        }
    return out;
}
 
SolverLDES& SolverLDES::init_from_solver(NetworkSolver& aux) {
    const qn::NetworkStruct<double> sn = model_->get_struct();
    const std::size_t M = sn.nstations, K = sn.nclasses;
    Matrix<double> P;
    if (aux.get_name() == "CTMC") {
        P = ldes::warm_start_placement_from_ctmc(sn, ctmc_options(aux.options()));
    } else {
        // The table drops all-zero rows, so the means are matched back by name.
        const AvgTable& t = aux.avg_table();
        Matrix<double> QN(M, K, 0.0);
        for (std::size_t row = 0; row < t.QLen.size(); ++row)
            for (std::size_t i = 0; i < M; ++i)
                for (std::size_t c = 0; c < K; ++c)
                    if (sn.stations[i].name == t.Station[row] && sn.classes[c].name == t.JobClass[row])
                        QN(i, c) = t.QLen[row];
        P = ldes::warm_start_placement_from_qlen(sn, QN);
    }
    // `initFromMarginal`: every state of each station with that many jobs, the first one
    // weighed, declared as the model initial state every state-driven solver starts from.
    for (std::size_t i = 1; i <= M; ++i) {
        if (sn.stations[i - 1].nodetype == lang::NodeType::Source) continue;
        std::vector<std::size_t> n(K, 0), ph(K, 1);
        for (std::size_t c = 0; c < K; ++c) {
            n[c] = static_cast<std::size_t>(P(i - 1, c) + 0.5);
            ph[c] = sn.phases_of(i, c + 1);
        }
        const std::vector<std::vector<double> > rows = qn::from_marginal(sn, i, n, ph);
        if (rows.empty())
            throw InputError("SolverLDES::init_from_solver: station '" + sn.stations[i - 1].name +
                             "' admits no state with the warm-start placement");
        Matrix<double> space(rows.size(), rows[0].size());
        for (std::size_t r = 0; r < rows.size(); ++r)
            for (std::size_t c = 0; c < rows[r].size(); ++c) space(r, c) = rows[r][c];
        std::vector<double> prior(rows.size(), 0.0);
        prior[0] = 1.0;
        model_->set_state_prior(sn.node_of_station(i), space, prior);
    }
    opts_.tranfilter = "fixed";
    opts_.warmupfrac = 0.0;
    reset();
    jobs_ready_ = false;
    return *this;
}
 
SolverLDES& SolverLDES::init_from_jobs(const Matrix<double>& jobs) {
    if (jobs.rows() == 0 || jobs.cols() < 7)
        throw InputError("SolverLDES::init_from_jobs needs one row per job with the columns [station, class, "
                         "remaining, elapsed, qage, sage, order, id], id optional");
    opts_.init_jobs = jobs;
    opts_.track_jobs = true;
    reset();
    jobs_ready_ = false;
    return *this;
}
 
void SolverLDES::run_tracked() {
    if (jobs_ready_ && opts_.cache) return;
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    ldes::LdesOptions opt = ldes_options(opts_);
    opt.track_jobs = true;
    const ldes::LdesResult r = ldes::solver_ldes(sn, opt);
    final_jobs_ = r.finalJobs;
    job_departures_ = r.jobDepartures;
    jobs_ready_ = true;
}
 
const Matrix<double>& SolverLDES::final_jobs() {
    run_tracked();
    return final_jobs_;
}
 
const Matrix<double>& SolverLDES::job_departures() {
    run_tracked();
    return job_departures_;
}
 
// ---------------------------------------------------------------------------
// SolverJMT
// ---------------------------------------------------------------------------
 
std::vector<std::vector<CdfCurve> > SolverJMT::cdf_respt() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    std::vector<std::vector<CdfCurve> > out(sn.nstations, std::vector<CdfCurve>(sn.nclasses));
    // `true` is the getCdfRespT contract: the model is first solved for its
    // steady-state queue lengths and the logged run starts preloaded at their
    // rounded values, which is what makes the seeded curve comparable to the
    // reference's for the same seed.
    const std::map<std::pair<std::size_t, std::size_t>, std::vector<std::pair<double, double> > >
        rd = jmt::jmt_get_cdf_resp_t(sn, jmt_options(opts_), true);
    for (std::map<std::pair<std::size_t, std::size_t>,
                  std::vector<std::pair<double, double> > >::const_iterator it = rd.begin();
         it != rd.end(); ++it) {
        const std::size_t i = it->first.first, c = it->first.second;
        if (i == 0 || i > out.size() || c == 0 || c > sn.nclasses) continue;
        // THE PAIRS ARE (F, X), as `ecdf` returns them, and the curve stores
        // the two apart: reading them the other way round would report a
        // probability as a time.
        for (std::size_t j = 0; j < it->second.size(); ++j) {
            out[i - 1][c - 1].F.push_back(it->second[j].first);
            out[i - 1][c - 1].t.push_back(it->second[j].second);
        }
    }
    return out;
}
 
TranAvg SolverJMT::tran_avg() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const jmt::JmtReplication<double> r = jmt::jmt_transient_replications(sn, jmt_options(opts_));
    TranAvg out;
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (std::size_t c = 0; c < sn.nclasses; ++c)
            out.label.push_back(sn.stations[i].name + "/" + sn.classes[c].name);
    out.t = r.t;
    for (std::size_t g = 0; g < r.t.size(); ++g) {
        std::vector<double> row;
        for (std::size_t i = 0; i < sn.nstations; ++i)
            for (std::size_t c = 0; c < sn.nclasses; ++c)
                row.push_back(i < r.QNt.size() && c < r.QNt[i].size() && g < r.QNt[i][c].size()
                                  ? r.QNt[i][c][g]
                                  : 0.0);
        out.QNt.push_back(row);
    }
    return out;
}
 
double SolverJMT::prob_aggr(std::size_t node, const std::vector<double>& state) {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    if (node == 0 || node > sn.nodes.size()) throw InputError("prob_aggr: node index out of range");
    const std::size_t ist = sn.nodes[node - 1].station;
    if (ist == 0) throw InputError("prob_aggr: the node is not a station");
    const jmt::JmtProbAggr r = jmt::jmt_prob_aggr(sn, jmt_options(opts_), ist, state);
    return r.station[ist - 1];
}
 
SamplePath SolverJMT::sample(std::size_t events, std::size_t node) {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    const jmt::JmtOptions opt = jmt_options(opts_);
    SamplePath out;
    if (node) {
        if (node > sn.nodes.size() || sn.nodes[node - 1].station == 0)
            throw InputError("sample: the node is not a station, so nothing is logged for it");
        const jmt::JmtNodeTrace<double> tr = jmt::jmt_sample_aggr(sn, node, events, opt);
        for (std::size_t c = 0; c < sn.nclasses; ++c) out.label.push_back(sn.classes[c].name);
        out.t = tr.t;
        out.state = tr.qlen;
        return out;
    }
    const jmt::JmtSysTrace<double> tr = jmt::jmt_sample_sys_aggr(sn, events, opt);
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (std::size_t c = 0; c < sn.nclasses; ++c)
            out.label.push_back(sn.stations[i].name + "/" + sn.classes[c].name);
    out.t = tr.t;
    for (std::size_t g = 0; g < tr.t.size(); ++g) {
        std::vector<double> row;
        for (std::size_t i = 0; i < sn.nstations; ++i)
            for (std::size_t c = 0; c < sn.nclasses; ++c)
                row.push_back(i < tr.state.size() && g < tr.state[i].size() &&
                                      c < tr.state[i][g].size()
                                  ? tr.state[i][g][c]
                                  : 0.0);
        out.state.push_back(row);
    }
    return out;
}
 
// ---------------------------------------------------------------------------
// SolverBA
// ---------------------------------------------------------------------------
 
BoundsTable SolverBA::bounds_table() {
    const qn::NetworkStruct<double>& sn = model_->get_struct();
    ba::BaOptions opt;
    if (!opts_.method.empty()) opt.method = opts_.method;
    opt.level = opts_.level;
    const ba::BaBounds<double> b = ba::ba_bounds(sn, opt);
    BoundsTable t;
    t.method = opt.method;
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (std::size_t c = 0; c < sn.nclasses; ++c) {
            if (i < b.keep.size() && c < b.keep[i].size() && !b.keep[i][c]) continue;
            t.Station.push_back(sn.stations[i].name);
            t.JobClass.push_back(sn.classes[c].name);
            t.Qlower.push_back(b.Qlower(i, c));
            t.Qupper.push_back(b.Qupper(i, c));
            t.Tlower.push_back(b.Tlower(i, c));
            t.Tupper.push_back(b.Tupper(i, c));
        }
    return t;
}
 
// ---------------------------------------------------------------------------
// SolverAUTO
// ---------------------------------------------------------------------------
 
namespace {
 
/** A ranking slot's facade name: "fld" is SolverFLD, the rest are their upper-case selves. */
std::string facade_name_of_slot(autosolver::AutoSolver s) {
    std::string name = autosolver::auto_solver_name(s);
    std::transform(name.begin(), name.end(), name.begin(), ::toupper);
    return name;
}
 
/** The facade name's method family, the inverse of `facade_name_of_family`. */
std::string family_of_facade_name(const std::string& tag) {
    std::string fam = tag;
    std::transform(fam.begin(), fam.end(), fam.begin(), ::tolower);
    return fam;
}
 
/**
 * Who answers `getter` under the configured method, in `delegate.m`'s order. A family names its
 * engine outright and is the whole list; an intent names a ranking, whose first slot alone is asked
 * with the method the ranking gated on ('exact' under the exact intent).
 */
struct AutoDelegatePlan {
    bool family = false;
    std::string family_name;
    std::vector<std::string> tags;
    std::string first_method;
};
 
AutoDelegatePlan auto_delegate_plan(const qn::NetworkStruct<double>& sn, const std::string& method,
                                    const std::string& getter) {
    AutoDelegatePlan plan;
    const autosolver::AutoToken tok = autosolver::auto_resolve_token(method, sn);
    if (!tok.is_intent) {
        plan.family = true;
        plan.family_name = tok.family;
        plan.tags.push_back(facade_name_of_family(tok.family));
        plan.first_method = tok.submethod;
        return plan;
    }
    const autosolver::AutoChoice c = autosolver::auto_choose_solver_mode(sn, getter, tok.mode);
    const std::vector<autosolver::AutoSolver> order = autosolver::auto_proposed_solvers(sn, getter, tok.mode);
    for (std::size_t i = 0; i < order.size(); ++i) plan.tags.push_back(facade_name_of_slot(order[i]));
    plan.first_method = c.method.empty() ? std::string("default") : c.method;
    return plan;
}
 
/**
 * `delegate.m` for a getter the facade answers on some solvers only. `answers` names the facades
 * that expose the call; a proposed solver outside it is passed over, since this port has nothing to
 * ask it. Input and internal errors propagate; any other moves on to the next candidate, and when
 * all fail the LAST error is reported under its own category, as the reference does.
 */
template <class R, class Call>
R auto_delegate(Network& model, const SolverOptions& o, const std::string& getter,
                const std::set<std::string>& answers, Call call, std::string& last) {
    const qn::NetworkStruct<double>& sn = model.get_struct();
    const AutoDelegatePlan plan = auto_delegate_plan(sn, o.method, getter);
    std::string last_error, last_cat, last_kind, proposed;
    bool tried = false;
    for (std::size_t i = 0; i < plan.tags.size(); ++i) {
        const std::string& tag = plan.tags[i];
        proposed += (proposed.empty() ? "" : ", ") + tag;
        if (!answers.count(tag)) continue;
        tried = true;
        SolverOptions oo = o;
        oo.method = i == 0 ? plan.first_method : std::string("default");
        try {
            R r = call(tag, model, oo);
            last = tag;
            return r;
        } catch (const Error& e) {
            const std::string cat = e.category();
            if (cat == "input" || cat == "internal") throw;
            last_error = tag + ": " + e.what();
            last_cat = cat;
            last_kind = e.kind();
        }
    }
    if (!tried) {
        if (plan.family)
            throw UnsupportedError(ErrorKind{"Method"},
                "SolverAUTO: method family '" + plan.family_name + "' does not answer " + getter +
                " in this port");
        throw UnsupportedError(ErrorKind{"Method"},
            "SolverAUTO: none of the solvers proposed for " + getter + " (" + proposed +
            ") answers it in this port");
    }
    const std::string msg = "SolverAUTO: All candidate solvers failed to execute method '" + getter +
                            "'. Last error: " + last_error;
    if (last_cat == "numerical") throw NumericError(ErrorKind{last_kind}, msg);
    if (last_cat == "external") throw ExternalError(ErrorKind{last_kind}, msg);
    throw UnsupportedError(ErrorKind{last_kind.empty() ? std::string("Solver") : last_kind}, msg);
}
 
/** True for a facade whose answer is an estimate from a simulated run under this method. */
bool facade_is_stochastic(const std::string& tag, const std::string& method) {
    if (tag == "SSA" || tag == "LDES") return true;
    if (tag == "JMT") return method.compare(0, 4, "jmva") != 0;
    return false;
}
 
}  // namespace
 
SolverAUTO::SolverAUTO(Network& m, const SolverOptions& o) : NetworkSolver(m, "AUTO", o) {
    // The constructor resolves the method as the reference's does, so an unknown name fails here
    // and not at the first getter. A bare method name needs the model to resolve, and is checked
    // against it only when the model-free resolution cannot place it.
    try {
        autosolver::auto_resolve_token(o.method);
    } catch (const InputError&) {
        autosolver::auto_resolve_token(o.method, m.get_struct());
    }
}
 
SolverAUTO::SolverAUTO(Network& m, const std::string& method)
    : SolverAUTO(m, SolverOptions().set_method(method)) {}
 
std::string SolverAUTO::selected_solver_name() {
    if (!last_.empty()) return last_;
    if (solved_) return table_.solver;
    try {
        return choose_solver("getAvg");
    } catch (const Error&) {
        // No feasible delegate: the caller learns this from the analysis itself.
        return std::string();
    }
}
 
std::string SolverAUTO::selected_method() {
    const std::string tag = selected_solver_name();
    return tag.empty() ? std::string() : family_of_facade_name(tag);
}
 
std::vector<std::string> SolverAUTO::candidate_solver_names() const {
    std::vector<std::string> out;
    const std::vector<autosolver::AutoSolver> c = autosolver::auto_candidates(model_->get_struct());
    for (std::size_t i = 0; i < c.size(); ++i) out.push_back(facade_name_of_slot(c[i]));
    return out;
}
 
std::string SolverAUTO::choose_solver(const std::string& getter) const {
    return auto_delegate_plan(model_->get_struct(), opts_.method, getter).tags.front();
}
 
std::vector<std::string> SolverAUTO::proposed_solver_names(const std::string& getter) const {
    return auto_delegate_plan(model_->get_struct(), opts_.method, getter).tags;
}
 
bool SolverAUTO::is_stochastic() {
    if (!last_.empty()) return facade_is_stochastic(last_, std::string());
    if (solved_) return facade_is_stochastic(table_.solver, table_.method);
    const std::vector<std::string> c = candidate_solver_names();
    for (std::size_t i = 0; i < c.size(); ++i)
        if (facade_is_stochastic(c[i], std::string())) return true;
    return false;
}
 
std::vector<std::string> SolverAUTO::list_all_methods() const {
    std::vector<std::string> out = selection_intents();
    out.push_back("auto");
    const std::vector<std::string> fams = family_names();
    for (std::size_t f = 0; f < fams.size(); ++f) {
        // A family contributes only when it can be instantiated on a Network here: ln, env and uq
        // take another model class, ldes needs its engine and lqns its qnsolver binary.
        std::string tag;
        try {
            tag = facade_name_of_family(fams[f]);
        } catch (const Error&) {
            continue;
        }
        if (tag == "LQNS" && !lqns::qnsolver_is_available()) continue;
        out.push_back(fams[f]);
        std::string key = fams[f];
        std::transform(key.begin(), key.end(), key.begin(), ::toupper);
        const std::vector<std::string> declared = method_registry::methods(key);
        for (std::size_t m = 0; m < declared.size(); ++m) out.push_back(fams[f] + "." + declared[m]);
    }
    std::sort(out.begin(), out.end());
    out.erase(std::unique(out.begin(), out.end()), out.end());
    return out;
}
 
std::vector<AutoCandidate> SolverAUTO::find_solver(const std::string& metric, bool show_all) const {
    const std::vector<autosolver::SolverCandidate> rows =
        autosolver::auto_find_solver(model_->get_struct(), metric, show_all);
    std::vector<AutoCandidate> out;
    for (std::size_t i = 0; i < rows.size(); ++i) {
        AutoCandidate c;
        c.solver = rows[i].solver;
        c.method = rows[i].method;
        c.runnable = rows[i].runnable;
        c.method_class = rows[i].method_class;
        c.metrics = rows[i].metrics;
        c.reason = rows[i].reason;
        out.push_back(c);
    }
    return out;
}
 
std::vector<std::vector<CdfCurve> > SolverAUTO::cdf_respt() {
    std::set<std::string> answers;
    answers.insert("CTMC");
    answers.insert("FLD");
    answers.insert("LDES");
    answers.insert("JMT");
    return auto_delegate<std::vector<std::vector<CdfCurve> > >(
        *model_, opts_, "getCdfRespT", answers,
        [](const std::string& tag, Network& m, const SolverOptions& oo) {
            if (tag == "CTMC") return SolverCTMC(m, oo).cdf_respt();
            if (tag == "FLD") return SolverFLD(m, oo).cdf_respt();
            if (tag == "LDES") return SolverLDES(m, oo).cdf_respt();
            return SolverJMT(m, oo).cdf_respt();
        },
        last_);
}
 
TranAvg SolverAUTO::tran_avg() {
    std::set<std::string> answers;
    answers.insert("FLD");
    answers.insert("JMT");
    return auto_delegate<TranAvg>(
        *model_, opts_, "getTranAvg", answers,
        [](const std::string& tag, Network& m, const SolverOptions& oo) {
            if (tag == "FLD") return SolverFLD(m, oo).tran_avg();
            return SolverJMT(m, oo).tran_avg();
        },
        last_);
}
 
double SolverAUTO::prob_aggr(std::size_t node, const std::vector<double>& state) {
    std::set<std::string> answers;
    answers.insert("CTMC");
    answers.insert("JMT");
    return auto_delegate<double>(
        *model_, opts_, "getProbAggr", answers,
        [node, &state](const std::string& tag, Network& m, const SolverOptions& oo) {
            if (tag == "CTMC") return SolverCTMC(m, oo).prob_aggr(node, state);
            return SolverJMT(m, oo).prob_aggr(node, state);
        },
        last_);
}
 
namespace {
 
/**
 * The options a `ctmcSolver()`-style accessor builds its solver with: the family's submethod when
 * the method names that family, its default otherwise (an intent names no CTMC method).
 */
SolverOptions family_solver_options(const SolverOptions& o, const std::string& family,
                                    const std::string& fallback) {
    SolverOptions oo = o;
    const std::size_t dot = o.method.find('.');
    const std::string head = dot == std::string::npos ? o.method : o.method.substr(0, dot);
    if (autosolver::auto_family_alias(head) == family && dot != std::string::npos)
        oo.method = o.method.substr(dot + 1);
    else
        oo.method = fallback;
    return oo;
}
 
}  // namespace
 
SolverCTMC SolverAUTO::ctmc_solver() const {
    return SolverCTMC(*model_, family_solver_options(opts_, "ctmc", "default"));
}
 
SolverFLD SolverAUTO::fld_solver() const {
    return SolverFLD(*model_, family_solver_options(opts_, "fld", "default"));
}
 
SolverNC SolverAUTO::nc_solver() const {
    return SolverNC(*model_, family_solver_options(opts_, "nc", "default"));
}
 
SolverBA SolverAUTO::ba_solver() const {
    // baSolver(): a bound method named outright is kept, every other name asks SolverBA for 'auto'.
    return SolverBA(*model_, family_solver_options(opts_, "ba", "auto"));
}
 
std::vector<std::string> SolverAUTO::selection_intents() {
    return {"default", "heur", "sim", "exact", "fast", "accurate", "bound"};
}
 
std::vector<std::string> SolverAUTO::family_names() {
    return {"mva", "nc", "ctmc", "fld", "mam", "ag", "ba", "ssa", "ldes", "jmt", "ln", "env", "lqns", "uq"};
}
 
std::vector<std::string> SolverAUTO::metric_groups() { return autosolver::auto_metric_groups(); }
 
std::string SolverAUTO::metric_group_of(const std::string& name) {
    return autosolver::auto_metric_group_of(name);
}
 
std::vector<std::string> SolverAUTO::family_metrics(const std::string& family) {
    return autosolver::auto_family_metrics(family);
}
 
namespace {
 
qn::FeatureSet auto_feature_union() {
    const qn::FeatureSet sets[] = {qn::mva_feature_set("default"),   qn::nc_feature_set("default"),
                                   qn::ctmc_feature_set("default"),  qn::fluid_class_feature_set(),
                                   qn::mam_feature_set("default"),   qn::ssa_feature_set("default"),
                                   qn::ldes_feature_set("default"),  qn::jmt_feature_set("default")};
    qn::FeatureSet out;
    for (std::size_t s = 0; s < sizeof(sets) / sizeof(*sets); ++s) {
        const std::vector<qn::Feature> fs = sets[s].list();
        for (std::size_t i = 0; i < fs.size(); ++i) out.set(fs[i]);
    }
    return out;
}
 
}  // namespace
 
std::vector<std::string> SolverAUTO::feature_set() {
    const std::vector<qn::Feature> fs = auto_feature_union().list();
    std::vector<std::string> out;
    for (std::size_t i = 0; i < fs.size(); ++i) out.push_back(qn::feature_name(fs[i]));
    return out;
}
 
bool SolverAUTO::supports(Network& model) {
    return qn::feature_set_supports("SolverAUTO", auto_feature_union(),
                                    qn::used_lang_features(model.get_struct()))
        .ok;
}
 
}  // namespace line