api_dispatch_mam.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 'mam' group (see api_dispatch_detail.h).
 *
 * Conventions shared by every entry below:
 *  - argument keys are the MATLAB parameter names of the reference .m file
 *    (matlab/lib/kpctoolbox/, matlab/lib/m3a/, matlab/src/api/{mam,me,mc}/),
 *    output keys its output names. A function with no MATLAB reference uses
 *    the parameter names of its C++ port;
 *  - a MAP is the list [D0, D1] and an MMAP the list [D0, D1, D11, ..., D1K],
 *    exactly the MATLAB cell arrays; a list of MAPs is a list of such lists;
 *  - an ABSENT optional argument selects the port's own shorter overload, or
 *    leaves a field of an options struct at the port's value. No default is
 *    restated here. Where the port only has a long form that takes several
 *    optionals together, giving one of them requires giving the others, and
 *    the refusal names the missing one;
 *  - indices that are 1-based in MATLAB (map_stochcomp retain_idx,
 *    smp_passage_moments target, me_* refstat) are 1-based here as well;
 *  - an infinite server count or buffer is the string "Inf", as MATLAB's Inf.
 */
 
#include <cmath>
#include <cstdint>
#include <initializer_list>
#include <random>
#include <string>
#include <utility>
#include <vector>
 
#include "api_dispatch_detail.h"
 
#include "line/api/mam/amap2_adjust_gamma.h"
#include "line/api/mam/amap2_assemble.h"
#include "line/api/mam/amap2_fit_gamma.h"
#include "line/api/mam/amap2_fitall_gamma.h"
#include "line/api/mam/aph2_adjust.h"
#include "line/api/mam/aph2_assemble.h"
#include "line/api/mam/aph2_fit.h"
#include "line/api/mam/aph2_fitall.h"
#include "line/api/mam/ldqbd.h"
#include "line/api/mam/ldqbd_censor.h"
#include "line/api/mam/ldqbd_mphc.h"
#include "line/api/mam/m3pp2m_fitc.h"
#include "line/api/mam/m3pp2m_fitc_approx.h"
#include "line/api/mam/m3pp2m_fitc_trace.h"
#include "line/api/mam/m3pp_superpos_fitc.h"
#include "line/api/mam/map2_fit.h"
#include "line/api/mam/map2_fit_idc.h"
#include "line/api/mam/map_acfc.h"
#include "line/api/mam/map_cdf.h"
#include "line/api/mam/map_count_mean.h"
#include "line/api/mam/map_count_moment.h"
#include "line/api/mam/map_count_var.h"
#include "line/api/mam/map_mark.h"
#include "line/api/mam/map_moment.h"
#include "line/api/mam/map_pdf.h"
#include "line/api/mam/map_rand.h"
#include "line/api/mam/map_transform.h"
#include "line/api/mam/map_varcount.h"
#include "line/api/mam/mfq_fluflu_sojourn.h"
#include "line/api/mam/mfq_ld_distr.h"
#include "line/api/mam/mfq_ld_mean.h"
#include "line/api/mam/mfq_ld_solve.h"
#include "line/api/mam/mfq_multiregime.h"
#include "line/api/mam/mfq_prio_queue.h"
#include "line/api/mam/mfq_sojourn.h"
#include "line/api/mam/mmap_compress.h"
#include "line/api/mam/mmap_count_var.h"
#include "line/api/mam/mmap_lambda.h"
#include "line/api/mam/mmpp2_fit.h"
#include "line/api/mam/mmpp2_fit1.h"
#include "line/api/mam/mmpp2_fit2.h"
#include "line/api/mam/mmpp2_fit3.h"
#include "line/api/mam/mmpp2_fit4.h"
#include "line/api/mam/mmpp2_fitc.h"
#include "line/api/mam/mmpp2_fitc_approx.h"
#include "line/api/mam/qbd_bmapbmap1.h"
#include "line/api/mam/qbd_mapmap1.h"
#include "line/api/mam/qbd_r.h"
#include "line/api/mam/qbd_rap.h"
#include "line/api/mam/qbd_setupdelayoff.h"
#include "line/api/mam/qbd_setupdelayoff_mmk.h"
#include "line/api/mc/ctmc_passage.h"
#include "line/api/me/me_cqn.h"
#include "line/api/me/me_gegec_mql.h"
#include "line/api/me/me_gegecn.h"
#include "line/api/me/me_gegecn_pb.h"
#include "line/api/me/me_mqn.h"
#include "line/api/me/me_oqn.h"
#include "line/api/me/me_oqn_blk.h"
 
