follow the smalles branch

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

correction in the sign of gomory_cut

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

fix in the gomory cut sign

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

try using lemmas of cut_solver as cuts

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

add find_cube() proposed by Nikolaj

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

restore m_int_branch_cut_solver to 8

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

accept empty lar_terms in theory_lra and also do not create empty lar_terms/lemmas

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

qflia_tactic

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

call find_feasible solution to recover for a failure in find_cube

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

do not tighten unused terms

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

get rid of inf_int_set

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

fix a bug with an accidental solution in cube

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

get rid of inf_int_set

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

bug fix with has_int_var() for lar_solver

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

fix in find_inf_int_base_column

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

src/smt/theory_lra.cpp

@@ -1095,7 +1095,7 @@ public:
     }
 
     void init_variable_values() {
-        if (m_solver.get() && th.get_num_vars() > 0) {
+        if (!m.canceled() && m_solver.get() && th.get_num_vars() > 0) {
             m_solver->get_model(m_variable_values);
         }
     }
@@ -1255,7 +1255,7 @@ public:
         return FC_GIVEUP;
     }
 
-    // create a bound atom representing term <= k
+    // create a bound atom representing term >= k is lower_bound is true, and term <= k if it is false
     app_ref mk_bound(lp::lar_term const& term, rational const& k, bool lower_bound) {
         app_ref t = mk_term(term, k.is_int());
         app_ref atom(m);
@@ -1289,7 +1289,7 @@ public:
         case lp::lia_move::ok:
             return l_true;
         case lp::lia_move::branch: {
-            app_ref b = mk_bound(term, k, upper);
+            app_ref b = mk_bound(term, k, !upper);
             // branch on term >= k + 1
             // branch on term <= k
             // TBD: ctx().force_phase(ctx().get_literal(b));
@@ -1300,7 +1300,7 @@ public:
         case lp::lia_move::cut: {
             ++m_stats.m_gomory_cuts;
             // m_explanation implies term <= k
-            app_ref b = mk_bound(term, k, upper);
+            app_ref b = mk_bound(term, k, false);
             m_eqs.reset();
             m_core.reset();
             m_params.reset();
@@ -2683,10 +2683,14 @@ public:
         if (!term.m_v.is_zero()) {
             args.push_back(a.mk_numeral(term.m_v, is_int));
         }
-        if (args.size() == 1) {
+        switch (args.size()) {
+        case 0:
+            return app_ref(a.mk_numeral(rational::zero(), is_int), m);
+        case 1:
             return app_ref(to_app(args[0].get()), m);
+        default:
+            return app_ref(a.mk_add(args.size(), args.c_ptr()), m);
         }
-        return app_ref(a.mk_add(args.size(), args.c_ptr()), m);
     }
 
     app_ref mk_obj(theory_var v) {
@@ -2806,12 +2810,14 @@ public:
         st.update("arith-make-feasible", m_solver->settings().st().m_make_feasible);
         st.update("arith-max-columns", m_solver->settings().st().m_max_cols);
         st.update("arith-max-rows", m_solver->settings().st().m_max_rows);
-        st.update("cut_solver-calls", m_solver->settings().st().m_cut_solver_calls);
+        st.update("cut-solver-calls", m_solver->settings().st().m_cut_solver_calls);
-        st.update("cut_solver-true", m_solver->settings().st().m_cut_solver_true);
+        st.update("cut-solver-true", m_solver->settings().st().m_cut_solver_true);
-        st.update("cut_solver-false", m_solver->settings().st().m_cut_solver_false);
+        st.update("cut-solver-false", m_solver->settings().st().m_cut_solver_false);
-        st.update("cut_solver-undef", m_solver->settings().st().m_cut_solver_undef);
+        st.update("cut-solver-undef", m_solver->settings().st().m_cut_solver_undef);
-        st.update("gcd_calls", m_solver->settings().st().m_gcd_calls);
+        st.update("gcd-calls", m_solver->settings().st().m_gcd_calls);
-        st.update("gcd_conflict", m_solver->settings().st().m_gcd_conflicts);
+        st.update("gcd-conflict", m_solver->settings().st().m_gcd_conflicts);
+        st.update("cube-calls", m_solver->settings().st().m_cube_calls);
+        st.update("cube-success", m_solver->settings().st().m_cube_success);
     }        
 };
     

src/util/lp/CMakeLists.txt

@@ -23,7 +23,6 @@ z3_add_component(lp
     matrix.cpp
     nra_solver.cpp
     permutation_matrix.cpp
-    quick_xplain.cpp
     row_eta_matrix.cpp
     scaler.cpp
     sparse_matrix.cpp

src/util/lp/cut_solver.h

@@ -33,7 +33,7 @@
 #include "util/lp/constraint.h"
 #include "util/lp/active_set.h"
 #include "util/lp/indexer_of_constraints.h"
-
+#include "util/lp/lar_term.h"
 namespace lp {
 class cut_solver; //forward definition
 
@@ -397,6 +397,7 @@ public:
     unsigned                                       m_pushes_to_trail;
     indexer_of_constraints                         m_indexer_of_constraints;
     
+    
     bool is_lower_bound(literal & l) const {
         return l.is_lower();
     }
@@ -2132,11 +2133,11 @@ public:
         }
         
         if (number_of_lemmas == m_lemmas_container.size()) {
-            if (lemma_has_been_modified) {
+            if (lemma_has_been_modified && lemma->poly().number_of_monomials() != 0) {
                 add_lemma_as_not_active(lemma);
             }
             else {
-                delete_constraint(lemma);;
+                delete_constraint(lemma);
             }
         }
     }
@@ -2420,6 +2421,7 @@ public:
     }
 
     void add_lemma_common(constraint* lemma) {
+		lp_assert(lemma->poly().number_of_monomials() > 0);
         m_lemmas_container.add_lemma(lemma);
         polynomial & p = lemma->poly();
         simplify_ineq(p);
@@ -2485,7 +2487,39 @@ public:
     }
 
     unsigned number_of_constraints() const { return m_asserts.size() + m_lemmas_container.size(); }
+
+    void copy_poly_coeffs_to_term(polynomial& poly, lar_term & t) {
+        for (auto & p :poly.m_coeffs)
+            t.add_monomial(p.coeff(), p.var());
+    }
     
+    bool try_getting_cut(lar_term& t, mpq &k, vector<impq>& x) {
+	// todo : create an efficient version based on var_info bounds
+
+        vector<constraint*> short_lemmas;
+        unsigned size = static_cast<unsigned>(-1);
+        for (constraint *c : m_lemmas_container.m_lemmas ) {
+            if (!c->is_ineq()) continue;
+            const auto & p = c->poly();
+            if (p.number_of_monomials() > size)
+                continue;
+            if (is_pos(c->poly().value(x))) {
+                if (p.number_of_monomials() < size) {
+                    size = p.number_of_monomials();
+                    // even is size == 1 it makes sense to look for a random cut
+                    short_lemmas.clear();
+                }
+                short_lemmas.push_back(c);
+            }
+        }
+        unsigned n = short_lemmas.size();
+        if (n == 0) return false;
+        constraint *c = short_lemmas[m_settings.random_next() % n];    
+        k = - c->poly().m_a;
+        copy_poly_coeffs_to_term(c->poly(), t);
+        TRACE("cut_solver_cuts", trace_print_constraint(tout, *c););
+        return true;
+    }
 };
 
 inline polynomial operator*(const mpq & a, polynomial & p) {

src/util/lp/int_solver.cpp

@@ -44,29 +44,35 @@ void int_solver::trace_inf_rows() const {
           );
 }
 
-int_set& int_solver::inf_int_set() {
-    return m_lar_solver->m_inf_int_set;
-}
-
-const int_set& int_solver::inf_int_set() const {
-    return m_lar_solver->m_inf_int_set;
-}
-
 bool int_solver::has_inf_int() const {
-    return !inf_int_set().is_empty();
+    return m_lar_solver->has_inf_int();
 }
 
 int int_solver::find_inf_int_base_column() {
-    if (inf_int_set().is_empty())
+    unsigned inf_int_count;
+    int j = find_inf_int_boxed_base_column_with_smallest_range(inf_int_count);
+	if (j != -1)
+		return j;
+    if (inf_int_count == 0)
         return -1;
-    int j = find_inf_int_boxed_base_column_with_smallest_range();
+    unsigned k = random() % inf_int_count;
-    if (j != -1)
+    return get_kth_inf_int(k);
-        return j;
-    unsigned k = settings().random_next() % inf_int_set().m_index.size();
-    return inf_int_set().m_index[k];
 }
 
-int int_solver::find_inf_int_boxed_base_column_with_smallest_range() {
+int int_solver::get_kth_inf_int(unsigned k) const {
+    unsigned inf_int_count = 0;
+    for (unsigned j : m_lar_solver->r_basis()) {
+        if (! column_is_int_inf(j) )
+            continue;
+        if (inf_int_count++ == k)
+            return j;
+    }
+    lp_assert(false);
+    return -1;
+}
+
+int int_solver::find_inf_int_boxed_base_column_with_smallest_range(unsigned & inf_int_count) {
+    inf_int_count = 0;
     int result = -1;
     mpq range;
     mpq new_range;
@@ -74,11 +80,13 @@ int int_solver::find_inf_int_boxed_base_column_with_smallest_range() {
     unsigned n;
     lar_core_solver & lcs = m_lar_solver->m_mpq_lar_core_solver;
 
-    for (int j : inf_int_set().m_index) {
+    for (unsigned j : m_lar_solver->r_basis()) {
-        lp_assert(is_base(j) && column_is_int_inf(j));
+        if (!column_is_int_inf(j))
-        lp_assert(!is_fixed(j));
+            continue;
+        inf_int_count++;
         if (!is_boxed(j))
             continue;
+        lp_assert(!is_fixed(j));
         new_range  = lcs.m_r_upper_bounds()[j].x - lcs.m_r_lower_bounds()[j].x;
         if (new_range > small_range_thresold) 
             continue;
@@ -325,7 +333,7 @@ lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsi
           tout << "inf_col = " << inf_col << std::endl;
           );
         
-    // gomory will be   t >= k
+    // gomory will be   t <= k and the current solution has a property t > k
     k = 1;
     mpq lcm_den(1);
     unsigned x_j;
@@ -356,7 +364,7 @@ lia_move int_solver::mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsi
 
     lp_assert(current_solution_is_inf_on_cut(t, k));
     m_lar_solver->subs_term_columns(t);
-    TRACE("gomory_cut", tout<<"precut:"; m_lar_solver->print_term(t, tout); tout << ">= " << k << std::endl;);
+    TRACE("gomory_cut", tout<<"precut:"; m_lar_solver->print_term(t, tout); tout << " <= " << k << std::endl;);
     return lia_move::cut;
 }
 
@@ -476,15 +484,69 @@ void int_solver::copy_values_from_cut_solver() {
 }
 
 void int_solver::catch_up_in_adding_constraints_to_cut_solver() {
-    lp_assert(m_cut_solver.number_of_asserts() <= m_lar_solver->constraints().size());
+	lp_assert(m_cut_solver.number_of_asserts() <= m_lar_solver->constraints().size());
     for (unsigned j = m_cut_solver.number_of_asserts(); j < m_lar_solver->constraints().size(); j++) {
         add_constraint_to_cut_solver(j, m_lar_solver->constraints()[j]);
     }
 }
 
+bool int_solver::tighten_term_for_cube(unsigned i) {
+    unsigned ti = i + m_lar_solver->terms_start_index();
+    if (!m_lar_solver->term_is_used_as_row(ti))
+        return true;
+    const lar_term* t = m_lar_solver->terms()[i];
+    mpq delta = zero_of_type<mpq>();
+    for (const auto & p : *t) {
+        delta += abs(p.coeff());
+    }
+    delta *= mpq(1, 2);
+    TRACE("cube", m_lar_solver->print_term_as_indices(*t, tout); tout << ", delta = " << delta;);
+
+    return m_lar_solver->tighten_term_bounds_by_delta(i, delta);
+}
+
+bool int_solver::tighten_terms_for_cube() {
+    for (unsigned i = 0; i < m_lar_solver->terms().size(); i++)
+        if (!tighten_term_for_cube(i)) {
+            TRACE("cube", tout << "cannot tighten";);
+            return false;
+        }
+    return true;
+}
+
+bool int_solver::find_cube() {
+	if (m_branch_cut_counter % settings().m_int_branch_find_cube != 0)
+        return false;
+    
+    settings().st().m_cube_calls++;
+    TRACE("cube",
+          for (unsigned j = 0; j < m_lar_solver->A_r().column_count(); j++)
+              display_column(tout, j);
+          m_lar_solver->print_terms(tout);
+          );
+    m_lar_solver->push();
+    if(!tighten_terms_for_cube()) {
+        m_lar_solver->pop();
+        return false;
+    }
+
+    lp_status st = m_lar_solver->find_feasible_solution();
+    if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
+        TRACE("cube", tout << "cannot find a feasiblie solution";);
+        m_lar_solver->pop();
+        move_non_basic_columns_to_bounds();
+        m_lar_solver->find_feasible_solution();
+        lp_assert(m_lar_solver->get_status() == lp_status::OPTIMAL);
+        // it can happen that we found an integer solution here
+        return !m_lar_solver->r_basis_has_inf_int();
+    }
+    m_lar_solver->round_to_integer_solution();
+    m_lar_solver->pop();
+    return true;
+}
+
 lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) {
     init_check_data();
-    lp_assert(inf_int_set_is_correct());
     // it is a reimplementation of 
     // final_check_status theory_arith<Ext>::check_int_feasibility()
     // from theory_arith_int.h with the addition of cut_solver
@@ -507,7 +569,14 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) {
     if (!has_inf_int())
         return lia_move::ok;
 
-    if ((++m_branch_cut_counter) % settings().m_int_branch_cut_solver == 0) {
+    ++m_branch_cut_counter;
+    if (find_cube()){
+        settings().st().m_cube_success++;
+        return lia_move::ok;
+    }
+    TRACE("cube", tout << "cube did not succeed";);
+    
+    if ((m_branch_cut_counter) % settings().m_int_branch_cut_solver == 0) {
         TRACE("check_main_int", tout<<"cut_solver";);
         catch_up_in_adding_constraints_to_cut_solver();
         auto check_res = m_cut_solver.check();
@@ -520,9 +589,18 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) {
         case cut_solver::lbool::l_true:
             settings().st().m_cut_solver_true++;
             copy_values_from_cut_solver();
+            lp_assert(m_lar_solver->all_constraints_hold());
             return lia_move::ok;
         case cut_solver::lbool::l_undef:
             settings().st().m_cut_solver_undef++;
+            if (m_cut_solver.try_getting_cut(t, k, m_lar_solver->m_mpq_lar_core_solver.m_r_x)) {
+                m_lar_solver->subs_term_columns(t);
+                TRACE("cut_solver_cuts",
+                      tout<<"precut from cut_solver:"; m_lar_solver->print_term(t, tout); tout << " <= " << k << std::endl;);
+
+
+                return lia_move::cut;
+            }
             break;
         default:
             return lia_move::give_up;
@@ -598,7 +676,6 @@ void int_solver::set_value_for_nbasic_column_ignore_old_values(unsigned j, const
     auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x[j];
     auto delta = new_val - x;
     x = new_val;
-    update_column_in_int_inf_set(j);
     m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta);
 }
 
@@ -612,7 +689,6 @@ void int_solver::set_value_for_nbasic_column(unsigned j, const impq & new_val) {
     }
     auto delta = new_val - x;
     x = new_val;
-    update_column_in_int_inf_set(j);
     m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta);
 }
 
@@ -677,7 +753,6 @@ void int_solver::patch_int_infeasible_nbasic_columns() {
         move_non_basic_columns_to_bounds();
         m_lar_solver->find_feasible_solution();
     }
-    lp_assert(is_feasible() && inf_int_set_is_correct());
 }
 
 mpq get_denominators_lcm(const row_strip<mpq> & row) {
@@ -990,26 +1065,9 @@ void int_solver::display_column(std::ostream & out, unsigned j) const {
     m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(j, out);
 }
 
-bool int_solver::inf_int_set_is_correct() const {
-    for (unsigned j = 0; j < m_lar_solver->A_r().column_count(); j++) {
-        if (inf_int_set().contains(j) != (is_int(j) && (!value_is_int(j)))) {
-            TRACE("arith_int", tout << "j= " << j << " inf_int_set().contains(j) = " << inf_int_set().contains(j) << ", is_int(j) = "   << is_int(j) << "\nvalue_is_int(j) = " << value_is_int(j) << ", val = " << get_value(j) << std::endl;); 
-            return false;
-        }
-    }
-    return true;
-}
-
 bool int_solver::column_is_int_inf(unsigned j) const {
     return is_int(j) && (!value_is_int(j));
 }
-    
-void int_solver::update_column_in_int_inf_set(unsigned j) {
-    if (is_int(j) && (!value_is_int(j)))
-        inf_int_set().insert(j);
-    else
-        inf_int_set().erase(j);
-}
 
 bool int_solver::is_base(unsigned j) const {
     return m_lar_solver->m_mpq_lar_core_solver.m_r_heading[j] >= 0;
@@ -1181,6 +1239,7 @@ const impq& int_solver::lower_bound(unsigned j) const {
 }
 
 lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column, bool & upper) {
+    TRACE("check_main_int", tout << "branching" << std::endl;);
     lp_assert(t.is_empty());
     lp_assert(j != -1);
     t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j));
@@ -1201,7 +1260,6 @@ lia_move int_solver::create_branch_on_column(int j, lar_term& t, mpq& k, bool fr
 }
 
 bool int_solver::left_branch_is_more_narrow_than_right(unsigned j) {
-    return settings().random_next() % 2;
     switch (m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_column_types[j] ) {
     case column_type::fixed:
         return false;
@@ -1222,16 +1280,7 @@ bool int_solver::left_branch_is_more_narrow_than_right(unsigned j) {
 const impq& int_solver::upper_bound(unsigned j) const {
     return m_lar_solver->column_upper_bound(j);
 }
-void int_solver::display_inf_or_int_inf_columns(std::ostream & out) const {
+
-    out << "int inf\n";
-    for (unsigned j : m_lar_solver->m_inf_int_set.m_index) {
-        display_column(out, j);
-    }
-    out << "regular inf\n";
-    for (unsigned j : m_lar_solver->m_mpq_lar_core_solver.m_r_solver.m_inf_set.m_index) {
-        display_column(out, j);
-    }
-}
 bool int_solver::is_term(unsigned j) const {
     return m_lar_solver->column_corresponds_to_term(j);
 }
@@ -1250,8 +1299,8 @@ void int_solver::pop(unsigned k) {
     m_cut_solver.pop_constraints();
 }
 
-void int_solver::push() {
+void int_solver::push() { m_cut_solver.push(); }
-    m_cut_solver.push();
+
-}
+unsigned int_solver::column_count() const  { return m_lar_solver->column_count(); }
 
 }

src/util/lp/int_solver.h

@@ -59,8 +59,7 @@ public:
     cut_solver m_cut_solver;
     // methods
     int_solver(lar_solver* lp);
-    int_set& inf_int_set();
+
-    const int_set& inf_int_set() const;
     // main function to check that the solution provided by lar_solver is valid for integral values,
     // or provide a way of how it can be adjusted.
     lia_move check(lar_term& t, mpq& k, explanation& ex, bool & upper);
@@ -115,13 +114,11 @@ private:
     void failed();
     bool is_feasible() const;
     const impq & get_value(unsigned j) const;
-    void display_column(std::ostream & out, unsigned j) const;
-    bool inf_int_set_is_correct() const;
-    void update_column_in_int_inf_set(unsigned j);
     bool column_is_int_inf(unsigned j) const;
     void trace_inf_rows() const;
     int find_inf_int_base_column();
-    int find_inf_int_boxed_base_column_with_smallest_range();
+    int find_inf_int_boxed_base_column_with_smallest_range(unsigned&);
+    int get_kth_inf_int(unsigned) const;
     lp_settings& settings();
     bool move_non_basic_columns_to_bounds();
     void branch_infeasible_int_var(unsigned);
@@ -143,6 +140,7 @@ private:
     }
 
 public:
+    void display_column(std::ostream & out, unsigned j) const;
     inline static
     mpq fractional_part(const impq & n) {
         lp_assert(is_rational(n));
@@ -164,7 +162,6 @@ private:
     lia_move create_branch_on_column(int j, lar_term& t, mpq& k, bool free_column, bool & upper);
     void catch_up_in_adding_constraints_to_cut_solver();
 public:
-    void display_inf_or_int_inf_columns(std::ostream & out) const;
     template <typename T>
     void fill_cut_solver_vars();
     template <typename T>
@@ -176,5 +173,9 @@ public:
     void push();
     void copy_values_from_cut_solver();
     bool left_branch_is_more_narrow_than_right(unsigned);
+    bool find_cube();
+    bool tighten_terms_for_cube();
+    bool tighten_term_for_cube(unsigned);
+    unsigned column_count() const;
 };
 }

src/util/lp/lar_core_solver.h

@@ -830,5 +830,6 @@ public:
         }
     }
 
