spn_net.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_NET_H
#define LINE_API_SPN_SPN_NET_H
 
/**
 * @file
 * @ingroup api_spn
 * Shared net model of the Petri-net property functions (spn_tinvariants,
 * spn_siphons, spn_traps, spn_structprops, spn_reachgraph, spn_reachprops,
 * spn_properties): the level-by-event arc matrices, the initial marking,
 * Farkas' algorithm, the siphon enumeration and the GSPN firing rule. The
 * MATLAB twin is matlab/src/api/spn/private/.
 *
 * Levels are the (place, class) pairs of `spn_sinvariants`, place-major
 * (pp * R + k, 0-based); events are the (transition, mode) pairs in node
 * order. Per event, pre = max(0, enabling), post = max(0, firing) and the
 * inhibitor threshold is inhibiting (infinite: no arc).
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <map>
#include <string>
#include <vector>
 
#include "line/api/spn/spn_sinvariants.h"
#include "line/lang/lang_types.h"
#include "line/lang/qn/network_struct.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace spn {
namespace detail {
 
/** Arc matrices of a net over place levels and events. */
struct SpnNet {
    /** 1-based node indices of the places, in level order. */
    std::vector<std::size_t> places;
    std::size_t R = 1;
    std::size_t n = 0;
    std::size_t E = 0;
    /** 1-based node index and 0-based mode index of each event. */
    std::vector<std::size_t> transitions, modes;
    /** pre[l][e], post[l][e], C = post - pre, inh[l][e] (infinite: no arc). */
    std::vector<std::vector<long long>> pre, post, C;
    std::vector<std::vector<double>> inh;
    std::vector<bool> immediate;
    std::vector<double> prio, weight;
    /** Features the structural checks ignore: inhibitor, immediate, priority. */
    std::vector<std::string> ignored;
 
    bool monotone() const { return ignored.empty(); }
 
    std::vector<bool> enabled(const std::vector<long long>& m) const {
        std::vector<bool> en(E, true);
        for (std::size_t e = 0; e < E; ++e)
            for (std::size_t l = 0; l < n && en[e]; ++l)
                if (pre[l][e] > m[l] || inh[l][e] <= static_cast<double>(m[l])) en[e] = false;
        return en;
    }
 
    /**
     * Events that may fire in m: the enabled immediate events of the highest
     * priority when any immediate event is enabled, else all enabled events.
     */
    std::vector<bool> fireable(const std::vector<long long>& m) const {
        std::vector<bool> en = enabled(m);
        bool anyimm = false;
        double best = -std::numeric_limits<double>::infinity();
        for (std::size_t e = 0; e < E; ++e)
            if (en[e] && immediate[e]) {
                anyimm = true;
                best = std::max(best, prio[e]);
            }
        if (!anyimm) return en;
        std::vector<bool> f(E, false);
        for (std::size_t e = 0; e < E; ++e) f[e] = en[e] && immediate[e] && prio[e] == best;
        return f;
    }
};
 
inline long long spn_net_integer(double x, const char* what, const std::string& caller) {
    const double r = std::floor(x + 0.5);
    if (std::fabs(x - r) > 1e-9)
        throw InputError(ErrorKind{"Model"}, caller + ": " + what +
                                                 " is not integral; a fractional arc multiplicity "
                                                 "has no Petri-net meaning");
    return static_cast<long long>(r);
}
 
