api_dispatch_misc.cpp

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 *
 * --api dispatch entries for the 'misc' group (see api_dispatch_detail.h):
 * polling, flow-equivalent servers (fes), caches, cache decomposition (da),
 * loss networks, retrieval, sum, MAP queueing-network bounds (mapqn),
 * workflow detectors, decision diagrams (mdd), stochastic Petri nets (spn)
 * and a handful of mam / pfqn / trace / dqsys helpers.
 *
 * CONVENTIONS:
 *  - argument keys are the MATLAB parameter names and result keys the MATLAB
 *    output names; a MATLAB struct argument or output (options, params, ch,
 *    mde, info, stats) is a nested object keyed by its MATLAB field names;
 *  - a MAP is the MATLAB cell {D0, D1}, written as a two-element list
 *    [D0, D1]; a list of MAPs is a list of such pairs, and in a LIST every
 *    matrix is written as an array of rows so the two shapes cannot be
 *    confused (a single MAP may write a 1x1 matrix as a bare number);
 *  - a station, item, list, node or place INDEX passes through in the C++
 *    port's own convention, which the port headers document: 0-based unless
 *    the header says otherwise (ljd_linearize and the mdd node ids are
 *    1-based because the reference's are). An index is not a measurement and
 *    is never shifted on the way out;
 *  - an argument MATLAB accepts but the port has no parameter for is refused
 *    by name (Args::unsupported) rather than dropped;
 *  - an absent optional argument calls the port's shorter overload, or leaves
 *    the port's options-struct default in place; no default is restated here.
 */
 
#include "api_dispatch_detail.h"
 
#include <cctype>
#include <cstdint>
#include <cstdlib>
#include <string>
#include <vector>
 
#include "line/api/mam/map_moment.h"
#include "line/api/polling/polling_qsys_1limited.h"
#include "line/api/polling/polling_qsys_exhaustive.h"
#include "line/api/polling/polling_qsys_klimited.h"
#include "line/api/polling/polling_qsys_table.h"
#include "line/lang/lang_types.h"
#include "line/api/fes/ljd_linearize.h"
#include "line/api/fes/fes_compute_throughputs.h"
#include "line/api/fes/fes_compute_metrics.h"
#include "line/api/fes/fes_build_isolated.h"
#include "line/api/fes/fes_map_interdeparture.h"
#include "line/api/fes/fes_map_interp.h"
#include "line/api/fes/fes_map_levels.h"
#include "line/api/fes/fes_map_moments.h"
#include "line/api/fes/fes_map_solve.h"
#include "line/api/fes/fes_map_aggregate.h"
#include "line/api/cache/cache_erec.h"
#include "line/api/cache/cache_cost.h"
#include "line/api/cache/cache_spm_size.h"
#include "line/api/cache/cache_gamma.h"
#include "line/api/cache/cache_gamma_lp.h"
#include "line/api/cache/cache_miss_asy.h"
#include "line/api/cache/cache_miss.h"
#include "line/api/cache/cache_miss_rmf.h"
#include "line/api/cache/cache_rrm_meanfield.h"
#include "line/api/cache/cache_rrm_meanfield_ode.h"
#include "line/api/cache/cache_miss_fpi.h"
#include "line/api/cache/cache_miss_spm.h"
#include "line/api/cache/cache_mva.h"
#include "line/api/cache/cache_mva_miss.h"
#include "line/api/cache/cache_prob_erec.h"
#include "line/api/cache/cache_prob_fpi.h"
#include "line/api/cache/cache_prob_spm.h"
#include "line/api/cache/cache_spm.h"
#include "line/api/cache/cache_t_hlru.h"
#include "line/api/cache/cache_ttl_hlru.h"
#include "line/api/cache/cache_xi_fp.h"
#include "line/api/cache/cache_xi_iter.h"
#include "line/api/cache/cache_is.h"
#include "line/api/cache/cache_miss_is.h"
#include "line/api/cache/cache_prob_is.h"
#include "line/api/cache/cache_lrum_map_levelstats.h"
#include "line/api/cache/cache_t_lrum_map.h"
#include "line/api/cache/cache_ttl_lrua.h"
#include "line/api/cache/cache_ttl_lrum_map.h"
#include "line/api/cache/cache_sojourn.h"
#include "line/api/da/da_cache_isolate.h"
#include "line/api/da/da_traffic_superpos.h"
#include "line/api/lossn/lossn_erlangfp.h"
#include "line/api/lossn/lossn_manjunath.h"
#include "line/api/lossn/lossn_mci.h"
#include "line/api/lossn/lossn_rec.h"
#include "line/api/lossn/lossn_asy.h"
#include "line/api/lossn/lossn_kr.h"
#include "line/api/lossn/lossn_links.h"
#include "line/api/lossn/lossn_revsens.h"
#include "line/api/dqsys/dqsys_bernoulli1.h"
#include "line/api/mam/map_transform.h"
#include "line/api/mam/map2mmpp.h"
#include "line/api/trace/autocov.h"
#include "line/api/mam/aph_fit.h"
#include "line/api/mam/m3pp22_fitc_cov.h"
#include "line/api/mam/ldqbd_mphc.h"
#include "line/api/mam/mmdp_isfeasible.h"
#include "line/api/pfqn/infradius_h.h"
#include "line/api/pfqn/pas_placement.h"
#include "line/api/mapqn/mapqn_amva.h"
#include "line/api/mapqn/mapqn_bnd_lr.h"
#include "line/api/mapqn/mapqn_bnd_lr_mva.h"
#include "line/api/mapqn/mapqn_bnd_qr.h"
#include "line/api/mapqn/mapqn_bnd_qr_ld.h"
#include "line/api/mapqn/mapqn_bnd_qr_delay.h"
#include "line/api/mapqn/mapqn_qr_bounds_bas.h"
#include "line/api/mapqn/mapqn_qr_bounds_rsrd.h"
#include "line/api/retrieval/retrieval_metrics.h"
#include "line/api/retrieval/retrieval_mva.h"
#include "line/api/retrieval/retrieval_fpi.h"
#include "line/api/retrieval/retrieval_fpi_latency.h"
#include "line/api/retrieval/retrieval_nc.h"
#include "line/api/retrieval/retrieval_rayint.h"
#include "line/api/sum/sum_closed.h"
#include "line/api/sum/sum_closing.h"
#include "line/api/wf/wf_branch_detector.h"
#include "line/api/wf/wf_loop_detector.h"
#include "line/api/wf/wf_parallel_detector.h"
#include "line/api/wf/wf_sequence_detector.h"
#include "line/api/mdd/mdd.h"
#include "line/api/mdd/mdd_rec.h"
#include "line/api/mdd/mdd_closedqn.h"
#include "line/api/spn/spn_conv.h"
#include "line/api/spn/spn_ispn.h"
#include "line/api/spn/spn_mdd.h"
#include "line/api/spn/spn_rec_enabled.h"
#include "line/api/spn/spn_sinvariants.h"
 
namespace line {
namespace reg {
namespace {
 
// ---------------------------------------------------------------------------
// Shared readers and encoders
// ---------------------------------------------------------------------------
 
std::string at(const std::string& where, std::size_t k) {
    return where + "[" + std::to_string(k) + "]";
}
 
/** One non-negative index or count from a JSON number or numeric string. */
std::size_t size_from_json(const Json& j, const std::string& where) {
    const std::string text = decimal_text(j, where);
    char* end = nullptr;
    const long long v = std::strtoll(text.c_str(), &end, 10);
    if (end == text.c_str() || *end != '\0')
        throw InputError(where + ": expected a non-negative integer, got " + text);
    if (v < 0) throw InputError(where + ": expected a non-negative integer, got " + text);
    return static_cast<std::size_t>(v);
}
 
/** A flat list of non-negative integers; a bare number is a one-element list. */
std::vector<std::size_t> sizes_from_json(const Json& j, const std::string& where) {
    std::vector<std::size_t> out;
    if (!j.is_array()) {
        out.push_back(size_from_json(j, where));
        return out;
    }
    for (std::size_t k = 0; k < j.size(); ++k) out.push_back(size_from_json(j[k], at(where, k)));
    return out;
}
 
Json encode_sizes(const std::vector<std::size_t>& v) {
    Json a = Json::array();
    for (std::size_t x : v) a.push_back(static_cast<std::uint64_t>(x));
    return a;
}
 
/** A flag vector: booleans or 0/1, a bare value being a one-element list. */
std::vector<bool> flags_from_json(const Json& j, const std::string& where) {
    std::vector<bool> out;
    const Json flat = j.is_array() ? j : Json::array({j});
    for (std::size_t k = 0; k < flat.size(); ++k) {
        const Json& e = flat[k];
        if (e.is_boolean()) {
            out.push_back(e.get<bool>());
            continue;
        }
        const std::string text = decimal_text(e, at(where, k));
        if (text != "0" && text != "1")
            throw InputError(at(where, k) + ": a flag must be true, false, 0 or 1, got " + text);
        out.push_back(text == "1");
    }
    return out;
}
 
/** A list of numeric vectors, e.g. a MATLAB cell of row vectors or a matrix read by rows. */
template <class T>
std::vector<std::vector<T> > vectors_from_json(const Json& j, const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected a list of vectors");
    std::vector<std::vector<T> > out;
    for (std::size_t k = 0; k < j.size(); ++k)
        out.push_back(vector_from_json<T>(j[k], at(where, k)));
    return out;
}
 
/** A list of matrices. */
template <class T>
std::vector<Matrix<T> > matrices_from_json(const Json& j, const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected a list of matrices");
    std::vector<Matrix<T> > out;
    for (std::size_t k = 0; k < j.size(); ++k)
        out.push_back(matrix_from_json<T>(j[k], at(where, k)));
    return out;
}
 
/** A list of lists of matrices (a MATLAB cell of cells, e.g. R{v}{i}). */
template <class T>
std::vector<std::vector<Matrix<T> > > matrix_lists_from_json(const Json& j,
                                                             const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected a list of lists of matrices");
    std::vector<std::vector<Matrix<T> > > out;
    for (std::size_t k = 0; k < j.size(); ++k)
        out.push_back(matrices_from_json<T>(j[k], at(where, k)));
    return out;
}
 
template <class T>
Json encode_matrix_lists(const std::vector<std::vector<Matrix<T> > >& v) {
    Json a = Json::array();
    for (const auto& x : v) a.push_back(encode_matrices(x));
    return a;
}
 
/** One MAP written [D0, D1]. */
template <class T>
mam::Map<T> map_from_json(const Json& j, const std::string& where) {
    if (!j.is_array() || j.size() != 2)
        throw InputError(where + ": a MAP is the two-element list [D0, D1]");
    mam::Map<T> m;
    m.D0 = matrix_from_json<T>(j[0], where + ".D0");
    m.D1 = matrix_from_json<T>(j[1], where + ".D1");
    if (m.D0.rows() != m.D0.cols() || m.D1.rows() != m.D0.rows() || m.D1.cols() != m.D0.cols())
        throw InputError(where + ": D0 and D1 must be square matrices of the same order");
    return m;
}
 
/**
 * True when j is a LIST of MAPs rather than one MAP: its first element is a
 * pair whose first matrix is written as an array of rows, i.e. j[0][0][0] is
 * itself an array. A single MAP has a number there.
 */
bool is_map_list(const Json& j) {
    return j.is_array() && !j.empty() && j[0].is_array() && !j[0].empty() &&
           j[0][0].is_array() && !j[0][0].empty() && j[0][0][0].is_array();
}
 
template <class T>
std::vector<mam::Map<T> > maps_from_json(const Json& j, const std::string& where) {
    if (!is_map_list(j))
        throw InputError(where + ": expected a list of MAPs [[D0, D1], ...], each matrix "
                                 "written as an array of rows");
    std::vector<mam::Map<T> > out;
    for (std::size_t k = 0; k < j.size(); ++k) out.push_back(map_from_json<T>(j[k], at(where, k)));
    return out;
}
 
template <class T>
Json encode_map(const mam::Map<T>& m) {
    Json a = Json::array();
    a.push_back(encode_matrix(m.D0));
    a.push_back(encode_matrix(m.D1));
    return a;
}
 
template <class T>
Json encode_maps(const std::vector<mam::Map<T> >& v) {
    Json a = Json::array();
    for (const auto& m : v) a.push_back(encode_map(m));
    return a;
}
 
std::string where_of(const char* fn, const char* key) { return std::string(fn) + ": " + key; }
 
// ---------------------------------------------------------------------------
// polling
// ---------------------------------------------------------------------------
 
/**
 * The moments the polling formulas consume, reduced exactly as MATLAB's
 * private/polling_moments.m does: each entry of arvMAPs, svcMAPs and
 * switchMAPs is either a MAP [D0, D1] or a numeric moment vector, namely
 * lambda_i, [E[S_i] E[S_i^2]] and [E[R_i] Var[R_i]] respectively.
 */
// per_queue_sw = false lets switchMAPs carry a different number of entries
// (one per table position for polling_qsys_table); r and delta2 then hold the
// moments of those entries in order.
template <class T>
polling::PollingMoments<T> polling_moments_from_json(const Json& arv, const Json& svc, const Json& sw,
                                                     const char* fn, bool per_queue_sw) {
    for (const Json* j : {&arv, &svc, &sw})
        if (!j->is_array()) throw InputError(std::string(fn) + ": arvMAPs, svcMAPs and switchMAPs are lists");
    const std::size_t n = arv.size();
    if (svc.size() != n || (per_queue_sw && sw.size() != n))
        throw InputError(std::string(fn) + ": arvMAPs, svcMAPs and switchMAPs must have one entry per queue");
    // An entry is a MAP when it is a pair whose elements are arrays; a scalar
    // or a flat list of numbers is a moment vector.
    const auto is_map = [](const Json& e) {
        return e.is_array() && e.size() == 2 && e[0].is_array() && e[1].is_array();
    };
    polling::PollingMoments<T> m;
    for (std::size_t i = 0; i < n; ++i) {
        const std::string wa = at(std::string(fn) + ": arvMAPs", i);
        const std::string ws = at(std::string(fn) + ": svcMAPs", i);
        if (is_map(arv[i])) {
            m.lambda.push_back(mam::map_lambda(map_from_json<T>(arv[i], wa)));
        } else {
            const std::vector<T> v = vector_from_json<T>(arv[i], wa);
            if (v.empty()) throw InputError(wa + ": empty moment vector");
            m.lambda.push_back(v[0]);
        }
        if (is_map(svc[i])) {
            const mam::Map<T> s = map_from_json<T>(svc[i], ws);
            m.b.push_back(mam::map_mean(s));
            m.b2.push_back(mam::map_moment(s, 2));
        } else {
            const std::vector<T> v = vector_from_json<T>(svc[i], ws);
            if (v.size() < 2) throw InputError(ws + ": expected [E[S] E[S^2]]");
            m.b.push_back(v[0]);
            m.b2.push_back(v[1]);
        }
    }
    for (std::size_t i = 0; i < sw.size(); ++i) {
        const std::string wr = at(std::string(fn) + ": switchMAPs", i);
        if (is_map(sw[i])) {
            const mam::Map<T> r = map_from_json<T>(sw[i], wr);
            m.r.push_back(mam::map_mean(r));
            m.delta2.push_back(mam::map_var(r));
        } else {
            const std::vector<T> v = vector_from_json<T>(sw[i], wr);
            if (v.size() < 2) throw InputError(wr + ": expected [E[R] Var[R]]");
            m.r.push_back(v[0]);
            m.delta2.push_back(v[1]);
        }
    }
    return m;
}
 
template <class T>
polling::PollingMoments<T> polling_moments_from_args(Args& a, const char* fn) {
    const Json& arv = a.get("arvMAPs");
    const Json& svc = a.get("svcMAPs");
    const Json& sw = a.get("switchMAPs");
    a.done();
    return polling_moments_from_json<T>(arv, svc, sw, fn, true);
}
 
#define LINE_POLLING_OP(Name, fn)                                                    \
    struct Name {                                                                    \
        template <class T>                                                           \
        static Json run(Args& a) {                                                   \
            const polling::PollingMoments<T> m = polling_moments_from_args<T>(a, #fn); \
            Json r;                                                                  \
            r["W"] = encode_vector(polling::fn(m));                                  \
            return r;                                                                \
        }                                                                            \
    };
LINE_POLLING_OP(OpPolling1Limited, polling_qsys_1limited)
LINE_POLLING_OP(OpPollingDecrementing, polling_qsys_decrementing)
LINE_POLLING_OP(OpPollingExhaustive, polling_qsys_exhaustive)
LINE_POLLING_OP(OpPollingGated, polling_qsys_gated)
#undef LINE_POLLING_OP
 
/**
 * polling_qsys_table(arvMAPs, svcMAPs, switchMAPs, table, discipline):
 * switchMAPs holds one entry per table position, table is 1-based and
 * discipline is "exhaustive" or "gated" (or the PollingType id of
 * EXHAUSTIVE or GATED). Returns W, EX (M x n) and EX2 (one n x n per position).
 */
struct OpPollingTable {
    template <class T>
    static Json run(Args& a) {
        const char* fn = "polling_qsys_table";
        const Json& arv = a.get("arvMAPs");
        const Json& svc = a.get("svcMAPs");
        const Json& sw = a.get("switchMAPs");
        const std::vector<int> tab = a.ints("table");
        const Json& disc = a.get("discipline");
        a.done();
        bool gated = false;
        if (disc.is_string()) {
            std::string d = disc.get<std::string>();
            for (auto& ch : d) ch = static_cast<char>(std::tolower(static_cast<unsigned char>(ch)));
            if (d == "exhaustive") gated = false;
            else if (d == "gated") gated = true;
            else throw InputError(ErrorKind{"Argument"}, std::string(fn) + ": The discipline must be 'exhaustive' or 'gated'.");
        } else {
            const std::vector<int> id = int_vector_from_json(disc, std::string(fn) + ": discipline");
            if (id.size() == 1 && id[0] == static_cast<int>(lang::PollingType::EXHAUSTIVE)) gated = false;
            else if (id.size() == 1 && id[0] == static_cast<int>(lang::PollingType::GATED)) gated = true;
            else throw UnsupportedError(ErrorKind{"SchedStrategy_POLLING"},
                                        std::string(fn) + ": polling_qsys_table is exact for exhaustive or gated service only.");
        }
        const polling::PollingMoments<T> m = polling_moments_from_json<T>(arv, svc, sw, fn, false);
        std::vector<std::size_t> table;
        for (int t : tab) {
            if (t < 1) throw InputError(ErrorKind{"Argument"}, std::string(fn) + ": The polling table must list every queue 1.." +
                                        std::to_string(m.lambda.size()) + " at least once.");
            table.push_back(static_cast<std::size_t>(t));
        }
        const polling::PollingTableResult<T> res =
            polling::polling_qsys_table_full(m.lambda, m.b, m.b2, table, m.r, m.delta2, gated);
        Json r;
        r["W"] = encode_vector(res.W);
        Matrix<T> EX(res.EX.size(), m.lambda.size());
        for (std::size_t p = 0; p < res.EX.size(); ++p)
            for (std::size_t i = 0; i < m.lambda.size(); ++i) EX(p, i) = res.EX[p][i];
        r["EX"] = encode_matrix(EX);
        Json ex2 = Json::array();
        for (const auto& S : res.EX2) ex2.push_back(encode_matrix(S));
        r["EX2"] = ex2;
        return r;
    }
};
 
/** polling_qsys_klimited(arvMAPs, svcMAPs, switchMAPs, K): K scalar or per queue. */
struct OpPollingKLimited {
    template <class T>
    static Json run(Args& a) {
        const char* fn = "polling_qsys_klimited";
        const Json& arv = a.get("arvMAPs");
        const Json& svc = a.get("svcMAPs");
        const Json& sw = a.get("switchMAPs");
        const std::vector<int> k = a.ints("K");
        a.done();
        const polling::PollingMoments<T> m = polling_moments_from_json<T>(arv, svc, sw, fn, true);
        Json r;
        r["W"] = encode_vector(polling::polling_qsys_klimited(m, std::vector<long>(k.begin(), k.end())));
        return r;
    }
};
 
/**
 * polling_pcl(arvMAPs, svcMAPs, switchMAPs, disc): disc is a string of
 * 'E','G','L','D' (one per queue, or one letter for all), or a PollingType id
 * (KLIMITED and GLIMITED stand for K=1). Returns w, rhs and pcl.
 */
struct OpPollingPcl {
    template <class T>
    static Json run(Args& a) {
        const char* fn = "polling_pcl";
        const Json& arv = a.get("arvMAPs");
        const Json& svc = a.get("svcMAPs");
        const Json& sw = a.get("switchMAPs");
        const Json& disc = a.get("disc");
        a.done();
        const polling::PollingMoments<T> m = polling_moments_from_json<T>(arv, svc, sw, fn, true);
        const std::size_t n = m.lambda.size();
        std::string code;
        if (disc.is_string()) {
            code = disc.get<std::string>();
            if (code.size() == 1) code = std::string(n, code[0]);
        } else {
            for (int id : int_vector_from_json(disc, std::string(fn) + ": disc")) {
                switch (static_cast<lang::PollingType>(id)) {
                    case lang::PollingType::EXHAUSTIVE: code += 'E'; break;
                    case lang::PollingType::GATED: code += 'G'; break;
                    case lang::PollingType::KLIMITED:
                    case lang::PollingType::GLIMITED: code += 'L'; break;
                    case lang::PollingType::DECREMENTING: code += 'D'; break;
                    default:
                        throw UnsupportedError(ErrorKind{"SchedStrategy_POLLING"},
                                               std::string(fn) + ": unsupported discipline.");
                }
            }
            if (code.size() == 1) code = std::string(n, code[0]);
        }
        const polling::PollingPcl<T> p = polling::polling_pcl(m, code);
        Json r;
        r["w"] = encode_vector(p.w);
        r["rhs"] = encode_scalar(p.rhs);
        r["pcl"] = encode_scalar(p.pcl);
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// fes
// ---------------------------------------------------------------------------
 
struct OpLjdLinearize {
    template <class T>
    static Json run(Args& a) {
        const std::vector<int> nvec = a.ints("nvec");
        const std::vector<int> cutoffs = a.ints("cutoffs");
        a.done();
        Json r;
        r["idx"] = encode_count(fes::ljd_linearize(nvec, cutoffs));  // 1-based, as MATLAB's
        return r;
    }
};
 
/** The (L, mi, isDelay, cutoffs) quadruple shared by fes_compute_throughputs and _metrics. */
template <class T>
struct FesLoad {
    Matrix<T> L;
    std::vector<int> mi;
    std::vector<bool> isDelay;
    std::vector<int> cutoffs;
};
 
template <class T>
FesLoad<T> read_fes_load(Args& a, const char* fn) {
    FesLoad<T> f;
    f.L = a.matrix<T>("L");
    f.mi = a.ints("mi");
    f.isDelay = flags_from_json(a.get("isDelay"), where_of(fn, "isDelay"));
    f.cutoffs = a.ints("cutoffs");
    a.unsupported("options", "its only field, verbose, prints progress and changes no result; "
                             "the port does not print.");
    a.done();
    return f;
}
 
struct OpFesComputeThroughputs {
    template <class T>
    static Json run(Args& a) {
        const FesLoad<T> f = read_fes_load<T>(a, "fes_compute_throughputs");
        Json r;
        r["scalingTable"] =
            encode_vectors(fes::fes_compute_throughputs(f.L, f.mi, f.isDelay, f.cutoffs));
        return r;
    }
};
 
struct OpFesComputeMetrics {
    template <class T>
    static Json run(Args& a) {
        const FesLoad<T> f = read_fes_load<T>(a, "fes_compute_metrics");
        const fes::FesConditionalMetrics<T> m =
            fes::fes_compute_metrics(f.L, f.mi, f.isDelay, f.cutoffs);
        Json r;
        r["Qtable"] = encode_matrices(m.QN);
        r["Utable"] = encode_matrices(m.UN);
        return r;
    }
};
 
/**
 * fes_build_isolated: MATLAB takes (sn, subsetIndices, stochCompS) and plucks
 * the subset's service rates out of sn; the port lifts that extraction into
 * its signature and takes the rates directly (fes_build_isolated.h). The
 * entry reads them under the port's name, `rates` (M_sub x K, the rows of
 * sn.rates at the subset), and refuses sn and subsetIndices. mi and isDelay,
 * which MATLAB returns only as a pass-through of sn, are not outputs here.
 */
struct OpFesBuildIsolated {
    template <class T>
    static Json run(Args& a) {
        a.unsupported("sn", "the port takes the subset's service rates as `rates` instead of "
                            "a NetworkStruct.");
        a.unsupported("subsetIndices", "the port takes the subset's rows of sn.rates as "
                                       "`rates`, so the subset is already selected.");
        const Matrix<T> rates = a.matrix<T>("rates");
        const Matrix<T> S = a.matrix<T>("stochCompS");
        a.done();
        const fes::FesIsolated<T> iso = fes::fes_build_isolated(rates, S);
        Json r;
        r["L"] = encode_matrix(iso.L);
        r["visits"] = encode_matrix(iso.visits);
        return r;
    }
};
 
/**
 * A flow-equivalent-server descriptor as MATLAB's fes_map_levels reads it:
 * a single MAP replicated (and, given mi, scaled by min(k, mi)) over n
 * levels, or a list of per-level MAPs returned as given.
 */
template <class T>
std::vector<mam::Map<T> > fes_levels_from_json(const Json& j, std::size_t n, const double* mi,
                                               const std::string& where) {
    if (is_map_list(j)) return fes::fes_map_levels(maps_from_json<T>(j, where), n);
    const mam::Map<T> m = map_from_json<T>(j, where);
    return mi ? fes::fes_map_levels(m, n, *mi) : fes::fes_map_levels(m, n);
}
 
struct OpFesMapLevels {
    template <class T>
    static Json run(Args& a) {
        const Json& F = a.get("FES");
        const std::size_t n = a.count("n");
        const bool hasMi = a.has("mi");
        const double mi = hasMi ? a.number<double>("mi") : 0.0;
        a.done();
        Json r;
        r["FESlev"] = encode_maps(
            fes_levels_from_json<T>(F, n, hasMi ? &mi : nullptr, "fes_map_levels: FES"));
        return r;
    }
};
 
struct OpFesMapInterdeparture {
    template <class T>
    static Json run(Args& a) {
        const Json& M = a.get("MAPs");
        const Json& F = a.get("FES");
        const std::size_t n = a.count("n");
        // mi = [mi_station mi_fes] scales a single MAP over the levels exactly
        // as MATLAB's fes_map_levels call inside fes_map_interdeparture does.
        const bool hasMi = a.has("mi");
        std::vector<double> mi;
        if (hasMi) {
            mi = a.vector<double>("mi");
            if (mi.size() != 2)
                throw InputError("fes_map_interdeparture: mi is [mi_station mi_fes]");
        }
        a.done();
        const std::vector<mam::Map<T> > maps = fes_levels_from_json<T>(
            M, n, hasMi ? &mi[0] : nullptr, "fes_map_interdeparture: MAPs");
        const std::vector<mam::Map<T> > fes = fes_levels_from_json<T>(
            F, n, hasMi ? &mi[1] : nullptr, "fes_map_interdeparture: FES");
        const mam::Map<T> out = fes::fes_map_interdeparture(maps, fes, n);
        Json r;
        r["T0"] = encode_matrix(out.D0);
        r["T1"] = encode_matrix(out.D1);
        return r;
    }
};
 
struct OpFesMapInterp {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> x = a.vector<T>("x");
        const std::vector<T> y = a.vector<T>("y");
        const std::vector<T> xq = a.vector<T>("xq");
        a.done();
        Json r;
        r["yq"] = encode_vector(fes::fes_map_interp(x, y, xq));
        return r;
    }
};
 
/**
 * fes_map_moments(T0, T1, options). The port's options are POSITIONAL
 * (method, step_safety, tol, iter_max), so an absent field is honoured by
 * calling with the leading ones only; a field given after an absent one would
 * need the absent one's default restated, and is refused instead.
 */
struct OpFesMapMoments {
    template <class T>
    static Json run(Args& a) {
        mam::Map<T> m;
        m.D0 = a.matrix<T>("T0");
        m.D1 = a.matrix<T>("T1");
        const Json* oj = a.opt("options");
        a.done();
        const Json empty = Json::object();
        Args o(oj ? *oj : empty, "fes_map_moments options");
        const char* order[] = {"method", "step_safety", "tol", "iter_max"};
        std::size_t given = 0;
        for (std::size_t k = 0; k < 4; ++k)
            if (o.has(order[k])) given = k + 1;
        for (std::size_t k = 0; k < given; ++k)
            if (!o.has(order[k]))
                throw InputError(std::string("fes_map_moments: options.") + order[given - 1] +
                                 " is given without options." + order[k] +
                                 "; the port takes these positionally (method, step_safety, "
                                 "tol, iter_max), so give every field before the last one used");
        const std::string method = given > 0 ? o.text("method", "") : std::string();
        const double step = given > 1 ? o.number<double>("step_safety") : 0.0;
        const double tol = given > 2 ? o.number<double>("tol") : 0.0;
        const std::size_t iters = given > 3 ? o.count("iter_max") : 0;
        o.done();
        fes::FesMapMoments<T> r0;
        switch (given) {
            case 0: r0 = fes::fes_map_moments(m); break;
            case 1: r0 = fes::fes_map_moments(m, method); break;
            case 2: r0 = fes::fes_map_moments(m, method, step); break;
            case 3: r0 = fes::fes_map_moments(m, method, step, tol); break;
            default: r0 = fes::fes_map_moments(m, method, step, tol, iters); break;
        }
        Json r;
        r["e1"] = encode_scalar(r0.e1);
        r["e2"] = encode_scalar(r0.e2);
        r["e3"] = encode_scalar(r0.e3);
        r["e11"] = encode_scalar(r0.e11);
        r["idc"] = encode_scalar(r0.idc);
        return r;
    }
};
 
struct OpFesMapSolve {
    template <class T>
    static Json run(Args& a) {
        const Json& F = a.get("FES");
        const mam::Map<T> think = map_from_json<T>(a.get("thinkMAP"), "fes_map_solve: thinkMAP");
        const std::size_t N = a.count("N");
        a.done();
        const std::vector<mam::Map<T> > fes = fes_levels_from_json<T>(F, N, nullptr, "fes_map_solve: FES");
        const fes::FesMapSolveResult<T> s = fes::fes_map_solve(fes, think, N);
        Json r;
        r["XN"] = encode_scalar(s.X);
        r["RN"] = encode_scalar(s.R);
        r["QN"] = encode_scalar(s.Q);
        r["pk"] = encode_vector(s.pk);
        return r;
    }
};
 
struct OpFesMapAggregate {
    template <class T>
    static Json run(Args& a) {
        const std::vector<mam::Map<T> > maps = maps_from_json<T>(a.get("maps"), "fes_map_aggregate: maps");
        const std::vector<double> servers = a.vector<double>("servers");
        const std::size_t n = a.count("n");
        const Json* oj = a.opt("options");
        a.done();
        const Json empty = Json::object();
        Args o(oj ? *oj : empty, "fes_map_aggregate options");
        o.unsupported("verbose", "it prints progress and changes no result; the port does not print.");
        const bool hasGrid = o.has("grid");
        const std::vector<std::size_t> grid =
            hasGrid ? sizes_from_json(o.get("grid"), "fes_map_aggregate: options.grid")
                    : std::vector<std::size_t>();
        o.done();
        const fes::FesMapAggregateResult<T> s =
            hasGrid ? fes::fes_map_aggregate(maps, servers, n, grid)
                    : fes::fes_map_aggregate(maps, servers, n);
        Json info;
        info["throughput"] = encode_vector(s.throughput);
        info["moments"] = encode_vectors(s.moments);
        info["status"] = encode_ints(s.status);
        info["grid"] = encode_sizes(s.grid);
        Json r;
        r["FES"] = encode_maps(s.fes);
        r["info"] = info;
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// cache
//
// gamma is the (n x h) access-factor matrix and m the (h) list capacities. A
// MATLAB request-rate array lambda(u, n, h+1) is read by the ports at its first
// page only, so `lambda` is that (u x n) page, except where a port takes the
// per-list rates (cache_gamma, cache_gamma_lp, cache_ttl_lrua), in which case it
// is the list of the u per-user (n x h+1) matrices and R the per-user list of
// per-item (h+1 x h+1) access graphs.
// ---------------------------------------------------------------------------
 
template <class T>
Json cache_miss_results(const T& M, const std::vector<T>& MU, const std::vector<T>& MI,
                        const std::vector<T>& pi0) {
    Json r;
    r["M"] = encode_scalar(M);
    r["MU"] = encode_vector(MU);
    r["MI"] = encode_vector(MI);
    r["pi0"] = encode_vector(pi0);
    return r;
}
 
/** The optional size-constrained pair (sigma, k): both or neither. */
struct SizeCaps {
    bool given = false;
    std::vector<int> sigma, k;
};
 
SizeCaps read_size_caps(Args& a, const char* fn, const char* capKey) {
    SizeCaps s;
    const bool hs = a.has("sigma"), hk = a.has(capKey);
    if (hs != hk)
        throw InputError(std::string(fn) + ": sigma and " + capKey +
                         " constrain the cache together; give both or neither");
    s.given = hs;
    s.sigma = a.ints_or_empty("sigma");
    s.k = a.ints_or_empty(capKey);
    return s;
}
 
struct OpCacheErec {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const SizeCaps s = read_size_caps(a, "cache_erec", "k");
        a.done();
        Json r;
        r["E"] = encode_scalar(s.given ? cache::cache_erec(gamma, m, s.sigma, s.k)
                                       : cache::cache_erec(gamma, m));
        return r;
    }
};
 
struct OpCacheProbErec {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const SizeCaps s = read_size_caps(a, "cache_prob_erec", "k");
        a.done();
        Json r;
        r["prob"] = encode_matrix(s.given ? cache::cache_prob_erec(gamma, m, s.sigma, s.k)
                                          : cache::cache_prob_erec(gamma, m));
        return r;
    }
};
 
/**
 * cache_cost(gamma, m, sigma[, k[, pij]]) -> K, pij. As in the reference, an
 * absent pij is the cache_prob_erec occupancy, and it is returned either way.
 */
struct OpCacheCost {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const std::vector<int> sigma = a.ints("sigma");
        const std::vector<int> k = a.ints_or_empty("k");
        const bool hasPij = a.has("pij");
        Matrix<T> pij = a.matrix_or_empty<T>("pij");
        a.done();
        if (!hasPij) pij = cache::cache_prob_erec(gamma, m, sigma, k);
        Json r;
        r["K"] = encode_vector(cache::cache_cost(gamma, m, sigma, k, pij));
        r["pij"] = encode_matrix(pij);
        return r;
    }
};
 
/** viol: one row [item, list, blockinglist] per blocked pair, 0-based like parent. */
struct OpCacheCostPathcheck {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> sigma = a.ints("sigma");
        const std::vector<int> k = a.ints("k");
        const std::vector<int> parent = a.ints("parent");  // 0-based, -1 = rooted in the miss list
        a.done();
        Json viol = Json::array();
        for (const cache::CacheBlockedPair& p : cache::cache_cost_pathcheck(gamma, sigma, k, parent))
            viol.push_back(encode_sizes({p.item, p.list, p.blocking_list}));
        Json r;
        r["viol"] = viol;
        return r;
    }
};
 
struct OpCacheSpmSize {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const std::vector<int> sigma = a.ints("sigma");
        const std::vector<int> k = a.ints("k");
        const bool hasMode = a.has("costmode");
        const std::string mode = a.text("costmode", "");
        a.done();
        cache::CacheSpmSizeResult<T> s;
        if (hasMode) {
            cache::CacheCostMode cm;
            if (mode == "atmost" || mode == "ATMOST" || mode == "AtMost")
                cm = cache::CacheCostMode::AtMost;
            else if (mode == "exact" || mode == "EXACT" || mode == "Exact")
                cm = cache::CacheCostMode::Exact;
            else
                throw InputError("cache_spm_size: costmode must be 'atmost' or 'exact', got '" +
                                 mode + "'");
            s = cache::cache_spm_size(gamma, m, sigma, k, cm);
        } else {
            s = cache::cache_spm_size(gamma, m, sigma, k);
        }
        Json out;
        out["xi"] = encode_vector(s.xi);
        out["zeta"] = encode_vector(s.zeta);
        Json b = Json::array();
        for (bool x : s.binding) b.push_back(x);
        out["binding"] = b;
        out["pij"] = encode_matrix(s.pij);
        out["K"] = encode_vector(s.k_mean);
        out["phi"] = encode_scalar(s.phi);
        out["logdetSigma"] = encode_scalar(s.logdet_sigma);
        out["iter"] = encode_count(s.iterations);
        out["relerrEst"] = s.relerr_est;
        out["method"] = s.method;
        out["span"] = s.span;
        Json r;
        r["E"] = encode_scalar(s.e);
        r["logE"] = encode_scalar(s.log_e);
        r["out"] = out;
        return r;
    }
};
 
struct OpCacheGamma {
    template <class T>
    static Json run(Args& a) {
        const std::vector<Matrix<T> > lambda = matrices_from_json<T>(a.get("lambda"), "cache_gamma: lambda");
        const auto R = matrix_lists_from_json<T>(a.get("R"), "cache_gamma: R");
        a.done();
        const cache::CacheGammaGraphResult<T> g = cache::cache_gamma(lambda, R);
        Json r;
        r["gamma"] = encode_matrix(g.gamma);
        r["u"] = encode_count(g.u);
        r["n"] = encode_count(g.n);
        r["h"] = encode_count(g.h);
        return r;
    }
};
 
struct OpCacheGammaLp {
    template <class T>
    static Json run(Args& a) {
        const std::vector<Matrix<T> > lambda = matrices_from_json<T>(a.get("lambda"), "cache_gamma_lp: lambda");
        const auto R = matrix_lists_from_json<T>(a.get("R"), "cache_gamma_lp: R");
        a.done();
        const cache::CacheGammaResult<T> g = cache::cache_gamma_lp(lambda, R);
        Json r;
        r["gamma"] = encode_matrix(g.gamma);
        r["u"] = encode_count(g.u);
        r["n"] = encode_count(g.n);
        r["h"] = encode_count(g.h);
        r["parent"] = encode_ints(g.parent);  // 0-based, -1 = rooted in the miss list
        return r;
    }
};
 
/** cache_miss_asy(gamma, m[, maxiter, tol]): the port takes the two together. */
struct OpCacheMissAsy {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const bool hi = a.has("maxiter"), ht = a.has("tol");
        if (hi != ht)
            throw InputError("cache_miss_asy: the port takes maxiter and tol together; give "
                             "both or neither");
        const int maxiter = a.integer("maxiter", 0);
        const T tol = a.scalar<T>("tol", T(0));
        a.done();
        Json r;
        r["missratio"] = encode_scalar(hi ? cache::cache_miss_asy(gamma, m, maxiter, tol)
                                          : cache::cache_miss_asy(gamma, m));
        return r;
    }
};
 
struct OpCacheMiss {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const bool hasLambda = a.has("lambda");
        const Matrix<T> lambda = a.matrix_or_empty<T>("lambda");
        const SizeCaps s = read_size_caps(a, "cache_miss", "cap");
        a.done();
        const cache::CacheMissResult<T> x =
            hasLambda ? (s.given ? cache::cache_miss(gamma, m, lambda, s.sigma, s.k)
                                 : cache::cache_miss(gamma, m, lambda))
                      : (s.given ? cache::cache_miss(gamma, m, s.sigma, s.k)
                                 : cache::cache_miss(gamma, m));
        return cache_miss_results(x.M, x.MU, x.MI, x.pi0);
    }
};
 
struct OpCacheMissFpi {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const Matrix<T> lambda = a.matrix<T>("lambda");
        a.done();
        const cache::CacheMissResult<T> x = cache::cache_miss_fpi(gamma, m, lambda);
        return cache_miss_results(x.M, x.MU, x.MI, x.pi0);
    }
};
 
