daq_graph_model.py

"""
Generic DAQ Data-Flow Graph Model
==================================
Models data flow, utilization, and bottlenecks as a directed graph of
Node and Edge objects.  The DAQ topology is expressed purely as a list
of node/edge definitions at the bottom of this file — no topology logic
is buried in the analysis code.

Concepts
--------
Node   — a processing stage (detector, buffer, event builder, …).
         Has a multiplicity (how many parallel instances exist) and
         zero or more resource constraints (e.g. disk write speed).

Edge   — a directed data-flow link between two nodes.
         Has a raw bandwidth capacity and an optional safety margin.
         Carries a 'flow_Bps' that is computed by the solver.

Flow solver
-----------
Each edge is labelled with a symbolic flow expression — a lambda that
takes the global parameter dict and returns bytes/s.  The solver
evaluates all flows, then computes utilization = flow / effective_capacity
for every edge and node-resource.  The bottleneck is the
edge/resource with the highest utilization.
"""

from __future__ import annotations
import math
from dataclasses import dataclass, field
from typing import Callable, Dict, List, Optional, Tuple
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.gridspec import GridSpec
import warnings
warnings.filterwarnings("ignore")


# ══════════════════════════════════════════════════════════════
# CORE DATA STRUCTURES
# ══════════════════════════════════════════════════════════════

@dataclass
class Resource:
    """A throughput constraint that lives on a Node (e.g. disk write speed)."""
    name: str                              # human label
    capacity_fn: Callable                  # fn(params) → raw capacity in Bps
    flow_fn:     Callable                  # fn(params) → actual flow in Bps
    safety_margin: Optional[float] = None  # overrides graph-level default if set
    # computed
    capacity_Bps: float  = field(default=0.0, init=False)
    flow_Bps:     float  = field(default=0.0, init=False)
    utilization:  float  = field(default=0.0, init=False)
    max_trigger_rate_hz: Optional[float] = field(default=None, init=False)
    # optional: fn(params) → max sustainable trigger rate (Hz)
    max_rate_fn:  Optional[Callable] = None

    def evaluate(self, params: dict, default_margin: float):
        sm = self.safety_margin if self.safety_margin is not None else default_margin
        self.capacity_Bps = self.capacity_fn(params) * sm
        self.flow_Bps     = self.flow_fn(params)
        self.utilization  = self.flow_Bps / self.capacity_Bps if self.capacity_Bps else float("inf")
        if self.max_rate_fn:
            self.max_trigger_rate_hz = self.max_rate_fn(params, sm)


@dataclass
class Node:
    """A processing stage in the DAQ pipeline."""
    id:           str
    label:        str                      # display name (may contain \n)
    sublabel_fn:  Callable                 # fn(params) → str, shown under label
    color:        str = "#1f6feb"
    multiplicity_fn: Callable = lambda p: 1   # fn(params) → int
    resources:    List[Resource] = field(default_factory=list)
    # layout hint: (x, y) in figure-space units
    pos:          Tuple[float, float] = (0.0, 0.0)
    # computed
    multiplicity: int = field(default=1, init=False)

    def evaluate(self, params: dict, default_margin: float):
        self.multiplicity = int(self.multiplicity_fn(params))
        for r in self.resources:
            r.evaluate(params, default_margin)

    @property
    def worst_resource(self) -> Optional[Resource]:
        if not self.resources:
            return None
        return max(self.resources, key=lambda r: r.utilization)


@dataclass
class Edge:
    """A directed data-flow link between two nodes (pure flow-rate arrow)."""
    src:          str                      # Node id
    dst:          str                      # Node id
    label:        str                      # human label for the link
    flow_fn:      Callable                 # fn(params) → actual avg flow in Bps
    dashed:       bool   = False           # draw dashed (control/metadata flows)
    label_frac:   float  = 0.5            # label position along edge (0=src, 1=dst)
    label_nudge:  float  = 1.0            # perpendicular nudge multiplier (neg = flip side)
    # computed
    flow_Bps:     float = field(default=0.0, init=False)

    def evaluate(self, params: dict, default_margin: float):
        self.flow_Bps = self.flow_fn(params)


@dataclass
class Annotation:
    """A free-floating text box on the flow diagram."""
    text_fn:  Callable          # fn(params) → str
    pos:      Tuple[float, float]
    color:    str = "#f85149"
    bg:       str = "#2a1010"
    fontsize: float = 7.5


# ══════════════════════════════════════════════════════════════
# GRAPH
# ══════════════════════════════════════════════════════════════

