From 92fe8c59688b89b165287e63cb560a43139e427e Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Fri, 3 Mar 2023 17:53:00 -0800 Subject: [PATCH] restore the previous state Signed-off-by: Lev Nachmanson --- src/math/lp/CMakeLists.txt | 3 + src/math/lp/indexed_value.h | 11 + src/math/lp/lar_core_solver.h | 30 +- src/math/lp/lp_dual_simplex.cpp | 24 + src/math/lp/lp_dual_simplex.h | 93 +++ src/math/lp/lp_dual_simplex_def.h | 376 +++++++++++ src/math/lp/lp_primal_core_solver.h | 1 + src/math/lp/lp_primal_simplex.cpp | 35 + src/math/lp/lp_primal_simplex.h | 106 +++ src/math/lp/lp_primal_simplex_def.h | 367 +++++++++++ src/math/lp/lp_solver.cpp | 55 ++ src/math/lp/lp_solver.h | 260 ++++++++ src/math/lp/lp_solver_def.h | 571 ++++++++++++++++ src/math/lp/mps_reader.h | 891 +++++++++++++++++++++++++ src/math/lp/static_matrix.cpp | 1 + src/sat/smt/arith_sls.h | 3 + src/sat/smt/arith_solver.h | 4 +- src/shell/CMakeLists.txt | 1 + src/shell/lp_frontend.cpp | 109 +++ src/shell/lp_frontend.h | 7 + src/shell/main.cpp | 10 +- src/smt/theory_lra.cpp | 3 + src/test/lp/lp.cpp | 982 ++++++++++++++++++++++++++-- src/test/lp/smt_reader.h | 4 + src/test/lp/test_file_reader.h | 1 + 25 files changed, 3878 insertions(+), 70 deletions(-) create mode 100644 src/math/lp/lp_dual_simplex.cpp create mode 100644 src/math/lp/lp_dual_simplex.h create mode 100644 src/math/lp/lp_dual_simplex_def.h create mode 100644 src/math/lp/lp_primal_simplex.cpp create mode 100644 src/math/lp/lp_primal_simplex.h create mode 100644 src/math/lp/lp_primal_simplex_def.h create mode 100644 src/math/lp/lp_solver.cpp create mode 100644 src/math/lp/lp_solver.h create mode 100644 src/math/lp/lp_solver_def.h create mode 100644 src/math/lp/mps_reader.h create mode 100644 src/shell/lp_frontend.cpp create mode 100644 src/shell/lp_frontend.h diff --git a/src/math/lp/CMakeLists.txt b/src/math/lp/CMakeLists.txt index 5719de44f..9f0fae6bc 100644 --- a/src/math/lp/CMakeLists.txt +++ b/src/math/lp/CMakeLists.txt @@ -20,8 +20,11 @@ z3_add_component(lp lar_core_solver.cpp lp_core_solver_base.cpp lp_dual_core_solver.cpp + lp_dual_simplex.cpp lp_primal_core_solver.cpp + lp_primal_simplex.cpp lp_settings.cpp + lp_solver.cpp lu.cpp lp_utils.cpp matrix.cpp diff --git a/src/math/lp/indexed_value.h b/src/math/lp/indexed_value.h index 92d8f2adf..c48376470 100644 --- a/src/math/lp/indexed_value.h +++ b/src/math/lp/indexed_value.h @@ -43,4 +43,15 @@ public: m_value = val; } }; +#ifdef Z3DEBUG +template +bool check_vector_for_small_values(indexed_vector & w, lp_settings & settings) { + for (unsigned i : w.m_index) { + const X & v = w[i]; + if ((!is_zero(v)) && settings.abs_val_is_smaller_than_drop_tolerance(v)) + return false; + } + return true; +} +#endif } diff --git a/src/math/lp/lar_core_solver.h b/src/math/lp/lar_core_solver.h index de8fe68ad..8a6c64ef0 100644 --- a/src/math/lp/lar_core_solver.h +++ b/src/math/lp/lar_core_solver.h @@ -651,7 +651,35 @@ public: } } - + void scale_problem_for_doubles( + static_matrix& A, + vector & lower_bounds, + vector & upper_bounds) { + vector column_scale_vector; + vector right_side_vector(A.column_count()); + settings().reps_in_scaler = 5; + scaler scaler(right_side_vector, + A, + settings().scaling_minimum, + settings().scaling_maximum, + column_scale_vector, + settings()); + if (! scaler.scale()) { + // the scale did not succeed, unscaling + A.clear(); + create_double_matrix(A); + } else { + for (unsigned j = 0; j < A.column_count(); j++) { + if (m_r_solver.column_has_upper_bound(j)) { + upper_bounds[j] /= column_scale_vector[j]; + } + if (m_r_solver.column_has_lower_bound(j)) { + lower_bounds[j] /= column_scale_vector[j]; + } + } + } + + } // returns the trace of basis changes vector find_solution_signature_with_doubles(lar_solution_signature & signature) { if (m_d_solver.m_factorization == nullptr || m_d_solver.m_factorization->get_status() != LU_status::OK) { diff --git a/src/math/lp/lp_dual_simplex.cpp b/src/math/lp/lp_dual_simplex.cpp new file mode 100644 index 000000000..aaf612f56 --- /dev/null +++ b/src/math/lp/lp_dual_simplex.cpp @@ -0,0 +1,24 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#include "math/lp/lp_dual_simplex_def.h" +template lp::mpq lp::lp_dual_simplex::get_current_cost() const; +template void lp::lp_dual_simplex::find_maximal_solution(); +template double lp::lp_dual_simplex::get_current_cost() const; +template void lp::lp_dual_simplex::find_maximal_solution(); diff --git a/src/math/lp/lp_dual_simplex.h b/src/math/lp/lp_dual_simplex.h new file mode 100644 index 000000000..75ef87492 --- /dev/null +++ b/src/math/lp/lp_dual_simplex.h @@ -0,0 +1,93 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once +#include "util/vector.h" +#include "math/lp/lp_utils.h" +#include "math/lp/lp_solver.h" +#include "math/lp/lp_dual_core_solver.h" +namespace lp { + +template +class lp_dual_simplex: public lp_solver { + lp_dual_core_solver * m_core_solver; + vector m_b_copy; + vector m_lower_bounds; // We don't have a convention here that all low bounds are zeros. At least it does not hold for the first stage solver + vector m_column_types_of_core_solver; + vector m_column_types_of_logicals; + vector m_can_enter_basis; +public: + ~lp_dual_simplex() override { + delete m_core_solver; + } + + lp_dual_simplex() : m_core_solver(nullptr) {} + + + void decide_on_status_after_stage1(); + + void fix_logical_for_stage2(unsigned j); + + void fix_structural_for_stage2(unsigned j); + + void unmark_boxed_and_fixed_columns_and_fix_structural_costs(); + + void restore_right_sides(); + + void solve_for_stage2(); + + void fill_x_with_zeros(); + + void stage1(); + + void stage2(); + + void fill_first_stage_solver_fields(); + + column_type get_column_type(unsigned j); + + void fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(unsigned j); + + void fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(unsigned j); + + void fill_costs_and_bounds_and_column_types_for_the_first_stage_solver(); + + void set_type_for_logical(unsigned j, column_type col_type) { + this->m_column_types_of_logicals[j - this->number_of_core_structurals()] = col_type; + } + + void fill_first_stage_solver_fields_for_row_slack_and_artificial(unsigned row, + unsigned & slack_var, + unsigned & artificial); + + void augment_matrix_A_and_fill_x_and_allocate_some_fields(); + + + + void copy_m_b_aside_and_set_it_to_zeros(); + + void find_maximal_solution() override; + + T get_column_value(unsigned column) const override { + return this->get_column_value_with_core_solver(column, m_core_solver); + } + + T get_current_cost() const override; +}; +} diff --git a/src/math/lp/lp_dual_simplex_def.h b/src/math/lp/lp_dual_simplex_def.h new file mode 100644 index 000000000..8af9d87c1 --- /dev/null +++ b/src/math/lp/lp_dual_simplex_def.h @@ -0,0 +1,376 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once + +#include "math/lp/lp_dual_simplex.h" +namespace lp{ + +template void lp_dual_simplex::decide_on_status_after_stage1() { + switch (m_core_solver->get_status()) { + case lp_status::OPTIMAL: + if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) { + this->m_status = lp_status::FEASIBLE; + } else { + this->m_status = lp_status::UNBOUNDED; + } + break; + case lp_status::DUAL_UNBOUNDED: + lp_unreachable(); + case lp_status::TIME_EXHAUSTED: + this->m_status = lp_status::TIME_EXHAUSTED; + break; + case lp_status::FLOATING_POINT_ERROR: + this->m_status = lp_status::FLOATING_POINT_ERROR; + break; + default: + lp_unreachable(); + } +} + +template void lp_dual_simplex::fix_logical_for_stage2(unsigned j) { + lp_assert(j >= this->number_of_core_structurals()); + switch (m_column_types_of_logicals[j - this->number_of_core_structurals()]) { + case column_type::lower_bound: + m_lower_bounds[j] = numeric_traits::zero(); + m_column_types_of_core_solver[j] = column_type::lower_bound; + m_can_enter_basis[j] = true; + break; + case column_type::fixed: + this->m_upper_bounds[j] = m_lower_bounds[j] = numeric_traits::zero(); + m_column_types_of_core_solver[j] = column_type::fixed; + m_can_enter_basis[j] = false; + break; + default: + lp_unreachable(); + } +} + +template void lp_dual_simplex::fix_structural_for_stage2(unsigned j) { + column_info * ci = this->m_map_from_var_index_to_column_info[this->m_core_solver_columns_to_external_columns[j]]; + switch (ci->get_column_type()) { + case column_type::lower_bound: + m_lower_bounds[j] = numeric_traits::zero(); + m_column_types_of_core_solver[j] = column_type::lower_bound; + m_can_enter_basis[j] = true; + break; + case column_type::fixed: + case column_type::upper_bound: + lp_unreachable(); + case column_type::boxed: + this->m_upper_bounds[j] = ci->get_adjusted_upper_bound() / this->m_column_scale[j]; + m_lower_bounds[j] = numeric_traits::zero(); + m_column_types_of_core_solver[j] = column_type::boxed; + m_can_enter_basis[j] = true; + break; + case column_type::free_column: + m_can_enter_basis[j] = true; + m_column_types_of_core_solver[j] = column_type::free_column; + break; + default: + lp_unreachable(); + } + // T cost_was = this->m_costs[j]; + this->set_scaled_cost(j); +} + +template void lp_dual_simplex::unmark_boxed_and_fixed_columns_and_fix_structural_costs() { + unsigned j = this->m_A->column_count(); + while (j-- > this->number_of_core_structurals()) { + fix_logical_for_stage2(j); + } + j = this->number_of_core_structurals(); + while (j--) { + fix_structural_for_stage2(j); + } +} + +template void lp_dual_simplex::restore_right_sides() { + unsigned i = this->m_A->row_count(); + while (i--) { + this->m_b[i] = m_b_copy[i]; + } +} + +template void lp_dual_simplex::solve_for_stage2() { + m_core_solver->restore_non_basis(); + m_core_solver->solve_yB(m_core_solver->m_y); + m_core_solver->fill_reduced_costs_from_m_y_by_rows(); + m_core_solver->start_with_initial_basis_and_make_it_dual_feasible(); + m_core_solver->set_status(lp_status::FEASIBLE); + m_core_solver->solve(); + switch (m_core_solver->get_status()) { + case lp_status::OPTIMAL: + this->m_status = lp_status::OPTIMAL; + break; + case lp_status::DUAL_UNBOUNDED: + this->m_status = lp_status::INFEASIBLE; + break; + case lp_status::TIME_EXHAUSTED: + this->m_status = lp_status::TIME_EXHAUSTED; + break; + case lp_status::FLOATING_POINT_ERROR: + this->m_status = lp_status::FLOATING_POINT_ERROR; + break; + default: + lp_unreachable(); + } + this->m_second_stage_iterations = m_core_solver->total_iterations(); + this->m_total_iterations = (this->m_first_stage_iterations + this->m_second_stage_iterations); +} + +template void lp_dual_simplex::fill_x_with_zeros() { + unsigned j = this->m_A->column_count(); + while (j--) { + this->m_x[j] = numeric_traits::zero(); + } +} + +template void lp_dual_simplex::stage1() { + lp_assert(m_core_solver == nullptr); + this->m_x.resize(this->m_A->column_count(), numeric_traits::zero()); + if (this->m_settings.get_message_ostream() != nullptr) + this->print_statistics_on_A(*this->m_settings.get_message_ostream()); + m_core_solver = new lp_dual_core_solver( + *this->m_A, + m_can_enter_basis, + this->m_b, // the right side vector + this->m_x, + this->m_basis, + this->m_nbasis, + this->m_heading, + this->m_costs, + this->m_column_types_of_core_solver, + this->m_lower_bounds, + this->m_upper_bounds, + this->m_settings, + *this); + m_core_solver->fill_reduced_costs_from_m_y_by_rows(); + m_core_solver->start_with_initial_basis_and_make_it_dual_feasible(); + if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) { + // skipping stage 1 + m_core_solver->set_status(lp_status::OPTIMAL); + m_core_solver->set_total_iterations(0); + } else { + m_core_solver->solve(); + } + decide_on_status_after_stage1(); + this->m_first_stage_iterations = m_core_solver->total_iterations(); +} + +template void lp_dual_simplex::stage2() { + unmark_boxed_and_fixed_columns_and_fix_structural_costs(); + restore_right_sides(); + solve_for_stage2(); +} + +template void lp_dual_simplex::fill_first_stage_solver_fields() { + unsigned slack_var = this->number_of_core_structurals(); + unsigned artificial = this->number_of_core_structurals() + this->m_slacks; + + for (unsigned row = 0; row < this->row_count(); row++) { + fill_first_stage_solver_fields_for_row_slack_and_artificial(row, slack_var, artificial); + } + fill_costs_and_bounds_and_column_types_for_the_first_stage_solver(); +} + +template column_type lp_dual_simplex::get_column_type(unsigned j) { + lp_assert(j < this->m_A->column_count()); + if (j >= this->number_of_core_structurals()) { + return m_column_types_of_logicals[j - this->number_of_core_structurals()]; + } + return this->m_map_from_var_index_to_column_info[this->m_core_solver_columns_to_external_columns[j]]->get_column_type(); +} + +template void lp_dual_simplex::fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(unsigned j) { + // see 4.7 in the dissertation of Achim Koberstein + lp_assert(this->m_core_solver_columns_to_external_columns.find(j) != + this->m_core_solver_columns_to_external_columns.end()); + + T free_bound = T(1e4); // see 4.8 + unsigned jj = this->m_core_solver_columns_to_external_columns[j]; + lp_assert(this->m_map_from_var_index_to_column_info.find(jj) != this->m_map_from_var_index_to_column_info.end()); + column_info * ci = this->m_map_from_var_index_to_column_info[jj]; + switch (ci->get_column_type()) { + case column_type::upper_bound: { + std::stringstream s; + s << "unexpected bound type " << j << " " + << column_type_to_string(get_column_type(j)); + throw_exception(s.str()); + break; + } + case column_type::lower_bound: { + m_can_enter_basis[j] = true; + this->set_scaled_cost(j); + this->m_lower_bounds[j] = numeric_traits::zero(); + this->m_upper_bounds[j] = numeric_traits::one(); + break; + } + case column_type::free_column: { + m_can_enter_basis[j] = true; + this->set_scaled_cost(j); + this->m_upper_bounds[j] = free_bound; + this->m_lower_bounds[j] = -free_bound; + break; + } + case column_type::boxed: + m_can_enter_basis[j] = false; + this->m_costs[j] = numeric_traits::zero(); + this->m_upper_bounds[j] = this->m_lower_bounds[j] = numeric_traits::zero(); // is it needed? + break; + default: + lp_unreachable(); + } + m_column_types_of_core_solver[j] = column_type::boxed; +} + +template void lp_dual_simplex::fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(unsigned j) { + this->m_costs[j] = 0; + lp_assert(get_column_type(j) != column_type::upper_bound); + if ((m_can_enter_basis[j] = (get_column_type(j) == column_type::lower_bound))) { + m_column_types_of_core_solver[j] = column_type::boxed; + this->m_lower_bounds[j] = numeric_traits::zero(); + this->m_upper_bounds[j] = numeric_traits::one(); + } else { + m_column_types_of_core_solver[j] = column_type::fixed; + this->m_lower_bounds[j] = numeric_traits::zero(); + this->m_upper_bounds[j] = numeric_traits::zero(); + } +} + +template void lp_dual_simplex::fill_costs_and_bounds_and_column_types_for_the_first_stage_solver() { + unsigned j = this->m_A->column_count(); + while (j-- > this->number_of_core_structurals()) { // go over logicals here + fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(j); + } + j = this->number_of_core_structurals(); + while (j--) { + fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(j); + } +} + +template void lp_dual_simplex::fill_first_stage_solver_fields_for_row_slack_and_artificial(unsigned row, + unsigned & slack_var, + unsigned & artificial) { + lp_assert(row < this->row_count()); + auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[row]]; + // we need to bring the program to the form Ax = b + T rs = this->m_b[row]; + switch (constraint.m_relation) { + case Equal: // no slack variable here + set_type_for_logical(artificial, column_type::fixed); + this->m_basis[row] = artificial; + this->m_costs[artificial] = numeric_traits::zero(); + (*this->m_A)(row, artificial) = numeric_traits::one(); + artificial++; + break; + + case Greater_or_equal: + set_type_for_logical(slack_var, column_type::lower_bound); + (*this->m_A)(row, slack_var) = - numeric_traits::one(); + if (rs > 0) { + // adding one artificial + set_type_for_logical(artificial, column_type::fixed); + (*this->m_A)(row, artificial) = numeric_traits::one(); + this->m_basis[row] = artificial; + this->m_costs[artificial] = numeric_traits::zero(); + artificial++; + } else { + // we can put a slack_var into the basis, and avoid adding an artificial variable + this->m_basis[row] = slack_var; + this->m_costs[slack_var] = numeric_traits::zero(); + } + slack_var++; + break; + case Less_or_equal: + // introduce a non-negative slack variable + set_type_for_logical(slack_var, column_type::lower_bound); + (*this->m_A)(row, slack_var) = numeric_traits::one(); + if (rs < 0) { + // adding one artificial + set_type_for_logical(artificial, column_type::fixed); + (*this->m_A)(row, artificial) = - numeric_traits::one(); + this->m_basis[row] = artificial; + this->m_costs[artificial] = numeric_traits::zero(); + artificial++; + } else { + // we can put slack_var into the basis, and avoid adding an artificial variable + this->m_basis[row] = slack_var; + this->m_costs[slack_var] = numeric_traits::zero(); + } + slack_var++; + break; + } +} + +template void lp_dual_simplex::augment_matrix_A_and_fill_x_and_allocate_some_fields() { + this->count_slacks_and_artificials(); + this->m_A->add_columns_at_the_end(this->m_slacks + this->m_artificials); + unsigned n = this->m_A->column_count(); + this->m_column_types_of_core_solver.resize(n); + m_column_types_of_logicals.resize(this->m_slacks + this->m_artificials); + this->m_costs.resize(n); + this->m_upper_bounds.resize(n); + this->m_lower_bounds.resize(n); + m_can_enter_basis.resize(n); + this->m_basis.resize(this->m_A->row_count()); +} + + + +template void lp_dual_simplex::copy_m_b_aside_and_set_it_to_zeros() { + for (unsigned i = 0; i < this->m_b.size(); i++) { + m_b_copy.push_back(this->m_b[i]); + this->m_b[i] = numeric_traits::zero(); // preparing for the first stage + } +} + +template void lp_dual_simplex::find_maximal_solution(){ + if (this->problem_is_empty()) { + this->m_status = lp_status::EMPTY; + return; + } + + this->flip_costs(); // do it for now, todo ( remove the flipping) + + this->cleanup(); + if (this->m_status == lp_status::INFEASIBLE) { + return; + } + this->fill_matrix_A_and_init_right_side(); + this->fill_m_b(); + this->scale(); + augment_matrix_A_and_fill_x_and_allocate_some_fields(); + fill_first_stage_solver_fields(); + copy_m_b_aside_and_set_it_to_zeros(); + stage1(); + if (this->m_status == lp_status::FEASIBLE) { + stage2(); + } +} + + +template T lp_dual_simplex::get_current_cost() const { + T ret = numeric_traits::zero(); + for (auto it : this->m_map_from_var_index_to_column_info) { + ret += this->get_column_cost_value(it.first, it.second); + } + return -ret; // we flip costs for now +} +} diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h index a60395ab0..4b6163df8 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -30,6 +30,7 @@ Revision History: #include #include #include "math/lp/lu.h" +#include "math/lp/lp_solver.h" #include "math/lp/static_matrix.h" #include "math/lp/core_solver_pretty_printer.h" #include "math/lp/lp_core_solver_base.h" diff --git a/src/math/lp/lp_primal_simplex.cpp b/src/math/lp/lp_primal_simplex.cpp new file mode 100644 index 000000000..634f52900 --- /dev/null +++ b/src/math/lp/lp_primal_simplex.cpp @@ -0,0 +1,35 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#include +#include +#include +#include "util/vector.h" +#include +#include "math/lp/lp_primal_simplex_def.h" +template bool lp::lp_primal_simplex::bounds_hold(std::unordered_map, std::equal_to, std::allocator > > const&); +template bool lp::lp_primal_simplex::row_constraints_hold(std::unordered_map, std::equal_to, std::allocator > > const&); +template double lp::lp_primal_simplex::get_current_cost() const; +template double lp::lp_primal_simplex::get_column_value(unsigned int) const; +template lp::lp_primal_simplex::~lp_primal_simplex(); +template lp::lp_primal_simplex::~lp_primal_simplex(); +template lp::mpq lp::lp_primal_simplex::get_current_cost() const; +template lp::mpq lp::lp_primal_simplex::get_column_value(unsigned int) const; +template void lp::lp_primal_simplex::find_maximal_solution(); +template void lp::lp_primal_simplex::find_maximal_solution(); diff --git a/src/math/lp/lp_primal_simplex.h b/src/math/lp/lp_primal_simplex.h new file mode 100644 index 000000000..77e12d088 --- /dev/null +++ b/src/math/lp/lp_primal_simplex.h @@ -0,0 +1,106 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once +#include "util/vector.h" +#include +#include +#include +#include "math/lp/lp_utils.h" +#include "math/lp/column_info.h" +#include "math/lp/lp_primal_core_solver.h" +#include "math/lp/lp_solver.h" +namespace lp { +template +class lp_primal_simplex: public lp_solver { + lp_primal_core_solver * m_core_solver; + vector m_lower_bounds; +private: + unsigned original_rows() { return this->m_external_rows_to_core_solver_rows.size(); } + + void fill_costs_and_x_for_first_stage_solver(unsigned original_number_of_columns); + + void init_buffer(unsigned k, vector & r); + + void refactor(); + + void set_scaled_costs(); +public: + lp_primal_simplex(): m_core_solver(nullptr) {} + + column_info * get_or_create_column_info(unsigned column); + + void set_status(lp_status status) { + this->m_status = status; + } + + lp_status get_status() { + return this->m_status; + } + + void fill_acceptable_values_for_x(); + + + void set_zero_bound(bool * bound_is_set, T * bounds, unsigned i); + + void fill_costs_and_x_for_first_stage_solver_for_row( + int row, + unsigned & slack_var, + unsigned & artificial); + + + + + void set_core_solver_bounds(); + + void find_maximal_solution() override; + + void fill_A_x_and_basis_for_stage_one_total_inf(); + + void fill_A_x_and_basis_for_stage_one_total_inf_for_row(unsigned row); + + void solve_with_total_inf(); + + + ~lp_primal_simplex() override; + + bool bounds_hold(std::unordered_map const & solution); + + T get_row_value(unsigned i, std::unordered_map const & solution, std::ostream * out); + + bool row_constraint_holds(unsigned i, std::unordered_map const & solution, std::ostream * out); + + bool row_constraints_hold(std::unordered_map const & solution); + + + T * get_array_from_map(std::unordered_map const & solution); + + bool solution_is_feasible(std::unordered_map const & solution) { + return bounds_hold(solution) && row_constraints_hold(solution); + } + + T get_column_value(unsigned column) const override { + return this->get_column_value_with_core_solver(column, m_core_solver); + } + + T get_current_cost() const override; + + +}; +} diff --git a/src/math/lp/lp_primal_simplex_def.h b/src/math/lp/lp_primal_simplex_def.h new file mode 100644 index 000000000..7ffe819b2 --- /dev/null +++ b/src/math/lp/lp_primal_simplex_def.h @@ -0,0 +1,367 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once + +#include +#include "util/vector.h" +#include "math/lp/lp_primal_simplex.h" + +namespace lp { +template void lp_primal_simplex::fill_costs_and_x_for_first_stage_solver(unsigned original_number_of_columns) { + unsigned slack_var = original_number_of_columns; + unsigned artificial = original_number_of_columns + this->m_slacks; + + for (unsigned row = 0; row < this->row_count(); row++) { + fill_costs_and_x_for_first_stage_solver_for_row(row, slack_var, artificial); + } +} + +template void lp_primal_simplex::init_buffer(unsigned k, vector & r) { + for (unsigned i = 0; i < k; i++) { + r[i] = 0; + } + r[k] = 1; + for (unsigned i = this->row_count() -1; i > k; i--) { + r[i] = 0; + } +} + +template void lp_primal_simplex::refactor() { + m_core_solver->init_lu(); + if (m_core_solver->factorization()->get_status() != LU_status::OK) { + throw_exception("cannot refactor"); + } +} + +template void lp_primal_simplex::set_scaled_costs() { + unsigned j = this->number_of_core_structurals(); + while (j-- > 0) { + this->set_scaled_cost(j); + } +} + +template column_info * lp_primal_simplex::get_or_create_column_info(unsigned column) { + auto it = this->m_columns.find(column); + return (it == this->m_columns.end())? ( this->m_columns[column] = new column_info) : it->second; +} + +template void lp_primal_simplex::fill_acceptable_values_for_x() { + for (auto t : this->m_core_solver_columns_to_external_columns) { + this->m_x[t.first] = numeric_traits::zero(); + } +} + + +template void lp_primal_simplex::set_zero_bound(bool * bound_is_set, T * bounds, unsigned i) { + bound_is_set[i] = true; + bounds[i] = numeric_traits::zero(); +} + +template void lp_primal_simplex::fill_costs_and_x_for_first_stage_solver_for_row( + int row, + unsigned & slack_var, + unsigned & artificial) { + lp_assert(row >= 0 && row < this->row_count()); + auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[row]]; + // we need to bring the program to the form Ax = b + T rs = this->m_b[row]; + T artificial_cost = - numeric_traits::one(); + switch (constraint.m_relation) { + case Equal: // no slack variable here + this->m_column_types[artificial] = column_type::lower_bound; + this->m_costs[artificial] = artificial_cost; // we are maximizing, so the artificial, which is non-negatiive, will be pushed to zero + this->m_basis[row] = artificial; + if (rs >= 0) { + (*this->m_A)(row, artificial) = numeric_traits::one(); + this->m_x[artificial] = rs; + } else { + (*this->m_A)(row, artificial) = - numeric_traits::one(); + this->m_x[artificial] = - rs; + } + artificial++; + break; + + case Greater_or_equal: + this->m_column_types[slack_var] = column_type::lower_bound; + (*this->m_A)(row, slack_var) = - numeric_traits::one(); + + if (rs > 0) { + lp_assert(numeric_traits::is_zero(this->m_x[slack_var])); + // adding one artificial + this->m_column_types[artificial] = column_type::lower_bound; + (*this->m_A)(row, artificial) = numeric_traits::one(); + this->m_costs[artificial] = artificial_cost; + this->m_basis[row] = artificial; + this->m_x[artificial] = rs; + artificial++; + } else { + // we can put a slack_var into the basis, and atemplate void lp_primal_simplex::adding an artificial variable + this->m_basis[row] = slack_var; + this->m_x[slack_var] = - rs; + } + slack_var++; + break; + case Less_or_equal: + // introduce a non-negative slack variable + this->m_column_types[slack_var] = column_type::lower_bound; + (*this->m_A)(row, slack_var) = numeric_traits::one(); + + if (rs < 0) { + // adding one artificial + lp_assert(numeric_traits::is_zero(this->m_x[slack_var])); + this->m_column_types[artificial] = column_type::lower_bound; + (*this->m_A)(row, artificial) = - numeric_traits::one(); + this->m_costs[artificial] = artificial_cost; + this->m_x[artificial] = - rs; + this->m_basis[row] = artificial++; + } else { + // we can put slack_var into the basis, and atemplate void lp_primal_simplex::adding an artificial variable + this->m_basis[row] = slack_var; + this->m_x[slack_var] = rs; + } + slack_var++; + break; + } +} + + + + + +template void lp_primal_simplex::set_core_solver_bounds() { + unsigned total_vars = this->m_A->column_count() + this->m_slacks + this->m_artificials; + this->m_column_types.resize(total_vars); + this->m_upper_bounds.resize(total_vars); + for (auto cit : this->m_map_from_var_index_to_column_info) { + column_info * ci = cit.second; + unsigned j = ci->get_column_index(); + if (!is_valid(j)) + continue; // the variable is not mapped to a column + switch (this->m_column_types[j] = ci->get_column_type()){ + case column_type::fixed: + this->m_upper_bounds[j] = numeric_traits::zero(); + break; + case column_type::boxed: + this->m_upper_bounds[j] = ci->get_adjusted_upper_bound() / this->m_column_scale[j]; + break; + + default: break; // do nothing + } + } +} + + +template void lp_primal_simplex::find_maximal_solution() { + if (this->problem_is_empty()) { + this->m_status = lp_status::EMPTY; + return; + } + + this->cleanup(); + this->fill_matrix_A_and_init_right_side(); + if (this->m_status == lp_status::INFEASIBLE) { + return; + } + this->m_x.resize(this->m_A->column_count()); + this->fill_m_b(); + this->scale(); + fill_acceptable_values_for_x(); + this->count_slacks_and_artificials(); + set_core_solver_bounds(); + solve_with_total_inf(); +} + +template void lp_primal_simplex::fill_A_x_and_basis_for_stage_one_total_inf() { + for (unsigned row = 0; row < this->row_count(); row++) + fill_A_x_and_basis_for_stage_one_total_inf_for_row(row); +} + +template void lp_primal_simplex::fill_A_x_and_basis_for_stage_one_total_inf_for_row(unsigned row) { + lp_assert(row < this->row_count()); + auto ext_row_it = this->m_core_solver_rows_to_external_rows.find(row); + lp_assert(ext_row_it != this->m_core_solver_rows_to_external_rows.end()); + unsigned ext_row = ext_row_it->second; + auto constr_it = this->m_constraints.find(ext_row); + lp_assert(constr_it != this->m_constraints.end()); + auto & constraint = constr_it->second; + unsigned j = this->m_A->column_count(); // j is a slack variable + this->m_A->add_column(); + // we need to bring the program to the form Ax = b + this->m_basis[row] = j; + switch (constraint.m_relation) { + case Equal: + this->m_x[j] = this->m_b[row]; + (*this->m_A)(row, j) = numeric_traits::one(); + this->m_column_types[j] = column_type::fixed; + this->m_upper_bounds[j] = m_lower_bounds[j] = zero_of_type(); + break; + + case Greater_or_equal: + this->m_x[j] = - this->m_b[row]; + (*this->m_A)(row, j) = - numeric_traits::one(); + this->m_column_types[j] = column_type::lower_bound; + this->m_upper_bounds[j] = zero_of_type(); + break; + case Less_or_equal: + this->m_x[j] = this->m_b[row]; + (*this->m_A)(row, j) = numeric_traits::one(); + this->m_column_types[j] = column_type::lower_bound; + this->m_upper_bounds[j] = m_lower_bounds[j] = zero_of_type(); + break; + default: + lp_unreachable(); + } +} + +template void lp_primal_simplex::solve_with_total_inf() { + int total_vars = this->m_A->column_count() + this->row_count(); + if (total_vars == 0) { + this->m_status = lp_status::OPTIMAL; + return; + } + m_lower_bounds.clear(); + m_lower_bounds.resize(total_vars, zero_of_type()); // low bounds are shifted ot zero + this->m_x.resize(total_vars, numeric_traits::zero()); + this->m_basis.resize(this->row_count()); + this->m_costs.clear(); + this->m_costs.resize(total_vars, zero_of_type()); + fill_A_x_and_basis_for_stage_one_total_inf(); + if (this->m_settings.get_message_ostream() != nullptr) + this->print_statistics_on_A(*this->m_settings.get_message_ostream()); + set_scaled_costs(); + + m_core_solver = new lp_primal_core_solver(*this->m_A, + this->m_b, + this->m_x, + this->m_basis, + this->m_nbasis, + this->m_heading, + this->m_costs, + this->m_column_types, + m_lower_bounds, + this->m_upper_bounds, + this->m_settings, *this); + m_core_solver->solve(); + this->set_status(m_core_solver->get_status()); + this->m_total_iterations = m_core_solver->total_iterations(); +} + + +template lp_primal_simplex::~lp_primal_simplex() { + delete m_core_solver; +} + +template bool lp_primal_simplex::bounds_hold(std::unordered_map const & solution) { + for (auto it : this->m_map_from_var_index_to_column_info) { + auto sol_it = solution.find(it.second->get_name()); + if (sol_it == solution.end()) { + std::stringstream s; + s << "cannot find column " << it.first << " in solution"; + throw_exception(s.str() ); + } + + if (!it.second->bounds_hold(sol_it->second)) { + it.second->bounds_hold(sol_it->second); + return false; + } + } + return true; +} + +template T lp_primal_simplex::get_row_value(unsigned i, std::unordered_map const & solution, std::ostream * out) { + auto it = this->m_A_values.find(i); + if (it == this->m_A_values.end()) { + std::stringstream s; + s << "cannot find row " << i; + throw_exception(s.str() ); + } + T ret = numeric_traits::zero(); + for (auto & pair : it->second) { + auto cit = this->m_map_from_var_index_to_column_info.find(pair.first); + lp_assert(cit != this->m_map_from_var_index_to_column_info.end()); + column_info * ci = cit->second; + auto sol_it = solution.find(ci->get_name()); + lp_assert(sol_it != solution.end()); + T column_val = sol_it->second; + if (out != nullptr) { + (*out) << pair.second << "(" << ci->get_name() << "=" << column_val << ") "; + } + ret += pair.second * column_val; + } + if (out != nullptr) { + (*out) << " = " << ret << std::endl; + } + return ret; +} + +template bool lp_primal_simplex::row_constraint_holds(unsigned i, std::unordered_map const & solution, std::ostream *out) { + T row_val = get_row_value(i, solution, out); + auto & constraint = this->m_constraints[i]; + T rs = constraint.m_rs; + bool print = out != nullptr; + switch (constraint.m_relation) { + case Equal: + if (fabs(numeric_traits::get_double(row_val - rs)) > 0.00001) { + if (print) { + (*out) << "should be = " << rs << std::endl; + } + return false; + } + return true; + case Greater_or_equal: + if (numeric_traits::get_double(row_val - rs) < -0.00001) { + if (print) { + (*out) << "should be >= " << rs << std::endl; + } + return false; + } + return true;; + + case Less_or_equal: + if (numeric_traits::get_double(row_val - rs) > 0.00001) { + if (print) { + (*out) << "should be <= " << rs << std::endl; + } + return false; + } + return true;; + } + lp_unreachable(); + return false; // it is unreachable +} + +template bool lp_primal_simplex::row_constraints_hold(std::unordered_map const & solution) { + for (auto it : this->m_A_values) { + if (!row_constraint_holds(it.first, solution, nullptr)) { + row_constraint_holds(it.first, solution, nullptr); + return false; + } + } + return true; +} + +template T lp_primal_simplex::get_current_cost() const { + T ret = numeric_traits::zero(); + for (auto it : this->m_map_from_var_index_to_column_info) { + ret += this->get_column_cost_value(it.first, it.second); + } + return ret; +} +} diff --git a/src/math/lp/lp_solver.cpp b/src/math/lp/lp_solver.cpp new file mode 100644 index 000000000..fc9514098 --- /dev/null +++ b/src/math/lp/lp_solver.cpp @@ -0,0 +1,55 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#include +#include "math/lp/lp_solver_def.h" +template void lp::lp_solver::add_constraint(lp::lp_relation, double, unsigned int); +template void lp::lp_solver::cleanup(); +template void lp::lp_solver::count_slacks_and_artificials(); +template void lp::lp_solver::fill_m_b(); +template void lp::lp_solver::fill_matrix_A_and_init_right_side(); +template void lp::lp_solver::flip_costs(); +template double lp::lp_solver::get_column_cost_value(unsigned int, lp::column_info*) const; +template int lp::lp_solver::get_column_index_by_name(std::string) const; +template double lp::lp_solver::get_column_value_with_core_solver(unsigned int, lp::lp_core_solver_base*) const; +template lp::column_info* lp::lp_solver::get_or_create_column_info(unsigned int); +template void lp::lp_solver::give_symbolic_name_to_column(std::string, unsigned int); +template void lp::lp_solver::print_statistics_on_A(std::ostream & out); +template bool lp::lp_solver::problem_is_empty(); +template void lp::lp_solver::scale(); +template void lp::lp_solver::set_scaled_cost(unsigned int); +template lp::lp_solver::~lp_solver(); +template void lp::lp_solver::add_constraint(lp::lp_relation, lp::mpq, unsigned int); +template void lp::lp_solver::cleanup(); +template void lp::lp_solver::count_slacks_and_artificials(); +template void lp::lp_solver::fill_m_b(); +template void lp::lp_solver::fill_matrix_A_and_init_right_side(); +template void lp::lp_solver::flip_costs(); +template lp::mpq lp::lp_solver::get_column_cost_value(unsigned int, lp::column_info*) const; +template int lp::lp_solver::get_column_index_by_name(std::string) const; +template lp::mpq lp::lp_solver::get_column_value_by_name(std::string) const; +template lp::mpq lp::lp_solver::get_column_value_with_core_solver(unsigned int, lp::lp_core_solver_base*) const; +template lp::column_info* lp::lp_solver::get_or_create_column_info(unsigned int); +template void lp::lp_solver::give_symbolic_name_to_column(std::string, unsigned int); +template void lp::lp_solver::print_statistics_on_A(std::ostream & out); +template bool lp::lp_solver::problem_is_empty(); +template void lp::lp_solver::scale(); +template void lp::lp_solver::set_scaled_cost(unsigned int); +template lp::lp_solver::~lp_solver(); +template double lp::lp_solver::get_column_value_by_name(std::string) const; diff --git a/src/math/lp/lp_solver.h b/src/math/lp/lp_solver.h new file mode 100644 index 000000000..ab16a686f --- /dev/null +++ b/src/math/lp/lp_solver.h @@ -0,0 +1,260 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ + +#pragma once +#include +#include +#include +#include "util/vector.h" +#include "math/lp/lp_settings.h" +#include "math/lp/column_info.h" +#include "math/lp/static_matrix.h" +#include "math/lp/lp_core_solver_base.h" +#include "math/lp/scaler.h" +#include "math/lp/bound_analyzer_on_row.h" +namespace lp { +enum lp_relation { + Less_or_equal, + Equal, + Greater_or_equal +}; + +template +struct lp_constraint { + X m_rs; // right side of the constraint + lp_relation m_relation; + lp_constraint() {} // empty constructor + lp_constraint(T rs, lp_relation relation): m_rs(rs), m_relation(relation) {} +}; + + +template +class lp_solver : public column_namer { + column_info * get_or_create_column_info(unsigned column); + +protected: + T get_column_cost_value(unsigned j, column_info * ci) const; +public: + unsigned m_total_iterations; + static_matrix* m_A; // this is the matrix of constraints + vector m_b; // the right side vector + unsigned m_first_stage_iterations; + unsigned m_second_stage_iterations; + std::unordered_map> m_constraints; + std::unordered_map*> m_map_from_var_index_to_column_info; + std::unordered_map > m_A_values; + std::unordered_map m_names_to_columns; // don't have to use it + std::unordered_map m_external_rows_to_core_solver_rows; + std::unordered_map m_core_solver_rows_to_external_rows; + std::unordered_map m_core_solver_columns_to_external_columns; + vector m_column_scale; + std::unordered_map m_name_map; + unsigned m_artificials; + unsigned m_slacks; + vector m_column_types; + vector m_costs; + vector m_x; + vector m_upper_bounds; + vector m_basis; + vector m_nbasis; + vector m_heading; + + + lp_status m_status; + + lp_settings m_settings; + lp_solver(): + m_A(nullptr), // this is the matrix of constraints + m_first_stage_iterations (0), + m_second_stage_iterations (0), + m_artificials (0), + m_slacks (0), + m_status(lp_status::UNKNOWN) + {} + + unsigned row_count() const { return this->m_A->row_count(); } + + void add_constraint(lp_relation relation, T right_side, unsigned row_index); + + void set_cost_for_column(unsigned column, T column_cost) { + get_or_create_column_info(column)->set_cost(column_cost); + } + std::string get_variable_name(unsigned j) const override; + + void set_row_column_coefficient(unsigned row, unsigned column, T const & val) { + m_A_values[row][column] = val; + } + // returns the current cost + virtual T get_current_cost() const = 0; + // do not have to call it + void give_symbolic_name_to_column(std::string name, unsigned column); + + virtual T get_column_value(unsigned column) const = 0; + + T get_column_value_by_name(std::string name) const; + + // returns -1 if not found + virtual int get_column_index_by_name(std::string name) const; + + void set_lower_bound(unsigned i, T bound) { + column_info *ci = get_or_create_column_info(i); + ci->set_lower_bound(bound); + } + + void set_upper_bound(unsigned i, T bound) { + column_info *ci = get_or_create_column_info(i); + ci->set_upper_bound(bound); + } + + void unset_lower_bound(unsigned i) { + get_or_create_column_info(i)->unset_lower_bound(); + } + + void unset_upper_bound(unsigned i) { + get_or_create_column_info(i)->unset_upper_bound(); + } + + void set_fixed_value(unsigned i, T val) { + column_info *ci = get_or_create_column_info(i); + ci->set_fixed_value(val); + } + + void unset_fixed_value(unsigned i) { + get_or_create_column_info(i)->unset_fixed(); + } + + lp_status get_status() const { + return m_status; + } + + void set_status(lp_status st) { + m_status = st; + } + + + ~lp_solver() override; + + void flip_costs(); + + virtual void find_maximal_solution() = 0; + void set_time_limit(unsigned time_limit_in_seconds) { + m_settings.time_limit = time_limit_in_seconds; + } + + +protected: + bool problem_is_empty(); + + void scale(); + + + void print_rows_scale_stats(std::ostream & out); + + void print_columns_scale_stats(std::ostream & out); + + void print_row_scale_stats(unsigned i, std::ostream & out); + + void print_column_scale_stats(unsigned j, std::ostream & out); + + void print_scale_stats(std::ostream & out); + + void get_max_abs_in_row(std::unordered_map & row_map); + + void pin_vars_down_on_row(std::unordered_map & row) { + pin_vars_on_row_with_sign(row, - numeric_traits::one()); + } + + void pin_vars_up_on_row(std::unordered_map & row) { + pin_vars_on_row_with_sign(row, numeric_traits::one()); + } + + void pin_vars_on_row_with_sign(std::unordered_map & row, T sign ); + + bool get_minimal_row_value(std::unordered_map & row, T & lower_bound); + + bool get_maximal_row_value(std::unordered_map & row, T & lower_bound); + + bool row_is_zero(std::unordered_map & row); + + bool row_e_is_obsolete(std::unordered_map & row, unsigned row_index); + + bool row_ge_is_obsolete(std::unordered_map & row, unsigned row_index); + + bool row_le_is_obsolete(std::unordered_map & row, unsigned row_index); + + // analyse possible max and min values that are derived from var boundaries + // Let us say that the we have a "ge" constraint, and the min value is equal to the rs. + // Then we know what values of the variables are. For each positive coeff of the row it has to be + // the low boundary of the var and for a negative - the upper. + + // this routing also pins the variables to the boundaries + bool row_is_obsolete(std::unordered_map & row, unsigned row_index ); + + void remove_fixed_or_zero_columns(); + + void remove_fixed_or_zero_columns_from_row(unsigned i, std::unordered_map & row); + + unsigned try_to_remove_some_rows(); + + void cleanup(); + + void map_external_rows_to_core_solver_rows(); + + void map_external_columns_to_core_solver_columns(); + + unsigned number_of_core_structurals() { + return static_cast(m_core_solver_columns_to_external_columns.size()); + } + + void restore_column_scales_to_one() { + for (unsigned i = 0; i < m_column_scale.size(); i++) m_column_scale[i] = numeric_traits::one(); + } + + void unscale(); + + void fill_A_from_A_values(); + + void fill_matrix_A_and_init_right_side(); + + void count_slacks_and_artificials(); + + void count_slacks_and_artificials_for_row(unsigned i); + + T lower_bound_shift_for_row(unsigned i); + + void fill_m_b(); + + T get_column_value_with_core_solver(unsigned column, lp_core_solver_base * core_solver) const; + void set_scaled_cost(unsigned j); + void print_statistics_on_A(std::ostream & out) { + out << "extended A[" << this->m_A->row_count() << "," << this->m_A->column_count() << "]" << std::endl; + } + +public: + lp_settings & settings() { return m_settings;} + void print_model(std::ostream & s) const { + s << "objective = " << get_current_cost() << std::endl; + s << "column values\n"; + for (auto & it : m_names_to_columns) { + s << it.first << " = " << get_column_value(it.second) << std::endl; + } + } +}; +} diff --git a/src/math/lp/lp_solver_def.h b/src/math/lp/lp_solver_def.h new file mode 100644 index 000000000..191832a24 --- /dev/null +++ b/src/math/lp/lp_solver_def.h @@ -0,0 +1,571 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once + +#include +#include +#include "util/vector.h" +#include "math/lp/lp_solver.h" +namespace lp { +template column_info * lp_solver::get_or_create_column_info(unsigned column) { + auto it = m_map_from_var_index_to_column_info.find(column); + return (it == m_map_from_var_index_to_column_info.end())? (m_map_from_var_index_to_column_info[column] = new column_info()) : it->second; +} + +template +std::string lp_solver::get_variable_name(unsigned j) const { // j here is the core solver index + if (!m_settings.print_external_var_name()) + return std::string("j")+T_to_string(j); + auto it = this->m_core_solver_columns_to_external_columns.find(j); + if (it == this->m_core_solver_columns_to_external_columns.end()) + return std::string("x")+T_to_string(j); + unsigned external_j = it->second; + auto t = this->m_map_from_var_index_to_column_info.find(external_j); + if (t == this->m_map_from_var_index_to_column_info.end()) { + return std::string("x") +T_to_string(external_j); + } + return t->second->get_name(); +} + +template T lp_solver::get_column_cost_value(unsigned j, column_info * ci) const { + if (ci->is_fixed()) { + return ci->get_cost() * ci->get_fixed_value(); + } + return ci->get_cost() * get_column_value(j); +} +template void lp_solver::add_constraint(lp_relation relation, T right_side, unsigned row_index) { + lp_assert(m_constraints.find(row_index) == m_constraints.end()); + lp_constraint cs(right_side, relation); + m_constraints[row_index] = cs; +} + +template void lp_solver::give_symbolic_name_to_column(std::string name, unsigned column) { + auto it = m_map_from_var_index_to_column_info.find(column); + column_info *ci; + if (it == m_map_from_var_index_to_column_info.end()){ + m_map_from_var_index_to_column_info[column] = ci = new column_info; + } else { + ci = it->second; + } + ci->set_name(name); + m_names_to_columns[name] = column; +} + + +template T lp_solver::get_column_value_by_name(std::string name) const { + auto it = m_names_to_columns.find(name); + if (it == m_names_to_columns.end()) { + std::stringstream s; + s << "get_column_value_by_name " << name; + throw_exception(s.str()); + } + return get_column_value(it -> second); +} + +// returns -1 if not found +template int lp_solver::get_column_index_by_name(std::string name) const { + auto t = m_names_to_columns.find(name); + if (t == m_names_to_columns.end()) { + return -1; + } + return t->second; +} + + +template lp_solver::~lp_solver(){ + delete m_A; + for (auto t : m_map_from_var_index_to_column_info) { + delete t.second; + } +} + +template void lp_solver::flip_costs() { + for (auto t : m_map_from_var_index_to_column_info) { + column_info *ci = t.second; + ci->set_cost(-ci->get_cost()); + } +} + +template bool lp_solver::problem_is_empty() { + for (auto & c : m_A_values) + if (!c.second.empty()) + return false; + return true; +} + +template void lp_solver::scale() { + if (numeric_traits::precise() || m_settings.use_scaling == false) { + m_column_scale.clear(); + m_column_scale.resize(m_A->column_count(), one_of_type()); + return; + } + + T smin = T(m_settings.scaling_minimum); + T smax = T(m_settings.scaling_maximum); + + scaler scaler(m_b, *m_A, smin, smax, m_column_scale, this->m_settings); + if (!scaler.scale()) { + unscale(); + } +} + + +template void lp_solver::print_rows_scale_stats(std::ostream & out) { + out << "rows max" << std::endl; + for (unsigned i = 0; i < m_A->row_count(); i++) { + print_row_scale_stats(i, out); + } + out << std::endl; +} + +template void lp_solver::print_columns_scale_stats(std::ostream & out) { + out << "columns max" << std::endl; + for (unsigned i = 0; i < m_A->column_count(); i++) { + print_column_scale_stats(i, out); + } + out << std::endl; +} + +template void lp_solver::print_row_scale_stats(unsigned i, std::ostream & out) { + out << "(" << std::min(m_A->get_min_abs_in_row(i), abs(m_b[i])) << " "; + out << std::max(m_A->get_max_abs_in_row(i), abs(m_b[i])) << ")"; +} + +template void lp_solver::print_column_scale_stats(unsigned j, std::ostream & out) { + out << "(" << m_A->get_min_abs_in_row(j) << " "; + out << m_A->get_max_abs_in_column(j) << ")"; +} + +template void lp_solver::print_scale_stats(std::ostream & out) { + print_rows_scale_stats(out); + print_columns_scale_stats(out); +} + +template void lp_solver::get_max_abs_in_row(std::unordered_map & row_map) { + T ret = numeric_traits::zero(); + for (auto jp : row_map) { + T ac = numeric_traits::abs(jp->second); + if (ac > ret) { + ret = ac; + } + } + return ret; +} + +template void lp_solver::pin_vars_on_row_with_sign(std::unordered_map & row, T sign ) { + for (auto t : row) { + unsigned j = t.first; + column_info * ci = m_map_from_var_index_to_column_info[j]; + T a = t.second; + if (a * sign > numeric_traits::zero()) { + lp_assert(ci->upper_bound_is_set()); + ci->set_fixed_value(ci->get_upper_bound()); + } else { + lp_assert(ci->lower_bound_is_set()); + ci->set_fixed_value(ci->get_lower_bound()); + } + } +} + +template bool lp_solver::get_minimal_row_value(std::unordered_map & row, T & lower_bound) { + lower_bound = numeric_traits::zero(); + for (auto & t : row) { + T a = t.second; + column_info * ci = m_map_from_var_index_to_column_info[t.first]; + if (a > numeric_traits::zero()) { + if (ci->lower_bound_is_set()) { + lower_bound += ci->get_lower_bound() * a; + } else { + return false; + } + } else { + if (ci->upper_bound_is_set()) { + lower_bound += ci->get_upper_bound() * a; + } else { + return false; + } + } + } + return true; +} + +template bool lp_solver::get_maximal_row_value(std::unordered_map & row, T & lower_bound) { + lower_bound = numeric_traits::zero(); + for (auto & t : row) { + T a = t.second; + column_info * ci = m_map_from_var_index_to_column_info[t.first]; + if (a < numeric_traits::zero()) { + if (ci->lower_bound_is_set()) { + lower_bound += ci->get_lower_bound() * a; + } else { + return false; + } + } else { + if (ci->upper_bound_is_set()) { + lower_bound += ci->get_upper_bound() * a; + } else { + return false; + } + } + } + return true; +} + +template bool lp_solver::row_is_zero(std::unordered_map & row) { + for (auto & t : row) { + if (!is_zero(t.second)) + return false; + } + return true; +} + +template bool lp_solver::row_e_is_obsolete(std::unordered_map & row, unsigned row_index) { + T rs = m_constraints[row_index].m_rs; + if (row_is_zero(row)) { + if (!is_zero(rs)) + m_status = lp_status::INFEASIBLE; + return true; + } + + T lower_bound; + bool lb = get_minimal_row_value(row, lower_bound); + if (lb) { + T diff = lower_bound - rs; + if (!val_is_smaller_than_eps(diff, m_settings.refactor_tolerance)){ + // lower_bound > rs + m_settings.refactor_epsilon + m_status = lp_status::INFEASIBLE; + return true; + } + if (val_is_smaller_than_eps(-diff, m_settings.refactor_tolerance)){ + pin_vars_down_on_row(row); + return true; + } + } + + T upper_bound; + bool ub = get_maximal_row_value(row, upper_bound); + if (ub) { + T diff = rs - upper_bound; + if (!val_is_smaller_than_eps(diff, m_settings.refactor_tolerance)) { + // upper_bound < rs - m_settings.refactor_tolerance + m_status = lp_status::INFEASIBLE; + return true; + } + if (val_is_smaller_than_eps(-diff, m_settings.refactor_tolerance)){ + pin_vars_up_on_row(row); + return true; + } + } + + return false; +} + +template bool lp_solver::row_ge_is_obsolete(std::unordered_map & row, unsigned row_index) { + T rs = m_constraints[row_index].m_rs; + if (row_is_zero(row)) { + if (rs > zero_of_type()) + m_status = lp_status::INFEASIBLE; + return true; + } + + T upper_bound; + if (get_maximal_row_value(row, upper_bound)) { + T diff = rs - upper_bound; + if (!val_is_smaller_than_eps(diff, m_settings.refactor_tolerance)) { + // upper_bound < rs - m_settings.refactor_tolerance + m_status = lp_status::INFEASIBLE; + return true; + } + if (val_is_smaller_than_eps(-diff, m_settings.refactor_tolerance)){ + pin_vars_up_on_row(row); + return true; + } + } + + return false; +} + +template bool lp_solver::row_le_is_obsolete(std::unordered_map & row, unsigned row_index) { + T lower_bound; + T rs = m_constraints[row_index].m_rs; + if (row_is_zero(row)) { + if (rs < zero_of_type()) + m_status = lp_status::INFEASIBLE; + return true; + } + + if (get_minimal_row_value(row, lower_bound)) { + T diff = lower_bound - rs; + if (!val_is_smaller_than_eps(diff, m_settings.refactor_tolerance)){ + // lower_bound > rs + m_settings.refactor_tolerance + m_status = lp_status::INFEASIBLE; + return true; + } + if (val_is_smaller_than_eps(-diff, m_settings.refactor_tolerance)){ + pin_vars_down_on_row(row); + return true; + } + } + + return false; +} + +// analyse possible max and min values that are derived from var boundaries +// Let us say that the we have a "ge" constraint, and the min value is equal to the rs. +// Then we know what values of the variables are. For each positive coeff of the row it has to be +// the low boundary of the var and for a negative - the upper. + +// this routing also pins the variables to the boundaries +template bool lp_solver::row_is_obsolete(std::unordered_map & row, unsigned row_index ) { + auto & constraint = m_constraints[row_index]; + switch (constraint.m_relation) { + case lp_relation::Equal: + return row_e_is_obsolete(row, row_index); + + case lp_relation::Greater_or_equal: + return row_ge_is_obsolete(row, row_index); + + case lp_relation::Less_or_equal: + return row_le_is_obsolete(row, row_index); + } + lp_unreachable(); + return false; // it is unreachable +} + +template void lp_solver::remove_fixed_or_zero_columns() { + for (auto & i_row : m_A_values) { + remove_fixed_or_zero_columns_from_row(i_row.first, i_row.second); + } +} + +template void lp_solver::remove_fixed_or_zero_columns_from_row(unsigned i, std::unordered_map & row) { + auto & constraint = m_constraints[i]; + vector removed; + for (auto & col : row) { + unsigned j = col.first; + lp_assert(m_map_from_var_index_to_column_info.find(j) != m_map_from_var_index_to_column_info.end()); + column_info * ci = m_map_from_var_index_to_column_info[j]; + if (ci->is_fixed()) { + removed.push_back(j); + T aj = col.second; + constraint.m_rs -= aj * ci->get_fixed_value(); + } else { + if (numeric_traits::is_zero(col.second)){ + removed.push_back(j); + } + } + } + + for (auto j : removed) { + row.erase(j); + } +} + +template unsigned lp_solver::try_to_remove_some_rows() { + vector rows_to_delete; + for (auto & t : m_A_values) { + if (row_is_obsolete(t.second, t.first)) { + rows_to_delete.push_back(t.first); + } + + if (m_status == lp_status::INFEASIBLE) { + return 0; + } + } + if (!rows_to_delete.empty()) { + for (unsigned k : rows_to_delete) { + m_A_values.erase(k); + } + } + remove_fixed_or_zero_columns(); + return static_cast(rows_to_delete.size()); +} + +template void lp_solver::cleanup() { + int n = 0; // number of deleted rows + int d; + while ((d = try_to_remove_some_rows()) > 0) + n += d; + + if (n == 1) { + LP_OUT(m_settings, "deleted one row" << std::endl); + } else if (n) { + LP_OUT(m_settings, "deleted " << n << " rows" << std::endl); + } +} + +template void lp_solver::map_external_rows_to_core_solver_rows() { + unsigned size = 0; + for (auto & row : m_A_values) { + m_external_rows_to_core_solver_rows[row.first] = size; + m_core_solver_rows_to_external_rows[size] = row.first; + size++; + } +} + +template void lp_solver::map_external_columns_to_core_solver_columns() { + unsigned size = 0; + for (auto & row : m_A_values) { + for (auto & col : row.second) { + if (col.second == numeric_traits::zero() || m_map_from_var_index_to_column_info[col.first]->is_fixed()) { + throw_exception("found fixed column"); + } + unsigned j = col.first; + auto column_info_it = m_map_from_var_index_to_column_info.find(j); + lp_assert(column_info_it != m_map_from_var_index_to_column_info.end()); + + auto j_column = column_info_it->second->get_column_index(); + if (!is_valid(j_column)) { // j is a newcomer + m_map_from_var_index_to_column_info[j]->set_column_index(size); + m_core_solver_columns_to_external_columns[size++] = j; + } + } + } +} + +template void lp_solver::unscale() { + delete m_A; + m_A = nullptr; + fill_A_from_A_values(); + restore_column_scales_to_one(); + fill_m_b(); +} + +template void lp_solver::fill_A_from_A_values() { + m_A = new static_matrix(static_cast(m_A_values.size()), number_of_core_structurals()); + for (auto & t : m_A_values) { + auto row_it = m_external_rows_to_core_solver_rows.find(t.first); + lp_assert(row_it != m_external_rows_to_core_solver_rows.end()); + unsigned row = row_it->second; + for (auto k : t.second) { + auto column_info_it = m_map_from_var_index_to_column_info.find(k.first); + lp_assert(column_info_it != m_map_from_var_index_to_column_info.end()); + column_info *ci = column_info_it->second; + unsigned col = ci->get_column_index(); + lp_assert(is_valid(col)); + bool col_is_flipped = m_map_from_var_index_to_column_info[k.first]->is_flipped(); + if (!col_is_flipped) { + (*m_A)(row, col) = k.second; + } else { + (*m_A)(row, col) = - k.second; + } + } + } +} + +template void lp_solver::fill_matrix_A_and_init_right_side() { + map_external_rows_to_core_solver_rows(); + map_external_columns_to_core_solver_columns(); + lp_assert(m_A == nullptr); + fill_A_from_A_values(); + m_b.resize(m_A->row_count()); +} + +template void lp_solver::count_slacks_and_artificials() { + for (int i = row_count() - 1; i >= 0; i--) { + count_slacks_and_artificials_for_row(i); + } +} + +template void lp_solver::count_slacks_and_artificials_for_row(unsigned i) { + lp_assert(this->m_constraints.find(this->m_core_solver_rows_to_external_rows[i]) != this->m_constraints.end()); + auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[i]]; + switch (constraint.m_relation) { + case Equal: + m_artificials++; + break; + case Greater_or_equal: + m_slacks++; + if (this->m_b[i] > 0) { + m_artificials++; + } + break; + case Less_or_equal: + m_slacks++; + if (this->m_b[i] < 0) { + m_artificials++; + } + break; + } +} + +template T lp_solver::lower_bound_shift_for_row(unsigned i) { + T ret = numeric_traits::zero(); + + auto row = this->m_A_values.find(i); + if (row == this->m_A_values.end()) { + throw_exception("cannot find row"); + } + for (auto col : row->second) { + ret += col.second * this->m_map_from_var_index_to_column_info[col.first]->get_shift(); + } + return ret; +} + +template void lp_solver::fill_m_b() { + for (int i = this->row_count() - 1; i >= 0; i--) { + lp_assert(this->m_constraints.find(this->m_core_solver_rows_to_external_rows[i]) != this->m_constraints.end()); + unsigned external_i = this->m_core_solver_rows_to_external_rows[i]; + auto & constraint = this->m_constraints[external_i]; + this->m_b[i] = constraint.m_rs - lower_bound_shift_for_row(external_i); + } +} + +template T lp_solver::get_column_value_with_core_solver(unsigned column, lp_core_solver_base * core_solver) const { + auto cit = this->m_map_from_var_index_to_column_info.find(column); + if (cit == this->m_map_from_var_index_to_column_info.end()) { + return numeric_traits::zero(); + } + + column_info * ci = cit->second; + + if (ci->is_fixed()) { + return ci->get_fixed_value(); + } + + unsigned cj = ci->get_column_index(); + if (cj != static_cast(-1)) { + T v = core_solver->get_var_value(cj) * this->m_column_scale[cj]; + if (ci->is_free()) { + return v; + } + if (!ci->is_flipped()) { + return v + ci->get_lower_bound(); + } + + // the flipped case when there is only upper bound + return -v + ci->get_upper_bound(); // + } + + return numeric_traits::zero(); // returns zero for out of boundary columns +} + +template void lp_solver::set_scaled_cost(unsigned j) { + // grab original costs but modify it with the column scales + lp_assert(j < this->m_column_scale.size()); + column_info * ci = this->m_map_from_var_index_to_column_info[this->m_core_solver_columns_to_external_columns[j]]; + T cost = ci->get_cost(); + if (ci->is_flipped()){ + cost *= T(-1); + } + lp_assert(ci->is_fixed() == false); + this->m_costs[j] = cost * this->m_column_scale[j]; +} +} diff --git a/src/math/lp/mps_reader.h b/src/math/lp/mps_reader.h new file mode 100644 index 000000000..8093954b1 --- /dev/null +++ b/src/math/lp/mps_reader.h @@ -0,0 +1,891 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ + +#pragma once + +// reads an MPS file representing a Mixed Integer Program +#include +#include +#include +#include "util/vector.h" +#include +#include +#include +#include +#include "math/lp/lp_primal_simplex.h" +#include "math/lp/lp_dual_simplex.h" +#include "math/lp/lar_solver.h" +#include "math/lp/lp_utils.h" +#include "math/lp/lp_solver.h" +namespace lp { +inline bool my_white_space(const char & a) { + return a == ' ' || a == '\t'; +} +inline size_t number_of_whites(const std::string & s) { + size_t i = 0; + for(;i < s.size(); i++) + if (!my_white_space(s[i])) return i; + return i; +} +inline size_t number_of_whites_from_end(const std::string & s) { + size_t ret = 0; + for(int i = static_cast(s.size()) - 1;i >= 0; i--) + if (my_white_space(s[i])) ret++;else break; + + return ret; +} + + + // trim from start +inline std::string <rim(std::string &s) { + s.erase(0, number_of_whites(s)); + return s; +} + + + + + // trim from end +inline std::string &rtrim(std::string &s) { + // s.erase(std::find_if(s.rbegin(), s.rend(), std::not1(std::ptr_fun(std::isspace))).base(), s.end()); + s.erase(s.end() - number_of_whites_from_end(s), s.end()); + return s; +} + // trim from both ends +inline std::string &trim(std::string &s) { + return ltrim(rtrim(s)); +} + +inline std::string trim(std::string const &r) { + std::string s = r; + return ltrim(rtrim(s)); +} + + +inline vector string_split(const std::string &source, const char *delimiter, bool keep_empty) { + vector results; + size_t prev = 0; + size_t next = 0; + while ((next = source.find_first_of(delimiter, prev)) != std::string::npos) { + if (keep_empty || (next - prev != 0)) { + results.push_back(source.substr(prev, next - prev)); + } + prev = next + 1; + } + if (prev < source.size()) { + results.push_back(source.substr(prev)); + } + return results; +} + +inline vector split_and_trim(const std::string &line) { + auto split = string_split(line, " \t", false); + vector ret; + for (auto s : split) { + ret.push_back(trim(s)); + } + return ret; +} + +template +class mps_reader { + enum row_type { Cost, Less_or_equal, Greater_or_equal, Equal }; + struct bound { + T m_low; + T m_upper; + bool m_low_is_set; + bool m_upper_is_set; + bool m_value_is_fixed; + T m_fixed_value; + bool m_free; + // constructor + bound() : m_low(numeric_traits::zero()), + m_low_is_set(true), + m_upper_is_set(false), + m_value_is_fixed(false), + m_free(false) {} // it seems all mps files I have seen have the default low value 0 on a variable + }; + + struct column { + std::string m_name; + bound * m_bound; + unsigned m_index; + column(const std::string &name, unsigned index): m_name(name), + m_bound(nullptr), + m_index(index) { + } + }; + + struct row { + row_type m_type; + std::string m_name; + std::unordered_map m_row_columns; + unsigned m_index; + T m_right_side; + T m_range; + row(row_type type, const std::string &name, unsigned index) : + m_type(type), + m_name(name), + m_index(index), + m_right_side(zero_of_type()), + m_range(zero_of_type()) + { + } + }; + + bool m_is_OK; + std::string m_file_name; + std::unordered_map m_rows; + std::unordered_map m_columns; + std::unordered_map m_names_to_var_index; + std::string m_line; + std::string m_name; + std::string m_cost_row_name; + std::ifstream m_file_stream; + // needed to adjust the index row + unsigned m_cost_line_count; + unsigned m_line_number; + std::ostream * m_message_stream; + + void set_m_ok_to_false() { + *m_message_stream << "setting m_is_OK to false" << std::endl; + m_is_OK = false; + } + + std::string get_string_from_position(unsigned offset) { + unsigned i = offset; + for (; i < m_line.size(); i++){ + if (m_line[i] == ' ') + break; + } + lp_assert(m_line.size() >= offset); + lp_assert(m_line.size() >> i); + lp_assert(i >= offset); + return m_line.substr(offset, i - offset); + } + + void set_boundary_for_column(unsigned col, bound * b, lp_solver * solver){ + if (b == nullptr) { + solver->set_lower_bound(col, numeric_traits::zero()); + return; + } + + if (b->m_free) { + return; + } + if (b->m_low_is_set) { + solver->set_lower_bound(col, b->m_low); + } + if (b->m_upper_is_set) { + solver->set_upper_bound(col, b->m_upper); + } + + if (b->m_value_is_fixed) { + solver->set_fixed_value(col, b->m_fixed_value); + } + } + + bool all_white_space() { + for (unsigned i = 0; i < m_line.size(); i++) { + char c = m_line[i]; + if (c != ' ' && c != '\t') { + return false; + } + } + return true; + } + + void read_line() { + while (m_is_OK) { + if (!getline(m_file_stream, m_line)) { + m_line_number++; + set_m_ok_to_false(); + *m_message_stream << "cannot read from file" << std::endl; + } + m_line_number++; + if (!m_line.empty() && m_line[0] != '*' && !all_white_space()) + break; + } + } + + void read_name() { + do { + read_line(); + if (m_line.find("NAME") != 0) { + continue; + } + m_line = m_line.substr(4); + m_name = trim(m_line); + break; + } while (m_is_OK); + } + + void read_rows() { + // look for start of the rows + read_line(); + do { + if (static_cast(m_line.find("ROWS")) >= 0) { + break; + } + } while (m_is_OK); + do { + read_line(); + if (m_line.find("COLUMNS") == 0) { + break; + } + add_row(); + } while (m_is_OK); + } + + void read_column_by_columns(const std::string & column_name, std::string column_data) { + // uph, let us try to work with columns + if (column_data.size() >= 22) { + std::string ss = column_data.substr(0, 8); + std::string row_name = trim(ss); + auto t = m_rows.find(row_name); + + if (t == m_rows.end()) { + *m_message_stream << "cannot find " << row_name << std::endl; + goto fail; + } else { + row * row = t->second; + row->m_row_columns[column_name] = numeric_traits::from_string(column_data.substr(8)); + if (column_data.size() > 24) { + column_data = column_data.substr(25); + if (column_data.size() >= 22) { + read_column_by_columns(column_name, column_data); + } + } + } + } else { + fail: + set_m_ok_to_false(); + *m_message_stream << "cannot understand this line\n" + "line = " << m_line << ", line number is " << m_line_number << std::endl; + return; + } + } + + void read_column(const std::string & column_name, const std::string & column_data){ + auto tokens = split_and_trim(column_data); + for (unsigned i = 0; i < tokens.size() - 1; i+= 2) { + auto row_name = tokens[i]; + if (row_name == "'MARKER'") return; // it is the integrality marker, no real data here + auto t = m_rows.find(row_name); + if (t == m_rows.end()) { + read_column_by_columns(column_name, column_data); + return; + } + row *r = t->second; + r->m_row_columns[column_name] = numeric_traits::from_string(tokens[i + 1]); + } + } + + void read_columns(){ + std::string column_name; + do { + read_line(); + if (m_line.find("RHS") == 0) { + break; + } + if (m_line.size() < 22) { + (*m_message_stream) << "line is too short for a column" << std::endl; + (*m_message_stream) << m_line << std::endl; + (*m_message_stream) << "line number is " << m_line_number << std::endl; + set_m_ok_to_false(); + return; + } + std::string column_name_tmp = trim(m_line.substr(4, 8)); + if (!column_name_tmp.empty()) { + column_name = column_name_tmp; + } + auto col_it = m_columns.find(column_name); + mps_reader::column * col; + if (col_it == m_columns.end()) { + col = new mps_reader::column(column_name, static_cast(m_columns.size())); + m_columns[column_name] = col; + // (*m_message_stream) << column_name << '[' << col->m_index << ']'<< std::endl; + } else { + col = col_it->second; + } + read_column(column_name, m_line.substr(14)); + } while (m_is_OK); + } + + void read_rhs() { + do { + read_line(); + if (m_line.find("BOUNDS") == 0 || m_line.find("ENDATA") == 0 || m_line.find("RANGES") == 0) { + break; + } + fill_rhs(); + } while (m_is_OK); + } + + + void fill_rhs_by_columns(std::string rhsides) { + // uph, let us try to work with columns + if (rhsides.size() >= 22) { + std::string ss = rhsides.substr(0, 8); + std::string row_name = trim(ss); + auto t = m_rows.find(row_name); + + if (t == m_rows.end()) { + (*m_message_stream) << "cannot find " << row_name << std::endl; + goto fail; + } else { + row * row = t->second; + row->m_right_side = numeric_traits::from_string(rhsides.substr(8)); + if (rhsides.size() > 24) { + rhsides = rhsides.substr(25); + if (rhsides.size() >= 22) { + fill_rhs_by_columns(rhsides); + } + } + } + } else { + fail: + set_m_ok_to_false(); + (*m_message_stream) << "cannot understand this line" << std::endl; + (*m_message_stream) << "line = " << m_line << ", line number is " << m_line_number << std::endl; + return; + } + } + + void fill_rhs() { + if (m_line.size() < 14) { + (*m_message_stream) << "line is too short" << std::endl; + (*m_message_stream) << m_line << std::endl; + (*m_message_stream) << "line number is " << m_line_number << std::endl; + set_m_ok_to_false(); + return; + } + std::string rhsides = m_line.substr(14); + vector splitted_line = split_and_trim(rhsides); + + for (unsigned i = 0; i < splitted_line.size() - 1; i += 2) { + auto t = m_rows.find(splitted_line[i]); + if (t == m_rows.end()) { + fill_rhs_by_columns(rhsides); + return; + } + row * row = t->second; + row->m_right_side = numeric_traits::from_string(splitted_line[i + 1]); + } + } + + void read_bounds() { + if (m_line.find("BOUNDS") != 0) { + return; + } + + do { + read_line(); + if (m_line[0] != ' ') { + break; + } + create_or_update_bound(); + } while (m_is_OK); + } + + void read_ranges() { + if (m_line.find("RANGES") != 0) { + return; + } + do { + read_line(); + auto sl = split_and_trim(m_line); + if (sl.size() < 2) { + break; + } + read_range(sl); + } while (m_is_OK); + } + + + void read_bound_by_columns(const std::string & colstr) { + if (colstr.size() < 14) { + (*m_message_stream) << "line is too short" << std::endl; + (*m_message_stream) << m_line << std::endl; + (*m_message_stream) << "line number is " << m_line_number << std::endl; + set_m_ok_to_false(); + return; + } + // uph, let us try to work with columns + if (colstr.size() >= 22) { + std::string ss = colstr.substr(0, 8); + std::string column_name = trim(ss); + auto t = m_columns.find(column_name); + + if (t == m_columns.end()) { + (*m_message_stream) << "cannot find " << column_name << std::endl; + goto fail; + } else { + vector bound_string; + bound_string.push_back(column_name); + if (colstr.size() > 14) { + bound_string.push_back(colstr.substr(14)); + } + mps_reader::column * col = t->second; + bound * b = col->m_bound; + if (b == nullptr) { + col->m_bound = b = new bound(); + } + update_bound(b, bound_string); + } + } else { + fail: + set_m_ok_to_false(); + (*m_message_stream) << "cannot understand this line" << std::endl; + (*m_message_stream) << "line = " << m_line << ", line number is " << m_line_number << std::endl; + return; + } + } + + void update_bound(bound * b, vector bound_string) { + /* + UP means an upper bound is applied to the variable. A bound of type LO means a lower bound is applied. A bound type of FX ("fixed") means that the variable has upper and lower bounds equal to a single value. A bound type of FR ("free") means the variable has neither lower nor upper bounds and so can take on negative values. A variation on that is MI for free negative, giving an upper bound of 0 but no lower bound. Bound type PL is for a free positive for zero to plus infinity, but as this is the normal default, it is seldom used. There are also bound types for use in MIP models - BV for binary, being 0 or 1. UI for upper integer and LI for lower integer. SC stands for semi-continuous and indicates that the variable may be zero, but if not must be equal to at least the value given. + */ + + std::string bound_type = get_string_from_position(1); + if (bound_type == "BV") { + b->m_upper_is_set = true; + b->m_upper = 1; + return; + } + + if (bound_type == "UP" || bound_type == "UI" || bound_type == "LIMITMAX") { + if (bound_string.size() <= 1){ + set_m_ok_to_false(); + return; + } + b->m_upper_is_set = true; + b->m_upper= numeric_traits::from_string(bound_string[1]); + } else if (bound_type == "LO" || bound_type == "LI") { + if (bound_string.size() <= 1){ + set_m_ok_to_false(); + return; + } + + b->m_low_is_set = true; + b->m_low = numeric_traits::from_string(bound_string[1]); + } else if (bound_type == "FR") { + b->m_free = true; + } else if (bound_type == "FX") { + if (bound_string.size() <= 1){ + set_m_ok_to_false(); + return; + } + + b->m_value_is_fixed = true; + b->m_fixed_value = numeric_traits::from_string(bound_string[1]); + } else if (bound_type == "PL") { + b->m_low_is_set = true; + b->m_low = 0; + } else if (bound_type == "MI") { + b->m_upper_is_set = true; + b->m_upper = 0; + } else { + (*m_message_stream) << "unexpected bound type " << bound_type << " at line " << m_line_number << std::endl; + set_m_ok_to_false(); + throw; + } + } + + void create_or_update_bound() { + const unsigned name_offset = 14; + lp_assert(m_line.size() >= 14); + vector bound_string = split_and_trim(m_line.substr(name_offset, m_line.size())); + + if (bound_string.empty()) { + set_m_ok_to_false(); + (*m_message_stream) << "error at line " << m_line_number << std::endl; + throw m_line; + } + + std::string name = bound_string[0]; + auto it = m_columns.find(name); + if (it == m_columns.end()){ + read_bound_by_columns(m_line.substr(14)); + return; + } + mps_reader::column * col = it->second; + bound * b = col->m_bound; + if (b == nullptr) { + col->m_bound = b = new bound(); + } + update_bound(b, bound_string); + } + + + + void read_range_by_columns(std::string rhsides) { + if (m_line.size() < 14) { + (*m_message_stream) << "line is too short" << std::endl; + (*m_message_stream) << m_line << std::endl; + (*m_message_stream) << "line number is " << m_line_number << std::endl; + set_m_ok_to_false(); + return; + } + // uph, let us try to work with columns + if (rhsides.size() >= 22) { + std::string ss = rhsides.substr(0, 8); + std::string row_name = trim(ss); + auto t = m_rows.find(row_name); + + if (t == m_rows.end()) { + (*m_message_stream) << "cannot find " << row_name << std::endl; + goto fail; + } else { + row * row = t->second; + row->m_range = numeric_traits::from_string(rhsides.substr(8)); + maybe_modify_current_row_and_add_row_for_range(row); + if (rhsides.size() > 24) { + rhsides = rhsides.substr(25); + if (rhsides.size() >= 22) { + read_range_by_columns(rhsides); + } + } + } + } else { + fail: + set_m_ok_to_false(); + (*m_message_stream) << "cannot understand this line" << std::endl; + (*m_message_stream) << "line = " << m_line << ", line number is " << m_line_number << std::endl; + return; + } + } + + + void read_range(vector & splitted_line){ + for (unsigned i = 1; i < splitted_line.size() - 1; i += 2) { + auto it = m_rows.find(splitted_line[i]); + if (it == m_rows.end()) { + read_range_by_columns(m_line.substr(14)); + return; + } + row * row = it->second; + row->m_range = numeric_traits::from_string(splitted_line[i + 1]); + maybe_modify_current_row_and_add_row_for_range(row); + } + } + + void maybe_modify_current_row_and_add_row_for_range(row * row_with_range) { + unsigned index= static_cast(m_rows.size() - m_cost_line_count); + std::string row_name = row_with_range->m_name + "_range"; + row * other_bound_range_row; + switch (row_with_range->m_type) { + case row_type::Greater_or_equal: + m_rows[row_name] = other_bound_range_row = new row(row_type::Less_or_equal, row_name, index); + other_bound_range_row->m_right_side = row_with_range->m_right_side + abs(row_with_range->m_range); + break; + case row_type::Less_or_equal: + m_rows[row_name] = other_bound_range_row = new row(row_type::Greater_or_equal, row_name, index); + other_bound_range_row->m_right_side = row_with_range->m_right_side - abs(row_with_range->m_range); + break; + case row_type::Equal: + if (row_with_range->m_range > 0) { + row_with_range->m_type = row_type::Greater_or_equal; // the existing row type change + m_rows[row_name] = other_bound_range_row = new row(row_type::Less_or_equal, row_name, index); + } else { // row->m_range < 0; + row_with_range->m_type = row_type::Less_or_equal; // the existing row type change + m_rows[row_name] = other_bound_range_row = new row(row_type::Greater_or_equal, row_name, index); + } + other_bound_range_row->m_right_side = row_with_range->m_right_side + row_with_range->m_range; + break; + default: + (*m_message_stream) << "unexpected bound type " << row_with_range->m_type << " at line " << m_line_number << std::endl; + set_m_ok_to_false(); + throw; + } + + for (auto s : row_with_range->m_row_columns) { + lp_assert(m_columns.find(s.first) != m_columns.end()); + other_bound_range_row->m_row_columns[s.first] = s.second; + } + } + + void add_row() { + if (m_line.length() < 2) { + return; + } + + m_line = trim(m_line); + char c = m_line[0]; + m_line = m_line.substr(1); + m_line = trim(m_line); + add_row(c); + } + + void add_row(char c) { + unsigned index= static_cast(m_rows.size() - m_cost_line_count); + switch (c) { + case 'E': + m_rows[m_line] = new row(row_type::Equal, m_line, index); + break; + case 'L': + m_rows[m_line] = new row(row_type::Less_or_equal, m_line, index); + break; + case 'G': + m_rows[m_line] = new row(row_type::Greater_or_equal, m_line, index); + break; + case 'N': + m_rows[m_line] = new row(row_type::Cost, m_line, index); + m_cost_row_name = m_line; + m_cost_line_count++; + break; + } + } + unsigned range_count() { + unsigned ret = 0; + for (auto s : m_rows) { + if (s.second->m_range != 0) { + ret++; + } + } + return ret; + } + + /* + If rhs is a constraint's right-hand-side value and range is the constraint's range value, then the range interval is defined according to the following table: + sense interval + G [rhs, rhs + |range|] + L [rhs - |range|, rhs] + E [rhs, rhs + |range|] if range > 0, + [rhs - |range|, rhs] if range < 0 + where |range| is range's absolute value. + */ + + lp_relation get_relation_from_row(row_type rt) { + switch (rt) { + case mps_reader::Less_or_equal: return lp_relation::Less_or_equal; + case mps_reader::Greater_or_equal: return lp_relation::Greater_or_equal; + case mps_reader::Equal: return lp_relation::Equal; + default: + (*m_message_stream) << "Unexpected rt " << rt << std::endl; + set_m_ok_to_false(); + throw; + } + } + + unsigned solver_row_count() { + return m_rows.size() - m_cost_line_count + range_count(); + } + + void fill_solver_on_row(row * row, lp_solver *solver) { + if (row->m_name != m_cost_row_name) { + solver->add_constraint(get_relation_from_row(row->m_type), row->m_right_side, row->m_index); + for (auto s : row->m_row_columns) { + lp_assert(m_columns.find(s.first) != m_columns.end()); + solver->set_row_column_coefficient(row->m_index, m_columns[s.first]->m_index, s.second); + } + } else { + set_solver_cost(row, solver); + } + } + + T abs(T & t) { return t < numeric_traits::zero() ? -t: t; } + + void fill_solver_on_rows(lp_solver * solver) { + for (auto row_it : m_rows) { + fill_solver_on_row(row_it.second, solver); + } + } + + + void fill_solver_on_columns(lp_solver * solver){ + for (auto s : m_columns) { + mps_reader::column * col = s.second; + unsigned index = col->m_index; + set_boundary_for_column(index, col->m_bound, solver); + // optional call + solver->give_symbolic_name_to_column(col->m_name, col->m_index); + } + } + + void fill_solver(lp_solver *solver) { + fill_solver_on_rows(solver); + fill_solver_on_columns(solver); + } + + void set_solver_cost(row * row, lp_solver *solver) { + for (auto s : row->m_row_columns) { + std::string name = s.first; + lp_assert(m_columns.find(name) != m_columns.end()); + mps_reader::column * col = m_columns[name]; + solver->set_cost_for_column(col->m_index, s.second); + } + } + +public: + + void set_message_stream(std::ostream * o) { + lp_assert(o != nullptr); + m_message_stream = o; + } + vector column_names() { + vector v; + for (auto s : m_columns) { + v.push_back(s.first); + } + return v; + } + + ~mps_reader() { + for (auto s : m_rows) { + delete s.second; + } + for (auto s : m_columns) { + auto col = s.second; + delete col->m_bound; + delete col; + } + } + + mps_reader(const std::string & file_name): + m_is_OK(true), + m_file_name(file_name), + m_file_stream(file_name), + m_cost_line_count(0), + m_line_number(0), + m_message_stream(& std::cout) {} + void read() { + if (!m_file_stream.is_open()){ + set_m_ok_to_false(); + return; + } + + read_name(); + read_rows(); + read_columns(); + read_rhs(); + if (m_line.find("BOUNDS") == 0) { + read_bounds(); + read_ranges(); + } else if (m_line.find("RANGES") == 0) { + read_ranges(); + read_bounds(); + } + } + + bool is_ok() { + return m_is_OK; + } + + lp_solver * create_solver(bool dual) { + lp_solver * solver = dual? (lp_solver*)new lp_dual_simplex() : new lp_primal_simplex(); + fill_solver(solver); + return solver; + } + + lconstraint_kind get_lar_relation_from_row(row_type rt) { + switch (rt) { + case Less_or_equal: return LE; + case Greater_or_equal: return GE; + case Equal: return EQ; + default: + (*m_message_stream) << "Unexpected rt " << rt << std::endl; + set_m_ok_to_false(); + throw; + } + } + + unsigned get_var_index(std::string s) { + auto it = m_names_to_var_index.find(s); + if (it != m_names_to_var_index.end()) + return it->second; + unsigned ret = static_cast(m_names_to_var_index.size()); + m_names_to_var_index[s] = ret; + return ret; + } + + void fill_lar_solver_on_row(row * row, lar_solver *solver, int row_index) { + if (row->m_name != m_cost_row_name) { + auto kind = get_lar_relation_from_row(row->m_type); + vector> ls; + for (auto s : row->m_row_columns) { + var_index i = solver->add_var(get_var_index(s.first), false); + ls.push_back(std::make_pair(s.second, i)); + } + unsigned j = solver->add_term(ls, row_index); + solver->add_var_bound(j, kind, row->m_right_side); + } else { + // ignore the cost row + } + } + + + void fill_lar_solver_on_rows(lar_solver * solver) { + int row_index = 0; + for (auto row_it : m_rows) { + fill_lar_solver_on_row(row_it.second, solver, row_index++); + } + } + + void create_low_constraint_for_var(column* col, bound * b, lar_solver *solver) { + var_index i = solver->add_var(col->m_index, false); + solver->add_var_bound(i, GE, b->m_low); + } + + void create_upper_constraint_for_var(column* col, bound * b, lar_solver *solver) { + var_index i = solver->add_var(col->m_index, false); + solver->add_var_bound(i, LE, b->m_upper); + } + + void create_equality_contraint_for_var(column* col, bound * b, lar_solver *solver) { + var_index i = solver->add_var(col->m_index, false); + solver->add_var_bound(i, LE, b->m_fixed_value); + solver->add_var_bound(i, GE, b->m_fixed_value); + } + + void fill_lar_solver_on_columns(lar_solver * solver) { + for (auto s : m_columns) { + mps_reader::column * col = s.second; + solver->add_var(col->m_index, false); + auto b = col->m_bound; + if (b == nullptr) return; + + if (b->m_free) continue; + + if (b->m_low_is_set) { + create_low_constraint_for_var(col, b, solver); + } + if (b->m_upper_is_set) { + create_upper_constraint_for_var(col, b, solver); + } + if (b->m_value_is_fixed) { + create_equality_contraint_for_var(col, b, solver); + } + } + } + + + void fill_lar_solver(lar_solver * solver) { + fill_lar_solver_on_columns(solver); + fill_lar_solver_on_rows(solver); + } + + lar_solver * create_lar_solver() { + lar_solver * solver = new lar_solver(); + fill_lar_solver(solver); + return solver; + } +}; +} diff --git a/src/math/lp/static_matrix.cpp b/src/math/lp/static_matrix.cpp index cde96277b..571e9b1d0 100644 --- a/src/math/lp/static_matrix.cpp +++ b/src/math/lp/static_matrix.cpp @@ -24,6 +24,7 @@ Revision History: #include "math/lp/static_matrix_def.h" #include "math/lp/lp_core_solver_base.h" #include "math/lp/lp_dual_core_solver.h" +#include "math/lp/lp_dual_simplex.h" #include "math/lp/lp_primal_core_solver.h" #include "math/lp/scaler.h" #include "math/lp/lar_solver.h" diff --git a/src/sat/smt/arith_sls.h b/src/sat/smt/arith_sls.h index 09a56c84e..af3a46234 100644 --- a/src/sat/smt/arith_sls.h +++ b/src/sat/smt/arith_sls.h @@ -19,6 +19,9 @@ Author: #include "util/obj_pair_set.h" #include "ast/ast_trail.h" #include "ast/arith_decl_plugin.h" +#include "math/lp/lp_solver.h" +#include "math/lp/lp_primal_simplex.h" +#include "math/lp/lp_dual_simplex.h" #include "math/lp/indexed_value.h" #include "math/lp/lar_solver.h" #include "math/lp/nla_solver.h" diff --git a/src/sat/smt/arith_solver.h b/src/sat/smt/arith_solver.h index 68d5f8025..9ae671c16 100644 --- a/src/sat/smt/arith_solver.h +++ b/src/sat/smt/arith_solver.h @@ -19,7 +19,9 @@ Author: #include "util/obj_pair_set.h" #include "ast/ast_trail.h" #include "ast/arith_decl_plugin.h" - +#include "math/lp/lp_solver.h" +#include "math/lp/lp_primal_simplex.h" +#include "math/lp/lp_dual_simplex.h" #include "math/lp/indexed_value.h" #include "math/lp/lar_solver.h" #include "math/lp/nla_solver.h" diff --git a/src/shell/CMakeLists.txt b/src/shell/CMakeLists.txt index c0e9c8505..d9b74f162 100644 --- a/src/shell/CMakeLists.txt +++ b/src/shell/CMakeLists.txt @@ -28,6 +28,7 @@ add_executable(shell opt_frontend.cpp smtlib_frontend.cpp z3_log_frontend.cpp + lp_frontend.cpp # FIXME: shell should really link against libz3 but it can't due to requiring # use of some hidden symbols. Also libz3 has the ``api_dll`` component which # we don't want (I think). diff --git a/src/shell/lp_frontend.cpp b/src/shell/lp_frontend.cpp new file mode 100644 index 000000000..8d6425533 --- /dev/null +++ b/src/shell/lp_frontend.cpp @@ -0,0 +1,109 @@ +/*++ +Copyright (c) 2016 Microsoft Corporation + +Author: + + Lev Nachmanson 2016-10-27 + +--*/ + +#include "math/lp/lp_settings.h" +#include "math/lp/mps_reader.h" +#include "util/timeout.h" +#include "util/cancel_eh.h" +#include "util/scoped_timer.h" +#include "util/rlimit.h" +#include "util/gparams.h" +#include "util/mutex.h" +#include +#include +#include "smt/params/smt_params_helper.hpp" + +namespace { +static mutex *display_stats_mux = new mutex; + +static lp::lp_solver* g_solver = nullptr; + +static void display_statistics() { + lock_guard lock(*display_stats_mux); + if (g_solver && g_solver->settings().print_statistics) { + // TBD display relevant information about statistics + } +} + +static void STD_CALL on_ctrl_c(int) { + signal (SIGINT, SIG_DFL); + display_statistics(); + raise(SIGINT); +} + +static void on_timeout() { + display_statistics(); + _Exit(0); +} + +struct front_end_resource_limit : public lp::lp_resource_limit { + reslimit& m_reslim; + + front_end_resource_limit(reslimit& lim): + m_reslim(lim) + {} + + bool get_cancel_flag() override { return !m_reslim.inc(); } +}; + +void run_solver(smt_params_helper & params, char const * mps_file_name) { + + reslimit rlim; + unsigned timeout = gparams::get_ref().get_uint("timeout", 0); + unsigned rlimit = gparams::get_ref().get_uint("rlimit", 0); + front_end_resource_limit lp_limit(rlim); + + scoped_rlimit _rlimit(rlim, rlimit); + cancel_eh eh(rlim); + scoped_timer timer(timeout, &eh); + + std::string fn(mps_file_name); + lp::mps_reader reader(fn); + reader.set_message_stream(&std::cout); // can be redirected + reader.read(); + if (!reader.is_ok()) { + std::cerr << "cannot process " << mps_file_name << std::endl; + return; + } + lp::lp_solver * solver = reader.create_solver(false); // false - to create the primal solver + solver->settings().set_resource_limit(lp_limit); + g_solver = solver; + if (params.arith_min()) { + solver->flip_costs(); + } + solver->settings().set_message_ostream(&std::cout); + solver->settings().report_frequency = params.arith_rep_freq(); + solver->settings().print_statistics = params.arith_print_stats(); + solver->settings().set_simplex_strategy(lp:: simplex_strategy_enum::lu); + + solver->find_maximal_solution(); + + *(solver->settings().get_message_ostream()) << "status is " << lp_status_to_string(solver->get_status()) << std::endl; + if (solver->get_status() == lp::lp_status::OPTIMAL) { + if (params.arith_min()) { + solver->flip_costs(); + } + solver->print_model(std::cout); + } + + display_statistics(); + g_solver = nullptr; + delete solver; +} +} + +unsigned read_mps_file(char const * mps_file_name) { + signal(SIGINT, on_ctrl_c); + register_on_timeout_proc(on_timeout); + smt_params_helper p; + param_descrs r; + p.collect_param_descrs(r); + run_solver(p, mps_file_name); + return 0; +} diff --git a/src/shell/lp_frontend.h b/src/shell/lp_frontend.h new file mode 100644 index 000000000..b24be811f --- /dev/null +++ b/src/shell/lp_frontend.h @@ -0,0 +1,7 @@ +/* + Copyright (c) 2013 Microsoft Corporation. All rights reserved. + + Author: Lev Nachmanson +*/ +#pragma once +unsigned read_mps_file(char const * mps_file_name); diff --git a/src/shell/main.cpp b/src/shell/main.cpp index 26325d3d0..4c26d91d9 100644 --- a/src/shell/main.cpp +++ b/src/shell/main.cpp @@ -37,13 +37,14 @@ Revision History: #include "util/gparams.h" #include "util/env_params.h" #include "util/file_path.h" +#include "shell/lp_frontend.h" #include "shell/drat_frontend.h" #if defined( _WINDOWS ) && defined( __MINGW32__ ) && ( defined( __GNUG__ ) || defined( __clang__ ) ) #include #endif -typedef enum { IN_UNSPECIFIED, IN_SMTLIB_2, IN_DATALOG, IN_DIMACS, IN_WCNF, IN_OPB, IN_LP, IN_Z3_LOG, IN_DRAT } input_kind; +typedef enum { IN_UNSPECIFIED, IN_SMTLIB_2, IN_DATALOG, IN_DIMACS, IN_WCNF, IN_OPB, IN_LP, IN_Z3_LOG, IN_MPS, IN_DRAT } input_kind; static char const * g_input_file = nullptr; static char const * g_drat_input_file = nullptr; @@ -376,6 +377,10 @@ int STD_CALL main(int argc, char ** argv) { else if (strcmp(ext, "smt2") == 0) { g_input_kind = IN_SMTLIB_2; } + else if (strcmp(ext, "mps") == 0 || strcmp(ext, "sif") == 0 || + strcmp(ext, "MPS") == 0 || strcmp(ext, "SIF") == 0) { + g_input_kind = IN_MPS; + } } } switch (g_input_kind) { @@ -401,6 +406,9 @@ int STD_CALL main(int argc, char ** argv) { case IN_Z3_LOG: replay_z3_log(g_input_file); break; + case IN_MPS: + return_value = read_mps_file(g_input_file); + break; case IN_DRAT: return_value = read_drat(g_drat_input_file); break; diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 2aa988282..d61910ff2 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -19,6 +19,9 @@ --*/ #include "util/stopwatch.h" +#include "math/lp/lp_solver.h" +#include "math/lp/lp_primal_simplex.h" +#include "math/lp/lp_dual_simplex.h" #include "math/lp/indexed_value.h" #include "math/lp/lar_solver.h" #include "math/lp/nla_solver.h" diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 78abf1a6f..c82cdd0a4 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -33,6 +33,8 @@ #include #include #include "math/lp/lp_utils.h" +#include "math/lp/lp_primal_simplex.h" +#include "math/lp/mps_reader.h" #include "test/lp/smt_reader.h" #include "math/lp/binary_heap_priority_queue.h" #include "test/lp/argument_parser.h" @@ -57,71 +59,6 @@ #include "math/lp/cross_nested.h" #include "math/lp/int_cube.h" #include "math/lp/emonics.h" - -bool my_white_space(const char & a) { - return a == ' ' || a == '\t'; -} -size_t number_of_whites(const std::string & s) { - size_t i = 0; - for(;i < s.size(); i++) - if (!my_white_space(s[i])) return i; - return i; -} -size_t number_of_whites_from_end(const std::string & s) { - size_t ret = 0; - for(int i = static_cast(s.size()) - 1;i >= 0; i--) - if (my_white_space(s[i])) ret++;else break; - - return ret; -} - - -std::string <rim(std::string &s) { - s.erase(0, number_of_whites(s)); - return s; -} - - - - - // trim from end -inline std::string &rtrim(std::string &s) { - // s.erase(std::find_if(s.rbegin(), s.rend(), std::not1(std::ptr_fun(std::isspace))).base(), s.end()); - s.erase(s.end() - number_of_whites_from_end(s), s.end()); - return s; -} - // trim from both ends -inline std::string &trim(std::string &s) { - return ltrim(rtrim(s)); -} - - -vector string_split(const std::string &source, const char *delimiter, bool keep_empty) { - vector results; - size_t prev = 0; - size_t next = 0; - while ((next = source.find_first_of(delimiter, prev)) != std::string::npos) { - if (keep_empty || (next - prev != 0)) { - results.push_back(source.substr(prev, next - prev)); - } - prev = next + 1; - } - if (prev < source.size()) { - results.push_back(source.substr(prev)); - } - return results; -} - -vector split_and_trim(const std::string &line) { - auto split = string_split(line, " \t", false); - vector ret; - for (auto s : split) { - ret.push_back(trim(s)); - } - return ret; -} - - namespace nla { void test_horner(); void test_monics(); @@ -1458,15 +1395,169 @@ void update_settings(argument_parser & args_parser, lp_settings& settings) { } } +template +void setup_solver(unsigned time_limit, bool look_for_min, argument_parser & args_parser, lp_solver * solver) { + if (time_limit > 0) + solver->set_time_limit(time_limit); + + if (look_for_min) + solver->flip_costs(); + + update_settings(args_parser, solver->settings()); +} bool values_are_one_percent_close(double a, double b); +void print_x(mps_reader & reader, lp_solver * solver) { + for (const auto & name : reader.column_names()) { + std::cout << name << "=" << solver->get_column_value_by_name(name) << ' '; + } + std::cout << std::endl; +} + +void compare_solutions(mps_reader & reader, lp_solver * solver, lp_solver * solver0) { + for (const auto & name : reader.column_names()) { + double a = solver->get_column_value_by_name(name); + double b = solver0->get_column_value_by_name(name); + if (!values_are_one_percent_close(a, b)) { + std::cout << "different values for " << name << ":" << a << " and " << b << std::endl; + } + } +} +void solve_mps_double(std::string file_name, bool look_for_min, unsigned time_limit, bool dual, bool compare_with_primal, argument_parser & args_parser) { + mps_reader reader(file_name); + reader.read(); + if (!reader.is_ok()) { + std::cout << "cannot process " << file_name << std::endl; + return; + } + + lp_solver * solver = reader.create_solver(dual); + setup_solver(time_limit, look_for_min, args_parser, solver); + stopwatch sw; + sw.start(); + if (dual) { + std::cout << "solving for dual" << std::endl; + } + solver->find_maximal_solution(); + sw.stop(); + double span = sw.get_seconds(); + std::cout << "Status: " << lp_status_to_string(solver->get_status()) << std::endl; + if (solver->get_status() == lp_status::OPTIMAL) { + if (reader.column_names().size() < 20) { + print_x(reader, solver); + } + double cost = solver->get_current_cost(); + if (look_for_min) { + cost = -cost; + } + std::cout << "cost = " << cost << std::endl; + } + std::cout << "processed in " << span / 1000.0 << " seconds, running for " << solver->m_total_iterations << " iterations" << " one iter for " << (double)span/solver->m_total_iterations << " ms" << std::endl; + if (compare_with_primal) { + auto * primal_solver = reader.create_solver(false); + setup_solver(time_limit, look_for_min, args_parser, primal_solver); + primal_solver->find_maximal_solution(); + if (solver->get_status() != primal_solver->get_status()) { + std::cout << "statuses are different: dual " << lp_status_to_string(solver->get_status()) << " primal = " << lp_status_to_string(primal_solver->get_status()) << std::endl; + } else { + if (solver->get_status() == lp_status::OPTIMAL) { + double cost = solver->get_current_cost(); + if (look_for_min) { + cost = -cost; + } + double primal_cost = primal_solver->get_current_cost(); + if (look_for_min) { + primal_cost = -primal_cost; + } + std::cout << "primal cost = " << primal_cost << std::endl; + if (!values_are_one_percent_close(cost, primal_cost)) { + compare_solutions(reader, primal_solver, solver); + print_x(reader, primal_solver); + std::cout << "dual cost is " << cost << ", but primal cost is " << primal_cost << std::endl; + lp_assert(false); + } + } + } + delete primal_solver; + } + delete solver; +} +void solve_mps_rational(std::string file_name, bool look_for_min, unsigned time_limit, bool dual, argument_parser & args_parser) { + mps_reader reader(file_name); + reader.read(); + if (reader.is_ok()) { + auto * solver = reader.create_solver(dual); + setup_solver(time_limit, look_for_min, args_parser, solver); + stopwatch sw; + sw.start(); + solver->find_maximal_solution(); + std::cout << "Status: " << lp_status_to_string(solver->get_status()) << std::endl; + + if (solver->get_status() == lp_status::OPTIMAL) { + // for (auto name: reader.column_names()) { + // std::cout << name << "=" << solver->get_column_value_by_name(name) << ' '; + // } + lp::mpq cost = solver->get_current_cost(); + if (look_for_min) { + cost = -cost; + } + std::cout << "cost = " << cost.get_double() << std::endl; + } + std::cout << "processed in " << sw.get_current_seconds() / 1000.0 << " seconds, running for " << solver->m_total_iterations << " iterations" << std::endl; + delete solver; + } else { + std::cout << "cannot process " << file_name << std::endl; + } +} void get_time_limit_and_max_iters_from_parser(argument_parser & args_parser, unsigned & time_limit); // forward definition +void solve_mps(std::string file_name, bool look_for_min, unsigned time_limit, bool solve_for_rational, bool dual, bool compare_with_primal, argument_parser & args_parser) { + if (!solve_for_rational) { + std::cout << "solving " << file_name << std::endl; + solve_mps_double(file_name, look_for_min, time_limit, dual, compare_with_primal, args_parser); + } + else { + std::cout << "solving " << file_name << " in rationals " << std::endl; + solve_mps_rational(file_name, look_for_min, time_limit, dual, args_parser); + } +} +void solve_mps(std::string file_name, argument_parser & args_parser) { + bool look_for_min = args_parser.option_is_used("--min"); + unsigned time_limit; + bool solve_for_rational = args_parser.option_is_used("--mpq"); + bool dual = args_parser.option_is_used("--dual"); + bool compare_with_primal = args_parser.option_is_used("--compare_with_primal"); + get_time_limit_and_max_iters_from_parser(args_parser, time_limit); + solve_mps(file_name, look_for_min, time_limit, solve_for_rational, dual, compare_with_primal, args_parser); +} + +void solve_mps_in_rational(std::string file_name, bool dual, argument_parser & /*args_parser*/) { + std::cout << "solving " << file_name << std::endl; + + mps_reader reader(file_name); + reader.read(); + if (reader.is_ok()) { + auto * solver = reader.create_solver(dual); + solver->find_maximal_solution(); + std::cout << "status is " << lp_status_to_string(solver->get_status()) << std::endl; + if (solver->get_status() == lp_status::OPTIMAL) { + if (reader.column_names().size() < 20) { + for (const auto & name : reader.column_names()) { + std::cout << name << "=" << solver->get_column_value_by_name(name).get_double() << ' '; + } + } + std::cout << std::endl << "cost = " << numeric_traits::get_double(solver->get_current_cost()) << std::endl; + } + delete solver; + } else { + std::cout << "cannot process " << file_name << std::endl; + } +} void test_upair_queue() { int n = 10; @@ -1535,6 +1626,53 @@ void test_binary_priority_queue() { std::cout << " done" << std::endl; } +bool solution_is_feasible(std::string file_name, const std::unordered_map & solution) { + mps_reader reader(file_name); + reader.read(); + if (reader.is_ok()) { + lp_primal_simplex * solver = static_cast *>(reader.create_solver(false)); + return solver->solution_is_feasible(solution); + } + return false; +} + + +void solve_mps_with_known_solution(std::string file_name, std::unordered_map * solution, lp_status status, bool dual) { + std::cout << "solving " << file_name << std::endl; + mps_reader reader(file_name); + reader.read(); + if (reader.is_ok()) { + auto * solver = reader.create_solver(dual); + solver->find_maximal_solution(); + std::cout << "status is " << lp_status_to_string(solver->get_status()) << std::endl; + if (status != solver->get_status()){ + std::cout << "status should be " << lp_status_to_string(status) << std::endl; + lp_assert(status == solver->get_status()); + throw "status is wrong"; + } + if (solver->get_status() == lp_status::OPTIMAL) { + std::cout << "cost = " << solver->get_current_cost() << std::endl; + if (solution != nullptr) { + for (auto it : *solution) { + if (fabs(it.second - solver->get_column_value_by_name(it.first)) >= 0.000001) { + std::cout << "expected:" << it.first << "=" << + it.second <<", got " << solver->get_column_value_by_name(it.first) << std::endl; + } + lp_assert(fabs(it.second - solver->get_column_value_by_name(it.first)) < 0.000001); + } + } + if (reader.column_names().size() < 20) { + for (const auto & name : reader.column_names()) { + std::cout << name << "=" << solver->get_column_value_by_name(name) << ' '; + } + std::cout << std::endl; + } + } + delete solver; + } else { + std::cout << "cannot process " << file_name << std::endl; + } +} int get_random_rows() { return 5 + my_random() % 2; @@ -1548,6 +1686,55 @@ int get_random_int() { return -1 + my_random() % 2; // (1.0 + RAND_MAX); } +void add_random_row(lp_primal_simplex * solver, int cols, int row) { + solver->add_constraint(lp_relation::Greater_or_equal, 1, row); + for (int i = 0; i < cols; i++) { + solver->set_row_column_coefficient(row, i, get_random_int()); + } +} + +void add_random_cost(lp_primal_simplex * solver, int cols) { + for (int i = 0; i < cols; i++) { + solver->set_cost_for_column(i, get_random_int()); + } +} + +lp_primal_simplex * generate_random_solver() { + int rows = get_random_rows(); + int cols = get_random_columns(); + auto * solver = new lp_primal_simplex(); + for (int i = 0; i < rows; i++) { + add_random_row(solver, cols, i); + } + add_random_cost(solver, cols); + return solver; +} + + + +void random_test_on_i(unsigned i) { + if (i % 1000 == 0) { + std::cout << "."; + } + srand(i); + auto *solver = generate_random_solver(); + solver->find_maximal_solution(); + // std::cout << lp_status_to_string(solver->get_status()) << std::endl; + delete solver; +} + +void random_test() { + for (unsigned i = 0; i < std::numeric_limits::max(); i++) { + try { + random_test_on_i(i); + } + catch (const char * error) { + std::cout << "i = " << i << ", throwing at ' " << error << "'" << std::endl; + break; + } + } +} + #ifndef _WINDOWS void fill_file_names(vector &file_names, std::set & minimums) { char *home_dir = getenv("HOME"); @@ -1709,9 +1896,140 @@ void find_dir_and_file_name(std::string a, std::string & dir, std::string& fn) { // std::cout << "fn = " << fn << std::endl; } +void process_test_file(std::string test_dir, std::string test_file_name, argument_parser & args_parser, std::string out_dir, unsigned max_iters, unsigned time_limit, unsigned & successes, unsigned & failures, unsigned & inconclusives); +void solve_some_mps(argument_parser & args_parser) { + unsigned max_iters = UINT_MAX, time_limit = UINT_MAX; + get_time_limit_and_max_iters_from_parser(args_parser, time_limit); + unsigned successes = 0; + unsigned failures = 0; + unsigned inconclusives = 0; + std::set minimums; + vector file_names; + fill_file_names(file_names, minimums); + bool solve_for_rational = args_parser.option_is_used("--mpq"); + bool dual = args_parser.option_is_used("--dual"); + bool compare_with_primal = args_parser.option_is_used("--compare_with_primal"); + bool compare_with_glpk = args_parser.option_is_used("--compare_with_glpk"); + if (compare_with_glpk) { + std::string out_dir = args_parser.get_option_value("--out_dir"); + if (out_dir.size() == 0) { + out_dir = "/tmp/test"; + } + test_out_dir(out_dir); + for (auto& a : file_names) { + try { + std::string file_dir; + std::string file_name; + find_dir_and_file_name(a, file_dir, file_name); + process_test_file(file_dir, file_name, args_parser, out_dir, max_iters, time_limit, successes, failures, inconclusives); + } + catch(const char *s){ + std::cout<< "exception: "<< s << std::endl; + } + } + std::cout << "comparing with glpk: successes " << successes << ", failures " << failures << ", inconclusives " << inconclusives << std::endl; + return; + } + if (!solve_for_rational) { + solve_mps(file_names[6], false, time_limit, false, dual, compare_with_primal, args_parser); + solve_mps_with_known_solution(file_names[3], nullptr, lp_status::INFEASIBLE, dual); // chvatal: 135(d) + std::unordered_map sol; + sol["X1"] = 0; + sol["X2"] = 6; + sol["X3"] = 0; + sol["X4"] = 15; + sol["X5"] = 2; + sol["X6"] = 1; + sol["X7"] = 1; + sol["X8"] = 0; + solve_mps_with_known_solution(file_names[9], &sol, lp_status::OPTIMAL, dual); + solve_mps_with_known_solution(file_names[0], &sol, lp_status::OPTIMAL, dual); + sol.clear(); + sol["X1"] = 25.0/14.0; + // sol["X2"] = 0; + // sol["X3"] = 0; + // sol["X4"] = 0; + // sol["X5"] = 0; + // sol["X6"] = 0; + // sol["X7"] = 9.0/14.0; + solve_mps_with_known_solution(file_names[5], &sol, lp_status::OPTIMAL, dual); // chvatal: 135(e) + solve_mps_with_known_solution(file_names[4], &sol, lp_status::OPTIMAL, dual); // chvatal: 135(e) + solve_mps_with_known_solution(file_names[2], nullptr, lp_status::UNBOUNDED, dual); // chvatal: 135(c) + solve_mps_with_known_solution(file_names[1], nullptr, lp_status::UNBOUNDED, dual); // chvatal: 135(b) + solve_mps(file_names[8], false, time_limit, false, dual, compare_with_primal, args_parser); + // return; + for (auto& s : file_names) { + try { + solve_mps(s, minimums.find(s) != minimums.end(), time_limit, false, dual, compare_with_primal, args_parser); + } + catch(const char *s){ + std::cout<< "exception: "<< s << std::endl; + } + } + } else { + // unsigned i = 0; + for (auto& s : file_names) { + // if (i++ > 9) return; + try { + solve_mps_in_rational(s, dual, args_parser); + } + catch(const char *s){ + std::cout<< "exception: "<< s << std::endl; + } + } + } +} #endif +void solve_rational() { + lp_primal_simplex solver; + solver.add_constraint(lp_relation::Equal, lp::mpq(7), 0); + solver.add_constraint(lp_relation::Equal, lp::mpq(-3), 1); + + // setting the cost + int cost[] = {-3, -1, -1, 2, -1, 1, 1, -4}; + std::string var_names[8] = {"x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"}; + + for (unsigned i = 0; i < 8; i++) { + solver.set_cost_for_column(i, lp::mpq(cost[i])); + solver.give_symbolic_name_to_column(var_names[i], i); + } + + int row0[] = {1, 0, 3, 1, -5, -2 , 4, -6}; + for (unsigned i = 0; i < 8; i++) { + solver.set_row_column_coefficient(0, i, lp::mpq(row0[i])); + } + + int row1[] = {0, 1, -2, -1, 4, 1, -3, 5}; + for (unsigned i = 0; i < 8; i++) { + solver.set_row_column_coefficient(1, i, lp::mpq(row1[i])); + } + + int bounds[] = {8, 6, 4, 15, 2, 10, 10, 3}; + for (unsigned i = 0; i < 8; i++) { + solver.set_lower_bound(i, lp::mpq(0)); + solver.set_upper_bound(i, lp::mpq(bounds[i])); + } + + std::unordered_map expected_sol; + expected_sol["x1"] = lp::mpq(0); + expected_sol["x2"] = lp::mpq(6); + expected_sol["x3"] = lp::mpq(0); + expected_sol["x4"] = lp::mpq(15); + expected_sol["x5"] = lp::mpq(2); + expected_sol["x6"] = lp::mpq(1); + expected_sol["x7"] = lp::mpq(1); + expected_sol["x8"] = lp::mpq(0); + solver.find_maximal_solution(); + lp_assert(solver.get_status() == lp_status::OPTIMAL); +#ifdef Z3DEBUG + for (const auto & it : expected_sol) { + (void)it; + lp_assert(it.second == solver.get_column_value_by_name(it.first)); + } +#endif +} std::string read_line(bool & end, std::ifstream & file) { @@ -1729,6 +2047,49 @@ bool contains(std::string const & s, char const * pattern) { } +std::unordered_map * get_solution_from_glpsol_output(std::string & file_name) { + std::ifstream file(file_name); + if (!file.is_open()){ + std::cerr << "cannot open " << file_name << std::endl; + return nullptr; + } + std::string s; + bool end; + do { + s = read_line(end, file); + if (end) { + std::cerr << "unexpected file end " << file_name << std::endl; + return nullptr; + } + if (contains(s, "Column name")){ + break; + } + } while (true); + + read_line(end, file); + if (end) { + std::cerr << "unexpected file end " << file_name << std::endl; + return nullptr; + } + + auto ret = new std::unordered_map(); + + do { + s = read_line(end, file); + if (end) { + std::cerr << "unexpected file end " << file_name << std::endl; + return nullptr; + } + auto split = string_split(s, " \t", false); + if (split.empty()) { + return ret; + } + + lp_assert(split.size() > 3); + (*ret)[split[1]] = atof(split[3].c_str()); + } while (true); +} + void test_init_U() { @@ -1828,6 +2189,7 @@ void setup_args_parser(argument_parser & parser) { parser.add_option_with_help_string("--test_lp_0", "solve a small lp"); parser.add_option_with_help_string("--solve_some_mps", "solves a list of mps problems"); parser.add_option_with_after_string_with_help("--test_file_directory", "loads files from the directory for testing"); + parser.add_option_with_help_string("--compare_with_glpk", "compares the results by running glpsol"); parser.add_option_with_after_string_with_help("--out_dir", "setting the output directory for tests, if not set /tmp is used"); parser.add_option_with_help_string("--dual", "using the dual simplex solver"); parser.add_option_with_help_string("--compare_with_primal", "using the primal simplex solver for comparison"); @@ -1998,9 +2360,237 @@ std::string get_status(std::string file_name) { throw 0; } +// returns true if the costs should be compared too +bool compare_statuses(std::string glpk_out_file_name, std::string lp_out_file_name, unsigned & successes, unsigned & failures) { + std::string glpk_status = get_status(glpk_out_file_name); + std::string lp_tst_status = get_status(lp_out_file_name); + + if (glpk_status != lp_tst_status) { + if (glpk_status == "UNDEFINED" && (lp_tst_status == "UNBOUNDED" || lp_tst_status == "INFEASIBLE")) { + successes++; + return false; + } else { + std::cout << "glpsol and lp_tst disagree: glpsol status is " << glpk_status; + std::cout << " but lp_tst status is " << lp_tst_status << std::endl; + failures++; + return false; + } + } + return lp_tst_status == "OPTIMAL"; +} + +double get_glpk_cost(std::string file_name) { + std::ifstream f(file_name); + if (!f.is_open()) { + std::cout << "cannot open " << file_name << std::endl; + throw 0; + } + std::string str; + while (getline(f, str)) { + if (str.find("Objective") != std::string::npos) { + vector tokens = split_and_trim(str); + if (tokens.size() != 5) { + std::cout << "unexpected Objective std::string " << str << std::endl; + throw 0; + } + return atof(tokens[3].c_str()); + } + } + std::cout << "cannot find the Objective line in " << file_name << std::endl; + throw 0; +} + +double get_lp_tst_cost(std::string file_name) { + std::ifstream f(file_name); + if (!f.is_open()) { + std::cout << "cannot open " << file_name << std::endl; + throw 0; + } + std::string str; + std::string cost_string; + while (getline(f, str)) { + if (str.find("cost") != std::string::npos) { + cost_string = str; + } + } + if (cost_string.empty()) { + std::cout << "cannot find the cost line in " << file_name << std::endl; + throw 0; + } + + vector tokens = split_and_trim(cost_string); + if (tokens.size() != 3) { + std::cout << "unexpected cost string " << cost_string << std::endl; + throw 0; + } + return atof(tokens[2].c_str()); +} + +bool values_are_one_percent_close(double a, double b) { + double maxval = std::max(fabs(a), fabs(b)); + if (maxval < 0.000001) { + return true; + } + + double one_percent = maxval / 100; + return fabs(a - b) <= one_percent; +} + +// returns true if both are optimal +void compare_costs(std::string glpk_out_file_name, + std::string lp_out_file_name, + unsigned & successes, + unsigned & failures) { + double a = get_glpk_cost(glpk_out_file_name); + double b = get_lp_tst_cost(lp_out_file_name); + + if (values_are_one_percent_close(a, b)) { + successes++; + } else { + failures++; + std::cout << "glpsol cost is " << a << " lp_tst cost is " << b << std::endl; + } +} +void compare_with_glpk(std::string glpk_out_file_name, std::string lp_out_file_name, unsigned & successes, unsigned & failures, std::string /*lp_file_name*/) { +#ifdef CHECK_GLPK_SOLUTION + std::unordered_map * solution_table = get_solution_from_glpsol_output(glpk_out_file_name); + if (solution_is_feasible(lp_file_name, *solution_table)) { + std::cout << "glpk solution is feasible" << std::endl; + } else { + std::cout << "glpk solution is infeasible" << std::endl; + } + delete solution_table; +#endif + if (compare_statuses(glpk_out_file_name, lp_out_file_name, successes, failures)) { + compare_costs(glpk_out_file_name, lp_out_file_name, successes, failures); + } +} +void test_lar_on_file(std::string file_name, argument_parser & args_parser); + +void process_test_file(std::string test_dir, std::string test_file_name, argument_parser & args_parser, std::string out_dir, unsigned max_iters, unsigned time_limit, unsigned & successes, unsigned & failures, unsigned & inconclusives) { + bool use_mpq = args_parser.option_is_used("--mpq"); + bool minimize = args_parser.option_is_used("--min"); + std::string full_lp_tst_out_name = out_dir + "/" + create_output_file_name(minimize, test_file_name, use_mpq); + + std::string input_file_name = test_dir + "/" + test_file_name; + if (input_file_name[input_file_name.size() - 1] == '~') { + // std::cout << "ignoring " << input_file_name << std::endl; + return; + } + std::cout <<"processing " << input_file_name << std::endl; + + std::ofstream out(full_lp_tst_out_name); + if (!out.is_open()) { + std::cout << "cannot open file " << full_lp_tst_out_name << std::endl; + throw 0; + } + std::streambuf *coutbuf = std::cout.rdbuf(); // save old buffer + std::cout.rdbuf(out.rdbuf()); // redirect std::cout to dir_entry->d_name! + bool dual = args_parser.option_is_used("--dual"); + try { + if (args_parser.option_is_used("--lar")) + test_lar_on_file(input_file_name, args_parser); + else + solve_mps(input_file_name, minimize, time_limit, use_mpq, dual, false, args_parser); + } + catch(...) { + std::cout << "catching the failure" << std::endl; + failures++; + std::cout.rdbuf(coutbuf); // reset to standard output again + return; + } + std::cout.rdbuf(coutbuf); // reset to standard output again + + if (args_parser.option_is_used("--compare_with_glpk")) { + std::string glpk_out_file_name = out_dir + "/" + create_output_file_name_for_glpsol(minimize, std::string(test_file_name)); + int glpk_exit_code = run_glpk(input_file_name, glpk_out_file_name, minimize, time_limit); + if (glpk_exit_code != 0) { + std::cout << "glpk failed" << std::endl; + inconclusives++; + } else { + compare_with_glpk(glpk_out_file_name, full_lp_tst_out_name, successes, failures, input_file_name); + } + } +} +/* + 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 + } +*/ +/* + 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; + + 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 { bool operator()(const std::pair &left, const std::pair &right) { return left.second < right.second; @@ -2008,11 +2598,121 @@ 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); + 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::unordered_map get_solution_map(lp_solver * lps, mps_reader & reader) { + std::unordered_map ret; + for (const auto & it : reader.column_names()) { + ret[it] = lps->get_column_value_by_name(it); + } + return ret; +} +void run_lar_solver(argument_parser & args_parser, lar_solver * solver, mps_reader * reader) { + std::string maxng = args_parser.get_option_value("--maxng"); + if (!maxng.empty()) { + solver->settings().max_number_of_iterations_with_no_improvements = atoi(maxng.c_str()); + } + if (args_parser.option_is_used("-pd")){ + solver->settings().presolve_with_double_solver_for_lar = true; + } + + if (args_parser.option_is_used("--compare_with_primal")){ + if (reader == nullptr) { + std::cout << "cannot compare with primal, the reader is null " << std::endl; + return; + } + auto * lps = reader->create_solver(false); + lps->find_maximal_solution(); + std::unordered_map sol = get_solution_map(lps, *reader); + std::cout << "status = " << lp_status_to_string(solver->get_status()) << std::endl; + return; + } + stopwatch sw; + sw.start(); + lp_status status = solver->solve(); + std::cout << "status is " << lp_status_to_string(status) << ", processed for " << sw.get_current_seconds() <<" seconds, and " << solver->get_total_iterations() << " iterations" << std::endl; + if (solver->get_status() == lp_status::INFEASIBLE) { + explanation evidence; + solver->get_infeasibility_explanation(evidence); + } + if (args_parser.option_is_used("--randomize_lar")) { + if (solver->get_status() != lp_status::OPTIMAL) { + std::cout << "cannot check randomize on an infeazible problem" << std::endl; + return; + } + std::cout << "checking randomize" << std::endl; + vector all_vars; + for (unsigned j = 0; j < solver->number_of_vars(); j++) + all_vars.push_back(j); + + unsigned m = all_vars.size(); + if (m > 100) + m = 100; + + var_index *vars = new var_index[m]; + for (unsigned i = 0; i < m; i++) + vars[i]=all_vars[i]; + + solver->random_update(m, vars); + delete []vars; + } +} +lar_solver * create_lar_solver_from_file(std::string file_name, argument_parser & args_parser) { + if (args_parser.option_is_used("--smt")) { + smt_reader reader(file_name); + reader.read(); + if (!reader.is_ok()){ + std::cout << "cannot process " << file_name << std::endl; + return nullptr; + } + return reader.create_lar_solver(); + } + mps_reader reader(file_name); + reader.read(); + if (!reader.is_ok()) { + std::cout << "cannot process " << file_name << std::endl; + return nullptr; + } + return reader.create_lar_solver(); +} +void test_lar_on_file(std::string file_name, argument_parser & args_parser) { + lar_solver * solver = create_lar_solver_from_file(file_name, args_parser); + mps_reader reader(file_name); + mps_reader * mps_reader = nullptr; + reader.read(); + if (reader.is_ok()) { + mps_reader = & reader; + run_lar_solver(args_parser, solver, mps_reader); + } + delete solver; +} vector get_file_names_from_file_list(std::string filelist) { std::ifstream file(filelist); @@ -2033,6 +2733,23 @@ vector get_file_names_from_file_list(std::string filelist) { return ret; } +void test_lar_solver(argument_parser & args_parser) { + + std::string file_name = args_parser.get_option_value("--file"); + if (!file_name.empty()) { + test_lar_on_file(file_name, args_parser); + return; + } + + std::string file_list = args_parser.get_option_value("--filelist"); + if (!file_list.empty()) { + for (const std::string & fn : get_file_names_from_file_list(file_list)) + test_lar_on_file(fn, args_parser); + return; + } + + std::cout << "give option --file or --filelist to test_lar_solver\n"; +} void test_numeric_pair() { numeric_pair a; @@ -3217,6 +3934,22 @@ void test_lp_local(int argn, char**argv) { return finalize(0); } + if (args_parser.option_is_used("--test_mpq")) { + test_rationals(); + return finalize(0); + } + + if (args_parser.option_is_used("--test_mpq_np")) { + test_rationals_no_numeric_pairs(); + return finalize(0); + } + + if (args_parser.option_is_used("--test_mpq_np_plus")) { + test_rationals_no_numeric_pairs_plus(); + return finalize(0); + } + + if (args_parser.option_is_used("--test_int_set")) { test_int_set(); @@ -3234,8 +3967,29 @@ void test_lp_local(int argn, char**argv) { return finalize(0); } - - +#ifdef Z3DEBUG + if (args_parser.option_is_used("--test_swaps")) { + square_sparse_matrix m(10, 0); + fill_matrix(m); + test_swap_rows_with_permutation(m); + test_swap_cols_with_permutation(m); + return finalize(0); + } +#endif + if (args_parser.option_is_used("--test_perm")) { + test_permutations(); + return finalize(0); + } + if (args_parser.option_is_used("--test_file_directory")) { + test_files_from_directory(args_parser.get_option_value("--test_file_directory"), args_parser); + return finalize(0); + } + std::string file_list = args_parser.get_option_value("--filelist"); + if (!file_list.empty()) { + for (const std::string & fn : get_file_names_from_file_list(file_list)) + solve_mps(fn, args_parser); + return finalize(0); + } if (args_parser.option_is_used("-tbq")) { test_binary_priority_queue(); @@ -3243,6 +3997,100 @@ void test_lp_local(int argn, char**argv) { return finalize(ret); } +#ifdef Z3DEBUG + lp_settings settings; + update_settings(args_parser, settings); + if (args_parser.option_is_used("--test_lu")) { + test_lu(settings); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--test_small_lu")) { + test_small_lu(settings); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--lar")){ + std::cout <<"calling test_lar_solver" << std::endl; + test_lar_solver(args_parser); + return finalize(0); + } + + + + if (args_parser.option_is_used("--test_larger_lu")) { + test_larger_lu(settings); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--test_larger_lu_with_holes")) { + test_larger_lu_with_holes(settings); + ret = 0; + return finalize(ret); + } +#endif + if (args_parser.option_is_used("--eti")) { + test_evidence_for_total_inf_simple(args_parser); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--maximize_term")) { + test_maximize_term(); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--test_lp_0")) { + test_lp_0(); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--smap")) { + test_stacked(); + ret = 0; + return finalize(ret); + } + if (args_parser.option_is_used("--term")) { + test_term(); + ret = 0; + return finalize(ret); + } + unsigned time_limit; + get_time_limit_and_max_iters_from_parser(args_parser, time_limit); + bool dual = args_parser.option_is_used("--dual"); + bool solve_for_rational = args_parser.option_is_used("--mpq"); + std::string file_name = args_parser.get_option_value("--file"); + if (!file_name.empty()) { + solve_mps(file_name, args_parser.option_is_used("--min"), time_limit, solve_for_rational, dual, args_parser.option_is_used("--compare_with_primal"), args_parser); + ret = 0; + return finalize(ret); + } + + if (args_parser.option_is_used("--solve_some_mps")) { +#ifndef _WINDOWS + solve_some_mps(args_parser); +#endif + ret = 0; + return finalize(ret); + } + // lp::ccc = 0; + return finalize(0); + test_init_U(); + test_replace_column(); +#ifdef Z3DEBUG + square_sparse_matrix_with_permutations_test(); + test_dense_matrix(); + test_swap_operations(); + test_permutations(); + test_pivot_like_swaps_and_pivot(); +#endif + tst1(); + std::cout << "done with LP tests\n"; return finalize(0); // has_violations() ? 1 : 0); } } diff --git a/src/test/lp/smt_reader.h b/src/test/lp/smt_reader.h index 75edb23b9..2ab0c1ea6 100644 --- a/src/test/lp/smt_reader.h +++ b/src/test/lp/smt_reader.h @@ -20,14 +20,18 @@ Revision History: #pragma once +// reads an MPS file representing a Mixed Integer Program #include #include #include +#include "math/lp/lp_primal_simplex.h" +#include "math/lp/lp_dual_simplex.h" #include "math/lp/lar_solver.h" #include #include #include #include +#include "math/lp/mps_reader.h" #include "math/lp/ul_pair.h" #include "math/lp/lar_constraints.h" #include diff --git a/src/test/lp/test_file_reader.h b/src/test/lp/test_file_reader.h index 36b273740..8f461ea1c 100644 --- a/src/test/lp/test_file_reader.h +++ b/src/test/lp/test_file_reader.h @@ -27,6 +27,7 @@ Revision History: #include #include #include "math/lp/lp_utils.h" +#include "math/lp/lp_solver.h" namespace lp {