Model.store_cls construction seam, Store.to_sparse(reason) for mixed merges (#1045)

doc/api.rst

@@ -691,10 +691,11 @@ Utilities
 ``options`` holds these keys:
 
 - ``display_max_rows``, ``display_max_terms``: size of the printed repr.
-- ``warn_on_densify``: emit a ``PerformanceWarning`` when a sparse (CSR)
-  backing is converted to dense. ``None`` (default) warns in a sparse model
-  on implicit conversions only, ``True`` always warns, ``False`` never
-  warns. See :doc:`sparse-models`.
+- ``warn_on_densify``: the conversion notice, a ``PerformanceWarning`` when
+  a sparse (CSR) backing is converted to dense or a dense operand is
+  converted to sparse for a merge with a sparse one. ``None`` (default)
+  warns in a sparse model on implicit conversions only, ``True`` always
+  warns, ``False`` never warns. See :doc:`sparse-models`.
 
 
 Warnings

doc/release_notes.rst

@@ -32,8 +32,13 @@ Version 1.0.0
   * ``Variables.print_labels``, ``Constraints.print_labels``, ``Model.print_infeasibilities`` and ``Model.compute_set_of_infeasible_constraints`` are removed; use ``format_labels``, ``Model.compute_infeasibilities`` and ``Model.format_infeasibilities``.
   * The unused ``Variables.dataset_names``, ``Constraints.dataset_names`` and ``LinearExpression.sanitize`` are removed.
 
-* ``linopy.csr.CSRLinearExpression`` is renamed to ``CSRStore`` without an alias. A ``LinearExpression`` now always holds a store from the new ``linopy.store`` module, a ``DenseStore`` or a ``CSRStore``, in place of the private ``_csr`` attribute. ``CSRStore.to_dense()`` returns a ``DenseStore`` instead of a ``LinearExpression``. Writing a model with a sparse objective to netcdf or copying the model no longer converts the objective to dense, and emits no densify notice. On both backends, a cell whose constant becomes NaN, e.g. by ``where``, ``reindex``, a join with ``fill_value=linopy.ABSENT`` or a constant operation yielding NaN such as ``0 / 0``, drops its terms, so its variables no longer reach ``sum()``, a constraint or the objective. The ``const`` setter and ``assign(const=...)`` of ``LinearExpression`` and ``QuadraticExpression`` raise ``ValueError`` when the assigned value holds a NaN or does not cover every coordinate of the expression. (`#1020 <https://github.com/PyPSA/linopy/issues/1020>`__)
+* ``linopy.csr.CSRLinearExpression`` is renamed to ``CSRStore`` without an alias. A ``LinearExpression`` now always holds a store from the new ``linopy.store`` module, a ``DenseStore`` or a ``CSRStore``, in place of the private ``_csr`` attribute. ``CSRStore.to_dense()`` returns a ``DenseStore`` instead of a ``LinearExpression``. Writing a model with a sparse objective to netcdf or copying the model no longer converts the objective to dense, and emits no densify notice. On both backends, a cell whose constant becomes NaN, e.g. by ``where``, ``reindex`` or a join with ``fill_value=linopy.ABSENT``, drops its terms, so its variables no longer reach ``sum()``, a constraint or the objective. The ``const`` setter and ``assign(const=...)`` of ``LinearExpression`` and ``QuadraticExpression`` raise ``ValueError`` when the assigned value holds a NaN or does not cover every coordinate of the expression. (`#1020 <https://github.com/PyPSA/linopy/issues/1020>`__)
 * An operation without a sparse path on a CSR-backed ``LinearExpression``, including reading ``.data``, now works on a dense copy and no longer converts the expression itself, which stays sparse; only assigning ``vars``, ``coeffs`` or an off-grid ``const`` replaces its store. The densify notice names the operation and the case as ``op(case)``, e.g. ``Sparse (CSR) backing densified: reindex(method=...)``, so filters matching the former reason texts need updating. A named sparse expression stays sparse through ``Model.add_expressions`` and ``Model.copy`` and keeps its name through assignments. ``densify_terms``, ``reset_const`` (so also ``Constraint.update(rhs=...)``) and ``Model.assign_coords`` keep a sparse expression sparse. ``LinearExpression.where`` and ``QuadraticExpression.where`` take ``drop`` as their only keyword argument instead of forwarding ``**kwargs`` to ``xarray.Dataset.where``, which accepts no other. (`#1020 <https://github.com/PyPSA/linopy/issues/1020>`__)
+* ``Model.store_cls`` names the store class a model builds expressions on, ``CSRStore`` in a sparse model and ``DenseStore`` otherwise. ``Variable.to_linexpr`` always returns an expression on it; the private ``_dense`` keyword is removed. (`#1045 <https://github.com/PyPSA/linopy/issues/1045>`__)
+* A dense operand in a merge with a CSR-backed one is converted up through the new ``Store.to_sparse(reason)``; a merge that stays sparse emits one ``Dense backing sparsified: merge with a sparse operand`` notice, under the same ``warn_on_densify`` option as the densify notice. A dense operand with coordinates the sparse grid cannot hold, e.g. two-dimensional auxiliary ones, now falls the whole merge back to dense with the reason ``merge(coordinates the Grid cannot hold: <names>)``. A merge of dense operands stays dense, also in a sparse model. (`#1045 <https://github.com/PyPSA/linopy/issues/1045>`__)
+* A named linear expression of a sparse model reads back from netcdf as a CSR-backed expression, unless it carries coordinates the sparse grid cannot hold, e.g. two-dimensional auxiliary ones; then it reads back dense with those coordinates. (`#1045 <https://github.com/PyPSA/linopy/issues/1045>`__)
+* A constant operation that produces NaN in a live cell, e.g. ``0 * inf``, ``0 / 0`` or ``inf - inf`` from multiplying, dividing, adding or subtracting a constant, now raises ``ValueError`` on both backends. Before, the cell kept its terms under a NaN constant. linopy trusts NaN only as absence, so mark cells to drop with ``.where(...)`` or ``mask=`` instead. Cells that were absent before stay absent. A sparse division by zero no longer emits numpy ``RuntimeWarning``\ s. (`#1045 <https://github.com/PyPSA/linopy/issues/1045>`__)
+* Cells with a NaN constant now drop their terms when a ``LinearExpression`` or ``QuadraticExpression`` is built from a dataset or its ``vars`` or ``coeffs`` are assigned, so neither a dataset with terms under a NaN constant nor a ``vars`` assignment can revive an absent cell. (`#1045 <https://github.com/PyPSA/linopy/issues/1045>`__)
 
 **Deprecations**
 

examples/sparse-models.ipynb

@@ -242,6 +242,12 @@
     "Explicit conversions are `CSRConstraint.mutable()`, `CSRConstraint.to_dense()`\n",
     "and `add_constraints(..., freeze=False)`. They stay quiet under `None`.\n",
     "\n",
+    "The option also governs the opposite conversion: a dense operand merged with\n",
+    "a sparse one, e.g. `sparse_expr + dense_expr`, is converted to sparse with one\n",
+    "`Dense backing sparsified: merge with a sparse operand` notice per merge. A\n",
+    "dense operand with coordinates the sparse grid cannot hold, such as\n",
+    "two-dimensional auxiliary ones, makes the whole merge dense instead.\n",
+    "\n",
     "To find the line that triggers a warning, turn it into an error with\n",
     "`warnings.simplefilter(\"error\", linopy.PerformanceWarning)`; the traceback\n",
     "then points to it.\n",

linopy/csr.py

@@ -53,7 +53,7 @@
 from linopy.common import coords_from_dataset, coords_to_dataset_vars
 from linopy.config import options
 from linopy.constants import HELPER_DIMS, TERM_DIM, PerformanceWarning
-from linopy.semantics import warn_outside_linopy
+from linopy.semantics import check_arithmetic_nan, warn_outside_linopy
 from linopy.store import DenseStore, Reason, Store, merged_grid
 
 if TYPE_CHECKING:
@@ -576,15 +576,27 @@ def scale(
         """
         Apply ``op`` with ``factor`` -- a scalar or one value per cell -- to
         every coefficient and to the constant. Explicit zeros are kept, as by
-        :meth:`add`; a cell made absent drops its terms.
+        :meth:`add`; a NaN it produces in a live cell raises.
         """
-        factor = np.broadcast_to(factor, (self.n_cells,))
-        per_term = np.repeat(factor, np.diff(self.csr.indptr))
+        with np.errstate(divide="ignore", invalid="ignore"):
+            const = op(self.const, factor)
+        check_arithmetic_nan(const, op, factor, self.const)
+        return self._scaled_terms(factor, op).assign_const(const)
+
+    def _scaled_terms(
+        self, factor: float | np.ndarray, op: Callable[[Any, Any], Any]
+    ) -> CSRStore:
+        """Apply ``op`` with a scalar or per-cell ``factor`` to every coefficient."""
+        per_term = factor
+        if np.ndim(factor):
+            per_term = np.repeat(factor, np.diff(self.csr.indptr))
+        with np.errstate(divide="ignore", invalid="ignore"):
+            data = op(self.csr.data, per_term)
+        check_arithmetic_nan(data, op, per_term)
         csr = scipy.sparse.csr_array(
-            (op(self.csr.data, per_term), self.csr.indices, self.csr.indptr),
-            shape=self.csr.shape,
+            (data, self.csr.indices, self.csr.indptr), shape=self.csr.shape
         )
-        return replace(self, csr=csr, const=op(self.const, factor)).absorbed()
+        return replace(self, csr=csr)
 
     def assign_const(self, const: np.ndarray) -> CSRStore:
         """
@@ -612,7 +624,7 @@ def combine(
         store = self.reindex(grid) if reindex else self
         if factor is not None:
             factor = factor.broadcast_like(const).transpose(*dims)
-            store = store.scale(factor.to_numpy().reshape(-1), op)
+            store = store._scaled_terms(factor.to_numpy().reshape(-1), op)
         return store.assign_const(const.to_numpy().reshape(-1)).with_grid(grid)
 
     def absorbed(self) -> CSRStore:
@@ -735,6 +747,16 @@ def add(self, *others: CSRStore) -> CSRStore:
         csr = _interleaved([p.csr for p in parts], shape, self.model)
         return replace(self, csr=csr, const=const, grid=merged_grid(parts)).absorbed()
 
+    @classmethod
+    def from_matmul(
+        cls,
+        source: CSRStore,
+        matrix: scipy.sparse.csr_array,
+        contracted_dims: Iterable[str],
+        new_indexes: Iterable[pd.Index],
+    ) -> CSRStore:
+        return source.dot(matrix, contracted_dims, new_indexes)
+
     def dot(
         self,
         matrix: scipy.sparse.csr_array,
@@ -983,18 +1005,26 @@ def csr_to_term_arrays(
     return vars_, coeffs
 
 
-def _densify_notice(reason: str, model: Model, explicit: bool = False) -> None:
+def _conversion_notice(message: str, model: Model, explicit: bool) -> None:
     """
-    Emit a :class:`~linopy.constants.PerformanceWarning` naming why a sparse
-    (CSR) backing is dropped, if ``options["warn_on_densify"]`` is set, or
-    left at its default ``None``, ``model`` is sparse and the conversion is
-    an implicit fallback rather than an ``explicit`` user request.
+    Emit ``message`` as a :class:`~linopy.constants.PerformanceWarning` if
+    ``options["warn_on_densify"]`` is set, or left at its default ``None``,
+    ``model`` is sparse and the conversion is an implicit fallback rather
+    than an ``explicit`` user request.
     """
     enabled = options["warn_on_densify"]
     if model.sparse and not explicit if enabled is None else enabled:
-        warn_outside_linopy(
-            f"Sparse (CSR) backing densified: {reason}.", PerformanceWarning
-        )
+        warn_outside_linopy(message, PerformanceWarning)
+
+
+def _densify_notice(reason: str, model: Model, explicit: bool = False) -> None:
+    """The conversion notice naming why a sparse (CSR) backing is dropped."""
+    _conversion_notice(f"Sparse (CSR) backing densified: {reason}.", model, explicit)
+
+
+def _sparsify_notice(reason: str, model: Model) -> None:
+    """The conversion notice naming why a dense backing is converted to sparse."""
+    _conversion_notice(f"Dense backing sparsified: {reason}.", model, False)
 
 
 def _readonly(values: np.ndarray, copy: bool = True) -> np.ndarray:

linopy/expressions.py

@@ -111,6 +111,7 @@
     FillValueLike,
     _shared_dim_mismatch_message,
     absorb_absence,
+    check_arithmetic_nan,
     check_join_fill_value,
     check_user_nan,
     enforce_aux_conflict,
@@ -556,11 +557,11 @@ def sum(
         multikey_frame = None if use_fallback else _multikey_value_frame(group, labels)
 
         store = self._store
-        sparse = self.model.sparse or isinstance(store, CSRStore)
+        store_cls = result_store_cls(store, self.model)
         if multikey_frame is not None and not observed:
             _warn_dense_grid(multikey_frame)
 
-        if sparse:
+        if store_cls is CSRStore:
             grouper = multikey_frame
             if grouper is None:
                 is_1d = isinstance(group, DataArray) and group.ndim == 1
@@ -578,7 +579,7 @@ def sum(
                     if isinstance(grouper, pd.DataFrame)
                     else str(grouper.name or "group")
                 )
-                res = CSRStore.from_grouper(
+                res = store_cls.from_grouper(
                     store, self.model, grouper, group_name, stacked, coord_dims
                 )
                 return LinearExpression._from_store(res, self.model)
@@ -1100,10 +1101,11 @@ def _combined_with_constant(
         """
         expr = self._reindexed_to(self_const, needs_data_reindex)
         fields = {"const": op(self_const, operand)}
+        check_arithmetic_nan(fields["const"], op, operand, self_const)
         if scale:
             fields["coeffs"] = op(expr.coeffs, operand)
-        data = absorb_absence(assign_fields(expr.data, **fields))
-        return self.__class__(data, self.model)
+            check_arithmetic_nan(fields["coeffs"], op, operand, labels=expr.vars)
+        return self.__class__(assign_fields(expr.data, **fields), self.model)
 
     def _apply_constant_op(
         self,
@@ -2118,7 +2120,7 @@ def _from_store(cls, store: Store, model: Model) -> Self:
         return obj
 
     def _set_data(self, data: Dataset) -> None:
-        self._store = DenseStore(data, self._model)
+        self._store = DenseStore(absorb_absence(data), self._model)
 
     def _dense_store(self, reason: Reason | None) -> DenseStore:
         """
@@ -2279,6 +2281,7 @@ def _combined_with_constant(
         scale: bool,
     ) -> Self:
         const = op(self_const, operand)
+        check_arithmetic_nan(const, op, operand, self_const)
         factor = operand if scale else None
         store = self._store.combine(const, op, factor, needs_data_reindex)
         if isinstance(store, str):
@@ -2571,7 +2574,7 @@ def __mul__(
 
         try:
             if isinstance(other, variables.Variable):
-                other = other.to_linexpr(_dense=True)
+                other = other._linexpr(DenseStore)
             elif isinstance(other, variables.ScalarVariable):
                 other = other.to_linexpr()
 
@@ -2647,9 +2650,10 @@ def _sparse_matmul(self, other: ConstantLike) -> LinearExpression | Reason:
         Contract against a constant as one sparse matrix product, skipping the
         dense broadcast intermediate of ``(self * other).sum(dim)``.
 
-        Returns the reason where the dense path owns the semantics: non-unique
-        grid labels, a zero-size grid, an operand sharing no dimension with
-        the grid, and an unlabelled output dimension. The result is the compact
+        Returns the reason where the dense path owns the semantics: a dense
+        expression with coordinates the grid cannot hold, non-unique grid
+        labels, a zero-size grid, an operand sharing no dimension with the
+        grid, and an unlabelled output dimension. The result is the compact
         canonical form of :meth:`CSRStore.dot`, so its term
         count may differ from the dense path's while the values agree; the
         result is CSR-backed when the input was or the model is sparse.
@@ -2659,6 +2663,8 @@ def _sparse_matmul(self, other: ConstantLike) -> LinearExpression | Reason:
         if not self.coord_dims:
             return "matmul(zero-dimensional expression)"
         csr = self._store.to_sparse()
+        if isinstance(csr, str):
+            return f"matmul({csr})"
         if not csr.grid.is_unique:
             return "matmul(non-unique labels)"
         if csr.grid.size == 0:
@@ -2673,10 +2679,9 @@ def _sparse_matmul(self, other: ConstantLike) -> LinearExpression | Reason:
             da, contracted, new_dims, csr.grid.indexes, csr.grid.aux
         )
         new_indexes = [da.indexes[d].rename(d) for d in new_dims]
-        res = csr.dot(matrix, contracted, new_indexes)
-        if self.is_sparse or self.model.sparse:
-            return type(self)._from_store(res, self.model)
-        return type(self)._from_store(res.to_dense(), self.model)
+        store_cls = result_store_cls(self._store, self.model)
+        res = store_cls.from_matmul(csr, matrix, contracted, new_indexes)
+        return type(self)._from_store(res, self.model)
 
     @property
     def flat(self) -> pd.DataFrame:
@@ -3157,7 +3162,7 @@ def __init__(self, data: Dataset | None, model: Model) -> None:
         self._data = data
 
     def _set_data(self, data: Dataset) -> None:
-        self._data = data
+        self._data = absorb_absence(data)
 
     @property
     def data(self) -> Dataset:
@@ -3407,6 +3412,11 @@ def as_expression(
         return LinearExpression(obj, model)
 
 
+def result_store_cls(store: Store, model: Model) -> type[Store]:
+    """The store class an operation on ``store`` builds its result on."""
+    return CSRStore if isinstance(store, CSRStore) else model.store_cls
+
+
 def _dense_copies(exprs: Sequence[Any], reason: Reason) -> list[Any]:
     """
     ``exprs`` with every CSR-backed LinearExpression replaced by a dense

linopy/io.py

@@ -31,8 +31,10 @@
     to_polars,
 )
 from linopy.constants import FACTOR_DIM, SOS_DIM_ATTR, SOS_TYPE_ATTR
+from linopy.csr import CSRStore
 from linopy.objective import Objective, linear_part
 from linopy.scaling import constraint_scaling_lookup, variable_scaling_lookup
+from linopy.store import DenseStore, expression_dataset
 
 if TYPE_CHECKING:
     from cuopt.linear_programming import DataModel as cuoptDataModel
@@ -1270,7 +1272,11 @@ def owner(k: str) -> str | None:
         ):
             expr = QuadraticExpression(expr_ds, m)
         else:
-            expr = LinearExpression(expr_ds, m)
+            dense = DenseStore(expression_dataset(expr_ds, m), m)
+            store = dense.to_sparse() if m.store_cls is CSRStore else dense
+            expr = LinearExpression._from_store(
+                dense if isinstance(store, str) else store, m
+            )
         expr.attrs["name"] = name
         expressions[name] = expr
 

linopy/model.py

@@ -60,7 +60,7 @@
     Constraints,
     CSRConstraint,
 )
-from linopy.csr import _densify_notice
+from linopy.csr import CSRStore, _densify_notice
 from linopy.dualization import dualize
 from linopy.expressions import (
     Expressions,
@@ -107,6 +107,7 @@
     sos_reformulation_context,
     undo_sos_reformulation,
 )
+from linopy.store import DenseStore, Store
 from linopy.types import (
     ConstantLike,
     ConstraintLike,
@@ -521,6 +522,18 @@ def sparse(self) -> bool:
         """
         return self._sparse
 
+    @property
+    def store_cls(self) -> type[Store]:
+        """
+        The store class the construction seam builds on: ``to_linexpr``,
+        ``groupby(...).sum()``, ``@`` and the netcdf reader yield a
+        :class:`~linopy.csr.CSRStore` in a sparse model and a
+        :class:`~linopy.store.DenseStore` otherwise. Quadratic products,
+        ``rolling``, ``cumsum`` and the ``LinearExpression`` constructor
+        build dense.
+        """
+        return CSRStore if self._sparse else DenseStore
+
     @property
     def set_names_in_solver_io(self) -> bool:
         """Whether direct solver exports include names by default."""