api_dispatch_mc.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 'mc' group (see api_dispatch_detail.h): the
 * Markov-chain functions not already in api_dispatch.cpp, the npfqn traffic and
 * bound functions, the discrete-time cycle constants dpfqn_nc / dpfqn_ncld and
 * the stochastic network calculus envelopes snc_*.
 *
 * CONVENTIONS, shared with the mc entries of api_dispatch.cpp:
 *  - a state, station or stage INDEX is 0-based on the way in and is reported
 *    unchanged on the way out (an index is not a measurement). The one
 *    exception is npfqn_sqd's stationToStateful, which the port documents as
 *    ONE-based because it is sn.stationToStateful copied verbatim;
 *  - a MATLAB struct output (ctx, info, aux, result, out) is a nested object
 *    keyed by the MATLAB field names;
 *  - an MMAP is the MATLAB cell {D0, D1, D1^(1), ..., D1^(R)} as a list of
 *    matrices, and a list of MMAPs is a list of such lists;
 *  - a complex number is a pair [re, im]; a plain number is read as real.
 *
 * WHERE MATLAB TAKES AN sn. npfqn_nonexp_approx and npfqn_sqd take `sn` in
 * MATLAB, but their ports take the handful of arrays the reference unpacks
 * from it as their first act (see the headers). Those arrays are plain data,
 * so the entries read them under the PORT's parameter names and refuse `sn`.
 */
 
#include "api_dispatch_detail.h"
 
#include <cmath>
#include <complex>
#include <cstdint>
#include <cstdlib>
#include <string>
#include <vector>
 
#include "line/api/dpfqn/dpfqn_nc.h"
#include "line/api/mc/ctmc_accumulated.h"
#include "line/api/mc/ctmc_classical.h"
#include "line/api/mc/ctmc_gth.h"
#include "line/api/mc/ctmc_ssd.h"
#include "line/api/mc/ctmc_isfeasible.h"
#include "line/api/mc/ctmc_memory_gate.h"
#include "line/api/mc/ctmc_passage.h"
#include "line/api/mc/ctmc_rand.h"
#include "line/api/mc/ctmc_simulate.h"
#include "line/api/mc/ctmc_uniformization.h"
#include "line/api/mc/dtmc_solve.h"
#include "line/api/mc/dtmc_transient.h"
#include "line/api/npfqn/npfqn_bnd_bgt.h"
#include "line/api/npfqn/npfqn_bnd_bpt.h"
#include "line/api/npfqn/npfqn_dps_morrison.h"
#include "line/api/npfqn/npfqn_feedback_elim.h"
#include "line/api/npfqn/npfqn_nonexp_approx.h"
#include "line/api/npfqn/npfqn_rqna_weight.h"
#include "line/api/npfqn/npfqn_sqd.h"
#include "line/api/npfqn/npfqn_traffic_idc.h"
#include "line/api/npfqn/npfqn_traffic_merge.h"
#include "line/api/npfqn/npfqn_traffic_merge_cs.h"
#include "line/api/npfqn/npfqn_traffic_rqt.h"
#include "line/api/npfqn/npfqn_traffic_split_cs.h"
#include "line/api/snc/snc_conv.h"
#include "line/api/snc/snc_env_cpoisson.h"
#include "line/api/snc/snc_env_map.h"
#include "line/api/snc/snc_env_poisson.h"
#include "line/api/snc/snc_env_tokenbucket.h"
#include "line/api/snc/snc_leftover.h"
#include "line/api/snc/snc_output.h"
#include "line/api/snc/snc_srv_exp.h"
#include "line/api/snc/snc_srv_rate.h"
 
namespace line {
namespace reg {
namespace {
 
// ---------------------------------------------------------------------------
// Shared readers and encoders
// ---------------------------------------------------------------------------
 
/** One non-negative index from a JSON number or numeric string. */
std::size_t index_from_json(const Json& j, const std::string& where) {
    const std::string text = decimal_text(j, where);
    char* end = nullptr;
    const long long v = std::strtoll(text.c_str(), &end, 10);
    if (end == text.c_str() || *end != '\0')
        throw InputError(where + ": an index must be an integer, got " + text);
    if (v < 0) throw InputError(where + ": an index cannot be negative, got " + text);
    return static_cast<std::size_t>(v);
}
 
/** A flat list of 0-based indices; a bare number is a one-element list. */
std::vector<std::size_t> indices_from_json(const Json& j, const std::string& where) {
    std::vector<std::size_t> out;
    if (!j.is_array()) {
        out.push_back(index_from_json(j, where));
        return out;
    }
    for (std::size_t k = 0; k < j.size(); ++k)
        out.push_back(index_from_json(j[k], where + "[" + std::to_string(k) + "]"));
    return out;
}
 
std::vector<std::size_t> read_indices(Args& a, const char* key) {
    return indices_from_json(a.get(key), a.get(key).is_null() ? key : std::string(key));
}
 
/** A list of index lists (sigma of npfqn_bnd_bgt). */
std::vector<std::vector<std::size_t> > read_index_lists(Args& a, const char* key) {
    const Json& j = a.get(key);
    const std::string where = std::string(key);
    if (!j.is_array()) throw InputError(where + ": expected a list of index lists");
    std::vector<std::vector<std::size_t> > out;
    for (std::size_t b = 0; b < j.size(); ++b)
        out.push_back(indices_from_json(j[b], where + "[" + std::to_string(b) + "]"));
    return out;
}
 
Json encode_indices(const std::vector<std::size_t>& v) {
    Json a = Json::array();
    for (std::size_t i = 0; i < v.size(); ++i) a.push_back(static_cast<std::uint64_t>(v[i]));
    return a;
}
 
/** A per-station flag vector: an array of booleans or of 0/1. */
std::vector<bool> read_flags(Args& a, const char* key) {
    const Json& j = a.get(key);
    const std::string where = std::string(key);
    std::vector<bool> out;
    const Json flat = j.is_array() ? j : Json::array({j});
    for (std::size_t k = 0; k < flat.size(); ++k) {
        const Json& e = flat[k];
        if (e.is_boolean()) {
            out.push_back(e.get<bool>());
        } else {
            const std::string text = decimal_text(e, where);
            if (text != "0" && text != "1")
                throw InputError(where + ": a flag must be true, false, 0 or 1, got " + text);
            out.push_back(text == "1");
        }
    }
    return out;
}
 
/** A list of numeric vectors (mu of npfqn_bnd_bgt). */
template <class T>
std::vector<std::vector<T> > read_vector_list(Args& a, const char* key) {
    const Json& j = a.get(key);
    const std::string where = std::string(key);
    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;
}
 
/** A matrix as the list of its rows (P of dpfqn_ncld, J x N). */
template <class T>
std::vector<std::vector<T> > rows_of(const Matrix<T>& m) {
    std::vector<std::vector<T> > out(m.rows(), std::vector<T>(m.cols()));
    for (std::size_t i = 0; i < m.rows(); ++i)
        for (std::size_t k = 0; k < m.cols(); ++k) out[i][k] = m(i, k);
    return out;
}
 
/** An MMAP cell {D0, D1, D1^(1), ...} as the mam::Mmap struct; [] is the empty flow. */
template <class T>
mam::Mmap<T> mmap_from_json(const Json& j, const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected an MMAP, a list of matrices");
    mam::Mmap<T> m;
    if (j.empty()) return m;
    if (j.size() < 2)
        throw InputError(where + ": an MMAP is {D0, D1, D1^(1), ...}, got one matrix");
    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) + "]"));
    return m;
}
 
