api_dispatch_qsys.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 'qsys' group (see api_dispatch_detail.h):
 * the remaining closed-form queueing systems, the fork-join (fj) family, the
 * age-of-information transforms, the moment conversions and the simulation
 * planning helpers.
 *
 * CONVENTIONS, shared with api_dispatch.cpp:
 *  - argument keys are the MATLAB parameter names, output keys the MATLAB
 *    output names (struct fields where MATLAB returns a struct);
 *  - an absent optional argument calls the port's shorter overload, or leaves
 *    the port's own default in place; no default is restated here. Where the
 *    port takes its optional arguments POSITIONALLY, a later one cannot be
 *    given without the earlier ones (prefix_count), and where it takes a group
 *    of tuning knobs only together, the group is all-or-none (all_or_none);
 *  - a function handle cannot cross JSON. Where the port offers a constructor
 *    for a NAMED law, the handle is replaced by that law (a model.json
 *    distribution record, or the aoi_lst_* constructors); otherwise the
 *    argument is refused by name, and a function whose handle is REQUIRED is
 *    listed as library-only below.
 */
 
#include "api_dispatch_detail.h"
 
#include <cmath>
#include <complex>
#include <initializer_list>
#include <string>
#include <vector>
 
#include "line/api/aoi/aoi_dist2ph.h"
#include "line/api/aoi/aoi_lst_det.h"
#include "line/api/aoi/aoi_lst_erlang.h"
#include "line/api/aoi/aoi_lst_exp.h"
#include "line/api/aoi/aoi_lst_ph.h"
#include "line/api/fj/fj_amva.h"
#include "line/api/fj/fj_branch_members.h"
#include "line/api/fj/fj_char_max_blom.h"
#include "line/api/fj/fj_char_max_discrete.h"
#include "line/api/fj/fj_cox_fit.h"
#include "line/api/fj/fj_dag_makespan.h"
#include "line/api/fj/fj_delay_opt.h"
#include "line/api/fj/fj_dispersion.h"
#include "line/api/fj/fj_dist2fj.h"
#include "line/api/fj/fj_ism_green.h"
#include "line/api/fj/fj_lst_max_het.h"
#include "line/api/fj/fj_mg1_respt_moments.h"
#include "line/api/fj/fj_qgb.h"
#include "line/api/fj/fj_respt_bulk.h"
#include "line/api/fj/fj_respt_closed.h"
#include "line/api/fj/fj_respt_nosplit.h"
#include "line/api/fj/fj_serialization.h"
#include "line/api/fj/fj_tail_forktail.h"
#include "line/api/fj/fj_tail_ordstat.h"
#include "line/api/fj/fj_tsm_capacity.h"
#include "line/api/fj/fj_xmax_coxian.h"
#include "line/api/fj/fj_xmax_het.h"
#include "line/api/fj/fj_xmax_hz.h"
#include "line/api/fj/fj_xmax_moments_het.h"
#include "line/api/moment/moment_cumulant.h"
#include "line/api/moment/moment_housematrix.h"
#include "line/api/moment/moment_joint.h"
#include "line/api/moment/moment_tail.h"
#include "line/api/moment/moment_tensor.h"
// qsys_bmapm1.h names lang::GlobalConstants without including the header that
// declares it, so lang_types.h must come first.
#include "line/lang/lang_types.h"
#include "line/api/qsys/qsys_bmapm1.h"
#include "line/api/qsys/qsys_hh1_lindley.h"
#include "line/api/qsys/qsys_ldps_workload.h"
#include "line/api/qsys/qsys_lindley_moment.h"
#include "line/api/qsys/qsys_mapg1k.h"
#include "line/api/qsys/qsys_mapg1k_perflow.h"
#include "line/api/qsys/qsys_mapmapc.h"
#include "line/api/qsys/qsys_maxima_twomoment.h"
#include "line/api/qsys/qsys_mg1_ps.h"
#include "line/api/qsys/qsys_ddckk.h"
#include "line/api/qsys/qsys_engset.h"
#include "line/api/qsys/qsys_gimc.h"
#include "line/api/qsys/qsys_md1.h"
#include "line/api/qsys/qsys_mer1.h"
#include "line/api/qsys/qsys_mg1_busy.h"
#include "line/api/qsys/qsys_mg1_conservation.h"
#include "line/api/qsys/qsys_mg1_exceptional.h"
#include "line/api/qsys/qsys_mg1_prio_pr.h"
#include "line/api/qsys/qsys_mg1_setup.h"
#include "line/api/qsys/qsys_mg1_sjf.h"
#include "line/api/qsys/qsys_mg1kk.h"
#include "line/api/qsys/qsys_mgc_ps.h"
#include "line/api/qsys/qsys_mh21.h"
#include "line/api/qsys/qsys_mmc_prio.h"
#include "line/api/qsys/qsys_mmc_prio_eta.h"
#include "line/api/qsys/qsys_mmc_prio_tardiness.h"
#include "line/api/qsys/qsys_mmcc.h"
#include "line/api/qsys/qsys_mmckk.h"
#include "line/api/qsys/qsys_mgisrgi_whitt.h"
#include "line/api/qsys/qsys_mm1_lindley.h"
#include "line/api/qsys/qsys_mm1_ps.h"
#include "line/api/qsys/qsys_mm1_tandem_lindley.h"
#include "line/api/qsys/qsys_mmapg1k.h"
#include "line/api/qsys/qsys_mmapgk1.h"
#include "line/api/qsys/qsys_mmapgk1_lcfspr.h"
#include "line/api/qsys/qsys_mmapmap1.h"
#include "line/api/qsys/qsys_tandem_lindley.h"
#include "line/api/qsys/qsys_tandem_ub_ciucu.h"
#include "line/api/sim/sim_quest_heuristic_ci.h"
#include "line/api/sim/sim_quest_options.h"
#include "line/api/sim/sim_runlength.h"
#include "line/io/network_reader.h"
#include "line/lang/distribution.h"
 
