mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-10 19:27:06 +00:00
Code Activity

before rm lu

Browse source 
Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2023-03-04 11:30:59 -08:00
parent f7c9c9ef72
commit ea16f6608c
25 changed files with 70 additions and 3878 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -20,11 +20,8 @@ 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

11
src/math/lp/indexed_value.h
View file

@ -43,15 +43,4 @@ public:
m_value = val;
}
};
#ifdef Z3DEBUG
template <typename X>
bool check_vector_for_small_values(indexed_vector<X> & 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
}

30
src/math/lp/lar_core_solver.h
View file

@ -651,35 +651,7 @@ public:
}
}
void scale_problem_for_doubles(
static_matrix<double, double>& A,
vector<double> & lower_bounds,
vector<double> & upper_bounds) {
vector<double> column_scale_vector;
vector<double> right_side_vector(A.column_count());
settings().reps_in_scaler = 5;
scaler<double, double > 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<unsigned> 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) {

24
src/math/lp/lp_dual_simplex.cpp
View file

@ -1,24 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "math/lp/lp_dual_simplex_def.h"
template lp::mpq lp::lp_dual_simplex<lp::mpq, lp::mpq>::get_current_cost() const;
template void lp::lp_dual_simplex<lp::mpq, lp::mpq>::find_maximal_solution();
template double lp::lp_dual_simplex<double, double>::get_current_cost() const;
template void lp::lp_dual_simplex<double, double>::find_maximal_solution();

93
src/math/lp/lp_dual_simplex.h
View file

@ -1,93 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<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 <typename T, typename X>
class lp_dual_simplex: public lp_solver<T, X> {
lp_dual_core_solver<T, X> * m_core_solver;
vector<T> m_b_copy;
vector<T> 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<column_type> m_column_types_of_core_solver;
vector<column_type> m_column_types_of_logicals;
vector<bool> 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;
};
}

376
src/math/lp/lp_dual_simplex_def.h
View file

@ -1,376 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "math/lp/lp_dual_simplex.h"
namespace lp{
template <typename T, typename X> void lp_dual_simplex<T, X>::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 <typename T, typename X> void lp_dual_simplex<T, X>::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<T>::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<T>::zero();
m_column_types_of_core_solver[j] = column_type::fixed;
m_can_enter_basis[j] = false;
break;
default:
lp_unreachable();
}
}
template <typename T, typename X> void lp_dual_simplex<T, X>::fix_structural_for_stage2(unsigned j) {
column_info<T> * 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<T>::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<T>::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 <typename T, typename X> void lp_dual_simplex<T, X>::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 <typename T, typename X> void lp_dual_simplex<T, X>::restore_right_sides() {
unsigned i = this->m_A->row_count();
while (i--) {
this->m_b[i] = m_b_copy[i];
}
}
template <typename T, typename X> void lp_dual_simplex<T, X>::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 <typename T, typename X> void lp_dual_simplex<T, X>::fill_x_with_zeros() {
unsigned j = this->m_A->column_count();
while (j--) {
this->m_x[j] = numeric_traits<T>::zero();
}
}
template <typename T, typename X> void lp_dual_simplex<T, X>::stage1() {
lp_assert(m_core_solver == nullptr);
this->m_x.resize(this->m_A->column_count(), numeric_traits<T>::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<T, X>(
*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 <typename T, typename X> void lp_dual_simplex<T, X>::stage2() {
unmark_boxed_and_fixed_columns_and_fix_structural_costs();
restore_right_sides();
solve_for_stage2();
}
template <typename T, typename X> void lp_dual_simplex<T, X>::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 <typename T, typename X> column_type lp_dual_simplex<T, X>::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 <typename T, typename X> void lp_dual_simplex<T, X>::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<T> * 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<T>::zero();
this->m_upper_bounds[j] = numeric_traits<T>::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<T>::zero();
this->m_upper_bounds[j] = this->m_lower_bounds[j] = numeric_traits<T>::zero(); // is it needed?
break;
default:
lp_unreachable();
}
m_column_types_of_core_solver[j] = column_type::boxed;
}
template <typename T, typename X> void lp_dual_simplex<T, X>::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<T>::zero();
this->m_upper_bounds[j] = numeric_traits<T>::one();
} else {
m_column_types_of_core_solver[j] = column_type::fixed;
this->m_lower_bounds[j] = numeric_traits<T>::zero();
this->m_upper_bounds[j] = numeric_traits<T>::zero();
}
}
template <typename T, typename X> void lp_dual_simplex<T, X>::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 <typename T, typename X> void lp_dual_simplex<T, X>::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<T>::zero();
(*this->m_A)(row, artificial) = numeric_traits<T>::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<T>::one();
if (rs > 0) {
// adding one artificial
set_type_for_logical(artificial, column_type::fixed);
(*this->m_A)(row, artificial) = numeric_traits<T>::one();
this->m_basis[row] = artificial;
this->m_costs[artificial] = numeric_traits<T>::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<T>::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<T>::one();
if (rs < 0) {
// adding one artificial
set_type_for_logical(artificial, column_type::fixed);
(*this->m_A)(row, artificial) = - numeric_traits<T>::one();
this->m_basis[row] = artificial;
this->m_costs[artificial] = numeric_traits<T>::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<T>::zero();
}
slack_var++;
break;
}
}
template <typename T, typename X> void lp_dual_simplex<T, X>::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 <typename T, typename X> void lp_dual_simplex<T, X>::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<T>::zero(); // preparing for the first stage
}
}
template <typename T, typename X> void lp_dual_simplex<T, X>::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 <typename T, typename X> T lp_dual_simplex<T, X>::get_current_cost() const {
T ret = numeric_traits<T>::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
}
}

