mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-14 21:08:46 +00:00
Code Activity

add var_register

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

fill the matrix A in hnf_cutter

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

fill the matrix A in hnf_cutter

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

first steps of hnf cutter

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

handle generated cases in hnf

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

call hnf only for a full rank matrix

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

get (H reversed) * b

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

finding the cut row randomly, exiting if is not there

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

produce first cuts with hnf

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

produce first cuts with hnf

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

define by lp_settings if to avoid calling hnf_cutter when the solution is not on the boundary

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

hnf

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

revert to the previous version

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2018-05-17 15:59:38 -07:00
parent 3b5337823a
commit 82eb80de6d
17 changed files with 482 additions and 129 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/smt/theory_lra.cpp
View file

@ -2822,6 +2822,8 @@ public:
st.update("cube-success", m_solver->settings().st().m_cube_success);
st.update("arith-patches", m_solver->settings().st().m_patches);
st.update("arith-patches-success", m_solver->settings().st().m_patches_success);
st.update("arith-hnf-calls", m_solver->settings().st().m_hnf_cutter_calls);
st.update("arith-hnf-cuts", m_solver->settings().st().m_hnf_cuts);
}
};

69
src/test/lp/lp.cpp
View file

@ -3410,6 +3410,8 @@ void test_gomory_cut_1() {
g.mk_gomory_cut(t, k, expl, inf_col, row);
}
void call_hnf(general_matrix & A);
void test_hnf_m_less_than_n() {
#ifdef Z3DEBUG
general_matrix A;
@ -3432,9 +3434,7 @@ void test_hnf_m_less_than_n() {
v.push_back(mpq(1));
v.push_back(mpq(5));
A.push_row(v);
unsigned r;
mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r);
hnf<general_matrix> h(A, d, r);
call_hnf(A);
#endif
}
void test_hnf_m_greater_than_n() {
@ -3456,9 +3456,7 @@ void test_hnf_m_greater_than_n() {
v.push_back(mpq(12));
v.push_back(mpq(55));
A.push_row(v);
unsigned r;
mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r);
hnf<general_matrix> h(A, d, r);
call_hnf(A);
#endif
}
@ -3508,11 +3506,11 @@ void test_determinant() {
M[0][0] = 3; M[0][1] = -1; M[0][2] = 1;
M[1][0] = 1; M[1][1] = 0; M[1][2] = 0;
M[2][0] = 0; M[2][1] = 1; M[2][2] = 4;
unsigned r;
svector<unsigned> r;
std::cout << "M = "; M.print(std::cout, 4); endl(std::cout);
mpq d = hnf_calc::determinant_of_rectangular_matrix(M, r);
std::cout << "det M = " << d << std::endl;
std::cout << "rank = " << r << std::endl;
std::cout << "rank = " << r.size() << std::endl;
}
{
auto M = general_matrix(4, 6);
@ -3520,11 +3518,11 @@ void test_determinant() {
M[1][0] = 1; M[1][1] = 0; M[1][2] = 0; M[1][3] = 0; M[1][4] = 2; M[1][5] = 7;
M[2][0] = 0; M[2][1] = 1; M[2][2] = 4; M[2][3] = 0; M[2][4] = 2; M[2][5] = 8;
M[3][0] = 6; M[3][1] = -2; M[3][2] = 2; M[3][3] = 2; M[3][4] = 6; M[3][5] = -2;
unsigned r;
svector<unsigned> r;
std::cout << "M = "; M.print(std::cout, 4); endl(std::cout);
mpq d = hnf_calc::determinant_of_rectangular_matrix(M, r);
std::cout << "det M = " << d << std::endl;
std::cout << "rank = " << r << std::endl;
std::cout << "rank = " << r.size() << std::endl;
}
#endif
}
@ -3539,11 +3537,11 @@ void fill_general_matrix(general_matrix & M) {
M[i][j] = mpq(static_cast<int>(my_random() % 13) - 6);
}
void call_hnf(general_matrix A) {
unsigned r;
void call_hnf(general_matrix& A) {
svector<unsigned> r;
mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r);
if (r == A.row_count())
hnf<general_matrix> h(A, d, r);
A.shrink_to_rank(r);
hnf<general_matrix> h(A, d);
}
@ -3559,9 +3557,7 @@ void test_hnf_1_2() {
v.push_back(mpq(5));
v.push_back(mpq(26));
A.push_row(v);
unsigned r;
mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r);
hnf<general_matrix> h(A, d, r);
call_hnf(A);
std::cout << "test_hnf_1_2 passed" << std::endl;
}
void test_hnf_2_2() {
@ -3575,9 +3571,7 @@ void test_hnf_2_2() {
v.push_back(mpq(2));
v.push_back(mpq(11));
A.push_row(v);
unsigned r;
mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r);
hnf<general_matrix> h(A, d, r);
call_hnf(A);
std::cout << "test_hnf_2_2 passed" << std::endl;
}
@ -3676,8 +3670,41 @@ void test_hnf_5_5() {
std::cout << "test_hnf_5_5 passed" << std::endl;
}
void test_small_generated_hnf() {
std::cout << "test_small_rank_hnf" << std::endl;
general_matrix A;
vector<mpq> v;
v.push_back(mpq(5));
v.push_back(mpq(26));
A.push_row(v);
v.clear();
v.push_back(zero_of_type<mpq>());
v.push_back(zero_of_type<mpq>());
A.push_row(v);
call_hnf(A);
std::cout << "test_small_rank_hnf passed" << std::endl;
}
void test_larger_generated_hnf() {
std::cout << "test_larger_generated_rank_hnf" << std::endl;
general_matrix A;
vector<mpq> v;
v.push_back(zero_of_type<mpq>());
v.push_back(zero_of_type<mpq>());
v.push_back(zero_of_type<mpq>());
A.push_row(v);
A.push_row(v);
v.clear();
v.push_back(mpq(5));
v.push_back(mpq(26));
v.push_back(mpq(3));
A.push_row(v);
call_hnf(A);
std::cout << "test_larger_generated_rank_hnf passed" << std::endl;
}
void test_hnf() {
test_determinant();
test_larger_generated_hnf();
test_small_generated_hnf();
test_hnf_1_2();
test_hnf_3_3();
test_hnf_4_4();

1
src/util/lp/dense_matrix_def.h
View file

@ -23,7 +23,6 @@ Revision History:
#include "util/lp/numeric_pair.h"
#include "util/lp/dense_matrix.h"
namespace lp {
template <typename T> void print_vector(const vector<T> & t, std::ostream & out);
template <typename T, typename X> dense_matrix<T, X>::dense_matrix(unsigned m, unsigned n) : m_m(m), m_n(n), m_values(m * n, numeric_traits<T>::zero()) {
}

31
src/util/lp/explanation.h Normal file
View file

@ -0,0 +1,31 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Nikolaj Bjorner (nbjorner)
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
namespace lp {
struct explanation {
vector<std::pair<mpq, constraint_index>> m_explanation;
void push_justification(constraint_index j, const mpq& v) {
m_explanation.push_back(std::make_pair(v, j));
}
void push_justification(constraint_index j) {
m_explanation.push_back(std::make_pair(one_of_type<mpq>(), j));
}
};
}

65
src/util/lp/general_matrix.h
View file

@ -50,7 +50,7 @@ public:
unsigned row_count() const { return m_data.size(); }
unsigned column_count() const { return m_data[0].size(); }
unsigned column_count() const { return m_data.size() > 0? m_data[0].size() : 0; }
class ref_row {
general_matrix& m_matrix;
@ -70,11 +70,15 @@ public:
ref_row operator[](unsigned i) { return ref_row(*this, m_data[adjust_row(i)]); }
ref_row_const operator[](unsigned i) const { return ref_row_const(*this, m_data[adjust_row(i)]); }
#ifdef Z3DEBUG
void print(std::ostream & out, unsigned blanks = 0) const {
print_matrix<mpq>(m_data, out, blanks);
}
void clear() { m_data.clear(); }
void print(std::ostream & out, const char * ss) const {
std::string s(ss);
out << s;
print(out, s.size());
}
void print_submatrix(std::ostream & out, unsigned k, unsigned blanks = 0) const {
general_matrix m(row_count() - k, column_count() - k);
@ -85,9 +89,29 @@ public:
print_matrix<mpq>(m.m_data, out, blanks);
}
#endif
void clear() { m_data.clear(); }
bool row_is_initialized_correctly(const vector<mpq>& row) {
lp_assert(row.size() == column_count());
for (unsigned j = 0; j < row.size(); j ++)
lp_assert(is_zero(row[j]));
return true;
}
template <typename T>
void init_row_from_container(int i, const T & c, std::function<unsigned (unsigned)> column_fix) {
auto & row = m_data[adjust_row(i)];
lp_assert(row_is_initialized_correctly(row));
for (const auto & p : c) {
unsigned j = adjust_column(column_fix(p.var()));
row[j] = p.coeff();
}
}
void copy_column_to_indexed_vector(unsigned entering, indexed_vector<mpq> &w ) const {
lp_assert(false); //
lp_assert(false); // not implemented
}
general_matrix operator*(const general_matrix & m) const {
lp_assert(m.row_count() == column_count());
@ -191,5 +215,38 @@ public:
v.resize(n);
}
}
void shrink_to_rank(const svector<unsigned>& basis_rows) {
if (basis_rows.size() == row_count()) return;
vector<vector<mpq>> data; // todo : not efficient code
for (unsigned i : basis_rows)
data.push_back(m_data[i]);
m_data = data;
}
// used for debug only
general_matrix take_first_n_columns(unsigned n) const {
lp_assert(n <= column_count());
if (n == column_count())
return *this;
general_matrix ret(row_count(), n);
for (unsigned i = 0; i < row_count(); i++)
for (unsigned j = 0; j < n; j++)
ret[i][j] = (*this)[i][j];
return ret;
}
inline
friend vector<mpq> operator*(const vector<mpq> & f, const general_matrix& a) {
vector<mpq> r(a.column_count());
for (unsigned j = 0; j < a.column_count(); j ++) {
mpq t = zero_of_type<mpq>();
for (unsigned i = 0; i < a.row_count(); i++) {
t += f[i] * a[i][j];
}
r[j] = t;
}
return r;
}
};
}

