spn_ispn.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_SPN_SPN_ISPN_H
#define LINE_API_SPN_SPN_ISPN_H
 
/**
 * @file
 * @ingroup api_spn
 * Infinite-state stochastic Petri nets (iSPN) by the matrix-geometric method.
 *
 * A net whose single unbounded place P0 is the level of a quasi-birth-death
 * process (B. R. Haverkort, "Performance of Computer Communication Systems: A
 * Model-Based Approach", Wiley, 1998, ch. 17; B. R. Haverkort, A. Ost, "Steady-
 * state analysis of infinite stochastic Petri nets: comparing the spectral
 * expansion and the matrix-geometric method", PNPM 1997, 36-45).
 *
 * The iSPN requirements of Haverkort sec. 17.2 are checked, and a violation is
 * refused with `LINE:unsupported:Structure` naming the condition: every firing
 * moves P0 by at most one token (also across a chain of immediate firings),
 * every other place is bounded, and above a threshold level K the levels repeat
 * with the same phase markings and level-independent rates. The threshold is
 * found by exploring a window of levels and comparing the generator blocks.
 *
 * Vanishing markings of immediate transitions are eliminated as
 * Q = Q_TT + Q_TV (I - P_VV)^-1 P_VT. R solves A0 + R A1 + R^2 A2 = 0 by
 * logarithmic reduction (G. Latouche, V. Ramaswami, J. Appl. Probab. 30(3),
 * 1993) or successive substitution, the drift condition pi A0 1 < pi A2 1
 * (Haverkort (17.7)) is checked, and the boundary levels 0..K solve x B = 0
 * with the R A2 correction and the normalisation x_K (I-R)^-1 1 = 1 - sum.
 *
 * NO TRUNCATION: a measure r of the marking that is affine in the level above
 * K, r(K+j) = a + (j-1) b for j >= 1, has the closed form
 *   E[r] = sum_{i<=K} x_i r_i + x_K R N a + x_K R^2 N^2 b,  N = (I-R)^-1,
 * which covers the mean marking of every place, the marginals of the bounded
 * places, the firing rates and the user rewards. A measure that is not affine
 * is refused rather than truncated.
 *
 * DOUBLE ONLY: R is a fixed point of a quadratic matrix equation, irrational in
 * general, and the stability test needs the spectral radius of R.
 *
 * Indices are 0-based: p0 is the row of the unbounded place in pre/post.
 * MATLAB twin: matlab/src/api/spn/spn_ispn.m (1-based P0); Java twin:
 * jline.api.spn.Spn_ispn; Python twin: line_solver.api.spn.ispn.
 */
 
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <functional>
#include <limits>
#include <map>
#include <string>
#include <vector>
 
#include "line/api/mam/qbd_r.h"
#include "line/util/eig.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace spn {
 
/** Options of spn_ispn. */
struct IspnOptions {
    int extra = 3;            ///< levels explored above the structural threshold (at least 3)
    int maxlevel = 200;       ///< largest explored window before the net is refused
    std::size_t maxstates = 200000;  ///< largest number of markings explored
    double tol = 1e-9;        ///< tolerance of the block-repetition and affinity tests
    double distmax = 1e5;     ///< largest level of the reported distribution of P0
    /** reward functions of the marking (length P), affine in the level above K */
    std::vector<std::function<double(const std::vector<double>&)>> rewards;
    std::string qbdmethod = "logred";  ///< "logred" or "successive"
};
 
/** Result of spn_ispn; fields as in the MATLAB help. */
struct IspnResult {
    int K = 0;                 ///< first repeating level
    std::size_t L = 0;         ///< phases per level
    std::vector<std::vector<double>> A0, A1, A2, R, B;
    std::vector<std::vector<double>> x;                       ///< stationary vectors of levels 0..K
    std::vector<std::vector<std::vector<double>>> levelmarkings;  ///< P x n markings of levels 0..K
    std::vector<double> drift;  ///< {up, down} mean drift above K
    double spR = 0;
    std::vector<double> meanmarking;
    double varlevel = 0;
    std::vector<double> leveldist;
    double leveltail = 0;
    std::vector<std::vector<double>> marginal;  ///< per place (P0: leveldist)
    std::vector<double> tput;
    std::vector<double> reward;
    std::size_t ntangible = 0;
    std::size_t nvanishing = 0;
};
 