@dataclass
class FlowGraph:
    nodes:       List[Node]
    edges:       List[Edge]
    params:      dict
    annotations: List[Annotation] = field(default_factory=list)

    def evaluate(self):
        sm = self.params.get("safety_margin", 1.0)
        for n in self.nodes:
            n.evaluate(self.params, sm)
        for e in self.edges:
            e.evaluate(self.params, sm)

    # ── analysis helpers ──────────────────────────────────────

    def node_by_id(self, nid: str) -> Node:
        return next(n for n in self.nodes if n.id == nid)

    def all_constraints(self) -> List[Tuple[str, float, Optional[float]]]:
        """Return list of (label, utilization, max_trigger_rate_hz) for every
        node resource, sorted by descending utilization."""
        rows = []
        for n in self.nodes:
            for r in n.resources:
                rows.append((f"{n.label.replace(chr(10),' ')} — {r.name}",
                             r.utilization, r.max_trigger_rate_hz))
        return sorted(rows, key=lambda x: -x[1])

    def bottleneck(self) -> Tuple[str, float, Optional[float]]:
        """Binding constraint on trigger rate (ignores fixed-rate links)."""
        rate_constraints = [(l, u, r) for l, u, r in self.all_constraints()
                            if r is not None]
        if rate_constraints:
            return min(rate_constraints, key=lambda x: x[2])
        return self.all_constraints()[0]

    def system_max_rate(self) -> float:
        rates = [r for _, _, r in self.all_constraints() if r is not None]
        return min(rates) if rates else float("nan")

    def print_summary(self):
        p  = self.params
        sm = p.get("safety_margin", 1.0)

        def bps(b):
            if b >= 1e12: return f"{b/1e12:.3f} Tb/s"
            if b >= 1e9:  return f"{b/1e9:.3f} Gb/s"
            if b >= 1e6:  return f"{b/1e6:.2f} Mb/s"
            return f"{b/1e3:.2f} kb/s"
        def Bps(b):
            if b >= 1e9:  return f"{b/1e9:.3f} GB/s"
            if b >= 1e6:  return f"{b/1e6:.2f} MB/s"
            return f"{b/1e3:.2f} KB/s"

        W = 70
        print("=" * W)
        print("  DAQ FLOW GRAPH — STEADY-STATE ANALYSIS")
        print("=" * W)
        print(f"  Safety margin: {sm*100:.0f}% usable  ({(1-sm)*100:.0f}% reserved)")
        print(f"  Trigger rate:  {p['trigger_rate_hz']} Hz")
        print(f"  Event window:  {p['event_window_s']*1e3:.2f} ms")
        print(f"  BDE event/node:{_bde_event_per_node_B(p)/1e6:.1f} MB  "
              f"({Bps(_bde_ro_node_ingress_Bps(p))} × {p['event_window_s']*1e3:.2f} ms)")
        print(f"  TDE event/node:{_tde_event_per_node_B(p)/1e6:.1f} MB  "
              f"({Bps(_tde_ro_node_ingress_Bps(p))} × {p['event_window_s']*1e3:.2f} ms)")
        print(f"  Total event:   {_event_total_B(p)/1e6:.1f} MB")

        print(f"\n{'─'*W}")
        print(f"  {'EDGE FLOWS'}")
        print(f"  {'EDGE / LINK':<42}  {'FLOW':>12}")
        print(f"  {'─'*45}")
        for e in self.edges:
            if e.dashed:
                continue
            print(f"  {e.label.replace(chr(10), ' '):<42}  {bps(e.flow_Bps*8):>12}")

        print(f"\n{'─'*W}")
        print(f"  {'NODE RESOURCE':<34}  {'FLOW':>12}  {'EFF.CAP':>12}  {'UTIL':>7}")
        print(f"{'─'*W}")
        for n in self.nodes:
            for r in n.resources:
                flag = "  ◀ OVER" if r.utilization > 1.0 else ""
                print(f"  {n.label.replace(chr(10),' ')[:16]+' — '+r.name:<34}  "
                      f"{Bps(r.flow_Bps):>12}  {Bps(r.capacity_Bps):>12}  "
                      f"{r.utilization*100:>6.1f}%{flag}")

        print(f"\n{'─'*W}")
        print(f"  BOTTLENECK RANKING")
        print(f"{'─'*W}")
        for i, (lbl, util, mr) in enumerate(self.all_constraints()):
            tag  = "  ◀ #1 BOTTLENECK" if i == 0 else f"  #{i+1}"
            mstr = f"  max={mr:.4f} Hz" if mr is not None else ""
            print(f"  {lbl:<40}  {util*100:6.1f}%{mstr}{tag}")

        smr = self.system_max_rate()
        bn  = self.bottleneck()
        print(f"\n  System max trigger rate : {smr:.4f} Hz  (limited by: {bn[0]})")
        print(f"  Current trigger rate    : {p['trigger_rate_hz']} Hz")
        print(f"  Headroom factor         : {smr/p['trigger_rate_hz']:.2f}×")

        num_crates = p["num_bde_crates"] + p["num_tde_crates"]
        event_per_crate_B = _event_total_B(p) / num_crates
        headroom_hz = smr - p["trigger_rate_hz"]
        per_crate_headroom_hz = headroom_hz * _event_total_B(p) / event_per_crate_B
        print(f"  Trigger rate headroom   : {headroom_hz:.4f} Hz  ({per_crate_headroom_hz:.2f} Hz per-crate equivalent)")

        tape_limit_B = 30e15  # 30 PB
        seconds_per_year = 365.25 * 24 * 3600
        yearly_events_B = _filtered_event_B(p) * p["trigger_rate_hz"] * seconds_per_year
        yearly_tp_B     = _tp_total_Bps(p) * seconds_per_year
        yearly_B        = yearly_events_B + yearly_tp_B
        tape_fraction   = yearly_B / tape_limit_B
        print(f"\n  Yearly triggered events : {yearly_events_B/1e15:.2f} PB/yr")
        print(f"  Yearly TP stream        : {yearly_tp_B/1e15:.2f} PB/yr")
        print(f"  Yearly total to tape    : {yearly_B/1e15:.2f} PB/yr  ({tape_fraction*100:.1f}% of 30 PB/yr limit)")

        # ── Buffer memory requirements ────────────────────────────────────
        ro_buffer_s  = 15.0
        eb_buffer_s  = 600.0

        bde_ro_mem_B = _bde_ro_node_ingress_Bps(p) * ro_buffer_s
        tde_ro_mem_B = _tde_ro_node_ingress_Bps(p) * ro_buffer_s
        ro_mem_total_B = (bde_ro_mem_B * p["num_bde_ro_nodes"]
                          + tde_ro_mem_B * p["num_tde_ro_nodes"])

        eb_trig_Bps  = _event_total_B(p) * p["trigger_rate_hz"]
        eb_tp_Bps    = _tp_total_Bps(p)
        eb_mem_B     = (eb_trig_Bps + eb_tp_Bps) * eb_buffer_s / p["num_eb_nodes"]
        eb_mem_total_B = eb_mem_B * p["num_eb_nodes"]

        print(f"\n{'─'*W}")
        print(f"  MEMORY REQUIREMENTS")
        print(f"{'─'*W}")
        print(f"  RO node buffer  ({ro_buffer_s:.0f} s of raw detector data per node):")
        print(f"    BDE RO node   : {fmt_size(bde_ro_mem_B)}  "
              f"({Bps(_bde_ro_node_ingress_Bps(p))} × {ro_buffer_s:.0f} s)")
        print(f"    TDE RO node   : {fmt_size(tde_ro_mem_B)}  "
              f"({Bps(_tde_ro_node_ingress_Bps(p))} × {ro_buffer_s:.0f} s)")
        print(f"    All RO nodes  : {fmt_size(ro_mem_total_B)}  total")
        print(f"  EB node buffer  ({eb_buffer_s:.0f} s of triggered events + TP stream per node):")
        print(f"    Triggered data: {Bps(eb_trig_Bps)}  ({Bps(eb_trig_Bps/p['num_eb_nodes'])} per node)")
        print(f"    TP stream     : {Bps(eb_tp_Bps)}  ({Bps(eb_tp_Bps/p['num_eb_nodes'])} per node)")
        print(f"    Per EB node   : {fmt_size(eb_mem_B)}")
        print(f"    All EB nodes  : {fmt_size(eb_mem_total_B)}  total")

        # ── 100s continuous drain scenario ────────────────────────────────
        drain_window_s = 100.0

        def _drain_report(label, volume_B, path):
            path_eff = {k: v * sm for k, v in path.items()}
            bn_lbl   = min(path_eff, key=path_eff.get)
            throughput_Bps = path_eff[bn_lbl]
            drain_time_s   = volume_B / throughput_Bps
            print(f"\n{'─'*W}")
            print(f"  {label}  ({fmt_size(volume_B)} total)")
            print(f"{'─'*W}")
            for lbl, cap in path.items():
                flag = "  ◀ BOTTLENECK" if lbl == bn_lbl else ""
                print(f"  {lbl:<32}  eff. cap {Bps(cap*sm):>10}{flag}")
            print(f"\n  Drain time : {drain_time_s:.0f} s  "
                  f"({drain_time_s/3600:.2f} hr)  at {Bps(throughput_Bps)}")

        raw_drain_path = {
            "RO egress (agg)":      p["ro_egress_bps"]         / 8 * _num_ro_nodes(p),
            "EB NIC ingress (agg)": p["eb_ingress_bps"]        / 8 * p["num_eb_nodes"],
            "EB → local storage":   p["eb_to_storage_Bps"]         * p["num_eb_nodes"],
            "Local storage write":  p["local_storage_write_Bps"]   * p["num_eb_nodes"],
            "Storage → offline":    p["offline_bps"]           / 8,
        }
        _drain_report(
            f"100s DRAIN — raw detector data (no filter)",
            _total_detector_Bps(p) * drain_window_s,
            raw_drain_path,
        )

        tp_drain_path = {
            "RO egress (agg)":      p["ro_egress_bps"]         / 8 * _num_ro_nodes(p),
            "EB NIC ingress (agg)": p["eb_ingress_bps"]        / 8 * p["num_eb_nodes"],
            "Local storage write":  p["local_storage_write_Bps"]   * p["num_eb_nodes"],
            "Storage → offline":    p["offline_bps"]           / 8,
        }
        _drain_report(
            f"100s DRAIN — TP stream only",
            _tp_total_Bps(p) * drain_window_s,
            tp_drain_path,
        )

        print("=" * W)