namespace line {
namespace reg {
 
namespace {
 
// ---------------------------------------------------------------------------
// Shared readers
// ---------------------------------------------------------------------------
 
/**
 * How many of the port's POSITIONAL optional arguments the caller gave. They
 * must form a prefix: naming a later one without the earlier would make this
 * boundary restate the earlier defaults, which it does not do.
 */
std::size_t prefix_count(const Args& a, std::initializer_list<const char*> keys,
                         const std::string& fn) {
    std::size_t n = 0;
    bool gap = false;
    const char* missing = nullptr;
    for (const char* k : keys) {
        if (a.has(k)) {
            if (gap)
                throw InputError(fn + ": '" + k + "' follows '" + missing +
                                 "' positionally in the port; give '" + missing + "' as well");
            ++n;
        } else if (!gap) {
            gap = true;
            missing = k;
        }
    }
    return n;
}
 
/** A group of tuning knobs the port takes only together: all of them or none. */
bool all_or_none(const Args& a, std::initializer_list<const char*> keys, const std::string& fn) {
    std::size_t n = 0;
    std::string names;
    for (const char* k : keys) {
        if (a.has(k)) ++n;
        if (!names.empty()) names += ", ";
        names += std::string("'") + k + "'";
    }
    if (n != 0 && n != keys.size())
        throw InputError(fn + ": " + names + " tune one computation together; give all of them "
                              "or none");
    return n != 0;
}
 
std::vector<unsigned> read_uints(Args& a, const char* key, const std::string& fn) {
    const std::vector<int> v = a.ints(key);
    std::vector<unsigned> out;
    for (std::size_t i = 0; i < v.size(); ++i) {
        if (v[i] < 0)
            throw InputError(fn + ": '" + key + "' holds counts and cannot be negative");
        out.push_back(static_cast<unsigned>(v[i]));
    }
    return out;
}
 
/** A required flag, which Args::boolean only offers with a fallback. */
bool required_flag(Args& a, const char* key, const std::string& fn) {
    if (!a.has(key)) throw InputError(fn + ": missing required argument '" + key + "'");
    return a.boolean(key, false);
}
 
/** A MAP as the (D0, D1) pair the reference passes around, under any two keys. */
template <class T>
mam::Map<T> read_map(Args& a, const char* d0, const char* d1) {
    mam::Map<T> m;
    m.D0 = a.matrix<T>(d0);
    m.D1 = a.matrix<T>(d1);
    return m;
}
 
/** A MAP given as MATLAB's cell {D0, D1}: a two-element list of matrices. */
template <class T>
mam::Map<T> map_from_pair(const Json& j, const std::string& where) {
    if (!j.is_array() || j.size() != 2)
        throw InputError(where + ": expected the MAP as a list [D0, D1] of two matrices");
    mam::Map<T> m;
    m.D0 = matrix_from_json<T>(j[0], where + "[0]");
    m.D1 = matrix_from_json<T>(j[1], where + "[1]");
    return m;
}
 
/** A list of MAPs, each a [D0, D1] pair. */
template <class T>
std::vector<mam::Map<T> > maps_from_json(const Json& j, const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected a list of [D0, D1] pairs");
    std::vector<mam::Map<T> > out;
    for (std::size_t i = 0; i < j.size(); ++i)
        out.push_back(map_from_pair<T>(j[i], where + "[" + std::to_string(i) + "]"));
    return out;
}
 
/**
 * A distribution given as its model.json record, e.g. {"type":"Exp",
 * "params":{"lambda":2}}. The record reader parses its numbers through double,
 * so at --arith real a law is as precise as a double literal, and a malformed
 * record is reported as the caller's input rather than as a JSON library fault.
 */
template <class T>
lang::Distrib<T> dist_record(const Json& j, const std::string& where) {
    if (!j.is_object() || !j.contains("type"))
        throw InputError(where + ": expected a distribution record {\"type\": ..., ...} as "
                                 "written to model.json");
    try {
        return io::detail::dist_from_json<T>(j);
    } catch (const Json::exception& e) {
        throw InputError(where + ": malformed distribution record: " + e.what());
    }
}
 
/**
 * The service law of the MAP/G/1/K family, MATLAB's svc struct: {"type":"det",
 * "d"}, {"type":"gamma","alpha","theta"} or {"type":"ph","alpha","T"}. The
 * 'density' type carries a function handle and has no JSON form.
 */
template <class T>
qsys::ServiceLaw<T> read_svc(Args& a, const std::string& fn) {
    Args sv(a.get("svc"), fn + ": svc");
    const std::string type = sv.text("type", "");
    qsys::ServiceLaw<T> svc;
    if (type == "det") {
        svc = qsys::ServiceLaw<T>::deterministic(sv.number<T>("d"));
    } else if (type == "gamma") {
        const T al = sv.number<T>("alpha");
        svc = qsys::ServiceLaw<T>::gamma(al, sv.number<T>("theta"));
    } else if (type == "ph") {
        const std::vector<T> al = sv.vector<T>("alpha");
        svc = qsys::ServiceLaw<T>::phase_type(al, sv.matrix<T>("T"));
    } else if (type == "density") {
        throw UnsupportedError(ErrorKind{"Api"},
            fn + ": a 'density' service law is a function handle and has no named-JSON form; "
                 "use 'det', 'gamma' or 'ph', or call it from the library");
    } else {
        throw UnsupportedError(ErrorKind{"Distribution"},
            fn + ": svc.type must be 'det', 'gamma' or 'ph', got '" + type + "'");
    }
    sv.done();
    return svc;
}
 
/** A vector output that MATLAB reshapes to the scalar it was given. */
template <class T>
Json encode_shaped(const std::vector<T>& v, bool scalar) {
    if (scalar && v.size() == 1) return encode_scalar(v[0]);
    return encode_vector(v);
}
 
/** The MATLAB message every exact refusal of the sim family shares. */
std::string sim_exact_refusal(const char* fn, const char* why) {
    return std::string(fn) + ": " + why +
           ", so exact rational arithmetic has nothing to preserve and is refused; rerun with "
           "--arith double or --arith real";
}
 
// ---------------------------------------------------------------------------
// Moment tensors: nested JSON arrays, outer index = dimension 1
// ---------------------------------------------------------------------------
 
template <class T>
void tensor_fill(const Json& node, std::size_t depth, std::size_t offset,
                 const std::vector<std::size_t>& sz, const std::vector<std::size_t>& stride,
                 std::vector<T>& data, const std::string& where) {
    if (depth == sz.size()) {
        if (node.is_array())
            throw InputError(where + ": the tensor is nested deeper than its first element");
        data[offset] = number_from_json<T>(node, where);
        return;
    }
    if (!node.is_array() || node.size() != sz[depth])
        throw InputError(where + ": ragged tensor, dimension " + std::to_string(depth + 1) +
                         " must have " + std::to_string(sz[depth]) + " entries throughout");
    for (std::size_t i = 0; i < sz[depth]; ++i)
        tensor_fill<T>(node[i], depth + 1, offset + i * stride[depth], sz, stride, data, where);
}
 
/**
 * A moment tensor from nested arrays: the nesting depth is the order and the
 * outermost index runs over dimension 1, so A[i][j] is MATLAB's A(i,j). A flat
 * array is an order-1 tensor, a bare number a one-entry one.
 */
template <class T>
moment::MomentTensor<T> tensor_from_json(const Json& j, const std::string& where) {
    std::vector<std::size_t> sz;
    if (!j.is_array()) {
        sz.push_back(1);
    } else {
        const Json* p = &j;
        while (p->is_array()) {
            sz.push_back(p->size());
            if (p->empty()) break;
            p = &(*p)[0];
        }
    }
    std::vector<std::size_t> stride(sz.size(), 1);
    for (std::size_t l = 1; l < sz.size(); ++l) stride[l] = stride[l - 1] * sz[l - 1];
    moment::MomentTensor<T> t(sz);
    if (!j.is_array())
        t.data[0] = number_from_json<T>(j, where);
    else if (t.numel() > 0)
        tensor_fill<T>(j, 0, 0, sz, stride, t.data, where);
    return t;
}
 
template <class T>
Json tensor_nest(const moment::MomentTensor<T>& t, std::size_t depth, std::size_t offset) {
    Json out = Json::array();
    const std::size_t st = t.stride(depth);
    for (std::size_t i = 0; i < t.sz[depth]; ++i) {
        if (depth + 1 == t.order())
            out.push_back(encode_scalar(t.data[offset + i * st]));
        else
            out.push_back(tensor_nest(t, depth + 1, offset + i * st));
    }
    return out;
}
 
/** The inverse of tensor_from_json; an order-1 tensor is a flat array. */
template <class T>
Json encode_tensor(const moment::MomentTensor<T>& t) {
    return tensor_nest(t, 0, 0);
}
 
Json encode_sizes(const std::vector<std::size_t>& v) {
    Json out = Json::array();
    for (std::size_t i = 0; i < v.size(); ++i) out.push_back(encode_count(v[i]));
    return out;
}
 
// ---------------------------------------------------------------------------
// moment: one-dimensional sequence conversions, vector in and vector out
// ---------------------------------------------------------------------------
 
#define LINE_MOMENT_VEC_OP(OPNAME, FN, IN, OUT)                                      \
    struct OPNAME {                                                                  \
        template <class T>                                                           \
        static Json run(Args& a) {                                                   \
            const std::vector<T> x = a.vector<T>(IN);                                \
            a.done();                                                                \
            Json out;                                                                \
            out[OUT] = encode_vector(moment::FN<T>(x));                              \
            return out;                                                              \
        }                                                                            \
    };
 
LINE_MOMENT_VEC_OP(OpMomentCumulantFromRaw, moment_cumulant_from_raw, "m", "kappa")
LINE_MOMENT_VEC_OP(OpMomentRawFromCumulant, moment_raw_from_cumulant, "kappa", "m")
LINE_MOMENT_VEC_OP(OpMomentFactcumulantFromFactorial, moment_factcumulant_from_factorial, "f",
                   "kappa")
LINE_MOMENT_VEC_OP(OpMomentFactorialFromFactcumulant, moment_factorial_from_factcumulant,
                   "kappa", "f")
LINE_MOMENT_VEC_OP(OpMomentBinomialFromTail, moment_binomial_from_tail, "t", "b")
LINE_MOMENT_VEC_OP(OpMomentTailFromBinomial, moment_tail_from_binomial, "b", "t")
#undef LINE_MOMENT_VEC_OP
 
// ---------------------------------------------------------------------------
// moment: joint (tensor) conversions, tensor in and tensor out
// ---------------------------------------------------------------------------
 
#define LINE_MOMENT_TENSOR_OP(OPNAME, FN, IN, OUT)                                   \
    struct OPNAME {                                                                  \
        template <class T>                                                           \
        static Json run(Args& a) {                                                   \
            const moment::MomentTensor<T> x =                                        \
                tensor_from_json<T>(a.get(IN), std::string(#FN ": ") + IN);          \
            a.done();                                                                \
            Json out;                                                                \
            out[OUT] = encode_tensor(moment::FN<T>(x));                              \
            return out;                                                              \
        }                                                                            \
    };
 
LINE_MOMENT_TENSOR_OP(OpMjBinomialFromFactorial, moment_joint_binomial_from_factorial, "f", "b")
LINE_MOMENT_TENSOR_OP(OpMjBinomialFromNegbinomial, moment_joint_binomial_from_negbinomial, "bm",
                      "b")
LINE_MOMENT_TENSOR_OP(OpMjBinomialFromTail, moment_joint_binomial_from_tail, "t", "b")
LINE_MOMENT_TENSOR_OP(OpMjCentralFromRaw, moment_joint_central_from_raw, "m", "mc")
LINE_MOMENT_TENSOR_OP(OpMjCentralFromTail, moment_joint_central_from_tail, "t", "mc")
LINE_MOMENT_TENSOR_OP(OpMjCumulantFromRaw, moment_joint_cumulant_from_raw, "m", "kappa")
LINE_MOMENT_TENSOR_OP(OpMjFactcumulantFromFactorial, moment_joint_factcumulant_from_factorial,
                      "f", "kappa")
LINE_MOMENT_TENSOR_OP(OpMjFactorialFromBinomial, moment_joint_factorial_from_binomial, "b", "f")
LINE_MOMENT_TENSOR_OP(OpMjFactorialFromFactcumulant, moment_joint_factorial_from_factcumulant,
                      "kappa", "f")
LINE_MOMENT_TENSOR_OP(OpMjFactorialFromRaw, moment_joint_factorial_from_raw, "m", "f")
LINE_MOMENT_TENSOR_OP(OpMjFactorialFromUpfactorial, moment_joint_factorial_from_upfactorial, "fp",
                      "f")
LINE_MOMENT_TENSOR_OP(OpMjNegbinomialFromBinomial, moment_joint_negbinomial_from_binomial, "b",
                      "bm")
LINE_MOMENT_TENSOR_OP(OpMjNegbinomialFromUpfactorial, moment_joint_negbinomial_from_upfactorial,
                      "fp", "bm")
LINE_MOMENT_TENSOR_OP(OpMjRawFromCumulant, moment_joint_raw_from_cumulant, "kappa", "m")
LINE_MOMENT_TENSOR_OP(OpMjRawFromFactorial, moment_joint_raw_from_factorial, "f", "m")
LINE_MOMENT_TENSOR_OP(OpMjRawFromUpfactorial, moment_joint_raw_from_upfactorial, "fp", "m")
LINE_MOMENT_TENSOR_OP(OpMjTailFromBinomial, moment_joint_tail_from_binomial, "b", "t")
LINE_MOMENT_TENSOR_OP(OpMjUpfactorialFromFactorial, moment_joint_upfactorial_from_factorial, "f",
                      "fp")
LINE_MOMENT_TENSOR_OP(OpMjUpfactorialFromNegbinomial, moment_joint_upfactorial_from_negbinomial,
                      "bm", "fp")
LINE_MOMENT_TENSOR_OP(OpMjUpfactorialFromRaw, moment_joint_upfactorial_from_raw, "m", "fp")
#undef LINE_MOMENT_TENSOR_OP
 
struct OpMjCentralFromRawMean {
    template <class T>
    static Json run(Args& a) {
        const moment::MomentTensor<T> m =
            tensor_from_json<T>(a.get("m"), "moment_joint_central_from_raw_mean: m");
        const std::vector<T> mu = a.vector<T>("mu");
        a.done();
        Json out;
        out["mc"] = encode_tensor(moment::moment_joint_central_from_raw_mean<T>(m, mu));
        return out;
    }
};
 
struct OpMjRawFromCentral {
    template <class T>
    static Json run(Args& a) {
        const moment::MomentTensor<T> mc =
            tensor_from_json<T>(a.get("mc"), "moment_joint_raw_from_central: mc");
        const std::vector<T> mu = a.vector<T>("mu");
        a.done();
        Json out;
        out["m"] = encode_tensor(moment::moment_joint_raw_from_central<T>(mc, mu));
        return out;
    }
};
 
struct OpMjAggregate {
    template <class T>
    static Json run(Args& a) {
        const moment::MomentTensor<T> F =
            tensor_from_json<T>(a.get("F"), "moment_joint_aggregate: F");
        a.done();
        Json out;
        out["f"] = encode_vector(moment::moment_joint_aggregate<T>(F));
        return out;
    }
};
 
struct OpMjMarking {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> f = a.vector<T>("f");
        const std::vector<T> p = a.vector<T>("p");
        const std::vector<unsigned> d = read_uints(a, "dims", "moment_joint_marking");
        a.done();
        const std::vector<std::size_t> dims(d.begin(), d.end());
        Json out;
        out["F"] = encode_tensor(moment::moment_joint_marking<T>(f, p, dims));
        return out;
    }
};
 
struct OpMomentHousematrix {
    template <class T>
    static Json run(Args& a) {
        const std::string edge = a.text("edge", "");
        if (edge.empty()) throw InputError("moment_housematrix: missing required argument 'edge'");
        const int n = a.required_integer("n");
        a.done();
        Json out;
        out["T"] = encode_matrix(moment::moment_housematrix<T>(edge, n));
        return out;
    }
};
 
struct OpMomentTensorsize {
    template <class T>
    static Json run(Args& a) {
        const moment::MomentTensor<T> A = tensor_from_json<T>(a.get("A"), "moment_tensorsize: A");
        a.done();
        Json out;
        out["sz"] = encode_sizes(moment::moment_tensorsize<T>(A));
        return out;
    }
};
 
struct OpMomentTensortrans {
    template <class T>
    static Json run(Args& a) {
        const moment::MomentTensor<T> A = tensor_from_json<T>(a.get("A"), "moment_tensortrans: A");
        const Matrix<T> Tm = a.matrix<T>("T");
        const std::size_t mode = a.count("mode");
        a.done();
        // MATLAB's mode is 1-based, the port's 0-based.
        if (mode < 1) throw InputError("moment_tensortrans: 'mode' is 1-based and must be >= 1");
        Json out;
        out["B"] = encode_tensor(moment::moment_tensortrans<T>(A, Tm, mode - 1));
        return out;
    }
};
 
struct OpMomentJointtrans {
    template <class T>
    static Json run(Args& a) {
        const moment::MomentTensor<T> A = tensor_from_json<T>(a.get("A"), "moment_jointtrans: A");
        const std::string edge = a.text("edge", "");
        if (edge.empty()) throw InputError("moment_jointtrans: missing required argument 'edge'");
        a.done();
        Json out;
        out["B"] = encode_tensor(moment::moment_jointtrans<T>(A, edge));
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// fj: fork-join response times and maxima
// ---------------------------------------------------------------------------
 
/**
 * fj_amva / fj_qgb: MATLAB defaults Z to 0, but the port has no overload
 * without Z, so Z is required here rather than restated.
 */
struct OpFjAmva {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> D = a.vector<T>("D");
        const std::vector<unsigned> P = read_uints(a, "P", "fj_amva");
        const unsigned M = a.uinteger("M");
        const T Z = a.number<T>("Z");
        a.done();
        const fj::FJAmvaResult<T> r = fj::fj_amva<T>(D, P, M, Z);
        Json out;
        out["R"] = encode_vector(r.R);
        out["Q"] = encode_vector(r.Q);
        out["X"] = encode_scalar(r.X);
        out["U"] = encode_vector(r.U);
        return out;
    }
};
 
struct OpFjQgb {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> D = a.vector<T>("D");
        const std::vector<unsigned> P = read_uints(a, "P", "fj_qgb");
        const unsigned M = a.uinteger("M");
        const T Z = a.number<T>("Z");
        a.done();
        const fj::FJQgbResult<T> r = fj::fj_qgb<T>(D, P, M, Z);
        Json out;
        out["Q"] = encode_vector(r.Q);
        out["y"] = encode_vector(r.y);
        return out;
    }
};
 
struct OpFjCharMaxBlom {
    template <class T>
    static Json run(Args& a) {
        const unsigned K = a.uinteger("K");
        a.unsupported("Finv", "it is a quantile function handle; omit it for the standard "
                              "normal the port selects");
        const std::size_t n = prefix_count(a, {"alpha", "beta"}, "fj_char_max_blom");
        const T alpha = n >= 1 ? a.number<T>("alpha") : T();
        const T beta = n >= 2 ? a.number<T>("beta") : T();
        a.done();
        const std::function<T(const T&)> none;
        const fj::FJCharMaxBlomResult<T> r =
            n == 0 ? fj::fj_char_max_blom<T>(K)
                   : (n == 1 ? fj::fj_char_max_blom<T>(K, none, alpha)
                             : fj::fj_char_max_blom<T>(K, none, alpha, beta));
        Json out;
        out["mK"] = encode_scalar(r.mK);
        // MATLAB returns NaN for the bracket when the port's guard says it does
        // not exist; the absence is reported by omitting it.
        if (r.bracket_available) {
            out["lo"] = encode_scalar(r.lo);
            out["hi"] = encode_scalar(r.hi);
        }
        return out;
    }
};
 
struct OpFjCharMaxDiscrete {
    template <class T>
    static Json run(Args& a) {
        const unsigned K = a.uinteger("K");
        std::string type = a.text("dist_type", "");
        const T par = a.number<T>("param");
        a.done();
        for (std::size_t i = 0; i < type.size(); ++i)
            type[i] = static_cast<char>(std::tolower(static_cast<unsigned char>(type[i])));
        fj::FJDiscreteDist dist;
        if (type == "geometric" || type == "geom")
            dist = fj::FJDiscreteDist::Geometric;
        else if (type == "poisson" || type == "pois")
            dist = fj::FJDiscreteDist::Poisson;
        else
            throw UnsupportedError(ErrorKind{"Distribution"},
                "fj_char_max_discrete: dist_type must be 'geometric' or 'poisson', got '" + type +
                    "'");
        const fj::FJCharMaxDiscreteResult<T> r = fj::fj_char_max_discrete<T>(K, dist, par);
        Json out;
        out["MK"] = encode_scalar(r.MK);
        out["mK"] = static_cast<long long>(r.mK);
        out["exact"] = encode_scalar(r.exact);
        return out;
    }
};
 
struct OpFjCoxFit {
    template <class T>
    static Json run(Args& a) {
        const T m1 = a.number<T>("m1");
        const T c2 = a.number<T>("c2");
        a.done();
        const fj::FJCoxFitResult<T> r = fj::fj_cox_fit<T>(m1, c2);
        Json out;
        out["mu1"] = encode_scalar(r.mu1);
        out["mu2"] = encode_scalar(r.mu2);
        out["q"] = encode_scalar(r.q);
        out["kmin"] = static_cast<long long>(r.kmin);
        out["kmax"] = static_cast<long long>(r.kmax);
        return out;
    }
};
 
struct OpFjDagMakespan {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> pred = a.matrix<T>("pred");
        Matrix<T> rate = a.matrix<T>("rate");
        a.done();
        // MATLAB also takes an n-vector of constant rates and repmats it to the
        // n by n concurrency table; the port takes only the table.
        const std::size_t n = pred.rows();
        if ((rate.rows() == 1 || rate.cols() == 1) && rate.rows() * rate.cols() == n && n > 1) {
            Matrix<T> full(n, n);
            for (std::size_t i = 0; i < n; ++i)
                for (std::size_t k = 0; k < n; ++k)
                    full(i, k) = rate.rows() == 1 ? rate(0, i) : rate(i, 0);
            rate = full;
        }
        const fj::FJDagMakespanResult<T> r = fj::fj_dag_makespan<T>(pred, rate);
        Json out;
        out["C"] = encode_scalar(r.C);
        out["I"] = encode_vector(r.I);
        out["Cend"] = encode_vector(r.Cend);
        out["E"] = encode_vector(r.E);
        return out;
    }
};
 
/** MATLAB's options.npoints counts Simpson NODES, the port's npanels panels. */
unsigned panels_from_npoints(Args& o, const std::string& fn) {
    const unsigned np = o.uinteger("npoints");
    if (np < 2) throw InputError(fn + ": options.npoints must be at least 2");
    return np - 1;
}
 
/**
 * fj_delay_opt(shape, rate, options). The options struct carries maxsweeps,
 * dtol and npoints, taken positionally by the port in that order. options.tol
 * is refused: the port fixes the dispersion tolerance and takes no argument
 * for it.
 */
struct OpFjDelayOpt {
    template <class T>
    static Json run(Args& a) {
        const std::vector<unsigned> shape = read_uints(a, "shape", "fj_delay_opt");
        const std::vector<T> rate = a.vector<T>("rate");
        std::size_t n = 0;
        unsigned maxsweeps = 0, npanels = 0;
        T dtol = T();
        if (const Json* oj = a.opt("options")) {
            Args o(*oj, "fj_delay_opt: options");
            o.unsupported("tol", "the port fixes the dispersion tolerance of the inner "
                                 "fj_dispersion calls and takes no argument for it");
            n = prefix_count(o, {"maxsweeps", "dtol", "npoints"}, "fj_delay_opt: options");
            if (n >= 1) maxsweeps = o.uinteger("maxsweeps");
            if (n >= 2) dtol = o.number<T>("dtol");
            if (n >= 3) npanels = panels_from_npoints(o, "fj_delay_opt");
            o.done();
        }
        a.done();
        fj::FJDelayOptResult<T> r;
        if (n == 0) r = fj::fj_delay_opt<T>(shape, rate);
        else if (n == 1) r = fj::fj_delay_opt<T>(shape, rate, maxsweeps);
        else if (n == 2) r = fj::fj_delay_opt<T>(shape, rate, maxsweeps, dtol);
        else r = fj::fj_delay_opt<T>(shape, rate, maxsweeps, dtol, npanels);
        Json out;
        out["d"] = encode_vector(r.d);
        out["Edisp"] = encode_scalar(r.Edisp);
        out["Emax"] = encode_scalar(r.Emax);
        return out;
    }
};
 
/**
 * fj_dispersion(shape, rate, d, options). Without d the port's two-argument
 * overload runs; options {tol, npoints} are positional after d, so they
 * require d.
 */
struct OpFjDispersion {
    template <class T>
    static Json run(Args& a) {
        const std::vector<unsigned> shape = read_uints(a, "shape", "fj_dispersion");
        const std::vector<T> rate = a.vector<T>("rate");
        const bool hasD = a.has("d");
        const std::vector<T> d = a.vector_or_empty<T>("d");
        std::size_t n = 0;
        unsigned npanels = 0;
        T tol = T();
        if (const Json* oj = a.opt("options")) {
            if (!hasD)
                throw InputError("fj_dispersion: 'options' follows 'd' positionally in the port; "
                                 "give 'd' as well");
            Args o(*oj, "fj_dispersion: options");
            n = prefix_count(o, {"tol", "npoints"}, "fj_dispersion: options");
            if (n >= 1) tol = o.number<T>("tol");
            if (n >= 2) npanels = panels_from_npoints(o, "fj_dispersion");
            o.done();
        }
        a.done();
        fj::FJDispersionResult<T> r;
        if (!hasD) r = fj::fj_dispersion<T>(shape, rate);
        else if (n == 0) r = fj::fj_dispersion<T>(shape, rate, d);
        else if (n == 1) r = fj::fj_dispersion<T>(shape, rate, d, tol);
        else r = fj::fj_dispersion<T>(shape, rate, d, tol, npanels);
        Json out;
        out["Edisp"] = encode_scalar(r.Edisp);
        out["Emax"] = encode_scalar(r.Emax);
        out["Emin"] = encode_scalar(r.Emin);
        return out;
    }
};
 
struct OpFjIsmGreen {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T mu = a.number<T>("mu");
        const unsigned s = a.uinteger("s");
        const std::vector<T> c = a.vector<T>("c");
        a.done();
        const fj::FJIsmGreenResult<T> r = fj::fj_ism_green<T>(lambda, mu, s, c);
        Json o;
        o["EB"] = encode_scalar(r.EB);
        o["EB2"] = encode_scalar(r.EB2);
        o["ED"] = encode_scalar(r.ED);
        o["ED2"] = encode_scalar(r.ED2);
        o["EQ"] = encode_scalar(r.EQ);
        o["EQbar"] = encode_scalar(r.EQbar);
        o["pq"] = encode_scalar(r.pq);
        o["pd"] = encode_scalar(r.pd);
        o["pi0"] = encode_scalar(r.pi0);
        o["rho"] = encode_scalar(r.rho);
        o["q"] = encode_vector(r.q);
        o["ES"] = encode_scalar(r.ES);
        Json out;
        out["W"] = encode_scalar(r.W);
        out["R"] = encode_scalar(r.R);
        out["out"] = o;
        return out;
    }
};
 
struct OpFjLstMaxHet {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const T s = a.number<T>("s");
        a.done();
        Json out;
        out["L"] = encode_scalar(fj::fj_lst_max_het<T>(lambda, s));
        return out;
    }
};
 
struct OpFjResptBulk {
    template <class T>
    static Json run(Args& a) {
        const unsigned K = a.uinteger("K");
        const T lambda = a.number<T>("lambda");
        const T mu = a.number<T>("mu");
        const unsigned c = a.uinteger("c");
        const bool hasN = a.has("nmax");
        const std::size_t nmax = a.count("nmax", 0);
        a.done();
        const fj::FJResptBulkResult<T> r = hasN ? fj::fj_respt_bulk<T>(K, lambda, mu, c, nmax)
                                                : fj::fj_respt_bulk<T>(K, lambda, mu, c);
        Json out;
        out["Rreq"] = encode_scalar(r.Rreq);
        out["Rtask"] = encode_scalar(r.Rtask);
        out["Q"] = encode_scalar(r.Q);
        out["p"] = encode_vector(r.p);
        return out;
    }
};
 
struct OpFjResptClosed {
    template <class T>
    static Json run(Args& a) {
        const unsigned K = a.uinteger("K");
        const T x = a.number<T>("x");
        const unsigned M = a.uinteger("M");
        const bool hasA = a.has("A");
        const T A = hasA ? a.number<T>("A") : T();
        a.done();
        const fj::FJResptClosedResult<T> r =
            hasA ? fj::fj_respt_closed<T>(K, x, M, A) : fj::fj_respt_closed<T>(K, x, M);
        Json out;
        out["R"] = encode_scalar(r.R);
        out["exact"] = r.exact;
        return out;
    }
};
 
struct OpFjResptNosplit {
    template <class T>
    static Json run(Args& a) {
        const unsigned K = a.uinteger("K");
        const T lambda = a.number<T>("lambda");
        const T mu = a.number<T>("mu");
        a.done();
        const fj::FJResptNosplitResult<T> r = fj::fj_respt_nosplit<T>(K, lambda, mu);
        Json out;
        out["R"] = encode_scalar(r.R);
        out["rho"] = encode_scalar(r.rho);
        return out;
    }
};
 
struct OpFjSerialization {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> Rs = a.vector<T>("Rs");
        const T R0 = a.number<T>("R0");
        const unsigned M = a.uinteger("M");
        const bool hasAlpha = a.has("alpha");
        const T alpha = hasAlpha ? a.number<T>("alpha") : T();
        a.done();
        const fj::FJSerializationResult<T> r = hasAlpha
                                                   ? fj::fj_serialization<T>(Rs, R0, M, alpha)
                                                   : fj::fj_serialization<T>(Rs, R0, M);
        Json out;
        out["P"] = encode_vector(r.P);
        out["delay"] = encode_vector(r.delay);
        out["Rtot"] = encode_scalar(r.Rtot);
        return out;
    }
};
 
