solver_mam_ispn.h

Go to the documentation of this file.

/**
 * @file solver_mam_ispn.h
 * The 'ispn' method of SolverMAM: steady state of a stochastic Petri net with
 * exactly one unbounded place, solved as a quasi-birth-death process by the
 * matrix-geometric method (Haverkort 1998, ch. 17; spn_ispn).
 *
 * Port of matlab/src/solvers/MAM/solver_mam_ispn.m (Java twin
 * jline.solvers.mam.analyzers.Solver_mam_ispn, Python twin
 * line_solver.api.solvers.mam.ispn).
 *
 * The net is flattened to (place, class) levels, place-major as in spn_mdd, so
 * level pp*R+k (0-based) is class k in place pp. An open class arriving from a
 * Source to a Place becomes a timed transition with no input place, firing at
 * the arrival rate times the routing probability; tokens a mode sends to a Sink
 * are dropped. The unbounded level is MamOptions::ispn_place when set (a place
 * name, or 'place.class' in a multiclass net), and otherwise the one level that
 * no non-negative S-invariant covers.
 *
 * The measures are those SolverCTMC reports for a Petri net: Q(ist,k) is the
 * mean number of class-k tokens in the place, U = Q (a Place is an INF
 * station), Tp(ist,k) the token rate out of the place (the sum over the modes
 * of their firing rate times the class-k input multiplicity), R = Q/Tp by
 * Little's law. The Source row carries the arrival rate.
 *
 * spn_ispn works at double; another arithmetic is converted to double on the
 * way in and back on the way out, since the matrix-geometric solve is a
 * floating-point iteration (the R matrix has no closed form in general).
 *
 * Copyright (c) 2012-2026, Imperial College London
 * All rights reserved.
 */
 
#ifndef LINE_SOLVERS_MAM_SOLVER_MAM_ISPN_H
#define LINE_SOLVERS_MAM_SOLVER_MAM_ISPN_H
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <map>
#include <string>
#include <vector>
 
#include "line/api/spn/spn_ispn.h"
#include "line/lang/distribution.h"
#include "line/lang/qn/network_struct.h"
#include "line/solvers/mam/mam_types.h"
#include "line/util/error.h"
 