namespace line {
namespace reg {
 
namespace {
 
// ---------------------------------------------------------------------------
// Shared helpers
// ---------------------------------------------------------------------------
 
std::string quoted(const char* key) { return std::string("'") + key + "'"; }
 
/**
 * The long form of a port whose short overload fixes some arguments itself:
 * selected as soon as any of `triggers` is given, and then every one of
 * `needed` must be given too, because the long form has no default for it.
 */
bool long_form(const Args& a, const char* fn, std::initializer_list<const char*> triggers,
               std::initializer_list<const char*> needed) {
    std::string trig;
    for (const char* t : triggers)
        if (a.has(t)) {
            trig = t;
            break;
        }
    if (trig.empty()) return false;
    for (const char* n : needed)
        if (!a.has(n))
            throw InputError(std::string(fn) + ": '" + trig + "' selects the long form of the C++ "
                             "port, which has no default for '" + n + "'; give '" + n +
                             "' as well, or drop '" + trig + "' to use the port's own defaults.");
    return true;
}
 
template <class T>
mam::Map<T> map_from_json(const Json& j, const std::string& where) {
    if (!j.is_array() || j.size() != 2)
        throw InputError(where + ": a MAP is a list [D0, D1] of two square 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]");
    if (m.D0.rows() != m.D0.cols() || m.D1.rows() != m.D0.rows() || m.D1.cols() != m.D0.cols())
        throw InputError(where + ": D0 and D1 must be square matrices of the same order");
    return m;
}
 
template <class T>
mam::Mmap<T> mmap_from_json(const Json& j, const std::string& where) {
    if (!j.is_array() || j.size() < 3)
        throw InputError(where + ": an MMAP is a list [D0, D1, D11, ..., D1K] of at least three "
                         "square matrices");
    mam::Mmap<T> m;
    m.D0 = matrix_from_json<T>(j[0], where + "[0]");
    m.D1 = matrix_from_json<T>(j[1], where + "[1]");
    for (std::size_t k = 2; k < j.size(); ++k)
        m.Dc.push_back(matrix_from_json<T>(j[k], where + "[" + std::to_string(k) + "]"));
    const std::size_t n = m.D0.rows();
    bool ok = m.D0.cols() == n && m.D1.rows() == n && m.D1.cols() == n;
    for (const Matrix<T>& D : m.Dc) ok = ok && D.rows() == n && D.cols() == n;
    if (!ok) throw InputError(where + ": every matrix of an MMAP must be square of the same order");
    return m;
}
 
template <class T>
mam::Map<T> read_map(Args& a, const char* key) {
    return map_from_json<T>(a.get(key), quoted(key));
}
 
template <class T>
mam::Mmap<T> read_mmap(Args& a, const char* key) {
    return mmap_from_json<T>(a.get(key), quoted(key));
}
 
template <class T>
std::vector<mam::Map<T> > read_maps(Args& a, const char* key) {
    const Json& j = a.get(key);
    if (!j.is_array() || j.empty())
        throw InputError(quoted(key) + ": expected a non-empty list of MAPs, each [D0, D1]");
    std::vector<mam::Map<T> > out;
    for (std::size_t i = 0; i < j.size(); ++i)
        out.push_back(map_from_json<T>(j[i], quoted(key) + "[" + std::to_string(i) + "]"));
    return out;
}
 
template <class T>
Json encode_map(const mam::Map<T>& m) {
    Json o = Json::array();
    o.push_back(encode_matrix(m.D0));
    o.push_back(encode_matrix(m.D1));
    return o;
}
 
template <class T>
Json encode_mmap(const mam::Mmap<T>& m) {
    Json o = Json::array();
    o.push_back(encode_matrix(m.D0));
    o.push_back(encode_matrix(m.D1));
    for (const Matrix<T>& D : m.Dc) o.push_back(encode_matrix(D));
    return o;
}
 
template <class T>
Json encode_maps(const std::vector<mam::Map<T> >& v) {
    Json o = Json::array();
    for (const mam::Map<T>& m : v) o.push_back(encode_map(m));
    return o;
}
 
template <class T>
Json encode_mmaps(const std::vector<mam::Mmap<T> >& v) {
    Json o = Json::array();
    for (const mam::Mmap<T>& m : v) o.push_back(encode_mmap(m));
    return o;
}
 
/** Non-negative integers, for lags, orders and powers. */
std::vector<unsigned> unsigneds(Args& a, const char* key) {
    std::vector<unsigned> out;
    for (int v : a.ints(key)) {
        if (v < 0)
            throw InputError(quoted(key) + ": entries must be non-negative integers, got " +
                             std::to_string(v));
        out.push_back(static_cast<unsigned>(v));
    }
    return out;
}
 
/** MATLAB 1-based indices, returned 0-based. */
std::vector<std::size_t> one_based(const Json& j, const std::string& where) {
    std::vector<std::size_t> out;
    for (int v : int_vector_from_json(j, where)) {
        if (v < 1)
            throw InputError(where + ": indices are 1-based, as in MATLAB; got " + std::to_string(v));
        out.push_back(static_cast<std::size_t>(v - 1));
    }
    return out;
}
 
/** A list of numeric vectors. */
template <class T>
std::vector<std::vector<T> > vectors_from_json(const Json& j, const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected a list of vectors");
    std::vector<std::vector<T> > out;
    for (std::size_t i = 0; i < j.size(); ++i)
        out.push_back(vector_from_json<T>(j[i], where + "[" + std::to_string(i) + "]"));
    return out;
}
 
template <class T>
std::vector<std::vector<T> > read_vectors(Args& a, const char* key) {
    return vectors_from_json<T>(a.get(key), quoted(key));
}
 
/**
 * A count per station in which MATLAB's Inf marks "unbounded" (an infinite
 * server, an infinite buffer); api/me spells that 0. A finite entry must be a
 * positive integer, so a literal 0 cannot be mistaken for either meaning.
 */
std::vector<long> counts_with_inf(Args& a, const char* key) {
    std::vector<long> out;
    for (double d : a.vector<double>(key)) {
        if (std::isinf(d) && d > 0) {
            out.push_back(0);
            continue;
        }
        if (d != std::floor(d) || d < 1)
            throw InputError(quoted(key) + ": entries must be positive integers or \"Inf\", got " +
                             shortest_decimal(d));
        out.push_back(static_cast<long>(d));
    }
    return out;
}
 
/** Class populations; with `open_inf` an Inf entry (an open class) reads as 0. */
std::vector<long> populations(Args& a, const char* key, bool open_inf) {
    std::vector<long> out;
    for (double d : a.vector<double>(key)) {
        if (open_inf && std::isinf(d) && d > 0) {
            out.push_back(0);
            continue;
        }
        if (d != std::floor(d) || d < 0)
            throw InputError(quoted(key) + ": a population must be a non-negative integer, got " +
                             shortest_decimal(d));
        out.push_back(static_cast<long>(d));
    }
    return out;
}
 
/** A per-entry flag vector (true/false or 0/1), as std::vector<char>. */
std::vector<char> flags(Args& a, const char* key) {
    const Json& j = a.get(key);
    std::vector<char> out;
    std::vector<Json> flat;
    if (j.is_array())
        for (const Json& e : j) flat.push_back(e);
    else
        flat.push_back(j);
    for (const Json& e : flat) {
        if (e.is_boolean()) {
            out.push_back(e.get<bool>() ? 1 : 0);
            continue;
        }
        const std::vector<int> v = int_vector_from_json(e, quoted(key));
        if (v.size() != 1 || (v[0] != 0 && v[0] != 1))
            throw InputError(quoted(key) + ": entries must be true, false, 0 or 1");
        out.push_back(static_cast<char>(v[0]));
    }
    return out;
}
 
/** Per-class routing: R matrices (M x M), or a single matrix when R == 1. */
template <class T>
std::vector<Matrix<T> > per_class_routing(Args& a, const char* key, std::size_t R) {
    const Json& j = a.get(key);
    const bool single = j.is_array() && !j.empty() && j[0].is_array() && !j[0].empty() &&
                        !j[0][0].is_array();
    if (single) {
        if (R != 1)
            throw InputError(quoted(key) + ": with R = " + std::to_string(R) +
                             " classes, give a list of R routing matrices (one M x M matrix per "
                             "class, P{r}(i,j) = MATLAB P(i,j,r))");
        return std::vector<Matrix<T> >(1, matrix_from_json<T>(j, quoted(key)));
    }
    return a.matrices<T>(key);
}
 
/** MATLAB refstat: 1-based per class, 0 for the default; api/me spells the default -1. */
std::vector<long> reference_stations(Args& a, const char* key, std::size_t R) {
    if (!a.has(key)) {
        a.opt(key);
        return std::vector<long>(R, -1);
    }
    std::vector<long> out;
    for (int v : a.ints(key)) {
        if (v < 0) throw InputError(quoted(key) + ": station indices are 1-based, 0 for the default");
        out.push_back(v == 0 ? -1L : static_cast<long>(v - 1));
    }
    return out;
}
 
/** The MATLAB options struct of me_cqn/me_mqn/me_oqn. */
me::MeOptions me_options(Args& a) {
    me::MeOptions o;
    const Json* j = a.opt("options");
    if (!j) return o;
    if (!j->is_object()) throw InputError("'options' must be an object {tol, maxiter, verbose}");
    for (auto it = j->begin(); it != j->end(); ++it) {
        const std::string k = it.key();
        if (k == "tol") {
            o.tol = number_from_json<double>(*it, "'options.tol'");
        } else if (k == "maxiter") {
            const std::vector<int> v = int_vector_from_json(*it, "'options.maxiter'");
            if (v.size() != 1 || v[0] < 1)
                throw InputError("'options.maxiter' must be a positive integer");
            o.maxiter = v[0];
        } else if (k == "verbose") {
            const bool on = it->is_boolean() ? it->get<bool>()
                                             : number_from_json<double>(*it, "'options.verbose'") != 0;
            if (on)
                throw UnsupportedError(ErrorKind{"JsonApi"},
                                       "'options.verbose' = true prints iteration traces to the "
                                       "console, which the JSON API has no channel for; omit it.");
        } else {
            throw InputError("'options': unknown field '" + k + "'; accepted: tol, maxiter, verbose");
        }
    }
    return o;
}
 
me::MeBlkOptions me_blk_options(Args& a) {
    me::MeBlkOptions o;
    const Json* j = a.opt("options");
    if (!j) return o;
    if (!j->is_object())
        throw InputError("'options' must be an object {tol, maxiter, verbose, damping}");
    for (auto it = j->begin(); it != j->end(); ++it) {
        const std::string k = it.key();
        if (k == "tol") {
            o.tol = number_from_json<double>(*it, "'options.tol'");
        } else if (k == "damping") {
            o.damping = number_from_json<double>(*it, "'options.damping'");
        } else if (k == "maxiter") {
            const std::vector<int> v = int_vector_from_json(*it, "'options.maxiter'");
            if (v.size() != 1 || v[0] < 1)
                throw InputError("'options.maxiter' must be a positive integer");
            o.maxiter = v[0];
        } else if (k == "verbose") {
            const bool on = it->is_boolean() ? it->get<bool>()
                                             : number_from_json<double>(*it, "'options.verbose'") != 0;
            if (on)
                throw UnsupportedError(ErrorKind{"JsonApi"},
                                       "'options.verbose' = true prints iteration traces to the "
                                       "console, which the JSON API has no channel for; omit it.");
        } else {
            throw InputError("'options': unknown field '" + k +
                             "'; accepted: tol, maxiter, verbose, damping");
        }
    }
    return o;
}
 
template <class T>
Json encode_me_result(const me::MeResult<T>& r, bool with_X) {
    Json o;
    o["L"] = encode_matrix(r.L);
    o["W"] = encode_matrix(r.W);
    o["Ca"] = encode_matrix(r.Ca);
    o["Cd"] = encode_matrix(r.Cd);
    o["lambda"] = encode_matrix(r.lambda);
    o["rho"] = encode_matrix(r.rho);
    if (with_X) o["X"] = encode_vector(r.X);
    o["iter"] = static_cast<std::int64_t>(r.iter);
    o["converged"] = r.converged;
    return o;
}
 
/** The building blocks returned by mfq_ld_solve, as mfq_ld_mean/distr take them. */
template <class T>
mam::LevelDependentFluidBlocks<T> read_ld_blocks(Args& a) {
    mam::LevelDependentFluidBlocks<T> b;
    b.masses = read_vectors<T>(a, "masses");
    b.iniF = read_vectors<T>(a, "iniF");
    b.KF = a.matrices<T>("KF");
    b.cloF = a.matrices<T>("cloF");
    b.iniB = read_vectors<T>(a, "iniB");
    b.KB = a.matrices<T>("KB");
    b.cloB = a.matrices<T>("cloB");
    b.Thr = a.vector<T>("T");
    return b;
}
 
/** A rate vector, given as MATLAB's diagonal matrix or as the vector itself. */
std::vector<double> diagonal_rates(const Json& j, const std::string& where) {
    const Matrix<double> m = matrix_from_json<double>(j, where);
    if (m.empty()) return std::vector<double>();
    if (m.rows() == 1 || m.cols() == 1) return vector_from_json<double>(j, where);
    if (m.rows() != m.cols())
        throw InputError(where + ": a drift-rate matrix must be square and diagonal");
    std::vector<double> d(m.rows());
    for (std::size_t i = 0; i < m.rows(); ++i)
        for (std::size_t k = 0; k < m.cols(); ++k) {
            if (i == k) {
                d[i] = m(i, i);
            } else if (m(i, k) != 0.0) {
                throw InputError(where + ": a drift-rate matrix must be diagonal (entry (" +
                                 std::to_string(i + 1) + "," + std::to_string(k + 1) +
                                 ") is nonzero)");
            }
        }
    return d;
}
 
std::vector<std::vector<double> > diagonal_rate_list(Args& a, const char* key, bool optional) {
    const Json* j = optional ? a.opt(key) : &a.get(key);
    std::vector<std::vector<double> > out;
    if (!j) return out;
    if (!j->is_array()) throw InputError(quoted(key) + ": expected a list of drift-rate matrices");
    for (std::size_t i = 0; i < j->size(); ++i)
        out.push_back(diagonal_rates((*j)[i], quoted(key) + "[" + std::to_string(i) + "]"));
    return out;
}
 
std::vector<mam::FluidBoundary> boundaries(Args& a, const char* key) {
    std::vector<mam::FluidBoundary> out;
    for (int v : a.ints_or_empty(key)) {
        if (v != 0 && v != 1)
            throw InputError(quoted(key) + ": entries are 0 (reflective) or 1 (absorbing)");
        out.push_back(v == 0 ? mam::FluidBoundary::Reflective : mam::FluidBoundary::Absorbing);
    }
    return out;
}
 
mam::FluidDistrKind distr_kind(Args& a) {
    const std::string w = a.text("what", "");
    if (!a.has("what")) a.get("what");
    if (w == "pdf") return mam::FluidDistrKind::Pdf;
    if (w == "pdfd") return mam::FluidDistrKind::Pdfd;
    if (w == "cdf") return mam::FluidDistrKind::Cdf;
    if (w == "cdfm") return mam::FluidDistrKind::Cdfm;
    throw InputError("'what' must be one of \"pdf\", \"pdfd\", \"cdf\" or \"cdfm\", got \"" + w + "\"");
}
 
mam::MmapCompressMethod compress_method(const std::string& s) {
    if (s == "default" || s == "mixture" || s == "mixture.order1")
        return mam::MmapCompressMethod::MixtureOrder1;
    if (s == "mixture.order2") return mam::MmapCompressMethod::MixtureOrder2;
    if (s == "mamap2") return mam::MmapCompressMethod::Mamap2;
    if (s == "mamap2.fb") return mam::MmapCompressMethod::Mamap2Fb;
    if (s == "m3pp.approx_cov") return mam::MmapCompressMethod::M3ppApproxCov;
    if (s == "m3pp.approx_ag") return mam::MmapCompressMethod::M3ppApproxAg;
    if (s == "m3pp.exact_delta") return mam::MmapCompressMethod::M3ppExactDelta;
    if (s == "m3pp.approx_delta") return mam::MmapCompressMethod::M3ppApproxDelta;
    throw InputError("'config.method': unknown compression method '" + s +
                     "'; accepted: default, mixture, mixture.order1, mixture.order2, mamap2, "
                     "mamap2.fb, m3pp.approx_cov, m3pp.approx_ag, m3pp.exact_delta, "
                     "m3pp.approx_delta");
}
 
// ---------------------------------------------------------------------------
// MAP descriptors (kpctoolbox/map)
// ---------------------------------------------------------------------------
 
#define LINE_MAM_MAP_UNARY(OpName, fn, outkey, enc)        \
    struct OpName {                                         \
        template <class T>                                  \
        static Json run(Args& a) {                          \
            const mam::Map<T> m = read_map<T>(a, "MAP");    \
            a.done();                                       \
            Json o;                                         \
            o[outkey] = enc(mam::fn(m));                    \
            return o;                                       \
        }                                                   \
    };
 
LINE_MAM_MAP_UNARY(OpMapInfgen, map_infgen, "Q", encode_matrix)
LINE_MAM_MAP_UNARY(OpMapProb, map_prob, "PIQ", encode_vector)
LINE_MAM_MAP_UNARY(OpMapLambda, map_lambda, "LAMBDA", encode_scalar)
LINE_MAM_MAP_UNARY(OpMapPie, map_pie, "PIE", encode_vector)
LINE_MAM_MAP_UNARY(OpMapMean, map_mean, "MEAN", encode_scalar)
LINE_MAM_MAP_UNARY(OpMapEmbedded, map_embedded, "P", encode_matrix)
LINE_MAM_MAP_UNARY(OpMapVar, map_var, "VAR", encode_scalar)
LINE_MAM_MAP_UNARY(OpMapScv, map_scv, "SCV", encode_scalar)
LINE_MAM_MAP_UNARY(OpMapIdc, map_idc, "I", encode_scalar)
LINE_MAM_MAP_UNARY(OpMapKurt, map_kurt, "KURT", encode_scalar)
LINE_MAM_MAP_UNARY(OpMapSkew, map_skew, "SKEWNESS", encode_scalar)
LINE_MAM_MAP_UNARY(OpMapNormalize, map_normalize, "MAP", encode_map)
#undef LINE_MAM_MAP_UNARY
 
struct OpMapMoment {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAP");
        const std::vector<unsigned> orders = unsigneds(a, "ORDERS");
        a.done();
        std::vector<T> out;
        for (unsigned k : orders) out.push_back(mam::map_moment(m, k));
        Json o;
        o["MOMENTS"] = encode_vector(out);
        return o;
    }
};
 
struct OpMapAcf {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAP");
        const std::vector<unsigned> lags = unsigneds(a, "LAGS");
        a.done();
        Json o;
        o["ACFCOEFFS"] = encode_vector(mam::map_acf(m, lags));
        return o;
    }
};
 
struct OpMapScale {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAPIN");
        const T mean = a.number<T>("NEWMEAN");
        a.done();
        Json o;
        o["MAPOUT"] = encode_map(mam::map_scale(m, mean));
        return o;
    }
};
 