struct OpFjTsmCapacity {
    template <class T>
    static Json run(Args& a) {
        const unsigned s = a.uinteger("s");
        const std::vector<T> f = a.vector<T>("f");
        const std::vector<unsigned> rr = read_uints(a, "r", "fj_tsm_capacity");
        const std::vector<T> x = a.vector<T>("x");
        a.done();
        const fj::FJTsmCapacityResult<T> r = fj::fj_tsm_capacity<T>(s, f, rr, x);
        Json out;
        out["Lmax"] = encode_scalar(r.Lmax);
        out["Llp"] = encode_scalar(r.Llp);
        // MATLAB's Lfcfs is NaN outside the two-class closed form; omitted then.
        if (r.fcfs_available) out["Lfcfs"] = encode_scalar(r.Lfcfs);
        Json states = Json::array();
        for (std::size_t i = 0; i < r.states.size(); ++i) {
            Json row = Json::array();
            for (std::size_t k = 0; k < r.states[i].size(); ++k)
                row.push_back(static_cast<long long>(r.states[i][k]));
            states.push_back(row);
        }
        Json p;
        p["states"] = states;
        p["prob"] = encode_vector(r.prob);
        out["p"] = p;
        return out;
    }
};
 
struct OpFjXmaxCoxian {
    template <class T>
    static Json run(Args& a) {
        const unsigned K = a.uinteger("K");
        const T mu1 = a.number<T>("mu1");
        const T mu2 = a.number<T>("mu2");
        const T q = a.number<T>("q");
        a.done();
        const fj::FJXmaxCoxianResult<T> r = fj::fj_xmax_coxian<T>(K, mu1, mu2, q);
        Json out;
        out["Xmax"] = encode_scalar(r.Xmax);
        out["m1"] = encode_scalar(r.m1);
        out["c2"] = encode_scalar(r.c2);
        return out;
    }
};
 
struct OpFjXmaxHet {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const bool hasN = a.has("n");
        const unsigned n = a.uinteger("n", 0);
        a.done();
        Json out;
        out["Mn"] = encode_scalar(hasN ? fj::fj_xmax_het<T>(lambda, n) : fj::fj_xmax_het<T>(lambda));
        return out;
    }
};
 
struct OpFjXmaxMomentsHet {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const bool hasN = a.has("n");
        const unsigned n = a.uinteger("n", 0);
        a.done();
        Json out;
        out["M"] = encode_vector(hasN ? fj::fj_xmax_moments_het<T>(lambda, n)
                                      : fj::fj_xmax_moments_het<T>(lambda));
        return out;
    }
};
 
struct OpFjXmaxHz {
    template <class T>
    static Json run(Args& a) {
        const T m1 = a.number<T>("m1");
        const T m2 = a.number<T>("m2");
        const unsigned K = a.uinteger("K");
        a.done();
        const fj::FJXmaxHzResult<T> r = fj::fj_xmax_hz<T>(m1, m2, K);
        Json out;
        out["Xmax"] = encode_scalar(r.Xmax);
        out["resid"] = encode_scalar(r.resid);
        return out;
    }
};
 
struct OpFjMg1ResptMoments {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T ES = a.number<T>("ES");
        const T ES2 = a.number<T>("ES2");
        const T ES3 = a.number<T>("ES3");
        a.done();
        const fj::Mg1ResptMoments<T> r = fj::fj_mg1_respt_moments<T>(lambda, ES, ES2, ES3);
        Json out;
        out["ET"] = encode_scalar(r.ET);
        out["VT"] = encode_scalar(r.VT);
        return out;
    }
};
 
template <class T>
Json forktail_results(const fj::ForkTailResult<T>& r) {
    Json out;
    out["xp"] = encode_scalar(r.xp);
    out["alpha"] = encode_vector(r.alpha);
    out["beta"] = encode_vector(r.beta);
    return out;
}
 
/** fj_tail_forktail(ET, VT, K, p, P): K, p and P are positional. */
struct OpFjTailForktail {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> ET = a.vector<T>("ET");
        const std::vector<T> VT = a.vector<T>("VT");
        const std::size_t n = prefix_count(a, {"K", "p", "P"}, "fj_tail_forktail");
        const std::vector<T> K = n >= 1 ? a.vector<T>("K") : std::vector<T>();
        const T p = n >= 2 ? a.number<T>("p") : T();
        const std::vector<T> P = n >= 3 ? a.vector<T>("P") : std::vector<T>();
        a.done();
        if (n == 0) return forktail_results(fj::fj_tail_forktail<T>(ET, VT));
        if (n == 1) return forktail_results(fj::fj_tail_forktail<T>(ET, VT, K));
        if (n == 2) return forktail_results(fj::fj_tail_forktail<T>(ET, VT, K, p));
        return forktail_results(fj::fj_tail_forktail<T>(ET, VT, K, p, P));
    }
};
 
/** fj_tail_ordstat(ET, VT, K, p, kreq): K, p and kreq are positional. */
struct OpFjTailOrdstat {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> ET = a.vector<T>("ET");
        const std::vector<T> VT = a.vector<T>("VT");
        const std::size_t n = prefix_count(a, {"K", "p", "kreq"}, "fj_tail_ordstat");
        const std::size_t K = n >= 1 ? a.count("K") : 0;
        const T p = n >= 2 ? a.number<T>("p") : T();
        const std::size_t kreq = n >= 3 ? a.count("kreq") : 0;
        a.done();
        if (n == 0) return forktail_results(fj::fj_tail_ordstat<T>(ET, VT));
        if (n == 1) return forktail_results(fj::fj_tail_ordstat<T>(ET, VT, K));
        if (n == 2) return forktail_results(fj::fj_tail_ordstat<T>(ET, VT, K, p));
        return forktail_results(fj::fj_tail_ordstat<T>(ET, VT, K, p, kreq));
    }
};
 
/**
 * fj_branch_members(lqn, joinaidx). Of the LayeredNetworkStruct only the four
 * fields the algorithm reads cross the boundary, as the object lqn; indices
 * stay 1-based, as in MATLAB.
 */
