z3-z3-4.13.0.src.math.hilbert.hilbert_basis.cpp Maven / Gradle / Ivy
The newest version!
/*++
Copyright (c) 2013 Microsoft Corporation
Module Name:
hilbert_basis.cpp
Abstract:
Basic Hilbert Basis computation.
Author:
Nikolaj Bjorner (nbjorner) 2013-02-09.
Revision History:
--*/
#include "math/hilbert/hilbert_basis.h"
#include "util/heap.h"
#include "util/map.h"
#include "math/hilbert/heap_trie.h"
#include "util/stopwatch.h"
typedef int_hashtable > int_table;
class hilbert_basis::value_index1 {
struct stats {
unsigned m_num_comparisons;
unsigned m_num_hit;
unsigned m_num_miss;
stats() { reset(); }
void reset() { memset(this, 0, sizeof(*this)); }
};
hilbert_basis& hb;
int_table m_table;
stats m_stats;
public:
value_index1(hilbert_basis& hb):
hb(hb)
{}
void insert(offset_t idx, values const& vs) {
m_table.insert(idx.m_offset);
}
void remove(offset_t idx, values const& vs) {
m_table.remove(idx.m_offset);
}
void reset() {
m_table.reset();
}
bool find(offset_t idx, values const& vs) {
// display_profile(idx, std::cout);
int_table::iterator it = m_table.begin(), end = m_table.end();
for (; it != end; ++it) {
offset_t offs(*it);
++m_stats.m_num_comparisons;
if (*it != static_cast(idx.m_offset) && hb.is_subsumed(idx, offs)) {
++m_stats.m_num_hit;
return true;
}
}
++m_stats.m_num_miss;
return false;
}
void collect_statistics(statistics& st) const {
st.update("hb.index.num_comparisons", m_stats.m_num_comparisons);
st.update("hb.index.num_hit", m_stats.m_num_hit);
st.update("hb.index.num_miss", m_stats.m_num_miss);
}
void reset_statistics() {
m_stats.reset();
}
unsigned size() const {
return m_table.size();
}
void display(std::ostream& out) const {
int_table::iterator it = m_table.begin(), end = m_table.end();
for (; it != end; ++it) {
offset_t offs(*it);
hb.display(out, offs);
}
display_profile(out);
}
private:
typedef hashtable numeral_set;
void display_profile(std::ostream& out) const {
vector weights;
weights.resize(hb.get_num_vars()+1);
int_table::iterator it = m_table.begin(), end = m_table.end();
for (; it != end; ++it) {
offset_t offs(*it);
values const& ws = hb.vec(offs);
weights[0].insert(ws.weight());
for (unsigned i = 0; i < hb.get_num_vars(); ++i) {
weights[i+1].insert(ws[i]);
}
}
out << "profile: ";
for (unsigned i = 0; i < weights.size(); ++i) {
out << weights[i].size() << " ";
}
out << "\n";
}
void display_profile(offset_t idx, std::ostream& out) {
unsigned_vector leq;
unsigned nv = hb.get_num_vars();
values const& vs = hb.vec(idx);
leq.resize(nv+1);
numeral maxw(0);
for (unsigned i = 0; i < nv; ++i) {
if (!hb.is_abs_geq(maxw, vs[i])) {
maxw = vs[i];
}
}
unsigned num_below_max = 0;
int_table::iterator it = m_table.begin(), end = m_table.end();
for (; it != end; ++it) {
offset_t offs(*it);
values const& ws = hb.vec(offs);
if (ws.weight() <= vs.weight()) {
leq[0]++;
}
bool filtered = false;
for (unsigned i = 0; i < nv; ++i) {
if (hb.is_abs_geq(vs[i], ws[i])) {
leq[i+1]++;
}
if (!hb.is_abs_geq(maxw, ws[i])) {
filtered = true;
}
}
if (!filtered) {
++num_below_max;
}
}
out << vs.weight() << ":" << leq[0] << " ";
for (unsigned i = 0; i < nv; ++i) {
out << vs[i] << ":" << leq[i+1] << " ";
}
out << " max<= " << num_below_max;
out << "\n";
}
};
class hilbert_basis::value_index2 {
struct key_le {
hilbert_basis& hb;
key_le(hilbert_basis& hb): hb(hb) {}
bool le(numeral const& n1, numeral const& n2) const {
return hb.is_abs_geq(n2, n1);
}
};
typedef heap_trie ht;
struct checker : public ht::check_value {
hilbert_basis* hb;
offset_t m_value;
checker(): hb(nullptr) {}
bool operator()(unsigned const& v) override {
if (m_value.m_offset != v) { // && hb->is_subsumed(m_value, offset_t(v))) {
return true;
}
else {
return false;
}
}
};
hilbert_basis& hb;
key_le m_le;
ht m_trie;
checker m_checker;
unsigned m_offset;
numeral const* get_keys(values const& vs) {
return vs()-m_offset;
}
public:
value_index2(hilbert_basis& hb):
hb(hb),
m_le(hb),
m_trie(m_le),
m_offset(1) {
m_checker.hb = &hb;
}
void insert(offset_t idx, values const& vs) {
m_trie.insert(get_keys(vs), idx.m_offset);
}
void remove(offset_t idx, values const& vs) {
m_trie.remove(get_keys(vs));
}
void reset(unsigned offset) {
m_offset = offset;
m_trie.reset(hb.get_num_vars()+m_offset);
}
bool find(offset_t idx, values const& vs) {
m_checker.m_value = idx;
return m_trie.find_le(get_keys(vs), m_checker);
}
void collect_statistics(statistics& st) const {
m_trie.collect_statistics(st);
}
void reset_statistics() {
m_trie.reset_statistics();
}
unsigned size() const {
return m_trie.size();
}
void display(std::ostream& out) const {
// m_trie.display(out);
}
};
class hilbert_basis::index {
// for each non-positive weight a separate index.
// for positive weights a shared value index.
// typedef value_index1 value_index;
typedef value_index2 value_index;
// typedef value_index3 value_index;
struct stats {
unsigned m_num_find;
unsigned m_num_insert;
stats() { reset(); }
void reset() { memset(this, 0, sizeof(*this)); }
};
template
class numeral_map : public map {};
typedef numeral_map value_map;
hilbert_basis& hb;
value_map m_neg;
value_index m_pos;
value_index m_zero;
stats m_stats;
unsigned m_num_ineqs;
public:
index(hilbert_basis& hb): hb(hb), m_pos(hb), m_zero(hb), m_num_ineqs(0) {}
void insert(offset_t idx, values const& vs) {
++m_stats.m_num_insert;
if (vs.weight().is_pos()) {
m_pos.insert(idx, vs);
}
else if (vs.weight().is_zero()) {
m_zero.insert(idx, vs);
}
else {
value_index* map = nullptr;
if (!m_neg.find(vs.weight(), map)) {
map = alloc(value_index, hb);
map->reset(m_num_ineqs);
m_neg.insert(vs.weight(), map);
}
map->insert(idx, vs);
}
}
void remove(offset_t idx, values const& vs) {
if (vs.weight().is_pos()) {
m_pos.remove(idx, vs);
}
else if (vs.weight().is_zero()) {
m_zero.remove(idx, vs);
}
else {
m_neg.find(vs.weight())->remove(idx, vs);
}
}
bool find(offset_t idx, values const& vs) {
++m_stats.m_num_find;
if (vs.weight().is_pos()) {
return m_pos.find(idx, vs);
}
else if (vs.weight().is_zero()) {
return m_zero.find(idx, vs);
}
else {
value_index* map;
return
m_neg.find(vs.weight(), map) &&
map->find(idx, vs);
}
}
void reset(unsigned num_ineqs) {
value_map::iterator it = m_neg.begin(), end = m_neg.end();
for (; it != end; ++it) {
dealloc(it->m_value);
}
m_pos.reset(num_ineqs);
m_zero.reset(num_ineqs);
m_num_ineqs = num_ineqs;
m_neg.reset();
}
void collect_statistics(statistics& st) const {
m_pos.collect_statistics(st);
m_zero.collect_statistics(st);
value_map::iterator it = m_neg.begin(), end = m_neg.end();
for (; it != end; ++it) {
it->m_value->collect_statistics(st);
}
st.update("hb.index.num_find", m_stats.m_num_find);
st.update("hb.index.num_insert", m_stats.m_num_insert);
st.update("hb.index.size", size());
}
void reset_statistics() {
m_pos.reset_statistics();
m_zero.reset_statistics();
value_map::iterator it = m_neg.begin(), end = m_neg.end();
for (; it != end; ++it) {
it->m_value->reset_statistics();
}
}
void display(std::ostream& out) const {
m_pos.display(out);
m_zero.display(out);
value_map::iterator it = m_neg.begin(), end = m_neg.end();
for (; it != end; ++it) {
it->m_value->display(out);
}
}
private:
unsigned size() const {
unsigned sz = m_pos.size();
sz += m_zero.size();
value_map::iterator it = m_neg.begin(), end = m_neg.end();
for (; it != end; ++it) {
sz += it->m_value->size();
}
return sz;
}
};
/**
\brief priority queue for passive list.
*/
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;
passive* m_this;
lt m_lt;
heap m_heap; // binary heap over weights
numeral sum_abs(offset_t idx) const {
numeral w(0);
unsigned nv = hb.get_num_vars();
for (unsigned i = 0; i < nv; ++i) {
w += abs(hb.vec(idx)[i]);
}
return w;
}
public:
passive(hilbert_basis& hb):
hb(hb),
m_lt(&m_this),
m_heap(10, m_lt)
{
m_this = this;
}
void reset() {
m_heap.reset();
m_free_list.reset();
m_passive.reset();
}
bool empty() const {
return m_heap.empty();
}
offset_t pop() {
SASSERT(!empty());
unsigned val = static_cast(m_heap.erase_min());
offset_t result = m_passive[val];
m_free_list.push_back(val);
m_passive[val] = mk_invalid_offset();
return result;
}
void insert(offset_t idx) {
unsigned v;
if (m_free_list.empty()) {
v = m_passive.size();
m_passive.push_back(idx);
m_heap.set_bounds(v+1);
}
else {
v = m_free_list.back();
m_free_list.pop_back();
m_passive[v] = idx;
}
m_heap.insert(v);
}
class iterator {
passive& p;
unsigned m_idx;
void fwd() {
while (m_idx < p.m_passive.size() &&
is_invalid_offset(p.m_passive[m_idx])) {
++m_idx;
}
}
public:
iterator(passive& p, unsigned i): p(p), m_idx(i) { fwd(); }
offset_t operator*() const { return p.m_passive[m_idx]; }
iterator& operator++() { ++m_idx; fwd(); return *this; }
iterator operator++(int) { iterator tmp = *this; ++*this; return tmp; }
bool operator==(iterator const& it) const {return m_idx == it.m_idx; }
bool operator!=(iterator const& it) const {return m_idx != it.m_idx; }
};
iterator begin() {
return iterator(*this, 0);
}
iterator end() {
return iterator(*this, m_passive.size());
}
bool operator()(int v1, int v2) const {
offset_t idx1 = m_passive[v1];
offset_t idx2 = m_passive[v2];
return sum_abs(idx1) < sum_abs(idx2);
}
};
class hilbert_basis::vector_lt_t {
hilbert_basis& hb;
public:
vector_lt_t(hilbert_basis& hb): hb(hb) {}
bool operator()(offset_t idx1, offset_t idx2) const {
return hb.vector_lt(idx1, idx2);
}
};
class hilbert_basis::passive2 {
struct lt {
passive2** p;
lt(passive2** p): p(p) {}
bool operator()(int v1, int v2) const {
return (**p)(v1, v2);
}
};
hilbert_basis& hb;
svector m_pos_sos;
svector m_neg_sos;
vector m_pos_sos_sum;
vector m_neg_sos_sum;
vector m_sum_abs;
unsigned_vector m_psos;
svector m_pas;
vector m_weight;
unsigned_vector m_free_list;
passive2* m_this;
lt m_lt;
heap m_heap; // binary heap over weights
numeral sum_abs(offset_t idx) const {
numeral w(0);
unsigned nv = hb.get_num_vars();
for (unsigned i = 0; i < nv; ++i) {
w += abs(hb.vec(idx)[i]);
}
return w;
}
public:
passive2(hilbert_basis& hb):
hb(hb),
m_lt(&m_this),
m_heap(10, m_lt)
{
m_this = this;
}
void init(svector const& I) {
for (unsigned i = 0; i < I.size(); ++i) {
numeral const& w = hb.vec(I[i]).weight();
if (w.is_pos()) {
m_pos_sos.push_back(I[i]);
m_pos_sos_sum.push_back(sum_abs(I[i]));
}
else {
m_neg_sos.push_back(I[i]);
m_neg_sos_sum.push_back(sum_abs(I[i]));
}
}
}
void reset() {
m_heap.reset();
m_free_list.reset();
m_psos.reset();
m_pas.reset();
m_sum_abs.reset();
m_pos_sos.reset();
m_neg_sos.reset();
m_pos_sos_sum.reset();
m_neg_sos_sum.reset();
m_weight.reset();
}
void insert(offset_t idx, unsigned offset) {
SASSERT(!m_pos_sos.empty());
unsigned v;
if (m_free_list.empty()) {
v = m_pas.size();
m_pas.push_back(idx);
m_psos.push_back(offset);
m_weight.push_back(numeral(0));
m_heap.set_bounds(v+1);
m_sum_abs.push_back(sum_abs(idx));
}
else {
v = m_free_list.back();
m_free_list.pop_back();
m_pas[v] = idx;
m_psos[v] = offset;
m_weight[v] = numeral(0);
m_sum_abs[v] = sum_abs(idx);
}
next_resolvable(hb.vec(idx).weight().is_pos(), v);
}
bool empty() const {
return m_heap.empty();
}
unsigned pop(offset_t& sos, offset_t& pas) {
SASSERT (!empty());
unsigned val = static_cast(m_heap.erase_min());
pas = m_pas[val];
numeral old_weight = hb.vec(pas).weight();
bool is_positive = old_weight.is_pos();
unsigned psos = m_psos[val];
sos = is_positive?m_neg_sos[psos]:m_pos_sos[psos];
m_psos[val]++;
next_resolvable(is_positive, val);
numeral new_weight = hb.vec(sos).weight() + old_weight;
if (new_weight.is_pos() != old_weight.is_pos()) {
psos = 0;
}
return psos;
}
bool operator()(int v1, int v2) const {
return m_weight[v1] < m_weight[v2];
}
class iterator {
passive2& p;
unsigned m_idx;
void fwd() {
while (m_idx < p.m_pas.size() &&
is_invalid_offset(p.m_pas[m_idx])) {
++m_idx;
}
}
public:
iterator(passive2& p, unsigned i): p(p), m_idx(i) { fwd(); }
offset_t pas() const { return p.m_pas[m_idx]; }
offset_t sos() const { return (p.hb.vec(pas()).weight().is_pos()?p.m_neg_sos:p.m_pos_sos)[p.m_psos[m_idx]]; }
iterator& operator++() { ++m_idx; fwd(); return *this; }
iterator operator++(int) { iterator tmp = *this; ++*this; return tmp; }
bool operator==(iterator const& it) const {return m_idx == it.m_idx; }
bool operator!=(iterator const& it) const {return m_idx != it.m_idx; }
};
iterator begin() {
return iterator(*this, 0);
}
iterator end() {
return iterator(*this, m_pas.size());
}
private:
void next_resolvable(bool is_positive, unsigned v) {
offset_t pas = m_pas[v];
svector const& soss = is_positive?m_neg_sos:m_pos_sos;
while (m_psos[v] < soss.size()) {
unsigned psos = m_psos[v];
offset_t sos = soss[psos];
if (hb.can_resolve(sos, pas, false)) {
m_weight[v] = m_sum_abs[v] + (is_positive?m_neg_sos_sum[psos]:m_pos_sos_sum[psos]);
m_heap.insert(v);
return;
}
++m_psos[v];
}
// add pas to free-list for hb if it is not in sos.
m_free_list.push_back(v);
m_psos[v] = UINT_MAX;
m_pas[v] = mk_invalid_offset();
}
};
hilbert_basis::hilbert_basis(reslimit& lim):
m_limit(lim),
m_use_support(true),
m_use_ordered_support(true),
m_use_ordered_subsumption(true)
{
m_index = alloc(index, *this);
m_passive = alloc(passive, *this);
m_passive2 = alloc(passive2, *this);
}
hilbert_basis::~hilbert_basis() {
dealloc(m_index);
dealloc(m_passive);
dealloc(m_passive2);
}
hilbert_basis::offset_t hilbert_basis::mk_invalid_offset() {
return offset_t(UINT_MAX);
}
bool hilbert_basis::is_invalid_offset(offset_t offs) {
return offs.m_offset == UINT_MAX;
}
void hilbert_basis::reset() {
m_ineqs.reset();
m_iseq.reset();
m_store.reset();
m_basis.reset();
m_free_list.reset();
m_sos.reset();
m_zero.reset();
m_active.reset();
if (m_passive) {
m_passive->reset();
}
if (m_passive2) {
m_passive2->reset();
}
if (m_index) {
m_index->reset(1);
}
m_ints.reset();
m_current_ineq = 0;
}
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);
}
void hilbert_basis::reset_statistics() {
m_stats.reset();
m_index->reset_statistics();
}
void hilbert_basis::add_ge(rational_vector const& v, rational const& b) {
SASSERT(m_ineqs.empty() || v.size() + 1 == m_ineqs.back().size());
num_vector w;
w.push_back(to_numeral(-b));
for (unsigned i = 0; i < v.size(); ++i) {
w.push_back(to_numeral(v[i]));
}
m_ineqs.push_back(w);
m_iseq.push_back(false);
}
void hilbert_basis::add_le(rational_vector const& v, rational const& b) {
rational_vector w(v);
for (unsigned i = 0; i < w.size(); ++i) {
w[i].neg();
}
add_ge(w, -b);
}
void hilbert_basis::add_eq(rational_vector const& v, rational const& b) {
SASSERT(m_ineqs.empty() || v.size() + 1 == m_ineqs.back().size());
num_vector w;
w.push_back(to_numeral(-b));
for (unsigned i = 0; i < v.size(); ++i) {
w.push_back(to_numeral(v[i]));
}
m_ineqs.push_back(w);
m_iseq.push_back(true);
}
void hilbert_basis::add_ge(rational_vector const& v) {
add_ge(v, rational(0));
}
void hilbert_basis::add_le(rational_vector const& v) {
add_le(v, rational(0));
}
void hilbert_basis::add_eq(rational_vector const& v) {
add_eq(v, rational(0));
}
void hilbert_basis::set_is_int(unsigned var_index) {
//
// The 0'th index is reserved for the constant
// coefficient. Shift indices by 1.
//
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;
}
else {
SASSERT(m_ineqs.back().size() > 1);
return m_ineqs.back().size();
}
}
hilbert_basis::values hilbert_basis::vec(offset_t offs) const {
return values(m_ineqs.size(), m_store.data() + offs.m_offset);
}
void hilbert_basis::init_basis() {
m_basis.reset();
m_store.reset();
m_free_list.reset();
unsigned nv = get_num_vars();
for (unsigned i = 0; i < nv; ++i) {
add_unit_vector(i, numeral(1));
}
for (unsigned i = 0; i < m_ints.size(); ++i) {
add_unit_vector(m_ints[i], numeral(-1));
}
}
void hilbert_basis::add_unit_vector(unsigned i, numeral const& e) {
unsigned num_vars = get_num_vars();
num_vector w(num_vars, numeral(0));
w[i] = e;
offset_t idx = alloc_vector();
values v = vec(idx);
for (unsigned j = 0; j < num_vars; ++j) {
v[j] = w[j];
}
m_basis.push_back(idx);
}
lbool hilbert_basis::saturate() {
init_basis();
m_current_ineq = 0;
while (checkpoint() && m_current_ineq < m_ineqs.size()) {
select_inequality();
stopwatch sw;
sw.start();
lbool r = saturate(m_ineqs[m_current_ineq], m_iseq[m_current_ineq]);
IF_VERBOSE(3,
{ statistics st;
collect_statistics(st);
st.display(verbose_stream());
sw.stop();
verbose_stream() << "time: " << sw.get_seconds() << "\n";
});
++m_stats.m_num_saturations;
if (r != l_true) {
return r;
}
++m_current_ineq;
}
if (!checkpoint()) {
return l_undef;
}
return l_true;
}
lbool hilbert_basis::saturate_orig(num_vector const& ineq, bool is_eq) {
m_active.reset();
m_passive->reset();
m_zero.reset();
m_index->reset(m_current_ineq+1);
int_table support;
TRACE("hilbert_basis", display_ineq(tout, ineq, is_eq););
iterator it = begin();
for (; it != end(); ++it) {
offset_t idx = *it;
values v = vec(idx);
v.weight() = get_weight(v, ineq);
for (unsigned k = 0; k < m_current_ineq; ++k) {
v.weight(k) = get_weight(v, m_ineqs[k]);
}
add_goal(idx);
if (m_use_support) {
support.insert(idx.m_offset);
}
}
TRACE("hilbert_basis", display(tout););
// resolve passive into active
offset_t j = alloc_vector();
while (!m_passive->empty()) {
if (!checkpoint()) {
return l_undef;
}
offset_t idx = m_passive->pop();
TRACE("hilbert_basis", display(tout););
if (is_subsumed(idx)) {
recycle(idx);
continue;
}
for (unsigned i = 0; checkpoint() && i < m_active.size(); ++i) {
if ((!m_use_support || support.contains(m_active[i].m_offset)) && can_resolve(idx, m_active[i], true)) {
resolve(idx, m_active[i], j);
if (add_goal(j)) {
j = alloc_vector();
}
}
}
m_active.push_back(idx);
}
m_free_list.push_back(j);
// Move positive from active and zeros to new basis.
m_basis.reset();
m_basis.append(m_zero);
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);
}
else {
m_free_list.push_back(idx);
}
}
m_active.reset();
m_passive->reset();
m_zero.reset();
TRACE("hilbert_basis", display(tout););
return m_basis.empty()?l_false:l_true;
}
bool hilbert_basis::vector_lt(offset_t idx1, offset_t idx2) const {
values v = vec(idx1);
values w = vec(idx2);
numeral a(0), b(0);
for (unsigned i = 0; i < get_num_vars(); ++i) {
a += abs(v[i]);
b += abs(w[i]);
}
return a < b;
}
lbool hilbert_basis::saturate(num_vector const& ineq, bool is_eq) {
m_zero.reset();
m_index->reset(m_current_ineq+1);
m_passive2->reset();
m_sos.reset();
TRACE("hilbert_basis", display_ineq(tout, ineq, is_eq););
unsigned init_basis_size = 0;
for (unsigned i = 0; i < m_basis.size(); ++i) {
offset_t idx = m_basis[i];
values v = vec(idx);
v.weight() = get_weight(v, ineq);
for (unsigned k = 0; k < m_current_ineq; ++k) {
v.weight(k) = get_weight(v, m_ineqs[k]);
}
m_index->insert(idx, v);
if (v.weight().is_zero()) {
m_zero.push_back(idx);
}
else {
if (v.weight().is_pos()) {
m_basis[init_basis_size++] = idx;
}
m_sos.push_back(idx);
}
}
m_basis.resize(init_basis_size);
m_passive2->init(m_sos);
// ASSERT basis is sorted by weight.
// initialize passive
for (unsigned i = 0; (init_basis_size > 0) && i < m_sos.size(); ++i) {
if (vec(m_sos[i]).weight().is_neg()) {
m_passive2->insert(m_sos[i], 0);
}
}
TRACE("hilbert_basis", display(tout););
// resolve passive into active
offset_t idx = alloc_vector();
while (checkpoint() && !m_passive2->empty()) {
offset_t sos, pas;
TRACE("hilbert_basis", display(tout); );
unsigned offset = m_passive2->pop(sos, pas);
SASSERT(can_resolve(sos, pas, true));
resolve(sos, pas, idx);
if (is_subsumed(idx)) {
continue;
}
values v = vec(idx);
m_index->insert(idx, v);
if (v.weight().is_zero()) {
m_zero.push_back(idx);
}
else {
if (!m_use_ordered_support) {
offset = 0;
}
m_passive2->insert(idx, offset);
if (v.weight().is_pos()) {
m_basis.push_back(idx);
}
}
idx = alloc_vector();
}
if (!checkpoint()) {
return l_undef;
}
m_free_list.push_back(idx);
// remove positive values from basis if we are looking for an equality.
while (is_eq && !m_basis.empty()) {
m_free_list.push_back(m_basis.back());
m_basis.pop_back();
}
m_basis.append(m_zero);
std::sort(m_basis.begin(), m_basis.end(), vector_lt_t(*this));
m_zero.reset();
TRACE("hilbert_basis", display(tout););
return m_basis.empty()?l_false:l_true;
}
void hilbert_basis::get_basis_solution(unsigned i, rational_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(to_rational(vec(offs)[i]));
}
is_initial = !vec(offs)[0].is_zero();
}
void hilbert_basis::get_ge(unsigned i, rational_vector& v, rational& b, bool& is_eq) {
v.reset();
for (unsigned j = 1; j < m_ineqs[i].size(); ++j) {
v.push_back(to_rational(m_ineqs[i][j]));
}
b = to_rational(-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 == 0) {
prod = prod2;
non_zeros = non_zeros2;
best = j;
break;
}
if (non_zeros2 < non_zeros || (non_zeros2 == non_zeros && prod2 < prod)) {
prod = prod2;
non_zeros = non_zeros2;
best = j;
}
}
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]);
}
}
unsigned hilbert_basis::get_num_nonzeros(num_vector const& ineq) {
unsigned count = 0;
for (unsigned i = 0; i < ineq.size(); ++i) {
if (!ineq[i].is_zero()) {
++count;
}
}
return count;
}
unsigned hilbert_basis::get_ineq_product(num_vector const& ineq) {
unsigned num_pos = 0, num_neg = 0;
iterator it = begin();
for (; it != end(); ++it) {
values v = vec(*it);
numeral w = get_weight(v, ineq);
if (w.is_pos()) {
++num_pos;
}
else if (w.is_neg()) {
++num_neg;
}
}
return num_pos * num_neg;
}
hilbert_basis::numeral hilbert_basis::get_ineq_diff(num_vector const& ineq) {
numeral max_pos(0), min_neg(0);
iterator it = begin();
for (; it != end(); ++it) {
values v = vec(*it);
numeral w = get_weight(v, ineq);
if (w > max_pos) {
max_pos = w;
}
else if (w < min_neg) {
min_neg = w;
}
}
return max_pos - min_neg;
}
void hilbert_basis::recycle(offset_t idx) {
m_index->remove(idx, vec(idx));
m_free_list.push_back(idx);
}
void hilbert_basis::resolve(offset_t i, offset_t j, offset_t r) {
++m_stats.m_num_resolves;
values v = vec(i);
values w = vec(j);
values u = vec(r);
unsigned nv = get_num_vars();
for (unsigned k = 0; k < nv; ++k) {
u[k] = v[k] + w[k];
}
u.weight() = v.weight() + w.weight();
for (unsigned k = 0; k < m_current_ineq; ++k) {
u.weight(k) = v.weight(k) + w.weight(k);
}
TRACE("hilbert_basis_verbose",
display(tout, i);
display(tout, j);
display(tout, r);
);
}
hilbert_basis::offset_t hilbert_basis::alloc_vector() {
if (m_free_list.empty()) {
unsigned sz = m_ineqs.size() + get_num_vars();
unsigned idx = m_store.size();
m_store.resize(idx + sz);
// std::cout << "alloc vector: " << idx << " " << sz << " " << m_store.c_ptr() + idx << " " << m_ineqs.size() << "\n";
return offset_t(idx);
}
else {
offset_t result = m_free_list.back();
m_free_list.pop_back();
return result;
}
}
bool hilbert_basis::checkpoint() {
return m_limit.inc();
}
bool hilbert_basis::add_goal(offset_t idx) {
TRACE("hilbert_basis", display(tout, idx););
values v = vec(idx);
if (is_subsumed(idx)) {
return false;
}
m_index->insert(idx, v);
if (v.weight().is_zero()) {
m_zero.push_back(idx);
}
else {
m_passive->insert(idx);
}
return true;
}
bool hilbert_basis::is_subsumed(offset_t idx) {
if (m_index->find(idx, vec(idx))) {
++m_stats.m_num_subsumptions;
return true;
}
return false;
}
bool hilbert_basis::can_resolve(offset_t i, offset_t j, bool check_sign) const {
if (check_sign && get_sign(i) == get_sign(j)) {
return false;
}
SASSERT(get_sign(i) != get_sign(j));
values const& v1 = vec(i);
values const& v2 = vec(j);
if (v1[0].is_one() && v2[0].is_one()) {
return false;
}
for (unsigned i = 0; i < m_ints.size(); ++i) {
unsigned j = m_ints[i];
if (v1[j].is_pos() && v2[j].is_neg()) {
return false;
}
if (v1[j].is_neg() && v2[j].is_pos()) {
return false;
}
}
return true;
}
hilbert_basis::sign_t hilbert_basis::get_sign(offset_t idx) const {
numeral const& val = vec(idx).weight();
if (val.is_pos()) {
return pos;
}
if (val.is_neg()) {
return neg;
}
return zero;
}
hilbert_basis::numeral hilbert_basis::get_weight(values const & 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];
}
return result;
}
void hilbert_basis::display(std::ostream& out) const {
out << "inequalities:\n";
for (unsigned i = 0; i < m_ineqs.size(); ++i) {
display_ineq(out, m_ineqs[i], m_iseq[i]);
}
if (!m_basis.empty()) {
out << "basis:\n";
for (iterator it = begin(); it != end(); ++it) {
display(out, *it);
}
}
if (!m_active.empty()) {
out << "active:\n";
for (unsigned i = 0; i < m_active.size(); ++i) {
display(out, m_active[i]);
}
}
if (!m_passive->empty()) {
passive::iterator it = m_passive->begin();
passive::iterator end = m_passive->end();
out << "passive:\n";
for (; it != end; ++it) {
display(out, *it);
}
}
if (!m_passive2->empty()) {
passive2::iterator it = m_passive2->begin();
passive2::iterator end = m_passive2->end();
out << "passive:\n";
for (; it != end; ++it) {
display(out << "sos:", it.sos());
display(out << "pas:", it.pas());
}
}
if (!m_zero.empty()) {
out << "zero:\n";
for (unsigned i = 0; i < m_zero.size(); ++i) {
display(out, m_zero[i]);
}
}
if (m_index) {
m_index->display(out);
}
}
void hilbert_basis::display(std::ostream& out, offset_t o) const {
display(out, vec(o));
out << " -> " << vec(o).weight() << "\n";
}
void hilbert_basis::display(std::ostream& out, values const& v) const {
unsigned nv = get_num_vars();
for (unsigned j = 0; j < nv; ++j) {
out << v[j] << " ";
}
}
void hilbert_basis::display_ineq(std::ostream& out, num_vector const& v, bool is_eq) const {
unsigned nv = v.size();
for (unsigned j = 1; j < nv; ++j) {
if (!v[j].is_zero()) {
if (j > 0) {
if (v[j].is_pos()) {
out << " + ";
}
else {
out << " - ";
}
}
else if (j == 0 && v[0].is_neg()) {
out << "-";
}
if (!v[j].is_one() && !v[j].is_minus_one()) {
out << abs(v[j]) << "*";
}
out << "x" << j;
}
}
if (is_eq) {
out << " = " << -v[0] << "\n";
}
else {
out << " >= " << -v[0] << "\n";
}
}
/**
Vector v is subsumed by vector w if
v[i] >= w[i] for each index i.
a*v >= a*w for the evaluation of vectors with respect to a.
. a*v < 0 => a*v = a*w
. a*v > 0 => a*w > 0
. a*v = 0 => a*w = 0
Justification:
let u := v - w, then
u[i] >= 0 for each index i
a*u = a*(v-w) >= 0
So v = u + w, where a*u >= 0, a*w >= 0.
If a*v >= a*w >= 0 then v and w are linear
solutions of e_i, and also v-w is a solution.
If a*v = a*w < 0, then a*(v-w) = 0, so v can be obtained from w + (v - w).
*/
bool hilbert_basis::is_subsumed(offset_t i, offset_t j) const {
values v = vec(i);
values w = vec(j);
numeral const& n = v.weight();
numeral const& m = w.weight();
bool r =
i.m_offset != j.m_offset &&
n >= m && (!m.is_neg() || n == m) &&
is_geq(v, w);
for (unsigned k = 0; r && k < m_current_ineq; ++k) {
r = v.weight(k) >= w.weight(k);
}
CTRACE("hilbert_basis", r,
display(tout, i);
tout << " <= \n";
display(tout, j);
tout << "\n";);
return r;
}
bool hilbert_basis::is_geq(values const& v, values const& w) const {
unsigned nv = get_num_vars();
for (unsigned i = 0; i < nv; ++i) {
if (!is_abs_geq(v[i], w[i])) {
return false;
}
}
return true;
}
bool hilbert_basis::is_abs_geq(numeral const& v, numeral const& w) const {
if (w.is_neg()) {
return v <= w;
}
else {
return v >= w;
}
}