struct OpMapRenewal {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAPIN");
        a.done();
        Json o;
        o["MAPOUT"] = encode_map(mam::map_renewal(m));
        return o;
    }
};
 
struct OpMapExponentialMean {
    template <class T>
    static Json run(Args& a) {
        const T mean = a.number<T>("mean");
        a.done();
        Json o;
        o["MAP"] = encode_map(mam::map_exponential_mean(mean));
        return o;
    }
};
 
struct OpMapErlang {
    template <class T>
    static Json run(Args& a) {
        const T mean = a.number<T>("MEAN");
        const unsigned k = a.uinteger("k");
        a.done();
        Json o;
        o["MAP"] = encode_map(mam::map_erlang(mean, k));
        return o;
    }
};
 
struct OpMapSum {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAP");
        const unsigned n = a.uinteger("n");
        a.done();
        Json o;
        o["MAP"] = encode_map(mam::map_sum(m, n));
        return o;
    }
};
 
struct OpMapSumind {
    template <class T>
    static Json run(Args& a) {
        const std::vector<mam::Map<T> > maps = read_maps<T>(a, "MAPs");
        a.done();
        Json o;
        o["MAP"] = encode_map(mam::map_sumind(maps));
        return o;
    }
};
 
struct OpMapMixture {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> alpha = a.vector<T>("alpha");
        const std::vector<mam::Map<T> > maps = read_maps<T>(a, "MAPs");
        a.done();
        Json o;
        o["MAP"] = encode_map(mam::map_mixture(alpha, maps));
        return o;
    }
};
 
struct OpPh2map {
    template <class T>
    static Json run(Args& a) {
        const Json& j = a.get("ph");
        if (!j.is_object())
            throw InputError("'ph' must be an object {\"alpha\": entry vector, \"subgen\": "
                             "sub-generator}");
        mam::PhType<T> ph;
        bool has_alpha = false, has_sub = false;
        for (auto it = j.begin(); it != j.end(); ++it) {
            if (it.key() == "alpha") {
                ph.alpha = vector_from_json<T>(*it, "'ph.alpha'");
                has_alpha = true;
            } else if (it.key() == "subgen") {
                ph.subgen = matrix_from_json<T>(*it, "'ph.subgen'");
                has_sub = true;
            } else {
                throw InputError("'ph': unknown field '" + it.key() + "'; accepted: alpha, subgen");
            }
        }
        if (!has_alpha || !has_sub) throw InputError("'ph' needs both 'alpha' and 'subgen'");
        a.done();
        Json o;
        o["MAP"] = encode_map(mam::ph2map(ph));
        return o;
    }
};
 
struct OpMapStochcomp {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAP");
        const std::vector<std::size_t> retain = one_based(a.get("retain_idx"), "'retain_idx'");
        a.done();
        Json o;
        o["MAPnew"] = encode_map(mam::map_stochcomp(m, retain));
        return o;
    }
};
 
struct OpMapJoint {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAP");
        const std::vector<unsigned> lags = unsigneds(a, "a");
        const std::vector<unsigned> pw = unsigneds(a, "i");
        a.done();
        Json o;
        o["JM"] = encode_scalar(mam::map_joint(m, lags, pw));
        return o;
    }
};
 
struct OpMapIsfeasible {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAP");
        const bool tol = a.has("TOL");
        const T t = tol ? a.number<T>("TOL") : T();
        a.opt("TOL");
        a.done();
        Json o;
        o["ISFEAS"] = tol ? mam::map_isfeasible(m, t) : mam::map_isfeasible(m);
        return o;
    }
};
 
struct OpMapHyperexp {
    template <class T>
    static Json run(Args& a) {
        const T mean = a.number<T>("MEAN");
        const T scv = a.number<T>("SCV");
        const bool hp = a.has("p");
        const T p = hp ? a.number<T>("p") : T();
        a.opt("p");
        a.done();
        Json o;
        o["MAP"] = encode_map(hp ? mam::map_hyperexp(mean, scv, p) : mam::map_hyperexp(mean, scv));
        return o;
    }
};
 
struct OpMapMark {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAP");
        const std::vector<T> prob = a.vector<T>("prob");
        a.done();
        Json o;
        o["MMAP"] = encode_mmap(mam::map_mark(m, prob));
        return o;
    }
};
 
struct OpMapAcfc {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAP");
        const std::vector<unsigned> kset = unsigneds(a, "kset");
        const T u = a.number<T>("u");
        a.done();
        Json o;
        o["acfc"] = encode_vector(mam::map_acfc(m, kset, u));
        return o;
    }
};
 
#define LINE_MAM_MAP_POINTS(OpName, fn, mapkey, ptkey, outkey)  \
    struct OpName {                                             \
        template <class T>                                      \
        static Json run(Args& a) {                              \
            const mam::Map<T> m = read_map<T>(a, mapkey);       \
            const std::vector<T> pts = a.vector<T>(ptkey);      \
            a.done();                                           \
            Json o;                                             \
            o[outkey] = encode_vector(mam::fn(m, pts));         \
            return o;                                           \
        }                                                       \
    };
 
LINE_MAM_MAP_POINTS(OpMapCdf, map_cdf, "MAP", "POINTS", "CDFVALS")
LINE_MAM_MAP_POINTS(OpMapPdf, map_pdf, "MAP", "tset", "f")
LINE_MAM_MAP_POINTS(OpMapVarcount, map_varcount, "MAP", "tset", "varc")
LINE_MAM_MAP_POINTS(OpMapCountMean, map_count_mean, "map", "t", "m")
LINE_MAM_MAP_POINTS(OpMapCountVar, map_count_var, "map", "t", "v")
#undef LINE_MAM_MAP_POINTS
 
struct OpMapCountMoment {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "MAP");
        const T t = a.number<T>("t");
        const std::vector<unsigned> orders = unsigneds(a, "order");
        a.done();
        Json o;
        o["M"] = encode_vector(mam::map_count_moment(m, t, orders));
        return o;
    }
};
 
// ---------------------------------------------------------------------------
// MMAP algebra (m3a/mmap)
// ---------------------------------------------------------------------------
 
struct OpMmapLambda {
    template <class T>
    static Json run(Args& a) {
        const mam::Mmap<T> m = read_mmap<T>(a, "MMAP");
        a.done();
        Json o;
        o["lambda"] = encode_vector(mam::mmap_lambda(m));
        return o;
    }
};
 
struct OpMmapCountLambda {
    template <class T>
    static Json run(Args& a) {
        const mam::Mmap<T> m = read_mmap<T>(a, "mmap");
        a.done();
        Json o;
        o["lk"] = encode_vector(mam::mmap_count_lambda(m));
        return o;
    }
};
 
struct OpMmapPc {
    template <class T>
    static Json run(Args& a) {
        const mam::Mmap<T> m = read_mmap<T>(a, "MMAP");
        a.done();
        Json o;
        o["pc"] = encode_vector(mam::mmap_pc(m));
        return o;
    }
};
 
struct OpMmapIsfeasible {
    template <class T>
    static Json run(Args& a) {
        const mam::Mmap<T> m = read_mmap<T>(a, "MMAP");
        const bool tol = a.has("TOL");
        const T t = tol ? a.number<T>("TOL") : T();
        a.opt("TOL");
        a.done();
        Json o;
        o["TF"] = tol ? mam::mmap_isfeasible_tol(m, t) : mam::mmap_isfeasible(m);
        return o;
    }
};
 
struct OpMmapNormalize {
    template <class T>
    static Json run(Args& a) {
        const mam::Mmap<T> m = read_mmap<T>(a, "MMAP");
        a.done();
        Json o;
        o["MMAP"] = encode_mmap(mam::mmap_normalize(m));
        return o;
    }
};
 
struct OpMmapScale {
    template <class T>
    static Json run(Args& a) {
        const mam::Mmap<T> m = read_mmap<T>(a, "MMAP");
        const Json& jm = a.get("M");
        const std::vector<T> M = vector_from_json<T>(jm, "'M'");
        if (M.size() != 1)
            throw UnsupportedError(ErrorKind{"JsonApi"},
                                   "mmap_scale: the C++ port scales to a single aggregate mean; "
                                   "the per-class form (M a vector of class means) is not ported.");
        a.unsupported("maxIter", "the MATLAB reference declares it but never reads it, so it "
                                 "has no effect to reproduce.");
        a.done();
        Json o;
        o["SCALED"] = encode_mmap(mam::mmap_scale(m, M[0]));
        return o;
    }
};
 
struct OpMmapSuper {
    template <class T>
    static Json run(Args& a) {
        const mam::Mmap<T> x = read_mmap<T>(a, "MMAPa");
        const mam::Mmap<T> y = read_mmap<T>(a, "MMAPb");
        const std::string opt = a.text("opt", "default");
        if (a.has("opt") && opt != "default")
            throw UnsupportedError(ErrorKind{"JsonApi"},
                                   "mmap_super: only opt = \"default\" is ported to C++; got \"" +
                                       opt + "\".");
        a.done();
        Json o;
        o["SUP"] = encode_mmap(mam::mmap_super(x, y));
        return o;
    }
};
 
struct OpMmapCountVar {
    template <class T>
    static Json run(Args& a) {
        const mam::Mmap<T> m = read_mmap<T>(a, "mmap");
        const T t = a.number<T>("t");
        a.done();
        Json o;
        o["vk"] = encode_vector(mam::mmap_count_var(m, t));
        return o;
    }
};
 
struct OpMmapCompress {
    template <class T>
    static Json run(Args& a) {
        const mam::Mmap<T> m = read_mmap<T>(a, "MMAP");
        const Json* c = a.opt("config");
        bool has_method = false;
        mam::MmapCompressMethod meth = mam::MmapCompressMethod::MixtureOrder1;
        if (c) {
            if (!c->is_object()) throw InputError("'config' must be an object {method}");
            for (auto it = c->begin(); it != c->end(); ++it) {
                if (it.key() != "method")
                    throw InputError("'config': unknown field '" + it.key() + "'; accepted: method");
                if (!it->is_string()) throw InputError("'config.method' must be a string");
                meth = compress_method(it->get<std::string>());
                has_method = true;
            }
        }
        a.done();
        Json o;
        o["MMAP"] = encode_mmap(has_method ? mam::mmap_compress(m, meth) : mam::mmap_compress(m));
        return o;
    }
};
 
// ---------------------------------------------------------------------------
// QBD
// ---------------------------------------------------------------------------
 
/** qbd_R and qbd_R_logred: MATLAB iter_max, plus the port's tol alongside it. */
template <bool Logred>
struct OpQbdR {
    template <class T>
    static Json run(Args& a) {
        const char* fn = Logred ? "qbd_R_logred" : "qbd_R";
        const Matrix<T> B = a.matrix<T>("B");
        const Matrix<T> L = a.matrix<T>("L");
        const Matrix<T> F = a.matrix<T>("F");
        const bool lf = long_form(a, fn, {"iter_max", "tol"}, {"iter_max", "tol"});
        const unsigned it = lf ? a.uinteger("iter_max") : 0u;
        const T tol = lf ? a.number<T>("tol") : T();
        a.opt("iter_max");
        a.opt("tol");
        a.done();
        Json o;
        if (Logred)
            o["R"] = encode_matrix(lf ? mam::qbd_R_logred(B, L, F, it, tol) : mam::qbd_R_logred(B, L, F));
        else
            o["R"] = encode_matrix(lf ? mam::qbd_R(B, L, F, it, tol) : mam::qbd_R(B, L, F));
        return o;
    }
};
 
