diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 47f5d3028..c240d243a 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -307,7 +307,7 @@ class theory_lra::imp { m_solver->settings().bound_propagation() = BP_NONE != propagation_mode(); m_solver->set_track_pivoted_rows(lp.bprop_on_pivoted_rows()); m_solver->settings().m_int_branch_cut_gomory_threshold = ctx().get_fparams().m_arith_branch_cut_ratio; - m_solver->settings().m_int_branch_cut_solver = std::max(4u, ctx().get_fparams().m_arith_branch_cut_ratio); + m_solver->settings().m_int_branch_cut_solver = std::max(8u, ctx().get_fparams().m_arith_branch_cut_ratio); m_solver->settings().m_run_gcd_test = ctx().get_fparams().m_arith_gcd_test; m_solver->settings().set_random_seed(ctx().get_fparams().m_random_seed); //m_solver->settings().set_ostream(0); @@ -672,6 +672,10 @@ class theory_lra::imp { } theory_var internalize_def(app* term, scoped_internalize_state& st) { + if (ctx().e_internalized(term)) { + IF_VERBOSE(0, verbose_stream() << "repeated term\n";); + return mk_var(term, false); + } linearize_term(term, st); if (is_unit_var(st)) { return st.vars()[0]; @@ -1252,17 +1256,21 @@ public: } // create a bound atom representing term <= k - app_ref mk_bound(lp::lar_term const& term, rational const& k) { + app_ref mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound) { app_ref t = mk_term(term, k.is_int()); - app_ref atom(a.mk_le(t, a.mk_numeral(k, k.is_int())), m); + app_ref atom(m); + if (lower_bound) { + atom = a.mk_ge(t, a.mk_numeral(k, k.is_int())); + } + else { + atom = a.mk_le(t, a.mk_numeral(k, k.is_int())); + } expr_ref atom1(m); proof_ref atomp(m); ctx().get_rewriter()(atom, atom1, atomp); atom = to_app(atom1); TRACE("arith", tout << atom << "\n"; - m_solver->print_term(term, tout << "bound atom: "); tout << " <= " << k << "\n"; - display(tout); - ); + m_solver->print_term(term, tout << "bound atom: "); tout << " <= " << k << "\n";); ctx().internalize(atom, true); ctx().mark_as_relevant(atom.get()); return atom; @@ -1276,12 +1284,15 @@ public: lp::lar_term term; lp::mpq k; lp::explanation ex; // TBD, this should be streamlined accross different explanations - switch(m_lia->check(term, k, ex)) { + bool upper; + switch(m_lia->check(term, k, ex, upper)) { case lp::lia_move::ok: return l_true; case lp::lia_move::branch: { - (void)mk_bound(term, k); + app_ref b = mk_bound(term, k, upper); + // branch on term >= k + 1 // branch on term <= k + // TBD: ctx().force_phase(ctx().get_literal(b)); // at this point we have a new unassigned atom that the // SAT core assigns a value to return l_false; @@ -1289,7 +1300,7 @@ public: case lp::lia_move::cut: { ++m_stats.m_gomory_cuts; // m_explanation implies term <= k - app_ref b = mk_bound(term, k); + app_ref b = mk_bound(term, k, upper); m_eqs.reset(); m_core.reset(); m_params.reset(); @@ -1298,7 +1309,11 @@ public: set_evidence(ev.second); } } - assign(literal(ctx().get_bool_var(b), false)); + literal lit(ctx().get_bool_var(b), false); + TRACE("arith", + ctx().display_lemma_as_smt_problem(tout << "new cut:\n", m_core.size(), m_core.c_ptr(), m_eqs.size(), m_eqs.c_ptr(), lit); + display(tout);); + assign(lit); return l_false; } case lp::lia_move::conflict: @@ -2725,6 +2740,8 @@ public: if (m_solver) { m_solver->print_constraints(out); m_solver->print_terms(out); + // auto pp = lp ::core_solver_pretty_printer(m_solver->m_mpq_lar_core_solver.m_r_solver, out); + // pp.print(); } unsigned nv = th.get_num_vars(); for (unsigned v = 0; v < nv; ++v) { @@ -2793,6 +2810,8 @@ public: st.update("cut_solver-true", m_solver->settings().st().m_cut_solver_true); st.update("cut_solver-false", m_solver->settings().st().m_cut_solver_false); st.update("cut_solver-undef", m_solver->settings().st().m_cut_solver_undef); + st.update("gcd_calls", m_solver->settings().st().m_gcd_calls); + st.update("gcd_conflict", m_solver->settings().st().m_gcd_conflicts); } }; diff --git a/src/test/lp/gomory_test.h b/src/test/lp/gomory_test.h new file mode 100644 index 000000000..03ad5b187 --- /dev/null +++ b/src/test/lp/gomory_test.h @@ -0,0 +1,236 @@ +namespace lp { +struct gomory_test { + gomory_test( + std::function name_function_p, + std::function get_value_p, + std::function at_low_p, + std::function at_upper_p, + std::function lower_bound_p, + std::function upper_bound_p, + std::function column_lower_bound_constraint_p, + std::function column_upper_bound_constraint_p + ) : + m_name_function(name_function_p), + get_value(get_value_p), + at_low(at_low_p), + at_upper(at_upper_p), + lower_bound(lower_bound_p), + upper_bound(upper_bound_p), + column_lower_bound_constraint(column_lower_bound_constraint_p), + column_upper_bound_constraint(column_upper_bound_constraint_p) + {} + + std::function m_name_function; + std::function get_value; + std::function at_low; + std::function at_upper; + std::function lower_bound; + std::function upper_bound; + std::function column_lower_bound_constraint; + std::function column_upper_bound_constraint; + + bool is_real(unsigned) { return false; } // todo: test real case + void print_row(std::ostream& out, vector> & row ) { + bool first = true; + for (const auto & it : row) { + auto val = it.first; + if (first) { + first = false; + } else { + if (numeric_traits::is_pos(val)) { + out << " + "; + } else { + out << " - "; + val = -val; + } + } + if (val == -numeric_traits::one()) + out << " - "; + else if (val != numeric_traits::one()) + out << T_to_string(val); + + out << m_name_function(it.second); + } + } + + void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl, const mpq& f_0, const mpq& one_minus_f_0) { + TRACE("gomory_cut_detail_real", tout << "real\n";); + mpq new_a; + if (at_low(x_j)) { + if (a.is_pos()) { + new_a = a / (1 - f_0); + } + else { + new_a = a / f_0; + new_a.neg(); + } + k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation than + // k += lower_bound(x_j).x * new_a; + expl.push_justification(column_lower_bound_constraint(x_j), new_a); + } + else { + lp_assert(at_upper(x_j)); + if (a.is_pos()) { + new_a = a / f_0; + new_a.neg(); // the upper terms are inverted. + } + else { + new_a = a / (mpq(1) - f_0); + } + k.addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a; + expl.push_justification(column_upper_bound_constraint(x_j), new_a); + } + TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << k << "\n";); + pol.add_monomial(new_a, x_j); + } + + void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & t, explanation& expl, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) { + lp_assert(is_int(x_j)); + lp_assert(!a.is_int()); + lp_assert(f_0 > zero_of_type() && f_0 < one_of_type()); + mpq f_j = int_solver::fractional_part(a); + TRACE("gomory_cut_detail", + tout << a << " x_j = " << x_j << ", k = " << k << "\n"; + tout << "f_j: " << f_j << "\n"; + tout << "f_0: " << f_0 << "\n"; + tout << "1 - f_0: " << one_minus_f_0 << "\n"; + tout << "at_low(" << x_j << ") = " << at_low(x_j) << std::endl; + ); + lp_assert (!f_j.is_zero()); + mpq new_a; + if (at_low(x_j)) { + if (f_j <= one_minus_f_0) { + new_a = f_j / one_minus_f_0; + } + else { + new_a = (1 - f_j) / f_0; + } + k.addmul(new_a, lower_bound(x_j).x); + expl.push_justification(column_lower_bound_constraint(x_j), new_a); + } + else { + lp_assert(at_upper(x_j)); + if (f_j <= f_0) { + new_a = f_j / f_0; + } + else { + new_a = (mpq(1) - f_j) / (one_minus_f_0); + } + new_a.neg(); // the upper terms are inverted + k.addmul(new_a, upper_bound(x_j).x); + expl.push_justification(column_upper_bound_constraint(x_j), new_a); + } + TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << k << "\n";); + t.add_monomial(new_a, x_j); + lcm_den = lcm(lcm_den, denominator(new_a)); + } + + + void report_conflict_from_gomory_cut(mpq &k) { + lp_assert(false); + } + + void adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq &lcm_den) { + lp_assert(!t.is_empty()); + auto pol = t.coeffs_as_vector(); + t.clear(); + if (pol.size() == 1) { + TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;); + unsigned v = pol[0].second; + lp_assert(is_int(v)); + const mpq& a = pol[0].first; + k /= a; + if (a.is_pos()) { // we have av >= k + if (!k.is_int()) + k = ceil(k); + // switch size + t.add_monomial(- mpq(1), v); + k.neg(); + } else { + if (!k.is_int()) + k = floor(k); + t.add_monomial(mpq(1), v); + } + } else { + TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;); + lcm_den = lcm(lcm_den, denominator(k)); + TRACE("gomory_cut_detail", tout << "k: " << k << " lcm_den: " << lcm_den << "\n"; + for (unsigned i = 0; i < pol.size(); i++) { + tout << pol[i].first << " " << pol[i].second << "\n"; + } + tout << "k: " << k << "\n";); + lp_assert(lcm_den.is_pos()); + if (!lcm_den.is_one()) { + // normalize coefficients of integer parameters to be integers. + for (auto & pi: pol) { + pi.first *= lcm_den; + SASSERT(!is_int(pi.second) || pi.first.is_int()); + } + k *= lcm_den; + } + TRACE("gomory_cut_detail", tout << "after *lcm\n"; + for (unsigned i = 0; i < pol.size(); i++) { + tout << pol[i].first << " * v" << pol[i].second << "\n"; + } + tout << "k: " << k << "\n";); + + // negate everything to return -pol <= -k + for (const auto & pi: pol) + t.add_monomial(-pi.first, pi.second); + k.neg(); + } + TRACE("gomory_cut_detail", tout << "k = " << k << std::endl;); + lp_assert(k.is_int()); + } + + void print_term(lar_term & t, std::ostream & out) { + lp_assert(is_zero(t.m_v)); + vector> row; + for (auto p : t.m_coeffs) + row.push_back(std::make_pair(p.second, p.first)); + print_row(out, row); + } + + void mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsigned inf_col, vector> & row) { + enable_trace("gomory_cut"); + enable_trace("gomory_cut_detail"); + + TRACE("gomory_cut", + tout << "applying cut at:\n"; print_row(tout, row); + tout << std::endl << "inf_col = " << inf_col << std::endl; + ); + + // gomory will be t >= k + k = 1; + mpq lcm_den(1); + unsigned x_j; + mpq a; + bool some_int_columns = false; + mpq f_0 = int_solver::fractional_part(get_value(inf_col)); + mpq one_min_f_0 = 1 - f_0; + for ( auto pp : row) { + a = pp.first; + x_j = pp.second; + if (x_j == inf_col) + continue; + // make the format compatible with the format used in: Integrating Simplex with DPLL(T) + a.neg(); + if (is_real(x_j)) + real_case_in_gomory_cut(a, x_j, k, t, expl, f_0, one_min_f_0); + else { + if (a.is_int()) continue; // f_j will be zero and no monomial will be added + some_int_columns = true; + int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den, f_0, one_min_f_0); + } + } + + if (t.is_empty()) + return report_conflict_from_gomory_cut(k); + if (some_int_columns) + adjust_term_and_k_for_some_ints_case_gomory(t, k, lcm_den); + + TRACE("gomory_cut", tout<<"new cut :"; print_term(t, tout); tout << " >= " << k << std::endl;); + + } +}; +} diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 4cc36b46a..3e6e93f62 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -50,6 +50,8 @@ Revision History: #include "util/lp/integer_domain.h" #include "util/lp/stacked_map.h" #include +#include "test/lp/gomory_test.h" + namespace lp { unsigned seed = 1; @@ -1883,6 +1885,7 @@ void test_replace_column() { void setup_args_parser(argument_parser & parser) { + parser.add_option_with_help_string("-gomory", "gomory"); parser.add_option_with_help_string("-intd", "test integer_domain"); parser.add_option_with_help_string("-xyz_sample", "run a small interactive scenario"); parser.add_option_with_after_string_with_help("--density", "the percentage of non-zeroes in the matrix below which it is not dense"); @@ -3255,6 +3258,153 @@ void test_resolve(cut_solver& cs, unsigned constraint_index, unsigned i0) { } +void test_gomory_cut_0() { + gomory_test g( + [](unsigned j) { return "v" + T_to_string(j);} // name_function_p + , + [](unsigned j) { //get_value_p + if (j == 1) + return mpq(2730, 1727); + if (j == 2) + return zero_of_type(); + if (j == 3) return mpq(3); + lp_assert(false); + return zero_of_type(); + }, + [](unsigned j) { // at_low_p + if (j == 1) + return false; + if (j == 2) + return true; + if (j == 3) + return true; + lp_assert(false); + return false; + }, + [](unsigned j) { // at_upper + if (j == 1) + return false; + if (j == 2) + return true; + if (j == 3) + return false; + lp_assert(false); + return false; + }, + [](unsigned j) { // lower_bound + if (j == 1) { + lp_assert(false); //unlimited from below + return 0; + } + if (j == 2) + return 0; + if (j == 3) + return 3; + lp_assert(false); + return 0; + }, + [](unsigned j) { // upper + if (j == 1) { + lp_assert(false); //unlimited from above + return 0; + } + if (j == 2) + return 0; + if (j == 3) + return 10; + lp_assert(false); + return 0; + }, + [] (unsigned) { return 0; }, + [] (unsigned) { return 0; } + ); + lar_term t; + mpq k; + explanation expl; + unsigned inf_col = 1; + vector> row; + row.push_back(std::make_pair(mpq(1), 1)); + row.push_back(std::make_pair(mpq(2731, 1727), 2)); + row.push_back(std::make_pair(mpq(-910, 1727), 3)); + g.mk_gomory_cut(t, k, expl, inf_col, row); +} + +void test_gomory_cut_1() { + gomory_test g( + [](unsigned j) { return "v" + T_to_string(j);} // name_function_p + , + [](unsigned j) { //get_value_p + if (j == 1) + return mpq(-2); + if (j == 2) + return mpq(4363334, 2730001); + if (j == 3) + return mpq(1); + lp_assert(false); + return zero_of_type(); + }, + [](unsigned j) { // at_low_p + if (j == 1) + return false; + if (j == 2) + return false; + if (j == 3) + return true; + lp_assert(false); + return false; + }, + [](unsigned j) { // at_upper + if (j == 1) + return true; + if (j == 2) + return false; + if (j == 3) + return true; + lp_assert(false); + return false; + }, + [](unsigned j) { // lower_bound + if (j == 1) { + lp_assert(false); //unlimited from below + return 0; + } + if (j == 2) + return 1; + if (j == 3) + return 1; + lp_assert(false); + return 0; + }, + [](unsigned j) { // upper + if (j == 1) { + return -2; + } + if (j == 2) + return 3333; + if (j == 3) + return 10000; + lp_assert(false); + return 0; + }, + [] (unsigned) { return 0; }, + [] (unsigned) { return 0; } + ); + lar_term t; + mpq k; + explanation expl; + unsigned inf_col = 2; + vector> row; + row.push_back(std::make_pair(mpq(1726667, 2730001), 1)); + row.push_back(std::make_pair(mpq(-910000, 2730001), 3)); + row.push_back(std::make_pair(mpq(1), 2)); + g.mk_gomory_cut(t, k, expl, inf_col, row); +} + +void test_gomory_cut() { + test_gomory_cut_0(); + test_gomory_cut_1(); +} + void test_lp_local(int argn, char**argv) { // initialize_util_module(); @@ -3270,7 +3420,10 @@ void test_lp_local(int argn, char**argv) { } args_parser.print(); - + if (args_parser.option_is_used("-gomory")) { + test_gomory_cut(); + return finalize(0); + } if (args_parser.option_is_used("-intd")) { test_integer_domain(); return finalize(0); diff --git a/src/util/lp/active_set.h b/src/util/lp/active_set.h new file mode 100644 index 000000000..587570559 --- /dev/null +++ b/src/util/lp/active_set.h @@ -0,0 +1,76 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: + + +--*/ + +#pragma once +#include "util/lp/binary_heap_priority_queue.h" +namespace lp { +class active_set { + std::unordered_set m_cs; + binary_heap_priority_queue m_q; + std::unordered_map m_id_to_constraint; +public: + std::unordered_set cs() const { return m_cs;} + + bool contains(const constraint* c) const { + return m_id_to_constraint.find(c->id()) != m_id_to_constraint.end(); + } + + bool is_empty() const { return m_cs.size() == 0; } + // low priority will be dequeued first + void add_constraint(constraint* c, int priority) { + if (contains(c)) + return; + m_cs.insert(c); + m_id_to_constraint[c->id()] = c; + m_q.enqueue(c->id(), priority); + } + + void clear() { + m_cs.clear(); + m_id_to_constraint.clear(); + m_q.clear(); + } + + + constraint* remove_constraint() { + if (m_cs.size() == 0) + return nullptr; + unsigned id = m_q.dequeue(); + auto it = m_id_to_constraint.find(id); + lp_assert(it != m_id_to_constraint.end()); + constraint* c = it->second; + m_cs.erase(c); + m_id_to_constraint.erase(it); + return c; + } + + unsigned size() const { + return static_cast(m_cs.size()); + } + + void remove_constraint(constraint * c) { + if (! contains(c)) return; + + m_cs.erase(c); + m_id_to_constraint.erase(c->id()); + m_q.remove(c->id()); + } +}; +} diff --git a/src/util/lp/binary_heap_priority_queue.h b/src/util/lp/binary_heap_priority_queue.h index f36373029..9a71fc01e 100644 --- a/src/util/lp/binary_heap_priority_queue.h +++ b/src/util/lp/binary_heap_priority_queue.h @@ -45,7 +45,7 @@ public: unsigned size() const { return m_heap_size; } binary_heap_priority_queue(): m_heap(1), m_heap_size(0) {} // the empty constructror // n is the initial queue capacity. - // The capacity will be enlarged two times automatically if needed + // The capacity will be enlarged each time twice if needed binary_heap_priority_queue(unsigned n); void clear() { diff --git a/src/util/lp/binary_heap_priority_queue_def.h b/src/util/lp/binary_heap_priority_queue_def.h index 5178aaf3e..0d70e5aca 100644 --- a/src/util/lp/binary_heap_priority_queue_def.h +++ b/src/util/lp/binary_heap_priority_queue_def.h @@ -130,8 +130,12 @@ template void binary_heap_priority_queue::enqueue_new(unsigned o // In this case the priority will be changed and the queue adjusted. template void binary_heap_priority_queue::enqueue(unsigned o, const T & priority) { if (o >= m_priorities.size()) { - resize(o << 1); // make the size twice larger + if (o == 0) + resize(2); + else + resize(o << 1); // make the size twice larger } + if (m_heap_inverse[o] == -1) enqueue_new(o, priority); else diff --git a/src/util/lp/bound_analyzer_on_row.h b/src/util/lp/bound_analyzer_on_row.h index 2c421da72..6df372426 100644 --- a/src/util/lp/bound_analyzer_on_row.h +++ b/src/util/lp/bound_analyzer_on_row.h @@ -19,7 +19,6 @@ Revision History: --*/ #pragma once #include "util/vector.h" -#include "util/lp/linear_combination_iterator.h" #include "implied_bound.h" #include "test_bound_analyzer.h" #include "util/lp/bound_propagator.h" @@ -28,27 +27,85 @@ Revision History: // In a loop we drive the partial sum down, denoting the variables of this process by _u. // In the same loop trying to pin variables by pushing the partial sum up, denoting the variable related to it by _l namespace lp { - +template // C plays a role of a container class bound_analyzer_on_row { - - linear_combination_iterator & m_it; + struct term_with_basis_col { + const C & m_row; + unsigned m_bj; + struct ival { + unsigned m_var; + const mpq & m_coeff; + ival(unsigned var, const mpq & val) : m_var(var), m_coeff(val) { + } + unsigned var() const { return m_var;} + const mpq & coeff() const { return m_coeff; } + }; + + term_with_basis_col(const C& row, unsigned bj) : m_row(row), m_bj(bj) {} + struct const_iterator { + // fields + typename C::const_iterator m_it; + unsigned m_bj; + + + //typedefs + + + typedef const_iterator self_type; + typedef ival value_type; + typedef ival reference; + typedef int difference_type; + typedef std::forward_iterator_tag iterator_category; + + reference operator*() const { + if (m_bj == static_cast(-1)) + return ival((*m_it).var(), (*m_it).coeff()); + return ival(m_bj, - 1); + } + self_type operator++() { self_type i = *this; operator++(1); return i; } + + self_type operator++(int) { + if (m_bj == static_cast(-1)) + m_it++; + else + m_bj = static_cast(-1); + return *this; + } + + // constructor + const_iterator(const typename C::const_iterator& it, unsigned bj) : + m_it(it), + m_bj(bj) + {} + bool operator==(const self_type &other) const { + return m_it == other.m_it && m_bj == other.m_bj ; + } + bool operator!=(const self_type &other) const { return !(*this == other); } + }; + const_iterator begin() const { + return const_iterator( m_row.begin(), m_bj); + } + const_iterator end() const { return const_iterator(m_row.end(), m_bj); } + }; + term_with_basis_col m_row; bound_propagator & m_bp; - unsigned m_row_or_term_index; - int m_column_of_u; // index of an unlimited from above monoid + unsigned m_row_or_term_index; + int m_column_of_u; // index of an unlimited from above monoid // -1 means that such a value is not found, -2 means that at least two of such monoids were found - int m_column_of_l; // index of an unlimited from below monoid - impq m_rs; + int m_column_of_l; // index of an unlimited from below monoid + impq m_rs; public : // constructor bound_analyzer_on_row( - linear_combination_iterator &it, + const C & it, + unsigned bj, // basis column for the row const numeric_pair& rs, unsigned row_or_term_index, bound_propagator & bp ) : - m_it(it), + m_row(it, bj), m_bp(bp), m_row_or_term_index(row_or_term_index), m_column_of_u(-1), @@ -59,11 +116,11 @@ public : unsigned j; void analyze() { - - mpq a; unsigned j; - while (((m_column_of_l != -2) || (m_column_of_u != -2)) && m_it.next(a, j)) - analyze_bound_on_var_on_coeff(j, a); - + for (auto c : m_row) { + if ((m_column_of_l == -2) && (m_column_of_u == -2)) + break; + analyze_bound_on_var_on_coeff(c.var(), c.coeff()); + } if (m_column_of_u >= 0) limit_monoid_u_from_below(); else if (m_column_of_u == -1) @@ -168,25 +225,23 @@ public : int strict = 0; mpq total; lp_assert(is_zero(total)); - m_it.reset(); - mpq a; unsigned j; - while (m_it.next(a, j)) { + for (auto p : m_row) { bool str; - total -= monoid_min(a, j, str); + total -= monoid_min(p.coeff(), p.var(), str); if (str) strict++; } - m_it.reset(); - while (m_it.next(a, j)) { + + for (auto p : m_row) { bool str; - bool a_is_pos = is_pos(a); - mpq bound = total / a + monoid_min_no_mult(a_is_pos, j, str); + bool a_is_pos = is_pos(p.coeff()); + mpq bound = total / p.coeff() + monoid_min_no_mult(a_is_pos, p.var(), str); if (a_is_pos) { - limit_j(j, bound, true, false, strict - static_cast(str) > 0); + limit_j(p.var(), bound, true, false, strict - static_cast(str) > 0); } else { - limit_j(j, bound, false, true, strict - static_cast(str) > 0); + limit_j(p.var(), bound, false, true, strict - static_cast(str) > 0); } } } @@ -195,25 +250,23 @@ public : int strict = 0; mpq total; lp_assert(is_zero(total)); - m_it.reset(); - mpq a; unsigned j; - while (m_it.next(a, j)) { + for (auto p : m_row) { bool str; - total -= monoid_max(a, j, str); + total -= monoid_max(p.coeff(), p.var(), str); if (str) strict++; } - m_it.reset(); - while (m_it.next(a, j)) { + + for (auto p : m_row) { bool str; - bool a_is_pos = is_pos(a); - mpq bound = total / a + monoid_max_no_mult(a_is_pos, j, str); + bool a_is_pos = is_pos(p.coeff()); + mpq bound = total / p.coeff() + monoid_max_no_mult(a_is_pos, p.var(), str); bool astrict = strict - static_cast(str) > 0; if (a_is_pos) { - limit_j(j, bound, true, true, astrict); + limit_j(p.var(), bound, true, true, astrict); } else { - limit_j(j, bound, false, false, astrict); + limit_j(p.var(), bound, false, false, astrict); } } } @@ -222,18 +275,18 @@ public : void limit_monoid_u_from_below() { // we are going to limit from below the monoid m_column_of_u, // every other monoid is impossible to limit from below - mpq u_coeff, a; + mpq u_coeff; unsigned j; mpq bound = -m_rs.x; - m_it.reset(); bool strict = false; - while (m_it.next(a, j)) { + for (auto p : m_row) { + j = p.var(); if (j == static_cast(m_column_of_u)) { - u_coeff = a; + u_coeff = p.coeff(); continue; } bool str; - bound -= monoid_max(a, j, str); + bound -= monoid_max(p.coeff(), j, str); if (str) strict = true; } @@ -251,19 +304,19 @@ public : void limit_monoid_l_from_above() { // we are going to limit from above the monoid m_column_of_l, // every other monoid is impossible to limit from above - mpq l_coeff, a; + mpq l_coeff; unsigned j; mpq bound = -m_rs.x; bool strict = false; - m_it.reset(); - while (m_it.next(a, j)) { + for (auto p : m_row) { + j = p.var(); if (j == static_cast(m_column_of_l)) { - l_coeff = a; + l_coeff = p.coeff(); continue; } bool str; - bound -= monoid_min(a, j, str); + bound -= monoid_min(p.coeff(), j, str); if (str) strict = true; } @@ -282,7 +335,7 @@ public : // bool lower_bound = be.m_lower_bound; // if (!coeff_is_pos) // lower_bound = !lower_bound; - // auto it = m_it.clone(); + // auto it = m_row.clone(); // mpq a; unsigned j; // while (it->next(a, j)) { // if (be.m_j == j) continue; @@ -336,14 +389,16 @@ public : } } - static void analyze_row(linear_combination_iterator &it, + static void analyze_row(const C & row, + unsigned bj, // basis column for the row const numeric_pair& rs, unsigned row_or_term_index, bound_propagator & bp ) { - bound_analyzer_on_row a(it, rs, row_or_term_index, bp); + bound_analyzer_on_row a(row, bj, rs, row_or_term_index, bp); a.analyze(); } }; } + diff --git a/src/util/lp/column_namer.h b/src/util/lp/column_namer.h index 7daf6676b..3eeff0e26 100644 --- a/src/util/lp/column_namer.h +++ b/src/util/lp/column_namer.h @@ -19,31 +19,19 @@ Revision History: --*/ #include -#include "util/lp/linear_combination_iterator.h" +#include "util/lp/static_matrix.h" namespace lp { class column_namer { public: virtual std::string get_column_name(unsigned j) const = 0; template - void print_linear_iterator(linear_combination_iterator* it, std::ostream & out) const { + void print_row(const row_strip & row, std::ostream & out) const { vector> coeff; - T a; - unsigned i; - while (it->next(a, i)) { - coeff.push_back(std::make_pair(a, i)); + for (auto p : row) { + coeff.push_back(std::make_pair(p.coeff(), p.var())); } print_linear_combination_of_column_indices(coeff, out); } - template - void print_linear_iterator_indices_only(linear_combination_iterator* it, std::ostream & out) const { - vector> coeff; - T a; - unsigned i; - while (it->next(a, i)) { - coeff.push_back(std::make_pair(a, i)); - } - print_linear_combination_of_column_indices_only(coeff, out); - } template void print_linear_combination_of_column_indices_only(const vector> & coeffs, std::ostream & out) const { diff --git a/src/util/lp/constraint.h b/src/util/lp/constraint.h new file mode 100644 index 000000000..84ec188d6 --- /dev/null +++ b/src/util/lp/constraint.h @@ -0,0 +1,99 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: + + +--*/ + +#pragma once +namespace lp { +class constraint; // forward definition +struct constraint_hash { + size_t operator() (const constraint* c) const; +}; + +struct constraint_equal { + bool operator() (const constraint * a, const constraint * b) const; +}; + +class constraint { // we only have less or equal for the inequality sign, which is enough for integral variables + int m_id; + bool m_is_ineq; + polynomial m_poly; + mpq m_d; // the divider for the case of a divisibility constraint + std::unordered_set m_assert_origins; // these indices come from the client and get collected during tightening +public : + unsigned id() const { return m_id; } + const polynomial & poly() const { return m_poly; } + polynomial & poly() { return m_poly; } + std::unordered_set & assert_origins() { return m_assert_origins;} + const std::unordered_set & assert_origins() const { return m_assert_origins;} + bool is_lemma() const { return !is_assert(); } + bool is_assert() const { return m_assert_origins.size() == 1; } + bool is_ineq() const { return m_is_ineq; } + const mpq & divider() const { return m_d; } +public: + constraint( + unsigned id, + constraint_index assert_origin, + const polynomial & p, + bool is_ineq): + m_id(id), + m_is_ineq(is_ineq), + m_poly(p) + { // creates an assert + m_assert_origins.insert(assert_origin); + } + constraint( + unsigned id, + const std::unordered_set& origins, + const polynomial & p, + bool is_ineq): + m_id(id), + m_is_ineq(is_ineq), + m_poly(p), + m_assert_origins(origins) + {} + + + + constraint( + unsigned id, + const polynomial & p, + bool is_ineq): + m_id(id), + m_is_ineq(is_ineq), + m_poly(p) { // creates a lemma + } + +public: + constraint() {} + + const mpq & coeff(var_index j) const { + return m_poly.coeff(j); + } + const vector& coeffs() const { return m_poly.m_coeffs;} + + bool is_tight(unsigned j) const { + const mpq & a = m_poly.coeff(j); + return a == 1 || a == -1; + } + void add_predecessor(const constraint* p) { + lp_assert(p != nullptr); + for (auto m : p->assert_origins()) + m_assert_origins.insert(m); } +}; +} diff --git a/src/util/lp/core_solver_pretty_printer_def.h b/src/util/lp/core_solver_pretty_printer_def.h index c8dac1b03..58ffbb481 100644 --- a/src/util/lp/core_solver_pretty_printer_def.h +++ b/src/util/lp/core_solver_pretty_printer_def.h @@ -105,6 +105,8 @@ template void core_solver_pretty_printer::init_m_ string name = m_core_solver.column_name(column); for (unsigned row = 0; row < nrows(); row ++) { + m_A[row].resize(ncols(), ""); + m_signs[row].resize(ncols(),""); set_coeff( m_A[row], m_signs[row], diff --git a/src/util/lp/cut_solver.cpp b/src/util/lp/cut_solver.cpp index 241759cb3..2bfc7c229 100644 --- a/src/util/lp/cut_solver.cpp +++ b/src/util/lp/cut_solver.cpp @@ -19,6 +19,5 @@ namespace lp { c.s.print_constraint(out, c.c); return out; } - } diff --git a/src/util/lp/cut_solver.h b/src/util/lp/cut_solver.h index 6a11da117..e39244006 100644 --- a/src/util/lp/cut_solver.h +++ b/src/util/lp/cut_solver.h @@ -1,7 +1,23 @@ -/* +/*++ Copyright (c) 2017 Microsoft Corporation - Author: Nikolaj Bjorner, Lev Nachmanson -*/ + + Module Name: + + + + Abstract: + + + + Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + + Revision History: + + + --*/ + #pragma once #include "util/vector.h" #include "util/trace.h" @@ -11,223 +27,15 @@ #include "util/lp/lp_utils.h" #include #include "util/lp/int_set.h" -#include "util/lp/linear_combination_iterator_on_std_vector.h" #include "util/lp/stacked_vector.h" +#include "util/lp/monomial.h" +#include "util/lp/polynomial.h" +#include "util/lp/constraint.h" +#include "util/lp/active_set.h" +#include "util/lp/indexer_of_constraints.h" namespace lp { -enum class lbool { l_false, l_true, l_undef }; -inline std::string lbool_to_string(lbool l) { - switch(l) { - case lbool::l_true: return std::string("true"); - case lbool::l_false: return std::string("false"); - case lbool::l_undef: return std::string("undef"); - default: - return std::string("what is it?"); - } -} - -class cut_solver; - -struct monomial { - mpq m_coeff; // the coefficient of the monomial - var_index m_var; // the variable index -public: - monomial(const mpq& coeff, var_index var) : m_coeff(coeff), m_var(var) {} - // copy constructor - monomial(const monomial& m) : monomial(m.coeff(), m.var()) {} - monomial(var_index var) : monomial(one_of_type(), var) {} - const mpq & coeff() const { return m_coeff; } - mpq & coeff() { return m_coeff; } - var_index var() const { return m_var; } - std::pair to_pair() const { return std::make_pair(coeff(), var());} -}; - -struct polynomial { - static mpq m_local_zero; - // the polynomial evaluates to m_coeffs + m_a - vector m_coeffs; - mpq m_a; // the free coefficient - polynomial(const vector& p, const mpq & a) : m_coeffs(p), m_a(a) {} - polynomial(const vector& p) : polynomial(p, zero_of_type()) {} - polynomial(): m_a(zero_of_type()) {} - polynomial(const polynomial & p) : m_coeffs(p.m_coeffs), m_a(p.m_a) {} - - const mpq & coeff(var_index j) const { - for (const auto & t : m_coeffs) { - if (j == t.var()) { - return t.coeff(); - } - } - return m_local_zero; - } - - polynomial & operator+=(const polynomial & p) { - m_a += p.m_a; - for (const auto & t: p.m_coeffs) - *this += monomial(t.coeff(), t.var()); - return *this; - } - - void add(const mpq & c, const polynomial &p) { - m_a += p.m_a * c; - - for (const auto & t: p.m_coeffs) - *this += monomial(c * t.coeff(), t.var()); - } - - void clear() { - m_coeffs.clear(); - m_a = zero_of_type(); - } - - bool is_empty() const { return m_coeffs.size() == 0 && numeric_traits::is_zero(m_a); } - - unsigned number_of_monomials() const { return m_coeffs.size();} - - void add(const monomial &m ){ - if (is_zero(m.coeff())) return; - for (unsigned k = 0; k < m_coeffs.size(); k++) { - auto & l = m_coeffs[k]; - if (m.var() == l.var()) { - l.coeff() += m.coeff(); - if (l.coeff() == 0) - m_coeffs.erase(m_coeffs.begin() + k); - return; - } - } - m_coeffs.push_back(m); - lp_assert(is_correct()); - } - - polynomial & operator+=(const monomial &m ){ - add(m); - return *this; - } - - polynomial & operator+=(const mpq &c ){ - m_a += c; - return *this; - } - - - bool is_correct() const { - std::unordered_set s; - for (auto & l : m_coeffs) { - if (l.coeff() == 0) - return false; - s.insert(l.var()); - } - return m_coeffs.size() == s.size(); - } - - bool var_coeff_is_unit(unsigned j) const { - const mpq & a = coeff(j); - return a == 1 || a == -1; - } - const vector & coeffs() const { return m_coeffs; } -}; - -class constraint; // forward definition -struct constraint_hash { - size_t operator() (const constraint* c) const; -}; - -struct constraint_equal { - bool operator() (const constraint * a, const constraint * b) const; -}; - -class constraint { // we only have less or equal for the inequality sign, which is enough for integral variables - int m_id; - bool m_is_ineq; - polynomial m_poly; - mpq m_d; // the divider for the case of a divisibility constraint - const svector m_assert_origins; // these indices come from the client - bool m_active; - std::unordered_set m_predecessors; // used in tight_constraints and lemmas -public : - void set_active_flag() {m_active = true;} - void remove_active_flag() { m_active = false; } - bool is_active() const { return m_active; } - unsigned id() const { return m_id; } - const polynomial & poly() const { return m_poly; } - polynomial & poly() { return m_poly; } - const svector & assert_origins() const { return m_assert_origins;} - bool is_lemma() const { return !is_assert(); } - bool is_assert() const { return m_predecessors.size() == 0; } - bool is_ineq() const { return m_is_ineq; } - const mpq & divider() const { return m_d; } -public : - static constraint * make_ineq_assert( - int id, - const vector& term, - const mpq& a, - const svector & origins) { - return new constraint(id, origins, polynomial(term, a), true); - } - - static constraint * make_ineq_constraint( - int id, - const polynomial & p, - std::unordered_set predecessors_set) { - auto c = new constraint(id, p, true); - c->predecessors() = predecessors_set; - return c; - } - - static constraint * make_ineq_lemma(unsigned id, const polynomial &p) { - return new constraint(id, p, true); - } - - // static constraint make_div_lemma(unsigned id, const polynomial &p, const mpq & div) { - // lp_assert(false); // not implemented - // // constraint * c = new constraint(id, p, true); - // // c->m_d = div; - // return nullptr; - // } -private: - constraint( - unsigned id, - const svector & assert_origins, - const polynomial & p, - bool is_ineq): - m_id(id), - m_is_ineq(is_ineq), - m_poly(p), - m_assert_origins(assert_origins), - m_active(false) { // creates an assert - } - - - - constraint( - unsigned id, - const polynomial & p, - bool is_ineq): - m_id(id), - m_is_ineq(is_ineq), - m_poly(p), - m_active(false) { // creates a lemma - } - -public: - constraint() {} - - const mpq & coeff(var_index j) const { - return m_poly.coeff(j); - } - const vector& coeffs() const { return m_poly.m_coeffs;} - - bool is_tight(unsigned j) const { - const mpq & a = m_poly.coeff(j); - return a == 1 || a == -1; - } - const std::unordered_set& predecessors() const { return m_predecessors; } - std::unordered_set& predecessors() { return m_predecessors; } - void add_predecessor(const constraint* p) { - lp_assert(p != nullptr); - m_predecessors.insert(p); } -}; - +class cut_solver; //forward definition struct pp_poly { cut_solver const& s; @@ -245,8 +53,18 @@ std::ostream& operator<<(std::ostream& out, pp_poly const& p); std::ostream& operator<<(std::ostream& out, pp_constraint const& p); class cut_solver : public column_namer { -public: // for debugging - +public: + enum class lbool { l_false, l_true, l_undef }; + inline std::string lbool_to_string(lbool l) { + switch(l) { + case lbool::l_true: return std::string("true"); + case lbool::l_false: return std::string("false"); + case lbool::l_undef: return std::string("undef"); + default: + return std::string("what is it?"); + } + } + typedef lp::polynomial polynomial; typedef lp::monomial monomial; @@ -334,7 +152,7 @@ public: // for debugging enum class bound_type { LOWER, UPPER, UNDEF - }; + }; struct bound_result { mpq m_bound; @@ -363,7 +181,7 @@ public: // for debugging unsigned m_internal_j; // it is just the index into m_var_infos of this var_info, if the var_info is not active then this value is set to (unsigned)-1 integer_domain m_domain; // the map of constraints using the var: bound_type = UNDEF stands for a div constraint - std::unordered_map m_dependent_constraints; + std::unordered_map m_dependent_constraints; // pair (c, true) means that c has a monomial (a, m_internal_j) for a > 0, and (c, false) means that a < 0 unsigned m_external_stack_level; unsigned m_number_of_bound_propagations; unsigned m_number_of_asserts; @@ -391,11 +209,11 @@ public: // for debugging void activate(unsigned internal_j) { m_internal_j = internal_j; } - void add_dependent_constraint(constraint* i, bound_type bt) { + void add_dependent_constraint(constraint* i, bool coeff_is_pos) { lp_assert(m_dependent_constraints.find(i) == m_dependent_constraints.end()); if (i->is_assert()) m_number_of_asserts++; - m_dependent_constraints[i] = bt; + m_dependent_constraints[i] = coeff_is_pos; } void remove_depended_constraint(constraint* i) { lp_assert(m_dependent_constraints.find(i) != m_dependent_constraints.end()); @@ -413,12 +231,16 @@ public: // for debugging bool intersect_with_lower_bound(const mpq & b, unsigned explanation, unsigned stack_level) { conditional_push(stack_level); - return m_domain.intersect_with_bound(b, true, explanation); + bool ret = m_domain.intersect_with_bound(b, true, explanation); + lp_assert(!m_domain.is_empty()); + return ret; } bool intersect_with_upper_bound(const mpq & b, unsigned explanation, unsigned external_level) { conditional_push(external_level); - return m_domain.intersect_with_bound(b, false, explanation); + bool ret = m_domain.intersect_with_bound(b, false, explanation); + lp_assert(!m_domain.is_empty()); + return ret; } bool is_fixed() const { return m_domain.is_fixed();} @@ -427,8 +249,8 @@ public: // for debugging bool get_upper_bound_with_expl(mpq & b, unsigned& expl) const { return m_domain.get_upper_bound_with_expl(b, expl); } bool get_lower_bound_with_expl(mpq & b, unsigned& expl) const { return m_domain.get_lower_bound_with_expl(b, expl); } void print_var_domain(std::ostream &out) const { m_domain.print(out); } - std::unordered_map & dependent_constraints() { return m_dependent_constraints; } - const std::unordered_map & dependent_constraints() const { return m_dependent_constraints; } + std::unordered_map & dependent_constraints() { return m_dependent_constraints; } + const std::unordered_map & dependent_constraints() const { return m_dependent_constraints; } int get_lower_bound_expl() const { return m_domain.get_lower_bound_expl();} int get_upper_bound_expl() const { return m_domain.get_upper_bound_expl();} public: @@ -444,6 +266,29 @@ public: // for debugging }; // end of var_info + struct lemmas_container { + std::unordered_set m_lemmas; + unsigned size() const { return m_lemmas.size(); } + void submit_assert_origin_for_delete(constraint_index ci) { m_deleted_assert_origins.insert(ci);} + vector remove_lemmas_depending_on_submitted_origins() { + vector r; + for (auto * l : m_lemmas) { + for (unsigned o : l->assert_origins()) + if (m_deleted_assert_origins.find(o) != m_deleted_assert_origins.end()) { + r.push_back(l); + break; // from the inner loop + } + } + for (auto * l : r) + m_lemmas.erase(l); + + m_deleted_assert_origins.clear(); + return r; + } + std::unordered_set m_deleted_assert_origins; + void add_lemma(constraint* l) { m_lemmas.insert(l); } + }; + bool lhs_is_int(const vector & lhs) const { for (auto & p : lhs) { if (numeric_traits::is_int(p.coeff()) == false) return false; @@ -492,11 +337,11 @@ public: if (coeff == one_of_type() || coeff == - one_of_type()) continue; if (!all_fixed_except_j(c->poly(), j)) continue; - if (p.second == bound_type::UPPER) { + if (!p.second) { upper_bounds++; upper = c; if (lower_bounds) return true; - } else if (p.second == bound_type::LOWER) { + } else { lower_bounds++; lower = c; if (upper_bounds) @@ -514,71 +359,21 @@ public: ~cut_solver() { for (constraint * c : m_asserts) delete c; - for (constraint * c : m_lemmas) + for (constraint * c : m_lemmas_container.m_lemmas) delete c; } - class active_set { - std::unordered_set m_cs; - public: - - std::unordered_set cs() const { return m_cs;} - - bool is_empty() const { return m_cs.size() == 0; } - - void add_constraint(constraint* c) { - if (c->is_active()) return; - m_cs.insert(c); - c->set_active_flag(); - } - - void clear() { - for (constraint * c: m_cs) { - c->remove_active_flag(); - } - m_cs.clear(); - } - - - constraint* remove_random_constraint(unsigned rand) { - if (m_cs.size() == 0) - return nullptr; - unsigned j = rand % m_cs.size(); - auto it = std::next(m_cs.begin(), j); - constraint * c = *it; - c->remove_active_flag(); - m_cs.erase(it); - return c; - } - - unsigned size() const { - return static_cast(m_cs.size()); - } - - void remove_constraint(constraint * c) { - m_cs.erase(c); - c->remove_active_flag(); - } - }; - struct scope { - unsigned m_asserts_size; - unsigned m_lemmas_size; unsigned m_trail_size; scope() {} - scope(unsigned asserts_size, - unsigned lemmas_size, - unsigned trail_size) : m_asserts_size(asserts_size), - m_lemmas_size(lemmas_size), - m_trail_size(trail_size) - {} - + scope( unsigned trail_size) : + m_trail_size(trail_size) {} }; // fields - vector m_var_infos; - svector m_asserts; - svector m_lemmas; + vector m_var_infos; + vector m_asserts; + lemmas_container m_lemmas_container; vector m_v; // the values of the variables std::function m_var_name_function; std::function m_print_constraint_function; @@ -587,7 +382,6 @@ public: active_set m_active_set; vector m_trail; lp_settings & m_settings; - unsigned m_max_constraint_id; std::set m_U; // the set of conflicting cores unsigned m_bounded_search_calls; unsigned m_number_of_conflicts; @@ -600,9 +394,9 @@ public: bool m_cancelled; // debug fields unsigned m_number_of_constraints_tried_for_propagaton; - unsigned m_number_of_pops; + unsigned m_pushes_to_trail; + indexer_of_constraints m_indexer_of_constraints; - bool is_lower_bound(literal & l) const { return l.is_lower(); } @@ -641,10 +435,9 @@ public: bool check_inconsistent() { if (at_base_lvl()) { - // the last added lemmas can give the contradiction - for (unsigned j = m_lemmas.size(); j--; ) { - if (lower_is_pos(m_lemmas[j])) { - TRACE("check_inconsistent_int", tout << pp_poly(*this, m_lemmas[j]->poly()) << "\n";); + for (constraint *c : m_lemmas_container.m_lemmas ) { + if (lower_is_pos(c)) { + TRACE("check_inconsistent_int", tout << pp_poly(*this, c->poly()) << "\n";); lp_assert(false); // not implemented return true; } @@ -671,7 +464,7 @@ public: return false; } } - for (auto c : m_lemmas) { + for (auto c : m_lemmas_container.m_lemmas) { if (is_pos(value(c->poly()))) { TRACE("all_constraints_hold_int", tout << "constraint does not hold\n"; print_constraint(tout, *c);); @@ -694,11 +487,6 @@ public: return lbool::l_true; } - bool resolve_conflict_core() { - // this is where the main action is. - return true; - } - void find_new_conflicting_cores_under_constraint(var_index j, const_constr* c) { // lp_assert(false); } @@ -736,11 +524,12 @@ public: } // 'j' is a variable that changed - void add_changed_var(unsigned j) { + void add_changed_var(unsigned j, bool lower_bound_got_tighter, int priority) { TRACE("add_changed_var_int", tout << "j = " << j << "\n";); for (auto & p: m_var_infos[j].dependent_constraints()) { TRACE("add_changed_var_int", tout << pp_constraint(*this, *p.first) << "\n";); - m_active_set.add_constraint(p.first); + if (p.second == lower_bound_got_tighter) + m_active_set.add_constraint(p.first, 2); // 2 is a priority } } @@ -762,6 +551,38 @@ public: test_bound_struct() :m_expl_lower(-1), m_expl_upper(-1) {} }; + bool state_is_correct() const { + return var_infos_are_correct() && constraint_indices_are_correct(); + } + + bool constraint_indices_are_correct() const { + std::unordered_set r; + unsigned n = m_indexer_of_constraints.max(); + for (auto c: m_asserts) { + if (c->id() >= n) { + std::cout << "c->id >= n\n"; + return false; + } + if (r.find(c->id()) != r.end()){ + std::cout << "id found\n"; + return false; + } + r.insert(c->id()); + } + + for (auto c: m_lemmas_container.m_lemmas) { + if (c->id() >= n) { + std::cout << "c->id >= n\n"; + return false; + } + if (r.find(c->id()) != r.end()){ + std::cout << "id found\n"; + return false; + } + r.insert(c->id()); + } + return true; + } bool var_infos_are_correct() const { vector bounds = get_bounds_from_trail(); @@ -885,7 +706,7 @@ public: if (!is_zero(c->coeff(j))) deps.insert(c); } - for (const auto c: m_lemmas) { + for (const auto c: m_lemmas_container.m_lemmas) { if (!is_zero(c->coeff(j))) deps.insert(c); } @@ -972,8 +793,9 @@ public: } void push_literal_to_trail(literal & l) { + m_pushes_to_trail ++; m_trail.push_back(l); - add_changed_var(l.var()); + add_changed_var(l.var(), l.is_lower(), 0); } unsigned find_large_enough_j(unsigned i) { @@ -1022,12 +844,8 @@ public: } - bool cycle_on_var_info(const var_info & vi) const { - return vi.number_of_bound_propagations() > m_settings.m_cut_solver_bound_propagation_factor * vi.number_of_asserts(); - } - bool too_many_propagations_for_var(const var_info& vi) const { - return vi.number_of_bound_propagations() > m_settings.m_cut_solver_bound_propagation_factor * vi.dependent_constraints().size(); + return vi.number_of_bound_propagations() > m_settings.m_cut_solver_cycle_on_var * vi.number_of_asserts(); } bool new_upper_bound_is_relevant(unsigned j, const mpq & v) const { @@ -1067,10 +885,8 @@ public: const mpq& lower_val, constraint* c) { unsigned j = p.var(); - if (cycle_on_var_info(m_var_infos[j])) { - return; - } - + if (m_var_infos[j].is_fixed()) + return; if (is_pos(p.coeff())) { mpq m; get_var_lower_bound(p.var(), m); @@ -1094,9 +910,7 @@ public: const mpq& rs, constraint *c) { unsigned j = p.var(); - if (cycle_on_var_info(m_var_infos[j])) { - return; - } + if (is_pos(p.coeff())) { mpq v = floor(rs / p.coeff()); if (new_upper_bound_is_relevant(j, v)) { @@ -1122,7 +936,7 @@ public: } if (vi.is_active()) out << ", active var "; print_var_domain(out, vi); - out << ", propagaions = " << vi.number_of_bound_propagations() << ", deps = " << vi.dependent_constraints().size(); + out << ", propagations = " << vi.number_of_bound_propagations() << ", deps = " << vi.dependent_constraints().size(); out << ", asserts = " << vi.number_of_asserts() << std::endl; // out << "external levels: "; // for (auto j : vi.external_stack_level()) @@ -1195,9 +1009,6 @@ public: for (auto j : c->assert_origins()) m_explanation.insert(j); - - for (const_constr* p : c->predecessors()) - add_premises(p, visited_constraints); } void fill_conflict_explanation(const constraint *confl) { @@ -1223,12 +1034,9 @@ public: void trace_print_constraint(std::ostream& out, const_constr* i) const { print_constraint(out, *i); out << "id = " << i->id() << ", "; - unsigned j; - auto pairs = to_pairs(i->poly().m_coeffs); - auto it = linear_combination_iterator_on_vector(pairs); - while (it.next(j)) { - out << "domain of " << var_name(j) << " = "; - print_var_domain(out, j); + for (auto p : i->poly().m_coeffs){ + out << "domain of " << var_name(p.var()) << " = "; + print_var_domain(out, p.var()); } if (i->assert_origins().size()) { out << (i->assert_origins().size() > 1?"origins: ":"origin: ") ; @@ -1524,7 +1332,7 @@ public: if (i.is_ineq()) { out << " <= 0"; } - out << " active = " << i.is_active() << " "; + out << " active = " << m_active_set.contains(&i) << " "; mpq b; if (lower(&i, b)) { out << ", lower = " << b; @@ -1739,20 +1547,14 @@ public: } void print_state_stats(std::ostream & out) const { - static bool one_time_print = true; if (m_bounded_search_calls % 10 ) return; out << "search_calls = " << m_bounded_search_calls << ", "; out << "vars = " << m_var_infos.size() << ","; out << "asserts = "<< m_asserts.size() << ", "; - out << "lemmas = " << m_lemmas.size() << ", "; + out << "lemmas = " << m_lemmas_container.size() << ", "; out << "trail = " << m_trail.size() << ", props = " << m_number_of_constraints_tried_for_propagaton << ", "; out << "cnfls = " << m_number_of_conflicts << ", "; - out << "pops = " << m_number_of_pops << std::endl; - if (one_time_print && m_lemmas.size() >= 500) { - print_state(out); - one_time_print = false; - } } @@ -1800,12 +1602,10 @@ public: m_bounded_search_calls = 0; m_stuck_state = false; m_cancelled = false; - for (constraint * c : m_asserts) - m_active_set.add_constraint(c); for (auto & vi : m_var_infos) vi.number_of_bound_propagations() = 0; m_number_of_constraints_tried_for_propagaton = 0; - m_number_of_pops = 0; + m_pushes_to_trail = 0; } constraint* propagate_constraint(constraint* c) { @@ -1825,7 +1625,9 @@ public: return c; } propagate_constraint_on_lower(c, b); - } else if (r == 1) { + } else if (r == 1 && + !too_many_propagations_for_var(c->poly().coeffs()[the_only_unlim].var())) { + // otherwise the new bound, even if it exists, will be rejected! lp_assert(!lower_is_pos(c->poly())); propagate_constraint_only_one_unlim(c, the_only_unlim); } @@ -1849,8 +1651,8 @@ public: print_constraint(out, *i); } out << "end of constraints\n"; - out << "lemmas total " << m_lemmas.size() << "\n"; - for (const auto i: m_lemmas) { + out << "lemmas total " << m_lemmas_container.size() << "\n"; + for (const auto i: m_lemmas_container.m_lemmas) { print_constraint(out, *i); } out << "end of constraints\n"; @@ -1920,10 +1722,21 @@ public: return false; } + bool cancel_when_propagation_speed_is_too_slow() { + if (m_number_of_constraints_tried_for_propagaton >= 10000) { + m_cancelled = m_pushes_to_trail <= 50; + m_pushes_to_trail = 0; + m_number_of_constraints_tried_for_propagaton = 0; + return m_cancelled; + } + return false; + } lbool propagate_and_backjump_step() { do { + if (cancel_when_propagation_speed_is_too_slow()) + return lbool::l_undef; constraint* c = propagate(); if (cancel()) return lbool::l_undef; @@ -2144,12 +1957,14 @@ public: l.tight_constr() = l.cnstr(); return; } - l.tight_constr() = constraint::make_ineq_constraint( m_max_constraint_id++, - c->poly(), - c->predecessors()); + l.tight_constr() = new constraint( get_new_constraint_id(), + c->assert_origins(), + c->poly(), + true); l.tight_constr()->add_predecessor(c); tighten(l.tight_constr(), j, a, index_of_literal); add_lemma_as_not_active(l.tight_constr()); + lp_assert(constraint_indices_are_correct()); } void backjump(unsigned index_of_literal, @@ -2167,12 +1982,13 @@ public: print_var_info(tout, j); tout << "p = " << pp_poly(*this, p) << "\n";); if (!improves(j, br)) { - CTRACE("int_backjump", lp_settings::ddd, br.print(tout); tout << "\nimproves is false\n"; - tout << "var info after pop = "; print_var_info(tout, j);); - if (lemma_has_been_modified) - add_lemma(lemma); - else { - m_active_set.add_constraint(orig_conflict); + TRACE("int_backjump", br.print(tout); tout << "\nimproves is false\n"; + tout << "var info after pop = "; print_var_info(tout, j);); + if (lemma_has_been_modified) { // flip_coin() gives + // priority 0 or 1, so it is be dequeued fast + add_lemma(lemma, flip_coin()); + } else { + m_active_set.add_constraint(orig_conflict, flip_coin()); } } else { constraint *c; @@ -2290,19 +2106,19 @@ public: return b; } + unsigned get_new_constraint_id() { return m_indexer_of_constraints.get_new_index(); } void resolve_conflict_for_inequality(constraint * c) { // Create a new constraint, almost copy of "c", that becomes a lemma and could be modified later - constraint *lemma = constraint::make_ineq_lemma(m_max_constraint_id++, c->poly()); - lemma->predecessors() = c->predecessors(); // copy predecessors - lemma->add_predecessor(c); // and add c + constraint *lemma = new constraint(get_new_constraint_id(), c->poly(), true); + lemma->add_predecessor(c); TRACE("resolve_conflict_for_inequality", print_constraint(tout, *lemma);); lp_assert(lower_is_pos(lemma->poly())); bool done = false; unsigned j = m_trail.size() - 1; bool lemma_has_been_modified = false; - unsigned number_of_lemmas = m_lemmas.size(); + unsigned number_of_lemmas = m_lemmas_container.size(); while (!done) { if (cancel()) { return; @@ -2315,15 +2131,21 @@ public: } } - if (number_of_lemmas == m_lemmas.size()) { - if (lemma_has_been_modified) + if (number_of_lemmas == m_lemmas_container.size()) { + if (lemma_has_been_modified) { add_lemma_as_not_active(lemma); + } else { - delete lemma; + delete_constraint(lemma);; } } } + void delete_constraint(constraint * c) { + m_indexer_of_constraints.release_index(c->id()); + delete c; + } + void resolve_conflict(constraint * i) { lp_assert(!at_base_lvl()); TRACE("int_resolve_confl", tout << "inconstistent_constraint = "; @@ -2338,58 +2160,45 @@ public: void print_scope(std::ostream& out) const { for (const scope & s : m_scopes) { - out << "asserts_size = " << s.m_asserts_size; - out << ", lemmas_size = " << s.m_lemmas_size << "\n"; out << ", trail_size = " << s.m_trail_size << "\n"; } } - void remove_active_flag_from_constraints_in_active_set() { - for (auto c : m_active_set.cs()) { - c->remove_active_flag(); - } - } - - void set_active_flag_for_constraints_in_active_set() { - for (auto c : m_active_set.cs()) { - c->set_active_flag(); - } - } - public: + unsigned number_of_asserts() const { return m_asserts.size(); } + void push() { - m_scopes.push_back(scope(m_asserts.size(), m_lemmas.size(), m_trail.size())); + m_scopes.push_back(scope(m_trail.size())); TRACE("pp_cs", tout << "level = " << m_scopes.size() << ", trail size = " << m_trail.size();); } - - void pop_constraints(unsigned n_asserts, unsigned n_lemmas) { - if (n_asserts >= m_asserts.size()) - return; // only shrink the lemmas if asserts are shrunk - while (m_asserts.size() > n_asserts) { - constraint * i = m_asserts.back();; - for (auto & p: i->poly().m_coeffs) { - m_var_infos[p.var()].remove_depended_constraint(i); - } - m_active_set.remove_constraint(i); - delete i; - m_asserts.pop_back(); + void pop_last_assert() { + constraint * i = m_asserts.back();; + for (auto & p: i->poly().m_coeffs) { + m_var_infos[p.var()].remove_depended_constraint(i); } - - while (m_lemmas.size() > n_lemmas) { - constraint * i = m_lemmas.back();; + lp_assert(i->assert_origins().size() == 1); + for (constraint_index ci : i->assert_origins()) + m_lemmas_container.submit_assert_origin_for_delete(ci); + m_active_set.remove_constraint(i); + delete_constraint(i); + m_asserts.pop_back(); + } + + void pop_constraints() { + vector gone_lemmas = m_lemmas_container.remove_lemmas_depending_on_submitted_origins(); + + for (constraint * i : gone_lemmas) { for (auto & p: i->poly().m_coeffs) { m_var_infos[p.var()].remove_depended_constraint(i); } m_active_set.remove_constraint(i); - delete i; - m_lemmas.pop_back(); + delete_constraint(i); } } public: - void pop(unsigned k) { - m_number_of_pops++; + void pop_trail(unsigned k) { unsigned new_scope_size = m_scopes.size() - k; scope s = m_scopes[new_scope_size]; m_scopes.shrink(new_scope_size); @@ -2397,12 +2206,17 @@ public: literal& lit = m_trail[j]; if (lit.is_decided()) m_decision_level--; + else + m_active_set.add_constraint(lit.cnstr(), 2); // 2 is a priority; this constraint will not be processed urgentlty var_info & vi = m_var_infos[lit.var()]; if (vi.external_stack_level() != lit.prev_var_level()) vi.pop(1, lit.prev_var_level()); m_trail.pop_back(); } - pop_constraints(s.m_asserts_size, s.m_lemmas_size); + } + void pop(unsigned k) { + pop_trail(k); + pop_constraints(); lp_assert(var_infos_are_correct()); } public: @@ -2418,7 +2232,6 @@ public: m_number_of_variables_function(number_of_variables_function), m_var_value_function(var_value_function), m_settings(settings), - m_max_constraint_id(0), m_decision_level(0) {} @@ -2456,15 +2269,15 @@ public: } } - constraint* find_constraint_to_propagate(unsigned rand) { + constraint* find_constraint_to_propagate() { handle_conflicting_cores(); - return m_active_set.remove_random_constraint(rand); + return m_active_set.remove_constraint(); } // returns nullptr if there is no conflict, or a conflict constraint otherwise constraint* propagate_constraints_on_active_set() { constraint *c; - while ((c = find_constraint_to_propagate(m_settings.random_next())) != nullptr) { + while ((c = find_constraint_to_propagate()) != nullptr) { c = propagate_constraint(c); if (cancel()) { return nullptr; @@ -2520,7 +2333,7 @@ public: m_v.resize(j + 1); m_v[j] = b; TRACE("decide_var_on_bound", tout<< "j="<< j<<" ";print_var_info(tout, j);); - add_changed_var(j); + add_changed_var(j, !decide_on_lower, 1); m_decision_level++; literal nl = literal::make_decided_literal(j, !decide_on_lower, b, decision_context_index, var_level); push_literal_to_trail(nl); @@ -2586,7 +2399,8 @@ public: void simplify_ineq(polynomial & p) const { TRACE("simplify_ineq_int", tout << "p = " << pp_poly(*this, p) << "\n";); auto & ms = p.m_coeffs; - lp_assert(ms.size()); + if (ms.size() == 0) + return; mpq g; if (ms.size() == 1) { g = abs(ms[0].coeff()); @@ -2606,11 +2420,11 @@ public: } void add_lemma_common(constraint* lemma) { - m_lemmas.push_back(lemma); + m_lemmas_container.add_lemma(lemma); polynomial & p = lemma->poly(); simplify_ineq(p); for (const auto & m : p.coeffs()) { - m_var_infos[m.var()].add_dependent_constraint(lemma, is_pos(m.coeff())? bound_type::UPPER: bound_type::LOWER); + m_var_infos[m.var()].add_dependent_constraint(lemma, is_pos(m.coeff())); } } @@ -2619,16 +2433,17 @@ public: TRACE("add_lemma_int", trace_print_constraint(tout, lemma);); } - void add_lemma(constraint * lemma) { + void add_lemma(constraint * lemma, int priority) { add_lemma_common(lemma); - m_active_set.add_constraint(lemma); + m_active_set.add_constraint(lemma, priority); + lp_assert(constraint_indices_are_correct()); TRACE("add_lemma_int", trace_print_constraint(tout, lemma);); } unsigned add_ineq(const vector & lhs, const mpq& free_coeff, - svector origins) { + constraint_index origin) { lp_assert(lhs_is_int(lhs)); lp_assert(is_int(free_coeff)); for (auto & p : lhs) { @@ -2643,19 +2458,16 @@ public: } } - constraint * c = constraint::make_ineq_assert(m_max_constraint_id++, lhs, free_coeff,origins); + constraint * c = new constraint(get_new_constraint_id(), origin, polynomial(lhs, free_coeff), true); + + lp_assert(constraint_indices_are_correct()); + m_asserts.push_back(c); add_constraint_to_dependend_for_its_vars(c); - m_active_set.add_constraint(c); + m_active_set.add_constraint(c, 2); // 2 is a priority - TRACE("add_ineq_int", - tout << "explanation :"; - for (auto i: origins) { - m_print_constraint_function(i, tout); - tout << "\n"; - }); - - TRACE("add_ineq_int", tout << "m_asserts[" << m_asserts.size() - 1 << "] = "; + TRACE("add_ineq_int",tout << "explanation :"; m_print_constraint_function(origin, tout); tout << "\n"; + tout << "m_asserts[" << m_asserts.size() - 1 << "] = "; print_constraint(tout, *m_asserts.back()); tout << "\n";); return m_asserts.size() - 1; @@ -2664,13 +2476,16 @@ public: void add_constraint_to_dependend_for_its_vars(constraint * c) { for (auto & p : c->poly().coeffs()) { - m_var_infos[p.var()].add_dependent_constraint(c, is_pos(p.coeff())? bound_type::UPPER : bound_type::LOWER); + m_var_infos[p.var()].add_dependent_constraint(c, is_pos(p.coeff())); } } bool var_has_no_bounds(const var_info& vi) const { return !lower_bound_exists(vi) && !upper_bound_exists(vi); } + + unsigned number_of_constraints() const { return m_asserts.size() + m_lemmas_container.size(); } + }; inline polynomial operator*(const mpq & a, polynomial & p) { diff --git a/src/util/lp/eta_matrix_def.h b/src/util/lp/eta_matrix_def.h index d5f371335..5c7661e24 100644 --- a/src/util/lp/eta_matrix_def.h +++ b/src/util/lp/eta_matrix_def.h @@ -74,7 +74,7 @@ void eta_matrix::apply_from_right(vector & w) { t += w[it.first] * it.second; } w[m_column_index] = t; -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(vectors_are_equal(clone_w, w, get_number_of_rows())); // delete clone_w; #endif @@ -114,7 +114,7 @@ void eta_matrix::apply_from_right(indexed_vector & w) { } } -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(w.is_OK()); // lp_assert(vectors_are_equal(wcopy, w.m_data)); #endif @@ -144,7 +144,7 @@ void eta_matrix::conjugate_by_permutation(permutation_matrix & p) { for (auto & pair : m_column_vector.m_data) { pair.first = p.get_rev(pair.first); } -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(deb == *this); #endif } diff --git a/src/util/lp/hash_helper.h b/src/util/lp/hash_helper.h deleted file mode 100644 index ab5fa844b..000000000 --- a/src/util/lp/hash_helper.h +++ /dev/null @@ -1,54 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once -#include -#include -#include "util/numerics/mpq.h" -#ifdef __CLANG__ -#pragma clang diagnostic push -#pragma clang diagnostic ignored "-Wmismatched-tags" -#endif -namespace std { -template<> -struct hash { - inline size_t operator()(const lp::mpq & v) const { - return v.hash(); - } -}; -} - -template -inline void hash_combine(std::size_t & seed, const T & v) { - seed ^= std::hash()(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2); -} - -namespace std { -template struct hash> { - inline size_t operator()(const pair & v) const { - size_t seed = 0; - hash_combine(seed, v.first); - hash_combine(seed, v.second); - return seed; - } -}; -} -#ifdef __CLANG__ -#pragma clang diagnostic pop -#endif diff --git a/src/util/lp/indexed_vector.h b/src/util/lp/indexed_vector.h index aabc886f0..f17eb76ec 100644 --- a/src/util/lp/indexed_vector.h +++ b/src/util/lp/indexed_vector.h @@ -166,6 +166,55 @@ public: } } } + + struct ival { + unsigned m_var; + const T & m_coeff; + ival(unsigned var, const T & val) : m_var(var), m_coeff(val) { + } + unsigned var() const { return m_var;} + const T & coeff() const { return m_coeff; } + }; + + struct const_iterator { + // fields + const unsigned *m_i; + const indexed_vector& m_v; + + //typedefs + + + typedef const_iterator self_type; + typedef ival value_type; + typedef const ival reference; + // typedef const column_cell* pointer; + typedef int difference_type; + typedef std::forward_iterator_tag iterator_category; + + reference operator*() const { + return ival(*m_i, m_v[*m_i]); + } + self_type operator++() { self_type i = *this; m_i++; return i; } + self_type operator++(int) { m_i++; return *this; } + + const_iterator(const unsigned* it, const indexed_vector& v) : + m_i(it), + m_v(v) + {} + bool operator==(const self_type &other) const { + return m_i == other.m_i; + } + bool operator!=(const self_type &other) const { return !(*this == other); } + }; + + const_iterator begin() const { + return const_iterator(m_index.begin(), *this); + } + + const_iterator end() const { + return const_iterator(m_index.end(), *this); + } + #ifdef Z3DEBUG bool is_OK() const; diff --git a/src/util/lp/indexer_of_constraints.h b/src/util/lp/indexer_of_constraints.h new file mode 100644 index 000000000..30976d496 --- /dev/null +++ b/src/util/lp/indexer_of_constraints.h @@ -0,0 +1,45 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: + + +--*/ + +#pragma once +#include "util/lp/binary_heap_priority_queue.h" +namespace lp { + +class indexer_of_constraints { + binary_heap_priority_queue m_queue_of_released_indices; + unsigned m_max; +public: + indexer_of_constraints() :m_max(0) {} + unsigned get_new_index() { + unsigned ret; + if (m_queue_of_released_indices.is_empty()) { + ret = m_max++; + } + else { + ret = m_queue_of_released_indices.dequeue(); + } + return ret; + }; + void release_index(unsigned i) { + m_queue_of_released_indices.enqueue(i, i); + }; + unsigned max() const { return m_max; } +}; +} diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index acb468058..1b948e2d6 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -31,17 +31,17 @@ void int_solver::trace_inf_rows() const { } } - num = 0; - for (unsigned i = 0; i < m_lar_solver->A_r().row_count(); i++) { - unsigned j = m_lar_solver->m_mpq_lar_core_solver.m_r_basis[i]; - if (column_is_int_inf(j)) { - num++; - iterator_on_row it(m_lar_solver->A_r().m_rows[i]); - m_lar_solver->print_linear_iterator(&it, tout); - tout << "\n"; - } - } - tout << "num of int infeasible: " << num << "\n"; + num = 0; + for (unsigned i = 0; i < m_lar_solver->A_r().row_count(); i++) { + unsigned j = m_lar_solver->m_mpq_lar_core_solver.m_r_basis[i]; + if (column_is_int_inf(j)) { + num++; + m_lar_solver->print_row(m_lar_solver->A_r().m_rows[i], tout); + tout << "\n"; + } + } + tout << "num of int infeasible: " << num << "\n"; + ); } int_set& int_solver::inf_int_set() { @@ -106,32 +106,29 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() { } -bool int_solver::is_gomory_cut_target(linear_combination_iterator &iter) { - unsigned j; - lp_assert(iter.is_reset()); +bool int_solver::is_gomory_cut_target(const row_strip& row) { // All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed). - while (iter.next(j)) { + unsigned j; + for (auto p : row) { + j = p.var(); if (is_base(j)) continue; if (!is_zero(get_value(j).y)) { TRACE("gomory_cut", tout << "row is not gomory cut target:\n"; display_column(tout, j); tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";); - iter.reset(); return false; } } - iter.reset(); return true; } -void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl, unsigned gomory_cut_inf_column) { +void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl, const mpq& f_0, const mpq& one_minus_f_0) { TRACE("gomory_cut_detail_real", tout << "real\n";); - mpq f_0 = fractional_part(get_value(gomory_cut_inf_column)); mpq new_a; if (at_low(x_j)) { if (a.is_pos()) { - new_a = a / (1 - f_0); + new_a = a / one_minus_f_0; } else { new_a = a / f_0; @@ -148,7 +145,7 @@ void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, l new_a.neg(); // the upper terms are inverted. } else { - new_a = a / (mpq(1) - f_0); + new_a = a / one_minus_f_0; } k.addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a; expl.push_justification(column_upper_bound_constraint(x_j), new_a); @@ -166,11 +163,9 @@ constraint_index int_solver::column_lower_bound_constraint(unsigned j) const { } -void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & t, explanation& expl, mpq & lcm_den, unsigned inf_column) { +void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & t, explanation& expl, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) { lp_assert(is_int(x_j)); lp_assert(!a.is_int()); - mpq f_0 = fractional_part(get_value(inf_column)); - lp_assert(f_0 > zero_of_type() && f_0 < one_of_type()); mpq f_j = fractional_part(a); TRACE("gomory_cut_detail", tout << a << " x_j" << x_j << " k = " << k << "\n"; @@ -182,7 +177,6 @@ void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, la lp_assert (!f_j.is_zero()); mpq new_a; if (at_low(x_j)) { - auto one_minus_f_0 = 1 - f_0; if (f_j <= one_minus_f_0) { new_a = f_j / one_minus_f_0; } @@ -198,7 +192,7 @@ void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, la new_a = f_j / f_0; } else { - new_a = (mpq(1) - f_j) / (1 - f_0); + new_a = (mpq(1) - f_j) / one_minus_f_0; } new_a.neg(); // the upper terms are inverted k.addmul(new_a, upper_bound(x_j).x); @@ -319,18 +313,15 @@ void int_solver::adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k -lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsigned inf_col, linear_combination_iterator& iter) { +lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsigned inf_col, const row_strip & row) { lp_assert(column_is_int_inf(inf_col)); TRACE("gomory_cut", - tout << "applying cut at:\n"; m_lar_solver->print_linear_iterator_indices_only(&iter, tout); tout << std::endl; - iter.reset(); - unsigned j; - while(iter.next(j)) { - m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(j, tout); + tout << "applying cut at:\n"; m_lar_solver->print_row(row, tout); tout << std::endl; + for (auto p : row) { + m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout); } - iter.reset(); tout << "inf_col = " << inf_col << std::endl; ); @@ -340,18 +331,21 @@ lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsi unsigned x_j; mpq a; bool some_int_columns = false; - lp_assert(iter.is_reset()); - while (iter.next(a, x_j)) { + mpq f_0 = int_solver::fractional_part(get_value(inf_col)); + mpq one_min_f_0 = 1 - f_0; + for (auto p : row) { + x_j = p.var(); if (x_j == inf_col) continue; // make the format compatible with the format used in: Integrating Simplex with DPLL(T) + a = p.coeff(); a.neg(); if (is_real(x_j)) - real_case_in_gomory_cut(a, x_j, k, t, expl, inf_col); + real_case_in_gomory_cut(a, x_j, k, t, expl, f_0, one_min_f_0); else { if (a.is_int()) continue; // f_j will be zero and no monomial will be added some_int_columns = true; - int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den, inf_col); + int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den, f_0, one_min_f_0); } } @@ -362,6 +356,7 @@ lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsi lp_assert(current_solution_is_inf_on_cut(t, k)); m_lar_solver->subs_term_columns(t); + TRACE("gomory_cut", tout<<"precut:"; m_lar_solver->print_term(t, tout); tout << ">= " << k << std::endl;); return lia_move::cut; } @@ -371,28 +366,30 @@ void int_solver::init_check_data() { m_old_values_data.resize(n); } -int int_solver::find_free_var_in_gomory_row(linear_combination_iterator& iter) { +int int_solver::find_free_var_in_gomory_row(const row_strip& row) { unsigned j; - while(iter.next(j)) { + for (auto p : row) { + j = p.var(); if (!is_base(j) && is_free(j)) return static_cast(j); } - iter.reset(); return -1; } -lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j) { +lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, bool & upper) { lia_move ret; - linear_combination_iterator* iter = m_lar_solver->get_iterator_on_row(row_of_basic_column(j)); - int free_j = find_free_var_in_gomory_row(*iter); + + const row_strip& row = m_lar_solver->get_row(row_of_basic_column(j)); + int free_j = find_free_var_in_gomory_row(row); if (free_j != -1) { - ret = create_branch_on_column(j, t, k, true); - } else if (!is_gomory_cut_target(*iter)) { - ret = create_branch_on_column(j, t, k, false); + ret = create_branch_on_column(j, t, k, true, upper); + } else if (!is_gomory_cut_target(row)) { + bool upper; + ret = create_branch_on_column(j, t, k, false, upper); } else { - ret = mk_gomory_cut(t, k, ex, j, *iter); + upper = false; + ret = mk_gomory_cut(t, k, ex, j, row); } - delete iter; return ret; } @@ -478,7 +475,14 @@ void int_solver::copy_values_from_cut_solver() { } } -lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) { +void int_solver::catch_up_in_adding_constraints_to_cut_solver() { + lp_assert(m_cut_solver.number_of_asserts() <= m_lar_solver->constraints().size()); + for (unsigned j = m_cut_solver.number_of_asserts(); j < m_lar_solver->constraints().size(); j++) { + add_constraint_to_cut_solver(j, m_lar_solver->constraints()[j]); + } +} + +lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) { init_check_data(); lp_assert(inf_int_set_is_correct()); // it is a reimplementation of @@ -486,9 +490,16 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) { // from theory_arith_int.h with the addition of cut_solver if (!has_inf_int()) return lia_move::ok; - if (settings().m_run_gcd_test) - if (!gcd_test(ex)) + if (settings().m_run_gcd_test) { + settings().st().m_gcd_calls++; + if (!gcd_test(ex)) { + TRACE("gcd_test", tout << "conflict";); + settings().st().m_gcd_conflicts++; return lia_move::conflict; + } + } else { + TRACE("gcd_test", tout << "no test";); + } pivoted_rows_tracking_control pc(m_lar_solver); /* if (m_params.m_arith_euclidean_solver) apply_euclidean_solver(); */ //m_lar_solver->pivot_fixed_vars_from_basis(); @@ -498,20 +509,20 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) { if ((++m_branch_cut_counter) % settings().m_int_branch_cut_solver == 0) { TRACE("check_main_int", tout<<"cut_solver";); + catch_up_in_adding_constraints_to_cut_solver(); auto check_res = m_cut_solver.check(); settings().st().m_cut_solver_calls++; switch (check_res) { - case lbool::l_false: + case cut_solver::lbool::l_false: copy_explanations_from_cut_solver(ex); settings().st().m_cut_solver_false++; return lia_move::conflict; - case lbool::l_true: + case cut_solver::lbool::l_true: settings().st().m_cut_solver_true++; copy_values_from_cut_solver(); return lia_move::ok; - case lbool::l_undef: + case cut_solver::lbool::l_undef: settings().st().m_cut_solver_undef++; - settings().m_int_branch_cut_solver *= (settings().m_int_branch_cut_solver); // take a square break; default: return lia_move::give_up; @@ -530,10 +541,11 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) { int j = find_inf_int_base_column(); if (j == -1) return lia_move::ok; TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";); - return proceed_with_gomory_cut(t, k, ex, j); + + return proceed_with_gomory_cut(t, k, ex, j, upper); } TRACE("check_main_int", tout << "branch"; ); - return create_branch_on_column(find_inf_int_base_column(), t, k, false); + return create_branch_on_column(find_inf_int_base_column(), t, k, false, upper); } bool int_solver::move_non_basic_column_to_bounds(unsigned j) { @@ -668,26 +680,24 @@ void int_solver::patch_int_infeasible_nbasic_columns() { lp_assert(is_feasible() && inf_int_set_is_correct()); } -mpq get_denominators_lcm(iterator_on_row &it) { +mpq get_denominators_lcm(const row_strip & row) { mpq r(1); - mpq a; - unsigned j; - while (it.next(a, j)) { - r = lcm(r, denominator(a)); + for (auto c : row) { + r = lcm(r, denominator(c.coeff())); } return r; } bool int_solver::gcd_test_for_row(static_matrix> & A, unsigned i, explanation & ex) { - iterator_on_row it(A.m_rows[i]); - mpq lcm_den = get_denominators_lcm(it); + mpq lcm_den = get_denominators_lcm(A.m_rows[i]); mpq consts(0); mpq gcds(0); mpq least_coeff(0); bool least_coeff_is_bounded = false; - mpq a; unsigned j; - while (it.next(a, j)) { + for (auto &c : A.m_rows[i]) { + j = c.var(); + const mpq& a = c.coeff(); if (m_lar_solver->column_is_fixed(j)) { mpq aux = lcm_den * a; consts += aux * m_lar_solver->column_lower_bound(j).x; @@ -725,8 +735,11 @@ bool int_solver::gcd_test_for_row(static_matrix> & A, uns return true; } - if (!(consts / gcds).is_int()) - fill_explanation_from_fixed_columns(it, ex); + if (!(consts / gcds).is_int()) { + TRACE("gcd_test", tout << "row failed the GCD test:\n"; display_row_info(tout, i);); + fill_explanation_from_fixed_columns(A.m_rows[i], ex); + return false; + } if (least_coeff.is_one() && !least_coeff_is_bounded) { SASSERT(gcds.is_one()); @@ -734,7 +747,7 @@ bool int_solver::gcd_test_for_row(static_matrix> & A, uns } if (least_coeff_is_bounded) { - return ext_gcd_test(it, least_coeff, lcm_den, consts, ex); + return ext_gcd_test(A.m_rows[i], least_coeff, lcm_den, consts, ex); } return true; } @@ -745,13 +758,11 @@ void int_solver::add_to_explanation_from_fixed_or_boxed_column(unsigned j, expla ex.m_explanation.push_back(std::make_pair(mpq(1), lc)); ex.m_explanation.push_back(std::make_pair(mpq(1), uc)); } -void int_solver::fill_explanation_from_fixed_columns(iterator_on_row & it, explanation & ex) { - it.reset(); - unsigned j; - while (it.next(j)) { - if (!m_lar_solver->column_is_fixed(j)) +void int_solver::fill_explanation_from_fixed_columns(const row_strip & row, explanation & ex) { + for (const auto & c : row) { + if (!m_lar_solver->column_is_fixed(c.var())) continue; - add_to_explanation_from_fixed_or_boxed_column(j, ex); + add_to_explanation_from_fixed_or_boxed_column(c.var(), ex); } } @@ -759,14 +770,13 @@ bool int_solver::gcd_test(explanation & ex) { auto & A = m_lar_solver->A_r(); // getting the matrix for (unsigned i = 0; i < A.row_count(); i++) if (!gcd_test_for_row(A, i, ex)) { - std::cout << "false from gcd_test\n" ; return false; } return true; } -bool int_solver::ext_gcd_test(iterator_on_row & it, +bool int_solver::ext_gcd_test(const row_strip & row, mpq const & least_coeff, mpq const & lcm_den, mpq const & consts, explanation& ex) { @@ -774,10 +784,11 @@ bool int_solver::ext_gcd_test(iterator_on_row & it, mpq l(consts); mpq u(consts); - it.reset(); mpq a; unsigned j; - while (it.next(a, j)) { + for (const auto & c : row) { + j = c.var(); + const mpq & a = c.coeff(); if (m_lar_solver->column_is_fixed(j)) continue; SASSERT(!m_lar_solver->column_is_real(j)); @@ -817,20 +828,20 @@ bool int_solver::ext_gcd_test(iterator_on_row & it, mpq u1 = floor(u/gcds); if (u1 < l1) { - fill_explanation_from_fixed_columns(it, ex); + fill_explanation_from_fixed_columns(row, ex); return false; } return true; } - +/* linear_combination_iterator * int_solver::get_column_iterator(unsigned j) { if (m_lar_solver->use_tableau()) return new iterator_on_column(m_lar_solver->A_r().m_columns[j], m_lar_solver->A_r()); return new iterator_on_indexed_vector(m_lar_solver->get_column_in_lu_mode(j)); } - +*/ int_solver::int_solver(lar_solver* lar_slv) : m_lar_solver(lar_slv), @@ -893,8 +904,7 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq return false; impq const & xj = get_value(j); - linear_combination_iterator *it = get_column_iterator(j); - + inf_l = true; inf_u = true; l = u = zero_of_type(); @@ -909,7 +919,12 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq mpq a; // the coefficient in the column unsigned row_index; - while (it->next(a, row_index)) { + lp_assert(settings().use_tableau()); + const auto & A = m_lar_solver->A_r(); + for (auto c : A.column(j)) { + row_index = c.var(); + const mpq & a = c.coeff(); + unsigned i = lcs.m_r_basis[row_index]; impq const & xi = get_value(i); if (is_int(i) && is_int(j) && !a.is_int()) @@ -930,7 +945,6 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq if (!inf_l && !inf_u && l == u) break;; } - delete it; TRACE("freedom_interval", tout << "freedom variable for:\n"; tout << m_lar_solver->get_column_name(j); @@ -1062,21 +1076,16 @@ lp_settings& int_solver::settings() { void int_solver::display_row_info(std::ostream & out, unsigned row_index) const { auto & rslv = m_lar_solver->m_mpq_lar_core_solver.m_r_solver; - auto it = m_lar_solver->get_iterator_on_row(row_index); - mpq a; - unsigned j; - while (it->next(a, j)) { - if (numeric_traits::is_pos(a)) + for (auto &c: rslv.m_A.m_rows[row_index]) { + if (numeric_traits::is_pos(c.coeff())) out << "+"; - out << a << rslv.column_name(j) << " "; + out << c.coeff() << rslv.column_name(c.var()) << " "; } - it->reset(); - while(it->next(j)) { - rslv.print_column_bound_info(j, out); + for (auto& c: rslv.m_A.m_rows[row_index]) { + rslv.print_column_bound_info(c.var(), out); } rslv.print_column_bound_info(rslv.m_basis[row_index], out); - delete it; } unsigned int_solver::random() { @@ -1171,11 +1180,18 @@ const impq& int_solver::lower_bound(unsigned j) const { return m_lar_solver->column_lower_bound(j); } -lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column) const { +lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column, bool & upper) { lp_assert(t.is_empty()); lp_assert(j != -1); t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j)); - k = free_column? mpq(0) : floor(get_value(j)); + if (free_column) { + upper = true; + k = mpq(0); + } else { + upper = left_branch_is_more_narrow_than_right(j); + k = upper? floor(get_value(j)) : ceil(get_value(j)); + } + TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n"; display_column(tout, j); tout << "k = " << k << std::endl; @@ -1184,6 +1200,25 @@ lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool fr } +bool int_solver::left_branch_is_more_narrow_than_right(unsigned j) { + return settings().random_next() % 2; + switch (m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_column_types[j] ) { + case column_type::fixed: + return false; + case column_type::boxed: + { + auto k = floor(get_value(j)); + return k - lower_bound(j).x < upper_bound(j).x - (k + mpq(1)); + } + case column_type::lower_bound: + return true; + case column_type::upper_bound: + return false; + default: + return false; + } +} + const impq& int_solver::upper_bound(unsigned j) const { return m_lar_solver->column_upper_bound(j); } @@ -1205,19 +1240,14 @@ void int_solver::add_constraint_to_cut_solver(unsigned ci, const lar_base_constr vector coeffs; mpq rs; get_int_coeffs_from_constraint(c, coeffs, rs); - svector explanation; - explanation.push_back(ci); - m_cut_solver.add_ineq(coeffs, -rs, explanation); -} - -void int_solver::notify_on_last_added_constraint() { - unsigned ci = m_lar_solver->constraints().size() - 1; - const lar_base_constraint* c = m_lar_solver->constraints()[ci]; - add_constraint_to_cut_solver(ci, c); + m_cut_solver.add_ineq(coeffs, -rs, ci); } void int_solver::pop(unsigned k) { - m_cut_solver.pop(k); + m_cut_solver.pop_trail(k); + while (m_cut_solver.number_of_asserts() > m_lar_solver->constraints().size()) + m_cut_solver.pop_last_assert(); + m_cut_solver.pop_constraints(); } void int_solver::push() { diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index 2c41220c9..badb55beb 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -1,11 +1,25 @@ -/* - Copyright (c) 2017 Microsoft Corporation - Author: Lev Nachmanson -*/ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: + + +--*/ #pragma once #include "util/lp/lp_settings.h" #include "util/lp/static_matrix.h" -#include "util/lp/iterator_on_row.h" #include "util/lp/int_set.h" #include "util/lp/lar_term.h" #include "util/lp/cut_solver.h" @@ -49,7 +63,7 @@ public: const int_set& inf_int_set() const; // main function to check that the solution provided by lar_solver is valid for integral values, // or provide a way of how it can be adjusted. - lia_move check(lar_term& t, mpq& k, explanation& ex); + lia_move check(lar_term& t, mpq& k, explanation& ex, bool & upper); bool move_non_basic_column_to_bounds(unsigned j); lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex); private: @@ -76,17 +90,16 @@ private: // creates a fresh inequality. bool branch(const lp_constraint & new_inequality); - bool ext_gcd_test(iterator_on_row & it, + bool ext_gcd_test(const row_strip& row, mpq const & least_coeff, mpq const & lcm_den, mpq const & consts, explanation & ex); - void fill_explanation_from_fixed_columns(iterator_on_row & it, explanation &); + void fill_explanation_from_fixed_columns(const row_strip & row, explanation &); void add_to_explanation_from_fixed_or_boxed_column(unsigned j, explanation &); void patch_int_infeasible_non_basic_column(unsigned j); void patch_int_infeasible_nbasic_columns(); bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m); - linear_combination_iterator * get_column_iterator(unsigned j); const impq & lower_bound(unsigned j) const; const impq & upper_bound(unsigned j) const; bool is_int(unsigned j) const; @@ -112,14 +125,13 @@ private: lp_settings& settings(); bool move_non_basic_columns_to_bounds(); void branch_infeasible_int_var(unsigned); - lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex, unsigned inf_col, linear_combination_iterator& iter); + lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex, unsigned inf_col, const row_strip& row); lia_move report_conflict_from_gomory_cut(mpq & k); void adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq& lcm_den); void init_check_data(); - bool constrain_free_vars(linear_combination_iterator * r); - lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j); - int find_free_var_in_gomory_row(linear_combination_iterator& iter); - bool is_gomory_cut_target(linear_combination_iterator &iter); + lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, bool & upper); + int find_free_var_in_gomory_row(const row_strip& ); + bool is_gomory_cut_target(const row_strip&); bool at_bound(unsigned j) const; bool at_low(unsigned j) const; bool at_upper(unsigned j) const; @@ -130,13 +142,15 @@ private: return is_zero(n.y); } +public: inline static mpq fractional_part(const impq & n) { lp_assert(is_rational(n)); return n.x - floor(n.x); } - void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex, unsigned inf_column); - void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den, unsigned inf_column); +private: + void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex, const mpq& f_0, const mpq& one_minus_f_0); + void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0); constraint_index column_upper_bound_constraint(unsigned j) const; constraint_index column_lower_bound_constraint(unsigned j) const; void display_row_info(std::ostream & out, unsigned row_index) const; @@ -147,7 +161,8 @@ public: private: unsigned random(); bool has_inf_int() const; - lia_move create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column) const; + lia_move create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column, bool & upper); + void catch_up_in_adding_constraints_to_cut_solver(); public: void display_inf_or_int_inf_columns(std::ostream & out) const; template @@ -155,11 +170,11 @@ public: template void get_int_coeffs_from_constraint(const lar_base_constraint* c, vector& coeff, T & rs); bool is_term(unsigned j) const; - void notify_on_last_added_constraint(); void add_constraint_to_cut_solver(unsigned,const lar_base_constraint*); void copy_explanations_from_cut_solver(explanation &); void pop(unsigned); void push(); void copy_values_from_cut_solver(); + bool left_branch_is_more_narrow_than_right(unsigned); }; } diff --git a/src/util/lp/iterator_on_column.h b/src/util/lp/iterator_on_column.h deleted file mode 100644 index 0d3b99ef5..000000000 --- a/src/util/lp/iterator_on_column.h +++ /dev/null @@ -1,67 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once -#include "util/lp/linear_combination_iterator.h" -#include "util/lp/static_matrix.h" -#include "util/lp/lar_term.h" -namespace lp { -template -struct iterator_on_column:linear_combination_iterator { - const vector& m_column; // the offset in term coeffs - const static_matrix & m_A; - int m_i; // the initial offset in the column - unsigned size() const override { return m_column.size(); } - iterator_on_column(const vector& column, const static_matrix & A) // the offset in term coeffs - : - m_column(column), - m_A(A), - m_i(-1) {} - - bool next(mpq & a, unsigned & i) override { - if (++m_i >= static_cast(m_column.size())) - return false; - - const column_cell& c = m_column[m_i]; - a = m_A.get_val(c); - i = c.m_i; - return true; - } - - bool next(unsigned & i) override { - if (++m_i >= static_cast(m_column.size())) - return false; - - const column_cell& c = m_column[m_i]; - i = c.m_i; - return true; - } - - void reset() override { - m_i = -1; - } - - bool is_reset() const { return m_i == -1;} - - linear_combination_iterator * clone() { - iterator_on_column * r = new iterator_on_column(m_column, m_A); - return r; - } -}; -} diff --git a/src/util/lp/iterator_on_indexed_vector.h b/src/util/lp/iterator_on_indexed_vector.h deleted file mode 100644 index b0ef87506..000000000 --- a/src/util/lp/iterator_on_indexed_vector.h +++ /dev/null @@ -1,56 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once -#include "util/lp/linear_combination_iterator.h" -namespace lp { -template -struct iterator_on_indexed_vector:linear_combination_iterator { - const indexed_vector & m_v; - unsigned m_offset; - iterator_on_indexed_vector(const indexed_vector & v) : - m_v(v), - m_offset(0) - {} - unsigned size() const override { return m_v.m_index.size(); } - bool next(T & a, unsigned & i) override { - if (m_offset >= m_v.m_index.size()) - return false; - i = m_v.m_index[m_offset++]; - a = m_v.m_data[i]; - return true; - } - - bool next(unsigned & i) override { - if (m_offset >= m_v.m_index.size()) - return false; - i = m_v.m_index[m_offset++]; - return true; - } - void reset() override { - m_offset = 0; - } - - bool is_reset() const { return m_offset == 0;} - - linear_combination_iterator* clone() { - return new iterator_on_indexed_vector(m_v); - } -}; -} diff --git a/src/util/lp/iterator_on_pivot_row.h b/src/util/lp/iterator_on_pivot_row.h deleted file mode 100644 index e20533c55..000000000 --- a/src/util/lp/iterator_on_pivot_row.h +++ /dev/null @@ -1,62 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once -#include "util/lp/iterator_on_indexed_vector.h" -namespace lp { -template -struct iterator_on_pivot_row:linear_combination_iterator { - bool m_basis_returned; - const indexed_vector & m_v; - unsigned m_basis_j; - iterator_on_indexed_vector m_it; - unsigned size() const override { return m_it.size(); } - iterator_on_pivot_row(const indexed_vector & v, unsigned basis_j) : - m_basis_returned(false), - m_v(v), m_basis_j(basis_j), m_it(v) {} - bool next(T & a, unsigned & i) override { - if (m_basis_returned == false) { - m_basis_returned = true; - a = one_of_type(); - i = m_basis_j; - return true; - } - return m_it.next(a, i); - } - bool next(unsigned & i) override { - if (m_basis_returned == false) { - m_basis_returned = true; - i = m_basis_j; - return true; - } - return m_it.next(i); - } - void reset() override { - m_basis_returned = false; - m_it.reset(); - } - - bool is_reset() const { return m_basis_returned == false && m_it.is_reset();} - - linear_combination_iterator * clone() { - iterator_on_pivot_row * r = new iterator_on_pivot_row(m_v, m_basis_j); - return r; - } -}; -} diff --git a/src/util/lp/iterator_on_row.h b/src/util/lp/iterator_on_row.h deleted file mode 100644 index 36b1d22ad..000000000 --- a/src/util/lp/iterator_on_row.h +++ /dev/null @@ -1,55 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once -#include "util/lp/linear_combination_iterator.h" -namespace lp { -template -struct iterator_on_row:linear_combination_iterator { - const vector> & m_row; - unsigned m_i; // offset - iterator_on_row(const vector> & row) : m_row(row), m_i(0) - {} - unsigned size() const override { return m_row.size(); } - bool next(T & a, unsigned & i) override { - if (m_i == m_row.size()) - return false; - auto &c = m_row[m_i++]; - i = c.m_j; - a = c.get_val(); - return true; - } - bool next(unsigned & i) override { - if (m_i == m_row.size()) - return false; - auto &c = m_row[m_i++]; - i = c.m_j; - return true; - } - void reset() override { - m_i = 0; - } - - bool is_reset() const { return m_i == 0;} - - linear_combination_iterator* clone() { - return new iterator_on_row(m_row); - } -}; -} diff --git a/src/util/lp/iterator_on_term_with_basis_var.h b/src/util/lp/iterator_on_term_with_basis_var.h deleted file mode 100644 index cf36046cb..000000000 --- a/src/util/lp/iterator_on_term_with_basis_var.h +++ /dev/null @@ -1,75 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once -#include "util/lp/linear_combination_iterator.h" -#include "util/lp/numeric_pair.h" -#include "util/lp/lar_term.h" -namespace lp { -struct iterator_on_term_with_basis_var:linear_combination_iterator { - const lar_term & m_term; - std::unordered_map::const_iterator m_i; // the offset in term coeffs - bool m_term_j_returned; - unsigned m_term_j; - unsigned size() const override {return static_cast(m_term.m_coeffs.size() + 1);} - iterator_on_term_with_basis_var(const lar_term & t, unsigned term_j) : - m_term(t), - m_i(t.m_coeffs.begin()), - m_term_j_returned(false), - m_term_j(term_j) {} - - bool next(mpq & a, unsigned & i) override { - if (m_term_j_returned == false) { - m_term_j_returned = true; - a = - one_of_type(); - i = m_term_j; - return true; - } - if (m_i == m_term.m_coeffs.end()) - return false; - i = m_i->first; - a = m_i->second; - m_i++; - return true; - } - bool next(unsigned & i) override { - if (m_term_j_returned == false) { - m_term_j_returned = true; - i = m_term_j; - return true; - } - if (m_i == m_term.m_coeffs.end()) - return false; - i = m_i->first; - m_i++; - return true; - } - - bool is_reset() const { return m_term_j_returned == false && m_i == m_term.m_coeffs.begin();} - - void reset() { - m_term_j_returned = false; - m_i = m_term.m_coeffs.begin(); - } - linear_combination_iterator * clone() override { - iterator_on_term_with_basis_var * r = new iterator_on_term_with_basis_var(m_term, m_term_j); - return r; - } -}; -} diff --git a/src/util/lp/lar_core_solver.h b/src/util/lp/lar_core_solver.h index c5b281f51..70277c848 100644 --- a/src/util/lp/lar_core_solver.h +++ b/src/util/lp/lar_core_solver.h @@ -30,8 +30,6 @@ Revision History: #include "util/lp/lp_primal_core_solver.h" #include "util/lp/stacked_vector.h" #include "util/lp/lar_solution_signature.h" -#include "util/lp/iterator_on_column.h" -#include "util/lp/iterator_on_indexed_vector.h" #include "util/lp/stacked_value.h" namespace lp { @@ -801,15 +799,6 @@ public: m_r_solver.init_column_row_non_zeroes(); } - linear_combination_iterator * get_column_iterator(unsigned j) { - if (settings().use_tableau()) { - return new iterator_on_column>(m_r_solver.m_A.m_columns[j], m_r_solver.m_A); - } else { - m_r_solver.solve_Bd(j); - return new iterator_on_indexed_vector(m_r_solver.m_ed); - } - } - bool column_is_fixed(unsigned j) const { return m_column_types()[j] == column_type::fixed || ( m_column_types()[j] == column_type::boxed && diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index fb975186c..f0c9cf09a 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -151,9 +151,10 @@ void lar_solver::analyze_new_bounds_on_row( unsigned row_index, bound_propagator & bp) { lp_assert(!use_tableau()); - iterator_on_pivot_row it(m_mpq_lar_core_solver.get_pivot_row(), m_mpq_lar_core_solver.m_r_basis[row_index]); - - bound_analyzer_on_row ra_pos(it, + unsigned j = m_mpq_lar_core_solver.m_r_basis[row_index]; // basis column for the row + bound_analyzer_on_row> + ra_pos(m_mpq_lar_core_solver.get_pivot_row(), + j, zero_of_type>(), row_index, bp @@ -163,14 +164,14 @@ void lar_solver::analyze_new_bounds_on_row( void lar_solver::analyze_new_bounds_on_row_tableau( unsigned row_index, - bound_propagator & bp - ) { + bound_propagator & bp ) { if (A_r().m_rows[row_index].size() > settings().max_row_length_for_bound_propagation) return; - iterator_on_row it(A_r().m_rows[row_index]); + lp_assert(use_tableau()); - bound_analyzer_on_row::analyze_row(it, + bound_analyzer_on_row>::analyze_row(A_r().m_rows[row_index], + static_cast(-1), zero_of_type>(), row_index, bp @@ -200,13 +201,6 @@ void lar_solver::calculate_implied_bounds_for_row(unsigned i, bound_propagator & } } - -linear_combination_iterator * lar_solver::create_new_iter_from_term(unsigned term_index) const { - lp_assert(false); // not implemented - return nullptr; - // new linear_combination_iterator_on_vector(m_terms[adjust_term_index(term_index)]->coeffs_as_vector()); -} - unsigned lar_solver::adjust_column_index_to_term_index(unsigned j) const { unsigned ext_var_or_term = m_columns_to_ext_vars_or_term_indices[j]; return ext_var_or_term < m_terms_start_index ? j : ext_var_or_term; @@ -407,6 +401,9 @@ void lar_solver::pop(unsigned k) { m_constraints.resize(m_constraint_count); m_term_count.pop(k); for (unsigned i = m_term_count; i < m_terms.size(); i++) { +#if Z3DEBUG_CHECK_UNIQUE_TERMS + m_set_of_terms.erase(m_terms[i]); +#endif delete m_terms[i]; } m_terms.resize(m_term_count); @@ -595,7 +592,8 @@ void lar_solver::substitute_terms_in_linear_expression(const vector>& explanation) const { -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG lconstraint_kind kind; lp_assert(the_relations_are_of_same_type(explanation, kind)); lp_assert(the_left_sides_sum_to_zero(explanation)); @@ -1048,7 +1046,7 @@ bool lar_solver::explanation_is_correct(const vector>& } bool lar_solver::inf_explanation_is_correct() const { -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG vector> explanation; get_infeasibility_explanation(explanation); return explanation_is_correct(explanation); @@ -1190,6 +1188,7 @@ void lar_solver::print_constraint(constraint_index ci, std::ostream & out) const } void lar_solver::print_constraints(std::ostream& out) const { + out << "number of constraints = " << m_constraints.size() << std::endl; for (auto c : m_constraints) { print_constraint(c, out); } @@ -1214,7 +1213,26 @@ void lar_solver::print_term(lar_term const& term, std::ostream & out) const { if (!numeric_traits::is_zero(term.m_v)) { out << term.m_v << " + "; } - print_linear_combination_of_column_indices(term.coeffs_as_vector(), out); + bool first = true; + for (const auto p : term) { + mpq val = p.coeff(); + if (first) { + first = false; + } else { + if (is_pos(val)) { + out << " + "; + } else { + out << " - "; + val = -val; + } + } + if (val == -numeric_traits::one()) + out << " - "; + else if (val != numeric_traits::one()) + out << T_to_string(val); + out << this->get_column_name(p.var()); + } + } void lar_solver::print_term_as_indices(lar_term const& term, std::ostream & out) const { @@ -1584,17 +1602,59 @@ void lar_solver::add_new_var_to_core_fields_for_mpq(bool register_in_basis) { var_index lar_solver::add_term_undecided(const vector> & coeffs, const mpq &m_v) { - m_terms.push_back(new lar_term(coeffs, m_v)); + push_and_register_term(new lar_term(coeffs, m_v)); return m_terms_start_index + m_terms.size() - 1; } +#if Z3DEBUG_CHECK_UNIQUE_TERMS +bool lar_solver::term_coeffs_are_ok(const vector> & coeffs, const mpq& v) { + if (coeffs.empty()) { + return is_zero(v); + } + + for (const auto & p : coeffs) { + if (column_is_real(p.second)) + return true; + } + + mpq g; + bool g_is_set = false; + for (const auto & p : coeffs) { + if (!p.first.is_int()) { + std::cout << "p.first = " << p.first << " is not an int\n"; + return false; + } + if (!g_is_set) { + g_is_set = true; + g = p.first; + } else { + g = gcd(g, p.first); + } + } + if (g == one_of_type()) + return true; + + std::cout << "term is not ok: g = " << g << std::endl; + this->print_linear_combination_of_column_indices_only(coeffs, std::cout); + std::cout << " rs = " << v << std::endl; + return false; +} +#endif +void lar_solver::push_and_register_term(lar_term* t) { +#if Z3DEBUG_CHECK_UNIQUE_TERMS + lp_assert(m_set_of_terms.find(t) == m_set_of_terms.end()); + m_set_of_terms.insert(t); +#endif + m_terms.push_back(t); +} + // terms var_index lar_solver::add_term(const vector> & 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)); + push_and_register_term(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()) { @@ -1607,6 +1667,7 @@ var_index lar_solver::add_term(const vector> & coeffs, } void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) { + TRACE("dump_terms", print_term(*term, tout); tout << std::endl;); register_new_ext_var_index(term_ext_index, term_is_int(term)); // j will be a new variable unsigned j = A_r().column_count(); @@ -1614,8 +1675,8 @@ void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned m_columns_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, + A_r().fill_last_row_with_pivoting(*term, + j, 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()); } @@ -1653,9 +1714,6 @@ constraint_index lar_solver::add_var_bound(var_index j, lconstraint_kind kind, c lp_assert(bound_is_integer_for_integer_column(j, right_side)); auto vc = new lar_var_constraint(j, kind, right_side); m_constraints.push_back(vc); - if (has_int_var()) { - m_int_solver->notify_on_last_added_constraint(); - } update_column_type_and_bound(j, kind, right_side, ci); } else { @@ -1697,8 +1755,6 @@ void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_k unsigned term_j = it->second.internal_j(); mpq rs = right_side - m_terms[adjusted_term_index]->m_v; m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side)); - if (has_int_var()) - m_int_solver->notify_on_last_added_constraint(); update_column_type_and_bound(term_j, kind, rs, ci); } else { diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index b2585dea8..662d0bb3d 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -37,15 +37,13 @@ Revision History: #include "util/lp/stacked_value.h" #include "util/lp/stacked_vector.h" #include "util/lp/stacked_unordered_set.h" -#include "util/lp/iterator_on_pivot_row.h" #include "util/lp/implied_bound.h" #include "util/lp/bound_analyzer_on_row.h" -#include "util/lp/iterator_on_term_with_basis_var.h" -#include "util/lp/iterator_on_row.h" #include "util/lp/quick_xplain.h" #include "util/lp/conversion_helper.h" #include "util/lp/int_solver.h" #include "util/lp/nra_solver.h" +#include "util/lp/bound_propagator.h" namespace lp { @@ -60,30 +58,71 @@ class lar_solver : public column_namer { unsigned internal_j() const { return m_internal_j;} bool is_integer() const {return m_is_integer;} }; +#if Z3DEBUG_CHECK_UNIQUE_TERMS + struct term_hasher { + std::size_t operator()(const lar_term *t) const + { + using std::size_t; + using std::hash; + using std::string; + size_t seed = 0; + for (const auto& p : t->m_coeffs) { + hash_combine(seed, p); + } + return seed; + } + }; + + struct term_ls_comparer { + bool operator()(const lar_term *a, const lar_term* b) const + { + // a is contained in b + for (auto & p : a->m_coeffs) { + auto t = b->m_coeffs.find(p.first); + if (t == b->m_coeffs.end()) + return false; + if (p.second != t->second) + return false; + } + // zz is contained in b + for (auto & p : b->m_coeffs) { + auto t = a->m_coeffs.find(p.first); + if (t == a->m_coeffs.end()) + return false; + if (p.second != t->second) + return false; + } + return true; + } + }; + std::unordered_set m_set_of_terms; +#endif + //////////////////// fields ////////////////////////// - lp_settings m_settings; - lp_status m_status; - stacked_value m_simplex_strategy; - std::unordered_map m_ext_vars_to_columns; - vector m_columns_to_ext_vars_or_term_indices; - stacked_vector m_columns_to_ul_pairs; - vector m_constraints; + lp_settings m_settings; + lp_status m_status; + stacked_value m_simplex_strategy; + std::unordered_map m_ext_vars_to_columns; + vector m_columns_to_ext_vars_or_term_indices; + stacked_vector m_columns_to_ul_pairs; + vector m_constraints; +private: + stacked_value m_constraint_count; + // the set of column indices j such that bounds have changed for j + int_set m_columns_with_changed_bound; + int_set m_rows_with_changed_bounds; + int_set m_basic_columns_with_changed_cost; + stacked_value m_infeasible_column_index; // such can be found at the initialization step + stacked_value m_term_count; + vector m_terms; + const var_index m_terms_start_index; + indexed_vector m_column_buffer; + + public : const vector& constraints() const { return m_constraints; } -private: - stacked_value m_constraint_count; - // the set of column indices j such that bounds have changed for j - int_set m_columns_with_changed_bound; - int_set m_rows_with_changed_bounds; - int_set m_basic_columns_with_changed_cost; - stacked_value m_infeasible_column_index; // such can be found at the initialization step - stacked_value m_term_count; - vector m_terms; - const var_index m_terms_start_index; - indexed_vector m_column_buffer; -public: lar_core_solver m_mpq_lar_core_solver; private: std::function m_tracker_of_x_change; @@ -138,13 +177,17 @@ public: void add_new_var_to_core_fields_for_mpq(bool register_in_basis); - var_index add_term_undecided(const vector> & coeffs, - const mpq &m_v); // terms var_index add_term(const vector> & coeffs, const mpq &m_v); + var_index add_term_undecided(const vector> & coeffs, + const mpq &m_v); + + bool term_coeffs_are_ok(const vector> & coeffs, const mpq& v); + void push_and_register_term(lar_term* t); + void add_row_for_term(const lar_term * term, unsigned term_ext_index); void add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index); @@ -218,9 +261,6 @@ public: void calculate_implied_bounds_for_row(unsigned i, bound_propagator & bp); - - linear_combination_iterator * create_new_iter_from_term(unsigned term_index) const; - unsigned adjust_column_index_to_term_index(unsigned j) const; void propagate_bounds_on_a_term(const lar_term& t, bound_propagator & bp, unsigned term_offset); @@ -301,8 +341,6 @@ public: void detect_rows_of_bound_change_column_for_nbasic_column(unsigned j); - - void detect_rows_of_bound_change_column_for_nbasic_column_tableau(unsigned j); bool use_tableau() const; @@ -479,10 +517,12 @@ public: } bool bound_is_integer_for_integer_column(unsigned j, const mpq & right_side) const; - linear_combination_iterator * get_iterator_on_row(unsigned i) { - return m_mpq_lar_core_solver.m_r_solver.get_iterator_on_row(i); + + const row_strip & get_row(unsigned i) { + return A_r().m_rows[i]; } + unsigned get_base_column_in_row(unsigned row_index) const { return m_mpq_lar_core_solver.m_r_solver.get_base_column_in_row(row_index); } @@ -495,17 +535,20 @@ public: return m_columns_to_ul_pairs()[j].lower_bound_witness(); } - void subs_terms_for_debugging(lar_term& t) { - vector> pol; + void subs_term_columns(lar_term& t) { + vector> columns_to_subs; for (const auto & m : t.m_coeffs) { - pol.push_back(std::make_pair(m.second, adjust_column_index_to_term_index(m.first))); + unsigned tj = adjust_column_index_to_term_index(m.first); + if (tj == m.first) continue; + columns_to_subs.push_back(std::make_pair(m.first, tj)); + } + for (const auto & p : columns_to_subs) { + auto it = t.m_coeffs.find(p.first); + lp_assert(it != t.m_coeffs.end()); + mpq v = it->second; + t.m_coeffs.erase(it); + t.m_coeffs[p.second] = v; } - mpq v = t.m_v; - - vector> pol_after_subs; - substitute_terms_in_linear_expression(pol, pol_after_subs, v); - t.clear(); - t = lar_term(pol_after_subs, v); } bool inf_int_set_is_correct_for_column(unsigned j) const { diff --git a/src/util/lp/lar_term.h b/src/util/lp/lar_term.h index a2aee5109..519847848 100644 --- a/src/util/lp/lar_term.h +++ b/src/util/lp/lar_term.h @@ -1,22 +1,22 @@ /*++ -Copyright (c) 2017 Microsoft Corporation + Copyright (c) 2017 Microsoft Corporation -Module Name: + Module Name: - + -Abstract: + Abstract: - + -Author: + Author: - Lev Nachmanson (levnach) + Lev Nachmanson (levnach) -Revision History: + Revision History: ---*/ + --*/ #pragma once #include "util/lp/indexed_vector.h" namespace lp { @@ -85,11 +85,55 @@ struct lar_term { t.second.neg(); } + template + T apply(const vector& x) const { + T ret = T(m_v); + for (const auto & t : m_coeffs) { + ret += t.second * x[t.first]; + } + return ret; + } + void clear() { m_coeffs.clear(); m_v = zero_of_type(); } + struct ival { + unsigned m_var; + const mpq & m_coeff; + ival(unsigned var, const mpq & val) : m_var(var), m_coeff(val) { + } + unsigned var() const { return m_var;} + const mpq & coeff() const { return m_coeff; } + }; + struct const_iterator { + //fields + std::unordered_map::const_iterator m_it; + + typedef const_iterator self_type; + typedef ival value_type; + typedef ival reference; + // typedef std::pair* pointer; + typedef int difference_type; + typedef std::forward_iterator_tag iterator_category; + + reference operator*() const { + return ival(m_it->first, m_it->second); + } + + self_type operator++() { self_type i = *this; m_it++; return i; } + self_type operator++(int) { m_it++; return *this; } + + const_iterator(std::unordered_map::const_iterator it) : m_it(it) {} + bool operator==(const self_type &other) const { + return m_it == other.m_it; + } + bool operator!=(const self_type &other) const { return !(*this == other); } + }; + + const_iterator begin() const { return m_coeffs.begin();} + const_iterator end() const { return m_coeffs.end(); } }; } diff --git a/src/util/lp/linear_combination_iterator.h b/src/util/lp/linear_combination_iterator.h deleted file mode 100644 index 3ba30bdd8..000000000 --- a/src/util/lp/linear_combination_iterator.h +++ /dev/null @@ -1,66 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once -namespace lp { -template -struct linear_combination_iterator { - virtual bool next(T & a, unsigned & i) = 0; - virtual bool next(unsigned & i) = 0; - virtual void reset() = 0; - virtual linear_combination_iterator * clone() = 0; - virtual ~linear_combination_iterator(){} - virtual unsigned size() const = 0; - virtual bool is_reset() const = 0; -}; -template -struct linear_combination_iterator_on_vector : linear_combination_iterator { - vector> & m_vector; - int m_offset; - bool next(T & a, unsigned & i) { - if(static_cast(m_offset) >= m_vector.size()) - return false; - auto & p = m_vector[m_offset]; - a = p.first; - i = p.second; - m_offset++; - return true; - } - - bool next(unsigned & i) { - if(static_cast(m_offset) >= m_vector.size()) - return false; - auto & p = m_vector[m_offset]; - i = p.second; - m_offset++; - return true; - } - bool is_reset() const { return m_offset == 0;} - void reset() {m_offset = 0;} - linear_combination_iterator * clone() { - return new linear_combination_iterator_on_vector(m_vector); - } - linear_combination_iterator_on_vector(vector> & vec): - m_vector(vec), - m_offset(0) - {} - unsigned size() const { return m_vector.size(); } -}; - -} diff --git a/src/util/lp/linear_combination_iterator_on_std_vector.h b/src/util/lp/linear_combination_iterator_on_std_vector.h deleted file mode 100644 index e69de29bb..000000000 diff --git a/src/util/lp/lp_bound_propagator.cpp b/src/util/lp/lp_bound_propagator.cpp deleted file mode 100644 index 53218fced..000000000 --- a/src/util/lp/lp_bound_propagator.cpp +++ /dev/null @@ -1,67 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#include "util/lp/lar_solver.h" -namespace lp { - lp_bound_propagator::lp_bound_propagator(lar_solver & ls): - m_lar_solver(ls) {} -column_type lp_bound_propagator::get_column_type(unsigned j) const { - return m_lar_solver.m_mpq_lar_core_solver.m_column_types()[j]; -} -const impq & lp_bound_propagator::get_low_bound(unsigned j) const { - return m_lar_solver.m_mpq_lar_core_solver.m_r_low_bounds()[j]; -} -const impq & lp_bound_propagator::get_upper_bound(unsigned j) const { - return m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j]; -} -void lp_bound_propagator::try_add_bound(const mpq & v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) { - unsigned term_j = m_lar_solver.adjust_column_index_to_term_index(j); - mpq w = v; - if (term_j != j) { - j = term_j; - w += m_lar_solver.get_term(term_j).m_v; // when terms are turned into the columns they "lose" the right side, at this moment they aquire it back - } - lconstraint_kind kind = is_low? GE : LE; - if (strict) - kind = static_cast(kind / 2); - - if (!bound_is_interesting(j, kind, w)) - return; - unsigned k; // index to ibounds - if (is_low) { - if (try_get_val(m_improved_low_bounds, j, k)) { - auto & found_bound = m_ibounds[k]; - if (w > found_bound.m_bound || (w == found_bound.m_bound && found_bound.m_strict == false && strict)) - found_bound = implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict); - } else { - m_improved_low_bounds[j] = m_ibounds.size(); - m_ibounds.push_back(implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict)); - } - } else { // the upper bound case - if (try_get_val(m_improved_upper_bounds, j, k)) { - auto & found_bound = m_ibounds[k]; - if (w < found_bound.m_bound || (w == found_bound.m_bound && found_bound.m_strict == false && strict)) - found_bound = implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict); - } else { - m_improved_upper_bounds[j] = m_ibounds.size(); - m_ibounds.push_back(implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict)); - } - } -} -} diff --git a/src/util/lp/lp_bound_propagator.h b/src/util/lp/lp_bound_propagator.h deleted file mode 100644 index f3deaac5e..000000000 --- a/src/util/lp/lp_bound_propagator.h +++ /dev/null @@ -1,42 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once -#include "util/lp/lp_settings.h" -namespace lp { -class lar_solver; -class bound_propagator { - std::unordered_map m_improved_lower_bounds; // these maps map a column index to the corresponding index in ibounds - std::unordered_map m_improved_upper_bounds; - lar_solver & m_lar_solver; -public: - vector m_ibounds; -public: - lp_bound_propagator(lar_solver & ls); - column_type get_column_type(unsigned) const; - const impq & get_lower_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"; } -}; -} diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index 72e15e362..f1a02b372 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -28,8 +28,6 @@ Revision History: #include "util/lp/lu.h" #include "util/lp/permutation_matrix.h" #include "util/lp/column_namer.h" -#include "util/lp/iterator_on_row.h" -#include "util/lp/iterator_on_pivot_row.h" namespace lp { @@ -751,13 +749,6 @@ public: return m_iters_with_no_cost_growing; } - linear_combination_iterator * get_iterator_on_row(unsigned i) { - if (m_settings.use_tableau()) - return new iterator_on_row(m_A.m_rows[i]); - calculate_pivot_row(i); - return new iterator_on_pivot_row(m_pivot_row, m_basis[i]); - } - void calculate_pivot_row(unsigned i); unsigned get_base_column_in_row(unsigned row_index) const { return m_basis[row_index]; diff --git a/src/util/lp/lp_core_solver_base_def.h b/src/util/lp/lp_core_solver_base_def.h index 44a0a0c74..abadabf32 100644 --- a/src/util/lp/lp_core_solver_base_def.h +++ b/src/util/lp/lp_core_solver_base_def.h @@ -159,7 +159,7 @@ solve_Bd(unsigned entering) { m_columns_nz[entering] = m_ed.m_index.size(); lp_assert(m_ed.is_OK()); lp_assert(m_w.is_OK()); -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // auto B = get_B(*m_factorization, m_basis); // vector a(m_m()); // m_A.copy_column_to_vector(entering, a); @@ -978,8 +978,8 @@ template void lp_core_solver_base::pivot_fixed_v if (get_column_type(basic_j) != column_type::fixed) continue; T a; unsigned j; - auto * it = get_iterator_on_row(i); - while (it->next(a, j)) { + for (auto &c : m_A.m_rows[i]) { + j = c.var(); if (j == basic_j) continue; if (get_column_type(j) != column_type::fixed) { @@ -987,7 +987,6 @@ template void lp_core_solver_base::pivot_fixed_v break; } } - delete it; } } diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index 3c95e8ee0..0af8f0af6 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -37,7 +37,6 @@ Revision History: #include "util/lp/breakpoint.h" #include "util/lp/binary_heap_priority_queue.h" #include "util/lp/int_set.h" -#include "util/lp/iterator_on_row.h" namespace lp { // This core solver solves (Ax=b, lower_bound_values \leq x \leq upper_bound_values, maximize costs*x ) diff --git a/src/util/lp/lp_primal_simplex.h b/src/util/lp/lp_primal_simplex.h index 10efdcb79..d0f5d7d7d 100644 --- a/src/util/lp/lp_primal_simplex.h +++ b/src/util/lp/lp_primal_simplex.h @@ -26,7 +26,6 @@ Revision History: #include "util/lp/column_info.h" #include "util/lp/lp_primal_core_solver.h" #include "util/lp/lp_solver.h" -#include "util/lp/iterator_on_row.h" namespace lp { template class lp_primal_simplex: public lp_solver { diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h index 0a7e3f6ea..460097c58 100644 --- a/src/util/lp/lp_settings.h +++ b/src/util/lp/lp_settings.h @@ -106,6 +106,8 @@ struct stats { unsigned m_cut_solver_true; unsigned m_cut_solver_false; unsigned m_cut_solver_undef; + unsigned m_gcd_calls; + unsigned m_gcd_conflicts; stats() { reset(); } void reset() { memset(this, 0, sizeof(*this)); } }; @@ -226,9 +228,8 @@ public: backup_costs(true), column_number_threshold_for_using_lu_in_lar_solver(4000), m_int_branch_cut_gomory_threshold(4), - m_int_branch_cut_solver(4), + m_int_branch_cut_solver(8), m_run_gcd_test(true), - m_cut_solver_bound_propagation_factor(5), m_cut_solver_cycle_on_var(10) {} @@ -339,7 +340,6 @@ public: unsigned m_int_branch_cut_gomory_threshold; unsigned m_int_branch_cut_solver; bool m_run_gcd_test; - unsigned m_cut_solver_bound_propagation_factor; unsigned m_cut_solver_cycle_on_var; }; // end of lp_settings class diff --git a/src/util/lp/lp_solver.h b/src/util/lp/lp_solver.h index fb1f92c24..f1fa15e00 100644 --- a/src/util/lp/lp_solver.h +++ b/src/util/lp/lp_solver.h @@ -28,7 +28,6 @@ Revision History: #include "util/lp/static_matrix.h" #include "util/lp/lp_core_solver_base.h" #include "util/lp/scaler.h" -#include "util/lp/linear_combination_iterator.h" #include "util/lp/bound_analyzer_on_row.h" namespace lp { enum lp_relation { diff --git a/src/util/lp/lp_utils.h b/src/util/lp/lp_utils.h index b7a896c11..fca4f782f 100644 --- a/src/util/lp/lp_utils.h +++ b/src/util/lp/lp_utils.h @@ -38,7 +38,7 @@ bool contains(const std::unordered_map & map, const A& key) { #ifdef lp_for_z3 #ifdef Z3DEBUG -#define LEAN_DEBUG 1 +#define Z3DEBUG 1 #endif namespace lp { diff --git a/src/util/lp/lu.cpp b/src/util/lp/lu.cpp index 58da60d81..4e543a73b 100644 --- a/src/util/lp/lu.cpp +++ b/src/util/lp/lu.cpp @@ -39,7 +39,7 @@ template lp::mpq lp::dot_product(vector const&, vect template void lp::init_factorization(lp::lu*&, lp::static_matrix&, vector&, lp::lp_settings&); template void lp::init_factorization(lp::lu*&, lp::static_matrix&, vector&, lp::lp_settings&); template void lp::init_factorization >(lp::lu >*&, lp::static_matrix >&, vector&, lp::lp_settings&); -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG template void lp::print_matrix(lp::sparse_matrix&, std::ostream & out); template void lp::print_matrix(lp::static_matrix&, std::ostream&); template void lp::print_matrix >(lp::static_matrix >&, std::ostream&); diff --git a/src/util/lp/lu.h b/src/util/lp/lu.h index 76a90a205..9713936a1 100644 --- a/src/util/lp/lu.h +++ b/src/util/lp/lu.h @@ -34,7 +34,7 @@ Revision History: #include "util/lp/square_dense_submatrix.h" #include "util/lp/dense_matrix.h" namespace lp { -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG template // print the nr x nc submatrix at the top left corner void print_submatrix(sparse_matrix & m, unsigned mr, unsigned nc); @@ -113,7 +113,7 @@ public: #endif m_i = p.apply_reverse(m_i); -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(*this == deb); #endif } diff --git a/src/util/lp/lu_def.h b/src/util/lp/lu_def.h index e42d2e050..f8abdb2a5 100644 --- a/src/util/lp/lu_def.h +++ b/src/util/lp/lu_def.h @@ -25,7 +25,7 @@ Revision History: #include "util/debug.h" #include "util/lp/lu.h" namespace lp { -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG template // print the nr x nc submatrix at the top left corner void print_submatrix(sparse_matrix & m, unsigned mr, unsigned nc, std::ostream & out) { vector> A; @@ -138,12 +138,12 @@ lu::lu(static_matrix const & A, m_row_eta_work_vector(A.row_count()), m_refactor_counter(0) { lp_assert(!(numeric_traits::precise() && settings.use_tableau())); -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG debug_test_of_basis(A, basis); #endif ++m_settings.st().m_num_factorizations; create_initial_factorization(); -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(check_correctness()); #endif } diff --git a/src/util/lp/monomial.h b/src/util/lp/monomial.h new file mode 100644 index 000000000..5aa1eeb0c --- /dev/null +++ b/src/util/lp/monomial.h @@ -0,0 +1,33 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once +namespace lp { +struct monomial { + mpq m_coeff; // the coefficient of the monomial + var_index m_var; // the variable index +public: + monomial(const mpq& coeff, var_index var) : m_coeff(coeff), m_var(var) {} + monomial(var_index var) : monomial(one_of_type(), var) {} + const mpq & coeff() const { return m_coeff; } + mpq & coeff() { return m_coeff; } + var_index var() const { return m_var; } + std::pair to_pair() const { return std::make_pair(coeff(), var());} +}; +} diff --git a/src/util/lp/permutation_matrix.h b/src/util/lp/permutation_matrix.h index 1639b0cf1..77fb14cc3 100644 --- a/src/util/lp/permutation_matrix.h +++ b/src/util/lp/permutation_matrix.h @@ -28,7 +28,7 @@ Revision History: #include "util/lp/matrix.h" #include "util/lp/tail_matrix.h" namespace lp { -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG inline bool is_even(int k) { return (k/2)*2 == k; } #endif diff --git a/src/util/lp/permutation_matrix_def.h b/src/util/lp/permutation_matrix_def.h index c6a9dcfcd..bb2ca6c2a 100644 --- a/src/util/lp/permutation_matrix_def.h +++ b/src/util/lp/permutation_matrix_def.h @@ -74,7 +74,7 @@ void permutation_matrix::apply_from_left(vector & w, lp_settings & ) { while (i-- > 0) { w[i] = m_X_buffer[i]; } -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(vectors_are_equal(deb_w, w, row_count())); // delete [] deb_w; #endif @@ -109,7 +109,7 @@ template void permutation_matrix::apply_from_righ for (unsigned i = 0; i < size(); i++) { w[i] = m_T_buffer[i]; } -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(vectors_are_equal(deb_w, w, row_count())); // delete [] deb_w; #endif @@ -133,7 +133,7 @@ template void permutation_matrix::apply_from_righ w.set_value(buffer[i], pj); } lp_assert(w.is_OK()); -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG lp_assert(vectors_are_equal(wcopy, w.m_data)); #endif } @@ -180,7 +180,7 @@ void permutation_matrix::apply_reverse_from_left(indexed_vector & w) { w[j] = t[i]; w.m_index[i] = j; } -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(vectors_are_equal(deb_w, w.m_data, row_count())); // delete [] deb_w; #endif diff --git a/src/util/lp/polynomial.h b/src/util/lp/polynomial.h new file mode 100644 index 000000000..e88cc4a9c --- /dev/null +++ b/src/util/lp/polynomial.h @@ -0,0 +1,106 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once +namespace lp { +struct polynomial { + static mpq m_local_zero; + // the polynomial evaluates to m_coeffs + m_a + vector m_coeffs; + mpq m_a; // the free coefficient + polynomial(const vector& p, const mpq & a) : m_coeffs(p), m_a(a) {} + polynomial(const vector& p) : polynomial(p, zero_of_type()) {} + polynomial(): m_a(zero_of_type()) {} + polynomial(const polynomial & p) : m_coeffs(p.m_coeffs), m_a(p.m_a) {} + + const mpq & coeff(var_index j) const { + for (const auto & t : m_coeffs) { + if (j == t.var()) { + return t.coeff(); + } + } + return m_local_zero; + } + + polynomial & operator+=(const polynomial & p) { + m_a += p.m_a; + for (const auto & t: p.m_coeffs) + *this += monomial(t.coeff(), t.var()); + return *this; + } + + void add(const mpq & c, const polynomial &p) { + m_a += p.m_a * c; + + for (const auto & t: p.m_coeffs) + *this += monomial(c * t.coeff(), t.var()); + } + + void clear() { + m_coeffs.clear(); + m_a = zero_of_type(); + } + + bool is_empty() const { return m_coeffs.size() == 0 && numeric_traits::is_zero(m_a); } + + unsigned number_of_monomials() const { return m_coeffs.size();} + + void add(const monomial &m ){ + if (is_zero(m.coeff())) return; + for (unsigned k = 0; k < m_coeffs.size(); k++) { + auto & l = m_coeffs[k]; + if (m.var() == l.var()) { + l.coeff() += m.coeff(); + if (l.coeff() == 0) + m_coeffs.erase(m_coeffs.begin() + k); + return; + } + } + m_coeffs.push_back(m); + lp_assert(is_correct()); + } + + polynomial & operator+=(const monomial &m ){ + add(m); + return *this; + } + + polynomial & operator+=(const mpq &c ){ + m_a += c; + return *this; + } + + + bool is_correct() const { + std::unordered_set s; + for (auto & l : m_coeffs) { + if (l.coeff() == 0) + return false; + s.insert(l.var()); + } + return m_coeffs.size() == s.size(); + } + + bool var_coeff_is_unit(unsigned j) const { + const mpq & a = coeff(j); + return a == 1 || a == -1; + } + const vector & coeffs() const { return m_coeffs; } +}; +} diff --git a/src/util/lp/random_updater.h b/src/util/lp/random_updater.h index 441038c00..e5067f2e5 100644 --- a/src/util/lp/random_updater.h +++ b/src/util/lp/random_updater.h @@ -25,7 +25,6 @@ Revision History: #include #include #include "util/lp/lp_settings.h" -#include "util/lp/linear_combination_iterator.h" // see http://research.microsoft.com/projects/z3/smt07.pdf // The class searches for a feasible solution with as many different values of variables as it can find namespace lp { diff --git a/src/util/lp/random_updater_def.h b/src/util/lp/random_updater_def.h index 664ad4eca..419335aa3 100644 --- a/src/util/lp/random_updater_def.h +++ b/src/util/lp/random_updater_def.h @@ -83,7 +83,7 @@ void random_updater::add_column_to_sets(unsigned j) { add_value(m_lar_solver.get_core_solver().m_r_x[j]); } else { unsigned row = m_lar_solver.get_core_solver().m_r_heading[j]; - for (auto row_c : m_lar_solver.get_core_solver().m_r_A.m_rows[row]) { + for (auto & row_c : m_lar_solver.get_core_solver().m_r_A.m_rows[row]) { unsigned cj = row_c.m_j; if (m_lar_solver.get_core_solver().m_r_heading[cj] < 0) { m_var_set.insert(cj); diff --git a/src/util/lp/row_eta_matrix_def.h b/src/util/lp/row_eta_matrix_def.h index 272b6313e..46b187a24 100644 --- a/src/util/lp/row_eta_matrix_def.h +++ b/src/util/lp/row_eta_matrix_def.h @@ -33,7 +33,7 @@ void row_eta_matrix::apply_from_left(vector & w, lp_settings &) { w_at_row += w[it.first] * it.second; } // w[m_row] = w_at_row; - // #ifdef LEAN_DEBUG + // #ifdef Z3DEBUG // lp_assert(vectors_are_equal(clone_w, w, m_dimension)); // delete [] clone_w; // #endif @@ -95,7 +95,7 @@ void row_eta_matrix::apply_from_right(vector & w) { for (auto & it : m_row_vector.m_data) { w[it.first] += w_row * it.second; } -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(vectors_are_equal(clone_w, w, m_dimension)); // delete clone_w; #endif @@ -144,7 +144,7 @@ void row_eta_matrix::apply_from_right(indexed_vector & w) { } } } -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(vectors_are_equal(wcopy, w.m_data)); #endif @@ -165,7 +165,7 @@ void row_eta_matrix::conjugate_by_permutation(permutation_matrix & p columns.push_back(it.first); for (unsigned i = static_cast(columns.size()); i-- > 0;) m_row_vector.m_data[i].first = p.get_rev(columns[i]); -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(deb == *this); #endif } diff --git a/src/util/lp/scaler_def.h b/src/util/lp/scaler_def.h index 4a0e51da6..3f21f1953 100644 --- a/src/util/lp/scaler_def.h +++ b/src/util/lp/scaler_def.h @@ -220,18 +220,18 @@ template void scaler::scale_row(unsigned i) { if (numeric_traits::is_zero(row_max)) { return; } - if (numeric_traits::get_double(row_max) < m_scaling_minimum) { + if (row_max < m_scaling_minimum) { do { alpha *= 2; row_max *= 2; - } while (numeric_traits::get_double(row_max) < m_scaling_minimum); + } while (row_max < m_scaling_minimum); m_A.multiply_row(i, alpha); m_b[i] *= alpha; - } else if (numeric_traits::get_double(row_max) > m_scaling_maximum) { + } else if (row_max > m_scaling_maximum) { do { alpha /= 2; row_max /= 2; - } while (numeric_traits::get_double(row_max) > m_scaling_maximum); + } while (row_max > m_scaling_maximum); m_A.multiply_row(i, alpha); m_b[i] *= alpha; } @@ -245,16 +245,16 @@ template void scaler::scale_column(unsigned i) return; // the column has zeros only } std::cout << "f"; - if (numeric_traits::get_double(column_max) < m_scaling_minimum) { + if (column_max < m_scaling_minimum) { do { alpha *= 2; column_max *= 2; - } while (numeric_traits::get_double(column_max) < m_scaling_minimum); - } else if (numeric_traits::get_double(column_max) > m_scaling_maximum) { + } while (column_max < m_scaling_minimum); + } else if (column_max > m_scaling_maximum) { do { alpha /= 2; column_max /= 2; - } while (numeric_traits::get_double(column_max) > m_scaling_maximum); + } while (column_max > m_scaling_maximum); } else { return; } diff --git a/src/util/lp/sparse_matrix.cpp b/src/util/lp/sparse_matrix.cpp index cdf255695..4ce8e64ef 100644 --- a/src/util/lp/sparse_matrix.cpp +++ b/src/util/lp/sparse_matrix.cpp @@ -82,7 +82,7 @@ template void sparse_matrix>::double_solve_U_y(index template void sparse_matrix >::double_solve_U_y >(indexed_vector>&, const lp_settings&); template void lp::sparse_matrix::solve_U_y_indexed_only(lp::indexed_vector&, const lp_settings&, vector &); template void lp::sparse_matrix::solve_U_y_indexed_only(lp::indexed_vector&, const lp_settings &, vector &); -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG template bool sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; template bool sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; template bool sparse_matrix >::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; diff --git a/src/util/lp/sparse_matrix_def.h b/src/util/lp/sparse_matrix_def.h index 8892e5b14..39159243c 100644 --- a/src/util/lp/sparse_matrix_def.h +++ b/src/util/lp/sparse_matrix_def.h @@ -490,7 +490,7 @@ void sparse_matrix::solve_y_U_indexed(indexed_vector & y, const lp_sett } lp_assert(y.is_OK()); -#if 0 && LEAN_DEBUG +#if 0 && Z3DEBUG if (numeric_traits::precise() == false) lp_assert(vectors_are_equal(ycopy, y.m_data)); #endif @@ -700,7 +700,7 @@ void sparse_matrix::solve_U_y_indexed_only(indexed_vector & y, const lp } lp_assert(y.is_OK()); -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // dense_matrix deb(this); // vector clone_y(y.m_data); // deb.apply_from_left(clone_y); diff --git a/src/util/lp/square_dense_submatrix_def.h b/src/util/lp/square_dense_submatrix_def.h index 5608c6e39..ca00285bb 100644 --- a/src/util/lp/square_dense_submatrix_def.h +++ b/src/util/lp/square_dense_submatrix_def.h @@ -300,7 +300,7 @@ void square_dense_submatrix::apply_from_left_to_vector(vector & w) { } template bool square_dense_submatrix::is_L_matrix() const { -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG lp_assert(m_row_permutation.is_identity()); for (unsigned i = 0; i < m_parent->dimension(); i++) { if (i < m_index_start) { @@ -341,7 +341,7 @@ template void square_dense_submatrix::apply_from_ t[adjust_column_inverse(j)] = column_by_vector_product(j, w); } w = t; -#ifdef LEAN_DEBUG +#ifdef Z3DEBUG // lp_assert(vector_are_equal(deb_w, w)); #endif } diff --git a/src/util/lp/static_matrix.h b/src/util/lp/static_matrix.h index a1c3b64a4..2f127b84a 100644 --- a/src/util/lp/static_matrix.h +++ b/src/util/lp/static_matrix.h @@ -26,7 +26,6 @@ Revision History: #include "util/lp/sparse_vector.h" #include "util/lp/indexed_vector.h" #include "util/lp/permutation_matrix.h" -#include "util/lp/linear_combination_iterator.h" #include namespace lp { @@ -35,20 +34,26 @@ struct column_cell { unsigned m_offset; // the offset of the element in the matrix row column_cell(unsigned i, unsigned offset) : m_i(i), m_offset(offset) { } + }; template struct row_cell { unsigned m_j; // points to the column unsigned m_offset; // offset in column + T m_value; row_cell(unsigned j, unsigned offset, T const & val) : m_j(j), m_offset(offset), m_value(val) { } const T & get_val() const { return m_value;} T & get_val() { return m_value;} void set_val(const T& v) { m_value = v; } - T m_value; + unsigned var() const { return m_j;} + const T & coeff() const { return m_value;} }; +template +using row_strip = vector>; + // each assignment for this matrix should be issued only once!!! template class static_matrix @@ -63,11 +68,10 @@ class static_matrix }; std::stack m_stack; public: - typedef vector> row_strip; typedef vector column_strip; vector m_vector_of_row_offsets; indexed_vector m_work_vector; - vector m_rows; + vector> m_rows; vector m_columns; // starting inner classes class ref { @@ -97,10 +101,6 @@ public: return m_rows[c.m_i][c.m_offset].get_val(); } - row_cell & get_row_cell(const column_cell & c) { - return m_rows[c.m_i][c.m_offset]; - } - column_cell & get_column_cell(const row_cell &rc) { return m_columns[rc.m_j][rc.m_offset]; } @@ -132,7 +132,7 @@ public: void add_columns_at_the_end(unsigned delta); void add_new_element(unsigned i, unsigned j, const T & v); - void add_row() {m_rows.push_back(row_strip());} + void add_row() {m_rows.push_back(row_strip());} void add_column() { m_columns.push_back(column_strip()); m_vector_of_row_offsets.push_back(-1); @@ -195,7 +195,7 @@ public: // void fix_row_indices_in_each_column_for_crossed_row(unsigned k); - void cross_out_row_from_columns(unsigned k, row_strip & row); + void cross_out_row_from_columns(unsigned k, row_strip & row); void cross_out_row_from_column(unsigned col, unsigned k); @@ -294,7 +294,7 @@ public: // pivot row i to row ii bool pivot_row_to_row_given_cell(unsigned i, column_cell& c, unsigned); - void scan_row_ii_to_offset_vector(const row_strip & rvals); + void scan_row_ii_to_offset_vector(const row_strip & rvals); void transpose_rows(unsigned i, unsigned ii) { auto t = m_rows[i]; @@ -313,46 +313,54 @@ public: } } + void fill_last_row_with_pivoting_loop_block(unsigned j, const vector & basis_heading) { + int row_index = basis_heading[j]; + if (row_index < 0) + return; + T & alpha = m_work_vector[j]; // the pivot alpha + if (is_zero(alpha)) + return; + + for (const auto & c : m_rows[row_index]) { + if (c.m_j == j) { + continue; + } + T & wv = m_work_vector.m_data[c.m_j]; + bool was_zero = is_zero(wv); + wv -= alpha * c.m_value; + if (was_zero) + m_work_vector.m_index.push_back(c.m_j); + else { + if (is_zero(wv)) { + m_work_vector.erase_from_index(c.m_j); + } + } + } + alpha = zero_of_type(); + m_work_vector.erase_from_index(j); + } - void fill_last_row_with_pivoting(linear_combination_iterator & it, const vector & basis_heading) { + + + template + void fill_last_row_with_pivoting(const term& row, + unsigned bj, // the index of the basis column + const vector & basis_heading) { lp_assert(numeric_traits::precise()); lp_assert(row_count() > 0); m_work_vector.resize(column_count()); T a; - unsigned j; - while (it.next(a, j)) { - m_work_vector.set_value(-a, j); // we use the form -it + 1 = 0 + // we use the form -it + 1 = 0 + m_work_vector.set_value(one_of_type(), bj); + for (auto p : row) { + m_work_vector.set_value(-p.coeff(), p.var()); // but take care of the basis 1 later } - it.reset(); - // not iterate with pivoting - while (it.next(j)) { - int row_index = basis_heading[j]; - if (row_index < 0) - continue; - - T & alpha = m_work_vector[j]; // the pivot alpha - if (is_zero(alpha)) - continue; - - for (const auto & c : m_rows[row_index]) { - if (c.m_j == j) { - continue; - } - T & wv = m_work_vector.m_data[c.m_j]; - bool was_zero = is_zero(wv); - wv -= alpha * c.m_value; - if (was_zero) - m_work_vector.m_index.push_back(c.m_j); - else { - if (is_zero(wv)) { - m_work_vector.erase_from_index(c.m_j); - } - } - } - alpha = zero_of_type(); - m_work_vector.erase_from_index(j); + // now iterate with pivoting + fill_last_row_with_pivoting_loop_block(bj, basis_heading); + for (auto p : row) { + fill_last_row_with_pivoting_loop_block(p.var(), basis_heading); } lp_assert(m_work_vector.is_OK()); unsigned last_row = row_count() - 1; @@ -382,5 +390,66 @@ public: } return ret; } + + struct column_cell_plus { + const column_cell & m_c; + const static_matrix& m_A; + // constructor + column_cell_plus(const column_cell & c, const static_matrix& A) : + m_c(c), m_A(A) {} + unsigned var() const { return m_c.m_i; } + const T & coeff() const { return m_A.m_rows[var()][m_c.m_offset].get_val(); } + + }; + + struct column_container { + unsigned m_j; // the column index + const static_matrix & m_A; + column_container(unsigned j, const static_matrix& A) : m_j(j), m_A(A) { + } + struct const_iterator { + // fields + const column_cell *m_c; + const static_matrix& m_A; + + //typedefs + + + typedef const_iterator self_type; + typedef column_cell_plus value_type; + typedef const column_cell_plus reference; + // typedef const column_cell* pointer; + typedef int difference_type; + typedef std::forward_iterator_tag iterator_category; + + reference operator*() const { + return column_cell_plus(*m_c, m_A); + } + self_type operator++() { self_type i = *this; m_c++; return i; } + self_type operator++(int) { m_c++; return *this; } + + const_iterator(const column_cell* it, const static_matrix& A) : + m_c(it), + m_A(A) + {} + bool operator==(const self_type &other) const { + return m_c == other.m_c; + } + bool operator!=(const self_type &other) const { return !(*this == other); } + }; + + const_iterator begin() const { + return const_iterator(m_A.m_columns[m_j].begin(), m_A); + } + + const_iterator end() const { + return const_iterator(m_A.m_columns[m_j].end(), m_A); + } + }; + + column_container column(unsigned j) const { + return column_container(j, *this); + } + }; } diff --git a/src/util/lp/static_matrix_def.h b/src/util/lp/static_matrix_def.h index 6797d265e..a246a5aea 100644 --- a/src/util/lp/static_matrix_def.h +++ b/src/util/lp/static_matrix_def.h @@ -27,7 +27,7 @@ template void static_matrix::init_row_columns(unsigned m, unsigned n) { lp_assert(m_rows.size() == 0 && m_columns.size() == 0); for (unsigned i = 0; i < m; i++){ - m_rows.push_back(row_strip()); + m_rows.push_back(row_strip()); } for (unsigned j = 0; j < n; j++){ m_columns.push_back(column_strip()); @@ -35,7 +35,7 @@ void static_matrix::init_row_columns(unsigned m, unsigned n) { } -template void static_matrix::scan_row_ii_to_offset_vector(const row_strip & rvals) { +template void static_matrix::scan_row_ii_to_offset_vector(const row_strip & rvals) { for (unsigned j = 0; j < rvals.size(); j++) m_vector_of_row_offsets[rvals[j].m_j] = j; } @@ -166,7 +166,7 @@ template std::set> static_matrix::get_domain() { std::set> ret; for (unsigned i = 0; i < m_rows.size(); i++) { - for (auto it : m_rows[i]) { + for (auto &it : m_rows[i]) { ret.insert(std::make_pair(i, it.m_j)); } } @@ -296,7 +296,7 @@ template void static_matrix::fix_row_indices_i } } -template void static_matrix::cross_out_row_from_columns(unsigned k, row_strip & row) { +template void static_matrix::cross_out_row_from_columns(unsigned k, row_strip & row) { for (auto & t : row) { cross_out_row_from_column(t.m_j, k); } diff --git a/src/util/lp/test_bound_analyzer.h b/src/util/lp/test_bound_analyzer.h index b11fe6d94..3765f480e 100644 --- a/src/util/lp/test_bound_analyzer.h +++ b/src/util/lp/test_bound_analyzer.h @@ -17,9 +17,9 @@ Revision History: --*/ +#if 0 #pragma once #include "util/vector.h" -#include "util/lp/linear_combination_iterator.h" #include "util/lp/implied_bound.h" #include "util/lp/lp_settings.h" #include @@ -273,3 +273,4 @@ public : }; } +#endif