progress in gomory cut

Signed-off-by: Lev Nachmanson <levnach@microsoft.com>

src/smt/theory_lra.cpp

@@ -105,6 +105,7 @@ namespace lp_api {
         unsigned m_bound_propagations1;
         unsigned m_bound_propagations2;
         unsigned m_assert_diseq;
+        unsigned m_gomory_cuts;
         stats() { reset(); }
         void reset() {
             memset(this, 0, sizeof(*this));
@@ -1253,6 +1254,7 @@ namespace smt {
                 return l_false;
             }
             case lp::lia_move::cut: {
+                ++m_stats.m_gomory_cuts;
                 // m_explanation implies term <= k
                 app_ref b = mk_bound(term, k);
                 m_eqs.reset();

src/util/lp/int_solver.cpp

@@ -5,6 +5,7 @@
 
 #include "util/lp/int_solver.h"
 #include "util/lp/lar_solver.h"
+#include "util/lp/antecedents.h"
 namespace lp {
 
 void int_solver::fix_non_base_columns() {
@@ -128,108 +129,107 @@ bool int_solver::is_gomory_cut_target() {
     return true;
 }
 
+
+void int_solver::is_real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, buffer<row_entry> & pol) {
+    mpq f_0  = fractional_part(get_value(m_gomory_cut_inf_column));
+    mpq new_a;
+    if (at_lower(x_j)) {
+        if (a.is_pos()) {
+            new_a  =  a / (mpq(1) - f_0);
+        }
+        else {
+            new_a  =  a / f_0;
+            new_a.neg();
+        }
+        k += lower_bound(x_j).x * k;  // k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation
+        //        lower(x_j)->push_justification(ante, new_a, coeffs_enabled());*/
+    }
+    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 / (mpq(1) - f_0); 
+        }
+        k += upper_bound(x_j).x * k; // k.addmul(new_a, upper_bound(x_j).get_rational());
+        //        upper(x_j)->push_justification(ante, new_a, coeffs_enabled());*/
+    }
+    TRACE("gomory_cut_detail", tout << a << "*v" << x_j << " k: " << k << "\n";);
+    pol.push_back(row_entry(new_a, x_j));
+}
+void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, buffer<row_entry> & pol) {
+    /*
+    ++num_ints;
+    SASSERT(is_int(x_j));
+    mpq f_j = Ext::fractional_part(a);
+    TRACE("gomory_cut_detail", 
+          tout << a << "*v" << x_j << "\n";
+          tout << "fractional_part: " << Ext::fractional_part(a) << "\n";
+          tout << "f_j: " << f_j << "\n";
+          tout << "f_0: " << f_0 << "\n";
+          tout << "one_minus_f_0: " << one_minus_f_0 << "\n";);
+    if (!f_j.is_zero()) {
+        mpq new_a;
+        if (at_lower(x_j)) {
+            if (f_j <= one_minus_f_0) {
+                new_a = f_j / one_minus_f_0;
+            }
+            else {
+                new_a = (mpq(1) - f_j) / f_0;
+            }
+            k.addmul(new_a, lower_bound(x_j).get_rational());
+            lower(x_j)->push_justification(ante, new_a, coeffs_enabled());
+        }
+        else {
+            SASSERT(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
+            k.addmul(new_a, upper_bound(x_j).get_rational());
+            upper(x_j)->push_justification(ante, new_a, coeffs_enabled());
+        }
+        TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << k << "\n";);
+        pol.push_back(row_entry(new_a, x_j));
+        lcm_den = lcm(lcm_den, denominator(new_a));
+        }*/
+}
+
 lia_move int_solver::mk_gomory_cut(explanation & ex) {
 
     lp_assert(column_is_int_inf(m_gomory_cut_inf_column));
 
-    return lia_move::give_up;
+    TRACE("gomory_cut", tout << "applying cut at:\n"; m_lar_solver->print_linear_iterator(m_iter_on_gomory_row, tout); tout << "\n";);
-    
-    /*
-    TRACE("gomory_cut", tout << "applying cut at:\n"; display_row_info(tout, r););
         
-    antecedents ante(*this);
+    antecedents ante();
-
-    m_stats.m_gomory_cuts++;
-
     // gomory will be   pol >= k
-    numeral k(1);
+    mpq k(1);
     buffer<row_entry> pol;
-        
+       
-    numeral f_0  = Ext::fractional_part(m_value[x_i]);
+    mpq f_0  = fractional_part(get_value(m_gomory_cut_inf_column));
-    numeral one_minus_f_0 = numeral(1) - f_0; 
+    mpq one_minus_f_0 = mpq(1) - f_0;
-    SASSERT(!f_0.is_zero());
+    lp_assert(!is_zero(f_0) && !is_zero(one_minus_f_0));
-    SASSERT(!one_minus_f_0.is_zero());
+    mpq lcm_den(1);
-        
-    numeral lcm_den(1);
     unsigned num_ints = 0;
+    unsigned x_j;
+    mpq a;
 
-    typename vector<row_entry>::const_iterator it  = r.begin_entries();
+    while (m_iter_on_gomory_row->next(a, x_j)) {
-    typename vector<row_entry>::const_iterator end = r.end_entries();
+        if (x_j == m_gomory_cut_inf_column)
-    for (; it != end; ++it) {
+            continue;
-        if (!it->is_dead() && it->m_var != x_i) {
+        // make the format compatible with the format used in: Integrating Simplex with DPLL(T)
-            theory_var x_j   = it->m_var;
+        a.neg();  
-            numeral a_ij = it->m_coeff;
+        if (is_real(x_j)) 
-            a_ij.neg();  // make the used format compatible with the format used in: Integrating Simplex with DPLL(T)
+            is_real_case_in_gomory_cut(a, x_j, k, pol);
-            if (is_real(x_j)) {
+        else
-                numeral new_a_ij;
+            int_case_in_gomory_cut(a, x_j, k, pol);
-                if (at_lower(x_j)) {
-                    if (a_ij.is_pos()) {
-                        new_a_ij  =  a_ij / one_minus_f_0;
-                    }
-                    else {
-                        new_a_ij  =  a_ij / f_0;
-                        new_a_ij.neg();
-                    }
-                    k.addmul(new_a_ij, lower_bound(x_j).get_rational());
-                    lower(x_j)->push_justification(ante, new_a_ij, coeffs_enabled());
-                }
-                else {
-                    SASSERT(at_upper(x_j));
-                    if (a_ij.is_pos()) {
-                        new_a_ij =   a_ij / f_0; 
-                        new_a_ij.neg(); // the upper terms are inverted.
-                    }
-                    else {
-                        new_a_ij =   a_ij / one_minus_f_0; 
-                    }
-                    k.addmul(new_a_ij, upper_bound(x_j).get_rational());
-                    upper(x_j)->push_justification(ante, new_a_ij, coeffs_enabled());
-                }
-                TRACE("gomory_cut_detail", tout << a_ij << "*v" << x_j << " k: " << k << "\n";);
-                pol.push_back(row_entry(new_a_ij, x_j));
-            }
-            else {
-                ++num_ints;
-                SASSERT(is_int(x_j));
-                numeral f_j = Ext::fractional_part(a_ij);
-                TRACE("gomory_cut_detail", 
-                      tout << a_ij << "*v" << x_j << "\n";
-                      tout << "fractional_part: " << Ext::fractional_part(a_ij) << "\n";
-                      tout << "f_j: " << f_j << "\n";
-                      tout << "f_0: " << f_0 << "\n";
-                      tout << "one_minus_f_0: " << one_minus_f_0 << "\n";);
-                if (!f_j.is_zero()) {
-                    numeral new_a_ij;
-                    if (at_lower(x_j)) {
-                        if (f_j <= one_minus_f_0) {
-                            new_a_ij = f_j / one_minus_f_0;
-                        }
-                        else {
-                            new_a_ij = (numeral(1) - f_j) / f_0;
-                        }
-                        k.addmul(new_a_ij, lower_bound(x_j).get_rational());
-                        lower(x_j)->push_justification(ante, new_a_ij, coeffs_enabled());
-                    }
-                    else {
-                        SASSERT(at_upper(x_j));
-                        if (f_j <= f_0) {
-                            new_a_ij = f_j / f_0;
-                        }
-                        else {
-                            new_a_ij = (numeral(1) - f_j) / one_minus_f_0;
-                        }
-                        new_a_ij.neg(); // the upper terms are inverted
-                        k.addmul(new_a_ij, upper_bound(x_j).get_rational());
-                        upper(x_j)->push_justification(ante, new_a_ij, coeffs_enabled());
-                    }
-                    TRACE("gomory_cut_detail", tout << "new_a_ij: " << new_a_ij << " k: " << k << "\n";);
-                    pol.push_back(row_entry(new_a_ij, x_j));
-                    lcm_den = lcm(lcm_den, denominator(new_a_ij));
-                }
-            }
-        }
     }
-
+    /*
     CTRACE("empty_pol", pol.empty(), display_row_info(tout, r););
 
     expr_ref bound(get_manager());
@@ -276,8 +276,8 @@ lia_move int_solver::mk_gomory_cut(explanation & ex) {
                   }
                   tout << "k: " << k << "\n";);
         }
-        mk_polynomial_ge(pol.size(), pol.c_ptr(), k.to_rational(), bound);            
+        mk_polynomial_ge(pol.size(), pol.c_ptr(), k.to_rational(), bound);            */
-    }
+    /*
     TRACE("gomory_cut", tout << "new cut:\n" << bound << "\n"; ante.display(tout););
     literal l     = null_literal;
     context & ctx = get_context();
@@ -292,6 +292,7 @@ lia_move int_solver::mk_gomory_cut(explanation & ex) {
                                                                 ante.eqs().size(), ante.eqs().c_ptr(), ante, l)));
     return true;
     */
+   return lia_move::give_up;
 }
 
 void int_solver::init_check_data() {
@@ -327,9 +328,10 @@ lia_move int_solver::proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& e
         return lia_move::continue_with_check;
     }
 