namespace ispn_detail {
 
typedef std::vector<double> Vec;
typedef std::vector<Vec> Mat;
 
inline std::string fmt(const char* f, double a = 0, double b = 0, double c = 0, double d = 0) {
    char buf[512];
    std::snprintf(buf, sizeof buf, f, a, b, c, d);
    return std::string(buf);
}
 
[[noreturn]] inline void refuse(const std::string& m) {
    throw UnsupportedError(ErrorKind{"Structure"}, "spn_ispn: " + m);
}
[[noreturn]] inline void unstable(const std::string& m) {
    throw NumericError(ErrorKind{"Instability"}, "spn_ispn: " + m);
}
[[noreturn]] inline void badarg(const std::string& m) {
    throw InputError(ErrorKind{"Argument"}, "spn_ispn: " + m);
}
 
inline Mat zeros(std::size_t n, std::size_t m) { return Mat(n, Vec(m, 0.0)); }
inline Mat eye(std::size_t n) {
    Mat I = zeros(n, n);
    for (std::size_t i = 0; i < n; ++i) I[i][i] = 1.0;
    return I;
}
inline std::size_t ncols(const Mat& A) { return A.empty() ? 0 : A[0].size(); }
inline Mat sub(const Mat& A, const std::vector<std::size_t>& r, const std::vector<std::size_t>& c) {
    Mat B = zeros(r.size(), c.size());
    for (std::size_t i = 0; i < r.size(); ++i)
        for (std::size_t j = 0; j < c.size(); ++j) B[i][j] = A[r[i]][c[j]];
    return B;
}
inline bool all_zero(const Mat& A) {
    for (const Vec& r : A)
        for (double v : r)
            if (v != 0) return false;
    return true;
}
inline double max_abs(const Vec& a) {
    double m = 0;
    for (double v : a) m = std::max(m, std::fabs(v));
    return m;
}
inline double max_abs(const Mat& A) {
    double m = 0;
    for (const Vec& r : A) m = std::max(m, max_abs(r));
    return m;
}
inline double max_abs_diff(const Mat& A, const Mat& B) {
    double m = 0;
    for (std::size_t i = 0; i < A.size(); ++i)
        for (std::size_t j = 0; j < A[i].size(); ++j) m = std::max(m, std::fabs(A[i][j] - B[i][j]));
    return m;
}
inline Mat plus(const Mat& A, const Mat& B) {
    Mat C = A;
    for (std::size_t i = 0; i < A.size(); ++i)
        for (std::size_t j = 0; j < A[i].size(); ++j) C[i][j] += B[i][j];
    return C;
}
inline Mat minus(const Mat& A, const Mat& B) {
    Mat C = A;
    for (std::size_t i = 0; i < A.size(); ++i)
        for (std::size_t j = 0; j < A[i].size(); ++j) C[i][j] -= B[i][j];
    return C;
}
inline Mat scale(const Mat& A, double s) {
    Mat C = A;
    for (Vec& r : C)
        for (double& v : r) v *= s;
    return C;
}
inline Mat mult(const Mat& A, const Mat& B) {
    const std::size_t n = A.size(), m = B.size(), p = ncols(B);
    Mat C = zeros(n, p);
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t k = 0; k < m; ++k) {
            const double a = A[i][k];
            if (a == 0) continue;
            for (std::size_t j = 0; j < p; ++j) C[i][j] += a * B[k][j];
        }
    return C;
}
inline Vec vec_mult(const Vec& x, const Mat& A) {
    const std::size_t p = ncols(A);
    Vec y(p, 0.0);
    for (std::size_t k = 0; k < x.size(); ++k)
        for (std::size_t j = 0; j < p; ++j) y[j] += x[k] * A[k][j];
    return y;
}
inline double sum(const Vec& a) {
    double s = 0;
    for (double v : a) s += v;
    return s;
}
inline Vec row_sums(const Mat& A) {
    Vec s(A.size());
    for (std::size_t i = 0; i < A.size(); ++i) s[i] = sum(A[i]);
    return s;
}
inline double dot(const Vec& a, const Vec& b) {
    double s = 0;
    for (std::size_t i = 0; i < a.size(); ++i) s += a[i] * b[i];
    return s;
}
 
/** LU with partial pivoting in place; returns the row permutation. */
inline std::vector<std::size_t> lu(Mat& A) {
    const std::size_t n = A.size();
    std::vector<std::size_t> piv(n);
    for (std::size_t i = 0; i < n; ++i) piv[i] = i;
    for (std::size_t k = 0; k < n; ++k) {
        std::size_t p = k;
        for (std::size_t i = k + 1; i < n; ++i)
            if (std::fabs(A[i][k]) > std::fabs(A[p][k])) p = i;
        if (p != k) {
            std::swap(A[p], A[k]);
            std::swap(piv[p], piv[k]);
        }
        const double d = A[k][k];
        if (d == 0) continue;
        for (std::size_t i = k + 1; i < n; ++i) {
            const double f = A[i][k] / d;
            A[i][k] = f;
            if (f == 0) continue;
            for (std::size_t j = k + 1; j < n; ++j) A[i][j] -= f * A[k][j];
        }
    }
    return piv;
}
 
/** Solves A X = B. */
inline Mat solve(const Mat& A0, const Mat& B) {
    const std::size_t n = A0.size(), m = ncols(B);
    Mat A = A0;
    const std::vector<std::size_t> piv = lu(A);
    Mat X = zeros(n, m);
    for (std::size_t c = 0; c < m; ++c) {
        Vec y(n, 0.0);
        for (std::size_t i = 0; i < n; ++i) {
            double s = B[piv[i]][c];
            for (std::size_t j = 0; j < i; ++j) s -= A[i][j] * y[j];
            y[i] = s;
        }
        for (std::size_t ii = n; ii-- > 0;) {
            double s = y[ii];
            for (std::size_t j = ii + 1; j < n; ++j) s -= A[ii][j] * X[j][c];
            X[ii][c] = s / A[ii][ii];
        }
    }
    return X;
}
inline Mat inv(const Mat& A) { return solve(A, eye(A.size())); }
 
