lqn_export_nlp.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_LQN_EXPORT_NLP_H
#define LINE_API_LQN_EXPORT_NLP_H
 
/**
 * @file
 * @ingroup api_lqn
 * A layered queueing network as ONE nonlinear program: linear constraints over
 * the flows, the service-time composition and the closures, plus a nonlinear
 * objective that is the sum of squares of the QD-AMVA congestion and Little laws
 * (Casale, Perez and Wang, PERFORMANCE 2015).
 *
 * Twin of python line_solver.api.lqn.lqn_export_nlp. The python exporter EMITS
 * a self-contained script; that emission is python-only by design. What is
 * ported is the PROGRAM the script builds, constructed in memory: the same
 * element numbering, the same variable layout
 *
 *     Xe(NE) Xa(NA) Xc(NC) W(NA) S(NE) R(NC) Q(NA) Qt(NC) BC(NC) BR(NR) GH(NGH) GT(NGT) SA(NE if forking)
 *
 * the same row order of AEQ/AUB, the same residual order, the same refusals
 * (`LINE:unsupported:NlpExport`) and the same solve: Newton on the square
 * system with a minimum-norm least-squares step, escalating the warm-start
 * sweeps 0, 5, 20, 100 before falling back to a constrained descent.
 *
 * THE FALLBACK DIFFERS FROM PYTHON. Python falls back to scipy's trust-constr;
 * this port runs a projected Levenberg-Marquardt on the merit
 * ||r||^2 + ||AEQ x - BEQ||^2 + ||max(0, AUB x - BUB)||^2 over x >= 0, identical
 * in MATLAB, the JAR and C++. Every model the reference suite solves is rooted
 * by the Newton ladder, which reaches the same root in every codebase.
 *
 * Rate functions r_i(n) = softmin(n, m_i) [* lld_i(n)] and the AND-fork closure
 * are differentiated by COMPLEX STEP, as in python. A declared lld multiplier is
 * an expression in `n` over a small grammar shared by the four codebases:
 *
 *     expr  := term (('+'|'-') term)*        term  := unary (('*'|'/') unary)*
 *     unary := ('+'|'-') unary | power       power := atom ('**' unary)?
 *     atom  := NUMBER | 'n' | FUNC '(' expr {',' expr} ')' | '(' expr ')'
 *     FUNC  := exp | log | sqrt | sin | cos | tanh | abs | min | max, optionally np. or math.
 *
 * `abs` is the real modulus and `min`/`max` reject a complex argument, so both
 * are caught by the analyticity check exactly as python catches them.
 */
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstdio>
#include <cstdlib>
#include <limits>
#include <map>
#include <memory>
#include <set>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
#include "line/api/lqn/lqn_balance_equations.h"
#include "line/lang/lqn/lqn_struct.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace lqn {
 
/**
 * A declared queue-dependent rate multiplier of one station: an expression in
 * `n`, or a tabulated lldscaling row, which is refused because its interpolant
 * has a corner at every knot.
 */
struct NlpLld {
    std::string expr;
    std::vector<double> table;
    bool tabulated = false;
    NlpLld() {}
    NlpLld(const std::string& e) : expr(e) {}  // NOLINT(runtime/explicit): a string IS an lld
    NlpLld(const char* e) : expr(e) {}         // NOLINT(runtime/explicit)
    static NlpLld tabulated_row(const std::vector<double>& v) {
        NlpLld l;
        l.table = v;
        l.tabulated = true;
        return l;
    }
};
 
/** Options of the program: python export_nlp(lld=, alpha=, fork=). */
struct NlpOptions {
    std::map<std::string, NlpLld> lld;  ///< station name -> multiplier
    double alpha = 20.0;                ///< softmin sharpness
    std::string fork = "refuse";        ///< 'refuse' or 'exp'
};
 
namespace nlp_detail {
 
using cplx = std::complex<double>;
 
inline std::string fmt(const char* f, double v) {
    char b[64];
    std::snprintf(b, sizeof(b), f, v);
    return b;
}
 
/** Python's repr of a float: the shortest round-tripping digits, in python's layout. */
inline std::string py_repr(double v) {
    if (std::isnan(v)) return "nan";
    if (std::isinf(v)) return v > 0 ? "inf" : "-inf";
    if (v == 0.0) return std::signbit(v) ? "-0.0" : "0.0";
    char b[64];
    int p = 1;
    for (; p <= 17; ++p) {
        std::snprintf(b, sizeof(b), "%.*e", p - 1, v);
        if (std::strtod(b, nullptr) == v) break;
    }
    const std::string s(b);
    const std::size_t epos = s.find('e');
    const int ex = std::atoi(s.c_str() + epos + 1);
    std::string mant = s.substr(0, epos);
    const bool neg = mant[0] == '-';
    if (neg) mant = mant.substr(1);
    std::string digits;
    for (char c : mant)
        if (c != '.') digits += c;
    while (digits.size() > 1 && digits.back() == '0') digits.pop_back();
    std::string out;
    if (ex >= -4 && ex < 16) {
        if (ex >= 0) {
            if (int(digits.size()) <= ex + 1) {
                out = digits + std::string(std::size_t(ex + 1 - int(digits.size())), '0') + ".0";
            } else {
                out = digits.substr(0, std::size_t(ex + 1)) + "." + digits.substr(std::size_t(ex + 1));
            }
        } else {
            out = "0." + std::string(std::size_t(-ex - 1), '0') + digits;
        }
    } else {
        out = digits.substr(0, 1);
        if (digits.size() > 1) out += "." + digits.substr(1);
        char e[16];
        std::snprintf(e, sizeof(e), "e%c%02d", ex < 0 ? '-' : '+', ex < 0 ? -ex : ex);
        out += e;
    }
    return neg ? "-" + out : out;
}
 
/** Python's repr of a value that may be complex. */
inline std::string py_repr(const cplx& z) {
    if (z.imag() == 0.0) return py_repr(z.real());
    std::string im = py_repr(z.imag());
    if (im.size() > 2 && im.compare(im.size() - 2, 2, ".0") == 0) im = im.substr(0, im.size() - 2);
    std::string re = py_repr(z.real());
    if (re.size() > 2 && re.compare(re.size() - 2, 2, ".0") == 0) re = re.substr(0, re.size() - 2);
    return "(" + re + (z.imag() < 0 || std::signbit(z.imag()) ? "" : "+") + im + "j)";
}
 
/** An evaluation failure, carrying the python exception name the message reports. */
struct LldEvalError : std::runtime_error {
    std::string kind;
    LldEvalError(const std::string& k, const std::string& m) : std::runtime_error(m), kind(k) {}
};
 
/** A parse failure: the text is not an expression of the grammar. */
struct LldSyntaxError : std::runtime_error {
    explicit LldSyntaxError(const std::string& m) : std::runtime_error(m) {}
};
 
/**
 * An lld expression, parsed once and evaluated over std::complex<double>.
 *
 * Each operation takes REAL semantics when its operands are real, as python
 * does on a float (log and sqrt of a negative number are nan, math.* raises),
 * and complex semantics otherwise, which is what the complex step reads.
 */
class LldExpr {
public:
    explicit LldExpr(const std::string& src) : src_(src) {
        pos_ = 0;
        next();
        root_ = parse_expr();
        if (tok_.kind != Tok::END) throw LldSyntaxError("invalid syntax");
    }
 
    cplx eval(const cplx& n) const { return eval_node(root_, n); }
 
private:
    struct Tok {
        enum Kind { END, NUM, ID, OP, LP, RP, COMMA } kind = END;
        std::string text;
        double num = 0.0;
    };
    struct Node {
        enum Kind { NUM, VAR, ID, NEG, POS, ADD, SUB, MUL, DIV, POW, CALL } kind = NUM;
        double num = 0.0;
        std::string name;
        std::vector<int> kids;
    };
 
    std::string src_;
    std::size_t pos_ = 0;
    Tok tok_;
    std::vector<Node> nodes_;
    int root_ = -1;
 
    static bool idstart(char c) { return std::isalpha(static_cast<unsigned char>(c)) || c == '_'; }
    static bool idchar(char c) { return std::isalnum(static_cast<unsigned char>(c)) || c == '_'; }
    static bool digit(char c) { return c >= '0' && c <= '9'; }
 
    void next() {
        while (pos_ < src_.size() && std::isspace(static_cast<unsigned char>(src_[pos_]))) ++pos_;
        tok_ = Tok();
        if (pos_ >= src_.size()) return;
        const char c = src_[pos_];
        if (digit(c) || (c == '.' && pos_ + 1 < src_.size() && digit(src_[pos_ + 1]))) {
            const std::size_t b = pos_;
            while (pos_ < src_.size() && digit(src_[pos_])) ++pos_;
            if (pos_ < src_.size() && src_[pos_] == '.') {
                ++pos_;
                while (pos_ < src_.size() && digit(src_[pos_])) ++pos_;
            }
            if (pos_ < src_.size() && (src_[pos_] == 'e' || src_[pos_] == 'E')) {
                std::size_t q = pos_ + 1;
                if (q < src_.size() && (src_[q] == '+' || src_[q] == '-')) ++q;
                if (q >= src_.size() || !digit(src_[q])) throw LldSyntaxError("invalid decimal literal");
                pos_ = q;
                while (pos_ < src_.size() && digit(src_[pos_])) ++pos_;
            }
            if (pos_ < src_.size() && idchar(src_[pos_])) throw LldSyntaxError("invalid decimal literal");
            tok_.kind = Tok::NUM;
            tok_.text = src_.substr(b, pos_ - b);
            tok_.num = std::strtod(tok_.text.c_str(), nullptr);
            return;
        }
        if (idstart(c)) {
            const std::size_t b = pos_;
            while (pos_ < src_.size() && idchar(src_[pos_])) ++pos_;
            while (pos_ + 1 < src_.size() && src_[pos_] == '.' && idstart(src_[pos_ + 1])) {
                ++pos_;
                while (pos_ < src_.size() && idchar(src_[pos_])) ++pos_;
            }
            tok_.kind = Tok::ID;
            tok_.text = src_.substr(b, pos_ - b);
            return;
        }
        ++pos_;
        if (c == '*' && pos_ < src_.size() && src_[pos_] == '*') {
            ++pos_;
            tok_.kind = Tok::OP;
            tok_.text = "**";
            return;
        }
        if (c == '+' || c == '-' || c == '*' || c == '/') {
            tok_.kind = Tok::OP;
            tok_.text = std::string(1, c);
            return;
        }
        if (c == '(') { tok_.kind = Tok::LP; return; }
        if (c == ')') { tok_.kind = Tok::RP; return; }
        if (c == ',') { tok_.kind = Tok::COMMA; return; }
        throw LldSyntaxError("invalid syntax");
    }
 
    int add(Node::Kind k, std::vector<int> kids, double num = 0.0, const std::string& name = "") {
        Node nd;
        nd.kind = k;
        nd.kids = std::move(kids);
        nd.num = num;
        nd.name = name;
        nodes_.push_back(nd);
        return int(nodes_.size()) - 1;
    }
 
    bool is_op(const char* s) const { return tok_.kind == Tok::OP && tok_.text == s; }
 
