Fix Xpress IIS mapping for masked constraints

doc/release_notes.rst

@@ -9,6 +9,7 @@ Upcoming Version
 * Add ``linopy.breakpoints()`` factory for convenient breakpoint construction from lists, dicts, or keyword arguments. Includes ``breakpoints.segments()`` for disjunctive formulations.
 * Add the `sphinx-copybutton` to the documentation
 * Add SOS1 and SOS2 reformulations for solvers not supporting them.
+* Fix Xpress IIS label mapping for masked constraints and add a regression test for matching infeasible coordinates.
 
 
 Version 0.6.5

linopy/model.py

@@ -241,12 +241,20 @@ def objective(self) -> Objective:
     @objective.setter
     def objective(
         self, obj: Objective | LinearExpression | QuadraticExpression
-    ) -> Objective:
+    ) -> None:
+        """
+        Set the objective function.
+
+        Parameters
+        ----------
+        obj : Objective, LinearExpression, or QuadraticExpression
+            The objective to assign to the model. If not an Objective instance,
+            it will be wrapped in an Objective.
+        """
         if not isinstance(obj, Objective):
             obj = Objective(obj, self)
 
         self._objective = obj
-        return self._objective
 
     @property
     def sense(self) -> str:
@@ -257,6 +265,9 @@ def sense(self) -> str:
 
     @sense.setter
     def sense(self, value: str) -> None:
+        """
+        Set the sense of the objective function.
+        """
         self.objective.sense = value
 
     @property
@@ -271,6 +282,9 @@ def parameters(self) -> Dataset:
 
     @parameters.setter
     def parameters(self, value: Dataset | Mapping) -> None:
+        """
+        Set the parameters of the model.
+        """
         self._parameters = Dataset(value)
 
     @property
@@ -296,6 +310,9 @@ def status(self) -> str:
 
     @status.setter
     def status(self, value: str) -> None:
+        """
+        Set the status of the model.
+        """
         self._status = ModelStatus[value].value
 
     @property
@@ -307,11 +324,13 @@ def termination_condition(self) -> str:
 
     @termination_condition.setter
     def termination_condition(self, value: str) -> None:
-        # TODO: remove if-clause, only kept for backward compatibility
-        if value:
-            self._termination_condition = TerminationCondition[value].value
-        else:
+        """
+        Set the termination condition of the model.
+        """
+        if value == "":
             self._termination_condition = value
+        else:
+            self._termination_condition = TerminationCondition[value].value
 
     @property
     def chunk(self) -> T_Chunks:
@@ -322,6 +341,9 @@ def chunk(self) -> T_Chunks:
 
     @chunk.setter
     def chunk(self, value: T_Chunks) -> None:
+        """
+        Set the chunk sizes of the model.
+        """
         self._chunk = value
 
     @property
@@ -339,6 +361,9 @@ def force_dim_names(self) -> bool:
 
     @force_dim_names.setter
     def force_dim_names(self, value: bool) -> None:
+        """
+        Set whether to force custom dimension names for variables and constraints.
+        """
         self._force_dim_names = bool(value)
 
     @property
@@ -351,6 +376,9 @@ def auto_mask(self) -> bool:
 
     @auto_mask.setter
     def auto_mask(self, value: bool) -> None:
+        """
+        Set whether to automatically mask variables and constraints with NaN values.
+        """
         self._auto_mask = bool(value)
 
     @property
@@ -362,6 +390,9 @@ def solver_dir(self) -> Path:
 
     @solver_dir.setter
     def solver_dir(self, value: str | Path) -> None:
+        """
+        Set the solver directory of the model.
+        """
         if not isinstance(value, str | Path):
             raise TypeError("'solver_dir' must path-like.")
         self._solver_dir = Path(value)
@@ -1638,7 +1669,14 @@ def _compute_infeasibilities_gurobi(self, solver_model: Any) -> list[int]:
         return labels
 
     def _compute_infeasibilities_xpress(self, solver_model: Any) -> list[int]:
-        """Compute infeasibilities for Xpress solver."""
+        """
+        Compute infeasibilities for Xpress solver.
+
+        This function correctly maps solver constraint positions to linopy
+        constraint labels, handling masked constraints where some labels may
+        be skipped (e.g., labels [0, 2, 4] with gaps instead of sequential
+        [0, 1, 2]).
+        """
         # Compute all IIS
         try:  # Try new API first
             solver_model.IISAll()
@@ -1652,20 +1690,21 @@ def _compute_infeasibilities_xpress(self, solver_model: Any) -> list[int]:
 
         labels = set()
 
