workflow_net.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_WFN_WORKFLOW_NET_H
#define LINE_API_WFN_WORKFLOW_NET_H
 
/**
 * @file
 * @ingroup api_wfn
 * A Workflow Net (WF-net) model for performance analysis.
 *
 * A WF-net (van der Aalst, 1997/98) is a Petri net with exactly one source
 * place i (no incoming arcs), exactly one sink place o (no outgoing arcs), and
 * every place and transition on a directed path from i to o. A case is one
 * token put in i; it is complete when one token reaches o.
 *
 * WorkflowNet has no solver of its own. It checks that the net is a WF-net and
 * that it is sound, and then either analyses one case exactly or compiles the
 * net into an ordinary Network:
 *
 *   validate()                 structural WF-net conditions
 *   check_soundness(k)         k-soundness on the reachability graph
 *   get_cycle_time()           exact phase-type law of the duration of one case
 *                              without contention (any sound net, PH timing)
 *   get_cycle_time_analysis()  moments, CDF, percentiles, expected firings per
 *                              transition and expected token-time per place
 *   get_resource_loads(lambda) per-resource load of an open stream of cases
 *   to_workflow()              equivalent Workflow, for block-structured nets
 *   to_network("forkjoin")     fork-join queueing network, block-structured nets
 *   to_network("spn")          stochastic Petri net with a case-counter place,
 *                              any sound net
 *   get_case_metrics(solver)   case throughput, WIP, cycle time and resource
 *                              utilisation from a solver of a compiled network
 *
 * TIMING. A timed transition fires after a phase-type delay (Exp, Erlang,
 * HyperExp, Coxian, Cox2, APH, PH) under race semantics with enabling memory: a
 * transition disabled by a conflicting firing loses its progress. An immediate
 * transition fires in zero time; among the enabled immediate transitions only
 * those of the highest priority fire, chosen in proportion to their weights.
 *
 * SEMANTICS OF THE TWO NETWORK TARGETS WITH MANY CASES. The fork-join target
 * keeps case identity, since LINE joins only synchronise siblings of the same
 * job. The SPN target uses uncoloured tokens, so an AND-join may pair tokens of
 * different cases; its case-counter place gives the mean number of cases in the
 * system and, by Little's law, the exact mean cycle time under that semantics.
 * Percentiles are not available from it.
 *
 *     wfn::WorkflowNet net("SimpleWFN");
 *     net.add_place("i"); net.add_place("p"); net.add_place("o");
 *     net.add_transition("Task1", Exp(1.0));
 *     net.add_transition("Task2", Exp(2.0));
 *     net.add_arc("i", "Task1"); net.add_arc("Task1", "p");
 *     net.add_arc("p", "Task2"); net.add_arc("Task2", "o");
 *     const wfn::SoundnessReport rep = net.check_soundness();
 *     const Dist ph = net.get_cycle_time();   // hypoexponential, mean 1.5
 *
 * References:
 *   W.M.P. van der Aalst, "The Application of Petri Nets to Workflow
 *   Management", J. Circuits, Systems and Computers 8(1), 1998.
 *
 * MATLAB twin (reference): matlab/src/api/wfn/WorkflowNet.m
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <deque>
#include <limits>
#include <map>
#include <memory>
#include <string>
#include <utility>
#include <vector>
 
#include "line/api/wfn/soundness_report.h"
#include "line/api/wfn/wfn_arc.h"
#include "line/api/wfn/wfn_place.h"
#include "line/api/wfn/wfn_resource.h"
#include "line/api/wfn/wfn_transition.h"
#include "line/io/pnml.h"
#include "line/lang/distributions.h"
#include "line/lang/lang_types.h"
#include "line/lang/workflow/workflow.h"
#include "line/solvers/solver.h"
#include "line/util/error.h"
#include "line/util/expm.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace wfn {
 
namespace detail {
 
/** `sprintf('%g', x)`. */
inline std::string num_g(double x) {
    char buf[64];
    std::snprintf(buf, sizeof buf, "%g", x);
    return buf;
}
 
/** `sprintf('%d', x)` for an integral double, which is how MATLAB prints a count held in a double. */
inline std::string num_d(double x) {
    if (std::isinf(x)) return x > 0 ? "Inf" : "-Inf";
    char buf[64];
    std::snprintf(buf, sizeof buf, "%.0f", x);
    return buf;
}
 
/** An LU factorisation reused across right-hand sides. */
struct Factor {
    Matrix<double> lu;
    std::vector<std::size_t> piv;
    explicit Factor(const Matrix<double>& A) : lu(A) { piv = lu_factor(lu); }
    std::vector<double> solve(std::vector<double> b) const {
        lu_solve(lu, piv, b);
        return b;
    }
};
 
inline Matrix<double> transpose(const Matrix<double>& A) {
    Matrix<double> B(A.cols(), A.rows(), 0.0);
    for (std::size_t i = 0; i < A.rows(); ++i)
        for (std::size_t j = 0; j < A.cols(); ++j) B(j, i) = A(i, j);
    return B;
}
 
/** P(duration <= t) of the possibly defective PH law (alpha, Q); the defect is an atom at zero. */
inline double ph_cdf(const std::vector<double>& alpha, const Matrix<double>& Q, double t) {
    double mass = 0.0;
    for (double a : alpha) mass += a;
    if (t <= 0) return 1.0 - mass;
    const Matrix<double> E = expm(Q, t);
    double surv = 0.0;
    for (std::size_t i = 0; i < alpha.size(); ++i)
        for (std::size_t j = 0; j < alpha.size(); ++j) surv += alpha[i] * E(i, j);
    return 1.0 - surv;
}
 
/** The q-quantile by bisection, bracketed by doubling from the mean MU. */
inline double ph_quantile(const std::vector<double>& alpha, const Matrix<double>& Q, double q, double mu) {
    if (!(q > 0 && q < 1))
        throw InputError(ErrorKind("Argument"), "Percentile levels must lie in (0,1).");
    if (ph_cdf(alpha, Q, 0) >= q) return 0.0;
    double hi = mu;
    while (ph_cdf(alpha, Q, hi) < q) hi = 2 * hi;
    double lo = 0.0;
    for (int it = 1; it <= 200; ++it) {
        const double mid = (lo + hi) / 2;
        if (ph_cdf(alpha, Q, mid) < q)
            lo = mid;
        else
            hi = mid;
        if (hi - lo <= 1e-12 * std::max(1.0, hi)) break;
    }
    return (lo + hi) / 2;
}
 
/** A process-tree node of a block-structured net, as `reduceToTree` builds it. */
struct TreeNode {
    enum Type { EPS, TASK, SEQ, AND, XOR, LOOP };
    Type type = EPS;
    std::size_t t = 0;              ///< transition index of a TASK
    std::vector<TreeNode> kids;     ///< SEQ/AND/XOR children; LOOP holds {do, redo}
    std::vector<double> probs;      ///< XOR branch probabilities
    double q = 0.0;                 ///< LOOP repeat probability
};
 
inline TreeNode tree_eps() { return TreeNode(); }
 
inline TreeNode tree_collapse(TreeNode::Type type, std::vector<TreeNode> kids) {
    if (kids.empty()) return tree_eps();
    if (kids.size() == 1) return kids[0];
    TreeNode n;
    n.type = type;
    n.kids = std::move(kids);
    return n;
}
 
/** Flatten nested blocks of the same TYPE and drop the zero-time parts. */
inline TreeNode tree_flat(TreeNode::Type type, const std::vector<TreeNode>& parts) {
    std::vector<TreeNode> kids;
    for (const TreeNode& p : parts) {
        if (p.type == TreeNode::EPS) continue;
        if (p.type == type)
            kids.insert(kids.end(), p.kids.begin(), p.kids.end());
        else
            kids.push_back(p);
    }
    return tree_collapse(type, kids);
}
 
inline TreeNode tree_seq(const std::vector<TreeNode>& parts) { return tree_flat(TreeNode::SEQ, parts); }
inline TreeNode tree_and(const std::vector<TreeNode>& parts) { return tree_flat(TreeNode::AND, parts); }
 
inline TreeNode tree_xor(const std::vector<TreeNode>& parts, const std::vector<double>& probs) {
    TreeNode n;
    n.type = TreeNode::XOR;
    for (std::size_t j = 0; j < parts.size(); ++j) {
        const TreeNode& p = parts[j];
        if (p.type == TreeNode::XOR) {
            n.kids.insert(n.kids.end(), p.kids.begin(), p.kids.end());
            for (double pr : p.probs) n.probs.push_back(probs[j] * pr);
        } else {
            n.kids.push_back(p);
            n.probs.push_back(probs[j]);
        }
    }
    bool all_eps = true;
    for (const TreeNode& k : n.kids) all_eps = all_eps && k.type == TreeNode::EPS;
    return all_eps ? tree_eps() : n;
}
 
inline TreeNode tree_loop(const TreeNode& do_part, const TreeNode& redo_part, double q) {
    TreeNode n;
    n.type = TreeNode::LOOP;
    n.kids.push_back(do_part);
    n.kids.push_back(redo_part);
    n.q = q;
    return n;
}
 
/** One row of an AvgTable; a pair the solver did not report is zero, as in MATLAB's getAvg. */
inline double table_cell(const AvgTable& tab, const std::vector<double>& col, const std::string& station,
                         const std::string& cls) {
    for (std::size_t i = 0; i < tab.Station.size(); ++i)
        if (tab.Station[i] == station && tab.JobClass[i] == cls && i < col.size())
            return std::isnan(col[i]) ? 0.0 : col[i];
    return 0.0;
}
 
/** `sum(M(station, :))` over the classes. */
inline double table_sum(const AvgTable& tab, const std::vector<double>& col, const std::string& station) {
    double s = 0.0;
    for (std::size_t i = 0; i < tab.Station.size(); ++i)
        if (tab.Station[i] == station && i < col.size() && !std::isnan(col[i])) s += col[i];
    return s;
}
 
}  // namespace detail
 
/** The full single-case analysis of `get_cycle_time_analysis`. */
struct CycleTimeAnalysis {
    /** PH law of the case duration. */
    lang::Distrib<double> ph;
    double mean = 0.0;
    double variance = 0.0;
    double scv = 0.0;
    /** The first three raw moments. */
    std::vector<double> moments;
    /** Probability that a case completes in zero time. */
    double zero_time_prob = 0.0;
    /** The (alpha, Q) pair behind `ph`; alpha may be defective by `zero_time_prob`. */
    std::vector<double> alpha;
    Matrix<double> Q;
    /** Number of tangible transient states. */
    std::size_t num_states = 0;
 
    struct Firings {
        std::vector<std::string> Transition;
        std::vector<double> ExpectedFirings;
    } firings;
    struct PlaceTime {
        std::vector<std::string> Place;
        std::vector<double> ExpectedTokenTime;
    } place_time;
 
    /** P(duration <= t). */
    double cdf(double t) const { return detail::ph_cdf(alpha, Q, t); }
    /** The q-quantile, q in (0,1). */
    double percentile(double q) const { return detail::ph_quantile(alpha, Q, q, mean); }
};
 
/** The table of `get_resource_loads`. */
struct ResourceLoads {
    std::vector<std::string> Resource;
    std::vector<double> Servers, DemandPerCase, Load;
    std::vector<bool> Stable;
    std::size_t size() const { return Resource.size(); }
};
 
/** Case-level metrics of `get_case_metrics`. */
struct CaseMetrics {
    std::string target;   ///< "forkjoin" or "spn"
    std::string mode;     ///< "open" or "closed"
    double throughput = 0.0;
    double wip = 0.0;
    double cycle_time = 0.0;
    struct Resources {
        std::vector<std::string> Resource;
        std::vector<double> Servers, Utilization;
    } resources;
    /** RespT at its station for the fork-join target, InService tokens for the SPN one. */
    struct Tasks {
        std::vector<std::string> Task;
        std::vector<double> RespT, InService;
    } tasks;
    /** Mean tokens per place, SPN target only. */
    struct Places {
        std::vector<std::string> Place;
        std::vector<double> QLen;
    } places;
};
 
/** A place, a transition, or the name of either, where an argument may be any of them. */
class ElementRef {
public:
    ElementRef(const WFNPlace& p) : place_(&p) {}
    ElementRef(const WFNTransition& t) : transition_(&t) {}
    ElementRef(const std::string& n) : name_(n), by_name_(true) {}
    ElementRef(const char* n) : name_(n), by_name_(true) {}
 
private:
    friend class WorkflowNet;
    const WFNPlace* place_ = nullptr;
    const WFNTransition* transition_ = nullptr;
    std::string name_;
    bool by_name_ = false;
};
 
class WorkflowNet {
public:
    typedef lang::Distrib<double> Dist;
 
    /** Names prefixed so that compiled networks never collide with user names. */
    static constexpr const char* ReservedPrefix = "wfn_";
 
    explicit WorkflowNet(const std::string& name) : name_(name) {}
 
    WorkflowNet(WorkflowNet&&) = default;
    WorkflowNet& operator=(WorkflowNet&&) = default;
    WorkflowNet(const WorkflowNet&) = delete;
    WorkflowNet& operator=(const WorkflowNet&) = delete;
 
    const std::string& name() const { return name_; }
    std::size_t num_places() const { return places_.size(); }
    std::size_t num_transitions() const { return transitions_.size(); }
    const WFNPlace& place(std::size_t i) const { return *places_.at(i); }
    WFNTransition& transition(std::size_t i) { return *transitions_.at(i); }
    const WFNTransition& transition(std::size_t i) const { return *transitions_.at(i); }
    const std::vector<WFNArc>& arcs() const { return arcs_; }
    std::size_t num_resources() const { return resources_.size(); }
    const WFNResource& resource(std::size_t r) const { return *resources_.at(r); }
    /** The source place, set by `validate` when not declared; nullptr before. */
    const WFNPlace* source_place() const { return source_ < 0 ? nullptr : places_[source_].get(); }
    /** The sink place, set by `validate` when not declared; nullptr before. */
    const WFNPlace* sink_place() const { return sink_ < 0 ? nullptr : places_[sink_].get(); }
    /** The interarrival law of cases (open analysis), or nullptr. */
    const Dist* arrival_process() const { return has_arrival_ ? &arrival_ : nullptr; }
    long closed_cases() const { return closed_cases_; }
    double state_cap() const { return state_cap_; }
 
    // ------------------------------------------------------------------
    // Construction
    // ------------------------------------------------------------------
 
