move hnf cut functionality

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

src/math/lp/CMakeLists.txt

@@ -10,6 +10,7 @@ z3_add_component(lp
     factorization.cpp
     factorization_factory_imp.cpp
     gomory.cpp
+    hnf_cutter.cpp
     horner.cpp
     indexed_vector.cpp
     int_branch.cpp

src/math/lp/hnf_cutter.cpp

@@ -0,0 +1,275 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    hnf_cutter.cpp
+
+Author:
+    Lev Nachmanson (levnach)
+
+--*/
+#include "math/lp/int_solver.h"
+#include "math/lp/lar_solver.h"
+#include "math/lp/hnf_cutter.h"
+
+namespace lp  {
+
+    hnf_cutter::hnf_cutter(int_solver& lia): 
+        lia(lia), 
+        lra(lia.lra),
+        m_settings(lia.settings()),
+        m_abs_max(zero_of_type<mpq>()),
+        m_var_register(0) {}
+    
+    bool hnf_cutter::is_full() const {
+        return
+            terms_count() >= lia.settings().limit_on_rows_for_hnf_cutter ||
+            vars().size() >= lia.settings().limit_on_columns_for_hnf_cutter;
+    }
+
+    void hnf_cutter::clear() {
+        // m_A will be filled from scratch in init_matrix_A
+        m_var_register.clear();
+        m_terms.clear();
+        m_terms_upper.clear();
+        m_constraints_for_explanation.clear();
+        m_right_sides.clear();
+        m_abs_max = zero_of_type<mpq>();
+        m_overflow = false;
+    }
+
+    void hnf_cutter::add_term(const lar_term* t, const mpq &rs, constraint_index ci, bool upper_bound) {
+        m_terms.push_back(t);
+        m_terms_upper.push_back(upper_bound);
+        if (upper_bound)
+            m_right_sides.push_back(rs);
+        else
+            m_right_sides.push_back(-rs);
+        
+        m_constraints_for_explanation.push_back(ci);
+       
+        for (const auto &p : *t) {
+            m_var_register.add_var(p.var());
+            mpq t = abs(ceil(p.coeff()));
+            if (t > m_abs_max)
+                m_abs_max = t;
+        }
+    }
+    
+    void hnf_cutter::print(std::ostream & out) {
+        out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl;
+    }
+
+    void hnf_cutter::initialize_row(unsigned i) {
+        mpq sign = m_terms_upper[i]? one_of_type<mpq>(): - one_of_type<mpq>();
+        m_A.init_row_from_container(i, * m_terms[i], [this](unsigned j) { return m_var_register.add_var(j);}, sign);
+    }
+
+    void hnf_cutter::init_matrix_A() {
+        m_A = general_matrix(terms_count(), vars().size());
+        for (unsigned i = 0; i < terms_count(); i++)
+            initialize_row(i);
+    }
+
+    // todo: as we need only one row i with non integral b[i] need to optimize later
+    void hnf_cutter::find_h_minus_1_b(const general_matrix& H, vector<mpq> & b) {
+        // the solution will be put into b
+        for (unsigned i = 0; i < H.row_count() ;i++) {
+            for (unsigned j = 0; j < i; j++) {
+                b[i] -= H[i][j]*b[j];
+            }
+            b[i] /= H[i][i];
+            // consider return from here if b[i] is not an integer and return i
+        }
+    }
+
+    vector<mpq> hnf_cutter::create_b(const svector<unsigned> & basis_rows) {
+        if (basis_rows.size() == m_right_sides.size())
+            return m_right_sides;
+        vector<mpq> b;
+        for (unsigned i : basis_rows) {
+            b.push_back(m_right_sides[i]);
+        }
+        return b;
+    }
+
+    int hnf_cutter::find_cut_row_index(const vector<mpq> & b) {
+        int ret = -1;
+        int n = 0;
+        for (int i = 0; i < static_cast<int>(b.size()); i++) {
+            if (is_integer(b[i])) continue;
+            if (n == 0 ) {
+                lp_assert(ret == -1);
+                n = 1;
+                ret = i;
+            } else {
+                if (m_settings.random_next() % (++n) == 0) {
+                    ret = i;
+                }
+            }
+        }
+        return ret;
+    }
+
+    // fills e_i*H_minus_1
+    void hnf_cutter::get_ei_H_minus_1(unsigned i, const general_matrix& H, vector<mpq> & row) {
+        // we solve x = ei * H_min_1
+        // or x * H = ei
+        unsigned m = H.row_count();
+        for (unsigned k = i + 1; k < m; k++) {
+            row[k] = zero_of_type<mpq>();
+        }
+        row[i] = one_of_type<mpq>() / H[i][i];
+        for(int k = i - 1; k >= 0; k--) {
+            mpq t = zero_of_type<mpq>();
+            for (unsigned l = k + 1; l <= i; l++) {
+                t += H[l][k]*row[l];
+            }
+            row[k] = -t / H[k][k]; 
+        }
+
+        // // test region
+        // vector<mpq> ei(H.row_count(), zero_of_type<mpq>());
+        // ei[i] = one_of_type<mpq>();
+        // vector<mpq> pr = row * H;
+        // pr.shrink(ei.size());
+        // lp_assert(ei == pr);
+        // // end test region
+        
+    }
+
+    void hnf_cutter::fill_term(const vector<mpq> & row, lar_term& t) {
+        for (unsigned j = 0; j < row.size(); j++) {
+            if (!is_zero(row[j]))
+                t.add_monomial(row[j], m_var_register.local_to_external(j));
+        }
+    }
+#ifdef Z3DEBUG
+    vector<mpq> hnf_cutter::transform_to_local_columns(const vector<impq> & x) const {
+        vector<mpq> ret;
+        for (unsigned j = 0; j < vars().size(); j++) {
+             ret.push_back(x[m_var_register.local_to_external(j)].x);
+        }
+        return ret;
+    }
+#endif
+    void hnf_cutter::shrink_explanation(const svector<unsigned>& basis_rows) {
+        svector<unsigned> new_expl;
+        for (unsigned i : basis_rows) {
+            new_expl.push_back(m_constraints_for_explanation[i]);
+        }
+        m_constraints_for_explanation = new_expl;
+    }
+
+    bool hnf_cutter::overflow() const { return m_overflow; }
+    
+    lia_move hnf_cutter::create_cut(lar_term& t, mpq& k, explanation* ex, bool & upper, const vector<mpq> & x0) {
+        // we suppose that x0 has at least one non integer element 
+        (void)x0;
+
+        init_matrix_A();
+        svector<unsigned> basis_rows;
+        mpq big_number = m_abs_max.expt(3);
+        mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows, big_number);
+        
+        if (d >= big_number) {
+            return lia_move::undef;
+        }
+        
+        if (m_settings.get_cancel_flag()) {
+            return lia_move::undef;
+        }
+
+        if (basis_rows.size() < m_A.row_count()) {
+            m_A.shrink_to_rank(basis_rows);
+            shrink_explanation(basis_rows);
+        }
+        
+        hnf<general_matrix> h(m_A, d);        
+        vector<mpq> b = create_b(basis_rows);
+        lp_assert(m_A * x0 == b);
+        find_h_minus_1_b(h.W(), b);
+        int cut_row = find_cut_row_index(b);
+
+        if (cut_row == -1) {
+            return lia_move::undef;
+        }
+
+        // the matrix is not square - we can get
+        // all integers in b's projection
+        
+        vector<mpq> row(m_A.column_count());
+        get_ei_H_minus_1(cut_row, h.W(), row);
+        vector<mpq> f = row * m_A;
+        fill_term(f, t);
+        k = floor(b[cut_row]);
+        upper = true;
+        return lia_move::cut;
+    }
+
+    svector<unsigned> hnf_cutter::vars() const { return m_var_register.vars(); }
+
+    void hnf_cutter::try_add_term_to_A_for_hnf(unsigned i) {
+        mpq rs;
+        const lar_term* t = lra.terms()[i];
+        constraint_index ci;
+        bool upper_bound;
+        if (!is_full() && lra.get_equality_and_right_side_for_term_on_current_x(i, rs, ci, upper_bound)) {
+            add_term(t, rs, ci, upper_bound);
+        }
+    }
+
+
+    bool hnf_cutter::hnf_has_var_with_non_integral_value() const {
+        for (unsigned j : vars())
+            if (!lia.get_value(j).is_int())
+                return true;
+        return false;
+    }
+
+    bool hnf_cutter::init_terms_for_hnf_cut() {
+        clear();
+        for (unsigned i = 0; i < lra.terms().size() && !is_full(); i++) {
+            try_add_term_to_A_for_hnf(i);
+        }
+        return hnf_has_var_with_non_integral_value();
+    }
+
+    lia_move hnf_cutter::make_hnf_cut() {
+        if (!init_terms_for_hnf_cut()) {
+            return lia_move::undef;
+        }
+        lia.settings().stats().m_hnf_cutter_calls++;
+        TRACE("hnf_cut", tout << "settings().stats().m_hnf_cutter_calls = " << lia.settings().stats().m_hnf_cutter_calls << "\n";
+              for (unsigned i : constraints_for_explanation()) {
+                  lra.constraints().display(tout, i);
+              }              
+              tout << lra.constraints();
+              );
+#ifdef Z3DEBUG
+        vector<mpq> x0 = transform_to_local_columns(lra.m_mpq_lar_core_solver.m_r_x);
+#else
+        vector<mpq> x0;
+#endif
+        lia_move r =  create_cut(lia.m_t, lia.m_k, lia.m_ex, lia.m_upper, x0);
+
+        if (r == lia_move::cut) {      
+            TRACE("hnf_cut",
+                  lra.print_term(lia.m_t, tout << "cut:"); 
+                  tout << " <= " << lia.m_k << std::endl;
+                  for (unsigned i : constraints_for_explanation()) {
+                      lra.constraints().display(tout, i);
+                  }              
+                  );
+            lp_assert(lia.current_solution_is_inf_on_cut());
+            lia.settings().stats().m_hnf_cuts++;
+            lia.m_ex->clear();        
+            for (unsigned i : constraints_for_explanation()) {
+                lia.m_ex->push_justification(i);
+            }
+        } 
+        return r;
+    }
+
+}

