# Complementarity semantics — the source contract **Status:** normative for the modeling layer, every complementarity lowering, and every model-rebuilding pass. **Issue:** [#1147](https://github.com/jkitchin/discopt/issues/1147) (slice 1 of the MCP/MPEC RFC, [#1123](https://github.com/jkitchin/discopt/issues/1123)), the companion of [disjunction semantics](disjunction_semantics.md) (#1124). This document fixes what a complementarity condition *is* at the source level, so a later consumer — a source residual, a KKT diagnostic, a native MCP backend — has a contract to read rather than a lowering to reverse-engineer. ## 1. Why the source relation has to survive `Model.complementarity(x, y)` lowers the condition on the spot: it becomes a selector binary, big-M rows and generated auxiliaries. Before #1147 the *relation* was dropped by every pass that rebuilds a `Model` — GDP lowering, integer-product expansion, factorable reformulation, binary-multilinear linearization — so after lowering the only surviving trace of a declared pair was its name, baked into identifiers like `_gdp_aux_disj_pair0_0_0`. That is not a cosmetic loss. The next slice ([#1148](https://github.com/jkitchin/discopt/issues/1148)) must report a **source** complementarity residual, `r = max_i min(f_i(x), g_i(x))` over the *declared* operands, separately from the lowered NLP's own residual. With the source operands gone, such a probe would have measured the relaxed row instead, printed a small number, and been believed — the instrument-that-measures-nothing failure the measurement discipline in `CLAUDE.md` §6–§7 exists to prevent. So: **the relation is the durable object; the rows are a lowering of it.** ## 2. The two forms ### NCP pair — `0 <= f ⊥ g >= 0` Both operands are declared nonnegative and at most one is nonzero (`f·g = 0`). This is what `Model.complementarity` and the `.nl` type-5 importer produce, and the only form any lowering in `discopt.mpec` encodes today (GDP disjunction, SOS1, Scholtes regularization). ### Box MCP — `F(z) ⊥ z in [l, u]` A *residual* paired with a *bounded* variable: | state of `z` | asserted | |---|---| | `z = l` | `F(z) >= 0` | | `l < z < u` | `F(z) = 0` | | `z = u` | `F(z) <= 0` | `l = 0, u = +inf` reduces exactly to the NCP pair, and `Model.mcp` records that case *as* an NCP pair and lowers it, so the two forms cannot diverge on the case they share. (Bounds at discopt's ±9.999e19 unbounded sentinel are normalized to infinities first; comparing `ub == inf` literally would misclassify `m.continuous("a", lb=0)` as a general box.) The general box form is **represented and not lowered** in this slice. It is a three-branch condition, so handing it to a two-branch lowering would encode a different feasible set and certify it. Both entry points refuse loudly instead: `mpec.reformulate_{gdp,sos1,scholtes}` raise `NotImplementedError`, and the **solver boundary** refuses any model carrying a relation no lowering emitted rows for — a declared condition solved as if it were absent is a false certificate, not a missing feature. That boundary is `solve_model`, not `Model.solve`. `Model.solve` checks too (so an alternate backend is refused before it dispatches), but the check must sit where every solve passes through: the differentiable-solve paths and the primal heuristics call `solve_model` directly, and with the guard only on `Model.solve` they certified `optimal` on a point the declared relation forbids. Use `mpec.require_all_relations_lowered(model, context=...)` from any new entry point. ## 3. What the relation carries `mpec.Complementarity` is compared and hashed **by identity** (`eq=False`): a relation is a node in the model's IR, not a value, so it stays a usable dict/set key as passes carry it from one `Model` to the next. (Value equality was never available: `Expression.__eq__` builds a `Constraint`.) | field | meaning | |---|---| | `f`, `g` | the **source** operands, never rewritten by a lowering or a pass | | `role` | `ComplementarityRole`: `NCP_PAIR`, `BOX_MCP`, `FROM_KKT`, `FROM_DISJUNCT` | | `f_bounds`, `g_bounds` | the bounds the *relation* declares — its semantics, not the operand's interval enclosure, which presolve moves | | `scale` | characteristic residual magnitude, for a meaningful tolerance when the operands carry unrelated physical magnitudes; `effective_scale` derives one from the declared bounds when it is not given | | `parent` | identity of the construct that generated the relation (the lower-level row for a KKT pair, the disjunction for a disjunct pair) | | `f_shape`, `g_shape`, `index`, `source` | shape/index information for vectorized relations: `elements(model)` yields one scalar relation per index, each with its multi-index and a back-pointer to the declared relation (a scalar relation yields `[self]`, so the identity key holds there too) | | (no lowering field) | which models carry the generated rows, and by which method, is **model** state — `Model._lowered_complementarities`, read via `pair.is_lowered_into(model)` / `pair.lowering_in(model)` | `role` is *provenance*, never a lowering switch. Every lowering and every bound rule branches on `Complementarity.is_symmetric_nonnegative` — the **declared bounds** — and never on `role`. This is not a stylistic preference: `Complementarity` is public and constructible directly, so a relation with `g_bounds=(-1, 1)` and the default `role=NCP_PAIR` reaches the lowerings. A `role`-gated refusal let exactly that through and emitted `g >= 0` against a declared `lb = -1`, certifying a point the relation forbids. A label can be wrong; the declared bounds *are* the semantics. `FROM_DISJUNCT` is reserved vocabulary with no in-tree producer yet — exactly as `DisjunctionSemantics.OR` was defined by #1124 before any lowering emitted it. The lowering marks — **and the method** — live on the model (`Model._lowered_complementarities`, an identity map from relation to method), not on the relation. The mark was weak references on the relation first, which made a model carrying one unpicklable and made `copy.deepcopy` silently drop it, so the clone refused to solve. The *method* then stayed behind as a plain field, which was the same mistake one level down: lowering a relation by GDP into `m1` and by SOS1 into `m2` left the field reading `"sos1"`, so `m1`'s provenance no longer described its own disjunctive rows. The rule this settles: **a relation is shared; what a given model did with it is that model's fact.** Ask the model that carries the rows you are reasoning about, never the relation. ## 4. Provenance is keyed to objects, never to indices or names Presolve/FBBT eliminate and renumber variables, so an index-keyed map is stale exactly when it becomes useful; and a name-keyed one matches a generated `_gdp_aux_disj_pair0_0_0` that is not the source operand at all. Therefore: * `resolve_source_variables(model, pair, context=...)` resolves by **object identity** and raises `ComplementarityProvenanceError` — naming the relation and the caller — when the model does not hold one of the source variables. * `flat_source_indices(model, pair)` derives backend-facing flat columns **on demand at the solver boundary**; no index is ever persisted on the relation. It is element-accurate for a statically indexed operand (`x[2] ⊥ y` returns the two columns it reads, not all of `x`), and widens to the whole variable for an index it cannot resolve statically — conservative, never narrower than the truth, so a consumer can rely on the set covering everything the relation reads. * Those columns are prefix sums over **the target model's own `_variables`**, not `Variable._index`. `_index` is the variable's position in the model it was *declared* on, and a relation's operands are shared objects: a model holding a two-element `x` and a scalar `y` as `[y, x]` has them at `[1, 2]` and `[0]`, while `_index` reports `[0, 1]` and `[1]`. `Model._flat_var_offset` is right for the rebuilding passes (they preserve declaration order) and wrong as a general provenance accessor. A provenance query must also never *write* shared state — no renumbering, no renaming — since other live models hold the same objects. * An expression node the walker cannot descend into raises rather than contributing nothing, so resolution can never report success without having looked at the operand. ## 5. What every rebuilding pass owes the relation A pass that constructs a fresh `Model` **must** call `mpec.carry_complementarities(src, dst, pass_name=...)`. It forwards the same relation objects (so identity-keyed provenance holds across any number of passes), checks each against `dst` first, and propagates the lowering mark with `src`'s own method — `dst` was rebuilt from `src`'s rows, so it already carries the generated ones and must not be re-lowered. The propagation test is `pair.is_lowered_into(src)` and nothing weaker. Keying it on "the relation has been lowered *somewhere*" let an unlowered `src` sharing a relation with another model hand `dst` a lowered mark, so `dst` claimed rows neither model had and walked straight past the solve guard of §2. Passes with a defensive `except Exception: return model` handler re-raise `ComplementarityProvenanceError` ahead of it. A broken provenance chain is not an out-of-scope model; degrading it to a silent abstain is the failure mode this contract exists to remove. One guard was *removed* under this contract rather than kept: `binary_multilinear_reform` used to abstain whenever `_complementarities` was non-empty — a guard keyed on state every earlier pass emptied, so its premise evaporated at a pass boundary. Forwarding the relation set replaces it. ## 6. Naming A relation's name is the base name of every row generated from it, so an unnamed relation must not reuse one. `_ensure_relation_name` assigns `compl{k}` unique against the model's existing relation and constraint names, at the single point every lowering funnels through (`_scalarize_pairs`) — not at a call site. The fallback previously counted *within the list handed to the lowering*, and `Model.complementarity` hands over exactly one relation, so every unnamed declaration on a model became `compl0` and emitted `compl0_f_nonneg` twice. ## 7. Source residuals — what a number means **Status:** normative for every consumer that reports a complementarity residual. **Issue:** [#1148](https://github.com/jkitchin/discopt/issues/1148) (slice 2), implemented in `discopt.mpec_report`. ### The two residuals are different numbers POUNCE's internal barrier/KKT complementarity is a property of the **generated NLP**. The MPCC complementarity `0 <= f ⊥ g >= 0` is a property of the **user's model**. They are not the same measurement and must never be reported as if they were. Measured on `max x + y s.t. x == y, 0 <= x ⊥ y >= 0`, after a Scholtes homotopy stopped at `t = 1e-8`: | quantity | value | |---|---| | max violation of the generated rows (`x >= 0`, `y >= 0`, `x·y <= t`) | `1.0e-8` | | source complementarity `min(x, y)` on the declared operands | `1.4e-4` | Four orders of magnitude, at the same point, because the *regularized* problem admits `x = y ≈ sqrt(t)` and every row it holds is satisfied there. A result that reports only the first number says the point is complementary. It is not: the true MPCC forces `x = y = 0`, and the returned objective (`-2.8e-4`) is *below* the true optimum (`0`). `SourceResidualReport` carries both, side by side. ### What is reported, and its definition Every residual is a `Residual` carrying `value`, the `definition` string that produced it, the `scale` it was measured against, and the algorithmic `floor` below which it cannot be driven. A number without its formula cannot be compared across solvers, which is why the definition is a field and not a docstring. | residual | definition | |---|---| | `source_complementarity` | `max_i min(f_i(x), g_i(x))` for a nonnegative pair; the normal map `\|z - mid(l, u, z - F(x))\|` for a box MCP. `product` and Fischer–Burmeister are selectable. | | `source_operand_bounds` | violation of the bounds the **relation** declares on each operand (nonnegativity for an NCP pair, the box for an MCP) | | `source_primal_feasibility` | max violation of the **source** constraint rows and the **original** declared variable bounds, on the source expression tree | | `selector_integrality` | `max_i y_i(1 - y_i)` over binary selectors, **after** checking `0 <= y_i <= 1` | | `lowered_row_residual` | max violation of the rows the lowering generated — the comparison quantity, never the source claim | Which definition is used is chosen from the relation's **declared bounds**, never from its `role` — the same rule §3 states for the lowerings, and for the same reason. ### Scale, and the residual the relaxation admits A uniform tolerance is meaningless when the two operands carry unrelated physical magnitudes (a multiplier in 1e-6 against a flow in 1e3), so a residual is reported both raw and divided by the relation's `effective_scale` (§3). Separately, a Scholtes regularization at `t` replaces the exact condition with `f, g >= 0` and `f·g <= t`, and that relaxed set **admits** points whose residual is nonzero. How large depends on which residual definition is being read: | definition | largest residual the row `f·g <= t` admits | | --- | --- | | `min`, `natural_map` | `sqrt(t)`, attained at `f = g = sqrt(t)` | | `product` | `t` — the row bounds this quantity directly | | `fischer_burmeister` | `(2 - sqrt(2))·sqrt(t)`, also attained at `f = g` | `Residual.admitted_scale` carries that number and `admitted_scale_definition` records the formula it came from, so the two can never be read against each other's units. `Residual.within_admitted_scale` compares the residual against that upper bound and nothing else. **This comparison alone does not establish feasibility of the relaxed rows.** Only one direction is implied: for a nonnegative pair, a `min` residual *above* `sqrt(t)` forces `f·g > t`, so being outside the scale does mean the relaxed row is violated. The converse fails — at `t = 1e-8`, `(f, g) = (1e-4, 1)` has `min(f, g) = 1e-4 = sqrt(t)` so the flag is true, while `f·g = 1e-4` exceeds `t` by four orders of magnitude. Feasibility of the generated rows is a separate number: `lowered_row_residual`. **This is an upper bound on the admitted residual, not a lower limit on attainable accuracy.** `(f, g) = (0, 1)` satisfies `f·g <= t` for every positive `t` with a residual of exactly `0`, and solvers routinely land far below the scale — so a residual much smaller than `sqrt(t)` is ordinary, not surprising. What the number buys is the other direction: a `min` residual of `1e-4` from a `t = 1e-8` continuation is no larger than the worst case the relaxation admits, which is a different statement from "the subsolver failed to converge" — and, per the previous paragraph, also a different statement from "the relaxed rows hold at this point". Whether the subsolver could have done better is an empirical question about the subsolver; the continuation trace's per-stage statuses are where that evidence lives, not here. POUNCE Gate 0 ([jkitchin/pounce#794](https://github.com/jkitchin/pounce/issues/794)) documents the `sqrt(t)` scale in exactly this sense. A report that omitted it would imply exact orthogonality. ### This field set is the shared benchmark schema [jkitchin/pounce#780](https://github.com/jkitchin/pounce/issues/780) lists "define the common benchmark metadata/result schema" as an immediate action. That schema is this one: `SourceResidualReport.as_dict()` is its serialization, and it is defined here rather than as a separate artifact so the two repositories do not end up with two vocabularies for the same numbers. Every entry carries `value`, `definition`, `scale`, `scaled_value`, `admitted_scale`, `admitted_scale_definition` and `within_admitted_scale`; the report adds `n_scalar_relations` (the executed-measurement count), the per-relation breakdown with `source_name`/`index` attribution, the continuation trace, and `stationarity` (always `null` — see §8). ### The instrument checks itself Operands are evaluated with the interval evaluator over a **degenerate** box (`[x, x]` per variable, keyed by object identity). A correct evaluation stays degenerate; the evaluator widens to `[-inf, +inf]` for an atom it cannot walk, so a widened result proves the graph was *not* evaluated and `evaluate_at_point` raises rather than reporting a midpoint. Likewise `relation_residuals` refuses when it is handed relations and measures none, and every constraint-list entry kind (`Constraint`, `_DisjunctiveConstraint`, `_SOSConstraint`) has an explicit residual or raises — a row silently skipped would make the maximum read as complete. This is CLAUDE.md §6 applied to the instrument itself. ## 8. Local results versus certificates — the contract **Status:** normative for every solver path that can return a local result. **Issue:** [#1148](https://github.com/jkitchin/discopt/issues/1148) §B/§C. Vocabulary in `discopt.status`; enforcement in `SolveResult.__post_init__`. ### A distinct terminal status, not a flag `discopt.status` defines three statuses that make **no global claim**: * `local_optimal` — a local stationary point of the problem actually solved; * `local_limit` — a *point*, from a local search that stopped on an iteration or step-size limit without establishing stationarity. Strictly weaker than `local_optimal`, and separate from it for the same reason `local_infeasible` is separate from `infeasible`: the weaker claim needs its own name, or it gets made under the stronger one. A Scholtes homotopy whose final accepted stage stalled reports the iterate here — it is still a usable warm start and still carries its residuals — rather than calling it stationary; * `local_infeasible` — a local solver found no point. This is **not** an infeasibility proof. A stalled MPEC continuation lands here and never on `infeasible`, which is a certificate in the other direction. The continuation trace keeps the same distinction, so the result and the trace cannot disagree: `ContinuationStage.certified` says whether the subsolver converged at that stage (as opposed to merely handing back the point it stopped at), `ContinuationTrace.reported_point_certified` says whether the stage behind the *reported* point did, and `converged` reads that rather than `any_stage_accepted`. It is a status and not a `certified=False` flag beside `"optimal"` because every consumer that pattern-matches on the *status* would inherit the bug. Measured on the benchmark harness: `is_solved` is `status == OPTIMAL`, and `proved_optimal_count` counts `gap is None and is_solved` as a proved optimum — a local stationary point reported as `"optimal"` is scored by the release gate as a **certificate**, which is CLAUDE.md §1's hard gate with no slack. An unmapped status maps to `UNKNOWN` and both counters skip the row, so a consumer that has never heard of the local mode fails closed. `DISCOPT_STATUS_MAP` names `local_optimal`/`local_limit`/`local_infeasible` explicitly so that skip is intentional and tested rather than accidental. The states the vocabulary keeps apart: certified optimal; certified bound + incumbent with a gap; feasible with no bound; local stationary with no global claim; a local point with no stationarity claim at all; and certified-infeasible versus local-search-failed. ### The asymmetry that makes a local mode safe > A local result **may become an incumbent only after independent feasibility > verification**, and **may never become a bound.** A feasible point is a valid upper bound on a minimum; a dual bound is the certificate and never comes from a local solve. That asymmetry is what lets a local mode feed the global solver without contaminating it. It is enforced structurally, not by convention: `SolveResult.__post_init__` **raises** if a local status carries a `bound` or a `root_bound`, and forces `gap_certified` to `False`. Refusing rather than silently clearing the field is deliberate — a caller that set one has a bug in what it believes it proved. `mpec_report.accept_local_incumbent` is the verification side: it runs the shared `validation.feasibility.verify_point` and returns `None` rather than vouching for a point it cannot verify. Two further rules: * **Validation is performed on the source model**, not only on the generated NLP. A smoothing tolerance can be small while a truth value near a strict predicate boundary has already flipped — §7's table is that case. * **Absent certification is interpreted as not certified.** `status.is_certified_status` is a membership test against a closed set, never `not is_local_status(...)`. ### Stationarity is not claimed C-/M-/S-stationarity classifications are reported **only** if discopt has actually checked the required conditions. It has not, so `SourceResidualReport.stationarity` is `None` and no code path sets it. Reporting a classification the solver did not verify would be exactly the §1 error this document exists to prevent. ### One return type `mpec.solve_mpec` returns a `modeling.core.SolveResult` for `scholtes`, `sos1` **and** `gdp`. Before #1148 switching one keyword changed the returned type, the type of `status` (enum versus str) and whether a certification field existed at all, so no caller could write one branch that read a result. The Scholtes arm carries its continuation trace — per-stage `t`, subsolver status, accept/reject reason and the achieved source residual — rather than discarding it, so "converged at 1e-8" and "stalled at 1e-2" are distinguishable from the result alone.