struct OpQbdFundmat {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> B = a.matrix<T>("B");
        const Matrix<T> L = a.matrix<T>("L");
        const Matrix<T> F = a.matrix<T>("F");
        const std::string which = a.text("matrices", "GR");
        if (which != "G" && which != "R" && which != "GR" && which != "RG")
            throw InputError("qbd_fundmat: 'matrices' must be \"G\", \"R\", \"GR\" or \"RG\"");
        const bool lf = long_form(a, "qbd_fundmat", {"precision", "maxNumIt"}, {"precision", "maxNumIt"});
        const T tol = lf ? a.number<T>("precision") : T();
        const unsigned it = lf ? a.uinteger("maxNumIt") : 0u;
        a.opt("precision");
        a.opt("maxNumIt");
        a.done();
        const mam::QbdFundMat<T> r = lf ? mam::qbd_fundmat(B, L, F, it, tol) : mam::qbd_fundmat(B, L, F);
        Json o;
        if (which.find('G') != std::string::npos) o["G"] = encode_matrix(r.G);
        if (which.find('R') != std::string::npos) o["R"] = encode_matrix(r.R);
        return o;
    }
};
 
struct OpQbdCaudal {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> R = a.matrix<T>("R");
        const bool lf = long_form(a, "qbd_caudal", {"iter_max", "tol"}, {"iter_max", "tol"});
        const unsigned it = lf ? a.uinteger("iter_max") : 0u;
        const T tol = lf ? a.number<T>("tol") : T();
        a.opt("iter_max");
        a.opt("tol");
        a.done();
        Json o;
        o["eta"] = encode_scalar(lf ? mam::qbd_caudal(R, it, tol) : mam::qbd_caudal(R));
        return o;
    }
};
 
struct OpQbdPi {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> B = a.matrix<T>("B");
        const Matrix<T> Lbar = a.matrix<T>("Lbar");
        const Matrix<T> R = a.matrix<T>("R");
        const bool lf = long_form(a, "qbd_pi", {"max_levels", "mass_tol"}, {"max_levels", "mass_tol"});
        const std::size_t ml = lf ? a.count("max_levels") : 0u;
        const T mt = lf ? a.number<T>("mass_tol") : T();
        a.opt("max_levels");
        a.opt("mass_tol");
        a.done();
        Json o;
        o["pi"] = encode_matrix(lf ? mam::qbd_pi(B, Lbar, R, ml, mt) : mam::qbd_pi(B, Lbar, R));
        return o;
    }
};
 
template <bool G>
struct OpQbdResidual {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> B = a.matrix<T>("B");
        const Matrix<T> L = a.matrix<T>("L");
        const Matrix<T> F = a.matrix<T>("F");
        const Matrix<T> X = a.matrix<T>(G ? "G" : "R");
        a.done();
        Json o;
        o["residual"] = encode_scalar(G ? mam::qbd_G_residual(B, L, F, X) : mam::qbd_R_residual(B, L, F, X));
        return o;
    }
};
 
struct OpQbdMapmap1Blocks {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> arr = read_map<T>(a, "arrival");
        const mam::Map<T> srv = read_map<T>(a, "service");
        a.done();
        const mam::QbdMapMap1Blocks<T> r = mam::qbd_mapmap1_blocks(arr, srv);
        Json o;
        o["B"] = encode_matrix(r.B);
        o["L"] = encode_matrix(r.L);
        o["F"] = encode_matrix(r.F);
        o["Lbar"] = encode_matrix(r.Lbar);
        return o;
    }
};
 
template <bool Fact>
struct OpQbdQlenMoment {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> pi0 = a.vector<T>("pi0");
        const Matrix<T> R = a.matrix<T>("R");
        const unsigned m = a.uinteger("m");
        a.done();
        Json o;
        if (Fact)
            o["factmoment"] = encode_scalar(mam::qbd_qlen_factmoment(pi0, R, m));
        else
            o["moment"] = encode_scalar(mam::qbd_qlen_moment(pi0, R, m));
        return o;
    }
};
 
struct OpQbdRg {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> arr = read_map<T>(a, "MAPa");
        const mam::Map<T> srv = read_map<T>(a, "MAPs");
        const bool hu = a.has("util");
        const T util = hu ? a.number<T>("util") : T();
        a.opt("util");
        a.done();
        const mam::QbdRg<T> r = hu ? mam::qbd_rg(arr, srv, util) : mam::qbd_rg(arr, srv);
        Json o;
        o["R"] = encode_matrix(r.R);
        o["G"] = encode_matrix(r.G);
        o["B"] = encode_matrix(r.B);
        o["L"] = encode_matrix(r.L);
        o["F"] = encode_matrix(r.F);
        o["U"] = encode_matrix(r.U);
        return o;
    }
};
 
struct OpQbdMapmap1 {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> arr = read_map<T>(a, "MAPa");
        const mam::Map<T> srv = read_map<T>(a, "MAPs");
        const bool lf = long_form(a, "qbd_mapmap1", {"util", "max_levels"}, {"util", "max_levels"});
        const T util = lf ? a.number<T>("util") : T();
        const std::size_t ml = lf ? a.count("max_levels") : 0u;
        a.opt("util");
        a.opt("max_levels");
        a.done();
        const mam::QbdMapMap1Result<T> r = lf ? mam::qbd_mapmap1(arr, srv, util, ml) : mam::qbd_mapmap1(arr, srv);
        Json o;
        o["XN"] = encode_scalar(r.XN);
        o["QN"] = encode_scalar(r.QN);
        o["UN"] = encode_scalar(r.UN);
        o["pqueue"] = encode_matrix(r.pqueue);
        o["R"] = encode_matrix(r.R);
        o["eta"] = encode_scalar(r.eta);
        o["G"] = encode_matrix(r.G);
        o["A_1"] = encode_matrix(r.B);
        o["A0"] = encode_matrix(r.L);
        o["A1"] = encode_matrix(r.F);
        o["U"] = encode_matrix(r.U);
        o["MAPs"] = encode_map(r.service);
        return o;
    }
};
 
struct OpQbdMapmap1QlenTruncated {
    template <class T>
    static Json run(Args& a) {
        const Json& j = a.get("res");
        if (!j.is_object() || j.find("pqueue") == j.end())
            throw InputError("'res' must be the result object of qbd_mapmap1, carrying at least "
                             "'pqueue' (the only field this function reads)");
        mam::QbdMapMap1Result<T> r;
        r.pqueue = matrix_from_json<T>(j["pqueue"], "'res.pqueue'");
        a.done();
        Json o;
        o["QN"] = encode_scalar(mam::qbd_mapmap1_qlen_truncated(r));
        return o;
    }
};
 
struct OpQbdRap {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> A0 = a.matrix<T>("A0");
        const Matrix<T> A1 = a.matrix<T>("A1");
        const Matrix<T> A2 = a.matrix<T>("A2");
        const bool lf = long_form(a, "qbd_rap", {"B0", "B1", "numLevels"}, {"B0", "B1", "numLevels"});
        const Matrix<T> B0 = lf ? a.matrix<T>("B0") : Matrix<T>();
        const Matrix<T> B1 = lf ? a.matrix<T>("B1") : Matrix<T>();
        const std::size_t nl = lf ? a.count("numLevels") : 0u;
        a.opt("B0");
        a.opt("B1");
        a.opt("numLevels");
        a.done();
        const mam::QbdRapResult<T> r = lf ? mam::qbd_rap(A0, A1, A2, B0, B1, nl) : mam::qbd_rap(A0, A1, A2);
        Json o;
        o["levelProb"] = encode_vector(r.levelProb);
        o["QN"] = encode_scalar(r.QN);
        o["R"] = encode_matrix(r.R);
        o["G"] = encode_matrix(r.G);
        o["U"] = encode_matrix(r.U);
        o["spr"] = r.spr;
        o["pqueue"] = encode_matrix(r.pqueue);
        o["pi0"] = encode_vector(r.pi0);
        return o;
    }
};
 
struct OpQbdBmapbmap1 {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> arr = read_map<T>(a, "MAPa");
        const std::vector<T> pb = a.vector<T>("pbatcha");
        const mam::Map<T> srv = read_map<T>(a, "MAPs");
        a.done();
        const mam::QbdBmapBmap1Blocks<T> r = mam::qbd_bmapbmap1(arr, pb, srv);
        Json o;
        o["A1"] = encode_matrices(r.A1);
        o["A0"] = encode_matrix(r.A0);
        o["A_1"] = encode_matrix(r.Am1);
        o["A0bar"] = encode_matrix(r.A0bar);
        o["B0"] = encode_matrix(r.B0);
        o["B1"] = encode_matrices(r.B1);
        return o;
    }
};
 
struct OpLdqbd {
    template <class T>
    static Json run(Args& a) {
        const std::vector<Matrix<T> > q0 = a.matrices<T>("Q0");
        const std::vector<Matrix<T> > q1 = a.matrices<T>("Q1");
        // MATLAB holds Q2{1..N} = the down-blocks of levels 1..N; the port indexes
        // them by level, 0..N, with slot 0 unused.
        std::vector<Matrix<T> > q2(1, Matrix<T>());
        for (const Matrix<T>& m : a.matrices<T>("Q2")) q2.push_back(m);
        a.unsupported("options", "its fields (epsilon, maxIter, verbose) are only printed by the "
                                 "MATLAB reference under verbose and never change the result; the "
                                 "recursion is exact and takes no tolerance.");
        a.done();
        const mam::LdqbdResult<T> r = mam::ldqbd(q0, q1, q2);
        std::vector<Matrix<T> > Rout;
        for (std::size_t n = 1; n < r.R.size(); ++n) Rout.push_back(r.R[n]);
        Json o;
        o["R"] = encode_matrices(Rout);
        o["pi"] = encode_vector(r.pi.pi);
        return o;
    }
};
 
struct OpLdqbdCensor {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        // phases: one vector of per-level phase counts per set (a MATLAB cell).
        std::vector<std::vector<std::size_t> > phases;
        for (const Matrix<double>& m : a.matrices<double>("phases")) {
            std::vector<std::size_t> P;
            for (std::size_t i = 0; i < m.rows(); ++i)
                for (std::size_t j = 0; j < m.cols(); ++j) {
                    const double v = m(i, j);
                    if (!(v >= 1) || v != static_cast<double>(static_cast<std::size_t>(v)))
                        throw InputError(ErrorKind{"Argument"},
                                         "ldqbd_censor: phases must hold positive integers");
                    P.push_back(static_cast<std::size_t>(v));
                }
            phases.push_back(P);
        }
        double tol = 1e-10;
        a.unsupported("options", "only options.tol is read by the reference and it is passed as 'tol'.");
        if (a.has("tol")) tol = a.number<double>("tol");
        a.done();
        const mam::LdqbdCensorResult<T> r = mam::ldqbd_censor(Q, phases, tol);
        Json cells = Json::array();
        for (const auto& set : r.pi_cells) cells.push_back(encode_vectors(set));
        Json o;
        o["pi"] = encode_vector(r.pi);
        o["pi_cell"] = cells;
        return o;
    }
};
 
