fluid_ddpp.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_FLUID_FLUID_DDPP_H
#define LINE_SOLVERS_FLUID_FLUID_DDPP_H
 
/**
 * @file
 * @ingroup line_solvers
 * Refined mean-field (O(1/N)) steady state of the matrix fluid drift,
 * `method = "ddpp"` (alias `fld.ddpp`). Port of `solver_fluid_ddpp.m` and
 * `fluid_ddpp_admits.m`.
 *
 * The model is read as a density-dependent population process whose size N is
 * the total closed population: the state is the job count of every (station,
 * class, phase) coordinate, a jump moving a job from a to b has rate
 * theta_a(x) * W(a,b), with theta_a = x_a * min(1, c_i / X_i) (x_a at an
 * infinite server). The family scales every population and finite server count
 * together, so the fixed point x* and the correction V are N-independent in job
 * counts. Gast and Van Houdt (POMACS 1(2), 2017) give E[X] = x* + V + O(1/N),
 * with J W + W J' + Q = 0 and V = -J^{-1} (1/2) H:W, both solved on an
 * orthonormal basis of the jump span (Allmeier and Gast, PER 49(4), 2022). The
 * rates are refined to the same order, so the corrected throughputs satisfy flow
 * balance exactly.
 *
 * The correction is nonzero only at a saturated MULTICLASS PS/FCFS share: a
 * single-class network gives V = 0 (piecewise-linear drift for exponential
 * service, vanishing Hessian contraction for phase-type service, so the result
 * is insensitive beyond the mean). The fixed point is found by the same
 * pseudo-transient continuation as the MATLAB, JAR and native Python ports.
 */
 
#include <algorithm>
#include <chrono>
#include <cmath>
#include <complex>
#include <cstddef>
#include <limits>
#include <string>
#include <vector>
 
#include "line/lang/qn/feature_set.h"
#include "line/lang/qn/network_struct.h"
#include "line/solvers/fluid/fluid_moments.h"
#include "line/solvers/fluid/fluid_nonhyperbolic.h"
#include "line/solvers/fluid/fluid_odes.h"
#include "line/solvers/fluid/solver_fluid.h"
#include "line/util/eig.h"
#include "line/util/error.h"
#include "line/util/linalg.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
#include "line/util/sylvester.h"
 
namespace line {
namespace fluid {
 
/** The verdict of `fluid_ddpp_admits`: ok, or the refusal and its identifier. */
struct FluidDdppAdmits {
    bool ok = false;
    std::string category;  ///< "unsupported" or "numerical" when refused
    std::string kind;      ///< the LINE:<category>:<kind> kind
    std::string reason;
 