1
src/math/lp/lp_primal_core_solver.h
View file

@ -30,7 +30,6 @@ Revision History:
#include <cstdlib>
#include <algorithm>
#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"

35
src/math/lp/lp_primal_simplex.cpp
View file

@ -1,35 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <utility>
#include <memory>
#include <string>
#include "util/vector.h"
#include <functional>
#include "math/lp/lp_primal_simplex_def.h"
template bool lp::lp_primal_simplex<double, double>::bounds_hold(std::unordered_map<std::string, double, std::hash<std::string>, std::equal_to<std::string>, std::allocator<std::pair<std::string const, double> > > const&);
template bool lp::lp_primal_simplex<double, double>::row_constraints_hold(std::unordered_map<std::string, double, std::hash<std::string>, std::equal_to<std::string>, std::allocator<std::pair<std::string const, double> > > const&);
template double lp::lp_primal_simplex<double, double>::get_current_cost() const;
template double lp::lp_primal_simplex<double, double>::get_column_value(unsigned int) const;
template lp::lp_primal_simplex<double, double>::~lp_primal_simplex();
template lp::lp_primal_simplex<lp::mpq, lp::mpq>::~lp_primal_simplex();
template lp::mpq lp::lp_primal_simplex<lp::mpq, lp::mpq>::get_current_cost() const;
template lp::mpq lp::lp_primal_simplex<lp::mpq, lp::mpq>::get_column_value(unsigned int) const;
template void lp::lp_primal_simplex<double, double>::find_maximal_solution();
template void lp::lp_primal_simplex<lp::mpq, lp::mpq>::find_maximal_solution();

106
src/math/lp/lp_primal_simplex.h
View file

@ -1,106 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
#include <unordered_map>
#include <string>
#include <algorithm>
#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 <typename T, typename X>
class lp_primal_simplex: public lp_solver<T, X> {
lp_primal_core_solver<T, X> * m_core_solver;
vector<X> 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<T> & r);
void refactor();
void set_scaled_costs();
public:
lp_primal_simplex(): m_core_solver(nullptr) {}
column_info<T> * 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<std::string, T> const & solution);
T get_row_value(unsigned i, std::unordered_map<std::string, T> const & solution, std::ostream * out);
bool row_constraint_holds(unsigned i, std::unordered_map<std::string, T> const & solution, std::ostream * out);
bool row_constraints_hold(std::unordered_map<std::string, T> const & solution);
T * get_array_from_map(std::unordered_map<std::string, T> const & solution);
bool solution_is_feasible(std::unordered_map<std::string, T> 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;
};
}

