Decision before v1: first-class MultiIndex support vs. flat dims + auxiliary level coords

[FBumann]

Note

AI-written analysis (Claude Code, prompted by @FBumann). This is a scope question for #717 — the decision window closes when v1 ships, because afterwards any change here is a breaking change with a deprecation cycle.

The question

Every difficult problem in the #732 → #737 → #742 → #717 chain traces back to one root: linopy supports stacked pd.MultiIndex dimensions as a first-class data model. Should v1 keep that — or replace it with a representation that needs none of the machinery?

This is the radical version of the "learn from the xarray community, which struggled a lot with MIs" argument that decided scenario B.

What MultiIndex support costs

Complexity that exists only for MI
_project_onto_multiindex_levels, _LevelProjection, _as_multiindex	~120 lines
MI branches in _coords_to_dict, validate_alignment, _broadcast_to_coords (expand-via-template)	~80 lines
assign_multiindex_safe (~20 call sites), get_dims_with_index_levels, MI serialization in coords_to_dataset_vars	~100 lines
§11's stacked-MI paragraph, the scenario A/B design discussion, pydata/xarray#11368 workarounds	weeks
TestMultiIndexProjection + MI tests across 6 files	~400 lines

Plus a permanent tax: every future feature must answer "and what about MultiIndex?".

The alternative: flat dim + auxiliary level coords

The same information, no pd.MultiIndex anywhere — and §11 already governs auxiliary coords:

# instead of:  snapshot = MultiIndex[(2020,t1), (2020,t2), (2030,t1), (2030,t2)]
snapshots = pd.RangeIndex(4, name="snapshot")
period    = xr.DataArray([2020, 2020, 2030, 2030], dims="snapshot", ...)
timestep  = xr.DataArray(["t1", "t2", "t1", "t2"], dims="snapshot", ...)

Verified against the current #717 branch (all snippets run):

x = m.add_variables(coords=[snapshots], name="x")              # ✅ works
expr = (1 * x).assign_coords(period=period)                    # ✅ aux coords attach

# per-period weighting — same explicit recipe as the MI case:
w = xr.DataArray(weights[period.values].values, dims="snapshot", ...)
expr * w                                                       # ✅ works, no projection machinery involved

# groupby a level:
expr.drop_vars("period").groupby(period.rename("inv_period")).sum()   # ⚠️ works, but needs the
                                                                      #    drop/rename dance — naming
                                                                      #    conflict otherwise (fixable)

Sparse indexes (not every combination exists) are naturally representable — that's what a flat list is. The stacked/unstacked round-trip, the projection, the coverage-gap concept: none of them exist in this representation.

Three options

	Internal complexity	PyPSA impact	Tuple .sel() / .unstack()
