3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-23 09:05:31 +00:00
Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
This commit is contained in:
Nikolaj Bjorner 2018-10-14 15:16:22 -07:00
commit e9d615e309
204 changed files with 4620 additions and 2435 deletions

View file

@ -1,11 +1,12 @@
if (EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/version.h")
message(FATAL_ERROR "\"${CMAKE_CURRENT_SOURCE_DIR}/version.h\""
if (EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/z3_version.h")
message(FATAL_ERROR "\"${CMAKE_CURRENT_SOURCE_DIR}/z3_version.h\""
${z3_polluted_tree_msg}
)
endif()
set(Z3_FULL_VERSION "\"${Z3_FULL_VERSION_STR}\"")
configure_file(version.h.cmake.in ${CMAKE_CURRENT_BINARY_DIR}/version.h)
configure_file(z3_version.h.cmake.in ${CMAKE_CURRENT_BINARY_DIR}/z3_version.h)
z3_add_component(util
SOURCES

View file

@ -45,7 +45,7 @@ Revision History:
// For SSE2, it is best to use compiler intrinsics because this makes it completely
// clear to the compiler what instructions should be used. E.g., for sqrt(), the Windows compiler selects
// the x87 FPU, even when /arch:SSE2 is on.
// Luckily, these are kind of standardized, at least for Windows/Linux/OSX.
// Luckily, these are kind of standardized, at least for Windows/Linux/macOS.
#ifdef __clang__
#undef USE_INTRINSICS
#endif
@ -75,14 +75,14 @@ hwf_manager::hwf_manager() :
#endif
#endif
#else
// OSX/Linux: Nothing.
// macOS/Linux: Nothing.
#endif
// We only set the precision of the FPU here in the constructor. At the moment, there are no
// other parts of the code that could overwrite this, and Windows takes care of context switches.
// CMW: I'm not sure what happens on CPUs with hyper-threading (since the FPU is shared).
// I have yet to discover whether Linux and OSX save the FPU state when switching context.
// I have yet to discover whether Linux and macOS save the FPU state when switching context.
// As long as we stick to using the SSE2 FPU though, there shouldn't be any problems with respect
// to the precision (not sure about the rounding modes though).
}
@ -279,7 +279,7 @@ void hwf_manager::fma(mpf_rounding_mode rm, hwf const & x, hwf const & y, hwf co
#endif
#endif
#else
// Linux, OSX
// Linux, macOS
o.value = ::fma(x.value, y.value, z.value);
#endif
#endif
@ -331,7 +331,7 @@ void hwf_manager::round_to_integral(mpf_rounding_mode rm, hwf const & x, hwf & o
}
#endif
#else
// Linux, OSX.
// Linux, macOS.
o.value = nearbyint(x.value);
#endif
}
@ -623,7 +623,7 @@ void hwf_manager::set_rounding_mode(mpf_rounding_mode rm)
UNREACHABLE(); // Note: MPF_ROUND_NEAREST_TAWAY is not supported by the hardware!
}
#endif
#else // OSX/Linux
#else // macOS/Linux
switch (rm) {
case MPF_ROUND_NEAREST_TEVEN:
SETRM(FE_TONEAREST);

View file

@ -2,10 +2,11 @@ z3_add_component(lp
SOURCES
binary_heap_priority_queue.cpp
binary_heap_upair_queue.cpp
bound_propagator.cpp
lp_bound_propagator.cpp
core_solver_pretty_printer.cpp
dense_matrix.cpp
eta_matrix.cpp
gomory.cpp
indexed_vector.cpp
int_solver.cpp
lar_solver.cpp

View file

@ -69,16 +69,6 @@ public:
m_column_index(static_cast<unsigned>(-1))
{}
column_info(unsigned column_index) :
m_lower_bound_is_set(false),
m_lower_bound_is_strict(false),
m_upper_bound_is_set (false),
m_upper_bound_is_strict (false),
m_is_fixed(false),
m_cost(numeric_traits<T>::zero()),
m_column_index(column_index) {
}
column_info(const column_info & ci) {
m_name = ci.m_name;
m_lower_bound_is_set = ci.m_lower_bound_is_set;

View file

@ -33,29 +33,6 @@ public:
print_linear_combination_of_column_indices(coeff, out);
}
template <typename T>
void print_linear_combination_of_column_indices_only(const vector<std::pair<T, unsigned>> & coeffs, std::ostream & out) const {
bool first = true;
for (const auto & it : coeffs) {
auto val = it.first;
if (first) {
first = false;
} else {
if (numeric_traits<T>::is_pos(val)) {
out << " + ";
} else {
out << " - ";
val = -val;
}
}
if (val == -numeric_traits<T>::one())
out << " - ";
else if (val != numeric_traits<T>::one())
out << T_to_string(val);
out << "v" << it.second;
}
}
template <typename T>

341
src/util/lp/gomory.cpp Normal file
View file

@ -0,0 +1,341 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Nikolaj Bjorner (nbjorner)
Lev Nachmanson (levnach)
Revision History:
--*/
#include "util/lp/gomory.h"
#include "util/lp/int_solver.h"
#include "util/lp/lar_solver.h"
#include "util/lp/lp_utils.h"
namespace lp {
class gomory::imp {
lar_term & m_t; // the term to return in the cut
mpq & m_k; // the right side of the cut
explanation& m_ex; // the conflict explanation
unsigned m_inf_col; // a basis column which has to be an integer but has a non integral value
const row_strip<mpq>& m_row;
const int_solver& m_int_solver;
mpq m_lcm_den;
mpq m_f;
mpq m_one_minus_f;
mpq m_fj;
mpq m_one_minus_fj;
const impq & get_value(unsigned j) const { return m_int_solver.get_value(j); }
bool is_real(unsigned j) const { return m_int_solver.is_real(j); }
bool at_lower(unsigned j) const { return m_int_solver.at_lower(j); }
bool at_upper(unsigned j) const { return m_int_solver.at_upper(j); }
const impq & lower_bound(unsigned j) const { return m_int_solver.lower_bound(j); }
const impq & upper_bound(unsigned j) const { return m_int_solver.upper_bound(j); }
constraint_index column_lower_bound_constraint(unsigned j) const { return m_int_solver.column_lower_bound_constraint(j); }
constraint_index column_upper_bound_constraint(unsigned j) const { return m_int_solver.column_upper_bound_constraint(j); }
bool column_is_fixed(unsigned j) const { return m_int_solver.m_lar_solver->column_is_fixed(j); }
void int_case_in_gomory_cut(unsigned j) {
lp_assert(is_int(j) && m_fj.is_pos());
TRACE("gomory_cut_detail",
tout << " k = " << m_k;
tout << ", fj: " << m_fj << ", ";
tout << (at_lower(j)?"at_lower":"at_upper")<< std::endl;
);
mpq new_a;
if (at_lower(j)) {
new_a = m_fj <= m_one_minus_f ? m_fj / m_one_minus_f : ((1 - m_fj) / m_f);
lp_assert(new_a.is_pos());
m_k.addmul(new_a, lower_bound(j).x);
m_ex.push_justification(column_lower_bound_constraint(j));
}
else {
lp_assert(at_upper(j));
// the upper terms are inverted: therefore we have the minus
new_a = - (m_fj <= m_f ? m_fj / m_f : ((1 - m_fj) / m_one_minus_f));
lp_assert(new_a.is_neg());
m_k.addmul(new_a, upper_bound(j).x);
m_ex.push_justification(column_upper_bound_constraint(j));
}
m_t.add_monomial(new_a, j);
m_lcm_den = lcm(m_lcm_den, denominator(new_a));
TRACE("gomory_cut_detail", tout << "v" << j << " new_a = " << new_a << ", k = " << m_k << ", m_lcm_den = " << m_lcm_den << "\n";);
}
void real_case_in_gomory_cut(const mpq & a, unsigned j) {
TRACE("gomory_cut_detail_real", tout << "real\n";);
mpq new_a;
if (at_lower(j)) {
if (a.is_pos()) {
new_a = a / m_one_minus_f;
}
else {
new_a = - a / m_f;
}
m_k.addmul(new_a, lower_bound(j).x); // is it a faster operation than
// k += lower_bound(j).x * new_a;
m_ex.push_justification(column_lower_bound_constraint(j));
}
else {
lp_assert(at_upper(j));
if (a.is_pos()) {
new_a = - a / m_f;
}
else {
new_a = a / m_one_minus_f;
}
m_k.addmul(new_a, upper_bound(j).x); // k += upper_bound(j).x * new_a;
m_ex.push_justification(column_upper_bound_constraint(j));
}
TRACE("gomory_cut_detail_real", tout << a << "*v" << j << " k: " << m_k << "\n";);
m_t.add_monomial(new_a, j);
}
lia_move report_conflict_from_gomory_cut() {
lp_assert(m_k.is_pos());
// conflict 0 >= k where k is positive
m_k.neg(); // returning 0 <= -k
return lia_move::conflict;
}
void adjust_term_and_k_for_some_ints_case_gomory() {
lp_assert(!m_t.is_empty());
// k = 1 + sum of m_t at bounds
auto pol = m_t.coeffs_as_vector();
m_t.clear();
if (pol.size() == 1) {
TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;);
unsigned v = pol[0].second;
lp_assert(is_int(v));
const mpq& a = pol[0].first;
m_k /= a;
if (a.is_pos()) { // we have av >= k
if (!m_k.is_int())
m_k = ceil(m_k);
// switch size
m_t.add_monomial(- mpq(1), v);
m_k.neg();
} else {
if (!m_k.is_int())
m_k = floor(m_k);
m_t.add_monomial(mpq(1), v);
}
} else {
m_lcm_den = lcm(m_lcm_den, denominator(m_k));
lp_assert(m_lcm_den.is_pos());
TRACE("gomory_cut_detail", tout << "pol.size() > 1 den: " << m_lcm_den << std::endl;);
if (!m_lcm_den.is_one()) {
// normalize coefficients of integer parameters to be integers.
for (auto & pi: pol) {
pi.first *= m_lcm_den;
SASSERT(!is_int(pi.second) || pi.first.is_int());
}
m_k *= m_lcm_den;
}
// negate everything to return -pol <= -m_k
for (const auto & pi: pol) {
TRACE("gomory_cut", tout << pi.first << "* " << "v" << pi.second << "\n";);
m_t.add_monomial(-pi.first, pi.second);
}
m_k.neg();
}
TRACE("gomory_cut_detail", tout << "k = " << m_k << std::endl;);
lp_assert(m_k.is_int());
}
std::string var_name(unsigned j) const {
return std::string("x") + std::to_string(j);
}
std::ostream& dump_coeff_val(std::ostream & out, const mpq & a) const {
if (a.is_int()) {
out << a;
}
else if ( a >= zero_of_type<mpq>())
out << "(/ " << numerator(a) << " " << denominator(a) << ")";
else {
out << "(- ( / " << numerator(-a) << " " << denominator(-a) << "))";
}
return out;
}
template <typename T>
void dump_coeff(std::ostream & out, const T& c) const {
out << "( * ";
dump_coeff_val(out, c.coeff());
out << " " << var_name(c.var()) << ")";
}
std::ostream& dump_row_coefficients(std::ostream & out) const {
mpq lc(1);
for (const auto& p : m_row) {
lc = lcm(lc, denominator(p.coeff()));
}
for (const auto& p : m_row) {
dump_coeff_val(out << " (* ", p.coeff()*lc) << " " << var_name(p.var()) << ")";
}
return out;
}
void dump_the_row(std::ostream& out) const {
out << "; the row, excluding fixed vars\n";
out << "(assert ( = ( +";
dump_row_coefficients(out) << ") 0))\n";
}
void dump_declaration(std::ostream& out, unsigned v) const {
out << "(declare-const " << var_name(v) << (is_int(v) ? " Int" : " Real") << ")\n";
}
void dump_declarations(std::ostream& out) const {
// for a column j the var name is vj
for (const auto & p : m_row) {
dump_declaration(out, p.var());
}
for (const auto& p : m_t) {
unsigned v = p.var();
if (m_int_solver.m_lar_solver->is_term(v)) {
dump_declaration(out, v);
}
}
}
void dump_lower_bound_expl(std::ostream & out, unsigned j) const {
out << "(assert (>= " << var_name(j) << " " << lower_bound(j).x << "))\n";
}
void dump_upper_bound_expl(std::ostream & out, unsigned j) const {
out << "(assert (<= " << var_name(j) << " " << upper_bound(j).x << "))\n";
}
void dump_explanations(std::ostream& out) const {
for (const auto & p : m_row) {
unsigned j = p.var();
if (j == m_inf_col || (!is_real(j) && p.coeff().is_int())) {
continue;
}
else if (at_lower(j)) {
dump_lower_bound_expl(out, j);
} else {
lp_assert(at_upper(j));
dump_upper_bound_expl(out, j);
}
}
}
std::ostream& dump_term_coefficients(std::ostream & out) const {
for (const auto& p : m_t) {
dump_coeff(out, p);
}
return out;
}
std::ostream& dump_term_sum(std::ostream & out) const {
return dump_term_coefficients(out << "(+ ") << ")";
}
std::ostream& dump_term_le_k(std::ostream & out) const {
return dump_term_sum(out << "(<= ") << " " << m_k << ")";
}
void dump_the_cut_assert(std::ostream & out) const {
dump_term_le_k(out << "(assert (not ") << "))\n";
}
void dump_cut_and_constraints_as_smt_lemma(std::ostream& out) const {
dump_declarations(out);
dump_the_row(out);
dump_explanations(out);
dump_the_cut_assert(out);
out << "(check-sat)\n";
}
public:
lia_move create_cut() {
TRACE("gomory_cut",
tout << "applying cut at:\n"; print_linear_combination_of_column_indices_only(m_row, tout); tout << std::endl;
for (auto & p : m_row) {
m_int_solver.m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout);
}
tout << "inf_col = " << m_inf_col << std::endl;
);
// gomory will be t <= k and the current solution has a property t > k
m_k = 1;
m_t.clear();
mpq m_lcm_den(1);
bool some_int_columns = false;
mpq m_f = fractional_part(get_value(m_inf_col));
TRACE("gomory_cut_detail", tout << "m_f: " << m_f << ", ";
tout << "1 - m_f: " << 1 - m_f << ", get_value(m_inf_col).x - m_f = " << get_value(m_inf_col).x - m_f;);
lp_assert(m_f.is_pos() && (get_value(m_inf_col).x - m_f).is_int());
mpq one_min_m_f = 1 - m_f;
for (const auto & p : m_row) {
unsigned j = p.var();
if (j == m_inf_col) {
lp_assert(p.coeff() == one_of_type<mpq>());
TRACE("gomory_cut_detail", tout << "seeing basic var";);
continue;
}
// use -p.coeff() to make the format compatible with the format used in: Integrating Simplex with DPLL(T)
if (is_real(j)) {
real_case_in_gomory_cut(- p.coeff(), j);
} else {
if (p.coeff().is_int()) {
// m_fj will be zero and no monomial will be added
continue;
}
some_int_columns = true;
m_fj = fractional_part(-p.coeff());
m_one_minus_fj = 1 - m_fj;
int_case_in_gomory_cut(j);
}
}
if (m_t.is_empty())
return report_conflict_from_gomory_cut();
if (some_int_columns)
adjust_term_and_k_for_some_ints_case_gomory();
lp_assert(m_int_solver.current_solution_is_inf_on_cut());
TRACE("gomory_cut_detail", dump_cut_and_constraints_as_smt_lemma(tout););
m_int_solver.m_lar_solver->subs_term_columns(m_t);
TRACE("gomory_cut", print_linear_combination_of_column_indices_only(m_t, tout << "gomory cut:"); tout << " <= " << m_k << std::endl;);
return lia_move::cut;
}
imp(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip<mpq>& row, const int_solver& int_slv ) :
m_t(t),
m_k(k),
m_ex(ex),
m_inf_col(basic_inf_int_j),
m_row(row),
m_int_solver(int_slv),
m_lcm_den(1),
m_f(fractional_part(get_value(basic_inf_int_j).x)),
m_one_minus_f(1 - m_f) {}
};
lia_move gomory::create_cut() {
return m_imp->create_cut();
}
gomory::gomory(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip<mpq>& row, const int_solver& s) {
m_imp = alloc(imp, t, k, ex, basic_inf_int_j, row, s);
}
gomory::~gomory() { dealloc(m_imp); }
}

36
src/util/lp/gomory.h Normal file
View file

@ -0,0 +1,36 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Nikolaj Bjorner (nbjorner)
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/lp/lar_term.h"
#include "util/lp/lia_move.h"
#include "util/lp/explanation.h"
#include "util/lp/static_matrix.h"
namespace lp {
class int_solver;
class gomory {
class imp;
imp *m_imp;
public :
gomory(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip<mpq>& row, const int_solver& s);
lia_move create_cut();
~gomory();
};
}

View file

@ -8,6 +8,7 @@
#include "util/lp/lp_utils.h"
#include <utility>
#include "util/lp/monomial.h"
#include "util/lp/gomory.h"
namespace lp {
@ -101,10 +102,7 @@ bool int_solver::is_gomory_cut_target(const row_strip<mpq>& row) {
unsigned j;
for (const auto & p : row) {
j = p.var();
if (is_base(j)) continue;
if (!at_bound(j))
return false;
if (!is_zero(get_value(j).y)) {
if (!is_base(j) && (!at_bound(j) || !is_zero(get_value(j).y))) {
TRACE("gomory_cut", tout << "row is not gomory cut target:\n";
display_column(tout, j);
tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";);
@ -115,36 +113,6 @@ bool int_solver::is_gomory_cut_target(const row_strip<mpq>& row) {
}
void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0) {
TRACE("gomory_cut_detail_real", tout << "real\n";);
mpq new_a;
if (at_low(x_j)) {
if (a.is_pos()) {
new_a = a / one_minus_f_0;
}
else {
new_a = a / f_0;
new_a.neg();
}
m_k->addmul(new_a, lower_bound(x_j).x); // is it a faster operation than
// k += lower_bound(x_j).x * new_a;
m_ex->push_justification(column_lower_bound_constraint(x_j), new_a);
}
else {
lp_assert(at_upper(x_j));
if (a.is_pos()) {
new_a = a / f_0;
new_a.neg(); // the upper terms are inverted.
}
else {
new_a = a / one_minus_f_0;
}
m_k->addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a;
m_ex->push_justification(column_upper_bound_constraint(x_j), new_a);
}
TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << *m_k << "\n";);
m_t->add_monomial(new_a, x_j);
}
constraint_index int_solver::column_upper_bound_constraint(unsigned j) const {
return m_lar_solver->get_column_upper_bound_witness(j);
@ -155,221 +123,31 @@ constraint_index int_solver::column_lower_bound_constraint(unsigned j) const {
}
void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j,
mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) {
lp_assert(is_int(x_j));
lp_assert(!a.is_int());
mpq f_j = fractional_part(a);
TRACE("gomory_cut_detail",
tout << a << " x_j" << x_j << " k = " << *m_k << "\n";
tout << "f_j: " << f_j << "\n";
tout << "f_0: " << f_0 << "\n";
tout << "1 - f_0: " << 1 - f_0 << "\n";
tout << "at_low(" << x_j << ") = " << at_low(x_j) << std::endl;
);
lp_assert (!f_j.is_zero());
mpq new_a;
if (at_low(x_j)) {
if (f_j <= one_minus_f_0) {
new_a = f_j / one_minus_f_0;
}
else {
new_a = (1 - f_j) / f_0;
}
m_k->addmul(new_a, lower_bound(x_j).x);
m_ex->push_justification(column_lower_bound_constraint(x_j), new_a);
}
else {
lp_assert(at_upper(x_j));
if (f_j <= f_0) {
new_a = f_j / f_0;
}
else {
new_a = (mpq(1) - f_j) / one_minus_f_0;
}
new_a.neg(); // the upper terms are inverted
m_k->addmul(new_a, upper_bound(x_j).x);
m_ex->push_justification(column_upper_bound_constraint(x_j), new_a);
}
TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << *m_k << "\n";);
m_t->add_monomial(new_a, x_j);
lcm_den = lcm(lcm_den, denominator(new_a));
}
lia_move int_solver::report_conflict_from_gomory_cut() {
TRACE("empty_pol",);
lp_assert(m_k->is_pos());
// conflict 0 >= k where k is positive
m_k->neg(); // returning 0 <= -k
return lia_move::conflict;
}
void int_solver::gomory_cut_adjust_t_and_k(vector<std::pair<mpq, unsigned>> & pol,
lar_term & t,
mpq &k,
bool some_ints,
mpq & lcm_den) {
if (!some_ints)
return;
t.clear();
if (pol.size() == 1) {
unsigned v = pol[0].second;
lp_assert(is_int(v));
bool k_is_int = k.is_int();
const mpq& a = pol[0].first;
k /= a;
if (a.is_pos()) { // we have av >= k
if (!k_is_int)
k = ceil(k);
// switch size
t.add_monomial(- mpq(1), v);
k.neg();
} else {
if (!k_is_int)
k = floor(k);
t.add_monomial(mpq(1), v);
}
} else if (some_ints) {
lcm_den = lcm(lcm_den, denominator(k));
lp_assert(lcm_den.is_pos());
if (!lcm_den.is_one()) {
// normalize coefficients of integer parameters to be integers.
for (auto & pi: pol) {
pi.first *= lcm_den;
SASSERT(!is_int(pi.second) || pi.first.is_int());
}
k *= lcm_den;
}
// negate everything to return -pol <= -k
for (const auto & pi: pol)
t.add_monomial(-pi.first, pi.second);
k.neg();
}
}
bool int_solver::current_solution_is_inf_on_cut() const {
const auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x;
impq v = m_t->apply(x);
mpq sign = *m_upper ? one_of_type<mpq>() : -one_of_type<mpq>();
CTRACE("current_solution_is_inf_on_cut", v * sign <= (*m_k) * sign,
tout << "m_upper = " << *m_upper << std::endl;
tout << "v = " << v << ", k = " << (*m_k) << std::endl;
impq v = m_t.apply(x);
mpq sign = m_upper ? one_of_type<mpq>() : -one_of_type<mpq>();
CTRACE("current_solution_is_inf_on_cut", v * sign <= m_k * sign,
tout << "m_upper = " << m_upper << std::endl;
tout << "v = " << v << ", k = " << m_k << std::endl;
);
return v * sign > (*m_k) * sign;
return v * sign > m_k * sign;
}
void int_solver::adjust_term_and_k_for_some_ints_case_gomory(mpq &lcm_den) {
lp_assert(!m_t->is_empty());
auto pol = m_t->coeffs_as_vector();
m_t->clear();
if (pol.size() == 1) {
TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;);
unsigned v = pol[0].second;
lp_assert(is_int(v));
const mpq& a = pol[0].first;
(*m_k) /= a;
if (a.is_pos()) { // we have av >= k
if (!(*m_k).is_int())
(*m_k) = ceil((*m_k));
// switch size
m_t->add_monomial(- mpq(1), v);
(*m_k).neg();
} else {
if (!(*m_k).is_int())
(*m_k) = floor((*m_k));
m_t->add_monomial(mpq(1), v);
}
} else {
TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;);
lcm_den = lcm(lcm_den, denominator((*m_k)));
lp_assert(lcm_den.is_pos());
if (!lcm_den.is_one()) {
// normalize coefficients of integer parameters to be integers.
for (auto & pi: pol) {
pi.first *= lcm_den;
SASSERT(!is_int(pi.second) || pi.first.is_int());
}
(*m_k) *= lcm_den;
}
// negate everything to return -pol <= -(*m_k)
for (const auto & pi: pol)
m_t->add_monomial(-pi.first, pi.second);
(*m_k).neg();
}
TRACE("gomory_cut_detail", tout << "k = " << (*m_k) << std::endl;);
lp_assert((*m_k).is_int());
}
lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip<mpq> & row) {
lp_assert(column_is_int_inf(inf_col));
TRACE("gomory_cut",
tout << "applying cut at:\n"; m_lar_solver->print_row(row, tout); tout << std::endl;
for (auto & p : row) {
m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout);
}
tout << "inf_col = " << inf_col << std::endl;
);
// gomory will be t <= k and the current solution has a property t > k
*m_k = 1;
mpq lcm_den(1);
unsigned x_j;
mpq a;
bool some_int_columns = false;
mpq f_0 = int_solver::fractional_part(get_value(inf_col));
mpq one_min_f_0 = 1 - f_0;
for (const auto & p : row) {
x_j = p.var();
if (x_j == inf_col)
continue;
// make the format compatible with the format used in: Integrating Simplex with DPLL(T)
a = p.coeff();
a.neg();
if (is_real(x_j))
real_case_in_gomory_cut(a, x_j, f_0, one_min_f_0);
else if (!a.is_int()) { // f_j will be zero and no monomial will be added
some_int_columns = true;
int_case_in_gomory_cut(a, x_j, lcm_den, f_0, one_min_f_0);
}
}
if (m_t->is_empty())
return report_conflict_from_gomory_cut();
if (some_int_columns)
adjust_term_and_k_for_some_ints_case_gomory(lcm_den);
lp_assert(current_solution_is_inf_on_cut());
m_lar_solver->subs_term_columns(*m_t);
TRACE("gomory_cut", tout<<"precut:"; m_lar_solver->print_term(*m_t, tout); tout << " <= " << *m_k << std::endl;);
return lia_move::cut;
}
int int_solver::find_free_var_in_gomory_row(const row_strip<mpq>& row) {
unsigned j;
for (const auto & p : row) {
j = p.var();
if (!is_base(j) && is_free(j))
return static_cast<int>(j);
}
return -1;
gomory gc(m_t, m_k, m_ex, inf_col, row, *this);
return gc.create_cut();
}
lia_move int_solver::proceed_with_gomory_cut(unsigned j) {
const row_strip<mpq>& row = m_lar_solver->get_row(row_of_basic_column(j));
if (-1 != find_free_var_in_gomory_row(row))
return lia_move::undef;
if (!is_gomory_cut_target(row))
return lia_move::undef;
return create_branch_on_column(j);
*m_upper = true;
m_upper = true;
return mk_gomory_cut(j, row);
}
@ -394,17 +172,19 @@ typedef monomial mono;
// this will allow to enable and disable tracking of the pivot rows
struct pivoted_rows_tracking_control {
lar_solver * m_lar_solver;
bool m_track_pivoted_rows;
pivoted_rows_tracking_control(lar_solver* ls) :
struct check_return_helper {
lar_solver * m_lar_solver;
const lia_move & m_r;
bool m_track_pivoted_rows;
check_return_helper(lar_solver* ls, const lia_move& r) :
m_lar_solver(ls),
m_r(r),
m_track_pivoted_rows(ls->get_track_pivoted_rows())
{
TRACE("pivoted_rows", tout << "pivoted rows = " << ls->m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows->size() << std::endl;);
m_lar_solver->set_track_pivoted_rows(false);
}
~pivoted_rows_tracking_control() {
~check_return_helper() {
TRACE("pivoted_rows", tout << "pivoted rows = " << m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows->size() << std::endl;);
m_lar_solver->set_track_pivoted_rows(m_track_pivoted_rows);
}
@ -589,21 +369,21 @@ lia_move int_solver::make_hnf_cut() {
#else
vector<mpq> x0;
#endif
lia_move r = m_hnf_cutter.create_cut(*m_t, *m_k, *m_ex, *m_upper, x0);
lia_move r = m_hnf_cutter.create_cut(m_t, m_k, m_ex, m_upper, x0);
if (r == lia_move::cut) {
TRACE("hnf_cut",
m_lar_solver->print_term(*m_t, tout << "cut:");
tout << " <= " << *m_k << std::endl;
m_lar_solver->print_term(m_t, tout << "cut:");
tout << " <= " << m_k << std::endl;
for (unsigned i : m_hnf_cutter.constraints_for_explanation()) {
m_lar_solver->print_constraint(i, tout);
}
);
lp_assert(current_solution_is_inf_on_cut());
settings().st().m_hnf_cuts++;
m_ex->clear();
m_ex.clear();
for (unsigned i : m_hnf_cutter.constraints_for_explanation()) {
m_ex->push_justification(i);
m_ex.push_justification(i);
}
}
return r;
@ -619,14 +399,17 @@ lia_move int_solver::hnf_cut() {
return lia_move::undef;
}
lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) {
lia_move int_solver::check() {
if (!has_inf_int()) return lia_move::sat;
m_t = &t; m_k = &k; m_ex = &ex; m_upper = &upper;
m_t.clear();
m_k.reset();
m_ex.clear();
m_upper = false;
lia_move r = run_gcd_test();
if (r != lia_move::undef) return r;
pivoted_rows_tracking_control pc(m_lar_solver);
check_return_helper pc(m_lar_solver, r);
if(settings().m_int_pivot_fixed_vars_from_basis)
m_lar_solver->pivot_fixed_vars_from_basis();
@ -862,8 +645,8 @@ bool int_solver::gcd_test_for_row(static_matrix<mpq, numeric_pair<mpq>> & A, uns
void int_solver::add_to_explanation_from_fixed_or_boxed_column(unsigned j) {
constraint_index lc, uc;
m_lar_solver->get_bound_constraint_witnesses_for_column(j, lc, uc);
m_ex->m_explanation.push_back(std::make_pair(mpq(1), lc));
m_ex->m_explanation.push_back(std::make_pair(mpq(1), uc));
m_ex.m_explanation.push_back(std::make_pair(mpq(1), lc));
m_ex.m_explanation.push_back(std::make_pair(mpq(1), uc));
}
void int_solver::fill_explanation_from_fixed_columns(const row_strip<mpq> & row) {
for (const auto & c : row) {
@ -1126,7 +909,7 @@ bool int_solver::at_bound(unsigned j) const {
}
}
bool int_solver::at_low(unsigned j) const {
bool int_solver::at_lower(unsigned j) const {
auto & mpq_solver = m_lar_solver->m_mpq_lar_core_solver.m_r_solver;
switch (mpq_solver.m_column_types[j] ) {
case column_type::fixed:
@ -1258,20 +1041,20 @@ const impq& int_solver::lower_bound(unsigned j) const {
lia_move int_solver::create_branch_on_column(int j) {
TRACE("check_main_int", tout << "branching" << std::endl;);
lp_assert(m_t->is_empty());
lp_assert(m_t.is_empty());
lp_assert(j != -1);
m_t->add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j));
m_t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j));
if (is_free(j)) {
*m_upper = true;
*m_k = mpq(0);
m_upper = true;
m_k = mpq(0);
} else {
*m_upper = left_branch_is_more_narrow_than_right(j);
*m_k = *m_upper? floor(get_value(j)) : ceil(get_value(j));
m_upper = left_branch_is_more_narrow_than_right(j);
m_k = m_upper? floor(get_value(j)) : ceil(get_value(j));
}
TRACE("arith_int", tout << "branching v" << j << " = " << get_value(j) << "\n";
display_column(tout, j);
tout << "k = " << *m_k << std::endl;
tout << "k = " << m_k << std::endl;
);
return lia_move::branch;

View file

@ -39,20 +39,31 @@ public:
// fields
lar_solver *m_lar_solver;
unsigned m_number_of_calls;
lar_term *m_t; // the term to return in the cut
mpq *m_k; // the right side of the cut
explanation *m_ex; // the conflict explanation
bool *m_upper; // we have a cut m_t*x <= k if m_upper is true nad m_t*x >= k otherwise
lar_term m_t; // the term to return in the cut
mpq m_k; // the right side of the cut
explanation m_ex; // the conflict explanation
bool m_upper; // we have a cut m_t*x <= k if m_upper is true nad m_t*x >= k otherwise
hnf_cutter m_hnf_cutter;
// methods
int_solver(lar_solver* lp);
// main function to check that the 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 & upper);
lia_move check();
lar_term const& get_term() const { return m_t; }
mpq const& get_offset() const { return m_k; }
explanation const& get_explanation() const { return m_ex; }
bool is_upper() const { return m_upper; }
bool move_non_basic_column_to_bounds(unsigned j);
lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex);
bool is_base(unsigned j) const;
bool is_real(unsigned j) const;
const impq & lower_bound(unsigned j) const;
const impq & upper_bound(unsigned j) const;
bool is_int(unsigned j) const;
const impq & get_value(unsigned j) const;
bool at_lower(unsigned j) const;
bool at_upper(unsigned j) const;
private:
@ -79,10 +90,7 @@ private:
void add_to_explanation_from_fixed_or_boxed_column(unsigned j);
lia_move patch_nbasic_columns();
bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m);
const impq & lower_bound(unsigned j) const;
const impq & upper_bound(unsigned j) const;
bool is_int(unsigned j) const;
bool is_real(unsigned j) const;
private:
bool is_boxed(unsigned j) const;
bool is_fixed(unsigned j) const;
bool is_free(unsigned j) const;
@ -91,7 +99,6 @@ private:
void set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val);
bool non_basic_columns_are_at_bounds() const;
bool is_feasible() const;
const impq & get_value(unsigned j) const;
bool column_is_int_inf(unsigned j) const;
void trace_inf_rows() const;
lia_move branch_or_sat();
@ -104,39 +111,22 @@ private:
bool move_non_basic_columns_to_bounds();
void branch_infeasible_int_var(unsigned);
lia_move mk_gomory_cut(unsigned inf_col, const row_strip<mpq>& row);
lia_move report_conflict_from_gomory_cut();
void adjust_term_and_k_for_some_ints_case_gomory(mpq& lcm_den);
lia_move proceed_with_gomory_cut(unsigned j);
int find_free_var_in_gomory_row(const row_strip<mpq>& );
bool is_gomory_cut_target(const row_strip<mpq>&);
bool at_bound(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);
}
public:
void display_column(std::ostream & out, unsigned j) const;
inline static
mpq fractional_part(const impq & n) {
lp_assert(is_rational(n));
return n.x - floor(n.x);
}
private:
void real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0);
void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0);
constraint_index column_upper_bound_constraint(unsigned j) const;
constraint_index column_lower_bound_constraint(unsigned j) const;
void display_row_info(std::ostream & out, unsigned row_index) const;
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;
public:
bool shift_var(unsigned j, unsigned range);
private:
void display_row_info(std::ostream & out, unsigned row_index) const;
unsigned random();
bool has_inf_int() const;
lia_move create_branch_on_column(int j);
@ -161,5 +151,5 @@ public:
bool hnf_has_var_with_non_integral_value() const;
bool hnf_cutter_is_full() const;
void patch_nbasic_column(unsigned j, bool patch_only_int_vals);
};
};
}

View file

@ -75,7 +75,7 @@ struct lar_term_constraint: public lar_base_constraint {
}
unsigned size() const override { return m_term->size();}
lar_term_constraint(const lar_term *t, lconstraint_kind kind, const mpq& right_side) : lar_base_constraint(kind, right_side), m_term(t) { }
mpq get_free_coeff_of_left_side() const override { return m_term->m_v;}
// mpq get_free_coeff_of_left_side() const override { return m_term->m_v;}
};

