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

chore(util/lp): no "using", indentation

Browse source
This commit is contained in:
Leonardo de Moura 2016-01-08 15:50:24 -08:00
parent 1841d17544
commit 8fbf66d01a
28 changed files with 3070 additions and 3316 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

297
src/tests/util/lp/lp.cpp
View file

@ -4,35 +4,34 @@ Released under Apache 2.0 license as described in the file LICENSE.
Author: Lev Nachmanson
*/
#include <dirent.h>
#include <algorithm>
#include <string>
#include <set>
#include <iostream>
#include <sys/types.h>
#include <sys/stat.h>
#include <cstdlib>
#include "util/lp/lp.h"
#include "util/lp/lp_primal_simplex.h"
#include <ctime>
#include <vector>
#include <stdlib.h>
#include <utility>
#include "util/pair.h"
#include "util/lp/lp.h"
#include "util/lp/lp_primal_simplex.h"
#include "tests/util/lp/mps_reader.h"
#include "tests/util/lp/smt_reader.h"
#include "util/numerics/mpq.h"
#include "util/lp/binary_heap_priority_queue.h"
#include "tests/util/lp/argument_parser.h"
#include "tests/util/lp/test_file_reader.h"
#include <string>
#include <set>
#include "util/lp/indexed_value.h"
#include <stdlib.h>
#include "tests/util/lp/init_module.h"
#include "util/numerics/init_module.h"
#include <dirent.h>
#include <algorithm>
#include "util/lp/lar_solver.h"
#include "util/lp/numeric_pair.h"
#include "util/lp/binary_heap_upair_queue.h"
#include <utility>
using namespace lean;
using namespace std;
unsigned seed = 1;
#ifdef LEAN_DEBUG
unsigned lp_settings::ddd = 0;
@ -136,7 +135,7 @@ void tst1() {
float test = m10by9(0, 1);
std::cout << "got " << test << std:: endl;
std::cout << "got " << test << std::endl;
m10by9.set(0, 8, 8);
@ -189,39 +188,39 @@ void test_small_lu(lp_settings & settings) {
lu<double, double> l(m, basis, heading, settings, non_basic_columns);
lean_assert(l.is_correct());
indexed_vector<double> w(m.row_count());
cout << "entering 2, leaving 0" << endl;
cout << "entering 2, leaving 0" << std::endl;
l.prepare_entering(2, w); // to init vector w
l.replace_column(0, 0, w);
l.change_basis(2, 0);
// #ifdef LEAN_DEBUG
// cout << "we were factoring " << endl;
// cout << "we were factoring " << std::endl;
// print_matrix(get_B(l));
// #endif
lean_assert(l.is_correct());
cout << "entering 4, leaving 3" << endl;
cout << "entering 4, leaving 3" << std::endl;
l.prepare_entering(4, w); // to init vector w
l.replace_column(3, 0, w);
l.change_basis(4, 3);
cout << "we were factoring " << endl;
cout << "we were factoring " << std::endl;
#ifdef LEAN_DEBUG
print_matrix(get_B(l));
#endif
lean_assert(l.is_correct());
cout << "entering 5, leaving 1" << endl;
cout << "entering 5, leaving 1" << std::endl;
l.prepare_entering(5, w); // to init vector w
l.replace_column(1, 0, w);
l.change_basis(5, 1);
cout << "we were factoring " << endl;
cout << "we were factoring " << std::endl;
#ifdef LEAN_DEBUG
print_matrix(get_B(l));
#endif
lean_assert(l.is_correct());
cout << "entering 3, leaving 2" << endl;
cout << "entering 3, leaving 2" << std::endl;
l.prepare_entering(3, w); // to init vector w
l.replace_column(2, 0, w);
l.change_basis(3, 2);
cout << "we were factoring " << endl;
cout << "we were factoring " << std::endl;
#ifdef LEAN_DEBUG
print_matrix(get_B(l));
#endif
@ -400,13 +399,13 @@ void test_larger_lu(lp_settings& settings) {
dense_matrix<double, double> left_side = l.get_left_side();
dense_matrix<double, double> right_side = l.get_right_side();
if (!(left_side == right_side)) {
cout << "left side" << endl;
cout << "left side" << std::endl;
print_matrix(left_side);
cout << "right side" << endl;
cout << "right side" << std::endl;
print_matrix(right_side);
std::cout << "different sides" << std::endl;
cout << "initial factorization is incorrect" << endl;
cout << "initial factorization is incorrect" << std::endl;
exit(1);
}
indexed_vector<double> w(m.row_count());
@ -532,7 +531,7 @@ void test_lp_primal_core_solver() {
#ifdef LEAN_DEBUG
template <typename T, typename X>
void test_swap_rows_with_permutation(sparse_matrix<T, X>& m){
cout << "testing swaps" << endl;
cout << "testing swaps" << std::endl;
unsigned dim = m.row_count();
dense_matrix<double, double> original(m);
permutation_matrix<double, double> q(dim);
@ -542,12 +541,12 @@ void test_swap_rows_with_permutation(sparse_matrix<T, X>& m){
unsigned row1 = lrand48() % dim;
unsigned row2 = lrand48() % dim;
if (row1 == row2) continue;
cout << "swap " << row1 << " " << row2 << endl;
cout << "swap " << row1 << " " << row2 << std::endl;
m.swap_rows(row1, row2);
q.transpose_from_left(row1, row2);
lean_assert(original == q * m);
print_matrix(m);
cout << endl;
cout << std::endl;
}
}
#endif
@ -579,14 +578,14 @@ void matrix_repro_test() {
auto l = d * m.m_column_permutation;
l = m.m_row_permutation * l;
cout << "matrix_repro_test" << endl;
cout << "matrix_repro_test" << std::endl;
print_matrix(l);
lean_assert(l == m);
}
template <typename T, typename X>
void test_swap_cols_with_permutation(sparse_matrix<T, X>& m){
cout << "testing swaps" << endl;
cout << "testing swaps" << std::endl;
unsigned dim = m.row_count();
dense_matrix<double, double> original(m);
permutation_matrix<double, double> q(dim);
@ -596,12 +595,12 @@ void test_swap_cols_with_permutation(sparse_matrix<T, X>& m){
unsigned row1 = lrand48() % dim;
unsigned row2 = lrand48() % dim;
if (row1 == row2) continue;
cout << "swap " << row1 << " " << row2 << endl;
cout << "swap " << row1 << " " << row2 << std::endl;
m.swap_rows(row1, row2);
q.transpose_from_right(row1, row2);
lean_assert(original == q * m);
print_matrix(m);
cout << endl;
cout << std::endl;
}
}
@ -697,7 +696,7 @@ void test_pivot_like_swaps_and_pivot(){
for (auto & t : row) {
cout << t << ",";
}
cout << endl;
cout << std::endl;
lp_settings settings;
m.pivot_row_to_row(pivot_row, alpha, target_row, settings);
m.pivot_row_to_row(pivot_row_0, beta, target_row, settings);
@ -963,7 +962,7 @@ void test_conjugate_eta_matrix() {
l.set_diagonal_element(10);
cout << "l" << endl;
cout << "l" << std::endl;
print_matrix(l);
dense_matrix<double, double> lcopy(l);
@ -972,12 +971,12 @@ void test_conjugate_eta_matrix() {
permutation_matrix<double, double> pr = p.get_inverse();
auto conj = lcopy * pr;
cout << "U*pr" << endl;
cout << "U*pr" << std::endl;
print_matrix(conj);
conj = p * conj;
cout << "conj " << endl;
cout << "conj " << std::endl;
print_matrix(conj);
cout << "l" << endl;
cout << "l" << std::endl;
print_matrix(l);
lean_assert(conj == l);
}
@ -1135,19 +1134,19 @@ void update_settings(argument_parser & args_parser, lp_settings& settings) {
if (get_int_from_args_parser("--percent_for_enter", args_parser, n))
settings.percent_of_entering_to_check = n;
if (get_int_from_args_parser("--partial_pivot", args_parser, n)) {
cout << "setting partial pivot constant to " << n << endl;
cout << "setting partial pivot constant to " << n << std::endl;
settings.c_partial_pivoting = n;
}
if (get_int_from_args_parser("--density", args_parser, n)) {
double density = static_cast<double>(n) / 100.0;
cout << "setting density to " << density << endl;
cout << "setting density to " << density << std::endl;
settings.density_threshold = density;
}
if (get_int_from_args_parser("--maxng", args_parser, n))
settings.max_number_of_iterations_with_no_improvements = n;
double d;
if (get_double_from_args_parser("--harris_toler", args_parser, d)) {
cout << "setting harris_feasibility_tolerance to " << d << endl;
cout << "setting harris_feasibility_tolerance to " << d << std::endl;
settings.harris_feasibility_tolerance = d;
}
}
@ -1172,7 +1171,7 @@ void print_x(mps_reader<double, double> & reader, lp_solver<double, double> * so
for (auto name : reader.column_names()) {
std::cout << name << "=" << solver->get_column_value_by_name(name) << ' ';
}
cout << endl;
cout << std::endl;
}
void compare_solutions(mps_reader<double, double> & reader, lp_solver<double, double> * solver, lp_solver<double, double> * solver0) {
@ -1180,7 +1179,7 @@ void compare_solutions(mps_reader<double, double> & reader, lp_solver<double, do
double a = solver->get_column_value_by_name(name);
double b = solver0->get_column_value_by_name(name);
if (!values_are_one_percent_close(a, b)) {
cout << "different values for " << name << ":" << a << " and " << b << endl;
cout << "different values for " << name << ":" << a << " and " << b << std::endl;
}
}
}
@ -1197,7 +1196,7 @@ void solve_mps_double(std::string file_name, bool look_for_min, unsigned max_ite
setup_solver(max_iterations, time_limit, look_for_min, args_parser, solver);
int begin = get_millisecond_count();
if (dual) {
cout << "solving for dual" << endl;
cout << "solving for dual" << std::endl;
}
solver->find_maximal_solution();
int span = get_millisecond_span(begin);
@ -1212,13 +1211,13 @@ void solve_mps_double(std::string file_name, bool look_for_min, unsigned max_ite
}
std::cout << "cost = " << cost << std::endl;
}
cout << "processed in " << span / 1000.0 << " seconds, running for " << solver->m_total_iterations << " iterations" << endl;
cout << "processed in " << span / 1000.0 << " seconds, running for " << solver->m_total_iterations << " iterations" << std::endl;
if (compare_with_primal) {
auto * primal_solver = reader.create_solver(false);
setup_solver(max_iterations, time_limit, look_for_min, args_parser, primal_solver);
primal_solver->find_maximal_solution();
if (solver->get_status() != primal_solver->get_status()) {
cout << "statuses are different: dual " << lp_status_to_string(solver->get_status()) << " primal = " << lp_status_to_string(primal_solver->get_status()) << endl;
cout << "statuses are different: dual " << lp_status_to_string(solver->get_status()) << " primal = " << lp_status_to_string(primal_solver->get_status()) << std::endl;
} else {
if (solver->get_status() == lp_status::OPTIMAL) {
double cost = solver->get_current_cost();
@ -1229,11 +1228,11 @@ void solve_mps_double(std::string file_name, bool look_for_min, unsigned max_ite
if (look_for_min) {
primal_cost = -primal_cost;
}
cout << "primal cost = " << primal_cost << endl;
cout << "primal cost = " << primal_cost << std::endl;
if (!values_are_one_percent_close(cost, primal_cost)) {
compare_solutions(reader, primal_solver, solver);
print_x(reader, primal_solver);
cout << "dual cost is " << cost << ", but primal cost is " << primal_cost << endl;
cout << "dual cost is " << cost << ", but primal cost is " << primal_cost << std::endl;
lean_assert(false);
}
}
@ -1272,7 +1271,7 @@ void solve_mps_rational(std::string file_name, bool look_for_min, unsigned max_i
}
std::cout << "cost = " << cost.get_double() << std::endl;
}
cout << "processed in " << get_millisecond_span(begin) / 1000.0 << " seconds, running for " << solver->m_total_iterations << " iterations" << endl;
cout << "processed in " << get_millisecond_span(begin) / 1000.0 << " seconds, running for " << solver->m_total_iterations << " iterations" << std::endl;
delete solver;
} else {
std::cout << "cannot process " << file_name << std::endl;
@ -1360,7 +1359,7 @@ void test_binary_priority_queue() {
for (unsigned i = 0; i < 10; i++) {
unsigned de = q.dequeue();
lean_assert(i == de);
cout << de << endl;
cout << de << std::endl;
}
q.enqueue(2, 2);
q.enqueue(1, 1);
@ -1384,11 +1383,11 @@ void test_binary_priority_queue() {
while (q.size() > 0) {
unsigned d =q.dequeue();
lean_assert(t++ == d);
cout << d << endl;
cout << d << std::endl;
}
test_upair_queue();
cout << " done" << endl;
cout << " done" << std::endl;
}
bool solution_is_feasible(std::string file_name, const std::unordered_map<string, double> & solution) {
@ -1411,7 +1410,7 @@ void solve_mps_with_known_solution(std::string file_name, std::unordered_map<str
solver->find_maximal_solution();
std::cout << "status is " << lp_status_to_string(solver->get_status()) << std::endl;
if (status != solver->get_status()){
cout << "status should be " << lp_status_to_string(status) << endl;
cout << "status should be " << lp_status_to_string(status) << std::endl;
lean_assert(status == solver->get_status());
throw "status is wrong";
}
@ -1430,7 +1429,7 @@ void solve_mps_with_known_solution(std::string file_name, std::unordered_map<str
for (auto name : reader.column_names()) {
std::cout << name << "=" << solver->get_column_value_by_name(name) << ' ';
}
cout << endl;
cout << std::endl;
}
}
delete solver;
@ -1484,7 +1483,7 @@ void random_test_on_i(unsigned i) {
srand(i);
auto *solver = generate_random_solver();
solver->find_maximal_solution();
// cout << lp_status_to_string(solver->get_status()) << endl;
// cout << lp_status_to_string(solver->get_status()) << std::endl;
delete solver;
}
@ -1494,7 +1493,7 @@ void random_test() {
random_test_on_i(i);
}
catch (const char * error) {
cout << "i = " << i << ", throwing at ' " << error << "'" << endl;
cout << "i = " << i << ", throwing at ' " << error << "'" << std::endl;
break;
}
}
@ -1634,10 +1633,10 @@ void fill_file_names(std::vector<std::string> &file_names, std::set<string> & m
void test_out_dir(string out_dir) {
DIR *out_dir_p = opendir(out_dir.c_str());
if (out_dir_p == nullptr) {
cout << "creating directory " << out_dir << endl;
cout << "creating directory " << out_dir << std::endl;
int res = mkdir(out_dir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
if (res) {
cout << "Cannot open output directory \"" << out_dir << "\"" << endl;
cout << "Cannot open output directory \"" << out_dir << "\"" << std::endl;
}
return;
}
@ -1648,13 +1647,13 @@ void find_dir_and_file_name(string a, string & dir, string& fn) {
// todo: make it system independent
size_t last_slash_pos = a.find_last_of("/");
if (last_slash_pos >= a.size()) {
cout << "cannot find file name in " << a << endl;
cout << "cannot find file name in " << a << std::endl;
throw;
}
dir = a.substr(0, last_slash_pos);
// cout << "dir = " << dir << endl;
// cout << "dir = " << dir << std::endl;
fn = a.substr(last_slash_pos + 1);
// cout << "fn = " << fn << endl;
// cout << "fn = " << fn << std::endl;
}
void process_test_file(string test_dir, string test_file_name, argument_parser & args_parser, string out_dir, unsigned max_iters, unsigned time_limit, unsigned & successes, unsigned & failures, unsigned & inconclusives);
@ -1689,7 +1688,7 @@ void solve_some_mps(argument_parser & args_parser) {
std::cout<< "exception: "<< s << std::endl;
}
}
cout << "comparing with glpk: successes " << successes << ", failures " << failures << ", inconclusives " << inconclusives << endl;
cout << "comparing with glpk: successes " << successes << ", failures " << failures << ", inconclusives " << inconclusives << std::endl;
return;
}
if (!solve_for_rational) {
@ -1807,7 +1806,7 @@ bool contains(string const & s, char const * pattern) {
unordered_map<string, double> * get_solution_from_glpsol_output(string & file_name) {
ifstream file(file_name);
if (!file.is_open()){
cerr << "cannot open " << file_name << endl;
cerr << "cannot open " << file_name << std::endl;
return nullptr;
}
string s;
@ -1815,7 +1814,7 @@ unordered_map<string, double> * get_solution_from_glpsol_output(string & file_na
do {
s = read_line(end, file);
if (end) {
cerr << "unexpected file end " << file_name << endl;
cerr << "unexpected file end " << file_name << std::endl;
return nullptr;
}
if (contains(s, "Column name")){
@ -1825,7 +1824,7 @@ unordered_map<string, double> * get_solution_from_glpsol_output(string & file_na
read_line(end, file);
if (end) {
cerr << "unexpected file end " << file_name << endl;
cerr << "unexpected file end " << file_name << std::endl;
return nullptr;
}
@ -1834,7 +1833,7 @@ unordered_map<string, double> * get_solution_from_glpsol_output(string & file_na
do {
s = read_line(end, file);
if (end) {
cerr << "unexpected file end " << file_name << endl;
cerr << "unexpected file end " << file_name << std::endl;
return nullptr;
}
auto split = string_split(s, " \t", false);
@ -1956,7 +1955,7 @@ void solve_test_flipped(bool dual) {
// solving a problem with a constraint xj <= c, a flipped constraint
char * home_dir = getenv("HOME");
if (home_dir == nullptr) {
cout << "cannot find home directory" << endl;
cout << "cannot find home directory" << std::endl;
return;
}
string file_name = string(home_dir) + "/projects/lean/src/tests/util/lp/l4.mps";
@ -1968,7 +1967,7 @@ void solve_test_flipped(bool dual) {
solver->find_maximal_solution();
lean_assert(solver->get_status() == OPTIMAL);
double x1_val = solver->get_column_value_by_name("X1");
cout << "X1 = " << x1_val << endl;
cout << "X1 = " << x1_val << std::endl;
mps_reader<double, double> reader_(file_name);
reader_.read();
auto solver_ = reader_.create_solver(dual);
@ -1977,7 +1976,7 @@ void solve_test_flipped(bool dual) {
solver_-> unset_low_bound(j);
solver_->set_upper_bound(j, x1_val + 1);
solver_->find_maximal_solution();
cout << "new X1 = " << solver_->get_column_value_by_name("X1") << endl;
cout << "new X1 = " << solver_->get_column_value_by_name("X1") << std::endl;
lean_assert(fabs(x1_val - solver_->get_column_value_by_name("X1")) < 1e-10);
delete solver;
delete solver_;
@ -1989,7 +1988,7 @@ void print_chunk(T * arr, unsigned len) {
for (unsigned i = 0; i < len; i++) {
cout << arr[i] << ", ";
}
cout << endl;
cout << std::endl;
}
struct mem_cpy_place_holder {
@ -2043,7 +2042,7 @@ int run_glpk(string file_name, string glpk_out_file_name, bool minimize, unsigne
string get_status(string file_name) {
std::ifstream f(file_name);
if (!f.is_open()) {
cout << "cannot open " << file_name << endl;
cout << "cannot open " << file_name << std::endl;
throw 0;
}
string str;
@ -2051,13 +2050,13 @@ string get_status(string file_name) {
if (str.find("Status") != string::npos) {
vector<string> tokens = split_and_trim(str);
if (tokens.size() != 2) {
cout << "unexpected Status string " << str << endl;
cout << "unexpected Status string " << str << std::endl;
throw 0;
}
return tokens[1];
}
}
cout << "cannot find the status line in " << file_name << endl;
cout << "cannot find the status line in " << file_name << std::endl;
throw 0;
}
@ -2072,7 +2071,7 @@ bool compare_statuses(string glpk_out_file_name, string lp_out_file_name, unsign
return false;
} else {
cout << "glpsol and lp_tst disagree: glpsol status is " << glpk_status;
cout << " but lp_tst status is " << lp_tst_status << endl;
cout << " but lp_tst status is " << lp_tst_status << std::endl;
failures++;
return false;
}
@ -2083,7 +2082,7 @@ bool compare_statuses(string glpk_out_file_name, string lp_out_file_name, unsign
double get_glpk_cost(string file_name) {
std::ifstream f(file_name);
if (!f.is_open()) {
cout << "cannot open " << file_name << endl;
cout << "cannot open " << file_name << std::endl;
throw 0;
}
string str;
@ -2091,20 +2090,20 @@ double get_glpk_cost(string file_name) {
if (str.find("Objective") != string::npos) {
vector<string> tokens = split_and_trim(str);
if (tokens.size() != 5) {
cout << "unexpected Objective string " << str << endl;
cout << "unexpected Objective string " << str << std::endl;
throw 0;
}
return atof(tokens[3].c_str());
}
}
cout << "cannot find the Objective line in " << file_name << endl;
cout << "cannot find the Objective line in " << file_name << std::endl;
throw 0;
}
double get_lp_tst_cost(string file_name) {
std::ifstream f(file_name);
if (!f.is_open()) {
cout << "cannot open " << file_name << endl;
cout << "cannot open " << file_name << std::endl;
throw 0;
}
string str;
@ -2115,13 +2114,13 @@ double get_lp_tst_cost(string file_name) {
}
}
if (cost_string.size() == 0) {
cout << "cannot find the cost line in " << file_name << endl;
cout << "cannot find the cost line in " << file_name << std::endl;
throw 0;
}
vector<string> tokens = split_and_trim(cost_string);
if (tokens.size() != 3) {
cout << "unexpected cost string " << cost_string << endl;
cout << "unexpected cost string " << cost_string << std::endl;
throw 0;
}
return atof(tokens[2].c_str());
@ -2149,7 +2148,7 @@ void compare_costs(string glpk_out_file_name,
successes++;
} else {
failures++;
cout << "glpsol cost is " << a << " lp_tst cost is " << b << endl;
cout << "glpsol cost is " << a << " lp_tst cost is " << b << std::endl;
}
}
@ -2158,9 +2157,9 @@ void compare_costs(string glpk_out_file_name,
void compare_with_glpk(string glpk_out_file_name, string lp_out_file_name, unsigned & successes, unsigned & failures, string lp_file_name) {
std::unordered_map<string, double> * solution_table = get_solution_from_glpsol_output(glpk_out_file_name);
if (solution_is_feasible(lp_file_name, *solution_table)) {
cout << "glpk solution is feasible" << endl;
cout << "glpk solution is feasible" << std::endl;
} else {
cout << "glpk solution is infeasible" << endl;
cout << "glpk solution is infeasible" << std::endl;
}
delete solution_table;
if (compare_statuses(glpk_out_file_name, lp_out_file_name, successes, failures)) {
@ -2176,14 +2175,14 @@ void process_test_file(string test_dir, string test_file_name, argument_parser &
string input_file_name = test_dir + "/" + test_file_name;
if (input_file_name[input_file_name.size() - 1] == '~') {
// cout << "ignoring " << input_file_name << endl;
// cout << "ignoring " << input_file_name << std::endl;
return;
}
cout <<"processing " << input_file_name << endl;
cout <<"processing " << input_file_name << std::endl;
std::ofstream out(full_lp_tst_out_name);
if (!out.is_open()) {
cout << "cannot open file " << full_lp_tst_out_name << endl;
cout << "cannot open file " << full_lp_tst_out_name << std::endl;
throw 0;
}
std::streambuf *coutbuf = std::cout.rdbuf(); // save old buffer
@ -2196,7 +2195,7 @@ void process_test_file(string test_dir, string test_file_name, argument_parser &
solve_mps(input_file_name, minimize, max_iters, time_limit, use_mpq, dual, false, args_parser);
}
catch(...) {
cout << "catching the failure" << endl;
cout << "catching the failure" << std::endl;
failures++;
std::cout.rdbuf(coutbuf); // reset to standard output again
return;
@ -2207,7 +2206,7 @@ void process_test_file(string test_dir, string test_file_name, argument_parser &
string glpk_out_file_name = out_dir + "/" + create_output_file_name_for_glpsol(minimize, string(test_file_name));
int glpk_exit_code = run_glpk(input_file_name, glpk_out_file_name, minimize, time_limit);
if (glpk_exit_code != 0) {
cout << "glpk failed" << endl;
cout << "glpk failed" << std::endl;
inconclusives++;
} else {
compare_with_glpk(glpk_out_file_name, full_lp_tst_out_name, successes, failures, input_file_name);
@ -2215,13 +2214,13 @@ void process_test_file(string test_dir, string test_file_name, argument_parser &
}
}
vector<pair<string, int>> get_file_list_of_dir(string test_file_dir) {
std::vector<std::pair<std::string, int>> get_file_list_of_dir(std::string test_file_dir) {
DIR *dir;
if ((dir = opendir(test_file_dir.c_str())) == nullptr) {
cout << "Cannot open directory " << test_file_dir << endl;
std::cout << "Cannot open directory " << test_file_dir << std::endl;
throw 0;
}
vector<pair<string, int>> ret;
std::vector<std::pair<std::string, int>> ret;
struct dirent entry;
struct dirent* result;
int return_code;
@ -2231,7 +2230,7 @@ vector<pair<string, int>> get_file_list_of_dir(string test_file_dir) {
DIR *tmp_dp = opendir(entry.d_name);
struct stat file_record;
if (tmp_dp == nullptr) {
string s = test_file_dir+ "/" + entry.d_name;
std::string s = test_file_dir+ "/" + entry.d_name;
int stat_ret = stat(s.c_str(), & file_record);
if (stat_ret!= -1) {
ret.push_back(make_pair(entry.d_name, file_record.st_size));
@ -2249,24 +2248,24 @@ vector<pair<string, int>> get_file_list_of_dir(string test_file_dir) {
struct file_size_comp {
unordered_map<string, int>& m_file_sizes;
file_size_comp(unordered_map<string, int>& fs) :m_file_sizes(fs) {}
int operator()(string a, string b) {
cout << m_file_sizes.size() << endl;
cout << a << endl;
cout << b << endl;
unordered_map<std::string, int>& m_file_sizes;
file_size_comp(unordered_map<std::string, int>& fs) :m_file_sizes(fs) {}
int operator()(std::string a, std::string b) {
std::cout << m_file_sizes.size() << std::endl;
std::cout << a << std::endl;
std::cout << b << std::endl;
auto ls = m_file_sizes.find(a);
cout << "fa" << endl;
std::cout << "fa" << std::endl;
auto rs = m_file_sizes.find(b);
cout << "fb" << endl;
std::cout << "fb" << std::endl;
if (ls != m_file_sizes.end() && rs != m_file_sizes.end()) {
cout << "fc " << endl;
std::cout << "fc " << std::endl;
int r = (*ls < *rs? -1: (*ls > *rs)? 1 : 0);
cout << "calc r " << endl;
std::cout << "calc r " << std::endl;
return r;
} else {
cout << "sc " << endl;
std::cout << "sc " << std::endl;
return 0;
}
}
@ -2274,25 +2273,25 @@ struct file_size_comp {
struct sort_pred {
bool operator()(const std::pair<string, int> &left, const std::pair<string, int> &right) {
bool operator()(const std::pair<std::string, int> &left, const std::pair<std::string, int> &right) {
return left.second < right.second;
}
};
void test_files_from_directory(string test_file_dir, argument_parser & args_parser) {
cout << "loading files from directory \"" << test_file_dir << "\"" << endl;
string out_dir = args_parser.get_option_value("--out_dir");
void test_files_from_directory(std::string test_file_dir, argument_parser & args_parser) {
std::cout << "loading files from directory \"" << test_file_dir << "\"" << std::endl;
std::string out_dir = args_parser.get_option_value("--out_dir");
if (out_dir.size() == 0) {
out_dir = "/tmp/test";
}
DIR *out_dir_p = opendir(out_dir.c_str());
if (out_dir_p == nullptr) {
cout << "Cannot open output directory \"" << out_dir << "\"" << endl;
std::cout << "Cannot open output directory \"" << out_dir << "\"" << std::endl;
return;
}
closedir(out_dir_p);
vector<pair<string, int>> files = get_file_list_of_dir(test_file_dir);
std::vector<std::pair<std::string, int>> files = get_file_list_of_dir(test_file_dir);
std::sort(files.begin(), files.end(), sort_pred());
unsigned max_iters, time_limit;
get_time_limit_and_max_iters_from_parser(args_parser, time_limit, max_iters);
@ -2300,12 +2299,12 @@ void test_files_from_directory(string test_file_dir, argument_parser & args_pars
for (auto & t : files) {
process_test_file(test_file_dir, t.first, args_parser, out_dir, max_iters, time_limit, successes, failures, inconclusives);
}
cout << "comparing with glpk: successes " << successes << ", failures " << failures << ", inconclusives " << inconclusives << endl;
std::cout << "comparing with glpk: successes " << successes << ", failures " << failures << ", inconclusives " << inconclusives << std::endl;
}
unordered_map<string, mpq> get_solution_map(lp_solver<mpq, mpq> * lps, mps_reader<mpq, mpq> & reader) {
unordered_map<string, mpq> ret;
unordered_map<std::string, mpq> get_solution_map(lp_solver<mpq, mpq> * lps, mps_reader<mpq, mpq> & reader) {
unordered_map<std::string, mpq> ret;
for (auto it : reader.column_names()) {
ret[it] = lps->get_column_value_by_name(it);
}
@ -2313,7 +2312,7 @@ unordered_map<string, mpq> get_solution_map(lp_solver<mpq, mpq> * lps, mps_reade
}
void run_lar_solver(argument_parser & args_parser, lar_solver * solver, mps_reader<mpq, mpq> * reader) {
string maxng = args_parser.get_option_value("--maxng");
std::string maxng = args_parser.get_option_value("--maxng");
if (maxng.size() > 0) {
solver->settings().max_number_of_iterations_with_no_improvements = atoi(maxng.c_str());
}
@ -2323,34 +2322,34 @@ void run_lar_solver(argument_parser & args_parser, lar_solver * solver, mps_read
if (args_parser.option_is_used("--mpq")) {
solver->settings().use_double_solver_for_lar = false;
}
string iter = args_parser.get_option_value("--max_iters");
std::string iter = args_parser.get_option_value("--max_iters");
if (iter.size() > 0) {
solver->settings().max_total_number_of_iterations = atoi(iter.c_str());
}
if (args_parser.option_is_used("--compare_with_primal")){
if (reader == nullptr) {
cout << "cannot compare with primal, the reader is null " << endl;
std::cout << "cannot compare with primal, the reader is null " << std::endl;
return;
}
auto * lps = reader->create_solver(false);
lps->find_maximal_solution();
unordered_map<string, mpq> sol = get_solution_map(lps, *reader);
unordered_map<std::string, mpq> sol = get_solution_map(lps, *reader);
mpq inf = solver->get_infeasibility_of_solution(sol);
cout << "inf with primal = " << inf << endl;
std::cout << "inf with primal = " << inf << std::endl;
return;
}
int begin = get_millisecond_count();
lp_status status = solver->check();
cout << "status is " << lp_status_to_string(status) << ", processed for " << get_millisecond_span(begin) / 1000.0 <<" seconds, and " << solver->get_total_iterations() << " iterations" << endl;
std::cout << "status is " << lp_status_to_string(status) << ", processed for " << get_millisecond_span(begin) / 1000.0 <<" seconds, and " << solver->get_total_iterations() << " iterations" << std::endl;
if (solver->get_status() == INFEASIBLE) {
buffer<pair<mpq, constraint_index>> evidence;
buffer<std::pair<mpq, constraint_index>> evidence;
solver->get_infeasibility_evidence(evidence);
}
}
void test_lar_on_file(string file_name, argument_parser & args_parser) {
void test_lar_on_file(std::string file_name, argument_parser & args_parser) {
lar_solver * solver = nullptr;
cout << "processing " << file_name << endl;
std::cout << "processing " << file_name << std::endl;
if (args_parser.option_is_used("--smt")) {
smt_reader reader(file_name);
reader.read();
@ -2374,16 +2373,16 @@ void test_lar_on_file(string file_name, argument_parser & args_parser) {
delete solver;
}
vector<string> get_file_names_from_file_list(string filelist) {
vector<std::string> get_file_names_from_file_list(std::string filelist) {
ifstream file(filelist);
if (!file.is_open()) {
cout << "cannot open " << filelist << endl;
return vector<string>();
std::cout << "cannot open " << filelist << std::endl;
return vector<std::string>();
}
vector<string> ret;
vector<std::string> ret;
bool end;
do {
string s = read_line(end, file);
std::string s = read_line(end, file);
if (end)
break;
if (s.size() == 0)
@ -2394,15 +2393,15 @@ vector<string> get_file_names_from_file_list(string filelist) {
}
void test_lar_solver(argument_parser & args_parser) {
string file_name = args_parser.get_option_value("--file");
std::string file_name = args_parser.get_option_value("--file");
if (file_name.size() > 0) {
test_lar_on_file(file_name, args_parser);
return;
}
string file_list = args_parser.get_option_value("--filelist");
std::string file_list = args_parser.get_option_value("--filelist");
if (file_list.size() > 0) {
for (string fn : get_file_names_from_file_list(file_list))
for (std::string fn : get_file_names_from_file_list(file_list))
test_lar_on_file(fn, args_parser);
return;
}
@ -2417,7 +2416,7 @@ void test_numeric_pair() {
a -= c;
lean_assert (a == b + c);
numeric_pair<mpq> d = a * 2;
cout << a << endl;
std::cout << a << std::endl;
lean_assert(b == b);
lean_assert(b < a);
lean_assert(b <= a);
@ -2429,15 +2428,15 @@ void test_numeric_pair() {
lean_assert(b + b > a);
lean_assert(mpq(2.1) * b + b > a);
lean_assert(-b * mpq(2.1) - b < mpq(0.99) * a);
cout << - b * mpq(2.1) - b << endl;
std::cout << - b * mpq(2.1) - b << std::endl;
lean_assert(-b *(mpq(2.1) + 1) == - b * mpq(2.1) - b);
}
void get_matrix_dimensions(ifstream & f, unsigned & m, unsigned & n) {
string line;
std::string line;
getline(f, line);
getline(f, line);
vector<string> r = split_and_trim(line);
vector<std::string> r = split_and_trim(line);
m = atoi(r[1].c_str());
getline(f, line);
r = split_and_trim(line);
@ -2446,7 +2445,7 @@ void get_matrix_dimensions(ifstream & f, unsigned & m, unsigned & n) {
void read_row_cols(unsigned i, static_matrix<double, double>& A, ifstream & f) {
do {
string line;
std::string line;
getline(f, line);
if (line== "row_end")
break;
@ -2459,13 +2458,13 @@ void read_row_cols(unsigned i, static_matrix<double, double>& A, ifstream & f) {
}
bool read_row(static_matrix<double, double> & A, ifstream & f) {
string line;
std::string line;
getline(f, line);
if (static_cast<int>(line.find("row")) == -1)
return false;
auto r = split_and_trim(line);
if (r[0] != "row")
cout << "wrong row line" << line << endl;
std::cout << "wrong row line" << line << std::endl;
unsigned i = atoi(r[1].c_str());
read_row_cols(i, A, f);
return true;
@ -2476,8 +2475,8 @@ void read_rows(static_matrix<double, double>& A, ifstream & f) {
}
void read_basis(vector<unsigned> & basis, ifstream & f) {
cout << "reading basis" << endl;
string line;
std::cout << "reading basis" << std::endl;
std::string line;
getline(f, line);
lean_assert(line == "basis_start");
do {
@ -2490,7 +2489,7 @@ void read_basis(vector<unsigned> & basis, ifstream & f) {
}
void read_indexed_vector(indexed_vector<double> & v, ifstream & f) {
string line;
std::string line;
getline(f, line);
lean_assert(line == "vector_start");
do {
@ -2500,18 +2499,18 @@ void read_indexed_vector(indexed_vector<double> & v, ifstream & f) {
unsigned i = atoi(r[0].c_str());
double val = atof(r[1].c_str());
v.set_value(val, i);
cout << "setting value " << i << " = " << val << endl;
std::cout << "setting value " << i << " = " << val << std::endl;
} while (true);
}
void check_lu_from_file(string lufile_name) {
void check_lu_from_file(std::string lufile_name) {
ifstream f(lufile_name);
if (!f.is_open()) {
cout << "cannot open file " << lufile_name << endl;
std::cout << "cannot open file " << lufile_name << std::endl;
}
unsigned m, n;
get_matrix_dimensions(f, m, n);
cout << "init matrix " << m << " by " << n << endl;
std::cout << "init matrix " << m << " by " << n << std::endl;
static_matrix<double, double> A(m, n);
read_rows(A, f);
vector<unsigned> basis;
@ -2540,7 +2539,7 @@ void check_lu_from_file(string lufile_name) {
}
void test_square_dense_submatrix() {
cout << "testing square_dense_submatrix" << endl;
std::cout << "testing square_dense_submatrix" << std::endl;
unsigned parent_dim = 7;
sparse_matrix<double, double> parent(parent_dim);
fill_matrix(parent);
@ -2567,7 +2566,7 @@ void test_square_dense_submatrix() {
m[i-index_start][j-index_start] = d[i][j];
print_matrix(m);
cout << endl;
std::cout << std::endl;
#endif
}
@ -2580,14 +2579,14 @@ int main(int argn, char * const * argv) {
setup_args_parser(args_parser);
if (!args_parser.parse()) {
cout << args_parser.m_error_message << endl;
cout << args_parser.usage_string();
std::cout << args_parser.m_error_message << std::endl;
std::cout << args_parser.usage_string();
ret = 1;
return finalize(ret);
}
cout << "the options are " << endl;
std::cout << "the options are " << std::endl;
args_parser.print();
string lufile = args_parser.get_option_value("--checklu");
std::string lufile = args_parser.get_option_value("--checklu");
if (lufile.size()) {
check_lu_from_file(lufile);
return finalize(0);
@ -2607,13 +2606,13 @@ int main(int argn, char * const * argv) {
return finalize(0);
}
if (args_parser.option_is_used("--lar")){
cout <<"calling test_lar_solver" << endl;
std::cout <<"calling test_lar_solver" << std::endl;
test_lar_solver(args_parser);
return finalize(0);
}
string file_list = args_parser.get_option_value("--filelist");
std::string file_list = args_parser.get_option_value("--filelist");
if (file_list.size() > 0) {
for (string fn : get_file_names_from_file_list(file_list))
for (std::string fn : get_file_names_from_file_list(file_list))
solve_mps(fn, args_parser);
return finalize(0);
}
@ -2658,7 +2657,7 @@ int main(int argn, char * const * argv) {
get_time_limit_and_max_iters_from_parser(args_parser, time_limit, max_iters);
bool dual = args_parser.option_is_used("--dual");
bool solve_for_rational = args_parser.option_is_used("--mpq");
string file_name = args_parser.get_option_value("--file");
std::string file_name = args_parser.get_option_value("--file");
if (file_name.size() > 0) {
solve_mps(file_name, args_parser.option_is_used("--min"), max_iters, time_limit, solve_for_rational, dual, args_parser.option_is_used("--compare_with_primal"), args_parser);
ret = 0;

76
src/util/lp/binary_heap_priority_queue.h
View file

@ -4,19 +4,17 @@
Author: Lev Nachmanson
*/
#pragma once
#include <vector>
namespace lean {
using std::vector;
// the elements with the smallest priority are dequeued first
// the elements with the smallest priority are dequeued first
template <typename T>
class binary_heap_priority_queue {
vector<T> m_priorities;
std::vector<T> m_priorities;
// indexing for A starts from 1
vector<unsigned> m_heap; // keeps the elements of the queue
vector<int> m_heap_inverse; // o = m_heap[m_heap_inverse[o]]
std::vector<unsigned> m_heap; // keeps the elements of the queue
std::vector<int> m_heap_inverse; // o = m_heap[m_heap_inverse[o]]
unsigned m_heap_size = 0;
// is is the child place in heap
@ -33,15 +31,15 @@ class binary_heap_priority_queue {
}
void decrease_priority(unsigned o, T newPriority) {
m_priorities[o] = newPriority;
int i = m_heap_inverse[o];
while (i > 1) {
if (m_priorities[m_heap[i]] < m_priorities[m_heap[i >> 1]])
swap_with_parent(i);
else
break;
i >>= 1;
}
m_priorities[o] = newPriority;
int i = m_heap_inverse[o];
while (i > 1) {
if (m_priorities[m_heap[i]] < m_priorities[m_heap[i >> 1]])
swap_with_parent(i);
else
break;
i >>= 1;
}
}
public:
@ -56,13 +54,13 @@ public:
for (int k = 0; k < 2; k++) {
if (ch > m_heap_size) break;
if (!(m_priorities[m_heap[i]] <= m_priorities[m_heap[ch]])){
cout << "m_heap_size = " << m_heap_size << endl;
cout << "i = " << i << endl;
cout << "m_heap[i] = " << m_heap[i] << endl;
cout << "ch = " << ch << endl;
cout << "m_heap[ch] = " << m_heap[ch] << endl;
cout << "m_priorities[m_heap[i]] = " << m_priorities[m_heap[i]] << endl;
cout << "m_priorities[m_heap[ch]] = " << m_priorities[m_heap[ch]]<< endl;
std::cout << "m_heap_size = " << m_heap_size << std::endl;
std::cout << "i = " << i << std::endl;
std::cout << "m_heap[i] = " << m_heap[i] << std::endl;
std::cout << "ch = " << ch << std::endl;
std::cout << "m_heap[ch] = " << m_heap[ch] << std::endl;
std::cout << "m_priorities[m_heap[i]] = " << m_priorities[m_heap[i]] << std::endl;
std::cout << "m_priorities[m_heap[ch]] = " << m_priorities[m_heap[ch]]<< std::endl;
return false;
}
ch++;
@ -109,7 +107,7 @@ public:
m_priorities(n),
m_heap(n + 1), // because the indexing for A starts from 1
m_heap_inverse(n, -1)
{ }
{ }
void clear() {
m_heap_size = 0;
@ -127,15 +125,15 @@ public:
}
void enqueue_new(unsigned o, const T& priority) {
m_heap_size++;
int i = m_heap_size;
lean_assert(o < m_priorities.size());
m_priorities[o] = priority;
put_at(i, o);
while (i > 1 && m_priorities[m_heap[i >> 1]] > priority) {
swap_with_parent(i);
i >>= 1;
}
m_heap_size++;
int i = m_heap_size;
lean_assert(o < m_priorities.size());
m_priorities[o] = priority;
put_at(i, o);
while (i > 1 && m_priorities[m_heap[i >> 1]] > priority) {
swap_with_parent(i);
i >>= 1;
}
}
// This method can work with an element that is already in the queue.
// In this case the priority will be changed and the queue adjusted.
@ -168,9 +166,7 @@ public:
/// return the first element of the queue and removes it from the queue
unsigned dequeue_and_get_priority(T & priority) {
if (m_heap_size == 0) {
throw "invalid size";
}
lean_assert(m_heap_size != 0);
int ret = m_heap[1];
priority = m_priorities[ret];
put_the_last_at_the_top_and_fix_the_heap();
@ -202,7 +198,7 @@ public:
m_heap_size--;
}
}
/// return the first element of the queue and removes it from the queue
/// return the first element of the queue and removes it from the queue
unsigned dequeue() {
lean_assert(m_heap_size);
int ret = m_heap[1];
@ -212,16 +208,16 @@ public:
}
void print() {
vector<int> index;
vector<T> prs;
std::vector<int> index;
std::vector<T> prs;
while (size()) {
T prior;
int j = dequeue_and_get_priority(prior);
index.push_back(j);
prs.push_back(prior);
cout << "(" << j << ", " << prior << ")";
std::cout << "(" << j << ", " << prior << ")";
}
cout << endl;
std::cout << std::endl;
// restore the queue
for (int i = 0; i < index.size(); i++)
enqueue(index[i], prs[i]);

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

@ -18,8 +18,8 @@ namespace lean {
template <typename T>
class binary_heap_upair_queue {
binary_heap_priority_queue<T> m_q;
unordered_map<upair, unsigned> m_pairs_to_index;
vector<upair> m_pairs; // inverse to index
std::unordered_map<upair, unsigned> m_pairs_to_index;
std::vector<upair> m_pairs; // inverse to index
std::vector<unsigned> m_available_spots;
public:
binary_heap_upair_queue(unsigned size) : m_q(size), m_pairs(size) {
@ -113,7 +113,8 @@ public:
unsigned j = p.second;
auto it = tmp.find(j);
if (it != tmp.end()) {
cout << "for pair (" << p.first.first << ", " << p.first.second << "), the index " << j << " is already inside " << endl;
std::cout << "for pair (" << p.first.first << ", " << p.first.second << "), the index " << j
<< " is already inside " << std::endl;
lean_assert(false);
} else {
tmp.insert(j);

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

@ -11,74 +11,74 @@
#include <algorithm>
namespace lean {
typedef unsigned var_index;
typedef unsigned constraint_index;
enum lconstraint_kind {
LE = -2, LT = -1 , GE = 2, GT = 1, EQ = 0
};
typedef unsigned var_index;
typedef unsigned constraint_index;
enum lconstraint_kind {
LE = -2, LT = -1 , GE = 2, GT = 1, EQ = 0
};
class lar_normalized_constraint; // forward definition
bool compare(const pair<mpq, var_index> & a, const pair<mpq, var_index> & b) {
return a.second < b.second;
class lar_normalized_constraint; // forward definition
bool compare(const pair<mpq, var_index> & a, const pair<mpq, var_index> & b) {
return a.second < b.second;
}
class canonic_left_side {
public:
int m_row_index = -1;
int m_column_index = -1; // this is the column of the left side variable in the matrix
std::vector<pair<mpq, var_index>> m_coeffs;
column_info<mpq> m_column_info;
lar_normalized_constraint * m_low_bound_witness = nullptr;
lar_normalized_constraint * m_upper_bound_witness = nullptr;
canonic_left_side(buffer<pair<mpq, var_index>> buffer) {
for (auto it : buffer) {
if (numeric_traits<mpq>::is_zero(it.first)) continue;
m_coeffs.push_back(it);
}
std::sort(m_coeffs.begin(), m_coeffs.end(), compare);
normalize();
}
class canonic_left_side {
public:
int m_row_index = -1;
int m_column_index = -1; // this is the column of the left side variable in the matrix
std::vector<pair<mpq, var_index>> m_coeffs;
column_info<mpq> m_column_info;
lar_normalized_constraint * m_low_bound_witness = nullptr;
lar_normalized_constraint * m_upper_bound_witness = nullptr;
void set_name(std::string name) {
m_column_info.set_name(name);
}
canonic_left_side(buffer<pair<mpq, var_index>> buffer) {
for (auto it : buffer) {
if (numeric_traits<mpq>::is_zero(it.first)) continue;
m_coeffs.push_back(it);
}
unsigned size() const { return m_coeffs.size(); }
std::sort(m_coeffs.begin(), m_coeffs.end(), compare);
normalize();
void normalize() {
if (m_coeffs.size() == 0) return;
auto t = m_coeffs[0].first;
for (auto & it : m_coeffs)
it.first /= t;
}
bool operator==(const canonic_left_side& a) const {
if (m_coeffs.size() != a.m_coeffs.size()) return false;
for (unsigned i = 0; i < m_coeffs.size(); i++) {
if (m_coeffs[i] != a.m_coeffs[i])
return false;
}
return true;
}
void set_name(string name) {
m_column_info.set_name(name);
std::size_t hash_of_ls() const {
std::size_t ret = 0;
std::hash<pair<mpq, var_index>> hash_fun;
for (auto v : m_coeffs) {
ret |= (hash_fun(v) << 2);
}
return ret;
}
};
unsigned size() const { return m_coeffs.size(); }
void normalize() {
if (m_coeffs.size() == 0) return;
auto t = m_coeffs[0].first;
for (auto & it : m_coeffs)
it.first /= t;
}
bool operator==(const canonic_left_side& a) const {
if (m_coeffs.size() != a.m_coeffs.size()) return false;
for (unsigned i = 0; i < m_coeffs.size(); i++) {
if (m_coeffs[i] != a.m_coeffs[i])
return false;
}
return true;
}
std::size_t hash_of_ls() const {
std::size_t ret = 0;
std::hash<pair<mpq, var_index>> hash_fun;
for (auto v : m_coeffs) {
ret |= (hash_fun(v) << 2);
}
return ret;
}
};
struct hash_and_equal_of_canonic_left_side_struct {
std::size_t operator() (const canonic_left_side* ls) const {
return ls->hash_of_ls();
}
bool operator() (const canonic_left_side* a, const canonic_left_side* b) const {
return (*a) == (*b);
}
};
struct hash_and_equal_of_canonic_left_side_struct {
std::size_t operator() (const canonic_left_side* ls) const {
return ls->hash_of_ls();
}
bool operator() (const canonic_left_side* a, const canonic_left_side* b) const {
return (*a) == (*b);
}
};
}

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

@ -13,12 +13,9 @@
#include <string>
#include <algorithm>
namespace lean {
using std::vector;
using std::tuple;
template <typename T>
class column_info {
string m_name;
std::string m_name;
bool m_low_bound_is_set = false;
bool m_low_bound_is_strict = false;
bool m_upper_bound_is_set = false;
@ -124,11 +121,11 @@ public:
m_cost = cost;
}
void set_name(string const & s) {
void set_name(std::string const & s) {
m_name = s;
}
string get_name() const {
std::string get_name() const {
return m_name;
}

52
src/util/lp/core_solver_pretty_printer.h
View file

@ -13,6 +13,8 @@ template <typename T, typename X> class lp_core_solver_base; // forward definiti
template <typename T, typename X>
class core_solver_pretty_printer {
template<typename A> using vector = std::vector<A>;
typedef std::string string;
lp_core_solver_base<T, X> & m_core_solver;
vector<unsigned> m_column_widths;
vector<vector<string>> m_A;
@ -207,7 +209,7 @@ public:
}
int blanks = m_title_width + 1 - m_x_title.size();
cout << m_x_title;
std::cout << m_x_title;
print_blanks(blanks);
auto bh = m_core_solver.m_x;
@ -215,9 +217,9 @@ public:
string s = T_to_string(bh[i]);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
cout << s << " "; // the column interval
std::cout << s << " "; // the column interval
}
cout << endl;
std::cout << std::endl;
}
std::string get_low_bound_string(unsigned j) {
@ -228,7 +230,7 @@ public:
if (m_core_solver.low_bounds_are_set())
return T_to_string(m_core_solver.low_bound_value(j));
else
return string("0");
return std::string("0");
break;
default:
return std::string();
@ -253,16 +255,16 @@ public:
return;
}
int blanks = m_title_width + 1 - m_low_bounds_title.size();
cout << m_low_bounds_title;
std::cout << m_low_bounds_title;
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_low_bound_string(i);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
cout << s << " "; // the column interval
std::cout << s << " "; // the column interval
}
cout << endl;
std::cout << std::endl;
}
void print_upps() {
@ -270,16 +272,16 @@ public:
return;
}
int blanks = m_title_width + 1 - m_upp_bounds_title.size();
cout << m_upp_bounds_title;
std::cout << m_upp_bounds_title;
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_upp_bound_string(i);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
cout << s << " "; // the column interval
std::cout << s << " "; // the column interval
}
cout << endl;
std::cout << std::endl;
}
string get_exact_column_norm_string(unsigned col) {
@ -288,28 +290,28 @@ public:
void print_exact_norms() {
int blanks = m_title_width + 1 - m_exact_norm_title.size();
cout << m_exact_norm_title;
std::cout << m_exact_norm_title;
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_exact_column_norm_string(i);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
cout << s << " ";
std::cout << s << " ";
}
cout << endl;
std::cout << std::endl;
}
void print_approx_norms() {
int blanks = m_title_width + 1 - m_approx_norm_title.size();
cout << m_approx_norm_title;
std::cout << m_approx_norm_title;
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(m_core_solver.m_column_norms[i]);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
cout << s << " ";
std::cout << s << " ";
}
cout << endl;
std::cout << std::endl;
}
void print() {
@ -324,12 +326,12 @@ public:
print_upps();
print_exact_norms();
print_approx_norms();
cout << endl;
std::cout << std::endl;
}
void print_basis_heading() {
int blanks = m_title_width + 1 - m_basis_heading_title.size();
cout << m_basis_heading_title;
std::cout << m_basis_heading_title;
print_blanks(blanks);
if (ncols() == 0) {
@ -340,9 +342,9 @@ public:
string s = T_to_string(bh[i]);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
cout << s << " "; // the column interval
std::cout << s << " "; // the column interval
}
cout << endl;
std::cout << std::endl;
}
void print_bottom_line() {
@ -351,7 +353,7 @@ public:
void print_cost() {
int blanks = m_title_width + 1 - m_cost_title.size();
cout << m_cost_title;
std::cout << m_cost_title;
print_blanks(blanks);
print_given_rows(m_costs, m_cost_signs, m_core_solver.get_cost());
}
@ -363,18 +365,18 @@ public:
int number_of_blanks = width - s.size();
lean_assert(number_of_blanks >= 0);
print_blanks(number_of_blanks);
cout << s << ' ';
std::cout << s << ' ';
if (col < row.size() - 1) {
cout << signs[col + 1] << ' ';
std::cout << signs[col + 1] << ' ';
}
}
cout << '=';
std::cout << '=';
string rs = T_to_string(rst);
int nb = m_rs_width - rs.size();
lean_assert(nb >= 0);
print_blanks(nb + 1);
cout << rs << std::endl;
std::cout << rs << std::endl;
}
void print_row(unsigned i){

259
src/util/lp/eta_matrix.h
View file

@ -8,172 +8,147 @@
#pragma once
namespace lean {
// This is the sum of a unit matrix and a one-column matrix
template <typename T, typename X>
class eta_matrix
: public tail_matrix<T, X> {
// This is the sum of a unit matrix and a one-column matrix
template <typename T, typename X>
class eta_matrix
: public tail_matrix<T, X> {
#ifdef LEAN_DEBUG
unsigned m_length;
unsigned m_length;
#endif
unsigned m_column_index;
sparse_vector<T> m_column_vector;
T m_diagonal_element;
public:
unsigned m_column_index;
sparse_vector<T> m_column_vector;
T m_diagonal_element;
public:
#ifdef LEAN_DEBUG
eta_matrix(unsigned column_index, unsigned length):
eta_matrix(unsigned column_index, unsigned length):
#else
eta_matrix(unsigned column_index):
eta_matrix(unsigned column_index):
#endif
#ifdef LEAN_DEBUG
m_length(length),
m_length(length),
#endif
m_column_index(column_index) {
m_column_index(column_index) {
}
void print() {
print_matrix(*this);
}
bool is_unit() {
return m_column_vector.size() == 0 && m_diagonal_element == 1;
}
void set_diagonal_element(T const & diagonal_element) {
lean_assert(!lp_settings::is_eps_small_general(diagonal_element, 1e-12));
m_diagonal_element = diagonal_element;
}
const T & get_diagonal_element() const {
return m_diagonal_element;
}
void apply_from_left(std::vector<X> & w, lp_settings & ) {
auto w_at_column_index = w[m_column_index];
w[m_column_index] /= m_diagonal_element;
for (auto it = sparse_vector_iterator<T>(m_column_vector); !it.done(); it.move()) {
w[it.index()] += w_at_column_index * it.value();
}
}
template <typename L>
void apply_from_left_local(indexed_vector<L> & w, lp_settings & settings) {
const L w_at_column_index = w[m_column_index];
if (is_zero(w_at_column_index)) return;
if (settings.abs_val_is_smaller_than_drop_tolerance(w[m_column_index] /= m_diagonal_element)) {
w[m_column_index] = zero_of_type<L>();
w.erase_from_index(m_column_index);
}
// T & get_column_element(unsigned j) {
// return m_column[j];
// }
void print() {
print_matrix(*this);
}
bool is_unit() {
return m_column_vector.size() == 0 && m_diagonal_element == 1;
}
void set_diagonal_element(T const & diagonal_element) {
lean_assert(!lp_settings::is_eps_small_general(diagonal_element, 1e-12));
m_diagonal_element = diagonal_element;
}
const T & get_diagonal_element() const {
return m_diagonal_element;
}
void apply_from_left(vector<X> & w, lp_settings & ) {
auto w_at_column_index = w[m_column_index];
w[m_column_index] /= m_diagonal_element;
for (auto it = sparse_vector_iterator<T>(m_column_vector); !it.done(); it.move()) {
w[it.index()] += w_at_column_index * it.value();
}
}
template <typename L>
void apply_from_left_local(indexed_vector<L> & w, lp_settings & settings) {
const L w_at_column_index = w[m_column_index];
if (is_zero(w_at_column_index)) return;
if (settings.abs_val_is_smaller_than_drop_tolerance(w[m_column_index] /= m_diagonal_element)) {
w[m_column_index] = zero_of_type<L>();
w.erase_from_index(m_column_index);
}
for (auto it = sparse_vector_iterator<T>(m_column_vector); !it.done(); it.move()) {
unsigned i = it.index();
if (is_zero(w[i])) {
L v = w[i] = w_at_column_index * it.value();
if (settings.abs_val_is_smaller_than_drop_tolerance(v)) {
w[i] = zero_of_type<L>();
continue;
}
w.m_index.push_back(i);
} else {
L v = w[i] += w_at_column_index * it.value(); // todo indexed_vector
if (settings.abs_val_is_smaller_than_drop_tolerance(v)) {
w[i] = zero_of_type<L>();
w.erase_from_index(i);
}
for (auto it = sparse_vector_iterator<T>(m_column_vector); !it.done(); it.move()) {
unsigned i = it.index();
if (is_zero(w[i])) {
L v = w[i] = w_at_column_index * it.value();
if (settings.abs_val_is_smaller_than_drop_tolerance(v)) {
w[i] = zero_of_type<L>();
continue;
}
w.m_index.push_back(i);
} else {
L v = w[i] += w_at_column_index * it.value(); // todo indexed_vector
if (settings.abs_val_is_smaller_than_drop_tolerance(v)) {
w[i] = zero_of_type<L>();
w.erase_from_index(i);
}
}
}
}
void apply_from_left_to_T(indexed_vector<T> & w, lp_settings & settings) {
apply_from_left_local(w, settings);
}
void apply_from_left_to_T(indexed_vector<T> & w, lp_settings & settings) {
apply_from_left_local(w, settings);
}
void push_back(unsigned row_index, T val ) {
lean_assert(row_index != m_column_index);
m_column_vector.push_back(row_index, val);
}
void push_back(unsigned row_index, T val ) {
lean_assert(row_index != m_column_index);
m_column_vector.push_back(row_index, val);
}
void apply_from_right(vector<T> & w) {
void apply_from_right(std::vector<T> & w) {
#ifdef LEAN_DEBUG
// dense_matrix<T, X> deb(*this);
// auto clone_w = clone_vector<T>(w, get_number_of_rows());
// deb.apply_from_right(clone_w);
// dense_matrix<T, X> deb(*this);
// auto clone_w = clone_vector<T>(w, get_number_of_rows());
// deb.apply_from_right(clone_w);
#endif
T t = w[m_column_index] / m_diagonal_element;
for (auto it = sparse_vector_iterator<T>(m_column_vector); !it.done(); it.move()) {
t += w[it.index()] * it.value();
T t = w[m_column_index] / m_diagonal_element;
for (auto it = sparse_vector_iterator<T>(m_column_vector); !it.done(); it.move()) {
t += w[it.index()] * it.value();
}
w[m_column_index] = t;
#ifdef LEAN_DEBUG
// lean_assert(vectors_are_equal<T>(clone_w, w, get_number_of_rows()));
// delete clone_w;
#endif
}
T get_elem(unsigned i, unsigned j) const {
if (j == m_column_index){
if (i == j) {
return 1 / m_diagonal_element;
}
w[m_column_index] = t;
#ifdef LEAN_DEBUG
// lean_assert(vectors_are_equal<T>(clone_w, w, get_number_of_rows()));
// delete clone_w;
#endif
return m_column_vector[i];
}
// void apply_from_right(indexed_vector<T> & w) {
// #ifdef LEAN_DEBUG
// // dense_matrix<T, X> deb(*this);
// // auto clone_w = clone_vector<T>(w, get_number_of_rows());
// // deb.apply_from_right(clone_w);
// #endif
// T t = w[m_column_index] / m_diagonal_element;
// for (auto it = sparse_vector_iterator<T>(m_column_vector); !it.done(); it.move()) {
// t += w[it.index()] * it.value();
// }
// if (numeric_traits<T>::is_zero(w[m_column_index])) {
// w.m_index.push_back(m_column_index);
// }
// w[m_column_index] = t; // we might get a zero here
// #ifdef LEAN_DEBUG
// // lean_assert(vectors_are_equal<T>(clone_w, w, get_number_of_rows()));
// // delete clone_w;
// #endif
// }
return i == j ? numeric_traits<T>::one() : numeric_traits<T>::zero();
}
#ifdef LEAN_DEBUG
unsigned row_count() const { return m_length; }
unsigned column_count() const { return m_length; }
void set_number_of_rows(unsigned /*m*/) { }
void set_number_of_columns(unsigned /*n*/) { }
#endif
sparse_vector_iterator<T> get_sparse_vector_iterator() {
return sparse_vector_iterator<T>(m_column_vector);
}
T get_elem(unsigned i, unsigned j) const {
if (j == m_column_index){
if (i == j) {
return 1 / m_diagonal_element;
}
return m_column_vector[i];
}
return i == j ? numeric_traits<T>::one() : numeric_traits<T>::zero();
void divide_by_diagonal_element() {
m_column_vector.divide(m_diagonal_element);
}
void conjugate_by_permutation(permutation_matrix<T, X> & p) {
// this = p * this * p(-1)
#ifdef LEAN_DEBUG
// auto rev = p.get_reverse();
// auto deb = ((*this) * rev);
// deb = p * deb;
#endif
m_column_index = p.get_rev(m_column_index);
for (auto & pair : m_column_vector.m_data) {
pair.first = p.get_rev(pair.first);
}
#ifdef LEAN_DEBUG
unsigned row_count() const { return m_length; }
unsigned column_count() const { return m_length; }
void set_number_of_rows(unsigned /*m*/) { }
void set_number_of_columns(unsigned /*n*/) { }
#endif
sparse_vector_iterator<T> get_sparse_vector_iterator() {
return sparse_vector_iterator<T>(m_column_vector);
}
void divide_by_diagonal_element() {
m_column_vector.divide(m_diagonal_element);
}
void conjugate_by_permutation(permutation_matrix<T, X> & p) {
// this = p * this * p(-1)
#ifdef LEAN_DEBUG
// auto rev = p.get_reverse();
// auto deb = ((*this) * rev);
// deb = p * deb;
#endif
m_column_index = p.get_rev(m_column_index);
for (auto & pair : m_column_vector.m_data) {
pair.first = p.get_rev(pair.first);
}
#ifdef LEAN_DEBUG
// lean_assert(deb == *this);
// lean_assert(deb == *this);
#endif
}
};
}
};
}

16
src/util/lp/indexed_value.h
View file

@ -46,13 +46,13 @@ public:
m_value = val;
}
};
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;
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;
}
}

21
src/util/lp/indexed_vector.h
View file

@ -20,11 +20,6 @@
#include "util/lp/sparse_vector.h"
#include <iomanip>
namespace lean {
using std::string;
using std::cout;
using std::endl;
using std::vector;
template <typename T>
void print_vector(const T * t, unsigned l) {
for (unsigned i = 0; i < l; i++)
@ -33,7 +28,7 @@ void print_vector(const T * t, unsigned l) {
}
template <typename T>
void print_vector(const vector<T> & t) {
void print_vector(const std::vector<T> & t) {
for (unsigned i = 0; i < t.size(); i++)
std::cout << t[i] << " ";
std::cout << std::endl;
@ -47,7 +42,7 @@ void print_vector(const buffer<T> & t) {
}
template <typename T>
void print_sparse_vector(const vector<T> & t) {
void print_sparse_vector(const std::vector<T> & t) {
for (unsigned i = 0; i < t.size(); i++) {
if (is_zero(t[i]))continue;
std::cout << "[" << i << "] = " << t[i] << ", ";
@ -55,7 +50,7 @@ void print_sparse_vector(const vector<T> & t) {
std::cout << std::endl;
}
void print_vector(const vector<mpq> & t) {
void print_vector(const std::vector<mpq> & t) {
for (unsigned i = 0; i < t.size(); i++)
std::cout << t[i].get_double() << std::setprecision(3) << " ";
std::cout << std::endl;
@ -66,7 +61,7 @@ class indexed_vector {
public:
// m_index points to non-zero elements of m_data
buffer<T> m_data;
vector<unsigned> m_index;
std::vector<unsigned> m_index;
indexed_vector(unsigned data_size) {
m_data.resize(data_size, numeric_traits<T>::zero());
}
@ -117,8 +112,6 @@ public:
if (it != m_index.end()) m_index.erase(it);
}
#ifdef LEAN_DEBUG
bool is_OK() const {
int size = 0;
@ -132,11 +125,11 @@ public:
return size == m_index.size();
}
void print() {
cout << "m_index " << endl;
std::cout << "m_index " << std::endl;
for (unsigned i = 0; i < m_index.size(); i++) {
cout << m_index[i] << " ";
std::cout << m_index[i] << " ";
}
cout << endl;
std::cout << std::endl;
print_vector(m_data);
}
#endif

181
src/util/lp/lar_constraints.h
View file

@ -16,101 +16,88 @@
#include <algorithm>
#include "util/lp/canonic_left_side.h"
namespace lean {
lconstraint_kind flip_kind(lconstraint_kind t) {
return static_cast<lconstraint_kind>( - static_cast<int>(t));
}
std::string lconstraint_kind_string(lconstraint_kind t) {
switch (t) {
case LE: return string("<=");
case LT: return string("<");
case GE: return string(">=");
case GT: return string(">");
case EQ: return string("=");
default:
throw "unexpected";
}
}
class lar_base_constraint {
public:
lconstraint_kind m_kind;
mpq m_right_side;
virtual buffer<pair<mpq, var_index>> get_left_side_coefficients() const = 0;
constraint_index m_index; // the index of constraint
lar_base_constraint() {}
lar_base_constraint(lconstraint_kind kind, mpq right_side, constraint_index index) :m_kind(kind), m_right_side(right_side), m_index(index) {}
virtual unsigned size() const = 0;
virtual ~lar_base_constraint(){}
};
class lar_constraint : public lar_base_constraint {
public:
std::unordered_map<var_index, mpq> m_left_side;
lar_constraint() {}
lar_constraint(const buffer<pair<mpq, var_index>> & left_side, lconstraint_kind kind, mpq right_side, constraint_index index) : lar_base_constraint(kind, right_side, index) {
for (auto & it : left_side) {
auto r = m_left_side.find(it.second);
if (r == m_left_side.end()) {
m_left_side[it.second] = it.first;
} else {
r->second += it.first;
}
}
}
lar_constraint(const lar_base_constraint & c): lar_base_constraint(c.m_kind, c.m_right_side, c.m_index) {
for (auto t : c.get_left_side_coefficients())
m_left_side.insert(std::make_pair(t.second, t.first));
}
unsigned size() const {
return m_left_side.size();
}
buffer<pair<mpq, var_index>> get_left_side_coefficients() const {
buffer<pair<mpq, var_index>> ret;
for (auto it : m_left_side) {
ret.push_back(std::make_pair(it.second, it.first));
}
return ret;
}
};
class lar_normalized_constraint : public lar_base_constraint {
public:
canonic_left_side* m_canonic_left_side;
mpq m_ratio_to_original; // by multiplying this constraint by m_ratio_to_original we get the original one
lar_constraint m_origin_constraint;
lar_normalized_constraint(canonic_left_side * ls, mpq ratio, lconstraint_kind kind, mpq right_side, const lar_constraint & origin):
lar_base_constraint(kind, right_side, origin.m_index),
m_canonic_left_side(ls),
m_ratio_to_original(ratio),
m_origin_constraint(origin) {
}
lar_normalized_constraint() {}
// lar_normalized_constraint & operator=(lar_normalized_constraint & other) {
// m_canonic_left_side = other.m_canonic_left_side;
// m_origin_constraint = other.m_origin_constraint;
// m_kind = other.m_kind;
// m_right_side = other.m_right_side;
// return *this;
// }
virtual buffer<pair<mpq, var_index>> get_left_side_coefficients() const {
buffer<pair<mpq, var_index>> ret;
for (auto t : m_canonic_left_side->m_coeffs) ret.push_back(t);
return ret;
}
virtual unsigned size() const {
return m_canonic_left_side->size();
}
};
lconstraint_kind flip_kind(lconstraint_kind t) {
return static_cast<lconstraint_kind>( - static_cast<int>(t));
}
std::string lconstraint_kind_string(lconstraint_kind t) {
switch (t) {
case LE: return std::string("<=");
case LT: return std::string("<");
case GE: return std::string(">=");
case GT: return std::string(">");
case EQ: return std::string("=");
}
lean_unreachable();
}
class lar_base_constraint {
public:
lconstraint_kind m_kind;
mpq m_right_side;
virtual buffer<pair<mpq, var_index>> get_left_side_coefficients() const = 0;
constraint_index m_index; // the index of constraint
lar_base_constraint() {}
lar_base_constraint(lconstraint_kind kind, mpq right_side, constraint_index index) :m_kind(kind), m_right_side(right_side), m_index(index) {}
virtual unsigned size() const = 0;
virtual ~lar_base_constraint(){}
};
class lar_constraint : public lar_base_constraint {
public:
std::unordered_map<var_index, mpq> m_left_side;
lar_constraint() {}
lar_constraint(const buffer<pair<mpq, var_index>> & left_side, lconstraint_kind kind, mpq right_side, constraint_index index) : lar_base_constraint(kind, right_side, index) {
for (auto & it : left_side) {
auto r = m_left_side.find(it.second);
if (r == m_left_side.end()) {
m_left_side[it.second] = it.first;
} else {
r->second += it.first;
}
}
}
lar_constraint(const lar_base_constraint & c): lar_base_constraint(c.m_kind, c.m_right_side, c.m_index) {
for (auto t : c.get_left_side_coefficients())
m_left_side.insert(std::make_pair(t.second, t.first));
}
unsigned size() const {
return m_left_side.size();
}
buffer<pair<mpq, var_index>> get_left_side_coefficients() const {
buffer<pair<mpq, var_index>> ret;
for (auto it : m_left_side) {
ret.push_back(std::make_pair(it.second, it.first));
}
return ret;
}
};
class lar_normalized_constraint : public lar_base_constraint {
public:
canonic_left_side* m_canonic_left_side;
mpq m_ratio_to_original; // by multiplying this constraint by m_ratio_to_original we get the original one
lar_constraint m_origin_constraint;
lar_normalized_constraint(canonic_left_side * ls, mpq ratio, lconstraint_kind kind, mpq right_side, const lar_constraint & origin):
lar_base_constraint(kind, right_side, origin.m_index),
m_canonic_left_side(ls),
m_ratio_to_original(ratio),
m_origin_constraint(origin) {
}
lar_normalized_constraint() {}
virtual buffer<pair<mpq, var_index>> get_left_side_coefficients() const {
buffer<pair<mpq, var_index>> ret;
for (auto t : m_canonic_left_side->m_coeffs) ret.push_back(t);
return ret;
}
virtual unsigned size() const {
return m_canonic_left_side->size();
}
};
}

1570
src/util/lp/lar_core_solver.h
View file

File diff suppressed because it is too large Load diff

8
src/util/lp/lar_solution_signature.h
View file

@ -11,8 +11,8 @@
#include "util/lp/lp_settings.h"
#include <unordered_map>
namespace lean {
struct lar_solution_signature {
std::unordered_map<unsigned, non_basic_column_value_position> non_basic_column_value_positions;
lp_status status;
};
struct lar_solution_signature {
std::unordered_map<unsigned, non_basic_column_value_position> non_basic_column_value_positions;
lp_status status;
};
}

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

@ -4,7 +4,6 @@
Author: Lev Nachmanson
*/
#pragma once
#include <vector>
#include "util/debug.h"
@ -23,77 +22,68 @@
#include "util/lp/scaler.h"
namespace std {
template<>
struct hash<lean::mpq> {
inline size_t operator()(const lean::mpq & v) const {
return v.hash();
template<>
struct hash<lean::mpq> {
inline size_t operator()(const lean::mpq & v) const {
return v.hash();
}
};
}
namespace lean {
using std::string;
using std::cout;
using std::endl;
using std::vector;
using std::unordered_map;
using std::unordered_set;
template <typename V>
struct conversion_helper {
static V get_low_bound(const column_info<mpq> & ci) {
return V(ci.get_low_bound(), ci.low_bound_is_strict()? 1 : 0);
}
template <typename V>
struct conversion_helper {
static V get_low_bound(const column_info<mpq> & ci) {
return V(ci.get_low_bound(), ci.low_bound_is_strict()? 1 : 0);
}
static V get_upper_bound(const column_info<mpq> & ci) {
return V(ci.get_upper_bound(), ci.upper_bound_is_strict()? -1 : 0);
}
};
static V get_upper_bound(const column_info<mpq> & ci) {
return V(ci.get_upper_bound(), ci.upper_bound_is_strict()? -1 : 0);
}
};
template<>
struct conversion_helper <double> {
static double get_low_bound(const column_info<mpq> & ci) {
if (!ci.low_bound_is_strict())
return ci.get_low_bound().get_double();
double eps = 0.00001;
if (!ci.upper_bound_is_set())
return ci.get_low_bound().get_double() + eps;
eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double()/1000, eps);
template<>
struct conversion_helper <double> {
static double get_low_bound(const column_info<mpq> & ci) {
if (!ci.low_bound_is_strict())
return ci.get_low_bound().get_double();
double eps = 0.00001;
if (!ci.upper_bound_is_set())
return ci.get_low_bound().get_double() + eps;
}
static double get_upper_bound(const column_info<mpq> & ci) {
if (!ci.upper_bound_is_strict())
return ci.get_upper_bound().get_double();
double eps = 0.00001;
if (!ci.low_bound_is_set())
return ci.get_upper_bound().get_double() - eps;
eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double()/1000, eps);
eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double()/1000, eps);
return ci.get_low_bound().get_double() + eps;
}
static double get_upper_bound(const column_info<mpq> & ci) {
if (!ci.upper_bound_is_strict())
return ci.get_upper_bound().get_double();
double eps = 0.00001;
if (!ci.low_bound_is_set())
return ci.get_upper_bound().get_double() - eps;
}
};
eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double()/1000, eps);
return ci.get_upper_bound().get_double() - eps;
}
};
class lar_solver {
unsigned m_available_var_index = 0;
unsigned m_available_constr_index = 0;
lp_status m_status = UNKNOWN;
unordered_map<string, var_index> m_var_names_to_var_index;
unordered_map<var_index, canonic_left_side*> m_map_from_var_index_to_left_side;
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;
unordered_map<unsigned, canonic_left_side*> m_column_indices_to_canonic_left_sides;
unordered_map<constraint_index, lar_normalized_constraint> m_normalized_constraints;
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_canonic_left_sides;
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;
static_matrix<mpq, numeric_pair<mpq>> m_A;
lar_core_solver<mpq, numeric_pair<mpq>> m_mpq_core_solver;
lp_settings m_settings;
vector<column_type> m_column_types; // this field is passed to the core solver
vector<numeric_pair<mpq>> m_low_bounds; // passed to the core solver
vector<numeric_pair<mpq>> m_upper_bounds; // passed to the core solver
vector<unsigned> m_basis;
vector<numeric_pair<mpq>> m_x; // the solution
unordered_map<unsigned, string> m_column_names;
vector<numeric_pair<mpq>> m_right_side_vector; // this vector will be all zeroes, it might change when the optimization with fixed variables will used
vector<mpq> m_costs;
std::vector<column_type> m_column_types; // this field is passed to the core solver
std::vector<numeric_pair<mpq>> m_low_bounds; // passed to the core solver
std::vector<numeric_pair<mpq>> m_upper_bounds; // passed to the core solver
std::vector<unsigned> m_basis;
std::vector<numeric_pair<mpq>> m_x; // the solution
std::unordered_map<unsigned, std::string> m_column_names;
std::vector<numeric_pair<mpq>> m_right_side_vector; // this vector will be all zeroes, it might change when the optimization with fixed variables will used
std::vector<mpq> m_costs;
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<pair<mpq, var_index>>& left_side_par) {
auto left_side = new canonic_left_side(left_side_par);
@ -125,7 +115,7 @@ class lar_solver {
return it->second;
}
void add_canonic_left_side_for_var(var_index i, string var_name) {
void add_canonic_left_side_for_var(var_index i, std::string var_name) {
buffer<pair<mpq, var_index>> b;
b.push_back(std::make_pair(numeric_traits<mpq>::one(), i));
auto can_ls = new canonic_left_side(b);
@ -265,11 +255,11 @@ class lar_solver {
break;
case EQ:
{
set_upper_bound_for_column_info(&norm_constr);
set_low_bound_for_column_info(&norm_constr);
}
break;
{
set_upper_bound_for_column_info(&norm_constr);
set_low_bound_for_column_info(&norm_constr);
}
break;
default:
throw "unexpected";
}
@ -290,9 +280,9 @@ class lar_solver {
auto & ci = t->m_column_info;
unsigned j = t->m_column_index;
lean_assert(valid_index(j));
string name = ci.get_name();
std::string name = ci.get_name();
if (name.size() == 0)
name = string("_s") + T_to_string(j);
name = std::string("_s") + T_to_string(j);
m_column_names[j] = name;
}
}
@ -309,7 +299,7 @@ class lar_solver {
}
template <typename V>
void fill_bounds_for_core_solver(vector<V> & lb, vector<V> & ub) {
void fill_bounds_for_core_solver(std::vector<V> & lb, std::vector<V> & ub) {
unsigned n = m_canonic_left_sides.size(); // this is the number of columns
lb.resize(n);
ub.resize(n);
@ -326,14 +316,14 @@ class lar_solver {
template <typename V>
void resize_x_and_init_with_zeros(vector<V> & x, unsigned n) {
void resize_x_and_init_with_zeros(std::vector<V> & x, unsigned n) {
x.clear();
x.resize(n, zero_of_type<V>()); // init with zeroes
}
template <typename V>
void resize_x_and_init_with_signature(vector<V> & x, vector<V> & low_bound,
vector<V> & upper_bound, const lar_solution_signature & signature) {
void resize_x_and_init_with_signature(std::vector<V> & x, std::vector<V> & low_bound,
std::vector<V> & upper_bound, const lar_solution_signature & signature) {
x.clear();
x.resize(low_bound.size());
for (auto & t : signature.non_basic_column_value_positions) {
@ -341,7 +331,7 @@ class lar_solver {
}
}
template <typename V> V get_column_val(vector<V> & low_bound, 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) {
switch (pos_type) {
case at_low_bound: return low_bound[j];
case at_fixed:
@ -354,7 +344,7 @@ class lar_solver {
public:
~lar_solver() {
vector<canonic_left_side*> to_delete;
std::vector<canonic_left_side*> to_delete;
for (auto it : m_canonic_left_sides)
to_delete.push_back(it);
for (auto t : to_delete)
@ -380,7 +370,7 @@ public:
}
var_index add_var(string s) {
var_index add_var(std::string s) {
auto got = m_var_names_to_var_index.find(s);
if (got != m_var_names_to_var_index.end()) return got->second;
@ -392,7 +382,7 @@ public:
constraint_index add_constraint(const buffer<pair<mpq, var_index>>& left_side, lconstraint_kind kind_par, mpq right_side_par) {
if (left_side.size() == 0) {
cout << "cannot add a constraint without left side" << endl;
std::cout << "cannot add a constraint without left side" << std::endl;
return (constraint_index)(-1);
}
constraint_index i = m_available_constr_index++;
@ -412,18 +402,18 @@ public:
}
bool all_constraints_hold() {
unordered_map<var_index, mpq> var_map;
std::unordered_map<var_index, mpq> var_map;
get_model(var_map);
for ( auto & it : m_normalized_constraints )
if (!constraint_holds(it.second.m_origin_constraint, var_map)) {
print_constraint(&it.second.m_origin_constraint);
cout << endl;
std::cout << std::endl;
return false;
}
return true;
}
bool constraint_holds(const lar_constraint & constr, unordered_map<var_index, mpq> & var_map) {
bool constraint_holds(const lar_constraint & constr, std::unordered_map<var_index, mpq> & var_map) {
mpq left_side_val = get_left_side_val(constr, var_map);
switch (constr.m_kind) {
case LE: return left_side_val <= constr.m_right_side;
@ -453,9 +443,9 @@ public:
for (auto & it : evidence) {
mpq coeff = it.first;
constraint_index con_ind = it.second;
cout << coeff.get_double() << endl;
std::cout << coeff.get_double() << std::endl;
lar_constraint & constr = m_normalized_constraints[con_ind].m_origin_constraint;
print_constraint(&constr); cout << endl;
print_constraint(&constr); std::cout << std::endl;
lconstraint_kind kind = coeff.is_pos()? constr.m_kind: flip_kind(constr.m_kind);
if (kind == GT || kind == LT)
@ -470,7 +460,7 @@ public:
return n_of_G == 0 || n_of_L == 0;
}
void register_in_map(unordered_map<var_index, mpq> & coeffs, lar_constraint & cn, const mpq & a) {
void register_in_map(std::unordered_map<var_index, mpq> & coeffs, lar_constraint & cn, const mpq & a) {
for (auto & it : cn.m_left_side) {
unsigned j = it.first;
auto p = coeffs.find(j);
@ -479,7 +469,7 @@ public:
}
}
bool the_left_sides_sum_to_zero(const buffer<pair<mpq, unsigned>> & evidence) {
unordered_map<var_index, mpq> coeff_map;
std::unordered_map<var_index, mpq> coeff_map;
for (auto & it : evidence) {
mpq coeff = it.first;
constraint_index con_ind = it.second;
@ -533,7 +523,7 @@ public:
}
template <typename V>
void init_right_sides_with_zeros(vector<V> & rs) {
void init_right_sides_with_zeros(std::vector<V> & rs) {
rs.clear();
rs.resize(m_basis.size(), zero_of_type<V>());
}
@ -556,10 +546,10 @@ public:
}
template <typename U, typename V>
void prepare_core_solver_fields(static_matrix<U, V> & A, vector<V> & x,
vector<V> & right_side_vector,
vector<V> & low_bound,
vector<V> & upper_bound) {
void prepare_core_solver_fields(static_matrix<U, V> & A, std::vector<V> & x,
std::vector<V> & right_side_vector,
std::vector<V> & low_bound,
std::vector<V> & upper_bound) {
create_matrix_A(A);
fill_bounds_for_core_solver(low_bound, upper_bound);
if (m_status == INFEASIBLE) {
@ -571,10 +561,10 @@ public:
}
template <typename U, typename V>
void prepare_core_solver_fields_with_signature(static_matrix<U, V> & A, vector<V> & x,
vector<V> & right_side_vector,
vector<V> & low_bound,
vector<V> & upper_bound, const lar_solution_signature & signature) {
void prepare_core_solver_fields_with_signature(static_matrix<U, V> & A, std::vector<V> & x,
std::vector<V> & right_side_vector,
std::vector<V> & low_bound,
std::vector<V> & upper_bound, const lar_solution_signature & signature) {
create_matrix_A(A);
fill_bounds_for_core_solver(low_bound, upper_bound);
if (m_status == INFEASIBLE) {
@ -586,9 +576,9 @@ public:
void find_solution_signature_with_doubles(lar_solution_signature & signature) {
static_matrix<double, double> A;
vector<double> x, right_side_vector, low_bounds, upper_bounds;
std::vector<double> x, right_side_vector, low_bounds, upper_bounds;
prepare_core_solver_fields<double, double>(A, x, right_side_vector, low_bounds, upper_bounds);
vector<double> column_scale_vector;
std::vector<double> column_scale_vector;
scaler<double, double > scaler(right_side_vector, A, m_settings.scaling_minimum, m_settings.scaling_maximum, column_scale_vector, this->m_settings);
if (!scaler.scale()) {
// the scale did not succeed, unscaling
@ -597,7 +587,7 @@ public:
for (auto & s : column_scale_vector)
s = one_of_type<double>();
}
vector<double> costs(A.column_count());
std::vector<double> costs(A.column_count());
auto core_solver = lp_primal_core_solver<double, double>(A,
right_side_vector,
x,
@ -643,7 +633,7 @@ public:
}
void get_infeasibility_evidence(buffer<pair<mpq, constraint_index>> & evidence){
if (!m_mpq_core_solver.get_infeasible_row_sign()) {
cout << "don't have the infeasibility evidence" << endl;
std::cout << "don't have the infeasibility evidence" << std::endl;
return;
}
// the infeasibility sign
@ -654,7 +644,7 @@ public:
}
void get_infeasibility_evidence_for_inf_sign(buffer<pair<mpq, constraint_index>> & evidence,
const vector<pair<mpq, unsigned>> & inf_row,
const std::vector<pair<mpq, unsigned>> & inf_row,
int inf_sign) {
for (auto & it : inf_row) {
mpq coeff = it.first;
@ -714,19 +704,19 @@ public:
}
}
void get_model(unordered_map<var_index, mpq> & variable_values){
void get_model(std::unordered_map<var_index, mpq> & variable_values){
lean_assert(m_status == OPTIMAL);
mpq delta = find_delta_for_strict_bounds();
for (auto & it : m_map_from_var_index_to_left_side) {
numeric_pair<mpq> & rp = m_x[it.second->m_column_index];
// cout << it.second->m_column_info.get_name() << " = " << rp << endl;
// std::cout << it.second->m_column_info.get_name() << " = " << rp << std::endl;
variable_values[it.first] = rp.x + delta * rp.y;
}
}
string get_variable_name(var_index vi) {
std::string get_variable_name(var_index vi) {
if (m_map_from_var_index_to_left_side.size() <= vi) {
string s="variable " + T_to_string(vi) + " is not found";
std::string s = "variable " + T_to_string(vi) + " is not found";
return s;
}
return m_map_from_var_index_to_left_side[vi]->m_column_info.get_name();
@ -735,8 +725,8 @@ public:
// ********** print region start
void print_constraint(constraint_index ci) {
if (m_normalized_constraints.size() <= ci) {
string s = "constraint " + T_to_string(ci) + " is not found";
cout << s << endl;
std::string s = "constraint " + T_to_string(ci) + " is not found";
std::cout << s << std::endl;
return;
}
@ -751,15 +741,15 @@ public:
first = false;
} else {
if (val.is_pos()) {
cout << " + ";
std::cout << " + ";
} else {
cout << " - ";
std::cout << " - ";
val = -val;
}
}
if (val != numeric_traits<mpq>::one())
cout << T_to_string(val);
cout << m_map_from_var_index_to_left_side[it.second]->m_column_info.get_name();
std::cout << T_to_string(val);
std::cout << m_map_from_var_index_to_left_side[it.second]->m_column_info.get_name();
}
}
@ -772,21 +762,20 @@ public:
first = false;
} else {
if (val.is_pos()) {
cout << " + ";
std::cout << " + ";
} else {
cout << " - ";
std::cout << " - ";
val = -val;
}
}
if (val != numeric_traits<mpq>::one())
cout << val;
cout << m_map_from_var_index_to_left_side[it.second]->m_column_info.get_name();
std::cout << val;
std::cout << m_map_from_var_index_to_left_side[it.second]->m_column_info.get_name();
}
}
numeric_pair<mpq> get_infeasibility_from_core_solver(unordered_map<string, mpq> & solution) {
numeric_pair<mpq> get_infeasibility_from_core_solver(std::unordered_map<std::string, mpq> & solution) {
prepare_independently_of_numeric_type();
prepare_core_solver_fields(m_A, m_x, m_right_side_vector, m_low_bounds, m_upper_bounds);
m_mpq_core_solver.prefix(); // just to fill the core solver
@ -803,34 +792,34 @@ public:
if (static_cast<unsigned>(ls->m_column_index) == j) {
auto & ci = ls->m_column_info;
if (ci.low_bound_is_set()) {
cout << "l = " << ci.get_low_bound();
std::cout << "l = " << ci.get_low_bound();
}
if (ci.upper_bound_is_set()) {
cout << "u = " << ci.get_upper_bound();
std::cout << "u = " << ci.get_upper_bound();
}
cout << endl;
std::cout << std::endl;
m_mpq_core_solver.print_column_info(j);
}
}
}
mpq get_infeasibility_of_solution(unordered_map<string, mpq> & solution) {
cout << "solution" << endl;
mpq get_infeasibility_of_solution(std::unordered_map<std::string, mpq> & solution) {
std::cout << "solution" << std::endl;
for (auto it : solution) {
cout << it.first << " = " << it.second.get_double() << endl;
std::cout << it.first << " = " << it.second.get_double() << std::endl;
}
mpq ret = numeric_traits<mpq>::zero();
for (auto it : m_normalized_constraints) {
ret += get_infeasibility_of_constraint(it.second, solution);
}
cout << "ret = " << ret.get_double() << endl;
std::cout << "ret = " << ret.get_double() << std::endl;
auto core_inf = get_infeasibility_from_core_solver(solution);
cout << "core inf = " << T_to_string(core_inf) << endl;
std::cout << "core inf = " << T_to_string(core_inf) << std::endl;
lean_assert(numeric_pair<mpq>(ret, 0) == core_inf);
return ret;
}
mpq get_infeasibility_of_constraint(const lar_normalized_constraint & norm_constr, unordered_map<string, mpq> & solution) {
mpq get_infeasibility_of_constraint(const lar_normalized_constraint & norm_constr, std::unordered_map<std::string, mpq> & solution) {
auto kind = norm_constr.m_kind;
mpq left_side_val = get_canonic_left_side_val(norm_constr.m_canonic_left_side, solution);
@ -849,14 +838,14 @@ public:
}
}
mpq get_canonic_left_side_val(canonic_left_side * ls, unordered_map<string, mpq> & solution) {
mpq get_canonic_left_side_val(canonic_left_side * ls, std::unordered_map<std::string, mpq> & solution) {
mpq ret = numeric_traits<mpq>::zero();
for (auto it : ls->m_coeffs) {
var_index j = it.second;
auto vi = m_map_from_var_index_to_left_side.find(j);
lean_assert(vi != m_map_from_var_index_to_left_side.end());
canonic_left_side * var_ls = vi->second;
string s = var_ls->m_column_info.get_name();
std::string s = var_ls->m_column_info.get_name();
auto t = solution.find(s);
lean_assert(t != solution.end());
ret += it.first * (t->second);
@ -864,7 +853,7 @@ public:
return ret;
}
mpq get_left_side_val(const lar_constraint & cns, const unordered_map<var_index, mpq> & var_map) {
mpq get_left_side_val(const lar_constraint & cns, const std::unordered_map<var_index, mpq> & var_map) {
mpq ret = numeric_traits<mpq>::zero();
for (auto it : cns.m_left_side) {
var_index j = it.first;
@ -877,7 +866,7 @@ public:
void print_constraint(const lar_base_constraint * c) {
print_left_side_of_constraint(c);
cout <<" " << lconstraint_kind_string(c->m_kind) << " " << c->m_right_side;
std::cout <<" " << lconstraint_kind_string(c->m_kind) << " " << c->m_right_side;
}
unsigned get_total_iterations() const { return m_mpq_core_solver.m_total_iterations; }
};

122
src/util/lp/lp_core_solver_base.h
View file

@ -4,23 +4,24 @@
Author: Lev Nachmanson
*/
#pragma once
#include <string>
#include "util/lp/lp.h"
#include "util/lp/core_solver_pretty_printer.h"
#include <set>
#include <vector>
#include <string>
#include "util/sstream.h"
#include "util/exception.h"
#include "util/lp/lp.h"
#include "util/lp/core_solver_pretty_printer.h"
#include "util/lp/numeric_pair.h"
namespace lean {
void init_basic_part_of_basis_heading(std::vector<unsigned> & basis, unsigned m, vector<int> & basis_heading) {
void init_basic_part_of_basis_heading(std::vector<unsigned> & basis, unsigned m, std::vector<int> & basis_heading) {
for (unsigned i = 0; i < m; i++) {
unsigned column = basis[i];
basis_heading[column] = i;
}
}
void init_non_basic_part_of_basis_heading(vector<int> & basis_heading, vector<unsigned> & non_basic_columns, unsigned n) {
void init_non_basic_part_of_basis_heading(std::vector<int> & basis_heading, std::vector<unsigned> & non_basic_columns, unsigned n) {
for (int j = n; j--;){
if (basis_heading[j] < 0) {
non_basic_columns.push_back(j);
@ -30,10 +31,10 @@ void init_non_basic_part_of_basis_heading(vector<int> & basis_heading, vector<un
}
}
void init_basis_heading_and_non_basic_columns_vector(std::vector<unsigned> & basis,
unsigned m,
vector<int> & basis_heading,
unsigned n,
vector<unsigned> & non_basic_columns) {
unsigned m,
std::vector<int> & basis_heading,
unsigned n,
std::vector<unsigned> & non_basic_columns) {
init_basic_part_of_basis_heading(basis, m, basis_heading);
init_non_basic_part_of_basis_heading(basis_heading, non_basic_columns, n);
}
@ -45,9 +46,9 @@ public:
unsigned m_n; // the number of columns in the matrix m_A
std::vector<T> m_pivot_row_of_B_1; // the pivot row of the reverse of B
std::vector<T> m_pivot_row; // this is the real pivot row of the simplex tableu
vector<unsigned> m_pivot_row_index;
std::vector<unsigned> m_pivot_row_index;
static_matrix<T, X> & m_A; // the matrix A
vector<X> & m_b; // the right side
std::vector<X> & m_b; // the right side
std::vector<unsigned> & m_basis;
std::vector<int> m_basis_heading;
std::vector<X> & m_x; // a feasible solution, the fist time set in the constructor
@ -57,31 +58,31 @@ public:
lp_status m_status;
// a device that is able to solve Bx=c, xB=d, and change the basis
lu<T, X> * m_factorization;
const unordered_map<unsigned, string> & m_column_names;
const std::unordered_map<unsigned, std::string> & m_column_names;
indexed_vector<T> m_w; // the vector featuring in 24.3 of the Chvatal book
std::vector<T> m_d; // the vector of reduced costs
std::vector<T> m_ed; // the solution of B*m_ed = a
vector<unsigned> m_index_of_ed;
std::vector<unsigned> m_index_of_ed;
unsigned m_total_iterations = 0;
int m_start_time;
unsigned m_iters_with_no_cost_growing = 0;
vector<unsigned> m_non_basic_columns;
vector<column_type> & m_column_type;
vector<X> & m_low_bound_values;
vector<X> & m_upper_bound_values;
vector<T> m_column_norms; // the approximate squares of column norms that help choosing a profitable column
vector<X> m_copy_of_xB;
std::vector<unsigned> m_non_basic_columns;
std::vector<column_type> & m_column_type;
std::vector<X> & m_low_bound_values;
std::vector<X> & m_upper_bound_values;
std::vector<T> m_column_norms; // the approximate squares of column norms that help choosing a profitable column
std::vector<X> m_copy_of_xB;
unsigned m_refactor_counter = 200;
lp_core_solver_base(static_matrix<T, X> & A,
vector<X> & b, // the right side vector
std::vector<X> & b, // the right side vector
std::vector<unsigned> & basis,
std::vector<X> & x,
std::vector<T> & costs,
lp_settings & settings,
const unordered_map<unsigned, string> & column_names,
vector<column_type> & column_types,
vector<X> & low_bound_values,
vector<X> & upper_bound_values):
const std::unordered_map<unsigned, std::string> & column_names,
std::vector<column_type> & column_types,
std::vector<X> & low_bound_values,
std::vector<X> & upper_bound_values):
m_m(A.row_count()),
m_n(A.column_count()),
m_pivot_row_of_B_1(m_m),
@ -130,7 +131,7 @@ public:
delete m_factorization;
}
vector<unsigned> & non_basis() {
std::vector<unsigned> & non_basis() {
return m_factorization->m_non_basic_columns;
}
@ -242,11 +243,11 @@ public:
X eps = feps * (one + T(0.1) * abs(m_b[i]));
if (delta >eps) {
cout << "x is off (";
cout << "m_b[" << i << "] = " << m_b[i] << " ";
cout << "left side = " << m_A.dot_product_with_row(i, m_x) << ' ';
cout << "delta = " << delta << ' ';
cout << "iters = " << m_total_iterations << ")" << endl;
std::cout << "x is off (";
std::cout << "m_b[" << i << "] = " << m_b[i] << " ";
std::cout << "left side = " << m_A.dot_product_with_row(i, m_x) << ' ';
std::cout << "delta = " << delta << ' ';
std::cout << "iters = " << m_total_iterations << ")" << std::endl;
return true;
}
}
@ -307,13 +308,13 @@ public:
}
void print_statistics(X cost) {
cout << "cost = " << T_to_string(cost) <<
", nonzeros = " << m_factorization->get_number_of_nonzeroes() << endl;
std::cout << "cost = " << T_to_string(cost) <<
", nonzeros = " << m_factorization->get_number_of_nonzeroes() << std::endl;
}
bool print_statistics_with_iterations_and_check_that_the_time_is_over(unsigned total_iterations) {
if (total_iterations % m_settings.report_frequency == 0) {
cout << "iterations = " << total_iterations << ", nonzeros = " << m_factorization->get_number_of_nonzeroes() << endl;
std::cout << "iterations = " << total_iterations << ", nonzeros = " << m_factorization->get_number_of_nonzeroes() << std::endl;
if (time_is_over()) {
return true;
}
@ -321,9 +322,9 @@ public:
return false;
}
bool print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over(string str, unsigned total_iterations) {
bool print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over(std::string str, unsigned total_iterations) {
if (total_iterations % m_settings.report_frequency == 0) {
cout << str << " iterations = " << total_iterations << " cost = " << T_to_string(get_cost()) <<", nonzeros = " << m_factorization->get_number_of_nonzeroes() << endl;
std::cout << str << " iterations = " << total_iterations << " cost = " << T_to_string(get_cost()) <<", nonzeros = " << m_factorization->get_number_of_nonzeroes() << std::endl;
if (time_is_over()) {
return true;
}
@ -333,7 +334,7 @@ public:
bool print_statistics_with_cost_and_check_that_the_time_is_over(unsigned total_iterations, X cost) {
if (total_iterations % m_settings.report_frequency == 0) {
cout << "iterations = " << total_iterations << ", ";
std::cout << "iterations = " << total_iterations << ", ";
print_statistics(cost);
if (time_is_over()) {
return true;
@ -344,7 +345,7 @@ public:
bool print_statistics_and_check_that_the_time_is_over(unsigned total_iterations) {
if (total_iterations % (numeric_traits<T>::precise()? static_cast<unsigned>(m_settings.report_frequency/10) : m_settings.report_frequency) == 0) {
cout << "iterations = " << total_iterations << ", ";
std::cout << "iterations = " << total_iterations << ", ";
if (time_is_over()) {
return true;
}
@ -414,16 +415,14 @@ public:
case low_bound:
return x_is_at_low_bound(j) && d_is_not_negative(j);
case upper_bound:
cout << "upper_bound type should be switched to low_bound" << endl;
std::cout << "upper_bound type should be switched to low_bound" << std::endl;
lean_assert(false); // impossible case
case free_column:
return numeric_traits<X>::is_zero(m_d[j]);
default:
cout << "column = " << j << endl;
cout << "unexpected column type = " << column_type_to_string(m_column_type[j]) << endl;
lean_assert(false);
throw "unexpected column type";
break;
std::cout << "column = " << j << std::endl;
std::cout << "unexpected column type = " << column_type_to_string(m_column_type[j]) << std::endl;
lean_unreachable();
}
}
bool d_is_not_negative(unsigned j) {
@ -444,13 +443,13 @@ public:
bool time_is_over() {
int span_in_mills = get_millisecond_span(m_start_time);
if (span_in_mills / 1000.0 > m_settings.time_limit) {
cout << "time is over" << endl;
std::cout << "time is over" << std::endl;
return true;
}
return false;
}
void rs_minus_Anx(vector<X> & rs) {
void rs_minus_Anx(std::vector<X> & rs) {
unsigned row = m_m;
while (row--) {
auto &rsv = rs[row] = m_b[row];
@ -467,9 +466,9 @@ public:
solve_Ax_eq_b();
bool ret= !A_mult_x_is_off();
if (ret)
cout << "find_x_by_solving succeeded" << endl;
std::cout << "find_x_by_solving succeeded" << std::endl;
else
cout << "find_x_by_solving did not succeed" << endl;
std::cout << "find_x_by_solving did not succeed" << std::endl;
return ret;
}
@ -487,8 +486,7 @@ public:
m_refactor_counter = 0;
m_iters_with_no_cost_growing++;
if (m_factorization->get_status() != LU_status::OK) {
cout << "failing refactor on off_result for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << m_total_iterations << endl;
throw "";
throw exception(sstream() << "failing refactor on off_result for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << m_total_iterations);
}
return false;
}
@ -509,11 +507,11 @@ public:
m_factorization->change_basis(entering, leaving);
init_factorization(m_factorization, m_A, m_basis, m_basis_heading, m_settings, m_non_basic_columns);
if (m_factorization->get_status() != LU_status::OK || A_mult_x_is_off()) {
cout << "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << m_total_iterations << endl;
std::cout << "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << m_total_iterations << std::endl;
restore_x_and_refactor(entering, leaving, tt);
lean_assert(!A_mult_x_is_off());
m_iters_with_no_cost_growing++;
cout << "rolled back after failing of init_factorization()" << endl;
std::cout << "rolled back after failing of init_factorization()" << std::endl;
m_status = UNSTABLE;
return false;
}
@ -571,12 +569,12 @@ public:
restore_x(entering, t);
init_factorization(m_factorization, m_A, m_basis, m_basis_heading, m_settings, m_non_basic_columns);
if (m_factorization->get_status() == LU_status::Degenerated) {
cout << "cannot refactor" << endl;
std::cout << "cannot refactor" << std::endl;
m_status = lp_status::FLOATING_POINT_ERROR;
}
// solve_Ax_eq_b();
if (A_mult_x_is_off()) {
cout << "cannot restore solution" << endl;
std::cout << "cannot restore solution" << std::endl;
m_status = lp_status::FLOATING_POINT_ERROR;
return;
}
@ -584,7 +582,7 @@ public:
void restore_x(unsigned entering, X const & t) {
if (is_zero(t)) return;
cout << "calling restore for entering " << entering << endl;
std::cout << "calling restore for entering " << entering << std::endl;
m_x[entering] -= t;
for (unsigned i : m_index_of_ed) {
m_x[m_basis[i]] = m_copy_of_xB[i];
@ -612,7 +610,7 @@ public:
}
}
void copy_rs_to_xB(vector<X> & rs) {
void copy_rs_to_xB(std::vector<X> & rs) {
unsigned j = m_m;
while (j--) {
m_x[m_basis[j]] = rs[j];
@ -636,26 +634,26 @@ public:
auto it = m_column_names.find(column);
if (it == m_column_names.end()) {
std::string name = T_to_string(column);
return std::string(string("u") + name);
return std::string(std::string("u") + name);
}
return it->second;
}
void copy_right_side(vector<X> & rs) {
void copy_right_side(std::vector<X> & rs) {
unsigned i = m_m;
while (i --) {
rs[i] = m_b[i];
}
}
void add_delta_to_xB(vector<X> & del) {
void add_delta_to_xB(std::vector<X> & del) {
unsigned i = m_m;
while (i--) {
this->m_x[this->m_basis[i]] -= del[i];
}
}
void find_error_in_BxB(vector<X>& rs){
void find_error_in_BxB(std::vector<X>& rs){
unsigned row = m_m;
while (row--) {
auto &rsv = rs[row];
@ -670,9 +668,9 @@ public:
// recalculates the projection of x to B, such that Ax = b, whereab is the right side
void solve_Ax_eq_b() {
vector<X> rs(m_m);
std::vector<X> rs(m_m);
rs_minus_Anx(rs);
vector<X> rrs = rs; // another copy of rs
std::vector<X> rrs = rs; // another copy of rs
m_factorization->solve_By(rs);
copy_rs_to_xB(rs);
find_error_in_BxB(rrs);

170
src/util/lp/lp_dual_core_solver.h
View file

@ -17,7 +17,7 @@ namespace lean {
template <typename T, typename X>
class lp_dual_core_solver:public lp_core_solver_base<T, X> {
public:
vector<bool> & m_can_enter_basis;
std::vector<bool> & m_can_enter_basis;
int m_r; // the row of the leaving column
int m_p; // leaving column; that is m_p = m_basis[m_r]
T m_delta; // the offset of the leaving basis variable
@ -29,21 +29,21 @@ public:
std::set<unsigned> m_breakpoint_set; // it is F in "Progress in the dual simplex method ..."
std::set<unsigned> m_flipped_boxed;
std::set<unsigned> m_tight_set; // it is the set of all breakpoints that become tight when m_q becomes tight
vector<T> m_a_wave;
vector<T> m_betas; // m_betas[i] is approximately a square of the norm of the i-th row of the reverse of B
std::vector<T> m_a_wave;
std::vector<T> m_betas; // m_betas[i] is approximately a square of the norm of the i-th row of the reverse of B
T m_harris_tolerance;
lp_dual_core_solver(static_matrix<T, X> & A,
vector<bool> & can_enter_basis,
vector<X> & b, // the right side vector
vector<X> & x, // the number of elements in x needs to be at least as large as the number of columns in A
vector<unsigned> & basis,
vector<T> & costs,
vector<column_type> & column_type_array,
vector<X> & low_bound_values,
vector<X> & upper_bound_values,
std::vector<bool> & can_enter_basis,
std::vector<X> & b, // the right side std::vector
std::vector<X> & x, // the number of elements in x needs to be at least as large as the number of columns in A
std::vector<unsigned> & basis,
std::vector<T> & costs,
std::vector<column_type> & column_type_array,
std::vector<X> & low_bound_values,
std::vector<X> & upper_bound_values,
lp_settings & settings,
unordered_map<unsigned, string> const & column_names):
std::unordered_map<unsigned, std::string> const & column_names):
lp_core_solver_base<T, X>(A,
b,
basis,
@ -82,11 +82,11 @@ public:
}
void print_nb() {
cout << "this is nb " << endl;
std::cout << "this is nb " << std::endl;
for (auto l : this->m_factorization->m_non_basic_columns) {
cout << l << " ";
std::cout << l << " ";
}
cout << endl;
std::cout << std::endl;
}
void restore_non_basis() {
@ -108,7 +108,7 @@ public:
if (!(this->m_refactor_counter++ >= 200)) {
this->m_factorization->replace_column(leaving, this->m_ed[this->m_factorization->basis_heading(leaving)], this->m_w);
if (this->m_factorization->get_status() != LU_status::OK) {
cout << "failed on replace_column( " << leaving << ", " << this->m_ed[this->m_factorization->basis_heading(leaving)] << ") total_iterations = " << this->m_total_iterations << endl;
std::cout << "failed on replace_column( " << leaving << ", " << this->m_ed[this->m_factorization->basis_heading(leaving)] << ") total_iterations = " << this->m_total_iterations << std::endl;
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
this->m_iters_with_no_cost_growing++;
this->m_status = UNSTABLE;
@ -119,7 +119,7 @@ public:
this->m_factorization->change_basis(entering, leaving);
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
if (this->m_factorization->get_status() != LU_status::OK) {
cout << "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << this->m_total_iterations << endl;
std::cout << "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << this->m_total_iterations << std::endl;
this->m_iters_with_no_cost_growing++;
return false;
}
@ -137,7 +137,7 @@ public:
this->fill_reduced_costs_from_m_y_by_rows();
}
vector<unsigned> & non_basis() { return this->m_factorization->m_non_basic_columns; }
std::vector<unsigned> & non_basis() { return this->m_factorization->m_non_basic_columns; }
void init_betas() {
// todo : look at page 194 of Progress in the dual simplex algorithm for solving large scale LP problems : techniques for a fast and stable implementation
@ -185,10 +185,10 @@ public:
}
// void print_x_and_low_bound(unsigned p) {
// cout << "x l[" << p << "] = " << this->m_x[p] << " " << this->m_low_bound_values[p] << endl;
// std::cout << "x l[" << p << "] = " << this->m_x[p] << " " << this->m_low_bound_values[p] << std::endl;
// }
// void print_x_and_upper_bound(unsigned p) {
// cout << "x u[" << p << "] = " << this->m_x[p] << " " << this->m_upper_bound_values[p] << endl;
// std::cout << "x u[" << p << "] = " << this->m_x[p] << " " << this->m_upper_bound_values[p] << std::endl;
// }
// returns the
@ -199,13 +199,13 @@ public:
case boxed:
if (this->x_below_low_bound(p)) {
T del = get_edge_steepness_for_low_bound(p);
// cout << "case boxed low_bound in pricing" << endl;
// std::cout << "case boxed low_bound in pricing" << std::endl;
// print_x_and_low_bound(p);
return del;
}
if (this->x_above_upper_bound(p)) {
T del = get_edge_steepness_for_upper_bound(p);
// cout << "case boxed at upper_bound in pricing" << endl;
// std::cout << "case boxed at upper_bound in pricing" << std::endl;
// print_x_and_upper_bound(p);
return del;
}
@ -213,7 +213,7 @@ public:
case low_bound:
if (this->x_below_low_bound(p)) {
T del = get_edge_steepness_for_low_bound(p);
// cout << "case low_bound in pricing" << endl;
// std::cout << "case low_bound in pricing" << std::endl;
// print_x_and_low_bound(p);
return del;
}
@ -222,7 +222,7 @@ public:
case upper_bound:
if (this->x_above_upper_bound(p)) {
T del = get_edge_steepness_for_upper_bound(p);
// cout << "case upper_bound in pricing" << endl;
// std::cout << "case upper_bound in pricing" << std::endl;
// print_x_and_upper_bound(p);
return del;
}
@ -232,11 +232,7 @@ public:
lean_assert(numeric_traits<T>::is_zero(this->m_d[p]));
return numeric_traits<T>::zero();
default:
cout << "row = " << i << " p = " << p << endl;
cout << "unexpected column type = " << column_type_to_string(this->m_column_type[p]) << endl;
lean_assert(false);
throw "unexpected column type";
break;
lean_unreachable();
}
}
@ -282,7 +278,7 @@ public:
rows_left = number_of_rows_to_try;
goto loop_start;
// if (this->m_total_iterations % 80 == 0) {
// cout << "m_delta = " << m_delta << endl;
// std::cout << "m_delta = " << m_delta << std::endl;
// }
}
}
@ -312,58 +308,22 @@ public:
if (this->x_above_upper_bound(m_p)) {
return 1;
}
lean_assert(false); // wrong choice of leaving row
throw "wrong choice of leaving row";
lean_unreachable();
case low_bound:
if (this->x_below_low_bound(m_p)) {
return -1;
}
lean_assert(false); // wrong choice of leaving row
throw "wrong choice of leaving row";
lean_unreachable();
case upper_bound:
if (this->x_above_upper_bound(m_p)) {
return 1;
}
lean_assert(false); // wrong choice of leaving row
throw "wrong choice of leaving row";
lean_unreachable();
default:
cout << column_type_to_string(this->m_column_type[m_p]) << endl;
lean_assert(false); // wrong choice of leaving row
throw "wrong choice of leaving row";
break;
lean_unreachable();
}
}
// void try_update_theta_D_and_q_for_boxed(unsigned j) {
// if (this->m_settings.abs_val_is_smaller_than_zero_tolerance(this->m_pivot_row[j])) {
// return;
// }
// T ratio = abs(this->m_d[j] / this->m_pivot_row[j]);
// if (ratio < m_theta_D) {
// m_theta_D = ratio;
// m_q = j;
// m_delta -= abs(bound_span(j) * this->m_pivot_row[j]);
// }
// }
// void try_update_theta_D_and_q(unsigned j) {
// if (this->m_settings.abs_val_is_smaller_than_zero_tolerance(this->m_pivot_row[j])) {
// return;
// cout << " investigate " << endl;
// throw "try_update_theta_D_and_q";
// }
// T ratio = this->m_d[j] / (m_sign_of_alpha_r * this->m_pivot_row[j]);
// if (ratio < m_theta_D) {
// m_theta_D = ratio;
// m_q = j;
// }
// }
bool can_be_breakpoint(unsigned j) {
if (this->pivot_row_element_is_too_small_for_ratio_test(j)) return false;
switch (this->m_column_type[j]) {
@ -416,35 +376,26 @@ public:
if (this->x_above_upper_bound(m_p)) {
return this->m_x[m_p] - this->m_upper_bound_values[m_p];;
}
cout << "incorrect m_p = " << m_p << endl;
lean_assert(false);
throw "incorrect m_p";
lean_unreachable();
case low_bound:
if (this->x_below_low_bound(m_p)) {
return this->m_x[m_p] - this->m_low_bound_values[m_p];
}
cout << "incorrect m_p = " << m_p << endl;
lean_assert(false);
throw "incorrect m_p";
lean_unreachable();
case upper_bound:
if (this->x_above_upper_bound(m_p)) {
return get_edge_steepness_for_upper_bound(m_p);
}
cout << "incorrect m_p = " << m_p << endl;
lean_assert(false);
throw "incorrect m_p";
lean_unreachable();
case fixed:
return this->m_x[m_p] - this->m_upper_bound_values[m_p];;
default:
cout << "unsigned column type " << this->m_column_type[m_p] << endl;
throw "unhandled column type";
break;
lean_unreachable();
}
}
void restore_d() {
cout << "restore_d" << endl;
std::cout << "restore_d" << std::endl;
this->m_d[m_p] = numeric_traits<T>::zero();
for (auto j : non_basis()) {
this->m_d[j] += m_theta_D * this->m_pivot_row[j];
@ -456,8 +407,8 @@ public:
for (auto j : non_basis()) {
T d = this->m_costs[j] - this->m_A.dot_product_with_column(this->m_y, j);
if (numeric_traits<T>::get_double(abs(d - this->m_d[j])) >= 0.001) {
cout << "m_total_iterations = " << this->m_total_iterations << endl;
cout << "d[" << j << "] = " << this->m_d[j] << " but should be " << d << endl;
std::cout << "m_total_iterations = " << this->m_total_iterations << std::endl;
std::cout << "d[" << j << "] = " << this->m_d[j] << " but should be " << d << std::endl;
return false;
}
}
@ -517,11 +468,7 @@ public:
case free_column:
break;
default:
cout << "column = " << j << endl;
cout << "unexpected column type = " << column_type_to_string(this->m_column_type[j]) << endl;
lean_assert(false);
throw "unexpected column type";
break;
lean_unreachable();
}
}
@ -532,7 +479,7 @@ public:
}
void init_beta_precisely(unsigned i) {
vector<T> vec(this->m_m, numeric_traits<T>::zero());
std::vector<T> vec(this->m_m, numeric_traits<T>::zero());
vec[i] = numeric_traits<T>::one();
this->m_factorization->solve_yB(vec);
T beta = numeric_traits<T>::zero();
@ -543,12 +490,12 @@ public:
}
void init_betas_precisely() {
cout << "init beta precisely..." << endl;
std::cout << "init beta precisely..." << std::endl;
unsigned i = this->m_m;
while (i--) {
init_beta_precisely(i);
}
cout << "done" << endl;
std::cout << "done" << std::endl;
}
// step 7 of the algorithm from Progress
@ -569,7 +516,7 @@ public:
}
void revert_to_previous_basis() {
cout << "recovering basis p = " << m_p << " q = " << m_q << endl;
std::cout << "recovering basis p = " << m_p << " q = " << m_q << std::endl;
this->m_factorization->change_basis(m_p, m_q);
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
if (this->m_factorization->get_status() != LU_status::OK) {
@ -590,15 +537,12 @@ public:
bool problem_is_dual_feasible() {
for (unsigned j : non_basis()){
if (!this->column_is_dual_feasible(j)) {
cout << "column " << j << " is not dual feasible" << endl;
cout << "m_d[" << j << "] = " << this->m_d[j] << endl;
cout << "x[" << j << "] = " << this->m_x[j] << endl;
cout << "type = " << column_type_to_string(this->m_column_type[j]) << endl;
cout << "bounds = " << this->m_low_bound_values[j] << "," << this->m_upper_bound_values[j] << endl;
cout << "m_total_iterations = " << this->m_total_iterations << endl;
// cout << "m_betas ";
// print_vector(m_betas);
// column_is_dual_feasible(j); // debug
std::cout << "column " << j << " is not dual feasible" << std::endl;
std::cout << "m_d[" << j << "] = " << this->m_d[j] << std::endl;
std::cout << "x[" << j << "] = " << this->m_x[j] << std::endl;
std::cout << "type = " << column_type_to_string(this->m_column_type[j]) << std::endl;
std::cout << "bounds = " << this->m_low_bound_values[j] << "," << this->m_upper_bound_values[j] << std::endl;
std::cout << "m_total_iterations = " << this->m_total_iterations << std::endl;
return false;
}
}
@ -628,14 +572,14 @@ public:
}
void set_status_to_tentative_dual_unbounded_or_dual_unbounded() {
cout << "cost = " << this->get_cost() << endl;
if (this->m_status == TENTATIVE_DUAL_UNBOUNDED) {
cout << "setting status to DUAL_UNBOUNDED" << endl;
this->m_status = DUAL_UNBOUNDED;
} else {
cout << "setting to TENTATIVE_DUAL_UNBOUNDED" << endl;
this->m_status = TENTATIVE_DUAL_UNBOUNDED;
}
std::cout << "cost = " << this->get_cost() << std::endl;
if (this->m_status == TENTATIVE_DUAL_UNBOUNDED) {
std::cout << "setting status to DUAL_UNBOUNDED" << std::endl;
this->m_status = DUAL_UNBOUNDED;
} else {
std::cout << "setting to TENTATIVE_DUAL_UNBOUNDED" << std::endl;
this->m_status = TENTATIVE_DUAL_UNBOUNDED;
}
}
// it is positive if going from low bound to upper bound and negative if going from upper bound to low bound
@ -817,7 +761,7 @@ public:
this->m_total_iterations = 0;
this->m_iters_with_no_cost_growing = 0;
do {
if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over(string(), this->m_total_iterations)){
if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over(std::string(), this->m_total_iterations)){
this->m_status = lp_status::TIME_EXHAUSTED;
return;
}

61
src/util/lp/lp_dual_simplex.h
View file

@ -4,17 +4,18 @@
Author: Lev Nachmanson
*/
#pragma once
#include <vector>
#include "util/sstream.h"
#include "util/exception.h"
#include "util/lp/lp_solver.h"
#include "util/lp/lp_dual_core_solver.h"
#include <vector>
namespace lean {
template <typename T, typename X>
class lp_dual_simplex: public lp_solver<T, X> {
lp_dual_core_solver<T, X> * m_core_solver = nullptr;
vector<T> m_b_copy;
std::vector<T> m_b_copy;
std::vector<T> m_low_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
std::vector<column_type> m_column_types_of_core_solver;
std::vector<column_type> m_column_types_of_logicals;
@ -31,33 +32,28 @@ public:
case OPTIMAL:
if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) {
this->m_status = FEASIBLE;
cout << "status is FEASIBLE" << endl;
std::cout << "status is FEASIBLE" << std::endl;
} else {
cout << "status is UNBOUNDED" << endl;
std::cout << "status is UNBOUNDED" << std::endl;
this->m_status = UNBOUNDED;
}
break;
case DUAL_UNBOUNDED:
cout << "the status cannot be DUAL_UNBOUNDED since the price is not positive!" << endl;
lean_assert(false);
throw "unexpected status DUAL_UNBOUNDED";
break;
lean_unreachable();
case ITERATIONS_EXHAUSTED:
cout << "status is ITERATIONS_EXHAUSTED" << endl;
std::cout << "status is ITERATIONS_EXHAUSTED" << std::endl;
this->m_status = ITERATIONS_EXHAUSTED;
break;
case TIME_EXHAUSTED:
cout << "status is TIME_EXHAUSTED" << endl;
std::cout << "status is TIME_EXHAUSTED" << std::endl;
this->m_status = TIME_EXHAUSTED;
break;
case FLOATING_POINT_ERROR:
cout << "status is FLOATING_POINT_ERROR" << endl;
std::cout << "status is FLOATING_POINT_ERROR" << std::endl;
this->m_status = FLOATING_POINT_ERROR;
break;
default:
cout << "the status is not expected: " << lp_status_to_string(m_core_solver->get_status()) << endl;
lean_assert(false);
throw "unexpected status";
lean_unreachable();
}
}
@ -75,8 +71,7 @@ public:
m_can_enter_basis[j] = false;
break;
default:
throw 1;
lean_assert(false);
lean_unreachable();
}
}
@ -90,9 +85,7 @@ public:
break;
case fixed:
case upper_bound:
lean_assert(false);
throw "unexpected bound type";
break;
lean_unreachable();
case boxed:
this->m_upper_bounds[j] = ci->get_adjusted_upper_bound() / this->m_column_scale[j];
m_low_bounds[j] = numeric_traits<T>::zero();
@ -104,14 +97,13 @@ public:
m_column_types_of_core_solver[j] = free_column;
break;
default:
lean_assert(false);
throw "unexpected bound type";
lean_unreachable();
}
T cost_was = this->m_costs[j];
this->set_scaled_cost(j);
bool in_basis = m_core_solver->m_factorization->m_basis_heading[j] >= 0;
if (in_basis && cost_was != this->m_costs[j]) {
cout << "cost change in basis" << endl;
std::cout << "cost change in basis" << std::endl;
}
}
@ -154,9 +146,7 @@ public:
this->m_status = FLOATING_POINT_ERROR;
break;
default:
cout << "status of core solver is " << lp_status_to_string(m_core_solver->get_status()) << endl;
lean_assert(false);
throw "unexpected status";
lean_unreachable();
}
this->m_second_stage_iterations = m_core_solver->m_total_iterations;
}
@ -187,11 +177,11 @@ public:
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())) {
cout << "skipping stage 1" << endl;
std::cout << "skipping stage 1" << std::endl;
m_core_solver->set_status(OPTIMAL);
m_core_solver->m_total_iterations = 0;
} else {
cout << "stage 1" << endl;
std::cout << "stage 1" << std::endl;
m_core_solver->solve();
}
decide_on_status_after_stage1();
@ -199,7 +189,7 @@ public:
}
void stage2() {
cout << "starting stage2" << endl;
std::cout << "starting stage2" << std::endl;
unmark_boxed_and_fixed_columns_and_fix_structural_costs();
restore_right_sides();
solve_for_stage2();
@ -230,15 +220,12 @@ public:
T free_bound = T(1e4); // see 4.8
unsigned jj = this->m_core_solver_columns_to_external_columns[j];
lean_assert(this->m_columns.find(jj) != this->m_columns.end());
lean_assert(this->m_columns.find(jj) != this->m_columns.end());
column_info<T> * ci = this->m_columns[jj];
switch (ci->get_column_type()) {
case upper_bound:
cout << j << endl;
cout << column_type_to_string(get_column_type(j)) << endl;
cout << "throwing upper_bound case " << endl;
throw "unexpected bound type"; // we flip columns so this case is impossible
break;
throw exception(sstream() << "unexpected bound type " << j << " "
<< column_type_to_string(get_column_type(j)));
case low_bound: {
m_can_enter_basis[j] = true;
this->set_scaled_cost(j);
@ -259,9 +246,7 @@ public:
this->m_upper_bounds[j] = this->m_low_bounds[j] = numeric_traits<T>::zero(); // is it needed?
break;
default:
lean_assert(false);
throw "unexpected column type";
break;
lean_unreachable();
}
m_column_types_of_core_solver[j] = boxed;
}

65
src/util/lp/lp_primal_core_solver.h
View file

@ -37,15 +37,15 @@ public: // todo : move public lower ; it is for debug only
T m_enter_price_eps;
X m_infeasibility = zero_of_type<X>();
int m_sign_of_entering_delta;
vector<breakpoint<X>> m_breakpoints;
std::vector<breakpoint<X>> m_breakpoints;
binary_heap_priority_queue<X> m_breakpoint_indices_queue;
bool m_using_inf_costs = false;
bool m_recalc_reduced_costs = false;
std::set<unsigned> m_forbidden_enterings;
vector<T> m_beta; // see Swietanowski working vector beta for column norms
std::vector<T> m_beta; // see Swietanowski working vector beta for column norms
T m_epsilon_of_reduced_cost = T(1e-7);
bool m_exit_on_feasible_solution = false;
vector<T> m_costs_backup;
std::vector<T> m_costs_backup;
bool m_current_x_is_feasible;
int choose_entering_column(unsigned number_of_benefitial_columns_to_go_over) { // at this moment m_y = cB * B(-1)
@ -82,9 +82,7 @@ public: // todo : move public lower ; it is for debug only
}
continue;
default:
cout << "unexpected column type " << column_type_to_string(this->m_column_type[j]) << endl;
throw "unexpected column_type";
break;
lean_unreachable();
}
// if we are here then j is a candidate to enter the basis
@ -196,11 +194,11 @@ public: // todo : move public lower ; it is for debug only
} while ( i != initial_i);
restore_harris_eps();
// if (m_j_uht.size() > 1) {
// cout << "under set " << endl;
// cout << "under set " << std::endl;
// for(auto j : m_j_uht)
// print_column(j);
// // }
// cout << "done with under set" << endl;
// cout << "done with under set" << std::endl;
return leaving; // debug
// return leaving_fixed != -1 ? leaving_fixed : leaving;
}
@ -443,8 +441,7 @@ public: // todo : move public lower ; it is for debug only
break;
default:
cout << column_type_to_string(this->m_column_type[j]) << endl;
throw "unexpected type";
lean_unreachable();
}
if (theta < zero_of_type<X>())
theta = zero_of_type<X>();
@ -515,9 +512,9 @@ public: // todo : move public lower ; it is for debug only
}
vector<T> m_low_bound_values_dummy; // needed for the base class only
std::vector<T> m_low_bound_values_dummy; // needed for the base class only
X get_max_bound(vector<X> & b) {
X get_max_bound(std::vector<X> & b) {
X ret = zero_of_type<X>();
for (auto & v : b) {
X a = abs(v);
@ -528,7 +525,7 @@ public: // todo : move public lower ; it is for debug only
// stage1 constructor
lp_primal_core_solver(static_matrix<T, X> & A,
vector<X> & b, // the right side vector
std::vector<X> & b, // the right side vector
std::vector<X> & x, // the number of elements in x needs to be at least as large as the number of columns in A
std::vector<unsigned> & basis,
std::vector<T> & costs,
@ -536,7 +533,7 @@ public: // todo : move public lower ; it is for debug only
std::vector<X> & low_bound_values,
std::vector<X> & upper_bound_values,
lp_settings & settings,
unordered_map<unsigned, string> const & column_names):
std::unordered_map<unsigned, std::string> const & column_names):
lp_core_solver_base<T, X>(A, b,
basis,
x,
@ -553,14 +550,14 @@ public: // todo : move public lower ; it is for debug only
// constructor
lp_primal_core_solver(static_matrix<T, X> & A,
vector<X> & b, // the right side vector
std::vector<X> & b, // the right side vector
std::vector<X> & x, // the number of elements in x needs to be at least as large as the number of columns in A
std::vector<unsigned> & basis,
std::vector<T> & costs,
std::vector<column_type> & column_type_array,
std::vector<X> & upper_bound_values,
lp_settings & settings,
unordered_map<unsigned, string> const & column_names):
std::unordered_map<unsigned, std::string> const & column_names):
lp_core_solver_base<T, X>(A, b,
basis,
x,
@ -593,19 +590,19 @@ public: // todo : move public lower ; it is for debug only
}
for (int j = 0; j < this->m_n; j++) {
if (column_has_low_bound(j) && this->m_x[j] < numeric_traits<T>::zero()) {
cout << "low bound for variable " << j << " does not hold: this->m_x[" << j << "] = " << this->m_x[j] << " is negative " << endl;
std::cout << "low bound for variable " << j << " does not hold: this->m_x[" << j << "] = " << this->m_x[j] << " is negative " << std::endl;
return false;
}
if (column_has_upper_bound(j) && this->m_x[j] > this->m_upper_bound_values[j]) {
cout << "upper bound for " << j << " does not hold: " << this->m_upper_bound_values[j] << ">" << this->m_x[j] << endl;
std::cout << "upper bound for " << j << " does not hold: " << this->m_upper_bound_values[j] << ">" << this->m_x[j] << std::endl;
return false;
}
if (basis_set.find(j) != basis_set.end()) continue;
if (this->m_column_type[j] == low_bound) {
if (numeric_traits<T>::zero() != this->m_x[j]) {
cout << "only low bound is set for " << j << " but low bound value " << numeric_traits<T>::zero() << " is not equal to " << this->m_x[j] << endl;
std::cout << "only low bound is set for " << j << " but low bound value " << numeric_traits<T>::zero() << " is not equal to " << this->m_x[j] << std::endl;
return false;
}
}
@ -754,7 +751,7 @@ public: // todo : move public lower ; it is for debug only
this->restore_x(entering, t * m_sign_of_entering_delta);
m_forbidden_enterings.insert(entering);
this->m_iters_with_no_cost_growing++;
cout << "failing in advance_on_entering_and_leaving for entering == leaving = " << leaving << endl;
std::cout << "failing in advance_on_entering_and_leaving for entering == leaving = " << leaving << std::endl;
return;
}
}
@ -827,7 +824,7 @@ public: // todo : move public lower ; it is for debug only
if (leaving == -1){
if (get_current_x_is_infeasible()) {
if (this->m_status == UNSTABLE) {
cout << "setting status to FLOATING_POINT_ERROR" << endl;
std::cout << "setting status to FLOATING_POINT_ERROR" << std::endl;
this->m_status = FLOATING_POINT_ERROR;
return;
}
@ -861,12 +858,12 @@ public: // todo : move public lower ; it is for debug only
}
void print_column_norms() {
cout << " column norms " << endl;
std::cout << " column norms " << std::endl;
for (unsigned j = 0; j < this->m_n; j++) {
cout << this->m_column_norms[j] << " ";
std::cout << this->m_column_norms[j] << " ";
}
cout << endl;
cout << endl;
std::cout << std::endl;
std::cout << std::endl;
}
void set_current_x_is_feasible() {
@ -987,7 +984,7 @@ public: // todo : move public lower ; it is for debug only
T calculate_column_norm_exactly(unsigned j) {
indexed_vector<T> w(this->m_m);
this->m_A.copy_column_to_vector(j, w);
vector<T> d(this->m_m);
std::vector<T> d(this->m_m);
this->m_factorization->solve_Bd_when_w_is_ready(d, w);
T ret = zero_of_type<T>();
for (auto v : d)
@ -1277,11 +1274,11 @@ public: // todo : move public lower ; it is for debug only
return true;
}
if (this->m_iters_with_no_cost_growing >= this->m_settings.max_number_of_iterations_with_no_improvements) {
cout << "m_iters_with_no_cost_growing = " << this->m_iters_with_no_cost_growing << endl;
std::cout << "m_iters_with_no_cost_growing = " << this->m_iters_with_no_cost_growing << std::endl;
this->m_status = ITERATIONS_EXHAUSTED; return true;
}
if (this->m_total_iterations >= this->m_settings.max_total_number_of_iterations) {
cout << "max_total_number_of_iterations " << this->m_total_iterations << " is reached " << endl;
std::cout << "max_total_number_of_iterations " << this->m_total_iterations << " is reached " << std::endl;
this->m_status = ITERATIONS_EXHAUSTED; return true;
}
return false;
@ -1301,22 +1298,22 @@ public: // todo : move public lower ; it is for debug only
}
void print_column(unsigned j) {
cout << this->column_name(j) << " " << j << " " << column_type_to_string(this->m_column_type[j]) << " x = " << this->m_x[j] << " " << "c = " << this->m_costs[j];;
std::cout << this->column_name(j) << " " << j << " " << column_type_to_string(this->m_column_type[j]) << " x = " << this->m_x[j] << " " << "c = " << this->m_costs[j];;
switch (this->m_column_type[j]) {
case fixed:
case boxed:
cout << "( " << this->m_low_bound_values[j] << " " << this->m_x[j] << " " << this->m_upper_bound_values[j] << ")" << endl;
std::cout << "( " << this->m_low_bound_values[j] << " " << this->m_x[j] << " " << this->m_upper_bound_values[j] << ")" << std::endl;
break;
case upper_bound:
cout << "( _" << this->m_x[j] << " " << this->m_upper_bound_values[j] << ")" << endl;
std::cout << "( _" << this->m_x[j] << " " << this->m_upper_bound_values[j] << ")" << std::endl;
break;
case low_bound:
cout << "( " << this->m_low_bound_values[j] << " " << this->m_x[j] << " " << "_ )" << endl;
std::cout << "( " << this->m_low_bound_values[j] << " " << this->m_x[j] << " " << "_ )" << std::endl;
break;
case free_column:
cout << "( _" << this->m_x[j] << "_)" << endl;
std::cout << "( _" << this->m_x[j] << "_)" << std::endl;
default:
throw "unexpected type";
lean_unreachable();
}
}

67
src/util/lp/lp_primal_simplex.h
View file

@ -4,24 +4,22 @@
Author: Lev Nachmanson
*/
#pragma once
#include "util/lp/lp_primal_core_solver.h"
#include "util/lp/lp_solver.h"
#include <vector>
#include <unordered_map>
#include <string>
#include <algorithm>
#include "util/exception.h"
#include "util/sstream.h"
#include "util/lp/column_info.h"
#include "util/lp/lp_primal_core_solver.h"
#include "util/lp/lp_solver.h"
namespace lean {
using std::vector;
using std::tuple;
template <typename T, typename X>
class lp_primal_simplex: public lp_solver<T, X> {
lp_primal_core_solver<T, X> * m_core_solver = nullptr;
vector<X> m_low_bounds;
std::vector<X> m_low_bounds;
private:
unsigned original_rows() { return this->m_external_rows_to_core_solver_rows.size(); }
@ -34,7 +32,7 @@ private:
}
}
void init_buffer(unsigned k, vector<T> & r) {
void init_buffer(unsigned k, std::vector<T> & r) {
for (unsigned i = 0; i < k; i++) {
r[i] = 0;
}
@ -47,8 +45,7 @@ private:
void refactor() {
m_core_solver->init_lu();
if (m_core_solver->factorization()->get_status() != LU_status::OK) {
cout << "cannot refactor \n";
throw "cannot refactor";
throw exception("cannot refactor");
}
}
@ -60,7 +57,7 @@ private:
}
void stage_two() {
cout << "starting stage 2" << endl;
std::cout << "starting stage 2" << std::endl;
lean_assert(!m_core_solver->A_mult_x_is_off());
int j = this->m_A->column_count() - 1;
unsigned core_solver_cols = this->number_of_core_structurals();
@ -73,7 +70,7 @@ private:
m_core_solver->set_status(lp_status::FEASIBLE);
this->m_second_stage_iterations = m_core_solver->solve();
this->m_status = m_core_solver->get_status();
// cout << "status is " << lp_status_to_string(this->m_status) << endl;
// std::cout << "status is " << lp_status_to_string(this->m_status) << std::endl;
}
public:
lp_primal_simplex() {}
@ -175,11 +172,11 @@ public:
string name_of_core_solver_column(unsigned j) { // j here is the core solver index
std::string name_of_core_solver_column(unsigned j) { // j here is the core solver index
unsigned external_j = this->m_core_solver_columns_to_external_columns[j];
auto t = this->m_columns.find(external_j);
if (t == this->m_columns.end()) {
return string("name_not_found");
return std::string("name_not_found");
}
return t->m_name;
}
@ -274,7 +271,7 @@ public:
}
void solve_with_total_inf() {
cout << "starting solve_with_total_inf()" << endl;
std::cout << "starting solve_with_total_inf()" << std::endl;
int total_vars = this->m_A->column_count() + this->row_count();
m_low_bounds.clear();
m_low_bounds.resize(total_vars, zero_of_type<X>()); // low bounds are shifted ot zero
@ -308,7 +305,7 @@ public:
// void stage_one_of_total_inf() {
// cout << "starting stage_one_of_total_inf()" << endl;
// std::cout << "starting stage_one_of_total_inf()" << std::endl;
// int total_vars = this->m_A->column_count() + this->row_count();
// m_low_bounds.clear();
// m_low_bounds.resize(total_vars, zero_of_type<X>()); // low bounds are shifted ot zero
@ -339,16 +336,16 @@ public:
}
}
bool bounds_hold(unordered_map<string, T> const & solution) {
bool bounds_hold(std::unordered_map<std::string, T> const & solution) {
for (auto it : this->m_columns) {
auto sol_it = solution.find(it.second->get_name());
if (sol_it == solution.end()) {
cout << "cannot find column " << it.first << " in solution " << endl;
std::cout << "cannot find column " << it.first << " in solution " << std::endl;
throw;
}
if (!it.second->bounds_hold(sol_it->second)) {
cout << "bounds do not hold for " << it.second->get_name() << endl;
std::cout << "bounds do not hold for " << it.second->get_name() << std::endl;
it.second->bounds_hold(sol_it->second);
return false;
}
@ -356,37 +353,37 @@ public:
return true;
}
T get_row_value(unsigned i, unordered_map<string, T> const & solution, bool print) {
T get_row_value(unsigned i, std::unordered_map<std::string, T> const & solution, bool print) {
auto it = this->m_A_values.find(i);
if (it == this->m_A_values.end()) {
cout << "cannot find row " << i << endl;
std::cout << "cannot find row " << i << std::endl;
throw "get_row_value";
}
T ret = numeric_traits<T>::zero();
for (auto & pair : it->second) {
auto cit = this->m_columns.find(pair.first);
if (cit == this->m_columns.end()){
cout << "cannot find column " << pair.first << endl;
std::cout << "cannot find column " << pair.first << std::endl;
}
column_info<T> * ci = cit->second;
auto sol_it = solution.find(ci->get_name());
if (sol_it == solution.end()) {
cout << "cannot find in the solution column " << ci->get_name() << endl;
std::cout << "cannot find in the solution column " << ci->get_name() << std::endl;
}
T column_val = sol_it->second;
if (print) {
cout << pair.second << "(" << ci->get_name() << "=" << column_val << ") ";
std::cout << pair.second << "(" << ci->get_name() << "=" << column_val << ") ";
}
ret += pair.second * column_val;
}
if (print) {
cout << " = " << ret << endl;
std::cout << " = " << ret << std::endl;
}
return ret;
}
bool row_constraint_holds(unsigned i, unordered_map<string, T> const & solution, bool print) {
bool row_constraint_holds(unsigned i, std::unordered_map<std::string, T> const & solution, bool print) {
T row_val = get_row_value(i, solution, print);
auto & constraint = this->m_constraints[i];
T rs = constraint.m_rs;
@ -394,7 +391,7 @@ public:
case Equal:
if (fabs(numeric_traits<T>::get_double(row_val - rs)) > 0.00001) {
if (print) {
cout << "should be = " << rs << endl;
std::cout << "should be = " << rs << std::endl;
}
return false;
}
@ -402,7 +399,7 @@ public:
case Greater_or_equal:
if (numeric_traits<T>::get_double(row_val - rs) < -0.00001) {
if (print) {
cout << "should be >= " << rs << endl;
std::cout << "should be >= " << rs << std::endl;
}
return false;
}
@ -411,17 +408,17 @@ public:
case Less_or_equal:
if (numeric_traits<T>::get_double(row_val - rs) > 0.00001) {
if (print) {
cout << "should be <= " << rs << endl;
std::cout << "should be <= " << rs << std::endl;
}
return false;
}
return true;;
}
cout << "throw in row_constraint_holds " << endl;
std::cout << "throw in row_constraint_holds " << std::endl;
throw "wrong case";
}
bool row_constraints_hold(unordered_map<string, T> const & solution) {
bool 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, false)) {
row_constraint_holds(it.first, solution, true);
@ -432,21 +429,19 @@ public:
}
T * get_array_from_map(unordered_map<string, T> const & solution) {
T * get_array_from_map(std::unordered_map<std::string, T> const & solution) {
T * t = new T[solution.size()];
for (auto it : solution) {
auto g = this->m_names_to_columns.find(it.first);
if (g == this->m_names_to_columns.end()) {
string s = string("cannot find name ") + " "+ it.first;
cout << "throw in get_array_from_map" << endl;
throw s;
throw exception(sstream() << "cannot find name " << it.first);
}
t[g->second] = it.second;
}
return t;
}
bool solution_is_feasible(unordered_map<string, T> const & solution) {
bool solution_is_feasible(std::unordered_map<std::string, T> const & solution) {
return bounds_hold(solution) && row_constraints_hold(solution);
}

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

@ -23,8 +23,6 @@ enum lp_relation {
Greater_or_equal
};
template <typename T, typename X>
struct lp_constraint {
X m_rs; // right side of the constraint
@ -50,26 +48,26 @@ protected:
public:
unsigned m_total_iterations;
static_matrix<T, X>* m_A = nullptr; // this is the matrix of constraints
vector<T> m_b; // the right side vector
std::vector<T> m_b; // the right side vector
unsigned m_first_stage_iterations = 0;
unsigned m_second_stage_iterations = 0;
unordered_map<unsigned, lp_constraint<T, X>> m_constraints;
unordered_map<unsigned, column_info<T>*> m_columns;
unordered_map<unsigned, unordered_map<unsigned, T> > m_A_values;
unordered_map<string, unsigned> m_names_to_columns; // don't have to use it
unordered_map<unsigned, unsigned> m_external_rows_to_core_solver_rows;
unordered_map<unsigned, unsigned> m_core_solver_rows_to_external_rows;
unordered_map<unsigned, unsigned> m_external_columns_to_core_solver_columns;
unordered_map<unsigned, unsigned> m_core_solver_columns_to_external_columns;
vector<T> m_column_scale;
unordered_map<unsigned, std::string> m_name_map;
std::unordered_map<unsigned, lp_constraint<T, X>> m_constraints;
std::unordered_map<unsigned, column_info<T>*> m_columns;
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_external_columns_to_core_solver_columns;
std::unordered_map<unsigned, unsigned> m_core_solver_columns_to_external_columns;
std::vector<T> m_column_scale;
std::unordered_map<unsigned, std::string> m_name_map;
unsigned m_artificials = 0;
unsigned m_slacks = 0;
vector<column_type> m_column_types;
vector<T> m_costs;
vector<T> m_x;
vector<T> m_upper_bounds;
vector<unsigned> m_basis;
std::vector<column_type> m_column_types;
std::vector<T> m_costs;
std::vector<T> m_x;
std::vector<T> m_upper_bounds;
std::vector<unsigned> m_basis;
lp_status m_status = lp_status::UNKNOWN;
@ -94,7 +92,7 @@ public:
// returns the current cost
virtual T get_current_cost() const = 0;
// do not have to call it
void give_symbolic_name_to_column(string name, unsigned column) {
void give_symbolic_name_to_column(std::string name, unsigned column) {
auto it = m_columns.find(column);
column_info<T> *ci;
if (it == m_columns.end()){
@ -111,7 +109,7 @@ public:
T get_column_value_by_name(std::string name) const {
auto it = m_names_to_columns.find(name);
if (it == m_names_to_columns.end()) {
cout << "throwing in get_column_value_by_name" << endl;
std::cout << "throwing in get_column_value_by_name" << std::endl;
throw "get_column_value_by_name";
}
return get_column_value(it -> second);
@ -211,29 +209,29 @@ protected:
void print_rows_scale_stats() {
cout << "rows max" << endl;
for (unsigned i = 0; i < m_A->row_count(); i++) {
std::cout << "rows max" << std::endl;
for (unsigned i = 0; i < m_A->row_count(); i++) {
print_row_scale_stats(i);
}
cout << endl;
std::cout << std::endl;
}
void print_columns_scale_stats() {
cout << "columns max" << endl;
std::cout << "columns max" << std::endl;
for (unsigned i = 0; i < m_A->column_count(); i++) {
print_column_scale_stats(i);
}
cout << endl;
std::cout << std::endl;
}
void print_row_scale_stats(unsigned i) {
cout << "(" << std::min(m_A->get_min_abs_in_row(i), abs(m_b[i])) << " ";
cout << std::max(m_A->get_max_abs_in_row(i), abs(m_b[i])) << ")";
std::cout << "(" << std::min(m_A->get_min_abs_in_row(i), abs(m_b[i])) << " ";
std::cout << std::max(m_A->get_max_abs_in_row(i), abs(m_b[i])) << ")";
}
void print_column_scale_stats(unsigned j) {
cout << "(" << m_A->get_min_abs_in_row(j) << " ";
cout << m_A->get_max_abs_in_column(j) << ")";
std::cout << "(" << m_A->get_min_abs_in_row(j) << " ";
std::cout << m_A->get_max_abs_in_column(j) << ")";
}
void print_scale_stats() {
@ -241,7 +239,7 @@ protected:
print_columns_scale_stats();
}
void get_max_abs_in_row(unordered_map<unsigned, T> & row_map) {
void 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);
@ -252,16 +250,16 @@ protected:
return ret;
}
void pin_vars_down_on_row(unordered_map<unsigned, T> & row) {
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(unordered_map<unsigned, T> & row) {
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(unordered_map<unsigned, T> & row, T sign ) {
unordered_map<unsigned, T> pinned;
void pin_vars_on_row_with_sign(std::unordered_map<unsigned, T> & row, T sign ) {
std::unordered_map<unsigned, T> pinned;
for (auto t : row) {
unsigned j = t.first;
column_info<T> * ci = m_columns[j];
@ -276,7 +274,7 @@ protected:
}
}
bool get_minimal_row_value(unordered_map<unsigned, T> & row, T & low_bound) {
bool get_minimal_row_value(std::unordered_map<unsigned, T> & row, T & low_bound) {
low_bound = numeric_traits<T>::zero();
for (auto & t : row) {
T a = t.second;
@ -298,7 +296,7 @@ protected:
return true;
}
bool get_maximal_row_value(unordered_map<unsigned, T> & row, T & low_bound) {
bool get_maximal_row_value(std::unordered_map<unsigned, T> & row, T & low_bound) {
low_bound = numeric_traits<T>::zero();
for (auto & t : row) {
T a = t.second;
@ -320,7 +318,7 @@ protected:
return true;
}
bool row_is_zero(unordered_map<unsigned, T> & row) {
bool row_is_zero(std::unordered_map<unsigned, T> & row) {
for (auto & t : row) {
if (!is_zero(t.second))
return false;
@ -328,7 +326,7 @@ protected:
return true;
}
bool row_e_is_obsolete(unordered_map<unsigned, T> & row, unsigned row_index) {
bool 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))
@ -369,7 +367,7 @@ protected:
return false;
}
int row_ge_is_obsolete(unordered_map<unsigned, T> & row, unsigned row_index) {
int 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>())
@ -394,7 +392,7 @@ protected:
return false;
}
bool row_le_is_obsolete(unordered_map<unsigned, T> & row, unsigned row_index) {
bool row_le_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index) {
T low_bound;
T rs = m_constraints[row_index].m_rs;
if (row_is_zero(row)) {
@ -425,7 +423,7 @@ protected:
// 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(unordered_map<unsigned, T> & row, unsigned row_index ) {
bool 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:
@ -448,9 +446,9 @@ protected:
}
}
void remove_fixed_or_zero_columns_from_row(unsigned i, unordered_map<unsigned, T> & row) {
void remove_fixed_or_zero_columns_from_row(unsigned i, std::unordered_map<unsigned, T> & row) {
auto & constraint = m_constraints[i];
vector<unsigned> removed;
std::vector<unsigned> removed;
for (auto & col : row) {
unsigned j = col.first;
lean_assert(m_columns.find(j) != m_columns.end());
@ -472,7 +470,7 @@ protected:
}
unsigned try_to_remove_some_rows() {
vector<unsigned> rows_to_delete;
std::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);
@ -498,9 +496,9 @@ protected:
n += d;
if (n == 1)
cout << "deleted one row" << endl;
std::cout << "deleted one row" << std::endl;
else if (n)
cout << "deleted " << n << " rows" << endl;
std::cout << "deleted " << n << " rows" << std::endl;
}
void map_external_rows_to_core_solver_rows() {
@ -517,7 +515,7 @@ protected:
for (auto & row : m_A_values) {
for (auto & col : row.second) {
if (col.second == numeric_traits<T>::zero() || m_columns[col.first]->is_fixed()) {
cout << "found fixed column " << endl;
std::cout << "found fixed column " << std::endl;
throw "map_external_columns_to_core_solver_columns";
}
unsigned j = col.first;
@ -675,10 +673,10 @@ protected:
this->m_costs[j] = cost * this->m_column_scale[j];
}
void print_statistics_on_A() {
cout << "extended A[" << this->m_A->row_count() << "," << this->m_A->column_count() << "]" << endl;
std::cout << "extended A[" << this->m_A->row_count() << "," << this->m_A->column_count() << "]" << std::endl;
// for (unsigned i = 0; i < this->m_A->row_count(); i++) {
// if (this->m_A->number_of_non_zeroes_in_row(i) <= 2 ) {
// cout << "m_p[" << i << "] = " << this->m_A->number_of_non_zeroes_in_row(i) << endl;
// std::cout << "m_p[" << i << "] = " << this->m_A->number_of_non_zeroes_in_row(i) << std::endl;
// }
// }
}

108
src/util/lp/lu.h
View file

@ -27,10 +27,10 @@ std::string T_to_string(const T & t); // forward definition
#ifdef LEAN_DEBUG
template <typename T, typename X> // print the nr x nc submatrix at the top left corner
void print_submatrix(sparse_matrix<T, X> & m, unsigned mr, unsigned nc) {
vector<vector<string>> A;
vector<unsigned> widths;
std::vector<std::vector<std::string>> A;
std::vector<unsigned> widths;
for (unsigned i = 0; i < m.row_count() && i < mr ; i++) {
A.push_back(vector<string>());
A.push_back(std::vector<std::string>());
for (unsigned j = 0; j < m.column_count() && j < nc; j++) {
A[i].push_back(T_to_string(static_cast<T>(m(i, j))));
}
@ -45,11 +45,11 @@ void print_submatrix(sparse_matrix<T, X> & m, unsigned mr, unsigned nc) {
template<typename T, typename X>
void print_matrix(static_matrix<T, X> &m) {
vector<vector<string>> A;
vector<unsigned> widths;
std::vector<std::vector<std::string>> A;
std::vector<unsigned> widths;
std::set<pair<unsigned, unsigned>> domain = m.get_domain();
for (unsigned i = 0; i < m.row_count(); i++) {
A.push_back(vector<string>());
A.push_back(std::vector<std::string>());
for (unsigned j = 0; j < m.column_count(); j++) {
A[i].push_back(T_to_string(static_cast<T>(m(i, j))));
}
@ -64,10 +64,10 @@ void print_matrix(static_matrix<T, X> &m) {
template <typename T, typename X>
void print_matrix(sparse_matrix<T, X>& m) {
vector<vector<string>> A;
vector<unsigned> widths;
std::vector<std::vector<std::string>> A;
std::vector<unsigned> widths;
for (unsigned i = 0; i < m.row_count(); i++) {
A.push_back(vector<string>());
A.push_back(std::vector<std::string>());
for (unsigned j = 0; j < m.column_count(); j++) {
A[i].push_back(T_to_string(static_cast<T>(m(i, j))));
}
@ -135,11 +135,11 @@ public:
unsigned row_count() const { return m_m; } // not defined }
unsigned column_count() const { return m_m; } // not defined }
#endif
void apply_from_left(vector<X> & w, lp_settings &) {
void apply_from_left(std::vector<X> & w, lp_settings &) {
w[m_i] /= m_val;
}
void apply_from_right(vector<T> & w) {
void apply_from_right(std::vector<T> & w) {
w[m_i] /= m_val;
}
@ -177,7 +177,7 @@ public:
// the fields
unsigned m_dim;
static_matrix<T, X> const &m_A;
vector<unsigned>& m_basis;
std::vector<unsigned>& m_basis;
permutation_matrix<T, X> m_Q;
permutation_matrix<T, X> m_R;
sparse_matrix<T, X> m_U;
@ -185,20 +185,20 @@ public:
// m_tail is composed of tail_matrices:
// one_off_diagonal_matrix, and transposition_matrix
vector<tail_matrix<T, X> *> m_tail;
std::vector<tail_matrix<T, X> *> m_tail;
lp_settings & m_settings;
vector<int> & m_basis_heading;
std::vector<int> & m_basis_heading;
bool m_failure = false;
vector<unsigned> & m_non_basic_columns;
std::vector<unsigned> & m_non_basic_columns;
indexed_vector<T> m_row_eta_work_vector;
// constructor
// if A is an m by n matrix then basis has length m and values in [0,n); the values are all different
// they represent the set of m columns
lu(static_matrix<T, X> const & A,
vector<unsigned>& basis,
vector<int> & basis_heading,
std::vector<unsigned>& basis,
std::vector<int> & basis_heading,
lp_settings & settings,
vector<unsigned> & non_basic_columns):
std::vector<unsigned> & non_basic_columns):
m_dim(A.row_count()),
m_A(A),
m_basis(basis),
@ -215,9 +215,9 @@ public:
create_initial_factorization();
if (get_status() != LU_status::OK) {
if (get_status() == LU_status::Degenerated) {
cout << "lu status is Degenerated" << endl;
std::cout << "lu status is Degenerated" << std::endl;
} else {
cout << "lu status is " <<static_cast<int>(get_status()) << endl;
std::cout << "lu status is " <<static_cast<int>(get_status()) << std::endl;
}
return;
}
@ -240,7 +240,7 @@ public:
return - m_basis_heading[j] - 1;
}
void solve_By(vector<X> & y) {
void solve_By(std::vector<X> & y) {
init_vector_y(y);
solve_By_when_y_is_ready(y);
}
@ -253,7 +253,7 @@ public:
}
template <typename L>
void solve_By_when_y_is_ready(vector<L> & y) {
void solve_By_when_y_is_ready(std::vector<L> & y) {
m_U.double_solve_U_y(y);
m_R.apply_reverse_from_left(y); // see 24.3 from Chvatal
if (precise<X>()) return;
@ -268,42 +268,42 @@ public:
void print_indexed_vector(indexed_vector<T> & w, std::ofstream & f) {
f << "vector_start" << endl;
f << "vector_start" << std::endl;
for (unsigned j : w.m_index) {
f << j << " " << w[j] << endl;
f << j << " " << w[j] << std::endl;
}
f << "vector_end" << endl;
f << "vector_end" << std::endl;
}
void print_basis(std::ofstream & f) {
f << "basis_start" << endl;
f << "basis_start" << std::endl;
for (unsigned j : m_basis)
f << j << endl;
f << "basis_end" << endl;
f << j << std::endl;
f << "basis_end" << std::endl;
}
void print_matrix_compact(std::ofstream & f) {
f << "matrix_start" << endl;
f << "nrows " << m_A.row_count() << endl;
f << "ncolumns " << m_A.column_count() << endl;
f << "matrix_start" << std::endl;
f << "nrows " << m_A.row_count() << std::endl;
f << "ncolumns " << m_A.column_count() << std::endl;
for (unsigned i = 0; i < m_A.row_count(); i++) {
auto & row = m_A.m_rows[i];
f << "row " << i << endl;
f << "row " << i << std::endl;
for (auto & t : row.m_cells) {
f << "column " << t.m_j << " value " << t.m_value << endl;
f << "column " << t.m_j << " value " << t.m_value << std::endl;
}
f << "row_end" << endl;
f << "row_end" << std::endl;
}
f << "matrix_end" << endl;
f << "matrix_end" << std::endl;
}
void print(indexed_vector<T> & w) {
string dump_file_name("/tmp/lu");
std::string dump_file_name("/tmp/lu");
remove(dump_file_name.c_str());
std::ofstream f(dump_file_name);
if (!f.is_open()) {
cout << "cannot open file " << dump_file_name << endl;
std::cout << "cannot open file " << dump_file_name << std::endl;
return;
}
cout << "writing lu dump to " << dump_file_name << endl;
std::cout << "writing lu dump to " << dump_file_name << std::endl;
print_matrix_compact(f);
print_basis(f);
print_indexed_vector(w, f);
@ -319,7 +319,7 @@ public:
return m_basis_heading[i] < 0;
}
void solve_yB_internal(vector<T>& y) {
void solve_yB_internal(std::vector<T>& y) {
// first solve yU = cb*R(-1)
m_R.apply_reverse_from_right(y); // got y = cb*R(-1)
m_U.solve_y_U(y); // got y*U=cb*R(-1)
@ -332,13 +332,13 @@ public:
}
}
void add_delta_to_solution(vector<T>& yc, vector<T>& y){
void add_delta_to_solution(std::vector<T>& yc, std::vector<T>& y){
unsigned i = y.size();
while (i--)
y[i]+=yc[i];
}
void find_error_of_yB(vector<T>& yc, const vector<T>& y) {
void find_error_of_yB(std::vector<T>& yc, const std::vector<T>& y) {
unsigned i = m_dim;
while (i--) {
yc[i] -= m_A.dot_product_with_column(y, m_basis[i]);
@ -348,8 +348,8 @@ public:
// solves y*B = y
// y is the input
void solve_yB(vector<T> & y) {
vector<T> yc(y); // copy y aside
void solve_yB(std::vector<T> & y) {
std::vector<T> yc(y); // copy y aside
solve_yB_internal(y);
find_error_of_yB(yc, y);
solve_yB_internal(yc);
@ -400,7 +400,7 @@ public:
return m_A(i, m_basis[j]);
}
void init_vector_y(vector<X> & y) {
void init_vector_y(std::vector<X> & y) {
apply_lp_lists_to_y(y);
m_Q.apply_reverse_from_left(y);
}
@ -426,7 +426,7 @@ public:
}
}
void apply_lp_lists_to_y(vector<X>& y) {
void apply_lp_lists_to_y(std::vector<X>& y) {
for (unsigned i = 0; i < m_tail.size(); i++) {
m_tail[i]->apply_from_left(y, m_settings);
}
@ -559,14 +559,14 @@ public:
unsigned pi, pj;
m_U.get_pivot_for_column(pi, pj, T(m_settings.c_partial_pivoting), j);
if (pi == -1) {
cout << "cannot find the pivot for column " << j << endl;
std::cout << "cannot find the pivot for column " << j << std::endl;
m_failure = true;
return;
}
swap_columns(j, pj);
swap_rows(j, pi);
if (!pivot_the_row(j)) {
cout << "pivot_the_row(" << j << ") failed" << endl;
std::cout << "pivot_the_row(" << j << ") failed" << std::endl;
m_failure = true;
}
}
@ -681,7 +681,7 @@ public:
return;
}
j++;
// cout << "switching to dense factoring for " << j << endl;
// std::cout << "switching to dense factoring for " << j << endl;
m_dense_LU = new square_dense_submatrix<T, X>(&m_U, j);
for (; j < m_dim; j++) {
pivot_in_dense_mode(j);
@ -715,7 +715,7 @@ public:
}
void scan_last_row_to_work_vector(unsigned lowest_row_of_the_bump) {
vector<indexed_value<T>> & last_row_vec = m_U.get_row_values(m_U.adjust_row(lowest_row_of_the_bump));
std::vector<indexed_value<T>> & last_row_vec = m_U.get_row_values(m_U.adjust_row(lowest_row_of_the_bump));
for (auto & iv : last_row_vec) {
if (is_zero(iv.m_value)) continue;
lean_assert(!m_settings.abs_val_is_smaller_than_drop_tolerance(iv.m_value));
@ -739,7 +739,7 @@ public:
T v = m_row_eta_work_vector[j];
if (numeric_traits<T>::is_zero(v)) continue; // this column does not contribute to the solution
unsigned aj = m_U.adjust_row(j);
vector<indexed_value<T>> & row = m_U.get_row_values(aj);
std::vector<indexed_value<T>> & row = m_U.get_row_values(aj);
for (auto & iv : row) {
unsigned col = m_U.adjust_column_inverse(iv.m_index);
lean_assert(col >= j || numeric_traits<T>::is_zero(iv.m_value));
@ -780,9 +780,9 @@ public:
!is_zero(pivot_elem_for_checking) &&
#endif
!m_settings.abs_val_is_smaller_than_pivot_tolerance((m_row_eta_work_vector[lowest_row_of_the_bump] - pivot_elem_for_checking) / denom)) {
// cout << "m_row_eta_work_vector[" << lowest_row_of_the_bump << "] = " << T_to_string(m_row_eta_work_vector[lowest_row_of_the_bump]) << ", but pivot = " << T_to_string(pivot_elem_for_checking) << endl;
// std::cout << "m_row_eta_work_vector[" << lowest_row_of_the_bump << "] = " << T_to_string(m_row_eta_work_vector[lowest_row_of_the_bump]) << ", but pivot = " << T_to_string(pivot_elem_for_checking) << endl;
set_status(LU_status::Degenerated);
// cout << "diagonal element is off" << endl;
// std::cout << "diagonal element is off" << endl;
return nullptr;
}
#ifdef LEAN_DEBUG
@ -861,13 +861,13 @@ public:
}
}; // end of lu
template <typename T, typename X>
void init_factorization(lu<T, X>* & factorization, static_matrix<T, X> & m_A, std::vector<unsigned> & m_basis, vector<int> & m_basis_heading, lp_settings &m_settings, vector<unsigned> & non_basic_columns) {
void init_factorization(lu<T, X>* & factorization, static_matrix<T, X> & m_A, std::vector<unsigned> & m_basis, std::vector<int> & m_basis_heading, lp_settings &m_settings, std::vector<unsigned> & non_basic_columns) {
if (factorization != nullptr) {
delete factorization;
}
factorization = new lu<T, X>(m_A, m_basis, m_basis_heading, m_settings, non_basic_columns);
if (factorization->get_status() != LU_status::OK) {
cout << "failing in init_factorization" << endl;
std::cout << "failing in init_factorization" << std::endl;
return;
}
}

6
src/util/lp/matrix_domain.h
View file

@ -36,15 +36,13 @@ template<typename S, typename T> struct hash<pair<S, T>> {
}
namespace lean {
using std::unordered_map;
template <typename T>
class matrix_domain {
vector<unordered_map<unsigned, void *>> m_domain;
std::vector<std::unordered_map<unsigned, void *>> m_domain;
public:
matrix_domain(unsigned rows) {
while (rows--) {
unordered_map<unsigned, void *> t;
std::unordered_map<unsigned, void *> t;
m_domain.push_back(t);
}
}

2167
src/util/lp/permutation_matrix.h
View file

File diff suppressed because it is too large Load diff

6
src/util/lp/row_eta_matrix.h
View file

@ -54,7 +54,7 @@ public:
return m_row_vector.m_data[m_row];
}
void apply_from_left(vector<X> & w, lp_settings &
void apply_from_left(std::vector<X> & w, lp_settings &
#ifdef LEAN_DEBUG
// settings
#endif
@ -114,7 +114,7 @@ public:
m_row_vector.push_back(row_index, val);
}
void apply_from_right(vector<T> & w) {
void apply_from_right(std::vector<T> & w) {
T w_row = w[m_row];
if (numeric_traits<T>::is_zero(w_row)) return;
#ifdef LEAN_DEBUG
@ -170,7 +170,7 @@ public:
#endif
m_row = p.apply_reverse(m_row);
// copy aside the column indices
vector<unsigned> columns;
std::vector<unsigned> columns;
for (auto it = sparse_vector_iterator<T>(m_row_vector); !it.done(); it.move()) {
columns.push_back(it.index());
}

473
src/util/lp/scaler.h
View file

@ -6,264 +6,263 @@
*/
#pragma once
#include "util/numerics/double.h"
#include <vector>
#include <math.h>
#include <algorithm>
#include <stdio.h> /* printf, fopen */
#include <stdlib.h> /* exit, EXIT_FAILURE */
#include "util/numerics/double.h"
namespace lean {
// for scaling an LP
template <typename T, typename X>
class scaler {
vector<X> & m_b; // right side
static_matrix<T, X> &m_A; // the constraint matrix
const T & m_scaling_minimum;
const T & m_scaling_maximum;
vector<T>& m_column_scale;
lp_settings & m_settings;
public:
// constructor
scaler(vector<X> & b, static_matrix<T, X> &A, const T & scaling_minimum, const T & scaling_maximum, vector<T> & column_scale,
lp_settings & settings):
m_b(b),
m_A(A),
m_scaling_minimum(scaling_minimum),
m_scaling_maximum(scaling_maximum),
m_column_scale(column_scale),
m_settings(settings) {
lean_assert(m_column_scale.size() == 0);
m_column_scale.resize(m_A.column_count(), numeric_traits<T>::one());
}
// for scaling an LP
template <typename T, typename X>
class scaler {
std::vector<X> & m_b; // right side
static_matrix<T, X> &m_A; // the constraint matrix
const T & m_scaling_minimum;
const T & m_scaling_maximum;
std::vector<T>& m_column_scale;
lp_settings & m_settings;
public:
// constructor
scaler(std::vector<X> & b, static_matrix<T, X> &A, const T & scaling_minimum, const T & scaling_maximum, std::vector<T> & column_scale,
lp_settings & settings):
m_b(b),
m_A(A),
m_scaling_minimum(scaling_minimum),
m_scaling_maximum(scaling_maximum),
m_column_scale(column_scale),
m_settings(settings) {
lean_assert(m_column_scale.size() == 0);
m_column_scale.resize(m_A.column_count(), numeric_traits<T>::one());
}
T right_side_balance() {
T ret = zero_of_type<T>();
unsigned i = m_A.row_count();
while (i--) {
T rs = abs(convert_struct<T, X>::convert(m_b[i]));
if (!is_zero<T>(rs)) {
numeric_traits<T>::log(rs);
ret += rs * rs;
}
}
return ret;
}
T get_balance() {
return m_A.get_balance();
}
T A_min() const {
T min = zero_of_type<T>();
for (unsigned i = 0; i < m_A.row_count(); i++) {
T t = m_A.get_min_abs_in_row(i);
min = i == 0 ? t : std::min(t, min);
}
return min;
}
T A_max() const {
T max = zero_of_type<T>();
for (unsigned i = 0; i < m_A.row_count(); i++) {
T t = m_A.get_max_abs_in_row(i);
max = i == 0? t : std::max(t, max);
}
return max;
}
T get_A_ratio() const {
T min = A_min();
T max = A_max();
T ratio = max / min;
return ratio;
}
T get_max_ratio_on_rows() const {
T ret = T(1);
unsigned i = m_A.row_count();
while (i--) {
T t = m_A.get_max_abs_in_row(i)/m_A.get_min_abs_in_row(i);
if (t > ret)
ret = t;
}
return ret;
}
T get_max_ratio_on_columns() const {
T ret = T(1);
unsigned i = m_A.column_count();
while (i--) {
T t = m_A.get_max_abs_in_column(i)/m_A.get_min_abs_in_column(i);
if (t > ret)
ret = t;
}
return ret;
}
void scale_rows_with_geometric_mean() {
unsigned i = m_A.row_count();
while (i--) {
T max = m_A.get_max_abs_in_row(i);
T min = m_A.get_min_abs_in_row(i);
lean_assert(max > zero_of_type<T>() && min > zero_of_type<T>());
T gm = T(sqrt(numeric_traits<T>::get_double(max*min)));
m_A.divide_row_by_constant(i, gm);
m_b[i] /= gm;
T right_side_balance() {
T ret = zero_of_type<T>();
unsigned i = m_A.row_count();
while (i--) {
T rs = abs(convert_struct<T, X>::convert(m_b[i]));
if (!is_zero<T>(rs)) {
numeric_traits<T>::log(rs);
ret += rs * rs;
}
}
return ret;
}
void scale_columns_with_geometric_mean() {
unsigned i = m_A.column_count();
while (i--) {
T max = m_A.get_max_abs_in_column(i);
T min = m_A.get_min_abs_in_column(i);
T gm = T(1)/T(sqrt(numeric_traits<T>::get_double(max*min)));
m_A.scale_column(i, gm);
m_column_scale[i]*=gm;
}
T get_balance() {
return m_A.get_balance();
}
T A_min() const {
T min = zero_of_type<T>();
for (unsigned i = 0; i < m_A.row_count(); i++) {
T t = m_A.get_min_abs_in_row(i);
min = i == 0 ? t : std::min(t, min);
}
return min;
}
void scale_once_for_ratio() {
T max_ratio_on_rows = get_max_ratio_on_rows();
T max_ratio_on_columns = get_max_ratio_on_columns();
bool scale_rows_first = max_ratio_on_rows > max_ratio_on_columns;
// if max_ratio_on_columns is the largerst then the rows are in worser shape then columns
if (scale_rows_first) {
scale_rows_with_geometric_mean();
scale_columns_with_geometric_mean();
T A_max() const {
T max = zero_of_type<T>();
for (unsigned i = 0; i < m_A.row_count(); i++) {
T t = m_A.get_max_abs_in_row(i);
max = i == 0? t : std::max(t, max);
}
return max;
}
T get_A_ratio() const {
T min = A_min();
T max = A_max();
T ratio = max / min;
return ratio;
}
T get_max_ratio_on_rows() const {
T ret = T(1);
unsigned i = m_A.row_count();
while (i--) {
T t = m_A.get_max_abs_in_row(i)/m_A.get_min_abs_in_row(i);
if (t > ret)
ret = t;
}
return ret;
}
T get_max_ratio_on_columns() const {
T ret = T(1);
unsigned i = m_A.column_count();
while (i--) {
T t = m_A.get_max_abs_in_column(i)/m_A.get_min_abs_in_column(i);
if (t > ret)
ret = t;
}
return ret;
}
void scale_rows_with_geometric_mean() {
unsigned i = m_A.row_count();
while (i--) {
T max = m_A.get_max_abs_in_row(i);
T min = m_A.get_min_abs_in_row(i);
lean_assert(max > zero_of_type<T>() && min > zero_of_type<T>());
T gm = T(sqrt(numeric_traits<T>::get_double(max*min)));
m_A.divide_row_by_constant(i, gm);
m_b[i] /= gm;
}
}
void scale_columns_with_geometric_mean() {
unsigned i = m_A.column_count();
while (i--) {
T max = m_A.get_max_abs_in_column(i);
T min = m_A.get_min_abs_in_column(i);
T gm = T(1)/T(sqrt(numeric_traits<T>::get_double(max*min)));
m_A.scale_column(i, gm);
m_column_scale[i]*=gm;
}
}
void scale_once_for_ratio() {
T max_ratio_on_rows = get_max_ratio_on_rows();
T max_ratio_on_columns = get_max_ratio_on_columns();
bool scale_rows_first = max_ratio_on_rows > max_ratio_on_columns;
// if max_ratio_on_columns is the largerst then the rows are in worser shape then columns
if (scale_rows_first) {
scale_rows_with_geometric_mean();
scale_columns_with_geometric_mean();
} else {
scale_columns_with_geometric_mean();
scale_rows_with_geometric_mean();
}
}
bool scale_with_ratio() {
T ratio = get_A_ratio();
// The ratio is greater than or equal to one. We would like to diminish it and bring it as close to 1 as possible
unsigned reps = 20;
do {
scale_once_for_ratio();
T new_r = get_A_ratio();
if (new_r >= T(0.9) * ratio)
break;
} while (reps--);
bring_rows_and_columns_maximums_to_one();
return true;
}
void bring_row_maximums_to_one() {
unsigned i = m_A.row_count();
while (i--) {
T t = m_A.get_max_abs_in_row(i);
m_A.divide_row_by_constant(i, t);
m_b[i] /= t;
}
}
void bring_column_maximums_to_one() {
unsigned i = m_A.column_count();
while (i--) {
T t = T(1) / m_A.get_max_abs_in_column(i);
m_A.scale_column(i, t);
m_column_scale[i] *= t;
}
}
void bring_rows_and_columns_maximums_to_one() {
if (get_max_ratio_on_rows() > get_max_ratio_on_columns()) {
bring_row_maximums_to_one();
bring_column_maximums_to_one();
} else {
bring_column_maximums_to_one();
bring_row_maximums_to_one();
}
}
bool scale_with_log_balance() {
T balance = get_balance();
T balance_before_scaling = balance;
// todo : analyze the scale order : rows-columns, or columns-rows. Iterate if needed
for (int i = 0; i < 10; i++) {
scale_rows();
scale_columns();
T nb = get_balance();
if (nb < T(0.9) * balance) {
balance = nb;
} else {
scale_columns_with_geometric_mean();
scale_rows_with_geometric_mean();
balance = nb;
break;
}
}
return balance <= balance_before_scaling;
}
// Returns true if and only if the scaling was successful.
// It is the caller responsibility to restore the matrix
bool scale() {
if (numeric_traits<T>::precise()) return true;
if (m_settings.scale_with_ratio)
return scale_with_ratio();
return scale_with_log_balance();
}
bool scale_with_ratio() {
T ratio = get_A_ratio();
// The ratio is greater than or equal to one. We would like to diminish it and bring it as close to 1 as possible
unsigned reps = 20;
void scale_rows() {
for (unsigned i = 0; i < m_A.row_count(); i++)
scale_row(i);
}
void scale_row(unsigned i) {
T row_max = std::max(m_A.get_max_abs_in_row(i), abs(convert_struct<T, X>::convert(m_b[i])));
T alpha = numeric_traits<T>::one();
if (numeric_traits<T>::is_zero(row_max)) {
return;
}
if (numeric_traits<T>::get_double(row_max) < m_scaling_minimum) {
do {
scale_once_for_ratio();
T new_r = get_A_ratio();
if (new_r >= T(0.9) * ratio)
break;
} while (reps--);
alpha *= 2;
row_max *= 2;
} while (numeric_traits<T>::get_double(row_max) < m_scaling_minimum);
m_A.scale_row(i, alpha);
m_b[i] *= alpha;
} else if (numeric_traits<T>::get_double(row_max) > m_scaling_maximum) {
do {
alpha /= 2;
row_max /= 2;
} while (numeric_traits<T>::get_double(row_max) > m_scaling_maximum);
m_A.scale_row(i, alpha);
m_b[i] *= alpha;
}
}
bring_rows_and_columns_maximums_to_one();
return true;
void scale_column(unsigned i){
T column_max = m_A.get_max_abs_in_column(i);
T alpha = numeric_traits<T>::one();
if (numeric_traits<T>::is_zero(column_max)){
return; // the column has zeros only
}
void bring_row_maximums_to_one() {
unsigned i = m_A.row_count();
while (i--) {
T t = m_A.get_max_abs_in_row(i);
m_A.divide_row_by_constant(i, t);
m_b[i] /= t;
}
if (numeric_traits<T>::get_double(column_max) < m_scaling_minimum) {
do {
alpha *= 2;
column_max *= 2;
} while (numeric_traits<T>::get_double(column_max) < m_scaling_minimum);
} else if (numeric_traits<T>::get_double(column_max) > m_scaling_maximum) {
do {
alpha /= 2;
column_max /= 2;
} while (numeric_traits<T>::get_double(column_max) > m_scaling_maximum);
} else {
return;
}
m_A.scale_column(i, alpha);
m_column_scale[i] = alpha;
}
void bring_column_maximums_to_one() {
unsigned i = m_A.column_count();
while (i--) {
T t = T(1) / m_A.get_max_abs_in_column(i);
m_A.scale_column(i, t);
m_column_scale[i] *= t;
}
void scale_columns() {
for (unsigned i = 0; i < m_A.column_count(); i++) {
scale_column(i);
}
void bring_rows_and_columns_maximums_to_one() {
if (get_max_ratio_on_rows() > get_max_ratio_on_columns()) {
bring_row_maximums_to_one();
bring_column_maximums_to_one();
} else {
bring_column_maximums_to_one();
bring_row_maximums_to_one();
}
}
bool scale_with_log_balance() {
T balance = get_balance();
T balance_before_scaling = balance;
// todo : analyze the scale order : rows-columns, or columns-rows. Iterate if needed
for (int i = 0; i < 10; i++) {
scale_rows();
scale_columns();
T nb = get_balance();
if (nb < T(0.9) * balance) {
balance = nb;
} else {
balance = nb;
break;
}
}
return balance <= balance_before_scaling;
}
// Returns true if and only if the scaling was successful.
// It is the caller responsibility to restore the matrix
bool scale() {
if (numeric_traits<T>::precise()) return true;
if (m_settings.scale_with_ratio)
return scale_with_ratio();
return scale_with_log_balance();
}
void scale_rows() {
for (unsigned i = 0; i < m_A.row_count(); i++)
scale_row(i);
}
void scale_row(unsigned i) {
T row_max = std::max(m_A.get_max_abs_in_row(i), abs(convert_struct<T, X>::convert(m_b[i])));
T alpha = numeric_traits<T>::one();
if (numeric_traits<T>::is_zero(row_max)) {
return;
}
if (numeric_traits<T>::get_double(row_max) < m_scaling_minimum) {
do {
alpha *= 2;
row_max *= 2;
} while (numeric_traits<T>::get_double(row_max) < m_scaling_minimum);
m_A.scale_row(i, alpha);
m_b[i] *= alpha;
} else if (numeric_traits<T>::get_double(row_max) > m_scaling_maximum) {
do {
alpha /= 2;
row_max /= 2;
} while (numeric_traits<T>::get_double(row_max) > m_scaling_maximum);
m_A.scale_row(i, alpha);
m_b[i] *= alpha;
}
}
void scale_column(unsigned i){
T column_max = m_A.get_max_abs_in_column(i);
T alpha = numeric_traits<T>::one();
if (numeric_traits<T>::is_zero(column_max)){
return; // the column has zeros only
}
if (numeric_traits<T>::get_double(column_max) < m_scaling_minimum) {
do {
alpha *= 2;
column_max *= 2;
} while (numeric_traits<T>::get_double(column_max) < m_scaling_minimum);
} else if (numeric_traits<T>::get_double(column_max) > m_scaling_maximum) {
do {
alpha /= 2;
column_max /= 2;
} while (numeric_traits<T>::get_double(column_max) > m_scaling_maximum);
} else {
return;
}
m_A.scale_column(i, alpha);
m_column_scale[i] = alpha;
}
void scale_columns() {
for (unsigned i = 0; i < m_A.column_count(); i++) {
scale_column(i);
}
}
};
}
};
}

146
src/util/lp/sparse_matrix.h
View file

@ -23,10 +23,7 @@
#include "util/lp/eta_matrix.h"
#include "util/lp/binary_heap_upair_queue.h"
namespace lean {
using std::vector;
using std::cout;
// it is a square matrix
// it is a square matrix
template <typename T, typename X>
class sparse_matrix
#ifdef LEAN_DEBUG
@ -35,7 +32,7 @@ class sparse_matrix
{
struct col_header {
unsigned m_shortened_markovitz = 0;
vector<indexed_value<T>> m_values; // the actual column values
std::vector<indexed_value<T>> m_values; // the actual column values
col_header() {}
@ -50,8 +47,8 @@ class sparse_matrix
binary_heap_upair_queue<unsigned> m_pivot_queue;
public:
vector<vector<indexed_value<T>>> m_rows;
vector<col_header> m_columns;
std::vector<std::vector<indexed_value<T>>> m_rows;
std::vector<col_header> m_columns;
permutation_matrix<T, X> m_row_permutation;
permutation_matrix<T, X> m_column_permutation;
indexed_vector<T> m_work_pivot_vector;
@ -80,13 +77,13 @@ public:
}
void copy_column_from_static_matrix(unsigned col, static_matrix<T, X> const &A, unsigned col_index_in_the_new_matrix) {
vector<column_cell<T>> const & A_col_vector = A.m_columns[col];
std::vector<column_cell<T>> const & A_col_vector = A.m_columns[col];
unsigned size = A_col_vector.size();
vector<indexed_value<T>> & new_column_vector = m_columns[col_index_in_the_new_matrix].m_values;
std::vector<indexed_value<T>> & new_column_vector = m_columns[col_index_in_the_new_matrix].m_values;
for (unsigned l = 0; l < size; l++) {
column_cell<T> const & col_cell = A_col_vector[l];
unsigned col_offset = new_column_vector.size();
vector<indexed_value<T>> & row_vector = m_rows[col_cell.m_i];
std::vector<indexed_value<T>> & row_vector = m_rows[col_cell.m_i];
unsigned row_offset = row_vector.size();
new_column_vector.push_back(indexed_value<T>(col_cell.m_value, col_cell.m_i, row_offset));
row_vector.push_back(indexed_value<T>(col_cell.m_value, col_index_in_the_new_matrix, col_offset));
@ -134,7 +131,7 @@ public:
};
void set_with_no_adjusting_for_row(unsigned row, unsigned col, T val) { // should not be used in efficient code
vector<indexed_value<T>> & row_vec = m_rows[row];
std::vector<indexed_value<T>> & row_vec = m_rows[row];
for (auto & iv : row_vec) {
if (iv.m_index == col) {
iv.set_value(val);
@ -146,7 +143,7 @@ public:
}
void set_with_no_adjusting_for_col(unsigned row, unsigned col, T val) { // should not be used in efficient code
vector<indexed_value<T>> & col_vec = m_columns[col].m_values;
std::vector<indexed_value<T>> & col_vec = m_columns[col].m_values;
for (auto & iv : col_vec) {
if (iv.m_index == row) {
iv.set_value(val);
@ -203,19 +200,19 @@ public:
// end of access for debugging helpers
vector<indexed_value<T>> & get_row_values(unsigned row) {
std::vector<indexed_value<T>> & get_row_values(unsigned row) {
return m_rows[row];
}
vector<indexed_value<T>> const & get_row_values(unsigned row) const {
std::vector<indexed_value<T>> const & get_row_values(unsigned row) const {
return m_rows[row];
}
vector<indexed_value<T>> & get_column_values(unsigned col) {
std::vector<indexed_value<T>> & get_column_values(unsigned col) {
return m_columns[col].m_values;
}
vector<indexed_value<T>> const & get_column_values(unsigned col) const {
std::vector<indexed_value<T>> const & get_column_values(unsigned col) const {
return m_columns[col].m_values;
}
@ -240,7 +237,7 @@ public:
void init_row_headers() {
for (unsigned l = 0; l < m_row_permutation.size(); l++) {
m_rows.push_back(vector<indexed_value<T>>());
m_rows.push_back(std::vector<indexed_value<T>>());
}
}
@ -284,7 +281,7 @@ public:
}
void remove_element(vector<indexed_value<T>> & row_vals, unsigned row_offset, vector<indexed_value<T>> & column_vals, unsigned column_offset) {
void remove_element(std::vector<indexed_value<T>> & row_vals, unsigned row_offset, std::vector<indexed_value<T>> & column_vals, unsigned column_offset) {
if (column_offset != column_vals.size() - 1) {
auto & column_iv = column_vals[column_offset] = column_vals.back(); // copy from the tail
column_iv_other(column_iv).m_other = column_offset;
@ -301,13 +298,13 @@ public:
row_vals.pop_back();
}
void remove_element(vector<indexed_value<T>> & row_chunk, indexed_value<T> & row_el_iv) {
void remove_element(std::vector<indexed_value<T>> & row_chunk, indexed_value<T> & row_el_iv) {
auto & column_chunk = get_column_values(row_el_iv.m_index);
indexed_value<T> & col_el_iv = column_chunk[row_el_iv.m_other];
remove_element(row_chunk, col_el_iv.m_other, column_chunk, row_el_iv.m_other);
}
void put_max_index_to_0(vector<indexed_value<T>> & row_vals, unsigned max_index) {
void put_max_index_to_0(std::vector<indexed_value<T>> & row_vals, unsigned max_index) {
if (max_index == 0) return;
indexed_value<T> * max_iv = & row_vals[max_index];
indexed_value<T> * start_iv = & row_vals[0];
@ -325,7 +322,7 @@ public:
set_max_in_row(m_rows[row]);
}
void set_max_in_row(vector<indexed_value<T>> & row_vals) {
void set_max_in_row(std::vector<indexed_value<T>> & row_vals) {
if (row_vals.size() == 0)
return;
T max_val = numeric_traits<T>::zero();
@ -535,11 +532,11 @@ public:
}
void swap_columns(unsigned a, unsigned b) {
// cout << "swaapoiiin" << endl;
// dense_matrix<T, X> d(*this);
// cout << "swaapoiiin" << std::endl;
// dense_matrix<T, X> d(*this);
m_column_permutation.transpose_from_left(a, b);
// d.swap_columns(a, b);
// lean_assert(*this == d);
// d.swap_columns(a, b);
// lean_assert(*this == d);
}
void swap_rows(unsigned a, unsigned b) {
@ -554,7 +551,7 @@ public:
for (auto & iv : m_rows[i]) {
T v = iv.m_value / t;
if (settings.abs_val_is_smaller_than_drop_tolerance(v)){
v = numeric_traits<T>::zero();
v = numeric_traits<T>::zero();
}
m_columns[iv.m_index].m_values[iv.m_other].set_value(iv.m_value = v);
}
@ -685,7 +682,7 @@ public:
}
template <typename L>
void find_error_in_solution_U_y(vector<L>& y_orig, vector<L> & y) {
void find_error_in_solution_U_y(std::vector<L>& y_orig, std::vector<L> & y) {
unsigned i = dimension();
while (i--) {
y_orig[i] -= dot_product_with_row(i, y);
@ -693,7 +690,7 @@ public:
}
template <typename L>
void add_delta_to_solution(const vector<L>& del, vector<L> & y) {
void add_delta_to_solution(const std::vector<L>& del, std::vector<L> & y) {
unsigned i = dimension();
while (i--) {
y[i] += del[i];
@ -701,8 +698,8 @@ public:
}
template <typename L>
void double_solve_U_y(vector<L>& y){
vector<L> y_orig(y); // copy y aside
void double_solve_U_y(std::vector<L>& y){
std::vector<L> y_orig(y); // copy y aside
solve_U_y(y);
find_error_in_solution_U_y(y_orig, y);
// y_orig contains the error now
@ -744,7 +741,7 @@ public:
virtual void set_number_of_columns(unsigned /*n*/) { }
#endif
template <typename L>
L dot_product_with_row (unsigned row, const vector<L> & y) const {
L dot_product_with_row (unsigned row, const std::vector<L> & y) const {
L ret = zero_of_type<L>();
auto & mc = get_row_values(adjust_row(row));
for (auto & c : mc) {
@ -781,7 +778,7 @@ public:
return false;
}
void remove_element_that_is_not_in_w(vector<indexed_value<T>> & column_vals, indexed_value<T> & col_el_iv) {
void remove_element_that_is_not_in_w(std::vector<indexed_value<T>> & column_vals, indexed_value<T> & col_el_iv) {
auto & row_chunk = m_rows[col_el_iv.m_index];
indexed_value<T> & row_el_iv = row_chunk[col_el_iv.m_other];
unsigned index_in_row = col_el_iv.m_other;
@ -865,11 +862,11 @@ public:
unsigned pivot_score(unsigned i, unsigned j) {
// It goes like this (rnz-1)(cnz-1) is the Markovitz number, that is the max number of
// new non zeroes we can obtain after the pivoting.
// In addition we will get another cnz - 1 elements in the eta matrix created for this pivot,
// which gives rnz(cnz-1). For example, is 0 for a column singleton, but not for
// a row singleton ( which is not a column singleton).
// It goes like this (rnz-1)(cnz-1) is the Markovitz number, that is the max number of
// new non zeroes we can obtain after the pivoting.
// In addition we will get another cnz - 1 elements in the eta matrix created for this pivot,
// which gives rnz(cnz-1). For example, is 0 for a column singleton, but not for
// a row singleton ( which is not a column singleton).
auto col_header = m_columns[j];
@ -906,7 +903,7 @@ public:
enqueue_domain_into_pivot_queue();
}
void recover_pivot_queue(vector<upair> & rejected_pivots) {
void recover_pivot_queue(std::vector<upair> & rejected_pivots) {
for (auto p : rejected_pivots) {
m_pivot_queue.enqueue(p.first, p.second, pivot_score(p.first, p.second));
}
@ -925,8 +922,6 @@ public:
return abs(iv.m_value) * c_partial_pivoting < max ? 1: 0;
}
return 2; // the element became zero but it still sits in the active pivots?
cout << "column " << j << " does not belong to row " << i << endl;
throw "cannot be here";
}
bool remove_row_from_active_pivots_and_shorten_columns(unsigned row) {
@ -994,7 +989,7 @@ public:
unsigned cnz = cols.size();
for (auto & iv : cols) {
if (adjust_row_inverse(iv.m_index) < k)
cnz--;
cnz--;
}
lean_assert(cnz > 0);
return m_rows[i].m_values.size() * (cnz - 1);
@ -1030,9 +1025,9 @@ public:
return true;
}
void print_active_matrix(unsigned k) {
cout << "active matrix for k = " << k << endl;
std::cout << "active matrix for k = " << k << std::endl;
if (k >= dimension()) {
cout << "empty" << endl;
std::cout << "empty" << std::endl;
return;
}
unsigned dim = dimension() - k;
@ -1057,7 +1052,7 @@ public:
unsigned arow = adjust_row(i);
for (auto & iv : m_rows[arow].m_values) {
lean_assert(pivot_score_without_shortened_counters(arow, iv.m_index, k + 1) ==
m_pivot_queue.get_priority(arow, iv.m_index));
m_pivot_queue.get_priority(arow, iv.m_index));
}
return true;
}
@ -1071,7 +1066,7 @@ public:
bool get_pivot_for_column(unsigned &i, unsigned &j, T const & c_partial_pivoting, unsigned k) {
vector<upair> pivots_candidates_that_are_too_small;
std::vector<upair> pivots_candidates_that_are_too_small;
while (!m_pivot_queue.is_empty()) {
m_pivot_queue.dequeue(i, j);
unsigned i_inv = adjust_row_inverse(i);
@ -1098,26 +1093,26 @@ public:
recover_pivot_queue(pivots_candidates_that_are_too_small);
return false;
/*
unsigned m = dimension();
unsigned markovitz_number_min = m * m;
unsigned iterations = 0;
for (unsigned k = 2; k <= m; k++){
if (markovitz_number_min < (k-1) * (k-1))
return true;
unsigned m = dimension();
unsigned markovitz_number_min = m * m;
unsigned iterations = 0;
for (unsigned k = 2; k <= m; k++){
if (markovitz_number_min < (k-1) * (k-1))
return true;
if (get_pivot_in_k_short_columns(k, markovitz_number_min, i, j, c_partial_pivoting, iterations)
||
get_pivot_in_k_short_rows(k, markovitz_number_min, i, j, c_partial_pivoting, iterations)) {
return true;
}
if (iterations > search_depth * 2 && markovitz_number_min < m * m)
return true;
}
return markovitz_number_min < m * m;
if (get_pivot_in_k_short_columns(k, markovitz_number_min, i, j, c_partial_pivoting, iterations)
||
get_pivot_in_k_short_rows(k, markovitz_number_min, i, j, c_partial_pivoting, iterations)) {
return true;
}
if (iterations > search_depth * 2 && markovitz_number_min < m * m)
return true;
}
return markovitz_number_min < m * m;
*/
}
bool elem_is_too_small(vector<indexed_value<T>> & row_chunk, indexed_value<T> & iv, T const & c_partial_pivoting) {
bool elem_is_too_small(std::vector<indexed_value<T>> & row_chunk, indexed_value<T> & iv, T const & c_partial_pivoting) {
if (&iv == &row_chunk[0]) {
return false; // the max element is at the head
}
@ -1137,7 +1132,7 @@ public:
bool shorten_columns_by_pivot_row(unsigned i, unsigned pivot_column) {
vector<indexed_value<T>> & row_chunk = get_row_values(i);
std::vector<indexed_value<T>> & row_chunk = get_row_values(i);
for (indexed_value<T> & iv : row_chunk) {
unsigned j = iv.m_index;
@ -1163,7 +1158,7 @@ public:
}
void fill_eta_matrix(unsigned j, eta_matrix<T, X> ** eta) {
vector<indexed_value<T>> & col_chunk = get_column_values(adjust_column(j));
std::vector<indexed_value<T>> & col_chunk = get_column_values(adjust_column(j));
bool is_unit = true;
for (auto & iv : col_chunk) {
unsigned i = adjust_row_inverse(iv.m_index);
@ -1203,7 +1198,7 @@ public:
bool is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings & settings) const {
for (unsigned i = 0; i < dimension(); i++) {
vector<indexed_value<T>> const & row_chunk = get_row_values(i);
std::vector<indexed_value<T>> const & row_chunk = get_row_values(i);
lean_assert(row_chunk.size());
T const & max = abs(row_chunk[0].m_value);
unsigned ai = adjust_row_inverse(i);
@ -1214,7 +1209,7 @@ public:
return false;
if (aj == ai) {
if (iv.m_value != 1) {
cout << "value at diagonal = " << iv.m_value << endl;
std::cout << "value at diagonal = " << iv.m_value << std::endl;
return false;
}
}
@ -1245,18 +1240,18 @@ public:
for (indexed_value<T> & column_iv : mc) {
indexed_value<T> & row_iv = column_iv_other(column_iv);
if (row_iv.m_index != col) {
cout << "m_other in row does not belong to column " << col << ", but to column " << row_iv.m_index << endl;
std::cout << "m_other in row does not belong to column " << col << ", but to column " << row_iv.m_index << std::endl;
lean_assert(false);
}
if (& row_iv_other(row_iv) != &column_iv) {
cout << "row and col do not point to each other" << endl;
std::cout << "row and col do not point to each other" << std::endl;
lean_assert(false);
}
if (row_iv.m_value != column_iv.m_value) {
cout << "the data from col " << col << " for row " << column_iv.m_index << " is different in the column " << endl;
cout << "in the col it is " << column_iv.m_value << ", but in the row it is " << row_iv.m_value << endl;
std::cout << "the data from col " << col << " for row " << column_iv.m_index << " is different in the column " << std::endl;
std::cout << "in the col it is " << column_iv.m_value << ", but in the row it is " << row_iv.m_value << std::endl;
lean_assert(false);
}
}
@ -1268,18 +1263,18 @@ public:
indexed_value<T> & column_iv = row_iv_other(row_iv);
if (column_iv.m_index != row) {
cout << "col_iv does not point to correct row " << row << " but to " << column_iv.m_index << endl;
std::cout << "col_iv does not point to correct row " << row << " but to " << column_iv.m_index << std::endl;
lean_assert(false);
}
if (& row_iv != & column_iv_other(column_iv)) {
cout << "row and col do not point to each other" << endl;
std::cout << "row and col do not point to each other" << std::endl;
lean_assert(false);
}
if (row_iv.m_value != column_iv.m_value) {
cout << "the data from col " << column_iv.m_index << " for row " << row << " is different in the column " << endl;
cout << "in the col it is " << column_iv.m_value << ", but in the row it is " << row_iv.m_value << endl;
std::cout << "the data from col " << column_iv.m_index << " for row " << row << " is different in the column " << std::endl;
std::cout << "in the col it is " << column_iv.m_value << ", but in the row it is " << row_iv.m_value << std::endl;
lean_assert(false);
}
}
@ -1297,10 +1292,11 @@ public:
}
}
void map_domain_to_vector(unordered_map<pair<unsigned, unsigned>, unsigned> & domain, vector<pair<unsigned, unsigned>> & vec) {
void map_domain_to_vector(std::unordered_map<pair<unsigned, unsigned>, unsigned> & domain,
std::vector<pair<unsigned, unsigned>> & vec) {
lean_assert(domain.size() == 0 && vec.size() == 0);
for (unsigned i = 0; i < dimension(); i++) {
vector<indexed_value<T>> & row = get_row_values(i);
std::vector<indexed_value<T>> & row = get_row_values(i);
for (auto & iv : row) {
unsigned j = iv.m_index;
unsigned nz = vec.size();

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

@ -8,6 +8,7 @@
#pragma once
#include <vector>
#include "util/debug.h"
#include "util/pair.h"
#include "util/numerics/numeric_traits.h"
#include "util/numerics/xnumeral.h"
#include "util/numerics/mpq.h"
@ -18,8 +19,6 @@
#include "util/numerics/mpfp.h"
#include "util/lp/lp_settings.h"
namespace lean {
using std::vector;
using std::pair;
template <typename T>
void zero_vector(T * t, unsigned size) {
while (size-- > 0) { // it can be made faster by copying big chunks
@ -37,7 +36,7 @@ class sparse_vector_iterator; // forward definition
template <typename T>
class sparse_vector {
public:
vector<pair<unsigned, T>> m_data;
std::vector<pair<unsigned, T>> m_data;
void push_back(unsigned index, T val) {
m_data.emplace_back(index, val);
}
@ -66,7 +65,7 @@ public:
template <typename T>
class sparse_vector_iterator {
typedef typename vector< pair<unsigned, T> >::iterator p_it;
typedef typename std::vector<pair<unsigned, T>>::iterator p_it;
p_it m_it;
p_it m_end;
public:

24
src/util/lp/square_dense_submatrix.h
View file

@ -23,8 +23,6 @@
#include "util/lp/eta_matrix.h"
#include "util/lp/binary_heap_upair_queue.h"
namespace lean {
using std::vector;
using std::cout;
template <typename T, typename X>
class square_dense_submatrix : public tail_matrix<T, X> {
// the submatrix uses the permutations of the parent matrix to access the elements
@ -45,7 +43,7 @@ class square_dense_submatrix : public tail_matrix<T, X> {
public: // debug
unsigned m_index_start;
unsigned m_dim;
vector<T> m_v;
std::vector<T> m_v;
sparse_matrix<T, X> * m_parent = nullptr;
permutation_matrix<T, X> m_row_permutation;
public:
@ -230,7 +228,7 @@ public:
}
}
template <typename L>
L row_by_vector_product(unsigned i, const vector<L> & v) {
L row_by_vector_product(unsigned i, const std::vector<L> & v) {
lean_assert(i >= m_index_start);
unsigned row_in_subm = i - m_index_start;
@ -242,7 +240,7 @@ public:
}
template <typename L>
L column_by_vector_product(unsigned j, const vector<L> & v) {
L column_by_vector_product(unsigned j, const std::vector<L> & v) {
lean_assert(j >= m_index_start);
unsigned offset = j - m_index_start;
@ -276,7 +274,7 @@ public:
#endif // use indexed vector here
#ifndef DO_NOT_USE_INDEX
vector<L> t(m_parent->dimension(), zero_of_type<L>());
std::vector<L> t(m_parent->dimension(), zero_of_type<L>());
for (auto k : w.m_index) {
unsigned j = adjust_column(k); // k-th element will contribute only to column j
if (j < m_index_start) {
@ -301,7 +299,7 @@ public:
}
}
#else
vector<L> t(m_parent->dimension());
std::vector<L> t(m_parent->dimension());
for (unsigned i = 0; i < m_index_start; i++) {
t[adjust_row_inverse(i)] = w[adjust_column_inverse(i)];
}
@ -327,7 +325,7 @@ public:
}
template <typename L>
void apply_from_left_to_vector(vector<L> & w) {
void apply_from_left_to_vector(std::vector<L> & w) {
// lp_settings & settings) {
// dense_matrix<T, L> deb(*this);
// vector<L> deb_w(w);
@ -338,7 +336,7 @@ public:
// // print_vector(w.m_data);
// // cout << "deb_w" << endl;
// // print_vector(deb_w);
vector<L> t(m_parent->dimension());
std::vector<L> t(m_parent->dimension());
for (unsigned i = 0; i < m_index_start; i++) {
t[adjust_row_inverse(i)] = w[adjust_column_inverse(i)];
}
@ -383,17 +381,17 @@ public:
void apply_from_right(indexed_vector<T> & w) {
lean_assert(false); // not implemented
}
void apply_from_left(vector<X> & w, lp_settings & /*settings*/) {
void apply_from_left(std::vector<X> & w, lp_settings & /*settings*/) {
apply_from_left_to_vector(w);// , settings);
}
void apply_from_right(vector<T> & w) {
void apply_from_right(std::vector<T> & w) {
#ifdef LEAN_DEBUG
// dense_matrix<T, X> deb(*this);
// vector<T> deb_w(w);
// deb.apply_from_right(deb_w);
#endif
vector<T> t(w.size());
std::vector<T> t(w.size());
for (unsigned j = 0; j < m_index_start; j++) {
t[adjust_column_inverse(j)] = w[adjust_row_inverse(j)];
}
@ -404,7 +402,7 @@ public:
lean_assert(w.size() == t.size());
w = t;
#ifdef LEAN_DEBUG
// lean_assert(vectors_are_equal<T>(deb_w, w));
// lean_assert(vector_are_equal<T>(deb_w, w));
#endif
}

19
src/util/lp/static_matrix.h
View file

@ -61,12 +61,12 @@ class static_matrix
std::set<pair<unsigned, unsigned>> m_domain;
#endif
public:
typedef vector<row_cell<T>> row_strip;
typedef vector<column_cell<T>> column_strip;
vector<T> m_work_pivot_vector;
typedef std::vector<row_cell<T>> row_strip;
typedef std::vector<column_cell<T>> column_strip;
std::vector<T> m_work_pivot_vector;
vector<row_strip> m_rows;
vector<column_strip> m_columns;
std::vector<row_strip> m_rows;
std::vector<column_strip> m_columns;
// starting inner classes
class ref {
static_matrix & m_matrix;
@ -283,7 +283,7 @@ public:
v.set_value(it.m_value, it.m_i);
}
}
void copy_column_to_vector (unsigned j, vector<T> & v) const {
void copy_column_to_vector (unsigned j, std::vector<T> & v) const {
v.resize(row_count(), numeric_traits<T>::zero());
for (auto & it : m_columns[j]) {
if (!is_zero(it.m_value))
@ -352,7 +352,7 @@ public:
#ifdef LEAN_DEBUG
void check_consistency() {
unordered_map<std::pair<unsigned, unsigned>, T> by_rows;
std::unordered_map<std::pair<unsigned, unsigned>, T> by_rows;
for (int i = 0; i < m_rows.size(); i++){
for (auto & t : m_rows[i]) {
pair<unsigned, unsigned> p(i, t.m_j);
@ -360,7 +360,7 @@ public:
by_rows[p] = t.get_val();
}
}
unordered_map<pair<unsigned, unsigned>, T> by_cols;
std::unordered_map<pair<unsigned, unsigned>, T> by_cols;
for (int i = 0; i < m_columns.size(); i++){
for (auto & t : m_columns[i]) {
pair<unsigned, unsigned> p(t.m_i, i);
@ -373,7 +373,8 @@ public:
for (auto & t : by_rows) {
auto ic = by_cols.find(t.first);
if (ic == by_cols.end()){
cout << "rows have pair (" << t.first.first <<"," << t.first.second << "), but columns don't " << endl;
std::cout << "rows have pair (" << t.first.first <<"," << t.first.second
<< "), but columns don't " << std::endl;
}
lean_assert(ic != by_cols.end());
lean_assert(t.second == ic->second);