+    const vector<unsigned>& r_basis() const { return m_r_basis; }
 };
 }

src/util/lp/lar_solver.cpp

@@ -31,11 +31,8 @@ lar_solver::lar_solver() : m_status(lp_status::OPTIMAL),
                            m_infeasible_column_index(-1),
                            m_terms_start_index(1000000),
                            m_mpq_lar_core_solver(m_settings, *this),
-                           m_tracker_of_x_change([&](unsigned j) {
+                           m_int_solver(nullptr),
-                                   call_assignment_tracker(j);
+                           m_has_int_var(false)
-                               }
-                               ),
-                           m_int_solver(nullptr)
 {}
     
 void lar_solver::set_track_pivoted_rows(bool v) {
@@ -277,11 +274,10 @@ lp_status lar_solver::get_status() const { return m_status;}
 void lar_solver::set_status(lp_status s) {m_status = s;}
 
 bool lar_solver::has_int_var() const {
-    return m_mpq_lar_core_solver.m_r_solver.m_tracker_of_x_change != nullptr;
+    return m_has_int_var;
 }
 
 lp_status lar_solver::find_feasible_solution() {
-    lp_assert(inf_int_set_is_correct());
     m_settings.st().m_make_feasible++;
     if (A_r().column_count() > m_settings.st().m_max_cols)
         m_settings.st().m_max_cols = A_r().column_count();
@@ -290,19 +286,12 @@ lp_status lar_solver::find_feasible_solution() {
     if (strategy_is_undecided())
         decide_on_strategy_and_adjust_initial_state();
 
-    if (has_int_var()) {
-        m_inf_int_set.resize(A_r().column_count());
-    }
-    
     m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = true;
     auto ret = solve();
-    TRACE("intinf", m_int_solver->display_inf_or_int_inf_columns(tout););
-    lp_assert(inf_int_set_is_correct());
     return ret;
 }
     
 lp_status lar_solver::solve() {
-    lp_assert(inf_int_set_is_correct());
     if (m_status == lp_status::INFEASIBLE) {
         return m_status;
     }
@@ -313,7 +302,6 @@ lp_status lar_solver::solve() {
     }
        
     m_columns_with_changed_bound.clear();
-    lp_assert(inf_int_set_is_correct());
     return m_status;
 }
 
@@ -363,7 +351,6 @@ void lar_solver::shrink_inf_set_after_pop(unsigned n, int_set & set) {
     
 void lar_solver::pop(unsigned k) {
     TRACE("arith_int", tout << "pop" << std::endl;);
-    lp_assert(inf_int_set_is_correct());
     TRACE("lar_solver", tout << "k = " << k << std::endl;);
 
     int n_was = static_cast<int>(m_ext_vars_to_columns.size());
@@ -381,11 +368,8 @@ void lar_solver::pop(unsigned k) {
 
     m_mpq_lar_core_solver.pop(k);
     clean_popped_elements(n, m_columns_with_changed_bound);
-    clean_popped_elements(n, m_inf_int_set);
     unsigned m = A_r().row_count();
-    lp_assert(inf_int_set_is_correct());
     clean_popped_elements(m, m_rows_with_changed_bounds);
-    lp_assert(inf_int_set_is_correct());
     clean_inf_set_of_r_solver_after_pop();
     m_status = m_mpq_lar_core_solver.m_r_solver.current_x_is_feasible()? lp_status::OPTIMAL: lp_status::UNKNOWN;
     lp_assert(m_settings.simplex_strategy() == simplex_strategy_enum::undecided ||
@@ -663,7 +647,6 @@ bool lar_solver::costs_are_used() const {
 }
     
 void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair<mpq> & delta) {
-    lp_assert(inf_int_set_is_correct());
     if (use_tableau()) {
         for (const auto & c : A_r().m_columns[j]) {
             unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.m_i];
@@ -685,7 +668,6 @@ void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j,
             m_mpq_lar_core_solver.m_r_solver.update_x_with_delta_and_track_feasibility(bj, -m_column_buffer[i] * delta); 
         }
     }
-    lp_assert(inf_int_set_is_correct());
 }
 
 void lar_solver::update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j) {
@@ -755,10 +737,8 @@ void lar_solver::solve_with_core_solver() {
         update_x_and_inf_costs_for_columns_with_changed_bounds_tableau();
     else 
         update_x_and_inf_costs_for_columns_with_changed_bounds();
-    TRACE("intinf", m_int_solver->display_inf_or_int_inf_columns(tout););
     m_mpq_lar_core_solver.solve();
     set_status(m_mpq_lar_core_solver.m_r_solver.get_status());
-    lp_assert(inf_int_set_is_correct());
     lp_assert(m_status != lp_status::OPTIMAL || all_constraints_hold());
 }
 
@@ -1140,8 +1120,8 @@ void lar_solver::get_infeasibility_explanation_for_inf_sign(
 }
 
 void lar_solver::get_model(std::unordered_map<var_index, mpq> & variable_values) const {
-    mpq delta = mpq(1, 2); // start from 0.5 to have less clashes
     lp_assert(m_status == lp_status::OPTIMAL);
+    mpq delta = mpq(1, 2); // start from 0.5 to have less clashes
     unsigned i;
     do {
         // different pairs have to produce different singleton values
@@ -1276,7 +1256,6 @@ void lar_solver::random_update(unsigned sz, var_index const * vars) {
     fill_var_set_for_random_update(sz, vars, column_list);
     random_updater ru(*this, column_list);
     ru.update();
-    lp_assert(inf_int_set_is_correct());
 }
 
 
@@ -1406,7 +1385,6 @@ void lar_solver::pop_tableau() {
 }
 
 void lar_solver::clean_inf_set_of_r_solver_after_pop() {
-    lp_assert(inf_int_set_is_correct());
     vector<unsigned> became_feas;
     clean_popped_elements(A_r().column_count(), m_mpq_lar_core_solver.m_r_solver.m_inf_set);
     std::unordered_set<unsigned> basic_columns_with_changed_cost;
@@ -1449,12 +1427,6 @@ void lar_solver::clean_inf_set_of_r_solver_after_pop() {
     }
 }
 
-void lar_solver::shrink_explanation_to_minimum(vector<std::pair<mpq, constraint_index>> & explanation) const {
-    // implementing quickXplain
-    quick_xplain::run(explanation, *this);
-    lp_assert(this->explanation_is_correct(explanation));
-}
-
 bool lar_solver::model_is_int_feasible() const {
     unsigned n = A_r().column_count();
     for (unsigned j = 0; j < n; j++) {
@@ -1495,7 +1467,7 @@ bool lar_solver::column_is_fixed(unsigned j) const {
 bool lar_solver::ext_var_is_int(var_index ext_var) const {
     auto it = m_ext_vars_to_columns.find(ext_var);
     lp_assert(it != m_ext_vars_to_columns.end());
-    return it == m_ext_vars_to_columns.end() || it->second.is_integer();
+    return it->second.is_integer();
 }
 
 // below is the initialization functionality of lar_solver
@@ -1511,6 +1483,8 @@ void lar_solver::catch_up_in_updating_int_solver() {
 }
 
 var_index lar_solver::add_var(unsigned ext_j, bool is_int) {
+    if (is_int)
+        m_has_int_var = true;
     if (is_int && !has_int_var())
         catch_up_in_updating_int_solver();
         
@@ -1529,10 +1503,6 @@ var_index lar_solver::add_var(unsigned ext_j, bool is_int) {
     m_columns_to_ul_pairs.push_back(ul_pair(static_cast<unsigned>(-1)));
     add_non_basic_var_to_core_fields(ext_j, is_int);
     lp_assert(sizes_are_correct());
-    if (is_int) {
-        m_mpq_lar_core_solver.m_r_solver.set_tracker_of_x(& m_tracker_of_x_change);
-    }
-    lp_assert(inf_int_set_is_correct());
     return i;
 }
 
@@ -1548,7 +1518,6 @@ void lar_solver::add_non_basic_var_to_core_fields(unsigned ext_j, bool is_int) {
     register_new_ext_var_index(ext_j, is_int);
     m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
     m_columns_with_changed_bound.increase_size_by_one();
-    m_inf_int_set.increase_size_by_one();
     add_new_var_to_core_fields_for_mpq(false);
     if (use_lu())
         add_new_var_to_core_fields_for_doubles(false);
@@ -1686,7 +1655,6 @@ void lar_solver::add_row_from_term_no_constraint(const lar_term * term, unsigned
     m_mpq_lar_core_solver.m_r_solver.update_x_and_call_tracker(j, get_basic_var_value_from_row_directly(A_r().row_count() - 1));
     if (use_lu())
         fill_last_row_of_A_d(A_d(), term);
-    lp_assert(inf_int_set_is_correct());
 }
 
 void lar_solver::add_basic_var_to_core_fields() {
@@ -1695,7 +1663,6 @@ void lar_solver::add_basic_var_to_core_fields() {
     m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
     m_columns_with_changed_bound.increase_size_by_one();
     m_rows_with_changed_bounds.increase_size_by_one();
-    m_inf_int_set.increase_size_by_one();
     add_new_var_to_core_fields_for_mpq(true);
     if (use_lu)
         add_new_var_to_core_fields_for_doubles(true);
@@ -1720,7 +1687,6 @@ constraint_index lar_solver::add_var_bound(var_index j, lconstraint_kind kind, c
         add_var_bound_on_constraint_for_term(j, kind, right_side, ci);
     }
     lp_assert(sizes_are_correct());
-    lp_assert(inf_int_set_is_correct());
     return ci;
 }
 
@@ -1749,7 +1715,7 @@ void lar_solver::update_column_type_and_bound(var_index j, lconstraint_kind kind
 void lar_solver::add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
     lp_assert(is_term(j));
     unsigned adjusted_term_index = adjust_term_index(j);
-    lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int());
+    //    lp_assert(!term_is_int(m_terms[adjusted_term_index]) || right_side.is_int());
     auto it = m_ext_vars_to_columns.find(j);
     if (it != m_ext_vars_to_columns.end()) {
         unsigned term_j = it->second.internal_j();
@@ -2164,6 +2130,87 @@ var_index lar_solver:: to_var_index(unsigned ext_j) const {
     lp_assert(it != m_ext_vars_to_columns.end());
     return it->second.internal_j();
 }
+
+bool lar_solver::tighten_term_bounds_by_delta(unsigned term_index, const mpq& delta) {
+    unsigned tj = term_index + m_terms_start_index;
+	auto it = m_ext_vars_to_columns.find(tj);
+	if (it == m_ext_vars_to_columns.end())
+		return true;
+	unsigned j = it->second.internal_j();
+    auto & slv = m_mpq_lar_core_solver.m_r_solver;
+    TRACE("cube", tout << "delta = " << delta << std::endl;
+          m_int_solver->display_column(tout, j); );
+    if (slv.column_has_upper_bound(j) && slv.column_has_lower_bound(j)) {
+        if (slv.m_upper_bounds[j].x - delta < slv.m_lower_bounds[j].x + delta) {
+            TRACE("cube", tout << "cannot tighten, delta = " << delta;);
+            return false;
+        }
+    }
+    TRACE("cube", tout << "can tighten";);
+    if (slv.column_has_upper_bound(j)) {
+        add_var_bound(tj, lconstraint_kind::LE, slv.m_upper_bounds[j].x - delta);
+    }
+    if (slv.column_has_lower_bound(j)) {
+        add_var_bound(tj, lconstraint_kind::GE, slv.m_lower_bounds[j].x + delta);
+    }
+    return true;
+}
+
+void lar_solver::update_delta_for_terms(const impq & delta, unsigned j, const vector<unsigned>& terms_of_var) {
+    for (unsigned i : terms_of_var) {
+        lar_term & t = *m_terms[i];
+        auto it = t.m_coeffs.find(j);
+        unsigned tj = to_var_index(i + m_terms_start_index);
+        TRACE("cube",
+              tout << "t.apply = " << t.apply(m_mpq_lar_core_solver.m_r_x) << ", m_mpq_lar_core_solver.m_r_x[tj]= " << m_mpq_lar_core_solver.m_r_x[tj];);
+        TRACE("cube", print_term_as_indices(t, tout);
+              tout << ", it->second = " << it->second;
+              tout << ", tj = " << tj << ", ";
+              m_int_solver->display_column(tout, tj);
+              );
+        
+        m_mpq_lar_core_solver.m_r_x[tj] += it->second * delta;
+        lp_assert(t.apply(m_mpq_lar_core_solver.m_r_x) == m_mpq_lar_core_solver.m_r_x[tj]);
+        TRACE("cube", m_int_solver->display_column(tout, tj); );
+    }
+}
+
+
+void lar_solver::fill_vars_to_terms(vector<vector<unsigned>> & vars_to_terms) {
+	for (unsigned j = 0; j < m_terms.size(); j++) {
+		if (!term_is_used_as_row(j + m_terms_start_index))
+			continue;
+		for (const auto & p : *m_terms[j]) {
+			if (p.var() >= vars_to_terms.size())
+				vars_to_terms.resize(p.var() + 1);
+			vars_to_terms[p.var()].push_back(j);
+		}
+	}
+}
+
+void lar_solver::round_to_integer_solution() {
+    vector<vector<unsigned>> vars_to_terms;
+    fill_vars_to_terms(vars_to_terms);
+
+    for (unsigned j = 0; j < column_count(); j++) {
+        if (column_corresponds_to_term(j)) continue;
+        TRACE("cube", m_int_solver->display_column(tout, j););
+        impq& v =  m_mpq_lar_core_solver.m_r_x[j];
+        if (v.is_int())
+            continue;
+        impq flv = floor(v);
+        auto del = flv - v; // del is negative
+        if (del < - mpq(1, 2)) {
+            del = impq(one_of_type<mpq>()) + del;
+            v = ceil(v);
+        } else {
+            v = flv;
+        }
+        TRACE("cube", m_int_solver->display_column(tout, j); tout << "v = " << v << " ,del = " << del;);
+        update_delta_for_terms(del, j, vars_to_terms[j]);
+    }
+}
+
 } // namespace lp
 
 

src/util/lp/lar_solver.h

@@ -39,7 +39,6 @@ Revision History:
 #include "util/lp/stacked_unordered_set.h"
 #include "util/lp/implied_bound.h"
 #include "util/lp/bound_analyzer_on_row.h"
-#include "util/lp/quick_xplain.h"
 #include "util/lp/conversion_helper.h"
 #include "util/lp/int_solver.h"
 #include "util/lp/nra_solver.h"
@@ -117,17 +116,19 @@ private:
     vector<lar_term*>                                   m_terms;
     const var_index                                     m_terms_start_index;
     indexed_vector<mpq>                                 m_column_buffer;
+public:
+    lar_core_solver                                     m_mpq_lar_core_solver;
+private:
+    int_solver *                                        m_int_solver;
+    bool                                                m_has_int_var;
 
     
 public :
+    unsigned terms_start_index() const { return m_terms_start_index; }
+    const vector<lar_term*> terms() const { return m_terms; }
     const vector<lar_base_constraint*>& constraints() const {
         return m_constraints;
     }
-    lar_core_solver m_mpq_lar_core_solver;
-private:
-    std::function<void (unsigned)> m_tracker_of_x_change;
-    int_solver * m_int_solver;
-public:
     void set_int_solver(int_solver * int_slv) {
         m_int_solver = int_slv;
     }
@@ -136,7 +137,6 @@ public:
     }
     unsigned constraint_count() const;
     const lar_base_constraint& get_constraint(unsigned ci) const;
-    int_set m_inf_int_set; 
     ////////////////// methods ////////////////////////////////
     static_matrix<mpq, numeric_pair<mpq>> & A_r();
     static_matrix<mpq, numeric_pair<mpq>> const & A_r() const;
@@ -551,43 +551,32 @@ public:
         }
     }
 
-    bool inf_int_set_is_correct_for_column(unsigned j) const {
-        if (m_inf_int_set.contains(j) != (column_is_int(j) && (!column_value_is_integer(j)))) {
-            TRACE("arith_int",
-                  tout << "j= " << j <<
-                  " inf_int_set().contains(j) = " << m_inf_int_set.contains(j) <<
-                  ", column_is_int(j) = "   << column_is_int(j) <<
-                  "\n column_value_is_integer(j) = " << column_value_is_integer(j) <<
-                  ", val = " << get_column_value(j) << std::endl;); 
-            return false;
-        }
-        return true;
-    }
-    
-    bool inf_int_set_is_correct() const {
-        if (!has_int_var())
-            return true;
-        for (unsigned j = 0; j < A_r().column_count(); j++) {
-            if (inf_int_set_is_correct_for_column(j) == false)
-                return false;
-        }
-        return true;
-    }
     bool has_int_var() const;
-    void call_assignment_tracker(unsigned j) {
+    bool has_inf_int() const {
-        if (!var_is_int(j)) {
+        for (unsigned j = 0; j < column_count(); j++) {
-            lp_assert(m_inf_int_set.contains(j) == false);
+            if (column_is_int(j) && ! column_value_is_int(j))
-            return;
+                return true;
         }
-        if (m_mpq_lar_core_solver.m_r_x[j].is_int())
+        return false;
-            m_inf_int_set.erase(j);
-        else
-            m_inf_int_set.insert(j);
     }
 
+    bool r_basis_has_inf_int() const {
+        for (unsigned j : r_basis()) {
+            if (column_is_int(j) && ! column_value_is_int(j))
+                return true;
+        }
+        return false;
+    }
+    
     lar_core_solver & get_core_solver() { return m_mpq_lar_core_solver; }
     bool column_corresponds_to_term(unsigned) const;
     void catch_up_in_updating_int_solver();
     var_index to_var_index(unsigned ext_j) const;
+    bool tighten_term_bounds_by_delta(unsigned, const mpq&);
+    void round_to_integer_solution();
+    void update_delta_for_terms(const impq & delta, unsigned j, const vector<unsigned>&);
+    void fill_vars_to_terms(vector<vector<unsigned>> & vars_to_terms);
+    unsigned column_count() const { return A_r().column_count(); }
+    const vector<unsigned> & r_basis() const { return m_mpq_lar_core_solver.r_basis(); }
 };
 }

src/util/lp/lp_core_solver_base.h

@@ -85,12 +85,7 @@ public:
     bool                  m_tracing_basis_changes;
     int_set*              m_pivoted_rows;
     bool                  m_look_for_feasible_solution_only;
-    std::function<void (unsigned)> * m_tracker_of_x_change;
 
-    void set_tracker_of_x(std::function<void (unsigned)>* tracker) {
-        m_tracker_of_x_change = tracker;
-    }
-    
     void start_tracing_basis_changes() {
         m_trace_of_basis_change_vector.resize(0);
         m_tracing_basis_changes = true;
@@ -705,14 +700,10 @@ public:
 
     void update_x_and_call_tracker(unsigned j, const X & v) {
         m_x[j] = v;
-        if (m_tracker_of_x_change != nullptr)
-            (*m_tracker_of_x_change)(j);
     }
 
     void add_delta_to_x_and_call_tracker(unsigned j, const X & delta) {
         m_x[j] += delta;
-        if (m_tracker_of_x_change != nullptr)
-            (*m_tracker_of_x_change)(j);
     }
     
     void track_column_feasibility(unsigned j) {
@@ -722,14 +713,10 @@ public:
             insert_column_into_inf_set(j);
     }
     void insert_column_into_inf_set(unsigned j) {
-        if (m_tracker_of_x_change != nullptr)
-            (*m_tracker_of_x_change)(j);
         m_inf_set.insert(j);
         lp_assert(!column_is_feasible(j));
     }
     void remove_column_from_inf_set(unsigned j) {
-        if (m_tracker_of_x_change != nullptr)
-            (*m_tracker_of_x_change)(j);
         m_inf_set.erase(j);
         lp_assert(column_is_feasible(j));
     }

src/util/lp/lp_core_solver_base_def.h

@@ -67,8 +67,7 @@ lp_core_solver_base(static_matrix<T, X> & A,
     m_steepest_edge_coefficients(A.column_count()),
     m_tracing_basis_changes(false),
     m_pivoted_rows(nullptr),
-    m_look_for_feasible_solution_only(false),
+    m_look_for_feasible_solution_only(false) {
-    m_tracker_of_x_change(nullptr) {
     lp_assert(bounds_for_boxed_are_set_correctly());    
     init();
     init_basis_heading_and_non_basic_columns_vector();

src/util/lp/lp_settings.h

@@ -108,6 +108,8 @@ struct stats {
     unsigned m_cut_solver_undef;
     unsigned m_gcd_calls;
     unsigned m_gcd_conflicts;
+    unsigned m_cube_calls;
+    unsigned m_cube_success;
     stats() { reset(); }
     void reset() { memset(this, 0, sizeof(*this)); }
 };
@@ -334,13 +336,14 @@ public:
     bool use_breakpoints_in_feasibility_search;
     unsigned random_next() { return m_rand(); }
     void set_random_seed(unsigned s) { m_rand.set_seed(s); }
-    unsigned max_row_length_for_bound_propagation;
+    unsigned         max_row_length_for_bound_propagation;
-    bool backup_costs;
+    bool             backup_costs;
-    unsigned column_number_threshold_for_using_lu_in_lar_solver;
+    unsigned         column_number_threshold_for_using_lu_in_lar_solver;
-    unsigned m_int_branch_cut_gomory_threshold;
+    unsigned         m_int_branch_cut_gomory_threshold;
-    unsigned m_int_branch_cut_solver;
+    unsigned         m_int_branch_cut_solver;
-    bool m_run_gcd_test;
+    unsigned         m_int_branch_find_cube;
-    unsigned m_cut_solver_cycle_on_var;
+    bool             m_run_gcd_test;
+    unsigned         m_cut_solver_cycle_on_var;
 }; // end of lp_settings class
 
 

src/util/lp/polynomial.h

@@ -101,6 +101,15 @@ struct polynomial {
         const mpq & a = coeff(j);
         return a == 1 || a == -1;
     }
+
+    template <typename c> // c plays a role of a map from indices to impq
+    mpq value(const c& v) const {
+        mpq r = m_a;
+        for (auto & p : m_coeffs)
+            r += v[p.var()].x * p.coeff();
+        return r;
+    }
+    
     const vector<monomial> & coeffs() const { return m_coeffs; }
 };
 }

src/util/lp/stacked_vector.h

@@ -52,6 +52,14 @@ public:
         bool operator==(B const& other) const {
             return m_vec.m_vector[m_i] == other;
         }
+        B& operator+=(B const &delta) {
+            // not tracking the change here!
+            return m_vec.m_vector[m_i] += delta;
+        }
+        B& operator-=(B const &delta) {
+            // not tracking the change here!
+            return m_vec.m_vector[m_i] -= delta;
+        }
     };
 
     class ref_const {