+    lia_move ret = mk_gomory_cut(ex);
     delete m_iter_on_gomory_row;
     m_iter_on_gomory_row = nullptr;
-    return mk_gomory_cut(ex);
+    return ret;
  }
 
 
@@ -792,6 +794,10 @@ bool int_solver::is_int(unsigned j) const {
     return m_lar_solver->column_is_int(j);
 }
 
+bool int_solver::is_real(unsigned j) const {
+    return !is_int(j);
+}
+
 bool int_solver::value_is_int(unsigned j) const {
     return m_lar_solver->m_mpq_lar_core_solver.m_r_x[j].is_int();
 }
@@ -856,6 +862,7 @@ bool int_solver::is_free(unsigned j) const {
 bool int_solver::at_bound(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:
     case column_type::boxed:
         return
             mpq_solver.m_low_bounds[j] == get_value(j) ||
@@ -869,6 +876,30 @@ bool int_solver::at_bound(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:
+    case column_type::boxed:
+    case column_type::low_bound:
+        return mpq_solver.m_low_bounds[j] == get_value(j);
+    default:
+        return false;
+    }
+}
+
+bool int_solver::at_upper(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:
+    case column_type::boxed:
+    case column_type::upper_bound:
+        return mpq_solver.m_upper_bounds[j] == get_value(j);
+    default:
+        return false;
+    }
+}
+
 
 
 lp_settings& int_solver::settings() {

src/util/lp/int_solver.h

@@ -27,6 +27,11 @@ struct explanation {
 };
 
 class int_solver {
+    struct row_entry {
+        mpq m_coeff;
+        unsigned m_var;
+        row_entry(const mpq & coeff, unsigned var) : m_coeff(coeff), m_var(var) {}
+    };
 public:
     // fields
     lar_solver *m_lar_solver;
@@ -82,6 +87,7 @@ private:
     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;
     bool is_base(unsigned j) const;
     bool is_boxed(unsigned j) const;
     bool is_free(unsigned j) const;
@@ -109,5 +115,19 @@ private:
     int find_next_free_var_in_gomory_row();
     bool is_gomory_cut_target();
     bool at_bound(unsigned j) const;
+    bool at_lower(unsigned j) const;
+    bool at_upper(unsigned j) const;
+
+    inline static bool is_rational(const impq & n) {
+        return is_zero(n.y);  
+    }
+
+    inline static
+    mpq fractional_part(const impq & n) {
+        lp_assert(is_rational);
+        return n.x - floor(n.x);
+    }
+    void is_real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, buffer<row_entry> & pol);
+    void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, buffer<row_entry> & pol);
 };
 }

src/util/lp/lar_solver.h

@@ -414,9 +414,7 @@ public:
         return v.is_int();
     }
 
-    
+    bool column_is_real(unsigned j) const {
-    
-     bool column_is_real(unsigned j) const {
         return !column_is_int(j);
     }