struct OpCacheMissSpm {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const Matrix<T> lambda = a.matrix<T>("lambda");
        a.done();
        const cache::CacheMissSpmResult<T> x = cache::cache_miss_spm(gamma, m, lambda);
        Json r = cache_miss_results(x.M, x.MU, x.MI, x.pi0);
        r["lE"] = encode_scalar(x.lE);
        return r;
    }
};
 
/**
 * cache_miss_rmf(gamma, m, lambda[, tspan, x0init, accost, wgen]).
 *
 * The steady state is the port's cache_miss_rmf; a tspan adds the transient
 * of cache_miss_rmf_transient over [tspan(1), tspan(end)], from x0init, with
 * the sojourn clock wgen = {A, c, phi0} when given. Two port gaps surface as
 * refusals: the transient port has no default initial occupancy, so tspan
 * needs x0init; and the accost overload takes the fixed-point horizon tmax
 * explicitly (the reference fixes it), so accost needs tmax. tmax is a
 * port-only key, optional otherwise.
 */
struct OpCacheMissRmf {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gm = a.matrix<T>("gamma");
        const std::vector<T> gamma(gm.data(), gm.data() + gm.size());
        const std::vector<int> m = a.ints("m");
        const Matrix<T> lambda = a.matrix<T>("lambda");
        const bool hasTmax = a.has("tmax");
        const T tmax = a.scalar<T>("tmax", T(0));
        const bool hasAccost = a.has("accost");
        const auto accost = hasAccost
                                ? matrix_lists_from_json<T>(a.get("accost"), "cache_miss_rmf: accost")
                                : std::vector<std::vector<Matrix<T> > >();
        const bool hasTspan = a.has("tspan");
        const std::vector<T> tspan = a.vector_or_empty<T>("tspan");
        const bool hasX0 = a.has("x0init");
        const std::vector<T> x0init = a.vector_or_empty<T>("x0init");
        const Json* wj = a.opt("wgen");
        a.done();
        if (hasAccost && !hasTmax)
            throw InputError("cache_miss_rmf: the port's access-graph overload takes the "
                             "fixed-point horizon explicitly; give tmax with accost");
        if (hasX0 && !hasTspan)
            throw InputError("cache_miss_rmf: x0init is the initial occupancy of the transient "
                             "and needs tspan");
        if (wj && !hasTspan)
            throw InputError("cache_miss_rmf: wgen weights the transient and needs tspan");
        if (hasTspan && !hasX0)
            throw InputError("cache_miss_rmf: the port's transient (cache_miss_rmf_transient) "
                             "has no default initial occupancy; give x0init with tspan");
        if (hasTspan && tspan.size() != 2)
            throw InputError("cache_miss_rmf: tspan is [t0 t1]; the port integrates between "
                             "the two end points only");
        cache::CacheSojournClock<T> clk;
        if (wj) {
            Args w(*wj, "cache_miss_rmf wgen");
            clk.A = w.matrix<T>("A");
            clk.c = w.vector<T>("c");
            clk.phi0 = w.vector<T>("phi0");
            w.done();
        }
        const cache::CacheMissRmfResult<T> s =
            hasAccost ? cache::cache_miss_rmf(gamma, m, lambda, tmax, accost)
                      : (hasTmax ? cache::cache_miss_rmf(gamma, m, lambda, tmax)
                                 : cache::cache_miss_rmf(gamma, m, lambda));
        Json r = cache_miss_results(s.M, s.MU, s.MI, s.pi0);
        r["xss"] = encode_vector(s.xss);
        r["Wcov"] = encode_matrix(s.Wcov);
        if (hasTspan) {
            const cache::CacheMissRmfResult<T> tr =
                hasAccost ? cache::cache_miss_rmf_transient(m, lambda, tspan[0], tspan[1], x0init,
                                                            accost, clk)
                          : cache::cache_miss_rmf_transient(m, lambda, tspan[0], tspan[1], x0init,
                                                            clk);
            r["tout"] = encode_vector(tr.tout);
            r["pi0_t"] = encode_matrix(tr.pi0_t);
            r["MU_t"] = encode_matrix(tr.MU_t);
            r["xtraj"] = encode_matrix(tr.xtraj);
            if (wj) {
                Json sj;
                sj["xbar"] = encode_vector(tr.sojourn.xbar);
                sj["wtot"] = encode_scalar(tr.sojourn.wtot);
                sj["MU"] = encode_vector(tr.sojourn.MU);
                r["sj"] = sj;
            }
        }
        return r;
    }
};
 
