# ADR-001 — Logit SUE: Dial's STOCH loading, MSA-SUE, and the fixed-point certificate

**Status:** accepted (implemented in v0.x)
**File:** `docs/design/adr-001-logit-sue-dial-certificate.md`

## Context

The roadmap (ARCHITECTURE §5, v0.x) calls for logit stochastic user equilibrium as the
first route-choice component, and T1 already promises "SUE fixed-point residual for SUE
tasks" (ARCHITECTURE §2). The hard question is not the solver — MSA with a stochastic
loader is textbook (Sheffi 1985, ch. 12) — but the **certificate**: P1 demands that the
harness recompute every scored metric from `(link_flows, scenario)` alone, for black
boxes included. For UE this is the externally computed relative gap. For SUE there is no
analogue via an objective: Fisk's (1980) equivalent convex program is

```
min  Σ_a ∫₀^{v_a} t_a(s) ds  +  (1/θ) Σ_rs Σ_p f_p^rs ln f_p^rs
```

whose entropy term is a function of **path** flows `f_p`, and a link-flow vector does not
determine its path-flow decomposition (nor, under Dial, even the admissible path set,
which moves with costs). So the Fisk objective cannot be evaluated from emitted link
flows and cannot be the certificate. What *can* be evaluated from link flows is the
defining fixed-point condition of SUE itself: `v = L(t(v), θ)`, where `L` is a pinned,
deterministic stochastic-loading map. This ADR pins that map (Dial's STOCH), the solver
(MSA), the certificate (an L1 fixed-point residual in AEC-like per-traveler units), and
where `θ` lives (an optional `Scenario` field that leaves every existing content hash
unchanged).

Terminology note: "stochastic" in SUE refers to traveler perception, not algorithmic
randomness. Dial's logit loading is a closed-form deterministic function of `(t, θ)` —
no RNG, no Monte Carlo — so this model runs on the *deterministic* track (M=1). P5's
"stochastic track: SUE" wording should be refined to "simulation-based SUE (e.g. probit,
v1)"; logit-SUE-via-Dial belongs with FW and MSA.

## Decision 1 — The loading map: Dial's STOCH, double-pass, origin-based efficiency

`L(t, θ)` is Dial's (1971) STOCH algorithm in the **double-pass form with the
origin-based efficient-link criterion `r(i) < r(j)`**, exactly as presented in Sheffi
(1985, §11.2 — the canonical free reference implementers actually use).

Dial's paper defines a "reasonable" (efficient) path with two labels — every link must
lead *away from the origin* (`r(i) < r(j)`, `r` = shortest time from the origin) *and
toward the destination* (`s(i) > s(j)`, `s` = shortest time to the destination) — but
his recommended double-pass algorithm, and Sheffi's textbook presentation of it, retain
only the origin criterion. We follow that choice because:

1. **Reproducibility.** Sheffi §11.2 is the presentation of record; matching it means a
   third party implementing "Dial per Sheffi" reproduces our certificate bit-for-bit in
   exact arithmetic. The two-label variant makes the efficient set destination-dependent,
   and implementations disagree on details.
2. **Complexity.** Origin-based efficiency lets one forward+backward sweep per origin
   serve *all* destinations (like AON): `O(|origins| · (m + n log n))`. The two-label
   criterion needs per-OD-pair sweeps: `O(|zones|²)` label trees.
3. **Task-definition honesty.** Whatever criterion we pick *defines* the benchmark's SUE
   task; models and harness share the same map, so the choice must be pinned, not
   optimal. We pin the canonical one and record it in the scenario card
   (`sue.loading: dial-stoch`).

### Exact algorithm (per origin `o` with a positive demand row)

All computation is on the **PathEngine expanded graph** (`models/_paths.py`): centroids
below `first_thru_node` keep their index in the tail role and get an arc-less shadow
index in the head role, so no efficient path can traverse a centroid — the distance
labels `r(·)` automatically honor centroid restrictions. Costs `t = Network.link_cost(v)`
are generalized (BPR time + toll/distance fixed cost), strictly positive and finite
(guarded exactly as in `PathEngine._graph`).

1. **Labels.** Dijkstra from `o`'s tail index: `r(i)` for every expanded node.
2. **Efficient links.** `E_o = { a=(i→j) : r(i), r(j) finite and r(i) < r(j) }`, strict
   float comparison, no tolerance (ties excluded; see Decision 3).
3. **Link likelihood** (Sheffi eq. 11.7-style):
   `L_a = exp(x_a)` with `x_a = θ·(r(j) − r(i) − t_a)` for `a ∈ E_o`, else `L_a = 0`.
   By shortest-path optimality `r(j) ≤ r(i) + t_a`, so `x_a ≤ 0` and `L_a ≤ 1`:
   likelihoods can never overflow, and `x_a = 0` exactly on shortest-path links.