    /** Throws the refusal as the error class its category names. */
    [[noreturn]] void raise() const {
        if (category == "numerical") throw NumericError(ErrorKind{kind}, "solver_fluid_ddpp: " + reason);
        throw UnsupportedError(ErrorKind{kind}, "solver_fluid_ddpp: " + reason);
    }
};
 
namespace detail {
 
template <class T>
bool fluid_ddpp_active(const qn::NetworkStruct<T>& sn, std::size_t i, std::size_t r) {
    if (sn.stations[i].sched == lang::SchedStrategy::EXT) return false;
    return fluid_service_defined(sn, i, r);
}
 
inline std::vector<double> fluid_ddpp_theta(const std::vector<double>& x,
                                            const std::vector<std::vector<std::size_t>>& members,
                                            const std::vector<double>& c) {
    std::vector<double> th = x;
    for (std::size_t i = 0; i < members.size(); ++i) {
        if (members[i].empty() || !std::isfinite(c[i])) continue;
        double X = 0.0;
        for (std::size_t a : members[i]) X += x[a];
        if (X > c[i])
            for (std::size_t a : members[i]) th[a] = c[i] * x[a] / X;
    }
    return th;
}
 
inline Matrix<double> fluid_ddpp_dtheta(const std::vector<double>& x,
                                        const std::vector<std::vector<std::size_t>>& members,
                                        const std::vector<double>& c) {
    Matrix<double> D = eye<double>(x.size());
    for (std::size_t i = 0; i < members.size(); ++i) {
        if (members[i].empty() || !std::isfinite(c[i])) continue;
        double X = 0.0;
        for (std::size_t a : members[i]) X += x[a];
        if (X > c[i])
            for (std::size_t a : members[i])
                for (std::size_t b : members[i])
                    D(a, b) = c[i] * ((a == b ? 1.0 / X : 0.0) - x[a] / (X * X));
    }
    return D;
}
 
/** (1/2) sum_{jk} d2theta_a/dx_j dx_k Sigma_jk, nonzero at saturated stations only. */
inline std::vector<double> fluid_ddpp_hesscontract(const std::vector<double>& x,
                                                   const std::vector<std::vector<std::size_t>>& members,
                                                   const std::vector<double>& c,
                                                   const Matrix<double>& Sigma) {
    std::vector<double> bh(x.size(), 0.0);
    for (std::size_t i = 0; i < members.size(); ++i) {
        if (members[i].empty() || !std::isfinite(c[i])) continue;
        double X = 0.0;
        for (std::size_t a : members[i]) X += x[a];
        if (!(X > c[i])) continue;
        std::vector<double> row(members[i].size(), 0.0);
        double tot = 0.0;
        for (std::size_t p = 0; p < members[i].size(); ++p) {
            for (std::size_t b : members[i]) row[p] += Sigma(members[i][p], b);
            tot += row[p];
        }
        for (std::size_t p = 0; p < members[i].size(); ++p) {
            const std::size_t a = members[i][p];
            bh[a] = c[i] * (-row[p] / (X * X) + x[a] * tot / (X * X * X));
        }
    }
    return bh;
}
 
/** Largest step in (0,1] that keeps the state nonnegative, backed off by 0.9. */
inline double fluid_ddpp_steplimit(const std::vector<double>& x, const std::vector<double>& dx) {
    double m = std::numeric_limits<double>::infinity();
    for (std::size_t a = 0; a < x.size(); ++a)
        if (dx[a] < 0.0 && x[a] + dx[a] < 0.0) m = std::min(m, -x[a] / dx[a]);
    return std::isfinite(m) ? std::min(1.0, 0.9 * m) : 1.0;
}
 
inline std::vector<double> fluid_ddpp_solve(const Matrix<double>& A, std::vector<double> b) {
    Matrix<double> LU = A;
    const std::vector<std::size_t> piv = lu_factor(LU);
    lu_solve(LU, piv, b);
    return b;
}
 
inline double fluid_ddpp_norminf(const std::vector<double>& v) {
    double m = 0.0;
    for (double z : v) m = std::max(m, std::fabs(z));
    return m;
}
 
inline double fluid_ddpp_norm2(const std::vector<double>& v) {
    double s = 0.0;
    for (double z : v) s += z * z;
    return std::sqrt(s);
}
 
/** 1-norm reciprocal condition estimate by explicit inverse; 0 when singular. */
inline double fluid_ddpp_rcond(const Matrix<double>& A) {
    const auto norm1 = [](const Matrix<double>& B) {
        double m = 0.0;
        for (std::size_t j = 0; j < B.cols(); ++j) {
            double s = 0.0;
            for (std::size_t i = 0; i < B.rows(); ++i) s += std::fabs(B(i, j));
            m = std::max(m, s);
        }
        return m;
    };
    try {
        const Matrix<double> Ai = inverse(A);
        const double na = norm1(A), ni = norm1(Ai);
        if (!std::isfinite(ni) || na == 0.0) return 0.0;
        return 1.0 / (na * ni);
    } catch (const Error&) {
        return 0.0;
    }
}
 
}  // namespace detail
 
/**
 * Port of `fluid_ddpp_admits.m`: the model class of 'ddpp', decided before the
 * fixed point is solved. Closed classes only, INF/PS/FCFS stations (FCFS with
 * exponential service at one rate for every class), renewal PH service, steady
 * state, no load dependence, and the moment_maxstate cap.
 */
template <class T>
FluidDdppAdmits fluid_ddpp_admits(const qn::NetworkStruct<T>& sn, const FluidOptions& opt) {
    FluidDdppAdmits v;
    const auto refuse = [&](const char* cat, const std::string& kind, const std::string& why) {
        v.ok = false;
        v.category = cat;
        v.kind = kind;
        v.reason = why;
        return v;
    };
    for (const auto& nd : sn.nodes)
        if (nd.nodetype == qn::NodeType::Transition || nd.nodetype == qn::NodeType::Place)
            return refuse("unsupported", "Transition",
                          "the 'ddpp' method has no drift for a stochastic Petri net. Use method 'dae'");
    for (const auto& nd : sn.nodes)
        if (nd.nodetype == qn::NodeType::Cache)
            return refuse("unsupported", "Cache",
                          "the 'ddpp' method does not support caching stations. Use method 'rmf'");
    for (const auto& nd : sn.nodes)
        if (nd.nodetype == qn::NodeType::Fork || nd.nodetype == qn::NodeType::Join)
            return refuse("unsupported", "Fork", "the 'ddpp' method does not support fork-join models");
    if (sn.has_open_classes())
        return refuse("unsupported", "OpenClass",
                      "the 'ddpp' method supports closed models only: its system size N is the total "
                      "closed population. Use method 'minnormal' for open or mixed models");
    if (std::isfinite(opt.timespan_end))
        return refuse("unsupported", "Analysis",
                      "the 'ddpp' method corrects the STATIONARY mean only, so it has no transient. "
                      "Leave timespan_end infinite, or use method 'matrix' for a trajectory");
    if (detail::fluid_has_time_varying_rates(opt))
        return refuse("unsupported", "Analysis",
                      "the 'ddpp' method needs an autonomous drift, but a rate schedule makes the "
                      "rates time-varying");
    const std::size_t M = sn.nstations, K = sn.nclasses;
    for (std::size_t i = 0; i < M; ++i)
        for (const T& s : sn.stations[i].lldscaling)
            if (std::fabs(num_traits<T>::to_double(s) - 1.0) > 1e-14)
                return refuse("unsupported", "LoadDependence",
                              "the 'ddpp' method does not evaluate load dependence. Use method 'closing'");
    std::size_t nstate = 0;
    for (std::size_t i = 0; i < M; ++i) {
        const lang::SchedStrategy sched = sn.stations[i].sched;
        if (!(sched == lang::SchedStrategy::INF || sched == lang::SchedStrategy::PS ||
              sched == lang::SchedStrategy::FCFS))
            return refuse("unsupported", qn::feature_name(qn::feature_of_sched(sched)),
                          "the 'ddpp' method supports INF, PS and FCFS stations only, and station " +
                              std::to_string(i + 1) + " uses " + lang::sched_to_text(sched));
        double fcfsRate = std::numeric_limits<double>::quiet_NaN();
        for (std::size_t r = 0; r < K; ++r) {
            if (!detail::fluid_service_defined(sn, i, r)) continue;
            const lang::Distrib<T>& d = sn.service[i][r];
            const lang::ProcessType pt = d.type;
            if (pt == lang::ProcessType::MAP || pt == lang::ProcessType::MMPP2 ||
                pt == lang::ProcessType::NHPP || pt == lang::ProcessType::MAPT ||
                pt == lang::ProcessType::PHT)
                return refuse("unsupported", "Distribution",
                              "the 'ddpp' method needs a renewal phase-type service, and class " +
                                  std::to_string(r + 1) + " at station " + std::to_string(i + 1) +
                                  " has a " + lang::process_to_text(pt) + " process");
            const std::size_t nph = d.D0.rows();
            nstate += nph;
            if (sched == lang::SchedStrategy::FCFS) {
                if (nph > 1)
                    return refuse("unsupported", "SchedStrategy_FCFS",
                                  "the 'ddpp' method takes an FCFS station only with exponential service, "
                                  "where its stationary counts are those of PS, and class " +
                                      std::to_string(r + 1) + " at station " + std::to_string(i + 1) +
                                      " has a " + std::to_string(nph) + "-phase service");
                const double rate = -num_traits<T>::to_double(d.D0(0, 0));
                if (std::isnan(fcfsRate)) {
                    fcfsRate = rate;
                } else if (std::fabs(rate - fcfsRate) > 1e-10 * std::max(1.0, fcfsRate)) {
                    return refuse("unsupported", "SchedStrategy_FCFS",
                                  "the 'ddpp' method takes an FCFS station only with one service rate for "
                                  "every class, where its stationary counts are those of PS, and station " +
                                      std::to_string(i + 1) + " serves at rates " +
                                      std::to_string(fcfsRate) + " and " + std::to_string(rate));
                }
            }
        }
    }
    const std::size_t maxstate = opt.moment_maxstate > 0 ? opt.moment_maxstate : 200;
    if (nstate > maxstate)
        return refuse("numerical", "StateSpace",
                      "the phase-resolved state has " + std::to_string(nstate) +
                          " coordinates, above the moment_maxstate limit of " + std::to_string(maxstate) +
                          " that bounds the Lyapunov solve");
    v.ok = true;
    return v;
}
 
/** What the expansion found beside the refined table: the fluid point and V. */
struct FluidDdppReport {
    std::vector<double> x;  ///< fluid fixed point x*
    std::vector<double> V;  ///< the O(1/N) correction, in job counts
    Matrix<double> QNfluid, UNfluid, TNfluid;
    double N = 0.0;         ///< total closed population
    double residual = 0.0;
};
 
/**
 * Port of `solver_fluid_ddpp.m`. Returns the refined QN, UN, RN, TN (M x K),
 * with `xvec` = x* + V. `report`, when given, receives the fluid point and V.
 * With `probe` it stops once the fixed point has passed the kink and stability
 * checks and returns an empty solution (see `fluid_ddpp_hyperbolic`).
 */
template <class T>
FluidSolution solver_fluid_ddpp(const qn::NetworkStruct<T>& sn, const FluidOptions& opt,
                                FluidDdppReport* report = nullptr, bool probe = false) {
    const FluidDdppAdmits adm = fluid_ddpp_admits(sn, opt);
    if (!adm.ok) adm.raise();
    const std::size_t M = sn.nstations, K = sn.nclasses;
    const double eps = std::numeric_limits<double>::epsilon();
 
    // state encoding: one coordinate per (station, class, phase)
    std::vector<std::vector<std::size_t>> start(M, std::vector<std::size_t>(K, 0));
    std::vector<std::vector<std::size_t>> len(M, std::vector<std::size_t>(K, 0));
    std::size_t n = 0;
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < K; ++r) {
            start[i][r] = n;
            if (detail::fluid_ddpp_active(sn, i, r)) {
                len[i][r] = sn.service[i][r].D0.rows();
                n += len[i][r];
            }
        }
    std::vector<std::size_t> st(n), cl(n);
    std::vector<double> exitRate(n, 0.0);
    std::vector<std::vector<std::vector<double>>> pie(M, std::vector<std::vector<double>>(K));
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < K; ++r) {
            if (len[i][r] == 0) continue;
            const lang::Distrib<T>& d = sn.service[i][r];
            pie[i][r] = detail::fluid_pie(d);
            for (std::size_t k = 0; k < len[i][r]; ++k) {
                const std::size_t a = start[i][r] + k;
                st[a] = i;
                cl[a] = r;
                double s = 0.0;
                for (std::size_t kk = 0; kk < d.D1.cols(); ++kk) s += num_traits<T>::to_double(d.D1(k, kk));
                exitRate[a] = s;
            }
        }
 