# ══════════════════════════════════════════════════════════════
# VISUALISATION
# ══════════════════════════════════════════════════════════════

BG, PANEL, BORDER = "#0d1117", "#161b22", "#30363d"
TXTMAIN, TXTSUB   = "#e6edf3", "#8b949e"

def util_color(u: float) -> tuple:
    u = float(np.clip(u, 0, 1))
    r = u * 2 if u < 0.5 else 1.0
    g = 1.0   if u < 0.5 else (1.0 - u) * 2
    return (0.12 + 0.85*r, 0.12 + 0.85*g, 0.12)

def fmt_bps(b: float) -> str:
    if b >= 1e12: return f"{b/1e12:.2f} Tb/s"
    if b >= 1e9:  return f"{b/1e9:.2f} Gb/s"
    if b >= 1e6:  return f"{b/1e6:.1f} Mb/s"
    return f"{b/1e3:.1f} kb/s"

def fmt_Bps(b: float) -> str:
    if b >= 1e9:  return f"{b/1e9:.2f} GB/s"
    if b >= 1e6:  return f"{b/1e6:.1f} MB/s"
    return f"{b/1e3:.1f} KB/s"

def fmt_size(b: float) -> str:
    if b >= 1e12: return f"{b/1e12:.2f} TB"
    if b >= 1e9:  return f"{b/1e9:.2f} GB"
    if b >= 1e6:  return f"{b/1e6:.1f} MB"
    return f"{b/1e3:.1f} KB"


