dlsim4gmns

13def fd_state(k, params):
14    """Map density k (scalar or array) to (v, q) on the triangular FD.
15
16    params: mapping/Series/row with v_f, q_cap, k_crit, k_jam.
17      free branch   k <= k_crit : v = v_f,          q = v_f * k
18      congested     k >  k_crit : q = q_cap*(k_jam-k)/(k_jam-k_crit), v = q/k
19    Returns (v, q) as numpy arrays with v <= v_f and q = k*v.
20    """
21    k = np.asarray(k, dtype=float)
22    vf = float(params["v_f"]); qc = float(params["q_cap"])
23    kc = float(params["k_crit"]); kj = float(params["k_jam"])
24    q = np.where(k <= kc, vf * k,
25                 qc * np.maximum(kj - k, 0.0) / max(kj - kc, 1e-9))
26    v = np.where(k > 1e-9, q / np.maximum(k, 1e-9), vf)
27    return np.minimum(v, vf), q

30def density_from_speed(v, params):
31    """Invert the congested branch: speed -> density (for speed->volume, Level 4).
32    Under-determined in free flow (v ~= v_f); returns the congested-branch density."""
33    v = np.asarray(v, dtype=float)
34    qc = float(params["q_cap"])
35    kc = float(params["k_crit"]); kj = float(params["k_jam"])
36    return qc * kj / (v * (kj - kc) + qc)

39def fd_from_attributes(free_speed_kmh, capacity_vph, lanes):
40    """Derive FD params from OSM/GMNS link attributes (geography-agnostic, Level 4)."""
41    vf = float(free_speed_kmh); qc = float(capacity_vph)
42    return dict(v_f=vf, q_cap=qc, k_crit=qc / max(vf, 1.0), k_jam=float(lanes) * 130.0)

11def read_network(net_dir: str) -> dict:
12    """Load a GMNS network directory. Tolerates both the competition schema
13    (free_speed_kmh/capacity_vph) and osm2gmns output (free_speed/capacity)."""
14    node = pd.read_csv(os.path.join(net_dir, "node.csv"))
15    link = pd.read_csv(os.path.join(net_dir, "link.csv"))
16    return {"nodes": node, "links": link, "fd": fd_table(link)}

26def fd_table(links: pd.DataFrame) -> pd.DataFrame:
27    """Per-link triangular FD params from link attributes."""
28    vf = pd.to_numeric(_col(links, "free_speed_kmh", "free_speed"), errors="coerce").fillna(100.0)
29    cap = pd.to_numeric(_col(links, "capacity_vph", "capacity"), errors="coerce")
30    lanes = pd.to_numeric(_col(links, "lanes"), errors="coerce").fillna(2).clip(lower=1)
31    cap = cap.where(cap > 0, lanes * 1900.0)
32    rows = [fd_from_attributes(v, c, n) for v, c, n in zip(vf, cap, lanes)]
33    fd = pd.DataFrame(rows)
34    fd.insert(0, "link_id", links["link_id"].values)
35    return fd

12def simulate(net: dict, tod_step_min: int = 15, peak_load: float = 1.15) -> pd.DataFrame:
13    """Produce a diurnal density-canonical state {speed, flow, density} per (link, tod).
14    Density follows a free-flow floor plus AM/PM Gaussian demand bumps, pushed past
15    k_crit on lower-capacity links (the bottlenecks). q = k*v holds by construction.
16    """
17    fd = net["fd"].set_index("link_id")
18    lanes = net["links"].set_index("link_id")
19    tods = [f"{h:02d}:{m:02d}" for h in range(24) for m in range(0, 60, tod_step_min)]
20    tmin = np.array([int(t[:2]) * 60 + int(t[3:]) for t in tods])
21    demand = np.exp(-((tmin - 450) / 90) ** 2) + np.exp(-((tmin - 1020) / 100) ** 2)
22    rows = []
23    for lid in fd.index:
24        r = fd.loc[lid]
25        ln = float(lanes.loc[lid].get("lanes", 3)) if lid in lanes.index else 3
26        load = peak_load if ln <= 2 else 0.85
27        k = 0.15 * r["k_crit"] + demand * load * r["k_crit"]
28        k = np.minimum(k, r["k_jam"])
29        v, q = fd_state(k, r)
30        for t, vv, qq, kk in zip(tods, v, q, k):
31            rows.append(dict(link_id=lid, tod=t, timestamp=f"2026-06-15T{t}:00Z",
32                             speed=vv, flow=qq, density=float(kk)))
33    return pd.DataFrame(rows)

