diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index 4715d34be..83a8619e3 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -565,40 +565,6 @@ bool lar_solver::remove_from_basis(unsigned j) { return m_mpq_lar_core_solver.m_r_solver.remove_from_basis(j); } -// val is the new value to be assigned to x[j] -// return true iff can find a new basic column that would be feasible -// after the pivoting -bool lar_solver::remove_from_basis(unsigned basic_j, const mpq& val) { - SASSERT(is_base(basic_j)); - impq del(0); - const auto& slv = m_mpq_lar_core_solver.m_r_solver; - bool grow = val < get_column_value(basic_j).x; // grow = true means that the monomial of the pivoted var has to grow - for (auto &c : A_r().m_rows[row_of_basic_column(basic_j)]) { - lpvar j = c.var(); - if (j == basic_j) { - SASSERT(c.coeff().is_one()); - continue; - } - - bool can_pivot = column_is_free(j); - if (!can_pivot && - ((grow && slv.monoid_can_increase(c))|| (!grow && slv.monoid_can_decrease(c)))) { - if (del.is_zero()) - del = impq(val) - get_column_value(basic_j); - - impq j_val = get_column_value(j) - del / c.coeff(); - if (inside_bounds(j, j_val)) - can_pivot = true; - } - - if (can_pivot) { - pivot_column_tableau(c.var(), row_of_basic_column(basic_j)); - return true; - } - } - return false; -} - lar_term lar_solver::get_term_to_maximize(unsigned j_or_term) const { if (tv::is_term(j_or_term)) { return get_term(j_or_term); @@ -1292,7 +1258,8 @@ void lar_solver::get_rid_of_inf_eps() { mpq delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(mpq(1)); for (unsigned j = 0; j < number_of_vars(); j++) { auto & r = m_mpq_lar_core_solver.m_r_x[j]; - r = impq(r.x + delta * r.y); + if (!r.y.is_zero()) + r = impq(r.x + delta * r.y); } } diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h index 5047486e7..a3fb26bbe 100644 --- a/src/math/lp/lar_solver.h +++ b/src/math/lp/lar_solver.h @@ -260,7 +260,6 @@ class lar_solver : public column_namer { void update_delta_for_terms(const impq & delta, unsigned j, const vector&); void fill_vars_to_terms(vector> & vars_to_terms); bool remove_from_basis(unsigned); - bool remove_from_basis(unsigned, const mpq&); lar_term get_term_to_maximize(unsigned ext_j) const; bool sum_first_coords(const lar_term& t, mpq & val) const; void collect_rounded_rows_to_fix(); @@ -362,14 +361,13 @@ public: const ChangeReport& change_report) { if (is_base(j)) { TRACE("nla_solver", get_int_solver()->display_row_info(tout, row_of_basic_column(j)) << "\n";); - if (!remove_from_basis(j, val)) - return false; + remove_from_basis(j); } impq ival(val); if (is_blocked(j, ival)) return false; - + TRACE("nla_solver", tout << "not blocked\n";); impq delta = get_column_value(j) - ival; for (const auto &c : A_r().column(j)) { unsigned row_index = c.var(); diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp index 0884eaf08..8c5cc5071 100644 --- a/src/math/lp/nla_core.cpp +++ b/src/math/lp/nla_core.cpp @@ -1215,7 +1215,7 @@ bool core::var_breaks_correct_monic_as_factor(lpvar j, const monic& m) const { if (!val(var(m)).is_zero()) return true; - if (!val(j).is_zero()) // j was not zero: the new value does not matter - m must have another zero product + if (!val(j).is_zero()) // j was not zero: the new value does not matter - m must have another zero factor return false; // do we have another zero in m? for (lpvar k : m) { @@ -1287,29 +1287,38 @@ bool core::has_real(const monic& m) const { return true; return false; } + // returns true if the patching is blocking -bool core::patch_is_blocked(lpvar u, const monic& m, const lp::impq& ival) const { - SASSERT(m_to_refine.contains(m.var())); - if (m_cautious_patching && !m_lar_solver.inside_bounds(u, ival)) +bool core::patch_is_blocked(lpvar u, const monic& patched_m, const lp::impq& ival) const { + if (m_cautious_patching && + (!m_lar_solver.inside_bounds(u, ival) || (var_is_int(u) && ival.is_int() == false))) { + TRACE("nla_solver", tout << "u = " << u << " blocked, for feas or integr\n";); return true; // blocking + } if (var_breaks_correct_monic(u)) { TRACE("nla_solver", tout << "u = " << u << " blocked as used in a correct monomial\n";); return true; } - bool ret = u == m.var() || (m.contains_var(u) && var_breaks_correct_monic_as_factor(u, m)); + bool ret = u == patched_m.var() || (patched_m.contains_var(u) && var_breaks_correct_monic_as_factor(u, patched_m)); - TRACE("nla_solver", tout << "u = " << u << ", m = "; print_monic(m, tout) << + TRACE("nla_solver", tout << "u = " << u << ", patched_m = "; print_monic(patched_m, tout) << "ret = " << ret << "\n";); return ret; } -bool core::try_to_patch(lpvar k, const rational& v, const monic & m) { - auto is_blocked = [this, k, m](lpvar u, const lp::impq& v) - { return u != k && patch_is_blocked(u, m, v); }; +bool core::try_to_patch(lpvar patched_var, const rational& v, const monic & m) { + auto is_blocked = [this, patched_var, m](lpvar u, const lp::impq& iv) + { + if (!m_lar_solver.inside_bounds(u, iv)) + return true; + if (u == patched_var) + return false; + return patch_is_blocked(u, m, iv); + }; auto change_report = [this](lpvar u) { update_to_refine_of_var(u); }; - return m_lar_solver.try_to_patch(k, v, is_blocked, change_report); + return m_lar_solver.try_to_patch(patched_var, v, is_blocked, change_report); } bool in_power(const svector& vs, unsigned l) { @@ -1348,31 +1357,36 @@ void core::patch_monomial(lpvar j) { return; } + + // Now we try patching the factor variables. + TRACE("nla_solver", tout << " trying squares\n";); // handle perfect squares if (m.vars().size() == 2 && m.vars()[0] == m.vars()[1]) { rational root; if (v.is_perfect_square(root)) { lpvar k = m.vars()[0]; - if (!var_is_int(k) && - !var_breaks_correct_monic(k) && - (try_to_patch(k, root, m) || try_to_patch(k, -root, m)) - ) { + if (!var_breaks_correct_monic(k) && (try_to_patch(k, root, m) || try_to_patch(k, -root, m))) { + TRACE("nla_solver", tout << "patched square\n";); + return; } } + TRACE("nla_solver", tout << " cannot patch\n";); return; } - // We have v != abc. Let us suppose we patch b. Then b should - // be equal to v/ac = v/(abc/b) = b(v/abc) + + // We have v != abc, but we need to have v = abc. + // If we patch b then b should be equal to v/ac = v/(abc/b) = b(v/abc) if (!v.is_zero()) { rational r = val(j) / v; SASSERT(m.is_sorted()); - + TRACE("nla_solver", tout << "r = " << r << ", v = " << v << "\n";); for (unsigned l = 0; l < m.size(); l++) { lpvar k = m.vars()[l]; if (!in_power(m.vars(), l) && !var_is_int(k) && !var_breaks_correct_monic(k) && try_to_patch(k, r * val(k), m)) { // r * val(k) gives the right value of k + TRACE("nla_solver", tout << "patched j " << j << "\n";); SASSERT(mul_val(m) == var_val(m)); erase_from_to_refine(j); break; @@ -1385,13 +1399,15 @@ void core::patch_monomials_on_to_refine() { auto to_refine = m_to_refine.index(); // the rest of the function might change m_to_refine, so have to copy unsigned sz = to_refine.size(); - TRACE("nla_solver", tout << "sz = " << sz << "\n";); + unsigned start = random(); for (unsigned i = 0; i < sz; i++) { patch_monomial(to_refine[(start + i) % sz]); if (m_to_refine.size() == 0) break; } + TRACE("nla_solver", tout << "sz = " << sz << ", m_to_refine = " << m_to_refine.size() << + (sz > m_to_refine.size()? " less" : "same" ) << "\n";); } void core::patch_monomials() { @@ -1400,7 +1416,8 @@ void core::patch_monomials() { if (m_to_refine.size() == 0 || !m_nla_settings.expensive_patching()) { return; } - m_cautious_patching = false; // + NOT_IMPLEMENTED_YET(); + m_cautious_patching = false; patch_monomials_on_to_refine(); m_lar_solver.push(); save_tableau();