    int parse_expr() {
        int l = parse_term();
        while (is_op("+") || is_op("-")) {
            const bool plus = is_op("+");
            next();
            const int r = parse_term();
            l = add(plus ? Node::ADD : Node::SUB, {l, r});
        }
        return l;
    }
    int parse_term() {
        int l = parse_unary();
        while (is_op("*") || is_op("/")) {
            const bool mul = is_op("*");
            next();
            const int r = parse_unary();
            l = add(mul ? Node::MUL : Node::DIV, {l, r});
        }
        return l;
    }
    int parse_unary() {
        if (is_op("+") || is_op("-")) {
            const bool neg = is_op("-");
            next();
            const int u = parse_unary();
            return add(neg ? Node::NEG : Node::POS, {u});
        }
        return parse_power();
    }
    int parse_power() {
        const int a = parse_atom();
        if (is_op("**")) {
            next();
            const int e = parse_unary();
            return add(Node::POW, {a, e});
        }
        return a;
    }
    int parse_atom() {
        if (tok_.kind == Tok::NUM) {
            const double v = tok_.num;
            next();
            return add(Node::NUM, {}, v);
        }
        if (tok_.kind == Tok::ID) {
            const std::string name = tok_.text;
            next();
            if (tok_.kind == Tok::LP) {
                next();
                std::vector<int> args;
                if (tok_.kind != Tok::RP) {
                    args.push_back(parse_expr());
                    while (tok_.kind == Tok::COMMA) {
                        next();
                        args.push_back(parse_expr());
                    }
                }
                if (tok_.kind != Tok::RP) throw LldSyntaxError("'(' was never closed");
                next();
                return add(Node::CALL, args, 0.0, name);
            }
            if (name == "n") return add(Node::VAR, {});
            if (name == "np.pi" || name == "math.pi") return add(Node::NUM, {}, 3.141592653589793);
            if (name == "np.e" || name == "math.e") return add(Node::NUM, {}, 2.718281828459045);
            return add(Node::ID, {}, 0.0, name);
        }
        if (tok_.kind == Tok::LP) {
            next();
            const int e = parse_expr();
            if (tok_.kind != Tok::RP) throw LldSyntaxError("'(' was never closed");
            next();
            return e;
        }
        throw LldSyntaxError("invalid syntax");
    }
 
    static bool real(const cplx& z) { return z.imag() == 0.0; }
 
    static cplx ipow(cplx b, long e) {
        // CPython's c_powi: repeated squaring for a small integral exponent
        const bool inv = e < 0;
        unsigned long m = static_cast<unsigned long>(inv ? -e : e);
        cplx r(1.0, 0.0);
        while (m) {
            if (m & 1UL) r *= b;
            b *= b;
            m >>= 1;
        }
        return inv ? cplx(1.0, 0.0) / r : r;
    }
 
    static cplx do_pow(const cplx& a, const cplx& b) {
        if (real(a) && real(b)) {
            const double x = a.real(), y = b.real();
            if (x == 0.0 && y < 0.0)
                throw LldEvalError("ZeroDivisionError", "0.0 cannot be raised to a negative power");
            if (x >= 0.0 || y == std::floor(y)) return cplx(std::pow(x, y), 0.0);
            return std::pow(a, b);  // a negative base to a fractional power is complex in python
        }
        if (a == cplx(0.0, 0.0)) {
            if (b.imag() != 0.0 || b.real() < 0.0)
                throw LldEvalError("ZeroDivisionError", "0.0 to a negative or complex power");
            return b.real() == 0.0 ? cplx(1.0, 0.0) : cplx(0.0, 0.0);
        }
        if (b.imag() == 0.0 && b.real() == std::floor(b.real()) && std::fabs(b.real()) <= 100.0)
            return ipow(a, static_cast<long>(b.real()));
        return std::pow(a, b);
    }
 
    static cplx call(const std::string& full, const std::vector<cplx>& v) {
        std::string prefix, base = full;
        const std::size_t dot = full.rfind('.');
        if (dot != std::string::npos) {
            prefix = full.substr(0, dot);
            base = full.substr(dot + 1);
        }
        if (!prefix.empty() && prefix != "np" && prefix != "math")
            throw LldEvalError("NameError", "name '" + prefix.substr(0, prefix.find('.')) +
                                                "' is not defined");
        const bool math = prefix == "math";
        static const char* unary[] = {"exp", "log", "sqrt", "sin", "cos", "tanh", "abs"};
        bool isunary = false;
        for (const char* u : unary)
            if (base == u) isunary = true;
        if (isunary) {
            if (v.size() != 1)
                throw LldEvalError("TypeError", base + "() takes exactly one argument (" +
                                                    std::to_string(v.size()) + " given)");
            const cplx z = v[0];
            if (base == "abs") return cplx(std::abs(z), 0.0);
            if (real(z)) {
                const double x = z.real();
                if (base == "exp") return cplx(std::exp(x), 0.0);
                if (base == "sin") return cplx(std::sin(x), 0.0);
                if (base == "cos") return cplx(std::cos(x), 0.0);
                if (base == "tanh") return cplx(std::tanh(x), 0.0);
                if (base == "log") {
                    if (math && x <= 0.0) throw LldEvalError("ValueError", "math domain error");
                    if (x < 0.0) return cplx(std::numeric_limits<double>::quiet_NaN(), 0.0);
                    return cplx(std::log(x), 0.0);
                }
                // sqrt
                if (math && x < 0.0) throw LldEvalError("ValueError", "math domain error");
                if (x < 0.0) return cplx(std::numeric_limits<double>::quiet_NaN(), 0.0);
                return cplx(std::sqrt(x), 0.0);
            }
            if (math && base != "abs")
                throw LldEvalError("TypeError", "must be real number, not complex");
            if (base == "exp") return std::exp(z);
            if (base == "sin") return std::sin(z);
            if (base == "cos") return std::cos(z);
            if (base == "tanh") return std::tanh(z);
            if (base == "log") return std::log(z);
            return std::sqrt(z);
        }
        if (base == "min" || base == "max") {
            if (v.size() < 2)
                throw LldEvalError("TypeError", base + " expected at least 2 arguments, got " +
                                                    std::to_string(v.size()));
            for (const cplx& z : v)
                if (!real(z))
                    throw LldEvalError("TypeError", std::string("'") + (base == "max" ? ">" : "<") +
                                                        "' not supported between instances of "
                                                        "'complex' and 'float'");
            double r = v[0].real();
            for (std::size_t i = 1; i < v.size(); ++i)
                r = base == "max" ? (v[i].real() > r ? v[i].real() : r)
                                  : (v[i].real() < r ? v[i].real() : r);
            return cplx(r, 0.0);
        }
        if (prefix == "np") throw LldEvalError("AttributeError", "module 'numpy' has no attribute '" + base + "'");
        if (prefix == "math") throw LldEvalError("AttributeError", "module 'math' has no attribute '" + base + "'");
        throw LldEvalError("NameError", "name '" + base + "' is not defined");
    }
 
    cplx eval_node(int i, const cplx& n) const {
        const Node& nd = nodes_[std::size_t(i)];
        switch (nd.kind) {
            case Node::NUM: return cplx(nd.num, 0.0);
            case Node::VAR: return n;
            case Node::ID: throw LldEvalError("NameError", "name '" + nd.name + "' is not defined");
            case Node::NEG: return -eval_node(nd.kids[0], n);
            case Node::POS: return eval_node(nd.kids[0], n);
            case Node::ADD: return eval_node(nd.kids[0], n) + eval_node(nd.kids[1], n);
            case Node::SUB: return eval_node(nd.kids[0], n) - eval_node(nd.kids[1], n);
            case Node::MUL: {
                const cplx a = eval_node(nd.kids[0], n), b = eval_node(nd.kids[1], n);
                if (real(a) && real(b)) return cplx(a.real() * b.real(), 0.0);
                return a * b;
            }
            case Node::DIV: {
                const cplx a = eval_node(nd.kids[0], n), b = eval_node(nd.kids[1], n);
                if (b == cplx(0.0, 0.0))
                    throw LldEvalError("ZeroDivisionError", real(a) ? "float division by zero"
                                                                    : "complex division by zero");
                if (real(a) && real(b)) return cplx(a.real() / b.real(), 0.0);
                return a / b;
            }
            case Node::POW: return do_pow(eval_node(nd.kids[0], n), eval_node(nd.kids[1], n));
            case Node::CALL: {
                std::vector<cplx> args;
                for (int k : nd.kids) args.push_back(eval_node(k, n));
                return call(nd.name, args);
            }
        }
        return cplx(0.0, 0.0);
    }
};
 
inline double re(double x) { return x; }
inline double re(const cplx& z) { return z.real(); }
 
/** Smooth stand-in for min(a, b), verbatim python _SOFTMIN_SRC; real or complex. */
template <class S>
S softmin(const S& a, const S& b, double alpha) {
    S lo = a, hi = b;
    if (!(re(a) <= re(b))) {
        lo = b;
        hi = a;
    }
    const S gap = hi - lo;
    if (re(gap) * alpha > 700.0) return lo;
    using std::exp;
    const S w = exp(-alpha * gap);
    return lo + gap * w / (1.0 + w);
}
 
/**
 * Why an lld expression cannot be used, or empty: python _check_lld verbatim.
 * Real and positive at n = 1, 2.5, 10, and analytic, which a complex step
 * compared against a central difference at n = 2.5 verifies.
 */
inline std::string check_lld(const std::string& src, std::shared_ptr<const LldExpr>& out) {
    std::shared_ptr<const LldExpr> e;
    try {
        e = std::make_shared<const LldExpr>(src);
    } catch (const LldSyntaxError& ex) {
        return std::string("not a Python expression (") + ex.what() + ")";
    }
    const double pts[3] = {1.0, 2.5, 10.0};
    for (double n : pts) {
        cplx v;
        try {
            v = e->eval(cplx(n, 0.0));
        } catch (const LldEvalError& ex) {
            return "fails at n=" + fmt("%g", n) + " (" + ex.kind + ": " + ex.what() + ")";
        }
        if (v.imag() != 0.0 || !std::isfinite(v.real()) || v.real() <= 0.0)
            return "must be real and positive, is " + py_repr(v) + " at n=" + fmt("%g", n);
    }
    const double n0 = 2.5, step = 1e-6;
    double dcs = 0.0;
    try {
        dcs = e->eval(cplx(n0, 1e-20)).imag() / 1e-20;
    } catch (const LldEvalError& ex) {
        return std::string("rejects a complex argument, so it cannot be differentiated by ") +
               "complex step (" + ex.kind + ": " + ex.what() + ")";
    }
    const double dfd = (e->eval(cplx(n0 + step, 0.0)).real() - e->eval(cplx(n0 - step, 0.0)).real()) /
                       (2.0 * step);
    if (std::fabs(dcs - dfd) > 1e-4 * (1.0 + std::fabs(dfd)))
        return "is not analytic at n=" + fmt("%g", n0) + ": its complex step gives " + fmt("%g", dcs) +
               " where a finite difference gives " + fmt("%g", dfd) +
               ". abs, min, max, a comparison or an interpolation will do this.";
    out = e;
    return "";
}
 
/** One-sided Jacobi SVD of A (m x n), returning the min-norm least-squares solution of A x = b. */
inline std::vector<double> lstsq_min_norm(const Matrix<double>& A, const std::vector<double>& b) {
    const std::size_t m = A.rows(), n = A.cols();
    std::vector<double> x(n, 0.0);
    if (m == 0 || n == 0) return x;
    // column-major working copy of A, rotated until its columns are orthogonal
    std::vector<std::vector<double>> W(n, std::vector<double>(m)), V(n, std::vector<double>(n, 0.0));
    for (std::size_t j = 0; j < n; ++j) {
        for (std::size_t i = 0; i < m; ++i) W[j][i] = A(i, j);
        V[j][j] = 1.0;
    }
    const double eps = std::numeric_limits<double>::epsilon();
    for (int sweep = 0; sweep < 100; ++sweep) {
        bool rotated = false;
        for (std::size_t p = 0; p + 1 < n; ++p) {
            for (std::size_t q = p + 1; q < n; ++q) {
                double alpha = 0.0, beta = 0.0, gamma = 0.0;
                for (std::size_t i = 0; i < m; ++i) {
                    alpha += W[p][i] * W[p][i];
                    beta += W[q][i] * W[q][i];
                    gamma += W[p][i] * W[q][i];
                }
                if (gamma == 0.0 || std::fabs(gamma) <= eps * std::sqrt(alpha * beta)) continue;
                rotated = true;
                const double zeta = (beta - alpha) / (2.0 * gamma);
                const double t = (zeta >= 0.0 ? 1.0 : -1.0) / (std::fabs(zeta) + std::sqrt(1.0 + zeta * zeta));
                const double c = 1.0 / std::sqrt(1.0 + t * t), s = c * t;
                for (std::size_t i = 0; i < m; ++i) {
                    const double wp = W[p][i], wq = W[q][i];
                    W[p][i] = c * wp - s * wq;
                    W[q][i] = s * wp + c * wq;
                }
                for (std::size_t i = 0; i < n; ++i) {
                    const double vp = V[p][i], vq = V[q][i];
                    V[p][i] = c * vp - s * vq;
                    V[q][i] = s * vp + c * vq;
                }
            }
        }
        if (!rotated) break;
    }
    std::vector<double> sig(n, 0.0);
    double smax = 0.0;
    for (std::size_t j = 0; j < n; ++j) {
        double s2 = 0.0;
        for (std::size_t i = 0; i < m; ++i) s2 += W[j][i] * W[j][i];
        sig[j] = std::sqrt(s2);
        smax = std::max(smax, sig[j]);
    }
    // numpy lstsq(rcond=None): singular values at or below eps*max(m,n)*smax count as zero
    const double cut = eps * double(std::max(m, n)) * smax;
    for (std::size_t j = 0; j < n; ++j) {
        if (!(sig[j] > cut)) continue;
        double ub = 0.0;  // (w_j . b) / sigma_j^2 = (u_j . b) / sigma_j
        for (std::size_t i = 0; i < m; ++i) ub += W[j][i] * b[i];
        const double f = ub / (sig[j] * sig[j]);
        for (std::size_t i = 0; i < n; ++i) x[i] += f * V[j][i];
    }
    return x;
}
 
inline double norm2(const std::vector<double>& v) {
    double s = 0.0;
    for (double x : v) s += x * x;
    return std::sqrt(s);
}
 
}  // namespace nlp_detail
 