    /**
     * `addPlace(name, initialMarking)`. The initial marking is informational:
     * the analyses always start from k tokens in the source place.
     */
    WFNPlace& add_place(const std::string& name, double initial_marking = 0.0) {
        check_new_name(name);
        places_.emplace_back(new WFNPlace(name, places_.size(), initial_marking));
        place_map_[name] = places_.size() - 1;
        return *places_.back();
    }
 
    /**
     * `addTransition(name, distribution, [], servers, priority, weight)`: the
     * timing follows the law, IMMEDIATE for an Immediate law and TIMED for any
     * other. SERVERS is the number of concurrent firings (default Inf, i.e. no
     * contention unless a resource is assigned); PRIORITY and WEIGHT rank and
     * share the enabled immediate transitions.
     */
    WFNTransition& add_transition(const std::string& name, const Dist& distribution = Exp(1.0),
                                  double servers = std::numeric_limits<double>::infinity(),
                                  double priority = 1.0, double weight = 1.0) {
        const lang::TimingStrategy timing = distribution.type == lang::ProcessType::IMMEDIATE
                                                ? lang::TimingStrategy::IMMEDIATE
                                                : lang::TimingStrategy::TIMED;
        return add_transition(name, distribution, timing, servers, priority, weight);
    }
 
    /** `addTransition(name, distribution, timing, servers, priority, weight)`. */
    WFNTransition& add_transition(const std::string& name, const Dist& distribution,
                                  lang::TimingStrategy timing,
                                  double servers = std::numeric_limits<double>::infinity(),
                                  double priority = 1.0, double weight = 1.0) {
        Dist d = distribution;
        if (timing == lang::TimingStrategy::IMMEDIATE) {
            d = Immediate();
        } else if (distribution.type == lang::ProcessType::IMMEDIATE) {
            throw InputError(ErrorKind("Argument"),
                             "Timed transition \"" + name +
                                 "\" cannot have an Immediate distribution; declare it TimingStrategy.IMMEDIATE.");
        }
        if (!(servers >= 1))
            throw InputError(ErrorKind("Argument"), "Number of servers must be at least 1.");
        if (!(weight > 0) || !(priority > 0))
            throw InputError(ErrorKind("Argument"), "Firing weight and priority must be positive.");
        check_new_name(name);
        transitions_.emplace_back(
            new WFNTransition(name, transitions_.size(), d, timing, servers, priority, weight));
        transition_map_[name] = transitions_.size() - 1;
        return *transitions_.back();
    }
 
    /**
     * `addArc(source, target, weight)`: from a place to a transition or from a
     * transition to a place, given as elements of this net or by name. WEIGHT
     * is a positive integer.
     */
    const WFNArc& add_arc(const ElementRef& source, const ElementRef& target, double weight = 1.0) {
        const ElementRef src = resolve(source);
        const ElementRef tgt = resolve(target);
        check_owned(src);
        check_owned(tgt);
        if ((src.place_ && tgt.place_) || (src.transition_ && tgt.transition_))
            throw InputError(ErrorKind("Model"),
                             "Arcs must connect a place to a transition or vice versa, not two places or two "
                             "transitions.");
        if (!(weight >= 1) || weight != std::round(weight))
            throw InputError(ErrorKind("Model"), "Arc " + element_name(src) + " -> " + element_name(tgt) +
                                                     " has weight " + detail::num_g(weight) +
                                                     "; arc weights must be positive integers.");
        if (src.place_)
            arcs_.push_back(WFNArc(true, src.place_, tgt.transition_, weight));
        else
            arcs_.push_back(WFNArc(false, tgt.place_, src.transition_, weight));
        return arcs_.back();
    }
 
    /** `setSource(place)`: declare the source place (checked by validate). */
    void set_source(const ElementRef& place) { source_ = declared_place(place, "Source"); }
 
    /** `setSink(place)`: declare the sink place (checked by validate). */
    void set_sink(const ElementRef& place) { sink_ = declared_place(place, "Sink"); }
 
    /**
     * `setResource(transition, name, servers, scheduling)`: assign a resource,
     * shared by every task that names it and by every case. An existing NAME is
     * reused; SERVERS and SCHEDULING, when given, update it.
     */
    WFNResource& set_resource(const ElementRef& transition, const std::string& name) {
        return set_resource_impl(transition, name, nullptr, nullptr);
    }
    WFNResource& set_resource(const ElementRef& transition, const std::string& name, double servers) {
        return set_resource_impl(transition, name, &servers, nullptr);
    }
    WFNResource& set_resource(const ElementRef& transition, const std::string& name, double servers,
                              lang::SchedStrategy scheduling) {
        return set_resource_impl(transition, name, &servers, &scheduling);
    }
 
    /** `setArrivalProcess(dist)`: any interarrival distribution for open analysis. */
    void set_arrival_process(const Dist& process) {
        arrival_ = process;
        has_arrival_ = true;
    }
    /** `setArrivalProcess('none')`, or a type that needs a parameter, which is then refused. */
    void set_arrival_process(const std::string& type) { set_arrival_typed(type, nullptr, nullptr); }
    /** `setArrivalProcess('poisson', rate)`. */
    void set_arrival_process(const std::string& type, double rate) { set_arrival_typed(type, &rate, nullptr); }
    /** `setArrivalProcess('map', MAP)`. */
    void set_arrival_process(const std::string& type, const Dist& map) { set_arrival_typed(type, nullptr, &map); }
 
    /** `setClosedCases(k)`: k concurrent cases (closed analysis); 0 = open. */
    void set_closed_cases(double k) {
        if (!(k >= 0) || k != std::round(k) || std::isinf(k))
            throw InputError(ErrorKind("Argument"), "The number of closed cases must be a non-negative integer.");
        closed_cases_ = static_cast<long>(k);
    }
 
    /** `setStateCap(cap)`: the largest state space explored by the analyses. */
    void set_state_cap(double cap) {
        if (!(cap >= 1)) throw InputError(ErrorKind("Argument"), "The state cap must be a positive number.");
        state_cap_ = cap;
    }
 
    /** `getPlace(name)`. */
    const WFNPlace& get_place(const std::string& name) const {
        const auto it = place_map_.find(name);
        if (it == place_map_.end()) throw InputError(ErrorKind("Model"), "Place \"" + name + "\" not found.");
        return *places_[it->second];
    }
 
    /** `getTransition(name)`. */
    WFNTransition& get_transition(const std::string& name) {
        const auto it = transition_map_.find(name);
        if (it == transition_map_.end())
            throw InputError(ErrorKind("Model"), "Transition \"" + name + "\" not found.");
        return *transitions_[it->second];
    }
 
    /** `getIncidence()`: input and output arc weights, places x transitions. */
    std::pair<Matrix<double>, Matrix<double>> get_incidence() const {
        const Inc inc = incidence();
        Matrix<double> Pre(inc.nP, inc.nT, 0.0), Post(inc.nP, inc.nT, 0.0);
        for (std::size_t p = 0; p < inc.nP; ++p)
            for (std::size_t t = 0; t < inc.nT; ++t) {
                Pre(p, t) = static_cast<double>(inc.pre[p][t]);
                Post(p, t) = static_cast<double>(inc.post[p][t]);
            }
        return std::make_pair(Pre, Post);
    }
 
    // ------------------------------------------------------------------
    // Structure and soundness
    // ------------------------------------------------------------------
 
    /**
     * `validate()`: exactly one place without incoming arcs (the source i),
     * exactly one without outgoing arcs (the sink o), and every place and
     * transition on a directed path from i to o. A declared source or sink must
     * be that place. On success the source and sink are set; on failure a
     * LINE:input:Model error names the offending nodes.
     */
    void validate() {
        if (places_.empty())
            throw InputError(ErrorKind("Model"), "A workflow net must have at least one place.");
        if (transitions_.empty())
            throw InputError(ErrorKind("Model"), "A workflow net must have at least one transition.");
        const Inc inc = incidence();
        const std::size_t nP = inc.nP, nT = inc.nT;
        for (std::size_t t = 0; t < nT; ++t) {
            bool in = false, out = false;
            for (std::size_t p = 0; p < nP; ++p) {
                in = in || inc.pre[p][t] != 0;
                out = out || inc.post[p][t] != 0;
            }
            if (!in)
                throw InputError(ErrorKind("Model"), "Transition \"" + transitions_[t]->name() +
                                                         "\" has no input place, so it is not on a path from the "
                                                         "source.");
            if (!out)
                throw InputError(ErrorKind("Model"), "Transition \"" + transitions_[t]->name() +
                                                         "\" has no output place, so it is not on a path to the "
                                                         "sink.");
        }
        std::vector<std::size_t> src_cand, snk_cand;
        for (std::size_t p = 0; p < nP; ++p) {
            long in = 0, out = 0;
            for (std::size_t t = 0; t < nT; ++t) {
                in += inc.post[p][t];
                out += inc.pre[p][t];
            }
            if (in == 0) src_cand.push_back(p);
            if (out == 0) snk_cand.push_back(p);
        }
        if (src_cand.size() != 1)
            throw InputError(ErrorKind("Model"),
                             "A WF-net has exactly one source place (a place without incoming arcs); found " +
                                 std::to_string(src_cand.size()) + ": " + place_list(src_cand) + ".");
        if (snk_cand.size() != 1)
            throw InputError(ErrorKind("Model"),
                             "A WF-net has exactly one sink place (a place without outgoing arcs); found " +
                                 std::to_string(snk_cand.size()) + ": " + place_list(snk_cand) + ".");
        const std::size_t s = src_cand[0], o = snk_cand[0];
        if (source_ >= 0 && static_cast<std::size_t>(source_) != s)
            throw InputError(ErrorKind("Model"), "Place \"" + places_[source_]->name() +
                                                     "\" was declared the source, but \"" + places_[s]->name() +
                                                     "\" is the only place without incoming arcs.");
        if (sink_ >= 0 && static_cast<std::size_t>(sink_) != o)
            throw InputError(ErrorKind("Model"), "Place \"" + places_[sink_]->name() +
                                                     "\" was declared the sink, but \"" + places_[o]->name() +
                                                     "\" is the only place without outgoing arcs.");
        // bipartite graph over places 0..nP-1 then transitions nP..nP+nT-1
        const std::size_t n = nP + nT;
        std::vector<std::vector<std::size_t>> fwd_adj(n), bwd_adj(n);
        for (std::size_t p = 0; p < nP; ++p)
            for (std::size_t t = 0; t < nT; ++t) {
                if (inc.pre[p][t] > 0) {
                    fwd_adj[p].push_back(nP + t);
                    bwd_adj[nP + t].push_back(p);
                }
                if (inc.post[p][t] > 0) {
                    fwd_adj[nP + t].push_back(p);
                    bwd_adj[p].push_back(nP + t);
                }
            }
        const std::vector<bool> fwd = bfs(fwd_adj, std::vector<std::size_t>(1, s));
        const std::vector<bool> bwd = bfs(bwd_adj, std::vector<std::size_t>(1, o));
        std::string off;
        for (std::size_t v = 0; v < n; ++v) {
            if (fwd[v] && bwd[v]) continue;
            if (!off.empty()) off += ", ";
            off += "\"" + (v < nP ? places_[v]->name() : transitions_[v - nP]->name()) + "\"";
        }
        if (!off.empty())
            throw InputError(ErrorKind("Model"), "Not a WF-net: " + off + " not on a path from source \"" +
                                                     places_[s]->name() + "\" to sink \"" + places_[o]->name() +
                                                     "\".");
        source_ = static_cast<long>(s);
        sink_ = static_cast<long>(o);
    }
 
    /**
     * `checkSoundness(k, stateCap)`: decide k-soundness on the reachability graph.
     *
     * From the marking with K tokens in the source the net is k-sound when (1)
     * the final marking, K tokens in the sink, is reachable from every reachable
     * marking (option to complete), (2) no reachable marking puts K tokens in
     * the sink together with other tokens (proper completion) and (3) every
     * transition fires in some reachable marking (no dead transitions). The
     * untimed net is explored, as soundness is a property of the net, not of its
     * timing. A marking that strictly covers one of its ancestors proves the net
     * unbounded, hence unsound. If STATECAP markings are explored first, the
     * verdict is "unknown", never "sound".
     */
    SoundnessReport check_soundness(double k = 1) { return check_soundness(k, state_cap_); }
 
    SoundnessReport check_soundness(double k, double state_cap) {
        if (!(k >= 1) || k != std::round(k) || std::isinf(k))
            throw InputError(ErrorKind("Argument"), "k must be a positive integer.");
        validate();
        const Inc inc = incidence();
        const std::size_t nP = inc.nP, nT = inc.nT;
        const long kk = static_cast<long>(k);
        Marking m0(nP, 0), mf(nP, 0);
        m0[source_] = kk;
        mf[sink_] = kk;
 
        SoundnessReport report;
        report.k = kk;
        report.verdict = "unknown";
 
        Reach rg = untimed_reachability(inc, m0, state_cap);
        report.state_space_size = rg.M.size();
        if (rg.status == Reach::UNBOUNDED) {
            report.verdict = "unsound";
            report.unbounded = true;
            report.unbounded_witness = rg.witness_cover;
            report.witnesses.push_back("unbounded: " + marking_string(rg.witness_cover) +
                                       " strictly covers its ancestor " + marking_string(rg.witness_ancestor));
            report.sound = false;
            return report;
        }
        if (rg.status == Reach::CAP) {
            report.verdict = "unknown";
            report.unbounded_or_too_large = true;
            report.witnesses.push_back("state cap of " + detail::num_d(state_cap) +
                                       " markings reached before the reachability graph was complete");
            report.sound = false;
            return report;
        }
 
        const std::size_t n = rg.M.size();
        std::vector<bool> is_final(n, false);
        bool any_final = false;
        for (std::size_t i = 0; i < n; ++i) {
            is_final[i] = rg.M[i] == mf;
            any_final = any_final || is_final[i];
        }
        std::vector<std::size_t> outdeg(n, 0);
        std::vector<std::vector<std::size_t>> radj(n);
        std::vector<bool> fired(nT, false);
        for (const ReachEdge& e : rg.E) {
            ++outdeg[e.from];
            radj[e.to].push_back(e.from);
            fired[e.t] = true;
        }
        for (std::size_t i = 0; i < n; ++i)
            if (outdeg[i] == 0 && !is_final[i]) {
                report.deadlocks.push_back(rg.M[i]);
                report.witnesses.push_back("deadlock: " + marking_string(rg.M[i]) + " enables no transition");
            }
        for (std::size_t i = 0; i < n; ++i)
            if (rg.M[i][sink_] >= kk && !is_final[i]) {
                report.improper_completion.push_back(rg.M[i]);
                report.witnesses.push_back("improper completion: " + marking_string(rg.M[i]));
            }
        std::vector<bool> coreach(n, false);
        if (any_final) {
            std::vector<std::size_t> finals;
            for (std::size_t i = 0; i < n; ++i)
                if (is_final[i]) finals.push_back(i);
            coreach = bfs(radj, finals);
        }
        for (std::size_t i = 0; i < n; ++i)
            if (!coreach[i]) {
                report.no_option_to_complete.push_back(rg.M[i]);
                if (outdeg[i] > 0)
                    report.witnesses.push_back("no option to complete: from " + marking_string(rg.M[i]) +
                                               " the final marking is unreachable");
            }
        for (std::size_t t = 0; t < nT; ++t)
            if (!fired[t]) {
                report.dead_transitions.push_back(transitions_[t]->name());
                report.witnesses.push_back("dead transition: \"" + transitions_[t]->name() + "\" never fires");
            }
        report.sound = any_final && report.deadlocks.empty() && report.improper_completion.empty() &&
                       report.no_option_to_complete.empty() && report.dead_transitions.empty();
        report.verdict = report.sound ? "sound" : "unsound";
        return report;
    }
 
