From b52e79b648c6e243109b355c6dee3283b9d54cc6 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Wed, 17 Apr 2019 10:41:14 -0700 Subject: [PATCH] emons branch Signed-off-by: Lev Nachmanson --- src/test/emonomials.cpp | 36 ++++ src/util/lp/CMakeLists.txt | 1 + src/util/lp/emonomials.cpp | 286 +++++++++++++++++++++++++++++++ src/util/lp/emonomials.h | 246 ++++++++++++++++++++++++++ src/util/lp/incremental_vector.h | 76 ++++++++ src/util/lp/monomial.h | 14 +- src/util/lp/nla_core.cpp | 65 ++----- src/util/lp/nla_core.h | 5 +- src/util/lp/nla_defs.h | 103 +++++++++++ src/util/lp/var_eqs.cpp | 282 +++++++++++++++--------------- src/util/lp/var_eqs.h | 77 +++++---- 11 files changed, 955 insertions(+), 236 deletions(-) create mode 100644 src/test/emonomials.cpp create mode 100644 src/util/lp/emonomials.cpp create mode 100644 src/util/lp/emonomials.h create mode 100644 src/util/lp/incremental_vector.h create mode 100644 src/util/lp/nla_defs.h diff --git a/src/test/emonomials.cpp b/src/test/emonomials.cpp new file mode 100644 index 000000000..f5d7e6e94 --- /dev/null +++ b/src/test/emonomials.cpp @@ -0,0 +1,36 @@ +#include "util/lp/emonomials.h" + +static void display_mons(nla::emonomials const& em) { + std::cout << em << "\n"; + for (auto const& m1 : em) { + std::cout << "factors: " << m1 << ": \n"; + for (auto const& m2 : em.get_factors_of(m1)) { + std::cout << m2 << "\n"; + } + } +} + +void tst_emonomials() { + nla::var_eqs ve; + nla::emonomials em(ve); + + unsigned x = 1; + unsigned y = 2; + unsigned z = 3; + unsigned u = 4; + unsigned v1 = 5; + unsigned v2 = 6; + unsigned v3 = 7; + em.add(v1, x, x); + em.add(v2, x, z); + em.add(v3, x, y, x); + display_mons(em); + em.push(); + ve.merge_plus(x, y, nla::eq_justification({})); + display_mons(em); + em.push(); + ve.merge_plus(x, z, nla::eq_justification({})); + display_mons(em); + em.pop(1); + display_mons(em); +} diff --git a/src/util/lp/CMakeLists.txt b/src/util/lp/CMakeLists.txt index 16cd21af6..bda497378 100644 --- a/src/util/lp/CMakeLists.txt +++ b/src/util/lp/CMakeLists.txt @@ -6,6 +6,7 @@ z3_add_component(lp core_solver_pretty_printer.cpp dense_matrix.cpp eta_matrix.cpp + emonomials.cpp factorization.cpp factorization_factory_imp.cpp gomory.cpp diff --git a/src/util/lp/emonomials.cpp b/src/util/lp/emonomials.cpp new file mode 100644 index 000000000..3b4532be6 --- /dev/null +++ b/src/util/lp/emonomials.cpp @@ -0,0 +1,286 @@ +/*++ + Copyright (c) 2019 Microsoft Corporation + + Module Name: + + emonomials.cpp + + Abstract: + + table that associate monomials to congruence class representatives modulo a union find structure. + + Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + + Revision History: + + to replace rooted_mons.h and rooted_mon, rooted_mon_tabled + + --*/ + +#include "util/lp/emonomials.h" +#include "util/lp/nla_defs.h" + +namespace nla { + + void emonomials::inc_canonical() { + ++m_canonical; + if (m_canonical == 0) { + for (auto& svt : m_canonized) { + svt.m_canonical = 0; + } + ++m_canonical; + } + } + + void emonomials::inc_visited() const { + ++m_visited; + if (m_visited == 0) { + for (auto& svt : m_canonized) { + svt.m_visited = 0; + } + ++m_visited; + } + } + + void emonomials::push() { + m_lim.push_back(m_monomials.size()); + m_region.push_scope(); + m_ve.push(); + } + + void emonomials::pop(unsigned n) { + m_ve.pop(n); + unsigned old_sz = m_lim[m_lim.size() - n]; + for (unsigned i = m_monomials.size(); i-- > old_sz; ) { + monomial const& m = m_monomials[i]; + m_var2index[m.var()] = UINT_MAX; + lpvar last_var = UINT_MAX; + for (lpvar v : m.vars()) { + if (v != last_var) { + remove_var2monomials(v, i); + last_var = v; + } + } + } + m_monomials.shrink(old_sz); + m_canonized.shrink(old_sz); + m_region.pop_scope(n); + m_lim.shrink(m_lim.size() - n); + } + + void emonomials::remove_var2monomials(lpvar v, unsigned mIndex) { + cell*& cur_head = m_use_lists[v].m_head; + cell*& cur_tail = m_use_lists[v].m_tail; + cell* old_head = cur_head->m_next; + if (old_head == cur_head) { + cur_head = nullptr; + cur_tail = nullptr; + } + else { + cur_head = old_head; + cur_tail->m_next = old_head; + } + } + + void emonomials::insert_var2monomials(lpvar v, unsigned mIndex) { + m_use_lists.reserve(v + 1); + cell*& cur_head = m_use_lists[v].m_head; + cell*& cur_tail = m_use_lists[v].m_tail; + cell* new_head = new (m_region) cell(mIndex, cur_head); + cur_head = new_head; + if (!cur_tail) cur_tail = new_head; + cur_tail->m_next = new_head; + } + + void emonomials::merge_var2monomials(lpvar root, lpvar other) { + if (root == other) return; + m_use_lists.reserve(std::max(root, other) + 1); + cell*& root_head = m_use_lists[root].m_head; + cell*& root_tail = m_use_lists[root].m_tail; + cell* other_head = m_use_lists[other].m_head; + cell* other_tail = m_use_lists[other].m_tail; + if (root_head == nullptr) { + root_head = other_head; + root_tail = other_tail; + } + else if (other_head) { + // other_head -> other_tail -> root_head --> root_tail -> other_head. + root_tail->m_next = other_head; + other_tail->m_next = root_head; + root_head = other_head; + } + else { + // other_head = other_tail = nullptr + } + } + + void emonomials::unmerge_var2monomials(lpvar root, lpvar other) { + if (root == other) return; + cell*& root_head = m_use_lists[root].m_head; + cell*& root_tail = m_use_lists[root].m_tail; + cell* other_head = m_use_lists[other].m_head; + cell* other_tail = m_use_lists[other].m_tail; + if (other_head == nullptr) { + // no-op + } + else if (root_tail == other_tail) { + root_head = nullptr; + root_tail = nullptr; + } + else { + root_head = other_tail->m_next; + root_tail->m_next = root_head; + other_tail->m_next = other_head; + } + } + + emonomials::cell* emonomials::head(lpvar v) const { + v = m_ve.find(v).var(); + m_use_lists.reserve(v + 1); + return m_use_lists[v].m_head; + } + + void emonomials::set_visited(monomial const& m) const { + m_canonized[m_var2index[m.var()]].m_visited = m_visited; + } + + bool emonomials::is_visited(monomial const& m) const { + return m_visited == m_canonized[m_var2index[m.var()]].m_visited; + } + + + /** + \brief insert a new monomial. + + Assume that the variables are canonical, that is, not equal in current + context so another variable. To support equal in current context we + could track sign information in monomials, or ensure that insert_var2monomials + works on the equivalence class roots. + */ + void emonomials::add(lpvar v, unsigned sz, lpvar const* vs) { + unsigned idx = m_monomials.size(); + m_monomials.push_back(monomial(v, sz, vs)); + m_canonized.push_back(signed_vars_ts()); + lpvar last_var = UINT_MAX; + for (unsigned i = 0; i < sz; ++i) { + lpvar w = vs[i]; + SASSERT(m_ve.is_root(w)); + if (w != last_var) { + insert_var2monomials(w, idx); + last_var = w; + } + } + SASSERT(m_ve.is_root(v)); + m_var2index.setx(v, idx, UINT_MAX); + } + + signed_vars const& emonomials::canonize(monomial const& mon) const { + unsigned index = m_var2index[mon.var()]; + if (m_canonized[index].m_canonical != m_canonical) { + m_vars.reset(); + bool sign = false; + for (lpvar v : mon) { + signed_var sv = m_ve.find(v); + m_vars.push_back(sv.var()); + sign ^= sv.sign(); + } + m_canonized[index].set_vars(m_vars.size(), m_vars.c_ptr()); + m_canonized[index].set_sign(sign); + m_canonized[index].set_var(mon.var()); + m_canonized[index].m_canonical = m_canonical; + } + return m_canonized[index]; + } + + bool emonomials::canonize_divides(monomial const& m1, monomial const& m2) const { + if (m1.size() > m2.size()) return false; + signed_vars const& s1 = canonize(m1); + signed_vars const& s2 = canonize(m2); + unsigned sz1 = s1.size(), sz2 = s2.size(); + unsigned i = 0, j = 0; + while (true) { + if (i == sz1) { + return true; + } + else if (j == sz2) { + return false; + } + else if (s1[i] == s2[j]) { + ++i; ++j; + } + else if (s1[i] < s2[j]) { + return false; + } + else { + ++j; + } + } + } + + void emonomials::explain_canonized(monomial const& m, lp::explanation& exp) { + for (lpvar v : m) { + signed_var w = m_ve.find(v); + m_ve.explain(signed_var(v, false), w, exp); + } + } + + // yes, assume that monomials are non-empty. + emonomials::pf_iterator::pf_iterator(emonomials const& m, monomial const& mon, bool at_end): + m(m), m_mon(mon), m_it(iterator(m, m.head(mon[0]), at_end)), m_end(iterator(m, m.head(mon[0]), true)) { + fast_forward(); + } + + void emonomials::pf_iterator::fast_forward() { + for (; m_it != m_end; ++m_it) { + if (m_mon.var() != (*m_it).var() && m.canonize_divides(m_mon, *m_it) && !m.is_visited(*m_it)) { + m.set_visited(*m_it); + break; + } + } + } + + void emonomials::merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) { + if (!r2.sign() && m_ve.find(~r2) != m_ve.find(r1)) { + merge_var2monomials(r2.var(), r1.var()); + } + inc_canonical(); + } + + void emonomials::after_merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) { + // skip + } + + void emonomials::unmerge_eh(signed_var r2, signed_var r1) { + if (!r2.sign() && m_ve.find(~r2) != m_ve.find(r1)) { + unmerge_var2monomials(r2.var(), r1.var()); + } + inc_canonical(); + } + + std::ostream& emonomials::display(std::ostream& out) const { + out << "monomials\n"; + unsigned idx = 0; + for (auto const& m : m_monomials) { + out << (idx++) << ": " << m << "\n"; + } + out << "use lists\n"; + idx = 0; + for (auto const& ht : m_use_lists) { + cell* c = ht.m_head; + if (c) { + out << "v" << idx << ": "; + do { + out << c->m_index << " "; + c = c->m_next; + } + while (c != ht.m_head); + out << "\n"; + } + ++idx; + } + return out; + } + +} diff --git a/src/util/lp/emonomials.h b/src/util/lp/emonomials.h new file mode 100644 index 000000000..8f1135104 --- /dev/null +++ b/src/util/lp/emonomials.h @@ -0,0 +1,246 @@ +/*++ + Copyright (c) 2019 Microsoft Corporation + + Module Name: + + emonomials.h + + Abstract: + + table that associate monomials to congruence class representatives modulo a union find structure. + + Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + + Revision History: + + to replace rooted_mons.h and rooted_mon, rooted_mon_tabled + + --*/ + +#pragma once +#include "util/lp/lp_utils.h" +#include "util/lp/var_eqs.h" +#include "util/lp/monomial.h" +#include "util/region.h" + +namespace nla { + + /** + \brief class used to summarize the coefficients to a monomial after + canonization with respect to current equalities. + */ + class signed_vars { + lpvar m_var; // variable representing original monomial + svector m_vars; + bool m_sign; + public: + signed_vars() : m_sign(false) {} + lpvar var() const { return m_var; } + svector const& vars() const { return m_vars; } + unsigned size() const { return m_vars.size(); } + lpvar operator[](unsigned i) const { return m_vars[i]; } + bool sign() const { return m_sign; } + void set_sign(bool s) { m_sign = s; } + void set_var(lpvar v) { m_var = v; } + void set_vars(unsigned n, lpvar const* vars) { + m_vars.reset(); + m_vars.append(n, vars); + std::sort(m_vars.begin(), m_vars.end()); + } + std::ostream& display(std::ostream& out) const { + out << "v" << var() << " := "; + if (sign()) out << "- "; + for (lpvar v : vars()) out << "v" << v << " "; + return out; + } + }; + + inline std::ostream& operator<<(std::ostream& out, signed_vars const& m) { return m.display(out); } + + + class emonomials : public var_eqs_merge_handler { + /** + \brief private fields used by emonomials for maintaining state of canonized monomials. + */ + class signed_vars_ts : public signed_vars { + public: + signed_vars_ts(): m_canonical(0), m_visited(0) {} + unsigned m_canonical; + unsigned m_visited; + }; + + /** + \brief singly-lined cyclic list of monomial indices where variable occurs. + Each variable points to the head and tail of the cyclic list. + Initially, head and tail are nullptr. + New elements are inserted in the beginning of the list. + Two lists are merged when equivalence class representatives are merged, + and the merge is undone when the representative variables are unmerged. + */ + struct cell { + cell(unsigned mIndex, cell* c): m_next(c), m_index(mIndex) {} + cell* m_next; + unsigned m_index; + }; + struct head_tail { + head_tail(): m_head(nullptr), m_tail(nullptr) {} + cell* m_head; + cell* m_tail; + }; + + var_eqs& m_ve; + vector m_monomials; // set of monomials + unsigned_vector m_lim; // backtracking point + unsigned_vector m_var2index; // var_mIndex -> mIndex + unsigned m_canonical; // timestamp of last merge + mutable unsigned m_visited; // timestamp of visited monomials during pf_iterator + region m_region; // region for allocating linked lists + mutable vector m_canonized; // canonized versions of signed variables + mutable svector m_vars; // temporary vector of variables + mutable svector m_use_lists; // use list of monomials where variables occur. + + void inc_canonical(); + void inc_visited() const; + + void remove_var2monomials(lpvar v, unsigned mIndex); + void insert_var2monomials(lpvar v, unsigned mIndex); + void merge_var2monomials(lpvar root, lpvar other); + void unmerge_var2monomials(lpvar root, lpvar other); + + cell* head(lpvar v) const; + void set_visited(monomial const& m) const; + bool is_visited(monomial const& m) const; + + public: + /** + \brief emonomials builds on top of var_eqs. + push and pop on emonomials calls push/pop on var_eqs, so no + other calls to push/pop to the var_eqs should take place. + */ + emonomials(var_eqs& ve): m_ve(ve), m_canonical(1), m_visited(0) { m_ve.set_merge_handler(this); } + + /** + \brief push/pop scopes. + The life-time of a merge is local within a scope. + */ + void push(); + + void pop(unsigned n); + + /** + \brief create a monomial from an equality v := vs + */ + void add(lpvar v, unsigned sz, lpvar const* vs); + void add(lpvar v, svector const& vs) { add(v, vs.size(), vs.c_ptr()); } + void add(lpvar v, lpvar x, lpvar y) { lpvar vs[2] = { x, y }; add(v, 2, vs); } + void add(lpvar v, lpvar x, lpvar y, lpvar z) { lpvar vs[3] = { x, y, z }; add(v, 3, vs); } + + /** + \brief retrieve monomial corresponding to variable v from definition v := vs + */ + monomial const* var2monomial(lpvar v) const { unsigned idx = m_var2index.get(v, UINT_MAX); return idx == UINT_MAX ? nullptr : &m_monomials[idx]; } + + /** + \brief obtain a canonized signed monomial + corresponding to current equivalence class. + */ + signed_vars const& canonize(monomial const& m) const; + + /** + \brief determine if m1 divides m2 over the canonization obtained from merged variables. + */ + bool canonize_divides(monomial const& m1, monomial const& m2) const; + + /** + \brief produce explanation for monomial canonization. + */ + void explain_canonized(monomial const& m, lp::explanation& exp); + + /** + \brief iterator over monomials that are declared. + */ + vector::const_iterator begin() const { return m_monomials.begin(); } + vector::const_iterator end() const { return m_monomials.end(); } + + /** + \brief iterators over monomials where an equivalent variable is used + */ + class iterator { + emonomials const& m; + cell* m_cell; + bool m_touched; + public: + iterator(emonomials const& m, cell* c, bool at_end): m(m), m_cell(c), m_touched(at_end || c == nullptr) {} + monomial const& operator*() { return m.m_monomials[m_cell->m_index]; } + iterator& operator++() { m_touched = true; m_cell = m_cell->m_next; return *this; } + iterator operator++(int) { iterator tmp = *this; ++*this; return tmp; } + bool operator==(iterator const& other) const { return m_cell == other.m_cell && m_touched == other.m_touched; } + bool operator!=(iterator const& other) const { return m_cell != other.m_cell || m_touched != other.m_touched; } + }; + + class use_list { + emonomials const& m; + lpvar m_var; + cell* head() { return m.head(m_var); } + public: + use_list(emonomials const& m, lpvar v): m(m), m_var(v) {} + iterator begin() { return iterator(m, head(), false); } + iterator end() { return iterator(m, head(), true); } + }; + + use_list get_use_list(lpvar v) const { return use_list(*this, v); } + + + /** + \brief retrieve monomials m' where m is a proper factor of modulo current equalities. + */ + class pf_iterator { + emonomials const& m; + monomial const& m_mon; // canonized monomial + iterator m_it; // iterator over the first variable occurs list, ++ filters out elements that are not factors. + iterator m_end; + + void fast_forward(); + public: + pf_iterator(emonomials const& m, monomial const& mon, bool at_end); + monomial const& operator*() { return *m_it; } + pf_iterator& operator++() { ++m_it; fast_forward(); return *this; } + pf_iterator operator++(int) { pf_iterator tmp = *this; ++*this; return tmp; } + bool operator==(pf_iterator const& other) const { return m_it == other.m_it; } + bool operator!=(pf_iterator const& other) const { return m_it != other.m_it; } + }; + + class factors_of { + emonomials const& m; + monomial const& mon; + public: + factors_of(emonomials const& m, monomial const& mon): m(m), mon(mon) {} + pf_iterator begin() { return pf_iterator(m, mon, false); } + pf_iterator end() { return pf_iterator(m, mon, true); } + }; + + factors_of get_factors_of(monomial const& m) const { inc_visited(); return factors_of(*this, m); } + + /** + \brief display state of emonomials + */ + std::ostream& display(std::ostream& out) const; + + /** + \brief + these are merge event handlers to interect the union-find handlers. + r2 becomes the new root. r2 is the root of v2, r1 is the old root of v1 + */ + void merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) override; + + void after_merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) override; + + void unmerge_eh(signed_var r2, signed_var r1) override; + + }; + + inline std::ostream& operator<<(std::ostream& out, emonomials const& m) { return m.display(out); } + +} diff --git a/src/util/lp/incremental_vector.h b/src/util/lp/incremental_vector.h new file mode 100644 index 000000000..55c76a1fa --- /dev/null +++ b/src/util/lp/incremental_vector.h @@ -0,0 +1,76 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ + +#pragma once +#include +#include +#include +#include "util/vector.h" +namespace lp { +template class incremental_vector { + vector m_stack_of_vector_sizes; + vector m_vector; +public: + const B & operator[](unsigned a) const { + return m_vector[a]; + } + + unsigned size() const { + return m_vector.size(); + } + + void push() { + m_stack_of_vector_sizes.push_back(m_vector.size()); + } + + void pop() { + pop(1); + } + + template + void pop_tail(vector & v, unsigned k) { + lp_assert(v.size() >= k); + v.shrink(v.size() - k); + } + + template + void resize(vector & v, unsigned new_size) { + v.resize(new_size); + } + + void pop(unsigned k) { + lp_assert(m_stack_of_vector_sizes.size() >= k); + lp_assert(k > 0); + m_vector.shrink(peek_size(k)); + unsigned new_st_size = m_stack_of_vector_sizes.size() - k; + m_stack_of_vector_sizes.shrink(k); + } + + void push_back(const B & b) { + m_vector.push_back(b); + } + + unsigned peek_size(unsigned k) const { + lp_assert(k > 0 && k <= m_stack_of_vector_sizes.size()); + return m_stack_of_vector_sizes[m_stack_of_vector_sizes.size() - k]; + } +}; +} + diff --git a/src/util/lp/monomial.h b/src/util/lp/monomial.h index d5c0bf24f..10f2e99e9 100644 --- a/src/util/lp/monomial.h +++ b/src/util/lp/monomial.h @@ -35,10 +35,22 @@ namespace nla { unsigned operator[](unsigned idx) const { return m_vs[idx]; } svector::const_iterator begin() const { return m_vs.begin(); } svector::const_iterator end() const { return m_vs.end(); } - const svector vars() const { return m_vs; } + const svector& vars() const { return m_vs; } bool empty() const { return m_vs.empty(); } + + std::ostream& display(std::ostream& out) const { + out << "v" << var() << " := "; + for (auto v : *this) { + out << "v" << v << " "; + } + return out; + } }; + inline std::ostream& operator<<(std::ostream& out, monomial const& m) { + return m.display(out); + } + typedef std::unordered_map variable_map_type; bool check_assignment(monomial const& m, variable_map_type & vars); diff --git a/src/util/lp/nla_core.cpp b/src/util/lp/nla_core.cpp index 0f862da28..27da6dd52 100644 --- a/src/util/lp/nla_core.cpp +++ b/src/util/lp/nla_core.cpp @@ -27,7 +27,8 @@ core::core(lp::lar_solver& s) : m_tangents(this), m_basics(this), m_order(this), - m_monotone(this) { + m_monotone(this), + m_emons(m_evars) { } bool core::compare_holds(const rational& ls, llc cmp, const rational& rs) const { @@ -1608,25 +1609,24 @@ void core::init_search() { register_monomials_in_tables(); } -#if 0 -void init_to_refine() { +void core::init_to_refine() { m_to_refine.clear(); for (auto const & m : m_emons) if (!check_monomial(m)) m_to_refine.push_back(m.var()); TRACE("nla_solver", - tout << m_to_refine.size() << " mons to refine:\n"); - for (unsigned v : m_to_refine) tout << m_emons.var2monomial(v) << "\n"; + tout << m_to_refine.size() << " mons to refine:\n"; + for (unsigned v : m_to_refine) tout << m_emons.var2monomial(v) << "\n";); } - -std::unordered_set collect_vars( const lemma& l) { + +std::unordered_set core::collect_vars(const lemma& l) const { std::unordered_set vars; auto insert_j = [&](lpvar j) { - vars.insert(j); - auto const* m = m_emons.var2monomial(j); - if (m) for (lpvar k : *m) vars.insert(k); - }; + vars.insert(j); + auto const* m = m_emons.var2monomial(j); + if (m) for (lpvar k : *m) vars.insert(k); + }; for (const auto& i : current_lemma().ineqs()) { for (const auto& p : i.term()) { @@ -1641,21 +1641,6 @@ std::unordered_set collect_vars( const lemma& l) { } return vars; } -#endif - -void core::init_to_refine() { - m_to_refine.clear(); - for (unsigned i = 0; i < m_monomials.size(); i++) { - if (!check_monomial(m_monomials[i])) - m_to_refine.push_back(i); - } - TRACE("nla_solver", - tout << m_to_refine.size() << " mons to refine:\n"; - for (unsigned i: m_to_refine) { - print_monomial_with_vars(m_monomials[i], tout); - } - ); -} bool core:: divide(const rooted_mon& bc, const factor& c, factor & b) const { svector c_vars = sorted_vars(c); @@ -1706,34 +1691,6 @@ void core::negate_factor_relation(const rational& a_sign, const factor& a, const mk_ineq(a_fs*a_sign, var(a), - b_fs*b_sign, var(b), cmp); } -std::unordered_set core::collect_vars(const lemma& l) const { - std::unordered_set vars; - for (const auto& i : current_lemma().ineqs()) { - for (const auto& p : i.term()) { - lpvar j = p.var(); - vars.insert(j); - auto it = m_var_to_its_monomial.find(j); - if (it != m_var_to_its_monomial.end()) { - for (lpvar k : m_monomials[it->second]) - vars.insert(k); - } - } - } - for (const auto& p : current_expl()) { - const auto& c = m_lar_solver.get_constraint(p.second); - for (const auto& r : c.coeffs()) { - lpvar j = r.second; - vars.insert(j); - auto it = m_var_to_its_monomial.find(j); - if (it != m_var_to_its_monomial.end()) { - for (lpvar k : m_monomials[it->second]) - vars.insert(k); - } - } - } - return vars; -} - std::ostream& core::print_lemma(std::ostream& out) const { print_specific_lemma(current_lemma(), out); return out; diff --git a/src/util/lp/nla_core.h b/src/util/lp/nla_core.h index d89dbab02..2272716a6 100644 --- a/src/util/lp/nla_core.h +++ b/src/util/lp/nla_core.h @@ -26,6 +26,7 @@ #include "util/lp/nla_basics_lemmas.h" #include "util/lp/nla_order_lemmas.h" #include "util/lp/nla_monotone_lemmas.h" +#include "util/lp/emonomials.h" namespace nla { template @@ -79,7 +80,6 @@ struct core { var_eqs m_evars; vector m_monomials; - rooted_mon_table m_rm_table; // this field is used for the push/pop operations unsigned_vector m_monomials_counts; lp::lar_solver& m_lar_solver; @@ -94,6 +94,7 @@ struct core { basics m_basics; order m_order; monotone m_monotone; + emonomials m_emons; // methods core(lp::lar_solver& s); @@ -115,7 +116,7 @@ struct core { rational vvr(const rooted_mon& rm) const { return vvr(m_monomials[rm.orig_index()].var()) * rm.orig_sign(); } - rational vvr(const factor& f) const { return f.is_var()? vvr(f.index()) : vvr(m_rm_table.rms()[f.index()]); } + rational vvr(const factor& f) const { return f.is_var()? vvr(f.index()) : vvr(m_rm_table.rms()[f.index()]); } lpvar var(const factor& f) const { return f.is_var()? f.index() : var(m_rm_table.rms()[f.index()]); } diff --git a/src/util/lp/nla_defs.h b/src/util/lp/nla_defs.h new file mode 100644 index 000000000..37c0085ed --- /dev/null +++ b/src/util/lp/nla_defs.h @@ -0,0 +1,103 @@ +/*++ + Copyright (c) 2017 Microsoft Corporation + + Module Name: + + + + Abstract: + + + + Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + + Revision History: + + + --*/ +#pragma once +#include "util/lp/lp_types.h" +#include "util/lp/column_info.h" +#include "util/lp/explanation.h" + +namespace nla { +typedef lp::constraint_index lpci; +typedef lp::lconstraint_kind llc; +typedef lp::constraint_index lpci; +typedef lp::explanation expl_set; +typedef lp::var_index lpvar; + +struct from_index_dummy{}; +class signed_var { + unsigned m_sv; +public: + // constructor, sign = true means minus + signed_var(lpvar v, bool sign): m_sv((v << 1) + (sign ? 1 : 0)) {} + // constructor + signed_var(unsigned sv, from_index_dummy): m_sv(sv) {} + bool sign() const { return 0 != (m_sv & 0x1); } + lpvar var() const { return m_sv >> 1; } + unsigned index() const { return m_sv; } + void neg() { m_sv = m_sv ^ 1; } + friend signed_var operator~(signed_var const& sv) { + return signed_var(sv.var(), !sv.sign()); + } + bool operator==(signed_var const& other) const { + return m_sv == other.m_sv; + } + bool operator!=(signed_var const& other) const { + return m_sv != other.m_sv; + } + rational rsign() const { return sign() ? rational::minus_one() : rational::one(); } + + std::ostream& display(std::ostream& out) const { + return out << (sign()?"-":"") << var(); + } +}; + +inline std::ostream& operator<<(std::ostream& out, signed_var const& sv) { return sv.display(out); } + +/* + * represents definition m_v = coeff* v1*v2*...*vn, + * where m_vs = [v1, v2, .., vn] + */ +class monomial_coeff { + svector m_vs; + rational m_coeff; +public: + monomial_coeff(const svector& vs, rational const& coeff): m_vs(vs), m_coeff(coeff) {} + rational const& coeff() const { return m_coeff; } + const svector & vars() const { return m_vs; } +}; +template bool has_zero(const T& product) { + for (const rational & t : product) { + if (t.is_zero()) + return true; + } + return false; +} +template +bool uniform_le(const T& a, const T& b, unsigned & strict_i) { + SASSERT(a.size() == b.size()); + strict_i = -1; + bool z_b = false; + + for (unsigned i = 0; i < a.size(); i++) { + if (a[i] > b[i]){ + return false; + } + + SASSERT(!a[i].is_neg()); + if (a[i] < b[i]){ + strict_i = i; + } else if (b[i].is_zero()) { + z_b = true; + } + } + if (z_b) {strict_i = -1;} + return true; +} + +} diff --git a/src/util/lp/var_eqs.cpp b/src/util/lp/var_eqs.cpp index 2dc2118e5..7fc711d25 100644 --- a/src/util/lp/var_eqs.cpp +++ b/src/util/lp/var_eqs.cpp @@ -23,164 +23,164 @@ namespace nla { - var_eqs::var_eqs(): m_uf(m_ufctx) {} +var_eqs::var_eqs(): m_merge_handler(nullptr), m_uf(*this), m_stack(*this) {} - void var_eqs::push() { - m_trail_lim.push_back(m_trail.size()); - m_ufctx.get_trail_stack().push_scope(); - } +void var_eqs::push() { + m_trail_lim.push_back(m_trail.size()); + get_trail_stack().push_scope(); +} - void var_eqs::pop(unsigned n) { - unsigned old_sz = m_trail_lim[m_trail_lim.size() - n]; - for (unsigned i = m_trail.size(); i-- > old_sz; ) { - auto const& sv = m_trail[i]; - m_eqs[sv.first.index()].pop_back(); - m_eqs[sv.second.index()].pop_back(); - m_eqs[(~sv.first).index()].pop_back(); - m_eqs[(~sv.second).index()].pop_back(); +void var_eqs::pop(unsigned n) { + unsigned old_sz = m_trail_lim[m_trail_lim.size() - n]; + for (unsigned i = m_trail.size(); i-- > old_sz; ) { + auto const& sv = m_trail[i]; + m_eqs[sv.first.index()].pop_back(); + m_eqs[sv.second.index()].pop_back(); + m_eqs[(~sv.first).index()].pop_back(); + m_eqs[(~sv.second).index()].pop_back(); + } + m_trail_lim.shrink(m_trail_lim.size() - n); + m_trail.shrink(old_sz); + get_trail_stack().pop_scope(n); +} + +void var_eqs::merge(signed_var v1, signed_var v2, eq_justification const& j) { + unsigned max_i = std::max(v1.index(), v2.index()) + 2; + m_eqs.reserve(max_i); + while (m_uf.get_num_vars() <= max_i) m_uf.mk_var(); + m_trail.push_back(std::make_pair(v1, v2)); + m_uf.merge(v1.index(), v2.index()); + m_uf.merge((~v1).index(), (~v2).index()); + m_eqs[v1.index()].push_back(justified_var(v2, j)); + m_eqs[v2.index()].push_back(justified_var(v1, j)); + m_eqs[(~v1).index()].push_back(justified_var(~v2, j)); + m_eqs[(~v2).index()].push_back(justified_var(~v1, j)); +} + +signed_var var_eqs::find(signed_var v) const { + if (v.index() >= m_uf.get_num_vars()) { + return v; + } + unsigned idx = m_uf.find(v.index()); + return signed_var(idx, from_index_dummy()); +} + +void var_eqs::explain_dfs(signed_var v1, signed_var v2, lp::explanation& e) const { + SASSERT(find(v1) == find(v2)); + if (v1 == v2) { + return; + } + m_todo.push_back(var_frame(v1, 0)); + m_jtrail.reset(); + m_marked.reserve(m_eqs.size(), false); + SASSERT(m_marked_trail.empty()); + m_marked[v1.index()] = true; + m_marked_trail.push_back(v1.index()); + while (true) { + SASSERT(!m_todo.empty()); + var_frame& f = m_todo.back(); + signed_var v = f.m_var; + if (v == v2) { + break; } - m_trail_lim.shrink(m_trail_lim.size() - n); - m_trail.shrink(old_sz); - m_ufctx.get_trail_stack().pop_scope(n); - } - - void var_eqs::merge(signed_var v1, signed_var v2, eq_justification const& j) { - unsigned max_i = std::max(v1.index(), v2.index()) + 2; - m_eqs.reserve(max_i); - while (m_uf.get_num_vars() <= max_i) m_uf.mk_var(); - m_trail.push_back(std::make_pair(v1, v2)); - m_uf.merge(v1.index(), v2.index()); - m_uf.merge((~v1).index(), (~v2).index()); - m_eqs[v1.index()].push_back(justified_var(v2, j)); - m_eqs[v2.index()].push_back(justified_var(v1, j)); - m_eqs[(~v1).index()].push_back(justified_var(~v2, j)); - m_eqs[(~v2).index()].push_back(justified_var(~v1, j)); - } - - signed_var var_eqs::find(signed_var v) const { - if (v.index() >= m_uf.get_num_vars()) { - return v; - } - unsigned idx = m_uf.find(v.index()); - return signed_var(idx, from_index()); - } - - void var_eqs::explain_dfs(signed_var v1, signed_var v2, lp::explanation& e) const { - SASSERT(find(v1) == find(v2)); - if (v1 == v2) { - return; - } - m_todo.push_back(var_frame(v1, 0)); - m_jtrail.reset(); - m_marked.reserve(m_eqs.size(), false); - SASSERT(m_marked_trail.empty()); - m_marked[v1.index()] = true; - m_marked_trail.push_back(v1.index()); - while (true) { - SASSERT(!m_todo.empty()); - var_frame& f = m_todo.back(); - signed_var v = f.m_var; - if (v == v2) { - break; - } - auto const& next = m_eqs[v.index()]; - bool seen_all = true; - unsigned sz = next.size(); - for (unsigned i = f.m_index; seen_all && i < sz; ++i) { - justified_var const& jv = next[i]; - signed_var v3 = jv.m_var; - if (!m_marked[v3.index()]) { - seen_all = false; - f.m_index = i + 1; - m_todo.push_back(var_frame(v3, 0)); - m_jtrail.push_back(jv.m_j); - m_marked_trail.push_back(v3.index()); - m_marked[v3.index()] = true; - } - } - if (seen_all) { - m_todo.pop_back(); - m_jtrail.pop_back(); + auto const& next = m_eqs[v.index()]; + bool seen_all = true; + unsigned sz = next.size(); + for (unsigned i = f.m_index; seen_all && i < sz; ++i) { + justified_var const& jv = next[i]; + signed_var v3 = jv.m_var; + if (!m_marked[v3.index()]) { + seen_all = false; + f.m_index = i + 1; + m_todo.push_back(var_frame(v3, 0)); + m_jtrail.push_back(jv.m_j); + m_marked_trail.push_back(v3.index()); + m_marked[v3.index()] = true; } } + if (seen_all) { + m_todo.pop_back(); + m_jtrail.pop_back(); + } + } - for (eq_justification const& j : m_jtrail) { - j.explain(e); - } - m_stats.m_num_explains += m_jtrail.size(); - m_stats.m_num_explain_calls++; - m_todo.reset(); - m_jtrail.reset(); - for (unsigned idx : m_marked_trail) { - m_marked[idx] = false; - } - m_marked_trail.reset(); - - // IF_VERBOSE(2, verbose_stream() << (double)m_stats.m_num_explains / m_stats.m_num_explain_calls << "\n"); + for (eq_justification const& j : m_jtrail) { + j.explain(e); } + m_stats.m_num_explains += m_jtrail.size(); + m_stats.m_num_explain_calls++; + m_todo.reset(); + m_jtrail.reset(); + for (unsigned idx : m_marked_trail) { + m_marked[idx] = false; + } + m_marked_trail.reset(); - void var_eqs::explain_bfs(signed_var v1, signed_var v2, lp::explanation& e) const { - SASSERT(find(v1) == find(v2)); - if (v1 == v2) { - return; + // IF_VERBOSE(2, verbose_stream() << (double)m_stats.m_num_explains / m_stats.m_num_explain_calls << "\n"); +} + +void var_eqs::explain_bfs(signed_var v1, signed_var v2, lp::explanation& e) const { + SASSERT(find(v1) == find(v2)); + if (v1 == v2) { + return; + } + m_todo.push_back(var_frame(v1, 0)); + m_jtrail.push_back(eq_justification({})); + m_marked.reserve(m_eqs.size(), false); + SASSERT(m_marked_trail.empty()); + m_marked[v1.index()] = true; + m_marked_trail.push_back(v1.index()); + unsigned head = 0; + for (; ; ++head) { + var_frame& f = m_todo[head]; + signed_var v = f.m_var; + if (v == v2) { + break; } - m_todo.push_back(var_frame(v1, 0)); - m_jtrail.push_back(eq_justification({})); - m_marked.reserve(m_eqs.size(), false); - SASSERT(m_marked_trail.empty()); - m_marked[v1.index()] = true; - m_marked_trail.push_back(v1.index()); - unsigned head = 0; - for (; ; ++head) { - var_frame& f = m_todo[head]; - signed_var v = f.m_var; - if (v == v2) { - break; - } - auto const& next = m_eqs[v.index()]; - unsigned sz = next.size(); - for (unsigned i = sz; i-- > 0; ) { - justified_var const& jv = next[i]; - signed_var v3 = jv.m_var; - if (!m_marked[v3.index()]) { - m_todo.push_back(var_frame(v3, head)); - m_jtrail.push_back(jv.m_j); - m_marked_trail.push_back(v3.index()); - m_marked[v3.index()] = true; - } + auto const& next = m_eqs[v.index()]; + unsigned sz = next.size(); + for (unsigned i = sz; i-- > 0; ) { + justified_var const& jv = next[i]; + signed_var v3 = jv.m_var; + if (!m_marked[v3.index()]) { + m_todo.push_back(var_frame(v3, head)); + m_jtrail.push_back(jv.m_j); + m_marked_trail.push_back(v3.index()); + m_marked[v3.index()] = true; } } + } - while (head != 0) { - m_jtrail[head].explain(e); - head = m_todo[head].m_index; - ++m_stats.m_num_explains; - } - ++m_stats.m_num_explain_calls; + while (head != 0) { + m_jtrail[head].explain(e); + head = m_todo[head].m_index; + ++m_stats.m_num_explains; + } + ++m_stats.m_num_explain_calls; - m_todo.reset(); - m_jtrail.reset(); - for (unsigned idx : m_marked_trail) { - m_marked[idx] = false; - } - m_marked_trail.reset(); - - // IF_VERBOSE(2, verbose_stream() << (double)m_stats.m_num_explains / m_stats.m_num_explain_calls << "\n"); + m_todo.reset(); + m_jtrail.reset(); + for (unsigned idx : m_marked_trail) { + m_marked[idx] = false; } + m_marked_trail.reset(); + + // IF_VERBOSE(2, verbose_stream() << (double)m_stats.m_num_explains / m_stats.m_num_explain_calls << "\n"); +} - std::ostream& var_eqs::display(std::ostream& out) const { - m_uf.display(out); - unsigned idx = 0; - for (auto const& edges : m_eqs) { - out << signed_var(idx, from_index()) << ": "; - for (auto const& jv : edges) { - out << jv.m_var << " "; - } - out << "\n"; - ++idx; +std::ostream& var_eqs::display(std::ostream& out) const { + m_uf.display(out); + unsigned idx = 0; + for (auto const& edges : m_eqs) { + out << signed_var(idx, from_index_dummy()) << ": "; + for (auto const& jv : edges) { + out << jv.m_var << " "; } - return out; + out << "\n"; + ++idx; } + return out; +} } diff --git a/src/util/lp/var_eqs.h b/src/util/lp/var_eqs.h index 43475de59..f7469ba54 100644 --- a/src/util/lp/var_eqs.h +++ b/src/util/lp/var_eqs.h @@ -24,41 +24,6 @@ #include "util/rational.h" #include "util/lp/explanation.h" namespace nla { - -typedef lp::var_index lpvar; -typedef lp::constraint_index lpci; -struct from_index{}; - -class signed_var { - unsigned m_sv; -public: - // constructor, sign = true means minus - signed_var(lpvar v, bool sign): m_sv((v << 1) + (sign ? 1 : 0)) {} - // constructor - signed_var(unsigned sv, from_index): m_sv(sv) {} - bool sign() const { return 0 != (m_sv & 0x1); } - lpvar var() const { return m_sv >> 1; } - unsigned index() const { return m_sv; } - void neg() { m_sv = m_sv ^ 1; } - friend signed_var operator~(signed_var const& sv) { - return signed_var(sv.var(), !sv.sign()); - } - bool operator==(signed_var const& other) const { - return m_sv == other.m_sv; - } - bool operator!=(signed_var const& other) const { - return m_sv != other.m_sv; - } - rational rsign() const { return sign() ? rational::minus_one() : rational::one(); } - - std::ostream& display(std::ostream& out) const { - return out << (sign()?"-":"") << var(); - } -}; - - inline std::ostream& operator<<(std::ostream& out, signed_var const& sv) { - return sv.display(out); - } class eq_justification { lpci m_cs[4]; @@ -80,6 +45,16 @@ public: } }; + class var_eqs_merge_handler { + public: + virtual void merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) = 0; + + virtual void after_merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) = 0; + + virtual void unmerge_eh(signed_var r2, signed_var r1) = 0;; + + }; + class var_eqs { struct justified_var { signed_var m_var; @@ -101,12 +76,14 @@ class var_eqs { typedef std::pair signed_var_pair; - union_find_default_ctx m_ufctx; - union_find<> m_uf; + var_eqs_merge_handler* m_merge_handler; + union_find m_uf; svector m_trail; unsigned_vector m_trail_lim; vector m_eqs; // signed_var-index -> justified_var corresponding to edges in a graph used to extract justifications. + typedef trail_stack trail_stack_t; + trail_stack_t m_stack; mutable svector m_todo; mutable svector m_marked; mutable unsigned_vector m_marked_trail; @@ -197,8 +174,32 @@ public: eq_class equiv_class(lpvar v) { return equiv_class(signed_var(v, false)); } + std::ostream& display(std::ostream& out) const; - + + // union find event handlers + void set_merge_handler(var_eqs_merge_handler* mh) { m_merge_handler = mh; } + + trail_stack_t& get_trail_stack() { return m_stack; } + + void unmerge_eh(unsigned i, unsigned j) { + if (m_merge_handler) { + m_merge_handler->unmerge_eh(signed_var(i, from_index_dummy()), signed_var(j, from_index_dummy())); + } + } + void merge_eh(unsigned r2, unsigned r1, unsigned v2, unsigned v1) { + if (m_merge_handler) { + m_merge_handler->merge_eh(signed_var(r2, from_index_dummy()), signed_var(r1, from_index_dummy()), signed_var(v2, from_index_dummy()), signed_var(v1, from_index_dummy())); + } + } + + void after_merge_eh(unsigned r2, unsigned r1, unsigned v2, unsigned v1) { + if (m_merge_handler) { + m_merge_handler->after_merge_eh(signed_var(r2, from_index_dummy()), signed_var(r1, from_index_dummy()), signed_var(v2, from_index_dummy()), signed_var(v1, from_index_dummy())); + } + } + + }; inline std::ostream& operator<<(var_eqs const& ve, std::ostream& out) { return ve.display(out); }