mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-05-02 21:37:02 +00:00
Code Activity

merge with master branch

Browse source 
Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
This commit is contained in:
Nikolaj Bjorner 2017-09-19 09:39:22 -07:00
commit 651587ce01
1602 changed files with 40496 additions and 27837 deletions
Download patch file Download diff file Expand all files Collapse all files

7
src/util/lp/CMakeLists.txt
View file

@ -3,13 +3,11 @@ z3_add_component(lp
lp_utils.cpp
binary_heap_priority_queue_instances.cpp
binary_heap_upair_queue_instances.cpp
bound_propagator.cpp
lp_bound_propagator.cpp
core_solver_pretty_printer_instances.cpp
dense_matrix_instances.cpp
eta_matrix_instances.cpp
indexed_vector_instances.cpp
int_solver.cpp
lar_solver_instances.cpp
lar_core_solver_instances.cpp
lp_core_solver_base_instances.cpp
lp_dual_core_solver_instances.cpp
@ -20,7 +18,6 @@ z3_add_component(lp
lp_solver_instances.cpp
lu_instances.cpp
matrix_instances.cpp
nra_solver.cpp
permutation_matrix_instances.cpp
quick_xplain.cpp
row_eta_matrix_instances.cpp
@ -31,8 +28,6 @@ z3_add_component(lp
random_updater_instances.cpp
COMPONENT_DEPENDENCIES
util
polynomial
nlsat
PYG_FILES
lp_params.pyg
)

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

@ -1,8 +1,23 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
#include "util/debug.h"
@ -22,7 +37,7 @@ class binary_heap_priority_queue {
void put_at(unsigned i, unsigned h);
void decrease_priority(unsigned o, T newPriority);
public:
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
bool is_consistent() const;
#endif
public:
@ -60,10 +75,10 @@ public:
/// return the first element of the queue and removes it from the queue
unsigned dequeue();
unsigned peek() const {
lp_assert(m_heap_size > 0);
SASSERT(m_heap_size > 0);
return m_heap[1];
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
void print(std::ostream & out);
#endif
};

46
src/util/lp/binary_heap_priority_queue.hpp
View file

@ -1,11 +1,26 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/vector.h"
#include "util/lp/binary_heap_priority_queue.h"
namespace lp {
// is is the child place in heap
// this is the child place in the heap
template <typename T> void binary_heap_priority_queue<T>::swap_with_parent(unsigned i) {
unsigned parent = m_heap[i >> 1];
put_at(i >> 1, m_heap[i]);
@ -29,12 +44,12 @@ template <typename T> void binary_heap_priority_queue<T>::decrease_priority(unsi
}
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T> bool binary_heap_priority_queue<T>::is_consistent() const {
for (int i = 0; i < m_heap_inverse.size(); i++) {
int i_index = m_heap_inverse[i];
lp_assert(i_index <= static_cast<int>(m_heap_size));
lp_assert(i_index == -1 || m_heap[i_index] == i);
SASSERT(i_index <= static_cast<int>(m_heap_size));
SASSERT(i_index == -1 || m_heap[i_index] == i);
}
for (unsigned i = 1; i < m_heap_size; i++) {
unsigned ch = i << 1;
@ -49,13 +64,14 @@ template <typename T> bool binary_heap_priority_queue<T>::is_consistent() const
return true;
}
#endif
template <typename T> void binary_heap_priority_queue<T>::remove(unsigned o) {
T priority_of_o = m_priorities[o];
int o_in_heap = m_heap_inverse[o];
if (o_in_heap == -1) {
return; // nothing to do
}
lp_assert(static_cast<unsigned>(o_in_heap) <= m_heap_size);
SASSERT(static_cast<unsigned>(o_in_heap) <= m_heap_size);
if (static_cast<unsigned>(o_in_heap) < m_heap_size) {
put_at(o_in_heap, m_heap[m_heap_size--]);
if (m_priorities[m_heap[o_in_heap]] > priority_of_o) {
@ -72,11 +88,11 @@ template <typename T> void binary_heap_priority_queue<T>::remove(unsigned o) {
}
}
} else {
lp_assert(static_cast<unsigned>(o_in_heap) == m_heap_size);
SASSERT(static_cast<unsigned>(o_in_heap) == m_heap_size);
m_heap_size--;
}
m_heap_inverse[o] = -1;
// lp_assert(is_consistent());
// SASSERT(is_consistent());
}
// n is the initial queue capacity.
// The capacity will be enlarged two times automatically if needed
@ -102,7 +118,7 @@ template <typename T> void binary_heap_priority_queue<T>::put_to_heap(unsigned i
template <typename T> void binary_heap_priority_queue<T>::enqueue_new(unsigned o, const T& priority) {
m_heap_size++;
int i = m_heap_size;
lp_assert(o < m_priorities.size());
SASSERT(o < m_priorities.size());
m_priorities[o] = priority;
put_at(i, o);
while (i > 1 && m_priorities[m_heap[i >> 1]] > priority) {
@ -134,7 +150,7 @@ template <typename T> void binary_heap_priority_queue<T>::change_priority_for_ex
/// return the first element of the queue and removes it from the queue
template <typename T> unsigned binary_heap_priority_queue<T>::dequeue_and_get_priority(T & priority) {
lp_assert(m_heap_size != 0);
SASSERT(m_heap_size != 0);
int ret = m_heap[1];
priority = m_priorities[ret];
put_the_last_at_the_top_and_fix_the_heap();
@ -168,13 +184,13 @@ 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() {
lp_assert(m_heap_size > 0);
SASSERT(m_heap_size > 0);
int ret = m_heap[1];
put_the_last_at_the_top_and_fix_the_heap();
m_heap_inverse[ret] = -1;
return ret;
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T> void binary_heap_priority_queue<T>::print(std::ostream & out) {
vector<int> index;
vector<T> prs;

23
src/util/lp/binary_heap_priority_queue_instances.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/numeric_pair.h"
#include "util/lp/binary_heap_priority_queue.hpp"
namespace lp {

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <unordered_set>
@ -38,7 +53,7 @@ public:
void enqueue(unsigned i, unsigned j, const T & priority);
void dequeue(unsigned & i, unsigned &j);
T get_priority(unsigned i, unsigned j) const;
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
bool pair_to_index_is_a_bijection() const;
bool available_spots_are_correct() const;
bool is_correct() const {

31
src/util/lp/binary_heap_upair_queue.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <set>
#include "util/lp/lp_utils.h"
@ -14,7 +29,7 @@ template <typename T> binary_heap_upair_queue<T>::binary_heap_upair_queue(unsign
template <typename T> unsigned
binary_heap_upair_queue<T>::dequeue_available_spot() {
lp_assert(m_available_spots.empty() == false);
SASSERT(m_available_spots.empty() == false);
unsigned ret = m_available_spots.back();
m_available_spots.pop_back();
return ret;
@ -54,7 +69,7 @@ template <typename T> void binary_heap_upair_queue<T>::enqueue(unsigned i, unsig
m_pairs.resize(new_size);
}
ij_index = dequeue_available_spot();
// lp_assert(ij_index<m_pairs.size() && ij_index_is_new(ij_index));
// SASSERT(ij_index<m_pairs.size() && ij_index_is_new(ij_index));
m_pairs[ij_index] = p;
m_pairs_to_index[p] = ij_index;
} else {
@ -64,7 +79,7 @@ template <typename T> void binary_heap_upair_queue<T>::enqueue(unsigned i, unsig
}
template <typename T> void binary_heap_upair_queue<T>::dequeue(unsigned & i, unsigned &j) {
lp_assert(!m_q.is_empty());
SASSERT(!m_q.is_empty());
unsigned ij_index = m_q.dequeue();
upair & p = m_pairs[ij_index];
i = p.first;
@ -81,7 +96,7 @@ template <typename T> T binary_heap_upair_queue<T>::get_priority(unsigned i, uns
return m_q.get_priority(it->second);
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T> bool binary_heap_upair_queue<T>::pair_to_index_is_a_bijection() const {
std::set<int> tmp;
for (auto p : m_pairs_to_index) {

23
src/util/lp/binary_heap_upair_queue_instances.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/binary_heap_upair_queue.hpp"
namespace lp {
template binary_heap_upair_queue<int>::binary_heap_upair_queue(unsigned int);

109
src/util/lp/bound_analyzer_on_row.h
View file

@ -1,14 +1,29 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
#include "util/lp/linear_combination_iterator.h"
#include "implied_bound.h"
#include "test_bound_analyzer.h"
#include "util/lp/implied_bound.h"
#include "util/lp/test_bound_analyzer.h"
#include <functional>
#include "util/lp/bound_propagator.h"
#include "util/lp/lp_bound_propagator.h"
// We have an equality : sum by j of row[j]*x[j] = rs
// We try to pin a var by pushing the total by using the variable bounds
// In a loop we drive the partial sum down, denoting the variables of this process by _u.
@ -16,9 +31,9 @@
namespace lp {
class bound_analyzer_on_row {
linear_combination_iterator<mpq> & m_it;
bound_propagator & m_bp;
lp_bound_propagator & m_bp;
unsigned m_row_or_term_index;
int m_column_of_u; // index of an unlimited from above monoid
// -1 means that such a value is not found, -2 means that at least two of such monoids were found
@ -31,7 +46,7 @@ public :
linear_combination_iterator<mpq> &it,
const numeric_pair<mpq>& rs,
unsigned row_or_term_index,
bound_propagator & bp
lp_bound_propagator & bp
)
:
m_it(it),
@ -45,7 +60,7 @@ public :
unsigned j;
void analyze() {
mpq a; unsigned j;
while (((m_column_of_l != -2) || (m_column_of_u != -2)) && m_it.next(a, j))
analyze_bound_on_var_on_coeff(j, a);
@ -91,11 +106,11 @@ public :
}
const impq & ub(unsigned j) const {
lp_assert(upper_bound_is_available(j));
SASSERT(upper_bound_is_available(j));
return m_bp.get_upper_bound(j);
}
const impq & lb(unsigned j) const {
lp_assert(low_bound_is_available(j));
SASSERT(low_bound_is_available(j));
return m_bp.get_low_bound(j);
}
@ -114,29 +129,29 @@ public :
}
return a * lb(j).x;
}
mpq monoid_max(const mpq & a, unsigned j, bool & strict) const {
if (is_pos(a)) {
strict = !is_zero(ub(j).y);
return a * ub(j).x;
}
strict = !is_zero(lb(j).y);
return a * lb(j).x;
}
const mpq & monoid_min_no_mult(bool a_is_pos, unsigned j, bool & strict) const {
if (!a_is_pos) {
strict = !is_zero(ub(j).y);
return ub(j).x;
}
strict = !is_zero(lb(j).y);
return lb(j).x;
}
mpq monoid_max(const mpq & a, unsigned j, bool & strict) const {
if (is_pos(a)) {
strict = !is_zero(ub(j).y);
return a * ub(j).x;
}
strict = !is_zero(lb(j).y);
return a * lb(j).x;
}
const mpq & monoid_min_no_mult(bool a_is_pos, unsigned j, bool & strict) const {
if (!a_is_pos) {
strict = !is_zero(ub(j).y);
return ub(j).x;
}
strict = !is_zero(lb(j).y);
return lb(j).x;
}
mpq monoid_min(const mpq & a, unsigned j, bool& strict) const {
if (is_neg(a)) {
strict = !is_zero(ub(j).y);
return a * ub(j).x;
}
strict = !is_zero(lb(j).y);
return a * lb(j).x;
}
@ -145,15 +160,15 @@ public :
if (is_neg(a)) {
return a * ub(j).x;
}
return a * lb(j).x;
}
void limit_all_monoids_from_above() {
int strict = 0;
mpq total;
lp_assert(is_zero(total));
SASSERT(is_zero(total));
m_it.reset();
mpq a; unsigned j;
while (m_it.next(a, j)) {
@ -166,7 +181,7 @@ public :
m_it.reset();
while (m_it.next(a, j)) {
bool str;
bool a_is_pos = is_pos(a);
bool a_is_pos = is_pos(a);
mpq bound = total / a + monoid_min_no_mult(a_is_pos, j, str);
if (a_is_pos) {
limit_j(j, bound, true, false, strict - static_cast<int>(str) > 0);
@ -180,7 +195,7 @@ public :
void limit_all_monoids_from_below() {
int strict = 0;
mpq total;
lp_assert(is_zero(total));
SASSERT(is_zero(total));
m_it.reset();
mpq a; unsigned j;
while (m_it.next(a, j)) {
@ -192,9 +207,9 @@ public :
m_it.reset();
while (m_it.next(a, j)) {
bool str;
bool a_is_pos = is_pos(a);
mpq bound = total / a + monoid_max_no_mult(a_is_pos, j, str);
bool astrict = strict - static_cast<int>(str) > 0;
bool a_is_pos = is_pos(a);
mpq bound = total / a + monoid_max_no_mult(a_is_pos, j, str);
bool astrict = strict - static_cast<int>(str) > 0;
if (a_is_pos) {
limit_j(j, bound, true, true, astrict);
}
@ -204,7 +219,7 @@ public :
}
}
void limit_monoid_u_from_below() {
// we are going to limit from below the monoid m_column_of_u,
// every other monoid is impossible to limit from below
@ -225,7 +240,7 @@ public :
}
bound /= u_coeff;
if (numeric_traits<impq>::is_pos(u_coeff)) {
limit_j(m_column_of_u, bound, true, true, strict);
} else {
@ -260,7 +275,7 @@ public :
limit_j(m_column_of_l, bound, false, true, strict);
}
}
// // it is the coefficent before the bounded column
// void provide_evidence(bool coeff_is_pos) {
// /*
@ -272,7 +287,7 @@ public :
// mpq a; unsigned j;
// while (it->next(a, j)) {
// if (be.m_j == j) continue;
// lp_assert(bound_is_available(j, is_neg(a) ? low_bound : !low_bound));
// SASSERT(bound_is_available(j, is_neg(a) ? low_bound : !low_bound));
// be.m_vector_of_bound_signatures.emplace_back(a, j, numeric_traits<impq>::
// is_neg(a)? low_bound: !low_bound);
// }
@ -284,27 +299,27 @@ public :
m_bp.try_add_bound(u, j, is_low_bound, coeff_before_j_is_pos, m_row_or_term_index, strict);
}
void advance_u(unsigned j) {
if (m_column_of_u == -1)
m_column_of_u = j;
else
m_column_of_u = -2;
}
void advance_l(unsigned j) {
if (m_column_of_l == -1)
m_column_of_l = j;
else
m_column_of_l = -2;
}
void analyze_bound_on_var_on_coeff(int j, const mpq &a) {
switch (m_bp.get_column_type(j)) {
case column_type::low_bound:
if (numeric_traits<mpq>::is_pos(a))
advance_u(j);
else
else
advance_l(j);
break;
case column_type::upper_bound:
@ -325,7 +340,7 @@ public :
static void analyze_row(linear_combination_iterator<mpq> &it,
const numeric_pair<mpq>& rs,
unsigned row_or_term_index,
bound_propagator & bp
lp_bound_propagator & bp
) {
bound_analyzer_on_row a(it, rs, row_or_term_index, bp);
a.analyze();

47
src/util/lp/bound_propagator.cpp
View file

@ -1,47 +0,0 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
#include "util/lp/lar_solver.h"
namespace lp {
bound_propagator::bound_propagator(lar_solver & ls):
m_lar_solver(ls) {}
column_type bound_propagator::get_column_type(unsigned j) const {
return m_lar_solver.m_mpq_lar_core_solver.m_column_types()[j];
}
const impq & bound_propagator::get_low_bound(unsigned j) const {
return m_lar_solver.m_mpq_lar_core_solver.m_r_low_bounds()[j];
}
const impq & bound_propagator::get_upper_bound(unsigned j) const {
return m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j];
}
void bound_propagator::try_add_bound(const mpq & v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) {
j = m_lar_solver.adjust_column_index_to_term_index(j);
lconstraint_kind kind = is_low? GE : LE;
if (strict)
kind = static_cast<lconstraint_kind>(kind / 2);
if (!bound_is_interesting(j, kind, v))
return;
unsigned k; // index to ibounds
if (is_low) {
if (try_get_val(m_improved_low_bounds, j, k)) {
auto & found_bound = m_ibounds[k];
if (v > found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict))
found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
} else {
m_improved_low_bounds[j] = m_ibounds.size();
m_ibounds.push_back(implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
}
} else { // the upper bound case
if (try_get_val(m_improved_upper_bounds, j, k)) {
auto & found_bound = m_ibounds[k];
if (v < found_bound.m_bound || (v == found_bound.m_bound && found_bound.m_strict == false && strict))
found_bound = implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
} else {
m_improved_upper_bounds[j] = m_ibounds.size();
m_ibounds.push_back(implied_bound(v, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
}
}
}
}

23
src/util/lp/breakpoint.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -100,11 +115,11 @@ public:
}
T get_low_bound() const {
lp_assert(m_low_bound_is_set);
SASSERT(m_low_bound_is_set);
return m_low_bound;
}
T get_upper_bound() const {
lp_assert(m_upper_bound_is_set);
SASSERT(m_upper_bound_is_set);
return m_upper_bound;
}
@ -156,7 +171,7 @@ public:
}
T get_fixed_value() const {
lp_assert(m_is_fixed);
SASSERT(m_is_fixed);
return m_fixed_value;
}

23
src/util/lp/column_namer.h
View file

@ -1,8 +1,23 @@
#pragma once
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <string>
#include "util/lp/linear_combination_iterator.h"
namespace lp {

23
src/util/lp/conversion_helper.h
View file

@ -1,8 +1,23 @@
/*
Copyright (c) 2013 Microsoft Corporation. All rights reserved.
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
Author: Lev Nachmanson
*/
#pragma once
namespace lp {
template <typename V>

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <limits>
#include <string>

29
src/util/lp/core_solver_pretty_printer.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <limits>
#include <string>
#include <algorithm>
@ -148,7 +163,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::adjust_
case column_type::free_column:
break;
default:
lp_assert(false);
SASSERT(false);
break;
}
}
@ -357,7 +372,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_g
unsigned width = m_column_widths[col];
string s = row[col];
int number_of_blanks = width - static_cast<unsigned>(s.size());
lp_assert(number_of_blanks >= 0);
SASSERT(number_of_blanks >= 0);
print_blanks(number_of_blanks, m_out);
m_out << s << ' ';
if (col < row.size() - 1) {
@ -368,7 +383,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_g
string rs = T_to_string(rst);
int nb = m_rs_width - static_cast<int>(rs.size());
lp_assert(nb >= 0);
SASSERT(nb >= 0);
print_blanks(nb + 1, m_out);
m_out << rs << std::endl;
}

23
src/util/lp/core_solver_pretty_printer_instances.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/numeric_pair.h"
#include "util/lp/core_solver_pretty_printer.hpp"
template lp::core_solver_pretty_printer<double, double>::core_solver_pretty_printer(lp::lp_core_solver_base<double, double> &, std::ostream & out);

27
src/util/lp/dense_matrix.h
View file

@ -1,9 +1,24 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
#include "util/vector.h"
#include "util/lp/matrix.h"
namespace lp {
@ -31,7 +46,7 @@ public:
dense_matrix(unsigned m, unsigned n);
dense_matrix operator*=(matrix<T, X> const & a) {
lp_assert(column_count() == a.row_count());
SASSERT(column_count() == a.row_count());
dense_matrix c(row_count(), a.column_count());
for (unsigned i = 0; i < row_count(); i++) {
for (unsigned j = 0; j < a.column_count(); j++) {

27
src/util/lp/dense_matrix.hpp
View file

@ -1,9 +1,24 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/lp_settings.h"
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
#include "util/vector.h"
#include "util/lp/numeric_pair.h"
#include "util/lp/dense_matrix.h"
@ -170,7 +185,7 @@ template <typename T, typename X> void dense_matrix<T, X>::multiply_row_by_const
template <typename T, typename X>
dense_matrix<T, X> operator* (matrix<T, X> & a, matrix<T, X> & b){
lp_assert(a.column_count() == b.row_count());
SASSERT(a.column_count() == b.row_count());
dense_matrix<T, X> ret(a.row_count(), b.column_count());
for (unsigned i = 0; i < ret.m_m; i++)
for (unsigned j = 0; j< ret.m_n; j++) {

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

@ -1,10 +1,25 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/lp_settings.h"
#include "util/lp/dense_matrix.hpp"
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
#include "util/vector.h"
template lp::dense_matrix<double, double> lp::operator*<double, double>(lp::matrix<double, double>&, lp::matrix<double, double>&);
template void lp::dense_matrix<double, double>::apply_from_left(vector<double> &);

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -13,20 +28,20 @@ namespace lp {
template <typename T, typename X>
class eta_matrix
: public tail_matrix<T, X> {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
unsigned m_length;
#endif
unsigned m_column_index;
public:
sparse_vector<T> m_column_vector;
T m_diagonal_element;
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
eta_matrix(unsigned column_index, unsigned length):
#else
eta_matrix(unsigned column_index):
#endif
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
m_length(length),
#endif
m_column_index(column_index) {}
@ -61,7 +76,7 @@ public:
void push_back(unsigned row_index, T val ) {
lp_assert(row_index != m_column_index);
SASSERT(row_index != m_column_index);
m_column_vector.push_back(row_index, val);
}
@ -69,7 +84,7 @@ public:
void apply_from_right(indexed_vector<T> & w);
T get_elem(unsigned i, unsigned j) const;
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
unsigned row_count() const { return m_length; }
unsigned column_count() const { return m_length; }
void set_number_of_rows(unsigned m) { m_length = m; }

45
src/util/lp/eta_matrix.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -49,7 +64,7 @@ apply_from_left_local(indexed_vector<L> & w, lp_settings & settings) {
}
template <typename T, typename X>
void eta_matrix<T, X>::apply_from_right(vector<T> & w) {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// dense_matrix<T, X> deb(*this);
// auto clone_w = clone_vector<T>(w, get_number_of_rows());
// deb.apply_from_right(clone_w);
@ -59,8 +74,8 @@ void eta_matrix<T, X>::apply_from_right(vector<T> & w) {
t += w[it.first] * it.second;
}
w[m_column_index] = t;
#ifdef LEAN_DEBUG
// lp_assert(vectors_are_equal<T>(clone_w, w, get_number_of_rows()));
#ifdef Z3DEBUG
// SASSERT(vectors_are_equal<T>(clone_w, w, get_number_of_rows()));
// delete clone_w;
#endif
}
@ -68,7 +83,7 @@ template <typename T, typename X>
void eta_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
if (w.m_index.size() == 0)
return;
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// vector<T> wcopy(w.m_data);
// apply_from_right(wcopy);
#endif
@ -99,12 +114,12 @@ void eta_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
}
}
#ifdef LEAN_DEBUG
// lp_assert(w.is_OK());
// lp_assert(vectors_are_equal<T>(wcopy, w.m_data));
#ifdef Z3DEBUG
// SASSERT(w.is_OK());
// SASSERT(vectors_are_equal<T>(wcopy, w.m_data));
#endif
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X>
T eta_matrix<T, X>::get_elem(unsigned i, unsigned j) const {
if (j == m_column_index){
@ -120,7 +135,7 @@ T eta_matrix<T, X>::get_elem(unsigned i, unsigned j) const {
template <typename T, typename X>
void eta_matrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & p) {
// this = p * this * p(-1)
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// auto rev = p.get_reverse();
// auto deb = ((*this) * rev);
// deb = p * deb;
@ -129,8 +144,8 @@ void eta_matrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & p) {
for (auto & pair : m_column_vector.m_data) {
pair.first = p.get_rev(pair.first);
}
#ifdef LEAN_DEBUG
// lp_assert(deb == *this);
#ifdef Z3DEBUG
// SASSERT(deb == *this);
#endif
}
}

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

@ -1,12 +1,27 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <memory>
#include "util/vector.h"
#include "util/lp/numeric_pair.h"
#include "util/lp/eta_matrix.hpp"
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template double lp::eta_matrix<double, double>::get_elem(unsigned int, unsigned int) const;
template lp::mpq lp::eta_matrix<lp::mpq, lp::mpq>::get_elem(unsigned int, unsigned int) const;
template lp::mpq lp::eta_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >::get_elem(unsigned int, unsigned int) const;

23
src/util/lp/hash_helper.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <utility>
#include <functional>

23
src/util/lp/implied_bound.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/lp_settings.h"
#include "util/lp/lar_constraints.h"

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
namespace lp {
@ -41,7 +56,7 @@ public:
m_value = val;
}
};
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename X>
bool check_vector_for_small_values(indexed_vector<X> & w, lp_settings & settings) {
for (unsigned i : w.m_index) {

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -75,7 +90,16 @@ public:
}
void set_value(const T& value, unsigned index);
void set_value_as_in_dictionary(unsigned index) {
SASSERT(index < m_data.size());
T & loc = m_data[index];
if (is_zero(loc)) {
m_index.push_back(index);
loc = one_of_type<T>(); // use as a characteristic function
}
}
void clear();
void clear_all();
const T& operator[] (unsigned i) const {
@ -152,7 +176,7 @@ public:
}
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
bool is_OK() const;
void print(std::ostream & out);
#endif

29
src/util/lp/indexed_vector.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/vector.h"
#include "util/lp/indexed_vector.h"
#include "util/lp/lp_settings.h"
@ -41,13 +56,13 @@ template <typename T>
void indexed_vector<T>::resize(unsigned data_size) {
clear();
m_data.resize(data_size, numeric_traits<T>::zero());
lp_assert(is_OK());
SASSERT(is_OK());
}
template <typename T>
void indexed_vector<T>::set_value(const T& value, unsigned index) {
m_data[index] = value;
lp_assert(std::find(m_index.begin(), m_index.end(), index) == m_index.end());
SASSERT(std::find(m_index.begin(), m_index.end(), index) == m_index.end());
m_index.push_back(index);
}
@ -70,7 +85,7 @@ void indexed_vector<T>::erase_from_index(unsigned j) {
m_index.erase(it);
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T>
bool indexed_vector<T>::is_OK() const {
return true;

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/vector.h"
#include "util/lp/indexed_vector.hpp"
namespace lp {
@ -17,7 +32,7 @@ template void indexed_vector<mpq>::resize(unsigned int);
template void indexed_vector<unsigned>::resize(unsigned int);
template void indexed_vector<mpq>::set_value(const mpq&, unsigned int);
template void indexed_vector<unsigned>::set_value(const unsigned&, unsigned int);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template bool indexed_vector<double>::is_OK() const;
template bool indexed_vector<mpq>::is_OK() const;
template bool indexed_vector<lp::numeric_pair<mpq> >::is_OK() const;

25
src/util/lp/int_set.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
#include "util/lp/indexed_vector.h"
@ -20,7 +35,7 @@ public:
return m_data[j] >= 0;
}
void insert(unsigned j) {
lp_assert(j < m_data.size());
SASSERT(j < m_data.size());
if (contains(j)) return;
m_data[j] = m_index.size();
m_index.push_back(j);

1007
src/util/lp/int_solver.cpp
View file

File diff suppressed because it is too large Load diff

62
src/util/lp/int_solver.h
View file

@ -8,6 +8,9 @@
#include "util/lp/iterator_on_row.h"
#include "util/lp/int_set.h"
#include "util/lp/lar_term.h"
#include "util/lp/cut_solver.h"
#include "util/lp/lar_constraints.h"
namespace lp {
class lar_solver;
template <typename T, typename X>
@ -35,9 +38,6 @@ public:
int_set m_old_values_set;
vector<impq> m_old_values_data;
unsigned m_branch_cut_counter;
linear_combination_iterator<mpq>* m_iter_on_gomory_row;
unsigned m_gomory_cut_inf_column;
bool m_found_free_var_in_gomory_row;
// methods
int_solver(lar_solver* lp);
@ -46,6 +46,8 @@ public:
// main function to check that solution provided by lar_solver is valid for integral values,
// or provide a way of how it can be adjusted.
lia_move check(lar_term& t, mpq& k, explanation& ex);
bool move_non_basic_column_to_bounds(unsigned j);
lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex);
private:
// how to tighten bounds for integer variables.
@ -77,22 +79,22 @@ private:
explanation & ex);
void fill_explanation_from_fixed_columns(iterator_on_row<mpq> & it, explanation &);
void add_to_explanation_from_fixed_or_boxed_column(unsigned j, explanation &);
void remove_fixed_vars_from_base();
void patch_int_infeasible_columns();
void patch_int_infeasible_non_basic_column(unsigned j);
void patch_int_infeasible_nbasic_columns();
bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m);
linear_combination_iterator<mpq> * get_column_iterator(unsigned j);
bool lower(unsigned j) const;
bool upper(unsigned j) const;
const impq & lower_bound(unsigned j) const;
const impq & low_bound(unsigned j) const;
const impq & upper_bound(unsigned j) const;
bool is_int(unsigned j) const;
bool is_real(unsigned j) const;
bool is_base(unsigned j) const;
bool is_boxed(unsigned j) const;
bool is_fixed(unsigned j) const;
bool is_free(unsigned j) const;
bool value_is_int(unsigned j) const;
void set_value_for_nbasic_column(unsigned j, const impq & new_val);
void fix_non_base_columns();
void set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val);
bool non_basic_columns_are_at_bounds() const;
void failed();
bool is_feasible() const;
const impq & get_value(unsigned j) const;
@ -104,36 +106,52 @@ private:
int find_inf_int_base_column();
int find_inf_int_boxed_base_column_with_smallest_range();
lp_settings& settings();
void move_non_base_vars_to_bounds();
bool move_non_basic_columns_to_bounds();
void branch_infeasible_int_var(unsigned);
lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex);
lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex, unsigned inf_col, linear_combination_iterator<mpq>& iter);
lia_move report_conflict_from_gomory_cut(mpq & k);
lia_move report_gomory_cut(lar_term& t, mpq& k, mpq& lcm_den, unsigned num_ints);
void adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq& lcm_den);
void init_check_data();
bool constrain_free_vars(linear_combination_iterator<mpq> * r);
lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex);
int find_next_free_var_in_gomory_row();
bool is_gomory_cut_target();
lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex, unsigned j);
int find_free_var_in_gomory_row(linear_combination_iterator<mpq>& iter);
bool is_gomory_cut_target(linear_combination_iterator<mpq> &iter);
bool at_bound(unsigned j) const;
bool at_lower(unsigned j) const;
bool at_low(unsigned j) const;
bool at_upper(unsigned j) const;
bool has_low(unsigned j) const;
bool has_upper(unsigned j) const;
unsigned row_of_basic_column(unsigned j) const;
inline static bool is_rational(const impq & n) {
return is_zero(n.y);
}
inline static
mpq fractional_part(const impq & n) {
lp_assert(is_rational);
lp_assert(is_rational(n));
return n.x - floor(n.x);
}
void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex);
void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den);
void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation & ex, unsigned inf_column);
void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& t, explanation& ex, mpq & lcm_den, unsigned inf_column);
constraint_index column_upper_bound_constraint(unsigned j) const;
constraint_index column_low_bound_constraint(unsigned j) const;
void display_row_info(std::ostream & out, unsigned row_index) const;
void gomory_cut_adjust_t_and_k_for_size_1(const vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k);
void gomory_cut_adjust_t_and_k_for_size_gt_1(vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k, unsigned num_ints, mpq &lcm_den);
void gomory_cut_adjust_t_and_k(vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k, bool num_ints, mpq &lcm_den);
bool current_solution_is_inf_on_cut(const lar_term& t, const mpq& k) const;
public:
bool shift_var(unsigned j, unsigned range);
private:
unsigned random();
bool has_inf_int() const;
lia_move create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column) const;
public:
void display_inf_or_int_inf_columns(std::ostream & out) const;
template <typename T>
void fill_cut_solver(cut_solver<T> & cs);
template <typename T>
void fill_cut_solver_for_constraint(const lar_base_constraint*, cut_solver<T>& );
template <typename T>
void get_int_coeffs_from_constraint(const lar_base_constraint* c, vector<std::pair<T, var_index>>& coeff, T & rs);
};
}

23
src/util/lp/iterator_on_column.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/linear_combination_iterator.h"
#include "util/lp/static_matrix.h"

23
src/util/lp/iterator_on_indexed_vector.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/linear_combination_iterator.h"
namespace lp {

23
src/util/lp/iterator_on_pivot_row.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/iterator_on_indexed_vector.h"
namespace lp {

23
src/util/lp/iterator_on_row.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/linear_combination_iterator.h"
namespace lp {

23
src/util/lp/iterator_on_term_with_basis_var.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/linear_combination_iterator.h"
#include "util/lp/numeric_pair.h"

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -25,7 +40,7 @@ inline std::string lconstraint_kind_string(lconstraint_kind t) {
case GT: return std::string(">");
case EQ: return std::string("=");
}
lp_unreachable();
SASSERT(false);
return std::string(); // it is unreachable
}
@ -74,7 +89,7 @@ public:
: lar_base_constraint(kind, right_side), m_coeffs(left_side) {}
lar_constraint(const lar_base_constraint & c) {
lp_assert(false); // should not be called : todo!
SASSERT(false); // should not be called : todo!
}
unsigned size() const {

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
#include <string>
@ -168,9 +183,9 @@ public:
}
void push() {
lp_assert(m_r_solver.basis_heading_is_correct());
lp_assert(!need_to_presolve_with_double_solver() || m_d_solver.basis_heading_is_correct());
lp_assert(m_column_types.size() == m_r_A.column_count());
SASSERT(m_r_solver.basis_heading_is_correct());
SASSERT(!need_to_presolve_with_double_solver() || m_d_solver.basis_heading_is_correct());
SASSERT(m_column_types.size() == m_r_A.column_count());
m_stacked_simplex_strategy = settings().simplex_strategy();
m_stacked_simplex_strategy.push();
m_column_types.push();
@ -192,7 +207,7 @@ public:
template <typename K>
void push_vector(stacked_vector<K> & pushed_vector, const vector<K> & vector) {
lp_assert(pushed_vector.size() <= vector.size());
SASSERT(pushed_vector.size() <= vector.size());
for (unsigned i = 0; i < vector.size();i++) {
if (i == pushed_vector.size()) {
pushed_vector.push_back(vector[i]);
@ -242,8 +257,8 @@ public:
pop_basis(k);
m_stacked_simplex_strategy.pop(k);
settings().simplex_strategy() = m_stacked_simplex_strategy;
lp_assert(m_r_solver.basis_heading_is_correct());
lp_assert(!need_to_presolve_with_double_solver() || m_d_solver.basis_heading_is_correct());
SASSERT(m_r_solver.basis_heading_is_correct());
SASSERT(!need_to_presolve_with_double_solver() || m_d_solver.basis_heading_is_correct());
}
bool need_to_presolve_with_double_solver() const {
@ -304,11 +319,11 @@ public:
break;
default:
lp_assert(false);
SASSERT(false);
}
break;
default:
lp_unreachable();
SASSERT(false);
}
m_r_solver.remove_column_from_inf_set(j);
return true;
@ -317,7 +332,7 @@ public:
void prepare_solver_x_with_signature_tableau(const lar_solution_signature & signature) {
lp_assert(m_r_solver.inf_set_is_correct());
SASSERT(m_r_solver.inf_set_is_correct());
for (auto &t : signature) {
unsigned j = t.first;
if (m_r_heading[j] >= 0)
@ -332,9 +347,9 @@ public:
m_r_solver.m_x[jb] -= delta * m_r_solver.m_A.get_val(cc);
m_r_solver.update_column_in_inf_set(jb);
}
lp_assert(m_r_solver.A_mult_x_is_off() == false);
SASSERT(m_r_solver.A_mult_x_is_off() == false);
}
lp_assert(m_r_solver.inf_set_is_correct());
SASSERT(m_r_solver.inf_set_is_correct());
}
@ -342,7 +357,7 @@ public:
void prepare_solver_x_with_signature(const lar_solution_signature & signature, lp_primal_core_solver<L,K> & s) {
for (auto &t : signature) {
unsigned j = t.first;
lp_assert(m_r_heading[j] < 0);
SASSERT(m_r_heading[j] < 0);
auto pos_type = t.second;
switch (pos_type) {
case at_low_bound:
@ -359,7 +374,7 @@ public:
case not_at_bound:
switch (m_column_types[j]) {
case column_type::free_column:
lp_assert(false); // unreachable
SASSERT(false); // unreachable
case column_type::upper_bound:
s.m_x[j] = s.m_upper_bounds[j];
break;
@ -377,15 +392,15 @@ public:
s.m_x[j] = s.m_low_bounds[j];
break;
default:
lp_assert(false);
SASSERT(false);
}
break;
default:
lp_unreachable();
SASSERT(false);
}
}
lp_assert(is_zero_vector(s.m_b));
SASSERT(is_zero_vector(s.m_b));
s.solve_Ax_eq_b();
}
@ -418,7 +433,7 @@ public:
// the queues of delayed indices
std::queue<unsigned> entr_q, leav_q;
auto * l = cs.m_factorization;
lp_assert(l->get_status() == LU_status::OK);
SASSERT(l->get_status() == LU_status::OK);
for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) {
unsigned entering = trace_of_basis_change[i];
unsigned leaving = trace_of_basis_change[i+1];
@ -446,8 +461,8 @@ public:
continue;
}
}
lp_assert(cs.m_basis_heading[entering] < 0);
lp_assert(cs.m_basis_heading[leaving] >= 0);
SASSERT(cs.m_basis_heading[entering] < 0);
SASSERT(cs.m_basis_heading[leaving] >= 0);
if (l->get_status() == LU_status::OK) {
l->prepare_entering(entering, w); // to init vector w
l->replace_column(zero_of_type<L>(), w, cs.m_basis_heading[leaving]);
@ -471,7 +486,7 @@ public:
void solve_on_signature_tableau(const lar_solution_signature & signature, const vector<unsigned> & changes_of_basis) {
r_basis_is_OK();
lp_assert(settings().use_tableau());
SASSERT(settings().use_tableau());
bool r = catch_up_in_lu_tableau(changes_of_basis, m_d_solver.m_basis_heading);
if (!r) { // it is the case where m_d_solver gives a degenerated basis
@ -490,10 +505,10 @@ public:
return;
m_r_solver.stop_tracing_basis_changes();
// and now catch up in the double solver
lp_assert(m_r_solver.total_iterations() >= m_r_solver.m_trace_of_basis_change_vector.size() /2);
SASSERT(m_r_solver.total_iterations() >= m_r_solver.m_trace_of_basis_change_vector.size() /2);
catch_up_in_lu(m_r_solver.m_trace_of_basis_change_vector, m_r_solver.m_basis_heading, m_d_solver);
}
lp_assert(r_basis_is_OK());
SASSERT(r_basis_is_OK());
}
bool adjust_x_of_column(unsigned j) {
@ -507,16 +522,16 @@ public:
}
m_r_solver.snap_column_to_bound_tableau(j);
lp_assert(m_r_solver.column_is_feasible(j));
SASSERT(m_r_solver.column_is_feasible(j));
m_r_solver.m_inf_set.erase(j);
*/
lp_assert(false);
SASSERT(false);
return true;
}
bool catch_up_in_lu_tableau(const vector<unsigned> & trace_of_basis_change, const vector<int> & basis_heading) {
lp_assert(r_basis_is_OK());
SASSERT(r_basis_is_OK());
// the queues of delayed indices
std::queue<unsigned> entr_q, leav_q;
for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) {
@ -546,47 +561,47 @@ public:
continue;
}
}
lp_assert(m_r_solver.m_basis_heading[entering] < 0);
lp_assert(m_r_solver.m_basis_heading[leaving] >= 0);
SASSERT(m_r_solver.m_basis_heading[entering] < 0);
SASSERT(m_r_solver.m_basis_heading[leaving] >= 0);
m_r_solver.change_basis_unconditionally(entering, leaving);
if(!m_r_solver.pivot_column_tableau(entering, m_r_solver.m_basis_heading[entering])) {
// unroll the last step
// unroll the last step
m_r_solver.change_basis_unconditionally(leaving, entering);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
bool t =
#endif
m_r_solver.pivot_column_tableau(leaving, m_r_solver.m_basis_heading[leaving]);
#ifdef LEAN_DEBUG
lp_assert(t);
#ifdef Z3DEBUG
SASSERT(t);
#endif
return false;
}
}
lp_assert(r_basis_is_OK());
SASSERT(r_basis_is_OK());
return true;
}
bool r_basis_is_OK() const {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
if (!m_r_solver.m_settings.use_tableau())
return true;
for (unsigned j : m_r_solver.m_basis) {
lp_assert(m_r_solver.m_A.m_columns[j].size() == 1);
lp_assert(m_r_solver.m_A.get_val(m_r_solver.m_A.m_columns[j][0]) == one_of_type<mpq>());
SASSERT(m_r_solver.m_A.m_columns[j].size() == 1);
SASSERT(m_r_solver.m_A.get_val(m_r_solver.m_A.m_columns[j][0]) == one_of_type<mpq>());
}
for (unsigned j =0; j < m_r_solver.m_basis_heading.size(); j++) {
if (m_r_solver.m_basis_heading[j] >= 0) continue;
if (m_r_solver.m_column_types[j] == column_type::fixed) continue;
lp_assert(static_cast<unsigned>(- m_r_solver.m_basis_heading[j] - 1) < m_r_solver.m_column_types.size());
lp_assert( m_r_solver.m_basis_heading[j] <= -1);
SASSERT(static_cast<unsigned>(- m_r_solver.m_basis_heading[j] - 1) < m_r_solver.m_column_types.size());
SASSERT( m_r_solver.m_basis_heading[j] <= -1);
}
#endif
return true;
}
void solve_on_signature(const lar_solution_signature & signature, const vector<unsigned> & changes_of_basis) {
lp_assert(!settings().use_tableau());
SASSERT(!settings().use_tableau());
if (m_r_solver.m_factorization == nullptr) {
for (unsigned j = 0; j < changes_of_basis.size(); j+=2) {
unsigned entering = changes_of_basis[j];
@ -615,7 +630,7 @@ public:
return;
m_r_solver.stop_tracing_basis_changes();
// and now catch up in the double solver
lp_assert(m_r_solver.total_iterations() >= m_r_solver.m_trace_of_basis_change_vector.size() /2);
SASSERT(m_r_solver.total_iterations() >= m_r_solver.m_trace_of_basis_change_vector.size() /2);
catch_up_in_lu(m_r_solver.m_trace_of_basis_change_vector, m_r_solver.m_basis_heading, m_d_solver);
}
}
@ -641,7 +656,7 @@ public:
template <typename L, typename K>
void extract_signature_from_lp_core_solver(const lp_primal_core_solver<L, K> & solver, lar_solution_signature & signature) {
signature.clear();
lp_assert(signature.size() == 0);
SASSERT(signature.size() == 0);
for (unsigned j = 0; j < solver.m_basis_heading.size(); j++) {
if (solver.m_basis_heading[j] < 0) {
signature[j] = solver.get_non_basic_column_value_position(j);
@ -664,7 +679,7 @@ public:
if (upper_bound_is_set(j)) {
const auto & ub = m_r_solver.m_upper_bounds[j];
m_d_upper_bounds[j] = ub.x.get_double() + delta * ub.y.get_double();
lp_assert(!low_bound_is_set(j) || (m_d_upper_bounds[j] >= m_d_low_bounds[j]));
SASSERT(!low_bound_is_set(j) || (m_d_upper_bounds[j] >= m_d_low_bounds[j]));
}
}
}
@ -729,7 +744,7 @@ public:
case column_type::fixed:
return true;
default:
lp_assert(false);
SASSERT(false);
}
return false;
}
@ -744,20 +759,20 @@ public:
case column_type::fixed:
return true;
default:
lp_assert(false);
SASSERT(false);
}
return false;
}
void update_delta(mpq& delta, numeric_pair<mpq> const& l, numeric_pair<mpq> const& u) const {
lp_assert(l <= u);
SASSERT(l <= u);
if (l.x < u.x && l.y > u.y) {
mpq delta1 = (u.x - l.x) / (l.y - u.y);
if (delta1 < delta) {
delta = delta1;
}
}
lp_assert(l.x + delta * l.y <= u.x + delta * u.y);
SASSERT(l.x + delta * l.y <= u.x + delta * u.y);
}
@ -796,37 +811,6 @@ public:
return new iterator_on_indexed_vector<mpq>(m_r_solver.m_ed);
}
}
bool column_is_fixed(unsigned j) const {
return m_column_types()[j] == column_type::fixed ||
( m_column_types()[j] == column_type::boxed &&
m_r_solver.m_low_bounds[j] == m_r_solver.m_upper_bounds[j]);
}
const impq & low_bound(unsigned j) const {
lp_assert(m_column_types()[j] == column_type::fixed ||
m_column_types()[j] == column_type::boxed ||
m_column_types()[j] == column_type::low_bound);
return m_r_low_bounds[j];
}
const impq & upper_bound(unsigned j) const {
lp_assert(m_column_types()[j] == column_type::fixed ||
m_column_types()[j] == column_type::boxed ||
m_column_types()[j] == column_type::upper_bound);
return m_r_upper_bounds[j];
}
const bool column_is_bounded(unsigned j) const {
switch(m_column_types()[j]) {
case column_type::fixed:
case column_type::boxed:
return true;
default:
return false;
}
}
};
}

110
src/util/lp/lar_core_solver.hpp
View file

@ -1,11 +1,41 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <string>
#include "util/vector.h"
#include "util/lp/lar_core_solver.h"
@ -42,9 +72,9 @@ lar_core_solver::lar_core_solver(
column_names){}
void lar_core_solver::init_costs(bool first_time) {
lp_assert(false); // should not be called
// lp_assert(this->m_x.size() >= this->m_n());
// lp_assert(this->m_column_types.size() >= this->m_n());
SASSERT(false); // should not be called
// SASSERT(this->m_x.size() >= this->m_n());
// SASSERT(this->m_column_types.size() >= this->m_n());
// if (first_time)
// this->m_costs.resize(this->m_n());
// X inf = this->m_infeasibility;
@ -54,7 +84,7 @@ void lar_core_solver::init_costs(bool first_time) {
// if (!(first_time || inf >= this->m_infeasibility)) {
// LP_OUT(this->m_settings, "iter = " << this->total_iterations() << std::endl);
// LP_OUT(this->m_settings, "inf was " << T_to_string(inf) << " and now " << T_to_string(this->m_infeasibility) << std::endl);
// lp_assert(false);
// SASSERT(false);
// }
// if (inf == this->m_infeasibility)
// this->m_iters_with_no_cost_growing++;
@ -105,7 +135,7 @@ void lar_core_solver::init_cost_for_column(unsigned j) {
this->m_costs[j] = numeric_traits<T>::zero();
break;
default:
lp_assert(false);
SASSERT(false);
break;
}*/
}
@ -138,14 +168,30 @@ int lar_core_solver::column_is_out_of_bounds(unsigned j) {
return 0;
break;
}*/
lp_assert(false);
SASSERT(false);
return true;
}
void lar_core_solver::calculate_pivot_row(unsigned i) {
m_r_solver.calculate_pivot_row(i);
SASSERT(!m_r_solver.use_tableau());
SASSERT(m_r_solver.m_pivot_row.is_OK());
m_r_solver.m_pivot_row_of_B_1.clear();
m_r_solver.m_pivot_row_of_B_1.resize(m_r_solver.m_m());
m_r_solver.m_pivot_row.clear();
m_r_solver.m_pivot_row.resize(m_r_solver.m_n());
if (m_r_solver.m_settings.use_tableau()) {
unsigned basis_j = m_r_solver.m_basis[i];
for (auto & c : m_r_solver.m_A.m_rows[i]) {
if (c.m_j != basis_j)
m_r_solver.m_pivot_row.set_value(c.get_val(), c.m_j);
}
return;
}
m_r_solver.calculate_pivot_row_of_B_1(i);
m_r_solver.calculate_pivot_row_when_pivot_row_of_B1_is_ready(i);
}
@ -192,7 +238,7 @@ void lar_core_solver::calculate_pivot_row(unsigned i) {
}
void lar_core_solver::fill_not_improvable_zero_sum_from_inf_row() {
lp_assert(m_r_solver.A_mult_x_is_off() == false);
SASSERT(m_r_solver.A_mult_x_is_off() == false);
unsigned bj = m_r_basis[m_r_solver.m_inf_row_index_for_tableau];
m_infeasible_sum_sign = m_r_solver.inf_sign_of_column(bj);
m_infeasible_linear_combination.clear();
@ -227,32 +273,32 @@ void lar_core_solver::fill_not_improvable_zero_sum() {
void lar_core_solver::solve() {
lp_assert(m_r_solver.non_basic_columns_are_set_correctly());
lp_assert(m_r_solver.inf_set_is_correct());
SASSERT(m_r_solver.non_basic_columns_are_set_correctly());
SASSERT(m_r_solver.inf_set_is_correct());
if (m_r_solver.current_x_is_feasible() && m_r_solver.m_look_for_feasible_solution_only) {
m_r_solver.set_status(lp_status::OPTIMAL);
m_r_solver.set_status(OPTIMAL);
return;
}
++settings().st().m_need_to_solve_inf;
lp_assert(!m_r_solver.A_mult_x_is_off());
lp_assert((!settings().use_tableau()) || r_basis_is_OK());
SASSERT(!m_r_solver.A_mult_x_is_off());
SASSERT((!settings().use_tableau()) || r_basis_is_OK());
if (need_to_presolve_with_double_solver()) {
prefix_d();
lar_solution_signature solution_signature;
vector<unsigned> changes_of_basis = find_solution_signature_with_doubles(solution_signature);
if (m_d_solver.get_status() == lp_status::TIME_EXHAUSTED) {
m_r_solver.set_status(lp_status::TIME_EXHAUSTED);
if (m_d_solver.get_status() == TIME_EXHAUSTED) {
m_r_solver.set_status(TIME_EXHAUSTED);
return;
}
if (settings().use_tableau())
solve_on_signature_tableau(solution_signature, changes_of_basis);
else
solve_on_signature(solution_signature, changes_of_basis);
lp_assert(!settings().use_tableau() || r_basis_is_OK());
SASSERT(!settings().use_tableau() || r_basis_is_OK());
} else {
if (!settings().use_tableau()) {
bool snapped = m_r_solver.snap_non_basic_x_to_bound();
lp_assert(m_r_solver.non_basic_columns_are_set_correctly());
SASSERT(m_r_solver.non_basic_columns_are_set_correctly());
if (snapped)
m_r_solver.solve_Ax_eq_b();
}
@ -260,16 +306,16 @@ void lar_core_solver::solve() {
m_r_solver.find_feasible_solution();
else
m_r_solver.solve();
lp_assert(!settings().use_tableau() || r_basis_is_OK());
SASSERT(!settings().use_tableau() || r_basis_is_OK());
}
if (m_r_solver.get_status() == lp_status::INFEASIBLE) {
if (m_r_solver.get_status() == INFEASIBLE) {
fill_not_improvable_zero_sum();
} else if (m_r_solver.get_status() != lp_status::UNBOUNDED) {
m_r_solver.set_status(lp_status::OPTIMAL);
} else if (m_r_solver.get_status() != UNBOUNDED) {
m_r_solver.set_status(OPTIMAL);
}
lp_assert(r_basis_is_OK());
lp_assert(m_r_solver.non_basic_columns_are_set_correctly());
lp_assert(m_r_solver.inf_set_is_correct());
SASSERT(r_basis_is_OK());
SASSERT(m_r_solver.non_basic_columns_are_set_correctly());
SASSERT(m_r_solver.inf_set_is_correct());
}

23
src/util/lp/lar_core_solver_instances.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <utility>
#include <memory>
#include <string>

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"

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

File diff suppressed because it is too large Load diff

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

File diff suppressed because it is too large Load diff

54
src/util/lp/lar_term.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/indexed_vector.h"
namespace lp {
@ -10,7 +25,7 @@ struct lar_term {
std::unordered_map<unsigned, mpq> m_coeffs;
mpq m_v;
lar_term() {}
void add_monoid(const mpq& c, unsigned j) {
void add_to_map(unsigned j, const mpq& c) {
auto it = m_coeffs.find(j);
if (it == m_coeffs.end()) {
m_coeffs.emplace(j, c);
@ -21,10 +36,6 @@ struct lar_term {
}
}
bool is_empty() const {
return m_coeffs.size() == 0 && is_zero(m_v);
}
unsigned size() const { return static_cast<unsigned>(m_coeffs.size()); }
const std::unordered_map<unsigned, mpq> & coeffs() const {
@ -34,7 +45,7 @@ struct lar_term {
lar_term(const vector<std::pair<mpq, unsigned>>& coeffs,
const mpq & v) : m_v(v) {
for (const auto & p : coeffs) {
add_monoid(p.first, p.second);
add_to_map(p.second, p.first);
}
}
bool operator==(const lar_term & a) const { return false; } // take care not to create identical terms
@ -56,7 +67,7 @@ struct lar_term {
if (it == m_coeffs.end()) return;
const mpq & b = it->second;
for (unsigned it_j :li.m_index) {
add_monoid(- b * li.m_data[it_j], it_j);
add_to_map(it_j, - b * li.m_data[it_j]);
}
m_coeffs.erase(it);
}
@ -64,26 +75,5 @@ struct lar_term {
bool contains(unsigned j) const {
return m_coeffs.find(j) != m_coeffs.end();
}
void negate() {
for (auto & t : m_coeffs)
t.second.neg();
}
template <typename T>
T apply(const vector<T>& x) const {
T ret = T(m_v);
for (const auto & t : m_coeffs) {
ret += t.second * x[t.first];
}
return ret;
}
void clear() {
m_coeffs.clear();
m_v = zero_of_type<mpq>();
}
};
}

23
src/util/lp/linear_combination_iterator.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
namespace lp {
template <typename T>

67
src/util/lp/lp_bound_propagator.cpp Normal file
View file

@ -0,0 +1,67 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/lar_solver.h"
namespace lp {
lp_bound_propagator::lp_bound_propagator(lar_solver & ls):
m_lar_solver(ls) {}
column_type lp_bound_propagator::get_column_type(unsigned j) const {
return m_lar_solver.m_mpq_lar_core_solver.m_column_types()[j];
}
const impq & lp_bound_propagator::get_low_bound(unsigned j) const {
return m_lar_solver.m_mpq_lar_core_solver.m_r_low_bounds()[j];
}
const impq & lp_bound_propagator::get_upper_bound(unsigned j) const {
return m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j];
}
void lp_bound_propagator::try_add_bound(const mpq & v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) {
unsigned term_j = m_lar_solver.adjust_column_index_to_term_index(j);
mpq w = v;
if (term_j != j) {
j = term_j;
w += m_lar_solver.get_term(term_j).m_v; // when terms are turned into the columns they "lose" the right side, at this moment they aquire it back
}
lconstraint_kind kind = is_low? GE : LE;
if (strict)
kind = static_cast<lconstraint_kind>(kind / 2);
if (!bound_is_interesting(j, kind, w))
return;
unsigned k; // index to ibounds
if (is_low) {
if (try_get_val(m_improved_low_bounds, j, k)) {
auto & found_bound = m_ibounds[k];
if (w > found_bound.m_bound || (w == found_bound.m_bound && found_bound.m_strict == false && strict))
found_bound = implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
} else {
m_improved_low_bounds[j] = m_ibounds.size();
m_ibounds.push_back(implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
}
} else { // the upper bound case
if (try_get_val(m_improved_upper_bounds, j, k)) {
auto & found_bound = m_ibounds[k];
if (w < found_bound.m_bound || (w == found_bound.m_bound && found_bound.m_strict == false && strict))
found_bound = implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict);
} else {
m_improved_upper_bounds[j] = m_ibounds.size();
m_ibounds.push_back(implied_bound(w, j, is_low, coeff_before_j_is_pos, row_or_term_index, strict));
}
}
}
}

27
src/util/lp/bound_propagator.h → src/util/lp/lp_bound_propagator.h
View file

@ -1,19 +1,34 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/lp_settings.h"
namespace lp {
class lar_solver;
class bound_propagator {
class lp_bound_propagator {
std::unordered_map<unsigned, unsigned> m_improved_low_bounds; // these maps map a column index to the corresponding index in ibounds
std::unordered_map<unsigned, unsigned> m_improved_upper_bounds;
lar_solver & m_lar_solver;
public:
vector<implied_bound> m_ibounds;
public:
bound_propagator(lar_solver & ls);
lp_bound_propagator(lar_solver & ls);
column_type get_column_type(unsigned) const;
const impq & get_low_bound(unsigned) const;
const impq & get_upper_bound(unsigned) const;

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <set>
#include "util/vector.h"
@ -13,9 +28,6 @@
#include "util/lp/lu.h"
#include "util/lp/permutation_matrix.h"
#include "util/lp/column_namer.h"
#include "util/lp/iterator_on_row.h"
#include "util/lp/iterator_on_pivot_row.h"
namespace lp {
template <typename T, typename X> // X represents the type of the x variable and the bounds
@ -26,14 +38,7 @@ class lp_core_solver_base {
private:
lp_status m_status;
public:
bool current_x_is_feasible() const {
TRACE("feas",
if (m_inf_set.size()) {
tout << "column " << m_inf_set.m_index[0] << " is infeasible" << std::endl;
}
);
return m_inf_set.size() == 0;
}
bool current_x_is_feasible() const { return m_inf_set.size() == 0; }
bool current_x_is_infeasible() const { return m_inf_set.size() != 0; }
int_set m_inf_set;
bool m_using_infeas_costs;
@ -69,12 +74,6 @@ public:
bool m_tracing_basis_changes;
int_set* m_pivoted_rows;
bool m_look_for_feasible_solution_only;
std::function<void (unsigned, const X &)> * m_tracker_of_x_change;
void set_tracker_of_x(std::function<void (unsigned, const X&)>* tracker) {
m_tracker_of_x_change = tracker;
}
void start_tracing_basis_changes() {
m_trace_of_basis_change_vector.resize(0);
m_tracing_basis_changes = true;
@ -198,11 +197,11 @@ public:
bool need_to_pivot_to_basis_tableau() const {
lp_assert(m_A.is_correct());
SASSERT(m_A.is_correct());
unsigned m = m_A.row_count();
for (unsigned i = 0; i < m; i++) {
unsigned bj = m_basis[i];
lp_assert(m_A.m_columns[bj].size() > 0);
SASSERT(m_A.m_columns[bj].size() > 0);
if (m_A.m_columns[bj].size() > 1 || m_A.get_val(m_A.m_columns[bj][0]) != one_of_type<mpq>()) return true;
}
return false;
@ -211,7 +210,7 @@ public:
bool reduced_costs_are_correct_tableau() const {
if (m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows)
return true;
lp_assert(m_A.is_correct());
SASSERT(m_A.is_correct());
if (m_using_infeas_costs) {
if (infeasibility_costs_are_correct() == false) {
std::cout << "infeasibility_costs_are_correct() does not hold" << std::endl;
@ -386,11 +385,11 @@ public:
}
bool make_column_feasible(unsigned j, numeric_pair<mpq> & delta) {
lp_assert(m_basis_heading[j] < 0);
SASSERT(m_basis_heading[j] < 0);
auto & x = m_x[j];
switch (m_column_types[j]) {
case column_type::fixed:
lp_assert(m_low_bounds[j] == m_upper_bounds[j]);
SASSERT(m_low_bounds[j] == m_upper_bounds[j]);
if (x != m_low_bounds[j]) {
delta = m_low_bounds[j] - x;
x = m_low_bounds[j];
@ -426,7 +425,7 @@ public:
case column_type::free_column:
break;
default:
lp_assert(false);
SASSERT(false);
break;
}
return false;
@ -445,7 +444,6 @@ public:
void init_lu();
int pivots_in_column_and_row_are_different(int entering, int leaving) const;
void pivot_fixed_vars_from_basis();
bool pivot_column_general(unsigned j, unsigned j_basic, indexed_vector<T> & w);
bool pivot_for_tableau_on_basis();
bool pivot_row_for_tableau_on_basis(unsigned row);
void init_basic_part_of_basis_heading() {
@ -475,7 +473,7 @@ public:
}
void change_basis_unconditionally(unsigned entering, unsigned leaving) {
lp_assert(m_basis_heading[entering] < 0);
SASSERT(m_basis_heading[entering] < 0);
int place_in_non_basis = -1 - m_basis_heading[entering];
if (static_cast<unsigned>(place_in_non_basis) >= m_nbasis.size()) {
// entering variable in not in m_nbasis, we need to put it back;
@ -494,8 +492,7 @@ public:
}
void change_basis(unsigned entering, unsigned leaving) {
lp_assert(m_basis_heading[entering] < 0);
lp_assert(m_basis_heading[leaving] >= 0);
SASSERT(m_basis_heading[entering] < 0);
int place_in_basis = m_basis_heading[leaving];
int place_in_non_basis = - m_basis_heading[entering] - 1;
@ -536,7 +533,7 @@ public:
case column_type::free_column:
break;
default:
lp_assert(false);
SASSERT(false);
break;
}
return true;
@ -584,10 +581,10 @@ public:
case column_type::free_column:
break;
default:
lp_assert(false);
SASSERT(false);
}
out << "basis heading = " << m_basis_heading[j] << std::endl;
out << "x = " << m_x[j] << std::endl;
std::cout << "basis heading = " << m_basis_heading[j] << std::endl;
std::cout << "x = " << m_x[j] << std::endl;
/*
std::cout << "cost = " << m_costs[j] << std::endl;
std:: cout << "m_d = " << m_d[j] << std::endl;*/
@ -676,28 +673,24 @@ public:
void update_column_in_inf_set(unsigned j) {
if (column_is_feasible(j)) {
remove_column_from_inf_set(j);
m_inf_set.erase(j);
} else {
insert_column_into_inf_set(j);
m_inf_set.insert(j);
}
}
void insert_column_into_inf_set(unsigned j) {
if (m_tracker_of_x_change != nullptr)
(*m_tracker_of_x_change)(j, m_x[j]);
m_inf_set.insert(j);
lp_assert(!column_is_feasible(j));
SASSERT(!column_is_feasible(j));
}
void remove_column_from_inf_set(unsigned j) {
if (m_tracker_of_x_change != nullptr)
(*m_tracker_of_x_change)(j, m_x[j]);
m_inf_set.erase(j);
lp_assert(column_is_feasible(j));
SASSERT(column_is_feasible(j));
}
bool costs_on_nbasis_are_zeros() const {
lp_assert(this->basis_heading_is_correct());
SASSERT(this->basis_heading_is_correct());
for (unsigned j = 0; j < this->m_n(); j++) {
if (this->m_basis_heading[j] < 0)
lp_assert(is_zero(this->m_costs[j]));
SASSERT(is_zero(this->m_costs[j]));
}
return true;
}
@ -708,17 +701,5 @@ public:
const unsigned & iters_with_no_cost_growing() const {
return m_iters_with_no_cost_growing;
}
linear_combination_iterator<T> * get_iterator_on_row(unsigned i) {
if (m_settings.use_tableau())
return new iterator_on_row<T>(m_A.m_rows[i]);
calculate_pivot_row(i);
return new iterator_on_pivot_row<T>(m_pivot_row, m_basis[i]);
}
void calculate_pivot_row(unsigned i);
unsigned get_base_column_in_row(unsigned row_index) const {
return m_basis[row_index];
}
};
}

193
src/util/lp/lp_core_solver_base.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <set>
#include <string>
#include "util/vector.h"
@ -24,7 +39,7 @@ lp_core_solver_base(static_matrix<T, X> & A,
const vector<X> & upper_bound_values):
m_total_iterations(0),
m_iters_with_no_cost_growing(0),
m_status(lp_status::FEASIBLE),
m_status(FEASIBLE),
m_inf_set(A.column_count()),
m_using_infeas_costs(false),
m_pivot_row_of_B_1(A.row_count()),
@ -52,9 +67,8 @@ lp_core_solver_base(static_matrix<T, X> & A,
m_steepest_edge_coefficients(A.column_count()),
m_tracing_basis_changes(false),
m_pivoted_rows(nullptr),
m_look_for_feasible_solution_only(false),
m_tracker_of_x_change(nullptr) {
lp_assert(bounds_for_boxed_are_set_correctly());
m_look_for_feasible_solution_only(false) {
SASSERT(bounds_for_boxed_are_set_correctly());
init();
init_basis_heading_and_non_basic_columns_vector();
}
@ -62,7 +76,7 @@ lp_core_solver_base(static_matrix<T, X> & A,
template <typename T, typename X> void lp_core_solver_base<T, X>::
allocate_basis_heading() { // the rest of initilization will be handled by the factorization class
init_basis_heading_and_non_basic_columns_vector();
lp_assert(basis_heading_is_correct());
SASSERT(basis_heading_is_correct());
}
template <typename T, typename X> void lp_core_solver_base<T, X>::
init() {
@ -128,7 +142,7 @@ solve_yB(vector<T> & y) {
// }
// }
template <typename T, typename X> void lp_core_solver_base<T, X>::solve_Bd(unsigned entering, indexed_vector<T> & column) {
lp_assert(!m_settings.use_tableau());
SASSERT(!m_settings.use_tableau());
if (m_factorization == nullptr) {
init_factorization(m_factorization, m_A, m_basis, m_settings);
}
@ -138,19 +152,19 @@ template <typename T, typename X> void lp_core_solver_base<T, X>::solve_Bd(unsig
template <typename T, typename X> void lp_core_solver_base<T, X>::
solve_Bd(unsigned entering) {
lp_assert(m_ed.is_OK());
SASSERT(m_ed.is_OK());
m_factorization->solve_Bd(entering, m_ed, m_w);
if (this->precise())
m_columns_nz[entering] = m_ed.m_index.size();
lp_assert(m_ed.is_OK());
lp_assert(m_w.is_OK());
#ifdef LEAN_DEBUG
SASSERT(m_ed.is_OK());
SASSERT(m_w.is_OK());
#ifdef Z3DEBUG
// auto B = get_B(*m_factorization, m_basis);
// vector<T> a(m_m());
// m_A.copy_column_to_vector(entering, a);
// vector<T> cd(m_ed.m_data);
// B.apply_from_left(cd, m_settings);
// lp_assert(vectors_are_equal(cd , a));
// SASSERT(vectors_are_equal(cd , a));
#endif
}
@ -209,7 +223,7 @@ restore_m_ed(T * buffer) {
template <typename T, typename X> bool lp_core_solver_base<T, X>::
A_mult_x_is_off() const {
lp_assert(m_x.size() == m_A.column_count());
SASSERT(m_x.size() == m_A.column_count());
if (numeric_traits<T>::precise()) {
for (unsigned i = 0; i < m_m(); i++) {
X delta = m_b[i] - m_A.dot_product_with_row(i, m_x);
@ -245,7 +259,7 @@ A_mult_x_is_off() const {
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::
A_mult_x_is_off_on_index(const vector<unsigned> & index) const {
lp_assert(m_x.size() == m_A.column_count());
SASSERT(m_x.size() == m_A.column_count());
if (numeric_traits<T>::precise()) return false;
#if RUN_A_MULT_X_IS_OFF_FOR_PRECESE
for (unsigned i : index) {
@ -285,13 +299,13 @@ A_mult_x_is_off_on_index(const vector<unsigned> & index) const {
// from page 182 of Istvan Maros's book
template <typename T, typename X> void lp_core_solver_base<T, X>::
calculate_pivot_row_of_B_1(unsigned pivot_row) {
lp_assert(! use_tableau());
lp_assert(m_pivot_row_of_B_1.is_OK());
SASSERT(! use_tableau());
SASSERT(m_pivot_row_of_B_1.is_OK());
m_pivot_row_of_B_1.clear();
m_pivot_row_of_B_1.set_value(numeric_traits<T>::one(), pivot_row);
lp_assert(m_pivot_row_of_B_1.is_OK());
SASSERT(m_pivot_row_of_B_1.is_OK());
m_factorization->solve_yB_with_error_check_indexed(m_pivot_row_of_B_1, m_basis_heading, m_basis, m_settings);
lp_assert(m_pivot_row_of_B_1.is_OK());
SASSERT(m_pivot_row_of_B_1.is_OK());
}
@ -381,11 +395,11 @@ set_non_basic_x_to_correct_bounds() {
break;
case column_type::low_bound:
m_x[j] = m_low_bounds[j];
lp_assert(column_is_dual_feasible(j));
SASSERT(column_is_dual_feasible(j));
break;
case column_type::upper_bound:
m_x[j] = m_upper_bounds[j];
lp_assert(column_is_dual_feasible(j));
SASSERT(column_is_dual_feasible(j));
break;
default:
break;
@ -403,15 +417,15 @@ column_is_dual_feasible(unsigned j) const {
return x_is_at_low_bound(j) && d_is_not_negative(j);
case column_type::upper_bound:
LP_OUT(m_settings, "upper_bound type should be switched to low_bound" << std::endl);
lp_assert(false); // impossible case
SASSERT(false); // impossible case
case column_type::free_column:
return numeric_traits<X>::is_zero(m_d[j]);
default:
LP_OUT(m_settings, "column = " << j << std::endl);
LP_OUT(m_settings, "unexpected column type = " << column_type_to_string(m_column_types[j]) << std::endl);
lp_unreachable();
SASSERT(false);
}
lp_unreachable();
SASSERT(false);
return false;
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::
@ -494,7 +508,7 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::column_is_feas
return true;
break;
default:
lp_unreachable();
SASSERT(false);
}
return false; // it is unreachable
}
@ -535,7 +549,7 @@ update_basis_and_x(int entering, int leaving, X const & tt) {
if (!find_x_by_solving()) {
restore_x(entering, tt);
if(A_mult_x_is_off()) {
m_status = lp_status::FLOATING_POINT_ERROR;
m_status = FLOATING_POINT_ERROR;
m_iters_with_no_cost_growing++;
return false;
}
@ -545,7 +559,7 @@ update_basis_and_x(int entering, int leaving, X const & tt) {
if (m_factorization->get_status() != LU_status::OK) {
std::stringstream s;
// s << "failing refactor on off_result for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << total_iterations();
m_status = lp_status::FLOATING_POINT_ERROR;
m_status = FLOATING_POINT_ERROR;
return false;
}
return false;
@ -567,19 +581,19 @@ update_basis_and_x(int entering, int leaving, X const & tt) {
init_lu();
if (m_factorization->get_status() != LU_status::OK) {
if (m_look_for_feasible_solution_only && !precise()) {
m_status = lp_status::UNSTABLE;
m_status = UNSTABLE;
delete m_factorization;
m_factorization = nullptr;
return false;
}
// LP_OUT(m_settings, "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << total_iterations() << std::endl);
restore_x_and_refactor(entering, leaving, tt);
if (m_status == lp_status::FLOATING_POINT_ERROR)
if (m_status == FLOATING_POINT_ERROR)
return false;
lp_assert(!A_mult_x_is_off());
SASSERT(!A_mult_x_is_off());
m_iters_with_no_cost_growing++;
// LP_OUT(m_settings, "rolled back after failing of init_factorization()" << std::endl);
m_status = lp_status::UNSTABLE;
m_status = UNSTABLE;
return false;
}
return true;
@ -588,7 +602,7 @@ update_basis_and_x(int entering, int leaving, X const & tt) {
template <typename T, typename X> bool lp_core_solver_base<T, X>::
divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col) {
lp_assert(numeric_traits<T>::precise());
SASSERT(numeric_traits<T>::precise());
int pivot_index = -1;
auto & row = m_A.m_rows[pivot_row];
unsigned size = row.size();
@ -629,7 +643,7 @@ pivot_column_tableau(unsigned j, unsigned piv_row_index) {
return false;
if (pivot_col_cell_index != 0) {
lp_assert(column.size() > 1);
SASSERT(column.size() > 1);
// swap the pivot column cell with the head cell
auto c = column[0];
column[0] = column[pivot_col_cell_index];
@ -640,7 +654,7 @@ pivot_column_tableau(unsigned j, unsigned piv_row_index) {
}
while (column.size() > 1) {
auto & c = column.back();
lp_assert(c.m_i != piv_row_index);
SASSERT(c.m_i != piv_row_index);
if(! m_A.pivot_row_to_row_given_cell(piv_row_index, c, j)) {
return false;
}
@ -688,7 +702,7 @@ non_basis_is_correctly_represented_in_heading() const {
}
for (unsigned j = 0; j < m_A.column_count(); j++) {
if (m_basis_heading[j] >= 0) {
lp_assert(static_cast<unsigned>(m_basis_heading[j]) < m_A.row_count() && m_basis[m_basis_heading[j]] == j);
SASSERT(static_cast<unsigned>(m_basis_heading[j]) < m_A.row_count() && m_basis[m_basis_heading[j]] == j);
}
}
return true;
@ -696,9 +710,9 @@ non_basis_is_correctly_represented_in_heading() const {
template <typename T, typename X> bool lp_core_solver_base<T, X>::
basis_heading_is_correct() const {
lp_assert(m_basis_heading.size() == m_A.column_count());
lp_assert(m_basis.size() == m_A.row_count());
lp_assert(m_nbasis.size() <= m_A.column_count() - m_A.row_count()); // for the dual the size of non basis can be smaller
SASSERT(m_basis_heading.size() == m_A.column_count());
SASSERT(m_basis.size() == m_A.row_count());
SASSERT(m_nbasis.size() <= m_A.column_count() - m_A.row_count()); // for the dual the size of non basis can be smaller
if (!basis_has_no_doubles()) {
// std::cout << "basis_has_no_doubles" << std::endl;
return false;
@ -842,7 +856,7 @@ solve_Ax_eq_b() {
template <typename T, typename X> void lp_core_solver_base<T, X>::
snap_non_basic_x_to_bound_and_free_to_zeroes() {
for (unsigned j : non_basis()) {
lp_assert(j < m_x.size());
SASSERT(j < m_x.size());
switch (m_column_types[j]) {
case column_type::fixed:
case column_type::boxed:
@ -893,9 +907,9 @@ get_non_basic_column_value_position(unsigned j) const {
case column_type::upper_bound:
return x_is_at_upper_bound(j)? at_upper_bound : not_at_bound;
default:
lp_unreachable();
SASSERT(false);
}
lp_unreachable();
SASSERT(false);
return at_low_bound;
}
@ -924,55 +938,42 @@ template <typename T, typename X> void lp_core_solver_base<T, X>::transpose_row
transpose_basis(i, j);
m_A.transpose_rows(i, j);
}
// j is the new basic column, j_basic - the leaving column
template <typename T, typename X> bool lp_core_solver_base<T, X>::pivot_column_general(unsigned j, unsigned j_basic, indexed_vector<T> & w) {
lp_assert(m_basis_heading[j] < 0);
lp_assert(m_basis_heading[j_basic] >= 0);
unsigned row_index = m_basis_heading[j_basic];
if (m_settings.m_simplex_strategy == simplex_strategy_enum::lu) {
if (m_factorization->need_to_refactor()) {
init_lu();
}
else {
m_factorization->prepare_entering(j, w); // to init vector w
m_factorization->replace_column(zero_of_type<T>(), w, row_index);
}
if (m_factorization->get_status() != LU_status::OK) {
init_lu();
return false;
}
else {
change_basis(j, j_basic);
}
}
else { // the tableau case
if (pivot_column_tableau(j, row_index))
change_basis(j, j_basic);
else return false;
}
return true;
}
template <typename T, typename X> void lp_core_solver_base<T, X>::pivot_fixed_vars_from_basis() {
// run over basis and non-basis at the same time
indexed_vector<T> w(m_basis.size()); // the buffer
unsigned i = 0; // points to basis
for (; i < m_basis.size(); i++) {
unsigned basic_j = m_basis[i];
unsigned j = 0; // points to nonbasis
for (; i < m_basis.size() && j < m_nbasis.size(); i++) {
unsigned ii = m_basis[i];
unsigned jj;
if (get_column_type(basic_j) != column_type::fixed) continue;
T a;
unsigned j;
auto * it = get_iterator_on_row(i);
while (it->next(a, j)) {
if (j == basic_j)
continue;
if (get_column_type(j) != column_type::fixed) {
if (pivot_column_general(j, basic_j, w))
if (get_column_type(ii) != column_type::fixed) continue;
while (j < m_nbasis.size()) {
for (; j < m_nbasis.size(); j++) {
jj = m_nbasis[j];
if (get_column_type(jj) != column_type::fixed)
break;
}
if (j >= m_nbasis.size())
break;
j++;
if (m_factorization->need_to_refactor()) {
change_basis(jj, ii);
init_lu();
} else {
m_factorization->prepare_entering(jj, w); // to init vector w
m_factorization->replace_column(zero_of_type<T>(), w, m_basis_heading[ii]);
change_basis(jj, ii);
}
if (m_factorization->get_status() != LU_status::OK) {
change_basis(ii, jj);
init_lu();
} else {
break;
}
}
delete it;
SASSERT(m_factorization->get_status()== LU_status::OK);
}
}
@ -980,7 +981,7 @@ template <typename T, typename X> bool
lp_core_solver_base<T, X>::infeasibility_costs_are_correct() const {
if (! this->m_using_infeas_costs)
return true;
lp_assert(costs_on_nbasis_are_zeros());
SASSERT(costs_on_nbasis_are_zeros());
for (unsigned j :this->m_basis) {
if (!infeasibility_cost_is_correct_for_column(j)) {
std::cout << "infeasibility_cost_is_correct_for_column does not hold\n";
@ -1025,31 +1026,9 @@ lp_core_solver_base<T, X>::infeasibility_cost_is_correct_for_column(unsigned j)
case column_type::free_column:
return is_zero(this->m_costs[j]);
default:
lp_assert(false);
SASSERT(false);
return true;
}
}
template <typename T, typename X>
void lp_core_solver_base<T, X>::calculate_pivot_row(unsigned i) {
lp_assert(!use_tableau());
lp_assert(m_pivot_row.is_OK());
m_pivot_row_of_B_1.clear();
m_pivot_row_of_B_1.resize(m_m());
m_pivot_row.clear();
m_pivot_row.resize(m_n());
if (m_settings.use_tableau()) {
unsigned basis_j = m_basis[i];
for (auto & c : m_A.m_rows[i]) {
if (c.m_j != basis_j)
m_pivot_row.set_value(c.get_val(), c.m_j);
}
return;
}
calculate_pivot_row_of_B_1(i);
calculate_pivot_row_when_pivot_row_of_B1_is_ready(i);
}
}

24
src/util/lp/lp_core_solver_base_instances.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <utility>
#include <memory>
#include <string>
@ -129,4 +144,3 @@ template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::inf_set_is_correct() co
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::infeasibility_costs_are_correct() const;
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq >::infeasibility_costs_are_correct() const;
template bool lp::lp_core_solver_base<double, double >::infeasibility_costs_are_correct() const;
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::calculate_pivot_row(unsigned int);

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/static_matrix.h"
#include "util/lp/lp_core_solver_base.h"

113
src/util/lp/lp_dual_core_solver.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <algorithm>
#include <string>
#include "util/vector.h"
@ -23,7 +38,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::restore_non_ba
while (j--) {
if (this->m_basis_heading[j] >= 0 ) continue;
if (m_can_enter_basis[j]) {
lp_assert(std::find(nb.begin(), nb.end(), j) == nb.end());
SASSERT(std::find(nb.begin(), nb.end(), j) == nb.end());
nb.push_back(j);
this->m_basis_heading[j] = - static_cast<int>(nb.size());
}
@ -82,25 +97,25 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::start_with_ini
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::done() {
if (this->get_status() == lp_status::OPTIMAL) {
if (this->get_status() == OPTIMAL) {
return true;
}
if (this->total_iterations() > this->m_settings.max_total_number_of_iterations) { // debug !!!!
this->set_status(lp_status::ITERATIONS_EXHAUSTED);
this->set_status(ITERATIONS_EXHAUSTED);
return true;
}
return false; // todo, need to be more cases
}
template <typename T, typename X> T lp_dual_core_solver<T, X>::get_edge_steepness_for_low_bound(unsigned p) {
lp_assert(this->m_basis_heading[p] >= 0 && static_cast<unsigned>(this->m_basis_heading[p]) < this->m_m());
SASSERT(this->m_basis_heading[p] >= 0 && static_cast<unsigned>(this->m_basis_heading[p]) < this->m_m());
T del = this->m_x[p] - this->m_low_bounds[p];
del *= del;
return del / this->m_betas[this->m_basis_heading[p]];
}
template <typename T, typename X> T lp_dual_core_solver<T, X>::get_edge_steepness_for_upper_bound(unsigned p) {
lp_assert(this->m_basis_heading[p] >= 0 && static_cast<unsigned>(this->m_basis_heading[p]) < this->m_m());
SASSERT(this->m_basis_heading[p] >= 0 && static_cast<unsigned>(this->m_basis_heading[p]) < this->m_m());
T del = this->m_x[p] - this->m_upper_bounds[p];
del *= del;
return del / this->m_betas[this->m_basis_heading[p]];
@ -135,12 +150,12 @@ template <typename T, typename X> T lp_dual_core_solver<T, X>::pricing_for_row(u
return numeric_traits<T>::zero();
break;
case column_type::free_column:
lp_assert(numeric_traits<T>::is_zero(this->m_d[p]));
SASSERT(numeric_traits<T>::is_zero(this->m_d[p]));
return numeric_traits<T>::zero();
default:
lp_unreachable();
SASSERT(false);
}
lp_unreachable();
SASSERT(false);
return numeric_traits<T>::zero();
}
@ -170,8 +185,8 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::pricing_loop(u
}
} while (i != initial_offset_in_rows && rows_left);
if (m_r == -1) {
if (this->get_status() != lp_status::UNSTABLE) {
this->set_status(lp_status::OPTIMAL);
if (this->get_status() != UNSTABLE) {
this->set_status(OPTIMAL);
}
} else {
m_p = this->m_basis[m_r];
@ -181,10 +196,10 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::pricing_loop(u
return;
}
// failure in advance_on_known_p
if (this->get_status() == lp_status::FLOATING_POINT_ERROR) {
if (this->get_status() == FLOATING_POINT_ERROR) {
return;
}
this->set_status(lp_status::UNSTABLE);
this->set_status(UNSTABLE);
m_forbidden_rows.insert(m_r);
}
}
@ -209,9 +224,9 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::advance_on_kno
int pivot_compare_result = this->pivots_in_column_and_row_are_different(m_q, m_p);
if (!pivot_compare_result){;}
else if (pivot_compare_result == 2) { // the sign is changed, cannot continue
lp_unreachable(); // not implemented yet
SASSERT(false); // not implemented yet
} else {
lp_assert(pivot_compare_result == 1);
SASSERT(pivot_compare_result == 1);
this->init_lu();
}
DSE_FTran();
@ -228,21 +243,21 @@ template <typename T, typename X> int lp_dual_core_solver<T, X>::define_sign_of_
if (this->x_above_upper_bound(m_p)) {
return 1;
}
lp_unreachable();
SASSERT(false);
case column_type::low_bound:
if (this->x_below_low_bound(m_p)) {
return -1;
}
lp_unreachable();
SASSERT(false);
case column_type::upper_bound:
if (this->x_above_upper_bound(m_p)) {
return 1;
}
lp_unreachable();
SASSERT(false);
default:
lp_unreachable();
SASSERT(false);
}
lp_unreachable();
SASSERT(false);
return 0;
}
@ -250,10 +265,10 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::can_be_breakpo
if (this->pivot_row_element_is_too_small_for_ratio_test(j)) return false;
switch (this->m_column_types[j]) {
case column_type::low_bound:
lp_assert(this->m_settings.abs_val_is_smaller_than_harris_tolerance(this->m_x[j] - this->m_low_bounds[j]));
SASSERT(this->m_settings.abs_val_is_smaller_than_harris_tolerance(this->m_x[j] - this->m_low_bounds[j]));
return m_sign_of_alpha_r * this->m_pivot_row[j] > 0;
case column_type::upper_bound:
lp_assert(this->m_settings.abs_val_is_smaller_than_harris_tolerance(this->m_x[j] - this->m_upper_bounds[j]));
SASSERT(this->m_settings.abs_val_is_smaller_than_harris_tolerance(this->m_x[j] - this->m_upper_bounds[j]));
return m_sign_of_alpha_r * this->m_pivot_row[j] < 0;
case column_type::boxed:
{
@ -292,23 +307,23 @@ template <typename T, typename X> T lp_dual_core_solver<T, X>::get_delta() {
if (this->x_above_upper_bound(m_p)) {
return this->m_x[m_p] - this->m_upper_bounds[m_p];
}
lp_unreachable();
SASSERT(false);
case column_type::low_bound:
if (this->x_below_low_bound(m_p)) {
return this->m_x[m_p] - this->m_low_bounds[m_p];
}
lp_unreachable();
SASSERT(false);
case column_type::upper_bound:
if (this->x_above_upper_bound(m_p)) {
return get_edge_steepness_for_upper_bound(m_p);
}
lp_unreachable();
SASSERT(false);
case column_type::fixed:
return this->m_x[m_p] - this->m_upper_bounds[m_p];
default:
lp_unreachable();
SASSERT(false);
}
lp_unreachable();
SASSERT(false);
return zero_of_type<T>();
}
@ -355,7 +370,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::update_betas()
template <typename T, typename X> void lp_dual_core_solver<T, X>::apply_flips() {
for (unsigned j : m_flipped_boxed) {
lp_assert(this->x_is_at_bound(j));
SASSERT(this->x_is_at_bound(j));
if (this->x_is_at_low_bound(j)) {
this->m_x[j] = this->m_upper_bounds[j];
} else {
@ -385,7 +400,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::snap_xN_column
case column_type::free_column:
break;
default:
lp_unreachable();
SASSERT(false);
}
}
@ -441,7 +456,7 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::basis_change_a
return false;
}
lp_assert(d_is_correct());
SASSERT(d_is_correct());
return true;
}
@ -457,7 +472,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::recover_leavin
case free_of_bounds:
this->m_x[m_q] = zero_of_type<X>();
default:
lp_unreachable();
SASSERT(false);
}
}
@ -466,12 +481,12 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::revert_to_prev
this->change_basis_unconditionally(m_p, m_q);
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_settings);
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(lp_status::FLOATING_POINT_ERROR); // complete failure
this->set_status(FLOATING_POINT_ERROR); // complete failure
return;
}
recover_leaving();
if (!this->find_x_by_solving()) {
this->set_status(lp_status::FLOATING_POINT_ERROR);
this->set_status(FLOATING_POINT_ERROR);
return;
}
recalculate_xB_and_d();
@ -551,10 +566,10 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::delta_keeps_th
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::set_status_to_tentative_dual_unbounded_or_dual_unbounded() {
if (this->get_status() == lp_status::TENTATIVE_DUAL_UNBOUNDED) {
this->set_status(lp_status::DUAL_UNBOUNDED);
if (this->get_status() == TENTATIVE_DUAL_UNBOUNDED) {
this->set_status(DUAL_UNBOUNDED);
} else {
this->set_status(lp_status::TENTATIVE_DUAL_UNBOUNDED);
this->set_status(TENTATIVE_DUAL_UNBOUNDED);
}
}
@ -584,7 +599,7 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::tight_breakpoi
template <typename T, typename X> T lp_dual_core_solver<T, X>::calculate_harris_delta_on_breakpoint_set() {
bool first_time = true;
T ret = zero_of_type<T>();
lp_assert(m_breakpoint_set.size() > 0);
SASSERT(m_breakpoint_set.size() > 0);
for (auto j : m_breakpoint_set) {
T t;
if (this->x_is_at_low_bound(j)) {
@ -633,7 +648,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::find_q_on_tigh
}
}
m_tight_set.erase(m_q);
lp_assert(m_q != -1);
SASSERT(m_q != -1);
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::find_q_and_tight_set() {
@ -660,7 +675,7 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::ratio_test() {
set_status_to_tentative_dual_unbounded_or_dual_unbounded();
return false;
}
this->set_status(lp_status::FEASIBLE);
this->set_status(FEASIBLE);
find_q_and_tight_set();
if (!tight_breakpoinst_are_all_boxed()) break;
T del = m_delta - delta_lost_on_flips_of_tight_breakpoints() * initial_delta_sign;
@ -716,19 +731,19 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::update_xb_afte
template <typename T, typename X> void lp_dual_core_solver<T, X>::one_iteration() {
unsigned number_of_rows_to_try = get_number_of_rows_to_try_for_leaving();
unsigned offset_in_rows = this->m_settings.random_next() % this->m_m();
if (this->get_status() == lp_status::TENTATIVE_DUAL_UNBOUNDED) {
if (this->get_status() == TENTATIVE_DUAL_UNBOUNDED) {
number_of_rows_to_try = this->m_m();
} else {
this->set_status(lp_status::FEASIBLE);
this->set_status(FEASIBLE);
}
pricing_loop(number_of_rows_to_try, offset_in_rows);
lp_assert(problem_is_dual_feasible());
SASSERT(problem_is_dual_feasible());
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::solve() { // see the page 35
lp_assert(d_is_correct());
lp_assert(problem_is_dual_feasible());
lp_assert(this->basis_heading_is_correct());
SASSERT(d_is_correct());
SASSERT(problem_is_dual_feasible());
SASSERT(this->basis_heading_is_correct());
this->set_total_iterations(0);
this->iters_with_no_cost_growing() = 0;
do {
@ -736,7 +751,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::solve() { // s
return;
}
one_iteration();
} while (this->get_status() != lp_status::FLOATING_POINT_ERROR && this->get_status() != lp_status::DUAL_UNBOUNDED && this->get_status() != lp_status::OPTIMAL &&
} while (this->get_status() != FLOATING_POINT_ERROR && this->get_status() != DUAL_UNBOUNDED && this->get_status() != OPTIMAL &&
this->iters_with_no_cost_growing() <= this->m_settings.max_number_of_iterations_with_no_improvements
&& this->total_iterations() <= this->m_settings.max_total_number_of_iterations);
}

23
src/util/lp/lp_dual_core_solver_instances.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <utility>
#include <memory>
#include <string>

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
#include "util/lp/lp_utils.h"

95
src/util/lp/lp_dual_simplex.hpp
View file

@ -1,37 +1,52 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/lp_dual_simplex.h"
namespace lp{
template <typename T, typename X> void lp_dual_simplex<T, X>::decide_on_status_after_stage1() {
switch (m_core_solver->get_status()) {
case lp_status::OPTIMAL:
case OPTIMAL:
if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) {
this->m_status = lp_status::FEASIBLE;
this->m_status = FEASIBLE;
} else {
this->m_status = lp_status::UNBOUNDED;
this->m_status = UNBOUNDED;
}
break;
case lp_status::DUAL_UNBOUNDED:
lp_unreachable();
case lp_status::ITERATIONS_EXHAUSTED:
this->m_status = lp_status::ITERATIONS_EXHAUSTED;
case DUAL_UNBOUNDED:
SASSERT(false);
case ITERATIONS_EXHAUSTED:
this->m_status = ITERATIONS_EXHAUSTED;
break;
case lp_status::TIME_EXHAUSTED:
this->m_status = lp_status::TIME_EXHAUSTED;
case TIME_EXHAUSTED:
this->m_status = TIME_EXHAUSTED;
break;
case lp_status::FLOATING_POINT_ERROR:
this->m_status = lp_status::FLOATING_POINT_ERROR;
case FLOATING_POINT_ERROR:
this->m_status = FLOATING_POINT_ERROR;
break;
default:
lp_unreachable();
SASSERT(false);
}
}
template <typename T, typename X> void lp_dual_simplex<T, X>::fix_logical_for_stage2(unsigned j) {
lp_assert(j >= this->number_of_core_structurals());
SASSERT(j >= this->number_of_core_structurals());
switch (m_column_types_of_logicals[j - this->number_of_core_structurals()]) {
case column_type::low_bound:
m_low_bounds[j] = numeric_traits<T>::zero();
@ -44,7 +59,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fix_logical_for_st
m_can_enter_basis[j] = false;
break;
default:
lp_unreachable();
SASSERT(false);
}
}
@ -58,7 +73,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fix_structural_for
break;
case column_type::fixed:
case column_type::upper_bound:
lp_unreachable();
SASSERT(false);
case column_type::boxed:
this->m_upper_bounds[j] = ci->get_adjusted_upper_bound() / this->m_column_scale[j];
m_low_bounds[j] = numeric_traits<T>::zero();
@ -70,7 +85,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fix_structural_for
m_column_types_of_core_solver[j] = column_type::free_column;
break;
default:
lp_unreachable();
SASSERT(false);
}
// T cost_was = this->m_costs[j];
this->set_scaled_cost(j);
@ -99,23 +114,23 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::solve_for_stage2()
m_core_solver->solve_yB(m_core_solver->m_y);
m_core_solver->fill_reduced_costs_from_m_y_by_rows();
m_core_solver->start_with_initial_basis_and_make_it_dual_feasible();
m_core_solver->set_status(lp_status::FEASIBLE);
m_core_solver->set_status(FEASIBLE);
m_core_solver->solve();
switch (m_core_solver->get_status()) {
case lp_status::OPTIMAL:
this->m_status = lp_status::OPTIMAL;
case OPTIMAL:
this->m_status = OPTIMAL;
break;
case lp_status::DUAL_UNBOUNDED:
this->m_status = lp_status::INFEASIBLE;
case DUAL_UNBOUNDED:
this->m_status = INFEASIBLE;
break;
case lp_status::TIME_EXHAUSTED:
this->m_status = lp_status::TIME_EXHAUSTED;
case TIME_EXHAUSTED:
this->m_status = TIME_EXHAUSTED;
break;
case lp_status::FLOATING_POINT_ERROR:
this->m_status = lp_status::FLOATING_POINT_ERROR;
case FLOATING_POINT_ERROR:
this->m_status = FLOATING_POINT_ERROR;
break;
default:
lp_unreachable();
SASSERT(false);
}
this->m_second_stage_iterations = m_core_solver->total_iterations();
this->m_total_iterations = (this->m_first_stage_iterations + this->m_second_stage_iterations);
@ -129,7 +144,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fill_x_with_zeros(
}
template <typename T, typename X> void lp_dual_simplex<T, X>::stage1() {
lp_assert(m_core_solver == nullptr);
SASSERT(m_core_solver == nullptr);
this->m_x.resize(this->m_A->column_count(), numeric_traits<T>::zero());
if (this->m_settings.get_message_ostream() != nullptr)
this->print_statistics_on_A(*this->m_settings.get_message_ostream());
@ -151,7 +166,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::stage1() {
m_core_solver->start_with_initial_basis_and_make_it_dual_feasible();
if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) {
// skipping stage 1
m_core_solver->set_status(lp_status::OPTIMAL);
m_core_solver->set_status(OPTIMAL);
m_core_solver->set_total_iterations(0);
} else {
m_core_solver->solve();
@ -177,7 +192,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_s
}
template <typename T, typename X> column_type lp_dual_simplex<T, X>::get_column_type(unsigned j) {
lp_assert(j < this->m_A->column_count());
SASSERT(j < this->m_A->column_count());
if (j >= this->number_of_core_structurals()) {
return m_column_types_of_logicals[j - this->number_of_core_structurals()];
}
@ -186,12 +201,12 @@ template <typename T, typename X> column_type lp_dual_simplex<T, X>::get_column_
template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(unsigned j) {
// see 4.7 in the dissertation of Achim Koberstein
lp_assert(this->m_core_solver_columns_to_external_columns.find(j) !=
SASSERT(this->m_core_solver_columns_to_external_columns.find(j) !=
this->m_core_solver_columns_to_external_columns.end());
T free_bound = T(1e4); // see 4.8
unsigned jj = this->m_core_solver_columns_to_external_columns[j];
lp_assert(this->m_map_from_var_index_to_column_info.find(jj) != this->m_map_from_var_index_to_column_info.end());
SASSERT(this->m_map_from_var_index_to_column_info.find(jj) != this->m_map_from_var_index_to_column_info.end());
column_info<T> * ci = this->m_map_from_var_index_to_column_info[jj];
switch (ci->get_column_type()) {
case column_type::upper_bound: {
@ -221,14 +236,14 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_bounds_
this->m_upper_bounds[j] = this->m_low_bounds[j] = numeric_traits<T>::zero(); // is it needed?
break;
default:
lp_unreachable();
SASSERT(false);
}
m_column_types_of_core_solver[j] = column_type::boxed;
}
template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(unsigned j) {
this->m_costs[j] = 0;
lp_assert(get_column_type(j) != column_type::upper_bound);
SASSERT(get_column_type(j) != column_type::upper_bound);
if ((m_can_enter_basis[j] = (get_column_type(j) == column_type::low_bound))) {
m_column_types_of_core_solver[j] = column_type::boxed;
this->m_low_bounds[j] = numeric_traits<T>::zero();
@ -254,7 +269,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_and_bou
template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_solver_fields_for_row_slack_and_artificial(unsigned row,
unsigned & slack_var,
unsigned & artificial) {
lp_assert(row < this->row_count());
SASSERT(row < this->row_count());
auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[row]];
// we need to bring the program to the form Ax = b
T rs = this->m_b[row];
@ -336,7 +351,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::find_maximal_solut
this->flip_costs(); // do it for now, todo ( remove the flipping)
this->cleanup();
if (this->m_status == lp_status::INFEASIBLE) {
if (this->m_status == INFEASIBLE) {
return;
}
this->fill_matrix_A_and_init_right_side();
@ -346,7 +361,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::find_maximal_solut
fill_first_stage_solver_fields();
copy_m_b_aside_and_set_it_to_zeros();
stage1();
if (this->m_status == lp_status::FEASIBLE) {
if (this->m_status == FEASIBLE) {
stage2();
}
}

23
src/util/lp/lp_dual_simplex_instances.cpp
View file

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

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <list>
@ -70,7 +85,7 @@ public:
// unsigned len = 100000000;
// for (unsigned j : this->m_inf_set.m_index) {
// int i = this->m_basis_heading[j];
// lp_assert(i >= 0);
// SASSERT(i >= 0);
// unsigned row_len = this->m_A.m_rows[i].size();
// if (row_len < len) {
// choices.clear();
@ -98,8 +113,8 @@ public:
bool column_is_benefitial_for_entering_basis_on_sign_row_strategy(unsigned j, int sign) const {
// sign = 1 means the x of the basis column of the row has to grow to become feasible, when the coeff before j is neg, or x - has to diminish when the coeff is pos
// we have xbj = -aj * xj
lp_assert(this->m_basis_heading[j] < 0);
lp_assert(this->column_is_feasible(j));
SASSERT(this->m_basis_heading[j] < 0);
SASSERT(this->column_is_feasible(j));
switch (this->m_column_types[j]) {
case column_type::free_column: return true;
case column_type::fixed: return false;
@ -117,13 +132,13 @@ public:
return !this->x_is_at_upper_bound(j);
}
lp_assert(false); // cannot be here
SASSERT(false); // cannot be here
return false;
}
bool needs_to_grow(unsigned bj) const {
lp_assert(!this->column_is_feasible(bj));
SASSERT(!this->column_is_feasible(bj));
switch(this->m_column_types[bj]) {
case column_type::free_column:
return false;
@ -134,12 +149,12 @@ public:
default:
return false;
}
lp_assert(false); // unreachable
SASSERT(false); // unreachable
return false;
}
int inf_sign_of_column(unsigned bj) const {
lp_assert(!this->column_is_feasible(bj));
SASSERT(!this->column_is_feasible(bj));
switch(this->m_column_types[bj]) {
case column_type::free_column:
return 0;
@ -151,7 +166,7 @@ public:
default:
return -1;
}
lp_assert(false); // unreachable
SASSERT(false); // unreachable
return 0;
}
@ -159,7 +174,7 @@ public:
bool monoid_can_decrease(const row_cell<T> & rc) const {
unsigned j = rc.m_j;
lp_assert(this->column_is_feasible(j));
SASSERT(this->column_is_feasible(j));
switch (this->m_column_types[j]) {
case column_type::free_column:
return true;
@ -186,13 +201,13 @@ public:
default:
return false;
}
lp_assert(false); // unreachable
SASSERT(false); // unreachable
return false;
}
bool monoid_can_increase(const row_cell<T> & rc) const {
unsigned j = rc.m_j;
lp_assert(this->column_is_feasible(j));
SASSERT(this->column_is_feasible(j));
switch (this->m_column_types[j]) {
case column_type::free_column:
return true;
@ -219,7 +234,7 @@ public:
default:
return false;
}
lp_assert(false); // unreachable
SASSERT(false); // unreachable
return false;
}
@ -329,24 +344,24 @@ public:
}
void limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound(unsigned j, const T & m, X & theta, bool & unlimited) {
lp_assert(m < 0 && this->m_column_types[j] == column_type::upper_bound);
SASSERT(m < 0 && this->m_column_types[j] == column_type::upper_bound);
limit_inf_on_upper_bound_m_neg(m, this->m_x[j], this->m_upper_bounds[j], theta, unlimited);
}
void limit_theta_on_basis_column_for_inf_case_m_neg_low_bound(unsigned j, const T & m, X & theta, bool & unlimited) {
lp_assert(m < 0 && this->m_column_types[j] == column_type::low_bound);
SASSERT(m < 0 && this->m_column_types[j] == column_type::low_bound);
limit_inf_on_bound_m_neg(m, this->m_x[j], this->m_low_bounds[j], theta, unlimited);
}
void limit_theta_on_basis_column_for_inf_case_m_pos_low_bound(unsigned j, const T & m, X & theta, bool & unlimited) {
lp_assert(m > 0 && this->m_column_types[j] == column_type::low_bound);
SASSERT(m > 0 && this->m_column_types[j] == column_type::low_bound);
limit_inf_on_low_bound_m_pos(m, this->m_x[j], this->m_low_bounds[j], theta, unlimited);
}
void limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound(unsigned j, const T & m, X & theta, bool & unlimited) {
lp_assert(m > 0 && this->m_column_types[j] == column_type::upper_bound);
SASSERT(m > 0 && this->m_column_types[j] == column_type::upper_bound);
limit_inf_on_bound_m_pos(m, this->m_x[j], this->m_upper_bounds[j], theta, unlimited);
};
@ -359,7 +374,7 @@ public:
X get_max_bound(vector<X> & b);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
void check_Ax_equal_b();
void check_the_bounds();
void check_bound(unsigned i);
@ -388,7 +403,7 @@ public:
bool need_to_switch_costs() const {
if (this->m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows)
return false;
// lp_assert(calc_current_x_is_feasible() == current_x_is_feasible());
// SASSERT(calc_current_x_is_feasible() == current_x_is_feasible());
return this->current_x_is_feasible() == this->m_using_infeas_costs;
}
@ -443,7 +458,7 @@ public:
if (j == -1)
return -1;
lp_assert(!this->column_is_feasible(j));
SASSERT(!this->column_is_feasible(j));
switch (this->m_column_types[j]) {
case column_type::fixed:
case column_type::upper_bound:
@ -459,7 +474,7 @@ public:
new_val_for_leaving = this->m_low_bounds[j];
break;
default:
lp_assert(false);
SASSERT(false);
new_val_for_leaving = numeric_traits<T>::zero(); // does not matter
}
return j;
@ -469,7 +484,7 @@ public:
X new_val_for_leaving;
int leaving = find_leaving_tableau_rows(new_val_for_leaving);
if (leaving == -1) {
this->set_status(lp_status::OPTIMAL);
this->set_status(OPTIMAL);
return;
}
@ -485,14 +500,14 @@ public:
T a_ent;
int entering = find_beneficial_column_in_row_tableau_rows(this->m_basis_heading[leaving], a_ent);
if (entering == -1) {
this->set_status(lp_status::INFEASIBLE);
this->set_status(INFEASIBLE);
return;
}
X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent;
advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta );
lp_assert(this->m_x[leaving] == new_val_for_leaving);
SASSERT(this->m_x[leaving] == new_val_for_leaving);
if (this->current_x_is_feasible())
this->set_status(lp_status::OPTIMAL);
this->set_status(OPTIMAL);
}
void fill_breakpoints_array(unsigned entering);
@ -507,13 +522,13 @@ public:
void update_basis_and_x_with_comparison(unsigned entering, unsigned leaving, X delta);
void decide_on_status_when_cannot_find_entering() {
lp_assert(!need_to_switch_costs());
this->set_status(this->current_x_is_feasible()? lp_status::OPTIMAL: lp_status::INFEASIBLE);
SASSERT(!need_to_switch_costs());
this->set_status(this->current_x_is_feasible()? OPTIMAL: INFEASIBLE);
}
// void limit_theta_on_basis_column_for_feas_case_m_neg(unsigned j, const T & m, X & theta) {
// lp_assert(m < 0);
// lp_assert(this->m_column_type[j] == low_bound || this->m_column_type[j] == boxed);
// SASSERT(m < 0);
// SASSERT(this->m_column_type[j] == low_bound || this->m_column_type[j] == boxed);
// const X & eps = harris_eps_for_bound(this->m_low_bounds[j]);
// if (this->above_bound(this->m_x[j], this->m_low_bounds[j])) {
// theta = std::min((this->m_low_bounds[j] -this->m_x[j] - eps) / m, theta);
@ -522,7 +537,7 @@ public:
// }
void limit_theta_on_basis_column_for_feas_case_m_neg_no_check(unsigned j, const T & m, X & theta, bool & unlimited) {
lp_assert(m < 0);
SASSERT(m < 0);
const X& eps = harris_eps_for_bound(this->m_low_bounds[j]);
limit_theta((this->m_low_bounds[j] - this->m_x[j] - eps) / m, theta, unlimited);
if (theta < zero_of_type<X>()) theta = zero_of_type<X>();
@ -530,7 +545,7 @@ public:
bool limit_inf_on_bound_m_neg(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) {
// x gets smaller
lp_assert(m < 0);
SASSERT(m < 0);
if (numeric_traits<T>::precise()) {
if (this->below_bound(x, bound)) return false;
if (this->above_bound(x, bound)) {
@ -554,7 +569,7 @@ public:
bool limit_inf_on_bound_m_pos(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) {
// x gets larger
lp_assert(m > 0);
SASSERT(m > 0);
if (numeric_traits<T>::precise()) {
if (this->above_bound(x, bound)) return false;
if (this->below_bound(x, bound)) {
@ -579,14 +594,14 @@ public:
void limit_inf_on_low_bound_m_pos(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) {
if (numeric_traits<T>::precise()) {
// x gets larger
lp_assert(m > 0);
SASSERT(m > 0);
if (this->below_bound(x, bound)) {
limit_theta((bound - x) / m, theta, unlimited);
}
}
else {
// x gets larger
lp_assert(m > 0);
SASSERT(m > 0);
const X& eps = harris_eps_for_bound(bound);
if (this->below_bound(x, bound)) {
limit_theta((bound - x + eps) / m, theta, unlimited);
@ -596,7 +611,7 @@ public:
void limit_inf_on_upper_bound_m_neg(const T & m, const X & x, const X & bound, X & theta, bool & unlimited) {
// x gets smaller
lp_assert(m < 0);
SASSERT(m < 0);
const X& eps = harris_eps_for_bound(bound);
if (this->above_bound(x, bound)) {
limit_theta((bound - x - eps) / m, theta, unlimited);
@ -604,7 +619,7 @@ public:
}
void limit_theta_on_basis_column_for_inf_case_m_pos_boxed(unsigned j, const T & m, X & theta, bool & unlimited) {
// lp_assert(m > 0 && this->m_column_type[j] == column_type::boxed);
// SASSERT(m > 0 && this->m_column_type[j] == column_type::boxed);
const X & x = this->m_x[j];
const X & lbound = this->m_low_bounds[j];
@ -624,7 +639,7 @@ public:
}
void limit_theta_on_basis_column_for_inf_case_m_neg_boxed(unsigned j, const T & m, X & theta, bool & unlimited) {
// lp_assert(m < 0 && this->m_column_type[j] == column_type::boxed);
// SASSERT(m < 0 && this->m_column_type[j] == column_type::boxed);
const X & x = this->m_x[j];
const X & ubound = this->m_upper_bounds[j];
if (this->above_bound(x, ubound)) {
@ -642,7 +657,7 @@ public:
}
}
void limit_theta_on_basis_column_for_feas_case_m_pos(unsigned j, const T & m, X & theta, bool & unlimited) {
lp_assert(m > 0);
SASSERT(m > 0);
const T& eps = harris_eps_for_bound(this->m_upper_bounds[j]);
if (this->below_bound(this->m_x[j], this->m_upper_bounds[j])) {
limit_theta((this->m_upper_bounds[j] - this->m_x[j] + eps) / m, theta, unlimited);
@ -654,7 +669,7 @@ public:
}
void limit_theta_on_basis_column_for_feas_case_m_pos_no_check(unsigned j, const T & m, X & theta, bool & unlimited ) {
lp_assert(m > 0);
SASSERT(m > 0);
const X& eps = harris_eps_for_bound(this->m_upper_bounds[j]);
limit_theta( (this->m_upper_bounds[j] - this->m_x[j] + eps) / m, theta, unlimited);
if (theta < zero_of_type<X>()) {
@ -720,7 +735,7 @@ public:
break;
default:
lp_unreachable();
SASSERT(false);
}
if (!unlimited && theta < zero_of_type<X>()) {
theta = zero_of_type<X>();
@ -779,7 +794,7 @@ public:
if (this->m_basis_heading[j] < 0)
continue;
if (!this->column_is_feasible(j))
this->insert_column_into_inf_set(j);
this->m_inf_set.insert(j);
}
}
@ -803,7 +818,7 @@ public:
case column_type::free_column:
return 0;
default:
lp_assert(false);
SASSERT(false);
}
return 0;
}
@ -838,7 +853,7 @@ public:
return -1;
break;
default:
lp_assert(false);
SASSERT(false);
}
return 0;
@ -864,7 +879,7 @@ public:
// the delta is between the old and the new cost (old - new)
void update_reduced_cost_for_basic_column_cost_change(const T & delta, unsigned j) {
lp_assert(this->m_basis_heading[j] >= 0);
SASSERT(this->m_basis_heading[j] >= 0);
unsigned i = static_cast<unsigned>(this->m_basis_heading[j]);
for (const row_cell<T> & rc : this->m_A.m_rows[i]) {
unsigned k = rc.m_j;
@ -915,7 +930,7 @@ public:
} else {
m_converted_harris_eps = zero_of_type<T>();
}
this->set_status(lp_status::UNKNOWN);
this->set_status(UNKNOWN);
}
// constructor
@ -943,10 +958,10 @@ public:
upper_bound_values),
m_beta(A.row_count()),
m_converted_harris_eps(convert_struct<T, double>::convert(this->m_settings.harris_feasibility_tolerance)) {
lp_assert(initial_x_is_correct());
SASSERT(initial_x_is_correct());
m_low_bounds_dummy.resize(A.column_count(), zero_of_type<T>());
m_enter_price_eps = numeric_traits<T>::precise() ? numeric_traits<T>::zero() : T(1e-5);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// check_correctness();
#endif
}

231
src/util/lp/lp_primal_core_solver.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <list>
#include "util/vector.h"
#include <fstream>
@ -15,7 +30,7 @@ namespace lp {
template <typename T, typename X>
void lp_primal_core_solver<T, X>::sort_non_basis_rational() {
lp_assert(numeric_traits<T>::precise());
SASSERT(numeric_traits<T>::precise());
if (this->m_settings.use_tableau()) {
std::sort(this->m_nbasis.begin(), this->m_nbasis.end(), [this](unsigned a, unsigned b) {
unsigned ca = this->m_A.number_of_non_zeroes_in_column(a);
@ -70,11 +85,11 @@ bool lp_primal_core_solver<T, X>::column_is_benefitial_for_entering_on_breakpoin
const T & d = this->m_d[j];
switch (this->m_column_types[j]) {
case column_type::low_bound:
lp_assert(this->x_is_at_low_bound(j));
SASSERT(this->x_is_at_low_bound(j));
ret = d < -m_epsilon_of_reduced_cost;
break;
case column_type::upper_bound:
lp_assert(this->x_is_at_upper_bound(j));
SASSERT(this->x_is_at_upper_bound(j));
ret = d > m_epsilon_of_reduced_cost;
break;
case column_type::fixed:
@ -83,7 +98,7 @@ bool lp_primal_core_solver<T, X>::column_is_benefitial_for_entering_on_breakpoin
case column_type::boxed:
{
bool low_bound = this->x_is_at_low_bound(j);
lp_assert(low_bound || this->x_is_at_upper_bound(j));
SASSERT(low_bound || this->x_is_at_upper_bound(j));
ret = (low_bound && d < -m_epsilon_of_reduced_cost) || ((!low_bound) && d > m_epsilon_of_reduced_cost);
}
break;
@ -91,7 +106,7 @@ bool lp_primal_core_solver<T, X>::column_is_benefitial_for_entering_on_breakpoin
ret = d > m_epsilon_of_reduced_cost || d < - m_epsilon_of_reduced_cost;
break;
default:
lp_unreachable();
SASSERT(false);
ret = false;
break;
}
@ -127,14 +142,14 @@ bool lp_primal_core_solver<T, X>::column_is_benefitial_for_entering_basis(unsign
}
break;
default:
lp_unreachable();
SASSERT(false);
break;
}
return false;
}
template <typename T, typename X>
bool lp_primal_core_solver<T, X>::column_is_benefitial_for_entering_basis_precise(unsigned j) const {
lp_assert (numeric_traits<T>::precise());
SASSERT (numeric_traits<T>::precise());
if (this->m_using_infeas_costs && this->m_settings.use_breakpoints_in_feasibility_search)
return column_is_benefitial_for_entering_on_breakpoints(j);
const T& dj = this->m_d[j];
@ -167,7 +182,7 @@ bool lp_primal_core_solver<T, X>::column_is_benefitial_for_entering_basis_precis
}
break;
default:
lp_unreachable();
SASSERT(false);
break;
}
return false;
@ -175,7 +190,7 @@ bool lp_primal_core_solver<T, X>::column_is_benefitial_for_entering_basis_precis
template <typename T, typename X>
int lp_primal_core_solver<T, X>::choose_entering_column_presize(unsigned number_of_benefitial_columns_to_go_over) { // at this moment m_y = cB * B(-1)
lp_assert(numeric_traits<T>::precise());
SASSERT(numeric_traits<T>::precise());
if (number_of_benefitial_columns_to_go_over == 0)
return -1;
if (this->m_basis_sort_counter == 0) {
@ -259,7 +274,7 @@ int lp_primal_core_solver<T, X>::choose_entering_column(unsigned number_of_benef
template <typename T, typename X> int lp_primal_core_solver<T, X>::advance_on_sorted_breakpoints(unsigned entering, X &t) {
T slope_at_entering = this->m_d[entering];
breakpoint<X> * last_bp = nullptr;
lp_assert(m_breakpoint_indices_queue.is_empty()==false);
SASSERT(m_breakpoint_indices_queue.is_empty()==false);
while (m_breakpoint_indices_queue.is_empty() == false) {
unsigned bi = m_breakpoint_indices_queue.dequeue();
breakpoint<X> *b = &m_breakpoints[bi];
@ -274,7 +289,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::advance_on_so
}
}
}
lp_assert (last_bp != nullptr);
SASSERT (last_bp != nullptr);
t = last_bp->m_delta;
return last_bp->m_j;
}
@ -282,13 +297,13 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::advance_on_so
template <typename T, typename X> int
lp_primal_core_solver<T, X>::find_leaving_and_t_with_breakpoints(unsigned entering, X & t){
lp_assert(this->precise() == false);
SASSERT(this->precise() == false);
fill_breakpoints_array(entering);
return advance_on_sorted_breakpoints(entering, t);
}
template <typename T, typename X> bool lp_primal_core_solver<T, X>::get_harris_theta(X & theta) {
lp_assert(this->m_ed.is_OK());
SASSERT(this->m_ed.is_OK());
bool unlimited = true;
for (unsigned i : this->m_ed.m_index) {
if (this->m_settings.abs_val_is_smaller_than_pivot_tolerance(this->m_ed[i])) continue;
@ -345,13 +360,13 @@ template <typename T, typename X> bool lp_primal_core_solver<T, X>::try_jump_to_
if (m_sign_of_entering_delta > 0) {
t = this->m_upper_bounds[entering] - this->m_x[entering];
if (unlimited || t <= theta){
lp_assert(t >= zero_of_type<X>());
SASSERT(t >= zero_of_type<X>());
return true;
}
} else { // m_sign_of_entering_delta == -1
t = this->m_x[entering] - this->m_low_bounds[entering];
if (unlimited || t <= theta) {
lp_assert(t >= zero_of_type<X>());
SASSERT(t >= zero_of_type<X>());
return true;
}
}
@ -360,7 +375,7 @@ template <typename T, typename X> bool lp_primal_core_solver<T, X>::try_jump_to_
if (m_sign_of_entering_delta > 0) {
t = this->m_upper_bounds[entering] - this->m_x[entering];
if (unlimited || t <= theta){
lp_assert(t >= zero_of_type<X>());
SASSERT(t >= zero_of_type<X>());
return true;
}
}
@ -369,7 +384,7 @@ template <typename T, typename X> bool lp_primal_core_solver<T, X>::try_jump_to_
if (m_sign_of_entering_delta < 0) {
t = this->m_x[entering] - this->m_low_bounds[entering];
if (unlimited || t <= theta) {
lp_assert(t >= zero_of_type<X>());
SASSERT(t >= zero_of_type<X>());
return true;
}
}
@ -405,7 +420,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
do {
unsigned i = this->m_ed.m_index[k];
const T & ed = this->m_ed[i];
lp_assert(!numeric_traits<T>::is_zero(ed));
SASSERT(!numeric_traits<T>::is_zero(ed));
unsigned j = this->m_basis[i];
limit_theta_on_basis_column(j, - ed * m_sign_of_entering_delta, t, unlimited);
if (!unlimited) {
@ -424,7 +439,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
while (k != initial_k) {
unsigned i = this->m_ed.m_index[k];
const T & ed = this->m_ed[i];
lp_assert(!numeric_traits<T>::is_zero(ed));
SASSERT(!numeric_traits<T>::is_zero(ed));
unsigned j = this->m_basis[i];
unlimited = true;
limit_theta_on_basis_column(j, -ed * m_sign_of_entering_delta, ratio, unlimited);
@ -464,7 +479,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leavi
return find_leaving_and_t_with_breakpoints(entering, t);
X theta;
bool unlimited = get_harris_theta(theta);
lp_assert(unlimited || theta >= zero_of_type<X>());
SASSERT(unlimited || theta >= zero_of_type<X>());
if (try_jump_to_another_bound_on_entering(entering, theta, t, unlimited)) return entering;
if (unlimited)
return -1;
@ -529,11 +544,11 @@ template <typename T, typename X> X lp_primal_core_solver<T, X>::get_max_boun
return ret;
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X> void lp_primal_core_solver<T, X>::check_Ax_equal_b() {
dense_matrix<T, X> d(this->m_A);
T * ls = d.apply_from_left_with_different_dims(this->m_x);
lp_assert(vectors_are_equal<T>(ls, this->m_b, this->m_m()));
SASSERT(vectors_are_equal<T>(ls, this->m_b, this->m_m()));
delete [] ls;
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::check_the_bounds() {
@ -543,8 +558,8 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::check_the
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::check_bound(unsigned i) {
lp_assert (!(this->column_has_low_bound(i) && (numeric_traits<T>::zero() > this->m_x[i])));
lp_assert (!(this->column_has_upper_bound(i) && (this->m_upper_bounds[i] < this->m_x[i])));
SASSERT (!(this->column_has_low_bound(i) && (numeric_traits<T>::zero() > this->m_x[i])));
SASSERT (!(this->column_has_upper_bound(i) && (this->m_upper_bounds[i] < this->m_x[i])));
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::check_correctness() {
@ -558,10 +573,10 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::check_cor
template <typename T, typename X>
void lp_primal_core_solver<T, X>::update_reduced_costs_from_pivot_row(unsigned entering, unsigned leaving) {
// the basis heading has changed already
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
auto & basis_heading = this->m_basis_heading;
lp_assert(basis_heading[entering] >= 0 && static_cast<unsigned>(basis_heading[entering]) < this->m_m());
lp_assert(basis_heading[leaving] < 0);
SASSERT(basis_heading[entering] >= 0 && static_cast<unsigned>(basis_heading[entering]) < this->m_m());
SASSERT(basis_heading[leaving] < 0);
#endif
T pivot = this->m_pivot_row[entering];
T dq = this->m_d[entering]/pivot;
@ -584,7 +599,7 @@ void lp_primal_core_solver<T, X>::update_reduced_costs_from_pivot_row(unsigned e
template <typename T, typename X> int lp_primal_core_solver<T, X>::refresh_reduced_cost_at_entering_and_check_that_it_is_off(unsigned entering) {
if (numeric_traits<T>::precise()) return 0;
T reduced_at_entering_was = this->m_d[entering]; // can benefit from going over non-zeros of m_ed
lp_assert(abs(reduced_at_entering_was) > m_epsilon_of_reduced_cost);
SASSERT(abs(reduced_at_entering_was) > m_epsilon_of_reduced_cost);
T refreshed_cost = this->m_costs[entering];
unsigned i = this->m_m();
while (i--) refreshed_cost -= this->m_costs[this->m_basis[i]] * this->m_ed[i];
@ -619,7 +634,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::backup_an
m_costs_backup = this->m_costs;
} else {
T cost_max = std::max(max_abs_in_vector(this->m_costs), T(1));
lp_assert(m_costs_backup.size() == 0);
SASSERT(m_costs_backup.size() == 0);
for (unsigned j = 0; j < this->m_costs.size(); j++)
m_costs_backup.push_back(this->m_costs[j] /= cost_max);
}
@ -649,16 +664,16 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run(
template <typename T, typename X> void lp_primal_core_solver<T, X>::calc_working_vector_beta_for_column_norms(){
lp_assert(numeric_traits<T>::precise() == false);
lp_assert(this->m_ed.is_OK());
lp_assert(m_beta.is_OK());
SASSERT(numeric_traits<T>::precise() == false);
SASSERT(this->m_ed.is_OK());
SASSERT(m_beta.is_OK());
m_beta = this->m_ed;
this->m_factorization->solve_yB_with_error_check_indexed(m_beta, this->m_basis_heading, this->m_basis, this->m_settings);
}
template <typename T, typename X>
void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving(int entering, X & t) {
lp_assert(!this->A_mult_x_is_off() );
SASSERT(!this->A_mult_x_is_off() );
this->update_x(entering, t * m_sign_of_entering_delta);
if (this->A_mult_x_is_off_on_index(this->m_ed.m_index) && !this->find_x_by_solving()) {
this->init_lu();
@ -670,7 +685,7 @@ void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving(int entering
}
}
if (this->m_using_infeas_costs) {
lp_assert(is_zero(this->m_costs[entering]));
SASSERT(is_zero(this->m_costs[entering]));
init_infeasibility_costs_for_changed_basis_only();
}
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
@ -683,10 +698,10 @@ void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving(int entering
}
template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_entering_and_leaving(int entering, int leaving, X & t) {
lp_assert(entering >= 0 && m_non_basis_list.back() == static_cast<unsigned>(entering));
lp_assert(this->m_using_infeas_costs || t >= zero_of_type<X>());
lp_assert(leaving >= 0 && entering >= 0);
lp_assert(entering != leaving || !is_zero(t)); // otherwise nothing changes
SASSERT(entering >= 0 && m_non_basis_list.back() == static_cast<unsigned>(entering));
SASSERT(this->m_using_infeas_costs || t >= zero_of_type<X>());
SASSERT(leaving >= 0 && entering >= 0);
SASSERT(entering != leaving || !is_zero(t)); // otherwise nothing changes
if (entering == leaving) {
advance_on_entering_equal_leaving(entering, t);
return;
@ -698,14 +713,14 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
int pivot_compare_result = this->pivots_in_column_and_row_are_different(entering, leaving);
if (!pivot_compare_result){;}
else if (pivot_compare_result == 2) { // the sign is changed, cannot continue
this->set_status(lp_status::UNSTABLE);
this->set_status(UNSTABLE);
this->iters_with_no_cost_growing()++;
return;
} else {
lp_assert(pivot_compare_result == 1);
SASSERT(pivot_compare_result == 1);
this->init_lu();
if (this->m_factorization == nullptr || this->m_factorization->get_status() != LU_status::OK) {
this->set_status(lp_status::UNSTABLE);
this->set_status(UNSTABLE);
this->iters_with_no_cost_growing()++;
return;
}
@ -717,10 +732,10 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
t = -t;
}
if (!this->update_basis_and_x(entering, leaving, t)) {
if (this->get_status() == lp_status::FLOATING_POINT_ERROR)
if (this->get_status() == FLOATING_POINT_ERROR)
return;
if (this->m_look_for_feasible_solution_only) {
this->set_status(lp_status::FLOATING_POINT_ERROR);
this->set_status(FLOATING_POINT_ERROR);
return;
}
init_reduced_costs();
@ -733,7 +748,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
}
if (this->current_x_is_feasible()) {
this->set_status(lp_status::FEASIBLE);
this->set_status(FEASIBLE);
if (this->m_look_for_feasible_solution_only)
return;
}
@ -746,7 +761,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
} else {
update_reduced_costs_from_pivot_row(entering, leaving);
}
lp_assert(!need_to_switch_costs());
SASSERT(!need_to_switch_costs());
std::list<unsigned>::iterator it = m_non_basis_list.end();
it--;
* it = static_cast<unsigned>(leaving);
@ -754,13 +769,13 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_entering_precise(int entering) {
lp_assert(numeric_traits<T>::precise());
lp_assert(entering > -1);
SASSERT(numeric_traits<T>::precise());
SASSERT(entering > -1);
this->solve_Bd(entering);
X t;
int leaving = find_leaving_and_t_precise(entering, t);
if (leaving == -1) {
this->set_status(lp_status::UNBOUNDED);
this->set_status(UNBOUNDED);
return;
}
advance_on_entering_and_leaving(entering, leaving, t);
@ -771,12 +786,12 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_e
advance_on_entering_precise(entering);
return;
}
lp_assert(entering > -1);
SASSERT(entering > -1);
this->solve_Bd(entering);
int refresh_result = refresh_reduced_cost_at_entering_and_check_that_it_is_off(entering);
if (refresh_result) {
if (this->m_look_for_feasible_solution_only) {
this->set_status(lp_status::FLOATING_POINT_ERROR);
this->set_status(FLOATING_POINT_ERROR);
return;
}
@ -791,7 +806,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_e
int leaving = find_leaving_and_t(entering, t);
if (leaving == -1){
if (!this->current_x_is_feasible()) {
lp_assert(!numeric_traits<T>::precise()); // we cannot have unbounded with inf costs
SASSERT(!numeric_traits<T>::precise()); // we cannot have unbounded with inf costs
// if (m_look_for_feasible_solution_only) {
// this->m_status = INFEASIBLE;
@ -799,19 +814,19 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_e
// }
if (this->get_status() == lp_status::UNSTABLE) {
this->set_status(lp_status::FLOATING_POINT_ERROR);
if (this->get_status() == UNSTABLE) {
this->set_status(FLOATING_POINT_ERROR);
return;
}
init_infeasibility_costs();
this->set_status(lp_status::UNSTABLE);
this->set_status(UNSTABLE);
return;
}
if (this->get_status() == lp_status::TENTATIVE_UNBOUNDED) {
this->set_status(lp_status::UNBOUNDED);
if (this->get_status() == TENTATIVE_UNBOUNDED) {
this->set_status(UNBOUNDED);
} else {
this->set_status(lp_status::TENTATIVE_UNBOUNDED);
this->set_status(TENTATIVE_UNBOUNDED);
}
return;
}
@ -825,7 +840,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 = static_cast<unsigned>(this->m_nbasis.size());
if (this->get_status() == lp_status::TENTATIVE_UNBOUNDED)
if (this->get_status() == TENTATIVE_UNBOUNDED)
return ret; // we really need to find entering with a large reduced cost
if (ret > 300) {
ret = (unsigned)(ret * this->m_settings.percent_of_entering_to_check / 100);
@ -852,12 +867,12 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
init_run();
if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only) {
this->set_status(lp_status::FEASIBLE);
this->set_status(FEASIBLE);
return 0;
}
if ((!numeric_traits<T>::precise()) && this->A_mult_x_is_off()) {
this->set_status(lp_status::FLOATING_POINT_ERROR);
this->set_status(FLOATING_POINT_ERROR);
return 0;
}
do {
@ -865,10 +880,10 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
return this->total_iterations();
}
one_iteration();
lp_assert(!this->m_using_infeas_costs || this->costs_on_nbasis_are_zeros());
SASSERT(!this->m_using_infeas_costs || this->costs_on_nbasis_are_zeros());
switch (this->get_status()) {
case lp_status::OPTIMAL: // double check that we are at optimum
case lp_status::INFEASIBLE:
case OPTIMAL: // double check that we are at optimum
case INFEASIBLE:
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
break;
if (!numeric_traits<T>::precise()) {
@ -877,7 +892,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status (lp_status::FLOATING_POINT_ERROR);
this->set_status (FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
@ -885,7 +900,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
decide_on_status_when_cannot_find_entering();
break;
}
this->set_status(lp_status::UNKNOWN);
this->set_status(UNKNOWN);
} else { // precise case
if (this->m_look_for_feasible_solution_only) { // todo: keep the reduced costs correct all the time!
init_reduced_costs();
@ -893,31 +908,31 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
decide_on_status_when_cannot_find_entering();
break;
}
this->set_status(lp_status::UNKNOWN);
this->set_status(UNKNOWN);
}
}
break;
case lp_status::TENTATIVE_UNBOUNDED:
case TENTATIVE_UNBOUNDED:
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(lp_status::FLOATING_POINT_ERROR);
this->set_status(FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
break;
case lp_status::UNBOUNDED:
case UNBOUNDED:
if (this->current_x_is_infeasible()) {
init_reduced_costs();
this->set_status(lp_status::UNKNOWN);
this->set_status(UNKNOWN);
}
break;
case lp_status::UNSTABLE:
lp_assert(! (numeric_traits<T>::precise()));
case UNSTABLE:
SASSERT(! (numeric_traits<T>::precise()));
this->init_lu();
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(lp_status::FLOATING_POINT_ERROR);
this->set_status(FLOATING_POINT_ERROR);
break;
}
init_reduced_costs();
@ -926,13 +941,13 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
default:
break; // do nothing
}
} while (this->get_status() != lp_status::FLOATING_POINT_ERROR
} while (this->get_status() != FLOATING_POINT_ERROR
&&
this->get_status() != lp_status::UNBOUNDED
this->get_status() != UNBOUNDED
&&
this->get_status() != lp_status::OPTIMAL
this->get_status() != OPTIMAL
&&
this->get_status() != lp_status::INFEASIBLE
this->get_status() != INFEASIBLE
&&
this->iters_with_no_cost_growing() <= this->m_settings.max_number_of_iterations_with_no_improvements
&&
@ -940,7 +955,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
&&
!(this->current_x_is_feasible() && this->m_look_for_feasible_solution_only));
lp_assert(this->get_status() == lp_status::FLOATING_POINT_ERROR
SASSERT(this->get_status() == FLOATING_POINT_ERROR
||
this->current_x_is_feasible() == false
||
@ -957,7 +972,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::delete_fa
// according to Swietanowski, " A new steepest edge approximation for the simplex method for linear programming"
template <typename T, typename X> void lp_primal_core_solver<T, X>::init_column_norms() {
lp_assert(numeric_traits<T>::precise() == false);
SASSERT(numeric_traits<T>::precise() == false);
for (unsigned j = 0; j < this->m_n(); j++) {
this->m_column_norms[j] = T(static_cast<int>(this->m_A.m_columns[j].size() + 1))
@ -967,7 +982,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_column_
// debug only
template <typename T, typename X> T lp_primal_core_solver<T, X>::calculate_column_norm_exactly(unsigned j) {
lp_assert(numeric_traits<T>::precise() == false);
SASSERT(numeric_traits<T>::precise() == false);
indexed_vector<T> w(this->m_m());
this->m_A.copy_column_to_vector(j, w);
vector<T> d(this->m_m());
@ -979,8 +994,8 @@ template <typename T, typename X> T lp_primal_core_solver<T, X>::calculate_colum
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::update_or_init_column_norms(unsigned entering, unsigned leaving) {
lp_assert(numeric_traits<T>::precise() == false);
lp_assert(m_column_norm_update_counter <= this->m_settings.column_norms_update_frequency);
SASSERT(numeric_traits<T>::precise() == false);
SASSERT(m_column_norm_update_counter <= this->m_settings.column_norms_update_frequency);
if (m_column_norm_update_counter == this->m_settings.column_norms_update_frequency) {
m_column_norm_update_counter = 0;
init_column_norms();
@ -992,7 +1007,7 @@ 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) {
lp_assert(numeric_traits<T>::precise() == false);
SASSERT(numeric_traits<T>::precise() == false);
T pivot = this->m_pivot_row[entering];
T g_ent = calculate_norm_of_entering_exactly() / pivot / pivot;
if (!numeric_traits<T>::precise()) {
@ -1027,8 +1042,8 @@ template <typename T, typename X> T lp_primal_core_solver<T, X>::calculate_no
// calling it stage1 is too cryptic
template <typename T, typename X> void lp_primal_core_solver<T, X>::find_feasible_solution() {
this->m_look_for_feasible_solution_only = true;
lp_assert(this->non_basic_columns_are_set_correctly());
this->set_status(lp_status::UNKNOWN);
SASSERT(this->non_basic_columns_are_set_correctly());
this->set_status(UNKNOWN);
solve();
}
@ -1072,15 +1087,15 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::fill_breakpo
template <typename T, typename X> bool lp_primal_core_solver<T, X>::done() {
if (this->get_status() == lp_status::OPTIMAL || this->get_status() == lp_status::FLOATING_POINT_ERROR) return true;
if (this->get_status() == lp_status::INFEASIBLE) {
if (this->get_status() == OPTIMAL || this->get_status() == FLOATING_POINT_ERROR) return true;
if (this->get_status() == INFEASIBLE) {
return true;
}
if (this->m_iters_with_no_cost_growing >= this->m_settings.max_number_of_iterations_with_no_improvements) {
this->get_status() = lp_status::ITERATIONS_EXHAUSTED; return true;
this->get_status() = ITERATIONS_EXHAUSTED; return true;
}
if (this->total_iterations() >= this->m_settings.max_total_number_of_iterations) {
this->get_status() = lp_status::ITERATIONS_EXHAUSTED; return true;
this->get_status() = ITERATIONS_EXHAUSTED; return true;
}
return false;
}
@ -1095,8 +1110,8 @@ void lp_primal_core_solver<T, X>::init_infeasibility_costs_for_changed_basis_onl
template <typename T, typename X>
void lp_primal_core_solver<T, X>::init_infeasibility_costs() {
lp_assert(this->m_x.size() >= this->m_n());
lp_assert(this->m_column_types.size() >= this->m_n());
SASSERT(this->m_x.size() >= this->m_n());
SASSERT(this->m_column_types.size() >= this->m_n());
for (unsigned j = this->m_n(); j--;)
init_infeasibility_cost_for_column(j);
this->m_using_infeas_costs = true;
@ -1138,7 +1153,7 @@ lp_primal_core_solver<T, X>::get_infeasibility_cost_for_column(unsigned j) const
ret = numeric_traits<T>::zero();
break;
default:
lp_assert(false);
SASSERT(false);
ret = numeric_traits<T>::zero(); // does not matter
break;
}
@ -1192,14 +1207,14 @@ lp_primal_core_solver<T, X>::init_infeasibility_cost_for_column(unsigned j) {
this->m_costs[j] = numeric_traits<T>::zero();
break;
default:
lp_assert(false);
SASSERT(false);
break;
}
if (numeric_traits<T>::is_zero(this->m_costs[j])) {
this->remove_column_from_inf_set(j);
this->m_inf_set.erase(j);
} else {
this->insert_column_into_inf_set(j);
this->m_inf_set.insert(j);
}
if (!this->m_settings.use_breakpoints_in_feasibility_search) {
this->m_costs[j] = - this->m_costs[j];
@ -1223,7 +1238,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::print_column
case column_type::free_column:
out << "( _" << this->m_x[j] << "_)" << std::endl;
default:
lp_unreachable();
SASSERT(false);
}
}
@ -1262,7 +1277,7 @@ template <typename T, typename X> std::string lp_primal_core_solver<T, X>::break
case upper_break: return "upper_break";
case fixed_break: return "fixed_break";
default:
lp_assert(false);
SASSERT(false);
break;
}
return "type is not found";
@ -1275,7 +1290,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::print_breakp
template <typename T, typename X>
void lp_primal_core_solver<T, X>::init_reduced_costs() {
lp_assert(!this->use_tableau());
SASSERT(!this->use_tableau());
if (this->current_x_is_infeasible() && !this->m_using_infeas_costs) {
init_infeasibility_costs();
} else if (this->current_x_is_feasible() && this->m_using_infeas_costs) {
@ -1290,12 +1305,12 @@ void lp_primal_core_solver<T, X>::init_reduced_costs() {
template <typename T, typename X> void lp_primal_core_solver<T, X>::change_slope_on_breakpoint(unsigned entering, breakpoint<X> * b, T & slope_at_entering) {
if (b->m_j == entering) {
lp_assert(b->m_type != fixed_break && (!is_zero(b->m_delta)));
SASSERT(b->m_type != fixed_break && (!is_zero(b->m_delta)));
slope_at_entering += m_sign_of_entering_delta;
return;
}
lp_assert(this->m_basis_heading[b->m_j] >= 0);
SASSERT(this->m_basis_heading[b->m_j] >= 0);
unsigned i_row = this->m_basis_heading[b->m_j];
const T & d = - this->m_ed[i_row];
if (numeric_traits<T>::is_zero(d)) return;
@ -1314,13 +1329,13 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::change_sl
slope_at_entering += delta;
break;
default:
lp_assert(false);
SASSERT(false);
}
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::try_add_breakpoint_in_row(unsigned i) {
lp_assert(i < this->m_m());
SASSERT(i < this->m_m());
const T & d = this->m_ed[i]; // the coefficient before m_entering in the i-th row
if (d == 0) return; // the change of x[m_entering] will not change the corresponding basis x
unsigned j = this->m_basis[i];
@ -1342,7 +1357,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::try_add_b
case column_type::free_column:
break;
default:
lp_assert(false);
SASSERT(false);
break;
}
}
@ -1366,7 +1381,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::print_bound_
out << "inf, inf" << std::endl;
break;
default:
lp_assert(false);
SASSERT(false);
break;
}
}

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <utility>
#include <memory>
#include <string>
@ -9,7 +24,7 @@
#include <functional>
#include "util/lp/lar_solver.h"
#include "util/lp/lp_primal_core_solver.hpp"
#include "util/lp/lp_primal_core_solver_tableau.hpp"
#include "util/lp/lp_primal_core_solver_tableau.h"
namespace lp {
template void lp_primal_core_solver<double, double>::find_feasible_solution();

21
src/util/lp/lp_primal_core_solver_tableau.hpp
View file

@ -97,8 +97,9 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over((this->m_using_infeas_costs? "inf t" : "feas t"), * this->m_settings.get_message_ostream())) {
return this->total_iterations();
}
if (this->m_settings.use_tableau_rows())
if (this->m_settings.use_tableau_rows()) {
one_iteration_tableau_rows();
}
else
one_iteration_tableau();
switch (this->get_status()) {
@ -184,7 +185,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
}
template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_entering_and_leaving_tableau(int entering, int leaving, X & t) {
lp_assert(this->A_mult_x_is_off() == false);
CASSERT("A_off", this->A_mult_x_is_off() == false);
lp_assert(leaving >= 0 && entering >= 0);
lp_assert((this->m_settings.simplex_strategy() ==
simplex_strategy_enum::tableau_rows) ||
@ -201,7 +202,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
t = -t;
}
this->update_basis_and_x_tableau(entering, leaving, t);
lp_assert(this->A_mult_x_is_off() == false);
CASSERT("A_off", this->A_mult_x_is_off() == false);
this->iters_with_no_cost_growing() = 0;
} else {
this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
@ -225,7 +226,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
template <typename T, typename X>
void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving_tableau(int entering, X & t) {
lp_assert(!this->A_mult_x_is_off() );
CASSERT("A_off", !this->A_mult_x_is_off() );
this->update_x_tableau(entering, t * m_sign_of_entering_delta);
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
return;
@ -298,7 +299,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tableau() {
// print_matrix(&(this->m_A), std::cout);
lp_assert(this->A_mult_x_is_off() == false);
CASSERT("A_off", this->A_mult_x_is_off() == false);
lp_assert(basis_columns_are_set_correctly());
this->m_basis_sort_counter = 0; // to initiate the sort of the basis
this->set_total_iterations(0);
@ -331,22 +332,20 @@ update_basis_and_x_tableau(int entering, int leaving, X const & tt) {
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::
update_x_tableau(unsigned entering, const X& delta) {
this->add_delta_to_x_and_call_tracker(entering, delta);
if (!this->m_using_infeas_costs) {
this->m_x[entering] += delta;
for (const auto & c : this->m_A.m_columns[entering]) {
unsigned i = c.m_i;
this->m_x[this->m_basis[i]] -= delta * this->m_A.get_val(c);
this->update_column_in_inf_set(this->m_basis[i]);
this->update_x_with_delta_and_track_feasibility(this->m_basis[i], - delta * this->m_A.get_val(c));
}
} else { // m_using_infeas_costs == true
this->m_x[entering] += delta;
lp_assert(this->column_is_feasible(entering));
lp_assert(this->m_costs[entering] == zero_of_type<T>());
// m_d[entering] can change because of the cost change for basic columns.
for (const auto & c : this->m_A.m_columns[entering]) {
unsigned i = c.m_i;
unsigned j = this->m_basis[i];
this->m_x[j] -= delta * this->m_A.get_val(c);
this->add_delta_to_x_and_call_tracker(j, -delta * this->m_A.get_val(c));
update_inf_cost_for_column_tableau(j);
if (is_zero(this->m_costs[j]))
this->remove_column_from_inf_set(j);
@ -354,7 +353,7 @@ update_x_tableau(unsigned entering, const X& delta) {
this->insert_column_into_inf_set(j);
}
}
lp_assert(this->A_mult_x_is_off() == false);
CASSERT("A_off", this->A_mult_x_is_off() == false);
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
#include <unordered_map>

45
src/util/lp/lp_primal_simplex.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <string>
#include "util/vector.h"
#include "util/lp/lp_primal_simplex.h"
@ -61,7 +76,7 @@ template <typename T, typename X> void lp_primal_simplex<T, X>::fill_costs_and_x
int row,
unsigned & slack_var,
unsigned & artificial) {
lp_assert(row >= 0 && row < this->row_count());
SASSERT(row >= 0 && row < this->row_count());
auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[row]];
// we need to bring the program to the form Ax = b
T rs = this->m_b[row];
@ -86,7 +101,7 @@ template <typename T, typename X> void lp_primal_simplex<T, X>::fill_costs_and_x
(*this->m_A)(row, slack_var) = - numeric_traits<T>::one();
if (rs > 0) {
lp_assert(numeric_traits<T>::is_zero(this->m_x[slack_var]));
SASSERT(numeric_traits<T>::is_zero(this->m_x[slack_var]));
// adding one artificial
this->m_column_types[artificial] = column_type::low_bound;
(*this->m_A)(row, artificial) = numeric_traits<T>::one();
@ -108,7 +123,7 @@ template <typename T, typename X> void lp_primal_simplex<T, X>::fill_costs_and_x
if (rs < 0) {
// adding one artificial
lp_assert(numeric_traits<T>::is_zero(this->m_x[slack_var]));
SASSERT(numeric_traits<T>::is_zero(this->m_x[slack_var]));
this->m_column_types[artificial] = column_type::low_bound;
(*this->m_A)(row, artificial) = - numeric_traits<T>::one();
this->m_costs[artificial] = artificial_cost;
@ -177,12 +192,12 @@ 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>::fill_A_x_and_basis_for_stage_one_total_inf_for_row(unsigned row) {
lp_assert(row < this->row_count());
SASSERT(row < this->row_count());
auto ext_row_it = this->m_core_solver_rows_to_external_rows.find(row);
lp_assert(ext_row_it != this->m_core_solver_rows_to_external_rows.end());
SASSERT(ext_row_it != this->m_core_solver_rows_to_external_rows.end());
unsigned ext_row = ext_row_it->second;
auto constr_it = this->m_constraints.find(ext_row);
lp_assert(constr_it != this->m_constraints.end());
SASSERT(constr_it != this->m_constraints.end());
auto & constraint = constr_it->second;
unsigned j = this->m_A->column_count(); // j is a slack variable
this->m_A->add_column();
@ -209,14 +224,14 @@ template <typename T, typename X> void lp_primal_simplex<T, X>::fill_A_x_and_bas
this->m_upper_bounds[j] = m_low_bounds[j] = zero_of_type<X>();
break;
default:
lp_unreachable();
SASSERT(false);
}
}
template <typename T, typename X> void lp_primal_simplex<T, X>::solve_with_total_inf() {
int total_vars = this->m_A->column_count() + this->row_count();
if (total_vars == 0) {
this->m_status = lp_status::OPTIMAL;
this->m_status = OPTIMAL;
return;
}
m_low_bounds.clear();
@ -281,10 +296,10 @@ template <typename T, typename X> T lp_primal_simplex<T, X>::get_row_value(unsig
T ret = numeric_traits<T>::zero();
for (auto & pair : it->second) {
auto cit = this->m_map_from_var_index_to_column_info.find(pair.first);
lp_assert(cit != this->m_map_from_var_index_to_column_info.end());
SASSERT(cit != this->m_map_from_var_index_to_column_info.end());
column_info<T> * ci = cit->second;
auto sol_it = solution.find(ci->get_name());
lp_assert(sol_it != solution.end());
SASSERT(sol_it != solution.end());
T column_val = sol_it->second;
if (out != nullptr) {
(*out) << pair.second << "(" << ci->get_name() << "=" << column_val << ") ";
@ -329,7 +344,7 @@ template <typename T, typename X> bool lp_primal_simplex<T, X>::row_constraint_h
}
return true;;
}
lp_unreachable();
SASSERT(false);
return false; // it is unreachable
}

23
src/util/lp/lp_primal_simplex_instances.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <utility>
#include <memory>
#include <string>

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -16,17 +31,13 @@ namespace lp {
typedef unsigned var_index;
typedef unsigned constraint_index;
typedef unsigned row_index;
typedef vector<std::pair<mpq, constraint_index>> explanation_t;
enum class column_type {
free_column = 0,
low_bound = 1,
upper_bound = 2,
boxed = 3,
fixed = 4
};
low_bound = 1,
upper_bound = 2,
boxed = 3,
fixed = 4
};
enum class simplex_strategy_enum {
undecided = 3,
@ -37,7 +48,7 @@ enum class simplex_strategy_enum {
std::string column_type_to_string(column_type t);
enum class lp_status {
enum lp_status {
UNKNOWN,
INFEASIBLE,
TENTATIVE_UNBOUNDED,
@ -79,14 +90,11 @@ public:
};
struct stats {
unsigned m_make_feasible;
unsigned m_total_iterations;
unsigned m_iters_with_no_cost_growing;
unsigned m_num_factorizations;
unsigned m_num_of_implied_bounds;
unsigned m_need_to_solve_inf;
unsigned m_max_cols;
unsigned m_max_rows;
stats() { reset(); }
void reset() { memset(this, 0, sizeof(*this)); }
};
@ -205,8 +213,7 @@ public:
use_breakpoints_in_feasibility_search(false),
max_row_length_for_bound_propagation(300),
backup_costs(true),
column_number_threshold_for_using_lu_in_lar_solver(4000),
m_int_branch_cut_threshold(100)
column_number_threshold_for_using_lu_in_lar_solver(4000)
{}
void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; }
@ -304,7 +311,7 @@ public:
unsigned column_norms_update_frequency;
bool scale_with_ratio;
double density_threshold; // need to tune it up, todo
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
static unsigned ddd; // used for debugging
#endif
bool use_breakpoints_in_feasibility_search;
@ -313,7 +320,6 @@ public:
unsigned max_row_length_for_bound_propagation;
bool backup_costs;
unsigned column_number_threshold_for_using_lu_in_lar_solver;
unsigned m_int_branch_cut_threshold;
}; // end of lp_settings class
@ -375,7 +381,7 @@ inline void print_blanks(int n, std::ostream & out) {
// after a push of the last element we ensure that the vector increases
// we also suppose that before the last push the vector was increasing
inline void ensure_increasing(vector<unsigned> & v) {
lp_assert(v.size() > 0);
SASSERT(v.size() > 0);
unsigned j = v.size() - 1;
for (; j > 0; j-- )
if (v[j] <= v[j - 1]) {
@ -390,7 +396,7 @@ inline void ensure_increasing(vector<unsigned> & v) {
#if LEAN_DEBUG
#if Z3DEBUG
bool D();
#endif
}

55
src/util/lp/lp_settings.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <cmath>
#include <string>
#include "util/vector.h"
@ -14,27 +29,27 @@ std::string column_type_to_string(column_type t) {
case column_type::low_bound: return "low_bound";
case column_type::upper_bound: return "upper_bound";
case column_type::free_column: return "free_column";
default: lp_unreachable();
default: SASSERT(false);
}
return "unknown"; // it is unreachable
}
const char* lp_status_to_string(lp_status status) {
switch (status) {
case lp_status::UNKNOWN: return "UNKNOWN";
case lp_status::INFEASIBLE: return "INFEASIBLE";
case lp_status::UNBOUNDED: return "UNBOUNDED";
case lp_status::TENTATIVE_DUAL_UNBOUNDED: return "TENTATIVE_DUAL_UNBOUNDED";
case lp_status::DUAL_UNBOUNDED: return "DUAL_UNBOUNDED";
case lp_status::OPTIMAL: return "OPTIMAL";
case lp_status::FEASIBLE: return "FEASIBLE";
case lp_status::FLOATING_POINT_ERROR: return "FLOATING_POINT_ERROR";
case lp_status::TIME_EXHAUSTED: return "TIME_EXHAUSTED";
case lp_status::ITERATIONS_EXHAUSTED: return "ITERATIONS_EXHAUSTED";
case lp_status::EMPTY: return "EMPTY";
case lp_status::UNSTABLE: return "UNSTABLE";
case UNKNOWN: return "UNKNOWN";
case INFEASIBLE: return "INFEASIBLE";
case UNBOUNDED: return "UNBOUNDED";
case TENTATIVE_DUAL_UNBOUNDED: return "TENTATIVE_DUAL_UNBOUNDED";
case DUAL_UNBOUNDED: return "DUAL_UNBOUNDED";
case OPTIMAL: return "OPTIMAL";
case FEASIBLE: return "FEASIBLE";
case FLOATING_POINT_ERROR: return "FLOATING_POINT_ERROR";
case TIME_EXHAUSTED: return "TIME_EXHAUSTED";
case ITERATIONS_EXHAUSTED: return "ITERATIONS_EXHAUSTED";
case EMPTY: return "EMPTY";
case UNSTABLE: return "UNSTABLE";
default:
lp_unreachable();
SASSERT(false);
}
return "UNKNOWN"; // it is unreachable
}
@ -49,7 +64,7 @@ lp_status lp_status_from_string(std::string status) {
if (status == "TIME_EXHAUSTED") return lp_status::TIME_EXHAUSTED;
if (status == "ITERATIONS_EXHAUSTED") return lp_status::ITERATIONS_EXHAUSTED;
if (status == "EMPTY") return lp_status::EMPTY;
lp_unreachable();
SASSERT(false);
return lp_status::UNKNOWN; // it is unreachable
}
@ -104,7 +119,7 @@ bool vectors_are_equal(const vector<T> & a, const vector<T> &b) {
}
return true;
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
unsigned lp_settings::ddd = 0;
#endif
}

23
src/util/lp/lp_settings_instances.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <memory>
#include "util/vector.h"
#include "util/lp/lp_settings.hpp"

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <string>

63
src/util/lp/lp_solver.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <string>
#include <algorithm>
#include "util/vector.h"
@ -32,7 +47,7 @@ template <typename T, typename X> T lp_solver<T, X>::get_column_cost_value(unsig
return ci->get_cost() * get_column_value(j);
}
template <typename T, typename X> void lp_solver<T, X>::add_constraint(lp_relation relation, T right_side, unsigned row_index) {
lp_assert(m_constraints.find(row_index) == m_constraints.end());
SASSERT(m_constraints.find(row_index) == m_constraints.end());
lp_constraint<T, X> cs(right_side, relation);
m_constraints[row_index] = cs;
}
@ -158,10 +173,10 @@ template <typename T, typename X> void lp_solver<T, X>::pin_vars_on_row_with_sig
column_info<T> * ci = m_map_from_var_index_to_column_info[j];
T a = t.second;
if (a * sign > numeric_traits<T>::zero()) {
lp_assert(ci->upper_bound_is_set());
SASSERT(ci->upper_bound_is_set());
ci->set_fixed_value(ci->get_upper_bound());
} else {
lp_assert(ci->low_bound_is_set());
SASSERT(ci->low_bound_is_set());
ci->set_fixed_value(ci->get_low_bound());
}
}
@ -223,7 +238,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_e_is_obsolete(std
T rs = m_constraints[row_index].m_rs;
if (row_is_zero(row)) {
if (!is_zero(rs))
m_status = lp_status::INFEASIBLE;
m_status = INFEASIBLE;
return true;
}
@ -233,7 +248,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_e_is_obsolete(std
T diff = low_bound - rs;
if (!val_is_smaller_than_eps(diff, m_settings.refactor_tolerance)){
// low_bound > rs + m_settings.refactor_epsilon
m_status = lp_status::INFEASIBLE;
m_status = INFEASIBLE;
return true;
}
if (val_is_smaller_than_eps(-diff, m_settings.refactor_tolerance)){
@ -248,7 +263,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_e_is_obsolete(std
T diff = rs - upper_bound;
if (!val_is_smaller_than_eps(diff, m_settings.refactor_tolerance)) {
// upper_bound < rs - m_settings.refactor_tolerance
m_status = lp_status::INFEASIBLE;
m_status = INFEASIBLE;
return true;
}
if (val_is_smaller_than_eps(-diff, m_settings.refactor_tolerance)){
@ -264,7 +279,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_ge_is_obsolete(std:
T rs = m_constraints[row_index].m_rs;
if (row_is_zero(row)) {
if (rs > zero_of_type<X>())
m_status = lp_status::INFEASIBLE;
m_status = INFEASIBLE;
return true;
}
@ -273,7 +288,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_ge_is_obsolete(std:
T diff = rs - upper_bound;
if (!val_is_smaller_than_eps(diff, m_settings.refactor_tolerance)) {
// upper_bound < rs - m_settings.refactor_tolerance
m_status = lp_status::INFEASIBLE;
m_status = INFEASIBLE;
return true;
}
if (val_is_smaller_than_eps(-diff, m_settings.refactor_tolerance)){
@ -290,7 +305,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_le_is_obsolete(std::
T rs = m_constraints[row_index].m_rs;
if (row_is_zero(row)) {
if (rs < zero_of_type<X>())
m_status = lp_status::INFEASIBLE;
m_status = INFEASIBLE;
return true;
}
@ -328,7 +343,7 @@ template <typename T, typename X> bool lp_solver<T, X>::row_is_obsolete(std::
case lp_relation::Less_or_equal:
return row_le_is_obsolete(row, row_index);
}
lp_unreachable();
SASSERT(false);
return false; // it is unreachable
}
@ -343,7 +358,7 @@ template <typename T, typename X> void lp_solver<T, X>::remove_fixed_or_zero_col
vector<unsigned> removed;
for (auto & col : row) {
unsigned j = col.first;
lp_assert(m_map_from_var_index_to_column_info.find(j) != m_map_from_var_index_to_column_info.end());
SASSERT(m_map_from_var_index_to_column_info.find(j) != m_map_from_var_index_to_column_info.end());
column_info<T> * ci = m_map_from_var_index_to_column_info[j];
if (ci->is_fixed()) {
removed.push_back(j);
@ -412,7 +427,7 @@ template <typename T, typename X> void lp_solver<T, X>::map_external_columns_to_
}
unsigned j = col.first;
auto column_info_it = m_map_from_var_index_to_column_info.find(j);
lp_assert(column_info_it != m_map_from_var_index_to_column_info.end());
SASSERT(column_info_it != m_map_from_var_index_to_column_info.end());
auto j_column = column_info_it->second->get_column_index();
if (!is_valid(j_column)) { // j is a newcomer
@ -435,14 +450,14 @@ template <typename T, typename X> void lp_solver<T, X>::fill_A_from_A_values() {
m_A = new static_matrix<T, X>(static_cast<unsigned>(m_A_values.size()), number_of_core_structurals());
for (auto & t : m_A_values) {
auto row_it = m_external_rows_to_core_solver_rows.find(t.first);
lp_assert(row_it != m_external_rows_to_core_solver_rows.end());
SASSERT(row_it != m_external_rows_to_core_solver_rows.end());
unsigned row = row_it->second;
for (auto k : t.second) {
auto column_info_it = m_map_from_var_index_to_column_info.find(k.first);
lp_assert(column_info_it != m_map_from_var_index_to_column_info.end());
SASSERT(column_info_it != m_map_from_var_index_to_column_info.end());
column_info<T> *ci = column_info_it->second;
unsigned col = ci->get_column_index();
lp_assert(is_valid(col));
SASSERT(is_valid(col));
bool col_is_flipped = m_map_from_var_index_to_column_info[k.first]->is_flipped();
if (!col_is_flipped) {
(*m_A)(row, col) = k.second;
@ -456,7 +471,7 @@ template <typename T, typename X> void lp_solver<T, X>::fill_A_from_A_values() {
template <typename T, typename X> void lp_solver<T, X>::fill_matrix_A_and_init_right_side() {
map_external_rows_to_core_solver_rows();
map_external_columns_to_core_solver_columns();
lp_assert(m_A == nullptr);
SASSERT(m_A == nullptr);
fill_A_from_A_values();
m_b.resize(m_A->row_count());
}
@ -468,7 +483,7 @@ 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) {
lp_assert(this->m_constraints.find(this->m_core_solver_rows_to_external_rows[i]) != this->m_constraints.end());
SASSERT(this->m_constraints.find(this->m_core_solver_rows_to_external_rows[i]) != this->m_constraints.end());
auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[i]];
switch (constraint.m_relation) {
case Equal:
@ -504,7 +519,7 @@ template <typename T, typename X> T lp_solver<T, X>::low_bound_shift_for_row(
template <typename T, typename X> void lp_solver<T, X>::fill_m_b() {
for (int i = this->row_count() - 1; i >= 0; i--) {
lp_assert(this->m_constraints.find(this->m_core_solver_rows_to_external_rows[i]) != this->m_constraints.end());
SASSERT(this->m_constraints.find(this->m_core_solver_rows_to_external_rows[i]) != this->m_constraints.end());
unsigned external_i = this->m_core_solver_rows_to_external_rows[i];
auto & constraint = this->m_constraints[external_i];
this->m_b[i] = constraint.m_rs - low_bound_shift_for_row(external_i);
@ -542,13 +557,13 @@ template <typename T, typename X> T lp_solver<T, X>::get_column_value_with_core_
template <typename T, typename X> void lp_solver<T, X>::set_scaled_cost(unsigned j) {
// grab original costs but modify it with the column scales
lp_assert(j < this->m_column_scale.size());
SASSERT(j < this->m_column_scale.size());
column_info<T> * ci = this->m_map_from_var_index_to_column_info[this->m_core_solver_columns_to_external_columns[j]];
T cost = ci->get_cost();
if (ci->is_flipped()){
cost *= -1;
}
lp_assert(ci->is_fixed() == false);
SASSERT(ci->is_fixed() == false);
this->m_costs[j] = cost * this->m_column_scale[j];
}
}

23
src/util/lp/lp_solver_instances.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <string>
#include "util/lp/lp_solver.hpp"
template void lp::lp_solver<double, double>::add_constraint(lp::lp_relation, double, unsigned int);

27
src/util/lp/lp_utils.cpp
View file

@ -1,11 +1,26 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/lp_utils.h"
#ifdef lp_for_z3
namespace lp {
double numeric_traits<double>::g_zero = 0.0;
double numeric_traits<double>::g_one = 1.0;
}
#endif

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

@ -1,8 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
This file should be present in z3 and in Lean.
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <string>
#include "util/lp/numeric_pair.h"
@ -18,23 +32,15 @@ bool try_get_val(const std::unordered_map<A,B> & map, const A& key, B & val) {
template <typename A, typename B>
bool contains(const std::unordered_map<A, B> & map, const A& key) {
return map.find(key) != map.end();
return map.find(key) != map.end();
}
#ifdef lp_for_z3
#ifdef Z3DEBUG
#define LEAN_DEBUG 1
#endif
namespace lp {
inline void throw_exception(const std::string & str) {
throw default_exception(str);
}
typedef z3_exception exception;
#define lp_assert(_x_) { SASSERT(_x_); }
inline void lp_unreachable() { lp_assert(false); }
template <typename X> inline X zero_of_type() { return numeric_traits<X>::zero(); }
template <typename X> inline X one_of_type() { return numeric_traits<X>::one(); }
template <typename X> inline bool is_zero(const X & v) { return numeric_traits<X>::is_zero(v); }
@ -78,64 +84,3 @@ struct hash<lp::numeric_pair<lp::mpq>> {
};
}
#else // else of #if lp_for_z3
#include <utility>
#include <functional>
//include "util/numerics/mpq.h"
//include "util/numerics/numeric_traits.h"
//include "util/numerics/double.h"
#ifdef __CLANG__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wmismatched-tags"
#endif
namespace std {
template<>
struct hash<lp::mpq> {
inline size_t operator()(const lp::mpq & v) const {
return v.hash();
}
};
}
namespace lp {
template <typename X> inline bool precise() { return numeric_traits<X>::precise();}
template <typename X> inline X one_of_type() { return numeric_traits<X>::one(); }
template <typename X> inline bool is_zero(const X & v) { return numeric_traits<X>::is_zero(v); }
template <typename X> inline double get_double(const X & v) { return numeric_traits<X>::get_double(v); }
template <typename T> inline T zero_of_type() {return numeric_traits<T>::zero();}
inline void throw_exception(std::string str) { throw exception(str); }
template <typename T> inline T from_string(std::string const & ) { lp_unreachable();}
template <> double inline from_string<double>(std::string const & str) { return atof(str.c_str());}
template <> mpq inline from_string<mpq>(std::string const & str) {
return mpq(atof(str.c_str()));
}
} // closing lp
template <class T>
inline void hash_combine(std::size_t & seed, const T & v) {
seed ^= std::hash<T>()(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
}
namespace std {
template<typename S, typename T> struct hash<pair<S, T>> {
inline size_t operator()(const pair<S, T> & v) const {
size_t seed = 0;
hash_combine(seed, v.first);
hash_combine(seed, v.second);
return seed;
}
};
template<>
struct hash<lp::numeric_pair<lp::mpq>> {
inline size_t operator()(const lp::numeric_pair<lp::mpq> & v) const {
size_t seed = 0;
hash_combine(seed, v.x);
hash_combine(seed, v.y);
return seed;
}
};
} // std
#ifdef __CLANG__
#pragma clang diagnostic pop
#endif
#endif

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
@ -19,7 +34,7 @@
#include "util/lp/square_dense_submatrix.h"
#include "util/lp/dense_matrix.h"
namespace lp {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X> // print the nr x nc submatrix at the top left corner
void print_submatrix(sparse_matrix<T, X> & m, unsigned mr, unsigned nc);
@ -32,7 +47,7 @@ void print_matrix(sparse_matrix<T, X>& m, std::ostream & out);
template <typename T, typename X>
X dot_product(const vector<T> & a, const vector<X> & b) {
lp_assert(a.size() == b.size());
SASSERT(a.size() == b.size());
auto r = zero_of_type<X>();
for (unsigned i = 0; i < a.size(); i++) {
r += a[i] * b[i];
@ -47,7 +62,7 @@ class one_elem_on_diag: public tail_matrix<T, X> {
T m_val;
public:
one_elem_on_diag(unsigned i, T val) : m_i(i), m_val(val) {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
m_one_over_val = numeric_traits<T>::one() / m_val;
#endif
}
@ -56,7 +71,7 @@ public:
one_elem_on_diag(const one_elem_on_diag & o);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
unsigned m_m;
unsigned m_n;
virtual void set_number_of_rows(unsigned m) { m_m = m; m_n = m; }
@ -91,15 +106,15 @@ public:
void conjugate_by_permutation(permutation_matrix<T, X> & p) {
// this = p * this * p(-1)
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// auto rev = p.get_reverse();
// auto deb = ((*this) * rev);
// deb = p * deb;
#endif
m_i = p.apply_reverse(m_i);
#ifdef LEAN_DEBUG
// lp_assert(*this == deb);
#ifdef Z3DEBUG
// SASSERT(*this == deb);
#endif
}
}; // end of one_elem_on_diag
@ -212,7 +227,7 @@ public:
// see page 407 of Chvatal
unsigned transform_U_to_V_by_replacing_column(indexed_vector<T> & w, unsigned leaving_column_of_U);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
void check_vector_w(unsigned entering);
void check_apply_matrix_to_vector(matrix<T, X> *lp, T *w);
@ -248,7 +263,7 @@ public:
bool is_correct(const vector<unsigned>& basis);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
dense_matrix<T, X> tail_product();
dense_matrix<T, X> get_left_side(const vector<unsigned>& basis);
@ -291,7 +306,7 @@ public:
bool need_to_refactor() { return m_refactor_counter >= 200; }
void adjust_dimension_with_matrix_A() {
lp_assert(m_A.row_count() >= m_dim);
SASSERT(m_A.row_count() >= m_dim);
m_dim = m_A.row_count();
m_U.resize(m_dim);
m_Q.resize(m_dim);
@ -305,7 +320,7 @@ public:
unsigned m = m_A.row_count();
unsigned m_prev = m_U.dimension();
lp_assert(m_A.column_count() == heading.size());
SASSERT(m_A.column_count() == heading.size());
for (unsigned i = m_prev; i < m; i++) {
for (const row_cell<T> & c : m_A.m_rows[i]) {
@ -321,14 +336,14 @@ public:
void add_last_rows_to_B(const vector<int> & heading, const std::unordered_set<unsigned> & columns_to_replace) {
unsigned m = m_A.row_count();
lp_assert(m_A.column_count() == heading.size());
SASSERT(m_A.column_count() == heading.size());
adjust_dimension_with_matrix_A();
m_w_for_extension.resize(m);
// At this moment the LU is correct
// for B extended by only by ones at the diagonal in the lower right corner
for (unsigned j :columns_to_replace) {
lp_assert(heading[j] >= 0);
SASSERT(heading[j] >= 0);
replace_column_with_only_change_at_last_rows(j, heading[j]);
if (get_status() == LU_status::Degenerated)
break;
@ -352,7 +367,7 @@ public:
template <typename T, typename X>
void init_factorization(lu<T, X>* & factorization, static_matrix<T, X> & m_A, vector<unsigned> & m_basis, lp_settings &m_settings);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X>
dense_matrix<T, X> get_B(lu<T, X>& f, const vector<unsigned>& basis);
#endif

135
src/util/lp/lu.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <string>
#include <algorithm>
#include <set>
@ -10,7 +25,7 @@
#include "util/debug.h"
#include "util/lp/lu.h"
namespace lp {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X> // print the nr x nc submatrix at the top left corner
void print_submatrix(sparse_matrix<T, X> & m, unsigned mr, unsigned nc, std::ostream & out) {
vector<vector<std::string>> A;
@ -72,13 +87,13 @@ template <typename T, typename X>
one_elem_on_diag<T, X>::one_elem_on_diag(const one_elem_on_diag & o) {
m_i = o.m_i;
m_val = o.m_val;
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
m_m = m_n = o.m_m;
m_one_over_val = numeric_traits<T>::one() / o.m_val;
#endif
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X>
T one_elem_on_diag<T, X>::get_elem(unsigned i, unsigned j) const {
if (i == j){
@ -122,29 +137,29 @@ lu<T, X>::lu(static_matrix<T, X> const & A,
m_failure(false),
m_row_eta_work_vector(A.row_count()),
m_refactor_counter(0) {
lp_assert(!(numeric_traits<T>::precise() && settings.use_tableau()));
#ifdef LEAN_DEBUG
SASSERT(!(numeric_traits<T>::precise() && settings.use_tableau()));
#ifdef Z3DEBUG
debug_test_of_basis(A, basis);
#endif
++m_settings.st().m_num_factorizations;
create_initial_factorization();
#ifdef LEAN_DEBUG
// lp_assert(check_correctness());
#ifdef Z3DEBUG
// SASSERT(check_correctness());
#endif
}
template <typename T, typename X>
void lu<T, X>::debug_test_of_basis(static_matrix<T, X> const & A, vector<unsigned> & basis) {
std::set<unsigned> set;
for (unsigned i = 0; i < A.row_count(); i++) {
lp_assert(basis[i]< A.column_count());
SASSERT(basis[i]< A.column_count());
set.insert(basis[i]);
}
lp_assert(set.size() == A.row_count());
SASSERT(set.size() == A.row_count());
}
template <typename T, typename X>
void lu<T, X>::solve_By(indexed_vector<X> & y) {
lp_assert(false); // not implemented
SASSERT(false); // not implemented
// init_vector_y(y);
// solve_By_when_y_is_ready(y);
}
@ -268,7 +283,7 @@ void lu<T, X>::solve_yB(vector<T>& y) {
m_U.solve_y_U(y); // got y*U=cb*R(-1)
m_Q.apply_reverse_from_right_to_T(y); //
for (auto e = m_tail.rbegin(); e != m_tail.rend(); ++e) {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
(*e)->set_number_of_columns(m_dim);
#endif
(*e)->apply_from_right(y);
@ -277,20 +292,20 @@ void lu<T, X>::solve_yB(vector<T>& y) {
template <typename T, typename X>
void lu<T, X>::solve_yB_indexed(indexed_vector<T>& y) {
lp_assert(y.is_OK());
SASSERT(y.is_OK());
// first solve yU = cb*R(-1)
m_R.apply_reverse_from_right_to_T(y); // got y = cb*R(-1)
lp_assert(y.is_OK());
SASSERT(y.is_OK());
m_U.solve_y_U_indexed(y, m_settings); // got y*U=cb*R(-1)
lp_assert(y.is_OK());
SASSERT(y.is_OK());
m_Q.apply_reverse_from_right_to_T(y);
lp_assert(y.is_OK());
SASSERT(y.is_OK());
for (auto e = m_tail.rbegin(); e != m_tail.rend(); ++e) {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
(*e)->set_number_of_columns(m_dim);
#endif
(*e)->apply_from_right(y);
lp_assert(y.is_OK());
SASSERT(y.is_OK());
}
}
@ -304,8 +319,8 @@ void lu<T, X>::add_delta_to_solution(const vector<T>& yc, vector<T>& y){
template <typename T, typename X>
void lu<T, X>::add_delta_to_solution_indexed(indexed_vector<T>& y) {
// the delta sits in m_y_copy, put result into y
lp_assert(y.is_OK());
lp_assert(m_y_copy.is_OK());
SASSERT(y.is_OK());
SASSERT(m_y_copy.is_OK());
m_ii.clear();
m_ii.resize(y.data_size());
for (unsigned i : y.m_index)
@ -315,7 +330,7 @@ void lu<T, X>::add_delta_to_solution_indexed(indexed_vector<T>& y) {
if (m_ii[i] == 0)
m_ii.set_value(1, i);
}
lp_assert(m_ii.is_OK());
SASSERT(m_ii.is_OK());
y.m_index.clear();
for (unsigned i : m_ii.m_index) {
@ -326,7 +341,7 @@ void lu<T, X>::add_delta_to_solution_indexed(indexed_vector<T>& y) {
v = zero_of_type<T>();
}
lp_assert(y.is_OK());
SASSERT(y.is_OK());
}
template <typename T, typename X>
@ -343,7 +358,7 @@ void lu<T, X>::find_error_of_yB_indexed(const indexed_vector<T>& y, const vector
// it is a non efficient version
indexed_vector<T> yc = m_y_copy;
yc.m_index.clear();
lp_assert(!numeric_traits<T>::precise());
SASSERT(!numeric_traits<T>::precise());
{
vector<unsigned> d_basis(y.m_data.size());
@ -364,10 +379,10 @@ void lu<T, X>::find_error_of_yB_indexed(const indexed_vector<T>& y, const vector
}
}
#endif
lp_assert(m_ii.is_OK());
SASSERT(m_ii.is_OK());
m_ii.clear();
m_ii.resize(y.data_size());
lp_assert(m_y_copy.is_OK());
SASSERT(m_y_copy.is_OK());
// put the error into m_y_copy
for (auto k : y.m_index) {
auto & row = m_A.m_rows[k];
@ -399,7 +414,7 @@ void lu<T, X>::find_error_of_yB_indexed(const indexed_vector<T>& y, const vector
m_y_copy.set_value(v, k);
}
}
lp_assert(m_y_copy.is_OK());
SASSERT(m_y_copy.is_OK());
}
@ -419,12 +434,12 @@ void lu<T, X>::solve_yB_with_error_check_indexed(indexed_vector<T> & y, const ve
}
return;
}
lp_assert(m_y_copy.is_OK());
lp_assert(y.is_OK());
SASSERT(m_y_copy.is_OK());
SASSERT(y.is_OK());
if (y.m_index.size() * ratio_of_index_size_to_all_size<T>() < m_A.column_count()) {
m_y_copy = y;
solve_yB_indexed(y);
lp_assert(y.is_OK());
SASSERT(y.is_OK());
if (y.m_index.size() * ratio_of_index_size_to_all_size<T>() >= m_A.column_count()) {
find_error_of_yB(m_y_copy.m_data, y.m_data, basis);
solve_yB(m_y_copy.m_data);
@ -436,7 +451,7 @@ void lu<T, X>::solve_yB_with_error_check_indexed(indexed_vector<T> & y, const ve
solve_yB_indexed(m_y_copy);
add_delta_to_solution_indexed(y);
}
lp_assert(m_y_copy.is_OK());
SASSERT(m_y_copy.is_OK());
} else {
solve_yB_with_error_check(y.m_data, basis);
y.restore_index_and_clean_from_data();
@ -489,7 +504,7 @@ template <typename T, typename X>
void lu<T, X>::perform_transformations_on_w(indexed_vector<T>& w) {
apply_lp_list_to_w(w);
m_Q.apply_reverse_from_left(w);
// TBD does not compile: lp_assert(numeric_traits<T>::precise() || check_vector_for_small_values(w, m_settings));
// TBD does not compile: SASSERT(numeric_traits<T>::precise() || check_vector_for_small_values(w, m_settings));
}
// see Chvatal 24.3
@ -503,7 +518,7 @@ template <typename T, typename X>
void lu<T, X>::apply_lp_list_to_w(indexed_vector<T> & w) {
for (unsigned i = 0; i < m_tail.size(); i++) {
m_tail[i]->apply_from_left_to_T(w, m_settings);
// TBD does not compile: lp_assert(check_vector_for_small_values(w, m_settings));
// TBD does not compile: SASSERT(check_vector_for_small_values(w, m_settings));
}
}
template <typename T, typename X>
@ -570,7 +585,7 @@ unsigned lu<T, X>::transform_U_to_V_by_replacing_column(indexed_vector<T> & w,
return column_to_replace;
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X>
void lu<T, X>::check_vector_w(unsigned entering) {
T * w = new T[m_dim];
@ -595,7 +610,7 @@ void lu<T, X>::check_apply_lp_lists_to_w(T * w) {
permutation_matrix<T, X> qr = m_Q.get_reverse();
apply_to_vector(qr, w);
for (int i = m_dim - 1; i >= 0; i--) {
lp_assert(abs(w[i] - w[i]) < 0.0000001);
SASSERT(abs(w[i] - w[i]) < 0.0000001);
}
}
@ -624,7 +639,7 @@ void lu<T, X>::process_column(int j) {
}
template <typename T, typename X>
bool lu<T, X>::is_correct(const vector<unsigned>& basis) {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
if (get_status() != LU_status::OK) {
return false;
}
@ -637,10 +652,10 @@ bool lu<T, X>::is_correct(const vector<unsigned>& basis) {
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X>
dense_matrix<T, X> lu<T, X>::tail_product() {
lp_assert(tail_size() > 0);
SASSERT(tail_size() > 0);
dense_matrix<T, X> left_side = permutation_matrix<T, X>(m_dim);
for (unsigned i = 0; i < tail_size(); i++) {
matrix<T, X>* lp = get_lp_matrix(i);
@ -690,8 +705,8 @@ template <typename T, typename X>
bool lu<T, X>::all_columns_and_rows_are_active() {
unsigned i = m_dim;
while (i--) {
lp_assert(m_U.col_is_active(i));
lp_assert(m_U.row_is_active(i));
SASSERT(m_U.col_is_active(i));
SASSERT(m_U.row_is_active(i));
}
return true;
}
@ -733,9 +748,9 @@ void lu<T, X>::create_initial_factorization(){
}
}
if (j == m_dim) {
// TBD does not compile: lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
// lp_assert(is_correct());
// lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
// TBD does not compile: SASSERT(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
// SASSERT(is_correct());
// SASSERT(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
return;
}
j++;
@ -748,12 +763,12 @@ void lu<T, X>::create_initial_factorization(){
}
}
m_dense_LU->update_parent_matrix(m_settings);
lp_assert(m_dense_LU->is_L_matrix());
SASSERT(m_dense_LU->is_L_matrix());
m_dense_LU->conjugate_by_permutation(m_Q);
push_matrix_to_tail(m_dense_LU);
m_refactor_counter = 0;
// lp_assert(is_correct());
// lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
// SASSERT(is_correct());
// SASSERT(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
}
template <typename T, typename X>
@ -780,7 +795,7 @@ void lu<T, X>::scan_last_row_to_work_vector(unsigned lowest_row_of_the_bump) {
vector<indexed_value<T>> & last_row_vec = m_U.get_row_values(m_U.adjust_row(lowest_row_of_the_bump));
for (auto & iv : last_row_vec) {
if (is_zero(iv.m_value)) continue;
lp_assert(!m_settings.abs_val_is_smaller_than_drop_tolerance(iv.m_value));
SASSERT(!m_settings.abs_val_is_smaller_than_drop_tolerance(iv.m_value));
unsigned adjusted_col = m_U.adjust_column_inverse(iv.m_index);
if (adjusted_col < lowest_row_of_the_bump) {
m_row_eta_work_vector.set_value(-iv.m_value, adjusted_col);
@ -801,14 +816,14 @@ void lu<T, X>::pivot_and_solve_the_system(unsigned replaced_column, unsigned low
vector<indexed_value<T>> & row = m_U.get_row_values(aj);
for (auto & iv : row) {
unsigned col = m_U.adjust_column_inverse(iv.m_index);
lp_assert(col >= j || numeric_traits<T>::is_zero(iv.m_value));
SASSERT(col >= j || numeric_traits<T>::is_zero(iv.m_value));
if (col == j) continue;
if (numeric_traits<T>::is_zero(iv.m_value)) {
continue;
}
// the -v is for solving the system ( to zero the last row), and +v is for pivoting
T delta = col < lowest_row_of_the_bump? -v * iv.m_value: v * iv.m_value;
lp_assert(numeric_traits<T>::is_zero(delta) == false);
SASSERT(numeric_traits<T>::is_zero(delta) == false);
@ -845,7 +860,7 @@ row_eta_matrix<T, X> *lu<T, X>::get_row_eta_matrix_and_set_row_vector(unsigned r
return nullptr;
}
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
auto ret = new row_eta_matrix<T, X>(replaced_column, lowest_row_of_the_bump, m_dim);
#else
auto ret = new row_eta_matrix<T, X>(replaced_column, lowest_row_of_the_bump);
@ -885,15 +900,15 @@ void lu<T, X>::replace_column(T pivot_elem_for_checking, indexed_vector<T> & w,
push_matrix_to_tail(row_eta);
}
calculate_Lwave_Pwave_for_bump(replaced_column, lowest_row_of_the_bump);
// lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
// lp_assert(w.is_OK() && m_row_eta_work_vector.is_OK());
// SASSERT(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
// SASSERT(w.is_OK() && m_row_eta_work_vector.is_OK());
}
template <typename T, typename X>
void lu<T, X>::calculate_Lwave_Pwave_for_bump(unsigned replaced_column, unsigned lowest_row_of_the_bump){
T diagonal_elem;
if (replaced_column < lowest_row_of_the_bump) {
diagonal_elem = m_row_eta_work_vector[lowest_row_of_the_bump];
// lp_assert(m_row_eta_work_vector.is_OK());
// SASSERT(m_row_eta_work_vector.is_OK());
m_U.set_row_from_work_vector_and_clean_work_vector_not_adjusted(m_U.adjust_row(lowest_row_of_the_bump), m_row_eta_work_vector, m_settings);
} else {
diagonal_elem = m_U(lowest_row_of_the_bump, lowest_row_of_the_bump); // todo - get it more efficiently
@ -904,13 +919,13 @@ void lu<T, X>::calculate_Lwave_Pwave_for_bump(unsigned replaced_column, unsigned
}
calculate_Lwave_Pwave_for_last_row(lowest_row_of_the_bump, diagonal_elem);
// lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
// SASSERT(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
}
template <typename T, typename X>
void lu<T, X>::calculate_Lwave_Pwave_for_last_row(unsigned lowest_row_of_the_bump, T diagonal_element) {
auto l = new one_elem_on_diag<T, X>(lowest_row_of_the_bump, diagonal_element);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
l->set_number_of_columns(m_dim);
#endif
push_matrix_to_tail(l);
@ -927,11 +942,11 @@ void init_factorization(lu<T, X>* & factorization, static_matrix<T, X> & m_A, ve
// LP_OUT(m_settings, "failing in init_factorization" << std::endl);
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X>
dense_matrix<T, X> get_B(lu<T, X>& f, const vector<unsigned>& basis) {
lp_assert(basis.size() == f.dimension());
lp_assert(basis.size() == f.m_U.dimension());
SASSERT(basis.size() == f.dimension());
SASSERT(basis.size() == f.m_U.dimension());
dense_matrix<T, X> B(f.dimension(), f.dimension());
for (unsigned i = 0; i < f.dimension(); i++)
for (unsigned j = 0; j < f.dimension(); j++)

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <utility>
#include <memory>
#include <string>
@ -24,7 +39,7 @@ template lp::mpq lp::dot_product<lp::mpq, lp::mpq>(vector<lp::mpq > const&, vect
template void lp::init_factorization<double, double>(lp::lu<double, double>*&, lp::static_matrix<double, double>&, vector<unsigned int>&, lp::lp_settings&);
template void lp::init_factorization<lp::mpq, lp::mpq>(lp::lu<lp::mpq, lp::mpq>*&, lp::static_matrix<lp::mpq, lp::mpq>&, vector<unsigned int>&, lp::lp_settings&);
template void lp::init_factorization<lp::mpq, lp::numeric_pair<lp::mpq> >(lp::lu<lp::mpq, lp::numeric_pair<lp::mpq> >*&, lp::static_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >&, vector<unsigned int>&, lp::lp_settings&);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template void lp::print_matrix<double, double>(lp::sparse_matrix<double, double>&, std::ostream & out);
template void lp::print_matrix<lp::mpq, lp::mpq>(lp::static_matrix<lp::mpq, lp::mpq>&, std::ostream&);
template void lp::print_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >(lp::static_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >&, std::ostream&);

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#ifdef Z3DEBUG
#pragma once
#include "util/lp/numeric_pair.h"

23
src/util/lp/matrix.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#ifdef Z3DEBUG
#include <cmath>

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

@ -1,9 +1,24 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/lp_settings.h"
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
#include "util/lp/matrix.hpp"
#include "util/lp/static_matrix.h"
#include <string>

49
src/util/lp/mps_reader.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
@ -160,9 +175,9 @@ class mps_reader {
if (m_line[i] == ' ')
break;
}
lp_assert(m_line.size() >= offset);
lp_assert(m_line.size() >> i);
lp_assert(i >= offset);
SASSERT(m_line.size() >= offset);
SASSERT(m_line.size() >> i);
SASSERT(i >= offset);
return m_line.substr(offset, i - offset);
}
@ -497,7 +512,7 @@ class mps_reader {
void create_or_update_bound() {
const unsigned name_offset = 14;
lp_assert(m_line.size() >= 14);
SASSERT(m_line.size() >= 14);
vector<std::string> bound_string = split_and_trim(m_line.substr(name_offset, m_line.size()));
if (bound_string.size() == 0) {
@ -603,7 +618,7 @@ class mps_reader {
}
for (auto s : row_with_range->m_row_columns) {
lp_assert(m_columns.find(s.first) != m_columns.end());
SASSERT(m_columns.find(s.first) != m_columns.end());
other_bound_range_row->m_row_columns[s.first] = s.second;
}
}
@ -679,7 +694,7 @@ class mps_reader {
if (row->m_name != m_cost_row_name) {
solver->add_constraint(get_relation_from_row(row->m_type), row->m_right_side, row->m_index);
for (auto s : row->m_row_columns) {
lp_assert(m_columns.find(s.first) != m_columns.end());
SASSERT(m_columns.find(s.first) != m_columns.end());
solver->set_row_column_coefficient(row->m_index, m_columns[s.first]->m_index, s.second);
}
} else {
@ -714,7 +729,7 @@ class mps_reader {
void set_solver_cost(row * row, lp_solver<T, X> *solver) {
for (auto s : row->m_row_columns) {
std::string name = s.first;
lp_assert(m_columns.find(name) != m_columns.end());
SASSERT(m_columns.find(name) != m_columns.end());
mps_reader::column * col = m_columns[name];
solver->set_cost_for_column(col->m_index, s.second);
}
@ -723,7 +738,7 @@ class mps_reader {
public:
void set_message_stream(std::ostream * o) {
lp_assert(o != nullptr);
SASSERT(o != nullptr);
m_message_stream = o;
}
vector<std::string> column_names() {
@ -810,7 +825,7 @@ public:
auto kind = get_lar_relation_from_row(row->m_type);
vector<std::pair<mpq, var_index>> ls;
for (auto s : row->m_row_columns) {
var_index i = solver->add_var(get_var_index(s.first), false);
var_index i = solver->add_var(get_var_index(s.first));
ls.push_back(std::make_pair(s.second, i));
}
solver->add_constraint(ls, kind, row->m_right_side);
@ -828,20 +843,20 @@ public:
void create_low_constraint_for_var(column* col, bound * b, lar_solver *solver) {
vector<std::pair<mpq, var_index>> ls;
var_index i = solver->add_var(col->m_index, false);
var_index i = solver->add_var(col->m_index);
ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
solver->add_constraint(ls, GE, b->m_low);
}
void create_upper_constraint_for_var(column* col, bound * b, lar_solver *solver) {
var_index i = solver->add_var(col->m_index, false);
var_index i = solver->add_var(col->m_index);
vector<std::pair<mpq, var_index>> ls;
ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
solver->add_constraint(ls, LE, b->m_upper);
}
void create_equality_contraint_for_var(column* col, bound * b, lar_solver *solver) {
var_index i = solver->add_var(col->m_index, false);
var_index i = solver->add_var(col->m_index);
vector<std::pair<mpq, var_index>> ls;
ls.push_back(std::make_pair(numeric_traits<T>::one(), i));
solver->add_constraint(ls, EQ, b->m_fixed_value);
@ -850,7 +865,7 @@ public:
void fill_lar_solver_on_columns(lar_solver * solver) {
for (auto s : m_columns) {
mps_reader::column * col = s.second;
solver->add_var(col->m_index, false);
solver->add_var(col->m_index);
auto b = col->m_bound;
if (b == nullptr) return;

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

@ -166,7 +166,6 @@ namespace nra {
polynomial::polynomial* ps[1] = { p };
bool is_even[1] = { false };
nlsat::literal lit;
nlsat::assumption a = this + idx;
switch (k) {
case lp::lconstraint_kind::LE:
@ -262,5 +261,4 @@ namespace nra {
return m_imp->am();
}
}

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

@ -1,33 +1,38 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
The idea is that it is only one different file in Lean and z3 source inside of LP
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#define lp_for_z3
#include <string>
#include <cmath>
#include <algorithm>
#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 lp {
#ifdef lp_for_z3 // rename rationals
typedef rational mpq;
#else
typedef lp::mpq mpq;
#endif
typedef rational mpq; // rename rationals
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<unsigned> {
@ -67,14 +72,13 @@ template <> class numeric_traits<double> {
static bool is_pos(const rational & d) {return d.is_pos();}
static bool is_neg(const rational & d) {return d.is_neg();}
};
#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(numeric_traits<X>::get_double(x)) < y; }
static bool below_bound_numeric(const X &, const X &, const Y &) { /*lp_unreachable();*/ return false;}
static bool above_bound_numeric(const X &, const X &, const Y &) { /*lp_unreachable();*/ return false; }
static bool below_bound_numeric(const X &, const X &, const Y &) { /*SASSERT(false);*/ return false;}
static bool above_bound_numeric(const X &, const X &, const Y &) { /*SASSERT(false);*/ return false; }
};
@ -104,9 +108,9 @@ struct numeric_pair {
template <typename X>
numeric_pair(const X & n) : x(n), y(0) {
}
numeric_pair(const numeric_pair<T> & n) : x(n.x), y(n.y) {}
template <typename X, typename Y>
numeric_pair(X xp, Y yp) : x(convert_struct<T, X>::convert(xp)), y(convert_struct<T, Y>::convert(yp)) {}
@ -144,16 +148,16 @@ struct numeric_pair {
}
numeric_pair operator/(const numeric_pair &) const {
// lp_unreachable();
// SASSERT(false);
}
numeric_pair operator+(const numeric_pair & a) const {
return numeric_pair(a.x + x, a.y + y);
}
numeric_pair operator*(const numeric_pair & /*a*/) const {
// lp_unreachable();
// SASSERT(false);
}
numeric_pair& operator+=(const numeric_pair & a) {
@ -195,15 +199,10 @@ struct numeric_pair {
bool is_pos() const { return x.is_pos() || (x.is_zero() && y.is_pos());}
bool is_neg() const { return x.is_neg() || (x.is_zero() && y.is_neg());}
std::string to_string() const {
return std::string("(") + T_to_string(x) + ", " + T_to_string(y) + ")";
}
bool is_int() const {
return x.is_int() && y.is_zero();
}
};
@ -230,7 +229,7 @@ numeric_pair<T> operator/(const numeric_pair<T> & r, const X & a) {
}
// template <numeric_pair, typename T> bool precise() { return numeric_traits<T>::precise();}
template <typename T> double get_double(const lp::numeric_pair<T> & ) { /* lp_unreachable(); */ return 0;}
template <typename T> double get_double(const lp::numeric_pair<T> & ) { /* SASSERT(false); */ return 0;}
template <typename T>
class numeric_traits<lp::numeric_pair<T>> {
public:
@ -238,7 +237,7 @@ class numeric_traits<lp::numeric_pair<T>> {
static lp::numeric_pair<T> zero() { return lp::numeric_pair<T>(numeric_traits<T>::zero(), numeric_traits<T>::zero()); }
static bool is_zero(const lp::numeric_pair<T> & v) { return numeric_traits<T>::is_zero(v.x) && numeric_traits<T>::is_zero(v.y); }
static double get_double(const lp::numeric_pair<T> & v){ return numeric_traits<T>::get_double(v.x); } // just return the double of the first coordinate
static double one() { /*lp_unreachable();*/ return 0;}
static double one() { /*SASSERT(false);*/ return 0;}
static bool is_pos(const numeric_pair<T> &p) {
return numeric_traits<T>::is_pos(p.x) ||
(numeric_traits<T>::is_zero(p.x) && numeric_traits<T>::is_pos(p.y));
@ -247,7 +246,7 @@ class numeric_traits<lp::numeric_pair<T>> {
return numeric_traits<T>::is_neg(p.x) ||
(numeric_traits<T>::is_zero(p.x) && numeric_traits<T>::is_neg(p.y));
}
};
template <>
@ -268,11 +267,11 @@ struct convert_struct<numeric_pair<T>, double> {
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 &) {
// lp_unreachable();
// SASSERT(false);
return false;
}
static bool above_bound_numeric(const numeric_pair<T> &, const numeric_pair<T> &, const double &) {
// lp_unreachable();
// SASSERT(false);
return false;
}
};
@ -329,26 +328,4 @@ struct convert_struct<double, double> {
template <typename X> bool is_epsilon_small(const X & v, const double &eps) { return convert_struct<X, double>::is_epsilon_small(v, eps);}
template <typename X> bool below_bound_numeric(const X & x, const X & bound, const double& eps) { return convert_struct<X, double>::below_bound_numeric(x, bound, eps);}
template <typename X> bool above_bound_numeric(const X & x, const X & bound, const double& eps) { return convert_struct<X, double>::above_bound_numeric(x, bound, eps);}
template <typename T> T floor(const numeric_pair<T> & r) {
if (r.x.is_int()) {
if (r.y.is_nonneg()) {
return r.x;
}
return r.x - mpq::one();
}
return floor(r.x);
}
template <typename T> T ceil(const numeric_pair<T> & r) {
if (r.x.is_int()) {
if (r.y.is_nonpos()) {
return r.x;
}
return r.x + mpq::one();
}
return ceil(r.x);
}
}

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
#include <algorithm>
@ -13,7 +28,7 @@
#include "util/lp/matrix.h"
#include "util/lp/tail_matrix.h"
namespace lp {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
inline bool is_even(int k) { return (k/2)*2 == k; }
#endif
@ -50,7 +65,7 @@ class permutation_matrix : public tail_matrix<T, X> {
void init(unsigned length);
unsigned get_rev(unsigned i) { return m_rev[i]; }
bool is_dense() const { return false; }
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
permutation_matrix get_inverse() const {
return permutation_matrix(size(), m_rev);
}
@ -86,14 +101,14 @@ class permutation_matrix : public tail_matrix<T, X> {
void apply_reverse_from_right_to_X(vector<X> & w);
void set_val(unsigned i, unsigned pi) {
lp_assert(i < size() && pi < size()); m_permutation[i] = pi; m_rev[pi] = i; }
SASSERT(i < size() && pi < size()); m_permutation[i] = pi; m_rev[pi] = i; }
void transpose_from_left(unsigned i, unsigned j);
unsigned apply_reverse(unsigned i) const { return m_rev[i]; }
void transpose_from_right(unsigned i, unsigned j);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
T get_elem(unsigned i, unsigned j) const{
return m_permutation[i] == j? numeric_traits<T>::one() : numeric_traits<T>::zero();
}

81
src/util/lp/permutation_matrix.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/vector.h"
#include "util/lp/permutation_matrix.h"
namespace lp {
@ -27,7 +42,7 @@ template <typename T, typename X> void permutation_matrix<T, X>::init(unsigned l
}
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X> void permutation_matrix<T, X>::print(std::ostream & out) const {
out << "[";
for (unsigned i = 0; i < size(); i++) {
@ -44,13 +59,13 @@ template <typename T, typename X> void permutation_matrix<T, X>::print(std::ostr
template <typename T, typename X>
void permutation_matrix<T, X>::apply_from_left(vector<X> & w, lp_settings & ) {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// dense_matrix<L, X> deb(*this);
// L * deb_w = clone_vector<L>(w, row_count());
// deb.apply_from_left(deb_w);
#endif
// std::cout << " apply_from_left " << std::endl;
lp_assert(m_X_buffer.size() == w.size());
SASSERT(m_X_buffer.size() == w.size());
unsigned i = size();
while (i-- > 0) {
m_X_buffer[i] = w[m_permutation[i]];
@ -59,8 +74,8 @@ void permutation_matrix<T, X>::apply_from_left(vector<X> & w, lp_settings & ) {
while (i-- > 0) {
w[i] = m_X_buffer[i];
}
#ifdef LEAN_DEBUG
// lp_assert(vectors_are_equal<L>(deb_w, w, row_count()));
#ifdef Z3DEBUG
// SASSERT(vectors_are_equal<L>(deb_w, w, row_count()));
// delete [] deb_w;
#endif
}
@ -81,12 +96,12 @@ void permutation_matrix<T, X>::apply_from_left_to_T(indexed_vector<T> & w, lp_se
}
template <typename T, typename X> void permutation_matrix<T, X>::apply_from_right(vector<T> & w) {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// dense_matrix<T, X> deb(*this);
// T * deb_w = clone_vector<T>(w, row_count());
// deb.apply_from_right(deb_w);
#endif
lp_assert(m_T_buffer.size() == w.size());
SASSERT(m_T_buffer.size() == w.size());
for (unsigned i = 0; i < size(); i++) {
m_T_buffer[i] = w[m_rev[i]];
}
@ -94,14 +109,14 @@ template <typename T, typename X> void permutation_matrix<T, X>::apply_from_righ
for (unsigned i = 0; i < size(); i++) {
w[i] = m_T_buffer[i];
}
#ifdef LEAN_DEBUG
// lp_assert(vectors_are_equal<T>(deb_w, w, row_count()));
#ifdef Z3DEBUG
// SASSERT(vectors_are_equal<T>(deb_w, w, row_count()));
// delete [] deb_w;
#endif
}
template <typename T, typename X> void permutation_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
vector<T> wcopy(w.m_data);
apply_from_right(wcopy);
#endif
@ -117,9 +132,9 @@ template <typename T, typename X> void permutation_matrix<T, X>::apply_from_righ
unsigned pj = m_permutation[j];
w.set_value(buffer[i], pj);
}
lp_assert(w.is_OK());
#ifdef LEAN_DEBUG
lp_assert(vectors_are_equal(wcopy, w.m_data));
SASSERT(w.is_OK());
#ifdef Z3DEBUG
SASSERT(vectors_are_equal(wcopy, w.m_data));
#endif
}
@ -147,7 +162,7 @@ void permutation_matrix<T, X>::clear_data(indexed_vector<L> & w) {
template <typename T, typename X>template <typename L>
void permutation_matrix<T, X>::apply_reverse_from_left(indexed_vector<L> & w) {
// the result will be w = p(-1) * w
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// dense_matrix<L, X> deb(get_reverse());
// L * deb_w = clone_vector<L>(w.m_data, row_count());
// deb.apply_from_left(deb_w);
@ -165,8 +180,8 @@ void permutation_matrix<T, X>::apply_reverse_from_left(indexed_vector<L> & w) {
w[j] = t[i];
w.m_index[i] = j;
}
#ifdef LEAN_DEBUG
// lp_assert(vectors_are_equal<L>(deb_w, w.m_data, row_count()));
#ifdef Z3DEBUG
// SASSERT(vectors_are_equal<L>(deb_w, w.m_data, row_count()));
// delete [] deb_w;
#endif
}
@ -174,7 +189,7 @@ void permutation_matrix<T, X>::apply_reverse_from_left(indexed_vector<L> & w) {
template <typename T, typename X>
void permutation_matrix<T, X>::apply_reverse_from_left_to_T(vector<T> & w) {
// the result will be w = p(-1) * w
lp_assert(m_T_buffer.size() == w.size());
SASSERT(m_T_buffer.size() == w.size());
unsigned i = size();
while (i-- > 0) {
m_T_buffer[m_permutation[i]] = w[i];
@ -187,7 +202,7 @@ void permutation_matrix<T, X>::apply_reverse_from_left_to_T(vector<T> & w) {
template <typename T, typename X>
void permutation_matrix<T, X>::apply_reverse_from_left_to_X(vector<X> & w) {
// the result will be w = p(-1) * w
lp_assert(m_X_buffer.size() == w.size());
SASSERT(m_X_buffer.size() == w.size());
unsigned i = size();
while (i-- > 0) {
m_X_buffer[m_permutation[i]] = w[i];
@ -201,7 +216,7 @@ void permutation_matrix<T, X>::apply_reverse_from_left_to_X(vector<X> & w) {
template <typename T, typename X>
void permutation_matrix<T, X>::apply_reverse_from_right_to_T(vector<T> & w) {
// the result will be w = w * p(-1)
lp_assert(m_T_buffer.size() == w.size());
SASSERT(m_T_buffer.size() == w.size());
unsigned i = size();
while (i-- > 0) {
m_T_buffer[i] = w[m_permutation[i]];
@ -215,11 +230,11 @@ void permutation_matrix<T, X>::apply_reverse_from_right_to_T(vector<T> & w) {
template <typename T, typename X>
void permutation_matrix<T, X>::apply_reverse_from_right_to_T(indexed_vector<T> & w) {
// the result will be w = w * p(-1)
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// vector<T> wcopy(w.m_data);
// apply_reverse_from_right_to_T(wcopy);
#endif
lp_assert(w.is_OK());
SASSERT(w.is_OK());
vector<T> tmp;
vector<unsigned> tmp_index(w.m_index);
for (auto i : w.m_index) {
@ -232,15 +247,15 @@ void permutation_matrix<T, X>::apply_reverse_from_right_to_T(indexed_vector<T> &
w.set_value(tmp[k], m_rev[j]);
}
// lp_assert(w.is_OK());
// lp_assert(vectors_are_equal(w.m_data, wcopy));
// SASSERT(w.is_OK());
// SASSERT(vectors_are_equal(w.m_data, wcopy));
}
template <typename T, typename X>
void permutation_matrix<T, X>::apply_reverse_from_right_to_X(vector<X> & w) {
// the result will be w = w * p(-1)
lp_assert(m_X_buffer.size() == w.size());
SASSERT(m_X_buffer.size() == w.size());
unsigned i = size();
while (i-- > 0) {
m_X_buffer[i] = w[m_permutation[i]];
@ -253,7 +268,7 @@ void permutation_matrix<T, X>::apply_reverse_from_right_to_X(vector<X> & w) {
template <typename T, typename X> void permutation_matrix<T, X>::transpose_from_left(unsigned i, unsigned j) {
// the result will be this = (i,j)*this
lp_assert(i < size() && j < size() && i != j);
SASSERT(i < size() && j < size() && i != j);
auto pi = m_rev[i];
auto pj = m_rev[j];
set_val(pi, j);
@ -262,7 +277,7 @@ template <typename T, typename X> void permutation_matrix<T, X>::transpose_from_
template <typename T, typename X> void permutation_matrix<T, X>::transpose_from_right(unsigned i, unsigned j) {
// the result will be this = this * (i,j)
lp_assert(i < size() && j < size() && i != j);
SASSERT(i < size() && j < size() && i != j);
auto pi = m_permutation[i];
auto pj = m_permutation[j];
set_val(i, pj);
@ -271,7 +286,7 @@ template <typename T, typename X> void permutation_matrix<T, X>::transpose_from_
template <typename T, typename X> void permutation_matrix<T, X>::multiply_by_permutation_from_left(permutation_matrix<T, X> & p) {
m_work_array = m_permutation;
lp_assert(p.size() == size());
SASSERT(p.size() == size());
unsigned i = size();
while (i-- > 0) {
set_val(i, m_work_array[p[i]]); // we have m(P)*m(Q) = m(QP), where m is the matrix of the permutation
@ -281,7 +296,7 @@ template <typename T, typename X> void permutation_matrix<T, X>::multiply_by_per
// this is multiplication in the matrix sense
template <typename T, typename X> void permutation_matrix<T, X>::multiply_by_permutation_from_right(permutation_matrix<T, X> & p) {
m_work_array = m_permutation;
lp_assert(p.size() == size());
SASSERT(p.size() == size());
unsigned i = size();
while (i-- > 0)
set_val(i, p[m_work_array[i]]); // we have m(P)*m(Q) = m(QP), where m is the matrix of the permutation
@ -289,7 +304,7 @@ template <typename T, typename X> void permutation_matrix<T, X>::multiply_by_per
}
template <typename T, typename X> void permutation_matrix<T, X>::multiply_by_reverse_from_right(permutation_matrix<T, X> & q){ // todo : condensed permutations ?
lp_assert(q.size() == size());
SASSERT(q.size() == size());
m_work_array = m_permutation;
// the result is this = this*q(-1)
unsigned i = size();

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <memory>
#include "util/vector.h"
#include "util/lp/permutation_matrix.hpp"

53
src/util/lp/quick_xplain.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/lar_solver.h"
namespace lp {
quick_xplain::quick_xplain(vector<std::pair<mpq, constraint_index>> & explanation, const lar_solver & ls, lar_solver & qsol) :
@ -15,12 +30,12 @@ void quick_xplain::add_constraint_to_qsol(unsigned j) {
auto ci = m_qsol.add_constraint(ls, lar_c.m_kind, lar_c.m_right_side);
m_local_ci_to_constraint_offsets[ci] = j;
}
void quick_xplain::copy_constraint_and_add_constraint_vars(const lar_constraint& lar_c) {
vector < std::pair<mpq, unsigned>> ls;
for (auto & p : lar_c.get_left_side_coefficients()) {
unsigned j = p.second;
unsigned lj = m_qsol.add_var(j, false);
unsigned lj = m_qsol.add_var(j);
ls.push_back(std::make_pair(p.first, lj));
}
m_constraints_in_local_vars.push_back(lar_constraint(ls, lar_c.m_kind, lar_c.m_right_side));
@ -29,7 +44,7 @@ void quick_xplain::copy_constraint_and_add_constraint_vars(const lar_constraint&
bool quick_xplain::infeasible() {
m_qsol.solve();
return m_qsol.get_status() == lp_status::INFEASIBLE;
return m_qsol.get_status() == INFEASIBLE;
}
// u - unexplored constraints
@ -56,9 +71,9 @@ void quick_xplain::minimize(const vector<unsigned>& u) {
}
}
if (m > 0) {
lp_assert(m_qsol.constraint_stack_size() >= initial_stack_size);
SASSERT(m_qsol.constraint_stack_size() >= initial_stack_size);
m_qsol.pop(m_qsol.constraint_stack_size() - initial_stack_size);
for (auto j : m_x)
for (auto j : m_x)
add_constraint_to_qsol(j);
if (!infeasible()) {
vector<unsigned> un;
@ -69,11 +84,11 @@ void quick_xplain::minimize(const vector<unsigned>& u) {
}
}
void quick_xplain::run(vector<std::pair<mpq, constraint_index>> & explanation, const lar_solver & ls){
if (explanation.size() <= 2) return;
lar_solver qsol;
lp_assert(ls.explanation_is_correct(explanation));
SASSERT(ls.explanation_is_correct(explanation));
quick_xplain q(explanation, ls, qsol);
q.solve();
}
@ -94,13 +109,13 @@ bool quick_xplain::is_feasible(const vector<unsigned> & x, unsigned k) const {
vector < std::pair<mpq, unsigned>> ls;
const lar_constraint & c = m_constraints_in_local_vars[i];
for (auto & p : c.get_left_side_coefficients()) {
unsigned lj = l.add_var(p.second, false);
unsigned lj = l.add_var(p.second);
ls.push_back(std::make_pair(p.first, lj));
}
l.add_constraint(ls, c.m_kind, c.m_right_side);
}
l.solve();
return l.get_status() != lp_status::INFEASIBLE;
return l.get_status() != INFEASIBLE;
}
bool quick_xplain::x_is_minimal() const {
@ -109,7 +124,7 @@ bool quick_xplain::x_is_minimal() const {
x.push_back(j);
for (unsigned k = 0; k < x.size(); k++) {
lp_assert(is_feasible(x, x[k]));
SASSERT(is_feasible(x, x[k]));
}
return true;
}
@ -117,8 +132,8 @@ bool quick_xplain::x_is_minimal() const {
void quick_xplain::solve() {
copy_constraints_to_local_constraints();
m_qsol.push();
lp_assert(m_qsol.constraint_count() == 0)
vector<unsigned> u;
SASSERT(m_qsol.constraint_count() == 0);
vector<unsigned> u;
for (unsigned k = 0; k < m_constraints_in_local_vars.size(); k++)
u.push_back(k);
minimize(u);
@ -127,10 +142,10 @@ void quick_xplain::solve() {
for (unsigned i : m_x)
add_constraint_to_qsol(i);
m_qsol.solve();
lp_assert(m_qsol.get_status() == lp_status::INFEASIBLE);
SASSERT(m_qsol.get_status() == INFEASIBLE);
m_qsol.get_infeasibility_explanation(m_explanation);
lp_assert(m_qsol.explanation_is_correct(m_explanation));
lp_assert(x_is_minimal());
SASSERT(m_qsol.explanation_is_correct(m_explanation));
SASSERT(x_is_minimal());
for (auto & p : m_explanation) {
p.second = this->m_local_constraint_offset_to_external_ci[m_local_ci_to_constraint_offsets[p.second]];
}

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

@ -1,7 +1,22 @@
/*
/*++
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"

22
src/util/lp/random_updater.h
View file

@ -1,7 +1,23 @@
/*
/*++
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <set>
#include "util/vector.h"

59
src/util/lp/random_updater.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/random_updater.h"
#include "util/lp/static_matrix.h"
#include "util/lp/lar_solver.h"
@ -36,7 +51,7 @@ random_updater::interval random_updater::get_interval_of_non_basic_var(unsigned
ret.set_upper_bound(m_core_solver.m_r_upper_bounds[j]);
break;
default:
lp_assert(false);
SASSERT(false);
}
return ret;
}
@ -44,15 +59,15 @@ random_updater::interval random_updater::get_interval_of_non_basic_var(unsigned
void random_updater::diminish_interval_for_basic_var(numeric_pair<mpq>& nb_x, unsigned j,
mpq & a,
interval & r) {
lp_assert(m_core_solver.m_r_heading[j] >= 0);
SASSERT(m_core_solver.m_r_heading[j] >= 0);
numeric_pair<mpq> delta;
lp_assert(a != zero_of_type<mpq>());
SASSERT(a != zero_of_type<mpq>());
switch (m_core_solver.get_column_type(j)) {
case column_type::free_column:
break;
case column_type::low_bound:
delta = m_core_solver.m_r_x[j] - m_core_solver.m_r_low_bounds[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
SASSERT(delta >= zero_of_type<numeric_pair<mpq>>());
if (a > 0) {
r.set_upper_bound(nb_x + delta / a);
} else {
@ -61,7 +76,7 @@ void random_updater::diminish_interval_for_basic_var(numeric_pair<mpq>& nb_x, un
break;
case column_type::upper_bound:
delta = m_core_solver.m_r_upper_bounds()[j] - m_core_solver.m_r_x[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
SASSERT(delta >= zero_of_type<numeric_pair<mpq>>());
if (a > 0) {
r.set_low_bound(nb_x - delta / a);
} else {
@ -71,17 +86,17 @@ void random_updater::diminish_interval_for_basic_var(numeric_pair<mpq>& nb_x, un
case column_type::boxed:
if (a > 0) {
delta = m_core_solver.m_r_x[j] - m_core_solver.m_r_low_bounds[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
SASSERT(delta >= zero_of_type<numeric_pair<mpq>>());
r.set_upper_bound(nb_x + delta / a);
delta = m_core_solver.m_r_upper_bounds()[j] - m_core_solver.m_r_x[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
SASSERT(delta >= zero_of_type<numeric_pair<mpq>>());
r.set_low_bound(nb_x - delta / a);
} else { // a < 0
delta = m_core_solver.m_r_upper_bounds()[j] - m_core_solver.m_r_x[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
SASSERT(delta >= zero_of_type<numeric_pair<mpq>>());
r.set_upper_bound(nb_x - delta / a);
delta = m_core_solver.m_r_x[j] - m_core_solver.m_r_low_bounds[j];
lp_assert(delta >= zero_of_type<numeric_pair<mpq>>());
SASSERT(delta >= zero_of_type<numeric_pair<mpq>>());
r.set_low_bound(nb_x + delta / a);
}
break;
@ -90,7 +105,7 @@ void random_updater::diminish_interval_for_basic_var(numeric_pair<mpq>& nb_x, un
r.set_upper_bound(nb_x);
break;
default:
lp_assert(false);
SASSERT(false);
}
}
@ -113,15 +128,15 @@ random_updater::interval random_updater::find_shift_interval(unsigned j) {
}
void random_updater::shift_var(unsigned j, interval & r) {
lp_assert(r.contains(m_core_solver.m_r_x[j]));
lp_assert(m_core_solver.m_r_solver.column_is_feasible(j));
SASSERT(r.contains(m_core_solver.m_r_x[j]));
SASSERT(m_core_solver.m_r_solver.column_is_feasible(j));
auto old_x = m_core_solver.m_r_x[j];
remove_value(old_x);
auto new_val = m_core_solver.m_r_x[j] = get_random_from_interval(r);
add_value(new_val);
lp_assert(r.contains(m_core_solver.m_r_x[j]));
lp_assert(m_core_solver.m_r_solver.column_is_feasible(j));
SASSERT(r.contains(m_core_solver.m_r_x[j]));
SASSERT(m_core_solver.m_r_solver.column_is_feasible(j));
auto delta = m_core_solver.m_r_x[j] - old_x;
unsigned i;
@ -130,9 +145,9 @@ void random_updater::shift_var(unsigned j, interval & r) {
while(m_column_j->next(a, i)) {
unsigned bj = m_core_solver.m_r_basis[i];
m_core_solver.m_r_x[bj] -= a * delta;
lp_assert(m_core_solver.m_r_solver.column_is_feasible(bj));
SASSERT(m_core_solver.m_r_solver.column_is_feasible(bj));
}
lp_assert(m_core_solver.m_r_solver.A_mult_x_is_off() == false);
SASSERT(m_core_solver.m_r_solver.A_mult_x_is_off() == false);
}
numeric_pair<mpq> random_updater::get_random_from_interval(interval & r) {
@ -143,7 +158,7 @@ numeric_pair<mpq> random_updater::get_random_from_interval(interval & r) {
return r.low_bound + numeric_pair<mpq>(rand % range, 0);
if ((!r.low_bound_is_set) && r.upper_bound_is_set)
return r.upper_bound - numeric_pair<mpq>(rand % range, 0);
lp_assert(r.low_bound_is_set && r.upper_bound_is_set);
SASSERT(r.low_bound_is_set && r.upper_bound_is_set);
return r.low_bound + (rand % range) * (r.upper_bound - r.low_bound)/ range;
}
@ -183,7 +198,7 @@ void random_updater::add_value(numeric_pair<mpq>& v) {
void random_updater::remove_value(numeric_pair<mpq>& v) {
std::unordered_map<numeric_pair<mpq>, unsigned>::iterator it = m_values.find(v);
lp_assert(it != m_values.end());
SASSERT(it != m_values.end());
it->second--;
if (it->second == 0)
m_values.erase((std::unordered_map<numeric_pair<mpq>, unsigned>::const_iterator)it);

23
src/util/lp/random_updater_instances.cpp
View file

@ -1,5 +1,20 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/random_updater.hpp"

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -16,20 +31,20 @@ namespace lp {
template <typename T, typename X>
class row_eta_matrix
: public tail_matrix<T, X> {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
unsigned m_dimension;
#endif
unsigned m_row_start;
unsigned m_row;
sparse_vector<T> m_row_vector;
public:
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
row_eta_matrix(unsigned row_start, unsigned row, unsigned dim):
#else
row_eta_matrix(unsigned row_start, unsigned row):
#endif
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
m_dimension(dim),
#endif
m_row_start(row_start), m_row(row) {
@ -55,7 +70,7 @@ public:
}
void push_back(unsigned row_index, T val ) {
lp_assert(row_index != m_row);
SASSERT(row_index != m_row);
m_row_vector.push_back(row_index, val);
}
@ -63,7 +78,7 @@ public:
void apply_from_right(indexed_vector<T> & w);
void conjugate_by_permutation(permutation_matrix<T, X> & p);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
T get_elem(unsigned row, unsigned col) const;
unsigned row_count() const { return m_dimension; }
unsigned column_count() const { return m_dimension; }

55
src/util/lp/row_eta_matrix.hpp
View file

@ -1,13 +1,28 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/vector.h"
#include "util/lp/row_eta_matrix.h"
namespace lp {
template <typename T, typename X>
void row_eta_matrix<T, X>::apply_from_left(vector<X> & w, lp_settings &) {
// #ifdef LEAN_DEBUG
// #ifdef Z3DEBUG
// dense_matrix<T> deb(*this);
// auto clone_w = clone_vector<T>(w, m_dimension);
// deb.apply_from_left(clone_w, settings);
@ -18,8 +33,8 @@ void row_eta_matrix<T, X>::apply_from_left(vector<X> & w, lp_settings &) {
w_at_row += w[it.first] * it.second;
}
// w[m_row] = w_at_row;
// #ifdef LEAN_DEBUG
// lp_assert(vectors_are_equal<T>(clone_w, w, m_dimension));
// #ifdef Z3DEBUG
// SASSERT(vectors_are_equal<T>(clone_w, w, m_dimension));
// delete [] clone_w;
// #endif
}
@ -43,7 +58,7 @@ void row_eta_matrix<T, X>::apply_from_left_local_to_T(indexed_vector<T> & w, lp_
auto it = std::find(w.m_index.begin(), w.m_index.end(), m_row);
w.m_index.erase(it);
}
// TBD: lp_assert(check_vector_for_small_values(w, settings));
// TBD: SASSERT(check_vector_for_small_values(w, settings));
}
template <typename T, typename X>
@ -65,14 +80,14 @@ void row_eta_matrix<T, X>::apply_from_left_local_to_X(indexed_vector<X> & w, lp_
auto it = std::find(w.m_index.begin(), w.m_index.end(), m_row);
w.m_index.erase(it);
}
// TBD: does not compile lp_assert(check_vector_for_small_values(w, settings));
// TBD: does not compile SASSERT(check_vector_for_small_values(w, settings));
}
template <typename T, typename X>
void row_eta_matrix<T, X>::apply_from_right(vector<T> & w) {
const T & w_row = w[m_row];
if (numeric_traits<T>::is_zero(w_row)) return;
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// dense_matrix<T> deb(*this);
// auto clone_w = clone_vector<T>(w, m_dimension);
// deb.apply_from_right(clone_w);
@ -80,18 +95,18 @@ void row_eta_matrix<T, X>::apply_from_right(vector<T> & w) {
for (auto & it : m_row_vector.m_data) {
w[it.first] += w_row * it.second;
}
#ifdef LEAN_DEBUG
// lp_assert(vectors_are_equal<T>(clone_w, w, m_dimension));
#ifdef Z3DEBUG
// SASSERT(vectors_are_equal<T>(clone_w, w, m_dimension));
// delete clone_w;
#endif
}
template <typename T, typename X>
void row_eta_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
lp_assert(w.is_OK());
SASSERT(w.is_OK());
const T & w_row = w[m_row];
if (numeric_traits<T>::is_zero(w_row)) return;
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// vector<T> wcopy(w.m_data);
// apply_from_right(wcopy);
#endif
@ -129,8 +144,8 @@ void row_eta_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
}
}
}
#ifdef LEAN_DEBUG
// lp_assert(vectors_are_equal(wcopy, w.m_data));
#ifdef Z3DEBUG
// SASSERT(vectors_are_equal(wcopy, w.m_data));
#endif
}
@ -138,7 +153,7 @@ void row_eta_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
template <typename T, typename X>
void row_eta_matrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & p) {
// this = p * this * p(-1)
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// auto rev = p.get_reverse();
// auto deb = ((*this) * rev);
// deb = p * deb;
@ -150,11 +165,11 @@ void row_eta_matrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & p
columns.push_back(it.first);
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
// lp_assert(deb == *this);
#ifdef Z3DEBUG
// SASSERT(deb == *this);
#endif
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X>
T row_eta_matrix<T, X>::get_elem(unsigned row, unsigned col) const {
if (row == m_row){

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <memory>
#include "util/vector.h"
#include "util/lp/row_eta_matrix.hpp"
@ -10,7 +25,7 @@ namespace lp {
template void row_eta_matrix<double, double>::conjugate_by_permutation(permutation_matrix<double, double>&);
template void row_eta_matrix<mpq, numeric_pair<mpq> >::conjugate_by_permutation(permutation_matrix<mpq, numeric_pair<mpq> >&);
template void row_eta_matrix<mpq, mpq>::conjugate_by_permutation(permutation_matrix<mpq, mpq>&);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template mpq row_eta_matrix<mpq, mpq>::get_elem(unsigned int, unsigned int) const;
template mpq row_eta_matrix<mpq, numeric_pair<mpq> >::get_elem(unsigned int, unsigned int) const;
template double row_eta_matrix<double, double>::get_elem(unsigned int, unsigned int) const;

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -31,7 +46,7 @@ public:
m_scaling_maximum(scaling_maximum),
m_column_scale(column_scale),
m_settings(settings) {
lp_assert(m_column_scale.size() == 0);
SASSERT(m_column_scale.size() == 0);
m_column_scale.resize(m_A.column_count(), numeric_traits<T>::one());
}

29
src/util/lp/scaler.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <algorithm>
#include "util/lp/scaler.h"
#include "util/lp/numeric_pair.h"
@ -41,7 +56,7 @@ template <typename T, typename X> T scaler<T, X>::A_max() const {
template <typename T, typename X> T scaler<T, X>::get_A_ratio() const {
T min = A_min();
T max = A_max();
lp_assert(!m_settings.abs_val_is_smaller_than_zero_tolerance(min));
SASSERT(!m_settings.abs_val_is_smaller_than_zero_tolerance(min));
T ratio = max / min;
return ratio;
}
@ -51,7 +66,7 @@ template <typename T, typename X> T scaler<T, X>::get_max_ratio_on_rows() con
unsigned i = m_A.row_count();
while (i--) {
T den = m_A.get_min_abs_in_row(i);
lp_assert(!m_settings.abs_val_is_smaller_than_zero_tolerance(den));
SASSERT(!m_settings.abs_val_is_smaller_than_zero_tolerance(den));
T t = m_A.get_max_abs_in_row(i)/ den;
if (t > ret)
ret = t;
@ -78,7 +93,7 @@ template <typename T, typename X> void scaler<T, X>::scale_rows_with_geometri
while (i--) {
T max = m_A.get_max_abs_in_row(i);
T min = m_A.get_min_abs_in_row(i);
lp_assert(max > zero_of_type<T>() && min > zero_of_type<T>());
SASSERT(max > zero_of_type<T>() && min > zero_of_type<T>());
if (is_zero(max) || is_zero(min))
continue;
T gm = T(sqrt(numeric_traits<T>::get_double(max*min)));

23
src/util/lp/scaler_instances.cpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/scaler.hpp"
template bool lp::scaler<double, double>::scale();
template bool lp::scaler<lp::mpq, lp::mpq>::scale();

23
src/util/lp/signature_bound_evidence.h
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/lp_settings.h"
#include "util/lp/lar_constraints.h"

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -25,7 +40,7 @@ namespace lp {
// it is a square matrix
template <typename T, typename X>
class sparse_matrix
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
: public matrix<T, X>
#endif
{
@ -57,7 +72,7 @@ public:
vector<bool> m_processed;
unsigned get_n_of_active_elems() const { return m_n_of_active_elems; }
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// dense_matrix<T> m_dense;
#endif
/*
@ -146,7 +161,7 @@ public:
unsigned dimension() const {return static_cast<unsigned>(m_row_permutation.size());}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
unsigned row_count() const {return dimension();}
unsigned column_count() const {return dimension();}
#endif
@ -206,19 +221,19 @@ public:
void multiply_from_right(permutation_matrix<T, X>& p) {
// m_dense = m_dense * p;
m_column_permutation.multiply_by_permutation_from_right(p);
// lp_assert(*this == m_dense);
// SASSERT(*this == m_dense);
}
void multiply_from_left(permutation_matrix<T, X>& p) {
// m_dense = p * m_dense;
m_row_permutation.multiply_by_permutation_from_left(p);
// lp_assert(*this == m_dense);
// SASSERT(*this == m_dense);
}
void multiply_from_left_with_reverse(permutation_matrix<T, X>& p) {
// m_dense = p * m_dense;
m_row_permutation.multiply_by_permutation_reverse_from_left(p);
// lp_assert(*this == m_dense);
// SASSERT(*this == m_dense);
}
// adding delta columns at the end of the matrix
@ -231,13 +246,13 @@ public:
// dense_matrix<T, X> d(*this);
m_column_permutation.transpose_from_left(a, b);
// d.swap_columns(a, b);
// lp_assert(*this == d);
// SASSERT(*this == d);
}
void swap_rows(unsigned a, unsigned b) {
m_row_permutation.transpose_from_right(a, b);
// m_dense.swap_rows(a, b);
// lp_assert(*this == m_dense);
// SASSERT(*this == m_dense);
}
void divide_row_by_constant(unsigned i, const T & t, lp_settings & settings);
@ -286,7 +301,7 @@ public:
template <typename L>
void solve_U_y_indexed_only(indexed_vector<L> & y, const lp_settings&, vector<unsigned> & sorted_active_rows );
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
T get_elem(unsigned i, unsigned j) const { return get(i, j); }
unsigned get_number_of_rows() const { return dimension(); }
unsigned get_number_of_columns() const { return dimension(); }
@ -341,7 +356,7 @@ public:
bool shorten_active_matrix(unsigned row, eta_matrix<T, X> *eta_matrix);
unsigned pivot_score_without_shortened_counters(unsigned i, unsigned j, unsigned k);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
bool can_improve_score_for_row(unsigned row, unsigned score, T const & c_partial_pivoting, unsigned k);
bool really_best_pivot(unsigned i, unsigned j, T const & c_partial_pivoting, unsigned k);
void print_active_matrix(unsigned k, std::ostream & out);
@ -373,7 +388,7 @@ public:
}
bool fill_eta_matrix(unsigned j, eta_matrix<T, X> ** eta);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
bool is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings & settings) const;
bool is_upper_triangular_until(unsigned k) const;
@ -393,7 +408,7 @@ public:
void process_index_recursively_for_y_U(unsigned j, vector<unsigned> & sorted_rows);
void resize(unsigned new_dim) {
unsigned old_dim = dimension();
lp_assert(new_dim >= old_dim);
SASSERT(new_dim >= old_dim);
for (unsigned j = old_dim; j < new_dim; j++) {
m_rows.push_back(vector<indexed_value<T>>());
m_columns.push_back(col_header());
@ -407,7 +422,7 @@ public:
add_new_element(j, j, numeric_traits<T>::one());
}
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
vector<T> get_full_row(unsigned i) const;
#endif
unsigned pivot_queue_size() const { return m_pivot_queue.size(); }

133
src/util/lp/sparse_matrix.hpp
View file

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/vector.h"
#include "util/lp/sparse_matrix.h"
@ -82,12 +97,12 @@ void sparse_matrix<T, X>::set_with_no_adjusting(unsigned row, unsigned col, T va
template <typename T, typename X>
void sparse_matrix<T, X>::set(unsigned row, unsigned col, T val) { // should not be used in efficient code
lp_assert(row < dimension() && col < dimension());
SASSERT(row < dimension() && col < dimension());
// m_dense.set_elem(row, col, val);
row = adjust_row(row);
col = adjust_column(col);
set_with_no_adjusting(row, col, val);
// lp_assert(*this == m_dense);
// SASSERT(*this == m_dense);
}
template <typename T, typename X>
@ -243,7 +258,7 @@ void sparse_matrix<T, X>::scan_row_to_work_vector_and_remove_pivot_column(unsign
}
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X>
vector<T> sparse_matrix<T, X>::get_full_row(unsigned i) const {
vector<T> r;
@ -261,8 +276,8 @@ vector<T> sparse_matrix<T, X>::get_full_row(unsigned i) const {
// Returns false if the resulting row is all zeroes, and true otherwise
template <typename T, typename X>
bool sparse_matrix<T, X>::pivot_row_to_row(unsigned i, const T& alpha, unsigned i0, lp_settings & settings ) {
lp_assert(i < dimension() && i0 < dimension());
lp_assert(i != i0);
SASSERT(i < dimension() && i0 < dimension());
SASSERT(i != i0);
unsigned pivot_col = adjust_column(i);
i = adjust_row(i);
i0 = adjust_row(i0);
@ -327,7 +342,7 @@ bool sparse_matrix<T, X>::set_row_from_work_vector_and_clean_work_vector_not_adj
if (numeric_traits<T>::is_zero(work_vec[j])) {
continue;
}
lp_assert(!settings.abs_val_is_smaller_than_drop_tolerance(work_vec[j]));
SASSERT(!settings.abs_val_is_smaller_than_drop_tolerance(work_vec[j]));
add_new_element(i0, adjust_column(j), work_vec[j]);
work_vec[j] = numeric_traits<T>::zero();
}
@ -372,7 +387,7 @@ void sparse_matrix<T, X>::remove_zero_elements_and_set_data_on_existing_elements
T val = work_vec[rj];
if (settings.abs_val_is_smaller_than_drop_tolerance(val)) {
remove_element(row_vals, row_el_iv);
lp_assert(numeric_traits<T>::is_zero(val));
SASSERT(numeric_traits<T>::is_zero(val));
} else {
m_columns[j].m_values[row_el_iv.m_other].set_value(row_el_iv.m_value = val);
work_vec[rj] = numeric_traits<T>::zero();
@ -393,7 +408,7 @@ void sparse_matrix<T, X>::add_columns_at_the_end(unsigned delta) {
template <typename T, typename X>
void sparse_matrix<T, X>::delete_column(int i) {
lp_assert(i < dimension());
SASSERT(i < dimension());
for (auto cell = m_columns[i].m_head; cell != nullptr;) {
auto next_cell = cell->m_down;
kill_cell(cell);
@ -403,7 +418,7 @@ void sparse_matrix<T, X>::delete_column(int i) {
template <typename T, typename X>
void sparse_matrix<T, X>::divide_row_by_constant(unsigned i, const T & t, lp_settings & settings) {
lp_assert(!settings.abs_val_is_smaller_than_zero_tolerance(t));
SASSERT(!settings.abs_val_is_smaller_than_zero_tolerance(t));
i = adjust_row(i);
for (auto & iv : m_rows[i]) {
T &v = iv.m_value;
@ -420,7 +435,7 @@ void sparse_matrix<T, X>::divide_row_by_constant(unsigned i, const T & t, lp_set
// the matrix here has to be upper triangular
template <typename T, typename X>
void sparse_matrix<T, X>::solve_y_U(vector<T> & y) const { // works by rows
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// T * rs = clone_vector<T>(y, dimension());
#endif
unsigned end = dimension();
@ -436,11 +451,11 @@ void sparse_matrix<T, X>::solve_y_U(vector<T> & y) const { // works by rows
}
}
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// dense_matrix<T> deb(*this);
// T * clone_y = clone_vector<T>(y, dimension());
// deb.apply_from_right(clone_y);
// lp_assert(vectors_are_equal(rs, clone_y, dimension()));
// SASSERT(vectors_are_equal(rs, clone_y, dimension()));
// delete [] clone_y;
// delete [] rs;
#endif
@ -450,7 +465,7 @@ void sparse_matrix<T, X>::solve_y_U(vector<T> & y) const { // works by rows
// the matrix here has to be upper triangular
template <typename T, typename X>
void sparse_matrix<T, X>::solve_y_U_indexed(indexed_vector<T> & y, const lp_settings & settings) {
#if 0 && LEAN_DEBUG
#if 0 && Z3DEBUG
vector<T> ycopy(y.m_data);
if (numeric_traits<T>::precise() == false)
solve_y_U(ycopy);
@ -474,10 +489,10 @@ void sparse_matrix<T, X>::solve_y_U_indexed(indexed_vector<T> & y, const lp_sett
y.m_data[j] = zero_of_type<T>();
}
lp_assert(y.is_OK());
#if 0 && LEAN_DEBUG
SASSERT(y.is_OK());
#if 0 && Z3DEBUG
if (numeric_traits<T>::precise() == false)
lp_assert(vectors_are_equal(ycopy, y.m_data));
SASSERT(vectors_are_equal(ycopy, y.m_data));
#endif
}
@ -537,8 +552,8 @@ void sparse_matrix<T, X>::add_delta_to_solution(const vector<L>& del, vector<L>
template <typename T, typename X>
template <typename L>
void sparse_matrix<T, X>::add_delta_to_solution(const indexed_vector<L>& del, indexed_vector<L> & y) {
// lp_assert(del.is_OK());
// lp_assert(y.is_OK());
// SASSERT(del.is_OK());
// SASSERT(y.is_OK());
for (auto i : del.m_index) {
y.add_value_at_index(i, del[i]);
}
@ -546,11 +561,11 @@ void sparse_matrix<T, X>::add_delta_to_solution(const indexed_vector<L>& del, in
template <typename T, typename X>
template <typename L>
void sparse_matrix<T, X>::double_solve_U_y(indexed_vector<L>& y, const lp_settings & settings){
lp_assert(y.is_OK());
SASSERT(y.is_OK());
indexed_vector<L> y_orig(y); // copy y aside
vector<unsigned> active_rows;
solve_U_y_indexed_only(y, settings, active_rows);
lp_assert(y.is_OK());
SASSERT(y.is_OK());
find_error_in_solution_U_y_indexed(y_orig, y, active_rows);
// y_orig contains the error now
if (y_orig.m_index.size() * ratio_of_index_size_to_all_size<T>() < 32 * dimension()) {
@ -563,7 +578,7 @@ void sparse_matrix<T, X>::double_solve_U_y(indexed_vector<L>& y, const lp_settin
add_delta_to_solution(y_orig.m_data, y.m_data);
y.restore_index_and_clean_from_data();
}
lp_assert(y.is_OK());
SASSERT(y.is_OK());
}
template <typename T, typename X>
template <typename L>
@ -581,7 +596,7 @@ void sparse_matrix<T, X>::double_solve_U_y(vector<L>& y){
template <typename T, typename X>
template <typename L>
void sparse_matrix<T, X>::solve_U_y(vector<L> & y) { // it is a column wise version
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// T * rs = clone_vector<T>(y, dimension());
#endif
@ -595,16 +610,16 @@ void sparse_matrix<T, X>::solve_U_y(vector<L> & y) { // it is a column wise vers
}
}
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// dense_matrix<T> deb(*this);
// T * clone_y = clone_vector<T>(y, dimension());
// deb.apply_from_left(clone_y);
// lp_assert(vectors_are_equal(rs, clone_y, dimension()));
// SASSERT(vectors_are_equal(rs, clone_y, dimension()));
#endif
}
template <typename T, typename X>
void sparse_matrix<T, X>::process_index_recursively_for_y_U(unsigned j, vector<unsigned> & sorted_active_rows) {
lp_assert(m_processed[j] == false);
SASSERT(m_processed[j] == false);
m_processed[j]=true;
auto & row = m_rows[adjust_row(j)];
for (auto & c : row) {
@ -619,7 +634,7 @@ void sparse_matrix<T, X>::process_index_recursively_for_y_U(unsigned j, vector<u
template <typename T, typename X>
void sparse_matrix<T, X>::process_column_recursively(unsigned j, vector<unsigned> & sorted_active_rows) {
lp_assert(m_processed[j] == false);
SASSERT(m_processed[j] == false);
auto & mc = m_columns[adjust_column(j)].m_values;
for (auto & iv : mc) {
unsigned i = adjust_row_inverse(iv.m_index);
@ -684,12 +699,12 @@ void sparse_matrix<T, X>::solve_U_y_indexed_only(indexed_vector<L> & y, const lp
y[j] = zero_of_type<L>();
}
lp_assert(y.is_OK());
#ifdef LEAN_DEBUG
SASSERT(y.is_OK());
#ifdef Z3DEBUG
// dense_matrix<T,X> deb(this);
// vector<T> clone_y(y.m_data);
// deb.apply_from_left(clone_y);
// lp_assert(vectors_are_equal(rs, clone_y));
// SASSERT(vectors_are_equal(rs, clone_y));
#endif
}
@ -802,7 +817,7 @@ void sparse_matrix<T, X>::add_new_elements_of_w_and_clear_w(unsigned column_to_r
unsigned ai = adjust_row(i);
add_new_element(ai, column_to_replace, w_at_i);
auto & row_chunk = m_rows[ai];
lp_assert(row_chunk.size() > 0);
SASSERT(row_chunk.size() > 0);
if (abs(w_at_i) > abs(row_chunk[0].m_value))
put_max_index_to_0(row_chunk, static_cast<unsigned>(row_chunk.size()) - 1);
}
@ -833,7 +848,7 @@ unsigned sparse_matrix<T, X>::pivot_score(unsigned i, unsigned j) {
template <typename T, typename X>
void sparse_matrix<T, X>::enqueue_domain_into_pivot_queue() {
lp_assert(m_pivot_queue.size() == 0);
SASSERT(m_pivot_queue.size() == 0);
for (unsigned i = 0; i < dimension(); i++) {
auto & rh = m_rows[i];
unsigned rnz = static_cast<unsigned>(rh.size());
@ -919,7 +934,7 @@ void sparse_matrix<T, X>::update_active_pivots(unsigned row) {
for (const auto & iv : m_rows[arow]) {
col_header & ch = m_columns[iv.m_index];
int cols = static_cast<int>(ch.m_values.size()) - ch.m_shortened_markovitz - 1;
lp_assert(cols >= 0);
SASSERT(cols >= 0);
for (const auto &ivc : ch.m_values) {
unsigned i = ivc.m_index;
if (adjust_row_inverse(i) <= row) continue; // the i is not an active row
@ -945,7 +960,7 @@ bool sparse_matrix<T, X>::shorten_active_matrix(unsigned row, eta_matrix<T, X> *
for (auto & iv : row_values) {
const col_header& ch = m_columns[iv.m_index];
int cnz = static_cast<int>(ch.m_values.size()) - ch.m_shortened_markovitz - 1;
lp_assert(cnz >= 0);
SASSERT(cnz >= 0);
m_pivot_queue.enqueue(row, iv.m_index, rnz * cnz);
}
}
@ -961,25 +976,25 @@ unsigned sparse_matrix<T, X>::pivot_score_without_shortened_counters(unsigned i,
if (adjust_row_inverse(iv.m_index) < k)
cnz--;
}
lp_assert(cnz > 0);
SASSERT(cnz > 0);
return m_rows[i].m_values.size() * (cnz - 1);
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X>
bool sparse_matrix<T, X>::can_improve_score_for_row(unsigned row, unsigned score, T const & c_partial_pivoting, unsigned k) {
unsigned arow = adjust_row(row);
auto & row_vals = m_rows[arow].m_values;
auto & begin_iv = row_vals[0];
T row_max = abs(begin_iv.m_value);
lp_assert(adjust_column_inverse(begin_iv.m_index) >= k);
SASSERT(adjust_column_inverse(begin_iv.m_index) >= k);
if (pivot_score_without_shortened_counters(arow, begin_iv.m_index, k) < score) {
print_active_matrix(k);
return true;
}
for (unsigned jj = 1; jj < row_vals.size(); jj++) {
auto & iv = row_vals[jj];
lp_assert(adjust_column_inverse(iv.m_index) >= k);
lp_assert(abs(iv.m_value) <= row_max);
SASSERT(adjust_column_inverse(iv.m_index) >= k);
SASSERT(abs(iv.m_value) <= row_max);
if (c_partial_pivoting * abs(iv.m_value) < row_max) continue;
if (pivot_score_without_shortened_counters(arow, iv.m_index, k) < score) {
print_active_matrix(k);
@ -993,7 +1008,7 @@ template <typename T, typename X>
bool sparse_matrix<T, X>::really_best_pivot(unsigned i, unsigned j, T const & c_partial_pivoting, unsigned k) {
unsigned queue_pivot_score = pivot_score_without_shortened_counters(i, j, k);
for (unsigned ii = k; ii < dimension(); ii++) {
lp_assert(!can_improve_score_for_row(ii, queue_pivot_score, c_partial_pivoting, k));
SASSERT(!can_improve_score_for_row(ii, queue_pivot_score, c_partial_pivoting, k));
}
return true;
}
@ -1026,7 +1041,7 @@ template <typename T, typename X>
bool sparse_matrix<T, X>::pivot_queue_is_correct_for_row(unsigned i, unsigned k) {
unsigned arow = adjust_row(i);
for (auto & iv : m_rows[arow].m_values) {
lp_assert(pivot_score_without_shortened_counters(arow, iv.m_index, k + 1) ==
SASSERT(pivot_score_without_shortened_counters(arow, iv.m_index, k + 1) ==
m_pivot_queue.get_priority(arow, iv.m_index));
}
return true;
@ -1035,8 +1050,8 @@ bool sparse_matrix<T, X>::pivot_queue_is_correct_for_row(unsigned i, unsigned k)
template <typename T, typename X>
bool sparse_matrix<T, X>::pivot_queue_is_correct_after_pivoting(int k) {
for (unsigned i = k + 1; i < dimension(); i++ )
lp_assert(pivot_queue_is_correct_for_row(i, k));
lp_assert(m_pivot_queue.is_correct());
SASSERT(pivot_queue_is_correct_for_row(i, k));
SASSERT(m_pivot_queue.is_correct());
return true;
}
#endif
@ -1052,10 +1067,10 @@ bool sparse_matrix<T, X>::get_pivot_for_column(unsigned &i, unsigned &j, int c_p
if (j_inv < k) continue;
int _small = elem_is_too_small(i, j, c_partial_pivoting);
if (!_small) {
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
// if (!really_best_pivot(i, j, c_partial_pivoting, k)) {
// print_active_matrix(k);
// lp_assert(false);
// SASSERT(false);
// }
#endif
recover_pivot_queue(pivots_candidates_that_are_too_small);
@ -1088,7 +1103,7 @@ bool sparse_matrix<T, X>::shorten_columns_by_pivot_row(unsigned i, unsigned pivo
for (indexed_value<T> & iv : row_chunk) {
unsigned j = iv.m_index;
if (j == pivot_column) {
lp_assert(!col_is_active(j));
SASSERT(!col_is_active(j));
continue;
}
m_columns[j].shorten_markovich_by_one();
@ -1121,7 +1136,7 @@ bool sparse_matrix<T, X>::fill_eta_matrix(unsigned j, eta_matrix<T, X> ** eta) {
return true;
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
*eta = new eta_matrix<T, X>(j, dimension());
#else
*eta = new eta_matrix<T, X>(j);
@ -1146,16 +1161,16 @@ bool sparse_matrix<T, X>::fill_eta_matrix(unsigned j, eta_matrix<T, X> ** eta) {
(*eta)->divide_by_diagonal_element();
return true;
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template <typename T, typename X>
bool sparse_matrix<T, X>::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings & settings) const {
for (unsigned i = 0; i < dimension(); i++) {
vector<indexed_value<T>> const & row_chunk = get_row_values(i);
lp_assert(row_chunk.size());
SASSERT(row_chunk.size());
T const & max = abs(row_chunk[0].m_value);
unsigned ai = adjust_row_inverse(i);
for (auto & iv : row_chunk) {
lp_assert(abs(iv.m_value) <= max);
SASSERT(abs(iv.m_value) <= max);
unsigned aj = adjust_column_inverse(iv.m_index);
if (!(ai <= aj || numeric_traits<T>::is_zero(iv.m_value)))
return false;
@ -1193,18 +1208,18 @@ void sparse_matrix<T, X>::check_column_vs_rows(unsigned col) {
indexed_value<T> & row_iv = column_iv_other(column_iv);
if (row_iv.m_index != col) {
// std::cout << "m_other in row does not belong to column " << col << ", but to column " << row_iv.m_index << std::endl;
lp_assert(false);
SASSERT(false);
}
if (& row_iv_other(row_iv) != &column_iv) {
// std::cout << "row and col do not point to each other" << std::endl;
lp_assert(false);
SASSERT(false);
}
if (row_iv.m_value != column_iv.m_value) {
// std::cout << "the data from col " << col << " for row " << column_iv.m_index << " is different in the column " << std::endl;
// std::cout << "in the col it is " << column_iv.m_value << ", but in the row it is " << row_iv.m_value << std::endl;
lp_assert(false);
SASSERT(false);
}
}
}
@ -1217,18 +1232,18 @@ void sparse_matrix<T, X>::check_row_vs_columns(unsigned row) {
if (column_iv.m_index != row) {
// std::cout << "col_iv does not point to correct row " << row << " but to " << column_iv.m_index << std::endl;
lp_assert(false);
SASSERT(false);
}
if (& row_iv != & column_iv_other(column_iv)) {
// std::cout << "row and col do not point to each other" << std::endl;
lp_assert(false);
SASSERT(false);
}
if (row_iv.m_value != column_iv.m_value) {
// std::cout << "the data from col " << column_iv.m_index << " for row " << row << " is different in the column " << std::endl;
// std::cout << "in the col it is " << column_iv.m_value << ", but in the row it is " << row_iv.m_value << std::endl;
lp_assert(false);
SASSERT(false);
}
}
}

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <memory>
#include "util/vector.h"
#include "util/lp/lp_settings.h"
@ -67,7 +82,7 @@ template void sparse_matrix<mpq, numeric_pair<mpq>>::double_solve_U_y<mpq>(index
template void sparse_matrix<mpq, numeric_pair<mpq> >::double_solve_U_y<numeric_pair<mpq> >(indexed_vector<numeric_pair<mpq>>&, const lp_settings&);
template void lp::sparse_matrix<double, double>::solve_U_y_indexed_only<double>(lp::indexed_vector<double>&, const lp_settings&, vector<unsigned> &);
template void lp::sparse_matrix<lp::mpq, lp::mpq>::solve_U_y_indexed_only<lp::mpq>(lp::indexed_vector<lp::mpq>&, const lp_settings &, vector<unsigned> &);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
template bool sparse_matrix<double, double>::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const;
template bool sparse_matrix<mpq, mpq>::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const;
template bool sparse_matrix<mpq, numeric_pair<mpq> >::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const;

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -18,7 +33,7 @@ public:
void push_back(unsigned index, T val) {
m_data.push_back(std::make_pair(index, val));
}
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
T operator[] (unsigned i) const {
for (auto &t : m_data) {
if (t.first == i) return t.second;
@ -27,7 +42,7 @@ public:
}
#endif
void divide(T const & a) {
lp_assert(!lp_settings::is_eps_small_general(a, 1e-12));
SASSERT(!lp_settings::is_eps_small_general(a, 1e-12));
for (auto & t : m_data) { t.second /= a; }
}

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

@ -1,7 +1,22 @@
/*
Copyright (c) 2017 Microsoft Corporation
Author: Lev Nachmanson
*/
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
@ -30,11 +45,11 @@ class square_dense_submatrix : public tail_matrix<T, X> {
ref(unsigned i, square_dense_submatrix & s) :
m_i_offset((i - s.m_index_start) * s.m_dim), m_s(s){}
T & operator[] (unsigned j) {
lp_assert(j >= m_s.m_index_start);
SASSERT(j >= m_s.m_index_start);
return m_s.m_v[m_i_offset + m_s.adjust_column(j) - m_s.m_index_start];
}
const T & operator[] (unsigned j) const {
lp_assert(j >= m_s.m_index_start);
SASSERT(j >= m_s.m_index_start);
return m_s.m_v[m_i_offset + m_s.adjust_column(j) - m_s.m_index_start];
}
};
@ -58,8 +73,8 @@ public:
bool is_dense() const { return true; }
ref operator[] (unsigned i) {
lp_assert(i >= m_index_start);
lp_assert(i < m_parent->dimension());
SASSERT(i >= m_index_start);
SASSERT(i < m_parent->dimension());
return ref(i, *this);
}
@ -148,7 +163,7 @@ public:
}
}
}
lp_assert(wcopy.is_OK());
SASSERT(wcopy.is_OK());
apply_from_right(w.m_data);
w.m_index.clear();
if (numeric_traits<T>::precise()) {
@ -167,11 +182,11 @@ public:
}
}
#else
lp_assert(w.is_OK());
lp_assert(m_work_vector.is_OK());
SASSERT(w.is_OK());
SASSERT(m_work_vector.is_OK());
m_work_vector.resize(w.data_size());
m_work_vector.clear();
lp_assert(m_work_vector.is_OK());
SASSERT(m_work_vector.is_OK());
unsigned end = m_index_start + m_dim;
for (unsigned k : w.m_index) {
// find j such that k = adjust_row_inverse(j)
@ -188,7 +203,7 @@ public:
}
}
m_work_vector.clean_up();
lp_assert(m_work_vector.is_OK());
SASSERT(m_work_vector.is_OK());
w = m_work_vector;
#endif
}
@ -198,7 +213,7 @@ public:
void apply_from_right(vector<T> & w);
#ifdef LEAN_DEBUG
#ifdef Z3DEBUG
T get_elem (unsigned i, unsigned j) const;
unsigned row_count() const { return m_parent->row_count();}
unsigned column_count() const { return row_count();}

Some files were not shown because too many files have changed in this diff Show more