    /** `formatReport(report)`: one-line summary of a soundness report. */
    static std::string format_report(const SoundnessReport& report) {
        if (report.verdict == "sound")
            return std::to_string(report.k) + "-sound, " + std::to_string(report.state_space_size) +
                   " reachable markings";
        std::string str;
        for (std::size_t i = 0; i < report.witnesses.size(); ++i) {
            if (i) str += "; ";
            str += report.witnesses[i];
        }
        return str.empty() ? report.verdict : str;
    }
 
    // ------------------------------------------------------------------
    // Exact single-case analysis
    // ------------------------------------------------------------------
 
    /**
     * `getCycleTime()`: exact case cycle time of one case without contention.
     *
     * The reachability graph of the timed net from one token in the source
     * becomes an absorbing CTMC whose absorbing state is the final marking;
     * immediate transitions are folded out of it, so the case duration is
     * exactly phase-type for any sound net with phase-type timing,
     * series-parallel or not. Requires 1-soundness.
     */
    Dist get_cycle_time() { return get_cycle_time_analysis().ph; }
 
    /** `getCycleTimeAnalysis()`: the full single-case analysis. */
    CycleTimeAnalysis get_cycle_time_analysis() {
        const Chain ch = single_case_chain();
        const std::size_t n = ch.Q.rows();
        double mass = 0.0;
        for (double a : ch.alpha) mass += a;
        if (n == 0 || mass <= lang::GlobalConstants::FineTol)
            throw UnsupportedError(ErrorKind("Structure"),
                                   "Every case of \"" + name_ +
                                       "\" completes in zero time, so its cycle time has no phase-type law.");
        Matrix<double> negQ(n, n, 0.0);
        for (std::size_t i = 0; i < n; ++i)
            for (std::size_t j = 0; j < n; ++j) negQ(i, j) = -ch.Q(i, j);
        // row-vector solves v / (-Q) are (-Q)' x = v'
        const detail::Factor fT(detail::transpose(negQ));
        const std::vector<double> Nrow = fT.solve(ch.alpha);
        std::vector<double> m(3, 0.0), v = ch.alpha;
        double fact = 1.0;
        for (int j = 1; j <= 3; ++j) {
            v = fT.solve(v);
            fact *= j;
            double s = 0.0;
            for (double x : v) s += x;
            m[j - 1] = fact * s;
        }
        CycleTimeAnalysis an;
        an.ph = PH(ch.alpha, ch.Q);
        an.mean = m[0];
        an.moments = m;
        an.variance = m[1] - m[0] * m[0];
        an.scv = an.variance / (m[0] * m[0]);
        an.zero_time_prob = 1.0 - mass;
        an.alpha = ch.alpha;
        an.Q = ch.Q;
        const std::size_t nT = transitions_.size(), nP = places_.size();
        for (std::size_t t = 0; t < nT; ++t) {
            double f = ch.firings0[t];
            for (std::size_t s = 0; s < n; ++s) f += Nrow[s] * (ch.Rfire(s, t) + ch.Gfire(s, t));
            an.firings.Transition.push_back(transitions_[t]->name());
            an.firings.ExpectedFirings.push_back(f);
        }
        for (std::size_t p = 0; p < nP; ++p) {
            double x = 0.0;
            for (std::size_t s = 0; s < n; ++s) x += Nrow[s] * static_cast<double>(ch.Mtang[s][p]);
            an.place_time.Place.push_back(places_[p]->name());
            an.place_time.ExpectedTokenTime.push_back(x);
        }
        an.num_states = n;
        return an;
    }
 
    /**
     * `getResourceLoads(lambda)`: one row per resource and per task with a
     * finite server count, the offered load LAMBDA * sum_t E[firings of t] *
     * E[time of t] divided by the servers, from the expected firings of one
     * case. The open system is stable only if every load is below 1.
     */
    ResourceLoads get_resource_loads(double lambda) {
        const CycleTimeAnalysis an = get_cycle_time_analysis();
        const std::vector<double>& f = an.firings.ExpectedFirings;
        ResourceLoads loads;
        for (std::size_t r = 0; r < resources_.size(); ++r) {
            double d = 0.0;
            for (std::size_t t = 0; t < transitions_.size(); ++t)
                if (transitions_[t]->resource_ == static_cast<long>(r))
                    d += f[t] * transitions_[t]->distribution().mean;
            loads.Resource.push_back(resources_[r]->name());
            loads.Servers.push_back(resources_[r]->servers());
            loads.DemandPerCase.push_back(d);
        }
        for (std::size_t t = 0; t < transitions_.size(); ++t) {
            const WFNTransition& tr = *transitions_[t];
            if (tr.is_timed() && std::isfinite(tr.servers())) {
                loads.Resource.push_back(tr.name());
                loads.Servers.push_back(tr.servers());
                loads.DemandPerCase.push_back(f[t] * tr.distribution().mean);
            }
        }
        for (std::size_t j = 0; j < loads.Resource.size(); ++j) {
            loads.Load.push_back(lambda * loads.DemandPerCase[j] / loads.Servers[j]);
            loads.Stable.push_back(loads.Load[j] < 1);
        }
        return loads;
    }
 
    /** `getResourceLoads()`: LAMBDA is the rate of the arrival process. */
    ResourceLoads get_resource_loads() {
        if (!has_arrival_) throw InputError(ErrorKind("Argument"), "Give LAMBDA or set an arrival process.");
        return get_resource_loads(1.0 / arrival_.mean);
    }
 
    // ------------------------------------------------------------------
    // Block structure and compilation
    // ------------------------------------------------------------------
 
    /**
     * `isBlockStructured()`: true if the net reduces to a process tree. Five
     * rules apply until nothing does: fold a single-input single-output
     * transition between two places into a place, fold a single-input
     * single-output place between two transitions into a transition, merge
     * parallel places (AND block), merge parallel transitions (XOR block,
     * choices resolved by immediate transitions of equal priority, in
     * proportion to their weights) and absorb a self-loop transition (geometric
     * loop). It is block-structured if it reduces to one place.
     */
    bool is_block_structured() {
        validate();
        detail::TreeNode tree;
        std::string reason;
        return reduce_to_tree(tree, reason);
    }
 
    /**
     * `toWorkflow()`: the equivalent Workflow of a block-structured net. Each
     * timed transition is an activity with its law and each block the Workflow
     * precedence of the same meaning. Immediate transitions disappear; where a
     * fork or join needs an anchor the net does not provide, an Immediate
     * activity wfn_aux_<n> is added. Workflow::Loop repeats one body a
     * geometric number of times >= 1, so a loop with work on its return path is
     * refused. Workflow models an Immediate activity at rate
     * GlobalConstants::Immediate, so its to_ph() agrees with get_cycle_time()
     * up to that rate, not to machine precision.
     */
    workflow::Workflow<double> to_workflow() {
        require_sound(1);
        detail::TreeNode tree;
        std::string reason;
        if (!reduce_to_tree(tree, reason))
            throw UnsupportedError(ErrorKind("Structure"),
                                   "toWorkflow requires a block-structured WF-net: " + reason);
        WfState st{workflow::Workflow<double>(name_ + "_workflow"), 0};
        std::string h, tl;
        wf_branch(tree, st, h, tl);
        return st.wf;
    }
 
    /**
     * `toNetwork(target, 'checkSoundness', tf)`: compile the net to a Network.
     *
     * "forkjoin": block-structured nets. AND-split -> Fork, AND-join -> Join,
     * XOR and loops -> probabilistic routing, a task -> Delay (no contention), a
     * Queue of its own (finite servers) or the Queue of its resource. Tasks
     * sharing a resource are told apart by class. Solvers: MVA (fork-join
     * approximations), JMT, LDES.
     *
     * "spn": any sound net, as a stochastic Petri net with a case-counter place.
     * A task on a resource becomes an immediate acquire followed by the timed
     * task, which releases the resource; waiting tasks take a free server in
     * random order, in proportion to their weights, whatever the resource
     * scheduling. Its transitions create and destroy tokens, which SolverSSA and
     * SolverLDES simulate.
     *
     * Open analysis (set_arrival_process) gives a Source and a Sink; closed
     * analysis (set_closed_cases(k)) recirculates k cases through a zero-time
     * restart. The soundness gate (k-soundness, k the number of closed cases or
     * 1) is on by default; pass false to skip it on nets too large to explore.
     */
    std::shared_ptr<Network> to_network(const std::string& target, bool check_soundness = true) {
        if (target.empty())
            throw InputError(ErrorKind("Argument"), "Target must be specified as 'forkjoin' or 'spn'.");
        if (closed_cases_ == 0 && !has_arrival_)
            throw InputError(ErrorKind("Model"), "Set an arrival process (open analysis) or setClosedCases(k) "
                                                 "(closed analysis) before compiling to a Network.");
        if (closed_cases_ > 0 && has_arrival_)
            throw InputError(ErrorKind("Model"),
                             "Both an arrival process and closed cases are set; clear one of them "
                             "(setArrivalProcess('none') or setClosedCases(0)).");
        validate();
        if (check_soundness) require_sound(std::max(1L, closed_cases_));
        std::string tg = target;
        std::transform(tg.begin(), tg.end(), tg.begin(), [](unsigned char c) { return std::tolower(c); });
        if (tg == "forkjoin") return compile_to_fork_join();
        if (tg == "spn") return compile_to_spn();
        throw InputError(ErrorKind("Argument"), "Unknown compilation target: " + target);
    }
 
    /**
     * `getCaseMetrics(solver)`: case-level metrics from a solver built on the
     * network returned by the last to_network call. The fork-join open cycle
     * time follows the process tree; the closed one, and the SPN one, are
     * Little's law on the case population.
     */
    CaseMetrics get_case_metrics(NetworkSolver& solver) {
        if (!compiled_ || &solver.model() != compiled_->network.get())
            throw InputError(ErrorKind("Argument"), "The solver is not built on the network returned by the last "
                                                    "toNetwork call of this WorkflowNet.");
        const Compiled& c = *compiled_;
        const AvgTable& tab = solver.avg_table();
        CaseMetrics metrics;
        metrics.target = c.target;
        metrics.mode = c.mode;
        if (c.target == "forkjoin") {
            double X;
            if (c.mode == "open") {
                X = c.lambda;
                std::size_t jidx = 0;
                metrics.cycle_time = fj_tree_time(c.tree, c, tab, jidx);
                metrics.wip = X * metrics.cycle_time;
            } else {
                X = detail::table_cell(tab, tab.Tput, c.restart, "Case");
                metrics.wip = static_cast<double>(c.k);
                metrics.cycle_time = static_cast<double>(c.k) / X;
            }
            metrics.throughput = X;
            for (const Compiled::Task& task : c.tasks) {
                metrics.tasks.Task.push_back(task.name);
                metrics.tasks.RespT.push_back(detail::table_cell(tab, tab.RespT, task.station, task.cls));
            }
            for (const Compiled::Res& res : c.resources) {
                metrics.resources.Resource.push_back(res.name);
                metrics.resources.Servers.push_back(res.servers);
                metrics.resources.Utilization.push_back(detail::table_sum(tab, tab.Util, res.station));
            }
        } else {
            metrics.wip = detail::table_sum(tab, tab.QLen, c.counter);
            const double X = c.mode == "open" ? c.lambda : detail::table_sum(tab, tab.Tput, c.counter);
            metrics.throughput = X;
            metrics.cycle_time = metrics.wip / X;
            for (const std::string& p : c.places) {
                metrics.places.Place.push_back(p);
                metrics.places.QLen.push_back(detail::table_sum(tab, tab.QLen, p));
            }
            for (const Compiled::Busy& b : c.busy) {
                metrics.tasks.Task.push_back(b.name);
                metrics.tasks.InService.push_back(detail::table_sum(tab, tab.QLen, b.place));
            }
            for (const Compiled::Res& res : c.resources) {
                double busy = 0.0;
                for (std::size_t b : res.busy) busy += detail::table_sum(tab, tab.QLen, c.busy[b].place);
                metrics.resources.Resource.push_back(res.name);
                metrics.resources.Servers.push_back(res.servers);
                metrics.resources.Utilization.push_back(busy / res.servers);
            }
        }
        return metrics;
    }
 
    // ------------------------------------------------------------------
    // PNML and conversion
    // ------------------------------------------------------------------
 
    /**
     * `savePNML(filename)`: write the plain Petri net (no case counter, no
     * resources) with one token in the source, so that from_pnml reads it back
     * with its timing.
     */
    void save_pnml(const std::string& filename) {
        validate();
        std::shared_ptr<Network> model = build_plain_network(1);
        io::pnml_save(model->get_struct(), filename);
    }
 
    /**
     * `WorkflowNet.fromPNML(filename, netName)`: read a net with pnml_load and
     * convert it with from_network; the source and sink are found from the
     * structure. A file from another tool is untimed and reads with every
     * transition exponential at rate 1.
     */
    static WorkflowNet from_pnml(const std::string& filename, const std::string& net_name = std::string()) {
        Network model = io::pnml_load<double>(filename, net_name);
        return from_network(model);
    }
 