52def validate(state: pd.DataFrame, fd: pd.DataFrame) -> tuple[bool, dict]:
53    """Run the gate. Returns (gate_passed, report). Gate-first: a failing gate means
54    the accuracy score is zero regardless of RMSE."""
55    checks = [check_schema(state)]
56    if checks[0]["passed"]:
57        checks.append(check_state(state, fd))
58    gate = all(c["passed"] for c in checks)
59    return gate, dict(gate_passed=gate, checks=checks)

20def check_schema(state: pd.DataFrame) -> dict:
21    """G0: required columns present; speed/flow/density finite and non-negative."""
22    missing = [c for c in REQUIRED if c not in state.columns]
23    if missing:
24        return dict(check="G0_schema", passed=False, detail={"missing": missing})
25    bad = 0
26    for c in ("speed", "flow", "density"):
27        v = state[c].astype(float)
28        bad += int((v < -1e-6).sum() + v.isna().sum())
29    return dict(check="G0_schema", passed=bad == 0, detail={"bad_or_missing": bad})

32def check_state(state: pd.DataFrame, fd: pd.DataFrame) -> dict:
33    """G1: q = k*v within tol, and v<=v_f, k<=k_jam, q<=q_cap (FD feasibility)."""
34    m = state.merge(fd, on="link_id", how="left")
35    v = m["speed"].astype(float).values
36    q = m["flow"].astype(float).values
37    k = m["density"].astype(float).values
38    qkv = np.abs(q - k * v) / np.maximum(np.abs(q), 1.0)
39    qkv_ok = _frac_ok(qkv <= QKV_TOL)
40    vf = m.get("v_f", pd.Series(np.full(len(m), 1e9))).astype(float).values
41    kj = m.get("k_jam", pd.Series(np.full(len(m), 1e9))).astype(float).values
42    cap = m.get("q_cap", pd.Series(np.full(len(m), 1e9))).astype(float).values
43    v_ok = _frac_ok(v <= vf * 1.05)
44    k_ok = _frac_ok(k <= kj * 1.05)
45    q_ok = _frac_ok(q <= cap * 1.15)
46    passed = all((1 - x) <= VIOL_FRAC for x in (qkv_ok, v_ok, k_ok, q_ok))
47    return dict(check="G1_state (q=k*v, FD-feasible)", passed=bool(passed),
48                detail=dict(qkv_ok=round(qkv_ok, 4), v_le_vf=round(v_ok, 4),
49                            k_le_kjam=round(k_ok, 4), q_le_cap=round(q_ok, 4)))

17def run(net_dir: str | None = None) -> dict:
18    """Run load -> simulate -> gate on a network directory (defaults to the sample)."""
19    net = read_network(net_dir or sample_network_dir())
20    state = simulate(net)
21    gate, report = validate(state, net["fd"])
22    return dict(net_dir=net_dir or "sample_corridor",
23                n_links=len(net["links"]), n_cells=len(state),
24                gate=gate, report=report, state=state)

12def sample_network_dir() -> str:
13    """Path to the bundled sample corridor shipped inside the package."""
14    return os.path.join(os.path.dirname(__file__), "data", "sample_corridor")

