z3-z3-4.13.0.src.test.hilbert_basis.cpp Maven / Gradle / Ivy
The newest version!
/*++
Copyright (c) 2015 Microsoft Corporation
--*/
#include "math/hilbert/hilbert_basis.h"
#include "ast/ast_pp.h"
#include "ast/reg_decl_plugins.h"
#include "ast/arith_decl_plugin.h"
#include "tactic/smtlogics/quant_tactics.h"
#include "tactic/tactic.h"
#include "solver/tactic2solver.h"
#include "solver/solver.h"
#include "util/rlimit.h"
#include
#include
#include
#include
static bool g_use_ordered_support = false;
static bool g_use_ordered_subsumption = false;
static bool g_use_support = false;
class hilbert_basis_validate {
ast_manager& m;
void validate_solution(hilbert_basis& hb, vector const& v, bool is_initial);
rational eval_ineq(hilbert_basis& hb, unsigned idx, vector const& v, bool& is_eq) {
vector w;
rational bound;
hb.get_ge(idx, w, bound, is_eq);
rational sum(0);
for (unsigned j = 0; j < v.size(); ++j) {
sum += w[j]*v[j];
}
sum -= bound;
return sum;
}
public:
hilbert_basis_validate(ast_manager& m);
expr_ref mk_validate(hilbert_basis& hb);
};
hilbert_basis_validate::hilbert_basis_validate(ast_manager& m):
m(m) {
}
void hilbert_basis_validate::validate_solution(hilbert_basis& hb, vector const& v, bool is_initial) {
unsigned sz = hb.get_num_ineqs();
rational bound;
for (unsigned i = 0; i < sz; ++i) {
bool is_eq;
vector w;
rational bound;
hb.get_ge(i, w, bound, is_eq);
rational sum(0);
for (unsigned j = 0; j < v.size(); ++j) {
sum += w[j]*v[j];
}
if (sum >= bound && !is_eq) {
continue;
}
if (sum == bound && is_eq) {
continue;
}
// homogeneous solutions should be non-negative.
if (!is_initial && sum.is_nonneg()) {
continue;
}
// validation failed.
std::cout << "validation failed\n";
std::cout << "constraint: ";
for (unsigned j = 0; j < v.size(); ++j) {
std::cout << v[j] << " ";
}
std::cout << (is_eq?" = ":" >= ") << bound << "\n";
std::cout << "vector: ";
for (unsigned j = 0; j < w.size(); ++j) {
std::cout << w[j] << " ";
}
std::cout << "\n";
std::cout << "sum: " << sum << "\n";
}
}
expr_ref hilbert_basis_validate::mk_validate(hilbert_basis& hb) {
arith_util a(m);
unsigned sz = hb.get_basis_size();
vector v;
// check that claimed solution really satisfies inequalities:
for (unsigned i = 0; i < sz; ++i) {
bool is_initial;
hb.get_basis_solution(i, v, is_initial);
validate_solution(hb, v, is_initial);
}
// check that solutions satisfying inequalities are in solution.
// build a formula that says solutions to linear inequalities
// coincide with linear combinations of basis.
vector offsets, increments;
expr_ref_vector xs(m), vars(m);
expr_ref var(m);
svector names;
sort_ref_vector sorts(m);
#define mk_mul(_r,_x) (_r.is_one()?((expr*)_x):((expr*)a.mk_mul(a.mk_numeral(_r,true),_x)))
for (unsigned i = 0; i < sz; ++i) {
bool is_initial;
hb.get_basis_solution(i, v, is_initial);
for (unsigned j = 0; xs.size() < v.size(); ++j) {
xs.push_back(m.mk_fresh_const("x", a.mk_int()));
}
if (is_initial) {
expr_ref_vector tmp(m);
for (unsigned j = 0; j < v.size(); ++j) {
tmp.push_back(a.mk_numeral(v[j], true));
}
offsets.push_back(tmp);
}
else {
var = m.mk_var(vars.size(), a.mk_int());
expr_ref_vector tmp(m);
for (unsigned j = 0; j < v.size(); ++j) {
tmp.push_back(mk_mul(v[j], var));
}
std::stringstream name;
name << "u" << i;
increments.push_back(tmp);
vars.push_back(var);
names.push_back(symbol(name.str()));
sorts.push_back(a.mk_int());
}
}
expr_ref_vector bounds(m);
for (unsigned i = 0; i < vars.size(); ++i) {
bounds.push_back(a.mk_ge(vars[i].get(), a.mk_numeral(rational(0), true)));
}
expr_ref_vector fmls(m);
expr_ref fml(m), fml1(m), fml2(m);
for (unsigned i = 0; i < offsets.size(); ++i) {
expr_ref_vector eqs(m);
eqs.append(bounds);
for (unsigned j = 0; j < xs.size(); ++j) {
expr_ref_vector sum(m);
sum.push_back(offsets[i][j].get());
for (unsigned k = 0; k < increments.size(); ++k) {
sum.push_back(increments[k][j].get());
}
eqs.push_back(m.mk_eq(xs[j].get(), a.mk_add(sum.size(), sum.data())));
}
fml = m.mk_and(eqs.size(), eqs.data());
if (!names.empty()) {
fml = m.mk_exists(names.size(), sorts.data(), names.data(), fml);
}
fmls.push_back(fml);
}
fml1 = m.mk_or(fmls.size(), fmls.data());
fmls.reset();
sz = hb.get_num_ineqs();
for (unsigned i = 0; i < sz; ++i) {
bool is_eq;
vector w;
rational bound;
hb.get_ge(i, w, bound, is_eq);
expr_ref_vector sum(m);
for (unsigned j = 0; j < w.size(); ++j) {
if (!w[j].is_zero()) {
sum.push_back(mk_mul(w[j], xs[j].get()));
}
}
expr_ref lhs(m), rhs(m);
lhs = a.mk_add(sum.size(), sum.data());
rhs = a.mk_numeral(bound, true);
if (is_eq) {
fmls.push_back(a.mk_eq(lhs, rhs));
}
else {
fmls.push_back(a.mk_ge(lhs, rhs));
}
}
fml2 = m.mk_and(fmls.size(), fmls.data());
fml = m.mk_eq(fml1, fml2);
bounds.reset();
for (unsigned i = 0; i < xs.size(); ++i) {
if (!hb.get_is_int(i)) {
bounds.push_back(a.mk_ge(xs[i].get(), a.mk_numeral(rational(0), true)));
}
}
if (!bounds.empty()) {
fml = m.mk_implies(m.mk_and(bounds.size(), bounds.data()), fml);
}
return fml;
}
hilbert_basis* g_hb = nullptr;
static double g_start_time;
static void display_statistics(hilbert_basis& hb) {
double time = static_cast(clock()) - g_start_time;
statistics st;
hb.collect_statistics(st);
st.display(std::cout);
std::cout << "time: " << (time / CLOCKS_PER_SEC) << " secs\n";
}
static void on_ctrl_c(int) {
signal (SIGINT, SIG_DFL);
display_statistics(*g_hb);
raise(SIGINT);
}
#if 0
static void validate_sat(hilbert_basis& hb) {
ast_manager m;
reg_decl_plugins(m);
hilbert_basis_validate val(m);
expr_ref fml = val.mk_validate(hb);
return;
std::cout << mk_pp(fml, m) << "\n";
fml = m.mk_not(fml);
params_ref p;
tactic_ref tac = mk_lra_tactic(m, p);
ref sol = mk_tactic2solver(m, tac.get(), p);
sol->assert_expr(fml);
lbool r = sol->check_sat(0,0);
std::cout << r << "\n";
}
#endif
static void saturate_basis(hilbert_basis& hb) {
signal(SIGINT, on_ctrl_c);
g_hb = &hb;
g_start_time = static_cast(clock());
hb.set_use_ordered_support(g_use_ordered_support);
hb.set_use_support(g_use_support);
hb.set_use_ordered_subsumption(g_use_ordered_subsumption);
lbool is_sat = hb.saturate();
switch(is_sat) {
case l_true:
std::cout << "sat\n";
hb.display(std::cout);
//validate_sat(hb);
break;
case l_false:
std::cout << "unsat\n";
break;
case l_undef:
std::cout << "undef\n";
break;
}
display_statistics(hb);
}
/**
n - number of variables.
k - subset of variables to be non-zero
bound - numeric value of upper and lower bound
num_ineqs - number of inequalities to create
*/
static void gorrila_test(unsigned seed, unsigned n, unsigned k, unsigned bound, unsigned num_ineqs) {
std::cout << "Gorrila test\n";
random_gen rand(seed);
reslimit rl;
hilbert_basis hb(rl);
ENSURE(0 < bound);
ENSURE(k <= n);
int ibound = static_cast(bound);
for (unsigned i = 0; i < num_ineqs; ++i) {
vector nv;
nv.resize(n);
rational a0;
unsigned num_selected = 0;
while (num_selected < k) {
unsigned s = rand(n);
if (nv[s].is_zero()) {
nv[s] = rational(ibound - static_cast(rand(2*bound+1)));
if (!nv[s].is_zero()) {
++num_selected;
}
}
}
a0 = rational(ibound - static_cast(rand(2*bound+1)));
hb.add_ge(nv, a0);
}
hb.display(std::cout << "Saturate\n");
saturate_basis(hb);
}
static vector vec(int i, int j) {
vector nv;
nv.resize(2);
nv[0] = rational(i);
nv[1] = rational(j);
return nv;
}
static vector vec(int i, int j, int k) {
vector nv;
nv.resize(3);
nv[0] = rational(i);
nv[1] = rational(j);
nv[2] = rational(k);
return nv;
}
static vector vec(int i, int j, int k, int l) {
vector nv;
nv.resize(4);
nv[0] = rational(i);
nv[1] = rational(j);
nv[2] = rational(k);
nv[3] = rational(l);
return nv;
}
static vector vec(int i, int j, int k, int l, int m) {
vector nv;
nv.resize(5);
nv[0] = rational(i);
nv[1] = rational(j);
nv[2] = rational(k);
nv[3] = rational(l);
nv[4] = rational(m);
return nv;
}
static vector vec(int i, int j, int k, int l, int x, int y, int z) {
vector nv;
nv.resize(7);
nv[0] = rational(i);
nv[1] = rational(j);
nv[2] = rational(k);
nv[3] = rational(l);
nv[4] = rational(x);
nv[5] = rational(y);
nv[6] = rational(z);
return nv;
}
// example 9, Ajili, Contenjean
// x + y - 2z = 0
// x - z = 0
// -y + z <= 0
static void tst1() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_eq(vec(1,1,-2));
hb.add_eq(vec(1,0,-1));
hb.add_le(vec(0,1,-1));
saturate_basis(hb);
}
// example 10, Ajili, Contenjean
// 23x - 12y - 9z <= 0
// x - 8y - 8z <= 0
void tst2() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_eq(vec(-23,12,9));
hb.add_eq(vec(-1,8,8));
saturate_basis(hb);
}
// example 6, Ajili, Contenjean
// 3x + 2y - z - 2u <= 0
static void tst3() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec(3,2,-1,-2));
saturate_basis(hb);
}
#define R rational
// Sigma_1, table 1, Ajili, Contejean
static void tst4() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec( 0,-2, 1, 3, 2,-2, 3), R(3));
hb.add_le(vec(-1, 7, 0, 1, 3, 5,-4), R(2));
hb.add_le(vec( 0,-1, 1,-1,-1, 0, 0), R(2));
hb.add_le(vec(-2, 0, 1, 4, 0, 0,-2), R(1));
hb.add_le(vec(-3, 2,-2, 2,-4,-1, 0), R(8));
hb.add_le(vec( 3,-2, 2,-2, 4, 1, 0), R(3));
hb.add_le(vec( 1, 0, 0,-1, 0, 1, 0), R(4));
hb.add_le(vec( 1,-2, 0, 0, 0, 0, 0), R(2));
hb.add_le(vec( 1, 1, 0, 0,-1, 0, 1), R(4));
hb.add_le(vec( 1, 0, 0, 0,-1, 0, 0), R(9));
saturate_basis(hb);
}
// Sigma_2 table 1, Ajili, Contejean
static void tst5() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec( 1, 2,-1, 1), R(3));
hb.add_le(vec( 2, 4, 1, 2), R(12));
hb.add_le(vec( 1, 4, 2, 1), R(9));
hb.add_le(vec( 1, 1, 0,-1), R(10));
hb.add_le(vec( 1, 1,-1, 0), R(6));
hb.add_le(vec( 1,-1, 0, 0), R(0));
hb.add_le(vec( 0, 0, 1,-1), R(2));
saturate_basis(hb);
}
// Sigma_3 table 1, Ajili, Contejean
static void tst6() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec( 4, 3, 0), R(6));
hb.add_le(vec(-3,-4, 0), R(-1));
hb.add_le(vec( 4, 0,-3), R(3));
hb.add_le(vec(-3, 0, 4), R(7));
hb.add_le(vec( 4, 0,-3), R(23));
hb.add_le(vec( 0,-3, 4), R(11));
saturate_basis(hb);
}
// Sigma_4 table 1, Ajili, Contejean
static void tst7() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_eq(vec( 1, 1, 1, 0), R(5));
hb.add_le(vec( 2, 1, 0, 1), R(6));
hb.add_le(vec( 1, 2, 1, 1), R(7));
hb.add_le(vec( 1, 3,-1, 2), R(8));
hb.add_le(vec( 1, 2,-9,-12), R(-11));
hb.add_le(vec( 0, 0,-1, 3), R(10));
saturate_basis(hb);
}
// Sigma_5 table 1, Ajili, Contejean
static void tst8() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec( 2, 1, 1), R(2));
hb.add_le(vec( 1, 2, 3), R(5));
hb.add_le(vec( 2, 2, 3), R(6));
hb.add_le(vec( 1,-1,-3), R(-2));
saturate_basis(hb);
}
// Sigma_6 table 1, Ajili, Contejean
static void tst9() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec( 1, 2, 3), R(11));
hb.add_le(vec( 2, 2, 5), R(13));
hb.add_le(vec( 1,-1,-11), R(3));
saturate_basis(hb);
}
// Sigma_7 table 1, Ajili, Contejean
static void tst10() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec( 1,-1,-1,-3), R(2));
hb.add_le(vec(-2, 3, 3,-5), R(3));
saturate_basis(hb);
}
// Sigma_8 table 1, Ajili, Contejean
static void tst11() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec( 7,-2,11, 3, -5), R(5));
saturate_basis(hb);
}
// Sigma_9 table 1, Ajili, Contejean
static void tst12() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_eq(vec( 1,-2,-3,4), R(0));
hb.add_le(vec(100,45,-78,-67), R(0));
saturate_basis(hb);
}
// Sigma_10 table 1, Ajili, Contejean
static void tst13() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec( 23, -56, -34, 12, 11), R(0));
saturate_basis(hb);
}
// Sigma_11 table 1, Ajili, Contejean
static void tst14() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_eq(vec(1, 0, -4, 8), R(2));
hb.add_le(vec(12,19,-11,-7), R(-7));
saturate_basis(hb);
}
static void tst15() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec(1, 0), R(1));
hb.add_le(vec(0, 1), R(1));
saturate_basis(hb);
}
static void tst16() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_le(vec(1, 0), R(100));
saturate_basis(hb);
}
static void tst17() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_eq(vec(1, 0), R(0));
hb.add_eq(vec(-1, 0), R(0));
hb.add_eq(vec(0, 2), R(0));
hb.add_eq(vec(0, -2), R(0));
saturate_basis(hb);
}
static void tst18() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_eq(vec(0, 1), R(0));
hb.add_eq(vec(1, -1), R(2));
saturate_basis(hb);
}
static void tst19() {
reslimit rl;
hilbert_basis hb(rl);
hb.add_eq(vec(0, 1, 0), R(0));
hb.add_eq(vec(1, -1, 0), R(2));
saturate_basis(hb);
}
static void test_A_5_5_3() {
reslimit rl;
hilbert_basis hb(rl);
for (unsigned i = 0; i < 15; ++i) {
vector v;
for (unsigned j = 0; j < 5; ++j) {
for (unsigned k = 0; k < 15; ++k) {
v.push_back(rational(k == i));
}
}
hb.add_ge(v, R(0));
}
for (unsigned i = 1; i <= 15; ++i) {
vector v;
for (unsigned k = 1; k <= 5; ++k) {
for (unsigned l = 1; l <= 5; ++l) {
for (unsigned j = 1; j <= 3; ++j) {
bool one = ((j*k <= i) && (((i - j) % 3) == 0)); // fixme
v.push_back(rational(one));
}
}
}
hb.add_ge(v, R(0));
}
// etc.
saturate_basis(hb);
}
void tst_hilbert_basis() {
std::cout << "hilbert basis test\n";
// tst3();
// return;
g_use_ordered_support = true;
test_A_5_5_3();
return;
tst18();
return;
tst19();
return;
tst17();
if (true) {
tst1();
tst2();
tst3();
tst4();
tst4();
tst4();
tst4();
tst4();
tst4();
tst5();
tst6();
tst7();
tst8();
tst9();
tst10();
tst11();
tst12();
tst13();
tst14();
tst15();
tst16();
gorrila_test(0, 4, 3, 20, 5);
gorrila_test(1, 4, 3, 20, 5);
//gorrila_test(2, 4, 3, 20, 5);
//gorrila_test(0, 4, 2, 20, 5);
//gorrila_test(0, 4, 2, 20, 5);
}
else {
gorrila_test(0, 10, 7, 20, 11);
}
return;
std::cout << "ordered support\n";
g_use_ordered_support = true;
tst4();
std::cout << "non-ordered support\n";
g_use_ordered_support = false;
tst4();
}