/** Solves x A = b for the row vector x. */
inline Vec solve_left(const Mat& A, const Vec& b) {
    const std::size_t n = A.size();
    Mat At = zeros(n, n), bb = zeros(n, 1);
    for (std::size_t i = 0; i < n; ++i) {
        bb[i][0] = b[i];
        for (std::size_t j = 0; j < n; ++j) At[j][i] = A[i][j];
    }
    const Mat X = solve(At, bb);
    Vec y(n);
    for (std::size_t i = 0; i < n; ++i) y[i] = X[i][0];
    return y;
}
inline double norm1(const Mat& A) {
    double m = 0;
    for (std::size_t j = 0; j < ncols(A); ++j) {
        double s = 0;
        for (std::size_t i = 0; i < A.size(); ++i) s += std::fabs(A[i][j]);
        m = std::max(m, s);
    }
    return m;
}
inline double cond1(const Mat& A) {
    const Mat Ai = inv(A);
    for (const Vec& r : Ai)
        for (double v : r)
            if (!std::isfinite(v)) return std::numeric_limits<double>::infinity();
    return norm1(A) * norm1(Ai);
}
inline Matrix<double> to_matrix(const Mat& A) {
    Matrix<double> M(A.size(), ncols(A));
    for (std::size_t i = 0; i < A.size(); ++i)
        for (std::size_t j = 0; j < A[i].size(); ++j) M(i, j) = A[i][j];
    return M;
}
inline Mat from_matrix(const Matrix<double>& M) {
    Mat A = zeros(M.rows(), M.cols());
    for (std::size_t i = 0; i < M.rows(); ++i)
        for (std::size_t j = 0; j < M.cols(); ++j) A[i][j] = M(i, j);
    return A;
}
 
struct Net {
    const Mat* pre;
    const Mat* post;
    const Mat* inh;
    const Vec* rate;
    const std::vector<bool>* imm;
    const Vec* srv;
    std::size_t p0, P, T;
    std::vector<std::size_t> other;
};
 
struct Arc {
    std::size_t src;  ///< 1-based source marking
    std::size_t dst;  ///< 1-based destination marking, 0 = escape above the window
    std::size_t trn;
    double val;
};
 
struct Walk {
    std::vector<Vec> S;
    std::vector<Arc> arcs;
    std::vector<bool> isvan;
};
 
struct Elim {
    Mat Q;                                       ///< nT x nT generator of the tangible markings
    Mat S;                                       ///< P x nT markings
    std::vector<std::vector<std::size_t>> levidx;  ///< tangible indices by level 0..Ktop+1
    Mat trate;                                   ///< nT x T firing rates
    std::size_t nvan = 0;
};
 
inline Vec enabled_weights(const Net& net, const Vec& m, std::vector<bool>& en) {
    Vec w(net.T, 0.0);
    for (std::size_t t = 0; t < net.T; ++t) {
        bool ok = true;
        for (std::size_t p = 0; p < net.P && ok; ++p)
            if (m[p] < (*net.pre)[p][t] || m[p] >= (*net.inh)[p][t]) ok = false;
        if (!ok) continue;
        double dg = std::numeric_limits<double>::infinity();
        for (std::size_t p = 0; p < net.P; ++p)
            if ((*net.pre)[p][t] > 0) dg = std::min(dg, std::floor(m[p] / (*net.pre)[p][t]));
        if (std::isinf(dg)) dg = 1;
        en[t] = true;
        w[t] = (*net.rate)[t] * std::min(dg, (*net.srv)[t]);
    }
    return w;
}
 
inline std::size_t lookup(Walk& W, std::map<std::vector<long long>, std::size_t>& keys, const Vec& m,
                          const Net& net, int cap, const IspnOptions& opt) {
    std::vector<long long> k(m.size());
    for (std::size_t i = 0; i < m.size(); ++i) k[i] = static_cast<long long>(m[i]);
    auto it = keys.find(k);
    if (it != keys.end()) return it->second;
    if (W.S.size() + 1 > opt.maxstates)
        refuse(fmt("iSPN requirement 1 violated: more than %.0f markings below level %.0f; a place other "
                   "than %.0f appears unbounded. Raise OPTIONS.maxstates if the net is just large",
                   static_cast<double>(opt.maxstates), cap, static_cast<double>(net.p0)));
    W.S.push_back(m);
    keys[k] = W.S.size();
    return W.S.size();
}
 
