Tracking: sparse/CSR path follow-ups after #961

[FabianHofmann]

Note

The following content was generated by AI.

Describe the feature you'd like to see

Tracking issue for the CSR follow-ups after #961 (umbrella #756 is closed). Each child names a place where a CSR-backed expression silently densifies, or a gap in the frozen-constraint API. This issue fixes the order in which to address them.

Key observations

Line references in this section point to the code before PR 1.

• LinearExpression.data (linopy/expressions.py:2358) is where a CSR-backed expression densifies on read. Other paths drop the backing without passing through it: @ returns dense when the input was dense and sparse_groupby is off (expressions.py:2609), CSRConstraint.to_dense()/mutable(), the rebuild at model.py:1345, and the return None fallbacks in _sparse_matmul, _try_csr_merge, _aligned and csr_rhs, which lose the reason.
• Grid holds only per-dim indexes. Aux coords live next to it in CSRLinearExpression.coords and are dropped by CSRConstraint.from_csr (constraints.py:1401), which is CSRConstraint drops auxiliary coordinates when a sparse expression becomes a constraint #941.
• Zero policy on the sparse path: every operation keeps explicit zeros (merge, scaling, sum, groupby, selection), so sparse and dense give the same terms up to order; only @/dot prunes them. Cell activeness is carried by const alone. Export culls zeros regardless; Zero culling forces model rebuilding #925 is an export/persistent-diff question and does not block this series.

Plan

Four pull requests, grouped by shared design rather than one per issue:

