feat: unlabeled-operand pairing by size (#736)

arithmetics-design/convention.md

@@ -122,10 +122,10 @@ it does when no dimension matches. The same goes for a 4×4 array against
 `(a: 4, b: 4)`: sizes cannot tell `(a, b)` from `(b, a)`. To name the
 dimensions, wrap the array in a DataArray.
 
-> **TODO — not yet implemented ([#736]).** Today an unlabeled array pairs
-> with the *leading* dimensions positionally, which silently guesses in the
-> ambiguous cases above. The pairing rule builds on the `as_dataarray` /
-> coords-as-truth seam ([#732], merged — now `linopy.alignment`).
+A scalar broadcasts over every dimension and so needs no pairing. A 0-d
+array is treated as a scalar; a Python `list` is read as a numpy array
+(it carries values, not labels). Implemented in `linopy.alignment`
+([#736]).
 
 ### §8. Shared dimensions must match exactly
 

arithmetics-design/legacy-removal.md

@@ -72,6 +72,11 @@ source tree.
   (`_project_onto_multiindex_levels`, `_LevelProjection`) stays:
   full-coverage full-level projections remain legal under v1 (they are
   the same coordinate spelled differently, §8).
+- `_dims_for_unlabeled_operand`: drop the legacy positional-pairing
+  fallback (the `warn_legacy(...)` branches plus the `return
+  list(candidates)`); the v1 size-pairing — the `is_v1()` block that
+  raises on ambiguity / no-match — becomes the whole function. The
+  `as_constant` / `_pair_axes_by_size` helpers stay (v1-clean).
 
 ### `linopy/piecewise.py` / `linopy/sos_reformulation.py`
 

linopy/examples.py

@@ -2,6 +2,7 @@
 This module contains examples of linear programming models using the linopy library.
 """
 
+import xarray as xr
 from numpy import arange
 
 from linopy import Model
@@ -73,7 +74,7 @@ def benchmark_model(n: int = 10, integerlabels: bool = False) -> Model:
         naxis, maxis = [arange(n, dtype=float), arange(n).astype(str)]
     x = m.add_variables(coords=[naxis, maxis])
     y = m.add_variables(coords=[naxis, maxis])
-    m.add_constraints(x - y >= naxis)
+    m.add_constraints(x - y >= xr.DataArray(naxis, dims=["dim_0"]))
     m.add_constraints(x + y >= 0)
     m.add_objective((2 * x).sum() + y.sum())
     return m

linopy/expressions.py

@@ -44,7 +44,13 @@
 from types import EllipsisType, NotImplementedType
 
 from linopy import constraints, variables
-from linopy.alignment import as_dataarray, broadcast_to_coords, fill_missing_coords
+from linopy.alignment import (
+    _matmul_operand_to_dataarray,
+    as_constant,
+    as_dataarray,
+    broadcast_to_coords,
+    fill_missing_coords,
+)
 from linopy.common import (
     EmptyDeprecationWrapper,
     LocIndexer,
@@ -680,9 +686,7 @@ def _add_constant_v1(
             if isinstance(other, float) and np.isnan(other):
                 check_user_nan()
             return self.assign(const=self.const + other)
-        da = broadcast_to_coords(
-            other, coords=self.coords, dims=self.coord_dims, strict=False
-        )
+        da = broadcast_to_coords(other, coords=self.coords, strict=False)
         if da.isnull().any():
             check_user_nan()
         self_const, da, needs_data_reindex = self._align_constant(
@@ -708,9 +712,7 @@ def _add_constant_legacy(
             if isinstance(other, float) and np.isnan(other):
                 check_user_nan()
             return self.assign(const=self.const.fillna(0) + other)
-        da = broadcast_to_coords(
-            other, coords=self.coords, dims=self.coord_dims, strict=False
-        )
+        da = broadcast_to_coords(other, coords=self.coords, strict=False)
         if da.isnull().any():
             check_user_nan()
         self_const, da, needs_data_reindex = self._align_constant(
@@ -752,9 +754,7 @@ def _apply_constant_op_v1(
         # §5: user NaN raised before we get here.
         if isinstance(other, float) and np.isnan(other):
             check_user_nan(op_kind=op_kind)
-        factor = broadcast_to_coords(
-            other, coords=self.coords, dims=self.coord_dims, strict=False
-        )
+        factor = broadcast_to_coords(other, coords=self.coords, strict=False)
         if factor.isnull().any():
             check_user_nan(op_kind=op_kind)
         self_const, factor, needs_data_reindex = self._align_constant(
@@ -785,9 +785,7 @@ def _apply_constant_op_legacy(
         # factor → fill_value (0 for mul, 1 for div), coeffs/const → 0.
         if isinstance(other, float) and np.isnan(other):
             check_user_nan(op_kind=op_kind)
-        factor = broadcast_to_coords(
-            other, coords=self.coords, dims=self.coord_dims, strict=False
-        )
+        factor = broadcast_to_coords(other, coords=self.coords, strict=False)
         if factor.isnull().any():
             check_user_nan(op_kind=op_kind)
         self_const, factor, needs_data_reindex = self._align_constant(
@@ -822,6 +820,7 @@ def _divide_by_constant(
         )
 
     def __div__(self: GenericExpression, other: SideLike) -> GenericExpression:
+        other = as_constant(other)
         try:
             if isinstance(other, SUPPORTED_EXPRESSION_TYPES):
                 raise TypeError(
@@ -1280,6 +1279,7 @@ def to_constraint(
         which are moved to the left-hand-side and constant values which are moved
         to the right-hand side.
         """
+        rhs = as_constant(rhs)
         if self.is_constant and is_constant(rhs):
             raise ValueError(
                 f"Both sides of the constraint are constant. At least one side must contain variables. {self} {rhs}"
@@ -1294,9 +1294,7 @@ def to_constraint(
             # through ``sub`` into the RHS and reaches downstream
             # auto-mask handling as "no constraint at this row" (§12).
             if isinstance(rhs, CONSTANT_TYPES):
-                rhs = broadcast_to_coords(
-                    rhs, coords=self.coords, dims=self.coord_dims, strict=False
-                )
+                rhs = broadcast_to_coords(rhs, coords=self.coords, strict=False)
                 extra_dims = set(rhs.dims) - set(self.coord_dims)
                 if extra_dims:
                     logger.warning(
@@ -1320,9 +1318,7 @@ def to_constraint(
         # part of normal arithmetic, so we restore the original NaN mask
         # afterward).
         if isinstance(rhs, CONSTANT_TYPES):
-            rhs = broadcast_to_coords(
-                rhs, coords=self.coords, dims=self.coord_dims, strict=False
-            )
+            rhs = broadcast_to_coords(rhs, coords=self.coords, strict=False)
 
             extra_dims = set(rhs.dims) - set(self.coord_dims)
             if extra_dims:
@@ -1844,6 +1840,7 @@ def __add__(
         Note: If other is a numpy array or pandas object without axes names,
         dimension names of self will be filled in other
         """
+        other = as_constant(other)
         if isinstance(other, QuadraticExpression):
             return other.__add__(self)
 
@@ -1879,6 +1876,7 @@ def __sub__(
         | LinearExpression
         | QuadraticExpression,
     ) -> LinearExpression | QuadraticExpression:
+        other = as_constant(other)
         try:
             return self.__add__(-other)
         except TypeError:
@@ -1903,6 +1901,7 @@ def __mul__(
         """
         Multiply the expr by a factor.
         """
+        other = as_constant(other)
         if isinstance(other, QuadraticExpression):
             return other.__rmul__(self)
 
@@ -1948,8 +1947,9 @@ def __matmul__(
         """
         Matrix multiplication with other, similar to xarray dot.
         """
+        other = as_constant(other)
         if not isinstance(other, LinearExpression | variables.Variable):
-            other = as_dataarray(other, coords=self.coords, dims=self.coord_dims)
+            other = _matmul_operand_to_dataarray(other, self.coords, self.coord_dims)
 
         common_dims = list(set(self.coord_dims).intersection(other.dims))
         return (self * other).sum(dim=common_dims)
@@ -2360,6 +2360,7 @@ def __mul__(self, other: SideLike) -> QuadraticExpression:
         """
         Multiply the expr by a factor.
         """
+        other = as_constant(other)
         if isinstance(other, SUPPORTED_EXPRESSION_TYPES):
             raise TypeError(
                 "unsupported operand type(s) for *: "
@@ -2381,6 +2382,7 @@ def __add__(self, other: SideLike) -> QuadraticExpression:
         Note: If other is a numpy array or pandas object without axes names,
         dimension names of self will be filled in other
         """
+        other = as_constant(other)
         try:
             if isinstance(other, CONSTANT_TYPES):
                 return self._add_constant(other)
@@ -2407,6 +2409,7 @@ def __sub__(self, other: SideLike) -> QuadraticExpression:
         Note: If other is a numpy array or pandas object without axes names,
         dimension names of self will be filled in other
         """
+        other = as_constant(other)
         try:
             return self.__add__(-other)
         except TypeError:
@@ -2430,12 +2433,13 @@ def __matmul__(
         """
         Matrix multiplication with other, similar to xarray dot.
         """
+        other = as_constant(other)
         if isinstance(other, SUPPORTED_EXPRESSION_TYPES):
             raise TypeError(
                 "Higher order non-linear expressions are not yet supported."
             )
 
-        other = as_dataarray(other, coords=self.coords, dims=self.coord_dims)
+        other = _matmul_operand_to_dataarray(other, self.coords, self.coord_dims)
         common_dims = list(set(self.coord_dims).intersection(other.dims))
         return (self * other).sum(dim=common_dims)
 

linopy/types.py

@@ -39,6 +39,10 @@
     | pl.Series
 )
 CONSTANT_TYPES: tuple[type, ...] = get_args(ConstantLike)
+
+UnlabeledLike: TypeAlias = numpy.ndarray | list | pl.Series
+UNLABELED_TYPES = get_args(UnlabeledLike)
+
 SignLike: TypeAlias = str | numpy.ndarray | DataArray | Series | DataFrame
 MaskLike: TypeAlias = numpy.ndarray | DataArray | Series | DataFrame
 PathLike: TypeAlias = str | Path

linopy/variables.py

@@ -333,9 +333,7 @@ def to_linexpr(
         linopy.LinearExpression
             Linear expression with the variables and coefficients.
         """
-        coefficient = broadcast_to_coords(
-            coefficient, coords=self.coords, dims=self.dims, strict=False
-        )
+        coefficient = broadcast_to_coords(coefficient, coords=self.coords, strict=False)
         # §5: user-supplied NaN in the coefficient must raise (v1) / warn
         # (legacy) — it's the multiplicative analogue of ``x + nan_data``
         # and otherwise enters the expression silently. The default

test/test_alignment.py

@@ -600,6 +600,25 @@ def test_extra_coord_entries_broadcast_in(self) -> None:
         assert list(da.coords["dim_0"].values) == ["a", "b"]
         assert list(da.coords["dim_2"].values) == ["A", "B"]
 
+    @pytest.mark.v1
+    def test_explicit_dims_bypass_size_pairing(self) -> None:
+        """
+        An explicit ``dims`` is honored positionally, like xarray — even when
+        size-pairing would otherwise be ambiguous (both coords dims size 4).
+        Pins the "infer order only when the user didn't name it" rule.
+        """
+        coords = {"a": [0, 1, 2, 3], "b": [4, 5, 6, 7]}  # both size 4
+
+        # No dims → size-pairing can't decide → raises.
+        with pytest.raises(ValueError, match=r"sizes alone cannot decide"):
+            broadcast_to_coords(np.arange(4), coords=coords, strict=False)
+
+        # Explicit dims=['a'] → the axis is labeled 'a' positionally, no
+        # pairing, no raise (matches xarray's positional dims assignment).
+        da = broadcast_to_coords(np.arange(4), coords=coords, dims=["a"], strict=False)
+        assert set(da.dims) == {"a", "b"}
+        assert (da.sel(b=4).values == np.arange(4)).all()
+
 
 # ---------------------------------------------------------------------------
 # Implicit MultiIndex-level projection — the legacy/v1 fork point

test/test_constraints.py

@@ -144,7 +144,6 @@ def test_constraint_assignment_with_reindex() -> None:
 @pytest.mark.parametrize(
     "rhs_factory",
     [
-        pytest.param(lambda m, v: v, id="numpy"),
         pytest.param(lambda m, v: xr.DataArray(v, dims=["dim_0"]), id="dataarray"),
         pytest.param(lambda m, v: pd.Series(v, index=v), id="series"),
         pytest.param(
@@ -168,10 +167,37 @@ def test_constraint_rhs_lower_dim(rhs_factory: Any) -> None:
     assert c.shape == (10, 10)
 
 
+@pytest.mark.legacy
+def test_constraint_rhs_unlabeled_lower_dim_legacy() -> None:
+    # An unlabeled array rhs pairs positionally with the leading dim under
+    # legacy; both dims are size 10 so the pairing is a size-coincidence.
+    m = Model()
+    naxis = np.arange(10, dtype=float)
+    maxis = np.arange(10).astype(str)
+    x = m.add_variables(coords=[naxis, maxis])
+    y = m.add_variables(coords=[naxis, maxis])
+
+    c = m.add_constraints(x - y >= naxis)
+    assert c.shape == (10, 10)
+
+
+@pytest.mark.v1
+def test_constraint_rhs_unlabeled_lower_dim_ambiguous_raises_v1() -> None:
+    # v1: both dims are size 10, so an unlabeled length-10 rhs cannot be
+    # paired by size — it raises (wrap in a DataArray to disambiguate).
+    m = Model()
+    naxis = np.arange(10, dtype=float)
+    maxis = np.arange(10).astype(str)
+    x = m.add_variables(coords=[naxis, maxis])
+    y = m.add_variables(coords=[naxis, maxis])
+
+    with pytest.raises(ValueError, match=r"sizes alone cannot decide"):
+        m.add_constraints(x - y >= naxis)
+
+
 @pytest.mark.parametrize(
     "rhs_factory",
     [
-        pytest.param(lambda m: np.ones((5, 3)), id="numpy"),
         pytest.param(lambda m: pd.DataFrame(np.ones((5, 3))), id="dataframe"),
     ],
 )
@@ -186,6 +212,29 @@ def test_constraint_rhs_higher_dim_constant_warns(
     assert "dimensions" in caplog.text
 
 
+@pytest.mark.legacy
+def test_constraint_rhs_unlabeled_higher_dim_warns_legacy(caplog: Any) -> None:
+    # Legacy: an unlabeled (5, 3) rhs pairs axis 0 with dim_0 positionally and
+    # broadcasts the extra axis, warning about the extra dimension.
+    m = Model()
+    x = m.add_variables(coords=[range(5)], name="x")
+
+    with caplog.at_level("WARNING", logger="linopy.expressions"):
+        m.add_constraints(x >= np.ones((5, 3)))
+    assert "dimensions" in caplog.text
+
+
+@pytest.mark.v1
+def test_constraint_rhs_unlabeled_higher_dim_raises_v1() -> None:
+    # v1: the (5, 3) rhs has an axis of length 3 matching no dim of x — it
+    # cannot be paired by size and raises.
+    m = Model()
+    x = m.add_variables(coords=[range(5)], name="x")
+
+    with pytest.raises(ValueError, match=r"no unambiguous dimension match"):
+        m.add_constraints(x >= np.ones((5, 3)))
+
+
 def test_constraint_rhs_higher_dim_dataarray_reindexes() -> None:
     """DataArray RHS with extra dims reindexes to expression coords (no raise)."""
     m = Model()