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

clean up int_solver

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

add a diagnostic method

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

white space change

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

cleanup in int_solver

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

some cleanup

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

remove m_became_zeros

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

start cut_solver, work on disjoint_intervals

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

start cut_solver, work on disjoint_intervals

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

work on disjoint_intervals

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

work on disjoint_intervals

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

work on disjoint_intervals

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

work on disjoint_intervals

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

work on disjoint_intervals

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

workin on disjoint_intervals

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

working on disjoint_intervals

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

working on disjoint_intervals

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

working on disjoint_intervals

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

work on disjoint_intervals

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

fix bugs in disjoint_intervals

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

fix bugs in gisjoint_intervals

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

fix bugs in gisjoint_intervals

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

fix bugs in disjoint_intervals

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

fix bugs in disjoint_intervals

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

fix bugs is disjoint intervals

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

bug fixes in disjoint_intervals

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

disjoint_intervals passes the test

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

test disjoint_intervals push(), pop()

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

cut_solver

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

work on cut_solver

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2017-08-19 20:38:49 -07:00
parent db8f01894f
commit 58ca4518e5
28 changed files with 2469 additions and 1307 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
src/smt/smt_context.cpp
View file

@ -3763,6 +3763,7 @@ namespace smt {
}
m_stats.m_num_final_checks++;
TRACE("final_check_stats", tout << "m_stats.m_num_final_checks = " << m_stats.m_num_final_checks << "\n";);
final_check_status ok = m_qmanager->final_check_eh(false);
if (ok != FC_DONE)

3
src/smt/theory_lra.cpp
View file

