diff --git a/CHANGELOG.md b/CHANGELOG.md index 0a2bc447..fea79a0c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,8 +14,57 @@ The release procedure that produces these entries is documented in ### Changed +- **`gap_certified` now requires BOTH ends of a gap** (`fix`, #875). A limit exit + with a valid dual bound but *no incumbent* used to report + `status="time_limit", objective=None, gap=None, gap_certified=True` — a certified + gap where no gap was ever formed. `SolveResult.__post_init__` already required a + finite bound; it now also requires an objective, for every status except + `infeasible` (whose certificate is not a gap). The downgrade is True → False only, + so it cannot manufacture a certificate, and the dual `bound` is left in place + because a bound with no incumbent is still a valid bound. + ### Fixed +- **Root setup no longer scales with `n_vars × n_constraints`, and no longer runs + unbudgeted** (`fix(performance)`, #875, successor to #863/#868). After #868, + `watercontamination0202` (106,711 vars / 107,209 rows) returned instead of hanging + but still took 579 s against a 30 s `time_limit` with `nodes=0`. Two distinct + defects, both entirely before the first branch-and-bound node: + - *Cost.* `_fix_single_var_equalities` was `O(n_constraints × n_vars)` — the affine + linearizer it calls allocated and zeroed a dense `n_vars` array per constraint, + and the caller then walked all `n_vars` entries in Python to find the single + nonzero — for bodies with **one leaf each**. That was ~460 s of the 579 s (23 of + 29 stack samples). Adds `_linearize_affine_expr_sparse` as the core, with the + dense `_linearize_affine_expr` as a view of it, and routes every per-constraint + scan through the sparse form (`_fix_single_var_equalities`, the one-hot group + scan, the monomial-factor fallback, the affine bound helpers). A new + `_any_linear_constraint_form` short-circuits the RLT-applicability probe, which + previously built and kept one dense row vector per constraint — 91 GB on that + instance — to evaluate a boolean. Measured on a synthetic probe at a fixed + constraint count: 0.85 s → 0.001 s at `n_vars=32,000`, and flat in `n_vars` + (18.8× → 1.0× across a 16× variable increase) instead of exactly linear. + - *Budget.* `tighten_nonlinear_bounds` (80.9 s over three root calls) had no + `deadline` parameter at all; it now takes one and polls at three granularities — + between rounds, between rules, and inside a rule's constraint scan — because a + single rule's sweep over 107k rows already exceeds a tight budget (#868's + `probing` lesson: a poll at the wrong granularity is not a poll). Rules declare + `row_scan_is_anytime`; the default is `False`, so a rule whose conclusion rests + on a row being *absent* (`PeriodicVariableBoundRule`) is skipped whole rather + than truncated. The convexity classifier's budget is a fraction of `time_limit` + recomputed per *model object*, so a reformulation restarted it and the fractions + added up (14.7 s over two runs); it is now additionally clamped to the absolute + `model._solve_deadline`. Every early exit does strictly less work — a looser box, + fewer rules, `convexity-unknown` routing to sound spatial branch-and-bound — + never a wrong one. +- **The QCP probe extractor's Hessian may be sparse** (`fix`, #875). + `_extract_quadratic_coefficients_from_values` still swept all `O(n²)` variable + pairs into a dense `(n, n)` array, per constraint — the one extractor #863/#868 + left dense. It now restricts off-diagonal probing to the evaluator's support (free: + the `O(n)` diagonal probes already identify it) and materialises through + `_materialise_Q`, matching the QP path. `_quadratic_row_has_terms` and the Gurobi + QCP entry point are made sparse-aware in step, so a sparse `Q` cannot reach a + consumer as the 0-d object array `np.asarray` silently produces. + ## [0.7.0] - 2026-07-24 ### Added diff --git a/python/discopt/_jax/mccormick_lp.py b/python/discopt/_jax/mccormick_lp.py index be815465..63d11ff8 100644 --- a/python/discopt/_jax/mccormick_lp.py +++ b/python/discopt/_jax/mccormick_lp.py @@ -29,8 +29,8 @@ from discopt._jax.milp_relaxation import ( _LIFT_MAX_CROSS_TERM_ARG_MAGNITUDE, _RELAX_NUMERIC_CAP, + _any_linear_constraint_form, _collect_affine_powers, - _linear_constraint_forms, _quadratic_constraint_forms, build_milp_relaxation, ) @@ -468,8 +468,12 @@ def __init__( # models (e.g. the indefinite integer QPs nvs17/nvs24, whose constraints # are all quadratic) from RLT entirely — exactly the dense-QP instances # whose McCormick bound is hopelessly loose without it. + # ``_any_linear_constraint_form`` short-circuits and stays sparse: the old + # ``bool(_linear_constraint_forms(...))`` linearized every row into a dense + # ``n_orig`` vector and kept them all just to test emptiness — O(rows * vars) + # in the relaxer constructor, on the pre-B&B root path (issue #875). self._rlt_applicable = (rlt_level1 or _tuning().rlt) and ( - bool(_linear_constraint_forms(model, self._n_orig)) + _any_linear_constraint_form(model, self._n_orig) or bool(_quadratic_constraint_forms(model, self._n_orig)) ) # Lever 1 (issue #194): solve each spatial-B&B node's relaxation as a pure diff --git a/python/discopt/_jax/milp_relaxation.py b/python/discopt/_jax/milp_relaxation.py index d9c86cc1..714daa3c 100644 --- a/python/discopt/_jax/milp_relaxation.py +++ b/python/discopt/_jax/milp_relaxation.py @@ -1237,24 +1237,75 @@ def _linear_constraint_forms(model: Model, n_vars: int) -> list[tuple[np.ndarray ``-g <= 0``. Nonlinear constraint bodies are skipped (they have no affine form). These are the factors level-1 RLT multiplies by variable bound factors to generate valid product cuts. + + Callers that only need to know *whether* any such factor exists must use + :func:`_any_linear_constraint_form` instead — materializing the dense arrays to + take ``bool(...)`` of the list is ``O(n_constraints * n_vars)`` in both time and + memory (issue #875). """ forms: list[tuple[np.ndarray, float]] = [] for constraint in model._constraints: - try: - coeff, const = _linearize_affine_expr(constraint.body, model, n_vars) - except ValueError: - continue # nonlinear body — not an affine factor - if not (np.all(np.isfinite(coeff)) and np.isfinite(const)): - continue sense = constraint.sense + if sense not in ("<=", "=="): + continue # no form is emitted below; skip the linearization entirely + terms = _linear_form_terms(constraint.body, model, n_vars) + if terms is None: + continue + coeff = np.zeros(n_vars, dtype=np.float64) + for j, c in terms[0].items(): + coeff[j] = c + const = terms[1] if sense == "<=": - forms.append((coeff, float(const))) - elif sense == "==": - forms.append((coeff, float(const))) - forms.append((-coeff, -float(const))) + forms.append((coeff, const)) + else: + forms.append((coeff, const)) + forms.append((-coeff, -const)) return forms +def _linear_form_terms( + body: Expression, model: Model, n_vars: int +) -> Optional[tuple[dict[int, float], float]]: + """Sparse ``(terms, const)`` for a linear constraint body, or ``None``. + + ``None`` means "not an affine factor": either the body is nonlinear (the + linearizer refuses) or a coefficient/constant is non-finite. Shared by + :func:`_linear_constraint_forms` and :func:`_any_linear_constraint_form` so both + apply the same acceptance test. + """ + try: + terms, const = _linearize_affine_expr_sparse(body, model, n_vars) + except ValueError: + return None # nonlinear body — not an affine factor + if not np.isfinite(const): + return None + for j, c in terms.items(): + # The dense form raised IndexError (uncaught) on an out-of-range index and + # dropped a non-finite coefficient; keep the drop, and treat an unindexable + # coefficient the same way rather than writing outside the array. + if not (0 <= j < n_vars) or not np.isfinite(c): + return None + return terms, float(const) + + +def _any_linear_constraint_form(model: Model, n_vars: int) -> bool: + """True when the model has at least one linear constraint factor for RLT. + + Short-circuits on the first hit, and never allocates an ``n_vars`` array. + ``bool(_linear_constraint_forms(model, n_vars))`` answers the same question but + linearizes every row into a dense vector first and *keeps them all* — on + ``watercontamination0202`` (106,711 vars / 107,209 rows) that is ~91 GB of + coefficient arrays to compute one boolean, and it sits in the + ``MccormickLPRelaxer`` constructor on the pre-B&B root path (issue #875). + """ + for constraint in model._constraints: + if constraint.sense not in ("<=", "=="): + continue + if _linear_form_terms(constraint.body, model, n_vars) is not None: + return True + return False + + # A quadratic constraint factor for level-1 RLT (issue #15, Phase 2): the body # ``g(x) = const + sum_i lin_i x_i + sum_{(k,l)} quad_{kl} x_k x_l`` with the sense # of the parent constraint. ``quad`` keys are sorted index pairs ``(k, l)``, @@ -1317,14 +1368,33 @@ def _quadratic_constraint_forms(model: Model, n_vars: int) -> list[_QuadForm]: return forms -def _linearize_affine_expr(expr: Expression, model: Model, n_vars: int) -> tuple[np.ndarray, float]: - """Linearize an affine expression over original variables. - - Raises ValueError when the expression contains nonlinear structure: only - affine arguments are soundly supported here, since any nonlinear structure - is relaxed by the uniform engine's atom envelopes rather than linearized. +def _linearize_affine_expr_sparse( + expr: Expression, model: Model, n_vars: int +) -> tuple[dict[int, float], float]: + """Sparse core of :func:`_linearize_affine_expr`: ``({flat_index: coeff}, const)``. + + Same recognition rules and same ``ValueError`` refusals as the dense wrapper — + only the accumulator differs. Callers that just need the *support* of an affine + body (how many variables it touches, and which) should use this: the dense + wrapper costs ``O(n_vars)`` per call to allocate and zero its array, so scanning + every constraint through it is ``O(n_constraints * n_vars)`` no matter how sparse + the bodies are. + + Measured (issue #875, ``watercontamination0202``: 106,711 vars / 107,209 rows): + ``_fix_single_var_equalities`` — a scan whose bodies have a SINGLE leaf each — + spent ~460 s of a 30 s solve budget entirely in that dense allocation plus the + caller's Python walk over all 106,711 entries. The synthetic scaling probe is + exactly linear in ``n_vars`` at a fixed constraint count (0.068 s at n=2,000 → + 3.650 s at n=128,000 for 400 rows), which is the signature of the dense array + rather than of the bodies. + + Keys are the flat original-variable indices; a key is emitted only for a + variable the body actually references, so the dict size is ``O(body leaves)``. + Indices are NOT range-checked here (the dense wrapper's array indexing is what + enforces ``n_vars``); a sparse caller that indexes arrays with them must bound + them itself. """ - coeff = np.zeros(n_vars, dtype=np.float64) + coeff: dict[int, float] = {} const_acc: list[float] = [0.0] def visit(e: Expression, scale: float) -> None: @@ -1335,7 +1405,7 @@ def visit(e: Expression, scale: float) -> None: if isinstance(e, Variable): offset = _compute_var_offset(e, model) if e.size == 1: - coeff[offset] += scale + coeff[offset] = coeff.get(offset, 0.0) + scale return raise ValueError(f"Cannot use array variable as scalar affine argument: {e}") @@ -1343,7 +1413,7 @@ def visit(e: Expression, scale: float) -> None: flat = _get_flat_index(e, model) if flat is None: raise ValueError(f"Cannot linearize IndexExpression: {e}") - coeff[flat] += scale + coeff[flat] = coeff.get(flat, 0.0) + scale return if isinstance(e, UnaryOp) and e.op == "neg": @@ -1388,7 +1458,7 @@ def visit(e: Expression, scale: float) -> None: if isinstance(op, Variable): offset = _compute_var_offset(op, model) for k in range(op.size): - coeff[offset + k] += scale + coeff[offset + k] = coeff.get(offset + k, 0.0) + scale return visit(op, scale) return @@ -1404,6 +1474,26 @@ def visit(e: Expression, scale: float) -> None: return coeff, const_acc[0] +def _linearize_affine_expr(expr: Expression, model: Model, n_vars: int) -> tuple[np.ndarray, float]: + """Linearize an affine expression over original variables. + + Raises ValueError when the expression contains nonlinear structure: only + affine arguments are soundly supported here, since any nonlinear structure + is relaxed by the uniform engine's atom envelopes rather than linearized. + + Dense view of :func:`_linearize_affine_expr_sparse` — identical coefficients, + identical constant, identical refusals, and (as before) an ``IndexError`` from + the array store when a body references a flat index at or beyond ``n_vars``. + Prefer the sparse core in any per-constraint scan: this allocates and zeroes an + ``n_vars`` array per call (issue #875). + """ + terms, const = _linearize_affine_expr_sparse(expr, model, n_vars) + coeff = np.zeros(n_vars, dtype=np.float64) + for j, c in terms.items(): + coeff[j] = c + return coeff, const + + def _match_scaled_constant_division( expr: Expression, scale: float, @@ -1545,16 +1635,19 @@ def visit(e: Expression, power: float) -> bool: # positive monomial factor. A nonzero constant or >1 variable breaks the # monomial form (reject); a non-positive coefficient breaks the log lift. try: - aff_coeffs, aff_const = _linearize_affine_expr(e, model, n_orig) + # Sparse: this runs per factor of every monomial body, and the dense + # form allocated an ``n_orig`` array plus a Python walk over it just to + # count the nonzeros (#875). + aff_terms, aff_const = _linearize_affine_expr_sparse(e, model, n_orig) except ValueError: return False if abs(float(aff_const)) > 1e-12: return False - nz = [(j, float(c)) for j, c in enumerate(aff_coeffs) if abs(float(c)) > 0.0] + nz = [(j, c) for j, c in aff_terms.items() if abs(c) > 0.0] if len(nz) != 1: return False j, c = nz[0] - if j >= n_orig or c <= 0.0: + if j < 0 or j >= n_orig or c <= 0.0: return False lo = float(flat_lb[j]) if not (lo > 1e-9) or not np.isfinite(lo): @@ -2285,18 +2378,25 @@ def _integer_affine_cos_lower_bound( """Return exact lower bound for scale*cos(integer-affine expr) on a small box.""" if not isinstance(expr, FunctionCall) or expr.func_name != "cos" or len(expr.args) != 1: return None + n_vars = len(flat_lb) try: - coeff, const = _linearize_affine_expr(expr.args[0], model, len(flat_lb)) + terms, const = _linearize_affine_expr_sparse(expr.args[0], model, n_vars) except ValueError: return None flat_types = _flat_variable_types(model) entries: list[tuple[float, range]] = [] n_values = 1 - for var_idx, c_i in enumerate(coeff): - c = float(c_i) + # Ascending index order, matching the dense ``enumerate(coeff)`` walk this + # replaces — the enumeration below is over a cartesian product, so the entry + # order is not observable in the result, but keeping it makes the two forms + # trivially comparable. Sparse: the dense walk was O(n_vars) per call (#875). + for var_idx in sorted(terms): + c = terms[var_idx] if abs(c) <= 1e-12: continue + if not (0 <= var_idx < n_vars): + return None # unindexable support: the dense form raised here values = _integer_domain_values(var_idx, flat_types, flat_lb, flat_ub) if values is None: return None @@ -2325,16 +2425,24 @@ def _scaled_affine_lower_bound( flat_lb: np.ndarray, flat_ub: np.ndarray, ) -> Optional[float]: + n_vars = len(flat_lb) try: - coeff, const = _linearize_affine_expr(expr, model, len(flat_lb)) + terms, const = _linearize_affine_expr_sparse(expr, model, n_vars) except ValueError: return None want_lower = scale >= 0.0 bound = float(const) - for c_i, lb_i, ub_i in zip(coeff, flat_lb, flat_ub): - c = float(c_i) + # Sparse walk over the body's support; the dense form zipped all ``n_vars`` + # coefficients against the box on every call (#875). Ascending index order keeps + # the floating-point accumulation identical to the dense walk's. + for var_idx in sorted(terms): + c = terms[var_idx] if abs(c) <= 1e-12: continue + if not (0 <= var_idx < n_vars): + return None # unindexable support: the dense form raised here + lb_i = flat_lb[var_idx] + ub_i = flat_ub[var_idx] chosen = float(lb_i) if (c >= 0.0) == want_lower else float(ub_i) if not _is_effectively_finite(chosen): return None diff --git a/python/discopt/_jax/nonlinear_bound_tightening.py b/python/discopt/_jax/nonlinear_bound_tightening.py index 75368fbc..f32bfe29 100644 --- a/python/discopt/_jax/nonlinear_bound_tightening.py +++ b/python/discopt/_jax/nonlinear_bound_tightening.py @@ -8,7 +8,9 @@ from __future__ import annotations +import inspect import logging +import time from dataclasses import dataclass from typing import Callable, NoReturn, Optional, Sequence, TypeVar @@ -174,6 +176,45 @@ class NonlinearBoundTighteningStats: applied_rules: tuple[str, ...] infeasible: bool = False infeasibility_reason: Optional[str] = None + #: True when the pass stopped early because its ``deadline`` passed, so the + #: returned box is a valid but *incomplete* tightening (fewer rules ran, or a + #: rule saw only a prefix of the constraints). Never set when no deadline was + #: given. Purely informational — the box itself is sound either way. + deadline_reached: bool = False + + +# Constraints between two ``time.perf_counter()`` reads inside a rule's row scan. +# The read is ~50 ns and a row's match attempt is orders of magnitude more, so the +# stride is free; it exists so the poll cannot be starved by a cheap-row model. +# Granularity matters more than the poll's existence: #868 found ``probing`` polling +# its deadline once per binary on an instance with 7 binaries and 107k constraints, +# which is a poll that never fires in time. +_DEADLINE_POLL_STRIDE = 64 + + +def _rows_until(model: Model, deadline: Optional[float]): + """Iterate ``model._constraints``, stopping once ``deadline`` passes. + + ``deadline`` is an absolute ``time.perf_counter()`` value, or ``None`` for the + unbounded legacy behaviour (the constraint list is then returned as-is, so a + budget-free call has *zero* added cost and is byte-identical to before). + + Every rule in this module only ever *intersects* the incoming box, and each + constraint's inference is independently valid, so consuming a prefix of the rows + yields a valid — merely looser — box. Stopping early is a "do less" bail, never + a "do wrong" one (issue #875; same contract as the Rust presolve passes' `_until` + variants from #868). + """ + if deadline is None: + return model._constraints + return _iter_rows_until(model._constraints, deadline) + + +def _iter_rows_until(rows, deadline: float): + for i, row in enumerate(rows): + if i % _DEADLINE_POLL_STRIDE == 0 and time.perf_counter() >= deadline: + return + yield row class NonlinearBoundTighteningInfeasible(ValueError): @@ -188,20 +229,66 @@ class _ReciprocalIntervalInfeasible: class NonlinearBoundTighteningRule: - """Base class for extensible nonlinear bound tightening rules.""" + """Base class for extensible nonlinear bound tightening rules. + + A rule MAY declare a trailing ``deadline: Optional[float] = None`` parameter on + :meth:`tighten`; :func:`tighten_nonlinear_bounds` detects that and passes the + absolute ``time.perf_counter()`` budget through, so the rule can stop scanning + rows (via :func:`_rows_until`) instead of running the whole constraint list. A + rule without the parameter keeps the four-argument contract exactly and is + simply called without it — external rules do not need to change. + """ name = "unnamed_rule" + #: True only when consuming a PREFIX of ``model._constraints`` yields a valid + #: (merely looser) box — i.e. every row's inference stands on its own and the + #: rule draws no conclusion from a row being *absent*. + #: + #: Default ``False``, deliberately: a rule that accumulates across rows and then + #: acts on the complement (``PeriodicVariableBoundRule`` restricts a variable to + #: one period precisely because nothing *else* in the model uses it) would cut + #: feasible points if a truncated scan missed the disqualifying row. Under a + #: budget such a rule is skipped whole rather than run on a prefix, so the safe + #: default is the one a new rule inherits without thinking about it. + row_scan_is_anytime: bool = False + def tighten( self, model: Model, flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: raise NotImplementedError +_RULE_ACCEPTS_DEADLINE: dict[type, bool] = {} + + +def _rule_accepts_deadline(rule: NonlinearBoundTighteningRule) -> bool: + """True when ``rule.tighten`` declares a ``deadline`` parameter (memoized). + + Signature inspection rather than a try/except around the call: a ``TypeError`` + raised from *inside* a rule must not be mistaken for "this rule has the old + signature" and silently retried unbudgeted. + """ + cls = type(rule) + hit = _RULE_ACCEPTS_DEADLINE.get(cls) + if hit is None: + try: + params = inspect.signature(cls.tighten).parameters + hit = "deadline" in params or any( + p.kind is inspect.Parameter.VAR_KEYWORD for p in params.values() + ) + except (TypeError, ValueError) as exc: + logger.debug("cannot inspect %s.tighten, assuming no deadline: %s", cls.__name__, exc) + hit = False + _RULE_ACCEPTS_DEADLINE[cls] = hit + return hit + + def _constant_value(expr) -> Optional[float]: if not isinstance(expr, Constant): return None @@ -623,6 +710,10 @@ class SumOfSquaresUpperBoundRule(NonlinearBoundTighteningRule): """Tighten bounds from constraints like sum(a_i * x_i^2) <= c.""" name = "sum_of_squares_upper_bound" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True def _match_scaled_square( self, @@ -638,11 +729,12 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) not in ("<=", "=="): continue @@ -717,6 +809,10 @@ class SqrtSumOfSquaresUpperBoundRule(NonlinearBoundTighteningRule): """Tighten bounds from constraints like sqrt(sum(a_i * x_i^2)) <= c.""" name = "sqrt_sum_of_squares_upper_bound" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True def _match_scaled_square( self, @@ -790,11 +886,12 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) not in ("<=", "=="): continue @@ -861,6 +958,10 @@ class SeparableQuadraticUpperBoundRule(NonlinearBoundTighteningRule): """Tighten bounds from separable convex quadratic constraints like x + y^2 <= c.""" name = "separable_quadratic_upper_bound" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True def _match_scaled_square( self, @@ -876,11 +977,12 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) not in ("<=", "=="): continue @@ -1083,6 +1185,10 @@ class MonotoneFunctionEqualityRule(NonlinearBoundTighteningRule): """Propagate equalities like y == exp(a*x + b) in both directions.""" name = "monotone_function_equality" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True def _match_scaled_function( self, @@ -1123,11 +1229,12 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) != "==": continue @@ -1266,6 +1373,10 @@ class MonotoneFunctionBoundsRule(NonlinearBoundTighteningRule): """Tighten affine arguments of monotone unary function constraints.""" name = "monotone_function_bounds" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True def _match_scaled_function( self, @@ -1306,11 +1417,12 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) not in ("<=", "=="): continue @@ -1399,6 +1511,10 @@ class QuadraticEqualityBoundsRule(NonlinearBoundTighteningRule): """Propagate equalities like x == y**2 in both directions.""" name = "quadratic_equality_bounds" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True @staticmethod def _square_interval(lb: float, ub: float) -> tuple[float, float]: @@ -1442,11 +1558,12 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) != "==": continue @@ -1566,6 +1683,10 @@ class SquareDifferenceLowerBoundRule(NonlinearBoundTighteningRule): """Tighten the leading square in equalities like ``a*x^2 = b*y^2 + c*z^2``.""" name = "square_difference_lower_bound" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True def _collect_square_sum( self, @@ -1626,11 +1747,12 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) != "==": continue @@ -1762,6 +1884,10 @@ class PositiveAffineReciprocalBoundsRule(NonlinearBoundTighteningRule): """Propagate ``c / positive_affine >= rhs`` into affine upper bounds.""" name = "positive_affine_reciprocal_bounds" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True def tighten( self, @@ -1769,12 +1895,13 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() n_vars = len(flat_lb) - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) not in ("<=", "=="): continue @@ -1831,6 +1958,10 @@ class NegativePowerBoundsRule(NonlinearBoundTighteningRule): """Infer strict positive lower bounds from ``x**p <= affine`` for ``p < 0``.""" name = "negative_power_bounds" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True def tighten( self, @@ -1838,12 +1969,13 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() n_vars = len(flat_lb) - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) not in ("<=", "=="): continue @@ -1924,6 +2056,10 @@ class ReciprocalBoundsRule(NonlinearBoundTighteningRule): """Tighten sign-stable affine denominators in simple reciprocal constraints.""" name = "reciprocal_bounds" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True def _match_scaled_reciprocal( self, @@ -1985,11 +2121,12 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) not in ("<=", "=="): continue @@ -2145,6 +2282,10 @@ class BilinearProductEqualityRule(NonlinearBoundTighteningRule): """ name = "bilinear_product_equality" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True def _match_product_term( self, term, metadata: FlatVariableMetadata @@ -2196,12 +2337,13 @@ def tighten( flat_lb: np.ndarray, flat_ub: np.ndarray, metadata: FlatVariableMetadata, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray]: tightened_lb = flat_lb.copy() tightened_ub = flat_ub.copy() eps = 1e-9 - for constraint in model._constraints: + for constraint in _rows_until(model, deadline): if getattr(constraint, "sense", None) != "==": continue @@ -2291,8 +2433,12 @@ class FunctionDomainBoundRule(NonlinearBoundTighteningRule): """ name = "function_domain_bound" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True - def tighten(self, model, flat_lb, flat_ub, metadata): + def tighten(self, model, flat_lb, flat_ub, metadata, deadline=None): lb = np.asarray(flat_lb, dtype=np.float64).copy() ub = np.asarray(flat_ub, dtype=np.float64).copy() @@ -2341,7 +2487,7 @@ def walk(expr): try: if model._objective is not None: walk(model._objective.expression) - for c in model._constraints: + for c in _rows_until(model, deadline): walk(c.body) except NonlinearBoundTighteningInfeasible: # An empty argument domain is a genuine infeasibility proof; let it @@ -2376,7 +2522,17 @@ class PeriodicVariableBoundRule(NonlinearBoundTighteningRule): name = "periodic_variable_bound" - def tighten(self, model, flat_lb, flat_ub, metadata): + # NOT anytime (the base-class default, restated because this rule is the reason + # the default exists). The reduction is justified by what the model does *not* + # contain: a variable stays in ``periodic - disqualified`` only while no scanned + # row uses it outside ``sin``/``cos``. A scan cut short at a deadline could miss + # exactly the disqualifying row and shrink a variable that is genuinely used + # elsewhere to one period — cutting feasible points. Under a budget this rule is + # therefore skipped whole; ``deadline`` is accepted and deliberately unused. + row_scan_is_anytime = False + + def tighten(self, model, flat_lb, flat_ub, metadata, deadline=None): + del deadline # see row_scan_is_anytime above: this scan must be complete lb = np.asarray(flat_lb, dtype=np.float64).copy() ub = np.asarray(flat_ub, dtype=np.float64).copy() try: @@ -2439,7 +2595,7 @@ def walk(expr): if model._objective is not None: walk(model._objective.expression) - for c in model._constraints: + for c in model._constraints: # must be complete; see row_scan_is_anytime walk(c.body) if not complete[0]: @@ -2514,8 +2670,12 @@ class DefinedVariableForwardRule(NonlinearBoundTighteningRule): """ name = "defined_variable_forward" + # Each row's inference stands alone (bounds derived from THIS constraint, + # intersected into the box), so a prefix of the rows is a looser-but-valid + # tightening -- safe to stop at a deadline mid-scan (#875). + row_scan_is_anytime = True - def tighten(self, model, flat_lb, flat_ub, metadata): + def tighten(self, model, flat_lb, flat_ub, metadata, deadline=None): from discopt._jax.convexity.interval import Interval from discopt._jax.convexity.interval_eval import evaluate_interval @@ -2535,7 +2695,7 @@ def tighten(self, model, flat_lb, flat_ub, metadata): if sz == 1: idx_to_var[off] = var - for con in model._constraints: + for con in _rows_until(model, deadline): if getattr(con, "sense", None) != "==": continue rhs = float(getattr(con, "rhs", 0.0) or 0.0) @@ -2667,8 +2827,38 @@ def tighten_nonlinear_bounds( flat_ub: np.ndarray, rules: Sequence[NonlinearBoundTighteningRule] = DEFAULT_NONLINEAR_BOUND_RULES, max_rounds: int = 5, + deadline: Optional[float] = None, ) -> tuple[np.ndarray, np.ndarray, NonlinearBoundTighteningStats]: - """Run registered nonlinear tightening rules on a variable box.""" + """Run registered nonlinear tightening rules on a variable box. + + ``deadline`` is an absolute ``time.perf_counter()`` value bounding the whole + pass, or ``None`` (the default) for the unbounded legacy behaviour, which is + byte-identical to before — the polls are the only difference, and none of them + can fire. + + Why it exists (issue #875). This is up to ``max_rounds`` sweeps over every rule, + and every rule walks every constraint; nothing downstream bounded it, because it + runs *before* the first branch-and-bound node exists. Measured on + ``watercontamination0202`` (106,711 vars / 107,209 rows), the three root-setup + calls cost **80.9 s against a 30 s ``time_limit``**. Same class as the LP-engine + loops of #858 and the Rust presolve passes of #868, one layer up each time. + + The budget is honoured at three granularities, coarsest first: between rounds, + between rules, and — for a rule that declares ``row_scan_is_anytime`` — inside + its constraint scan (:func:`_rows_until`). The innermost one is what makes the + poll real rather than nominal: a single rule's sweep over 107k rows already + exceeds a tight budget, so a between-rules-only poll would be #868's ``probing`` + all over again. + + Soundness. Every early exit does strictly *less* work: rules only intersect the + incoming box, so a dropped round, a skipped rule, or an unscanned row leaves the + box **looser**, never wrong. A rule whose conclusion depends on a row being + *absent* is skipped whole rather than truncated (see + :attr:`NonlinearBoundTighteningRule.row_scan_is_anytime`). The one thing a + budget can cost is an infeasibility *proof* the full pass would have found — + ``infeasible`` stays unset and branch-and-bound discovers it instead, which is + weaker, never false. + """ tightened_lb = np.asarray(flat_lb, dtype=np.float64).copy() tightened_ub = np.asarray(flat_ub, dtype=np.float64).copy() initial_lb = tightened_lb.copy() @@ -2708,13 +2898,35 @@ def _count_tightened(lb: np.ndarray, ub: np.ndarray) -> int: # so they are explored and validated by the node NLP rather than pruned. _snap_tolerant_crossovers(tightened_lb, tightened_ub) + deadline_reached = False + + def _out_of_budget() -> bool: + return deadline is not None and time.perf_counter() >= deadline + for _ in range(max(1, int(max_rounds))): + if _out_of_budget(): + deadline_reached = True + break round_changed = False for rule in rules: + # Between-rules poll. A rule that is not ``row_scan_is_anytime`` cannot + # be cut mid-scan without risking an unsound conclusion, so this is the + # only place it can be declined — never started rather than truncated, + # which bounds the overrun to one in-flight rule. + if _out_of_budget(): + deadline_reached = True + break prev_lb = tightened_lb.copy() prev_ub = tightened_ub.copy() + rule_kwargs = {} + if ( + deadline is not None + and getattr(rule, "row_scan_is_anytime", False) + and _rule_accepts_deadline(rule) + ): + rule_kwargs["deadline"] = deadline try: - cand_lb, cand_ub = rule.tighten(model, prev_lb, prev_ub, metadata) + cand_lb, cand_ub = rule.tighten(model, prev_lb, prev_ub, metadata, **rule_kwargs) except NonlinearBoundTighteningInfeasible as exc: _mark_rule(rule.name) return ( @@ -2725,6 +2937,7 @@ def _count_tightened(lb: np.ndarray, ub: np.ndarray) -> int: applied_rules=tuple(applied_rules), infeasible=True, infeasibility_reason=str(exc), + deadline_reached=deadline_reached, ), ) @@ -2746,6 +2959,7 @@ def _count_tightened(lb: np.ndarray, ub: np.ndarray) -> int: infeasibility_reason=( f"{rule.name} returned an empty interval for flat variable {first_idx}" ), + deadline_reached=deadline_reached, ), ) # Snap sub-tolerance crossovers introduced by intersecting this rule's @@ -2759,7 +2973,7 @@ def _count_tightened(lb: np.ndarray, ub: np.ndarray) -> int: if n_changed > 0: _mark_rule(rule.name) round_changed = True - if not round_changed: + if deadline_reached or not round_changed: break return ( @@ -2768,5 +2982,6 @@ def _count_tightened(lb: np.ndarray, ub: np.ndarray) -> int: NonlinearBoundTighteningStats( n_tightened=_count_tightened(tightened_lb, tightened_ub), applied_rules=tuple(applied_rules), + deadline_reached=deadline_reached, ), ) diff --git a/python/discopt/_jax/primal_heuristics.py b/python/discopt/_jax/primal_heuristics.py index f2eb22d8..ba95db5c 100644 --- a/python/discopt/_jax/primal_heuristics.py +++ b/python/discopt/_jax/primal_heuristics.py @@ -1951,7 +1951,7 @@ def _detect_one_hot_groups(model: Model, binary_mask: np.ndarray, n_vars: int) - Returns the list of groups (each a sorted list of flat binary indices, one entry per slot), or ``[]`` when no such structure is present. """ - from discopt._jax.milp_relaxation import _linearize_affine_expr + from discopt._jax.milp_relaxation import _linearize_affine_expr_sparse groups: list[list[int]] = [] seen: set[int] = set() @@ -1959,18 +1959,24 @@ def _detect_one_hot_groups(model: Model, binary_mask: np.ndarray, n_vars: int) - if getattr(c, "sense", None) != "==": continue try: - coeff, const = _linearize_affine_expr(c.body, model, n_vars) + # Sparse: this scan touches EVERY constraint, and the dense + # linearization costs O(n_vars) per row to allocate and zero — the + # #875 shape (~460 s of root setup on a 106,711-var instance from the + # identical pattern in ``_fix_single_var_equalities``). + terms, const = _linearize_affine_expr_sparse(c.body, model, n_vars) except Exception as exc: # noqa: BLE001 - a non-affine row is not a one-hot row logger.debug("one-hot row scan skipped a body: %s: %s", type(exc).__name__, exc) continue if not np.isfinite(const) or abs(float(const) + 1.0) > 1e-9: continue # not ``... == 1`` - nz = np.nonzero(np.abs(coeff) > 1e-9)[0] - if nz.size < 2: + nz = sorted(j for j, v in terms.items() if abs(v) > 1e-9) + if len(nz) < 2: continue - if not np.all(np.abs(coeff[nz] - 1.0) <= 1e-9): + if nz[0] < 0 or nz[-1] >= n_vars: + continue # out of the flat range the dense array bounded by raising + if any(abs(terms[j] - 1.0) > 1e-9 for j in nz): continue # non-unit coefficients — not a plain one-hot sum - if nz.max() >= binary_mask.size or not np.all(binary_mask[nz]): + if nz[-1] >= binary_mask.size or not np.all(binary_mask[nz]): continue # support is not entirely binary g = [int(i) for i in nz] if seen.intersection(g): diff --git a/python/discopt/_jax/problem_classifier.py b/python/discopt/_jax/problem_classifier.py index 130dad65..a494f285 100644 --- a/python/discopt/_jax/problem_classifier.py +++ b/python/discopt/_jax/problem_classifier.py @@ -953,6 +953,15 @@ def _extract_constraints_algebraic(model: Model, n_orig: int): def _quadratic_row_has_terms(Q: np.ndarray, tol: float = 1e-12) -> bool: + """True when ``Q`` holds a nonzero entry above ``tol``. + + Sparse-aware (#875): ``np.abs`` on a scipy sparse matrix returns a sparse matrix, + and ``np.any`` on one does not mean what it does on an ndarray — the stored + values are the only candidates, so test those directly. Explicit zeros can be + stored, hence the ``> tol`` test rather than ``nnz``. + """ + if _sp_issparse(Q): + return bool(np.any(np.abs(Q.data) > tol)) return bool(np.any(np.abs(Q) > tol)) @@ -1512,7 +1521,34 @@ def _extract_qp_data_from_repr(model: Model) -> QPData: def _extract_quadratic_coefficients_from_values(evaluate, n_vars: int): - """Extract 0.5*x'Q*x + c'x + d from a quadratic scalar evaluator.""" + """Extract 0.5*x'Q*x + c'x + d from a quadratic scalar evaluator. + + Same probe identities as :func:`_extract_qp_data_from_repr`: + - ``d = f(0)`` + - ``Q[j,j] = f(e_j) + f(-e_j) - 2d`` + - ``Q[i,j] = f(e_i+e_j) - f(e_i) - f(e_j) + d`` (``i != j``) + - ``c_j = f(e_j) - d - 0.5*Q[j,j]`` + + and now the same two #863 economies, which this (the QCP/QCQP probe path) was + left out of because ``watercontamination0202`` does not route here — #868 + declined to widen speculatively, and #875 closes the gap rather than leave two + extractors of one shape with two different cost profiles: + + * **Off-diagonal probing is restricted to the evaluator's support.** A variable + absent from this row has ``f(e_j) == f(-e_j) == d`` and a zero diagonal, so + every product involving it is identically zero. The support falls out of the + ``O(n)`` diagonal probes already taken, so the restriction is free, and it + takes the pair sweep from ``O(n^2)`` probes to ``O(|support|^2)``. + * **Q is materialised through :func:`_materialise_Q`** — dense while ``(n, n)`` + float64 fits the budget (bit-identical to the ``np.zeros((n, n))`` it + replaces, same entries written into the same zeros), scipy CSR beyond it. A + dense ``(n, n)`` is 91 GB at ``n = 106,711``, and this function was called + once per constraint. + + Both are pure cost reductions: the entries not probed are provably zero, and an + entry absent from the accumulator densifies back to the 0.0 the dense array + already held. + """ x_zero = np.zeros(n_vars, dtype=np.float64) d = float(evaluate(x_zero)) @@ -1526,23 +1562,26 @@ def _extract_quadratic_coefficients_from_values(evaluate, n_vars: int): ej[j] = -1.0 f_neg_ej[j] = float(evaluate(ej)) - Q = np.zeros((n_vars, n_vars), dtype=np.float64) + diag = f_ej + f_neg_ej - 2.0 * d + + terms: dict[tuple[int, int], float] = {} for j in range(n_vars): - Q[j, j] = f_ej[j] + f_neg_ej[j] - 2.0 * d + if diag[j] != 0.0: + terms[(j, j)] = float(diag[j]) - for i in range(n_vars): - for j in range(i + 1, n_vars): + support = [j for j in range(n_vars) if f_ej[j] != d or f_neg_ej[j] != d or diag[j] != 0.0] + for _si, i in enumerate(support): + for j in support[_si + 1 :]: eij = np.zeros(n_vars, dtype=np.float64) eij[i] = 1.0 eij[j] = 1.0 - f_eij = float(evaluate(eij)) - qij = f_eij - f_ej[i] - f_ej[j] + d - Q[i, j] = qij - Q[j, i] = qij + qij = float(evaluate(eij)) - f_ej[i] - f_ej[j] + d + if qij != 0.0: + terms[(i, j)] = qij + terms[(j, i)] = qij - c_vec = np.zeros(n_vars, dtype=np.float64) - for j in range(n_vars): - c_vec[j] = f_ej[j] - d - 0.5 * Q[j, j] + Q = _materialise_Q(terms, n_vars) + c_vec = f_ej - d - 0.5 * diag return Q, c_vec, d @@ -1578,7 +1617,9 @@ def _extract_qcp_data_from_repr(model: Model) -> QCPData: if _quadratic_row_has_terms(row_Q): q_rows.append( QuadraticConstraintData( - Q=np.asarray(row_Q), # type: ignore[arg-type] + # Preserve sparsity (see dense_Q / #863): np.asarray() on a scipy + # sparse matrix returns a 0-d object array instead of raising. + Q=row_Q if _sp_issparse(row_Q) else np.asarray(row_Q), # type: ignore[arg-type] c=np.asarray(row_c), # type: ignore[arg-type] sense=sense, rhs=rhs, @@ -1607,7 +1648,8 @@ def _extract_qcp_data_from_repr(model: Model) -> QCPData: obj_const = -obj_const return QCPData( - Q=np.asarray(Q), # type: ignore[arg-type] + # Preserve sparsity (see dense_Q / #863, #875). + Q=Q if _sp_issparse(Q) else np.asarray(Q), # type: ignore[arg-type] c=np.asarray(c_vec), # type: ignore[arg-type] # Preserve sparsity (see dense_A / #863). A_ub=A_ub if _sp_issparse(A_ub) else np.asarray(A_ub), # type: ignore[arg-type] diff --git a/python/discopt/_jax/uniform_relax.py b/python/discopt/_jax/uniform_relax.py index 429d5215..c94521db 100644 --- a/python/discopt/_jax/uniform_relax.py +++ b/python/discopt/_jax/uniform_relax.py @@ -67,7 +67,7 @@ class of the R1.1 canonical DAG (``canonical_expr.atomize`` taxonomy) soundly an _expression_lower_bound_for_lift, _flat_variable_types, _integer_domain_values, - _linearize_affine_expr, + _linearize_affine_expr_sparse, ) from discopt._jax.model_utils import flat_variable_bounds from discopt.modeling.core import Model, ObjectiveSense, UnaryOp @@ -2481,6 +2481,16 @@ def _fix_single_var_equalities( when the pinned value lies inside the current box; an out-of-box value is left for the equality's own LP rows to expose as infeasible (never silently widen or empty the box here). Returns fresh arrays; inputs are untouched. + + Cost (issue #875). This scan is ``O(sum of equality body sizes)``, NOT + ``O(n_constraints * n_vars)``: it reads the *sparse* affine linearization, whose + dict carries one entry per referenced variable. The dense + ``_linearize_affine_expr`` allocates and zeroes an ``n_vars`` array per call, and + walking that array in Python to find the single nonzero costs ``n_vars`` again — + on ``watercontamination0202`` (106,711 vars / 107,209 rows) that was **~460 s of + a 30 s solve budget**, 23 of 29 stack samples, for bodies with ONE leaf each. + The pass runs before the relaxation build and before any node exists, so nothing + downstream bounded it. """ lb = np.array(flat_lb, dtype=np.float64) ub = np.array(flat_ub, dtype=np.float64) @@ -2489,13 +2499,19 @@ def _fix_single_var_equalities( if getattr(con, "sense", None) != "==": continue try: - coeff, const = _linearize_affine_expr(con.body, model, n) + terms, const = _linearize_affine_expr_sparse(con.body, model, n) except (ValueError, TypeError): continue - nz = [(j, float(c)) for j, c in enumerate(np.asarray(coeff)) if abs(float(c)) > 1e-12] + nz = [(j, c) for j, c in terms.items() if abs(c) > 1e-12] if len(nz) != 1: continue j, c = nz[0] + # The sparse core does not range-check its keys; the dense wrapper's array + # store did, by raising an (uncaught) IndexError. Both are unreachable for a + # box covering the model's own flat variables, but a pin we cannot index is + # a pin we must not apply — skip rather than write outside the box. + if not (0 <= j < n): + continue val = (float(con.rhs) - float(const)) / c if not math.isfinite(val): continue diff --git a/python/discopt/modeling/core.py b/python/discopt/modeling/core.py index 277f6a13..18c80059 100644 --- a/python/discopt/modeling/core.py +++ b/python/discopt/modeling/core.py @@ -1673,6 +1673,14 @@ def __repr__(self): # Solve Result # ───────────────────────────────────────────────────────────── +#: The one status whose ``gap_certified=True`` certifies something OTHER than an +#: optimality gap — an infeasibility proof legitimately carries neither a dual bound +#: nor an incumbent. Every other status must produce both ends of a gap to claim one; +#: see :meth:`SolveResult.__post_init__`. Deliberately not widened to ``unbounded``: +#: that would *weaken* an existing guard, and an unbounded minimize already carries +#: ``bound=-inf``, so it is decertified by the dual-side check regardless. +_NON_GAP_CERTIFICATE_STATUSES = frozenset({"infeasible"}) + @dataclass class SolveResult: @@ -1854,15 +1862,40 @@ def __post_init__(self) -> None: # global lower bound at ``-inf``. Reporting ``gap_certified=True`` there # is a false certification (the benchmark gate would miscount it as a # solved/certified instance), so we downgrade it and clear the - # meaningless bound/gap. Infeasibility certificates are exempt: - # ``status="infeasible"`` with ``gap_certified=True`` certifies - # infeasibility, not a gap, and legitimately carries ``bound=None``. - if self.gap_certified and self.status != "infeasible": + # meaningless bound/gap. Infeasibility certificates are exempt + # (``_NON_GAP_CERTIFICATE_STATUSES``): ``status="infeasible"`` with + # ``gap_certified=True`` certifies infeasibility, not a gap, and + # legitimately carries ``bound=None``. + if self.gap_certified and self.status not in _NON_GAP_CERTIFICATE_STATUSES: if self.bound is None or not np.isfinite(self.bound): self.gap_certified = False self.bound = None self.gap = None + # Same guard, other end of the gap (#875). A gap has two ends, and the check + # above only requires the dual one: a ``time_limit`` exit that never found an + # incumbent came back ``objective=None, gap=None, gap_certified=True``, which + # claims a certified gap where no gap was ever formed. Harmless where nothing + # is reported at all, but the flag is exactly what a consumer checks *before* + # reading the values — the graduation panels count a ``gap_certified=True`` + # instance as certified, and the "no certification regression" rule compares + # that flag across a flag flip. Reproduced with a 0.5 s budget on a 300-var + # bilinear model: ``status=time_limit objective=None bound=-7497.0 gap=None + # gap_certified=True nodes=0``. + # + # This only ever downgrades True -> False, so it cannot manufacture a + # certificate; and it leaves ``bound`` alone, because a dual bound with no + # incumbent is still a perfectly valid dual bound worth reporting — it is the + # *gap* claim that was unfounded. Same exemption as above: a status whose + # certificate is not a gap (``infeasible``) is not asked to produce one. + if ( + self.gap_certified + and self.objective is None + and self.status not in _NON_GAP_CERTIFICATE_STATUSES + ): + self.gap_certified = False + self.gap = None + def value(self, var: Variable) -> np.ndarray: """ Get the optimal value of a variable. diff --git a/python/discopt/solver.py b/python/discopt/solver.py index 4d8513d9..fca9639b 100644 --- a/python/discopt/solver.py +++ b/python/discopt/solver.py @@ -2372,7 +2372,25 @@ def _apply_nonlinear_tightening_with_status( try: from discopt._jax.nonlinear_bound_tightening import tighten_nonlinear_bounds - tightened_lb, tightened_ub, stats = tighten_nonlinear_bounds(model, lb, ub) + # #875: this was the ONE nonlinear-tightening entry point with no budget. The + # other three all pass a deadline — the root declared-box pass, the + # periodic/domain pass, and AMP — so this one ran the pass to completion on + # every invocation. On watercontamination0202 (107k rows) that is ~23 s per + # call and it fired 3x inside a 30 s ``time_limit``, i.e. ~70 s of a 126 s + # profile: the single largest remaining overrun after the sparse-linearizer + # fix, and the reason a deadline on ``tighten_nonlinear_bounds`` alone did not + # bound the solve. + # + # Bound it by the solve's own ABSOLUTE deadline rather than a fresh per-call + # fraction. A per-call fraction is what let the convexity classifier's budget + # multiply across model objects (each reformulation restarted it), and this + # helper is called per node as well as at the root, so a fraction here would + # compound the same way. ``None`` (no time limit) keeps the unbounded pass, + # which is the current behavior for an untimed solve. + _nbt_deadline = getattr(model, "_solve_deadline", None) + tightened_lb, tightened_ub, stats = tighten_nonlinear_bounds( + model, lb, ub, deadline=_nbt_deadline + ) except Exception as exc: logger.debug("Skipping nonlinear tightening after error: %s", exc) return lb, ub, False @@ -3315,13 +3333,20 @@ def _format_bad_bound_entries( return bad_vars -def _declared_box_tightening(model: Model): +def _declared_box_tightening(model: Model, deadline: Optional[float] = None): """Run ``tighten_nonlinear_bounds`` **once** on the model's declared box. Returns ``(tightened_lb, tightened_ub, stats)``, or ``None`` when the pass raised (both consumers below degrade to "no information", exactly as they did when each caught the exception itself). + ``deadline`` (an absolute ``time.perf_counter()``) bounds the pass; ``None`` + leaves it unbudgeted, which is what a standalone caller outside ``solve_model`` + gets. Sharing the result (#863) halved this phase but did not bound it: the + remaining single call still cost ~27 s against a 30 s ``time_limit`` on + ``watercontamination0202``, because the pass had no budget parameter at all + (#875). A truncated pass returns a looser — never wrong — box. + Why this exists (#863): ``_check_finite_bounds`` and ``_detect_nonlinear_bound_infeasibility`` are called back-to-back in ``solve_model`` on an *unmodified* model, and each ran the whole pass over the @@ -3346,7 +3371,7 @@ def _declared_box_tightening(model: Model): try: from discopt._jax.nonlinear_bound_tightening import tighten_nonlinear_bounds - return tighten_nonlinear_bounds(model, raw_lb, raw_ub) + return tighten_nonlinear_bounds(model, raw_lb, raw_ub, deadline=deadline) except Exception as exc: logger.debug("Nonlinear bound tightening on the declared box failed: %s", exc) return None @@ -3585,6 +3610,16 @@ def _classify_model_convexity( overrun that budget it is abandoned and the model is reported as convexity-unknown, which routes to the sound spatial Branch and Bound. This keeps a tight ``time_limit`` from being blown by classification alone. + + That per-call budget alone does not bound the *solve*, because the memo is keyed + to a model object and a reformulation produces a new one: each classification + then gets a fresh ``0.2 * time_limit``, and they add up (#875 measured two runs + at 14.7 s under a 30 s ``time_limit`` on ``watercontamination0202``, after which + half the budget was already gone). So the deadline is additionally clamped to the + absolute ``model._solve_deadline`` — classification, however many times it runs, + can never on its own carry the solve past ``time_limit``. Hitting the clamp is + sound: an abandoned classification reports convexity-unknown, which is the + conservative verdict that routes to spatial Branch and Bound. """ cached = getattr(model, "_convexity_classification_cache", None) if cached is not None: @@ -3596,6 +3631,9 @@ def _classify_model_convexity( # solve_model (e.g. on a model produced by factorable reformulation). budget = getattr(model, "_convexity_time_budget", 15.0) deadline = (time.perf_counter() + budget) if budget else None + solve_deadline = getattr(model, "_solve_deadline", None) + if solve_deadline is not None: + deadline = solve_deadline if deadline is None else min(deadline, float(solve_deadline)) result: tuple[bool, bool, list[bool] | None] try: from discopt._jax.convexity import classify_model as _classify_convexity @@ -3718,7 +3756,7 @@ def _apply_auto_cut_policy(model: "Model", relaxer) -> None: ``_psd_cuts`` / ``_rlt_cuts`` flags; purely a performance choice — every cut family is sound, so this never affects correctness. """ - from discopt._jax.milp_relaxation import _linear_constraint_forms + from discopt._jax.milp_relaxation import _any_linear_constraint_form try: n = sum(v.size for v in model._variables) @@ -3727,7 +3765,7 @@ def _apply_auto_cut_policy(model: "Model", relaxer) -> None: # product structure is sparse, past the raw variable-count gate. if n > _AUTO_CUTS_MAX_VARS and not _rlt_sparse_admit(model, n): return # size gate: leave cuts off - has_linear_constraints = bool(_linear_constraint_forms(model, n)) + has_linear_constraints = _any_linear_constraint_form(model, n) if has_linear_constraints: relaxer._rlt_cuts = True relaxer._psd_cuts = False @@ -5475,11 +5513,18 @@ def _note_ignored_gp_minlp(name: str, should_warn: bool) -> None: # to one period, and clamp log/sqrt arguments to their natural domain so # the local NLP never wanders into the undefined region (issue #265's # false-infeasible from a free log argument). + # Budgeted like the other root-setup passes (#875): both rules walk every + # constraint body, and this runs before a single node exists. + # ``PeriodicVariableBoundRule`` is not row-anytime (its conclusion rests on a + # variable being absent elsewhere), so the budget can only decline to start + # it — never truncate it — which is exactly what keeps the reduction sound. + _per_budget_s = min(min(max(0.05 * float(time_limit), 1.0), 10.0), _remaining_budget()) _per_lb, _per_ub, _per_stats = tighten_nonlinear_bounds( model, _origin_lb_chk, _origin_ub_chk, rules=(PeriodicVariableBoundRule(), FunctionDomainBoundRule()), + deadline=time.perf_counter() + _per_budget_s, ) if _per_stats.n_tightened > 0: from discopt.solvers.amp import _apply_flat_bounds_to_model @@ -5570,6 +5615,10 @@ def _note_ignored_gp_minlp(name: str, should_warn: bool) -> None: model = _bml model._convexity_classification_cache = None model._convexity_time_budget = _convexity_time_budget + # Carry the absolute deadline onto the reformulated model too, + # or its classification (and the engines that read this stash) + # would run against no wall clamp at all (#654 / #875). + model._solve_deadline = _solve_t0 + float(time_limit) # Route like the integer-bilinear reform: skip the (slow, # redundant) FBBT root presolve on the lifted rows and use # the monolithic Rust simplex MILP engine, unless the @@ -6199,7 +6248,15 @@ class _DcbSkip(Exception): # on an unmodified model, so run it once and share it (#863): measured on # watercontamination0202 the second run was 39.78 s of bit-identical repeat work # against a 30 s budget. See _declared_box_tightening. - _declared_tightening = _declared_box_tightening(model) + # + # #875: sharing halved the phase but left it unbounded — the surviving call was + # still ~27 s of a 30 s time_limit, because the pass had no budget parameter. + # Same shape as the root presolve budget below: a share of ``time_limit``, + # clamped to what is actually left. Truncation only weakens the box. + _nbt_budget_s = min(min(max(0.15 * float(time_limit), 2.0), 30.0), _remaining_budget()) + _declared_tightening = _declared_box_tightening( + model, deadline=time.perf_counter() + _nbt_budget_s + ) _check_finite_bounds(model, _declared_tightening) nonlinear_infeasibility = _detect_nonlinear_bound_infeasibility(model, _declared_tightening) if nonlinear_infeasibility is not None: @@ -14580,7 +14637,11 @@ def _solve_qcp_gurobi( integrality = int_arr result = _gurobi_solve_qcp( - Q=np.asarray(qcp_data.Q[:n_orig, :n_orig]), + # ``_dense_Q``, not ``np.asarray``: the QCP extractor may now emit a scipy + # sparse Q (#875, matching the QP extractor since #863), and ``np.asarray`` + # on one returns a 0-d object array rather than raising — it would smuggle + # garbage into the Gurobi call instead of failing loudly. + Q=_dense_Q(qcp_data.Q)[:n_orig, :n_orig], c=np.asarray(qcp_data.c[:n_orig]), A_ub=A_ub_arg, b_ub=b_ub_arg, diff --git a/python/discopt/solvers/amp.py b/python/discopt/solvers/amp.py index 691846b3..00144e2a 100644 --- a/python/discopt/solvers/amp.py +++ b/python/discopt/solvers/amp.py @@ -2403,8 +2403,16 @@ def _from_minimization_space(value: float) -> float: if root_changed: logger.info("AMP: root FBBT tightened variable bounds before relaxation") + # Budgeted root-setup pass (#875): every rule walks every constraint, and this + # runs before AMP's first relaxation. A share of ``time_limit`` clamped to the + # wall actually left; truncation drops tightenings, so the box comes back looser + # but never wrong. + _nbt_budget_s = min( + min(max(0.15 * float(time_limit), 2.0), 30.0), + max(0.0, float(time_limit) - (time.perf_counter() - t_start)), + ) tightened_lb, tightened_ub, nonlinear_bt_stats = tighten_nonlinear_bounds( - model, flat_lb, flat_ub + model, flat_lb, flat_ub, deadline=time.perf_counter() + _nbt_budget_s ) if nonlinear_bt_stats.infeasible: logger.info( diff --git a/python/tests/test_863_shared_declared_box_tightening.py b/python/tests/test_863_shared_declared_box_tightening.py index 4fda78f1..ef534842 100644 --- a/python/tests/test_863_shared_declared_box_tightening.py +++ b/python/tests/test_863_shared_declared_box_tightening.py @@ -158,9 +158,12 @@ def test_solve_runs_the_declared_box_pass_once(monkeypatch): real = _solver._declared_box_tightening calls: list[int] = [] - def _counting(model): + # ``*a, **kw`` rather than ``(model)``: ``solve_model`` now also passes the pass's + # ``deadline`` (#875), and a fixed-arity spy would turn a plumbing change into a + # TypeError inside the solve rather than a count mismatch here. + def _counting(*a, **kw): calls.append(1) - return real(model) + return real(*a, **kw) monkeypatch.setattr(_solver, "_declared_box_tightening", _counting) with pytest.warns(UserWarning, match="very large or infinite declared bounds"): diff --git a/python/tests/test_875_nbt_call_site_budget.py b/python/tests/test_875_nbt_call_site_budget.py new file mode 100644 index 00000000..8e05d27d --- /dev/null +++ b/python/tests/test_875_nbt_call_site_budget.py @@ -0,0 +1,157 @@ +"""#875: every ``tighten_nonlinear_bounds`` call site must carry a budget. + +``tighten_nonlinear_bounds`` grew a ``deadline`` parameter, and three of its four +call sites pass one — the root declared-box pass, the periodic/domain pass, and +AMP. ``_apply_nonlinear_tightening_with_status`` did not, so it ran the pass to +completion on every invocation. + +Measured on ``watercontamination0202`` (106,711 vars / 107,209 rows), the instance +#875 is about: that unbudgeted site cost ~23 s per call and fired 3x inside a 30 s +``time_limit`` — about 70 s of a 126 s profile, the single largest remaining overrun +after the sparse-linearizer fix, and the reason a ``deadline`` on +``tighten_nonlinear_bounds`` alone did not bound the solve. End to end: + +====================================== ========== ========== +stage T=30 T=60 +====================================== ========== ========== +before #878 579.3 s 620.8 s +#878 (sparse linearizer + nbt poll) 88.5 s 118.0 s +plus this call site 47.2 s 76.3 s +====================================== ========== ========== + +This test pins the wiring rather than the wall-clock, so it cannot rot into a +machine-speed assertion: it asserts the deadline actually *arrives* at +``tighten_nonlinear_bounds``, and that it equals the solve's absolute deadline +rather than a fresh per-call fraction (a fraction recomputed per call is how the +convexity classifier's budget used to multiply across model objects, and this +helper runs per node as well as at the root). +""" + +from __future__ import annotations + +import discopt.modeling as dm +import numpy as np +import pytest +from discopt import solver as solver_mod + + +def _tiny_nonlinear_model(): + m = dm.Model("nbt_budget") + x = m.continuous("x", lb=0.5, ub=3.0) + y = m.continuous("y", lb=0.5, ub=3.0) + m.minimize(x + y) + m.subject_to(x * y >= 1.0) + return m + + +def test_call_site_forwards_the_solve_deadline(monkeypatch): + """The deadline reaches ``tighten_nonlinear_bounds``, and it is the absolute one.""" + seen: list[float | None] = [] + real = solver_mod._apply_nonlinear_tightening_with_status + + import discopt._jax.nonlinear_bound_tightening as nbt + + orig = nbt.tighten_nonlinear_bounds + + def spy(model, lb, ub, *args, **kwargs): + seen.append(kwargs.get("deadline", "MISSING")) # type: ignore[arg-type] + return orig(model, lb, ub, *args, **kwargs) + + monkeypatch.setattr(nbt, "tighten_nonlinear_bounds", spy) + + m = _tiny_nonlinear_model() + sentinel = 12345.0 + m._solve_deadline = sentinel + lb = np.array([0.5, 0.5]) + ub = np.array([3.0, 3.0]) + real(m, lb, ub) + + assert seen, "tighten_nonlinear_bounds was never called — the test is vacuous" + assert "MISSING" not in seen, ( + "_apply_nonlinear_tightening_with_status called tighten_nonlinear_bounds with no " + "deadline= argument; this call site is unbudgeted (#875)" + ) + assert seen[0] == sentinel, ( + f"expected the solve's absolute deadline {sentinel!r}, got {seen[0]!r} — a per-call " + "fraction here compounds across the per-node invocations" + ) + + +def test_absent_solve_deadline_leaves_the_pass_unbounded(monkeypatch): + """No ``time_limit`` must keep today's behavior: an unbounded pass, not a zero budget. + + Guards the direction of the fix: reading a missing attribute as ``0.0`` would + make every untimed solve skip tightening entirely. + """ + seen: list[object] = [] + import discopt._jax.nonlinear_bound_tightening as nbt + + orig = nbt.tighten_nonlinear_bounds + + def spy(model, lb, ub, *args, **kwargs): + seen.append(kwargs.get("deadline", "MISSING")) + return orig(model, lb, ub, *args, **kwargs) + + monkeypatch.setattr(nbt, "tighten_nonlinear_bounds", spy) + + m = _tiny_nonlinear_model() + assert not hasattr(m, "_solve_deadline") + solver_mod._apply_nonlinear_tightening_with_status( + m, np.array([0.5, 0.5]), np.array([3.0, 3.0]) + ) + assert seen, "vacuous: tighten_nonlinear_bounds never called" + assert seen[0] is None, f"expected deadline=None for an untimed solve, got {seen[0]!r}" + + +def test_tightening_still_tightens_under_a_generous_budget(): + """The budget must not cost the tightening it is meant to bound. + + ``x*y >= 1`` on ``[0.5,3]^2`` is tightenable; with a far-future deadline the pass + must still return a box no wider than it was given (soundness is one-directional + here: tightening may only shrink). + """ + m = _tiny_nonlinear_model() + import time + + m._solve_deadline = time.perf_counter() + 3600.0 + lb = np.array([0.5, 0.5]) + ub = np.array([3.0, 3.0]) + new_lb, new_ub, infeasible = solver_mod._apply_nonlinear_tightening_with_status(m, lb, ub) + assert not infeasible + assert np.all(new_lb >= lb - 1e-12), "tightening widened a lower bound" + assert np.all(new_ub <= ub + 1e-12), "tightening widened an upper bound" + assert new_lb.shape == lb.shape and new_ub.shape == ub.shape + + +def test_expired_budget_is_sound_and_does_not_widen(): + """An already-expired deadline must return a box that is still valid, never wider.""" + m = _tiny_nonlinear_model() + m._solve_deadline = 1.0 # far in the past for perf_counter + lb = np.array([0.5, 0.5]) + ub = np.array([3.0, 3.0]) + new_lb, new_ub, infeasible = solver_mod._apply_nonlinear_tightening_with_status(m, lb, ub) + assert not infeasible, "an expired budget must not manufacture an infeasibility" + assert np.all(new_lb >= lb - 1e-12) and np.all(new_ub <= ub + 1e-12) + + +@pytest.mark.parametrize("attr", ["_apply_nonlinear_tightening_with_status"]) +def test_no_other_unbudgeted_call_site_regresses(attr): + """Source-level guard: no ``tighten_nonlinear_bounds(`` call in solver.py may omit + ``deadline``. Cheap, and it catches a new call site added without a budget.""" + import inspect + import re + + src = inspect.getsource(solver_mod) + calls = [m.start() for m in re.finditer(r"tighten_nonlinear_bounds\(", src)] + assert calls, "vacuous: no call sites found in solver.py" + checked = 0 + for pos in calls: + window = src[pos : pos + 400] + # the import line itself is not a call + if window.startswith("tighten_nonlinear_bounds(") and "deadline" not in window: + raise AssertionError( + f"a tighten_nonlinear_bounds call in solver.py has no deadline= within 400 " + f"chars: ...{window[:120]}..." + ) + checked += 1 + assert checked == len(calls) diff --git a/python/tests/test_875_root_setup_budget.py b/python/tests/test_875_root_setup_budget.py new file mode 100644 index 00000000..31773318 --- /dev/null +++ b/python/tests/test_875_root_setup_budget.py @@ -0,0 +1,488 @@ +"""#875: root setup must be *cheap* and *bounded* before the first B&B node. + +Successor to #863/#868. After PR #868, ``watercontamination0202`` (106,711 vars / +107,209 rows, 7 binaries) returned instead of hanging, but ``solve(time_limit=30)`` +took **579.3 s** — 19.3x — with ``nodes=0``. The time was fully attributed: + + phase cost note + _fix_single_var_equalities ~460 s dominant, 23/29 stack samples + tighten_nonlinear_bounds x3 80.9 s no deadline awareness at all + _classify_model_convexity x2 14.7 s per-model budget, not per-solve + presolve + load + classification ~12 s already capped by #868 + +Every item is pre-B&B root setup — the class #858 fixed in the LP engine and #868 +fixed in the Rust presolve passes, one layer up each time. Two distinct defects hide +under that one heading, and the tests below separate them: + +* a **cost** bug — ``_fix_single_var_equalities`` was ``O(n_constraints * n_vars)`` + because the affine linearizer it calls allocates and zeroes a dense ``n_vars`` + array per call, and the caller then walked all ``n_vars`` entries in Python to find + the single nonzero. That is not a budget problem; no deadline should have been + needed for a scan over one-leaf bodies. Fixed by a sparse linearizer core. +* a **budget** bug — ``tighten_nonlinear_bounds`` had no ``deadline`` parameter, and + the convexity classifier's budget is a fraction of ``time_limit`` recomputed per + *model object*, so a reformulation restarts it and the fractions add up. + +The in-repo corpus is far too small to show either as a wall-clock overrun, so these +tests target the mechanism directly: the cost bug via a scaling probe that fails +before the fix, the budget bugs via an expired deadline. +""" + +from __future__ import annotations + +import os + +os.environ.setdefault("JAX_PLATFORMS", "cpu") +os.environ.setdefault("JAX_ENABLE_X64", "1") + +import time # noqa: E402 +from pathlib import Path # noqa: E402 + +import discopt.modeling as dm # noqa: E402 +import numpy as np # noqa: E402 +import pytest # noqa: E402 +from discopt._jax.milp_relaxation import ( # noqa: E402 + _any_linear_constraint_form, + _linear_constraint_forms, + _linearize_affine_expr, + _linearize_affine_expr_sparse, +) +from discopt._jax.model_utils import flat_variable_bounds # noqa: E402 +from discopt._jax.nonlinear_bound_tightening import ( # noqa: E402 + DEFAULT_NONLINEAR_BOUND_RULES, + FunctionDomainBoundRule, + NonlinearBoundTighteningRule, + PeriodicVariableBoundRule, + tighten_nonlinear_bounds, +) +from discopt._jax.uniform_relax import _fix_single_var_equalities # noqa: E402 +from discopt.modeling.core import SolveResult # noqa: E402 + +# -------------------------------------------------------------------------- +# The sparse affine linearizer: same answer, no O(n_vars) per call +# -------------------------------------------------------------------------- + + +def _affine_model(n_vars: int, n_eq: int): + """``n_eq`` single-variable equalities over ``n_vars`` variables. + + Every equality body has exactly ONE leaf, so an honest scan is O(n_eq); the dense + linearizer made it O(n_eq * n_vars) regardless. + """ + m = dm.Model(f"aff{n_vars}") + x = m.continuous("x", shape=(n_vars,), lb=-10.0, ub=10.0) + for k in range(n_eq): + m.subject_to(x[k % n_vars] == 1.0) + m.minimize(x[0]) + return m + + +def test_sparse_linearizer_agrees_with_the_dense_one(): + """The dense wrapper is now a view of the sparse core; they must not diverge.""" + m = dm.Model("agree") + x = m.continuous("x", shape=(6,), lb=-1.0, ub=1.0) + y = m.continuous("y", lb=-1.0, ub=1.0) + m.minimize(x[0]) + + m.subject_to(2.0 * x[0] + 3.0 * x[3] - 4.0 * y == 1.0) + m.subject_to(x[1] - x[1] + x[5] <= 2.0) # cancelling terms -> an explicit zero + m.subject_to(dm.sum(x) + y >= 0.0) + + n = 7 + for con in m._constraints: + dense, dense_const = _linearize_affine_expr(con.body, m, n) + terms, sparse_const = _linearize_affine_expr_sparse(con.body, m, n) + assert sparse_const == dense_const + rebuilt = np.zeros(n, dtype=np.float64) + for j, c in terms.items(): + rebuilt[j] = c + assert np.array_equal(rebuilt, dense) + + +def test_sparse_linearizer_refuses_exactly_what_the_dense_one_refused(): + m = dm.Model("refuse") + x = m.continuous("x", shape=(2,), lb=-1.0, ub=1.0) + m.minimize(x[0]) + m.subject_to(x[0] * x[1] <= 1.0) # nonlinear body + body = m._constraints[0].body + with pytest.raises(ValueError): + _linearize_affine_expr(body, m, 2) + with pytest.raises(ValueError): + _linearize_affine_expr_sparse(body, m, 2) + + +def test_fix_single_var_equalities_is_flat_in_the_variable_count(): + """The decisive cost test, and the one that FAILS before the fix. + + At a fixed constraint count the pass must not get slower as variables are added: + the bodies do not change. Measured before the fix, 400 equalities: + + n_vars=2,000 0.068 s + n_vars=8,000 0.240 s (3.55x for 4x n_vars) + n_vars=32,000 0.951 s (3.96x) + n_vars=128,000 3.650 s (3.84x) + + i.e. exactly linear in ``n_vars`` — which extrapolates to the ~460 s the issue + measured. After: 0.001 s at n_vars=32,000 (0.002 s at 128,000), and flat. + + Two assertions, and the ABSOLUTE one is primary. A ratio of two post-fix timings + divides ~1 ms by ~1 ms, which on a shared CI runner is noise over noise — one + scheduler hiccup flips it. The absolute ceiling has ~250x headroom for the sparse + scan and still fails the dense one on a runner 3x slower than the machine these + numbers came from, so it measures the implementation rather than the box. The + ratio is kept as the complexity-class signal but only evaluated when the baseline + is far enough above the timer floor to mean anything. + + ``min`` over repetitions, not mean: the fastest run is the one least polluted by + whatever else the machine was doing, which is the whole question here. + """ + n_eq = 300 + reps = 3 + walls = {} + for n_vars in (2_000, 32_000): + best = float("inf") + for _ in range(reps): + m = _affine_model(n_vars, n_eq) + lb, ub = flat_variable_bounds(m) + t0 = time.perf_counter() + out_lb, out_ub = _fix_single_var_equalities(m, lb, ub) + best = min(best, time.perf_counter() - t0) + # the pass must still do its job: every pinned variable collapsed to a point + pinned = np.flatnonzero(out_lb == out_ub) + assert pinned.size >= min(n_eq, n_vars) + assert np.all(out_lb[pinned] == 1.0) + walls[n_vars] = best + + # Primary: dense was 0.851 s here at n_vars=32,000, sparse is ~0.001 s. + assert walls[32_000] < 0.25, ( + f"the scan is still paying O(n_vars) per row: {walls[32_000]:.3f}s at " + f"n_vars=32,000 with only {n_eq} rows (dense measured 0.851s, sparse 0.001s)" + ) + # Secondary: 16x the variables for the same rows must not cost 16x the wall. + # Skipped when the baseline is at the timer floor, where the quotient is noise. + if walls[2_000] > 5e-3: + ratio = walls[32_000] / walls[2_000] + assert ratio < 4.0, ( + f"cost still scales with n_vars (16x vars -> {ratio:.1f}x wall): {walls}" + ) + + +def test_fix_single_var_equalities_still_pins_and_still_refuses(): + """Behaviour parity on the cases the docstring promises.""" + m = dm.Model("pins") + x = m.continuous("x", shape=(4,), lb=-10.0, ub=10.0) + m.minimize(x[0]) + m.subject_to(2.0 * x[1] == 6.0) # pins x1 = 3 + m.subject_to(x[2] == 99.0) # OUTSIDE the box -> left for the LP rows + m.subject_to(x[0] + x[3] == 1.0) # two variables -> not a pin + lb, ub = flat_variable_bounds(m) + out_lb, out_ub = _fix_single_var_equalities(m, lb, ub) + + assert out_lb[1] == out_ub[1] == pytest.approx(3.0) + assert (out_lb[2], out_ub[2]) == (-10.0, 10.0), "an out-of-box pin must not be applied" + assert (out_lb[0], out_ub[0]) == (-10.0, 10.0) + assert (out_lb[3], out_ub[3]) == (-10.0, 10.0) + # inputs untouched + assert lb[1] == -10.0 and ub[1] == 10.0 + + +def test_any_linear_constraint_form_agrees_with_the_list(monkeypatch): + """The boolean probe must answer exactly what ``bool(_linear_constraint_forms())`` + answered, on both a model that has linear factors and one that has none.""" + m = _affine_model(4_000, 200) + assert _any_linear_constraint_form(m, 4_000) is True + assert bool(_linear_constraint_forms(m, 4_000)) is True + + nonlinear = dm.Model("nl") + z = nonlinear.continuous("z", shape=(3,), lb=0.1, ub=2.0) + nonlinear.minimize(z[0]) + nonlinear.subject_to(z[0] * z[1] <= 1.0) + nonlinear.subject_to(dm.log(z[2]) <= 1.0) + assert _any_linear_constraint_form(nonlinear, 3) is False + assert bool(_linear_constraint_forms(nonlinear, 3)) is False + + +def test_any_linear_constraint_form_short_circuits(monkeypatch): + """Counted, not timed: the probe must linearize ONE row when the first row is + linear, rather than every row. + + Counting the linearizations is the actual claim (``bool()`` of a fully + materialised list is what this replaced), and unlike a wall-clock comparison it + cannot flake on a loaded runner — the failure mode that cost real time on #863 + and briefly on this branch. + """ + import discopt._jax.milp_relaxation as mr + + calls = [] + real = mr._linearize_affine_expr_sparse + + def _counting(expr, model, n_vars): + calls.append(1) + return real(expr, model, n_vars) + + monkeypatch.setattr(mr, "_linearize_affine_expr_sparse", _counting) + + m = _affine_model(4_000, 200) + calls.clear() + assert _any_linear_constraint_form(m, 4_000) is True + short_circuit_calls = len(calls) + + calls.clear() + _linear_constraint_forms(m, 4_000) + full_calls = len(calls) + + assert short_circuit_calls == 1, ( + f"the boolean probe linearized {short_circuit_calls} rows; it must stop at " + f"the first linear one" + ) + assert full_calls == 200, f"expected one linearization per row, got {full_calls}" + + +# -------------------------------------------------------------------------- +# tighten_nonlinear_bounds honours a deadline +# -------------------------------------------------------------------------- + + +def _nbt_model(n: int = 40): + m = dm.Model(f"nbt{n}") + x = m.continuous("x", shape=(n,), lb=-100.0, ub=100.0) + y = m.continuous("y", shape=(n,)) # free: the rules exist to bound these + m.minimize(x[0]) + for k in range(n): + m.subject_to(x[k] * x[k] <= 4.0) + m.subject_to(y[k] - x[k] * x[k] == 0.0) + return m + + +def test_tightening_with_an_expired_deadline_does_nothing_and_says_so(): + m = _nbt_model() + lb, ub = flat_variable_bounds(m) + out_lb, out_ub, stats = tighten_nonlinear_bounds(m, lb, ub, deadline=time.perf_counter()) + assert stats.deadline_reached is True + assert stats.n_tightened == 0 + assert stats.applied_rules == () + assert np.array_equal(out_lb, lb) and np.array_equal(out_ub, ub) + assert stats.infeasible is False + + +def test_a_future_deadline_changes_nothing(): + """The no-regression half: when the budget is not binding — the normal case — the + poll must be the ONLY difference. Same box, same stats, byte for byte.""" + m = _nbt_model() + lb, ub = flat_variable_bounds(m) + base_lb, base_ub, base_stats = tighten_nonlinear_bounds(m, lb, ub) + dl_lb, dl_ub, dl_stats = tighten_nonlinear_bounds( + m, lb, ub, deadline=time.perf_counter() + 3600.0 + ) + assert np.array_equal(base_lb, dl_lb) + assert np.array_equal(base_ub, dl_ub) + assert base_stats.n_tightened == dl_stats.n_tightened + assert base_stats.applied_rules == dl_stats.applied_rules + assert base_stats.infeasible == dl_stats.infeasible + assert dl_stats.deadline_reached is False + assert base_stats.deadline_reached is False + + +def test_a_truncated_pass_only_ever_loosens(): + """Anytime contract: whatever a budgeted pass returns must be a SUPERSET of the + box an unbudgeted pass returns (looser or equal, never tighter, never wrong).""" + m = _nbt_model(60) + lb, ub = flat_variable_bounds(m) + full_lb, full_ub, _ = tighten_nonlinear_bounds(m, lb, ub) + for budget in (0.0, 1e-4, 1e-3, 1e-2): + cut_lb, cut_ub, _stats = tighten_nonlinear_bounds( + m, lb, ub, deadline=time.perf_counter() + budget + ) + assert np.all(cut_lb <= full_lb + 1e-12), "a budgeted pass tightened PAST the full pass" + assert np.all(cut_ub >= full_ub - 1e-12), "a budgeted pass tightened PAST the full pass" + assert np.all(cut_lb >= lb - 1e-12) and np.all(cut_ub <= ub + 1e-12) + + +def test_a_row_scan_that_needs_completeness_is_never_truncated(): + """``PeriodicVariableBoundRule`` concludes from what the model does NOT contain: a + variable is restricted to one period only while no row uses it outside sin/cos. A + truncated scan could miss the disqualifying row and cut feasible points, so the + rule must be skipped whole under a budget, never run on a prefix. + + Guarded structurally (the flag) and behaviourally (a model whose disqualifying + use sits in the LAST row). + """ + assert PeriodicVariableBoundRule.row_scan_is_anytime is False + assert NonlinearBoundTighteningRule.row_scan_is_anytime is False, ( + "the base-class default must stay the SAFE one, so a new accumulating rule " + "does not silently inherit truncation" + ) + assert FunctionDomainBoundRule.row_scan_is_anytime is True + + from discopt._jax.nonlinear_bound_tightening import _cached_flat_metadata + + def _build(disqualify: bool): + """``t`` free and used inside ``cos``; optionally ALSO used bare, in the very + last row, which is what makes the period reduction invalid.""" + m = dm.Model(f"periodic{int(disqualify)}") + m.continuous("t") # flat index 0 + pad = m.continuous("pad", shape=(400,), lb=0.0, ub=1.0) + m.minimize(dm.cos(m._variables[0])) + for k in range(400): + m.subject_to(pad[k] <= 1.0) + if disqualify: + m.subject_to(m._variables[0] + pad[0] <= 1e6) + return m + + expired = time.perf_counter() + rule = PeriodicVariableBoundRule() + + # Positive control: with nothing disqualifying it, the rule DOES reduce ``t`` to + # one period even under an expired deadline — so the scan really did run whole. + ok = _build(disqualify=False) + lb, ub = flat_variable_bounds(ok) + out_lb, out_ub = rule.tighten( + ok, lb.copy(), ub.copy(), _cached_flat_metadata(ok), deadline=expired + ) + assert (out_lb[0], out_ub[0]) == pytest.approx((-np.pi, np.pi)) + + # The case that matters: the disqualifying use is the LAST row, so any truncation + # would miss it and wrongly shrink ``t``. + bad = _build(disqualify=True) + lb, ub = flat_variable_bounds(bad) + out_lb, out_ub = rule.tighten( + bad, lb.copy(), ub.copy(), _cached_flat_metadata(bad), deadline=expired + ) + assert (out_lb[0], out_ub[0]) == (lb[0], ub[0]), ( + "t was restricted to one period despite a non-periodic use in the final row" + ) + + +def test_every_default_rule_declares_its_row_scan_contract(): + """A rule added without an explicit contract inherits ``False`` (skip-whole), + which is safe; this test exists so the choice is visible rather than accidental.""" + for rule in DEFAULT_NONLINEAR_BOUND_RULES: + assert isinstance(rule.row_scan_is_anytime, bool), rule.name + + +def test_a_legacy_four_argument_rule_still_works(): + """External rules keep the old signature; the runner must not pass ``deadline`` + to a rule that does not declare it.""" + + class LegacyRule(NonlinearBoundTighteningRule): + name = "legacy" + row_scan_is_anytime = True # even so: it cannot accept the kwarg + + def tighten(self, model, flat_lb, flat_ub, metadata): + del model, metadata + out_ub = flat_ub.copy() + out_ub[0] = min(float(out_ub[0]), 2.0) + return flat_lb.copy(), out_ub + + m = dm.Model("legacy") + x = m.continuous("x", lb=0.0, ub=10.0) + m.minimize(x) + lb = np.array([0.0]) + ub = np.array([10.0]) + out_lb, out_ub, stats = tighten_nonlinear_bounds( + m, lb, ub, rules=(LegacyRule(),), deadline=time.perf_counter() + 3600.0 + ) + assert out_ub[0] == pytest.approx(2.0) + assert stats.applied_rules == ("legacy",) + assert out_lb[0] == pytest.approx(0.0) + + +# -------------------------------------------------------------------------- +# gap_certified needs BOTH ends of the gap +# -------------------------------------------------------------------------- + + +def test_a_limit_exit_with_no_incumbent_is_not_certified(): + """``status=time_limit, objective=None, bound=, gap_certified=True`` claims + a certified gap where no gap was ever formed. The dual bound is kept — it is + valid — but the certification claim is dropped.""" + r = SolveResult(status="time_limit", objective=None, bound=-7497.0, gap_certified=True) + assert r.gap_certified is False + assert r.bound == pytest.approx(-7497.0), "a valid dual bound must survive the downgrade" + assert r.gap is None + + r2 = SolveResult(status="time_limit", objective=None, bound=None, gap_certified=True) + assert r2.gap_certified is False + + +def test_an_infeasibility_certificate_is_still_exempt(): + r = SolveResult(status="infeasible", objective=None, bound=None, gap_certified=True) + assert r.gap_certified is True + + +def test_a_real_certified_optimum_is_untouched(): + r = SolveResult(status="optimal", objective=1.5, bound=1.5, gap=0.0, gap_certified=True) + assert r.gap_certified is True + assert r.bound == pytest.approx(1.5) + assert r.gap == pytest.approx(0.0) + + +# -------------------------------------------------------------------------- +# end to end +# -------------------------------------------------------------------------- + + +@pytest.mark.slow +def test_a_wide_model_reaches_the_solver_inside_its_budget(): + """The issue's shape in miniature: many variables, many single-variable + equalities, a nonlinear core. Before the fix the root setup alone scaled with + ``n_vars * n_constraints`` and the solve never reached a node.""" + n = 4_000 + m = dm.Model("wide") + x = m.continuous("x", shape=(n,), lb=-5.0, ub=5.0) + b = m.binary("b", shape=(4,)) + for k in range(0, n, 2): + m.subject_to(x[k] == 0.5) + for k in range(0, 200, 2): + m.subject_to(x[k] * x[k + 1] <= 4.0) + m.subject_to(sum(b[j] for j in range(4)) == 2) + m.minimize(sum(x[k] for k in range(200)) + sum(b)) + + t0 = time.perf_counter() + res = m.solve(time_limit=20.0) + wall = time.perf_counter() - t0 + + assert wall < 60.0, f"root setup still dominates: {wall:.1f}s against a 20 s limit" + assert res.status in {"optimal", "feasible", "time_limit"} + if res.gap_certified: + assert res.objective is not None and res.bound is not None + assert res.bound <= res.objective + 1e-6 + + +_BIG = Path( + os.path.expanduser( + "~/Dropbox/projects/discopt-minlp-benchmark/minlplib/nl/watercontamination0202.nl" + ) +) + + +@pytest.mark.slow +@pytest.mark.skipif(not _BIG.exists(), reason="needs the full MINLPLib snapshot") +@pytest.mark.parametrize("budget", [30.0, 60.0]) +def test_watercontamination0202_honours_its_time_limit(budget): + """The issue's definition of done, on the instance that exposed the class. + + Before: ``solve(time_limit=30)`` = 579.3 s (19.3x) and ``solve(time_limit=60)`` = + 620.8 s (10.4x), both with ``nodes=0`` — the whole budget spent in root setup + before a single node existed. Required: within ~1.25x, with a sound result. + + Needs the full MINLPLib snapshot (the in-repo 61-file corpus has nothing near + this size), so it skips without it. The in-repo evidence for the same fixes is + the mechanism-level tests above: the scaling probe for the cost bug and the + expired-deadline tests for the budget bug. + """ + m = dm.from_nl(str(_BIG)) + t0 = time.perf_counter() + res = m.solve(time_limit=budget) + wall = time.perf_counter() - t0 + + assert wall < 1.25 * budget, ( + f"solve took {wall:.1f}s against a {budget:.0f}s time_limit " + f"({wall / budget:.1f}x); root setup is still unbounded" + ) + # Sound, whatever it managed to prove. + assert res.status != "infeasible", "FALSE-INFEASIBLE on a feasible instance" + if res.objective is not None and res.bound is not None: + assert res.bound <= res.objective + 1e-6, "UNSOUND CERT (bound > incumbent)" + if res.gap_certified: + assert res.objective is not None and res.bound is not None diff --git a/python/tests/test_875_sparse_qcp_probe_extraction.py b/python/tests/test_875_sparse_qcp_probe_extraction.py new file mode 100644 index 00000000..b19fb8f4 --- /dev/null +++ b/python/tests/test_875_sparse_qcp_probe_extraction.py @@ -0,0 +1,152 @@ +"""#875: the QCP *probe* extractor's Hessian may be sparse, and must agree exactly. + +``_extract_quadratic_coefficients_from_values`` is the numeric-probe counterpart of +``_extract_qp_data_from_repr`` and the sibling of the algebraic walk covered by +``test_863_sparse_algebraic_Q.py``. #863 made both of those support-restricted and +sparse; this one was left dense because ``watercontamination0202`` does not route +through it, and #868 declined to widen speculatively. It is the same shape, and it is +worse: ``_extract_qcp_data_from_repr`` calls it **once per constraint**, so the dense +``(n, n)`` and the all-pairs sweep are paid per row. + +Two independent economies, each tested here: + +* **support restriction** — a variable absent from the evaluator has + ``f(e_j) == f(-e_j) == d``, so every product involving it is identically zero. The + ``O(n)`` diagonal probes already identify the support, so the pair sweep drops from + ``O(n^2)`` to ``O(|support|^2)`` for free. +* **sparse materialisation** — through ``_materialise_Q``, dense while ``(n, n)`` + float64 fits the budget and scipy CSR beyond it. + +Every forced-sparse arm asserts ``sp.issparse``: without that, a sparse branch that +quietly failed would leave the comparison dense-against-dense and prove nothing — the +trap recorded in ``test_863_sparse_algebraic_Q.py``'s header. +""" + +from __future__ import annotations + +import os + +os.environ.setdefault("JAX_PLATFORMS", "cpu") +os.environ.setdefault("JAX_ENABLE_X64", "1") + +import discopt._jax.problem_classifier as pc # noqa: E402 +import numpy as np # noqa: E402 +import pytest # noqa: E402 +import scipy.sparse as sp # noqa: E402 +from discopt._jax.problem_classifier import ( # noqa: E402 + _extract_quadratic_coefficients_from_values, + _quadratic_row_has_terms, + dense_Q, +) + + +@pytest.fixture +def q_budget(monkeypatch): + def _set(nbytes): + monkeypatch.setattr(pc, "_QP_DENSE_Q_MAX_BYTES", nbytes) + + return _set + + +def _quadratic_evaluator(n: int, support: int = 5): + """``0.5 x'Qx + c'x + d`` touching only the first ``support`` variables. + + Returned as a plain callable so the test exercises the probe arithmetic rather + than the model plumbing; ``probes`` counts evaluator calls. + """ + rng = np.random.default_rng(7) + Q = np.zeros((n, n), dtype=np.float64) + block = rng.standard_normal((support, support)) + Q[:support, :support] = block + block.T # symmetric + c = np.zeros(n, dtype=np.float64) + c[:support] = rng.standard_normal(support) + d = 1.25 + probes = [0] + + def evaluate(x): + probes[0] += 1 + x = np.asarray(x, dtype=np.float64) + return 0.5 * float(x @ (Q @ x)) + float(c @ x) + d + + return evaluate, probes, Q, c, d + + +@pytest.mark.parametrize("n", [30, 60]) +def test_probe_extraction_recovers_the_quadratic(q_budget, n): + q_budget(10**12) + evaluate, _probes, Q, c, d = _quadratic_evaluator(n) + Q_out, c_out, d_out = _extract_quadratic_coefficients_from_values(evaluate, n) + assert not sp.issparse(Q_out), "a 10^12-byte budget should have stayed dense" + assert d_out == pytest.approx(d) + assert np.allclose(dense_Q(Q_out), Q, atol=1e-9) + assert np.allclose(np.asarray(c_out), c, atol=1e-9) + + +@pytest.mark.parametrize("n", [30, 60]) +def test_forced_sparse_arm_equals_the_dense_arm(q_budget, n): + """The whole safety argument: flipping the representation must not move an entry.""" + q_budget(10**12) + evaluate, _p, _Q, _c, _d = _quadratic_evaluator(n) + Q_dense, c_dense, d_dense = _extract_quadratic_coefficients_from_values(evaluate, n) + assert not sp.issparse(Q_dense) + + q_budget(1) + evaluate2, _p2, _Q2, _c2, _d2 = _quadratic_evaluator(n) + Q_sparse, c_sparse, d_sparse = _extract_quadratic_coefficients_from_values(evaluate2, n) + assert sp.issparse(Q_sparse), "a 1-byte budget should have forced a sparse Q" + + assert np.array_equal(dense_Q(Q_dense), dense_Q(Q_sparse)) + assert np.array_equal(np.asarray(c_dense), np.asarray(c_sparse)) + assert d_dense == d_sparse + + +def test_off_diagonal_probing_is_restricted_to_the_support(): + """Probe COUNT, not wall clock: the all-pairs sweep is ``n(n-1)/2`` off-diagonal + probes; restricted to a support of ``s`` it is ``s(s-1)/2``. This is the test that + fails before the fix (n=120, support=5: 7,140 probes vs 10).""" + n, support = 120, 5 + evaluate, probes, _Q, _c, _d = _quadratic_evaluator(n, support=support) + _extract_quadratic_coefficients_from_values(evaluate, n) + + fixed = 1 + 2 * n # f(0) plus f(e_j) / f(-e_j) for every j + off_diagonal = probes[0] - fixed + all_pairs = n * (n - 1) // 2 + support_pairs = support * (support - 1) // 2 + assert off_diagonal <= support_pairs, ( + f"off-diagonal probing is not support-restricted: {off_diagonal} probes " + f"(support allows {support_pairs}, all pairs would be {all_pairs})" + ) + + +def test_a_row_with_no_quadratic_terms_is_still_seen_as_linear(q_budget): + """``_quadratic_row_has_terms`` decides the linear/quadratic split for every QCP + row. It used ``np.any(np.abs(Q) > tol)``, which does not mean what it looks like + on a scipy sparse matrix — the split must not move when Q sparsifies.""" + n = 40 + + def linear(x): + return float(np.asarray(x, dtype=np.float64)[:3].sum()) + 2.0 + + for budget in (10**12, 1): + q_budget(budget) + Q, c, d = _extract_quadratic_coefficients_from_values(linear, n) + assert _quadratic_row_has_terms(Q) is False + assert d == pytest.approx(2.0) + assert np.allclose(np.asarray(c)[:3], 1.0) + assert np.allclose(np.asarray(c)[3:], 0.0) + + q_budget(1) + evaluate, _p, _Q, _c, _d = _quadratic_evaluator(n) + Q_q, _c_q, _d_q = _extract_quadratic_coefficients_from_values(evaluate, n) + assert sp.issparse(Q_q) + assert _quadratic_row_has_terms(Q_q) is True + + +def test_quadratic_row_has_terms_respects_its_tolerance_when_sparse(q_budget): + """A stored value below ``tol`` must not count, sparse or dense — otherwise the + sparse arm would classify a numerically-linear row as quadratic.""" + tiny = sp.csr_matrix(([1e-15, -1e-15], ([0, 1], [1, 0])), shape=(3, 3)) + assert _quadratic_row_has_terms(tiny) is False + real = sp.csr_matrix(([2.0], ([0], [0])), shape=(3, 3)) + assert _quadratic_row_has_terms(real) is True + assert _quadratic_row_has_terms(np.zeros((3, 3))) is False diff --git a/python/tests/test_cert_finite_bound_soundness.py b/python/tests/test_cert_finite_bound_soundness.py index 1bfd2555..7b12ef17 100644 --- a/python/tests/test_cert_finite_bound_soundness.py +++ b/python/tests/test_cert_finite_bound_soundness.py @@ -68,9 +68,17 @@ def test_optimal_with_matching_bound_is_preserved(self): assert r.gap_certified is True assert r.bound == pytest.approx(31.0) - def test_finite_zero_bound_is_a_valid_certificate(self): + def test_finite_zero_bound_survives_the_guard(self): # bound == 0.0 is finite (a weak but sound lower bound, e.g. ex1252): - # the guard must NOT mistake a falsy-but-finite value for "no bound". + # the guard must NOT mistake a falsy-but-finite value for "no bound" and + # clear it. That is what this test is for, and the surviving ``bound`` is + # the assertion that checks it. + # + # ``gap_certified`` is a separate question and is False here for a reason + # unrelated to the bound: there is no incumbent, so no gap was ever formed + # and none can be certified (#875 added that second-end check; this case + # used to report a certified gap with ``objective=None, gap=inf``). The + # companion case with an incumbent is below. r = SolveResult( status="time_limit", objective=None, @@ -78,8 +86,23 @@ def test_finite_zero_bound_is_a_valid_certificate(self): gap=float("inf"), gap_certified=True, ) + assert r.bound == pytest.approx(0.0), "a falsy-but-finite bound must not be cleared" + assert r.gap_certified is False + assert r.gap is None + + def test_finite_zero_bound_is_a_valid_certificate_with_an_incumbent(self): + # The same falsy-but-finite bound, now with both ends of a gap present: + # certification stands. + r = SolveResult( + status="time_limit", + objective=5.0, + bound=0.0, + gap=1.0, + gap_certified=True, + ) assert r.gap_certified is True assert r.bound == pytest.approx(0.0) + assert r.gap == pytest.approx(1.0) def test_infeasibility_certificate_is_exempt(self): # An infeasibility certificate certifies infeasibility, not a gap, and