template <class T>
SpnNet spn_build_net(const qn::NetworkStruct<T>& sn, const std::string& caller) {
    std::vector<std::size_t> places, trans;
    for (std::size_t i = 1; i <= sn.nodes.size(); ++i) {
        if (sn.nodes[i - 1].nodetype == lang::NodeType::Place) places.push_back(i);
        else if (sn.nodes[i - 1].nodetype == lang::NodeType::Transition) trans.push_back(i);
    }
    if (places.empty() || trans.empty())
        throw InputError(ErrorKind{"Model"}, caller + ": the model holds no Place or no Transition node");
    const std::size_t R = sn.nclasses;
    for (std::size_t k = 0; k < R && k < sn.classes.size(); ++k)
        if (!std::isfinite(sn.classes[k].population))
            throw UnsupportedError(ErrorKind{"OpenClass"},
                                   caller + ": class " + std::to_string(k + 1) +
                                       " is open: its arrivals enter a place through a Source "
                                       "routing rather than a transition, so the incidence "
                                       "matrix does not describe the token game");
    SpnNet net;
    net.places = places;
    net.R = R;
    const std::size_t P = places.size();
    net.n = P * R;
    const std::size_t n = net.n;
    net.pre.assign(n, std::vector<long long>());
    net.post.assign(n, std::vector<long long>());
    net.inh.assign(n, std::vector<double>());
    const double inf = std::numeric_limits<double>::infinity();
    for (std::size_t t = 0; t < trans.size(); ++t) {
        const typename std::map<std::size_t, qn::TransitionParam<T>>::const_iterator it =
            sn.transparam.find(trans[t]);
        if (it == sn.transparam.end()) continue;
        const qn::TransitionParam<T>& tp = it->second;
        for (std::size_t m = 0; m < tp.nmodes; ++m) {
            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;
            for (std::size_t pp = 0; pp < P; ++pp) {
                const std::size_t q = places[pp] - 1;
                for (std::size_t k = 0; k < R; ++k) {
                    const std::size_t l = pp * R + k;
                    double a = 0, b = 0, h = inf;
                    if (en && q < en->rows() && k < en->cols())
                        a = std::max(0.0, num_traits<T>::to_double((*en)(q, k)));
                    if (fi && q < fi->rows() && k < fi->cols())
                        b = std::max(0.0, num_traits<T>::to_double((*fi)(q, k)));
                    if (ih && q < ih->rows() && k < ih->cols())
                        h = num_traits<T>::to_double((*ih)(q, k));
                    net.pre[l].push_back(spn_net_integer(a, "an enabling arc", caller));
                    net.post[l].push_back(spn_net_integer(b, "a firing arc", caller));
                    net.inh[l].push_back(h);
                }
            }
            net.transitions.push_back(trans[t]);
            net.modes.push_back(m);
            net.immediate.push_back(m < tp.timing.size() &&
                                    tp.timing[m] == lang::TimingStrategy::IMMEDIATE);
            net.prio.push_back(m < tp.firingprio.size() ? tp.firingprio[m] : 1.0);
            net.weight.push_back(m < tp.fireweight.size()
                                     ? num_traits<T>::to_double(tp.fireweight[m])
                                     : 1.0);
        }
    }
    net.E = net.transitions.size();
    if (net.E == 0) throw InputError(ErrorKind{"Model"}, caller + ": the net has no firing mode");
    net.C.assign(n, std::vector<long long>(net.E, 0));
    bool anyinh = false;
    for (std::size_t l = 0; l < n; ++l)
        for (std::size_t e = 0; e < net.E; ++e) {
            net.C[l][e] = net.post[l][e] - net.pre[l][e];
            anyinh = anyinh || std::isfinite(net.inh[l][e]);
        }
    bool anyimm = false, anytimed = false;
    std::vector<double> ip;
    for (std::size_t e = 0; e < net.E; ++e) {
        if (net.immediate[e]) {
            anyimm = true;
            if (std::find(ip.begin(), ip.end(), net.prio[e]) == ip.end()) ip.push_back(net.prio[e]);
        } else {
            anytimed = true;
        }
    }
    if (anyinh) net.ignored.push_back("inhibitor");
    if (anyimm && anytimed) net.ignored.push_back("immediate");
    if (ip.size() > 1) net.ignored.push_back("priority");
    return net;
}
 
inline long long spn_marking_integer(double x, const std::string& caller) {
    const double r = std::floor(x + 0.5);
    if (std::fabs(x - r) > 1e-9 || r < 0)
        throw InputError(ErrorKind{"Argument"},
                         caller + ": an initial marking must be a non-negative integer");
    return static_cast<long long>(r);
}
 
/** init if given, else the declared Place state if any, else the reference station populations. */
template <class T>
std::vector<long long> spn_initmarking(const qn::NetworkStruct<T>& sn, const SpnNet& net,
                                       const std::vector<double>& init, const std::string& caller) {
    const std::size_t n = net.n, R = net.R, P = net.places.size();
    std::vector<long long> m0(n, 0);
    if (!init.empty()) {
        if (init.size() != n)
            throw InputError(ErrorKind{"Argument"},
                             caller + ": init must hold one token count per place level");
        for (std::size_t i = 0; i < n; ++i) m0[i] = spn_marking_integer(init[i], caller);
        return m0;
    }
    bool anyset = false;
    std::vector<double> dec(n, 0.0);
    for (std::size_t pp = 0; pp < P; ++pp) {
        const typename std::map<std::size_t, std::vector<T>>::const_iterator it =
            sn.initmarking.find(net.places[pp]);
        if (it == sn.initmarking.end()) continue;
        for (std::size_t k = 0; k < R && k < it->second.size(); ++k) {
            dec[pp * R + k] = num_traits<T>::to_double(it->second[k]);
            anyset = anyset || dec[pp * R + k] > 0;
        }
    }
    if (anyset) {
        for (std::size_t i = 0; i < n; ++i) m0[i] = spn_marking_integer(dec[i], caller);
        return m0;
    }
    for (std::size_t k = 0; k < R; ++k) {
        const std::size_t rs = sn.classes[k].refstat;
        if (rs < 1 || rs > sn.station_to_node.size()) continue;
        const std::size_t ref_node = sn.station_to_node[rs - 1];
        for (std::size_t pp = 0; pp < P; ++pp)
            if (net.places[pp] == ref_node)
                m0[pp * R + k] += spn_marking_integer(sn.classes[k].population, caller);
    }
    return m0;
}
 
