api_dispatch_trace.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 'trace' group (see api_dispatch_detail.h):
 * the trace, mtrace and infer domains.
 *
 * ARGUMENT AND RESULT NAMES are the MATLAB ones: a trace is a JSON vector, a
 * marked trace a vector plus a label vector, and the results carry the MATLAB
 * output names. Where the port leaves a statistic UNDEFINED (a class with no
 * sample, which MATLAB reports as 0/0 = NaN) the entry is written as JSON null
 * at every arithmetic, rather than as whatever the port left in the slot.
 *
 * POSITIONAL OPTIONALS. A MATLAB positional optional maps onto a C++ default
 * argument. An absent optional calls the shorter overload, so no default is
 * restated here; a later optional given without the earlier ones is refused,
 * because reaching it would mean inventing the value of the earlier one.
 */
 
#include "api_dispatch_detail.h"
 
#include <array>
#include <cctype>
#include <cmath>
#include <cstddef>
#include <functional>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/infer/infer_compute_ql_at_arrival.h"
#include "line/api/infer/infer_get_qlen_arrival.h"
#include "line/api/infer/infer_gibbs.h"
#include "line/api/infer/infer_lqn_findbyname.h"
#include "line/api/infer/infer_lqn_getobs.h"
#include "line/api/infer/infer_minps_setup.h"
#include "line/api/infer/infer_mlps.h"
#include "line/api/infer/infer_nhpp_ks.h"
#include "line/api/infer/infer_qmle.h"
#include "line/api/infer/infer_quick_model.h"
#include "line/api/infer/infer_rps.h"
#include "line/api/infer/infer_smp.h"
#include "line/api/infer/infer_variational.h"
#include "line/api/trace/mtrace_backward_moment.h"
#include "line/api/trace/mtrace_count.h"
#include "line/api/trace/mtrace_cov.h"
#include "line/api/trace/mtrace_cross_moment.h"
#include "line/api/trace/mtrace_forward_moment.h"
#include "line/api/trace/mtrace_iat2counts.h"
#include "line/api/trace/mtrace_joint.h"
#include "line/api/trace/mtrace_mean.h"
#include "line/api/trace/mtrace_merge.h"
#include "line/api/trace/mtrace_moment.h"
#include "line/api/trace/mtrace_moment_simple.h"
#include "line/api/trace/mtrace_pc.h"
#include "line/api/trace/mtrace_sigma.h"
#include "line/api/trace/mtrace_sigma2.h"
#include "line/api/trace/mtrace_split.h"
#include "line/api/trace/mtrace_summary.h"
#include "line/api/trace/trace_acf.h"
#include "line/api/trace/trace_bicov.h"
#include "line/api/trace/trace_gamma.h"
#include "line/api/trace/trace_iat2bins.h"
#include "line/api/trace/trace_iat2counts.h"
#include "line/api/trace/trace_idc.h"
#include "line/api/trace/trace_idi.h"
#include "line/api/trace/trace_joint.h"
#include "line/api/trace/trace_mean.h"
#include "line/api/trace/trace_pmf.h"
#include "line/api/trace/trace_scv.h"
#include "line/api/trace/trace_skew.h"
#include "line/api/trace/trace_summary.h"
#include "line/api/trace/trace_var.h"
#include "line/io/network_reader.h"
#include "line/io/network_writer.h"
 