struct OpLdqbdMphc {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> D0 = a.matrix<T>("D0");
        const Matrix<T> D1 = a.matrix<T>("D1");
        const std::vector<T> alpha = a.vector<T>("alpha");
        const double c = a.number<double>("c");
        const std::vector<T> arr = a.vector<T>("arrRate");
        const std::vector<T> sf = a.vector_or_empty<T>("sf");
        const bool hc = a.has("carry");
        const bool carry = a.boolean("carry", false);
        a.done();
        const mam::LdqbdMphcBlocks<T> r =
            hc ? mam::ldqbd_mphc(D0, D1, alpha, c, arr, sf, carry) : mam::ldqbd_mphc(D0, D1, alpha, c, arr, sf);
        // MATLAB's Q2 holds the down-blocks of levels 1..Nlev; the port's slot 0 is unused.
        std::vector<Matrix<T> > q2;
        for (std::size_t n = 1; n < r.Q2.size(); ++n) q2.push_back(r.Q2[n]);
        Json o;
        o["Q0"] = encode_matrices(r.Q0);
        o["Q1"] = encode_matrices(r.Q1);
        o["Q2"] = encode_matrices(q2);
        return o;
    }
};
 
// ---------------------------------------------------------------------------
// Setup / delay-off queues
// ---------------------------------------------------------------------------
 
struct OpQbdSetupdelayoff {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T mu = a.number<T>("mu");
        const T ar = a.number<T>("alpharate");
        const T as = a.number<T>("alphascv");
        const T br = a.number<T>("betarate");
        const T bs = a.number<T>("betascv");
        a.done();
        Json o;
        o["QN"] = encode_scalar(mam::qbd_setupdelayoff(lambda, mu, ar, as, br, bs));
        return o;
    }
};
 
struct OpQbdSetupdelayoffClosed {
    template <class T>
    static Json run(Args& a) {
        const T N = a.number<T>("N");
        const T Z = a.number<T>("Z");
        const T mu = a.number<T>("mu");
        const T ar = a.number<T>("alpharate");
        const T as = a.number<T>("alphascv");
        const T br = a.number<T>("betarate");
        const T bs = a.number<T>("betascv");
        a.done();
        const mam::SetupDelayoffClosed<T> r = mam::qbd_setupdelayoff_closed(N, Z, mu, ar, as, br, bs);
        Json o;
        o["QN"] = encode_scalar(r.QN);
        o["XN"] = encode_scalar(r.XN);
        return o;
    }
};
 
struct OpQbdSetupdelayoffMmk {
    template <class T>
    static Json run(Args& a) {
        const int k = a.required_integer("k");
        const double N = a.number<double>("N");
        const T arv = a.number<T>("arv");
        const T mu = a.number<T>("mu");
        const T ar = a.number<T>("alpharate");
        const T as = a.number<T>("alphascv");
        const T br = a.number<T>("betarate");
        const T bs = a.number<T>("betascv");
        a.done();
        const mam::SetupDelayoffMMk<T> r = mam::qbd_setupdelayoff_mmk(k, N, arv, mu, ar, as, br, bs);
        Json o;
        o["QN"] = encode_scalar(r.QN);
        o["XN"] = encode_scalar(r.XN);
        return o;
    }
};
 
// ---------------------------------------------------------------------------
// AMAP(2) / APH(2)
// ---------------------------------------------------------------------------
 
struct OpAmap2AdjustGamma {
    template <class T>
    static Json run(Args& a) {
        const T M1 = a.number<T>("M1");
        const T M2 = a.number<T>("M2");
        const T M3 = a.number<T>("M3");
        const T G = a.number<T>("GAMMA");
        const bool full = long_form(a, "amap2_adjust_gamma", {"weights", "constraints", "tol"},
                                    {"weights", "method", "constraints", "tol"});
        const bool hm = a.has("method");
        mam::Amap2AdjustGammaResult<T> r;
        if (full) {
            const std::vector<T> w = a.vector<T>("weights");
            const int method = a.required_integer("method");
            const int cons = a.required_integer("constraints");
            const T tol = a.number<T>("tol");
            a.done();
            r = mam::amap2_adjust_gamma(M1, M2, M3, G, w, method, cons, tol);
        } else if (hm) {
            const int method = a.required_integer("method");
            a.opt("weights");
            a.opt("constraints");
            a.opt("tol");
            a.done();
            r = mam::amap2_adjust_gamma(M1, M2, M3, G, method);
        } else {
            a.opt("weights");
            a.opt("method");
            a.opt("constraints");
            a.opt("tol");
            a.done();
            r = mam::amap2_adjust_gamma(M1, M2, M3, G);
        }
        Json o;
        o["M2a"] = encode_scalar(r.M2a);
        o["M3a"] = encode_scalar(r.M3a);
        o["GAMMAa"] = encode_scalar(r.GAMMAa);
        o["objective"] = encode_scalar(r.objective);
        o["feasible"] = r.feasible;
        return o;
    }
};
 
struct OpAmap2Assemble {
    template <class T>
    static Json run(Args& a) {
        const T l1 = a.number<T>("l1");
        const T l2 = a.number<T>("l2");
        const T p1 = a.number<T>("p1");
        const T p2 = a.number<T>("p2");
        const int form = a.required_integer("form");
        a.done();
        Json o;
        o["AMAP"] = encode_map(mam::amap2_assemble(l1, l2, p1, p2, form));
        return o;
    }
};
 
struct OpAmap2FitGamma {
    template <class T>
    static Json run(Args& a) {
        const T M1 = a.number<T>("M1");
        const T M2 = a.number<T>("M2");
        const T M3 = a.number<T>("M3");
        const T G = a.number<T>("GAMMA");
        const bool ht = a.has("cvtol");
        const T cvtol = ht ? a.number<T>("cvtol") : T();
        a.opt("cvtol");
        a.done();
        const mam::Amap2FitGammaResult<T> r =
            ht ? mam::amap2_fit_gamma(M1, M2, M3, G, cvtol) : mam::amap2_fit_gamma(M1, M2, M3, G);
        Json o;
        o["AMAP"] = encode_map(r.amap);
        o["AMAPS"] = encode_maps(r.amaps);
        o["adjusted"] = r.adjusted;
        o["poisson_fallback"] = r.poisson_fallback;
        return o;
    }
};
 
struct OpAmap2FitallGamma {
    template <class T>
    static Json run(Args& a) {
        const T M1 = a.number<T>("M1");
        const T M2 = a.number<T>("M2");
        const T M3 = a.number<T>("M3");
        const T G = a.number<T>("GAMMA");
        const bool lf = long_form(a, "amap2_fitall_gamma", {"degentol", "r12tol"}, {"degentol", "r12tol"});
        const T dt = lf ? a.number<T>("degentol") : T();
        const T rt = lf ? a.number<T>("r12tol") : T();
        a.opt("degentol");
        a.opt("r12tol");
        a.done();
        Json o;
        o["AMAPS"] = encode_maps(lf ? mam::amap2_fitall_gamma(M1, M2, M3, G, dt, rt)
                                    : mam::amap2_fitall_gamma(M1, M2, M3, G));
        return o;
    }
};
 
struct OpAph2Adjust {
    template <class T>
    static Json run(Args& a) {
        const T M1 = a.number<T>("M1");
        const T M2 = a.number<T>("M2");
        const T M3 = a.number<T>("M3");
        const std::string method = a.text("method", "simple");
        if (a.has("method") && method != "simple")
            throw UnsupportedError(ErrorKind{"JsonApi"},
                                   "aph2_adjust: only method \"simple\" is ported to C++; got \"" +
                                       method + "\".");
        const bool ht = a.has("tol");
        const T tol = ht ? a.number<T>("tol") : T();
        a.opt("tol");
        a.done();
        const mam::Aph2AdjustResult<T> r = ht ? mam::aph2_adjust(M1, M2, M3, tol) : mam::aph2_adjust(M1, M2, M3);
        Json o;
        o["M2a"] = encode_scalar(r.M2a);
        o["M3a"] = encode_scalar(r.M3a);
        return o;
    }
};
 
struct OpAph2Assemble {
    template <class T>
    static Json run(Args& a) {
        const T l1 = a.number<T>("l1");
        const T l2 = a.number<T>("l2");
        const T p1 = a.number<T>("p1");
        a.done();
        Json o;
        o["APH"] = encode_map(mam::aph2_assemble(l1, l2, p1));
        return o;
    }
};
 
struct OpAph2Fit {
    template <class T>
    static Json run(Args& a) {
        const T M1 = a.number<T>("M1");
        const T M2 = a.number<T>("M2");
        const T M3 = a.number<T>("M3");
        a.done();
        const mam::Aph2FitResult<T> r = mam::aph2_fit(M1, M2, M3);
        Json o;
        o["APH"] = encode_map(r.aph);
        o["APHS"] = encode_maps(r.aphs);
        o["adjusted"] = r.adjusted;
        return o;
    }
};
 
struct OpAph2Fitall {
    template <class T>
    static Json run(Args& a) {
        const T M1 = a.number<T>("M1");
        const T M2 = a.number<T>("M2");
        const T M3 = a.number<T>("M3");
        const bool ht = a.has("degentol");
        const T dt = ht ? a.number<T>("degentol") : T();
        a.opt("degentol");
        a.done();
        Json o;
        o["APHS"] = encode_maps(ht ? mam::aph2_fitall(M1, M2, M3, dt) : mam::aph2_fitall(M1, M2, M3));
        return o;
    }
};
 
// ---------------------------------------------------------------------------
// MAP(2) / MMPP(2) fitting
// ---------------------------------------------------------------------------
 
struct OpMap2Fit {
    template <class T>
    static Json run(Args& a) {
        const T e1 = a.number<T>("e1");
        const T e2 = a.number<T>("e2");
        const bool h3 = a.has("e3");
        const T e3 = h3 ? a.number<T>("e3") : T();
        a.opt("e3");
        const T g2 = a.number<T>("g2");
        a.done();
        const mam::Map2FitResult<T> r = h3 ? mam::map2_fit(e1, e2, e3, g2) : mam::map2_fit(e1, e2, g2);
        Json o;
        o["MAP"] = r.has_map ? encode_map(r.map) : Json::array();
        o["ERR"] = r.err;
        return o;
    }
};
 
struct OpMap2FitIdc {
    template <class T>
    static Json run(Args& a) {
        const T e1 = a.number<T>("e1");
        const T e2 = a.number<T>("e2");
        const T e3 = a.number<T>("e3");
        const T I = a.number<T>("I");
        a.done();
        const mam::Map2FitIdcResult<T> r = mam::map2_fit_idc(e1, e2, e3, I);
        Json o;
        o["MAP"] = encode_map(r.map);
        o["status"] = r.status;
        return o;
    }
};
 
struct OpMmpp2Fit {
    template <class T>
    static Json run(Args& a) {
        const T E1 = a.number<T>("E1");
        const T E2 = a.number<T>("E2");
        const T E3 = a.number<T>("E3");
        const T g = a.number<T>("ACFLAG1");
        a.done();
        Json o;
        o["MAP"] = encode_map(mam::mmpp2_fit(E1, E2, E3, g));
        return o;
    }
};
 
struct OpMmpp2Fit1 {
    template <class T>
    static Json run(Args& a) {
        const T mean = a.number<T>("mean");
        const T scv = a.number<T>("scv");
        const T skew = a.number<T>("skew");
        const T idc = a.number<T>("idc");
        a.done();
        const mam::Map2FitResult<T> r = mam::mmpp2_fit1(mean, scv, skew, idc);
        Json o;
        o["MAP"] = r.has_map ? encode_map(r.map) : Json::array();
        return o;
    }
};
 