/**
 * E[max] of K independent EXPONENTIAL branch completion times of means m, with
 * the reciprocals cleared so a branch of zero mean contributes zero:
 * sum over nonempty S of (-1)^(|S|+1) prod_{i in S} m_i / sum_{i in S} prod_{j in S, j != i} m_j.
 */
template <class S>
S fj_max_exp(const std::vector<S>& m) {
    const std::size_t k = m.size();
    S tot = S(0.0);
    for (std::size_t mask = 1; mask < (std::size_t(1) << k); ++mask) {
        S num = S(1.0), den = S(0.0);
        std::size_t card = 0;
        for (std::size_t i = 0; i < k; ++i) {
            if (!((mask >> i) & 1U)) continue;
            ++card;
            num = num * m[i];
            S q = S(1.0);
            for (std::size_t j = 0; j < k; ++j)
                if (j != i && ((mask >> j) & 1U)) q = q * m[j];
            den = den + q;
        }
        tot = tot + (card % 2 ? 1.0 : -1.0) * num / den;
    }
    return tot;
}
 
/** Gradient of fj_max_exp by complex step, one coordinate at a time. */
inline std::vector<double> fj_max_exp_grad(const std::vector<double>& u) {
    std::vector<double> g(u.size(), 0.0);
    for (std::size_t j = 0; j < u.size(); ++j) {
        std::vector<nlp_detail::cplx> z(u.begin(), u.end());
        z[j] += nlp_detail::cplx(0.0, 1e-20);
        g[j] = fj_max_exp(z).imag() / 1e-20;
    }
    return g;
}
 
/** The queue-dependent rate of one station, r(n) = softmin(n, m) [* lld(n)], or n at a delay. */
struct NlpRate {
    bool delay = false;
    double mult = 1.0;
    double alpha = 20.0;
    std::shared_ptr<const nlp_detail::LldExpr> lld;
 
    double operator()(double n) const {
        if (delay) return n;
        double s = nlp_detail::softmin<double>(n, mult, alpha);
        if (lld) s = s * lld->eval(nlp_detail::cplx(n, 0.0)).real();
        return s;
    }
    nlp_detail::cplx operator()(const nlp_detail::cplx& n) const {
        if (delay) return n;
        nlp_detail::cplx s = nlp_detail::softmin<nlp_detail::cplx>(n, nlp_detail::cplx(mult, 0.0), alpha);
        if (lld) s = s * lld->eval(n);
        return s;
    }
    /** Derivative by complex step: imag(f(n + 1e-20 i)) / 1e-20. */
    double deriv(double n) const { return (*this)(nlp_detail::cplx(n, 1e-20)).imag() / 1e-20; }
};
 
/** What NlpProgram::solve returns: python's Result. */
struct NlpSolveResult {
    std::vector<double> x;
    std::string message;
    int niter = 0;          ///< fallback iterations; 0 when the Newton ladder rooted it
    bool newton = true;     ///< the Newton ladder reached a feasible root
};
 
/**
 * The program of one layered model. Every member is named after the python
 * script's global of the same meaning (lower case); indices are 0-based in
 * declaration order, so element i of a block is the LQN's absolute index
 * shift + i + 1.
 */
class NlpProgram {
public:
    enum Block { XE = 0, XA, XC, W, S, R, Q, QT, BC, BR, GH, GT, SA, NBLOCKS };
 
    // ---- model parameters ------------------------------------------------
    std::string name;
    std::vector<std::string> hosts, tasks, entries, acts, calls;
    std::vector<double> host_mult, host_maxmult;
    std::vector<bool> host_inf;
    std::vector<std::size_t> task_host;
    std::vector<double> task_mult, task_think;
    std::vector<bool> task_inf, task_isref;
    std::vector<std::size_t> entry_task, entry_order;
    std::vector<std::size_t> act_task, act_entry;
    std::vector<double> act_demand, act_visits;
    std::vector<std::size_t> call_act, call_entry, call_caller_task;
    std::vector<double> call_mean;
    std::vector<std::size_t> fork_entry;
    std::vector<std::vector<std::vector<std::size_t>>> fork_branches;
    std::vector<std::string> fork_label;
    double alpha = 20.0;
    std::vector<std::string> host_lld, task_lld;
    std::vector<bool> host_rate_capped;
    std::vector<NlpRate> host_rate, task_rate;
 
    // ---- derived index sets ----------------------------------------------
    std::size_t NH = 0, NT = 0, NE = 0, NA = 0, NC = 0, NFORK = 0, NR = 0, NGH = 0, NGT = 0;
    std::vector<std::vector<std::size_t>> acts_of_host, acts_of_entry, calls_of_act,
        calls_into_entry, calls_into_task, entries_of_task, blocks_of_entry, seq_acts_of_entry,
        sibling_acts, sibling_calls;
    std::set<std::size_t> fork_acts;
    std::vector<std::size_t> ref_tasks, ref_entry, qd_hosts, qd_tasks;
    std::map<std::size_t, std::size_t> gh_of, gt_of;
    std::vector<bool> task_delay, host_delay;
    std::vector<double> task_pop, host_ntot, task_ntot, host_delta, task_delta;
 
    // ---- layout ----------------------------------------------------------
    std::size_t offset[NBLOCKS] = {0};
    std::size_t blocksize[NBLOCKS] = {0};
    std::size_t NVAR = 0, NRES = 0;
 
    // ---- linear block ----------------------------------------------------
    Matrix<double> AEQ, AUB;
    std::vector<double> BEQ, BUB;
 
    // ---- residuals: const + lin.x + bilinear + rate terms + fork terms ----
    struct Bil { std::size_t k, i, j; double c; };
    struct Nlt { std::size_t k; const NlpRate* f; double s, c0; std::vector<double> a; };
    struct Mnl { std::size_t k; double s; std::vector<double> c; Matrix<double> a; };
    Matrix<double> res_lin;
    std::vector<double> res_const;
    std::vector<Bil> bil;
    std::vector<Nlt> nlt;
    std::vector<Mnl> mnl;
 
    NlpProgram() {}
    NlpProgram(const NlpProgram& o) { *this = o; }
    NlpProgram& operator=(const NlpProgram& o) {
        if (this == &o) return *this;
        copy_fields(o);
        // the rate terms point into this program's own rate tables
        rebind_rates(o);
        return *this;
    }
 
    std::size_t col(Block b, std::size_t i) const { return offset[b] + i; }
    std::size_t sent(std::size_t e) const { return NFORK ? col(SA, e) : col(S, e); }
 
    /** The block's slice of an iterate. */
    std::vector<double> block(const std::vector<double>& x, Block b) const {
        return std::vector<double>(x.begin() + long(offset[b]), x.begin() + long(offset[b] + blocksize[b]));
    }
 
    std::vector<std::string> var_names() const {
        static const char* tag[NBLOCKS] = {"Xe", "Xa", "Xc", "W", "S", "R", "Q", "Qt", "BC", "BR", "GH", "GT", "SA"};
        std::vector<std::string> out;
        for (int b = 0; b < NBLOCKS; ++b) {
            for (std::size_t i = 0; i < blocksize[b]; ++i) {
                std::string lab;
                switch (Block(b)) {
                    case XE: case S: case SA: lab = entries[i]; break;
                    case XA: case W: case Q: lab = acts[i]; break;
                    case XC: case R: case QT: case BC: lab = calls[i]; break;
                    case BR: lab = tasks[ref_tasks[i]]; break;
                    case GH: lab = hosts[qd_hosts[i]]; break;
                    case GT: lab = tasks[qd_tasks[i]]; break;
                    default: break;
                }
                out.push_back(std::string(tag[b]) + "[" + lab + "]");
            }
        }
        return out;
    }
 
    // ---- QD-AMVA arrival-instant estimate --------------------------------
    std::vector<std::pair<std::size_t, double>> qbar_host_terms(std::size_t h, std::size_t a, double g) const {
        std::vector<std::pair<std::size_t, double>> out;
        for (std::size_t ap : acts_of_host[h]) out.push_back({col(Q, ap), g});
        const double pop = task_pop[act_task[a]];
        if (!std::isinf(pop) && pop > 0)
            for (std::size_t ap : sibling_acts[a]) out.push_back({col(Q, ap), -g / pop});
        return out;
    }
    std::vector<std::pair<std::size_t, double>> qbar_task_terms(std::size_t t, std::size_t c, double g) const {
        std::vector<std::pair<std::size_t, double>> out;
        for (std::size_t cp : calls_into_task[t]) out.push_back({col(QT, cp), g});
        const double pop = task_pop[call_caller_task[c]];
        if (!std::isinf(pop) && pop > 0)
            for (std::size_t cp : sibling_calls[c]) out.push_back({col(QT, cp), -g / pop});
        return out;
    }
    double qd_pop_host(const std::vector<double>& Qv, std::size_t h) const {
        double s = 0.0;
        for (std::size_t a : acts_of_host[h]) s += Qv[a];
        return 1.0 + host_delta[h] * s;
    }
    double qd_pop_task(const std::vector<double>& Qtv, std::size_t t) const {
        double s = 0.0;
        for (std::size_t c : calls_into_task[t]) s += Qtv[c];
        return 1.0 + task_delta[t] * s;
    }
    double qbar_host(const std::vector<double>& Qv, std::size_t h, std::size_t a) const {
        const double pop = task_pop[act_task[a]];
        double tot = 0.0, own = 0.0;
        for (std::size_t ap : acts_of_host[h]) tot += Qv[ap];
        for (std::size_t ap : sibling_acts[a]) own += Qv[ap];
        return tot - ((!std::isinf(pop) && pop > 0) ? own / pop : 0.0);
    }
    double qbar_task(const std::vector<double>& Qtv, std::size_t t, std::size_t c) const {
        const double pop = task_pop[call_caller_task[c]];
        double tot = 0.0, own = 0.0;
        for (std::size_t cp : calls_into_task[t]) tot += Qtv[cp];
        for (std::size_t cp : sibling_calls[c]) own += Qtv[cp];
        return tot - ((!std::isinf(pop) && pop > 0) ? own / pop : 0.0);
    }
 