namespace line {
namespace reg {
 
namespace {
 
// ---------------------------------------------------------------------------
// Shared readers and encoders
// ---------------------------------------------------------------------------
 
Json encode_longs(const std::vector<long>& v) {
    Json a = Json::array();
    for (long x : v) a.push_back(static_cast<long long>(x));
    return a;
}
 
Json encode_long_matrix(const Matrix<long>& m) {
    Json a = Json::array();
    for (std::size_t i = 0; i < m.rows(); ++i) {
        Json row = Json::array();
        for (std::size_t k = 0; k < m.cols(); ++k) row.push_back(static_cast<long long>(m(i, k)));
        a.push_back(row);
    }
    return a;
}
 
Json encode_sizes(const std::vector<std::size_t>& v) {
    Json a = Json::array();
    for (std::size_t x : v) a.push_back(encode_count(x));
    return a;
}
 
/** An undefined statistic (no sample behind it) is null, at every arithmetic. */
template <class T>
Json encode_defined(const T& v, long count) {
    return count == 0 ? Json(nullptr) : encode_scalar(v);
}
 
template <class T>
Json encode_vector_defined(const std::vector<T>& v, const std::vector<long>& count) {
    Json a = Json::array();
    for (std::size_t k = 0; k < v.size(); ++k) a.push_back(encode_defined(v[k], count[k]));
    return a;
}
 
template <class T>
Json encode_matrix_defined(const Matrix<T>& m, const Matrix<long>& count) {
    Json a = Json::array();
    for (std::size_t i = 0; i < m.rows(); ++i) {
        Json row = Json::array();
        for (std::size_t k = 0; k < m.cols(); ++k) row.push_back(encode_defined(m(i, k), count(i, k)));
        a.push_back(row);
    }
    return a;
}
 
std::string where_of(const std::string& fn, const char* key) { return fn + ": " + key; }
 
std::vector<long> to_longs(const std::vector<int>& v) { return std::vector<long>(v.begin(), v.end()); }
 
/** Nonnegative integers, e.g. moment orders or 0-based-or-sink state labels. */
std::vector<std::size_t> nonneg_sizes(const std::vector<int>& v, const std::string& where) {
    std::vector<std::size_t> out;
    out.reserve(v.size());
    for (int x : v) {
        if (x < 0) throw InputError(where + ": entries must be nonnegative, got " + std::to_string(x));
        out.push_back(static_cast<std::size_t>(x));
    }
    return out;
}
 
std::vector<unsigned> nonneg_unsigned(const std::vector<int>& v, const std::string& where) {
    const std::vector<std::size_t> s = nonneg_sizes(v, where);
    return std::vector<unsigned>(s.begin(), s.end());
}
 
/** MATLAB 1-based class labels to the port's 0-based ones. */
std::vector<std::size_t> one_based_to_zero(const std::vector<int>& v, const std::string& where) {
    std::vector<std::size_t> out;
    out.reserve(v.size());
    for (int x : v) {
        if (x < 1) throw InputError(where + ": class labels are 1-based, got " + std::to_string(x));
        out.push_back(static_cast<std::size_t>(x - 1));
    }
    return out;
}
 
std::vector<std::size_t> one_based_sizes(const std::vector<int>& v, const std::string& where) {
    std::vector<std::size_t> out;
    out.reserve(v.size());
    for (int x : v) {
        if (x < 1) throw InputError(where + ": class labels are 1-based, got " + std::to_string(x));
        out.push_back(static_cast<std::size_t>(x));
    }
    return out;
}
 
/** A number where JSON null stands for `fill` (MATLAB's NaN or Inf placeholder). */
template <class T>
T number_or_fill(const Json& j, const std::string& where, const T& fill) {
    return j.is_null() ? fill : number_from_json<T>(j, where);
}
 
/**
 * A matrix like matrix_from_json, but with JSON null read as `fill`. The JSON
 * encoder writes a double NaN as null, so this is how a NaN placeholder
 * (MATLAB's 'not observed', 'not applicable') crosses the boundary.
 */
template <class T>
Matrix<T> matrix_or_fill(const Json& j, const std::string& where, const T& fill) {
    if (!j.is_array()) {
        Matrix<T> m(1, 1);
        m(0, 0) = number_or_fill<T>(j, where, fill);
        return m;
    }
    if (j.empty()) return Matrix<T>();
    if (j[0].is_array()) {
        const std::size_t rows = j.size(), cols = j[0].size();
        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) + " is not an array of " +
                                 std::to_string(cols) + " entries");
        if (cols == 0) return Matrix<T>();
        Matrix<T> m(rows, cols);
        for (std::size_t i = 0; i < rows; ++i)
            for (std::size_t k = 0; k < cols; ++k) m(i, k) = number_or_fill<T>(j[i][k], where, fill);
        return m;
    }
    Matrix<T> m(1, j.size());
    for (std::size_t k = 0; k < j.size(); ++k) m(0, k) = number_or_fill<T>(j[k], where, fill);
    return m;
}
 
/** MATLAB's reshape(X, 1, numel(X)): column-major flattening of any shape. */
template <class T>
std::vector<T> column_major(const Matrix<T>& m) {
    std::vector<T> v;
    v.reserve(m.rows() * m.cols());
    for (std::size_t c = 0; c < m.cols(); ++c)
        for (std::size_t r = 0; r < m.rows(); ++r) v.push_back(m(r, c));
    return v;
}
 
std::string lower(std::string s) {
    for (char& c : s) c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
    return s;
}
 
std::string string_of(const Json& j, const std::string& where) {
    if (!j.is_string())
        throw InputError(where + ": expected a string, got " + std::string(j.type_name()));
    return j.get<std::string>();
}
 
std::vector<std::string> strings_of(const Json& j, const std::string& where) {
    if (!j.is_array()) throw InputError(where + ": expected a list of strings");
    std::vector<std::string> out;
    for (const Json& e : j) out.push_back(string_of(e, where));
    return out;
}
 
/** A flag stored as a JSON boolean or 0/1. */
bool flag_of(const Json& j, const std::string& where) {
    if (j.is_boolean()) return j.get<bool>();
    const std::vector<int> v = int_vector_from_json(j, where);
    if (v.size() != 1 || (v[0] != 0 && v[0] != 1))
        throw InputError(where + ": expected true, false, 0 or 1");
    return v[0] != 0;
}
 
/**
 * The MATLAB 6 x (K+1) cell array `data` of the inference traces, of which the
 * references read row 3 (arrival times, ms) and row 4 (response times) of the
 * first K columns. Row-major JSON: data[row][column], each cell a vector.
 */
template <class T>
void read_trace_cells(const Json& data, const std::string& where, std::vector<std::vector<T>>& at_ms,
                      std::vector<std::vector<T>>& rt) {
    if (!data.is_array() || data.size() < 4)
        throw InputError(where + ": expected the 6 x (K+1) cell array, as a list of rows; rows 3 "
                                 "(arrival times, ms) and 4 (response times) are read");
    const std::size_t cols = data[0].is_array() ? data[0].size() : 0;
    for (std::size_t r = 0; r < data.size(); ++r)
        if (!data[r].is_array() || data[r].size() != cols)
            throw InputError(where + ": row " + std::to_string(r + 1) + " must have " +
                             std::to_string(cols) + " cells, like row 1");
    if (cols < 2) throw InputError(where + ": at least one class column plus the trailing column");
    const std::size_t K = cols - 1;
    at_ms.assign(K, std::vector<T>());
    rt.assign(K, std::vector<T>());
    for (std::size_t k = 0; k < K; ++k) {
        at_ms[k] = vector_from_json<T>(data[2][k], where + "{3," + std::to_string(k + 1) + "}");
        rt[k] = vector_from_json<T>(data[3][k], where + "{4," + std::to_string(k + 1) + "}");
    }
}
 
// ---------------------------------------------------------------------------
// trace: one trace in, one scalar out
// ---------------------------------------------------------------------------
 
#define LINE_TRACE_SCALAR_OP(Name, fn, outkey)                          \
    struct Name {                                                       \
        template <class T>                                              \
        static Json run(Args& a) {                                      \
            const std::vector<T> S = a.vector<T>("S");                  \
            a.done();                                                   \
            Json out;                                                   \
            out[outkey] = encode_scalar(line::trace::fn<T>(S));         \
            return out;                                                 \
        }                                                               \
    };
 
LINE_TRACE_SCALAR_OP(OpTraceMean, trace_mean, "m")
LINE_TRACE_SCALAR_OP(OpTraceVar, trace_var, "v")
LINE_TRACE_SCALAR_OP(OpTraceScv, trace_scv, "scv")
LINE_TRACE_SCALAR_OP(OpTraceSkew, trace_skew, "skew")
LINE_TRACE_SCALAR_OP(OpTraceIdc, trace_idc, "IDC")
#undef LINE_TRACE_SCALAR_OP
 
struct OpTraceAcf {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> S = a.vector<T>("S");
        const bool has_lags = a.has("lags");
        const std::vector<int> lags = a.ints_or_empty("lags");
        a.done();
        Json out;
        out["rho"] = encode_vector(has_lags ? line::trace::trace_acf<T>(S, lags)
                                            : line::trace::trace_acf<T>(S));
        return out;
    }
};
 
struct OpTraceBicov {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> S = a.vector<T>("S");
        const std::vector<int> grid = a.ints("GRID");
        a.done();
        const line::trace::TraceBicovResult<T> r = line::trace::trace_bicov<T>(S, grid);
        Json lags = Json::array();
        for (const std::vector<int>& l : r.lags) lags.push_back(encode_ints(l));
        Json out;
        out["BiCov"] = encode_vector(r.bicov);
        out["BiCovLags"] = lags;
        return out;
    }
};
 
struct OpTraceGamma {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> S = a.vector<T>("T");
        const bool has_limit = a.has("limit");
        const int limit = has_limit ? a.required_integer("limit") : 0;
        a.opt("limit");
        a.done();
        const line::trace::TraceGammaResult<T> r =
            has_limit ? line::trace::trace_gamma<T>(S, static_cast<long>(limit))
                      : line::trace::trace_gamma<T>(S);
        Json out;
        out["GAMMA"] = encode_scalar(r.gamma);
        out["RHO0"] = encode_scalar(r.rho0);
        out["RESIDUALS"] = encode_scalar(r.residuals);
        return out;
    }
};
 
struct OpTraceIat2bins {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> S = a.vector<T>("S");
        const T scale = a.number<T>("scale");
        a.done();
        const line::trace::TraceBinsResult r = line::trace::trace_iat2bins<T>(S, scale);
        Json out;
        out["C"] = encode_longs(r.counts);
        out["bC"] = encode_longs(r.membership);
        return out;
    }
};
 
struct OpTraceIat2counts {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> S = a.vector<T>("S");
        const T scale = a.number<T>("scale");
        a.done();
        Json out;
        out["C"] = encode_longs(line::trace::trace_iat2counts<T>(S, scale));
        return out;
    }
};
 
/**
 * trace_idi(S, kset, OPTION, n). OPTION 'aggregate' with n is the port's
 * aggregate_n; 'aggregate-mix' is not ported (see trace_idi.h). The kset
 * default of the reference (min(1000, ceil(n/30))) has no overload in the port,
 * so kset is required rather than restated here.
 */
struct OpTraceIdi {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> S = a.vector<T>("S");
        const std::vector<long> kset = to_longs(a.ints("kset"));
        const bool has_option = a.has("OPTION");
        const std::string option = lower(a.text("OPTION", std::string()));
        const bool has_n = a.has("n");
        a.opt("n");
        line::trace::TraceIdiResult<T> r;
        if (!has_option) {
            if (has_n)
                throw InputError("trace_idi: 'n' is the aggregation count of OPTION 'aggregate'; "
                                 "give OPTION too");
            a.done();
            r = line::trace::trace_idi<T>(S, kset);
        } else if (option == "aggregate") {
            if (!has_n) throw InputError("trace_idi: OPTION 'aggregate' needs 'n'");
            const int n = a.required_integer("n");
            if (n < 1) throw InputError("trace_idi: 'n' must be a positive integer");
            a.done();
            r = line::trace::trace_idi<T>(S, kset, static_cast<long>(n));
        } else if (option == "aggregate-mix") {
            throw UnsupportedError(ErrorKind{"JsonApi"},
                "trace_idi: OPTION 'aggregate-mix' is not ported; the reference branch keeps "
                "only the statistic of its last partition (see trace_idi.h)");
        } else {
            throw InputError("trace_idi: OPTION must be 'aggregate', got '" + option + "'");
        }
        Json out;
        out["IDIk"] = encode_vector(r.idi);
        out["support"] = encode_longs(r.support);
        return out;
    }
};
 