template <bool Four>
struct OpMmpp2Fit24 {
    template <class T>
    static Json run(Args& a) {
        const T mean = a.number<T>("mean");
        const T scv = a.number<T>("scv");
        const T skew = a.number<T>("skew");
        const T g = a.number<T>(Four ? "acf1" : "g2");
        a.done();
        const mam::Mmpp2FitResult<T> r = Four ? mam::mmpp2_fit4(mean, scv, skew, g) : mam::mmpp2_fit2(mean, scv, skew, g);
        Json o;
        o["MAP"] = encode_map(r.map);
        o["feasible"] = r.feasible;
        return o;
    }
};
 
struct OpMmpp2Fit3 {
    template <class T>
    static Json run(Args& a) {
        const T E1 = a.number<T>("E1");
        const T E2 = a.number<T>("E2");
        const T E3 = a.number<T>("E3");
        const T G2 = a.number<T>("G2");
        const bool ht = a.has("g2tol");
        const T tol = ht ? a.number<T>("g2tol") : T();
        a.opt("g2tol");
        a.done();
        Json o;
        o["MAP"] = encode_map(ht ? mam::mmpp2_fit3(E1, E2, E3, G2, tol) : mam::mmpp2_fit3(E1, E2, E3, G2));
        return o;
    }
};
 
struct OpMmpp2Fitc {
    template <class T>
    static Json run(Args& a) {
        const T mu = a.number<T>("mu");
        const T bt1 = a.number<T>("bt1");
        const T bt2 = a.number<T>("bt2");
        const T binf = a.number<T>("binf");
        const T m3t2 = a.number<T>("m3t2");
        const T t1 = a.number<T>("t1");
        const T t2 = a.number<T>("t2");
        a.done();
        const mam::Mmpp2FitcResult<T> r = mam::mmpp2_fitc(mu, bt1, bt2, binf, m3t2, t1, t2);
        Json o;
        o["mmpp"] = encode_map(r.map);
        o["degenerate"] = r.degenerate;
        o["third_moment_ok"] = r.third_moment_ok;
        return o;
    }
};
 
struct OpMmpp2FitcApprox {
    template <class T>
    static Json run(Args& a) {
        const T x = a.number<T>("a");
        const T bt1 = a.number<T>("bt1");
        const T bt2 = a.number<T>("bt2");
        const T binf = a.number<T>("binf");
        const T m3t2 = a.number<T>("m3t2");
        const T t1 = a.number<T>("t1");
        const T t2 = a.number<T>("t2");
        a.done();
        const mam::Mmpp2FitcApproxResult<T> r = mam::mmpp2_fitc_approx(x, bt1, bt2, binf, m3t2, t1, t2);
        Json o;
        o["FIT"] = encode_map(r.map);
        o["objective"] = encode_scalar(r.objective);
        o["violation"] = encode_scalar(r.violation);
        o["converged"] = r.converged;
        return o;
    }
};
 
struct OpMmpp2FitcTheoretical {
    template <class T>
    static Json run(Args& a) {
        const mam::Map<T> m = read_map<T>(a, "map");
        const bool lf = long_form(a, "mmpp2_fitc_theoretical", {"t1", "t2", "tinf"}, {"t1", "t2", "tinf"});
        const T t1 = lf ? a.number<T>("t1") : T();
        const T t2 = lf ? a.number<T>("t2") : T();
        const T ti = lf ? a.number<T>("tinf") : T();
        a.opt("t1");
        a.opt("t2");
        a.opt("tinf");
        a.done();
        const mam::Mmpp2FitcResult<T> r = lf ? mam::mmpp2_fitc_theoretical(m, t1, t2, ti) : mam::mmpp2_fitc_theoretical(m);
        Json o;
        o["fit"] = encode_map(r.map);
        o["degenerate"] = r.degenerate;
        o["third_moment_ok"] = r.third_moment_ok;
        return o;
    }
};
 
// ---------------------------------------------------------------------------
// M3PP
// ---------------------------------------------------------------------------
 
struct OpM3ppSuperposFitc {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> av = a.vector<T>("av");
        const std::vector<T> btv = a.vector<T>("btv");
        const std::vector<T> binfv = a.vector<T>("binfv");
        const std::vector<T> m3tv = a.vector<T>("m3tv");
        const T t = a.number<T>("t");
        const T tinf = a.number<T>("tinf");
        a.done();
        const mam::M3ppSuperposResult<T> r = mam::m3pp_superpos_fitc(av, btv, binfv, m3tv, t, tinf);
        Json o;
        o["fit"] = encode_mmap(r.mmap);
        o["m3pps"] = encode_mmaps(r.parts);
        return o;
    }
};
 
/** The ten characteristics m3pp2m_fitc and its approximate variants share. */
template <class T>
struct M3pp2mChars {
    T a, bt1, bt2, binf, m3t2, t1, t2, t3;
    std::vector<T> ai, third;
};
 
template <class T>
M3pp2mChars<T> read_m3pp2m(Args& a, const char* third) {
    M3pp2mChars<T> c;
    c.a = a.number<T>("a");
    c.bt1 = a.number<T>("bt1");
    c.bt2 = a.number<T>("bt2");
    c.binf = a.number<T>("binf");
    c.m3t2 = a.number<T>("m3t2");
    c.t1 = a.number<T>("t1");
    c.t2 = a.number<T>("t2");
    c.ai = a.vector<T>("ai");
    c.third = a.vector<T>(third);
    c.t3 = a.number<T>("t3");
    a.done();
    return c;
}
 
struct OpM3pp2mFitc {
    template <class T>
    static Json run(Args& a) {
        const M3pp2mChars<T> c = read_m3pp2m<T>(a, "dvt3");
        const mam::M3pp2mFitcResult<T> r =
            mam::m3pp2m_fitc(c.a, c.bt1, c.bt2, c.binf, c.m3t2, c.t1, c.t2, c.ai, c.third, c.t3);
        Json o;
        o["mmmpp"] = encode_mmap(r.mmap);
        o["degenerate"] = r.degenerate;
        return o;
    }
};
 
template <bool Ag>
struct OpM3pp2mFitcApprox {
    template <class T>
    static Json run(Args& a) {
        const M3pp2mChars<T> c = read_m3pp2m<T>(a, Ag ? "gt3" : "dvt3");
        const mam::M3pp2mFitcApproxResult<T> r =
            Ag ? mam::m3pp2m_fitc_approx_ag(c.a, c.bt1, c.bt2, c.binf, c.m3t2, c.t1, c.t2, c.ai, c.third, c.t3)
               : mam::m3pp2m_fitc_approx(c.a, c.bt1, c.bt2, c.binf, c.m3t2, c.t1, c.t2, c.ai, c.third, c.t3);
        Json o;
        o["FIT"] = encode_mmap(r.mmap);
        o["mmpp_objective"] = encode_scalar(r.mmpp_objective);
        o["class_objective"] = encode_scalar(r.class_objective);
        o["class_violation"] = encode_scalar(r.class_violation);
        o["degenerate"] = r.degenerate;
        o["feasible"] = r.feasible;
        return o;
    }
};
 
struct OpM3pp2mFitcTrace {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> Tv = a.vector<T>("T");
        const std::vector<int> A = a.ints("A");
        const Json& jm = a.get("method");
        if (!jm.is_string()) throw InputError("'method' must be a string");
        const std::string method = jm.get<std::string>();
        const bool lf = long_form(a, "m3pp2m_fitc_trace", {"t1", "tinf"}, {"t1", "tinf"});
        const T t1 = lf ? a.number<T>("t1") : T();
        const T ti = lf ? a.number<T>("tinf") : T();
        a.opt("t1");
        a.opt("tinf");
        a.done();
        const mam::M3pp2mFitcTraceResult<T> r =
            lf ? mam::m3pp2m_fitc_trace(Tv, A, method, t1, ti) : mam::m3pp2m_fitc_trace(Tv, A, method);
        Json o;
        o["fit"] = encode_mmap(r.mmap);
        return o;
    }
};
 
// ---------------------------------------------------------------------------
// Fluid queues
// ---------------------------------------------------------------------------
 
struct OpMfqSojourn {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const Matrix<T> Rin = a.matrix<T>("Rin");
        const Matrix<T> Rout = a.matrix<T>("Rout");
        const bool lf = long_form(a, "mfq_sojourn", {"Q0", "transToPH", "prec"}, {"transToPH", "prec"});
        const Matrix<T> Q0 = a.matrix_or_empty<T>("Q0");
        const bool tp = a.boolean("transToPH", false);
        const T prec = lf ? a.number<T>("prec") : T();
        a.opt("prec");
        a.done();
        const mam::MeRepresentation<T> r = lf ? mam::mfq_sojourn(Q, Rin, Rout, Q0, tp, prec) : mam::mfq_sojourn(Q, Rin, Rout);
        Json o;
        o["alpha"] = encode_vector(r.alpha);
        o["A"] = encode_matrix(r.A);
        return o;
    }
};
 
struct OpMfqFlufluSojourn {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Qin = a.matrix<T>("Qin");
        const Matrix<T> Rin = a.matrix<T>("Rin");
        const Matrix<T> Qout = a.matrix<T>("Qout");
        const Matrix<T> Rout = a.matrix<T>("Rout");
        if (!a.has("srv0stop")) a.get("srv0stop");
        const bool s0 = a.boolean("srv0stop", false);
        const bool lf = long_form(a, "mfq_fluflu_sojourn", {"transToPH", "prec"}, {"transToPH", "prec"});
        const bool tp = a.boolean("transToPH", false);
        const T prec = lf ? a.number<T>("prec") : T();
        a.opt("prec");
        a.done();
        const mam::MeRepresentation<T> r = lf ? mam::mfq_fluflu_sojourn(Qin, Rin, Qout, Rout, s0, tp, prec)
                                              : mam::mfq_fluflu_sojourn(Qin, Rin, Qout, Rout, s0);
        Json o;
        o["alpha"] = encode_vector(r.alpha);
        o["A"] = encode_matrix(r.A);
        return o;
    }
};
 
struct OpMfqLdMean {
    template <class T>
    static Json run(Args& a) {
        const mam::LevelDependentFluidBlocks<T> b = read_ld_blocks<T>(a);
        a.done();
        Json o;
        o["res"] = encode_scalar(mam::mfq_ld_mean(b));
        return o;
    }
};
 
struct OpMfqLdDistr {
    template <class T>
    static Json run(Args& a) {
        const mam::LevelDependentFluidBlocks<T> b = read_ld_blocks<T>(a);
        const mam::FluidDistrKind what = distr_kind(a);
        const std::vector<T> pts = a.vector<T>("points");
        a.done();
        Json o;
        o["res"] = encode_vectors(mam::mfq_ld_distr(b, what, pts));
        return o;
    }
};
 
struct OpMfqLdSolve {
    template <class T>
    static Json run(Args& a) {
        const std::vector<Matrix<T> > Q = a.matrices<T>("Q");
        const std::vector<Matrix<T> > R = a.matrices<T>("R");
        const std::vector<Matrix<T> > S = a.matrices<T>("S");
        const std::vector<T> Thr = a.vector<T>("T");
        const bool lf = long_form(a, "mfq_ld_solve", {"boundaryL", "boundaryU", "Qt", "prec"}, {"prec"});
        const std::vector<mam::FluidBoundary> bL = boundaries(a, "boundaryL");
        const std::vector<mam::FluidBoundary> bU = boundaries(a, "boundaryU");
        const std::vector<Matrix<T> > Qt = a.matrices_or_empty<T>("Qt");
        const T prec = lf ? a.number<T>("prec") : T();
        a.opt("prec");
        a.done();
        const mam::LevelDependentFluidBlocks<T> r =
            lf ? mam::mfq_ld_solve(Q, R, S, Thr, bL, bU, Qt, prec) : mam::mfq_ld_solve(Q, R, S, Thr);
        Json o;
        o["masses"] = encode_vectors(r.masses);
        o["iniF"] = encode_vectors(r.iniF);
        o["KF"] = encode_matrices(r.KF);
        o["cloF"] = encode_matrices(r.cloF);
        o["iniB"] = encode_vectors(r.iniB);
        o["KB"] = encode_matrices(r.KB);
        o["cloB"] = encode_matrices(r.cloB);
        return o;
    }
};
 
