feat(csr): one read-only Model(sparse=True) key for the sparse path

benchmarks/patterns/nodal_balance.py

@@ -9,7 +9,7 @@
 current (dense) kernel.
 
 ``nodal_balance_sparse`` builds the identical constraint through the
-CSR-backed path (``sum(sparse=True)`` + ``freeze=True`` under v1): the
+CSR-backed path (``Model(sparse=True)`` under v1): the
 grouped sum never materializes the padded rectangle and is realized directly
 as a CSRConstraint, so its peak memory should stay flat across the severity
 sweep — the pair makes the padding cost visible.
@@ -62,15 +62,15 @@ def _build(severity: int, sparse: bool) -> linopy.Model:
     buses = pd.RangeIndex(N_BUS, name="bus")
     rng = np.random.default_rng(1)
 
-    m = linopy.Model()
+    m = linopy.Model(sparse=sparse)
     gen = m.add_variables(lower=0, coords=[gens, time], name="gen")
 
     bus_of_gen = pd.Series(_bus_of_gen(severity), index=gens, name="bus")
-    supply = (1 * gen).groupby(bus_of_gen).sum(sparse=sparse)
+    supply = (1 * gen).groupby(bus_of_gen).sum()
     demand = xr.DataArray(
         rng.uniform(10.0, 100.0, size=(N_BUS, N_TIME)), coords=[buses, time]
     )
-    m.add_constraints(supply == demand, name="balance", freeze=sparse)
+    m.add_constraints(supply == demand, name="balance")
     m.add_objective(gen.sum())
     return m
 

doc/api.rst

@@ -411,7 +411,7 @@ CSRConstraint
 =============
 
 Memory-efficient, immutable constraint representation backed by a scipy
-CSR sparse matrix. Opt in via ``Model(freeze_constraints=True)`` or
+CSR sparse matrix. Opt in via ``Model(sparse=True)`` or
 ``Model.add_constraints(..., freeze=True)``. See the
 :doc:`creating-constraints` guide for usage.
 

doc/release_notes.rst

@@ -35,6 +35,11 @@ Upcoming Version
 * A frozen ``CSRConstraint`` is softened natively, without densifying: the slack terms are appended to its sparse rows in place, so ``add_constraints(..., penalty=..., freeze=True)`` works as well. The slack is kept by ``to_dense()``/``mutable()``, ``Constraint.freeze()`` and the netcdf round trip. (`#975 <https://github.com/PyPSA/linopy/issues/975>`__)
 * The slack variable(s) created by ``soften()`` can be retrieved afterwards via the new ``slack`` property of ``Constraint`` and ``CSRConstraint``.
 
+*Sparse model*
+
+* New read-only model key ``Model(sparse=True)`` turns on the whole sparse (CSR) path for one model: ``groupby(...).sum()`` and ``@``/``dot`` against a constant return CSR-backed expressions, and ``add_constraints`` freezes every constraint unless ``freeze=False`` is passed. It requires the v1 semantics: ``Model(sparse=True)`` raises under legacy, and so do ``groupby(...).sum()`` and ``@`` on a sparse model, and ``read_netcdf`` of one, once the semantics are switched back to legacy. It rejects ``chunk``. The key is kept by ``Model.copy`` and the netcdf round trip. ``@`` no longer reads ``linopy.options["sparse_groupby"]``. (`#976 <https://github.com/PyPSA/linopy/issues/976>`__)
+* Deprecated in favour of ``Model(sparse=True)``, each with a ``FutureWarning`` and to be removed with the legacy semantics: ``Model(freeze_constraints=...)`` and the ``Model.freeze_constraints`` setter, ``groupby(...).sum(sparse=...)`` and ``linopy.options["sparse_groupby"]``. They keep their current behaviour until then, except that ``@`` ignores ``sparse_groupby``, and netcdf files that store ``freeze_constraints`` still load. (`#976 <https://github.com/PyPSA/linopy/issues/976>`__)
+
 *Other*
 
 * New method :meth:`linopy.Model.assign_coords` reassigns coordinate values across an existing model — variables, constraints (dense and CSR-backed), expressions and parameters — without changing the model's shape: ``m.assign_coords(snapshot=new_snapshots)``. Values-only: the new values must match the length of the dimension's full-index container, and containers holding subsets of the dimension are mapped by label, preserving the subset relation. Dataset variable order is preserved. Under v1 semantics, ``Model.solve()`` raises when containers carry labels on a shared dimension that are neither equal nor subsets of one another. Typical use is advancing the window in rolling-horizon optimization with the persistent solver interface. (https://github.com/PyPSA/linopy/issues/767)

examples/creating-constraints.ipynb

@@ -294,9 +294,9 @@
    "id": "27",
    "metadata": {},
    "source": [
-    "### Freezing all constraints globally\n",
+    "### Building a sparse model\n",
     "\n",
-    "Set `freeze_constraints=True` on the `Model` to automatically freeze every constraint added via `add_constraints`:"
+    "Set `sparse=True` on the `Model` to build it on the sparse path: every constraint added via `add_constraints` is frozen, and `groupby(...).sum()` and `@` return CSR-backed expressions. Pass `freeze=False` to keep a single constraint mutable. The setting is fixed at creation and requires the v1 semantics:"
    ]
   },
   {
@@ -306,7 +306,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "m3 = Model(freeze_constraints=True)\n",
+    "m3 = Model(sparse=True)\n",
     "z = m3.add_variables(coords=[np.arange(50)], name=\"z\")\n",
     "m3.add_constraints(z >= 0, name=\"lower\")\n",
     "m3.add_constraints(z <= 100, name=\"upper\")\n",
@@ -382,7 +382,7 @@
     "Additionally, if you don't need variable and constraint names in the solver (e.g. for batch solves), you can disable name export for extra speed:\n",
     "\n",
     "```python\n",
-    "m = Model(freeze_constraints=True, set_names_in_solver_io=False)\n",
+    "m = Model(sparse=True, set_names_in_solver_io=False)\n",
     "```"
    ]
   }