template <class T>
std::vector<mam::Mmap<T> > read_mmap_list(Args& a, const char* key) {
    const Json& j = a.get(key);
    const std::string where = std::string(key);
    if (!j.is_array()) throw InputError(where + ": expected a list of MMAPs");
    std::vector<mam::Mmap<T> > out;
    for (std::size_t i = 0; i < j.size(); ++i)
        out.push_back(mmap_from_json<T>(j[i], where + "[" + std::to_string(i) + "]"));
    return out;
}
 
template <class T>
Json encode_mmap(const mam::Mmap<T>& m) {
    std::vector<Matrix<T> > cell;
    cell.push_back(m.D0);
    cell.push_back(m.D1);
    for (const Matrix<T>& d : m.Dc) cell.push_back(d);
    return encode_matrices(cell);
}
 
/** Complex points: a number is real, a pair [re, im] is complex. */
std::vector<std::complex<double> > read_complex(Args& a, const char* key) {
    const Json& j = a.get(key);
    const std::string where = std::string(key);
    const Json flat = j.is_array() ? j : Json::array({j});
    std::vector<std::complex<double> > out;
    for (std::size_t k = 0; k < flat.size(); ++k) {
        const Json& e = flat[k];
        const std::string w = where + "[" + std::to_string(k) + "]";
        if (e.is_array()) {
            if (e.size() != 2) throw InputError(w + ": a complex number is a pair [re, im]");
            out.emplace_back(number_from_json<double>(e[0], w), number_from_json<double>(e[1], w));
        } else {
            out.emplace_back(number_from_json<double>(e, w), 0.0);
        }
    }
    return out;
}
 
Json encode_complex(const std::vector<std::complex<double> >& v) {
    Json a = Json::array();
    for (const std::complex<double>& z : v) a.push_back(Json::array({z.real(), z.imag()}));
    return a;
}
 
/** A nested struct argument (config, options, corrections), or null when absent. */
const Json* sub_object(Args& a, const char* key) {
    const Json* j = a.opt(key);
    if (j && !j->is_object())
        throw InputError(std::string(key) + ": expected an object with named fields");
    return j;
}
 
/** The C++ default arguments of options.method / options.lti_method. */
struct PassageOptions {
    bool hasMethod = false;
    std::string method;
    bool hasLti = false;
    std::string lti;
};
 
/**
 * MATLAB's `options` of ctmc_passage_time / ctmc_passage_taboo. lti_method
 * picks the Laplace inversion of method 'lt' and has no effect under any other
 * method, so naming it without method 'lt' is refused rather than ignored.
 */
PassageOptions read_passage_options(Args& a, const std::string& fn) {
    PassageOptions o;
    const Json* j = sub_object(a, "options");
    if (!j) return o;
    Args op(*j, fn + ": options");
    o.hasMethod = op.has("method");
    if (o.hasMethod) o.method = op.text("method", "");
    o.hasLti = op.has("lti_method");
    if (o.hasLti) o.lti = op.text("lti_method", "");
    op.done();
    if (o.hasLti && !(o.hasMethod && o.method == "lt"))
        throw InputError(fn + ": options.lti_method selects the inversion of method 'lt' and "
                         "has no effect otherwise; set options.method to 'lt' or drop it");
    return o;
}
 
// ---------------------------------------------------------------------------
// mc
// ---------------------------------------------------------------------------
 
/** [S,Q11,Q12,Q21,Q22,T,X] = ctmc_stochcomp(Q, I); I is 0-based and optional. */
struct OpCtmcStochcomp {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const bool hasI = a.has("I");
        const std::vector<std::size_t> I = hasI ? read_indices(a, "I") : std::vector<std::size_t>();
        a.done();
        const mc::StochCompResult<T> r = hasI ? mc::ctmc_stochcomp(Q, I) : mc::ctmc_stochcomp(Q);
        Json out;
        out["S"] = encode_matrix(r.S);
        out["Q11"] = encode_matrix(r.Q11);
        out["Q12"] = encode_matrix(r.Q12);
        out["Q21"] = encode_matrix(r.Q21);
        out["Q22"] = encode_matrix(r.Q22);
        out["T"] = encode_matrix(r.T12);
        out["X"] = encode_matrix(r.X);
        return out;
    }
};
 
/**
 * tol / maxiter of ctmc_uniformization and ctmc_timeaverage. They are
 * positional in MATLAB, so maxiter cannot be named without tol there; the same
 * holds here, which keeps every call one of the port's three arities.
 */
template <class R, class F>
R call_tol_maxiter(Args& a, const std::string& fn, F f) {
    const bool hasTol = a.has("tol"), hasMax = a.has("maxiter");
    if (hasMax && !hasTol)
        throw InputError(fn + ": 'maxiter' follows 'tol' in the MATLAB signature; name both");
    const double tol = hasTol ? a.number<double>("tol") : 0.0;
    const long maxiter = hasMax ? static_cast<long>(a.required_integer("maxiter")) : 0;
    a.done();
    if (hasMax) return f(3, tol, maxiter);
    if (hasTol) return f(2, tol, maxiter);
    return f(1, tol, maxiter);
}
 
/** [pi,kmax] = ctmc_uniformization(pi0,Q,t,tol,maxiter). */
struct OpCtmcUniformization {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> pi0 = a.vector<T>("pi0");
        const Matrix<T> Q = a.matrix<T>("Q");
        const T t = a.number<T>("t");
        const mc::UniformizationResult<T> r = call_tol_maxiter<mc::UniformizationResult<T> >(
            a, "ctmc_uniformization", [&](int n, double tol, long maxiter) {
                if (n == 3) return mc::ctmc_uniformization(pi0, Q, t, tol, maxiter);
                if (n == 2) return mc::ctmc_uniformization(pi0, Q, t, tol);
                return mc::ctmc_uniformization(pi0, Q, t);
            });
        Json out;
        out["pi"] = encode_vector(r.pi);
        out["kmax"] = encode_count(r.kmax);
        return out;
    }
};
 
/** [piTimeAvg,piExit,kmax] = ctmc_timeaverage(pi0,Q,t,tol,maxiter). */
struct OpCtmcTimeaverage {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> pi0 = a.vector<T>("pi0");
        const Matrix<T> Q = a.matrix<T>("Q");
        const T t = a.number<T>("t");
        const mc::TimeAverageResult<T> r = call_tol_maxiter<mc::TimeAverageResult<T> >(
            a, "ctmc_timeaverage", [&](int n, double tol, long maxiter) {
                if (n == 3) return mc::ctmc_timeaverage(pi0, Q, t, tol, maxiter);
                if (n == 2) return mc::ctmc_timeaverage(pi0, Q, t, tol);
                return mc::ctmc_timeaverage(pi0, Q, t);
            });
        Json out;
        out["piTimeAvg"] = encode_vector(r.piTimeAvg);
        out["piExit"] = encode_vector(r.piExit);
        out["kmax"] = encode_count(r.kmax);
        return out;
    }
};
 
/**
 * [p,flag,relres,iter] = ctmc_power / ctmc_jacobi / ctmc_gaussseidel(Q,tol,maxiter,p0)
 * and [p,flag,relres,iter,omega] = ctmc_sor(Q,tol,maxiter,omega,p0). Every
 * argument but Q is optional and named, a nonpositive tol / maxiter / omega
 * selecting the default (omega <= 0 is the adaptive SOR rule).
 */
