Structural analysis: infer the coupling graph from the IR (separable dims, linking constraints)

[FBumann]

Where this comes from

Prior-art scan of declarative optimisation languages surfaced GBOML (U. Liège, JOSS 2022) — a stand-alone MILP language whose organising abstraction is a hierarchical hypergraph: #NODE blocks own internal variables and local constraints, #HYPEREDGE blocks carry the coupling constraints, and coupling may only reference the external: variables of incident nodes. Their JOSS paper diagnoses a gap we recognise — algebraic modeling languages "rarely provide language constructs to represent the block structure in a model and typically fail to take advantage of it", while framework-style tools (they name PyPSA, Calliope, Balmorel) "lack expressiveness and adding components is often cumbersome". They claim four payoffs from carrying that structure: simpler encoding, model re-use, parallelised model generation, and structure-exploiting solvers.

The observation that motivates this issue: we have the same structural information, and we do not have to ask the user for it. GBOML makes the modeller draw the hypergraph by hand. In our IR the coupling is already marked syntactically, and can be recovered by a static pass over Program.

This issue is only the analysis-and-reporting half. Consuming the structure (block-annotated sinks, decomposition solvers) is explicitly out of scope — see the bottom.

The claim

relational/ir.py Program is a declaration-level object: a handful of ConstraintDecls, each with dims and an affine Expr tree. Nothing is flattened, so the coupling pattern survives — and the expression nodes that couple coords are exactly the ones ARCHITECTURE.md already classifies for executor locality:

IR node	coupling
Add, Mul, Div, Neg, Var, Param, Cmp masks	none — pointwise, row (g1, t5) touches only variables at (g1, t5)
Shift(dim=d, n=k)	banded along d with bandwidth k (bounded halo)
Sum(over=d)	collapses d — a linking row across all of d
GroupSum(mapping=m, into=i)	collapses the mapped dim into i — linking, but partitioned by the mapping data

From which the rule: a dim that appears in some dims: but never inside a Sum.over, a GroupSum, or a Shift.dim anywhere in the program is separable — the model is block-diagonal along it. This is sound by construction (every coupling route is one of the nodes above) and conservative in the right direction: it can under-report separability, never over-report it.

Caveat worth encoding in the output: Sum and Shift are statically decidable, but GroupSum coupling is data-dependent. Whether gen_bus links everything into one component or partitions into five islands lives in the mapping table. So the pass has two tiers — a static verdict from the Program alone, and an optional refinement that runs connected components over the mapping table(s). The refinement is a query over one small parameter, not over the constraint matrix.

What to build

1. A pure function over Program — relational/structure.py, no engine imports, mirroring ir.py's "frozen dataclasses, no execution logic" discipline — returning:
   • separable dims (block-diagonal axes) with the block count where statically known,
   • linking constraint declarations, and which dim each one collapses,
   • time-coupling bandwidth per dim (max |n| over Shift nodes on that dim),
   • the declaration-level bipartite incidence: which constraint families touch which variable families.
2. Optional data-tier refinement for GroupSum: connected components over the mapping parameter, giving the real block count.
3. Surface it. check() is the natural home — it is already the "tell me about this file without building it" entry point, and this needs no data. Overlaps with Build/solve observability: executor stats #34 (executor stats); structure is a pre-build property, so it should be reportable from check() alone and echoed in stats when a build happens.

Output should read as a diagnostic a user acts on, e.g. "block-diagonal along region (14 blocks); linking constraints: co2_budget (sums region), transmission_balance (group_sum line_region → region); snapshot coupled with bandwidth 1 (storage_balance)."

Why this is worth it independent of decomposition

Diagnostics alone justify it — knowing your 14 regions are independent except for one CO₂ budget row is a real modelling fact, and today nothing tells you.

But the reason to build it now is that separability is the prerequisite for sequential per-block build-and-solve. If a dim is provably separable, the engine can build and solve one slice at a time, holding a single block in memory. That is an out-of-core scaling story that falls straight out of hard rule 4 rather than straining it, needs no external solver, and is a strictly better use of the structure than exporting it.

Explicitly out of scope

• Block-annotated sinks (SMPS-with-blocks, DSP format). Structurally this is just another sink and the COO batches are already there — but it bets on a small, mostly academic structure-exploiting-solver ecosystem. Separate issue if ever.
• Implementing decomposition (Benders, Dantzig-Wolfe). These are solve algorithms — cut management, convergence, subproblem lifecycle — requiring a Python driver loop around highspy. That is Python-side algorithmic machinery of exactly the kind hard rule 6 exists to keep out. Not a goal.

Definition of done

check() reports the coupling structure of a model with no data bound; tests cover a block-diagonal program, a program with one linking Sum, a Shift-banded program, and a GroupSum program where the static verdict says "coupled" and the data-tier refinement finds N components. ARCHITECTURE.md's locality section notes that the same classification also yields the coupling graph.

🤖 Generated with Claude Code

[FBumann]

Parking this rather than scheduling it:

• the description scopes it to the analysis-and-reporting half, and no consumer for the structure exists yet
• the natural consumers are Namespacing: qualified names for component libraries #29 / Compose-then-build: schema merge as a pure step before one lower/stream pass #30 (component libraries), neither of which has resolved

Deprioritisation is expressed by the absence of the roadmap label rather than by a new one, so this stays open as backlog. The GBOML prior-art reasoning is the part worth keeping.

[FBumann]

ROADMAP Track 2d (#378) names this as the block-detection step for decomposition: the master/sub split is a GROUP BY over the A frame, which is what this analysis produces. Sequenced with #204.

[FBumann]

Path fix: this issue is written against relational/ir.py and Program. The module is relational/plan.py and the type is plan.Program; the frozen-dataclass, no-execution-logic discipline it cites is intact under the new name.

The node table also needs one substitution — Shift(dim, n) is now plan.Translate, one node for both surface spellings differing only in wrap. That strengthens the claim rather than weakening it: bandwidth is |by| off a single node kind instead of two, and wrap=True is the case where the band is cyclic, which the coupling report would want to say.