mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-23 17:15:31 +00:00
Code Activity

working on stand-alone simplex

Browse source 
Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
This commit is contained in:
Nikolaj Bjorner 2014-01-26 19:46:42 -08:00
parent f68eff3276
commit c14c65465a
9 changed files with 332 additions and 50 deletions
Download patch file Download diff file Expand all files Collapse all files

28
src/math/simplex/simplex.h
View file

@ -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<var_info> m_vars;
svector<var_t> 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<bool is_below>
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;
};
};

308
src/math/simplex/simplex_def.h
View file

@ -44,7 +44,8 @@ namespace simplex {
template<typename Ext>
void simplex<Ext>::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<typename Ext>
lbool simplex<Ext>::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<typename Ext>
void simplex<Ext>::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<typename Ext>
void simplex<Ext>::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<typename Ext>
bool simplex<Ext>::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<typename Ext>
bool simplex<Ext>::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<typename Ext>
bool simplex<Ext>::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 Ext>
typename simplex<Ext>::var_t
simplex<Ext>::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<true>(x_i, out_a_ij);
}
else {
return select_pivot_core<false>(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<typename Ext>
template<bool is_below>
typename simplex<Ext>::var_t
simplex<Ext>::select_pivot_core(var_t x_i, scoped_numeral & out_a_ij) {
simplex<Ext>::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 Ext>
typename simplex<Ext>::var_t
simplex<Ext>::select_blands_pivot(var_t x_i, bool is_below, scoped_numeral & out_a_ij) {
simplex<Ext>::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<typename Ext>
lbool simplex<Ext>::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<Ext>::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<typename Ext>
void simplex<Ext>::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<typename Ext>
void simplex<Ext>::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 Ext>
typename simplex<Ext>::var_t
simplex<Ext>::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<typename Ext>
void simplex<Ext>::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 Ext>
typename simplex<Ext>::pivot_strategy_t
simplex<Ext>::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<typename Ext>
bool simplex<Ext>::is_feasible() const {
for (unsigned i = 0; i < m_vars.size(); ++i) {
if (below_lower(i) || above_upper(i)) return false;
}
return true;
}
};

2
src/math/simplex/sparse_matrix.h
View file

@ -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;

6
src/math/simplex/sparse_matrix_def.h
View file

@ -448,7 +448,7 @@ namespace simplex {
*/
template<typename Ext>
void sparse_matrix<Ext>::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<typename Ext>
void sparse_matrix<Ext>::display(std::ostream& out) const {
void sparse_matrix<Ext>::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";