struct OpFjBranchMembers {
    template <class T>
    static Json run(Args& a) {
        Args l(a.get("lqn"), "fj_branch_members: lqn");
        fj::LqnBranchView<T> v;
        v.graph = l.matrix<T>("graph");
        v.ashift = l.count("ashift");
        v.nacts = l.count("nacts");
        v.actposttype = l.ints("actposttype");
        l.done();
        const std::size_t joinaidx = a.count("joinaidx");
        a.done();
        const std::vector<std::vector<std::size_t> > m = fj::fj_branch_members<T>(v, joinaidx);
        Json members = Json::array();
        for (std::size_t i = 0; i < m.size(); ++i) members.push_back(encode_sizes(m[i]));
        Json out;
        out["members"] = members;
        return out;
    }
};
 
/**
 * fj_dist2fj(lineMAP, distType, sn, ist, r). The reference reads ONE field of
 * sn, procid(ist, r), so sn crosses as {"procid": [[..]]} and nothing else;
 * ist and r are 1-based. The codes are MATLAB's ProcessType values.
 */
struct OpFjDist2fj {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = map_from_pair<T>(a.get("lineMAP"), "fj_dist2fj: lineMAP");
        const std::string distType = a.text("distType", "");
        Args sn(a.get("sn"), "fj_dist2fj: sn");
        const Matrix<double> procid = sn.matrix<double>("procid");
        sn.done();
        const std::size_t ist = a.count("ist");
        const std::size_t r = a.count("r");
        a.done();
        fj::FjDistKind kind;
        if (distType == "arrival")
            kind = fj::FjDistKind::Arrival;
        else if (distType == "service")
            kind = fj::FjDistKind::Service;
        else
            throw InputError("fj_dist2fj: distType must be 'arrival' or 'service', got '" +
                             distType + "'");
        if (ist < 1 || r < 1 || ist > procid.rows() || r > procid.cols())
            throw InputError("fj_dist2fj: (ist, r) are 1-based and must index sn.procid");
        const double code = procid(ist - 1, r - 1);
        fj::FjProcType pt;
        if (code == 0.0) pt = fj::FjProcType::Exp;
        else if (code == 1.0) pt = fj::FjProcType::Erlang;
        else if (code == 2.0) pt = fj::FjProcType::HyperExp;
        else if (code == 5.0) pt = fj::FjProcType::Map;
        else
            throw UnsupportedError(ErrorKind{"Distribution"},
                "fj_dist2fj: sn.procid(ist, r) = " + std::to_string(code) +
                    " is not a process FJ_codes takes (EXP=0, ERLANG=1, HYPEREXP=2, MAP=5)");
        const fj::FjDist<T> d = fj::fj_dist2fj<T>(m, kind, pt);
        Json out;
        if (kind == fj::FjDistKind::Arrival) {
            out["lambda"] = encode_scalar(d.lambda);
            out["lambda0"] = encode_matrix(d.lambda0);
            out["lambda1"] = encode_matrix(d.lambda1);
            out["ma"] = encode_count(d.ma);
            out["Ia"] = encode_matrix(d.Ia);
            out["ArrChoice"] = d.choice;
        } else {
            out["mu"] = encode_scalar(d.mu);
            out["ST"] = encode_matrix(d.ST);
            out["St"] = encode_vector(d.St);
            out["tau_st"] = encode_vector(d.tau_st);
            out["SerChoice"] = d.choice;
        }
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// aoi: the service-law transforms, evaluated where the caller asks
// ---------------------------------------------------------------------------
//
// aoi_lst_* return a function handle in MATLAB. A handle cannot cross JSON, so
// the boundary takes the points "s" at which to evaluate it and returns the
// values as "lst". The law is the port's own constructor, unchanged.
 
template <class T>
Json lst_values(Args& a, const aoi::Lst<T>& L) {
    const std::vector<T> s = a.vector<T>("s");
    a.done();
    std::vector<T> v(s.size());
    for (std::size_t i = 0; i < s.size(); ++i) v[i] = L(s[i]);
    Json out;
    out["lst"] = encode_vector(v);
    return out;
}
 
struct OpAoiLstDet {
    template <class T>
    static Json run(Args& a) {
        const T d = a.number<T>("d");
        return lst_values<T>(a, aoi::aoi_lst_det<T>(d));
    }
};
 
struct OpAoiLstErlang {
    template <class T>
    static Json run(Args& a) {
        const unsigned k = a.uinteger("k");
        const T mu = a.number<T>("mu");
        return lst_values<T>(a, aoi::aoi_lst_erlang<T>(k, mu));
    }
};
 
struct OpAoiLstExp {
    template <class T>
    static Json run(Args& a) {
        const T mu = a.number<T>("mu");
        return lst_values<T>(a, aoi::aoi_lst_exp<T>(mu));
    }
};
 
struct OpAoiLstPh {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> alpha = a.vector<T>("alpha");
        const Matrix<T> Tm = a.matrix<T>("T");
        return lst_values<T>(a, aoi::aoi_lst_ph<T>(alpha, Tm));
    }
};
 
/** aoi_dist2ph(proc): proc is MATLAB's cell {D0, D1}. */
struct OpAoiDist2ph {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = map_from_pair<T>(a.get("proc"), "aoi_dist2ph: proc");
        a.done();
        const aoi::AoiPh<T> r = aoi::aoi_dist2ph<T>(m);
        Json out;
        out["alpha"] = encode_vector(r.alpha);
        out["T"] = encode_matrix(r.Tmat);
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// qsys: MAP/G/1/K and its marked and per-flow variants
// ---------------------------------------------------------------------------
 
/** The 'tol' and 'nmax' pair of the MAP/G/1/K family, which tune one truncation. */
template <class T>
struct Truncation {
    bool given = false;
    T tol = T();
    std::size_t nmax = 0;
};
 
template <class T>
Truncation<T> read_truncation(Args& a, const std::string& fn) {
    Truncation<T> tr;
    tr.given = all_or_none(a, {"tol", "nmax"}, fn);
    if (tr.given) {
        tr.tol = a.number<T>("tol");
        tr.nmax = a.count("nmax");
    } else {
        a.opt("tol");
        a.opt("nmax");
    }
    return tr;
}
 
struct OpQsysMapg1k {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> arrival = read_map<T>(a, "D0", "D1");
        const qsys::ServiceLaw<T> svc = read_svc<T>(a, "qsys_mapg1k");
        const std::size_t K = a.count("K");
        const Truncation<T> tr = read_truncation<T>(a, "qsys_mapg1k");
        a.done();
        const qsys::MapG1kResult<T> r = tr.given
                                            ? qsys::qsys_mapg1k<T>(arrival, svc, K, tr.tol, tr.nmax)
                                            : qsys::qsys_mapg1k<T>(arrival, svc, K);
        Json out;
        out["p0"] = encode_scalar(r.p0);
        out["pK"] = encode_scalar(r.pK);
        out["lossProbability"] = encode_scalar(r.lossProbability);
        out["throughput"] = encode_scalar(r.throughput);
        out["lambda"] = encode_scalar(r.lambda);
        out["meanServiceTime"] = encode_scalar(r.meanServiceTime);
        out["utilization"] = encode_scalar(r.utilization);
        out["rho"] = encode_scalar(r.rho);
        out["nmax"] = encode_count(r.nmax);
        out["sigma"] = encode_vector(r.sigma);
        out["pKvec"] = encode_vector(r.pKvec);
        out["p0vec"] = encode_vector(r.p0vec);
        out["plevel"] = encode_vector(r.plevel);
        out["meanQueueLength"] = encode_scalar(r.meanQueueLength);
        out["countingResidual"] = encode_scalar(r.countingResidual);
        return out;
    }
};
 
struct OpQsysMmapg1k {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> D0 = a.matrix<T>("D0");
        const std::vector<Matrix<T> > D1c = a.matrices<T>("D1c");
        const qsys::ServiceLaw<T> svc = read_svc<T>(a, "qsys_mmapg1k");
        const std::size_t K = a.count("K");
        const Truncation<T> tr = read_truncation<T>(a, "qsys_mmapg1k");
        a.done();
        const qsys::MmapG1kResult<T> r =
            tr.given ? qsys::qsys_mmapg1k<T>(D0, D1c, svc, K, tr.tol, tr.nmax)
                     : qsys::qsys_mmapg1k<T>(D0, D1c, svc, K);
        Json out;
        out["throughput"] = encode_vector(r.throughput);
        out["lossRatio"] = encode_vector(r.lossRatio);
        out["lambda"] = encode_vector(r.lambda);
        out["lambdaAggregate"] = encode_scalar(r.lambdaAggregate);
        out["throughputAggregate"] = encode_scalar(r.throughputAggregate);
        out["lossAggregate"] = encode_scalar(r.lossAggregate);
        out["p0"] = encode_scalar(r.p0);
        out["pK"] = encode_scalar(r.pK);
        out["pKvec"] = encode_vector(r.pKvec);
        out["plevel"] = encode_vector(r.plevel);
        out["meanQueueLength"] = encode_scalar(r.meanQueueLength);
        out["meanServiceTime"] = encode_scalar(r.meanServiceTime);
        out["utilization"] = encode_scalar(r.utilization);
        out["rho"] = encode_scalar(r.rho);
        return out;
    }
};
 
/** qsys_mapg1k_perflow(MAPS, svc, K): MAPS is a list of [D0, D1] pairs. */
struct OpQsysMapg1kPerflow {
    template <class T>
    static Json run(Args& a) {
        const std::vector<mam::Map<T> > flows =
            maps_from_json<T>(a.get("MAPS"), "qsys_mapg1k_perflow: MAPS");
        const qsys::ServiceLaw<T> svc = read_svc<T>(a, "qsys_mapg1k_perflow");
        const std::size_t K = a.count("K");
        const Truncation<T> tr = read_truncation<T>(a, "qsys_mapg1k_perflow");
        a.done();
        const qsys::MapG1kPerflowResult<T> r =
            tr.given ? qsys::qsys_mapg1k_perflow<T>(flows, svc, K, tr.tol, tr.nmax)
                     : qsys::qsys_mapg1k_perflow<T>(flows, svc, K);
        Json out;
        out["throughput"] = encode_vector(r.throughput);
        out["lossRatio"] = encode_vector(r.lossRatio);
        out["lambda"] = encode_vector(r.lambda);
        out["lambdaAggregate"] = encode_scalar(r.lambdaAggregate);
        out["throughputAggregate"] = encode_scalar(r.throughputAggregate);
        out["lossAggregate"] = encode_scalar(r.lossAggregate);
        out["p0"] = encode_vector(r.p0);
        out["pK"] = encode_vector(r.pK);
        out["meanServiceTime"] = encode_scalar(r.meanServiceTime);
        out["rho"] = encode_scalar(r.rho);
        return out;
    }
};
 