    /**
     * `WorkflowNet.fromNetwork(model)`: copy the places, transitions and arcs of
     * a single-class Network of Place and Transition nodes. A transition must
     * have exactly one mode and no inhibitor arcs or firing-rate dependence,
     * since a WF-net transition has one firing rule. The result is validated.
     */
    static WorkflowNet from_network(Network& model) {
        const qn::NetworkStruct<double>& sn = model.get_struct();
        if (sn.classes.size() != 1)
            throw UnsupportedError(ErrorKind("MultiClass"),
                                   "A WF-net has uncoloured tokens; the Network must have exactly one class.");
        WorkflowNet net(sn.name + "_wfn");
        const std::size_t nn = sn.nodes.size();
        for (std::size_t i = 0; i < nn; ++i) {
            const lang::NodeType nt = sn.nodes[i].nodetype;
            if (nt == lang::NodeType::Place) {
                double st = 0.0;
                const auto it = sn.initmarking.find(i + 1);
                if (it != sn.initmarking.end())
                    for (double x : it->second) st += x;
                net.add_place(sn.nodes[i].name, st);
            } else if (nt != lang::NodeType::Transition) {
                throw UnsupportedError(ErrorKind("Node"), "Node \"" + sn.nodes[i].name + "\" is a " +
                                                              lang::node_type_to_text(nt) +
                                                              "; a WF-net holds only places and transitions.");
            }
        }
        for (std::size_t i = 0; i < nn; ++i) {
            if (sn.nodes[i].nodetype != lang::NodeType::Transition) continue;
            const std::string& tname = sn.nodes[i].name;
            const auto it = sn.transparam.find(i + 1);
            const std::size_t nmodes = it == sn.transparam.end() ? 0 : it->second.nmodes;
            if (nmodes != 1)
                throw UnsupportedError(ErrorKind("Transition"), "Transition \"" + tname + "\" has " +
                                                                    std::to_string(nmodes) +
                                                                    " modes; a WF-net transition has exactly one "
                                                                    "firing rule.");
            const qn::TransitionParam<double>& tp = it->second;
            const Matrix<double>& Ein = tp.enabling[0];
            const Matrix<double>& Eout = tp.firing[0];
            const Matrix<double> Inh = tp.inhibiting.empty() ? Matrix<double>() : tp.inhibiting[0];
            for (std::size_t r = 0; r < Inh.rows(); ++r)
                for (std::size_t c = 0; c < Inh.cols(); ++c)
                    if (std::isfinite(Inh(r, c)) && Inh(r, c) > 0)
                        throw UnsupportedError(ErrorKind("Transition"), "Transition \"" + tname +
                                                                            "\" has inhibitor arcs, which WF-nets do "
                                                                            "not have.");
            if (!tp.firingdep.empty() && tp.firingdep[0])
                throw UnsupportedError(ErrorKind("Transition"),
                                       "Transition \"" + tname +
                                           "\" has a marking-dependent firing rate, which WF-net analysis does not "
                                           "support.");
            const lang::TimingStrategy timing = tp.timing[0];
            const Dist dist = timing == lang::TimingStrategy::IMMEDIATE ? Dist(Immediate()) : tp.firingproc[0];
            net.add_transition(tname, dist, timing, tp.nmodeservers[0], tp.firingprio[0], tp.fireweight[0]);
            for (std::size_t j = 0; j < Ein.rows() && j < nn; ++j)
                if (Ein(j, 0) > 0) net.add_arc(sn.nodes[j].name, tname, Ein(j, 0));
            for (std::size_t j = 0; j < Eout.rows() && j < nn; ++j) {
                if (!(Eout(j, 0) > 0)) continue;
                if (sn.nodes[j].nodetype != lang::NodeType::Place)
                    throw UnsupportedError(ErrorKind("Node"), "Transition \"" + tname + "\" fires into \"" +
                                                                  sn.nodes[j].name + "\", which is not a place.");
                net.add_arc(tname, sn.nodes[j].name, Eout(j, 0));
            }
        }
        net.validate();
        return net;
    }
 
private:
    typedef std::vector<long> Marking;
 
    struct Inc {
        std::size_t nP = 0, nT = 0;
        std::vector<std::vector<long>> pre, post;   ///< [place][transition]
    };
 
    struct ReachEdge {
        std::size_t from, to, t;
    };
 
    struct Reach {
        enum Status { COMPLETE, CAP, UNBOUNDED };
        std::vector<Marking> M;
        std::vector<ReachEdge> E;
        Status status = COMPLETE;
        Marking witness_cover, witness_ancestor;
    };
 
    /** The single-case chain of `singleCaseChain`, zero-time states folded out. */
    struct Chain {
        std::vector<double> alpha;
        Matrix<double> Q;
        Matrix<double> Rfire, Gfire;           ///< tangible x transitions
        std::vector<Marking> Mtang;            ///< marking of each tangible state
        std::vector<double> firings0;          ///< immediate firings before the first tangible state
    };
 
    struct WfState {
        workflow::Workflow<double> wf;
        std::size_t naux;
    };
 
    /** The element map of the last network built by `to_network`. */
    struct Compiled {
        std::shared_ptr<Network> network;
        std::string target, mode;
        long k = 0;
        double lambda = std::numeric_limits<double>::quiet_NaN();
        std::string restart;                    ///< closed fork-join: the restart Delay
        struct Task {
            std::string name, station, cls;
        };
        std::vector<Task> tasks;
        detail::TreeNode tree;
        std::vector<std::string> joins;         ///< Join names, in the preorder fjEmit created them
        std::vector<std::string> task_station;  ///< per transition, empty when it has none
        std::vector<std::string> task_class;
        struct Res {
            std::string name;
            double servers;
            std::string station;
            std::vector<std::size_t> busy;
        };
        std::vector<Res> resources;
        std::string counter;                    ///< SPN: the wfn_cases place
        std::vector<std::string> places;
        struct Busy {
            std::string name, place;
        };
        std::vector<Busy> busy;
    };
 
    std::string name_;
    std::vector<std::unique_ptr<WFNPlace>> places_;
    std::vector<std::unique_ptr<WFNTransition>> transitions_;
    std::vector<WFNArc> arcs_;
    std::vector<std::unique_ptr<WFNResource>> resources_;
    std::map<std::string, std::size_t> place_map_, transition_map_;
    long source_ = -1, sink_ = -1;
    Dist arrival_;
    bool has_arrival_ = false;
    long closed_cases_ = 0;
    double state_cap_ = 1e5;
    std::shared_ptr<Compiled> compiled_;
 
    // ------------------------------------------------------------------
    // Helpers of the public API
    // ------------------------------------------------------------------
 
    void check_new_name(const std::string& name) const {
        if (place_map_.count(name) || transition_map_.count(name))
            throw InputError(ErrorKind("Model"), "A place or transition named \"" + name +
                                                     "\" already exists in workflow net \"" + name_ + "\".");
        const std::string pre = ReservedPrefix;
        if (name.compare(0, pre.size(), pre) == 0)
            throw InputError(ErrorKind("Model"), "Names starting with \"" + pre +
                                                     "\" are reserved for compiled networks: \"" + name + "\".");
    }
 
    /** A name resolved to the place or transition it denotes. */
    ElementRef resolve(const ElementRef& ref) const {
        if (!ref.by_name_) return ref;
        const auto ip = place_map_.find(ref.name_);
        if (ip != place_map_.end()) return ElementRef(*places_[ip->second]);
        const auto it = transition_map_.find(ref.name_);
        if (it != transition_map_.end()) return ElementRef(*transitions_[it->second]);
        throw InputError(ErrorKind("Model"), "No place or transition named \"" + ref.name_ + "\" found.");
    }
 
    void check_owned(const ElementRef& ref) const {
        bool ok = false;
        if (ref.place_)
            ok = ref.place_->index() < places_.size() && places_[ref.place_->index()].get() == ref.place_;
        else if (ref.transition_)
            ok = ref.transition_->index() < transitions_.size() &&
                 transitions_[ref.transition_->index()].get() == ref.transition_;
        if (!ok)
            throw InputError(ErrorKind("Model"),
                             "Element is not a place or transition of workflow net \"" + name_ + "\".");
    }
 
    static std::string element_name(const ElementRef& ref) {
        return ref.place_ ? ref.place_->name() : ref.transition_->name();
    }
 
    long declared_place(const ElementRef& ref, const char* what) const {
        ElementRef r = ref;
        if (ref.by_name_) r = ElementRef(get_place(ref.name_));
        if (!r.place_) throw InputError(ErrorKind("Model"), std::string(what) + " must be a WFNPlace object.");
        check_owned(r);
        return static_cast<long>(r.place_->index());
    }
 
    WFNResource& set_resource_impl(const ElementRef& transition, const std::string& name, const double* servers,
                                   const lang::SchedStrategy* scheduling) {
        ElementRef r = transition;
        if (transition.by_name_) {
            const auto it = transition_map_.find(transition.name_);
            if (it == transition_map_.end())
                throw InputError(ErrorKind("Model"), "Transition \"" + transition.name_ + "\" not found.");
            r = ElementRef(*transitions_[it->second]);
        }
        if (!r.transition_) throw InputError(ErrorKind("Argument"), "First argument must be a transition.");
        check_owned(r);
        WFNTransition& tr = *transitions_[r.transition_->index()];
        if (tr.is_immediate())
            throw InputError(ErrorKind("Argument"), "Immediate transition \"" + tr.name() +
                                                        "\" takes no time and cannot hold a resource.");
        if (std::isfinite(tr.servers()))
            throw InputError(ErrorKind("Argument"),
                             "Transition \"" + tr.name() + "\" already limits its own concurrency to " +
                                 detail::num_d(tr.servers()) +
                                 " servers; a task is constrained either by its own servers or by a shared "
                                 "resource, not both.");
        long idx = -1;
        for (std::size_t i = 0; i < resources_.size(); ++i)
            if (resources_[i]->name() == name) {
                idx = static_cast<long>(i);
                break;
            }
        if (idx < 0) {
            resources_.emplace_back(new WFNResource(name, servers ? *servers : 1.0,
                                                    scheduling ? *scheduling : lang::SchedStrategy::FCFS));
            idx = static_cast<long>(resources_.size()) - 1;
        } else {
            if (servers) resources_[idx]->set_servers(*servers);
            if (scheduling) resources_[idx]->set_scheduling(*scheduling);
        }
        tr.resource_ = idx;
        return *resources_[idx];
    }
 
    void set_arrival_typed(const std::string& type, const double* rate, const Dist* map) {
        std::string tp = type;
        std::transform(tp.begin(), tp.end(), tp.begin(), [](unsigned char c) { return std::tolower(c); });
        if (tp == "none") {
            has_arrival_ = false;
            arrival_ = Dist();
        } else if (tp == "poisson") {
            if (!rate || !(*rate > 0))
                throw InputError(ErrorKind("Argument"), "A Poisson stream needs a positive rate.");
            arrival_ = Exp(*rate);
            has_arrival_ = true;
        } else if (tp == "map") {
            const bool modulated =
                map && (map->type == lang::ProcessType::MAP || map->type == lang::ProcessType::MMPP2 ||
                        map->type == lang::ProcessType::MMAP || map->type == lang::ProcessType::BMAP);
            if (!modulated) throw InputError(ErrorKind("Argument"), "A MAP stream needs a MAP object.");
            arrival_ = *map;
            has_arrival_ = true;
        } else {
            throw InputError(ErrorKind("Argument"), "Unknown arrival process type: " + type);
        }
    }
 
    std::string place_list(const std::vector<std::size_t>& idx) const {
        std::string s;
        for (std::size_t j = 0; j < idx.size(); ++j) s += (j ? ", \"" : "\"") + places_[idx[j]]->name() + "\"";
        return s.empty() ? "none" : s;
    }
 
    std::string transition_list(const std::vector<std::size_t>& idx) const {
        std::string s;
        for (std::size_t j = 0; j < idx.size(); ++j)
            s += (j ? ", \"" : "\"") + transitions_[idx[j]]->name() + "\"";
        return s.empty() ? "none" : s;
    }
 
    std::string marking_string(const Marking& m) const {
        std::string s = "[";
        bool first = true;
        for (std::size_t p = 0; p < m.size(); ++p) {
            if (m[p] == 0) continue;
            if (!first) s += ", ";
            s += places_[p]->name() + ":" + std::to_string(m[p]);
            first = false;
        }
        return s + "]";
    }
 
    Inc incidence() const {
        Inc inc;
        inc.nP = places_.size();
        inc.nT = transitions_.size();
        inc.pre.assign(inc.nP, std::vector<long>(inc.nT, 0));
        inc.post.assign(inc.nP, std::vector<long>(inc.nT, 0));
        for (const WFNArc& a : arcs_) {
            const long w = static_cast<long>(std::llround(a.weight()));
            if (a.from_place())
                inc.pre[a.place().index()][a.transition().index()] += w;
            else
                inc.post[a.place().index()][a.transition().index()] += w;
        }
        return inc;
    }
 
    static std::vector<bool> bfs(const std::vector<std::vector<std::size_t>>& adj,
                                 const std::vector<std::size_t>& starts) {
        std::vector<bool> visited(adj.size(), false);
        std::deque<std::size_t> queue;
        for (std::size_t s : starts) {
            visited[s] = true;
            queue.push_back(s);
        }
        while (!queue.empty()) {
            const std::size_t cur = queue.front();
            queue.pop_front();
            for (std::size_t nxt : adj[cur])
                if (!visited[nxt]) {
                    visited[nxt] = true;
                    queue.push_back(nxt);
                }
        }
        return visited;
    }
 
    static std::vector<double> enabling_degree(const Inc& inc, const Marking& m) {
        std::vector<double> deg(inc.nT, std::numeric_limits<double>::infinity());
        for (std::size_t t = 0; t < inc.nT; ++t) {
            for (std::size_t p = 0; p < inc.nP; ++p)
                if (inc.pre[p][t] > 0)
                    deg[t] = std::min(deg[t], std::floor(static_cast<double>(m[p]) / inc.pre[p][t]));
        }
        return deg;
    }
 
