mirror of
https://github.com/Z3Prover/z3
synced 2025-04-27 02:45:51 +00:00
add stand-alone simplex
Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
This commit is contained in:
Nikolaj Bjorner
parent
c6a9dae00a
commit
26a3d2ca31
9 changed files with 1560 additions and 20 deletions
|
|
@ -15,6 +15,7 @@ def init_project_def():
|
|||
add_lib('sat', ['util'])
|
||||
add_lib('nlsat', ['polynomial', 'sat'])
|
||||
add_lib('hilbert', ['util'], 'math/hilbert')
|
||||
add_lib('simplex', ['util'], 'math/simplex')
|
||||
add_lib('interval', ['util'], 'math/interval')
|
||||
add_lib('realclosure', ['interval'], 'math/realclosure')
|
||||
add_lib('subpaving', ['interval'], 'math/subpaving')
|
||||
|
|
@ -75,7 +76,7 @@ def init_project_def():
|
|||
add_lib('api', ['portfolio', 'user_plugin', 'smtparser', 'realclosure','interp','opt'],
|
||||
includes2install=['z3.h', 'z3_v1.h', 'z3_macros.h'] + API_files)
|
||||
add_exe('shell', ['api', 'sat', 'extra_cmds','opt'], exe_name='z3')
|
||||
add_exe('test', ['api', 'fuzzing'], exe_name='test-z3', install=False)
|
||||
add_exe('test', ['api', 'fuzzing', 'simplex'], exe_name='test-z3', install=False)
|
||||
add_dll('api_dll', ['api', 'sat', 'extra_cmds'], 'api/dll',
|
||||
reexports=['api'],
|
||||
dll_name='libz3',
|
||||
|
|
|
|||
26
src/math/simplex/simplex.cpp
Normal file
26
src/math/simplex/simplex.cpp
Normal file
|
|
@ -0,0 +1,26 @@
|
|||
/*++
|
||||
Copyright (c) 2014 Microsoft Corporation
|
||||
|
||||
Module Name:
|
||||
|
||||
simplex.h
|
||||
|
||||
Abstract:
|
||||
|
||||
Multi-precision simplex tableau.
|
||||
|
||||
Author:
|
||||
|
||||
Nikolaj Bjorner (nbjorner) 2014-01-15
|
||||
|
||||
Notes:
|
||||
|
||||
--*/
|
||||
|
||||
#include"simplex.h"
|
||||
#include"sparse_matrix_def.h"
|
||||
#include"simplex_def.h"
|
||||
namespace simplex {
|
||||
template class simplex<mpz_ext>;
|
||||
template class simplex<mpq_ext>;
|
||||
};
|
||||
157
src/math/simplex/simplex.h
Normal file
157
src/math/simplex/simplex.h
Normal file
|
|
@ -0,0 +1,157 @@
|
|||
/*++
|
||||
Copyright (c) 2014 Microsoft Corporation
|
||||
|
||||
Module Name:
|
||||
|
||||
simplex.h
|
||||
|
||||
Abstract:
|
||||
|
||||
Multi-precision simplex tableau.
|
||||
|
||||
- It uses code from theory_arith where applicable.
|
||||
|
||||
- It is detached from the theory class and ASTs.
|
||||
|
||||
- It uses non-shared mpz/mpq's avoiding global locks and operations on rationals.
|
||||
|
||||
- It follows the same sparse tableau layout (no LU yet).
|
||||
|
||||
- It does not include features for non-linear arithmetic.
|
||||
|
||||
- Branch/bound/cuts is external.
|
||||
|
||||
Author:
|
||||
|
||||
Nikolaj Bjorner (nbjorner) 2014-01-15
|
||||
|
||||
Notes:
|
||||
|
||||
--*/
|
||||
|
||||
#ifndef _SIMPLEX_H_
|
||||
#define _SIMPLEX_H_
|
||||
|
||||
#include "sparse_matrix.h"
|
||||
#include "mpq_inf.h"
|
||||
#include "heap.h"
|
||||
#include "lbool.h"
|
||||
|
||||
namespace simplex {
|
||||
|
||||
template<typename Ext>
|
||||
class simplex {
|
||||
|
||||
typedef unsigned var_t;
|
||||
typedef typename Ext::eps_numeral eps_numeral;
|
||||
typedef typename Ext::numeral numeral;
|
||||
typedef typename Ext::manager manager;
|
||||
typedef typename Ext::eps_manager eps_manager;
|
||||
typedef typename Ext::scoped_numeral scoped_numeral;
|
||||
typedef _scoped_numeral<eps_manager> scoped_eps_numeral;
|
||||
|
||||
typedef typename _scoped_numeral_vector<eps_manager> scoped_eps_numeral_vector;
|
||||
|
||||
typedef sparse_matrix<Ext> matrix;
|
||||
|
||||
struct var_lt {
|
||||
bool operator()(var_t v1, var_t v2) const { return v1 < v2; }
|
||||
};
|
||||
typedef heap<var_lt> var_heap;
|
||||
|
||||
enum pivot_strategy_t {
|
||||
S_BLAND,
|
||||
S_GREATEST_ERROR,
|
||||
S_LEAST_ERROR,
|
||||
S_DEFAULT
|
||||
};
|
||||
|
||||
struct var_info {
|
||||
unsigned m_base2row:29;
|
||||
unsigned m_is_base:1;
|
||||
unsigned m_lower_valid:1;
|
||||
unsigned m_upper_valid:1;
|
||||
eps_numeral m_value;
|
||||
eps_numeral m_lower;
|
||||
eps_numeral m_upper;
|
||||
numeral m_base_coeff;
|
||||
var_info():
|
||||
m_base2row(0),
|
||||
m_is_base(false),
|
||||
m_lower_valid(false),
|
||||
m_upper_valid(false)
|
||||
{}
|
||||
};
|
||||
|
||||
static const var_t null_var = UINT_MAX;
|
||||
matrix M;
|
||||
unsigned m_max_iterations;
|
||||
volatile bool m_cancel;
|
||||
var_heap m_to_patch;
|
||||
mutable manager m;
|
||||
mutable eps_manager em;
|
||||
vector<var_info> m_vars;
|
||||
svector<var_t> m_row2base;
|
||||
bool m_bland;
|
||||
random_gen m_random;
|
||||
|
||||
public:
|
||||
simplex():
|
||||
m_max_iterations(UINT_MAX),
|
||||
m_cancel(false),
|
||||
m_to_patch(1024),
|
||||
m_bland(false) {}
|
||||
|
||||
typedef typename matrix::row row;
|
||||
typedef typename matrix::row_iterator row_iterator;
|
||||
typedef typename matrix::col_iterator col_iterator;
|
||||
|
||||
void ensure_var(var_t v);
|
||||
row add_row(unsigned num_vars, var_t base, var_t const* vars, numeral const* coeffs);
|
||||
void del_row(row const& r);
|
||||
void set_lower(var_t var, eps_numeral const& b);
|
||||
void set_upper(var_t var, eps_numeral const& b);
|
||||
void unset_lower(var_t var);
|
||||
void unset_upper(var_t var);
|
||||
void set_value(var_t var, eps_numeral const& b);
|
||||
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);
|
||||
eps_numeral const& get_value(var_t v);
|
||||
void display(std::ostream& out) const;
|
||||
|
||||
private:
|
||||
|
||||
var_t select_var_to_fix();
|
||||
pivot_strategy_t pivot_strategy();
|
||||
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) { }
|
||||
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);
|
||||
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);
|
||||
int get_num_non_free_dep_vars(var_t x_j, int best_so_far);
|
||||
|
||||
bool below_lower(var_t v) const;
|
||||
bool above_upper(var_t v) const;
|
||||
bool above_lower(var_t v) const;
|
||||
bool below_upper(var_t v) const;
|
||||
bool outside_bounds(var_t v) const { return below_lower(v) || above_upper(v); }
|
||||
bool is_free(var_t v) const { return !m_vars[v].m_lower_valid && !m_vars[v].m_upper_valid; }
|
||||
bool is_non_free(var_t v) const { return !is_free(v); }
|
||||
unsigned get_num_vars() const { return m_vars.size(); }
|
||||
bool is_base(var_t x) const { return m_vars[x].m_is_base; }
|
||||
|
||||
bool well_formed() const;
|
||||
};
|
||||
|
||||
};
|
||||
|
||||
#endif
|
||||
480
src/math/simplex/simplex_def.h
Normal file
480
src/math/simplex/simplex_def.h
Normal file
|
|
@ -0,0 +1,480 @@
|
|||
/*++
|
||||
Copyright (c) 2014 Microsoft Corporation
|
||||
|
||||
Module Name:
|
||||
|
||||
simplex_def.h
|
||||
|
||||
Abstract:
|
||||
|
||||
Author:
|
||||
|
||||
Nikolaj Bjorner (nbjorner) 2014-01-15
|
||||
|
||||
Notes:
|
||||
|
||||
--*/
|
||||
|
||||
#ifndef _SIMPLEX_DEF_H_
|
||||
#define _SIMPLEX_DEF_H_
|
||||
|
||||
|
||||
namespace simplex {
|
||||
|
||||
template<typename Ext>
|
||||
typename simplex<Ext>::row
|
||||
simplex<Ext>::add_row(var_t base, unsigned num_vars, var_t const* vars, numeral const* coeffs) {
|
||||
DEBUG_CODE(
|
||||
bool found = false;
|
||||
for (unsigned i = 0; !found && i < num_vars; ++i) found = vars[i] == base;
|
||||
SASSERT(found);
|
||||
);
|
||||
row r = M.mk_row();
|
||||
for (unsigned i = 0; i < num_vars; ++i) {
|
||||
M.add(r, coeffs[i], vars[i]);
|
||||
}
|
||||
while (m_row2base.size() <= r.id()) {
|
||||
m_row2base.push_back(null_var);
|
||||
}
|
||||
m_row2base[r.id()] = base;
|
||||
m_vars[base].m_base2row = r.id();
|
||||
m_vars[base].m_is_base = true;
|
||||
return r;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::del_row(row const& r) {
|
||||
m_is_base[m_row2base[r.id()]] = false;
|
||||
M.del(r);
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::set_lower(var_t var, eps_numeral const& b) {
|
||||
var_info& vi = m_vars[var];
|
||||
em.set(vi.m_lower, b);
|
||||
vi.m_lower_valid = true;
|
||||
if (em.lt(vi.m_value, b)) {
|
||||
scoped_eps_numeral delta(em);
|
||||
em.sub(b, vi.m_value, delta);
|
||||
update_value(var, delta);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::set_upper(var_t var, eps_numeral const& b) {
|
||||
var_info& vi = m_vars[var];
|
||||
em.set(vi.m_upper, b);
|
||||
vi.m_upper_valid = true;
|
||||
if (em.gt(vi.m_value, b)) {
|
||||
scoped_eps_numeral delta(em);
|
||||
em.sub(b, vi.m_value, delta);
|
||||
update_value(var, delta);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::unset_lower(var_t var) {
|
||||
m_vars[var].m_lower_valid = false;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::unset_upper(var_t var) {
|
||||
m_vars[var].m_upper_valid = false;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::set_value(var_t var, eps_numeral const& b) {
|
||||
scoped_eps_numeral delta(em);
|
||||
em.sub(b, m_vars[var].m_value, delta);
|
||||
update_value(var, delta);
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
typename simplex<Ext>::eps_numeral const&
|
||||
simplex<Ext>::get_value(var_t v) {
|
||||
return m_vars[v].m_value;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::display(std::ostream& out) const {
|
||||
M.display(out);
|
||||
for (unsigned i = 0; i < m_vars.size(); ++i) {
|
||||
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 << " [";
|
||||
if (vi.m_lower_valid) em.display(out, vi.m_lower) else out << "-oo";
|
||||
out << ":";
|
||||
if (vi.m_upper_valid) em.display(out, vi.m_upper) else out << "oo";
|
||||
out << "]";
|
||||
out << "\n";
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::ensure_var(var_t v) {
|
||||
while (v >= m_vars.size()) {
|
||||
M.ensure_var(m_vars.size());
|
||||
m_vars.push_back(var_info());
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
lbool simplex<Ext>::make_feasible() {
|
||||
unsigned num_iterations = 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);
|
||||
if (!make_var_feasible(v)) {
|
||||
return l_false;
|
||||
}
|
||||
++num_iterations;
|
||||
SASSERT(well_formed());
|
||||
}
|
||||
return l_true;
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Make x_j the new base variable for row of x_i.
|
||||
x_i is assumed base variable of row r_i.
|
||||
x_j is assumed to have coefficient a_ij in r_i.
|
||||
|
||||
a_ii*x_i + a_ij*x_j + r_i = 0
|
||||
|
||||
current value of x_i is v_i
|
||||
new value of x_i is new_value
|
||||
a_ii*(x_i + new_value - x_i) + a_ij*((x_i - new_value)*a_ii/a_ij + x_j) + r_i = 0
|
||||
|
||||
Let r_k be a row where x_j has coefficient x_kj != 0.
|
||||
r_k <- r_k * a_ij - r_i * a_kj
|
||||
*/
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::update_and_pivot(var_t x_i, var_t x_j, numeral const& a_ij, eps_numeral const& new_value) {
|
||||
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;
|
||||
numeral const& a_ii = x_iI.m_base_coeff;
|
||||
em.mul(theta, a_ii, theta);
|
||||
em.div(theta, a_ij, theta);
|
||||
update_value(x_j, theta);
|
||||
SASSERT(em.eq(x_iI.m_value, new_value));
|
||||
unsigned r_i = x_iI.m_base2row;
|
||||
x_jI.m_base2row = r_i;
|
||||
m_row2base[r_i] = x_j;
|
||||
x_jI.m_is_base = true;
|
||||
x_iI.m_is_base = false;
|
||||
if (outside_bounds(x_j)) {
|
||||
m_to_patch.insert(x_j);
|
||||
}
|
||||
col_iterator it = M.col_begin(x_j), end = M.col_end(x_j);
|
||||
scoped_numeral a_kj(m), g(m);
|
||||
for (; it != end; ++it) {
|
||||
row r_k = it.get_row();
|
||||
if (r_k != r_i) {
|
||||
a_kj = it.get_row_entry().m_coeff;
|
||||
a_kj.neg();
|
||||
M.mul(r_k, a_ij);
|
||||
M.add(r_k, a_kj, r_i);
|
||||
var_t s = m_row2base[r_k.id()];
|
||||
numeral& coeff = m_vars[s].m_base_coeff;
|
||||
m.mul(coeff, a_ij, coeff);
|
||||
M.gcd_normalize(r_k, g);
|
||||
if (!m.is_one(g)) {
|
||||
m.div(coeff, g, coeff);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::update_value(var_t v, eps_numeral const& delta) {
|
||||
update_value_core(v, delta);
|
||||
col_iterator it = M.col_begin(v), end = M.col_end(v);
|
||||
|
||||
// v <- v + delta
|
||||
// s*s_coeff + v*v_coeff + R = 0
|
||||
// ->
|
||||
// (v + delta)*v_coeff + (s - delta*v_coeff/s_coeff)*v + R = 0
|
||||
for (; it != end; ++it) {
|
||||
row r = it.get_row();
|
||||
var_t s = m_row2base[r.id()];
|
||||
var_info& si = m_vars[s];
|
||||
scoped_eps_numeral delta2(em);
|
||||
numeral const& coeff = it.get_row_entry().m_coeff;
|
||||
em.mul(delta, coeff, delta2);
|
||||
em.div(delta2, si.m_base_coeff, delta2);
|
||||
delta2.neg();
|
||||
update_value_core(s, delta2);
|
||||
// TBD m.add(si.m_base_coeff, delta2, si.m_base_coeff);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void simplex<Ext>::update_value_core(var_t v, eps_numeral const& delta) {
|
||||
eps_numeral& val = m_vars[v].m_value;
|
||||
em.add(val, delta, val);
|
||||
if (is_base(v) && outside_bounds(v)) {
|
||||
m_to_patch.insert(v);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
bool simplex<Ext>::below_lower(var_t v) const {
|
||||
var_info const& vi = m_vars[v];
|
||||
return vi.m_lower_valid && em.lt(vi.m_value, vi.m_lower);
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
bool simplex<Ext>::above_upper(var_t v) const {
|
||||
var_info const& vi = m_vars[v];
|
||||
return vi.m_upper_valid && em.gt(vi.m_value, vi.m_upper);
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
bool simplex<Ext>::above_lower(var_t v) const {
|
||||
var_info const& vi = m_vars[v];
|
||||
return !vi.m_lower_valid || em.gt(vi.m_value, vi.m_lower);
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
bool simplex<Ext>::below_upper(var_t v) const {
|
||||
var_info const& vi = m_vars[v];
|
||||
return !vi.m_upper_valid || em.lt(vi.m_value, vi.m_upper);
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
bool simplex<Ext>::make_var_feasible(var_t x_i) {
|
||||
scoped_numeral a_ij(m);
|
||||
scoped_eps_numeral value(em);
|
||||
bool is_below;
|
||||
if (below_lower(x_i)) {
|
||||
is_below = true;
|
||||
value = m_vars[x_i].m_lower;
|
||||
}
|
||||
else if (above_upper(x_i)) {
|
||||
is_below = false;
|
||||
value = m_vars[x_i].m_upper;
|
||||
}
|
||||
else {
|
||||
// x_i is already feasible
|
||||
return true;
|
||||
}
|
||||
var_t x_j = select_pivot(x_i, is_below, a_ij);
|
||||
if (x_j != null_var) {
|
||||
update_and_pivot(x_i, x_j, a_ij, value);
|
||||
}
|
||||
return x_j != null_var;
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Wrapper for select_blands_pivot_core 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);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
\brief 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
|
||||
if there is no x_j. Otherwise, return x_j and store its coefficient
|
||||
in out_a_ij.
|
||||
|
||||
The argument is_below is true (false) if x_i is below its lower
|
||||
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) {
|
||||
SASSERT(is_base(x_i));
|
||||
var_t max = get_num_vars();
|
||||
var_t result = max;
|
||||
row r = row(m_vars[x_i].m_base2row);
|
||||
int n;
|
||||
unsigned best_col_sz = UINT_MAX;
|
||||
int best_so_far = INT_MAX;
|
||||
|
||||
row_iterator it = M.row_begin(r), end = M.row_end(r);
|
||||
|
||||
for (; it != end; ++it) {
|
||||
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)))) {
|
||||
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)) {
|
||||
result = x_j;
|
||||
out_a_ij = a_ij;
|
||||
best_so_far = num;
|
||||
best_col_sz = col_sz;
|
||||
n = 1;
|
||||
}
|
||||
else if (num == best_so_far && col_sz == best_col_sz) {
|
||||
n++;
|
||||
if (m_random()%n == 0) {
|
||||
result = x_j;
|
||||
out_a_ij = a_ij;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return result < max ? result : null_var;
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Return the number of base variables that are non free and are v dependent.
|
||||
The function adds 1 to the result if v is non free.
|
||||
The function returns with a partial result r if r > best_so_far.
|
||||
This function is used to select the pivot variable.
|
||||
*/
|
||||
template<typename Ext>
|
||||
int simplex<Ext>::get_num_non_free_dep_vars(var_t x_j, int best_so_far) {
|
||||
int result = is_non_free(x_j);
|
||||
col_iterator it = M.col_begin(x_j), end = M.col_end(x_j);
|
||||
for (; it != end; ++it) {
|
||||
var_t s = m_row2base[it.get_row().id()];
|
||||
result += is_non_free(s);
|
||||
if (result > best_so_far)
|
||||
return result;
|
||||
}
|
||||
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
|
||||
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) {
|
||||
SASSERT(is_base(x_i));
|
||||
unsigned max = get_num_vars();
|
||||
var_t result = max;
|
||||
row r(m_vars[x_i].m_base2row);
|
||||
row_iterator it = M.row_begin(r), end = M.row_end(r);
|
||||
for (; it != end; ++it) {
|
||||
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)))) {
|
||||
SASSERT(!is_base(x_j));
|
||||
if (x_j < result) {
|
||||
result = x_j;
|
||||
out_a_ij = a_ij;
|
||||
}
|
||||
}
|
||||
}
|
||||
return result < max ? result : null_var;
|
||||
}
|
||||
|
||||
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.
|
||||
return l_true;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
typename simplex<Ext>::pivot_strategy_t
|
||||
simplex<Ext>::pivot_strategy() {
|
||||
if (m_bland) {
|
||||
return S_BLAND;
|
||||
}
|
||||
return S_DEFAULT;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
typename simplex<Ext>::var_t
|
||||
simplex<Ext>::select_var_to_fix() {
|
||||
switch (pivot_strategy()) {
|
||||
case S_BLAND:
|
||||
return select_smallest_var();
|
||||
case S_GREATEST_ERROR:
|
||||
return select_error_var(false);
|
||||
case S_LEAST_ERROR:
|
||||
return select_error_var(true);
|
||||
default:
|
||||
return select_smallest_var();
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
typename simplex<Ext>::var_t
|
||||
simplex<Ext>::select_error_var(bool least) {
|
||||
var_t best = null_var;
|
||||
scoped_eps_numeral best_error(em);
|
||||
scoped_eps_numeral curr_error(em);
|
||||
typename var_heap::iterator it = m_to_patch.begin();
|
||||
typename var_heap::iterator end = m_to_patch.end();
|
||||
for (; it != end; ++it) {
|
||||
var_t v = *it;
|
||||
var_info const& vi = m_vars[v];
|
||||
if (below_lower(v))
|
||||
em.sub(vi.m_lower, vi.m_value, curr_error);
|
||||
else if (above_upper(v))
|
||||
em.sub(vi.m_value, vi.m_lower, curr_error);
|
||||
else
|
||||
continue;
|
||||
SASSERT(is_pos(curr_error));
|
||||
if ((best == null_var) ||
|
||||
(!least && curr_error > best_error) ||
|
||||
(least && curr_error < best_error)) {
|
||||
best = v;
|
||||
best_error = curr_error;
|
||||
}
|
||||
}
|
||||
if (best == null_var)
|
||||
m_to_patch.clear(); // all variables are satisfied
|
||||
else
|
||||
m_to_patch.erase(best);
|
||||
return best;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
bool simplex<Ext>::well_formed() const {
|
||||
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.
|
||||
//
|
||||
// TBD: extract coefficient of base variable and compare
|
||||
// with m_vars[s].m_base_coeff;
|
||||
//
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
254
src/math/simplex/sparse_matrix.h
Normal file
254
src/math/simplex/sparse_matrix.h
Normal file
|
|
@ -0,0 +1,254 @@
|
|||
/*++
|
||||
Copyright (c) 2014 Microsoft Corporation
|
||||
|
||||
Module Name:
|
||||
|
||||
sparse_matrix.h
|
||||
|
||||
Abstract:
|
||||
|
||||
|
||||
Author:
|
||||
|
||||
Nikolaj Bjorner (nbjorner) 2014-01-15
|
||||
|
||||
Notes:
|
||||
|
||||
--*/
|
||||
|
||||
#ifndef _SPARSE_MATRIX_H_
|
||||
#define _SPARSE_MATRIX_H_
|
||||
|
||||
#include "mpq_inf.h"
|
||||
|
||||
namespace simplex {
|
||||
|
||||
template<typename Ext>
|
||||
class sparse_matrix {
|
||||
public:
|
||||
typedef typename Ext::numeral numeral;
|
||||
typedef typename Ext::scoped_numeral scoped_numeral;
|
||||
typedef typename Ext::manager manager;
|
||||
typedef unsigned var_t;
|
||||
|
||||
struct row_entry {
|
||||
numeral m_coeff;
|
||||
var_t m_var;
|
||||
row_entry(numeral& c, var_t v): m_coeff(c), m_var(v) {}
|
||||
};
|
||||
|
||||
private:
|
||||
|
||||
static const int dead_id = -1;
|
||||
|
||||
/**
|
||||
\brief A row_entry is: m_var*m_coeff
|
||||
|
||||
m_col_idx points to the place in the
|
||||
column where the variable occurs.
|
||||
*/
|
||||
struct _row_entry : public row_entry {
|
||||
union {
|
||||
int m_col_idx;
|
||||
int m_next_free_row_entry_idx;
|
||||
};
|
||||
_row_entry(numeral const & c, var_t v): row_entry(c, v), m_col_idx(0) {}
|
||||
_row_entry() : row_entry(numeral(), dead_id), m_col_idx(0) {}
|
||||
bool is_dead() const { return m_var == dead_id; }
|
||||
};
|
||||
|
||||
/**
|
||||
\brief A column entry points to the row and the row_entry within the row
|
||||
that has a non-zero coefficient on the variable associated
|
||||
with the column entry.
|
||||
*/
|
||||
struct col_entry {
|
||||
int m_row_id;
|
||||
union {
|
||||
int m_row_idx;
|
||||
int m_next_free_col_entry_idx;
|
||||
};
|
||||
col_entry(int r, int i): m_row_id(r), m_row_idx(i) {}
|
||||
col_entry(): m_row_id(0), m_row_idx(0) {}
|
||||
bool is_dead() const { return m_row_id == dead_id; }
|
||||
};
|
||||
|
||||
struct column;
|
||||
|
||||
/**
|
||||
\brief A row contains a base variable and set of
|
||||
row_entries. The base variable must occur in the set of
|
||||
row_entries with coefficient 1.
|
||||
*/
|
||||
struct _row {
|
||||
vector<_row_entry> m_entries;
|
||||
unsigned m_size; // the real size, m_entries contains dead row_entries.
|
||||
int m_first_free_idx; // first available position.
|
||||
_row();
|
||||
unsigned size() const { return m_size; }
|
||||
unsigned num_entries() const { return m_entries.size(); }
|
||||
void reset(manager& m);
|
||||
_row_entry & add_row_entry(unsigned & pos_idx);
|
||||
void del_row_entry(unsigned idx);
|
||||
void compress(manager& m, vector<column> & cols);
|
||||
void compress_if_needed(manager& m, vector<column> & cols);
|
||||
void save_var_pos(svector<int> & result_map) const;
|
||||
void reset_var_pos(svector<int> & result_map) const;
|
||||
bool is_coeff_of(var_t v, numeral const & expected) const;
|
||||
int get_idx_of(var_t v) const;
|
||||
};
|
||||
|
||||
/**
|
||||
\brief A column stores in which rows a variable occurs.
|
||||
The column may have free/dead entries. The field m_first_free_idx
|
||||
is a reference to the first free/dead entry.
|
||||
*/
|
||||
struct column {
|
||||
svector<col_entry> m_entries;
|
||||
unsigned m_size;
|
||||
int m_first_free_idx;
|
||||
|
||||
column():m_size(0), m_first_free_idx(-1) {}
|
||||
unsigned size() const { return m_size; }
|
||||
unsigned num_entries() const { return m_entries.size(); }
|
||||
void reset();
|
||||
void compress(vector<_row> & rows);
|
||||
void compress_if_needed(vector<_row> & rows);
|
||||
void compress_singleton(vector<_row> & rows, unsigned singleton_pos);
|
||||
col_entry const * get_first_col_entry() const;
|
||||
col_entry & add_col_entry(int & pos_idx);
|
||||
void del_col_entry(unsigned idx);
|
||||
};
|
||||
|
||||
manager m;
|
||||
vector<_row> m_rows;
|
||||
svector<unsigned> m_dead_rows; // rows to recycle
|
||||
vector<column> m_columns; // per var
|
||||
svector<int> m_var_pos; // temporary map from variables to positions in row
|
||||
|
||||
bool well_formed_row(unsigned row_id) const;
|
||||
bool well_formed_column(unsigned column_id) const;
|
||||
void del_row_entry(_row& r, unsigned pos);
|
||||
|
||||
public:
|
||||
|
||||
sparse_matrix() {}
|
||||
~sparse_matrix();
|
||||
|
||||
class row {
|
||||
unsigned m_id;
|
||||
public:
|
||||
row(unsigned r):m_id(r) {}
|
||||
row():m_id(UINT_MAX) {}
|
||||
bool operator!=(row const& other) const {
|
||||
return m_id != other.m_id;
|
||||
}
|
||||
unsigned id() const { return m_id; }
|
||||
};
|
||||
|
||||
void ensure_var(var_t v);
|
||||
|
||||
row mk_row();
|
||||
void add(row r, numeral const& n, var_t var);
|
||||
void add(row r, numeral const& n, row src);
|
||||
void mul(row r, numeral const& n);
|
||||
void neg(row r);
|
||||
void del(row r);
|
||||
|
||||
void gcd_normalize(row const& r, scoped_numeral& g);
|
||||
|
||||
class row_iterator {
|
||||
friend sparse_matrix;
|
||||
unsigned m_curr;
|
||||
_row & m_row;
|
||||
void move_to_used() {
|
||||
while (m_curr < m_row.m_entries.size() &&
|
||||
m_row.m_entries[m_curr].is_dead()) {
|
||||
++m_curr;
|
||||
}
|
||||
}
|
||||
row_iterator(_row & r, bool begin):
|
||||
m_curr(0), m_row(r) {
|
||||
if (begin) {
|
||||
move_to_used();
|
||||
}
|
||||
else {
|
||||
m_curr = m_row.m_entries.size();
|
||||
}
|
||||
}
|
||||
public:
|
||||
row_entry & operator*() const { return m_row.m_entries[m_curr]; }
|
||||
row_entry * operator->() const { return &(operator*()); }
|
||||
row_iterator & operator++() { ++m_curr; move_to_used(); return *this; }
|
||||
row_iterator operator++(int) { row_iterator tmp = *this; ++*this; return tmp; }
|
||||
bool operator==(row_iterator const & it) const { return m_curr == it.m_curr; }
|
||||
bool operator!=(row_iterator const & it) const { return m_curr != it.m_curr; }
|
||||
};
|
||||
|
||||
row_iterator row_begin(row const& r) { return row_iterator(m_rows[r.id()], true); }
|
||||
row_iterator row_end(row const& r) { return row_iterator(m_rows[r.id()], false); }
|
||||
|
||||
unsigned column_size(var_t v) const { return m_columns[v].size(); }
|
||||
|
||||
class col_iterator {
|
||||
friend sparse_matrix;
|
||||
unsigned m_curr;
|
||||
column const& m_col;
|
||||
vector<_row> const& m_rows;
|
||||
void move_to_used() {
|
||||
while (m_curr < m_col.m_entries.size() && m_col.m_entries[m_curr].is_dead()) {
|
||||
++m_curr;
|
||||
}
|
||||
}
|
||||
col_iterator(column const& c, vector<_row> const& r, bool begin):
|
||||
m_curr(0), m_col(c), m_rows(r) {
|
||||
if (begin) {
|
||||
move_to_used();
|
||||
}
|
||||
else {
|
||||
m_curr = m_col.m_entries.size();
|
||||
}
|
||||
}
|
||||
public:
|
||||
row get_row() { return row(m_col.m_entries[m_curr].m_row_id); }
|
||||
row_entry const& get_row_entry() {
|
||||
col_entry const& c = m_col.m_entries[m_curr];
|
||||
int row_id = c.m_row_id;
|
||||
return m_rows[row_id].m_entries[c.m_row_idx];
|
||||
}
|
||||
|
||||
col_iterator & operator++() { ++m_curr; move_to_used(); return *this; }
|
||||
col_iterator operator++(int) { col_iterator tmp = *this; ++*this; return tmp; }
|
||||
bool operator==(col_iterator const & it) const { return m_curr == it.m_curr; }
|
||||
bool operator!=(col_iterator const & it) const { return m_curr != it.m_curr; }
|
||||
};
|
||||
|
||||
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;
|
||||
bool well_formed() const;
|
||||
|
||||
|
||||
};
|
||||
|
||||
struct mpz_ext {
|
||||
typedef mpz numeral;
|
||||
typedef scoped_mpz scoped_numeral;
|
||||
typedef unsynch_mpz_manager manager;
|
||||
typedef mpq_inf eps_numeral;
|
||||
typedef unsynch_mpq_inf_manager eps_manager;
|
||||
};
|
||||
|
||||
struct mpq_ext {
|
||||
typedef mpq numeral;
|
||||
typedef scoped_mpq scoped_numeral;
|
||||
typedef unsynch_mpq_manager manager;
|
||||
typedef mpq_inf eps_numeral;
|
||||
typedef unsynch_mpq_inf_manager eps_manager;
|
||||
};
|
||||
|
||||
};
|
||||
|
||||
|
||||
#endif
|
||||
580
src/math/simplex/sparse_matrix_def.h
Normal file
580
src/math/simplex/sparse_matrix_def.h
Normal file
|
|
@ -0,0 +1,580 @@
|
|||
/*++
|
||||
Copyright (c) 2014 Microsoft Corporation
|
||||
|
||||
Module Name:
|
||||
|
||||
sparse_matrix_def.h
|
||||
|
||||
Abstract:
|
||||
|
||||
|
||||
Author:
|
||||
|
||||
Nikolaj Bjorner (nbjorner) 2014-01-15
|
||||
|
||||
Notes:
|
||||
|
||||
mainly hoisted from theory_arith.h and theory_arith_aux.h
|
||||
|
||||
--*/
|
||||
|
||||
#ifndef _SPARSE_MATRIX_DEF_H_
|
||||
#define _SPARSE_MATRIX_DEF_H_
|
||||
|
||||
#include "sparse_matrix.h"
|
||||
#include "uint_set.h"
|
||||
|
||||
namespace simplex {
|
||||
|
||||
// -----------------------------------
|
||||
//
|
||||
// Rows
|
||||
//
|
||||
// -----------------------------------
|
||||
|
||||
template<typename Ext>
|
||||
sparse_matrix<Ext>::_row::_row():
|
||||
m_size(0),
|
||||
m_first_free_idx(-1) {
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::_row::reset(manager& m) {
|
||||
for (unsigned i = 0; i < m_entries.size(); ++i) {
|
||||
m.reset(m_entries[i].m_coeff);
|
||||
}
|
||||
m_entries.reset();
|
||||
m_size = 0;
|
||||
m_first_free_idx = -1;
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Add a new row_entry. The result is a reference to the new row_entry.
|
||||
The position of the new row_entry in the
|
||||
row is stored in pos_idx.
|
||||
*/
|
||||
template<typename Ext>
|
||||
typename sparse_matrix<Ext>::_row_entry &
|
||||
sparse_matrix<Ext>::_row::add_row_entry(unsigned & pos_idx) {
|
||||
m_size++;
|
||||
if (m_first_free_idx == -1) {
|
||||
pos_idx = m_entries.size();
|
||||
m_entries.push_back(_row_entry());
|
||||
return m_entries.back();
|
||||
}
|
||||
else {
|
||||
SASSERT(m_first_free_idx >= 0);
|
||||
pos_idx = static_cast<unsigned>(m_first_free_idx);
|
||||
_row_entry & result = m_entries[pos_idx];
|
||||
SASSERT(result.is_dead());
|
||||
m_first_free_idx = result.m_next_free_row_entry_idx;
|
||||
return result;
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Delete row_entry at position idx.
|
||||
*/
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::_row::del_row_entry(unsigned idx) {
|
||||
_row_entry & t = m_entries[idx];
|
||||
SASSERT(!t.is_dead());
|
||||
t.m_next_free_row_entry_idx = m_first_free_idx;
|
||||
t.m_var = dead_id;
|
||||
m_size--;
|
||||
m_first_free_idx = idx;
|
||||
SASSERT(t.is_dead());
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Remove holes (i.e., dead entries) from the row.
|
||||
*/
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::_row::compress(manager& m, vector<column> & cols) {
|
||||
unsigned i = 0;
|
||||
unsigned j = 0;
|
||||
unsigned sz = m_entries.size();
|
||||
for (; i < sz; i++) {
|
||||
_row_entry & t1 = m_entries[i];
|
||||
if (!t1.is_dead()) {
|
||||
if (i != j) {
|
||||
_row_entry & t2 = m_entries[j];
|
||||
t2.m_coeff.swap(t1.m_coeff);
|
||||
t2.m_var = t1.m_var;
|
||||
t2.m_col_idx = t1.m_col_idx;
|
||||
SASSERT(!t2.is_dead());
|
||||
column & col = cols[t2.m_var];
|
||||
col.m_entries[t2.m_col_idx].m_row_idx = j;
|
||||
}
|
||||
j++;
|
||||
}
|
||||
}
|
||||
SASSERT(j == m_size);
|
||||
//
|
||||
// alternative: update the free-list to point to the
|
||||
// tail and avoid shrinking.
|
||||
// if m.does not allocate memory (for wrapper around
|
||||
// double), also bypass this step.
|
||||
//
|
||||
for (unsigned i = m_size; i < m_entries.size(); ++i) {
|
||||
m.reset(m_entries[i].m_coeff);
|
||||
}
|
||||
m_entries.shrink(m_size);
|
||||
m_first_free_idx = -1;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::_row::compress_if_needed(manager& m, vector<column> & cols) {
|
||||
if (size() *2 < num_entries()) {
|
||||
compress(m, cols);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Fill the map var -> pos/idx
|
||||
*/
|
||||
template<typename Ext>
|
||||
inline void sparse_matrix<Ext>::_row::save_var_pos(svector<int> & result_map) const {
|
||||
typename vector<_row_entry>::const_iterator it = m_entries.begin();
|
||||
typename vector<_row_entry>::const_iterator end = m_entries.end();
|
||||
unsigned idx = 0;
|
||||
for (; it != end; ++it, ++idx) {
|
||||
if (!it->is_dead()) {
|
||||
result_map[it->m_var] = idx;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Reset the map var -> pos/idx. That is for all variables v in the row, set result[v] = -1
|
||||
This method can be viewed as the "inverse" of save_var_pos.
|
||||
*/
|
||||
template<typename Ext>
|
||||
inline void sparse_matrix<Ext>::_row::reset_var_pos(svector<int> & result_map) const {
|
||||
typename vector<_row_entry>::const_iterator it = m_entries.begin();
|
||||
typename vector<_row_entry>::const_iterator end = m_entries.end();
|
||||
for (; it != end; ++it) {
|
||||
if (!it->is_dead()) {
|
||||
result_map[it->m_var] = -1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
int sparse_matrix<Ext>::_row::get_idx_of(var_t v) const {
|
||||
typename vector<_row_entry>::const_iterator it = m_entries.begin();
|
||||
typename vector<_row_entry>::const_iterator end = m_entries.end();
|
||||
for (unsigned idx = 0; it != end; ++it, ++idx) {
|
||||
if (!it->is_dead() && it->m_var == v)
|
||||
return idx;
|
||||
}
|
||||
return -1;
|
||||
}
|
||||
|
||||
// -----------------------------------
|
||||
//
|
||||
// Columns
|
||||
//
|
||||
// -----------------------------------
|
||||
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::column::reset() {
|
||||
m_entries.reset();
|
||||
m_size = 0;
|
||||
m_first_free_idx = -1;
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Remove holes (i.e., dead entries) from the column.
|
||||
*/
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::column::compress(vector<_row> & rows) {
|
||||
unsigned i = 0;
|
||||
unsigned j = 0;
|
||||
unsigned sz = m_entries.size();
|
||||
for (; i < sz; i++) {
|
||||
col_entry & e1 = m_entries[i];
|
||||
if (!e1.is_dead()) {
|
||||
if (i != j) {
|
||||
m_entries[j] = e1;
|
||||
_row & r = rows[e1.m_row_id];
|
||||
r.m_entries[e1.m_row_idx].m_col_idx = j;
|
||||
}
|
||||
j++;
|
||||
}
|
||||
}
|
||||
SASSERT(j == m_size);
|
||||
m_entries.shrink(m_size);
|
||||
m_first_free_idx = -1;
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Invoke compress if the column contains too many holes (i.e., dead entries).
|
||||
*/
|
||||
template<typename Ext>
|
||||
inline void sparse_matrix<Ext>::column::compress_if_needed(vector<_row> & rows) {
|
||||
if (size() * 2 < num_entries()) {
|
||||
compress(rows);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Special version of compress, that is used when the column contain
|
||||
only one entry located at position singleton_pos.
|
||||
*/
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::column::compress_singleton(vector<_row> & rows, unsigned singleton_pos) {
|
||||
SASSERT(m_size == 1);
|
||||
if (singleton_pos != 0) {
|
||||
col_entry & s = m_entries[singleton_pos];
|
||||
m_entries[0] = s;
|
||||
row & r = rows[s.m_row_id];
|
||||
r[s.m_row_idx].m_col_idx = 0;
|
||||
}
|
||||
m_first_free_idx = -1;
|
||||
m_entries.shrink(1);
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
const typename sparse_matrix<Ext>::col_entry *
|
||||
sparse_matrix<Ext>::column::get_first_col_entry() const {
|
||||
typename svector<col_entry>::const_iterator it = m_entries.begin();
|
||||
typename svector<col_entry>::const_iterator end = m_entries.end();
|
||||
for (; it != end; ++it) {
|
||||
if (!it->is_dead()) {
|
||||
return it;
|
||||
}
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
typename sparse_matrix<Ext>::col_entry &
|
||||
sparse_matrix<Ext>::column::add_col_entry(int & pos_idx) {
|
||||
m_size++;
|
||||
if (m_first_free_idx == -1) {
|
||||
pos_idx = m_entries.size();
|
||||
m_entries.push_back(col_entry());
|
||||
return m_entries.back();
|
||||
}
|
||||
else {
|
||||
pos_idx = m_first_free_idx;
|
||||
col_entry & result = m_entries[pos_idx];
|
||||
SASSERT(result.is_dead());
|
||||
m_first_free_idx = result.m_next_free_col_entry_idx;
|
||||
return result;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::column::del_col_entry(unsigned idx) {
|
||||
col_entry & c = m_entries[idx];
|
||||
SASSERT(!c.is_dead());
|
||||
c.m_row_id = dead_id;
|
||||
c.m_next_free_col_entry_idx = m_first_free_idx;
|
||||
m_first_free_idx = idx;
|
||||
m_size--;
|
||||
}
|
||||
|
||||
// -----------------------------------
|
||||
//
|
||||
// Matrix
|
||||
//
|
||||
// -----------------------------------
|
||||
|
||||
template<typename Ext>
|
||||
sparse_matrix<Ext>::~sparse_matrix() {
|
||||
for (unsigned i = 0; i < m_rows.size(); ++i) {
|
||||
_row& r = m_rows[i];
|
||||
for (unsigned j = 0; j < r.m_entries.size(); ++j) {
|
||||
m.reset(r.m_entries[j].m_coeff);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::ensure_var(var_t v) {
|
||||
while (m_columns.size() <= v) {
|
||||
m_columns.push_back(column());
|
||||
m_var_pos.push_back(-1);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
typename sparse_matrix<Ext>::row
|
||||
sparse_matrix<Ext>::mk_row() {
|
||||
if (m_dead_rows.empty()) {
|
||||
row r(m_rows.size());
|
||||
m_rows.push_back(_row());
|
||||
return r;
|
||||
}
|
||||
else {
|
||||
row r(m_dead_rows.back());
|
||||
m_dead_rows.pop_back();
|
||||
return r;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::add(row dst, numeral const& n, var_t v) {
|
||||
_row& r = m_rows[dst.id()];
|
||||
column& c = m_columns[v];
|
||||
unsigned r_idx;
|
||||
int c_idx;
|
||||
_row_entry & r_entry = r.add_row_entry(r_idx);
|
||||
col_entry& c_entry = c.add_col_entry(c_idx);
|
||||
r_entry.m_var = v;
|
||||
m.set(r_entry.m_coeff, n);
|
||||
r_entry.m_col_idx = c_idx;
|
||||
c_entry.m_row_id = dst.id();
|
||||
c_entry.m_row_idx = r_idx;
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Set row1 <- row1 + row2 * n
|
||||
*/
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::add(row row1, numeral const& n, row row2) {
|
||||
// m_stats.m_add_rows++;
|
||||
_row & r1 = m_rows[row1.id()];
|
||||
_row & r2 = m_rows[row2.id()];
|
||||
|
||||
r1.save_var_pos(m_var_pos);
|
||||
|
||||
//
|
||||
// loop over variables in row2,
|
||||
// add terms in row2 to row1.
|
||||
//
|
||||
|
||||
#define ADD_ROW(_SET_COEFF_, _ADD_COEFF_) \
|
||||
row_iterator it = row_begin(row2); \
|
||||
row_iterator end = row_end(row2); \
|
||||
for (; it != end; ++it) { \
|
||||
var_t v = it->m_var; \
|
||||
int pos = m_var_pos[v]; \
|
||||
if (pos == -1) { \
|
||||
/* variable v is not in row1 */ \
|
||||
unsigned row_idx; \
|
||||
_row_entry & r_entry = r1.add_row_entry(row_idx); \
|
||||
r_entry.m_var = v; \
|
||||
m.set(r_entry.m_coeff, it->m_coeff); \
|
||||
_SET_COEFF_; \
|
||||
column & c = m_columns[v]; \
|
||||
int col_idx; \
|
||||
col_entry & c_entry = c.add_col_entry(col_idx); \
|
||||
r_entry.m_col_idx = col_idx; \
|
||||
c_entry.m_row_id = row1.id(); \
|
||||
c_entry.m_row_idx = row_idx; \
|
||||
} \
|
||||
else { \
|
||||
/* variable v is in row1 */ \
|
||||
_row_entry & r_entry = r1.m_entries[pos]; \
|
||||
SASSERT(r_entry.m_var == v); \
|
||||
_ADD_COEFF_; \
|
||||
if (m.is_zero(r_entry.m_coeff)) { \
|
||||
del_row_entry(r1, pos); \
|
||||
} \
|
||||
} \
|
||||
} \
|
||||
|
||||
((void) 0)
|
||||
|
||||
if (m.is_one(n)) {
|
||||
ADD_ROW({},
|
||||
m.add(r_entry.m_coeff, it->m_coeff, r_entry.m_coeff));
|
||||
}
|
||||
else if (m.is_minus_one(n)) {
|
||||
ADD_ROW(m.neg(r_entry.m_coeff),
|
||||
m.sub(r_entry.m_coeff, it->m_coeff, r_entry.m_coeff));
|
||||
}
|
||||
else {
|
||||
scoped_numeral tmp(m);
|
||||
ADD_ROW(m.mul(r_entry.m_coeff, n, r_entry.m_coeff),
|
||||
m.mul(it->m_coeff, n, tmp);
|
||||
m.add(r_entry.m_coeff, tmp, r_entry.m_coeff));
|
||||
}
|
||||
|
||||
r1.reset_var_pos(m_var_pos);
|
||||
r1.compress_if_needed(m, m_columns);
|
||||
}
|
||||
|
||||
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::del_row_entry(_row& r, unsigned pos) {
|
||||
_row_entry & r_entry = r.m_entries[pos];
|
||||
var_t v = r_entry.m_var;
|
||||
int col_idx = r_entry.m_col_idx;
|
||||
r.del_row_entry(pos);
|
||||
column & c = m_columns[v];
|
||||
c.del_col_entry(col_idx);
|
||||
c.compress_if_needed(m_rows);
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Set row <- -row
|
||||
*/
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::neg(row r) {
|
||||
row_iterator it = row_begin(r);
|
||||
row_iterator end = row_end(r);
|
||||
for (; it != end; ++it) {
|
||||
m.neg(it->m_coeff);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Set row <- n*row
|
||||
*/
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::mul(row r, numeral const& n) {
|
||||
SASSERT(!m.is_zero(n));
|
||||
if (m.is_one(n)) {
|
||||
// no op
|
||||
}
|
||||
else if (m.is_minus_one(n)) {
|
||||
neg(r);
|
||||
}
|
||||
else {
|
||||
row_iterator it = row_begin(r);
|
||||
row_iterator end = row_end(r);
|
||||
for (; it != end; ++it) {
|
||||
m.mul(it->m_coeff, n, it->m_coeff);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
\brief Delete row.
|
||||
*/
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::del(row r) {
|
||||
_row& rw = r.id();
|
||||
for (unsigned i = 0; i < rw.m_entries.size(); ++i) {
|
||||
_row_entry& e = rw.m_entries[i];
|
||||
if (!e.is_dead()) {
|
||||
del_row_entry(rw, i);
|
||||
}
|
||||
}
|
||||
SASSERT(rw.m_size == 0);
|
||||
m_dead_rows.push_back(r.id());
|
||||
}
|
||||
|
||||
/**
|
||||
\brief normalize coefficients by dividing with they coefficients.
|
||||
return the gcd.
|
||||
*/
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::gcd_normalize(row const& r, scoped_numeral& g) {
|
||||
g.reset();
|
||||
row_iterator it = row_begin(r), end = row_end(r);
|
||||
for (; it != end && !m.is_one(g); ++it) {
|
||||
if (m.is_zero(g)) g = it->m_coeff;
|
||||
else m.gcd(g, it->m_coeff, g);
|
||||
}
|
||||
if (m.is_zero(g)) {
|
||||
g = numeral(1);
|
||||
}
|
||||
if (!m.is_one(g)) {
|
||||
row_iterator it2 = row_begin(r);
|
||||
for (; it2 != end; ++it2) {
|
||||
m.div(it2->m_coeff, g, it2->m_coeff);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
\brief well_formed check
|
||||
*/
|
||||
template<typename Ext>
|
||||
bool sparse_matrix<Ext>::well_formed() const {
|
||||
for (unsigned i = 0; i < m_rows.size(); ++i) {
|
||||
well_formed_row(i);
|
||||
}
|
||||
for (unsigned i = 0; i < m_columns.size(); ++i) {
|
||||
well_formed_column(i);
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
/**
|
||||
\brief well_formed row check
|
||||
*/
|
||||
template<typename Ext>
|
||||
bool sparse_matrix<Ext>::well_formed_row(unsigned row_id) const {
|
||||
uint_set vars, dead;
|
||||
_row const& r = m_rows[row_id];
|
||||
for (unsigned i = 0; i < r.num_entries(); ++i) {
|
||||
_row_entry const& e = r.m_entries[i];
|
||||
if (e.is_dead()) {
|
||||
dead.insert(i);
|
||||
continue;
|
||||
}
|
||||
SASSERT(!vars.contains(e.m_var));
|
||||
SASSERT(!m.is_zero(e.m_coeff));
|
||||
SASSERT(e.m_var != dead_id);
|
||||
col_entry const& c = m_columns[e.m_var].m_entries[e.m_col_idx];
|
||||
SASSERT(c.m_row_id == row_id);
|
||||
SASSERT(c.m_row_idx == i);
|
||||
vars.insert(e.m_var);
|
||||
}
|
||||
int idx = r.m_first_free_idx;
|
||||
while (idx != -1) {
|
||||
SASSERT(dead.contains(idx));
|
||||
dead.remove(idx);
|
||||
idx = r.m_entries[idx].m_next_free_row_entry_idx;
|
||||
}
|
||||
SASSERT(dead.empty());
|
||||
return true;
|
||||
}
|
||||
|
||||
/**
|
||||
\brief well_formed column check
|
||||
*/
|
||||
template<typename Ext>
|
||||
bool sparse_matrix<Ext>::well_formed_column(var_t v) const {
|
||||
uint_set rows, dead;
|
||||
column const& col = m_columns[v];
|
||||
for (unsigned i = 0; i < col.num_entries(); ++i) {
|
||||
col_entry const& c = col.m_entries[i];
|
||||
if (c.is_dead()) {
|
||||
dead.insert(i);
|
||||
continue;
|
||||
}
|
||||
SASSERT(!rows.contains(c.m_row_id));
|
||||
_row const& row = m_rows[c.m_row_id];
|
||||
_row_entry const& r = row.m_entries[c.m_row_idx];
|
||||
SASSERT(r.m_var == v);
|
||||
SASSERT(r.m_col_idx == i);
|
||||
rows.insert(c.m_row_id);
|
||||
}
|
||||
int idx = col.m_first_free_idx;
|
||||
while (idx != -1) {
|
||||
SASSERT(dead.contains(idx));
|
||||
dead.remove(idx);
|
||||
idx = col.m_entries[idx].m_next_free_col_entry_idx;
|
||||
}
|
||||
SASSERT(dead.empty());
|
||||
return true;
|
||||
}
|
||||
|
||||
/**
|
||||
\brief display method
|
||||
*/
|
||||
template<typename Ext>
|
||||
void sparse_matrix<Ext>::display(std::ostream& out) const {
|
||||
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)
|
||||
out << "*v" << it->m_var << " ";
|
||||
}
|
||||
out << "\n";
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
#endif
|
||||
24
src/test/simplex.cpp
Normal file
24
src/test/simplex.cpp
Normal file
|
|
@ -0,0 +1,24 @@
|
|||
#include "sparse_matrix.h"
|
||||
#include "sparse_matrix_def.h"
|
||||
#include "simplex.h"
|
||||
#include "simplex_def.h"
|
||||
|
||||
typedef simplex::simplex<simplex::mpz_ext> Simplex;
|
||||
|
||||
void tst_simplex() {
|
||||
simplex::sparse_matrix<simplex::mpz_ext> M;
|
||||
Simplex S;
|
||||
|
||||
S.make_feasible();
|
||||
|
||||
unsynch_mpz_manager m;
|
||||
scoped_mpz_vector coeffs(m);
|
||||
svector<unsigned> vars;
|
||||
for (unsigned i = 0; i < 5; ++i) {
|
||||
S.ensure_var(i);
|
||||
vars.push_back(i);
|
||||
coeffs.push_back(mpz(i+1));
|
||||
}
|
||||
|
||||
Simplex::row r = S.add_row(1, coeffs.size(), vars.c_ptr(), coeffs.c_ptr());
|
||||
}
|
||||
|
|
@ -19,7 +19,7 @@ Revision History:
|
|||
#include"mpq_inf.h"
|
||||
|
||||
template<bool SYNCH>
|
||||
std::string mpq_inf_manager<SYNCH>::to_string(mpq_inf const & a) const {
|
||||
std::string mpq_inf_manager<SYNCH>::to_string(mpq_inf const & a) {
|
||||
if (m.is_zero(a.second))
|
||||
return m.to_string(a.first);
|
||||
|
||||
|
|
|
|||
|
|
@ -26,12 +26,12 @@ typedef std::pair<mpq, mpq> mpq_inf;
|
|||
|
||||
template<bool SYNCH = true>
|
||||
class mpq_inf_manager {
|
||||
mpq_manager<SYNCH> & m;
|
||||
mpq_manager<SYNCH> m;
|
||||
double m_inf;
|
||||
public:
|
||||
typedef mpq_inf numeral;
|
||||
|
||||
mpq_inf_manager(mpq_manager<SYNCH> & _m, double inf = 0.0001):m(_m) {
|
||||
mpq_inf_manager(double inf = 0.0001) {
|
||||
set_inf(inf);
|
||||
}
|
||||
|
||||
|
|
@ -83,6 +83,14 @@ public:
|
|||
}
|
||||
|
||||
bool is_int(mpq_inf const & a) const { return m.is_int(a.first) && m.is_zero(a.second); }
|
||||
|
||||
bool is_pos(mpq_inf const & a) const {
|
||||
return m.is_pos(a.first) || (m.is_zero(a.first) && m.is_pos(a.second));
|
||||
}
|
||||
|
||||
bool is_neg(mpq_inf const & a) const {
|
||||
return m.is_neg(a.first) || (m.is_zero(a.first) && m.is_neg(a.second));
|
||||
}
|
||||
|
||||
bool is_rational(mpq_inf const & a) const { return m.is_zero(a.second); }
|
||||
|
||||
|
|
@ -104,15 +112,15 @@ public:
|
|||
return m.is_zero(a.first) && m.is_zero(a.second);
|
||||
}
|
||||
|
||||
bool eq(mpq_inf const & a, mpq_inf const & b) const {
|
||||
bool eq(mpq_inf const & a, mpq_inf const & b) {
|
||||
return m.eq(a.first, b.first) && m.eq(a.second, b.second);
|
||||
}
|
||||
|
||||
bool eq(mpq_inf const & a, mpq const & b) const {
|
||||
bool eq(mpq_inf const & a, mpq const & b) {
|
||||
return m.eq(a.first, b) && m.is_zero(a.second);
|
||||
}
|
||||
|
||||
bool eq(mpq_inf const & a, mpq const & b, inf_kind k) const {
|
||||
bool eq(mpq_inf const & a, mpq const & b, inf_kind k) {
|
||||
if (!m.eq(a.first, b))
|
||||
return false;
|
||||
switch (k) {
|
||||
|
|
@ -124,15 +132,15 @@ public:
|
|||
return false;
|
||||
}
|
||||
|
||||
bool lt(mpq_inf const & a, mpq_inf const & b) const {
|
||||
bool lt(mpq_inf const & a, mpq_inf const & b) {
|
||||
return m.lt(a.first, b.first) || (m.lt(a.second, b.second) && m.eq(a.first, b.first));
|
||||
}
|
||||
|
||||
bool lt(mpq_inf const & a, mpq const & b) const {
|
||||
bool lt(mpq_inf const & a, mpq const & b) {
|
||||
return m.lt(a.first, b) || (m.is_neg(a.second) && m.eq(a.first, b));
|
||||
}
|
||||
|
||||
bool lt(mpq_inf const & a, mpq const & b, inf_kind k) const {
|
||||
bool lt(mpq_inf const & a, mpq const & b, inf_kind k) {
|
||||
if (m.lt(a.first, b))
|
||||
return true;
|
||||
if (m.eq(a.first, b)) {
|
||||
|
|
@ -146,13 +154,13 @@ public:
|
|||
return false;
|
||||
}
|
||||
|
||||
bool gt(mpq_inf const & a, mpq_inf const & b) const { return lt(b, a); }
|
||||
bool gt(mpq_inf const & a, mpq_inf const & b) { return lt(b, a); }
|
||||
|
||||
bool gt(mpq_inf const & a, mpq const & b) const {
|
||||
bool gt(mpq_inf const & a, mpq const & b) {
|
||||
return m.gt(a.first, b) || (m.is_pos(a.second) && m.eq(a.first, b));
|
||||
}
|
||||
|
||||
bool gt(mpq_inf const & a, mpq const & b, inf_kind k) const {
|
||||
bool gt(mpq_inf const & a, mpq const & b, inf_kind k) {
|
||||
if (m.gt(a.first, b))
|
||||
return true;
|
||||
if (m.eq(a.first, b)) {
|
||||
|
|
@ -166,17 +174,17 @@ public:
|
|||
return false;
|
||||
}
|
||||
|
||||
bool le(mpq_inf const & a, mpq_inf const & b) const { return !gt(a, b); }
|
||||
bool le(mpq_inf const & a, mpq_inf const & b) { return !gt(a, b); }
|
||||
|
||||
bool le(mpq_inf const & a, mpq const & b) const { return !gt(a, b); }
|
||||
bool le(mpq_inf const & a, mpq const & b) { return !gt(a, b); }
|
||||
|
||||
bool le(mpq_inf const & a, mpq const & b, inf_kind k) const { return !gt(a, b, k); }
|
||||
bool le(mpq_inf const & a, mpq const & b, inf_kind k) { return !gt(a, b, k); }
|
||||
|
||||
bool ge(mpq_inf const & a, mpq_inf const & b) const { return !lt(a, b); }
|
||||
bool ge(mpq_inf const & a, mpq_inf const & b) { return !lt(a, b); }
|
||||
|
||||
bool ge(mpq_inf const & a, mpq const & b) const { return !lt(a, b); }
|
||||
bool ge(mpq_inf const & a, mpq const & b) { return !lt(a, b); }
|
||||
|
||||
bool ge(mpq_inf const & a, mpq const & b, inf_kind k) const { return !lt(a, b, k); }
|
||||
bool ge(mpq_inf const & a, mpq const & b, inf_kind k) { return !lt(a, b, k); }
|
||||
|
||||
void add(mpq_inf const & a, mpq_inf const & b, mpq_inf & c) {
|
||||
m.add(a.first, b.first, c.first);
|
||||
|
|
@ -208,6 +216,16 @@ public:
|
|||
m.mul(b, a.second, c.second);
|
||||
}
|
||||
|
||||
void div(mpq_inf const & a, mpq const & b, mpq_inf & c) {
|
||||
m.div(a.first, b, c.first);
|
||||
m.div(a.second, b, c.second);
|
||||
}
|
||||
|
||||
void div(mpq_inf const & a, mpz const & b, mpq_inf & c) {
|
||||
m.div(a.first, b, c.first);
|
||||
m.div(a.second, b, c.second);
|
||||
}
|
||||
|
||||
void inc(mpq_inf & a) {
|
||||
m.inc(a.first);
|
||||
}
|
||||
|
|
@ -246,7 +264,7 @@ public:
|
|||
}
|
||||
}
|
||||
|
||||
std::string to_string(mpq_inf const & a) const;
|
||||
std::string to_string(mpq_inf const & a);
|
||||
};
|
||||
|
||||
typedef mpq_inf_manager<true> synch_mpq_inf_manager;
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue