fix: broadcast pandas/DataArray bounds in coords and preserve dim order

doc/release_notes.rst

@@ -52,6 +52,8 @@ Most users should keep calling ``model.solve(...)``. If you want more control, y
 
 **Bug Fixes**
 
+* ``Model.add_variables``: 0.7.0 made ``coords`` (dims, order, and values) the source of truth for ``DataArray`` bounds; this release closes the two remaining gaps. Pandas ``Series`` / ``DataFrame`` bounds missing a dimension are broadcast to ``coords`` instead of being silently dropped (`#709 <https://github.com/PyPSA/linopy/issues/709>`__), and the variable's dimension order always follows ``coords`` regardless of bound type (`#706 <https://github.com/PyPSA/linopy/issues/706>`__).
+* ``add_piecewise_formulation`` now produces a reproducible dimension order in the broadcast breakpoint array. The previous set-based expansion gave a hash-randomized order that varied between processes.
 * SOS constraints on masked variables no longer cause solver-specific failures (Gurobi ``IndexError``, Xpress ``?404 Invalid column number``, LP parse errors, silent set corruption). ``Model.solve()`` and ``Model.to_file()`` now raise a clear ``NotImplementedError`` referring users to `#688 <https://github.com/PyPSA/linopy/issues/688>`__; pass ``reformulate_sos=True`` as a workaround.
 * ``Model.solve(..., reformulate_sos=True)`` now actually reformulates SOS constraints even when the solver supports them natively. Previously it was silently ignored with a warning.
 

linopy/common.py

@@ -9,7 +9,7 @@
 
 import operator
 import os
-from collections.abc import Callable, Generator, Hashable, Iterable, Sequence
+from collections.abc import Callable, Generator, Hashable, Iterable, Mapping, Sequence
 from functools import cached_property, partial, reduce, wraps
 from pathlib import Path
 from typing import TYPE_CHECKING, Any, Generic, TypeVar, overload