inline Walk explore(const Net& net, const Vec& m0, int cap, const IspnOptions& opt) {
    Walk W;
    std::map<std::vector<long long>, std::size_t> keys;
    if (m0[net.p0] > cap) badarg("the initial marking of P0 exceeds the explored window");
    lookup(W, keys, m0, net, cap, opt);
    std::size_t head = 0;
    while (head < W.S.size()) {
        const Vec m = W.S[head];
        ++head;
        std::vector<bool> en(net.T, false);
        Vec w = enabled_weights(net, m, en);
        bool v = false;
        for (std::size_t t = 0; t < net.T; ++t) v = v || (en[t] && (*net.imm)[t]);
        W.isvan.push_back(v);
        std::vector<std::size_t> fire;
        double tot = 0;
        for (std::size_t t = 0; t < net.T; ++t)
            if (en[t] && (*net.imm)[t] == v && w[t] > 0) {
                fire.push_back(t);
                tot += w[t];
            }
        if (v) {
            if (tot <= 0) badarg("an enabled immediate transition has zero weight");
            for (double& x : w) x /= tot;
        }
        for (std::size_t t : fire) {
            Vec mn(net.P);
            for (std::size_t p = 0; p < net.P; ++p) mn[p] = m[p] - (*net.pre)[p][t] + (*net.post)[p][t];
            const std::size_t id = mn[net.p0] > cap ? 0 : lookup(W, keys, mn, net, cap, opt);
            W.arcs.push_back(Arc{head, id, t, w[t]});
        }
    }
    return W;
}
 
inline Elim eliminate(const Walk& W, const Net& net, int Ktop) {
    const std::size_t n = W.S.size(), p0 = net.p0, T = net.T;
    std::vector<std::size_t> tgl, vnl;
    for (std::size_t i = 0; i < n; ++i) (W.isvan[i] ? vnl : tgl).push_back(i);
    const std::vector<std::size_t>& other = net.other;
    std::stable_sort(tgl.begin(), tgl.end(), [&](std::size_t a, std::size_t b) {
        if (W.S[a][p0] != W.S[b][p0]) return W.S[a][p0] < W.S[b][p0];
        for (std::size_t p : other)
            if (W.S[a][p] != W.S[b][p]) return W.S[a][p] < W.S[b][p];
        return false;
    });
    const std::size_t nT = tgl.size(), nV = vnl.size();
    std::vector<std::size_t> mapT(n, 0), mapV(n, 0);
    for (std::size_t i = 0; i < nT; ++i) mapT[tgl[i]] = i;
    for (std::size_t i = 0; i < nV; ++i) mapV[vnl[i]] = i;
    const std::size_t esc = nT;
    Mat QTT = zeros(nT, nT + 1), QTV = zeros(nT, nV), PVT = zeros(nV, nT + 1), PVV = zeros(nV, nV);
    Mat trate = zeros(nT, T), Fimm = zeros(nV, T);
    for (const Arc& arc : W.arcs) {
        const std::size_t s0 = arc.src - 1, d = arc.dst, t = arc.trn;
        const double v = arc.val;
        if (!W.isvan[s0]) {
            const std::size_t i = mapT[s0];
            trate[i][t] += v;
            if (d == 0) QTT[i][esc] += v;
            else if (W.isvan[d - 1]) QTV[i][mapV[d - 1]] += v;
            else QTT[i][mapT[d - 1]] += v;
        } else {
            const std::size_t i = mapV[s0];
            Fimm[i][t] += v;
            if (d == 0) PVT[i][esc] += v;
            else if (W.isvan[d - 1]) PVV[i][mapV[d - 1]] += v;
            else PVT[i][mapT[d - 1]] += v;
        }
    }
    Mat Q;
    if (nV > 0) {
        const Mat IP = minus(eye(nV), PVV);
        if (cond1(IP) > 1e14)
            refuse("the immediate transitions contain a loop of vanishing markings that is never left "
                   "(a timeless trap)");
        const Mat IPinv = inv(IP);
        Q = plus(QTT, mult(QTV, mult(IPinv, PVT)));
        trate = plus(trate, mult(QTV, mult(IPinv, Fimm)));
    } else {
        Q = QTT;
    }
    Vec levT(nT + 1);
    for (std::size_t i = 0; i < nT; ++i) levT[i] = W.S[tgl[i]][p0];
    levT[nT] = Ktop + 3;
    for (std::size_t a = 0; a < nT; ++a)
        for (std::size_t b = 0; b <= nT; ++b)
            if (std::fabs(Q[a][b]) > 1e-14 && std::fabs(levT[b] - levT[a]) > 1 && levT[a] <= Ktop)
                refuse(fmt("iSPN requirement 2 violated: a timed firing followed by immediate firings moves "
                           "the unbounded place from level %.0f to level %.0f, but a level may change by at "
                           "most one",
                           levT[a], levT[b]));
    Elim E;
    E.Q = zeros(nT, nT);
    for (std::size_t i = 0; i < nT; ++i) {
        double rs = 0;
        for (std::size_t j = 0; j <= nT; ++j)
            if (j != i) rs += Q[i][j];
        for (std::size_t j = 0; j < nT; ++j)
            if (j != i) E.Q[i][j] = Q[i][j];
        E.Q[i][i] = -rs;
    }
    E.S = zeros(net.P, nT);
    for (std::size_t p = 0; p < net.P; ++p)
        for (std::size_t i = 0; i < nT; ++i) E.S[p][i] = W.S[tgl[i]][p];
    E.levidx.assign(static_cast<std::size_t>(Ktop + 2), std::vector<std::size_t>());
    for (std::size_t i = 0; i < nT; ++i) {
        const double l = levT[i];
        if (l >= 0 && l < Ktop + 2) E.levidx[static_cast<std::size_t>(l)].push_back(i);
    }
    E.trate = trate;
    E.nvan = nV;
    return E;
}
 
