diff --git a/src/smt/smt_context.cpp b/src/smt/smt_context.cpp index fd9b27069..926ab653a 100644 --- a/src/smt/smt_context.cpp +++ b/src/smt/smt_context.cpp @@ -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) diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index e8ab2b88e..1dccc11a4 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -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()); diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 911c127a1..f95dba98a 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -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 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 void test_matrix(sparse_matrix & a) { auto m = a.dimension(); -// copy a to b in the reversed order + // copy a to b in the reversed order sparse_matrix 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 & basis_heading, vector(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(non_basic_columns.size()); } } } @@ -554,7 +557,7 @@ void test_lp_0() { vector nbasis; vector heading; - lp_primal_core_solver lpsolver(m_, b, x_star, basis, nbasis, heading, costs, column_types, upper_bound_values, settings, cn); + lp_primal_core_solver 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& 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& 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 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 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 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 & 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 * 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> 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> 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> 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> 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& m_file_sizes; - file_size_comp(unordered_map& 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& m_file_sizes; + file_size_comp(unordered_map& 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> 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> 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 non_basic_columns; lu lsuhl(A, basis, settings); - indexed_vector d(A.row_count()); + indexed_vector 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 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 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 & 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 & 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 & 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 d; + vector> 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 cs([](unsigned i) + { + if (i == 0) return std::string("x"); + if (i == 1) return std::string("y"); + return std::to_string(i); + }); + vector> 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); } diff --git a/src/util/lp/CMakeLists.txt b/src/util/lp/CMakeLists.txt index efe683b42..681c18463 100644 --- a/src/util/lp/CMakeLists.txt +++ b/src/util/lp/CMakeLists.txt @@ -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 diff --git a/src/util/lp/bound_propagator.cpp b/src/util/lp/bound_propagator.cpp index aac57af48..8b106c521 100644 --- a/src/util/lp/bound_propagator.cpp +++ b/src/util/lp/bound_propagator.cpp @@ -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); diff --git a/src/util/lp/cut_solver.cpp b/src/util/lp/cut_solver.cpp new file mode 100644 index 000000000..2b527aa01 --- /dev/null +++ b/src/util/lp/cut_solver.cpp @@ -0,0 +1,8 @@ +/* + Copyright (c) 2017 Microsoft Corporation + Author: Nikolaj Bjorner, Lev Nachmanson +*/ +#include "util/lp/cut_solver.h" +namespace lp { +} + diff --git a/src/util/lp/cut_solver.h b/src/util/lp/cut_solver.h new file mode 100644 index 000000000..6fcf041b6 --- /dev/null +++ b/src/util/lp/cut_solver.h @@ -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 +namespace lp { +template +class cut_solver : public column_namer { +public: // for debugging + struct polynomial { + // the polynom evaluates to m_term + m_a + vector> m_term; + mpq m_a; // the free coefficient + polynomial(vector>& p, const mpq & a) : m_term(p), m_a(a) { + } + polynomial(vector>& p) : polynomial(p, 0) { + } + + }; + + struct ineq { // we only have less or equal, which is enough for integral variables + polynomial m_poly; + ineq(vector>& term, const mpq& a): m_poly(term, a) { + } + }; + + vector 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 m_literals; // point to m_trail + integer_domain m_domain; + }; + + vector m_var_infos; + + bool lhs_is_int(const vector> & lhs) const { + for (auto & p : lhs) { + if (numeric_traits::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> & lhs, const mpq& free_coeff) { + lp_assert(lhs_is_int(lhs)); + lp_assert(free_coeff.is_int()); + vector> 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 m_var_name_function; + bool m_inconsistent; // tracks if state is consistent + unsigned m_scope_lvl; // tracks the number of case splits + + svector 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 m_scopes; + std::unordered_map 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 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 & 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; + } + } + } +}; +} diff --git a/src/util/lp/disjoint_intervals.h b/src/util/lp/disjoint_intervals.h deleted file mode 100644 index 5f4f31af6..000000000 --- a/src/util/lp/disjoint_intervals.h +++ /dev/null @@ -1,334 +0,0 @@ -/* - Copyright (c) 2017 Microsoft Corporation - Author: Lev Nachmanson -*/ -#pragma once -#include -namespace lp { -// represents the set of disjoint intervals of integer number -struct disjoint_intervals { - std::map 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::iterator & it) const { - return is_start(it->second); - } - bool is_start(const std::map::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::iterator & it) const { - return is_end(it->second); - } - bool is_end(const std::map::reverse_iterator & it) const { - return is_end(it->second); - } - - int pos(const std::map::iterator & it) const { - return it->first; - } - int pos(const std::map::reverse_iterator & it) const { - return it->first; - } - - int bound_kind(const std::map::iterator & it) const { - return it->second; - } - - int bound_kind(const std::map::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::iterator & it) const { - return is_proper_end(it->second); - } - bool is_proper_end(const std::map::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::iterator & it) const { - return is_one_point_interval(it->second); - } - bool is_one_point_interval(const std::map::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"; - } - - -}; -} diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index d717c4047..e1beb73ec 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -5,24 +5,10 @@ #include "util/lp/int_solver.h" #include "util/lp/lar_solver.h" +#include "util/lp/cut_solver.h" +#include 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& iter) { @@ -394,117 +380,166 @@ int int_solver::find_free_var_in_gomory_row(linear_combination_iterator& it return -1; } -lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, linear_combination_iterator& 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* 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(); - 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 +void int_solver::fill_cut_solver(cut_solver & cs) { + for (lar_base_constraint * c : m_lar_solver->constraints()) + fill_cut_solver_for_constraint(c, cs); +} + +template +void int_solver::fill_cut_solver_for_constraint(const lar_base_constraint* c, cut_solver & cs) { + vector> coeffs; + T rs; + get_int_coeffs_from_constraint(c, coeffs, rs); + cs.add_ineq(coeffs, rs); +} + + +// it produces an inequality coeff*x <= rs +template +void int_solver::get_int_coeffs_from_constraint(const lar_base_constraint* c, vector>& 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> 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::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 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 &it) { @@ -595,55 +645,55 @@ bool int_solver::gcd_test_for_row(static_matrix> & 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 & 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 & it, @@ -735,7 +785,7 @@ bool int_solver::ext_gcd_test(iterator_on_row & it, linear_combination_iterator * int_solver::get_column_iterator(unsigned j) { if (m_lar_solver->use_tableau()) return new iterator_on_column(m_lar_solver->A_r().m_columns[j], m_lar_solver->A_r()); - return new iterator_on_indexed_vector(m_lar_solver->get_column_in_lu_mode(j)); + return new iterator_on_indexed_vector(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); + } +} } diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index 12723effc..640ed9b3e 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -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 @@ -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 & 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 * 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& 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 * r); - lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j, linear_combination_iterator& 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& iter); bool is_gomory_cut_target(linear_combination_iterator &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 + void fill_cut_solver(cut_solver & cs); + template + void fill_cut_solver_for_constraint(const lar_base_constraint*, cut_solver& ); + template + void get_int_coeffs_from_constraint(const lar_base_constraint* c, vector>& coeff, T & rs); + }; } diff --git a/src/util/lp/integer_domain.h b/src/util/lp/integer_domain.h new file mode 100644 index 000000000..961a34da3 --- /dev/null +++ b/src/util/lp/integer_domain.h @@ -0,0 +1,876 @@ +/* + Copyright (c) 2017 Microsoft Corporation + Author: Lev Nachmanson +*/ +#pragma once +#include +#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 +class integer_domain { +#ifdef Z3DEBUG + std::set m_domain; +#endif + typedef typename std::map::iterator iter; + typedef typename std::map::const_iterator const_iter; + typedef typename std::map::reverse_iterator riter; + stacked_map m_endpoints; // 0 means start, 1 means end, 2 means both - for a point interval + stacked_value 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) : +#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) { + 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; + } +}; +} diff --git a/src/util/lp/lar_constraints.h b/src/util/lp/lar_constraints.h index 280b24f20..ac15028bb 100644 --- a/src/util/lp/lar_constraints.h +++ b/src/util/lp/lar_constraints.h @@ -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> get_left_side_coefficients() const = 0; diff --git a/src/util/lp/lar_core_solver.h b/src/util/lp/lar_core_solver.h index 3daf2a3c1..136db9119 100644 --- a/src/util/lp/lar_core_solver.h +++ b/src/util/lp/lar_core_solver.h @@ -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; } diff --git a/src/util/lp/lar_core_solver.hpp b/src/util/lp/lar_core_solver.hpp index 19833f5eb..0943d6977 100644 --- a/src/util/lp/lar_core_solver.hpp +++ b/src/util/lp/lar_core_solver.hpp @@ -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(); diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 191d1ab6f..70c322c4a 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -148,8 +148,8 @@ bool lar_solver::implied_bound_is_correctly_explained(implied_bound const & be, void lar_solver::analyze_new_bounds_on_row( - unsigned row_index, - bound_propagator & bp) { + unsigned row_index, + bound_propagator & bp) { lp_assert(!use_tableau()); iterator_on_pivot_row it(m_mpq_lar_core_solver.get_pivot_row(), m_mpq_lar_core_solver.m_r_basis[row_index]); @@ -162,8 +162,8 @@ void lar_solver::analyze_new_bounds_on_row( } void lar_solver::analyze_new_bounds_on_row_tableau( - unsigned row_index, - bound_propagator & bp + unsigned row_index, + bound_propagator & bp ) { if (A_r().m_rows[row_index].size() > settings().max_row_length_for_bound_propagation) @@ -285,7 +285,6 @@ bool lar_solver::has_int_var() const { } lp_status lar_solver::find_feasible_solution() { - TRACE("arith_int", ); lp_assert(inf_int_set_is_correct()); m_settings.st().m_make_feasible++; if (A_r().column_count() > m_settings.st().m_max_cols) @@ -301,6 +300,7 @@ lp_status lar_solver::find_feasible_solution() { m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = true; auto ret = solve(); + TRACE("intinf", m_int_solver->display_inf_or_int_inf_columns(tout);); lp_assert(inf_int_set_is_correct()); return ret; } @@ -379,10 +379,7 @@ void lar_solver::pop(unsigned k) { if (m_settings.use_tableau()) { pop_tableau(); } - TRACE("arith_int", tout << "pop" << std::endl;); - lp_assert(inf_int_set_is_correct()); - TRACE("arith_int", tout << "pop" << std::endl;); - + lp_assert(A_r().column_count() == n); m_columns_to_ul_pairs.pop(k); m_mpq_lar_core_solver.pop(k); @@ -395,7 +392,7 @@ void lar_solver::pop(unsigned k) { clean_inf_set_of_r_solver_after_pop(); m_status = m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()? lp_status::OPTIMAL: lp_status::UNKNOWN; lp_assert(m_settings.simplex_strategy() == simplex_strategy_enum::undecided || - (!use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); + (!use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()); lp_assert(ax_is_correct()); @@ -755,6 +752,7 @@ void lar_solver::solve_with_core_solver() { update_x_and_inf_costs_for_columns_with_changed_bounds_tableau(); else update_x_and_inf_costs_for_columns_with_changed_bounds(); + TRACE("intinf", m_int_solver->display_inf_or_int_inf_columns(tout);); m_mpq_lar_core_solver.solve(); set_status(m_mpq_lar_core_solver.m_r_solver.get_status()); lp_assert(inf_int_set_is_correct()); @@ -1121,9 +1119,9 @@ void lar_solver::get_infeasibility_explanation(vector> & explanation, - const vector> & inf_row, - int inf_sign) const { + vector> & explanation, + const vector> & inf_row, + int inf_sign) const { for (auto & it : inf_row) { mpq coeff = it.first; @@ -1143,7 +1141,6 @@ void lar_solver::get_model(std::unordered_map & variable_values) lp_assert(m_status == lp_status::OPTIMAL); unsigned i; do { - // different pairs have to produce different singleton values std::unordered_set set_of_different_pairs; std::unordered_set set_of_different_singles; @@ -1316,10 +1313,9 @@ void lar_solver::remove_last_row_and_column_from_tableau(unsigned j) { slv.m_b.pop_back(); } -void lar_solver::remove_last_column_from_tableau(unsigned j) { - lp_assert(j == A_r().column_count() - 1); +void lar_solver::remove_last_column_from_A() { // the last column has to be empty - lp_assert(A_r().m_columns[j].size() == 0); + lp_assert(A_r().m_columns.back().size() == 0); A_r().m_columns.pop_back(); } @@ -1351,9 +1347,9 @@ void lar_solver::remove_last_column_from_basis_tableau(unsigned j) { lp_assert(rslv.basis_heading_is_correct()); } -void lar_solver::remove_column_from_tableau(unsigned j) { +void lar_solver::remove_last_column_from_tableau() { auto& rslv = m_mpq_lar_core_solver.m_r_solver; - lp_assert(j == A_r().column_count() - 1); + unsigned j = A_r().column_count() - 1; lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); if (column_represents_row_in_tableau(j)) { remove_last_row_and_column_from_tableau(j); @@ -1361,7 +1357,7 @@ void lar_solver::remove_column_from_tableau(unsigned j) { rslv.change_basis_unconditionally(j, rslv.m_basis[A_r().row_count()]); // A_r().row_count() is the index of the last row in the basis still } else { - remove_last_column_from_tableau(j); + remove_last_column_from_A(); } rslv.m_x.pop_back(); rslv.m_d.pop_back(); @@ -1379,8 +1375,8 @@ void lar_solver::pop_tableau() { lp_assert(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct()); // We remove last variables starting from m_column_names.size() to m_vec_of_canonic_left_sides.size(). // At this moment m_column_names is already popped - for (unsigned j = A_r().column_count(); j-- > m_columns_to_ext_vars_or_term_indices.size();) - remove_column_from_tableau(j); + while (A_r().column_count() > m_columns_to_ext_vars_or_term_indices.size()) + remove_last_column_from_tableau(); lp_assert(m_mpq_lar_core_solver.m_r_solver.m_costs.size() == A_r().column_count()); lp_assert(m_mpq_lar_core_solver.m_r_solver.m_basis.size() == A_r().row_count()); lp_assert(m_mpq_lar_core_solver.m_r_solver.basis_heading_is_correct()); @@ -1439,10 +1435,10 @@ void lar_solver::shrink_explanation_to_minimum(vectorsecond.is_integer(); } - // below is the initialization functionality of lar_solver +// below is the initialization functionality of lar_solver bool lar_solver::strategy_is_undecided() const { - return m_settings.simplex_strategy() == simplex_strategy_enum::undecided; + return m_settings.simplex_strategy() == simplex_strategy_enum::undecided; } var_index lar_solver::add_var(unsigned ext_j, bool is_int) { TRACE("add_var", tout << "adding var " << ext_j << (is_int? " int" : " nonint") << std::endl;); - var_index i; - lp_assert(ext_j < m_terms_start_index); + var_index i; + lp_assert(ext_j < m_terms_start_index); - if (ext_j >= m_terms_start_index) - throw 0; // todo : what is the right way to exit? - auto it = m_ext_vars_to_columns.find(ext_j); - if (it != m_ext_vars_to_columns.end()) { - return it->second.ext_j(); - } - lp_assert(m_columns_to_ul_pairs.size() == A_r().column_count()); - i = A_r().column_count(); - m_columns_to_ul_pairs.push_back(ul_pair(static_cast(-1))); - add_non_basic_var_to_core_fields(ext_j, is_int); - lp_assert(sizes_are_correct()); + if (ext_j >= m_terms_start_index) + throw 0; // todo : what is the right way to exit? + auto it = m_ext_vars_to_columns.find(ext_j); + if (it != m_ext_vars_to_columns.end()) { + return it->second.ext_j(); + } + lp_assert(m_columns_to_ul_pairs.size() == A_r().column_count()); + i = A_r().column_count(); + m_columns_to_ul_pairs.push_back(ul_pair(static_cast(-1))); + add_non_basic_var_to_core_fields(ext_j, is_int); + lp_assert(sizes_are_correct()); if (is_int) { m_mpq_lar_core_solver.m_r_solver.set_tracker_of_x(& m_tracker_of_x_change); } lp_assert(inf_int_set_is_correct()); - return i; + return i; } void lar_solver::register_new_ext_var_index(unsigned ext_v, bool is_int) { - lp_assert(!contains(m_ext_vars_to_columns, ext_v)); - unsigned j = static_cast(m_ext_vars_to_columns.size()); - m_ext_vars_to_columns.insert(std::make_pair(ext_v, ext_var_info(j, is_int))); - lp_assert(m_columns_to_ext_vars_or_term_indices.size() == j); - m_columns_to_ext_vars_or_term_indices.push_back(ext_v); + lp_assert(!contains(m_ext_vars_to_columns, ext_v)); + unsigned j = static_cast(m_ext_vars_to_columns.size()); + m_ext_vars_to_columns.insert(std::make_pair(ext_v, ext_var_info(j, is_int))); + lp_assert(m_columns_to_ext_vars_or_term_indices.size() == j); + m_columns_to_ext_vars_or_term_indices.push_back(ext_v); } void lar_solver::add_non_basic_var_to_core_fields(unsigned ext_j, bool is_int) { - register_new_ext_var_index(ext_j, is_int); - m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column); - m_columns_with_changed_bound.increase_size_by_one(); + register_new_ext_var_index(ext_j, is_int); + m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column); + m_columns_with_changed_bound.increase_size_by_one(); m_inf_int_set.increase_size_by_one(); - add_new_var_to_core_fields_for_mpq(false); - if (use_lu()) - add_new_var_to_core_fields_for_doubles(false); + add_new_var_to_core_fields_for_mpq(false); + if (use_lu()) + add_new_var_to_core_fields_for_doubles(false); } void lar_solver::add_new_var_to_core_fields_for_doubles(bool register_in_basis) { - unsigned j = A_d().column_count(); - A_d().add_column(); - lp_assert(m_mpq_lar_core_solver.m_d_x.size() == j); - // lp_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j); // restore later - m_mpq_lar_core_solver.m_d_x.resize(j + 1); - m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1); - m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1); - lp_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method - if (register_in_basis) { - A_d().add_row(); - m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size()); - m_mpq_lar_core_solver.m_d_basis.push_back(j); - } - else { - m_mpq_lar_core_solver.m_d_heading.push_back(-static_cast(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1); - m_mpq_lar_core_solver.m_d_nbasis.push_back(j); - } + unsigned j = A_d().column_count(); + A_d().add_column(); + lp_assert(m_mpq_lar_core_solver.m_d_x.size() == j); + // lp_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j); // restore later + m_mpq_lar_core_solver.m_d_x.resize(j + 1); + m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1); + m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1); + lp_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method + if (register_in_basis) { + A_d().add_row(); + m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size()); + m_mpq_lar_core_solver.m_d_basis.push_back(j); + } + else { + m_mpq_lar_core_solver.m_d_heading.push_back(-static_cast(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1); + m_mpq_lar_core_solver.m_d_nbasis.push_back(j); + } } void lar_solver::add_new_var_to_core_fields_for_mpq(bool register_in_basis) { - unsigned j = A_r().column_count(); - A_r().add_column(); - lp_assert(m_mpq_lar_core_solver.m_r_x.size() == j); - // lp_assert(m_mpq_lar_core_solver.m_r_low_bounds.size() == j && m_mpq_lar_core_solver.m_r_upper_bounds.size() == j); // restore later - m_mpq_lar_core_solver.m_r_x.resize(j + 1); - m_mpq_lar_core_solver.m_r_low_bounds.increase_size_by_one(); - m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one(); - m_mpq_lar_core_solver.m_r_solver.m_inf_set.increase_size_by_one(); - m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1); - m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1); - lp_assert(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method - if (register_in_basis) { - A_r().add_row(); - m_mpq_lar_core_solver.m_r_heading.push_back(m_mpq_lar_core_solver.m_r_basis.size()); - m_mpq_lar_core_solver.m_r_basis.push_back(j); - if (m_settings.bound_propagation()) - m_rows_with_changed_bounds.insert(A_r().row_count() - 1); - } - else { - m_mpq_lar_core_solver.m_r_heading.push_back(-static_cast(m_mpq_lar_core_solver.m_r_nbasis.size()) - 1); - m_mpq_lar_core_solver.m_r_nbasis.push_back(j); - } + unsigned j = A_r().column_count(); + A_r().add_column(); + lp_assert(m_mpq_lar_core_solver.m_r_x.size() == j); + // lp_assert(m_mpq_lar_core_solver.m_r_low_bounds.size() == j && m_mpq_lar_core_solver.m_r_upper_bounds.size() == j); // restore later + m_mpq_lar_core_solver.m_r_x.resize(j + 1); + m_mpq_lar_core_solver.m_r_low_bounds.increase_size_by_one(); + m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one(); + m_mpq_lar_core_solver.m_r_solver.m_inf_set.increase_size_by_one(); + m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1); + m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1); + lp_assert(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method + if (register_in_basis) { + A_r().add_row(); + m_mpq_lar_core_solver.m_r_heading.push_back(m_mpq_lar_core_solver.m_r_basis.size()); + m_mpq_lar_core_solver.m_r_basis.push_back(j); + if (m_settings.bound_propagation()) + m_rows_with_changed_bounds.insert(A_r().row_count() - 1); + } + else { + m_mpq_lar_core_solver.m_r_heading.push_back(-static_cast(m_mpq_lar_core_solver.m_r_nbasis.size()) - 1); + m_mpq_lar_core_solver.m_r_nbasis.push_back(j); + } } var_index lar_solver::add_term_undecided(const vector> & coeffs, - const mpq &m_v) { - m_terms.push_back(new lar_term(coeffs, m_v)); - return m_terms_start_index + m_terms.size() - 1; + const mpq &m_v) { + m_terms.push_back(new lar_term(coeffs, m_v)); + return m_terms_start_index + m_terms.size() - 1; } // terms var_index lar_solver::add_term(const vector> & coeffs, - const mpq &m_v) { - if (strategy_is_undecided()) - return add_term_undecided(coeffs, m_v); + const mpq &m_v) { + if (strategy_is_undecided()) + return add_term_undecided(coeffs, m_v); - m_terms.push_back(new lar_term(coeffs, m_v)); - unsigned adjusted_term_index = m_terms.size() - 1; - var_index ret = m_terms_start_index + adjusted_term_index; - if (use_tableau() && !coeffs.empty()) { - add_row_from_term_no_constraint(m_terms.back(), ret); - if (m_settings.bound_propagation()) - m_rows_with_changed_bounds.insert(A_r().row_count() - 1); - } - lp_assert(m_ext_vars_to_columns.size() == A_r().column_count()); - return ret; + m_terms.push_back(new lar_term(coeffs, m_v)); + unsigned adjusted_term_index = m_terms.size() - 1; + var_index ret = m_terms_start_index + adjusted_term_index; + if (use_tableau() && !coeffs.empty()) { + add_row_from_term_no_constraint(m_terms.back(), ret); + if (m_settings.bound_propagation()) + m_rows_with_changed_bounds.insert(A_r().row_count() - 1); + } + lp_assert(m_ext_vars_to_columns.size() == A_r().column_count()); + return ret; } void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) { - register_new_ext_var_index(term_ext_index, term_is_int(term)); - // j will be a new variable - unsigned j = A_r().column_count(); - ul_pair ul(j); - m_columns_to_ul_pairs.push_back(ul); - add_basic_var_to_core_fields(); - if (use_tableau()) { - auto it = iterator_on_term_with_basis_var(*term, j); - A_r().fill_last_row_with_pivoting(it, - m_mpq_lar_core_solver.m_r_solver.m_basis_heading); - m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type()); - } - else { - fill_last_row_of_A_r(A_r(), term); - } + register_new_ext_var_index(term_ext_index, term_is_int(term)); + // j will be a new variable + unsigned j = A_r().column_count(); + ul_pair ul(j); + m_columns_to_ul_pairs.push_back(ul); + add_basic_var_to_core_fields(); + if (use_tableau()) { + auto it = iterator_on_term_with_basis_var(*term, j); + A_r().fill_last_row_with_pivoting(it, + m_mpq_lar_core_solver.m_r_solver.m_basis_heading); + m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type()); + } + else { + fill_last_row_of_A_r(A_r(), term); + } m_mpq_lar_core_solver.m_r_solver.update_x_and_call_tracker(j, get_basic_var_value_from_row_directly(A_r().row_count() - 1)); - if (use_lu()) - fill_last_row_of_A_d(A_d(), term); + if (use_lu()) + fill_last_row_of_A_d(A_d(), term); lp_assert(inf_int_set_is_correct()); } void lar_solver::add_basic_var_to_core_fields() { - bool use_lu = m_mpq_lar_core_solver.need_to_presolve_with_double_solver(); - lp_assert(!use_lu || A_r().column_count() == A_d().column_count()); - m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column); - m_columns_with_changed_bound.increase_size_by_one(); - m_rows_with_changed_bounds.increase_size_by_one(); + bool use_lu = m_mpq_lar_core_solver.need_to_presolve_with_double_solver(); + lp_assert(!use_lu || A_r().column_count() == A_d().column_count()); + m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column); + m_columns_with_changed_bound.increase_size_by_one(); + m_rows_with_changed_bounds.increase_size_by_one(); m_inf_int_set.increase_size_by_one(); - add_new_var_to_core_fields_for_mpq(true); - if (use_lu) - add_new_var_to_core_fields_for_doubles(true); + add_new_var_to_core_fields_for_mpq(true); + if (use_lu) + add_new_var_to_core_fields_for_doubles(true); } bool lar_solver::bound_is_integer_if_needed(unsigned j, const mpq & right_side) const { @@ -1638,450 +1634,450 @@ bool lar_solver::bound_is_integer_if_needed(unsigned j, const mpq & right_side) constraint_index lar_solver::add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side) { TRACE("lar_solver", tout << "j = " << j << std::endl;); - constraint_index ci = m_constraints.size(); - if (!is_term(j)) { // j is a var + constraint_index ci = m_constraints.size(); + if (!is_term(j)) { // j is a var lp_assert(bound_is_integer_if_needed(j, right_side)); - auto vc = new lar_var_constraint(j, kind, right_side); - m_constraints.push_back(vc); - update_column_type_and_bound(j, kind, right_side, ci); - } - else { - add_var_bound_on_constraint_for_term(j, kind, right_side, ci); - } - lp_assert(sizes_are_correct()); + auto vc = new lar_var_constraint(j, kind, right_side); + m_constraints.push_back(vc); + update_column_type_and_bound(j, kind, right_side, ci); + } + else { + add_var_bound_on_constraint_for_term(j, kind, right_side, ci); + } + lp_assert(sizes_are_correct()); lp_assert(inf_int_set_is_correct()); - return ci; + return ci; } void lar_solver::update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) { - switch (m_mpq_lar_core_solver.m_column_types[j]) { - case column_type::free_column: - update_free_column_type_and_bound(j, kind, right_side, constr_index); - break; - case column_type::boxed: - update_boxed_column_type_and_bound(j, kind, right_side, constr_index); - break; - case column_type::low_bound: - update_low_bound_column_type_and_bound(j, kind, right_side, constr_index); - break; - case column_type::upper_bound: - update_upper_bound_column_type_and_bound(j, kind, right_side, constr_index); - break; - case column_type::fixed: - update_fixed_column_type_and_bound(j, kind, right_side, constr_index); - break; - default: - lp_assert(false); // cannot be here - } + switch (m_mpq_lar_core_solver.m_column_types[j]) { + case column_type::free_column: + update_free_column_type_and_bound(j, kind, right_side, constr_index); + break; + case column_type::boxed: + update_boxed_column_type_and_bound(j, kind, right_side, constr_index); + break; + case column_type::low_bound: + update_low_bound_column_type_and_bound(j, kind, right_side, constr_index); + break; + case column_type::upper_bound: + update_upper_bound_column_type_and_bound(j, kind, right_side, constr_index); + break; + case column_type::fixed: + update_fixed_column_type_and_bound(j, kind, right_side, constr_index); + break; + default: + lp_assert(false); // cannot be here + } } void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) { - lp_assert(is_term(j)); - unsigned adjusted_term_index = adjust_term_index(j); + lp_assert(is_term(j)); + unsigned adjusted_term_index = adjust_term_index(j); lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int()); - auto it = m_ext_vars_to_columns.find(j); - if (it != m_ext_vars_to_columns.end()) { - unsigned term_j = it->second.ext_j(); - mpq rs = right_side - m_terms[adjusted_term_index]->m_v; - m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side)); - update_column_type_and_bound(term_j, kind, rs, ci); - } - else { - add_constraint_from_term_and_create_new_column_row(j, m_terms[adjusted_term_index], kind, right_side); - } + auto it = m_ext_vars_to_columns.find(j); + if (it != m_ext_vars_to_columns.end()) { + unsigned term_j = it->second.ext_j(); + mpq rs = right_side - m_terms[adjusted_term_index]->m_v; + m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side)); + update_column_type_and_bound(term_j, kind, rs, ci); + } + else { + add_constraint_from_term_and_create_new_column_row(j, m_terms[adjusted_term_index], kind, right_side); + } } constraint_index lar_solver::add_constraint(const vector>& left_side_with_terms, lconstraint_kind kind_par, const mpq& right_side_parm) { - vector> left_side; - mpq rs = -right_side_parm; - substitute_terms_in_linear_expression(left_side_with_terms, left_side, rs); - unsigned term_index = add_term(left_side, zero_of_type()); - constraint_index ci = m_constraints.size(); - add_var_bound_on_constraint_for_term(term_index, kind_par, -rs, ci); - return ci; + vector> left_side; + mpq rs = -right_side_parm; + substitute_terms_in_linear_expression(left_side_with_terms, left_side, rs); + unsigned term_index = add_term(left_side, zero_of_type()); + constraint_index ci = m_constraints.size(); + add_var_bound_on_constraint_for_term(term_index, kind_par, -rs, ci); + return ci; } void lar_solver::add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term, - lconstraint_kind kind, const mpq & right_side) { + lconstraint_kind kind, const mpq & right_side) { - add_row_from_term_no_constraint(term, term_j); - unsigned j = A_r().column_count() - 1; - update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size()); - m_constraints.push_back(new lar_term_constraint(term, kind, right_side)); - lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); + add_row_from_term_no_constraint(term, term_j); + unsigned j = A_r().column_count() - 1; + update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size()); + m_constraints.push_back(new lar_term_constraint(term, kind, right_side)); + lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size()); } void lar_solver::decide_on_strategy_and_adjust_initial_state() { - lp_assert(strategy_is_undecided()); - if (m_columns_to_ul_pairs.size() > m_settings.column_number_threshold_for_using_lu_in_lar_solver) { - m_settings.simplex_strategy() = simplex_strategy_enum::lu; - } - else { - m_settings.simplex_strategy() = simplex_strategy_enum::tableau_rows; // todo: when to switch to tableau_costs? - } - adjust_initial_state(); + lp_assert(strategy_is_undecided()); + if (m_columns_to_ul_pairs.size() > m_settings.column_number_threshold_for_using_lu_in_lar_solver) { + m_settings.simplex_strategy() = simplex_strategy_enum::lu; + } + else { + m_settings.simplex_strategy() = simplex_strategy_enum::tableau_rows; // todo: when to switch to tableau_costs? + } + adjust_initial_state(); } void lar_solver::adjust_initial_state() { - switch (m_settings.simplex_strategy()) { - case simplex_strategy_enum::lu: - adjust_initial_state_for_lu(); - break; - case simplex_strategy_enum::tableau_rows: - adjust_initial_state_for_tableau_rows(); - break; - case simplex_strategy_enum::tableau_costs: - lp_assert(false); // not implemented - case simplex_strategy_enum::undecided: - adjust_initial_state_for_tableau_rows(); - break; - } + switch (m_settings.simplex_strategy()) { + case simplex_strategy_enum::lu: + adjust_initial_state_for_lu(); + break; + case simplex_strategy_enum::tableau_rows: + adjust_initial_state_for_tableau_rows(); + break; + case simplex_strategy_enum::tableau_costs: + lp_assert(false); // not implemented + case simplex_strategy_enum::undecided: + adjust_initial_state_for_tableau_rows(); + break; + } } void lar_solver::adjust_initial_state_for_lu() { - copy_from_mpq_matrix(A_d()); - unsigned n = A_d().column_count(); - m_mpq_lar_core_solver.m_d_x.resize(n); - m_mpq_lar_core_solver.m_d_low_bounds.resize(n); - m_mpq_lar_core_solver.m_d_upper_bounds.resize(n); - m_mpq_lar_core_solver.m_d_heading = m_mpq_lar_core_solver.m_r_heading; - m_mpq_lar_core_solver.m_d_basis = m_mpq_lar_core_solver.m_r_basis; + copy_from_mpq_matrix(A_d()); + unsigned n = A_d().column_count(); + m_mpq_lar_core_solver.m_d_x.resize(n); + m_mpq_lar_core_solver.m_d_low_bounds.resize(n); + m_mpq_lar_core_solver.m_d_upper_bounds.resize(n); + m_mpq_lar_core_solver.m_d_heading = m_mpq_lar_core_solver.m_r_heading; + m_mpq_lar_core_solver.m_d_basis = m_mpq_lar_core_solver.m_r_basis; - /* - unsigned j = A_d().column_count(); - A_d().add_column(); - lp_assert(m_mpq_lar_core_solver.m_d_x.size() == j); - // lp_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j); // restore later - m_mpq_lar_core_solver.m_d_x.resize(j + 1 ); - m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1); - m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1); - lp_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method - if (register_in_basis) { - A_d().add_row(); - m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size()); - m_mpq_lar_core_solver.m_d_basis.push_back(j); - }else { - m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1); - m_mpq_lar_core_solver.m_d_nbasis.push_back(j); - }*/ + /* + unsigned j = A_d().column_count(); + A_d().add_column(); + lp_assert(m_mpq_lar_core_solver.m_d_x.size() == j); + // lp_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j); // restore later + m_mpq_lar_core_solver.m_d_x.resize(j + 1 ); + m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1); + m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1); + lp_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method + if (register_in_basis) { + A_d().add_row(); + m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size()); + m_mpq_lar_core_solver.m_d_basis.push_back(j); + }else { + m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1); + m_mpq_lar_core_solver.m_d_nbasis.push_back(j); + }*/ } void lar_solver::adjust_initial_state_for_tableau_rows() { - for (unsigned j = 0; j < m_terms.size(); j++) { - if (contains(m_ext_vars_to_columns, j + m_terms_start_index)) - continue; - add_row_from_term_no_constraint(m_terms[j], j + m_terms_start_index); - } + for (unsigned j = 0; j < m_terms.size(); j++) { + if (contains(m_ext_vars_to_columns, j + m_terms_start_index)) + continue; + add_row_from_term_no_constraint(m_terms[j], j + m_terms_start_index); + } } // this fills the last row of A_d and sets the basis column: -1 in the last column of the row void lar_solver::fill_last_row_of_A_d(static_matrix & A, const lar_term* ls) { - lp_assert(A.row_count() > 0); - lp_assert(A.column_count() > 0); - unsigned last_row = A.row_count() - 1; - lp_assert(A.m_rows[last_row].empty()); + lp_assert(A.row_count() > 0); + lp_assert(A.column_count() > 0); + unsigned last_row = A.row_count() - 1; + lp_assert(A.m_rows[last_row].empty()); - for (auto & t : ls->m_coeffs) { - lp_assert(!is_zero(t.second)); - var_index j = t.first; - A.set(last_row, j, -t.second.get_double()); - } + for (auto & t : ls->m_coeffs) { + lp_assert(!is_zero(t.second)); + var_index j = t.first; + A.set(last_row, j, -t.second.get_double()); + } - unsigned basis_j = A.column_count() - 1; - A.set(last_row, basis_j, -1); + unsigned basis_j = A.column_count() - 1; + A.set(last_row, basis_j, -1); } void lar_solver::update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) { - mpq y_of_bound(0); - switch (kind) { - case LT: - y_of_bound = -1; - case LE: - m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound; - lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound); - lp_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j); - { - auto up = numeric_pair(right_side, y_of_bound); - m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - } - set_upper_bound_witness(j, constr_ind); - break; - case GT: - y_of_bound = 1; - case GE: - m_mpq_lar_core_solver.m_column_types[j] = column_type::low_bound; - lp_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j); - { - auto low = numeric_pair(right_side, y_of_bound); - m_mpq_lar_core_solver.m_r_low_bounds[j] = low; - } - set_low_bound_witness(j, constr_ind); - break; - case EQ: - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = numeric_pair(right_side, zero_of_type()); - set_upper_bound_witness(j, constr_ind); - set_low_bound_witness(j, constr_ind); - break; + mpq y_of_bound(0); + switch (kind) { + case LT: + y_of_bound = -1; + case LE: + m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound; + lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound); + lp_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j); + { + auto up = numeric_pair(right_side, y_of_bound); + m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; + } + set_upper_bound_witness(j, constr_ind); + break; + case GT: + y_of_bound = 1; + case GE: + m_mpq_lar_core_solver.m_column_types[j] = column_type::low_bound; + lp_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j); + { + auto low = numeric_pair(right_side, y_of_bound); + m_mpq_lar_core_solver.m_r_low_bounds[j] = low; + } + set_low_bound_witness(j, constr_ind); + break; + case EQ: + m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; + m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = numeric_pair(right_side, zero_of_type()); + set_upper_bound_witness(j, constr_ind); + set_low_bound_witness(j, constr_ind); + break; - default: - lp_unreachable(); + default: + lp_unreachable(); - } - m_columns_with_changed_bound.insert(j); + } + m_columns_with_changed_bound.insert(j); } void lar_solver::update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) { - lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound); - mpq y_of_bound(0); - switch (kind) { - case LT: - y_of_bound = -1; - case LE: + lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound); + mpq y_of_bound(0); + switch (kind) { + case LT: + y_of_bound = -1; + case LE: { - auto up = numeric_pair(right_side, y_of_bound); - if (up < m_mpq_lar_core_solver.m_r_upper_bounds()[j]) { - m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - set_upper_bound_witness(j, ci); - m_columns_with_changed_bound.insert(j); - } + auto up = numeric_pair(right_side, y_of_bound); + if (up < m_mpq_lar_core_solver.m_r_upper_bounds()[j]) { + m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; + set_upper_bound_witness(j, ci); + m_columns_with_changed_bound.insert(j); + } } break; - case GT: - y_of_bound = 1; - case GE: - m_mpq_lar_core_solver.m_column_types[j] = column_type::boxed; - { - auto low = numeric_pair(right_side, y_of_bound); - m_mpq_lar_core_solver.m_r_low_bounds[j] = low; - set_low_bound_witness(j, ci); - m_columns_with_changed_bound.insert(j); - if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { + case GT: + y_of_bound = 1; + case GE: + m_mpq_lar_core_solver.m_column_types[j] = column_type::boxed; + { + auto low = numeric_pair(right_side, y_of_bound); + m_mpq_lar_core_solver.m_r_low_bounds[j] = low; + set_low_bound_witness(j, ci); + m_columns_with_changed_bound.insert(j); + if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - } - else { - m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j] ? column_type::boxed : column_type::fixed; - } - } - break; - case EQ: + m_infeasible_column_index = j; + } + else { + m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j] ? column_type::boxed : column_type::fixed; + } + } + break; + case EQ: { - auto v = numeric_pair(right_side, zero_of_type()); - if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_status = lp_status::INFEASIBLE; - set_low_bound_witness(j, ci); - m_infeasible_column_index = j; - } - else { - m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v; - m_columns_with_changed_bound.insert(j); - set_low_bound_witness(j, ci); - set_upper_bound_witness(j, ci); - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - } - break; + auto v = numeric_pair(right_side, zero_of_type()); + if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { + m_status = lp_status::INFEASIBLE; + set_low_bound_witness(j, ci); + m_infeasible_column_index = j; + } + else { + m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v; + m_columns_with_changed_bound.insert(j); + set_low_bound_witness(j, ci); + set_upper_bound_witness(j, ci); + m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; + } + break; } break; - default: - lp_unreachable(); + default: + lp_unreachable(); - } + } } void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) { - lp_assert(m_status == lp_status::INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j])); - mpq y_of_bound(0); - switch (kind) { - case LT: - y_of_bound = -1; - case LE: + lp_assert(m_status == lp_status::INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j])); + mpq y_of_bound(0); + switch (kind) { + case LT: + y_of_bound = -1; + case LE: { - auto up = numeric_pair(right_side, y_of_bound); - if (up < m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - set_upper_bound_witness(j, ci); - m_columns_with_changed_bound.insert(j); - } + auto up = numeric_pair(right_side, y_of_bound); + if (up < m_mpq_lar_core_solver.m_r_upper_bounds[j]) { + m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; + set_upper_bound_witness(j, ci); + m_columns_with_changed_bound.insert(j); + } - if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - lp_assert(false); - m_infeasible_column_index = j; - } - else { - if (m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]) - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - } + if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) { + m_status = lp_status::INFEASIBLE; + lp_assert(false); + m_infeasible_column_index = j; + } + else { + if (m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]) + m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; + } } break; - case GT: - y_of_bound = 1; - case GE: + case GT: + y_of_bound = 1; + case GE: { - auto low = numeric_pair(right_side, y_of_bound); - if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_mpq_lar_core_solver.m_r_low_bounds[j] = low; - m_columns_with_changed_bound.insert(j); - set_low_bound_witness(j, ci); - } - if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - } - else if (low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - } + auto low = numeric_pair(right_side, y_of_bound); + if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) { + m_mpq_lar_core_solver.m_r_low_bounds[j] = low; + m_columns_with_changed_bound.insert(j); + set_low_bound_witness(j, ci); + } + if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + } + else if (low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) { + m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; + } } break; - case EQ: + case EQ: { - auto v = numeric_pair(right_side, zero_of_type()); - if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_upper_bound_witness(j, ci); - } - else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_low_bound_witness(j, ci); - } - else { - m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v; - set_low_bound_witness(j, ci); - set_upper_bound_witness(j, ci); - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - m_columns_with_changed_bound.insert(j); - } + auto v = numeric_pair(right_side, zero_of_type()); + if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + set_upper_bound_witness(j, ci); + } + else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + set_low_bound_witness(j, ci); + } + else { + m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v; + set_low_bound_witness(j, ci); + set_upper_bound_witness(j, ci); + m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; + m_columns_with_changed_bound.insert(j); + } - break; + break; } - default: - lp_unreachable(); + default: + lp_unreachable(); - } + } } void lar_solver::update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) { - lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::low_bound); - mpq y_of_bound(0); - switch (kind) { - case LT: - y_of_bound = -1; - case LE: + lp_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::low_bound); + mpq y_of_bound(0); + switch (kind) { + case LT: + y_of_bound = -1; + case LE: { - auto up = numeric_pair(right_side, y_of_bound); - m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; - set_upper_bound_witness(j, ci); - m_columns_with_changed_bound.insert(j); + auto up = numeric_pair(right_side, y_of_bound); + m_mpq_lar_core_solver.m_r_upper_bounds[j] = up; + set_upper_bound_witness(j, ci); + m_columns_with_changed_bound.insert(j); - if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - } - else { - m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j] ? column_type::boxed : column_type::fixed; - } + if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + } + else { + m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j] ? column_type::boxed : column_type::fixed; + } } break; - case GT: - y_of_bound = 1; - case GE: + case GT: + y_of_bound = 1; + case GE: { - auto low = numeric_pair(right_side, y_of_bound); - if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_mpq_lar_core_solver.m_r_low_bounds[j] = low; - m_columns_with_changed_bound.insert(j); - set_low_bound_witness(j, ci); - } + auto low = numeric_pair(right_side, y_of_bound); + if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) { + m_mpq_lar_core_solver.m_r_low_bounds[j] = low; + m_columns_with_changed_bound.insert(j); + set_low_bound_witness(j, ci); + } } break; - case EQ: + case EQ: { - auto v = numeric_pair(right_side, zero_of_type()); - if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_upper_bound_witness(j, ci); - } - else { - m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v; - set_low_bound_witness(j, ci); - set_upper_bound_witness(j, ci); - m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; - } - m_columns_with_changed_bound.insert(j); - break; + auto v = numeric_pair(right_side, zero_of_type()); + if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + set_upper_bound_witness(j, ci); + } + else { + m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v; + set_low_bound_witness(j, ci); + set_upper_bound_witness(j, ci); + m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed; + } + m_columns_with_changed_bound.insert(j); + break; } - default: - lp_unreachable(); + default: + lp_unreachable(); - } + } } void lar_solver::update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) { - lp_assert(m_status == lp_status::INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j])); - lp_assert(m_status == lp_status::INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero())); - auto v = numeric_pair(right_side, mpq(0)); + lp_assert(m_status == lp_status::INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j])); + lp_assert(m_status == lp_status::INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero())); + auto v = numeric_pair(right_side, mpq(0)); - mpq y_of_bound(0); - switch (kind) { - case LT: - if (v <= m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_upper_bound_witness(j, ci); - } - break; - case LE: + mpq y_of_bound(0); + switch (kind) { + case LT: + if (v <= m_mpq_lar_core_solver.m_r_low_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + set_upper_bound_witness(j, ci); + } + break; + case LE: { - if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_upper_bound_witness(j, ci); - } + if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + set_upper_bound_witness(j, ci); + } } break; - case GT: + case GT: { - if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_low_bound_witness(j, ci); - } + if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + set_low_bound_witness(j, ci); + } } break; - case GE: + case GE: { - if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_low_bound_witness(j, ci); - } + if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + set_low_bound_witness(j, ci); + } } break; - case EQ: + case EQ: { - if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_upper_bound_witness(j, ci); - } - else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { - m_status = lp_status::INFEASIBLE; - m_infeasible_column_index = j; - set_low_bound_witness(j, ci); - } - break; + if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + set_upper_bound_witness(j, ci); + } + else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) { + m_status = lp_status::INFEASIBLE; + m_infeasible_column_index = j; + set_low_bound_witness(j, ci); + } + break; } - default: - lp_unreachable(); + default: + lp_unreachable(); - } + } } diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index 47a1936e1..aecc337cf 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -34,7 +34,6 @@ Revision History: #include "util/lp/lp_primal_core_solver.h" #include "util/lp/random_updater.h" #include -#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 m_columns_to_ext_vars_or_term_indices; stacked_vector m_columns_to_ul_pairs; vector m_constraints; +public : + const vector& constraints() const { + return m_constraints; + } +private: stacked_value 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 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>(), - 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 it(A_r().m_rows[row_index]); - SASSERT(use_tableau()); - bound_analyzer_on_row::analyze_row(it, - zero_of_type>(), - 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>& left_side_with_terms, - vector> &left_side, mpq & free_coeff) const; + vector> &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> & explanation, - const vector> & inf_row, - int inf_sign) const; + vector> & explanation, + const vector> & 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> & 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; } + }; } diff --git a/src/util/lp/lar_term.h b/src/util/lp/lar_term.h index 3f83a3fba..a2aee5109 100644 --- a/src/util/lp/lar_term.h +++ b/src/util/lp/lar_term.h @@ -89,6 +89,7 @@ struct lar_term { m_coeffs.clear(); m_v = zero_of_type(); } + }; } diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index 4b239e5de..1ed6e9a3e 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -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;*/ diff --git a/src/util/lp/lp_core_solver_base.hpp b/src/util/lp/lp_core_solver_base.hpp index 51f86a8b5..d54b7d363 100644 --- a/src/util/lp/lp_core_solver_base.hpp +++ b/src/util/lp/lp_core_solver_base.hpp @@ -226,7 +226,7 @@ template bool lp_core_solver_base:: A_mult_x_is_off() const { lp_assert(m_x.size() == m_A.column_count()); if (numeric_traits::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::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 bool lp_core_solver_base:: 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; diff --git a/src/util/lp/lp_primal_core_solver.hpp b/src/util/lp/lp_primal_core_solver.hpp index eda85175b..436a13c5b 100644 --- a/src/util/lp/lp_primal_core_solver.hpp +++ b/src/util/lp/lp_primal_core_solver.hpp @@ -673,7 +673,7 @@ template void lp_primal_core_solver::calc_work template void lp_primal_core_solver::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(); diff --git a/src/util/lp/lp_primal_core_solver_tableau.h b/src/util/lp/lp_primal_core_solver_tableau.h index 77cebc821..cd47abb47 100644 --- a/src/util/lp/lp_primal_core_solver_tableau.h +++ b/src/util/lp/lp_primal_core_solver_tableau.h @@ -206,7 +206,7 @@ unsigned lp_primal_core_solver::solve_with_tableau() { } template void lp_primal_core_solver::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 void lp_primal_core_solver::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 void lp_primal_core_solver::advance_on_en template void lp_primal_core_solver::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 int lp_primal_core_solver::find_leaving_ } template void lp_primal_core_solver::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 void lp_primal_core_solver:: diff --git a/src/util/lp/lp_utils.h b/src/util/lp/lp_utils.h index fb20959e4..f6c98645a 100644 --- a/src/util/lp/lp_utils.h +++ b/src/util/lp/lp_utils.h @@ -23,7 +23,7 @@ Revision History: #include "util/debug.h" #include template -bool try_get_val(const std::unordered_map & map, const A& key, B & val) { +bool try_get_value(const std::unordered_map & map, const A& key, B & val) { const auto it = map.find(key); if (it == map.end()) return false; val = it->second; diff --git a/src/util/lp/numeric_pair.h b/src/util/lp/numeric_pair.h index c5829f03e..ce8dba3a0 100644 --- a/src/util/lp/numeric_pair.h +++ b/src/util/lp/numeric_pair.h @@ -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 { +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 { public: static bool precise() { return false; } @@ -79,6 +91,7 @@ template <> class numeric_traits { 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 diff --git a/src/util/lp/stacked_map.h b/src/util/lp/stacked_map.h index 5ab5ff869..471aa41a0 100644 --- a/src/util/lp/stacked_map.h +++ b/src/util/lp/stacked_map.h @@ -21,29 +21,28 @@ Revision History: #pragma once // this class implements a map with some stack functionality #include +#include #include #include namespace lp { - - -template , - typename KeyEqual = std::equal_to, - typename Allocator = std::allocator< std::pair > > +template, + class _Alloc = std::allocator > > class stacked_map { struct delta { - std::unordered_set m_new; - std::unordered_map m_original_changed; + std::unordered_set m_new; + std::unordered_map m_original_changed; // std::unordered_map m_deb_copy; }; - std::unordered_map m_map; + std::map m_map; std::stack m_stack; public: class ref { - stacked_map & m_map; + stacked_map & m_map; const A & m_key; public: - ref(stacked_map & m, const A & key) :m_map(m), m_key(key) {} + ref(stacked_map & 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& operator()() const { return m_map;} + bool empty() const { return m_map.empty(); } + + const std::map& operator()() const { return m_map;} }; } diff --git a/src/util/lp/static_matrix.h b/src/util/lp/static_matrix.h index 6833352ba..a1c3b64a4 100644 --- a/src/util/lp/static_matrix.h +++ b/src/util/lp/static_matrix.h @@ -62,7 +62,6 @@ class static_matrix dim(unsigned m, unsigned n) :m_m(m), m_n(n) {} }; std::stack m_stack; - vector m_became_zeros; // the row indices that became zeroes during the pivoting public: typedef vector> row_strip; typedef vector 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 & get_row_cell(const column_cell & c) { + return m_rows[c.m_i][c.m_offset]; + } + + column_cell & get_column_cell(const row_cell &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]; diff --git a/src/util/lp/static_matrix.hpp b/src/util/lp/static_matrix.hpp index 4e05c669c..6797d265e 100644 --- a/src/util/lp/static_matrix.hpp +++ b/src/util/lp/static_matrix.hpp @@ -35,57 +35,46 @@ void static_matrix::init_row_columns(unsigned m, unsigned n) { } -template void static_matrix::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 void static_matrix::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 bool static_matrix::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 T static_matrix::get_row_balance(unsi } template bool static_matrix::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 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 bool static_matrix::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 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); diff --git a/src/util/lp/ul_pair.h b/src/util/lp/ul_pair.h index cbf511d90..e0edaadcc 100644 --- a/src/util/lp/ul_pair.h +++ b/src/util/lp/ul_pair.h @@ -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 & a, const std::pair & b) { return a.second < b.second; @@ -69,7 +72,7 @@ public: m_low_bound_witness(static_cast(-1)), m_upper_bound_witness(static_cast(-1)), m_i(static_cast(-1)) -{} + {} ul_pair(row_index ri) : m_low_bound_witness(static_cast(-1)), m_upper_bound_witness(static_cast(-1)), diff --git a/src/util/rational.h b/src/util/rational.h index dc7b79419..930d79010 100644 --- a/src/util/rational.h +++ b/src/util/rational.h @@ -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); }