struct OpQsysSvclawCoeffs {
    template <class T>
    static Json run(Args& a) {
        const qsys::ServiceLaw<T> svc = read_svc<T>(a, "qsys_svclaw_coeffs");
        const T theta = a.number<T>("theta");
        const T tol = a.number<T>("tol");
        const std::size_t nmaxCap = a.count("nmaxCap");
        a.done();
        const qsys::SvclawCoeffsResult<T> r = qsys::qsys_svclaw_coeffs<T>(svc, theta, tol, nmaxCap);
        Json out;
        out["cn"] = encode_vector(r.cn);
        out["Smean"] = encode_scalar(r.Smean);
        out["residual"] = encode_scalar(r.residual);
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// qsys: laws given as distribution records
// ---------------------------------------------------------------------------
 
/**
 * qsys_ldps_workload(lambda, B, alpha, N, t, ngrid). B is MATLAB's
 * Distribution object; it crosses as its model.json record, and the port is
 * handed exactly the three things the reference reads from it: getMean,
 * getSCV and evalCDF. An empty or absent t selects the automatic grid, as in
 * MATLAB; the port takes t and ngrid together, so ngrid is required whenever
 * either is named.
 */
struct OpQsysLdpsWorkload {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const lang::Distrib<T> law = dist_record<T>(a.get("B"), "qsys_ldps_workload: B");
        const std::vector<T> alpha = a.vector<T>("alpha");
        const std::size_t N = a.count("N");
        const bool tuned = a.has("t") || a.has("ngrid");
        if (tuned && !a.has("ngrid"))
            throw InputError("qsys_ldps_workload: the port takes 't' and 'ngrid' together; give "
                             "'ngrid' with 't'");
        const std::vector<T> t = a.vector_or_empty<T>("t");
        const std::size_t ngrid = a.count("ngrid", 0);
        a.done();
        qsys::WorkloadServiceLaw<T> B;
        B.mean = lang::dist_moment(law, 1u);
        const T m2 = lang::dist_moment(law, 2u);
        B.scv = T(m2 / (B.mean * B.mean) - num_traits<T>::from_int(1));
        B.cdf = [law](const T& x) { return lang::dist_cdf(law, x); };
        const qsys::LdpsWorkloadResult<T> r =
            tuned ? qsys::qsys_ldps_workload<T>(lambda, B, alpha, N, t, ngrid)
                  : qsys::qsys_ldps_workload<T>(lambda, B, alpha, N);
        Json out;
        out["F"] = encode_vector(r.F);
        out["t"] = encode_vector(r.t);
        out["p"] = encode_vector(r.p);
        return out;
    }
};
 
/**
 * qsys_mmapgk1(MMAP, svc, ...). MMAP is the LINE list [D0, D1, D^(1), ...].
 * Each svc entry is a [D0s, D1s] pair (a phase-type MAP) or a distribution
 * record; MATLAB's third form, a struct holding an lst handle, cannot cross.
 * The port takes wPoints, numWMoms, tol and the iteration cap together (MATLAB
 * fixes the cap internally, the port exposes it as iterMax), and quadNodes
 * only after them.
 */
struct OpQsysMmapgk1 {
    template <class T>
    static Json run(Args& a) {
        const std::vector<Matrix<T> > MMAP = a.matrices<T>("MMAP");
        const Json& sj = a.get("svc");
        if (!sj.is_array()) throw InputError("qsys_mmapgk1: svc must be a list, one law per type");
        std::vector<qsys::MmapGk1Service<T> > svc;
        for (std::size_t k = 0; k < sj.size(); ++k) {
            const std::string where = "qsys_mmapgk1: svc[" + std::to_string(k) + "]";
            if (sj[k].is_array()) {
                svc.push_back(qsys::MmapGk1Service<T>(map_from_pair<T>(sj[k], where)));
            } else if (sj[k].is_object() && sj[k].contains("lst")) {
                throw UnsupportedError(ErrorKind{"JsonApi"},
                    where + ": a transform struct carries an lst function handle, which cannot "
                            "cross the JSON API boundary; give the law as a distribution record "
                            "or a [D0, D1] pair");
            } else {
                svc.push_back(qsys::MmapGk1Service<T>(dist_record<T>(sj[k], where)));
            }
        }
        const bool tuned = all_or_none(a, {"wPoints", "numWMoms", "tol", "iterMax"}, "qsys_mmapgk1");
        if (a.has("quadNodes") && !tuned)
            throw InputError("qsys_mmapgk1: 'quadNodes' follows 'wPoints', 'numWMoms', 'tol' and "
                             "'iterMax' positionally in the port; give them as well");
        const std::vector<T> w = a.vector_or_empty<T>("wPoints");
        const std::size_t nm = a.count("numWMoms", 0);
        const double tol = a.scalar<double>("tol", 0.0);
        const std::size_t iterMax = a.count("iterMax", 0);
        const bool hasQ = a.has("quadNodes");
        const std::size_t q = a.count("quadNodes", 0);
        a.done();
        qsys::MmapGk1Result<T> r;
        if (!tuned) r = qsys::qsys_mmapgk1<T>(MMAP, svc);
        else if (!hasQ) r = qsys::qsys_mmapgk1<T>(MMAP, svc, w, nm, tol, iterMax);
        else r = qsys::qsys_mmapgk1<T>(MMAP, svc, w, nm, tol, iterMax, q);
        Json out;
        out["lambda"] = encode_vector(r.lambdas);
        out["arrivalRate"] = encode_scalar(r.arrivalRate);
        out["utilization"] = encode_scalar(r.utilization);
        out["idleVector"] = encode_vector(r.idleVector);
        out["waitMoments"] = encode_vectors(r.waitMoments);
        out["meanWaitingTime"] = encode_vector(r.meanWaitingTime);
        out["meanSojournTime"] = encode_vector(r.meanSojournTime);
        out["meanQueueLength"] = encode_scalar(r.meanQueueLength);
        out["waitCDF"] = encode_vectors(r.waitCDF);
        out["waitPoints"] = encode_vector(r.waitPoints);
        return out;
    }
};
 
/**
 * qsys_mmapgk1_lcfspr(MMAP, svc, ...). MMAP and svc as for qsys_mmapgk1; the
 * MATLAB name-value options nMax, wPoints, quadNodes and tol are each optional.
 * respLST is a function handle and cannot cross JSON; as for aoi_lst_*, the
 * optional points "s" (real) say where to evaluate it, and "respLST" returns
 * its values there, [type][point]. respCDF carries the CDF values.
 */
struct OpQsysMmapgk1Lcfspr {
    template <class T>
    static Json run(Args& a) {
        const std::vector<Matrix<T> > MMAP = a.matrices<T>("MMAP");
        const Json& sj = a.get("svc");
        if (!sj.is_array())
            throw InputError("qsys_mmapgk1_lcfspr: svc must be a list, one law per type");
        std::vector<qsys::MmapGk1Service<T> > svc;
        for (std::size_t k = 0; k < sj.size(); ++k) {
            const std::string where = "qsys_mmapgk1_lcfspr: svc[" + std::to_string(k) + "]";
            if (sj[k].is_array()) {
                svc.push_back(qsys::MmapGk1Service<T>(map_from_pair<T>(sj[k], where)));
            } else if (sj[k].is_object() && sj[k].contains("lst")) {
                throw UnsupportedError(ErrorKind{"JsonApi"},
                    where + ": a transform struct carries an lst function handle, which cannot "
                            "cross the JSON API boundary; give the law as a distribution record "
                            "or a [D0, D1] pair");
            } else {
                svc.push_back(qsys::MmapGk1Service<T>(dist_record<T>(sj[k], where)));
            }
        }
        qsys::MmapGk1LcfsprOptions<T> opt;
        if (a.has("nMax")) opt.nMax = static_cast<long>(a.count("nMax"));
        opt.wPoints = a.vector_or_empty<T>("wPoints");
        opt.quadNodes = a.count("quadNodes", 120);
        opt.tol = a.scalar<double>("tol", 1e-12);
        const bool wantLst = a.has("s");
        const std::vector<T> sPts = wantLst ? a.vector<T>("s") : std::vector<T>();
        a.done();
        const qsys::MmapGk1LcfsprResult<T> r = qsys::qsys_mmapgk1_lcfspr<T>(MMAP, svc, opt);
        Json out;
        out["lambda"] = encode_vector(r.lambdas);
        out["arrivalRate"] = encode_scalar(r.arrivalRate);
        out["utilization"] = encode_scalar(r.utilization);
        out["Q"] = encode_matrix(r.Q);
        out["kappa"] = encode_vector(r.kappa);
        out["Rk"] = encode_matrices(r.Rk);
        out["R"] = encode_matrix(r.R);
        out["classMeanQueueLength"] = encode_vector(r.classMeanQueueLength);
        out["meanQueueLength"] = encode_scalar(r.meanQueueLength);
        out["meanSojournTime"] = encode_vector(r.meanSojournTime);
        out["queueLengthDistr"] = encode_vector(r.queueLengthDistr);
        out["classQueueLengthDistr"] = encode_vectors(r.classQueueLengthDistr);
        out["respCDF"] = encode_vectors(r.respCDF);
        out["respPoints"] = encode_vector(r.respPoints);
        if (wantLst) {
            const std::size_t K = r.lambdas.size();
            std::vector<std::vector<T> > lv(K, std::vector<T>(sPts.size()));
            for (std::size_t i = 0; i < sPts.size(); ++i) {
                const std::vector<std::complex<double> > w =
                    r.respLST(std::complex<double>(num_traits<T>::to_double(sPts[i]), 0.0));
                for (std::size_t k = 0; k < K; ++k) lv[k][i] = num_traits<T>::from_double(w[k].real());
            }
            out["respLST"] = encode_vectors(lv);
        }
        out["analyzer"] = r.analyzer;
        return out;
    }
};
 
/**
 * qsys_tandem_ub_ciucu(x, lst, p, mu, dlst). lst is the interarrival
 * Laplace-Stieltjes transform; it crosses as the distribution record of the
 * interarrival law, whose transform the port evaluates. dlst is refused: the
 * port differentiates lst itself when it is absent.
 */
struct OpQsysTandemUbCiucu {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> x = a.vector<T>("x");
        const lang::Distrib<T> law = dist_record<T>(a.get("lst"), "qsys_tandem_ub_ciucu: lst");
        const std::vector<T> p = a.vector<T>("p");
        const std::vector<T> mu = a.vector<T>("mu");
        a.unsupported("dlst", "it is a function handle; omit it and the port differentiates the "
                              "transform of lst numerically");
        a.done();
        const std::function<T(const T&)> lst = [law](const T& s) { return lang::dist_lst(law, s); };
        const qsys::TandemUbResult<T> r = qsys::qsys_tandem_ub_ciucu<T>(x, lst, p, mu);
        Json out;
        out["S"] = encode_vector(r.S);
        out["W"] = encode_vector(r.W);
        out["theta"] = encode_scalar(r.theta);
        out["alpha"] = encode_scalar(r.alpha);
        out["A"] = encode_scalar(r.A);
        out["B"] = encode_scalar(r.B);
        out["C"] = encode_scalar(r.C);
        out["D"] = encode_scalar(r.D);
        out["analyzer"] = r.analyzer;
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// qsys: Lindley recursions
// ---------------------------------------------------------------------------
 
template <class T>
Json lindley_results(const qsys::LindleyResult<T>& r, bool scalar) {
    Json out;
    out["mean"] = encode_shaped(r.mean, scalar);
    out["var"] = encode_shaped(r.var, scalar);
    out["moments"] = encode_matrix(r.moments);
    out["mmax"] = static_cast<long long>(r.mmax);
    out["analyzer"] = r.analyzer;
    return out;
}
 
struct OpQsysHh1Lindley {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const std::vector<T> pa = a.vector<T>("pa");
        const std::vector<T> mu = a.vector<T>("mu");
        const std::vector<T> ps = a.vector<T>("ps");
        const bool scalar = !a.get("Wn").is_array();
        const std::vector<T> Wn = a.vector<T>("Wn");
        const bool hasM = a.has("mmax");
        const unsigned mmax = a.uinteger("mmax", 0);
        a.done();
        if (scalar)
            return lindley_results(
                hasM ? qsys::qsys_hh1_lindley<T>(lambda, pa, mu, ps, Wn[0], mmax)
                     : qsys::qsys_hh1_lindley<T>(lambda, pa, mu, ps, Wn[0]),
                true);
        return lindley_results(hasM ? qsys::qsys_hh1_lindley<T>(lambda, pa, mu, ps, Wn, mmax)
                                    : qsys::qsys_hh1_lindley<T>(lambda, pa, mu, ps, Wn),
                               false);
    }
};
 
struct OpQsysMm1Lindley {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T mu = a.number<T>("mu");
        const bool scalar = !a.get("Wn").is_array();
        const std::vector<T> Wn = a.vector<T>("Wn");
        const bool hasM = a.has("mmax");
        const unsigned mmax = a.uinteger("mmax", 0);
        a.done();
        if (scalar)
            return lindley_results(hasM ? qsys::qsys_mm1_lindley<T>(lambda, mu, Wn[0], mmax)
                                        : qsys::qsys_mm1_lindley<T>(lambda, mu, Wn[0]),
                                   true);
        return lindley_results(hasM ? qsys::qsys_mm1_lindley<T>(lambda, mu, Wn, mmax)
                                    : qsys::qsys_mm1_lindley<T>(lambda, mu, Wn),
                               false);
    }
};
 
struct OpQsysLindleyMoment {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T mu = a.number<T>("mu");
        const bool scalar = !a.get("Wn").is_array();
        const std::vector<T> Wn = a.vector<T>("Wn");
        const unsigned m = a.uinteger("m");
        a.done();
        Json out;
        if (scalar)
            out["value"] = encode_scalar(qsys::qsys_lindley_moment<T>(lambda, mu, Wn[0], m));
        else
            out["value"] = encode_vector(qsys::qsys_lindley_moment<T>(lambda, mu, Wn, m));
        return out;
    }
};
 
struct OpQsysMm1TandemLindley {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T mu1 = a.number<T>("mu1");
        const T mu2 = a.number<T>("mu2");
        const bool s1 = !a.get("Wk").is_array();
        const bool s2 = !a.get("Wk1").is_array();
        const std::vector<T> Wk = a.vector<T>("Wk");
        const std::vector<T> Wk1 = a.vector<T>("Wk1");
        a.done();
        const bool scalar = s1 && s2;
        const qsys::Mm1TandemLindleyResult<T> r =
            scalar ? qsys::qsys_mm1_tandem_lindley<T>(lambda, mu1, mu2, Wk[0], Wk1[0])
                   : qsys::qsys_mm1_tandem_lindley<T>(lambda, mu1, mu2, Wk, Wk1);
        Json out;
        out["mean"] = encode_shaped(r.mean, scalar);
        out["interdepMean"] = encode_shaped(r.interdepMean, scalar);
        out["idleProb"] = encode_shaped(r.idleProb, scalar);
        out["analyzer"] = r.analyzer;
        return out;
    }
};
 
struct OpQsysTandemLindley {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> A = a.vector<T>("A");
        const Matrix<T> S = a.matrix<T>("S");
        const bool hasW0 = a.has("W0");
        const std::vector<T> W0 = a.vector_or_empty<T>("W0");
        a.done();
        const qsys::TandemLindleyResult<T> r =
            hasW0 ? qsys::qsys_tandem_lindley<T>(A, S, W0) : qsys::qsys_tandem_lindley<T>(A, S);
        Json out;
        out["W"] = encode_matrix(r.W);
        out["G"] = encode_matrix(r.G);
        out["T"] = encode_matrix(r.T_);
        out["departure"] = encode_matrix(r.departure);
        out["analyzer"] = r.analyzer;
        return out;
    }
};
 