inline int find_K(const Elim& E, int Ktop, const IspnOptions& opt, std::size_t p0) {
    const auto& lv = E.levidx;
    const Mat& S = E.S;
    for (int c = 0; c < Ktop - 2; ++c) {
        const std::vector<std::size_t>& ref = lv[c];
        if (ref.empty()) continue;
        bool ok = true;
        for (int k = c + 1; k < Ktop + 2 && ok; ++k) {
            if (lv[k].size() != ref.size()) {
                ok = false;
                break;
            }
            for (std::size_t p = 0; p < S.size() && ok; ++p) {
                if (p == p0) continue;
                for (std::size_t j = 0; j < ref.size(); ++j)
                    if (S[p][lv[k][j]] != S[p][ref[j]]) {
                        ok = false;
                        break;
                    }
            }
        }
        if (!ok) continue;
        const Mat F0 = sub(E.Q, lv[c], lv[c + 1]);
        const Mat L1 = sub(E.Q, lv[c + 1], lv[c + 1]);
        const Mat D1 = sub(E.Q, lv[c + 1], lv[c]);
        const double sc = std::max(1.0, std::max(max_abs(F0), std::max(max_abs(L1), max_abs(D1))));
        for (int k = c + 1; k <= Ktop && ok; ++k)
            if (max_abs_diff(sub(E.Q, lv[k], lv[k + 1]), F0) > opt.tol * sc ||
                max_abs_diff(sub(E.Q, lv[k], lv[k]), L1) > opt.tol * sc ||
                max_abs_diff(sub(E.Q, lv[k], lv[k - 1]), D1) > opt.tol * sc)
                ok = false;
        if (ok) return c;
    }
    return -1;
}
 
inline double expect(const Vec& r, const std::string& what, const Mat& x,
                     const std::vector<std::vector<std::size_t>>& lv, int K, int Ktop, const Vec& c1,
                     const Vec& c2, double tol) {
    double y = 0;
    for (int i = 0; i <= K; ++i)
        for (std::size_t j = 0; j < x[i].size(); ++j) y += x[i][j] * r[lv[i][j]];
    const std::size_t n = lv[K + 1].size();
    Vec a(n), b(n);
    for (std::size_t j = 0; j < n; ++j) {
        a[j] = r[lv[K + 1][j]];
        b[j] = r[lv[K + 2][j]] - a[j];
    }
    const double sc = std::max(1.0, std::max(max_abs(a), max_abs(b)));
    for (int j = 3; j <= Ktop - K; ++j)
        for (std::size_t i = 0; i < n; ++i)
            if (std::fabs(r[lv[K + j][i]] - a[i] - (j - 1) * b[i]) > tol * sc * j)
                badarg("the " + what + " is not affine in the level above K = " + std::to_string(K) +
                       ", so it has no closed form");
    return y + dot(c1, a) + dot(c2, b);
}
 
inline long long gcd_ll(long long a, long long b) {
    while (b != 0) {
        const long long t = a % b;
        a = b;
        b = t;
    }
    return a;
}
 
inline std::vector<Vec> minimal(const std::vector<Vec>& D, std::size_t T, std::size_t P) {
    const std::size_t n = D.size();
    if (n <= 1) return D;
    std::vector<std::vector<bool>> supp(n, std::vector<bool>(P));
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t p = 0; p < P; ++p) supp[i][p] = D[i][T + p] > 1e-12;
    std::vector<bool> drop(n, false);
    for (std::size_t i = 0; i < n; ++i) {
        if (drop[i]) continue;
        for (std::size_t j = 0; j < n; ++j) {
            if (i == j || drop[j]) continue;
            bool le = true, lt = false;
            for (std::size_t p = 0; p < P; ++p) {
                if (supp[j][p] && !supp[i][p]) le = false;
                if (!supp[j][p] && supp[i][p]) lt = true;
            }
            if (le && (lt || j < i)) {
                drop[i] = true;
                break;
            }
        }
    }
    std::vector<Vec> out;
    for (std::size_t i = 0; i < n; ++i)
        if (!drop[i]) out.push_back(D[i]);
    return out;
}
 
}  // namespace ispn_detail
 
/**
 * Places covered by a non-negative S-invariant (P-semiflow) y' C = 0, by
 * Farkas' algorithm on [C | I] keeping minimal supports. The uncovered places
 * are the candidates for the unbounded place of an iSPN.
 *
 * @param C P x T incidence matrix POST - PRE
 */
