Merge pull request #1834 from levnach/gomory

Gomory

src/smt/smt_arith_value.h

@@ -17,7 +17,7 @@ Author:
 Revision History:
 
 --*/
-#pragma once;
+#pragma once
 
 #include "ast/arith_decl_plugin.h"
 #include "smt/smt_context.h"

src/test/lp/gomory_test.h

@@ -1,4 +1,5 @@
 namespace lp {
+#include "util/lp/lp_utils.h"
 struct gomory_test {
     gomory_test(
         std::function<std::string (unsigned)> name_function_p,
@@ -88,7 +89,7 @@ struct gomory_test {
         lp_assert(is_int(x_j));
         lp_assert(!a.is_int());
              lp_assert(f_0 > zero_of_type<mpq>() && f_0 < one_of_type<mpq>());
-        mpq f_j =  int_solver::fractional_part(a);
+        mpq f_j =  fractional_part(a);
         TRACE("gomory_cut_detail", 
               tout << a << " x_j = " << x_j << ", k = " << k << "\n";
               tout << "f_j: " << f_j << "\n";
@@ -206,7 +207,7 @@ struct gomory_test {
         unsigned x_j;
         mpq a;
         bool some_int_columns = false;
-        mpq f_0  = int_solver::fractional_part(get_value(inf_col));
+        mpq f_0  = fractional_part(get_value(inf_col));
         mpq one_min_f_0 = 1 - f_0;
         for ( auto pp : row) {
             a = pp.first;

src/util/lp/column_namer.h

@@ -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>

src/util/lp/gomory.cpp

@@ -20,6 +20,7 @@
 #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 {
@@ -43,32 +44,31 @@ class gomory::imp {
     void int_case_in_gomory_cut(const mpq & a, unsigned j,
                                             mpq & lcm_den, const mpq& f0, const mpq& one_minus_f0) {
         lp_assert(is_int(j) && !a.is_int());
-        mpq fj =  int_solver::fractional_part(a);
-        lp_assert(fj.is_pos());
+        mpq fj =  fractional_part(a);
         TRACE("gomory_cut_detail", 
               tout << a << " j=" << j << " k = " << m_k;
               tout << ", fj: " << fj << ", ";
-              tout << "f0: " << f0 << ", ";
-              tout << "1 - f0: " << 1 - f0 << ", ";
+              tout << "a - fj = " << a - fj <<  ", ";
               tout << (at_lower(j)?"at_lower":"at_upper")<< std::endl;
               );
+        lp_assert(fj.is_pos() && (a - fj).is_int());
         mpq new_a;
         mpq one_minus_fj = 1 - fj;
         if (at_lower(j)) {
-            new_a = fj < one_minus_f0? fj / one_minus_f0 : one_minus_fj / f0;
+            new_a = fj < one_minus_f0? fj / one_minus_f0 : (- one_minus_fj / f0);
             m_k.addmul(new_a, lower_bound(j).x);
             m_ex.push_justification(column_lower_bound_constraint(j), new_a);
         }
         else {
             lp_assert(at_upper(j));
             // the upper terms are inverted: therefore we have the minus
-            new_a = - (fj < f0? fj / f0 : one_minus_fj / one_minus_f0);
+            new_a = fj < f0? (- fj / f0 ) : (one_minus_fj / one_minus_f0);
             m_k.addmul(new_a, upper_bound(j).x);
             m_ex.push_justification(column_upper_bound_constraint(j), new_a);
         }
-        TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << m_k << "\n";);
         m_t.add_monomial(new_a, j);
         lcm_den = lcm(lcm_den, denominator(new_a));
+        TRACE("gomory_cut_detail", tout << "new_a = " << new_a << ", k = " << m_k << ", lcm_den = " << lcm_den << "\n";);
     }
 
     void real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f0, const mpq& one_minus_f0) {
@@ -150,10 +150,120 @@ class gomory::imp {
         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);
+    }
+
+    void dump_coeff_val(std::ostream & out, const mpq & a) const {
+        if (a.is_int()) {
+            if ( a >= zero_of_type<mpq>())
+                out << a;
+            else {
+                out << "( - " <<  - a << ") ";
+            }
+        } else {
+            if ( a >= zero_of_type<mpq>())
+                out << "(div " << numerator(a) << " " << denominator(a) << ")";
+            else {
+                out << "(- ( div " <<   numerator(-a) << " " << denominator(-a) << "))";
+            }
+            
+        }
+    }
+    
+    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()) << ")";
+    }
+    
+    void dump_row_coefficients(std::ostream & out) const {
+        for (const auto& p : m_row) {
+            if (!column_is_fixed(p.var()))
+                dump_coeff(out, p);
+        }
+    }
+    
+    void dump_the_row(std::ostream& out) const {
+        out << "; the row, excluding fixed vars\n";
+        out << "(assert ( = ( + ";
+        dump_row_coefficients(out);
+        out << ") 0))\n";
+    }
+    
+    void dump_declarations(std::ostream& out) const {
+        // for a column j the var name is vj
+        for (const auto & p : m_row) {
+            if (column_is_fixed(p.var())) continue;
+            out << "(declare-fun " << var_name(p.var()) << " () "
+                << (is_int(p.var())? "Int" : "Real") << ")\n";
+        }
+    }
+
+    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 (column_is_fixed(j)) continue;
+            if (j == m_inf_col) {
+                continue;
+            }
+
+            if (column_is_fixed(j)) {
+                dump_lower_bound_expl(out, j);
+                dump_upper_bound_expl(out, j);
+                continue;
+            }
+
+            if (at_lower(j)) {
+                dump_lower_bound_expl(out, j);
+            } else {
+                dump_upper_bound_expl(out, j);
+            }
+        }
+    }
+
+    void dump_terms_coefficients(std::ostream & out) const {
+        for (const auto& p : m_t) {
+            dump_coeff(out, p);
+        }
+    }
+    
+    void dump_term_sum(std::ostream & out) const {
+        out << "( + ";
+        dump_terms_coefficients(out);
+        out << ")";
+    }
+    
+    void dump_term_le_k(std::ostream & out) const {
+        out << "( <= ";
+        dump_term_sum(out);
+        out << m_k << ")";
+    }
+    void dump_the_cut_assert(std::ostream & out) const {
+        out <<"(assert (not ";
+        dump_term_le_k(out);
+        out << "))\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"; m_int_solver.m_lar_solver->print_row(m_row, tout); tout << std::endl;
+              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);
               }
@@ -164,7 +274,11 @@ public:
         m_k = 1;
         mpq lcm_den(1);
         bool some_int_columns = false;
-        mpq f0  = int_solver::fractional_part(get_value(m_inf_col));
+        mpq f0  = fractional_part(get_value(m_inf_col));
+        TRACE("gomory_cut_detail", tout << "f0: " << f0 << ", ";
+              tout << "1 - f0: " << 1 - f0 << ", get_value(m_inf_col).x - f0 = " << get_value(m_inf_col).x - f0;);
+        lp_assert(f0.is_pos() && (get_value(m_inf_col).x - f0).is_int());  
+
         mpq one_min_f0 = 1 - f0;
         for (const auto & p : m_row) {
             unsigned j = p.var();
@@ -194,7 +308,8 @@ public:
             adjust_term_and_k_for_some_ints_case_gomory(lcm_den);
         lp_assert(m_int_solver.current_solution_is_inf_on_cut());
         m_int_solver.m_lar_solver->subs_term_columns(m_t, m_k);
-        TRACE("gomory_cut", tout<<"gomory cut:"; m_int_solver.m_lar_solver->print_term(m_t, tout); tout << " <= " << m_k << std::endl;);
+        TRACE("gomory_cut", tout<<"gomory cut:"; print_linear_combination_of_column_indices_only(m_t, tout); tout << " <= " << m_k << std::endl;);
+        TRACE("gomory_cut_detail", dump_cut_and_constraints_as_smt_lemma(tout););
         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 ) :

src/util/lp/int_solver.h

@@ -113,17 +113,9 @@ private:
     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);
-    }
     constraint_index column_upper_bound_constraint(unsigned j) const;
     constraint_index column_lower_bound_constraint(unsigned j) const;
     bool current_solution_is_inf_on_cut() const;

src/util/lp/lar_solver.cpp

@@ -1273,7 +1273,7 @@ std::ostream& lar_solver::print_term_as_indices(lar_term const& term, std::ostre
     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);
+    print_linear_combination_of_column_indices_only(term, out);
     return out;
 }
 

src/util/lp/lp_settings.h

@@ -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();

src/util/lp/lp_utils.h

@@ -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); }