struct OpTraceJoint {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> S = a.vector<T>("S");
        const std::vector<int> lag = a.ints("lag");
        const std::vector<unsigned> order = nonneg_unsigned(a.ints("order"), "trace_joint: order");
        a.done();
        Json out;
        out["JM"] = encode_scalar(line::trace::trace_joint<T>(S, lag, order));
        return out;
    }
};
 
struct OpTracePmf {
    template <class T>
    static Json run(Args& a) {
        const std::vector<int> X = a.ints("X");
        a.done();
        const line::trace::TracePmfResult<T> r = line::trace::trace_pmf<T>(X);
        Json out;
        out["pmf"] = encode_vector(r.pmf);
        out["px"] = encode_ints(r.values);
        return out;
    }
};
 
/**
 * trace_summary(m). The port's skewness is bias-corrected and its percentiles
 * are type 7, its ACF covers lags 1..4 and its tail is the 95th percentile
 * alone; those are the values reported (see trace_summary.h). KURT is MATLAB's
 * kurtosis, i.e. the port's excess kurtosis plus 3.
 */
struct OpTraceSummary {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> S = a.vector<T>("m");
        a.unsupported("fid", "it is a MATLAB file identifier for the printed report; the API "
                             "returns the values instead of printing them");
        a.done();
        const line::trace::TraceSummary<T> r = line::trace::trace_summary<T>(S);
        Json out;
        out["MEAN"] = encode_scalar(r.mean);
        out["SCV"] = encode_scalar(r.scv);
        out["MAD"] = encode_scalar(r.mad);
        out["SKEW"] = encode_scalar(r.skew);
        out["KURT"] = encode_scalar(T(r.kurt_excess + num_traits<T>::from_int(3)));
        out["QUART"] = encode_vector(std::vector<T>{r.q25, r.q50, r.q75});
        out["TAILS1PERC"] = encode_vector(std::vector<T>{r.p95});
        out["MINMAX"] = encode_vector(std::vector<T>{r.min, r.max});
        out["IQR"] = encode_scalar(r.iqr);
        out["ACF"] = encode_vector(r.acf);
        out["IDC"] = encode_scalar(r.idc);
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// mtrace: marked traces
// ---------------------------------------------------------------------------
 
struct OpMtraceMoment {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> Tv = a.vector<T>("T");
        const std::vector<int> A = a.ints("A");
        const std::vector<unsigned> orders = nonneg_unsigned(a.ints("ORDERS"), "mtrace_moment: ORDERS");
        const bool has_after = a.has("AFTER"), has_norm = a.has("NORM");
        if (has_norm && !has_after)
            throw InputError("mtrace_moment: NORM follows AFTER positionally; give AFTER too");
        const bool after = has_after ? a.boolean("AFTER", false) : false;
        const bool norm = has_norm ? a.boolean("NORM", false) : false;
        a.opt("AFTER");
        a.opt("NORM");
        a.done();
        Matrix<T> M;
        if (has_norm)
            M = line::trace::mtrace_moment<T>(Tv, A, orders, after, norm);
        else if (has_after)
            M = line::trace::mtrace_moment<T>(Tv, A, orders, after);
        else
            M = line::trace::mtrace_moment<T>(Tv, A, orders);
        Json out;
        out["M"] = encode_matrix(M);
        return out;
    }
};
 
#define LINE_MTRACE_DIRECTED_MOMENT_OP(Name, fn)                                          \
    struct Name {                                                                         \
        template <class T>                                                                \
        static Json run(Args& a) {                                                        \
            const std::vector<T> Tv = a.vector<T>("T");                                   \
            const std::vector<int> A = a.ints("A");                                       \
            const std::vector<unsigned> orders =                                          \
                nonneg_unsigned(a.ints("ORDERS"), #fn ": ORDERS");                        \
            const bool has_norm = a.has("NORM");                                          \
            const bool norm = has_norm ? a.boolean("NORM", false) : false;                \
            a.opt("NORM");                                                                \
            a.done();                                                                     \
            Json out;                                                                     \
            out["M"] = encode_matrix(has_norm ? line::trace::fn<T>(Tv, A, orders, norm)   \
                                              : line::trace::fn<T>(Tv, A, orders));       \
            return out;                                                                   \
        }                                                                                 \
    };
 
LINE_MTRACE_DIRECTED_MOMENT_OP(OpMtraceForwardMoment, mtrace_forward_moment)
LINE_MTRACE_DIRECTED_MOMENT_OP(OpMtraceBackwardMoment, mtrace_backward_moment)
#undef LINE_MTRACE_DIRECTED_MOMENT_OP
 
#define LINE_MTRACE_CROSS_MOMENT_OP(Name, fn)                                      \
    struct Name {                                                                  \
        template <class T>                                                         \
        static Json run(Args& a) {                                                 \
            const std::vector<T> Tv = a.vector<T>("T");                            \
            const std::vector<int> L = a.ints("L");                                \
            const unsigned k = a.uinteger("k");                                    \
            a.done();                                                              \
            const line::trace::MtraceCrossMomentResult<T> r =                      \
                line::trace::fn<T>(Tv, L, k);                                      \
            Json out;                                                              \
            out["MC"] = encode_matrix_defined(r.mc, r.count);                      \
            return out;                                                            \
        }                                                                          \
    };
 
LINE_MTRACE_CROSS_MOMENT_OP(OpMtraceCrossMoment, mtrace_cross_moment)
LINE_MTRACE_CROSS_MOMENT_OP(OpMtraceMomentSimple, mtrace_moment_simple)
#undef LINE_MTRACE_CROSS_MOMENT_OP
 
struct OpMtraceCount {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> Tv = a.vector<T>("T");
        const std::vector<int> A = a.ints("A");
        const T t = a.number<T>("t");
        a.done();
        Json out;
        out["N"] = encode_long_matrix(line::trace::mtrace_count<T>(Tv, A, t).counts);
        return out;
    }
};
 
struct OpMtraceCov {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> Tv = a.vector<T>("T");
        const std::vector<int> A = a.ints("A");
        a.done();
        const std::vector<std::vector<Matrix<T>>> cov = line::trace::mtrace_cov<T>(Tv, A);
        Json c = Json::array();
        for (const std::vector<Matrix<T>>& row : cov) c.push_back(encode_matrices(row));
        Json out;
        out["COV"] = c;
        return out;
    }
};
 
struct OpMtraceIat2counts {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> Tv = a.vector<T>("T");
        const std::vector<int> A = a.ints("A");
        const T scale = a.number<T>("scale");
        a.unsupported("no_mex", "it selects between MATLAB's MEX and script implementations of "
                                "the same counts; the port has one implementation");
        a.done();
        Json out;
        out["C"] = encode_long_matrix(line::trace::mtrace_iat2counts<T>(Tv, A, scale).counts);
        return out;
    }
};
 
