hilbert validation

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

src/muz_qe/hilbert_basis.cpp

@@ -410,6 +410,7 @@ void hilbert_basis::collect_statistics(statistics& st) const {
     st.update("hb.num_subsumptions", m_stats.m_num_subsumptions);
     st.update("hb.num_resolves", m_stats.m_num_resolves);
     st.update("hb.num_saturations", m_stats.m_num_saturations);
+    st.update("hb.basis_size", get_basis_size());
     m_index->collect_statistics(st);
 }
 
@@ -464,6 +465,10 @@ void hilbert_basis::set_is_int(unsigned var_index) {
     m_ints.push_back(var_index+1);
 }
 
+bool hilbert_basis::get_is_int(unsigned var_index) const {    
+    return m_ints.contains(var_index+1);
+}
+
 unsigned hilbert_basis::get_num_vars() const {
     if (m_ineqs.empty()) {
         return 0;
@@ -579,6 +584,23 @@ lbool hilbert_basis::saturate(num_vector const& ineq, bool is_eq) {
     return l_true;
 }
 
+void hilbert_basis::get_basis_solution(unsigned i, num_vector& v, bool& is_initial) {
+    offset_t offs = m_basis[i];
+    v.reset();
+    for (unsigned i = 1; i < get_num_vars(); ++i) {
+        v.push_back(vec(offs)[i]);
+    }
+    is_initial = !vec(offs)[0].is_zero();
+}
+
+void hilbert_basis::get_ge(unsigned i, num_vector& v, numeral& b, bool& is_eq) {
+    v.reset();
+    v.append(get_num_vars()-1, m_ineqs[i].c_ptr() + 1);
+    b = -m_ineqs[i][0];
+    is_eq = m_iseq[i];
+}
+
+
 void hilbert_basis::select_inequality() {
     SASSERT(m_current_ineq < m_ineqs.size());
     unsigned best = m_current_ineq;
@@ -806,19 +828,6 @@ void hilbert_basis::display_ineq(std::ostream& out, num_vector const& v, bool is
 }
 
 
-void hilbert_isl_basis::add_le(num_vector const& v, numeral bound) {
-    unsigned sz = v.size();
-    num_vector w;
-    w.push_back(-bound);
-    w.push_back(bound);
-    for (unsigned i = 0; i < sz; ++i) {
-        w.push_back(v[i]);
-        w.push_back(-v[i]);
-    }
-    m_basis.add_le(w);
-}
-
-
 /**
    Vector v is subsumed by vector w if
        

src/muz_qe/hilbert_basis.h

@@ -11,8 +11,6 @@ Abstract:
 
     hilbert_basis computes a Hilbert basis for linear
     homogeneous inequalities over naturals.
-    hilbert_sl_basis  computes a semi-linear set over naturals.
-    hilbert_isl_basis computes semi-linear sets over integers.
 
 Author:
 
@@ -140,33 +138,24 @@ public:
     void add_eq(num_vector const& v, numeral const& b);
 
     void set_is_int(unsigned var_index);
+    bool get_is_int(unsigned var_index) const;
 
     lbool saturate();
 
+    unsigned get_basis_size() const { return m_basis.size(); }
+    void get_basis_solution(unsigned i, num_vector& v, bool& is_initial);
+
+    unsigned get_num_ineqs() const { return m_ineqs.size(); }
+    void get_ge(unsigned i, num_vector& v, numeral& b, bool& is_eq);    
+
     void set_cancel(bool f) { m_cancel = f; }
 
     void display(std::ostream& out) const;
 
     void collect_statistics(statistics& st) const;
     void reset_statistics();     
+
 };
 
 
-class hilbert_isl_basis {
-public:
-    typedef rational        numeral;
-    typedef vector<numeral> num_vector;
-private:
-    hilbert_basis m_basis;    
-public:
-    hilbert_isl_basis() {}    
-    void reset() { m_basis.reset(); }
-
-    // add inequality v*x >= bound, x ranges over integers
-    void add_le(num_vector const& v, numeral bound);
-    lbool saturate() { return m_basis.saturate(); }
-    void set_cancel(bool f) { m_basis.set_cancel(f); }
-    void display(std::ostream& out) const { m_basis.display(out); }    
-};
-
 #endif 

src/muz_qe/hilbert_basis_validate.cpp

@@ -0,0 +1,178 @@
+/*++
+Copyright (c) 2013 Microsoft Corporation
+
+Module Name:
+
+    hilbert_basis_validate.cpp
+
+Abstract:
+
+    Basic Hilbert Basis validation.
+
+    hilbert_basis computes a Hilbert basis for linear
+    homogeneous inequalities over naturals.
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2013-02-15.
+
+Revision History:
+
+--*/
+
+#include "hilbert_basis_validate.h"
+#include "arith_decl_plugin.h"
+#include "ast_pp.h"
+#include <sstream>
+
+
+hilbert_basis_validate::hilbert_basis_validate(ast_manager& m): 
+    m(m) {
+}
+
+void hilbert_basis_validate::validate_solution(hilbert_basis& hb, vector<rational> const& v, bool is_initial) {
+    unsigned sz = hb.get_num_ineqs();
+    rational bound;
+    for (unsigned i = 0; i < sz; ++i) {
+        bool is_eq;
+        vector<rational> w;
+        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 (bound > sum || 
+            (is_eq && bound != sum)) {
+            // validation failed.
+            std::cout << "validation failed for inequality\n";
+            for (unsigned j = 0; j < v.size(); ++j) {
+                std::cout << v[j] << " ";
+            }
+            std::cout << "\n";
+            for (unsigned j = 0; j < w.size(); ++j) {
+                std::cout << w[j] << " ";
+            }
+            std::cout << (is_eq?" = ":" >= ") << bound << "\n";
+            std::cout << "is initial: " << (is_initial?"true":"false") << "\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<rational> v;
+    bool is_initial;
+
+    // check that claimed solution really satisfies inequalities:        
+    for (unsigned i = 0; i < sz; ++i) {
+        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<expr_ref_vector> offsets, increments;
+    expr_ref_vector xs(m), vars(m);
+    expr_ref var(m);
+    svector<symbol> 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) {
+        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().c_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.c_ptr())));
+        }
+        fml = m.mk_and(eqs.size(), eqs.c_ptr());
+        if (!names.empty()) {
+            fml = m.mk_exists(names.size(), sorts.c_ptr(), names.c_ptr(), fml);
+        }
+        fmls.push_back(fml);
+    }
+    fml1 = m.mk_or(fmls.size(), fmls.c_ptr());
+    fmls.reset();
+    
+    sz = hb.get_num_ineqs();
+    for (unsigned i = 0; i < sz; ++i) {
+        bool is_eq;
+        vector<rational> 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.c_ptr());
+        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.c_ptr());
+    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.c_ptr()), fml);
+    }
+    return fml;
+
+}
+