struct OpQsysMm1Ps {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const std::vector<T> mu = a.vector<T>("mu");
        a.done();
        const qsys::Mm1PsResult<T> r = qsys::qsys_mm1_ps<T>(lambda, mu);
        Json out;
        out["W"] = encode_vector(r.W);
        out["W2"] = encode_vector(r.W2);
        out["alpha"] = encode_scalar(r.alpha);
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// qsys: the remaining stand-alone systems
// ---------------------------------------------------------------------------
 
/**
 * qsys_mg1_ps(lambda, svc, svcparam, ...), the phase-type path: svc is the
 * initial vector alpha and svcparam the subgenerator T. The transform-handle
 * path (svc a handle, svcparam the mean) and the 'pdf' option are handles and
 * cannot cross. The handle-valued outputs are not reported; their tabulations
 * on x, s and t are.
 */
struct OpQsysMg1Ps {
    static Json run(Args& a) {
        const double lambda = a.number<double>("lambda");
        if (a.get("svc").is_string())
            throw UnsupportedError(ErrorKind{"JsonApi"},
                "qsys_mg1_ps: argument 'svc' cannot cross the JSON API boundary as a transform "
                "handle; give the phase-type initial vector alpha, with svcparam its "
                "subgenerator T");
        const std::vector<double> alpha = a.vector<double>("svc");
        const Matrix<double> Tm = a.matrix<double>("svcparam");
        qsys::Mg1PsOptions opt;
        if (a.has("x")) opt.x = a.vector<double>("x");
        if (a.has("s")) opt.s = a.vector<double>("s");
        if (a.has("t")) opt.t = a.vector<double>("t");
        opt.nterms = a.count("nterms", opt.nterms);
        a.unsupported("pdf", "it is a density function handle; on the phase-type path the port "
                             "derives the density itself");
        a.done();
        const qsys::Mg1PsResult r = qsys::qsys_mg1_ps(lambda, alpha, Tm, opt);
        Json out;
        out["rho"] = r.rho;
        out["m1"] = r.m1;
        out["m2"] = r.m2;
        out["x"] = encode_vector(r.x);
        out["s"] = encode_vector(r.s);
        out["t"] = encode_vector(r.t);
        out["lstCondVal"] = encode_matrix(r.lstCondVal);
        out["lstUncondVal"] = encode_vector(r.lstUncondVal);
        out["atomCond"] = encode_vector(r.atomCond);
        out["atomUncond"] = r.atomUncond;
        out["meanCond"] = encode_vector(r.meanCond);
        out["m2Cond"] = encode_vector(r.m2Cond);
        out["varCond"] = encode_vector(r.varCond);
        out["meanUncond"] = r.meanUncond;
        out["m2Uncond"] = r.m2Uncond;
        out["varUncond"] = r.varUncond;
        out["pdfCond"] = encode_matrix(r.pdfCond);
        out["cdfCond"] = encode_matrix(r.cdfCond);
        out["pdfUncond"] = encode_vector(r.pdfUncond);
        out["cdfUncond"] = encode_vector(r.cdfUncond);
        return out;
    }
};
 
/**
 * qsys_bmapm1(D, mu, ...). The port's tuned overload takes all five options;
 * Uniformization and MaxLevel have the port's own "select automatically"
 * values (<= 0 and 0), which are what MATLAB's empty defaults mean, while
 * MaxIter, Tolerance and TailTolerance have numeric defaults that are not
 * restated: naming any option therefore requires those three.
 */
struct OpQsysBmapm1 {
    template <class T>
    static Json run(Args& a) {
        const std::vector<Matrix<T> > D = a.matrices<T>("D");
        const T mu = a.number<T>("mu");
        const bool tuned = a.has("Uniformization") || a.has("MaxLevel") || a.has("MaxIter") ||
                           a.has("Tolerance") || a.has("TailTolerance");
        if (tuned && !(a.has("MaxIter") && a.has("Tolerance") && a.has("TailTolerance")))
            throw InputError("qsys_bmapm1: the port takes its options together; with any of them "
                             "give 'MaxIter', 'Tolerance' and 'TailTolerance'");
        const T q = a.scalar<T>("Uniformization", num_traits<T>::from_int(0));
        const std::size_t maxLevel = a.count("MaxLevel", 0);
        const unsigned maxIter = a.uinteger("MaxIter", 0);
        const T tol = a.scalar<T>("Tolerance", num_traits<T>::from_int(0));
        const double tailTol = a.scalar<double>("TailTolerance", 0.0);
        a.done();
        const qsys::BmapM1Result<T> r =
            tuned ? qsys::qsys_bmapm1<T>(D, mu, q, maxLevel, maxIter, tol, tailTol)
                  : qsys::qsys_bmapm1<T>(D, mu);
        Json out;
        out["theta"] = encode_vector(r.theta);
        out["lambda"] = encode_scalar(r.lambda);
        out["rho"] = encode_scalar(r.rho);
        out["q"] = encode_scalar(r.q);
        out["A0"] = encode_matrix(r.A0);
        out["A1"] = encode_matrix(r.A1);
        out["B0"] = encode_matrix(r.B0);
        out["Bk"] = encode_matrices(r.Bk);
        out["A"] = encode_matrix(r.A);
        out["alpha"] = encode_vector(r.alpha);
        out["G"] = encode_matrix(r.G);
        out["drift"] = encode_scalar(r.drift);
        out["decayRate"] = r.decayRate;
        out["levelProb"] = encode_matrix(r.levelProb);
        out["pi0"] = encode_scalar(r.pi0);
        out["meanQueueLength"] = encode_scalar(r.meanQueueLength);
        out["utilization"] = encode_scalar(r.utilization);
        out["throughput"] = encode_scalar(r.throughput);
        out["truncLevel"] = encode_count(r.truncLevel);
        out["truncError"] = r.truncError;
        out["gConverged"] = r.gConverged;
        return out;
    }
};
 
/** qsys_maxima_twomoment(n, mean_, cs2, 'q', 'exactFitted'), positional in the port. */
struct OpQsysMaximaTwomoment {
    template <class T>
    static Json run(Args& a) {
        const std::size_t n = a.count("n");
        const T mean = a.number<T>("mean_");
        const T cs2 = a.number<T>("cs2");
        const std::size_t k = prefix_count(a, {"q", "exactFitted"}, "qsys_maxima_twomoment");
        const T q = k >= 1 ? a.number<T>("q") : T();
        const bool ef = k >= 2 ? a.boolean("exactFitted", false) : false;
        a.done();
        const qsys::MaximaResult<T> r =
            k == 0 ? qsys::qsys_maxima_twomoment<T>(n, mean, cs2)
                   : (k == 1 ? qsys::qsys_maxima_twomoment<T>(n, mean, cs2, q)
                             : qsys::qsys_maxima_twomoment<T>(n, mean, cs2, q, ef));
        Json out;
        out["value"] = encode_scalar(r.value);
        out["slope"] = encode_scalar(r.slope);
        out["eta"] = encode_scalar(r.eta);
        out["threshold"] = encode_scalar(r.threshold);
        out["reliable"] = r.reliable;
        out["family"] = r.family;
        if (r.hasFitted) out["exactFittedValue"] = encode_scalar(r.exactFittedValue);
        return out;
    }
};
 
/**
 * qsys_mgisrgi_whitt(lambda, mu, s, r, patience, ...). r may be the string
 * "Inf". A numeric patience is the exponential rate theta, as in MATLAB; the
 * hazard and ccdf forms are function handles and cannot cross. The options
 * keep the port's MgisrgiOptions defaults unless named.
 */
struct OpQsysMgisrgiWhitt {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T mu = a.number<T>("mu");
        const unsigned s = a.uinteger("s");
        const double r = a.number<double>("r");
        const Json& pj = a.get("patience");
        if (pj.is_object() || pj.is_string())
            throw UnsupportedError(ErrorKind{"JsonApi"},
                "qsys_mgisrgi_whitt: argument 'patience' cannot cross the JSON API boundary as a "
                "hazard or ccdf function handle; give the exponential rate theta as a number");
        const qsys::Patience<T> patience =
            qsys::Patience<T>::exponential(a.number<T>("patience"));
        qsys::MgisrgiOptions opt;
        if (a.has("wPoints")) opt.wPoints = a.vector<double>("wPoints");
        opt.maxQueue = a.count("maxQueue", opt.maxQueue);
        opt.tol = a.scalar<double>("tol", opt.tol);
        opt.invMethod = a.text("invMethod", opt.invMethod);
        opt.invN = a.count("invN", opt.invN);
        a.done();
        const qsys::QsysAbandonResult<T> res =
            qsys::qsys_mgisrgi_whitt<T>(lambda, mu, s, r, patience, opt);
        Json out;
        out["queueLengthDist"] = encode_vector(res.queueLengthDist);
        out["probLoss"] = encode_scalar(res.probLoss);
        out["probNoWait"] = encode_scalar(res.probNoWait);
        out["probServed"] = encode_scalar(res.probServed);
        out["probAbandon"] = encode_scalar(res.probAbandon);
        out["meanNumber"] = encode_scalar(res.meanNumber);
        out["varNumber"] = encode_scalar(res.varNumber);
        out["meanQueueLength"] = encode_scalar(res.meanQueueLength);
        out["varQueueLength"] = encode_scalar(res.varQueueLength);
        out["utilization"] = encode_scalar(res.utilization);
        out["throughput"] = encode_scalar(res.throughput);
        out["abandonRate"] = encode_scalar(res.abandonRate);
        out["meanWaitServed"] = encode_scalar(res.meanWaitServed);
        out["varWaitServed"] = encode_scalar(res.varWaitServed);
        out["meanWaitAbandon"] = encode_scalar(res.meanWaitAbandon);
        out["varWaitAbandon"] = encode_scalar(res.varWaitAbandon);
        out["meanWait"] = encode_scalar(res.meanWait);
        out["secondMomentWait"] = encode_scalar(res.secondMomentWait);
        out["abandonRates"] = encode_vector(res.abandonRates);
        out["totalAbandonRates"] = encode_vector(res.totalAbandonRates);
        out["numWaitingSpaces"] = encode_count(res.numWaitingSpaces);
        out["exponentialPatience"] = res.exponentialPatience;
        out["patienceRate"] = encode_scalar(res.patienceRate);
        out["waitPoints"] = encode_vector(res.waitPoints);
        out["cdfWaitServed"] = encode_vector(res.cdfWaitServed);
        out["cdfWaitAbandon"] = encode_vector(res.cdfWaitAbandon);
        out["cdfWait"] = encode_vector(res.cdfWait);
        return out;
    }
};
 
/** qsys_mapmapc(D0, D1, S0, S1, c, 'maxNumComp', 'numWMoms'), positional in the port. */
struct OpQsysMapmapc {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> arrival = read_map<T>(a, "D0", "D1");
        const mam::Map<T> service = read_map<T>(a, "S0", "S1");
        const unsigned c = a.uinteger("c");
        const std::size_t k = prefix_count(a, {"maxNumComp", "numWMoms"}, "qsys_mapmapc");
        const std::size_t dsz = k >= 1 ? a.count("maxNumComp") : 0;
        const std::size_t nw = k >= 2 ? a.count("numWMoms") : 0;
        a.done();
        const qsys::MapPhcResult<T> r =
            k == 0 ? qsys::qsys_mapmapc<T>(arrival, service, c)
                   : (k == 1 ? qsys::qsys_mapmapc<T>(arrival, service, c, dsz)
                             : qsys::qsys_mapmapc<T>(arrival, service, c, dsz, nw));
        Json out;
        out["meanQueueLength"] = encode_scalar(r.meanQueueLength);
        out["meanWaitingTime"] = encode_scalar(r.meanWaitingTime);
        out["meanSojournTime"] = encode_scalar(r.meanSojournTime);
        out["utilization"] = encode_scalar(r.utilization);
        out["queueLengthDist"] = encode_vector(r.queueLengthDist);
        out["waitingTimeMoments"] = encode_vector(r.waitingTimeMoments);
        out["waitingTimeCCDF"] = encode_vector(r.waitingTimeCCDF);
        out["waitingTimePoints"] = encode_vector(r.waitingTimePoints);
        out["probWait"] = encode_scalar(r.probWait);
        out["phaseCount"] = encode_count(r.phaseCount);
        return out;
    }
};
 
/**
 * qsys_mmapmap1(D, svc, 'carry'). D is MATLAB's BUTools-layout cell
 * {D0, D1_1, ..., D1_K}: the per-mark blocks, whose sum is the aggregate D1.
 * svc is a list of [S0, S1] pairs, one per mark.
 */