linopy/config.py

@@ -12,6 +12,10 @@
 LEGACY_SEMANTICS = "legacy"
 V1_SEMANTICS = "v1"
 VALID_SEMANTICS = {LEGACY_SEMANTICS, V1_SEMANTICS}
+SPARSE_DEPRECATION = (
+    "is deprecated and will be removed with the legacy semantics; use "
+    "Model(sparse=True) instead"
+)
 
 
 class LinopySemanticsWarning(FutureWarning):
@@ -50,6 +54,13 @@ def set_value(self, **kwargs: Any) -> None:
                     f"Invalid semantics: {v!r}. "
                     f"Must be one of {sorted(VALID_SEMANTICS)}."
                 )
+            if k == "sparse_groupby" and v:
+                from linopy.semantics import warn_outside_linopy
+
+                warn_outside_linopy(
+                    f'linopy.options["sparse_groupby"] {SPARSE_DEPRECATION}.',
+                    FutureWarning,
+                )
             self._current_values[k] = v
 
     def get_value(self, name: str) -> Any:

linopy/csr.py

@@ -1,8 +1,8 @@
 """
 The sparse backing of a LinearExpression: ``A @ x + c`` in CSR form.
 
-Under v1 semantics, ``expr.groupby(g).sum(sparse=True)`` (or
-``linopy.options["sparse_groupby"]``) returns an ordinary
+In a sparse model (``Model(sparse=True)``, v1 semantics only),
+``expr.groupby(g).sum()`` returns an ordinary
 :class:`~linopy.expressions.LinearExpression` backed by a
 :class:`CSRLinearExpression`: same public type, different backing, akin to
 dask-backed xarray objects. The CSR form is canonical (duplicate variables

linopy/expressions.py

@@ -94,6 +94,7 @@
     to_polars,
 )
 from linopy.config import (
+    SPARSE_DEPRECATION,
     options,
 )
 from linopy.constants import (
@@ -127,6 +128,7 @@
     join_fill,
     reindex_like_if_needed,
     warn_legacy,
+    warn_outside_linopy,
 )
 from linopy.types import (
     CONSTANT_TYPES,
@@ -571,6 +573,7 @@ def sum(
             default one. Kept as an escape hatch. Leave at False unless the
             default misbehaves. Defaults to False.
         sparse : bool, optional
+            Deprecated, build the model with ``Model(sparse=True)`` instead.
             Build the grouped sum in CSR form behind the ordinary
             LinearExpression type — no group-size padding; a still-sparse
             lhs reaching ``Model.add_constraints`` with ``freeze=True``
@@ -579,7 +582,8 @@ def sum(
             Series/DataFrame, 1-D DataArray and coordinate-name-list groupers
             over an existing dimension; with a name list and ``observed=True``
             the CSR result stays compact over the observed key combinations.
-            Requires v1 semantics. Defaults to
+            Requires v1 semantics. Defaults to True for a sparse model, a
+            CSR-backed expression or the deprecated
             ``linopy.options["sparse_groupby"]``. See :mod:`linopy.csr`.
         observed : bool
             Only applies when grouping by a list of coordinate names. If True,
@@ -607,15 +611,22 @@ def sum(
 
         multikey_frame = None if use_fallback else _multikey_value_frame(group, labels)
 
+        self.model._check_sparse_semantics()
         csr = self._csr
         explicit_sparse = sparse is True
         if sparse is None:
-            sparse = is_v1() and (options["sparse_groupby"] or csr is not None)
-        elif sparse and not is_v1():
-            raise ValueError(
-                "sparse groupby-sum requires v1 semantics; opt in with "
-                "linopy.options['semantics'] = 'v1'."
+            sparse = is_v1() and (
+                self.model.sparse or options["sparse_groupby"] or csr is not None
+            )
+        else:
+            warn_outside_linopy(
+                f"groupby(...).sum(sparse=...) {SPARSE_DEPRECATION}.", FutureWarning
             )
+            if sparse and not is_v1():
+                raise ValueError(
+                    "sparse groupby-sum requires v1 semantics; opt in with "
+                    "linopy.options['semantics'] = 'v1'."
+                )
         if multikey_frame is not None and not observed:
             _warn_dense_grid(multikey_frame)
 
@@ -653,6 +664,8 @@ def sum(
                     "DataArray or list of coordinate names as grouper over an "
                     "existing dimension, without use_fallback."
                 )
+            if csr is None:
+                _densify_notice("groupby-sum with a grouper without a sparse path")
 
         if multikey_frame is not None:
             group = multikey_frame
@@ -1586,7 +1599,8 @@ def dot(self, other: SideLike) -> Self | QuadraticExpression:
 
         Identical to ``@``. For a :class:`LinearExpression` that includes the
         sparse contraction under v1; :class:`QuadraticExpression` always takes
-        the dense path. There is no per-call ``sparse=`` flag.
+        the dense path. The result is CSR-backed for a sparse model
+        (``Model(sparse=True)``) or a CSR-backed input.
         """
         return self.__matmul__(other)
 