    // ---- evaluation ------------------------------------------------------
    std::vector<double> residuals(const std::vector<double>& x) const {
        std::vector<double> r(NRES, 0.0);
        for (std::size_t k = 0; k < NRES; ++k) {
            double s = res_const[k];
            for (std::size_t j = 0; j < NVAR; ++j) s += res_lin(k, j) * x[j];
            r[k] = s;
        }
        for (const Bil& b : bil) r[b.k] += b.c * x[b.i] * x[b.j];
        for (const Nlt& t : nlt) r[t.k] += t.s * (*t.f)(t.c0 + dot(t.a, x));
        for (const Mnl& t : mnl) r[t.k] += t.s * fj_max_exp(mnl_args(t, x));
        return r;
    }
 
    Matrix<double> jacobian(const std::vector<double>& x) const {
        Matrix<double> J = res_lin;
        for (const Bil& b : bil) {
            J(b.k, b.i) += b.c * x[b.j];
            J(b.k, b.j) += b.c * x[b.i];
        }
        for (const Nlt& t : nlt) {
            const double d = t.s * t.f->deriv(t.c0 + dot(t.a, x));
            for (std::size_t j = 0; j < NVAR; ++j) J(t.k, j) += d * t.a[j];
        }
        for (const Mnl& t : mnl) {
            const std::vector<double> g = fj_max_exp_grad(mnl_args(t, x));
            for (std::size_t j = 0; j < NVAR; ++j) {
                double s = 0.0;
                for (std::size_t i = 0; i < g.size(); ++i) s += g[i] * t.a(i, j);
                J(t.k, j) += t.s * s;
            }
        }
        return J;
    }
 
    double objective(const std::vector<double>& x) const {
        const std::vector<double> r = residuals(x);
        double s = 0.0;
        for (double v : r) s += v * v;
        return s;
    }
 
    std::vector<double> aeq_resid(const std::vector<double>& x) const { return affine(AEQ, BEQ, x); }
    std::vector<double> aub_resid(const std::vector<double>& x) const { return affine(AUB, BUB, x); }
 
    // ---- starting point --------------------------------------------------
    double compose(const std::vector<std::size_t>& as, const std::vector<double>& Wv,
                   const std::vector<double>& Rv) const {
        double s = 0.0;
        for (std::size_t a : as) {
            s += act_visits[a] * Wv[a];
            for (std::size_t c : calls_of_act[a]) s += act_visits[a] * call_mean[c] * Rv[c];
        }
        return s;
    }
 
    struct Sweep { std::vector<double> W, S, SA, R, Xe, Xa, Xc; };
 
    /** One pass of the relations at the given queue lengths: python _sweep. */
    Sweep sweep(const std::vector<double>& Qv, const std::vector<double>& Qtv) const {
        Sweep o;
        o.W.assign(NA, 0.0);
        for (std::size_t h = 0; h < NH; ++h)
            for (std::size_t a : acts_of_host[h])
                o.W[a] = host_delay[h] ? act_demand[a]
                                       : act_demand[a] * (1.0 + qbar_host(Qv, h, a)) / host_rate[h](qd_pop_host(Qv, h));
        o.S.assign(NE, 0.0);
        o.SA.assign(NE, 0.0);
        o.R.assign(NC, 0.0);
        for (std::size_t e : entry_order) {
            o.SA[e] = compose(acts_of_entry[e], o.W, o.R);
            if (!blocks_of_entry[e].empty()) {
                double s = compose(seq_acts_of_entry[e], o.W, o.R);
                for (std::size_t b : blocks_of_entry[e]) {
                    std::vector<double> m;
                    for (const std::vector<std::size_t>& br : fork_branches[b]) m.push_back(compose(br, o.W, o.R));
                    s += fj_max_exp(m);
                }
                o.S[e] = s;
            } else {
                o.S[e] = o.SA[e];
            }
            const std::size_t t = entry_task[e];
            for (std::size_t c : calls_into_entry[e])
                o.R[c] = task_delay[t] ? o.S[e]
                                       : o.S[e] * (1.0 + qbar_task(Qtv, t, c)) / task_rate[t](qd_pop_task(Qtv, t));
        }
        o.Xe.assign(NE, 0.0);
        o.Xa.assign(NA, 0.0);
        o.Xc.assign(NC, 0.0);
        for (std::size_t ri = 0; ri < NR; ++ri) {
            const std::size_t t = ref_tasks[ri], e = ref_entry[ri];
            const double cyc = task_think[t] + o.S[e];
            o.Xe[e] = cyc > 0 ? task_mult[t] / cyc : 0.0;
        }
        for (std::size_t q = entry_order.size(); q-- > 0;) {
            const std::size_t e = entry_order[q];
            if (!task_isref[entry_task[e]]) {
                double s = 0.0;
                for (std::size_t c : calls_into_entry[e]) s += o.Xc[c];
                o.Xe[e] = s;
            }
            for (std::size_t a : acts_of_entry[e]) {
                o.Xa[a] = act_visits[a] * o.Xe[e];
                for (std::size_t c : calls_of_act[a]) o.Xc[c] = call_mean[c] * o.Xa[a];
            }
        }
        return o;
    }
 
    /** A point satisfying every linear relation, after WARM fixed-point sweeps. */
    std::vector<double> initial_point(int warm = 0) const {
        std::vector<double> Qv(NA, 0.0), Qtv(NC, 0.0);
        for (int it = 0; it < std::max(0, warm) + 1; ++it) {
            const Sweep s = sweep(Qv, Qtv);
            for (std::size_t a = 0; a < NA; ++a) Qv[a] = s.Xa[a] * s.W[a];
            for (std::size_t c = 0; c < NC; ++c) Qtv[c] = s.Xc[c] * s.R[c];
        }
        const Sweep s = sweep(Qv, Qtv);
        std::vector<double> x(NVAR, 0.0);
        put(x, XE, s.Xe);
        put(x, XA, s.Xa);
        put(x, XC, s.Xc);
        put(x, W, s.W);
        put(x, S, s.S);
        if (NFORK) put(x, SA, s.SA);
        put(x, R, s.R);
        if (warm) {
            put(x, Q, Qv);
            put(x, QT, Qtv);
        }
        for (std::size_t c = 0; c < NC; ++c) x[col(BC, c)] = s.Xc[c] * s.SA[call_entry[c]];
        for (std::size_t ri = 0; ri < NR; ++ri) {
            const std::size_t e = ref_entry[ri];
            x[col(BR, ri)] = s.Xe[e] * s.S[e];
        }
        const std::vector<double> Qs = block(x, Q), Qts = block(x, QT);
        for (std::size_t h : qd_hosts) x[col(GH, gh_of.at(h))] = host_rate[h](qd_pop_host(Qs, h));
        for (std::size_t t : qd_tasks) x[col(GT, gt_of.at(t))] = task_rate[t](qd_pop_task(Qts, t));
        return x;
    }
 
    // ---- solve -----------------------------------------------------------
    /** Newton on the square system [residuals; AEQ x - BEQ], min-norm step: python polish. */
    std::vector<double> polish(std::vector<double> x, int iters = 50) const {
        std::vector<double> best = x;
        double bobj = objective(x);
        const double scale = 1.0 + maxabs(x);
        for (int it = 0; it < iters; ++it) {
            const std::vector<double> r = residuals(x), e = aeq_resid(x);
            const Matrix<double> J = jacobian(x);
            Matrix<double> JF(NRES + AEQ.rows(), NVAR, 0.0);
            std::vector<double> F(NRES + AEQ.rows(), 0.0);
            for (std::size_t k = 0; k < NRES; ++k) {
                F[k] = -r[k];
                for (std::size_t j = 0; j < NVAR; ++j) JF(k, j) = J(k, j);
            }
            for (std::size_t k = 0; k < AEQ.rows(); ++k) {
                F[NRES + k] = -e[k];
                for (std::size_t j = 0; j < NVAR; ++j) JF(NRES + k, j) = AEQ(k, j);
            }
            const std::vector<double> dx = nlp_detail::lstsq_min_norm(JF, F);
            for (std::size_t j = 0; j < NVAR; ++j) x[j] += dx[j];
            if (minval(x) < -1e-6 * scale) break;
            if (AUB.rows() && maxval(aub_resid(x)) > 1e-6 * scale) break;
            const double obj = objective(x);
            if (obj < bobj) {
                best = x;
                bobj = obj;
            }
            if (nlp_detail::norm2(dx) <= 1e-14 * (1.0 + nlp_detail::norm2(x))) break;
        }
        for (double& v : best) v = std::max(v, 0.0);
        return best;
    }
 
    /** A root of the square system, feasible to the bounds and the inequalities. */
    bool solved(const std::vector<double>& x) const {
        const double scale = 1.0 + (x.empty() ? 0.0 : maxabs(x));
        if (!x.empty() && minval(x) < -1e-9 * scale) return false;
        if (AUB.rows() && maxval(aub_resid(x)) > 1e-7 * scale) return false;
        if (AEQ.rows()) {
            const std::vector<double> e = aeq_resid(x);
            if (maxabs(e) > 1e-9 * scale) return false;
        }
        return NRES ? maxabs(residuals(x)) <= 1e-10 * scale : true;
    }
 
    /** The merit the fallback descends: ||r||^2 + ||AEQ x - BEQ||^2 + ||max(0, AUB x - BUB)||^2. */
    double merit(const std::vector<double>& x) const {
        double s = objective(x);
        for (double v : aeq_resid(x)) s += v * v;
        for (double v : aub_resid(x))
            if (v > 0) s += v * v;
        return s;
    }
 
    /**
     * The projected Levenberg-Marquardt fallback, shared by MATLAB, the JAR and
     * C++ in place of python's trust-constr. `merits`, when given, receives the
     * merit after each ACCEPTED step, which is how a test sees it descend.
     */
    std::vector<double> lm_fallback(std::vector<double> x, int maxiter, int* iters = nullptr,
                                    std::vector<double>* merits = nullptr) const {
        for (double& v : x) v = std::max(v, 0.0);
        double phi = merit(x);
        double mu = -1.0;
        int it = 0;
        while (it < maxiter && phi > 1e-30) {
            const std::vector<double> r = residuals(x), e = aeq_resid(x), u = aub_resid(x);
            const Matrix<double> J = jacobian(x);
            std::vector<std::vector<double>> rows;
            std::vector<double> G;
            for (std::size_t k = 0; k < NRES; ++k) {
                rows.push_back(std::vector<double>(&J(k, 0), &J(k, 0) + NVAR));
                G.push_back(r[k]);
            }
            for (std::size_t k = 0; k < AEQ.rows(); ++k) {
                rows.push_back(std::vector<double>(&AEQ(k, 0), &AEQ(k, 0) + NVAR));
                G.push_back(e[k]);
            }
            for (std::size_t k = 0; k < AUB.rows(); ++k)
                if (u[k] > 0) {
                    rows.push_back(std::vector<double>(&AUB(k, 0), &AUB(k, 0) + NVAR));
                    G.push_back(u[k]);
                }
            Matrix<double> JtJ(NVAR, NVAR, 0.0);
            std::vector<double> JtG(NVAR, 0.0);
            for (std::size_t q = 0; q < rows.size(); ++q)
                for (std::size_t i = 0; i < NVAR; ++i) {
                    if (rows[q][i] == 0.0) continue;
                    JtG[i] += rows[q][i] * G[q];
                    for (std::size_t j = 0; j < NVAR; ++j) JtJ(i, j) += rows[q][i] * rows[q][j];
                }
            if (mu < 0) {
                double dmax = 0.0;
                for (std::size_t i = 0; i < NVAR; ++i) dmax = std::max(dmax, JtJ(i, i));
                mu = 1e-3 * dmax;
            }
            bool stop = false;
            while (it < maxiter) {
                ++it;
                Matrix<double> M = JtJ;
                for (std::size_t i = 0; i < NVAR; ++i) M(i, i) += mu;
                std::vector<double> rhs(NVAR);
                for (std::size_t i = 0; i < NVAR; ++i) rhs[i] = -JtG[i];
                const std::vector<double> d = line::solve(M, rhs);
                std::vector<double> xn(NVAR);
                for (std::size_t i = 0; i < NVAR; ++i) xn[i] = std::max(x[i] + d[i], 0.0);
                const double pn = merit(xn);
                if (nlp_detail::norm2(d) <= 1e-14 * (1.0 + nlp_detail::norm2(x))) stop = true;
                if (pn < phi) {
                    x = xn;
                    phi = pn;
                    mu /= 3.0;
                    if (merits) merits->push_back(phi);
                    break;
                }
                mu *= 2.0;
                if (stop) break;
            }
            if (stop) break;
        }
        if (iters) *iters = it;
        return x;
    }
 