/**
 * Minimal-support non-negative integer rows x with x' A = 0 (Farkas on
 * [A | I], gcd-normalised, support-minimality filter at every step).
 */
inline std::vector<std::vector<long long>> spn_farkas(const std::vector<std::vector<long long>>& A,
                                                      std::size_t nc) {
    const std::size_t nr = A.size();
    const std::size_t w = nc + nr;
    std::vector<std::vector<long long>> rows(nr, std::vector<long long>(w, 0));
    for (std::size_t p = 0; p < nr; ++p) {
        for (std::size_t c = 0; c < nc; ++c) rows[p][c] = A[p][c];
        rows[p][nc + p] = 1;
    }
    for (std::size_t c = 0; c < nc; ++c) {
        std::vector<std::vector<long long>> nxt;
        for (std::size_t r = 0; r < rows.size(); ++r)
            if (rows[r][c] == 0) nxt.push_back(rows[r]);
        for (std::size_t a = 0; a < rows.size(); ++a) {
            if (rows[a][c] <= 0) continue;
            for (std::size_t b = 0; b < rows.size(); ++b) {
                if (rows[b][c] >= 0) continue;
                const long long pa = rows[a][c], nb = -rows[b][c];
                const long long d = spn_gcd(pa, nb);
                std::vector<long long> combo(w, 0);
                long long g = 0;
                for (std::size_t k = 0; k < w; ++k) {
                    combo[k] = (nb / d) * rows[a][k] + (pa / d) * rows[b][k];
                    g = spn_gcd(g, combo[k]);
                }
                if (g > 1)
                    for (std::size_t k = 0; k < w; ++k) combo[k] /= g;
                bool nz = false;
                for (std::size_t k = 0; k < nr; ++k) nz = nz || combo[nc + k] != 0;
                if (nz) nxt.push_back(combo);
            }
        }
        std::vector<bool> keep(nxt.size(), true);
        for (std::size_t r = 0; r < nxt.size(); ++r)
            for (std::size_t s = 0; s < nxt.size(); ++s) {
                if (s == r || !keep[r]) continue;
                if (!spn_support_subset(nxt[s], nxt[r], nc, nr)) continue;
                if (!spn_support_equal(nxt[s], nxt[r], nc, nr) || s < r) keep[r] = false;
            }
        rows.clear();
        for (std::size_t r = 0; r < nxt.size(); ++r)
            if (keep[r]) rows.push_back(nxt[r]);
    }
    std::vector<std::vector<long long>> X;
    for (std::size_t r = 0; r < rows.size(); ++r) {
        std::vector<long long> x(rows[r].begin() + static_cast<std::ptrdiff_t>(nc), rows[r].end());
        bool nonneg = true;
        for (std::size_t k = 0; k < nr; ++k) nonneg = nonneg && x[k] >= 0;
        if (nonneg) X.push_back(x);
    }
    return X;
}
 
/** Largest siphon inside the level set X (drop every level fed by an event taking nothing from X). */
inline std::vector<bool> spn_maxsiphon(const std::vector<std::vector<long long>>& pre,
                                       const std::vector<std::vector<long long>>& post,
                                       std::vector<bool> X) {
    const std::size_t n = pre.size();
    const std::size_t E = n == 0 ? 0 : pre[0].size();
    for (;;) {
        bool any = false;
        for (std::size_t l = 0; l < n; ++l) any = any || X[l];
        if (!any) return X;
        std::vector<bool> fr(E, true);
        for (std::size_t e = 0; e < E; ++e)
            for (std::size_t l = 0; l < n && fr[e]; ++l)
                if (X[l] && pre[l][e] > 0) fr[e] = false;
        bool drop = false;
        for (std::size_t l = 0; l < n; ++l) {
            if (!X[l]) continue;
            for (std::size_t e = 0; e < E; ++e)
                if (fr[e] && post[l][e] > 0) {
                    X[l] = false;
                    drop = true;
                    break;
                }
        }
        if (!drop) return X;
    }
}
 
