api_dispatch_pfqn_a.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 'pfqn_a' group (see api_dispatch_detail.h).
 *
 * Conventions shared by every entry below:
 *  - argument keys are the MATLAB parameter names of matlab/src/api/pfqn/<fn>.m,
 *    output keys its output names (or the C++ field name where they agree);
 *  - an ABSENT optional argument selects the port's own shorter overload, or
 *    leaves a trailing C++ default argument unwritten. No default value is
 *    restated here. When the port takes optional arguments POSITIONALLY and the
 *    caller sets a later one while leaving an earlier one out, the earlier one
 *    has no way to stay at its default and the call is refused, naming it
 *    (`depth` below). An empty matrix or vector is the port's own "not given"
 *    for the arguments marked fillable, so those never block a later one;
 *  - node, station, class and path indices are 0-based and passed through
 *    unchanged, as read_partition/encode_indices do in api_dispatch.cpp;
 *  - a scheduling vector is a list of names, "PS", "FCFS" or "INF";
 *  - a convergence flag is reported whenever the port returns one.
 */
 
#include <cmath>
#include <cstdint>
#include <initializer_list>
#include <string>
#include <utility>
#include <vector>
 
#include "api_dispatch_detail.h"
 
#include "line/api/pfqn/pfqn_ab_amva.h"
#include "line/api/pfqn/pfqn_aghq.h"
#include "line/api/pfqn/pfqn_aql.h"
#include "line/api/pfqn/pfqn_bk.h"
#include "line/api/pfqn/pfqn_bkt.h"
#include "line/api/pfqn/pfqn_ble.h"
#include "line/api/pfqn/pfqn_bs.h"
#include "line/api/pfqn/pfqn_busyp.h"
#include "line/api/pfqn/pfqn_busyp_clw.h"
#include "line/api/pfqn/pfqn_busyp_multiclass.h"
#include "line/api/pfqn/pfqn_cbh.h"
#include "line/api/pfqn/pfqn_cftp.h"
#include "line/api/pfqn/pfqn_chow.h"
#include "line/api/pfqn/pfqn_clust.h"
#include "line/api/pfqn/pfqn_clw.h"
#include "line/api/pfqn/pfqn_cntol.h"
#include "line/api/pfqn/pfqn_comomrm_ld.h"
#include "line/api/pfqn/pfqn_conwayms.h"
#include "line/api/pfqn/pfqn_cub.h"
#include "line/api/pfqn/pfqn_cub_evals.h"
#include "line/api/pfqn/pfqn_cyclet_ofree.h"
#include "line/api/pfqn/pfqn_dac.h"
#include "line/api/pfqn/pfqn_divdiff_ld.h"
#include "line/api/pfqn/pfqn_dmlin.h"
#include "line/api/pfqn/pfqn_dnc.h"
#include "line/api/pfqn/pfqn_egflinearizer.h"
#include "line/api/pfqn/pfqn_expand.h"
#include "line/api/pfqn/pfqn_explicit.h"
#include "line/api/pfqn/pfqn_explicit_ld.h"
#include "line/api/pfqn/pfqn_fnc.h"
#include "line/api/pfqn/pfqn_gerasimov.h"
#include "line/api/pfqn/pfqn_gflinearizer.h"
#include "line/api/pfqn/pfqn_gldsingle.h"
#include "line/api/pfqn/pfqn_grnmol.h"
#include "line/api/pfqn/pfqn_harel_bounds.h"
#include "line/api/pfqn/pfqn_hst.h"
#include "line/api/pfqn/pfqn_is.h"
#include "line/api/pfqn/pfqn_joint.h"
#include "line/api/pfqn/pfqn_jointmarg.h"
#include "line/api/pfqn/pfqn_kt.h"
#include "line/api/pfqn/pfqn_lap.h"
#include "line/api/pfqn/pfqn_lcfsqn_ca.h"
#include "line/api/pfqn/pfqn_lcfsqn_mva.h"
#include "line/api/pfqn/pfqn_lcfsqn_nc.h"
#include "line/api/pfqn/pfqn_lcp.h"
#include "line/api/pfqn/pfqn_ld_is.h"
#include "line/api/pfqn/pfqn_ldbcmp.h"
#include "line/api/pfqn/pfqn_ldmx_ec.h"
#include "line/api/pfqn/pfqn_le.h"
#include "line/api/pfqn/pfqn_lekt.h"
#include "line/api/pfqn/pfqn_linearizer.h"
#include "line/api/pfqn/pfqn_linearizerms.h"
#include "line/api/pfqn/pfqn_linearizermx.h"
#include "line/api/pfqn/pfqn_lld.h"
#include "line/api/pfqn/pfqn_lldfun.h"
#include "line/api/pfqn/pfqn_lldsingle.h"
#include "line/api/pfqn/pfqn_looping.h"
#include "line/api/pfqn/pfqn_ls.h"
 
namespace line {
namespace reg {
 
namespace {
 
// ---------------------------------------------------------------------------
// Shared helpers
// ---------------------------------------------------------------------------
 
/** A positional optional argument: its key, and whether an empty value is the
 *  port's own "not given" (so leaving it out never blocks a later one). */
struct Pos {
    const char* key;
    bool fillable;
};
 
/**
 * How many of the trailing positional optionals to hand the port: one past the
 * last one present. An absent one before it that the port cannot take as empty
 * is refused, since writing its default here would restate the port's value.
 */
std::size_t depth(const Args& a, const char* fn, std::initializer_list<Pos> keys) {
    const std::vector<Pos> k(keys);
    std::size_t d = 0;
    for (std::size_t i = 0; i < k.size(); ++i)
        if (a.has(k[i].key)) d = i + 1;
    for (std::size_t i = 0; i + 1 < d; ++i) {
        if (k[i].fillable || a.has(k[i].key)) continue;
        std::string later;
        for (std::size_t j = d; j-- > i + 1;)
            if (a.has(k[j].key)) later = k[j].key;
        throw InputError(std::string(fn) + ": '" + later + "' is given but '" + k[i].key +
                         "' is not; the C++ port takes '" + k[i].key + "' positionally before '" +
                         later + "', so it cannot be left to its default while '" + later +
                         "' is set. Give '" + k[i].key + "' as well.");
    }
    return d;
}
 
/**
 * 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;
}
 
/** A scheduling name, case-insensitive. */
int sched_code(const Json& e, const std::string& where) {
    if (!e.is_string())
        throw InputError(where + ": a scheduling entry must be one of \"PS\", \"FCFS\" or \"INF\", got " +
                         std::string(e.type_name()));
    std::string s = e.get<std::string>();
    for (char& c : s) c = static_cast<char>(std::toupper(static_cast<unsigned char>(c)));
    if (s == "PS") return 0;
    if (s == "FCFS") return 1;
    if (s == "INF") return 2;
    throw InputError(where + ": unknown scheduling strategy '" + e.get<std::string>() +
                     "'; expected \"PS\", \"FCFS\" or \"INF\"");
}
 
/** A per-station scheduling vector; absent reads as empty, the ports' "not given". */
template <class E>
std::vector<E> sched_vector(Args& a, const char* key, E ps, E fcfs, E inf) {
    std::vector<E> out;
    const Json* j = a.opt(key);
    if (!j) return out;
    const std::string where = std::string("'") + key + "'";
    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) {
        const int c = sched_code(e, where);
        out.push_back(c == 0 ? ps : (c == 1 ? fcfs : inf));
    }
    return out;
}
 
std::vector<pfqn::AmvaSched> amva_sched(Args& a, const char* key) {
    return sched_vector(a, key, pfqn::AmvaSched::PS, pfqn::AmvaSched::FCFS, pfqn::AmvaSched::INF);
}
 
std::vector<pfqn::SchedStrategy> lin_sched(Args& a, const char* key) {
    return sched_vector(a, key, pfqn::SchedStrategy::PS, pfqn::SchedStrategy::FCFS,
                        pfqn::SchedStrategy::INF);
}
 
/** A double option; with `cn` the string "cn" selects the Chandy-Neuse test (NaN). */
double real_option(Args& a, const char* key, bool cn) {
    const Json& j = a.get(key);
    if (cn && j.is_string()) {
        std::string s = j.get<std::string>();
        for (char& c : s) c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
        if (s == "cn") return std::nan("");
    }
    return number_from_json<double>(j, std::string("'") + key + "'");
}
 
/** A scalar integer that must be given. */
int int_arg(Args& a, const char* key) { return a.required_integer(key); }
 
/** A list of 0-based indices, passed through unchanged. */
std::vector<std::size_t> index_list(const Json& j, const std::string& where) {
    const std::vector<int> v = int_vector_from_json(j, where);
    std::vector<std::size_t> out;
    out.reserve(v.size());
    for (int x : v) {
        if (x < 0) throw InputError(where + ": an index cannot be negative, got " + std::to_string(x));
        out.push_back(static_cast<std::size_t>(x));
    }
    return out;
}
 
std::vector<std::size_t> indices(Args& a, const char* key) {
    return index_list(a.get(key), std::string("'") + key + "'");
}
 
std::vector<std::size_t> indices_or_empty(Args& a, const char* key) {
    const Json* j = a.opt(key);
    return j ? index_list(*j, std::string("'") + key + "'") : std::vector<std::size_t>();
}
 
/** A list of index lists; a flat list is a single one when `flat_is_one`. */
std::vector<std::vector<std::size_t> > index_lists(const Json& j, const std::string& where,
                                                    bool flat_is_one) {
    std::vector<std::vector<std::size_t> > out;
    if (!j.is_array()) throw InputError(where + ": expected an array of index arrays");
    bool nested = false;
    for (const Json& e : j)
        if (e.is_array()) nested = true;
    if (!nested && !j.empty()) {
        if (!flat_is_one) throw InputError(where + ": expected an array of index arrays");
        out.push_back(index_list(j, where));
        return out;
    }
    for (std::size_t b = 0; b < j.size(); ++b) {
        if (!j[b].is_array())
            throw InputError(where + ": entry " + std::to_string(b) + " must be an array");
        out.push_back(index_list(j[b], where + "[" + std::to_string(b) + "]"));
    }
    return out;
}
 
/** An integer matrix (2-D array as rows, 1-D array as one row). */
Matrix<int> int_matrix(Args& a, const char* key) {
    const Json& j = a.get(key);
    const std::string where = std::string("'") + key + "'";
    if (j.is_array() && !j.empty() && j[0].is_array()) {
        const std::size_t rows = j.size(), cols = j[0].size();
        Matrix<int> m(rows, cols, 0);
        for (std::size_t i = 0; i < rows; ++i) {
            if (!j[i].is_array() || j[i].size() != cols)
                throw InputError(where + ": row " + std::to_string(i) + " has the wrong length");
            const std::vector<int> r = int_vector_from_json(j[i], where);
            for (std::size_t k = 0; k < cols; ++k) m(i, k) = r[k];
        }
        return m;
    }
    const std::vector<int> r = int_vector_from_json(j, where);
    Matrix<int> m(1, r.size(), 0);
    for (std::size_t k = 0; k < r.size(); ++k) m(0, k) = r[k];
    return m;
}
 
/** A population vector in which an infinite entry marks an open class. */
std::vector<int> mixed_population(Args& a, const char* key) {
    const Json& j = a.get(key);
    const std::string where = std::string("'") + key + "'";
    std::vector<int> out;
    for (double d : vector_from_json<double>(j, where)) {
        if (std::isinf(d) && d > 0) {
            out.push_back(pfqn::kOpenClass);
            continue;
        }
        if (d != std::floor(d) || d < 0)
            throw InputError(where + ": a closed-class population must be a non-negative integer");
        out.push_back(static_cast<int>(d));
    }
    return out;
}
 
/** A boolean vector, from true/false or 0/1 entries. */
std::vector<bool> bool_vector(const Json& j, const std::string& where) {
    std::vector<bool> 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>());
            continue;
        }
        const std::vector<int> v = int_vector_from_json(e, where);
        if (v.size() != 1 || (v[0] != 0 && v[0] != 1))
            throw InputError(where + ": entries must be true, false, 0 or 1");
        out.push_back(v[0] != 0);
    }
    return out;
}
 
