refac: introduce consistent convention for linopy operations with subsets and supersets

CLAUDE.md

@@ -110,27 +110,6 @@ When modifying the codebase, maintain consistency with these patterns and ensure
 * Always create a feature branch for new features or bug fixes.
 * Use the github cli (gh) to interact with the Github repository.
 
-### GitHub Claude Code Integration
-
-This repository includes Claude Code GitHub Actions for automated assistance:
-
-1. **Automated PR Reviews** (`claude-code-review.yml`):
-   - Automatically reviews PRs only when first created (opened)
-   - Subsequent reviews require manual `@claude` mention
-   - Focuses on Python best practices, xarray patterns, and optimization correctness
-   - Can run tests and linting as part of the review
-   - **Skip initial review by**: Adding `[skip-review]` or `[WIP]` to PR title, or using draft PRs
-
-2. **Manual Claude Assistance** (`claude.yml`):
-   - Trigger by mentioning `@claude` in any:
-     - Issue comments
-     - Pull request comments
-     - Pull request reviews
-     - New issue body or title
-   - Claude can help with bug fixes, feature implementation, code explanations, etc.
-
-**Note**: Both workflows require the `ANTHROPIC_API_KEY` secret to be configured in the repository settings.
-
 
 ## Development Guidelines
 
@@ -140,3 +119,4 @@ This repository includes Claude Code GitHub Actions for automated assistance:
 4. Use type hints and mypy for type checking.
 5. Always write tests into the `test` directory, following the naming convention `test_*.py`.
 6. Always write temporary and non git-tracked code in the `dev-scripts` directory.
+7. In test scripts use linopy assertions from the testing.py module where useful (assert_linequal, assert_varequal, etc.)

doc/release_notes.rst

@@ -6,6 +6,12 @@ Upcoming Version
 
 * Fix docs (pick highs solver)
 * Add the `sphinx-copybutton` to the documentation