def draw_graph(graph: FlowGraph, out_path: str = "daq_flow_diagram.png"):
    graph.evaluate()
    p  = graph.params
    sm = p.get("safety_margin", 1.0)

    fig = plt.figure(figsize=(22, 15), facecolor=BG)
    fig.suptitle("Physics Detector DAQ  —  Data Flow & Bottleneck Model",
                 fontsize=16, color=TXTMAIN, fontweight="bold", y=0.988,
                 fontfamily="monospace")

    gs = GridSpec(2, 2, figure=fig,
                  hspace=0.44, wspace=0.30,
                  left=0.03, right=0.97, top=0.95, bottom=0.05)
    ax_flow = fig.add_subplot(gs[0, :])
    ax_util = fig.add_subplot(gs[1, 0])
    ax_rate = fig.add_subplot(gs[1, 1])

    for ax in (ax_flow, ax_util, ax_rate):
        ax.set_facecolor(PANEL)
        for sp in ax.spines.values():
            sp.set_color(BORDER); sp.set_linewidth(0.8)

    # ── node id → pixel position lookup ─────────────────────
    node_pos = {n.id: n.pos for n in graph.nodes}

    # ── bounding box for axis limits ────────────────────────
    xs = [p[0] for p in node_pos.values()]
    ys = [p[1] for p in node_pos.values()]
    pad = 1.5
    ax_flow.set_xlim(min(xs)-pad, max(xs)+pad)
    ax_flow.set_ylim(min(ys)-pad*0.6, max(ys)+pad*1.2)
    ax_flow.axis("off")
    ax_flow.set_title(
        f"Data Flow Topology  —  safety margin {sm*100:.0f}%"
        f"  (utilization shown vs. effective capacity)",
        color=TXTSUB, fontsize=10, pad=4, loc="left")

    # ── draw edges ───────────────────────────────────────────
    for e in graph.edges:
        x1, y1 = node_pos[e.src]
        x2, y2 = node_pos[e.dst]
        col = "#4d9de0" if not e.dashed else "#7c5cdb"
        lw  = 1.2 if e.dashed else 2.0
        ax_flow.annotate("",
            xy=(x2, y2), xytext=(x1, y1),
            arrowprops=dict(arrowstyle="-|>", color=col, lw=lw,
                            linestyle="dashed" if e.dashed else "solid"),
            zorder=2)
        if e.label and not e.dashed:
            mx = x1 + e.label_frac * (x2 - x1)
            my = y1 + e.label_frac * (y2 - y1)
            dx, dy = x2-x1, y2-y1
            length = math.hypot(dx, dy) or 1
            nx = -dy/length * 0.32 * e.label_nudge
            ny =  dx/length * 0.32 * e.label_nudge
            ax_flow.text(mx+nx, my+ny, f"{e.label}\n{fmt_bps(e.flow_Bps*8)}",
                ha="center", va="center", color=col,
                fontsize=6.8, fontweight="bold", zorder=5,
                bbox=dict(facecolor=PANEL, edgecolor="none", pad=1.2, alpha=0.88))

    # ── draw nodes ───────────────────────────────────────────
    for n in graph.nodes:
        x, y  = n.pos
        w, h  = 2.1, 1.0
        worst = n.worst_resource
        ec    = "#f85149" if (worst and worst.utilization > 1.0) else "#58a6ff"
        ax_flow.add_patch(mpatches.FancyBboxPatch(
            (x-w/2, y-h/2), w, h,
            boxstyle="round,pad=0.08", facecolor=n.color,
            edgecolor=ec, linewidth=1.4, zorder=3, alpha=0.93))
        ax_flow.text(x, y+h*0.14, n.label,
            ha="center", va="center",
            color=TXTMAIN, fontsize=8.0, fontweight="bold", zorder=4)
        sub = n.sublabel_fn(p)
        ax_flow.text(x, y-h*0.24, sub,
            ha="center", va="center",
            color=TXTSUB, fontsize=6.5, zorder=4, linespacing=1.3)
        # worst-resource utilization badge
        if worst:
            badge_col = util_color(min(worst.utilization, 1.0))
            ax_flow.text(x+w/2-0.05, y+h/2-0.05,
                f"{worst.utilization*100:.0f}%",
                ha="right", va="top", color=badge_col,
                fontsize=6.5, fontweight="bold", zorder=6)

    # ── annotations ─────────────────────────────────────────
    for ann in graph.annotations:
        ax_flow.text(ann.pos[0], ann.pos[1], ann.text_fn(p),
            ha="center", va="bottom", color=ann.color,
            fontsize=ann.fontsize, fontweight="bold", zorder=6,
            bbox=dict(facecolor=ann.bg, edgecolor=ann.color,
                      boxstyle="round,pad=0.38", linewidth=1.1))

    # ── utilization bar chart ────────────────────────────────
    ax = ax_util
    ax.set_title(
        f"Node Resource Utilization  @  {p['trigger_rate_hz']} Hz"
        f"  |  safety margin {sm*100:.0f}%",
        color=TXTSUB, fontsize=9.5, pad=6)

    bar_items = []
    for n in graph.nodes:
        for r in n.resources:
            bar_items.append((
                f"{n.label.replace(chr(10),' ')[:14]} — {r.name}",
                r.utilization))
    bar_items.sort(key=lambda x: -x[1])

    yi   = np.arange(len(bar_items))
    vals = [v for _, v in bar_items]
    cols = [util_color(min(v, 1.0)) for v in vals]

    ax.barh(yi, [min(v, 1.35)*100 for v in vals],
            color=cols, height=0.6, zorder=3)
    # hatch: headroom to raw hardware limit
    for i, v in enumerate(vals):
        raw_pct = min(v, 1.0) * 100 / sm
        ax.barh(i, raw_pct - min(v*100, 100),
                left=min(v*100, 100),
                color="none", edgecolor="#444", linewidth=0.6,
                height=0.6, zorder=2, hatch="////")

    ax.set_yticks(yi)
    ax.set_yticklabels([lbl for lbl, _ in bar_items],
                       fontsize=8.0, color=TXTMAIN)
    ax.set_xlabel("Utilization (% of effective capacity)",
                  color=TXTSUB, fontsize=9)
    ax.tick_params(colors=TXTSUB)
    ax.set_xlim(0, 148)
    ax.axvline(100, color="#f85149", lw=1.5, ls="--", alpha=0.8, zorder=4)
    ax.text(101, -0.8, "eff. cap", color="#f85149", fontsize=7)
    ax.grid(axis="x", color=BORDER, alpha=0.5, zorder=0)
    for i, v in enumerate(vals):
        label = f"{v*100:.1f}%" + ("  ⚠" if v > 1.0 else "")
        ax.text(min(v*100, 134) + 1.0, i,
                label, va="center",
                color="#f85149" if v > 1.0 else TXTMAIN, fontsize=8)

    # ── max rate chart ───────────────────────────────────────
    ax = ax_rate
    ax.set_title(
        "Max Sustainable Trigger Rate per Constraint\n(effective capacities)",
        color=TXTSUB, fontsize=9.5, pad=6)

    rate_items = [(lbl, mr) for lbl, _, mr in graph.all_constraints()
                  if mr is not None]
    rate_items.sort(key=lambda x: x[1])

    if rate_items:
        yi2   = np.arange(len(rate_items))
        rates = [r for _, r in rate_items]
        rcols = [util_color(p["trigger_rate_hz"] / r) for r in rates]
        ax.barh(yi2, rates, color=rcols, height=0.6, zorder=3)
        ax.set_yticks(yi2)
        ax.set_yticklabels([lbl for lbl, _ in rate_items],
                           fontsize=8.0, color=TXTMAIN)
        ax.set_xlabel("Max trigger rate  (Hz)", color=TXTSUB, fontsize=9)
        ax.tick_params(colors=TXTSUB)
        ax.grid(axis="x", color=BORDER, alpha=0.5, zorder=0)
        ax.axvline(p["trigger_rate_hz"], color="#58a6ff", lw=2.0,
                   ls="--", zorder=5)
        ax.text(p["trigger_rate_hz"]*1.04, -0.7,
                f"Current\n{p['trigger_rate_hz']} Hz",
                color="#58a6ff", fontsize=7.5)
        for i, r in enumerate(rates):
            ax.text(r*1.02, i, f"{r:.4f} Hz",
                    va="center", color=TXTMAIN, fontsize=8.5)

    # shared colorbar
    from matplotlib.cm import ScalarMappable
    from matplotlib.colors import LinearSegmentedColormap
    cmap = LinearSegmentedColormap.from_list("u",
               ["#22c55e", "#f59e0b", "#ef4444"])
    sm_obj = ScalarMappable(cmap=cmap, norm=plt.Normalize(0, 100))
    sm_obj.set_array([])
    cbar = fig.colorbar(sm_obj, ax=[ax_util, ax_rate],
                        orientation="vertical",
                        fraction=0.014, pad=0.01, shrink=0.65)
    cbar.set_label("Utilization %", color=TXTSUB, fontsize=8)
    cbar.ax.yaxis.set_tick_params(color=TXTSUB)
    plt.setp(cbar.ax.yaxis.get_ticklabels(), color=TXTSUB, fontsize=7)

    plt.savefig(out_path, dpi=150, bbox_inches="tight", facecolor=BG)
    print(f"Figure saved → {out_path}")
    plt.close()


# ══════════════════════════════════════════════════════════════
# DAQ TOPOLOGY DEFINITION
# ══════════════════════════════════════════════════════════════
# Everything below here is the *description* of the specific DAQ
# system.  The Node/Edge/Annotation classes above are generic and
# can be reused for any flow system.
# ══════════════════════════════════════════════════════════════

# ── Global parameters ─────────────────────────────────────────
PARAMS = {
    # BDE front-end (12 × 10G links/crate, 4 streams/link)
    "bde_frame_bytes":          7_200,
    "bde_frame_interval_ns":    2048 * 16,      # 32 768 ns
    "bde_links_per_crate":      12,
    "bde_streams_per_link":     4,
    "bde_link_bps":             10e9,
    "num_bde_crates":           80,

    # TDE front-end (4 × 25G links/crate, 12 streams/link)
    "tde_frame_bytes":          7_232,
    "tde_frame_interval_ns":    2000 * 16,      # 32 000 ns
    "tde_links_per_crate":      4,
    "tde_streams_per_link":     12,
    "tde_link_bps":             25e9,
    "num_tde_crates":           80,

    "num_bde_switches":         10,
    "num_tde_switches":         10,
    "num_bde_ro_nodes":         20,
    "num_tde_ro_nodes":         20,

    # Links (raw hardware, bits/s)
    "ro_switch_bps":            400e9,          # switch → RO node bundle
    "ro_egress_bps":            25e9,
    "eb_ingress_bps":           100e9,
    "offline_bps":              100e9,

    # Trigger primitives (50 MHz × 24 B, 1 ms buckets, even across all RO nodes)
    "tp_rate_hz":               50e6,
    "tp_size_bytes":            24,
    "tp_bucket_ms":             1,

    # Event selection
    "trigger_rate_hz":          0.1,
    "event_window_s":           8.5e-3,

    # Event builders
    "num_eb_nodes":             6,
    "eb_to_storage_Bps":        10e9,           # per node, NVMe-like handoff (not a bottleneck)

    # Data filter cluster
    "num_filter_nodes":         10,
    "filter_link_bps":          25e9,           # per filter node (ingress and egress)
    "filter_reduction":         2.0,            # data reduction factor
    "filter_latency_s":         60.0,           # processing time per trigger record per core
    "filter_cores_per_node":    24,             # parallel cores per filter node

    "local_storage_write_Bps":  1000e6,          # per node, sustained disk write rate

    # Safety margin (applied to ALL edges/resources unless overridden)
    "safety_margin":            0.80,
}

