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

dev(lp): integrate with z3

Browse source 
Signed-off-by: Lev Nachmanson <levnach@microsoft.com>
This commit is contained in:
Lev Nachmanson 2016-05-09 15:24:47 -07:00 committed by Leonardo de Moura
parent b58eac5013
commit d823349d2e
52 changed files with 418 additions and 385 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -87,7 +87,7 @@ void test_matrix(sparse_matrix<T, X> & a) {
void tst1() {
std::cout << "testing the minimial matrix with 1 row and 1 column" << std::endl;
sparse_matrix<float, float> m0(1);
sparse_matrix<double, double> m0(1);
m0.set(0, 0, 1);
// print_matrix(m0);
m0.set(0, 0, 0);
@ -95,17 +95,17 @@ void tst1() {
test_matrix(m0);
unsigned rows = 2;
sparse_matrix<float, float> m(rows);
sparse_matrix<double, double> m(rows);
std::cout << "setting m(0,1)=" << std::endl;
m.set(0, 1, 11);
m.set(0, 0, 12);
// print_matrix<float>(m);
// print_matrix<double>(m);
test_matrix(m);
sparse_matrix<float, float> m1(2);
sparse_matrix<double, double> m1(2);
m1.set(0, 0, 2);
m1.set(1, 0, 3);
// print_matrix(m1);
@ -118,7 +118,7 @@ void tst1() {
std::cout << "printing zero matrix 3 by 1" << std::endl;
sparse_matrix<float, float> m2(3);
sparse_matrix<double, double> m2(3);
// print_matrix(m2);
m2.set(0, 0, 1);
@ -128,12 +128,12 @@ void tst1() {
test_matrix(m2);
sparse_matrix<float, float> m10by9(10);
sparse_matrix<double, double> m10by9(10);
m10by9.set(0, 1, 1);
m10by9(0, 1) = 4;
float test = m10by9(0, 1);
double test = m10by9(0, 1);
std::cout << "got " << test << std::endl;
@ -629,7 +629,7 @@ void fill_matrix(sparse_matrix<T, X>& m){
}
void test_pivot_like_swaps_and_pivot(){
sparse_matrix<float, float> m(10);
sparse_matrix<double, double> m(10);
fill_matrix(m);
// print_matrix(m);
// pivot at 2,7
@ -648,12 +648,12 @@ void test_pivot_like_swaps_and_pivot(){
m.swap_rows(1, 3);
// print_matrix(m);
vector<float> row;
float alpha = 2.33;
vector<double> row;
double alpha = 2.33;
unsigned pivot_row = 1;
unsigned target_row = 2;
unsigned pivot_row_0 = 3;
float beta = 3.1;
double beta = 3.1;
m(target_row, 3) = 0;
m(target_row, 5) = 0;
m(pivot_row, 6) = 0;
@ -680,7 +680,7 @@ void test_pivot_like_swaps_and_pivot(){
#ifdef LEAN_DEBUG
void test_swap_rows() {
sparse_matrix<float, float> m(10);
sparse_matrix<double, double> m(10);
fill_matrix(m);
// print_matrix(m);
test_swap_rows(m, 3, 5);
@ -701,7 +701,7 @@ void test_swap_rows() {
test_swap_rows(m, 0, 7);
// go over some corner cases
sparse_matrix<float, float> m0(2);
sparse_matrix<double, double> m0(2);
test_swap_rows(m0, 0, 1);
m0(0, 0) = 3;
test_swap_rows(m0, 0, 1);
@ -709,7 +709,7 @@ void test_swap_rows() {
test_swap_rows(m0, 0, 1);
sparse_matrix<float, float> m1(10);
sparse_matrix<double, double> m1(10);
test_swap_rows(m1, 0, 1);
m1(0, 0) = 3;
test_swap_rows(m1, 0, 1);
@ -721,7 +721,7 @@ void test_swap_rows() {
test_swap_rows(m1, 0, 1);
sparse_matrix<float, float> m2(3);
sparse_matrix<double, double> m2(3);
test_swap_rows(m2, 0, 1);
m2(0, 0) = 3;
test_swap_rows(m2, 0, 1);
@ -833,7 +833,7 @@ void sparse_matrix_with_permutaions_test() {
}
void test_swap_columns() {
sparse_matrix<float, float> m(10);
sparse_matrix<double, double> m(10);
fill_matrix(m);
// print_matrix(m);
@ -853,7 +853,7 @@ void test_swap_columns() {
test_swap_columns(m, 0, 7);
// go over some corner cases
sparse_matrix<float, float> m0(2);
sparse_matrix<double, double> m0(2);
test_swap_columns(m0, 0, 1);
m0(0, 0) = 3;
test_swap_columns(m0, 0, 1);
@ -861,7 +861,7 @@ void test_swap_columns() {
test_swap_columns(m0, 0, 1);
sparse_matrix<float, float> m1(10);
sparse_matrix<double, double> m1(10);
test_swap_columns(m1, 0, 1);
m1(0, 0) = 3;
test_swap_columns(m1, 0, 1);
@ -873,7 +873,7 @@ void test_swap_columns() {
test_swap_columns(m1, 0, 1);
sparse_matrix<float, float> m2(3);
sparse_matrix<double, double> m2(3);
test_swap_columns(m2, 0, 1);
m2(0, 0) = 3;
test_swap_columns(m2, 0, 1);
@ -1015,6 +1015,13 @@ void update_settings(argument_parser & args_parser, lp_settings& settings) {
cout << "setting partial pivot constant to " << n << std::endl;
settings.c_partial_pivoting = n;
}
#if LEAN_DEBUG
if (get_int_from_args_parser("--ddd", args_parser, n)) {
cout << "setting ddd to " << n << std::endl;
settings.ddd = n;
}
#endif
if (get_int_from_args_parser("--density", args_parser, n)) {
double density = static_cast<double>(n) / 100.0;
cout << "setting density to " << density << std::endl;
@ -1029,8 +1036,8 @@ void update_settings(argument_parser & args_parser, lp_settings& settings) {
}
}
void setup_solver(unsigned max_iterations, unsigned time_limit, bool look_for_min, argument_parser & args_parser, lp_solver<double, double> * solver) {
template <typename T, typename X>
void setup_solver(unsigned max_iterations, unsigned time_limit, bool look_for_min, argument_parser & args_parser, lp_solver<T, X> * solver) {
if (max_iterations > 0)
solver->set_max_iterations_per_stage(max_iterations);
@ -1120,7 +1127,7 @@ void solve_mps_double(std::string file_name, bool look_for_min, unsigned max_ite
delete solver;
}
void solve_mps_rational(std::string file_name, bool look_for_min, unsigned max_iterations, unsigned time_limit, bool dual, argument_parser & /*args_parser*/) {
void solve_mps_rational(std::string file_name, bool look_for_min, unsigned max_iterations, unsigned time_limit, bool dual, argument_parser & args_parser) {
mps_reader<mpq, mpq> reader(file_name);
reader.read();
if (reader.is_ok()) {
@ -1128,14 +1135,8 @@ void solve_mps_rational(std::string file_name, bool look_for_min, unsigned max_i
if (look_for_min) {
solver->flip_costs();
}
setup_solver(max_iterations, time_limit, look_for_min, args_parser, solver);
int begin = get_millisecond_count();
if (max_iterations > 0) {
solver->set_max_iterations_per_stage(max_iterations);
}
if (time_limit > 0) {
solver->set_time_limit(time_limit);
}
solver->find_maximal_solution();
std::cout << "Status: " << lp_status_to_string(solver->get_status()) << std::endl;
@ -1756,7 +1757,7 @@ void test_init_U() {
}
void test_replace_column() {
sparse_matrix<float, float> m(10);
sparse_matrix<double, double> m(10);
fill_matrix(m);
m.swap_columns(0, 7);
m.swap_columns(6, 3);
@ -1765,7 +1766,7 @@ void test_replace_column() {
m(i, 0) = 0;
}
indexed_vector<float> w(m.dimension());
indexed_vector<double> w(m.dimension());
for (unsigned i = 0; i < m.dimension(); i++) {
w.set_value(i % 3, i);
}
@ -1783,6 +1784,7 @@ void test_replace_column() {
void setup_args_parser(argument_parser & parser) {
parser.add_option_with_after_string_with_help("--density", "the percentage of non-zeroes in the matrix below which it is not dense");
parser.add_option_with_after_string_with_help("--ddd", "the debug var");
parser.add_option_with_after_string_with_help("--harris_toler", "harris tolerance");
parser.add_option_with_help_string("--test_swaps", "test row swaps with a permutation");
parser.add_option_with_after_string_with_help("--checklu", "the file name for lu checking");
@ -1862,8 +1864,8 @@ struct mem_cpy_place_holder {
};
int finalize(unsigned ret) {
finalize_util_module();
finalize_numerics_module();
lean::finalize_util_module();
lean::finalize_numerics_module();
return ret;
}
@ -2450,8 +2452,8 @@ void test_square_dense_submatrix() {
}
int main(int argn, char * const * argv) {
initialize_util_module();
initialize_numerics_module();
lean::initialize_util_module();
lean::initialize_numerics_module();
int ret;
argument_parser args_parser(argn, argv);
setup_args_parser(args_parser);
@ -2560,5 +2562,5 @@ int main(int argn, char * const * argv) {
#endif
tst1();
std::cout<< "done with LP tests\n";
return finalize(has_violations() ? 1 : 0);
return finalize(lean::has_violations() ? 1 : 0);
}

99
src/tests/util/lp/mps_reader.h
View file

@ -8,22 +8,42 @@
#pragma once
// reads an MPS file reperesenting a Mixed Integer Program
#include <functional>
#include <algorithm>
#include <string>
#include <vector>
#include <unordered_map>
#include <iostream>
#include <fstream>
#include <locale>
#include "util/lp/lp_primal_simplex.h"
#include "util/lp/lp_dual_simplex.h"
#include "util/lp/lar_solver.h"
#include <iostream>
#include <fstream>
#include <util/numerics/mpq.h>
#include <functional>
#include <algorithm>
#include "util/lp/lp_utils.h"
#include "util/lp/lp_solver.h"
namespace lean {
using namespace std;
inline bool my_white_space(const char & a) {
return a == ' ' || a == '\t';
}
inline size_t number_of_whites(const string & s) {
size_t i = 0;
for (; i < s.size(); i++)
if (!my_white_space(s[i])) return i;
return i;
}
inline size_t number_of_whites_from_end(const string & s) {
size_t ret = 0;
for (int i = static_cast<int>(s.size()) - 1; i >= 0; i--)
if (my_white_space(s[i])) ret++; else break;
return ret;
}
// trim from start
inline string &ltrim(std::string &s) {
s.erase(s.begin(), find_if(s.begin(), s.end(), std::not1(std::ptr_fun<int, int>(std::isspace))));
inline std::string &ltrim(std::string &s) {
s.erase(0, number_of_whites(s));
return s;
}
@ -32,7 +52,8 @@ inline string &ltrim(std::string &s) {
// trim from end
inline string &rtrim(std::string &s) {
s.erase(std::find_if(s.rbegin(), s.rend(), std::not1(std::ptr_fun<int, int>(std::isspace))).base(), s.end());
// s.erase(std::find_if(s.rbegin(), s.rend(), std::not1(std::ptr_fun<int, int>(std::isspace))).base(), s.end());
s.erase(s.end() - number_of_whites_from_end(s), s.end());
return s;
}
// trim from both ends
@ -46,8 +67,8 @@ inline string trim(std::string const &r) {
}
inline vector<string> string_split(const string &source, const char *delimiter, bool keep_empty) {
vector<string> results;
inline std::vector<string> string_split(const string &source, const char *delimiter, bool keep_empty) {
std::vector<string> results;
size_t prev = 0;
size_t next = 0;
while ((next = source.find_first_of(delimiter, prev)) != std::string::npos) {
@ -62,9 +83,9 @@ inline vector<string> string_split(const string &source, const char *delimiter,
return results;
}
inline vector<string> split_and_trim(string line) {
inline std::vector<string> split_and_trim(string line) {
auto split = string_split(line, " \t", false);
vector<string> ret;
std::vector<string> ret;
for (auto s : split) {
ret.push_back(trim(s));
}
@ -220,7 +241,7 @@ class mps_reader {
goto fail;
} else {
row * row = t->second;
row->m_row_columns[column_name] = atof(column_data.substr(8).c_str());
row->m_row_columns[column_name] = from_string<T>(column_data.substr(8));
if (column_data.size() > 24) {
column_data = column_data.substr(25);
if (column_data.size() >= 22) {
@ -248,7 +269,7 @@ class mps_reader {
return;
}
row *r = t->second;
r->m_row_columns[column_name] = atof(tokens[i + 1].c_str());
r->m_row_columns[column_name] = from_string<T>(tokens[i + 1]);
}
}
@ -307,7 +328,7 @@ class mps_reader {
goto fail;
} else {
row * row = t->second;
row->m_right_side = atof(rhsides.substr(8).c_str());
row->m_right_side = from_string<T>(rhsides.substr(8));
if (rhsides.size() > 24) {
rhsides = rhsides.substr(25);
if (rhsides.size() >= 22) {
@ -333,7 +354,7 @@ class mps_reader {
return;
}
string rhsides = m_line.substr(14);
vector<string> splitted_line = split_and_trim(rhsides);
std::vector<string> splitted_line = split_and_trim(rhsides);
for (unsigned i = 0; i < splitted_line.size() - 1; i += 2) {
auto t = m_rows.find(splitted_line[i]);
@ -342,7 +363,7 @@ class mps_reader {
return;
}
row * row = t->second;
row->m_right_side = atof(splitted_line[i + 1].c_str());
row->m_right_side = from_string<X>(splitted_line[i + 1]);
}
}
@ -393,7 +414,7 @@ class mps_reader {
cout << "cannot find " << column_name << endl;
goto fail;
} else {
vector<string> bound_string;
std::vector<string> bound_string;
bound_string.push_back(column_name);
if (colstr.size() > 14) {
bound_string.push_back(colstr.substr(14));
@ -414,7 +435,7 @@ class mps_reader {
}
}
void update_bound(bound * b, vector<string> bound_string) {
void update_bound(bound * b, std::vector<string> bound_string) {
/*
UP means an upper bound is applied to the variable. A bound of type LO means a lower bound is applied. A bound type of FX ("fixed") means that the variable has upper and lower bounds equal to a single value. A bound type of FR ("free") means the variable has neither lower nor upper bounds and so can take on negative values. A variation on that is MI for free negative, giving an upper bound of 0 but no lower bound. Bound type PL is for a free positive for zero to plus infinity, but as this is the normal default, it is seldom used. There are also bound types for use in MIP models - BV for binary, being 0 or 1. UI for upper integer and LI for lower integer. SC stands for semi-continuous and indicates that the variable may be zero, but if not must be equal to at least the value given.
*/
@ -431,19 +452,16 @@ class mps_reader {
set_m_ok_to_false();
return;
}
double val = atof(bound_string[1].c_str());
b->m_upper_is_set = true;
b->m_upper= val;
b->m_upper= from_string<T>(bound_string[1]);
} else if (bound_type == "LO" || bound_type == "LI") {
if (bound_string.size() <= 1){
set_m_ok_to_false();
return;
}
double val = atof(bound_string[1].c_str());
b->m_low_is_set = true;
b->m_low = val;
b->m_low = from_string<T>(bound_string[1]);
} else if (bound_type == "FR") {
b->m_free = true;
} else if (bound_type == "FX") {
@ -452,9 +470,8 @@ class mps_reader {
return;
}
double val = atof(bound_string[1].c_str());
b->m_value_is_fixed = true;
b->m_fixed_value = val;
b->m_fixed_value = from_string<T>(bound_string[1]);
} else if (bound_type == "PL") {
b->m_low_is_set = true;
b->m_low = 0;
@ -471,7 +488,7 @@ class mps_reader {
void create_or_update_bound() {
const unsigned name_offset = 14;
lean_assert(m_line.size() >= 14);
vector<string> bound_string = split_and_trim(m_line.substr(name_offset, m_line.size()));
std::vector<string> bound_string = split_and_trim(m_line.substr(name_offset, m_line.size()));
if (bound_string.size() == 0) {
set_m_ok_to_false();
@ -514,7 +531,7 @@ class mps_reader {
goto fail;
} else {
row * row = t->second;
row->m_range = atof(rhsides.substr(8).c_str());
row->m_range = from_string<T>(rhsides.substr(8));
maybe_modify_current_row_and_add_row_for_range(row);
if (rhsides.size() > 24) {
rhsides = rhsides.substr(25);
@ -533,7 +550,7 @@ class mps_reader {
}
void read_range(vector<string> & splitted_line){
void read_range(std::vector<string> & splitted_line){
for (unsigned i = 1; i < splitted_line.size() - 1; i += 2) {
auto it = m_rows.find(splitted_line[i]);
if (it == m_rows.end()) {
@ -541,7 +558,7 @@ class mps_reader {
return;
}
row * row = it->second;
row->m_range = atof(splitted_line[i + 1].c_str());
row->m_range = from_string<X>(splitted_line[i + 1]);
maybe_modify_current_row_and_add_row_for_range(row);
}
}
@ -632,13 +649,13 @@ class mps_reader {
where |range| is range's absolute value.
*/
lp_relation get_relation_from_row(mps_reader::row_type row_type) {
switch (row_type) {
case mps_reader::row_type::Less_or_equal: return lp_relation::Less_or_equal;
case mps_reader::row_type::Greater_or_equal: return lp_relation::Greater_or_equal;
case mps_reader::row_type::Equal: return lp_relation::Equal;
lp_relation get_relation_from_row(row_type rt) {
switch (rt) {
case mps_reader::Less_or_equal: return lp_relation::Less_or_equal;
case mps_reader::Greater_or_equal: return lp_relation::Greater_or_equal;
case mps_reader::Equal: return lp_relation::Equal;
default:
std::cout << "Unexpected row_type " << row_type << endl;
std::cout << "Unexpected rt " << rt << endl;
set_m_ok_to_false();
throw;
}
@ -694,8 +711,8 @@ class mps_reader {
}
public:
vector<string> column_names() {
vector<string> v;
std::vector<string> column_names() {
std::vector<string> v;
for (auto s : m_columns) {
v.push_back(s.first);
}
@ -748,13 +765,13 @@ public:
return solver;
}
lconstraint_kind get_lar_relation_from_row(mps_reader::row_type row_type) {
switch (row_type) {
lconstraint_kind get_lar_relation_from_row(row_type rt) {
switch (rt) {
case Less_or_equal: return LE;
case Greater_or_equal: return GE;
case Equal: return EQ;
default:
std::cout << "Unexpected row_type " << row_type << endl;
std::cout << "Unexpected rt " << rt << endl;
set_m_ok_to_false();
throw;
}

32
src/tests/util/lp/smt_reader.h
View file

@ -16,29 +16,20 @@
#include "util/lp/lar_solver.h"
#include <iostream>
#include <fstream>
#include <util/numerics/mpq.h>
#include <functional>
#include <algorithm>
#include "tests/util/lp/mps_reader.h"
#include "util/lp/canonic_left_side.h"
#include "util/lp/lar_constraints.h"
#include <sstream>
#include <cstdlib>
namespace lean {
using namespace std;
template<typename T>
T from_string(const std::string& str) {
std::istringstream ss(str);
T ret;
ss >> ret;
return ret;
}
class smt_reader {
public:
struct lisp_elem {
string m_head;
vector<lisp_elem> m_elems;
std::vector<lisp_elem> m_elems;
void print() {
if (m_elems.size()) {
cout << '(';
@ -51,12 +42,12 @@ namespace lean {
cout << " " << m_head;
}
}
unsigned size() const { return m_elems.size(); }
unsigned size() const { return static_cast<unsigned>(m_elems.size()); }
bool is_simple() const { return size() == 0; }
};
struct formula_constraint {
lconstraint_kind m_kind;
vector<pair<mpq, string>> m_coeffs;
std::vector<pair<mpq, string>> m_coeffs;
mpq m_right_side = numeric_traits<mpq>::zero();
void add_pair(mpq c, string name) {
m_coeffs.push_back(make_pair(c, name));
@ -65,13 +56,12 @@ namespace lean {
lisp_elem m_formula_lisp_elem;
vector<formula_constraint> m_constraints;
std::vector<formula_constraint> m_constraints;
string m_file_name;
ifstream m_file_stream;
string m_line;
bool m_is_OK = true;
unsigned m_line_number = 0;
smt_reader(string file_name):
m_file_name(file_name), m_file_stream(file_name) {
}
@ -90,9 +80,9 @@ namespace lean {
}
int first_separator() {
unsigned blank_pos = m_line.find(' ');
unsigned br_pos = m_line.find('(');
unsigned reverse_br_pos = m_line.find(')');
unsigned blank_pos = static_cast<unsigned>(m_line.find(' '));
unsigned br_pos = static_cast<unsigned>(m_line.find('('));
unsigned reverse_br_pos = static_cast<unsigned>(m_line.find(')'));
return min(blank_pos, min(br_pos, reverse_br_pos));
}
@ -150,7 +140,7 @@ namespace lean {
void parse_line() {
if (m_line.find(":formula") == 0) {
int first_br = m_line.find('(');
int first_br = static_cast<int>(m_line.find('('));
if (first_br == -1) {
cout << "empty formula" << endl;
return;
@ -214,7 +204,7 @@ namespace lean {
}
void add_mult_elem(formula_constraint & c, vector<lisp_elem> & els) {
void add_mult_elem(formula_constraint & c, std::vector<lisp_elem> & els) {
lean_assert(els.size() == 2);
mpq coeff = get_coeff(els[0]);
string col_name = get_name(els[1]);
@ -249,7 +239,7 @@ namespace lean {
}
}
void add_sum(formula_constraint & c, vector<lisp_elem> & sum_els) {
void add_sum(formula_constraint & c, std::vector<lisp_elem> & sum_els) {
for (auto & el : sum_els)
add_sum_elem(c, el);
}

2
src/util/lp/binary_heap_priority_queue.cpp
View file

@ -169,7 +169,7 @@ template <typename T> void binary_heap_priority_queue<T>::put_the_last_at_the_to
}
/// return the first element of the queue and removes it from the queue
template <typename T> unsigned binary_heap_priority_queue<T>::dequeue() {
lean_assert(m_heap_size);
lean_assert(m_heap_size > 0);
int ret = m_heap[1];
put_the_last_at_the_top_and_fix_the_heap();
m_heap_inverse[ret] = -1;

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

@ -1,3 +1,4 @@
/*
Copyright (c) 2013 Microsoft Corporation. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.
@ -7,6 +8,7 @@
#pragma once
#include <vector>
#include "util/debug.h"
#include "util/lp/lp_utils.h"
namespace lean {
// the elements with the smallest priority are dequeued first
template <typename T>

5
src/util/lp/binary_heap_upair_queue.cpp
View file

@ -6,10 +6,11 @@
*/
#include <set>
#include "util/lp/lp_utils.h"
#include "util/lp/binary_heap_upair_queue.h"
namespace lean {
template <typename T> binary_heap_upair_queue<T>::binary_heap_upair_queue(unsigned size) : m_q(size), m_pairs(size) {
lean_assert(size);
lean_assert(size > 0);
for (unsigned i = 0; i < size; i++)
m_available_spots.push_back(i);
}
@ -49,7 +50,7 @@ template <typename T> void binary_heap_upair_queue<T>::enqueue(unsigned i, unsig
// it is a new pair, let us find a spot for it
if (m_available_spots.empty()) {
// we ran out of empty spots
unsigned size_was = m_pairs.size();
unsigned size_was = static_cast<unsigned>(m_pairs.size());
unsigned new_size = size_was << 1;
for (unsigned i = size_was; i < new_size; i++)
m_available_spots.push_back(i);

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

@ -13,7 +13,7 @@
#include <set>
#include <utility>
#include "util/lp/binary_heap_priority_queue.h"
#include "util/lp/hash_helper.h"
typedef std::pair<unsigned, unsigned> upair;

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

@ -10,9 +10,7 @@
#include <string>
#include <algorithm>
#include <utility>
#include "util/numerics/mpq.h"
#include "util/lp/column_info.h"
#include "util/lp/hash_helper.h"
namespace lean {
typedef unsigned var_index;
typedef unsigned constraint_index;
@ -48,7 +46,7 @@ public:
m_column_info.set_name(name);
}
unsigned size() const { return m_coeffs.size(); }
unsigned size() const { return static_cast<unsigned>(m_coeffs.size()); }
void normalize() {
if (m_coeffs.size() == 0) return;

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

@ -100,7 +100,7 @@ public:
}
T get_fixed_value() const {
lean_assert(m_is_fixed)
lean_assert(m_is_fixed);
return m_fixed_value;
}

44
src/util/lp/core_solver_pretty_printer.cpp
View file

@ -7,8 +7,10 @@
#include <limits>
#include <string>
#include <algorithm>
#include "util/lp/lp_utils.h"
#include "util/lp/lp_core_solver_base.h"
#include "util/lp/core_solver_pretty_printer.h"
#include "util/lp/numeric_pair.h"
namespace lean {
@ -32,7 +34,7 @@ core_solver_pretty_printer<T, X>::core_solver_pretty_printer(lp_core_solver_base
m_cost_title = "costs";
m_basis_heading_title = "heading";
m_x_title = "x*";
m_title_width = std::max(std::max(m_cost_title.size(), std::max(m_basis_heading_title.size(), m_x_title.size())), m_approx_norm_title.size());
m_title_width = static_cast<unsigned>(std::max(std::max(m_cost_title.size(), std::max(m_basis_heading_title.size(), m_x_title.size())), m_approx_norm_title.size()));
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_costs() {
@ -52,9 +54,9 @@ template <typename T, typename X> core_solver_pretty_printer<T, X>::~core_solver
delete [] m_ed_buff;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_rs_width() {
m_rs_width = T_to_string(m_core_solver.get_cost()).size();
m_rs_width = static_cast<unsigned>(T_to_string(m_core_solver.get_cost()).size());
for (unsigned i = 0; i < nrows(); i++) {
auto wt = T_to_string(m_rs[i]).size();
unsigned wt = static_cast<unsigned>(T_to_string(m_rs[i]).size());
if (wt > m_rs_width) {
m_rs_width = wt;
}
@ -81,7 +83,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_m_
name);
m_rs[row] += m_core_solver.m_ed[row] * m_core_solver.m_x[column];
}
m_exact_column_norms.push_back(current_column_norm() + 1);
m_exact_column_norms.push_back(current_column_norm() + T(1)); // a conversion missing 1 -> T
}
}
@ -122,11 +124,11 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::adjust_
template <typename T, typename X> unsigned core_solver_pretty_printer<T, X>:: get_column_width(unsigned column) {
unsigned w = std::max(m_costs[column].size(), T_to_string(m_core_solver.m_x[column]).size());
unsigned w = static_cast<unsigned>(std::max(m_costs[column].size(), T_to_string(m_core_solver.m_x[column]).size()));
adjust_width_with_bounds(column, w);
adjust_width_with_basis_heading(column, w);
for (unsigned i = 0; i < nrows(); i++) {
unsigned cellw = m_A[i][column].size();
unsigned cellw = static_cast<unsigned>(m_A[i][column].size());
if (cellw > w) {
w = cellw;
}
@ -171,14 +173,14 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_x
return;
}
int blanks = m_title_width + 1 - m_x_title.size();
int blanks = m_title_width + 1 - static_cast<int>(m_x_title.size());
m_out << m_x_title;
print_blanks(blanks, m_out);
auto bh = m_core_solver.m_x;
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(bh[i]);
int blanks = m_column_widths[i] - s.size();
int blanks = m_column_widths[i] - static_cast<int>(s.size());
print_blanks(blanks, m_out);
m_out << s << " "; // the column interval
}
@ -217,13 +219,13 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_l
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - m_low_bounds_title.size();
int blanks = m_title_width + 1 - static_cast<unsigned>(m_low_bounds_title.size());
m_out << m_low_bounds_title;
print_blanks(blanks, m_out);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_low_bound_string(i);
int blanks = m_column_widths[i] - s.size();
int blanks = m_column_widths[i] - static_cast<unsigned>(s.size());
print_blanks(blanks, m_out);
m_out << s << " "; // the column interval
}
@ -234,13 +236,13 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_u
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - m_upp_bounds_title.size();
int blanks = m_title_width + 1 - static_cast<unsigned>(m_upp_bounds_title.size());
m_out << m_upp_bounds_title;
print_blanks(blanks, m_out);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_upp_bound_string(i);
int blanks = m_column_widths[i] - s.size();
int blanks = m_column_widths[i] - static_cast<unsigned>(s.size());
print_blanks(blanks, m_out);
m_out << s << " "; // the column interval
}
@ -248,12 +250,12 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_u
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_exact_norms() {
int blanks = m_title_width + 1 - m_exact_norm_title.size();
int blanks = m_title_width + 1 - static_cast<int>(m_exact_norm_title.size());
m_out << m_exact_norm_title;
print_blanks(blanks, m_out);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_exact_column_norm_string(i);
int blanks = m_column_widths[i] - s.size();
int blanks = m_column_widths[i] - static_cast<int>(s.size());
print_blanks(blanks, m_out);
m_out << s << " ";
}
@ -261,12 +263,12 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_e
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_approx_norms() {
int blanks = m_title_width + 1 - m_approx_norm_title.size();
int blanks = m_title_width + 1 - static_cast<int>(m_approx_norm_title.size());
m_out << m_approx_norm_title;
print_blanks(blanks, m_out);
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();
int blanks = m_column_widths[i] - static_cast<int>(s.size());
print_blanks(blanks, m_out);
m_out << s << " ";
}
@ -289,7 +291,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print()
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_basis_heading() {
int blanks = m_title_width + 1 - m_basis_heading_title.size();
int blanks = m_title_width + 1 - static_cast<int>(m_basis_heading_title.size());
m_out << m_basis_heading_title;
print_blanks(blanks, m_out);
@ -299,7 +301,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_b
auto bh = m_core_solver.m_basis_heading;
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(bh[i]);
int blanks = m_column_widths[i] - s.size();
int blanks = m_column_widths[i] - static_cast<unsigned>(s.size());
print_blanks(blanks, m_out);
m_out << s << " "; // the column interval
}
@ -307,7 +309,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_b
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_cost() {
int blanks = m_title_width + 1 - m_cost_title.size();
int blanks = m_title_width + 1 - static_cast<int>(m_cost_title.size());
m_out << m_cost_title;
print_blanks(blanks, m_out);
print_given_rows(m_costs, m_cost_signs, m_core_solver.get_cost());
@ -317,7 +319,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_g
for (unsigned col = 0; col < row.size(); col++) {
unsigned width = m_column_widths[col];
string s = row[col];
int number_of_blanks = width - s.size();
int number_of_blanks = width - static_cast<unsigned>(s.size());
lean_assert(number_of_blanks >= 0);
print_blanks(number_of_blanks, m_out);
m_out << s << ' ';
@ -328,7 +330,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_g
m_out << '=';
string rs = T_to_string(rst);
int nb = m_rs_width - rs.size();
int nb = m_rs_width - static_cast<int>(rs.size());
lean_assert(nb >= 0);
print_blanks(nb + 1, m_out);
m_out << rs << std::endl;

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

@ -5,9 +5,6 @@
Author: Lev Nachmanson
*/
#pragma once
#include "util/numerics/float.h"
#include "util/numerics/double.h"
#include "util/numerics/mpq.h"
namespace lean {
template <typename T>

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

@ -8,17 +8,9 @@
#pragma once
#include <vector>
#include "util/debug.h"
#include "util/numerics/numeric_traits.h"
#include "util/numerics/xnumeral.h"
#include "util/numerics/mpq.h"
#include "util/numerics/mpz.h"
#include "util/numerics/mpbq.h"
#include "util/numerics/double.h"
#include "util/numerics/float.h"
#include "util/numerics/mpfp.h"
#include <string>
#include "util/lp/sparse_vector.h"
#include <iomanip>
#include "util/lp/lp_utils.h"
namespace lean {
template <typename T> void print_vector(const std::vector<T> & t, std::ostream & out);

3
src/util/lp/indexed_vector_instances.cpp
View file

@ -11,7 +11,6 @@ template void indexed_vector<double>::clear();
template void indexed_vector<double>::clear_all();
template void indexed_vector<double>::erase_from_index(unsigned int);
template void indexed_vector<double>::set_value(double, unsigned int);
template void indexed_vector<float>::set_value(float, unsigned int);
template void indexed_vector<mpq>::clear();
template void indexed_vector<mpq>::clear_all();
template void indexed_vector<mpq>::erase_from_index(unsigned int);
@ -19,6 +18,6 @@ template void indexed_vector<mpq>::resize(unsigned int);
template void indexed_vector<mpq>::set_value(mpq, unsigned int);
#ifdef LEAN_DEBUG
template bool indexed_vector<double>::is_OK() const;
template bool lean::indexed_vector<mpq>::is_OK() const;
template bool indexed_vector<mpq>::is_OK() const;
#endif
}

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

@ -8,13 +8,10 @@
#pragma once
#include <vector>
#include <utility>
#include "util/debug.h"
#include "util/buffer.h"
#include "util/numerics/numeric_traits.h"
#include "util/numerics/xnumeral.h"
#include <unordered_map>
#include <string>
#include <algorithm>
#include "util/lp/lp_utils.h"
#include "util/lp/canonic_left_side.h"
namespace lean {
inline lconstraint_kind flip_kind(lconstraint_kind t) {
@ -30,6 +27,7 @@ inline std::string lconstraint_kind_string(lconstraint_kind t) {
case EQ: return std::string("=");
}
lean_unreachable();
return std::string(); // it is unreachable
}
class lar_base_constraint {
@ -54,7 +52,7 @@ public:
lar_constraint(const lar_base_constraint & c);
unsigned size() const {
return m_left_side.size();
return static_cast<unsigned>(m_left_side.size());
}
buffer<std::pair<mpq, var_index>> get_left_side_coefficients() const;

3
src/util/lp/lar_core_solver.cpp
View file

@ -647,6 +647,7 @@ template <typename T, typename X> bool lar_core_solver<T, X>::improves_pivot_
default:
lean_unreachable();
}
return false; // it is unreachable
}
template <typename T, typename X> int lar_core_solver<T, X>::choose_entering_column_for_row_inf_strategy(int inf_sign) {
@ -744,7 +745,7 @@ template <typename T, typename X> void lar_core_solver<T, X>::advance_on_infe
int entering_delta_sign;
std::vector<unsigned> leaving_candidates;
X delta = find_initial_delta_and_its_sign(inf_row, entering, inf_sign, entering_delta_sign, leaving_candidates);
lean_assert(leaving_candidates.size());
lean_assert(leaving_candidates.size() > 0);
lean_assert(delta > zero_of_type<X>());
unsigned row = my_random() % this->m_m;
unsigned initial_row = row;

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

@ -7,6 +7,7 @@
#pragma once
#include <vector>
#include <string>
#include <utility>
#include "util/lp/lp_core_solver_base.h"
#include <algorithm>
#include "util/lp/indexed_vector.h"
@ -23,7 +24,7 @@ class lar_core_solver : public lp_core_solver_base<T, X> {
std::vector<unsigned> m_tight_basic_columns;
std::vector<breakpoint<X>> m_breakpoints;
binary_heap_priority_queue<X> m_breakpoint_indices_queue;
std::vector<pair<mpq, unsigned>> m_infeasible_row;
std::vector<std::pair<mpq, unsigned>> m_infeasible_row;
int m_infeasible_row_sign = 0;
// with a breakpoint at this delta
public:
@ -38,7 +39,7 @@ public:
int get_infeasible_row_sign() const { return m_infeasible_row_sign; }
const std::vector<pair<mpq, unsigned>> & get_infeasibility_info(int & inf_sign) const {
const std::vector<std::pair<mpq, unsigned>> & get_infeasibility_info(int & inf_sign) const {
inf_sign = m_infeasible_row_sign;
return m_infeasible_row;
}
@ -155,7 +156,7 @@ public:
std::vector<unsigned> & leaving_candidates);
unsigned find_leaving_for_inf_row_strategy(std::vector<unsigned> & leaving_candidates) {
lean_assert(leaving_candidates.size());
lean_assert(leaving_candidates.size() > 0);
return leaving_candidates[my_random() % leaving_candidates.size()]; // more randomness
}

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

@ -7,6 +7,7 @@
#include <string>
#include <algorithm>
#include <vector>
#include <utility>
#include "util/lp/lar_solver.h"
namespace lean {
double conversion_helper <double>::get_low_bound(const column_info<mpq> & ci) {
@ -29,7 +30,7 @@ double conversion_helper <double>::get_upper_bound(const column_info<mpq> & ci)
return ci.get_upper_bound().get_double() - eps;
}
canonic_left_side * lar_solver::create_or_fetch_existing_left_side(const buffer<pair<mpq, var_index>>& left_side_par) {
canonic_left_side * lar_solver::create_or_fetch_existing_left_side(const buffer<std::pair<mpq, var_index>>& left_side_par) {
auto left_side = new canonic_left_side(left_side_par);
lean_assert(left_side->size() > 0);
auto it = m_canonic_left_sides.find(left_side);
@ -60,7 +61,7 @@ mpq lar_solver::find_ratio(canonic_left_side * ls, const lar_constraint & constr
}
void lar_solver::add_canonic_left_side_for_var(var_index i, std::string var_name) {
buffer<pair<mpq, var_index>> b;
buffer<std::pair<mpq, var_index>> b;
b.push_back(std::make_pair(numeric_traits<mpq>::one(), i));
auto can_ls = new canonic_left_side(b);
can_ls->set_name(var_name);
@ -74,7 +75,7 @@ void lar_solver::map_left_side_to_column_of_A(canonic_left_side* left_side, uns
lean_assert(m_column_indices_to_canonic_left_sides.find(j) == m_column_indices_to_canonic_left_sides.end());
left_side->m_column_index = j; // assigning this index does not change the hash of canonic_left_side
if (left_side->size() > 1) { // if size is one we will not create a row for this left side
left_side->m_row_index = m_basis.size();
left_side->m_row_index = static_cast<int>(m_basis.size());
m_basis.push_back(j); // j will be a basis column, so we put it into the basis as well
}
m_column_indices_to_canonic_left_sides[j++] = left_side; // pointing from the column to the left side
@ -105,8 +106,8 @@ void lar_solver::add_row_to_A(static_matrix<U, V> & A, unsigned i, canonic_left_
template <typename U, typename V>
void lar_solver::create_matrix_A(static_matrix<U, V> & A) {
unsigned n = m_column_indices_to_canonic_left_sides.size();
unsigned m = m_basis.size();
unsigned n = static_cast<unsigned>(m_column_indices_to_canonic_left_sides.size());
unsigned m = static_cast<unsigned>(m_basis.size());
A.init_empty_matrix(m, n);
unsigned i = 0;
for (auto & t : m_column_indices_to_canonic_left_sides)
@ -244,7 +245,7 @@ void lar_solver::fill_column_types() {
template <typename V>
void lar_solver::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
unsigned n = static_cast<unsigned>(m_canonic_left_sides.size()); // this is the number of columns
lb.resize(n);
ub.resize(n);
for (auto t : m_canonic_left_sides) {
@ -284,6 +285,7 @@ template <typename V> V lar_solver::get_column_val(std::vector<V> & low_bound, s
default:
lean_unreachable();
}
return zero_of_type<V>(); // it is unreachable
}
lar_solver::~lar_solver() {
@ -307,8 +309,8 @@ var_index lar_solver::add_var(std::string s) {
return i;
}
constraint_index lar_solver::add_constraint(const buffer<pair<mpq, var_index>>& left_side, lconstraint_kind kind_par, mpq right_side_par) {
lean_assert(left_side.size());
constraint_index lar_solver::add_constraint(const buffer<std::pair<mpq, var_index>>& left_side, lconstraint_kind kind_par, mpq right_side_par) {
lean_assert(left_side.size() > 0);
constraint_index i = m_available_constr_index++;
lar_constraint original_constr(left_side, kind_par, right_side_par, i);
canonic_left_side * ls = create_or_fetch_existing_left_side(left_side);
@ -341,6 +343,7 @@ bool lar_solver::constraint_holds(const lar_constraint & constr, std::unordered_
default:
lean_unreachable();
}
return false; // it is unreachable
}
void lar_solver::solve_with_core_solver() {
@ -352,7 +355,7 @@ void lar_solver::solve_with_core_solver() {
#endif
}
bool lar_solver::the_relations_are_of_same_type(const buffer<pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum) {
bool lar_solver::the_relations_are_of_same_type(const buffer<std::pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum) {
unsigned n_of_G = 0, n_of_L = 0;
bool strict = false;
for (auto & it : evidence) {
@ -381,7 +384,7 @@ void lar_solver::register_in_map(std::unordered_map<var_index, mpq> & coeffs, la
else p->second += it.second * a;
}
}
bool lar_solver::the_left_sides_sum_to_zero(const buffer<pair<mpq, unsigned>> & evidence) {
bool lar_solver::the_left_sides_sum_to_zero(const buffer<std::pair<mpq, unsigned>> & evidence) {
std::unordered_map<var_index, mpq> coeff_map;
for (auto & it : evidence) {
mpq coeff = it.first;
@ -395,7 +398,7 @@ bool lar_solver::the_left_sides_sum_to_zero(const buffer<pair<mpq, unsigned>> &
return true;
}
bool lar_solver::the_righ_sides_do_not_sum_to_zero(const buffer<pair<mpq, unsigned>> & evidence) {
bool lar_solver::the_righ_sides_do_not_sum_to_zero(const buffer<std::pair<mpq, unsigned>> & evidence) {
mpq ret = numeric_traits<mpq>::zero();
for (auto & it : evidence) {
mpq coeff = it.first;
@ -407,7 +410,7 @@ bool lar_solver::the_righ_sides_do_not_sum_to_zero(const buffer<pair<mpq, unsign
}
#ifdef LEAN_DEBUG
bool lar_solver::the_evidence_is_correct() {
buffer<pair<mpq, unsigned>> evidence;
buffer<std::pair<mpq, unsigned>> evidence;
get_infeasibility_evidence(evidence);
lconstraint_kind kind;
lean_assert(the_relations_are_of_same_type(evidence, kind));
@ -442,7 +445,7 @@ void lar_solver::init_right_sides_with_zeros(std::vector<V> & rs) {
rs.resize(m_basis.size(), zero_of_type<V>());
}
mpq lar_solver::sum_of_right_sides_of_evidence(const buffer<pair<mpq, unsigned>> & evidence) {
mpq lar_solver::sum_of_right_sides_of_evidence(const buffer<std::pair<mpq, unsigned>> & evidence) {
mpq ret = numeric_traits<mpq>::zero();
for (auto & it : evidence) {
mpq coeff = it.first;
@ -499,7 +502,7 @@ void lar_solver::find_solution_signature_with_doubles(lar_solution_signature & s
A.clear();
create_matrix_A(A);
for (auto & s : column_scale_vector)
s = one_of_type<double>();
s = 1.0;
}
std::vector<double> costs(A.column_count());
auto core_solver = lp_primal_core_solver<double, double>(A,
@ -545,7 +548,7 @@ lp_status lar_solver::check() {
solve();
return m_status;
}
void lar_solver::get_infeasibility_evidence(buffer<pair<mpq, constraint_index>> & evidence){
void lar_solver::get_infeasibility_evidence(buffer<std::pair<mpq, constraint_index>> & evidence){
if (!m_mpq_core_solver.get_infeasible_row_sign()) {
return;
}
@ -556,8 +559,8 @@ void lar_solver::get_infeasibility_evidence(buffer<pair<mpq, constraint_index>>
get_infeasibility_evidence_for_inf_sign(evidence, inf_row, inf_sign);
}
void lar_solver::get_infeasibility_evidence_for_inf_sign(buffer<pair<mpq, constraint_index>> & evidence,
const std::vector<pair<mpq, unsigned>> & inf_row,
void lar_solver::get_infeasibility_evidence_for_inf_sign(buffer<std::pair<mpq, constraint_index>> & evidence,
const std::vector<std::pair<mpq, unsigned>> & inf_row,
int inf_sign) {
for (auto & it : inf_row) {
mpq coeff = it.first;
@ -727,6 +730,7 @@ mpq lar_solver::get_infeasibility_of_constraint(const lar_normalized_constraint
default:
lean_unreachable();
}
return mpq(0); // it is unreachable
}
mpq lar_solver::get_canonic_left_side_val(canonic_left_side * ls, std::unordered_map<std::string, mpq> & solution) {

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

@ -6,10 +6,9 @@
*/
#pragma once
#include <vector>
#include <utility>
#include "util/debug.h"
#include "util/buffer.h"
#include "util/numerics/numeric_traits.h"
#include "util/numerics/xnumeral.h"
#include <unordered_map>
#include <unordered_set>
#include <string>
@ -61,7 +60,7 @@ class lar_solver {
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);
canonic_left_side * create_or_fetch_existing_left_side(const buffer<std::pair<mpq, var_index>>& left_side_par);
bool var_is_fixed(unsigned j);
@ -148,7 +147,7 @@ public:
var_index add_var(std::string s);
constraint_index add_constraint(const buffer<pair<mpq, var_index>>& left_side, lconstraint_kind kind_par, mpq right_side_par);
constraint_index add_constraint(const buffer<std::pair<mpq, var_index>>& left_side, lconstraint_kind kind_par, mpq right_side_par);
bool is_infeasible(const column_info<mpq> & ci) {
return ci.low_bound_is_set() && ci.upper_bound_is_set() && ci.get_low_bound() > ci.get_upper_bound();
@ -162,12 +161,12 @@ public:
void solve_with_core_solver();
bool the_relations_are_of_same_type(const buffer<pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum);
bool the_relations_are_of_same_type(const buffer<std::pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum);
void register_in_map(std::unordered_map<var_index, mpq> & coeffs, lar_constraint & cn, const mpq & a);
bool the_left_sides_sum_to_zero(const buffer<pair<mpq, unsigned>> & evidence);
bool the_left_sides_sum_to_zero(const buffer<std::pair<mpq, unsigned>> & evidence);
bool the_righ_sides_do_not_sum_to_zero(const buffer<pair<mpq, unsigned>> & evidence);
bool the_righ_sides_do_not_sum_to_zero(const buffer<std::pair<mpq, unsigned>> & evidence);
bool the_evidence_is_correct();
@ -176,7 +175,7 @@ public:
template <typename V>
void init_right_sides_with_zeros(std::vector<V> & rs);
mpq sum_of_right_sides_of_evidence(const buffer<pair<mpq, unsigned>> & evidence);
mpq sum_of_right_sides_of_evidence(const buffer<std::pair<mpq, unsigned>> & evidence);
void prepare_independently_of_numeric_type();
template <typename U, typename V>
@ -201,10 +200,10 @@ public:
void solve();
lp_status check();
void get_infeasibility_evidence(buffer<pair<mpq, constraint_index>> & evidence);
void get_infeasibility_evidence(buffer<std::pair<mpq, constraint_index>> & evidence);
void get_infeasibility_evidence_for_inf_sign(buffer<pair<mpq, constraint_index>> & evidence,
const std::vector<pair<mpq, unsigned>> & inf_row,
void get_infeasibility_evidence_for_inf_sign(buffer<std::pair<mpq, constraint_index>> & evidence,
const std::vector<std::pair<mpq, unsigned>> & inf_row,
int inf_sign);

2
src/util/lp/lar_solver_instances.cpp
View file

@ -4,5 +4,5 @@
Author: Lev Nachmanson
*/
#
#include "util/lp/lar_solver.cpp"

17
src/util/lp/lp_core_solver_base.cpp
View file

@ -7,6 +7,7 @@
#include <set>
#include <string>
#include <vector>
#include "util/lp/lp_utils.h"
#include "util/lp/lp_core_solver_base.h"
namespace lean {
void init_basic_part_of_basis_heading(std::vector<unsigned> & basis, unsigned m, std::vector<int> & basis_heading) {
@ -21,7 +22,7 @@ void init_non_basic_part_of_basis_heading(std::vector<int> & basis_heading, std:
if (basis_heading[j] < 0) {
non_basic_columns.push_back(j);
// the index of column j in m_non_basic_columns is (- basis_heading[j] - 1)
basis_heading[j] = - non_basic_columns.size();
basis_heading[j] = - static_cast<int>(non_basic_columns.size());
}
}
}
@ -116,7 +117,7 @@ solve_yB(std::vector<T> & y) {
template <typename T, typename X> void lp_core_solver_base<T, X>::
update_index_of_ed() {
m_index_of_ed.clear();
unsigned i = m_ed.size();
unsigned i = static_cast<unsigned>(m_ed.size());
while (i--) {
if (!is_zero(m_ed[i]))
m_index_of_ed.push_back(i);
@ -248,7 +249,7 @@ calculate_pivot_row_when_pivot_row_of_B1_is_ready() {
}
}
unsigned j = m_pivot_row.size();
unsigned j = static_cast<unsigned>(m_pivot_row.size());
while (j--) {
if (!is_zero(m_pivot_row[j]))
m_pivot_row_index.push_back(j);
@ -274,7 +275,7 @@ print_statistics(X cost) {
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::
print_statistics_with_iterations_and_check_that_the_time_is_over(unsigned total_iterations) {
if (total_iterations % m_settings.report_frequency == 0) {
if (m_settings.print_statistics && total_iterations % m_settings.report_frequency == 0) {
std::cout << "iterations = " << total_iterations << ", nonzeros = " << m_factorization->get_number_of_nonzeroes() << std::endl;
if (time_is_over()) {
return true;
@ -356,6 +357,8 @@ column_is_dual_feasible(unsigned j) const {
std::cout << "unexpected column type = " << column_type_to_string(m_column_type[j]) << std::endl;
lean_unreachable();
}
lean_unreachable();
return false;
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::
d_is_not_negative(unsigned j) const {
@ -419,7 +422,9 @@ update_basis_and_x(int entering, int leaving, X const & tt) {
m_refactor_counter = 0;
m_iters_with_no_cost_growing++;
if (m_factorization->get_status() != LU_status::OK) {
throw exception(sstream() << "failing refactor on off_result for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << m_total_iterations);
std::stringstream s;
s << "failing refactor on off_result for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << m_total_iterations;
throw_exception(s.str());
}
return false;
}
@ -691,6 +696,8 @@ get_non_basic_column_value_position(unsigned j) {
default:
lean_unreachable();
}
lean_unreachable();
return at_low_bound;
}
template <typename T, typename X> void lp_core_solver_base<T, X>::init_lu() {

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

@ -8,8 +8,7 @@
#include <set>
#include <vector>
#include <string>
#include "util/sstream.h"
#include "util/exception.h"
#include "util/lp/lp_utils.h"
#include "util/lp/core_solver_pretty_printer.h"
#include "util/lp/numeric_pair.h"
#include "util/lp/static_matrix.h"

14
src/util/lp/lp_dual_core_solver.cpp
View file

@ -55,7 +55,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::fill_non_basis
while (j--) {
if (this->m_basis_heading[j] >= 0 || !m_can_enter_basis[j]) continue;
nb.push_back(j);
this->m_basis_heading[j] = - nb.size();
this->m_basis_heading[j] = - static_cast<int>(nb.size());
}
}
@ -68,7 +68,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::restore_non_ba
if (m_can_enter_basis[j]) {
lean_assert(std::find(nb.begin(), nb.end(), j) == nb.end());
nb.push_back(j);
this->m_basis_heading[j] = - nb.size();
this->m_basis_heading[j] = - static_cast<int>(nb.size());
}
}
}
@ -183,6 +183,8 @@ template <typename T, typename X> T lp_dual_core_solver<T, X>::pricing_for_row(u
default:
lean_unreachable();
}
lean_unreachable();
return numeric_traits<T>::zero();
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::pricing_loop(unsigned number_of_rows_to_try, unsigned offset_in_rows) {
@ -283,6 +285,8 @@ template <typename T, typename X> int lp_dual_core_solver<T, X>::define_sign_of_
default:
lean_unreachable();
}
lean_unreachable();
return 0;
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::can_be_breakpoint(unsigned j) {
@ -329,7 +333,7 @@ template <typename T, typename X> T lp_dual_core_solver<T, X>::get_delta() {
return this->m_x[m_p] - this->m_low_bound_values[m_p];
}
if (this->x_above_upper_bound(m_p)) {
return this->m_x[m_p] - this->m_upper_bound_values[m_p];;
return this->m_x[m_p] - this->m_upper_bound_values[m_p];
}
lean_unreachable();
case low_bound:
@ -343,10 +347,12 @@ template <typename T, typename X> T lp_dual_core_solver<T, X>::get_delta() {
}
lean_unreachable();
case fixed:
return this->m_x[m_p] - this->m_upper_bound_values[m_p];;
return this->m_x[m_p] - this->m_upper_bound_values[m_p];
default:
lean_unreachable();
}
lean_unreachable();
return zero_of_type<T>();
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::restore_d() {

12
src/util/lp/lp_dual_simplex.cpp
View file

@ -5,7 +5,7 @@
Author: Lev Nachmanson
*/
#include "util/lp/lp_dual_simplex.h"
namespace lean {
namespace lean{
template <typename T, typename X> void lp_dual_simplex<T, X>::decide_on_status_after_stage1() {
switch (m_core_solver->get_status()) {
@ -193,9 +193,13 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_bounds_
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:
throw exception(sstream() << "unexpected bound type " << j << " "
<< column_type_to_string(get_column_type(j)));
case upper_bound: {
std::stringstream s;
s << "unexpected bound type " << j << " "
<< column_type_to_string(get_column_type(j));
throw_exception(s.str());
break;
}
case low_bound: {
m_can_enter_basis[j] = true;
this->set_scaled_cost(j);

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

@ -6,8 +6,7 @@
*/
#pragma once
#include <vector>
#include "util/sstream.h"
#include "util/exception.h"
#include "util/lp/lp_utils.h"
#include "util/lp/lp_solver.h"
#include "util/lp/lp_dual_core_solver.h"
namespace lean {

37
src/util/lp/lp_primal_core_solver.cpp
View file

@ -30,7 +30,7 @@ void lp_primal_core_solver<T, X>::sort_non_basis() {
// reinit m_basis_heading
for (unsigned j = 0; j < this->m_non_basic_columns.size(); j++) {
unsigned col = this->m_non_basic_columns[j];
this->m_basis_heading[col] = - j - 1;
this->m_basis_heading[col] = - static_cast<int>(j) - 1;
m_non_basis_list.push_back(col);
}
}
@ -38,7 +38,7 @@ void lp_primal_core_solver<T, X>::sort_non_basis() {
template <typename T, typename X>
int lp_primal_core_solver<T, X>::choose_entering_column(unsigned number_of_benefitial_columns_to_go_over) { // at this moment m_y = cB * B(-1)
lean_assert(number_of_benefitial_columns_to_go_over);
lean_assert(number_of_benefitial_columns_to_go_over > 0);
if (this->m_sort_counter == 0) {
sort_non_basis();
this->m_sort_counter = 20;
@ -162,7 +162,8 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leavi
}
if (++i == this->m_m) i = 0;
} while ( i != initial_i);
restore_harris_eps();
if (!precise<T>())
restore_harris_eps();
return leaving;
}
@ -189,7 +190,7 @@ template <typename T, typename X> bool lp_primal_core_solver<T, X>::try_jump_
}
template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_and_t(unsigned entering, X & t){
if (numeric_traits<T>::precise()) return find_leaving_and_t_precisely(entering, t);
// if (numeric_traits<T>::precise()) return find_leaving_and_t_precisely(entering, t);
X theta = get_harris_theta();
lean_assert(theta >= zero_of_type<X>());
if (try_jump_to_another_bound_on_entering(entering, theta, t)) return entering;
@ -251,7 +252,7 @@ template <typename T, typename X> lp_primal_core_solver<T, X>::lp_primal_core_so
std::vector<X> & upper_bound_values,
lp_settings & settings,
std::unordered_map<unsigned, std::string> const & column_names):
lp_core_solver_base<T, X>(A, b,
lp_core_solver_base<T, X>(A, b,
basis,
x,
costs,
@ -260,8 +261,12 @@ template <typename T, typename X> lp_primal_core_solver<T, X>::lp_primal_core_so
column_type_array,
low_bound_values,
upper_bound_values),
m_beta(A.row_count()),
m_converted_harris_eps(convert_struct<T, double>::convert(this->m_settings.harris_feasibility_tolerance)) {
m_beta(A.row_count()) {
if (!(numeric_traits<T>::precise())) {
m_converted_harris_eps = convert_struct<T, double>::convert(this->m_settings.harris_feasibility_tolerance);
} else {
m_converted_harris_eps = zero_of_type<T>();
}
this->m_status = UNKNOWN;
this->m_column_norm_update_counter = settings.column_norms_update_frequency;
}
@ -504,7 +509,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
// this->init_reduced_costs_for_one_iteration();
// }
template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_entering(int entering) {
template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_entering(int entering) {
lean_assert(entering > -1);
this->solve_Bd(entering);
int refresh_result = refresh_reduced_cost_at_entering_and_check_that_it_is_off(entering);
@ -545,7 +550,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::push_forw
}
template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::get_number_of_non_basic_column_to_try_for_enter() {
unsigned ret = this->m_non_basic_columns.size();
unsigned ret = static_cast<unsigned>(this->m_non_basic_columns.size());
if (this->m_status == TENTATIVE_UNBOUNDED)
return ret; // we really need to find entering with a large reduced cost
if (ret > 300) {
@ -663,7 +668,11 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::update_or
// following Swietanowski - A new steepest ...
template <typename T, typename X> void lp_primal_core_solver<T, X>::update_column_norms(unsigned entering, unsigned leaving) {
T pivot = this->m_pivot_row[entering];
T g_ent = std::max(calculate_norm_of_entering_exactly() / pivot /pivot, T(0.000001));
T g_ent = calculate_norm_of_entering_exactly() / pivot / pivot;
if (!numeric_traits<T>::precise()) {
if (g_ent < T(0.000001))
g_ent = T(0.000001);
}
this->m_column_norms[leaving] = g_ent;
for (unsigned j : this->m_pivot_row_index) {
@ -733,7 +742,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::add_colum
}
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::one_iteration() {
template <typename T, typename X> void lp_primal_core_solver<T, X>::one_iteration() {
this->m_total_iterations++;
unsigned number_of_benefitial_columns_to_go_over = get_number_of_non_basic_column_to_try_for_enter();
int entering = choose_entering_column(number_of_benefitial_columns_to_go_over);
@ -743,7 +752,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::one_itera
advance_on_entering(entering);
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::fill_breakpoints_array(unsigned entering) {
template <typename T, typename X> void lp_primal_core_solver<T, X>::fill_breakpoints_array(unsigned entering) {
clear_breakpoints();
for (unsigned i = this->m_m; i--;)
try_add_breakpoint_in_row(i);
@ -846,7 +855,8 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::update_ba
else
this->update_x(entering, delta);
}
template <typename T, typename X> bool lp_primal_core_solver<T, X>::column_is_feasible(unsigned j) const {
template <typename T, typename X> bool lp_primal_core_solver<T, X>::column_is_feasible(unsigned j) const {
const X& x = this->m_x[j];
switch (this->m_column_type[j]) {
case fixed:
@ -879,6 +889,7 @@ template <typename T, typename X> bool lp_primal_core_solver<T, X>::column_is
default:
lean_unreachable();
}
return false; // it is unreachable
}
template <typename T, typename X> bool lp_primal_core_solver<T, X>::calc_current_x_is_feasible() const {

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

@ -16,7 +16,6 @@
#include <math.h>
#include <cstdlib>
#include <algorithm>
#include "util/numerics/double.h"
#include "util/lp/lu.h"
#include "util/lp/lp_solver.h"
#include "util/lp/static_matrix.h"
@ -43,7 +42,7 @@ public:
bool m_recalc_reduced_costs = false;
std::set<unsigned> m_forbidden_enterings;
std::vector<T> m_beta; // see Swietanowski working vector beta for column norms
T m_epsilon_of_reduced_cost = T(1e-7);
T m_epsilon_of_reduced_cost = T(1)/T(10000000);
bool m_exit_on_feasible_solution = false;
std::vector<T> m_costs_backup;
bool m_current_x_is_feasible;
@ -56,7 +55,7 @@ public:
int find_leaving_and_t_precisely(unsigned entering, X & t);
static X positive_infinity() {
return convert_struct<X, double>::convert(std::numeric_limits<unsigned>::max());
return convert_struct<X, unsigned>::convert(std::numeric_limits<unsigned>::max());
}
X get_harris_theta();
@ -91,7 +90,7 @@ public:
limit_inf_on_bound_m_pos(m, this->m_x[j], this->m_upper_bound_values[j], theta);
};
X harris_eps_for_bound(const X & bound) const { return ( convert_struct<X, double>::convert(1) + abs(bound)/10) * m_converted_harris_eps/3;
X harris_eps_for_bound(const X & bound) const { return ( convert_struct<X, int>::convert(1) + abs(bound)/10) * m_converted_harris_eps/3;
}

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

@ -31,7 +31,7 @@ template <typename T, typename X> void lp_primal_simplex<T, X>::init_buffer(unsi
template <typename T, typename X> void lp_primal_simplex<T, X>::refactor() {
m_core_solver->init_lu();
if (m_core_solver->factorization()->get_status() != LU_status::OK) {
throw exception("cannot refactor");
throw_exception("cannot refactor");
}
}
@ -226,6 +226,10 @@ template <typename T, typename X> void lp_primal_simplex<T, X>::fill_A_x_and_bas
template <typename T, typename X> void lp_primal_simplex<T, X>::solve_with_total_inf() {
int total_vars = this->m_A->column_count() + this->row_count();
if (total_vars == 0) {
this->m_status = OPTIMAL;
return;
}
m_low_bounds.clear();
m_low_bounds.resize(total_vars, zero_of_type<X>()); // low bounds are shifted ot zero
this->m_x.resize(total_vars, numeric_traits<T>::zero());
@ -237,10 +241,9 @@ template <typename T, typename X> void lp_primal_simplex<T, X>::solve_with_total
this->fill_column_names_for_core_solver();
unsigned j = this->m_A->column_count() - 1;
unsigned core_solver_cols = this->number_of_core_structurals();
while (j >= core_solver_cols)
while (j >= core_solver_cols) {
this->m_costs[j--] = numeric_traits<T>::zero();
}
set_scaled_costs();
m_core_solver = new lp_primal_core_solver<T, X>(*this->m_A,
this->m_b,
@ -267,7 +270,9 @@ template <typename T, typename X> bool lp_primal_simplex<T, X>::bounds_hold(std:
for (auto it : this->m_columns) {
auto sol_it = solution.find(it.second->get_name());
if (sol_it == solution.end()) {
throw exception(sstream() << "cannot find column " << it.first << " in solutio");
std::stringstream s;
s << "cannot find column " << it.first << " in solution";
throw_exception(s.str() );
}
if (!it.second->bounds_hold(sol_it->second)) {
@ -282,7 +287,9 @@ template <typename T, typename X> bool lp_primal_simplex<T, X>::bounds_hold(std:
template <typename T, typename X> T lp_primal_simplex<T, X>::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()) {
throw exception(sstream() << "cannot find row " << i);
std::stringstream s;
s << "cannot find row " << i;
throw_exception(s.str() );
}
T ret = numeric_traits<T>::zero();
for (auto & pair : it->second) {
@ -340,6 +347,7 @@ template <typename T, typename X> bool lp_primal_simplex<T, X>::row_constraint_h
return true;;
}
lean_unreachable();
return false; // it is unreachable
}
template <typename T, typename X> bool lp_primal_simplex<T, X>::row_constraints_hold(std::unordered_map<std::string, T> const & solution) {

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

@ -9,8 +9,7 @@
#include <unordered_map>
#include <string>
#include <algorithm>
#include "util/exception.h"
#include "util/sstream.h"
#include "util/lp/lp_utils.h"
#include "util/lp/column_info.h"
#include "util/lp/lp_primal_core_solver.h"
#include "util/lp/lp_solver.h"

10
src/util/lp/lp_settings.cpp
View file

@ -7,7 +7,6 @@
#include <cmath>
#include <string>
#include <vector>
#include "util/numerics/double.h"
#include "util/lp/lp_settings.h"
namespace lean {
std::string column_type_to_string(column_type t) {
@ -25,6 +24,7 @@ std::string column_type_to_string(column_type t) {
default:
lean_unreachable();
}
return "unknown"; // it is unreachable
}
const char* lp_status_to_string(lp_status status) {
@ -44,6 +44,7 @@ const char* lp_status_to_string(lp_status status) {
default:
lean_unreachable();
}
return "UNKNOWN"; // it is unreachable
}
lp_status lp_status_from_string(std::string status) {
@ -57,6 +58,7 @@ lp_status lp_status_from_string(std::string status) {
if (status == "ITERATIONS_EXHAUSTED") return lp_status::ITERATIONS_EXHAUSTED;
if (status == "EMPTY") return lp_status::EMPTY;
lean_unreachable();
return lp_status::UNKNOWN; // it is unreachable
}
int get_millisecond_count() {
timeb tb;
@ -71,14 +73,14 @@ int get_millisecond_span(int start_time) {
return span;
}
void my_random_init(unsigned * seed) {
#ifdef LEAN_WINDOWS
#if defined(LEAN_WINDOWS) || defined(WIN32)
srand(*seed);
#else
rand_r(seed);
#endif
}
unsigned my_random() {
#ifdef LEAN_WINDOWS
#if defined(LEAN_WINDOWS) || defined(WIN32)
return rand(); // NOLINT
#else
return lrand48();
@ -129,7 +131,7 @@ bool vectors_are_equal(const std::vector<T> & a, const buffer<T> &b) {
template <typename T>
bool vectors_are_equal(const std::vector<T> & a, const std::vector<T> &b) {
unsigned n = a.size();
unsigned n = static_cast<unsigned>(a.size());
if (n != b.size()) return false;
if (numeric_traits<T>::precise()) {
for (unsigned i = 0; i < n; i ++){

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

@ -11,10 +11,7 @@
#include <algorithm>
#include <limits>
#include <sys/timeb.h>
#include "util/debug.h"
#include "util/numerics/numeric_traits.h"
#include "util/sstream.h"
#include "util/numerics/mpq.h"
#include "util/lp/lp_utils.h"
namespace lean {
@ -51,11 +48,6 @@ lp_status lp_status_from_string(std::string status);
enum non_basic_column_value_position { at_low_bound, at_upper_bound, at_fixed, free_of_bounds };
template <typename X> bool is_epsilon_small(const X & v, const double& eps); // forward definition
template <typename X> bool precise() { return numeric_traits<X>::precise();}
template <typename X> X zero_of_type() { return numeric_traits<X>::zero(); }
template <typename X> X one_of_type() { return numeric_traits<X>::one(); }
template <typename X> bool is_zero(const X & v) { return numeric_traits<X>::is_zero(v); }
template <typename X> double get_double(const X & v) { return numeric_traits<X>::get_double(v); }
struct lp_settings {
// when the absolute value of an element is less than pivot_epsilon
@ -66,7 +58,7 @@ struct lp_settings {
// a quatation "if some choice of the entering vairable leads to an eta matrix
// whose diagonal element in the eta column is less than e2 (entering_diag_epsilon) in magnitude, the this choice is rejected ...
double entering_diag_epsilon = 1e-8;
double c_partial_pivoting = 10; // this is the constant c from page 410
int c_partial_pivoting = 10; // this is the constant c from page 410
unsigned depth_of_rook_search = 4;
bool using_partial_pivoting = true;
// dissertation of Achim Koberstein
@ -151,7 +143,8 @@ struct lp_settings {
// the method of lar solver to use
bool row_feasibility = true;
bool use_double_solver_for_lar = true;
int report_frequency = 10000000;
int report_frequency = 1000;
bool print_statistics = false;
unsigned column_norms_update_frequency = 1000;
bool scale_with_ratio = true;
double density_threshold = 0.7; // need to tune it up, todo
@ -210,4 +203,3 @@ inline void print_blanks(int n, std::ostream & out) {
while (n--) {out << ' '; }
}
}

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

@ -43,7 +43,9 @@ template <typename T, typename X> void lp_solver<T, X>::give_symbolic_name_to_co
template <typename T, typename X> T lp_solver<T, X>::get_column_value_by_name(std::string name) const {
auto it = m_names_to_columns.find(name);
if (it == m_names_to_columns.end()) {
throw exception(sstream() << "get_column_value_by_name " << name);
std::stringstream s;
s << "get_column_value_by_name " << name;
throw_exception(s.str());
}
return get_column_value(it -> second);
}
@ -249,7 +251,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_e_is_obsolete(std
return false;
}
template <typename T, typename X> int lp_solver<T, X>::row_ge_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index) {
template <typename T, typename X> bool lp_solver<T, X>::row_ge_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index) {
T rs = m_constraints[row_index].m_rs;
if (row_is_zero(row)) {
if (rs > zero_of_type<X>())
@ -274,7 +276,7 @@ template <typename T, typename X> int lp_solver<T, X>::row_ge_is_obsolete(std
return false;
}
template <typename T, typename X> bool lp_solver<T, X>::row_le_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index) {
template <typename T, typename X> bool lp_solver<T, X>::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)) {
@ -318,6 +320,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_is_obsolete(std::
return row_le_is_obsolete(row, row_index);
}
lean_unreachable();
return false; // it is unreachable
}
template <typename T, typename X> void lp_solver<T, X>::remove_fixed_or_zero_columns() {
@ -349,7 +352,7 @@ template <typename T, typename X> void lp_solver<T, X>::remove_fixed_or_zero_col
}
}
template <typename T, typename X> unsigned lp_solver<T, X>::try_to_remove_some_rows() {
template <typename T, typename X> unsigned lp_solver<T, X>::try_to_remove_some_rows() {
std::vector<unsigned> rows_to_delete;
for (auto & t : m_A_values) {
if (row_is_obsolete(t.second, t.first)) {
@ -366,13 +369,13 @@ template <typename T, typename X> unsigned lp_solver<T, X>::try_to_remove_som
}
}
remove_fixed_or_zero_columns();
return rows_to_delete.size();
return static_cast<unsigned>(rows_to_delete.size());
}
template <typename T, typename X> void lp_solver<T, X>::cleanup() {
int n = 0; // number of deleted rows
int d;
while ((d = try_to_remove_some_rows()))
while ((d = try_to_remove_some_rows() > 0))
n += d;
// if (n == 1)
@ -395,7 +398,7 @@ template <typename T, typename X> void lp_solver<T, X>::map_external_columns_to_
for (auto & row : m_A_values) {
for (auto & col : row.second) {
if (col.second == numeric_traits<T>::zero() || m_columns[col.first]->is_fixed()) {
throw exception("found fixed column");
throw_exception("found fixed column");
}
unsigned j = col.first;
auto j_place = m_external_columns_to_core_solver_columns.find(j);
@ -426,7 +429,7 @@ template <typename T, typename X> void lp_solver<T, X>::unscale() {
}
template <typename T, typename X> void lp_solver<T, X>::fill_A_from_A_values() {
m_A = new static_matrix<T, X>(m_A_values.size(), number_of_core_structurals());
m_A = new static_matrix<T, X>(static_cast<unsigned>(m_A_values.size()), number_of_core_structurals());
for (auto & t : m_A_values) {
lean_assert(m_external_rows_to_core_solver_rows.find(t.first) != m_external_rows_to_core_solver_rows.end());
unsigned row = m_external_rows_to_core_solver_rows[t.first];
@ -461,22 +464,19 @@ template <typename T, typename X> void lp_solver<T, X>::count_slacks_and_artific
template <typename T, typename X> void lp_solver<T, X>::count_slacks_and_artificials_for_row(unsigned i) {
lean_assert(this->m_constraints.find(this->m_core_solver_rows_to_external_rows[i]) != this->m_constraints.end());
auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[i]];
T rs;
switch (constraint.m_relation) {
case Equal:
m_artificials++;
break;
case Greater_or_equal:
m_slacks++;
rs = this->m_b[i];
if (rs > 0) {
if (this->m_b[i] > 0) {
m_artificials++;
}
break;
case Less_or_equal:
m_slacks++;
rs = this->m_b[i];
if (rs < 0) {
if (this->m_b[i] < 0) {
m_artificials++;
}
break;
@ -488,7 +488,7 @@ template <typename T, typename X> T lp_solver<T, X>::low_bound_shift_for_row(
auto row = this->m_A_values.find(i);
if (row == this->m_A_values.end()) {
throw exception("cannot find row");
throw_exception("cannot find row");
}
for (auto col : row->second) {
ret += col.second * this->m_columns[col.first]->get_shift();

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

@ -10,8 +10,6 @@
#include <unordered_map>
#include <algorithm>
#include <vector>
#include "util/exception.h"
#include "util/sstream.h"
#include "util/lp/lp_settings.h"
#include "util/lp/column_info.h"
#include "util/lp/static_matrix.h"
@ -173,7 +171,7 @@ protected:
bool row_e_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index);
int row_ge_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index);
bool row_ge_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index);
bool row_le_is_obsolete(std::unordered_map<unsigned, T> & row, unsigned row_index);
@ -199,7 +197,7 @@ protected:
void fill_column_names_for_core_solver();
unsigned number_of_core_structurals() { return m_external_columns_to_core_solver_columns.size(); }
unsigned number_of_core_structurals() { return static_cast<unsigned>(m_external_columns_to_core_solver_columns.size()); }
void restore_column_scales_to_one() {

19
src/util/lp/lu.cpp
View file

@ -8,6 +8,7 @@
#include <algorithm>
#include <set>
#include <vector>
#include <utility>
#include "util/lp/lu.h"
namespace lean {
#ifdef LEAN_DEBUG
@ -33,7 +34,7 @@ template<typename T, typename X>
void print_matrix(static_matrix<T, X> &m, std::ostream & out) {
std::vector<std::vector<std::string>> A;
std::vector<unsigned> widths;
std::set<pair<unsigned, unsigned>> domain = m.get_domain();
std::set<std::pair<unsigned, unsigned>> domain = m.get_domain();
for (unsigned i = 0; i < m.row_count(); i++) {
A.push_back(std::vector<std::string>());
for (unsigned j = 0; j < m.column_count(); j++) {
@ -169,9 +170,13 @@ void lu<T, X>::solve_Bd_when_w_is_ready(std::vector<T> & d, indexed_vector<T>& w
template <typename T, typename X>
template <typename L>
void lu<T, X>::solve_By_when_y_is_ready(std::vector<L> & y) {
if (numeric_traits<T>::precise()) {
m_U.solve_U_y(y);
m_R.apply_reverse_from_left(y); // see 24.3 from Chvatal
return;
}
m_U.double_solve_U_y(y);
m_R.apply_reverse_from_left(y); // see 24.3 from Chvatal
if (precise<X>()) return;
unsigned i = m_dim;
while (i--) {
if (is_zero(y[i])) continue;
@ -232,7 +237,7 @@ void lu<T, X>::solve_yB_internal(std::vector<T>& y) {
}
template <typename T, typename X>
void lu<T, X>::add_delta_to_solution(std::vector<T>& yc, std::vector<T>& y){
unsigned i = y.size();
unsigned i = static_cast<unsigned>(y.size());
while (i--)
y[i]+=yc[i];
}
@ -249,6 +254,10 @@ void lu<T, X>::find_error_of_yB(std::vector<T>& yc, const std::vector<T>& y) {
// y is the input
template <typename T, typename X>
void lu<T, X>::solve_yB(std::vector<T> & y) {
if (numeric_traits<T>::precise()) {
solve_yB_internal(y);
return;
}
std::vector<T> yc(y); // copy y aside
solve_yB_internal(y);
find_error_of_yB(yc, y);
@ -427,7 +436,7 @@ void lu<T, X>::check_apply_lp_lists_to_w(T * w) {
template <typename T, typename X>
void lu<T, X>::process_column(int j) {
unsigned pi, pj;
m_U.get_pivot_for_column(pi, pj, T(m_settings.c_partial_pivoting), j);
m_U.get_pivot_for_column(pi, pj, m_settings.c_partial_pivoting, j);
if (static_cast<int>(pi) == -1) {
std::cout << "cannot find the pivot for column " << j << std::endl;
m_failure = true;
@ -693,7 +702,7 @@ void lu<T, X>::replace_column(unsigned leaving, T pivot_elem, indexed_vector<T>
push_matrix_to_tail(row_eta);
}
calculate_Lwave_Pwave_for_bump(replaced_column, lowest_row_of_the_bump);
lean_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
lean_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
lean_assert(w.is_OK() && m_row_eta_work_vector.is_OK());
}
template <typename T, typename X>

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

@ -9,8 +9,6 @@
#include <vector>
#include "util/debug.h"
#include "util/numerics/numeric_traits.h"
#include "util/numerics/xnumeral.h"
#include <algorithm>
#include <set>
#include "util/lp/sparse_matrix.h"

1
src/util/lp/lu_instances.cpp
View file

@ -33,7 +33,6 @@ template void lean::init_factorization<lean::mpq, lean::mpq>(lean::lu<lean::mpq,
template void lean::init_factorization<lean::mpq, lean::numeric_pair<lean::mpq> >(lean::lu<lean::mpq, lean::numeric_pair<lean::mpq> >*&, lean::static_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >&, std::vector<unsigned int, std::allocator<unsigned int> >&, std::vector<int, std::allocator<int> >&, lean::lp_settings&, std::vector<unsigned int, std::allocator<unsigned int> >&);
#ifdef LEAN_DEBUG
template void lean::print_matrix<double, double>(lean::sparse_matrix<double, double>&, std::ostream & out);
template void lean::print_matrix<float, float>(lean::sparse_matrix<float, float>&, std::ostream & out);
template void lean::print_matrix<double, double>(lean::static_matrix<double, double>&, std::ostream & out);
template bool lean::lu<double, double>::is_correct();
template lean::dense_matrix<double, double> lean::get_B<double, double>(lean::lu<double, double>&);

2
src/util/lp/matrix.h
View file

@ -7,7 +7,7 @@
#ifdef LEAN_DEBUG
#pragma once
#include "util/numerics/numeric_traits.h"
#include "util/numerics/double.h"
// #include "util/numerics/double.h"
#include <vector>
#include <string>
#include "util/lp/lp_settings.h"

93
src/util/lp/numeric_pair.h
View file

@ -3,25 +3,67 @@
Released under Apache 2.0 license as described in the file LICENSE.
Author: Lev Nachmanson
The idea is that it is only one different file in Lean and z3 source inside of LP
*/
#pragma once
// #define lp_for_z3
#include <string>
#include <cmath>
#include <algorithm>
#include "util/numerics/mpq.h"
#include "util/numerics/double.h"
#include "util/lp/lp_settings.h"
#ifdef lp_for_z3
#include "../rational.h"
#include "../sstream.h"
#include "../z3_exception.h"
#else
#include "util/numerics/mpq.h"
#include "util/numerics/numeric_traits.h"
#endif
namespace lean {
#ifdef lp_for_z3 // rename rationals
typedef rational mpq;
#else
typedef lean::mpq mpq;
#endif
template <typename T>
std::string T_to_string(const T & t); // forward definition
#ifdef lp_for_z3
template <typename T> class numeric_traits {};
template <> class numeric_traits<double> {
public:
static bool precise() { return false; }
static double g_zero;
static double const &zero() { return g_zero; }
static double g_one;
static double const &one() { return g_one; }
static bool is_zero(double v) { return v == 0.0; }
static double const & get_double(double const & d) { return d;}
static double log(double const & d) { NOT_IMPLEMENTED_YET(); return d;}
static double from_string(std::string const & str) { return atof(str.c_str()); }
};
template<>
class numeric_traits<rational> {
public:
static bool precise() { return true; }
static rational const & zero() { return rational::zero(); }
static rational const & one() { return rational::one(); }
static bool is_zero(const rational & v) { return v.is_zero(); }
static double const get_double(const rational & d) { return d.get_double();}
static rational log(rational const& r) { UNREACHABLE(); return r; }
static rational from_string(std::string const & str) { return rational(str.c_str()); }
};
#endif
template <typename X, typename Y>
struct convert_struct {
static X convert(const Y & y){ return X(y);}
static bool is_epsilon_small(const X & x, const double & y) { return std::abs(get_double(x)) < y; }
static bool below_bound_numeric(const X &, const X &, const Y &) { lean_unreachable(); }
static bool above_bound_numeric(const X &, const X &, const Y &) { lean_unreachable(); }
static bool is_epsilon_small(const X & x, const double & y) { return std::abs(numeric_traits<X>::get_double(x)) < y; }
static bool below_bound_numeric(const X &, const X &, const Y &) { /*lean_unreachable();*/ return false;}
static bool above_bound_numeric(const X &, const X &, const Y &) { /*lean_unreachable();*/ return false; }
};
@ -75,7 +117,7 @@ struct numeric_pair {
}
numeric_pair operator/(const numeric_pair &) const {
lean_unreachable();
// lean_unreachable();
}
@ -84,7 +126,7 @@ struct numeric_pair {
}
numeric_pair operator*(const numeric_pair & /*a*/) const {
lean_unreachable();
// lean_unreachable();
}
numeric_pair& operator+=(const numeric_pair & a) {
@ -119,7 +161,7 @@ struct numeric_pair {
return numeric_pair(-x, -y);
}
static bool precize() { return numeric_traits<T>::precize();}
static bool precize() { return lean::numeric_traits<T>::precize();}
std::string to_string() const { return std::string("(") + T_to_string(x) + ", " + T_to_string(y) + ")"; }
};
@ -141,20 +183,23 @@ numeric_pair<T> operator*(const numeric_pair<T> & r, const X & a) {
return numeric_pair<T>(a * r.x, a * r.y);
}
} // close namespace lean
namespace lean {
// template <numeric_pair, typename T> bool precise() { return numeric_traits<T>::precise();}
template <typename T> double get_double(const numeric_pair<T> & ) { lean_unreachable(); }
template <typename T> double get_double(const lean::numeric_pair<T> & ) { /* lean_unreachable(); */ return 0;}
template <typename T>
class numeric_traits<numeric_pair<T>> {
class numeric_traits<lean::numeric_pair<T>> {
public:
static bool precise() { return numeric_traits<T>::precise();}
static numeric_pair<T> zero() { return numeric_pair<T>(numeric_traits<T>::zero(), numeric_traits<T>::zero()); }
static bool is_zero(const numeric_pair<T> & v) { return numeric_traits<T>::is_zero(v.x) && numeric_traits<T>::is_zero(v.y); }
static double get_double(const numeric_pair<T> & v){ return numeric_traits<T>::get_double(v.x); } // just return the double of the first coordinate
static double one() { lean_unreachable(); }
static lean::numeric_pair<T> zero() { return lean::numeric_pair<T>(numeric_traits<T>::zero(), numeric_traits<T>::zero()); }
static bool is_zero(const lean::numeric_pair<T> & v) { return numeric_traits<T>::is_zero(v.x) && numeric_traits<T>::is_zero(v.y); }
static double get_double(const lean::numeric_pair<T> & v){ return numeric_traits<T>::get_double(v.x); } // just return the double of the first coordinate
static double one() { /*lean_unreachable();*/ return 0;}
};
} // close namespace lean
namespace lean {
template <>
struct convert_struct<double, numeric_pair<double>> {
static double convert(const numeric_pair<double> & q) {return q.x;}
@ -165,19 +210,20 @@ template <typename X> bool is_epsilon_small(const X & v, const double& eps); /
template <typename T>
struct convert_struct<numeric_pair<T>, double> {
static numeric_pair<T> convert(const double & q) {
return numeric_pair<T>(convert_struct<T, double>::convert(q), zero_of_type<T>());
return numeric_pair<T>(convert_struct<T, double>::convert(q), numeric_traits<T>::zero());
}
static bool is_epsilon_small(const numeric_pair<T> & p, const double & eps) {
return convert_struct<T, double>::is_epsilon_small(p.x, eps) && convert_struct<T, double>::is_epsilon_small(p.y, eps);
}
static bool below_bound_numeric(const numeric_pair<T> &, const numeric_pair<T> &, const double &) {
lean_unreachable();
// lean_unreachable();
return false;
}
static bool above_bound_numeric(const numeric_pair<T> &, const numeric_pair<T> &, const double &) {
lean_unreachable();
// lean_unreachable();
return false;
}
};
template <>
struct convert_struct<numeric_pair<double>, double> {
static numeric_pair<double> convert(const double & q) {
@ -188,7 +234,6 @@ struct convert_struct<numeric_pair<double>, double> {
}
static int compare_on_coord(const double & x, const double & bound, const double eps) {
lean_assert(eps > 0);
if (bound == 0) return (x < - eps)? -1: (x > eps? 1 : 0); // it is an important special case
double relative = (bound > 0)? - eps: eps;
return (x < bound * (1.0 + relative) - eps)? -1 : ((x > bound * (1.0 - relative) + eps)? 1 : 0);
@ -199,7 +244,6 @@ struct convert_struct<numeric_pair<double>, double> {
if (r == 1) return false;
if (r == -1) return true;
// the first coordinates are almost the same
lean_assert(r == 0);
return compare_on_coord(x.y, bound.y, eps) == -1;
}
@ -208,7 +252,6 @@ struct convert_struct<numeric_pair<double>, double> {
if (r == -1) return false;
if (r == 1) return true;
// the first coordinates are almost the same
lean_assert(r == 0);
return compare_on_coord(x.y, bound.y, eps) == 1;
}
};
@ -220,13 +263,11 @@ struct convert_struct<double, double> {
}
static double convert(const double & y){ return y;}
static bool below_bound_numeric(const double & x, const double & bound, const double & eps) {
lean_assert(eps > 0);
if (bound == 0) return x < - eps;
double relative = (bound > 0)? - eps: eps;
return x < bound * (1.0 + relative) - eps;
}
static bool above_bound_numeric(const double & x, const double & bound, const double & eps) {
lean_assert(eps > 0);
if (bound == 0) return x > eps;
double relative = (bound > 0)? eps: - eps;
return x > bound * (1.0 + relative) + eps;

8
src/util/lp/permutation_matrix.cpp
View file

@ -77,7 +77,7 @@ void permutation_matrix<T, X>::apply_from_left_perm(indexed_vector<L> & w, lp_se
copy_aside(t, tmp_index, w);
clear_data(w);
// set the new values
for (unsigned i = t.size(); i > 0;) {
for (unsigned i = static_cast<unsigned>(t.size()); i > 0;) {
i--;
unsigned j = m_rev[tmp_index[i]];
w[j] = t[i];
@ -112,7 +112,7 @@ template <typename T, typename X> void permutation_matrix<T, X>::apply_from_righ
template <typename T, typename X> template <typename L>
void permutation_matrix<T, X>::copy_aside(std::vector<L> & t, std::vector<unsigned> & tmp_index, indexed_vector<L> & w) {
for (unsigned i = t.size(); i > 0;) {
for (unsigned i = static_cast<unsigned>(t.size()); i > 0;) {
i--;
unsigned j = w.m_index[i];
t[i] = w[j]; // copy aside all non-zeroes
@ -123,7 +123,7 @@ void permutation_matrix<T, X>::copy_aside(std::vector<L> & t, std::vector<unsign
template <typename T, typename X> template <typename L>
void permutation_matrix<T, X>::clear_data(indexed_vector<L> & w) {
// clear old non-zeroes
for (unsigned i = w.m_index.size(); i > 0;) {
for (unsigned i = static_cast<unsigned>(w.m_index.size()); i > 0;) {
i--;
unsigned j = w.m_index[i];
w[j] = zero_of_type<L>();
@ -145,7 +145,7 @@ void permutation_matrix<T, X>::apply_reverse_from_left(indexed_vector<L> & w) {
clear_data(w);
// set the new values
for (unsigned i = t.size(); i > 0;) {
for (unsigned i = static_cast<unsigned>(t.size()); i > 0;) {
i--;
unsigned j = m_permutation[tmp_index[i]];
w[j] = t[i];

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

@ -8,14 +8,6 @@
#include <vector>
#include <algorithm>
#include "util/debug.h"
#include "util/numerics/numeric_traits.h"
#include "util/numerics/xnumeral.h"
#include "util/numerics/mpq.h"
#include "util/numerics/mpz.h"
#include "util/numerics/mpbq.h"
#include "util/numerics/double.h"
#include "util/numerics/float.h"
#include "util/numerics/mpfp.h"
#include <string>
#include "util/lp/sparse_vector.h"
#include "util/lp/indexed_vector.h"
@ -129,7 +121,7 @@ namespace lean {
bool is_identity() const;
unsigned size() const { return m_rev.size(); }
unsigned size() const { return static_cast<unsigned>(m_rev.size()); }
unsigned * values() const { return m_permutation; }
}; // end of the permutation class

7
src/util/lp/permutation_matrix_instances.cpp
View file

@ -16,10 +16,7 @@ template void lean::permutation_matrix<double, double>::multiply_by_permutation_
template void lean::permutation_matrix<double, double>::multiply_by_reverse_from_right(lean::permutation_matrix<double, double>&);
template lean::permutation_matrix<double, double>::permutation_matrix(unsigned int, std::vector<unsigned int, std::allocator<unsigned int> > const&);
template void lean::permutation_matrix<double, double>::transpose_from_left(unsigned int, unsigned int);
template void lean::permutation_matrix<float, float>::apply_from_right(std::vector<float, std::allocator<float> >&);
template lean::permutation_matrix<float, float>::permutation_matrix(unsigned int);
template void lean::permutation_matrix<float, float>::transpose_from_left(unsigned int, unsigned int);
template void lean::permutation_matrix<float, float>::transpose_from_right(unsigned int, unsigned int);
template void lean::permutation_matrix<lean::mpq, lean::mpq>::apply_from_right(std::vector<lean::mpq, std::allocator<lean::mpq> >&);
template bool lean::permutation_matrix<lean::mpq, lean::mpq>::is_identity() const;
template void lean::permutation_matrix<lean::mpq, lean::mpq>::multiply_by_permutation_from_left(lean::permutation_matrix<lean::mpq, lean::mpq>&);
@ -44,8 +41,6 @@ template void lean::permutation_matrix<double, double>::apply_reverse_from_left<
template void lean::permutation_matrix<double, double>::apply_reverse_from_left<double>(std::vector<double, std::allocator<double> >&);
template void lean::permutation_matrix<double, double>::apply_reverse_from_right<double>(std::vector<double, std::allocator<double> >&);
template void lean::permutation_matrix<double, double>::transpose_from_right(unsigned int, unsigned int);
template void lean::permutation_matrix<float, float>::apply_from_left_perm<float>(lean::indexed_vector<float>&, lean::lp_settings&);
template void lean::permutation_matrix<float, float>::apply_from_left_perm<float>(std::vector<float, std::allocator<float> >&);
template void lean::permutation_matrix<lean::mpq, lean::mpq>::apply_from_left_perm<lean::mpq>(lean::indexed_vector<lean::mpq>&, lean::lp_settings&);
template void lean::permutation_matrix<lean::mpq, lean::mpq>::apply_from_left_perm<lean::mpq>(std::vector<lean::mpq, std::allocator<lean::mpq> >&);
template void lean::permutation_matrix<lean::mpq, lean::mpq>::apply_reverse_from_left<lean::mpq>(lean::indexed_vector<lean::mpq>&);

2
src/util/lp/row_eta_matrix.cpp
View file

@ -131,7 +131,7 @@ void row_eta_matrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & p
std::vector<unsigned> columns;
for (auto & it : m_row_vector.m_data)
columns.push_back(it.first);
for (unsigned i = columns.size(); i-- > 0;)
for (unsigned i = static_cast<unsigned>(columns.size()); i-- > 0;)
m_row_vector.m_data[i].first = p.get_rev(columns[i]);
#ifdef LEAN_DEBUG
// lean_assert(deb == *this);

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

@ -8,14 +8,6 @@
#pragma once
#include <vector>
#include "util/debug.h"
#include "util/numerics/numeric_traits.h"
#include "util/numerics/xnumeral.h"
#include "util/numerics/mpq.h"
#include "util/numerics/mpz.h"
#include "util/numerics/mpbq.h"
#include "util/numerics/double.h"
#include "util/numerics/float.h"
#include "util/numerics/mpfp.h"
#include <string>
#include "util/lp/sparse_vector.h"
#include "util/lp/indexed_vector.h"

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

@ -11,7 +11,7 @@
#include <algorithm>
#include <stdio.h> /* printf, fopen */
#include <stdlib.h> /* exit, EXIT_FAILURE */
#include "util/numerics/double.h"
#include "util/lp/lp_utils.h"
#include "util/lp/static_matrix.h"
namespace lean {
// for scaling an LP

40
src/util/lp/sparse_matrix.cpp
View file

@ -11,13 +11,13 @@ namespace lean {
template <typename T, typename X>
void sparse_matrix<T, X>::copy_column_from_static_matrix(unsigned col, static_matrix<T, X> const &A, unsigned col_index_in_the_new_matrix) {
std::vector<column_cell<T>> const & A_col_vector = A.m_columns[col];
unsigned size = A_col_vector.size();
unsigned size = static_cast<unsigned>(A_col_vector.size());
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();
unsigned col_offset = static_cast<unsigned>(new_column_vector.size());
std::vector<indexed_value<T>> & row_vector = m_rows[col_cell.m_i];
unsigned row_offset = row_vector.size();
unsigned row_offset = static_cast<unsigned>(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));
}
@ -193,7 +193,7 @@ void sparse_matrix<T, X>::set_max_in_row(std::vector<indexed_value<T>> & row_val
return;
T max_val = numeric_traits<T>::zero();
int max_index = 0;
for (unsigned i = row_vals.size(); i-- > 0;) {
for (unsigned i = static_cast<unsigned>(row_vals.size()); i-- > 0;) {
T iabs = abs(row_vals[i].m_value);
if (iabs > max_val) {
max_val = iabs;
@ -324,7 +324,7 @@ bool sparse_matrix<T, X>::set_row_from_work_vector_and_clean_work_vector(unsigne
template <typename T, typename X>
void sparse_matrix<T, X>::remove_zero_elements_and_set_data_on_existing_elements(unsigned row) {
auto & row_vals = m_rows[row];
for (unsigned k = row_vals.size(); k-- > 0;) { // we cannot simply run the iterator since we are removing
for (unsigned k = static_cast<unsigned>(row_vals.size()); k-- > 0;) { // we cannot simply run the iterator since we are removing
// elements from row_vals
auto & row_el_iv = row_vals[k];
unsigned j = row_el_iv.m_index;
@ -342,7 +342,7 @@ void sparse_matrix<T, X>::remove_zero_elements_and_set_data_on_existing_elements
template <typename T, typename X>
void sparse_matrix<T, X>::remove_zero_elements_and_set_data_on_existing_elements_not_adjusted(unsigned row, indexed_vector<T> & work_vec, lp_settings & settings) {
auto & row_vals = m_rows[row];
for (unsigned k = row_vals.size(); k-- > 0;) { // we cannot simply run the iterator since we are removing
for (unsigned k = static_cast<unsigned>(row_vals.size()); k-- > 0;) { // we cannot simply run the iterator since we are removing
// elements from row_vals
auto & row_el_iv = row_vals[k];
unsigned j = row_el_iv.m_index;
@ -619,7 +619,7 @@ void sparse_matrix<T, X>::remove_elements_that_are_not_in_w_and_update_common_el
// --------------------------------
// column_vals represents the old column
auto & column_vals = m_columns[column_to_replace].m_values;
for (unsigned k = column_vals.size(); k-- > 0;) {
for (unsigned k = static_cast<unsigned>(column_vals.size()); k-- > 0;) {
indexed_value<T> & col_el_iv = column_vals[k];
unsigned i = col_el_iv.m_index;
T w_data_at_i = w[adjust_row_inverse(i)];
@ -647,8 +647,8 @@ template <typename T, typename X>
void sparse_matrix<T, X>::add_new_element(unsigned row, unsigned col, T val) {
auto & row_vals = m_rows[row];
auto & col_vals = m_columns[col].m_values;
unsigned row_el_offs = row_vals.size();
unsigned col_el_offs = col_vals.size();
unsigned row_el_offs = static_cast<unsigned>(row_vals.size());
unsigned col_el_offs = static_cast<unsigned>(col_vals.size());
row_vals.push_back(indexed_value<T>(val, col, col_el_offs));
col_vals.push_back(indexed_value<T>(val, row, row_el_offs));
}
@ -666,7 +666,7 @@ void sparse_matrix<T, X>::add_new_elements_of_w_and_clear_w(unsigned column_to_r
auto & row_chunk = m_rows[ai];
lean_assert(row_chunk.size() > 0);
if (abs(w_at_i) > abs(row_chunk[0].m_value))
put_max_index_to_0(row_chunk, row_chunk.size() - 1);
put_max_index_to_0(row_chunk, static_cast<unsigned>(row_chunk.size()) - 1);
}
w[i] = numeric_traits<T>::zero();
}
@ -690,7 +690,7 @@ unsigned sparse_matrix<T, X>::pivot_score(unsigned i, unsigned j) {
auto col_header = m_columns[j];
return get_row_values(i).size() * (col_header.m_values.size() - col_header.m_shortened_markovitz - 1);
return static_cast<unsigned>(get_row_values(i).size() * (col_header.m_values.size() - col_header.m_shortened_markovitz - 1));
}
template <typename T, typename X>
@ -698,10 +698,10 @@ void sparse_matrix<T, X>::enqueue_domain_into_pivot_queue() {
lean_assert(m_pivot_queue.size() == 0);
for (unsigned i = 0; i < dimension(); i++) {
auto & rh = m_rows[i];
unsigned rnz = rh.size();
unsigned rnz = static_cast<unsigned>(rh.size());
for (auto iv : rh) {
unsigned j = iv.m_index;
m_pivot_queue.enqueue(i, j, rnz * (m_columns[j].m_values.size() - 1));
m_pivot_queue.enqueue(i, j, rnz * (static_cast<unsigned>(m_columns[j].m_values.size()) - 1));
}
}
}
@ -735,7 +735,7 @@ void sparse_matrix<T, X>::recover_pivot_queue(std::vector<upair> & rejected_pivo
}
template <typename T, typename X>
int sparse_matrix<T, X>::elem_is_too_small(unsigned i, unsigned j, const T & c_partial_pivoting) {
int sparse_matrix<T, X>::elem_is_too_small(unsigned i, unsigned j, int c_partial_pivoting) {
auto & row_chunk = m_rows[i];
if (j == row_chunk[0].m_index) {
@ -779,12 +779,12 @@ void sparse_matrix<T, X>::update_active_pivots(unsigned row) {
unsigned arow = adjust_row(row);
for (auto & iv : m_rows[arow]) {
col_header & ch = m_columns[iv.m_index];
int cols = ch.m_values.size() - ch.m_shortened_markovitz - 1;
int cols = static_cast<int>(ch.m_values.size()) - ch.m_shortened_markovitz - 1;
lean_assert(cols >= 0);
for (auto &ivc : ch.m_values) {
unsigned i = ivc.m_index;
if (adjust_row_inverse(i) <= row) continue; // the i is not an active row
m_pivot_queue.enqueue(i, iv.m_index, m_rows[i].size()*cols);
m_pivot_queue.enqueue(i, iv.m_index, static_cast<unsigned>(m_rows[i].size())*cols);
}
}
}
@ -802,10 +802,10 @@ bool sparse_matrix<T, X>::shorten_active_matrix(unsigned row, eta_matrix<T, X> *
for (auto & it : eta_matrix->m_column_vector.m_data) {
unsigned row = adjust_row(it.first);
auto & row_values = m_rows[row];
unsigned rnz = row_values.size();
unsigned rnz = static_cast<unsigned>(row_values.size());
for (auto & iv : row_values) {
col_header& ch = m_columns[iv.m_index];
int cnz = ch.m_values.size() - ch.m_shortened_markovitz - 1;
int cnz = static_cast<int>(ch.m_values.size()) - ch.m_shortened_markovitz - 1;
lean_assert(cnz >= 0);
m_pivot_queue.enqueue(row, iv.m_index, rnz * cnz);
}
@ -903,7 +903,7 @@ bool sparse_matrix<T, X>::pivot_queue_is_correct_after_pivoting(int k) {
#endif
template <typename T, typename X>
bool sparse_matrix<T, X>::get_pivot_for_column(unsigned &i, unsigned &j, T const & c_partial_pivoting, unsigned k) {
bool sparse_matrix<T, X>::get_pivot_for_column(unsigned &i, unsigned &j, int c_partial_pivoting, unsigned k) {
std::vector<upair> pivots_candidates_that_are_too_small;
while (!m_pivot_queue.is_empty()) {
m_pivot_queue.dequeue(i, j);
@ -933,7 +933,7 @@ bool sparse_matrix<T, X>::get_pivot_for_column(unsigned &i, unsigned &j, T const
}
template <typename T, typename X>
bool sparse_matrix<T, X>::elem_is_too_small(std::vector<indexed_value<T>> & row_chunk, indexed_value<T> & iv, T const & c_partial_pivoting) {
bool sparse_matrix<T, X>::elem_is_too_small(std::vector<indexed_value<T>> & row_chunk, indexed_value<T> & iv, int c_partial_pivoting) {
if (&iv == &row_chunk[0]) {
return false; // the max element is at the head
}

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

@ -7,7 +7,6 @@
#pragma once
#include <vector>
#include "util/numerics/float.h"
#include "util/lp/permutation_matrix.h"
#include <unordered_map>
#include "util/lp/static_matrix.h"
@ -140,7 +139,7 @@ public:
unsigned dimension() const {return m_row_permutation.size();}
unsigned dimension() const {return static_cast<unsigned>(m_row_permutation.size());}
#ifdef LEAN_DEBUG
unsigned row_count() const {return dimension();}
@ -316,7 +315,7 @@ public:
void recover_pivot_queue(std::vector<upair> & rejected_pivots);
int elem_is_too_small(unsigned i, unsigned j, const T & c_partial_pivoting);
int elem_is_too_small(unsigned i, unsigned j, int c_partial_pivoting);
bool remove_row_from_active_pivots_and_shorten_columns(unsigned row);
@ -336,12 +335,12 @@ public:
bool pivot_queue_is_correct_after_pivoting(int k);
bool get_pivot_for_column(unsigned &i, unsigned &j, T const & c_partial_pivoting, unsigned k);
bool get_pivot_for_column(unsigned &i, unsigned &j, int c_partial_pivoting, unsigned k);
bool elem_is_too_small(std::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, int c_partial_pivoting);
unsigned number_of_non_zeroes_in_row(unsigned row) const {
return m_rows[row].size();
return static_cast<unsigned>(m_rows[row].size());
}
unsigned number_of_non_zeroes_in_column(unsigned col) const {

11
src/util/lp/sparse_matrix_instances.cpp
View file

@ -16,7 +16,7 @@ template void sparse_matrix<double, double>::fill_eta_matrix(unsigned int, eta_m
template const double & sparse_matrix<double, double>::get(unsigned int, unsigned int) const;
template unsigned sparse_matrix<double, double>::get_number_of_nonzeroes() const;
template unsigned sparse_matrix<double, double>::get_number_of_nonzeroes_below_row(unsigned int) const;
template bool sparse_matrix<double, double>::get_pivot_for_column(unsigned int&, unsigned int&, double const&, unsigned int);
template bool sparse_matrix<double, double>::get_pivot_for_column(unsigned int&, unsigned int&, int, unsigned int);
template unsigned sparse_matrix<double, double>::lowest_row_in_column(unsigned int);
template bool sparse_matrix<double, double>::pivot_row_to_row(unsigned int, double, unsigned int, lp_settings&);
template bool sparse_matrix<double, double>::pivot_with_eta(unsigned int, eta_matrix<double, double>*, lp_settings&);
@ -30,18 +30,13 @@ template bool sparse_matrix<double, double>::shorten_active_matrix(unsigned int,
template void sparse_matrix<double, double>::solve_y_U(std::vector<double>&) const;
template sparse_matrix<double, double>::sparse_matrix(static_matrix<double, double> const&, std::vector<unsigned int>&);
template sparse_matrix<double, double>::sparse_matrix(unsigned int);
template float const & sparse_matrix<float, float>::get(unsigned int, unsigned int) const;
template bool sparse_matrix<float, float>::pivot_row_to_row(unsigned int, float, unsigned int, lp_settings&);
template void sparse_matrix<float, float>::replace_column(unsigned int, indexed_vector<float>&, lp_settings&);
template void sparse_matrix<float, float>::set(unsigned int, unsigned int, float);
template sparse_matrix<float, float>::sparse_matrix(unsigned int);
template void sparse_matrix<mpq, mpq>::add_new_element(unsigned int, unsigned int, mpq);
template void sparse_matrix<mpq, mpq>::divide_row_by_constant(unsigned int, mpq&, lp_settings&);
template void sparse_matrix<mpq, mpq>::fill_eta_matrix(unsigned int, eta_matrix<mpq, mpq>**);
template mpq const & sparse_matrix<mpq, mpq>::get(unsigned int, unsigned int) const;
template unsigned sparse_matrix<mpq, mpq>::get_number_of_nonzeroes() const;
template unsigned sparse_matrix<mpq, mpq>::get_number_of_nonzeroes_below_row(unsigned int) const;
template bool sparse_matrix<mpq, mpq>::get_pivot_for_column(unsigned int&, unsigned int&, mpq const&, unsigned int);
template bool sparse_matrix<mpq, mpq>::get_pivot_for_column(unsigned int&, unsigned int&, int, unsigned int);
template unsigned sparse_matrix<mpq, mpq>::lowest_row_in_column(unsigned int);
template bool sparse_matrix<mpq, mpq>::pivot_with_eta(unsigned int, eta_matrix<mpq, mpq>*, lp_settings&);
template void sparse_matrix<mpq, mpq>::prepare_for_factorization();
@ -58,7 +53,7 @@ template void sparse_matrix<mpq, numeric_pair<mpq>>::fill_eta_matrix(unsigne
template const mpq & sparse_matrix<mpq, numeric_pair<mpq>>::get(unsigned int, unsigned int) const;
template unsigned sparse_matrix<mpq, numeric_pair<mpq>>::get_number_of_nonzeroes() const;
template unsigned sparse_matrix<mpq, numeric_pair<mpq>>::get_number_of_nonzeroes_below_row(unsigned int) const;
template bool sparse_matrix<mpq, numeric_pair<mpq>>::get_pivot_for_column(unsigned int&, unsigned int&, mpq const&, unsigned int);
template bool sparse_matrix<mpq, numeric_pair<mpq>>::get_pivot_for_column(unsigned int&, unsigned int&, int, unsigned int);
template unsigned sparse_matrix<mpq, numeric_pair<mpq>>::lowest_row_in_column(unsigned int);
template bool sparse_matrix<mpq, numeric_pair<mpq>>::pivot_with_eta(unsigned int, eta_matrix<mpq, numeric_pair<mpq> >*, lp_settings&);
template void sparse_matrix<mpq, numeric_pair<mpq>>::prepare_for_factorization();

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

@ -7,23 +7,16 @@
#pragma once
#include <vector>
#include <utility>
#include "util/debug.h"
#include "util/pair.h"
#include "util/numerics/numeric_traits.h"
#include "util/numerics/xnumeral.h"
#include "util/numerics/mpq.h"
#include "util/numerics/mpz.h"
#include "util/numerics/mpbq.h"
#include "util/numerics/double.h"
#include "util/numerics/float.h"
#include "util/numerics/mpfp.h"
#include "util/lp/lp_utils.h"
#include "util/lp/lp_settings.h"
namespace lean {
template <typename T>
class sparse_vector {
public:
std::vector<pair<unsigned, T>> m_data;
std::vector<std::pair<unsigned, T>> m_data;
void push_back(unsigned index, T val) {
m_data.emplace_back(index, val);
}

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

@ -7,7 +7,6 @@
#pragma once
#include <vector>
#include "util/numerics/float.h"
#include "util/lp/permutation_matrix.h"
#include <unordered_map>
#include "util/lp/static_matrix.h"

16
src/util/lp/static_matrix.cpp
View file

@ -153,19 +153,19 @@ template <typename T, typename X> void static_matrix<T, X>::set(unsigned row,
if (numeric_traits<T>::is_zero(val)) return;
lean_assert(row < row_count() && col < column_count());
#ifdef LEAN_DEBUG
pair<unsigned, unsigned> p(row, col);
std::pair<unsigned, unsigned> p(row, col);
lean_assert(m_domain.find(p) == m_domain.end());
m_domain.insert(p);
#endif
auto & r = m_rows[row];
unsigned offs_in_cols = m_columns[col].size();
m_columns[col].push_back(make_column_cell(row, r.size(), val));
unsigned offs_in_cols = static_cast<unsigned>(m_columns[col].size());
m_columns[col].push_back(make_column_cell(row, static_cast<unsigned>(r.size()), val));
r.push_back(make_row_cell(col, offs_in_cols, val));
}
template <typename T, typename X>
std::set<pair<unsigned, unsigned>> static_matrix<T, X>::get_domain() {
std::set<pair<unsigned, unsigned>> ret;
std::set<std::pair<unsigned, unsigned>> static_matrix<T, X>::get_domain() {
std::set<std::pair<unsigned, unsigned>> ret;
for (unsigned i = 0; i < m_rows.size(); i++) {
for (auto it : m_rows[i]) {
ret.insert(std::make_pair(i, it.m_j));
@ -255,15 +255,15 @@ template <typename T, typename X> void static_matrix<T, X>::check_consistency
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);
std::pair<unsigned, unsigned> p(i, t.m_j);
lean_assert(by_rows.find(p) == by_rows.end());
by_rows[p] = t.get_val();
}
}
std::unordered_map<pair<unsigned, unsigned>, T> by_cols;
std::unordered_map<std::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);
std::pair<unsigned, unsigned> p(t.m_i, i);
lean_assert(by_cols.find(p) == by_cols.end());
by_cols[p] = get_value_of_column_cell(t);
}

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

@ -10,14 +10,6 @@
#include <set>
#include <unordered_map>
#include <utility>
#include "util/numerics/numeric_traits.h"
#include "util/numerics/xnumeral.h"
#include "util/numerics/mpq.h"
#include "util/numerics/mpz.h"
#include "util/numerics/mpbq.h"
#include "util/numerics/double.h"
#include "util/numerics/float.h"
#include "util/numerics/mpfp.h"
#include "util/lp/sparse_vector.h"
#include "util/lp/indexed_vector.h"
#include "util/lp/permutation_matrix.h"
@ -53,7 +45,7 @@ class static_matrix
#endif
{
#ifdef LEAN_DEBUG
std::set<pair<unsigned, unsigned>> m_domain;
std::set<std::pair<unsigned, unsigned>> m_domain;
#endif
public:
typedef std::vector<row_cell<T>> row_strip;
@ -103,9 +95,9 @@ public:
void init_empty_matrix(unsigned m, unsigned n);
unsigned row_count() const { return m_rows.size(); }
unsigned row_count() const { return static_cast<unsigned>(m_rows.size()); }
unsigned column_count() const { return m_columns.size(); }
unsigned column_count() const { return static_cast<unsigned>(m_columns.size()); }
template <typename L>
L dot_product_with_row(unsigned row, const std::vector<L> & w);;
@ -144,7 +136,7 @@ public:
ref operator()(unsigned row, unsigned col) { return ref(*this, row, col); }
std::set<pair<unsigned, unsigned>> get_domain();
std::set<std::pair<unsigned, unsigned>> get_domain();
void copy_column_to_vector (unsigned j, indexed_vector<T> & v) const;

3
src/util/lp/static_matrix_instances.cpp
View file

@ -7,6 +7,7 @@
#include <vector>
#include <memory>
#include <set>
#include <utility>
#include "util/lp/static_matrix.cpp"
#include "util/lp/lp_core_solver_base.h"
#include "util/lp/lp_dual_core_solver.h"
@ -29,7 +30,7 @@ template void static_matrix<double, double>::copy_column_to_vector(unsigned int,
template void static_matrix<double, double>::divide_row_by_constant(unsigned int, double const&);
template double static_matrix<double, double>::dot_product_with_column(std::vector<double, std::allocator<double> > const&, unsigned int) const;
template double static_matrix<double, double>::get_balance() const;
template std::set<pair<unsigned, unsigned>> static_matrix<double, double>::get_domain();
template std::set<std::pair<unsigned, unsigned>> static_matrix<double, double>::get_domain();
template double static_matrix<double, double>::get_elem(unsigned int, unsigned int) const;
template double static_matrix<double, double>::get_max_abs_in_column(unsigned int) const;
template double static_matrix<double, double>::get_min_abs_in_column(unsigned int) const;

1
src/util/lp/tail_matrix.h
View file

@ -9,6 +9,7 @@
#include <vector>
#include "util/lp/indexed_vector.h"
#include "util/lp/matrix.h"
#include "util/lp/lp_settings.h"
// These matrices appear at the end of the list
namespace lean {