    /**
     * Breadth-first reachability graph of the P/T net. A new marking that
     * strictly covers a marking on its parent chain proves unboundedness; the
     * state cap stops the search with status CAP.
     */
    static Reach untimed_reachability(const Inc& inc, const Marking& m0, double cap) {
        Reach rg;
        std::vector<std::size_t> parent;
        std::map<Marking, std::size_t> keys;
        rg.M.push_back(m0);
        parent.push_back(npos_parent());
        keys[m0] = 0;
        for (std::size_t head = 0; head < rg.M.size(); ++head) {
            const Marking m = rg.M[head];
            for (std::size_t t = 0; t < inc.nT; ++t) {
                bool en = true;
                for (std::size_t p = 0; p < inc.nP && en; ++p) en = m[p] >= inc.pre[p][t];
                if (!en) continue;
                Marking m2(inc.nP);
                for (std::size_t p = 0; p < inc.nP; ++p) m2[p] = m[p] + inc.post[p][t] - inc.pre[p][t];
                std::size_t to;
                const auto it = keys.find(m2);
                if (it != keys.end()) {
                    to = it->second;
                } else {
                    for (std::size_t a = head; a != npos_parent(); a = parent[a]) {
                        bool ge = true, gt = false;
                        for (std::size_t p = 0; p < inc.nP; ++p) {
                            ge = ge && m2[p] >= rg.M[a][p];
                            gt = gt || m2[p] > rg.M[a][p];
                        }
                        if (ge && gt) {
                            rg.status = Reach::UNBOUNDED;
                            rg.witness_cover = m2;
                            rg.witness_ancestor = rg.M[a];
                            return rg;
                        }
                    }
                    if (static_cast<double>(rg.M.size()) >= cap) {
                        rg.status = Reach::CAP;
                        return rg;
                    }
                    rg.M.push_back(m2);
                    parent.push_back(head);
                    to = rg.M.size() - 1;
                    keys[m2] = to;
                }
                rg.E.push_back(ReachEdge{head, to, t});
            }
        }
        return rg;
    }
 
    static std::size_t npos_parent() { return static_cast<std::size_t>(-1); }
 
    void require_sound(long k) {
        const SoundnessReport report = check_soundness(static_cast<double>(k));
        if (report.verdict == "unknown")
            throw NumericError(ErrorKind("StateSpace"), "Soundness of \"" + name_ + "\" (k=" + std::to_string(k) +
                                                            ") is undecided: " + format_report(report) +
                                                            ". Raise stateCap, or skip the gate where the method "
                                                            "allows it.");
        if (!report.sound)
            throw UnsupportedError(ErrorKind("Structure"), "The WF-net \"" + name_ + "\" is not " +
                                                               std::to_string(k) +
                                                               "-sound, so the case cycle time is undefined: " +
                                                               format_report(report));
    }
 
    /** `phRepresentation`: the (alpha, S) pair of a phase-type firing law. */
    static void ph_representation(const Dist& d, const std::string& name, std::vector<double>& alpha,
                                  Matrix<double>& S) {
        const lang::ProcessType ty = d.type;
        const bool ph_family = ty == lang::ProcessType::EXP || ty == lang::ProcessType::ERLANG ||
                               ty == lang::ProcessType::HYPEREXP || ty == lang::ProcessType::PH ||
                               ty == lang::ProcessType::APH || ty == lang::ProcessType::COXIAN ||
                               ty == lang::ProcessType::COX2;
        if (!ph_family || d.disabled || !d.has_map())
            throw UnsupportedError(ErrorKind("Distribution"),
                                   "Transition \"" + name + "\" has a " + lang::process_to_text(ty) +
                                       " firing law; the exact single-case analysis needs a phase-type law (Exp, "
                                       "Erlang, HyperExp, Coxian, Cox2, APH, PH). Compile with toNetwork('spn') and "
                                       "use a simulator (SolverSSA, SolverLDES) for other laws.");
        const std::size_t n = d.D0.rows();
        S = d.D0;
        alpha.assign(n, 0.0);
        if ((ty == lang::ProcessType::PH || ty == lang::ProcessType::APH) && d.params.size() == n) {
            alpha.assign(d.params.begin(), d.params.end());
        } else {
            // D1 = (-D0 e) alpha for a renewal law, so any row with a positive exit recovers alpha
            for (std::size_t i = 0; i < n; ++i) {
                double out = 0.0;
                for (std::size_t j = 0; j < n; ++j) out += d.D1(i, j);
                if (!(out > 0)) continue;
                for (std::size_t j = 0; j < n; ++j) alpha[j] = d.D1(i, j) / out;
                break;
            }
        }
        double sum = 0.0;
        bool negative = false;
        for (double a : alpha) {
            sum += a;
            negative = negative || a < -lang::GlobalConstants::FineTol;
        }
        if (std::abs(sum - 1) > lang::GlobalConstants::CoarseTol || negative)
            throw UnsupportedError(ErrorKind("Distribution"), "The firing law of transition \"" + name +
                                                                  "\" has a defective or negative initial vector.");
    }
 
    /** `afterFiring`: marking and phases after T fires, under enabling memory. */
    static std::vector<long> after_firing(const Inc& inc, const std::vector<long>& x, std::size_t t,
                                          const std::vector<std::size_t>& ord) {
        const std::size_t nP = inc.nP, nT = inc.nT;
        Marking m2(nP);
        for (std::size_t p = 0; p < nP; ++p) m2[p] = x[p] + inc.post[p][t] - inc.pre[p][t];
        const std::vector<double> deg = enabling_degree(inc, m2);
        std::vector<long> y(nP + nT);
        for (std::size_t p = 0; p < nP; ++p) y[p] = m2[p];
        for (std::size_t u = 0; u < nT; ++u) {
            long ph = x[nP + u];
            if (deg[u] == 0 || u == t || ord[u] <= 1) ph = 0;
            y[nP + u] = ph;
        }
        return y;
    }
 
    /**
     * `singleCaseChain`: the absorbing CTMC of one case. A state is the marking
     * plus the phase of every enabled phase-type transition with more than one
     * phase. Zero-time states (an immediate transition is enabled, or a newly
     * enabled transition still has to draw its initial phase) are eliminated by
     * solving for their absorption probabilities into tangible states, and the
     * expected immediate firings are carried along.
     */
    Chain single_case_chain() {
        require_sound(1);
        const Inc inc = incidence();
        const std::size_t nP = inc.nP, nT = inc.nT;
        std::vector<bool> is_imm(nT, false);
        std::vector<double> prio(nT), wgt(nT), srv(nT);
        std::vector<std::vector<double>> pha(nT);
        std::vector<Matrix<double>> phS(nT);
        std::vector<std::size_t> ord(nT, 0);
        for (std::size_t t = 0; t < nT; ++t) {
            const WFNTransition& tr = *transitions_[t];
            is_imm[t] = tr.is_immediate();
            prio[t] = tr.priority();
            wgt[t] = tr.weight();
            srv[t] = tr.servers();
            if (!is_imm[t]) {
                ph_representation(tr.distribution(), tr.name(), pha[t], phS[t]);
                ord[t] = pha[t].size();
            }
        }
        Marking mf(nP, 0);
        mf[sink_] = 1;
        std::vector<long> x0(nP + nT, 0);
        x0[source_] = 1;
 
        struct Edge {
            std::size_t from, to;
            double val;
            long lab;   ///< firing transition, -1 for a phase move
        };
        std::map<std::vector<long>, std::size_t> keys;
        std::vector<std::vector<long>> S;
        std::vector<int> kind;   // 1 = zero-time, 2 = tangible, 3 = absorbing
        std::vector<Edge> edges;
        S.push_back(x0);
        keys[x0] = 0;
        for (std::size_t head = 0; head < S.size(); ++head) {
            const std::vector<long> x = S[head];
            const Marking m(x.begin(), x.begin() + static_cast<long>(nP));
            const std::vector<double> deg = enabling_degree(inc, m);
            std::vector<std::vector<long>> succ;
            std::vector<double> val;
            std::vector<long> lab;
            bool imm_enabled = false;
            for (std::size_t t = 0; t < nT; ++t) imm_enabled = imm_enabled || (is_imm[t] && deg[t] > 0);
            int k;
            if (m == mf) {
                k = 3;
            } else if (imm_enabled) {
                k = 1;
                double pmax = -std::numeric_limits<double>::infinity();
                for (std::size_t t = 0; t < nT; ++t)
                    if (is_imm[t] && deg[t] > 0) pmax = std::max(pmax, prio[t]);
                std::vector<std::size_t> en;
                double wsum = 0.0;
                for (std::size_t t = 0; t < nT; ++t)
                    if (is_imm[t] && deg[t] > 0 && prio[t] == pmax) {
                        en.push_back(t);
                        wsum += wgt[t];
                    }
                for (std::size_t t : en) {
                    succ.push_back(after_firing(inc, x, t, ord));
                    val.push_back(wgt[t] / wsum);
                    lab.push_back(static_cast<long>(t));
                }
            } else {
                std::vector<double> act(nT);
                for (std::size_t t = 0; t < nT; ++t) act[t] = std::min(deg[t], srv[t]);
                for (std::size_t t = 0; t < nT; ++t)
                    if (!is_imm[t] && ord[t] > 1 && act[t] > 1)
                        throw UnsupportedError(
                            ErrorKind("Distribution"),
                            "Transition \"" + transitions_[t]->name() + "\" is enabled " + detail::num_d(act[t]) +
                                " times at once with a " + std::to_string(ord[t]) +
                                "-phase firing law; the exact single-case analysis keeps one phase per transition, "
                                "so it supports concurrent firings only for exponential laws.");
                std::size_t need = nT;
                for (std::size_t t = 0; t < nT && need == nT; ++t)
                    if (!is_imm[t] && ord[t] > 1 && act[t] > 0 && x[nP + t] == 0) need = t;
                if (need < nT) {
                    k = 1;
                    for (std::size_t j = 0; j < pha[need].size(); ++j)
                        if (pha[need][j] > 0) {
                            std::vector<long> y = x;
                            y[nP + need] = static_cast<long>(j + 1);
                            succ.push_back(y);
                            val.push_back(pha[need][j]);
                            lab.push_back(-1);
                        }
                } else {
                    k = 2;
                    bool any = false;
                    for (std::size_t t = 0; t < nT; ++t) {
                        if (is_imm[t] || !(act[t] > 0)) continue;
                        any = true;
                        if (ord[t] == 1) {
                            succ.push_back(after_firing(inc, x, t, ord));
                            val.push_back(-phS[t](0, 0) * act[t]);
                            lab.push_back(static_cast<long>(t));
                        } else {
                            const std::size_t j = static_cast<std::size_t>(x[nP + t] - 1);
                            for (std::size_t j2 = 0; j2 < ord[t]; ++j2)
                                if (j2 != j && phS[t](j, j2) > 0) {
                                    std::vector<long> y = x;
                                    y[nP + t] = static_cast<long>(j2 + 1);
                                    succ.push_back(y);
                                    val.push_back(phS[t](j, j2));
                                    lab.push_back(-1);
                                }
                            double exit_rate = 0.0;
                            for (std::size_t j2 = 0; j2 < ord[t]; ++j2) exit_rate -= phS[t](j, j2);
                            if (exit_rate > lang::GlobalConstants::FineTol) {
                                succ.push_back(after_firing(inc, x, t, ord));
                                val.push_back(exit_rate);
                                lab.push_back(static_cast<long>(t));
                            }
                        }
                    }
                    if (!any)
                        throw InternalError(ErrorKind("Invariant"), "Sound net \"" + name_ +
                                                                        "\" reached a dead marking " +
                                                                        marking_string(m) + " in the timed chain.");
                }
            }
            kind.push_back(k);
            for (std::size_t e = 0; e < succ.size(); ++e) {
                std::size_t to;
                const auto it = keys.find(succ[e]);
                if (it != keys.end()) {
                    to = it->second;
                } else {
                    to = S.size();
                    if (static_cast<double>(to + 1) > state_cap_)
                        throw NumericError(ErrorKind("StateSpace"), "The single-case chain of \"" + name_ +
                                                                        "\" exceeds the state cap of " +
                                                                        detail::num_d(state_cap_) + " states.");
                    S.push_back(succ[e]);
                    keys[succ[e]] = to;
                }
                edges.push_back(Edge{head, to, val[e], lab[e]});
            }
        }
 
        const std::size_t n = S.size();
        std::vector<std::size_t> pos(n, 0);
        std::size_t nZ = 0, nTT = 0;
        for (std::size_t i = 0; i < n; ++i) {
            if (kind[i] == 1) pos[i] = nZ++;
            else if (kind[i] == 2) pos[i] = nTT++;
        }
        // zero-time part: absorption probabilities into tangible states, and the firings on the way
        Matrix<double> IZ(nZ, nZ, 0.0), PZT(nZ, nTT, 0.0), CZ(nZ, nT, 0.0);
        for (std::size_t z = 0; z < nZ; ++z) IZ(z, z) = 1.0;
        for (const Edge& e : edges) {
            if (kind[e.from] != 1) continue;
            const std::size_t z = pos[e.from];
            if (kind[e.to] == 1) IZ(z, pos[e.to]) -= e.val;
            else if (kind[e.to] == 2) PZT(z, pos[e.to]) += e.val;
            if (e.lab >= 0) CZ(z, static_cast<std::size_t>(e.lab)) += e.val;
        }
        Matrix<double> BT(nZ, nTT, 0.0), FZ(nZ, nT, 0.0);
        if (nZ > 0) {
            const detail::Factor f(IZ);
            for (std::size_t c = 0; c < nTT; ++c) {
                std::vector<double> b(nZ);
                for (std::size_t z = 0; z < nZ; ++z) b[z] = PZT(z, c);
                b = f.solve(b);
                for (std::size_t z = 0; z < nZ; ++z) BT(z, c) = b[z];
            }
            for (std::size_t c = 0; c < nT; ++c) {
                std::vector<double> b(nZ);
                for (std::size_t z = 0; z < nZ; ++z) b[z] = CZ(z, c);
                b = f.solve(b);
                for (std::size_t z = 0; z < nZ; ++z) FZ(z, c) = b[z];
            }
        }
        // tangible part: rates, with each zero-time successor replaced by where it lands
        Chain ch;
        ch.Q = Matrix<double>(nTT, nTT, 0.0);
        ch.Rfire = Matrix<double>(nTT, nT, 0.0);
        ch.Gfire = Matrix<double>(nTT, nT, 0.0);
        for (const Edge& e : edges) {
            if (kind[e.from] != 2) continue;
            const std::size_t s = pos[e.from];
            const double r = e.val;
            ch.Q(s, s) -= r;
            if (e.lab >= 0) ch.Rfire(s, static_cast<std::size_t>(e.lab)) += r;
            if (kind[e.to] == 2) {
                ch.Q(s, pos[e.to]) += r;
            } else if (kind[e.to] == 1) {
                const std::size_t z = pos[e.to];
                for (std::size_t j = 0; j < nTT; ++j) ch.Q(s, j) += r * BT(z, j);
                for (std::size_t j = 0; j < nT; ++j) ch.Gfire(s, j) += r * FZ(z, j);
            }
        }
        for (std::size_t i = 0; i < n; ++i)
            if (kind[i] == 2) ch.Mtang.push_back(Marking(S[i].begin(), S[i].begin() + static_cast<long>(nP)));
        ch.firings0.assign(nT, 0.0);
        if (kind[0] == 2) {
            ch.alpha.assign(nTT, 0.0);
            ch.alpha[pos[0]] = 1.0;
        } else if (kind[0] == 1) {
            ch.alpha.assign(nTT, 0.0);
            for (std::size_t j = 0; j < nTT; ++j) ch.alpha[j] = BT(pos[0], j);
            for (std::size_t j = 0; j < nT; ++j) ch.firings0[j] = FZ(pos[0], j);
        }
        return ch;
    }
 