# ── Shorthand helpers used in lambdas ─────────────────────────
def _bde_stream_Bps(p):
    return p["bde_frame_bytes"] / (p["bde_frame_interval_ns"] * 1e-9)

def _tde_stream_Bps(p):
    return p["tde_frame_bytes"] / (p["tde_frame_interval_ns"] * 1e-9)

def _bde_streams_per_crate(p):
    return p["bde_links_per_crate"] * p["bde_streams_per_link"]

def _tde_streams_per_crate(p):
    return p["tde_links_per_crate"] * p["tde_streams_per_link"]

def _bde_crate_Bps(p):
    return _bde_stream_Bps(p) * _bde_streams_per_crate(p)

def _tde_crate_Bps(p):
    return _tde_stream_Bps(p) * _tde_streams_per_crate(p)

def _bde_total_Bps(p):
    return _bde_crate_Bps(p) * p["num_bde_crates"]

def _tde_total_Bps(p):
    return _tde_crate_Bps(p) * p["num_tde_crates"]

def _total_detector_Bps(p):
    return _bde_total_Bps(p) + _tde_total_Bps(p)

def _num_ro_nodes(p):
    return p["num_bde_ro_nodes"] + p["num_tde_ro_nodes"]

def _bde_ro_node_ingress_Bps(p):
    return _bde_total_Bps(p) / p["num_bde_ro_nodes"]

def _tde_ro_node_ingress_Bps(p):
    return _tde_total_Bps(p) / p["num_tde_ro_nodes"]

def _bde_event_per_node_B(p):
    return _bde_ro_node_ingress_Bps(p) * p["event_window_s"]

def _tde_event_per_node_B(p):
    return _tde_ro_node_ingress_Bps(p) * p["event_window_s"]

def _event_total_B(p):
    return (_bde_total_Bps(p) + _tde_total_Bps(p)) * p["event_window_s"]

def _tp_total_Bps(p):
    return p["tp_rate_hz"] * p["tp_size_bytes"]

def _tp_per_node_Bps(p):
    return _tp_total_Bps(p) / _num_ro_nodes(p)

def _filtered_event_B(p):
    return _event_total_B(p) / p["filter_reduction"]

def _bde_ro_trigger_avg_Bps(p):
    return _bde_event_per_node_B(p) * p["trigger_rate_hz"]

def _tde_ro_trigger_avg_Bps(p):
    return _tde_event_per_node_B(p) * p["trigger_rate_hz"]

def _eb_receive_time_s(p, sm=1.0):
    return _event_total_B(p) / (p["eb_ingress_bps"] * sm / 8)

