diff --git a/src/muz_qe/hilbert_basis.cpp b/src/muz_qe/hilbert_basis.cpp index b452fbc77..bc31572ea 100644 --- a/src/muz_qe/hilbert_basis.cpp +++ b/src/muz_qe/hilbert_basis.cpp @@ -21,8 +21,6 @@ Revision History: #include "heap.h" #include "map.h" -typedef u_map offset_refs_t; - template class rational_map : public map {}; @@ -132,7 +130,6 @@ private: - class hilbert_basis::index { // for each non-positive weight a separate index. // for positive weights a shared value index. @@ -142,7 +139,8 @@ class hilbert_basis::index { unsigned m_num_insert; stats() { reset(); } void reset() { memset(this, 0, sizeof(*this)); } - }; + }; + typedef rational_map value_map; hilbert_basis& hb; value_map m_neg; @@ -249,16 +247,15 @@ class hilbert_basis::passive { struct lt { passive& p; lt(passive& p): p(p) {} - bool operator()(int v1, int v2) const { return p(v1, v2); } }; - hilbert_basis& hb; - svector m_passive; - unsigned_vector m_free_list; - lt m_lt; - heap m_heap; // binary heap over weights + hilbert_basis& hb; + svector m_passive; + unsigned_vector m_free_list; + lt m_lt; + heap m_heap; // binary heap over weights numeral get_value(offset_t idx) const { numeral w(0); @@ -275,7 +272,8 @@ public: hb(hb) , m_lt(*this), m_heap(10, m_lt) - {} + { + } void reset() { m_heap.reset(); @@ -412,6 +410,7 @@ void hilbert_basis::collect_statistics(statistics& st) const { st.update("hb.num_subsumptions", m_stats.m_num_subsumptions); st.update("hb.num_resolves", m_stats.m_num_resolves); st.update("hb.num_saturations", m_stats.m_num_saturations); + st.update("hb.basis_size", get_basis_size()); m_index->collect_statistics(st); } @@ -426,6 +425,7 @@ void hilbert_basis::add_ge(num_vector const& v, numeral const& b) { w.push_back(-b); w.append(v); m_ineqs.push_back(w); + m_iseq.push_back(false); } void hilbert_basis::add_le(num_vector const& v, numeral const& b) { @@ -437,8 +437,12 @@ void hilbert_basis::add_le(num_vector const& v, numeral const& b) { } void hilbert_basis::add_eq(num_vector const& v, numeral const& b) { - add_le(v, b); - add_ge(v, b); + SASSERT(m_ineqs.empty() || v.size() + 1 == get_num_vars()); + num_vector w; + w.push_back(-b); + w.append(v); + m_ineqs.push_back(w); + m_iseq.push_back(true); } void hilbert_basis::add_ge(num_vector const& v) { @@ -450,8 +454,7 @@ void hilbert_basis::add_le(num_vector const& v) { } void hilbert_basis::add_eq(num_vector const& v) { - add_le(v); - add_ge(v); + add_eq(v, numeral(0)); } void hilbert_basis::set_is_int(unsigned var_index) { @@ -462,6 +465,10 @@ void hilbert_basis::set_is_int(unsigned var_index) { m_ints.push_back(var_index+1); } +bool hilbert_basis::get_is_int(unsigned var_index) const { + return m_ints.contains(var_index+1); +} + unsigned hilbert_basis::get_num_vars() const { if (m_ineqs.empty()) { return 0; @@ -501,14 +508,16 @@ void hilbert_basis::add_unit_vector(unsigned i, numeral const& e) { } lbool hilbert_basis::saturate() { - init_basis(); - for (unsigned i = 0; !m_cancel && i < m_ineqs.size(); ++i) { - select_inequality(i); - lbool r = saturate(m_ineqs[i]); + init_basis(); + m_current_ineq = 0; + while (!m_cancel && m_current_ineq < m_ineqs.size()) { + select_inequality(); + lbool r = saturate(m_ineqs[m_current_ineq], m_iseq[m_current_ineq]); ++m_stats.m_num_saturations; if (r != l_true) { return r; } + ++m_current_ineq; } if (m_cancel) { return l_undef; @@ -516,17 +525,17 @@ lbool hilbert_basis::saturate() { return l_true; } -lbool hilbert_basis::saturate(num_vector const& ineq) { +lbool hilbert_basis::saturate(num_vector const& ineq, bool is_eq) { m_active.reset(); m_passive->reset(); m_zero.reset(); m_index->reset(); - TRACE("hilbert_basis", display_ineq(tout, ineq);); + TRACE("hilbert_basis", display_ineq(tout, ineq, is_eq);); bool has_non_negative = false; iterator it = begin(); for (; it != end(); ++it) { values v = vec(*it); - set_eval(v, ineq); + v.weight() = get_weight(v, ineq); add_goal(*it); if (v.weight().is_nonneg()) { has_non_negative = true; @@ -559,7 +568,7 @@ lbool hilbert_basis::saturate(num_vector const& ineq) { // Move positive from active and zeros to new basis. m_basis.reset(); m_basis.append(m_zero); - for (unsigned i = 0; i < m_active.size(); ++i) { + for (unsigned i = 0; !is_eq && i < m_active.size(); ++i) { offset_t idx = m_active[i]; if (vec(idx).weight().is_pos()) { m_basis.push_back(idx); @@ -575,12 +584,29 @@ lbool hilbert_basis::saturate(num_vector const& ineq) { return l_true; } -void hilbert_basis::select_inequality(unsigned i) { - SASSERT(i < m_ineqs.size()); - unsigned best = i; - unsigned non_zeros = get_num_nonzeros(m_ineqs[i]); - unsigned prod = get_ineq_product(m_ineqs[i]); - for (unsigned j = i+1; prod != 0 && j < m_ineqs.size(); ++j) { +void hilbert_basis::get_basis_solution(unsigned i, num_vector& v, bool& is_initial) { + offset_t offs = m_basis[i]; + v.reset(); + for (unsigned i = 1; i < get_num_vars(); ++i) { + v.push_back(vec(offs)[i]); + } + is_initial = !vec(offs)[0].is_zero(); +} + +void hilbert_basis::get_ge(unsigned i, num_vector& v, numeral& b, bool& is_eq) { + v.reset(); + v.append(get_num_vars()-1, m_ineqs[i].c_ptr() + 1); + b = -m_ineqs[i][0]; + is_eq = m_iseq[i]; +} + + +void hilbert_basis::select_inequality() { + SASSERT(m_current_ineq < m_ineqs.size()); + unsigned best = m_current_ineq; + unsigned non_zeros = get_num_nonzeros(m_ineqs[best]); + unsigned prod = get_ineq_product(m_ineqs[best]); + for (unsigned j = best+1; prod != 0 && j < m_ineqs.size(); ++j) { unsigned non_zeros2 = get_num_nonzeros(m_ineqs[j]); unsigned prod2 = get_ineq_product(m_ineqs[j]); if (prod2 < prod || (prod2 == prod && non_zeros2 < non_zeros)) { @@ -589,8 +615,9 @@ void hilbert_basis::select_inequality(unsigned i) { best = j; } } - if (best != i) { - std::swap(m_ineqs[i], m_ineqs[best]); + if (best != m_current_ineq) { + std::swap(m_ineqs[m_current_ineq], m_ineqs[best]); + std::swap(m_iseq[m_current_ineq], m_iseq[best]); } } @@ -609,11 +636,11 @@ unsigned hilbert_basis::get_ineq_product(num_vector const& ineq) { iterator it = begin(); for (; it != end(); ++it) { values v = vec(*it); - set_eval(v, ineq); - if (v.weight().is_pos()) { + numeral w = get_weight(v, ineq); + if (w.is_pos()) { ++num_pos; } - else if (v.weight().is_neg()) { + else if (w.is_neg()) { ++num_neg; } } @@ -716,20 +743,20 @@ hilbert_basis::sign_t hilbert_basis::get_sign(offset_t idx) const { return zero; } -void hilbert_basis::set_eval(values& val, num_vector const& ineq) const { +hilbert_basis::numeral hilbert_basis::get_weight(values& val, num_vector const& ineq) const { numeral result(0); unsigned num_vars = get_num_vars(); for (unsigned i = 0; i < num_vars; ++i) { result += val[i]*ineq[i]; } - val.weight() = result; + return result; } void hilbert_basis::display(std::ostream& out) const { unsigned nv = get_num_vars(); out << "inequalities:\n"; for (unsigned i = 0; i < m_ineqs.size(); ++i) { - display_ineq(out, m_ineqs[i]); + display_ineq(out, m_ineqs[i], m_iseq[i]); } if (!m_basis.empty()) { out << "basis:\n"; @@ -757,7 +784,6 @@ void hilbert_basis::display(std::ostream& out) const { display(out, m_zero[i]); } } - } void hilbert_basis::display(std::ostream& out, offset_t o) const { @@ -772,7 +798,7 @@ void hilbert_basis::display(std::ostream& out, values const& v) const { } } -void hilbert_basis::display_ineq(std::ostream& out, num_vector const& v) const { +void hilbert_basis::display_ineq(std::ostream& out, num_vector const& v, bool is_eq) const { unsigned nv = get_num_vars(); for (unsigned j = 0; j < nv; ++j) { if (!v[j].is_zero()) { @@ -793,20 +819,12 @@ void hilbert_basis::display_ineq(std::ostream& out, num_vector const& v) const { out << "x" << j; } } - out << " >= 0\n"; -} - - -void hilbert_isl_basis::add_le(num_vector const& v, numeral bound) { - unsigned sz = v.size(); - num_vector w; - w.push_back(-bound); - w.push_back(bound); - for (unsigned i = 0; i < sz; ++i) { - w.push_back(v[i]); - w.push_back(-v[i]); + if (is_eq) { + out << " = 0\n"; + } + else { + out << " >= 0\n"; } - m_basis.add_le(w); } @@ -847,6 +865,9 @@ bool hilbert_basis::is_subsumed(offset_t i, offset_t j) const { i.m_offset != j.m_offset && n >= m && (!m.is_nonpos() || n == m) && is_geq(v, w); + for (unsigned k = 0; r && k < m_current_ineq; ++k) { + r = get_weight(vec(i), m_ineqs[k]) >= get_weight(vec(j), m_ineqs[k]); + } CTRACE("hilbert_basis", r, display(tout, i); tout << " <= \n"; diff --git a/src/muz_qe/hilbert_basis.h b/src/muz_qe/hilbert_basis.h index 000f27ec8..00b854271 100644 --- a/src/muz_qe/hilbert_basis.h +++ b/src/muz_qe/hilbert_basis.h @@ -11,8 +11,6 @@ Abstract: hilbert_basis computes a Hilbert basis for linear homogeneous inequalities over naturals. - hilbert_sl_basis computes a semi-linear set over naturals. - hilbert_isl_basis computes semi-linear sets over integers. Author: @@ -64,6 +62,7 @@ private: }; vector m_ineqs; // set of asserted inequalities + svector m_iseq; // inequalities that are equalities num_vector m_store; // store of vectors svector m_basis; // vector of current basis svector m_free_list; // free list of unused storage @@ -74,6 +73,7 @@ private: stats m_stats; index* m_index; // index of generated vectors unsigned_vector m_ints; // indices that can be both positive and negative + unsigned m_current_ineq; class iterator { hilbert_basis const& hb; unsigned m_idx; @@ -88,15 +88,15 @@ private: static offset_t mk_invalid_offset(); static bool is_invalid_offset(offset_t offs); - lbool saturate(num_vector const& ineq); + lbool saturate(num_vector const& ineq, bool is_eq); void init_basis(); - void select_inequality(unsigned i); + void select_inequality(); unsigned get_num_nonzeros(num_vector const& ineq); unsigned get_ineq_product(num_vector const& ineq); void add_unit_vector(unsigned i, numeral const& e); unsigned get_num_vars() const; - void set_eval(values& val, num_vector const& ineq) const; + numeral get_weight(values& val, num_vector const& ineq) const; bool is_geq(values const& v, values const& w) const; bool is_abs_geq(numeral const& v, numeral const& w) const; bool is_subsumed(offset_t idx); @@ -114,7 +114,7 @@ private: void display(std::ostream& out, offset_t o) const; void display(std::ostream& out, values const & v) const; - void display_ineq(std::ostream& out, num_vector const& v) const; + void display_ineq(std::ostream& out, num_vector const& v, bool is_eq) const; public: @@ -138,33 +138,24 @@ public: void add_eq(num_vector const& v, numeral const& b); void set_is_int(unsigned var_index); + bool get_is_int(unsigned var_index) const; lbool saturate(); + unsigned get_basis_size() const { return m_basis.size(); } + void get_basis_solution(unsigned i, num_vector& v, bool& is_initial); + + unsigned get_num_ineqs() const { return m_ineqs.size(); } + void get_ge(unsigned i, num_vector& v, numeral& b, bool& is_eq); + void set_cancel(bool f) { m_cancel = f; } void display(std::ostream& out) const; void collect_statistics(statistics& st) const; void reset_statistics(); + }; -class hilbert_isl_basis { -public: - typedef rational numeral; - typedef vector num_vector; -private: - hilbert_basis m_basis; -public: - hilbert_isl_basis() {} - void reset() { m_basis.reset(); } - - // add inequality v*x >= bound, x ranges over integers - void add_le(num_vector const& v, numeral bound); - lbool saturate() { return m_basis.saturate(); } - void set_cancel(bool f) { m_basis.set_cancel(f); } - void display(std::ostream& out) const { m_basis.display(out); } -}; - #endif diff --git a/src/muz_qe/hilbert_basis_validate.cpp b/src/muz_qe/hilbert_basis_validate.cpp new file mode 100644 index 000000000..df65146f1 --- /dev/null +++ b/src/muz_qe/hilbert_basis_validate.cpp @@ -0,0 +1,178 @@ +/*++ +Copyright (c) 2013 Microsoft Corporation + +Module Name: + + hilbert_basis_validate.cpp + +Abstract: + + Basic Hilbert Basis validation. + + hilbert_basis computes a Hilbert basis for linear + homogeneous inequalities over naturals. + +Author: + + Nikolaj Bjorner (nbjorner) 2013-02-15. + +Revision History: + +--*/ + +#include "hilbert_basis_validate.h" +#include "arith_decl_plugin.h" +#include "ast_pp.h" +#include + + +hilbert_basis_validate::hilbert_basis_validate(ast_manager& m): + m(m) { +} + +void hilbert_basis_validate::validate_solution(hilbert_basis& hb, vector const& v, bool is_initial) { + unsigned sz = hb.get_num_ineqs(); + rational bound; + for (unsigned i = 0; i < sz; ++i) { + bool is_eq; + vector w; + hb.get_ge(i, w, bound, is_eq); + rational sum(0); + for (unsigned j = 0; j < v.size(); ++j) { + sum += w[j]*v[j]; + } + if (bound > sum || + (is_eq && bound != sum)) { + // validation failed. + std::cout << "validation failed for inequality\n"; + for (unsigned j = 0; j < v.size(); ++j) { + std::cout << v[j] << " "; + } + std::cout << "\n"; + for (unsigned j = 0; j < w.size(); ++j) { + std::cout << w[j] << " "; + } + std::cout << (is_eq?" = ":" >= ") << bound << "\n"; + std::cout << "is initial: " << (is_initial?"true":"false") << "\n"; + std::cout << "sum: " << sum << "\n"; + } + } +} + +expr_ref hilbert_basis_validate::mk_validate(hilbert_basis& hb) { + arith_util a(m); + unsigned sz = hb.get_basis_size(); + vector v; + bool is_initial; + + // check that claimed solution really satisfies inequalities: + for (unsigned i = 0; i < sz; ++i) { + hb.get_basis_solution(i, v, is_initial); + validate_solution(hb, v, is_initial); + } + + // check that solutions satisfying inequalities are in solution. + // build a formula that says solutions to linear inequalities + // coincide with linear combinations of basis. + vector offsets, increments; + expr_ref_vector xs(m), vars(m); + expr_ref var(m); + svector names; + sort_ref_vector sorts(m); + +#define mk_mul(_r,_x) (_r.is_one()?((expr*)_x):((expr*)a.mk_mul(a.mk_numeral(_r,true),_x))) + + + for (unsigned i = 0; i < sz; ++i) { + hb.get_basis_solution(i, v, is_initial); + + for (unsigned j = 0; xs.size() < v.size(); ++j) { + xs.push_back(m.mk_fresh_const("x", a.mk_int())); + } + + if (is_initial) { + expr_ref_vector tmp(m); + for (unsigned j = 0; j < v.size(); ++j) { + tmp.push_back(a.mk_numeral(v[j], true)); + } + offsets.push_back(tmp); + } + else { + var = m.mk_var(vars.size(), a.mk_int()); + expr_ref_vector tmp(m); + for (unsigned j = 0; j < v.size(); ++j) { + tmp.push_back(mk_mul(v[j], var)); + } + std::stringstream name; + name << "u" << i; + increments.push_back(tmp); + vars.push_back(var); + names.push_back(symbol(name.str().c_str())); + sorts.push_back(a.mk_int()); + } + } + + expr_ref_vector bounds(m); + for (unsigned i = 0; i < vars.size(); ++i) { + bounds.push_back(a.mk_ge(vars[i].get(), a.mk_numeral(rational(0), true))); + } + expr_ref_vector fmls(m); + expr_ref fml(m), fml1(m), fml2(m); + for (unsigned i = 0; i < offsets.size(); ++i) { + expr_ref_vector eqs(m); + eqs.append(bounds); + for (unsigned j = 0; j < xs.size(); ++j) { + expr_ref_vector sum(m); + sum.push_back(offsets[i][j].get()); + for (unsigned k = 0; k < increments.size(); ++k) { + sum.push_back(increments[k][j].get()); + } + eqs.push_back(m.mk_eq(xs[j].get(), a.mk_add(sum.size(), sum.c_ptr()))); + } + fml = m.mk_and(eqs.size(), eqs.c_ptr()); + if (!names.empty()) { + fml = m.mk_exists(names.size(), sorts.c_ptr(), names.c_ptr(), fml); + } + fmls.push_back(fml); + } + fml1 = m.mk_or(fmls.size(), fmls.c_ptr()); + fmls.reset(); + + sz = hb.get_num_ineqs(); + for (unsigned i = 0; i < sz; ++i) { + bool is_eq; + vector w; + rational bound; + hb.get_ge(i, w, bound, is_eq); + expr_ref_vector sum(m); + for (unsigned j = 0; j < w.size(); ++j) { + if (!w[j].is_zero()) { + sum.push_back(mk_mul(w[j], xs[j].get())); + } + } + expr_ref lhs(m), rhs(m); + lhs = a.mk_add(sum.size(), sum.c_ptr()); + rhs = a.mk_numeral(bound, true); + if (is_eq) { + fmls.push_back(a.mk_eq(lhs, rhs)); + } + else { + fmls.push_back(a.mk_ge(lhs, rhs)); + } + } + fml2 = m.mk_and(fmls.size(), fmls.c_ptr()); + fml = m.mk_eq(fml1, fml2); + + bounds.reset(); + for (unsigned i = 0; i < xs.size(); ++i) { + if (!hb.get_is_int(i)) { + bounds.push_back(a.mk_ge(xs[i].get(), a.mk_numeral(rational(0), true))); + } + } + if (!bounds.empty()) { + fml = m.mk_implies(m.mk_and(bounds.size(), bounds.c_ptr()), fml); + } + return fml; + +} + diff --git a/src/muz_qe/hilbert_basis_validate.h b/src/muz_qe/hilbert_basis_validate.h new file mode 100644 index 000000000..defa5805c --- /dev/null +++ b/src/muz_qe/hilbert_basis_validate.h @@ -0,0 +1,43 @@ +/*++ +Copyright (c) 2013 Microsoft Corporation + +Module Name: + + hilbert_basis_validate.h + +Abstract: + + Basic Hilbert Basis validation. + + hilbert_basis computes a Hilbert basis for linear + homogeneous inequalities over naturals. + +Author: + + Nikolaj Bjorner (nbjorner) 2013-02-15. + +Revision History: + +--*/ + +#ifndef _HILBERT_BASIS_VALIDATE_H_ +#define _HILBERT_BASIS_VALIDATE_H_ + +#include "hilbert_basis.h" +#include "ast.h" + +class hilbert_basis_validate { + ast_manager& m; + + void validate_solution(hilbert_basis& hb, vector const& v, bool is_initial); + +public: + + hilbert_basis_validate(ast_manager& m); + + expr_ref mk_validate(hilbert_basis& hb); + +}; + + +#endif diff --git a/src/test/hilbert_basis.cpp b/src/test/hilbert_basis.cpp index 8629e50be..11faad819 100644 --- a/src/test/hilbert_basis.cpp +++ b/src/test/hilbert_basis.cpp @@ -1,4 +1,12 @@ #include "hilbert_basis.h" +#include "hilbert_basis_validate.h" +#include "ast_pp.h" +#include "reg_decl_plugins.h" +#include "quant_tactics.h" +#include "tactic.h" +#include "tactic2solver.h" +#include "solver.h" + #include #include @@ -19,6 +27,24 @@ static void on_ctrl_c(int) { raise(SIGINT); } +static void validate_sat(hilbert_basis& hb) { + ast_manager m; + reg_decl_plugins(m); + hilbert_basis_validate val(m); + + expr_ref fml = val.mk_validate(hb); + + std::cout << mk_pp(fml, m) << "\n"; + + fml = m.mk_not(fml); + params_ref p; + tactic_ref tac = mk_lra_tactic(m, p); + ref sol = mk_tactic2solver(m, tac.get(), p); + sol->assert_expr(fml); + lbool r = sol->check_sat(0,0); + std::cout << r << "\n"; +} + static void saturate_basis(hilbert_basis& hb) { signal(SIGINT, on_ctrl_c); g_hb = &hb; @@ -29,6 +55,7 @@ static void saturate_basis(hilbert_basis& hb) { case l_true: std::cout << "sat\n"; hb.display(std::cout); + // validate_sat(hb); break; case l_false: std::cout << "unsat\n"; @@ -40,6 +67,7 @@ static void saturate_basis(hilbert_basis& hb) { display_statistics(hb); } + /** n - number of variables. k - subset of variables to be non-zero @@ -74,6 +102,14 @@ static void gorrila_test(unsigned seed, unsigned n, unsigned k, unsigned bound, saturate_basis(hb); } +static vector vec(int i, int j) { + vector nv; + nv.resize(2); + nv[0] = rational(i); + nv[1] = rational(j); + return nv; +} + static vector vec(int i, int j, int k) { vector nv; nv.resize(3); @@ -243,16 +279,44 @@ static void tst11() { saturate_basis(hb); } +static void tst12() { + hilbert_basis hb; + hb.add_le(vec(1, 0), R(1)); + hb.add_le(vec(0, 1), R(1)); + saturate_basis(hb); +} + +// Sigma_9 table 1, Ajili, Contejean +static void tst13() { + hilbert_basis hb; + hb.add_eq(vec( 1,-2,-4,4), R(0)); + hb.add_le(vec(100,45,-78,-67), R(0)); + saturate_basis(hb); +} + +// Sigma_10 table 1, Ajili, Contejean +static void tst14() { + hilbert_basis hb; + hb.add_le(vec( 23, -56, -34, 12, 11), R(0)); + saturate_basis(hb); +} + +// Sigma_11 table 1, Ajili, Contejean +static void tst15() { +// hilbert_basis hb; +// hb.add_le(vec( 23, -56, -34, 12, 11), R(0)); +// saturate_basis(hb); +} + + void tst_hilbert_basis() { std::cout << "hilbert basis test\n"; - tst4(); - return; if (true) { tst1(); tst2(); tst3(); - tst4(); + // tst4(); tst5(); tst6(); tst7(); @@ -260,11 +324,15 @@ void tst_hilbert_basis() { tst9(); tst10(); tst11(); + tst12(); + tst13(); + tst14(); + tst15(); gorrila_test(0, 4, 3, 20, 5); gorrila_test(1, 4, 3, 20, 5); - gorrila_test(2, 4, 3, 20, 5); - gorrila_test(0, 4, 2, 20, 5); - gorrila_test(0, 4, 2, 20, 5); + //gorrila_test(2, 4, 3, 20, 5); + //gorrila_test(0, 4, 2, 20, 5); + //gorrila_test(0, 4, 2, 20, 5); } else { gorrila_test(0, 10, 7, 20, 11);