@@ -275,6 +275,170 @@ def as_dataarray(
     return arr
 
 
+def _coords_to_dict(
+    coords: Sequence[Sequence | pd.Index] | Mapping,
+) -> dict[str, Any]:
+    """
+    Normalize coords to a dict mapping dim names to coordinate values.
+
+    Entries must be ``pd.Index`` (named or not) or unnamed sequences
+    (``list`` / ``tuple`` / ``range`` / ``np.ndarray``). Other types —
+    notably ``xarray.DataArray`` — raise ``TypeError`` rather than
+    being silently dropped: callers should convert via
+    ``variable.indexes[<dim>]`` (or ``pd.Index(...)``) first.
+    """
+    if isinstance(coords, Mapping):
+        return dict(coords)
+    result: dict[str, Any] = {}
+    for c in coords:
+        if isinstance(c, pd.Index):
+            if c.name:
+                result[c.name] = c
+        elif isinstance(c, list | tuple | range | np.ndarray):
+            pass  # unnamed sequence contributes no named dim
+        else:
+            raise TypeError(
+                f"coords entries must be pd.Index or an unnamed sequence "
+                f"(list / tuple / range / numpy.ndarray); got "
+                f"{type(c).__name__}. For an xarray DataArray coord, pass "
+                f"`variable.indexes[<dim>]` (a pd.Index) instead."
+            )
+    return result
+
+
+def _named_pandas_to_dataarray(arr: pd.Series | pd.DataFrame) -> DataArray | None:
+    """
+    Convert a pandas Series or DataFrame with fully named axes to a DataArray.
+
+    DataFrame columns (and column-MultiIndex levels) are stacked into the row
+    MultiIndex so each axis name becomes its own dimension. Returns ``None``
+    if any axis (or MultiIndex level) is unnamed, so the caller can fall back
+    to ``as_dataarray``.
+    """
+    names = list(arr.index.names)
+    if isinstance(arr, pd.DataFrame):
+        names += list(arr.columns.names)
+    # pd.Index.names entries can be any hashable (tuples, ints, ...). Only
+    # strings map cleanly to xarray dim names; everything else falls through.
+    if any(not isinstance(n, str) for n in names):
+        return None
+
+    if isinstance(arr, pd.DataFrame):
+        arr = arr.stack(list(range(arr.columns.nlevels)), future_stack=True)
+
+    return arr.to_xarray()
+
+
+def as_dataarray_in_coords(arr: Any, coords: Any, **kwargs: Any) -> DataArray:
+    """
+    Coerce ``arr`` into a DataArray that matches ``coords``.
+
+    Strict-coords counterpart to ``as_dataarray``: ``coords`` is the
+    source of truth, so the returned DataArray has the dimensions,
+    dimension order, and coordinate values of ``coords``, regardless
+    of the input type. Pandas inputs with fully named axes are
+    converted via ``to_xarray`` so their axis names map to dimensions;
+    scalars, numpy arrays, and unnamed pandas go through
+    ``as_dataarray``. The result is then validated, expanded over
+    missing dims, and transposed; ``expand_dims`` and ``transpose``
+    are no-ops when the array already matches.
+
+    - Raises ``ValueError`` if the input has dimensions not in
+      ``coords``.
+    - Raises ``ValueError`` if shared dimension coordinates differ in
+      values. Same-values-different-order coordinates are reindexed.
+    """
+    if coords is None:
+        return as_dataarray(arr, coords, **kwargs)
+
+    expected = _coords_to_dict(coords)
+    if not expected:
+        return as_dataarray(arr, coords, **kwargs)
+
+    orig_type_name = type(arr).__name__
+
+    if isinstance(arr, pd.Series | pd.DataFrame):
+        converted = _named_pandas_to_dataarray(arr)
+        if converted is not None:
+            arr = converted
+
+    if not isinstance(arr, DataArray):
+        # numpy/polars/unnamed-pandas inputs are positional — their only
+        # meaningful information is the values; any axis labels are
+        # auto-generated. Default dims to coords' keys so as_dataarray
+        # labels axes correctly (instead of dim_0/dim_1), then re-assign
+        # coords from expected so positional inputs align to coords by
+        # position. A shape mismatch surfaces here as a clear xarray
+        # "conflicting sizes" error rather than a confusing
+        # "coordinates do not match" further down.
+        kwargs.setdefault("dims", list(expected))
+        arr = as_dataarray(arr, coords, **kwargs)
+        # Skip MultiIndex dims — re-assigning a PandasMultiIndex coord emits
+        # a FutureWarning and isn't needed (as_dataarray already used it).
+        arr = arr.assign_coords(
+            {
+                d: expected[d]
+                for d in arr.dims
+                if d in expected and not isinstance(arr.indexes.get(d), pd.MultiIndex)
+            }
+        )
+
+    extra = set(arr.dims) - set(expected)
+    if extra:
+        raise ValueError(
+            f"{orig_type_name} has extra dimensions not in coords: {extra}"
+        )
+
+    for dim, coord_values in expected.items():
+        if dim not in arr.dims:
+            continue
+        if isinstance(arr.indexes.get(dim), pd.MultiIndex):
+            continue
+        expected_idx = (
+            coord_values
+            if isinstance(coord_values, pd.Index)
+            else pd.Index(coord_values)
+        )
+        actual_idx = arr.coords[dim].to_index()
+        if not actual_idx.equals(expected_idx):
+            # Same values, different order → reindex to match expected order
+            if len(actual_idx) == len(expected_idx) and set(actual_idx) == set(
+                expected_idx
+            ):
+                arr = arr.reindex({dim: expected_idx})
+            else:
+                raise ValueError(
+                    f"Coordinates for dimension '{dim}' do not match: "
+                    f"expected {expected_idx.tolist()}, got {actual_idx.tolist()}"
+                )
+
+    # expand_dims prepends new dimensions and their coordinate variables;
+    # the subsequent transpose restores coords order. Both are no-ops when
+    # the array already matches. Reconstruct so the DataArray's coords
+    # iteration order also follows coords (a Dataset built from this picks
+    # up its dim order from coord insertion).
+    expand = {k: v for k, v in expected.items() if k not in arr.dims}
+    if expand:
+        arr = arr.expand_dims(expand)
+
+    target_dims = tuple(d for d in expected if d in arr.dims) + tuple(
+        d for d in arr.dims if d not in expected
+    )
+    arr = arr.transpose(*target_dims)
+
+    coord_order = [c for c in target_dims if c in arr.coords] + [
+        c for c in arr.coords if c not in target_dims
+    ]
+    if list(arr.coords) != coord_order:
+        arr = DataArray(
+            arr.variable,
+            coords={c: arr.coords[c] for c in coord_order},
+            name=arr.name,
+        )
+
+    return arr
+
+
 def broadcast_mask(mask: DataArray, labels: DataArray) -> DataArray:
     """
     Broadcast a boolean mask to match the shape of labels.

linopy/expressions.py

@@ -1844,7 +1844,7 @@ def from_rule(
         cls,
         model: Model,
         rule: Callable,
-        coords: Sequence[Sequence | pd.Index | DataArray] | Mapping | None = None,
+        coords: Sequence[Sequence | pd.Index] | Mapping | None = None,
     ) -> LinearExpression:
         """
         Create a linear expression from a rule and a set of coordinates.

linopy/model.py

@@ -29,6 +29,7 @@
 from linopy import solvers
 from linopy.common import (
     as_dataarray,
+    as_dataarray_in_coords,
     assign_multiindex_safe,
     best_int,
     broadcast_mask,
@@ -112,73 +113,6 @@
 logger = logging.getLogger(__name__)
 
 
-def _coords_to_dict(
-    coords: Sequence[Sequence | pd.Index | DataArray] | Mapping,
-) -> dict[str, Any]:
-    """Normalize coords to a dict mapping dim names to coordinate values."""
-    if isinstance(coords, Mapping):
-        return dict(coords)
-    # Sequence of indexes
-    result: dict[str, Any] = {}
-    for c in coords:
-        if isinstance(c, pd.Index) and c.name:
-            result[c.name] = c
-    return result
-
-
-def _validate_dataarray_bounds(arr: Any, coords: Any) -> Any:
-    """
-    Validate and expand DataArray bounds against explicit coords.
-
-    If ``arr`` is not a DataArray, return it unchanged (``as_dataarray``
-    will handle conversion). For DataArray inputs:
-
-    - Raises ``ValueError`` if the array has dimensions not in coords.
-    - Raises ``ValueError`` if shared dimension coordinates don't match.
-    - Expands missing dimensions via ``expand_dims``.
-    """
-    if not isinstance(arr, DataArray):
-        return arr
-
-    expected = _coords_to_dict(coords)
-    if not expected:
-        return arr
-
-    extra = set(arr.dims) - set(expected)
-    if extra:
-        raise ValueError(f"DataArray has extra dimensions not in coords: {extra}")
-
-    for dim, coord_values in expected.items():
-        if dim not in arr.dims:
-            continue
-        if isinstance(arr.indexes.get(dim), pd.MultiIndex):
-            continue
-        expected_idx = (
-            coord_values
-            if isinstance(coord_values, pd.Index)
-            else pd.Index(coord_values)
-        )
-        actual_idx = arr.coords[dim].to_index()
-        if not actual_idx.equals(expected_idx):
-            # Same values, different order → reindex to match expected order
-            if len(actual_idx) == len(expected_idx) and set(actual_idx) == set(
-                expected_idx
-            ):
-                arr = arr.reindex({dim: expected_idx})
-            else:
-                raise ValueError(
-                    f"Coordinates for dimension '{dim}' do not match: "
-                    f"expected {expected_idx.tolist()}, got {actual_idx.tolist()}"
-                )
-
-    # Expand missing dimensions
-    expand = {k: v for k, v in expected.items() if k not in arr.dims}
-    if expand:
-        arr = arr.expand_dims(expand)
-
-    return arr
-
-
 class Model:
     """
     Linear optimization model.
@@ -657,7 +591,7 @@ def add_variables(
         self,
         lower: Any = -inf,
         upper: Any = inf,
-        coords: Sequence[Sequence | pd.Index | DataArray] | Mapping | None = None,
+        coords: Sequence[Sequence | pd.Index] | Mapping | None = None,
         name: str | None = None,
         mask: DataArray | ndarray | Series | None = None,
         binary: bool = False,
@@ -683,11 +617,13 @@ def add_variables(
             Upper bound of the variable(s). Ignored if `binary` is True.
             The default is inf.
         coords : list/xarray.Coordinates, optional
-            The coords of the variable array.
-            These are directly passed to the DataArray creation of
-            `lower` and `upper`. For every single combination of
-            coordinates a optimization variable is added to the model.
-            The default is None.
+            The coords of the variable array. When provided, ``coords``
+            is the source of truth for the variable's dimensions,
+            dimension order, and coordinate values; ``lower`` and
+            ``upper`` are broadcast and aligned to match. One
+            optimization variable is added per combination of
+            coordinates. The default is None, in which case the shape
+            is inferred from the bounds.
         name : str, optional
             Reference name of the added variables. The default None results in
             a name like "var1", "var2" etc.
@@ -765,14 +701,10 @@ def add_variables(
                     "Semi-continuous variables require a positive scalar lower bound."
                 )
 
-        if coords is not None:
-            lower = _validate_dataarray_bounds(lower, coords)
-            upper = _validate_dataarray_bounds(upper, coords)
-
         data = Dataset(
             {
-                "lower": as_dataarray(lower, coords, **kwargs),
-                "upper": as_dataarray(upper, coords, **kwargs),
+                "lower": as_dataarray_in_coords(lower, coords, **kwargs),
+                "upper": as_dataarray_in_coords(upper, coords, **kwargs),
                 "labels": -1,
             }
         )
@@ -891,7 +823,7 @@ def add_constraints(
         sign: SignLike | None = ...,
         rhs: ConstantLike | VariableLike | ExpressionLike | None = ...,
         name: str | None = ...,
-        coords: Sequence[Sequence | pd.Index | DataArray] | Mapping | None = ...,
+        coords: Sequence[Sequence | pd.Index] | Mapping | None = ...,
         mask: MaskLike | None = ...,
         freeze: Literal[False] = ...,
     ) -> Constraint: ...
@@ -907,7 +839,7 @@ def add_constraints(
         sign: SignLike | None = ...,
         rhs: ConstantLike | VariableLike | ExpressionLike | None = ...,
         name: str | None = ...,
-        coords: Sequence[Sequence | pd.Index | DataArray] | Mapping | None = ...,
+        coords: Sequence[Sequence | pd.Index] | Mapping | None = ...,
         mask: MaskLike | None = ...,
         freeze: Literal[True] = ...,
     ) -> CSRConstraint: ...
@@ -922,7 +854,7 @@ def add_constraints(
         sign: SignLike | None = None,
         rhs: ConstantLike | VariableLike | ExpressionLike | None = None,
         name: str | None = None,
-        coords: Sequence[Sequence | pd.Index | DataArray] | Mapping | None = None,
+        coords: Sequence[Sequence | pd.Index] | Mapping | None = None,
         mask: MaskLike | None = None,
         freeze: bool | None = None,
     ) -> ConstraintBase:
@@ -1428,7 +1360,7 @@ def calculate_block_maps(self) -> None:
 
     @overload
     def linexpr(
-        self, *args: Sequence[Sequence | pd.Index | DataArray] | Mapping
+        self, *args: Sequence[Sequence | pd.Index] | Mapping
     ) -> LinearExpression: ...
 
     @overload
@@ -1441,7 +1373,7 @@ def linexpr(
         *args: tuple[ConstantLike, str | Variable | ScalarVariable]
         | ConstantLike
         | Callable
-        | Sequence[Sequence | pd.Index | DataArray]
+        | Sequence[Sequence | pd.Index]
         | Mapping,
     ) -> LinearExpression:
         """

linopy/piecewise.py

@@ -1006,20 +1006,18 @@ def _broadcast_points(
 
     lin_exprs = [_to_linexpr(e) for e in exprs]
 
-    target_dims: set[str] = set()
-    for le in lin_exprs:
-        target_dims.update(str(d) for d in le.coord_dims)
-
-    missing = target_dims - skip - {str(d) for d in points.dims}
-    if not missing:
-        return points
+    point_dims = {str(d) for d in points.dims}
 
+    # Iterate exprs/dims in order; a set would give a hash-dependent,
+    # run-varying expanded dimension order.
     expand_map: dict[str, list] = {}
-    for d in missing:
-        for le in lin_exprs:
+    for le in lin_exprs:
+        for dim in le.coord_dims:
+            d = str(dim)
+            if d in skip or d in point_dims or d in expand_map:
+                continue
             if d in le.coords:
-                expand_map[str(d)] = list(le.coords[d].values)
-                break
+                expand_map[d] = list(le.coords[d].values)
 
     if expand_map:
         points = points.expand_dims(expand_map)