Json encode_index_list(const std::vector<std::size_t>& v) {
    Json a = Json::array();
    for (std::size_t x : v) a.push_back(static_cast<std::uint64_t>(x));
    return a;
}
 
/** exp of a log normalizing constant, at an arithmetic that has it. */
template <class T>
T exp_of_log(Args& a, const char* key) {
    if constexpr (num_traits<T>::has_transcendental) {
        using std::exp;
        const T l = a.number<T>(key);
        return T(exp(l));
    } else {
        throw UnsupportedError(ErrorKind{"Arithmetic"},
                               std::string("'") + key + "' is a logarithm, and exp() is not in the exact "
                               "rational field; omit it to have the port compute G(N) exactly, or rerun "
                               "with --arith double or --arith real:<digits>.");
    }
}
 
// ---------------------------------------------------------------------------
// Result shapes
// ---------------------------------------------------------------------------
 
template <class T>
Json amva_results(const pfqn::AmvaResult<T>& r, const char* itkey) {
    Json out;
    out["XN"] = encode_vector(r.XN);
    out["QN"] = encode_matrix(r.QN);
    out["UN"] = encode_matrix(r.UN);
    out["RN"] = encode_matrix(r.RN);
    out[itkey] = encode_count(r.iterations);
    out["converged"] = r.converged;
    return out;
}
 
/** The six Linearizer outputs, under the names the MATLAB function gives them. */
template <class T>
Json lin_results(const pfqn::LinearizerResult<T>& r, const char* q, const char* u, const char* w,
                 const char* c, const char* x) {
    Json out;
    out[q] = encode_matrix(r.Q);
    out[u] = encode_matrix(r.U);
    out[w] = encode_matrix(r.W);
    out[c] = encode_vector(r.C);
    out[x] = encode_vector(r.X);
    out["totiter"] = r.totiter;
    return out;
}
 
template <class T>
Json lin_results(const pfqn::LinearizerResult<T>& r) {
    return lin_results(r, "Q", "U", "W", "C", "X");
}
 
template <class G, class L>
Json g_lg(const G& g, const L& lg, const char* gk = "G", const char* lk = "lG") {
    Json out;
    out[gk] = encode_scalar(g);
    out[lk] = encode_scalar(lg);
    return out;
}
 
// ---------------------------------------------------------------------------
// Entries
// ---------------------------------------------------------------------------
 
struct OpCycletOfree {
    static Json run(Args& a) {
        const std::vector<double> v = a.vector<double>("v");
        const std::vector<double> mu = a.vector<double>("mu");
        const std::size_t N = a.count("N");
        const std::vector<std::vector<std::size_t> > paths =
            index_lists(a.get("path"), "'path'", true);
        const std::vector<double> tset = a.vector<double>("tset");
        pfqn::CycletResult r;
        const Json* oj = a.opt("options");
        a.done();
        if (!oj) {
            r = pfqn::pfqn_cyclet_ofree(v, mu, N, paths, tset);
        } else {
            Args o(*oj, "pfqn_cyclet_ofree options");
            const std::size_t d = depth(o, "pfqn_cyclet_ofree options",
                                        {{"method", false}, {"nmom", false}, {"pathprob", true},
                                         {"lti_method", false}, {"tol", false}});
            const std::string method = o.has("method") ? o.text("method", "") : std::string();
            const std::size_t nmom = o.has("nmom") ? o.count("nmom") : 0;
            const std::vector<double> pathprob = o.vector_or_empty<double>("pathprob");
            const std::string lti = o.has("lti_method") ? o.text("lti_method", "") : std::string();
            const double tol = o.has("tol") ? o.number<double>("tol") : 0.0;
            o.done();
            switch (d) {
                case 0: r = pfqn::pfqn_cyclet_ofree(v, mu, N, paths, tset); break;
                case 1: r = pfqn::pfqn_cyclet_ofree(v, mu, N, paths, tset, method); break;
                case 2: r = pfqn::pfqn_cyclet_ofree(v, mu, N, paths, tset, method, nmom); break;
                case 3:
                    r = pfqn::pfqn_cyclet_ofree(v, mu, N, paths, tset, method, nmom, pathprob);
                    break;
                case 4:
                    r = pfqn::pfqn_cyclet_ofree(v, mu, N, paths, tset, method, nmom, pathprob, lti);
                    break;
                default:
                    r = pfqn::pfqn_cyclet_ofree(v, mu, N, paths, tset, method, nmom, pathprob, lti,
                                                tol);
            }
        }
        Json out;
        out["f"] = encode_vector(r.f);
        out["F"] = encode_vector(r.F);
        out["mom"] = encode_vector(r.mom);
        Json info = Json::array();
        for (const pfqn::CycletPathInfo& p : r.info) {
            Json e;
            e["method"] = p.method;
            e["lG"] = p.lG;
            e["path"] = encode_index_list(p.path);
            info.push_back(e);
        }
        out["out"] = info;
        return out;
    }
};
 
struct OpExpand {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> QN = a.matrix<T>("QN");
        const Matrix<T> UN = a.matrix<T>("UN");
        const Matrix<T> CN = a.matrix<T>("CN");
        std::vector<std::size_t> mapping = indices(a, "mapping");
        const std::size_t M = a.count("M_original");
        a.done();
        // MATLAB fills rows 1..M_original from mapping(1..M_original).
        if (mapping.size() < M)
            throw InputError("pfqn_expand: mapping has " + std::to_string(mapping.size()) +
                             " entries, fewer than M_original = " + std::to_string(M));
        mapping.resize(M);
        const pfqn::ExpandResult<T> r = pfqn::pfqn_expand(QN, UN, CN, mapping);
        Json out;
        out["QN_full"] = encode_matrix(r.QN);
        out["UN_full"] = encode_matrix(r.UN);
        out["CN_full"] = encode_matrix(r.CN);
        return out;
    }
};
 
template <class T>
Json explicit_results(const pfqn::ExplicitResult<T>& r) {
    Json out;
    out["lG"] = encode_scalar(r.lG);
    out["G"] = encode_scalar(r.G);
    out["method"] = r.method;
    out["lossDigits"] = r.lossDigits;
    out["valid"] = r.valid;
    return out;
}
 