src/math/lp/hnf_cutter.h

@@ -3,19 +3,17 @@ Copyright (c) 2017 Microsoft Corporation
 
 Module Name:
 
-    <name>
+    hnf_cutter.h
 
 Abstract:
 
-    <abstract>
+    Cuts (branches) from Hermite matrices
 
 Author:
     Lev Nachmanson (levnach)
 
-Revision History:
-
-
 --*/
+
 #pragma once
 #include "math/lp/lar_term.h"
 #include "math/lp/hnf.h"
@@ -25,212 +23,64 @@ Revision History:
 #include "math/lp/explanation.h"
 
 namespace lp  {
+class int_solver;
+class lar_solver;
+
 class hnf_cutter {
+    int_solver&                lia;
+    lar_solver&                lra;
+    lp_settings &              m_settings;
     general_matrix             m_A;
     vector<const lar_term*>    m_terms;
     vector<bool>               m_terms_upper;
     svector<constraint_index>  m_constraints_for_explanation;
     vector<mpq>                m_right_sides;
-    lp_settings &              m_settings;
     mpq                        m_abs_max;
     bool                       m_overflow;
     var_register               m_var_register;
+
 public:
 
+    
+    hnf_cutter(int_solver& lia);
+
+    lia_move hnf_cutter::make_hnf_cut();
+
+    bool hnf_cutter::init_terms_for_hnf_cut();
+    bool hnf_cutter::hnf_has_var_with_non_integral_value() const;
+    void hnf_cutter::try_add_term_to_A_for_hnf(unsigned i);
+
+    unsigned terms_count() const { return m_terms.size();  }
     const mpq & abs_max() const { return m_abs_max; }
-    
-    hnf_cutter(lp_settings & settings) : m_settings(settings),
-                                         m_abs_max(zero_of_type<mpq>()),
-                                         m_var_register(0) {}
-
-    unsigned terms_count() const {
-        return m_terms.size();
-    }
-
     const vector<const lar_term*>& terms() const { return m_terms; }
-    const svector<unsigned>& constraints_for_explanation() const {
+    const svector<unsigned>& constraints_for_explanation() const { return m_constraints_for_explanation; }
-        return m_constraints_for_explanation;
-    }
     const vector<mpq> & right_sides() const { return m_right_sides; }
-    void clear() {
-        // m_A will be filled from scratch in init_matrix_A
-        m_var_register.clear();
-        m_terms.clear();
-        m_terms_upper.clear();
-        m_constraints_for_explanation.clear();
-        m_right_sides.clear();
-        m_abs_max = zero_of_type<mpq>();
-        m_overflow = false;
-    }
-    void add_term(const lar_term* t, const mpq &rs, constraint_index ci, bool upper_bound) {
-        m_terms.push_back(t);
-        m_terms_upper.push_back(upper_bound);
-        if (upper_bound)
-            m_right_sides.push_back(rs);
-        else
-            m_right_sides.push_back(-rs);
 
-        m_constraints_for_explanation.push_back(ci);
+    bool is_full() const;
 
-        for (const auto &p : *t) {
+    void clear();
-            m_var_register.add_var(p.var());
+    void add_term(const lar_term* t, const mpq &rs, constraint_index ci, bool upper_bound);
-            mpq t = abs(ceil(p.coeff()));
-            if (t > m_abs_max)
-                m_abs_max = t;
-        }
-    }
     
-    void print(std::ostream & out) {
+    void print(std::ostream & out);
-        out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl;
-    }
 
-    void initialize_row(unsigned i) {
+    void initialize_row(unsigned i);
-        mpq sign = m_terms_upper[i]? one_of_type<mpq>(): - one_of_type<mpq>();
+    void init_matrix_A();
-        m_A.init_row_from_container(i, * m_terms[i], [this](unsigned j) { return m_var_register.add_var(j);}, sign);
+    void find_h_minus_1_b(const general_matrix& H, vector<mpq> & b);
-    }
 
-    void init_matrix_A() {
+    vector<mpq> create_b(const svector<unsigned> & basis_rows);
-        m_A = general_matrix(terms_count(), vars().size());
-        for (unsigned i = 0; i < terms_count(); i++)
-            initialize_row(i);
-    }
 
-    // todo: as we need only one row i with non integral b[i] need to optimize later
+    int find_cut_row_index(const vector<mpq> & b);
-    void find_h_minus_1_b(const general_matrix& H, vector<mpq> & b) {
-        // the solution will be put into b
-        for (unsigned i = 0; i < H.row_count() ;i++) {
-            for (unsigned j = 0; j < i; j++) {
-                b[i] -= H[i][j]*b[j];
-            }
-            b[i] /= H[i][i];
-            // consider return from here if b[i] is not an integer and return i
-        }
-    }
-
-    vector<mpq> create_b(const svector<unsigned> & basis_rows) {
-        if (basis_rows.size() == m_right_sides.size())
-            return m_right_sides;
-        vector<mpq> b;
-        for (unsigned i : basis_rows) {
-            b.push_back(m_right_sides[i]);
-        }
-        return b;
-    }
-
-    int find_cut_row_index(const vector<mpq> & b) {
-        int ret = -1;
-        int n = 0;
-        for (int i = 0; i < static_cast<int>(b.size()); i++) {
-            if (is_integer(b[i])) continue;
-            if (n == 0 ) {
-                lp_assert(ret == -1);
-                n = 1;
-                ret = i;
-            } else {
-                if (m_settings.random_next() % (++n) == 0) {
-                    ret = i;
-                }
-            }
-        }
-        return ret;
-    }
 
     // fills e_i*H_minus_1
-    void get_ei_H_minus_1(unsigned i, const general_matrix& H, vector<mpq> & row) {
+    void get_ei_H_minus_1(unsigned i, const general_matrix& H, vector<mpq> & row);
-        // we solve x = ei * H_min_1
-        // or x * H = ei
-        unsigned m = H.row_count();
-        for (unsigned k = i + 1; k < m; k++) {
-            row[k] = zero_of_type<mpq>();
-        }
-        row[i] = one_of_type<mpq>() / H[i][i];
-        for(int k = i - 1; k >= 0; k--) {
-            mpq t = zero_of_type<mpq>();
-            for (unsigned l = k + 1; l <= i; l++) {
-                t += H[l][k]*row[l];
-            }
-            row[k] = -t / H[k][k]; 
-        }
 
-        // // test region
+    void fill_term(const vector<mpq> & row, lar_term& t);
-        // vector<mpq> ei(H.row_count(), zero_of_type<mpq>());
-        // ei[i] = one_of_type<mpq>();
-        // vector<mpq> pr = row * H;
-        // pr.shrink(ei.size());
-        // lp_assert(ei == pr);
-        // // end test region
-        
-    }
-
-    void fill_term(const vector<mpq> & row, lar_term& t) {
-        for (unsigned j = 0; j < row.size(); j++) {
-            if (!is_zero(row[j]))
-                t.add_monomial(row[j], m_var_register.local_to_external(j));
-        }
-    }
 #ifdef Z3DEBUG
-    vector<mpq> transform_to_local_columns(const vector<impq> & x) const {
+    vector<mpq> transform_to_local_columns(const vector<impq> & x) const;
-        vector<mpq> ret;
-        for (unsigned j = 0; j < vars().size(); j++) {
-             ret.push_back(x[m_var_register.local_to_external(j)].x);
-        }
-        return ret;
-    }
 #endif
-    void shrink_explanation(const svector<unsigned>& basis_rows) {
+    void shrink_explanation(const svector<unsigned>& basis_rows);
-        svector<unsigned> new_expl;
+    bool overflow() const;    
-        for (unsigned i : basis_rows) {
+    lia_move create_cut(lar_term& t, mpq& k, explanation* ex, bool & upper, const vector<mpq> & x0);
-            new_expl.push_back(m_constraints_for_explanation[i]);
+    svector<unsigned> vars() const;
-        }
-        m_constraints_for_explanation = new_expl;
-    }
-
-    bool overflow() const { return m_overflow; }
-    
-    lia_move create_cut(lar_term& t, mpq& k, explanation* ex, bool & upper, const vector<mpq> & x0) {
-        // we suppose that x0 has at least one non integer element 
-        (void)x0;
-
-        init_matrix_A();
-        svector<unsigned> basis_rows;
-        mpq big_number = m_abs_max.expt(3);
-        mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows, big_number);
-        
-        if (d >= big_number) {
-            return lia_move::undef;
-        }
-        
-        if (m_settings.get_cancel_flag()) {
-            return lia_move::undef;
-        }
-
-        if (basis_rows.size() < m_A.row_count()) {
-            m_A.shrink_to_rank(basis_rows);
-            shrink_explanation(basis_rows);
-        }
-        
-        hnf<general_matrix> h(m_A, d);        
-        vector<mpq> b = create_b(basis_rows);
-        lp_assert(m_A * x0 == b);
-        find_h_minus_1_b(h.W(), b);
-        int cut_row = find_cut_row_index(b);
-
-        if (cut_row == -1) {
-            return lia_move::undef;
-        }
-
-        // the matrix is not square - we can get
-        // all integers in b's projection
-        
-        vector<mpq> row(m_A.column_count());
-        get_ei_H_minus_1(cut_row, h.W(), row);
-        vector<mpq> f = row * m_A;
-        fill_term(f, t);
-        k = floor(b[cut_row]);
-        upper = true;
-        return lia_move::cut;
-    }
-
-    svector<unsigned> vars() const { return m_var_register.vars(); }
 };
 }