template <mc::detail::SweepKind K>
struct OpCtmcClassical {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const double tol = a.scalar<double>("tol", mc::CLASSICAL_DEFAULT_TOL);
        const long maxiter = static_cast<long>(a.integer("maxiter", static_cast<int>(mc::CLASSICAL_DEFAULT_MAXITER)));
        const double omega = (K == mc::detail::SweepKind::Sor) ? a.scalar<double>("omega", 0.0)
                             : (K == mc::detail::SweepKind::GaussSeidel) ? 1.0 : 0.0;
        const std::vector<T> p0 = a.vector_or_empty<T>("p0");
        a.done();
        const mc::ClassicalResult<T> r = mc::detail::stationary_sweep_iterate(K, Q, tol, maxiter, omega, p0);
        Json out;
        out["p"] = encode_vector(r.p);
        out["flag"] = r.flag;
        out["relres"] = encode_scalar(r.relres);
        out["iter"] = static_cast<std::int64_t>(r.iter);
        if (K == mc::detail::SweepKind::Sor) out["omega"] = encode_scalar(r.omega);
        return out;
    }
};
 
/**
 * [p,flag] = ctmc_gth(Q) / dtmc_gth(P): GTH elimination at every arithmetic
 * (no subtraction, so exact at --arith exact). With several closed classes p is
 * NaN and flag 1, as the MATLAB two-output form returns.
 */
template <bool Dtmc>
struct OpGth {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>(Dtmc ? "P" : "Q");
        a.done();
        const mc::GthResult<T> r = mc::ctmc_gth_flag(Q);
        Json out;
        if (r.flag == 0) {
            out["p"] = encode_vector(r.p);
        } else {
            Json nan = Json::array();
            for (std::size_t i = 0; i < Q.rows(); ++i) nan.push_back("NaN");
            out["p"] = nan;
        }
        out["flag"] = r.flag;
        return out;
    }
};
 
/** J = ctmc_cesaro(Q) / dtmc_cesaro(P): limit matrix of any finite chain. */
template <bool Dtmc>
struct OpCesaro {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>(Dtmc ? "P" : "Q");
        a.done();
        Json out;
        out["J"] = encode_matrix(mc::ctmc_cesaro(Q));
        return out;
    }
};
 
/** p = ctmc_ssd(C,R) / dtmc_ssd(C,R); an absent edge is R = "Inf" or C = 0. */
struct OpSsd {
    static Json run(Args& a) {
        const Matrix<double> C = a.matrix<double>("C");
        const Matrix<double> R = a.matrix<double>("R");
        a.done();
        Json out;
        out["p"] = encode_vector(mc::ctmc_ssd(C, R));
        return out;
    }
};
 
/** [L,piT,kmax] = ctmc_accumulated(pi0,Q,t,tol,maxiter); t may be a vector. */
struct OpCtmcAccumulated {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> pi0 = a.vector<T>("pi0");
        const Matrix<T> Q = a.matrix<T>("Q");
        const std::vector<T> t = a.vector<T>("t");
        const mc::AccumulatedResult<T> r = call_tol_maxiter<mc::AccumulatedResult<T> >(
            a, "ctmc_accumulated", [&](int n, double tol, long maxiter) {
                if (n == 3) return mc::ctmc_accumulated(pi0, Q, t, tol, maxiter);
                if (n == 2) return mc::ctmc_accumulated(pi0, Q, t, tol);
                return mc::ctmc_accumulated(pi0, Q, t);
            });
        Json out;
        out["L"] = encode_matrix(r.L);
        out["piT"] = encode_matrix(r.piT);
        out["kmax"] = encode_count(r.kmax);
        return out;
    }
};
 
/**
 * [ok,msg] = ctmc_memory_gate(logNstates,force,verbose,safetyFraction). The
 * kernel is double-valued whatever the arithmetic; 'verbose' is refused because
 * it prints the calibration on stdout, the stream that carries the JSON reply.
 */
struct OpCtmcMemoryGate {
    template <class T>
    static Json run(Args& a) {
        const double logN = a.number<double>("logNstates");
        a.unsupported("verbose",
                      "ctmc_memory_gate: 'verbose' prints to stdout, which carries the JSON "
                      "reply; the decision and its message are returned as ok and msg");
        const bool hasForce = a.has("force"), hasSafety = a.has("safetyFraction");
        const bool force = hasForce ? a.boolean("force", false) : false;
        const double safety = hasSafety ? a.number<double>("safetyFraction") : 0.0;
        a.done();
        mc::CtmcGateResult r;
        if (hasSafety)
            r = mc::ctmc_memory_gate(logN, force, safety);
        else if (hasForce)
            r = mc::ctmc_memory_gate(logN, force);
        else
            r = mc::ctmc_memory_gate(logN);
        Json out;
        out["ok"] = r.ok;
        out["msg"] = r.message;
        return out;
    }
};
 
/** bool = ctmc_isfeasible(Q, tol). */
struct OpCtmcIsfeasible {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const bool hasTol = a.has("tol");
        const T tol = hasTol ? a.number<T>("tol") : T();
        a.done();
        Json out;
        out["bool"] = hasTol ? mc::ctmc_isfeasible(Q, tol) : mc::ctmc_isfeasible(Q);
        return out;
    }
};
 
/** h = ctmc_hitting_time(Q, targetStates) and h = dtmc_hitting_time(P, targetStates). */
struct OpCtmcHittingTime {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const std::vector<std::size_t> target = read_indices(a, "targetStates");
        a.done();
        Json out;
        out["h"] = encode_vector(mc::ctmc_hitting_time(Q, target));
        return out;
    }
};
 
struct OpDtmcHittingTime {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> P = a.matrix<T>("P");
        const std::vector<std::size_t> target = read_indices(a, "targetStates");
        a.done();
        Json out;
        out["h"] = encode_vector(mc::dtmc_hitting_time(P, target));
        return out;
    }
};
 
/**
 * pi_t = dtmc_transient(P, pi0, steps). An absent pi0 is the uniform law (the
 * port's empty-vector form). The port has no (P, pi0) overload, so a call that
 * names pi0 but not steps passes MATLAB's steps = 1 explicitly.
 */
struct OpDtmcTransient {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> P = a.matrix<T>("P");
        const bool hasPi0 = a.has("pi0"), hasSteps = a.has("steps");
        const std::vector<T> pi0 = a.vector_or_empty<T>("pi0");
        const std::size_t steps = hasSteps ? a.count("steps") : 1;
        a.done();
        Json out;
        out["pi_t"] = encode_matrix((hasPi0 || hasSteps) ? mc::dtmc_transient(P, pi0, steps)
                                                         : mc::dtmc_transient(P));
        return out;
    }
};
 
/** [alpha,S,s0,keep,atom] = ctmc_passage_ph(Q,pi0,target); an absent pi0 is the stationary law. */
struct OpCtmcPassagePh {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const std::vector<T> pi0 = a.vector_or_empty<T>("pi0");
        const std::vector<std::size_t> target = read_indices(a, "target");
        a.done();
        const mc::PassagePh<T> r = mc::ctmc_passage_ph(Q, pi0, target);
        Json out;
        out["alpha"] = encode_vector(r.alpha);
        out["S"] = encode_matrix(r.S);
        out["s0"] = encode_vector(r.s0);
        out["keep"] = encode_indices(r.keep);
        out["atom"] = encode_scalar(r.atom);
        return out;
    }
};
 