/** (n*(h+1)) flat occupancy x[k + s*n] as the (n x h+1) matrix MATLAB reshapes it to. */
template <class T>
Matrix<T> occupancy_matrix(const std::vector<T>& x, std::size_t n) {
    const std::size_t c = n ? x.size() / n : 0;
    Matrix<T> X(n, c);
    for (std::size_t s = 0; s < c; ++s)
        for (std::size_t k = 0; k < n; ++k) X(k, s) = x[k + s * n];
    return X;
}
 
/**
 * cache_rrm_meanfield is a SCRIPT in MATLAB with its data fixed inline; the
 * port takes that data as (lambda, m[, tmax]), and the entry reads it under
 * those names. x is returned as the (n x h+1) matrix the script reshapes.
 */
struct OpCacheRrmMeanfield {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const std::vector<int> m = a.ints("m");
        const bool hasTmax = a.has("tmax");
        const T tmax = a.scalar<T>("tmax", T(0));
        a.done();
        const cache::CacheRrmMeanfieldResult<T> s =
            hasTmax ? cache::cache_rrm_meanfield(lambda, m, tmax) : cache::cache_rrm_meanfield(lambda, m);
        Json r;
        r["x"] = encode_matrix(occupancy_matrix(s.x, lambda.size()));
        r["missrate"] = encode_scalar(s.missrate);
        r["missratio"] = encode_scalar(s.missratio);
        return r;
    }
};
 
/**
 * cache_rrm_meanfield_ode(t, x, lambda, m, n, h). The drift is autonomous, so
 * the time argument (unnamed in MATLAB) is not taken. The port derives n and h
 * from lambda and m; when given they are checked against those sizes.
 */
struct OpCacheRrmMeanfieldOde {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> X = a.matrix<T>("x");
        const std::vector<T> lambda = a.vector<T>("lambda");
        const std::vector<int> m = a.ints("m");
        const bool hn = a.has("n"), hh = a.has("h");
        const std::size_t n = a.count("n", 0), h = a.count("h", 0);
        a.done();
        if (hn && n != lambda.size())
            throw InputError("cache_rrm_meanfield_ode: n = " + std::to_string(n) +
                             " but lambda has " + std::to_string(lambda.size()) + " items");
        if (hh && h != m.size())
            throw InputError("cache_rrm_meanfield_ode: h = " + std::to_string(h) + " but m has " +
                             std::to_string(m.size()) + " lists");
        // x is the flat column x(:) of the (n x h+1) occupancy, or that matrix itself.
        std::vector<T> x;
        if (X.rows() == lambda.size() && X.cols() == m.size() + 1) {
            for (std::size_t s = 0; s < X.cols(); ++s)
                for (std::size_t k = 0; k < X.rows(); ++k) x.push_back(X(k, s));
        } else {
            x.assign(X.data(), X.data() + X.size());
        }
        Json r;
        r["dxdt"] = encode_vector(cache::cache_rrm_meanfield_ode(x, lambda, m));
        return r;
    }
};
 
struct OpCacheMva {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        a.done();
        const cache::CacheMvaResult<T> s = cache::cache_mva(gamma, m);
        Json r;
        r["pi"] = encode_vector(s.pi);
        r["pi0"] = encode_vector(s.pi0);
        r["pij"] = encode_matrix(s.pij);
        r["x"] = encode_vector(s.x);
        r["u"] = encode_matrix(s.u);
        r["E"] = encode_scalar(s.E);
        return r;
    }
};
 
struct OpCacheMvaMiss {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> p = a.vector<T>("p");
        const std::vector<int> m = a.ints("m");
        const Matrix<T> R = a.matrix<T>("R");
        a.done();
        const cache::CacheMvaMissResult<T> s = cache::cache_mva_miss(p, m, R);
        Json r;
        r["M"] = encode_scalar(s.M);
        r["Mk"] = encode_vector(s.Mk);
        return r;
    }
};
 
/** The (gamma, m) -> prob family. */
#define LINE_CACHE_GM_MATRIX_OP(Name, fn, out, ...)              \
    struct Name {                                                \
        template <class T>                                       \
        static Json run(Args& a) {                               \
            const Matrix<T> gamma = a.matrix<T>("gamma");        \
            const std::vector<int> m = a.ints("m");              \
            __VA_ARGS__;                                         \
            a.done();                                            \
            Json r;                                              \
            r[out] = encode_matrix(cache::fn(gamma, m));         \
            return r;                                            \
        }                                                        \
    };
LINE_CACHE_GM_MATRIX_OP(OpCacheProbFpi, cache_prob_fpi, "prob", (void)0)
LINE_CACHE_GM_MATRIX_OP(OpCacheProbSpm, cache_prob_spm, "prob",
                        a.unsupported("lE", "the port computes the log normalizing constant "
                                            "itself and takes no precomputed value."))