4. **Forward pass (node weights, log domain), ascending `r`.** Sheffi's recursion is
   `W(o)=1`, `W(j) = Σ_{(i→j)∈E_o} L_ij · W(i)`. `W(j)` sums `exp(θ(r(j) − C_P))` over
   efficient paths `P: o→j` (telescoping product), and the number of efficient paths can
   be exponential — `W` can overflow float64 on large networks at small θ. We therefore
   compute `b(j) = log W(j)` directly: `b(o)=0`,
   `b(j) = logsumexp_{(i→j)∈E_o} ( x_ij + b(i) )` (max-shifted). This is algebraically
   identical to Sheffi's recursion and overflow-free. Invariant: every reachable node's
   shortest path consists of efficient links with `x=0`, so `W(j) ≥ 1`, `b(j) ≥ 0` —
   no zero-weight node, no 0/0 in the backward pass, ever (guaranteed by the Network
   validation `free_flow_time > 0`).
5. **Backward pass (assignment), descending `r`.** Initialize node volume
   `V(head_index(d)) = q_od` for each destination `d ≠ o` with `q_od > 0` (intrazonal
   demand skipped; unreachable positive-demand destination raises `RuntimeError`,
   mirroring `PathEngine.all_or_nothing`). Sweep nodes in descending `r`; at node `j`
   with `V(j) > 0`, split over incoming efficient links with fractions
   `φ_ij = exp(x_ij + b(i) − b(j))` — each exponent ≤ 0 by definition of `b(j)`, so
   `φ ∈ [0,1]` with `Σφ = 1` in exact arithmetic; renormalize by the float sum so link
   flows conserve demand to machine precision — then `flow_ij += V(j)·φ_ij`,
   `V(i) += V(j)·φ_ij`.

The resulting path shares are exactly the logit model over Dial-efficient paths:
`Pr(P) = exp(−θ C_P) / Σ_{P'∈E_o} exp(−θ C_{P'})`.

**Processing order and ties.** Nodes with equal `r` cannot be joined by an efficient
link (strict inequality), so any order among equal labels is correct; a stable argsort
on `r` (as in PathEngine) keeps runs bit-deterministic. Underflow of `exp(x)` for deep
detours silently rounds negligible path shares to zero — benign, since shortest-path
terms are `exp(0)=1` and never underflow.

**θ units.** θ multiplies generalized-cost differences, so its unit is 1/(native cost
unit) and is **network-specific** (P9): Sioux Falls times are 0.01-hour units, so
`θ = 1.0` there means 100 h⁻¹. Scenario cards must state the unit next to the value.