src/math/lp/int_solver.cpp

@@ -18,7 +18,7 @@ namespace lp {
 int_solver::int_solver(lar_solver& lar_slv) :
     lra(lar_slv),
     m_number_of_calls(0),
-    m_hnf_cutter(settings()),
+    m_hnf_cutter(*this),
     m_hnf_cut_period(settings().hnf_cut_period()) {
     lra.set_int_solver(this);
 }
@@ -172,80 +172,12 @@ bool int_solver::should_gomory_cut() {
     return m_number_of_calls % settings().m_int_gomory_cut_period == 0;
 }
 
-void int_solver::try_add_term_to_A_for_hnf(unsigned i) {
-    mpq rs;
-    const lar_term* t = lra.terms()[i];
-    constraint_index ci;
-    bool upper_bound;
-    if (!hnf_cutter_is_full() && lra.get_equality_and_right_side_for_term_on_current_x(i, rs, ci, upper_bound)) {
-        m_hnf_cutter.add_term(t, rs, ci, upper_bound);
-    }
-}
-
-bool int_solver::hnf_cutter_is_full() const {
-    return
-        m_hnf_cutter.terms_count() >= settings().limit_on_rows_for_hnf_cutter
-                                     ||
-        m_hnf_cutter.vars().size() >= settings().limit_on_columns_for_hnf_cutter;
-}
-
-bool int_solver::hnf_has_var_with_non_integral_value() const {
-    for (unsigned j : m_hnf_cutter.vars())
-        if (!get_value(j).is_int())
-            return true;
-    return false;
-}
-
-bool int_solver::init_terms_for_hnf_cut() {
-    m_hnf_cutter.clear();
-    for (unsigned i = 0; i < lra.terms().size() && !hnf_cutter_is_full(); i++) {
-        try_add_term_to_A_for_hnf(i);
-    }
-    return hnf_has_var_with_non_integral_value();
-}
-
-lia_move int_solver::make_hnf_cut() {
-    if (!init_terms_for_hnf_cut()) {
-        return lia_move::undef;
-    }
-    settings().stats().m_hnf_cutter_calls++;
-    TRACE("hnf_cut", tout << "settings().stats().m_hnf_cutter_calls = " << settings().stats().m_hnf_cutter_calls << "\n";
-          for (unsigned i : m_hnf_cutter.constraints_for_explanation()) {
-              lra.constraints().display(tout, i);
-          }              
-          tout << lra.constraints();
-          );
-#ifdef Z3DEBUG
-    vector<mpq> x0 = m_hnf_cutter.transform_to_local_columns(lra.m_mpq_lar_core_solver.m_r_x);
-#else
-    vector<mpq> x0;
-#endif
-    lia_move r =  m_hnf_cutter.create_cut(m_t, m_k, m_ex, m_upper, x0);
-
-    if (r == lia_move::cut) {      
-        TRACE("hnf_cut",
-              lra.print_term(m_t, tout << "cut:"); 
-              tout << " <= " << m_k << std::endl;
-              for (unsigned i : m_hnf_cutter.constraints_for_explanation()) {
-                  lra.constraints().display(tout, i);
-              }              
-              );
-        lp_assert(current_solution_is_inf_on_cut());
-        settings().stats().m_hnf_cuts++;
-        m_ex->clear();        
-        for (unsigned i : m_hnf_cutter.constraints_for_explanation()) {
-             m_ex->push_justification(i);
-        }
-    } 
-    return r;
-}
-
 bool int_solver::should_hnf_cut() {
     return settings().m_enable_hnf && m_number_of_calls % m_hnf_cut_period == 0;
 }
 
 lia_move int_solver::hnf_cut() {
-    lia_move r = make_hnf_cut();
+    lia_move r = m_hnf_cutter.make_hnf_cut();
     if (r == lia_move::undef) {
         m_hnf_cut_period *= 2;
     }

src/math/lp/int_solver.h

@@ -36,6 +36,7 @@ class int_solver {
     friend class int_cube;
     friend class int_branch;
     friend class int_gcd_test;
+    friend class hnf_cutter;
 
     lar_solver&         lra;
     unsigned            m_number_of_calls;
@@ -55,7 +56,6 @@ public:
     lar_term const& get_term() const { return m_t; }
     mpq const& get_offset() const { return m_k; }
     bool is_upper() const { return m_upper; }
-    //lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex);    
     bool is_base(unsigned j) const;
     bool is_real(unsigned j) const;
     const impq & lower_bound(unsigned j) const;
@@ -106,12 +106,6 @@ public:
     bool all_columns_are_bounded() const;
     void find_feasible_solution();
     lia_move hnf_cut();
-    lia_move make_hnf_cut();
-    bool init_terms_for_hnf_cut();
-    bool hnf_matrix_is_empty() const;
-    void try_add_term_to_A_for_hnf(unsigned term_index);
-    bool hnf_has_var_with_non_integral_value() const;
-    bool hnf_cutter_is_full() const;
     void patch_nbasic_column(unsigned j);
   };
 }