struct OpExplicit {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const std::size_t d =
            depth(a, "pfqn_explicit", {{"tol", false}, {"method", false}, {"maxloss", false}});
        const double tol = a.has("tol") ? a.number<double>("tol") : 0.0;
        const std::string method = a.text("method", "");
        const double maxloss = a.has("maxloss") ? a.number<double>("maxloss") : 0.0;
        a.done();
        switch (d) {
            case 0: return explicit_results(pfqn::pfqn_explicit(L, N));
            case 1: return explicit_results(pfqn::pfqn_explicit(L, N, tol));
            case 2: return explicit_results(pfqn::pfqn_explicit(L, N, tol, method));
            default: return explicit_results(pfqn::pfqn_explicit(L, N, tol, method, maxloss));
        }
    }
};
 
struct OpExplicitLd {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const Matrix<T> mu = a.matrix<T>("mu");
        const std::size_t d =
            depth(a, "pfqn_explicit_ld", {{"tol", false}, {"method", false}, {"maxloss", false}});
        const double tol = a.has("tol") ? a.number<double>("tol") : 0.0;
        const std::string method = a.text("method", "");
        const double maxloss = a.has("maxloss") ? a.number<double>("maxloss") : 0.0;
        a.done();
        switch (d) {
            case 0: return explicit_results(pfqn::pfqn_explicit_ld(L, N, mu));
            case 1: return explicit_results(pfqn::pfqn_explicit_ld(L, N, mu, tol));
            case 2: return explicit_results(pfqn::pfqn_explicit_ld(L, N, mu, tol, method));
            default:
                return explicit_results(pfqn::pfqn_explicit_ld(L, N, mu, tol, method, maxloss));
        }
    }
};
 
struct OpDivdiffLd {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const bool hasZ = a.has("Z"), hasMu = a.has("mu");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const Matrix<T> mu = a.matrix_or_empty<T>("mu");
        a.unsupported("options", "the C++ port of pfqn_divdiff_ld takes no solver options");
        a.done();
        pfqn::DivdiffLdResult<T> r = hasMu ? pfqn::pfqn_divdiff_ld(L, N, Z, mu)
                                           : (hasZ ? pfqn::pfqn_divdiff_ld(L, N, Z)
                                                   : pfqn::pfqn_divdiff_ld(L, N));
        Json out;
        out["lG"] = encode_scalar(r.lG);
        out["G"] = encode_scalar(r.G);
        out["lossDigits"] = r.lossDigits;
        return out;
    }
};
 
/** pfqn_bs, pfqn_lcp and pfqn_chow share MATLAB's (L,N,Z,tol,maxiter,QN0,type). */
template <class T>
struct AmvaArgs {
    Matrix<T> L;
    std::vector<T> N, Z;
    bool hasZ;
    std::vector<pfqn::AmvaSched> type;
    double tol;
    std::size_t maxiter;
    Matrix<T> QN0;
    std::size_t d;  // how many of tol, maxiter, QN0 to pass
    bool longf;
};
 
template <class T>
AmvaArgs<T> read_amva(Args& a, const char* fn, bool cn, std::initializer_list<Pos> extra) {
    AmvaArgs<T> s;
    s.L = a.matrix<T>("L");
    s.N = a.vector<T>("N");
    s.hasZ = a.has("Z");
    s.Z = a.vector_or_empty<T>("Z");
    std::vector<Pos> keys = {{"tol", false}, {"maxiter", false}, {"QN0", true}};
    keys.insert(keys.end(), extra.begin(), extra.end());
    // depth() takes an initializer_list; spell out the two shapes used here.
    if (keys.size() == 3)
        s.d = depth(a, fn, {{"tol", false}, {"maxiter", false}, {"QN0", true}});
    else
        s.d = depth(a, fn, {{"tol", false}, {"maxiter", false}, {"QN0", true}, keys[3]});
    s.longf = s.d > 0 || a.has("type");
    s.type = amva_sched(a, "type");
    s.tol = a.has("tol") ? real_option(a, "tol", cn) : 0.0;
    s.maxiter = a.has("maxiter") ? a.count("maxiter") : 0;
    s.QN0 = a.matrix_or_empty<T>("QN0");
    return s;
}
 
struct OpBs {
    template <class T>
    static Json run(Args& a) {
        AmvaArgs<T> s = read_amva<T>(a, "pfqn_bs", true, {});
        a.done();
        pfqn::AmvaResult<T> r;
        if (!s.longf)
            r = s.hasZ ? pfqn::pfqn_bs(s.L, s.N, s.Z) : pfqn::pfqn_bs(s.L, s.N);
        else if (s.d == 0)
            r = pfqn::pfqn_bs(s.L, s.N, s.Z, s.type);
        else if (s.d == 1)
            r = pfqn::pfqn_bs(s.L, s.N, s.Z, s.type, s.tol);
        else if (s.d == 2)
            r = pfqn::pfqn_bs(s.L, s.N, s.Z, s.type, s.tol, s.maxiter);
        else
            r = pfqn::pfqn_bs(s.L, s.N, s.Z, s.type, s.tol, s.maxiter, s.QN0);
        return amva_results(r, "it");
    }
};
 
struct OpLcp {
    template <class T>
    static Json run(Args& a) {
        AmvaArgs<T> s = read_amva<T>(a, "pfqn_lcp", false, {});
        a.done();
        pfqn::AmvaResult<T> r;
        if (!s.longf)
            r = s.hasZ ? pfqn::pfqn_lcp(s.L, s.N, s.Z) : pfqn::pfqn_lcp(s.L, s.N);
        else if (s.d == 0)
            r = pfqn::pfqn_lcp(s.L, s.N, s.Z, s.type);
        else if (s.d == 1)
            r = pfqn::pfqn_lcp(s.L, s.N, s.Z, s.type, s.tol);
        else if (s.d == 2)
            r = pfqn::pfqn_lcp(s.L, s.N, s.Z, s.type, s.tol, s.maxiter);
        else
            r = pfqn::pfqn_lcp(s.L, s.N, s.Z, s.type, s.tol, s.maxiter, s.QN0);
        return amva_results(r, "it");
    }
};
 
pfqn::ChowVariant chow_variant(Args& a) {
    std::string v = a.text("variant", "");
    for (char& c : v) c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
    if (v == "forward") return pfqn::ChowVariant::Forward;
    if (v == "backward") return pfqn::ChowVariant::Backward;
    throw InputError("pfqn_chow: unknown variant '" + v + "'; expected \"forward\" or \"backward\"");
}
 
struct OpChow {
    template <class T>
    static Json run(Args& a) {
        AmvaArgs<T> s = read_amva<T>(a, "pfqn_chow", false, {{"variant", false}});
        const bool hasVariant = a.has("variant");
        const pfqn::ChowVariant variant = hasVariant ? chow_variant(a) : pfqn::ChowVariant::Forward;
        a.done();
        pfqn::AmvaResult<T> r;
        if (!s.longf && !hasVariant)
            r = s.hasZ ? pfqn::pfqn_chow(s.L, s.N, s.Z) : pfqn::pfqn_chow(s.L, s.N);
        else if (s.d == 0)
            r = pfqn::pfqn_chow(s.L, s.N, s.Z, s.type);
        else if (s.d == 1)
            r = pfqn::pfqn_chow(s.L, s.N, s.Z, s.type, s.tol);
        else if (s.d == 2)
            r = pfqn::pfqn_chow(s.L, s.N, s.Z, s.type, s.tol, s.maxiter);
        else if (s.d == 3)
            r = pfqn::pfqn_chow(s.L, s.N, s.Z, s.type, s.tol, s.maxiter, s.QN0);
        else
            r = pfqn::pfqn_chow(s.L, s.N, s.Z, s.type, s.tol, s.maxiter, s.QN0, variant);
        return amva_results(r, "it");
    }
};
 
struct OpCntol {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> N = a.vector<T>("N");
        a.done();
        Json out;
        out["tol"] = pfqn::pfqn_cntol(N);
        return out;
    }
};
 
struct OpAql {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<T> N = a.vector<T>("N");
        const bool hasZ = a.has("Z");
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        const std::size_t d = depth(a, "pfqn_aql", {{"TOL", false}, {"MAXITER", false}});
        const double tol = a.has("TOL") ? a.number<double>("TOL") : 0.0;
        const std::size_t maxiter = a.has("MAXITER") ? a.count("MAXITER") : 0;
        a.unsupported("QN0", "the C++ port of pfqn_aql takes no warm start; omit QN0");
        a.done();
        pfqn::AmvaResult<T> r;
        if (d == 0)
            r = hasZ ? pfqn::pfqn_aql(L, N, Z) : pfqn::pfqn_aql(L, N);
        else if (d == 1)
            r = pfqn::pfqn_aql(L, N, Z, tol);
        else
            r = pfqn::pfqn_aql(L, N, Z, tol, maxiter);
        return amva_results(r, "numIters");
    }
};
 
/** The long-form arguments of the Linearizer family. */
template <class T>
struct LinArgs {
    std::vector<pfqn::SchedStrategy> type;
    double tol = 0.0;
    int maxiter = 0;
    Matrix<T> QN0;
};
 