**θ → ∞ limit.** Shares concentrate on minimum-cost efficient paths, split *evenly*
across exact ties (unlike AON's arbitrary single-tree tie-breaking) — documented, and
irrelevant off ties.

## Decision 2 — The reference solver: MSA-SUE (`sue-msa`)

`models/sue_logit.py::DialSUEModel`, registered as `sue-msa`, mirrors `msa.py`:

- `v_1 = L(t(0), θ)` (stochastic loading at free-flow costs); then for `k = 1, 2, …`:
  `y_k = L(t(v_k), θ)`; checkpoint `v_k` with self-report
  `sue_fixed_point_residual = ‖y_k − v_k‖₁ / D`; step `v_{k+1} = v_k + (y_k − v_k)/k`
  (the repo's MSA indexing; identical to `v_{k+1} = v_k + (1/(k+1))(L(t(v_k)) − v_k)`
  up to an index shift). Powell & Sheffi (1982) prove convergence for predetermined
  Blum step sizes (`Σα_k = ∞`, `Σα_k² < ∞`) given a continuous loading map; Dial's
  cost-dependent efficient sets make `L` only piecewise continuous, so the guarantee is
  local (efficient set constant near the fixed point) — a documented caveat, not a
  blocker (Decision 3).
- θ is read from `scenario.sue_theta`; `None` raises `ValueError` — θ is task data, not
  a model factor, so no solver may tune it (P7).
- Budget convention: one full Dial load = 1 `sp_call` (the same one-batched-sweep unit
  as one `all_or_nothing` in FW/MSA; the Dijkstra sweep dominates both).
- `Capabilities(paradigm="sue", deterministic=True, provides_gap=True, seedable=True)` —
  `"sue"` is already in `PARADIGMS`; `deterministic=True` routes it to M=1 in the runner.
- The self-reported residual at checkpoint `v_k` is byte-comparable to the harness value
  (both are `‖L(t(v_k),θ) − v_k‖₁/D`), giving the FW-style honesty regression for free.

**Why the Fisk objective is not reported:** it needs path flows (entropy term) that link
flows do not determine; under Dial the path set additionally moves with costs. Stated in
the module docstring so nobody "fixes" this later.

## Decision 3 — The certificate: SUE fixed-point residual (extends P1)

For a scenario with `sue_theta = θ`, the harness computes, from emitted flows `v` that
pass the existing demand-aware feasibility audit:

```
t = network.link_cost(v)                    # same costs as the UE certification
v̂ = L(t, θ)                                 # ONE Dial load, deterministic
sue_fixed_point_residual = ‖v − v̂‖₁ / D     # D = total demand
```

**Norm and normalization.** Both `v` (audited) and `v̂` (constructed) route the same OD
totals, so `v − v̂` is a pure redistribution of flow among routes; the L1 norm measures
the total redistributed volume in vehicle·link units. Dividing by total demand `D`
yields *misallocated link-traversals per traveler* — the same intensive, per-traveler
convention as AEC (= excess time per traveler), zero exactly at equilibrium, and
comparable across ladder rungs (Braess `D=6` vs Chicago `D≈10⁶`). L2 has no flow
interpretation and is dominated by high-volume links; L∞ ignores diffuse misallocation;
normalizing by `‖v‖₁` instead of `D` would conflate residual with average trip length.

**Discontinuity caveat.** `t ↦ L(t, θ)` is discontinuous exactly where an efficient set
changes — a label tie `r(i) = r(j)`, a measure-zero event in cost space but reachable at
symmetric networks. If a scenario's fixed point sits *on* a tie, the residual need not
approach 0 and MSA can oscillate between the two branch maps. Within a fixed efficient
set, `L` is smooth and the residual is locally Lipschitz
(`≈ ‖(I − ∇L·∇t)(v − v*)‖₁/D` near `v*`). Consequences: the strict-inequality rule is
pinned (no tolerance parameter to disagree about); SUE scenario cards must verify —
analytically or empirically — that the fixed point is tie-free, recording exceptions in
the known-defects style (P9).

**Expected magnitudes.** At an exact analytic fixed point the residual is pure float
noise (cost evaluation + Dijkstra + logsumexp + brentq root tolerance): observed
`< 1e-12` on the two-route scenario; tests assert `< 1e-8`. For MSA-after-k the harness
residual at checkpoint `v_k` *is* the MSA direction norm; measured on the two-route
network: `2.8e-3` (k=10), `6.2e-6` (k=100), `1.1e-7` (k=500); on real networks expect
markedly slower, roughly `O(1/k)`-like decay (Powell & Sheffi prove no rate).

**Why this is P1-valid for black boxes.** The residual is a function of
`(link_flows, scenario, θ)` only — no path flows, no iterates, no model internals — and
costs one deterministic Dial load per checkpoint (≈ one AON; certification cost for SUE
scenarios roughly doubles). Censoring precedes certification: flows failing the demand
audit get `feasible = 0` and a NaN residual, so an all-zero emission is censored, never
scored. UE metrics (relative gap, AEC, Beckmann) remain **descriptive columns** on SUE
tasks — they are strictly positive at SUE by design — and the SUE task ranks by the
residual alone. Fit-vs-gap stays a pair, never a scalar (ARCHITECTURE §2).

## Decision 4 — Where θ lives: `Scenario.sue_theta`, hash-append-when-set

```python
@dataclass(frozen=True)
class Scenario:
    ...
    reference: ReferenceSolution | None = None
    family: str = field(default="")
    sue_theta: float | None = None      # NEW: logit dispersion, 1/(native cost unit)
```

`__post_init__` gains: `sue_theta` must be `None` or finite and `> 0`.
`content_hash()` gains, **after** the existing array loop and before `hexdigest()`:

```python
if self.sue_theta is not None:
    h.update(f"sue_theta={float(self.sue_theta)!r};".encode())
```

Because the update is conditional and appended last, the byte stream fed to SHA-256 when
`sue_theta is None` is identical to today's — every existing hash (Braess
`cf00f411cdccec88…`, Sioux Falls pins in manifests, `trained_on` lineage tokens) is
preserved, enforced by a new golden-hash regression test. When θ is set the hash
changes, which P2 *requires*: θ parametrizes a scored metric, so two scenarios differing
only in θ are different benchmark instances and must never share a hash. `!r` on the
float matches the existing `tw=…!r;dw=…!r` convention (exact round-trip repr).