# ── Node definitions ──────────────────────────────────────────
NODES = [
    Node(
        id    = "bde",
        label = "BDE\nFront-ends",
        sublabel_fn = lambda p: (
            f"{p['num_bde_crates']} crates  |  "
            f"{p['bde_links_per_crate']}×10G / {p['bde_streams_per_link']} str/link\n"
            f"{fmt_bps(_bde_crate_Bps(p)*8)}/crate  →  {fmt_bps(_bde_total_Bps(p)*8)} total"
        ),
        color = "#0e3320",
        multiplicity_fn = lambda p: p["num_bde_crates"],
        pos   = (1.2, 2.8),
        resources = [
            Resource(
                name        = "link egress",
                capacity_fn = lambda p: p["bde_link_bps"] * p["bde_links_per_crate"] * p["num_bde_crates"] / 8,
                flow_fn     = lambda p: _bde_total_Bps(p),
                safety_margin = 1.0,
                max_rate_fn = None,
            ),
        ],
    ),
    Node(
        id    = "tde",
        label = "TDE\nFront-ends",
        sublabel_fn = lambda p: (
            f"{p['num_tde_crates']} crates  |  "
            f"{p['tde_links_per_crate']}×25G / {p['tde_streams_per_link']} str/link\n"
            f"{fmt_bps(_tde_crate_Bps(p)*8)}/crate  →  {fmt_bps(_tde_total_Bps(p)*8)} total"
        ),
        color = "#0e2a3a",
        multiplicity_fn = lambda p: p["num_tde_crates"],
        pos   = (1.2, 1.2),
        resources = [
            Resource(
                name        = "link egress",
                capacity_fn = lambda p: p["tde_link_bps"] * p["tde_links_per_crate"] * p["num_tde_crates"] / 8,
                flow_fn     = lambda p: _tde_total_Bps(p),
                safety_margin = 1.0,
                max_rate_fn = None,
            ),
        ],
    ),
    Node(
        id    = "bde_switch",
        label = "BDE Readout\nSwitches",
        sublabel_fn = lambda p: f"×{p['num_bde_switches']}  |  aggregation",
        color = "#1c2128",
        multiplicity_fn = lambda p: p["num_bde_switches"],
        pos   = (3.6, 2.8),
    ),
    Node(
        id    = "tde_switch",
        label = "TDE Readout\nSwitches",
        sublabel_fn = lambda p: f"×{p['num_tde_switches']}  |  aggregation",
        color = "#1c2128",
        multiplicity_fn = lambda p: p["num_tde_switches"],
        pos   = (3.6, 1.2),
    ),
    Node(
        id    = "bde_ro_nodes",
        label = "BDE Readout\nNodes",
        sublabel_fn = lambda p: (
            f"×{p['num_bde_ro_nodes']}  |  400G in / 25G out\n"
            f"buffer + TP processing"
        ),
        color = "#0c2a52",
        multiplicity_fn = lambda p: p["num_bde_ro_nodes"],
        pos   = (6.2, 2.8),
        resources = [
            Resource(
                name        = "switch ingress",
                capacity_fn = lambda p: p["ro_switch_bps"] / 8 * p["num_bde_ro_nodes"],
                flow_fn     = lambda p: _bde_total_Bps(p),
                safety_margin = 1.0,
                max_rate_fn = None,
            ),
            Resource(
                name        = "NIC egress",
                capacity_fn = lambda p: p["ro_egress_bps"] / 8 * p["num_bde_ro_nodes"],
                flow_fn     = lambda p: (_bde_ro_trigger_avg_Bps(p) + 2 * _tp_per_node_Bps(p)) * p["num_bde_ro_nodes"],
                max_rate_fn = lambda p, sm: (
                    (p["ro_egress_bps"] * sm / 8 - 2 * _tp_per_node_Bps(p))
                    / _bde_event_per_node_B(p)
                ),
            ),
            Resource(
                name        = "buffer (net fill)",
                capacity_fn = lambda p: float("inf"),
                flow_fn     = lambda p: _bde_ro_node_ingress_Bps(p) - p["ro_egress_bps"] / 8,
                max_rate_fn = None,
            ),
        ],
    ),
    Node(
        id    = "tde_ro_nodes",
        label = "TDE Readout\nNodes",
        sublabel_fn = lambda p: (
            f"×{p['num_tde_ro_nodes']}  |  400G in / 25G out\n"
            f"buffer + TP processing"
        ),
        color = "#0c2a52",
        multiplicity_fn = lambda p: p["num_tde_ro_nodes"],
        pos   = (6.2, 1.2),
        resources = [
            Resource(
                name        = "switch ingress",
                capacity_fn = lambda p: p["ro_switch_bps"] / 8 * p["num_tde_ro_nodes"],
                flow_fn     = lambda p: _tde_total_Bps(p),
                safety_margin = 1.0,
                max_rate_fn = None,
            ),
            Resource(
                name        = "NIC egress",
                capacity_fn = lambda p: p["ro_egress_bps"] / 8 * p["num_tde_ro_nodes"],
                flow_fn     = lambda p: (_tde_ro_trigger_avg_Bps(p) + 2 * _tp_per_node_Bps(p)) * p["num_tde_ro_nodes"],
                max_rate_fn = lambda p, sm: (
                    (p["ro_egress_bps"] * sm / 8 - 2 * _tp_per_node_Bps(p))
                    / _tde_event_per_node_B(p)
                ),
            ),
            Resource(
                name        = "buffer (net fill)",
                capacity_fn = lambda p: float("inf"),
                flow_fn     = lambda p: _tde_ro_node_ingress_Bps(p) - p["ro_egress_bps"] / 8,
                max_rate_fn = None,
            ),
        ],
    ),
    Node(
        id    = "trigger_nodes",
        label = "Trigger\nNodes",
        sublabel_fn = lambda p: (
            f"data selection\n"
            f"{fmt_bps(_tp_total_Bps(p)*8)} TP in"
        ),
        color = "#2b1969",
        pos   = (9.0, 3.4),
    ),
    Node(
        id    = "eb_nodes",
        label = "Event Builders",
        sublabel_fn = lambda p: (
            f"×{p['num_eb_nodes']}  |  100G in\n"
            f"assemble + pass to storage"
        ),
        color = "#4a1800",
        multiplicity_fn = lambda p: p["num_eb_nodes"],
        pos   = (9.0, 0.7),
        resources = [
            Resource(
                name        = "NIC ingress",
                capacity_fn = lambda p: p["eb_ingress_bps"] * p["num_eb_nodes"] / 8,
                flow_fn     = lambda p: _event_total_B(p) * p["trigger_rate_hz"] + _tp_total_Bps(p),
                max_rate_fn = lambda p, sm: p["num_eb_nodes"] / _eb_receive_time_s(p, sm),
            ),
        ],
    ),
    Node(
        id    = "data_filter",
        label = "Data Filter\nCluster",
        sublabel_fn = lambda p: (
            f"×{p['num_filter_nodes']}  |  25G links\n"
            f"reduction ×{p['filter_reduction']:.0f}  |  {p['filter_cores_per_node']} cores × {p['filter_latency_s']:.0f}s/record"
        ),
        color = "#1a1040",
        multiplicity_fn = lambda p: p["num_filter_nodes"],
        pos   = (12.1, 2.0),
        resources = [
            Resource(
                name        = "filter NIC ingress",
                capacity_fn = lambda p: p["filter_link_bps"] * p["num_filter_nodes"] / 8,
                flow_fn     = lambda p: _event_total_B(p) * p["trigger_rate_hz"],
                max_rate_fn = lambda p, sm: (
                    p["filter_link_bps"] * sm * p["num_filter_nodes"] / 8 / _event_total_B(p)
                ),
            ),
            Resource(
                name        = "NIC egress",
                capacity_fn = lambda p: p["filter_link_bps"] * p["num_filter_nodes"] / 8,
                flow_fn     = lambda p: _filtered_event_B(p) * p["trigger_rate_hz"],
                max_rate_fn = lambda p, sm: (
                    p["filter_link_bps"] * sm * p["num_filter_nodes"] / 8 / _filtered_event_B(p)
                ),
            ),
            Resource(
                name        = "processing throughput",
                # capacity expressed as equivalent Bps: max records/s × event size
                # each node has filter_cores_per_node parallel cores, each handling one record
                capacity_fn = lambda p: (
                    p["num_filter_nodes"] * p["filter_cores_per_node"]
                    / p["filter_latency_s"] * _event_total_B(p)
                ),
                flow_fn     = lambda p: _event_total_B(p) * p["trigger_rate_hz"],
                max_rate_fn = lambda p, sm: (
                    p["num_filter_nodes"] * p["filter_cores_per_node"] * sm / p["filter_latency_s"]
                ),
            ),
        ],
    ),
    Node(
        id    = "local_storage",
        label = "Local Storage\n(EB disks)",
        sublabel_fn = lambda p: (
            f"×{p['num_eb_nodes']}  |  {fmt_Bps(p['local_storage_write_Bps'])}/s write\n"
            f"= {fmt_Bps(p['local_storage_write_Bps']*p['num_eb_nodes'])} total"
        ),
        color = "#272000",
        pos   = (12.1, 0.7),
        resources = [
            Resource(
                name        = "disk write",
                capacity_fn = lambda p: p["local_storage_write_Bps"] * p["num_eb_nodes"],
                # disk write serves both filtered events (bursty) and continuous TP archive
                flow_fn     = lambda p: (
                    _filtered_event_B(p) * p["trigger_rate_hz"] + _tp_total_Bps(p)
                ),
                max_rate_fn = lambda p, sm: (
                    (p["local_storage_write_Bps"] * sm * p["num_eb_nodes"] - _tp_total_Bps(p))
                    / _filtered_event_B(p)
                ),
            ),
        ],
    ),
    Node(
        id    = "offline",
        label = "Offline /\nTape Staging",
        sublabel_fn = lambda p: "100G link",
        color = "#0e2a1a",
        pos   = (15.0, 0.7),
        resources = [
            Resource(
                name        = "link ingress",
                capacity_fn = lambda p: p["offline_bps"] / 8,
                flow_fn     = lambda p: _filtered_event_B(p) * p["trigger_rate_hz"] + _tp_total_Bps(p),
                max_rate_fn = lambda p, sm: (
                    (p["offline_bps"] * sm / 8 - _tp_total_Bps(p)) / _filtered_event_B(p)
                ),
            ),
        ],
    ),
]