src/muz_qe/hilbert_basis_validate.h

@@ -0,0 +1,43 @@
+/*++
+Copyright (c) 2013 Microsoft Corporation
+
+Module Name:
+
+    hilbert_basis_validate.h
+
+Abstract:
+
+    Basic Hilbert Basis validation.
+
+    hilbert_basis computes a Hilbert basis for linear
+    homogeneous inequalities over naturals.
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2013-02-15.
+
+Revision History:
+
+--*/
+
+#ifndef _HILBERT_BASIS_VALIDATE_H_
+#define _HILBERT_BASIS_VALIDATE_H_
+
+#include "hilbert_basis.h"
+#include "ast.h"
+
+class hilbert_basis_validate {
+    ast_manager& m;
+
+    void validate_solution(hilbert_basis& hb, vector<rational> const& v, bool is_initial);
+
+public:
+        
+    hilbert_basis_validate(ast_manager& m);
+
+    expr_ref mk_validate(hilbert_basis& hb);
+
+};
+
+
+#endif 

src/test/hilbert_basis.cpp

@@ -1,4 +1,12 @@
 #include "hilbert_basis.h"
+#include "hilbert_basis_validate.h"
+#include "ast_pp.h"
+#include "reg_decl_plugins.h"
+#include "quant_tactics.h"
+#include "tactic.h"
+#include "tactic2solver.h"
+#include "solver.h"
+
 #include<signal.h>
 #include<time.h>
 
@@ -19,6 +27,24 @@ static void on_ctrl_c(int) {
     raise(SIGINT);
 }
 
+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);
+
+    std::cout << mk_pp(fml, m) << "\n";
+
+    fml = m.mk_not(fml);
+    params_ref p;
+    tactic_ref tac = mk_lra_tactic(m, p);
+    ref<solver> sol = mk_tactic2solver(m, tac.get(), p);
+    sol->assert_expr(fml);
+    lbool r = sol->check_sat(0,0);
+    std::cout << r << "\n";
+}
+
 static void saturate_basis(hilbert_basis& hb) {
     signal(SIGINT, on_ctrl_c);
     g_hb = &hb;
@@ -29,6 +55,7 @@ static void saturate_basis(hilbert_basis& hb) {
     case l_true:  
         std::cout << "sat\n"; 
         hb.display(std::cout);
+        // validate_sat(hb);
         break;
     case l_false: 
         std::cout << "unsat\n"; 
@@ -40,6 +67,7 @@ static void saturate_basis(hilbert_basis& hb) {
     display_statistics(hb);
 }
 
+
 /**
    n         - number of variables.
    k         - subset of variables to be non-zero
@@ -258,16 +286,37 @@ static void tst12() {
     saturate_basis(hb);
 }
 
+// Sigma_9 table 1,  Ajili, Contejean
+static void tst13() {
+    hilbert_basis hb;
+    hb.add_eq(vec( 1,-2,-4,4), R(0));
+    hb.add_le(vec(100,45,-78,-67), R(0));
+    saturate_basis(hb);
+}
+
+// Sigma_10 table 1,  Ajili, Contejean
+static void tst14() {
+    hilbert_basis hb;
+    hb.add_le(vec( 23, -56, -34, 12, 11), R(0));
+    saturate_basis(hb);
+}
+
+// Sigma_11 table 1,  Ajili, Contejean
+static void tst15() {
+//    hilbert_basis hb;
+//    hb.add_le(vec( 23, -56, -34, 12, 11), R(0));
+//    saturate_basis(hb);
+}
+
+
 void tst_hilbert_basis() {
     std::cout << "hilbert basis test\n";
-    tst12();
-    return;
 
     if (true) {
         tst1();
         tst2();
         tst3();
-        tst4();
+        // tst4();
         tst5();
         tst6();
         tst7();
@@ -275,11 +324,15 @@ void tst_hilbert_basis() {
         tst9();
         tst10();
         tst11();
+        tst12();
+        tst13();
+        tst14();
+        tst15();
         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);
+        //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);