struct OpMtraceJoint {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> Tv = a.vector<T>("T");
        const std::vector<int> A = a.ints("A");
        const std::vector<unsigned> i = nonneg_unsigned(a.ints("i"), "mtrace_joint: i");
        if (i.size() != 2) throw InputError("mtrace_joint: 'i' must hold the two orders [i1, i2]");
        a.done();
        const line::trace::MtraceJointResult<T> r = line::trace::mtrace_joint<T>(Tv, A, i[0], i[1]);
        Json out;
        out["JM"] = encode_vector_defined(r.jm, r.count);
        return out;
    }
};
 
struct OpMtraceMean {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> Tv = a.vector<T>("trace");
        const int ntypes = a.required_integer("ntypes");
        const std::vector<int> type = a.ints("type");
        a.done();
        const line::trace::MtraceMeanResult<T> r =
            line::trace::mtrace_mean<T>(Tv, static_cast<long>(ntypes), type);
        Json out;
        out["m"] = encode_vector_defined(r.mean, r.count);
        return out;
    }
};
 
struct OpMtraceMerge {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> t1 = a.vector<T>("t1");
        const std::vector<T> t2 = a.vector<T>("t2");
        a.done();
        const line::trace::MtraceMergeResult<T> r = line::trace::mtrace_merge<T>(t1, t2);
        Json out;
        out["T"] = encode_vector(r.times);
        out["L"] = encode_ints(r.labels);
        return out;
    }
};
 
/** The reference's first argument T is never read: the statistic is of the labels. */
void refuse_unused_times(Args& a, const char* fn) {
    if (a.has("T"))
        throw InputError(std::string(fn) + ": 'T' is not read by the reference, whose result "
                                           "depends on the labels L alone; omit it");
    a.opt("T");
}
 
struct OpMtracePc {
    template <class T>
    static Json run(Args& a) {
        const std::vector<int> C = a.ints("C");
        a.done();
        Json out;
        out["pc"] = encode_vector(line::trace::mtrace_pc<T>(C));
        return out;
    }
};
 
struct OpMtraceSigma {
    template <class T>
    static Json run(Args& a) {
        refuse_unused_times(a, "mtrace_sigma");
        const std::vector<int> L = a.ints("L");
        a.done();
        Json out;
        out["sigma"] = encode_matrix(line::trace::mtrace_sigma<T>(L));
        return out;
    }
};
 
/** sigma(i,j,h) as a C x C x C nested array, from the port's (i, j*C+h) layout. */
struct OpMtraceSigma2 {
    template <class T>
    static Json run(Args& a) {
        refuse_unused_times(a, "mtrace_sigma2");
        const std::vector<int> L = a.ints("L");
        a.done();
        const Matrix<T> s = line::trace::mtrace_sigma2<T>(L);
        const std::size_t C = s.rows();
        Json cube = Json::array();
        for (std::size_t i = 0; i < C; ++i) {
            Json plane = Json::array();
            for (std::size_t j = 0; j < C; ++j) {
                Json row = Json::array();
                for (std::size_t h = 0; h < C; ++h) row.push_back(encode_scalar(s(i, j * C + h)));
                plane.push_back(row);
            }
            cube.push_back(plane);
        }
        Json out;
        out["sigma"] = cube;
        return out;
    }
};
 
struct OpMtraceSplit {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> Tv = a.vector<T>("T");
        const std::vector<int> L = a.ints("L");
        a.done();
        Json out;
        out["TL"] = encode_vectors(line::trace::mtrace_split<T>(Tv, L).traces);
        return out;
    }
};
 
struct OpMtraceSummary {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> Tv = a.vector<T>("T");
        const std::vector<int> C = a.ints("C");
        a.done();
        const line::trace::MtraceSummary<T> r = line::trace::mtrace_summary<T>(Tv, C);
        Json s;
        s["M"] = encode_vector(r.M);
        s["ACF"] = encode_vector(r.acf);
        s["F1"] = encode_matrix(r.F1);
        s["F2"] = encode_matrix(r.F2);
        s["B1"] = encode_matrix(r.B1);
        s["B2"] = encode_matrix(r.B2);
        s["C1"] = encode_matrix(r.C1);
        s["C2"] = encode_matrix(r.C2);
        s["Pc"] = encode_vector(r.Pc);
        s["Pab"] = encode_matrix(r.Pab);
        Json out;
        out["summary"] = s;
        return out;
    }
};
 
// ---------------------------------------------------------------------------
// infer: queue lengths and demands from traces
// ---------------------------------------------------------------------------
 
struct OpInferComputeQlAtArrival {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> at = a.vector<T>("at");
        const std::vector<long> at_jobid = to_longs(a.ints("at_jobid"));
        const std::vector<T> rt = a.vector<T>("rt");
        const std::vector<long> rt_jobid = to_longs(a.ints("rt_jobid"));
        const std::vector<std::size_t> cls =
            one_based_to_zero(a.ints("class"), "infer_compute_ql_at_arrival: class");
        const std::size_t R = a.count("R");
        a.done();
        Json out;
        out["ql"] = encode_matrix(
            line::infer::infer_compute_ql_at_arrival<T>(at, at_jobid, rt, rt_jobid, cls, R));
        return out;
    }
};
 
struct OpInferGetQlenArrival {
    template <class T>
    static Json run(Args& a) {
        std::vector<std::vector<T>> at_ms, rt;
        read_trace_cells<T>(a.get("data"), "infer_get_qlen_arrival: data", at_ms, rt);
        a.done();
        Json out;
        out["ql"] = encode_matrices(line::infer::infer_get_qlen_arrival<T>(at_ms, rt));
        return out;
    }
};
 
struct OpInferLqnFindbyname {
    template <class T>
    static Json run(Args& a) {
        const Json& container = a.get("container");
        const std::string name = string_of(a.get("name"), "infer_lqn_findbyname: name");
        a.done();
        if (!container.is_array())
            throw InputError("infer_lqn_findbyname: 'container' must be a list of elements");
        std::vector<std::string> names;
        for (std::size_t k = 0; k < container.size(); ++k) {
            const Json& e = container[k];
            if (!e.is_object() || e.find("name") == e.end())
                throw InputError("infer_lqn_findbyname: container entry " + std::to_string(k + 1) +
                                 " has no 'name' field");
            names.push_back(string_of(e["name"], "infer_lqn_findbyname: container.name"));
        }
        const std::size_t idx = line::infer::infer_lqn_findbyname(names, name);
        Json out;
        out["el"] = idx == line::infer::INFER_LQN_NOT_FOUND ? Json(nullptr) : container[idx];
        return out;
    }
};
 
struct OpInferLqnGetobs {
    template <class T>
    static Json run(Args& a) {
        const std::vector<std::string> names = strings_of(a.get("names"), "infer_lqn_getobs: names");
        const Json& mj = a.get("metrics");
        const Json& sj = a.get("obsSpec");
        a.done();
        line::infer::LqnMetrics<T> m;
        {
            Args ma(mj, "infer_lqn_getobs: metrics");
            m.QLen = ma.vector_or_empty<T>("QLen");
            m.Util = ma.vector_or_empty<T>("Util");
            m.RespT = ma.vector_or_empty<T>("RespT");
            m.Tput = ma.vector_or_empty<T>("Tput");
            ma.done();
        }
        if (!sj.is_array()) throw InputError("infer_lqn_getobs: 'obsSpec' must be a list of {metric, name}");
        std::vector<line::infer::LqnObsSpec> spec;
        for (const Json& e : sj) {
            Args ea(e, "infer_lqn_getobs: obsSpec");
            const std::string metric = lower(string_of(ea.get("metric"), "infer_lqn_getobs: obsSpec.metric"));
            line::infer::LqnObsSpec s;
            s.name = string_of(ea.get("name"), "infer_lqn_getobs: obsSpec.name");
            ea.done();
            if (metric == "respt")
                s.metric = line::infer::LqnMetric::RespT;
            else if (metric == "util")
                s.metric = line::infer::LqnMetric::Util;
            else if (metric == "tput")
                s.metric = line::infer::LqnMetric::Tput;
            else if (metric == "qlen")
                s.metric = line::infer::LqnMetric::QLen;
            else
                throw InputError("infer_lqn_getobs: Unknown metric '" + metric + "'.");
            spec.push_back(s);
        }
        Json out;
        out["z"] = encode_vector(line::infer::infer_lqn_getobs<T>(names, m, spec));
        return out;
    }
};
 