@@ -2829,8 +2843,9 @@ def __matmul__(
         form -- duplicate variables summed, terms label-ordered, explicit
         zeros pruned -- so its term count may differ from the dense path's
         while the values agree. Returns a CSR-backed result when ``self`` is
-        CSR-backed.
+        CSR-backed or the model is sparse.
         """
+        self.model._check_sparse_semantics()
         other = as_constant(other)
         other_is_const = not isinstance(other, LinearExpression | variables.Variable)
         if other_is_const and is_v1() and type(self) is LinearExpression:
@@ -2857,7 +2872,7 @@ def _sparse_matmul(self, other: ConstantLike) -> LinearExpression | None:
         grid, and an unlabelled output dimension. The result is the compact
         canonical form of :meth:`CSRLinearExpression.contracted`, so its term
         count may differ from the dense path's while the values agree; the
-        result is CSR-backed when the input was.
+        result is CSR-backed when the input was or the model is sparse.
         """
         if is_nan_scalar(other):
             check_user_nan(op_kind="mul")
@@ -2890,7 +2905,7 @@ def _sparse_matmul(self, other: ConstantLike) -> LinearExpression | None:
         )
         new_indexes = [da.indexes[d].rename(d) for d in new_dims]
         res = csr.contracted(matrix, contracted, new_indexes)
-        if self._csr is not None or options["sparse_groupby"]:
+        if self._csr is not None or self.model.sparse:
             return type(self)._from_csr(res, self.model)
         return res.to_dense()
 

linopy/io.py

@@ -1099,6 +1099,7 @@ def with_prefix(ds: xr.Dataset, prefix: str) -> xr.Dataset:
     params = [with_prefix(m.parameters, "parameters")]
 
     scalars = {k: getattr(m, k) for k in m.scalar_attrs}
+    scalars |= {"sparse": m.sparse, "freeze_constraints": m._freeze_constraints}
     ds = xr.merge(vars + cons + exprs + obj + params, combine_attrs="drop_conflicts")
     ds = ds.assign_attrs(scalars)
     ds.attrs[NETCDF_VERSION_ATTR] = version("linopy")
@@ -1169,8 +1170,8 @@ def read_netcdf(path: Path | str, **kwargs: Any) -> Model:
     if isinstance(path, str):
         path = Path(path)
 
-    m = Model()
     ds = xr.load_dataset(path, **kwargs)
+    m = Model(sparse=bool(ds.attrs.get("sparse", False)))
 
     def has_prefix(k: str, prefix: str) -> bool:
         return k.rsplit("-", 1)[0] == prefix
@@ -1279,6 +1280,7 @@ def container_names(kind: str) -> list[str]:
     for k in m.scalar_attrs:
         if k in ds.attrs:
             setattr(m, k, ds.attrs[k])
+    m._freeze_constraints = bool(ds.attrs.get("freeze_constraints", False))
 
     if max(m._xCounter, m._cCounter) > np.iinfo(np.int32).max:
         m._dtypes["labels"] = np.int64
@@ -1354,9 +1356,9 @@ def copy(m: Model, include_solution: bool = False, deep: bool = True) -> Model:
         chunk=m._chunk,
         force_dim_names=m._force_dim_names,
         auto_mask=m._auto_mask,
-        freeze_constraints=m.freeze_constraints,
         set_names_in_solver_io=m.set_names_in_solver_io,
         solver_dir=str(m._solver_dir),
+        sparse=m.sparse,
     )
 
     new_model._variables = Variables(
@@ -1438,6 +1440,7 @@ def _copy_con(name: str, con: ConstraintBase) -> ConstraintBase:
     for attr in m.scalar_attrs:
         if include_solution or attr not in SOLVE_STATE_ATTRS:
             setattr(new_model, attr, getattr(m, attr))
+    new_model._freeze_constraints = m._freeze_constraints
 
     if m._sos_reformulation_state is not None:
         new_model._sos_reformulation_state = _copy.deepcopy(m._sos_reformulation_state)