-        # Create constraint mapping for efficient lookups
-        constraint_to_index = {
-            constraint: idx
-            for idx, constraint in enumerate(solver_model.getConstraint())
-        }
+        clabels = self.matrices.clabels
+        constraint_position_map = {}
+        for position, constraint_obj in enumerate(solver_model.getConstraint()):
+            if 0 <= position < len(clabels):
+                constraint_label = clabels[position]
+                if constraint_label >= 0:
+                    constraint_position_map[constraint_obj] = constraint_label
 
         # Retrieve each IIS
         for iis_num in range(1, num_iis + 1):
             iis_constraints = self._extract_iis_constraints(solver_model, iis_num)
 
-            # Convert constraint objects to indices
             for constraint_obj in iis_constraints:
-                if constraint_obj in constraint_to_index:
-                    labels.add(constraint_to_index[constraint_obj])
+                if constraint_obj in constraint_position_map:
+                    labels.add(constraint_position_map[constraint_obj])
                 # Note: Silently skip constraints not found in mapping
                 # This can happen if the model structure changed after solving
 

linopy/variables.py

@@ -291,9 +291,15 @@ def at(self) -> AtIndexer:
 
     @property
     def loc(self) -> LocIndexer:
+        """
+        Indexing the variable using coordinates.
+        """
         return LocIndexer(self)
 
     def to_pandas(self) -> pd.Series:
+        """
+        Convert the variable labels to a pandas Series.
+        """
         return self.labels.to_pandas()
 
     def to_linexpr(
@@ -734,10 +740,16 @@ def type(self) -> str:
 
     @property
     def coord_dims(self) -> tuple[Hashable, ...]:
+        """
+        Get the coordinate dimensions of the variable.
+        """
         return tuple(k for k in self.dims if k not in HELPER_DIMS)
 
     @property
     def coord_sizes(self) -> dict[Hashable, int]:
+        """
+        Get the coordinate sizes of the variable.
+        """
         return {k: v for k, v in self.sizes.items() if k not in HELPER_DIMS}
 
     @property
@@ -1111,6 +1123,19 @@ def sanitize(self) -> Variable:
         return self
 
     def equals(self, other: Variable) -> bool:
+        """
+        Check if this Variable is equal to another.
+
+        Parameters
+        ----------
+        other : Variable
+            The Variable to compare with.
+
+        Returns
+        -------
+        bool
+            True if the variables have equal labels, False otherwise.
+        """
         return self.labels.equals(other.labels)
 
     # Wrapped function which would convert variable to dataarray

test/test_infeasibility.py

@@ -3,6 +3,8 @@
 Test infeasibility detection for different solvers.
 """
 
+from typing import cast
+
 import pandas as pd
 import pytest
 
@@ -242,3 +244,58 @@ def test_deprecated_method(
 
         # Check that it contains constraint labels
         assert len(subset) > 0
+
+    @pytest.mark.parametrize("solver", ["gurobi", "xpress"])
+    def test_masked_constraint_infeasibility(
+        self, solver: str, capsys: pytest.CaptureFixture[str]
+    ) -> None:
+        """
+        Test infeasibility detection with masked constraints.
+
+        This test verifies that the solver correctly maps constraint positions
+        to constraint labels when constraints are masked (some rows skipped).
+        The enumeration creates positions [0, 1, 2, ...] that should correspond
+        to the actual constraint labels which may have gaps like [0, 2, 4, 6].
+        """
+        if solver not in available_solvers:
+            pytest.skip(f"{solver} not available")
+
+        m = Model()
+
+        time = pd.RangeIndex(8, name="time")
+        x = m.add_variables(lower=0, upper=5, coords=[time], name="x")
+        y = m.add_variables(lower=0, upper=5, coords=[time], name="y")
+
+        # Create a mask that keeps only even time indices (0, 2, 4, 6)
+        mask = pd.Series([i % 2 == 0 for i in range(len(time))])
+        m.add_constraints(x + y >= 10, name="sum_lower", mask=mask)
+
+        mask = pd.Series([False] * (len(time) // 2) + [True] * (len(time) // 2))
+        m.add_constraints(x <= 4, name="x_upper", mask=mask)
+
+        m.add_objective(x.sum() + y.sum())
+        status, condition = m.solve(solver_name=solver)
+
+        assert status == "warning"
+        assert "infeasible" in condition
+
+        labels = m.compute_infeasibilities()
+        assert labels
+
+        positions = [
+            cast(tuple[str, dict[str, int]], m.constraints.get_label_position(label))
+            for label in labels
+        ]
+        grouped_coords: dict[str, set[int]] = {"sum_lower": set(), "x_upper": set()}
+        for name, coord in positions:
+            assert name in grouped_coords
+            grouped_coords[name].add(coord["time"])
+
+        assert grouped_coords["sum_lower"]
+        assert grouped_coords["sum_lower"] == grouped_coords["x_upper"]
+
+        m.print_infeasibilities()
+        output = capsys.readouterr().out
+        for time_coord in grouped_coords["sum_lower"]:
+            assert f"sum_lower[{time_coord}]" in output
+            assert f"x_upper[{time_coord}]" in output