"""Scenario objects: frozen, content-hashed, declarative problem instances.
Design principle P2 (docs/ARCHITECTURE.md): a scenario is data, never code.
Node ids follow the TNTP convention: 1-based, zones are nodes ``1..n_zones``,
and nodes with id below ``first_thru_node`` are zone centroids that cannot
carry through traffic.
"""
from __future__ import annotations
import hashlib
from dataclasses import dataclass, field
import numpy as np
__all__ = [
"Network",
"Demand",
"ElasticDemand",
"CombinedDemand",
"MulticlassDemand",
"ReferenceSolution",
"Scenario",
]
def _as_f64(x: np.ndarray) -> np.ndarray:
return np.ascontiguousarray(x, dtype=np.float64)
[docs]
@dataclass(frozen=True)
class Network:
"""Directed road network with BPR-family link performance functions.
The generalized link cost (in time units) is
t_a(v) = fft_a * (1 + b_a * (v / cap_a)^power_a)
+ toll_weight * toll_a + distance_weight * length_a
where the toll and distance weights are per-network metadata (they are
NOT stored in ``.tntp`` files; see docs/ARCHITECTURE.md P9).
``units`` records per-network unit conventions (e.g. Sioux Falls free-flow
times are 0.01 hours) as documentation; costs are always computed in the
network's native units.
"""
name: str
n_nodes: int
n_zones: int
first_thru_node: int
init_node: np.ndarray # (n_links,) int64, 1-based tail node ids
term_node: np.ndarray # (n_links,) int64, 1-based head node ids
capacity: np.ndarray
length: np.ndarray
free_flow_time: np.ndarray
b: np.ndarray # BPR "alpha"
power: np.ndarray # BPR "beta"
toll: np.ndarray
link_type: np.ndarray
toll_weight: float = 0.0
distance_weight: float = 0.0
units: tuple[tuple[str, str], ...] = ()
def __post_init__(self) -> None:
n = self.n_links
for name in ("term_node", "capacity", "length", "free_flow_time", "b", "power", "toll"):
if len(getattr(self, name)) != n:
raise ValueError(f"Network '{self.name}': column {name!r} has wrong length")
pairs = set(zip(self.init_node.tolist(), self.term_node.tolist(), strict=True))
if len(pairs) != n:
raise ValueError(
f"Network '{self.name}' contains parallel (duplicate) links; "
"TABenchmark v0 requires unique (init_node, term_node) pairs "
"(see the known-defects registry for affected TNTP networks)."
)
if np.any(self.free_flow_time <= 0):
raise ValueError(
f"Network '{self.name}': free_flow_time must be strictly positive "
"(zero-cost links break shortest-path sparsity)."
)
if np.any(self.capacity <= 0):
raise ValueError(f"Network '{self.name}': capacity must be strictly positive")
if np.any(self.b < 0):
raise ValueError(
f"Network '{self.name}': BPR coefficient b (alpha) must be nonnegative "
"(negative b makes link costs decrease in flow, voiding Beckmann convexity)"
)
if np.any(self.power < 0):
raise ValueError(
f"Network '{self.name}': BPR exponent power (beta) must be nonnegative"
)
if np.any((self.init_node < 1) | (self.init_node > self.n_nodes)) or np.any(
(self.term_node < 1) | (self.term_node > self.n_nodes)
):
raise ValueError(f"Network '{self.name}': node ids out of range 1..n_nodes")
@property
def n_links(self) -> int:
return len(self.init_node)
@property
def fixed_cost(self) -> np.ndarray:
"""Flow-independent generalized-cost component per link."""
return self.toll_weight * self.toll + self.distance_weight * self.length
[docs]
def link_cost(self, link_flows: np.ndarray) -> np.ndarray:
"""Generalized link travel cost t_a(v_a) at the given flows."""
v = np.maximum(np.asarray(link_flows, dtype=np.float64), 0.0)
ratio = v / self.capacity
return self.free_flow_time * (1.0 + self.b * ratio**self.power) + self.fixed_cost
[docs]
def link_cost_integral(self, link_flows: np.ndarray) -> np.ndarray:
"""Per-link Beckmann integral: integral of t_a(s) ds from 0 to v_a."""
v = np.maximum(np.asarray(link_flows, dtype=np.float64), 0.0)
ratio = v / self.capacity
variable = self.free_flow_time * (v + self.b * v * ratio**self.power / (self.power + 1.0))
return variable + self.fixed_cost * v
[docs]
def link_cost_derivative(self, link_flows: np.ndarray) -> np.ndarray:
"""Per-link derivative dt_a/dv_a — the diagonal Hessian of the Beckmann
objective, used by conjugate-direction Frank-Wolfe variants.
For BPR, t'(v) = fft * b * p * (v/cap)^(p-1) / cap. Edge cases: p = 0
or b = 0 gives exactly 0 (constant cost); p = 1 gives the constant
fft*b/cap (numpy 0**0 = 1); 0 < p < 1 has an unbounded derivative at
v = 0, where H is defined as 0 (it is only a conjugacy scale) and
0**(negative) is never evaluated.
"""
v = np.maximum(np.asarray(link_flows, dtype=np.float64), 0.0)
ratio = v / self.capacity
exponent = self.power - 1.0
coeff = self.free_flow_time * self.b * self.power / self.capacity
# Zero wherever the analytic value is 0 (p or b zero) or defined as 0
# (v = 0 with p < 1) — without ever evaluating 0**(negative). For
# 0 < p < 1 at subnormal positive flows the analytic value exceeds
# float64 range: clamp to the largest finite float (H is only a
# conjugacy scale) so no inf/RuntimeWarning ever escapes.
zero = (coeff == 0.0) | ((exponent < 0) & (ratio <= 0.0))
safe = np.where(zero, 1.0, ratio)
with np.errstate(over="ignore"):
h = np.minimum(coeff * safe**exponent, np.finfo(np.float64).max)
return np.where(zero, 0.0, h)
[docs]
@dataclass(frozen=True)
class Demand:
"""Origin-destination demand. ``matrix[i, j]`` is flow from zone i+1 to j+1."""
matrix: np.ndarray # (n_zones, n_zones) float64
def __post_init__(self) -> None:
m = self.matrix
if m.ndim != 2 or m.shape[0] != m.shape[1]:
raise ValueError("Demand matrix must be square (n_zones, n_zones)")
if np.any(m < 0):
raise ValueError("Demand must be nonnegative")
@property
def n_zones(self) -> int:
return self.matrix.shape[0]
@property
def total(self) -> float:
return float(self.matrix.sum())
[docs]
@dataclass(frozen=True)
class ElasticDemand:
"""Elastic (variable) demand law: the demand between an OD pair is a
strictly decreasing function of that pair's equilibrium travel cost,
``d_rs = D_rs(u_rs)``.
The *reference* demand ``d0_rs = D_rs(0)`` (demand at zero cost — an upper
bound) is carried by the scenario's ordinary :class:`Demand` matrix; this
object adds only the decay law and its single shared parameter. Two forms:
* ``"linear"`` ``D_rs(u) = d0_rs * max(0, 1 - u/param)`` (``param`` = u0)
inverse on ``[0, d0]``: ``D_rs^{-1}(d) = param * (1 - d/d0)`` — bounded,
singularity-free (the excess-demand arc cost ``param * e/d0`` is linear).
* ``"exponential"`` ``D_rs(u) = d0_rs * exp(-param * u)`` (``param`` = beta)
inverse on ``(0, d0]``: ``D_rs^{-1}(d) = ln(d0/d) / param`` — the
cited/logit-consistent form, unbounded as ``d -> 0`` (guarded below).
Objective ``min Sum_a integral t_a - Sum_rs integral D_rs^{-1}`` and the
excess-demand ("dummy arc") transformation follow Sheffi (1985) ch. 6 and
Boyles, Lownes & Unnikrishnan, *Transportation Network Analysis* sec. 9.1;
the transformation itself is due to Gartner (1980), and the seminal
computational treatment is Florian & Nguyen (1974). See docs/design/adr-005.
``param`` is a scalar shared across OD pairs (per-OD curves are a future
extension). It is task data, content-hashed with the scenario when set.
"""
form: str
param: float
_FORMS = ("linear", "exponential")
def __post_init__(self) -> None:
if self.form not in self._FORMS:
raise ValueError(
f"ElasticDemand form must be one of {self._FORMS}, got {self.form!r}"
)
if not (np.isfinite(self.param) and self.param > 0):
raise ValueError(
f"ElasticDemand param must be finite and > 0, got {self.param!r}"
)
[docs]
def realized_demand(self, d0: np.ndarray, cost: np.ndarray) -> np.ndarray:
"""``D_rs(u)`` elementwise: realized demand at OD cost ``cost``, clamped
to ``[0, d0]``. Both arguments broadcast (typically ``(n_zones,
n_zones)`` reference-demand and OD shortest-path-cost matrices)."""
d0 = np.asarray(d0, dtype=np.float64)
u = np.asarray(cost, dtype=np.float64)
if self.form == "linear":
d = d0 * (1.0 - u / self.param)
else: # exponential
d = d0 * np.exp(-self.param * u)
return np.clip(d, 0.0, d0)
[docs]
def inverse_demand(self, d0: np.ndarray, d: np.ndarray) -> np.ndarray:
"""``D_rs^{-1}(d)``: the OD cost at which realized demand equals ``d``
(``d`` in ``[0, d0]``). Prices the excess-demand arc. For the
exponential form the ``d -> 0`` log-singularity is floored so the
transformed solver never evaluates ``ln(inf)``; ``d0 == 0`` pairs
(no demand) return 0."""
d0 = np.asarray(d0, dtype=np.float64)
d = np.asarray(d, dtype=np.float64)
positive = d0 > 0
if self.form == "linear":
return self.param * np.where(positive, 1.0 - d / np.where(positive, d0, 1.0), 0.0)
floor = np.maximum(d0, 1.0) * 1e-12
d = np.maximum(d, floor)
return np.where(positive, np.log(np.where(positive, d0, 1.0) / d) / self.param, 0.0)
[docs]
def excess_arc_cost(self, d0: np.ndarray, excess: np.ndarray) -> np.ndarray:
"""Cost of the excess-demand (Gartner dummy) arc carrying unmet flow
``excess``: ``W(e) = D_rs^{-1}(d0 - e)`` — increasing in ``e``."""
d0 = np.asarray(d0, dtype=np.float64)
return self.inverse_demand(d0, d0 - np.asarray(excess, dtype=np.float64))
[docs]
@dataclass(frozen=True)
class CombinedDemand:
"""Combined trip-distribution + assignment task data (Evans 1976).
Instead of a fixed OD matrix, the demand is *endogenous*: only the
trip-end margins are fixed — the productions ``O_i`` (trips originating in
zone ``i``) and attractions ``D_j`` (trips destined for zone ``j``) — and
the OD matrix ``d_ij`` is distributed across pairs by a **doubly-constrained
gravity model** at the equilibrium travel costs. The combined convex
program (Evans 1976; Sheffi, *Urban Transportation Networks* 1985 ch. 6;
Boyles, Lownes & Unnikrishnan, *Transportation Network Analysis* §6) is
min_{x,d} Σ_a ∫_0^{x_a} t_a(w) dw + (1/β) Σ_ij d_ij (ln d_ij − 1)
s.t. Σ_j d_ij = O_i, Σ_i d_ij = D_j, d ≥ 0, x = assign(d)
whose stationarity conditions give the doubly-constrained gravity
``d_ij = A_i B_j exp(−β u_ij)`` (``u_ij`` the equilibrium OD cost, ``A``/``B``
the balancing factors) *and* Wardrop route equilibrium of ``d``. The single
parameter ``beta`` is the gravity dispersion / impedance sensitivity (large
``beta`` → trips concentrate on cheap destinations; ``beta → 0`` → the
cost-free maximum-entropy distribution ``O_i D_j / T``).
Intrazonal trips never enter the network, so the gravity is distributed over
interzonal pairs only (``i ≠ j`` with ``O_i > 0`` and ``D_j > 0``). This is
task data, content-hashed with the scenario (docs/design/adr-007); the
reference :class:`Demand` matrix carries the free-flow gravity distribution
(a deterministic reference — the uncongested equilibrium demand). See
docs/design/adr-007-combined-distribution-assignment.md.
"""
productions: np.ndarray # (n_zones,) O_i, trips originating at each zone
attractions: np.ndarray # (n_zones,) D_j, trips attracted to each zone
beta: float
def __post_init__(self) -> None:
o = _as_f64(self.productions)
d = _as_f64(self.attractions)
object.__setattr__(self, "productions", o)
object.__setattr__(self, "attractions", d)
if o.ndim != 1 or d.ndim != 1 or o.shape != d.shape:
raise ValueError("CombinedDemand productions/attractions must be 1-D, equal length")
if np.any(o < 0) or np.any(d < 0):
raise ValueError("CombinedDemand productions/attractions must be nonnegative")
if not (np.isfinite(o).all() and np.isfinite(d).all()):
raise ValueError("CombinedDemand productions/attractions must be finite")
if not (np.isfinite(self.beta) and self.beta > 0):
raise ValueError(f"CombinedDemand beta must be finite and > 0, got {self.beta!r}")
# Doubly-constrained feasibility: total productions must equal total
# attractions, else no OD matrix reproduces both margins.
total_o, total_d = float(o.sum()), float(d.sum())
if not np.isclose(total_o, total_d, rtol=1e-9, atol=1e-9):
raise ValueError(
"CombinedDemand is infeasible: total productions "
f"({total_o}) must equal total attractions ({total_d})"
)
@property
def n_zones(self) -> int:
return self.productions.shape[0]
@property
def total(self) -> float:
return float(self.productions.sum())
[docs]
def support(self) -> np.ndarray:
"""Boolean ``(n_zones, n_zones)`` mask of distributable OD pairs:
interzonal (``i ≠ j``) pairs with positive production and attraction."""
o = self.productions
d = self.attractions
n = o.shape[0]
return (o[:, None] > 0) & (d[None, :] > 0) & ~np.eye(n, dtype=bool)
[docs]
def gravity(
self, od_cost: np.ndarray, *, tol: float = 1e-13, max_iters: int = 10000
) -> np.ndarray:
"""Doubly-constrained gravity demand ``d_ij = A_i B_j exp(−β u_ij)`` at
OD costs ``u = od_cost`` (a ``(n_zones, n_zones)`` matrix; only support
entries are read), balanced to the productions/attractions by the
classic Furness / iterative-proportional-fitting recursion
A_i = O_i / Σ_j B_j f_ij, B_j = D_j / Σ_i A_i f_ij,
f_ij = exp(−β u_ij) on the interzonal support (0 elsewhere).
Deterministic (fixed tolerance/iteration cap) so the harness certificate
and the solver recompute byte-identical demand from the same costs (P1).
A final row rescale makes the production margins exact; the attraction
margins converge to ``tol``. Returns a dense matrix, zero off support.
"""
o = self.productions
d = self.attractions
support = self.support()
u = np.asarray(od_cost, dtype=np.float64)
# f_ij = exp(-beta u_ij) on support. u >= fft > 0 there, so the exponent
# is strictly negative and exp never overflows; off support f = 0. When
# beta * u is so large (>~708) that exp underflows, every support entry
# of a producing row could otherwise collapse to 0, silently dropping
# that row's production margin (and the scored realized_demand). Floor
# the on-support deterrence at the smallest normal float so no producing
# row can lose its margin; the floor is ~1e-308 and never binds in the
# normal regime, so it changes nothing there (it only degrades
# gracefully to a near-uniform split when all relative weights underflow).
tiny = np.finfo(np.float64).tiny
f = np.where(support, np.maximum(np.exp(-self.beta * np.where(support, u, 0.0)), tiny), 0.0)
active_o = o > 0
active_d = d > 0
a = np.ones(o.shape[0], dtype=np.float64)
b = np.ones(d.shape[0], dtype=np.float64)
for _ in range(max_iters):
row = f @ b
a = np.where(active_o & (row > 0), o / np.where(row > 0, row, 1.0), 0.0)
col = f.T @ a
b_new = np.where(active_d & (col > 0), d / np.where(col > 0, col, 1.0), 0.0)
if np.max(np.abs(b_new - b)) <= tol:
b = b_new
break
b = b_new
# Final row solve makes Σ_j d_ij = O_i exact (attractions within tol).
row = f @ b
a = np.where(active_o & (row > 0), o / np.where(row > 0, row, 1.0), 0.0)
return (a[:, None] * b[None, :]) * f
[docs]
@dataclass(frozen=True)
class MulticlassDemand:
"""Multiclass-user demand + linear class-interaction (Dafermos 1972; adr-013).
``K >= 2`` user classes share one physical network. Each class ``i`` has its
own OD demand matrix ``matrices[i]`` (so the classes couple only through the
link costs, never the demand constraints — the feasible set is a Cartesian
product of per-class demand polytopes) and perceives a class-specific link
cost that is the shared BPR latency of the **total** link flow plus a linear
class interaction
t_a^i(V) = t_a^BPR(v_a) + sum_j interaction[i, j] * v_a^j,
v_a = sum_j v_a^j (total flow on link a), v_a^j = class j's flow on a.
``interaction`` is a ``(K, K)`` matrix applied per link (units of cost per
unit class-``j`` link flow). It is what makes the split among classes
well-defined: with ``interaction = 0`` every class sees the identical cost
``t_a^BPR(v_a)`` and the per-class split is arbitrary (the model degenerates
to a single-class UE on the summed demand). A **symmetric** ``interaction``
(``interaction[i,j] == interaction[j,i]``) makes the equilibrium the
minimizer of a convex multiclass-Beckmann potential (the integrable case
Dafermos 1972 characterizes); an **asymmetric** one makes it a genuine
variational inequality with no equivalent optimization (Smith 1979; Dafermos
1980) — the block-structured generalization of the single-class ``vi-asym``
cost ``t(v) = t_BPR(v) + C v``.
This is task data, content-hashed with the scenario when set. No sign
constraint is imposed on ``interaction`` (a pathological entry that drives a
cost non-positive is caught by the solver and censored by the certificate,
exactly as for ``vi-asym``)."""
matrices: np.ndarray # (K, n_zones, n_zones) float64, per-class OD demand
interaction: np.ndarray # (K, K) float64, linear class-coupling applied per link
def __post_init__(self) -> None:
m = _as_f64(self.matrices)
c = _as_f64(self.interaction)
object.__setattr__(self, "matrices", m)
object.__setattr__(self, "interaction", c)
if m.ndim != 3 or m.shape[1] != m.shape[2]:
raise ValueError(
"MulticlassDemand matrices must have shape (n_classes, n_zones, n_zones)"
)
if m.shape[0] < 2:
raise ValueError(
"MulticlassDemand needs >= 2 classes (a single class is an ordinary "
"UE; use Demand); got n_classes="
f"{m.shape[0]}"
)
if np.any(m < 0) or not np.all(np.isfinite(m)):
raise ValueError("MulticlassDemand matrices must be finite and nonnegative")
if c.shape != (m.shape[0], m.shape[0]):
raise ValueError(
f"MulticlassDemand interaction must have shape ({m.shape[0]}, "
f"{m.shape[0]}), got {c.shape}"
)
if not np.all(np.isfinite(c)):
raise ValueError("MulticlassDemand interaction must be finite")
@property
def n_classes(self) -> int:
return self.matrices.shape[0]
@property
def n_zones(self) -> int:
return self.matrices.shape[1]
@property
def total_matrix(self) -> np.ndarray:
"""Class-summed OD matrix — the aggregate demand the network carries."""
return self.matrices.sum(axis=0)
@property
def total(self) -> float:
return float(self.matrices.sum())
[docs]
def symmetric(self, tol: float = 1e-12) -> bool:
"""True iff the class interaction is symmetric (the integrable case)."""
return bool(np.allclose(self.interaction, self.interaction.T, atol=tol, rtol=0.0))
[docs]
@dataclass(frozen=True)
class ReferenceSolution:
"""Best-known solution used as a regression oracle (with provenance)."""
link_flows: np.ndarray
source: str
note: str = ""
[docs]
@dataclass(frozen=True)
class Scenario:
"""A frozen benchmark instance: network + demand (+ optional oracle).
``family`` names the scenario lineage used by the fairness gate (P7):
a learned model whose ``trained_on`` includes this family is refused
evaluation on this scenario.
``sue_theta`` (optional) makes this an SUE task: the dispersion dial of the
pinned loading map, in 1/(native cost unit) — so its meaning is
network-specific (P9; scenario cards state the unit). It is task data,
never a model factor, and it is content-hashed when set (two scenarios
differing only in theta are different benchmark instances).
``sue_family`` selects the choice-model family of that SUE task —
``"logit"`` (Dial-STOCH closed-form loading, docs/design/adr-001) or
``"probit"`` (Monte Carlo loading + pinned certificate, adr-003). For
probit, ``sue_theta`` carries ``beta``, the perception variance per unit
free-flow time. This is unrelated to ``family``, which is *data lineage*
for the ``trained_on`` fairness gate (P7); ``sue_family`` names the task's
equilibrium definition. It is hashed only when non-default so every logit
scenario keeps the byte-identical hash it had before this field existed.
``elastic_demand`` (optional) makes this a variable-demand UE task: the
:class:`Demand` matrix becomes the *reference* demand ``d0 = D(0)`` and this
:class:`ElasticDemand` supplies the decay law (docs/design/adr-005). Like
``sue_theta`` it is task data content-hashed only when set, so every
fixed-demand scenario keeps its prior hash. An elastic task is a
deterministic UE with endogenous demand and so is mutually exclusive with
the SUE fields.
``combined_demand`` (optional) makes this a combined trip-distribution +
assignment task (Evans 1976; docs/design/adr-007): the OD matrix is
endogenous, distributed by a doubly-constrained gravity model from the fixed
trip-end margins in this :class:`CombinedDemand`, and the :class:`Demand`
matrix carries the free-flow gravity *reference*. Like the other optional
task fields it is content-hashed only when set. It is mutually exclusive
with both the SUE fields and ``elastic_demand`` (all three make demand
non-fixed in incompatible ways).
``br_epsilon`` (optional) makes this a **boundedly-rational** user-equilibrium
task (Mahmassani & Chang 1976/1987; docs/design/adr-008): a flow is acceptable
if every *used* route lies within an absolute indifference band ``epsilon``
(native cost units) of its OD's minimum route cost, ``c_pi <= kappa_rs +
epsilon``. This relaxes Wardrop's equality to a one-sided band, so the
equilibrium is a *set*, not a point (``epsilon=0`` recovers UE). It is task
data, content-hashed only when set, and mutually exclusive with the SUE /
elastic / combined fields (BR-UE is a deterministic fixed-demand route
equilibrium; those make the demand non-fixed). The band unit is
network-specific (scenario cards state it, P9).
``side_capacities`` (optional) makes this a **side-constrained** UE task
(Larsson & Patriksson 1995; docs/design/adr-009): a per-link hard capacity
``u_a`` on the link flow, ``v_a <= u_a``, imposed *in addition to* the BPR
latency (these are physical throughput limits, distinct from the BPR
``Network.capacity`` reference volume). The equilibrium is Wardrop UE on the
capacity-augmented cost ``t_a(v_a) + beta_a`` with ``beta_a >= 0`` a queueing
delay / toll that is zero off the binding set. It is per-link task data,
content-hashed only when set, and mutually exclusive with the SUE / elastic /
combined / BR fields.
``multiclass`` (optional) makes this a **multiclass-user** equilibrium task
(Dafermos 1972; docs/design/adr-013): ``K >= 2`` user classes share the
network, each with its own OD demand and a class-specific link cost coupling
the classes (:class:`MulticlassDemand`). The aggregate ``demand`` must equal
the class sum (so single-class consumers of ``demand`` stay consistent). It
is task data, content-hashed only when set, and mutually exclusive with the
SUE / elastic / combined / BR / side-constrained / link-interaction fields
(``link_interaction`` is its single-class special case). Certifying it needs
the model's per-class link flows (``FlowState.class_link_flows``); the scored
quantity is the class-summed VI residual.
"""
name: str
network: Network
demand: Demand
reference: ReferenceSolution | None = None
family: str = field(default="")
sue_theta: float | None = None
sue_family: str = "logit"
elastic_demand: ElasticDemand | None = None
combined_demand: CombinedDemand | None = None
br_epsilon: float | None = None
side_capacities: np.ndarray | None = None
link_interaction: np.ndarray | None = None
multiclass: MulticlassDemand | None = None
def __post_init__(self) -> None:
if self.demand.n_zones != self.network.n_zones:
raise ValueError(
f"Scenario '{self.name}': demand has {self.demand.n_zones} zones, "
f"network declares {self.network.n_zones}"
)
if not self.family:
object.__setattr__(self, "family", self.name)
if self.sue_theta is not None and not (
np.isfinite(self.sue_theta) and self.sue_theta > 0
):
raise ValueError(
f"Scenario '{self.name}': sue_theta must be finite and > 0, "
f"got {self.sue_theta!r}"
)
if self.sue_family not in ("logit", "probit"):
raise ValueError(
f"Scenario '{self.name}': sue_family must be 'logit' or 'probit', "
f"got {self.sue_family!r}"
)
if self.sue_family == "probit" and self.sue_theta is None:
raise ValueError(
f"Scenario '{self.name}': sue_family='probit' requires sue_theta "
"(beta, the perception variance per unit free-flow time)"
)
if self.elastic_demand is not None and self.sue_theta is not None:
raise ValueError(
f"Scenario '{self.name}': elastic_demand and sue_theta are "
"mutually exclusive (elastic UE is deterministic with endogenous "
"demand; SUE is stochastic with fixed demand)"
)
if self.combined_demand is not None:
if self.combined_demand.n_zones != self.network.n_zones:
raise ValueError(
f"Scenario '{self.name}': combined_demand has "
f"{self.combined_demand.n_zones} zones, network declares "
f"{self.network.n_zones}"
)
if self.sue_theta is not None or self.elastic_demand is not None:
raise ValueError(
f"Scenario '{self.name}': combined_demand is mutually exclusive "
"with sue_theta and elastic_demand (each makes the OD demand "
"non-fixed in an incompatible way)"
)
if self.br_epsilon is not None:
if not (np.isfinite(self.br_epsilon) and self.br_epsilon > 0):
raise ValueError(
f"Scenario '{self.name}': br_epsilon must be finite and > 0, "
f"got {self.br_epsilon!r}"
)
if (
self.sue_theta is not None
or self.elastic_demand is not None
or self.combined_demand is not None
):
raise ValueError(
f"Scenario '{self.name}': br_epsilon is mutually exclusive with "
"sue_theta, elastic_demand and combined_demand (BR-UE is a "
"deterministic fixed-demand route equilibrium)"
)
if self.side_capacities is not None:
u = _as_f64(self.side_capacities)
object.__setattr__(self, "side_capacities", u)
if u.shape != (self.network.n_links,):
raise ValueError(
f"Scenario '{self.name}': side_capacities must have shape "
f"({self.network.n_links},), got {u.shape}"
)
if np.any(u <= 0) or not np.all(np.isfinite(u)):
raise ValueError(
f"Scenario '{self.name}': side_capacities must be finite and > 0"
)
if (
self.sue_theta is not None
or self.elastic_demand is not None
or self.combined_demand is not None
or self.br_epsilon is not None
):
raise ValueError(
f"Scenario '{self.name}': side_capacities is mutually exclusive "
"with sue_theta, elastic_demand, combined_demand and br_epsilon"
)
if self.link_interaction is not None:
c = _as_f64(self.link_interaction)
object.__setattr__(self, "link_interaction", c)
m = self.network.n_links
if c.shape != (m, m):
raise ValueError(
f"Scenario '{self.name}': link_interaction must have shape "
f"({m}, {m}), got {c.shape}"
)
if not np.all(np.isfinite(c)):
raise ValueError(
f"Scenario '{self.name}': link_interaction must be finite"
)
if (
self.sue_theta is not None
or self.elastic_demand is not None
or self.combined_demand is not None
or self.br_epsilon is not None
or self.side_capacities is not None
):
raise ValueError(
f"Scenario '{self.name}': link_interaction is mutually exclusive "
"with sue_theta, elastic_demand, combined_demand, br_epsilon and "
"side_capacities (it is a deterministic fixed-demand VI with "
"non-separable link costs)"
)
if self.multiclass is not None:
mc = self.multiclass
if mc.n_zones != self.network.n_zones:
raise ValueError(
f"Scenario '{self.name}': multiclass demand has {mc.n_zones} "
f"zones, network declares {self.network.n_zones}"
)
# The aggregate `demand` must equal the class sum, so the network's
# total loading and every single-class consumer of `demand` (feasibility
# scale, fairness gate) stay consistent with the per-class task.
if not np.allclose(mc.total_matrix, self.demand.matrix, rtol=1e-9, atol=1e-9):
raise ValueError(
f"Scenario '{self.name}': demand matrix must equal the multiclass "
"class sum (multiclass.matrices.sum(axis=0)); pass "
"Demand(multiclass.total_matrix) as the aggregate demand"
)
if (
self.sue_theta is not None
or self.elastic_demand is not None
or self.combined_demand is not None
or self.br_epsilon is not None
or self.side_capacities is not None
or self.link_interaction is not None
):
raise ValueError(
f"Scenario '{self.name}': multiclass is mutually exclusive with "
"sue_theta, elastic_demand, combined_demand, br_epsilon, "
"side_capacities and link_interaction (multiclass demand is a "
"distinct per-class equilibrium; link_interaction is its "
"single-class special case)"
)
[docs]
def content_hash(self) -> str:
"""SHA-256 over the canonical serialization of all scored content (P2)."""
h = hashlib.sha256()
net = self.network
h.update(
f"nodes={net.n_nodes};zones={net.n_zones};ftn={net.first_thru_node};"
f"tw={net.toll_weight!r};dw={net.distance_weight!r};".encode()
)
for label, arr in (
("init", net.init_node),
("term", net.term_node),
("cap", net.capacity),
("len", net.length),
("fft", net.free_flow_time),
("b", net.b),
("pow", net.power),
("toll", net.toll),
("od", self.demand.matrix),
):
h.update(label.encode())
h.update(_as_f64(arr).tobytes())
# Conditional, appended last: scenarios without a theta hash exactly
# as before this field existed (pinned by a golden-hash test).
if self.sue_theta is not None:
h.update(f"sue_theta={float(self.sue_theta)!r};".encode())
# Appended after theta and only when non-default, so every logit
# scenario hashes exactly as before this field existed (golden test):
# a probit task can never collide with the logit task at the same theta.
if self.sue_family != "logit":
h.update(f"sue_family={self.sue_family};".encode())
# Appended last and only when set: every fixed-demand scenario hashes
# exactly as before this field existed (golden-hash test). The
# reference demand d0 is already covered by the ("od", matrix) bytes
# above; only the decay law and its parameter are new content.
if self.elastic_demand is not None:
ed = self.elastic_demand
h.update(f"elastic_form={ed.form};elastic_param={float(ed.param)!r};".encode())
# Appended last and only when set: every scenario without a combined
# task hashes exactly as before this field existed (golden-hash test).
# The reference free-flow-gravity matrix is already covered by the
# ("od", matrix) bytes above; only the fixed margins and dispersion are
# new scored content.
if self.combined_demand is not None:
cd = self.combined_demand
h.update(f"combined_beta={float(cd.beta)!r};".encode())
h.update(b"combined_prod")
h.update(_as_f64(cd.productions).tobytes())
h.update(b"combined_attr")
h.update(_as_f64(cd.attractions).tobytes())
# Appended last and only when set: every non-BR scenario hashes exactly as
# before this field existed (golden-hash test). Two scenarios differing
# only in the indifference band are different benchmark instances.
if self.br_epsilon is not None:
h.update(f"br_epsilon={float(self.br_epsilon)!r};".encode())
# Per-link hard capacities are scored instance content (a change in any
# u_a is a different benchmark instance); hashed only when set.
if self.side_capacities is not None:
h.update(b"side_cap")
h.update(_as_f64(self.side_capacities).tobytes())
# Non-separable link-cost interaction operator C (t(v) = t_BPR(v) + C v);
# a possibly-asymmetric C is scored instance content, hashed only when set,
# so every scenario without it hashes exactly as before (golden-hash test).
if self.link_interaction is not None:
h.update(b"link_interaction")
h.update(_as_f64(self.link_interaction).tobytes())
# Multiclass per-class demand + class-interaction is scored instance
# content; appended last and only when set, so every scenario without it
# (the aggregate `demand` bytes above are the class sum) hashes exactly as
# before this field existed (golden-hash test).
if self.multiclass is not None:
h.update(b"multiclass_matrices")
h.update(_as_f64(self.multiclass.matrices).tobytes())
h.update(b"multiclass_interaction")
h.update(_as_f64(self.multiclass.interaction).tobytes())
return h.hexdigest()