3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-10 03:07:07 +00:00

use heap to track infeasible columns. (#6771)

* use heap to track infeasible columns

* fix the formatting

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

---------

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2023-06-19 13:50:14 -07:00 committed by GitHub
parent 4d44e60c33
commit 32ec02778e
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
12 changed files with 116 additions and 102 deletions

View file

@ -279,9 +279,9 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print()
print_row(i);
}
m_out << std::endl;
if (m_core_solver.inf_set().size()) {
if (m_core_solver.inf_heap().size()) {
m_out << "inf columns: ";
print_vector(m_core_solver.inf_set(), m_out);
print_vector(m_core_solver.inf_heap(), m_out);
m_out << std::endl;
}
}

View file

@ -85,8 +85,8 @@ unsigned lar_core_solver::get_number_of_non_ints() const {
void lar_core_solver::solve() {
TRACE("lar_solver", tout << m_r_solver.get_status() << "\n";);
lp_assert(m_r_solver.non_basic_columns_are_set_correctly());
lp_assert(m_r_solver.inf_set_is_correct());
TRACE("find_feas_stats", tout << "infeasibles = " << m_r_solver.inf_set_size() << ", int_infs = " << get_number_of_non_ints() << std::endl;);
lp_assert(m_r_solver.inf_heap_is_correct());
TRACE("find_feas_stats", tout << "infeasibles = " << m_r_solver.inf_heap_size() << ", int_infs = " << get_number_of_non_ints() << std::endl;);
if (m_r_solver.current_x_is_feasible() && m_r_solver.m_look_for_feasible_solution_only) {
m_r_solver.set_status(lp_status::OPTIMAL);
TRACE("lar_solver", tout << m_r_solver.get_status() << "\n";);
@ -117,11 +117,9 @@ void lar_core_solver::solve() {
}
lp_assert(r_basis_is_OK());
lp_assert(m_r_solver.non_basic_columns_are_set_correctly());
lp_assert(m_r_solver.inf_set_is_correct());
TRACE("lar_solver", tout << m_r_solver.get_status() << "\n";);
}
lp_assert(m_r_solver.inf_heap_is_correct());
TRACE("lar_solver", tout << m_r_solver.get_status() << "\n";);
}
} // namespace lp

View file