# ── Edge definitions ──────────────────────────────────────────
EDGES = [
    Edge(
        src   = "bde",
        dst   = "bde_switch",
        label = "BDE → Switch",
        flow_fn = lambda p: _bde_total_Bps(p),
    ),
    Edge(
        src   = "tde",
        dst   = "tde_switch",
        label = "TDE → Switch",
        flow_fn = lambda p: _tde_total_Bps(p),
    ),
    Edge(
        src   = "bde_switch",
        dst   = "bde_ro_nodes",
        label = "BDE Switch → RO\n400G/node",
        flow_fn = lambda p: _bde_total_Bps(p),
    ),
    Edge(
        src   = "tde_switch",
        dst   = "tde_ro_nodes",
        label = "TDE Switch → RO\n400G/node",
        flow_fn = lambda p: _tde_total_Bps(p),
    ),
    Edge(
        src       = "bde_ro_nodes",
        dst       = "trigger_nodes",
        label     = "BDE RO → Trig\nTP primitives",
        flow_fn   = lambda p: _tp_per_node_Bps(p) * p["num_bde_ro_nodes"],
        label_frac = 0.35,
    ),
    Edge(
        src       = "tde_ro_nodes",
        dst       = "trigger_nodes",
        label     = "TDE RO → Trig\nTP primitives",
        flow_fn   = lambda p: _tp_per_node_Bps(p) * p["num_tde_ro_nodes"],
        label_frac = 0.65,
    ),
    Edge(
        src       = "bde_ro_nodes",
        dst       = "eb_nodes",
        label     = "BDE RO → EB\ntriggered readout",
        flow_fn   = lambda p: _bde_ro_trigger_avg_Bps(p) * p["num_bde_ro_nodes"],
        label_frac = 0.25,
    ),
    Edge(
        src        = "bde_ro_nodes",
        dst        = "eb_nodes",
        label      = "BDE TP → EB\nfor recording",
        flow_fn    = lambda p: _tp_total_Bps(p) / 2,
        label_frac  = 0.55,
        label_nudge = -1.5,
    ),
    Edge(
        src       = "tde_ro_nodes",
        dst       = "eb_nodes",
        label     = "TDE RO → EB\ntriggered readout",
        flow_fn   = lambda p: _tde_ro_trigger_avg_Bps(p) * p["num_tde_ro_nodes"],
        label_frac = 0.25,
    ),
    Edge(
        src        = "tde_ro_nodes",
        dst        = "eb_nodes",
        label      = "TDE TP → EB\nfor recording",
        flow_fn    = lambda p: _tp_total_Bps(p) / 2,
        label_frac  = 0.55,
        label_nudge = -1.5,
    ),
    Edge(
        src       = "eb_nodes",
        dst       = "data_filter",
        label     = "EB → Filter\nfull events",
        flow_fn   = lambda p: _event_total_B(p) * p["trigger_rate_hz"],
        label_frac = 0.45,
    ),
    Edge(
        src        = "data_filter",
        dst        = "eb_nodes",
        label      = "Filter → EB\nfiltered events",
        flow_fn    = lambda p: _filtered_event_B(p) * p["trigger_rate_hz"],
        label_frac  = 0.55,
        label_nudge = -1.5,
    ),
    Edge(
        src       = "eb_nodes",
        dst       = "local_storage",
        label     = "EB → Storage\nfiltered handoff",
        flow_fn   = lambda p: _filtered_event_B(p) * p["trigger_rate_hz"],
        label_frac = 0.2,
    ),
    Edge(
        src       = "eb_nodes",
        dst       = "local_storage",
        label     = "EB → Storage\nTPStream",
        flow_fn   = lambda p: _tp_total_Bps(p),
        label_frac = 0.8,
    ),
    Edge(
        src     = "local_storage",
        dst     = "offline",
        label   = "Storage → Offline\n100G link",
        flow_fn = lambda p: _filtered_event_B(p) * p["trigger_rate_hz"] + _tp_total_Bps(p),
    ),
]

# ── Annotations ───────────────────────────────────────────────
ANNOTATIONS = [
    Annotation(
        text_fn = lambda p: (
            f"⚡  {fmt_size(_event_total_B(p))} per event  →  1 EB  "
            f"over {_eb_receive_time_s(p, p['safety_margin']):.2f}s  "
            f"(eff. {fmt_bps(p['eb_ingress_bps']*p['safety_margin'])} limited)\n"
            f"{_num_ro_nodes(p)} nodes × 25G = {_num_ro_nodes(p)*25//10/10:.0f} Tb/s offered  →  EB NIC is drain bottleneck"
        ),
        pos   = (9.0, 0.0),
        color = "#f85149",
        bg    = "#2a1010",
    ),
    Annotation(
        text_fn = lambda p: (
            f"⚠  BDE buffer grows at "
            f"{fmt_Bps(_bde_ro_node_ingress_Bps(p) - p['ro_egress_bps']/8)}/node"
        ),
        pos   = (6.2, 3.8),
        color = "#e3b341",
        bg    = "#261d00",
        fontsize = 8.0,
    ),
    Annotation(
        text_fn = lambda p: (
            f"⚠  TDE buffer grows at "
            f"{fmt_Bps(_tde_ro_node_ingress_Bps(p) - p['ro_egress_bps']/8)}/node"
        ),
        pos   = (6.2, 0.2),
        color = "#e3b341",
        bg    = "#261d00",
        fontsize = 8.0,
    ),
]

# ══════════════════════════════════════════════════════════════
# DISCRETE-EVENT SIMULATION
# ══════════════════════════════════════════════════════════════

def _simulate_stage(job_ready_times: np.ndarray,
                    n_servers: int,
                    svc_time: float):
    """
    Simulate a G/D/c queue (general arrivals, deterministic service,
    c parallel servers).

    Parameters
    ----------
    job_ready_times : sorted array of times at which each job becomes
                      available to this stage (output of the previous stage,
                      or trigger arrival times for the first stage).
    n_servers       : number of parallel servers (cores / nodes).
    svc_time        : deterministic service time per job (seconds).

    Returns
    -------
    start_times  : when each job enters service
    end_times    : when each job leaves this stage
    wait_times   : start_times - job_ready_times  (queueing delay)
    queue_depths : number of jobs waiting (not yet in service) when
                   each job arrives at this stage
    """
    n = len(job_ready_times)
    server_free_at = np.zeros(n_servers)
    start_times    = np.empty(n)
    end_times      = np.empty(n)

    for i, ready in enumerate(job_ready_times):
        s              = int(np.argmin(server_free_at))
        start          = max(ready, server_free_at[s])
        server_free_at[s] = start + svc_time
        start_times[i] = start
        end_times[i]   = start + svc_time

    wait_times = start_times - job_ready_times

    # Queue depth when job i arrives = jobs that arrived before i
    # but have not yet started service (still waiting in the queue).
    queue_depths = np.array([
        int(np.sum(start_times[:i] > job_ready_times[i]))
        for i in range(n)
    ], dtype=int)

    return start_times, end_times, wait_times, queue_depths