struct OpQsysMmapmap1 {
    template <class T>
    static Json run(Args& a) {
        const std::vector<Matrix<T> > D = a.matrices<T>("D");
        const std::vector<mam::Map<T> > svc = maps_from_json<T>(a.get("svc"), "qsys_mmapmap1: svc");
        const bool hasCarry = a.has("carry");
        const std::vector<int> ci = a.ints_or_empty("carry");
        a.done();
        if (D.size() < 2)
            throw InputError("qsys_mmapmap1: D must be {D0, D1_1, ..., D1_K} with K >= 1");
        mam::Mmap<T> mm;
        mm.D0 = D[0];
        mm.D1 = D[1];
        mm.Dc.push_back(D[1]);
        for (std::size_t k = 2; k < D.size(); ++k) {
            if (D[k].rows() != mm.D1.rows() || D[k].cols() != mm.D1.cols())
                throw InputError("qsys_mmapmap1: the blocks of D must share one size");
            for (std::size_t i = 0; i < mm.D1.rows(); ++i)
                for (std::size_t j = 0; j < mm.D1.cols(); ++j)
                    mm.D1(i, j) = T(mm.D1(i, j) + D[k](i, j));
            mm.Dc.push_back(D[k]);
        }
        std::vector<bool> carry;
        for (std::size_t i = 0; i < ci.size(); ++i) {
            if (ci[i] != 0 && ci[i] != 1)
                throw InputError("qsys_mmapmap1: 'carry' holds flags, 0 or 1");
            carry.push_back(ci[i] != 0);
        }
        const qsys::MmapMap1Result<T> r = hasCarry ? qsys::qsys_mmapmap1<T>(mm, svc, carry)
                                                   : qsys::qsys_mmapmap1<T>(mm, svc);
        Json out;
        out["meanQueueLength"] = encode_vector(r.meanQueueLength);
        out["meanSojournTime"] = encode_vector(r.meanSojournTime);
        out["arrivalRate"] = encode_vector(r.arrivalRate);
        Json cj = Json::array();
        for (std::size_t i = 0; i < r.carry.size(); ++i) cj.push_back(static_cast<bool>(r.carry[i]));
        out["carry"] = cj;
        out["markCount"] = encode_count(r.markCount);
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// sim: planning and option helpers
// ---------------------------------------------------------------------------
 
/** sim_quest_options(options): the MATLAB options struct, validated and completed. */
struct OpSimQuestOptions {
    static Json run(Args& a) {
        sim::QuestOptions in;
        const Json* oj = a.opt("options");
        if (oj) {
            Args o(*oj, "sim_quest_options: options");
            in.b0 = static_cast<long>(o.integer("b0", static_cast<int>(in.b0)));
            in.m0 = static_cast<long>(o.integer("m0", static_cast<int>(in.m0)));
            if (o.has("s")) {
                const std::vector<int> s = o.ints("s");
                in.s.assign(s.begin(), s.end());
            }
            in.beta = o.scalar<double>("beta", in.beta);
            in.eta = o.scalar<double>("eta", in.eta);
            in.theta = o.scalar<double>("theta", in.theta);
            in.weight = o.scalar<double>("weight", in.weight);
            in.force = o.boolean("force", in.force);
            o.done();
        }
        a.done();
        const sim::QuestOptions r = oj ? sim::sim_quest_options(in) : sim::sim_quest_options();
        Json opt;
        opt["b0"] = static_cast<long long>(r.b0);
        opt["m0"] = static_cast<long long>(r.m0);
        Json s = Json::array();
        for (std::size_t i = 0; i < r.s.size(); ++i) s.push_back(static_cast<long long>(r.s[i]));
        opt["s"] = s;
        opt["beta"] = r.beta;
        opt["eta"] = r.eta;
        opt["theta"] = r.theta;
        opt["weight"] = r.weight;
        opt["force"] = r.force;
        Json out;
        out["opt"] = opt;
        return out;
    }
};
 
struct OpSimQuestHeuristicCi {
    template <class T>
    static Json run(Args& a) {
        if constexpr (!num_traits<T>::has_transcendental) {
            throw UnsupportedError(ErrorKind{"Arithmetic"}, sim_exact_refusal(
                "sim_quest_heuristic_ci", "the interval is a quantile times a square root"));
        } else {
            const std::vector<T> bqe = a.vector<T>("bqe");
            const T centre = a.number<T>("centre");
            const T Ap = a.number<T>("Ap");
            const T Np = a.number<T>("Np");
            const std::size_t nstar = a.count("nstar");
            const double alpha = a.number<double>("alpha");
            const bool ac = required_flag(a, "useAutocorr", "sim_quest_heuristic_ci");
            a.done();
            const sim::QuestInterval<T> r =
                sim::sim_quest_heuristic_ci<T>(bqe, centre, Ap, Np, nstar, alpha, ac);
            Json out;
            out["lower"] = encode_scalar(r.lower);
            out["upper"] = encode_scalar(r.upper);
            return out;
        }
    }
};
 
/** sim_runlength(mean_, asymVar, 'relPrecision', 'confidence', 'runLength'), positional. */
struct OpSimRunlength {
    template <class T>
    static Json run(Args& a) {
        if constexpr (!num_traits<T>::has_transcendental) {
            throw UnsupportedError(ErrorKind{"Arithmetic"}, sim_exact_refusal(
                "sim_runlength", "the run length is set by a normal quantile"));
        } else {
            const T mean = a.number<T>("mean_");
            const T av = a.number<T>("asymVar");
            const std::size_t k =
                prefix_count(a, {"relPrecision", "confidence", "runLength"}, "sim_runlength");
            const T rp = k >= 1 ? a.number<T>("relPrecision") : T();
            const T cf = k >= 2 ? a.number<T>("confidence") : T();
            const T rl = k >= 3 ? a.number<T>("runLength") : T();
            a.done();
            sim::RunLengthResult<T> r;
            if (k == 0) r = sim::sim_runlength<T>(mean, av);
            else if (k == 1) r = sim::sim_runlength<T>(mean, av, rp);
            else if (k == 2) r = sim::sim_runlength<T>(mean, av, rp, cf);
            else r = sim::sim_runlength<T>(mean, av, rp, cf, rl);
            Json out;
            out["requiredRunLength"] = encode_scalar(r.requiredRunLength);
            out["z"] = encode_scalar(r.z);
            if (r.hasRun) {
                out["halfWidth"] = encode_scalar(r.halfWidth);
                out["achievedRelPrecision"] = encode_scalar(r.achievedRelPrecision);
            }
            return out;
        }
    }
};
 
// ---------------------------------------------------------------------------
// Exact closed forms: finite-source, loss, priority, busy-period, setup and
// renewal-arrival queues. A law argument (svc, LA, EU, lst) is the model.json
// distribution record; MATLAB's LST-handle form has no JSON spelling.
// ---------------------------------------------------------------------------
 
/** The refusal of an LST function handle where the reference also takes one. */
[[noreturn]] void refuse_lst_handle(const std::string& fn, const char* key) {
    throw UnsupportedError(ErrorKind{"JsonApi"},
        fn + ": argument '" + key + "' cannot cross the JSON API boundary as an LST function "
             "handle; give it as a distribution record {\"type\": ..., ...}");
}
 
/** qsys_mmckk(alpha, mu, c, K, 'tPoints', 'percentiles'): M/M/c/K/K. */
struct OpQsysMmckk {
    static Json run(Args& a) {
        const double alpha = a.number<double>("alpha");
        const double mu = a.number<double>("mu");
        const unsigned c = a.uinteger("c");
        const unsigned K = a.uinteger("K");
        const qsys::QsysResOptions o = read_qsys_res_options(a, false);
        a.done();
        return encode_qsys_res(qsys::qsys_mmckk(alpha, mu, c, K, o));
    }
};
 
/** qsys_engset(alpha, mu, c, K): the Engset loss system. */
struct OpQsysEngset {
    static Json run(Args& a) {
        const double alpha = a.number<double>("alpha");
        const double mu = a.number<double>("mu");
        const unsigned c = a.uinteger("c");
        const unsigned K = a.uinteger("K");
        a.done();
        return encode_qsys_res(qsys::qsys_engset(alpha, mu, c, K));
    }
};
 
/** qsys_mg1kk(alpha, K, svc[, ES]): the finite-source M/G/1//K queue (Takacs). */
struct OpQsysMg1kk {
    static Json run(Args& a) {
        const double alpha = a.number<double>("alpha");
        const unsigned K = a.uinteger("K");
        if (a.get("svc").is_string()) refuse_lst_handle("qsys_mg1kk", "svc");
        const lang::Distrib<double> svc = dist_record<double>(a.get("svc"), "qsys_mg1kk: svc");
        // With a Distribution the reference reads the mean from the law and ignores ES.
        if (a.has("ES")) (void)a.number<double>("ES");
        a.done();
        return encode_qsys_res(qsys::qsys_mg1kk(alpha, K, svc));
    }
};
 
/** qsys_ddckk(EO, ES, c, K): the deterministic finite-source D/D/c/K/K queue. */
struct OpQsysDdckk {
    template <class T>
    static Json run(Args& a) {
        const T EO = a.number<T>("EO");
        const T ES = a.number<T>("ES");
        const unsigned c = a.uinteger("c");
        const unsigned K = a.uinteger("K");
        a.done();
        return encode_qsys_res(qsys::qsys_ddckk(EO, ES, c, K));
    }
};
 
/** qsys_mg1_setup(lambda, mu, cs, EU[, EU2]): M/G/1 with setup times; EU a number or a law. */
struct OpQsysMg1Setup {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T mu = a.number<T>("mu");
        const T cs = a.number<T>("cs");
        qsys::Mg1SetupResult<T> r;
        if (a.get("EU").is_object()) {
            const lang::Distrib<T> U = dist_record<T>(a.get("EU"), "qsys_mg1_setup: EU");
            a.done();
            r = qsys::qsys_mg1_setup(lambda, mu, cs, U);
        } else {
            const T EU = a.number<T>("EU");
            if (!a.has("EU2"))
                throw InputError(ErrorKind("Argument"),
                    "qsys_mg1_setup: qsys_mg1_setup needs the second moment EU2 of the setup time "
                    "after EU, or a setup Distribution in place of EU");
            const T EU2 = a.number<T>("EU2");
            a.done();
            r = qsys::qsys_mg1_setup(lambda, mu, cs, EU, EU2);
        }
        Json out;
        out["W"] = encode_scalar(r.W);
        out["Wq"] = encode_scalar(r.Wq);
        out["U"] = encode_scalar(r.U);
        out["Q"] = encode_scalar(r.Q);
        out["res"] = encode_qsys_res(r.res);
        return out;
    }
};
 
/** qsys_mg1_exceptional(lambda, mu, cs, muf, csf): M/G/1 with exceptional first service. */
struct OpQsysMg1Exceptional {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T mu = a.number<T>("mu");
        const T cs = a.number<T>("cs");
        const T muf = a.number<T>("muf");
        const T csf = a.number<T>("csf");
        a.done();
        const qsys::Mg1ExceptionalResult<T> r = qsys::qsys_mg1_exceptional(lambda, mu, cs, muf, csf);
        Json out;
        out["W"] = encode_scalar(r.W);
        out["Wq"] = encode_scalar(r.Wq);
        out["U"] = encode_scalar(r.U);
        out["Q"] = encode_scalar(r.Q);
        out["res"] = encode_qsys_res(r.res);
        return out;
    }
};
 
/** qsys_mmcc(lambda, mu, c): the Erlang-B loss system. */
struct OpQsysMmcc {
    static Json run(Args& a) {
        const double lambda = a.number<double>("lambda");
        const double mu = a.number<double>("mu");
        const unsigned c = a.uinteger("c");
        a.done();
        return encode_qsys_res(qsys::qsys_mmcc(lambda, mu, c));
    }
};
 
/** qsys_mgc_ps(lambda, mu, c): the insensitive M/G/c processor-sharing queue. */
struct OpQsysMgcPs {
    static Json run(Args& a) {
        const double lambda = a.number<double>("lambda");
        const double mu = a.number<double>("mu");
        const unsigned c = a.uinteger("c");
        a.done();
        return encode_qsys_res(qsys::qsys_mgc_ps(lambda, mu, c));
    }
};
 
/** qsys_mg1_prio_pr(lambda, mu, cs, 'ES3'): preemptive-resume priority M/G/1. */
struct OpQsysMg1PrioPr {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const std::vector<T> mu = a.vector<T>("mu");
        const std::vector<T> cs = a.vector<T>("cs");
        const std::vector<T> ES3 = a.vector_or_empty<T>("ES3");
        a.done();
        const qsys::Mg1PrioPrResult<T> r = qsys::qsys_mg1_prio_pr(lambda, mu, cs, ES3);
        Json out;
        out["W"] = encode_vector(r.W);
        out["rhohat"] = encode_scalar(r.rho);
        out["res"] = encode_qsys_res(r.res);
        return out;
    }
};
 
/** qsys_mmc_prio(lambda, mu, c): non-preemptive priority M/M/c, common service rate. */
struct OpQsysMmcPrio {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const T mu = a.number<T>("mu");
        const unsigned c = a.uinteger("c");
        a.done();
        const qsys::MmcPrioResult<T> r = qsys::qsys_mmc_prio(lambda, mu, c);
        Json out;
        out["W"] = encode_vector(r.W);
        out["rhohat"] = encode_scalar(r.rho);
        out["res"] = encode_qsys_res(r.res);
        return out;
    }
};
 
/** qsys_mmc_prio_tardiness(lambda, mu, c, d[, type, method]): exact priority M/M/c tardiness (HGK 2013). */
struct OpQsysMmcPrioTardiness {
    static Json run(Args& a) {
        const std::vector<double> lambda = a.vector<double>("lambda");
        const double mu = a.number<double>("mu");
        const unsigned c = a.uinteger("c");
        const std::vector<double> d = a.vector<double>("d");
        const std::string type = a.text("type", "np");
        const std::string method = a.text("method", "euler");
        a.done();
        const qsys::MmcPrioTardinessResult r = qsys::qsys_mmc_prio_tardiness(lambda, mu, c, d, type, method);
        Json out;
        out["Tard"] = encode_vector(r.Tard);
        out["Ptard"] = encode_vector(r.Ptard);
        out["res"] = encode_qsys_res(r.res);
        return out;
    }
};
 
/** qsys_mmc_prio_eta(lambda, mu, c, r, v, q, d[, type, method]): expected tardiness at arrival (HGK 2013). */
struct OpQsysMmcPrioEta {
    static Json run(Args& a) {
        const std::vector<double> lambda = a.vector<double>("lambda");
        const double mu = a.number<double>("mu");
        const unsigned c = a.uinteger("c");
        const unsigned r = a.uinteger("r");
        const unsigned v = a.uinteger("v");
        const unsigned q = a.uinteger("q");
        const std::vector<double> d = a.vector<double>("d");
        const std::string type = a.text("type", "np");
        const std::string method = a.text("method", "euler");
        a.done();
        const qsys::MmcPrioEtaResult e = qsys::qsys_mmc_prio_eta(lambda, mu, c, r, v, q, d, type, method);
        Json out;
        out["Tard"] = encode_vector(e.Tard);
        out["Ptard"] = encode_vector(e.Ptard);
        out["S"] = encode_scalar(e.S);
        return out;
    }
};
 
/** qsys_mg1_conservation(lambda, mu, cs[, W]): Kleinrock's conservation law. */
struct OpQsysMg1Conservation {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const std::vector<T> mu = a.vector<T>("mu");
        const std::vector<T> cs = a.vector<T>("cs");
        const std::vector<T> W = a.vector_or_empty<T>("W");
        a.done();
        const qsys::Mg1ConservationResult<T> r = qsys::qsys_mg1_conservation(lambda, mu, cs, W);
        Json res = Json::object();
        res["invariant"] = encode_scalar(r.V);
        res["fcfsWaitingTime"] = encode_scalar(r.fcfsWaitingTime);
        // MATLAB's NaN when W is not given; JSON has no NaN, so null.
        res["residual"] = r.hasResidual ? encode_scalar(r.residual) : Json();
        Json out;
        out["V"] = encode_scalar(r.V);
        out["res"] = res;
        return out;
    }
};
 
/** qsys_mg1_sjf(lambda, mu, cs): non-preemptive shortest-job-first M/G/1. */
struct OpQsysMg1Sjf {
    static Json run(Args& a) {
        const std::vector<double> lambda = a.vector<double>("lambda");
        const std::vector<double> mu = a.vector<double>("mu");
        const std::vector<double> cs = a.vector<double>("cs");
        a.done();
        const qsys::Mg1SjfResult r = qsys::qsys_mg1_sjf(lambda, mu, cs);
        Json out;
        out["W"] = encode_vector(r.W);
        out["rhohat"] = r.rhohat;
        out["res"] = encode_qsys_res(r.res);
        return out;
    }
};
 
/** qsys_mg1_busy(lambda, mu, cs, 'ES3', 'lst', 'sPoints'): the M/G/1 busy period. */
struct OpQsysMg1Busy {
    static Json run(Args& a) {
        const double lambda = a.number<double>("lambda");
        const double mu = a.number<double>("mu");
        const double cs = a.number<double>("cs");
        qsys::Mg1BusyOptions o;
        if (a.has("ES3")) o.ES3 = a.number<double>("ES3");
        if (a.has("lst")) {
            if (a.get("lst").is_string()) refuse_lst_handle("qsys_mg1_busy", "lst");
            const lang::Distrib<double> d = dist_record<double>(a.get("lst"), "qsys_mg1_busy: lst");
            o.lst = [d](double s) { return lang::dist_lst(d, s); };
        }
        if (a.has("sPoints")) o.sPoints = a.vector<double>("sPoints");
        a.done();
        return encode_qsys_res(qsys::qsys_mg1_busy(lambda, mu, cs, o));
    }
};
 