@ -244,6 +244,14 @@ namespace lp {
set.erase(j);
}
void lar_solver::clean_popped_elements_for_heap(unsigned n, lpvar_heap& heap) {
vector<int> to_remove;
for (unsigned j : heap)
if (j >= n)
to_remove.push_back(j);
for (unsigned j : to_remove)
heap.erase(j);
}
void lar_solver::pop(unsigned k) {
@ -271,7 +279,7 @@ namespace lp {
unsigned m = A_r().row_count();
clean_popped_elements(m, m_rows_with_changed_bounds);
clean_inf_set_of_r_solver_after_pop();
clean_inf_heap_of_r_solver_after_pop();
lp_assert(
m_settings.simplex_strategy() == simplex_strategy_enum::undecided ||
m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
@ -328,7 +336,7 @@ namespace lp {
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.inf_set(); // hijack this set that should be empty right now
auto& jset = m_mpq_lar_core_solver.m_r_solver.inf_heap(); // hijack this set that should be empty right now
lp_assert(jset.empty());
for (lar_term::ival p : term) {
@ -674,9 +682,9 @@ namespace lp {
void lar_solver::update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j) {
if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) {
if (costs_are_used()) {
bool was_infeas = m_mpq_lar_core_solver.m_r_solver.inf_set_contains(j);
bool was_infeas = m_mpq_lar_core_solver.m_r_solver.inf_heap_contains(j);
m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j);
if (was_infeas != m_mpq_lar_core_solver.m_r_solver.inf_set_contains(j))
if (was_infeas != m_mpq_lar_core_solver.m_r_solver.inf_heap_contains(j))
m_basic_columns_with_changed_cost.insert(j);
}
else {
@ -1293,12 +1301,12 @@ namespace lp {
lp_assert(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct());
}
void lar_solver::clean_inf_set_of_r_solver_after_pop() {
void lar_solver::clean_inf_heap_of_r_solver_after_pop() {
vector<unsigned> became_feas;
clean_popped_elements(A_r().column_count(), m_mpq_lar_core_solver.m_r_solver.inf_set());
clean_popped_elements_for_heap(A_r().column_count(), m_mpq_lar_core_solver.m_r_solver.inf_heap());
std::unordered_set<unsigned> basic_columns_with_changed_cost;
m_inf_index_copy.reset();
for (auto j : m_mpq_lar_core_solver.m_r_solver.inf_set())
for (auto j : m_mpq_lar_core_solver.m_r_solver.inf_heap())
m_inf_index_copy.push_back(j);
for (auto j : m_inf_index_copy) {
if (m_mpq_lar_core_solver.m_r_heading[j] >= 0) {
@ -1316,16 +1324,16 @@ namespace lp {
lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis_heading[j] < 0);
m_mpq_lar_core_solver.m_r_solver.m_d[j] -= m_mpq_lar_core_solver.m_r_solver.m_costs[j];
m_mpq_lar_core_solver.m_r_solver.m_costs[j] = zero_of_type<mpq>();
m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_set(j);
m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_heap(j);
}
became_feas.clear();
for (unsigned j : m_mpq_lar_core_solver.m_r_solver.inf_set()) {
for (unsigned j : m_mpq_lar_core_solver.m_r_solver.inf_heap()) {
lp_assert(m_mpq_lar_core_solver.m_r_heading[j] >= 0);
if (m_mpq_lar_core_solver.m_r_solver.column_is_feasible(j))
became_feas.push_back(j);
}
for (unsigned j : became_feas)
m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_set(j);
m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_heap(j);
}
@ -1504,7 +1512,7 @@ namespace lp {
m_mpq_lar_core_solver.m_r_x.resize(j + 1);
m_mpq_lar_core_solver.m_r_lower_bounds.increase_size_by_one();
m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one();
m_mpq_lar_core_solver.m_r_solver.inf_set_increase_size_by_one();
m_mpq_lar_core_solver.m_r_solver.inf_heap_increase_size_by_one();
m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1);
m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1);
lp_assert(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method

View file

@ -94,7 +94,7 @@ class lar_solver : public column_namer {
// these are basic columns with the value changed, so the corresponding row in the tableau
// does not sum to zero anymore
u_set m_incorrect_columns;
// copy of m_r_solver.inf_set()
// copy of m_r_solver.inf_heap()
unsigned_vector m_inf_index_copy;
stacked_value<unsigned> m_term_count;
vector<lar_term*> m_terms;
@ -178,7 +178,9 @@ class lar_solver : public column_namer {
bp);
}
static void clean_popped_elements_for_heap(unsigned n, lpvar_heap& set);
static void clean_popped_elements(unsigned n, u_set& set);
bool maximize_term_on_tableau(const lar_term & term,
impq &term_max);
bool costs_are_zeros_for_r_solver() const;
@ -230,7 +232,7 @@ class lar_solver : public column_namer {
void remove_last_column_from_basis_tableau(unsigned j);
void remove_last_column_from_tableau();
void pop_tableau();
void clean_inf_set_of_r_solver_after_pop();
void clean_inf_heap_of_r_solver_after_pop();
inline bool column_value_is_integer(unsigned j) const { return get_column_value(j).is_int(); }
bool model_is_int_feasible() const;
@ -386,14 +388,15 @@ public:
return tv::is_term(idx)?
m_term_register.local_to_external(idx) : m_var_register.local_to_external(idx);
}
bool column_corresponds_to_term(unsigned) const;
inline unsigned row_count() const { return A_r().row_count(); }
bool var_is_registered(var_index vj) const;
void clear_inf_set() {
m_mpq_lar_core_solver.m_r_solver.inf_set().clear();
void clear_inf_heap() {
m_mpq_lar_core_solver.m_r_solver.inf_heap().clear();
}
inline void remove_column_from_inf_set(unsigned j) {
m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_set(j);
inline void remove_column_from_inf_heap(unsigned j) {
m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_heap(j);
}
template <typename ChangeReport>
void change_basic_columns_dependend_on_a_given_nb_column_report(unsigned j,

View file

@ -63,8 +63,8 @@ template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::colu
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq>>::pivot_column_tableau(unsigned int, unsigned int);
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::pivot_column_tableau(unsigned int, unsigned int);
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::transpose_rows_tableau(unsigned int, unsigned int);
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::inf_set_is_correct() const;
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::inf_set_is_correct() const;
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::inf_heap_is_correct() const;
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::inf_heap_is_correct() const;
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::remove_from_basis(unsigned int);

View file

@ -28,9 +28,13 @@ Revision History:
#include "math/lp/permutation_matrix.h"
#include "math/lp/column_namer.h"
#include "math/lp/u_set.h"
#include "util/heap.h"
namespace lp {
struct lpvar_lt {
bool operator()(lpvar v1, lpvar v2) const { return v1 < v2; }
};
typedef heap<lpvar_lt> lpvar_heap;
template <typename T, typename X>
X dot_product(const vector<T> & a, const vector<X> & b) {
lp_assert(a.size() == b.size());
@ -51,24 +55,24 @@ private:
public:
bool current_x_is_feasible() const {
TRACE("feas",
if (m_inf_set.size()) {
tout << "column " << *m_inf_set.begin() << " is infeasible" << std::endl;
print_column_info(*m_inf_set.begin(), tout);
if (m_inf_heap.size()) {
tout << "column " << *m_inf_heap.begin() << " is infeasible" << std::endl;
print_column_info(*m_inf_heap.begin(), tout);
} else {
tout << "x is feasible\n";
}
);
return m_inf_set.size() == 0;
return m_inf_heap.empty();
}
bool current_x_is_infeasible() const { return m_inf_set.size() != 0; }
bool current_x_is_infeasible() const { return m_inf_heap.size() != 0; }
private:
u_set m_inf_set;
lpvar_heap m_inf_heap;
public:
const u_set& inf_set() const { return m_inf_set; }
u_set& inf_set() { return m_inf_set; }
void inf_set_increase_size_by_one() { m_inf_set.increase_size_by_one(); }
bool inf_set_contains(unsigned j) const { return m_inf_set.contains(j); }
unsigned inf_set_size() const { return m_inf_set.size(); }
const lpvar_heap& inf_heap() const { return m_inf_heap; }
lpvar_heap& inf_heap() { return m_inf_heap; }
void inf_heap_increase_size_by_one() { m_inf_heap.reserve(m_inf_heap.size() + 1); }
bool inf_heap_contains(unsigned j) const { return m_inf_heap.contains(j); }
unsigned inf_heap_size() const { return m_inf_heap.size(); }
indexed_vector<T> m_pivot_row; // this is the real pivot row of the simplex tableu
static_matrix<T, X> & m_A; // the matrix A
// vector<X> const & m_b; // the right side
@ -255,7 +259,7 @@ public:
bool calc_current_x_is_feasible_include_non_basis() const;
bool inf_set_is_correct() const;
bool inf_heap_is_correct() const;
bool column_is_dual_feasible(unsigned j) const;
@ -526,8 +530,8 @@ public:
swap(m_basis_heading, m_basis[i], m_basis[ii]);
}
bool column_is_in_inf_set(unsigned j) const {
return m_inf_set.contains(j);
bool column_is_in_inf_heap(unsigned j) const {
return m_inf_heap.contains(j);
}
bool column_is_base(unsigned j) const {
@ -560,29 +564,26 @@ public:
void track_column_feasibility(unsigned j) {
if (column_is_feasible(j))
remove_column_from_inf_set(j);
remove_column_from_inf_heap(j);
else
insert_column_into_inf_set(j);
insert_column_into_inf_heap(j);
}
void insert_column_into_inf_set(unsigned j) {
void insert_column_into_inf_heap(unsigned j) {
TRACE("lar_solver", tout << "j = " << j << "\n";);
m_inf_set.insert(j);
if (!m_inf_heap.contains(j))
m_inf_heap.insert(j);
lp_assert(!column_is_feasible(j));
}
void remove_column_from_inf_set(unsigned j) {
void remove_column_from_inf_heap(unsigned j) {
TRACE("lar_solver", tout << "j = " << j << "\n";);
m_inf_set.erase(j);
if (m_inf_heap.contains(j))
m_inf_heap.erase(j);
lp_assert(column_is_feasible(j));
}
void resize_inf_set(unsigned size) {
void clear_inf_heap() {
TRACE("lar_solver",);
m_inf_set.resize(size);
}
void clear_inf_set() {
TRACE("lar_solver",);
m_inf_set.clear();
m_inf_heap.clear();
}
bool costs_on_nbasis_are_zeros() const {

View file

@ -42,7 +42,7 @@ lp_core_solver_base(static_matrix<T, X> & A,
m_total_iterations(0),
m_iters_with_no_cost_growing(0),
m_status(lp_status::FEASIBLE),
m_inf_set(A.column_count()),
m_inf_heap(std::max(static_cast<unsigned>(1024), A.column_count())),
m_pivot_row(A.column_count()),
m_A(A),
m_basis(basis),
@ -250,9 +250,9 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::calc_current_x
return true;
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::inf_set_is_correct() const {
template <typename T, typename X> bool lp_core_solver_base<T, X>::inf_heap_is_correct() const {
for (unsigned j = 0; j < this->m_n(); j++) {
bool belongs_to_set = m_inf_set.contains(j);
bool belongs_to_set = m_inf_heap.contains(j);
bool is_feas = column_is_feasible(j);
if (is_feas == belongs_to_set) {
TRACE("lp_core", tout << "incorrectly set column in inf set "; print_column_info(j, tout) << "\n";);

View file

@ -33,6 +33,8 @@ Revision History:
#include <sstream>
#include <string>
#include <unordered_map>
#include "util/heap.h"
namespace lp {
// This core solver solves (Ax=b, lower_bound_values \leq x \leq
@ -152,19 +154,28 @@ public:
UNREACHABLE(); // unreachable
return false;
}
unsigned get_number_of_basic_vars_that_might_become_inf(
unsigned j) const { // consider looking at the signs here: todo
/**
* Return the number of base non-free variables depending on the column j,
* different from bj,
* but take the min with the (bound+1).
* This function is used to select the pivot variable.
*/
unsigned get_num_of_not_free_basic_dependent_vars(
unsigned j, unsigned bound, unsigned bj) const { // consider looking at the signs here: todo
unsigned r = 0;
for (const auto &cc : this->m_A.m_columns[j]) {
unsigned k = this->m_basis[cc.var()];
if (this->m_column_types[k] != column_type::free_column)
r++;
unsigned basic_for_row = this->m_basis[cc.var()];
if (basic_for_row == bj)
continue;
// std::cout << this->m_A.m_rows[cc.var()] << std::endl;
if (this->m_column_types[basic_for_row] != column_type::free_column)
if (r++ > bound) return r;
}
return r;
}
int find_beneficial_column_in_row_tableau_rows_bland_mode(int i, T &a_ent) {
int find_beneficial_entering_in_row_tableau_rows_bland_mode(int i, T &a_ent) {
int j = -1;
unsigned bj = this->m_basis[i];
bool bj_needs_to_grow = needs_to_grow(bj);
@ -190,15 +201,15 @@ public:
return j;
}
int find_beneficial_column_in_row_tableau_rows(int i, T &a_ent) {
int find_beneficial_entering_tableau_rows(int i, T &a_ent) {
if (m_bland_mode_tableau)
return find_beneficial_column_in_row_tableau_rows_bland_mode(i, a_ent);
return find_beneficial_entering_in_row_tableau_rows_bland_mode(i, a_ent);
// a short row produces short infeasibility explanation and benefits at
// least one pivot operation
int choice = -1;
int nchoices = 0;
unsigned num_of_non_free_basics = 1000000;
unsigned len = 100000000;
unsigned min_non_free_so_far = -1;
unsigned best_col_sz = -1;
unsigned bj = this->m_basis[i];
bool bj_needs_to_grow = needs_to_grow(bj);
for (unsigned k = 0; k < this->m_A.m_rows[i].size(); k++) {
@ -213,17 +224,18 @@ public:
if (!monoid_can_increase(rc))
continue;
}
unsigned damage = get_number_of_basic_vars_that_might_become_inf(j);
if (damage < num_of_non_free_basics) {
num_of_non_free_basics = damage;
len = this->m_A.m_columns[j].size();
unsigned not_free = get_num_of_not_free_basic_dependent_vars(j, min_non_free_so_far, bj);
unsigned col_sz = this->m_A.m_columns[j].size();
if (not_free < min_non_free_so_far || (not_free == min_non_free_so_far && col_sz < best_col_sz)) {
min_non_free_so_far = not_free;
best_col_sz = this->m_A.m_columns[j].size();
choice = k;
nchoices = 1;
} else if (damage == num_of_non_free_basics &&
this->m_A.m_columns[j].size() <= len &&
(this->m_settings.random_next() % (++nchoices))) {
choice = k;
len = this->m_A.m_columns[j].size();
} else if (not_free == min_non_free_so_far &&
col_sz == best_col_sz) {
if (this->m_settings.random_next(++nchoices) == 0){
choice = k;
}
}
}
@ -282,7 +294,7 @@ public:
unsigned j, vector<unsigned> &leavings, T m, X &t,
T &abs_of_d_of_leaving);
X get_max_bound(vector<X> &b);
#ifdef Z3DEBUG
void check_Ax_equal_b();
@ -291,14 +303,6 @@ public:
void check_correctness();
#endif
// from page 183 of Istvan Maros's book
// the basis structures have not changed yet
void update_reduced_costs_from_pivot_row(unsigned entering, unsigned leaving);
// return 0 if the reduced cost at entering is close enough to the refreshed
// 1 if it is way off, and 2 if it is unprofitable
int refresh_reduced_cost_at_entering_and_check_that_it_is_off(
unsigned entering);
void backup_and_normalize_costs();
@ -356,15 +360,13 @@ public:
}
int find_smallest_inf_column() {
int j = -1;
for (unsigned k : this->inf_set()) {
if (k < static_cast<unsigned>(j)) {
j = k;
}
}
return j;
if (this->inf_heap().empty())
return -1;
return this->inf_heap().min_value();
}
const X &get_val_for_leaving(unsigned j) const {
lp_assert(!this->column_is_feasible(j));
switch (this->m_column_types[j]) {
@ -405,7 +407,7 @@ public:
}
}
T a_ent;
int entering = find_beneficial_column_in_row_tableau_rows(
int entering = find_beneficial_entering_tableau_rows(
this->m_basis_heading[leaving], a_ent);
if (entering == -1) {
this->set_status(lp_status::INFEASIBLE);
@ -414,7 +416,7 @@ public:
const X &new_val_for_leaving = get_val_for_leaving(leaving);
X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent;
this->m_x[leaving] = new_val_for_leaving;
this->remove_column_from_inf_set(leaving);
this->remove_column_from_inf_heap(leaving);
advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta);
if (this->current_x_is_feasible())
this->set_status(lp_status::OPTIMAL);

View file

@ -30,7 +30,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::one_iteratio
else {
advance_on_entering_tableau(entering);
}
lp_assert(this->inf_set_is_correct());
lp_assert(this->inf_heap_is_correct());
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_entering_tableau(int entering) {
@ -256,7 +256,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tab
this->m_basis_sort_counter = 0; // to initiate the sort of the basis
// this->set_total_iterations(0);
this->iters_with_no_cost_growing() = 0;
lp_assert(this->inf_set_is_correct());
lp_assert(this->inf_heap_is_correct());
if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only)
return;
if (this->m_settings.backup_costs)

View file

@ -218,6 +218,8 @@ public:
unsigned hnf_cut_period() const { return m_hnf_cut_period; }
void set_hnf_cut_period(unsigned period) { m_hnf_cut_period = period; }
unsigned random_next() { return m_rand(); }
unsigned random_next(unsigned u ) { return m_rand(u); }
void set_random_seed(unsigned s) { m_rand.set_seed(s); }
bool bound_progation() const {

View file

@ -1436,7 +1436,7 @@ void core::patch_monomials() {
} else {
m_lar_solver.pop();
restore_tableau();
m_lar_solver.clear_inf_set();
m_lar_solver.clear_inf_heap();
}
SASSERT(m_lar_solver.ax_is_correct());
}

View file

@ -75,8 +75,8 @@ private:
c().negate_relation(lemma, m_jy, m_y.rat_sign()*pl.y);
#if Z3DEBUG
SASSERT(c().val(m_x) == m_xy.x && c().val(m_y) == m_xy.y);
int mult_sign = nla::rat_sign(pl.x - m_xy.x)*nla::rat_sign(pl.y - m_xy.y);
SASSERT((mult_sign == 1) == m_below);
// int mult_sign = nla::rat_sign(pl.x - m_xy.x)*nla::rat_sign(pl.y - m_xy.y);
SASSERT((nla::rat_sign(pl.x - m_xy.x)*nla::rat_sign(pl.y - m_xy.y) == 1) == m_below);
// If "mult_sign is 1" then (a - x)(b-y) > 0 and ab - bx - ay + xy > 0
// or -ab + bx + ay < xy or -ay - bx + xy > -ab
// val(j) stands for xy. So, finally we have -ay - bx + j > - ab