#undef LINE_CACHE_GM_MATRIX_OP
 
struct OpCacheTtlHlru {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> lambda = a.matrix<T>("lambda");
        const std::vector<int> m = a.ints("m");
        a.done();
        Json r;
        r["pij"] = encode_matrix(cache::cache_ttl_hlru(lambda, m));
        return r;
    }
};
 
struct OpCacheTHlru {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        a.done();
        Json r;
        r["t"] = encode_vector(cache::cache_t_hlru(gamma, m));
        return r;
    }
};
 
struct OpCacheSpm {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        a.unsupported("xi0", "the port starts the saddle-point iteration from its own initial "
                             "point and takes none.");
        a.done();
        const cache::CacheSpmResult<T> s = cache::cache_spm(gamma, m);
        Json r;
        r["Z"] = encode_scalar(s.Z);
        r["lZ"] = encode_scalar(s.lZ);
        r["xi"] = encode_vector(s.xi);
        return r;
    }
};
 
struct OpCacheXiFp {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        a.unsupported("xi", "the port starts the fixed point from its own initial point and "
                            "takes none.");
        a.done();
        const cache::CacheXiFpResult<T> s = cache::cache_xi_fp(gamma, m);
        Json r;
        r["xi"] = encode_vector(s.xi);
        r["pi0"] = encode_vector(s.pi0);
        r["pij"] = encode_matrix(s.pij);
        r["it"] = s.it;
        return r;
    }
};
 
struct OpCacheXiIter {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        a.unsupported("tmax", "the port has no iteration-count argument.");
        a.done();
        Json r;
        r["z"] = encode_vector(cache::cache_xi_iter(gamma, m));
        return r;
    }
};
 
/**
 * Sampling controls of cache_is / cache_miss_is / cache_prob_is: samples and
 * the port-only seed come together (the port's overload takes both), and the
 * size constraint (sigma, k) needs them, since the port's constrained
 * overload is the one taking all four.
 */
struct IsControls {
    bool sampled = false;
    std::size_t samples = 0;
    std::uint64_t seed = 0;
    SizeCaps caps;
};
 
IsControls read_is_controls(Args& a, const char* fn) {
    IsControls c;
    const bool hs = a.has("samples"), hd = a.has("seed");
    if (hs != hd)
        throw InputError(std::string(fn) + ": the port takes samples and seed together; give "
                                           "both or neither");
    c.sampled = hs;
    c.samples = a.count("samples", 0);
    c.seed = static_cast<std::uint64_t>(a.count("seed", 0));
    c.caps = read_size_caps(a, fn, "k");
    if (c.caps.given && !c.sampled)
        throw InputError(std::string(fn) + ": the port's size-constrained overload also takes "
                                           "samples and seed; give them with sigma and k");
    return c;
}
 
struct OpCacheIs {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const IsControls c = read_is_controls(a, "cache_is");
        a.done();
        const cache::CacheIsResult<T> s =
            c.caps.given ? cache::cache_is(gamma, m, c.samples, c.seed, c.caps.sigma, c.caps.k)
                         : (c.sampled ? cache::cache_is(gamma, m, c.samples, c.seed)
                                      : cache::cache_is(gamma, m));
        Json r;
        r["E"] = encode_scalar(s.E);
        r["lE"] = encode_scalar(s.lE);
        return r;
    }
};
 
struct OpCacheProbIs {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const IsControls c = read_is_controls(a, "cache_prob_is");
        a.done();
        Json r;
        r["prob"] = encode_matrix(
            c.caps.given ? cache::cache_prob_is(gamma, m, c.samples, c.seed, c.caps.sigma, c.caps.k)
                         : (c.sampled ? cache::cache_prob_is(gamma, m, c.samples, c.seed)
                                      : cache::cache_prob_is(gamma, m)));
        return r;
    }
};
 
struct OpCacheMissIs {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        const Matrix<T> lambda = a.matrix<T>("lambda");
        const IsControls c = read_is_controls(a, "cache_miss_is");
        a.done();
        const cache::CacheMissIsResult<T> s =
            c.caps.given
                ? cache::cache_miss_is(gamma, m, lambda, c.samples, c.seed, c.caps.sigma, c.caps.k)
                : (c.sampled ? cache::cache_miss_is(gamma, m, lambda, c.samples, c.seed)
                             : cache::cache_miss_is(gamma, m, lambda));
        Json r = cache_miss_results(s.M, s.MU, s.MI, s.pi0);
        r["lE"] = encode_scalar(s.lE);
        return r;
    }
};
 
struct OpCacheLrumMapLevelstats {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> D0 = a.matrix<T>("D0");
        const Matrix<T> D1 = a.matrix<T>("D1");
        const std::vector<T> Tv = a.vector<T>("T");
        a.done();
        const cache::CacheLrumMapLevelStats<T> s = cache::cache_lrum_map_levelstats(D0, D1, Tv);
        Json r;
        r["prob"] = encode_vector(s.prob);
        r["occ"] = encode_vector(s.occ);
        r["hitfrac"] = encode_vector(s.hitfrac);
        return r;
    }
};
 
/** Per-item MAPs from MATLAB's two cells D0c{k}, D1c{k}. */
template <class T>
std::vector<mam::Map<T> > item_maps(Args& a, const char* fn) {
    const std::vector<Matrix<T> > D0 = matrices_from_json<T>(a.get("D0c"), where_of(fn, "D0c"));
    const std::vector<Matrix<T> > D1 = matrices_from_json<T>(a.get("D1c"), where_of(fn, "D1c"));
    if (D0.size() != D1.size())
        throw InputError(std::string(fn) + ": D0c and D1c must have one matrix per item");
    std::vector<mam::Map<T> > items(D0.size());
    for (std::size_t k = 0; k < D0.size(); ++k) {
        items[k].D0 = D0[k];
        items[k].D1 = D1[k];
    }
    return items;
}
 
/**
 * The LRU(m) MAP characteristic-time solvers. MATLAB stops on fsolve's
 * defaults; the port iterates to a relative tolerance it takes as a required
 * argument (internal callers pass GlobalConstants::FineTol), so `tol` is a
 * required port-only key here, and `maxswp` an optional one.
 */
struct OpCacheTLrumMap {
    template <class T>
    static Json run(Args& a) {
        const std::vector<mam::Map<T> > items = item_maps<T>(a, "cache_t_lrum_map");
        const std::vector<T> m = a.vector<T>("m");
        const T tol = a.number<T>("tol");
        const bool hasSw = a.has("maxswp");
        const unsigned maxswp = a.uinteger("maxswp", 0);
        a.done();
        Json r;
        r["t"] = encode_vector(hasSw ? cache::cache_t_lrum_map(items, m, tol, maxswp)
                                     : cache::cache_t_lrum_map(items, m, tol));
        return r;
    }
};
 
struct OpCacheTtlLrumMap {
    template <class T>
    static Json run(Args& a) {
        const std::vector<mam::Map<T> > items = item_maps<T>(a, "cache_ttl_lrum_map");
        const std::vector<T> m = a.vector<T>("m");
        const T tol = a.number<T>("tol");
        a.done();
        const cache::CacheTtlLrumMapResult<T> s = cache::cache_ttl_lrum_map(items, m, tol);
        Json r;
        r["pij"] = encode_matrix(s.pij);
        r["pijtime"] = encode_matrix(s.pijtime);
        return r;
    }
};
 
struct OpCacheTtlLrua {
    template <class T>
    static Json run(Args& a) {
        const std::vector<Matrix<T> > lambda = matrices_from_json<T>(a.get("lambda"), "cache_ttl_lrua: lambda");
        const auto R = matrix_lists_from_json<T>(a.get("R"), "cache_ttl_lrua: R");
        const std::vector<T> m = a.vector<T>("m");
        a.unsupported("seed", "the port solves for the characteristic times by deterministic "
                              "bisection and draws nothing.");
        const T tol = a.number<T>("tol");
        const bool hasSw = a.has("maxswp");
        const unsigned maxswp = a.uinteger("maxswp", 0);
        a.done();
        Json r;
        r["prob"] = encode_matrix(hasSw ? cache::cache_ttl_lrua(lambda, R, m, tol, maxswp)
                                        : cache::cache_ttl_lrua(lambda, R, m, tol));
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// da: decomposition-aggregation helpers
// ---------------------------------------------------------------------------
 
/**
 * da_cache_isolate(ch, lambda). ch carries MATLAB's sn.nodeparam fields:
 * itemcap (h), nitems, pread (one row per class; null or [] for MATLAB's NaN)
 * and accost (u x n of (h+1) x (h+1), or [] for the default linear routing).
 */
struct OpDaCacheIsolate {
    template <class T>
    static Json run(Args& a) {
        const Json& chj = a.get("ch");
        const std::vector<T> lambda = a.vector<T>("lambda");
        a.done();
        if (!chj.is_object()) throw InputError("da_cache_isolate: ch must be an object");
        Args ch(chj, "da_cache_isolate ch");
        da::CacheParam<T> p;
        p.itemcap = ch.ints("itemcap");
        p.nitems = ch.count("nitems");
        const Json& pr = ch.get("pread");
        if (!pr.is_array()) throw InputError("da_cache_isolate: ch.pread must be a list, one row per class");
        for (std::size_t v = 0; v < pr.size(); ++v) {
            if (pr[v].is_null() || (pr[v].is_array() && pr[v].empty())) {
                p.pread.push_back(std::vector<T>());
                continue;
            }
            p.pread.push_back(vector_from_json<T>(pr[v], at("da_cache_isolate: ch.pread", v)));
        }
        if (ch.has("accost")) {
            const Json& ac = ch.get("accost");
            if (!(ac.is_array() && ac.empty()))
                p.accost = matrix_lists_from_json<T>(ac, "da_cache_isolate: ch.accost");
        }
        ch.done();
        const da::CacheIsolateResult<T> s = da::da_cache_isolate(p, lambda);
        Json r;
        r["gamma"] = encode_matrix(s.gamma);
        r["lambda_cache"] = encode_matrices(s.lambda_cache);  // per class, n x (h+1)
        r["Rcost"] = encode_matrix_lists(s.Rcost);
        r["parent"] = encode_ints(s.parent);
        return r;
    }
};
 
struct OpDaTrafficSuperpos {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const std::vector<T> a2 = a.vector<T>("a2");
        a.done();
        Json r;
        r["d2"] = encode_scalar(da::da_traffic_superpos(lambda, a2));
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// lossn: loss networks
// ---------------------------------------------------------------------------
 
struct OpLossnErlangfp {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> nu = a.vector<T>("nu");
        const Matrix<T> A = a.matrix<T>("A");
        const std::vector<T> C = a.vector<T>("C");
        a.done();
        const lossn::ErlangFpResult<T> s = lossn::lossn_erlangfp(nu, A, C);
        Json r;
        r["QLen"] = encode_vector(s.QLen);
        r["Loss"] = encode_vector(s.Loss);
        r["E"] = encode_matrix(s.E);
        r["niter"] = encode_count(s.iterations);
        return r;
    }
};
 
/** lossn_manjunath(nu, A, C); options.max_live_states is the port's state cap. */
struct OpLossnManjunath {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> nu = a.vector<T>("nu");
        const Matrix<T> A = a.matrix<T>("A");
        const std::vector<T> C = a.vector<T>("C");
        lossn::LossnManjunathOptions o;
        if (a.has("options")) {
            Args op(*a.opt("options"), "lossn_manjunath options");
            if (op.has("max_live_states")) o.max_live_states = op.count("max_live_states");
            op.done();
        }
        a.done();
        const lossn::LossnManjunathResult<T> s = lossn::lossn_manjunath(nu, A, C, o);
        Json r;
        r["QLen"] = encode_vector(s.QLen);
        r["Loss"] = encode_vector(s.Loss);
        r["lG"] = s.lG;
        r["niter"] = encode_count(s.iterations);
        return r;
    }
};
 
/**
 * lossn_mci(nu, A, C, options): options {samples, alpha, gamma, seed}. The
 * port draws from a seeded generator only, so options.seed is required: MATLAB
 * falls back to the session RNG state, which has no JSON form.
 */
struct OpLossnMci {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> nu = a.vector<T>("nu");
        const Matrix<T> A = a.matrix<T>("A");
        const std::vector<T> C = a.vector<T>("C");
        const Json* oj = a.opt("options");
        a.done();
        if (!oj)
            throw InputError("lossn_mci: options.seed is required; the port draws from a seeded "
                             "generator and has no session RNG to fall back on");
        Args op(*oj, "lossn_mci options");
        lossn::LossnMciOptions<T> o;
        if (op.has("samples")) o.samples = op.count("samples");
        if (op.has("alpha")) o.alpha = op.scalar<double>("alpha", 0.0);
        if (op.has("gamma")) o.gamma = op.vector<T>("gamma");
        if (op.has("revenue")) o.revenue = op.vector<T>("revenue");
        if (!op.has("seed"))
            throw InputError("lossn_mci: options.seed is required; the port draws from a seeded "
                             "generator and has no session RNG to fall back on");
        const std::uint64_t seed = static_cast<std::uint64_t>(op.uinteger("seed"));
        op.done();
        const lossn::LossnMciResult<T> s = lossn::lossn_mci(nu, A, C, o, seed);
        Json ci;
        ci["accept"] = encode_matrix(s.acceptCI);
        ci["loss"] = encode_matrix(s.lossCI);
        ci["acceptPoint"] = encode_vector(s.acceptPoint);
        ci["lossPoint"] = encode_vector(s.lossPoint);
        ci["level"] = encode_scalar(s.level);
        ci["linkLossPoint"] = encode_matrix(s.linkLossPoint);
        if (s.hasRevenue) {
            ci["revenue"] = encode_scalar(s.revenue);
            ci["revsens"] = encode_vector(s.revsens);
            ci["revsensCI"] = encode_matrix(s.revsensCI);
            ci["impliedCost"] = encode_vector(s.impliedCost);
            ci["linkImpliedCost"] = encode_matrix(s.linkImpliedCost);
        }
        Json r;
        r["QLen"] = encode_vector(s.QLen);
        r["Loss"] = encode_vector(s.Loss);
        r["lG"] = s.lG;
        r["ci"] = ci;
        r["nsamples"] = encode_count(s.nsamples);
        return r;
    }
};
 
struct OpLossnRec {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> nu = a.vector<T>("nu");
        const Matrix<T> A = a.matrix<T>("A");
        const std::vector<T> C = a.vector<T>("C");
        a.done();
        const lossn::LossnRecResult<T> s = lossn::lossn_rec(nu, A, C);
        Json r;
        r["QLen"] = encode_vector(s.QLen);
        r["Loss"] = encode_vector(s.Loss);
        r["lG"] = s.lG;
        r["niter"] = s.iterations;
        return r;
    }
};
 
/** lossn_kr(rho, b, C): occupancy q, blocking B, carried load QLen. */
struct OpLossnKr {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> rho = a.vector<T>("rho");
        const std::vector<T> b = a.vector<T>("b");
        const T C = a.number<T>("C");
        a.done();
        const lossn::KrResult<T> s = lossn::lossn_kr(rho, b, C);
        Json r;
        r["q"] = encode_vector(s.q);
        r["B"] = encode_vector(s.B);
        r["QLen"] = encode_vector(s.QLen);
        return r;
    }
};
 
/** lossn_kr_dim(rho, b, Btarget); a scalar Btarget applies to every class. */
struct OpLossnKrDim {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> rho = a.vector<T>("rho");
        const std::vector<T> b = a.vector<T>("b");
        std::vector<T> Bt = a.vector<T>("Btarget");
        a.done();
        if (Bt.size() == 1) Bt.assign(rho.size(), Bt[0]);
        const lossn::KrDimResult<T> s = lossn::lossn_kr_dim(rho, b, Bt);
        Json r;
        r["C"] = s.C;
        r["B"] = encode_vector(s.B);
        return r;
    }
};
 
/** lossn_kr_asy(rho, b, C): critical-loading blocking, alpha, sigma, delta. */
struct OpLossnKrAsy {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> rho = a.vector<T>("rho");
        const std::vector<T> b = a.vector<T>("b");
        const T C = a.number<T>("C");
        a.done();
        const lossn::KrAsyResult<T> s = lossn::lossn_kr_asy(rho, b, C);
        Json r;
        r["B"] = encode_vector(s.B);
        r["alpha"] = encode_scalar(s.alpha);
        r["sigma"] = encode_scalar(s.sigma);
        r["delta"] = encode_scalar(s.delta);
        return r;
    }
};
 