**Alternatives rejected.** *Task object:* the right home eventually, but no `Task` class
exists in code (tasks are YAML strings), and θ must be content-hashed *today*; the
Scenario field is forward-compatible — a v1 Task layer can construct/point to the hashed
scenario. *Evaluator argument:* violates P2 — two runs with different θ would share a
`scenario_hash`, a silently different problem masquerading as the same instance, and
manifests would not pin θ. *Model factor:* worst — models tuning the task definition
violates P7 outright.

## Decision 5 — Analytic anchor: the two-route scenario

`Network` validation forbids parallel links, so the classic two-parallel-links example
is realized as two disjoint 2-link routes. `data/builtin.py::two_route_scenario(demand=4.0,
sue_theta=0.5)`: nodes 1 (origin zone), 2 (destination zone), 3, 4; `first_thru_node=1`;
links, in order, with linear costs via the existing `bpr_linear` device:

| link | cost | role |
|---|---|---|
| 1→3 | `t = 1` (constant) | route A first leg |
| 3→2 | `t = 1 + f` | route A, `c_A(f_A) = 2 + f_A` |
| 1→4 | `t = 1` (constant) | route B first leg |
| 4→2 | `t = 0.5 + 2f` | route B, `c_B(f_B) = 1.5 + 2 f_B` |

Since `r(3) = r(4) = 1 < min(c_A, c_B) ∈ [1.5, ∞)` at every nonnegative flow, both
routes are Dial-efficient at *all* costs encountered, so `L` reduces exactly to a binary
logit and the SUE fixed point to the scalar equation
`f_A = D / (1 + exp(θ·(c_A(f_A) − c_B(D − f_A))))`, solvable by `brentq` inside the test
(no hard-coded truth digits). Key verified numbers (θ=0.5, D=4): free-flow Dial load
`(1.7512939965, …)`; fixed point `f_A* = 2.2990959494`; UE `f_A = 2.5`; θ-ladder
`|f_A*(θ) − 2.5| = 0.2009, 0.0724, 0.0318, 0.0034` for θ = 0.5, 2, 5, 50; certificate
`≈ 1.165` after shifting 0.5 units A→B. Full test list in the test-plan section below.

## Decision 6 — Integration

The loader lives in **`src/tabench/models/_stoch.py`** as a `StochEngine` class
(mirroring `_paths.py::PathEngine`), because both the SUE model and the Evaluator need
it and `metrics/gaps.py` *already* imports `..models._paths.PathEngine` — the precedent
direction is `metrics → models → core`, acyclic. `_stoch.py` imports only
`core.scenario` (+ scipy) and reuses PathEngine's expansion. A top-level shared
`routing/` package (moving `_paths.py` too) is the cleaner long-term layering but a
larger refactor; deferred, noted for v1.

Touch points: `Evaluator` gains the residual column when `scenario.sue_theta` is set;
runner CSV gains `sue_fixed_point_residual` and `self_sue_residual`;
`models/sue_logit.py` + registry export; `scenarios/0tworoute-sue.yaml` card with an
explicit `sue: {theta, theta_units, loading}` block; ARCHITECTURE metric-definitions
bullet + P5 wording refinement; README model/metric rows. Exact diffs in the
implementation spec accompanying this ADR.

## Consequences

- Anything emitting feasible link flows on a SUE scenario — a 1971 STOCH heir or a 2025
  GNN — receives an externally certified SUE residual, in per-traveler units, from one
  deterministic Dial load. P1 extends to SUE with no new trust assumptions.
- The efficient-set discontinuity is accepted (canonical fidelity over smoothness) and
  managed by scenario-card vetting rather than by a hidden tie tolerance.
- Existing content hashes, manifests, and lineage tokens are untouched; a golden-hash
  test makes that a regression guarantee.
- Evaluator cost on SUE scenarios roughly doubles per checkpoint (AON + Dial load);
  acceptable, with a shared-Dijkstra optimization available later.

## References

- Dial, R.B. (1971). A probabilistic multipath traffic assignment model which obviates
  path enumeration. *Transportation Research* 5(2), 83–111.
- Sheffi, Y. (1985). *Urban Transportation Networks*, ch. 11 (STOCH), ch. 12 (MSA-SUE).
  Free PDF at sheffi.mit.edu.
- Fisk, C. (1980). Some developments in equilibrium traffic assignment.
  *Transportation Research Part B* 14(3), 243–255.
- Powell, W.B. & Sheffi, Y. (1982). The convergence of equilibrium algorithms with
  predetermined step sizes. *Transportation Science* 16(1), 45–55.
- Daganzo, C.F. & Sheffi, Y. (1977). On stochastic models of traffic assignment.
  *Transportation Science* 11(3), 253–274.