367
src/math/lp/lp_primal_simplex_def.h
View file

@ -1,367 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <string>
#include "util/vector.h"
#include "math/lp/lp_primal_simplex.h"
namespace lp {
template <typename T, typename X> void lp_primal_simplex<T, X>::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 <typename T, typename X> void lp_primal_simplex<T, X>::init_buffer(unsigned k, vector<T> & 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 <typename T, typename X> void lp_primal_simplex<T, X>::refactor() {
m_core_solver->init_lu();
if (m_core_solver->factorization()->get_status() != LU_status::OK) {
throw_exception("cannot refactor");
}
}
template <typename T, typename X> void lp_primal_simplex<T, X>::set_scaled_costs() {
unsigned j = this->number_of_core_structurals();
while (j-- > 0) {
this->set_scaled_cost(j);
}
}
template <typename T, typename X> column_info<T> * lp_primal_simplex<T, X>::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<T>) : it->second;
}
template <typename T, typename X> void lp_primal_simplex<T, X>::fill_acceptable_values_for_x() {
for (auto t : this->m_core_solver_columns_to_external_columns) {
this->m_x[t.first] = numeric_traits<T>::zero();
}
}
template <typename T, typename X> void lp_primal_simplex<T, X>::set_zero_bound(bool * bound_is_set, T * bounds, unsigned i) {
bound_is_set[i] = true;
bounds[i] = numeric_traits<T>::zero();
}
template <typename T, typename X> void lp_primal_simplex<T, X>::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<T>::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<T>::one();
this->m_x[artificial] = rs;
} else {
(*this->m_A)(row, artificial) = - numeric_traits<T>::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<T>::one();
if (rs > 0) {
lp_assert(numeric_traits<T>::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<T>::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 <typename T, typename X> void lp_primal_simplex<T, X>::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<T>::one();
if (rs < 0) {
// adding one artificial
lp_assert(numeric_traits<T>::is_zero(this->m_x[slack_var]));
this->m_column_types[artificial] = column_type::lower_bound;
(*this->m_A)(row, artificial) = - numeric_traits<T>::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 <typename T, typename X> void lp_primal_simplex<T, X>::adding an artificial variable
this->m_basis[row] = slack_var;
this->m_x[slack_var] = rs;
}
slack_var++;
break;
}
}
template <typename T, typename X> void lp_primal_simplex<T, X>::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<T> * 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<T>::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 <typename T, typename X> void lp_primal_simplex<T, X>::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 <typename T, typename X> void lp_primal_simplex<T, X>::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 <typename T, typename X> void lp_primal_simplex<T, X>::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<T>::one();
this->m_column_types[j] = column_type::fixed;
this->m_upper_bounds[j] = m_lower_bounds[j] = zero_of_type<X>();
break;
case Greater_or_equal:
this->m_x[j] = - this->m_b[row];
(*this->m_A)(row, j) = - numeric_traits<T>::one();
this->m_column_types[j] = column_type::lower_bound;
this->m_upper_bounds[j] = zero_of_type<X>();
break;
case Less_or_equal:
this->m_x[j] = this->m_b[row];
(*this->m_A)(row, j) = numeric_traits<T>::one();
this->m_column_types[j] = column_type::lower_bound;
this->m_upper_bounds[j] = m_lower_bounds[j] = zero_of_type<X>();
break;
default:
lp_unreachable();
}
}
template <typename T, typename X> void lp_primal_simplex<T, X>::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<X>()); // low bounds are shifted ot zero
this->m_x.resize(total_vars, numeric_traits<T>::zero());
this->m_basis.resize(this->row_count());
this->m_costs.clear();
this->m_costs.resize(total_vars, zero_of_type<T>());
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<T, X>(*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 <typename T, typename X> lp_primal_simplex<T, X>::~lp_primal_simplex() {
delete m_core_solver;
}
template <typename T, typename X> bool lp_primal_simplex<T, X>::bounds_hold(std::unordered_map<std::string, T> 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 <typename T, typename X> T lp_primal_simplex<T, X>::get_row_value(unsigned i, std::unordered_map<std::string, T> 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<T>::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<T> * 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 <typename T, typename X> bool lp_primal_simplex<T, X>::row_constraint_holds(unsigned i, std::unordered_map<std::string, T> 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<T>::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<T>::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<T>::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 <typename T, typename X> bool lp_primal_simplex<T, X>::row_constraints_hold(std::unordered_map<std::string, T> 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 <typename T, typename X> T lp_primal_simplex<T, X>::get_current_cost() const {
T ret = numeric_traits<T>::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;
}
}

55
src/math/lp/lp_solver.cpp
View file

@ -1,55 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <string>
#include "math/lp/lp_solver_def.h"
template void lp::lp_solver<double, double>::add_constraint(lp::lp_relation, double, unsigned int);
template void lp::lp_solver<double, double>::cleanup();
template void lp::lp_solver<double, double>::count_slacks_and_artificials();
template void lp::lp_solver<double, double>::fill_m_b();
template void lp::lp_solver<double, double>::fill_matrix_A_and_init_right_side();
template void lp::lp_solver<double, double>::flip_costs();
template double lp::lp_solver<double, double>::get_column_cost_value(unsigned int, lp::column_info<double>*) const;
template int lp::lp_solver<double, double>::get_column_index_by_name(std::string) const;
template double lp::lp_solver<double, double>::get_column_value_with_core_solver(unsigned int, lp::lp_core_solver_base<double, double>*) const;
template lp::column_info<double>* lp::lp_solver<double, double>::get_or_create_column_info(unsigned int);
template void lp::lp_solver<double, double>::give_symbolic_name_to_column(std::string, unsigned int);
template void lp::lp_solver<double, double>::print_statistics_on_A(std::ostream & out);
template bool lp::lp_solver<double, double>::problem_is_empty();
template void lp::lp_solver<double, double>::scale();
template void lp::lp_solver<double, double>::set_scaled_cost(unsigned int);
template lp::lp_solver<double, double>::~lp_solver();
template void lp::lp_solver<lp::mpq, lp::mpq>::add_constraint(lp::lp_relation, lp::mpq, unsigned int);
template void lp::lp_solver<lp::mpq, lp::mpq>::cleanup();
template void lp::lp_solver<lp::mpq, lp::mpq>::count_slacks_and_artificials();
template void lp::lp_solver<lp::mpq, lp::mpq>::fill_m_b();
template void lp::lp_solver<lp::mpq, lp::mpq>::fill_matrix_A_and_init_right_side();
template void lp::lp_solver<lp::mpq, lp::mpq>::flip_costs();
template lp::mpq lp::lp_solver<lp::mpq, lp::mpq>::get_column_cost_value(unsigned int, lp::column_info<lp::mpq>*) const;
template int lp::lp_solver<lp::mpq, lp::mpq>::get_column_index_by_name(std::string) const;
template lp::mpq lp::lp_solver<lp::mpq, lp::mpq>::get_column_value_by_name(std::string) const;
template lp::mpq lp::lp_solver<lp::mpq, lp::mpq>::get_column_value_with_core_solver(unsigned int, lp::lp_core_solver_base<lp::mpq, lp::mpq>*) const;
template lp::column_info<lp::mpq>* lp::lp_solver<lp::mpq, lp::mpq>::get_or_create_column_info(unsigned int);
template void lp::lp_solver<lp::mpq, lp::mpq>::give_symbolic_name_to_column(std::string, unsigned int);
template void lp::lp_solver<lp::mpq, lp::mpq>::print_statistics_on_A(std::ostream & out);
template bool lp::lp_solver<lp::mpq, lp::mpq>::problem_is_empty();
template void lp::lp_solver<lp::mpq, lp::mpq>::scale();
template void lp::lp_solver<lp::mpq, lp::mpq>::set_scaled_cost(unsigned int);
template lp::lp_solver<lp::mpq, lp::mpq>::~lp_solver();
template double lp::lp_solver<double, double>::get_column_value_by_name(std::string) const;

260
src/math/lp/lp_solver.h
View file

@ -1,260 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <string>
#include <unordered_map>
#include <algorithm>
#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 <typename T, typename X>
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 <typename T, typename X>
class lp_solver : public column_namer {
column_info<T> * get_or_create_column_info(unsigned column);
protected:
T get_column_cost_value(unsigned j, column_info<T> * ci) const;
public:
unsigned m_total_iterations;
static_matrix<T, X>* m_A; // this is the matrix of constraints
vector<T> m_b; // the right side vector
unsigned m_first_stage_iterations;
unsigned m_second_stage_iterations;
std::unordered_map<unsigned, lp_constraint<T, X>> m_constraints;
std::unordered_map<var_index, column_info<T>*> m_map_from_var_index_to_column_info;
std::unordered_map<unsigned, std::unordered_map<unsigned, T> > m_A_values;
std::unordered_map<std::string, unsigned> m_names_to_columns; // don't have to use it
std::unordered_map<unsigned, unsigned> m_external_rows_to_core_solver_rows;
std::unordered_map<unsigned, unsigned> m_core_solver_rows_to_external_rows;
std::unordered_map<unsigned, unsigned> m_core_solver_columns_to_external_columns;
vector<T> m_column_scale;
std::unordered_map<unsigned, std::string> m_name_map;
unsigned m_artificials;
unsigned m_slacks;
vector<column_type> m_column_types;
vector<T> m_costs;
vector<T> m_x;
vector<T> m_upper_bounds;
vector<unsigned> m_basis;
vector<unsigned> m_nbasis;
vector<int> 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<T> *ci = get_or_create_column_info(i);
ci->set_lower_bound(bound);
}
void set_upper_bound(unsigned i, T bound) {
column_info<T> *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<T> *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<unsigned, T> & row_map);
void pin_vars_down_on_row(std::unordered_map<unsigned, T> & row) {
pin_vars_on_row_with_sign(row, - numeric_traits<T>::one());
}
void pin_vars_up_on_row(std::unordered_map<unsigned, T> & row) {
pin_vars_on_row_with_sign(row, numeric_traits<T>::one());
}
void pin_vars_on_row_with_sign(std::unordered_map<unsigned, T> & row, T sign );
bool get_minimal_row_value(std::unordered_map<unsigned, T> & row, T & lower_bound);
bool get_maximal_row_value(std::unordered_map<unsigned, T> & row, T & lower_bound);
bool row_is_zero(std::unordered_map<unsigned, T> & row);
bool row_e_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index);
bool row_ge_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index);
bool row_le_is_obsolete(std::unordered_map<unsigned, T> & 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<unsigned, T> & row, unsigned row_index );
void remove_fixed_or_zero_columns();
void remove_fixed_or_zero_columns_from_row(unsigned i, std::unordered_map<unsigned, T> & 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<unsigned>(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<T>::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<T, X> * 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;
}
}
};
}

