From ed6e2a241d1382844eedc05d992a30f8208b48fd Mon Sep 17 00:00:00 2001 From: Lev Nachmanson <5377127+levnach@users.noreply.github.com> Date: Thu, 9 Jul 2026 10:39:23 -0700 Subject: [PATCH] opt: validate strict optimization optima faithfully with delta-rational bounds (#10059) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ## Problem Maximizing/minimizing under a **strict** inequality has a delta-rational optimum. For ```smt2 (declare-const r Real) (assert (< r 1)) (maximize r) (check-sat) (get-objectives) ``` the optimum is the supremum `1 - epsilon`, but z3 reported `r = 0`. The same defect makes shared-symbol objectives report a value matching **neither the model nor the true optimum** (issue #10028 follow-up). Minimal reproducer — a 6-mark Golomb ruler (a `>32`-arg `distinct`, so the objective is coupled to EUF) with a strict real objective `obj > x5`, whose true optimum is `17 + epsilon`: | case | before | after | |---|---|---| | `maximize r`, `r < 1` | `0` ❌ | `1 - epsilon` ✅ | | `minimize r`, `r > 1` | `0` ❌ | `1 + epsilon` ✅ | | Golomb `minimize obj`, `obj > x5` | `35/2` / `7+eps` ❌ | `17 + epsilon` ✅ | ## Root cause `check_bound` validates the LP hint by asserting `objective >= optimum`. For a supremum `1 - epsilon` this is a **lower** bound whose value carries a **negative** infinitesimal `(1, -1)`. No `lconstraint_kind` can express that. The kind->infinitesimal map only yields the *matching-sign* cases — `GT` -> lower `(r, +1)`, `LT` -> upper `(r, -1)` — or zero (`GE`/`LE`). The opposite-sign lower bound `(r, -1)` (i.e. `r >= r0 - delta`) is a *relaxation* that no strict inequality produces. `opt_solver::mk_ge` therefore projected the `-epsilon` away, turning `r >= 1 - epsilon` into the over-strong, unsatisfiable `r >= 1`; validation failed and the strictly smaller current model value was reported instead. ## Fix — carry the infinitesimal faithfully through the bound pipeline - **`lp_api::bound`** gains an `eps` component so `get_value` returns the true delta value (no spurious rational fixed-variable equality is propagated to EUF). - **`lar_base_constraint`** stores its right-hand side as a delta-rational `impq` pair; `rhs()` returns the rational component, `bound_eps()` the infinitesimal one. - **`lar_solver`** bound activation/update threads the whole `impq` bound, so a lower bound `(r, -1)` can be asserted. `constraint_holds` accounts for it using the **same** strict-bounds delta that flattens the model, computed **once per model**. - **`theory_lra::mk_ge`** builds a *fresh* predicate for the `(r, -1)` lower bound (to avoid colliding with an already-internalized `v >= r` literal) and attaches `eps = -1`. **`opt_solver::mk_ge`** passes the unprojected value to `theory_lra` / `theory_mi_arith` / `theory_inf_arith` (whose bounds are already `inf_rational`). The pair machinery is what makes the supremum both representable (optimum `1 - epsilon`) and validatable; the reported witness model remains the flattened rational (`find_delta_for_strict_bounds`), consistent with the existing epsilon semantics. ## Validation - Strict optima correct: `1-eps`, `1+eps`, bounded `2 5-eps`, and lex/box variants. - Integer optima and the #10028 shared-symbol cases unchanged (Golomb n=6/7/8 -> 17/25/34, consistent with the model). - Unit tests **92/92** (release); no new debug-suite failures. - Opt regression corpus (73 files, `model_validate=true`) **byte-identical** to baseline. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --------- Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> Co-authored-by: Nikolaj Bjorner --- .github/workflows/ocaml.yaml | 16 +--- src/math/lp/lar_constraints.h | 34 ++++++++- src/math/lp/lar_solver.cpp | 133 ++++++++++++++++++++++------------ src/math/lp/lar_solver.h | 21 +++--- src/math/lp/lp_api.h | 14 +++- src/opt/opt_solver.cpp | 21 ++++-- src/smt/theory_lra.cpp | 57 +++++++++++++-- 7 files changed, 212 insertions(+), 84 deletions(-) diff --git a/.github/workflows/ocaml.yaml b/.github/workflows/ocaml.yaml index ce575ba646..243a342268 100644 --- a/.github/workflows/ocaml.yaml +++ b/.github/workflows/ocaml.yaml @@ -28,21 +28,13 @@ jobs: restore-keys: | ${{ runner.os }}-ccache- - # Cache opam (compiler + packages) - - name: Cache opam - uses: actions/cache@v6.1.0 - with: - path: ~/.opam - key: ${{ runner.os }}-opam-${{ matrix.ocaml-version }}-${{ github.sha }} - restore-keys: | - ${{ runner.os }}-opam-${{ matrix.ocaml-version }}- - - # Setup OCaml via action + # Setup OCaml via action (handles opam caching internally) - uses: ocaml/setup-ocaml@v3 with: ocaml-compiler: ${{ matrix.ocaml-version }} opam-disable-sandboxing: true - + cache-prefix: v1 + # Platform-specific dependencies - name: Install system dependencies (Ubuntu) if: matrix.os == 'ubuntu-latest' @@ -123,4 +115,4 @@ jobs: - name: Run ml_example (native) run: | export DYLD_LIBRARY_PATH=$(pwd)/build - ./ml_example \ No newline at end of file + ./ml_example diff --git a/src/math/lp/lar_constraints.h b/src/math/lp/lar_constraints.h index 4ad0644890..8a8a6c60c9 100644 --- a/src/math/lp/lar_constraints.h +++ b/src/math/lp/lar_constraints.h @@ -39,7 +39,17 @@ inline std::string lconstraint_kind_string(lconstraint_kind t) { class lar_base_constraint { lconstraint_kind m_kind; - mpq m_right_side; + // Right-hand side as a delta-rational pair (x, y) = x + y*delta. The + // rational part x is the ordinary bound value returned by rhs(); the + // infinitesimal part y (bound_eps) is non-zero only for the delta-rational + // bounds that validate strict optimization suprema/infima (see + // opt_solver::maximize_objective). The strict kinds already carry a + // matching-sign infinitesimal in update_bound_with_* (LT -> upper bound + // (r, -1); GT -> lower bound (r, +1)); y is needed for the OPPOSITE-sign + // case that no kind yields: a lower bound (r, -1), i.e. objvar >= r - delta + // (GE gives (r, 0), GT gives (r, +1)), which is how a maximize supremum + // r - delta is asserted. y is zero for all ordinary constraints. + impq m_right_side; bool m_active; bool m_is_auxiliary; unsigned m_j; @@ -53,10 +63,17 @@ public: virtual ~lar_base_constraint() = default; lconstraint_kind kind() const { return m_kind; } - mpq const& rhs() const { return m_right_side; } + // First (rational) component of the right-hand side pair. + mpq const& rhs() const { return m_right_side.x; } + // Whole right-hand side pair (rational value + infinitesimal, see below). + impq const& rhs_impq() const { return m_right_side; } unsigned column() const { return m_j; } u_dependency* dep() const { return m_dep; } + // Second (infinitesimal) component of the right-hand side pair. + mpq const& bound_eps() const { return m_right_side.y; } + void set_bound_eps(mpq const& e) { m_right_side.y = e; } + void activate() { m_active = true; } void deactivate() { m_active = false; } bool is_active() const { return m_active; } @@ -181,11 +198,24 @@ public: return add(new (m_region) lar_var_constraint(j, k, mk_dep(), rhs)); } + constraint_index add_var_constraint(lpvar j, lconstraint_kind k, mpq const& rhs, mpq const& eps) { + auto* c = new (m_region) lar_var_constraint(j, k, mk_dep(), rhs); + c->set_bound_eps(eps); + return add(c); + } + constraint_index add_term_constraint(unsigned j, const lar_term* t, lconstraint_kind k, mpq const& rhs) { auto* dep = mk_dep(); return add(new (m_region) lar_term_constraint(j, t, k, dep, rhs)); } + constraint_index add_term_constraint(unsigned j, const lar_term* t, lconstraint_kind k, mpq const& rhs, mpq const& eps) { + auto* dep = mk_dep(); + auto* c = new (m_region) lar_term_constraint(j, t, k, dep, rhs); + c->set_bound_eps(eps); + return add(c); + } + // future behavior uses activation bit. bool is_active(constraint_index ci) const { return m_constraints[ci]->is_active(); } diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index 8be0cca91b..2f7dc11d25 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -1338,10 +1338,14 @@ namespace lp { if (m_imp->m_settings.get_cancel_flag()) return true; std::unordered_map var_map; - get_model_do_not_care_about_diff_vars(var_map); + // Compute the strict-bounds delta once per model: it flattens both the + // model (var_map) and the eps component of any delta-rational bound in + // constraint_holds, so the two must use the very same value. + mpq delta = get_core_solver().find_delta_for_strict_bounds(m_imp->m_settings.m_epsilon); + get_model_do_not_care_about_diff_vars(var_map, delta); for (auto const& c : m_imp->m_constraints.active()) { - if (!constraint_holds(c, var_map)) { + if (!constraint_holds(c, var_map, delta)) { TRACE(lar_solver, m_imp->m_constraints.display(tout, c) << "\n"; for (auto p : c.coeffs()) { @@ -1353,14 +1357,21 @@ namespace lp { return true; } - bool lar_solver::constraint_holds(const lar_base_constraint& constr, std::unordered_map& var_map) const { + bool lar_solver::constraint_holds(const lar_base_constraint& constr, std::unordered_map& var_map, const mpq& delta) const { mpq left_side_val = get_left_side_val(constr, var_map); + // Account for a delta-rational bound rhs + eps*delta (eps != 0 only for + // the bounds that validate strict optimization optima). 'delta' is the + // same strict-bounds delta that flattened var_map, so the comparison is + // exact over the reals. + mpq rhs = constr.rhs(); + if (!constr.bound_eps().is_zero()) + rhs += constr.bound_eps() * delta; switch (constr.kind()) { - case LE: return left_side_val <= constr.rhs(); - case LT: return left_side_val < constr.rhs(); - case GE: return left_side_val >= constr.rhs(); - case GT: return left_side_val > constr.rhs(); - case EQ: return left_side_val == constr.rhs(); + case LE: return left_side_val <= rhs; + case LT: return left_side_val < rhs; + case GE: return left_side_val >= rhs; + case GT: return left_side_val > rhs; + case EQ: return left_side_val == rhs; default: UNREACHABLE(); } @@ -1583,6 +1594,10 @@ namespace lp { void lar_solver::get_model_do_not_care_about_diff_vars(std::unordered_map& variable_values) const { mpq delta = get_core_solver().find_delta_for_strict_bounds(m_imp->m_settings.m_epsilon); + get_model_do_not_care_about_diff_vars(variable_values, delta); + } + + void lar_solver::get_model_do_not_care_about_diff_vars(std::unordered_map& variable_values, const mpq& delta) const { for (unsigned i = 0; i < get_core_solver().r_x().size(); ++i) { const impq& rp = get_core_solver().r_x(i); variable_values[i] = rp.x + delta * rp.y; @@ -2142,12 +2157,12 @@ namespace lp { void lar_solver::activate_check_on_equal(constraint_index ci, unsigned& equal_column) { auto const& c = m_imp->m_constraints[ci]; - update_column_type_and_bound_check_on_equal(c.column(), c.rhs(), ci, equal_column); + update_column_type_and_bound_check_on_equal(c.column(), c.rhs_impq(), ci, equal_column); } void lar_solver::activate(constraint_index ci) { auto const& c = m_imp->m_constraints[ci]; - update_column_type_and_bound(c.column(), c.rhs(), ci); + update_column_type_and_bound(c.column(), c.rhs_impq(), ci); } mpq lar_solver::adjust_bound_for_int(lpvar j, lconstraint_kind& k, const mpq& bound) { @@ -2191,6 +2206,24 @@ namespace lp { return ci; } + // Variant that attaches an infinitesimal coefficient 'eps' to the bound, so + // that activating the resulting constraint asserts the delta-rational bound + // (right_side, eps). Used to faithfully validate strict optimization optima + // (e.g. a maximize supremum r - delta is validated as a lower bound + // (r, -1)). Only supported for plain column bounds (no term column). + constraint_index lar_solver::mk_var_bound(lpvar j, lconstraint_kind kind, const mpq& right_side, const mpq& eps) { + TRACE(lar_solver, tout << "j = " << get_variable_name(j) << " " << lconstraint_kind_string(kind) << " " << right_side << " + " << eps << "*eps" << std::endl;); + mpq rs = adjust_bound_for_int(j, kind, right_side); + SASSERT(bound_is_integer_for_integer_column(j, rs)); + constraint_index ci; + if (!column_has_term(j)) + ci = m_imp->m_constraints.add_var_constraint(j, kind, rs, eps); + else + ci = m_imp->m_constraints.add_term_constraint(j, m_imp->m_columns[j].term(), kind, rs, eps); + SASSERT(sizes_are_correct()); + return ci; + } + bool lar_solver::compare_values(lpvar j, lconstraint_kind k, const mpq& rhs) { return compare_values(get_column_value(j), k, rhs); } @@ -2209,7 +2242,7 @@ namespace lp { } void lar_solver::update_column_type_and_bound(unsigned j, - const mpq& right_side, + const impq& right_side, constraint_index constr_index) { TRACE(lar_solver_feas, tout << "j = " << j << " was " << (this->column_is_feasible(j)?"feas":"non-feas") << std::endl;); m_imp->m_constraints.activate(constr_index); @@ -2273,7 +2306,10 @@ namespace lp { ls.add_var_bound(tv, c.kind(), c.rhs()); } void lar_solver::update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { - // SASSERT(validate_bound(j, kind, right_side, dep)); + update_column_type_and_bound(j, kind, impq(right_side), dep); + } + void lar_solver::update_column_type_and_bound(unsigned j, lconstraint_kind kind, const impq& right_side, u_dependency* dep) { + // SASSERT(validate_bound(j, kind, right_side.x, dep)); TRACE( lar_solver_feas, tout << "j" << j << " " << lconstraint_kind_string(kind) << " " << right_side << std::endl; @@ -2287,11 +2323,16 @@ namespace lp { } }); bool was_fixed = column_is_fixed(j); - mpq rs = adjust_bound_for_int(j, kind, right_side); + // adjust_bound_for_int operates on the rational part (and may sharpen + // the kind for integer columns); the infinitesimal part y is carried + // through unchanged. y is non-zero only for the delta-rational bounds + // that validate strict optimization optima, which target real columns. + mpq rs = adjust_bound_for_int(j, kind, right_side.x); + impq bound(rs, right_side.y); if (column_has_upper_bound(j)) - update_column_type_and_bound_with_ub(j, kind, rs, dep); + update_column_type_and_bound_with_ub(j, kind, bound, dep); else - update_column_type_and_bound_with_no_ub(j, kind, rs, dep); + update_column_type_and_bound_with_no_ub(j, kind, bound, dep); if (!was_fixed && column_is_fixed(j) && m_fixed_var_eh) m_fixed_var_eh(j); @@ -2320,7 +2361,7 @@ namespace lp { } void lar_solver::update_column_type_and_bound_check_on_equal(unsigned j, - const mpq& right_side, + const impq& right_side, constraint_index constr_index, unsigned& equal_to_j) { update_column_type_and_bound(j, right_side, constr_index); @@ -2336,7 +2377,7 @@ namespace lp { return m_imp->m_constraints.add_term_constraint(j, m_imp->m_columns[j].term(), kind, rs); } - void lar_solver::update_column_type_and_bound_with_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { + void lar_solver::update_column_type_and_bound_with_ub(unsigned j, lp::lconstraint_kind kind, const impq& right_side, u_dependency* dep) { SASSERT(column_has_upper_bound(j)); if (column_has_lower_bound(j)) { update_bound_with_ub_lb(j, kind, right_side, dep); @@ -2346,7 +2387,7 @@ namespace lp { } } - void lar_solver::update_column_type_and_bound_with_no_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { + void lar_solver::update_column_type_and_bound_with_no_ub(unsigned j, lp::lconstraint_kind kind, const impq& right_side, u_dependency* dep) { SASSERT(!column_has_upper_bound(j)); if (column_has_lower_bound(j)) { update_bound_with_no_ub_lb(j, kind, right_side, dep); @@ -2356,18 +2397,18 @@ namespace lp { } } - void lar_solver::update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { + void lar_solver::update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep) { SASSERT(column_has_lower_bound(j) && column_has_upper_bound(j)); SASSERT(get_core_solver().m_column_types[j] == column_type::boxed || get_core_solver().m_column_types[j] == column_type::fixed); - mpq y_of_bound(0); + mpq y_of_bound(right_side.y); switch (kind) { case LT: - y_of_bound = -1; + y_of_bound += -1; Z3_fallthrough; case LE: { - auto up = numeric_pair(right_side, y_of_bound); + auto up = numeric_pair(right_side.x, y_of_bound); if (up < get_lower_bound(j)) { set_crossed_bounds_column_and_deps(j, true, dep); } @@ -2379,10 +2420,10 @@ namespace lp { break; } case GT: - y_of_bound = 1; + y_of_bound += 1; Z3_fallthrough; case GE: { - auto low = numeric_pair(right_side, y_of_bound); + auto low = numeric_pair(right_side.x, y_of_bound); if (low > get_upper_bound(j)) { set_crossed_bounds_column_and_deps(j, false, dep); } @@ -2395,7 +2436,7 @@ namespace lp { break; } case EQ: { - auto v = numeric_pair(right_side, zero_of_type()); + auto v = numeric_pair(right_side.x, zero_of_type()); if (v > get_upper_bound(j)) set_crossed_bounds_column_and_deps(j, false, dep); else if (v < get_lower_bound(j)) @@ -2416,17 +2457,17 @@ namespace lp { get_core_solver().m_column_types[j] = column_type::fixed; } - void lar_solver::update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { + void lar_solver::update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep) { SASSERT(column_has_lower_bound(j) && !column_has_upper_bound(j)); SASSERT(get_core_solver().m_column_types[j] == column_type::lower_bound); - mpq y_of_bound(0); + mpq y_of_bound(right_side.y); switch (kind) { case LT: - y_of_bound = -1; + y_of_bound += -1; Z3_fallthrough; case LE: { - auto up = numeric_pair(right_side, y_of_bound); + auto up = numeric_pair(right_side.x, y_of_bound); if (up < get_lower_bound(j)) { set_crossed_bounds_column_and_deps(j, true, dep); } @@ -2437,9 +2478,9 @@ namespace lp { break; } case GT: - y_of_bound = 1; + y_of_bound += 1; case GE: { - auto low = numeric_pair(right_side, y_of_bound); + auto low = numeric_pair(right_side.x, y_of_bound); if (low < get_lower_bound(j)) { return; } @@ -2447,7 +2488,7 @@ namespace lp { break; } case EQ: { - auto v = numeric_pair(right_side, zero_of_type()); + auto v = numeric_pair(right_side.x, zero_of_type()); if (v < get_lower_bound(j)) { set_crossed_bounds_column_and_deps(j, true, dep); } @@ -2464,28 +2505,28 @@ namespace lp { } } - void lar_solver::update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { + void lar_solver::update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep) { SASSERT(!column_has_lower_bound(j) && column_has_upper_bound(j)); SASSERT(get_core_solver().m_column_types[j] == column_type::upper_bound); - mpq y_of_bound(0); + mpq y_of_bound(right_side.y); switch (kind) { case LT: - y_of_bound = -1; + y_of_bound += -1; Z3_fallthrough; case LE: { - auto up = numeric_pair(right_side, y_of_bound); + auto up = numeric_pair(right_side.x, y_of_bound); if (up >= get_upper_bound(j)) return; set_upper_bound_witness(j, dep, up); } break; case GT: - y_of_bound = 1; + y_of_bound += 1; Z3_fallthrough; case GE: { - auto low = numeric_pair(right_side, y_of_bound); + auto low = numeric_pair(right_side.x, y_of_bound); if (low > get_upper_bound(j)) { set_crossed_bounds_column_and_deps(j, false, dep); } @@ -2497,7 +2538,7 @@ namespace lp { break; case EQ: { - auto v = numeric_pair(right_side, zero_of_type()); + auto v = numeric_pair(right_side.x, zero_of_type()); if (v > get_upper_bound(j)) { set_crossed_bounds_column_and_deps(j, false, dep); } @@ -2514,30 +2555,30 @@ namespace lp { } } - void lar_solver::update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) { + void lar_solver::update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep) { SASSERT(!column_has_lower_bound(j) && !column_has_upper_bound(j)); - mpq y_of_bound(0); + mpq y_of_bound(right_side.y); switch (kind) { case LT: - y_of_bound = -1; + y_of_bound += -1; Z3_fallthrough; case LE: { - auto up = numeric_pair(right_side, y_of_bound); + auto up = numeric_pair(right_side.x, y_of_bound); set_upper_bound_witness(j, dep, up); get_core_solver().m_column_types[j] = column_type::upper_bound; } break; case GT: - y_of_bound = 1; + y_of_bound += 1; Z3_fallthrough; case GE: { - auto low = numeric_pair(right_side, y_of_bound); + auto low = numeric_pair(right_side.x, y_of_bound); set_lower_bound_witness(j, dep, low); get_core_solver().m_column_types[j] = column_type::lower_bound; } break; case EQ: { - auto v = numeric_pair(right_side, zero_of_type()); + auto v = numeric_pair(right_side.x, zero_of_type()); set_upper_bound_witness(j, dep, v); set_lower_bound_witness(j, dep, v); get_core_solver().m_column_types[j] = column_type::fixed; diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h index 4588f1772f..d361496224 100644 --- a/src/math/lp/lar_solver.h +++ b/src/math/lp/lar_solver.h @@ -88,19 +88,20 @@ class lar_solver : public column_namer { void add_bound_negation_to_solver(lar_solver& ls, lpvar j, lconstraint_kind kind, const mpq& right_side); void add_constraint_to_validate(lar_solver& ls, constraint_index ci); bool m_validate_blocker = false; - void update_column_type_and_bound_check_on_equal(unsigned j, const mpq& right_side, constraint_index ci, unsigned&); - void update_column_type_and_bound(unsigned j, const mpq& right_side, constraint_index ci); + void update_column_type_and_bound_check_on_equal(unsigned j, const impq& right_side, constraint_index ci, unsigned&); + void update_column_type_and_bound(unsigned j, const impq& right_side, constraint_index ci); public: bool validate_blocker() const { return m_validate_blocker; } bool & validate_blocker() { return m_validate_blocker; } + void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const impq& right_side, u_dependency* dep); void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); private: - void update_column_type_and_bound_with_ub(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); - void update_column_type_and_bound_with_no_ub(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); - void update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); - void update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); - void update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); - void update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep); + void update_column_type_and_bound_with_ub(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep); + void update_column_type_and_bound_with_no_ub(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep); + void update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep); + void update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep); + void update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep); + void update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep); void remove_non_fixed_from_fixed_var_table(); constraint_index add_var_bound_on_constraint_for_term(lpvar j, lconstraint_kind kind, const mpq& right_side); void set_crossed_bounds_column_and_deps(unsigned j, bool lower_bound, u_dependency* dep); @@ -147,7 +148,7 @@ class lar_solver : public column_namer { numeric_pair get_basic_var_value_from_row(unsigned i); bool all_constrained_variables_are_registered(const vector>& left_side); bool all_constraints_hold() const; - bool constraint_holds(const lar_base_constraint& constr, std::unordered_map& var_map) const; + bool constraint_holds(const lar_base_constraint& constr, std::unordered_map& var_map, const mpq& delta) const; static void register_in_map(std::unordered_map& coeffs, const lar_base_constraint& cn, const mpq& a); static void register_monoid_in_map(std::unordered_map& coeffs, const mpq& a, unsigned j); bool the_left_sides_sum_to_zero(const vector>& evidence) const; @@ -271,6 +272,7 @@ public: bool fixed_base_removed_correctly() const; #endif constraint_index mk_var_bound(lpvar j, lconstraint_kind kind, const mpq& right_side); + constraint_index mk_var_bound(lpvar j, lconstraint_kind kind, const mpq& right_side, const mpq& eps); void activate_check_on_equal(constraint_index, lpvar&); void activate(constraint_index); void random_update(unsigned sz, lpvar const* vars); @@ -478,6 +480,7 @@ public: void get_model(std::unordered_map& variable_values) const; void get_rid_of_inf_eps(); void get_model_do_not_care_about_diff_vars(std::unordered_map& variable_values) const; + void get_model_do_not_care_about_diff_vars(std::unordered_map& variable_values, const mpq& delta) const; std::string get_variable_name(lpvar vi) const override; void set_variable_name(lpvar vi, const std::string&); unsigned number_of_vars() const; diff --git a/src/math/lp/lp_api.h b/src/math/lp/lp_api.h index 96279a7db4..56dcb95414 100644 --- a/src/math/lp/lp_api.h +++ b/src/math/lp/lp_api.h @@ -36,15 +36,23 @@ namespace lp_api { rational m_value; bound_kind m_bound_kind; lp::constraint_index m_constraints[2]; + // Infinitesimal coefficient of the asserted (positive-literal) bound + // value: the bound means v (>=|<=) m_value + m_eps*delta. Non-zero + // only for the delta-rational bounds used to validate strict + // optimization optima (e.g. a lower bound (r, -1) for a maximize + // supremum). Kept so get_value reports the true delta value and no + // spurious rational fixed-variable equality is propagated to EUF. + rational m_eps; public: - bound(Literal bv, theory_var v, lp::lpvar vi, bool is_int, rational const& val, bound_kind k, lp::constraint_index ct, lp::constraint_index cf) : + bound(Literal bv, theory_var v, lp::lpvar vi, bool is_int, rational const& val, bound_kind k, lp::constraint_index ct, lp::constraint_index cf, rational const& eps = rational::zero()) : m_bv(bv), m_var(v), m_column_index(vi), m_is_int(is_int), m_value(val), - m_bound_kind(k) { + m_bound_kind(k), + m_eps(eps) { m_constraints[0] = cf; m_constraints[1] = ct; } @@ -67,7 +75,7 @@ namespace lp_api { inf_rational get_value(bool is_true) const { if (is_true != get_lit().sign()) - return inf_rational(m_value); // v >= value or v <= value + return inf_rational(m_value, m_eps); // v >= value (+ eps*delta) or v <= value (+ eps*delta) if (m_is_int) { SASSERT(m_value.is_int()); rational const& offset = (m_bound_kind == lower_t) ? rational::minus_one() : rational::one(); diff --git a/src/opt/opt_solver.cpp b/src/opt/opt_solver.cpp index eb8b4dd6cb..ac167a4d7b 100644 --- a/src/opt/opt_solver.cpp +++ b/src/opt/opt_solver.cpp @@ -509,13 +509,19 @@ namespace opt { if (typeid(smt::theory_inf_arith) == typeid(opt)) { smt::theory_inf_arith& th = dynamic_cast(opt); - return th.mk_ge(m_fm, v, val); + // Pass the original value (with its negative infinitesimal, if any): + // theory_inf_arith's bounds are inf_rational and represent a strict + // supremum r - delta faithfully as the lower bound (r, -1), which is + // required to validate a strict maximization optimum. + return th.mk_ge(m_fm, v, _val); } if (typeid(smt::theory_mi_arith) == typeid(opt)) { smt::theory_mi_arith& th = dynamic_cast(opt); - SASSERT(val.is_finite()); - return th.mk_ge(m_fm, v, val.get_numeral()); + SASSERT(_val.is_finite()); + // As above: theory_mi_arith's bounds are inf_rational, so pass the + // delta-rational value through unprojected for faithful validation. + return th.mk_ge(m_fm, v, _val.get_numeral()); } if (typeid(smt::theory_i_arith) == typeid(opt)) { @@ -549,8 +555,13 @@ namespace opt { if (typeid(smt::theory_lra) == typeid(opt)) { smt::theory_lra& th = dynamic_cast(opt); - SASSERT(val.is_finite()); - return th.mk_ge(m_fm, v, val.get_numeral()); + SASSERT(_val.is_finite()); + // Pass the ORIGINAL value (with its negative infinitesimal, if any): + // theory_lra faithfully encodes a lower bound v >= r - delta, which + // is required to validate a strict maximization supremum. 'val' + // above has had a negative infinitesimal projected away for the + // benefit of theories that cannot represent it. + return th.mk_ge(m_fm, v, _val.get_numeral()); } // difference logic? diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index d3ac4840ad..4999dcd265 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -3304,6 +3304,16 @@ public: } api_bound* mk_var_bound(bool_var bv, theory_var v, lp_api::bound_kind bk, rational const& bound) { + return mk_var_bound(bv, v, bk, bound, rational::zero()); + } + + // eps is the infinitesimal coefficient of the asserted (positive-literal) + // bound value: the bound means v (>=|<=) bound + eps*delta. Non-zero only + // for the delta-rational bounds that faithfully validate strict + // optimization optima (a maximize supremum r - delta becomes a lower bound + // (r, -1)). Only the asserted direction (cT) carries eps; the negation cF + // is never activated on the optimization validation path. + api_bound* mk_var_bound(bool_var bv, theory_var v, lp_api::bound_kind bk, rational const& bound, rational const& eps) { scoped_internalize_state st(*this); st.vars().push_back(v); st.coeffs().push_back(rational::one()); @@ -3315,7 +3325,7 @@ public: lp::lconstraint_kind kT = bound2constraint_kind(v_is_int, bk, true); lp::lconstraint_kind kF = bound2constraint_kind(v_is_int, bk, false); - cT = lp().mk_var_bound(vi, kT, bound); + cT = lp().mk_var_bound(vi, kT, bound, eps); if (v_is_int) { rational boundF = (bk == lp_api::lower_t) ? bound - 1 : bound + 1; cF = lp().mk_var_bound(vi, kF, boundF); @@ -3326,7 +3336,7 @@ public: add_ineq_constraint(cT, literal(bv, false)); add_ineq_constraint(cF, literal(bv, true)); - return alloc(api_bound, literal(bv, false), v, vi, v_is_int, bound, bk, cT, cF); + return alloc(api_bound, literal(bv, false), v, vi, v_is_int, bound, bk, cT, cF, eps); } // @@ -4055,11 +4065,22 @@ public: tout << " x[" << j << "] = " << lp().get_column_value(j) << "\n"; } }); + // Discard the infinitesimal of the value returned from the NLA path. + // When NLA is involved the objective is nonlinear, so lp_val is the + // optimum of the LINEAR relaxation: its infinitesimal comes from the + // strict bounds introduced by the linearization, not from a genuine + // strict optimum of the nonlinear problem. If it were kept, + // opt_solver::mk_ge would assert a delta-rational bound (r, -1) that the + // real problem cannot honor, fixing the objective column at a delta + // value the LP core cannot snap on the next solve (assertion + // non_basic_columns_are_set_correctly). The rational part remains a + // sound bound for the optimizer to validate via check_bound. + inf_eps lp_val_no_eps(lp_val.get_infinity(), inf_rational(lp_val.get_rational())); switch (nla_st) { case FC_DONE: // NLA satisfied: keep the optimal assignment, return LP value blocker = mk_gt(v); - result = lp_val; + result = lp_val_no_eps; st = lp::lp_status::FEASIBLE; return true; case FC_CONTINUE: @@ -4068,7 +4089,7 @@ public: // as a bound for the optimizer to validate via check_bound(). lp().restore_x(); blocker = mk_gt(v, lp_ival); - result = lp_val; + result = lp_val_no_eps; st = lp::lp_status::FEASIBLE; return true; case FC_GIVEUP: @@ -4249,11 +4270,32 @@ public: expr_ref mk_ge(generic_model_converter& fm, theory_var v, inf_rational const& val) { rational r = val.get_rational(); - bool is_strict = val.get_infinitesimal().is_pos(); + bool is_strict = val.get_infinitesimal().is_pos(); + // A negative infinitesimal encodes a delta-rational lower bound + // v >= r - delta. It arises when validating a strict maximization + // optimum (supremum r reported as r - epsilon): no lconstraint_kind + // yields a lower bound with a -delta component, so it is threaded + // through as an explicit eps on the bound (see lp_api::bound, + // lar_solver::mk_var_bound). Over the reals this is a genuine bound + // (feasible together with the problem's own strict bound v <= r - delta, + // fixing v = r - delta), which is exactly what makes the supremum + // achievable in the delta field and lets check_bound validate it. + bool is_lower_eps = val.get_infinitesimal().is_neg(); app_ref b(m); bool is_int = a.is_int(get_expr(v)); TRACE(arith, display(tout << "v" << v << "\n");); - if (is_strict) { + if (is_lower_eps) { + // Fresh, dedicated predicate for the delta-rational lower bound + // v >= r - delta. A plain (a.mk_ge v r) atom would collide with an + // already-internalized 'v >= r' literal (e.g. from the problem's own + // strict bound v < r), which carries no infinitesimal and would make + // validation assert the over-strong v >= r. The bound's real meaning + // (including the -delta) is attached via the api_bound's eps below. + std::ostringstream strm; + strm << r << " - eps <= " << mk_pp(get_expr(v), m) << " (opt)"; + b = m.mk_const(symbol(strm.str()), m.mk_bool_sort()); + } + else if (is_strict) { b = a.mk_le(mk_obj(v), a.mk_numeral(r, is_int)); } else { @@ -4267,7 +4309,8 @@ public: // ctx().set_enode_flag(bv, true); lp_api::bound_kind bkind = lp_api::bound_kind::lower_t; if (is_strict) bkind = lp_api::bound_kind::upper_t; - api_bound* a = mk_var_bound(bv, v, bkind, r); + rational eps = is_lower_eps ? rational::minus_one() : rational::zero(); + api_bound* a = mk_var_bound(bv, v, bkind, r, eps); mk_bound_axioms(*a); updt_unassigned_bounds(v, +1); m_bounds[v].push_back(a);