    /** `choiceResolvable`: a choice the block structure can express. */
    bool choice_resolvable(const std::vector<std::size_t>& comp, const std::vector<bool>& head_imm,
                           const std::vector<double>& head_prio, const std::string& place_name,
                           std::string& why) const {
        bool ok = true;
        for (std::size_t c : comp) ok = ok && head_imm[c] && head_prio[c] == head_prio[comp[0]];
        why.clear();
        if (!ok)
            why = "the choice at place \"" + place_name +
                  "\" is not resolved by immediate transitions of equal priority, so it is a race or a priority "
                  "choice that the block structure cannot express";
        return ok;
    }
 
    /**
     * `reduceToTree`: reduce a block-structured WF-net to a process tree. A
     * LOOP node carries {do, redo} and the repeat probability q: the case runs
     * do, then with probability q runs redo and starts over. Immediate
     * transitions become EPS.
     */
    bool reduce_to_tree(detail::TreeNode& tree, std::string& reason) const {
        using detail::TreeNode;
        const Inc inc = incidence();
        reason.clear();
        const std::size_t nP = inc.nP, nT = inc.nT;
        for (std::size_t p = 0; p < nP; ++p)
            for (std::size_t t = 0; t < nT; ++t)
                if (inc.pre[p][t] > 1 || inc.post[p][t] > 1) {
                    reason = "an arc has weight greater than one";
                    return false;
                }
        std::vector<std::vector<bool>> pre(nP, std::vector<bool>(nT)), post(nP, std::vector<bool>(nT));
        for (std::size_t p = 0; p < nP; ++p)
            for (std::size_t t = 0; t < nT; ++t) {
                pre[p][t] = inc.pre[p][t] > 0;
                post[p][t] = inc.post[p][t] > 0;
            }
        std::vector<bool> alive_p(nP, true), alive_t(nT, true);
        const std::size_t src = static_cast<std::size_t>(source_);
        std::size_t snk = static_cast<std::size_t>(sink_);
        std::vector<TreeNode> ptree(nP, detail::tree_eps()), ttree(nT);
        std::vector<bool> head_imm(nT, false);
        std::vector<double> head_w(nT), head_prio(nT);
        for (std::size_t t = 0; t < nT; ++t) {
            const WFNTransition& tr = *transitions_[t];
            if (tr.is_immediate()) {
                ttree[t] = detail::tree_eps();
                head_imm[t] = true;
            } else {
                ttree[t].type = TreeNode::TASK;
                ttree[t].t = t;
            }
            head_w[t] = tr.weight();
            head_prio[t] = tr.priority();
        }
        auto col = [&](const std::vector<std::vector<bool>>& A, std::size_t t) {
            std::vector<std::size_t> out;
            for (std::size_t p = 0; p < nP; ++p)
                if (A[p][t]) out.push_back(p);
            return out;
        };
        auto row = [&](const std::vector<std::vector<bool>>& A, std::size_t p) {
            std::vector<std::size_t> out;
            for (std::size_t t = 0; t < nT; ++t)
                if (A[p][t]) out.push_back(t);
            return out;
        };
        auto alive = [](const std::vector<bool>& a) {
            std::vector<std::size_t> out;
            for (std::size_t i = 0; i < a.size(); ++i)
                if (a[i]) out.push_back(i);
            return out;
        };
        auto one = [](std::size_t x) { return std::vector<std::size_t>(1, x); };
        bool changed = true;
        while (changed) {
            changed = false;
            // R1: fold transition t between p1 and p2 into place p1
            for (std::size_t t : alive(alive_t)) {
                const std::vector<std::size_t> ins = col(pre, t), outs = col(post, t);
                if (ins.size() == 1 && outs.size() == 1 && ins[0] != outs[0]) {
                    const std::size_t p1 = ins[0], p2 = outs[0];
                    if (row(pre, p1) == one(t) && row(post, p2) == one(t)) {
                        ptree[p1] = detail::tree_seq({ptree[p1], ttree[t], ptree[p2]});
                        pre[p1] = pre[p2];
                        pre[p2].assign(nT, false);
                        post[p2].assign(nT, false);
                        for (std::size_t p = 0; p < nP; ++p) pre[p][t] = post[p][t] = false;
                        alive_t[t] = false;
                        alive_p[p2] = false;
                        if (p2 == snk) snk = p1;
                        changed = true;
                    }
                }
            }
            // R2: fold place p between t1 and t2 into transition t1
            for (std::size_t p : alive(alive_p)) {
                if (p == src || p == snk) continue;
                const std::vector<std::size_t> prod = row(post, p), cons = row(pre, p);
                if (prod.size() == 1 && cons.size() == 1 && prod[0] != cons[0]) {
                    const std::size_t t1 = prod[0], t2 = cons[0];
                    if (col(post, t1) == one(p) && col(pre, t2) == one(p)) {
                        ttree[t1] = detail::tree_seq({ttree[t1], ptree[p], ttree[t2]});
                        for (std::size_t q = 0; q < nP; ++q) {
                            post[q][t1] = post[q][t2];
                            pre[q][t2] = post[q][t2] = false;
                        }
                        pre[p].assign(nT, false);
                        post[p].assign(nT, false);
                        alive_t[t2] = false;
                        alive_p[p] = false;
                        changed = true;
                    }
                }
            }
            // R3: AND block, parallel places with the same single producer and consumer
            const std::vector<std::size_t> ap = alive(alive_p);
            for (std::size_t a : ap) {
                if (!alive_p[a] || a == src || a == snk) continue;
                const std::vector<std::size_t> pa = row(post, a), ca = row(pre, a);
                if (pa.size() != 1 || ca.size() != 1) continue;
                for (std::size_t b : ap) {
                    if (b <= a || !alive_p[b] || b == src || b == snk) continue;
                    if (row(post, b) == pa && row(pre, b) == ca) {
                        ptree[a] = detail::tree_and({ptree[a], ptree[b]});
                        pre[b].assign(nT, false);
                        post[b].assign(nT, false);
                        alive_p[b] = false;
                        changed = true;
                    }
                }
            }
            // R4: XOR block, parallel transitions with the same single input and output
            const std::vector<std::size_t> at = alive(alive_t);
            for (std::size_t a : at) {
                if (!alive_t[a]) continue;
                const std::vector<std::size_t> ia = col(pre, a), oa = col(post, a);
                if (ia.size() != 1 || oa.size() != 1) continue;
                for (std::size_t b : at) {
                    if (b <= a || !alive_t[b]) continue;
                    if (col(pre, b) == ia && col(post, b) == oa) {
                        const std::vector<std::size_t> comp = row(pre, ia[0]);
                        std::string why;
                        if (!choice_resolvable(comp, head_imm, head_prio, places_[ia[0]]->name(), why)) {
                            reason = why;
                            continue;
                        }
                        const double wa = head_w[a], wb = head_w[b];
                        ttree[a] = detail::tree_xor({ttree[a], ttree[b]}, {wa / (wa + wb), wb / (wa + wb)});
                        head_w[a] = wa + wb;
                        for (std::size_t p = 0; p < nP; ++p) pre[p][b] = post[p][b] = false;
                        alive_t[b] = false;
                        changed = true;
                    }
                }
            }
            // R5: loop, a transition whose only input and only output is p
            for (std::size_t t : alive(alive_t)) {
                const std::vector<std::size_t> ins = col(pre, t), outs = col(post, t);
                if (ins.size() == 1 && ins == outs) {
                    const std::size_t p = ins[0];
                    const std::vector<std::size_t> comp = row(pre, p);
                    if (comp.size() < 2) continue;
                    std::string why;
                    if (!choice_resolvable(comp, head_imm, head_prio, places_[p]->name(), why)) {
                        reason = why;
                        continue;
                    }
                    double wsum = 0.0;
                    for (std::size_t c : comp) wsum += head_w[c];
                    ptree[p] = detail::tree_loop(ptree[p], ttree[t], head_w[t] / wsum);
                    for (std::size_t q = 0; q < nP; ++q) pre[q][t] = post[q][t] = false;
                    alive_t[t] = false;
                    changed = true;
                }
            }
        }
        const std::vector<std::size_t> left_t = alive(alive_t), left_p = alive(alive_p);
        const bool ok = left_t.empty() && left_p.size() == 1;
        if (ok) {
            tree = ptree[left_p[0]];
        } else if (reason.empty()) {
            reason = "the reduction to series, AND, XOR and loop blocks stops with transitions " +
                     transition_list(left_t) + " left, so the net is not block-structured";
        }
        return ok;
    }
 
    // ------------------------------------------------------------------
    // to_workflow
    // ------------------------------------------------------------------
 
    std::string wf_task(std::size_t t, WfState& st) const {
        const WFNTransition& tr = *transitions_[t];
        if (st.wf.activity_index(tr.name()) == workflow::Workflow<double>::npos)
            st.wf.add_activity(tr.name(), tr.distribution());
        return tr.name();
    }
 
    static std::string wf_aux(WfState& st) {
        ++st.naux;
        const std::string name = std::string(ReservedPrefix) + "aux_" + std::to_string(st.naux);
        st.wf.add_activity(name, Immediate());
        return name;
    }
 
    /** `wfBranch`: emit NODE as a self-contained fragment, returning its head and tail activities. */
    void wf_branch(const detail::TreeNode& node, WfState& st, std::string& h, std::string& tl) const {
        using detail::TreeNode;
        switch (node.type) {
            case TreeNode::TASK:
                h = wf_task(node.t, st);
                tl = h;
                break;
            case TreeNode::SEQ: {
                std::vector<TreeNode> kids = node.kids;
                if (kids.front().type == TreeNode::TASK) {
                    h = wf_task(kids.front().t, st);
                    kids.erase(kids.begin());
                } else {
                    h = wf_aux(st);
                }
                if (!kids.empty() && kids.back().type == TreeNode::TASK) {
                    tl = wf_task(kids.back().t, st);
                    kids.pop_back();
                } else if (kids.empty()) {
                    tl = h;
                } else {
                    tl = wf_aux(st);
                }
                if (h != tl) wf_emit(detail::tree_seq(kids), h, tl, st);
                break;
            }
            case TreeNode::EPS:
                h = wf_aux(st);
                tl = h;
                break;
            default:
                h = wf_aux(st);
                tl = wf_aux(st);
                wf_emit(node, h, tl, st);
        }
    }
 
    /** `wfEmit`: emit the precedences leading from activity PRE through NODE to activity POST. */
    void wf_emit(const detail::TreeNode& node, const std::string& pre, const std::string& post,
                 WfState& st) const {
        using detail::TreeNode;
        typedef workflow::Workflow<double> Wf;
        switch (node.type) {
            case TreeNode::EPS:
                st.wf.add_precedence(Wf::Serial(pre, post));
                break;
            case TreeNode::TASK: {
                const std::string a = wf_task(node.t, st);
                st.wf.add_precedence(Wf::Serial(pre, a));
                st.wf.add_precedence(Wf::Serial(a, post));
                break;
            }
            case TreeNode::SEQ: {
                const std::vector<TreeNode>& kids = node.kids;
                std::string cur = pre;
                std::size_t i = 0;
                while (i < kids.size()) {
                    const TreeNode& k = kids[i];
                    if (k.type == TreeNode::TASK) {
                        const std::string a = wf_task(k.t, st);
                        st.wf.add_precedence(Wf::Serial(cur, a));
                        cur = a;
                        ++i;
                    } else {
                        std::string nxt;
                        std::size_t step;
                        if (i + 1 < kids.size() && kids[i + 1].type == TreeNode::TASK) {
                            nxt = wf_task(kids[i + 1].t, st);
                            step = 2;
                        } else if (i + 1 == kids.size()) {
                            nxt = post;
                            step = 1;
                        } else {
                            nxt = wf_aux(st);
                            step = 1;
                        }
                        wf_emit(k, cur, nxt, st);
                        cur = nxt;
                        i += step;
                    }
                }
                if (cur != post) st.wf.add_precedence(Wf::Serial(cur, post));
                break;
            }
            case TreeNode::AND:
            case TreeNode::XOR: {
                std::vector<std::string> heads(node.kids.size()), tails(node.kids.size());
                for (std::size_t j = 0; j < node.kids.size(); ++j) wf_branch(node.kids[j], st, heads[j], tails[j]);
                if (node.type == TreeNode::AND) {
                    st.wf.add_precedence(Wf::AndFork(pre, heads));
                    st.wf.add_precedence(Wf::AndJoin(tails, post));
                } else {
                    st.wf.add_precedence(Wf::OrFork(pre, heads, node.probs));
                    st.wf.add_precedence(Wf::OrJoin(tails, post));
                }
                break;
            }
            case TreeNode::LOOP: {
                if (node.kids[1].type != TreeNode::EPS)
                    throw UnsupportedError(ErrorKind("Structure"),
                                           "toWorkflow cannot express a loop with work on its return path: "
                                           "Workflow.Loop repeats a single body a geometric number of times. "
                                           "getCycleTime and toNetwork handle this net.");
                std::string h, tl;
                wf_branch(node.kids[0], st, h, tl);
                st.wf.add_precedence(Wf::Loop(pre, {h, post}, 1.0 / (1.0 - node.q)));
                break;
            }
        }
    }
 
