From 8989e10e7141d52d9696b5c1a39100e16fb0e6e1 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Fri, 3 Mar 2023 15:41:30 -0800 Subject: [PATCH] rm lp_dual_simplex --- 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/sat/smt/arith_solver.h | 1 - 4 files changed, 494 deletions(-) delete mode 100644 src/math/lp/lp_dual_simplex.cpp delete mode 100644 src/math/lp/lp_dual_simplex.h delete mode 100644 src/math/lp/lp_dual_simplex_def.h diff --git a/src/math/lp/lp_dual_simplex.cpp b/src/math/lp/lp_dual_simplex.cpp deleted file mode 100644 index aaf612f56..000000000 --- a/src/math/lp/lp_dual_simplex.cpp +++ /dev/null @@ -1,24 +0,0 @@ -/*++ -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 deleted file mode 100644 index 75ef87492..000000000 --- a/src/math/lp/lp_dual_simplex.h +++ /dev/null @@ -1,93 +0,0 @@ -/*++ -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 deleted file mode 100644 index 8af9d87c1..000000000 --- a/src/math/lp/lp_dual_simplex_def.h +++ /dev/null @@ -1,376 +0,0 @@ -/*++ -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/sat/smt/arith_solver.h b/src/sat/smt/arith_solver.h index 9ae671c16..eac73c719 100644 --- a/src/sat/smt/arith_solver.h +++ b/src/sat/smt/arith_solver.h @@ -21,7 +21,6 @@ Author: #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"