/** L = ctmc_passage_lst(Q,pi0,target,s); s and L are lists of [re, im]. */
struct OpCtmcPassageLst {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const std::vector<T> pi0 = a.vector_or_empty<T>("pi0");
        const std::vector<std::size_t> target = read_indices(a, "target");
        const std::vector<std::complex<double> > s = read_complex(a, "s");
        a.done();
        Json out;
        out["L"] = encode_complex(mc::ctmc_passage_lst(Q, pi0, target, s));
        return out;
    }
};
 
/** [mall,m] = ctmc_passage_moments(Q,pi0,target,nmax). */
struct OpCtmcPassageMoments {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const std::vector<T> pi0 = a.vector_or_empty<T>("pi0");
        const std::vector<std::size_t> target = read_indices(a, "target");
        const bool hasN = a.has("nmax");
        const std::size_t nmax = hasN ? a.count("nmax") : 0;
        a.done();
        const mc::PassageMoments<T> r = hasN ? mc::ctmc_passage_moments(Q, pi0, target, nmax)
                                             : mc::ctmc_passage_moments(Q, pi0, target);
        Json out;
        out["mall"] = encode_matrix(r.mall);
        out["m"] = encode_vector(r.m);
        return out;
    }
};
 
/**
 * [F,f,out] = ctmc_passage_time(Q,pi0,target,tset,options). out carries atom
 * and the phase-type representation (alpha, S, s0, keep) the curve is computed
 * from, which is ctmc_passage_ph of the same arguments.
 */
struct OpCtmcPassageTime {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const std::vector<T> pi0 = a.vector_or_empty<T>("pi0");
        const std::vector<std::size_t> target = read_indices(a, "target");
        const std::vector<double> tset = a.vector<double>("tset");
        const PassageOptions o = read_passage_options(a, "ctmc_passage_time");
        a.done();
        const mc::PassageCurve<T> c =
            o.hasLti      ? mc::ctmc_passage_time(Q, pi0, target, tset, o.method, o.lti)
            : o.hasMethod ? mc::ctmc_passage_time(Q, pi0, target, tset, o.method)
                          : mc::ctmc_passage_time(Q, pi0, target, tset);
        const mc::PassagePh<T> ph = mc::ctmc_passage_ph(Q, pi0, target);
        Json out;
        out["F"] = encode_vector(c.F);
        out["f"] = encode_vector(c.f);
        Json o2;
        o2["atom"] = encode_scalar(c.atom);
        o2["alpha"] = encode_vector(ph.alpha);
        o2["S"] = encode_matrix(ph.S);
        o2["s0"] = encode_vector(ph.s0);
        o2["keep"] = encode_indices(ph.keep);
        out["out"] = o2;
        return out;
    }
};
 
/** [F,f,out] = ctmc_passage_taboo(Q,pi0,target,taboo,tset,options). */
struct OpCtmcPassageTaboo {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const std::vector<T> pi0 = a.vector_or_empty<T>("pi0");
        const std::vector<std::size_t> target = read_indices(a, "target");
        const std::vector<std::size_t> taboo = read_indices(a, "taboo");
        const bool hasT = a.has("tset");
        const std::vector<double> tset = a.vector_or_empty<double>("tset");
        const PassageOptions o = read_passage_options(a, "ctmc_passage_taboo");
        a.done();
        const mc::PassageTaboo<T> r =
            o.hasLti      ? mc::ctmc_passage_taboo(Q, pi0, target, taboo, tset, o.method, o.lti)
            : o.hasMethod ? mc::ctmc_passage_taboo(Q, pi0, target, taboo, tset, o.method)
            : hasT        ? mc::ctmc_passage_taboo(Q, pi0, target, taboo, tset)
                          : mc::ctmc_passage_taboo(Q, pi0, target, taboo);
        Json out;
        out["F"] = encode_vector(r.F);
        out["f"] = encode_vector(r.f);
        Json o2;
        o2["atom"] = encode_scalar(r.atom);
        o2["reachprob"] = encode_scalar(r.reachprob);
        o2["tabooprob"] = encode_scalar(r.tabooprob);
        o2["pB"] = encode_vector(r.pB);
        o2["pC"] = encode_vector(r.pC);
        o2["alpha"] = encode_vector(r.alpha);
        o2["S"] = encode_matrix(r.S);
        o2["sB"] = encode_vector(r.sB);
        o2["sC"] = encode_vector(r.sC);
        o2["keep"] = encode_indices(r.keep);
        o2["target"] = encode_indices(r.target);
        o2["taboo"] = encode_indices(r.taboo);
        o2["method"] = r.method;
        out["out"] = o2;
        return out;
    }
};
 
/** [mall,m,pB] = ctmc_passage_taboo_moments(Q,pi0,target,taboo,nmax). */
struct OpCtmcPassageTabooMoments {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const std::vector<T> pi0 = a.vector_or_empty<T>("pi0");
        const std::vector<std::size_t> target = read_indices(a, "target");
        const std::vector<std::size_t> taboo = read_indices(a, "taboo");
        const bool hasN = a.has("nmax");
        const std::size_t nmax = hasN ? a.count("nmax") : 0;
        a.done();
        const mc::PassageTabooMoments<T> r =
            hasN ? mc::ctmc_passage_taboo_moments(Q, pi0, target, taboo, nmax)
                 : mc::ctmc_passage_taboo_moments(Q, pi0, target, taboo);
        Json out;
        out["mall"] = encode_matrix(r.mall);
        out["m"] = encode_vector(r.m);
        out["pB"] = encode_vector(r.pB);
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// dpfqn
// ---------------------------------------------------------------------------
 
/** [lG,G,G1] = dpfqn_nc(p, N). */
struct OpDpfqnNc {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> p = a.vector<T>("p");
        const std::size_t N = a.count("N");
        a.done();
        const dpfqn::DtNcResult<T> r = dpfqn::dpfqn_nc(p, N);
        Json out;
        out["lG"] = encode_scalar(r.lG);
        out["G"] = encode_scalar(r.G);
        out["G1"] = encode_vector(r.G1);
        return out;
    }
};
 