    // ------------------------------------------------------------------
    // to_network("forkjoin")
    // ------------------------------------------------------------------
 
    /** A vertex of the flow graph fjEmit builds; VIRT vertices are eliminated by fj_resolve. */
    struct FjVertex {
        std::string kind;       ///< endpoint, sink, virt, station, fork, join, router
        std::size_t node = 0;   ///< 1-based network node, 0 for a virtual vertex
        std::size_t task = 0;   ///< 1-based transition whose class the vertex carries, 0 = case class
    };
    struct FjEdge {
        std::size_t u, v;
        double p;
    };
    struct FjGraph {
        std::vector<FjVertex> g;
        std::vector<FjEdge> E;
        std::size_t nfork = 0;
        std::vector<std::size_t> task_vertex;    ///< per transition, npos when it has none
        std::vector<std::size_t> res_node;       ///< per resource, 0 until created
        std::vector<std::size_t> task_node;      ///< per transition with its own Queue, else 0
        std::vector<std::string> joins;
        std::map<std::size_t, std::string> names;   ///< station node -> its name, for the AvgTable lookups
    };
 
    static std::size_t fj_vertex(FjGraph& fj, const std::string& kind, std::size_t node, std::size_t task) {
        FjVertex v;
        v.kind = kind;
        v.node = node;
        v.task = task;
        fj.g.push_back(v);
        return fj.g.size() - 1;
    }
 
    /** `fjEmit`: add the flow graph of NODE between the virtual vertices VIN and VOUT. */
    void fj_emit(const detail::TreeNode& node, std::size_t vin, std::size_t vout, FjGraph& fj, Network& net) const {
        using detail::TreeNode;
        switch (node.type) {
            case TreeNode::EPS:
                fj.E.push_back(FjEdge{vin, vout, 1.0});
                break;
            case TreeNode::TASK: {
                const std::size_t t = node.t;
                const WFNTransition& tr = *transitions_[t];
                std::size_t nd;
                if (tr.resource_ >= 0) {
                    const std::size_t r = static_cast<std::size_t>(tr.resource_);
                    if (fj.res_node[r] == 0) {
                        const WFNResource& res = *resources_[r];
                        fj.res_node[r] = net.add_queue(res.name(), res.scheduling());
                        net.set_number_of_servers(fj.res_node[r], res.servers());
                        fj.names[fj.res_node[r]] = res.name();
                    }
                    nd = fj.res_node[r];
                } else if (std::isfinite(tr.servers())) {
                    nd = net.add_queue(tr.name(), lang::SchedStrategy::FCFS);
                    net.set_number_of_servers(nd, tr.servers());
                    fj.task_node[t] = nd;
                    fj.names[nd] = tr.name();
                } else {
                    nd = net.add_delay(tr.name());
                    fj.names[nd] = tr.name();
                }
                const std::size_t v = fj_vertex(fj, "station", nd, t + 1);
                fj.task_vertex[t] = v;
                fj.E.push_back(FjEdge{vin, v, 1.0});
                fj.E.push_back(FjEdge{v, vout, 1.0});
                break;
            }
            case TreeNode::SEQ: {
                std::size_t cur = vin;
                for (std::size_t j = 0; j < node.kids.size(); ++j) {
                    const std::size_t nxt = j + 1 == node.kids.size() ? vout : fj_vertex(fj, "virt", 0, 0);
                    fj_emit(node.kids[j], cur, nxt, fj, net);
                    cur = nxt;
                }
                break;
            }
            case TreeNode::XOR:
                for (std::size_t j = 0; j < node.kids.size(); ++j) {
                    const std::size_t b = fj_vertex(fj, "virt", 0, 0);
                    fj.E.push_back(FjEdge{vin, b, node.probs[j]});
                    fj_emit(node.kids[j], b, vout, fj, net);
                }
                break;
            case TreeNode::LOOP: {
                const std::size_t v1 = fj_vertex(fj, "virt", 0, 0);
                const std::size_t v2 = fj_vertex(fj, "virt", 0, 0);
                const std::size_t v3 = fj_vertex(fj, "virt", 0, 0);
                fj.E.push_back(FjEdge{vin, v1, 1.0});
                fj_emit(node.kids[0], v1, v2, fj, net);
                fj.E.push_back(FjEdge{v2, v3, node.q});
                fj_emit(node.kids[1], v3, v1, fj, net);
                fj.E.push_back(FjEdge{v2, vout, 1.0 - node.q});
                break;
            }
            case TreeNode::AND: {
                ++fj.nfork;
                const std::string pre = ReservedPrefix;
                const std::size_t fk = net.add_fork(pre + "fork" + std::to_string(fj.nfork));
                const std::string jname = pre + "join" + std::to_string(fj.nfork);
                const std::size_t jn = net.add_join(jname, fk);
                fj.joins.push_back(jname);
                const std::size_t vf = fj_vertex(fj, "fork", fk, 0);
                const std::size_t vj = fj_vertex(fj, "join", jn, 0);
                fj.E.push_back(FjEdge{vin, vf, 1.0});
                for (std::size_t j = 0; j < node.kids.size(); ++j) {
                    const std::size_t b = fj_vertex(fj, "virt", 0, 0);
                    const std::size_t bo = fj_vertex(fj, "virt", 0, 0);
                    fj.E.push_back(FjEdge{vf, b, 1.0});
                    fj_emit(node.kids[j], b, bo, fj, net);
                    fj.E.push_back(FjEdge{bo, vj, 1.0});
                }
                fj.E.push_back(FjEdge{vj, vout, 1.0});
                break;
            }
        }
    }
 
    /**
     * `fjResolve`: eliminate the virtual vertices. A fork branch that is not
     * deterministic is routed through a Router, since a Fork sends one job down
     * every outgoing link.
     */
    static std::vector<FjEdge> fj_resolve(FjGraph& fj, Network& net) {
        const std::size_t nv = fj.g.size();
        std::vector<bool> is_virt(nv);
        std::vector<std::size_t> iv, ip, vpos(nv, 0), ppos(nv, 0);
        for (std::size_t v = 0; v < nv; ++v) {
            is_virt[v] = fj.g[v].kind == "virt";
            if (is_virt[v]) {
                vpos[v] = iv.size();
                iv.push_back(v);
            } else {
                ppos[v] = ip.size();
                ip.push_back(v);
            }
        }
        const std::size_t ni = iv.size(), np = ip.size();
        Matrix<double> Iv(ni, ni, 0.0), Pvp(ni, np, 0.0);
        for (std::size_t i = 0; i < ni; ++i) Iv(i, i) = 1.0;
        for (const FjEdge& e : fj.E) {
            if (!is_virt[e.u]) continue;
            if (is_virt[e.v])
                Iv(vpos[e.u], vpos[e.v]) -= e.p;
            else
                Pvp(vpos[e.u], ppos[e.v]) += e.p;
        }
        const std::string loop_msg =
            "The net has a loop of zero-time choices that repeats with probability one, so a case never leaves it.";
        Matrix<double> Res(ni, np, 0.0);
        if (ni > 0) {
            Matrix<double> inv(ni, ni, 0.0);
            try {
                const detail::Factor f(Iv);
                for (std::size_t c = 0; c < ni; ++c) {
                    std::vector<double> b(ni, 0.0);
                    b[c] = 1.0;
                    b = f.solve(b);
                    for (std::size_t r = 0; r < ni; ++r) inv(r, c) = b[r];
                }
            } catch (const NumericError&) {
                throw UnsupportedError(ErrorKind("Structure"), loop_msg);
            }
            // condest(Iv) as the exact 1-norm condition number
            double n1 = 0.0, n1i = 0.0;
            for (std::size_t c = 0; c < ni; ++c) {
                double s = 0.0, si = 0.0;
                for (std::size_t r = 0; r < ni; ++r) {
                    s += std::abs(Iv(r, c));
                    si += std::abs(inv(r, c));
                }
                n1 = std::max(n1, s);
                n1i = std::max(n1i, si);
            }
            if (!(n1 * n1i <= 1.0 / lang::GlobalConstants::FineTol))
                throw UnsupportedError(ErrorKind("Structure"), loop_msg);
            for (std::size_t r = 0; r < ni; ++r)
                for (std::size_t c = 0; c < np; ++c) {
                    double s = 0.0;
                    for (std::size_t k = 0; k < ni; ++k) s += inv(r, k) * Pvp(k, c);
                    Res(r, c) = s;
                }
        }
        const double tol = lang::GlobalConstants::FineTol;
        std::vector<FjEdge> dist;
        std::size_t nrouter = 0;
        for (std::size_t u : ip) {
            if (fj.g[u].kind == "sink") continue;
            std::vector<double> d;
            if (fj.g[u].kind == "fork") {
                for (std::size_t ei = 0; ei < fj.E.size(); ++ei) {
                    const FjEdge e = fj.E[ei];
                    if (e.u != u) continue;
                    d.assign(np, 0.0);
                    if (is_virt[e.v])
                        for (std::size_t j = 0; j < np; ++j) d[j] = Res(vpos[e.v], j);
                    else
                        d[ppos[e.v]] = 1.0;
                    std::vector<std::size_t> tgt;
                    for (std::size_t j = 0; j < np; ++j)
                        if (d[j] > tol) tgt.push_back(j);
                    if (tgt.size() == 1 && std::abs(d[tgt[0]] - 1) < tol) {
                        dist.push_back(FjEdge{u, ip[tgt[0]], 1.0});
                    } else {
                        ++nrouter;
                        const std::size_t rt =
                            net.add_router(std::string(ReservedPrefix) + "router" + std::to_string(nrouter));
                        const std::size_t r = fj_vertex(fj, "router", rt, 0);
                        dist.push_back(FjEdge{u, r, 1.0});
                        for (std::size_t j : tgt) dist.push_back(FjEdge{r, ip[j], d[j]});
                    }
                }
            } else {
                d.assign(np, 0.0);
                for (const FjEdge& e : fj.E) {
                    if (e.u != u) continue;
                    if (is_virt[e.v])
                        for (std::size_t j = 0; j < np; ++j) d[j] += e.p * Res(vpos[e.v], j);
                    else
                        d[ppos[e.v]] += e.p;
                }
                for (std::size_t j = 0; j < np; ++j)
                    if (d[j] > tol) dist.push_back(FjEdge{u, ip[j], d[j]});
            }
        }
        return dist;
    }
 
    /** `compileToForkJoin`: the fork-join queueing network of a block-structured net. */
    std::shared_ptr<Network> compile_to_fork_join() {
        detail::TreeNode tree;
        std::string reason;
        if (!reduce_to_tree(tree, reason))
            throw UnsupportedError(ErrorKind("Structure"), "toNetwork('forkjoin') requires a block-structured "
                                                           "WF-net: " + reason + ". Use toNetwork('spn').");
        const bool closed = closed_cases_ > 0;
        const std::string pre = ReservedPrefix;
        std::shared_ptr<Network> network = std::make_shared<Network>(name_ + "_forkjoin");
        Network& net = *network;
        const std::size_t nT = transitions_.size(), nR = resources_.size();
        // tasks on a resource shared with another task need their own class
        std::vector<std::size_t> use_count(nR, 0);
        for (std::size_t t = 0; t < nT; ++t)
            if (transitions_[t]->resource_ >= 0) ++use_count[static_cast<std::size_t>(transitions_[t]->resource_)];
        FjGraph fj;
        fj.task_vertex.assign(nT, static_cast<std::size_t>(-1));
        fj.res_node.assign(nR, 0);
        fj.task_node.assign(nT, 0);
        if (closed) {
            fj_vertex(fj, "endpoint", net.add_delay(pre + "restart"), 0);
        } else {
            fj_vertex(fj, "endpoint", net.add_source(pre + "source"), 0);
            fj_vertex(fj, "sink", net.add_sink(pre + "sink"), 0);
        }
        const std::size_t vin = fj_vertex(fj, "virt", 0, 0);
        const std::size_t vout = fj_vertex(fj, "virt", 0, 0);
        fj.E.push_back(FjEdge{0, vin, 1.0});
        fj.E.push_back(FjEdge{vout, closed ? 0u : 1u, 1.0});
        fj_emit(tree, vin, vout, fj, net);
        const std::vector<FjEdge> dist = fj_resolve(fj, net);
        // classes
        const std::size_t endpoint = fj.g[0].node;
        const std::size_t case_class = closed
                                           ? net.add_closed_class("Case", static_cast<double>(closed_cases_), endpoint, 0)
                                           : net.add_open_class("Case", 0);
        std::vector<std::size_t> cls_of(nT, case_class);
        std::vector<std::string> cls_name(nT, "Case");
        for (std::size_t t = 0; t < nT; ++t) {
            const long r = transitions_[t]->resource_;
            if (r >= 0 && use_count[static_cast<std::size_t>(r)] > 1) {
                cls_name[t] = "Case_" + transitions_[t]->name();
                cls_of[t] = closed ? net.add_closed_class(cls_name[t], 0.0, endpoint, 0)
                                   : net.add_open_class(cls_name[t], 0);
            }
        }
        std::vector<std::size_t> vcls(fj.g.size(), case_class);
        for (std::size_t v = 0; v < fj.g.size(); ++v)
            if (fj.g[v].task > 0) vcls[v] = cls_of[fj.g[v].task - 1];
        // service
        if (closed)
            net.set_service(endpoint, case_class, Immediate());
        else
            net.set_arrival(endpoint, case_class, arrival_);
        std::shared_ptr<Compiled> c = std::make_shared<Compiled>();
        c->task_station.assign(nT, std::string());
        for (std::size_t t = 0; t < nT; ++t) {
            if (fj.task_vertex[t] == static_cast<std::size_t>(-1)) continue;
            const WFNTransition& tr = *transitions_[t];
            const std::size_t nd = fj.g[fj.task_vertex[t]].node;
            net.set_service(nd, cls_of[t], tr.distribution());
            const std::string station = fj.names.at(nd);
            c->tasks.push_back(Compiled::Task{tr.name(), station, cls_name[t]});
            c->task_station[t] = station;
        }
        // routing
        RoutingMatrix P;
        for (const FjEdge& e : dist) {
            const std::size_t r = vcls[e.u];
            const std::size_t s = fj.g[e.v].kind == "sink" ? r : vcls[e.v];
            const std::size_t nu = fj.g[e.u].node, nv = fj.g[e.v].node;
            P.set(r, s, nu, nv, P.get(r, s, nu, nv) + e.p);
        }
        net.link(P);
        for (std::size_t r = 0; r < nR; ++r)
            if (fj.res_node[r] != 0)
                c->resources.push_back(Compiled::Res{resources_[r]->name(), resources_[r]->servers(),
                                                     fj.names.at(fj.res_node[r]), {}});
        for (std::size_t t = 0; t < nT; ++t)
            if (fj.task_node[t] != 0)
                c->resources.push_back(Compiled::Res{transitions_[t]->name(), transitions_[t]->servers(),
                                                     fj.names.at(fj.task_node[t]), {}});
        c->network = network;
        c->target = "forkjoin";
        c->k = closed_cases_;
        if (closed) {
            c->mode = "closed";
            c->restart = pre + "restart";
        } else {
            c->mode = "open";
            c->lambda = 1.0 / arrival_.mean;
        }
        c->tree = tree;
        c->joins = fj.joins;
        c->task_class = cls_name;
        compiled_ = c;
        return network;
    }
 