1. Expression kernel (parts A, C, B, D), feat(csr): keep CSR backing through metadata, constant arithmetic, sum/groupby and selection #973: observability, aux coords in Grid, and all operations that keep the backing. Merged.
2. Boundaries (parts E, F), feat(csr): constraint and export boundaries of the sparse path (PR 2 of #972) #974: constraint side and export. Merged.
3. Sparse soften (part S), feat(csr): native soften for frozen constraints #979 (closes Native sparse soften for frozen constraints (penalty= with freeze=True) #975): native soften / penalty= on frozen constraints, so the model key in PR 4 needs no special case for softening. Merged.
4. Control (part G): the read-only model key from One Model-level switch for the sparse path instead of three independent opt-ins #976, feat(csr): one read-only Model(sparse=True) key for the sparse path #985. Merged. The docs (part H) moved to PR 9 in phase 2.

Part	Issues	Content	Depends on
A: observability	#962, #969 (parts 1+2)	Serve shape/sizes/coords/dims/isnull from the CSR store. .is_sparse, repr marker. Opt-in warn_on_densify through one _densify_notice(reason) helper, called at data, at each return None fallback and in CSRConstraint.to_dense. data stays a one-way conversion (no cached dense copy next to _csr, the setters would make it stale). Move aux coords into Grid.	–
C	#965	scaled_by / shifted for non-scalar constants, aligned through _matmul_operand_to_matrix	A
B	#964	sum(dim) by merging rows directly (CSRLinearExpression.aggregated), keeping explicit zeros. Chained groupby from a CSR source; fix the gate at expressions.py:613 to use coord_dims instead of self.data.dims.	A
D	#966, mask= part of #970	Row masking, sel through reindexed, isel as row gather. Keep _try_csr_merge sparse when aux coords differ across grids.	A
E: constraint side	#941, #963, rest of #970	Aux coords on CSRConstraint (trivial once in Grid), incl. netcdf round trip (io.py:1260). Explanatory errors for loc/update/from_rule naming .mutable(). Expression rhs moved to lhs.	D, #806
F: export	#968, then #967	flat/to_polars from the CSR store, modelled on CSRConstraint.to_polars. Sparse objective across objective.py, matrices.py, io.py, persistent/diff.py; decide where the objective name lives (CSRLinearExpression has no attrs).	B
S: sparse soften	#975, PR #979	Soften a CSRConstraint natively (one slack term per active row), drop the penalty + freeze error	E, F
G: control	#969 (part 3), #976	One read-only model key instead of the three opt-ins (semantics="v1", sparse=/sparse_groupby, freeze_constraints), as proposed in #976 (under discussion). Replaces the earlier plan of per-call sparse= on @/dot/merge; removes the coupling of @ to options["sparse_groupby"].	S, #976
H: docs	#971	API entries, user-guide section, list of operations that keep the backing	G

Ordering notes:

• Within PR 1, commit in the order A, C, B, D. PR 1 stays out of constraints.py so fix: support adding/removing variables after adding frozen constraint #806 can land before PR 2.
• Tests extend the existing helpers in test/test_csr.py (assert_frozen_equal, assert_terms_equal, assert_cells_equal, assert_contracted_equal) with a "backing kept" assertion instead of adding a new fixture.
• Out of scope for now: DataArray sign in to_constraint and non-scalar fillna still densify. groupby-sum pads every group to the largest group size #745 is likely covered by sparse groupby; re-check and close or narrow.

Phase 2: performance and remaining densification

Note

The following section was generated by AI.

Phase 1 made the sparse path complete and observable. A PyPSA-like benchmark on f665a260 shows where the time and memory now go. It has 200 buses x 720 snapshots, 2.09M variables, 1.73M constraints and 4.2M nonzeros.

phase	sparse	dense	dense, frozen
supply + storage + flow @ incidence	69 ms / 137 MB	2225 ms / 2965 MB	743 ms / 2965 MB
add_constraints(balance)	4 ms	33 ms	602 ms / 1739 MB
m.matrices	106 ms / 254 MB	1337 ms / 1832 MB	110 ms / 254 MB
to_highspy + to_file (LP)	1179 ms	3040 ms	1094 ms
total build + export	1.9-2.2 s	7.4 s	3.3 s

Findings:

• The sparse gain comes from three places: groupby with uneven group sizes, @, and merges of their results. From add_constraints on, a sparse model and a frozen dense model run the same code.
• Export now takes about 65% of the sparse build time. A large share of it is spent inside HiGHS.
• These operations still densify a CSR-backed expression, silently by default: shift, roll, diff, merge(dim=...), multiplication by a DataArray with a new dimension, rolling, cumsum and quadratic products. For an uneven groupby result, this costs about 250 ms and 1.4 GB.
• A CSR backing for expressions built straight from variables gives no gain. They have one term per cell and no padding, and a chain of CSR operations was not faster than dense.

Plan: four pull requests, then the docs.

PR	Issue	Content	Depends on
5: export path	#1008	Cache Model.matrices when all constraints are frozen, skip scaling lookups without scaling, skip eliminate_zeros for frozen blocks. Add the benchmark under benchmark/.	–
6: fast freeze	#1009	Direct rectangle-to-CSR in CSRConstraint.from_dense, sort only when nterm > 1	5
7: CSR kernel	#1010	Merge all operands in one COO pass, row gather for pure reorderings, no int64 intermediate in aggregated	–
8: densification gaps	#1011, #969	Sparse shift/roll/diff and merge(dim=...) (see the SnapshotWindow.merge comment below), warn_on_densify on by default in a sparse model	7
9: docs	#971	Part H, after PR 8, when the list of operations that keep the backing is final	8

Ordering notes:

• PRs 5 and 7 are independent. PRs 5 and 6 help dense models too.
• Rejected: CSR backing for variables and one-term expressions, G @ csr for sum (same speed as COO), and a lazy expression graph (it does not reduce the export costs, and it adds a second set of semantics).
• Out of scope: multiplication by a DataArray with a new dimension, rolling, cumsum and quadratic products. Revisit them if users report them.

Full benchmark output (sparse)

=== sparse  N_BUS=200 N_T=720 N_GEN=2000 N_LINE=300  nvars=2088000 ncons=1728000 nnz=4175200
vars                                58 ms   peak     44.2 MB
expr: groupby nodal supply          48 ms   peak     88.7 MB
expr: sto groupby                   66 ms   peak     19.3 MB
expr: flow @ incidence              43 ms   peak     35.3 MB
expr: balance sum                   69 ms   peak    137.0 MB
cons: balance                        4 ms   peak      4.6 MB
expr: soc shift                    128 ms   peak     24.1 MB
cons: soc                           48 ms   peak     23.6 MB
expr+cons: p <= pmax               120 ms   peak    152.7 MB
objective                           71 ms   peak     65.2 MB
matrices                           106 ms   peak    254.0 MB
to_highspy                         436 ms   peak    215.9 MB
to_file(lp)                        743 ms   peak    187.1 MB
TOTAL                             1940 ms

Timings are from one machine with tracemalloc on, so they are higher than without it. The memory peaks are stable across runs.

Children

• repr() silently and permanently densifies a CSR-backed LinearExpression #962
• Make the sparse path observable and controllable (.is_sparse, densify warning, per-call sparse=) #969
• Keep CSR backing for elementwise ops with non-scalar constants #965
• Sparse sum(dim=...) and chained groupby on CSR-backed expressions #964
• Keep CSR backing through where/mask/sel/isel/loc #966
• Widen the CSRConstraint fast path in Model.add_constraints #970
• CSRConstraint drops auxiliary coordinates when a sparse expression becomes a constraint #941
• CSRConstraint lacks loc/update/from_rule and raises bare AttributeError #963
• Emit flat/to_polars directly from the CSR store #968
• Sparse path for the objective #967
• Native sparse soften for frozen constraints (penalty= with freeze=True) #975 (PR feat(csr): native soften for frozen constraints #979)
• One Model-level switch for the sparse path instead of three independent opt-ins #976
• Document the sparse expression path (CSRLinearExpression, sparse=True, sparse_groupby) #971
• Cache Model.matrices and trim the matrix export path #1008
• Faster freeze of dense constraints into CSRConstraint #1009
• Faster CSR kernel: merge in one pass, row gather for reorderings #1010
• Keep the CSR backing through shift, roll, diff and merge(dim=...) #1011

[MaykThewessen]

Note

The following experiment results and reproduction were generated with AI assistance.

BirdFlow validation of the sparse metadata fix merged in #1003, plus one remaining fallback observed in the workload:

• has_terms and auxiliary-coordinate drop_vars caused three densifications in a 168-hour PyPSA core-model build. After the fix, one remains at SnapshotWindow.merge.
• Three fresh-process repetitions per mode: median whole-process peak RSS was 1.959 GiB with pinned linopy 0.9.1 dense and 1.740 GiB with the fixed sparse branch; matrix export was 0.120 s versus 0.030 s. This includes imports, network loading, build/export and diagnostics, not solver peak memory. Build timings were close. Source version and arithmetic mode also differ, so these numbers do not isolate the metadata fix itself.
• A separate production-chunk experiment solved the same frozen 1008-hour network as six independent weeks in each mode. All twelve solves were optimal. Dense legacy and sparse v1 gave exactly identical objective, generator dispatch, Line/Link/Transformer flows, storage SOC and nodal prices; every compared output was finite. NTC runtime state and operative-limit sidecars were restored in both modes.

Remaining fallback: SnapshotWindow.merge calls linopy.merge(..., dim="snapshot", join="outer"). _try_csr_merge currently supports term-wise merging, so coordinate concatenation densifies both operands and result. The reproduction below matches dense numerically. A sparse implementation must stack rows and preserve input order. Replacing concatenation with term-wise addition would change the operation. Overlapping labels, absent cells, constants and auxiliary coordinates need explicit coverage.

Standalone fallback reproduction and solve provenance

import warnings
import linopy
import pandas as pd
from linopy import Model, LinearExpression
from linopy.csr import CSRLinearExpression
from linopy.testing import assert_linequal

linopy.options["semantics"] = "v1"
linopy.options["warn_on_densify"] = True
m = Model(sparse=True)
x = m.add_variables(coords=[pd.Index([0, 1, 2], name="snapshot")], name="x")
parts = [(2 * x).sel(snapshot=[0, 1]), (3 * x).sel(snapshot=[2])]
expected = linopy.merge(parts, dim="snapshot", join="outer")
sparse = [LinearExpression._from_csr(CSRLinearExpression.from_dense(e.data, m), m) for e in parts]
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    result = linopy.merge(sparse, dim="snapshot", join="outer")
assert_linequal(result, expected)
print([str(w.message) for w in caught])
print(result._csr is not None, [e._csr is not None for e in sparse])
# One coordinate-merge densification notice; False, [False, False].

The private conversion only constructs this reproduction's sparse operands.

Solve provenance: linopy fa770b7cea1ffc89b0f01952775711e5230e20b9 for both modes, installed PyPSA 1.3.0, HiGHS 1.15.1 (7e5e6228c84d7778bd632520de39899ab8f18b07), HiPO with crossover, 18 threads, blueprint reuse disabled. Frozen BirdFlow source: 324e2aa26f0f4f0a797716522c8df6493cfd3b19. Source and input hashes were checked; solved networks and scorecards retained locally.

Physical balance/cycle/thermal checks passed identically. The existing scenario health gate failed identically: 23.020316 GWh localized slack versus a 1 GWh threshold. This is numerical equivalence evidence, not healthy-scenario or full-year certification. Each week has independent storage boundaries. Production dependency pins were unchanged by this experiment.

Standalone SVG. The unchanged slack is attributable to thin-line limits in the local grid model and is unrelated to this linopy issue.

Note: AI-assisted (written with model identity unavailable).