gomory cut worked on a toy example

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

src/util/lp/int_solver.cpp

@@ -213,7 +213,82 @@ void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, la
     }
 }
 
-lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl ) {
+lia_move int_solver::report_conflict_from_gomory_cut(mpq & k) {
+    TRACE("empty_pol",
+          display_row_info(tout,
+                           m_lar_solver->m_mpq_lar_core_solver.m_r_heading[m_gomory_cut_inf_column]););
+    lp_assert(k.is_pos());
+    // conflict 0 >= k where k is positive
+    k.neg(); // returning 0 <= -k
+    return lia_move::conflict;
+}
+
+void int_solver::gomory_cut_adjust_t_and_k_for_size_gt_1(
+                                                         vector<std::pair<mpq, unsigned>> & pol,
+                                                         lar_term & t,
+                                                         mpq &k,
+                                                         unsigned num_ints,
+                                                         mpq & lcm_den) {
+        if (num_ints > 0) {
+            lcm_den = lcm(lcm_den, denominator(k));
+            TRACE("gomory_cut_detail", tout << "k: " << k << " lcm_den: " << lcm_den << "\n";
+                  linear_combination_iterator_on_vector<mpq> pi(pol);
+                  m_lar_solver->print_linear_iterator(&pi, tout);
+                  tout << "\nk: " << k << "\n";);
+            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;
+            }
+            TRACE("gomory_cut_detail", tout << "after *lcm_den\n";
+                  for (unsigned i = 0; i < pol.size(); i++) {
+                      tout << pol[i].first << " * v" << pol[i].second << "\n";
+                  }
+                  tout << "k: " << k << "\n";);
+        }
+        t.clear();
+        // negate everything to return -pol <= -k
+        for (const auto & pi: pol)
+            t.add_monoid(-pi.first, pi.second);
+        k.neg();
+}
+
+
+void int_solver::gomory_cut_adjust_t_and_k_for_size_1(const vector<std::pair<mpq, unsigned>> & pol, lar_term& t, mpq &k) {
+    lp_assert(pol.size() == 1);
+    unsigned j = pol[0].second;
+    k /= pol[0].first;
+    bool is_lower = pol[0].first.is_pos();
+    if (is_int(j) && !k.is_int()) {
+        k = is_lower?ceil(k):floor(k);
+    }
+    if (is_lower) { // returning -t <= -k which is equivalent to t >= k
+        k.neg();
+        t.negate();
+    }
+}
+
+
+
+lia_move int_solver::report_gomory_cut(lar_term& t, mpq& k, mpq &lcm_den, unsigned num_ints) {
+    lp_assert(!t.is_empty());
+    auto pol = t.coeffs_as_vector();
+    if (pol.size() == 1)
+        gomory_cut_adjust_t_and_k_for_size_1(pol, t, k);
+    else
+        gomory_cut_adjust_t_and_k_for_size_gt_1(pol, t, k, num_ints, lcm_den);
+    m_lar_solver->subs_terms_for_debugging(t); // todo: remove later
+    return lia_move::cut;
+}
+
+
+
+
+lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl) {
 
     lp_assert(column_is_int_inf(m_gomory_cut_inf_column));
 
@@ -238,59 +313,11 @@ lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl ) {
         }
     }
 
-    if (t.is_empty()) {
+    if (t.is_empty())
-        TRACE("empty_pol",
+        return report_conflict_from_gomory_cut(k);
-              display_row_info(tout,
-                               m_lar_solver->m_mpq_lar_core_solver.m_r_heading[m_gomory_cut_inf_column]););
-        lp_assert(k.is_pos());
-        // conflict 0 >= k where k is positive
-        k.neg(); // returning 0 <= -k
-        return lia_move::conflict;
-    }
 
-    auto pol = t.coeffs_as_vector();
+    return report_gomory_cut(t, k, lcm_den, num_ints);
-    if (pol.size() == 1) {
+    
-        unsigned j = pol[0].second;
-        k /= pol[0].first;
-        bool is_lower = pol[0].first.is_pos();
-        if (is_int(j) && !k.is_int()) {
-            k = is_lower?ceil(k):floor(k);
-        }
-        if (is_lower) { // returning -t <= -k which is equivalent to t >= k
-            k.neg();
-            t.negate();
-        }
-    } else { 
-        if (num_ints > 0) {
-            lcm_den = lcm(lcm_den, denominator(k));
-            TRACE("gomory_cut_detail", tout << "k: " << k << " lcm_den: " << lcm_den << "\n";
-                  linear_combination_iterator_on_vector<mpq> pi(pol);
-                  m_lar_solver->print_linear_iterator(&pi, tout);
-                  tout << "\nk: " << k << "\n";);
-            SASSERT(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;
-            }
-            TRACE("gomory_cut_detail", tout << "after *lcm\n";
-                  for (unsigned i = 0; i < pol.size(); i++) {
-                      tout << pol[i].first << " * v" << pol[i].second << "\n";
-                  }
-                  tout << "k: " << k << "\n";);
-            t.clear();
-            // negate everything to return -pol <= -k
-            for (const auto & pi: pol)
-                t.add_monoid(-pi.first, pi.second);
-            k.neg();
-        } else {
-            lp_assert(false); // not sure what happens here
-        }
-    }
-    return lia_move::cut;
 }
 
 void int_solver::init_check_data() {
@@ -370,6 +397,7 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex) {
         if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
             return lia_move::give_up;
         }
