From b47e0cd4fb072646ba637f2a4a5523aed17d3d34 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Thu, 18 Apr 2019 15:21:40 -0700 Subject: [PATCH] fix tableau's Ax=b after cube heuristic Signed-off-by: Lev Nachmanson --- src/util/lp/int_solver.cpp | 13 ++++--- src/util/lp/lar_solver.cpp | 36 ++++++++++++++++--- src/util/lp/lar_solver.h | 18 +++++----- src/util/lp/lp_core_solver_base.h | 6 ++-- .../lp/lp_primal_core_solver_tableau_def.h | 4 +-- src/util/lp/nla_core.cpp | 2 +- 6 files changed, 55 insertions(+), 24 deletions(-) diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index 39d743258..b539947c8 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -248,25 +248,30 @@ lia_move int_solver::find_cube() { TRACE("cube", for (unsigned j = 0; j < m_lar_solver->A_r().column_count(); j++) display_column(tout, j); - m_lar_solver->print_terms(tout); + m_lar_solver->print_constraints(tout); ); - lar_solver::scoped_push _sp(*m_lar_solver); + m_lar_solver->push(); if (!tighten_terms_for_cube()) { + m_lar_solver->pop(); return lia_move::undef; } lp_status st = m_lar_solver->find_feasible_solution(); if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { TRACE("cube", tout << "cannot find a feasiblie solution";); - _sp.pop(); + m_lar_solver->pop(); m_lar_solver->move_non_basic_columns_to_bounds(); find_feasible_solution(); // it can happen that we found an integer solution here return !m_lar_solver->r_basis_has_inf_int()? lia_move::sat: lia_move::undef; } - _sp.pop(); + m_lar_solver->x_shifted_in_cube() = true; + m_lar_solver->pop(); m_lar_solver->round_to_integer_solution(); + m_lar_solver->set_status(lp_status::FEASIBLE); + lp_assert(is_feasible()); + TRACE("cube", tout << "success";); settings().st().m_cube_success++; TRACE("cube", tout << "sat with cube\n";); return lia_move::sat; diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index dc84855c7..4b2c717d9 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -27,7 +27,8 @@ void clear() {lp_assert(false); // not implemented } -lar_solver::lar_solver() : m_status(lp_status::UNKNOWN), +lar_solver::lar_solver() : m_x_shifted_in_cube(false), + m_status(lp_status::UNKNOWN), m_infeasible_column(-1), m_mpq_lar_core_solver(m_settings, *this), m_int_solver(nullptr), @@ -357,7 +358,6 @@ void lar_solver::shrink_inf_set_after_pop(unsigned n, int_set & set) { void lar_solver::pop(unsigned k) { - TRACE("int_solver", tout << "pop" << std::endl;); TRACE("lar_solver", tout << "k = " << k << std::endl;); restore_rounded_columns(); // if it is not done now, the basis changes and restore_rounded_columns would now work m_infeasible_column.pop(k); @@ -398,6 +398,7 @@ void lar_solver::pop(unsigned k) { m_settings.simplex_strategy() = m_simplex_strategy; lp_assert(sizes_are_correct()); lp_assert((!m_settings.use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); + lp_assert(x_shifted_in_cube() || ax_is_correct()); set_status(lp_status::UNKNOWN); } @@ -784,8 +785,9 @@ unsigned lar_solver::row_of_basic_column(unsigned j) const { bool lar_solver::ax_is_correct() const { for (unsigned i = 0; i < A_r().row_count(); i++) { - if (!row_is_correct(i)) + if (!row_is_correct(i)) { return false; + } } return true; } @@ -873,8 +875,27 @@ void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds_tableau( } } +void lar_solver::fix_Ax_b_on_row(unsigned i) { + unsigned bj = m_mpq_lar_core_solver.m_r_basis[i]; + auto& v = m_mpq_lar_core_solver.m_r_x[bj]; + v = zero_of_type>(); + for (const auto & c : A_r().m_rows[i]) { + if (c.var() != bj) + v -= c.coeff() * m_mpq_lar_core_solver.m_r_x[c.var()]; + } +} +void lar_solver::fix_Ax_b() { + for (unsigned i = 0; i < A_r().row_count(); i++) { + fix_Ax_b_on_row(i); + } + lp_assert(ax_is_correct()); +} void lar_solver::solve_with_core_solver() { + if (x_shifted_in_cube()) { + fix_Ax_b(); + x_shifted_in_cube() = false; + } if (!use_tableau()) add_last_rows_to_lu(m_mpq_lar_core_solver.m_r_solver); if (m_mpq_lar_core_solver.need_to_presolve_with_double_solver()) { @@ -1026,6 +1047,12 @@ bool lar_solver::all_constraints_hold() const { for (unsigned i = 0; i < m_constraints.size(); i++) { if (!constraint_holds(*m_constraints[i], var_map)) { + TRACE("lar_solver", + tout << "i = " << i << "; "; + print_constraint(m_constraints[i], tout) << "\n"; + for(auto p : m_constraints[i]->coeffs()) { + m_mpq_lar_core_solver.m_r_solver.print_column_info(p.second, tout); + }); return false; } } @@ -2232,11 +2259,11 @@ void lar_solver::round_to_integer_solution() { for (unsigned j = 0; j < column_count(); j++) { if (!column_is_int(j)) continue; if (column_corresponds_to_term(j)) continue; - TRACE("cube", m_int_solver->display_column(tout, j);); impq& v = m_mpq_lar_core_solver.m_r_x[j]; if (v.is_int()) continue; SASSERT(is_base(j)); + TRACE("cube", m_int_solver->display_column(tout, j);); impq flv = impq(floor(v)); auto del = flv - v; // del is negative if (del < - impq(mpq(1, 2))) { @@ -2246,6 +2273,7 @@ void lar_solver::round_to_integer_solution() { v = flv; } m_incorrect_columns.insert(j); + TRACE("cube", tout << "new val = " << v << "\n";); } } // return true if all y coords are zeroes diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index b8fe67ba2..2eee63430 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -88,6 +88,9 @@ class lar_solver : public column_namer { #endif //////////////////// fields ////////////////////////// + // it is set to true when during the cube heuristic the column values got shifted + // and in this case the tableau becames invalid + bool m_x_shifted_in_cube; lp_settings m_settings; lp_status m_status; stacked_value m_simplex_strategy; @@ -114,6 +117,10 @@ public: stacked_value m_term_count; vector m_terms; indexed_vector m_column_buffer; + // end of fields + + bool x_shifted_in_cube() const { return m_x_shifted_in_cube; } + bool& x_shifted_in_cube() { return m_x_shifted_in_cube; } unsigned terms_start_index() const { return m_terms_start_index; } const vector & terms() const { return m_terms; } @@ -347,15 +354,6 @@ public: void pop(unsigned k); - class scoped_push { - lar_solver& m_solver; - bool m_pop; - public: - scoped_push(lar_solver& s):m_solver(s), m_pop(true) { s.push(); } - ~scoped_push() { if (m_pop) m_solver.pop(); } - void pop() { SASSERT(m_pop); m_solver.pop(); m_pop = false; } - }; - vector get_all_constraint_indices() const; bool maximize_term_on_tableau(const lar_term & term, @@ -644,5 +642,7 @@ public: lar_term get_term_to_maximize(unsigned ext_j) const; void set_cut_strategy(unsigned cut_frequency); bool sum_first_coords(const lar_term& t, mpq & val) const; + void fix_Ax_b(); + void fix_Ax_b_on_row(unsigned); }; } diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index 7e909f3e6..64f6f5696 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -252,13 +252,11 @@ public: } bool below_bound(const X & x, const X & bound) const { - if (precise()) return x < bound; - return below_bound_numeric(x, bound, m_settings.primal_feasibility_tolerance); + return precise()? x < bound : below_bound_numeric(x, bound, m_settings.primal_feasibility_tolerance); } bool above_bound(const X & x, const X & bound) const { - if (precise()) return x > bound; - return above_bound_numeric(x, bound, m_settings.primal_feasibility_tolerance); + return precise()? x > bound : above_bound_numeric(x, bound, m_settings.primal_feasibility_tolerance); } bool x_below_low_bound(unsigned p) const { diff --git a/src/util/lp/lp_primal_core_solver_tableau_def.h b/src/util/lp/lp_primal_core_solver_tableau_def.h index 0d6f7ac92..a2f8367b9 100644 --- a/src/util/lp/lp_primal_core_solver_tableau_def.h +++ b/src/util/lp/lp_primal_core_solver_tableau_def.h @@ -115,9 +115,9 @@ unsigned lp_primal_core_solver::solve_with_tableau() { } if (this->m_settings.use_tableau_rows()) { one_iteration_tableau_rows(); - } - else + } else { one_iteration_tableau(); + } TRACE("lar_solver", tout << "one iteration tableau " << this->get_status() << "\n";); switch (this->get_status()) { case lp_status::OPTIMAL: // double check that we are at optimum diff --git a/src/util/lp/nla_core.cpp b/src/util/lp/nla_core.cpp index d25d1f09d..910414414 100644 --- a/src/util/lp/nla_core.cpp +++ b/src/util/lp/nla_core.cpp @@ -2165,7 +2165,7 @@ lbool core:: check(vector& l_vec) { TRACE("nla_solver", tout << "calls = " << settings().st().m_nla_calls << "\n";); m_lemma_vec = &l_vec; if (!(m_lar_solver.get_status() == lp::lp_status::OPTIMAL || m_lar_solver.get_status() == lp::lp_status::FEASIBLE )) { - TRACE("nla_solver", tout << "unknown because of the m_lar_solver.m_status = " << lp_status_to_string(m_lar_solver.get_status()) << "\n";); + TRACE("nla_solver", tout << "unknown because of the m_lar_solver.m_status = " << m_lar_solver.get_status() << "\n";); return l_undef; }