michael/lean2
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

dev(lp): simplify the design of lar_solver

Browse source 
Signed-off-by: Lev Nachmanson <levnach@microsoft.com>
This commit is contained in:
Lev Nachmanson 2016-05-19 14:59:32 -07:00 committed by Leonardo de Moura
parent 1d3c46e712
commit 0bb21914e6
8 changed files with 229 additions and 178 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

13
src/util/lp/canonic_left_side.h
View file

@ -12,8 +12,7 @@
#include <utility> #include <utility>
#include "util/lp/column_info.h" #include "util/lp/column_info.h"
namespace lean { namespace lean {
typedef unsigned var_index;
typedef unsigned constraint_index;
enum lconstraint_kind { enum lconstraint_kind {
LE = -2, LT = -1 , GE = 2, GT = 1, EQ = 0 LE = -2, LT = -1 , GE = 2, GT = 1, EQ = 0
}; };
@ -25,10 +24,10 @@ inline bool compare(const std::pair<mpq, var_index> & a, const std::pair<mpq,
class canonic_left_side { class canonic_left_side {
public: public:
int m_row_index = -1; // this index is exposed to the user, this is not the index that points to the row
int m_column_index = -1; // this is the column of the left side variable in the matrix unsigned m_row_index = static_cast<unsigned>(-1);
var_index m_additional_var_index = static_cast<var_index>(-1); // this is the index of the additional variable created for this constraint
std::vector<std::pair<mpq, var_index>> m_coeffs; std::vector<std::pair<mpq, var_index>> m_coeffs;
column_info<mpq> m_column_info;
lar_normalized_constraint * m_low_bound_witness = nullptr; lar_normalized_constraint * m_low_bound_witness = nullptr;
lar_normalized_constraint * m_upper_bound_witness = nullptr; lar_normalized_constraint * m_upper_bound_witness = nullptr;
@ -42,10 +41,6 @@ public:
normalize(); normalize();
} }
void set_name(std::string name) {
m_column_info.set_name(name);
}
unsigned size() const { return static_cast<unsigned>(m_coeffs.size()); } unsigned size() const { return static_cast<unsigned>(m_coeffs.size()); }
void normalize() { void normalize() {

35
src/util/lp/column_info.h
View file

@ -12,6 +12,10 @@
#include <algorithm> #include <algorithm>
#include "util/lp/lp_settings.h" #include "util/lp/lp_settings.h"
namespace lean { namespace lean {
inline bool is_valid(unsigned j) { return static_cast<int>(j) >= 0;}
typedef unsigned var_index;
typedef unsigned constraint_index;
template <typename T> template <typename T>
class column_info { class column_info {
std::string m_name; std::string m_name;
@ -24,13 +28,40 @@ class column_info {
T m_cost = numeric_traits<T>::zero(); T m_cost = numeric_traits<T>::zero();
T m_fixed_value; T m_fixed_value;
bool m_is_fixed = false; bool m_is_fixed = false;
unsigned m_column_index;
public: public:
void set_column_index(unsigned j) {
m_column_index = j;
}
// the default constructor
column_info() {}
column_info(unsigned column_index) : m_column_index(column_index) {
}
column_info(const column_info & ci) {
m_name = ci.m_name;
m_low_bound_is_set = ci.m_low_bound_is_set;
m_low_bound_is_strict = ci.m_low_bound_is_strict;
m_upper_bound_is_set = ci.m_upper_bound_is_set;
m_upper_bound_is_strict = ci.m_upper_bound_is_strict;
m_low_bound = ci.m_low_bound;
m_upper_bound = ci.m_upper_bound;
m_cost = ci.m_cost;
m_fixed_value = ci.m_fixed_value;
m_is_fixed = ci.m_is_fixed;
m_column_index = ci.m_column_index;
}
unsigned get_column_index() const {
return m_column_index;
}
column_type get_column_type() const { column_type get_column_type() const {
return m_is_fixed? fixed : (m_low_bound_is_set? (m_upper_bound_is_set? boxed : low_bound) : (m_upper_bound_is_set? upper_bound: free_column)); return m_is_fixed? fixed : (m_low_bound_is_set? (m_upper_bound_is_set? boxed : low_bound) : (m_upper_bound_is_set? upper_bound: free_column));
} }
column_type get_column_type_no_flipping() { column_type get_column_type_no_flipping() {
if (m_is_fixed) { if (m_is_fixed) {
return column_type::fixed; return column_type::fixed;
} }

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

@ -33,9 +33,15 @@ double conversion_helper <double>::get_upper_bound(const column_info<mpq> & ci)
canonic_left_side * lar_solver::create_or_fetch_existing_left_side(const buffer<std::pair<mpq, var_index>>& left_side_par) { canonic_left_side * lar_solver::create_or_fetch_existing_left_side(const buffer<std::pair<mpq, var_index>>& left_side_par) {
auto left_side = new canonic_left_side(left_side_par); auto left_side = new canonic_left_side(left_side_par);
lean_assert(left_side->size() > 0); lean_assert(left_side->size() > 0);
auto it = m_canonic_left_sides.find(left_side); auto it = m_set_of_canonic_left_sides.find(left_side);
if (it == m_canonic_left_sides.end()) { if (it == m_set_of_canonic_left_sides.end()) {
m_canonic_left_sides.insert(left_side); m_set_of_canonic_left_sides.insert(left_side);
lean_assert(m_map_from_var_index_to_column_info_with_cls.find(m_available_var_index) == m_map_from_var_index_to_column_info_with_cls.end());
unsigned vj = m_available_var_index;
m_map_from_var_index_to_column_info_with_cls[vj] = column_info_with_cls(left_side);
m_map_from_var_index_to_column_info_with_cls[vj].m_column_info.set_name("_s" + T_to_string(vj));
left_side->m_additional_var_index = vj;
m_available_var_index++;
} else { } else {
delete left_side; delete left_side;
left_side = *it; left_side = *it;
@ -43,14 +49,7 @@ canonic_left_side * lar_solver::create_or_fetch_existing_left_side(const buffer<
return left_side; return left_side;
} }
mpq lar_solver::find_ratio_of_original_constraint_to_normalized(canonic_left_side * ls, const lar_constraint & constraint) {
bool lar_solver::var_is_fixed(unsigned j) {
auto it = m_map_from_var_index_to_left_side.find(j);
lean_assert(it != m_map_from_var_index_to_left_side.end());
return it->second->m_column_info.is_fixed();
}
mpq lar_solver::find_ratio(canonic_left_side * ls, const lar_constraint & constraint) {
lean_assert(ls->m_coeffs.size() > 0); lean_assert(ls->m_coeffs.size() > 0);
auto first_pair = ls->m_coeffs[0]; auto first_pair = ls->m_coeffs[0];
lean_assert(first_pair.first == numeric_traits<mpq>::one()); lean_assert(first_pair.first == numeric_traits<mpq>::one());
@ -60,81 +59,68 @@ mpq lar_solver::find_ratio(canonic_left_side * ls, const lar_constraint & constr
return it->second; return it->second;
} }
void lar_solver::add_canonic_left_side_for_var(var_index i, std::string var_name) { void lar_solver::map_left_side_to_A_of_core_solver(canonic_left_side* left_side, unsigned j) {
buffer<std::pair<mpq, var_index>> b; var_index additional_var = left_side->m_additional_var_index;
b.push_back(std::make_pair(numeric_traits<mpq>::one(), i)); lean_assert(valid_index(additional_var));
auto can_ls = new canonic_left_side(b); auto it = m_map_from_var_index_to_column_info_with_cls.find(additional_var);
can_ls->set_name(var_name); lean_assert(it != m_map_from_var_index_to_column_info_with_cls.end());
column_info<mpq> & ci = it->second.m_column_info;
lean_assert(m_canonic_left_sides.find(can_ls) == m_canonic_left_sides.end()); lean_assert(!is_valid(ci.get_column_index()));
m_canonic_left_sides.insert(can_ls); lean_assert(left_side->size() > 0); // if size is zero we have an empty row
m_map_from_var_index_to_left_side[i] = can_ls; left_side->m_row_index = m_basis.size();
} m_basis.push_back(j); // j will be a basis column, so we put it into the basis as well
lean_assert(m_map_from_column_indices_to_var_index.find(j) == m_map_from_column_indices_to_var_index.end());
void lar_solver::map_left_side_to_column_of_A(canonic_left_side* left_side, unsigned & j) { ci.set_column_index(j);
lean_assert(m_column_indices_to_canonic_left_sides.find(j) == m_column_indices_to_canonic_left_sides.end()); m_map_from_column_indices_to_var_index[j] = additional_var;
left_side->m_column_index = j; // assigning this index does not change the hash of canonic_left_side
if (left_side->size() > 1) { // if size is one we will not create a row for this left side
left_side->m_row_index = static_cast<int>(m_basis.size());
m_basis.push_back(j); // j will be a basis column, so we put it into the basis as well
}
m_column_indices_to_canonic_left_sides[j++] = left_side; // pointing from the column to the left side
} }
void lar_solver::map_left_sides_to_columns_of_A() { void lar_solver::map_left_sides_to_A_of_core_solver() {
unsigned j = 0; unsigned j = m_map_from_column_indices_to_var_index.size();
for (auto it : m_canonic_left_sides) for (auto it : m_set_of_canonic_left_sides)
map_left_side_to_column_of_A(it, j); map_left_side_to_A_of_core_solver(it, j++);
} }
bool valid_index(unsigned j) { return static_cast<int>(j) >= 0;}
// this adds a row to A // this adds a row to A
template <typename U, typename V> template <typename U, typename V>
void lar_solver::add_row_to_A(static_matrix<U, V> & A, unsigned i, canonic_left_side * ls) { void lar_solver::fill_row_of_A(static_matrix<U, V> & A, unsigned i, canonic_left_side * ls) {
for (auto & t : ls->m_coeffs) { for (auto & t : ls->m_coeffs) {
var_index vi = t.second; var_index vi = t.second;
canonic_left_side *var_left_side = m_map_from_var_index_to_left_side[vi]; unsigned column = get_column_index_from_var_index(vi);
unsigned column = var_left_side->m_column_index; lean_assert(is_valid(column));
lean_assert(valid_index(column));
A.set(i, column, convert_struct<U, mpq>::convert(t.first)); A.set(i, column, convert_struct<U, mpq>::convert(t.first));
} }
A.set(i, ls->m_column_index, - one_of_type<U>()); unsigned additional_column = get_column_index_from_var_index(ls->m_additional_var_index);
lean_assert(is_valid(additional_column));
A.set(i, additional_column, - one_of_type<U>());
}
void lar_solver::fill_set_of_active_var_indices() {
for (auto & t : m_set_of_canonic_left_sides)
for (auto & a : t->m_coeffs)
m_set_of_active_var_indices.insert(a.second);
} }
template <typename U, typename V> template <typename U, typename V>
void lar_solver::create_matrix_A(static_matrix<U, V> & A) { void lar_solver::create_matrix_A(static_matrix<U, V> & A) {
unsigned n = static_cast<unsigned>(m_column_indices_to_canonic_left_sides.size()); unsigned m = m_set_of_canonic_left_sides.size();
unsigned m = static_cast<unsigned>(m_basis.size()); fill_set_of_active_var_indices();
unsigned n = m_set_of_active_var_indices.size() + m;
A.init_empty_matrix(m, n); A.init_empty_matrix(m, n);
unsigned i = 0; unsigned i = 0;
for (auto & t : m_column_indices_to_canonic_left_sides) for (auto t : m_set_of_canonic_left_sides) {
if (valid_index(t.second->m_row_index)) lean_assert(t->size() > 0);
add_row_to_A(A, i++, t.second); if (is_valid(t->m_row_index))
#ifdef LEAN_DEBUG fill_row_of_A(A, i++, t);
// print_matrix(m_A); }
#endif
} }
// void lar_solver::fill_column_info_names() {
// for (unsigned j = 0; j < m_A.column_count(); j++) {
// column_info<mpq> t;
// m_column_infos.push_back(t);
// if (j < m_map_from_var_index_to_name_left_side_pair.size()) {
// m_column_infos.back().set_name(m_map_from_var_index_to_name_left_side_pair[j]);
// } else {
// string pref("_s_");
// m_column_infos.back().set_name(pref + T_to_string(j));
// }
// m_column_names
// }
// }
void lar_solver::set_upper_bound_for_column_info(lar_normalized_constraint * norm_constr) { void lar_solver::set_upper_bound_for_column_info(lar_normalized_constraint * norm_constr) {
const mpq & v = norm_constr->m_right_side; const mpq & v = norm_constr->m_right_side;
canonic_left_side * ls = norm_constr->m_canonic_left_side; canonic_left_side * ls = norm_constr->m_canonic_left_side;
column_info<mpq> & ci = ls->m_column_info; var_index additional_var_index = ls->m_additional_var_index;
lean_assert(is_valid(additional_var_index));
column_info<mpq> & ci = get_column_info_from_var_index(additional_var_index);
lean_assert(norm_constr->m_kind == LE || norm_constr->m_kind == LT || norm_constr->m_kind == EQ); lean_assert(norm_constr->m_kind == LE || norm_constr->m_kind == LT || norm_constr->m_kind == EQ);
bool strict = norm_constr->m_kind == LT; bool strict = norm_constr->m_kind == LT;
if (!ci.upper_bound_is_set()) { if (!ci.upper_bound_is_set()) {
@ -165,7 +151,7 @@ bool lar_solver::try_to_set_fixed(column_info<mpq> & ci) {
void lar_solver::set_low_bound_for_column_info(lar_normalized_constraint * norm_constr) { void lar_solver::set_low_bound_for_column_info(lar_normalized_constraint * norm_constr) {
const mpq & v = norm_constr->m_right_side; const mpq & v = norm_constr->m_right_side;
canonic_left_side * ls = norm_constr->m_canonic_left_side; canonic_left_side * ls = norm_constr->m_canonic_left_side;
column_info<mpq> & ci = ls->m_column_info; column_info<mpq> & ci = get_column_info_from_var_index(ls->m_additional_var_index);
lean_assert(norm_constr->m_kind == GE || norm_constr->m_kind == GT || norm_constr->m_kind == EQ); lean_assert(norm_constr->m_kind == GE || norm_constr->m_kind == GT || norm_constr->m_kind == EQ);
bool strict = norm_constr->m_kind == GT; bool strict = norm_constr->m_kind == GT;
if (!ci.low_bound_is_set()) { if (!ci.low_bound_is_set()) {
@ -221,10 +207,10 @@ column_type lar_solver::get_column_type(column_info<mpq> & ci) {
void lar_solver::fill_column_names() { void lar_solver::fill_column_names() {
m_column_names.clear(); m_column_names.clear();
for (auto t : m_canonic_left_sides) { for (auto & t : m_map_from_var_index_to_column_info_with_cls) {
auto & ci = t->m_column_info; column_info<mpq> & ci = t.second.m_column_info;
unsigned j = t->m_column_index; unsigned j = ci.get_column_index();
lean_assert(valid_index(j)); lean_assert(is_valid(j));
std::string name = ci.get_name(); std::string name = ci.get_name();
if (name.size() == 0) if (name.size() == 0)
name = std::string("_s") + T_to_string(j); name = std::string("_s") + T_to_string(j);
@ -234,24 +220,25 @@ void lar_solver::fill_column_names() {
void lar_solver::fill_column_types() { void lar_solver::fill_column_types() {
m_column_types.clear(); m_column_types.clear();
m_column_types.resize(m_canonic_left_sides.size(), free_column); m_column_types.resize(m_map_from_var_index_to_column_info_with_cls.size(), free_column);
for (auto t : m_canonic_left_sides) { for (auto t : m_set_of_canonic_left_sides) {
auto & ci = t->m_column_info; var_index additional_vj = t->m_additional_var_index;
unsigned j = t->m_column_index; unsigned j = get_column_index_from_var_index(additional_vj);
lean_assert(valid_index(j)); lean_assert(is_valid(j));
m_column_types[j] = get_column_type(ci); m_column_types[j] = get_column_type(get_column_info_from_var_index(additional_vj));
} }
} }
template <typename V> template <typename V>
void lar_solver::fill_bounds_for_core_solver(std::vector<V> & lb, std::vector<V> & ub) { void lar_solver::fill_bounds_for_core_solver(std::vector<V> & lb, std::vector<V> & ub) {
unsigned n = static_cast<unsigned>(m_canonic_left_sides.size()); // this is the number of columns unsigned n = static_cast<unsigned>(m_map_from_var_index_to_column_info_with_cls.size()); // this is the number of columns
lb.resize(n); lb.resize(n);
ub.resize(n); ub.resize(n);
for (auto t : m_canonic_left_sides) { for (auto t : m_set_of_canonic_left_sides) {
auto & ci = t->m_column_info; auto & ci = get_column_info_from_var_index(t->m_additional_var_index);
unsigned j = t->m_column_index; unsigned j = ci.get_column_index();
lean_assert(valid_index(j)); lean_assert(is_valid(j));
lean_assert(j < n);
if (ci.low_bound_is_set()) if (ci.low_bound_is_set())
lb[j] = conversion_helper<V>::get_low_bound(ci); lb[j] = conversion_helper<V>::get_low_bound(ci);
if (ci.upper_bound_is_set()) if (ci.upper_bound_is_set())
@ -290,31 +277,32 @@ template <typename V> V lar_solver::get_column_val(std::vector<V> & low_bound, s
lar_solver::~lar_solver() { lar_solver::~lar_solver() {
std::vector<canonic_left_side*> to_delete; std::vector<canonic_left_side*> to_delete;
for (auto it : m_canonic_left_sides) for (auto it : m_set_of_canonic_left_sides)
to_delete.push_back(it); to_delete.push_back(it);
for (auto t : to_delete) for (auto t : to_delete)
delete t; delete t;
} }
var_index lar_solver::add_var(std::string s) { var_index lar_solver::add_var(std::string s) {
auto got = m_var_names_to_var_index.find(s); auto it = m_var_names_to_var_index.find(s);
if (got != m_var_names_to_var_index.end()) return got->second; if (it != m_var_names_to_var_index.end())
return it->second;
var_index i = m_available_var_index++; var_index i = m_available_var_index++;
lean_assert(m_map_from_var_index_to_column_info_with_cls.find(i) == m_map_from_var_index_to_column_info_with_cls.end());
auto ci_with_cls = column_info_with_cls();
ci_with_cls.m_column_info.set_name(s);
m_map_from_var_index_to_column_info_with_cls[i] = ci_with_cls;
m_var_names_to_var_index[s] = i; m_var_names_to_var_index[s] = i;
add_canonic_left_side_for_var(i, s);
return i; return i;
} }
constraint_index lar_solver::add_constraint(const buffer<std::pair<mpq, var_index>>& left_side, lconstraint_kind kind_par, mpq right_side_par) { constraint_index lar_solver::add_constraint(const buffer<std::pair<mpq, var_index>>& left_side, lconstraint_kind kind_par, mpq right_side_par) {
lean_assert(left_side.size() > 0); lean_assert(left_side.size() > 0);
constraint_index i = m_available_constr_index++; constraint_index i = m_available_constr_index++;
lean_assert(m_normalized_constraints.find(i) == m_normalized_constraints.end());
lar_constraint original_constr(left_side, kind_par, right_side_par, i); lar_constraint original_constr(left_side, kind_par, right_side_par, i);
canonic_left_side * ls = create_or_fetch_existing_left_side(left_side); canonic_left_side * ls = create_or_fetch_existing_left_side(left_side);
mpq ratio = find_ratio(ls, original_constr); mpq ratio = find_ratio_of_original_constraint_to_normalized(ls, original_constr);
auto kind = ratio.is_neg()? flip_kind(kind_par): kind_par; auto kind = ratio.is_neg()? flip_kind(kind_par): kind_par;
mpq right_side = right_side_par / ratio; mpq right_side = right_side_par / ratio;
lar_normalized_constraint normalized_constraint(ls, ratio, kind, right_side, original_constr); lar_normalized_constraint normalized_constraint(ls, ratio, kind, right_side, original_constr);
@ -455,9 +443,21 @@ mpq lar_solver::sum_of_right_sides_of_evidence(const buffer<std::pair<mpq, unsig
} }
return ret; return ret;
} }
// do not touch additional vars here
void lar_solver::map_var_indices_to_columns_of_A() {
int i = 0;
for (auto & it : m_map_from_var_index_to_column_info_with_cls) {
if (it.second.m_canonic_left_side != nullptr) continue;
lean_assert(m_map_from_column_indices_to_var_index.find(i) == m_map_from_column_indices_to_var_index.end());
it.second.m_column_info.set_column_index(i);
m_map_from_column_indices_to_var_index[i] = it.first;
i++;
}
}
void lar_solver::prepare_independently_of_numeric_type() { void lar_solver::prepare_independently_of_numeric_type() {
update_column_info_of_normalized_constraints(); update_column_info_of_normalized_constraints();
map_left_sides_to_columns_of_A(); map_var_indices_to_columns_of_A();
map_left_sides_to_A_of_core_solver();
fill_column_names(); fill_column_names();
fill_column_types(); fill_column_types();
} }
@ -565,7 +565,12 @@ void lar_solver::get_infeasibility_evidence_for_inf_sign(buffer<std::pair<mpq, c
for (auto & it : inf_row) { for (auto & it : inf_row) {
mpq coeff = it.first; mpq coeff = it.first;
unsigned j = it.second; unsigned j = it.second;
canonic_left_side * ls = m_column_indices_to_canonic_left_sides[j]; auto it1 = m_map_from_column_indices_to_var_index.find(j);
lean_assert(it1 != m_map_from_column_indices_to_var_index.end());
unsigned var_j = it1->second;
auto it2 = m_map_from_var_index_to_column_info_with_cls.find(var_j);
lean_assert(it2 != m_map_from_var_index_to_column_info_with_cls.end());
canonic_left_side * ls = it2->second.m_canonic_left_side;
int adj_sign = coeff.is_pos() ? inf_sign : -inf_sign; int adj_sign = coeff.is_pos() ? inf_sign : -inf_sign;
lar_normalized_constraint * bound_constr = adj_sign < 0? ls->m_upper_bound_witness : ls->m_low_bound_witness; lar_normalized_constraint * bound_constr = adj_sign < 0? ls->m_upper_bound_witness : ls->m_low_bound_witness;
@ -577,10 +582,10 @@ void lar_solver::get_infeasibility_evidence_for_inf_sign(buffer<std::pair<mpq, c
mpq lar_solver::find_delta_for_strict_bounds() { mpq lar_solver::find_delta_for_strict_bounds() {
mpq delta = numeric_traits<mpq>::one(); mpq delta = numeric_traits<mpq>::one();
for (auto t : m_canonic_left_sides) { for (auto t : m_set_of_canonic_left_sides) {
auto & ci = t->m_column_info; auto & ci = get_column_info_from_var_index(t->m_additional_var_index);
unsigned j = t->m_column_index; unsigned j = ci.get_column_index();
lean_assert (valid_index(j)); lean_assert (is_valid(j));
if (ci.low_bound_is_set()) if (ci.low_bound_is_set())
restrict_delta_on_low_bound_column(delta, j); restrict_delta_on_low_bound_column(delta, j);
if (ci.upper_bound_is_set()) if (ci.upper_bound_is_set())
@ -623,18 +628,21 @@ void lar_solver::restrict_delta_on_upper_bound(mpq& delta, unsigned j) {
void lar_solver::get_model(std::unordered_map<var_index, mpq> & variable_values){ void lar_solver::get_model(std::unordered_map<var_index, mpq> & variable_values){
lean_assert(m_status == OPTIMAL); lean_assert(m_status == OPTIMAL);
mpq delta = find_delta_for_strict_bounds(); mpq delta = find_delta_for_strict_bounds();
for (auto & it : m_map_from_var_index_to_left_side) { for (auto & it : m_map_from_var_index_to_column_info_with_cls) {
numeric_pair<mpq> & rp = m_x[it.second->m_column_index]; column_info<mpq> & ci = it.second.m_column_info;
unsigned j = ci.get_column_index();
numeric_pair<mpq> & rp = m_x[j];
variable_values[it.first] = rp.x + delta * rp.y; variable_values[it.first] = rp.x + delta * rp.y;
} }
} }
std::string lar_solver::get_variable_name(var_index vi) { std::string lar_solver::get_variable_name(var_index vi) const {
if (m_map_from_var_index_to_left_side.size() <= vi) { auto it = m_map_from_var_index_to_column_info_with_cls.find(vi);
if (it == m_map_from_var_index_to_column_info_with_cls.end()) {
std::string s = "variable " + T_to_string(vi) + " is not found"; std::string s = "variable " + T_to_string(vi) + " is not found";
return s; return s;
} }
return m_map_from_var_index_to_left_side[vi]->m_column_info.get_name(); return it->second.m_column_info.get_name();
} }
// ********** print region start // ********** print region start
@ -650,7 +658,7 @@ void lar_solver::print_constraint(constraint_index ci, std::ostream & out) {
void lar_solver::print_canonic_left_side(const canonic_left_side & c, std::ostream & out) { void lar_solver::print_canonic_left_side(const canonic_left_side & c, std::ostream & out) {
bool first = true; bool first = true;
for (auto it : c.m_coeffs) { for (auto & it : c.m_coeffs) {
auto val = it.first; auto val = it.first;
if (first) { if (first) {
first = false; first = false;
@ -664,13 +672,13 @@ void lar_solver::print_canonic_left_side(const canonic_left_side & c, std::ostre
} }
if (val != numeric_traits<mpq>::one()) if (val != numeric_traits<mpq>::one())
out << T_to_string(val); out << T_to_string(val);
out << m_map_from_var_index_to_left_side[it.second]->m_column_info.get_name(); out << get_variable_name(it.second);
} }
} }
void lar_solver::print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) { void lar_solver::print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) {
bool first = true; bool first = true;
for (auto it : c->get_left_side_coefficients()) { for (auto & it : c->get_left_side_coefficients()) {
auto val = it.first; auto val = it.first;
if (numeric_traits<mpq>::is_zero(val)) continue; if (numeric_traits<mpq>::is_zero(val)) continue;
if (first) { if (first) {
@ -686,29 +694,29 @@ void lar_solver::print_left_side_of_constraint(const lar_base_constraint * c, st
if (val != numeric_traits<mpq>::one()) if (val != numeric_traits<mpq>::one())
out << val; out << val;
out << m_map_from_var_index_to_left_side[it.second]->m_column_info.get_name(); out << get_variable_name(it.second);
} }
} }
void lar_solver::print_info_on_column(unsigned j, std::ostream & out) { // void lar_solver::print_info_on_column(unsigned j, std::ostream & out) {
for (auto ls : m_canonic_left_sides) { // for (auto ls : m_set_of_canonic_left_sides) {
if (static_cast<unsigned>(ls->m_column_index) == j) { // if (static_cast<unsigned>(ls->get_ci()) == j) {
auto & ci = ls->m_column_info; // auto & ci = ls->m_column_info;
if (ci.low_bound_is_set()) { // if (ci.low_bound_is_set()) {
out << "l = " << ci.get_low_bound(); // out << "l = " << ci.get_low_bound();
} // }
if (ci.upper_bound_is_set()) { // if (ci.upper_bound_is_set()) {
out << "u = " << ci.get_upper_bound(); // out << "u = " << ci.get_upper_bound();
} // }
out << std::endl; // out << std::endl;
m_mpq_core_solver.print_column_info(j, out); // m_mpq_core_solver.print_column_info(j, out);
} // }
} // }
} // }
mpq lar_solver::get_infeasibility_of_solution(std::unordered_map<std::string, mpq> & solution) { mpq lar_solver::get_infeasibility_of_solution(std::unordered_map<std::string, mpq> & solution) {
mpq ret = numeric_traits<mpq>::zero(); mpq ret = numeric_traits<mpq>::zero();
for (auto it : m_normalized_constraints) { for (auto & it : m_normalized_constraints) {
ret += get_infeasibility_of_constraint(it.second, solution); ret += get_infeasibility_of_constraint(it.second, solution);
} }
return ret; return ret;
@ -735,12 +743,9 @@ mpq lar_solver::get_infeasibility_of_constraint(const lar_normalized_constraint
mpq lar_solver::get_canonic_left_side_val(canonic_left_side * ls, std::unordered_map<std::string, mpq> & solution) { mpq lar_solver::get_canonic_left_side_val(canonic_left_side * ls, std::unordered_map<std::string, mpq> & solution) {
mpq ret = numeric_traits<mpq>::zero(); mpq ret = numeric_traits<mpq>::zero();
for (auto it : ls->m_coeffs) { for (auto & it : ls->m_coeffs) {
var_index j = it.second; var_index vi = it.second;
auto vi = m_map_from_var_index_to_left_side.find(j); std::string s = get_variable_name(vi);
lean_assert(vi != m_map_from_var_index_to_left_side.end());
canonic_left_side * var_ls = vi->second;
std::string s = var_ls->m_column_info.get_name();
auto t = solution.find(s); auto t = solution.find(s);
lean_assert(t != solution.end()); lean_assert(t != solution.end());
ret += it.first * (t->second); ret += it.first * (t->second);
@ -750,7 +755,7 @@ mpq lar_solver::get_canonic_left_side_val(canonic_left_side * ls, std::unordered
mpq lar_solver::get_left_side_val(const lar_constraint & cns, const std::unordered_map<var_index, mpq> & var_map) { mpq lar_solver::get_left_side_val(const lar_constraint & cns, const std::unordered_map<var_index, mpq> & var_map) {
mpq ret = numeric_traits<mpq>::zero(); mpq ret = numeric_traits<mpq>::zero();
for (auto it : cns.m_left_side) { for (auto & it : cns.m_left_side) {
var_index j = it.first; var_index j = it.first;
auto vi = var_map.find(j); auto vi = var_map.find(j);
lean_assert(vi != var_map.end()); lean_assert(vi != var_map.end());
@ -763,5 +768,21 @@ void lar_solver::print_constraint(const lar_base_constraint * c, std::ostream &
print_left_side_of_constraint(c, out); print_left_side_of_constraint(c, out);
out <<" " << lconstraint_kind_string(c->m_kind) << " " << c->m_right_side; out <<" " << lconstraint_kind_string(c->m_kind) << " " << c->m_right_side;
} }
unsigned lar_solver::get_column_index_from_var_index(var_index vi) const{
auto it = m_map_from_var_index_to_column_info_with_cls.find(vi);
if (it == m_map_from_var_index_to_column_info_with_cls.end())
return static_cast<unsigned>(-1);
return it->second.m_column_info.get_column_index();
}
column_info<mpq> & lar_solver::get_column_info_from_var_index(var_index vi) {
auto it = m_map_from_var_index_to_column_info_with_cls.find(vi);
if (it == m_map_from_var_index_to_column_info_with_cls.end()) {
lean_assert(false);
throw 0;
}
return it->second.m_column_info;
}
} }

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

@ -33,6 +33,14 @@ struct conversion_helper {
} }
}; };
struct column_info_with_cls { // column_info_with canonic_left_side
canonic_left_side * m_canonic_left_side;
column_info<mpq> m_column_info;
column_info_with_cls(): m_canonic_left_side(nullptr), m_column_info(static_cast<unsigned>(-1)) {}
column_info_with_cls(canonic_left_side * cls) : m_canonic_left_side(cls), m_column_info(static_cast<unsigned>(-1)) {}
};
template<> template<>
struct conversion_helper <double> { struct conversion_helper <double> {
static double get_low_bound(const column_info<mpq> & ci); static double get_low_bound(const column_info<mpq> & ci);
@ -43,11 +51,12 @@ class lar_solver {
unsigned m_available_var_index = 0; unsigned m_available_var_index = 0;
unsigned m_available_constr_index = 0; unsigned m_available_constr_index = 0;
lp_status m_status = UNKNOWN; lp_status m_status = UNKNOWN;
std::unordered_set<var_index> m_set_of_active_var_indices; // a variable is active if it is referenced in a left side
std::unordered_map<std::string, var_index> m_var_names_to_var_index; std::unordered_map<std::string, var_index> m_var_names_to_var_index;
std::unordered_map<var_index, canonic_left_side*> m_map_from_var_index_to_left_side; std::unordered_set<canonic_left_side*, hash_and_equal_of_canonic_left_side_struct, hash_and_equal_of_canonic_left_side_struct> m_set_of_canonic_left_sides;
std::unordered_set<canonic_left_side*, hash_and_equal_of_canonic_left_side_struct, hash_and_equal_of_canonic_left_side_struct> m_canonic_left_sides; std::unordered_map<unsigned, var_index> m_map_from_column_indices_to_var_index;
std::unordered_map<unsigned, canonic_left_side*> m_column_indices_to_canonic_left_sides;
std::unordered_map<constraint_index, lar_normalized_constraint> m_normalized_constraints; std::unordered_map<constraint_index, lar_normalized_constraint> m_normalized_constraints;
std::unordered_map<var_index, column_info_with_cls> m_map_from_var_index_to_column_info_with_cls;
static_matrix<mpq, numeric_pair<mpq>> m_A; static_matrix<mpq, numeric_pair<mpq>> m_A;
lar_core_solver<mpq, numeric_pair<mpq>> m_mpq_core_solver; lar_core_solver<mpq, numeric_pair<mpq>> m_mpq_core_solver;
lp_settings m_settings; lp_settings m_settings;
@ -62,24 +71,20 @@ class lar_solver {
canonic_left_side * m_infeasible_canonic_left_side = nullptr; // such can be found at the initialization step canonic_left_side * m_infeasible_canonic_left_side = nullptr; // such can be found at the initialization step
canonic_left_side * create_or_fetch_existing_left_side(const buffer<std::pair<mpq, var_index>>& left_side_par); canonic_left_side * create_or_fetch_existing_left_side(const buffer<std::pair<mpq, var_index>>& left_side_par);
mpq find_ratio_of_original_constraint_to_normalized(canonic_left_side * ls, const lar_constraint & constraint);
bool var_is_fixed(unsigned j);
mpq find_ratio(canonic_left_side * ls, const lar_constraint & constraint);
void add_canonic_left_side_for_var(var_index i, std::string var_name); void add_canonic_left_side_for_var(var_index i, std::string var_name);
void map_left_side_to_column_of_A(canonic_left_side* left_side, unsigned & j); void map_left_side_to_A_of_core_solver(canonic_left_side* left_side, var_index vj);
void map_left_sides_to_A_of_core_solver();
void map_left_sides_to_columns_of_A();
bool valid_index(unsigned j) { return static_cast<int>(j) >= 0;} bool valid_index(unsigned j) { return static_cast<int>(j) >= 0;}
// this adds a row to A // this adds a row to A
template <typename U, typename V> template <typename U, typename V>
void add_row_to_A(static_matrix<U, V> & A, unsigned i, canonic_left_side * ls); void fill_row_of_A(static_matrix<U, V> & A, unsigned i, canonic_left_side * ls);
template <typename U, typename V> template <typename U, typename V>
void create_matrix_A(static_matrix<U, V> & A); void create_matrix_A(static_matrix<U, V> & A);
@ -123,6 +128,12 @@ class lar_solver {
std::vector<V> & upper_bound, const lar_solution_signature & signature); std::vector<V> & upper_bound, const lar_solution_signature & signature);
template <typename V> V get_column_val(std::vector<V> & low_bound, std::vector<V> & upper_bound, non_basic_column_value_position pos_type, unsigned j); template <typename V> V get_column_val(std::vector<V> & low_bound, std::vector<V> & upper_bound, non_basic_column_value_position pos_type, unsigned j);
void map_var_indices_to_columns_of_A();
void register_in_map(std::unordered_map<var_index, mpq> & coeffs, lar_constraint & cn, const mpq & a);
unsigned get_column_index_from_var_index(var_index vi) const;
column_info<mpq> & get_column_info_from_var_index(var_index vi);
void fill_set_of_active_var_indices();
public: public:
~lar_solver(); ~lar_solver();
@ -149,10 +160,6 @@ public:
constraint_index add_constraint(const buffer<std::pair<mpq, var_index>>& left_side, lconstraint_kind kind_par, mpq right_side_par); constraint_index add_constraint(const buffer<std::pair<mpq, var_index>>& left_side, lconstraint_kind kind_par, mpq right_side_par);
bool is_infeasible(const column_info<mpq> & ci) {
return ci.low_bound_is_set() && ci.upper_bound_is_set() && ci.get_low_bound() > ci.get_upper_bound();
}
bool all_constraints_hold(); bool all_constraints_hold();
bool constraint_holds(const lar_constraint & constr, std::unordered_map<var_index, mpq> & var_map); bool constraint_holds(const lar_constraint & constr, std::unordered_map<var_index, mpq> & var_map);
@ -163,7 +170,6 @@ public:
bool the_relations_are_of_same_type(const buffer<std::pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum); bool the_relations_are_of_same_type(const buffer<std::pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum);
void register_in_map(std::unordered_map<var_index, mpq> & coeffs, lar_constraint & cn, const mpq & a);
bool the_left_sides_sum_to_zero(const buffer<std::pair<mpq, unsigned>> & evidence); bool the_left_sides_sum_to_zero(const buffer<std::pair<mpq, unsigned>> & evidence);
bool the_righ_sides_do_not_sum_to_zero(const buffer<std::pair<mpq, unsigned>> & evidence); bool the_righ_sides_do_not_sum_to_zero(const buffer<std::pair<mpq, unsigned>> & evidence);
@ -214,7 +220,7 @@ public:
void get_model(std::unordered_map<var_index, mpq> & variable_values); void get_model(std::unordered_map<var_index, mpq> & variable_values);
std::string get_variable_name(var_index vi); std::string get_variable_name(var_index vi) const;
void print_constraint(constraint_index ci, std::ostream & out); void print_constraint(constraint_index ci, std::ostream & out);
@ -224,7 +230,7 @@ public:
numeric_pair<mpq> get_infeasibility_from_core_solver(std::unordered_map<std::string, mpq> & solution); numeric_pair<mpq> get_infeasibility_from_core_solver(std::unordered_map<std::string, mpq> & solution);
void print_info_on_column(unsigned j, std::ostream & out); // void print_info_on_column(unsigned j, std::ostream & out);
mpq get_infeasibility_of_solution(std::unordered_map<std::string, mpq> & solution); mpq get_infeasibility_of_solution(std::unordered_map<std::string, mpq> & solution);

6
src/util/lp/lp_primal_simplex.cpp
View file

@ -145,9 +145,9 @@ template <typename T, typename X> void lp_primal_simplex<T, X>::set_core_solver_
this->m_upper_bounds.resize(total_vars); this->m_upper_bounds.resize(total_vars);
for (auto cit : this->m_columns) { for (auto cit : this->m_columns) {
column_info<T> * ci = cit.second; column_info<T> * ci = cit.second;
auto p = this->m_external_columns_to_core_solver_columns.find(cit.first); unsigned j = ci->get_column_index();
if (p == this->m_external_columns_to_core_solver_columns.end()) continue; if (!is_valid(j))
unsigned j = p->second; continue; // the variable is not mapped to a column
switch (this->m_column_types[j] = ci->get_column_type()){ switch (this->m_column_types[j] = ci->get_column_type()){
case fixed: case fixed:
this->m_upper_bounds[j] = numeric_traits<T>::zero(); this->m_upper_bounds[j] = numeric_traits<T>::zero();

25
src/util/lp/lp_solver.cpp
View file

@ -401,9 +401,9 @@ template <typename T, typename X> void lp_solver<T, X>::map_external_columns_to_
throw_exception("found fixed column"); throw_exception("found fixed column");
} }
unsigned j = col.first; unsigned j = col.first;
auto j_place = m_external_columns_to_core_solver_columns.find(j); auto j_column = m_columns[j]->get_column_index();
if (j_place == m_external_columns_to_core_solver_columns.end()) { // j is a newcomer if (!is_valid(j_column)) { // j is a newcomer
m_external_columns_to_core_solver_columns[j] = size; m_columns[j]->set_column_index(size);
m_core_solver_columns_to_external_columns[size++] = j; m_core_solver_columns_to_external_columns[size++] = j;
} }
} }
@ -412,10 +412,9 @@ template <typename T, typename X> void lp_solver<T, X>::map_external_columns_to_
template <typename T, typename X> void lp_solver<T, X>::fill_column_names_for_core_solver() { template <typename T, typename X> void lp_solver<T, X>::fill_column_names_for_core_solver() {
for (auto it : this->m_columns) { for (auto it : this->m_columns) {
auto p = this->m_external_columns_to_core_solver_columns.find(it.first); unsigned j = it.second->get_column_index();
if (p != this->m_external_columns_to_core_solver_columns.end()) { if (is_valid(j))
this->m_name_map[p->second] = it.second->get_name(); this->m_name_map[j] = it.second->get_name();
}
} }
} }
@ -434,10 +433,9 @@ template <typename T, typename X> void lp_solver<T, X>::fill_A_from_A_values() {
lean_assert(m_external_rows_to_core_solver_rows.find(t.first) != m_external_rows_to_core_solver_rows.end()); lean_assert(m_external_rows_to_core_solver_rows.find(t.first) != m_external_rows_to_core_solver_rows.end());
unsigned row = m_external_rows_to_core_solver_rows[t.first]; unsigned row = m_external_rows_to_core_solver_rows[t.first];
for (auto k : t.second) { for (auto k : t.second) {
lean_assert(m_external_columns_to_core_solver_columns.find(k.first) != m_external_columns_to_core_solver_columns.end()); unsigned col = m_columns[k.first]->get_column_index();
unsigned col = m_external_columns_to_core_solver_columns[k.first]; lean_assert(is_valid(col));
bool col_is_flipped = m_columns[k.first]->is_flipped(); bool col_is_flipped = m_columns[k.first]->is_flipped();
if (!col_is_flipped) { if (!col_is_flipped) {
(*m_A)(row, col) = k.second; (*m_A)(row, col) = k.second;
} else { } else {
@ -505,7 +503,7 @@ template <typename T, typename X> void lp_solver<T, X>::fill_m_b() {
} }
} }
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 { 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_columns.find(column); auto cit = this->m_columns.find(column);
if (cit == this->m_columns.end()) { if (cit == this->m_columns.end()) {
return numeric_traits<T>::zero(); return numeric_traits<T>::zero();
@ -517,9 +515,8 @@ template <typename T, typename X> T lp_solver<T, X>::get_column_value_with_co
return ci->get_fixed_value(); return ci->get_fixed_value();
} }
auto t = this->m_external_columns_to_core_solver_columns.find(column); unsigned cj = ci->get_column_index();
if (t != this->m_external_columns_to_core_solver_columns.end()){ if (cj != static_cast<unsigned>(-1)) {
unsigned cj = t->second;
T v = core_solver->get_var_value(cj) * this->m_column_scale[cj]; T v = core_solver->get_var_value(cj) * this->m_column_scale[cj];
if (ci->is_free()) { if (ci->is_free()) {
return v; return v;

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

@ -43,12 +43,11 @@ public:
unsigned m_first_stage_iterations = 0; unsigned m_first_stage_iterations = 0;
unsigned m_second_stage_iterations = 0; unsigned m_second_stage_iterations = 0;
std::unordered_map<unsigned, lp_constraint<T, X>> m_constraints; std::unordered_map<unsigned, lp_constraint<T, X>> m_constraints;
std::unordered_map<unsigned, column_info<T>*> m_columns; std::unordered_map<var_index, column_info<T>*> m_columns;
std::unordered_map<unsigned, std::unordered_map<unsigned, T> > m_A_values; 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<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_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_rows_to_external_rows;
std::unordered_map<unsigned, unsigned> m_external_columns_to_core_solver_columns;
std::unordered_map<unsigned, unsigned> m_core_solver_columns_to_external_columns; std::unordered_map<unsigned, unsigned> m_core_solver_columns_to_external_columns;
std::vector<T> m_column_scale; std::vector<T> m_column_scale;
std::unordered_map<unsigned, std::string> m_name_map; std::unordered_map<unsigned, std::string> m_name_map;
@ -197,7 +196,9 @@ protected:
void fill_column_names_for_core_solver(); void fill_column_names_for_core_solver();
unsigned number_of_core_structurals() { return static_cast<unsigned>(m_external_columns_to_core_solver_columns.size()); } unsigned number_of_core_structurals() {
return static_cast<unsigned>(m_core_solver_columns_to_external_columns.size());
}
void restore_column_scales_to_one() { void restore_column_scales_to_one() {

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

@ -18,7 +18,7 @@ namespace lean {
inline void lean_assert(bool b) {} inline void lean_assert(bool b) {}
#else #else
#define lean_assert(_x_) {} #define lean_assert(_x_) {}
#endif #endif
inline void lean_unreachable() { lean_assert(false); } inline void lean_unreachable() { lean_assert(false); }
template <typename X> inline X zero_of_type() { return numeric_traits<X>::zero(); } template <typename X> inline X zero_of_type() { return numeric_traits<X>::zero(); }
template <typename X> inline X one_of_type() { return numeric_traits<X>::one(); } template <typename X> inline X one_of_type() { return numeric_traits<X>::one(); }
@ -75,7 +75,7 @@ template <typename X> inline bool is_zero(const X & v) { return numeric_traits<X
template <typename X> inline double get_double(const X & v) { return numeric_traits<X>::get_double(v); } template <typename X> inline double get_double(const X & v) { return numeric_traits<X>::get_double(v); }
template <typename T> inline T zero_of_type() {return numeric_traits<T>::zero();} template <typename T> inline T zero_of_type() {return numeric_traits<T>::zero();}
inline void throw_exception(std::string str) { throw exception(str); } inline void throw_exception(std::string str) { throw exception(str); }
template <typename T> inline T from_string(std::string const & str) { lean_unreachable();} template <typename T> inline T from_string(std::string const & ) { lean_unreachable();}
template <> double inline from_string<double>(std::string const & str) { return atof(str.c_str());} template <> double inline from_string<double>(std::string const & str) { return atof(str.c_str());}
template <> mpq inline from_string<mpq>(std::string const & str) { template <> mpq inline from_string<mpq>(std::string const & str) {
return mpq(atof(str.c_str())); return mpq(atof(str.c_str()));