/** lossn_asy(nu, A, C): the Kelly limit. */
struct OpLossnAsy {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> nu = a.vector<T>("nu");
        const Matrix<T> A = a.matrix<T>("A");
        const std::vector<T> C = a.vector<T>("C");
        a.done();
        const lossn::LossnAsyResult<T> s = lossn::lossn_asy(nu, A, C);
        Json r;
        r["QLen"] = encode_vector(s.QLen);
        r["Loss"] = encode_vector(s.Loss);
        r["a"] = encode_vector(s.a);
        r["niter"] = s.niter;
        return r;
    }
};
 
/**
 * The options struct lossn_links and lossn_revsens forward to lossn_mci:
 * {samples, alpha, gamma, seed}. The seed is required only when the method is
 * 'mci', for the reason OpLossnMci gives.
 */
template <class T>
void lossn_parse_mci_options(const Json* oj, const std::string& fn, bool isMci,
                             lossn::LossnMciOptions<T>& o, std::uint64_t& seed) {
    bool haveSeed = false;
    if (oj) {
        Args op(*oj, fn + " options");
        if (op.has("samples")) o.samples = op.count("samples");
        if (op.has("alpha")) o.alpha = op.scalar<double>("alpha", 0.0);
        if (op.has("gamma")) o.gamma = op.vector<T>("gamma");
        if (op.has("seed")) {
            seed = static_cast<std::uint64_t>(op.uinteger("seed"));
            haveSeed = true;
        }
        op.done();
    }
    if (isMci && !haveSeed)
        throw InputError(fn + ": options.seed is required for method 'mci'; the port draws from a "
                              "seeded generator and has no session RNG to fall back on");
}
 
/** lossn_links(nu, A, C, method, options): per-link blocking E and occupancy Occ. */
struct OpLossnLinks {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> nu = a.vector<T>("nu");
        const Matrix<T> A = a.matrix<T>("A");
        const std::vector<T> C = a.vector<T>("C");
        const std::string method = a.text("method", "exact");
        const Json* oj = a.opt("options");
        a.done();
        lossn::LossnMciOptions<T> o;
        std::uint64_t seed = 0;
        lossn_parse_mci_options<T>(oj, "lossn_links", lossn::detail::lossn_lower(method) == "mci", o,
                                   seed);
        const lossn::LossnLinksResult<T> s = lossn::lossn_links(nu, A, C, method, o, seed);
        Json r;
        r["E"] = encode_matrix(s.E);
        r["Occ"] = encode_vector(s.Occ);
        r["Loss"] = encode_vector(s.Loss);
        r["QLen"] = encode_vector(s.QLen);
        return r;
    }
};
 
/** lossn_revsens(nu, A, C, w, method, options): revenue sensitivities and implied costs. */
struct OpLossnRevsens {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> nu = a.vector<T>("nu");
        const Matrix<T> A = a.matrix<T>("A");
        const std::vector<T> C = a.vector<T>("C");
        const std::vector<T> w = a.vector<T>("w");
        const std::string method = a.text("method", "exact");
        const Json* oj = a.opt("options");
        a.done();
        lossn::LossnMciOptions<T> o;
        std::uint64_t seed = 0;
        lossn_parse_mci_options<T>(oj, "lossn_revsens", lossn::detail::lossn_lower(method) == "mci",
                                   o, seed);
        const lossn::LossnRevsensResult<T> s = lossn::lossn_revsens(nu, A, C, w, method, o, seed);
        Json r;
        r["dW"] = encode_vector(s.dW);
        r["c"] = encode_vector(s.c);
        r["W"] = encode_scalar(s.W);
        r["clink"] = encode_matrix(s.clink);
        if (s.hasCI) {
            Json ci;
            ci["revsens"] = encode_vector(s.revsens);
            ci["revsensCI"] = encode_matrix(s.revsensCI);
            ci["level"] = encode_scalar(s.level);
            r["ci"] = ci;
        } else {
            r["ci"] = Json::array();
        }
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// dqsys_bernoulli1
// ---------------------------------------------------------------------------
 
/**
 * dqsys_bernoulli1(b, p[, L]). Without L the buffer is unbounded and b, p are
 * scalars; MATLAB then returns pmf and arrivalPmf as function handles, which
 * have no JSON form and are left out (the closed-form law is Geo/Geo/1, see
 * dqsys_geogeo1). With L every field is a plain value.
 */
struct OpDqsysBernoulli1 {
    template <class T>
    static Json run(Args& a) {
        Json r;
        if (!a.has("L")) {
            const T b = a.number<T>("b");
            const T p = a.number<T>("p");
            a.done();
            const dqsys::Bernoulli1InfiniteResult<T> s = dqsys::dqsys_bernoulli1(b, p);
            r["capacity"] = "Inf";
            r["arrivalProb"] = encode_scalar(s.arrivalProb);
            r["serviceProb"] = encode_scalar(s.serviceProb);
            r["emptyProb"] = encode_scalar(s.emptyProb);
            r["utilization"] = encode_scalar(s.utilization);
            r["throughput"] = encode_scalar(s.throughput);
            r["lossProb"] = encode_scalar(s.lossProb);
            r["meanQueueLength"] = encode_scalar(s.meanQueueLength);
            r["meanWaitingQueue"] = encode_scalar(s.meanWaitingQueue);
            r["meanSojournTime"] = encode_scalar(s.meanSojournTime);
            r["meanWaitingTime"] = encode_scalar(s.meanWaitingTime);
            r["normConst"] = encode_scalar(s.normConst);
            r["analyzer"] = s.analyzer;
            return r;
        }
        const std::vector<T> b = a.vector<T>("b");
        const std::vector<T> p = a.vector<T>("p");
        const std::size_t L = a.count("L");
        a.done();
        const dqsys::Bernoulli1Result<T> s = dqsys::dqsys_bernoulli1(b, p, L);
        r["capacity"] = encode_count(s.capacity);
        r["arrivalProb"] = encode_vector(s.arrivalProb);
        r["serviceProb"] = encode_vector(s.serviceProb);
        r["pmf"] = encode_vector(s.pmf);
        r["arrivalPmf"] = encode_vector(s.arrivalPmf);
        r["emptyProb"] = encode_scalar(s.emptyProb);
        r["utilization"] = encode_scalar(s.utilization);
        r["throughput"] = encode_scalar(s.throughput);
        r["lossProb"] = encode_scalar(s.lossProb);
        r["meanQueueLength"] = encode_scalar(s.meanQueueLength);
        r["meanWaitingQueue"] = encode_scalar(s.meanWaitingQueue);
        r["meanSojournTime"] = encode_scalar(s.meanSojournTime);
        r["meanWaitingTime"] = encode_scalar(s.meanWaitingTime);
        r["normConst"] = encode_scalar(s.normConst);
        r["analyzer"] = s.analyzer;
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// MAP and phase-type helpers
// ---------------------------------------------------------------------------
 
struct OpMap2ph {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = map_from_json<T>(a.get("MAPIN"), "map2ph: MAPIN");
        a.done();
        const mam::PhType<T> ph = mam::map2ph(m);
        Json r;
        r["alpha"] = encode_vector(ph.alpha);
        r["T"] = encode_matrix(ph.subgen);
        return r;
    }
};
 
struct OpMap2mmpp {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = map_from_json<T>(a.get("MAP"), "map2mmpp: MAP");
        a.done();
        const mam::Map2mmppResult<T> s = mam::map2mmpp(m);
        Json r;
        r["Q"] = encode_matrix(s.Q);
        r["LAMBDA"] = encode_matrix(s.LAMBDA);
        return r;
    }
};
 
struct OpAutocov {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> X = a.vector<T>("X");
        a.done();
        Json r;
        r["acv"] = encode_vector(trace::autocov(X));
        return r;
    }
};
 
/** aph_fit(e1, e2, e3[, nmax]); tol is the port's moment-matching tolerance and needs nmax. */
struct OpAphFit {
    template <class T>
    static Json run(Args& a) {
        const T e1 = a.number<T>("e1");
        const T e2 = a.number<T>("e2");
        const T e3 = a.number<T>("e3");
        const bool hn = a.has("nmax"), ht = a.has("tol");
        if (ht && !hn)
            throw InputError("aph_fit: tol follows nmax in the port's overloads; give nmax too");
        mam::AphFitResult<T> s;
        if (ht) {
            const unsigned nmax = static_cast<unsigned>(a.count("nmax"));
            const T tol = a.number<T>("tol");
            a.done();
            s = mam::aph_fit(e1, e2, e3, nmax, tol);
        } else if (hn) {
            const unsigned nmax = static_cast<unsigned>(a.count("nmax"));
            a.done();
            s = mam::aph_fit(e1, e2, e3, nmax);
        } else {
            a.done();
            s = mam::aph_fit(e1, e2, e3);
        }
        Json r;
        r["APH"] = encode_map(s.aph);
        r["isexact"] = s.isexact;
        return r;
    }
};
 
/** m3pp22_fitc_approx_cov: FIT is the MMAP cell {D0, D1, D1_1, D1_2}. */
struct OpM3pp22FitcApproxCov {
    template <class T>
    static Json run(Args& a) {
        const T av = a.number<T>("a");
        const T bt1 = a.number<T>("bt1");
        const T bt2 = a.number<T>("bt2");
        const T binf = a.number<T>("binf");
        const T m3t2 = a.number<T>("m3t2");
        const T t1 = a.number<T>("t1");
        const T t2 = a.number<T>("t2");
        const std::vector<T> ai = a.vector<T>("ai");
        const T st3 = a.number<T>("st3");
        const T t3 = a.number<T>("t3");
        a.done();
        const mam::M3pp22FitcCovResult<T> s =
            mam::m3pp22_fitc_approx_cov(av, bt1, bt2, binf, m3t2, t1, t2, ai, st3, t3);
        Json fit = Json::array();
        fit.push_back(encode_matrix(s.mmap.D0));
        fit.push_back(encode_matrix(s.mmap.D1));
        for (const auto& d : s.mmap.Dc) fit.push_back(encode_matrix(d));
        Json r;
        r["FIT"] = fit;
        return r;
    }
};
 
struct OpPhMultisets {
    template <class T>
    static Json run(Args& a) {
        const std::size_t p = a.count("p");
        const std::size_t k = a.count("k");
        a.done();
        const std::vector<std::vector<int> > M = mam::ph_multisets(p, k);
        Json rows = Json::array();
        for (const auto& row : M) rows.push_back(encode_ints(row));
        Json r;
        r["M"] = rows;
        return r;
    }
};
 
struct OpMmdpIsfeasible {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const Matrix<T> R = a.matrix<T>("R");
        a.done();
        Json r;
        r["bool"] = mam::mmdp_isfeasible(Q, R);
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// pfqn helpers
// ---------------------------------------------------------------------------
 
struct OpInfradiusH {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> x = a.vector<T>("x");
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<T> N = a.vector<T>("N");
        const Matrix<T> alpha = a.matrix<T>("alpha");
        a.done();
        Json r;
        r["y"] = encode_scalar(pfqn::infradius_h(x, L, N, alpha));
        return r;
    }
};
 
/**
 * pas_placement(H) -> P. MATLAB's second output is the closure placeable(x),
 * which has no JSON form; give the remaining per-class counts x to have it
 * evaluated, as 0-based class indices.
 */
struct OpPasPlacement {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> H = a.matrix<T>("H");
        const bool hx = a.has("x");
        const std::vector<T> x = hx ? a.vector<T>("x") : std::vector<T>();
        a.done();
        const pfqn::PasPlacement<T> s = pfqn::pas_placement(H);
        Json r;
        r["P"] = encode_matrix(s.P);
        if (hx) r["placeable"] = encode_sizes(s.placeable(x));
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// mapqn: MAP queueing networks
// ---------------------------------------------------------------------------
 
struct OpMapqnAmva {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> mu = a.vector<T>("mu");
        const std::vector<Matrix<T> > D0s = matrices_from_json<T>(a.get("D0s"), "mapqn_amva: D0s");
        const std::vector<Matrix<T> > D1s = matrices_from_json<T>(a.get("D1s"), "mapqn_amva: D1s");
        const std::vector<int> N = a.ints("N");
        a.done();
        const mapqn::MapqnAmvaResult<T> s = mapqn::mapqn_amva(mu, D0s, D1s, N);
        Json r;
        r["X"] = encode_vector(s.X);
        r["Qq"] = encode_vector(s.Qq);
        r["U"] = encode_vector(s.U);
        r["pik"] = encode_vector(s.pi);
        r["ES"] = encode_vector(s.ES);
        return r;
    }
};
 
/** sense is 'max' or 'min'; absent, the port's own default overload is called. */
mapqn::MapqnSense mapqn_sense(const std::string& s, const char* fn) {
    if (s == "max" || s == "Max" || s == "MAX") return mapqn::MapqnSense::Max;
    if (s == "min" || s == "Min" || s == "MIN") return mapqn::MapqnSense::Min;
    throw InputError(std::string(fn) + ": sense must be 'max' or 'min', got '" + s + "'");
}
 
/** The MATLAB params fields every mapqn bound shares: M, N, K, mu, v, r. */
template <class T, class P>
void read_mapqn_core(Args& p, P& out, const char* fn) {
    out.M = p.required_integer("M");
    out.N = p.required_integer("N");
    out.K = p.ints("K");
    out.mu = matrices_from_json<T>(p.get("mu"), where_of(fn, "params.mu"));
    out.v = matrices_from_json<T>(p.get("v"), where_of(fn, "params.v"));
    out.r = p.matrix<T>("r");
    p.unsupported("verbose", "the port prints nothing, so there is no verbosity to set.");
    p.unsupported("lpAlgorithm", "the port solves the LP with its own simplex and offers no "
                                 "algorithm choice.");
}
 
const Json& params_object(Args& a, const char* fn) {
    const Json& j = a.get("params");
    if (!j.is_object()) throw InputError(std::string(fn) + ": params must be an object");
    return j;
}
 
template <class R>
Json mapqn_lp_status(const R& s) {
    Json r;
    r["objective"] = encode_scalar(s.objective);
    r["ok"] = s.ok;
    r["status"] = s.status;
    return r;
}
 
/**
 * mapqn_bnd_lr and mapqn_bnd_qr share the signature
 * (params, objective_queue, objective_phase[, sense]) and the result
 * {objective, U, IT, Q}; queues and phases are 0-based. The port has no
 * default for objective_queue and objective_phase, so both are required.
 */