def simulate_flow(graph: "FlowGraph",
                  n_triggers: int = 10_000,
                  seed: int = 42) -> dict:
    """
    Run a discrete-event simulation of the DAQ trigger pipeline.

    Triggers arrive as a Poisson process.  Each trigger is treated as
    a job that passes sequentially through:
        1. EB assembly        (num_eb_nodes parallel servers)
        2. Filter processing  (num_filter_nodes × filter_cores_per_node servers)
        3. Storage write      (num_eb_nodes parallel write channels)

    Returns a dict with per-stage statistics and end-to-end latency.
    """
    graph.evaluate()
    p  = graph.params
    sm = p.get("safety_margin", 1.0)
    rng = np.random.default_rng(seed)

    # Poisson trigger arrivals
    inter_arrivals = rng.exponential(1.0 / p["trigger_rate_hz"], n_triggers)
    arrival_times  = np.cumsum(inter_arrivals)

    stages = [
        dict(
            name       = "EB assembly",
            n_servers  = p["num_eb_nodes"],
            svc_time   = _eb_receive_time_s(p, sm),
            job_size_B = _event_total_B(p),
        ),
        dict(
            name       = "Filter processing",
            n_servers  = p["num_filter_nodes"] * p["filter_cores_per_node"],
            svc_time   = p["filter_latency_s"],
            job_size_B = _event_total_B(p),
        ),
        dict(
            name       = "Storage write",
            n_servers  = p["num_eb_nodes"],
            svc_time   = _filtered_event_B(p) / p["local_storage_write_Bps"],
            job_size_B = _filtered_event_B(p),
        ),
    ]

    results      = {"stages": {}, "arrival_times": arrival_times}
    job_times    = arrival_times.copy()

    for st in stages:
        starts, ends, waits, qdepths = _simulate_stage(
            job_times, st["n_servers"], st["svc_time"])
        results["stages"][st["name"]] = {
            **st,
            "start_times":  starts,
            "end_times":    ends,
            "wait_times":   waits,
            "queue_depths": qdepths,
        }
        job_times = ends

    results["total_latency"] = job_times - arrival_times
    return results


def print_simulation_summary(results: dict, params: dict):
    p  = params
    W  = 70
    n  = len(results["arrival_times"])
    sim_duration_s = results["arrival_times"][-1]

    print("=" * W)
    print(f"  DISCRETE-EVENT SIMULATION")
    print(f"  {n} triggers  |  Poisson @ {p['trigger_rate_hz']} Hz  |  "
          f"{sim_duration_s/3600:.1f} hr simulated")
    print("=" * W)

    pcts = [50, 95, 99, 99.9]

    for name, r in results["stages"].items():
        waits   = r["wait_times"]
        qdepths = r["queue_depths"]
        svc     = r["svc_time"]
        n_srv   = r["n_servers"]
        job_B   = r["job_size_B"]
        rho     = p["trigger_rate_hz"] * svc / n_srv  # per-server utilisation

        w_pct = np.percentile(waits,   pcts)
        q_pct = np.percentile(qdepths, pcts)

        print(f"\n  {name}")
        print(f"  {'─'*60}")
        print(f"  Servers: {n_srv}   service time: {svc:.3f} s   "
              f"ρ (per-server util): {rho:.4f}")
        print(f"  Wait time (s)  : "
              + "  ".join(f"p{q:.0f}={v:.2f}" for q, v in zip(pcts, w_pct)))
        print(f"  Queue depth    : "
              + "  ".join(f"p{q:.0f}={v:.0f}" for q, v in zip(pcts, q_pct)))
        buf_p99 = int(np.percentile(qdepths, 99)) * job_B
        buf_max = int(qdepths.max())               * job_B
        print(f"  Buffer (bytes) : p99={fmt_size(buf_p99)}  "
              f"max={fmt_size(buf_max)}  (@ {fmt_size(job_B)}/job)")

    tl   = results["total_latency"]
    tl_p = np.percentile(tl, pcts)
    print(f"\n  {'─'*60}")
    print(f"  End-to-end latency  (trigger arrival → written to disk)")
    print(f"  " + "  ".join(f"p{q:.0f}={v:.1f}s" for q, v in zip(pcts, tl_p))
          + f"  max={tl.max():.1f}s")
    print("=" * W)


def draw_simulation(results: dict, params: dict,
                    out_path: str = "daq_simulation.png"):
    """
    Plot queue depth over time and wait-time distributions for each
    pipeline stage.
    """
    stages    = list(results["stages"].items())
    n_stages  = len(stages)
    arrivals  = results["arrival_times"]

    fig, axes = plt.subplots(n_stages, 2,
                             figsize=(14, 4 * n_stages),
                             facecolor=BG)
    fig.suptitle("DAQ Trigger Pipeline — Discrete-Event Simulation",
                 fontsize=13, color=TXTMAIN, fontweight="bold", y=0.995)

    for ax in axes.flat:
        ax.set_facecolor(PANEL)
        for sp in ax.spines.values():
            sp.set_color(BORDER); sp.set_linewidth(0.8)
        ax.tick_params(colors=TXTSUB)

    for row, (name, r) in enumerate(stages):
        qdepths = r["queue_depths"]
        waits   = r["wait_times"]
        job_B   = r["job_size_B"]

        # ── left: queue depth over (simulated) time ──────────────
        ax = axes[row, 0]
        ax.plot(arrivals, qdepths, color="#4d9de0", lw=0.6, alpha=0.7)
        p99_q = np.percentile(qdepths, 99)
        ax.axhline(p99_q, color="#f85149", lw=1.2, ls="--",
                   label=f"p99 = {p99_q:.0f} jobs")
        ax.set_title(f"{name} — queue depth over time",
                     color=TXTSUB, fontsize=9)
        ax.set_xlabel("Simulated time (s)", color=TXTSUB, fontsize=8)
        ax.set_ylabel("Jobs waiting",       color=TXTSUB, fontsize=8)
        ax.legend(fontsize=8, facecolor=PANEL,
                  labelcolor=TXTMAIN, edgecolor=BORDER)

        # annotate buffer size at p99
        buf_p99 = int(p99_q) * job_B
        ax.text(0.98, 0.97, f"p99 buffer ≈ {fmt_size(buf_p99)}",
                transform=ax.transAxes, ha="right", va="top",
                color="#f85149", fontsize=8, fontweight="bold")

        # ── right: wait time histogram ────────────────────────────
        ax = axes[row, 1]
        nonzero = waits[waits > 0]
        if len(nonzero):
            ax.hist(nonzero, bins=60, color="#4d9de0", alpha=0.8,
                    edgecolor="none")
        p99_w = np.percentile(waits, 99)
        ax.axvline(p99_w, color="#f85149", lw=1.4, ls="--",
                   label=f"p99 = {p99_w:.2f} s")
        ax.set_title(f"{name} — wait time distribution",
                     color=TXTSUB, fontsize=9)
        ax.set_xlabel("Wait time (s)",       color=TXTSUB, fontsize=8)
        ax.set_ylabel("Trigger count",       color=TXTSUB, fontsize=8)
        ax.legend(fontsize=8, facecolor=PANEL,
                  labelcolor=TXTMAIN, edgecolor=BORDER)

    plt.tight_layout()
    plt.savefig(out_path, dpi=150, bbox_inches="tight", facecolor=BG)
    print(f"Simulation figure saved → {out_path}")
    plt.close()


# ── Build and run ─────────────────────────────────────────────
if __name__ == "__main__":
    graph = FlowGraph(
        nodes       = NODES,
        edges       = EDGES,
        params      = PARAMS,
        annotations = ANNOTATIONS,
    )
    graph.evaluate()
    graph.print_summary()
    draw_graph(graph, out_path="daq_flow_diagram.png")

    sim = simulate_flow(graph, n_triggers=10_000)
    print_simulation_summary(sim, PARAMS)
    draw_simulation(sim, PARAMS, out_path="daq_simulation.png")