22def plot_fd(params, ax=None):
23    """The triangular fundamental diagram: flow q(k) and speed v(k) vs density."""
24    plt = _mpl()
25    kj = float(params["k_jam"])
26    k = np.linspace(0, kj, 200)
27    v, q = fd_state(k, params)
28    if ax is None:
29        fig, ax = plt.subplots(figsize=(5, 3.4), dpi=120)
30    else:
31        fig = ax.figure
32    ax.plot(k, q, color="#8c1d40", lw=2, label="flow $q(k)$")
33    ax.axvline(float(params["k_crit"]), ls="--", color="#888", lw=1)
34    ax.set_xlabel("density $k$ (veh/km)"); ax.set_ylabel("flow $q$ (veh/h)")
35    ax.set_title("Fundamental diagram ($q=k\\,v$)")
36    ax2 = ax.twinx(); ax2.plot(k, v, color="#0a6", lw=1.2, alpha=.7, label="speed $v(k)$")
37    ax2.set_ylabel("speed $v$ (km/h)")
38    fig.tight_layout()
39    return fig

42def spacetime_heatmap(state: pd.DataFrame, links: pd.DataFrame | None = None,
43                      value: str = "speed", vmax: float | None = None):
44    """Space-time heatmap: link position (y) x time-of-day (x), colored by `value`.
45    Links are ordered by `links` (if given, by row order) else by link_id."""
46    plt = _mpl()
47    order = list(links["link_id"]) if links is not None else sorted(state["link_id"].unique())
48    tods = sorted(state["tod"].unique())
49    pv = state.pivot_table(value, "link_id", "tod", aggfunc="mean").reindex(index=order, columns=tods)
50    fig, ax = plt.subplots(figsize=(7, 3.6), dpi=120)
51    vmx = vmax if vmax is not None else (110 if value == "speed" else float(np.nanmax(pv.values)))
52    im = ax.imshow(pv.values, aspect="auto", cmap=_SPEED_CMAP, vmin=0, vmax=vmx, origin="lower",
53                   extent=[0, 24, 0, len(order)])
54    ax.set_xlabel("hour of day"); ax.set_ylabel("link (upstream $\\to$ downstream)")
55    ax.set_xticks(range(0, 25, 4)); ax.set_title(f"Space-time {value}")
56    fig.colorbar(im, ax=ax, label=value, shrink=.9); fig.tight_layout()
57    return fig

 78def network_map(net: dict, state: pd.DataFrame | None = None, value: str = "speed"):
 79    """Network geometry (nodes/links); links colored by mean `value` if a state is given."""
 80    plt = _mpl()
 81    nd, lk = net["nodes"], net["links"]
 82    pos = {r.node_id: (float(r.x_coord), float(r.y_coord)) for r in nd.itertuples()}
 83    col = None
 84    if state is not None:
 85        col = state.groupby("link_id")[value].mean()
 86    fig, ax = plt.subplots(figsize=(6.5, 4), dpi=120)
 87    import matplotlib
 88    from matplotlib.colors import Normalize
 89    norm = Normalize(0, 110); cmap = matplotlib.colormaps[_SPEED_CMAP]
 90    for r in lk.itertuples():
 91        a, b = pos.get(r.from_node_id), pos.get(r.to_node_id)
 92        if not a or not b:
 93            continue
 94        c = cmap(norm(col.get(r.link_id, np.nan))) if col is not None else "#0366d6"
 95        ax.plot([a[0], b[0]], [a[1], b[1]], "-", color=c, lw=2 + 0.6 * getattr(r, "lanes", 2))
 96    xs = [p[0] for p in pos.values()]; ys = [p[1] for p in pos.values()]
 97    ax.scatter(xs, ys, s=8, c="#333", zorder=3)
 98    ax.set_xlabel("lon"); ax.set_ylabel("lat"); ax.set_aspect("equal")
 99    ax.set_title("Network" + (f" (mean {value})" if col is not None else "")); fig.tight_layout()
100    return fig