struct OpInferQmle {
    template <class T>
    static Json run(Args& a) {
        const Matrix<T> Q = a.matrix<T>("Q");
        const std::vector<T> N = a.vector<T>("N");
        const std::vector<T> Z = a.vector<T>("Z");
        a.done();
        Json out;
        out["D"] = encode_matrix(line::infer::infer_qmle<T>(Q, N, Z));
        return out;
    }
};
 
struct OpInferRps {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> rt = a.vector<T>("rt");
        const std::vector<std::size_t> cls = one_based_to_zero(a.ints("class"), "infer_rps: class");
        const Matrix<T> ql = a.matrix<T>("ql");
        const int V = a.required_integer("V");
        a.done();
        Json out;
        out["demandEst"] = encode_vector(line::infer::infer_rps<T>(rt, cls, ql, static_cast<long>(V)));
        return out;
    }
};
 
/**
 * infer_quick_model(open, stations, classes, servers, jobs, routing). The
 * routing of a closed model is one M x M matrix per class (MATLAB's cell
 * form), or a single M x M matrix that every class uses.
 */
struct OpInferQuickModel {
    template <class T>
    static Json run(Args& a) {
        const bool is_open = flag_of(a.get("open"), "infer_quick_model: open");
        std::vector<lang::SchedStrategy> stations;
        for (const std::string& s : strings_of(a.get("stations"), "infer_quick_model: stations"))
            stations.push_back(io::detail::sched_from_json(s));
        const Matrix<T> demands = a.matrix<T>("classes");
        const bool has_servers = a.has("servers"), has_jobs = a.has("jobs"),
                   has_routing = a.has("routing");
        if (has_jobs && !has_servers)
            throw InputError("infer_quick_model: 'jobs' follows 'servers' positionally; give servers too");
        if (has_routing && !(has_servers && has_jobs))
            throw InputError("infer_quick_model: 'routing' follows 'servers' and 'jobs' "
                             "positionally; give both");
        if (has_routing && is_open)
            throw InputError("infer_quick_model: 'routing' applies to a closed model; an open "
                             "model is routed serially from the source to the sink");
        const std::vector<double> servers = a.vector_or_empty<double>("servers");
        const std::vector<double> jobs = a.vector_or_empty<double>("jobs");
        std::vector<Matrix<T>> routing;
        if (has_routing) {
            const Json& rj = a.get("routing");
            const bool cell_form = rj.is_array() && !rj.empty() && rj[0].is_array() &&
                                   !rj[0].empty() && rj[0][0].is_array();
            if (cell_form) {
                routing = a.matrices<T>("routing");
            } else {
                const Matrix<T> P = a.matrix<T>("routing");
                routing.assign(demands.rows(), P);
            }
        }
        a.done();
        qn::Network<T> model =
            has_routing ? line::infer::infer_quick_model<T>(is_open, stations, demands, servers, jobs, routing)
            : has_jobs  ? line::infer::infer_quick_model<T>(is_open, stations, demands, servers, jobs)
            : has_servers ? line::infer::infer_quick_model<T>(is_open, stations, demands, servers)
                          : line::infer::infer_quick_model<T>(is_open, stations, demands);
        Json out;
        out["model"] = io::network_json_envelope(model.get_struct());
        return out;
    }
};
 
/**
 * infer_nhpp_ks(times, T, 'cumRate', f, 'method', m, 'T0', t0). The port's
 * options are positional (cumRate, method, T0); an absent cumRate is the
 * constant rate in both, so method is reached through the empty callable.
 */
struct OpInferNhppKs {
    template <class T>
    static Json run(Args& a) {
        const std::vector<T> times = a.vector<T>("times");
        const T Tend = a.number<T>("T");
        a.unsupported("cumRate", "it is a function handle; without it the rate is taken "
                                 "constant over the interval");
        const bool has_method = a.has("method"), has_t0 = a.has("T0");
        const std::string m = lower(a.text("method", std::string()));
        if (has_t0 && !has_method)
            throw InputError("infer_nhpp_ks: the port takes T0 after method; give 'method' too");
        line::infer::NhppKsResult<T> r;
        if (!has_method) {
            a.done();
            r = line::infer::infer_nhpp_ks<T>(times, Tend);
        } else {
            line::infer::NhppKsMethod method;
            if (m == "cu")
                method = line::infer::NhppKsMethod::Cu;
            else if (m == "lewis")
                method = line::infer::NhppKsMethod::Lewis;
            else
                throw InputError("infer_nhpp_ks: method must be 'cu' or 'lewis'");
            const std::function<T(const T&)> constant;
            if (has_t0) {
                const T t0 = a.number<T>("T0");
                a.done();
                r = line::infer::infer_nhpp_ks<T>(times, Tend, constant, method, t0);
            } else {
                a.opt("T0");
                a.done();
                r = line::infer::infer_nhpp_ks<T>(times, Tend, constant, method);
            }
        }
        Json res;
        res["statistic"] = encode_scalar(r.statistic);
        res["pvalue"] = encode_scalar(r.pvalue);
        res["n"] = encode_count(r.n);
        res["uniforms"] = encode_vector(r.uniforms);
        res["transformed"] = encode_vector(r.transformed);
        Json out;
        out["result"] = res;
        return out;
    }
};
 
// ---- semi-Markov inference (double only) ----------------------------------
 
Json encode_smp_fit(const line::infer::SmpFit& f) {
    Json arcs = Json::array();
    for (const line::infer::SmpArc& c : f.arcs) {
        Json j;
        j["from"] = encode_count(c.from);
        j["to"] = encode_count(c.to);
        j["family"] = c.family;
        j["theta"] = encode_vector(c.theta);
        j["params"] = encode_vector(c.params);
        j["beta"] = encode_vector(c.beta);
        j["se"] = encode_vector(c.se);
        j["prob"] = c.prob;
        j["nevents"] = encode_count(c.nevents);
        arcs.push_back(j);
    }
    Json j;
    j["approach"] = f.approach;
    j["nstates"] = encode_count(f.nstates);
    j["ncov"] = encode_count(f.ncov);
    j["arcs"] = arcs;
    j["loglik"] = f.loglik;
    j["nparams"] = encode_count(f.nparams);
    j["aic"] = f.aic;
    j["P"] = encode_matrix(f.P);
    return j;
}
 
std::vector<double> doubles_nan(const Json& j, const std::string& where) {
    const Matrix<double> m = matrix_or_fill<double>(j, where, std::numeric_limits<double>::quiet_NaN());
    return column_major(m);
}
 