A. Disallow MI (TypeError → point to flat+aux)	deleted	n.snapshots API migration	gone
B. Convert at the boundary (MI accepted as input, stored flat+aux, re-stacked on output)	mostly deleted	thin adapter at model build + solution extraction	gone on linopy objects
C. Status quo (#717 as-is)	stays forever	none	works

What B would mean concretely

• coords=[multiindex] still works — linopy decomposes it into flat dim + level aux-coords on entry.
• The §11 stacked-MI paragraph reduces to one sentence ("MultiIndex coords are stored as a flat dimension with level coords; §11 governs the levels").
• solution / dual come back flat-indexed with level coords as columns; PyPSA re-stacks for its users (one set_index(levels) call).
• The groupby naming rough edge needs fixing (small).
• Cannot reindex onto a stacked MultiIndex via indexers — only reindex_like works pydata/xarray#11368 stops mattering to linopy entirely.

Why this needs deciding now

#717 currently implements C — including the machinery that the legacy-removal checklist says survives 1.0. If the answer is A or B, that machinery (and the spec section, and the tests) should not ship in v1 at all. After v1 ships with C, moving to A/B is a user-facing breaking change.

No position is taken here — the trade-off is real on both sides (PyPSA's n.snapshots API is the crux). But it should be a decision, not a default.

[FBumann]

@FabianHofmann I think this could really simplify things...

[FabianHofmann]

good points and right timing to raise that. I have to do some experiments with pypsa (there is another MI for stochastic optimization), but I have the strong feeling we should go with A or B.

[FBumann]

Same

[FBumann]

Note

Written with an agent (Claude Code, prompted by @FBumann). All snippets were run against the current branch and PyPSA 1.0.1.

Follow-up findings

Two things from digging into the code that sharpen the A-vs-B-vs-C call.

1. Arithmetic does not anchor C — v1 already drops implicit MI projection

A natural worry about A is "won't arithmetic break if it stops accepting MultiIndexes?" It won't, because v1 has already decided to remove the only MI-specific behaviour arithmetic has. expr * <per-period weight> today routes through _project_onto_multiindex_levels and emits:

multiindex-projection: implicitly broadcasting level subset ['period'] onto MultiIndex dimension 'snapshot'. This is deprecated and will raise under the v1 convention; project the input onto the dimension explicitly.

So under v1 the implicit-projection path raises regardless of A/B/C. Once it's gone, arithmetic with MI decomposes into:

case	under flat+aux (A)
operands share the full index	plain dim alignment — xarray-native, zero linopy machinery
operand is a level subset (per-period weight)	user projects explicitly — the same code v1 already requires under C

The projection engine's only remaining job in arithmetic under v1 is to emit the deprecation and then raise. That's not a feature v1 endorses — so giving it up costs nothing the spec wasn't already taking away.

2. PyPSA already bypasses linopy's MI groupby, and selects on levels mid-build

PyPSA's hardest MI paths don't use the C machinery anyway. The storage-SOC and store-energy per-period constraints hand-roll the rolling because linopy/xarray MI groupby is broken (pydata/xarray#6836):

# pypsa/optimization/constraints.py — "Normally we should use groupby, but it's broken for multi-index"
previous = [soc.data.sel(snapshot=(p, slice(None))).roll(snapshot=1) for p in sns.unique("period")]

This matters for option B: that sel(snapshot=(period, slice)) is a tuple selection on the stored MI, done mid-build, not at the I/O boundary. Verified:

	x.data.sel(snapshot=(2020, slice(None)))
C (MI stored)	✅
flat+aux	❌ KeyError → rewrite to .where(period==2020, drop=True)

So B's "MI in at input, flat+aux internally, re-stack on output" does not spare PyPSA the migration — these call sites break under B exactly as under A. B's transparency advantage is narrower than it looks; the clean split is A vs C.

What A requires in user code

Three concrete migrations (all have a verified one-liner equivalent):

# (a) per-level weighting — was implicit, now explicit (v1 requires this under C too)
expr * weights_by_period                          # was: implicit projection (deprecated)
expr * xr.DataArray(weights[period.values], dims="snapshot")   # A

# (b) tuple selection on a level
obj.sel(snapshot=(2020, slice(None)))             # C
obj.where(period == 2020, drop=True)              # A  (or .sel on the flat dim by mask)

# (c) groupby a level
expr.groupby(period).sum()                        # rough edge today (KeyError)
expr.drop_vars("period").groupby(period.rename("inv_period")).sum()   # A, works now; naming edge fixable

For PyPSA specifically: n.snapshots as a pd.MultiIndex is the crux — it's public API for multi-investment-period and stochastic networks, so A is a coordinated PyPSA-side change (migrate n.snapshots to a flat dim + period/scenario level coords), not a pure linopy decision.

Suggestion

If input ergonomics are the only reason to keep MI, the lean option is A + input-only sugar: accept coords=[multiindex], decompose to flat dim + level coords on entry, and never round-trip back. That deletes the §3 machinery (the ~40 assign_multiindex_safe call sites, the projection engine) while keeping coords=[mi] working — without B's permanent conversion adapter or its false promise of transparency.

[FBumann]

@FabianHofmann This is the one blocking v1 semantics imo. But its not trivial. Ill open a PR documenting a discussion around it, including reasoning and tests.