60def plot_odme_fit(y_obs, y_est, ax=None):
61    """ODME goodness-of-fit: observed vs estimated link counts (the 45-degree line)."""
62    plt = _mpl()
63    y_obs = np.asarray(y_obs, float); y_est = np.asarray(y_est, float)
64    if ax is None:
65        fig, ax = plt.subplots(figsize=(4.2, 4.2), dpi=120)
66    else:
67        fig = ax.figure
68    hi = max(float(y_obs.max()), float(y_est.max())) * 1.05
69    ax.plot([0, hi], [0, hi], "--", color="#888", lw=1)
70    ax.scatter(y_obs, y_est, s=28, color="#8c1d40", alpha=.8)
71    mape = float(np.mean(np.abs(y_est - y_obs)) / max(y_obs.mean(), 1) * 100)
72    ax.set_xlabel("observed link count"); ax.set_ylabel("ODME estimated")
73    ax.set_title(f"ODME fit (link MAPE {mape:.1f}%)"); ax.set_aspect("equal")
74    fig.tight_layout()
75    return fig

12def gls_odme(B, v_hat, A_OD=None, q_hat=None, lam=0.05, w=None):
13    """Non-negative, prior-regularised least squares.
14
15    B     : (n_sensor_links x P)  sensor-link x path incidence
16    v_hat : (n_sensor_links,)     observed link counts
17    A_OD  : (n_od x P)            OD x path incidence (optional prior term)
18    q_hat : (n_od,)               OD prior volumes (optional)
19    Returns dict(x=path flow, y=fitted link flow, d=OD, link_mape).
20    """
21    from scipy.optimize import lsq_linear
22    B = np.asarray(B, float); v_hat = np.asarray(v_hat, float)
23    W = np.ones(len(v_hat)) if w is None else np.asarray(w, float)
24    blocks = [np.sqrt(W)[:, None] * B]
25    targets = [np.sqrt(W) * v_hat]
26    if A_OD is not None and q_hat is not None and lam > 0:
27        A_OD = np.asarray(A_OD, float)
28        blocks.append(np.sqrt(lam) * A_OD)
29        targets.append(np.sqrt(lam) * np.asarray(q_hat, float))
30    M = np.vstack(blocks); t = np.concatenate(targets)
31    x = lsq_linear(M, t, bounds=(0, np.inf), max_iter=2000).x
32    y = B @ x
33    d = A_OD @ x if A_OD is not None else None
34    mape = float(np.mean(np.abs(y - v_hat)) / max(v_hat.mean(), 1) * 100)
35    return dict(x=x, y=y, d=d, link_mape=mape)

13def accumulate(lam, mu, N0=0.0, dt=1.0):
14    """Integrate net flow to the vehicle-number trajectory N(t).
15    For T net-flow steps returns T+1 values: N[0]=N0, N[i+1]=N[i]+dt*(lam[i]-mu[i])."""
16    lam = np.asarray(lam, float); mu = np.asarray(mu, float)
17    return np.concatenate([[N0], N0 + np.cumsum(dt * (lam - mu))])

20def conservation_residual(N, lam, mu, dt=1.0):
21    """Mean absolute violation of N(t+dt) - N(t) = dt*(lambda - mu). 0.0 == conserved.
22    N has length len(lam)+1 (the trajectory); lam, mu are the per-step net flows."""
23    N = np.asarray(N, float); lam = np.asarray(lam, float); mu = np.asarray(mu, float)
24    dN = np.diff(N)
25    expected = dt * (lam - mu)
26    return float(np.mean(np.abs(dN - expected))) if len(dN) else 0.0