inline bool spn_isminimal(const std::vector<std::vector<long long>>& pre,
                          const std::vector<std::vector<long long>>& post,
                          const std::vector<bool>& S) {
    for (std::size_t p = 0; p < S.size(); ++p) {
        if (!S[p]) continue;
        std::vector<bool> X = S;
        X[p] = false;
        const std::vector<bool> Y = spn_maxsiphon(pre, post, X);
        for (std::size_t l = 0; l < Y.size(); ++l)
            if (Y[l]) return false;
    }
    return true;
}
 
/**
 * All minimal siphons, by the branch and bound documented in spn_siphons.h,
 * sorted by size, then lexicographically by their level indices.
 */
inline std::vector<std::vector<bool>> spn_minsiphons(const std::vector<std::vector<long long>>& pre,
                                                     const std::vector<std::vector<long long>>& post) {
    const std::size_t n = pre.size();
    std::vector<std::vector<bool>> found;
    std::vector<std::vector<bool>> stI(1, std::vector<bool>(n, false));
    std::vector<std::vector<bool>> stO(1, std::vector<bool>(n, false));
    while (!stI.empty()) {
        const std::vector<bool> I = stI.back();
        const std::vector<bool> O = stO.back();
        stI.pop_back();
        stO.pop_back();
        std::vector<bool> notO(n);
        for (std::size_t l = 0; l < n; ++l) notO[l] = !O[l];
        const std::vector<bool> S = spn_maxsiphon(pre, post, notO);
        bool anyS = false, miss = false;
        for (std::size_t l = 0; l < n; ++l) {
            anyS = anyS || S[l];
            miss = miss || (I[l] && !S[l]);
        }
        if (!anyS || miss) continue;
        bool inside = false;
        for (std::size_t f = 0; f < found.size() && !inside; ++f) {
            bool sub = true, nz = false;
            for (std::size_t l = 0; l < n; ++l) {
                nz = nz || found[f][l];
                if (found[f][l] && !I[l]) sub = false;
            }
            inside = sub && nz;
        }
        if (inside) continue;
        if (spn_isminimal(pre, post, S)) {
            found.push_back(S);
            continue;
        }
        std::size_t q = n;
        for (std::size_t l = 0; l < n; ++l)
            if (S[l] && !I[l]) {
                q = l;
                break;
            }
        if (q == n) continue;
        std::vector<bool> Oq = O, Iq = I;
        Oq[q] = true;
        Iq[q] = true;
        stI.push_back(Iq);
        stO.push_back(O);
        stI.push_back(I);
        stO.push_back(Oq);
    }
    std::stable_sort(found.begin(), found.end(),
                     [](const std::vector<bool>& a, const std::vector<bool>& b) {
                         std::vector<std::size_t> ia, ib;
                         for (std::size_t l = 0; l < a.size(); ++l) {
                             if (a[l]) ia.push_back(l);
                             if (b[l]) ib.push_back(l);
                         }
                         if (ia.size() != ib.size()) return ia.size() < ib.size();
                         return ia < ib;
                     });
    return found;
}
 
/** Rows of logical sets as 0/1 integers. */
inline std::vector<std::vector<int>> spn_sets_to_int(const std::vector<std::vector<bool>>& L) {
    std::vector<std::vector<int>> out(L.size());
    for (std::size_t i = 0; i < L.size(); ++i) {
        out[i].resize(L[i].size());
        for (std::size_t l = 0; l < L[i].size(); ++l) out[i][l] = L[i][l] ? 1 : 0;
    }
    return out;
}
 
/** The transpose of an nr x nc integer matrix. */
inline std::vector<std::vector<long long>> spn_transpose(const std::vector<std::vector<long long>>& A,
                                                         std::size_t nc) {
    std::vector<std::vector<long long>> B(nc, std::vector<long long>(A.size(), 0));
    for (std::size_t i = 0; i < A.size(); ++i)
        for (std::size_t j = 0; j < nc; ++j) B[j][i] = A[i][j];
    return B;
}
 
}  // namespace detail
}  // namespace spn
}  // namespace line
 
#endif  // LINE_API_SPN_SPN_NET_H