+        init_inf_int_set(); // todo - can we avoid this call?
         int j = find_inf_int_base_column();
         if (j != -1) {
             // setup the call for gomory cut

src/util/lp/int_solver.h

@@ -106,6 +106,8 @@ private:
     void move_non_base_vars_to_bounds();
     void branch_infeasible_int_var(unsigned);
     lia_move mk_gomory_cut(lar_term& t, mpq& k,explanation & ex);
+    lia_move report_conflict_from_gomory_cut(mpq & k);
+    lia_move report_gomory_cut(lar_term& t, mpq& k, mpq& lcm_den, unsigned num_ints);
 	void init_check_data();
     bool constrain_free_vars(linear_combination_iterator<mpq> *  r);
     lia_move proceed_with_gomory_cut(lar_term& t, mpq& k, explanation& ex);
@@ -129,5 +131,7 @@ private:
     constraint_index column_upper_bound_constraint(unsigned j) const;
     constraint_index column_low_bound_constraint(unsigned j) const;
     void display_row_info(std::ostream & out, unsigned row_index) const;
+    void gomory_cut_adjust_t_and_k_for_size_1(const vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k);
+    void gomory_cut_adjust_t_and_k_for_size_gt_1(vector<std::pair<mpq, unsigned>> & pol, lar_term & t, mpq &k, unsigned num_ints, mpq &lcm_den);
 };
 }

src/util/lp/lar_solver.cpp

@@ -542,7 +542,7 @@ void lar_solver::set_low_bound_witness(var_index j, constraint_index ci) {
     m_columns_to_ul_pairs[j] = ul;
 }
 
-void lar_solver::register_one_coeff_in_map(std::unordered_map<var_index, mpq> & coeffs, const mpq & a, unsigned j) {
+void lar_solver::register_monoid_in_map(std::unordered_map<var_index, mpq> & coeffs, const mpq & a, unsigned j) {
     auto it = coeffs.find(j);
     if (it == coeffs.end()) {
         coeffs[j] = a;
@@ -553,18 +553,18 @@ void lar_solver::register_one_coeff_in_map(std::unordered_map<var_index, mpq> &
 
 
 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 & right_side) const {
+                                                       vector<std::pair<mpq, var_index>> &left_side, mpq & free_coeff) const {
     std::unordered_map<var_index, mpq> coeffs;
     for (auto & t : left_side_with_terms) {
         unsigned j = t.second;
         if (!is_term(j)) {
-            register_one_coeff_in_map(coeffs, t.first, j);
+            register_monoid_in_map(coeffs, t.first, j);
         } else {
             const lar_term & term = * m_terms[adjust_term_index(t.second)];
             for (auto & p : term.coeffs()){
-                register_one_coeff_in_map(coeffs, t.first * p.second , p.first);
+                register_monoid_in_map(coeffs, t.first * p.second , p.first);
             }
-            right_side += t.first * term.m_v;
+            free_coeff += t.first * term.m_v;
         }
     }
 

src/util/lp/lar_solver.h

@@ -149,8 +149,6 @@ public:
     lp_settings const & settings() const;
 
     void clear();
-
-
     lar_solver();
     void set_propagate_bounds_on_pivoted_rows_mode(bool v);
 
@@ -263,7 +261,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 & right_side) const;
+                          vector<std::pair<mpq, var_index>> &left_side, mpq & free_coeff) const;
 
 
     void detect_rows_of_bound_change_column_for_nbasic_column(unsigned j);
@@ -343,7 +341,7 @@ public:
     bool the_relations_are_of_same_type(const vector<std::pair<mpq, unsigned>> & evidence, lconstraint_kind & the_kind_of_sum) const;
 
     static void register_in_map(std::unordered_map<var_index, mpq> & coeffs, const lar_base_constraint & cn, const mpq & a);
-    static void register_one_coeff_in_map(std::unordered_map<var_index, mpq> & coeffs, const mpq & a, unsigned j);
+    static void register_monoid_in_map(std::unordered_map<var_index, mpq> & coeffs, const mpq & a, unsigned j);
 
 
     bool the_left_sides_sum_to_zero(const vector<std::pair<mpq, unsigned>> & evidence) const;
@@ -463,5 +461,17 @@ public:
         return m_columns_to_ul_pairs()[j].low_bound_witness();
     }
 
+    void subs_terms_for_debugging(lar_term& t) {
+        vector<std::pair<mpq, unsigned>> pol;
+        for (const auto & m : t.m_coeffs) {
+            pol.push_back(std::make_pair(m.second, adjust_column_index_to_term_index(m.first)));
+        }
+        mpq v = t.m_v;
+        
+        vector<std::pair<mpq, unsigned>> pol_after_subs;
+        substitute_terms_in_linear_expression(pol, pol_after_subs, v);
+        t.clear();
+        t = lar_term(pol_after_subs, v);
+    }
 };
 }