template <class T>
LinArgs<T> read_lin(Args& a, bool cn) {
    LinArgs<T> s;
    s.type = lin_sched(a, "type");
    if (a.has("tol")) s.tol = real_option(a, "tol", cn);
    if (a.has("maxiter")) s.maxiter = static_cast<int>(a.count("maxiter"));
    s.QN0 = a.matrix_or_empty<T>("QN0");
    return s;
}
 
struct OpLinearizer {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const bool hasZ = a.has("Z");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const bool lf = long_form(a, "pfqn_linearizer", {"type", "tol", "maxiter", "QN0"},
                                  {"tol", "maxiter"});
        const LinArgs<T> s = read_lin<T>(a, true);
        a.done();
        if (lf) return lin_results(pfqn::pfqn_linearizer(L, N, Z, s.type, s.tol, s.maxiter, s.QN0));
        return lin_results(hasZ ? pfqn::pfqn_linearizer(L, N, Z) : pfqn::pfqn_linearizer(L, N));
    }
};
 
struct OpLinearizerms {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const std::vector<int> nservers = a.ints("nservers");
        const bool lf = long_form(a, "pfqn_linearizerms", {"type", "tol", "maxiter", "QN0"},
                                  {"tol", "maxiter"});
        const LinArgs<T> s = read_lin<T>(a, false);
        a.done();
        const pfqn::LinearizerResult<T> r =
            lf ? pfqn::pfqn_linearizerms(L, N, Z, nservers, s.type, s.tol, s.maxiter, s.QN0)
               : pfqn::pfqn_linearizerms(L, N, Z, nservers);
        return lin_results(r, "Q", "U", "R", "C", "X");
    }
};
 
struct OpConwayms {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const std::vector<int> nservers = a.ints("nservers");
        const bool lf = long_form(a, "pfqn_conwayms", {"type", "tol", "maxiter", "QN0"},
                                  {"tol", "maxiter"});
        const LinArgs<T> s = read_lin<T>(a, false);
        a.done();
        const pfqn::LinearizerResult<T> r =
            lf ? pfqn::pfqn_conwayms(L, N, Z, nservers, s.type, s.tol, s.maxiter, s.QN0)
               : pfqn::pfqn_conwayms(L, N, Z, nservers);
        return lin_results(r, "Q", "U", "R", "C", "X");
    }
};
 
pfqn::LinearizerMxMethod mx_method(Args& a) {
    std::string m = a.text("method", "");
    for (char& c : m) c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
    if (m == "lin") return pfqn::LinearizerMxMethod::Lin;
    if (m == "gflin") return pfqn::LinearizerMxMethod::Gflin;
    if (m == "egflin") return pfqn::LinearizerMxMethod::Egflin;
    throw InputError("pfqn_linearizermx: unknown method '" + m +
                     "'; expected \"lin\", \"gflin\" or \"egflin\"");
}
 
struct OpLinearizermx {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = mixed_population(a, "N");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const std::vector<int> nservers = a.ints("nservers");
        const bool lf = long_form(a, "pfqn_linearizermx", {"type", "tol", "maxiter", "QN0"},
                                  {"tol", "maxiter", "method"});
        const LinArgs<T> s = read_lin<T>(a, false);
        const bool hasMethod = a.has("method");
        const pfqn::LinearizerMxMethod method =
            hasMethod ? mx_method(a) : pfqn::LinearizerMxMethod::Egflin;
        a.done();
        pfqn::LinearizerResult<T> r;
        if (lf)
            r = pfqn::pfqn_linearizermx(lambda, L, N, Z, nservers, s.type, s.tol, s.maxiter, method,
                                        s.QN0);
        else if (hasMethod)
            r = pfqn::pfqn_linearizermx(lambda, L, N, Z, nservers, method);
        else
            r = pfqn::pfqn_linearizermx(lambda, L, N, Z, nservers);
        return lin_results(r, "QN", "UN", "WN", "CN", "XN");
    }
};
 
struct OpGflinearizer {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const bool lf = long_form(a, "pfqn_gflinearizer", {"type", "tol", "maxiter", "QN0"},
                                  {"tol", "maxiter"});
        const LinArgs<T> s = read_lin<T>(a, true);
        const T alpha = a.number<T>("alpha");
        a.done();
        if (lf)
            return lin_results(
                pfqn::pfqn_gflinearizer(L, N, Z, s.type, s.tol, s.maxiter, alpha, s.QN0));
        return lin_results(pfqn::pfqn_gflinearizer(L, N, Z, alpha));
    }
};
 
struct OpEgflinearizer {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const bool lf = long_form(a, "pfqn_egflinearizer",
                                  {"type", "tol", "maxiter", "QN0", "npasses"}, {"tol", "maxiter"});
        const LinArgs<T> s = read_lin<T>(a, true);
        const std::vector<T> alpha = a.vector<T>("alpha");
        const bool hasNp = a.has("npasses");
        const int npasses = hasNp ? int_arg(a, "npasses") : 0;
        a.done();
        if (!lf) return lin_results(pfqn::pfqn_egflinearizer(L, N, Z, alpha));
        if (hasNp)
            return lin_results(pfqn::pfqn_egflinearizer(L, N, Z, s.type, s.tol, s.maxiter, alpha,
                                                        s.QN0, npasses));
        return lin_results(
            pfqn::pfqn_egflinearizer(L, N, Z, s.type, s.tol, s.maxiter, alpha, s.QN0));
    }
};
 
struct OpDmlin {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const bool hasZ = a.has("Z");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const bool lf = long_form(a, "pfqn_dmlin", {"type", "tol", "maxiter", "QN0", "npasses"},
                                  {"tol", "maxiter"});
        const LinArgs<T> s = read_lin<T>(a, false);
        const bool hasNp = a.has("npasses");
        const int npasses = hasNp ? int_arg(a, "npasses") : 0;
        a.done();
        if (!lf) return lin_results(hasZ ? pfqn::pfqn_dmlin(L, N, Z) : pfqn::pfqn_dmlin(L, N));
        if (hasNp)
            return lin_results(
                pfqn::pfqn_dmlin(L, N, Z, s.type, s.tol, s.maxiter, s.QN0, npasses));
        return lin_results(pfqn::pfqn_dmlin(L, N, Z, s.type, s.tol, s.maxiter, s.QN0));
    }
};
 
struct OpLldfun {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> n = a.vector<T>("n");
        const Matrix<T> lldscaling = a.matrix_or_empty<T>("lldscaling");
        const bool hasC = a.has("nservers");
        const std::vector<double> nservers = a.vector_or_empty<double>("nservers");
        a.done();
        Json out;
        out["r"] = encode_vector(hasC ? pfqn::pfqn_lldfun(n, lldscaling, nservers)
                                      : pfqn::pfqn_lldfun(n, lldscaling));
        return out;
    }
};
 
struct OpComomrmLd {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const Matrix<T> mu = a.matrix<T>("mu");
        a.unsupported("options", "the C++ port of pfqn_comomrm_ld takes no solver options "
                                 "(MATLAB reads options.tol there)");
        a.done();
        const pfqn::ComomRmResult<T> r = pfqn::pfqn_comomrm_ld(L, N, Z, mu);
        Json out = g_lg(r.G, r.lG);
        out["prob"] = encode_vector(r.prob);
        return out;
    }
};
 
struct OpGldsingle {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const int N = int_arg(a, "N");
        const Matrix<T> mu = a.matrix<T>("mu");
        a.unsupported("options", "the C++ port of pfqn_gldsingle takes no solver options");
        a.done();
        const pfqn::NcResult<T> r = pfqn::pfqn_gldsingle(L, N, mu);
        return g_lg(r.G, r.lG);
    }
};
 
struct OpLldsingle {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const int N = int_arg(a, "N");
        const Matrix<T> mu = a.matrix<T>("mu");
        a.unsupported("options", "the C++ port of pfqn_lldsingle takes no solver options");
        a.done();
        const pfqn::NcResult<T> r = pfqn::pfqn_lldsingle(L, N, mu);
        return g_lg(r.G, r.lG);
    }
};
 
struct OpLld {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const bool hasMu = a.has("mu");
        const Matrix<T> mu = a.matrix_or_empty<T>("mu");
        a.unsupported("options", "the C++ port of pfqn_lld takes no solver options");
        a.done();
        const pfqn::LldResult<T> r = hasMu ? pfqn::pfqn_lld(L, N, mu) : pfqn::pfqn_lld(L, N);
        Json out = g_lg(r.G, r.lG);
        out["s"] = encode_ints(r.s);
        return out;
    }
};
 
struct OpJoint {
    template <class T>
    static Json run(Args& a) {
        const Matrix<int> n = int_matrix(a, "n");
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const bool hasG = a.has("lGN");
        const T G = hasG ? exp_of_log<T>(a, "lGN") : num_traits<T>::from_int(0);
        a.done();
        T p;
        // MATLAB: a single column is the per-station total form, R columns the per-class form.
        if (n.cols() == 1) {
            std::vector<int> m(n.rows());
            for (std::size_t i = 0; i < n.rows(); ++i) m[i] = n(i, 0);
            p = hasG ? pfqn::pfqn_joint_total(m, L, N, Z, G) : pfqn::pfqn_joint_total(m, L, N, Z);
        } else {
            p = hasG ? pfqn::pfqn_joint(n, L, N, Z, G) : pfqn::pfqn_joint(n, L, N, Z);
        }
        Json out;
        out["pjoint"] = encode_scalar(p);
        return out;
    }
};
 