/** qsys_gimc(LA, mu, c, 'lambda', 'tPoints', 'percentiles', 'nmax'): the GI/M/c queue. */
struct OpQsysGimc {
    static Json run(Args& a) {
        if (a.get("LA").is_string()) refuse_lst_handle("qsys_gimc", "LA");
        const lang::Distrib<double> LA = dist_record<double>(a.get("LA"), "qsys_gimc: LA");
        const double mu = a.number<double>("mu");
        const unsigned c = a.uinteger("c");
        const double lambda = a.has("lambda") ? a.number<double>("lambda") : -1.0;
        const qsys::QsysResOptions o = read_qsys_res_options(a);
        a.done();
        return encode_qsys_res(qsys::qsys_gimc(LA, mu, c, o, lambda));
    }
};
 
/** qsys_md1(lambda, mu, 'nmax'): the M/D/1 queue. */
struct OpQsysMd1 {
    static Json run(Args& a) {
        const double lambda = a.number<double>("lambda");
        const double mu = a.number<double>("mu");
        const long nmax = a.has("nmax") ? static_cast<long>(a.count("nmax")) : -1L;
        a.done();
        return encode_qsys_res(qsys::qsys_md1(lambda, mu, nmax));
    }
};
 
/** The complex roots and coefficients, as [re, im] pairs. */
Json encode_cplx_pairs(const std::vector<std::complex<double> >& v) {
    Json out = Json::array();
    for (const std::complex<double>& z : v) out.push_back(Json::array({z.real(), z.imag()}));
    return out;
}
 
/** qsys_mer1(lambda, mu, r, 'tPoints', 'percentiles', 'nmax'): the M/Er/1 queue. */
struct OpQsysMer1 {
    static Json run(Args& a) {
        const double lambda = a.number<double>("lambda");
        const double mu = a.number<double>("mu");
        const unsigned r = a.uinteger("r");
        const qsys::QsysResOptions o = read_qsys_res_options(a);
        a.done();
        const qsys::Mer1Result m = qsys::qsys_mer1(lambda, mu, r, o);
        Json out = encode_qsys_res(m.res);
        // MATLAB's res.roots and res.coeffs are complex; a complex number is a pair [re, im].
        out["roots"] = encode_cplx_pairs(m.roots);
        out["coeffs"] = encode_cplx_pairs(m.coeffs);
        return out;
    }
};
 
/** qsys_mh21(lambda, p1, mu1, mu2, 'tPoints', 'percentiles', 'nmax'): the M/H2/1 queue. */
struct OpQsysMh21 {
    static Json run(Args& a) {
        const double lambda = a.number<double>("lambda");
        const double p1 = a.number<double>("p1");
        const double mu1 = a.number<double>("mu1");
        const double mu2 = a.number<double>("mu2");
        const qsys::QsysResOptions o = read_qsys_res_options(a);
        a.done();
        return encode_qsys_res(qsys::qsys_mh21(lambda, p1, mu1, mu2, o));
    }
};
 
}  // namespace
 
void register_api_qsys(DispatchTable& t) {
    // moment: all three arithmetics
    t["moment_cumulant_from_raw"] = make<OpMomentCumulantFromRaw>();
    t["moment_raw_from_cumulant"] = make<OpMomentRawFromCumulant>();
    t["moment_factcumulant_from_factorial"] = make<OpMomentFactcumulantFromFactorial>();
    t["moment_factorial_from_factcumulant"] = make<OpMomentFactorialFromFactcumulant>();
    t["moment_binomial_from_tail"] = make<OpMomentBinomialFromTail>();
    t["moment_tail_from_binomial"] = make<OpMomentTailFromBinomial>();
    t["moment_housematrix"] = make<OpMomentHousematrix>();
    t["moment_tensorsize"] = make<OpMomentTensorsize>();
    t["moment_tensortrans"] = make<OpMomentTensortrans>();
    t["moment_jointtrans"] = make<OpMomentJointtrans>();
    t["moment_joint_aggregate"] = make<OpMjAggregate>();
    t["moment_joint_marking"] = make<OpMjMarking>();
    t["moment_joint_central_from_raw_mean"] = make<OpMjCentralFromRawMean>();
    t["moment_joint_raw_from_central"] = make<OpMjRawFromCentral>();
    t["moment_joint_binomial_from_factorial"] = make<OpMjBinomialFromFactorial>();
    t["moment_joint_binomial_from_negbinomial"] = make<OpMjBinomialFromNegbinomial>();
    t["moment_joint_binomial_from_tail"] = make<OpMjBinomialFromTail>();
    t["moment_joint_central_from_raw"] = make<OpMjCentralFromRaw>();
    t["moment_joint_central_from_tail"] = make<OpMjCentralFromTail>();
    t["moment_joint_cumulant_from_raw"] = make<OpMjCumulantFromRaw>();
    t["moment_joint_factcumulant_from_factorial"] = make<OpMjFactcumulantFromFactorial>();
    t["moment_joint_factorial_from_binomial"] = make<OpMjFactorialFromBinomial>();
    t["moment_joint_factorial_from_factcumulant"] = make<OpMjFactorialFromFactcumulant>();
    t["moment_joint_factorial_from_raw"] = make<OpMjFactorialFromRaw>();
    t["moment_joint_factorial_from_upfactorial"] = make<OpMjFactorialFromUpfactorial>();
    t["moment_joint_negbinomial_from_binomial"] = make<OpMjNegbinomialFromBinomial>();
    t["moment_joint_negbinomial_from_upfactorial"] = make<OpMjNegbinomialFromUpfactorial>();
    t["moment_joint_raw_from_cumulant"] = make<OpMjRawFromCumulant>();
    t["moment_joint_raw_from_factorial"] = make<OpMjRawFromFactorial>();
    t["moment_joint_raw_from_upfactorial"] = make<OpMjRawFromUpfactorial>();
    t["moment_joint_tail_from_binomial"] = make<OpMjTailFromBinomial>();
    t["moment_joint_upfactorial_from_factorial"] = make<OpMjUpfactorialFromFactorial>();
    t["moment_joint_upfactorial_from_negbinomial"] = make<OpMjUpfactorialFromNegbinomial>();
    t["moment_joint_upfactorial_from_raw"] = make<OpMjUpfactorialFromRaw>();
 
    // fj: all three arithmetics
    t["fj_amva"] = make<OpFjAmva>();
    t["fj_char_max_blom"] = make<OpFjCharMaxBlom>();
    t["fj_cox_fit"] = make<OpFjCoxFit>();
    t["fj_dag_makespan"] = make<OpFjDagMakespan>();
    t["fj_ism_green"] = make<OpFjIsmGreen>();
    t["fj_lst_max_het"] = make<OpFjLstMaxHet>();
    t["fj_qgb"] = make<OpFjQgb>();
    t["fj_respt_bulk"] = make<OpFjResptBulk>();
    t["fj_respt_closed"] = make<OpFjResptClosed>();
    t["fj_respt_nosplit"] = make<OpFjResptNosplit>();
    t["fj_serialization"] = make<OpFjSerialization>();
    t["fj_tsm_capacity"] = make<OpFjTsmCapacity>();
    t["fj_xmax_coxian"] = make<OpFjXmaxCoxian>();
    t["fj_xmax_het"] = make<OpFjXmaxHet>();
    t["fj_xmax_hz"] = make<OpFjXmaxHz>();
    t["fj_xmax_moments_het"] = make<OpFjXmaxMomentsHet>();
    t["fj_mg1_respt_moments"] = make<OpFjMg1ResptMoments>();
    t["fj_branch_members"] = make<OpFjBranchMembers>();
    t["fj_dist2fj"] = make<OpFjDist2fj>();
    // fj: double and real
    t["fj_char_max_discrete"] = make<NeedsTranscendental<OpFjCharMaxDiscrete> >();
    t["fj_delay_opt"] = make<NeedsTranscendental<OpFjDelayOpt> >();
    t["fj_dispersion"] = make<NeedsTranscendental<OpFjDispersion> >();
    t["fj_tail_forktail"] = make<NeedsTranscendental<OpFjTailForktail> >();
    t["fj_tail_ordstat"] = make<NeedsTranscendental<OpFjTailOrdstat> >();
 
    // aoi
    t["aoi_lst_det"] = make<NeedsTranscendental<OpAoiLstDet> >();
    t["aoi_lst_erlang"] = make<OpAoiLstErlang>();
    t["aoi_lst_exp"] = make<OpAoiLstExp>();
    t["aoi_lst_ph"] = make<OpAoiLstPh>();
    t["aoi_dist2ph"] = make<OpAoiDist2ph>();
    // qsys: the registry lists all three arithmetics for these six, but the hh1,
    // lindley_moment, mm1 and mm1_tandem ports static_assert a transcendental
    // scalar, so an exact call is refused here with the Arithmetic reason.
    t["qsys_hh1_lindley"] = make<NeedsTranscendental<OpQsysHh1Lindley>>();
    t["qsys_lindley_moment"] = make<NeedsTranscendental<OpQsysLindleyMoment>>();
    t["qsys_mm1_lindley"] = make<NeedsTranscendental<OpQsysMm1Lindley>>();
    t["qsys_mm1_ps"] = make<OpQsysMm1Ps>();
    t["qsys_mm1_tandem_lindley"] = make<NeedsTranscendental<OpQsysMm1TandemLindley>>();
    t["qsys_tandem_lindley"] = make<OpQsysTandemLindley>();
    // qsys: double and real
    t["qsys_mapg1k"] = make<NeedsTranscendental<OpQsysMapg1k> >();
    t["qsys_mmapg1k"] = make<NeedsTranscendental<OpQsysMmapg1k> >();
    t["qsys_mapg1k_perflow"] = make<NeedsTranscendental<OpQsysMapg1kPerflow> >();
    t["qsys_svclaw_coeffs"] = make<NeedsTranscendental<OpQsysSvclawCoeffs> >();
    t["qsys_ldps_workload"] = make<NeedsTranscendental<OpQsysLdpsWorkload> >();
    t["qsys_mmapgk1"] = make<NeedsTranscendental<OpQsysMmapgk1> >();
    t["qsys_mmapgk1_lcfspr"] = make<NeedsTranscendental<OpQsysMmapgk1Lcfspr> >();
    t["qsys_tandem_ub_ciucu"] = make<NeedsTranscendental<OpQsysTandemUbCiucu> >();
    t["qsys_bmapm1"] = make<NeedsTranscendental<OpQsysBmapm1> >();
    t["qsys_maxima_twomoment"] = make<NeedsTranscendental<OpQsysMaximaTwomoment> >();
    t["qsys_mgisrgi_whitt"] = make<NeedsTranscendental<OpQsysMgisrgiWhitt> >();
    t["qsys_mapmapc"] = make<NeedsTranscendental<OpQsysMapmapc> >();
    t["qsys_mmapmap1"] = make<NeedsTranscendental<OpQsysMmapmap1> >();
    // qsys: double only
    t["qsys_mg1_ps"] = make<DoubleOnly<OpQsysMg1Ps> >();
    t["qsys_mmckk"] = make<DoubleOnly<OpQsysMmckk> >();
    t["qsys_engset"] = make<DoubleOnly<OpQsysEngset> >();
    t["qsys_mg1kk"] = make<DoubleOnly<OpQsysMg1kk> >();
    t["qsys_mmcc"] = make<DoubleOnly<OpQsysMmcc> >();
    t["qsys_mgc_ps"] = make<DoubleOnly<OpQsysMgcPs> >();
    t["qsys_mg1_sjf"] = make<DoubleOnly<OpQsysMg1Sjf> >();
    t["qsys_mg1_busy"] = make<DoubleOnly<OpQsysMg1Busy> >();
    t["qsys_gimc"] = make<DoubleOnly<OpQsysGimc> >();
    t["qsys_md1"] = make<DoubleOnly<OpQsysMd1> >();
    t["qsys_mer1"] = make<DoubleOnly<OpQsysMer1> >();
    t["qsys_mh21"] = make<DoubleOnly<OpQsysMh21> >();
    // qsys exact closed forms: all three arithmetics
    t["qsys_ddckk"] = make<OpQsysDdckk>();
    t["qsys_mg1_setup"] = make<OpQsysMg1Setup>();
    t["qsys_mg1_exceptional"] = make<OpQsysMg1Exceptional>();
    t["qsys_mg1_prio_pr"] = make<OpQsysMg1PrioPr>();
    t["qsys_mmc_prio"] = make<OpQsysMmcPrio>();
    t["qsys_mmc_prio_tardiness"] = make<DoubleOnly<OpQsysMmcPrioTardiness> >();
    t["qsys_mmc_prio_eta"] = make<DoubleOnly<OpQsysMmcPrioEta> >();
    t["qsys_mg1_conservation"] = make<OpQsysMg1Conservation>();
 
    // sim
    t["sim_quest_options"] = make<DoubleOnly<OpSimQuestOptions> >();
    t["sim_quest_heuristic_ci"] = make<OpSimQuestHeuristicCi>();
    t["sim_runlength"] = make<OpSimRunlength>();
}
 
void register_library_only_qsys(LibraryOnlyTable& t) {
    t["qsys_gig1_rq"] =
        "its interarrival law IaFun is a function handle (a CDF or transform evaluated by the "
        "port), and the port offers no constructor that builds it from a named law.";
    t["qsys_mg1k_loss"] =
        "its service law svc_density is a density function handle the port integrates, and the "
        "port takes no named-law form for it.";
    t["qsys_ggingi_tga"] =
        "its patience law (patienceCcdf, patiencePdf) and service law (serviceCcdf) are function "
        "handles the port integrates, with no named-law form.";
    t["qsys_ggisgi_fluid"] =
        "its patience and service laws (patienceCcdf, servingCcdf) are function handles the port "
        "integrates, with no named-law form.";
    t["qsys_gtmtst_fluid"] =
        "its time-varying arrival rate, staffing and service rate (lambdaFun, sFun, muFun) and "
        "its patience law (patienceCcdf) are function handles of time.";
    t["qsys_mtginf"] =
        "its time-varying arrival rate lambdaFun is a function handle of time and its service law "
        "serviceCcdf a function handle, with no named-law form.";
    t["qsys_mtgs0_mol"] =
        "its time-varying arrival rate lambdaFun is a function handle of time and its service law "
        "serviceCcdf a function handle, with no named-law form.";
    t["fj_order_stat"] =
        "its branch law F_X is a CDF function handle the port evaluates, with no named-law form.";
    t["fj_xmax_hz_het"] =
        "its branch laws are a list of CDF function handles the port integrates, with no "
        "named-law form.";
}
 
}  // namespace reg
}  // namespace line