feat(csr): keep CSR backing through metadata, constant arithmetic, sum/groupby and selection

doc/release_notes.rst

@@ -44,6 +44,7 @@ Upcoming Version
 *Internal*
 
 * The sparse backing of a ``LinearExpression`` moved from ``linopy.sparse_expression`` to ``linopy.csr`` and the class ``CSRExpression`` was renamed to ``CSRLinearExpression``. Dense/sparse conversion is now spelled the same way on both CSR types: ``CSRLinearExpression.from_dense`` / ``.to_dense`` and ``CSRConstraint.from_dense`` (previously ``CSRConstraint.from_mutable``) / ``.to_dense``. ``Constraint.freeze()`` and ``CSRConstraint.mutable()`` are unchanged.
+* The first constructor argument of ``LinearExpressionGroupby`` is now named ``obj`` (previously ``data``), as it holds a ``Dataset`` or an expression; ``LinearExpressionGroupby(data=...)`` must be spelled ``LinearExpressionGroupby(obj=...)``.
 
 *Documentation*
 
@@ -54,6 +55,12 @@ Upcoming Version
 * ``@``/``dot`` against a constant matrix that holds zeros no longer densifies the result to one term per contracted member. The zero-coefficient terms are dropped, so the term dimension shrinks to the widest non-zero cell. On PyPSA's Kirchhoff Voltage Law constraint (a cycle matrix with ~3 branches per cycle) this cuts the expression from 852 to 3 terms — 284x fewer cells — which in turn shrinks the downstream ``merge``. A constant without zeros is unaffected. (`#748 <https://github.com/PyPSA/linopy/issues/748>`__)
 * ``densify_terms`` (used by ``sum(drop_zeros=True)`` and the sparse ``@`` path) is now fully vectorised. It previously counted the non-zero positions with a Python loop that scaled quadratically in the number of non-zero terms — 127 s for a (2000 x 60) expression, now 3 ms — and allocated the compacted output at the full original term width. It now allocates only the compacted width and returns the expression unchanged when it holds no zeros.
 * Under v1, ``@``/``dot`` against a constant now runs as sparse linear algebra instead of building the dense broadcast intermediate (``self * other`` then ``.sum()``): peak memory scales with ``nnz(C) x nterm`` rather than with the full broadcast shape. A ``CSRLinearExpression``-backed operand (from ``groupby(...).sum(sparse=True)``) stays CSR-backed through ``@``, and the result is the compact canonical form (duplicate variables summed, terms label-ordered, explicit zeros pruned — cell activeness is carried by ``const`` alone). (`#748 <https://github.com/PyPSA/linopy/issues/748>`__, `#756 <https://github.com/PyPSA/linopy/issues/756>`__, `#925 <https://github.com/PyPSA/linopy/issues/925>`__)
+* Reading metadata of a CSR-backed ``LinearExpression`` (``repr()``, ``shape``, ``sizes``, ``dims``, ``coords``, ``indexes``, ``isnull()``) no longer converts it to dense; it is served from the sparse backing, auxiliary coordinates included. (`#962 <https://github.com/PyPSA/linopy/issues/962>`__)
+* New ``LinearExpression.is_sparse`` tells whether an expression is CSR-backed, and its repr header reads ``LinearExpression (sparse)``. The opt-in option ``linopy.options["warn_on_densify"]`` (default ``False``) emits a ``PerformanceWarning`` naming the reason whenever a sparse backing is dropped, e.g. on ``.data``, an operation without a sparse path, or ``CSRConstraint.mutable()``. (`#969 <https://github.com/PyPSA/linopy/issues/969>`__)
+* Multiplying, dividing, adding or subtracting a non-scalar constant (``DataArray``, ``Series``, ``ndarray``, …) keeps a CSR-backed ``LinearExpression`` sparse: the coefficients are scaled per cell and only the constant changes, instead of expanding the dense rectangle. The constant is aligned exactly as on the dense path (NaN, label mismatch, joins and auxiliary coordinates); an operand that adds dimensions still densifies. ``expr.const`` is now served from the sparse backing as well. (`#965 <https://github.com/PyPSA/linopy/issues/965>`__)
+* ``sum(dim=...)`` and a further ``groupby(...).sum()`` keep a CSR-backed ``LinearExpression`` sparse: summed rows are merged in the sparse backing instead of stacking the summed dimensions into the term dimension of the dense rectangle. Auxiliary coordinates, absent cells and the constant follow the dense path; ``groupby(...).sum(sparse=False)`` still densifies. (`#964 <https://github.com/PyPSA/linopy/issues/964>`__)
+* ``where``, ``sel``, ``isel``, ``loc`` and ``[]`` keep a CSR-backed ``LinearExpression`` sparse: the selection is evaluated on the grid's row numbers with the dense (xarray) semantics, including ``drop=``, ``method=``, scalar coordinates left by a scalar selection and auxiliary coordinates, and the selected rows are gathered from the sparse backing; masked cells become absent. A condition or indexer that introduces dimensions still densifies. A merge of differing grids now stays sparse when operands carry auxiliary coordinates, which are checked for conflicts as on the dense path. (`#966 <https://github.com/PyPSA/linopy/issues/966>`__)
+* ``Model.add_constraints(..., mask=..., freeze=True)`` with a CSR-backed left-hand side applies the mask sparsely and returns a ``CSRConstraint`` instead of expanding the dense rectangle. (`#970 <https://github.com/PyPSA/linopy/issues/970>`__)
 * Persistent snapshots of tz-aware ``DatetimeIndex`` coordinates no longer materialise an object array of ``Timestamp`` per container per capture and diff. Coordinates are stored as UTC-ns arrays with the timezone identity carried alongside, making snapshot capture ~24x and warm-start diffs ~33x faster on tz-aware models, while naive and tz-aware coordinates — and differing timezones — stay correctly unequal. (`#960 <https://github.com/PyPSA/linopy/pull/960>`__)
 
 **Bug fixes**