struct OpJointmarg {
    template <class T>
    static Json run(Args& a) {
        const std::vector<int> n = a.ints("n");
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const std::vector<std::size_t> infset = indices_or_empty(a, "infset");
        const bool hasG = a.has("lGN");
        const T G = hasG ? exp_of_log<T>(a, "lGN") : num_traits<T>::from_int(0);
        const bool hasE = a.has("engine");
        const std::string engine = a.text("engine", "");
        a.done();
        T p;
        if (hasG)
            p = hasE ? pfqn::pfqn_jointmarg(n, L, N, infset, G, engine)
                     : pfqn::pfqn_jointmarg(n, L, N, infset, G);
        else
            p = hasE ? pfqn::pfqn_jointmarg(n, L, N, infset, engine)
                     : pfqn::pfqn_jointmarg(n, L, N, infset);
        Json out;
        out["pjoint"] = encode_scalar(p);
        if constexpr (num_traits<T>::has_transcendental) {
            using std::log;
            out["lpjoint"] = encode_scalar(T(log(p)));
        } else {
            out["lpjoint"] = std::log(num_traits<T>::to_double(p));
        }
        return out;
    }
};
 
struct OpLcfsqnCa {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> alpha = a.vector<T>("alpha");
        const std::vector<T> beta = a.vector<T>("beta");
        const std::vector<int> N = a.ints("N");
        a.done();
        const pfqn::LcfsQnResult<T> r = pfqn::pfqn_lcfsqn_ca(alpha, beta, N);
        Json out;
        out["G"] = encode_scalar(r.G);
        out["V"] = encode_scalar(r.V);
        return out;
    }
};
 
struct OpLcfsqnMva {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> alpha = a.vector<T>("alpha");
        const std::vector<T> beta = a.vector<T>("beta");
        const std::vector<int> N = a.ints("N");
        a.done();
        const pfqn::LcfsMvaResult<T> r = pfqn::pfqn_lcfsqn_mva(alpha, beta, N);
        Json out;
        out["T"] = encode_vector(r.T_);
        out["Q"] = encode_matrix(r.Q);
        out["U"] = encode_matrix(r.U);
        out["B"] = encode_matrix(r.B);
        return out;
    }
};
 
struct OpLcfsqnNc {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> alpha = a.vector<T>("alpha");
        const std::vector<T> beta = a.vector<T>("beta");
        const std::vector<int> N = a.ints("N");
        a.done();
        Json out;
        out["G"] = encode_scalar(pfqn::pfqn_lcfsqn_nc(alpha, beta, N));
        return out;
    }
};
 
struct OpLdmxEc {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> lambda = a.vector<T>("lambda");
        const Matrix<T> D = a.matrix<T>("D");
        const Matrix<T> mu = a.matrix<T>("mu");
        a.done();
        const pfqn::LdmxEcResult<T> r = pfqn::pfqn_ldmx_ec(lambda, D, mu);
        Json out;
        out["EC"] = encode_matrix(r.EC);
        out["E"] = encode_matrix(r.E);
        out["Eprime"] = encode_matrix(r.Eprime);
        out["Lo"] = encode_vector(r.Lo);
        return out;
    }
};
 
struct OpCbh {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> L = a.vector<T>("L");
        const int N = int_arg(a, "N");
        const T Z = a.number<T>("Z");
        const bool hasLevel = a.has("level");
        const int level = hasLevel ? int_arg(a, "level") : 0;
        a.done();
        const pfqn::CbhBounds<T> r = hasLevel ? pfqn::pfqn_cbh(L, N, Z, level) : pfqn::pfqn_cbh(L, N, Z);
        Json out;
        out["Xlo"] = encode_scalar(r.Xlo);
        out["Xhi"] = encode_scalar(r.Xhi);
        return out;
    }
};
 
struct OpCub {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        const bool lf = long_form(a, "pfqn_cub", {"order", "atol"}, {"order", "atol"});
        const int order = lf ? int_arg(a, "order") : 0;
        const T atol = lf ? a.number<T>("atol") : num_traits<T>::from_int(0);
        a.done();
        const pfqn::CubResult<T> r = lf ? pfqn::pfqn_cub(L, N, Z, order, atol) : pfqn::pfqn_cub(L, N, Z);
        return g_lg(r.G, r.lG, "Gn", "lGn");
    }
};
 
struct OpDac {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        const bool hasMu = a.has("mu");
        const Matrix<T> mu = a.matrix_or_empty<T>("mu");
        a.done();
        const pfqn::DacResult<T> r = hasMu ? pfqn::pfqn_dac(L, N, Z, mu) : pfqn::pfqn_dac(L, N, Z);
        Json out;
        out["Pjoint"] = encode_vector(r.Pjoint);
        Json st = Json::array();
        for (const std::vector<int>& s : r.states) st.push_back(encode_ints(s));
        out["states"] = st;
        out["XN"] = encode_vector(r.XN);
        out["QN"] = encode_matrix(r.QN);
        out["UN"] = encode_vector(r.UN);
        out["CN"] = encode_vector(r.CN);
        out["pi"] = encode_matrix(r.PI);
        return out;
    }
};
 
struct OpFnc {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> alpha = a.matrix<T>("alpha");
        const bool hasC = a.has("c");
        const std::vector<T> c = a.vector_or_empty<T>("c");
        a.done();
        const pfqn::FncResult<T> r = hasC ? pfqn::pfqn_fnc(alpha, c) : pfqn::pfqn_fnc(alpha);
        Json out;
        out["mu"] = encode_matrix(r.mu);
        out["c"] = encode_vector(r.c);
        return out;
    }
};
 
struct OpGrnmol {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        a.done();
        Json out;
        out["G"] = encode_scalar(pfqn::pfqn_grnmol(L, N));
        return out;
    }
};
 
/** The (L, N[, Z]) shape of the asymptotic expansions, returning G and lG. */
#define LINE_PFQN_A_LNZ(Name, fn, gk, lk)                                          \
    struct Name {                                                                  \
        template <class T>                                                         \
        static Json run(Args& a) {                                                 \
            const Matrix<T> L = a.matrix<T>("L");                                  \
            const std::vector<T> N = a.vector<T>("N");                             \
            const bool hasZ = a.has("Z");                                          \
            const std::vector<T> Z = a.vector_or_empty<T>("Z");                    \
            a.done();                                                              \
            const auto r = hasZ ? pfqn::fn(L, N, Z) : pfqn::fn(L, N);              \
            return g_lg(r.G, r.lG, gk, lk);                                        \
        }                                                                          \
    };
 
LINE_PFQN_A_LNZ(OpBkt, pfqn_bkt, "Gn", "lGn")
LINE_PFQN_A_LNZ(OpLe, pfqn_le, "Gn", "lGn")
LINE_PFQN_A_LNZ(OpBle, pfqn_ble, "Gn", "lGn")
#undef LINE_PFQN_A_LNZ
 
struct OpKt {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<T> N = a.vector<T>("N");
        const bool hasZ = a.has("Z");
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        a.done();
        const pfqn::KtResult<T> r = hasZ ? pfqn::pfqn_kt(L, N, Z) : pfqn::pfqn_kt(L, N);
        Json out = g_lg(r.G, r.lG);
        out["X"] = encode_vector(r.X);
        out["Q"] = encode_matrix(r.Q);
        return out;
    }
};
 
struct OpLekt {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<T> N = a.vector<T>("N");
        const bool hasZ = a.has("Z");
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        a.done();
        const pfqn::LektResult<T> r = hasZ ? pfqn::pfqn_lekt(L, N, Z) : pfqn::pfqn_lekt(L, N);
        Json out = g_lg(r.G, r.lG, "Gn", "lGn");
        out["route"] = r.route;
        return out;
    }
};
 
struct OpAghq {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<T> N = a.vector<T>("N");
        const std::size_t d = depth(a, "pfqn_aghq", {{"Z", true}, {"q", false}});
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        const std::size_t q = a.has("q") ? a.count("q") : 0;
        a.done();
        pfqn::LeResult<T> r;
        if (d == 0)
            r = pfqn::pfqn_aghq(L, N);
        else if (d == 1)
            r = pfqn::pfqn_aghq(L, N, Z);
        else
            r = pfqn::pfqn_aghq(L, N, Z, q);
        return g_lg(r.G, r.lG, "Gn", "lGn");
    }
};
 
struct OpLap {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> L = a.vector<T>("L");
        const std::vector<T> N = a.vector<T>("N");
        const std::vector<T> Z = a.vector<T>("Z");
        a.done();
        Json out;
        out["logI"] = encode_scalar(pfqn::pfqn_lap(L, N, Z));
        return out;
    }
};
 