inline std::vector<bool> spn_ispn_covered(const std::vector<std::vector<double>>& C) {
    using namespace ispn_detail;
    const std::size_t P = C.size(), T = P > 0 ? C[0].size() : 0;
    std::vector<Vec> D;
    for (std::size_t p = 0; p < P; ++p) {
        if (C[p].size() != T) throw InputError(ErrorKind{"Argument"}, "spn_ispn_covered: C must be P x T");
        Vec row(T + P, 0.0);
        for (std::size_t t = 0; t < T; ++t) row[t] = C[p][t];
        row[T + p] = 1.0;
        D.push_back(row);
    }
    for (std::size_t t = 0; t < T; ++t) {
        std::vector<Vec> nd, pos, neg;
        for (const Vec& r : D) {
            if (std::fabs(r[t]) < 1e-12) nd.push_back(r);
            else if (r[t] > 0) pos.push_back(r);
            else neg.push_back(r);
        }
        for (const Vec& a : pos)
            for (const Vec& b : neg) {
                Vec row(T + P);
                long long g = 0;
                for (std::size_t j = 0; j < T + P; ++j) {
                    row[j] = -b[t] * a[j] + a[t] * b[j];
                    const long long rv = std::llabs(std::llround(row[j]));
                    if (rv != 0) g = gcd_ll(g, rv);
                }
                if (g > 1)
                    for (double& v : row) v /= static_cast<double>(g);
                nd.push_back(row);
            }
        D = minimal(nd, T, P);
    }
    std::vector<bool> covered(P, false);
    for (const Vec& r : D)
        for (std::size_t p = 0; p < P; ++p)
            if (r[T + p] > 1e-12) covered[p] = true;
    return covered;
}
 
/**
 * Steady state of an infinite-state SPN.
 *
 * @param pre  P x T input multiplicities
 * @param post P x T output multiplicities
 * @param inh  P x T inhibitor thresholds (t disabled when m(p) >= inh), +inf for none; empty: none
 * @param rate length-T timed rate, or the weight of an immediate transition
 * @param imm  length-T, true for an immediate transition; empty: all timed
 * @param srv  length-T servers (+inf allowed); empty: all 1. A transition with no input fires at rate(t)
 * @param m0   length-P initial marking
 * @param p0   0-based index of the unbounded place
 */