@ -1202,7 +1202,6 @@ namespace smt {
}
final_check_status final_check_eh() {
TRACE("lar_solver", tout << "ddd=" <<++lp::lp_settings::ddd << std::endl;);
m_use_nra_model = false;
lbool is_sat = l_true;
if (m_solver->get_status() != lp::lp_status::OPTIMAL) {
@ -2450,7 +2449,7 @@ namespace smt {
m_nra->am().set(r, 0);
while (!m_todo_terms.empty()) {
rational wcoeff = m_todo_terms.back().second;
lp::var_index wi = m_todo_terms.back().first; // todo : got a warning "wi is not used"
// lp::var_index wi = m_todo_terms.back().first; // todo : got a warning "wi is not used"
m_todo_terms.pop_back();
lp::lar_term const& term = m_solver->get_term(vi);
scoped_anum r1(m_nra->am());

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

@ -47,13 +47,16 @@ Revision History:
#include "util/lp/stacked_unordered_set.h"
#include "util/lp/int_set.h"
#include "util/stopwatch.h"
#include "util/lp/integer_domain.h"
#include "util/lp/stacked_map.h"
#include <cstdlib>
namespace lp {
unsigned seed = 1;
random_gen g_rand;
static unsigned my_random() {
return g_rand();
}
random_gen g_rand;
static unsigned my_random() {
return g_rand();
}
struct simple_column_namer:public column_namer
{
std::string get_column_name(unsigned j) const override {
@ -66,7 +69,7 @@ template <typename T, typename X>
void test_matrix(sparse_matrix<T, X> & a) {
auto m = a.dimension();
// copy a to b in the reversed order
// copy a to b in the reversed order
sparse_matrix<T, X> b(m);
std::cout << "copy b to a"<< std::endl;
for (int row = m - 1; row >= 0; row--)
@ -198,9 +201,9 @@ void init_non_basic_part_of_basis_heading(vector<int> & basis_heading, vector<un
non_basic_columns.clear();
for (int j = basis_heading.size(); j--;){
if (basis_heading[j] < 0) {
non_basic_columns.push_back(j);
// the index of column j in m_nbasis is (- basis_heading[j] - 1)
basis_heading[j] = - static_cast<int>(non_basic_columns.size());
non_basic_columns.push_back(j);
// the index of column j in m_nbasis is (- basis_heading[j] - 1)
basis_heading[j] = - static_cast<int>(non_basic_columns.size());
}
}
}
@ -554,7 +557,7 @@ void test_lp_0() {
vector<unsigned> nbasis;
vector<int> heading;
lp_primal_core_solver<double, double> lpsolver(m_, b, x_star, basis, nbasis, heading, costs, column_types, upper_bound_values, settings, cn);
lp_primal_core_solver<double, double> lpsolver(m_, b, x_star, basis, nbasis, heading, costs, column_types, upper_bound_values, settings, cn);
lpsolver.solve();
}
@ -674,7 +677,7 @@ void test_swap_rows(sparse_matrix<T, X>& m, unsigned i0, unsigned i1){
for (unsigned i = 0; i < m.dimension(); i++)
for (unsigned j = 0; j < m.dimension(); j++) {
mcopy(i, j)= m(i, j);
}
}
std::cout << "swapping rows "<< i0 << "," << i1 << std::endl;
m.swap_rows(i0, i1);
@ -690,7 +693,7 @@ void test_swap_columns(sparse_matrix<T, X>& m, unsigned i0, unsigned i1){
for (unsigned i = 0; i < m.dimension(); i++)
for (unsigned j = 0; j < m.dimension(); j++) {
mcopy(i, j)= m(i, j);
}
}
m.swap_columns(i0, i1);
for (unsigned j = 0; j < m.dimension(); j++) {
@ -722,7 +725,7 @@ void test_pivot_like_swaps_and_pivot(){
sparse_matrix<double, double> m(10);
fill_matrix(m);
// print_matrix(m);
// pivot at 2,7
// pivot at 2,7
m.swap_columns(0, 7);
// print_matrix(m);
m.swap_rows(2, 0);
@ -732,7 +735,7 @@ void test_pivot_like_swaps_and_pivot(){
}
// print_matrix(m);
// say pivot at 3,4
// say pivot at 3,4
m.swap_columns(1, 4);
// print_matrix(m);
m.swap_rows(1, 3);
@ -790,7 +793,7 @@ void test_swap_rows() {
// print_matrix(m);
test_swap_rows(m, 0, 7);
// go over some corner cases
// go over some corner cases
sparse_matrix<double, double> m0(2);
test_swap_rows(m0, 0, 1);
m0(0, 0) = 3;
@ -942,7 +945,7 @@ void test_swap_columns() {
test_swap_columns(m, 0, 7);
// go over some corner cases
// go over some corner cases
sparse_matrix<double, double> m0(2);
test_swap_columns(m0, 0, 1);
m0(0, 0) = 3;
@ -1042,7 +1045,7 @@ void test_apply_reverse_from_right_to_perm(permutation_matrix<double, double> &
auto rs = pclone * rev;
lp_assert(p == rs)
#endif
}
}
void test_apply_reverse_from_right() {
auto vec = vector_of_permutaions();
@ -1818,7 +1821,7 @@ std::unordered_map<std::string, double> * get_solution_from_glpsol_output(std::s
}
auto split = string_split(s, " \t", false);
if (split.size() == 0) {
return ret;
return ret;
}
lp_assert(split.size() > 3);
@ -1880,6 +1883,8 @@ void test_replace_column() {
void setup_args_parser(argument_parser & parser) {
parser.add_option_with_help_string("-dji", "test integer_domain");
parser.add_option_with_help_string("-cs", "test cut_solver");
parser.add_option_with_help_string("-xyz_sample", "run a small interactive scenario");
parser.add_option_with_after_string_with_help("--density", "the percentage of non-zeroes in the matrix below which it is not dense");
parser.add_option_with_after_string_with_help("--harris_toler", "harris tolerance");
@ -2058,15 +2063,15 @@ void test_stacked() {
}
char * find_home_dir() {
#ifdef _WINDOWS
#else
#ifdef _WINDOWS
#else
char * home_dir = getenv("HOME");
if (home_dir == nullptr) {
if (home_dir == nullptr) {
std::cout << "cannot find home directory" << std::endl;
return nullptr;
}
#endif
return nullptr;
#endif
return nullptr;
}
@ -2088,8 +2093,8 @@ struct mem_cpy_place_holder {
void finalize(unsigned ret) {
/*
finalize_util_module();
finalize_numerics_module();
finalize_util_module();
finalize_numerics_module();
*/
// return ret;
}
@ -2227,7 +2232,7 @@ bool values_are_one_percent_close(double a, double b) {
// returns true if both are optimal
void compare_costs(std::string glpk_out_file_name,
std::string lp_out_file_name,
std::string lp_out_file_name,
unsigned & successes,
unsigned & failures) {
double a = get_glpk_cost(glpk_out_file_name);
@ -2306,78 +2311,78 @@ void process_test_file(std::string test_dir, std::string test_file_name, argumen
}
/*
int my_readdir(DIR *dirp, struct dirent *
#ifndef LEAN_WINDOWS
entry
#endif
, struct dirent **result) {
#ifdef LEAN_WINDOWS
*result = readdir(dirp); // NOLINT
return *result != nullptr? 0 : 1;
#else
return readdir_r(dirp, entry, result);
#endif
}
#ifndef LEAN_WINDOWS
entry
#endif
, struct dirent **result) {
#ifdef LEAN_WINDOWS
*result = readdir(dirp); // NOLINT
return *result != nullptr? 0 : 1;
#else
return readdir_r(dirp, entry, result);
#endif
}
*/
/*
vector<std::pair<std::string, int>> get_file_list_of_dir(std::string test_file_dir) {
DIR *dir;
if ((dir = opendir(test_file_dir.c_str())) == nullptr) {
std::cout << "Cannot open directory " << test_file_dir << std::endl;
throw 0;
}
vector<std::pair<std::string, int>> ret;
struct dirent entry;
struct dirent* result;
int return_code;
for (return_code = my_readdir(dir, &entry, &result);
#ifndef LEAN_WINDOWS
result != nullptr &&
#endif
return_code == 0;
return_code = my_readdir(dir, &entry, &result)) {
DIR *tmp_dp = opendir(result->d_name);
struct stat file_record;
if (tmp_dp == nullptr) {
std::string s = test_file_dir+ "/" + result->d_name;
int stat_ret = stat(s.c_str(), & file_record);
if (stat_ret!= -1) {
ret.push_back(make_pair(result->d_name, file_record.st_size));
} else {
perror("stat");
exit(1);
}
} else {
closedir(tmp_dp);
}
}
closedir(dir);
return ret;
}
vector<std::pair<std::string, int>> get_file_list_of_dir(std::string test_file_dir) {
DIR *dir;
if ((dir = opendir(test_file_dir.c_str())) == nullptr) {
std::cout << "Cannot open directory " << test_file_dir << std::endl;
throw 0;
}
vector<std::pair<std::string, int>> ret;
struct dirent entry;
struct dirent* result;
int return_code;
for (return_code = my_readdir(dir, &entry, &result);
#ifndef LEAN_WINDOWS
result != nullptr &&
#endif
return_code == 0;
return_code = my_readdir(dir, &entry, &result)) {
DIR *tmp_dp = opendir(result->d_name);
struct stat file_record;
if (tmp_dp == nullptr) {
std::string s = test_file_dir+ "/" + result->d_name;
int stat_ret = stat(s.c_str(), & file_record);
if (stat_ret!= -1) {
ret.push_back(make_pair(result->d_name, file_record.st_size));
} else {
perror("stat");
exit(1);
}
} else {
closedir(tmp_dp);
}
}
closedir(dir);
return ret;
}
*/
/*
struct file_size_comp {
unordered_map<std::string, int>& m_file_sizes;
file_size_comp(unordered_map<std::string, int>& fs) :m_file_sizes(fs) {}
int operator()(std::string a, std::string b) {
std::cout << m_file_sizes.size() << std::endl;
std::cout << a << std::endl;
std::cout << b << std::endl;
struct file_size_comp {
unordered_map<std::string, int>& m_file_sizes;
file_size_comp(unordered_map<std::string, int>& fs) :m_file_sizes(fs) {}
int operator()(std::string a, std::string b) {
std::cout << m_file_sizes.size() << std::endl;
std::cout << a << std::endl;
std::cout << b << std::endl;
auto ls = m_file_sizes.find(a);
std::cout << "fa" << std::endl;
auto rs = m_file_sizes.find(b);
std::cout << "fb" << std::endl;
if (ls != m_file_sizes.end() && rs != m_file_sizes.end()) {
std::cout << "fc " << std::endl;
int r = (*ls < *rs? -1: (*ls > *rs)? 1 : 0);
std::cout << "calc r " << std::endl;
return r;
} else {
std::cout << "sc " << std::endl;
return 0;
}
}
};
auto ls = m_file_sizes.find(a);
std::cout << "fa" << std::endl;
auto rs = m_file_sizes.find(b);
std::cout << "fb" << std::endl;
if (ls != m_file_sizes.end() && rs != m_file_sizes.end()) {
std::cout << "fc " << std::endl;
int r = (*ls < *rs? -1: (*ls > *rs)? 1 : 0);
std::cout << "calc r " << std::endl;
return r;
} else {
std::cout << "sc " << std::endl;
return 0;
}
}
};
*/
struct sort_pred {
@ -2389,26 +2394,26 @@ struct sort_pred {
void test_files_from_directory(std::string test_file_dir, argument_parser & args_parser) {
/*
std::cout << "loading files from directory \"" << test_file_dir << "\"" << std::endl;
std::string out_dir = args_parser.get_option_value("--out_dir");
if (out_dir.size() == 0) {
out_dir = "/tmp/test";
}
DIR *out_dir_p = opendir(out_dir.c_str());
if (out_dir_p == nullptr) {
std::cout << "Cannot open output directory \"" << out_dir << "\"" << std::endl;
return;
}
closedir(out_dir_p);
vector<std::pair<std::string, int>> files = get_file_list_of_dir(test_file_dir);
std::sort(files.begin(), files.end(), sort_pred());
unsigned max_iters, time_limit;
get_time_limit_and_max_iters_from_parser(args_parser, time_limit, max_iters);
unsigned successes = 0, failures = 0, inconclusives = 0;
for (auto & t : files) {
process_test_file(test_file_dir, t.first, args_parser, out_dir, max_iters, time_limit, successes, failures, inconclusives);
}
std::cout << "comparing with glpk: successes " << successes << ", failures " << failures << ", inconclusives " << inconclusives << std::endl;
std::cout << "loading files from directory \"" << test_file_dir << "\"" << std::endl;
std::string out_dir = args_parser.get_option_value("--out_dir");
if (out_dir.size() == 0) {
out_dir = "/tmp/test";
}
DIR *out_dir_p = opendir(out_dir.c_str());
if (out_dir_p == nullptr) {
std::cout << "Cannot open output directory \"" << out_dir << "\"" << std::endl;
return;
}
closedir(out_dir_p);
vector<std::pair<std::string, int>> files = get_file_list_of_dir(test_file_dir);
std::sort(files.begin(), files.end(), sort_pred());
unsigned max_iters, time_limit;
get_time_limit_and_max_iters_from_parser(args_parser, time_limit, max_iters);
unsigned successes = 0, failures = 0, inconclusives = 0;
for (auto & t : files) {
process_test_file(test_file_dir, t.first, args_parser, out_dir, max_iters, time_limit, successes, failures, inconclusives);
}
std::cout << "comparing with glpk: successes " << successes << ", failures " << failures << ", inconclusives " << inconclusives << std::endl;
*/
}
@ -2654,7 +2659,7 @@ void check_lu_from_file(std::string lufile_name) {
lp_settings settings;
vector<unsigned> non_basic_columns;
lu<double, double> lsuhl(A, basis, settings);
indexed_vector<double> d(A.row_count());
indexed_vector<double> d(A.row_count());
unsigned entering = 26;
lsuhl.solve_Bd(entering, d, v);
#ifdef Z3DEBUG
@ -2765,21 +2770,21 @@ void test_evidence_for_total_inf_simple(argument_parser & args_parser) {
}
void test_bound_propagation_one_small_sample1() {
/*
(<= (+ a (* (- 1.0) b)) 0.0)
(<= (+ b (* (- 1.0) x_13)) 0.0)
--> (<= (+ a (* (- 1.0) c)) 0.0)
(<= (+ a (* (- 1.0) b)) 0.0)
(<= (+ b (* (- 1.0) x_13)) 0.0)
--> (<= (+ a (* (- 1.0) c)) 0.0)
the inequality on (<= a c) is obtained from a triangle inequality (<= a b) (<= b c).
If b becomes basic variable, then it is likely the old solver ends up with a row that implies (<= a c).
a - b <= 0.0
b - c <= 0.0
the inequality on (<= a c) is obtained from a triangle inequality (<= a b) (<= b c).
If b becomes basic variable, then it is likely the old solver ends up with a row that implies (<= a c).
a - b <= 0.0
b - c <= 0.0
got to get a <= c
got to get a <= c
*/
std::function<bool (unsigned, bool, bool, const mpq & )> bound_is_relevant =
[&](unsigned j, bool is_low_bound, bool strict, const rational& bound_val) {
return true;
};
return true;
};
lar_solver ls;
unsigned a = ls.add_var(0, false);
unsigned b = ls.add_var(1, false);
@ -2815,34 +2820,34 @@ void test_bound_propagation_one_small_samples() {
test_bound_propagation_one_small_sample1();
/*
(>= x_46 0.0)
(<= x_29 0.0)
(not (<= x_68 0.0))
(<= (+ (* (/ 1001.0 1998.0) x_10) (* (- 1.0) x_151) x_68) (- (/ 1001.0 999.0)))
(<= (+ (* (/ 1001.0 999.0) x_9)
(* (- 1.0) x_152)
(* (/ 1001.0 999.0) x_151)
(* (/ 1001.0 999.0) x_68))
(- (/ 1502501.0 999000.0)))
(not (<= (+ (* (/ 999.0 2.0) x_10) (* (- 1.0) x_152) (* (- (/ 999.0 2.0)) x_151))
(/ 1001.0 2.0)))
(not (<= x_153 0.0))z
(>= (+ x_9 (* (- (/ 1001.0 999.0)) x_10) (* (- 1.0) x_153) (* (- 1.0) x_68))
(/ 5003.0 1998.0))
--> (not (<= (+ x_10 x_46 (* (- 1.0) x_29)) 0.0))
(<= x_29 0.0)
(not (<= x_68 0.0))
(<= (+ (* (/ 1001.0 1998.0) x_10) (* (- 1.0) x_151) x_68) (- (/ 1001.0 999.0)))
(<= (+ (* (/ 1001.0 999.0) x_9)
(* (- 1.0) x_152)
(* (/ 1001.0 999.0) x_151)
(* (/ 1001.0 999.0) x_68))
(- (/ 1502501.0 999000.0)))
(not (<= (+ (* (/ 999.0 2.0) x_10) (* (- 1.0) x_152) (* (- (/ 999.0 2.0)) x_151))
(/ 1001.0 2.0)))
(not (<= x_153 0.0))z
(>= (+ x_9 (* (- (/ 1001.0 999.0)) x_10) (* (- 1.0) x_153) (* (- 1.0) x_68))
(/ 5003.0 1998.0))
--> (not (<= (+ x_10 x_46 (* (- 1.0) x_29)) 0.0))
and
and
(<= (+ a (* (- 1.0) b)) 0.0)
(<= (+ b (* (- 1.0) x_13)) 0.0)
--> (<= (+ a (* (- 1.0) x_13)) 0.0)
(<= (+ a (* (- 1.0) b)) 0.0)
(<= (+ b (* (- 1.0) x_13)) 0.0)
--> (<= (+ a (* (- 1.0) x_13)) 0.0)
In the first case, there typically are no atomic formulas for bounding x_10. So there is never some
basic lemma of the form (>= x46 0), (<= x29 0), (>= x10 0) -> (not (<= (+ x10 x46 (- x29)) 0)).
Instead the bound on x_10 falls out from a bigger blob of constraints.
In the first case, there typically are no atomic formulas for bounding x_10. So there is never some
basic lemma of the form (>= x46 0), (<= x29 0), (>= x10 0) -> (not (<= (+ x10 x46 (- x29)) 0)).
Instead the bound on x_10 falls out from a bigger blob of constraints.
In the second case, the inequality on (<= x19 x13) is obtained from a triangle inequality (<= x19 x9) (<= x9 x13).
If x9 becomes basic variable, then it is likely the old solver ends up with a row that implies (<= x19 x13).
*/
In the second case, the inequality on (<= x19 x13) is obtained from a triangle inequality (<= x19 x9) (<= x9 x13).
If x9 becomes basic variable, then it is likely the old solver ends up with a row that implies (<= x19 x13).
*/
}
void test_bound_propagation_one_row() {
lar_solver ls;
@ -3083,11 +3088,134 @@ void test_rationals() {
std::cout << T_to_string(r) << std::endl;
}
void test_lp_local(int argn, char**argv) {
std::cout << "resize\n";
vector<mpq> r;
r.resize(1);
void get_random_interval(bool& neg_inf, bool& pos_inf, int& x, int &y) {
int i = my_random() % 10;
if (i == 0) {
neg_inf = true;
} else {
neg_inf = false;
x = my_random() % 100;
}
i = my_random() % 10;
if (i == 0) {
pos_inf = true;
} else {
pos_inf = false;
if (!neg_inf) {
y = x + my_random() % (101 - x);
lp_assert(y >= x);
}
else {
y = my_random() % 100;
}
}
lp_assert((neg_inf || (0 <= x && x <= 100)) && (pos_inf || (0 <= y && y <= 100)));
}
void test_integer_domain_intersection(integer_domain<int> & d) {
int x, y; bool neg_inf, pos_inf;
get_random_interval(neg_inf, pos_inf, x, y);
if (neg_inf) {
if (!pos_inf) {
d.intersect_with_upper_bound(y);
}
}
else if (pos_inf)
d.intersect_with_lower_bound(x);
else
d.intersect_with_interval(x, y);
}
void test_integer_domain_union(integer_domain<int> & d) {
int x, y; bool neg_inf, pos_inf;
get_random_interval(neg_inf, pos_inf, x, y);
if (neg_inf) {
if (!pos_inf) {
d.unite_with_interval_neg_inf_x(y);
}
else
d.init_to_contain_all();
}
else if (pos_inf)
d.unite_with_interval_x_pos_inf(x);
else
d.unite_with_interval(x, y);
lp_assert(d.is_correct());
}
void test_integer_domain_randomly(integer_domain<int> & d) {
int i = my_random() % 10;
if (i == 0)
test_integer_domain_intersection(d);
else
test_integer_domain_union(d);
}
void test_integer_domain() {
integer_domain<int> d;
vector<integer_domain<int>> stack;
for (int i = 0; i < 10000; i++) {
test_integer_domain_randomly(d);
stack.push_back(d);
d.push();
if (i > 0 && i%100 == 0) {
if (stack.size() == 0) continue;
unsigned k = my_random() % stack.size();
if (k == 0)
k = 1;
d.pop(k);
d.restore_domain();
for (unsigned j = 0; j + 1 < k; j++) {
stack.pop_back();
}
std::cout<<"comparing i = " << i << std::endl;
lp_assert(d == *stack.rbegin());
stack.pop_back();
}
//d.print(std::cout);
}
}
void test_cut_solver() {
cut_solver<int> cs([](unsigned i)
{
if (i == 0) return std::string("x");
if (i == 1) return std::string("y");
return std::to_string(i);
});
vector<std::pair<int, unsigned>> term;
unsigned x = 0;
unsigned y = 1;
term.push_back(std::make_pair(2, x));
term.push_back(std::make_pair(-3, y));
unsigned ineq_index = cs.add_ineq(term, mpq(3));
cs.print_ineq(ineq_index, std::cout);
mpq l;
auto ineq = cs.m_ineqs[ineq_index];
cs.add_lower_bound_for_user_var(x, 1);
cs.add_lower_bound_for_user_var(y, 1);
bool has_lower = cs.lower(ineq.m_poly, l);
if (has_lower) {
std::cout << "lower = " << l << std::endl;
} else {
std::cout << "no lower" << std::endl;
}
cs.add_upper_bound_for_user_var(y, 1);
has_lower = cs.lower(ineq.m_poly, l);
if (has_lower) {
std::cout << "lower = " << l << std::endl;
} else {
std::cout << "no lower" << std::endl;
}
}
void test_lp_local(int argn, char**argv) {
// initialize_util_module();
// initialize_numerics_module();
int ret;
@ -3102,6 +3230,14 @@ void test_lp_local(int argn, char**argv) {
args_parser.print();
if (args_parser.option_is_used("-dji")) {
test_integer_domain();
return finalize(0);
}
if (args_parser.option_is_used("-cs")) {
test_cut_solver();
return finalize(0);
}
if (args_parser.option_is_used("--test_mpq")) {
test_rationals();
return finalize(0);
@ -3112,7 +3248,7 @@ void test_lp_local(int argn, char**argv) {
return finalize(0);
}
if (args_parser.option_is_used("--test_mpq_np_plus")) {
if (args_parser.option_is_used("--test_mpq_np_plus")) {
test_rationals_no_numeric_pairs_plus();
return finalize(0);
}

3
src/util/lp/CMakeLists.txt
View file

@ -3,7 +3,8 @@ z3_add_component(lp
lp_utils.cpp
binary_heap_priority_queue_instances.cpp
binary_heap_upair_queue_instances.cpp
lp_bound_propagator.cpp
bound_propagator.cpp
cut_solver.cpp
core_solver_pretty_printer_instances.cpp
dense_matrix_instances.cpp
eta_matrix_instances.cpp

4
src/util/lp/bound_propagator.cpp
View file

@ -30,7 +30,7 @@ void bound_propagator::try_add_bound(mpq v, unsigned j, bool is_low, bool coeff
return;
unsigned k; // index to ibounds
if (is_low) {
if (try_get_val(m_improved_low_bounds, j, k)) {
if (try_get_value(m_improved_low_bounds, j, k)) {
auto & found_bound = m_ibounds[k];
if (v > found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict)) {
found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
@ -42,7 +42,7 @@ void bound_propagator::try_add_bound(mpq v, unsigned j, bool is_low, bool coeff
TRACE("try_add_bound", m_lar_solver.print_implied_bound(m_ibounds.back(), tout););
}
} else { // the upper bound case
if (try_get_val(m_improved_upper_bounds, j, k)) {
if (try_get_value(m_improved_upper_bounds, j, k)) {
auto & found_bound = m_ibounds[k];
if (v < found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict)) {
found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);

8
src/util/lp/cut_solver.cpp Normal file
View file

@ -0,0 +1,8 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Nikolaj Bjorner, Lev Nachmanson
*/
#include "util/lp/cut_solver.h"
namespace lp {
}

359
src/util/lp/cut_solver.h Normal file
View file

@ -0,0 +1,359 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Nikolaj Bjorner, Lev Nachmanson
*/
#pragma once
#include "util/vector.h"
#include "util/trace.h"
#include "util/lp/lp_settings.h"
#include "util/lp/column_namer.h"
#include "util/lp/integer_domain.h"
#include <functional>
namespace lp {
template <typename T>
class cut_solver : public column_namer {
public: // for debugging
struct polynomial {
// the polynom evaluates to m_term + m_a
vector<std::pair<T, var_index>> m_term;
mpq m_a; // the free coefficient
polynomial(vector<std::pair<T, var_index>>& p, const mpq & a) : m_term(p), m_a(a) {
}
polynomial(vector<std::pair<T, var_index>>& p) : polynomial(p, 0) {
}
};
struct ineq { // we only have less or equal, which is enough for integral variables
polynomial m_poly;
ineq(vector<std::pair<T, var_index>>& term, const mpq& a): m_poly(term, a) {
}
};
vector<ineq> m_ineqs;
enum class lbool {
l_false, // false
l_true, // true
l_undef // undef
};
enum class literal_type {
BOOL,
INEQ,
BOUND
};
struct literal {
literal_type m_tag;
bool m_sign; // true means the pointed inequality is negated, or bound is negated, or boolean value is negated
bool m_is_lower;
unsigned m_id;
unsigned m_index_of_ineq; // index into m_ineqs
bool m_bool_val; // used if m_tag is equal to BOOL
mpq m_bound; // used if m_tag is BOUND
literal(bool sign, bool val): m_tag(literal_type::BOOL),
m_bool_val(val){
}
literal(bool sign, unsigned index_of_ineq) : m_tag(literal_type::INEQ), m_index_of_ineq(index_of_ineq) {}
};
struct var_info {
unsigned m_user_var_index;
var_info(unsigned user_var_index) : m_user_var_index(user_var_index) {}
vector<unsigned> m_literals; // point to m_trail
integer_domain<T> m_domain;
};
vector<var_info> m_var_infos;
bool lhs_is_int(const vector<std::pair<T, var_index>> & lhs) const {
for (auto & p : lhs) {
if (numeric_traits<T>::is_int(p.first) == false) return false;
}
return true;
}
public:
std::string get_column_name(unsigned j) const {
return m_var_name_function(m_var_infos[j].m_user_var_index);
}
unsigned add_ineq(vector<std::pair<T, var_index>> & lhs, const mpq& free_coeff) {
lp_assert(lhs_is_int(lhs));
lp_assert(free_coeff.is_int());
vector<std::pair<T, var_index>> local_lhs;
for (auto & p : lhs)
local_lhs.push_back(std::make_pair(p.first, add_var(p.second)));
m_ineqs.push_back(ineq(local_lhs, free_coeff));
return m_ineqs.size() - 1;
}
std::function<std::string (unsigned)> m_var_name_function;
bool m_inconsistent; // tracks if state is consistent
unsigned m_scope_lvl; // tracks the number of case splits
svector<literal> m_trail;
// backtracking state from the SAT solver:
struct scope {
unsigned m_trail_lim; // pointer into assignment stack
unsigned m_clauses_to_reinit_lim; // ignore for now
bool m_inconsistent; // really needed?
};
svector<scope> m_scopes;
std::unordered_map<unsigned, unsigned> m_user_vars_to_cut_solver_vars;
unsigned add_var(unsigned user_var_index) {
unsigned ret;
if (try_get_value(m_user_vars_to_cut_solver_vars, user_var_index, ret))
return ret;
unsigned j = m_var_infos.size();
m_var_infos.push_back(var_info(user_var_index));
return m_user_vars_to_cut_solver_vars[user_var_index] = j;
}
bool is_lower_bound(literal & l) const {
if (l.m_tag != literal_type::BOUND || l.m_is_lower)
return false;
l = l.m_bound;
return true;
}
bool lower_for_var(unsigned j, T & lower) const {
bool ret = false;
for (unsigned i : m_var_infos[j].m_literals)
if (is_lower_bound(m_trail[i])) {
if (ret == false) {
ret = true;
lower = get_bound(m_trail[i]);
} else {
lower = std::max(lower, get_bound(m_trail[i]));
}
}
return ret;
}
bool is_upper_bound(literal & l) const {
if (l.m_tag != literal_type::BOUND || l.m_is_upper)
return false;
l = l.m_bound;
return true;
}
bool upper_for_var(unsigned j, T & upper) const {
bool ret = false;
for (unsigned i : m_var_infos[j].m_literals)
if (is_upper_bound(m_trail[i])) {
if (ret == false) {
ret = true;
upper = get_bound(m_trail[i]);
} else {
upper = std::min(upper, get_bound(m_trail[i]));
}
}
return ret;
}
// used for testing only
void add_lower_bound_for_user_var(unsigned user_var_index, const T& bound) {
unsigned j = m_user_vars_to_cut_solver_vars[user_var_index];
auto & vi = m_var_infos[j];
vi.m_domain.intersect_with_lower_bound(bound);
}
// used for testing only
void add_upper_bound_for_user_var(unsigned user_var_index, const T& bound) {
unsigned j = m_user_vars_to_cut_solver_vars[user_var_index];
auto & vi = m_var_infos[j];
vi.m_domain.intersect_with_upper_bound(bound);
}
bool at_base_lvl() const { return m_scope_lvl == 0; }
lbool check() {
init_search();
propagate();
while (true) {
lbool r = bounded_search();
if (r != lbool::l_undef)
return r;
restart();
simplify_problem();
if (check_inconsistent()) return lbool::l_false;
gc();
}
}
cut_solver(std::function<std::string (unsigned)> var_name_function) : m_var_name_function(var_name_function) {
}
void init_search() {
// TBD
// initialize data-structures
}
void simplify_problem() {
// no-op
}
void gc() {
// no-op
}
void restart() {
// no-op for now
}
bool check_inconsistent() {
// TBD
return false;
}
lbool bounded_search() {
while (true) {
checkpoint();
bool done = false;
while (!done) {
lbool is_sat = propagate_and_backjump_step(done);
if (is_sat != lbool::l_true) return is_sat;
}
gc();
if (!decide()) {
lbool is_sat = final_check();
if (is_sat != lbool::l_undef) {
return is_sat;
}
}
}
}
void checkpoint() {
// check for cancelation
}
void cleanup() {
}
lbool propagate_and_backjump_step(bool& done) {
done = true;
propagate();
if (!inconsistent())
return lbool::l_true;
if (!resolve_conflict())
return lbool::l_false;
if (at_base_lvl()) {
cleanup(); // cleaner may propagate frozen clauses
if (inconsistent()) {
TRACE("sat", tout << "conflict at level 0\n";);
return lbool::l_false;
}
gc();
}
done = false;
return lbool::l_true;
}
lbool final_check() {
// there are no more case splits, and all clauses are satisfied.
// prepare the model for external consumption.
return lbool::l_true;
}
bool resolve_conflict() {
while (true) {
bool r = resolve_conflict_core();
// after pop, clauses are reinitialized,
// this may trigger another conflict.
if (!r)
return false;
if (!inconsistent())
return true;
}
}
bool resolve_conflict_core() {
// this is where the main action is.
return true;
}
void propagate() {
// this is where the main action is.
}
bool decide() {
// this is where the main action is.
// pick the next variable and bound or value on the variable.
// return false if all variables have been assigned.
return false;
}
bool inconsistent() const { return m_inconsistent; }
bool get_var_low_bound(var_index i, T & bound) const {
const var_info & v = m_var_infos[i];
return v.m_domain.get_lower_bound(bound);
}
bool get_var_upper_bound(var_index i, T & bound) const {
const var_info & v = m_var_infos[i];
return v.m_domain.get_upper_bound(bound);
}
// finds the lower bound of the monomial,
// otherwise returns false
bool lower_monomial(const std::pair<T, var_index> & p, mpq & lb) const {
lp_assert(p.first != 0);
T var_bound;
if (p.first > 0) {
if (!get_var_low_bound(p.second, var_bound))
return false;
lb = p.first * var_bound;
}
else {
if (!get_var_upper_bound(p.second, var_bound))
return false;
lb = p.first * var_bound;
}
return true;
}
// returns false if not limited from below
// otherwise the answer is put into lp
bool lower(const polynomial & f, mpq & lb) const {
lb = 0;
mpq lm;
for (const auto & p : f.m_term) {
if (lower_monomial(p, lm)) {
lb += lm;
} else {
return false;
}
}
lb += f.m_a;
return true;
}
void print_ineq(unsigned i, std::ostream & out) {
print_polynomial(m_ineqs[i].m_poly, out);
out << " <= 0" << std::endl;
}
void print_polynomial(const polynomial & p, std::ostream & out) {
this->print_linear_combination_of_column_indices(p.m_term, out);
if (!is_zero(p.m_a)) {
if (p.m_a < 0) {
out << " - " << -p.m_a;
} else {
out << " + " << p.m_a;
}
}
}
};
}

334
src/util/lp/disjoint_intervals.h
View file

@ -1,334 +0,0 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include <map>
namespace lp {
// represents the set of disjoint intervals of integer number
struct disjoint_intervals {
std::map<int, short> m_endpoints; // 0 means start, 1 means end, 2 means both - for a point interval
bool m_empty;
// constructors create an interval containing all integer numbers or an empty interval
disjoint_intervals() : m_empty(false) {}
disjoint_intervals(bool is_empty) : m_empty(is_empty) {}
bool is_start(short x) const { return x == 0 || x == 2; }
bool is_start(const std::map<int, short>::iterator & it) const {
return is_start(it->second);
}
bool is_start(const std::map<int, short>::reverse_iterator & it) const {
return is_start(it->second);
}
bool is_end(short x) const { return x == 1 || x == 2; }
bool is_end(const std::map<int, short>::iterator & it) const {
return is_end(it->second);
}
bool is_end(const std::map<int, short>::reverse_iterator & it) const {
return is_end(it->second);
}
int pos(const std::map<int, short>::iterator & it) const {
return it->first;
}
int pos(const std::map<int, short>::reverse_iterator & it) const {
return it->first;
}
int bound_kind(const std::map<int, short>::iterator & it) const {
return it->second;
}
int bound_kind(const std::map<int, short>::reverse_iterator & it) const {
return it->second;
}
bool is_proper_start(short x) const { return x == 0; }
bool is_proper_end(short x) const { return x == 1; }
bool is_proper_end(const std::map<int, short>::iterator & it) const {
return is_proper_end(it->second);
}
bool is_proper_end(const std::map<int, short>::reverse_iterator & it) const {
return is_proper_end(it->second);
}
bool is_one_point_interval(short x) const { return x == 2; }
bool is_one_point_interval(const std::map<int, short>::iterator & it) const {
return is_one_point_interval(it->second);
}
bool is_one_point_interval(const std::map<int, short>::reverse_iterator & it) const {
return is_one_point_interval(it->second);
}
void erase(int x) {
m_endpoints.erase(x);
}
void set_one_point_segment(int x) {
m_endpoints[x] = 2;
}
void set_start(int x) {
m_endpoints[x] = 0;
}
void set_end(int x) {
m_endpoints[x] = 1;
}
void remove_all_endpoints_below(int x) {
while (m_endpoints.begin() != m_endpoints.end() && m_endpoints.begin()->first < x)
m_endpoints.erase(m_endpoints.begin());
}
// we intersect the existing set with the half open to the right interval
void intersect_with_lower_bound(int x) {
if (m_empty)
return;
if (m_endpoints.empty()) {
set_start(x);
return;
}
bool pos_inf = has_pos_inf();
auto it = m_endpoints.begin();
while (it != m_endpoints.end() && pos(it) < x) {
m_endpoints.erase(it);
it = m_endpoints.begin();
}
if (m_endpoints.empty()) {
if (!pos_inf) {
m_empty = true;
return;
}
set_start(x);
return;
}
lp_assert(pos(it) >= x);
if (pos(it) == x) {
if (is_proper_end(it))
set_one_point_segment(x);
}
else { // x(it) > x
if (is_proper_end(it)) {
set_start(x);
}
}
lp_assert(is_correct());
}
// we intersect the existing set with the half open interval
void intersect_with_upper_bound(int x) {
if (m_empty)
return;
if (m_endpoints.empty()) {
set_end(x);
return;
}
bool neg_inf = has_neg_inf();
auto it = m_endpoints.rbegin();
while (!m_endpoints.empty() && pos(it) > x) {
m_endpoints.erase(std::prev(m_endpoints.end()));
it = m_endpoints.rbegin();
}
if (m_endpoints.empty()) {
if (!neg_inf) {
m_empty = true;
return;
}
set_end(x);
}
lp_assert(pos(it) <= x);
if (pos(it) == x) {
if (is_one_point_interval(it)) {}
else if (is_proper_end(it)) {}
else {// is_proper_start(it->second)
set_one_point_segment(x);
}
}
else { // pos(it) < x}
if (is_start(it))
set_end(x);
}
lp_assert(is_correct());
}
bool has_pos_inf() const {
if (m_empty)
return false;
if (m_endpoints.empty())
return true;
lp_assert(m_endpoints.rbegin() != m_endpoints.rend());
return m_endpoints.rbegin()->second == 0;
}
bool has_neg_inf() const {
if (m_empty)
return false;
if (m_endpoints.empty())
return true;
auto it = m_endpoints.begin();
return is_proper_end(it->second);//m_endpoints.begin());
}
// we are intersecting
void intersect_with_interval(int x, int y) {
if (m_empty)
return;
lp_assert(x <= y);
intersect_with_lower_bound(x);
intersect_with_upper_bound(y);
}
// add an intervar [x, inf]
void unite_with_interval_x_pos_inf(int x) {
if (m_empty) {
set_start(x);
m_empty = false;
return;
}
while (!m_endpoints.empty() && pos(m_endpoints.rbegin()) > x) {
m_endpoints.erase(std::prev(m_endpoints.end()));
}
if (m_endpoints.empty()) {
set_start(x);
return;
}
auto it = m_endpoints.rbegin();
lp_assert(pos(it) <= x);
if (pos(it) == x) {
if (is_end(it)) {
m_endpoints.erase(x);
} else {
set_start(x);
}
} else if (pos(it) == x - 1 && is_end(it)) {
m_endpoints.erase(x - 1); // closing the gap
} else {
if (!has_pos_inf())
set_start(x);
}
}
// add an interval [-inf, x]
void unite_with_interval_neg_inf_x(int x) {
if (m_empty) {
set_end(x);
m_empty = false;
return;
}
auto it = m_endpoints.upper_bound(x);
if (it == m_endpoints.end()) {
bool pos_inf = has_pos_inf();
m_endpoints.clear();
// it could be the case where x is inside of the last infinite interval with pos inf
if (!pos_inf)
set_end(x);
return;
}
lp_assert(pos(it) > x);
if (is_one_point_interval(pos(it))) {
set_end(it->second);
} else {
if (is_start(it->second)) {
set_end(x);
}
}
while (!m_endpoints.empty() && m_endpoints.begin()->first < x) {
m_endpoints.erase(m_endpoints.begin());
}
lp_assert(is_correct());
}
void unite_with_interval(int x, int y) {
lp_assert(false); // not implemented
}
bool is_correct() const {
if (m_empty) {
if (m_endpoints.size() > 0) {
std::cout << "is empty is true but m_endpoints.size() = " << m_endpoints.size() << std::endl;
return false;
}
return true;
}
bool expect_end;
bool prev = false;
int prev_x;
for (auto t : m_endpoints) {
if (prev && (expect_end != t.second > 0)) {
std::cout << "x = " << t.first << "\n";
if (expect_end) {
std::cout << "expecting an interval end\n";
} else {
std::cout << "expecting an interval start\n";
}
return false;
}
if (t.second == 2) {
expect_end = false; // swallow a point interval
} else {
if (prev)
expect_end = !expect_end;
else
expect_end = is_start(t.second);
}
if (prev) {
if (t.first - prev_x <= 1) {
std::cout << "the sequence is not increasing or the gap is too small: " << prev_x << ", " << t.first << std::endl;
return false;
}
}
prev = true;
prev_x = t.first;
}
return true;
}
void print(std::ostream & out) const {
if (m_empty) {
out << "empty\n";
return;
}
if (m_endpoints.empty()){
out << "[-oo,oo]\n";
return;
}
bool first = true;
for (auto t : m_endpoints) {
if (first) {
if (t.second == 1) {
out << "[-oo," << t.first << "]";
}
else if (t.second == 0)
out << "[" << t.first << ",";
else if (t.second == 2)
out << "[" << t.first << "]";
first = false;
} else {
if (t.second==0)
out << "[" << t.first << ",";
else if (t.second == 1)
out << t.first << "]";
else if (t.second == 2)
out << "[" << t.first << "]";
}
}
if (has_pos_inf())
out << "oo]";
out << "\n";
}
};
}

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

@ -5,24 +5,10 @@
#include "util/lp/int_solver.h"
#include "util/lp/lar_solver.h"
#include "util/lp/cut_solver.h"
#include <utility>
namespace lp {
void int_solver::fix_non_base_columns() {
auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
bool change = false;
for (unsigned j : lcs.m_r_nbasis) {
if (column_is_int_inf(j)) {
change = true;
set_value_for_nbasic_column(j, floor(lcs.m_r_x[j].x));
}
}
if (!change)
return;
if (m_lar_solver->find_feasible_solution() == lp_status::INFEASIBLE)
failed();
lp_assert(is_feasible() && inf_int_set_is_correct());
}
void int_solver::failed() {
auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
@ -278,10 +264,10 @@ bool int_solver::current_solution_is_inf_on_cut(const lar_term& t, const mpq& k)
const auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x;
impq v = t.apply(x);
TRACE(
"current_solution_is_inf_on_cut", tout << "v = " << v << " k = " << k << std::endl;
if (v <=k) {
tout << "v <= k - it should not happen!\n";
}
"current_solution_is_inf_on_cut", tout << "v = " << v << " k = " << k << std::endl;
if (v <=k) {
tout << "v <= k - it should not happen!\n";
}
);
return v > k;
@ -379,9 +365,9 @@ lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsi
}
void int_solver::init_check_data() {
unsigned n = m_lar_solver->A_r().column_count();
m_old_values_set.resize(n);
m_old_values_data.resize(n);
unsigned n = m_lar_solver->A_r().column_count();
m_old_values_set.resize(n);
m_old_values_data.resize(n);
}
int int_solver::find_free_var_in_gomory_row(linear_combination_iterator<mpq>& iter) {
@ -394,117 +380,166 @@ int int_solver::find_free_var_in_gomory_row(linear_combination_iterator<mpq>& it
return -1;
}
lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, linear_combination_iterator<mpq>& iter) {
int free_j = find_free_var_in_gomory_row(iter);
lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j) {
lia_move ret;
linear_combination_iterator<mpq>* iter = m_lar_solver->get_iterator_on_row(row_of_basic_column(j));
int free_j = find_free_var_in_gomory_row(*iter);
if (free_j != -1) {
lp_assert(t.is_empty());
t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(free_j));
k = zero_of_type<mpq>();
return lia_move::branch; // branch on a free column
ret = create_branch_on_column(j, t, k, true);
} else if (!is_gomory_cut_target(*iter)) {
ret = create_branch_on_column(j, t, k, false);
} else {
ret = mk_gomory_cut(t, k, ex, j, *iter);
}
if (!is_gomory_cut_target(iter))
return create_branch_on_column(j, t, k);
return mk_gomory_cut(t, k, ex, j, iter);
}
delete iter;
return ret;
}
unsigned int_solver::row_of_basic_column(unsigned j) const {
return m_lar_solver->m_mpq_lar_core_solver.m_r_heading[j];
}
template <typename T>
void int_solver::fill_cut_solver(cut_solver<T> & cs) {
for (lar_base_constraint * c : m_lar_solver->constraints())
fill_cut_solver_for_constraint(c, cs);
}
template <typename T>
void int_solver::fill_cut_solver_for_constraint(const lar_base_constraint* c, cut_solver<T> & cs) {
vector<std::pair<T, var_index>> coeffs;
T rs;
get_int_coeffs_from_constraint(c, coeffs, rs);
cs.add_ineq(coeffs, rs);
}
// it produces an inequality coeff*x <= rs
template <typename T>
void int_solver::get_int_coeffs_from_constraint(const lar_base_constraint* c, vector<std::pair<T, var_index>>& coeffs, T & rs) {
lp_assert(c->m_kind != EQ); // it is not implemented, we need to create two inequalities in this case
int sign = ((int)c->m_kind > 0) ? -1 : 1;
vector<std::pair<T, var_index>> lhs = c->get_left_side_coefficients();
T den = denominator(c->m_right_side);
for (auto & kv : lhs) {
den = lcm(den, denominator(kv.first));
}
lp_assert(den > 0);
for (auto& kv : lhs) {
coeffs.push_back(std::make_pair(den * kv.first * sign, kv.second));
}
rs = den * c->m_right_side * sign;
if (kind_is_strict(c->m_kind))
rs--;
}
// this will allow to enable and disable tracking of the pivot rows
struct pivoted_rows_tracking_control {
lar_solver * m_lar_solver;
bool m_track_pivoted_rows;
pivoted_rows_tracking_control(lar_solver* ls) :
m_lar_solver(ls),
m_track_pivoted_rows(ls->get_track_pivoted_rows())
{
TRACE("pivoted_rows", tout << "pivoted rows = " << ls->m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows->size() << std::endl;);
m_lar_solver->set_track_pivoted_rows(false);
}
~pivoted_rows_tracking_control() {
TRACE("pivoted_rows", tout << "pivoted rows = " << m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows->size() << std::endl;);
m_lar_solver->set_track_pivoted_rows(m_track_pivoted_rows);
}
};
lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
init_check_data();
init_check_data();
lp_assert(inf_int_set_is_correct());
// it is mostly a reimplementation of
// final_check_status theory_arith<Ext>::check_int_feasibility()
// from theory_arith_int.h
if (!has_inf_int())
return lia_move::ok;
if (!has_inf_int())
return lia_move::ok;
if (settings().m_run_gcd_test)
if (!gcd_test(ex))
return lia_move::conflict;
/*
if (m_params.m_arith_euclidean_solver)
apply_euclidean_solver();
*/
bool track_pivoted_rows = m_lar_solver->get_track_pivoted_rows();
m_lar_solver->set_track_pivoted_rows(false);
m_lar_solver->pivot_fixed_vars_from_basis();
lean_assert(m_lar_solver->m_mpq_lar_core_solver.r_basis_is_OK());
patch_int_infeasible_columns();
lean_assert(m_lar_solver->m_mpq_lar_core_solver.r_basis_is_OK());
fix_non_base_columns();
lean_assert(m_lar_solver->m_mpq_lar_core_solver.r_basis_is_OK());
lean_assert(is_feasible());
TRACE("arith_int_rows", trace_inf_rows(););
if (!has_inf_int()) {
m_lar_solver->set_track_pivoted_rows(track_pivoted_rows);
pivoted_rows_tracking_control pc(m_lar_solver);
/* if (m_params.m_arith_euclidean_solver) apply_euclidean_solver(); */
//m_lar_solver->pivot_fixed_vars_from_basis();
patch_int_infeasible_nbasic_columns();
if (!has_inf_int())
return lia_move::ok;
}
TRACE("gomory_cut", tout << m_branch_cut_counter+1 << ", " << settings().m_int_branch_cut_gomory_threshold << std::endl;);
if ((++m_branch_cut_counter) % settings().m_int_branch_cut_gomory_threshold == 0) {
if (move_non_base_vars_to_bounds()) {
lp_status st = m_lar_solver->find_feasible_solution();
lp_assert(non_basic_columns_are_at_bounds());
if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
TRACE("arith_int", tout << "give_up\n";);
m_lar_solver->set_track_pivoted_rows(track_pivoted_rows);
return lia_move::give_up;
}
}
int j = find_inf_int_base_column();
lp_assert(j != -1);
TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";);
auto iter_on_gomory_row = m_lar_solver->get_iterator_on_row(row_of_basic_column(j));
lia_move ret = proceed_with_gomory_cut(t, k, ex, j, *iter_on_gomory_row);
delete iter_on_gomory_row;
m_lar_solver->set_track_pivoted_rows(track_pivoted_rows);
return ret;
}
m_lar_solver->set_track_pivoted_rows(track_pivoted_rows);
return create_branch_on_column(find_inf_int_base_column(), t, k);
// lp_assert(non_basic_columns_are_at_bounds());
TRACE("gomory_cut", tout << m_branch_cut_counter+1 << ", " << settings().m_int_branch_cut_gomory_threshold << std::endl;);
if (++m_branch_cut_counter > 0) { // testing cut_solver
cut_solver<mpq> cs([this](unsigned j) {return m_lar_solver->get_column_name(j);});
fill_cut_solver(cs);
} else
if ((++m_branch_cut_counter) % settings().m_int_branch_cut_gomory_threshold == 0) {
if (move_non_basic_columns_to_bounds()) {
lp_status st = m_lar_solver->find_feasible_solution();
lp_assert(non_basic_columns_are_at_bounds());
if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
TRACE("arith_int", tout << "give_up\n";);
return lia_move::give_up;
}
}
int j = find_inf_int_base_column();
if (j == -1) return lia_move::ok;
TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";);
return proceed_with_gomory_cut(t, k, ex, j);
}
return create_branch_on_column(find_inf_int_base_column(), t, k, false);
}
bool int_solver::move_non_base_vars_to_bounds() {
bool int_solver::move_non_basic_column_to_bounds(unsigned j) {
auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
auto & val = lcs.m_r_x[j];
switch (lcs.m_column_types()[j]) {
case column_type::boxed:
if (val != lcs.m_r_low_bounds()[j] && val != lcs.m_r_upper_bounds()[j]) {
if (random() % 2 == 0)
set_value_for_nbasic_column(j, lcs.m_r_low_bounds()[j]);
else
set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
return true;
}
break;
case column_type::low_bound:
if (val != lcs.m_r_low_bounds()[j]) {
set_value_for_nbasic_column(j, lcs.m_r_low_bounds()[j]);
return true;
}
break;
case column_type::upper_bound:
if (val != lcs.m_r_upper_bounds()[j]) {
set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
return true;
}
break;
default:
if (is_int(j) && !val.is_int()) {
set_value_for_nbasic_column(j, impq(floor(val)));
return true;
}
break;
}
return false;
}
bool int_solver::move_non_basic_columns_to_bounds() {
auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
bool change = false;
for (unsigned j : lcs.m_r_nbasis) {
auto & val = lcs.m_r_x[j];
switch (lcs.m_column_types()[j]) {
case column_type::boxed:
if (val != lcs.m_r_low_bounds()[j] && val != lcs.m_r_upper_bounds()[j]) {
set_value_for_nbasic_column(j, lcs.m_r_low_bounds()[j]);
change = true;
}
break;
case column_type::low_bound:
if (val != lcs.m_r_low_bounds()[j]) {
set_value_for_nbasic_column(j, lcs.m_r_low_bounds()[j]);
change = true;
}
break;
case column_type::upper_bound:
if (val != lcs.m_r_upper_bounds()[j]) {
set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
change = true;
}
break;
default:
if (is_int(j) && !val.is_int()) {
set_value_for_nbasic_column(j, impq(floor(val)));
change = true;
}
}
if (move_non_basic_column_to_bounds(j))
change = true;
}
return change;
}
void int_solver::set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val) {
lp_assert(!is_base(j));
lp_assert(!is_base(j));
auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x[j];
auto delta = new_val - x;
x = new_val;
@ -514,7 +549,7 @@ void int_solver::set_value_for_nbasic_column_ignore_old_values(unsigned j, const
void int_solver::set_value_for_nbasic_column(unsigned j, const impq & new_val) {
lp_assert(!is_base(j));
lp_assert(!is_base(j));
auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x[j];
if (m_lar_solver->has_int_var() && !m_old_values_set.contains(j)) {
m_old_values_set.insert(j);
@ -526,53 +561,68 @@ void int_solver::set_value_for_nbasic_column(unsigned j, const impq & new_val) {
m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta);
}
void int_solver::patch_int_infeasible_columns() {
void int_solver::patch_int_infeasible_non_basic_column(unsigned j) {
if (!is_int(j)) return;
bool inf_l, inf_u;
impq l, u;
mpq m;
auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
for (unsigned j : lcs.m_r_nbasis) {
if (!is_int(j))
continue;
get_freedom_interval_for_column(j, inf_l, l, inf_u, u, m);
impq & val = lcs.m_r_x[j];
bool val_is_int = val.is_int();
bool m_is_one = m.is_one();
if (m.is_one() && val_is_int)
continue;
// check whether value of j is already a multiple of m.
if (val_is_int && (val.x / m).is_int())
continue;
TRACE("patch_int",
tout << "TARGET j" << j << " -> [";
if (inf_l) tout << "-oo"; else tout << l;
tout << ", ";
if (inf_u) tout << "oo"; else tout << u;
tout << "]";
tout << ", m: " << m << ", val: " << val << ", is_int: " << m_lar_solver->column_is_int(j) << "\n";);
if (!inf_l) {
l = m_is_one? ceil(l) : m * ceil(l / m);
if (inf_u || l <= u) {
TRACE("patch_int",
tout << "patching with l: " << l << '\n';);
set_value(j, l);
} else {
TRACE("patch_int",
tout << "not patching " << l << "\n";);
}
} else if (!inf_u) {
u = m_is_one? floor(u) : m * floor(u / m);
set_value(j, u);
TRACE("patch_int",
tout << "patching with u: " << u << '\n';);
} else {
set_value(j, impq(0));
TRACE("patch_int",
tout << "patching with 0\n";);
}
lp_assert(is_feasible() && inf_int_set_is_correct());
if (!get_value(j).is_int() || !get_freedom_interval_for_column(j, inf_l, l, inf_u, u, m)) {
move_non_basic_column_to_bounds(j);
return;
}
auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
impq & val = lcs.m_r_x[j];
bool val_is_int = val.is_int();
bool m_is_one = m.is_one();
if (m.is_one() && val_is_int)
return;
// check whether value of j is already a multiple of m.
if (val_is_int && (val.x / m).is_int())
return;
TRACE("patch_int",
tout << "TARGET j" << j << " -> [";
if (inf_l) tout << "-oo"; else tout << l;
tout << ", ";
if (inf_u) tout << "oo"; else tout << u;
tout << "]";
tout << ", m: " << m << ", val: " << val << ", is_int: " << m_lar_solver->column_is_int(j) << "\n";);
if (!inf_l) {
l = m_is_one ? ceil(l) : m * ceil(l / m);
if (inf_u || l <= u) {
TRACE("patch_int",
tout << "patching with l: " << l << '\n';);
set_value_for_nbasic_column(j, l);
}
else {
TRACE("patch_int",
tout << "not patching " << l << "\n";);
}
}
else if (!inf_u) {
u = m_is_one ? floor(u) : m * floor(u / m);
set_value_for_nbasic_column(j, u);
TRACE("patch_int",
tout << "patching with u: " << u << '\n';);
}
else {
set_value_for_nbasic_column(j, impq(0));
TRACE("patch_int",
tout << "patching with 0\n";);
}
}
void int_solver::patch_int_infeasible_nbasic_columns() {
lp_assert(is_feasible());
for (unsigned j : m_lar_solver->m_mpq_lar_core_solver.m_r_nbasis) {
patch_int_infeasible_non_basic_column(j);
if (!is_feasible())
break;
}
if (!is_feasible()) {
move_non_basic_columns_to_bounds();
m_lar_solver->find_feasible_solution();
}
lp_assert(is_feasible() && inf_int_set_is_correct());
}
mpq get_denominators_lcm(iterator_on_row<mpq> &it) {
@ -595,55 +645,55 @@ bool int_solver::gcd_test_for_row(static_matrix<mpq, numeric_pair<mpq>> & A, uns
mpq a;
unsigned j;
while (it.next(a, j)) {
if (m_lar_solver->column_is_fixed(j)) {
mpq aux = lcm_den * a;
consts += aux * m_lar_solver->column_low_bound(j).x;
}
else if (m_lar_solver->column_is_real(j)) {
return true;
}
else if (gcds.is_zero()) {
gcds = abs(lcm_den * a);
least_coeff = gcds;
if (m_lar_solver->column_is_fixed(j)) {
mpq aux = lcm_den * a;
consts += aux * m_lar_solver->column_low_bound(j).x;
}
else if (m_lar_solver->column_is_real(j)) {
return true;
}
else if (gcds.is_zero()) {
gcds = abs(lcm_den * a);
least_coeff = gcds;
least_coeff_is_bounded = m_lar_solver->column_is_bounded(j);
}
else {
mpq aux = abs(lcm_den * a);
gcds = gcd(gcds, aux);
if (aux < least_coeff) {
least_coeff = aux;
least_coeff_is_bounded = m_lar_solver->column_is_bounded(j);
}
else {
mpq aux = abs(lcm_den * a);
gcds = gcd(gcds, aux);
if (aux < least_coeff) {
least_coeff = aux;
least_coeff_is_bounded = m_lar_solver->column_is_bounded(j);
}
else if (least_coeff_is_bounded && aux == least_coeff) {
least_coeff_is_bounded = m_lar_solver->column_is_bounded(j);
}
else if (least_coeff_is_bounded && aux == least_coeff) {
least_coeff_is_bounded = m_lar_solver->column_is_bounded(j);
}
SASSERT(gcds.is_int());
SASSERT(least_coeff.is_int());
TRACE("gcd_test_bug", tout << "coeff: " << a << ", gcds: " << gcds
<< " least_coeff: " << least_coeff << " consts: " << consts << "\n";);
}
SASSERT(gcds.is_int());
SASSERT(least_coeff.is_int());
TRACE("gcd_test_bug", tout << "coeff: " << a << ", gcds: " << gcds
<< " least_coeff: " << least_coeff << " consts: " << consts << "\n";);
}
if (gcds.is_zero()) {
// All variables are fixed.
// This theory guarantees that the assignment satisfies each row, and
// fixed integer variables are assigned to integer values.
return true;
}
if (!(consts / gcds).is_int())
fill_explanation_from_fixed_columns(it, ex);
if (least_coeff.is_one() && !least_coeff_is_bounded) {
SASSERT(gcds.is_one());
return true;
}
if (least_coeff_is_bounded) {
return ext_gcd_test(it, least_coeff, lcm_den, consts, ex);
}
if (gcds.is_zero()) {
// All variables are fixed.
// This theory guarantees that the assignment satisfies each row, and
// fixed integer variables are assigned to integer values.
return true;
}
if (!(consts / gcds).is_int())
fill_explanation_from_fixed_columns(it, ex);
if (least_coeff.is_one() && !least_coeff_is_bounded) {
SASSERT(gcds.is_one());
return true;
}
if (least_coeff_is_bounded) {
return ext_gcd_test(it, least_coeff, lcm_den, consts, ex);
}
return true;
}
void int_solver::add_to_explanation_from_fixed_or_boxed_column(unsigned j, explanation & ex) {
@ -663,14 +713,14 @@ void int_solver::fill_explanation_from_fixed_columns(iterator_on_row<mpq> & it,
}
bool int_solver::gcd_test(explanation & ex) {
auto & A = m_lar_solver->A_r(); // getting the matrix
for (unsigned i = 0; i < A.row_count(); i++)
auto & A = m_lar_solver->A_r(); // getting the matrix
for (unsigned i = 0; i < A.row_count(); i++)
if (!gcd_test_for_row(A, i, ex)) {
std::cout << "false from gcd_test\n" ;
return false;
}
return true;
return true;
}
bool int_solver::ext_gcd_test(iterator_on_row<mpq> & it,
@ -735,7 +785,7 @@ bool int_solver::ext_gcd_test(iterator_on_row<mpq> & it,
linear_combination_iterator<mpq> * int_solver::get_column_iterator(unsigned j) {
if (m_lar_solver->use_tableau())
return new iterator_on_column<mpq, impq>(m_lar_solver->A_r().m_columns[j], m_lar_solver->A_r());
return new iterator_on_indexed_vector<mpq>(m_lar_solver->get_column_in_lu_mode(j));
return new iterator_on_indexed_vector<mpq>(m_lar_solver->get_column_in_lu_mode(j));
}
@ -866,8 +916,8 @@ bool int_solver::value_is_int(unsigned j) const {
bool int_solver::is_feasible() const {
auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
lp_assert(
lcs.m_r_solver.calc_current_x_is_feasible_include_non_basis() ==
lcs.m_r_solver.current_x_is_feasible());
lcs.m_r_solver.calc_current_x_is_feasible_include_non_basis() ==
lcs.m_r_solver.current_x_is_feasible());
return lcs.m_r_solver.current_x_is_feasible();
}
const impq & int_solver::get_value(unsigned j) const {
@ -1073,17 +1123,15 @@ const impq& int_solver::low_bound(unsigned j) const {
return m_lar_solver->column_low_bound(j);
}
lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k) const {
lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column) const {
lp_assert(t.is_empty());
lp_assert(j != -1);
t.add_monomial(mpq(1), j);
k = floor(get_value(j));
t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j));
k = free_column? mpq(0) : floor(get_value(j));
TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n";
display_column(tout, j);
tout << "k = " << k << std::endl;
);
lp_assert(current_solution_is_inf_on_cut(t, k));
m_lar_solver->subs_term_columns(t);
return lia_move::branch;
}
@ -1091,4 +1139,14 @@ lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k) const {
const impq& int_solver::upper_bound(unsigned j) const {
return m_lar_solver->column_upper_bound(j);
}
void int_solver::display_inf_or_int_inf_columns(std::ostream & out) const {
out << "int inf\n";
for (unsigned j : m_lar_solver->m_inf_int_set.m_index) {
display_column(out, j);
}
out << "regular inf\n";
for (unsigned j : m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_inf_set.m_index) {
display_column(out, j);
}
}
}

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

@ -8,6 +8,9 @@
#include "util/lp/iterator_on_row.h"
#include "util/lp/int_set.h"
#include "util/lp/lar_term.h"
#include "util/lp/cut_solver.h"
#include "util/lp/lar_constraints.h"
namespace lp {
class lar_solver;
template <typename T, typename X>
@ -42,6 +45,8 @@ public:
// main function to check that solution provided by lar_solver is valid for integral values,
// or provide a way of how it can be adjusted.
lia_move check(lar_term& t, mpq& k, explanation& ex);
bool move_non_basic_column_to_bounds(unsigned j);
lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex);
private:
// how to tighten bounds for integer variables.
@ -73,8 +78,8 @@ private:
explanation & ex);
void fill_explanation_from_fixed_columns(iterator_on_row<mpq> & it, explanation &);
void add_to_explanation_from_fixed_or_boxed_column(unsigned j, explanation &);
void remove_fixed_vars_from_base();
void patch_int_infeasible_columns();
void patch_int_infeasible_non_basic_column(unsigned j);
void patch_int_infeasible_nbasic_columns();
bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m);
linear_combination_iterator<mpq> * get_column_iterator(unsigned j);
const impq & low_bound(unsigned j) const;
@ -88,7 +93,7 @@ private:
bool value_is_int(unsigned j) const;
void set_value_for_nbasic_column(unsigned j, const impq & new_val);
void set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val);
void fix_non_base_columns();
bool non_basic_columns_are_at_bounds() const;
void failed();
bool is_feasible() const;
const impq & get_value(unsigned j) const;
@ -100,14 +105,14 @@ private:
int find_inf_int_base_column();
int find_inf_int_boxed_base_column_with_smallest_range();
lp_settings& settings();
bool move_non_base_vars_to_bounds();
bool move_non_basic_columns_to_bounds();
void branch_infeasible_int_var(unsigned);
lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex, unsigned inf_col, linear_combination_iterator<mpq>& iter);
lia_move report_conflict_from_gomory_cut(mpq & k);
void adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq& lcm_den);
void init_check_data();
bool constrain_free_vars(linear_combination_iterator<mpq> * r);
lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, linear_combination_iterator<mpq>& iter);
lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j);
int find_free_var_in_gomory_row(linear_combination_iterator<mpq>& iter);
bool is_gomory_cut_target(linear_combination_iterator<mpq> &iter);
bool at_bound(unsigned j) const;
@ -136,8 +141,16 @@ public:
bool shift_var(unsigned j, unsigned range);
private:
unsigned random();
bool non_basic_columns_are_at_bounds() const;
bool has_inf_int() const;
lia_move create_branch_on_column(int j, lar_term& t, mpq& k) const;
lia_move create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column) const;
public:
void display_inf_or_int_inf_columns(std::ostream & out) const;
template <typename T>
void fill_cut_solver(cut_solver<T> & cs);
template <typename T>
void fill_cut_solver_for_constraint(const lar_base_constraint*, cut_solver<T>& );
template <typename T>
void get_int_coeffs_from_constraint(const lar_base_constraint* c, vector<std::pair<T, var_index>>& coeff, T & rs);
};
}

876
src/util/lp/integer_domain.h Normal file
View file

@ -0,0 +1,876 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#pragma once
#include <map>
#include "util/trace.h"
#include "util/lp/stacked_value.h"
#include "util/lp/stacked_map.h"
namespace lp {
// represents the set of disjoint intervals of integer number
template <typename T>
class integer_domain {
#ifdef Z3DEBUG
std::set<int> m_domain;
#endif
typedef typename std::map<T, char>::iterator iter;
typedef typename std::map<T, char>::const_iterator const_iter;
typedef typename std::map<T, char>::reverse_iterator riter;
stacked_map<T, char> m_endpoints; // 0 means start, 1 means end, 2 means both - for a point interval
stacked_value<bool> m_empty;
public:
// the default constructor creates a set containing all integer numbers
integer_domain() : integer_domain(false) {}
// if is_empty = false then the constructor creates a set containing all integer numbers,
// otherwise it creates an empty set
integer_domain(bool is_empty) : m_empty(is_empty) {
#if Z3DEBUG
if (!is_empty) {
for (int i = 0; i <= 100; i++)
m_domain.insert(i);
}
#endif
}
void init_to_contain_all() {
m_empty = false;
m_endpoints.clear();
#if Z3DEBUG
for (int i = 0; i <= 100; i++)
m_domain.insert(i);
#endif
}
// copy constructor
integer_domain(const integer_domain<T> & t) :
#if Z3DEBUG
m_domain(t.m_domain),
#endif
m_endpoints(t.m_endpoints),
m_empty(t.m_empty)
{
}
// needed for debug only
void restore_domain() {
#if Z3DEBUG
for (int i = 0; i <= 100; i++)
if (contains(i))
m_domain.insert(i);
else
m_domain.erase(i);
#endif
}
bool operator==(const integer_domain<T> & t) {
return m_empty == t.m_empty && m_endpoints() == t.m_endpoints();
}
bool contains_all() const {
return !m_empty && m_endpoints.empty();
}
bool contains(const T & x) const {
if (contains_all())
return true;
bool neg_inf;
const_iter l;
bool found_left_point = get_left_point(x, neg_inf, l);
if (!found_left_point)
return has_neg_inf();
if (neg_inf)
return true;
if (pos(l) == x)
return true;
return is_proper_start(l);
}
void handle_right_point_in_union(const_iter &r, const T &y) {
if (pos(r) == y) {
if (is_proper_start(r))
erase(r);
else
set_end(y);
}
else if (pos(r) == y + 1) {
if (is_proper_start(r)) {
erase(r);
}
else {
set_end(r->first);
}
}
else if (!is_proper_end(r))
set_end(y);
lp_assert(is_correct());
}
void handle_left_point_in_union(const_iter& l, const T &x, const T & y) {
if (pos(l) == x || pos(l) + 1 == x) {
if (is_proper_end(l)) {
l++;
erase(std::prev(l));
}
else {
if (is_one_point_interval(l)) {
set_start(pos(l));
}
}
}
else {
if (!is_proper_start(l)) {
set_start(x);
}
}
while (l!= m_endpoints.end() && pos(l) <= x)
l++;
remove_from_the_left(y, l);
}
void unite_with_interval(const T& x, const T& y) {
TRACE("disj_intervals", tout << "unite_with_interval(" << x << ", " << y << ")\n";);
#if Z3DEBUG
for (int i = std::max(x, 0); i <= std::min(100, y); i++)
m_domain.insert(i);
#endif
lp_assert(x <= y);
if (x == y) {
unite_with_one_point_interval(x);
return;
}
const_iter l, r;
bool neg_inf, pos_inf;
bool found_left_point = get_left_point(x, neg_inf, l);
bool found_right_point = get_right_point(y, pos_inf, r);
m_empty = false;
if (!found_left_point) {
if (!found_right_point) {
m_endpoints.clear();
set_start(x);
set_end(y);
return;
}
// found_right_point is true
if (pos_inf) {
m_endpoints.clear();
set_start(x);
return;
}
remove_from_the_left(y);
if (pos(m_endpoints.begin()) == y || pos(m_endpoints.begin()) == y + 1) {
if (is_proper_start(m_endpoints.begin()))
m_endpoints.erase(m_endpoints.begin());
else
set_end(pos(m_endpoints.begin()));
set_start(x);
}
else {
lp_assert(pos(m_endpoints.begin()) > y + 1);
if (is_start(m_endpoints.begin()))
set_end(y);
set_start(x);
}
return;
}
lp_assert(found_left_point);
if (!found_right_point) {
bool neg_inf = has_neg_inf();
remove_from_the_right(x);
if (m_endpoints.empty()) {
if (!neg_inf)
set_start_end(x, y);
else
set_end(y);
return;
}
if (pos(m_endpoints.rbegin()) == x || pos(m_endpoints.rbegin()) == x - 1) {
if (is_proper_end(m_endpoints.rbegin())) {
m_endpoints.erase(m_endpoints.rbegin());
}
else if (is_one_point_interval(m_endpoints.rbegin())) {
set_start(pos(m_endpoints.rbegin()));
}
set_end(y);
}
else {
if (is_end(m_endpoints.rbegin())) {
set_start(x);
set_end(y);
}
else {
set_end(y);
}
}
return;
}
// found_right_point and found_left_point
if (!neg_inf)
handle_left_point_in_union(l, x, y);
else {
remove_from_the_left(y);
}
if (!pos_inf)
handle_right_point_in_union(r, y);
else
remove_from_the_right(x);
}
bool has_pos_inf() const {
if (m_empty)
return false;
if (m_endpoints.empty())
return true;
lp_assert(m_endpoints.rbegin() != m_endpoints.rend());
return m_endpoints.rbegin()->second == 0;
}
bool has_neg_inf() const {
if (m_empty)
return false;
if (m_endpoints.empty())
return true;
auto it = m_endpoints.begin();
return is_proper_end(it->second);//m_endpoints.begin());
}
bool is_correct() const {
if (m_empty) {
if (m_endpoints.size() > 0) {
TRACE("disj_intervals", tout << "is empty is true but m_endpoints.size() = " << m_endpoints.size() << std::endl;);
return false;
}
return true;
}
bool expect_end;
bool prev = false;
T prev_x;
for (auto t : m_endpoints()) {
if (prev && ((expect_end && !is_end(t.second)) || (!expect_end && !is_start(t.second)))) {
TRACE("disj_intervals", tout << "x = " << t.first << "\n";);
if (expect_end) {
TRACE("disj_intervals", tout << "expecting an interval end\n";);
} else {
TRACE("disj_intervals", tout << "expecting an interval start\n";);
}
return false;
}
if (prev) {
if (t.first - prev_x <= 1 && !expect_end) {
TRACE("disj_intervals", tout << "the sequence is not increasing or the gap is too small: " << prev_x << ", " << t.first << std::endl;);
return false;
}
}
if (t.second == 2) {
expect_end = false; // swallow a point interval
}
else {
if (prev)
expect_end = !expect_end;
else
expect_end = is_start(t.second);
}
prev = true;
prev_x = t.first;
}
#if Z3DEBUG
for (int i = 0; i <= 100; i++ ) {
if ( (m_domain.find(i) != m_domain.end()) != contains(i)) {
TRACE("disj_intervals", tout << "incorrect value of contains(" << i << ") is = " << contains(i) << std::endl;);
return false;
}
}
#endif
return true;
}
public:
void print(std::ostream & out) const {
if (m_empty) {
out << "empty\n";
return;
}
if (m_endpoints.empty()){
out << "[-oo,oo]\n";
return;
}
bool first = true;
for (auto t : m_endpoints()) {
if (first) {
if (t.second == 1) {
out << "[-oo," << t.first << "]";
}
else if (t.second == 0)
out << "[" << t.first << ",";
else if (t.second == 2)
out << "[" << t.first << "]";
first = false;
} else {
if (t.second==0)
out << "[" << t.first << ",";
else if (t.second == 1)
out << t.first << "]";
else if (t.second == 2)
out << "[" << t.first << "]";
}
}
if (has_pos_inf())
out << "oo]";
out << "\n";
}
void push() { m_endpoints.push(); m_empty.push(); }
void pop() { m_endpoints.pop(); m_empty.pop(); }
void pop(unsigned k) { while(k--) pop(); }
// we intersect the existing set with the half open to the right interval
void intersect_with_lower_bound(const T& x) {
#ifdef Z3DEBUG
for (int i = 0; i < x; i++)
m_domain.erase(i);
#endif
TRACE("disj_intervals", tout << "intersect_with_lower_bound(" << x << ")\n";);
if (m_empty)
return;
if (m_endpoints.empty()) {
set_start(x);
return;
}
bool pos_inf = has_pos_inf();
auto it = m_endpoints.begin();
while (it != m_endpoints.end() && pos(it) < x) {
m_endpoints.erase(it);
it = m_endpoints.begin();
}
if (m_endpoints.empty()) {
if (!pos_inf) {
m_empty = true;
return;
}
set_start(x);
return;
}
lp_assert(pos(it) >= x);
if (pos(it) == x) {
if (is_proper_end(it))
set_one_point_interval(x);
}
else { // x(it) > x
if (is_proper_end(it)) {
set_start(x);
}
}
}
// we intersect the existing set with the half open interval
void intersect_with_upper_bound(const T& x) {
#ifdef Z3DEBUG
for (int i = 100; i > x; i--)
m_domain.erase(i);
#endif
TRACE("disj_intervals", tout << "intersect_with_upper_bound(" << x << ")\n";);
if (m_empty)
return;
if (m_endpoints.empty()) {
set_end(x);
return;
}
bool neg_inf = has_neg_inf();
auto it = m_endpoints.rbegin();
while (!m_endpoints.empty() && pos(it) > x) {
m_endpoints.erase(std::prev(m_endpoints.end()));
it = m_endpoints.rbegin();
}
if (m_endpoints.empty()) {
if (!neg_inf) {
m_empty = true;
return;
}
set_end(x);
}
lp_assert(pos(it) <= x);
if (pos(it) == x) {
if (is_one_point_interval(it)) {}
else if (is_proper_end(it)) {}
else {// is_proper_start(it->second)
set_one_point_interval(x);
}
}
else { // pos(it) < x}
if (is_proper_start(it))
set_end(x);
}
lp_assert(is_correct());
}
public:
void intersect_with_interval(const T& x, const T & y) {
#ifdef Z3DEBUG
for (int i = 0; i <= 100; i++)
if (i < x || i > y)
m_domain.erase(i);
#endif
TRACE("disj_intervals", tout << "intersect_with_interval(" << x << ", " << y <<")\n";);
if (m_empty)
return;
lp_assert(x <= y);
intersect_with_lower_bound(x);
intersect_with_upper_bound(y);
}
// add an intervar [x, inf]
void unite_with_interval_x_pos_inf(const T& x) {
if (contains_all())
return;
#if Z3DEBUG
for (int i = x; i <= 100; i++)
m_domain.insert(i);
#endif
TRACE("disj_intervals", tout << "unite_with_interval_x_pos_inf(" << x << ")\n";);
if (m_empty) {
set_start(x);
m_empty = false;
return;
}
bool neg_inf = has_neg_inf();
remove_from_the_right(x);
if (m_endpoints.empty()) {
if (!neg_inf)
set_start(x);
return;
}
auto it = m_endpoints.rbegin();
lp_assert(pos(it) <= x);
if (pos(it) == x) {
if (is_proper_end(it)) {
m_endpoints.erase(x);
} else {
set_start(x);
}
} else if (pos(it) == x - 1 && is_end(it)) {
if (is_proper_start(it)) {
// do nothing
}
else if (is_proper_end(it)) {
m_endpoints.erase(it);
}
else {
lp_assert(is_one_point_interval(it));
set_start(it);
}
} else {
if (!has_pos_inf())
set_start(x);
}
}
// add an interval [-inf, x]
void unite_with_interval_neg_inf_x(const T& x) {
#if Z3DEBUG
for (int i = 0; i <= x; i++)
m_domain.insert(i);
#endif
TRACE("disj_intervals", tout << "unite_with_interval_neg_inf_x(" << x << ")\n";);
if (m_empty) {
set_end(x);
m_empty = false;
return;
}
bool pos_inf;
const_iter r;
bool found_right_point = get_right_point(x, pos_inf, r);
if (!found_right_point) {
m_endpoints.clear();
set_end(x);
return;
}
if (pos_inf) {
m_endpoints.clear();
return;
}
lp_assert(pos(r) >= x);
if (pos(r) == x || pos(r) == x + 1) {
if (is_proper_start(r))
erase(r);
else if (is_one_point_interval(r)) {
set_end(pos(r));
} // do nothing for the proper end
} else {
if (!is_proper_end(r))
set_end(x);
}
while (!m_endpoints.empty() && m_endpoints.begin()->first < x) {
m_endpoints.erase(m_endpoints.begin());
}
lp_assert(is_correct());
}
private:
bool is_start(char x) const { return x == 0 || x == 2; }
bool is_start(const iter & it) const { return is_start(it->second); }
bool is_start(const const_iter & it) const { return is_start(it->second); }
bool is_start(const riter & it) const {
return is_start(it->second);
}
bool is_end(char x) const { return x == 1 || x == 2; }
bool is_end(const iter & it) const {
return is_end(it->second);
}
bool is_end(const const_iter & it) const {
return is_end(it->second);
}
bool is_end(const riter & it) const {
return is_end(it->second);
}
T pos(const iter & it) const {
return it->first;
}
T pos(const const_iter & it) const {
return it->first;
}
T pos(const riter & it) const {
return it->first;
}
T bound_kind(iter & it) const {
return it->second;
}
T bound_kind(riter & it) const {
return it->second;
}
bool is_proper_start(char x) const { return x == 0; }
bool is_proper_start(const riter &x) const { return is_proper_start(x->second);}
bool is_proper_start(const iter &x) const { return is_proper_start(x->second);}
bool is_proper_start(const const_iter &x) const { return is_proper_start(x->second);}
bool is_proper_end(char x) const { return x == 1; }
bool is_proper_end(const iter & it) const { return is_proper_end(it->second); }
bool is_proper_end(const const_iter & it) const { return is_proper_end(it->second); }
bool is_proper_end(const riter & it) const {
return is_proper_end(it->second);
}
bool is_one_point_interval(char x) const { return x == 2; }
bool is_one_point_interval(const iter & it) const {
return is_one_point_interval(it->second);
}
bool is_one_point_interval(const const_iter & it) const {
return is_one_point_interval(it->second);
}
bool is_one_point_interval(const riter & it) const {
return is_one_point_interval(it->second);
}
void erase(const iter& it) {
m_endpoints.erase(it);
}
void erase(const const_iter& it) {
m_endpoints.erase(it);
}
void erase(const riter& it) {
m_endpoints.erase(it);
}
void erase(const T& x) {
m_endpoints.erase(x);
}
void set_one_point_interval(const T& x) {
m_endpoints[x] = 2;
}
void set_start(const T& x) {
m_endpoints[x] = 0;
}
void set_start(const iter &t ) {
t->second = 0;
}
void set_start(riter &t ) {
t->second = 0;
}
void set_end(const T& x) {
m_endpoints[x] = 1;
}
void set_end(const iter& t) {
t->second = 1;
}
void set_end(riter& t) {
t->second = 1;
}
private:
void set_start_end(const T& x, const T & y) {
set_start(x);
set_end(y);
}
void unite_with_one_point_interval(const T &x) {
TRACE("disj_intervals", tout << "unite_with_one_point_interval(" << x << ")\n";);
if (m_empty) {
m_empty = false;
set_one_point_interval(x);
return;
}
bool has_left_neig, has_right_neig;
bool neg_inf, pos_inf;
const_iter l, r;
has_left_neig = get_left_point(x, neg_inf, l);
has_right_neig = get_right_point(x, pos_inf, r);
if (!has_left_neig) {
if (!has_right_neig) {
set_one_point_interval(x);
return;
}
lp_assert(!neg_inf);
// if (pos(r) == x ) nothing happens
if (pos(r) == x + 1) {
if (is_one_point_interval(r)) {
set_start(x);
set_end(x + 1);
} else {
lp_assert(is_proper_start(r));
erase(r);
set_start(x);
}
} else {
lp_assert(pos(r) > x);
set_one_point_interval(x);
}
return;
}
// has_left_neig
if (!has_right_neig) {
if (pos_inf)
return;
// if (pos(l) == x nothing happens
if (pos(l) == x - 1) {
if (is_proper_end(l)) {
erase(l);
set_end(x);
} else {
lp_assert(is_one_point_interval(l));
set_start(l->first);
set_end(x);
}
} else {
lp_assert(pos(l) < x - 1);
set_one_point_interval(x);
}
return;
}
// has both neighbors
if (neg_inf || pos_inf)
return;
if (pos(l) == x || pos(r) == x)
return;
// now the cases pos(l) == pos(r) or pos(l) + 1 == pos(r) are impossible
if (is_proper_start(l)) {
lp_assert(is_proper_end(r));
return;
}
if (pos(l) + 2 < pos(r)) { // we can glue only to one neighbor
if (pos(l) + 1 == x) {
if (is_proper_end(l)) {
erase(l);
set_end(x);
}
else {
lp_assert(!is_proper_start(l));
set_start(pos(l));
set_end(x);
}
}
else if (x + 1 == pos(r)) {
if (is_proper_start(r)) {
erase(r);
}
else if (is_one_point_interval(r)) {
set_end(pos(r));
}
set_start(x);
}
else {
set_one_point_interval(x);
}
} else {
lp_assert(pos(l) + 2 == pos(r));
lp_assert(pos(l) + 1 == x); // x is just in between l and r
if (is_proper_end(l)) {
erase(l);
} else {
lp_assert(is_one_point_interval(l));
set_start(l->first);
}
if (is_proper_start(r)) {
erase(r);
} else {
lp_assert(is_one_point_interval(r));
set_end(r->first);
}
}
}
// return false if there are no points y in m_endpoints with pos(y) <= x
// and negative infiniti is not true
bool get_left_point(const T& x, bool &neg_inf, const_iter & l) const {
if (m_empty) {
neg_inf = false;
return false;
}
if (m_endpoints.empty()) {
neg_inf = true;
return true;
}
l = m_endpoints.lower_bound(x);
if (l == m_endpoints.end()) {
l--;
neg_inf = false;
return true;
}
if (pos(l) == x) {
neg_inf = false;
return true;
}
if (l == m_endpoints.begin())
return neg_inf = has_neg_inf();
l--;
neg_inf = false;
return true;
}
// return false iff there are no points y in m_endpoints with pos(y) >= x, and positive infinity is not true
bool get_right_point(const T& x, bool &pos_inf, const_iter & r) const {
if (m_empty) {
pos_inf = false;
return false;
}
if (m_endpoints.empty()) {
pos_inf = true;
return true;
}
r = m_endpoints.lower_bound(x);
if (r == m_endpoints.end()) {
return pos_inf = has_pos_inf();;
}
pos_inf = false;
return true;
}
// inserts x, if needed and return a pointer to the leftmost point that is a candidate to removal
/*
iter insert_start_for_unite_with_interval(const T& x) {
auto lbx = m_endpoints.lower_bound(x);
if (lbx == m_endpoints.end()) {
if (!has_pos_inf()) {
T lastx = pos(m_endpoints.rbegin());
if (lastx + 1 < x)
set_start_end(x, y);
else {
m_endpoints.erase(lastx);
set_end(y);
}
}
return;
}
if (pos(lbx) == x) {
if(is_proper_end(lbx)) {
m_endpoints.erase(x);
} else { set_start(x);}
}
if (pos(lbx) > x) {
if (pos(lbx) > y + 1) {
if (is_end(lbx))
return;
set_start_end(x, y);
return;
}
if (pos(lpx) == y + 1) {
if (is_end(lbx))
}
}
lp_assert(false); // not implemented
}
*/
void remove_from_the_left(const T & y ) {
while (!m_endpoints.empty() && pos(m_endpoints.begin()) < y) {
m_endpoints.erase(m_endpoints.begin());
}
}
void remove_from_the_left(const T & y, const_iter& l ) {
while (l!= m_endpoints.end() && pos(l) < y) {
l++;
erase(std::prev(l));
}
}
void remove_from_the_right(const T & x) {
while (!m_endpoints.empty() && pos(m_endpoints.rbegin()) > x) {
m_endpoints.erase(m_endpoints.rbegin());
}
}
void remove_from_the_right(const T & x, riter & r) {
while (!m_endpoints.empty() && pos(r) > x) {
r--;
m_endpoints.erase(std::next(r));
}
}
public:
bool get_lower_bound(T& b) const {
if (m_empty)
return false;
if (has_neg_inf())
return false;
b = pos(m_endpoints.begin());
return true;
}
bool get_upper_bound(T& b) const {
if (m_empty)
return false;
if (has_pos_inf())
return false;
b = m_endpoints.rbegin()->first;
return true;
}
};
}

3
src/util/lp/lar_constraints.h
View file

@ -44,8 +44,7 @@ inline std::string lconstraint_kind_string(lconstraint_kind t) {
return std::string(); // it is unreachable
}
class lar_base_constraint {
public:
struct lar_base_constraint {
lconstraint_kind m_kind;
mpq m_right_side;
virtual vector<std::pair<mpq, var_index>> get_left_side_coefficients() const = 0;

3
src/util/lp/lar_core_solver.h
View file

@ -346,7 +346,7 @@ public:
unsigned jb = m_r_solver.m_basis[i];
m_r_solver.update_x_with_delta_and_track_feasibility(jb, - delta * m_r_solver.m_A.get_val(cc));
}
lp_assert(m_r_solver.A_mult_x_is_off() == false);
CASSERT("A_off", m_r_solver.A_mult_x_is_off() == false);
}
lp_assert(m_r_solver.inf_set_is_correct());
}
@ -781,7 +781,6 @@ public:
if (m_column_types()[j] != column_type::boxed)
continue;
update_delta(delta, m_r_low_bounds[j], m_r_upper_bounds[j]);
}
return delta;
}

4
src/util/lp/lar_core_solver.hpp
View file

@ -222,7 +222,7 @@ void lar_core_solver::calculate_pivot_row(unsigned i) {
}
void lar_core_solver::fill_not_improvable_zero_sum_from_inf_row() {
lp_assert(m_r_solver.A_mult_x_is_off() == false);
CASSERT("A_off", m_r_solver.A_mult_x_is_off() == false);
unsigned bj = m_r_basis[m_r_solver.m_inf_row_index_for_tableau];
m_infeasible_sum_sign = m_r_solver.inf_sign_of_column(bj);
m_infeasible_linear_combination.clear();
@ -273,7 +273,7 @@ void lar_core_solver::solve() {
return;
}
++settings().st().m_need_to_solve_inf;
lp_assert(!m_r_solver.A_mult_x_is_off());
CASSERT("A_off", !m_r_solver.A_mult_x_is_off());
lp_assert((!settings().use_tableau()) || r_basis_is_OK());
if (need_to_presolve_with_double_solver()) {
prefix_d();

954
src/util/lp/lar_solver.cpp
View file

File diff suppressed because it is too large Load diff

75
src/util/lp/lar_solver.h
View file

@ -34,7 +34,6 @@ Revision History:
#include "util/lp/lp_primal_core_solver.h"
#include "util/lp/random_updater.h"
#include <stack>
#include "util/lp/stacked_map.h"
#include "util/lp/stacked_value.h"
#include "util/lp/stacked_vector.h"
#include "util/lp/stacked_unordered_set.h"
@ -67,6 +66,11 @@ class lar_solver : public column_namer {
vector<unsigned> m_columns_to_ext_vars_or_term_indices;
stacked_vector<ul_pair> m_columns_to_ul_pairs;
vector<lar_base_constraint*> m_constraints;
public :
const vector<lar_base_constraint*>& constraints() const {
return m_constraints;
}
private:
stacked_value<unsigned> m_constraint_count;
// the set of column indices j such that bounds have changed for j
@ -160,33 +164,12 @@ public:
unsigned adjust_term_index(unsigned j) const;
void analyze_new_bounds_on_row(
unsigned row_index,
lp_bound_propagator & bp) {
SASSERT(!use_tableau());
iterator_on_pivot_row<mpq> it(m_mpq_lar_core_solver.get_pivot_row(), m_mpq_lar_core_solver.m_r_basis[row_index]);
bound_analyzer_on_row ra_pos(it,
zero_of_type<numeric_pair<mpq>>(),
row_index,
bp
);
ra_pos.analyze();
}
unsigned row_index,
bound_propagator & bp);
void analyze_new_bounds_on_row_tableau(
unsigned row_index,
lp_bound_propagator & bp
) {
if (A_r().m_rows[row_index].size() > settings().max_row_length_for_bound_propagation)
return;
iterator_on_row<mpq> it(A_r().m_rows[row_index]);
SASSERT(use_tableau());
bound_analyzer_on_row::analyze_row(it,
zero_of_type<numeric_pair<mpq>>(),
row_index,
bp
);
unsigned row_index,
bound_propagator & bp);
void substitute_basis_var_in_terms_for_row(unsigned i) {
// todo : create a map from term basic vars to the rows where they are used
@ -588,7 +571,7 @@ public:
void substitute_terms_in_linear_expression( const vector<std::pair<mpq, var_index>>& left_side_with_terms,
vector<std::pair<mpq, var_index>> &left_side, mpq & free_coeff) const;
vector<std::pair<mpq, var_index>> &left_side, mpq & free_coeff) const;
void detect_rows_of_bound_change_column_for_nbasic_column(unsigned j) {
@ -1080,9 +1063,9 @@ public:
}
void get_infeasibility_explanation_for_inf_sign(
vector<std::pair<mpq, constraint_index>> & explanation,
const vector<std::pair<mpq, unsigned>> & inf_row,
int inf_sign) const;
vector<std::pair<mpq, constraint_index>> & explanation,
const vector<std::pair<mpq, unsigned>> & inf_row,
int inf_sign) const;
int adj_sign = coeff.is_pos() ? inf_sign : -inf_sign;
@ -1173,10 +1156,10 @@ public:
bool column_represents_row_in_tableau(unsigned j);
void make_sure_that_the_bottom_right_elem_not_zero_in_tableau(unsigned i, unsigned j);
void remove_last_row_and_column_from_tableau(unsigned j);
void remove_last_column_from_tableau(unsigned j);
void remove_last_column_from_A();
void remove_last_column_from_basis_tableau(unsigned j);
void remove_column_from_tableau(unsigned j);
void remove_last_column_from_tableau();
void pop_tableau();
void clean_inf_set_of_r_solver_after_pop();
void shrink_explanation_to_minimum(vector<std::pair<mpq, constraint_index>> & explanation) const;
@ -1189,7 +1172,7 @@ public:
return !column_is_int(j);
}
bool model_is_int_feasible() const;
bool model_is_int_feasible() const;
void remove_last_row_and_column_from_tableau(unsigned j) {
@ -1363,16 +1346,16 @@ bool model_is_int_feasible() const;
}
bool inf_int_set_is_correct_for_column(unsigned j) const {
if (m_inf_int_set.contains(j) != (column_is_int(j) && (!column_value_is_integer(j)))) {
TRACE("arith_int",
tout << "j= " << j <<
" inf_int_set().contains(j) = " << m_inf_int_set.contains(j) <<
", column_is_int(j) = " << column_is_int(j) <<
"\n column_value_is_integer(j) = " << column_value_is_integer(j) <<
", val = " << get_column_value(j) << std::endl;);
return false;
}
return true;
if (m_inf_int_set.contains(j) != (column_is_int(j) && (!column_value_is_integer(j)))) {
TRACE("arith_int",
tout << "j= " << j <<
" inf_int_set().contains(j) = " << m_inf_int_set.contains(j) <<
", column_is_int(j) = " << column_is_int(j) <<
"\n column_value_is_integer(j) = " << column_value_is_integer(j) <<
", val = " << get_column_value(j) << std::endl;);
return false;
}
return true;
}
bool inf_int_set_is_correct() const {
@ -1383,10 +1366,7 @@ bool model_is_int_feasible() const;
return false;
}
return true;
}
}
bool has_int_var() const;
void call_assignment_tracker(unsigned j) {
if (!var_is_int(j)) {
@ -1400,5 +1380,6 @@ bool model_is_int_feasible() const;
}
lar_core_solver & get_core_solver() { return m_mpq_lar_core_solver; }
};
}

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

@ -89,6 +89,7 @@ struct lar_term {
m_coeffs.clear();
m_v = zero_of_type<mpq>();
}
};
}

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

@ -603,7 +603,12 @@ public:
lp_assert(false);
}
// out << "basis heading = " << m_basis_heading[j] << std::endl;
out << " x = " << m_x[j] << std::endl;
out << " x = " << m_x[j];
if (m_basis_heading[j] >= 0)
out << " base\n";
else
out << " nbas\n";
/*
std::cout << "cost = " << m_costs[j] << std::endl;
std:: cout << "m_d = " << m_d[j] << std::endl;*/

6
src/util/lp/lp_core_solver_base.hpp
View file

@ -226,7 +226,7 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::
A_mult_x_is_off() const {
lp_assert(m_x.size() == m_A.column_count());
if (numeric_traits<T>::precise()) {
for (unsigned i = 0; i < m_m(); i++) {
for (unsigned i = 0; i < m_m(); i++) {
X delta = m_b[i] - m_A.dot_product_with_row(i, m_x);
if (delta != numeric_traits<X>::zero()) {
std::cout << "precise x is off (";
@ -591,7 +591,7 @@ update_basis_and_x(int entering, int leaving, X const & tt) {
restore_x_and_refactor(entering, leaving, tt);
if (m_status == lp_status::FLOATING_POINT_ERROR)
return false;
lp_assert(!A_mult_x_is_off());
CASSERT("A_off", !A_mult_x_is_off());
m_iters_with_no_cost_growing++;
// LP_OUT(m_settings, "rolled back after failing of init_factorization()" << std::endl);
m_status = lp_status::UNSTABLE;
@ -630,7 +630,7 @@ divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col) {
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::
pivot_column_tableau(unsigned j, unsigned piv_row_index) {
if (!divide_row_by_pivot(piv_row_index, j))
if (!divide_row_by_pivot(piv_row_index, j))
return false;
auto &column = m_A.m_columns[j];
int pivot_col_cell_index = -1;

2
src/util/lp/lp_primal_core_solver.hpp
View file

@ -673,7 +673,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::calc_work
template <typename T, typename X>
void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving(int entering, X & t) {
lp_assert(!this->A_mult_x_is_off() );
CASSERT("A_off", !this->A_mult_x_is_off() );
this->update_x(entering, t * m_sign_of_entering_delta);
if (this->A_mult_x_is_off_on_index(this->m_ed.m_index) && !this->find_x_by_solving()) {
this->init_lu();

10
src/util/lp/lp_primal_core_solver_tableau.h
View file

@ -206,7 +206,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
}
template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_entering_and_leaving_tableau(int entering, int leaving, X & t) {
lp_assert(this->A_mult_x_is_off() == false);
CASSERT("A_off", this->A_mult_x_is_off() == false);
lp_assert(leaving >= 0 && entering >= 0);
lp_assert((this->m_settings.simplex_strategy() ==
simplex_strategy_enum::tableau_rows) ||
@ -223,7 +223,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
t = -t;
}
this->update_basis_and_x_tableau(entering, leaving, t);
lp_assert(this->A_mult_x_is_off() == false);
CASSERT("A_off", this->A_mult_x_is_off() == false);
this->iters_with_no_cost_growing() = 0;
} else {
this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
@ -247,7 +247,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
template <typename T, typename X>
void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving_tableau(int entering, X & t) {
lp_assert(!this->A_mult_x_is_off() );
CASSERT("A_off", !this->A_mult_x_is_off() );
this->update_x_tableau(entering, t * m_sign_of_entering_delta);
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
return;
@ -320,7 +320,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tableau() {
// print_matrix(&(this->m_A), std::cout);
lp_assert(this->A_mult_x_is_off() == false);
CASSERT("A_off", this->A_mult_x_is_off() == false);
lp_assert(basis_columns_are_set_correctly());
this->m_basis_sort_counter = 0; // to initiate the sort of the basis
this->set_total_iterations(0);
@ -374,7 +374,7 @@ update_x_tableau(unsigned entering, const X& delta) {
this->insert_column_into_inf_set(j);
}
}
lp_assert(this->A_mult_x_is_off() == false);
CASSERT("A_off", this->A_mult_x_is_off() == false);
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::

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

@ -23,7 +23,7 @@ Revision History:
#include "util/debug.h"
#include <unordered_map>
template <typename A, typename B>
bool try_get_val(const std::unordered_map<A,B> & map, const A& key, B & val) {
bool try_get_value(const std::unordered_map<A,B> & map, const A& key, B & val) {
const auto it = map.find(key);
if (it == map.end()) return false;
val = it->second;

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

@ -25,7 +25,6 @@ Revision History:
#include "../rational.h"
#include "../sstream.h"
#include "../z3_exception.h"
#else
// include "util/numerics/mpq.h"
// include "util/numerics/numeric_traits.h"
@ -50,8 +49,21 @@ public:
static unsigned one() { return 1; }
static bool is_zero(unsigned v) { return v == 0; }
static double get_double(unsigned const & d) { return d; }
static bool is_int(unsigned) {return true;}
};
template <> class numeric_traits<int> {
public:
static bool precise() { return true; }
static int const zero() { return 0; }
static int const one() { return 1; }
static bool is_zero(int v) { return v == 0; }
static double const get_double(int const & d) { return d; }
static bool is_int(int) {return true;}
static bool is_pos(int j) {return j > 0;}
};
template <> class numeric_traits<double> {
public:
static bool precise() { return false; }
@ -79,6 +91,7 @@ template <> class numeric_traits<double> {
static rational from_string(std::string const & str) { return rational(str.c_str()); }
static bool is_pos(const rational & d) {return d.is_pos();}
static bool is_neg(const rational & d) {return d.is_neg();}
static bool is_int(const rational & d) {return d.is_int();}
};
template <typename X, typename Y>

73
src/util/lp/stacked_map.h
View file

@ -21,29 +21,28 @@ Revision History:
#pragma once
// this class implements a map with some stack functionality
#include <unordered_map>
#include <unordered_set>
#include <set>
#include <stack>
namespace lp {
template <typename A, typename B,
typename Hash = std::hash<A>,
typename KeyEqual = std::equal_to<A>,
typename Allocator = std::allocator< std::pair<const A, B> > >
template<class A,
class B,
class _Pr = std::less<A>,
class _Alloc = std::allocator<std::pair<const A, B> > >
class stacked_map {
struct delta {
std::unordered_set<A, Hash, KeyEqual> m_new;
std::unordered_map<A, B, Hash, KeyEqual, Allocator> m_original_changed;
std::unordered_set<A> m_new;
std::unordered_map<A, B> m_original_changed;
// std::unordered_map<A,B, Hash, KeyEqual, Allocator > m_deb_copy;
};
std::unordered_map<A,B,Hash, KeyEqual, Allocator> m_map;
std::map<A,B,_Pr, _Alloc> m_map;
std::stack<delta> m_stack;
public:
class ref {
stacked_map<A,B,Hash, KeyEqual, Allocator> & m_map;
stacked_map<A,B,_Pr, _Alloc> & m_map;
const A & m_key;
public:
ref(stacked_map<A,B,Hash, KeyEqual, Allocator> & m, const A & key) :m_map(m), m_key(key) {}
ref(stacked_map<A,B,_Pr, _Alloc> & m, const A & key) :m_map(m), m_key(key) {}
ref & operator=(const B & b) {
m_map.emplace_replace(m_key, b);
return *this;
@ -55,8 +54,44 @@ public:
return it->second;
}
};
typename std::map < A, B, _Pr, _Alloc>::iterator upper_bound(const A& k) {
return m_map.upper_bound(k);
}
typename std::map < A, B, _Pr, _Alloc>::const_iterator upper_bound(const A& k) const {
return m_map.upper_bound(k);
}
typename std::map < A, B, _Pr, _Alloc>::iterator lower_bound(const A& k) {
return m_map.lower_bound(k);
}
typename std::map < A, B, _Pr, _Alloc>::const_iterator lower_bound(const A& k) const {
return m_map.lower_bound(k);
}
typename std::map < A, B, _Pr, _Alloc>::iterator end() {
return m_map.end();
}
typename std::map < A, B, _Pr, _Alloc>::const_iterator end() const {
return m_map.end();
}
typename std::map < A, B, _Pr, _Alloc>::reverse_iterator rend() {
return m_map.rend();
}
typename std::map < A, B, _Pr, _Alloc>::const_reverse_iterator rend() const {
return m_map.rend();
}
typename std::map < A, B, _Pr, _Alloc>::iterator begin() {
return m_map.begin();
}
typename std::map < A, B, _Pr, _Alloc>::const_iterator begin() const {
return m_map.begin();
}
typename std::map < A, B, _Pr, _Alloc>::reverse_iterator rbegin() {
return m_map.rbegin();
}
typename std::map < A, B, _Pr, _Alloc>::const_reverse_iterator rbegin() const {
return m_map.rbegin();
}
private:
void emplace_replace(const A & a, const B & b) {
void emplace_replace(const A & a, const B & b) {
if (!m_stack.empty()) {
delta & d = m_stack.top();
auto it = m_map.find(a);
@ -147,7 +182,15 @@ public:
m_stack.pop();
}
}
void erase(const typename std::map < A, B, _Pr, _Alloc>::iterator & it) {
erase(it->first);
}
void erase(const typename std::map < A, B, _Pr, _Alloc>::const_iterator & it) {
erase(it->first);
}
void erase(const typename std::map < A, B, _Pr, _Alloc>::reverse_iterator & it) {
erase(it->first);
}
void erase(const A & key) {
if (m_stack.empty()) {
m_map.erase(key);
@ -189,6 +232,8 @@ public:
m_map.clear();
}
const std::unordered_map<A, B,Hash, KeyEqual, Allocator>& operator()() const { return m_map;}
bool empty() const { return m_map.empty(); }
const std::map<A, B, _Pr, _Alloc>& operator()() const { return m_map;}
};
}

11
src/util/lp/static_matrix.h
View file

@ -62,7 +62,6 @@ class static_matrix
dim(unsigned m, unsigned n) :m_m(m), m_n(n) {}
};
std::stack<dim> m_stack;
vector<unsigned> m_became_zeros; // the row indices that became zeroes during the pivoting
public:
typedef vector<row_cell<T>> row_strip;
typedef vector<column_cell> column_strip;
@ -97,6 +96,14 @@ public:
const T & get_val(const column_cell & c) const {
return m_rows[c.m_i][c.m_offset].get_val();
}
row_cell<T> & get_row_cell(const column_cell & c) {
return m_rows[c.m_i][c.m_offset];
}
column_cell & get_column_cell(const row_cell<T> &rc) {
return m_columns[rc.m_j][rc.m_offset];
}
void init_row_columns(unsigned m, unsigned n);
@ -287,7 +294,7 @@ public:
// pivot row i to row ii
bool pivot_row_to_row_given_cell(unsigned i, column_cell& c, unsigned);
void scan_row_ii_to_offset_vector(unsigned ii);
void scan_row_ii_to_offset_vector(const row_strip & rvals);
void transpose_rows(unsigned i, unsigned ii) {
auto t = m_rows[i];

57
src/util/lp/static_matrix.hpp
View file

@ -35,57 +35,46 @@ void static_matrix<T, X>::init_row_columns(unsigned m, unsigned n) {
}
template <typename T, typename X> void static_matrix<T, X>::scan_row_ii_to_offset_vector(unsigned ii) {
auto & rvals = m_rows[ii];
unsigned size = rvals.size();
for (unsigned j = 0; j < size; j++)
template <typename T, typename X> void static_matrix<T, X>::scan_row_ii_to_offset_vector(const row_strip & rvals) {
for (unsigned j = 0; j < rvals.size(); j++)
m_vector_of_row_offsets[rvals[j].m_j] = j;
}
template <typename T, typename X> bool static_matrix<T, X>::pivot_row_to_row_given_cell(unsigned i, column_cell & c, unsigned pivot_col) {
unsigned ii = c.m_i;
lp_assert(i < row_count() && ii < column_count());
lp_assert(i != ii);
m_became_zeros.reset();
lp_assert(i < row_count() && ii < column_count() && i != ii);
T alpha = -get_val(c);
lp_assert(!is_zero(alpha));
auto & ii_row_vals = m_rows[ii];
remove_element(ii_row_vals, ii_row_vals[c.m_offset]);
scan_row_ii_to_offset_vector(ii);
lp_assert(!is_zero(alpha));
unsigned prev_size_ii = ii_row_vals.size();
auto & rowii = m_rows[ii];
remove_element(rowii, rowii[c.m_offset]);
scan_row_ii_to_offset_vector(rowii);
unsigned prev_size_ii = rowii.size();
// run over the pivot row and update row ii
for (const auto & iv : m_rows[i]) {
unsigned j = iv.m_j;
if (j == pivot_col) continue;
T alv = alpha * iv.m_value;
lp_assert(!is_zero(iv.m_value));
lp_assert(!is_zero(iv.m_value));
int j_offs = m_vector_of_row_offsets[j];
if (j_offs == -1) { // it is a new element
add_new_element(ii, j, alv);
}
else {
auto & row_el_iv = ii_row_vals[j_offs];
row_el_iv.m_value += alv;
if (is_zero(row_el_iv.m_value)) {
m_became_zeros.push_back(j_offs);
ensure_increasing(m_became_zeros);
}
rowii[j_offs].m_value += alv;
}
}
// clean the work vector
for (unsigned k = 0; k < prev_size_ii; k++) {
m_vector_of_row_offsets[ii_row_vals[k].m_j] = -1;
m_vector_of_row_offsets[rowii[k].m_j] = -1;
}
for (unsigned k = m_became_zeros.size(); k-- > 0; ) {
unsigned j = m_became_zeros[k];
remove_element(ii_row_vals, ii_row_vals[j]);
// remove zeroes
for (unsigned k = rowii.size(); k-- > 0; ) {
if (is_zero(rowii[k].m_value))
remove_element(rowii, rowii[k]);
}
return !ii_row_vals.empty();
return !rowii.empty();
}
@ -353,11 +342,12 @@ template <typename T, typename X> T static_matrix<T, X>::get_row_balance(unsi
}
template <typename T, typename X> bool static_matrix<T, X>::is_correct() const {
for (auto & r : m_rows) {
for (unsigned i = 0; i < m_rows.size(); i++) {
auto &r = m_rows[i];
std::unordered_set<unsigned> s;
for (auto & rc : r) {
if (s.find(rc.m_j) != s.end()) {
std::cout << "found column " << rc.m_j << " twice in a row" << std::endl;
std::cout << "found column " << rc.m_j << " twice in a row " << i << std::endl;
return false;
}
s.insert(rc.m_j);
@ -367,13 +357,20 @@ template <typename T, typename X> bool static_matrix<T, X>::is_correct() const {
return false;
if (rc.get_val() != get_val(m_columns[rc.m_j][rc.m_offset]))
return false;
if (is_zero(rc.get_val())) {
std::cout << "found zero column " << rc.m_j << " in row " << i << std::endl;
return false;
}
}
}
for (auto & c : m_columns) {
for (unsigned j = 0; j < m_columns.size(); j++) {
auto & c = m_columns[j];
std::unordered_set<unsigned> s;
for (auto & cc : c) {
if (s.find(cc.m_i) != s.end()) {
std::cout << "found row " << cc.m_i << " twice in a column" << std::endl;
std::cout << "found row " << cc.m_i << " twice in a column " << j << std::endl;
return false;
}
s.insert(cc.m_i);

31
src/util/lp/ul_pair.h
View file

@ -27,20 +27,23 @@ Revision History:
namespace lp {
enum lconstraint_kind {
LE = -2, LT = -1 , GE = 2, GT = 1, EQ = 0
};
inline std::ostream& operator<<(std::ostream& out, lconstraint_kind k) {
switch (k) {
case LE: return out << "<=";
case LT: return out << "<";
case GE: return out << ">=";
case GT: return out << ">";
case EQ: return out << "=";
}
return out << "??";
enum lconstraint_kind {
LE = -2, LT = -1 , GE = 2, GT = 1, EQ = 0
};
inline bool kind_is_strict(lconstraint_kind kind) { return kind == LT || kind == GT;}
inline std::ostream& operator<<(std::ostream& out, lconstraint_kind k) {
switch (k) {
case LE: return out << "<=";
case LT: return out << "<";
case GE: return out << ">=";
case GT: return out << ">";
case EQ: return out << "=";
}
return out << "??";
}
inline bool compare(const std::pair<mpq, var_index> & a, const std::pair<mpq, var_index> & b) {
return a.second < b.second;
@ -69,7 +72,7 @@ public:
m_low_bound_witness(static_cast<constraint_index>(-1)),
m_upper_bound_witness(static_cast<constraint_index>(-1)),
m_i(static_cast<row_index>(-1))
{}
{}
ul_pair(row_index ri) :
m_low_bound_witness(static_cast<constraint_index>(-1)),
m_upper_bound_witness(static_cast<constraint_index>(-1)),

1
src/util/rational.h
View file

@ -54,7 +54,6 @@ public:
rational(mpz const & z) { m().set(m_val, z); }
rational(double z) { UNREACHABLE(); }
explicit rational(char const * v) { m().set(m_val, v); }