fix(semantics): close gaps and false positives in legacy v1 warnings

doc/release_notes.rst

@@ -16,9 +16,11 @@ Upcoming Version
   * conflicting auxiliary coordinates raise, instead of being silently dropped.
   * a reindexing ``join=`` fills the positions it creates per side: the linopy operand contributes the zero expression, the constant operand contributes ``fill_value=`` (new on ``.add`` / ``.sub`` / ``.mul`` / ``.div``), which defaults to "does not apply here" — so a missing divisor zeroes the term instead of leaving it unscaled, and a missing numerator no longer comes out as ``1 / divisor``. Absence carried in by an operand is untouched and still propagates. Pass ``fill_value=linopy.ABSENT`` to opt out of the fill altogether — the join's coordinates are kept but every position it creates stays absent, on constant and expression operands alike (``linopy.merge`` included). (https://github.com/PyPSA/linopy/issues/890)
   * a first-class ``pd.MultiIndex`` dimension is rejected in favour of a flat dimension with the levels as auxiliary coordinates.
+  * a quadratic product (``x * y``, ``expr * expr``) over shared dimensions with a differing label set or a reorder raises, instead of being silently outer-joined by label; align the factors with ``.sel`` / ``.reindex`` / ``linopy.align`` first, as the product takes no ``join=``.
   * a ``groupby`` grouper aligns to the grouped dimension by label — a grouper whose labels differ in set or order from the dimension raises instead of being matched positionally, and a multi-key grouper yields a flat ``group`` dimension (keys as aux coords) rather than a stacked ``group`` MultiIndex. (https://github.com/PyPSA/linopy/issues/827)
 
 * Every operation whose result changes under v1 emits a ``LinopySemanticsWarning`` under legacy, naming the fix — so a model can be migrated incrementally before opting in. The full rules are specified in :doc:`the arithmetic convention <design/convention>`.
+* The legacy warnings now fire where the legacy and v1 results actually diverge, once per operation and pointing at the user's code. Arithmetic that fills an expression's absent slots with 0 (``x.diff("t") >= 0``, ``expr.shift(t=1) + y``, ``expr.where(mask) + 5``, ``expr.reindex(...) + 1``) warns, as does a constraint of a masked variable or expression against a constant (``x.where(mask) <= 3``, which legacy keeps as a labelled empty row and v1 drops), a reindexing ``join=`` that creates divisor positions on ``.div`` (filled with 1 under legacy, ``inf`` under v1) and a quadratic product over mismatched or reordered coordinates, which keeps integer labels under legacy. Reductions (``sum``, ``groupby(...).sum()``) of a masked variable no longer warn, as their results agree. The coordinate-mismatch messages describe the alignment legacy actually applies: positional when sizes match, otherwise an outer join by label for a merge; by label for a quadratic product, a constraint RHS (ignoring an explicit ``join=``, which v1 applies) and ``Variable * constant``.
 
 
 *Numerical scaling*

linopy/alignment.py

@@ -791,8 +791,7 @@ def _reindex_reordered_dims(
                 warn_legacy(
                     _legacy_const_reorder_message(
                         str(dim), expected_idx.values, actual_idx.values
-                    ),
-                    stacklevel=6,
+                    )
                 )
             arr = arr.reindex({dim: expected_idx})
     return arr

linopy/expressions.py

@@ -116,6 +116,7 @@
     FillValueLike,
     _legacy_coord_mismatch_message,
     _legacy_group_multiindex_message,
+    _legacy_join_divisor_fill_message,
     _legacy_nan_rhs_constraint_message,
     _shared_dim_mismatch_message,
     absorb_absence,
@@ -127,8 +128,10 @@
     is_nan_scalar,
     is_v1,
     join_fill,
+    legacy_warnings_muted,
     reindex_like_if_needed,
     warn_legacy,
+    warn_legacy_absence,
     warn_outside_linopy,
 )
 from linopy.types import (
@@ -1073,13 +1076,25 @@ def _multiply_by_linear_expression(
         # multiplication: (v1 + c1) * (v2 + c2) = v1 * v2 + c1 * v2 + c2 * v1 + c1 * c2
         # with v being the variables and c the constants
         # merge on factor dimension only returns v1 * v2 + c1 * c2
+        enforce_merge_dims(
+            [self.data, other.data],
+            concat_dim=FACTOR_DIM,
+            context="quadratic product",
+            alignment="by label, as an outer join",
+            joinable=False,
+        )
+        warn_legacy_absence(self.const, other.const)
         ds = other.data[["coeffs", "vars"]].sel(_term=0).broadcast_like(self.data)
         ds = assign_multiindex_safe(ds, const=other.const)
-        res = merge([self, ds], dim=FACTOR_DIM, cls=QuadraticExpression)
-        if self.has_constant:
-            res = res + other.reset_const() * self.const.fillna(0)
-        if other.has_constant:
-            res = res + self.reset_const() * other.const.fillna(0)
+        join: JoinOptions = (
+            "override" if _coord_sizes_agree([self.data, ds]) else "outer"
+        )
+        res = merge([self, ds], dim=FACTOR_DIM, cls=QuadraticExpression, join=join)
+        with legacy_warnings_muted():
+            if self.has_constant:
+                res = res + other.reset_const() * self.const.fillna(0)
+            if other.has_constant:
+                res = res + self.reset_const() * other.const.fillna(0)
         return cast(QuadraticExpression, res)
 
     def _absorb_join_absence(self, fill_value: FillValueLike) -> Self:
@@ -1132,7 +1147,7 @@ def _align_constant(
         rather than by ``xr.align(join="exact")``, whose wording is not
         API-stable across xarray releases.
         """
-        enforce_aux_conflict([self.const, other], stacklevel=4)
+        enforce_aux_conflict([self.const, other])
         other_fill = {other.name: join_fill(fill_value, 0)}
         if join is None:
             if is_v1():
@@ -1143,8 +1158,7 @@ def _align_constant(
                     warn_legacy(
                         _legacy_coord_mismatch_message(
                             "this operator's constant operand", *mismatch
-                        ),
-                        stacklevel=4,
+                        )
                     )
                 if other.sizes == self.const.sizes:
                     return self.const, other.assign_coords(coords=self.coords), False
@@ -1258,10 +1272,12 @@ def _add_constant_legacy(
         if np.isscalar(other) and join is None:
             if is_nan_scalar(other):
                 check_user_nan()
+            warn_legacy_absence(self.const)
             return self.assign(const=self.const.fillna(0) + other)
         self_const, da, needs_reindex = self._broadcast_and_align(
             other, fill_value, join
         )
+        warn_legacy_absence(self_const)
         expr = self._reindexed_to(self_const, needs_reindex)
         return expr.assign(const=self_const.fillna(0) + da.fillna(0))
 
@@ -1315,6 +1331,7 @@ def _apply_constant_op_legacy(
         self_const, factor, needs_reindex = self._broadcast_and_align(
             other, factor_fill, join, op_kind
         )
+        warn_legacy_absence(self_const)
         factor = factor.fillna(nan_fill)
         self_const = self_const.fillna(0)
         expr = self._reindexed_to(self_const, needs_reindex)
@@ -1343,8 +1360,14 @@ def _divide_by_constant(
         join: JoinOptions | None = None,
         fill_value: FillValueLike = None,
     ) -> Self:
-        if fill_value is None:
-            fill_value = np.inf if is_v1() else 1
+        if fill_value is None and is_v1():
+            fill_value = np.inf
+        elif fill_value is None:  # LEGACY: remove at 1.0
+            if join in ("left", "outer") and not _labels_within(
+                self.const, as_dataarray(other)
+            ):
+                warn_legacy(_legacy_join_divisor_fill_message())
+            fill_value = 1
         return self._apply_constant_op(
             other,
             operator.truediv,
@@ -1912,19 +1935,14 @@ def to_constraint(
                 rhs = self._broadcast_rhs(rhs)
                 if rhs.isnull().any():
                     check_user_nan()
-        else:  # LEGACY: remove at 1.0 — see doc/design/legacy-removal.rst.
-            if isinstance(rhs, CONSTANT_TYPES):
-                rhs = self._broadcast_rhs(rhs)
-                mismatch = first_mismatched_dim(self.const, rhs)
-                if mismatch is not None:
-                    warn_legacy(
-                        _legacy_coord_mismatch_message("constraint RHS", *mismatch)
-                    )
-                if bool(rhs.isnull().any()):
-                    warn_legacy(_legacy_nan_rhs_constraint_message())
-                rhs = rhs.reindex_like(self.const, fill_value=np.nan)
-            if isinstance(rhs, DataArray):
-                rhs_nan_mask = rhs.isnull()
+        elif isinstance(rhs, CONSTANT_TYPES):  # LEGACY: remove at 1.0
+            rhs = self._broadcast_rhs(rhs)
+            self._warn_legacy_rhs_alignment(rhs, join)
+            rhs_nan = rhs.isnull()
+            if bool(rhs_nan.any()):
+                warn_legacy(_legacy_nan_rhs_constraint_message())
+            rhs_nan_mask = rhs_nan.reindex_like(self.const, fill_value=True)
+            rhs = rhs.fillna(0).reindex_like(self.const, fill_value=0)
 
         all_to_lhs = self.sub(rhs, join=join).data
         computed_rhs = -all_to_lhs.const
@@ -1935,6 +1953,40 @@ def to_constraint(
         )
         return constraints.Constraint(data, model=self.model)
 
+    # LEGACY: remove at 1.0
+    def _warn_legacy_rhs_alignment(
+        self, rhs: DataArray, join: JoinOptions | None
+    ) -> None:
+        """
+        Legacy reindexes a constant RHS onto the expression by label, whatever
+        ``join``; warn unless the RHS covers every label and ``join`` would not
+        add, drop or reorder any.
+        """
+        mismatch = first_mismatched_dim(self.const, rhs)
+        if mismatch is None:
+            return
+        alignment = (
+            "by label onto the expression, no constraint where the RHS lacks a label"
+        )
+        if join is None:
+            warn_legacy(
+                _legacy_coord_mismatch_message(
+                    "constraint RHS", *mismatch, alignment=alignment
+                )
+            )
+            return
+        covered = join in ("left", "inner", "outer") and _labels_within(self.const, rhs)
+        if covered and (join != "outer" or _labels_within(rhs, self.const)):
+            return
+        warn_legacy(
+            _legacy_coord_mismatch_message(
+                "constraint RHS",
+                *mismatch,
+                alignment=f"{alignment}, ignoring `join={join!r}`",
+                v1=f"Under v1 `join={join!r}` is applied instead.",
+            )
+        )
+
     def _broadcast_rhs(self, rhs: ConstantLike) -> DataArray:
         """Broadcast a constant RHS onto the expression's coordinates."""
         rhs = broadcast_to_coords(
@@ -3429,11 +3481,11 @@ class QuadraticExpression(BaseExpression):
     _fill_value = {"vars": -1, "coeffs": np.nan, "const": np.nan}
 
     def __init__(self, data: Dataset | None, model: Model) -> None:
-        super().__init__(data, model)
-
         if data is None:
             da = xr.DataArray([[], []], dims=[FACTOR_DIM, TERM_DIM])
-            data = Dataset({"coeffs": da, "vars": da, "const": 0})
+            data = Dataset({"coeffs": da, "vars": da, "const": 0.0})
+        super().__init__(data, model)
+        data = self._data
         if FACTOR_DIM not in data.vars.dims:
             raise ValueError(f"Data does not include dimension {FACTOR_DIM}")
         elif data.sizes[FACTOR_DIM] != 2:
@@ -3709,6 +3761,18 @@ def _densify_all(exprs: Iterable[Any], reason: str) -> None:
         e._csr = None
 
 
+def _coord_sizes_agree(objs: Sequence[Any]) -> bool:
+    """Whether all ``objs`` carry the same non-helper dimension sizes."""
+    sizes = [{k: v for k, v in o.sizes.items() if k not in HELPER_DIMS} for o in objs]
+    return check_common_keys_values(sizes)
+
+
+def _labels_within(a: DataArray, b: DataArray) -> bool:
+    """Whether every label of ``a`` on a dim shared with ``b`` is also in ``b``."""
+    shared = set(a.indexes) & set(b.indexes)
+    return all(a.indexes[d].isin(b.indexes[d]).all() for d in shared)
+
+
 def _aligned(
     csrs: list[CSRLinearExpression], join: JoinOptions | None, fill: float
 ) -> list[CSRLinearExpression] | str:
@@ -3929,10 +3993,7 @@ def merge(
     if join is not None:
         override = join == "override"
     elif issubclass(cls, linopy_types) and dim in HELPER_DIMS:
-        coord_dims = [
-            {k: v for k, v in e.sizes.items() if k not in HELPER_DIMS} for e in exprs
-        ]
-        override = check_common_keys_values(coord_dims)  # type: ignore
+        override = _coord_sizes_agree(exprs)
     else:
         override = False
 
@@ -3955,6 +4016,7 @@ def merge(
 
     skipna = not is_v1()
     if dim == TERM_DIM:
+        warn_legacy_absence(*(d["const"] for d in data))
         ds = xr.concat([d[["coeffs", "vars"]] for d in data], dim, **kwargs)
         subkwargs = {**kwargs, "fill_value": {"const": join_fill(fill_value, 0)}}
         const = xr.concat([d[["const"]] for d in data], dim, **subkwargs)["const"]