    /**
     * Root the program: the Newton ladder over warm starts max(0,warm), 5, 20,
     * 100, then the fallback from the ladder start of smallest objective.
     */
    NlpSolveResult solve(int warm = 0, int maxiter = 3000) const {
        std::vector<int> ladder{std::max(0, warm)};
        for (int w : {5, 20, 100})
            if (w > warm) ladder.push_back(w);
        std::vector<double> best_x0;
        double best_obj = std::numeric_limits<double>::infinity();
        for (int w : ladder) {
            const std::vector<double> x0 = initial_point(w);
            const std::vector<double> x = polish(x0);
            if (solved(x)) {
                NlpSolveResult res;
                res.x = x;
                res.message = "Newton polish reached a feasible root from " + std::to_string(w) +
                              " warm sweep" + (w == 1 ? "" : "s") + ".";
                return res;
            }
            const double obj = objective(x0);
            if (obj < best_obj) {
                best_x0 = x0;
                best_obj = obj;
            }
        }
        if (best_x0.empty()) best_x0 = initial_point(ladder.front());
        NlpSolveResult res;
        res.newton = false;
        res.x = polish(lm_fallback(best_x0, maxiter, &res.niter));
        res.message = "projected Levenberg-Marquardt fallback, " + std::to_string(res.niter) +
                      " iterations, then Newton polish; merit " + nlp_detail::fmt("%.3e", merit(res.x)) + ".";
        return res;
    }
 
    // ---- helpers ---------------------------------------------------------
    static double dot(const std::vector<double>& a, const std::vector<double>& x) {
        double s = 0.0;
        for (std::size_t j = 0; j < a.size(); ++j) s += a[j] * x[j];
        return s;
    }
    static double maxabs(const std::vector<double>& v) {
        double m = 0.0;
        for (double x : v) m = std::max(m, std::fabs(x));
        return m;
    }
    static double minval(const std::vector<double>& v) {
        double m = std::numeric_limits<double>::infinity();
        for (double x : v) m = std::min(m, x);
        return m;
    }
    static double maxval(const std::vector<double>& v) {
        double m = -std::numeric_limits<double>::infinity();
        for (double x : v) m = std::max(m, x);
        return m;
    }
 
private:
    std::vector<double> mnl_args(const Mnl& t, const std::vector<double>& x) const {
        std::vector<double> u(t.c);
        for (std::size_t i = 0; i < u.size(); ++i)
            for (std::size_t j = 0; j < NVAR; ++j) u[i] += t.a(i, j) * x[j];
        return u;
    }
    std::vector<double> affine(const Matrix<double>& A, const std::vector<double>& b,
                               const std::vector<double>& x) const {
        std::vector<double> out(A.rows(), 0.0);
        for (std::size_t k = 0; k < A.rows(); ++k) {
            double s = 0.0;
            for (std::size_t j = 0; j < NVAR; ++j) s += A(k, j) * x[j];
            out[k] = s - b[k];
        }
        return out;
    }
    void put(std::vector<double>& x, Block b, const std::vector<double>& v) const {
        for (std::size_t i = 0; i < blocksize[b]; ++i) x[offset[b] + i] = v[i];
    }
    void copy_fields(const NlpProgram& o);
    void rebind_rates(const NlpProgram& o) {
        for (Nlt& t : nlt) {
            for (std::size_t h = 0; h < o.host_rate.size(); ++h)
                if (t.f == &o.host_rate[h]) t.f = &host_rate[h];
            for (std::size_t q = 0; q < o.task_rate.size(); ++q)
                if (t.f == &o.task_rate[q]) t.f = &task_rate[q];
        }
    }
};
 
inline void NlpProgram::copy_fields(const NlpProgram& o) {
    name = o.name; hosts = o.hosts; tasks = o.tasks; entries = o.entries; acts = o.acts; calls = o.calls;
    host_mult = o.host_mult; host_maxmult = o.host_maxmult; host_inf = o.host_inf;
    task_host = o.task_host; task_mult = o.task_mult; task_think = o.task_think;
    task_inf = o.task_inf; task_isref = o.task_isref; entry_task = o.entry_task;
    entry_order = o.entry_order; act_task = o.act_task; act_entry = o.act_entry;
    act_demand = o.act_demand; act_visits = o.act_visits; call_act = o.call_act;
    call_entry = o.call_entry; call_caller_task = o.call_caller_task; call_mean = o.call_mean;
    fork_entry = o.fork_entry; fork_branches = o.fork_branches; fork_label = o.fork_label;
    alpha = o.alpha; host_lld = o.host_lld; task_lld = o.task_lld;
    host_rate_capped = o.host_rate_capped; host_rate = o.host_rate; task_rate = o.task_rate;
    NH = o.NH; NT = o.NT; NE = o.NE; NA = o.NA; NC = o.NC; NFORK = o.NFORK; NR = o.NR;
    NGH = o.NGH; NGT = o.NGT;
    acts_of_host = o.acts_of_host; acts_of_entry = o.acts_of_entry; calls_of_act = o.calls_of_act;
    calls_into_entry = o.calls_into_entry; calls_into_task = o.calls_into_task;
    entries_of_task = o.entries_of_task; blocks_of_entry = o.blocks_of_entry;
    seq_acts_of_entry = o.seq_acts_of_entry; sibling_acts = o.sibling_acts;
    sibling_calls = o.sibling_calls; fork_acts = o.fork_acts; ref_tasks = o.ref_tasks;
    ref_entry = o.ref_entry; qd_hosts = o.qd_hosts; qd_tasks = o.qd_tasks; gh_of = o.gh_of;
    gt_of = o.gt_of; task_delay = o.task_delay; host_delay = o.host_delay; task_pop = o.task_pop;
    host_ntot = o.host_ntot; task_ntot = o.task_ntot; host_delta = o.host_delta;
    task_delta = o.task_delta;
    for (int b = 0; b < NBLOCKS; ++b) {
        offset[b] = o.offset[b];
        blocksize[b] = o.blocksize[b];
    }
    NVAR = o.NVAR; NRES = o.NRES; AEQ = o.AEQ; AUB = o.AUB; BEQ = o.BEQ; BUB = o.BUB;
    res_lin = o.res_lin; res_const = o.res_const; bil = o.bil; nlt = o.nlt; mnl = o.mnl;
}
 