/** [lG,G,W,Gc,Wa] = dpfqn_ncld(P, N), P a J-by-N matrix with P(j,n) = p_j(n). */
struct OpDpfqnNcld {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> P = a.matrix<T>("P");
        const std::size_t N = a.count("N");
        a.done();
        const dpfqn::DtNcLdResult<T> r = dpfqn::dpfqn_ncld(rows_of(P), N);
        Json out;
        out["lG"] = encode_scalar(r.lG);
        out["G"] = encode_vector(r.G);
        out["W"] = encode_vectors(r.W);
        out["Gc"] = encode_vectors(r.Gc);
        out["Wa"] = encode_vectors(r.Wa);
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// npfqn
// ---------------------------------------------------------------------------
 
/**
 * ctx = npfqn_traffic_idc(lambda0,P,c2a0,a0IdcFun,mu,cs2,sIdcFun,corrections).
 * The two IDC callables are stored in MATLAB's ctx only for npfqn_traffic_idc_at
 * (library-only) and play no part in ctx itself, so they are refused by name.
 */
struct OpNpfqnTrafficIdc {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda0 = a.vector<T>("lambda0");
        const Matrix<T> P = a.matrix<T>("P");
        const std::vector<T> c2a0 = a.vector<T>("c2a0");
        const std::vector<T> mu = a.vector<T>("mu");
        const std::vector<T> cs2 = a.vector<T>("cs2");
        a.unsupported("a0IdcFun", "npfqn_traffic_idc: a0IdcFun is a function handle with no "
                                  "JSON form; it is used only by npfqn_traffic_idc_at");
        a.unsupported("sIdcFun", "npfqn_traffic_idc: sIdcFun is a function handle with no "
                                 "JSON form; it is used only by npfqn_traffic_idc_at");
        const Json* cj = sub_object(a, "corrections");
        npfqn::TrafficIdcCorrections corr;
        if (cj) {
            Args c(*cj, "npfqn_traffic_idc: corrections");
            if (c.has("alpha")) corr.alpha = c.boolean("alpha", corr.alpha);
            if (c.has("beta")) corr.beta = c.boolean("beta", corr.beta);
            c.done();
        }
        a.done();
        const npfqn::TrafficIdcContext<T> x = cj ? npfqn::npfqn_traffic_idc(lambda0, P, c2a0, mu, cs2, corr)
                                                 : npfqn::npfqn_traffic_idc(lambda0, P, c2a0, mu, cs2);
        Json ctx;
        ctx["K"] = encode_count(x.K);
        ctx["lambda"] = encode_vector(x.lambda);
        ctx["lambda0"] = encode_vector(x.lambda0);
        ctx["mu"] = encode_vector(x.mu);
        ctx["rho"] = encode_vector(x.rho);
        ctx["cs2"] = encode_vector(x.cs2);
        ctx["c2a0"] = encode_vector(x.c2a0);
        ctx["c2a"] = encode_vector(x.c2a);
        ctx["c2d"] = encode_vector(x.c2d);
        ctx["c2x"] = encode_vector(x.c2x);
        ctx["P"] = encode_matrix(x.P);
        ctx["Xi"] = encode_matrix(x.Xi);
        ctx["lam_ji"] = encode_matrix(x.lam_ji);
        ctx["c2aij"] = encode_matrix(x.c2aij);
        ctx["c2alpha"] = encode_matrix(x.c2alpha);
        ctx["zetaAll"] = encode_matrices(x.zetaAll);
        Json out;
        out["ctx"] = ctx;
        return out;
    }
};
 
/**
 * SMMAP = npfqn_traffic_split_cs(MMAP, P, config). The reference never reads
 * config, so naming it is refused rather than silently dropped.
 */
struct OpNpfqnTrafficSplitCs {
    template <class T>
    static Json run(Args& a) {
        const std::vector<Matrix<T> > M = a.matrices<T>("MMAP");
        const Matrix<T> P = a.matrix<T>("P");
        a.unsupported("config", "npfqn_traffic_split_cs: the reference accepts config but never "
                                "reads it; there is nothing for it to select");
        a.done();
        const std::vector<npfqn::Mmap<T> > s = npfqn::npfqn_traffic_split_cs(M, P);
        Json out = Json::array();
        for (const npfqn::Mmap<T>& m : s) out.push_back(encode_matrices(m));
        Json r;
        r["SMMAP"] = out;
        return r;
    }
};
 
npfqn::Merge merge_from_text(const std::string& s, const std::string& fn) {
    if (s == "default") return npfqn::Merge::Default;
    if (s == "super") return npfqn::Merge::Super;
    if (s == "mixture") return npfqn::Merge::Mixture;
    if (s == "interpos") return npfqn::Merge::Interpos;
    throw InputError(fn + ": config.merge must be 'default', 'super', 'mixture' or 'interpos', got '" +
                     s + "'");
}
 
npfqn::Compress compress_from_text(const std::string& s, const std::string& fn) {
    if (s == "default") return npfqn::Compress::Default;
    if (s == "none") return npfqn::Compress::None;
    throw InputError(fn + ": config.compress must be 'default' or 'none', got '" + s + "'");
}
 
/** SMMAP = npfqn_traffic_merge_cs(MMAP, prob, config); config.merge only. */
struct OpNpfqnTrafficMergeCs {
    template <class T>
    static Json run(Args& a) {
        const std::vector<mam::Mmap<T> > flows = read_mmap_list<T>(a, "MMAP");
        const Matrix<T> prob = a.matrix<T>("prob");
        const Json* cj = sub_object(a, "config");
        bool hasMerge = false;
        npfqn::Merge merge = npfqn::Merge::Default;
        if (cj) {
            Args c(*cj, "npfqn_traffic_merge_cs: config");
            hasMerge = c.has("merge");
            if (hasMerge) merge = merge_from_text(c.text("merge", ""), "npfqn_traffic_merge_cs");
            c.done();
        }
        a.done();
        const mam::Mmap<T> s = hasMerge ? npfqn::npfqn_traffic_merge_cs(flows, prob, merge)
                                        : npfqn::npfqn_traffic_merge_cs(flows, prob);
        Json out;
        out["SMMAP"] = encode_mmap(s);
        return out;
    }
};
 
/** SMMAP = npfqn_traffic_merge(MMAP, config); config.merge and config.compress. */
struct OpNpfqnTrafficMerge {
    template <class T>
    static Json run(Args& a) {
        const std::vector<mam::Mmap<T> > flows = read_mmap_list<T>(a, "MMAP");
        const Json* cj = sub_object(a, "config");
        npfqn::MergeConfig cfg;
        if (cj) {
            Args c(*cj, "npfqn_traffic_merge: config");
            if (c.has("merge")) cfg.merge = merge_from_text(c.text("merge", ""), "npfqn_traffic_merge");
            if (c.has("compress"))
                cfg.compress = compress_from_text(c.text("compress", ""), "npfqn_traffic_merge");
            c.done();
        }
        a.done();
        const mam::Mmap<T> s = cj ? npfqn::npfqn_traffic_merge(flows, cfg) : npfqn::npfqn_traffic_merge(flows);
        Json out;
        out["SMMAP"] = encode_mmap(s);
        return out;
    }
};
 
/** w = npfqn_rqna_weight(t): a scalar t gives a scalar w, an array an array. */
struct OpNpfqnRqnaWeight {
    template <class T>
    static Json run(Args& a) {
        const bool vec = a.get("t").is_array();
        Json out;
        if (vec) {
            const std::vector<T> t = a.vector<T>("t");
            a.done();
            out["w"] = encode_vector(npfqn::npfqn_rqna_weight(t));
        } else {
            const T t = a.number<T>("t");
            a.done();
            out["w"] = encode_scalar(npfqn::npfqn_rqna_weight(t));
        }
        return out;
    }
};
 
/** [lambda,Gamma,alpha] = npfqn_traffic_rqt(lambda0,Gamma0,alpha0,F). */
struct OpNpfqnTrafficRqt {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda0 = a.vector<T>("lambda0");
        const std::vector<T> Gamma0 = a.vector<T>("Gamma0");
        const std::vector<T> alpha0 = a.vector<T>("alpha0");
        const Matrix<T> F = a.matrix<T>("F");
        a.done();
        const npfqn::TrafficRqt<T> r = npfqn::npfqn_traffic_rqt(lambda0, Gamma0, alpha0, F);
        Json out;
        out["lambda"] = encode_vector(r.lambda);
        out["Gamma"] = encode_vector(r.Gamma);
        out["alpha"] = encode_vector(r.alpha);
        return out;
    }
};
 
/**
 * [zlb,x,info] = npfqn_bnd_bpt(lambda0,mu,P,stationOf,c); stationOf is 0-based
 * and an absent c is MATLAB's empty c (all ones).
 */