struct OpMfqMultiregime {
    static Json run(Args& a) {
        const std::vector<Matrix<double> > Q = a.matrices<double>("Q");
        const std::vector<std::vector<double> > R = diagonal_rate_list(a, "R", false);
        const std::vector<Matrix<double> > Qt = a.matrices_or_empty<double>("Qt");
        const std::vector<std::vector<double> > Rt = diagonal_rate_list(a, "Rt", true);
        const std::vector<double> Thr = a.vector<double>("T");
        const std::vector<double> pdfp = a.vector_or_empty<double>("pdfpoints");
        const std::vector<double> cdfp = a.vector_or_empty<double>("cdfpoints");
        a.done();
        const mam::MultiRegimeResult r = mam::mfq_multiregime(Q, R, Qt, Rt, Thr, pdfp, cdfp);
        Json o;
        o["pdf"] = encode_vectors(r.pdf);
        o["pdfd"] = encode_vectors(r.pdfd);
        o["cdf"] = encode_vectors(r.cdf);
        o["cdfm"] = encode_vectors(r.cdfm);
        return o;
    }
};
 
struct OpMfqPrioQueue {
    static Json run(Args& a) {
        const Matrix<double> Q = a.matrix<double>("Q");
        const Matrix<double> R = a.matrix<double>("R");
        const double d = a.number<double>("d");
        mam::FluidPrioOptions opt;
        const bool fm = a.has("flMoms"), fd = a.has("flDistr"), sm = a.has("stMoms"), sd = a.has("stDistr");
        if (fm) opt.flMoms = a.count("flMoms");
        if (fd) opt.flDistr = a.vector<double>("flDistr");
        if (sm) opt.stMoms = a.count("stMoms");
        if (sd) opt.stDistr = a.vector<double>("stDistr");
        a.opt("flMoms");
        a.opt("flDistr");
        a.opt("stMoms");
        a.opt("stDistr");
        if (a.has("prec")) opt.prec = a.number<double>("prec");
        a.opt("prec");
        if (a.has("erlMaxOrder")) opt.erlMaxOrder = a.count("erlMaxOrder");
        a.opt("erlMaxOrder");
        if (a.has("classes")) {
            for (int c : a.ints("classes")) {
                if (c < 1) throw InputError("'classes' are 1-based, as in MATLAB");
                opt.classes.push_back(static_cast<std::size_t>(c));
            }
        }
        a.opt("classes");
        a.done();
        if (!fm && !fd && !sm && !sd)
            throw InputError("mfq_prio_queue: ask for at least one measure: flMoms, flDistr, "
                             "stMoms or stDistr");
        const mam::FluidPrioResult r = mam::mfq_prio_queue(Q, R, d, opt);
        Json o;
        Json cls = Json::array();
        for (std::size_t c : r.classes) cls.push_back(static_cast<std::uint64_t>(c));
        o["classes"] = cls;
        if (fm) o["flMoms"] = encode_vectors(r.flMoms);
        if (fd) o["flDistr"] = encode_vectors(r.flDistr);
        if (sm) o["stMoms"] = encode_vectors(r.stMoms);
        if (sd) o["stDistr"] = encode_vectors(r.stDistr);
        return o;
    }
};
 
// ---------------------------------------------------------------------------
// Semi-Markov passage times
// ---------------------------------------------------------------------------
 
struct OpSmpPassageMoments {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> P = a.matrix<T>("P");
        const Json& jh = a.get("hmom");
        const std::vector<T> pi0 = a.vector_or_empty<T>("pi0");
        const std::vector<std::size_t> target = one_based(a.get("target"), "'target'");
        const bool hn = a.has("nmax");
        const std::size_t nmax = a.count("nmax", 1);
        a.done();
        // hmom is either (nstates x nmax), or the per-pair cell {i,k} of moment vectors.
        const bool cell = jh.is_array() && !jh.empty() && jh[0].is_array() && !jh[0].empty() &&
                          jh[0][0].is_array();
        mc::PassageMoments<T> r;
        if (cell) {
            std::vector<std::vector<std::vector<T> > > h;
            for (std::size_t i = 0; i < jh.size(); ++i)
                h.push_back(vectors_from_json<T>(jh[i], "'hmom'[" + std::to_string(i) + "]"));
            r = hn ? mc::smp_passage_moments(P, h, pi0, target, nmax) : mc::smp_passage_moments(P, h, pi0, target);
        } else {
            const Matrix<T> h = matrix_from_json<T>(jh, "'hmom'");
            r = hn ? mc::smp_passage_moments(P, h, pi0, target, nmax) : mc::smp_passage_moments(P, h, pi0, target);
        }
        Json o;
        o["Mall"] = encode_matrix(r.mall);
        o["m"] = encode_vector(r.m);
        return o;
    }
};
 
// ---------------------------------------------------------------------------
// Maximum entropy (api/me)
// ---------------------------------------------------------------------------
 
struct OpMeCqn {
    template <class T>
    static Json run(Args& a) {
        const std::size_t M = a.count("M");
        const std::size_t R = a.count("R");
        const std::vector<long> N = populations(a, "N", false);
        const Matrix<T> mu = a.matrix<T>("mu");
        const Matrix<T> Cs = a.matrix<T>("Cs");
        const std::vector<Matrix<T> > P = per_class_routing<T>(a, "P", R);
        const std::vector<long> c = counts_with_inf(a, "c");
        const std::vector<long> ref = reference_stations(a, "refstat", R);
        const std::vector<char> ins = flags(a, "insens");
        const bool ho = a.has("options");
        const me::MeOptions opt = me_options(a);
        a.done();
        const me::MeResult<T> r = ho ? me::me_cqn(M, R, N, mu, Cs, P, c, ref, ins, opt)
                                     : me::me_cqn(M, R, N, mu, Cs, P, c, ref, ins);
        return encode_me_result(r, true);
    }
};
 
struct OpMeMqn {
    template <class T>
    static Json run(Args& a) {
        const std::size_t M = a.count("M");
        const std::size_t R = a.count("R");
        const std::vector<char> oc = flags(a, "openClasses");
        const Matrix<T> lambda0 = a.matrix<T>("lambda0");
        const Matrix<T> Ca0 = a.matrix<T>("Ca0");
        const std::vector<long> N = populations(a, "N", true);
        const Matrix<T> mu = a.matrix<T>("mu");
        const Matrix<T> Cs = a.matrix<T>("Cs");
        const std::vector<Matrix<T> > P = per_class_routing<T>(a, "P", R);
        const std::vector<long> c = counts_with_inf(a, "c");
        const std::vector<long> ref = reference_stations(a, "refstat", R);
        const std::vector<char> ins = flags(a, "insens");
        const bool ho = a.has("options");
        const me::MeOptions opt = me_options(a);
        a.done();
        const me::MeResult<T> r = ho ? me::me_mqn(M, R, oc, lambda0, Ca0, N, mu, Cs, P, c, ref, ins, opt)
                                     : me::me_mqn(M, R, oc, lambda0, Ca0, N, mu, Cs, P, c, ref, ins);
        return encode_me_result(r, true);
    }
};
 
struct OpMeOqn {
    template <class T>
    static Json run(Args& a) {
        const std::size_t M = a.count("M");
        const std::size_t R = a.count("R");
        const Matrix<T> lambda0 = a.matrix<T>("lambda0");
        const Matrix<T> Ca0 = a.matrix<T>("Ca0");
        const Matrix<T> mu = a.matrix<T>("mu");
        const Matrix<T> Cs = a.matrix<T>("Cs");
        const std::vector<Matrix<T> > P = per_class_routing<T>(a, "P", R);
        const std::vector<long> c = counts_with_inf(a, "c");
        const std::vector<char> ins = flags(a, "insens");
        const bool ho = a.has("options");
        const me::MeOptions opt = me_options(a);
        a.done();
        const me::MeResult<T> r = ho ? me::me_oqn(M, R, lambda0, Ca0, mu, Cs, P, c, ins, opt)
                                     : me::me_oqn(M, R, lambda0, Ca0, mu, Cs, P, c, ins);
        return encode_me_result(r, false);
    }
};
 
struct OpMeOqnBlk {
    template <class T>
    static Json run(Args& a) {
        const std::size_t M = a.count("M");
        const std::vector<T> lambda0 = a.vector<T>("lambda0");
        const std::vector<T> Ca0 = a.vector<T>("Ca0");
        const std::vector<T> mu = a.vector<T>("mu");
        const std::vector<T> Cs = a.vector<T>("Cs");
        const Matrix<T> P = a.matrix<T>("P");
        const std::vector<long> c = counts_with_inf(a, "c");
        const std::vector<long> N = counts_with_inf(a, "N");
        const std::vector<int> br = a.ints("blockrule");
        const bool ho = a.has("options");
        const me::MeBlkOptions opt = me_blk_options(a);
        a.done();
        const me::MeBlkResult<T> r = ho ? me::me_oqn_blk(M, lambda0, Ca0, mu, Cs, P, c, N, br, opt)
                                        : me::me_oqn_blk(M, lambda0, Ca0, mu, Cs, P, c, N, br);
        Json o;
        o["Q"] = encode_vector(r.Q);
        o["W"] = encode_vector(r.W);
        o["T"] = encode_vector(r.T_);
        o["U"] = encode_vector(r.U);
        o["Ca"] = encode_vector(r.Ca);
        o["Cd"] = encode_vector(r.Cd);
        o["PBa"] = encode_vector(r.PBa);
        o["lambda"] = encode_vector(r.lambda);
        o["iter"] = static_cast<std::int64_t>(r.iter);
        o["converged"] = r.converged;
        return o;
    }
};
 
struct OpMeGegecMql {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T Ca = a.number<T>("Ca");
        const T mu = a.number<T>("mu");
        const T Cs = a.number<T>("Cs");
        const long c = a.required_integer("c");
        a.done();
        Json o;
        o["L"] = encode_scalar(me::me_gegec_mql(lambda, Ca, mu, Cs, c));
        return o;
    }
};
 
struct OpMeGegecn {
    template <class T>
    static Json run(Args& a) {
        const T lambda = a.number<T>("lambda");
        const T Ca = a.number<T>("Ca");
        const T mu = a.number<T>("mu");
        const T Cs = a.number<T>("Cs");
        const long c = a.required_integer("c");
        const long K = a.required_integer("K");
        const long N = a.required_integer("N");
        a.done();
        const me::GegecnResult<T> r = me::me_gegecn(lambda, Ca, mu, Cs, c, K, N);
        Json o;
        o["p"] = encode_vector(r.p);
        o["L"] = encode_scalar(r.L);
        o["U"] = encode_scalar(r.U);
        o["PB"] = encode_scalar(r.PB);
        o["Lq"] = encode_scalar(r.Lq);
        return o;
    }
};
 
