Using SOS on masked variables breaks

[FBumann]

Summary

Direct-API SOS construction is broken for both solvers that implement it (Gurobi and Xpress, after #684). Both pass linopy variable labels as if they were solver column positions; the two are only equal when no variable in the model is masked. With a single masked entry anywhere, the SOS set is built against the wrong columns — and the solvers happen to raise rather than silently mis-build, so the bug has been invisible until you try a masked SOS variable.

Repro

# dev-scripts/repro_sos_minus1_both.py
import pandas as pd
from linopy import Model


def build():
    m = Model()
    coords = pd.Index([0, 1, 2, 3], name="i")
    mask = pd.Series([True, True, False, True], index=coords)  # one slot masked
    sos_var = m.add_variables(
        lower=0, upper=1, coords=[coords], mask=mask, name="sos_var"
    )
    m.add_sos_constraints(sos_var, sos_type=1, sos_dim="i")
    m.add_objective(-sos_var.sum())
    return m


for solver in ("gurobi", "xpress"):
    print(f"\n=== {solver} ===")
    try:
        build().solve(solver_name=solver, io_api="direct")
    except Exception as e:
        print(f"  RAISED {type(e).__name__}: {e}")

Output (against the xpress-direct branch from #684, on master for Gurobi):

=== gurobi ===
  RAISED IndexError: index 3 is out of bounds for axis 0 with size 3

=== xpress ===
  ?404 Error: Invalid column number passed to XPRSaddsets.
  Element 2 of your array has invalid column number 3
  RAISED SolverError

Root cause

In both _build_solver_model implementations, the SOS code passes var.labels values straight into the solver's addSOS as if they were column indices:

• Gurobi (linopy/solvers.py:1530-1532): gm.addSOS(sos_type, x[indices].tolist(), weights) — x is gm.addMVar(M.vlabels.shape, ...), indexed by position. The code passes linopy labels as positions.
• Xpress (feat(xpress): add direct API support via loadproblem #684, linopy/solvers.py:2168-2170): problem.addSOS(indices, weights, type=sos_type) — problem was built via loadproblem, where column position in the solver matches position in M.vlabels. Same mistake.

Whenever any variable in the model is masked, M.vlabels becomes non-contiguous (e.g. [0, 2, 3] instead of [0, 1, 2, 3]), and linopy label ≠ solver column position. The labels that worked correctly under the no-mask happy path now over-index the solver model.

The Xpress branch additionally filters out -1 labels (Xpress lines 2161-2167), which the Gurobi branch doesn't — that catches one shape of the bug (the masked entry inside the SOS variable itself) but not the broader label-vs-position mismatch caused by masks anywhere in the model.

Fix direction

The SOS preparation is the same logic in both solvers up to the vendor call. Worth pulling it into a single shared step that:

1. drops -1 (missing) entries from var.labels and the matching weights
2. resolves the remaining linopy labels to solver column positions (e.g. via model.variables.label_index)
3. returns (indices, weights) ready for the vendor addSOS

Each solver then keeps just its native call. Where the helper lives is open.

Needs regression coverage with a masked SOS variable in both test_optimization.py (end-to-end solve) and test_sos_constraints.py (model construction).

[FBumann]

I found that this also fails on file based io!

[FBumann]

@FabianHofmann Im digging deeper into the sos stuff, and im seeing some other inconsistencies between Model.solve() and the new Solver.from_model().solve().

Model.solve() vs Solver.from_model(...).solve()

Both paths run the same Solver subclass underneath — Model.solve() calls Solver.from_name() at model.py:1755. The tables below catalog what the high-level wrapper does around the build/run.

Each note is tagged with one of four categories; the categories drive the choice at the bottom.

• shape — whether/how the model can be solved (validation, sanitization)
• transform — mutates the model going in, restores going out
• orchestration — how the user wants the solve invoked
• wiring — copying the result back into the model

Legend: ✅ done for the user · (✅) achievable manually · ❌ missing · n/a not applicable. (M = Model.solve(), S = Solver.from_model().solve().)

Result handling

Feature	M	S	Notes
Primal → model.variables[k].solution	✅	(✅)	wiring — model.assign_result(result) (model.py:1782)
Duals → model.constraints[k].dual	✅	(✅)	wiring — covered by assign_result()
model.objective._value	✅	(✅)	wiring — covered by assign_result()
model.status / termination_condition	✅	(✅)	wiring — covered by assign_result()
model.solver / solver_model / solver_name	✅	(✅)	wiring — caller must model.solver = solver for compute_infeasibilities() / IIS
Close previous solver before re-assign	✅	❌	wiring — lifecycle owned by the Model
Return type	tuple	Result	different ergonomics

Pre-solve validation

Feature	M	S	Notes
Raise on empty objective	✅	❌	shape — clear Python error instead of opaque solver complaint (:1616)
Raise if picked solver not installed	✅ (any-installed check, :1659)	✅ (this-installed check, solvers.py:409)	both check, different scope — Model.solve() asserts ≥1 solver available, Solver.__post_init__ rejects the specific subclass
Default solver_name = available_solvers[0]	✅	n/a	caller picks subclass
Quadratic without QUADRATIC_OBJECTIVE	✅	❌	shape — fail-fast at submit (:1700)
Semi-continuous without SEMI_CONTINUOUS_VARIABLES	✅	❌	shape — same fail-fast (:1729)
SOS without SOS_CONSTRAINTS + no reformulation	✅	❌	shape — coupled to reformulate_sos (:1723)

Model preparation

Feature	M	S	Notes
sanitize_zeros	✅	❌ for io_api ∈ {lp,mps,lp-polars} · (✅) for io_api="direct"	shape — marks zero-coeff terms with -1. LP/MPS writer would otherwise emit literal 0 x1 terms; direct API's _matrix_export_data (+ solver pruning) makes it effectively a no-op (:1694)
sanitize_infinities	✅	❌ for io_api ∈ {lp,mps,lp-polars} · (✅) for io_api="direct"	shape — marks trivial <= inf / >= -inf rows; same split — biggest win is shrinking the LP/MPS file the solver has to parse, direct API ignores them anyway (:1697)
reset_solution() before solve	✅	❌	transform — clears stale solution if next solve fails (:1669)
reformulate_sos (True/False/auto)	✅	❌	transform — lets non-SOS solvers handle SOS via binary + linking (:1713-1727)
Undo SOS reformulation after solve	✅	❌	transform — without undo, synthetic vars/constraints leak (:1779)

Orchestration / convenience

Feature	M	S	Notes
mock_solve=True	✅	❌	orchestration — bypasses Solver entirely (:1611, :1817)
remote=RemoteHandler/OetcHandler	✅	❌	orchestration — ships model elsewhere (:1628-1657)
keep_files auto-cleanup	✅	❌	orchestration — tempfile lifecycle under model.solver_dir (:1772-1774)
Default solution_fn	✅	❌	orchestration — via model.get_solution_file() (:1684-1692)
Default problem_fn	✅	✅	Solver path uses the same helper
Info logs (solver, options, log path)	✅	❌	orchestration — cosmetic (:1665, :1671, :1674-1678)

Two options

1. Keep the Solver path as a lower-level building block

Status quo, formalised. Model.solve() stays the canonical path. Callers on Solver.from_model() compose the exposed helpers (constraints.sanitize_zeros(), reformulate_sos_constraints(), model.assign_result()) themselves. No churn; the asymmetry stays, documented.

2. Push shape + transform + wiring down to Solver, keep orchestration on Model

Map straight off the category tags above:

• shape → new kwargs on Solver.from_model() (built-in safety checks + opt-in sanitize_*). _build() stays private.
• transform → shared context manager (with reformulate_sos_if_needed(model, solver_cls):) reused by both paths.
• wiring → extend model.assign_result() to also stash model.solver = solver; document it on Solver.solve().
• orchestration → stays on Model.solve().

Note: kwarg placement should be clean — since from_model() always runs before solve(), kwargs that today appear on both (env, log_fn) should live on from_model() only. Final split:

• Solver.from_model() (build-time): io_api, sanitize_*, reformulate_sos, explicit_coordinate_names, set_names, slice_size, progress, problem_fn, log_fn, env.
• Solver.solve() (run-time only): solution_fn, warmstart_fn, basis_fn.

Cost: Solver.from_model() grows; need to decide on-by-default vs. opt-in per behavior (sanitization probably on; reformulate_sos probably opt-in to preserve current semantics). Benefit: Solver path becomes first-class. Canonical low-level pattern collapses to:

solver = Solver.from_name("gurobi", model, io_api="direct")
result = solver.solve()
model.assign_result(result)