    // jump rates per unit theta, self-loops dropped; station routing by stochastic complement
    const std::size_t S = sn.nof_stateful();
    const bool have_rt = sn.rt.rows() == S * K;
    Matrix<double> P(M * K, M * K, 0.0);
    if (have_rt) {
        std::vector<std::size_t> keep;
        keep.reserve(M * K);
        for (std::size_t i = 0; i < M; ++i) {
            const std::size_t isf = sn.stateful_of_station(i + 1) - 1;
            for (std::size_t r = 0; r < K; ++r) keep.push_back(isf * K + r);
        }
        Matrix<double> rt_full(S * K, S * K, 0.0);
        for (std::size_t a = 0; a < S * K; ++a)
            for (std::size_t b = 0; b < S * K; ++b) rt_full(a, b) = num_traits<T>::to_double(sn.rt(a, b));
        P = mc::dtmc_stochcomp(rt_full, keep);
    }
    Matrix<double> Woff(n, n, 0.0);
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < K; ++r) {
            if (len[i][r] == 0) continue;
            const lang::Distrib<T>& d = sn.service[i][r];
            for (std::size_t k = 0; k < len[i][r]; ++k) {
                const std::size_t a = start[i][r] + k;
                for (std::size_t kk = 0; kk < len[i][r]; ++kk) {
                    const double v = num_traits<T>::to_double(d.D0(k, kk));
                    if (kk != k && v > 0.0) Woff(a, start[i][r] + kk) += v;
                }
                if (exitRate[a] <= 0.0) continue;
                for (std::size_t j = 0; j < M; ++j)
                    for (std::size_t s = 0; s < K; ++s) {
                        const double p = P(i * K + r, j * K + s);
                        if (!(p > 0.0) || len[j][s] == 0) continue;
                        for (std::size_t kk = 0; kk < len[j][s]; ++kk) {
                            const std::size_t b = start[j][s] + kk;
                            const double pj = kk < pie[j][s].size() ? pie[j][s][kk] : 0.0;
                            if (b != a && pj > 0.0) Woff(a, b) += exitRate[a] * p * pj;
                        }
                    }
            }
        }
    Matrix<double> WdT(n, n, 0.0);  // transpose of Wd = Woff - diag(rowsum)
    double rateScale = 1e-12;
    for (std::size_t a = 0; a < n; ++a) {
        double rs = 0.0;
        for (std::size_t b = 0; b < n; ++b) {
            WdT(b, a) = Woff(a, b);
            rs += Woff(a, b);
        }
        WdT(a, a) -= rs;
        for (std::size_t b = 0; b < n; ++b) rateScale = std::max(rateScale, std::fabs(WdT(b, a)));
    }
 
    std::vector<double> c(M);
    std::vector<std::vector<std::size_t>> members(M);
    for (std::size_t i = 0; i < M; ++i) c[i] = sn.stations[i].nservers;
    for (std::size_t a = 0; a < n; ++a) members[st[a]].push_back(a);
 
    // reduced basis: orthonormal span of the jump directions (column-major find order)
    std::size_t nj = 0;
    for (std::size_t b = 0; b < n; ++b)
        for (std::size_t a = 0; a < n; ++a)
            if (Woff(a, b) > 0.0) ++nj;
    Matrix<double> Djump(n, nj, 0.0);
    {
        std::size_t e = 0;
        for (std::size_t b = 0; b < n; ++b)
            for (std::size_t a = 0; a < n; ++a)
                if (Woff(a, b) > 0.0) {
                    Djump(a, e) = -1.0;
                    Djump(b, e) += 1.0;
                    ++e;
                }
    }
    const Matrix<double> Pb = detail::fluid_orth(Djump);
    const std::size_t d = Pb.cols();
    const Matrix<double> PbT = Pb.transpose();
 
    // initial state: each class population on its reference station, entry phases
    std::vector<double> x(n, 0.0);
    double Ntot = 0.0;
    const std::vector<double> njobs = sn.njobs();
    for (std::size_t r = 0; r < K; ++r) {
        const double Nr = njobs[r];
        if (!(std::isfinite(Nr) && Nr > 0.0)) continue;
        const std::size_t ref = sn.classes[r].refstat - 1;
        if (ref >= M || len[ref][r] == 0)
            throw InputError(ErrorKind{"Model"},
                             "solver_fluid_ddpp: class " + std::to_string(r + 1) +
                                 " has no service process at its reference station, so its population "
                                 "has no coordinate to start from");
        for (std::size_t k = 0; k < len[ref][r]; ++k)
            x[start[ref][r] + k] += Nr * (k < pie[ref][r].size() ? pie[ref][r][k] : 0.0);
        Ntot += Nr;
    }
 
    // fixed point: pseudo-transient continuation, then Newton on the reduced system
    const std::size_t iterMax = std::max<std::size_t>(100, opt.iter_max);
    const double tolF = 1e-12 * rateScale * std::max(1.0, Ntot);
    double h = 1.0 / rateScale;
    const double hmax = 1e10 / rateScale;
    std::vector<double> F = mulvec(WdT, detail::fluid_ddpp_theta(x, members, c));
    double nF = detail::fluid_ddpp_norminf(F);
    std::size_t iters = 0;
    while (nF > tolF && iters < iterMax) {
        ++iters;
        Matrix<double> A = matmul(WdT, detail::fluid_ddpp_dtheta(x, members, c));
        for (std::size_t a = 0; a < n; ++a)
            for (std::size_t b = 0; b < n; ++b) A(a, b) = -A(a, b);
        for (std::size_t a = 0; a < n; ++a) A(a, a) += 1.0 / h;
        const std::vector<double> dx = detail::fluid_ddpp_solve(A, F);
        const double alpha = detail::fluid_ddpp_steplimit(x, dx);
        for (std::size_t a = 0; a < n; ++a) x[a] = std::max(0.0, x[a] + alpha * dx[a]);
        std::vector<double> Fnew = mulvec(WdT, detail::fluid_ddpp_theta(x, members, c));
        const double nFnew = detail::fluid_ddpp_norminf(Fnew);
        const double ratio = nF / std::max(nFnew, std::numeric_limits<double>::min());
        h = std::min(hmax, h * std::min(10.0, std::max(0.5, ratio)));
        F = Fnew;
        nF = nFnew;
    }
    for (int polish = 0; polish < 3 && d > 0; ++polish) {
        const Matrix<double> Jr = matmul(matmul(PbT, matmul(WdT, detail::fluid_ddpp_dtheta(x, members, c))), Pb);
        if (detail::fluid_ddpp_rcond(Jr) < eps) break;
        const std::vector<double> rdx = detail::fluid_ddpp_solve(Jr, mulvec(PbT, F));
        std::vector<double> dx = mulvec(Pb, rdx);
        for (double& z : dx) z = -z;
        const double alpha = detail::fluid_ddpp_steplimit(x, dx);
        for (std::size_t a = 0; a < n; ++a) x[a] = std::max(0.0, x[a] + alpha * dx[a]);
        F = mulvec(WdT, detail::fluid_ddpp_theta(x, members, c));
    }
    nF = detail::fluid_ddpp_norminf(F);
    if (!(nF <= 1e3 * tolF))
        throw NumericError(ErrorKind{"NonConvergence"},
                           "solver_fluid_ddpp: the fluid fixed point was not found in " +
                               std::to_string(iters) + " continuation steps (drift residual " +
                               std::to_string(nF) + " against a tolerance of " +
                               std::to_string(1e3 * tolF) + ")");
 
    // the expansion needs a smooth, hyperbolic fixed point
    for (std::size_t i = 0; i < M; ++i) {
        if (!std::isfinite(c[i]) || members[i].empty()) continue;
        double Xi = 0.0;
        for (std::size_t a : members[i]) Xi += x[a];
        if (std::fabs(Xi - c[i]) <= 1e-6 * std::max(1.0, c[i]))
            throw FluidNonHyperbolicError(
                "solver_fluid_ddpp: the fluid fixed point puts " + std::to_string(Xi) +
                " jobs at station " + std::to_string(i + 1) + ", which has " + std::to_string(c[i]) +
                " servers: the drift has a kink there, so the 1/N expansion has no Jacobian to expand "
                "around. Use method 'minnormal' or an exact solver");
    }
    const Matrix<double> Dth = detail::fluid_ddpp_dtheta(x, members, c);
    const Matrix<double> J = matmul(WdT, Dth);
    Matrix<double> Jr(0, 0, 0.0);
    if (d > 0) {
        Jr = matmul(matmul(PbT, J), Pb);
        double maxRe = -std::numeric_limits<double>::infinity();
        for (const std::complex<double>& z : eig_values(Jr)) maxRe = std::max(maxRe, z.real());
        if (maxRe >= -std::sqrt(eps) * rateScale)
            throw FluidNonHyperbolicError(
                "solver_fluid_ddpp: the fluid fixed point is not exponentially stable on the reachable "
                "subspace (largest Jacobian eigenvalue has real part " + std::to_string(maxRe) +
                "), so the stationary covariance and the 1/N correction do not exist");
    }
    if (probe) return FluidSolution();
 
    // noise covariance and its Lyapunov solve on the jump span
    const std::vector<double> th0 = detail::fluid_ddpp_theta(x, members, c);
    Matrix<double> Qn(n, n, 0.0);
    for (std::size_t a = 0; a < n; ++a)
        for (std::size_t b = 0; b < n; ++b) {
            const double rj = th0[a] * Woff(a, b);
            if (rj == 0.0) continue;
            Qn(a, a) += rj;
            Qn(b, b) += rj;
            Qn(a, b) -= rj;
            Qn(b, a) -= rj;
        }
    Matrix<double> Sigma(n, n, 0.0);
    if (d > 0) {
        Matrix<double> Qr = matmul(matmul(PbT, Qn), Pb);
        detail::fluid_symmetrize(Qr);
        Matrix<double> negQr(d, d, 0.0);
        for (std::size_t a = 0; a < d; ++a)
            for (std::size_t b = 0; b < d; ++b) negQr(a, b) = -Qr(a, b);
        Matrix<double> Wr = sylvester_schur(Jr, Jr.transpose(), negQr);
        detail::fluid_symmetrize(Wr);
        Sigma = matmul(matmul(Pb, Wr), PbT);
        detail::fluid_symmetrize(Sigma);
    }
 
    // Hessian contraction and the correction
    const std::vector<double> bh = detail::fluid_ddpp_hesscontract(x, members, c, Sigma);
    std::vector<double> V(n, 0.0);
    if (d > 0) {
        const std::vector<double> rv = detail::fluid_ddpp_solve(Jr, mulvec(PbT, mulvec(WdT, bh)));
        const std::vector<double> pv = mulvec(Pb, rv);
        for (std::size_t a = 0; a < n; ++a) V[a] = -pv[a];
    }
    const double normV = detail::fluid_ddpp_norm2(V), normX = detail::fluid_ddpp_norm2(x);
    if (normV > 0.5 * std::max(normX, std::sqrt(eps)))
        throw NumericError(ErrorKind{"Accuracy"},
                           "solver_fluid_ddpp: the 1/N correction (norm " + std::to_string(normV) +
                               ") is not small against the fluid fixed point (norm " + std::to_string(normX) +
                               "), so the expansion is outside its range of validity at this population");
    const std::vector<double> dthV = mulvec(Dth, V);
    std::vector<double> thR(n), xR(n);
    for (std::size_t a = 0; a < n; ++a) {
        thR[a] = th0[a] + dthV[a] + bh[a];
        xR[a] = x[a] + V[a];
    }
 
    // metrics
    const auto metrics = [&](const std::vector<double>& xs, const std::vector<double>& th,
                             Matrix<double>& QN, Matrix<double>& UN, Matrix<double>& TN) {
        QN = Matrix<double>(M, K, 0.0);
        UN = Matrix<double>(M, K, 0.0);
        TN = Matrix<double>(M, K, 0.0);
        for (std::size_t a = 0; a < n; ++a) {
            QN(st[a], cl[a]) += xs[a];
            TN(st[a], cl[a]) += exitRate[a] * th[a];
            if (std::isfinite(c[st[a]])) UN(st[a], cl[a]) += th[a] / c[st[a]];
        }
        for (std::size_t i = 0; i < M; ++i)
            if (sn.stations[i].sched == lang::SchedStrategy::INF)
                for (std::size_t r = 0; r < K; ++r) UN(i, r) = QN(i, r);
    };
    FluidSolution out;
    metrics(xR, thR, out.QN, out.UN, out.TN);
    out.RN = Matrix<double>(M, K, 0.0);
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < K; ++r)
            if (out.TN(i, r) > 0.0) out.RN(i, r) = out.QN(i, r) / out.TN(i, r);
    out.xvec = xR;
    out.iters = iters;
    out.method = "ddpp";
    out.XN.assign(K, 0.0);
    out.CN.assign(K, 0.0);
    for (std::size_t k = 0; k < K; ++k) {
        const std::size_t rs = sn.classes[k].refstat;
        if (rs >= 1 && rs <= M) out.XN[k] = out.TN(rs - 1, k);
        double q = 0.0;
        for (std::size_t i = 0; i < M; ++i) q += out.QN(i, k);
        if (out.XN[k] > 0.0) out.CN[k] = q / out.XN[k];
    }
    if (report) {
        report->x = x;
        report->V = V;
        metrics(x, th0, report->QNfluid, report->UNfluid, report->TNfluid);
        report->N = Ntot;
        report->residual = nF;
    }
    return out;
}
 