struct OpNpfqnBndBpt {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda0 = a.vector<T>("lambda0");
        const std::vector<T> mu = a.vector<T>("mu");
        const Matrix<T> P = a.matrix<T>("P");
        const std::vector<std::size_t> stationOf = read_indices(a, "stationOf");
        const std::vector<T> c = a.vector_or_empty<T>("c");
        a.done();
        const npfqn::BndBpt<T> r = npfqn::npfqn_bnd_bpt(lambda0, mu, P, stationOf, c);
        Json out;
        out["zlb"] = encode_scalar(r.zlb);
        out["x"] = encode_vector(r.x);
        Json info;
        info["lambda"] = encode_vector(r.lambda);
        info["rho"] = encode_vector(r.rho);
        info["rhoStation"] = encode_vector(r.rhoStation);
        info["nvars"] = encode_count(r.nvars);
        info["nrows"] = encode_count(r.nrows);
        out["info"] = info;
        return out;
    }
};
 
/**
 * [Qub,info] = npfqn_bnd_bgt(lambda,mu,sigma,J). mu and sigma are lists with
 * one entry per type; sigma holds 0-based stations. An absent J is derived from
 * sigma as MATLAB derives it (one past the largest station index); the port has
 * no shorter overload. info.finite marks where Qub is a bound at all (MATLAB
 * writes Inf there, which JSON cannot carry).
 */
struct OpNpfqnBndBgt {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const std::vector<std::vector<T> > mu = read_vector_list<T>(a, "mu");
        const std::vector<std::vector<std::size_t> > sigma = read_index_lists(a, "sigma");
        std::size_t J = 0;
        if (a.has("J")) {
            J = a.count("J");
        } else {
            for (const std::vector<std::size_t>& s : sigma)
                for (std::size_t j : s) J = std::max(J, j + 1);
        }
        a.done();
        const npfqn::BndBgt<T> r = npfqn::npfqn_bnd_bgt(lambda, mu, sigma, J);
        Json out;
        out["Qub"] = encode_vectors(r.Qub);
        Json info;
        Json fin = Json::array();
        for (bool b : r.finite) fin.push_back(b);
        info["finite"] = fin;
        info["gamma"] = encode_scalar(r.gamma);
        info["Lmax"] = encode_scalar(r.Lmax);
        info["L"] = encode_vectors(r.L);
        info["V"] = encode_vector(r.V);
        info["B"] = encode_scalar(r.B);
        info["U"] = encode_scalar(r.U);
        info["tailRatio"] = encode_scalar(r.tailRatio);
        info["tailStep"] = encode_scalar(r.tailStep);
        info["rho"] = encode_vector(r.rho);
        info["rhoStation"] = encode_vector(r.rhoStation);
        info["scale"] = encode_scalar(r.scale);
        info["classType"] = encode_indices(r.classType);
        info["classStage"] = encode_indices(r.classStage);
        info["classStation"] = encode_indices(r.classStation);
        out["info"] = info;
        return out;
    }
};
 
/**
 * [ST,gamma,nservers,rho,scva,scvs,eta] =
 *     npfqn_nonexp_approx(method,sn,ST,V,SCV,T,U,gamma,nservers).
 * The two fields the reference reads from sn are named as the port names them:
 * isFCFS (sn.sched(i) == FCFS) and rates (sn.rates).
 */
struct OpNpfqnNonexpApprox {
    template <class T>
    static Json run(Args& a) {
        a.unsupported("sn", "npfqn_nonexp_approx: pass the two sn fields the routine reads, "
                            "isFCFS (sn.sched == FCFS per station) and rates (sn.rates)");
        const std::string method = a.text("method", "");
        if (method.empty()) throw InputError("npfqn_nonexp_approx: missing required argument 'method'");
        const std::vector<bool> isFCFS = read_flags(a, "isFCFS");
        const Matrix<T> rates = a.matrix<T>("rates");
        const Matrix<T> ST = a.matrix<T>("ST");
        const Matrix<T> V = a.matrix<T>("V");
        const Matrix<T> SCV = a.matrix<T>("SCV");
        const Matrix<T> Tput = a.matrix<T>("T");
        const Matrix<T> U = a.matrix<T>("U");
        const std::vector<T> gamma = a.vector<T>("gamma");
        const std::vector<T> nservers = a.vector<T>("nservers");
        a.done();
        const npfqn::NonexpApproxResult<T> r =
            npfqn::npfqn_nonexp_approx(method, isFCFS, rates, ST, V, SCV, Tput, U, gamma, nservers);
        Json out;
        out["ST"] = encode_matrix(r.ST);
        out["gamma"] = encode_vector(r.gamma);
        out["nservers"] = encode_vector(r.nservers);
        out["rho"] = encode_vector(r.rho);
        out["scva"] = encode_vector(r.scva);
        out["scvs"] = encode_vector(r.scvs);
        out["eta"] = encode_vector(r.eta);
        return out;
    }
};
 
/**
 * [X,Q,U,R] = npfqn_sqd(sn,N,calibrationMode,serverBlockingTime,neighborMode,
 * v1Policy,initialV1). The sn is replaced by the arrays the port takes (see
 * npfqn_sqd.h): ST, V, cap, isDelay, rt, stationToStateful (ONE-based, as in
 * sn) and nclasses. N is required: MATLAB's fallback reads sn.nclosedjobs.
 */
struct OpNpfqnSqd {
    template <class T>
    static Json run(Args& a) {
        a.unsupported("sn", "npfqn_sqd: pass the arrays the routine unpacks from sn: ST, V, cap, "
                            "isDelay, rt, stationToStateful and nclasses (see npfqn_sqd.h)");
        const std::vector<T> ST = a.vector<T>("ST");
        const std::vector<T> V = a.vector<T>("V");
        const std::vector<T> cap = a.vector<T>("cap");
        const std::vector<bool> isDelay = read_flags(a, "isDelay");
        const Matrix<T> rt = a.matrix<T>("rt");
        const std::vector<std::size_t> s2s = read_indices(a, "stationToStateful");
        const std::size_t nclasses = a.count("nclasses");
        const int N = a.required_integer("N");
        npfqn::SqdOptions<T> opt;
        bool any = false;
        if (a.has("calibrationMode")) {
            opt.calibrationMode = a.required_integer("calibrationMode");
            any = true;
        }
        if (a.has("serverBlockingTime")) {
            opt.serverBlockingTime = a.boolean("serverBlockingTime", opt.serverBlockingTime);
            any = true;
        }
        if (a.has("neighborMode")) {
            const std::string m = a.text("neighborMode", "");
            if (m != "downstream" && m != "ownserver")
                throw InputError("npfqn_sqd: neighborMode must be 'downstream' or 'ownserver', got '" + m + "'");
            opt.ownServerNeighbor = (m == "ownserver");
            any = true;
        }
        if (a.has("v1Policy")) {
            const std::string m = a.text("v1Policy", "");
            if (m != "compound" && m != "fresh")
                throw InputError("npfqn_sqd: v1Policy must be 'compound' or 'fresh', got '" + m + "'");
            opt.freshV1 = (m == "fresh");
            any = true;
        }
        if (a.has("initialV1")) {
            opt.initialV1 = a.vector_or_empty<T>("initialV1");
            any = true;
        }
        a.done();
        const npfqn::SqdResult<T> r =
            any ? npfqn::npfqn_sqd(ST, V, cap, isDelay, rt, s2s, nclasses, N, opt)
                : npfqn::npfqn_sqd(ST, V, cap, isDelay, rt, s2s, nclasses, N);
        Json out;
        out["X"] = encode_vector(r.X);
        out["Q"] = encode_vector(r.Q);
        out["U"] = encode_vector(r.U);
        out["R"] = encode_vector(r.R);
        return out;
    }
};
 