@@ -67,6 +74,7 @@ Upcoming Version
 * ``linopy.read_netcdf`` now restores variables, expressions and constraints in their original insertion order instead of alphabetically, so ``matrices.A`` of a round-tripped model is no longer a row/column permutation of the original. The order is stored in the file; files written by earlier versions still load in sorted order. ``linopy.testing.assert_model_equal`` now also compares container order. (`#934 <https://github.com/PyPSA/linopy/issues/934>`__)
 * Adding or removing variables after a constraint was added with ``freeze=True`` no longer breaks ``model.matrices``. The frozen constraint stored dense variable positions as its matrix columns, so blocks frozen at different times disagreed on their width and stacking them raised ``ValueError: inconsistent shapes``. Raw variable labels are stored instead and mapped to positions when the matrix is assembled. Frozen constraints in netCDF files written by earlier versions are read as before. (`#926 <https://github.com/PyPSA/linopy/issues/926>`__)
 
+* ``linopy.merge`` with ``join="outer"``, ``"left"`` or ``"right"`` no longer raises an xarray ``MergeError`` when only one operand carries an auxiliary coordinate on a joined dimension. The constant was reindexed with a fill value of ``0`` that was also written into the auxiliary coordinate, so it no longer matched the coefficients' copy. The coordinate now follows its rows onto the joined grid and is NaN where no operand provides it, under both semantics. The same holds for ``.add`` / ``.sub`` / ``.mul`` / ``.div`` with a constant and a reindexing ``join=``.
 * A frozen constraint caches the label-to-position mapping of its matrix columns by weak reference and only rebuilds it when the constraint or the set of variables changes. While a persistent snapshot holds the arrays, repeated matrix assembly on an unchanged model returns the same objects, so the snapshot diff can again skip the comparison of untouched frozen constraints by object identity; one-off exports such as ``to_file`` retain no extra memory. (`#933 <https://github.com/PyPSA/linopy/issues/933>`__)
 * ``Solver.close()`` no longer leaves dangling native handles behind. The solver model is now dropped before the environment that owns it, instead of after. And the COPT and MindOpt file interfaces no longer hand back a model they already disposed: after a file-based COPT or MindOpt solve, ``model.solver_model`` is ``None`` rather than a handle into freed memory. (`#899 <https://github.com/PyPSA/linopy/pull/899>`__)
 

linopy/alignment.py

@@ -77,6 +77,7 @@ def as_constant(other: Any) -> Any:
 
 
 if TYPE_CHECKING:
+    from linopy.csr import AuxCoords
     from linopy.expressions import LinearExpression, QuadraticExpression
     from linopy.variables import Variable
 
@@ -696,7 +697,7 @@ def _matmul_operand_to_matrix(
     contracted: Sequence[str],
     new_dims: Sequence[str],
     indexes: Mapping[str, pd.Index],
-    aux_coords: Mapping[str, tuple[str, np.ndarray]],
+    aux_coords: AuxCoords,
 ) -> scipy.sparse.csr_array:
     """
     Flatten a ``@`` constant to the ``(contracted, new)`` matrix of the sparse

linopy/config.py

@@ -85,4 +85,5 @@ def __repr__(self) -> str:
     display_max_terms=6,
     semantics=LEGACY_SEMANTICS,
     sparse_groupby=False,
+    warn_on_densify=False,
 )

linopy/constraints.py

@@ -79,7 +79,7 @@
     PerformanceWarning,
     SIGNS_pretty,
 )
-from linopy.csr import Grid, csr_nterm, csr_to_term_arrays
+from linopy.csr import Grid, _densify_notice, csr_nterm, csr_to_term_arrays
 from linopy.scaling import ensure_scaling, validate_scaling
 from linopy.semantics import check_user_nan
 from linopy.types import (
@@ -1262,6 +1262,7 @@ def freeze(self) -> CSRConstraint:
 
     def to_dense(self) -> Constraint:
         """Convert to a Constraint."""
+        _densify_notice("frozen constraint converted by `to_dense()`/`mutable()`")
         return Constraint(self.data, self._model, self._name)
 
     def mutable(self) -> Constraint:
@@ -1419,25 +1420,26 @@ def from_csr(
             active,
             rhs_flat[active],
             sign,
-            grid=expr.grid,
+            grid=Grid(expr.grid.indexes),
             model=expr.model,
         )
 
 
-def csr_rhs(expr: CSRLinearExpression, rhs: Any) -> DataArray | None:
+def csr_rhs(expr: CSRLinearExpression, rhs: Any) -> DataArray | str:
     """
-    Return ``rhs`` as a DataArray on the expression grid, or None if the sparse
-    path cannot take it: a non-constant rhs, one that is no DataArray-like, or
-    one with helper dims or dims outside the grid falls back to the dense path.
+    Return ``rhs`` as a DataArray on the expression grid, or the reason the
+    sparse path cannot take it: a non-constant rhs, one that is no
+    DataArray-like, or one with helper dims or dims outside the grid falls
+    back to the dense path.
     """
     if not is_constant(rhs):
-        return None
+        return "constraint with a non-constant rhs"
     try:
         da = as_dataarray(rhs)
     except (TypeError, ValueError):
-        return None
+        return "constraint with an rhs that is not array-like"
     if set(da.dims) & set(HELPER_DIMS) or not set(da.dims) <= set(expr.grid.dims):
-        return None
+        return "constraint with an rhs over dimensions outside the grid"
     return da