10def volume_from_speed(v, params, cap_drop=0.70):
11    """Regime-aware volume estimate + 1-sigma uncertainty from a speed field.
12
13      congested (v < 0.75 v_f): near capacity, q ~= cap_drop * C   (tight sigma)
14      transition (0.75-0.9 v_f): FD congested branch                (loose sigma)
15      free (v >= 0.9 v_f):       UNDER-DETERMINED -> (nan, inf)     (defer to a prior)
16
17    Returns (q_hat, sigma) as arrays; free-flow cells carry (nan, inf) so an ODME step
18    weights them out (weight = 1/sigma^2 = 0) and lets the OD/diurnal prior fill them.
19    """
20    v = np.asarray(v, float)
21    vf = float(params["v_f"]); qc = float(params["q_cap"])
22    kc = float(params["k_crit"]); kj = float(params["k_jam"])
23    r = v / max(vf, 1e-6)
24    q = np.full_like(v, np.nan); sig = np.full_like(v, np.inf)
25    cong = r < 0.75
26    q = np.where(cong, cap_drop * qc, q)
27    sig = np.where(cong, 0.30 * cap_drop * qc, sig)
28    tr = (r >= 0.75) & (r < 0.9)
29    kt = qc * kj / (v * (kj - kc) + qc)          # congested-branch density
30    q = np.where(tr, kt * v, q)
31    sig = np.where(tr, 0.60 * np.maximum(kt * v, 1.0), sig)
32    return q, sig                                 # free flow stays (nan, inf)

21def newell_corridor(length, lanes, vf, w, kjam, qmax, inflow, dt_sec: float = 6.0):
22    length = np.ascontiguousarray(length, float)
23    lanes = np.ascontiguousarray(lanes, float)
24    vf = np.ascontiguousarray(vf, float)
25    w = np.ascontiguousarray(w, float)
26    kjam = np.ascontiguousarray(kjam, float)
27    inflow = np.ascontiguousarray(inflow, float)
28    L, T = len(length), len(inflow)
29    qmax = np.ascontiguousarray(np.asarray(qmax, float).reshape(L, T))
30    dT_h = dt_sec / 3600.0
31
32    # per-link free-flow / backward-wave lags (ticks) and storage (veh) -- exactly the C++
33    fftt = np.maximum(1, (length / vf / dT_h + 0.5).astype(int))
34    bwtt = np.maximum(1, (length / w / dT_h + 0.5).astype(int))
35    storage = kjam * length * lanes
36
37    A = np.zeros((L, T)); D = np.zeros((L, T))
38    A[0, 1:] = np.cumsum(inflow[1:] * dT_h)              # cumulative demand into link 0
39
40    for t in range(1, T):
41        Vt = np.empty(L); capin = np.empty(L)
42        for l in range(L):                               # Step 2: states at t
43            Vt[l] = A[l, t - fftt[l]] if t - fftt[l] >= 0 else 0.0        # forward wave
44            D_bw = D[l, t - bwtt[l]] if t - bwtt[l] >= 0 else 0.0
45            capin[l] = max(0.0, D_bw + storage[l] - A[l, t - 1])         # backward wave + storage
46        for l in range(L):                               # Step 3: capacity-constrained transfer
47            capout = min(qmax[l, t] * dT_h, capin[l + 1]) if l < L - 1 else qmax[l, t] * dT_h
48            Dt = min(Vt[l], D[l, t - 1] + capout)
49            D[l, t] = Dt
50            if l < L - 1:
51                A[l + 1, t] = Dt                         # A(l+1,t) = D(l,t)
52    return A, D

 77def newell_merge(length, lanes, vf, w, kjam, qmax, downstream, inflow, dt_sec: float = 6.0,
 78                 off_down=None, off_frac=None):
 79    """Newell KW loading on a freeway-junction topology (merge + diverge).
 80
 81    Movement-based: each link sends a mainline movement to `downstream[l]` and, if
 82    `off_frac[l] > 0`, an off-ramp movement (fraction `off_frac[l]` of its outflow) to
 83    `off_down[l]`. A downstream link's inflow capacity is split across the movements entering
 84    it **proportional to the # of incoming lanes** (DTALite merge). Movements are
 85    independent, so an off-ramp can never block the mainline (off-ramp spillback disabled --
 86    see docs/assumptions.md). `downstream`/`off_down` = -1 means the movement exits the
 87    network. Generalizes `newell_corridor` and the pure merge.
 88
 89    inflow : (L, T) external demand entering each origin link (0 for interior links); unmet
 90             demand waits in a per-origin queue off-network. Returns (A, D), each (L, T).
 91    """
 92    length = np.asarray(length, float)
 93    lanes = np.maximum(np.asarray(lanes, float), 1e-6)   # guard the lane-proportional split
 94    vf = np.asarray(vf, float); w = np.asarray(w, float); kjam = np.asarray(kjam, float)
 95    downstream = np.asarray(downstream, int)
 96    L = len(length)
 97    off_down = np.full(L, -1, int) if off_down is None else np.asarray(off_down, int)
 98    off_frac = np.zeros(L) if off_frac is None else np.asarray(off_frac, float)
 99    inflow = np.ascontiguousarray(np.asarray(inflow, float).reshape(L, -1))