/** [Q,R,X,aux] = npfqn_dps_morrison(N,Z,S,w). */
struct OpNpfqnDpsMorrison {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> N = a.vector<T>("N");
        const std::vector<T> Z = a.vector<T>("Z");
        const std::vector<T> S = a.vector<T>("S");
        const std::vector<T> w = a.vector<T>("w");
        a.done();
        const npfqn::DpsMorrisonResult<T> r = npfqn::npfqn_dps_morrison(N, Z, S, w);
        Json out;
        out["Q"] = encode_vector(r.Q);
        out["R"] = encode_vector(r.R);
        out["X"] = encode_vector(r.X);
        Json aux;
        aux["rho"] = encode_scalar(r.rho);
        aux["a"] = encode_scalar(r.a);
#define LINE_MC_AUX(f) aux[#f] = encode_scalar(r.f);
        LINE_MC_AUX(cB) LINE_MC_AUX(cC) LINE_MC_AUX(cD) LINE_MC_AUX(cH) LINE_MC_AUX(cI)
        LINE_MC_AUX(cJ) LINE_MC_AUX(cK) LINE_MC_AUX(cL) LINE_MC_AUX(cM) LINE_MC_AUX(cQ)
        LINE_MC_AUX(delta) LINE_MC_AUX(cR) LINE_MC_AUX(cS) LINE_MC_AUX(cU) LINE_MC_AUX(cA)
        LINE_MC_AUX(cV)
#undef LINE_MC_AUX
        aux["sigma"] = encode_vector(r.sigma);
        aux["W"] = encode_vector(r.W);
        aux["Qlead"] = encode_vector(r.Qlead);
        aux["Rlead"] = encode_vector(r.Rlead);
        out["aux"] = aux;
        return out;
    }
};
 
/**
 * result = npfqn_feedback_elim(P, rho, 'cs2', .., 'lambda', .., 'immediateOnly', ..).
 * An absent cs2 or lambda is MATLAB's empty option, which skips that output.
 */
struct OpNpfqnFeedbackElim {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> P = a.matrix<T>("P");
        const std::vector<T> rho = a.vector<T>("rho");
        const std::vector<T> cs2 = a.vector_or_empty<T>("cs2");
        const std::vector<T> lambda = a.vector_or_empty<T>("lambda");
        const bool hasImm = a.has("immediateOnly");
        const bool imm = hasImm ? a.boolean("immediateOnly", false) : false;
        a.done();
        const npfqn::FeedbackElimResult<T> r = hasImm ? npfqn::npfqn_feedback_elim(P, rho, cs2, lambda, imm)
                                                      : npfqn::npfqn_feedback_elim(P, rho, cs2, lambda);
        Json res;
        res["feedbackProb"] = encode_vector(r.feedbackProb);
        res["visitInflation"] = encode_vector(r.visitInflation);
        if (!r.modifiedScv.empty()) res["modifiedScv"] = encode_vector(r.modifiedScv);
        if (!r.modifiedRates.empty()) res["modifiedRates"] = encode_vector(r.modifiedRates);
        res["modifiedRouting"] = encode_matrix(r.modifiedRouting);
        res["reductionExact"] = r.reductionExact;
        Json out;
        out["result"] = res;
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// snc: every envelope function returns [sigma, rho]
// ---------------------------------------------------------------------------
 
Json env_json(const snc::Env& e) {
    Json out;
    out["sigma"] = e.sigma;
    out["rho"] = e.rho;
    return out;
}
 
/** The (sigma, rho) pair of an envelope passed as two scalars. */
snc::Env read_env(Args& a, const char* s, const char* r) {
    snc::Env e;
    e.sigma = a.number<double>(s);
    e.rho = a.number<double>(r);
    return e;
}
 
struct OpSncEnvPoisson {
    static Json run(Args& a) {
        const double lambda = a.number<double>("lambda"), theta = a.number<double>("theta");
        a.done();
        return env_json(snc::snc_env_poisson(lambda, theta));
    }
};
 
struct OpSncEnvCpoisson {
    static Json run(Args& a) {
        const double lambda = a.number<double>("lambda"), mu = a.number<double>("mu");
        const double theta = a.number<double>("theta");
        a.done();
        return env_json(snc::snc_env_cpoisson(lambda, mu, theta));
    }
};
 
/** [sigma,rho] = snc_env_tokenbucket(b, r, theta); theta is optional and only validated. */
struct OpSncEnvTokenbucket {
    static Json run(Args& a) {
        const double b = a.number<double>("b"), r = a.number<double>("r");
        const bool hasTheta = a.has("theta");
        const double theta = hasTheta ? a.number<double>("theta") : 0.0;
        a.done();
        return env_json(hasTheta ? snc::snc_env_tokenbucket(b, r, theta) : snc::snc_env_tokenbucket(b, r));
    }
};
 
struct OpSncEnvMap {
    static Json run(Args& a) {
        const Matrix<double> D0 = a.matrix<double>("D0"), D1 = a.matrix<double>("D1");
        const double theta = a.number<double>("theta");
        a.done();
        return env_json(snc::snc_env_map(D0, D1, theta));
    }
};
 
/** [sigma,rho] = snc_srv_rate(C, theta); theta is optional and only validated. */
struct OpSncSrvRate {
    static Json run(Args& a) {
        const double C = a.number<double>("C");
        const bool hasTheta = a.has("theta");
        const double theta = hasTheta ? a.number<double>("theta") : 0.0;
        a.done();
        return env_json(hasTheta ? snc::snc_srv_rate(C, theta) : snc::snc_srv_rate(C));
    }
};
 
struct OpSncSrvExp {
    static Json run(Args& a) {
        const double mu = a.number<double>("mu"), theta = a.number<double>("theta");
        a.done();
        return env_json(snc::snc_srv_exp(mu, theta));
    }
};
 
struct OpSncLeftover {
    static Json run(Args& a) {
        const snc::Env srv = read_env(a, "sigmaS", "rhoS");
        const snc::Env cross = read_env(a, "sigmaX", "rhoX");
        a.done();
        return env_json(snc::snc_leftover(srv, cross));
    }
};
 
/** [sigma,rho] = snc_conv(sigma1,rho1,sigma2,rho2,theta,delta). */
struct OpSncConv {
    static Json run(Args& a) {
        const snc::Env s1 = read_env(a, "sigma1", "rho1");
        const snc::Env s2 = read_env(a, "sigma2", "rho2");
        const double theta = a.number<double>("theta");
        const bool hasDelta = a.has("delta");
        const double delta = hasDelta ? a.number<double>("delta") : 0.0;
        a.done();
        return env_json(hasDelta ? snc::snc_conv(s1, s2, theta, delta) : snc::snc_conv(s1, s2, theta));
    }
};
 
struct OpSncOutput {
    static Json run(Args& a) {
        const snc::Env arv = read_env(a, "sigmaA", "rhoA");
        const snc::Env srv = read_env(a, "sigmaS", "rhoS");
        const double theta = a.number<double>("theta");
        a.done();
        return env_json(snc::snc_output(arv, srv, theta));
    }
};
 
 
/**
 * Q = ctmc_rand(n, seed). The uniforms come from mc::LcgUniform seeded by the
 * caller: its variates are dyadic rationals, so the exact arm returns the same
 * generator the double arm does, entry for entry. MATLAB draws from the session
 * RNG instead, so the matrices agree with it in law, not draw for draw.
 */
struct OpCtmcRand {
    template <class T>
    static Json run(Args& a) {
        const std::size_t n = a.count("n");
        mc::LcgUniform gen(required_seed(a, "ctmc_rand"));
        a.done();
        Json out;
        out["Q"] = encode_matrix(mc::ctmc_rand<T>(n, gen));
        return out;
    }
};
 
/**
 * [soujt, sts] = ctmc_simulate(Q, pi0, n, seed). pi0 may be empty, which draws
 * the initial state the way the reference does; sts is 0-based.
 */
struct OpCtmcSimulate {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const std::vector<T> pi0 = a.vector_or_empty<T>("pi0");
        const std::size_t n = a.count("n");
        pfqn::McRng rng(required_seed(a, "ctmc_simulate"));
        a.done();
        const mc::CtmcPath<T> p = mc::ctmc_simulate(Q, pi0, n, rng);
        Json out;
        out["soujt"] = encode_vector(p.sojourn);
        out["sts"] = encode_indices(p.states);
        return out;
    }
};
 
}  // namespace
 
