From 89eec4cc803b6bcabf52c2a99ec29f476d1f439f Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Mon, 20 Jan 2025 07:31:30 -1000 Subject: [PATCH] test dio Signed-off-by: Lev Nachmanson --- src/math/lp/dioph_eq.cpp | 378 +++++++++++++++++++++++-------------- src/math/lp/lar_solver.cpp | 5 +- src/math/lp/lar_solver.h | 3 +- 3 files changed, 246 insertions(+), 140 deletions(-) diff --git a/src/math/lp/dioph_eq.cpp b/src/math/lp/dioph_eq.cpp index 199ea0d9d..a3f0c5252 100644 --- a/src/math/lp/dioph_eq.cpp +++ b/src/math/lp/dioph_eq.cpp @@ -9,6 +9,7 @@ #include "math/lp/int_solver.h" #include "math/lp/lar_solver.h" #include "math/lp/lp_utils.h" +#include "math/lp/var_register_dio.h" /* Following paper: "A Practical Approach to Satisfiability Modulo Linear Integer Arithmetic" by Alberto Griggio(griggio@fbk.eu). @@ -212,7 +213,7 @@ namespace lp { entry_status m_entry_status; }; - var_register m_var_register; + var_register_dio m_var_register; 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, @@ -282,62 +283,44 @@ namespace lp { } }; - struct undo_add_term : public trail { - imp& m_s; - const lar_term* m_t; - unsigned m_l_column_count; - unsigned m_row_count; - unsigned m_e_column_count; - undo_add_term(imp& s, const lar_term *t): - m_s(s), - m_t(t), - m_l_column_count(m_s.m_l_matrix.column_count()), - m_row_count(m_s.m_l_matrix.row_count()), - m_e_column_count(m_s.m_e_matrix.column_count()) { - TRACE("dioph_eq", m_s.print_lar_term_L(*t, tout); tout << "t->j()=" << t->j() << std::endl;); - } - void undo () { - SASSERT(m_s.entries_are_ok()); - TRACE("dioph_eq", m_s.lra.print_term(*m_t, tout); tout << ", m_t->j() =" << m_t->j() << std::endl;); - if (!contains(m_s.m_active_terms, m_t)) { - for (int i = m_s.m_added_terms.size() - 1; i >= 0; --i) { - if (m_s.m_added_terms[i] != m_t) continue; - // the address is the same - if (i != m_s.m_added_terms.size() -1) - m_s.m_added_terms[i] = m_s.m_added_terms.back(); - m_s.m_added_terms.pop_back(); - SASSERT(m_s.m_l_matrix.row_count() == m_row_count && m_s.m_l_matrix.column_count() == m_l_column_count - && m_s.m_e_matrix.column_count() == m_e_column_count); - return; // all is done since the term has not made it to entries, etc - } + void remove_term_callback(const lar_term *t) { + TRACE("d_undo", tout << "t:"<< t<<", t->j():"<< t->j() << std::endl;); + TRACE("dioph_eq", lra.print_term(*t, tout); tout << ", t->j() =" << t->j() << std::endl;); + if (!contains(m_active_terms, t)) { + for (int i = m_added_terms.size() - 1; i >= 0; --i) { + if (m_added_terms[i] != t) continue; + if (i != m_added_terms.size() -1) + m_added_terms[i] = m_added_terms.back(); + m_added_terms.pop_back(); + break; // all is done since the term has not made it to m_active_terms } - // deregister the term that has been activated - for (const auto & p: m_t->ext_coeffs()) { - auto it = m_s.m_columns_to_terms.find(p.var()); - SASSERT(it != m_s.m_columns_to_terms.end()); - it->second.erase(m_t->j()); - if (it->second.size() == 0) { - m_s.m_columns_to_terms.erase(it); - } - } - SASSERT(std::find(m_s.m_added_terms.begin(), m_s.m_added_terms.end(), m_t) == m_s.m_added_terms.end()); - SASSERT(contains(m_s.m_active_terms, m_t)); - m_s.m_active_terms.erase(m_t); - std::cout << ++ lp_settings::ddd << std::endl; - TRACE("dioph_eq", - tout << "the deleted term column in m_l_matrix" << std::endl; - for (auto p: m_s.m_l_matrix.column(m_t->j())) { - tout << "p.coeff():" << p.coeff()<< ", row " << p.var() << std::endl; - } - tout << "m_l_matrix has " << m_s.m_l_matrix.column_count() << " columns"<< std::endl; - tout << "and " << m_s.m_l_matrix.row_count() << " rows" << std::endl; - m_s.print_lar_term_L(*m_t, tout); tout << "; m_t->j()=" << m_t->j() << std::endl; - ); - m_s.shrink_L_to_size(m_row_count, m_l_column_count, m_e_column_count, m_t); + return; } - }; - - + // deregister the term that has been activated + for (const auto & p: t->ext_coeffs()) { + TRACE("dio_reg", tout << "derigister p.var():" << p.var() <<"->" << t->j() << std::endl; ); + auto it = m_columns_to_terms.find(p.var()); + SASSERT(it != m_columns_to_terms.end()); + it->second.erase(t->j()); + if (it->second.size() == 0) { + m_columns_to_terms.erase(it); + } + } + SASSERT(std::find(m_added_terms.begin(), m_added_terms.end(), t) == m_added_terms.end()); + SASSERT(contains(m_active_terms, t)); + m_active_terms.erase(t); + TRACE("dioph_eq", + tout << "the deleted term column in m_l_matrix" << std::endl; + for (auto p: m_l_matrix.column(t->j())) { + tout << "p.coeff():" << p.coeff()<< ", row " << p.var() << std::endl; + } + tout << "m_l_matrix has " << m_l_matrix.column_count() << " columns"<< std::endl; + tout << "and " << m_l_matrix.row_count() << " rows" << std::endl; + print_lar_term_L(*t, tout); tout << "; t->j()=" << t->j() << std::endl; + ); + shrink_L_to_sizes(); + } + struct undo_fixed_column : public trail { imp& m_imp; unsigned m_j; // the column that has been added @@ -350,7 +333,6 @@ namespace lp { m_imp.add_changed_column(m_j); } }; - void remove_last_entry() { unsigned ei = m_entries.size() - 1; @@ -360,100 +342,118 @@ namespace lp { } else { remove_entry_index(this->m_s, ei); } - } - void shrink_L_to_size(unsigned row_count, unsigned l_column_count, unsigned e_column_count, const lar_term* t) { - while (row_count < m_l_matrix.row_count()) { - shrink_L_one_size(t); - } - SASSERT(entries_are_ok()); - SASSERT(l_column_count == m_l_matrix.column_count() && e_column_count == m_e_matrix.column_count()); + m_entries.pop_back(); } - - void remove_last_row_column_in_matrix(static_matrix& m) { - auto & last_row = m.m_rows.back(); - for (unsigned k = last_row.size(); k-- > 0;) { - m.remove_element(last_row, last_row[k]); - } - SASSERT(m.m_columns.back().size() == 0); - m.m_rows.pop_back(); - m.m_columns.pop_back(); - } - - void shrink_L_one_size(const lar_term*t) { - unsigned j = m_l_matrix.column_count() - 1; // the last column - unsigned i = m_l_matrix.row_count() - 1; // the last row - auto last_col = m_l_matrix.column(j); - TRACE("dioph_eq", tout << "the last entry: i:" << i <j()="<< t->j() << " is mapped to " << lar_solver_to_local(t->j()) << std::endl; - ); - SASSERT(t->j() == j); - - bool row_i_found = false; - // prepare for pivoting - std_vector entries_to_change; - - for (const auto p : last_col ) { - entries_to_change.push_back(p.var()); ; - if (p.var() == i) { - row_i_found = true; - } - } - SASSERT(row_i_found); - - for (unsigned ei : entries_to_change) - move_entry_from_s_to_f(ei); - - // find the coefficient before the term variable, k, in the last row of m_e_matrix - unsigned k = lar_solver_to_local(j); - auto &last_e_row = m_e_matrix.m_rows.back(); + void eliminate_last_term_column() { + unsigned j = m_l_matrix.column_count() - 1; + make_sure_j_is_in_the_last_row_of_l_matrix(); + const auto &last_e_row = m_l_matrix.m_rows.back(); mpq alpha; for (const auto p: last_e_row) { - if (p.var() == k) { + if (p.var() == j) { alpha = p.coeff(); break; } } - this->m_e_matrix.divide_row(i, alpha); - this->m_l_matrix.divide_row(i, alpha); - this->m_entries[i].m_c /= alpha; - SASSERT(entry_invariant(i)); - eliminate_var_in_f(i, k, 1); - TRACE("dioph_eq", tout << "last row before removal\n" ; - print_entry(m_entries.size() -1, tout) << std::endl); - remove_last_row_column_in_matrix(m_l_matrix); - remove_last_row_column_in_matrix(m_e_matrix); - m_entries.pop_back(); - remove_entry_index(m_f, i); - remove_entry_index(m_s, i); + unsigned last_row_index= m_l_matrix.row_count() - 1; + m_l_matrix.divide_row(last_row_index, alpha); // divide the last row by alpha + std_vector rows_to_change; + + auto &column = m_l_matrix.m_columns[j]; + int pivot_col_cell_index = -1; + for (unsigned k = 0; k < column.size(); k++) { + if (column[k].var() == last_row_index) { + pivot_col_cell_index = k; + break; + } + } + SASSERT(pivot_col_cell_index >= 0); + + if (pivot_col_cell_index != 0) { + SASSERT(column.size() > 1); + // swap the pivot column cell with the head cell + auto c = column[0]; + column[0] = column[pivot_col_cell_index]; + column[pivot_col_cell_index] = c; + + m_l_matrix.m_rows[last_row_index][column[0].offset()].offset() = 0; + m_l_matrix.m_rows[c.var()][c.offset()].offset() = pivot_col_cell_index; + } + while (column.size() > 1) { + auto & c = column.back(); + SASSERT(c.var() != last_row_index); + m_l_matrix.pivot_row_to_row_given_cell(last_row_index, c, j); + rows_to_change.push_back(c.var()); + } + + for (unsigned i : rows_to_change) { + recalculate_entry(i); + } + } + + void make_sure_j_is_in_the_last_row_of_l_matrix() { + unsigned j = m_l_matrix.column_count() - 1; + const auto &last_e_row = m_l_matrix.m_rows.back(); + mpq alpha; + for (const auto p: last_e_row) { + if (p.var() == j) { + return; + } + } + SASSERT(m_l_matrix.m_columns.back().size()); + unsigned i = m_l_matrix.m_columns[j][0].var(); + m_l_matrix.add_rows(mpq(1), i, m_l_matrix.row_count() - 1); + } + + void shrink_L_to_sizes() { + unsigned i = m_l_matrix.row_count() - 1; + eliminate_last_term_column(); + remove_last_row_in_matrix(m_l_matrix); + remove_last_row_in_matrix(m_e_matrix); + while(m_l_matrix.column_count() && m_l_matrix.m_columns.back().size() == 0) { + m_l_matrix.m_columns.pop_back(); + } + while(m_e_matrix.column_count() && m_e_matrix.m_columns.back().size() == 0) { + m_e_matrix.m_columns.pop_back(); + } + m_var_register.shrink(m_e_matrix.column_count()); + auto it = std::find_if(m_fresh_definitions.begin(), m_fresh_definitions.end(), [i](auto const& fe) { return fe.m_origin == i; }); - if (it != m_fresh_definitions.end()) - m_fresh_definitions.erase(it); + if (it != m_fresh_definitions.end()) + *it = fresh_definition(-1, -1); for (unsigned k = 0; k < m_k2s.size() ; k++) { if (m_k2s[k] == i) { m_k2s[k] = -1; + break; } } - - m_var_register.shrink(m_e_matrix.column_count()); - } + remove_entry_index(m_f, i); + remove_entry_index(m_s, i); + m_entries.pop_back(); + } + + + void remove_last_row_in_matrix(static_matrix& m) { + auto & last_row = m.m_rows.back(); + for (unsigned k = last_row.size(); k-- > 0;) { + m.remove_element(last_row, last_row[k]); + } + m.m_rows.pop_back(); + } void remove_entry_index(std::list & l, unsigned ei) { auto it = std::find(l.begin(), l.end(), ei); if (it != l.end()) l.erase(it); } - - - + std::unordered_set m_changed_columns; // m_column_to_terms[j] is the set of all k such lra.get_term(k) depends on j std::unordered_map> m_columns_to_terms; - friend undo_fixed_column; void add_changed_column(unsigned j) { TRACE("dioph_eq", lra.print_column_info(j, tout);); m_changed_columns.insert(j); @@ -463,7 +463,7 @@ namespace lp { std_vector m_branch_stats; std_vector m_branch_stack; std_vector m_explanation_of_branches; - void add_term_delegate(const lar_term* t) { + void add_term_callback(const lar_term* t) { unsigned j = t->j(); TRACE("dioph_eq", tout << "term column t->j():" << j << std::endl; lra.print_term(*t, tout) << std::endl; ); if (!lra.column_is_int(j)) { @@ -477,12 +477,10 @@ namespace lp { 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); + m_added_terms.push_back(t); } - void update_column_bound_delegate(unsigned j) { + void update_column_bound_callback(unsigned j) { if (!lra.column_is_int(j) || !lra.column_is_fixed(j)) return; m_changed_columns.insert(j); @@ -492,8 +490,9 @@ namespace lp { public: imp(int_solver& lia, lar_solver& lra) : lia(lia), lra(lra) { - lra.m_add_term_callback=[this](const lar_term*t){add_term_delegate(t);}; - lra.m_update_column_bound_callback = [this](unsigned j){update_column_bound_delegate(j);}; + lra.m_add_term_callback=[this](const lar_term*t){add_term_callback(t);}; + lra.m_remove_term_callback = [this](const lar_term*t){remove_term_callback(t);}; + lra.m_update_column_bound_callback = [this](unsigned j){update_column_bound_callback(j);}; } term_o get_term_from_entry(unsigned i) const { term_o t; @@ -517,6 +516,8 @@ namespace lp { TRACE("dioph_eq", tout << "register term:"; lra.print_term(t, tout); tout << ", t.j()=" << t.j() << std::endl;); for (const auto &p: t.ext_coeffs()) { auto it = m_columns_to_terms.find(p.var()); + TRACE("dio_reg", tout << "register p.var():" << p.var() <<"->" << t.j() << std::endl; ); + if (it != m_columns_to_terms.end()) { it->second.insert(t.j()); } @@ -528,7 +529,7 @@ namespace lp { } } // the term has form sum(a_i*x_i) - t.j() = 0, - void fill_entry(const lar_term& t) { + void fill_entry(const lar_term& t) { TRACE("dioph_eq", print_lar_term_L(t, tout) << std::endl;); entry te = { mpq(0), entry_status::F}; unsigned entry_index = (unsigned) m_entries.size(); @@ -585,10 +586,19 @@ namespace lp { c = mpq(0); open_l_term_to_work_vector(ei, c); clear_e_row(ei); + mpq denom(1); for (const auto & p: m_indexed_work_vector) { unsigned lj = add_var(p.var()); m_e_matrix.add_columns_up_to(lj); m_e_matrix.add_new_element(ei, lj, p.coeff()); + if (!denominator(p.coeff()).is_one()) { + denom = lcm(denom, denominator(p.coeff())); + } + } + if (!denom.is_one()) { + c *= denom; + m_l_matrix.multiply_row(ei, denom); + m_e_matrix.multiply_row(ei, denom); } SASSERT(entry_invariant(ei)); } @@ -641,17 +651,92 @@ namespace lp { } tout << std::endl; ); - if (fresh_entries_to_remove.size()) { - NOT_IMPLEMENTED_YET(); + while (fresh_entries_to_remove.size()) { +// enable_trace("d_once"); + unsigned xt = fresh_entries_to_remove.back(); + fresh_entries_to_remove.pop_back(); + const fresh_definition & fd = m_fresh_definitions[xt]; + TRACE("d_once", print_entry(fd.m_ei, tout) << std::endl; tout << "xt:" << xt << std::endl;); + unsigned last_ei = m_entries.size() - 1; + if (fd.m_ei != last_ei) { // not the last entry + transpose_entries(fd.m_ei, m_entries.size() -1 ); + // we are not going to recalculate fd.m_ei + // but we might need to recalculate last_ei, which becomes fd.m_ei + if (contains(entries_to_recalculate, last_ei )) { + entries_to_recalculate.erase(last_ei), + entries_to_recalculate.insert(fd.m_ei); + } + } + for (const auto p: m_e_matrix.m_columns[xt]) { + entries_to_recalculate.insert(p.var()); + } + + m_fresh_definitions[xt] = fresh_definition(-1,-1); + remove_last_entry(); + remove_last_row_in_matrix(m_l_matrix); + remove_last_row_in_matrix(m_e_matrix); } for(unsigned k : entries_to_recalculate) { + if (k >= m_entries.size()) + continue;; recalculate_entry(k); - } - move_recalculated_to_F(entries_to_recalculate); + move_entry_from_s_to_f(k); + if (m_e_matrix.m_columns.back().size() == 0) { + m_e_matrix.m_columns.pop_back(); + m_var_register.shrink(m_e_matrix.column_count()); + } + if (m_l_matrix.m_columns.back().size() == 0) { + m_l_matrix.m_columns.pop_back(); + } + } m_changed_columns.clear(); } + void transpose_entries(unsigned i, unsigned k) { + SASSERT(i != k); + m_l_matrix.transpose_rows(i, k); + m_e_matrix.transpose_rows(i, k); + remove_entry_from_lists(i); + remove_entry_from_lists(k); + std::swap(m_entries[i], m_entries[k]); + add_entry_to_lists(i); + add_entry_to_lists(k); + // transpose fresh definitions + for (auto & fd: m_fresh_definitions) { + if (fd.m_ei == i) + fd.m_ei = k; + else if (fd.m_ei == k) + fd.m_ei = i; + } + //transpose m_k2s + for (unsigned& t : m_k2s) { + if (t == i) + t = k; + else if (t == k) + t = i; + } + } + + void remove_entry_from_lists(unsigned ei) { + auto status = m_entries[ei].m_entry_status; + if (status == entry_status::F) { + m_f.remove(ei); + } else { + m_s.remove(ei); + } + } + + void add_entry_to_lists(unsigned ei) { + auto status = m_entries[ei].m_entry_status; + if (status == entry_status::F) { + m_f.push_back(ei); + } else { + m_s.push_back(ei); + } + } + + void move_recalculated_to_F(const std::unordered_set &entries_to_recalculate) { for (auto it = this->m_s.begin(); it != m_s.end(); ) { if (contains(entries_to_recalculate, *it)) { @@ -1563,10 +1648,18 @@ namespace lp { std::ostream& print_e_row(unsigned i, std::ostream& out) { return print_term_o(get_term_from_entry(i), out); } + bool j_sign_is_correct(unsigned ei, unsigned j, int j_sign) { + const auto& row = m_e_matrix.m_rows[ei]; + auto it = std::find_if (row.begin(), row.end(), [j](const auto& p) {return p.var() == j;} ); + if (it == row.end()) return false; + return it->coeff() == mpq(1)&& j_sign == 1 || + it->coeff() == mpq(-1) && j_sign == -1; + } // j is the variable to eliminate, it appears in row ei of m_e_matrix with // a coefficient equal to j_sign which is +-1 void eliminate_var_in_f(unsigned ei, unsigned j, int j_sign) { const auto & e = m_entries[ei]; + SASSERT(j_sign_is_correct(ei, j, j_sign)); TRACE("dioph_eq", tout << "eliminate var:" << j << " by using:"; print_entry(ei, tout) << std::endl;); auto& column = m_e_matrix.m_columns[j]; @@ -1628,6 +1721,12 @@ namespace lp { } bool entry_invariant(unsigned ei) const { + for (const auto &p: m_e_matrix.m_rows[ei]) { + if (!p.coeff().is_int()) { + return false; + } + } + bool ret = term_to_lar_solver(remove_fresh_vars(get_term_from_entry(ei))) == fix_vars(open_ml(m_l_matrix.m_rows[ei])); @@ -1711,6 +1810,7 @@ namespace lp { // it clears the row, and puts the term in the work vector void move_row_to_work_vector(unsigned ei) { + m_indexed_work_vector.clear(); m_indexed_work_vector.resize(m_e_matrix.column_count()); auto& row = m_e_matrix.m_rows[ei]; for (const auto& cell : row) @@ -1720,8 +1820,8 @@ namespace lp { // k is the variable to substitute void fresh_var_step(unsigned h, unsigned k, const mpq& ahk) { - move_row_to_work_vector( - h); // it clears the row, and puts the term in the work vector + move_row_to_work_vector(h); // it clears the row, and puts the term in the work vector + // step 7 from the paper // xt is the fresh variable unsigned xt = add_var(UINT_MAX); @@ -1756,6 +1856,8 @@ namespace lp { m_e_matrix.add_new_element(fresh_row, i, q); } + m_l_matrix.add_row(); + m_k2s.resize(k + 1, -1); m_k2s[k] = fresh_row; fresh_definition fd(-1, -1); diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index 0c96a673e..a70b54c9d 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -1558,8 +1558,11 @@ namespace lp { if (col.term() != nullptr) { if (s.m_need_register_terms) s.deregister_normalized_term(*col.term()); + if (s.m_remove_term_callback) + s.m_remove_term_callback(col.term()); delete col.term(); - s.m_terms.pop_back(); + s.m_terms.pop_back(); + } s.remove_last_column_from_tableau(); s.m_columns.pop_back(); diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h index b3f697e61..f3a4daa6f 100644 --- a/src/math/lp/lar_solver.h +++ b/src/math/lp/lar_solver.h @@ -408,7 +408,8 @@ public: public: std::function m_add_term_callback; - std::function m_update_column_bound_callback; + std::function m_remove_term_callback; + std::function m_update_column_bound_callback; bool external_is_used(unsigned) const; void pop(unsigned k); unsigned num_scopes() const { return m_trail.get_num_scopes(); }