100    T = inflow.shape[1]
101    qmax = np.ascontiguousarray(np.asarray(qmax, float).reshape(L, T))
102    dT_h = dt_sec / 3600.0
103
104    fftt = np.maximum(1, (length / vf / dT_h + 0.5).astype(int))
105    bwtt = np.maximum(1, (length / w / dT_h + 0.5).astype(int))
106    storage = kjam * length * lanes
107    # incoming movements per downstream link: (from_link, is_offramp)
108    in_moves = [[] for _ in range(L)]
109    for l in range(L):
110        if downstream[l] >= 0:
111            in_moves[downstream[l]].append((l, False))
112        if off_down[l] >= 0 and off_frac[l] > 0:
113            in_moves[off_down[l]].append((l, True))
114    is_origin = [not in_moves[l] for l in range(L)]
115
116    A = np.zeros((L, T)); D = np.zeros((L, T)); oq = np.zeros(L)
117    for l in range(L):
118        if is_origin[l]:
119            oq[l] += inflow[l, 0] * dT_h                # inject t=0 demand (loop starts at t=1)
120    for t in range(1, T):
121        S = np.zeros(L); capin = np.zeros(L)
122        for l in range(L):
123            Vl = A[l, t - fftt[l]] if t - fftt[l] >= 0 else 0.0
124            S[l] = min(max(0.0, Vl - D[l, t - 1]), qmax[l, t] * dT_h)
125            Dbw = D[l, t - bwtt[l]] if t - bwtt[l] >= 0 else 0.0
126            capin[l] = max(0.0, Dbw + storage[l] - A[l, t - 1])
127        main_dem = S * (1.0 - off_frac); off_dem = S * off_frac
128        alloc = np.zeros(L); arrivals = np.zeros(L)
129        for l in range(L):                              # movements that exit the network
130            if downstream[l] < 0:
131                alloc[l] += main_dem[l]
132            if off_down[l] < 0 and off_frac[l] > 0:
133                alloc[l] += off_dem[l]
134        for d in range(L):                              # merge over movements entering d
135            moves = in_moves[d]
136            if not moves:
137                continue
138            dem = [(off_dem[fl] if is_off else main_dem[fl]) for (fl, is_off) in moves]
139            wt = [lanes[fl] for (fl, _) in moves]
140            acc = _cap_split(dem, wt, range(len(moves)), capin[d])
141            for k, (fl, _) in enumerate(moves):
142                alloc[fl] += acc[k]; arrivals[d] += acc[k]
143        for l in range(L):                              # admit external demand at origins
144            if is_origin[l]:
145                oq[l] += inflow[l, t] * dT_h
146                a = min(oq[l], capin[l]); arrivals[l] += a; oq[l] -= a
147        for l in range(L):
148            D[l, t] = D[l, t - 1] + alloc[l]
149            A[l, t] = A[l, t - 1] + arrivals[l]
150    return A, D

153def corridor_engine(prefer_native: bool = True):
154    """Return (fn, backend_name): the native corridor loader if built, else pure-Python.
155    Both take the same signature and return identical (A, D)."""
156    if prefer_native:
157        try:
158            from . import native
159            if native.available():
160                return native.newell_corridor, "native-c++"
161            native.note_unavailable()
162        except Exception:
163            pass
164    return newell_corridor, "python"