571
src/math/lp/lp_solver_def.h
View file

@ -1,571 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <string>
#include <algorithm>
#include "util/vector.h"
#include "math/lp/lp_solver.h"
namespace lp {
template <typename T, typename X> column_info<T> * lp_solver<T, X>::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<T>()) : it->second;
}
template <typename T, typename X>
std::string lp_solver<T, X>::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 <typename T, typename X> T lp_solver<T, X>::get_column_cost_value(unsigned j, column_info<T> * ci) const {
if (ci->is_fixed()) {
return ci->get_cost() * ci->get_fixed_value();
}
return ci->get_cost() * get_column_value(j);
}
template <typename T, typename X> void lp_solver<T, X>::add_constraint(lp_relation relation, T right_side, unsigned row_index) {
lp_assert(m_constraints.find(row_index) == m_constraints.end());
lp_constraint<T, X> cs(right_side, relation);
m_constraints[row_index] = cs;
}
template <typename T, typename X> void lp_solver<T, X>::give_symbolic_name_to_column(std::string name, unsigned column) {
auto it = m_map_from_var_index_to_column_info.find(column);
column_info<T> *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<T>;
} else {
ci = it->second;
}
ci->set_name(name);
m_names_to_columns[name] = column;
}
template <typename T, typename X> T lp_solver<T, X>::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 <typename T, typename X> int lp_solver<T, X>::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 <typename T, typename X> lp_solver<T, X>::~lp_solver(){
delete m_A;
for (auto t : m_map_from_var_index_to_column_info) {
delete t.second;
}
}
template <typename T, typename X> void lp_solver<T, X>::flip_costs() {
for (auto t : m_map_from_var_index_to_column_info) {
column_info<T> *ci = t.second;
ci->set_cost(-ci->get_cost());
}
}
template <typename T, typename X> bool lp_solver<T, X>::problem_is_empty() {
for (auto & c : m_A_values)
if (!c.second.empty())
return false;
return true;
}
template <typename T, typename X> void lp_solver<T, X>::scale() {
if (numeric_traits<T>::precise() || m_settings.use_scaling == false) {
m_column_scale.clear();
m_column_scale.resize(m_A->column_count(), one_of_type<T>());
return;
}
T smin = T(m_settings.scaling_minimum);
T smax = T(m_settings.scaling_maximum);
scaler<T, X> scaler(m_b, *m_A, smin, smax, m_column_scale, this->m_settings);
if (!scaler.scale()) {
unscale();
}
}
template <typename T, typename X> void lp_solver<T, X>::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 <typename T, typename X> void lp_solver<T, X>::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 <typename T, typename X> void lp_solver<T, X>::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 <typename T, typename X> void lp_solver<T, X>::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 <typename T, typename X> void lp_solver<T, X>::print_scale_stats(std::ostream & out) {
print_rows_scale_stats(out);
print_columns_scale_stats(out);
}
template <typename T, typename X> void lp_solver<T, X>::get_max_abs_in_row(std::unordered_map<unsigned, T> & row_map) {
T ret = numeric_traits<T>::zero();
for (auto jp : row_map) {
T ac = numeric_traits<T>::abs(jp->second);
if (ac > ret) {
ret = ac;
}
}
return ret;
}
template <typename T, typename X> void lp_solver<T, X>::pin_vars_on_row_with_sign(std::unordered_map<unsigned, T> & row, T sign ) {
for (auto t : row) {
unsigned j = t.first;
column_info<T> * ci = m_map_from_var_index_to_column_info[j];
T a = t.second;
if (a * sign > numeric_traits<T>::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 <typename T, typename X> bool lp_solver<T, X>::get_minimal_row_value(std::unordered_map<unsigned, T> & row, T & lower_bound) {
lower_bound = numeric_traits<T>::zero();
for (auto & t : row) {
T a = t.second;
column_info<T> * ci = m_map_from_var_index_to_column_info[t.first];
if (a > numeric_traits<T>::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 <typename T, typename X> bool lp_solver<T, X>::get_maximal_row_value(std::unordered_map<unsigned, T> & row, T & lower_bound) {
lower_bound = numeric_traits<T>::zero();
for (auto & t : row) {
T a = t.second;
column_info<T> * ci = m_map_from_var_index_to_column_info[t.first];
if (a < numeric_traits<T>::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 <typename T, typename X> bool lp_solver<T, X>::row_is_zero(std::unordered_map<unsigned, T> & row) {
for (auto & t : row) {
if (!is_zero(t.second))
return false;
}
return true;
}
template <typename T, typename X> bool lp_solver<T, X>::row_e_is_obsolete(std::unordered_map<unsigned, T> & 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 <typename T, typename X> bool lp_solver<T, X>::row_ge_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index) {
T rs = m_constraints[row_index].m_rs;
if (row_is_zero(row)) {
if (rs > zero_of_type<X>())
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 <typename T, typename X> bool lp_solver<T, X>::row_le_is_obsolete(std::unordered_map<unsigned, T> & 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<X>())
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 <typename T, typename X> bool lp_solver<T, X>::row_is_obsolete(std::unordered_map<unsigned, T> & 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 <typename T, typename X> void lp_solver<T, X>::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 <typename T, typename X> void lp_solver<T, X>::remove_fixed_or_zero_columns_from_row(unsigned i, std::unordered_map<unsigned, T> & row) {
auto & constraint = m_constraints[i];
vector<unsigned> 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<T> * 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<T>::is_zero(col.second)){
removed.push_back(j);
}
}
}
for (auto j : removed) {
row.erase(j);
}
}
template <typename T, typename X> unsigned lp_solver<T, X>::try_to_remove_some_rows() {
vector<unsigned> 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<unsigned>(rows_to_delete.size());
}
template <typename T, typename X> void lp_solver<T, X>::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 <typename T, typename X> void lp_solver<T, X>::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 <typename T, typename X> void lp_solver<T, X>::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<T>::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 <typename T, typename X> void lp_solver<T, X>::unscale() {
delete m_A;
m_A = nullptr;
fill_A_from_A_values();
restore_column_scales_to_one();
fill_m_b();
}
template <typename T, typename X> void lp_solver<T, X>::fill_A_from_A_values() {
m_A = new static_matrix<T, X>(static_cast<unsigned>(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<T> *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 <typename T, typename X> void lp_solver<T, X>::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 <typename T, typename X> void lp_solver<T, X>::count_slacks_and_artificials() {
for (int i = row_count() - 1; i >= 0; i--) {
count_slacks_and_artificials_for_row(i);
}
}
template <typename T, typename X> void lp_solver<T, X>::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 <typename T, typename X> T lp_solver<T, X>::lower_bound_shift_for_row(unsigned i) {
T ret = numeric_traits<T>::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 <typename T, typename X> void lp_solver<T, X>::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 <typename T, typename X> T lp_solver<T, X>::get_column_value_with_core_solver(unsigned column, lp_core_solver_base<T, X> * 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<T>::zero();
}
column_info<T> * ci = cit->second;
if (ci->is_fixed()) {
return ci->get_fixed_value();
}
unsigned cj = ci->get_column_index();
if (cj != static_cast<unsigned>(-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<T>::zero(); // returns zero for out of boundary columns
}
template <typename T, typename X> void lp_solver<T, X>::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<T> * 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];
}
}

891
src/math/lp/mps_reader.h
View file

@ -1,891 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
// reads an MPS file representing a Mixed Integer Program
#include <functional>
#include <algorithm>
#include <string>
#include "util/vector.h"
#include <unordered_map>
#include <ostream>
#include <fstream>
#include <locale>
#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<int>(s.size()) - 1;i >= 0; i--)
if (my_white_space(s[i])) ret++;else break;
return ret;
}
// trim from start
inline std::string &ltrim(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<int, int>(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<std::string> string_split(const std::string &source, const char *delimiter, bool keep_empty) {
vector<std::string> 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<std::string> split_and_trim(const std::string &line) {
auto split = string_split(line, " \t", false);
vector<std::string> ret;
for (auto s : split) {
ret.push_back(trim(s));
}
return ret;
}
template <typename T, typename X>
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<T>::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<std::string, T> 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<T>()),
m_range(zero_of_type<T>())
{
}
};
bool m_is_OK;
std::string m_file_name;
std::unordered_map<std::string, row *> m_rows;
std::unordered_map<std::string, column *> m_columns;
std::unordered_map<std::string, unsigned> 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<T, X> * solver){
if (b == nullptr) {
solver->set_lower_bound(col, numeric_traits<T>::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<int>(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<T>::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<T>::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<unsigned>(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<T>::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<std::string> 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<T>::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<std::string> 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<std::string> 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<T>::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<T>::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<T>::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<std::string> 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<T>::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<std::string> & 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<T>::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<unsigned>(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<unsigned>(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<T, X> *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<T>::zero() ? -t: t; }
void fill_solver_on_rows(lp_solver<T, X> * solver) {
for (auto row_it : m_rows) {
fill_solver_on_row(row_it.second, solver);
}
}
void fill_solver_on_columns(lp_solver<T, X> * 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<T, X> *solver) {
fill_solver_on_rows(solver);
fill_solver_on_columns(solver);
}
void set_solver_cost(row * row, lp_solver<T, X> *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<std::string> column_names() {
vector<std::string> 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<T, X> * create_solver(bool dual) {
lp_solver<T, X> * solver = dual? (lp_solver<T, X>*)new lp_dual_simplex<T, X>() : new lp_primal_simplex<T, X>();
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<unsigned>(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<std::pair<mpq, var_index>> 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;
}
};
}

1
src/math/lp/static_matrix.cpp
View file

@ -24,7 +24,6 @@ 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"

3
src/sat/smt/arith_sls.h
View file

@ -19,9 +19,6 @@ 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"

4
src/sat/smt/arith_solver.h
View file

@ -19,9 +19,7 @@ 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"

1
src/shell/CMakeLists.txt
View file

@ -28,7 +28,6 @@ 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).

109
src/shell/lp_frontend.cpp
View file

@ -1,109 +0,0 @@
/*++
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 <iostream>
#include <signal.h>
#include "smt/params/smt_params_helper.hpp"
namespace {
static mutex *display_stats_mux = new mutex;
static lp::lp_solver<double, double>* 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<reslimit> eh(rlim);
scoped_timer timer(timeout, &eh);
std::string fn(mps_file_name);
lp::mps_reader<double, double> 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<double, double> * 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;
}

7
src/shell/lp_frontend.h
View file

@ -1,7 +0,0 @@
/*
Copyright (c) 2013 Microsoft Corporation. All rights reserved.
Author: Lev Nachmanson
*/
#pragma once
unsigned read_mps_file(char const * mps_file_name);

10
src/shell/main.cpp
View file

@ -37,14 +37,13 @@ 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 <crtdbg.h>
#endif
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;
typedef enum { IN_UNSPECIFIED, IN_SMTLIB_2, IN_DATALOG, IN_DIMACS, IN_WCNF, IN_OPB, IN_LP, IN_Z3_LOG, IN_DRAT } input_kind;
static char const * g_input_file = nullptr;
static char const * g_drat_input_file = nullptr;
@ -377,10 +376,6 @@ 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) {
@ -406,9 +401,6 @@ 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;

3
src/smt/theory_lra.cpp
View file

@ -19,9 +19,6 @@
--*/
#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"

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

File diff suppressed because it is too large Load diff

4
src/test/lp/smt_reader.h
View file

@ -20,18 +20,14 @@ Revision History:
#pragma once
// reads an MPS file representing a Mixed Integer Program
#include <string>
#include <vector>
#include <unordered_map>
#include "math/lp/lp_primal_simplex.h"
#include "math/lp/lp_dual_simplex.h"
#include "math/lp/lar_solver.h"
#include <iostream>
#include <fstream>
#include <functional>
#include <algorithm>
#include "math/lp/mps_reader.h"
#include "math/lp/ul_pair.h"
#include "math/lp/lar_constraints.h"
#include <sstream>

1
src/test/lp/test_file_reader.h
View file

@ -27,7 +27,6 @@ Revision History:
#include <iostream>
#include <fstream>
#include "math/lp/lp_utils.h"
#include "math/lp/lp_solver.h"
namespace lp {