namespace nlp_detail {
 
inline UnsupportedError nlp_error(const std::string& msg) {
    return UnsupportedError(ErrorKind{"NlpExport"}, msg);
}
 
template <class T>
std::string elem(const LqnStruct<T>& lqn, std::size_t idx) {
    if (idx < lqn.names.size() && !lqn.names[idx].empty()) return lqn.names[idx];
    return "#" + std::to_string(idx);
}
 
template <class T>
std::string callnm(const LqnStruct<T>& lqn, std::size_t c) {
    if (c < lqn.callnames.size() && !lqn.callnames[c].empty()) return lqn.callnames[c];
    return "C" + std::to_string(c);
}
 
template <class T>
double hostdem_mean(const LqnStruct<T>& lqn, std::size_t a) {
    if (a >= lqn.hostdem.size() || lqn.hostdem[a].disabled) return 0.0;
    const double d = num_traits<T>::to_double(lqn.hostdem[a].mean);
    return std::isnan(d) ? 0.0 : d;
}
 
/** One AND-fork/join block, in absolute LQN indices. */
struct ForkBlock {
    std::size_t entry = 0, join = 0;
    std::vector<std::size_t> tails;
    std::vector<std::vector<std::size_t>> branches;
    std::string label;
};
 
template <class T>
bool is_post_and(const LqnStruct<T>& lqn, std::size_t a) {
    return a < lqn.actposttype.size() && lqn.actposttype[a] == lang::PrecedenceType::POST_AND;
}
template <class T>
bool is_pre_and(const LqnStruct<T>& lqn, std::size_t a) {
    return a < lqn.actpretype.size() && lqn.actpretype[a] == lang::PrecedenceType::PRE_AND;
}
template <class T>
std::size_t quorum_of(const LqnStruct<T>& lqn, std::size_t a) {
    return a < lqn.actquorum.size() ? lqn.actquorum[a] : 0;
}
 
/** python _read_block: one candidate block, or the reason it is not one. */
template <class T>
std::string read_block(const LqnStruct<T>& lqn, const std::vector<std::size_t>& A,
                       const std::vector<std::size_t>& H, const std::vector<std::size_t>& srcs,
                       const std::map<std::size_t, std::vector<std::size_t>>& succ,
                       const std::map<std::size_t, std::vector<std::size_t>>& prd, ForkBlock& blk) {
    auto nm = [&](std::size_t a) { return elem(lqn, a); };
    auto reach = [&](std::size_t start) {
        std::set<std::size_t> seen;
        std::vector<std::size_t> stack{start};
        while (!stack.empty()) {
            const std::size_t x = stack.back();
            stack.pop_back();
            if (seen.count(x)) continue;
            seen.insert(x);
            for (std::size_t y : succ.at(x)) stack.push_back(y);
        }
        return seen;
    };
    std::vector<std::size_t> Hs = H;
    std::sort(Hs.begin(), Hs.end());
    std::string lbl;
    for (std::size_t i = 0; i < Hs.size(); ++i) lbl += (i ? "+" : "") + nm(Hs[i]);
    if (H.size() < 2) return "AND-fork at " + lbl + " has one branch";
 
    std::vector<std::set<std::size_t>> Rs;
    for (std::size_t h : H) Rs.push_back(reach(h));
    std::set<std::size_t> common = Rs[0];
    for (std::size_t i = 1; i < Rs.size(); ++i) {
        std::set<std::size_t> c;
        for (std::size_t x : common)
            if (Rs[i].count(x)) c.insert(x);
        common = c;
    }
    if (common.empty())
        return "the AND-fork branches at " + lbl + " never reconverge, so the entry has no single completion";
    std::vector<std::size_t> cands;
    for (std::size_t x : common) {
        if (quorum_of(lqn, x) == 0) continue;
        const std::set<std::size_t> rx = reach(x);
        bool sub = true;
        for (std::size_t y : common)
            if (!rx.count(y)) sub = false;
        if (sub) cands.push_back(x);
    }
    if (cands.size() != 1) return "the AND-fork branches at " + lbl + " do not meet at one join target";
    const std::size_t j = cands[0];
 
    std::vector<std::size_t> tails;
    const std::vector<std::size_t>& pj = prd.at(j);
    for (std::size_t a : A)
        if (is_pre_and(lqn, a) && std::find(pj.begin(), pj.end(), a) != pj.end()) tails.push_back(a);
    if (tails.size() != H.size())
        return "the join " + nm(j) + " takes " + std::to_string(tails.size()) +
               " AND predecessors where the fork at " + lbl + " opens " + std::to_string(H.size()) +
               " branches";
    const std::size_t q = quorum_of(lqn, j);
    if (q != H.size())
        return "quorum join at " + nm(j) + " fires on " + std::to_string(q) + " of " +
               std::to_string(H.size()) + " branches; the k-th order statistic is a different closure";
 
    std::vector<std::vector<std::size_t>> branches;
    for (const std::set<std::size_t>& r : Rs) {
        std::vector<std::size_t> br;
        for (std::size_t x : r)
            if (!common.count(x)) br.push_back(x);
        branches.push_back(br);
    }
    std::set<std::size_t> seen;
    for (std::size_t i = 0; i < branches.size(); ++i) {
        if (branches[i].empty()) return "branch " + nm(H[i]) + " of the AND-fork at " + lbl + " is empty";
        for (std::size_t a : branches[i])
            if (seen.count(a)) return "the AND-fork branches at " + lbl + " overlap, so they are not independent";
        for (std::size_t a : branches[i]) seen.insert(a);
    }
    for (std::size_t i = 0; i < branches.size(); ++i) {
        const std::vector<std::size_t>& br = branches[i];
        for (std::size_t a : br) {
            for (std::size_t b : prd.at(a))
                if (std::find(br.begin(), br.end(), b) == br.end() &&
                    std::find(srcs.begin(), srcs.end(), b) == srcs.end())
                    return "activity " + nm(a) + " is inside a branch of the AND-fork at " + lbl +
                           " but is also entered from " + nm(b) + ", so the branch is not a block";
            if (is_post_and(lqn, a) && a != H[i])
                return "AND-fork at " + nm(a) + " nested inside a branch of the one at " + lbl;
            if (is_pre_and(lqn, a) && std::find(tails.begin(), tails.end(), a) == tails.end())
                return "AND-join at " + nm(a) + " nested inside a branch of the fork at " + lbl;
        }
        bool work = false;
        for (std::size_t a : br) {
            if (hostdem_mean(lqn, a) > 0) work = true;
            for (std::size_t c = 1; c <= lqn.ncalls; ++c)
                if (lqn.callpair_src[c] == a) work = true;
        }
        if (!work)
            return "branch " + nm(H[i]) + " of the AND-fork at " + lbl +
                   " has no host demand and no call, so its completion time is identically zero";
    }
    blk.join = j;
    blk.tails = tails;
    blk.branches = branches;
    return "";
}
 
/** python _fork_blocks: the well-formed blocks, and every AND marking outside one. */
template <class T>
std::vector<ForkBlock> fork_blocks(const LqnStruct<T>& lqn, const std::string& fork,
                                   std::vector<std::string>& bad) {
    std::vector<ForkBlock> blocks;
    const std::size_t a0 = lqn.ashift + 1, a1 = lqn.ashift + lqn.nacts;
    std::vector<std::size_t> marked;
    for (std::size_t a = a0; a <= a1; ++a)
        if (is_post_and(lqn, a) || is_pre_and(lqn, a)) marked.push_back(a);
    if (marked.empty()) return blocks;
    if (fork != "exp") {
        for (std::size_t a : marked)
            bad.push_back("AND-fork or AND-join at activity " + elem(lqn, a) +
                          " (pass fork='exp' to carry it under the exponential-branch closure)");
        return blocks;
    }
    const auto G = [&](std::size_t i, std::size_t j) {
        return i <= lqn.graph.n && j <= lqn.graph.n && num_traits<T>::to_double(lqn.graph.get(i, j)) != 0.0;
    };
    std::set<std::size_t> placed;
    for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
        const std::size_t tidx = lqn.tshift + t;
        std::vector<std::size_t> ents = tidx < lqn.entriesof.size() ? lqn.entriesof[tidx] : std::vector<std::size_t>();
        std::sort(ents.begin(), ents.end());
        for (std::size_t eidx : ents) {
            std::set<std::size_t> As;
            if (eidx < lqn.actsof.size())
                for (std::size_t a : lqn.actsof[eidx])
                    if (a >= a0 && a <= a1) As.insert(a);
            const std::vector<std::size_t> A(As.begin(), As.end());
            if (A.empty()) continue;
            std::map<std::size_t, std::vector<std::size_t>> succ, prd;
            for (std::size_t a : A) {
                succ[a];
                prd[a];
                for (std::size_t b : A) {
                    if (G(a, b)) succ[a].push_back(b);
                    if (G(b, a)) prd[a].push_back(b);
                }
            }
            std::map<std::vector<std::size_t>, std::vector<std::size_t>> groups;
            for (std::size_t a : A)
                if (is_post_and(lqn, a)) {
                    std::vector<std::size_t> key = prd[a];
                    std::sort(key.begin(), key.end());
                    groups[key].push_back(a);
                }
            for (const auto& g : groups) {
                ForkBlock blk;
                const std::string why = read_block(lqn, A, g.second, g.first, succ, prd, blk);
                if (!why.empty()) {
                    bad.push_back(why);
                    continue;
                }
                blk.entry = eidx;
                std::string src;
                for (std::size_t i = 0; i < g.first.size(); ++i) src += (i ? "+" : "") + elem(lqn, g.first[i]);
                if (src.empty()) src = elem(lqn, eidx);
                blk.label = src + " -> " + elem(lqn, blk.join);
                blocks.push_back(blk);
                for (std::size_t h : g.second) placed.insert(h);
                for (std::size_t tl : blk.tails) placed.insert(tl);
            }
        }
    }
    for (std::size_t a : marked)
        if (!placed.count(a))
            bad.push_back("AND marking at activity " + elem(lqn, a) + " outside any well-formed fork/join block");
    return blocks;
}
 
/** python _scope_check: the blocks, or LINE:unsupported:NlpExport naming every offender. */
template <class T>
std::vector<ForkBlock> scope_check(const LqnStruct<T>& lqn, const std::string& fork) {
    std::vector<std::string> bad;
    for (std::size_t c = 1; c <= lqn.ncalls; ++c) {
        const CallType ct = c < lqn.calltype.size() ? lqn.calltype[c] : CallType::SYNC;
        if (ct != CallType::SYNC)
            bad.push_back(std::string(ct == CallType::ASYNC ? "asynchronous" : "forwarding") +
                          " call: " + callnm(lqn, c));
    }
    for (std::size_t i = 1; i <= lqn.nidx; ++i)
        if (i < lqn.has_arrival.size() && lqn.has_arrival[i])
            bad.push_back("open arrival at entry " + elem(lqn, i));
    for (std::size_t i = 1; i <= lqn.nidx; ++i) {
        const double r = i < lqn.repl.size() && !std::isnan(lqn.repl[i]) ? lqn.repl[i] : 1.0;
        if (r > 1) bad.push_back("replication of " + elem(lqn, i) + " (repl=" + fmt("%g", r) + ")");
    }
    for (std::size_t a = 1; a <= lqn.nacts; ++a)
        if (a < lqn.actphase.size() && lqn.actphase[a] > 1)
            bad.push_back("phase-2 activity " + elem(lqn, lqn.ashift + a));
    for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
        const std::size_t tidx = lqn.tshift + t;
        if (tidx < lqn.hassetup.size() && lqn.hassetup[tidx]) bad.push_back("setup time on task " + elem(lqn, tidx));
    }
    for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
        const std::size_t tidx = lqn.tshift + t;
        if (tidx < lqn.iscache.size() && lqn.iscache[tidx]) bad.push_back("cache task " + elem(lqn, tidx));
    }
    const std::vector<ForkBlock> blocks = fork_blocks(lqn, fork, bad);
    for (std::size_t i = 1; i < lqn.lincon_A.size(); ++i)
        if (lqn.lincon_A[i].rows() > 0) bad.push_back("admission region on " + elem(lqn, i));
    for (std::size_t a = 1; a <= lqn.nacts; ++a) {
        const std::size_t aidx = lqn.ashift + a;
        const std::size_t tidx = lqn.parent[aidx];
        std::vector<std::size_t> owners;
        if (tidx < lqn.entriesof.size())
            for (std::size_t e : lqn.entriesof[tidx])
                if (e < lqn.actsof.size() &&
                    std::find(lqn.actsof[e].begin(), lqn.actsof[e].end(), aidx) != lqn.actsof[e].end())
                    owners.push_back(e);
        if (owners.size() > 1) {
            std::string s;
            for (std::size_t i = 0; i < owners.size(); ++i) s += (i ? ", " : "") + elem(lqn, owners[i]);
            bad.push_back("activity " + elem(lqn, aidx) + " shared by entries " + s);
        } else if (owners.empty()) {
            bad.push_back("activity " + elem(lqn, aidx) + " bound to no entry");
        }
    }
    std::vector<std::size_t> refs;
    for (std::size_t t = 1; t <= lqn.ntasks; ++t)
        if (lqn.tshift + t < lqn.isref.size() && lqn.isref[lqn.tshift + t]) refs.push_back(lqn.tshift + t);
    if (refs.empty()) bad.push_back("no reference task: nothing drives the model");
    for (std::size_t tidx : refs) {
        const std::size_t ne = tidx < lqn.entriesof.size() ? lqn.entriesof[tidx].size() : 0;
        if (ne != 1)
            bad.push_back("reference task " + elem(lqn, tidx) + " has " + std::to_string(ne) +
                          " entries, expected 1");
    }
    if (!bad.empty()) {
        std::string msg = "cannot express this model as linear constraints plus a nonlinear objective:";
        for (const std::string& b : bad) msg += "\n  " + b;
        throw nlp_error(msg);
    }
    return blocks;
}
 
inline double tot(const std::vector<double>& pops) {
    double s = 0.0;
    for (double n : pops) {
        if (std::isinf(n)) return std::numeric_limits<double>::infinity();
        s += n;
    }
    return s;
}
 
inline double delta(double t) {
    if (std::isinf(t)) return 1.0;
    return t > 0 ? (t - 1.0) / t : 0.0;
}
 
inline std::string strip(const std::string& s) {
    const std::size_t b = s.find_first_not_of(" \t\n\r\f\v");
    if (b == std::string::npos) return "";
    const std::size_t e = s.find_last_not_of(" \t\n\r\f\v");
    return s.substr(b, e - b + 1);
}
 
}  // namespace nlp_detail
 
/**
 * Build the nonlinear program of a layered model: python export_nlp without the
 * emission, i.e. _collect + _attach_lld + the script's index sets, layout,
 * build_linear and build_residuals.
 *
 * @throws UnsupportedError LINE:unsupported:NlpExport naming every feature the
 *         form cannot express, a bad fork mode or alpha, or an unusable lld
 */
