From f508854fe5c6f890928c397ecbe44240a15ab3f8 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson <5377127+levnach@users.noreply.github.com> Date: Sun, 14 Jun 2026 16:25:21 -0700 Subject: [PATCH] Lcube (#9858) Implemented the largest cube heuristic from Bromberger and Weidenbach's paper on cubes. Also fixes an overflow bug in mzp. Use vswhere to find the visual studio version on windows in the build's ymls. --------- Signed-off-by: Lev Nachmanson Co-authored-by: Claude Fable 5 Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com> --- .github/workflows/Windows.yml | 10 +- .github/workflows/nightly.yml | 9 +- .github/workflows/nuget-build.yml | 9 +- .github/workflows/release.yml | 9 +- src/math/lp/int_cube.cpp | 249 ++++++++++++++++++++++++++++ src/math/lp/int_cube.h | 29 +++- src/math/lp/int_solver.cpp | 26 ++- src/math/lp/lar_solver.cpp | 7 +- src/math/lp/lar_solver.h | 4 +- src/math/lp/lp_params_helper.pyg | 4 +- src/math/lp/lp_settings.cpp | 2 + src/math/lp/lp_settings.h | 10 ++ src/smt/theory_lra.cpp | 2 +- src/test/CMakeLists.txt | 1 + src/test/lcube.cpp | 261 ++++++++++++++++++++++++++++++ src/test/lp/lp.cpp | 2 +- src/test/main.cpp | 3 +- src/util/mpz.cpp | 7 +- 18 files changed, 620 insertions(+), 24 deletions(-) create mode 100644 src/test/lcube.cpp diff --git a/.github/workflows/Windows.yml b/.github/workflows/Windows.yml index afd36f7d1..e35f79fa6 100644 --- a/.github/workflows/Windows.yml +++ b/.github/workflows/Windows.yml @@ -37,9 +37,10 @@ jobs: ${{ matrix.cmd1 }} ${{ matrix.cmd2 }} ${{ matrix.cmd3 }} - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" ${{ matrix.arch }} - cmake ${{ matrix.bindings }} -G "NMake Makefiles" ../ - nmake + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" ${{ matrix.arch }} || exit /b 1 + cmake ${{ matrix.bindings }} -G "NMake Makefiles" ../ || exit /b 1 + nmake || exit /b 1 cd .. shell: cmd - name: Run Regressions @@ -52,7 +53,8 @@ jobs: if: ${{ matrix.test }} run: | pushd build - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" ${{ matrix.arch }} + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" ${{ matrix.arch }} || exit /b 1 pushd build\python python z3test.py z3 python z3test.py z3num diff --git a/.github/workflows/nightly.yml b/.github/workflows/nightly.yml index 26a512715..b6c8126a7 100644 --- a/.github/workflows/nightly.yml +++ b/.github/workflows/nightly.yml @@ -552,7 +552,8 @@ jobs: - name: Build shell: cmd run: | - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x64 + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x64 || exit /b 1 python scripts\mk_win_dist.py --x64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --zip - name: Upload artifact @@ -578,7 +579,8 @@ jobs: - name: Build shell: cmd run: | - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x86 + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x86 || exit /b 1 python scripts\mk_win_dist.py --x86-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --zip - name: Upload artifact @@ -604,7 +606,8 @@ jobs: - name: Build shell: cmd run: | - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64 + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64 || exit /b 1 python scripts\mk_win_dist_cmake.py --arm64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --assembly-version=${{ env.MAJOR }}.${{ env.MINOR }}.${{ env.PATCH }} --zip - name: Upload artifact diff --git a/.github/workflows/nuget-build.yml b/.github/workflows/nuget-build.yml index 6cf089243..ca890aeba 100644 --- a/.github/workflows/nuget-build.yml +++ b/.github/workflows/nuget-build.yml @@ -30,7 +30,8 @@ jobs: - name: Build Windows x64 shell: cmd run: | - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x64 + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x64 || exit /b 1 python scripts\mk_win_dist.py --x64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --assembly-version=${{ github.event.inputs.version || '4.17.0' }} --zip - name: Upload Windows x64 artifact @@ -54,7 +55,8 @@ jobs: - name: Build Windows x86 shell: cmd run: | - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x86 + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x86 || exit /b 1 python scripts\mk_win_dist.py --x86-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --assembly-version=${{ github.event.inputs.version || '4.17.0' }} --zip - name: Upload Windows x86 artifact @@ -78,7 +80,8 @@ jobs: - name: Build Windows ARM64 shell: cmd run: | - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64 + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64 || exit /b 1 python scripts\mk_win_dist_cmake.py --arm64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --assembly-version=${{ github.event.inputs.version || '4.17.0' }} --zip - name: Upload Windows ARM64 artifact diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 112b47960..f5911d28d 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -562,7 +562,8 @@ jobs: - name: Build shell: cmd run: | - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x64 + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x64 || exit /b 1 python scripts\mk_win_dist.py --x64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --zip - name: Upload artifact @@ -588,7 +589,8 @@ jobs: - name: Build shell: cmd run: | - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" x86 + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" x86 || exit /b 1 python scripts\mk_win_dist.py --x86-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --zip - name: Upload artifact @@ -614,7 +616,8 @@ jobs: - name: Build shell: cmd run: | - call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64 + for /f "usebackq delims=" %%i in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -latest -prerelease -products * -requires Microsoft.VisualStudio.Component.VC.Tools.x86.x64 -property installationPath`) do set "VSPATH=%%i" + call "%VSPATH%\VC\Auxiliary\Build\vcvarsall.bat" amd64_arm64 || exit /b 1 python scripts\mk_win_dist_cmake.py --arm64-only --dotnet-key=%GITHUB_WORKSPACE%\resources\z3.snk --assembly-version=${{ env.RELEASE_VERSION }} --zip - name: Upload artifact diff --git a/src/math/lp/int_cube.cpp b/src/math/lp/int_cube.cpp index aad7fae6e..37e046239 100644 --- a/src/math/lp/int_cube.cpp +++ b/src/math/lp/int_cube.cpp @@ -16,6 +16,9 @@ Author: Revision History: --*/ +#include +#include + #include "math/lp/int_solver.h" #include "math/lp/lar_solver.h" #include "math/lp/int_cube.h" @@ -81,6 +84,252 @@ namespace lp { SASSERT(lp_status::OPTIMAL == lra.get_status() || lp_status::FEASIBLE == lra.get_status()); } + // The largest cube test of Bromberger and Weidenbach: + // maximize x_e subject to Ax + a'(x_e/2) <= b, x_e >= 0, where a'_i = ||a_i||_1, + // with the 1-norm taken over the integer variables of the row. + // The solution is the center z of a largest cube contained in the polyhedron. + // If the maximal edge length is at least 1, then the rounding of z is + // an integer solution; otherwise the rounding is checked, and possibly repaired, + // against the original constraints. + lia_move int_cube::find_largest_cube() { + lia.settings().stats().m_lcube_calls++; + TRACE(cube, + for (unsigned j = 0; j < lra.number_of_vars(); ++j) + lia.display_column(tout, j); + tout << lra.constraints(); + ); + + lra.push(); + // The edge rows are ephemeral: suppress the add-term callback, + // dioph_eq's reaction to it is not undone by pop(). + auto add_term_cb = lra.m_add_term_callback; + lra.m_add_term_callback = nullptr; + unsigned x_e = lra.add_var(UINT_MAX, false); // the edge length of the cube + lra.add_var_bound(x_e, lconstraint_kind::GE, mpq(0)); + bool ok = add_cube_edge_rows(x_e); + lra.m_add_term_callback = add_term_cb; + if (!ok) { + lra.pop(); + lra.set_status(lp_status::OPTIMAL); + return lia_move::undef; + } + + lp_status st = lra.find_feasible_solution(); + if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { + TRACE(cube, tout << "cannot find a feasible solution";); + lra.pop(); + lra.move_non_basic_columns_to_bounds(); + // it can happen that we found an integer solution here + return !lra.r_basis_has_inf_int()? lia_move::sat: lia_move::undef; + } + + impq e; // the maximal edge length + st = lra.maximize_term(x_e, e, /*fix_int_cols*/ false); + if (lia.settings().get_cancel_flag()) { + lra.pop(); + return lia_move::undef; + } + if (st == lp_status::UNBOUNDED) { + // infinite lattice width: the polyhedron contains cubes of arbitrary edge length + lra.add_var_bound(x_e, lconstraint_kind::GE, mpq(1)); + st = lra.find_feasible_solution(); + if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { + lra.pop(); + return lia_move::undef; + } + lra.pop(); + return sat_after_rounding(); + } + TRACE(cube, tout << "max edge length = " << e << "\n";); + if (e >= impq(mpq(1))) { + lra.pop(); + return sat_after_rounding(); + } + // the largest cube is smaller than the unit cube: + // the rounded center is only a candidate + lra.pop(); + return round_and_repair(); + } + + bool int_cube::add_cube_edge_rows(unsigned x_e) { + // snapshot the term columns: add_edge_rows_for_term appends to lra.terms() + svector term_columns; + for (const lar_term* t : lra.terms()) + term_columns.push_back(t->j()); + for (unsigned j : term_columns) + if (!add_edge_rows_for_term(j, x_e)) { + TRACE(cube, tout << "cannot add the edge rows";); + return false; + } + return true; + } + + // i is the column index having the term + bool int_cube::add_edge_rows_for_term(unsigned i, unsigned x_e) { + if (!lra.column_associated_with_row(i)) + return true; + const lar_term& t = lra.get_term(i); + impq delta = get_cube_delta_for_term(t); + TRACE(cube, lra.print_term_as_indices(t, tout); tout << ", delta = " << delta << "\n";); + if (is_zero(delta)) + return true; + if (!is_zero(delta.y)) + // the infinitesimal delta does not scale with x_e: tighten statically, + // it is sound for any edge length + return lra.tighten_term_bounds_by_delta(i, delta); + if (lra.column_has_upper_bound(i)) { + impq u = lra.get_upper_bound(i); // copy: add_term invalidates bound references + vector> coeffs; + coeffs.push_back(std::make_pair(mpq(1), i)); + coeffs.push_back(std::make_pair(delta.x, x_e)); + unsigned s = lra.add_term(coeffs, UINT_MAX); + lra.add_var_bound(s, is_zero(u.y) ? lconstraint_kind::LE : lconstraint_kind::LT, u.x); + } + if (lra.column_has_lower_bound(i)) { + impq l = lra.get_lower_bound(i); // copy: add_term invalidates bound references + vector> coeffs; + coeffs.push_back(std::make_pair(mpq(1), i)); + coeffs.push_back(std::make_pair(-delta.x, x_e)); + unsigned s = lra.add_term(coeffs, UINT_MAX); + lra.add_var_bound(s, is_zero(l.y) ? lconstraint_kind::GE : lconstraint_kind::GT, l.x); + } + return true; + } + + lia_move int_cube::sat_after_rounding() { + lra.round_to_integer_solution(); + lra.set_status(lp_status::FEASIBLE); + SASSERT(lia.settings().get_cancel_flag() || lia.is_feasible()); + TRACE(cube, tout << "largest cube success";); + lia.settings().stats().m_lcube_success++; + return lia_move::sat; + } + + lia_move int_cube::round_and_repair() { + lra.backup_x(); // remember the cube center + vector flips; + for (unsigned j = 0; j < lra.column_count(); ++j) { + if (!lra.column_is_int(j) || lra.column_has_term(j)) + continue; + const impq& v = lra.get_column_value(j); + if (v.is_int()) + continue; + flip_candidate f; + f.m_j = j; + f.m_lo = floor(v); + flips.push_back(f); + } + lra.round_to_integer_solution(); + for (auto& f : flips) + f.m_at_hi = lra.get_column_value(f.m_j).x > f.m_lo; + if (repair_rounded_candidate(flips)) { + lra.set_status(lp_status::FEASIBLE); + SASSERT(lia.settings().get_cancel_flag() || lia.is_feasible()); + TRACE(cube, tout << "largest cube success";); + lia.settings().stats().m_lcube_success++; + return lia_move::sat; + } + // return to the cube center: an interior point of the polyhedron + lra.restore_x(); + lra.set_status(lp_status::FEASIBLE); + return lia_move::undef; + } + + // Checks the rounded center against the original constraints. On failure + // searches the vertices of the lattice cell around the center greedily: + // flip a coordinate between floor and ceiling to maximally decrease the + // total bound violation, within a budget. + bool int_cube::repair_rounded_candidate(vector& flips) { + vector rows; + for (const lar_term* t : lra.terms()) { + unsigned j = t->j(); + if (!lra.column_associated_with_row(j)) + continue; + if (!lra.column_has_upper_bound(j) && !lra.column_has_lower_bound(j)) + continue; + bounded_row r; + r.m_j = j; + r.m_val = t->apply(lra.r_x()); + rows.push_back(r); + } + auto row_violation = [&](unsigned ri, const impq& v) { + impq w; + unsigned j = rows[ri].m_j; + if (lra.column_has_upper_bound(j) && v > lra.get_upper_bound(j)) + w += v - lra.get_upper_bound(j); + if (lra.column_has_lower_bound(j) && v < lra.get_lower_bound(j)) + w += lra.get_lower_bound(j) - v; + return w; + }; + impq violation; + for (unsigned ri = 0; ri < rows.size(); ++ri) + violation += row_violation(ri, rows[ri].m_val); + if (is_zero(violation)) + return true; // the rounded center fits as it is + if (flips.empty()) + return false; + + std::unordered_map flip_of_var; + for (unsigned fi = 0; fi < flips.size(); ++fi) + flip_of_var[flips[fi].m_j] = fi; + // occurrences of the flip candidates in the bounded rows + vector>> occs; + occs.resize(flips.size()); + for (unsigned ri = 0; ri < rows.size(); ++ri) { + const lar_term& t = lra.get_term(rows[ri].m_j); + for (lar_term::ival p : t) { + auto it = flip_of_var.find(p.j()); + if (it != flip_of_var.end()) + occs[it->second].push_back(std::make_pair(ri, p.coeff())); + } + } + + unsigned budget = std::min(2 * flips.size(), lia.settings().lcube_flips()); + bool flipped = false; + while (!is_zero(violation) && budget-- > 0) { + unsigned best_fi = UINT_MAX; + impq best_gain; + for (unsigned fi = 0; fi < flips.size(); ++fi) { + if (occs[fi].empty()) + continue; + mpq step = flips[fi].m_at_hi ? mpq(-1) : mpq(1); + impq gain; + for (const auto& o : occs[fi]) { + const impq& v = rows[o.first].m_val; + gain += row_violation(o.first, v + impq(step * o.second)) - row_violation(o.first, v); + } + if (gain < best_gain) { + best_gain = gain; + best_fi = fi; + } + } + if (best_fi == UINT_MAX) + return false; // no flip decreases the violation + mpq step = flips[best_fi].m_at_hi ? mpq(-1) : mpq(1); + for (const auto& o : occs[best_fi]) + rows[o.first].m_val += impq(step * o.second); + flips[best_fi].m_at_hi = !flips[best_fi].m_at_hi; + violation += best_gain; + flipped = true; + TRACE(cube, tout << "flipped column " << flips[best_fi].m_j << ", violation = " << violation << "\n";); + } + if (!is_zero(violation)) + return false; + + // apply the repaired candidate + for (const auto& f : flips) + lra.set_column_value(f.m_j, impq(f.m_at_hi ? f.m_lo + 1 : f.m_lo)); + for (const lar_term* t : lra.terms()) { + unsigned j = t->j(); + if (!lra.column_associated_with_row(j)) + continue; + lra.set_column_value(j, t->apply(lra.r_x())); + } + if (flipped) + lia.settings().stats().m_lcube_flip_success++; + return true; + } + impq int_cube::get_cube_delta_for_term(const lar_term& t) const { if (t.size() == 2) { bool seen_minus = false; diff --git a/src/math/lp/int_cube.h b/src/math/lp/int_cube.h index 4addc096b..a7d672a52 100644 --- a/src/math/lp/int_cube.h +++ b/src/math/lp/int_cube.h @@ -10,9 +10,15 @@ Abstract: Cube finder This routine attempts to find a feasible integer solution - by tightnening bounds and running an LRA solver on the + by tightnening bounds and running an LRA solver on the tighter system. + find_largest_cube() implements the largest cube test of + Bromberger and Weidenbach (Fast Cube Tests for LIA Constraint + Solving, IJCAR 2016): a fresh variable x_e for the cube edge + length is introduced and maximized; the center of the largest + cube is rounded to a candidate integer solution. + Author: Nikolaj Bjorner (nbjorner) Lev Nachmanson (levnach) @@ -21,7 +27,10 @@ Revision History: --*/ #pragma once +#include "util/vector.h" #include "math/lp/lia_move.h" +#include "math/lp/numeric_pair.h" +#include "math/lp/lar_term.h" namespace lp { class int_solver; @@ -29,12 +38,30 @@ namespace lp { class int_cube { class int_solver& lia; class lar_solver& lra; + // a fractional integer coordinate of the cube center: + // the candidate value is m_lo or m_lo + 1 + struct flip_candidate { + unsigned m_j = 0; + mpq m_lo; + bool m_at_hi = false; + }; + // a term column with at least one bound, tracked during the repair + struct bounded_row { + unsigned m_j = 0; + impq m_val; + }; bool tighten_term_for_cube(unsigned i); bool tighten_terms_for_cube(); void find_feasible_solution(); impq get_cube_delta_for_term(const lar_term& t) const; + bool add_edge_rows_for_term(unsigned i, unsigned x_e); + bool add_cube_edge_rows(unsigned x_e); + lia_move sat_after_rounding(); + lia_move round_and_repair(); + bool repair_rounded_candidate(vector& flips); public: int_cube(int_solver& lia); lia_move operator()(); + lia_move find_largest_cube(); }; } diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp index cf0878597..aec275b77 100644 --- a/src/math/lp/int_solver.cpp +++ b/src/math/lp/int_solver.cpp @@ -44,7 +44,8 @@ namespace lp { dioph_eq m_dio; int_gcd_test m_gcd; unsigned m_initial_dio_calls_period; - + unsigned m_lcube_period; + bool column_is_int_inf(unsigned j) const { return lra.column_is_int(j) && (!lia.value_is_int(j)); } @@ -52,7 +53,8 @@ namespace lp { imp(int_solver& lia): lia(lia), lra(lia.lra), lrac(lia.lrac), m_hnf_cutter(lia), m_dio(lia), m_gcd(lia) { m_hnf_cut_period = settings().hnf_cut_period(); m_initial_dio_calls_period = settings().dio_calls_period(); - } + m_lcube_period = settings().m_int_find_cube_period; + } bool has_lower(unsigned j) const { switch (lrac.m_column_types()[j]) { @@ -196,6 +198,25 @@ namespace lp { return m_number_of_calls % settings().m_int_find_cube_period == 0; } + // The largest cube test is throttled exponentially: when the polyhedron + // does not contain a large enough cube it is unlikely to contain one + // later, after more constraints are added, so each failure doubles the + // period and a success resets it. + bool should_find_lcube() { + return settings().lcube() && m_number_of_calls % m_lcube_period == 0; + } + + lia_move find_lcube() { + lia_move r = int_cube(lia).find_largest_cube(); + if (r == lia_move::undef) { + if (m_lcube_period < (1u << 30)) + m_lcube_period *= 2; + } + else + m_lcube_period = settings().m_int_find_cube_period; + return r; + } + bool should_gomory_cut() { bool dio_allows_gomory = !settings().dio() || settings().dio_enable_gomory_cuts() || m_dio.some_terms_are_ignored(); @@ -246,6 +267,7 @@ namespace lp { ++m_number_of_calls; if (r == lia_move::undef) r = patch_basic_columns(); if (r == lia_move::undef && should_find_cube()) r = int_cube(lia)(); + if (r == lia_move::undef && should_find_lcube()) r = find_lcube(); if (r == lia_move::undef) lra.move_non_basic_columns_to_bounds(); if (r == lia_move::undef && should_hnf_cut()) r = hnf_cut(); if (r == lia_move::undef && should_gomory_cut()) r = gomory(lia).get_gomory_cuts(2); diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index ff5ff973e..2e7b2ef8c 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -864,7 +864,7 @@ namespace lp { } lp_status lar_solver::maximize_term(unsigned j, - impq& term_max) { + impq& term_max, bool fix_int_cols) { TRACE(lar_solver, print_values(tout);); SASSERT(get_core_solver().m_r_solver.calc_current_x_is_feasible_include_non_basis()); lar_term term = get_term_to_maximize(j); @@ -879,6 +879,11 @@ namespace lp { return lp_status::UNBOUNDED; } + if (!fix_int_cols) { + set_status(lp_status::OPTIMAL); + return lp_status::OPTIMAL; + } + impq opt_val = term_max; bool change = false; diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h index a5d49dd33..79c729e5d 100644 --- a/src/math/lp/lar_solver.h +++ b/src/math/lp/lar_solver.h @@ -205,7 +205,9 @@ public: set_column_value(j, v); } - lp_status maximize_term(unsigned j_or_term, impq& term_max); + // fix_int_cols: after maximizing try to move the integer columns to integer values; + // pass false to keep the optimal (possibly fractional) vertex intact, e.g., for the largest cube test + lp_status maximize_term(unsigned j_or_term, impq& term_max, bool fix_int_cols); core_solver_pretty_printer pp(std::ostream& out) const; diff --git a/src/math/lp/lp_params_helper.pyg b/src/math/lp/lp_params_helper.pyg index c3d50ab86..e9604e04b 100644 --- a/src/math/lp/lp_params_helper.pyg +++ b/src/math/lp/lp_params_helper.pyg @@ -8,6 +8,8 @@ def_module_params(module_name='lp', ('dio_cuts_enable_hnf', BOOL, True, 'enable hnf cuts together with Diophantine cuts, only relevant when dioph_eq is true'), ('dio_ignore_big_nums', BOOL, True, 'Ignore the terms with big numbers in the Diophantine handler, only relevant when dioph_eq is true'), ('dio_calls_period', UINT, 1, 'Period of calling the Diophantine handler in the final_check()'), - ('dio_run_gcd', BOOL, False, 'Run the GCD heuristic if dio is on, if dio is disabled the option is not used'), + ('dio_run_gcd', BOOL, False, 'Run the GCD heuristic if dio is on, if dio is disabled the option is not used'), + ('lcube', BOOL, True, 'use the largest cube test for integer feasibility'), + ('lcube_flips', UINT, 16, 'maximal number of coordinate flips when repairing the rounded largest cube center, only relevant when lcube is true'), )) diff --git a/src/math/lp/lp_settings.cpp b/src/math/lp/lp_settings.cpp index 42d6d8ef6..63836aafa 100644 --- a/src/math/lp/lp_settings.cpp +++ b/src/math/lp/lp_settings.cpp @@ -43,5 +43,7 @@ void lp::lp_settings::updt_params(params_ref const& _p) { m_dio_ignore_big_nums = lp_p.dio_ignore_big_nums(); m_dio_calls_period = lp_p.dio_calls_period(); m_dio_run_gcd = lp_p.dio_run_gcd(); + m_lcube = lp_p.lcube(); + m_lcube_flips = lp_p.lcube_flips(); m_max_conflicts = p.max_conflicts(); } diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h index d2e79ea90..cb27b1628 100644 --- a/src/math/lp/lp_settings.h +++ b/src/math/lp/lp_settings.h @@ -112,6 +112,9 @@ struct statistics { unsigned m_gcd_conflicts = 0; unsigned m_cube_calls = 0; unsigned m_cube_success = 0; + unsigned m_lcube_calls = 0; + unsigned m_lcube_success = 0; + unsigned m_lcube_flip_success = 0; unsigned m_patches = 0; unsigned m_patches_success = 0; unsigned m_hnf_cutter_calls = 0; @@ -152,6 +155,9 @@ struct statistics { st.update("arith-gcd-conflict", m_gcd_conflicts); st.update("arith-cube-calls", m_cube_calls); st.update("arith-cube-success", m_cube_success); + st.update("arith-lcube-calls", m_lcube_calls); + st.update("arith-lcube-success", m_lcube_success); + st.update("arith-lcube-flip-success", m_lcube_flip_success); st.update("arith-patches", m_patches); st.update("arith-patches-success", m_patches_success); st.update("arith-hnf-calls", m_hnf_cutter_calls); @@ -258,7 +264,11 @@ private: bool m_dio_ignore_big_nums = false; unsigned m_dio_calls_period = 4; bool m_dio_run_gcd = true; + bool m_lcube = true; + unsigned m_lcube_flips = 16; public: + bool lcube() const { return m_lcube; } + unsigned lcube_flips() const { return m_lcube_flips; } unsigned dio_calls_period() const { return m_dio_calls_period; } unsigned & dio_calls_period() { return m_dio_calls_period; } bool print_external_var_name() const { return m_print_external_var_name; } diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index faf367396..d3ac4840a 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -4029,7 +4029,7 @@ public: if (!lp().is_feasible() || lp().has_changed_columns()) make_feasible(); vi = get_lpvar(v); - auto st = lp().maximize_term(vi, term_max); + auto st = lp().maximize_term(vi, term_max, /*fix_int_cols*/ true); if (has_int() && lp().has_inf_int()) { st = lp::lp_status::FEASIBLE; lp().restore_x(); diff --git a/src/test/CMakeLists.txt b/src/test/CMakeLists.txt index 39050620a..ef5902c0c 100644 --- a/src/test/CMakeLists.txt +++ b/src/test/CMakeLists.txt @@ -79,6 +79,7 @@ add_executable(test-z3 "${CMAKE_CURRENT_BINARY_DIR}/install_tactic.cpp" interval.cpp karr.cpp + lcube.cpp list.cpp main.cpp map.cpp diff --git a/src/test/lcube.cpp b/src/test/lcube.cpp new file mode 100644 index 000000000..06acbf0fd --- /dev/null +++ b/src/test/lcube.cpp @@ -0,0 +1,261 @@ +/*++ + Copyright (c) 2020 Microsoft Corporation + + Module Name: + + lcube.cpp + + Abstract: + + Tests for the largest cube test of Bromberger and Weidenbach + (Fast Cube Tests for LIA Constraint Solving, IJCAR 2016), + implemented in int_cube::find_largest_cube(). + + This file lives directly under src/test (not src/test/lp) so that the + scripts/mk_make.py build, which only compiles the top-level test + directory, links tst_lcube(). + + Author: + + Lev Nachmanson (levnach) + + --*/ + +#include +#include +#include + +#include "util/debug.h" +#include "util/params.h" +#include "math/lp/int_cube.h" +#include "math/lp/int_solver.h" +#include "math/lp/lar_solver.h" +#include "math/lp/numeric_pair.h" + +namespace lp { + +// tests for the largest cube test of Bromberger and Weidenbach +namespace lcube_test { + + struct ineq { + vector> m_coeffs; + lconstraint_kind m_kind; + mpq m_rs; + }; + + // builds for every inequality a term with the bound and solves + static void setup(lar_solver& solver, const vector& ineqs, svector* term_columns = nullptr) { + unsigned term_ext = 1000; + for (const auto& in : ineqs) { + unsigned t = solver.add_term(in.m_coeffs, term_ext++); + solver.add_var_bound(t, in.m_kind, in.m_rs); + if (term_columns) + term_columns->push_back(t); + } + auto st = solver.solve(); + VERIFY(st == lp_status::OPTIMAL || st == lp_status::FEASIBLE); + } + + static void verify_model(const lar_solver& solver, const vector& ineqs) { + for (const auto& in : ineqs) { + impq v; + for (const auto& p : in.m_coeffs) + v += solver.get_column_value(p.second) * p.first; + switch (in.m_kind) { + case lconstraint_kind::LE: VERIFY(v <= impq(in.m_rs)); break; + case lconstraint_kind::GE: VERIFY(v >= impq(in.m_rs)); break; + default: VERIFY(false); + } + } + } + + static void verify_int_values(const lar_solver& solver, std::initializer_list vars) { + for (unsigned j : vars) + VERIFY(solver.get_column_value(j).is_int()); + } + + // The example of Bromberger and Weidenbach: 3x1 - x2 <= 0, -2x1 - x2 <= -2, -2x1 + x2 <= 1. + // The largest cube is smaller than the unit cube, the rounded center is not a solution, + // no coordinate flip repairs it (the only integer solution lies outside the lattice cell + // of the center): expect undef and an intact solver state. + static void test_paper_example_undef() { + std::cout << "lcube: paper example, expecting undef\n"; + lar_solver solver; + unsigned x1 = solver.add_named_var(0, true, "x1"); + unsigned x2 = solver.add_named_var(1, true, "x2"); + vector ineqs; + ineqs.push_back(ineq{{{mpq(3), x1}, {mpq(-1), x2}}, lconstraint_kind::LE, mpq(0)}); + ineqs.push_back(ineq{{{mpq(-2), x1}, {mpq(-1), x2}}, lconstraint_kind::LE, mpq(-2)}); + ineqs.push_back(ineq{{{mpq(-2), x1}, {mpq(1), x2}}, lconstraint_kind::LE, mpq(1)}); + setup(solver, ineqs); + int_solver i_s(solver); + solver.set_int_solver(&i_s); + lia_move m = int_cube(i_s).find_largest_cube(); + std::cout << "lcube returned " << lia_move_to_string(m) << "\n"; + VERIFY(m == lia_move::undef); + VERIFY(solver.ax_is_correct()); + } + + // 3x1 - x2 <= 0, -2x1 - x2 <= -1, -2x1 + x2 <= 1: the largest cube has + // edge 4/17 with the center (3/17, 1) that rounds to the solution (0, 1), + // while the unit cube test fails: the largest cube test is stronger here. + static void test_beats_unit_cube() { + std::cout << "lcube: beating the unit cube test\n"; + lar_solver solver; + unsigned x1 = solver.add_named_var(0, true, "x1"); + unsigned x2 = solver.add_named_var(1, true, "x2"); + vector ineqs; + ineqs.push_back(ineq{{{mpq(3), x1}, {mpq(-1), x2}}, lconstraint_kind::LE, mpq(0)}); + ineqs.push_back(ineq{{{mpq(-2), x1}, {mpq(-1), x2}}, lconstraint_kind::LE, mpq(-1)}); + ineqs.push_back(ineq{{{mpq(-2), x1}, {mpq(1), x2}}, lconstraint_kind::LE, mpq(1)}); + svector tcols; + setup(solver, ineqs, &tcols); + // move the solution to a fractional feasible point: the cube tests only + // run when the current solution has fractional integer variables + solver.set_column_value_test(x1, impq(mpq(1, 4))); + solver.set_column_value_test(x2, impq(mpq(5, 4))); + solver.set_column_value_test(tcols[0], impq(mpq(-1, 2))); + solver.set_column_value_test(tcols[1], impq(mpq(-7, 4))); + solver.set_column_value_test(tcols[2], impq(mpq(3, 4))); + int_solver i_s(solver); + solver.set_int_solver(&i_s); + lia_move m = int_cube(i_s)(); + std::cout << "unit cube returned " << lia_move_to_string(m) << "\n"; + VERIFY(m == lia_move::undef); + m = int_cube(i_s).find_largest_cube(); + std::cout << "lcube returned " << lia_move_to_string(m) << "\n"; + VERIFY(m == lia_move::sat); + verify_int_values(solver, {x1, x2}); + verify_model(solver, ineqs); + } + + // 9/10 <= x + y + r <= 11/10, -11/10 <= x - y + r <= 1/10, 0 <= r <= 1/10, + // x, y integer, r real. The real variable keeps the terms and their bounds + // non-integer. A fractional center, e.g. (1/2, 1/2), rounds to an infeasible + // point that is repaired by flipping one coordinate: expect sat. + static void test_flip_repair() { + std::cout << "lcube: rounding repair\n"; + lar_solver solver; + unsigned x = solver.add_named_var(0, true, "x"); + unsigned y = solver.add_named_var(1, true, "y"); + unsigned r = solver.add_named_var(2, false, "r"); + solver.add_var_bound(r, lconstraint_kind::GE, mpq(0)); + solver.add_var_bound(r, lconstraint_kind::LE, mpq(1, 10)); + vector ineqs; + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(1), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(9, 10)}); + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(1), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(11, 10)}); + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-1), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(-11, 10)}); + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-1), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(1, 10)}); + setup(solver, ineqs); + int_solver i_s(solver); + solver.set_int_solver(&i_s); + lia_move m = int_cube(i_s).find_largest_cube(); + std::cout << "lcube returned " << lia_move_to_string(m) + << ", flip successes: " << solver.settings().stats().m_lcube_flip_success << "\n"; + VERIFY(m == lia_move::sat); + verify_int_values(solver, {x, y}); + verify_model(solver, ineqs); + } + + // 3x + 5y >= 7 alone has infinite lattice width: the edge length is + // unbounded and any cube center of edge >= 1 rounds to a solution. + static void test_infinite_lattice_width() { + std::cout << "lcube: infinite lattice width\n"; + lar_solver solver; + unsigned x = solver.add_named_var(0, true, "x"); + unsigned y = solver.add_named_var(1, true, "y"); + vector ineqs; + ineqs.push_back(ineq{{{mpq(3), x}, {mpq(5), y}}, lconstraint_kind::GE, mpq(7)}); + setup(solver, ineqs); + int_solver i_s(solver); + solver.set_int_solver(&i_s); + lia_move m = int_cube(i_s).find_largest_cube(); + std::cout << "lcube returned " << lia_move_to_string(m) << "\n"; + VERIFY(m == lia_move::sat); + verify_int_values(solver, {x, y}); + verify_model(solver, ineqs); + } + + // 0 <= x + 2y + r <= 8, 0 <= x - 2y + r <= 8: the maximal edge length is + // 8/3 >= 1, so the rounded center is guaranteed to be a solution. + static void test_edge_at_least_one() { + std::cout << "lcube: edge length at least 1\n"; + lar_solver solver; + unsigned x = solver.add_named_var(0, true, "x"); + unsigned y = solver.add_named_var(1, true, "y"); + unsigned r = solver.add_named_var(2, false, "r"); + solver.add_var_bound(r, lconstraint_kind::GE, mpq(0)); + solver.add_var_bound(r, lconstraint_kind::LE, mpq(1, 10)); + vector ineqs; + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(2), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(0)}); + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(2), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(8)}); + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-2), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(0)}); + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-2), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(8)}); + setup(solver, ineqs); + int_solver i_s(solver); + solver.set_int_solver(&i_s); + lia_move m = int_cube(i_s).find_largest_cube(); + std::cout << "lcube returned " << lia_move_to_string(m) << "\n"; + VERIFY(m == lia_move::sat); + verify_int_values(solver, {x, y}); + verify_model(solver, ineqs); + } + + // runs the flip-repair instance through int_solver::check() with the + // lp.lcube parameter set and the cube period lowered to 1: verifies the + // dispatch and the parameter plumbing + static void test_dispatch() { + std::cout << "lcube: dispatch through int_solver::check\n"; + lar_solver solver; + params_ref p; + p.set_bool("lcube", true); + solver.settings().updt_params(p); + VERIFY(solver.settings().lcube()); + solver.settings().m_int_find_cube_period = 1; + unsigned x = solver.add_named_var(0, true, "x"); + unsigned y = solver.add_named_var(1, true, "y"); + unsigned r = solver.add_named_var(2, false, "r"); + solver.add_var_bound(r, lconstraint_kind::GE, mpq(0)); + solver.add_var_bound(r, lconstraint_kind::LE, mpq(1, 10)); + vector ineqs; + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(1), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(9, 10)}); + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(1), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(11, 10)}); + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-1), y}, {mpq(1), r}}, lconstraint_kind::GE, mpq(-11, 10)}); + ineqs.push_back(ineq{{{mpq(1), x}, {mpq(-1), y}, {mpq(1), r}}, lconstraint_kind::LE, mpq(1, 10)}); + svector tcols; + setup(solver, ineqs, &tcols); + // a fractional feasible point, so that check() does not return sat right away + solver.set_column_value_test(x, impq(mpq(1, 2))); + solver.set_column_value_test(y, impq(mpq(1, 2))); + solver.set_column_value_test(r, impq(mpq(1, 10))); + solver.set_column_value_test(tcols[0], impq(mpq(11, 10))); + solver.set_column_value_test(tcols[1], impq(mpq(11, 10))); + solver.set_column_value_test(tcols[2], impq(mpq(1, 10))); + solver.set_column_value_test(tcols[3], impq(mpq(1, 10))); + int_solver i_s(solver); + solver.set_int_solver(&i_s); + explanation ex; + lia_move m = i_s.check(&ex); + std::cout << "check returned " << lia_move_to_string(m) + << ", lcube calls: " << solver.settings().stats().m_lcube_calls << "\n"; + VERIFY(m == lia_move::sat); + VERIFY(solver.settings().stats().m_lcube_calls >= 1); + verify_int_values(solver, {x, y}); + verify_model(solver, ineqs); + } + + static void run() { + test_paper_example_undef(); + test_beats_unit_cube(); + test_flip_repair(); + test_infinite_lattice_width(); + test_edge_at_least_one(); + test_dispatch(); + std::cout << "lcube tests passed\n"; + } +} // namespace lcube_test +} // namespace lp + +void tst_lcube() { + lp::lcube_test::run(); +} diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 7a2b27155..7eaef8543 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -1645,7 +1645,7 @@ void test_maximize_term() { lia_move lm = i_solver.check(&ex); VERIFY(lm == lia_move::sat); impq term_max; - lp_status st = solver.maximize_term(term_2x_pl_2y, term_max); + lp_status st = solver.maximize_term(term_2x_pl_2y, term_max, /*fix_int_cols*/ true); std::cout << "status = " << lp_status_to_string(st) << std::endl; std::cout << "term_max = " << term_max << std::endl; diff --git a/src/test/main.cpp b/src/test/main.cpp index 1727cb9dc..4a5239f19 100644 --- a/src/test/main.cpp +++ b/src/test/main.cpp @@ -193,7 +193,8 @@ X(ho_matcher) \ X(finite_set) \ X(finite_set_rewriter) \ - X(fpa) + X(fpa) \ + X(lcube) #define FOR_EACH_TEST(X, X_ARGV) \ FOR_EACH_ALL_TEST(X, X_ARGV) \ diff --git a/src/util/mpz.cpp b/src/util/mpz.cpp index 10c345b0c..35de2deb0 100644 --- a/src/util/mpz.cpp +++ b/src/util/mpz.cpp @@ -1904,8 +1904,11 @@ std::string mpz_manager::to_string(mpz const & a) const { template unsigned mpz_manager::hash(mpz const & a) { - if (is_small(a)) - return ::abs(a.m_val); + if (is_small(a)) { + // compute abs in unsigned arithmetic: ::abs(INT_MIN) is undefined + unsigned u = static_cast(a.m_val); + return a.m_val < 0 ? 0u - u : u; + } #ifndef _MP_GMP unsigned sz = size(a); if (sz == 1)