Merge pull request #2455 from levnach/fix

fix a bug in lar_solver in querying if a column is int

src/test/lp/gomory_test.h

@@ -86,7 +86,7 @@ struct gomory_test {
     }
     
     void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & t, explanation& expl, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) {
-        lp_assert(is_int(x_j));
+        lp_assert(is_integer(x_j));
         lp_assert(!a.is_int());
              lp_assert(f_0 > zero_of_type<mpq>() && f_0 < one_of_type<mpq>());
         mpq f_j =  fractional_part(a);
@@ -138,7 +138,7 @@ struct gomory_test {
         if (pol.size() == 1) {
             TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;);
             unsigned v = pol[0].second;
-            lp_assert(is_int(v));
+            lp_assert(is_integer(v));
             const mpq& a = pol[0].first;
             k /= a;
             if (a.is_pos()) { // we have av >= k
@@ -165,7 +165,7 @@ struct gomory_test {
                 // 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());
+                    SASSERT(!is_integer(pi.second) || pi.first.is_int());
                 }
                 k *= lcm_den;
             }

src/util/lp/gomory.cpp

@@ -47,7 +47,7 @@ class gomory::imp {
     bool column_is_fixed(unsigned j) const { return m_int_solver.m_lar_solver->column_is_fixed(j); }
 
     void int_case_in_gomory_cut(unsigned j) {
-        lp_assert(is_int(j) && m_fj.is_pos());
+        lp_assert(m_int_solver.column_is_int(j) && m_fj.is_pos());
         TRACE("gomory_cut_detail", 
               tout << " k = " << m_k;
               tout << ", fj: " << m_fj << ", ";
@@ -117,7 +117,7 @@ class gomory::imp {
         if (pol.size() == 1) {
             TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;);
             unsigned v = pol[0].second;
-            lp_assert(is_int(v));
+            lp_assert(m_int_solver.column_is_int(v));
             const mpq& a = pol[0].first;
             m_k /= a;
             if (a.is_pos()) { // we have av >= k
@@ -139,7 +139,7 @@ class gomory::imp {
                 // normalize coefficients of integer parameters to be integers.
                 for (auto & pi: pol) {
                     pi.first *= m_lcm_den;
-                    SASSERT(!is_int(pi.second) || pi.first.is_int());
+                    SASSERT(!m_int_solver.column_is_int(pi.second) || pi.first.is_int());
                 }
                 m_k *= m_lcm_den;
             }
@@ -195,7 +195,7 @@ class gomory::imp {
     }
 
     void dump_declaration(std::ostream& out, unsigned v) const {
-        out << "(declare-const " << var_name(v) << (is_int(v) ? " Int" : " Real") << ")\n";
+        out << "(declare-const " << var_name(v) << (m_int_solver.column_is_int(v) ? " Int" : " Real") << ")\n";
     }
     
     void dump_declarations(std::ostream& out) const {

src/util/lp/hnf.h

@@ -137,7 +137,7 @@ void pivot_column_non_fractional(M &m, unsigned r, bool & overflow, const mpq &
                 overflow = true;
                 return;
             }
-            lp_assert(is_int(m[i][j]));
+            lp_assert(is_integer(m[i][j]));
         }
     }
 }
@@ -577,7 +577,7 @@ private:
             process_row_modulo();
             lp_assert(is_pos(m_W[m_i][m_i]));
             m_R /= m_W[m_i][m_i];
-            lp_assert(is_int(m_R));
+            lp_assert(is_integer(m_R));
             m_half_R = floor(m_R / 2);
         }
     }

src/util/lp/hnf_cutter.h

@@ -120,7 +120,7 @@ public:
         int ret = -1;
         int n = 0;
         for (int i = 0; i < static_cast<int>(b.size()); i++) {
-            if (is_int(b[i])) continue;
+            if (is_integer(b[i])) continue;
             if (n == 0 ) {
                 lp_assert(ret == -1);
                 n = 1;

src/util/lp/int_solver.cpp

@@ -16,7 +16,7 @@ void int_solver::trace_inf_rows() const {
     TRACE("arith_int_rows",
           unsigned num = m_lar_solver->A_r().column_count();
           for (unsigned v = 0; v < num; v++) {
-              if (is_int(v) && !get_value(v).is_int()) {
+              if (column_is_int(v) && !get_value(v).is_int()) {
                   display_column(tout, v);
               }
           }
@@ -197,7 +197,7 @@ impq int_solver::get_cube_delta_for_term(const lar_term& t) const {
         bool seen_minus = false;
         bool seen_plus = false;
         for(const auto & p : t) {
-            if (!is_int(p.var()))
+            if (!column_is_int(p.var()))
                 goto usual_delta;
             const mpq & c = p.coeff();
             if (c == one_of_type<mpq>()) {
@@ -215,7 +215,7 @@ impq int_solver::get_cube_delta_for_term(const lar_term& t) const {
  usual_delta:
     mpq delta = zero_of_type<mpq>();
     for (const auto & p : t)
-        if (is_int(p.var()))
+        if (column_is_int(p.var()))
             delta += abs(p.coeff());
     
     delta *= mpq(1, 2);
@@ -759,7 +759,7 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq
         
         unsigned i = lcs.m_r_basis[row_index];
         impq const & xi = get_value(i);
-        if (is_int(i) && is_int(j) && !a.is_int())
+        if (column_is_int(i) && column_is_int(j) && !a.is_int())
             m = lcm(m, denominator(a));
         if (a.is_neg()) {
             if (has_low(i))
@@ -791,12 +791,12 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq
     return (inf_l || inf_u || l <= u);
 }
 
-bool int_solver::is_int(unsigned j) const {
+bool int_solver::column_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);
+    return !column_is_int(j);
 }
 
 bool int_solver::value_is_int(unsigned j) const {
@@ -821,7 +821,7 @@ void int_solver::display_column(std::ostream & out, unsigned j) const {
 }
 
 bool int_solver::column_is_int_inf(unsigned j) const {
-    return is_int(j) && (!value_is_int(j));
+    return column_is_int(j) && (!value_is_int(j));
 }
 
 bool int_solver::is_base(unsigned j) const {
@@ -912,7 +912,7 @@ bool int_solver::shift_var(unsigned j, unsigned range) {
         set_value_for_nbasic_column_ignore_old_values(j, new_val);
         return true;
     }
-    if (is_int(j)) {
+    if (column_is_int(j)) {
         if (!inf_l) {
             l = ceil(l);
             if (!m.is_one())
@@ -940,7 +940,7 @@ bool int_solver::shift_var(unsigned j, unsigned range) {
         set_value_for_nbasic_column_ignore_old_values(j, new_val);
         return true;
     }
-    if (!is_int(j)) {
+    if (!column_is_int(j)) {
         SASSERT(!inf_l && !inf_u);
         mpq delta       = mpq(random() % (range + 1));
         impq new_val = l + ((delta * (u - l)) / mpq(range)); 
@@ -975,7 +975,7 @@ bool int_solver::non_basic_columns_are_at_bounds() const {
                 return false;
             break;
         default:
-            if (is_int(j) && !val.is_int()) {
+            if (column_is_int(j) && !val.is_int()) {
                 return false;
             }
         }

src/util/lp/int_solver.h

@@ -59,7 +59,7 @@ public:
     bool is_real(unsigned j) const;
     const impq & lower_bound(unsigned j) const;
     const impq & upper_bound(unsigned j) const;
-    bool is_int(unsigned j) const;
+    bool column_is_int(unsigned j) const;
     const impq & get_value(unsigned j) const;
     bool at_lower(unsigned j) const;
     bool at_upper(unsigned j) const;

src/util/lp/lar_solver.cpp

@@ -512,7 +512,7 @@ bool lar_solver::move_non_basic_column_to_bounds(unsigned j) {
         }
         break;
     default:
-        if (is_int(j) && !val.is_int()) {
+        if (column_is_int(j) && !val.is_int()) {
             set_value_for_nbasic_column(j, impq(floor(val)));
             return true;
         }
@@ -526,6 +526,7 @@ void lar_solver::set_value_for_nbasic_column(unsigned j, const impq & new_val) {
     auto & x = m_mpq_lar_core_solver.m_r_x[j];
     auto delta = new_val - x;
     x = new_val;
+    m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j);
     change_basic_columns_dependend_on_a_given_nb_column(j, delta);
 }
 

src/util/lp/lp_core_solver_base.h

@@ -550,6 +550,7 @@ public:
     bool non_basic_columns_are_set_correctly() const {
         for (unsigned j : this->m_nbasis)
             if (!column_is_feasible(j)) {
+                TRACE("lp_core", tout << "inf col "; print_column_info(j, tout) << "\n";);
                 return false;
             }
         return true;
@@ -573,7 +574,7 @@ public:
         }
     }
 
-    void print_column_info(unsigned j, std::ostream & out) const {
+    std::ostream& print_column_info(unsigned j, std::ostream & out) const {
         out << "j = " << j << ",\tname = "<< column_name(j) << "\t";
         switch (m_column_types[j]) {
         case column_type::fixed:
@@ -598,6 +599,7 @@ public:
             out << " base\n";
         else
            out << " \n";
+        return out;
     }
 
     bool column_is_free(unsigned j) const { return this->m_column_type[j] == free; }

src/util/lp/lp_core_solver_base_def.h

@@ -511,12 +511,11 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::calc_current_x
 }
 
 template <typename T, typename X> bool lp_core_solver_base<T, X>::inf_set_is_correct() const {
-    unsigned j = this->m_n();
+    for (unsigned j = 0; j < this->m_n(); j++) {
-    while (j--) {
         bool belongs_to_set = m_inf_set.contains(j);
         bool is_feas = column_is_feasible(j);
-    
         if (is_feas == belongs_to_set) {
+            TRACE("lp_core", tout << "incorrectly set column in inf set "; print_column_info(j, tout) << "\n";);
             return false;
         }
     }

src/util/lp/lp_utils.h

@@ -86,7 +86,7 @@ template <typename X> inline X one_of_type() { return numeric_traits<X>::one();
 template <typename X> inline bool is_zero(const X & v) { return numeric_traits<X>::is_zero(v); }
 template <typename X> inline bool is_pos(const X & v) { return numeric_traits<X>::is_pos(v); }
 template <typename X> inline bool is_neg(const X & v) { return numeric_traits<X>::is_neg(v); }
-template <typename X> inline bool is_int(const X & v) { return numeric_traits<X>::is_int(v); }
+template <typename X> inline bool is_integer(const X & v) { return numeric_traits<X>::is_int(v); }
 
 template <typename X> inline X ceil_ratio(const X & a, const X & b) { return numeric_traits<X>::ceil_ratio(a, b); }
 template <typename X> inline X floor_ratio(const X & a, const X & b) { return numeric_traits<X>::floor_ratio(a, b); }