#define LINE_MAPQN_UIQ_OP(OPNAME, FN, RESULT)                                             \
    struct OPNAME {                                                                       \
        template <class T>                                                                \
        static Json run(Args& a) {                                                        \
            Args p(params_object(a, #FN), #FN " params");                                 \
            mapqn::MapqnParams<T> prm;                                                    \
            read_mapqn_core<T>(p, prm, #FN);                                              \
            p.done();                                                                     \
            const int q = a.required_integer("objective_queue");                          \
            const int ph = a.required_integer("objective_phase");                         \
            const bool hs = a.has("sense");                                               \
            const std::string sense = a.text("sense", "");                                \
            a.done();                                                                     \
            const mapqn::RESULT<T> s = hs ? mapqn::FN(prm, q, ph, mapqn_sense(sense, #FN))\
                                          : mapqn::FN(prm, q, ph);                        \
            Json res = mapqn_lp_status(s);                                                \
            res["U"] = encode_matrix(s.U);                                                \
            res["IT"] = encode_matrix(s.IT);                                              \
            res["Q"] = encode_matrix(s.Q);                                                \
            Json r;                                                                       \
            r["result"] = res;                                                            \
            r["x"] = encode_vector(s.x);                                                  \
            r["fval"] = encode_scalar(s.objective);                                       \
            return r;                                                                     \
        }                                                                                 \
    };
 
LINE_MAPQN_UIQ_OP(OpMapqnBndLr, mapqn_bnd_lr, MapqnBndLrResult)
LINE_MAPQN_UIQ_OP(OpMapqnBndQr, mapqn_bnd_qr, MapqnBndQrResult)
#undef LINE_MAPQN_UIQ_OP
 
/**
 * mapqn_bnd_qr_ld and mapqn_bnd_qr_delay: (params, objective_queue,
 * objective_phase, objective_n[, sense]) with the load-dependent scalings
 * params.alpha (M x N); the delay variant also reads params.Z and params.D1.
 */
template <class T>
Json mapqn_qr_ld_like(Args& a, const char* fn, bool delay) {
    Args p(params_object(a, fn), std::string(fn) + " params");
    mapqn::MapqnParams<T> prm;
    read_mapqn_core<T>(p, prm, fn);
    prm.alpha = p.matrix<T>("alpha");
    if (delay) {
        prm.Z = p.number<T>("Z");
        prm.D1 = p.number<T>("D1");
    }
    p.done();
    const int q = a.required_integer("objective_queue");
    const int ph = a.required_integer("objective_phase");
    const int n = a.required_integer("objective_n");
    const bool hs = a.has("sense");
    const std::string sense = a.text("sense", "");
    a.done();
    mapqn::MapqnQrResult<T> s;
    if (delay)
        s = hs ? mapqn::mapqn_bnd_qr_delay(prm, q, ph, n, mapqn_sense(sense, fn))
               : mapqn::mapqn_bnd_qr_delay(prm, q, ph, n);
    else
        s = hs ? mapqn::mapqn_bnd_qr_ld(prm, q, ph, n, mapqn_sense(sense, fn))
               : mapqn::mapqn_bnd_qr_ld(prm, q, ph, n);
    Json res = mapqn_lp_status(s);
    res["p2marginals"] = encode_matrices(s.p2marginals);
    Json r;
    r["result"] = res;
    r["x"] = encode_vector(s.x);
    r["fval"] = encode_scalar(s.objective);
    return r;
}
 
struct OpMapqnBndQrLd {
    template <class T>
    static Json run(Args& a) { return mapqn_qr_ld_like<T>(a, "mapqn_bnd_qr_ld", false); }
};
 
struct OpMapqnBndQrDelay {
    template <class T>
    static Json run(Args& a) { return mapqn_qr_ld_like<T>(a, "mapqn_bnd_qr_delay", true); }
};
 
/**
 * mapqn_bnd_lr_mva(params, objective_queue, objective_level[, sense[, objective_var]]):
 * params {M, N, K, muM, muMAP, v, r}; objective_level -1 is the sum over levels;
 * objective_var 'UN' or 'QN'.
 */
struct OpMapqnBndLrMva {
    template <class T>
    static Json run(Args& a) {
        const char* fn = "mapqn_bnd_lr_mva";
        Args p(params_object(a, fn), "mapqn_bnd_lr_mva params");
        mapqn::LrMvaParams<T> prm;
        prm.M = p.required_integer("M");
        prm.N = p.required_integer("N");
        prm.K = p.required_integer("K");
        prm.muM = p.vector<T>("muM");
        prm.muMAP = p.matrix<T>("muMAP");
        prm.v = p.matrix<T>("v");
        prm.r = p.matrix<T>("r");
        p.unsupported("verbose", "the port prints nothing, so there is no verbosity to set.");
        p.unsupported("lpAlgorithm", "the port solves the LP with its own simplex and offers no "
                                     "algorithm choice.");
        p.done();
        const int q = a.required_integer("objective_queue");
        const int lev = a.required_integer("objective_level");
        const bool hs = a.has("sense"), hv = a.has("objective_var");
        if (hv && !hs)
            throw InputError("mapqn_bnd_lr_mva: objective_var follows sense in the port's "
                             "signature; give sense too");
        const std::string sense = a.text("sense", "");
        const std::string var = a.text("objective_var", "");
        a.done();
        mapqn::MapqnBndLrMvaResult<T> s;
        if (hv) {
            mapqn::MapqnObjectiveVar ov;
            if (var == "UN") ov = mapqn::MapqnObjectiveVar::UN;
            else if (var == "QN") ov = mapqn::MapqnObjectiveVar::QN;
            else throw InputError("mapqn_bnd_lr_mva: objective_var must be 'UN' or 'QN', got '" + var + "'");
            s = mapqn::mapqn_bnd_lr_mva(prm, q, lev, mapqn_sense(sense, fn), ov);
        } else if (hs) {
            s = mapqn::mapqn_bnd_lr_mva(prm, q, lev, mapqn_sense(sense, fn));
        } else {
            s = mapqn::mapqn_bnd_lr_mva(prm, q, lev);
        }
        Json res = mapqn_lp_status(s);
        res["UN"] = encode_matrix(s.UN);
        res["QN"] = encode_matrix(s.QN);
        Json r;
        r["result"] = res;
        return r;
    }
};
 
std::vector<std::vector<int> > int_rows_from_json(const Json& j, const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected a list of integer rows");
    std::vector<std::vector<int> > out;
    for (std::size_t k = 0; k < j.size(); ++k) out.push_back(int_vector_from_json(j[k], at(where, k)));
    return out;
}
 
/**
 * mapqn_qr_bounds_bas(params, objective_queue[, sense]). There is no MATLAB
 * function; the keys are the port's QrBasParams fields (those of the JAR and
 * native Python parameter objects), all indices 0-based.
 */
struct OpMapqnQrBoundsBas {
    template <class T>
    static Json run(Args& a) {
        const char* fn = "mapqn_qr_bounds_bas";
        Args p(params_object(a, fn), "mapqn_qr_bounds_bas params");
        mapqn::QrBasParams<T> prm;
        prm.M = p.required_integer("M");
        prm.N = p.required_integer("N");
        prm.f = p.required_integer("f");
        prm.F = p.ints("F");
        prm.K = p.ints("K");
        prm.mu = matrices_from_json<T>(p.get("mu"), where_of(fn, "params.mu"));
        prm.v = matrices_from_json<T>(p.get("v"), where_of(fn, "params.v"));
        prm.r = p.matrix<T>("r");
        prm.MR = p.required_integer("MR");
        prm.BB = int_rows_from_json(p.get("BB"), where_of(fn, "params.BB"));
        prm.MM = int_rows_from_json(p.get("MM"), where_of(fn, "params.MM"));
        prm.ZZ = p.ints("ZZ");
        prm.ZM = p.required_integer("ZM");
        prm.MM1 = int_rows_from_json(p.get("MM1"), where_of(fn, "params.MM1"));
        p.done();
        const int q = a.required_integer("objective_queue");
        const bool hs = a.has("sense");
        const std::string sense = a.text("sense", "");
        a.done();
        const mapqn::QrBasResult<T> s = hs ? mapqn::mapqn_qr_bounds_bas(prm, q, mapqn_sense(sense, fn))
                                           : mapqn::mapqn_qr_bounds_bas(prm, q);
        Json r = mapqn_lp_status(s);
        r["U"] = encode_vector(s.U);
        r["occupancy"] = encode_vector(s.occupancy);
        r["e"] = encode_matrix(s.e);
        r["x"] = encode_vector(s.x);
        return r;
    }
};
 
/** mapqn_qr_bounds_rsrd(params, objective_queue[, sense]); params as QrRsrdParams, alpha optional. */
struct OpMapqnQrBoundsRsrd {
    template <class T>
    static Json run(Args& a) {
        const char* fn = "mapqn_qr_bounds_rsrd";
        Args p(params_object(a, fn), "mapqn_qr_bounds_rsrd params");
        mapqn::QrRsrdParams<T> prm;
        prm.M = p.required_integer("M");
        prm.N = p.required_integer("N");
        prm.F = p.ints("F");
        prm.K = p.ints("K");
        prm.mu = matrices_from_json<T>(p.get("mu"), where_of(fn, "params.mu"));
        prm.v = matrices_from_json<T>(p.get("v"), where_of(fn, "params.v"));
        prm.r = p.matrix<T>("r");
        if (p.has("alpha")) prm.alpha = vectors_from_json<T>(p.get("alpha"), where_of(fn, "params.alpha"));
        p.done();
        const int q = a.required_integer("objective_queue");
        const bool hs = a.has("sense");
        const std::string sense = a.text("sense", "");
        a.done();
        const mapqn::QrRsrdResult<T> s = hs ? mapqn::mapqn_qr_bounds_rsrd(prm, q, mapqn_sense(sense, fn))
                                            : mapqn::mapqn_qr_bounds_rsrd(prm, q);
        Json r = mapqn_lp_status(s);
        r["U"] = encode_vector(s.U);
        r["Ueff"] = encode_vector(s.Ueff);
        r["pb"] = encode_vector(s.pb);
        r["x"] = encode_vector(s.x);
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// retrieval: caches with a retrieval (fetching) system
// ---------------------------------------------------------------------------
 
/**
 * retrieval_metrics and retrieval_mva share (m, lambda, eta, gamma) and the
 * [pmiss, phit, pdh] return list.
 */
#define LINE_RETRIEVAL_EXACT_OP(OPNAME, FN)                                       \
    struct OPNAME {                                                               \
        template <class T>                                                        \
        static Json run(Args& a) {                                                \
            const std::vector<int> m = a.ints("m");                               \
            const std::vector<T> lambda = a.vector<T>("lambda");                  \
            const Matrix<T> eta = a.matrix<T>("eta");                             \
            const Matrix<T> gamma = a.matrix<T>("gamma");                         \
            a.done();                                                             \
            const auto s = retrieval::FN(m, lambda, eta, gamma);                  \
            Json r;                                                               \
            r["pmiss"] = encode_vector(s.pmiss);                                  \
            r["phit"] = encode_matrix(s.phit);                                    \
            r["pdh"] = encode_matrix(s.pdh);                                      \
            return r;                                                             \
        }                                                                         \
    };
 
LINE_RETRIEVAL_EXACT_OP(OpRetrievalMetrics, retrieval_metrics)
LINE_RETRIEVAL_EXACT_OP(OpRetrievalMva, retrieval_mva)
#undef LINE_RETRIEVAL_EXACT_OP
 
/** The trailing (max_iter, tol) of the MATLAB FPI functions; each key sets only its own field. */
retrieval::FpiOptions read_fpi_options(Args& a) {
    retrieval::FpiOptions o;
    if (a.has("max_iter")) o.max_iter = a.count("max_iter");
    if (a.has("tol")) o.tol = a.number<double>("tol");
    return o;
}
 
struct OpRetrievalFpi {
    template <class T>
    static Json run(Args& a) {
        const std::vector<int> m = a.ints("m");
        const std::vector<T> lambda = a.vector<T>("lambda");
        const Matrix<T> eta = a.matrix<T>("eta");
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const retrieval::FpiOptions o = read_fpi_options(a);
        a.done();
        const retrieval::RetrievalFpiResult<T> s = retrieval::retrieval_fpi(m, lambda, eta, gamma, o);
        Json r;
        r["pmiss"] = encode_vector(s.pmiss);
        r["phit"] = encode_matrix(s.phit);
        r["pdh"] = encode_matrix(s.pdh);
        return r;
    }
};
 
retrieval::RetrievalStationType retrieval_station_type(const Json& j, const std::string& where) {
    if (!j.is_string()) throw InputError(where + ": a station type is a string");
    const std::string s = j.get<std::string>();
    if (s == "IS") return retrieval::RetrievalStationType::IS;
    if (s == "PS") return retrieval::RetrievalStationType::PS;
    if (s == "SIRO") return retrieval::RetrievalStationType::SIRO;
    if (s == "FCFS") return retrieval::RetrievalStationType::FCFS;
    if (s == "LCFSPR") return retrieval::RetrievalStationType::LCFSPR;
    throw InputError(where + ": station type must be IS, PS, SIRO, FCFS or LCFSPR, got '" + s + "'");
}
 
/**
 * retrieval_fpi_latency(m, lambda, gamma, alpha, T, R, station_type).
 * alpha{s} is the (n x f) matrix of per-item entry vectors of station s, T{s}
 * the list of its n per-item subgenerators, and R the list of n routing
 * matrices, (S+1) x (S+1) with index 0 = outside (MATLAB's third dimension
 * becomes the list index). max_iter and tol are the port's FpiOptions; MATLAB
 * fixes them inside the call.
 */
struct OpRetrievalFpiLatency {
    template <class T>
    static Json run(Args& a) {
        const char* fn = "retrieval_fpi_latency";
        const std::vector<int> m = a.ints("m");
        const std::vector<T> lambda = a.vector<T>("lambda");
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<Matrix<T> > alpha = matrices_from_json<T>(a.get("alpha"), where_of(fn, "alpha"));
        const std::vector<std::vector<Matrix<T> > > Tsub =
            matrix_lists_from_json<T>(a.get("T"), where_of(fn, "T"));
        const std::vector<Matrix<T> > R = matrices_from_json<T>(a.get("R"), where_of(fn, "R"));
        const Json& ty = a.get("station_type");
        const retrieval::FpiOptions o = read_fpi_options(a);
        a.done();
        const std::size_t S = alpha.size();
        if (Tsub.size() != S)
            throw InputError("retrieval_fpi_latency: alpha and T disagree on the station count");
        std::vector<retrieval::RetrievalStationType> types;
        if (ty.is_string()) {
            types.push_back(retrieval_station_type(ty, where_of(fn, "station_type")));
        } else if (ty.is_array()) {
            for (std::size_t s = 0; s < ty.size(); ++s)
                types.push_back(retrieval_station_type(ty[s], at(where_of(fn, "station_type"), s)));
        } else {
            throw InputError("retrieval_fpi_latency: station_type is a list of strings");
        }
        if (types.size() != S)
            throw InputError("retrieval_fpi_latency: station_type and alpha disagree on the station count");
        std::vector<retrieval::RetrievalStationPH<T> > st(S);
        for (std::size_t s = 0; s < S; ++s) {
            st[s].alpha = alpha[s];
            st[s].sub = Tsub[s];
            st[s].type = types[s];
        }
        const retrieval::RetrievalFpiLatencyResult<T> res =
            retrieval::retrieval_fpi_latency(m, lambda, gamma, st, R, o);
        Json r;
        r["Z"] = encode_scalar(res.Z);
        r["d"] = encode_vector(res.d);
        r["phi"] = encode_vector(res.phi);
        r["pi0"] = encode_vector(res.pi0);
        return r;
    }
};
 
struct OpRetrievalNc {
    template <class T>
    static Json run(Args& a) {
        const std::vector<int> v = a.ints("v");
        const std::vector<int> m = a.ints("m");
        const std::vector<T> lambda = a.vector<T>("lambda");
        const Matrix<T> eta = a.matrix<T>("eta");
        const Matrix<T> gamma = a.matrix<T>("gamma");
        a.done();
        Json r;
        r["E"] = encode_scalar(retrieval::retrieval_nc(v, m, lambda, eta, gamma));
        return r;
    }
};
 
/** retrieval_rayint(gamma, m): the matrix form. n and nquad belong to the gfun form, whose gamma is a callable. */
struct OpRetrievalRayint {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> gamma = a.matrix<T>("gamma");
        const std::vector<int> m = a.ints("m");
        a.unsupported("n", "it selects the form whose gamma is a function handle gfun(t), which "
                           "has no JSON form; give gamma as an (n x h) matrix instead.");
        a.unsupported("nquad", "it is the quadrature size of the gfun form, whose gamma is a "
                               "function handle with no JSON form.");
        a.done();
        const retrieval::RetrievalRayintResult<T> s = retrieval::retrieval_rayint(gamma, m);
        Json out;
        out["xi"] = encode_vector(s.xi);
        out["logdetSigma"] = encode_scalar(s.logdet_sigma);
        out["iter"] = encode_count(s.iterations);
        out["relerrEst"] = s.relerr_est;
        out["method"] = s.method;
        Json r;
        r["E"] = encode_scalar(s.e);
        r["logE"] = encode_scalar(s.log_e);
        r["out"] = out;
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// sum: the summation method
// ---------------------------------------------------------------------------
 
/** mi(i): a server count, or "Inf" for a delay station. */
std::vector<sum::Servers> servers_from_json(const Json& j, const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected a list of server counts");
    std::vector<sum::Servers> out;
    for (std::size_t k = 0; k < j.size(); ++k) {
        const Json& e = j[k];
        if (e.is_string()) {
            const std::string s = e.get<std::string>();
            if (s == "Inf" || s == "inf") {
                out.push_back(sum::Servers::inf());
                continue;
            }
        }
        const long m = static_cast<long>(size_from_json(e, at(where, k)));
        if (m < 1) throw InputError(at(where, k) + ": a station needs at least one server");
        out.push_back(sum::Servers::of(m));
    }
    return out;
}
 
/** N(r) as a long; "Inf" is accepted only where allowInf (an open class in sum_closing). */
std::vector<long> populations_from_json(const Json& j, const std::string& where,
                                        const std::vector<bool>* allowInf) {
    if (!j.is_array()) throw InputError(where + ": expected a list of populations");
    std::vector<long> out;
    for (std::size_t k = 0; k < j.size(); ++k) {
        const Json& e = j[k];
        if (e.is_string() && (e.get<std::string>() == "Inf" || e.get<std::string>() == "inf")) {
            if (!allowInf || k >= allowInf->size() || !(*allowInf)[k])
                throw InputError(at(where, k) + ": Inf marks an open class (lambda0 > 0) only");
            out.push_back(0);  // the port ignores the N entry of an open class
            continue;
        }
        out.push_back(static_cast<long>(size_from_json(e, at(where, k))));
    }
    return out;
}
 
sum::SumOptions read_sum_options(Args& a) {
    sum::SumOptions o;
    if (a.has("tol")) o.tol = a.number<double>("tol");
    if (a.has("maxiter")) o.maxiter = a.count("maxiter");
    return o;
}
 
/**
 * sum_closed(L, N, Z, mi, scv, tol, maxiter). MATLAB defaults Z, mi and scv;
 * the port takes all three as required parameters, so they are required here.
 */
struct OpSumClosed {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<long> N = populations_from_json(a.get("N"), "sum_closed: N", nullptr);
        const std::vector<T> Z = a.vector<T>("Z");
        const std::vector<sum::Servers> mi = servers_from_json(a.get("mi"), "sum_closed: mi");
        const Matrix<T> scv = a.matrix<T>("scv");
        const sum::SumOptions o = read_sum_options(a);
        a.done();
        const sum::SumClosedResult<T> s = sum::sum_closed(L, N, Z, mi, scv, o);
        Json r;
        r["XN"] = encode_vector(s.XN);
        r["QN"] = encode_matrix(s.QN);
        r["UN"] = encode_matrix(s.UN);
        r["RN"] = encode_matrix(s.RN);
        r["it"] = encode_count(s.it);
        return r;
    }
};
 
/**
 * sum_closing(lambda0, scva, L, mi, scv, N, Z, Kclosed, tol, maxiter). As for
 * sum_closed, every vector the port takes is required. N(r) may be "Inf" only
 * for an open class, whose entry the port ignores.
 */
struct OpSumClosing {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda0 = a.vector<T>("lambda0");
        const std::vector<T> scva = a.vector<T>("scva");
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<sum::Servers> mi = servers_from_json(a.get("mi"), "sum_closing: mi");
        const Matrix<T> scv = a.matrix<T>("scv");
        std::vector<bool> open(lambda0.size());
        for (std::size_t r = 0; r < lambda0.size(); ++r)
            open[r] = num_traits<T>::to_double(lambda0[r]) > 0.0;
        const std::vector<long> N = populations_from_json(a.get("N"), "sum_closing: N", &open);
        const std::vector<T> Z = a.vector<T>("Z");
        sum::ClosingOptions o;
        if (a.has("Kclosed")) o.Kclosed = static_cast<long>(a.count("Kclosed"));
        o.sum = read_sum_options(a);
        a.done();
        const sum::SumClosingResult<T> s = sum::sum_closing(lambda0, scva, L, mi, scv, N, Z, o);
        Json r;
        r["XN"] = encode_vector(s.XN);
        r["QN"] = encode_matrix(s.QN);
        r["UN"] = encode_matrix(s.UN);
        r["RN"] = encode_matrix(s.RN);
        r["TN"] = encode_vector(s.TN);
        r["it"] = encode_count(s.it);
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// wf: workflow pattern detectors (JAR reference; no MATLAB counterpart)
// ---------------------------------------------------------------------------
//
// Each detector takes the (nedges x 3) linkMatrix [from, to, probability] and
// the node-id lists the JAR method takes (serviceNodes, joinNodes, routerNodes,
// forkNodes; node ids as they appear in the link matrix), and returns the
// detected patterns together with the statistics map the JAR's get*Stats
// method computes over them, keyed as the JAR keys it.
 
struct OpWfBranchDetector {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> link = a.matrix<T>("linkMatrix");
        const std::vector<int> serviceNodes = a.ints("serviceNodes");
        const std::vector<int> joinNodes = a.ints("joinNodes");
        a.done();
        const std::vector<wf::BranchPattern<T> > p = wf::detect_branches(link, serviceNodes, joinNodes);
        Json pats = Json::array();
        for (std::size_t k = 0; k < p.size(); ++k) {
            Json e;
            e["branchNodes"] = encode_ints(p[k].branchNodes);
            e["probabilities"] = encode_vector(p[k].probabilities);
            e["forkNode"] = p[k].forkNode;
            e["joinNode"] = p[k].hasJoinNode ? Json(p[k].joinNode) : Json(nullptr);
            const wf::BranchDiversity<T> d = wf::calculate_branch_diversity(p[k]);
            Json dj;
            dj["entropy"] = encode_scalar(d.entropy);
            dj["normalizedEntropy"] = encode_scalar(d.normalizedEntropy);
            dj["gini"] = encode_scalar(d.gini);
            dj["balance"] = encode_scalar(d.balance);
            e["diversity"] = dj;
            pats.push_back(e);
        }
        const wf::BranchStats<T> s = wf::get_branch_stats(p);
        Json st;
        st["numPatterns"] = encode_count(s.numPatterns);
        st["totalBranchNodes"] = encode_count(s.totalBranchNodes);
        st["avgBranches"] = encode_scalar(s.avgBranches);
        st["maxBranches"] = encode_count(s.maxBranches);
        st["minBranches"] = encode_count(s.minBranches);
        st["avgEntropy"] = encode_scalar(s.avgEntropy);
        st["avgBalance"] = encode_scalar(s.avgBalance);
        Json r;
        r["patterns"] = pats;
        r["stats"] = st;
        return r;
    }
};
 
/** routerNodes is required; joinNodes absent calls the JAR's three-argument detectLoops. */
struct OpWfLoopDetector {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> link = a.matrix<T>("linkMatrix");
        const std::vector<int> serviceNodes = a.ints("serviceNodes");
        const std::vector<int> routerNodes = a.ints("routerNodes");
        const bool hj = a.has("joinNodes");
        const std::vector<int> joinNodes = a.ints_or_empty("joinNodes");
        a.done();
        const std::vector<int> loops = hj ? wf::detect_loops(link, serviceNodes, routerNodes, joinNodes)
                                          : wf::detect_loops(link, serviceNodes, routerNodes);
        Json lj = Json::array();
        for (std::size_t k = 0; k < loops.size(); ++k) {
            Json e;
            e["node"] = loops[k];
            const T p = wf::get_loop_probability(loops[k], link, routerNodes);
            e["probability"] = encode_scalar(p);
            const wf::ExpectedIterations<T> it = wf::get_expected_loop_iterations(p);
            e["expectedIterations"] = it.infinite ? Json("Inf") : encode_scalar(it.value);
            lj.push_back(e);
        }
        const wf::LoopStats<T> s = wf::get_loop_stats(loops, link, routerNodes);
        Json st;
        st["numLoops"] = encode_count(s.numLoops);
        st["avgLoopProbability"] = encode_scalar(s.avgLoopProbability);
        st["maxLoopProbability"] = encode_scalar(s.maxLoopProbability);
        st["minLoopProbability"] = encode_scalar(s.minLoopProbability);
        st["avgExpectedIterations"] = encode_scalar(s.avgExpectedIterations);
        st["maxExpectedIterations"] = encode_scalar(s.maxExpectedIterations);
        Json r;
        r["loops"] = lj;
        r["stats"] = st;
        return r;
    }
};
 
Json encode_int_rows(const std::vector<std::vector<int> >& v) {
    Json j = Json::array();
    for (std::size_t k = 0; k < v.size(); ++k) j.push_back(encode_ints(v[k]));
    return j;
}
 
struct OpWfParallelDetector {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> link = a.matrix<T>("linkMatrix");
        const std::vector<int> serviceNodes = a.ints("serviceNodes");
        const std::vector<int> forkNodes = a.ints("forkNodes");
        const std::vector<int> joinNodes = a.ints("joinNodes");
        a.done();
        const std::vector<std::vector<int> > p = wf::detect_parallel(link, serviceNodes, forkNodes, joinNodes);
        const wf::ParallelStats<T> s = wf::get_parallel_stats<T>(p);
        Json st;
        st["numPatterns"] = encode_count(s.numPatterns);
        st["totalParallelNodes"] = encode_count(s.totalParallelNodes);
        st["avgParallelism"] = encode_scalar(s.avgParallelism);
        st["maxParallelism"] = encode_count(s.maxParallelism);
        Json r;
        r["patterns"] = encode_int_rows(p);
        r["stats"] = st;
        return r;
    }
};
 
struct OpWfSequenceDetector {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> link = a.matrix<T>("linkMatrix");
        const std::vector<int> serviceNodes = a.ints("serviceNodes");
        a.done();
        const std::vector<std::vector<int> > p = wf::detect_sequences(link, serviceNodes);
        const wf::SequenceStats<T> s = wf::get_sequence_stats<T>(p);
        Json st;
        st["numSequences"] = encode_count(s.numSequences);
        st["totalNodes"] = encode_count(s.totalNodes);
        st["avgLength"] = encode_scalar(s.avgLength);
        st["maxLength"] = encode_count(s.maxLength);
        st["minLength"] = encode_count(s.minLength);
        Json r;
        r["sequences"] = encode_int_rows(p);
        r["stats"] = st;
        return r;
    }
};
 
// ---------------------------------------------------------------------------
// mdd / spn: decision-diagram product forms and stochastic Petri nets
// ---------------------------------------------------------------------------
 
/** An MDD.toStruct: {K, domain, root, nnodes, node}, node{k}{p}{v} the child of arc v of node id p+1. */
mdd::MddStruct mdd_struct_from_json(const Json& j, const std::string& where) {
    if (!j.is_object()) throw InputError(where + ": expected the MDD.toStruct object {K, domain, root, nnodes, node}");
    Args s(j, where);
    mdd::MddStruct m;
    m.K = s.count("K");
    m.domain = s.ints("domain");
    m.root = s.required_integer("root");
    m.nnodes = s.ints("nnodes");
    const Json& node = s.get("node");
    s.done();
    if (!node.is_array()) throw InputError(where + ".node: expected one list of nodes per level");
    for (std::size_t k = 0; k < node.size(); ++k)
        m.node.push_back(int_rows_from_json(node[k], at(where + ".node", k)));
    if (m.domain.size() != m.K || m.nnodes.size() != m.K || m.node.size() != m.K)
        throw InputError(where + ": domain, nnodes and node need one entry per level (K)");
    for (std::size_t k = 0; k < m.K; ++k)
        if (static_cast<int>(m.node[k].size()) != m.nnodes[k])
            throw InputError(at(where + ".node", k) + ": holds a different node count from nnodes");
    return m;
}
 
Json encode_mdd_struct(const mdd::MddStruct& m) {
    Json j;
    j["K"] = encode_count(m.K);
    j["domain"] = encode_ints(m.domain);
    j["root"] = m.root;
    j["nnodes"] = encode_ints(m.nnodes);
    Json node = Json::array();
    for (std::size_t k = 0; k < m.node.size(); ++k) node.push_back(encode_int_rows(m.node[k]));
    j["node"] = node;
    return j;
}
 
/** mdd_rec(mdds, g): g{l}(v+1) is g_l(v). */
struct OpMddRec {
    template <class T>
    static Json run(Args& a) {
        const mdd::MddStruct m = mdd_struct_from_json(a.get("mdds"), "mdd_rec: mdds");
        const std::vector<std::vector<T> > g = vectors_from_json<T>(a.get("g"), "mdd_rec: g");
        a.done();
        Json r;
        r["G"] = encode_scalar(mdd::mdd_rec(m, g));
        return r;
    }
};
 
/** mdd_closedqn(mu, P, servers, N[, options]); servers(i) may be "Inf" for a delay. */
struct OpMddClosedqn {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> mu = a.vector<T>("mu");
        const Matrix<T> P = a.matrix<T>("P");
        const std::vector<double> servers = a.vector<double>("servers");
        const int N = a.required_integer("N");
        if (const Json* o = a.opt("options")) {
            Args op(*o, "mdd_closedqn options");
            op.unsupported("verbose", "the port prints no storage summary; the same figures are "
                                      "returned in stats.");
            op.unsupported("mdd", "it is a live MDD object to reuse, which has no JSON form; the "
                                  "reachable set is rebuilt instead.");
            op.done();
        }
        a.done();
        const mdd::MddClosedQnResult<T> s = mdd::mdd_closedqn(mu, P, servers, N);
        const mdd::MddStats st = s.mdd.stats();
        Json sj;
        sj["levels"] = encode_count(st.levels);
        sj["nodesPerLevel"] = encode_ints(st.nodes_per_level);
        sj["numNodes"] = st.num_nodes;
        sj["liveNodes"] = st.num_nodes;
        sj["tableNodes"] = st.table_nodes;
        sj["numStates"] = st.num_states;
        sj["mddInts"] = st.mdd_ints;
        sj["explicitInts"] = st.explicit_ints;
        sj["compression"] = st.compression();
        Json tj;
        tj["reach"] = s.time_reach;
        tj["gen"] = s.time_gen;
        tj["solve"] = s.time_solve;
        tj["metrics"] = s.time_metrics;
        Json r;
        r["mdd"] = encode_mdd_struct(s.mdd.to_struct());
        r["Q"] = encode_matrix(s.Q);
        r["pi"] = encode_vector(s.pi);
        r["states"] = encode_int_rows(s.states);
        r["QLen"] = encode_vector(s.QLen);
        r["U"] = encode_vector(s.U);
        r["X"] = encode_vector(s.X);
        r["stats"] = sj;
        r["times"] = tj;
        return r;
    }
};
 
std::vector<long long> ll_vector_from_json(const Json& j, const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected a list of integers");
    std::vector<long long> out;
    for (std::size_t k = 0; k < j.size(); ++k) {
        const Json& e = j[k];
        if (!e.is_number_integer() && !e.is_number_unsigned())
            throw InputError(at(where, k) + ": expected an integer");
        out.push_back(e.get<long long>());
    }
    return out;
}
 
std::vector<std::vector<long long> > ll_rows_from_json(const Json& j, const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected a list of integer rows");
    std::vector<std::vector<long long> > out;
    for (std::size_t k = 0; k < j.size(); ++k) out.push_back(ll_vector_from_json(j[k], at(where, k)));
    return out;
}
 
/**
 * spn_conv(S, V, g), or spn_conv(INV, g) with S the struct spn_sinvariants
 * returns ({places, S, V, m0}). In the struct form MATLAB binds g to its
 * second parameter; by name it is always "g".
 */
struct OpSpnConv {
    template <class T>
    static Json run(Args& a) {
        const Json& Sj = a.get("S");
        const std::vector<std::vector<T> > g = vectors_from_json<T>(a.get("g"), "spn_conv: g");
        T G;
        if (Sj.is_object()) {
            if (a.has("V"))
                throw InputError("spn_conv: S is the spn_sinvariants struct, which carries V; do not give V too");
            Args s(Sj, "spn_conv S");
            spn::SpnInvariants inv;
            inv.S = ll_rows_from_json(s.get("S"), "spn_conv: S.S");
            inv.V = ll_vector_from_json(s.get("V"), "spn_conv: S.V");
            if (const Json* p = s.opt("places")) {
                const std::vector<int> pl = int_vector_from_json(*p, "spn_conv: S.places");
                inv.places.assign(pl.begin(), pl.end());
            }
            if (const Json* m0 = s.opt("m0")) inv.m0 = ll_vector_from_json(*m0, "spn_conv: S.m0");
            s.done();
            a.done();
            G = spn::spn_conv(inv, g);
        } else {
            const std::vector<std::vector<long long> > S = ll_rows_from_json(Sj, "spn_conv: S");
            const std::vector<long long> V = ll_vector_from_json(a.get("V"), "spn_conv: V");
            a.done();
            G = spn::spn_conv(S, V, g);
        }
        Json r;
        r["G"] = encode_scalar(G);
        return r;
    }
};
 
/**
 * spn_ispn(pre, post, inh, rate, imm, srv, m0, p0, options): p0 is 1-based as
 * in MATLAB. inh, imm and srv may be empty (no inhibitors, all timed, one
 * server). options carries extra, maxlevel, maxstates, tol, distmax and
 * qbdmethod; rewards are callables and have no JSON form.
 */
struct OpSpnIspn {
    static std::vector<std::vector<double> > rows(const Matrix<double>& m) {
        std::vector<std::vector<double> > r(m.rows(), std::vector<double>(m.cols()));
        for (std::size_t i = 0; i < m.rows(); ++i)
            for (std::size_t j = 0; j < m.cols(); ++j) r[i][j] = m(i, j);
        return r;
    }
    static Json run(Args& a) {
        const std::vector<std::vector<double> > pre = rows(a.matrix<double>("pre"));
        const std::vector<std::vector<double> > post = rows(a.matrix<double>("post"));
        const std::vector<std::vector<double> > inh = rows(a.matrix_or_empty<double>("inh"));
        const std::vector<double> rate = a.vector<double>("rate");
        const std::vector<double> immd = a.vector_or_empty<double>("imm");
        const std::vector<double> srv = a.vector_or_empty<double>("srv");
        const std::vector<double> m0 = a.vector<double>("m0");
        const std::size_t P0 = a.count("p0");
        spn::IspnOptions o;
        if (const Json* oj = a.opt("options")) {
            Args op(*oj, "spn_ispn options");
            if (op.has("extra")) o.extra = static_cast<int>(op.count("extra"));
            if (op.has("maxlevel")) o.maxlevel = static_cast<int>(op.count("maxlevel"));
            if (op.has("maxstates")) o.maxstates = op.count("maxstates");
            if (op.has("tol")) o.tol = op.number<double>("tol");
            if (op.has("distmax")) o.distmax = op.number<double>("distmax");
            o.qbdmethod = op.text("qbdmethod", o.qbdmethod);
            op.unsupported("rewards", "they are functions of the marking, callables with no JSON form.");
            op.done();
        }
        a.done();
        if (P0 < 1) throw InputError(ErrorKind{"Argument"}, "spn_ispn: p0 is a 1-based place index");
        std::vector<bool> imm(immd.size());
        for (std::size_t t = 0; t < immd.size(); ++t) imm[t] = immd[t] != 0;
        const spn::IspnResult r = spn::spn_ispn(pre, post, inh, rate, imm, srv, m0, P0 - 1, o);
        Json res;
        res["K"] = encode_count(static_cast<std::size_t>(r.K));
        res["L"] = encode_count(r.L);
        res["A0"] = encode_vectors(r.A0);
        res["A1"] = encode_vectors(r.A1);
        res["A2"] = encode_vectors(r.A2);
        res["R"] = encode_vectors(r.R);
        res["B"] = encode_vectors(r.B);
        res["x"] = encode_vectors(r.x);
        Json lm = Json::array();
        for (const auto& m : r.levelmarkings) lm.push_back(encode_vectors(m));
        res["levelmarkings"] = lm;
        res["drift"] = encode_vector(r.drift);
        res["spR"] = encode_scalar(r.spR);
        res["meanmarking"] = encode_vector(r.meanmarking);
        res["varlevel"] = encode_scalar(r.varlevel);
        res["leveldist"] = encode_vector(r.leveldist);
        res["leveltail"] = encode_scalar(r.leveltail);
        res["marginal"] = encode_vectors(r.marginal);
        res["tput"] = encode_vector(r.tput);
        res["ntangible"] = encode_count(r.ntangible);
        res["nvanishing"] = encode_count(r.nvanishing);
        Json out;
        out["res"] = res;
        return out;
    }
};
 
/** spn_ispn_covered(C) with C = POST - PRE, P x T. */
struct OpSpnIspnCovered {
    static Json run(Args& a) {
        const Matrix<double> C = a.matrix<double>("C");
        a.done();
        const std::vector<bool> cov = spn::spn_ispn_covered(OpSpnIspn::rows(C));
        Json v = Json::array();
        for (bool b : cov) v.push_back(b);
        Json out;
        out["covered"] = v;
        return out;
    }
};
 
/**
 * spn_rec_enabled(mdds, g, mde, nplacelevels). mde is one entry of
 * INFO.modes from spn_mdd: enab and inhib (an entry may be "Inf") are what the
 * recursion reads; trans, mode, fire, D0, D1, pie, nph and srv are carried
 * into the port's SpnMode as given. dep is a callable.
 */
struct OpSpnRecEnabled {
    template <class T>
    static Json run(Args& a) {
        const mdd::MddStruct m = mdd_struct_from_json(a.get("mdds"), "spn_rec_enabled: mdds");
        const std::vector<std::vector<T> > g = vectors_from_json<T>(a.get("g"), "spn_rec_enabled: g");
        const Json& mj = a.get("mde");
        const std::size_t L = a.count("nplacelevels");
        a.done();
        if (!mj.is_object()) throw InputError("spn_rec_enabled: mde is one entry of INFO.modes, an object");
        Args md(mj, "spn_rec_enabled mde");
        spn::SpnMode<T> mode;
        mode.enab = md.vector<double>("enab");
        mode.inhib = md.vector<double>("inhib");
        if (md.has("fire")) mode.fire = md.vector<double>("fire");
        if (md.has("trans")) mode.trans = md.count("trans");
        if (md.has("mode")) mode.mode = md.count("mode");
        if (md.has("D0")) mode.D0 = md.matrix<T>("D0");
        if (md.has("D1")) mode.D1 = md.matrix<T>("D1");
        if (md.has("pie")) mode.pie = md.vector<T>("pie");
        if (md.has("nph")) mode.nph = md.count("nph");
        if (md.has("srv")) mode.srv = md.number<double>("srv");
        md.unsupported("dep", "it is a marking-dependent rate function, a callable with no JSON form.");
        md.done();
        if (mode.enab.size() < L || mode.inhib.size() < L)
            throw InputError("spn_rec_enabled: mde.enab and mde.inhib need one entry per place level");
        const spn::SpnEnabling<T> en = spn::spn_rec_enabled(m, g, mode, L);
        Json r;
        Json ej;
        ej["ge"] = encode_vector(en.ge);
        ej["eq"] = encode_vector(en.eq);
        ej["maxDegree"] = encode_count(en.max_degree);
        r["en"] = ej;
        return r;
    }
};
 
 
}  // namespace
 
void register_api_misc(DispatchTable& t) {
    t["polling_qsys_1limited"] = make<OpPolling1Limited>();
    t["polling_qsys_decrementing"] = make<OpPollingDecrementing>();
    t["polling_pcl"] = make<OpPollingPcl>();
    t["polling_qsys_exhaustive"] = make<OpPollingExhaustive>();
    t["polling_qsys_gated"] = make<OpPollingGated>();
    t["polling_qsys_klimited"] = make<OpPollingKLimited>();
    t["polling_qsys_table"] = make<OpPollingTable>();
    t["ljd_linearize"] = make<OpLjdLinearize>();
    t["fes_compute_throughputs"] = make<OpFesComputeThroughputs>();
    t["fes_compute_metrics"] = make<OpFesComputeMetrics>();
    t["fes_build_isolated"] = make<OpFesBuildIsolated>();
    t["fes_map_levels"] = make<OpFesMapLevels>();
    t["fes_map_interdeparture"] = make<OpFesMapInterdeparture>();
    t["fes_map_interp"] = make<OpFesMapInterp>();
    t["fes_map_moments"] = make<OpFesMapMoments>();
    t["fes_map_solve"] = make<OpFesMapSolve>();
    t["fes_map_aggregate"] = make<NeedsTranscendental<OpFesMapAggregate> >();
    t["cache_erec"] = make<OpCacheErec>();
    t["cache_prob_erec"] = make<OpCacheProbErec>();
    t["cache_cost"] = make<OpCacheCost>();
    t["cache_cost_pathcheck"] = make<OpCacheCostPathcheck>();
    t["cache_gamma"] = make<OpCacheGamma>();
    t["cache_gamma_lp"] = make<OpCacheGammaLp>();
    t["cache_miss"] = make<OpCacheMiss>();
    t["cache_mva"] = make<OpCacheMva>();
    t["cache_mva_miss"] = make<OpCacheMvaMiss>();
    t["cache_rrm_meanfield_ode"] = make<OpCacheRrmMeanfieldOde>();
    t["cache_spm_size"] = make<NeedsTranscendental<OpCacheSpmSize> >();
    t["cache_miss_asy"] = make<NeedsTranscendental<OpCacheMissAsy> >();
    t["cache_miss_fpi"] = make<NeedsTranscendental<OpCacheMissFpi> >();
    t["cache_miss_spm"] = make<NeedsTranscendental<OpCacheMissSpm> >();
    t["cache_miss_rmf"] = make<NeedsTranscendental<OpCacheMissRmf> >();
    t["cache_rrm_meanfield"] = make<NeedsTranscendental<OpCacheRrmMeanfield> >();
    t["cache_prob_fpi"] = make<NeedsTranscendental<OpCacheProbFpi> >();
    t["cache_prob_spm"] = make<NeedsTranscendental<OpCacheProbSpm> >();
    t["cache_ttl_hlru"] = make<NeedsTranscendental<OpCacheTtlHlru> >();
    t["cache_t_hlru"] = make<NeedsTranscendental<OpCacheTHlru> >();
    t["cache_spm"] = make<NeedsTranscendental<OpCacheSpm> >();
    t["cache_xi_fp"] = make<NeedsTranscendental<OpCacheXiFp> >();
    t["cache_xi_iter"] = make<NeedsTranscendental<OpCacheXiIter> >();
    t["cache_is"] = make<NeedsTranscendental<OpCacheIs> >();
    t["cache_prob_is"] = make<NeedsTranscendental<OpCacheProbIs> >();
    t["cache_miss_is"] = make<NeedsTranscendental<OpCacheMissIs> >();
    t["cache_lrum_map_levelstats"] = make<NeedsTranscendental<OpCacheLrumMapLevelstats> >();
    t["cache_t_lrum_map"] = make<NeedsTranscendental<OpCacheTLrumMap> >();
    t["cache_ttl_lrum_map"] = make<NeedsTranscendental<OpCacheTtlLrumMap> >();
    t["cache_ttl_lrua"] = make<NeedsTranscendental<OpCacheTtlLrua> >();
    t["da_cache_isolate"] = make<OpDaCacheIsolate>();
    t["da_traffic_superpos"] = make<OpDaTrafficSuperpos>();
    t["lossn_manjunath"] = make<OpLossnManjunath>();
    t["lossn_rec"] = make<OpLossnRec>();
    t["lossn_kr"] = make<OpLossnKr>();
    t["lossn_kr_dim"] = make<OpLossnKrDim>();
    t["dqsys_bernoulli1"] = make<OpDqsysBernoulli1>();
    t["map2ph"] = make<OpMap2ph>();
    t["map2mmpp"] = make<OpMap2mmpp>();
    t["autocov"] = make<OpAutocov>();
    t["ph_multisets"] = make<OpPhMultisets>();
    t["mmdp_isfeasible"] = make<OpMmdpIsfeasible>();
    t["pas_placement"] = make<OpPasPlacement>();
    t["lossn_erlangfp"] = make<NeedsTranscendental<OpLossnErlangfp> >();
    t["lossn_mci"] = make<NeedsTranscendental<OpLossnMci> >();
    t["lossn_kr_asy"] = make<NeedsTranscendental<OpLossnKrAsy> >();
    t["lossn_asy"] = make<NeedsTranscendental<OpLossnAsy> >();
    t["lossn_links"] = make<NeedsTranscendental<OpLossnLinks> >();
    t["lossn_revsens"] = make<NeedsTranscendental<OpLossnRevsens> >();
    t["aph_fit"] = make<NeedsTranscendental<OpAphFit> >();
    t["m3pp22_fitc_approx_cov"] = make<NeedsTranscendental<OpM3pp22FitcApproxCov> >();
    t["infradius_h"] = make<NeedsTranscendental<OpInfradiusH> >();
    t["mapqn_amva"] = make<OpMapqnAmva>();
    t["mapqn_bnd_lr"] = make<OpMapqnBndLr>();
    t["mapqn_bnd_lr_mva"] = make<OpMapqnBndLrMva>();
    t["mapqn_bnd_qr"] = make<OpMapqnBndQr>();
    t["mapqn_bnd_qr_ld"] = make<OpMapqnBndQrLd>();
    t["mapqn_bnd_qr_delay"] = make<OpMapqnBndQrDelay>();
    t["mapqn_qr_bounds_bas"] = make<OpMapqnQrBoundsBas>();
    t["mapqn_qr_bounds_rsrd"] = make<OpMapqnQrBoundsRsrd>();
    t["retrieval_metrics"] = make<OpRetrievalMetrics>();
    t["retrieval_mva"] = make<OpRetrievalMva>();
    t["retrieval_nc"] = make<OpRetrievalNc>();
    t["wf_loop_detector"] = make<OpWfLoopDetector>();
    t["wf_parallel_detector"] = make<OpWfParallelDetector>();
    t["wf_sequence_detector"] = make<OpWfSequenceDetector>();
    t["mdd_rec"] = make<OpMddRec>();
    t["mdd_closedqn"] = make<OpMddClosedqn>();
    t["spn_conv"] = make<OpSpnConv>();
    t["spn_ispn"] = make<DoubleOnly<OpSpnIspn> >();
    t["spn_ispn_covered"] = make<DoubleOnly<OpSpnIspnCovered> >();
    t["spn_rec_enabled"] = make<OpSpnRecEnabled>();
    t["retrieval_fpi"] = make<NeedsTranscendental<OpRetrievalFpi> >();
    t["retrieval_fpi_latency"] = make<NeedsTranscendental<OpRetrievalFpiLatency> >();
    t["retrieval_rayint"] = make<NeedsTranscendental<OpRetrievalRayint> >();
    t["sum_closed"] = make<NeedsTranscendental<OpSumClosed> >();
    t["sum_closing"] = make<NeedsTranscendental<OpSumClosing> >();
    t["wf_branch_detector"] = make<NeedsTranscendental<OpWfBranchDetector> >();
}
 
void register_library_only_misc(LibraryOnlyTable& t) {
    const std::string fn_arg =
        "it takes a C++ callable as an argument, which has no JSON form.";
    t["cd_peak_scaling"] = fn_arg;
    t["da_fpi"] = fn_arg;
    t["hyperexp_fit_longtail"] = fn_arg;
    t["laplace_invert"] = fn_arg;
    t["laplaceapprox"] = fn_arg;
    t["pas_swap2order"] = fn_arg;
    t["mdd_reachset"] = fn_arg;
    t["fes_beta_handle"] =
        "it returns a C++ callable (the per-class rate function of the aggregate), which has "
        "no JSON form.";
    const std::string mdd_desc =
        "it returns a Kronecker rate descriptor whose successor function is a C++ callable, "
        "which has no JSON form; run line-cli -s ctmc -a <analysis> on the model instead.";
    t["mdd_descriptor"] = mdd_desc;
    t["mdd_ps"] = mdd_desc;
    t["mdd_mcd"] =
        "it takes a Kronecker rate descriptor whose successor function is a C++ callable, "
        "which has no JSON form; run line-cli -s ctmc -a <analysis> on the model instead.";
    const std::string sn_mva =
        "it takes a NetworkStruct; build the model and run line-cli -s mva -a <analysis> "
        "on it instead.";
    t["cache_retrieval_inputs"] = sn_mva;
    t["fes_validate"] = sn_mva;
    t["da_cacheqn"] = sn_mva;
    t["da_cacheqn_retrieval"] = sn_mva;
    t["spn_lpbnd"] =
        "it takes a NetworkStruct; build the model and run line-cli -s ba -a <analysis> "
        "on it instead.";
    t["spn_mdd"] =
        "it takes a NetworkStruct; build the model and run line-cli -s ctmc -a <analysis> "
        "on it instead.";
    const std::string sn_nc =
        "it takes a NetworkStruct; build the model and run line-cli -s nc -a <analysis> "
        "on it instead.";
    t["spn_pf"] = sn_nc;
    t["spn_sinvariants"] = sn_nc;
    const std::string sn_props =
        "it takes a NetworkStruct and returns the Petri-net property structs; build the "
        "model and call it from C++ (line/api/spn/spn_properties.h) instead.";
    // the shared net model the property functions build; same carrier, same reason
    t["spn_net"] = sn_props;
    t["spn_properties"] = sn_props;
    t["spn_reachgraph"] = sn_props;
    t["spn_reachprops"] = sn_props;
    t["spn_siphons"] = sn_props;
    t["spn_structprops"] = sn_props;
    t["spn_tinvariants"] = sn_props;
    t["spn_traps"] = sn_props;
    t["spn_metrics"] =
        "it takes the SpnInfo that spn_mdd returns, which holds an MDD object and callable "
        "marking-dependence laws with no JSON form; build the net and run line-cli -s nc "
        "-a <analysis> on it instead.";
}
 
}  // namespace reg
}  // namespace line