    /**
     * `fjTreeTime`: mean case time through NODE from the per-visit response
     * times of a fork-join network. An AND block takes the mean of its branch
     * times plus the response time of its join: a join reports the
     * synchronisation wait averaged over the sibling tasks, one per branch, and
     * every sibling leaves the block when the last one arrives. Joins are
     * visited in the preorder fjEmit created them in.
     */
    static double fj_tree_time(const detail::TreeNode& node, const Compiled& c, const AvgTable& tab,
                               std::size_t& jidx) {
        using detail::TreeNode;
        switch (node.type) {
            case TreeNode::EPS:
                return 0.0;
            case TreeNode::TASK:
                return detail::table_cell(tab, tab.RespT, c.task_station[node.t], c.task_class[node.t]);
            case TreeNode::SEQ: {
                double x = 0.0;
                for (const TreeNode& k : node.kids) x += fj_tree_time(k, c, tab, jidx);
                return x;
            }
            case TreeNode::XOR: {
                double x = 0.0;
                for (std::size_t j = 0; j < node.kids.size(); ++j)
                    x += node.probs[j] * fj_tree_time(node.kids[j], c, tab, jidx);
                return x;
            }
            case TreeNode::LOOP: {
                const double d = fj_tree_time(node.kids[0], c, tab, jidx);
                const double r = fj_tree_time(node.kids[1], c, tab, jidx);
                return (d + node.q * r) / (1.0 - node.q);
            }
            case TreeNode::AND: {
                const std::string& jn = c.joins.at(jidx++);
                double y = 0.0;
                for (const TreeNode& k : node.kids) y += fj_tree_time(k, c, tab, jidx);
                // fjEmit routes every task into a join in the case class
                return y / static_cast<double>(node.kids.size()) + detail::table_cell(tab, tab.RespT, jn, "Case");
            }
        }
        return 0.0;
    }
 
    // ------------------------------------------------------------------
    // to_network("spn") and the plain net
    // ------------------------------------------------------------------
 
    /** `setTransitionTiming`: one mode of a LINE Transition from a WF-net transition. */
    static void set_transition_timing(Network& net, std::size_t node, std::size_t mode, const WFNTransition& w,
                                      double servers) {
        if (w.is_immediate()) {
            net.set_mode_timing(node, mode, lang::TimingStrategy::IMMEDIATE);
            net.set_firing_weight(node, mode, w.weight());
            net.set_firing_priority(node, mode, w.priority());
        } else {
            net.set_mode_timing(node, mode, lang::TimingStrategy::TIMED);
            net.set_mode_distribution(node, mode, w.distribution());
            net.set_mode_servers(node, mode, servers);
            net.set_firing_priority(node, mode, w.priority());
        }
    }
 
    /**
     * `compileToSPN`: a stochastic Petri net with a case-counter place. Arrivals
     * (open) or the k recirculating cases (closed) wait in wfn_arrivals; the
     * immediate wfn_start moves one into the source and adds a token to
     * wfn_cases, and the immediate wfn_end takes the token in the sink together
     * with one of wfn_cases. The mean of wfn_cases is therefore the mean number
     * of cases in the system.
     */
    std::shared_ptr<Network> compile_to_spn() {
        const bool closed = closed_cases_ > 0;
        const Inc inc = incidence();
        const std::size_t nP = inc.nP, nT = inc.nT, nR = resources_.size();
        for (std::size_t t = 0; t < nT; ++t) {
            const WFNTransition& tr = *transitions_[t];
            if (tr.resource_ < 0) continue;
            for (std::size_t p = 0; p < nP; ++p) {
                if (inc.pre[p][t] == 0) continue;
                std::vector<std::size_t> rivals;
                for (std::size_t u = 0; u < nT; ++u)
                    if (u != t && inc.pre[p][u] > 0) rivals.push_back(u);
                if (!rivals.empty())
                    throw UnsupportedError(ErrorKind("Structure"),
                                           "Task \"" + tr.name() + "\" holds resource \"" +
                                               resources_[static_cast<std::size_t>(tr.resource_)]->name() +
                                               "\" and competes for place \"" + places_[p]->name() + "\" with " +
                                               transition_list(rivals) +
                                               "; in the SPN the resource is acquired by an immediate transition, "
                                               "which would always win that race.");
            }
        }
        const std::string pre = ReservedPrefix;
        std::shared_ptr<Network> network = std::make_shared<Network>(name_ + "_spn");
        Network& net = *network;
        std::size_t source = 0, sink = 0;
        if (!closed) {
            source = net.add_source(pre + "source");
            sink = net.add_sink(pre + "sink");
        }
        std::vector<std::size_t> pl(nP);
        for (std::size_t p = 0; p < nP; ++p) pl[p] = net.add_place(places_[p]->name());
        const std::size_t p_arr = net.add_place(pre + "arrivals");
        const std::size_t p_cnt = net.add_place(pre + "cases");
        std::vector<std::size_t> p_res(nR);
        for (std::size_t r = 0; r < nR; ++r) p_res[r] = net.add_place(pre + "res_" + resources_[r]->name());
        std::vector<std::size_t> tr_(nT), acq(nT, 0), p_busy(nT, 0);
        for (std::size_t t = 0; t < nT; ++t) {
            tr_[t] = net.add_transition(transitions_[t]->name());
            if (transitions_[t]->resource_ >= 0) {
                acq[t] = net.add_transition(pre + "acq_" + transitions_[t]->name());
                p_busy[t] = net.add_place(pre + "busy_" + transitions_[t]->name());
            }
        }
        const std::size_t t_start = net.add_transition(pre + "start");
        const std::size_t t_end = net.add_transition(pre + "end");
        std::size_t jc;
        if (closed) {
            // resource tokens are tokens of the same uncoloured class, so they count in the population
            double n_res = 0.0;
            for (std::size_t r = 0; r < nR; ++r) n_res += resources_[r]->servers();
            jc = net.add_closed_class("Case", static_cast<double>(closed_cases_) + n_res, p_arr, 0);
        } else {
            jc = net.add_open_class("Case", 0);
            net.set_arrival(source, jc, arrival_);
        }
        RoutingMatrix R;
        auto arc = [&R](std::size_t a, std::size_t b) { R.set(1, 1, a, b, 1.0); };
        // start and end
        std::size_t m = net.add_mode(t_start, "start");
        net.set_mode_timing(t_start, m, lang::TimingStrategy::IMMEDIATE);
        net.set_enabling_conditions(t_start, m, jc, p_arr, 1);
        net.set_firing_outcome(t_start, m, jc, pl[source_], 1);
        net.set_firing_outcome(t_start, m, jc, p_cnt, 1);
        arc(p_arr, t_start);
        arc(t_start, pl[source_]);
        arc(t_start, p_cnt);
        m = net.add_mode(t_end, "end");
        net.set_mode_timing(t_end, m, lang::TimingStrategy::IMMEDIATE);
        net.set_enabling_conditions(t_end, m, jc, pl[sink_], 1);
        net.set_enabling_conditions(t_end, m, jc, p_cnt, 1);
        arc(pl[sink_], t_end);
        arc(p_cnt, t_end);
        if (closed) {
            net.set_firing_outcome(t_end, m, jc, p_arr, 1);
            arc(t_end, p_arr);
        } else {
            net.set_firing_outcome(t_end, m, jc, sink, 1);
            arc(t_end, sink);
            arc(source, p_arr);
        }
        std::shared_ptr<Compiled> c = std::make_shared<Compiled>();
        for (std::size_t r = 0; r < nR; ++r)
            c->resources.push_back(Compiled::Res{resources_[r]->name(), resources_[r]->servers(), std::string(), {}});
        for (std::size_t t = 0; t < nT; ++t) {
            const WFNTransition& w = *transitions_[t];
            if (w.resource_ < 0) {
                m = net.add_mode(tr_[t], w.name());
                set_transition_timing(net, tr_[t], m, w, w.servers());
                for (std::size_t p = 0; p < nP; ++p)
                    if (inc.pre[p][t] > 0) {
                        net.set_enabling_conditions(tr_[t], m, jc, pl[p], static_cast<double>(inc.pre[p][t]));
                        arc(pl[p], tr_[t]);
                    }
                for (std::size_t p = 0; p < nP; ++p)
                    if (inc.post[p][t] > 0) {
                        net.set_firing_outcome(tr_[t], m, jc, pl[p], static_cast<double>(inc.post[p][t]));
                        arc(tr_[t], pl[p]);
                    }
            } else {
                const std::size_t r = static_cast<std::size_t>(w.resource_);
                const std::size_t ma = net.add_mode(acq[t], "acq_" + w.name());
                net.set_mode_timing(acq[t], ma, lang::TimingStrategy::IMMEDIATE);
                net.set_firing_weight(acq[t], ma, w.weight());
                net.set_firing_priority(acq[t], ma, w.priority());
                for (std::size_t p = 0; p < nP; ++p)
                    if (inc.pre[p][t] > 0) {
                        net.set_enabling_conditions(acq[t], ma, jc, pl[p], static_cast<double>(inc.pre[p][t]));
                        arc(pl[p], acq[t]);
                    }
                net.set_enabling_conditions(acq[t], ma, jc, p_res[r], 1);
                arc(p_res[r], acq[t]);
                net.set_firing_outcome(acq[t], ma, jc, p_busy[t], 1);
                arc(acq[t], p_busy[t]);
                m = net.add_mode(tr_[t], w.name());
                set_transition_timing(net, tr_[t], m, w, std::numeric_limits<double>::infinity());
                net.set_enabling_conditions(tr_[t], m, jc, p_busy[t], 1);
                arc(p_busy[t], tr_[t]);
                for (std::size_t p = 0; p < nP; ++p)
                    if (inc.post[p][t] > 0) {
                        net.set_firing_outcome(tr_[t], m, jc, pl[p], static_cast<double>(inc.post[p][t]));
                        arc(tr_[t], pl[p]);
                    }
                net.set_firing_outcome(tr_[t], m, jc, p_res[r], 1);
                arc(tr_[t], p_res[r]);
                c->busy.push_back(Compiled::Busy{w.name(), pre + "busy_" + w.name()});
                c->resources[r].busy.push_back(c->busy.size() - 1);
            }
        }
        net.link(R);
        for (std::size_t p = 0; p < nP; ++p) net.set_initial_marking(pl[p], {0.0});
        net.set_initial_marking(p_cnt, {0.0});
        for (std::size_t t = 0; t < nT; ++t)
            if (p_busy[t] != 0) net.set_initial_marking(p_busy[t], {0.0});
        for (std::size_t r = 0; r < nR; ++r) net.set_initial_marking(p_res[r], {resources_[r]->servers()});
        net.set_initial_marking(p_arr, {closed ? static_cast<double>(closed_cases_) : 0.0});
        c->network = network;
        c->target = "spn";
        c->k = closed_cases_;
        if (closed) {
            c->mode = "closed";
        } else {
            c->mode = "open";
            c->lambda = 1.0 / arrival_.mean;
        }
        c->counter = pre + "cases";
        for (std::size_t p = 0; p < nP; ++p) c->places.push_back(places_[p]->name());
        std::vector<Compiled::Res> kept;
        for (const Compiled::Res& res : c->resources)
            if (res.servers > 0 && !res.busy.empty()) kept.push_back(res);
        c->resources = kept;
        compiled_ = c;
        return network;
    }
 
    /** `buildPlainNetwork(k)`: the net itself as a closed LINE Petri net, k tokens in the source. */
    std::shared_ptr<Network> build_plain_network(long k) const {
        const Inc inc = incidence();
        const std::size_t nP = inc.nP, nT = inc.nT;
        std::shared_ptr<Network> model = std::make_shared<Network>(name_);
        Network& net = *model;
        std::vector<std::size_t> pl(nP), tr(nT);
        for (std::size_t p = 0; p < nP; ++p) pl[p] = net.add_place(places_[p]->name());
        for (std::size_t t = 0; t < nT; ++t) tr[t] = net.add_transition(transitions_[t]->name());
        const std::size_t jc = net.add_closed_class("Case", static_cast<double>(k), pl[source_], 0);
        RoutingMatrix R;
        for (std::size_t t = 0; t < nT; ++t) {
            const WFNTransition& w = *transitions_[t];
            const std::size_t m = net.add_mode(tr[t], w.name());
            set_transition_timing(net, tr[t], m, w, w.servers());
            for (std::size_t p = 0; p < nP; ++p)
                if (inc.pre[p][t] > 0) {
                    net.set_enabling_conditions(tr[t], m, jc, pl[p], static_cast<double>(inc.pre[p][t]));
                    R.set(1, 1, pl[p], tr[t], 1.0);
                }
            for (std::size_t p = 0; p < nP; ++p)
                if (inc.post[p][t] > 0) {
                    net.set_firing_outcome(tr[t], m, jc, pl[p], static_cast<double>(inc.post[p][t]));
                    R.set(1, 1, tr[t], pl[p], 1.0);
                }
        }
        net.link(R);
        for (std::size_t p = 0; p < nP; ++p)
            net.set_initial_marking(pl[p], {static_cast<double>(static_cast<long>(p) == source_ ? k : 0)});
        return model;
    }
};
 
}  // namespace wfn
}  // namespace line
 
#endif  // LINE_API_WFN_WORKFLOW_NET_H