struct OpMeGegecnPb {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> p = a.vector<T>("p");
        const long K = a.required_integer("K");
        const long N = a.required_integer("N");
        const long c = a.required_integer("c");
        const T Cs = a.number<T>("Cs");
        const T Ca = a.number<T>("Ca");
        a.done();
        Json o;
        o["PB"] = encode_scalar(me::me_gegecn_pb(p, K, N, c, Cs, Ca));
        return o;
    }
};
 
 
/**
 * MMAP = m3pp_rand(order, classes, seed): a random M3PP. MATLAB's signature
 * has no seed and draws from the session RNG; the port takes its engine by
 * reference, so the JSON form requires the seed of a std::mt19937_64. The
 * result is the [D0, D1, D1_1, ..., D1_C] cell the other mam entries return.
 */
struct OpM3ppRand {
    template <class T>
    static Json run(Args& a) {
        const std::size_t K = a.count("order");
        const std::size_t classes = a.count("classes");
        std::mt19937_64 gen(required_seed(a, "m3pp_rand"));
        a.done();
        Json out;
        out["MMAP"] = encode_mmap(mam::m3pp_rand<T>(K, classes, gen));
        return out;
    }
};
 
}  // namespace
 
void register_api_mam(DispatchTable& t) {
    t["m3pp_rand"] = make<NeedsTranscendental<OpM3ppRand> >();
    // Double, Exact and Real.
    t["map_infgen"] = make<OpMapInfgen>();
    t["map_prob"] = make<OpMapProb>();
    t["map_lambda"] = make<OpMapLambda>();
    t["map_pie"] = make<OpMapPie>();
    t["map_mean"] = make<OpMapMean>();
    t["map_embedded"] = make<OpMapEmbedded>();
    t["map_moment"] = make<OpMapMoment>();
    t["map_var"] = make<OpMapVar>();
    t["map_scv"] = make<OpMapScv>();
    t["map_acf"] = make<OpMapAcf>();
    t["map_idc"] = make<OpMapIdc>();
    t["mmap_lambda"] = make<OpMmapLambda>();
    t["mmap_count_lambda"] = make<OpMmapCountLambda>();
    t["mmap_pc"] = make<OpMmapPc>();
    t["mmap_isfeasible"] = make<OpMmapIsfeasible>();
    t["mmap_normalize"] = make<OpMmapNormalize>();
    t["mmap_scale"] = make<OpMmapScale>();
    t["mmap_super"] = make<OpMmapSuper>();
    t["map_normalize"] = make<OpMapNormalize>();
    t["map_scale"] = make<OpMapScale>();
    t["map_exponential_mean"] = make<OpMapExponentialMean>();
    t["map_erlang"] = make<OpMapErlang>();
    t["map_sum"] = make<OpMapSum>();
    t["map_sumind"] = make<OpMapSumind>();
    t["map_mixture"] = make<OpMapMixture>();
    t["map_renewal"] = make<OpMapRenewal>();
    t["ph2map"] = make<OpPh2map>();
    t["map_stochcomp"] = make<OpMapStochcomp>();
    t["map_kurt"] = make<OpMapKurt>();
    t["map_joint"] = make<OpMapJoint>();
    t["map_isfeasible"] = make<OpMapIsfeasible>();
    t["qbd_pi"] = make<OpQbdPi>();
    t["qbd_R_residual"] = make<OpQbdResidual<false> >();
    t["qbd_G_residual"] = make<OpQbdResidual<true> >();
    t["qbd_mapmap1_blocks"] = make<OpQbdMapmap1Blocks>();
    t["qbd_qlen_factmoment"] = make<OpQbdQlenMoment<true> >();
    t["qbd_qlen_moment"] = make<OpQbdQlenMoment<false> >();
    t["amap2_assemble"] = make<OpAmap2Assemble>();
    t["aph2_assemble"] = make<OpAph2Assemble>();
    t["map_mark"] = make<OpMapMark>();
    t["map_count_mean"] = make<OpMapCountMean>();
    t["smp_passage_moments"] = make<OpSmpPassageMoments>();
    t["ldqbd"] = make<OpLdqbd>();
    t["ldqbd_censor"] = make<OpLdqbdCensor>();
    t["ldqbd_mphc"] = make<OpLdqbdMphc>();
    t["qbd_bmapbmap1"] = make<OpQbdBmapbmap1>();
 
    // Double and Real.
    t["me_cqn"] = make<NeedsTranscendental<OpMeCqn> >();
    t["me_mqn"] = make<NeedsTranscendental<OpMeMqn> >();
    t["me_oqn"] = make<NeedsTranscendental<OpMeOqn> >();
    t["qbd_R"] = make<NeedsTranscendental<OpQbdR<false> > >();
    t["qbd_R_logred"] = make<NeedsTranscendental<OpQbdR<true> > >();
    t["qbd_fundmat"] = make<NeedsTranscendental<OpQbdFundmat> >();
    t["qbd_caudal"] = make<NeedsTranscendental<OpQbdCaudal> >();
    t["qbd_rg"] = make<NeedsTranscendental<OpQbdRg> >();
    t["qbd_mapmap1"] = make<NeedsTranscendental<OpQbdMapmap1> >();
    t["qbd_mapmap1_qlen_truncated"] = make<NeedsTranscendental<OpQbdMapmap1QlenTruncated> >();
    t["map_hyperexp"] = make<NeedsTranscendental<OpMapHyperexp> >();
    t["map_skew"] = make<NeedsTranscendental<OpMapSkew> >();
    t["amap2_adjust_gamma"] = make<NeedsTranscendental<OpAmap2AdjustGamma> >();
    t["amap2_fit_gamma"] = make<NeedsTranscendental<OpAmap2FitGamma> >();
    t["amap2_fitall_gamma"] = make<NeedsTranscendental<OpAmap2FitallGamma> >();
    t["aph2_adjust"] = make<NeedsTranscendental<OpAph2Adjust> >();
    t["aph2_fit"] = make<NeedsTranscendental<OpAph2Fit> >();
    t["aph2_fitall"] = make<NeedsTranscendental<OpAph2Fitall> >();
    t["m3pp_superpos_fitc"] = make<NeedsTranscendental<OpM3ppSuperposFitc> >();
    t["m3pp2m_fitc"] = make<NeedsTranscendental<OpM3pp2mFitc> >();
    t["m3pp2m_fitc_approx"] = make<NeedsTranscendental<OpM3pp2mFitcApprox<false> > >();
    t["m3pp2m_fitc_approx_ag"] = make<NeedsTranscendental<OpM3pp2mFitcApprox<true> > >();
    t["m3pp2m_fitc_trace"] = make<NeedsTranscendental<OpM3pp2mFitcTrace> >();
    t["map2_fit"] = make<NeedsTranscendental<OpMap2Fit> >();
    t["mmpp2_fit"] = make<NeedsTranscendental<OpMmpp2Fit> >();
    t["mmpp2_fit1"] = make<NeedsTranscendental<OpMmpp2Fit1> >();
    t["mmpp2_fit2"] = make<NeedsTranscendental<OpMmpp2Fit24<false> > >();
    t["mmpp2_fit3"] = make<NeedsTranscendental<OpMmpp2Fit3> >();
    t["mmpp2_fit4"] = make<NeedsTranscendental<OpMmpp2Fit24<true> > >();
    t["mmpp2_fitc"] = make<NeedsTranscendental<OpMmpp2Fitc> >();
    t["mmpp2_fitc_approx"] = make<NeedsTranscendental<OpMmpp2FitcApprox> >();
    t["mmpp2_fitc_theoretical"] = make<NeedsTranscendental<OpMmpp2FitcTheoretical> >();
    t["map_acfc"] = make<NeedsTranscendental<OpMapAcfc> >();
    t["map_cdf"] = make<NeedsTranscendental<OpMapCdf> >();
    t["map_count_moment"] = make<NeedsTranscendental<OpMapCountMoment> >();
    t["map_count_var"] = make<NeedsTranscendental<OpMapCountVar> >();
    t["map_pdf"] = make<NeedsTranscendental<OpMapPdf> >();
    t["map_varcount"] = make<NeedsTranscendental<OpMapVarcount> >();
    t["mmap_count_var"] = make<NeedsTranscendental<OpMmapCountVar> >();
    t["mfq_sojourn"] = make<NeedsTranscendental<OpMfqSojourn> >();
    t["mfq_fluflu_sojourn"] = make<NeedsTranscendental<OpMfqFlufluSojourn> >();
    t["mfq_ld_mean"] = make<NeedsTranscendental<OpMfqLdMean> >();
    t["mfq_ld_distr"] = make<NeedsTranscendental<OpMfqLdDistr> >();
    t["map2_fit_idc"] = make<NeedsTranscendental<OpMap2FitIdc> >();
    t["me_gegec_mql"] = make<NeedsTranscendental<OpMeGegecMql> >();
    t["me_gegecn"] = make<NeedsTranscendental<OpMeGegecn> >();
    t["me_gegecn_pb"] = make<NeedsTranscendental<OpMeGegecnPb> >();
    t["me_oqn_blk"] = make<NeedsTranscendental<OpMeOqnBlk> >();
    t["mfq_ld_solve"] = make<NeedsTranscendental<OpMfqLdSolve> >();
    t["mmap_compress"] = make<NeedsTranscendental<OpMmapCompress> >();
    t["qbd_rap"] = make<NeedsTranscendental<OpQbdRap> >();
    t["qbd_setupdelayoff"] = make<NeedsTranscendental<OpQbdSetupdelayoff> >();
    t["qbd_setupdelayoff_closed"] = make<NeedsTranscendental<OpQbdSetupdelayoffClosed> >();
    t["qbd_setupdelayoff_mmk"] = make<NeedsTranscendental<OpQbdSetupdelayoffMmk> >();
 
    // Double only.
    t["mfq_multiregime"] = make<DoubleOnly<OpMfqMultiregime> >();
    t["mfq_prio_queue"] = make<DoubleOnly<OpMfqPrioQueue> >();
}
 
void register_library_only_mam(LibraryOnlyTable& t) {
    const std::string callables =
        "its holding-time transforms hlst are callables (std::function of a complex argument, "
        "one per state or per state pair), which have no JSON representation.";
    t["smp_passage_lst"] = "smp_passage_lst: " + callables;
    t["smp_passage_time"] = "smp_passage_time: " + callables;
    const std::string distribs =
        "its arcs are lang::Distrib objects (arbitrary distribution instances with their own "
        "sampling and transform methods), which have no JSON representation in the --api "
        "argument vocabulary.";
    t["smp_stage_cif"] = "smp_stage_cif: " + distribs;
    t["smp_stage_quad"] = "smp_stage_quad: " + distribs;
    t["mmap_mark"] =
        "mmap_mark: the C++ port mmap_mark(Map base, weights n x C) is a different function from "
        "the MATLAB reference mmap_mark(MMAP, prob K x R): it marks the arrivals of a MAP by "
        "phase-dependent weights, whereas the reference re-marks the K classes of an MMAP into R "
        "classes. Neither signature can be served from the other, so the MATLAB name is not "
        "exposed rather than returning a different quantity under it.";
    t["map_transform"] =
        "map_transform: the C++ library has no function of this name; it is the header "
        "map_transform.h, whose functions (map_normalize, map_scale, map_erlang, map_sum, "
        "map_sumind, map_mixture, map_renewal, map_stochcomp, map_hyperexp, ...) are exposed "
        "under their own names.";
    t["qbd_r"] =
        "qbd_r: the C++ library has no function of this name; it is the header qbd_r.h, whose "
        "functions (qbd_R, qbd_R_logred, qbd_fundmat, qbd_caudal, qbd_pi, qbd_R_residual, "
        "qbd_G_residual) are exposed under their own names.";
}
 
}  // namespace reg
}  // namespace line