struct OpLdbcmp {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> L = a.vector<T>("L");
        const T N = a.number<T>("N");
        const T Z = a.number<T>("Z");
        const bool lf = long_form(a, "pfqn_ldbcmp", {"c", "tol"}, {"c", "tol"});
        const std::vector<T> c = a.vector_or_empty<T>("c");
        const T tol = lf ? a.number<T>("tol") : num_traits<T>::from_int(0);
        a.done();
        const pfqn::LdBcmpBound<T> r = lf ? pfqn::pfqn_ldbcmp(L, N, Z, c, tol) : pfqn::pfqn_ldbcmp(L, N, Z);
        Json out;
        out["Xlo"] = encode_scalar(r.Xlo);
        out["Rhi"] = encode_scalar(r.Rhi);
        out["Qhat"] = encode_scalar(r.Qhat);
        out["applicable"] = r.applicable;
        return out;
    }
};
 
struct OpAbAmva {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> D = a.matrix<T>("D");
        const std::vector<int> N = a.ints("N");
        const Matrix<T> V = a.matrix<T>("V");
        const std::vector<int> nservers = a.ints("nservers");
        if (!a.has("sched")) (void)a.get("sched");  // required: the port has no default
        const std::vector<pfqn::SchedStrategy> sched = lin_sched(a, "sched");
        const bool lf = long_form(a, "pfqn_ab_amva", {"fcfsSchmidt", "marginalProbMethod"},
                                  {"fcfsSchmidt", "marginalProbMethod"});
        const bool schmidt = a.boolean("fcfsSchmidt", false);
        pfqn::AbMarginalMethod method = pfqn::AbMarginalMethod::Ab;
        if (lf) {
            std::string m = a.text("marginalProbMethod", "");
            for (char& ch : m) ch = static_cast<char>(std::tolower(static_cast<unsigned char>(ch)));
            if (m == "ab")
                method = pfqn::AbMarginalMethod::Ab;
            else if (m == "scat")
                method = pfqn::AbMarginalMethod::Scat;
            else
                throw InputError("pfqn_ab_amva: unknown marginalProbMethod '" + m +
                                 "'; expected \"ab\" or \"scat\"");
        }
        a.done();
        const pfqn::AbAmvaResult<T> r =
            lf ? pfqn::pfqn_ab_amva(D, N, V, nservers, sched, schmidt, method)
               : pfqn::pfqn_ab_amva(D, N, V, nservers, sched);
        Json out;
        out["QN"] = encode_matrix(r.QN);
        out["UN"] = encode_matrix(r.UN);
        out["RN"] = encode_matrix(r.RN);
        out["CN"] = encode_vector(r.CN);
        out["XN"] = encode_vector(r.XN);
        out["totiter"] = encode_count(r.totiter);
        return out;
    }
};
 
struct OpCubEvals {
    template <class T>
    static Json run(Args& a) {
        const int M = int_arg(a, "M");
        const int order = int_arg(a, "order");
        const bool hasZ = a.has("Z");
        const Matrix<double> Z = a.matrix_or_empty<double>("Z");
        a.done();
        double Zsum = 0.0;  // MATLAB reads sum(Z(:))
        for (std::size_t i = 0; i < Z.size(); ++i) Zsum += Z[i];
        Json out;
        out["nevals"] = hasZ ? pfqn::pfqn_cub_evals(M, order, Zsum) : pfqn::pfqn_cub_evals(M, order);
        return out;
    }
};
 
struct OpDnc {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> L = a.vector<T>("L");
        const T N = a.number<T>("N");
        a.done();
        const pfqn::DncResult<T> r = pfqn::pfqn_dnc(L, N);
        Json out = g_lg(r.G, r.lG);
        out["X"] = encode_scalar(r.X);
        return out;
    }
};
 
struct OpHst {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> L = a.vector<T>("L");
        const int N = int_arg(a, "N");
        const std::size_t d = depth(a, "pfqn_hst", {{"Z", false}, {"ist", false}});
        const T Z = a.has("Z") ? a.number<T>("Z") : num_traits<T>::from_int(0);
        const std::size_t ist = a.has("ist") ? a.count("ist") : 0;
        a.done();
        pfqn::HstResult<T> r;
        if (d == 0)
            r = pfqn::pfqn_hst(L, N);
        else if (d == 1)
            r = pfqn::pfqn_hst(L, N, Z);
        else
            r = pfqn::pfqn_hst(L, N, Z, ist);
        Json s;
        s["station"] = encode_count(r.station);
        s["X"] = encode_scalar(r.X);
        s["U"] = encode_scalar(r.U);
        s["Q"] = encode_scalar(r.Q);
        s["Pgeq"] = encode_vector(r.Pgeq);
        s["p"] = encode_vector(r.p);
        s["c"] = encode_vector(r.c);
        s["total"] = encode_scalar(r.total);
        s["worst"] = encode_scalar(r.worst);
        s["astar"] = encode_vector(r.astar);
        Json out;
        out["sens"] = s;
        return out;
    }
};
 
pfqn::ClustInner clust_inner(Args& a) {
    std::string m = a.text("inner", "");
    for (char& c : m) c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
    if (m == "lin") return pfqn::ClustInner::Linearizer;
    if (m == "bs") return pfqn::ClustInner::ProportionalEstimation;
    throw InputError("pfqn_clust: unknown inner method '" + m + "'; expected \"lin\" or \"bs\"");
}
 
struct OpClust {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<T> N = a.vector<T>("N");
        const bool hasZ = a.has("Z");
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        const std::size_t d = depth(a, "pfqn_clust",
                                    {{"subnets", true}, {"localclasses", true}, {"inner", false},
                                     {"tol", false}, {"maxiter", false}});
        std::vector<std::vector<std::size_t> > subnets, localclasses;
        if (const Json* j = a.opt("subnets")) subnets = index_lists(*j, "'subnets'", false);
        if (const Json* j = a.opt("localclasses"))
            localclasses = index_lists(*j, "'localclasses'", false);
        const pfqn::ClustInner inner = a.has("inner") ? clust_inner(a) : pfqn::ClustInner::Linearizer;
        const double tol = a.has("tol") ? a.number<double>("tol") : 0.0;
        const std::size_t maxiter = a.has("maxiter") ? a.count("maxiter") : 0;
        a.done();
        pfqn::AmvaResult<T> r;
        if (d == 0)
            r = hasZ ? pfqn::pfqn_clust(L, N, Z) : pfqn::pfqn_clust(L, N);
        else if (d <= 2)
            r = pfqn::pfqn_clust(L, N, Z, subnets, localclasses);
        else if (d == 3)
            r = pfqn::pfqn_clust(L, N, Z, subnets, localclasses, inner);
        else if (d == 4)
            r = pfqn::pfqn_clust(L, N, Z, subnets, localclasses, inner, tol);
        else
            r = pfqn::pfqn_clust(L, N, Z, subnets, localclasses, inner, tol, maxiter);
        return amva_results(r, "it");
    }
};
 
struct OpGerasimov {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const bool hasZ = a.has("Z");
        const Matrix<T> Z = a.matrix_or_empty<T>("Z");
        const Json* oj = a.opt("options");
        a.done();
        pfqn::NcResult<T> r;
        if (!oj) {
            r = hasZ ? pfqn::pfqn_gerasimov(L, N, Z) : pfqn::pfqn_gerasimov(L, N);
        } else {
            Args o(*oj, "pfqn_gerasimov options");
            const std::size_t d =
                depth(o, "pfqn_gerasimov options", {{"tol", false}, {"maxterms", false}});
            const double tol = o.has("tol") ? o.number<double>("tol") : 0.0;
            const std::size_t maxterms = o.has("maxterms") ? o.count("maxterms") : 0;
            o.done();
            if (d == 0)
                r = pfqn::pfqn_gerasimov(L, N, Z);
            else if (d == 1)
                r = pfqn::pfqn_gerasimov(L, N, Z, tol);
            else
                r = pfqn::pfqn_gerasimov(L, N, Z, tol, maxterms);
        }
        return g_lg(r.G, r.lG);
    }
};
 
struct OpLooping {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<T> N = a.vector<T>("N");
        const std::size_t d =
            depth(a, "pfqn_looping", {{"Z", true}, {"tol", false}, {"maxiter", false}});
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        const double tol = a.has("tol") ? a.number<double>("tol") : 0.0;
        const std::size_t maxiter = a.has("maxiter") ? a.count("maxiter") : 0;
        a.done();
        pfqn::LoopingBounds<T> r;
        if (d == 0)
            r = pfqn::pfqn_looping(L, N);
        else if (d == 1)
            r = pfqn::pfqn_looping(L, N, Z);
        else if (d == 2)
            r = pfqn::pfqn_looping(L, N, Z, tol);
        else
            r = pfqn::pfqn_looping(L, N, Z, tol, maxiter);
        Json out;
        out["Xlo"] = encode_vector(r.Xlo);
        out["Xup"] = encode_vector(r.Xup);
        out["QN"] = encode_matrix(r.Q);
        out["RN"] = encode_matrix(r.R);
        out["it"] = encode_count(r.iterations);
        out["converged"] = r.converged;
        return out;
    }
};
 
std::vector<std::size_t> busy_orders(Args& a) {
    const std::vector<int> v = a.ints("n");
    std::vector<std::size_t> out;
    for (int x : v) {
        if (x < 0) throw InputError("'n': a busy period order cannot be negative");
        out.push_back(static_cast<std::size_t>(x));
    }
    return out;
}
 
