mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-28 11:25:51 +00:00
Code Activity

updates to sorting networks

Browse source 
Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
This commit is contained in:
Nikolaj Bjorner 2017-09-23 22:36:19 -05:00
parent 3c4ac9aee5
commit edb3569599
18 changed files with 2070 additions and 170 deletions
Download patch file Download diff file Expand all files Collapse all files

35
src/ast/rewriter/pb2bv_rewriter.cpp
View file

@ -43,6 +43,7 @@ struct pb2bv_rewriter::imp {
struct card2bv_rewriter {
typedef expr* literal;
typedef ptr_vector<expr> literal_vector;
sorting_network_config m_cfg;
psort_nw<card2bv_rewriter> m_sort;
ast_manager& m;
imp& m_imp;
@ -570,8 +571,8 @@ struct pb2bv_rewriter::imp {
public:
card2bv_rewriter(imp& i, ast_manager& m):
m_sort(*this),
card2bv_rewriter(imp& i, ast_manager& m):
m_sort(*this, m_cfg),
m(m),
m_imp(i),
au(m),
@ -760,8 +761,8 @@ struct pb2bv_rewriter::imp {
m_trail.push_back(l);
return l;
}
literal fresh() {
expr_ref fr(m.mk_fresh_const("sn", m.mk_bool_sort()), m);
literal fresh(char const* n) {
expr_ref fr(m.mk_fresh_const(n, m.mk_bool_sort()), m);
m_imp.m_fresh.push_back(to_app(fr)->get_decl());
return trail(fr);
}
@ -785,6 +786,8 @@ struct pb2bv_rewriter::imp {
void pb_totalizer(bool f) {
m_pb_totalizer = f;
}
void set_at_most1(sorting_network_encoding enc) { m_cfg.m_encoding = enc; }
};
struct card2bv_rewriter_cfg : public default_rewriter_cfg {
@ -800,6 +803,8 @@ struct pb2bv_rewriter::imp {
void keep_pb_constraints(bool f) { m_r.keep_pb_constraints(f); }
void pb_num_system(bool f) { m_r.pb_num_system(f); }
void pb_totalizer(bool f) { m_r.pb_totalizer(f); }
void set_at_most1(sorting_network_encoding enc) { m_r.set_at_most1(enc); }
};
class card_pb_rewriter : public rewriter_tpl<card2bv_rewriter_cfg> {
@ -812,6 +817,7 @@ struct pb2bv_rewriter::imp {
void keep_pb_constraints(bool f) { m_cfg.keep_pb_constraints(f); }
void pb_num_system(bool f) { m_cfg.pb_num_system(f); }
void pb_totalizer(bool f) { m_cfg.pb_totalizer(f); }
void set_at_most1(sorting_network_encoding e) { m_cfg.set_at_most1(e); }
};
card_pb_rewriter m_rw;
@ -844,15 +850,25 @@ struct pb2bv_rewriter::imp {
gparams::get_module("sat").get_sym("pb.solver", symbol()) == symbol("totalizer");
}
sorting_network_encoding atmost1_encoding() const {
symbol enc = m_params.get_sym("atmost1_encoding", enc);
if (enc == symbol()) {
enc = gparams::get_module("sat").get_sym("atmost1_encoding", symbol());
}
if (enc == symbol("grouped")) return sorting_network_encoding::grouped_at_most_1;
if (enc == symbol("bimander")) return sorting_network_encoding::bimander_at_most_1;
if (enc == symbol("ordered")) return sorting_network_encoding::ordered_at_most_1;
return grouped_at_most_1;
}
imp(ast_manager& m, params_ref const& p):
m(m), m_params(p), m_lemmas(m),
m_fresh(m),
m_num_translated(0),
m_rw(*this, m) {
m_rw.keep_cardinality_constraints(keep_cardinality());
m_rw.keep_pb_constraints(keep_pb());
m_rw.pb_num_system(pb_num_system());
m_rw.pb_totalizer(pb_totalizer());
updt_params(p);
}
void updt_params(params_ref const & p) {
@ -860,7 +876,8 @@ struct pb2bv_rewriter::imp {
m_rw.keep_cardinality_constraints(keep_cardinality());
m_rw.keep_pb_constraints(keep_pb());
m_rw.pb_num_system(pb_num_system());
m_rw.pb_totalizer(pb_totalizer());
m_rw.pb_totalizer(pb_totalizer());
m_rw.set_at_most1(atmost1_encoding());
}
void collect_param_descrs(param_descrs& r) const {
r.insert("keep_cardinality_constraints", CPK_BOOL, "(default: true) retain cardinality constraints (don't bit-blast them) and use built-in cardinality solver");

7
src/opt/sortmax.cpp
View file

@ -32,12 +32,13 @@ namespace opt {
public:
typedef expr* literal;
typedef ptr_vector<expr> literal_vector;
sorting_network_config m_cfg;
psort_nw<sortmax> m_sort;
expr_ref_vector m_trail;
func_decl_ref_vector m_fresh;
ref<filter_model_converter> m_filter;
sortmax(maxsat_context& c, weights_t& ws, expr_ref_vector const& soft):
maxsmt_solver_base(c, ws, soft), m_sort(*this), m_trail(m), m_fresh(m) {}
maxsmt_solver_base(c, ws, soft), m_sort(*this, m_cfg), m_trail(m), m_fresh(m) {}
virtual ~sortmax() {}
@ -138,8 +139,8 @@ namespace opt {
m_trail.push_back(l);
return l;
}
literal fresh() {
expr_ref fr(m.mk_fresh_const("sn", m.mk_bool_sort()), m);
literal fresh(char const* n) {
expr_ref fr(m.mk_fresh_const(n, m.mk_bool_sort()), m);
func_decl* f = to_app(fr)->get_decl();
m_fresh.push_back(f);
m_filter->insert(f);

3
src/sat/sat_config.cpp
View file

@ -87,6 +87,7 @@ namespace sat {
m_local_search_threads = p.local_search_threads();
m_lookahead_simplify = p.lookahead_simplify();
m_lookahead_search = p.lookahead_search();
m_lookahead_reward = p.lookahead_reward();
m_ccc = p.ccc();
// These parameters are not exposed
@ -163,7 +164,7 @@ namespace sat {
m_pb_solver = PB_SOLVER;
}
else {
throw sat_param_exception("invalid PB solver");
throw sat_param_exception("invalid PB solver: solver, totalizer, circuit, sorting");
}
}

1
src/sat/sat_config.h
View file

@ -75,6 +75,7 @@ namespace sat {
bool m_local_search;
bool m_lookahead_search;
bool m_lookahead_simplify;
symbol m_lookahead_reward;
bool m_ccc;
unsigned m_simplify_mult1;

228
src/sat/sat_lookahead.cpp
View file

@ -41,7 +41,6 @@ namespace sat {
}
}
void lookahead::flip_prefix() {
if (m_trail_lim.size() < 64) {
uint64 mask = (1ull << m_trail_lim.size());
@ -228,11 +227,6 @@ namespace sat {
if (is_sat()) {
return false;
}
if (newbies) {
enable_trace("sat");
TRACE("sat", display(tout););
TRACE("sat", tout << sum << "\n";);
}
}
SASSERT(!m_candidates.empty());
// cut number of candidates down to max_num_cand.
@ -292,9 +286,8 @@ namespace sat {
double lookahead::init_candidates(unsigned level, bool newbies) {
m_candidates.reset();
double sum = 0;
for (bool_var const* it = m_freevars.begin(), * end = m_freevars.end(); it != end; ++it) {
SASSERT(is_undef(*it));
bool_var x = *it;
for (bool_var x : m_freevars) {
SASSERT(is_undef(x));
if (!m_select_lookahead_vars.empty()) {
if (m_select_lookahead_vars.contains(x)) {
m_candidates.push_back(candidate(x, m_rating[x]));
@ -333,20 +326,20 @@ namespace sat {
}
bool lookahead::is_sat() const {
for (bool_var const* it = m_freevars.begin(), * end = m_freevars.end(); it != end; ++it) {
literal l(*it, false);
for (bool_var x : m_freevars) {
literal l(x, false);
literal_vector const& lits1 = m_binary[l.index()];
for (unsigned i = 0; i < lits1.size(); ++i) {
if (!is_true(lits1[i])) {
TRACE("sat", tout << l << " " << lits1[i] << "\n";);
for (literal lit1 : lits1) {
if (!is_true(lit1)) {
TRACE("sat", tout << l << " " << lit1 << "\n";);
return false;
}
}
l.neg();
literal_vector const& lits2 = m_binary[l.index()];
for (unsigned i = 0; i < lits2.size(); ++i) {
if (!is_true(lits2[i])) {
TRACE("sat", tout << l << " " << lits2[i] << "\n";);
for (literal lit2 : lits2) {
if (!is_true(lit2)) {
TRACE("sat", tout << l << " " << lit2 << "\n";);
return false;
}
}
@ -389,12 +382,114 @@ namespace sat {
break;
}
case heule_schur_reward:
heule_schur_scores();
break;
case heule_unit_reward:
heule_unit_scores();
break;
case unit_literal_reward:
heule_schur_scores();
break;
}
}
static unsigned counter = 0;
void lookahead::heule_schur_scores() {
if (counter % 10 != 0) return;
++counter;
for (bool_var x : m_freevars) {
literal l(x, false);
m_rating[l.var()] = heule_schur_score(l) * heule_schur_score(~l);
}
}
double lookahead::heule_schur_score(literal l) {
double sum = 0;
for (literal lit : m_binary[l.index()]) {
if (is_undef(lit)) sum += literal_occs(lit) / 4.0;
}
watch_list& wlist = m_watches[l.index()];
for (auto & w : wlist) {
switch (w.get_kind()) {
case watched::BINARY:
UNREACHABLE();
break;
case watched::TERNARY: {
literal l1 = w.get_literal1();
literal l2 = w.get_literal2();
if (is_undef(l1) && is_undef(l2)) {
sum += (literal_occs(l1) + literal_occs(l2)) / 8.0;
}
break;
}
case watched::CLAUSE: {
clause_offset cls_off = w.get_clause_offset();
clause & c = *(m_cls_allocator.get_clause(cls_off));
unsigned sz = 0;
double to_add = 0;
for (literal lit : c) {
if (is_undef(lit) && lit != ~l) {
to_add += literal_occs(lit);
++sz;
}
}
sum += pow(0.5, sz) * to_add / sz;
break;
}
default:
break;
}
}
return sum;
}
void lookahead::heule_unit_scores() {
if (counter % 10 != 0) return;
++counter;
for (bool_var x : m_freevars) {
literal l(x, false);
m_rating[l.var()] = heule_unit_score(l) * heule_unit_score(~l);
}
}
double lookahead::heule_unit_score(literal l) {
double sum = 0;
for (literal lit : m_binary[l.index()]) {
if (is_undef(lit)) sum += 0.25;
}
watch_list& wlist = m_watches[l.index()];
for (auto & w : wlist) {
switch (w.get_kind()) {
case watched::BINARY:
UNREACHABLE();
break;
case watched::TERNARY: {
literal l1 = w.get_literal1();
literal l2 = w.get_literal2();
if (is_undef(l1) && is_undef(l2)) {
sum += 0.25;
}
break;
}
case watched::CLAUSE: {
clause_offset cls_off = w.get_clause_offset();
clause & c = *(m_cls_allocator.get_clause(cls_off));
unsigned sz = 0;
for (literal lit : c) {
if (is_undef(lit) && lit != ~l) {
++sz;
}
}
sum += pow(0.5, sz);
break;
}
default:
break;
}
}
return sum;
}
void lookahead::ensure_H(unsigned level) {
while (m_H.size() <= level) {
m_H.push_back(svector<double>());
@ -404,16 +499,16 @@ namespace sat {
void lookahead::h_scores(svector<double>& h, svector<double>& hp) {
double sum = 0;
for (bool_var const* it = m_freevars.begin(), * end = m_freevars.end(); it != end; ++it) {
literal l(*it, false);
for (bool_var x : m_freevars) {
literal l(x, false);
sum += h[l.index()] + h[(~l).index()];
}
if (sum == 0) sum = 0.0001;
double factor = 2 * m_freevars.size() / sum;
double sqfactor = factor * factor;
double afactor = factor * m_config.m_alpha;
for (bool_var const* it = m_freevars.begin(), * end = m_freevars.end(); it != end; ++it) {
literal l(*it, false);
for (bool_var x : m_freevars) {
literal l(x, false);
double pos = l_score(l, h, factor, sqfactor, afactor);
double neg = l_score(~l, h, factor, sqfactor, afactor);
hp[l.index()] = pos;
@ -425,28 +520,25 @@ namespace sat {
double lookahead::l_score(literal l, svector<double> const& h, double factor, double sqfactor, double afactor) {
double sum = 0, tsum = 0;
literal_vector::iterator it = m_binary[l.index()].begin(), end = m_binary[l.index()].end();
for (; it != end; ++it) {
bool_var v = it->var();
if (is_undef(*it)) sum += h[it->index()];
// if (m_freevars.contains(it->var())) sum += h[it->index()];
for (literal lit : m_binary[l.index()]) {
if (is_undef(lit)) sum += h[lit.index()];
// if (m_freevars.contains(lit.var())) sum += h[lit.index()];
}
watch_list& wlist = m_watches[l.index()];
watch_list::iterator wit = wlist.begin(), wend = wlist.end();
for (; wit != wend; ++wit) {
switch (wit->get_kind()) {
for (auto & w : wlist) {
switch (w.get_kind()) {
case watched::BINARY:
UNREACHABLE();
break;
case watched::TERNARY: {
literal l1 = wit->get_literal1();
literal l2 = wit->get_literal2();
literal l1 = w.get_literal1();
literal l2 = w.get_literal2();
// if (is_undef(l1) && is_undef(l2))
tsum += h[l1.index()] * h[l2.index()];
break;
}
case watched::CLAUSE: {
clause_offset cls_off = wit->get_clause_offset();
clause_offset cls_off = w.get_clause_offset();
clause & c = *(m_cls_allocator.get_clause(cls_off));
// approximation compared to ternary clause case:
// we pick two other literals from the clause.
@ -865,8 +957,6 @@ namespace sat {
copy_clauses(m_s.m_clauses);
copy_clauses(m_s.m_learned);
// m_config.m_use_ternary_reward &= !m_s.m_ext;
// copy units
unsigned trail_sz = m_s.init_trail_size();
for (unsigned i = 0; i < trail_sz; ++i) {
@ -995,15 +1085,17 @@ namespace sat {
return unsat;
}
void lookahead::push_lookahead1(literal lit, unsigned level) {
unsigned lookahead::push_lookahead1(literal lit, unsigned level) {
SASSERT(m_search_mode == lookahead_mode::searching);
m_search_mode = lookahead_mode::lookahead1;
scoped_level _sl(*this, level);
unsigned old_sz = m_trail.size();
assign(lit);
propagate();
return m_trail.size() - old_sz;
}
void lookahead::pop_lookahead1(literal lit) {
void lookahead::pop_lookahead1(literal lit, unsigned num_units) {
bool unsat = inconsistent();
SASSERT(m_search_mode == lookahead_mode::lookahead1);
m_inconsistent = false;
@ -1025,8 +1117,15 @@ namespace sat {
}
m_stats.m_windfall_binaries += m_wstack.size();
}
if (m_config.m_reward_type == unit_literal_reward) {
m_lookahead_reward += m_wstack.size();
switch (m_config.m_reward_type) {
case unit_literal_reward:
m_lookahead_reward += num_units;
break;
case heule_unit_reward:
case heule_schur_reward:
break;
default:
break;
}
m_wstack.reset();
}
@ -1226,6 +1325,9 @@ namespace sat {
case heule_schur_reward:
m_lookahead_reward += (literal_occs(l1) + literal_occs(l2)) / 8.0;
break;
case heule_unit_reward:
m_lookahead_reward += 0.25;
break;
case unit_literal_reward:
break;
}
@ -1253,6 +1355,9 @@ namespace sat {
m_lookahead_reward += pow(0.5, sz) * to_add / sz;
break;
}
case heule_unit_reward:
m_lookahead_reward += pow(0.5, sz);
break;
case ternary_reward:
m_lookahead_reward = (double)0.001;
break;
@ -1326,11 +1431,13 @@ namespace sat {
IF_VERBOSE(30, verbose_stream() << scope_lvl() << " " << lit << " binary: " << m_binary_trail.size() << " trail: " << m_trail_lim.back() << "\n";);
}
TRACE("sat", tout << "lookahead: " << lit << " @ " << m_lookahead[i].m_offset << "\n";);
unsigned old_trail_sz = m_trail.size();
reset_lookahead_reward(lit);
push_lookahead1(lit, level);
if (!first) do_double(lit, base);
bool unsat = inconsistent();
pop_lookahead1(lit);
bool unsat = inconsistent();
unsigned num_units = m_trail.size() - old_trail_sz;
pop_lookahead1(lit, num_units);
if (unsat) {
TRACE("sat", tout << "backtracking and settting " << ~lit << "\n";);
reset_lookahead_reward();
@ -1407,6 +1514,7 @@ namespace sat {
switch (m_config.m_reward_type) {
case ternary_reward: return l + r + (1 << 10) * l * r;
case heule_schur_reward: return l * r;
case heule_unit_reward: return l * r;
case unit_literal_reward: return l * r;
default: UNREACHABLE(); return l * r;
}
@ -1489,7 +1597,7 @@ namespace sat {
}
}
}
else {
else {
inc_lookahead_reward(l, m_lookahead_reward);
}
}
@ -1625,7 +1733,13 @@ namespace sat {
}
TRACE("sat", tout << "choose: " << l << " " << trail << "\n";);
++m_stats.m_decisions;
IF_VERBOSE(1, verbose_stream() << "select " << pp_prefix(m_prefix, m_trail_lim.size()) << ": " << l << " " << m_trail.size() << "\n";);
IF_VERBOSE(1, printf("\r");
std::stringstream strm;
strm << pp_prefix(m_prefix, m_trail_lim.size());
for (unsigned i = 0; i < 50; ++i) strm << " ";
printf(strm.str().c_str());
fflush(stdout);
);
push(l, c_fixed_truth);
trail.push_back(l);
SASSERT(inconsistent() || !is_unsat());
@ -1851,16 +1965,20 @@ namespace sat {
m_lookahead.reset();
}
std::ostream& lookahead::display(std::ostream& out) const {
std::ostream& lookahead::display_summary(std::ostream& out) const {
out << "Prefix: " << pp_prefix(m_prefix, m_trail_lim.size()) << "\n";
out << "Level: " << m_level << "\n";
out << "Free vars: " << m_freevars.size() << "\n";
return out;
}
std::ostream& lookahead::display(std::ostream& out) const {
display_summary(out);
display_values(out);
display_binary(out);
display_clauses(out);
out << "free vars: ";
for (bool_var const* it = m_freevars.begin(), * end = m_freevars.end(); it != end; ++it) {
out << *it << " ";
}
for (bool_var b : m_freevars) out << b << " ";
out << "\n";
for (unsigned i = 0; i < m_watches.size(); ++i) {
watch_list const& wl = m_watches[i];
@ -1879,6 +1997,24 @@ namespace sat {
return m_model;
}
void lookahead::init_config() {
if (m_s.m_config.m_lookahead_reward == symbol("hs")) {
m_config.m_reward_type = heule_schur_reward;
}
else if (m_s.m_config.m_lookahead_reward == symbol("heuleu")) {
m_config.m_reward_type = heule_unit_reward;
}
else if (m_s.m_config.m_lookahead_reward == symbol("ternary")) {
m_config.m_reward_type = ternary_reward;
}
else if (m_s.m_config.m_lookahead_reward == symbol("unit")) {
m_config.m_reward_type = unit_literal_reward;
}
else {
warning_msg("Reward type not recognized");
}
}
void lookahead::collect_statistics(statistics& st) const {
st.update("lh bool var", m_vprefix.size());
st.update("lh clauses", m_clauses.size());

15
src/sat/sat_lookahead.h
View file

@ -69,7 +69,8 @@ namespace sat {
enum reward_t {
ternary_reward,
unit_literal_reward,
heule_schur_reward
heule_schur_reward,
heule_unit_reward
};
struct config {
@ -277,6 +278,10 @@ namespace sat {
void init_pre_selection(unsigned level);
void ensure_H(unsigned level);
void h_scores(svector<double>& h, svector<double>& hp);
void heule_schur_scores();
double heule_schur_score(literal l);
void heule_unit_scores();
double heule_unit_score(literal l);
double l_score(literal l, svector<double> const& h, double factor, double sqfactor, double afactor);
// ------------------------------------
@ -393,8 +398,8 @@ namespace sat {
void push(literal lit, unsigned level);
void pop();
bool push_lookahead2(literal lit, unsigned level);
void push_lookahead1(literal lit, unsigned level);
void pop_lookahead1(literal lit);
unsigned push_lookahead1(literal lit, unsigned level);
void pop_lookahead1(literal lit, unsigned num_units);
double mix_diff(double l, double r) const;
clause const& get_clause(watch_list::iterator it) const;
bool is_nary_propagation(clause const& c, literal l) const;
@ -444,6 +449,8 @@ namespace sat {
void init_search();
void checkpoint();
void init_config();
public:
lookahead(solver& s) :
m_s(s),
@ -453,6 +460,7 @@ namespace sat {
m_level(2),
m_prefix(0) {
m_s.rlimit().push_child(&m_rlimit);
init_config();
}
~lookahead() {
@ -488,6 +496,7 @@ namespace sat {
void scc();
std::ostream& display(std::ostream& out) const;
std::ostream& display_summary(std::ostream& out) const;
model const& get_model();
void collect_statistics(statistics& st) const;

2
src/sat/sat_params.pyg
View file

@ -31,9 +31,11 @@ def_module_params('sat',
('cardinality.solver', BOOL, False, 'use cardinality solver'),
('pb.solver', SYMBOL, 'circuit', 'method for handling Pseudo-Boolean constraints: circuit (arithmetical circuit), sorting (sorting circuit), totalizer (use totalizer encoding), solver (use SMT solver)'),
('xor.solver', BOOL, False, 'use xor solver'),
('atmost1_encoding', SYMBOL, 'grouped', 'encoding used for at-most-1 constraints grouped, bimander, ordered'),
('local_search_threads', UINT, 0, 'number of local search threads to find satisfiable solution'),
('local_search', BOOL, False, 'use local search instead of CDCL'),
('lookahead_search', BOOL, False, 'use lookahead solver'),
('lookahead_simplify', BOOL, False, 'use lookahead solver during simplification'),
('lookahead.reward', SYMBOL, 'heuleu', 'select lookahead heuristic: ternary, hs (Heule Schur), heuleu (Heule Unit), or unit'),
('ccc', BOOL, False, 'use Concurrent Cube and Conquer solver')
))

8
src/smt/theory_pb.cpp
View file

@ -1396,7 +1396,7 @@ namespace smt {
th(th),
pb(m) {}
literal fresh() {
literal fresh(char const* ) {
app_ref y(m);
y = pb.mk_fresh_bool();
return literal(ctx.mk_bool_var(y));
@ -1441,7 +1441,8 @@ namespace smt {
theory_pb_params p;
theory_pb th(ctx.get_manager(), p);
psort_expr ps(ctx, th);
psort_nw<psort_expr> sort(ps);
sorting_network_config cfg;
psort_nw<psort_expr> sort(ps, cfg);
return sort.ge(false, k, n, xs);
}
@ -1577,7 +1578,8 @@ namespace smt {
psort_expr ps(ctx, *this);
psort_nw<psort_expr> sortnw(ps);
sorting_network_config cfg;
psort_nw<psort_expr> sortnw(ps, cfg);
sortnw.m_stats.reset();
if (ctx.get_assignment(thl) == l_true &&

1
src/test/main.cpp
View file

@ -239,7 +239,6 @@ int main(int argc, char ** argv) {
TST(pdr);
TST_ARGV(ddnf);
TST(model_evaluator);
TST_ARGV(lp);
TST(get_consequences);
TST(pb2bv);
TST_ARGV(sat_lookahead);

10
src/test/sat_local_search.cpp
View file

@ -1,8 +1,8 @@
#include "sat_local_search.h"
#include "sat_solver.h"
#include "cancel_eh.h"
#include "scoped_ctrl_c.h"
#include "scoped_timer.h"
#include "sat/sat_local_search.h"
#include "sat/sat_solver.h"
#include "util/cancel_eh.h"
#include "util/scoped_ctrl_c.h"
#include "util/scoped_timer.h"
static bool build_instance(char const * filename, sat::solver& s, sat::local_search& local_search)
{

10
src/test/sat_lookahead.cpp
View file

@ -1,8 +1,8 @@
#include "sat_solver.h"
#include "sat_watched.h"
#include "statistics.h"
#include "sat_lookahead.h"
#include "dimacs.h"
#include "sat/sat_solver.h"
#include "sat/sat_watched.h"
#include "util/statistics.h"
#include "sat/sat_lookahead.h"
#include "sat/dimacs.h"
static void display_model(sat::model const & m) {
for (unsigned i = 1; i < m.size(); i++) {

157
src/test/sorting_network.cpp
View file

@ -6,14 +6,14 @@ Copyright (c) 2015 Microsoft Corporation
#include "util/trace.h"
#include "util/vector.h"
#include "util/sorting_network.h"
#include "ast/ast.h"
#include "ast/ast_pp.h"
#include "ast/reg_decl_plugins.h"
#include "util/sorting_network.h"
#include "smt/smt_kernel.h"
#include "model/model_smt2_pp.h"
#include "smt/params/smt_params.h"
#include "ast/ast_util.h"
#include "model/model_smt2_pp.h"
#include "smt/smt_kernel.h"
#include "smt/params/smt_params.h"
@ -174,16 +174,61 @@ struct ast_ext2 {
std::ostream& pp(std::ostream& out, literal lit) {
return out << mk_pp(lit, m);
}
literal fresh() {
return trail(m.mk_fresh_const("x", m.mk_bool_sort()));
literal fresh(char const* n) {
return trail(m.mk_fresh_const(n, m.mk_bool_sort()));
}
void mk_clause(unsigned n, literal const* lits) {
m_clauses.push_back(mk_or(m, n, lits));
}
};
static void test_eq1(unsigned n, sorting_network_encoding enc) {
//std::cout << "test eq1 " << n << " for encoding: " << enc << "\n";
ast_manager m;
reg_decl_plugins(m);
ast_ext2 ext(m);
expr_ref_vector in(m), out(m);
for (unsigned i = 0; i < n; ++i) {
in.push_back(m.mk_fresh_const("a",m.mk_bool_sort()));
}
smt_params fp;
smt::kernel solver(m, fp);
sorting_network_config cfg;
cfg.m_encoding = enc;
psort_nw<ast_ext2> sn(ext, cfg);
expr_ref result1(m), result2(m);
static void test_sorting_eq(unsigned n, unsigned k) {
// equality:
solver.push();
result1 = sn.eq(true, 1, in.size(), in.c_ptr());
for (expr* cl : ext.m_clauses) {
solver.assert_expr(cl);
}
expr_ref_vector ors(m);
for (unsigned i = 0; i < n; ++i) {
expr_ref_vector ands(m);
for (unsigned j = 0; j < n; ++j) {
ands.push_back(j == i ? in[j].get() : m.mk_not(in[j].get()));
}
ors.push_back(mk_and(ands));
}
result2 = mk_or(ors);
solver.assert_expr(m.mk_not(m.mk_eq(result1, result2)));
//std::cout << ext.m_clauses << "\n";
//std::cout << result1 << "\n";
//std::cout << result2 << "\n";
lbool res = solver.check();
if (res == l_true) {
model_ref model;
solver.get_model(model);
model_smt2_pp(std::cout, m, *model, 0);
TRACE("pb", model_smt2_pp(tout, m, *model, 0););
}
ENSURE(l_false == res);
ext.m_clauses.reset();
}
static void test_sorting_eq(unsigned n, unsigned k, sorting_network_encoding enc) {
ENSURE(k < n);
ast_manager m;
reg_decl_plugins(m);
@ -194,16 +239,19 @@ static void test_sorting_eq(unsigned n, unsigned k) {
}
smt_params fp;
smt::kernel solver(m, fp);
psort_nw<ast_ext2> sn(ext);
sorting_network_config cfg;
cfg.m_encoding = enc;
psort_nw<ast_ext2> sn(ext, cfg);
expr_ref result(m);
// equality:
std::cout << "eq " << k << "\n";
std::cout << "eq " << k << " out of " << n << " for encoding " << enc << "\n";
solver.push();
result = sn.eq(true, k, in.size(), in.c_ptr());
result = sn.eq(false, k, in.size(), in.c_ptr());
std::cout << result << "\n" << ext.m_clauses << "\n";
solver.assert_expr(result);
for (unsigned i = 0; i < ext.m_clauses.size(); ++i) {
solver.assert_expr(ext.m_clauses[i].get());
for (expr* cl : ext.m_clauses) {
solver.assert_expr(cl);
}
lbool res = solver.check();
ENSURE(res == l_true);
@ -232,7 +280,7 @@ static void test_sorting_eq(unsigned n, unsigned k) {
ext.m_clauses.reset();
}
static void test_sorting_le(unsigned n, unsigned k) {
static void test_sorting_le(unsigned n, unsigned k, sorting_network_encoding enc) {
ast_manager m;
reg_decl_plugins(m);
ast_ext2 ext(m);
@ -242,7 +290,9 @@ static void test_sorting_le(unsigned n, unsigned k) {
}
smt_params fp;
smt::kernel solver(m, fp);
psort_nw<ast_ext2> sn(ext);
sorting_network_config cfg;
cfg.m_encoding = enc;
psort_nw<ast_ext2> sn(ext, cfg);
expr_ref result(m);
// B <= k
std::cout << "le " << k << "\n";
@ -279,7 +329,7 @@ static void test_sorting_le(unsigned n, unsigned k) {
}
void test_sorting_ge(unsigned n, unsigned k) {
void test_sorting_ge(unsigned n, unsigned k, sorting_network_encoding enc) {
ast_manager m;
reg_decl_plugins(m);
ast_ext2 ext(m);
@ -289,7 +339,9 @@ void test_sorting_ge(unsigned n, unsigned k) {
}
smt_params fp;
smt::kernel solver(m, fp);
psort_nw<ast_ext2> sn(ext);
sorting_network_config cfg;
cfg.m_encoding = enc;
psort_nw<ast_ext2> sn(ext, cfg);
expr_ref result(m);
// k <= B
std::cout << "ge " << k << "\n";
@ -325,11 +377,11 @@ void test_sorting_ge(unsigned n, unsigned k) {
solver.pop(1);
}
void test_sorting5(unsigned n, unsigned k) {
void test_sorting5(unsigned n, unsigned k, sorting_network_encoding enc) {
std::cout << "n: " << n << " k: " << k << "\n";
test_sorting_le(n, k);
test_sorting_eq(n, k);
test_sorting_ge(n, k);
test_sorting_le(n, k, enc);
test_sorting_eq(n, k, enc);
test_sorting_ge(n, k, enc);
}
expr_ref naive_at_most1(expr_ref_vector const& xs) {
@ -343,7 +395,7 @@ expr_ref naive_at_most1(expr_ref_vector const& xs) {
return mk_and(clauses);
}
void test_at_most_1(unsigned n, bool full) {
void test_at_most_1(unsigned n, bool full, sorting_network_encoding enc) {
ast_manager m;
reg_decl_plugins(m);
expr_ref_vector in(m), out(m);
@ -352,12 +404,17 @@ void test_at_most_1(unsigned n, bool full) {
}
ast_ext2 ext(m);
psort_nw<ast_ext2> sn(ext);
sorting_network_config cfg;
cfg.m_encoding = enc;
psort_nw<ast_ext2> sn(ext, cfg);
expr_ref result1(m), result2(m);
result1 = sn.le(full, 1, in.size(), in.c_ptr());
result2 = naive_at_most1(in);
std::cout << "clauses: " << ext.m_clauses << "\n-----\n";
std::cout << "encoded: " << result1 << "\n";
std::cout << "naive: " << result2 << "\n";
smt_params fp;
smt::kernel solver(m, fp);
@ -369,15 +426,27 @@ void test_at_most_1(unsigned n, bool full) {
solver.assert_expr(m.mk_not(m.mk_eq(result1, result2)));
std::cout << result1 << "\n";
lbool res = solver.check();
if (res == l_true) {
model_ref model;
solver.get_model(model);
model_smt2_pp(std::cout, m, *model, 0);
}
VERIFY(l_false == solver.check());
VERIFY(l_false == res);
solver.pop(1);
}
if (n >= 9) return;
if (n <= 1) return;
for (unsigned i = 0; i < static_cast<unsigned>(1 << n); ++i) {
std::cout << "checking: " << n << ": " << i << "\n";
std::cout << "checking n: " << n << " bits: ";
for (unsigned j = 0; j < n; ++j) {
bool is_true = (i & (1 << j)) != 0;
std::cout << (is_true?"1":"0");
}
std::cout << "\n";
solver.push();
unsigned k = 0;
for (unsigned j = 0; j < n; ++j) {
@ -388,7 +457,6 @@ void test_at_most_1(unsigned n, bool full) {
std::cout << atom << "\n";
if (is_true) ++k;
}
VERIFY(l_false == solver.check());
if (k > 1) {
solver.assert_expr(result1);
}
@ -405,7 +473,7 @@ void test_at_most_1(unsigned n, bool full) {
}
static void test_at_most1() {
static void test_at_most1(sorting_network_encoding enc) {
ast_manager m;
reg_decl_plugins(m);
expr_ref_vector in(m), out(m);
@ -415,28 +483,45 @@ static void test_at_most1() {
in[4] = in[3].get();
ast_ext2 ext(m);
psort_nw<ast_ext2> sn(ext);
sorting_network_config cfg;
cfg.m_encoding = enc;
psort_nw<ast_ext2> sn(ext, cfg);
expr_ref result(m);
result = sn.le(true, 1, in.size(), in.c_ptr());
std::cout << result << "\n";
std::cout << ext.m_clauses << "\n";
}
void tst_sorting_network() {
for (unsigned i = 1; i < 17; ++i) {
test_at_most_1(i, true);
test_at_most_1(i, false);
}
test_at_most1();
test_sorting_eq(11,7);
static void test_sorting5(sorting_network_encoding enc) {
test_sorting_eq(11,7, enc);
for (unsigned n = 3; n < 20; n += 2) {
for (unsigned k = 1; k < n; ++k) {
test_sorting5(n, k);
test_sorting5(n, k, enc);
}
}
}
static void tst_sorting_network(sorting_network_encoding enc) {
for (unsigned i = 1; i < 17; ++i) {
test_at_most_1(i, true, enc);
test_at_most_1(i, false, enc);
}
for (unsigned n = 2; n < 20; ++n) {
std::cout << "verify eq-1 out of " << n << "\n";
test_sorting_eq(n, 1, enc);
test_eq1(n, enc);
}
test_at_most1(enc);
test_sorting5(enc);
}
void tst_sorting_network() {
tst_sorting_network(sorting_network_encoding::ordered_at_most_1);
tst_sorting_network(sorting_network_encoding::grouped_at_most_1);
tst_sorting_network(sorting_network_encoding::bimander_at_most_1);
test_sorting1();
test_sorting2();
test_sorting3();
test_sorting4();
}

27
src/util/lp/bound_propagator.h Normal file
View file

@ -0,0 +1,27 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include "util/lp/lp_settings.h"
namespace lp {
class lar_solver;
class bound_propagator {
std::unordered_map<unsigned, unsigned> m_improved_low_bounds; // these maps map a column index to the corresponding index in ibounds
std::unordered_map<unsigned, unsigned> m_improved_upper_bounds;
lar_solver & m_lar_solver;
public:
vector<implied_bound> m_ibounds;
public:
bound_propagator(lar_solver & ls);
column_type get_column_type(unsigned) const;
const impq & get_low_bound(unsigned) const;
const impq & get_upper_bound(unsigned) const;
void try_add_bound(mpq v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict);
virtual bool bound_is_interesting(unsigned vi,
lp::lconstraint_kind kind,
const rational & bval) {return true;}
unsigned number_of_found_bounds() const { return m_ibounds.size(); }
virtual void consume(mpq const& v, unsigned j) { std::cout << "doh\n"; }
};
}

201
src/util/lp/cut_solver.h Normal file
View file

@ -0,0 +1,201 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Nikolaj Bjorner, Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include "util/trace.h"
#include "util/lp/lp_settings.h"
namespace lp {
template <typename T>
class cut_solver {
struct ineq { // we only have less or equal, which is enough for integral variables
mpq m_bound;
vector<std::pair<T, var_index>> m_term;
ineq(vector<std::pair<T, var_index>>& term, mpq bound):m_bound(bound),m_term(term) {
}
};
vector<ineq> m_ineqs;
enum class lbool {
l_false, // false
l_true, // true
l_undef // undef
};
enum class literal_type {
BOOL,
INEQ,
BOUND
};
struct literal {
literal_type m_tag;
bool m_sign; // true means the pointed inequality is negated, or bound is negated, or boolean value is negated
unsigned m_id;
unsigned m_index_of_ineq; // index into m_ineqs
bool m_bool_val; // used if m_tag is equal to BOOL
mpq m_bound; // used if m_tag is BOUND
literal(bool sign, bool val): m_tag(literal_type::BOOL),
m_bool_val(val){
}
literal(bool sign, unsigned index_of_ineq) : m_tag(literal_type::INEQ), m_index_of_ineq(index_of_ineq) {}
};
bool lhs_is_int(const vector<std::pair<mpq, var_index>> & lhs) const {
for (auto & p : lhs)
if (p.first.is_int() == false) return false;
return true;
}
public:
void add_ineq(vector<std::pair<mpq, var_index>> & lhs, mpq rhs) {
lp_assert(lhs_is_int(lhs));
lp_assert(rhs.is_int());
m_ineqs.push_back(ineq(lhs, rhs));
}
bool m_inconsistent; // tracks if state is consistent
unsigned m_scope_lvl; // tracks the number of case splits
svector<literal> m_trail;
// backtracking state from the SAT solver:
struct scope {
unsigned m_trail_lim; // pointer into assignment stack
unsigned m_clauses_to_reinit_lim; // ignore for now
bool m_inconsistent; // really needed?
};
svector<scope> m_scopes;
bool at_base_lvl() const { return m_scope_lvl == 0; }
lbool check() {
init_search();
propagate();
while (true) {
lbool r = bounded_search();
if (r != lbool::l_undef)
return r;
restart();
simplify_problem();
if (check_inconsistent()) return lbool::l_false;
gc();
}
}
cut_solver() {
}
void init_search() {
// TBD
// initialize data-structures
}
void simplify_problem() {
// no-op
}
void gc() {
// no-op
}
void restart() {
// no-op for now
}
bool check_inconsistent() {
// TBD
return false;
}
lbool bounded_search() {
while (true) {
checkpoint();
bool done = false;
while (!done) {
lbool is_sat = propagate_and_backjump_step(done);
if (is_sat != lbool::l_true) return is_sat;
}
gc();
if (!decide()) {
lbool is_sat = final_check();
if (is_sat != lbool::l_undef) {
return is_sat;
}
}
}
}
void checkpoint() {
// check for cancelation
}
void cleanup() {
}
lbool propagate_and_backjump_step(bool& done) {
done = true;
propagate();
if (!inconsistent())
return lbool::l_true;
if (!resolve_conflict())
return lbool::l_false;
if (at_base_lvl()) {
cleanup(); // cleaner may propagate frozen clauses
if (inconsistent()) {
TRACE("sat", tout << "conflict at level 0\n";);
return lbool::l_false;
}
gc();
}
done = false;
return lbool::l_true;
}
lbool final_check() {
// there are no more case splits, and all clauses are satisfied.
// prepare the model for external consumption.
return lbool::l_true;
}
bool resolve_conflict() {
while (true) {
bool r = resolve_conflict_core();
// after pop, clauses are reinitialized,
// this may trigger another conflict.
if (!r)
return false;
if (!inconsistent())
return true;
}
}
bool resolve_conflict_core() {
// this is where the main action is.
return true;
}
void propagate() {
// this is where the main action is.
}
bool decide() {
// this is where the main action is.
// pick the next variable and bound or value on the variable.
// return false if all variables have been assigned.
return false;
}
bool inconsistent() const { return m_inconsistent; }
};
}

334
src/util/lp/disjoint_intervals.h Normal file
View file

@ -0,0 +1,334 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include <map>
namespace lp {
// represents the set of disjoint intervals of integer number
struct disjoint_intervals {
std::map<int, short> m_endpoints; // 0 means start, 1 means end, 2 means both - for a point interval
bool m_empty;
// constructors create an interval containing all integer numbers or an empty interval
disjoint_intervals() : m_empty(false) {}
disjoint_intervals(bool is_empty) : m_empty(is_empty) {}
bool is_start(short x) const { return x == 0 || x == 2; }
bool is_start(const std::map<int, short>::iterator & it) const {
return is_start(it->second);
}
bool is_start(const std::map<int, short>::reverse_iterator & it) const {
return is_start(it->second);
}
bool is_end(short x) const { return x == 1 || x == 2; }
bool is_end(const std::map<int, short>::iterator & it) const {
return is_end(it->second);
}
bool is_end(const std::map<int, short>::reverse_iterator & it) const {
return is_end(it->second);
}
int pos(const std::map<int, short>::iterator & it) const {
return it->first;
}
int pos(const std::map<int, short>::reverse_iterator & it) const {
return it->first;
}
int bound_kind(const std::map<int, short>::iterator & it) const {
return it->second;
}
int bound_kind(const std::map<int, short>::reverse_iterator & it) const {
return it->second;
}
bool is_proper_start(short x) const { return x == 0; }
bool is_proper_end(short x) const { return x == 1; }
bool is_proper_end(const std::map<int, short>::iterator & it) const {
return is_proper_end(it->second);
}
bool is_proper_end(const std::map<int, short>::reverse_iterator & it) const {
return is_proper_end(it->second);
}
bool is_one_point_interval(short x) const { return x == 2; }
bool is_one_point_interval(const std::map<int, short>::iterator & it) const {
return is_one_point_interval(it->second);
}
bool is_one_point_interval(const std::map<int, short>::reverse_iterator & it) const {
return is_one_point_interval(it->second);
}
void erase(int x) {
m_endpoints.erase(x);
}
void set_one_point_segment(int x) {
m_endpoints[x] = 2;
}
void set_start(int x) {
m_endpoints[x] = 0;
}
void set_end(int x) {
m_endpoints[x] = 1;
}
void remove_all_endpoints_below(int x) {
while (m_endpoints.begin() != m_endpoints.end() && m_endpoints.begin()->first < x)
m_endpoints.erase(m_endpoints.begin());
}
// we intersect the existing set with the half open to the right interval
void intersect_with_lower_bound(int x) {
if (m_empty)
return;
if (m_endpoints.empty()) {
set_start(x);
return;
}
bool pos_inf = has_pos_inf();
auto it = m_endpoints.begin();
while (it != m_endpoints.end() && pos(it) < x) {
m_endpoints.erase(it);
it = m_endpoints.begin();
}
if (m_endpoints.empty()) {
if (!pos_inf) {
m_empty = true;
return;
}
set_start(x);
return;
}
lp_assert(pos(it) >= x);
if (pos(it) == x) {
if (is_proper_end(it))
set_one_point_segment(x);
}
else { // x(it) > x
if (is_proper_end(it)) {
set_start(x);
}
}
lp_assert(is_correct());
}
// we intersect the existing set with the half open interval
void intersect_with_upper_bound(int x) {
if (m_empty)
return;
if (m_endpoints.empty()) {
set_end(x);
return;
}
bool neg_inf = has_neg_inf();
auto it = m_endpoints.rbegin();
while (!m_endpoints.empty() && pos(it) > x) {
m_endpoints.erase(std::prev(m_endpoints.end()));
it = m_endpoints.rbegin();
}
if (m_endpoints.empty()) {
if (!neg_inf) {
m_empty = true;
return;
}
set_end(x);
}
lp_assert(pos(it) <= x);
if (pos(it) == x) {
if (is_one_point_interval(it)) {}
else if (is_proper_end(it)) {}
else {// is_proper_start(it->second)
set_one_point_segment(x);
}
}
else { // pos(it) < x}
if (is_start(it))
set_end(x);
}
lp_assert(is_correct());
}
bool has_pos_inf() const {
if (m_empty)
return false;
if (m_endpoints.empty())
return true;
lp_assert(m_endpoints.rbegin() != m_endpoints.rend());
return m_endpoints.rbegin()->second == 0;
}
bool has_neg_inf() const {
if (m_empty)
return false;
if (m_endpoints.empty())
return true;
auto it = m_endpoints.begin();
return is_proper_end(it->second);//m_endpoints.begin());
}
// we are intersecting
void intersect_with_interval(int x, int y) {
if (m_empty)
return;
lp_assert(x <= y);
intersect_with_lower_bound(x);
intersect_with_upper_bound(y);
}
// add an intervar [x, inf]
void unite_with_interval_x_pos_inf(int x) {
if (m_empty) {
set_start(x);
m_empty = false;
return;
}
while (!m_endpoints.empty() && pos(m_endpoints.rbegin()) > x) {
m_endpoints.erase(std::prev(m_endpoints.end()));
}
if (m_endpoints.empty()) {
set_start(x);
return;
}
auto it = m_endpoints.rbegin();
lp_assert(pos(it) <= x);
if (pos(it) == x) {
if (is_end(it)) {
m_endpoints.erase(x);
} else {
set_start(x);
}
} else if (pos(it) == x - 1 && is_end(it)) {
m_endpoints.erase(x - 1); // closing the gap
} else {
if (!has_pos_inf())
set_start(x);
}
}
// add an interval [-inf, x]
void unite_with_interval_neg_inf_x(int x) {
if (m_empty) {
set_end(x);
m_empty = false;
return;
}
auto it = m_endpoints.upper_bound(x);
if (it == m_endpoints.end()) {
bool pos_inf = has_pos_inf();
m_endpoints.clear();
// it could be the case where x is inside of the last infinite interval with pos inf
if (!pos_inf)
set_end(x);
return;
}
lp_assert(pos(it) > x);
if (is_one_point_interval(pos(it))) {
set_end(it->second);
} else {
if (is_start(it->second)) {
set_end(x);
}
}
while (!m_endpoints.empty() && m_endpoints.begin()->first < x) {
m_endpoints.erase(m_endpoints.begin());
}
lp_assert(is_correct());
}
void unite_with_interval(int x, int y) {
lp_assert(false); // not implemented
}
bool is_correct() const {
if (m_empty) {
if (m_endpoints.size() > 0) {
std::cout << "is empty is true but m_endpoints.size() = " << m_endpoints.size() << std::endl;
return false;
}
return true;
}
bool expect_end;
bool prev = false;
int prev_x;
for (auto t : m_endpoints) {
if (prev && (expect_end != t.second > 0)) {
std::cout << "x = " << t.first << "\n";
if (expect_end) {
std::cout << "expecting an interval end\n";
} else {
std::cout << "expecting an interval start\n";
}
return false;
}
if (t.second == 2) {
expect_end = false; // swallow a point interval
} else {
if (prev)
expect_end = !expect_end;
else
expect_end = is_start(t.second);
}
if (prev) {
if (t.first - prev_x <= 1) {
std::cout << "the sequence is not increasing or the gap is too small: " << prev_x << ", " << t.first << std::endl;
return false;
}
}
prev = true;
prev_x = t.first;
}
return true;
}
void print(std::ostream & out) const {
if (m_empty) {
out << "empty\n";
return;
}
if (m_endpoints.empty()){
out << "[-oo,oo]\n";
return;
}
bool first = true;
for (auto t : m_endpoints) {
if (first) {
if (t.second == 1) {
out << "[-oo," << t.first << "]";
}
else if (t.second == 0)
out << "[" << t.first << ",";
else if (t.second == 2)
out << "[" << t.first << "]";
first = false;
} else {
if (t.second==0)
out << "[" << t.first << ",";
else if (t.second == 1)
out << t.first << "]";
else if (t.second == 2)
out << "[" << t.first << "]";
}
}
if (has_pos_inf())
out << "oo]";
out << "\n";
}
};
}

591
src/util/lp/init_lar_solver.h Normal file
View file

@ -0,0 +1,591 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
// here we are inside lp::lar_solver class
bool strategy_is_undecided() const {
return m_settings.simplex_strategy() == simplex_strategy_enum::undecided;
}
var_index add_var(unsigned ext_j) {
var_index i;
SASSERT (ext_j < m_terms_start_index);
if (ext_j >= m_terms_start_index)
throw 0; // todo : what is the right way to exit?
if (try_get_val(m_ext_vars_to_columns, ext_j, i)) {
return i;
}
SASSERT(m_vars_to_ul_pairs.size() == A_r().column_count());
i = A_r().column_count();
m_vars_to_ul_pairs.push_back (ul_pair(static_cast<unsigned>(-1)));
add_non_basic_var_to_core_fields(ext_j);
SASSERT(sizes_are_correct());
return i;
}
void register_new_ext_var_index(unsigned ext_v) {
SASSERT(!contains(m_ext_vars_to_columns, ext_v));
unsigned j = static_cast<unsigned>(m_ext_vars_to_columns.size());
m_ext_vars_to_columns[ext_v] = j;
SASSERT(m_columns_to_ext_vars_or_term_indices.size() == j);
m_columns_to_ext_vars_or_term_indices.push_back(ext_v);
}
void add_non_basic_var_to_core_fields(unsigned ext_j) {
register_new_ext_var_index(ext_j);
m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
m_columns_with_changed_bound.increase_size_by_one();
add_new_var_to_core_fields_for_mpq(false);
if (use_lu())
add_new_var_to_core_fields_for_doubles(false);
}
void add_new_var_to_core_fields_for_doubles(bool register_in_basis) {
unsigned j = A_d().column_count();
A_d().add_column();
SASSERT(m_mpq_lar_core_solver.m_d_x.size() == j);
// SASSERT(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j); // restore later
m_mpq_lar_core_solver.m_d_x.resize(j + 1 );
m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
SASSERT(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
if (register_in_basis) {
A_d().add_row();
m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
m_mpq_lar_core_solver.m_d_basis.push_back(j);
}else {
m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
}
}
void add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
unsigned j = A_r().column_count();
A_r().add_column();
SASSERT(m_mpq_lar_core_solver.m_r_x.size() == j);
// SASSERT(m_mpq_lar_core_solver.m_r_low_bounds.size() == j && m_mpq_lar_core_solver.m_r_upper_bounds.size() == j); // restore later
m_mpq_lar_core_solver.m_r_x.resize(j + 1);
m_mpq_lar_core_solver.m_r_low_bounds.increase_size_by_one();
m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one();
m_mpq_lar_core_solver.m_r_solver.m_inf_set.increase_size_by_one();
m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1);
m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1);
SASSERT(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method
if (register_in_basis) {
A_r().add_row();
m_mpq_lar_core_solver.m_r_heading.push_back(m_mpq_lar_core_solver.m_r_basis.size());
m_mpq_lar_core_solver.m_r_basis.push_back(j);
if (m_settings.bound_propagation())
m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
} else {
m_mpq_lar_core_solver.m_r_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_r_nbasis.size()) - 1);
m_mpq_lar_core_solver.m_r_nbasis.push_back(j);
}
}
var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
const mpq &m_v) {
m_terms.push_back(new lar_term(coeffs, m_v));
m_orig_terms.push_back(new lar_term(coeffs, m_v));
return m_terms_start_index + m_terms.size() - 1;
}
// terms
var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
const mpq &m_v) {
if (strategy_is_undecided())
return add_term_undecided(coeffs, m_v);
m_terms.push_back(new lar_term(coeffs, m_v));
m_orig_terms.push_back(new lar_term(coeffs, m_v));
unsigned adjusted_term_index = m_terms.size() - 1;
var_index ret = m_terms_start_index + adjusted_term_index;
if (use_tableau() && !coeffs.empty()) {
add_row_for_term(m_orig_terms.back(), ret);
if (m_settings.bound_propagation())
m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
}
SASSERT(m_ext_vars_to_columns.size() == A_r().column_count());
return ret;
}
void add_row_for_term(const lar_term * term, unsigned term_ext_index) {
SASSERT(sizes_are_correct());
add_row_from_term_no_constraint(term, term_ext_index);
SASSERT(sizes_are_correct());
}
void add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) {
register_new_ext_var_index(term_ext_index);
// j will be a new variable
unsigned j = A_r().column_count();
ul_pair ul(j);
m_vars_to_ul_pairs.push_back(ul);
add_basic_var_to_core_fields();
if (use_tableau()) {
auto it = iterator_on_term_with_basis_var(*term, j);
A_r().fill_last_row_with_pivoting(it,
m_mpq_lar_core_solver.m_r_solver.m_basis_heading);
m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type<mpq>());
} else {
fill_last_row_of_A_r(A_r(), term);
}
m_mpq_lar_core_solver.m_r_x[j] = get_basic_var_value_from_row_directly(A_r().row_count() - 1);
if (use_lu())
fill_last_row_of_A_d(A_d(), term);
}
void add_basic_var_to_core_fields() {
bool use_lu = m_mpq_lar_core_solver.need_to_presolve_with_double_solver();
SASSERT(!use_lu || A_r().column_count() == A_d().column_count());
m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
m_columns_with_changed_bound.increase_size_by_one();
m_rows_with_changed_bounds.increase_size_by_one();
add_new_var_to_core_fields_for_mpq(true);
if (use_lu)
add_new_var_to_core_fields_for_doubles(true);
}
constraint_index add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side) {
constraint_index ci = m_constraints.size();
if (!is_term(j)) { // j is a var
auto vc = new lar_var_constraint(j, kind, right_side);
m_constraints.push_back(vc);
update_column_type_and_bound(j, kind, right_side, ci);
} else {
add_var_bound_on_constraint_for_term(j, kind, right_side, ci);
}
SASSERT(sizes_are_correct());
return ci;
}
void update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) {
switch(m_mpq_lar_core_solver.m_column_types[j]) {
case column_type::free_column:
update_free_column_type_and_bound(j, kind, right_side, constr_index);
break;
case column_type::boxed:
update_boxed_column_type_and_bound(j, kind, right_side, constr_index);
break;
case column_type::low_bound:
update_low_bound_column_type_and_bound(j, kind, right_side, constr_index);
break;
case column_type::upper_bound:
update_upper_bound_column_type_and_bound(j, kind, right_side, constr_index);
break;
case column_type::fixed:
update_fixed_column_type_and_bound(j, kind, right_side, constr_index);
break;
default:
SASSERT(false); // cannot be here
}
}
void add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
SASSERT(is_term(j));
unsigned adjusted_term_index = adjust_term_index(j);
unsigned term_j;
if (try_get_val(m_ext_vars_to_columns, j, term_j)) {
mpq rs = right_side - m_orig_terms[adjusted_term_index]->m_v;
m_constraints.push_back(new lar_term_constraint(m_orig_terms[adjusted_term_index], kind, right_side));
update_column_type_and_bound(term_j, kind, rs, ci);
}
else {
add_constraint_from_term_and_create_new_column_row(j, m_orig_terms[adjusted_term_index], kind, right_side);
}
}
void add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
lconstraint_kind kind, const mpq & right_side) {
add_row_from_term_no_constraint(term, term_j);
unsigned j = A_r().column_count() - 1;
update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size());
m_constraints.push_back(new lar_term_constraint(term, kind, right_side));
SASSERT(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
}
void decide_on_strategy_and_adjust_initial_state() {
SASSERT(strategy_is_undecided());
if (m_vars_to_ul_pairs.size() > m_settings.column_number_threshold_for_using_lu_in_lar_solver) {
m_settings.simplex_strategy() = simplex_strategy_enum::lu;
} else {
m_settings.simplex_strategy() = simplex_strategy_enum::tableau_rows; // todo: when to switch to tableau_costs?
}
adjust_initial_state();
}
void adjust_initial_state() {
switch (m_settings.simplex_strategy()) {
case simplex_strategy_enum::lu:
adjust_initial_state_for_lu();
break;
case simplex_strategy_enum::tableau_rows:
adjust_initial_state_for_tableau_rows();
break;
case simplex_strategy_enum::tableau_costs:
SASSERT(false); // not implemented
case simplex_strategy_enum::undecided:
adjust_initial_state_for_tableau_rows();
break;
}
}
void adjust_initial_state_for_lu() {
copy_from_mpq_matrix(A_d());
unsigned n = A_d().column_count();
m_mpq_lar_core_solver.m_d_x.resize(n);
m_mpq_lar_core_solver.m_d_low_bounds.resize(n);
m_mpq_lar_core_solver.m_d_upper_bounds.resize(n);
m_mpq_lar_core_solver.m_d_heading = m_mpq_lar_core_solver.m_r_heading;
m_mpq_lar_core_solver.m_d_basis = m_mpq_lar_core_solver.m_r_basis;
/*
unsigned j = A_d().column_count();
A_d().add_column();
SASSERT(m_mpq_lar_core_solver.m_d_x.size() == j);
// SASSERT(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j); // restore later
m_mpq_lar_core_solver.m_d_x.resize(j + 1 );
m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
SASSERT(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
if (register_in_basis) {
A_d().add_row();
m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
m_mpq_lar_core_solver.m_d_basis.push_back(j);
}else {
m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
}*/
}
void adjust_initial_state_for_tableau_rows() {
for (unsigned j = 0; j < m_terms.size(); j++) {
if (contains(m_ext_vars_to_columns, j + m_terms_start_index))
continue;
add_row_from_term_no_constraint(m_terms[j], j + m_terms_start_index);
}
}
// this fills the last row of A_d and sets the basis column: -1 in the last column of the row
void fill_last_row_of_A_d(static_matrix<double, double> & A, const lar_term* ls) {
SASSERT(A.row_count() > 0);
SASSERT(A.column_count() > 0);
unsigned last_row = A.row_count() - 1;
SASSERT(A.m_rows[last_row].empty());
for (auto & t : ls->m_coeffs) {
SASSERT(!is_zero(t.second));
var_index j = t.first;
A.set(last_row, j, - t.second.get_double());
}
unsigned basis_j = A.column_count() - 1;
A.set(last_row, basis_j, - 1 );
}
void update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) {
mpq y_of_bound(0);
switch (kind) {
case LT:
y_of_bound = -1;
case LE:
m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound;
SASSERT(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
SASSERT(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
{
auto up = numeric_pair<mpq>(right_side, y_of_bound);
m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
}
set_upper_bound_witness(j, constr_ind);
break;
case GT:
y_of_bound = 1;
case GE:
m_mpq_lar_core_solver.m_column_types[j] = column_type::low_bound;
SASSERT(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
{
auto low = numeric_pair<mpq>(right_side, y_of_bound);
m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
}
set_low_bound_witness(j, constr_ind);
break;
case EQ:
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
set_upper_bound_witness(j, constr_ind);
set_low_bound_witness(j, constr_ind);
break;
default:
SASSERT(false);
}
m_columns_with_changed_bound.insert(j);
}
void update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
SASSERT(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
mpq y_of_bound(0);
switch (kind) {
case LT:
y_of_bound = -1;
case LE:
{
auto up = numeric_pair<mpq>(right_side, y_of_bound);
if (up < m_mpq_lar_core_solver.m_r_upper_bounds()[j]) {
m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
set_upper_bound_witness(j, ci);
m_columns_with_changed_bound.insert(j);
}
}
break;
case GT:
y_of_bound = 1;
case GE:
m_mpq_lar_core_solver.m_column_types[j] = column_type::boxed;
{
auto low = numeric_pair<mpq>(right_side, y_of_bound);
m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
set_low_bound_witness(j, ci);
m_columns_with_changed_bound.insert(j);
if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
} else {
m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
}
}
break;
case EQ:
{
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
set_low_bound_witness(j, ci);
m_infeasible_column_index = j;
} else {
m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
m_columns_with_changed_bound.insert(j);
set_low_bound_witness(j, ci);
set_upper_bound_witness(j, ci);
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
}
break;
}
break;
default:
SASSERT(false);
}
}
void update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
SASSERT(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
mpq y_of_bound(0);
switch (kind) {
case LT:
y_of_bound = -1;
case LE:
{
auto up = numeric_pair<mpq>(right_side, y_of_bound);
if (up < m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
set_upper_bound_witness(j, ci);
m_columns_with_changed_bound.insert(j);
}
if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
SASSERT(false);
m_infeasible_column_index = j;
} else {
if (m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j])
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
}
}
break;
case GT:
y_of_bound = 1;
case GE:
{
auto low = numeric_pair<mpq>(right_side, y_of_bound);
if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
m_columns_with_changed_bound.insert(j);
set_low_bound_witness(j, ci);
}
if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
} else if ( low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
}
}
break;
case EQ:
{
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
} else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_low_bound_witness(j, ci);
} else {
m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
set_low_bound_witness(j, ci);
set_upper_bound_witness(j, ci);
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
m_columns_with_changed_bound.insert(j);
}
break;
}
default:
SASSERT(false);
}
}
void update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
SASSERT(m_mpq_lar_core_solver.m_column_types()[j] == column_type::low_bound);
mpq y_of_bound(0);
switch (kind) {
case LT:
y_of_bound = -1;
case LE:
{
auto up = numeric_pair<mpq>(right_side, y_of_bound);
m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
set_upper_bound_witness(j, ci);
m_columns_with_changed_bound.insert(j);
if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
} else {
m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
}
}
break;
case GT:
y_of_bound = 1;
case GE:
{
auto low = numeric_pair<mpq>(right_side, y_of_bound);
if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
m_columns_with_changed_bound.insert(j);
set_low_bound_witness(j, ci);
}
}
break;
case EQ:
{
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
} else {
m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
set_low_bound_witness(j, ci);
set_upper_bound_witness(j, ci);
m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
}
m_columns_with_changed_bound.insert(j);
break;
}
default:
SASSERT(false);
}
}
void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
SASSERT(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
SASSERT(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero()));
auto v = numeric_pair<mpq>(right_side, mpq(0));
mpq y_of_bound(0);
switch (kind) {
case LT:
if (v <= m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
}
break;
case LE:
{
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
}
}
break;
case GT:
{
if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index =j;
set_low_bound_witness(j, ci);
}
}
break;
case GE:
{
if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_low_bound_witness(j, ci);
}
}
break;
case EQ:
{
if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_upper_bound_witness(j, ci);
} else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
m_status = INFEASIBLE;
m_infeasible_column_index = j;
set_low_bound_witness(j, ci);
}
break;
}
default:
SASSERT(false);
}
}

408
src/util/lp/lp_primal_core_solver_tableau.h Normal file
View file

@ -0,0 +1,408 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
// this is a part of lp_primal_core_solver that deals with the tableau
#include "util/lp/lp_primal_core_solver.h"
namespace lp {
template <typename T, typename X> void lp_primal_core_solver<T, X>::one_iteration_tableau() {
int entering = choose_entering_column_tableau();
if (entering == -1) {
decide_on_status_when_cannot_find_entering();
}
else {
advance_on_entering_tableau(entering);
}
SASSERT(this->inf_set_is_correct());
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_entering_tableau(int entering) {
X t;
int leaving = find_leaving_and_t_tableau(entering, t);
if (leaving == -1) {
this->set_status(UNBOUNDED);
return;
}
advance_on_entering_and_leaving_tableau(entering, leaving, t);
}
/*
template <typename T, typename X> int lp_primal_core_solver<T, X>::choose_entering_column_tableau_rows() {
int i = find_inf_row();
if (i == -1)
return -1;
return find_shortest_beneficial_column_in_row(i);
}
*/
template <typename T, typename X> int lp_primal_core_solver<T, X>::choose_entering_column_tableau() {
//this moment m_y = cB * B(-1)
unsigned number_of_benefitial_columns_to_go_over = get_number_of_non_basic_column_to_try_for_enter();
SASSERT(numeric_traits<T>::precise());
if (number_of_benefitial_columns_to_go_over == 0)
return -1;
if (this->m_basis_sort_counter == 0) {
sort_non_basis();
this->m_basis_sort_counter = 20;
}
else {
this->m_basis_sort_counter--;
}
unsigned j_nz = this->m_m() + 1; // this number is greater than the max column size
std::list<unsigned>::iterator entering_iter = m_non_basis_list.end();
for (auto non_basis_iter = m_non_basis_list.begin(); number_of_benefitial_columns_to_go_over && non_basis_iter != m_non_basis_list.end(); ++non_basis_iter) {
unsigned j = *non_basis_iter;
if (!column_is_benefitial_for_entering_basis(j))
continue;
// if we are here then j is a candidate to enter the basis
unsigned t = this->m_A.number_of_non_zeroes_in_column(j);
if (t < j_nz) {
j_nz = t;
entering_iter = non_basis_iter;
if (number_of_benefitial_columns_to_go_over)
number_of_benefitial_columns_to_go_over--;
}
else if (t == j_nz && this->m_settings.random_next() % 2 == 0) {
entering_iter = non_basis_iter;
}
}// while (number_of_benefitial_columns_to_go_over && initial_offset_in_non_basis != offset_in_nb);
if (entering_iter == m_non_basis_list.end())
return -1;
unsigned entering = *entering_iter;
m_sign_of_entering_delta = this->m_d[entering] > 0 ? 1 : -1;
if (this->m_using_infeas_costs && this->m_settings.use_breakpoints_in_feasibility_search)
m_sign_of_entering_delta = -m_sign_of_entering_delta;
m_non_basis_list.erase(entering_iter);
m_non_basis_list.push_back(entering);
return entering;
}
template <typename T, typename X>
unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
init_run_tableau();
if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only) {
this->set_status(FEASIBLE);
return 0;
}
if ((!numeric_traits<T>::precise()) && this->A_mult_x_is_off()) {
this->set_status(FLOATING_POINT_ERROR);
return 0;
}
do {
if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over((this->m_using_infeas_costs? "inf t" : "feas t"), * this->m_settings.get_message_ostream())) {
return this->total_iterations();
}
if (this->m_settings.use_tableau_rows())
one_iteration_tableau_rows();
else
one_iteration_tableau();
switch (this->get_status()) {
case OPTIMAL: // double check that we are at optimum
case INFEASIBLE:
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
break;
if (!numeric_traits<T>::precise()) {
if(this->m_look_for_feasible_solution_only)
break;
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
if (choose_entering_column(1) == -1) {
decide_on_status_when_cannot_find_entering();
break;
}
this->set_status(UNKNOWN);
} else { // precise case
if ((!this->infeasibility_costs_are_correct())) {
init_reduced_costs_tableau(); // forcing recalc
if (choose_entering_column_tableau() == -1) {
decide_on_status_when_cannot_find_entering();
break;
}
this->set_status(UNKNOWN);
}
}
break;
case TENTATIVE_UNBOUNDED:
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
break;
case UNBOUNDED:
if (this->current_x_is_infeasible()) {
init_reduced_costs();
this->set_status(UNKNOWN);
}
break;
case UNSTABLE:
SASSERT(! (numeric_traits<T>::precise()));
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
break;
default:
break; // do nothing
}
} while (this->get_status() != FLOATING_POINT_ERROR
&&
this->get_status() != UNBOUNDED
&&
this->get_status() != OPTIMAL
&&
this->get_status() != INFEASIBLE
&&
this->iters_with_no_cost_growing() <= this->m_settings.max_number_of_iterations_with_no_improvements
&&
this->total_iterations() <= this->m_settings.max_total_number_of_iterations
&&
!(this->current_x_is_feasible() && this->m_look_for_feasible_solution_only));
SASSERT(this->get_status() == FLOATING_POINT_ERROR
||
this->current_x_is_feasible() == false
||
this->calc_current_x_is_feasible_include_non_basis());
return this->total_iterations();
}
template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_entering_and_leaving_tableau(int entering, int leaving, X & t) {
SASSERT(this->A_mult_x_is_off() == false);
SASSERT(leaving >= 0 && entering >= 0);
SASSERT((this->m_settings.simplex_strategy() ==
simplex_strategy_enum::tableau_rows) ||
m_non_basis_list.back() == static_cast<unsigned>(entering));
SASSERT(this->m_using_infeas_costs || !is_neg(t));
SASSERT(entering != leaving || !is_zero(t)); // otherwise nothing changes
if (entering == leaving) {
advance_on_entering_equal_leaving_tableau(entering, t);
return;
}
if (!is_zero(t)) {
if (this->current_x_is_feasible() || !this->m_settings.use_breakpoints_in_feasibility_search ) {
if (m_sign_of_entering_delta == -1)
t = -t;
}
this->update_basis_and_x_tableau(entering, leaving, t);
SASSERT(this->A_mult_x_is_off() == false);
this->iters_with_no_cost_growing() = 0;
} else {
this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
this->change_basis(entering, leaving);
}
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
return;
if (this->m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows) {
if (need_to_switch_costs()) {
this->init_reduced_costs_tableau();
}
SASSERT(!need_to_switch_costs());
std::list<unsigned>::iterator it = m_non_basis_list.end();
it--;
* it = static_cast<unsigned>(leaving);
}
}
template <typename T, typename X>
void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving_tableau(int entering, X & t) {
SASSERT(!this->A_mult_x_is_off() );
this->update_x_tableau(entering, t * m_sign_of_entering_delta);
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
return;
if (need_to_switch_costs()) {
init_reduced_costs_tableau();
}
this->iters_with_no_cost_growing() = 0;
}
template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_and_t_tableau(unsigned entering, X & t) {
unsigned k = 0;
bool unlimited = true;
unsigned row_min_nz = this->m_n() + 1;
m_leaving_candidates.clear();
auto & col = this->m_A.m_columns[entering];
unsigned col_size = col.size();
for (;k < col_size && unlimited; k++) {
const column_cell & c = col[k];
unsigned i = c.m_i;
const T & ed = this->m_A.get_val(c);
SASSERT(!numeric_traits<T>::is_zero(ed));
unsigned j = this->m_basis[i];
limit_theta_on_basis_column(j, - ed * m_sign_of_entering_delta, t, unlimited);
if (!unlimited) {
m_leaving_candidates.push_back(j);
row_min_nz = this->m_A.m_rows[i].size();
}
}
if (unlimited) {
if (try_jump_to_another_bound_on_entering_unlimited(entering, t))
return entering;
return -1;
}
X ratio;
for (;k < col_size; k++) {
const column_cell & c = col[k];
unsigned i = c.m_i;
const T & ed = this->m_A.get_val(c);
SASSERT(!numeric_traits<T>::is_zero(ed));
unsigned j = this->m_basis[i];
unlimited = true;
limit_theta_on_basis_column(j, -ed * m_sign_of_entering_delta, ratio, unlimited);
if (unlimited) continue;
unsigned i_nz = this->m_A.m_rows[i].size();
if (ratio < t) {
t = ratio;
m_leaving_candidates.clear();
m_leaving_candidates.push_back(j);
row_min_nz = i_nz;
} else if (ratio == t && i_nz < row_min_nz) {
m_leaving_candidates.clear();
m_leaving_candidates.push_back(j);
row_min_nz = this->m_A.m_rows[i].size();
} else if (ratio == t && i_nz == row_min_nz) {
m_leaving_candidates.push_back(j);
}
}
ratio = t;
unlimited = false;
if (try_jump_to_another_bound_on_entering(entering, t, ratio, unlimited)) {
t = ratio;
return entering;
}
if (m_leaving_candidates.size() == 1)
return m_leaving_candidates[0];
k = this->m_settings.random_next() % m_leaving_candidates.size();
return m_leaving_candidates[k];
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tableau() {
// print_matrix(&(this->m_A), std::cout);
SASSERT(this->A_mult_x_is_off() == false);
SASSERT(basis_columns_are_set_correctly());
this->m_basis_sort_counter = 0; // to initiate the sort of the basis
this->set_total_iterations(0);
this->iters_with_no_cost_growing() = 0;
SASSERT(this->inf_set_is_correct());
if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only)
return;
if (this->m_settings.backup_costs)
backup_and_normalize_costs();
m_epsilon_of_reduced_cost = numeric_traits<X>::precise() ? zero_of_type<T>() : T(1) / T(10000000);
if (this->m_settings.use_breakpoints_in_feasibility_search)
m_breakpoint_indices_queue.resize(this->m_n());
if (!numeric_traits<X>::precise()) {
this->m_column_norm_update_counter = 0;
init_column_norms();
}
if (this->m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows)
init_tableau_rows();
SASSERT(this->reduced_costs_are_correct_tableau());
SASSERT(!this->need_to_pivot_to_basis_tableau());
}
template <typename T, typename X> bool lp_primal_core_solver<T, X>::
update_basis_and_x_tableau(int entering, int leaving, X const & tt) {
SASSERT(this->use_tableau());
update_x_tableau(entering, tt);
this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
this->change_basis(entering, leaving);
return true;
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::
update_x_tableau(unsigned entering, const X& delta) {
if (!this->m_using_infeas_costs) {
this->m_x[entering] += delta;
for (const auto & c : this->m_A.m_columns[entering]) {
unsigned i = c.m_i;
this->m_x[this->m_basis[i]] -= delta * this->m_A.get_val(c);
this->update_column_in_inf_set(this->m_basis[i]);
}
} else { // m_using_infeas_costs == true
this->m_x[entering] += delta;
SASSERT(this->column_is_feasible(entering));
SASSERT(this->m_costs[entering] == zero_of_type<T>());
// m_d[entering] can change because of the cost change for basic columns.
for (const auto & c : this->m_A.m_columns[entering]) {
unsigned i = c.m_i;
unsigned j = this->m_basis[i];
this->m_x[j] -= delta * this->m_A.get_val(c);
update_inf_cost_for_column_tableau(j);
if (is_zero(this->m_costs[j]))
this->m_inf_set.erase(j);
else
this->m_inf_set.insert(j);
}
}
SASSERT(this->A_mult_x_is_off() == false);
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::
update_inf_cost_for_column_tableau(unsigned j) {
SASSERT(this->m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows);
SASSERT(this->m_using_infeas_costs);
T new_cost = get_infeasibility_cost_for_column(j);
T delta = this->m_costs[j] - new_cost;
if (is_zero(delta))
return;
this->m_costs[j] = new_cost;
update_reduced_cost_for_basic_column_cost_change(delta, j);
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::init_reduced_costs_tableau() {
if (this->current_x_is_infeasible() && !this->m_using_infeas_costs) {
init_infeasibility_costs();
} else if (this->current_x_is_feasible() && this->m_using_infeas_costs) {
if (this->m_look_for_feasible_solution_only)
return;
this->m_costs = m_costs_backup;
this->m_using_infeas_costs = false;
}
unsigned size = this->m_basis_heading.size();
for (unsigned j = 0; j < size; j++) {
if (this->m_basis_heading[j] >= 0)
this->m_d[j] = zero_of_type<T>();
else {
T& d = this->m_d[j] = this->m_costs[j];
for (auto & cc : this->m_A.m_columns[j]) {
d -= this->m_costs[this->m_basis[cc.m_i]] * this->m_A.get_val(cc);
}
}
}
}
}

202
src/util/sorting_network.h
View file

@ -24,6 +24,28 @@ Notes:
#ifndef SORTING_NETWORK_H_
#define SORTING_NETWORK_H_
enum sorting_network_encoding {
grouped_at_most_1,
bimander_at_most_1,
ordered_at_most_1
};
inline std::ostream& operator<<(std::ostream& out, sorting_network_encoding enc) {
switch (enc) {
case grouped_at_most_1: return out << "grouped";
case bimander_at_most_1: return out << "bimander";
case ordered_at_most_1: return out << "ordered";
}
return out << "???";
}
struct sorting_network_config {
sorting_network_encoding m_encoding;
sorting_network_config() {
m_encoding = grouped_at_most_1;
}
};
template <typename Ext>
class sorting_network {
typedef typename Ext::vector vect;
@ -88,7 +110,7 @@ Notes:
}
public:
sorting_network(Ext& ext):
sorting_network(Ext& ext):
m_ext(ext),
m_current(&m_currentv),
m_next(&m_nextv)
@ -121,6 +143,7 @@ Notes:
class psort_nw {
typedef typename psort_expr::literal literal;
typedef typename psort_expr::literal_vector literal_vector;
sorting_network_config m_cfg;
class vc {
unsigned v; // number of vertices
@ -185,7 +208,7 @@ Notes:
}
};
psort_nw(psort_expr& c): ctx(c) {}
psort_nw(psort_expr& c, sorting_network_config const& cfg): ctx(c), m_cfg(cfg) {}
literal ge(bool full, unsigned k, unsigned n, literal const* xs) {
if (k > n) {
@ -220,7 +243,17 @@ Notes:
else if (k == 1) {
literal_vector ors;
// scoped_stats _ss(m_stats, k, n);
return mk_at_most_1(full, n, xs, ors, false);
switch (m_cfg.m_encoding) {
case grouped_at_most_1:
return mk_at_most_1(full, n, xs, ors, false);
case bimander_at_most_1:
return mk_at_most_1_bimander(full, n, xs, ors);
case ordered_at_most_1:
return mk_ordered_atmost_1(full, n, xs);
default:
UNREACHABLE();
return xs[0];
}
}
else {
SASSERT(2*k <= n);
@ -278,7 +311,7 @@ Notes:
if (n == 1) {
return ors[0];
}
literal result = fresh();
literal result = fresh("or");
add_implies_or(result, n, ors);
add_or_implies(result, n, ors);
return result;
@ -317,15 +350,26 @@ Notes:
if (ands.size() == 1) {
return ands[0];
}
literal result = fresh();
literal result = fresh("and");
add_implies_and(result, ands);
add_and_implies(result, ands);
return result;
}
literal mk_exactly_1(bool full, unsigned n, literal const* xs) {
TRACE("pb", tout << "exactly 1 with " << n << " arguments " << (full?"full":"not full") << "\n";);
literal_vector ors;
literal r1 = mk_at_most_1(full, n, xs, ors, true);
literal r1;
switch (m_cfg.m_encoding) {
case grouped_at_most_1:
r1 = mk_at_most_1(full, n, xs, ors, true);
break;
case bimander_at_most_1:
r1 = mk_at_most_1_bimander(full, n, xs, ors);
break;
case ordered_at_most_1:
return mk_ordered_exactly_1(full, n, xs);
}
if (full) {
r1 = mk_and(r1, mk_or(ors));
@ -340,12 +384,8 @@ Notes:
TRACE("pb_verbose", tout << (full?"full":"partial") << " ";
for (unsigned i = 0; i < n; ++i) tout << xs[i] << " ";
tout << "\n";);
if (n >= 4 && false) {
return mk_at_most_1_bimander(full, n, xs, ors);
}
literal_vector in(n, xs);
literal result = fresh();
literal result = fresh("at-most-1");
unsigned inc_size = 4;
literal_vector ands;
ands.push_back(result);
@ -387,7 +427,7 @@ Notes:
// xs[0] + ... + xs[n-1] <= 1 => and_x
if (full) {
literal and_i = fresh();
literal and_i = fresh("and");
for (unsigned i = 0; i < n; ++i) {
literal_vector lits;
lits.push_back(and_i);
@ -407,29 +447,6 @@ Notes:
}
#if 0
literal result = fresh();
// result => xs[0] + ... + xs[n-1] <= 1
for (unsigned i = 0; i < n; ++i) {
for (unsigned j = i + 1; j < n; ++j) {
add_clause(ctx.mk_not(result), ctx.mk_not(xs[i]), ctx.mk_not(xs[j]));
}
}
// xs[0] + ... + xs[n-1] <= 1 => result
for (unsigned i = 0; i < n; ++i) {
literal_vector lits;
lits.push_back(result);
for (unsigned j = 0; j < n; ++j) {
if (j != i) lits.push_back(xs[j]);
}
add_clause(lits);
}
return result;
#endif
#if 1
// r <=> and( or(!xi,!xj))
//
literal_vector ands;
@ -439,30 +456,100 @@ Notes:
}
}
return mk_and(ands);
#else
// r <=> or (and !x_{j != i})
literal_vector ors;
for (unsigned i = 0; i < n; ++i) {
literal_vector ands;
for (unsigned j = 0; j < n; ++j) {
if (j != i) {
ands.push_back(ctx.mk_not(xs[j]));
}
}
ors.push_back(mk_and(ands));
}
return mk_or(ors);
#endif
}
literal mk_ordered_exactly_1(bool full, unsigned n, literal const* xs) {
return mk_ordered_1(full, true, n, xs);
}
literal mk_ordered_atmost_1(bool full, unsigned n, literal const* xs) {
return mk_ordered_1(full, false, n, xs);
}
literal mk_ordered_1(bool full, bool is_eq, unsigned n, literal const* xs) {
if (n <= 1 && !is_eq) return ctx.mk_true();
if (n == 0) {
return ctx.mk_false();
}
if (n == 1) {
return xs[0];
}
// y0 -> y1
// x0 -> y0
// x1 -> y1
// r, y0 -> ~x1
// r, y1 -> ~x2
// r -> x3 | y1
// r -> ~x3 | ~y1
// x0,x1,x2, .., x_{n-1}, x_n
// y0,y1,y2, .., y_{n-1}
// y_i -> y_{i+1} i = 0, ..., n - 2
// x_i -> y_i i = 0, ..., n - 1
// r, y_i -> ~x_{i+1} i = 0, ..., n - 1
// exactly 1:
// r -> x_n | y_{n-1}
// full (exactly 1):
// two_i -> y_i & x_{i+1}
// zero -> ~x_n
// zero -> ~y_{n-1}
// r | zero | two_0 | ... | two_{n-1}
// full atmost 1:
// r | two | two_0 | ... | two_{n-1}
literal r = fresh("ordered");
literal_vector ys;
for (unsigned i = 0; i + 1 < n; ++i) {
ys.push_back(fresh("y"));
}
for (unsigned i = 0; i + 2 < n; ++i) {
add_clause(ctx.mk_not(ys[i]), ys[i + 1]);
}
for (unsigned i = 0; i + 1 < n; ++i) {
add_clause(ctx.mk_not(xs[i]), ys[i]);
add_clause(ctx.mk_not(r), ctx.mk_not(ys[i]), ctx.mk_not(xs[i + 1]));
}
add_clause(ctx.mk_not(r), ys[n-2], xs[n-1]);
add_clause(ctx.mk_not(r), ctx.mk_not(ys[n-2]), ctx.mk_not(xs[n-1]));
for (unsigned i = 1; i < n - 1; ++i) {
add_clause(ctx.mk_not(ys[i]), xs[i], ys[i-1]);
}
add_clause(ctx.mk_not(ys[0]), xs[0]);
if (full) {
literal_vector twos;
for (unsigned i = 0; i < n - 1; ++i) {
twos.push_back(fresh("two"));
}
add_clause(ctx.mk_not(twos[0]), ys[0]);
add_clause(ctx.mk_not(twos[0]), xs[1]);
for (unsigned i = 1; i < n - 1; ++i) {
add_clause(ctx.mk_not(twos[i]), ys[i], twos[i-1]);
add_clause(ctx.mk_not(twos[i]), xs[i + 1], twos[i-1]);
}
if (is_eq) {
literal zero = fresh("zero");
add_clause(ctx.mk_not(zero), ctx.mk_not(xs[n-1]));
add_clause(ctx.mk_not(zero), ctx.mk_not(ys[n-2]));
add_clause(r, zero, twos.back());
}
else {
add_clause(r, twos.back());
}
}
return r;
}
//
literal mk_at_most_1_bimander(bool full, unsigned n, literal const* xs, literal_vector& ors) {
if (full) {
return mk_at_most_1(full, n, xs, ors, true);
}
literal_vector in(n, xs);
literal result = fresh();
literal result = fresh("bimander");
unsigned inc_size = 2;
literal_vector ands;
for (unsigned i = 0; i < n; i += inc_size) {
@ -477,7 +564,7 @@ Notes:
}
literal_vector bits;
for (unsigned k = 0; k < nbits; ++k) {
bits.push_back(fresh());
bits.push_back(fresh("bit"));
}
for (unsigned i = 0; i < ors.size(); ++i) {
for (unsigned k = 0; k < nbits; ++k) {
@ -539,9 +626,9 @@ Notes:
return ctx.mk_min(a, b);
}
literal fresh() {
literal fresh(char const* n) {
m_stats.m_num_compiled_vars++;
return ctx.fresh();
return ctx.fresh(n);
}
void add_clause(literal l1, literal l2, literal l3) {
literal lits[3] = { l1, l2, l3 };
@ -558,7 +645,6 @@ Notes:
m_stats.m_num_compiled_clauses++;
m_stats.m_num_clause_vars += n;
literal_vector tmp(n, ls);
TRACE("pb_verbose", for (unsigned i = 0; i < n; ++i) tout << ls[i] << " "; tout << "\n";);
ctx.mk_clause(n, tmp.c_ptr());
}
@ -925,7 +1011,7 @@ Notes:
SASSERT(b <= c);
SASSERT(a + b >= c);
for (unsigned i = 0; i < c; ++i) {
out.push_back(fresh());
out.push_back(fresh("dsmerge"));
}
if (m_t != GE) {
for (unsigned i = 0; i < a; ++i) {
@ -983,7 +1069,7 @@ Notes:
SASSERT(m <= n);
literal_vector lits;
for (unsigned i = 0; i < m; ++i) {
out.push_back(fresh());
out.push_back(fresh("dsort"));
}
if (m_t != GE) {
for (unsigned k = 1; k <= m; ++k) {