From df18885f97e8e364a81b664cc3fedb313c04b1df Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Sat, 24 Aug 2024 13:03:25 -1000 Subject: [PATCH] debug dio --- src/math/lp/dioph_eq.cpp | 218 ++++++++++++++++++++----------------- src/math/lp/int_solver.cpp | 11 +- src/math/lp/int_solver.h | 4 + src/math/lp/lar_term.h | 17 ++- src/test/lp/lp.cpp | 66 +++++++++++ 5 files changed, 206 insertions(+), 110 deletions(-) diff --git a/src/math/lp/dioph_eq.cpp b/src/math/lp/dioph_eq.cpp index cc0dd2519..1e6c25fce 100644 --- a/src/math/lp/dioph_eq.cpp +++ b/src/math/lp/dioph_eq.cpp @@ -23,20 +23,16 @@ namespace lp { const mpq& c() const { return m_c; } mpq& c() { return m_c; } term_o():m_c(0) {} - term_o& operator*=(mpq const& k) { - for (auto & t : m_coeffs) - t.m_value *= k; - return *this; - } void substitute(const term_o& t, unsigned term_column) { + SASSERT(!t.contains(term_column)); auto it = this->m_coeffs.find_core(term_column); if (it == nullptr) return; - mpq a = it->get_data().m_value; - this->m_coeffs.erase(term_column); + const mpq& a = it->get_data().m_value; for (auto p : t) { this->add_monomial(a * p.coeff(), p.j()); } this->c() += a * t.c(); + this->m_coeffs.erase(term_column); } }; @@ -114,39 +110,38 @@ namespace lp { imp(int_solver& lia, lar_solver& lra): lia(lia), lra(lra) {} void init() { - unsigned n_of_rows = static_cast(lra.r_basis().size()); + unsigned n_of_rows = lra.A_r().row_count(); unsigned skipped = 0; unsigned fn = 0; + + auto all_vars_are_int = [this](const auto& row) { + for (const auto& p : row) { + if (!lia.column_is_int(p.var())) { + return false; + } + } + return true; + }; for (unsigned i = 0; i < n_of_rows; i++) { auto & row = lra.get_row(i); + if (!all_vars_are_int(row)) continue; + const auto lcm = get_denominators_lcm(row); term_o t; - bool all_vars_are_int = true; - - for (const auto & p : row) { - if (!lia.column_is_int(p.var())){ - all_vars_are_int = false; - break; - } + for (const auto & p: row) + if (lia.is_fixed(p.var())) + t.c() += lia.lower_bound(p.var()).x; + else + t.add_monomial(lcm * p.coeff(), p.var()); + t.j() = i; //highjack this field to point to the original tableau row, TODO - is it used? + TRACE("dioph_eq", tout << "row = "; + lra.print_row(row, tout) << std::endl; + tout << "t:"; print_term_o(t, tout) << std::endl;); - } - - if (all_vars_are_int) { - term_o t; - const auto lcm = get_denominators_lcm(row); - for (const auto & p: row) - if (lia.is_fixed(p.var())) - t.c() += lia.lower_bound(p.var()).x; - else - t.add_monomial(lcm * p.coeff(), p.var()); - - t.j() = i; //hijach this field to point to the original tableau row - - unsigned lvar = static_cast(m_f.size()); - lar_term lt = lar_term(lvar); - eprime_pair pair = {t, lt}; - m_eprime.push_back(pair); - m_f.push_back(lvar); - } + unsigned lvar = static_cast(m_f.size()); + lar_term lt = lar_term(lvar); + eprime_pair pair = {t, lt}; + m_eprime.push_back(pair); + m_f.push_back(lvar); } } @@ -169,9 +164,9 @@ namespace lp { } if (g.is_one()) return true; + TRACE("dioph_eq", tout << "g:" << g << std::endl;); mpq new_c = ep.m_e.c() / g; if (new_c.is_int() == false) { - // conlict: todo - collect the conflict TRACE("dioph_eq", print_term_o(ep.m_e, tout << "conflict m_e:") << std::endl; tout << "g:" << g << std::endl; @@ -192,6 +187,10 @@ namespace lp { } ep.m_e.c() = new_c; ep.m_l *= (1/g); + TRACE("dioph_eq", + tout << "ep_m_e:"; print_term_o(ep.m_e, tout) << std::endl; + tout << "ep.ml:"; + print_lar_term_L(ep.m_l, tout ) << std::endl;); } return true; @@ -221,14 +220,14 @@ namespace lp { } void substitute_var_on_f(unsigned k, int k_sign, const term_o& k_subst, const lar_term& l_term) { - TRACE("dioph_eq", print_term_o(k_subst, tout<< k<< "->") << std::endl;); + TRACE("dioph_eq", print_term_o(k_subst, tout<< "x" << k<< " -> ") << std::endl;); for (unsigned e_index: m_f) { term_o& e = m_eprime[e_index].m_e; if (!e.contains(k)) continue; const mpq& k_coeff = e.get_coeff(k); TRACE("dioph_eq", print_term_o(e, tout << "before:") << std::endl; tout << "k_coeff:" << k_coeff << std::endl;); - m_eprime[e_index].m_l = k_sign*k_coeff*l_term + lar_term(e_index); // m_l is set to k_sign*e + e_sub + m_eprime[e_index].m_l += mpq(k_sign)*k_coeff*l_term; e.substitute(k_subst, k); TRACE("dioph_eq", print_term_o(e, tout << "after :") << std::endl; print_lar_term_L(m_eprime[e_index].m_l, tout) << std::endl;); @@ -261,10 +260,56 @@ namespace lp { if (p.j() == k) continue; t.add_monomial(-k_sign*p.coeff(), p.j()); } - t.c() = eh.c(); + t.c() = -k_sign*eh.c(); TRACE("dioph_eq", tout << "subst_term:"; print_term_o(t, tout) << std::endl;); return t; } + void eliminate_var(eprime_pair& e_pair, unsigned k, int k_sign) { + term_o t = get_term_to_subst(e_pair.m_e, k, k_sign); + substitute_var_on_f(k, k_sign, t, e_pair.m_l) ; + } + + void fresh_var_step(eprime_pair & e_pair, mpq& ahk, unsigned k) { +// step 7 from the paper + auto & eh = e_pair.m_e; + // xt is the fresh variable + unsigned xt = m_last_fresh_x_var--; + TRACE("dioph_eq", tout << "introduce fresh xt:" << "~x"<< fresh_index(xt) << std::endl; + tout << "eh:"; print_term_o(eh,tout) << std::endl;); + /* Let eh = sum (a_i*x_i) + c = 0. + For each i != k, let a_i = a_qi*ahk + a_ri, and let c = c_q * ahk + c_r + eh = ahk * (x_k + sum {a_qi*x_i) + c_q | i != k }) + sum {a_ri*x_i | i != k} + c_r = 0 + x_t = x_k + sum {a_qi*x_i) + c_q | i != k }, it will be xt_term + Then x_k -> - sum {a_qi*x_i) + c_q | i != k } + x_t, it will be subs_k + eh = ahk*x_t + ... + */ + term_o xt_term; + term_o k_subs; + mpq r, q = machine_div_rem(eh.c(), ahk, r); + xt_term.add_var(k); + xt_term.c() = q; + xt_term.add_monomial(mpq(-1), xt); + k_subs.add_var(xt); + k_subs.c() = - q; + + for (const auto & p: eh) { + if (p.j() == k) continue; + q = machine_div_rem(p.coeff(), ahk, r); + xt_term.add_monomial(q, p.j()); + k_subs.add_monomial(-q, p.j()); + } + TRACE("dioph_eq", tout << "xt_term:"; print_term_o(xt_term, tout) << std::endl; + tout << "k_subs:"; print_term_o(k_subs, tout) << std::endl;); + eprime_pair x_t_entry = {xt_term, lar_term()}; // 0 for m_l field + m_eprime.push_back(x_t_entry); + substitute_var_on_f(k, 1, k_subs, e_pair.m_l); + + // the term to eliminate the fresh variable + term_o xt_subs = xt_term.clone(); + xt_subs.add_monomial(mpq(1), xt); + TRACE("dioph_eq", tout << "xt_subs: x~"<< fresh_index(xt) << " -> "; print_term_o(xt_subs, tout) << std::endl;); + m_sigma.insert(xt, xt_subs); + } // this is the step 6 or 7 of the algorithm void rewrite_eqs() { auto eh_it = pick_eh(); @@ -273,61 +318,19 @@ namespace lp { print_term_o(m_eprime[*eh_it].m_e, tout << "\tm_e:") << "," << std::endl; print_lar_term_L(m_eprime[*eh_it].m_l, tout<< "\tm_l:") << "\n}"<< std::endl;); - term_o& eh = m_eprime[*eh_it].m_e; + term_o& eh = eprime_entry.m_e; auto [ahk, k, k_sign] = find_minimal_abs_coeff(eh); TRACE("dioph_eq", tout << "ahk:" << ahk << ", k:" << k << ", k_sign:" << k_sign << std::endl;); if (ahk.is_one()) { - // step 6 m_s.push_back(*eh_it); m_f.erase(eh_it); - term_o t = get_term_to_subst(eh, k, k_sign); - m_sigma.insert(k, t); - substitute_var_on_f(k, k_sign, t, eprime_entry.m_l) ; + eliminate_var(eprime_entry, k , k_sign); } else { - // step 7 - // the fresh variable - unsigned xt = m_last_fresh_x_var--; // xt is a fresh variable - ahk *= k_sign; // get the original value of ahk - TRACE("dioph_eq", tout << "introduce fresh xt:" << "~x"<< UINT_MAX -xt << std::endl; - tout << "eh:"; print_term_o(eh,tout) << std::endl;); - /* Let eh = sum (a_i*x_i) + c = 0. - For each i != k, let a_i = a_qi*ahk + a_ri, and let c = c_q * ahk + c_r - eh = ahk * (x_k + sum {a_qi*x_i) + c_q | i != k }) + sum {a_ri*x_i | i != k} + c_r = 0 - x_t = x_k + sum {a_qi*x_i) + c_q | i != k }, it will be x_t_term - Then x_k -> - sum {a_qi*x_i) + c_q | i != k } + x_t, it will be subs_k - eh = ahk*x_t + ... - */ - term_o x_t_term; - term_o k_subs; - mpq c_r; - mpq c_q = machine_div_rem(eh.c(), ahk, c_r); - x_t_term.add_var(k); - x_t_term.c() = c_q; - x_t_term.add_monomial(mpq(-1), xt); - k_subs.add_var(xt); - k_subs.c() = - c_q; - - for (const auto & p: eh) { - mpq a_q, a_r; - if (p.j() == k) continue; - a_q = machine_div_rem(p.coeff(), ahk, a_r); - x_t_term.add_monomial(a_q, p.j()); - k_subs.add_monomial(-a_q, p.j()); - } - TRACE("dioph_eq", tout << "x_t_term:"; print_term_o(x_t_term, tout) << std::endl; - tout << "k_subs:"; print_term_o(k_subs, tout) << std::endl;); - eprime_pair x_t_entry = {x_t_term, lar_term()}; // 0 for m_l field - m_eprime.push_back(x_t_entry); - substitute_var_on_f(k, 1, k_subs, eprime_entry.m_l); - - term_o x_t_subs = x_t_term.clone(); - x_t_subs.add_monomial(mpq(1), xt); - - TRACE("dioph_eq", tout << "x_t_subs:"; print_term_o(x_t_subs, tout) << std::endl;); - m_sigma.insert(xt, x_t_subs); + ahk *= k_sign; // restore the original coefficient + fresh_var_step(eprime_entry, ahk, k); } - } + void explain(lp::explanation& ex) { SASSERT(ex.empty()); auto & ep = m_eprime[m_conflict_index]; @@ -342,28 +345,40 @@ namespace lp { } TRACE("dioph_eq", lra.print_expl(tout, ex);); } + unsigned fresh_index(unsigned j) const {return UINT_MAX - j;} + struct fresh_lt { + bool operator()(unsigned v1, unsigned v2) const { return v1 > v2; } + }; + typedef heap fresh_queue; + void remove_fresh_variables(term_o& t) { - std::set seen_fresh_vars; + fresh_queue q(fresh_index(m_last_fresh_x_var)+1); + std::set processed; for (auto p : t) { auto j = p.j(); - if (is_fresh_var(j)) - seen_fresh_vars.insert(j); + if (is_fresh_var(j)) { + q.insert(fresh_index(j)); + processed.insert(j); + } } - CTRACE("dioph_eq_sub_terms", seen_fresh_vars.empty() == false, print_term_o(t, tout)<< std::endl); - while (!seen_fresh_vars.empty()) { - unsigned j = *seen_fresh_vars.begin(); - seen_fresh_vars.erase(j); - const term_o& ot = m_sigma.find(j); - TRACE("dioph_eq_sub_terms", tout << "ot:"; print_term_o(ot, tout) << std::endl;); - for (auto p : ot) - if (is_fresh_var(p.j())) - seen_fresh_vars.insert(p.j()); - t.substitute(ot, j); - - TRACE("dioph_eq_sub_terms", tout << "ot:"; print_term_o(ot, tout) << std::endl; - tout << "after sub:";print_term_o(t, tout) << std::endl;); + CTRACE("dioph_eq_fresh", q.empty() == false, print_term_o(t, tout)<< std::endl); + while (!q.empty()) { + unsigned j = q.erase_min(); + j = fresh_index(j); + const term_o& sub_t = m_sigma.find(j); + TRACE("dioph_eq_fresh", tout << "x~" << fresh_index(j) << "; sub_t:"; print_term_o(sub_t, tout) << std::endl;); + t.substitute(sub_t, j); + TRACE("dioph_eq_fresh", tout << "sub_t:"; print_term_o(sub_t, tout) << std::endl; + tout << "after sub:";print_term_o(t, tout) << std::endl;); + for (auto p : sub_t) + if (is_fresh_var(p.j())) { + if (processed.find(p.j()) != processed.end()) continue; + q.insert(fresh_index(p.j())); + processed.insert(j); + } } } + bool is_fresh_var(unsigned j) const { return j > lra.column_count(); } @@ -379,6 +394,7 @@ namespace lp { lia_move dioph_eq::check() { return m_imp->check(); } + void dioph_eq::explain(lp::explanation& ex) { m_imp->explain(ex); } diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp index fcf1e292f..999ac8b0b 100644 --- a/src/math/lp/int_solver.cpp +++ b/src/math/lp/int_solver.cpp @@ -228,9 +228,9 @@ namespace lp { return lia_move::undef; ++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_solve_dioph_eq()) r = solve_dioph_eq(); + // 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 && (true||should_solve_dioph_eq())) r = solve_dioph_eq(); if (r == lia_move::undef) lra.move_non_basic_columns_to_bounds(); if (r == lia_move::undef && should_hnf_cut()) r = hnf_cut(); @@ -879,6 +879,7 @@ namespace lp { bool int_solver::is_upper() const { return m_imp->m_upper; } bool& int_solver::is_upper() { return m_imp->m_upper; } explanation* int_solver::expl() { return m_imp->m_ex; } + void int_solver::set_expl(lp::explanation * ex) { m_imp->m_ex = ex; } bool int_solver::column_is_int_inf(unsigned j) const { return m_imp->column_is_int_inf(j); } @@ -895,5 +896,7 @@ namespace lp { bool int_solver::current_solution_is_inf_on_cut() const { return m_imp->current_solution_is_inf_on_cut(); } const impq & int_solver::lower_bound(unsigned j) const { return m_imp->lower_bound(j);} const impq & int_solver::upper_bound(unsigned j) const { return m_imp->upper_bound(j);} - + #if Z3DEBUG + lia_move int_solver::dio_test() {return m_imp->solve_dioph_eq();} + #endif } diff --git a/src/math/lp/int_solver.h b/src/math/lp/int_solver.h index 14d5c8e4e..e72445ca5 100644 --- a/src/math/lp/int_solver.h +++ b/src/math/lp/int_solver.h @@ -93,6 +93,10 @@ public: bool is_term(unsigned j) const; unsigned column_count() const; int select_int_infeasible_var(); + void set_expl(lp::explanation * ex); explanation * expl(); + #if Z3DEBUG + lia_move dio_test(); + #endif }; } diff --git a/src/math/lp/lar_term.h b/src/math/lp/lar_term.h index 5f5cdc2d2..b316daa9f 100644 --- a/src/math/lp/lar_term.h +++ b/src/math/lp/lar_term.h @@ -171,12 +171,19 @@ public: friend lar_term operator*(const mpq& k, const lar_term& term) { - lar_term result; - for (const auto& p : term) { - result.add_monomial(p.coeff()*k, p.j()); - } - return result; + lar_term r; + for (const auto& p : term) { + r.add_monomial(p.coeff()*k, p.j()); } + return r; + } + + lar_term& operator+=(const lar_term& a) { + for (const auto& p : a) { + add_monomial(p.coeff(), p.j()); + } + return *this; + } lar_term& operator*=(mpq const& k) { for (auto & t : m_coeffs) diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 84077b98e..811e8482a 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -491,6 +491,7 @@ void setup_args_parser(argument_parser &parser) { "-nla_blnt_fm", "test_basic_lemma_for_mon_neutral_from_factors_to_monomial"); parser.add_option_with_help_string("-hnf", "test hermite normal form"); + parser.add_option_with_help_string("-dio", "dioph equalities"); parser.add_option_with_help_string("-gomory", "gomory"); parser.add_option_with_help_string("-intd", "test integer_domain"); parser.add_option_with_help_string("-xyz_sample", @@ -1605,6 +1606,7 @@ void test_larger_generated_hnf() { std::cout << "test_larger_generated_rank_hnf passed" << std::endl; } #endif + void test_maximize_term() { std::cout << "test_maximize_term\n"; lar_solver solver; @@ -1646,6 +1648,64 @@ void test_maximize_term() { std::cout << "v[" << p.first << "] = " << p.second << std::endl; } } + +void test_dio() { + std::cout << "test dio\n"; + lar_solver solver; + int_solver i_solver(solver); + lp::explanation exp; + i_solver.set_expl(&exp); + unsigned _x1 = 0; + unsigned _x2 = 1; + unsigned _x3 = 2; + unsigned _fx_7 = 3; + unsigned _fx_17 = 4; +/* + 3x1 + 3x2 + 14x3 − 7 = 0 + 7x1 + 12x2 + 31x3 − 17 = 0 +*/ + lpvar x1 = solver.add_var(_x1, true); + lpvar x2 = solver.add_var(_x2, true); + lpvar x3 = solver.add_var(_x3, true); + lpvar fx_7 = solver.add_var(_fx_7, true); + lpvar fx_17 = solver.add_var(_fx_17, true); + vector> term_ls; + /* 3x1 + 3x2 + 14x3 − 7 */ + term_ls.push_back(std::pair(mpq(3), x1)); + term_ls.push_back(std::pair(mpq(3), x2)); + term_ls.push_back(std::pair(mpq(14), x3)); + term_ls.push_back(std::pair(mpq(-1), fx_7)); + for (auto & p: term_ls) { + p.first = -p.first; + } + unsigned t0 = solver.add_term(term_ls, 10); + term_ls.clear(); + /* 7x1 + 12x2 + 31x3 − 17 = 0*/ + term_ls.push_back(std::pair(mpq(7), x1)); + term_ls.push_back(std::pair(mpq(12), x2)); + term_ls.push_back(std::pair(mpq(31), x3)); + term_ls.push_back(std::pair(mpq(-1), fx_17)); + + for (auto & p: term_ls) { + p.first = -p.first; + } + unsigned t1 = solver.add_term(term_ls, 11); + + solver.add_var_bound(fx_7, LE, mpq(-7)); + solver.add_var_bound(fx_7, GE, mpq(-7)); + solver.add_var_bound(fx_17, LE, mpq(-17)); + solver.add_var_bound(fx_17, GE, mpq(-17)); + solver.add_var_bound(t0, LE, mpq(0)); + solver.add_var_bound(t0, GE, mpq(0)); + solver.add_var_bound(t1, LE, mpq(0)); + solver.add_var_bound(t1, GE, mpq(0)); +// solver.find_feasible_solution(); + //lp_assert(solver.get_status() == lp_status::OPTIMAL); + enable_trace("dioph_eq"); + enable_trace("dioph_eq_fresh"); + auto r = i_solver.dio_test(); + +} #ifdef Z3DEBUG void test_hnf() { test_larger_generated_hnf(); @@ -1782,6 +1842,12 @@ void test_lp_local(int argn, char **argv) { return finalize(0); } + if (args_parser.option_is_used("-dio")) { + test_dio(); + return finalize(0); + } + + if (args_parser.option_is_used("-gomory")) { test_gomory_cut(); return finalize(0);