struct OpBusyp {
    static Json run(Args& a) {
        const std::vector<double> alpha = a.vector<double>("alpha");
        const Matrix<double> mu = a.matrix<double>("mu");
        const Matrix<double> P = a.matrix<double>("P");
        const double N = a.number<double>("N");
        const std::vector<std::size_t> subnet = indices(a, "subnet");
        const std::vector<std::size_t> n = busy_orders(a);
        const std::size_t d = depth(a, "pfqn_busyp", {{"gamma", true}, {"tol", false}});
        const std::vector<double> gamma = a.vector_or_empty<double>("gamma");
        const double tol = a.has("tol") ? a.number<double>("tol") : 0.0;
        a.done();
        pfqn::BusyPeriodResult r;
        if (d == 0)
            r = pfqn::pfqn_busyp(alpha, mu, P, N, subnet, n);
        else if (d == 1)
            r = pfqn::pfqn_busyp(alpha, mu, P, N, subnet, n, gamma);
        else
            r = pfqn::pfqn_busyp(alpha, mu, P, N, subnet, n, gamma, tol);
        Json out;
        out["b"] = encode_vector(r.b);
        out["lG"] = encode_vector(r.lG);
        out["lH"] = encode_vector(r.lH);
        return out;
    }
};
 
struct OpBusypClw {
    static Json run(Args& a) {
        const Matrix<double> alpha = a.matrix<double>("alpha");
        const Matrix<double> mu = a.matrix<double>("mu");
        const std::vector<Matrix<double> > P = a.matrices<double>("P");
        const std::vector<double> N = a.vector<double>("N");
        const std::vector<std::size_t> subnet = indices(a, "subnet");
        const std::vector<std::size_t> n = busy_orders(a);
        const std::size_t d = depth(a, "pfqn_busyp_clw",
                                    {{"gamma", true}, {"isdelay", true}, {"method", false}});
        const Matrix<double> gamma = a.matrix_or_empty<double>("gamma");
        std::vector<bool> isdelay;
        if (const Json* j = a.opt("isdelay")) isdelay = bool_vector(*j, "'isdelay'");
        const std::string method = a.text("method", "");
        a.done();
        Json out;
        if (d == 0)
            out["b"] = encode_vector(pfqn::pfqn_busyp_clw(alpha, mu, P, N, subnet, n));
        else if (d == 1)
            out["b"] = encode_vector(pfqn::pfqn_busyp_clw(alpha, mu, P, N, subnet, n, gamma));
        else if (d == 2)
            out["b"] = encode_vector(pfqn::pfqn_busyp_clw(alpha, mu, P, N, subnet, n, gamma, isdelay));
        else
            out["b"] = encode_vector(
                pfqn::pfqn_busyp_clw(alpha, mu, P, N, subnet, n, gamma, isdelay, method));
        return out;
    }
};
 
struct OpBusypMulticlass {
    static Json run(Args& a) {
        const Matrix<double> alpha = a.matrix<double>("alpha");
        const Matrix<double> mu = a.matrix<double>("mu");
        const std::vector<Matrix<double> > P = a.matrices<double>("P");
        const std::vector<double> N = a.vector<double>("N");
        const std::vector<std::size_t> subnet = indices(a, "subnet");
        const std::vector<std::size_t> n = busy_orders(a);
        const std::size_t d =
            depth(a, "pfqn_busyp_multiclass",
                  {{"gamma", true}, {"phi", true}, {"tol", false}, {"jobclass", false}});
        const Matrix<double> gamma = a.matrix_or_empty<double>("gamma");
        const Matrix<double> phi = a.matrix_or_empty<double>("phi");
        const double tol = a.has("tol") ? a.number<double>("tol") : 0.0;
        const int jobclass = a.has("jobclass") ? int_arg(a, "jobclass") : 0;
        a.done();
        std::vector<double> b;
        if (d == 0)
            b = pfqn::pfqn_busyp_multiclass(alpha, mu, P, N, subnet, n);
        else if (d == 1)
            b = pfqn::pfqn_busyp_multiclass(alpha, mu, P, N, subnet, n, gamma);
        else if (d == 2)
            b = pfqn::pfqn_busyp_multiclass(alpha, mu, P, N, subnet, n, gamma, phi);
        else if (d == 3)
            b = pfqn::pfqn_busyp_multiclass(alpha, mu, P, N, subnet, n, gamma, phi, tol);
        else
            b = pfqn::pfqn_busyp_multiclass(alpha, mu, P, N, subnet, n, gamma, phi, tol, jobclass);
        Json out;
        out["b"] = encode_vector(b);
        return out;
    }
};
 
struct OpBk {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<T> N = a.vector<T>("N");
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        a.done();
        const pfqn::BkResult<T> r = pfqn::pfqn_bk(L, N, Z);
        Json out = g_lg(r.G, r.lG);
        out["X"] = encode_vector(r.X);
        out["U"] = encode_matrix(r.U);
        out["A"] = encode_index_list(r.A);
        out["B"] = encode_index_list(r.B);
        return out;
    }
};
 
struct OpClw {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const std::vector<T> Z = a.vector<T>("Z");
        const bool hasM = a.has("m");
        std::vector<long> m;
        for (int x : a.ints_or_empty("m")) m.push_back(static_cast<long>(x));
        const Json* oj = a.opt("options");
        a.done();
        pfqn::ClwResult<T> r;
        if (oj) {
            // Every field the caller leaves out stays at the port's ClwOptions default.
            pfqn::ClwOptions opt;
            Args o(*oj, "pfqn_clw options");
            if (o.has("l")) opt.l = o.ints("l");
            if (o.has("gamma")) opt.gamma = o.vector<double>("gamma");
            if (o.has("euler")) opt.euler = o.boolean("euler", opt.euler);
            if (o.has("euler_n")) opt.eulerN = int_arg(o, "euler_n");
            if (o.has("euler_m")) opt.eulerM = int_arg(o, "euler_m");
            if (o.has("euler_tol")) opt.eulerTol = o.number<double>("euler_tol");
            if (o.has("euler_maxm")) opt.eulerMaxM = int_arg(o, "euler_maxm");
            if (o.has("dimred")) opt.dimred = o.boolean("dimred", opt.dimred);
            if (o.has("dimred_maxd")) opt.dimredMaxD = int_arg(o, "dimred_maxd");
            if (o.has("beta")) opt.beta = o.vector<double>("beta");
            o.done();
            r = pfqn::pfqn_clw(L, N, Z, m, opt);
        } else {
            r = hasM ? pfqn::pfqn_clw(L, N, Z, m) : pfqn::pfqn_clw(L, N, Z);
        }
        return g_lg(r.G, r.lG);
    }
};
 
struct OpHarelBounds {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> rho = a.vector<T>("rho");
        const int N = int_arg(a, "N");
        const bool lf = a.has("Z") || a.has("maxUB");
        const T Z = a.has("Z") ? a.number<T>("Z") : num_traits<T>::from_int(0);
        // maxUB <= 0 is the port's documented "not given" (min(N,7)), as [] is in MATLAB.
        const int maxUB = a.has("maxUB") ? int_arg(a, "maxUB") : 0;
        a.done();
        if (lf && !a.has("Z"))
            throw InputError("pfqn_harel_bounds: 'maxUB' is given but 'Z' is not; the C++ port "
                             "takes 'Z' positionally before 'maxUB'. Give 'Z' as well.");
        const pfqn::HarelBoundsResult<T> r =
            lf ? pfqn::pfqn_harel_bounds(rho, N, Z, maxUB) : pfqn::pfqn_harel_bounds(rho, N);
        // The port indexes UB and TH by population, leaving slot 0 unused; MATLAB's
        // UB(n) and TH(n) run over n = 1..maxUB, so slot 0 is dropped.
        const std::vector<T> UB(r.UB.begin() + 1, r.UB.end());
        const std::vector<T> TH(r.TH.begin() + 1, r.TH.end());
        Json out;
        out["LB"] = encode_scalar(r.LB);
        out["UB"] = encode_vector(UB);
        out["TH"] = encode_vector(TH);
        return out;
    }
};
 
 
// ---------------------------------------------------------------------------
// Seeded Monte Carlo estimators. Each port draws from a caller-owned
// pfqn::McRng; the JSON form names it by a required seed (required_seed).
// ---------------------------------------------------------------------------
 
/**
 * pfqn_cftp(L, N, S, nsamples, method, seed): MATLAB's [Q, X, T]. S lists the
 * servers per station, Inf for a delay; method is 'cftp' (exact coupling from
 * the past) or 'approx'.
 */