46
src/util/lp/hnf.h
View file

@ -104,7 +104,7 @@ void extended_gcd_minimal_uv(const mpq & a, const mpq & b, mpq & d, mpq & u, mpq
template <typename M> bool prepare_pivot_for_lower_triangle(M &m, unsigned r) {
template <typename M> bool prepare_pivot_for_lower_triangle(M &m, unsigned r, svector<unsigned> & basis_rows) {
lp_assert(m.row_count() <= m.column_count());
for (unsigned i = r; i < m.row_count(); i++) {
for (unsigned j = r; j < m.column_count(); j++) {
@ -112,6 +112,7 @@ template <typename M> bool prepare_pivot_for_lower_triangle(M &m, unsigned r) {
if (i != r) {
m.transpose_rows(i, r);
}
basis_rows.push_back(i);
if (j != r) {
m.transpose_columns(j, r);
}
@ -137,12 +138,12 @@ template <typename M> void pivot_column_non_fractional(M &m, unsigned & r) {
}
// returns the rank of the matrix
template <typename M> unsigned to_lower_triangle_non_fractional(M &m) {
template <typename M> void to_lower_triangle_non_fractional(M &m, svector<unsigned> & basis_rows ) {
lp_assert(m.row_count() <= m.column_count());
unsigned i = 0;
for (; i < m.row_count() - 1; i++) {
if (!prepare_pivot_for_lower_triangle(m, i)) {
return i;
if (!prepare_pivot_for_lower_triangle(m, i, basis_rows)) {
return;
}
pivot_column_non_fractional(m, i);
}
@ -150,10 +151,10 @@ template <typename M> unsigned to_lower_triangle_non_fractional(M &m) {
// go over the last row and try to find a non-zero in the row to the right of diagonal
for (unsigned j = i; j < m.column_count(); j++) {
if (!is_zero(m[i][j])) {
return m.row_count();
basis_rows.push_back(i);
break;
}
}
return i;
}
template <typename M>
mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) {
@ -173,8 +174,8 @@ mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) {
}
// it returns "r" - the rank of the matrix and the gcd of r-minors
template <typename M> mpq determinant_of_rectangular_matrix(const M& m, unsigned & r) {
// it fills "r" - the basic rows of m
template <typename M> mpq determinant_of_rectangular_matrix(const M& m, svector<unsigned> & basis_rows) {
if (m.column_count() < m.row_count())
throw "not implemented"; // consider working with the transposed m, or create a "transposed" version if needed
// The plan is to transform m to the lower triangular form by using non-fractional Gaussian Elimination by columns.
@ -182,19 +183,18 @@ template <typename M> mpq determinant_of_rectangular_matrix(const M& m, unsigned
// m[r-1][r-1], m[r-1][r], ..., m[r-1]m[m.column_count() - 1] give the determinants of all minors of rank r.
// The gcd of these minors is the return value
auto mc = m;
r = to_lower_triangle_non_fractional(mc);
if (r == 0)
to_lower_triangle_non_fractional(mc, basis_rows);
if (basis_rows.size() == 0)
return one_of_type<mpq>();
lp_assert(!is_zero(gcd_of_row_starting_from_diagonal(mc, r - 1)));
return gcd_of_row_starting_from_diagonal(mc, r - 1);
return gcd_of_row_starting_from_diagonal(mc, basis_rows.size() - 1);
}
template <typename M> mpq determinant(const M& m) {
lp_assert(m.row_count() == m.column_count());
auto mc = m;
unsigned r;
mpq d = determinant_of_rectangular_matrix(mc, r);
return r < m.row_count() ? zero_of_type<mpq>() : d;
svector<unsigned> basis_rows;
mpq d = determinant_of_rectangular_matrix(mc, basis_rows);
return basis_rows.size() < m.row_count() ? zero_of_type<mpq>() : d;
}
} // end of namespace hnf_calc
@ -204,7 +204,7 @@ class hnf {
// fields
#ifdef Z3DEBUG
M & m_H;
M m_H;
M m_U;
M m_U_reverse;
#endif
@ -215,7 +215,7 @@ class hnf {
unsigned m_m;
unsigned m_n;
mpq m_d; // it is a positive number and a multiple of gcd of r-minors of m_A_orig, where r is the rank of m_A_orig
mpq m_r; // the rank of m_A
// we suppose that the rank of m_A is equal to row_count(), and that row_count() <= column_count(), that is m_A has the full rank
unsigned m_i;
unsigned m_j;
mpq m_R;
@ -391,7 +391,7 @@ class hnf {
void work_on_columns_less_than_i_in_the_triangle(unsigned i) {
const mpq & mii = m_H[i][i];
lp_assert(is_pos(mii));
if (is_zero(mii)) return;
for (unsigned j = 0; j < i; j++) {
const mpq & mij = m_H[i][j];
if (!is_pos(mij) && - mij < mii)
@ -444,7 +444,7 @@ class hnf {
const mpq & hij = m_H[i][j];
if (is_pos(hij))
return false;
if (- hij >= hii)
if (!is_zero(hii) && - hij >= hii)
return false;
}
@ -594,7 +594,7 @@ private:
}
public:
hnf(M & A, const mpq & d, unsigned r) :
hnf(M & A, const mpq & d) :
#ifdef Z3DEBUG
m_H(A),
#endif
@ -604,12 +604,10 @@ public:
m_m(A.row_count()),
m_n(A.column_count()),
m_d(d),
m_r(r),
m_R(m_d),
m_half_R(floor(m_R / 2))
{
lp_assert(m_m > 0 && m_n > 0);
if (is_zero(m_d))
if (m_m == 0 || m_n == 0 || is_zero(m_d))
return;
#ifdef Z3DEBUG
prepare_U_and_U_reverse();
@ -627,7 +625,7 @@ public:
#endif
}
const M & W() const { return m_W; }
};

157
src/util/lp/hnf_cutter.h
View file

@ -22,25 +22,170 @@ Revision History:
#include "util/lp/hnf.h"
#include "util/lp/general_matrix.h"
#include "util/lp/var_register.h"
#include "util/lp/lia_move.h"
#include "util/lp/explanation.h"
namespace lp {
class hnf_cutter {
var_register m_var_register;
general_matrix m_A;
vector<const lar_term*> m_terms;
var_register m_var_register;
general_matrix m_A;
vector<const lar_term*> m_terms;
vector<mpq> m_right_sides;
unsigned m_row_count;
unsigned m_column_count;
std::function<unsigned ()> m_random_next;
public:
hnf_cutter(std::function<unsigned()> random) : m_random_next(random) {}
void clear() {
m_A.clear();
m_var_register.clear();
m_terms.clear();
m_row_count = m_column_count = 0;
}
void add_term_to_A_for_hnf(const lar_term* t, const mpq &) {
void add_term(const lar_term* t, const mpq &rs) {
m_terms.push_back(t);
for (const auto &p : *t) {
m_var_register.register_user_var(p.var());
m_var_register.add_var(p.var());
}
m_right_sides.push_back(rs);
if (m_terms.size() <= m_var_register.size()) { // capture the maximal acceptable sizes
m_row_count = m_terms.size();
m_column_count = m_var_register.size();
}
}
void print(std::ostream & out) {
out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl;
}
void initialize_row(unsigned i) {
m_A.init_row_from_container(i, * m_terms[i], [this](unsigned j) { return m_var_register.add_var(j);});
}
void init_matrix_A() {
m_A = general_matrix(m_row_count, m_column_count);
// use the last suitable counts to make the number
// of rows less than or equal to the number of columns
for (unsigned i = 0; i < m_row_count; i++)
initialize_row(i);
}
// todo: as we need only one row i with non integral b[i] need to optimize later
void find_h_minus_1_b(const general_matrix& H, vector<mpq> & b) {
// the solution will be put into b
for (unsigned i = 0; i < H.row_count() ;i++) {
for (unsigned j = 0; j < i; j++) {
b[i] -= H[i][j]*b[j];
}
b[i] /= H[i][i];
// consider return from here if b[i] is not an integer and return i
}
}
vector<mpq> create_b(const svector<unsigned> & basis_rows) {
if (basis_rows.size() == m_right_sides.size())
return m_right_sides;
vector<mpq> b;
for (unsigned i : basis_rows)
b.push_back(m_right_sides[i]);
return b;
}
int find_cut_row_index(const vector<mpq> & b) {
int ret = -1;
int n = 0;
for (int i = 0; i < static_cast<int>(b.size()); i++) {
if (!is_int(b[i])) {
if (n == 0 ) {
lp_assert(ret == -1);
n = 1;
ret = i;
} else {
if (m_random_next() % (++n) == 0) {
ret = i;
}
}
}
}
return ret;
}
// fills e_i*H_minus_1
void get_ei_H_minus_1(unsigned i, const general_matrix& H, vector<mpq> & row) {
// we solve x = ei * H_min_1
// or x * H = ei
unsigned m = H.row_count();
for (unsigned k = i + 1; k < m; k++) {
row[k] = zero_of_type<mpq>();
}
row[i] = one_of_type<mpq>() / H[i][i];
for(int k = i - 1; k >= 0; k--) {
mpq t = zero_of_type<mpq>();
for (unsigned l = k + 1; l <= i; l++) {
t += H[l][k]*row[l];
}
row[k] = -t / H[k][k];
}
// test region
vector<mpq> ei(H.row_count(), zero_of_type<mpq>());
ei[i] = one_of_type<mpq>();
vector<mpq> pr = row * H;
pr.shrink(ei.size());
lp_assert(ei == pr);
// end test region
}
void fill_term(const vector<mpq> & row, lar_term& t) {
for (unsigned j = 0; j < row.size(); j++) {
if (!is_zero(row[j]))
t.add_monomial(row[j], m_var_register.local_var_to_user_var(j));
}
}
lia_move create_cut(lar_term& t, mpq& k, explanation& ex, bool & upper) {
init_matrix_A();
svector<unsigned> basis_rows;
mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows);
if (basis_rows.size() < m_A.row_count())
m_A.shrink_to_rank(basis_rows);
hnf<general_matrix> h(m_A, d);
// general_matrix A_orig = m_A;
vector<mpq> b = create_b(basis_rows);
vector<mpq> bcopy = b;
find_h_minus_1_b(h.W(), b);
lp_assert(bcopy == h.W().take_first_n_columns(b.size()) * b);
int cut_row = find_cut_row_index(b);
if (cut_row == -1) {
return lia_move::undef;
}
// test region
/*
general_matrix U(m_A.column_count());
vector<mpq> rt(m_A.column_count());
for (unsigned i = 0; i < U.row_count(); i++) {
get_ei_H_minus_1(i, h.W(), rt);
vector<mpq> ft = rt * A_orig;
for (unsigned j = 0; j < ft.size(); j++)
U[i][j] = ft[j];
}
std::cout << "U reverse = "; U.print(std::cout, 12); std::cout << std::endl;
*/
// end test region
vector<mpq> row(m_A.column_count());
get_ei_H_minus_1(cut_row, h.W(), row);
vector<mpq> f = row * m_A;
fill_term(f, t);
k = floor(b[cut_row]);
upper = true;
return lia_move::cut;
}
};
}

12
src/util/lp/indexed_vector.cpp
View file

@ -42,12 +42,12 @@ template void lp::indexed_vector<double>::print(std::basic_ostream<char,struct s
template void lp::indexed_vector<lp::numeric_pair<lp::mpq> >::print(std::ostream&);
#endif
}
template void lp::print_vector<double>(vector<double> const&, std::ostream&);
template void lp::print_vector<unsigned int>(vector<unsigned int> const&, std::ostream&);
template void lp::print_vector<std::string>(vector<std::string> const&, std::ostream&);
template void lp::print_vector<lp::numeric_pair<lp::mpq> >(vector<lp::numeric_pair<lp::mpq>> const&, std::ostream&);
// template void lp::print_vector<double, vectro>(vector<double> const&, std::ostream&);
// template void lp::print_vector<unsigned int>(vector<unsigned int> const&, std::ostream&);
// template void lp::print_vector<std::string>(vector<std::string> const&, std::ostream&);
// template void lp::print_vector<lp::numeric_pair<lp::mpq> >(vector<lp::numeric_pair<lp::mpq>> const&, std::ostream&);
template void lp::indexed_vector<double>::resize(unsigned int);
template void lp::print_vector< lp::mpq>(vector< lp::mpq> const &, std::basic_ostream<char, std::char_traits<char> > &);
template void lp::print_vector<std::pair<lp::mpq, unsigned int> >(vector<std::pair<lp::mpq, unsigned int>> const&, std::ostream&);
// template void lp::print_vector< lp::mpq>(vector< lp::mpq> const &, std::basic_ostream<char, std::char_traits<char> > &);
// template void lp::print_vector<std::pair<lp::mpq, unsigned int> >(vector<std::pair<lp::mpq, unsigned int>> const&, std::ostream&);
template void lp::indexed_vector<lp::numeric_pair<lp::mpq> >::erase_from_index(unsigned int);

2
src/util/lp/indexed_vector.h
View file

@ -28,8 +28,6 @@ Revision History:
#include <unordered_set>
namespace lp {
template <typename T> void print_vector(const vector<T> & t, std::ostream & out);
template <typename T> void print_vector(const buffer<T> & t, std::ostream & out);
template <typename T> void print_sparse_vector(const vector<T> & t, std::ostream & out);
void print_vector_as_doubles(const vector<mpq> & t, std::ostream & out);

15
src/util/lp/indexed_vector_def.h
View file

@ -22,21 +22,6 @@ Revision History:
#include "util/lp/lp_settings.h"
namespace lp {
template <typename T>
void print_vector(const vector<T> & t, std::ostream & out) {
for (unsigned i = 0; i < t.size(); i++)
out << t[i] << " ";
out << std::endl;
}
template <typename T>
void print_vector(const buffer<T> & t, std::ostream & out) {
for (unsigned i = 0; i < t.size(); i++)
out << t[i] << " ";
out << std::endl;
}
template <typename T>
void print_sparse_vector(const vector<T> & t, std::ostream & out) {
for (unsigned i = 0; i < t.size(); i++) {

63
src/util/lp/int_solver.cpp
View file

@ -272,14 +272,20 @@ void int_solver::gomory_cut_adjust_t_and_k(vector<std::pair<mpq, unsigned>> & po
bool int_solver::current_solution_is_inf_on_cut() const {
const auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x;
impq v = m_t->apply(x);
TRACE(
"current_solution_is_inf_on_cut", tout << "v = " << v << " k = " << (*m_k) << std::endl;
if (v <=(*m_k)) {
tout << "v <= k - it should not happen!\n";
}
mpq sign = !(*m_upper) ? one_of_type<mpq>() : -one_of_type<mpq>();
TRACE("current_solution_is_inf_on_cut",
if (is_pos(sign)) {
tout << "v = " << v << " k = " << (*m_k) << std::endl;
if (v <=(*m_k)) {
tout << "v <= k - it should not happen!\n";
}
} else {
if (v >= (*m_k)) {
tout << "v > k - it should not happen!\n";
}
}
);
return v > (*m_k);
return v * sign > (*m_k) * sign;
}
void int_solver::adjust_term_and_k_for_some_ints_case_gomory(mpq &lcm_den) {
@ -635,33 +641,44 @@ lia_move int_solver::gomory_cut() {
}
void int_solver::try_add_term_to_A_for_hnf(unsigned i) {
bool int_solver::try_add_term_to_A_for_hnf(unsigned i) {
mpq rs;
const lar_term* t = m_lar_solver->terms()[i];
for (const auto & p : *t) {
if (!is_int(p.var()))
return; // todo : the mix case!
return false; // todo : the mix case!
}
if (m_lar_solver->get_equality_for_term_on_corrent_x(i, rs)) {
m_hnf_cutter.add_term(t, rs);
return true;
} else {
return false;
}
if (!m_lar_solver->get_equality_for_term_on_corrent_x(i, rs))
return;
m_hnf_cutter.add_term_to_A_for_hnf(t, rs);
}
bool int_solver::hnf_matrix_is_empty() const { return true; }
bool int_solver::prepare_matrix_A_for_hnf_cut() {
m_hnf_cutter.clear();
for (unsigned i = 0; i < m_lar_solver->terms().size(); i++)
try_add_term_to_A_for_hnf(i);
m_hnf_cutter.print(std::cout);
return ! hnf_matrix_is_empty();
for (unsigned i = 0; i < m_lar_solver->terms().size(); i++) {
bool r = try_add_term_to_A_for_hnf(i);
if (!r && settings().hnf_cutter_exit_if_x_is_not_on_bound_or_mixed )
return false;
}
return true;
}
lia_move int_solver::make_hnf_cut() {
if( !prepare_matrix_A_for_hnf_cut())
if (!prepare_matrix_A_for_hnf_cut()) {
return lia_move::undef;
return lia_move::undef;
}
settings().st().m_hnf_cutter_calls++;
lia_move r = m_hnf_cutter.create_cut(*m_t, *m_k, *m_ex, *m_upper);
CTRACE("hnf_cut", r == lia_move::cut, tout<< "cut:"; m_lar_solver->print_term(*m_t, tout); tout << " <= " << *m_k << std::endl;);
if (r == lia_move::cut)
settings().st().m_hnf_cuts++;
return r;
}
lia_move int_solver::hnf_cut() {
@ -1009,10 +1026,12 @@ int_solver::int_solver(lar_solver* lar_slv) :
m_lar_solver(lar_slv),
m_branch_cut_counter(0),
m_chase_cut_solver([this](unsigned j) {return m_lar_solver->get_column_name(j);},
[this](unsigned j, std::ostream &o) {m_lar_solver->print_constraint(j, o);},
[this]() {return m_lar_solver->A_r().column_count();},
[this](unsigned j) {return get_value(j);},
settings()) {
[this](unsigned j, std::ostream &o) {m_lar_solver->print_constraint(j, o);},
[this]() {return m_lar_solver->A_r().column_count();},
[this](unsigned j) {return get_value(j);},
settings()),
m_hnf_cutter([this](){ return settings().random_next(); })
{
m_lar_solver->set_int_solver(this);
}

23
src/util/lp/int_solver.h
View file

@ -25,30 +25,15 @@ Revision History:
#include "util/lp/chase_cut_solver.h"
#include "util/lp/lar_constraints.h"
#include "util/lp/hnf_cutter.h"
#include "util/lp/lia_move.h"
#include "util/lp/explanation.h"
namespace lp {
class lar_solver;
template <typename T, typename X>
struct lp_constraint;
enum class lia_move {
sat,
branch,
cut,
conflict,
continue_with_check,
undef,
unsat
};
struct explanation {
vector<std::pair<mpq, constraint_index>> m_explanation;
void push_justification(constraint_index j, const mpq& v) {
m_explanation.push_back(std::make_pair(v, j));
}
void push_justification(constraint_index j) {
m_explanation.push_back(std::make_pair(one_of_type<mpq>(), j));
}
};
class int_solver {
public:
@ -182,6 +167,6 @@ public:
lia_move make_hnf_cut();
bool prepare_matrix_A_for_hnf_cut();
bool hnf_matrix_is_empty() const;
void try_add_term_to_A_for_hnf(unsigned term_index);
bool try_add_term_to_A_for_hnf(unsigned term_index);
};
}

31
src/util/lp/lia_move.h Normal file
View file

@ -0,0 +1,31 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Nikolaj Bjorner (nbjorner)
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
namespace lp {
enum class lia_move {
sat,
branch,
cut,
conflict,
continue_with_check,
undef,
unsat
};
}

21
src/util/lp/lp_settings.h
View file

@ -112,6 +112,8 @@ struct stats {
unsigned m_cube_success;
unsigned m_patches;
unsigned m_patches_success;
unsigned m_hnf_cutter_calls;
unsigned m_hnf_cuts;
stats() { reset(); }
void reset() { memset(this, 0, sizeof(*this)); }
};
@ -141,22 +143,22 @@ private:
random_gen m_rand;
public:
unsigned reps_in_scaler;
unsigned reps_in_scaler;
// when the absolute value of an element is less than pivot_epsilon
// in pivoting, we treat it as a zero
double pivot_epsilon;
double pivot_epsilon;
// see Chatal, page 115
double positive_price_epsilon;
double positive_price_epsilon;
// a quatation "if some choice of the entering vairable leads to an eta matrix
// whose diagonal element in the eta column is less than e2 (entering_diag_epsilon) in magnitude, the this choice is rejected ...
double entering_diag_epsilon;
int c_partial_pivoting; // this is the constant c from page 410
unsigned depth_of_rook_search;
bool using_partial_pivoting;
double entering_diag_epsilon;
int c_partial_pivoting; // this is the constant c from page 410
unsigned depth_of_rook_search;
bool using_partial_pivoting;
// dissertation of Achim Koberstein
// if Bx - b is different at any component more that refactor_epsilon then we refactor
double refactor_tolerance;
double pivot_tolerance;
double refactor_tolerance;
double pivot_tolerance;
double zero_tolerance;
double drop_tolerance;
double tolerance_for_artificials;
@ -175,6 +177,7 @@ public:
double dual_feasibility_tolerance; // // page 71 of the PhD thesis of Achim Koberstein
double primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein
double relative_primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein
bool hnf_cutter_exit_if_x_is_not_on_bound_or_mixed = true;
bool m_bound_propagation;

7
src/util/lp/lp_utils.h
View file

@ -22,6 +22,13 @@ Revision History:
#include "util/lp/numeric_pair.h"
#include "util/debug.h"
#include <unordered_map>
template <typename C>
void print_vector(const C & t, std::ostream & out) {
for (const auto & p : t)
out << p << " ";
out << std::endl;
}
template <typename A, typename B>
bool try_get_value(const std::unordered_map<A,B> & map, const A& key, B & val) {
const auto it = map.find(key);

48
src/util/lp/var_register.h Normal file
View file

@ -0,0 +1,48 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
namespace lp {
class var_register {
svector<unsigned> m_local_vars_to_external;
std::unordered_map<unsigned, unsigned> m_external_vars_to_local;
public:
unsigned add_var(unsigned user_var) {
auto t = m_external_vars_to_local.find(user_var);
if (t != m_external_vars_to_local.end()) {
return t->second;
}
unsigned ret = size();
m_external_vars_to_local[user_var] = ret;
m_local_vars_to_external.push_back(user_var);
return ret;
}
unsigned local_var_to_user_var(unsigned local_var) const {
return m_local_vars_to_external[local_var];
}
unsigned size() const {
return m_local_vars_to_external.size();
}
void clear() {
m_local_vars_to_external.clear();
m_external_vars_to_local.clear();
}
};
}

18
src/util/vector.h
View file

@ -188,6 +188,24 @@ public:
m_data = nullptr;
}
bool operator==(vector const & other) const {
if (this == &other) {
return true;
}
if (size() != other.size())
return false;
for (unsigned i = 0; i < size(); i++) {
if ((*this)[i] != other[i])
return false;
}
return true;
}
bool operator!=(vector const & other) const {
return !(*this == other);
}
vector & operator=(vector const & source) {
if (this == &source) {
return *this;