/** The `fit` struct infer_smp returns, read back field for field. */
line::infer::SmpFit read_smp_fit(const Json& j) {
    const std::string w = "infer_smp_hazard: fit";
    Args f(j, w);
    line::infer::SmpFit fit;
    fit.approach = string_of(f.get("approach"), w + ".approach");
    fit.nstates = f.count("nstates");
    fit.ncov = f.count("ncov");
    fit.loglik = number_or_fill<double>(f.get("loglik"), w + ".loglik",
                                        std::numeric_limits<double>::quiet_NaN());
    fit.nparams = f.count("nparams");
    fit.aic = number_or_fill<double>(f.get("aic"), w + ".aic", std::numeric_limits<double>::quiet_NaN());
    fit.P = matrix_or_fill<double>(f.get("P"), w + ".P", std::numeric_limits<double>::quiet_NaN());
    const Json& arcs = f.get("arcs");
    f.done();
    if (!arcs.is_array()) throw InputError(w + ".arcs: expected a list of arcs");
    for (const Json& e : arcs) {
        Args c(e, w + ".arcs");
        line::infer::SmpArc arc;
        arc.from = c.count("from");
        arc.to = c.count("to");
        arc.family = string_of(c.get("family"), w + ".arcs.family");
        arc.theta = doubles_nan(c.get("theta"), w + ".arcs.theta");
        arc.params = doubles_nan(c.get("params"), w + ".arcs.params");
        arc.beta = doubles_nan(c.get("beta"), w + ".arcs.beta");
        arc.se = doubles_nan(c.get("se"), w + ".arcs.se");
        arc.prob = number_or_fill<double>(c.get("prob"), w + ".arcs.prob",
                                          std::numeric_limits<double>::quiet_NaN());
        arc.nevents = c.count("nevents");
        c.done();
        fit.arcs.push_back(arc);
    }
    return fit;
}
 
struct OpInferSmp {
    static Json run(Args& a) {
        const std::string w = "infer_smp";
        const std::vector<std::size_t> from = nonneg_sizes(a.ints("from"), w + ": from");
        const std::vector<std::size_t> to = nonneg_sizes(a.ints("to"), w + ": to");
        const std::vector<double> dur = a.vector<double>("dur");
        line::infer::SmpOptions opt;
        if (a.has("approach")) opt.approach = lower(a.text("approach", std::string()));
        if (const Json* d = a.opt("dist")) {
            if (d->is_string()) {
                opt.dist = lower(d->get<std::string>());
            } else {
                if (!d->is_array()) throw InputError(w + ": 'dist' must be a family name or an E x E table of them");
                for (const Json& row : *d) {
                    if (!row.is_array()) throw InputError(w + ": 'dist' table rows must be lists");
                    std::vector<std::string> r;
                    for (const Json& e : row) r.push_back(e.is_null() ? std::string() : lower(string_of(e, w + ": dist")));
                    opt.dist_table.push_back(r);
                }
            }
        }
        if (a.has("covariates")) opt.covariates = a.matrix<double>("covariates");
        if (const Json* cm = a.opt("covmask")) {
            if (!cm->is_array()) throw InputError(w + ": 'covmask' must be an E x E x p array");
            for (const Json& plane : *cm) {
                if (!plane.is_array()) throw InputError(w + ": 'covmask' must be an E x E x p array");
                std::vector<std::vector<bool>> p;
                for (const Json& row : plane) {
                    if (!row.is_array()) throw InputError(w + ": 'covmask' must be an E x E x p array");
                    std::vector<bool> r;
                    for (const Json& e : row) r.push_back(flag_of(e, w + ": covmask"));
                    p.push_back(r);
                }
                opt.covmask.push_back(p);
            }
        }
        if (a.has("nstates")) opt.nstates = a.count("nstates");
        if (a.has("tolx")) opt.tolx = a.number<double>("tolx");
        if (a.has("tolf")) opt.tolf = a.number<double>("tolf");
        if (a.has("maxiter")) opt.maxiter = a.count("maxiter");
        a.opt("approach");
        a.opt("covariates");
        a.opt("nstates");
        a.opt("tolx");
        a.opt("tolf");
        a.opt("maxiter");
        a.done();
        Json out;
        out["fit"] = encode_smp_fit(line::infer::infer_smp(from, to, dur, opt));
        return out;
    }
};
 
struct OpInferSmpHazard {
    static Json run(Args& a) {
        const line::infer::SmpFit fit = read_smp_fit(a.get("fit"));
        const std::vector<double> t = a.vector<double>("t");
        const bool has_z = a.has("z");
        const std::vector<double> z = a.vector_or_empty<double>("z");
        a.done();
        const std::vector<line::infer::SmpHazard> h =
            has_z ? line::infer::infer_smp_hazard(fit, t, z) : line::infer::infer_smp_hazard(fit, t);
        Json hz = Json::array();
        for (const line::infer::SmpHazard& e : h) {
            Json j;
            j["from"] = encode_count(e.from);
            j["to"] = encode_count(e.to);
            j["intensity"] = encode_vector(e.intensity);
            j["sojourn"] = encode_vector(e.sojourn);
            j["cumulative"] = encode_vector(e.cumulative);
            j["cif"] = encode_vector(e.cif);
            hz.push_back(j);
        }
        Json out;
        out["hz"] = hz;
        return out;
    }
};
 
struct OpInferSmpNelsonAalen {
    static Json run(Args& a) {
        const std::string w = "infer_smp_nelson_aalen";
        const std::vector<std::size_t> from = nonneg_sizes(a.ints("from"), w + ": from");
        const std::vector<std::size_t> to = nonneg_sizes(a.ints("to"), w + ": to");
        const std::vector<double> dur = a.vector<double>("dur");
        const bool has_e = a.has("nstates");
        const std::size_t E = has_e ? a.count("nstates") : 0;
        a.opt("nstates");
        a.done();
        const std::vector<line::infer::SmpNelsonAalen> r =
            has_e ? line::infer::infer_smp_nelson_aalen(from, to, dur, E)
                  : line::infer::infer_smp_nelson_aalen(from, to, dur);
        Json na = Json::array();
        for (const line::infer::SmpNelsonAalen& e : r) {
            Json j;
            j["from"] = encode_count(e.from);
            j["to"] = encode_count(e.to);
            j["time"] = encode_vector(e.time);
            j["A"] = encode_vector(e.A);
            j["var"] = encode_vector(e.var);
            na.push_back(j);
        }
        Json out;
        out["na"] = na;
        return out;
    }
};
 
/**
 * infer_smp_family(family, mode, ...), one entry for the five modes, with the
 * mode's positional arguments under the names of the reference's help text.
 */
struct OpInferSmpFamily {
    static Json run(Args& a) {
        const std::string w = "infer_smp_family";
        const std::string family = string_of(a.get("family"), w + ": family");
        const std::string mode = string_of(a.get("mode"), w + ": mode");
        // Validates the family name before any parameter vector is indexed.
        const std::size_t npar = line::infer::infer_smp_family_npar(family);
        auto theta = [&]() {
            const std::vector<double> th = a.vector<double>("THETA");
            if (th.size() != npar)
                throw InputError(w + ": THETA must have " + std::to_string(npar) + " entries for '" +
                                 family + "'");
            return th;
        };
        Json out;
        if (mode == "npar") {
            a.done();
            out["N"] = encode_count(npar);
        } else if (mode == "logfs") {
            const std::vector<double> th = theta();
            const std::vector<double> t = a.vector<double>("T");
            const std::vector<double> eta = a.vector_or_empty<double>("ETA");
            if (eta.size() > 1 && eta.size() != t.size())
                throw InputError(w + ": ETA must be a scalar or have one entry per instant of T");
            a.done();
            std::vector<double> logf, logS;
            line::infer::infer_smp_family_logfs(family, th, t, eta, logf, logS);
            out["LOGF"] = encode_vector(logf);
            out["LOGS"] = encode_vector(logS);
        } else if (mode == "natural") {
            const std::vector<double> th = theta();
            a.done();
            out["P"] = encode_vector(line::infer::infer_smp_family_natural(family, th));
        } else if (mode == "theta0") {
            const std::vector<double> tev = a.vector<double>("TEV");
            const double texp = a.number<double>("TEXP");
            a.done();
            out["THETA"] = encode_vector(line::infer::infer_smp_family_theta0(family, tev, texp));
        } else if (mode == "dist") {
            const std::vector<double> th = theta();
            const double eta = a.number<double>("ETA");
            a.done();
            out["D"] = io::detail::dist_to_json(line::infer::infer_smp_family_dist(family, th, eta));
        } else {
            throw InputError(w + ": Unknown mode '" + mode +
                             "': use 'npar', 'logfs', 'natural', 'theta0' or 'dist'.");
        }
        return out;
    }
};
 