/**
 * Does the fluid fixed point admit the 1/N expansion? `fluid_ddpp_admits` is
 * static, but smoothness and exponential stability of the fixed point are a
 * property of the numbers: a closed Delay -> FCFS repairman with N = 3, think
 * rate 1 and service rate 2 puts exactly one job at the single server, on the
 * kink of min(1, c/X). The support gate asks this after the static predicate,
 * so findSolver does not offer a run that would refuse. Runs the solver up to
 * its hyperbolicity checks. Port of `fluid_ddpp_hyperbolic.m`.
 */
template <class T>
FluidDdppAdmits fluid_ddpp_hyperbolic(const qn::NetworkStruct<T>& sn, const FluidOptions& opt) {
    FluidDdppAdmits v;
    try {
        solver_fluid_ddpp(sn, opt, nullptr, true);
    } catch (const FluidNonHyperbolicError& e) {
        v.category = "numerical";
        v.kind = "FluidNonHyperbolic";
        v.reason = e.what();
        return v;
    } catch (const NumericError& e) {
        v.category = "numerical";
        v.kind = "NonConvergence";
        v.reason = e.what();
        return v;
    }
    v.ok = true;
    return v;
}
 
}  // namespace fluid
}  // namespace line
 
#endif  // LINE_SOLVERS_FLUID_FLUID_DDPP_H