From 862dc91cb218ad5ae782a766af170e2ecf00b9c9 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Mon, 23 Dec 2024 18:01:57 -1000 Subject: [PATCH] work on incremental dio Signed-off-by: Lev Nachmanson --- src/math/lp/dioph_eq.cpp | 137 +++++++++++++++++++++++--------- src/math/lp/lar_solver.cpp | 11 +++ src/math/lp/lar_solver.h | 7 ++ src/math/lp/static_matrix.cpp | 2 + src/math/lp/static_matrix.h | 7 +- src/math/lp/static_matrix_def.h | 39 ++++++++- src/test/lp/lp.cpp | 38 +++++++-- 7 files changed, 189 insertions(+), 52 deletions(-) diff --git a/src/math/lp/dioph_eq.cpp b/src/math/lp/dioph_eq.cpp index 9ed838649..3e18a75a4 100644 --- a/src/math/lp/dioph_eq.cpp +++ b/src/math/lp/dioph_eq.cpp @@ -125,7 +125,7 @@ namespace lp { m_c += t.c(); return *this; } - }; + }; std::ostream& print_S(std::ostream& out) { out << "S:\n"; @@ -204,7 +204,7 @@ namespace lp { NO_S_NO_F }; struct entry { - lar_term m_l; + //lar_term m_l; the term is taken from matrix m_l_matrix of the index entry mpq m_c; // the constant of the term, the term is taken from the row of // m_e_matrix with the same index as the entry entry_status m_entry_status; @@ -214,6 +214,7 @@ namespace lp { std_vector m_entries; // the terms are stored in m_A and m_c static_matrix m_e_matrix; // the rows of the matrix are the terms, + static_matrix m_l_matrix; // the rows of the matrix are the l_terms : providing the certificate to the entries // without the constant part int_solver& lia; lar_solver& lra; @@ -234,12 +235,12 @@ namespace lp { // {coeff*lar.get_term(k)}) std_vector m_k2s; - std_vector m_fresh_definitions; // seems only needed in the debug - // version in remove_fresh_vars + std_vector m_fresh_definitions; unsigned m_conflict_index = -1; // m_entries[m_conflict_index] gives the conflict unsigned m_max_number_of_iterations = 100; unsigned m_number_of_iterations; + std_vector> m_columns_to_entries; // m_columnn_to_terms[j] is the set of of all k such that m_entry[k].m_l depends on j struct branch { unsigned m_j = UINT_MAX; mpq m_rs; @@ -277,12 +278,57 @@ namespace lp { } }; + struct undo_add_term : public trail { + imp& m_s; + lar_term m_t; + undo_add_term(imp& s, const lar_term &t): m_s(s), m_t(t) { + + } + void undo () { + NOT_IMPLEMENTED_YET(); + } + }; + + + struct undo_add_column_bound : public trail { + imp& m_s; + unsigned m_j; // the column that has been added + undo_add_column_bound(imp& s, unsigned j) : m_s(s), m_j(j) {} + + void undo() override { + NOT_IMPLEMENTED_YET(); + } + }; + + std_vector m_added_terms; + std_vector m_branch_stats; std_vector m_branch_stack; std_vector m_explanation_of_branches; + void add_term_delegate(const lar_term* t) { + unsigned j = t->j(); + if (!lra.column_is_int(j) || !lra.column_has_term(j)) + return; + if (!all_vars_are_int(*t)) { + TRACE("dioph_eq", tout << "not all vars are integrall\n";); + return; + } + m_added_terms.push_back(t); + auto undo = undo_add_term(*this, *t); + lra.trail().push(undo); + } + + void add_column_bound_delegate(unsigned j) { + if (!lra.column_is_int(j)) + return; + auto undo = undo_add_column_bound(*this, j); + lra.trail().push(undo) ; + } public: - imp(int_solver& lia, lar_solver& lra) : lia(lia), lra(lra) {} + imp(int_solver& lia, lar_solver& lra) : lia(lia), lra(lra) { + lra.register_add_term_delegate([this](const lar_term*t){add_term_delegate(t);}); + } term_o get_term_from_entry(unsigned i) const { term_o t; for (const auto& p : m_e_matrix.m_rows[i]) { @@ -304,11 +350,16 @@ namespace lp { // the term has form sum(a_i*x_i) - t.j() = 0, void fill_entry(const lar_term& t) { TRACE("dioph_eq", print_lar_term_L(t, tout) << std::endl;); - entry te = {lar_term(t.j()), mpq(0), entry_status::F}; + entry te = { mpq(0), entry_status::F}; unsigned entry_index = (unsigned) m_entries.size(); m_f.push_back(entry_index); m_entries.push_back(te); entry& e = m_entries.back(); + SASSERT(m_l_matrix.row_count() == m_e_matrix.row_count()); +// fill m_l_matrix row + m_l_matrix.add_row(); + m_l_matrix.add_new_element(entry_index, t.j(), mpq(1)); +// fill E-entry m_e_matrix.add_row(); SASSERT(m_e_matrix.row_count() == m_entries.size()); @@ -356,6 +407,11 @@ namespace lp { m_number_of_iterations = 0; m_branch_stack.clear(); m_lra_level = 0; + for (const lar_term* t: m_added_terms) { + fill_entry(*t); + } + m_added_terms.clear(); + /* for (unsigned j = 0; j < lra.column_count(); j++) { if (!lra.column_is_int(j) || !lra.column_has_term(j)) continue; @@ -365,7 +421,7 @@ namespace lp { continue; } fill_entry(t); - } + }*/ } template @@ -437,7 +493,11 @@ namespace lp { p.coeff() /= g; } m_entries[ei].m_c = c_g; - e.m_l *= (1 / g); + // e.m_l *= (1 / g); + for (auto& p : m_l_matrix.m_rows[ei]) { + p.coeff() /= g; + } + TRACE("dioph_eq", tout << "ep_m_e:"; print_entry(ei, tout) << std::endl;); SASSERT(entry_invariant(ei)); @@ -511,13 +571,22 @@ namespace lp { q.push(j); } m_c += coeff * e.m_c; - m_tmp_l += coeff * e.m_l; + + m_tmp_l += coeff * l_term_from_row(k); // improve later TRACE("dioph_eq", tout << "after subs k:" << k << "\n"; print_term_o(create_term_from_ind_c(), tout) << std::endl; tout << "m_tmp_l:{"; print_lar_term_L(m_tmp_l, tout); tout << "}, opened:"; print_ml(m_tmp_l, tout) << std::endl;); } + lar_term l_term_from_row(unsigned k) const { + lar_term ret; + for (const auto & p: m_l_matrix.m_rows[k]) + ret.add_monomial(p.coeff(), p.var()); + + return ret; + } + term_o create_term_from_l(const lar_term& l) { term_o ret; for (const auto& p : l) { @@ -718,8 +787,8 @@ namespace lp { // returns true only on a conflict. // m_indexed_work_vector contains the coefficients of the term // m_c contains the constant term - // m_tmp_l is the linear combination of the equations that removs the - // substituted variablse returns true iff the conflict is found + // m_tmp_l is the linear combination of the equations that removes the + // substituted variables, returns true iff the conflict is found bool tighten_bounds_for_non_trivial_gcd(const mpq& g, unsigned j, bool is_upper) { mpq rs; @@ -781,7 +850,7 @@ namespace lp { return true; } - u_dependency* explain_fixed_in_meta_term(const lar_term& t) { + template u_dependency* explain_fixed_in_meta_term (const T& t) { return explain_fixed(open_ml(t)); } @@ -831,16 +900,6 @@ namespace lp { return lia_move::undef; } - struct undo_fixed : public trail { - lar_solver& m_lra; - unsigned m_j; // the variable that became fixed - undo_fixed(lar_solver& lra, unsigned j) : m_lra(lra), m_j(j) {} - - void undo() override { - NOT_IMPLEMENTED_YET(); - } - }; - void collect_evidence() { lra.get_infeasibility_explanation(m_infeas_explanation); for (const auto& p : m_infeas_explanation) { @@ -894,7 +953,7 @@ namespace lp { tout << "fixed j:" << j <<", was substited by "; print_entry(m_k2s[j], tout);); if (check_fixing(j) == lia_move::conflict) { auto& ep = m_entries[m_k2s[j]]; - for (auto ci : lra.flatten(explain_fixed_in_meta_term(ep.m_l))) { + for (auto ci : lra.flatten(explain_fixed_in_meta_term(m_l_matrix.m_rows[m_k2s[j]]))) { m_explanation_of_branches.push_back(ci); } return lia_move::conflict; @@ -929,11 +988,11 @@ namespace lp { } TRACE("dio_br", lra.print_column_info(b.m_j, tout) <<"add bound" << std::endl;); if (lra.column_is_fixed(b.m_j)) { - unsigned local_mj; - if (! m_var_register.external_is_used(b.m_j, local_mj)) + unsigned local_bj; + if (! m_var_register.external_is_used(b.m_j, local_bj)) return lia_move::undef; - if (fix_var(lar_solver_to_local(b.m_j)) == lia_move::conflict) { + if (fix_var(local_bj) == lia_move::conflict) { TRACE("dio_br", tout << "conflict in fix_var" << std::endl;) ; return lia_move::conflict; } @@ -1197,14 +1256,15 @@ namespace lp { print_entry(i, tout) << std::endl;); m_entries[i].m_c -= j_sign * coeff * e.m_c; m_e_matrix.pivot_row_to_row_given_cell_with_sign(ei, c, j, j_sign); - m_entries[i].m_l -= j_sign * coeff * e.m_l; + //m_entries[i].m_l -= j_sign * coeff * e.m_l; + m_l_matrix.add_rows( -j_sign*coeff, ei, i); TRACE("dioph_eq", tout << "after pivoting c_row:"; print_entry(i, tout);); CTRACE( "dioph_eq", !entry_invariant(i), tout << "invariant delta:"; { const auto& e = m_entries[i]; print_term_o(get_term_from_entry(ei) - - fix_vars(open_ml(e.m_l)), + fix_vars(open_ml(m_l_matrix.m_rows[ei])), tout) << std::endl; }); @@ -1226,7 +1286,7 @@ namespace lp { const auto& e = m_entries[ei]; bool ret = term_to_lar_solver(remove_fresh_vars(get_term_from_entry(ei))) == - fix_vars(open_ml(e.m_l)); + fix_vars(open_ml(m_l_matrix.m_rows[ei])); if (ret) return true; TRACE("dioph_eq", tout << "get_term_from_entry(" << ei << "):"; @@ -1234,11 +1294,12 @@ namespace lp { tout << "remove_fresh_vars:"; print_term_o(remove_fresh_vars(get_term_from_entry(ei)), tout) << std::endl; - tout << "e.m_l:"; print_lar_term_L(e.m_l, tout) << std::endl; + tout << "e.m_l:"; print_lar_term_L(l_term_from_row(ei), tout) << std::endl; tout << "open_ml(e.m_l):"; - print_term_o(open_ml(e.m_l), tout) << std::endl; + print_term_o(open_ml(l_term_from_row(ei)), tout) << std::endl; tout << "fix_vars(open_ml(e.m_l)):"; - print_term_o(fix_vars(open_ml(e.m_l)), tout) << std::endl;); + print_term_o(fix_vars(open_ml(l_term_from_row(ei))), tout) << std::endl; + ); return false; } @@ -1276,7 +1337,7 @@ namespace lp { return print_term_o(opened_ml, out); } - term_o open_ml(const lar_term& ml) const { + template term_o open_ml(const T& ml) const { term_o r; for (const auto& p : ml) { r += p.coeff() * (lra.get_term(p.var()) - lar_term(p.var())); @@ -1319,7 +1380,7 @@ namespace lp { e.m_c = r; m_e_matrix.add_new_element(h, xt, ahk); - m_entries.push_back({lar_term(), q, entry_status::NO_S_NO_F}); + m_entries.push_back({q, entry_status::NO_S_NO_F}); m_e_matrix.add_new_element(fresh_row, xt, -mpq(1)); m_e_matrix.add_new_element(fresh_row, k, mpq(1)); for (unsigned i : m_indexed_work_vector.m_index) { @@ -1359,9 +1420,9 @@ namespace lp { // out << "\tstatus:" << (int)e.m_entry_status; if (need_print_dep) { out << "\tm_l:{"; - print_lar_term_L(e.m_l, out) << "}, "; - print_ml(e.m_l, out << "\n\tfixed expl of m_l:\n") << "\n"; - print_dep(out, explain_fixed_in_meta_term(e.m_l)) << ","; + print_lar_term_L(l_term_from_row(ei), out) << "}, "; + print_ml(l_term_from_row(ei), out << "\n\tfixed expl of m_l:\n") << "\n"; + print_dep(out, explain_fixed_in_meta_term(l_term_from_row(ei))) << ","; } switch (e.m_entry_status) { case entry_status::F: @@ -1436,7 +1497,7 @@ namespace lp { TRACE("dioph_eq", tout << "conflict:"; print_entry(m_conflict_index, tout, true) << std::endl;); auto& ep = m_entries[m_conflict_index]; - for (auto ci : lra.flatten(explain_fixed_in_meta_term(ep.m_l))) { + for (auto ci : lra.flatten(explain_fixed_in_meta_term(m_l_matrix.m_rows[m_conflict_index]))) { ex.push_back(ci); } TRACE("dioph_eq", lra.print_expl(tout, ex);); diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index 5818641f1..3ddc00334 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -4,6 +4,7 @@ */ #include "math/lp/lar_solver.h" #include "smt/params/smt_params_helper.hpp" +#include "lar_solver.h" namespace lp { @@ -1598,6 +1599,12 @@ namespace lp { bool lar_solver::external_is_used(unsigned v) const { return m_var_register.external_is_used(v); } + void lar_solver::register_add_term_delegate(const std::function& f) { + this->m_add_term_delegates.push_back(f); + } + void lar_solver::register_add_column_bound_delegate(const std::function& f) { + this->m_add_column_bound_delegates.push_back(f); + } void lar_solver::add_non_basic_var_to_core_fields(unsigned ext_j, bool is_int) { register_new_external_var(ext_j, is_int); @@ -1695,6 +1702,8 @@ namespace lp { lp_assert(m_var_register.size() == A_r().column_count()); if (m_need_register_terms) register_normalized_term(*t, A_r().column_count() - 1); + for (const auto & f: m_add_term_delegates) + f(t); return ret; } @@ -1742,6 +1751,8 @@ namespace lp { constraint_index lar_solver::add_var_bound(lpvar j, lconstraint_kind kind, const mpq& right_side) { constraint_index ci = mk_var_bound(j, kind, right_side); activate(ci); + for (const auto & f: m_add_column_bound_delegates) + f(j); return ci; } diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h index 2f39022c5..9e3480788 100644 --- a/src/math/lp/lar_solver.h +++ b/src/math/lp/lar_solver.h @@ -405,6 +405,13 @@ public: } } + void register_add_term_delegate(const std::function&); + void register_add_column_bound_delegate(const std::function&); + + private: + std_vector> m_add_term_delegates; + std_vector> m_add_column_bound_delegates; + public: bool external_is_used(unsigned) const; void pop(unsigned k); unsigned num_scopes() const { return m_trail.get_num_scopes(); } diff --git a/src/math/lp/static_matrix.cpp b/src/math/lp/static_matrix.cpp index f5c121eaa..83a18416c 100644 --- a/src/math/lp/static_matrix.cpp +++ b/src/math/lp/static_matrix.cpp @@ -56,5 +56,7 @@ template bool lp::static_matrix >::pivot_row_ template void lp::static_matrix >::pivot_row_to_row_given_cell_with_sign(unsigned int, column_cell&, unsigned int, int); template void lp::static_matrix >::remove_element(vector, true, unsigned int>&, lp::row_cell&); template void lp::static_matrix::pivot_row_to_row_given_cell_with_sign(unsigned int, lp::row_cell&, unsigned int, int); +template void lp::static_matrix >::add_rows(mpq const&, unsigned int, unsigned int); +template void lp::static_matrix::add_rows(class rational const &,unsigned int,unsigned int); } diff --git a/src/math/lp/static_matrix.h b/src/math/lp/static_matrix.h index 51c4a1d4c..3504312b0 100644 --- a/src/math/lp/static_matrix.h +++ b/src/math/lp/static_matrix.h @@ -139,6 +139,8 @@ public: void add_new_element(unsigned i, unsigned j, const T & v); + // adds row i muliplied by coeff to row k + void add_rows(const mpq& coeff, unsigned i, unsigned k); void add_row() {m_rows.push_back(row_strip());} void add_column() { m_columns.push_back(column_strip()); @@ -218,8 +220,7 @@ public: return ret; } - void scan_row_to_work_vector(unsigned i); - + void scan_row_strip_to_work_vector(const row_strip & rvals); void clean_row_work_vector(unsigned i); @@ -300,8 +301,6 @@ public: // pivot row i to row ii bool pivot_row_to_row_given_cell(unsigned i, column_cell& c, unsigned j); void pivot_row_to_row_given_cell_with_sign(unsigned piv_row_index, column_cell& c, unsigned j, int j_sign); - void scan_row_ii_to_offset_vector(const row_strip & rvals); - void transpose_rows(unsigned i, unsigned ii) { auto t = m_rows[i]; m_rows[i] = m_rows[ii]; diff --git a/src/math/lp/static_matrix_def.h b/src/math/lp/static_matrix_def.h index a3b946782..9791900f5 100644 --- a/src/math/lp/static_matrix_def.h +++ b/src/math/lp/static_matrix_def.h @@ -42,7 +42,7 @@ void static_matrix::init_row_columns(unsigned m, unsigned n) { } -template void static_matrix::scan_row_ii_to_offset_vector(const row_strip & rvals) { +template void static_matrix:: scan_row_strip_to_work_vector(const row_strip & rvals) { for (unsigned j = 0; j < rvals.size(); j++) m_work_vector_of_row_offsets[rvals[j].var()] = j; } @@ -56,7 +56,7 @@ column_cell & c, unsigned pivot_col) { SASSERT(!is_zero(alpha)); auto & rowii = m_rows[ii]; remove_element(rowii, rowii[c.offset()]); - scan_row_ii_to_offset_vector(rowii); + scan_row_strip_to_work_vector(rowii); unsigned prev_size_ii = rowii.size(); // run over the pivot row and update row ii for (const auto & iv : m_rows[i]) { @@ -85,6 +85,37 @@ column_cell & c, unsigned pivot_col) { return !rowii.empty(); } + + template void static_matrix::add_rows(const mpq& alpha, unsigned i, unsigned k) { + lp_assert(i < row_count() && k < column_count() && i != k); + auto & rowk = m_rows[k]; + scan_row_strip_to_work_vector(rowk); + unsigned prev_size_k = rowk.size(); + // run over the pivot row and update row k + for (const auto & iv : m_rows[i]) { + unsigned j = iv.var(); + int j_offs = m_work_vector_of_row_offsets[j]; + if (j_offs == -1) { // it is a new element + T alv = alpha * iv.coeff(); + add_new_element(k, j, alv); + } + else { + addmul(rowk[j_offs].coeff(), iv.coeff(), alpha); + } + } + // clean the work vector + for (unsigned k = 0; k < prev_size_k; k++) { + m_work_vector_of_row_offsets[rowk[k].var()] = -1; + } + + // remove zeroes + for (unsigned k = rowk.size(); k-- > 0; ) { + if (is_zero(rowk[k].coeff())) + remove_element(rowk, rowk[k]); + } +} + + template inline void static_matrix::pivot_row_to_row_given_cell_with_sign(unsigned piv_row_index, column_cell& c, unsigned pivot_col, int pivot_sign) { @@ -94,7 +125,7 @@ column_cell& c, unsigned pivot_col, int pivot_sign) { SASSERT(!is_zero(alpha)); auto & rowii = m_rows[ii]; remove_element(rowii, rowii[c.offset()]); - scan_row_ii_to_offset_vector(rowii); + scan_row_strip_to_work_vector(rowii); unsigned prev_size_ii = rowii.size(); // run over the pivot row and update row ii for (const auto & iv : m_rows[piv_row_index]) { @@ -145,8 +176,8 @@ template void static_matrix::init_vector_of_row_o } template void static_matrix::init_empty_matrix(unsigned m, unsigned n) { - init_vector_of_row_offsets(); init_row_columns(m, n); + init_vector_of_row_offsets(); } template unsigned static_matrix::lowest_row_in_column(unsigned col) { diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 811e8482a..abf738b53 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -468,6 +468,7 @@ bool contains(std::string const &s, char const *pattern) { } void setup_args_parser(argument_parser &parser) { + parser.add_option_with_help_string("-add_rows", "test add_rows of static matrix"); parser.add_option_with_help_string("-monics", "test emonics"); parser.add_option_with_help_string("-nex_order", "test nex order"); parser.add_option_with_help_string("-nla_cn", "test cross nornmal form"); @@ -1728,6 +1729,33 @@ void test_gomory_cut() { test_gomory_cut_1(); } + void test_add_rows() { +// Create a static_matrix object + lp::static_matrix matrix; + matrix.init_empty_matrix(3, 3); + + // Populate the matrix with initial values + matrix.set(0, 0, mpq(1)); + matrix.set(0, 1, mpq(2)); + matrix.set(1, 0, mpq(3)); + matrix.set(1, 2, mpq(4)); + matrix.set(2, 1, mpq(5)); + matrix.set(2, 2, mpq(6)); + + // Perform add_rows operation + matrix.add_rows(mpq(2), 0, 1); // row 1 = row 1 + 2 * row 0 + + // Verify the results + SASSERT(matrix.get_elem(1, 0) == 5); // 3 + 2*1 + SASSERT(matrix.get_elem(1, 1) == 4); // 0 + 2*2 + SASSERT(matrix.get_elem(1, 2) == 4); // unchanged + + matrix.add_rows(mpq(-2), 0, 1); + SASSERT(matrix.get_elem(1, 0) == 3); // 5 - 2*1 + SASSERT(matrix.get_elem(1, 1) == 0); // 4 - 2*2 + SASSERT(matrix.get_elem(1, 2) == 4); // unchanged + } + void test_nla_order_lemma() { nla::test_order_lemma(); } void test_lp_local(int argn, char **argv) { @@ -1744,6 +1772,10 @@ void test_lp_local(int argn, char **argv) { } args_parser.print(); + if (args_parser.option_is_used("-add_rows")) { + test_add_rows(); + return finalize(0); + } if (args_parser.option_is_used("-monics")) { nla::test_monics(); return finalize(0); @@ -1842,12 +1874,6 @@ 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);