From ff4e6ab0116f04b701ee45f901d739d45baa2daa Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Sun, 22 Feb 2026 11:10:41 -0600 Subject: [PATCH 01/12] Add indicator constraints --- doc/release_notes.rst | 2 +- linopy/constants.py | 4 + linopy/io.py | 111 ++++++++++++++++ linopy/model.py | 139 +++++++++++++++++++ linopy/solver_capabilities.py | 3 + test/test_indicator_constraints.py | 205 +++++++++++++++++++++++++++++ 6 files changed, 463 insertions(+), 1 deletion(-) create mode 100644 test/test_indicator_constraints.py diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 979b2263..72e5a637 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -5,7 +5,7 @@ Upcoming Version ---------------- * Add SOS1 and SOS2 reformulations for solvers not supporting them. - +* Add indicator constraints for solvers that support them. Version 0.6.4 -------------- diff --git a/linopy/constants.py b/linopy/constants.py index 2e1ef47a..fa6b07d8 100644 --- a/linopy/constants.py +++ b/linopy/constants.py @@ -54,6 +54,10 @@ SOS_DIM_ATTR = "sos_dim" SOS_BIG_M_ATTR = "big_m_upper" +# Indicator constraint attribute keys +INDICATOR_BINARY_VAR_ATTR = "indicator_binary_var" +INDICATOR_BINARY_VAL_ATTR = "indicator_binary_val" + class ModelStatus(Enum): """ diff --git a/linopy/io.py b/linopy/io.py index 54090e87..5863736b 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -399,6 +399,66 @@ def sos_to_file( _format_and_write(df, columns, f) +def indicator_constraints_to_file( + m: Model, + f: BufferedWriter, + explicit_coordinate_names: bool = False, +) -> None: + """ + Write indicator constraints to the s.t. section of an LP file. + + Indicator constraints appear in the Subject To section with the format: + ``ic0: b = 1 -> +1.0 x <= 5.0`` + """ + if not m.indicator_constraints: + return + + # If no regular constraints were written, we need the s.t. header + if not len(m.constraints): + f.write(b"\n\ns.t.\n\n") + + print_variable_scalar, _ = get_printers_scalar( + m, explicit_coordinate_names=explicit_coordinate_names + ) + + for ic_name, ic_data in m.indicator_constraints.items(): + labels_flat = ic_data.labels.values.flatten() + binary_var_flat = ic_data.binary_var.values.flatten() + binary_val_flat = np.broadcast_to( + ic_data.binary_val.values, labels_flat.shape + ).flatten() + coeffs_flat = ic_data.coeffs.values.reshape(len(labels_flat), -1) + vars_flat = ic_data.vars.values.reshape(len(labels_flat), -1) + sign_flat = np.broadcast_to(ic_data.sign.values, labels_flat.shape).flatten() + rhs_flat = np.broadcast_to(ic_data.rhs.values, labels_flat.shape).flatten() + + for i in range(len(labels_flat)): + if labels_flat[i] == -1: + continue + + bvar_name = print_variable_scalar(int(binary_var_flat[i])) + bval = int(binary_val_flat[i]) + + # Build LHS string from coeffs and vars + terms = [] + for c, v in zip(coeffs_flat[i], vars_flat[i]): + if v == -1: + continue + var_name = print_variable_scalar(int(v)) + coeff = float(c) + if coeff >= 0: + terms.append(f"+{coeff} {var_name}") + else: + terms.append(f"{coeff} {var_name}") + + lhs_str = " ".join(terms) + sign_str = str(sign_flat[i]) + rhs_val = float(rhs_flat[i]) + + line = f"ic{labels_flat[i]}: {bvar_name} = {bval} -> {lhs_str} {sign_str} {rhs_val}\n" + f.write(line.encode()) + + def constraints_to_file( m: Model, f: BufferedWriter, @@ -487,6 +547,11 @@ def to_lp_file( slice_size=slice_size, explicit_coordinate_names=explicit_coordinate_names, ) + indicator_constraints_to_file( + m, + f=f, + explicit_coordinate_names=explicit_coordinate_names, + ) bounds_to_file( m, f=f, @@ -757,6 +822,46 @@ def add_sos(s: xr.DataArray, sos_type: int, sos_dim: str) -> None: for _, s in stacked.groupby("_sos_group"): add_sos(s.unstack("_sos_group"), sos_type, sos_dim) + if m.indicator_constraints: + import gurobipy + + sense_map = { + "<=": gurobipy.GRB.LESS_EQUAL, + ">=": gurobipy.GRB.GREATER_EQUAL, + "=": gurobipy.GRB.EQUAL, + } + + for ic_name, ic_data in m.indicator_constraints.items(): + labels_flat = ic_data.labels.values.flatten() + binary_var_flat = ic_data.binary_var.values.flatten() + binary_val_flat = np.broadcast_to( + ic_data.binary_val.values, labels_flat.shape + ).flatten() + coeffs_flat = ic_data.coeffs.values.reshape(len(labels_flat), -1) + vars_flat = ic_data.vars.values.reshape(len(labels_flat), -1) + sign_flat = np.broadcast_to( + ic_data.sign.values, labels_flat.shape + ).flatten() + rhs_flat = np.broadcast_to(ic_data.rhs.values, labels_flat.shape).flatten() + + x_list = x.tolist() + for i in range(len(labels_flat)): + if labels_flat[i] == -1: + continue + + lhs = gurobipy.LinExpr() + for c, v in zip(coeffs_flat[i], vars_flat[i]): + if v != -1: + lhs.add(x_list[v], float(c)) + + model.addGenConstrIndicator( + x_list[int(binary_var_flat[i])], + bool(binary_val_flat[i]), + lhs, + sense_map[str(sign_flat[i])], + float(rhs_flat[i]), + ) + model.update() return model @@ -784,6 +889,12 @@ def to_highspy(m: Model, explicit_coordinate_names: bool = False) -> Highs: "Use io_api='lp' instead." ) + if m.indicator_constraints: + raise NotImplementedError( + "Indicator constraints are not supported by the HiGHS direct API. " + "Use a solver that supports them (gurobi, cplex)." + ) + import highspy print_variable, print_constraint = get_printers_scalar( diff --git a/linopy/model.py b/linopy/model.py index e72b3efa..8f6b6803 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -126,6 +126,7 @@ class Model: # containers "_variables", "_constraints", + "_indicator_constraints", "_objective", "_parameters", "_solution", @@ -136,8 +137,10 @@ class Model: # TODO: move counters to Variables and Constraints class "_xCounter", "_cCounter", + "_icCounter", "_varnameCounter", "_connameCounter", + "_icnameCounter", "_blocks", # TODO: check if these should not be mutable "_chunk", @@ -185,6 +188,7 @@ def __init__( """ self._variables: Variables = Variables({}, model=self) self._constraints: Constraints = Constraints({}, model=self) + self._indicator_constraints: dict[str, Dataset] = {} self._objective: Objective = Objective(LinearExpression(None, self), self) self._parameters: Dataset = Dataset() @@ -192,8 +196,10 @@ def __init__( self._termination_condition: str = "" self._xCounter: int = 0 self._cCounter: int = 0 + self._icCounter: int = 0 self._varnameCounter: int = 0 self._connameCounter: int = 0 + self._icnameCounter: int = 0 self._blocks: DataArray | None = None self._chunk: T_Chunks = chunk @@ -219,6 +225,16 @@ def constraints(self) -> Constraints: """ return self._constraints + @property + def indicator_constraints(self) -> dict[str, Dataset]: + """ + Indicator constraints assigned to the model. + + Each entry is a Dataset with fields: coeffs, vars, sign, rhs, + labels, binary_var, binary_val. + """ + return self._indicator_constraints + @property def objective(self) -> Objective: """ @@ -814,6 +830,122 @@ def add_constraints( self.constraints.add(constraint) return constraint + def add_indicator_constraints( + self, + binary_var: Variable, + binary_val: int, + lhs: ConstraintLike | ExpressionLike | VariableLike, + sign: SignLike | None = None, + rhs: ConstantLike | None = None, + name: str | None = None, + ) -> Dataset: + """ + Add indicator constraints to the model. + + An indicator constraint has the form: + (binary_var == binary_val) => (linear_constraint) + + The linear constraint is only enforced when binary_var equals + binary_val. These constraints are handled natively by solvers + like Gurobi and CPLEX via general constraints. + + Parameters + ---------- + binary_var : linopy.Variable + Binary variable serving as the indicator. Must have binary=True. + binary_val : int + Triggering value, must be 0 or 1. + lhs : linopy.Constraint, linopy.LinearExpression, or linopy.Variable + The conditionally enforced constraint. If a LinearExpression or + Variable is passed, ``sign`` and ``rhs`` must also be provided. + sign : str, optional + Constraint sign ('<=', '>=', '='). Required when ``lhs`` is an + expression. + rhs : numeric, optional + Right-hand side. Required when ``lhs`` is an expression. + name : str, optional + Name for the indicator constraint group. + + Returns + ------- + xarray.Dataset + The stored indicator constraint data. + """ + if not binary_var.attrs.get("binary", False): + raise ValueError( + "Indicator variable must be binary. " + f"Variable '{binary_var.name}' is not binary." + ) + + if binary_val not in (0, 1): + raise ValueError(f"binary_val must be 0 or 1, got {binary_val}.") + + if name in self._indicator_constraints: + raise ValueError( + f"Indicator constraint '{name}' already assigned to model." + ) + if name is None: + name = f"indcon{self._icnameCounter}" + self._icnameCounter += 1 + + # Build constraint data from lhs + if isinstance(lhs, Constraint): + if sign is not None or rhs is not None: + raise ValueError("sign and rhs must be None when lhs is a Constraint.") + data = lhs.data + elif isinstance(lhs, LinearExpression): + if sign is None or rhs is None: + raise ValueError( + "sign and rhs are required when lhs is a LinearExpression." + ) + sign = maybe_replace_signs(as_dataarray(sign)) + data = lhs.to_constraint(sign, rhs).data + elif isinstance(lhs, Variable | ScalarVariable | ScalarLinearExpression): + if sign is None or rhs is None: + raise ValueError("sign and rhs are required when lhs is a Variable.") + sign = maybe_replace_signs(as_dataarray(sign)) + data = lhs.to_linexpr().to_constraint(sign, rhs).data + elif isinstance(lhs, AnonymousScalarConstraint): + if sign is not None or rhs is not None: + raise ValueError("sign and rhs must be None when lhs is a Constraint.") + data = lhs.to_constraint().data + else: + raise TypeError( + f"lhs must be a Constraint, LinearExpression, or Variable, " + f"got {type(lhs)}." + ) + + # Add binary variable labels and value + binary_var_labels = binary_var.labels + data["binary_var"] = binary_var_labels + data["binary_val"] = binary_val + + # Broadcast all fields together + (data,) = xr.broadcast(data, exclude=[TERM_DIM]) + + # Assign unique labels with the same shape as rhs (non-term dims) + labels = DataArray( + np.full(data.rhs.shape, -1, dtype=int), + coords=data.rhs.coords, + dims=data.rhs.dims, + ) + data["labels"] = labels + start = self._icCounter + end = start + data.labels.size + data.labels.values = np.arange(start, end).reshape(data.labels.shape) + self._icCounter += data.labels.size + + data = data.assign_attrs(label_range=(start, end), name=name) + + self._indicator_constraints[name] = data + return data + + def remove_indicator_constraints(self, name: str) -> None: + """ + Remove indicator constraint by name. + """ + del self._indicator_constraints[name] + def add_objective( self, expr: Variable @@ -1408,6 +1540,13 @@ def solve( "Use reformulate_sos=True or a solver that supports SOS (gurobi, cplex)." ) + if self.indicator_constraints: + if not solver_supports(solver_name, SolverFeature.INDICATOR_CONSTRAINTS): + raise ValueError( + f"Solver {solver_name} does not support indicator constraints. " + "Use a solver that supports them (gurobi, cplex)." + ) + try: solver_class = getattr(solvers, f"{solvers.SolverName(solver_name).name}") # initialize the solver as object of solver subclass diff --git a/linopy/solver_capabilities.py b/linopy/solver_capabilities.py index f0507317..9e543861 100644 --- a/linopy/solver_capabilities.py +++ b/linopy/solver_capabilities.py @@ -49,6 +49,7 @@ class SolverFeature(Enum): # Special constraint types SOS_CONSTRAINTS = auto() # Special Ordered Sets (SOS1/SOS2) constraints + INDICATOR_CONSTRAINTS = auto() # Indicator (conditional) constraints # Solver-specific SOLVER_ATTRIBUTE_ACCESS = auto() # Direct access to solver variable attributes @@ -86,6 +87,7 @@ def supports(self, feature: SolverFeature) -> bool: SolverFeature.SOLUTION_FILE_NOT_NEEDED, SolverFeature.IIS_COMPUTATION, SolverFeature.SOS_CONSTRAINTS, + SolverFeature.INDICATOR_CONSTRAINTS, SolverFeature.SOLVER_ATTRIBUTE_ACCESS, } ), @@ -134,6 +136,7 @@ def supports(self, feature: SolverFeature) -> bool: SolverFeature.LP_FILE_NAMES, SolverFeature.READ_MODEL_FROM_FILE, SolverFeature.SOS_CONSTRAINTS, + SolverFeature.INDICATOR_CONSTRAINTS, } ), ), diff --git a/test/test_indicator_constraints.py b/test/test_indicator_constraints.py new file mode 100644 index 00000000..2b404168 --- /dev/null +++ b/test/test_indicator_constraints.py @@ -0,0 +1,205 @@ +"""Tests for indicator constraint support.""" + +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +from linopy import Model, available_solvers + + +def test_add_indicator_constraints_basic() -> None: + """Indicator constraint is created with correct fields.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + assert "ic0" in m.indicator_constraints + ic = m.indicator_constraints["ic0"] + assert "coeffs" in ic + assert "vars" in ic + assert "sign" in ic + assert "rhs" in ic + assert "binary_var" in ic + assert "binary_val" in ic + assert "labels" in ic + + +def test_add_indicator_constraints_from_constraint() -> None: + """Indicator constraint accepts a Constraint object as lhs.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + con = x <= 5 + m.add_indicator_constraints(b, 1, con, name="ic0") + assert "ic0" in m.indicator_constraints + + +def test_indicator_constraints_validation_non_binary() -> None: + """Non-binary variable raises ValueError.""" + m = Model() + y = m.add_variables(lower=0, upper=1, name="y") + x = m.add_variables(lower=0, upper=10, name="x") + with pytest.raises(ValueError, match="must be binary"): + m.add_indicator_constraints(y, 1, x, "<=", 5) + + +def test_indicator_constraints_validation_bad_value() -> None: + """Invalid binary_val raises ValueError.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + with pytest.raises(ValueError, match="must be 0 or 1"): + m.add_indicator_constraints(b, 2, x, "<=", 5) + + +def test_indicator_constraints_duplicate_name() -> None: + """Duplicate name raises ValueError.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + with pytest.raises(ValueError, match="already assigned"): + m.add_indicator_constraints(b, 0, x, ">=", 0, name="ic0") + + +def test_indicator_constraints_not_in_regular_constraints() -> None: + """Indicator constraints are separate from regular constraints.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_constraints(x >= 0, name="regular") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + assert "regular" in m.constraints + assert "ic0" not in m.constraints + assert "ic0" in m.indicator_constraints + + +def test_remove_indicator_constraints() -> None: + """Indicator constraints can be removed.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + assert "ic0" in m.indicator_constraints + m.remove_indicator_constraints("ic0") + assert "ic0" not in m.indicator_constraints + + +def test_indicator_constraints_lp_file(tmp_path: Path) -> None: + """LP file contains general constraints section.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_constraints(x >= 0, name="dummy") + m.add_objective(x) + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + fn = tmp_path / "test.lp" + m.to_file(fn) + content = fn.read_text() + assert "= 1 ->" in content + assert "ic0:" in content + + +@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") +def test_indicator_constraints_solve_gurobi_active() -> None: + """ + Indicator constraint enforced when binary is at trigger value. + + Maximize x subject to: + b = 1 (fixed) + (b == 1) => (x <= 5) + x in [0, 10] + Since b=1 and the indicator enforces x<=5, optimal x=5. + """ + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + # Force b = 1 + m.add_constraints(b >= 1, name="fix_b") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + m.add_objective(x, sense="max") + m.solve(solver_name="gurobi") + assert m.objective.value is not None + assert np.isclose(m.objective.value, 5, atol=1e-6) + + +@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") +def test_indicator_constraints_solve_gurobi_inactive() -> None: + """ + Indicator constraint NOT enforced when binary differs from trigger. + + Maximize x subject to: + b = 1 (fixed) + (b == 0) => (x <= 5) + x in [0, 10] + Since b=1 but trigger is 0, the constraint is inactive. Optimal x=10. + """ + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + # Force b = 1 + m.add_constraints(b >= 1, name="fix_b") + m.add_indicator_constraints(b, 0, x, "<=", 5, name="ic0") + m.add_objective(x, sense="max") + m.solve(solver_name="gurobi") + assert m.objective.value is not None + assert np.isclose(m.objective.value, 10, atol=1e-6) + + +@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") +def test_indicator_constraints_solve_gurobi_multiple() -> None: + """ + Multiple indicator constraints work together. + + b1=1, b2=1 (forced), x in [0, 20]. + (b1 == 1) => (x <= 10) + (b2 == 1) => (x <= 5) + Maximize x. + Both constraints active, so x <= min(10, 5) = 5. Optimal x=5. + """ + m = Model() + b1 = m.add_variables(name="b1", binary=True) + b2 = m.add_variables(name="b2", binary=True) + x = m.add_variables(lower=0, upper=20, name="x") + m.add_constraints(b1 >= 1, name="fix_b1") + m.add_constraints(b2 >= 1, name="fix_b2") + m.add_indicator_constraints(b1, 1, x, "<=", 10, name="ic1") + m.add_indicator_constraints(b2, 1, x, "<=", 5, name="ic2") + m.add_objective(x, sense="max") + m.solve(solver_name="gurobi") + assert m.objective.value is not None + assert np.isclose(m.objective.value, 5, atol=1e-6) + + +def test_unsupported_solver_raises() -> None: + """Solvers without indicator support raise ValueError.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_constraints(x >= 0, name="dummy") + m.add_objective(x) + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + for solver in ["glpk", "highs", "mosek", "mindopt"]: + if solver in available_solvers: + with pytest.raises( + ValueError, match="does not support indicator constraints" + ): + m.solve(solver_name=solver) + + +def test_indicator_constraints_with_coords() -> None: + """Indicator constraints work with multi-dimensional coords.""" + m = Model() + idx = pd.RangeIndex(3, name="i") + b = m.add_variables(coords=[idx], name="b", binary=True) + x = m.add_variables(coords=[idx], lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + ic = m.indicator_constraints["ic0"] + # Should have 3 indicator constraints (one per index) + assert ic.labels.size == 3 From b395e656ced5dc51974f62db6eec24c973f55a72 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 2 Jun 2026 10:41:49 +0200 Subject: [PATCH 02/12] Phase A: indicator fields in unified constraint data model Add is_indicator/binary_var/binary_val to ConstraintBase, Constraint, and CSRConstraint (slots, init, from_mutable, data, netcdf round-trip). Add Constraints.indicator/.regular container accessors. --- linopy/constraints.py | 111 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 109 insertions(+), 2 deletions(-) diff --git a/linopy/constraints.py b/linopy/constraints.py index b74dee5c..952437d3 100644 --- a/linopy/constraints.py +++ b/linopy/constraints.py @@ -179,6 +179,21 @@ def dual(self) -> DataArray: def dual(self, value: DataArray) -> None: """Set the dual values DataArray.""" + @property + @abstractmethod + def is_indicator(self) -> bool: + """Whether the constraint is an indicator constraint.""" + + @property + @abstractmethod + def binary_var(self) -> DataArray | None: + """Get the indicator binary variable labels, or None.""" + + @property + @abstractmethod + def binary_val(self) -> int | np.ndarray | None: + """Get the indicator triggering value(s), or None.""" + @abstractmethod def has_variable(self, variable: variables.Variable) -> bool: """Check if the constraint references any of the given variable labels.""" @@ -509,6 +524,8 @@ class CSRConstraint(ConstraintBase): "_name", "_cindex", "_dual", + "_binvar_labels", + "_binval", ) def __init__( @@ -522,6 +539,8 @@ def __init__( name: str = "", cindex: int | None = None, dual: np.ndarray | None = None, + binvar_labels: np.ndarray | None = None, + binval: int | np.ndarray | None = None, ) -> None: self._csr = csr self._con_labels = con_labels @@ -532,6 +551,8 @@ def __init__( self._name = name self._cindex = cindex self._dual = dual + self._binvar_labels = binvar_labels + self._binval = binval @property def model(self) -> Model: @@ -686,6 +707,20 @@ def dual(self, value: DataArray) -> None: vals = np.asarray(value).ravel() self._dual = vals[self.active_positions] + @property + def is_indicator(self) -> bool: + return self._binvar_labels is not None + + @property + def binary_var(self) -> DataArray | None: + if self._binvar_labels is None: + return None + return self._active_to_dataarray(self._binvar_labels, fill=-1) + + @property + def binary_val(self) -> int | np.ndarray | None: + return self._binval + def _to_dataset(self, nterm: int) -> Dataset: """ Reconstruct labels/coeffs/vars Dataset from the CSR matrix. @@ -743,9 +778,18 @@ def _to_dataset(self, nterm: int) -> Dataset: def data(self) -> Dataset: """Reconstruct the xarray Dataset from the CSR representation.""" ds = self._to_dataset(self.nterm) - extra = {"sign": self.sign, "rhs": self.rhs} + extra: dict[str, Any] = {"sign": self.sign, "rhs": self.rhs} if self._dual is not None: extra["dual"] = self._active_to_dataarray(self._dual, fill=np.nan) + if self._binvar_labels is not None: + extra["binary_var"] = self._active_to_dataarray( + self._binvar_labels, fill=-1 + ) + binval = self._binval + if isinstance(binval, np.ndarray) and binval.ndim > 0: + extra["binary_val"] = self._active_to_dataarray(binval, fill=-1) + else: + extra["binary_val"] = binval return assign_multiindex_safe(ds, **extra).assign_attrs(self.attrs) def __repr__(self) -> str: @@ -833,6 +877,14 @@ def to_netcdf_ds(self) -> Dataset: } if isinstance(self._sign, str): attrs["sign"] = self._sign + if self._binvar_labels is not None: + attrs["is_indicator"] = True + data_vars["_binvar_labels"] = DataArray(self._binvar_labels, dims=["_flat"]) + binval = self._binval + if isinstance(binval, np.ndarray) and binval.ndim > 0: + data_vars["_binval"] = DataArray(binval, dims=["_flat"]) + else: + attrs["binval"] = int(binval) # type: ignore[arg-type] return Dataset(data_vars, attrs=attrs) @classmethod @@ -859,8 +911,23 @@ def from_netcdf_ds(cls, ds: Dataset, model: Model, name: str) -> CSRConstraint: con_labels = np.arange(cindex, cindex + len(rhs), dtype=np.intp) else: con_labels = np.arange(len(rhs), dtype=np.intp) + binvar_labels: np.ndarray | None = None + binval: int | np.ndarray | None = None + if "_binvar_labels" in ds: + binvar_labels = ds["_binvar_labels"].values + binval = ds["_binval"].values if "_binval" in ds else attrs["binval"] return cls( - csr, con_labels, rhs, sign, coords, model, name, cindex=cindex, dual=dual + csr, + con_labels, + rhs, + sign, + coords, + model, + name, + cindex=cindex, + dual=dual, + binvar_labels=binvar_labels, + binval=binval, ) def has_variable(self, variable: variables.Variable) -> bool: @@ -1023,6 +1090,16 @@ def from_mutable( dual = ( con.data["dual"].values.ravel()[active_mask] if "dual" in con.data else None ) + binvar_labels: np.ndarray | None = None + binval: int | np.ndarray | None = None + binary_var = con.binary_var + if binary_var is not None: + binvar_labels = binary_var.values.ravel()[active_mask] + bv = con.binary_val + if isinstance(bv, np.ndarray) and bv.ndim > 0: + binval = bv.ravel()[active_mask] + else: + binval = bv return cls( csr, con_labels, @@ -1033,6 +1110,8 @@ def from_mutable( con.name, cindex=cindex, dual=dual, + binvar_labels=binvar_labels, + binval=binval, ) @@ -1157,6 +1236,20 @@ def rhs(self, value: ExpressionLike) -> None: self.lhs = self.lhs - value.reset_const() self._data = assign_multiindex_safe(self.data, rhs=value.const) + @property + def is_indicator(self) -> bool: + return "binary_var" in self._data + + @property + def binary_var(self) -> DataArray | None: + return self._data.get("binary_var") + + @property + def binary_val(self) -> int | np.ndarray | None: + if "binary_val" in self._data: + return self._data["binary_val"].values + return None + @property def lhs(self) -> expressions.LinearExpression: data = self.data[["coeffs", "vars"]].rename({self.term_dim: TERM_DIM}) @@ -1677,6 +1770,20 @@ def equalities(self) -> Constraints: """ return self[[n for n, s in self.items() if (s.sign == EQUAL).all()]] + @property + def indicator(self) -> Constraints: + """ + Get the subset of constraints which are indicator constraints. + """ + return self[[n for n, c in self.items() if c.is_indicator]] + + @property + def regular(self) -> Constraints: + """ + Get the subset of constraints which are not indicator constraints. + """ + return self[[n for n, c in self.items() if not c.is_indicator]] + def sanitize_zeros(self) -> None: """ Filter out terms with zero and close-to-zero coefficient. From 59c5390f368678642c824883fa2386323a7a4904 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 2 Jun 2026 10:46:43 +0200 Subject: [PATCH 03/12] feat(indicator): unify indicator constraints into normal container (model+common+matrices) --- linopy/common.py | 6 +++- linopy/matrices.py | 70 ++++++++++++++++++++++++++++++++-------------- linopy/model.py | 58 +++++++++++++++----------------------- 3 files changed, 76 insertions(+), 58 deletions(-) diff --git a/linopy/common.py b/linopy/common.py index 831b6bc8..3b471e7f 100644 --- a/linopy/common.py +++ b/linopy/common.py @@ -1014,7 +1014,11 @@ def __init__(self, constraints: Any) -> None: @cached_property def clabels(self) -> np.ndarray: """Active constraint labels in build order, shape (n_active_cons,).""" - label_lists = [c.active_labels() for c in self._constraints.data.values()] + label_lists = [ + c.active_labels() + for c in self._constraints.data.values() + if not c.is_indicator + ] return ( np.concatenate(label_lists) if label_lists else np.array([], dtype=np.intp) ) diff --git a/linopy/matrices.py b/linopy/matrices.py index e694a720..bc1d9328 100644 --- a/linopy/matrices.py +++ b/linopy/matrices.py @@ -71,29 +71,55 @@ def _build_vars(self) -> None: def _build_cons(self) -> None: m = self._parent - - if not len(m.constraints): - self.clabels: ndarray = np.array([], dtype=np.intp) - self.b: ndarray = np.array([]) - self.sense: ndarray = np.array([], dtype=object) - self.A: scipy.sparse.csr_array | None = None - return - label_index = m.variables.label_index - csrs = [] - b_list = [] - sense_list = [] + label_to_pos = label_index.label_to_pos + + reg_csrs, reg_b, reg_sense = [], [], [] + ind_csrs, ind_b, ind_sense, ind_binvar, ind_binval = [], [], [], [], [] for c in m.constraints.data.values(): - csr, _, b, sense = c.to_matrix_with_rhs(label_index) - csrs.append(csr) - b_list.append(b) - sense_list.append(sense) - - self.A = cast(scipy.sparse.csr_array, scipy.sparse.vstack(csrs, format="csr")) - self.clabels = m.constraints.label_index.clabels - self.b = np.concatenate(b_list) if b_list else np.array([]) - self.sense = ( - np.concatenate(sense_list) if sense_list else np.array([], dtype=object) + if c.is_indicator: + cc = c if isinstance(c, CSRConstraint) else c.freeze() + csr, _, b, sense = cc.to_matrix_with_rhs(label_index) + ind_csrs.append(csr) + ind_b.append(b) + ind_sense.append(sense) + ind_binvar.append(label_to_pos[cc._binvar_labels]) + binval = cc._binval + n = len(b) + if np.ndim(binval) == 0: + ind_binval.append(np.full(n, int(binval), dtype=np.intp)) + else: + ind_binval.append(np.asarray(binval, dtype=np.intp).ravel()) + else: + csr, _, b, sense = c.to_matrix_with_rhs(label_index) + reg_csrs.append(csr) + reg_b.append(b) + reg_sense.append(sense) + + self.clabels: ndarray = m.constraints.label_index.clabels + self.A: scipy.sparse.csr_array | None = ( + cast(scipy.sparse.csr_array, scipy.sparse.vstack(reg_csrs, format="csr")) + if reg_csrs + else None + ) + self.b: ndarray = np.concatenate(reg_b) if reg_b else np.array([]) + self.sense: ndarray = ( + np.concatenate(reg_sense) if reg_sense else np.array([], dtype=object) + ) + self.indicator_A: scipy.sparse.csr_array | None = ( + cast(scipy.sparse.csr_array, scipy.sparse.vstack(ind_csrs, format="csr")) + if ind_csrs + else None + ) + self.indicator_b: ndarray = np.concatenate(ind_b) if ind_b else np.array([]) + self.indicator_sense: ndarray = ( + np.concatenate(ind_sense) if ind_sense else np.array([], dtype=object) + ) + self.indicator_binvar: ndarray = ( + np.concatenate(ind_binvar) if ind_binvar else np.array([], dtype=np.intp) + ) + self.indicator_binval: ndarray = ( + np.concatenate(ind_binval) if ind_binval else np.array([], dtype=np.intp) ) @property @@ -156,6 +182,8 @@ def dual(self) -> ndarray: dual_list = [] has_dual = False for c in m.constraints.data.values(): + if c.is_indicator: + continue if isinstance(c, CSRConstraint): # _dual is active-only if c._dual is not None: diff --git a/linopy/model.py b/linopy/model.py index 23d9d876..e086b99d 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -219,7 +219,6 @@ class Model: # containers "_variables", "_constraints", - "_indicator_constraints", "_objective", "_parameters", "_solution", @@ -230,10 +229,8 @@ class Model: # TODO: move counters to Variables and Constraints class "_xCounter", "_cCounter", - "_icCounter", "_varnameCounter", "_connameCounter", - "_icnameCounter", "_pwlCounter", "_blocks", # TODO: check if these should not be mutable @@ -294,7 +291,6 @@ def __init__( """ self._variables: Variables = Variables({}, model=self) self._constraints: Constraints = Constraints({}, model=self) - self._indicator_constraints: dict[str, Dataset] = {} self._objective: Objective = Objective(LinearExpression(None, self), self) self._parameters: Dataset = Dataset() @@ -302,10 +298,8 @@ def __init__( self._termination_condition: str = "" self._xCounter: int = 0 self._cCounter: int = 0 - self._icCounter: int = 0 self._varnameCounter: int = 0 self._connameCounter: int = 0 - self._icnameCounter: int = 0 self._pwlCounter: int = 0 self._blocks: DataArray | None = None @@ -372,14 +366,14 @@ def constraints(self) -> Constraints: return self._constraints @property - def indicator_constraints(self) -> dict[str, Dataset]: + def indicator_constraints(self) -> Constraints: """ Indicator constraints assigned to the model. - Each entry is a Dataset with fields: coeffs, vars, sign, rhs, - labels, binary_var, binary_val. + Returns the subset of ``model.constraints`` for which + ``is_indicator`` is True. """ - return self._indicator_constraints + return self.constraints.indicator @property def objective(self) -> Objective: @@ -1114,7 +1108,7 @@ def add_indicator_constraints( sign: SignLike | None = None, rhs: ConstantLike | None = None, name: str | None = None, - ) -> Dataset: + ) -> Constraint: """ Add indicator constraints to the model. @@ -1144,8 +1138,8 @@ def add_indicator_constraints( Returns ------- - xarray.Dataset - The stored indicator constraint data. + linopy.Constraint + The added indicator constraint. """ if not binary_var.attrs.get("binary", False): raise ValueError( @@ -1156,13 +1150,11 @@ def add_indicator_constraints( if binary_val not in (0, 1): raise ValueError(f"binary_val must be 0 or 1, got {binary_val}.") - if name in self._indicator_constraints: - raise ValueError( - f"Indicator constraint '{name}' already assigned to model." - ) - if name is None: - name = f"indcon{self._icnameCounter}" - self._icnameCounter += 1 + if name in list(self.constraints): + raise ValueError(f"Constraint '{name}' already assigned to model.") + elif name is None: + name = f"indcon{self._connameCounter}" + self._connameCounter += 1 # Build constraint data from lhs if isinstance(lhs, Constraint): @@ -1191,36 +1183,28 @@ def add_indicator_constraints( f"got {type(lhs)}." ) - # Add binary variable labels and value - binary_var_labels = binary_var.labels - data["binary_var"] = binary_var_labels + data["binary_var"] = binary_var.labels data["binary_val"] = binary_val - # Broadcast all fields together + data["labels"] = -1 (data,) = xr.broadcast(data, exclude=[TERM_DIM]) - # Assign unique labels with the same shape as rhs (non-term dims) - labels = DataArray( - np.full(data.rhs.shape, -1, dtype=int), - coords=data.rhs.coords, - dims=data.rhs.dims, - ) - data["labels"] = labels - start = self._icCounter + start = self._cCounter end = start + data.labels.size data.labels.values = np.arange(start, end).reshape(data.labels.shape) - self._icCounter += data.labels.size + self._cCounter += data.labels.size data = data.assign_attrs(label_range=(start, end), name=name) - self._indicator_constraints[name] = data - return data + con = Constraint(data, name=name, model=self, skip_broadcast=True) + freeze = self.freeze_constraints + return self.constraints.add(con, freeze=freeze and not self.chunk) def remove_indicator_constraints(self, name: str) -> None: """ Remove indicator constraint by name. """ - del self._indicator_constraints[name] + self.constraints.remove(name) def add_objective( self, @@ -2013,6 +1997,8 @@ def assign_result( if len(result.solution.dual): dual = result.solution.dual for _, con in self.constraints.items(): + if con.is_indicator: + continue start, end = con.range coords = {dim: con.coords[dim] for dim in con.coord_dims} con.dual = xr.DataArray( From 80c7ec98dd6134969a1ec6666d2d6a606f67b4ee Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 2 Jun 2026 10:48:03 +0200 Subject: [PATCH 04/12] feat(indicator): route LP write and model copy through unified container (io) --- linopy/io.py | 35 ++++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/linopy/io.py b/linopy/io.py index 0c1bfe45..7e1e99cf 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -470,17 +470,18 @@ def indicator_constraints_to_file( Indicator constraints appear in the Subject To section with the format: ``ic0: x0 = 1 -> +1.0 x1 <= 5.0`` """ - if not m.indicator_constraints: + if not len(m.constraints.indicator): return - if not len(m.constraints): + if not len(m.constraints.regular): f.write(b"\n\ns.t.\n\n") print_variable_scalar, _ = get_printers_scalar( m, explicit_coordinate_names=explicit_coordinate_names ) - for ic_data in m.indicator_constraints.values(): + for con in m.constraints.indicator.data.values(): + ic_data = con.data labels_flat = ic_data.labels.values.flatten() binary_var_flat = ic_data.binary_var.values.flatten() binary_val_flat = np.broadcast_to( @@ -521,7 +522,8 @@ def constraints_to_file( slice_size: int = 2_000_000, explicit_coordinate_names: bool = False, ) -> None: - if not len(m.constraints): + regular = m.constraints.regular + if not len(regular): return print_variable, print_constraint = get_printers( @@ -529,10 +531,10 @@ def constraints_to_file( ) f.write(b"\n\ns.t.\n\n") - names = m.constraints + names = list(regular) if progress: names = tqdm( - list(names), + names, desc="Writing constraints.", colour=TQDM_COLOR, ) @@ -540,7 +542,7 @@ def constraints_to_file( # to make this even faster, we can use polars expression # https://docs.pola.rs/user-guide/expressions/plugins/#output-data-types for name in names: - con = m.constraints[name] + con = regular[name] for con_slice in con.iterate_slices(slice_size): df = con_slice.to_polars() @@ -1151,7 +1153,7 @@ def copy(m: Model, include_solution: bool = False, deep: bool = True) -> Model: Model A deep or shallow copy of the model. """ - from linopy.constraints import Constraint, Constraints + from linopy.constraints import Constraint, ConstraintBase, Constraints from linopy.expressions import LinearExpression from linopy.model import Model, Objective from linopy.variables import Variable, Variables @@ -1181,15 +1183,18 @@ def copy(m: Model, include_solution: bool = False, deep: bool = True) -> Model: new_model, ) + def _copy_con_data(con: ConstraintBase) -> xr.Dataset: + d = con.mutable().data + if include_solution: + return d.copy(deep=deep) + cols = m.constraints.dataset_attrs + ( + ["binary_var", "binary_val"] if con.is_indicator else [] + ) + return d[cols].copy(deep=deep) + new_model._constraints = Constraints( { - name: Constraint( - con.mutable().data.copy(deep=deep) - if include_solution - else con.mutable().data[m.constraints.dataset_attrs].copy(deep=deep), - new_model, - name, - ) + name: Constraint(_copy_con_data(con), new_model, name) for name, con in m.constraints.items() }, new_model, From 5e36b19bc4fb49a410ad7195fe051266fb249d22 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 2 Jun 2026 10:48:46 +0200 Subject: [PATCH 05/12] feat(indicator): consume MatrixAccessor indicator arrays in Gurobi (solvers) --- linopy/solvers.py | 46 +++++++++++++++------------------------------- 1 file changed, 15 insertions(+), 31 deletions(-) diff --git a/linopy/solvers.py b/linopy/solvers.py index 8d9321b9..4d627c09 100644 --- a/linopy/solvers.py +++ b/linopy/solvers.py @@ -1601,8 +1601,7 @@ def _build_solver_model( if model.objective.sense == "max": gm.ModelSense = -1 - if len(model.constraints): - assert M.A is not None + if M.A is not None: c = gm.addMConstr(M.A, x, M.sense, M.b) if set_names: names = print_constraints(M.clabels) @@ -1611,41 +1610,26 @@ def _build_solver_model( for sos_type, positions, weights in _iter_sos_sets(model): gm.addSOS(sos_type, x[positions.tolist()].tolist(), weights.tolist()) - if model.indicator_constraints: + if M.indicator_A is not None: sense_map = { "<=": gurobipy.GRB.LESS_EQUAL, ">=": gurobipy.GRB.GREATER_EQUAL, "=": gurobipy.GRB.EQUAL, } x_list = x.tolist() - for ic_data in model.indicator_constraints.values(): - labels_flat = ic_data.labels.values.flatten() - binary_var_flat = ic_data.binary_var.values.flatten() - binary_val_flat = np.broadcast_to( - ic_data.binary_val.values, labels_flat.shape - ).flatten() - coeffs_flat = ic_data.coeffs.values.reshape(len(labels_flat), -1) - vars_flat = ic_data.vars.values.reshape(len(labels_flat), -1) - sign_flat = np.broadcast_to( - ic_data.sign.values, labels_flat.shape - ).flatten() - rhs_flat = np.broadcast_to( - ic_data.rhs.values, labels_flat.shape - ).flatten() - for i in range(len(labels_flat)): - if labels_flat[i] == -1: - continue - lhs = gurobipy.LinExpr() - for c, v in zip(coeffs_flat[i], vars_flat[i]): - if v != -1: - lhs.add(x_list[v], float(c)) - gm.addGenConstrIndicator( - x_list[int(binary_var_flat[i])], - bool(binary_val_flat[i]), - lhs, - sense_map[str(sign_flat[i])], - float(rhs_flat[i]), - ) + A = M.indicator_A + for i in range(A.shape[0]): + lhs = gurobipy.LinExpr() + start, end = A.indptr[i], A.indptr[i + 1] + for col, coeff in zip(A.indices[start:end], A.data[start:end]): + lhs.add(x_list[int(col)], float(coeff)) + gm.addGenConstrIndicator( + x_list[int(M.indicator_binvar[i])], + bool(M.indicator_binval[i]), + lhs, + sense_map[str(M.indicator_sense[i])], + float(M.indicator_b[i]), + ) gm.update() return gm From 36c18639ed97571fbadbd582b873c5b4fad88c2f Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 2 Jun 2026 10:49:24 +0200 Subject: [PATCH 06/12] test(indicator): repair assertions for unified container semantics --- test/test_indicator_constraints.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/test/test_indicator_constraints.py b/test/test_indicator_constraints.py index 2b404168..0cecf7e1 100644 --- a/test/test_indicator_constraints.py +++ b/test/test_indicator_constraints.py @@ -74,7 +74,10 @@ def test_indicator_constraints_not_in_regular_constraints() -> None: m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") assert "regular" in m.constraints - assert "ic0" not in m.constraints + assert "regular" in m.constraints.regular + assert "ic0" in m.constraints + assert "ic0" in m.constraints.indicator + assert "ic0" not in m.constraints.regular assert "ic0" in m.indicator_constraints @@ -102,7 +105,8 @@ def test_indicator_constraints_lp_file(tmp_path: Path) -> None: m.to_file(fn) content = fn.read_text() assert "= 1 ->" in content - assert "ic0:" in content + label = int(m.indicator_constraints["ic0"].labels.item()) + assert f"ic{label}:" in content @pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") From 33f5a24bf47364068574238294cfa9010675fb96 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 2 Jun 2026 10:50:23 +0200 Subject: [PATCH 07/12] fix(indicator): satisfy mypy on return type and binval cast --- linopy/matrices.py | 2 +- linopy/model.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/linopy/matrices.py b/linopy/matrices.py index bc1d9328..c4202731 100644 --- a/linopy/matrices.py +++ b/linopy/matrices.py @@ -84,7 +84,7 @@ def _build_cons(self) -> None: ind_b.append(b) ind_sense.append(sense) ind_binvar.append(label_to_pos[cc._binvar_labels]) - binval = cc._binval + binval = cast("int | np.ndarray", cc._binval) n = len(b) if np.ndim(binval) == 0: ind_binval.append(np.full(n, int(binval), dtype=np.intp)) diff --git a/linopy/model.py b/linopy/model.py index e086b99d..7f45f6a1 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -1108,7 +1108,7 @@ def add_indicator_constraints( sign: SignLike | None = None, rhs: ConstantLike | None = None, name: str | None = None, - ) -> Constraint: + ) -> ConstraintBase: """ Add indicator constraints to the model. @@ -1138,7 +1138,7 @@ def add_indicator_constraints( Returns ------- - linopy.Constraint + linopy.constraints.ConstraintBase The added indicator constraint. """ if not binary_var.attrs.get("binary", False): From f468c1e6087e6ffdce635ee0dc66c95363a227af Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 2 Jun 2026 10:52:59 +0200 Subject: [PATCH 08/12] fix(gurobi): key indicator sense_map on single-char senses MatrixAccessor.indicator_sense holds single-char senses (<, >, =) from to_matrix_with_rhs; the direct-API sense_map was keyed on two-char strings, raising KeyError on the direct path (LP-file path was fine). --- linopy/solvers.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/linopy/solvers.py b/linopy/solvers.py index 4d627c09..fbd7815e 100644 --- a/linopy/solvers.py +++ b/linopy/solvers.py @@ -1612,8 +1612,8 @@ def _build_solver_model( if M.indicator_A is not None: sense_map = { - "<=": gurobipy.GRB.LESS_EQUAL, - ">=": gurobipy.GRB.GREATER_EQUAL, + "<": gurobipy.GRB.LESS_EQUAL, + ">": gurobipy.GRB.GREATER_EQUAL, "=": gurobipy.GRB.EQUAL, } x_list = x.tolist() From 0552c4ef13fca96126938aec87ca1f182157f487 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 2 Jun 2026 10:55:20 +0200 Subject: [PATCH 09/12] Phase C: tests for unified indicator constraints Add freeze round-trip, netCDF round-trip (mutable + frozen), matrix split, copy preservation, and direct-API Gurobi solve coverage. --- test/test_indicator_constraints.py | 117 +++++++++++++++++++++++++++++ 1 file changed, 117 insertions(+) diff --git a/test/test_indicator_constraints.py b/test/test_indicator_constraints.py index 0cecf7e1..fad586c4 100644 --- a/test/test_indicator_constraints.py +++ b/test/test_indicator_constraints.py @@ -6,7 +6,9 @@ import pandas as pd import pytest +import linopy from linopy import Model, available_solvers +from linopy.constraints import Constraint, CSRConstraint def test_add_indicator_constraints_basic() -> None: @@ -207,3 +209,118 @@ def test_indicator_constraints_with_coords() -> None: ic = m.indicator_constraints["ic0"] # Should have 3 indicator constraints (one per index) assert ic.labels.size == 3 + + +def test_indicator_constraint_is_in_unified_container() -> None: + """Indicator constraint lives in the unified container, not in regular.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_constraints(x >= 0, name="regular") + ic = m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + assert isinstance(ic, Constraint) + assert ic.is_indicator + assert "ic0" in m.constraints + assert "ic0" in m.constraints.indicator + assert "ic0" not in m.constraints.regular + assert "ic0" in m.indicator_constraints + assert "regular" in m.constraints.regular + assert "regular" not in m.constraints.indicator + + +def test_indicator_constraint_freeze_roundtrip() -> None: + """is_indicator and binary fields survive freeze and freeze->mutable.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + frozen = m.constraints["ic0"].freeze() + assert isinstance(frozen, CSRConstraint) + assert frozen.is_indicator + assert frozen.binary_var is not None + assert frozen.binary_val == 1 + + mutable = frozen.mutable() + assert isinstance(mutable, Constraint) + assert mutable.is_indicator + assert mutable.binary_var is not None + assert np.all(mutable.binary_val == 1) + + +@pytest.mark.parametrize("freeze_constraints", [False, True]) +def test_indicator_constraint_netcdf_roundtrip( + tmp_path: Path, freeze_constraints: bool +) -> None: + """is_indicator and binary fields survive a netCDF round-trip.""" + m = Model(freeze_constraints=freeze_constraints) + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_constraints(x >= 0, name="regular") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + fn = tmp_path / "model.nc" + m.to_netcdf(fn) + m2 = linopy.read_netcdf(fn) + + ic = m2.constraints["ic0"] + assert ic.is_indicator + assert ic.binary_var is not None + assert np.all(ic.binary_val == 1) + assert "regular" not in m2.constraints.indicator + + +def test_indicator_constraint_matrix_split() -> None: + """Regular A excludes indicator rows; indicator arrays carry them.""" + m = Model() + idx = pd.RangeIndex(3, name="i") + b = m.add_variables(coords=[idx], name="b", binary=True) + x = m.add_variables(coords=[idx], lower=0, upper=10, name="x") + m.add_constraints(x >= 0, name="regular") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + n = idx.size + assert m.matrices.A.shape[0] == n + assert m.matrices.indicator_A.shape[0] == n + assert len(m.matrices.clabels) == n + + np.testing.assert_array_equal(m.matrices.indicator_binval, np.full(n, 1)) + np.testing.assert_array_equal(m.matrices.indicator_b, np.full(n, 5.0)) + np.testing.assert_array_equal(m.matrices.indicator_sense, np.full(n, "<")) + assert m.matrices.indicator_binvar.shape == (n,) + + +def test_indicator_constraint_copy_preserves() -> None: + """Model.copy preserves is_indicator and the binary fields.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + mc = m.copy() + ic = mc.constraints["ic0"] + assert ic.is_indicator + assert ic.binary_var is not None + assert np.all(ic.binary_val == 1) + + +@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") +@pytest.mark.parametrize( + ("trigger", "expected"), + [(1, 5), (0, 10)], + ids=["active", "inactive"], +) +def test_indicator_constraints_solve_gurobi_direct( + trigger: int, expected: float +) -> None: + """Direct-API Gurobi solve enforces the indicator only at the trigger.""" + m = Model() + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_constraints(b >= 1, name="fix_b") + m.add_indicator_constraints(b, trigger, x, "<=", 5, name="ic0") + m.add_objective(x, sense="max") + m.solve(solver_name="gurobi", io_api="direct") + assert m.objective.value is not None + assert np.isclose(m.objective.value, expected, atol=1e-6) From 2f036d4c89a1b657f6eb2c5de765a00d69f73183 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 2 Jun 2026 11:05:09 +0200 Subject: [PATCH 10/12] =?UTF-8?q?Exclude=20indicators=20from=20Constraints?= =?UTF-8?q?.ncons=20and=20to=5Fmatrix=20(plan=20=C2=A74)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ncons feeds model.ncons which the Mosek backend uses to size the (regular-only) matrix M.A; keep them aligned. Also exclude indicators from the consolidated to_matrix. Add ncons assertion and release note. --- doc/release_notes.rst | 2 +- linopy/constraints.py | 6 +++++- test/test_indicator_constraints.py | 1 + 3 files changed, 7 insertions(+), 2 deletions(-) diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 61472f13..f741e57b 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -5,7 +5,7 @@ Upcoming Version ---------------- * Add documentation about `LinearExpression.where` with `drop=True`. Add `BaseExpression.variable_names` property. -* Add indicator constraints for solvers that support them. +* Add indicator constraints for solvers that support them. They are part of the unified constraints container: ``model.add_indicator_constraints`` returns a ``Constraint`` and the constraint is stored in ``model.constraints`` (filterable via ``model.constraints.indicator`` / ``model.constraints.regular``), so it round-trips through netCDF and ``model.copy()``. **Features** diff --git a/linopy/constraints.py b/linopy/constraints.py index 952437d3..2b76a072 100644 --- a/linopy/constraints.py +++ b/linopy/constraints.py @@ -1754,7 +1754,7 @@ def ncons(self) -> int: This excludes constraints with missing labels or where all variables are masked (vars == -1). """ - return sum(con.ncons for con in self.data.values()) + return sum(con.ncons for con in self.data.values() if not con.is_indicator) @property def inequalities(self) -> Constraints: @@ -1972,9 +1972,13 @@ def to_matrix(self) -> tuple[scipy.sparse.csr_array, np.ndarray]: csrs = [] con_labels_list = [] for c in self.data.values(): + if c.is_indicator: + continue csr, con_labels = c.to_matrix(label_index) csrs.append(csr) con_labels_list.append(con_labels) + if not csrs: + raise ValueError("No constraints available to convert to matrix.") return scipy.sparse.vstack(csrs, format="csr"), np.concatenate(con_labels_list) def reset_dual(self) -> None: diff --git a/test/test_indicator_constraints.py b/test/test_indicator_constraints.py index fad586c4..7ec67382 100644 --- a/test/test_indicator_constraints.py +++ b/test/test_indicator_constraints.py @@ -284,6 +284,7 @@ def test_indicator_constraint_matrix_split() -> None: assert m.matrices.A.shape[0] == n assert m.matrices.indicator_A.shape[0] == n assert len(m.matrices.clabels) == n + assert m.ncons == n np.testing.assert_array_equal(m.matrices.indicator_binval, np.full(n, 1)) np.testing.assert_array_equal(m.matrices.indicator_b, np.full(n, 5.0)) From 170c693029084b2a264bfe960de8c06a955f0a7a Mon Sep 17 00:00:00 2001 From: Fabian Date: Wed, 3 Jun 2026 10:24:42 +0200 Subject: [PATCH 11/12] fix mypy in new tests --- test/test_indicator_constraints.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_indicator_constraints.py b/test/test_indicator_constraints.py index 7ec67382..e5348ff6 100644 --- a/test/test_indicator_constraints.py +++ b/test/test_indicator_constraints.py @@ -281,6 +281,8 @@ def test_indicator_constraint_matrix_split() -> None: m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") n = idx.size + assert m.matrices.A is not None + assert m.matrices.indicator_A is not None assert m.matrices.A.shape[0] == n assert m.matrices.indicator_A.shape[0] == n assert len(m.matrices.clabels) == n From b320c4ffce590a9f12fbfd5a59ab0236009213d4 Mon Sep 17 00:00:00 2001 From: Fabian Date: Wed, 3 Jun 2026 12:37:58 +0200 Subject: [PATCH 12/12] refactor tests --- test/test_indicator_constraints.py | 708 ++++++++++++++++------------- 1 file changed, 402 insertions(+), 306 deletions(-) diff --git a/test/test_indicator_constraints.py b/test/test_indicator_constraints.py index e5348ff6..bba6498a 100644 --- a/test/test_indicator_constraints.py +++ b/test/test_indicator_constraints.py @@ -9,321 +9,417 @@ import linopy from linopy import Model, available_solvers from linopy.constraints import Constraint, CSRConstraint +from linopy.variables import Variable - -def test_add_indicator_constraints_basic() -> None: - """Indicator constraint is created with correct fields.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - - assert "ic0" in m.indicator_constraints - ic = m.indicator_constraints["ic0"] - assert "coeffs" in ic - assert "vars" in ic - assert "sign" in ic - assert "rhs" in ic - assert "binary_var" in ic - assert "binary_val" in ic - assert "labels" in ic - - -def test_add_indicator_constraints_from_constraint() -> None: - """Indicator constraint accepts a Constraint object as lhs.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - con = x <= 5 - m.add_indicator_constraints(b, 1, con, name="ic0") - assert "ic0" in m.indicator_constraints - - -def test_indicator_constraints_validation_non_binary() -> None: - """Non-binary variable raises ValueError.""" - m = Model() - y = m.add_variables(lower=0, upper=1, name="y") - x = m.add_variables(lower=0, upper=10, name="x") - with pytest.raises(ValueError, match="must be binary"): - m.add_indicator_constraints(y, 1, x, "<=", 5) - - -def test_indicator_constraints_validation_bad_value() -> None: - """Invalid binary_val raises ValueError.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - with pytest.raises(ValueError, match="must be 0 or 1"): - m.add_indicator_constraints(b, 2, x, "<=", 5) - - -def test_indicator_constraints_duplicate_name() -> None: - """Duplicate name raises ValueError.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - with pytest.raises(ValueError, match="already assigned"): - m.add_indicator_constraints(b, 0, x, ">=", 0, name="ic0") - - -def test_indicator_constraints_not_in_regular_constraints() -> None: - """Indicator constraints are separate from regular constraints.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - m.add_constraints(x >= 0, name="regular") - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - - assert "regular" in m.constraints - assert "regular" in m.constraints.regular - assert "ic0" in m.constraints - assert "ic0" in m.constraints.indicator - assert "ic0" not in m.constraints.regular - assert "ic0" in m.indicator_constraints - - -def test_remove_indicator_constraints() -> None: - """Indicator constraints can be removed.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - assert "ic0" in m.indicator_constraints - m.remove_indicator_constraints("ic0") - assert "ic0" not in m.indicator_constraints - - -def test_indicator_constraints_lp_file(tmp_path: Path) -> None: - """LP file contains general constraints section.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - m.add_constraints(x >= 0, name="dummy") - m.add_objective(x) - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - - fn = tmp_path / "test.lp" - m.to_file(fn) - content = fn.read_text() - assert "= 1 ->" in content - label = int(m.indicator_constraints["ic0"].labels.item()) - assert f"ic{label}:" in content - - -@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") -def test_indicator_constraints_solve_gurobi_active() -> None: - """ - Indicator constraint enforced when binary is at trigger value. - - Maximize x subject to: - b = 1 (fixed) - (b == 1) => (x <= 5) - x in [0, 10] - Since b=1 and the indicator enforces x<=5, optimal x=5. - """ - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - # Force b = 1 - m.add_constraints(b >= 1, name="fix_b") - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - m.add_objective(x, sense="max") - m.solve(solver_name="gurobi") - assert m.objective.value is not None - assert np.isclose(m.objective.value, 5, atol=1e-6) - - -@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") -def test_indicator_constraints_solve_gurobi_inactive() -> None: - """ - Indicator constraint NOT enforced when binary differs from trigger. - - Maximize x subject to: - b = 1 (fixed) - (b == 0) => (x <= 5) - x in [0, 10] - Since b=1 but trigger is 0, the constraint is inactive. Optimal x=10. - """ - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - # Force b = 1 - m.add_constraints(b >= 1, name="fix_b") - m.add_indicator_constraints(b, 0, x, "<=", 5, name="ic0") - m.add_objective(x, sense="max") - m.solve(solver_name="gurobi") - assert m.objective.value is not None - assert np.isclose(m.objective.value, 10, atol=1e-6) - - -@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") -def test_indicator_constraints_solve_gurobi_multiple() -> None: - """ - Multiple indicator constraints work together. - - b1=1, b2=1 (forced), x in [0, 20]. - (b1 == 1) => (x <= 10) - (b2 == 1) => (x <= 5) - Maximize x. - Both constraints active, so x <= min(10, 5) = 5. Optimal x=5. - """ - m = Model() - b1 = m.add_variables(name="b1", binary=True) - b2 = m.add_variables(name="b2", binary=True) - x = m.add_variables(lower=0, upper=20, name="x") - m.add_constraints(b1 >= 1, name="fix_b1") - m.add_constraints(b2 >= 1, name="fix_b2") - m.add_indicator_constraints(b1, 1, x, "<=", 10, name="ic1") - m.add_indicator_constraints(b2, 1, x, "<=", 5, name="ic2") - m.add_objective(x, sense="max") - m.solve(solver_name="gurobi") - assert m.objective.value is not None - assert np.isclose(m.objective.value, 5, atol=1e-6) - - -def test_unsupported_solver_raises() -> None: - """Solvers without indicator support raise ValueError.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - m.add_constraints(x >= 0, name="dummy") - m.add_objective(x) - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - - for solver in ["glpk", "highs", "mosek", "mindopt"]: - if solver in available_solvers: - with pytest.raises( - ValueError, match="does not support indicator constraints" - ): - m.solve(solver_name=solver) - - -def test_indicator_constraints_with_coords() -> None: - """Indicator constraints work with multi-dimensional coords.""" - m = Model() - idx = pd.RangeIndex(3, name="i") - b = m.add_variables(coords=[idx], name="b", binary=True) - x = m.add_variables(coords=[idx], lower=0, upper=10, name="x") - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - ic = m.indicator_constraints["ic0"] - # Should have 3 indicator constraints (one per index) - assert ic.labels.size == 3 - - -def test_indicator_constraint_is_in_unified_container() -> None: - """Indicator constraint lives in the unified container, not in regular.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - m.add_constraints(x >= 0, name="regular") - ic = m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - - assert isinstance(ic, Constraint) - assert ic.is_indicator - assert "ic0" in m.constraints - assert "ic0" in m.constraints.indicator - assert "ic0" not in m.constraints.regular - assert "ic0" in m.indicator_constraints - assert "regular" in m.constraints.regular - assert "regular" not in m.constraints.indicator +requires_gurobi = pytest.mark.skipif( + "gurobi" not in available_solvers, reason="Gurobi not installed" +) -def test_indicator_constraint_freeze_roundtrip() -> None: - """is_indicator and binary fields survive freeze and freeze->mutable.""" +@pytest.fixture +def mbx() -> tuple[Model, Variable, Variable]: + """Model with a scalar binary ``b`` and continuous ``x`` in [0, 10].""" m = Model() b = m.add_variables(name="b", binary=True) x = m.add_variables(lower=0, upper=10, name="x") - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - - frozen = m.constraints["ic0"].freeze() - assert isinstance(frozen, CSRConstraint) - assert frozen.is_indicator - assert frozen.binary_var is not None - assert frozen.binary_val == 1 - - mutable = frozen.mutable() - assert isinstance(mutable, Constraint) - assert mutable.is_indicator - assert mutable.binary_var is not None - assert np.all(mutable.binary_val == 1) - - -@pytest.mark.parametrize("freeze_constraints", [False, True]) -def test_indicator_constraint_netcdf_roundtrip( - tmp_path: Path, freeze_constraints: bool -) -> None: - """is_indicator and binary fields survive a netCDF round-trip.""" - m = Model(freeze_constraints=freeze_constraints) - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - m.add_constraints(x >= 0, name="regular") - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - - fn = tmp_path / "model.nc" - m.to_netcdf(fn) - m2 = linopy.read_netcdf(fn) + return m, b, x - ic = m2.constraints["ic0"] - assert ic.is_indicator - assert ic.binary_var is not None - assert np.all(ic.binary_val == 1) - assert "regular" not in m2.constraints.indicator - -def test_indicator_constraint_matrix_split() -> None: - """Regular A excludes indicator rows; indicator arrays carry them.""" +@pytest.fixture +def coords_mbx() -> tuple[Model, Variable, Variable, pd.RangeIndex]: + """Model with index-aligned binary ``b`` and continuous ``x`` over a length-3 index.""" m = Model() idx = pd.RangeIndex(3, name="i") b = m.add_variables(coords=[idx], name="b", binary=True) x = m.add_variables(coords=[idx], lower=0, upper=10, name="x") - m.add_constraints(x >= 0, name="regular") - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - - n = idx.size - assert m.matrices.A is not None - assert m.matrices.indicator_A is not None - assert m.matrices.A.shape[0] == n - assert m.matrices.indicator_A.shape[0] == n - assert len(m.matrices.clabels) == n - assert m.ncons == n - - np.testing.assert_array_equal(m.matrices.indicator_binval, np.full(n, 1)) - np.testing.assert_array_equal(m.matrices.indicator_b, np.full(n, 5.0)) - np.testing.assert_array_equal(m.matrices.indicator_sense, np.full(n, "<")) - assert m.matrices.indicator_binvar.shape == (n,) - - -def test_indicator_constraint_copy_preserves() -> None: - """Model.copy preserves is_indicator and the binary fields.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") - - mc = m.copy() - ic = mc.constraints["ic0"] - assert ic.is_indicator - assert ic.binary_var is not None - assert np.all(ic.binary_val == 1) - - -@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobi not installed") -@pytest.mark.parametrize( - ("trigger", "expected"), - [(1, 5), (0, 10)], - ids=["active", "inactive"], -) -def test_indicator_constraints_solve_gurobi_direct( - trigger: int, expected: float -) -> None: - """Direct-API Gurobi solve enforces the indicator only at the trigger.""" - m = Model() - b = m.add_variables(name="b", binary=True) - x = m.add_variables(lower=0, upper=10, name="x") - m.add_constraints(b >= 1, name="fix_b") - m.add_indicator_constraints(b, trigger, x, "<=", 5, name="ic0") - m.add_objective(x, sense="max") - m.solve(solver_name="gurobi", io_api="direct") - assert m.objective.value is not None - assert np.isclose(m.objective.value, expected, atol=1e-6) + return m, b, x, idx + + +class TestConstruction: + """Creating indicator constraints and validating their arguments.""" + + def test_basic_fields(self, mbx: tuple[Model, Variable, Variable]) -> None: + """Indicator constraint is created with correct fields.""" + m, b, x = mbx + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + assert "ic0" in m.indicator_constraints + ic = m.indicator_constraints["ic0"] + for field in ( + "coeffs", + "vars", + "sign", + "rhs", + "binary_var", + "binary_val", + "labels", + ): + assert field in ic + + def test_auto_name(self, mbx: tuple[Model, Variable, Variable]) -> None: + """Omitting ``name`` auto-generates an ``indcon`` name.""" + m, b, x = mbx + con = m.add_indicator_constraints(b, 1, x, "<=", 5) + assert con.name.startswith("indcon") + assert con.name in m.indicator_constraints + + @pytest.mark.parametrize( + ("build_lhs", "kwargs"), + [ + (lambda x: x <= 5, {}), + (lambda x: 2 * x, {"sign": "<=", "rhs": 5}), + (lambda x: 1 * x.at[()] <= 5, {}), + ], + ids=["constraint", "linexpr", "scalarcon"], + ) + def test_valid_lhs( + self, mbx: tuple[Model, Variable, Variable], build_lhs: object, kwargs: dict + ) -> None: + """Each accepted lhs form (constraint, expression, scalar) builds an indicator.""" + m, b, x = mbx + ic = m.add_indicator_constraints(b, 1, build_lhs(x), name="ic0", **kwargs) # type: ignore[operator] + assert ic.is_indicator + assert "ic0" in m.indicator_constraints + + @pytest.mark.parametrize( + ("build_lhs", "kwargs", "exc", "match"), + [ + (lambda x: x <= 5, {"sign": "<=", "rhs": 5}, ValueError, "must be None"), + ( + lambda x: 2 * x, + {}, + ValueError, + "required when lhs is a LinearExpression", + ), + (lambda x: x, {}, ValueError, "required when lhs is a Variable"), + ( + lambda x: 1 * x.at[()] <= 5, + {"sign": "<=", "rhs": 5}, + ValueError, + "must be None", + ), + (lambda x: "not-a-constraint", {}, TypeError, "must be a Constraint"), + ], + ids=[ + "con+sign", + "linexpr-no-sign", + "var-no-sign", + "scalarcon+sign", + "bad-type", + ], + ) + def test_invalid_lhs( + self, + mbx: tuple[Model, Variable, Variable], + build_lhs: object, + kwargs: dict, + exc: type[Exception], + match: str, + ) -> None: + """Invalid lhs / sign / rhs combinations raise with a clear message.""" + m, b, x = mbx + with pytest.raises(exc, match=match): + m.add_indicator_constraints(b, 1, build_lhs(x), name="ic0", **kwargs) # type: ignore[operator] + + def test_non_binary_var_raises(self, mbx: tuple[Model, Variable, Variable]) -> None: + """Non-binary variable raises ValueError.""" + m, _, x = mbx + y = m.add_variables(lower=0, upper=1, name="y") + with pytest.raises(ValueError, match="must be binary"): + m.add_indicator_constraints(y, 1, x, "<=", 5) + + def test_bad_binary_val_raises(self, mbx: tuple[Model, Variable, Variable]) -> None: + """Invalid binary_val raises ValueError.""" + m, b, x = mbx + with pytest.raises(ValueError, match="must be 0 or 1"): + m.add_indicator_constraints(b, 2, x, "<=", 5) + + def test_duplicate_name_raises(self, mbx: tuple[Model, Variable, Variable]) -> None: + """Duplicate name raises ValueError.""" + m, b, x = mbx + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + with pytest.raises(ValueError, match="already assigned"): + m.add_indicator_constraints(b, 0, x, ">=", 0, name="ic0") + + +class TestContainer: + """How indicator constraints live in (and stay separate within) the container.""" + + def test_separate_from_regular(self, mbx: tuple[Model, Variable, Variable]) -> None: + """Indicator constraints are separate from regular constraints.""" + m, b, x = mbx + m.add_constraints(x >= 0, name="regular") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + assert "regular" in m.constraints + assert "regular" in m.constraints.regular + assert "ic0" in m.constraints + assert "ic0" in m.constraints.indicator + assert "ic0" not in m.constraints.regular + assert "ic0" in m.indicator_constraints + + def test_in_unified_container(self, mbx: tuple[Model, Variable, Variable]) -> None: + """Indicator constraint lives in the unified container, not in regular.""" + m, b, x = mbx + m.add_constraints(x >= 0, name="regular") + ic = m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + assert isinstance(ic, Constraint) + assert ic.is_indicator + assert "ic0" in m.constraints + assert "ic0" in m.constraints.indicator + assert "ic0" not in m.constraints.regular + assert "ic0" in m.indicator_constraints + assert "regular" in m.constraints.regular + assert "regular" not in m.constraints.indicator + + def test_remove(self, mbx: tuple[Model, Variable, Variable]) -> None: + """Indicator constraints can be removed.""" + m, b, x = mbx + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + assert "ic0" in m.indicator_constraints + m.remove_indicator_constraints("ic0") + assert "ic0" not in m.indicator_constraints + + def test_regular_constraint_has_no_indicator_fields( + self, mbx: tuple[Model, Variable, Variable] + ) -> None: + """Regular constraints report no indicator metadata, frozen or not.""" + m, _, x = mbx + m.add_constraints(x <= 5, name="c") + con = m.constraints["c"] + assert not con.is_indicator + assert con.binary_var is None + assert con.binary_val is None + + frozen = con.freeze() + assert not frozen.is_indicator + assert frozen.binary_var is None + assert frozen.binary_val is None + + def test_with_coords( + self, coords_mbx: tuple[Model, Variable, Variable, pd.RangeIndex] + ) -> None: + """Indicator constraints work with multi-dimensional coords.""" + m, b, x, idx = coords_mbx + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + assert m.indicator_constraints["ic0"].labels.size == idx.size + + +class TestPersistence: + """Indicator metadata survives freeze, copy, and netCDF round-trips.""" + + def test_freeze_roundtrip(self, mbx: tuple[Model, Variable, Variable]) -> None: + """is_indicator and binary fields survive freeze and freeze->mutable.""" + m, b, x = mbx + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + frozen = m.constraints["ic0"].freeze() + assert isinstance(frozen, CSRConstraint) + assert frozen.is_indicator + assert frozen.binary_var is not None + assert frozen.binary_val == 1 + + mutable = frozen.mutable() + assert isinstance(mutable, Constraint) + assert mutable.is_indicator + assert mutable.binary_var is not None + assert np.all(mutable.binary_val == 1) + + def test_copy_preserves(self, mbx: tuple[Model, Variable, Variable]) -> None: + """Model.copy preserves is_indicator and the binary fields.""" + m, b, x = mbx + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + ic = m.copy().constraints["ic0"] + assert ic.is_indicator + assert ic.binary_var is not None + assert np.all(ic.binary_val == 1) + + @pytest.mark.parametrize("freeze_constraints", [False, True]) + def test_netcdf_roundtrip(self, tmp_path: Path, freeze_constraints: bool) -> None: + """is_indicator and binary fields survive a netCDF round-trip.""" + m = Model(freeze_constraints=freeze_constraints) + b = m.add_variables(name="b", binary=True) + x = m.add_variables(lower=0, upper=10, name="x") + m.add_constraints(x >= 0, name="regular") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + fn = tmp_path / "model.nc" + m.to_netcdf(fn) + ic = linopy.read_netcdf(fn).constraints["ic0"] + + assert ic.is_indicator + assert ic.binary_var is not None + assert np.all(ic.binary_val == 1) + + def test_array_binval_roundtrip(self, tmp_path: Path) -> None: + """A coords-based indicator has per-element binary_val that round-trips.""" + m = Model(freeze_constraints=True) + idx = pd.RangeIndex(3, name="i") + b = m.add_variables(coords=[idx], name="b", binary=True) + x = m.add_variables(coords=[idx], lower=0, upper=10, name="x") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + frozen = m.constraints["ic0"] + assert isinstance(frozen, CSRConstraint) + assert isinstance(frozen.binary_val, np.ndarray) and frozen.binary_val.ndim == 1 + assert "binary_val" in frozen.data + + fn = tmp_path / "model.nc" + m.to_netcdf(fn) + ic = linopy.read_netcdf(fn).constraints["ic0"] + assert ic.is_indicator + assert np.all(ic.binary_val == 1) + + +class TestMatrices: + """Indicator rows are split out of the regular constraint matrix.""" + + def test_matrix_split( + self, coords_mbx: tuple[Model, Variable, Variable, pd.RangeIndex] + ) -> None: + """Regular A excludes indicator rows; indicator arrays carry them.""" + m, b, x, idx = coords_mbx + m.add_constraints(x >= 0, name="regular") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + n = idx.size + assert m.matrices.A is not None + assert m.matrices.indicator_A is not None + assert m.matrices.A.shape[0] == n + assert m.matrices.indicator_A.shape[0] == n + assert len(m.matrices.clabels) == n + assert m.ncons == n + + np.testing.assert_array_equal(m.matrices.indicator_binval, np.full(n, 1)) + np.testing.assert_array_equal(m.matrices.indicator_b, np.full(n, 5.0)) + np.testing.assert_array_equal(m.matrices.indicator_sense, np.full(n, "<")) + assert m.matrices.indicator_binvar.shape == (n,) + + def test_to_matrix_skips_indicator( + self, coords_mbx: tuple[Model, Variable, Variable, pd.RangeIndex] + ) -> None: + """Constraints.to_matrix drops indicator rows; an indicator-only set raises.""" + m, b, x, idx = coords_mbx + m.add_constraints(x <= 5, name="regular") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + A, labels = m.constraints.to_matrix() + assert A.shape[0] == idx.size + + with pytest.raises(ValueError, match="No constraints available"): + m.constraints.indicator.to_matrix() + + +class TestLPFile: + """LP export of indicator constraints.""" + + def test_with_regular_constraints( + self, mbx: tuple[Model, Variable, Variable], tmp_path: Path + ) -> None: + """LP file contains general constraints section.""" + m, b, x = mbx + m.add_constraints(x >= 0, name="dummy") + m.add_objective(x) + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + fn = tmp_path / "test.lp" + m.to_file(fn) + content = fn.read_text() + assert "= 1 ->" in content + label = int(m.indicator_constraints["ic0"].labels.item()) + assert f"ic{label}:" in content + + def test_indicator_only( + self, mbx: tuple[Model, Variable, Variable], tmp_path: Path + ) -> None: + """An LP file with no regular constraints still writes the indicator section.""" + m, b, x = mbx + m.add_objective(x) + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + fn = tmp_path / "test.lp" + m.to_file(fn) + content = fn.read_text() + assert "s.t." in content + assert "= 1 ->" in content + + +class TestSolve: + """Solving models with indicator constraints.""" + + @requires_gurobi + @pytest.mark.parametrize( + ("io_api", "trigger", "expected"), + [ + (None, 1, 5), + (None, 0, 10), + ("direct", 1, 5), + ("direct", 0, 10), + ], + ids=["lp-active", "lp-inactive", "direct-active", "direct-inactive"], + ) + def test_gurobi_enforces_at_trigger( + self, + mbx: tuple[Model, Variable, Variable], + io_api: str | None, + trigger: int, + expected: float, + ) -> None: + """ + The indicator is enforced (x<=5) only when b matches the trigger value. + + With b fixed to 1, an active trigger caps x at 5; an inactive one leaves + x free to its upper bound of 10. + """ + m, b, x = mbx + m.add_constraints(b >= 1, name="fix_b") + m.add_indicator_constraints(b, trigger, x, "<=", 5, name="ic0") + m.add_objective(x, sense="max") + m.solve(solver_name="gurobi", io_api=io_api) + assert m.objective.value is not None + assert np.isclose(m.objective.value, expected, atol=1e-6) + + @requires_gurobi + def test_gurobi_multiple(self) -> None: + """Multiple indicators combine: x <= min(10, 5) = 5.""" + m = Model() + b1 = m.add_variables(name="b1", binary=True) + b2 = m.add_variables(name="b2", binary=True) + x = m.add_variables(lower=0, upper=20, name="x") + m.add_constraints(b1 >= 1, name="fix_b1") + m.add_constraints(b2 >= 1, name="fix_b2") + m.add_indicator_constraints(b1, 1, x, "<=", 10, name="ic1") + m.add_indicator_constraints(b2, 1, x, "<=", 5, name="ic2") + m.add_objective(x, sense="max") + m.solve(solver_name="gurobi") + assert m.objective.value is not None + assert np.isclose(m.objective.value, 5, atol=1e-6) + + def test_unsupported_solver_raises( + self, mbx: tuple[Model, Variable, Variable] + ) -> None: + """Solvers without indicator support raise ValueError.""" + m, b, x = mbx + m.add_constraints(x >= 0, name="dummy") + m.add_objective(x) + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + + for solver in ["glpk", "highs", "mosek", "mindopt"]: + if solver in available_solvers: + with pytest.raises( + ValueError, match="does not support indicator constraints" + ): + m.solve(solver_name=solver) + + @requires_gurobi + def test_mip_has_no_duals(self, mbx: tuple[Model, Variable, Variable]) -> None: + """Indicator rows are skipped when collecting duals; MIPs expose none.""" + m, b, x = mbx + m.add_constraints(b >= 1, name="fix_b") + m.add_indicator_constraints(b, 1, x, "<=", 5, name="ic0") + m.add_objective(x, sense="max") + m.solve(solver_name="gurobi") + with pytest.raises(AttributeError, match="dual"): + _ = m.matrices.dual