diff --git a/src/math/simplex/simplex.h b/src/math/simplex/simplex.h index 0cb4852f2..974aa7cfb 100644 --- a/src/math/simplex/simplex.h +++ b/src/math/simplex/simplex.h @@ -36,6 +36,7 @@ Notes: #include "mpq_inf.h" #include "heap.h" #include "lbool.h" +#include "uint_set.h" namespace simplex { @@ -84,7 +85,7 @@ namespace simplex { }; static const var_t null_var = UINT_MAX; - matrix M; + mutable matrix M; unsigned m_max_iterations; volatile bool m_cancel; var_heap m_to_patch; @@ -93,14 +94,17 @@ namespace simplex { vector m_vars; svector m_row2base; bool m_bland; + unsigned m_blands_rule_threshold; random_gen m_random; + uint_set m_left_basis; public: simplex(): m_max_iterations(UINT_MAX), m_cancel(false), m_to_patch(1024), - m_bland(false) {} + m_bland(false), + m_blands_rule_threshold(1000) {} typedef typename matrix::row row; typedef typename matrix::row_iterator row_iterator; @@ -117,7 +121,7 @@ namespace simplex { void set_cancel(bool f) { m_cancel = f; } void set_max_iterations(unsigned m) { m_max_iterations = m; } lbool make_feasible(); - lbool optimize(var_t var); + lbool minimize(var_t var); eps_numeral const& get_value(var_t v); void display(std::ostream& out) const; @@ -128,17 +132,26 @@ namespace simplex { var_t select_smallest_var() { return m_to_patch.empty()?null_var:m_to_patch.erase_min(); } var_t select_error_var(bool least); // row get_infeasible_row() { } - void check_blands_rule(var_t v) { } + void check_blands_rule(var_t v, unsigned& num_repeated); bool make_var_feasible(var_t x_i); void update_and_pivot(var_t x_i, var_t x_j, numeral const& a_ij, eps_numeral const& new_value); void update_value(var_t v, eps_numeral const& delta); void update_value_core(var_t v, eps_numeral const& delta); + void pivot(var_t x_i, var_t x_j, numeral const& a_ij); + void move_to_bound(var_t x, bool to_lower); var_t select_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij); - var_t select_blands_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij); - template - var_t select_pivot_core(var_t x_i, scoped_numeral& out_a_ij); + var_t select_pivot_blands(var_t x_i, bool is_below, scoped_numeral& out_a_ij); + var_t select_pivot_core(var_t x_i, bool is_below, scoped_numeral& out_a_ij); int get_num_non_free_dep_vars(var_t x_j, int best_so_far); + var_t pick_var_to_leave(var_t x_j, bool inc, scoped_eps_numeral& gain, scoped_numeral& new_a_ij); + + + void select_pivot_primal(var_t v, var_t& x_i, var_t& x_j, scoped_numeral& a_ij, bool& inc); + + + bool at_lower(var_t v) const; + bool at_upper(var_t v) const; bool below_lower(var_t v) const; bool above_upper(var_t v) const; bool above_lower(var_t v) const; @@ -150,6 +163,7 @@ namespace simplex { bool is_base(var_t x) const { return m_vars[x].m_is_base; } bool well_formed() const; + bool is_feasible() const; }; }; diff --git a/src/math/simplex/simplex_def.h b/src/math/simplex/simplex_def.h index 0e05de291..26c7e2e03 100644 --- a/src/math/simplex/simplex_def.h +++ b/src/math/simplex/simplex_def.h @@ -44,7 +44,8 @@ namespace simplex { template void simplex::del_row(row const& r) { - m_is_base[m_row2base[r.id()]] = false; + m_vars[m_row2base[r.id()]].m_is_base = false; + m_row2base[r.id()] = null_var; M.del(r); } @@ -102,11 +103,11 @@ namespace simplex { var_info const& vi = m_vars[i]; out << "v" << i << " "; if (vi.m_is_base) out << "b:" << vi.m_base2row << " "; - em.display(out, vi.m_value); + out << em.to_string(vi.m_value); out << " ["; - if (vi.m_lower_valid) em.display(out, vi.m_lower) else out << "-oo"; + if (vi.m_lower_valid) out << em.to_string(vi.m_lower); else out << "-oo"; out << ":"; - if (vi.m_upper_valid) em.display(out, vi.m_upper) else out << "oo"; + if (vi.m_upper_valid) out << em.to_string(vi.m_upper); else out << "oo"; out << "]"; out << "\n"; } @@ -122,13 +123,15 @@ namespace simplex { template lbool simplex::make_feasible() { + m_left_basis.reset(); unsigned num_iterations = 0; + unsigned num_repeated = 0; var_t v = null_var; while ((v = select_var_to_fix()) != null_var) { if (m_cancel || num_iterations > m_max_iterations) { return l_undef; } - check_blands_rule(v); + check_blands_rule(v, num_repeated); if (!make_var_feasible(v)) { return l_false; } @@ -157,7 +160,6 @@ namespace simplex { SASSERT(is_base(x_i)); SASSERT(!is_base(x_j)); var_info& x_iI = m_vars[x_i]; - var_info& x_jI = m_vars[x_j]; scoped_eps_numeral theta(em); theta = x_iI.m_value; theta -= new_value; @@ -166,6 +168,13 @@ namespace simplex { em.div(theta, a_ij, theta); update_value(x_j, theta); SASSERT(em.eq(x_iI.m_value, new_value)); + pivot(x_i, x_j, a_ij); + } + + template + void simplex::pivot(var_t x_i, var_t x_j, numeral const& a_ij) { + var_info& x_iI = m_vars[x_i]; + var_info& x_jI = m_vars[x_j]; unsigned r_i = x_iI.m_base2row; x_jI.m_base2row = r_i; m_row2base[r_i] = x_j; @@ -196,6 +205,9 @@ namespace simplex { template void simplex::update_value(var_t v, eps_numeral const& delta) { + if (em.is_zero(delta)) { + return; + } update_value_core(v, delta); col_iterator it = M.col_begin(v), end = M.col_end(v); @@ -249,6 +261,18 @@ namespace simplex { return !vi.m_upper_valid || em.lt(vi.m_value, vi.m_upper); } + template + bool simplex::at_lower(var_t v) const { + var_info const& vi = m_vars[v]; + return vi.m_lower_valid && em.eq(vi.m_value, vi.m_lower); + } + + template + bool simplex::at_upper(var_t v) const { + var_info const& vi = m_vars[v]; + return vi.m_upper_valid && em.eq(vi.m_value, vi.m_upper); + } + template bool simplex::make_var_feasible(var_t x_i) { scoped_numeral a_ij(m); @@ -274,20 +298,15 @@ namespace simplex { } /** - \brief Wrapper for select_blands_pivot_core and select_pivot_core + \brief Wrapper for select_pivot_blands and select_pivot_core */ template typename simplex::var_t simplex::select_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij) { if (m_bland) { - return select_blands_pivot(x_i, is_below, out_a_ij); - } - if (is_below) { - return select_pivot_core(x_i, out_a_ij); - } - else { - return select_pivot_core(x_i, out_a_ij); + return select_pivot_blands(x_i, is_below, out_a_ij); } + return select_pivot_core(x_i, is_below, out_a_ij); } /** @@ -300,9 +319,8 @@ namespace simplex { bound (above its upper bound). */ template - template typename simplex::var_t - simplex::select_pivot_core(var_t x_i, scoped_numeral & out_a_ij) { + simplex::select_pivot_core(var_t x_i, bool is_below, scoped_numeral & out_a_ij) { SASSERT(is_base(x_i)); var_t max = get_num_vars(); var_t result = max; @@ -315,11 +333,13 @@ namespace simplex { for (; it != end; ++it) { var_t x_j = it->m_var; + if (x_i == x_j) continue; numeral const & a_ij = it->m_coeff; bool is_neg = is_below ? m.is_neg(a_ij) : m.is_pos(a_ij); bool is_pos = !is_neg; - if (x_i != x_j && ((is_pos && above_lower(x_j)) || (is_neg && below_upper(x_j)))) { + bool can_pivot = ((is_pos && above_lower(x_j)) || (is_neg && below_upper(x_j))); + if (can_pivot) { int num = get_num_non_free_dep_vars(x_j, best_so_far); unsigned col_sz = M.column_size(x_j); if (num < best_so_far || (num == best_so_far && col_sz < best_col_sz)) { @@ -360,16 +380,15 @@ namespace simplex { return result; } - /** \brief Using Bland's rule, select a variable x_j in the row r defining the base var x_i, - s.t. x_j can be used to patch the error in x_i. Return null_theory_var + s.t. x_j can be used to patch the error in x_i. Return null_var if there is no x_j. Otherwise, return x_j and store its coefficient in out_a_ij. */ template typename simplex::var_t - simplex::select_blands_pivot(var_t x_i, bool is_below, scoped_numeral & out_a_ij) { + simplex::select_pivot_blands(var_t x_i, bool is_below, scoped_numeral & out_a_ij) { SASSERT(is_base(x_i)); unsigned max = get_num_vars(); var_t result = max; @@ -379,8 +398,7 @@ namespace simplex { var_t x_j = it->m_var; numeral const & a_ij = it->m_coeff; bool is_neg = is_below ? m.is_neg(a_ij) : m.is_pos(a_ij); - bool is_pos = !is_neg; - if (x_i != x_j && ((is_pos && above_lower(x_j)) || (is_neg && below_upper(x_j)))) { + if (x_i != x_j && ((!is_neg && above_lower(x_j)) || (is_neg && below_upper(x_j)))) { SASSERT(!is_base(x_j)); if (x_j < result) { result = x_j; @@ -392,16 +410,230 @@ namespace simplex { } template - lbool simplex::optimize(var_t v) { - // TBD SASSERT(is_feasible()); - // pick row for v and check if primal - // bounds are slack. - // return l_false for unbounded. - // return l_undef for canceled. - // return l_true for optimal. + lbool simplex::minimize(var_t v) { + + // minimize v, such that tableau is feasible. + // Assume there are no bounds on v. + // Let k*v + c*x = 0 e.g, maximize c*x over + // tableau constraints: + // + // max { c*x | A*x = 0 and l <= x <= u } + // + // start with feasible assigment + // A*x0 = 0 and l <= x0 <= u + // + // Identify pivot: i, j: such that x_i is base, + // there is a row k1*x_i + k2*x_j + R = 0 + // and a delta such that: + // + // x_i' <- x_i + delta + // x_j' <- x_j - delta*k1/k2 + // l_i <= x_i' <= u_i + // l_j <= x_j' <= u_j + // and c*x' > c*x + // e.g., c*x := c_i*x_i + c_j*x_j + ... + // and c_i*delta > c_j*delta*k1/k2 + // and x_i < u_i (if delta > 0), l_i < x_i (if delta < 0) + // and l_j < x_j (if delta > 0), x_j < u_j (if delta < 0) + // + // update all rows, including c*x, using the pivot. + // + // If there is c_i*x_i in c*x such that c_i > 0 + // and upper_i = oo and complementary lower_j = -oo + // then the objective is unbounded. + // + // There is a singularity if there is a pivot such that + // c_i*delta == c_j*delta*k1/k2, e.g., nothing is improved, + // pivot, but use bland's rule to ensure + // convergence in the limit. + // + + SASSERT(is_feasible()); + scoped_eps_numeral delta(em); + scoped_numeral a_ij(m); + var_t x_i, x_j; + bool inc; + + while (true) { + if (m_cancel) { + return l_undef; + } + select_pivot_primal(v, x_i, x_j, a_ij, inc); + if (x_j == null_var) { + // optimal + return l_true; + } + var_info& vj = m_vars[x_j]; + if (x_i == null_var) { + if (inc && vj.m_upper_valid) { + delta = vj.m_upper; + delta -= vj.m_value; + update_value(x_j, delta); + } + else if (!inc && vj.m_lower_valid) { + delta = vj.m_lower; + delta -= vj.m_value; + update_value(x_j, delta); + } + else { + // unbounded + return l_false; + } + continue; + } + + // TBD: Change the value of x_j directly without pivoting: + // + // if (!vj.is_fixed() && vj.bounded() && gain >= upper - lower) { + // + // } + // + + pivot(x_i, x_j, a_ij); + move_to_bound(x_i, inc == m.is_pos(a_ij)); + } return l_true; } + template + void simplex::move_to_bound(var_t x, bool to_lower) { + scoped_eps_numeral delta(em), delta2(em); + var_info& vi = m_vars[x]; + if (to_lower) { + em.sub(vi.m_value, vi.m_lower, delta); + } + else { + em.sub(vi.m_upper, vi.m_value, delta); + } + col_iterator it = M.col_begin(x), end = M.col_end(x); + for (; it != end && is_pos(delta); ++it) { + var_t s = m_row2base[it.get_row().id()]; + var_info& vs = m_vars[s]; + numeral const& coeff = it.get_row_entry().m_coeff; + SASSERT(!m.is_zero(coeff)); + bool inc_s = (m.is_pos(coeff) == to_lower); + eps_numeral const* bound = 0; + if (inc_s && vs.m_upper_valid) { + bound = &vs.m_upper; + } + else if (!inc_s && vs.m_lower_valid) { + bound = &vs.m_lower; + } + if (bound) { + em.sub(*bound, vs.m_value, delta2); + em.div(delta2, coeff, delta2); + abs(delta2); + if (delta2 < delta) { + delta = delta2; + } + } + } + if (to_lower) { + delta.neg(); + } + update_value(x, delta); + } + + template + void simplex::select_pivot_primal(var_t v, var_t& x_i, var_t& x_j, scoped_numeral& a_ij, bool& inc) { + row r(m_vars[v].m_base2row); + row_iterator it = M.row_begin(r), end = M.row_end(r); + + scoped_eps_numeral gain(em), new_gain(em); + scoped_numeral new_a_ij(m); + x_i = null_var; + x_j = null_var; + inc = false; + + for (; it != end; ++it) { + var_t x = it->m_var; + if (x == v) continue; + bool is_pos = m.is_pos(it->m_coeff); + if ((is_pos && at_upper(x)) || (!is_pos && at_lower(x))) { + continue; // variable cannot be used for improving bounds. + // TBD check? + } + var_t y = pick_var_to_leave(x, is_pos, new_gain, new_a_ij); + if (y == null_var) { + // unbounded. + x_i = y; + x_j = x; + inc = is_pos; + a_ij = new_a_ij; + break; + } + bool better = + (new_gain > gain) || + ((is_zero(new_gain) && is_zero(gain) && (x_i == null_var || y < x_i))); + + if (better) { + x_i = y; + x_j = x; + inc = is_pos; + gain = new_gain; + a_ij = new_a_ij; + } + } + } + + template + typename simplex::var_t + simplex::pick_var_to_leave( + var_t x_j, bool inc, + scoped_eps_numeral& gain, scoped_numeral& new_a_ij) { + var_t x_i = null_var; + gain.reset(); + scoped_eps_numeral curr_gain(em); + col_iterator it = M.col_begin(x_j), end = M.col_end(x_j); + for (; it != end; ++it) { + row r = it.get_row(); + var_t s = m_row2base[r.id()]; + var_info& vi = m_vars[s]; + numeral const& a_ij = it.get_row_entry().m_coeff; + numeral const& a_ii = vi.m_base_coeff; + bool inc_s = m.is_neg(a_ij) ? inc : !inc; + if ((inc_s && !vi.m_upper_valid) || (!inc_s && !vi.m_lower_valid)) { + continue; + } + // + // current gain: (value(x_i)-bound)*a_ii/a_ij + // + curr_gain = vi.m_value; + curr_gain -= inc_s?vi.m_upper:vi.m_lower; + em.mul(curr_gain, a_ii, curr_gain); + em.div(curr_gain, a_ij, curr_gain); + if (is_neg(curr_gain)) { + curr_gain.neg(); + } + if (x_i == null_var || (curr_gain < gain) || + (is_zero(gain) && is_zero(curr_gain) && s < x_i)) { + x_i = s; + gain = curr_gain; + new_a_ij = a_ij; + } + } + TRACE("simplex", tout << "x_i v" << x_i << "\n";); + return x_i; + } + + template + void simplex::check_blands_rule(var_t v, unsigned& num_repeated) { + if (m_bland) + return; + if (m_left_basis.contains(v)) { + num_repeated++; + if (num_repeated > m_blands_rule_threshold) { + TRACE("simplex", tout << "using blands rule, " << num_repeated << "\n";); + // std::cerr << "BLANDS RULE...\n"; + m_bland = true; + } + } + else { + m_left_basis.insert(v); + } + } + + template typename simplex::pivot_strategy_t simplex::pivot_strategy() { @@ -463,15 +695,31 @@ namespace simplex { SASSERT(M.well_formed()); for (unsigned i = 0; i < m_row2base.size(); ++i) { var_t s = m_row2base[i]; - SASSERT(i == m_vars[s].m_base2row); // unless row is deleted. + if (s == null_var) continue; + SASSERT(i == m_vars[s].m_base2row); // // TBD: extract coefficient of base variable and compare // with m_vars[s].m_base_coeff; // + // check that sum of assignments add up to 0 for every row. + row_iterator it = M.row_begin(row(i)), end = M.row_end(row(i)); + scoped_eps_numeral sum(em), tmp(em); + for (; it != end; ++it) { + em.mul(m_vars[it->m_var].m_value, it->m_coeff, tmp); + sum += tmp; + } + SASSERT(em.is_zero(sum)); } return true; } + template + bool simplex::is_feasible() const { + for (unsigned i = 0; i < m_vars.size(); ++i) { + if (below_lower(i) || above_upper(i)) return false; + } + return true; + } }; diff --git a/src/math/simplex/sparse_matrix.h b/src/math/simplex/sparse_matrix.h index cecd1c44a..6a85d69f0 100644 --- a/src/math/simplex/sparse_matrix.h +++ b/src/math/simplex/sparse_matrix.h @@ -226,7 +226,7 @@ namespace simplex { col_iterator col_begin(int v) const { return col_iterator(m_columns[v], m_rows, true); } col_iterator col_end(int v) const { return col_iterator(m_columns[v], m_rows, false); } - void display(std::ostream& out) const; + void display(std::ostream& out); bool well_formed() const; diff --git a/src/math/simplex/sparse_matrix_def.h b/src/math/simplex/sparse_matrix_def.h index 4e46c823d..20e9d676f 100644 --- a/src/math/simplex/sparse_matrix_def.h +++ b/src/math/simplex/sparse_matrix_def.h @@ -448,7 +448,7 @@ namespace simplex { */ template void sparse_matrix::del(row r) { - _row& rw = r.id(); + _row& rw = m_rows[r.id()]; for (unsigned i = 0; i < rw.m_entries.size(); ++i) { _row_entry& e = rw.m_entries[i]; if (!e.is_dead()) { @@ -561,12 +561,12 @@ namespace simplex { \brief display method */ template - void sparse_matrix::display(std::ostream& out) const { + void sparse_matrix::display(std::ostream& out) { for (unsigned i = 0; i < m_rows.size(); ++i) { if (m_rows[i].size() == 0) continue; row_iterator it = row_begin(row(i)), end = row_end(row(i)); for (; it != end; ++it) { - m.display(out, it->m_coeff) + m.display(out, it->m_coeff); out << "*v" << it->m_var << " "; } out << "\n"; diff --git a/src/opt/weighted_maxsat.cpp b/src/opt/weighted_maxsat.cpp index b5c53c3be..36047d434 100644 --- a/src/opt/weighted_maxsat.cpp +++ b/src/opt/weighted_maxsat.cpp @@ -867,6 +867,7 @@ namespace opt { m_upper += w; if (wmax < w) wmax = w; b = m.mk_fresh_const("b", m.mk_bool_sort()); + block.push_back(b); expr* bb = b; s.mc().insert(b->get_decl()); a = m.mk_fresh_const("a", m.mk_bool_sort()); diff --git a/src/smt/theory_arith_aux.h b/src/smt/theory_arith_aux.h index 77845fee6..e430b6515 100644 --- a/src/smt/theory_arith_aux.h +++ b/src/smt/theory_arith_aux.h @@ -1359,8 +1359,7 @@ namespace smt { SASSERT(valid_row_assignment()); SASSERT(satisfy_bounds()); result = skipped_row?BEST_EFFORT:OPTIMIZED; - break; - } + break; } if (x_i == null_theory_var) { // can increase/decrease x_j as much as we want. diff --git a/src/tactic/core/pb_preprocess_tactic.cpp b/src/tactic/core/pb_preprocess_tactic.cpp index 166cecc97..51e92ef11 100644 --- a/src/tactic/core/pb_preprocess_tactic.cpp +++ b/src/tactic/core/pb_preprocess_tactic.cpp @@ -244,6 +244,7 @@ public: expr* arg = args1[0].get(); bool neg = m.is_not(arg, arg); if (!is_uninterp_const(arg)) continue; + if (!m_vars.contains(to_app(arg))) continue; rec const& r = m_vars.find(to_app(arg)); unsigned_vector const& pos = neg?r.neg:r.pos; for (unsigned j = 0; j < pos.size(); ++j) { @@ -251,7 +252,8 @@ public: if (k == i) continue; if (!to_ge(g->form(k), args2, coeffs2, k2)) continue; if (subsumes(args1, coeffs1, k1, args2, coeffs2, k2)) { - g->update(k, m.mk_true()); + IF_VERBOSE(3, verbose_stream() << "replace " << mk_pp(g->form(k), m) << "\n";); + g->update(k, m.mk_true()); m_progress = true; } } @@ -518,6 +520,7 @@ private: expr_ref tmp1(m), tmp2(m); expr* fml1 = g->form(idx1); expr* fml2 = g->form(idx2); + TRACE("pb", tout << mk_pp(fml1, m) << " " << mk_pp(fml2, m) << "\n";); expr_ref_vector args1(m), args2(m); vector coeffs1, coeffs2; rational k1, k2; @@ -533,6 +536,7 @@ private: m.is_not(x, x); if (!is_app(x)) return; if (!m_vars.contains(to_app(x))) return; + TRACE("pb", tout << mk_pp(x, m) << "\n";); rec const& r = m_vars.find(to_app(x)); if (r.pos.size() != 1 || r.neg.size() != 1) return; if (r.pos[0] != idx2 && r.neg[0] != idx2) return; @@ -632,13 +636,17 @@ private: m_r.set_substitution(&sub); for (unsigned i = 0; i < positions.size(); ++i) { unsigned idx = positions[i]; - expr* f = g->form(idx); + expr_ref f(m); + f = g->form(idx); if (!m.is_true(f)) { m_r(f, tmp); - TRACE("pb", tout << mk_pp(f, m) << " -> " << tmp - << " by " << mk_pp(e, m) << " |-> " << mk_pp(v, m) << "\n";); - g->update(idx, tmp); // proof & dependencies. - m_progress = true; + if (tmp != f) { + TRACE("pb", tout << mk_pp(f, m) << " -> " << tmp + << " by " << mk_pp(e, m) << " |-> " << mk_pp(v, m) << "\n";); + IF_VERBOSE(3, verbose_stream() << "replace " << mk_pp(f, m) << " -> " << tmp << "\n";); + g->update(idx, tmp); // proof & dependencies. + m_progress = true; + } } } m_r.set_substitution(0); @@ -653,7 +661,7 @@ private: bool found = false; for (unsigned j = 0; !found && j < args2.size(); ++j) { if (args1[i] == args2[j]) { - if (coeffs1[1] > coeffs2[j]) return false; + if (coeffs1[i] > coeffs2[j]) return false; found = true; } } diff --git a/src/util/mpq_inf.h b/src/util/mpq_inf.h index eef066a6b..4b6a7d33d 100644 --- a/src/util/mpq_inf.h +++ b/src/util/mpq_inf.h @@ -239,6 +239,12 @@ public: m.neg(a.second); } + void abs(mpq_inf & a) { + if (is_neg(a)) { + neg(a); + } + } + void ceil(mpq_inf const & a, mpq & b) { if (m.is_int(a.first)) { // special cases for k - delta*epsilon where k is an integer diff --git a/src/util/scoped_numeral.h b/src/util/scoped_numeral.h index e03b56336..54d55f827 100644 --- a/src/util/scoped_numeral.h +++ b/src/util/scoped_numeral.h @@ -133,6 +133,12 @@ public: return a.m().is_nonpos(a); } + friend _scoped_numeral abs(_scoped_numeral const& a) { + _scoped_numeral res(a); + a.m().abs(res); + return res; + } + void neg() { m().neg(m_num); }