struct OpCftp {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> L = a.vector<T>("L");
        const int N = int_arg(a, "N");
        std::vector<int> S;
        if (const Json* sj = a.opt("S")) {
            const std::vector<double> sd = vector_from_json<double>(*sj, "pfqn_cftp: S");
            for (double v : sd) {
                if (std::isinf(v) && v > 0) {
                    S.push_back(pfqn::cftp_inf_servers);
                } else if (!(v >= 1) || v != std::floor(v)) {
                    throw InputError("pfqn_cftp: S must hold positive integer server counts or "
                                     "Inf, got " + shortest_decimal(v));
                } else {
                    S.push_back(static_cast<int>(v));
                }
            }
        }
        const std::size_t nsamples = a.count("nsamples", 1);
        const std::string m = a.text("method", "cftp");
        pfqn::CftpMethod method;
        if (m == "cftp") method = pfqn::CftpMethod::Cftp;
        else if (m == "approx") method = pfqn::CftpMethod::Approx;
        else throw InputError("pfqn_cftp: method must be 'cftp' or 'approx', got '" + m + "'");
        pfqn::McRng rng(required_seed(a, "pfqn_cftp"));
        a.done();
        const pfqn::CftpResult<T> r = pfqn::pfqn_cftp(L, N, S, nsamples, method, rng);
        Matrix<double> X(r.X.rows(), r.X.cols());
        for (std::size_t i = 0; i < X.rows(); ++i)
            for (std::size_t j = 0; j < X.cols(); ++j) X(i, j) = r.X(i, j);
        std::vector<double> horizon(r.horizon.begin(), r.horizon.end());
        Json out;
        out["Q"] = encode_vector(r.Q);
        out["X"] = encode_matrix(X);
        out["T"] = encode_vector(horizon);
        return out;
    }
};
 
/**
 * The options object of pfqn_is and pfqn_ld_is: {samples (1e4), seed}. The
 * seed is required there as everywhere else in this family.
 */
struct IsOptions {
    std::size_t samples = 10000;
    std::uint64_t seed = 0;
};
 
inline IsOptions read_is_options(Args& a, const std::string& fn) {
    const Json* oj = a.opt("options");
    if (!oj)
        throw InputError(fn + ": options.seed is required; the port draws from a seeded "
                              "generator and has no session RNG to fall back on");
    Args op(*oj, fn + " options");
    IsOptions o;
    o.samples = op.count("samples", o.samples);
    o.seed = required_seed(op, fn, "options.seed");
    op.done();
    return o;
}
 
/** pfqn_is(L, N, Z, options): importance-sampling estimate of [G, lG]. */
struct OpIs {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        const IsOptions o = read_is_options(a, "pfqn_is");
        a.done();
        pfqn::McRng rng(o.seed);
        const pfqn::NcResult<T> r = pfqn::pfqn_is(L, N, Z, o.samples, rng);
        return g_lg(r.G, r.lG);
    }
};
 
/** pfqn_ld_is(L, N, Z, mu, options): the load-dependent estimator. */
struct OpLdIs {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<int> N = a.ints("N");
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        const Matrix<T> mu = a.matrix_or_empty<T>("mu");
        const IsOptions o = read_is_options(a, "pfqn_ld_is");
        a.done();
        pfqn::McRng rng(o.seed);
        const pfqn::NcResult<T> r = pfqn::pfqn_ld_is(L, N, Z, mu, o.samples, rng);
        return g_lg(r.G, r.lG);
    }
};
 
/** pfqn_ls(L, N, Z, I, seed): logistic-sampling estimate, MATLAB's [Gn, lGn]. */
struct OpLs {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> L = a.matrix<T>("L");
        const std::vector<T> N = a.vector<T>("N");
        const std::vector<T> Z = a.vector_or_empty<T>("Z");
        const std::size_t I = a.count("I", 100000);
        pfqn::McRng rng(required_seed(a, "pfqn_ls"));
        a.done();
        const pfqn::LsResult<T> r = pfqn::pfqn_ls(L, N, Z, I, rng);
        return g_lg(r.G, r.lG, "Gn", "lGn");
    }
};
 
}  // namespace
 
void register_api_pfqn_a(DispatchTable& t) {
    // double only
    t["pfqn_cyclet_ofree"] = make<DoubleOnly<OpCycletOfree> >();
    t["pfqn_busyp"] = make<DoubleOnly<OpBusyp> >();
    t["pfqn_busyp_clw"] = make<DoubleOnly<OpBusypClw> >();
    t["pfqn_busyp_multiclass"] = make<DoubleOnly<OpBusypMulticlass> >();
    // double, exact and real
    t["pfqn_expand"] = make<OpExpand>();
    // The registry lists exact, but the port static_asserts transcendental arithmetic.
    t["pfqn_lldfun"] = make<NeedsTranscendental<OpLldfun> >();
    t["pfqn_comomrm_ld"] = make<OpComomrmLd>();
    t["pfqn_gldsingle"] = make<OpGldsingle>();
    t["pfqn_lld"] = make<OpLld>();
    t["pfqn_joint"] = make<OpJoint>();
    t["pfqn_lcfsqn_ca"] = make<OpLcfsqnCa>();
    t["pfqn_lcfsqn_mva"] = make<OpLcfsqnMva>();
    t["pfqn_lcfsqn_nc"] = make<OpLcfsqnNc>();
    t["pfqn_ldmx_ec"] = make<OpLdmxEc>();
    t["pfqn_lldsingle"] = make<OpLldsingle>();
    t["pfqn_cbh"] = make<OpCbh>();
    t["pfqn_dac"] = make<OpDac>();
    t["pfqn_fnc"] = make<OpFnc>();
    t["pfqn_grnmol"] = make<OpGrnmol>();
    t["pfqn_ldbcmp"] = make<OpLdbcmp>();
    t["pfqn_cub_evals"] = make<OpCubEvals>();
    t["pfqn_gerasimov"] = make<OpGerasimov>();
    t["pfqn_jointmarg"] = make<OpJointmarg>();
    // double and real
    t["pfqn_explicit"] = make<NeedsTranscendental<OpExplicit> >();
    t["pfqn_explicit_ld"] = make<NeedsTranscendental<OpExplicitLd> >();
    t["pfqn_divdiff_ld"] = make<NeedsTranscendental<OpDivdiffLd> >();
    t["pfqn_bs"] = make<NeedsTranscendental<OpBs> >();
    t["pfqn_cntol"] = make<NeedsTranscendental<OpCntol> >();
    t["pfqn_aql"] = make<NeedsTranscendental<OpAql> >();
    t["pfqn_linearizer"] = make<NeedsTranscendental<OpLinearizer> >();
    t["pfqn_linearizerms"] = make<NeedsTranscendental<OpLinearizerms> >();
    t["pfqn_linearizermx"] = make<NeedsTranscendental<OpLinearizermx> >();
    t["pfqn_gflinearizer"] = make<NeedsTranscendental<OpGflinearizer> >();
    t["pfqn_egflinearizer"] = make<NeedsTranscendental<OpEgflinearizer> >();
    t["pfqn_conwayms"] = make<NeedsTranscendental<OpConwayms> >();
    t["pfqn_cub"] = make<NeedsTranscendental<OpCub> >();
    t["pfqn_kt"] = make<NeedsTranscendental<OpKt> >();
    t["pfqn_bkt"] = make<NeedsTranscendental<OpBkt> >();
    t["pfqn_lekt"] = make<NeedsTranscendental<OpLekt> >();
    t["pfqn_lap"] = make<NeedsTranscendental<OpLap> >();
    t["pfqn_le"] = make<NeedsTranscendental<OpLe> >();
    t["pfqn_ble"] = make<NeedsTranscendental<OpBle> >();
    t["pfqn_aghq"] = make<NeedsTranscendental<OpAghq> >();
    t["pfqn_ab_amva"] = make<NeedsTranscendental<OpAbAmva> >();
    t["pfqn_dnc"] = make<NeedsTranscendental<OpDnc> >();
    t["pfqn_hst"] = make<NeedsTranscendental<OpHst> >();
    t["pfqn_chow"] = make<NeedsTranscendental<OpChow> >();
    t["pfqn_clust"] = make<NeedsTranscendental<OpClust> >();
    t["pfqn_dmlin"] = make<NeedsTranscendental<OpDmlin> >();
    t["pfqn_lcp"] = make<NeedsTranscendental<OpLcp> >();
    t["pfqn_looping"] = make<NeedsTranscendental<OpLooping> >();
    t["pfqn_bk"] = make<NeedsTranscendental<OpBk> >();
    t["pfqn_clw"] = make<NeedsTranscendental<OpClw> >();
    t["pfqn_harel_bounds"] = make<NeedsTranscendental<OpHarelBounds> >();
    t["pfqn_cftp"] = make<NeedsTranscendental<OpCftp> >();
    t["pfqn_is"] = make<NeedsTranscendental<OpIs> >();
    t["pfqn_ld_is"] = make<NeedsTranscendental<OpLdIs> >();
    t["pfqn_ls"] = make<NeedsTranscendental<OpLs> >();
}
 
void register_library_only_pfqn_a(LibraryOnlyTable& t) {
    const std::string callables =
        "its scaling argument is a vector of callables (one std::function per station), which "
        "has no JSON representation.";
    t["pfqn_cdfun"] = "pfqn_cdfun: " + callables;
    t["pfqn_jdfun"] = "pfqn_jdfun: " + callables;
    const std::string oirate =
        "its rate argument mu is a vector of OiRate callables (std::function of the occupancy "
        "vector), which has no JSON representation.";
    t["pfqn_clwjd"] = "pfqn_clwjd: " + oirate;
    t["pfqn_clwoi"] = "pfqn_clwoi: " + oirate;
}
 
}  // namespace reg
}  // namespace line