inline IspnResult spn_ispn(const std::vector<std::vector<double>>& pre,
                           const std::vector<std::vector<double>>& post,
                           std::vector<std::vector<double>> inh, const std::vector<double>& rate,
                           std::vector<bool> imm, std::vector<double> srv, const std::vector<double>& m0,
                           std::size_t p0, const IspnOptions& opt = IspnOptions()) {
    using namespace ispn_detail;
    const double Inf = std::numeric_limits<double>::infinity();
    if (opt.extra < 3) badarg("OPTIONS.extra must be at least 3");
    const std::size_t P = pre.size(), T = P > 0 ? pre[0].size() : 0;
    if (post.size() != P) badarg("PRE and POST must both be P x T");
    for (std::size_t p = 0; p < P; ++p)
        if (pre[p].size() != T || post[p].size() != T) badarg("PRE and POST must both be P x T");
    if (inh.empty()) inh.assign(P, Vec(T, Inf));
    if (imm.empty()) imm.assign(T, false);
    if (srv.empty()) srv.assign(T, 1.0);
    bool shape = inh.size() == P && rate.size() == T && imm.size() == T && srv.size() == T && m0.size() == P;
    for (std::size_t p = 0; shape && p < P; ++p) shape = inh[p].size() == T;
    if (!shape) badarg("INH must be P x T, RATE, IMM and SRV of length T and M0 of length P");
    for (std::size_t p = 0; p < P; ++p) {
        if (m0[p] < 0 || m0[p] != std::rint(m0[p]))
            badarg("arc multiplicities and the initial marking must be non-negative integers");
        for (std::size_t t = 0; t < T; ++t)
            if (pre[p][t] < 0 || post[p][t] < 0 || pre[p][t] != std::rint(pre[p][t]) ||
                post[p][t] != std::rint(post[p][t]))
                badarg("arc multiplicities and the initial marking must be non-negative integers");
    }
    for (std::size_t t = 0; t < T; ++t)
        if (rate[t] < 0 || srv[t] < 1)
            badarg("arc multiplicities, the initial marking and the rates must be non-negative, servers >= 1");
    if (p0 >= P) badarg("P0 must be the index of a place");
 
    for (std::size_t t = 0; t < T; ++t) {
        const double d = post[p0][t] - pre[p0][t];
        if (std::fabs(d) > 1)
            refuse(fmt("iSPN requirement 2 violated: transition %.0f changes the unbounded place %.0f by %.0f "
                       "tokens, but a firing may move the level by at most one",
                       static_cast<double>(t), static_cast<double>(p0), d));
    }
    std::vector<std::size_t> other;
    for (std::size_t p = 0; p < P; ++p)
        if (p != p0) other.push_back(p);
 
    int Kstr = 0;
    std::vector<bool> infdeg(T, false);
    for (std::size_t t = 0; t < T; ++t) {
        if (pre[p0][t] > 0) {
            Kstr = std::max(Kstr, static_cast<int>(pre[p0][t]));
            if (!std::isinf(srv[t])) {
                Kstr = std::max(Kstr, static_cast<int>(pre[p0][t] * srv[t]));
            } else {
                bool any_other = false;
                for (std::size_t p : other) any_other = any_other || pre[p][t] > 0;
                if (any_other)
                    infdeg[t] = true;
                else
                    refuse(fmt("iSPN requirement 3 violated: transition %.0f is infinite-server and consumes "
                               "only from the unbounded place %.0f, so its rate grows with the level and never "
                               "becomes level independent",
                               static_cast<double>(t), static_cast<double>(p0)));
            }
        }
        if (!std::isinf(inh[p0][t])) Kstr = std::max(Kstr, static_cast<int>(inh[p0][t]));
    }
 
    Net net;
    net.pre = &pre;
    net.post = &post;
    net.inh = &inh;
    net.rate = &rate;
    net.imm = &imm;
    net.srv = &srv;
    net.p0 = p0;
    net.P = P;
    net.T = T;
    net.other = other;
 
    int Kneed = Kstr;
    int Ktop = Kneed + opt.extra;
    Elim E;
    int K;
    while (true) {
        if (Ktop > opt.maxlevel)
            refuse(fmt("iSPN requirement 1/3 violated: the levels of place %.0f do not become identical (same "
                       "phase markings and level-independent rates) up to level %.0f; either another place is "
                       "unbounded or the rates depend on the level. Raise OPTIONS.maxlevel if the threshold is "
                       "higher",
                       static_cast<double>(p0), opt.maxlevel));
        const Walk W = explore(net, m0, Ktop + 2, opt);
        for (std::size_t t = 0; t < T; ++t) {
            if (!infdeg[t]) continue;
            double mx = 0;
            for (const Vec& m : W.S) {
                double dg = Inf;
                for (std::size_t p : other)
                    if (pre[p][t] > 0) dg = std::min(dg, std::floor(m[p] / pre[p][t]));
                mx = std::max(mx, dg);
            }
            Kneed = std::max(Kneed, static_cast<int>(pre[p0][t] * mx));
        }
        if (Ktop < Kneed + opt.extra) {
            Ktop = Kneed + opt.extra;
            continue;
        }
        E = eliminate(W, net, Ktop);
        K = find_K(E, Ktop, opt, p0);
        if (K >= 0) break;
        Ktop = 2 * Ktop + 1;
    }
 
    const auto& lv = E.levidx;
    const Mat& Qt = E.Q;
    const Mat A0 = sub(Qt, lv[K], lv[K + 1]);
    const Mat A1 = sub(Qt, lv[K + 1], lv[K + 1]);
    const Mat A2 = sub(Qt, lv[K + 1], lv[K]);
    const std::size_t nL = lv[K].size();
    if (all_zero(A0))
        refuse(fmt("place %.0f is bounded: no transition raises its marking above level %.0f. The net has a "
                   "finite state space; solve it with SolverCTMC",
                   static_cast<double>(p0), K));
    if (all_zero(A2))
        unstable(fmt("the iSPN is unstable: no transition lowers the marking of place %.0f above level %.0f, "
                     "so it grows without bound",
                     static_cast<double>(p0), K));
 
    // drift (Haverkort (17.7))
    Mat Aa = plus(plus(A0, A1), A2);
    for (std::size_t i = 0; i < nL; ++i) Aa[i][nL - 1] = 1.0;
    Vec eLast(nL, 0.0);
    eLast[nL - 1] = 1.0;
    const Vec piA = solve_left(Aa, eLast);
    const double up = dot(piA, row_sums(A0));
    const double down = dot(piA, row_sums(A2));
    if (!(up < down * (1 - 1e-12)))
        unstable(fmt("the iSPN is unstable: the mean drift of place %.0f above level %.0f is up %.6g >= down "
                     "%.6g (Haverkort (17.7)); no steady state exists",
                     static_cast<double>(p0), K, up, down));
 
    std::string qm = opt.qbdmethod;
    std::transform(qm.begin(), qm.end(), qm.begin(), [](unsigned char c) { return std::tolower(c); });
    Matrix<double> Rmat;
    if (qm == "logred")
        Rmat = mam::qbd_R_logred(to_matrix(A2), to_matrix(A1), to_matrix(A0), 100000u, 1e-12);
    else if (qm == "successive")
        Rmat = mam::qbd_R(to_matrix(A2), to_matrix(A1), to_matrix(A0));
    else
        badarg("unknown OPTIONS.qbdmethod '" + opt.qbdmethod + "'; use 'logred' or 'successive'");
    const Mat Rm = from_matrix(Rmat);
    const double spR = spectral_radius(Rmat);
    if (!(spR < 1)) unstable(fmt("the iSPN is unstable: the spectral radius of R is %.6g >= 1", spR));
    const Mat I = eye(nL);
    const Mat Ni = inv(minus(I, Rm));
 
    // boundary
    std::vector<std::size_t> bidx;
    for (int i = 0; i <= K; ++i) bidx.insert(bidx.end(), lv[i].begin(), lv[i].end());
    const std::size_t nb = bidx.size();
    const Mat B = sub(Qt, bidx, bidx);
    Mat Bm = B;
    const Mat RA2 = mult(Rm, A2);
    for (std::size_t i = 0; i < nL; ++i)
        for (std::size_t j = 0; j < nL; ++j) Bm[nb - nL + i][nb - nL + j] += RA2[i][j];
    const Vec NiOne = row_sums(Ni);
    for (std::size_t i = 0; i < nb; ++i) Bm[i][nb - 1] = 1.0;
    for (std::size_t i = 0; i < nL; ++i) Bm[nb - nL + i][nb - 1] = NiOne[i];
    Vec rhs(nb, 0.0);
    rhs[nb - 1] = 1.0;
    Vec xb = solve_left(Bm, rhs);
    for (double& v : xb)
        if (std::fabs(v) < 1e-300) v = 0.0;
    IspnResult res;
    for (int i = 0, o = 0; i <= K; ++i) {
        res.x.push_back(Vec(xb.begin() + o, xb.begin() + o + static_cast<long>(lv[i].size())));
        o += static_cast<int>(lv[i].size());
    }
    const Vec& xK = res.x[K];
    const Mat RN = mult(Rm, Ni);
    const Vec c1 = vec_mult(xK, RN);
    const Vec c2 = vec_mult(xK, mult(mult(Rm, Rm), mult(Ni, Ni)));
 
    const Mat& S = E.S;
    const std::size_t nT = S[0].size();
    Vec mm(P);
    for (std::size_t p = 0; p < P; ++p)
        mm[p] = expect(S[p], "marking of place " + std::to_string(p), res.x, lv, K, Ktop, c1, c2, opt.tol);
    double EL2 = 0;
    for (int i = 0; i <= K; ++i) EL2 += static_cast<double>(i) * i * sum(res.x[i]);
    const Mat N2 = mult(Ni, Ni);
    const Mat N3 = mult(N2, Ni);
    const Mat term = plus(plus(scale(RN, static_cast<double>(K) * K), scale(mult(Rm, N2), 2.0 * K)),
                          mult(mult(Rm, plus(I, Rm)), N3));
    EL2 += dot(xK, row_sums(term));
    res.varlevel = std::max(EL2 - mm[p0] * mm[p0], 0.0);
 
    Vec dist;
    double tail = 1.0;
    for (int i = 0; i <= K; ++i) {
        dist.push_back(sum(res.x[i]));
        tail -= dist.back();
    }
    Vec v = vec_mult(xK, Rm);
    double lvl = K;
    while (lvl < opt.distmax && tail > 1e-12) {
        lvl += 1;
        const double s = sum(v);
        dist.push_back(s);
        tail -= s;
        v = vec_mult(v, Rm);
    }
    if (!std::isinf(opt.distmax))
        while (static_cast<double>(dist.size()) > std::floor(opt.distmax) + 1) dist.pop_back();
    res.leveldist = dist;
    res.leveltail = std::max(1.0 - sum(dist), 0.0);
 
    for (std::size_t p = 0; p < P; ++p) {
        if (p == p0) {
            res.marginal.push_back(res.leveldist);
            continue;
        }
        std::size_t nmax = 0;
        for (std::size_t s = 0; s < nT; ++s) nmax = std::max(nmax, static_cast<std::size_t>(S[p][s]));
        Vec mg(nmax + 1);
        for (std::size_t n = 0; n <= nmax; ++n) {
            Vec ind(nT);
            for (std::size_t s = 0; s < nT; ++s) ind[s] = S[p][s] == static_cast<double>(n) ? 1.0 : 0.0;
            mg[n] = expect(ind, "marginal of place " + std::to_string(p), res.x, lv, K, Ktop, c1, c2, opt.tol);
        }
        res.marginal.push_back(mg);
    }
 
    res.tput.assign(T, 0.0);
    for (std::size_t t = 0; t < T; ++t) {
        Vec col(nT);
        for (std::size_t s = 0; s < nT; ++s) col[s] = E.trate[s][t];
        res.tput[t] = expect(col, "firing rate of transition " + std::to_string(t), res.x, lv, K, Ktop, c1, c2,
                             opt.tol);
    }
 
    res.reward.assign(opt.rewards.size(), 0.0);
    for (std::size_t j = 0; j < opt.rewards.size(); ++j) {
        Vec rv(nT);
        for (std::size_t s = 0; s < nT; ++s) {
            Vec m(P);
            for (std::size_t p = 0; p < P; ++p) m[p] = S[p][s];
            rv[s] = opt.rewards[j](m);
        }
        res.reward[j] = expect(rv, "reward " + std::to_string(j), res.x, lv, K, Ktop, c1, c2, opt.tol);
    }
 
    res.K = K;
    res.L = nL;
    res.A0 = A0;
    res.A1 = A1;
    res.A2 = A2;
    res.R = Rm;
    res.B = B;
    for (int i = 0; i <= K; ++i) {
        Mat lm = zeros(P, lv[i].size());
        for (std::size_t p = 0; p < P; ++p)
            for (std::size_t j = 0; j < lv[i].size(); ++j) lm[p][j] = S[p][lv[i][j]];
        res.levelmarkings.push_back(lm);
    }
    res.drift = {up, down};
    res.spR = spR;
    res.meanmarking = mm;
    res.ntangible = nT;
    res.nvanishing = E.nvan;
    return res;
}
 
}  // namespace spn
}  // namespace line
 
#endif  // LINE_API_SPN_SPN_ISPN_H