// ---- variational inference --------------------------------------------------
 
/**
 * infer_variational(spec, options). The spec is MATLAB's struct; NaN entries
 * (arcrate where the arc is governed by a parameter, obsData where nothing was
 * observed) are written as null, and so is an Inf capacity. obsRange and
 * capacity are flattened column-major, as the reference's reshape does.
 */
struct OpInferVariational {
    template <class T>
    static Json run(Args& a) {
        const std::string w = "infer_variational";
        const Json& sj = a.get("spec");
        const Json* oj = a.opt("options");
        a.done();
        const T nan = num_traits<T>::from_double(std::numeric_limits<double>::quiet_NaN());
        const T inf = num_traits<T>::from_double(std::numeric_limits<double>::infinity());
 
        line::infer::VariationalSpec<T> spec;
        {
            Args s(sj, w + ": spec");
            const Matrix<T> arcs = s.matrix<T>("arcs");
            if (!arcs.empty() && arcs.cols() != 3)
                throw InputError(w + ": spec.arcs must have three columns [i j c].");
            const std::vector<int> arcint = int_vector_from_json(s.get("arcs"), w + ": spec.arcs");
            const std::vector<std::size_t> arcsz = nonneg_sizes(arcint, w + ": spec.arcs");
            for (std::size_t e = 0; e < arcs.rows(); ++e)
                spec.arcs.push_back(std::array<std::size_t, 3>{arcsz[3 * e], arcsz[3 * e + 1], arcsz[3 * e + 2]});
            spec.x0 = s.matrix<T>("x0");
            spec.sched = s.ints("sched");
            spec.nservers = s.vector<T>("nservers");
            spec.routeprob = s.vector<T>("routeprob");
            spec.arcparam = nonneg_sizes(s.ints("arcparam"), w + ": spec.arcparam");
            spec.arcrate = column_major(matrix_or_fill<T>(s.get("arcrate"), w + ": spec.arcrate", nan));
            spec.alpha0 = s.vector<T>("alpha0");
            spec.beta0 = s.vector<T>("beta0");
            spec.obsTimes = s.vector<T>("obsTimes");
            spec.obsData = matrix_or_fill<T>(s.get("obsData"), w + ": spec.obsData",
                                             line::infer::VariationalSpec<T>::unobserved());
            spec.obsRange = column_major(s.matrix<T>("obsRange"));
            spec.epsilon = s.number<T>("epsilon");
            if (const Json* c = s.opt("capacity"))
                spec.capacity = column_major(matrix_or_fill<T>(*c, w + ": spec.capacity", inf));
            s.done();
        }
        line::infer::VariationalOptions<T> opt;
        if (oj) {
            Args o(*oj, w + ": options");
            if (o.has("verbose")) opt.verbose = o.required_integer("verbose");
            if (o.has("iter_max")) opt.iter_max = o.count("iter_max");
            if (o.has("tol")) opt.tol = o.number<double>("tol");
            if (o.has("nsamples")) opt.nsamples = o.count("nsamples");
            if (o.has("delta")) opt.delta = o.number<double>("delta");
            if (o.has("floor")) opt.floor = o.number<double>("floor");
            if (o.has("rate_max")) opt.rate_max = o.number<double>("rate_max");
            if (o.has("rate_cap_factor")) opt.rate_cap_factor = o.number<double>("rate_cap_factor");
            if (o.has("unifmax")) opt.unifmax = o.number<double>("unifmax");
            if (o.has("unif_tol")) opt.unif_tol = o.number<double>("unif_tol");
            if (o.has("unif_max_terms")) opt.unif_max_terms = o.count("unif_max_terms");
            if (o.has("tmax")) opt.tmax = o.number<double>("tmax");
            if (o.has("dt")) opt.dt = o.number<double>("dt");
            if (o.has("ngrid")) opt.ngrid = static_cast<long>(o.count("ngrid"));
            if (o.has("ymax")) opt.ymax = static_cast<long>(o.count("ymax"));
            for (const char* k : {"verbose", "iter_max", "tol", "nsamples", "delta", "floor",
                                  "rate_max", "rate_cap_factor", "unifmax", "unif_tol",
                                  "unif_max_terms", "tmax", "dt", "ngrid", "ymax"})
                o.opt(k);
            o.done();
        }
        const line::infer::VariationalResult<T> r =
            oj ? line::infer::infer_variational<T>(spec, opt) : line::infer::infer_variational<T>(spec);
        Json res;
        res["alpha"] = encode_vector(r.alpha);
        res["beta"] = encode_vector(r.beta);
        res["rates"] = encode_vector(r.rates);
        res["meanServiceTime"] = encode_vector(r.mean_service_time);
        res["bound"] = encode_vector(r.bound);
        res["alphaTrace"] = encode_matrix(r.alpha_trace);
        res["betaTrace"] = encode_matrix(r.beta_trace);
        res["Y"] = encode_matrices(r.Y);
        res["nu"] = encode_matrices(r.nu);
        res["tgrid"] = encode_vector(r.tgrid);
        res["qlen"] = encode_matrix(r.qlen);
        res["iter"] = encode_count(r.iter);
        res["converged"] = r.converged;
        res["tailmass"] = encode_scalar(r.tailmass);
        Json out;
        out["out"] = res;
        return out;
    }
};
 
// ---- MLPS / MINPS (double only) --------------------------------------------
 
const char* const kMlpsModelWhy =
    "the reference reads the delay rates and the server count off the model; pass them as "
    "'muZ' (per-class think rates) and 'nCores' instead";
 
std::vector<line::api::MlpsSample> read_mlps_samples(Args& a, const std::string& w) {
    const std::vector<double> rt = a.vector<double>("rt");
    const std::vector<std::size_t> cls = one_based_sizes(a.ints("class"), w + ": class");
    const Matrix<double> ql = a.matrix<double>("ql");
    if (cls.size() != rt.size() || ql.rows() != rt.size())
        throw InputError(w + ": rt, class and ql must have one row per sample");
    std::vector<line::api::MlpsSample> samples(rt.size());
    for (std::size_t i = 0; i < rt.size(); ++i) {
        samples[i].rt = rt[i];
        samples[i].cls = cls[i];
        samples[i].ql.resize(ql.cols());
        for (std::size_t c = 0; c < ql.cols(); ++c) samples[i].ql[c] = ql(i, c);
    }
    return samples;
}
 
struct OpInferMlps {
    static Json run(Args& a) {
        a.unsupported("model", kMlpsModelWhy);
        a.unsupported("node", kMlpsModelWhy);
        const std::vector<double> muZ = a.vector<double>("muZ");
        const double nCores = a.number<double>("nCores");
        const std::vector<line::api::MlpsSample> samples = read_mlps_samples(a, "infer_mlps");
        a.done();
        Json out;
        out["demandEst"] = encode_vector(line::api::infer_mlps(muZ, nCores, samples));
        return out;
    }
};
 
/**
 * infer_minps_setup(data, initSample, sampleSize, V, model, node): prepare the
 * window, then run MINPS on it, as the reference's last line does. The port
 * also reports what the preparation decided (lambda, Wexp, classMap).
 */