30def load_scenario(scn_dir: str) -> dict:
31    """Read a scenario directory (node/link/demand csv + scenario.json)."""
32    net = read_network(scn_dir)
33    demand = pd.read_csv(os.path.join(scn_dir, "demand.csv"))
34    with open(os.path.join(scn_dir, "scenario.json")) as f:
35        cfg = json.load(f)
36    return {"dir": scn_dir, "net": net, "demand": demand, "cfg": cfg}

142def run_scenario(scn: dict, demand_scale: float = 1.0, bin_min: int = 5,
143                 prefer_native: bool = True) -> dict:
144    """Run merge dynamic loading and derive the density-canonical state, queue profile,
145    bottleneck report, and physics gate. Returns a rich result dict."""
146    p = _build_inputs(scn, demand_scale); topo = p["topo"]
147    engine, backend = merge_engine(prefer_native)
148    A, D = engine(p["length"], p["lanes"], p["vf"], p["w"], p["kjam_per_lane"],
149                  p["qmax"], np.asarray(topo["downstream"]), p["inflow"], float(p["tick"]),
150                  off_down=p["off_down"], off_frac=p["off_frac"])
151    A = np.asarray(A); D = np.asarray(D)
152    L, T = A.shape
153    on = A - D
154    dt_h = p["tick"] / 3600.0
155    bin_ticks = max(1, int(bin_min * 60 // p["tick"]))
156    fd = scn["net"]["fd"].set_index("link_id")
157
158    state_rows, queue_rows = [], []
159    for i, lid in enumerate(topo["lids"]):
160        length_i, vf_i = p["length"][i], p["vf"][i]
161        for b0 in range(0, T - 1, bin_ticks):
162            b1 = min(b0 + bin_ticks, T - 1)
163            hrs = (b1 - b0) * dt_h
164            if hrs <= 0:
165                continue
166            flow = (D[i, b1] - D[i, b0]) / hrs            # veh/h discharged in the bin
167            n_avg = float(on[i, b0:b1].mean())            # mean vehicles on link
168            dens = n_avg / max(length_i, 1e-6)            # total veh/mi
169            spd = flow / dens if dens > 1e-6 else vf_i
170            if spd > vf_i:                                # free-flow discretization overshoot:
171                spd, flow = vf_i, dens * vf_i             # cap v at v_f, keep q = k*v exact
172            tod = "%02d:%02d" % (int(p["base"] + b0 * dt_h * 60) // 60,
173                                 int(p["base"] + b0 * dt_h * 60) % 60)
174            tt_min = length_i / max(spd, 1e-6) * 60.0
175            state_rows.append(dict(link_id=lid, tod=tod, speed=spd, flow=flow,
176                                   density=dens, queue=max(0.0, n_avg), travel_time_min=tt_min))
177        # congestion episode = the contiguous over-capacity run containing the peak queue
178        thresh = float(fd.loc[lid]["k_crit"]) * length_i    # ~free-flow content
179        above = on[i] > thresh
180        onset = clear = None
181        if above.any():
182            peak = int(on[i].argmax())
183            s = peak
184            while s > 0 and above[s - 1]:
185                s -= 1
186            e = peak
187            while e < T - 1 and above[e + 1]:
188                e += 1
189            onset, clear = s, (e + 1 if e < T - 1 else None)
190        to_min = lambda t: None if t is None else round(p["base"] + t * dt_h * 60, 1)
191        queue_rows.append(dict(link_id=lid, max_on_link=float(on[i].max()),
192                               onset_min=to_min(onset), clear_min=to_min(clear),
193                               duration_min=None if (onset is None or clear is None) else round((clear - onset) * dt_h * 60, 1)))
194
195    state = pd.DataFrame(state_rows)
196    gate, report = validate(state[["link_id", "speed", "flow", "density"]], scn["net"]["fd"])
197    bottleneck = (pd.DataFrame(queue_rows).sort_values("max_on_link", ascending=False)
198                  .reset_index(drop=True))
199    released = float(sum(A[topo["idx"][lid], -1] for lid in topo["origins"]))
200    on_end = float(on[:, -1].sum())                     # vehicles still on the network
201    throughput = released - on_end                      # left via a mainline exit or an off-ramp
202    conserved = on_end < max(5.0, 0.001 * max(released, 1.0))   # network cleared (no residual)
203    return dict(backend=backend, A=A, D=D, on=on, topo=topo, tmin=p["tmin"],
204                inflow=p["inflow"], qmax=p["qmax"], state=state, gate=gate,
205                gate_report=report, bottleneck=bottleneck, conserved=bool(conserved),
206                throughput=throughput, released=released, notes=topo["notes"])

209def recover_demand_scale(scn: dict, observed_throughput: float,
210                         lo: float = 0.5, hi: float = 1.6, iters: int = 24) -> dict:
211    """Illustrative 1-parameter ODME: recover the demand multiplier whose simulated
212    corridor throughput matches an observed total, by bisection on the monotone map
213    scale -> throughput. The full multi-target ODME is `dlsim4gmns.odme` (FTT/GLS)."""
214    def err(s):
215        return run_scenario(scn, demand_scale=s)["throughput"] - observed_throughput
216    e_lo, e_hi = err(lo), err(hi)
217    if e_lo * e_hi > 0:
218        best = lo if abs(e_lo) < abs(e_hi) else hi
219        return dict(scale=best, matched=False,
220                    throughput=observed_throughput + (e_lo if best == lo else e_hi))
221    for _ in range(iters):
222        mid = 0.5 * (lo + hi)
223        if err(mid) * e_lo > 0:
224            lo = mid
225        else:
226            hi = mid
227    scale = 0.5 * (lo + hi)
228    return dict(scale=scale, matched=True,
229                throughput=run_scenario(scn, demand_scale=scale)["throughput"])

39def topology(net: dict) -> dict:
40    """Infer junction topology from GMNS: downstream[l] = the MAINLINE link l feeds into (the
41    outgoing link with the most lanes at its end node, or -1 at an exit), the origin links (no
42    incoming link), and each origin's zone. Off-ramp branches at a diverge node are activated
43    by declaring them in scenario.json 'diverges' (see docs/assumptions.md)."""
44    links, nodes = net["links"], net["nodes"]
45    lids = list(links["link_id"])
46    idx = {lid: i for i, lid in enumerate(lids)}
47    frm = dict(zip(links["link_id"], links["from_node_id"]))
48    to = dict(zip(links["link_id"], links["to_node_id"]))
49    lanes_of = dict(zip(links["link_id"], pd.to_numeric(
50        links.get("lanes", pd.Series([2] * len(links))), errors="coerce").fillna(1)))
51    out_by_node, in_nodes, notes = {}, set(), []
52    for lid in lids:
53        out_by_node.setdefault(frm[lid], []).append(lid)
54        in_nodes.add(to[lid])
55    downstream = []
56    for lid in lids:
57        outs = out_by_node.get(to[lid], [])
58        if not outs:
59            downstream.append(-1)
60            continue
61        main = max(outs, key=lambda x: lanes_of.get(x, 1))    # mainline = most lanes
62        downstream.append(idx[main])
63        if len(outs) > 1:
64            offs = [o for o in outs if o != main]
65            notes.append(f"node {to[lid]}: diverge (mainline={main}, off-ramp(s)={offs}); "
66                         f"declare in scenario.json 'diverges' to route off-ramp flow")
67    origins = [lid for lid in lids if frm[lid] not in in_nodes]
68    zone = dict(zip(nodes["node_id"], nodes.get("zone_id", pd.Series([None] * len(nodes)))))
69    origin_zone = {lid: zone.get(frm[lid]) for lid in origins}
70    return dict(lids=lids, idx=idx, downstream=downstream, origins=origins,
71                origin_zone=origin_zone, exits=[i for i, d in enumerate(downstream) if d < 0],
72                notes=notes)