template <class T>
NlpProgram lqn_nlp_program(const LqnStruct<T>& lqn, const NlpOptions& opt = NlpOptions(),
                           const std::string& name = "") {
    using namespace nlp_detail;
    if (opt.fork != "refuse" && opt.fork != "exp")
        throw nlp_error("fork takes 'refuse' or 'exp', got '" + opt.fork + "'");
    const std::vector<ForkBlock> blocks = scope_check(lqn, opt.fork);
 
    NlpProgram p;
    p.name = name.empty() ? "lqn" : name;
    const double inf = std::numeric_limits<double>::infinity();
    auto todbl = [](const T& v) { return num_traits<T>::to_double(v); };
    auto mult = [&](std::size_t i) { return i < lqn.mult.size() ? lqn.mult[i] : std::nan(""); };
    auto isinf_elem = [&](std::size_t i) {
        return std::isinf(mult(i)) || (i < lqn.sched.size() && lqn.sched[i] == SchedStrategy::INF);
    };
    const std::size_t hs = lqn.hshift, ts = lqn.tshift, es = lqn.eshift, as = lqn.ashift;
    p.NH = lqn.nhosts;
    p.NT = lqn.ntasks;
    p.NE = lqn.nentries;
    p.NA = lqn.nacts;
    p.NC = lqn.ncalls;
    const std::vector<T> visits = detail::lqn_act_visits<T>(lqn);
 
    for (std::size_t h = 0; h < p.NH; ++h) {
        const std::size_t i = hs + h + 1;
        p.hosts.push_back(elem(lqn, i));
        p.host_mult.push_back(mult(i));
        const double mm = i < lqn.maxmult.size() ? lqn.maxmult[i] : std::nan("");
        p.host_maxmult.push_back((std::isnan(mm) || mm <= 0) ? mult(i) : mm);
        p.host_inf.push_back(isinf_elem(i));
    }
    for (std::size_t t = 0; t < p.NT; ++t) {
        const std::size_t i = ts + t + 1;
        p.tasks.push_back(elem(lqn, i));
        p.task_host.push_back(lqn.parent[i] - hs - 1);
        p.task_mult.push_back(mult(i));
        p.task_inf.push_back(isinf_elem(i));
        p.task_isref.push_back(i < lqn.isref.size() && lqn.isref[i]);
        p.task_think.push_back(todbl(detail::lqn_ref_thinktime(lqn, i)));
    }
    for (std::size_t e = 0; e < p.NE; ++e) {
        const std::size_t i = es + e + 1;
        p.entries.push_back(elem(lqn, i));
        p.entry_task.push_back(lqn.parent[i] - ts - 1);
    }
    for (std::size_t a = 0; a < p.NA; ++a) {
        const std::size_t i = as + a + 1;
        p.acts.push_back(elem(lqn, i));
        p.act_task.push_back(lqn.parent[i] - ts - 1);
        p.act_entry.push_back(detail::lqn_entry_of_activity(lqn, i) - es - 1);
        p.act_demand.push_back(hostdem_mean(lqn, i));
        p.act_visits.push_back(todbl(visits[i]));
    }
    for (std::size_t c = 0; c < p.NC; ++c) {
        p.calls.push_back(callnm(lqn, c + 1));
        p.call_act.push_back(lqn.callpair_src[c + 1] - as - 1);
        p.call_entry.push_back(lqn.callpair_dst[c + 1] - es - 1);
        const double y = c + 1 < lqn.callproc_mean.size() ? todbl(lqn.callproc_mean[c + 1]) : 1.0;
        if (!std::isfinite(y) || y < 0)
            throw nlp_error("call " + callnm(lqn, c + 1) + " has a non-finite mean count " + py_repr(y));
        p.call_mean.push_back(y);
        p.call_caller_task.push_back(p.act_task[p.call_act[c]]);
    }
    for (const ForkBlock& b : blocks) {
        p.fork_entry.push_back(b.entry - es - 1);
        std::vector<std::vector<std::size_t>> brs;
        for (const std::vector<std::size_t>& br : b.branches) {
            std::vector<std::size_t> v;
            for (std::size_t a : br) v.push_back(a - as - 1);
            brs.push_back(v);
        }
        p.fork_branches.push_back(brs);
        p.fork_label.push_back(b.label);
    }
 
    // python _entry_order: callees before callers; a cycle is a call-path deadlock
    {
        std::vector<std::set<std::size_t>> dep(p.NE);
        for (std::size_t c = 0; c < p.NC; ++c) dep[p.act_entry[p.call_act[c]]].insert(p.call_entry[c]);
        std::vector<bool> done(p.NE, false);
        while (p.entry_order.size() < p.NE) {
            bool progressed = false;
            for (std::size_t e = 0; e < p.NE; ++e) {
                if (done[e]) continue;
                bool ok = true;
                for (std::size_t d : dep[e])
                    if (!done[d]) ok = false;
                if (ok) {
                    done[e] = true;
                    p.entry_order.push_back(e);
                    progressed = true;
                }
            }
            if (!progressed) {
                std::string stuck;
                for (std::size_t e = 0; e < p.NE; ++e)
                    if (!done[e]) stuck += (stuck.empty() ? "" : ", ") + p.entries[e];
                throw nlp_error("cyclic call graph over entries: " + stuck);
            }
        }
    }
 
    // python _attach_lld
    if (!std::isfinite(opt.alpha) || opt.alpha <= 0)
        throw nlp_error("the softmin sharpness alpha must be finite and positive, got " + py_repr(opt.alpha));
    p.alpha = opt.alpha;
    p.host_lld.assign(p.NH, "1.0");
    p.task_lld.assign(p.NT, "1.0");
    std::vector<std::shared_ptr<const LldExpr>> host_expr(p.NH), task_expr(p.NT);
    {
        std::vector<std::string> bad;
        for (const auto& kv : opt.lld) {
            const std::string& key = kv.first;
            if (kv.second.tabulated) {
                bad.push_back(key + ": a tabulated lldscaling row is not smooth -- its interpolant has a "
                                    "corner at every knot. Pass an expression in n instead.");
                continue;
            }
            std::vector<std::pair<int, std::size_t>> where;
            for (std::size_t i = 0; i < p.NH; ++i)
                if (p.hosts[i] == key) where.push_back({0, i});
            for (std::size_t i = 0; i < p.NT; ++i)
                if (p.tasks[i] == key) where.push_back({1, i});
            if (where.empty()) {
                bad.push_back(key + ": no processor or task by that name");
                continue;
            }
            if (where.size() > 1) {
                bad.push_back(key + ": names both a processor and a task");
                continue;
            }
            std::shared_ptr<const LldExpr> ex;
            const std::string why = check_lld(kv.second.expr, ex);
            if (!why.empty()) {
                bad.push_back(key + ": " + why);
                continue;
            }
            if (where[0].first == 0) {
                p.host_lld[where[0].second] = kv.second.expr;
                host_expr[where[0].second] = ex;
            } else {
                p.task_lld[where[0].second] = kv.second.expr;
                task_expr[where[0].second] = ex;
            }
        }
        if (!bad.empty()) {
            std::string msg = "cannot use these queue-dependent rates:";
            for (const std::string& b : bad) msg += "\n  " + b;
            throw nlp_error(msg);
        }
    }
    for (std::size_t h = 0; h < p.NH; ++h) {
        const bool capped = strip(p.host_lld[h]) == "1.0";
        p.host_rate_capped.push_back(capped);
        NlpRate r;
        r.delay = p.host_inf[h] || std::isinf(p.host_mult[h]);
        r.mult = p.host_mult[h];
        r.alpha = p.alpha;
        if (!capped) r.lld = host_expr[h];
        p.host_rate.push_back(r);
    }
    for (std::size_t t = 0; t < p.NT; ++t) {
        NlpRate r;
        r.delay = p.task_inf[t] || std::isinf(p.task_mult[t]);
        r.mult = p.task_mult[t];
        r.alpha = p.alpha;
        if (strip(p.task_lld[t]) != "1.0") r.lld = task_expr[t];
        p.task_rate.push_back(r);
    }
 
    // ---- derived index sets ----------------------------------------------
    const std::size_t NH = p.NH, NT = p.NT, NE = p.NE, NA = p.NA, NC = p.NC;
    p.acts_of_host.assign(NH, {});
    for (std::size_t h = 0; h < NH; ++h)
        for (std::size_t a = 0; a < NA; ++a)
            if (p.task_host[p.act_task[a]] == h) p.acts_of_host[h].push_back(a);
    p.acts_of_entry.assign(NE, {});
    for (std::size_t e = 0; e < NE; ++e)
        for (std::size_t a = 0; a < NA; ++a)
            if (p.act_entry[a] == e) p.acts_of_entry[e].push_back(a);
    p.calls_of_act.assign(NA, {});
    for (std::size_t a = 0; a < NA; ++a)
        for (std::size_t c = 0; c < NC; ++c)
            if (p.call_act[c] == a) p.calls_of_act[a].push_back(c);
    p.calls_into_entry.assign(NE, {});
    for (std::size_t e = 0; e < NE; ++e)
        for (std::size_t c = 0; c < NC; ++c)
            if (p.call_entry[c] == e) p.calls_into_entry[e].push_back(c);
    p.calls_into_task.assign(NT, {});
    for (std::size_t t = 0; t < NT; ++t)
        for (std::size_t c = 0; c < NC; ++c)
            if (p.entry_task[p.call_entry[c]] == t) p.calls_into_task[t].push_back(c);
    p.entries_of_task.assign(NT, {});
    for (std::size_t t = 0; t < NT; ++t)
        for (std::size_t e = 0; e < NE; ++e)
            if (p.entry_task[e] == t) p.entries_of_task[t].push_back(e);
 
    p.NFORK = p.fork_entry.size();
    for (const auto& blk : p.fork_branches)
        for (const auto& br : blk)
            for (std::size_t a : br) p.fork_acts.insert(a);
    p.blocks_of_entry.assign(NE, {});
    for (std::size_t e = 0; e < NE; ++e)
        for (std::size_t b = 0; b < p.NFORK; ++b)
            if (p.fork_entry[b] == e) p.blocks_of_entry[e].push_back(b);
    p.seq_acts_of_entry.assign(NE, {});
    for (std::size_t e = 0; e < NE; ++e)
        for (std::size_t a : p.acts_of_entry[e])
            if (!p.fork_acts.count(a)) p.seq_acts_of_entry[e].push_back(a);
 
    for (std::size_t t = 0; t < NT; ++t)
        if (p.task_isref[t]) p.ref_tasks.push_back(t);
    p.NR = p.ref_tasks.size();
    for (std::size_t t : p.ref_tasks) p.ref_entry.push_back(p.entries_of_task[t][0]);
 
    for (std::size_t t = 0; t < NT; ++t) p.task_delay.push_back(p.task_inf[t] || std::isinf(p.task_mult[t]));
    for (std::size_t h = 0; h < NH; ++h) p.host_delay.push_back(p.host_inf[h] || std::isinf(p.host_mult[h]));
    p.task_pop = p.task_mult;
 
    p.sibling_acts.assign(NA, {});
    for (std::size_t a = 0; a < NA; ++a)
        for (std::size_t ap : p.acts_of_host[p.task_host[p.act_task[a]]])
            if (p.act_task[ap] == p.act_task[a]) p.sibling_acts[a].push_back(ap);
    p.sibling_calls.assign(NC, {});
    for (std::size_t c = 0; c < NC; ++c)
        for (std::size_t cp : p.calls_into_task[p.entry_task[p.call_entry[c]]])
            if (p.call_caller_task[cp] == p.call_caller_task[c]) p.sibling_calls[c].push_back(cp);
 
    for (std::size_t h = 0; h < NH; ++h)
        if (!p.host_delay[h] && !p.acts_of_host[h].empty()) p.qd_hosts.push_back(h);
    for (std::size_t t = 0; t < NT; ++t)
        if (!p.task_delay[t] && !p.calls_into_task[t].empty()) p.qd_tasks.push_back(t);
    for (std::size_t i = 0; i < p.qd_hosts.size(); ++i) p.gh_of[p.qd_hosts[i]] = i;
    for (std::size_t i = 0; i < p.qd_tasks.size(); ++i) p.gt_of[p.qd_tasks[i]] = i;
    p.NGH = p.qd_hosts.size();
    p.NGT = p.qd_tasks.size();
 
    for (std::size_t h = 0; h < NH; ++h) {
        std::set<std::size_t> ts2;
        for (std::size_t a : p.acts_of_host[h]) ts2.insert(p.act_task[a]);
        std::vector<double> pops;
        for (std::size_t t : ts2) pops.push_back(p.task_pop[t]);
        p.host_ntot.push_back(tot(pops));
        p.host_delta.push_back(delta(p.host_ntot.back()));
    }
    for (std::size_t t = 0; t < NT; ++t) {
        std::set<std::size_t> ts2;
        for (std::size_t c : p.calls_into_task[t]) ts2.insert(p.call_caller_task[c]);
        std::vector<double> pops;
        for (std::size_t q : ts2) pops.push_back(p.task_pop[q]);
        p.task_ntot.push_back(tot(pops));
        p.task_delta.push_back(delta(p.task_ntot.back()));
    }
 
    // ---- layout ----------------------------------------------------------
    const std::size_t sizes[NlpProgram::NBLOCKS] = {NE, NA, NC, NA, NE, NC, NA, NC, NC, p.NR, p.NGH, p.NGT,
                                                    p.NFORK ? NE : 0};
    std::size_t off = 0;
    for (int b = 0; b < NlpProgram::NBLOCKS; ++b) {
        p.offset[b] = off;
        p.blocksize[b] = sizes[b];
        off += sizes[b];
    }
    p.NVAR = off;
    const std::size_t NV = p.NVAR;
    using B = NlpProgram;
 
    // ---- linear constraints: python build_linear, row for row -------------
    std::vector<std::vector<double>> aeq, aub;
    std::vector<double> beq, bub;
    auto row = [&]() { return std::vector<double>(NV, 0.0); };
    for (std::size_t c = 0; c < NC; ++c) {
        auto r = row();
        r[p.col(B::XC, c)] = 1.0;
        r[p.col(B::XA, p.call_act[c])] -= p.call_mean[c];
        aeq.push_back(r);
        beq.push_back(0.0);
    }
    for (std::size_t e = 0; e < NE; ++e) {
        if (p.task_isref[p.entry_task[e]]) continue;
        auto r = row();
        r[p.col(B::XE, e)] = 1.0;
        for (std::size_t c : p.calls_into_entry[e]) r[p.col(B::XC, c)] -= 1.0;
        aeq.push_back(r);
        beq.push_back(0.0);
    }
    for (std::size_t a = 0; a < NA; ++a) {
        auto r = row();
        r[p.col(B::XA, a)] = 1.0;
        r[p.col(B::XE, p.act_entry[a])] -= p.act_visits[a];
        aeq.push_back(r);
        beq.push_back(0.0);
    }
    for (std::size_t h = 0; h < NH; ++h) {
        if (!p.host_delay[h]) continue;
        for (std::size_t a : p.acts_of_host[h]) {
            auto r = row();
            r[p.col(B::W, a)] = 1.0;
            aeq.push_back(r);
            beq.push_back(p.act_demand[a]);
        }
    }
    auto composition = [&](std::vector<double>& r, std::size_t e) {
        for (std::size_t a : p.acts_of_entry[e]) {
            r[p.col(B::W, a)] -= p.act_visits[a];
            for (std::size_t c : p.calls_of_act[a]) r[p.col(B::R, c)] -= p.act_visits[a] * p.call_mean[c];
        }
    };
    for (std::size_t e = 0; e < NE; ++e) {
        if (!p.blocks_of_entry[e].empty()) continue;
        auto r = row();
        r[p.col(B::S, e)] = 1.0;
        composition(r, e);
        aeq.push_back(r);
        beq.push_back(0.0);
    }
    if (p.NFORK)
        for (std::size_t e = 0; e < NE; ++e) {
            auto r = row();
            r[p.col(B::SA, e)] = 1.0;
            composition(r, e);
            aeq.push_back(r);
            beq.push_back(0.0);
        }
    for (std::size_t c = 0; c < NC; ++c)
        if (p.task_delay[p.entry_task[p.call_entry[c]]]) {
            auto r = row();
            r[p.col(B::R, c)] = 1.0;
            r[p.col(B::S, p.call_entry[c])] -= 1.0;
            aeq.push_back(r);
            beq.push_back(0.0);
        }
    for (std::size_t ri = 0; ri < p.NR; ++ri) {
        const std::size_t t = p.ref_tasks[ri];
        auto r = row();
        r[p.col(B::XE, p.ref_entry[ri])] = p.task_think[t];
        r[p.col(B::BR, ri)] = 1.0;
        aeq.push_back(r);
        beq.push_back(p.task_mult[t]);
    }
    for (std::size_t h = 0; h < NH; ++h) {
        if (p.host_delay[h] || !p.host_rate_capped[h]) continue;
        auto r = row();
        for (std::size_t a : p.acts_of_host[h]) r[p.col(B::XA, a)] += p.act_demand[a];
        aub.push_back(r);
        bub.push_back(p.host_mult[h]);
    }
    for (std::size_t t = 0; t < NT; ++t) {
        if (p.task_isref[t] || p.task_delay[t] || p.calls_into_task[t].empty()) continue;
        auto r = row();
        for (std::size_t c : p.calls_into_task[t]) r[p.col(B::BC, c)] = 1.0;
        aub.push_back(r);
        bub.push_back(p.task_mult[t]);
    }
    for (std::size_t h = 0; h < NH; ++h) {
        if (p.host_delay[h] || std::isinf(p.host_ntot[h]) || p.acts_of_host[h].empty()) continue;
        auto r = row();
        for (std::size_t a : p.acts_of_host[h]) r[p.col(B::Q, a)] = 1.0;
        aub.push_back(r);
        bub.push_back(p.host_ntot[h]);
    }
    for (std::size_t t = 0; t < NT; ++t) {
        if (p.task_delay[t] || std::isinf(p.task_ntot[t]) || p.calls_into_task[t].empty()) continue;
        auto r = row();
        for (std::size_t c : p.calls_into_task[t]) r[p.col(B::QT, c)] = 1.0;
        aub.push_back(r);
        bub.push_back(p.task_ntot[t]);
    }
    auto tomat = [&](const std::vector<std::vector<double>>& rows) {
        Matrix<double> M(rows.size(), NV, 0.0);
        for (std::size_t i = 0; i < rows.size(); ++i)
            for (std::size_t j = 0; j < NV; ++j) M(i, j) = rows[i][j];
        return M;
    };
    p.AEQ = tomat(aeq);
    p.AUB = tomat(aub);
    p.BEQ = beq;
    p.BUB = bub;
 
    // ---- residuals: python build_residuals, in its order ------------------
    std::vector<std::vector<double>> lin;
    std::size_t k = 0;
    for (std::size_t a = 0; a < NA; ++a) {
        auto r = row();
        r[p.col(B::Q, a)] = 1.0;
        lin.push_back(r);
        p.res_const.push_back(0.0);
        p.bil.push_back({k++, p.col(B::XA, a), p.col(B::W, a), -1.0});
    }
    for (std::size_t c = 0; c < NC; ++c) {
        auto r = row();
        r[p.col(B::QT, c)] = 1.0;
        lin.push_back(r);
        p.res_const.push_back(0.0);
        p.bil.push_back({k++, p.col(B::XC, c), p.col(B::R, c), -1.0});
    }
    for (std::size_t c = 0; c < NC; ++c) {
        auto r = row();
        r[p.col(B::BC, c)] = 1.0;
        lin.push_back(r);
        p.res_const.push_back(0.0);
        p.bil.push_back({k++, p.col(B::XC, c), p.sent(p.call_entry[c]), -1.0});
    }
    for (std::size_t ri = 0; ri < p.NR; ++ri) {
        const std::size_t e = p.ref_entry[ri];
        auto r = row();
        r[p.col(B::BR, ri)] = 1.0;
        lin.push_back(r);
        p.res_const.push_back(0.0);
        p.bil.push_back({k++, p.col(B::XE, e), p.col(B::S, e), -1.0});
    }
    for (std::size_t e = 0; e < NE; ++e) {
        if (p.blocks_of_entry[e].empty()) continue;
        auto r = row();
        r[p.col(B::S, e)] = 1.0;
        for (std::size_t a : p.seq_acts_of_entry[e]) {
            r[p.col(B::W, a)] -= p.act_visits[a];
            for (std::size_t c : p.calls_of_act[a]) r[p.col(B::R, c)] -= p.act_visits[a] * p.call_mean[c];
        }
        lin.push_back(r);
        p.res_const.push_back(0.0);
        for (std::size_t b : p.blocks_of_entry[e]) {
            NlpProgram::Mnl t;
            t.k = k;
            t.s = -1.0;
            t.c.assign(p.fork_branches[b].size(), 0.0);
            t.a = Matrix<double>(p.fork_branches[b].size(), NV, 0.0);
            for (std::size_t i = 0; i < p.fork_branches[b].size(); ++i)
                for (std::size_t a : p.fork_branches[b][i]) {
                    t.a(i, p.col(B::W, a)) += p.act_visits[a];
                    for (std::size_t c : p.calls_of_act[a]) t.a(i, p.col(B::R, c)) += p.act_visits[a] * p.call_mean[c];
                }
            p.mnl.push_back(t);
        }
        ++k;
    }
    for (std::size_t h = 0; h < NH; ++h) {
        if (p.host_delay[h]) continue;
        const std::size_t g = p.col(B::GH, p.gh_of.at(h));
        for (std::size_t a : p.acts_of_host[h]) {
            auto r = row();
            for (const auto& jw : p.qbar_host_terms(h, a, -p.act_demand[a])) r[jw.first] += jw.second;
            lin.push_back(r);
            p.res_const.push_back(-p.act_demand[a]);
            p.bil.push_back({k, p.col(B::W, a), g, 1.0});
            ++k;
        }
    }
    for (std::size_t c = 0; c < NC; ++c) {
        const std::size_t e = p.call_entry[c], t = p.entry_task[e];
        if (p.task_delay[t]) continue;
        auto r = row();
        r[p.col(B::S, e)] = -1.0;
        lin.push_back(r);
        p.res_const.push_back(0.0);
        p.bil.push_back({k, p.col(B::R, c), p.col(B::GT, p.gt_of.at(t)), 1.0});
        for (const auto& jw : p.qbar_task_terms(t, c, -1.0)) p.bil.push_back({k, p.col(B::S, e), jw.first, jw.second});
        ++k;
    }
    for (std::size_t h : p.qd_hosts) {
        auto r = row();
        r[p.col(B::GH, p.gh_of.at(h))] = 1.0;
        lin.push_back(r);
        p.res_const.push_back(0.0);
        NlpProgram::Nlt t;
        t.k = k++;
        t.f = &p.host_rate[h];
        t.s = -1.0;
        t.c0 = 1.0;
        t.a.assign(NV, 0.0);
        for (std::size_t a : p.acts_of_host[h]) t.a[p.col(B::Q, a)] += p.host_delta[h];
        p.nlt.push_back(t);
    }
    for (std::size_t tq : p.qd_tasks) {
        auto r = row();
        r[p.col(B::GT, p.gt_of.at(tq))] = 1.0;
        lin.push_back(r);
        p.res_const.push_back(0.0);
        NlpProgram::Nlt t;
        t.k = k++;
        t.f = &p.task_rate[tq];
        t.s = -1.0;
        t.c0 = 1.0;
        t.a.assign(NV, 0.0);
        for (std::size_t c : p.calls_into_task[tq]) t.a[p.col(B::QT, c)] += p.task_delta[tq];
        p.nlt.push_back(t);
    }
    p.res_lin = tomat(lin);
    p.NRES = lin.size();
    return p;
}
 
}  // namespace lqn
}  // namespace line
 
#endif  // LINE_API_LQN_EXPORT_NLP_H