struct OpInferMinpsSetup {
    static Json run(Args& a) {
        const std::string w = "infer_minps_setup";
        a.unsupported("model", kMlpsModelWhy);
        a.unsupported("node", kMlpsModelWhy);
        if (a.has("V"))
            throw InputError(w + ": 'V' is not read by the reference, which takes the server count "
                                 "from the node; pass it as 'nCores'");
        a.opt("V");
        std::vector<std::vector<double>> at_ms, rt;
        read_trace_cells<double>(a.get("data"), w + ": data", at_ms, rt);
        std::vector<line::api::MinpsClassTrace> traces(at_ms.size());
        for (std::size_t k = 0; k < at_ms.size(); ++k) {
            traces[k].arrival_ms = at_ms[k];
            traces[k].rt = rt[k];
        }
        const std::size_t initSample = a.count("initSample");
        const std::size_t sampleSize = a.count("sampleSize");
        const std::vector<double> muZ = a.vector<double>("muZ");
        const double nCores = a.number<double>("nCores");
        a.done();
        const line::api::MinpsSetup st = line::api::infer_minps_setup(traces, initSample, sampleSize);
        if (muZ.size() != st.classMap.size())
            throw InputError(w + ": 'muZ' must have one entry per class with samples (" +
                             std::to_string(st.classMap.size()) + "), in the order of classMap");
        Json out;
        out["demandEst"] = encode_vector(line::api::infer_minps(muZ, nCores, st.samples));
        out["lambda"] = encode_vector(st.lambda);
        out["Wexp"] = st.threads;
        out["classMap"] = encode_sizes(st.classMap);
        return out;
    }
};
 
 
/**
 * demand = infer_gibbs(data, nbCores, tol, seed). data is the 6 x (K+1) cell
 * array, of which rows 3 (arrival times, ms), 4 (response times, s) and 6
 * (the think-time observations) of the first K columns are read; tol defaults
 * to 1e-3 as in the reference. The test-set draw and the slice sampler take a
 * pfqn::McRng named by the required seed: MATLAB uses the session RNG, so its
 * estimates agree with these in law, not draw for draw.
 */
struct OpInferGibbs {
    static Json run(Args& a) {
        const std::string w = "infer_gibbs";
        const Json& data = a.get("data");
        std::vector<std::vector<double>> at_ms, rt;
        read_trace_cells<double>(data, w + ": data", at_ms, rt);
        if (data.size() < 6)
            throw InputError(w + ": data needs row 6 (the think-time observations) as well as "
                                 "rows 3 and 4");
        std::vector<line::infer::GibbsTrace<double>> traces(at_ms.size());
        for (std::size_t k = 0; k < at_ms.size(); ++k) {
            traces[k].arrival_ms = at_ms[k];
            traces[k].respt_s = rt[k];
            traces[k].think_obs =
                vector_from_json<double>(data[5][k], w + ": data{6," + std::to_string(k + 1) + "}");
        }
        const double nbCores = a.number<double>("nbCores");
        const double tol = a.scalar<double>("tol", 1e-3);
        pfqn::McRng rng(required_seed(a, w));
        a.done();
        Json out;
        out["demand"] = encode_vector(line::infer::infer_gibbs<double>(traces, nbCores, tol, rng));
        return out;
    }
};
 
}  // namespace
 
void register_api_trace(DispatchTable& t) {
    // trace
    t["trace_mean"] = make<OpTraceMean>();
    t["trace_var"] = make<OpTraceVar>();
    t["trace_scv"] = make<OpTraceScv>();
    t["trace_skew"] = make<NeedsTranscendental<OpTraceSkew>>();
    t["trace_idc"] = make<OpTraceIdc>();
    t["trace_acf"] = make<OpTraceAcf>();
    t["trace_bicov"] = make<OpTraceBicov>();
    t["trace_gamma"] = make<OpTraceGamma>();
    t["trace_iat2bins"] = make<OpTraceIat2bins>();
    t["trace_iat2counts"] = make<OpTraceIat2counts>();
    t["trace_idi"] = make<OpTraceIdi>();
    t["trace_joint"] = make<OpTraceJoint>();
    t["trace_pmf"] = make<OpTracePmf>();
    t["trace_summary"] = make<NeedsTranscendental<OpTraceSummary>>();
    // mtrace
    t["mtrace_moment"] = make<OpMtraceMoment>();
    t["mtrace_forward_moment"] = make<OpMtraceForwardMoment>();
    t["mtrace_backward_moment"] = make<OpMtraceBackwardMoment>();
    t["mtrace_cross_moment"] = make<OpMtraceCrossMoment>();
    t["mtrace_moment_simple"] = make<OpMtraceMomentSimple>();
    t["mtrace_count"] = make<OpMtraceCount>();
    t["mtrace_cov"] = make<OpMtraceCov>();
    t["mtrace_iat2counts"] = make<OpMtraceIat2counts>();
    t["mtrace_joint"] = make<OpMtraceJoint>();
    t["mtrace_mean"] = make<OpMtraceMean>();
    t["mtrace_merge"] = make<OpMtraceMerge>();
    t["mtrace_pc"] = make<OpMtracePc>();
    t["mtrace_sigma"] = make<OpMtraceSigma>();
    t["mtrace_sigma2"] = make<OpMtraceSigma2>();
    t["mtrace_split"] = make<OpMtraceSplit>();
    t["mtrace_summary"] = make<OpMtraceSummary>();
    // infer
    t["infer_compute_ql_at_arrival"] = make<OpInferComputeQlAtArrival>();
    t["infer_get_qlen_arrival"] = make<OpInferGetQlenArrival>();
    t["infer_lqn_findbyname"] = make<OpInferLqnFindbyname>();
    t["infer_lqn_getobs"] = make<OpInferLqnGetobs>();
    t["infer_qmle"] = make<OpInferQmle>();
    t["infer_rps"] = make<OpInferRps>();
    t["infer_quick_model"] = make<OpInferQuickModel>();
    t["infer_nhpp_ks"] = make<NeedsTranscendental<OpInferNhppKs>>();
    t["infer_variational"] = make<NeedsTranscendental<OpInferVariational>>();
    t["infer_smp"] = make<DoubleOnly<OpInferSmp>>();
    t["infer_smp_hazard"] = make<DoubleOnly<OpInferSmpHazard>>();
    t["infer_smp_nelson_aalen"] = make<DoubleOnly<OpInferSmpNelsonAalen>>();
    t["infer_smp_family"] = make<DoubleOnly<OpInferSmpFamily>>();
    t["infer_mlps"] = make<DoubleOnly<OpInferMlps>>();
    t["infer_minps_setup"] = make<DoubleOnly<OpInferMinpsSetup>>();
    t["infer_gibbs"] = make<DoubleOnly<OpInferGibbs>>();
}
 
void register_library_only_trace(LibraryOnlyTable& t) {
    t["infer_lqn_jacobian"] =
        "its first argument hfun is a function handle (the LQN observation map), which JSON "
        "cannot carry.";
    t["infer_lqn_ekf"] =
        "its first argument hfun is a function handle (the LQN observation map), which JSON "
        "cannot carry.";
    t["infer_lqn"] =
        "it calibrates a LayeredNetwork model in place through a solver callable; no line-cli "
        "analysis serves LQN inference yet.";
    t["infer_lqn_setparams"] =
        "it writes parameters into a LayeredNetwork model object in place; there is no data "
        "result to return.";
    t["infer_fmlps"] =
        "it reads the Network model and the station node (NetworkStruct), not a data argument; "
        "no line-cli analysis serves demand inference yet.";
    t["infer_smp_env"] =
        "it takes the per-stage Network models and returns an Environment model object; solve "
        "the resulting environment with line-cli -s env.";
}
 
}  // namespace reg
}  // namespace line