namespace line {
namespace mam {
 
/** True when the model holds a Place node, i.e. it is a stochastic Petri net. */
template <class T>
bool ispn_has_place(const qn::NetworkStruct<T>& L) {
    for (const qn::NodeDef& nd : L.nodes)
        if (nd.nodetype == qn::NodeType::Place) return true;
    return false;
}
 
/**
 * Refuses the SolverMAM entry point ENTRY on a stochastic Petri net model with
 * LINE:unsupported:Place. The ispn route returns steady-state averages only, so
 * no state probability, response time distribution, transient or MAM result
 * entry point can answer on a Place model. MATLAB solver_mam_ispn_refuse.m, Java
 * SolverMAM.refusePlace and Python _refuse_place raise the same id and wording.
 */
template <class T>
void ispn_refuse_entry(const qn::NetworkStruct<T>& L, const std::string& entry) {
    if (ispn_has_place(L))
        throw UnsupportedError(ErrorKind{"Place"},
            "SolverMAM: " + entry + " has no stochastic Petri net analysis; the ispn method "
            "reports steady-state averages only");
}
 
/**
 * Why the iSPN route cannot take this model under METHOD, or "" when it can or
 * when the model is not a net. 'ispn' needs a Place model; a Place model under
 * 'default' or 'ispn' must be a pure net (Place, Transition, Source, Sink); any
 * other method on a Place model is refused, since no queueing analyzer of
 * SolverMAM reads a Place. Mirrors MATLAB SolverMAM.supportsModelMethod and
 * Java SolverMAM.ispnRefusal.
 */
template <class T>
std::string ispn_refusal(const qn::NetworkStruct<T>& L, const std::string& method) {
    const bool hp = ispn_has_place(L);
    if (method == "ispn" && !hp)
        return "SolverMAM: the ispn method solves stochastic Petri nets and the model holds no "
               "Place node. Use the 'default' method of SolverMAM for a queueing network.";
    if (!hp) return std::string();
    if (method != "default" && method != "ispn")
        return "SolverMAM: the " + method + " method does not support Place nodes: on a "
               "stochastic Petri net SolverMAM offers only the 'ispn' method, which 'default' "
               "selects.";
    for (const qn::NodeDef& nd : L.nodes)
        if (nd.nodetype != qn::NodeType::Place && nd.nodetype != qn::NodeType::Transition &&
            nd.nodetype != qn::NodeType::Source && nd.nodetype != qn::NodeType::Sink)
            return "SolverMAM: node " + nd.name + " is neither a Place, a Transition, a Source "
                   "nor a Sink; the ispn method of SolverMAM solves pure stochastic Petri nets. "
                   "Use SolverCTMC for a mixed queueing and Petri net model.";
    return std::string();
}
 
/** The analyzer. Returns the MvaSolution layout with actualmethod "ispn". */
template <class T>
MamSolution<T> solver_mam_ispn(const qn::NetworkStruct<T>& L, const MamOptions& opt) {
    using Vec = std::vector<double>;
    using Mat = std::vector<Vec>;
    const double Inf = std::numeric_limits<double>::infinity();
    auto D = [](const T& v) { return num_traits<T>::to_double(v); };
    auto structure = [](const std::string& m) {
        throw UnsupportedError(ErrorKind{"Structure"}, "solver_mam_ispn: " + m);
    };
 
    const std::size_t R = L.nclasses;
    const std::size_t nn = L.nodes.size();
    std::vector<std::size_t> places, transitions;  // 1-based node indices
    for (std::size_t i = 1; i <= nn; ++i) {
        const qn::NodeType t = L.nodes[i - 1].nodetype;
        if (t == qn::NodeType::Place) places.push_back(i);
        else if (t == qn::NodeType::Transition) transitions.push_back(i);
        else if (t != qn::NodeType::Source && t != qn::NodeType::Sink)
            structure("node " + L.nodes[i - 1].name + " is neither a Place, a Transition, a "
                      "Source nor a Sink; the ispn method solves pure stochastic Petri nets");
    }
    if (places.empty())
        throw InputError(ErrorKind{"Argument"},
                         "solver_mam_ispn: the ispn method needs a model holding Place nodes");
    const std::size_t NP = places.size();
    const std::size_t PL = NP * R;
    auto lvl = [R](std::size_t pp, std::size_t k) { return pp * R + k; };
    auto place_of = [&places](std::size_t node) -> long {
        for (std::size_t pp = 0; pp < places.size(); ++pp)
            if (places[pp] == node) return static_cast<long>(pp);
        return -1;
    };
 
    // ---- modes of the transitions, one column each
    Mat pre(PL), post(PL), inh(PL);
    Vec rate, srv;
    std::vector<bool> imm;
    std::size_t ncol = 0;
    auto add_col = [&]() {
        for (std::size_t l = 0; l < PL; ++l) {
            pre[l].push_back(0.0);
            post[l].push_back(0.0);
            inh[l].push_back(Inf);
        }
        return ncol++;
    };
    for (const std::size_t ind : transitions) {
        const auto it = L.transparam.find(ind);
        if (it == L.transparam.end()) continue;
        const qn::TransitionParam<T>& tp = it->second;
        const std::string& tname = L.nodes[ind - 1].name;
        for (std::size_t m = 0; m < tp.nmodes; ++m) {
            if (m < tp.firingdep.size() && tp.firingdep[m])
                structure("mode " + std::to_string(m + 1) + " of transition " + tname +
                          " has a marking-dependent firing rate, which the ispn method does "
                          "not support");
            const bool isimm = m < tp.timing.size() && tp.timing[m] == lang::TimingStrategy::IMMEDIATE;
            double r = 1.0;
            if (isimm) {
                if (m < tp.fireweight.size()) r = D(tp.fireweight[m]);
            } else {
                bool expo = m < tp.firingproc.size() && !tp.firingproc[m].disabled;
                if (expo) {
                    const mam::Map<T> proc = lang::dist_to_map(tp.firingproc[m]);
                    expo = proc.order() == 1;
                    if (expo) r = D(proc.D1(0, 0));
                }
                if (!expo)
                    throw UnsupportedError(ErrorKind{"NonExponential"},
                        "solver_mam_ispn: mode " + std::to_string(m + 1) + " of transition " +
                        tname + " has a non-exponential firing time; the ispn method needs "
                        "exponential or immediate transitions");
            }
            const Matrix<T>* en = m < tp.enabling.size() ? &tp.enabling[m] : nullptr;
            const Matrix<T>* ih = m < tp.inhibiting.size() ? &tp.inhibiting[m] : nullptr;
            const Matrix<T>* fi = m < tp.firing.size() ? &tp.firing[m] : nullptr;
            const std::size_t c = add_col();
            for (std::size_t pp = 0; pp < NP; ++pp) {
                const std::size_t q = places[pp] - 1;  // arcs are indexed by node
                for (std::size_t k = 0; k < R; ++k) {
                    const std::size_t l = lvl(pp, k);
                    if (en && q < en->rows() && k < en->cols()) pre[l][c] = std::max(0.0, D((*en)(q, k)));
                    if (fi && q < fi->rows() && k < fi->cols()) post[l][c] = std::max(0.0, D((*fi)(q, k)));
                    if (ih && q < ih->rows() && k < ih->cols()) {
                        const double v = D((*ih)(q, k));
                        // a zero threshold is the "no arc" filler of older models
                        inh[l][c] = (v <= 0 || std::isnan(v)) ? Inf : v;
                    }
                }
            }
            double ns = m < tp.nmodeservers.size() ? tp.nmodeservers[m] : 1.0;
            if (ns >= 1e9) ns = Inf;  // GlobalConstants.MaxInt stands for infinite servers
            rate.push_back(r);
            imm.push_back(isimm);
            srv.push_back(ns);
        }
    }
    const std::size_t nmodes = ncol;
 
    // ---- arrivals from Sources, as transitions with no input place
    Vec lambda(R, 0.0);
    std::vector<std::size_t> srcst(R, 0);  // 1-based station, 0 = none
    for (std::size_t ind = 1; ind <= nn; ++ind) {
        if (L.nodes[ind - 1].nodetype != qn::NodeType::Source) continue;
        const std::size_t ist = L.nodes[ind - 1].station;
        if (ist == 0) continue;
        for (std::size_t k = 0; k < R; ++k) {
            if (L.disabled[ist - 1][k]) continue;
            const double lk = D(L.rates(ist - 1, k));
            if (std::isnan(lk) || lk <= 0) continue;
            if (L.procid(ist, k + 1) != lang::ProcessType::EXP)
                throw UnsupportedError(ErrorKind{"NonExponential"},
                    "solver_mam_ispn: class " + L.classes[k].name + " arrives at the Source with "
                    "a non-exponential process; the ispn method needs Poisson arrivals");
            lambda[k] += lk;
            srcst[k] = ist;
            for (std::size_t jnd = 1; jnd <= nn; ++jnd)
                for (std::size_t r = 0; r < R; ++r) {
                    const std::size_t row = (ind - 1) * R + k, col = (jnd - 1) * R + r;
                    if (row >= L.rtnodes.rows() || col >= L.rtnodes.cols()) continue;
                    const double pr = D(L.rtnodes(row, col));
                    if (!(pr > 0)) continue;
                    const long pp = place_of(jnd);
                    if (pp < 0)
                        structure("the Source routes class " + L.classes[k].name + " to " +
                                  L.nodes[jnd - 1].name + ", which is not a Place");
                    const std::size_t c = add_col();
                    post[lvl(static_cast<std::size_t>(pp), r)][c] = 1.0;
                    rate.push_back(lk * pr);
                    imm.push_back(false);
                    srv.push_back(1.0);
                }
        }
    }
 
    // ---- initial marking: the declared place state, else njobs at the reference station
    Vec m0(PL, 0.0);
    bool anyset = false;
    for (std::size_t pp = 0; pp < NP; ++pp) {
        const auto it = L.initmarking.find(places[pp]);
        if (it == L.initmarking.end()) continue;
        for (std::size_t k = 0; k < R && k < it->second.size(); ++k) {
            const double v = D(it->second[k]);
            m0[lvl(pp, k)] = v;
            anyset = anyset || v != 0;
        }
    }
    if (!anyset) {
        std::fill(m0.begin(), m0.end(), 0.0);
        for (std::size_t k = 0; k < R; ++k) {
            const double nj = L.classes[k].population;
            if (!std::isfinite(nj) || nj <= 0) continue;
            const std::size_t rs = L.classes[k].refstat;
            if (rs == 0 || rs > L.nstations) continue;
            for (std::size_t pp = 0; pp < NP; ++pp)
                if (L.nodes[places[pp] - 1].station == rs) {
                    m0[lvl(pp, k)] = nj;
                    break;
                }
        }
    }
 
    // ---- level names and the unbounded level
    std::vector<std::string> lname(PL);
    for (std::size_t pp = 0; pp < NP; ++pp)
        for (std::size_t k = 0; k < R; ++k)
            lname[lvl(pp, k)] = R == 1 ? L.nodes[places[pp] - 1].name
                                       : L.nodes[places[pp] - 1].name + "." + L.classes[k].name;
    auto join = [&lname](const std::vector<std::size_t>& ix) {
        std::string s;
        for (std::size_t i = 0; i < ix.size(); ++i) s += (i ? ", " : "") + lname[ix[i]];
        return s;
    };
    auto used_level = [&](std::size_t l) {
        for (std::size_t c = 0; c < ncol; ++c)
            if (pre[l][c] != 0 || post[l][c] != 0) return true;
        return false;
    };
    std::size_t p0 = PL;
    if (!opt.ispn_place.empty()) {
        const std::string& want = opt.ispn_place;
        for (std::size_t l = 0; l < PL && p0 == PL; ++l)
            if (lname[l] == want) p0 = l;
        if (p0 == PL && R > 1) {
            for (std::size_t pp = 0; pp < NP; ++pp) {
                if (L.nodes[places[pp] - 1].name != want) continue;
                std::vector<std::size_t> used;
                for (std::size_t k = 0; k < R; ++k)
                    if (used_level(lvl(pp, k))) used.push_back(lvl(pp, k));
                if (used.size() == 1) p0 = used[0];
                break;
            }
        }
        if (p0 == PL) {
            std::vector<std::size_t> all(PL);
            for (std::size_t l = 0; l < PL; ++l) all[l] = l;
            throw InputError(ErrorKind{"Argument"},
                "solver_mam_ispn: options.config.ispn_place = '" + want + "' names no place "
                "level; the levels are " + join(all));
        }
    } else {
        Mat C(PL, Vec(ncol, 0.0));
        for (std::size_t l = 0; l < PL; ++l)
            for (std::size_t c = 0; c < ncol; ++c) C[l][c] = post[l][c] - pre[l][c];
        const std::vector<bool> covered = spn::spn_ispn_covered(C);
        std::vector<std::size_t> cand;
        for (std::size_t l = 0; l < PL; ++l)
            if (!covered[l] && (used_level(l) || m0[l] > 0)) cand.push_back(l);
        if (cand.empty())
            structure("every place is covered by a non-negative S-invariant, so the net is "
                      "bounded and has a finite state space; solve it with SolverCTMC");
        if (cand.size() > 1)
            structure("the iSPN needs exactly one unbounded place, but " +
                      std::to_string(cand.size()) + " places are not covered by an S-invariant (" +
                      join(cand) + "); name the unbounded one with options.config.ispn_place");
        p0 = cand[0];
    }
 
    // ---- solve
    spn::IspnOptions iopt;
    if (opt.ispn_maxlevel > 0) iopt.maxlevel = static_cast<int>(opt.ispn_maxlevel);
    const spn::IspnResult res = spn::spn_ispn(pre, post, inh, rate, imm, srv, m0, p0, iopt);
 
    // ---- pack the analyzer contract
    const std::size_t M = L.nstations;
    const T zero = num_traits<T>::from_int(0);
    auto F = [](double v) { return num_traits<T>::from_double(v); };
    MamSolution<T> out;
    mva::MvaSolution<T>& s = out.sol;
    s.Q = Matrix<T>(M, R, zero);
    s.U = Matrix<T>(M, R, zero);
    s.R = Matrix<T>(M, R, zero);
    s.Tp = Matrix<T>(M, R, zero);
    Matrix<double> Q(M, R, 0.0), Tp(M, R, 0.0);
    for (std::size_t pp = 0; pp < NP; ++pp) {
        const std::size_t ist = L.nodes[places[pp] - 1].station;
        if (ist == 0) continue;
        for (std::size_t k = 0; k < R; ++k) {
            const std::size_t l = lvl(pp, k);
            double tk = 0;
            for (std::size_t c = 0; c < nmodes && c < res.tput.size(); ++c) tk += res.tput[c] * pre[l][c];
            Q(ist - 1, k) = res.meanmarking[l];
            Tp(ist - 1, k) = tk;
        }
    }
    Vec X(R, 0.0), Cv(R, 0.0);
    for (std::size_t k = 0; k < R; ++k) {
        if (srcst[k] > 0) {
            Tp(srcst[k] - 1, k) = lambda[k];
            X[k] = lambda[k];
        } else {
            const std::size_t rs = L.classes[k].refstat;
            if (rs >= 1 && rs <= M) X[k] = Tp(rs - 1, k);
        }
        double qs = 0;
        for (std::size_t i = 0; i < M; ++i) qs += Q(i, k);
        if (X[k] > 0) Cv[k] = qs / X[k];
    }
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t k = 0; k < R; ++k) {
            s.Q(i, k) = F(Q(i, k));
            s.U(i, k) = F(Q(i, k));
            s.Tp(i, k) = F(Tp(i, k));
            if (Tp(i, k) > 0) s.R(i, k) = F(Q(i, k) / Tp(i, k));
        }
    s.X.assign(R, zero);
    s.C.assign(R, zero);
    for (std::size_t k = 0; k < R; ++k) {
        s.X[k] = F(X[k]);
        s.C[k] = F(Cv[k]);
    }
    s.iter = 1;
    s.method = "ispn";
    out.actualmethod = "ispn";
    return out;
}
 
}  // namespace mam
}  // namespace line
 
#endif  // LINE_SOLVERS_MAM_SOLVER_MAM_ISPN_H