From 26a3d2ca31072116a8f5c2cd7c0447f43acdc87e Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Tue, 21 Jan 2014 08:40:28 -0800 Subject: [PATCH 1/5] add stand-alone simplex Signed-off-by: Nikolaj Bjorner --- scripts/mk_project.py | 3 +- src/math/simplex/simplex.cpp | 26 ++ src/math/simplex/simplex.h | 157 ++++++++ src/math/simplex/simplex_def.h | 480 ++++++++++++++++++++++ src/math/simplex/sparse_matrix.h | 254 ++++++++++++ src/math/simplex/sparse_matrix_def.h | 580 +++++++++++++++++++++++++++ src/test/simplex.cpp | 24 ++ src/util/mpq_inf.cpp | 2 +- src/util/mpq_inf.h | 54 ++- 9 files changed, 1560 insertions(+), 20 deletions(-) create mode 100644 src/math/simplex/simplex.cpp create mode 100644 src/math/simplex/simplex.h create mode 100644 src/math/simplex/simplex_def.h create mode 100644 src/math/simplex/sparse_matrix.h create mode 100644 src/math/simplex/sparse_matrix_def.h create mode 100644 src/test/simplex.cpp diff --git a/scripts/mk_project.py b/scripts/mk_project.py index 6215e0fa0..952c22c20 100644 --- a/scripts/mk_project.py +++ b/scripts/mk_project.py @@ -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', diff --git a/src/math/simplex/simplex.cpp b/src/math/simplex/simplex.cpp new file mode 100644 index 000000000..494c0b6bb --- /dev/null +++ b/src/math/simplex/simplex.cpp @@ -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; + template class simplex; +}; diff --git a/src/math/simplex/simplex.h b/src/math/simplex/simplex.h new file mode 100644 index 000000000..0cb4852f2 --- /dev/null +++ b/src/math/simplex/simplex.h @@ -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 + 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 scoped_eps_numeral; + + typedef typename _scoped_numeral_vector scoped_eps_numeral_vector; + + typedef sparse_matrix matrix; + + struct var_lt { + bool operator()(var_t v1, var_t v2) const { return v1 < v2; } + }; + typedef heap 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 m_vars; + svector 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 + 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 diff --git a/src/math/simplex/simplex_def.h b/src/math/simplex/simplex_def.h new file mode 100644 index 000000000..b92ec7551 --- /dev/null +++ b/src/math/simplex/simplex_def.h @@ -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 simplex::row + simplex::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 + void simplex::del_row(row const& r) { + m_is_base[m_row2base[r.id()]] = false; + M.del(r); + } + + template + void simplex::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 + void simplex::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 + void simplex::unset_lower(var_t var) { + m_vars[var].m_lower_valid = false; + } + + template + void simplex::unset_upper(var_t var) { + m_vars[var].m_upper_valid = false; + } + + template + void simplex::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 simplex::eps_numeral const& + simplex::get_value(var_t v) { + return m_vars[v].m_value; + } + + template + void simplex::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 + void simplex::ensure_var(var_t v) { + while (v >= m_vars.size()) { + M.ensure_var(m_vars.size()); + m_vars.push_back(var_info()); + } + } + + template + lbool simplex::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 + void simplex::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 + void simplex::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 + void simplex::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 + bool simplex::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 + bool simplex::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 + bool simplex::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 + bool simplex::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 + bool simplex::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 simplex::var_t + simplex::select_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij) { + if (m_bland) { + return select_blands_pivot(x_i, is_below, out_a_ij); + } + if (is_below) { + return select_pivot_core(x_i, out_a_ij); + } + else { + return select_pivot_core(x_i, out_a_ij); + } + } + + /** + \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 + template + typename simplex::var_t + simplex::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 + int simplex::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 simplex::var_t + simplex::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 + lbool simplex::optimize(var_t v) { + // TBD SASSERT(is_feasible()); + // pick row for v and check if primal + // bounds are slack. + // return l_false for unbounded. + // return l_undef for canceled. + // return l_true for optimal. + return l_true; + } + + template + typename simplex::pivot_strategy_t + simplex::pivot_strategy() { + if (m_bland) { + return S_BLAND; + } + return S_DEFAULT; + } + + template + typename simplex::var_t + simplex::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 simplex::var_t + simplex::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 + bool simplex::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 + diff --git a/src/math/simplex/sparse_matrix.h b/src/math/simplex/sparse_matrix.h new file mode 100644 index 000000000..cecd1c44a --- /dev/null +++ b/src/math/simplex/sparse_matrix.h @@ -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 + 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 & cols); + void compress_if_needed(manager& m, vector & cols); + void save_var_pos(svector & result_map) const; + void reset_var_pos(svector & 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 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 m_dead_rows; // rows to recycle + vector m_columns; // per var + svector 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 diff --git a/src/math/simplex/sparse_matrix_def.h b/src/math/simplex/sparse_matrix_def.h new file mode 100644 index 000000000..4e46c823d --- /dev/null +++ b/src/math/simplex/sparse_matrix_def.h @@ -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 + sparse_matrix::_row::_row(): + m_size(0), + m_first_free_idx(-1) { + } + + template + void sparse_matrix::_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 sparse_matrix::_row_entry & + sparse_matrix::_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(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 + void sparse_matrix::_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 + void sparse_matrix::_row::compress(manager& m, vector & 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 + void sparse_matrix::_row::compress_if_needed(manager& m, vector & cols) { + if (size() *2 < num_entries()) { + compress(m, cols); + } + } + + /** + \brief Fill the map var -> pos/idx + */ + template + inline void sparse_matrix::_row::save_var_pos(svector & 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 + inline void sparse_matrix::_row::reset_var_pos(svector & 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 + int sparse_matrix::_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 + void sparse_matrix::column::reset() { + m_entries.reset(); + m_size = 0; + m_first_free_idx = -1; + } + + /** + \brief Remove holes (i.e., dead entries) from the column. + */ + template + void sparse_matrix::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 + inline void sparse_matrix::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 + void sparse_matrix::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 + const typename sparse_matrix::col_entry * + sparse_matrix::column::get_first_col_entry() const { + typename svector::const_iterator it = m_entries.begin(); + typename svector::const_iterator end = m_entries.end(); + for (; it != end; ++it) { + if (!it->is_dead()) { + return it; + } + } + return 0; + } + + template + typename sparse_matrix::col_entry & + sparse_matrix::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 + void sparse_matrix::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 + sparse_matrix::~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 + void sparse_matrix::ensure_var(var_t v) { + while (m_columns.size() <= v) { + m_columns.push_back(column()); + m_var_pos.push_back(-1); + } + } + + template + typename sparse_matrix::row + sparse_matrix::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 + void sparse_matrix::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 + void sparse_matrix::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 + void sparse_matrix::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 + void sparse_matrix::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 + void sparse_matrix::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 + void sparse_matrix::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 + void sparse_matrix::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 + bool sparse_matrix::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 + bool sparse_matrix::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 + bool sparse_matrix::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 + void sparse_matrix::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 diff --git a/src/test/simplex.cpp b/src/test/simplex.cpp new file mode 100644 index 000000000..ca69cec1f --- /dev/null +++ b/src/test/simplex.cpp @@ -0,0 +1,24 @@ +#include "sparse_matrix.h" +#include "sparse_matrix_def.h" +#include "simplex.h" +#include "simplex_def.h" + +typedef simplex::simplex Simplex; + +void tst_simplex() { + simplex::sparse_matrix M; + Simplex S; + + S.make_feasible(); + + unsynch_mpz_manager m; + scoped_mpz_vector coeffs(m); + svector 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()); +} diff --git a/src/util/mpq_inf.cpp b/src/util/mpq_inf.cpp index 6f6d805f2..cb3232d6c 100644 --- a/src/util/mpq_inf.cpp +++ b/src/util/mpq_inf.cpp @@ -19,7 +19,7 @@ Revision History: #include"mpq_inf.h" template -std::string mpq_inf_manager::to_string(mpq_inf const & a) const { +std::string mpq_inf_manager::to_string(mpq_inf const & a) { if (m.is_zero(a.second)) return m.to_string(a.first); diff --git a/src/util/mpq_inf.h b/src/util/mpq_inf.h index 217b58567..eef066a6b 100644 --- a/src/util/mpq_inf.h +++ b/src/util/mpq_inf.h @@ -26,12 +26,12 @@ typedef std::pair mpq_inf; template class mpq_inf_manager { - mpq_manager & m; + mpq_manager m; double m_inf; public: typedef mpq_inf numeral; - mpq_inf_manager(mpq_manager & _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 synch_mpq_inf_manager; From f6fd426c2808442ad9ca5ff3c7c4c3bdabcd5a26 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Tue, 21 Jan 2014 08:46:02 -0800 Subject: [PATCH 2/5] moved network flow Signed-off-by: Nikolaj Bjorner --- src/{smt => math/simplex}/network_flow.h | 0 src/{smt => math/simplex}/network_flow_def.h | 0 src/math/simplex/simplex_def.h | 1 - src/smt/theory_diff_logic.h | 1 - src/smt/theory_diff_logic_def.h | 9 ++++++++- 5 files changed, 8 insertions(+), 3 deletions(-) rename src/{smt => math/simplex}/network_flow.h (100%) rename src/{smt => math/simplex}/network_flow_def.h (100%) diff --git a/src/smt/network_flow.h b/src/math/simplex/network_flow.h similarity index 100% rename from src/smt/network_flow.h rename to src/math/simplex/network_flow.h diff --git a/src/smt/network_flow_def.h b/src/math/simplex/network_flow_def.h similarity index 100% rename from src/smt/network_flow_def.h rename to src/math/simplex/network_flow_def.h diff --git a/src/math/simplex/simplex_def.h b/src/math/simplex/simplex_def.h index b92ec7551..0e05de291 100644 --- a/src/math/simplex/simplex_def.h +++ b/src/math/simplex/simplex_def.h @@ -213,7 +213,6 @@ namespace simplex { 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); } } diff --git a/src/smt/theory_diff_logic.h b/src/smt/theory_diff_logic.h index ed52b76fc..dff1b718c 100644 --- a/src/smt/theory_diff_logic.h +++ b/src/smt/theory_diff_logic.h @@ -38,7 +38,6 @@ Revision History: #include"numeral_factory.h" #include"smt_clause.h" #include"theory_opt.h" -#include"network_flow.h" // The DL theory can represent term such as n + k, where n is an enode and k is a numeral. namespace smt { diff --git a/src/smt/theory_diff_logic_def.h b/src/smt/theory_diff_logic_def.h index 1982f2eb7..3dc9d3a7f 100644 --- a/src/smt/theory_diff_logic_def.h +++ b/src/smt/theory_diff_logic_def.h @@ -28,7 +28,6 @@ Revision History: #include"ast_pp.h" #include"warning.h" #include"smt_model_generator.h" -#include"network_flow_def.h" #include"model_implicant.h" using namespace smt; @@ -998,6 +997,10 @@ void theory_diff_logic::get_implied_bound_antecedents(edge_id bridge_edge, template inf_eps_rational theory_diff_logic::maximize(theory_var v) { + +#ifdef 0 + // disabled until fixed. + objective_term const& objective = m_objectives[v]; IF_VERBOSE(1, @@ -1058,6 +1061,10 @@ inf_eps_rational theory_diff_logic::maximize(theory_var v) { IF_VERBOSE(1, verbose_stream() << "Unbounded objective" << std::endl;); return inf_eps_rational::infinity(); } + +#endif + return inf_eps_rational::infinity(); + } template From f68eff32763a2811154b05aa38db07018d19fba6 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Tue, 21 Jan 2014 09:06:30 -0800 Subject: [PATCH 3/5] move network flow code Signed-off-by: Nikolaj Bjorner --- src/smt/theory_diff_logic_def.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smt/theory_diff_logic_def.h b/src/smt/theory_diff_logic_def.h index 3dc9d3a7f..8e3818a51 100644 --- a/src/smt/theory_diff_logic_def.h +++ b/src/smt/theory_diff_logic_def.h @@ -998,7 +998,7 @@ void theory_diff_logic::get_implied_bound_antecedents(edge_id bridge_edge, template inf_eps_rational theory_diff_logic::maximize(theory_var v) { -#ifdef 0 +#if 0 // disabled until fixed. objective_term const& objective = m_objectives[v]; From c14c65465a2602896ba3462ced9dadbfdbddaf27 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Sun, 26 Jan 2014 19:46:42 -0800 Subject: [PATCH 4/5] working on stand-alone simplex Signed-off-by: Nikolaj Bjorner --- src/math/simplex/simplex.h | 28 ++- src/math/simplex/simplex_def.h | 308 ++++++++++++++++++++--- src/math/simplex/sparse_matrix.h | 2 +- src/math/simplex/sparse_matrix_def.h | 6 +- src/opt/weighted_maxsat.cpp | 1 + src/smt/theory_arith_aux.h | 3 +- src/tactic/core/pb_preprocess_tactic.cpp | 22 +- src/util/mpq_inf.h | 6 + src/util/scoped_numeral.h | 6 + 9 files changed, 332 insertions(+), 50 deletions(-) diff --git a/src/math/simplex/simplex.h b/src/math/simplex/simplex.h index 0cb4852f2..974aa7cfb 100644 --- a/src/math/simplex/simplex.h +++ b/src/math/simplex/simplex.h @@ -36,6 +36,7 @@ Notes: #include "mpq_inf.h" #include "heap.h" #include "lbool.h" +#include "uint_set.h" namespace simplex { @@ -84,7 +85,7 @@ namespace simplex { }; static const var_t null_var = UINT_MAX; - matrix M; + mutable matrix M; unsigned m_max_iterations; volatile bool m_cancel; var_heap m_to_patch; @@ -93,14 +94,17 @@ namespace simplex { vector m_vars; svector m_row2base; bool m_bland; + unsigned m_blands_rule_threshold; random_gen m_random; + uint_set m_left_basis; public: simplex(): m_max_iterations(UINT_MAX), m_cancel(false), m_to_patch(1024), - m_bland(false) {} + m_bland(false), + m_blands_rule_threshold(1000) {} typedef typename matrix::row row; typedef typename matrix::row_iterator row_iterator; @@ -117,7 +121,7 @@ namespace simplex { void set_cancel(bool f) { m_cancel = f; } void set_max_iterations(unsigned m) { m_max_iterations = m; } lbool make_feasible(); - lbool optimize(var_t var); + lbool minimize(var_t var); eps_numeral const& get_value(var_t v); void display(std::ostream& out) const; @@ -128,17 +132,26 @@ namespace simplex { var_t select_smallest_var() { return m_to_patch.empty()?null_var:m_to_patch.erase_min(); } var_t select_error_var(bool least); // row get_infeasible_row() { } - void check_blands_rule(var_t v) { } + void check_blands_rule(var_t v, unsigned& num_repeated); bool make_var_feasible(var_t x_i); void update_and_pivot(var_t x_i, var_t x_j, numeral const& a_ij, eps_numeral const& new_value); void update_value(var_t v, eps_numeral const& delta); void update_value_core(var_t v, eps_numeral const& delta); + void pivot(var_t x_i, var_t x_j, numeral const& a_ij); + void move_to_bound(var_t x, bool to_lower); var_t select_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij); - var_t select_blands_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij); - template - var_t select_pivot_core(var_t x_i, scoped_numeral& out_a_ij); + var_t select_pivot_blands(var_t x_i, bool is_below, scoped_numeral& out_a_ij); + var_t select_pivot_core(var_t x_i, bool is_below, scoped_numeral& out_a_ij); int get_num_non_free_dep_vars(var_t x_j, int best_so_far); + var_t pick_var_to_leave(var_t x_j, bool inc, scoped_eps_numeral& gain, scoped_numeral& new_a_ij); + + + void select_pivot_primal(var_t v, var_t& x_i, var_t& x_j, scoped_numeral& a_ij, bool& inc); + + + bool at_lower(var_t v) const; + bool at_upper(var_t v) const; bool below_lower(var_t v) const; bool above_upper(var_t v) const; bool above_lower(var_t v) const; @@ -150,6 +163,7 @@ namespace simplex { bool is_base(var_t x) const { return m_vars[x].m_is_base; } bool well_formed() const; + bool is_feasible() const; }; }; diff --git a/src/math/simplex/simplex_def.h b/src/math/simplex/simplex_def.h index 0e05de291..26c7e2e03 100644 --- a/src/math/simplex/simplex_def.h +++ b/src/math/simplex/simplex_def.h @@ -44,7 +44,8 @@ namespace simplex { template void simplex::del_row(row const& r) { - m_is_base[m_row2base[r.id()]] = false; + m_vars[m_row2base[r.id()]].m_is_base = false; + m_row2base[r.id()] = null_var; M.del(r); } @@ -102,11 +103,11 @@ namespace simplex { var_info const& vi = m_vars[i]; out << "v" << i << " "; if (vi.m_is_base) out << "b:" << vi.m_base2row << " "; - em.display(out, vi.m_value); + out << em.to_string(vi.m_value); out << " ["; - if (vi.m_lower_valid) em.display(out, vi.m_lower) else out << "-oo"; + if (vi.m_lower_valid) out << em.to_string(vi.m_lower); else out << "-oo"; out << ":"; - if (vi.m_upper_valid) em.display(out, vi.m_upper) else out << "oo"; + if (vi.m_upper_valid) out << em.to_string(vi.m_upper); else out << "oo"; out << "]"; out << "\n"; } @@ -122,13 +123,15 @@ namespace simplex { template lbool simplex::make_feasible() { + m_left_basis.reset(); unsigned num_iterations = 0; + unsigned num_repeated = 0; var_t v = null_var; while ((v = select_var_to_fix()) != null_var) { if (m_cancel || num_iterations > m_max_iterations) { return l_undef; } - check_blands_rule(v); + check_blands_rule(v, num_repeated); if (!make_var_feasible(v)) { return l_false; } @@ -157,7 +160,6 @@ namespace simplex { SASSERT(is_base(x_i)); SASSERT(!is_base(x_j)); var_info& x_iI = m_vars[x_i]; - var_info& x_jI = m_vars[x_j]; scoped_eps_numeral theta(em); theta = x_iI.m_value; theta -= new_value; @@ -166,6 +168,13 @@ namespace simplex { em.div(theta, a_ij, theta); update_value(x_j, theta); SASSERT(em.eq(x_iI.m_value, new_value)); + pivot(x_i, x_j, a_ij); + } + + template + void simplex::pivot(var_t x_i, var_t x_j, numeral const& a_ij) { + var_info& x_iI = m_vars[x_i]; + var_info& x_jI = m_vars[x_j]; unsigned r_i = x_iI.m_base2row; x_jI.m_base2row = r_i; m_row2base[r_i] = x_j; @@ -196,6 +205,9 @@ namespace simplex { template void simplex::update_value(var_t v, eps_numeral const& delta) { + if (em.is_zero(delta)) { + return; + } update_value_core(v, delta); col_iterator it = M.col_begin(v), end = M.col_end(v); @@ -249,6 +261,18 @@ namespace simplex { return !vi.m_upper_valid || em.lt(vi.m_value, vi.m_upper); } + template + bool simplex::at_lower(var_t v) const { + var_info const& vi = m_vars[v]; + return vi.m_lower_valid && em.eq(vi.m_value, vi.m_lower); + } + + template + bool simplex::at_upper(var_t v) const { + var_info const& vi = m_vars[v]; + return vi.m_upper_valid && em.eq(vi.m_value, vi.m_upper); + } + template bool simplex::make_var_feasible(var_t x_i) { scoped_numeral a_ij(m); @@ -274,20 +298,15 @@ namespace simplex { } /** - \brief Wrapper for select_blands_pivot_core and select_pivot_core + \brief Wrapper for select_pivot_blands and select_pivot_core */ template typename simplex::var_t simplex::select_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij) { if (m_bland) { - return select_blands_pivot(x_i, is_below, out_a_ij); - } - if (is_below) { - return select_pivot_core(x_i, out_a_ij); - } - else { - return select_pivot_core(x_i, out_a_ij); + return select_pivot_blands(x_i, is_below, out_a_ij); } + return select_pivot_core(x_i, is_below, out_a_ij); } /** @@ -300,9 +319,8 @@ namespace simplex { bound (above its upper bound). */ template - template typename simplex::var_t - simplex::select_pivot_core(var_t x_i, scoped_numeral & out_a_ij) { + simplex::select_pivot_core(var_t x_i, bool is_below, scoped_numeral & out_a_ij) { SASSERT(is_base(x_i)); var_t max = get_num_vars(); var_t result = max; @@ -315,11 +333,13 @@ namespace simplex { for (; it != end; ++it) { var_t x_j = it->m_var; + if (x_i == x_j) continue; numeral const & a_ij = it->m_coeff; bool is_neg = is_below ? m.is_neg(a_ij) : m.is_pos(a_ij); bool is_pos = !is_neg; - if (x_i != x_j && ((is_pos && above_lower(x_j)) || (is_neg && below_upper(x_j)))) { + bool can_pivot = ((is_pos && above_lower(x_j)) || (is_neg && below_upper(x_j))); + if (can_pivot) { int num = get_num_non_free_dep_vars(x_j, best_so_far); unsigned col_sz = M.column_size(x_j); if (num < best_so_far || (num == best_so_far && col_sz < best_col_sz)) { @@ -360,16 +380,15 @@ namespace simplex { return result; } - /** \brief Using Bland's rule, select a variable x_j in the row r defining the base var x_i, - s.t. x_j can be used to patch the error in x_i. Return null_theory_var + s.t. x_j can be used to patch the error in x_i. Return null_var if there is no x_j. Otherwise, return x_j and store its coefficient in out_a_ij. */ template typename simplex::var_t - simplex::select_blands_pivot(var_t x_i, bool is_below, scoped_numeral & out_a_ij) { + simplex::select_pivot_blands(var_t x_i, bool is_below, scoped_numeral & out_a_ij) { SASSERT(is_base(x_i)); unsigned max = get_num_vars(); var_t result = max; @@ -379,8 +398,7 @@ namespace simplex { var_t x_j = it->m_var; numeral const & a_ij = it->m_coeff; bool is_neg = is_below ? m.is_neg(a_ij) : m.is_pos(a_ij); - bool is_pos = !is_neg; - if (x_i != x_j && ((is_pos && above_lower(x_j)) || (is_neg && below_upper(x_j)))) { + if (x_i != x_j && ((!is_neg && above_lower(x_j)) || (is_neg && below_upper(x_j)))) { SASSERT(!is_base(x_j)); if (x_j < result) { result = x_j; @@ -392,16 +410,230 @@ namespace simplex { } template - lbool simplex::optimize(var_t v) { - // TBD SASSERT(is_feasible()); - // pick row for v and check if primal - // bounds are slack. - // return l_false for unbounded. - // return l_undef for canceled. - // return l_true for optimal. + lbool simplex::minimize(var_t v) { + + // minimize v, such that tableau is feasible. + // Assume there are no bounds on v. + // Let k*v + c*x = 0 e.g, maximize c*x over + // tableau constraints: + // + // max { c*x | A*x = 0 and l <= x <= u } + // + // start with feasible assigment + // A*x0 = 0 and l <= x0 <= u + // + // Identify pivot: i, j: such that x_i is base, + // there is a row k1*x_i + k2*x_j + R = 0 + // and a delta such that: + // + // x_i' <- x_i + delta + // x_j' <- x_j - delta*k1/k2 + // l_i <= x_i' <= u_i + // l_j <= x_j' <= u_j + // and c*x' > c*x + // e.g., c*x := c_i*x_i + c_j*x_j + ... + // and c_i*delta > c_j*delta*k1/k2 + // and x_i < u_i (if delta > 0), l_i < x_i (if delta < 0) + // and l_j < x_j (if delta > 0), x_j < u_j (if delta < 0) + // + // update all rows, including c*x, using the pivot. + // + // If there is c_i*x_i in c*x such that c_i > 0 + // and upper_i = oo and complementary lower_j = -oo + // then the objective is unbounded. + // + // There is a singularity if there is a pivot such that + // c_i*delta == c_j*delta*k1/k2, e.g., nothing is improved, + // pivot, but use bland's rule to ensure + // convergence in the limit. + // + + SASSERT(is_feasible()); + scoped_eps_numeral delta(em); + scoped_numeral a_ij(m); + var_t x_i, x_j; + bool inc; + + while (true) { + if (m_cancel) { + return l_undef; + } + select_pivot_primal(v, x_i, x_j, a_ij, inc); + if (x_j == null_var) { + // optimal + return l_true; + } + var_info& vj = m_vars[x_j]; + if (x_i == null_var) { + if (inc && vj.m_upper_valid) { + delta = vj.m_upper; + delta -= vj.m_value; + update_value(x_j, delta); + } + else if (!inc && vj.m_lower_valid) { + delta = vj.m_lower; + delta -= vj.m_value; + update_value(x_j, delta); + } + else { + // unbounded + return l_false; + } + continue; + } + + // TBD: Change the value of x_j directly without pivoting: + // + // if (!vj.is_fixed() && vj.bounded() && gain >= upper - lower) { + // + // } + // + + pivot(x_i, x_j, a_ij); + move_to_bound(x_i, inc == m.is_pos(a_ij)); + } return l_true; } + template + void simplex::move_to_bound(var_t x, bool to_lower) { + scoped_eps_numeral delta(em), delta2(em); + var_info& vi = m_vars[x]; + if (to_lower) { + em.sub(vi.m_value, vi.m_lower, delta); + } + else { + em.sub(vi.m_upper, vi.m_value, delta); + } + col_iterator it = M.col_begin(x), end = M.col_end(x); + for (; it != end && is_pos(delta); ++it) { + var_t s = m_row2base[it.get_row().id()]; + var_info& vs = m_vars[s]; + numeral const& coeff = it.get_row_entry().m_coeff; + SASSERT(!m.is_zero(coeff)); + bool inc_s = (m.is_pos(coeff) == to_lower); + eps_numeral const* bound = 0; + if (inc_s && vs.m_upper_valid) { + bound = &vs.m_upper; + } + else if (!inc_s && vs.m_lower_valid) { + bound = &vs.m_lower; + } + if (bound) { + em.sub(*bound, vs.m_value, delta2); + em.div(delta2, coeff, delta2); + abs(delta2); + if (delta2 < delta) { + delta = delta2; + } + } + } + if (to_lower) { + delta.neg(); + } + update_value(x, delta); + } + + template + void simplex::select_pivot_primal(var_t v, var_t& x_i, var_t& x_j, scoped_numeral& a_ij, bool& inc) { + row r(m_vars[v].m_base2row); + row_iterator it = M.row_begin(r), end = M.row_end(r); + + scoped_eps_numeral gain(em), new_gain(em); + scoped_numeral new_a_ij(m); + x_i = null_var; + x_j = null_var; + inc = false; + + for (; it != end; ++it) { + var_t x = it->m_var; + if (x == v) continue; + bool is_pos = m.is_pos(it->m_coeff); + if ((is_pos && at_upper(x)) || (!is_pos && at_lower(x))) { + continue; // variable cannot be used for improving bounds. + // TBD check? + } + var_t y = pick_var_to_leave(x, is_pos, new_gain, new_a_ij); + if (y == null_var) { + // unbounded. + x_i = y; + x_j = x; + inc = is_pos; + a_ij = new_a_ij; + break; + } + bool better = + (new_gain > gain) || + ((is_zero(new_gain) && is_zero(gain) && (x_i == null_var || y < x_i))); + + if (better) { + x_i = y; + x_j = x; + inc = is_pos; + gain = new_gain; + a_ij = new_a_ij; + } + } + } + + template + typename simplex::var_t + simplex::pick_var_to_leave( + var_t x_j, bool inc, + scoped_eps_numeral& gain, scoped_numeral& new_a_ij) { + var_t x_i = null_var; + gain.reset(); + scoped_eps_numeral curr_gain(em); + col_iterator it = M.col_begin(x_j), end = M.col_end(x_j); + for (; it != end; ++it) { + row r = it.get_row(); + var_t s = m_row2base[r.id()]; + var_info& vi = m_vars[s]; + numeral const& a_ij = it.get_row_entry().m_coeff; + numeral const& a_ii = vi.m_base_coeff; + bool inc_s = m.is_neg(a_ij) ? inc : !inc; + if ((inc_s && !vi.m_upper_valid) || (!inc_s && !vi.m_lower_valid)) { + continue; + } + // + // current gain: (value(x_i)-bound)*a_ii/a_ij + // + curr_gain = vi.m_value; + curr_gain -= inc_s?vi.m_upper:vi.m_lower; + em.mul(curr_gain, a_ii, curr_gain); + em.div(curr_gain, a_ij, curr_gain); + if (is_neg(curr_gain)) { + curr_gain.neg(); + } + if (x_i == null_var || (curr_gain < gain) || + (is_zero(gain) && is_zero(curr_gain) && s < x_i)) { + x_i = s; + gain = curr_gain; + new_a_ij = a_ij; + } + } + TRACE("simplex", tout << "x_i v" << x_i << "\n";); + return x_i; + } + + template + void simplex::check_blands_rule(var_t v, unsigned& num_repeated) { + if (m_bland) + return; + if (m_left_basis.contains(v)) { + num_repeated++; + if (num_repeated > m_blands_rule_threshold) { + TRACE("simplex", tout << "using blands rule, " << num_repeated << "\n";); + // std::cerr << "BLANDS RULE...\n"; + m_bland = true; + } + } + else { + m_left_basis.insert(v); + } + } + + template typename simplex::pivot_strategy_t simplex::pivot_strategy() { @@ -463,15 +695,31 @@ namespace simplex { SASSERT(M.well_formed()); for (unsigned i = 0; i < m_row2base.size(); ++i) { var_t s = m_row2base[i]; - SASSERT(i == m_vars[s].m_base2row); // unless row is deleted. + if (s == null_var) continue; + SASSERT(i == m_vars[s].m_base2row); // // TBD: extract coefficient of base variable and compare // with m_vars[s].m_base_coeff; // + // check that sum of assignments add up to 0 for every row. + row_iterator it = M.row_begin(row(i)), end = M.row_end(row(i)); + scoped_eps_numeral sum(em), tmp(em); + for (; it != end; ++it) { + em.mul(m_vars[it->m_var].m_value, it->m_coeff, tmp); + sum += tmp; + } + SASSERT(em.is_zero(sum)); } return true; } + template + bool simplex::is_feasible() const { + for (unsigned i = 0; i < m_vars.size(); ++i) { + if (below_lower(i) || above_upper(i)) return false; + } + return true; + } }; diff --git a/src/math/simplex/sparse_matrix.h b/src/math/simplex/sparse_matrix.h index cecd1c44a..6a85d69f0 100644 --- a/src/math/simplex/sparse_matrix.h +++ b/src/math/simplex/sparse_matrix.h @@ -226,7 +226,7 @@ namespace simplex { col_iterator col_begin(int v) const { return col_iterator(m_columns[v], m_rows, true); } col_iterator col_end(int v) const { return col_iterator(m_columns[v], m_rows, false); } - void display(std::ostream& out) const; + void display(std::ostream& out); bool well_formed() const; diff --git a/src/math/simplex/sparse_matrix_def.h b/src/math/simplex/sparse_matrix_def.h index 4e46c823d..20e9d676f 100644 --- a/src/math/simplex/sparse_matrix_def.h +++ b/src/math/simplex/sparse_matrix_def.h @@ -448,7 +448,7 @@ namespace simplex { */ template void sparse_matrix::del(row r) { - _row& rw = r.id(); + _row& rw = m_rows[r.id()]; for (unsigned i = 0; i < rw.m_entries.size(); ++i) { _row_entry& e = rw.m_entries[i]; if (!e.is_dead()) { @@ -561,12 +561,12 @@ namespace simplex { \brief display method */ template - void sparse_matrix::display(std::ostream& out) const { + void sparse_matrix::display(std::ostream& out) { for (unsigned i = 0; i < m_rows.size(); ++i) { if (m_rows[i].size() == 0) continue; row_iterator it = row_begin(row(i)), end = row_end(row(i)); for (; it != end; ++it) { - m.display(out, it->m_coeff) + m.display(out, it->m_coeff); out << "*v" << it->m_var << " "; } out << "\n"; diff --git a/src/opt/weighted_maxsat.cpp b/src/opt/weighted_maxsat.cpp index b5c53c3be..36047d434 100644 --- a/src/opt/weighted_maxsat.cpp +++ b/src/opt/weighted_maxsat.cpp @@ -867,6 +867,7 @@ namespace opt { m_upper += w; if (wmax < w) wmax = w; b = m.mk_fresh_const("b", m.mk_bool_sort()); + block.push_back(b); expr* bb = b; s.mc().insert(b->get_decl()); a = m.mk_fresh_const("a", m.mk_bool_sort()); diff --git a/src/smt/theory_arith_aux.h b/src/smt/theory_arith_aux.h index 77845fee6..e430b6515 100644 --- a/src/smt/theory_arith_aux.h +++ b/src/smt/theory_arith_aux.h @@ -1359,8 +1359,7 @@ namespace smt { SASSERT(valid_row_assignment()); SASSERT(satisfy_bounds()); result = skipped_row?BEST_EFFORT:OPTIMIZED; - break; - } + break; } if (x_i == null_theory_var) { // can increase/decrease x_j as much as we want. diff --git a/src/tactic/core/pb_preprocess_tactic.cpp b/src/tactic/core/pb_preprocess_tactic.cpp index 166cecc97..51e92ef11 100644 --- a/src/tactic/core/pb_preprocess_tactic.cpp +++ b/src/tactic/core/pb_preprocess_tactic.cpp @@ -244,6 +244,7 @@ public: expr* arg = args1[0].get(); bool neg = m.is_not(arg, arg); if (!is_uninterp_const(arg)) continue; + if (!m_vars.contains(to_app(arg))) continue; rec const& r = m_vars.find(to_app(arg)); unsigned_vector const& pos = neg?r.neg:r.pos; for (unsigned j = 0; j < pos.size(); ++j) { @@ -251,7 +252,8 @@ public: if (k == i) continue; if (!to_ge(g->form(k), args2, coeffs2, k2)) continue; if (subsumes(args1, coeffs1, k1, args2, coeffs2, k2)) { - g->update(k, m.mk_true()); + IF_VERBOSE(3, verbose_stream() << "replace " << mk_pp(g->form(k), m) << "\n";); + g->update(k, m.mk_true()); m_progress = true; } } @@ -518,6 +520,7 @@ private: expr_ref tmp1(m), tmp2(m); expr* fml1 = g->form(idx1); expr* fml2 = g->form(idx2); + TRACE("pb", tout << mk_pp(fml1, m) << " " << mk_pp(fml2, m) << "\n";); expr_ref_vector args1(m), args2(m); vector coeffs1, coeffs2; rational k1, k2; @@ -533,6 +536,7 @@ private: m.is_not(x, x); if (!is_app(x)) return; if (!m_vars.contains(to_app(x))) return; + TRACE("pb", tout << mk_pp(x, m) << "\n";); rec const& r = m_vars.find(to_app(x)); if (r.pos.size() != 1 || r.neg.size() != 1) return; if (r.pos[0] != idx2 && r.neg[0] != idx2) return; @@ -632,13 +636,17 @@ private: m_r.set_substitution(&sub); for (unsigned i = 0; i < positions.size(); ++i) { unsigned idx = positions[i]; - expr* f = g->form(idx); + expr_ref f(m); + f = g->form(idx); if (!m.is_true(f)) { m_r(f, tmp); - TRACE("pb", tout << mk_pp(f, m) << " -> " << tmp - << " by " << mk_pp(e, m) << " |-> " << mk_pp(v, m) << "\n";); - g->update(idx, tmp); // proof & dependencies. - m_progress = true; + if (tmp != f) { + TRACE("pb", tout << mk_pp(f, m) << " -> " << tmp + << " by " << mk_pp(e, m) << " |-> " << mk_pp(v, m) << "\n";); + IF_VERBOSE(3, verbose_stream() << "replace " << mk_pp(f, m) << " -> " << tmp << "\n";); + g->update(idx, tmp); // proof & dependencies. + m_progress = true; + } } } m_r.set_substitution(0); @@ -653,7 +661,7 @@ private: bool found = false; for (unsigned j = 0; !found && j < args2.size(); ++j) { if (args1[i] == args2[j]) { - if (coeffs1[1] > coeffs2[j]) return false; + if (coeffs1[i] > coeffs2[j]) return false; found = true; } } diff --git a/src/util/mpq_inf.h b/src/util/mpq_inf.h index eef066a6b..4b6a7d33d 100644 --- a/src/util/mpq_inf.h +++ b/src/util/mpq_inf.h @@ -239,6 +239,12 @@ public: m.neg(a.second); } + void abs(mpq_inf & a) { + if (is_neg(a)) { + neg(a); + } + } + void ceil(mpq_inf const & a, mpq & b) { if (m.is_int(a.first)) { // special cases for k - delta*epsilon where k is an integer diff --git a/src/util/scoped_numeral.h b/src/util/scoped_numeral.h index e03b56336..54d55f827 100644 --- a/src/util/scoped_numeral.h +++ b/src/util/scoped_numeral.h @@ -133,6 +133,12 @@ public: return a.m().is_nonpos(a); } + friend _scoped_numeral abs(_scoped_numeral const& a) { + _scoped_numeral res(a); + a.m().abs(res); + return res; + } + void neg() { m().neg(m_num); } From 363af825c06a5af8f03d842b48131d54c2080582 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Sun, 26 Jan 2014 20:25:36 -0800 Subject: [PATCH 5/5] working on stand-alone simplex Signed-off-by: Nikolaj Bjorner --- src/math/simplex/simplex_def.h | 16 ++++++++++++---- src/smt/theory_arith_aux.h | 3 ++- src/test/main.cpp | 1 + src/test/simplex.cpp | 18 +++++++++++++++++- 4 files changed, 32 insertions(+), 6 deletions(-) diff --git a/src/math/simplex/simplex_def.h b/src/math/simplex/simplex_def.h index 26c7e2e03..46915695f 100644 --- a/src/math/simplex/simplex_def.h +++ b/src/math/simplex/simplex_def.h @@ -29,9 +29,13 @@ namespace simplex { for (unsigned i = 0; !found && i < num_vars; ++i) found = vars[i] == base; SASSERT(found); ); + scoped_numeral base_coeff(m); row r = M.mk_row(); for (unsigned i = 0; i < num_vars; ++i) { M.add(r, coeffs[i], vars[i]); + if (vars[i] == base) { + m.set(base_coeff, coeffs[i]); + } } while (m_row2base.size() <= r.id()) { m_row2base.push_back(null_var); @@ -39,6 +43,7 @@ namespace simplex { m_row2base[r.id()] = base; m_vars[base].m_base2row = r.id(); m_vars[base].m_is_base = true; + m.set(m_vars[base].m_base_coeff, base_coeff); return r; } @@ -54,7 +59,8 @@ namespace simplex { var_info& vi = m_vars[var]; em.set(vi.m_lower, b); vi.m_lower_valid = true; - if (em.lt(vi.m_value, b)) { + SASSERT(!vi.m_upper_valid || em.le(b, vi.m_upper)); + if (!vi.m_is_base && em.lt(vi.m_value, b)) { scoped_eps_numeral delta(em); em.sub(b, vi.m_value, delta); update_value(var, delta); @@ -66,7 +72,8 @@ namespace simplex { var_info& vi = m_vars[var]; em.set(vi.m_upper, b); vi.m_upper_valid = true; - if (em.gt(vi.m_value, b)) { + SASSERT(!vi.m_lower_valid || em.le(vi.m_lower, b)); + if (!vi.m_is_base && em.gt(vi.m_value, b)) { scoped_eps_numeral delta(em); em.sub(b, vi.m_value, delta); update_value(var, delta); @@ -102,13 +109,13 @@ namespace simplex { 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 << " "; out << em.to_string(vi.m_value); out << " ["; if (vi.m_lower_valid) out << em.to_string(vi.m_lower); else out << "-oo"; out << ":"; if (vi.m_upper_valid) out << em.to_string(vi.m_upper); else out << "oo"; out << "]"; + if (vi.m_is_base) out << " b:" << vi.m_base2row; out << "\n"; } } @@ -176,8 +183,9 @@ namespace simplex { 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; + x_jI.m_base2row = r_i; + m.set(x_jI.m_base_coeff, a_ij); x_jI.m_is_base = true; x_iI.m_is_base = false; if (outside_bounds(x_j)) { diff --git a/src/smt/theory_arith_aux.h b/src/smt/theory_arith_aux.h index e430b6515..eeafb1aa5 100644 --- a/src/smt/theory_arith_aux.h +++ b/src/smt/theory_arith_aux.h @@ -1359,7 +1359,8 @@ namespace smt { SASSERT(valid_row_assignment()); SASSERT(satisfy_bounds()); result = skipped_row?BEST_EFFORT:OPTIMIZED; - break; } + break; + } if (x_i == null_theory_var) { // can increase/decrease x_j as much as we want. diff --git a/src/test/main.cpp b/src/test/main.cpp index f6a4eab29..4912916a1 100644 --- a/src/test/main.cpp +++ b/src/test/main.cpp @@ -218,6 +218,7 @@ int main(int argc, char ** argv) { TST(expr_substitution); TST(sorting_network); TST(theory_pb); + TST(simplex); } void initialize_mam() {} diff --git a/src/test/simplex.cpp b/src/test/simplex.cpp index ca69cec1f..def6f8be3 100644 --- a/src/test/simplex.cpp +++ b/src/test/simplex.cpp @@ -2,6 +2,7 @@ #include "sparse_matrix_def.h" #include "simplex.h" #include "simplex_def.h" +#include "mpq_inf.h" typedef simplex::simplex Simplex; @@ -9,9 +10,13 @@ void tst_simplex() { simplex::sparse_matrix M; Simplex S; - S.make_feasible(); + std::cout << "simplex\n"; + + lbool is_sat = S.make_feasible(); + std::cout << "feasible: " << is_sat << "\n"; unsynch_mpz_manager m; + unsynch_mpq_inf_manager em; scoped_mpz_vector coeffs(m); svector vars; for (unsigned i = 0; i < 5; ++i) { @@ -21,4 +26,15 @@ void tst_simplex() { } Simplex::row r = S.add_row(1, coeffs.size(), vars.c_ptr(), coeffs.c_ptr()); + is_sat = S.make_feasible(); + std::cout << "feasible: " << is_sat << "\n"; + S.display(std::cout); + _scoped_numeral num(em); + num = std::make_pair(mpq(1), mpq(0)); + S.set_lower(0, num); + S.set_upper(0, num); + + is_sat = S.make_feasible(); + std::cout << "feasible: " << is_sat << "\n"; + S.display(std::cout); }