+* Harmonize coordinate alignment for operations with subset/superset objects:
+  - Multiplication and division fill missing coords with 0 (variable doesn't participate)
+  - Addition and subtraction of constants fill missing coords with 0 (identity element) and pin result to LHS coords
+  - Comparison operators (``==``, ``<=``, ``>=``) fill missing RHS coords with NaN (no constraint created)
+  - Fixes crash on ``subset + var`` / ``subset + expr`` reverse addition
+  - Fixes superset DataArrays expanding result coords beyond the variable's coordinate space
 * Add ``auto_mask`` parameter to ``Model`` class that automatically masks variables and constraints where bounds, coefficients, or RHS values contain NaN. This eliminates the need to manually create mask arrays when working with sparse or incomplete data.
 * Speed up LP file writing by 2-2.7x on large models through Polars streaming engine, join-based constraint assembly, and reduced per-constraint overhead
 * Fix multiplication of constant-only ``LinearExpression`` with other expressions

linopy/expressions.py

@@ -95,6 +95,7 @@
 
 SUPPORTED_CONSTANT_TYPES = (
     np.number,
+    np.bool_,
     int,
     float,
     DataArray,
@@ -533,15 +534,33 @@ def _multiply_by_linear_expression(
             res = res + self.reset_const() * other.const
         return res
 
+    def _add_constant(
+        self: GenericExpression, other: ConstantLike
+    ) -> GenericExpression:
+        if np.isscalar(other):
+            return self.assign(const=self.const + other)
+        da = as_dataarray(other, coords=self.coords, dims=self.coord_dims)
+        da = da.reindex_like(self.const, fill_value=0)
+        return self.assign(const=self.const + da)
+
     def _multiply_by_constant(
         self: GenericExpression, other: ConstantLike
     ) -> GenericExpression:
         multiplier = as_dataarray(other, coords=self.coords, dims=self.coord_dims)
+        multiplier = multiplier.reindex_like(self.const, fill_value=0)
         coeffs = self.coeffs * multiplier
-        assert all(coeffs.sizes[d] == s for d, s in self.coeffs.sizes.items())
         const = self.const * multiplier
         return self.assign(coeffs=coeffs, const=const)
 
+    def _divide_by_constant(
+        self: GenericExpression, other: ConstantLike
+    ) -> GenericExpression:
+        divisor = as_dataarray(other, coords=self.coords, dims=self.coord_dims)
+        divisor = divisor.reindex_like(self.const, fill_value=1)
+        coeffs = self.coeffs / divisor
+        const = self.const / divisor
+        return self.assign(coeffs=coeffs, const=const)
+
     def __div__(self: GenericExpression, other: SideLike) -> GenericExpression:
         try:
             if isinstance(
@@ -557,7 +576,7 @@ def __div__(self: GenericExpression, other: SideLike) -> GenericExpression:
                     f"{type(self)} and {type(other)}"
                     "Non-linear expressions are not yet supported."
                 )
-            return self._multiply_by_constant(other=1 / other)
+            return self._divide_by_constant(other)
         except TypeError:
             return NotImplemented
 
@@ -863,7 +882,10 @@ def to_constraint(self, sign: SignLike, rhs: SideLike) -> Constraint:
         sign : str, array-like
             Sign(s) of the constraints.
         rhs : constant, Variable, LinearExpression
-            Right-hand side of the constraint.
+            Right-hand side of the constraint. If a DataArray, it is
+            reindexed to match expression coordinates (fill_value=np.nan).
+            Extra dimensions in the RHS not present in the expression
+            raise a ValueError. NaN entries in the RHS mean "no constraint".
 
         Returns
         -------
@@ -876,6 +898,15 @@ def to_constraint(self, sign: SignLike, rhs: SideLike) -> Constraint:
                 f"Both sides of the constraint are constant. At least one side must contain variables. {self} {rhs}"
             )
 
+        if isinstance(rhs, DataArray):
+            extra_dims = set(rhs.dims) - set(self.coord_dims)
+            if extra_dims:
+                raise ValueError(
+                    f"RHS DataArray has dimensions {extra_dims} not present "
+                    f"in the expression. Cannot create constraint."
+                )
+            rhs = rhs.reindex_like(self.const, fill_value=np.nan)
+
         all_to_lhs = (self - rhs).data
         data = assign_multiindex_safe(
             all_to_lhs[["coeffs", "vars"]], sign=sign, rhs=-all_to_lhs.const
@@ -1312,11 +1343,11 @@ def __add__(
             return other.__add__(self)
 
         try:
-            if np.isscalar(other):
-                return self.assign(const=self.const + other)
-
-            other = as_expression(other, model=self.model, dims=self.coord_dims)
-            return merge([self, other], cls=self.__class__)
+            if isinstance(other, SUPPORTED_CONSTANT_TYPES):
+                return self._add_constant(other)
+            else:
+                other = as_expression(other, model=self.model, dims=self.coord_dims)
+                return merge([self, other], cls=self.__class__)
         except TypeError:
             return NotImplemented
 
@@ -1852,15 +1883,15 @@ def __add__(self, other: SideLike) -> QuadraticExpression:
         dimension names of self will be filled in other
         """
         try:
-            if np.isscalar(other):
-                return self.assign(const=self.const + other)
-
-            other = as_expression(other, model=self.model, dims=self.coord_dims)
+            if isinstance(other, SUPPORTED_CONSTANT_TYPES):
+                return self._add_constant(other)
+            else:
+                other = as_expression(other, model=self.model, dims=self.coord_dims)
 
-            if isinstance(other, LinearExpression):
-                other = other.to_quadexpr()
+                if isinstance(other, LinearExpression):
+                    other = other.to_quadexpr()
 
-            return merge([self, other], cls=self.__class__)
+                return merge([self, other], cls=self.__class__)
         except TypeError:
             return NotImplemented
 
@@ -1878,13 +1909,7 @@ def __sub__(self, other: SideLike) -> QuadraticExpression:
         dimension names of self will be filled in other
         """
         try:
-            if np.isscalar(other):
-                return self.assign(const=self.const - other)
-
-            other = as_expression(other, model=self.model, dims=self.coord_dims)
-            if type(other) is LinearExpression:
-                other = other.to_quadexpr()
-            return merge([self, -other], cls=self.__class__)
+            return self.__add__(-other)
         except TypeError:
             return NotImplemented
 

linopy/model.py

@@ -741,6 +741,16 @@ def add_constraints(
             # TODO: add a warning here, routines should be safe against this
             data = data.drop_vars(drop_dims)
 
+        rhs_nan = data.rhs.isnull()
+        if rhs_nan.any():
+            data["rhs"] = data.rhs.fillna(0)
+            rhs_mask = ~rhs_nan
+            mask = (
+                rhs_mask
+                if mask is None
+                else (as_dataarray(mask).astype(bool) & rhs_mask)
+            )
+
         data["labels"] = -1
         (data,) = xr.broadcast(data, exclude=[TERM_DIM])
 

linopy/monkey_patch_xarray.py

@@ -1,37 +1,45 @@
 from __future__ import annotations
 
 from collections.abc import Callable
-from functools import partialmethod, update_wrapper
-from types import NotImplementedType
+from functools import update_wrapper
 from typing import Any
 
 from xarray import DataArray
 
 from linopy import expressions, variables
 
-
-def monkey_patch(cls: type[DataArray], pass_unpatched_method: bool = False) -> Callable:
-    def deco(func: Callable) -> Callable:
-        func_name = func.__name__
-        wrapped = getattr(cls, func_name)
-        update_wrapper(func, wrapped)
-        if pass_unpatched_method:
-            func = partialmethod(func, unpatched_method=wrapped)  # type: ignore
-        setattr(cls, func_name, func)
-        return func
-
-    return deco
-
-
-@monkey_patch(DataArray, pass_unpatched_method=True)
-def __mul__(
-    da: DataArray, other: Any, unpatched_method: Callable
-) -> DataArray | NotImplementedType:
-    if isinstance(
-        other,
-        variables.Variable
-        | expressions.LinearExpression
-        | expressions.QuadraticExpression,
-    ):
-        return NotImplemented
-    return unpatched_method(da, other)
+_LINOPY_TYPES = (
+    variables.Variable,
+    variables.ScalarVariable,
+    expressions.LinearExpression,
+    expressions.ScalarLinearExpression,
+    expressions.QuadraticExpression,
+)
+
+
+def _make_patched_op(op_name: str) -> None:
+    """Patch a DataArray operator to return NotImplemented for linopy types, enabling reflected operators."""
+    original = getattr(DataArray, op_name)
+
+    def patched(
+        da: DataArray, other: Any, unpatched_method: Callable = original
+    ) -> Any:
+        if isinstance(other, _LINOPY_TYPES):
+            return NotImplemented
+        return unpatched_method(da, other)
+
+    update_wrapper(patched, original)
+    setattr(DataArray, op_name, patched)
+
+
+for _op in (
+    "__mul__",
+    "__add__",
+    "__sub__",
+    "__truediv__",
+    "__le__",
+    "__ge__",
+    "__eq__",
+):
+    _make_patched_op(_op)
+del _op

linopy/variables.py

@@ -315,6 +315,7 @@ def to_linexpr(
             Linear expression with the variables and coefficients.
         """
         coefficient = as_dataarray(coefficient, coords=self.coords, dims=self.dims)
+        coefficient = coefficient.reindex_like(self.labels, fill_value=0)
         ds = Dataset({"coeffs": coefficient, "vars": self.labels}).expand_dims(
             TERM_DIM, -1
         )
@@ -454,7 +455,7 @@ def __div__(
                 f"{type(self)} and {type(other)}. "
                 "Non-linear expressions are not yet supported."
             )
-        return self.to_linexpr(1 / other)
+        return self.to_linexpr()._divide_by_constant(other)
 
     def __truediv__(
         self, coefficient: float | int | LinearExpression | Variable