void register_api_mc(DispatchTable& t) {
    // mc
    t["ctmc_stochcomp"] = make<OpCtmcStochcomp>();
    t["ctmc_uniformization"] = make<NeedsTranscendental<OpCtmcUniformization> >();
    t["ctmc_timeaverage"] = make<NeedsTranscendental<OpCtmcTimeaverage> >();
    t["ctmc_accumulated"] = make<NeedsTranscendental<OpCtmcAccumulated> >();
    t["ctmc_power"] = make<NeedsTranscendental<OpCtmcClassical<mc::detail::SweepKind::Power> > >();
    t["ctmc_jacobi"] = make<NeedsTranscendental<OpCtmcClassical<mc::detail::SweepKind::Jacobi> > >();
    t["ctmc_gaussseidel"] = make<NeedsTranscendental<OpCtmcClassical<mc::detail::SweepKind::GaussSeidel> > >();
    t["ctmc_sor"] = make<NeedsTranscendental<OpCtmcClassical<mc::detail::SweepKind::Sor> > >();
    t["ctmc_gth"] = make<OpGth<false> >();
    t["dtmc_gth"] = make<OpGth<true> >();
    t["ctmc_cesaro"] = make<OpCesaro<false> >();
    t["dtmc_cesaro"] = make<OpCesaro<true> >();
    t["ctmc_ssd"] = make<DoubleOnly<OpSsd> >();
    t["dtmc_ssd"] = make<DoubleOnly<OpSsd> >();
    t["ctmc_memory_gate"] = make<NeedsTranscendental<OpCtmcMemoryGate> >();
    t["ctmc_isfeasible"] = make<OpCtmcIsfeasible>();
    t["ctmc_hitting_time"] = make<NeedsTranscendental<OpCtmcHittingTime> >();
    // The registry lists Exact for dtmc_hitting_time, but the port static_asserts
    // an arithmetic with an infinity (an unreachable state's time is Inf), so the
    // exact arm is refused with that reason instead of failing to compile.
    t["dtmc_hitting_time"] = make<NeedsTranscendental<OpDtmcHittingTime> >();
    t["dtmc_transient"] = make<OpDtmcTransient>();
    t["ctmc_passage_ph"] = make<OpCtmcPassagePh>();
    t["ctmc_passage_lst"] = make<NeedsTranscendental<OpCtmcPassageLst> >();
    t["ctmc_passage_moments"] = make<NeedsTranscendental<OpCtmcPassageMoments> >();
    t["ctmc_passage_time"] = make<NeedsTranscendental<OpCtmcPassageTime> >();
    t["ctmc_passage_taboo"] = make<NeedsTranscendental<OpCtmcPassageTaboo> >();
    t["ctmc_passage_taboo_moments"] = make<NeedsTranscendental<OpCtmcPassageTabooMoments> >();
    // dpfqn
    t["dpfqn_nc"] = make<OpDpfqnNc>();
    t["dpfqn_ncld"] = make<OpDpfqnNcld>();
    // npfqn
    t["npfqn_traffic_idc"] = make<OpNpfqnTrafficIdc>();
    t["npfqn_traffic_split_cs"] = make<OpNpfqnTrafficSplitCs>();
    t["npfqn_traffic_merge_cs"] = make<OpNpfqnTrafficMergeCs>();
    t["npfqn_traffic_merge"] = make<NeedsTranscendental<OpNpfqnTrafficMerge> >();
    t["npfqn_rqna_weight"] = make<NeedsTranscendental<OpNpfqnRqnaWeight> >();
    t["npfqn_traffic_rqt"] = make<NeedsTranscendental<OpNpfqnTrafficRqt> >();
    t["npfqn_bnd_bpt"] = make<OpNpfqnBndBpt>();
    t["npfqn_bnd_bgt"] = make<OpNpfqnBndBgt>();
    t["npfqn_nonexp_approx"] = make<NeedsTranscendental<OpNpfqnNonexpApprox> >();
    t["npfqn_sqd"] = make<NeedsTranscendental<OpNpfqnSqd> >();
    t["npfqn_dps_morrison"] = make<NeedsTranscendental<OpNpfqnDpsMorrison> >();
    t["npfqn_feedback_elim"] = make<OpNpfqnFeedbackElim>();
    // snc (double-only kernels)
    t["snc_env_poisson"] = make<DoubleOnly<OpSncEnvPoisson> >();
    t["snc_env_cpoisson"] = make<DoubleOnly<OpSncEnvCpoisson> >();
    t["snc_env_tokenbucket"] = make<DoubleOnly<OpSncEnvTokenbucket> >();
    t["snc_env_map"] = make<DoubleOnly<OpSncEnvMap> >();
    t["snc_srv_rate"] = make<DoubleOnly<OpSncSrvRate> >();
    t["snc_srv_exp"] = make<DoubleOnly<OpSncSrvExp> >();
    t["snc_leftover"] = make<DoubleOnly<OpSncLeftover> >();
    t["snc_conv"] = make<DoubleOnly<OpSncConv> >();
    t["snc_output"] = make<DoubleOnly<OpSncOutput> >();
    t["ctmc_rand"] = make<OpCtmcRand>();
    t["ctmc_simulate"] = make<NeedsTranscendental<OpCtmcSimulate> >();
}
 
void register_library_only_mc(LibraryOnlyTable& t) {
    const std::string sn =
        "it takes a NetworkStruct; build the model and run line-cli -s ctmc -a <analysis> "
        "on it instead.";
    t["ctmc_ssg"] = sn;
    t["ctmc_ssg_reachability"] = sn;
    t["ctmc_state_space_logsize"] = sn;
    t["npfqn_traffic_split_rr"] =
        "it takes a NetworkStruct; build the model and run line-cli -s <solver> -a <analysis> "
        "on it instead.";
    t["npfqn_traffic_idc_at"] =
        "it evaluates two caller-supplied IDC functions (a0IdcFun and sIdcFun), which are "
        "callables with no JSON form.";
    const std::string env =
        "its arrival and service envelopes are functions of theta (e.g. snc_env_poisson as a "
        "handle), callables with no JSON form; the envelope values at one theta are exposed "
        "through the snc_env_* and snc_srv_* entries.";
    t["snc_bound_backlog"] = env;
    t["snc_bound_delay"] = env;
    t["snc_perc_backlog"] = env;
    t["snc_perc_delay"] = env;
    t["snc_mean_backlog"] = env;
    t["snc_mean_delay"] = env;
    t["snc_thetaopt"] =
        "its objective is a function of theta passed as a callable, which has no JSON form.";
}
 
}  // namespace reg
}  // namespace line