From eeaca949e029effdf9e1938132bb213013a3bed6 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Fri, 8 Jun 2018 15:05:08 -0700 Subject: [PATCH] cleanup in hnf.h Signed-off-by: Lev Nachmanson fixex in maximize_term Signed-off-by: Lev Nachmanson prepare for LIRA (#76) * merge Signed-off-by: Nikolaj Bjorner * simple mixed integer example Signed-off-by: Nikolaj Bjorner Nikolaj's fixes in theory_lra.cpp Signed-off-by: Lev Nachmanson fixes in maximize_term, disable find_cube for mixed case Signed-off-by: Lev Nachmanson fix cube heuristic for the mixed case Signed-off-by: Lev Nachmanson fix the build Signed-off-by: Lev Nachmanson fix in find cube delta's calculation Signed-off-by: Lev Nachmanson remove a printout Signed-off-by: Lev Nachmanson remove a blank line Signed-off-by: Lev Nachmanson test credentials Signed-off-by: Lev Nachmanson avoid double checks to add terms in hnf_cutter Signed-off-by: Lev Nachmanson fixes in add terms to hnf_cutter Signed-off-by: Lev Nachmanson remove a variable used for debug only Signed-off-by: Lev Nachmanson remove a field Signed-off-by: Lev Nachmanson remove a warning and hide m_A_orig under debug Signed-off-by: Lev Nachmanson use var_register to deal with mapping between external and local variables Signed-off-by: Lev Nachmanson tighten the loop in explain_implied_bound Signed-off-by: Lev Nachmanson fixes in theory_lra and relaxing debug checks in pop(), for the case when push() has been done from not fully initialized state Signed-off-by: Lev Nachmanson suppress hnf cutter when the numbers become too large Signed-off-by: Lev Nachmanson remove some debug code Signed-off-by: Lev Nachmanson adjusting hnf Signed-off-by: Lev Nachmanson --- src/smt/theory_lra.cpp | 128 ++++++++++++------ src/test/lp/lp.cpp | 61 ++++++++- src/util/lp/hnf.h | 36 ++--- src/util/lp/hnf_cutter.h | 71 +++++----- src/util/lp/int_solver.cpp | 57 ++++---- src/util/lp/int_solver.h | 7 +- src/util/lp/lar_solver.cpp | 185 ++++++++++++-------------- src/util/lp/lar_solver.h | 27 ++-- src/util/lp/lp_core_solver_base_def.h | 2 + src/util/lp/lp_primal_core_solver.h | 3 +- src/util/lp/lp_settings.h | 12 +- src/util/lp/numeric_pair.h | 4 +- src/util/lp/var_register.h | 91 ++++++++++--- 13 files changed, 424 insertions(+), 260 deletions(-) diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 7e2ee2a7e..512bf7787 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -319,7 +319,8 @@ class theory_lra::imp { void ensure_nra() { if (!m_nra) { m_nra = alloc(nra::solver, *m_solver.get(), m.limit(), ctx().get_params()); - for (unsigned i = 0; i < m_scopes.size(); ++i) { + for (auto const& _s : m_scopes) { + (void)_s; m_nra->push(); } } @@ -418,14 +419,63 @@ class theory_lra::imp { terms[index] = n1; } else if (is_app(n) && a.get_family_id() == to_app(n)->get_family_id()) { + bool is_first = !ctx().e_internalized(n); app* t = to_app(n); - found_not_handled(n); internalize_args(t); mk_enode(t); + theory_var v = mk_var(n); coeffs[vars.size()] = coeffs[index]; vars.push_back(v); ++index; + if (!is_first) { + // skip recursive internalization + } + else if (a.is_to_int(n, n1)) { + if (!ctx().relevancy()) + mk_to_int_axiom(t); + } + else if (a.is_to_real(n, n1)) { + theory_var v1 = mk_var(n1); + internalize_eq(v, v1); + } + else if (a.is_idiv(n, n1, n2)) { + if (!a.is_numeral(n2, r) || r.is_zero()) found_not_handled(n); + app * mod = a.mk_mod(n1, n2); + ctx().internalize(mod, false); + if (ctx().relevancy()) ctx().add_relevancy_dependency(n, mod); + } + else if (a.is_mod(n, n1, n2)) { + bool is_num = a.is_numeral(n2, r) && !r.is_zero(); + if (!is_num) { + found_not_handled(n); + } + else { + app_ref div(a.mk_idiv(n1, n2), m); + mk_enode(div); + theory_var w = mk_var(div); + theory_var u = mk_var(n1); + // add axioms: + // u = v + r*w + // abs(r) > v >= 0 + assert_idiv_mod_axioms(u, v, w, r); + } + if (!ctx().relevancy() && !is_num) mk_idiv_mod_axioms(n1, n2); + } + else if (a.is_rem(n, n1, n2)) { + if (!a.is_numeral(n2, r) || r.is_zero()) found_not_handled(n); + if (!ctx().relevancy()) mk_rem_axiom(n1, n2); + } + else if (a.is_div(n, n1, n2)) { + if (!a.is_numeral(n2, r) || r.is_zero()) found_not_handled(n); + if (!ctx().relevancy()) mk_div_axiom(n1, n2); + } + else if (a.is_power(n)) { + found_not_handled(n); + } + else { + found_not_handled(n); + } } else { if (is_app(n)) { @@ -629,6 +679,11 @@ class theory_lra::imp { ++m_stats.m_add_rows; TRACE("arith", m_solver->print_constraint(index, tout); tout << "\n";); } + + void add_def_constraint(lp::constraint_index index) { + m_constraint_sources.setx(index, definition_source, null_source); + ++m_stats.m_add_rows; + } void add_def_constraint(lp::constraint_index index, theory_var v) { m_constraint_sources.setx(index, definition_source, null_source); @@ -639,13 +694,12 @@ class theory_lra::imp { void internalize_eq(theory_var v1, theory_var v2) { enode* n1 = get_enode(v1); enode* n2 = get_enode(v2); - scoped_internalize_state st(*this); - st.vars().push_back(v1); - st.vars().push_back(v2); - st.coeffs().push_back(rational::one()); - st.coeffs().push_back(rational::minus_one()); - init_left_side(st); - add_eq_constraint(m_solver->add_constraint(m_left_side, lp::EQ, rational::zero()), n1, n2); + expr* o1 = n1->get_owner(); + expr* o2 = n2->get_owner(); + app_ref term(a.mk_sub(o1, o2), m); + theory_var z = internalize_def(term); + lp::var_index vi = get_var_index(z); + add_eq_constraint(m_solver->add_var_bound(vi, lp::EQ, rational::zero()), n1, n2); TRACE("arith", tout << "v" << v1 << " = " << "v" << v2 << ": " << mk_pp(n1->get_owner(), m) << " = " << mk_pp(n2->get_owner(), m) << "\n";); @@ -947,6 +1001,23 @@ public: mk_axiom(is_int, ~eq); } + // create axiom for + // u = v + r*w, + /// abs(r) > r >= 0 + void assert_idiv_mod_axioms(theory_var u, theory_var v, theory_var w, rational const& r) { + app_ref term(m); + term = a.mk_sub(get_enode(u)->get_owner(), + a.mk_add(get_enode(v)->get_owner(), + a.mk_mul(a.mk_numeral(r, true), + get_enode(w)->get_owner()))); + theory_var z = internalize_def(term); + lp::var_index vi = get_var_index(z); + add_def_constraint(m_solver->add_var_bound(vi, lp::GE, rational::zero())); + add_def_constraint(m_solver->add_var_bound(vi, lp::LE, rational::zero())); + add_def_constraint(m_solver->add_var_bound(get_var_index(v), lp::GE, rational::zero())); + add_def_constraint(m_solver->add_var_bound(get_var_index(v), lp::LT, abs(r))); + } + void mk_idiv_mod_axioms(expr * p, expr * q) { if (a.is_zero(q)) { return; @@ -966,6 +1037,7 @@ public: // q <= 0 or (p mod q) >= 0 // q <= 0 or (p mod q) < q // q >= 0 or (p mod q) < -q + // enable_trace("mk_bool_var"); mk_axiom(q_ge_0, eq); mk_axiom(q_le_0, eq); mk_axiom(q_ge_0, mod_ge_0); @@ -1001,26 +1073,22 @@ public: ctx().mk_th_axiom(get_id(), l1, l2); if (ctx().relevancy()) { ctx().mark_as_relevant(l1); - expr_ref e(m); - ctx().literal2expr(l2, e); - ctx().add_rel_watch(~l1, e); + ctx().add_rel_watch(~l1, ctx().bool_var2expr(l2.var())); } } void mk_axiom(literal l1, literal l2, literal l3) { ctx().mk_th_axiom(get_id(), l1, l2, l3); if (ctx().relevancy()) { - expr_ref e(m); ctx().mark_as_relevant(l1); - ctx().literal2expr(l2, e); - ctx().add_rel_watch(~l1, e); - ctx().literal2expr(l3, e); - ctx().add_rel_watch(~l2, e); + ctx().add_rel_watch(~l1, ctx().bool_var2expr(l2.var())); + ctx().add_rel_watch(~l2, ctx().bool_var2expr(l3.var())); } } literal mk_literal(expr* e) { expr_ref pinned(e, m); + TRACE("mk_bool_var", tout << pinned << " " << pinned->get_id() << "\n";); if (!ctx().e_internalized(e)) { ctx().internalize(e, false); } @@ -2606,33 +2674,15 @@ public: } theory_lra::inf_eps maximize(theory_var v, expr_ref& blocker, bool& has_shared) { - lp::var_index vi = m_theory_var2var_index.get(v, UINT_MAX); - vector > coeffs; - rational coeff(0); - // - // TBD: change API for maximize_term to take a proper term as input. - // - if (m_solver->is_term(vi)) { - const lp::lar_term& term = m_solver->get_term(vi); - for (auto & ti : term.m_coeffs) { - SASSERT(!m_solver->is_term(ti.first)); - coeffs.push_back(std::make_pair(ti.second, ti.first)); - } - coeff = term.m_v; - } - else { - coeffs.push_back(std::make_pair(rational::one(), vi)); - coeff = rational::zero(); - } lp::impq term_max; - lp::lp_status st = m_solver->maximize_term(coeffs, term_max); + lp::lp_status st = m_solver->maximize_term(v, term_max); + inf_rational val(term_max.x, term_max.y); if (st == lp::lp_status::OPTIMAL) { blocker = mk_gt(v); - inf_rational val(term_max.x + coeff, term_max.y); return inf_eps(rational::zero(), val); } if (st == lp::lp_status::FEASIBLE) { - lp_assert(false); // not implemented - // todo: what to return? + // todo , TODO , not sure what happens here + return inf_eps(rational::zero(), val); } else { SASSERT(st == lp::lp_status::UNBOUNDED); diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index e762398a8..76f9bd7f2 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -1943,6 +1943,7 @@ void setup_args_parser(argument_parser & parser) { parser.add_option_with_help_string("--test_mpq", "test rationals"); parser.add_option_with_help_string("--test_mpq_np", "test rationals"); parser.add_option_with_help_string("--test_mpq_np_plus", "test rationals using plus instead of +="); + parser.add_option_with_help_string("--maximize_term", "test maximize_term()"); } struct fff { int a; int b;}; @@ -2491,6 +2492,8 @@ void test_lar_solver(argument_parser & args_parser) { test_lar_on_file(fn, args_parser); return; } + + std::cout << "give option --file or --filelist to test_lar_solver\n"; } void test_numeric_pair() { @@ -3302,7 +3305,7 @@ void test_determinant() { M[2][0] = 0; M[2][1] = 1; M[2][2] = 4; svector r; std::cout << "M = "; M.print(std::cout, 4); endl(std::cout); - mpq d = hnf_calc::determinant_of_rectangular_matrix(M, r); + mpq d = hnf_calc::determinant_of_rectangular_matrix(M, r, mpq((int)100000)); std::cout << "det M = " << d << std::endl; std::cout << "rank = " << r.size() << std::endl; } @@ -3314,7 +3317,7 @@ void test_determinant() { M[3][0] = 6; M[3][1] = -2; M[3][2] = 2; M[3][3] = 2; M[3][4] = 6; M[3][5] = -2; svector r; std::cout << "M = "; M.print(std::cout, 4); endl(std::cout); - mpq d = hnf_calc::determinant_of_rectangular_matrix(M, r); + mpq d = hnf_calc::determinant_of_rectangular_matrix(M, r, mpq((int)1000000)); std::cout << "det M = " << d << std::endl; std::cout << "rank = " << r.size() << std::endl; } @@ -3333,7 +3336,7 @@ void fill_general_matrix(general_matrix & M) { void call_hnf(general_matrix& A) { svector r; - mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r); + mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r, mpq((int)1000000000)); A.shrink_to_rank(r); hnf h(A, d); } @@ -3509,7 +3512,50 @@ void test_larger_generated_hnf() { call_hnf(A); std::cout << "test_larger_generated_rank_hnf passed" << std::endl; } - +#endif +void test_maximize_term() { + std::cout << "test_maximize_term\n"; + lar_solver solver; + int_solver i_solver(&solver); // have to create it too + unsigned _x = 0; + unsigned _y = 1; + var_index x = solver.add_var(_x, false); + var_index y = solver.add_var(_y, true); + vector> term_ls; + term_ls.push_back(std::pair((int)1, x)); + term_ls.push_back(std::pair((int)-1, y)); + unsigned term_x_min_y = solver.add_term(term_ls, mpq(0)); + term_ls.clear(); + term_ls.push_back(std::pair((int)2, x)); + term_ls.push_back(std::pair((int)2, y)); + + unsigned term_2x_pl_2y = solver.add_term(term_ls, mpq(0)); + solver.add_var_bound(term_x_min_y, LE, zero_of_type()); + solver.add_var_bound(term_2x_pl_2y, LE, mpq((int)5)); + solver.find_feasible_solution(); + lp_assert(solver.get_status() == lp_status::OPTIMAL); + solver.print_constraints(std::cout); + std::unordered_map model; + solver.get_model(model); + for (auto p : model) { + std::cout<< "v[" << p.first << "] = " << p.second << std::endl; + } + std::cout << "calling int_solver\n"; + lar_term t; mpq k; explanation ex; bool upper; + lia_move lm = i_solver.check(t, k, ex, upper); + lp_assert(lm == lia_move::sat); + impq term_max; + lp_status st = solver.maximize_term(term_2x_pl_2y, term_max); + + std::cout << "status = " << lp_status_to_string(st) << std::endl; + std::cout << "term_max = " << term_max << std::endl; + solver.get_model(model); + for (auto p : model) { + std::cout<< "v[" << p.first << "] = " << p.second << std::endl; + } + +} +#ifdef Z3DEBUG void test_hnf() { test_larger_generated_hnf(); test_small_generated_hnf(); @@ -3662,7 +3708,12 @@ void test_lp_local(int argn, char**argv) { ret = 0; return finalize(ret); } - + + if (args_parser.option_is_used("--maximize_term")) { + test_maximize_term(); + ret = 0; + return finalize(ret); + } if (args_parser.option_is_used("--test_lp_0")) { test_lp_0(); diff --git a/src/util/lp/hnf.h b/src/util/lp/hnf.h index e0d7b0f20..0e14972da 100644 --- a/src/util/lp/hnf.h +++ b/src/util/lp/hnf.h @@ -126,33 +126,35 @@ bool prepare_pivot_for_lower_triangle(M &m, unsigned r) { } template -void pivot_column_non_fractional(M &m, unsigned r) { +void pivot_column_non_fractional(M &m, unsigned r, bool & overflow, const mpq & big_number) { lp_assert(!is_zero(m[r][r])); for (unsigned j = r + 1; j < m.column_count(); j++) { - for (unsigned i = r + 1; i < m.row_count(); i++) { + for (unsigned i = r + 1; i < m.row_count(); i++) { m[i][j] = (r > 0) ? (m[r][r]*m[i][j] - m[i][r]*m[r][j]) / m[r-1][r-1] : (m[r][r]*m[i][j] - m[i][r]*m[r][j]); + if (m[i][j] >= big_number) { + overflow = true; + return; + } + lp_assert(is_int(m[i][j])); } } - // debug - for (unsigned k = r + 1; k < m.row_count(); k++) { - m[k][r] = zero_of_type(); - } - } // returns the rank of the matrix template -unsigned to_lower_triangle_non_fractional(M &m) { +unsigned to_lower_triangle_non_fractional(M &m, bool & overflow, const mpq& big_number) { unsigned i = 0; for (; i < m.row_count(); i++) { if (!prepare_pivot_for_lower_triangle(m, i)) { return i; } - pivot_column_non_fractional(m, i); + pivot_column_non_fractional(m, i, overflow, big_number); + if (overflow) + return 0; } lp_assert(i == m.row_count()); return i; @@ -184,9 +186,12 @@ mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) { // The gcd of these minors is the return value. template -mpq determinant_of_rectangular_matrix(const M& m, svector & basis_rows) { +mpq determinant_of_rectangular_matrix(const M& m, svector & basis_rows, const mpq& big_number) { auto m_copy = m; - unsigned rank = to_lower_triangle_non_fractional(m_copy); + bool overflow = false; + unsigned rank = to_lower_triangle_non_fractional(m_copy, overflow, big_number); + if (overflow) + return big_number; if (rank == 0) return one_of_type(); @@ -216,9 +221,8 @@ class hnf { M m_H; M m_U; M m_U_reverse; -#endif - M m_A_orig; +#endif M m_W; vector m_buffer; unsigned m_m; @@ -529,14 +533,14 @@ private: // with some changes related to that we create a low triangle matrix // with non-positive elements under the diagonal void process_column_in_row_modulo() { - mpq& aii = m_W[m_i][m_i]; + const mpq& aii = m_W[m_i][m_i]; const mpq& aij = m_W[m_i][m_j]; mpq d, p,q; hnf_calc::extended_gcd_minimal_uv(aii, aij, d, p, q); if (is_zero(d)) return; - mpq aii_over_d = aii / d; - mpq aij_over_d = aij / d; + mpq aii_over_d = mod_R(aii / d); + mpq aij_over_d = mod_R(aij / d); buffer_p_col_i_plus_q_col_j_W_modulo(p, q); pivot_column_i_to_column_j_W_modulo(- aij_over_d, aii_over_d); copy_buffer_to_col_i_W_modulo(); diff --git a/src/util/lp/hnf_cutter.h b/src/util/lp/hnf_cutter.h index 1035183d7..90cdd5a6d 100644 --- a/src/util/lp/hnf_cutter.h +++ b/src/util/lp/hnf_cutter.h @@ -31,15 +31,18 @@ class hnf_cutter { vector m_terms; svector m_constraints_for_explanation; vector m_right_sides; - unsigned m_row_count; - unsigned m_column_count; lp_settings & m_settings; - + mpq m_abs_max; + bool m_overflow; public: - hnf_cutter(lp_settings & settings) : m_row_count(0), m_column_count(0), m_settings(settings) {} - unsigned row_count() const { - return m_row_count; + const mpq & abs_max() const { return m_abs_max; } + + hnf_cutter(lp_settings & settings) : m_settings(settings), + m_abs_max(zero_of_type()) {} + + unsigned terms_count() const { + return m_terms.size(); } const vector& terms() const { return m_terms; } @@ -52,8 +55,9 @@ public: m_var_register.clear(); m_terms.clear(); m_constraints_for_explanation.clear(); - m_right_sides.clear(); - m_row_count = m_column_count = 0; + m_right_sides.clear(); + m_abs_max = zero_of_type(); + m_overflow = false; } void add_term(const lar_term* t, const mpq &rs, constraint_index ci) { m_terms.push_back(t); @@ -61,12 +65,12 @@ public: m_constraints_for_explanation.push_back(ci); for (const auto &p : *t) { m_var_register.add_var(p.var()); - } - if (m_terms.size() <= m_var_register.size()) { // capture the maximal acceptable sizes - m_row_count = m_terms.size(); - m_column_count = m_var_register.size(); + mpq t = abs(ceil(p.coeff())); + if (t > m_abs_max) + m_abs_max = t; } } + void print(std::ostream & out) { out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl; } @@ -76,10 +80,8 @@ public: } void init_matrix_A() { - m_A = general_matrix(m_row_count, m_column_count); - // use the last suitable counts to make the number - // of rows less than or equal to the number of columns - for (unsigned i = 0; i < m_row_count; i++) + m_A = general_matrix(terms_count(), vars().size()); + for (unsigned i = 0; i < terms_count(); i++) initialize_row(i); } @@ -123,7 +125,6 @@ public: return ret; } - // fills e_i*H_minus_1 void get_ei_H_minus_1(unsigned i, const general_matrix& H, vector & row) { // we solve x = ei * H_min_1 @@ -141,29 +142,28 @@ public: row[k] = -t / H[k][k]; } - // test region - vector ei(H.row_count(), zero_of_type()); - ei[i] = one_of_type(); - vector pr = row * H; - pr.shrink(ei.size()); - lp_assert(ei == pr); - // end test region + // // test region + // vector ei(H.row_count(), zero_of_type()); + // ei[i] = one_of_type(); + // vector pr = row * H; + // pr.shrink(ei.size()); + // lp_assert(ei == pr); + // // end test region } void fill_term(const vector & row, lar_term& t) { for (unsigned j = 0; j < row.size(); j++) { if (!is_zero(row[j])) - t.add_monomial(row[j], m_var_register.local_var_to_user_var(j)); + t.add_monomial(row[j], m_var_register.local_to_external(j)); } } #ifdef Z3DEBUG vector transform_to_local_columns(const vector & x) const { vector ret; - lp_assert(m_column_count <= m_var_register.size()); - for (unsigned j = 0; j < m_column_count;j++) { - lp_assert(is_zero(x[m_var_register.local_var_to_user_var(j)].y)); - ret.push_back(x[m_var_register.local_var_to_user_var(j)].x); + for (unsigned j = 0; j < vars().size(); j++) { + lp_assert(is_zero(x[m_var_register.local_to_external(j)].y)); + ret.push_back(x[m_var_register.local_to_external(j)].x); } return ret; } @@ -176,6 +176,8 @@ public: m_constraints_for_explanation = new_expl; } + bool overflow() const { return m_overflow; } + lia_move create_cut(lar_term& t, mpq& k, explanation& ex, bool & upper #ifdef Z3DEBUG , @@ -185,7 +187,16 @@ public: // we suppose that x0 has at least one non integer element init_matrix_A(); svector basis_rows; - mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows); + mpq big_number = m_abs_max.expt(3); + mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows, big_number); + + // std::cout << "max = " << m_abs_max << ", d = " << d << ", d/max = " << ceil (d /m_abs_max) << std::endl; + //std::cout << "max cube " << m_abs_max * m_abs_max * m_abs_max << std::endl; + + if (d >= big_number) { + return lia_move::undef; + } + if (m_settings.get_cancel_flag()) return lia_move::undef; if (basis_rows.size() < m_A.row_count()) { diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index 020b24171..5033f2852 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -416,6 +416,8 @@ impq int_solver::get_cube_delta_for_term(const lar_term& t) const { bool seen_minus = false; bool seen_plus = false; for(const auto & p : t) { + if (!is_int(p.var())) + goto usual_delta; const mpq & c = p.coeff(); if (c == one_of_type()) { seen_plus = true; @@ -431,9 +433,10 @@ impq int_solver::get_cube_delta_for_term(const lar_term& t) const { } usual_delta: mpq delta = zero_of_type(); - for (const auto & p : t) { - delta += abs(p.coeff()); - } + for (const auto & p : t) + if (is_int(p.var())) + delta += abs(p.coeff()); + delta *= mpq(1, 2); return impq(delta); } @@ -443,7 +446,6 @@ bool int_solver::tighten_term_for_cube(unsigned i) { if (!m_lar_solver->term_is_used_as_row(ti)) return true; const lar_term* t = m_lar_solver->terms()[i]; - impq delta = get_cube_delta_for_term(*t); TRACE("cube", m_lar_solver->print_term_as_indices(*t, tout); tout << ", delta = " << delta;); if (is_zero(delta)) @@ -533,20 +535,28 @@ lia_move int_solver::gomory_cut() { void int_solver::try_add_term_to_A_for_hnf(unsigned i) { mpq rs; const lar_term* t = m_lar_solver->terms()[i]; -#if Z3DEBUG - for (const auto & p : *t) { - if (!is_int(p.var())) { - lp_assert(false); - } - } -#endif constraint_index ci; - if (m_lar_solver->get_equality_and_right_side_for_term_on_current_x(i, rs, ci)) { + if (!hnf_cutter_is_full() && m_lar_solver->get_equality_and_right_side_for_term_on_current_x(i, rs, ci)) { m_hnf_cutter.add_term(t, rs, ci); } } -bool int_solver::hnf_has_non_integral_var() const { +bool int_solver::hnf_cutter_is_full() const { + return + m_hnf_cutter.terms_count() >= settings().limit_on_rows_for_hnf_cutter + || + m_hnf_cutter.vars().size() >= settings().limit_on_columns_for_hnf_cutter; +} + +lp_settings& int_solver::settings() { + return m_lar_solver->settings(); +} + +const lp_settings& int_solver::settings() const { + return m_lar_solver->settings(); +} + +bool int_solver::hnf_has_var_with_non_integral_value() const { for (unsigned j : m_hnf_cutter.vars()) if (get_value(j).is_int() == false) return true; @@ -555,10 +565,10 @@ bool int_solver::hnf_has_non_integral_var() const { bool int_solver::init_terms_for_hnf_cut() { m_hnf_cutter.clear(); - for (unsigned i = 0; i < m_lar_solver->terms().size() && m_hnf_cutter.row_count() < settings().limit_on_rows_for_hnf_cutter; i++) { + for (unsigned i = 0; i < m_lar_solver->terms().size() && !hnf_cutter_is_full(); i++) { try_add_term_to_A_for_hnf(i); } - return hnf_has_non_integral_var(); + return hnf_has_var_with_non_integral_value(); } lia_move int_solver::make_hnf_cut() { @@ -583,7 +593,7 @@ lia_move int_solver::make_hnf_cut() { for (unsigned i : m_hnf_cutter.constraints_for_explanation()) { m_ex->push_justification(i); } - } + } return r; } @@ -609,7 +619,6 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) { r = patch_nbasic_columns(); if (r != lia_move::undef) return r; ++m_branch_cut_counter; - r = find_cube(); if (r != lia_move::undef) return r; @@ -697,14 +706,12 @@ void int_solver::set_value_for_nbasic_column(unsigned j, const impq & new_val) { m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta); } -void int_solver::patch_nbasic_column(unsigned j) { +void int_solver::patch_nbasic_column(unsigned j, bool patch_only_int_vals) { auto & lcs = m_lar_solver->m_mpq_lar_core_solver; impq & val = lcs.m_r_x[j]; bool val_is_int = val.is_int(); - if (settings().m_int_patch_only_integer_values) { - if (!val_is_int) + if (patch_only_int_vals && !val_is_int) return; - } bool inf_l, inf_u; impq l, u; @@ -754,7 +761,7 @@ lia_move int_solver::patch_nbasic_columns() { settings().st().m_patches++; lp_assert(is_feasible()); for (unsigned j : m_lar_solver->m_mpq_lar_core_solver.m_r_nbasis) { - patch_nbasic_column(j); + patch_nbasic_column(j, settings().m_int_patch_only_integer_values); } lp_assert(is_feasible()); if (!has_inf_int()) { @@ -1128,12 +1135,6 @@ bool int_solver::at_upper(unsigned j) const { } } - - -lp_settings& int_solver::settings() { - return m_lar_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; for (auto &c: rslv.m_A.m_rows[row_index]) { diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index 0bc7d5ef0..13db7470e 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -50,6 +50,7 @@ public: // 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, bool & upper); + 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); bool is_base(unsigned j) const; @@ -100,6 +101,7 @@ private: int find_inf_int_boxed_base_column_with_smallest_range(unsigned&); int get_kth_inf_int(unsigned) const; lp_settings& settings(); + const lp_settings& settings() const; bool move_non_basic_columns_to_bounds(); void branch_infeasible_int_var(unsigned); lia_move mk_gomory_cut(unsigned inf_col, const row_strip& row); @@ -157,7 +159,8 @@ public: bool init_terms_for_hnf_cut(); bool hnf_matrix_is_empty() const; void try_add_term_to_A_for_hnf(unsigned term_index); - bool hnf_has_non_integral_var() const; - void patch_nbasic_column(unsigned j); + bool hnf_has_var_with_non_integral_value() const; + bool hnf_cutter_is_full() const; + void patch_nbasic_column(unsigned j, bool patch_only_int_vals); }; } diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index def7edb15..1fd124ab2 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -31,8 +31,7 @@ lar_solver::lar_solver() : m_status(lp_status::OPTIMAL), m_infeasible_column_index(-1), m_terms_start_index(1000000), m_mpq_lar_core_solver(m_settings, *this), - m_int_solver(nullptr), - m_has_int_var(false) + m_int_solver(nullptr) {} void lar_solver::set_track_pivoted_rows(bool v) { @@ -199,7 +198,7 @@ void lar_solver::calculate_implied_bounds_for_row(unsigned i, bound_propagator & } 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]; + unsigned ext_var_or_term = m_var_register.local_to_external(j); return ext_var_or_term < m_terms_start_index ? j : ext_var_or_term; } @@ -212,21 +211,14 @@ void lar_solver::explain_implied_bound(implied_bound & ib, bound_propagator & bp unsigned i = ib.m_row_or_term_index; int bound_sign = ib.m_is_lower_bound? 1: -1; int j_sign = (ib.m_coeff_before_j_is_pos ? 1 :-1) * bound_sign; - unsigned m_j = ib.m_j; - if (is_term(m_j)) { - auto it = m_ext_vars_to_columns.find(m_j); - lp_assert(it != m_ext_vars_to_columns.end()); - m_j = it->second.internal_j(); + unsigned bound_j = ib.m_j; + if (is_term(bound_j)) { + bound_j = m_var_register.external_to_local(bound_j); } for (auto const& r : A_r().m_rows[i]) { unsigned j = r.m_j; - mpq const& a = r.get_val(); if (j == m_j) continue; - if (is_term(j)) { - auto it = m_ext_vars_to_columns.find(j); - lp_assert(it != m_ext_vars_to_columns.end()); - j = it->second.internal_j(); - } + mpq const& a = r.get_val(); int a_sign = is_pos(a)? 1: -1; int sign = j_sign * a_sign; const ul_pair & ul = m_columns_to_ul_pairs[j]; @@ -240,7 +232,7 @@ void lar_solver::explain_implied_bound(implied_bound & ib, bound_propagator & bp bool lar_solver::term_is_used_as_row(unsigned term) const { lp_assert(is_term(term)); - return contains(m_ext_vars_to_columns, term); + return m_var_register.external_is_used(term); } void lar_solver::propagate_bounds_on_terms(bound_propagator & bp) { @@ -273,10 +265,6 @@ lp_status lar_solver::get_status() const { return m_status;} void lar_solver::set_status(lp_status s) {m_status = s;} -bool lar_solver::has_int_var() const { - return m_has_int_var; -} - lp_status lar_solver::find_feasible_solution() { m_settings.st().m_make_feasible++; if (A_r().column_count() > m_settings.st().m_max_cols) @@ -352,12 +340,9 @@ void lar_solver::pop(unsigned k) { TRACE("arith_int", tout << "pop" << std::endl;); TRACE("lar_solver", tout << "k = " << k << std::endl;); - int n_was = static_cast(m_ext_vars_to_columns.size()); m_infeasible_column_index.pop(k); unsigned n = m_columns_to_ul_pairs.peek_size(k); - for (unsigned j = n_was; j-- > n;) - m_ext_vars_to_columns.erase(m_columns_to_ext_vars_or_term_indices[j]); - m_columns_to_ext_vars_or_term_indices.resize(n); + m_var_register.shrink(n); TRACE("arith_int", tout << "pop" << std::endl;); if (m_settings.use_tableau()) { pop_tableau(); @@ -375,7 +360,6 @@ void lar_solver::pop(unsigned k) { lp_assert(ax_is_correct()); - lp_assert(m_mpq_lar_core_solver.m_r_solver.inf_set_is_correct()); m_constraint_count.pop(k); for (unsigned i = m_constraint_count; i < m_constraints.size(); i++) delete m_constraints[i]; @@ -404,18 +388,17 @@ vector lar_solver::get_all_constraint_indices() const { return ret; } -bool lar_solver::maximize_term_on_tableau(const vector> & term, +bool lar_solver::maximize_term_on_tableau(const lar_term & term, impq &term_max) { if (settings().simplex_strategy() == simplex_strategy_enum::undecided) decide_on_strategy_and_adjust_initial_state(); - + + m_mpq_lar_core_solver.m_r_solver.set_status(lp_status::FEASIBLE); m_mpq_lar_core_solver.solve(); if (m_mpq_lar_core_solver.m_r_solver.get_status() == lp_status::UNBOUNDED) return false; - term_max = 0; - for (auto & p : term) - term_max += p.first * m_mpq_lar_core_solver.m_r_x[p.second]; + term_max = term.apply(m_mpq_lar_core_solver.m_r_x); return true; } @@ -433,13 +416,13 @@ bool lar_solver::reduced_costs_are_zeroes_for_r_solver() const { return true; } -void lar_solver::set_costs_to_zero(const vector> & term) { +void lar_solver::set_costs_to_zero(const lar_term& term) { auto & rslv = m_mpq_lar_core_solver.m_r_solver; auto & jset = m_mpq_lar_core_solver.m_r_solver.m_inf_set; // hijack this set that should be empty right now lp_assert(jset.m_index.size()==0); - for (auto & p : term) { - unsigned j = p.second; + for (const auto & p : term) { + unsigned j = p.var(); rslv.m_costs[j] = zero_of_type(); int i = rslv.m_basis_heading[j]; if (i < 0) @@ -459,25 +442,25 @@ void lar_solver::set_costs_to_zero(const vector> & ter lp_assert(costs_are_zeros_for_r_solver()); } -void lar_solver::prepare_costs_for_r_solver(const vector> & term) { +void lar_solver::prepare_costs_for_r_solver(const lar_term & term) { auto & rslv = m_mpq_lar_core_solver.m_r_solver; rslv.m_using_infeas_costs = false; lp_assert(costs_are_zeros_for_r_solver()); lp_assert(reduced_costs_are_zeroes_for_r_solver()); rslv.m_costs.resize(A_r().column_count(), zero_of_type()); - for (auto & p : term) { - unsigned j = p.second; - rslv.m_costs[j] = p.first; + for (const auto & p : term) { + unsigned j = p.var(); + rslv.m_costs[j] = p.coeff(); if (rslv.m_basis_heading[j] < 0) - rslv.m_d[j] += p.first; + rslv.m_d[j] += p.coeff(); else - rslv.update_reduced_cost_for_basic_column_cost_change(- p.first, j); + rslv.update_reduced_cost_for_basic_column_cost_change(- p.coeff(), j); } lp_assert(rslv.reduced_costs_are_correct_tableau()); } -bool lar_solver::maximize_term_on_corrected_r_solver(const vector> & term, +bool lar_solver::maximize_term_on_corrected_r_solver(lar_term & term, impq &term_max) { settings().backup_costs = false; switch (settings().simplex_strategy()) { @@ -514,16 +497,35 @@ bool lar_solver::remove_from_basis(unsigned j) { return m_mpq_lar_core_solver.m_r_solver.remove_from_basis(j); } +lar_term lar_solver::get_term_to_maximize(unsigned ext_j) const { + unsigned local_j; + if (m_var_register.external_is_used(ext_j, local_j)) { + lar_term r; + r. add_monomial(one_of_type(), local_j); + return r; + } + + return get_term(ext_j); +} -// starting from a given feasible state look for the maximum of the term -// return true if found and false if unbounded -lp_status lar_solver::maximize_term(const vector> & term, +lp_status lar_solver::maximize_term(unsigned ext_j, impq &term_max) { - lp_assert(m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()); + bool was_feasible = m_mpq_lar_core_solver.m_r_solver.calc_current_x_is_feasible_include_non_basis(); + impq prev_value; + lar_term term = get_term_to_maximize(ext_j); + auto backup = m_mpq_lar_core_solver.m_r_x; + if (was_feasible) { + prev_value = term.apply(m_mpq_lar_core_solver.m_r_x); + } + m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = false; - if (!maximize_term_on_corrected_r_solver(term, term_max)) + if (!maximize_term_on_corrected_r_solver(term, term_max)) { + m_mpq_lar_core_solver.m_r_x = backup; return lp_status::UNBOUNDED; + } + impq opt_val = term_max; + bool change = false; for (unsigned j = 0; j < m_mpq_lar_core_solver.m_r_x.size(); j++) { if (!column_is_int(j)) @@ -534,18 +536,19 @@ lp_status lar_solver::maximize_term(const vector> & te if (!remove_from_basis(j)) // consider a special version of remove_from_basis that would not remove inf_int columns return lp_status::FEASIBLE; // it should not happen } - m_int_solver->patch_nbasic_column(j); + m_int_solver->patch_nbasic_column(j, false); if (!column_value_is_integer(j)) return lp_status::FEASIBLE; change = true; } if (change) { - term_max = zero_of_type(); - for (const auto& p : term) { - term_max += p.first * m_mpq_lar_core_solver.m_r_x[p.second]; - } + term_max = term.apply(m_mpq_lar_core_solver.m_r_x); } - return lp_status::OPTIMAL; + if (was_feasible && term_max < prev_value) { + term_max = prev_value; + m_mpq_lar_core_solver.m_r_x = backup; + } + return term_max == opt_val? lp_status::OPTIMAL :lp_status::FEASIBLE; } @@ -907,10 +910,10 @@ column_type lar_solver::get_column_type(const column_info & ci) { std::string lar_solver::get_column_name(unsigned j) const { if (j >= m_terms_start_index) return std::string("_t") + T_to_string(j); - if (j >= m_columns_to_ext_vars_or_term_indices.size()) + if (j >= m_var_register.size()) return std::string("_s") + T_to_string(j); - return std::string("v") + T_to_string(m_columns_to_ext_vars_or_term_indices[j]); + return std::string("v") + T_to_string(m_var_register.local_to_external(j)); } bool lar_solver::all_constrained_variables_are_registered(const vector>& left_side) { @@ -1392,7 +1395,8 @@ void lar_solver::pop_tableau() { lp_assert(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct()); // We remove last variables starting from m_column_names.size() to m_vec_of_canonic_left_sides.size(). // At this moment m_column_names is already popped - while (A_r().column_count() > m_columns_to_ext_vars_or_term_indices.size()) + unsigned size = m_var_register.size(); + while (A_r().column_count() > size) remove_last_column_from_tableau(); lp_assert(m_mpq_lar_core_solver.m_r_solver.m_costs.size() == A_r().column_count()); lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis.size() == A_r().row_count()); @@ -1469,9 +1473,7 @@ bool lar_solver::var_is_int(var_index v) const { } bool lar_solver::column_is_int(unsigned j) const { - unsigned ext_var = m_columns_to_ext_vars_or_term_indices[j]; - lp_assert(contains(m_ext_vars_to_columns, ext_var)); - return m_ext_vars_to_columns.find(ext_var)->second.is_integer(); + return m_var_register.local_is_int(j); } bool lar_solver::column_is_fixed(unsigned j) const { @@ -1480,9 +1482,7 @@ bool lar_solver::column_is_fixed(unsigned j) const { bool lar_solver::ext_var_is_int(var_index ext_var) const { - auto it = m_ext_vars_to_columns.find(ext_var); - lp_assert(it != m_ext_vars_to_columns.end()); - return it->second.is_integer(); + return m_var_register.external_is_int(ext_var); } // below is the initialization functionality of lar_solver @@ -1492,30 +1492,22 @@ bool lar_solver::strategy_is_undecided() const { } var_index lar_solver::add_var(unsigned ext_j, bool is_int) { - if (is_int) - m_has_int_var = true; - TRACE("add_var", tout << "adding var " << ext_j << (is_int? " int" : " nonint") << std::endl;); - var_index i; + var_index local_j; lp_assert(ext_j < m_terms_start_index); - auto it = m_ext_vars_to_columns.find(ext_j); - if (it != m_ext_vars_to_columns.end()) { - return it->second.internal_j(); - } + if (m_var_register.external_is_used(ext_j, local_j)) + return local_j; lp_assert(m_columns_to_ul_pairs.size() == A_r().column_count()); - i = A_r().column_count(); + local_j = A_r().column_count(); m_columns_to_ul_pairs.push_back(ul_pair(static_cast(-1))); add_non_basic_var_to_core_fields(ext_j, is_int); lp_assert(sizes_are_correct()); - return i; + return local_j; } void lar_solver::register_new_ext_var_index(unsigned ext_v, bool is_int) { - lp_assert(!contains(m_ext_vars_to_columns, ext_v)); - unsigned j = static_cast(m_ext_vars_to_columns.size()); - m_ext_vars_to_columns.insert(std::make_pair(ext_v, ext_var_info(j, is_int))); - lp_assert(m_columns_to_ext_vars_or_term_indices.size() == j); - m_columns_to_ext_vars_or_term_indices.push_back(ext_v); + lp_assert(!m_var_register.external_is_used(ext_v)); + m_var_register.add_var(ext_v, is_int); } void lar_solver::add_non_basic_var_to_core_fields(unsigned ext_j, bool is_int) { @@ -1631,7 +1623,7 @@ var_index lar_solver::add_term(const vector> & coeffs, if (m_settings.bound_propagation()) m_rows_with_changed_bounds.insert(A_r().row_count() - 1); } - lp_assert(m_ext_vars_to_columns.size() == A_r().column_count()); + lp_assert(m_var_register.size() == A_r().column_count()); return ret; } @@ -1716,9 +1708,8 @@ void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_k lp_assert(is_term(j)); unsigned adjusted_term_index = adjust_term_index(j); // lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int()); - auto it = m_ext_vars_to_columns.find(j); - if (it != m_ext_vars_to_columns.end()) { - unsigned term_j = it->second.internal_j(); + unsigned term_j; + if (m_var_register.external_is_used(j, term_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)); update_column_type_and_bound(term_j, kind, rs, ci); @@ -1804,10 +1795,10 @@ void lar_solver::adjust_initial_state_for_lu() { } void lar_solver::adjust_initial_state_for_tableau_rows() { - for (unsigned j = 0; j < m_terms.size(); j++) { - if (contains(m_ext_vars_to_columns, j + m_terms_start_index)) + for (unsigned i = 0; i < m_terms.size(); i++) { + if (m_var_register.external_is_used(i + m_terms_start_index)) continue; - add_row_from_term_no_constraint(m_terms[j], j + m_terms_start_index); + add_row_from_term_no_constraint(m_terms[i], i + m_terms_start_index); } } @@ -2122,21 +2113,18 @@ void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kin } bool lar_solver::column_corresponds_to_term(unsigned j) const { - return m_columns_to_ext_vars_or_term_indices[j] >= m_terms_start_index; + return m_var_register.local_to_external(j) >= m_terms_start_index; } -var_index lar_solver:: to_var_index(unsigned ext_j) const { - auto it = m_ext_vars_to_columns.find(ext_j); - lp_assert(it != m_ext_vars_to_columns.end()); - return it->second.internal_j(); +var_index lar_solver::to_column(unsigned ext_j) const { + return m_var_register.external_to_local(ext_j); } bool lar_solver::tighten_term_bounds_by_delta(unsigned term_index, const impq& delta) { unsigned tj = term_index + m_terms_start_index; - auto it = m_ext_vars_to_columns.find(tj); - if (it == m_ext_vars_to_columns.end()) - return true; - unsigned j = it->second.internal_j(); + unsigned j; + if (m_var_register.external_is_used(tj, j) == false) + return true; // the term is not a column so it has no bounds auto & slv = m_mpq_lar_core_solver.m_r_solver; TRACE("cube", tout << "delta = " << delta << std::endl; m_int_solver->display_column(tout, j); ); @@ -2166,7 +2154,7 @@ void lar_solver::update_delta_for_terms(const impq & delta, unsigned j, const ve for (unsigned i : terms_of_var) { lar_term & t = *m_terms[i]; auto it = t.m_coeffs.find(j); - unsigned tj = to_var_index(i + m_terms_start_index); + unsigned tj = to_column(i + m_terms_start_index); TRACE("cube", tout << "t.apply = " << t.apply(m_mpq_lar_core_solver.m_r_x) << ", m_mpq_lar_core_solver.m_r_x[tj]= " << m_mpq_lar_core_solver.m_r_x[tj];); TRACE("cube", print_term_as_indices(t, tout); @@ -2199,6 +2187,7 @@ void lar_solver::round_to_integer_solution() { fill_vars_to_terms(vars_to_terms); for (unsigned j = 0; j < column_count(); j++) { + if (column_is_int(j)) continue; if (column_corresponds_to_term(j)) continue; TRACE("cube", m_int_solver->display_column(tout, j);); impq& v = m_mpq_lar_core_solver.m_r_x[j]; @@ -2219,16 +2208,17 @@ void lar_solver::round_to_integer_solution() { bool lar_solver::get_equality_and_right_side_for_term_on_current_x(unsigned term_index, mpq & rs, constraint_index& ci) const { unsigned tj = term_index + m_terms_start_index; - auto it = m_ext_vars_to_columns.find(tj); - if (it == m_ext_vars_to_columns.end()) + unsigned j; + bool is_int; + if (m_var_register.external_is_used(tj, j, is_int) == false) + return false; // the term does not have bound because it does not correspond to a column + if (!is_int) // todo - allow for the next version of hnf return false; - unsigned j = it->second.internal_j(); impq term_val; bool term_val_ready = false; mpq b; bool is_strict; - if (has_upper_bound(j, ci, b, is_strict)) { - lp_assert(!is_strict); + if (has_upper_bound(j, ci, b, is_strict) && !is_strict) { lp_assert(b.is_int()); term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x); term_val_ready = true; @@ -2237,8 +2227,7 @@ bool lar_solver::get_equality_and_right_side_for_term_on_current_x(unsigned term return true; } } - if (has_lower_bound(j, ci, b, is_strict)) { - lp_assert(!is_strict); + if (has_lower_bound(j, ci, b, is_strict) && !is_strict) { if (!term_val_ready) term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x); lp_assert(b.is_int()); diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index ce975738f..d7f0fb076 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -47,15 +47,6 @@ namespace lp { class lar_solver : public column_namer { - class ext_var_info { - unsigned m_internal_j; // the internal index - bool m_is_integer; - public: - ext_var_info(unsigned j, var_index internal_j): ext_var_info(j, false) {} - ext_var_info(unsigned j , bool is_int) : m_internal_j(j), m_is_integer(is_int) {} - 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 @@ -100,8 +91,7 @@ class lar_solver : public column_namer { 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; + var_register m_var_register; stacked_vector m_columns_to_ul_pairs; vector m_constraints; private: @@ -119,8 +109,6 @@ public: lar_core_solver m_mpq_lar_core_solver; private: int_solver * m_int_solver; - bool m_has_int_var; - public : unsigned terms_start_index() const { return m_terms_start_index; } @@ -311,21 +299,21 @@ public: vector get_all_constraint_indices() const; - bool maximize_term_on_tableau(const vector> & term, + bool maximize_term_on_tableau(const lar_term & term, impq &term_max); bool costs_are_zeros_for_r_solver() const; bool reduced_costs_are_zeroes_for_r_solver() const; - void set_costs_to_zero(const vector> & term); + void set_costs_to_zero(const lar_term & term); - void prepare_costs_for_r_solver(const vector> & term); + void prepare_costs_for_r_solver(const lar_term & term); - bool maximize_term_on_corrected_r_solver(const vector> & term, + bool maximize_term_on_corrected_r_solver(lar_term & term, impq &term_max); // starting from a given feasible state look for the maximum of the term // return true if found and false if unbounded - lp_status maximize_term(const vector> & term, + lp_status maximize_term(unsigned ext_j , impq &term_max); @@ -578,7 +566,7 @@ public: lar_core_solver & get_core_solver() { return m_mpq_lar_core_solver; } bool column_corresponds_to_term(unsigned) const; void catch_up_in_updating_int_solver(); - var_index to_var_index(unsigned ext_j) const; + var_index to_column(unsigned ext_j) const; bool tighten_term_bounds_by_delta(unsigned, const impq&); void round_to_integer_solution(); void update_delta_for_terms(const impq & delta, unsigned j, const vector&); @@ -588,5 +576,6 @@ public: const vector & r_nbasis() const { return m_mpq_lar_core_solver.r_nbasis(); } bool get_equality_and_right_side_for_term_on_current_x(unsigned i, mpq &rs, constraint_index& ci) const; bool remove_from_basis(unsigned); + lar_term get_term_to_maximize(unsigned ext_j) const; }; } diff --git a/src/util/lp/lp_core_solver_base_def.h b/src/util/lp/lp_core_solver_base_def.h index dd3de910f..f20eaf042 100644 --- a/src/util/lp/lp_core_solver_base_def.h +++ b/src/util/lp/lp_core_solver_base_def.h @@ -973,6 +973,8 @@ template bool lp_core_solver_base::remove_from_ba indexed_vector w(m_basis.size()); // the buffer unsigned i = m_basis_heading[basic_j]; for (auto &c : m_A.m_rows[i]) { + if (c.var() == basic_j) + continue; if (pivot_column_general(c.var(), basic_j, w)) return true; } diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index 424823a7c..26fc550b3 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -504,8 +504,7 @@ public: } X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent; advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta ); - X xleaving = this->m_x[leaving]; - lp_assert(xleaving == new_val_for_leaving); + lp_assert(this->m_x[leaving] == new_val_for_leaving); if (this->current_x_is_feasible()) this->set_status(lp_status::OPTIMAL); } diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h index 553e58c62..0ef385cab 100644 --- a/src/util/lp/lp_settings.h +++ b/src/util/lp/lp_settings.h @@ -200,6 +200,7 @@ public: bool m_int_pivot_fixed_vars_from_basis; bool m_int_patch_only_integer_values; unsigned limit_on_rows_for_hnf_cutter; + unsigned limit_on_columns_for_hnf_cutter; unsigned random_next() { return m_rand(); } void set_random_seed(unsigned s) { m_rand.set_seed(s); } @@ -219,10 +220,10 @@ public: reps_in_scaler(20), pivot_epsilon(0.00000001), positive_price_epsilon(1e-7), - entering_diag_epsilon ( 1e-8), - c_partial_pivoting ( 10), // this is the constant c from page 410 - depth_of_rook_search ( 4), - using_partial_pivoting ( true), + entering_diag_epsilon (1e-8), + c_partial_pivoting (10), // this is the constant c from page 410 + depth_of_rook_search (4), + using_partial_pivoting (true), // dissertation of Achim Koberstein // if Bx - b is different at any component more that refactor_epsilon then we refactor refactor_tolerance ( 1e-4), @@ -264,7 +265,8 @@ public: m_chase_cut_solver_cycle_on_var(10), m_int_pivot_fixed_vars_from_basis(false), m_int_patch_only_integer_values(true), - limit_on_rows_for_hnf_cutter(75) + limit_on_rows_for_hnf_cutter(75), + limit_on_columns_for_hnf_cutter(150) {} void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; } diff --git a/src/util/lp/numeric_pair.h b/src/util/lp/numeric_pair.h index baff91dac..e98d76cbb 100644 --- a/src/util/lp/numeric_pair.h +++ b/src/util/lp/numeric_pair.h @@ -285,7 +285,9 @@ class numeric_traits> { return numeric_traits::is_neg(p.x) || (numeric_traits::is_zero(p.x) && numeric_traits::is_neg(p.y)); } - + static bool is_int(const numeric_pair & p) { + return numeric_traits::is_int(p.x) && numeric_traits::is_zero(p.y); + } }; template <> diff --git a/src/util/lp/var_register.h b/src/util/lp/var_register.h index 7528e8d49..86c238e12 100644 --- a/src/util/lp/var_register.h +++ b/src/util/lp/var_register.h @@ -18,33 +18,94 @@ Revision History: --*/ #pragma once namespace lp { +class ext_var_info { + unsigned m_internal_j; // the internal index + bool m_is_integer; +public: + ext_var_info() {} + ext_var_info(unsigned j): ext_var_info(j, true) {} + ext_var_info(unsigned j , bool is_int) : m_internal_j(j), m_is_integer(is_int) {} + unsigned internal_j() const { return m_internal_j;} + bool is_integer() const {return m_is_integer;} +}; + class var_register { - svector m_local_vars_to_external; - std::unordered_map m_external_vars_to_local; + svector m_local_to_external; + std::unordered_map m_external_to_local; public: unsigned add_var(unsigned user_var) { - auto t = m_external_vars_to_local.find(user_var); - if (t != m_external_vars_to_local.end()) { - return t->second; + return add_var(user_var, true); + } + unsigned add_var(unsigned user_var, bool is_int) { + auto t = m_external_to_local.find(user_var); + if (t != m_external_to_local.end()) { + return t->second.internal_j(); } - unsigned ret = size(); - m_external_vars_to_local[user_var] = ret; - m_local_vars_to_external.push_back(user_var); - return ret; + unsigned j = size(); + m_external_to_local[user_var] = ext_var_info(j, is_int); + m_local_to_external.push_back(user_var); + return j; } - const svector & vars() const { return m_local_vars_to_external; } + const svector & vars() const { return m_local_to_external; } - unsigned local_var_to_user_var(unsigned local_var) const { - return m_local_vars_to_external[local_var]; + unsigned local_to_external(unsigned local_var) const { + return m_local_to_external[local_var]; } unsigned size() const { - return m_local_vars_to_external.size(); + return m_local_to_external.size(); } void clear() { - m_local_vars_to_external.clear(); - m_external_vars_to_local.clear(); + m_local_to_external.clear(); + m_external_to_local.clear(); } + + unsigned external_to_local(unsigned j) const { + auto it = m_external_to_local.find(j); + lp_assert(it != m_external_to_local.end()); + return it->second.internal_j(); + } + + bool external_is_int(unsigned j) const { + auto it = m_external_to_local.find(j); + lp_assert(it != m_external_to_local.end()); + return it->second.is_integer(); + } + + bool external_is_used(unsigned ext_j) const { + auto it = m_external_to_local.find(ext_j); + return it != m_external_to_local.end(); + } + + bool external_is_used(unsigned ext_j, unsigned & local_j ) const { + auto it = m_external_to_local.find(ext_j); + if ( it == m_external_to_local.end()) + return false; + local_j = it->second.internal_j(); + return true; + } + + bool external_is_used(unsigned ext_j, unsigned & local_j, bool & is_int ) const { + auto it = m_external_to_local.find(ext_j); + if ( it == m_external_to_local.end()) + return false; + local_j = it->second.internal_j(); + is_int = it->second.is_integer(); + return true; + } + + + bool local_is_int(unsigned j) const { + return external_is_int(m_local_to_external[j]); + } + + void shrink(unsigned shrunk_size) { + for (unsigned j = size(); j-- > shrunk_size;) { + m_external_to_local.erase(m_local_to_external[j]); + } + m_local_to_external.resize(shrunk_size); + } + }; }