View file

@ -817,7 +817,7 @@ public:
}
const bool column_is_bounded(unsigned j) const {
bool column_is_bounded(unsigned j) const {
switch(m_column_types()[j]) {
case column_type::fixed:
case column_type::boxed:

View file

@ -27,7 +27,7 @@ void clear() {lp_assert(false); // not implemented
}
lar_solver::lar_solver() : m_status(lp_status::OPTIMAL),
lar_solver::lar_solver() : m_status(lp_status::UNKNOWN),
m_infeasible_column_index(-1),
m_terms_start_index(1000000),
m_mpq_lar_core_solver(m_settings, *this),
@ -71,7 +71,7 @@ bool lar_solver::sizes_are_correct() const {
}
void lar_solver::print_implied_bound(const implied_bound& be, std::ostream & out) const {
std::ostream& lar_solver::print_implied_bound(const implied_bound& be, std::ostream & out) const {
out << "implied bound\n";
unsigned v = be.m_j;
if (is_term(v)) {
@ -83,6 +83,7 @@ void lar_solver::print_implied_bound(const implied_bound& be, std::ostream & out
}
out << " " << lconstraint_kind_string(be.kind()) << " " << be.m_bound << std::endl;
out << "end of implied bound" << std::endl;
return out;
}
bool lar_solver::implied_bound_is_correctly_explained(implied_bound const & be, const vector<std::pair<mpq, unsigned>> & explanation) const {
@ -136,7 +137,7 @@ bool lar_solver::implied_bound_is_correctly_explained(implied_bound const & be,
kind = static_cast<lconstraint_kind>(-kind);
}
rs_of_evidence /= ratio;
rs_of_evidence += t->m_v * ratio;
// rs_of_evidence += t->m_v * ratio;
}
return kind == be.kind() && rs_of_evidence == be.m_bound;
@ -376,7 +377,7 @@ void lar_solver::pop(unsigned k) {
m_settings.simplex_strategy() = m_simplex_strategy;
lp_assert(sizes_are_correct());
lp_assert((!m_settings.use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
m_status = m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()? lp_status::OPTIMAL: lp_status::UNKNOWN;
set_status(lp_status::UNKNOWN);
}
vector<constraint_index> lar_solver::get_all_constraint_indices() const {
@ -601,7 +602,7 @@ void lar_solver::register_monoid_in_map(std::unordered_map<var_index, mpq> & coe
void lar_solver::substitute_terms_in_linear_expression(const vector<std::pair<mpq, var_index>>& left_side_with_terms,
vector<std::pair<mpq, var_index>> &left_side, mpq & free_coeff) const {
vector<std::pair<mpq, var_index>> &left_side) const {
std::unordered_map<var_index, mpq> coeffs;
for (auto & t : left_side_with_terms) {
unsigned j = t.second;
@ -612,7 +613,6 @@ void lar_solver::substitute_terms_in_linear_expression(const vector<std::pair<mp
for (auto & p : term.coeffs()){
register_monoid_in_map(coeffs, t.first * p.second , p.first);
}
free_coeff += t.first * term.m_v;
}
}
@ -781,7 +781,7 @@ void lar_solver::solve_with_core_solver() {
update_x_and_inf_costs_for_columns_with_changed_bounds();
m_mpq_lar_core_solver.solve();
set_status(m_mpq_lar_core_solver.m_r_solver.get_status());
lp_assert(m_status != lp_status::OPTIMAL || all_constraints_hold());
lp_assert((m_settings.random_next() % 100) != 0 || m_status != lp_status::OPTIMAL || all_constraints_hold());
}
@ -909,13 +909,8 @@ bool lar_solver::try_to_set_fixed(column_info<mpq> & ci) {
return false;
}
column_type lar_solver::get_column_type(const column_info<mpq> & ci) {
auto ret = ci.get_column_type_no_flipping();
if (ret == column_type::boxed) { // changing boxed to fixed because of the no span
if (ci.get_lower_bound() == ci.get_upper_bound())
ret = column_type::fixed;
}
return ret;
column_type lar_solver::get_column_type(unsigned j) const{
return m_mpq_lar_core_solver.m_column_types[j];
}
std::string lar_solver::get_column_name(unsigned j) const {
@ -954,9 +949,9 @@ bool lar_solver::constraint_holds(const lar_base_constraint & constr, std::unord
mpq left_side_val = get_left_side_val(constr, var_map);
switch (constr.m_kind) {
case LE: return left_side_val <= constr.m_right_side;
case LT: return left_side_val < constr.m_right_side;
case LT: return left_side_val < constr.m_right_side;
case GE: return left_side_val >= constr.m_right_side;
case GT: return left_side_val > constr.m_right_side;
case GT: return left_side_val > constr.m_right_side;
case EQ: return left_side_val == constr.m_right_side;
default:
lp_unreachable();
@ -975,8 +970,10 @@ bool lar_solver::the_relations_are_of_same_type(const vector<std::pair<mpq, unsi
flip_kind(m_constraints[con_ind]->m_kind);
if (kind == GT || kind == LT)
strict = true;
if (kind == GE || kind == GT) n_of_G++;
else if (kind == LE || kind == LT) n_of_L++;
if (kind == GE || kind == GT)
n_of_G++;
else if (kind == LE || kind == LT)
n_of_L++;
}
the_kind_of_sum = n_of_G ? GE : (n_of_L ? LE : EQ);
if (strict)
@ -1116,7 +1113,7 @@ bool lar_solver::has_upper_bound(var_index var, constraint_index& ci, mpq& value
bool lar_solver::has_value(var_index var, mpq& value) const {
if (is_term(var)) {
lar_term const& t = get_term(var);
value = t.m_v;
value = 0;
for (auto const& cv : t) {
impq const& r = get_column_value(cv.var());
if (!numeric_traits<mpq>::is_zero(r.y)) return false;
@ -1177,6 +1174,7 @@ void lar_solver::get_model(std::unordered_map<var_index, mpq> & variable_values)
std::unordered_set<impq> set_of_different_pairs;
std::unordered_set<mpq> set_of_different_singles;
delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(delta);
TRACE("get_model", tout << "delta=" << delta << "size = " << m_mpq_lar_core_solver.m_r_x.size() << std::endl;);
for (i = 0; i < m_mpq_lar_core_solver.m_r_x.size(); i++ ) {
const numeric_pair<mpq> & rp = m_mpq_lar_core_solver.m_r_x[i];
set_of_different_pairs.insert(rp);
@ -1208,42 +1206,40 @@ std::string lar_solver::get_variable_name(var_index vi) const {
}
// ********** print region start
void lar_solver::print_constraint(constraint_index ci, std::ostream & out) const {
std::ostream& lar_solver::print_constraint(constraint_index ci, std::ostream & out) const {
if (ci >= m_constraints.size()) {
out << "constraint " << T_to_string(ci) << " is not found";
out << std::endl;
return;
return out;
}
print_constraint(m_constraints[ci], out);
return print_constraint(m_constraints[ci], out);
}
void lar_solver::print_constraints(std::ostream& out) const {
std::ostream& lar_solver::print_constraints(std::ostream& out) const {
out << "number of constraints = " << m_constraints.size() << std::endl;
for (auto c : m_constraints) {
print_constraint(c, out);
}
return out;
}
void lar_solver::print_terms(std::ostream& out) const {
std::ostream& lar_solver::print_terms(std::ostream& out) const {
for (auto it : m_terms) {
print_term(*it, out);
out << "\n";
print_term(*it, out) << "\n";
}
return out;
}
void lar_solver::print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const {
std::ostream& lar_solver::print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const {
print_linear_combination_of_column_indices(c->get_left_side_coefficients(), out);
mpq free_coeff = c->get_free_coeff_of_left_side();
if (!is_zero(free_coeff))
out << " + " << free_coeff;
return out;
}
void lar_solver::print_term(lar_term const& term, std::ostream & out) const {
if (!numeric_traits<mpq>::is_zero(term.m_v)) {
out << term.m_v << " + ";
}
std::ostream& lar_solver::print_term(lar_term const& term, std::ostream & out) const {
bool first = true;
for (const auto p : term) {
mpq val = p.coeff();
@ -1263,14 +1259,12 @@ void lar_solver::print_term(lar_term const& term, std::ostream & out) const {
out << T_to_string(val);
out << this->get_column_name(p.var());
}
return out;
}
void lar_solver::print_term_as_indices(lar_term const& term, std::ostream & out) const {
if (!numeric_traits<mpq>::is_zero(term.m_v)) {
out << term.m_v << " + ";
}
print_linear_combination_of_column_indices_only(term.coeffs_as_vector(), out);
std::ostream& lar_solver::print_term_as_indices(lar_term const& term, std::ostream & out) const {
print_linear_combination_of_column_indices_only(term, out);
return out;
}
mpq lar_solver::get_left_side_val(const lar_base_constraint & cns, const std::unordered_map<var_index, mpq> & var_map) const {
@ -1284,9 +1278,10 @@ mpq lar_solver::get_left_side_val(const lar_base_constraint & cns, const std::u
return ret;
}
void lar_solver::print_constraint(const lar_base_constraint * c, std::ostream & out) const {
std::ostream& lar_solver::print_constraint(const lar_base_constraint * c, std::ostream & out) const {
print_left_side_of_constraint(c, out);
out << " " << lconstraint_kind_string(c->m_kind) << " " << c->m_right_side << std::endl;
return out;
}
void lar_solver::fill_var_set_for_random_update(unsigned sz, var_index const * vars, vector<unsigned>& column_list) {
@ -1492,7 +1487,7 @@ bool lar_solver::term_is_int(const lar_term * t) const {
for (auto const & p : t->m_coeffs)
if (! (column_is_int(p.first) && p.second.is_int()))
return false;
return t->m_v.is_int();
return true;
}
bool lar_solver::var_is_int(var_index v) const {
@ -1593,17 +1588,13 @@ void lar_solver::add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
}
var_index lar_solver::add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
const mpq &m_v) {
push_and_register_term(new lar_term(coeffs, m_v));
var_index lar_solver::add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs) {
push_and_register_term(new lar_term(coeffs));
return m_terms_start_index + m_terms.size() - 1;
}
#if Z3DEBUG_CHECK_UNIQUE_TERMS
bool lar_solver::term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs, const mpq& v) {
if (coeffs.empty()) {
return is_zero(v);
}
bool lar_solver::term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs) {
for (const auto & p : coeffs) {
if (column_is_real(p.second))
@ -1638,12 +1629,11 @@ void lar_solver::push_and_register_term(lar_term* t) {
}
// terms
var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
const mpq &m_v) {
var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs) {
if (strategy_is_undecided())
return add_term_undecided(coeffs, m_v);
return add_term_undecided(coeffs);
push_and_register_term(new lar_term(coeffs, m_v));
push_and_register_term(new lar_term(coeffs));
unsigned adjusted_term_index = m_terms.size() - 1;
var_index ret = m_terms_start_index + adjusted_term_index;
if (use_tableau() && !coeffs.empty()) {
@ -1656,7 +1646,7 @@ var_index lar_solver::add_term(const vector<std::pair<mpq, var_index>> & coeffs,
}
void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) {
TRACE("dump_terms", print_term(*term, tout); tout << std::endl;);
TRACE("dump_terms", print_term(*term, tout) << std::endl;);
register_new_ext_var_index(term_ext_index, term_is_int(term));
// j will be a new variable
unsigned j = A_r().column_count();
@ -1738,9 +1728,8 @@ void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_k
// lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int());
unsigned term_j;
if (m_var_register.external_is_used(j, term_j)) {
mpq rs = right_side - m_terms[adjusted_term_index]->m_v;
m_constraints.push_back(new lar_term_constraint(m_terms[adjusted_term_index], kind, right_side));
update_column_type_and_bound(term_j, kind, rs, ci);
update_column_type_and_bound(term_j, kind, right_side, ci);
}
else {
add_constraint_from_term_and_create_new_column_row(j, m_terms[adjusted_term_index], kind, right_side);
@ -1749,11 +1738,10 @@ void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_k
constraint_index lar_solver::add_constraint(const vector<std::pair<mpq, var_index>>& left_side_with_terms, lconstraint_kind kind_par, const mpq& right_side_parm) {
vector<std::pair<mpq, var_index>> left_side;
mpq rs = -right_side_parm;
substitute_terms_in_linear_expression(left_side_with_terms, left_side, rs);
unsigned term_index = add_term(left_side, zero_of_type<mpq>());
substitute_terms_in_linear_expression(left_side_with_terms, left_side);
unsigned term_index = add_term(left_side);
constraint_index ci = m_constraints.size();
add_var_bound_on_constraint_for_term(term_index, kind_par, -rs, ci);
add_var_bound_on_constraint_for_term(term_index, kind_par, right_side_parm, ci);
return ci;
}
@ -1762,7 +1750,7 @@ void lar_solver::add_constraint_from_term_and_create_new_column_row(unsigned ter
add_row_from_term_no_constraint(term, term_j);
unsigned j = A_r().column_count() - 1;
update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size());
update_column_type_and_bound(j, kind, right_side, m_constraints.size());
m_constraints.push_back(new lar_term_constraint(term, kind, right_side));
lp_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
}
@ -1964,7 +1952,6 @@ void lar_solver::update_boxed_column_type_and_bound(var_index j, lconstraint_kin
if (up < m_mpq_lar_core_solver.m_r_lower_bounds[j]) {
m_status = lp_status::INFEASIBLE;
lp_assert(false);
m_infeasible_column_index = j;
}
else {
@ -2261,6 +2248,7 @@ void lar_solver::set_cut_strategy(unsigned cut_frequency) {
}
}
} // namespace lp

View file

@ -94,7 +94,6 @@ class lar_solver : public column_namer {
var_register m_var_register;
stacked_vector<ul_pair> m_columns_to_ul_pairs;
vector<lar_base_constraint*> m_constraints;
private:
stacked_value<unsigned> m_constraint_count;
// the set of column indices j such that bounds have changed for j
int_set m_columns_with_changed_bound;
@ -164,13 +163,11 @@ public:
// terms
var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
const mpq &m_v);
var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs);
var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
const mpq &m_v);
var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs);
bool term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs, const mpq& v);
bool term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & coeffs);
void push_and_register_term(lar_term* t);
void add_row_for_term(const lar_term * term, unsigned term_ext_index);
@ -208,7 +205,10 @@ public:
void update_lower_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci);
void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci);
//end of init region
lp_settings & settings();
lp_settings const & settings() const;
@ -227,9 +227,7 @@ public:
bool use_lu() const;
bool sizes_are_correct() const;
void print_implied_bound(const implied_bound& be, std::ostream & out) const;
bool implied_bound_is_correctly_explained(implied_bound const & be, const vector<std::pair<mpq, unsigned>> & explanation) const;
void analyze_new_bounds_on_row(
@ -238,8 +236,7 @@ public:
void analyze_new_bounds_on_row_tableau(
unsigned row_index,
bound_propagator & bp
);
bound_propagator & bp);
void substitute_basis_var_in_terms_for_row(unsigned i);
@ -331,7 +328,7 @@ public:
void substitute_terms_in_linear_expression( const vector<std::pair<mpq, var_index>>& left_side_with_terms,
vector<std::pair<mpq, var_index>> &left_side, mpq & free_coeff) const;
vector<std::pair<mpq, var_index>> &left_side) const;
void detect_rows_of_bound_change_column_for_nbasic_column(unsigned j);
@ -397,9 +394,9 @@ public:
bool try_to_set_fixed(column_info<mpq> & ci);
column_type get_column_type(const column_info<mpq> & ci);
column_type get_column_type(unsigned j) const;
std::string get_column_name(unsigned j) const;
std::string get_column_name(unsigned j) const override;
bool all_constrained_variables_are_registered(const vector<std::pair<mpq, var_index>>& left_side);
@ -436,30 +433,33 @@ public:
int inf_sign) const;
void get_model(std::unordered_map<var_index, mpq> & variable_values) const;
void get_model_do_not_care_about_diff_vars(std::unordered_map<var_index, mpq> & variable_values) const;
std::string get_variable_name(var_index vi) const;
// ********** print region start
void print_constraint(constraint_index ci, std::ostream & out) const;
// print utilities
void print_constraints(std::ostream& out) const ;
std::ostream& print_constraint(constraint_index ci, std::ostream & out) const;
void print_terms(std::ostream& out) const;
std::ostream& print_constraints(std::ostream& out) const ;
void print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const;
std::ostream& print_terms(std::ostream& out) const;
void print_term(lar_term const& term, std::ostream & out) const;
std::ostream& print_left_side_of_constraint(const lar_base_constraint * c, std::ostream & out) const;
void print_term_as_indices(lar_term const& term, std::ostream & out) const;
std::ostream& print_term(lar_term const& term, std::ostream & out) const;
std::ostream& print_term_as_indices(lar_term const& term, std::ostream & out) const;
std::ostream& print_constraint(const lar_base_constraint * c, std::ostream & out) const;
std::ostream& print_implied_bound(const implied_bound& be, std::ostream & out) const;
mpq get_left_side_val(const lar_base_constraint & cns, const std::unordered_map<var_index, mpq> & var_map) const;
void print_constraint(const lar_base_constraint * c, std::ostream & out) const;
void fill_var_set_for_random_update(unsigned sz, var_index const * vars, vector<unsigned>& column_list);
void random_update(unsigned sz, var_index const * vars);

View file

@ -21,9 +21,9 @@
#include "util/lp/indexed_vector.h"
namespace lp {
struct lar_term {
// the term evaluates to sum of m_coeffs + m_v
// the term evaluates to sum of m_coeffs
std::unordered_map<unsigned, mpq> m_coeffs;
mpq m_v;
// mpq m_v;
lar_term() {}
void add_monomial(const mpq& c, unsigned j) {
auto it = m_coeffs.find(j);
@ -37,7 +37,7 @@ struct lar_term {
}
bool is_empty() const {
return m_coeffs.size() == 0 && is_zero(m_v);
return m_coeffs.size() == 0; // && is_zero(m_v);
}
unsigned size() const { return static_cast<unsigned>(m_coeffs.size()); }
@ -46,8 +46,7 @@ struct lar_term {
return m_coeffs;
}
lar_term(const vector<std::pair<mpq, unsigned>>& coeffs,
const mpq & v) : m_v(v) {
lar_term(const vector<std::pair<mpq, unsigned>>& coeffs) {
for (const auto & p : coeffs) {
add_monomial(p.first, p.second);
}
@ -87,7 +86,7 @@ struct lar_term {
template <typename T>
T apply(const vector<T>& x) const {
T ret = T(m_v);
T ret(0);
for (const auto & t : m_coeffs) {
ret += t.second * x[t.first];
}
@ -96,7 +95,6 @@ struct lar_term {
void clear() {
m_coeffs.clear();
m_v = zero_of_type<mpq>();
}
struct ival {

View file

@ -17,10 +17,6 @@ const impq & bound_propagator::get_upper_bound(unsigned j) const {
}
void bound_propagator::try_add_bound(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);
if (m_lar_solver.is_term(j)) {
// lp treats terms as not having a free coefficient, restoring it below for the outside consumption
v += m_lar_solver.get_term(j).m_v;
}
lconstraint_kind kind = is_low? GE : LE;
if (strict)

View file

@ -84,6 +84,7 @@ template lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::lp_core_s
const vector<lp::numeric_pair<lp::mpq> >&);
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::print_statistics_with_cost_and_check_that_the_time_is_over(lp::numeric_pair<lp::mpq>, std::ostream&);
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::snap_xN_to_bounds_and_fill_xB();
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::solve_Ax_eq_b();
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::solve_Bd(unsigned int);
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::update_basis_and_x(int, int, lp::numeric_pair<lp::mpq> const&);
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::update_x(unsigned int, const lp::numeric_pair<lp::mpq>&);

View file

@ -577,7 +577,7 @@ public:
}
void print_column_info(unsigned j, std::ostream & out) const {
out << "j = " << j << ", name = "<< column_name(j);
out << "j = " << j << ",\tname = "<< column_name(j) << "\t";
switch (m_column_types[j]) {
case column_type::fixed:
case column_type::boxed:
@ -596,11 +596,11 @@ public:
lp_assert(false);
}
// out << "basis heading = " << m_basis_heading[j] << std::endl;
out << " x = " << m_x[j];
out << "\tx = " << m_x[j];
if (m_basis_heading[j] >= 0)
out << " base\n";
else
out << " nbas\n";
out << " \n";
}
bool column_is_free(unsigned j) const { return this->m_column_type[j] == free; }

View file

@ -1238,6 +1238,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::print_column
break;
case column_type::free_column:
out << "( _" << this->m_x[j] << "_)" << std::endl;
break;
default:
lp_unreachable();
}

View file

@ -357,7 +357,7 @@ public:
}
#ifdef Z3DEBUG
static unsigned ddd; // used for debugging
static unsigned ddd; // used for debugging
#endif
}; // end of lp_settings class
@ -433,7 +433,15 @@ inline void ensure_increasing(vector<unsigned> & v) {
}
}
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(n));
return n.x - floor(n.x);
}
inline static mpq fractional_part(const mpq & n) {
return n - floor(n);
}
#if Z3DEBUG
bool D();

View file

@ -24,7 +24,7 @@ Revision History:
namespace lp {
template <typename T, typename X> column_info<T> * lp_solver<T, X>::get_or_create_column_info(unsigned column) {
auto it = m_map_from_var_index_to_column_info.find(column);
return (it == m_map_from_var_index_to_column_info.end())? (m_map_from_var_index_to_column_info[column] = new column_info<T>(static_cast<unsigned>(-1))) : it->second;
return (it == m_map_from_var_index_to_column_info.end())? (m_map_from_var_index_to_column_info[column] = new column_info<T>()) : it->second;
}
template <typename T, typename X>

View file

@ -50,11 +50,34 @@ bool contains(const std::unordered_map<A, B> & map, const A& key) {
namespace lp {
inline void throw_exception(std::string && str) {
throw default_exception(std::move(str));
template <typename T>
void print_linear_combination_of_column_indices_only(const T & coeffs, std::ostream & out) {
bool first = true;
for (const auto & it : coeffs) {
auto val = it.coeff();
if (first) {
first = false;
} else {
if (val.is_pos()) {
out << " + ";
} else {
out << " - ";
val = -val;
}
}
if (val == 1)
out << " ";
else
out << T_to_string(val);
out << "x" << it.var();
}
typedef z3_exception exception;
}
inline void throw_exception(std::string && str) {
throw default_exception(std::move(str));
}
typedef z3_exception exception;
#define lp_assert(_x_) { SASSERT(_x_); }
inline void lp_unreachable() { lp_assert(false); }

View file

@ -154,7 +154,7 @@ namespace nra {
polynomial::polynomial* ps[1] = { p };
bool even[1] = { false };
nlsat::literal lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, even);
m_nlsat->mk_clause(1, &lit, 0);
m_nlsat->mk_clause(1, &lit, nullptr);
}
void add_constraint(unsigned idx) {

View file

@ -57,10 +57,10 @@ public:
template <> class numeric_traits<int> {
public:
static bool precise() { return true; }
static int const zero() { return 0; }
static int const one() { return 1; }
static int zero() { return 0; }
static int one() { return 1; }
static bool is_zero(int v) { return v == 0; }
static double const get_double(int const & d) { return d; }
static double get_double(int const & d) { return d; }
static bool is_int(int) {return true;}
static bool is_pos(int j) {return j > 0;}
static bool is_neg(int j) {return j < 0;}

View file

@ -27,8 +27,9 @@ Revision History:
# define __has_builtin(x) 0
#endif
#ifdef __GNUC__
# if (__GNUC__ * 100 + __GNUC_MINOR__) >= 409 || __has_builtin(returns_nonnull)
# if ((__GNUC__ * 100 + __GNUC_MINOR__) >= 409 || __has_builtin(returns_nonnull)) && !defined(__INTEL_COMPILER)
# define GCC_RET_NON_NULL __attribute__((returns_nonnull))
# else
# define GCC_RET_NON_NULL

View file

@ -152,7 +152,7 @@ class mpz_manager {
// make sure that n is a big number and has capacity equal to at least c.
void allocate_if_needed(mpz & n, unsigned c) {
c = std::max(c, m_init_cell_capacity);
if (m_init_cell_capacity > c) c = m_init_cell_capacity;
if (n.m_ptr == nullptr || capacity(n) < c) {
deallocate(n);
n.m_val = 1;

View file

@ -115,7 +115,7 @@ public:
*/
T * steal() {
T * r = m_obj;
m_obj = 0;
m_obj = nullptr;
return r;
}
};

View file

@ -27,7 +27,7 @@ Revision History:
// Windows
#include<windows.h>
#elif defined(__APPLE__) && defined(__MACH__)
// Mac OS X
// macOS
#include<mach/mach.h>
#include<mach/clock.h>
#include<sys/time.h>
@ -59,7 +59,7 @@ struct scoped_timer::imp {
HANDLE m_timer;
bool m_first;
#elif defined(__APPLE__) && defined(__MACH__)
// Mac OS X
// macOS
pthread_t m_thread_id;
pthread_attr_t m_attributes;
unsigned m_interval;
@ -89,7 +89,7 @@ struct scoped_timer::imp {
}
}
#elif defined(__APPLE__) && defined(__MACH__)
// Mac OS X
// macOS
static void * thread_func(void * arg) {
scoped_timer::imp * st = static_cast<scoped_timer::imp*>(arg);
@ -153,7 +153,7 @@ struct scoped_timer::imp {
ms,
WT_EXECUTEINTIMERTHREAD);
#elif defined(__APPLE__) && defined(__MACH__)
// Mac OS X
// macOS
m_interval = ms?ms:0xFFFFFFFF;
if (pthread_attr_init(&m_attributes) != 0)
throw default_exception("failed to initialize timer thread attributes");
@ -194,7 +194,7 @@ struct scoped_timer::imp {
m_timer,
INVALID_HANDLE_VALUE);
#elif defined(__APPLE__) && defined(__MACH__)
// Mac OS X
// macOS
// If the waiting-thread is not up and waiting yet,
// we can make sure that it finishes quickly by

View file

@ -23,8 +23,6 @@ Revision History:
#if defined(_WINDOWS) || defined(_CYGWIN)
// Does this redefinition work?
#define ARRAYSIZE_TEMP ARRAYSIZE
#undef ARRAYSIZE
#include <windows.h>
@ -68,13 +66,11 @@ public:
}
};
#undef ARRAYSIZE
#define ARRAYSIZE ARRAYSIZE_TEMP
#undef max
#undef min
#elif defined(__APPLE__) && defined (__MACH__) // Mac OS X
#elif defined(__APPLE__) && defined (__MACH__) // macOS
#include<mach/mach.h>
#include<mach/clock.h>

View file

@ -44,6 +44,8 @@ class string_buffer {
m_buffer = new_buffer;
}
static const unsigned c_buffer_size = 24;
public:
string_buffer():
m_buffer(m_initial_buffer),
@ -80,29 +82,21 @@ public:
}
void append(int n) {
char buffer[24];
#ifdef _WINDOWS
sprintf_s(buffer, ARRAYSIZE(buffer), "%d", n);
#else
sprintf(buffer, "%d", n);
#endif
char buffer[c_buffer_size];
SPRINTF_D(buffer, n);
append(buffer);
}
void append(unsigned n) {
char buffer[24];
#ifdef _WINDOWS
sprintf_s(buffer, ARRAYSIZE(buffer), "%d", n);
#else
sprintf(buffer, "%d", n);
#endif
char buffer[c_buffer_size];
SPRINTF_U(buffer, n);
append(buffer);
}
void append(long n) {
char buffer[24];
char buffer[c_buffer_size];
#ifdef _WINDOWS
sprintf_s(buffer, ARRAYSIZE(buffer), "%ld", n);
sprintf_s(buffer, Z3_ARRAYSIZE(buffer), "%ld", n);
#else
sprintf(buffer, "%ld", n);
#endif

View file

@ -122,8 +122,10 @@ public:
TRACE("union_find", tout << "merging " << r1 << " " << r2 << "\n";);
if (r1 == r2)
return;
if (m_size[r1] > m_size[r2])
if (m_size[r1] > m_size[r2]) {
std::swap(r1, r2);
std::swap(v1, v2);
}
m_ctx.merge_eh(r2, r1, v2, v1);
m_find[r1] = r2;
m_size[r2] += m_size[r1];

View file

@ -48,13 +48,19 @@ static_assert(sizeof(int64_t) == 8, "64 bits");
#ifdef _WINDOWS
#define SSCANF sscanf_s
#define SPRINTF sprintf_s
// #define SPRINTF sprintf_s
#define SPRINTF_D(_buffer_, _i_) sprintf_s(_buffer_, Z3_ARRAYSIZE(_buffer_), "%d", _i_)
#define SPRINTF_U(_buffer_, _u_) sprintf_s(_buffer_, Z3_ARRAYSIZE(_buffer_), "%u", _u_)
#define _Exit exit
#else
#define SSCANF sscanf
#define SPRINTF sprintf
// #define SPRINTF sprintf
#define SPRINTF_D(_buffer_, _i_) sprintf(_buffer_, "%d", _i_)
#define SPRINTF_U(_buffer_, _u_) sprintf(_buffer_, "%u", _u_)
#endif
#define VEC2PTR(_x_) ((_x_).size() ? &(_x_)[0] : 0)
#ifdef _WINDOWS
@ -142,9 +148,7 @@ static inline uint64_t shift_left(uint64_t x, uint64_t y) {
template<class T, size_t N> char (*ArraySizer(T (&)[N]))[N];
// For determining the length of an array. See ARRAYSIZE() macro. This function is never actually called.
#ifndef ARRAYSIZE
#define ARRAYSIZE(a) sizeof(*ArraySizer(a))
#endif
#define Z3_ARRAYSIZE(a) sizeof(*ArraySizer(a))
template<typename IT>
void display(std::ostream & out, const IT & begin, const IT & end, const char * sep, bool & first) {