1 /*++
2 Copyright (c) 2017 Microsoft Corporation
3 
4 Module Name:
5 
6     hnf_cutter.cpp
7 
8 Author:
9     Lev Nachmanson (levnach)
10 
11 --*/
12 #include "math/lp/int_solver.h"
13 #include "math/lp/lar_solver.h"
14 #include "math/lp/hnf_cutter.h"
15 
16 namespace lp  {
17 
hnf_cutter(int_solver & lia)18     hnf_cutter::hnf_cutter(int_solver& lia):
19         lia(lia),
20         lra(lia.lra),
21         m_settings(lia.settings()),
22         m_abs_max(zero_of_type<mpq>()),
23         m_var_register(0) {}
24 
is_full() const25     bool hnf_cutter::is_full() const {
26         return
27             terms_count() >= lia.settings().limit_on_rows_for_hnf_cutter ||
28             vars().size() >= lia.settings().limit_on_columns_for_hnf_cutter;
29     }
30 
clear()31     void hnf_cutter::clear() {
32         // m_A will be filled from scratch in init_matrix_A
33         m_var_register.clear();
34         m_terms.clear();
35         m_terms_upper.clear();
36         m_constraints_for_explanation.clear();
37         m_right_sides.clear();
38         m_abs_max = zero_of_type<mpq>();
39         m_overflow = false;
40     }
41 
add_term(const lar_term * t,const mpq & rs,constraint_index ci,bool upper_bound)42     void hnf_cutter::add_term(const lar_term* t, const mpq &rs, constraint_index ci, bool upper_bound) {
43         m_terms.push_back(t);
44         m_terms_upper.push_back(upper_bound);
45         if (upper_bound)
46             m_right_sides.push_back(rs);
47         else
48             m_right_sides.push_back(-rs);
49 
50         m_constraints_for_explanation.push_back(ci);
51 
52         for (lar_term::ival p : *t) {
53             m_var_register.add_var(p.column().index(), true); // hnf only deals with integral variables for now
54             mpq t = abs(ceil(p.coeff()));
55             if (t > m_abs_max)
56                 m_abs_max = t;
57         }
58     }
59 
print(std::ostream & out)60     void hnf_cutter::print(std::ostream & out) {
61         out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl;
62     }
63 
initialize_row(unsigned i)64     void hnf_cutter::initialize_row(unsigned i) {
65         mpq sign = m_terms_upper[i]? one_of_type<mpq>(): - one_of_type<mpq>();
66         m_A.init_row_from_container(i, * m_terms[i], [this](unsigned j) { return m_var_register.add_var(j, true);}, sign);// hnf only deals with integral variables for now
67     }
68 
init_matrix_A()69     void hnf_cutter::init_matrix_A() {
70         m_A = general_matrix(terms_count(), vars().size());
71         for (unsigned i = 0; i < terms_count(); i++)
72             initialize_row(i);
73     }
74 
75     // todo: as we need only one row i with non integral b[i] need to optimize later
find_h_minus_1_b(const general_matrix & H,vector<mpq> & b)76     void hnf_cutter::find_h_minus_1_b(const general_matrix& H, vector<mpq> & b) {
77         // the solution will be put into b
78         for (unsigned i = 0; i < H.row_count() ;i++) {
79             for (unsigned j = 0; j < i; j++) {
80                 b[i] -= H[i][j]*b[j];
81             }
82             b[i] /= H[i][i];
83             // consider return from here if b[i] is not an integer and return i
84         }
85     }
86 
create_b(const svector<unsigned> & basis_rows)87     vector<mpq> hnf_cutter::create_b(const svector<unsigned> & basis_rows) {
88         if (basis_rows.size() == m_right_sides.size())
89             return m_right_sides;
90         vector<mpq> b;
91         for (unsigned i : basis_rows) {
92             b.push_back(m_right_sides[i]);
93         }
94         return b;
95     }
96 
find_cut_row_index(const vector<mpq> & b)97     int hnf_cutter::find_cut_row_index(const vector<mpq> & b) {
98         int ret = -1;
99         int n = 0;
100         for (int i = 0; i < static_cast<int>(b.size()); i++) {
101             if (is_integer(b[i])) continue;
102             if (n == 0 ) {
103                 lp_assert(ret == -1);
104                 n = 1;
105                 ret = i;
106             } else {
107                 if (m_settings.random_next() % (++n) == 0) {
108                     ret = i;
109                 }
110             }
111         }
112         return ret;
113     }
114 
115     // fills e_i*H_minus_1
get_ei_H_minus_1(unsigned i,const general_matrix & H,vector<mpq> & row)116     void hnf_cutter::get_ei_H_minus_1(unsigned i, const general_matrix& H, vector<mpq> & row) {
117         // we solve x = ei * H_min_1
118         // or x * H = ei
119         unsigned m = H.row_count();
120         for (unsigned k = i + 1; k < m; k++) {
121             row[k] = zero_of_type<mpq>();
122         }
123         row[i] = one_of_type<mpq>() / H[i][i];
124         for(int k = i - 1; k >= 0; k--) {
125             mpq t = zero_of_type<mpq>();
126             for (unsigned l = k + 1; l <= i; l++) {
127                 t += H[l][k]*row[l];
128             }
129             row[k] = -t / H[k][k];
130         }
131     }
132 
fill_term(const vector<mpq> & row,lar_term & t)133     void hnf_cutter::fill_term(const vector<mpq> & row, lar_term& t) {
134         for (unsigned j = 0; j < row.size(); j++) {
135             if (!is_zero(row[j]))
136                 t.add_monomial(row[j], m_var_register.local_to_external(j));
137         }
138     }
139 #ifdef Z3DEBUG
transform_to_local_columns(const vector<impq> & x) const140     vector<mpq> hnf_cutter::transform_to_local_columns(const vector<impq> & x) const {
141         vector<mpq> ret;
142         for (unsigned j = 0; j < vars().size(); j++) {
143              ret.push_back(x[m_var_register.local_to_external(j)].x);
144         }
145         return ret;
146     }
147 #endif
shrink_explanation(const svector<unsigned> & basis_rows)148     void hnf_cutter::shrink_explanation(const svector<unsigned>& basis_rows) {
149         svector<unsigned> new_expl;
150         for (unsigned i : basis_rows) {
151             new_expl.push_back(m_constraints_for_explanation[i]);
152         }
153         m_constraints_for_explanation = new_expl;
154     }
155 
overflow() const156     bool hnf_cutter::overflow() const { return m_overflow; }
157 /*
158   Here is the citation from "Cutting the Mix" by Jürgen Christ and Jochen Hoenicke.
159 
160   The algorithm is based on the Simplex algorithm. The solution space
161   forms a polyhedron in Q^n . If the solution space is non-empty, the
162   Simplex algorithm returns a solution of Ax <= b . We further assume
163   that the returned solution x0 is a vertex of the polyhedron, i. e.,
164   there is a nonsingular square submatrix A′ and a corresponding
165   vector b′ , such that A′x0=b′ . We call A′x <=b′ the defining
166   constraints of the vertex. If the returned solution is not on a
167   vertex we introduce artificial branches on input variables into A
168   and use these branches as defining constraints. These branches are
169   rarely needed in practise.
170 
171 The main idea is to bring the constraint system A′x<=b′ into a Hermite
172 normal form H and to compute the unimodular matrix U with A′U=H . The
173 Hermite normal form is uniquely defined. The constraint system A′x<=b′
174 is equivalent to Hy <=b′ with y:=(U−1)x . Since the solution x0 of
175 A′x0=b′ is not integral, the corresponding vector y0=(U−1)x0 is not
176 integral, either. The cuts from proofs algorithm creates an extended
177 branch on one of the components y_i of y , i. e., y_i <= floor(y0_i) or
178 y_i>=ceil(y0_i).  Further on in the paper there is a lemma showing that
179 branch y_i >= ceil(y0_i) is impossible.
180  */
create_cut(lar_term & t,mpq & k,explanation * ex,bool & upper,const vector<mpq> & x0)181     lia_move hnf_cutter::create_cut(lar_term& t, mpq& k, explanation* ex, bool & upper
182 #ifdef Z3DEBUG
183                                     , const vector<mpq> & x0
184                                     // we suppose that x0 has at least one non integer element
185 #endif
186                                     ) {
187         init_matrix_A();
188         svector<unsigned> basis_rows;
189         mpq big_number = m_abs_max.expt(3);
190         mpq d = hnf_calc::determinant_of_rectangular_matrix(m_A, basis_rows, big_number);
191 
192         if (d >= big_number) {
193             return lia_move::undef;
194         }
195 
196         if (m_settings.get_cancel_flag()) {
197             return lia_move::undef;
198         }
199 
200         if (basis_rows.size() < m_A.row_count()) {
201             m_A.shrink_to_rank(basis_rows);
202             shrink_explanation(basis_rows);
203         }
204 
205         hnf<general_matrix> h(m_A, d);
206         vector<mpq> b = create_b(basis_rows);
207 #ifdef Z3DEBUG
208         lp_assert(m_A * x0 == b);
209 #endif
210 
211         find_h_minus_1_b(h.W(), b);
212         int cut_row = find_cut_row_index(b);
213 
214         if (cut_row == -1) {
215             return lia_move::undef;
216         }
217 
218         // the matrix is not square - we can get
219         // all integers in b's projection
220 
221         vector<mpq> row(m_A.column_count());
222         get_ei_H_minus_1(cut_row, h.W(), row);
223         vector<mpq> f = row * m_A;
224         fill_term(f, t);
225         k = floor(b[cut_row]);
226         upper = true;
227         return lia_move::cut;
228     }
229 
vars() const230     svector<unsigned> hnf_cutter::vars() const { return m_var_register.vars(); }
231 
try_add_term_to_A_for_hnf(tv const & i)232     void hnf_cutter::try_add_term_to_A_for_hnf(tv const &i) {
233         mpq rs;
234         const lar_term& t = lra.get_term(i);
235         constraint_index ci;
236         bool upper_bound;
237         if (!is_full() && lra.get_equality_and_right_side_for_term_on_current_x(i, rs, ci, upper_bound)) {
238             add_term(&t, rs, ci, upper_bound);
239         }
240     }
241 
hnf_has_var_with_non_integral_value() const242     bool hnf_cutter::hnf_has_var_with_non_integral_value() const {
243         for (unsigned j : vars())
244             if (!lia.get_value(j).is_int())
245                 return true;
246         return false;
247     }
248 
init_terms_for_hnf_cut()249     bool hnf_cutter::init_terms_for_hnf_cut() {
250         clear();
251         for (unsigned i = 0; i < lra.terms().size() && !is_full(); i++) {
252             try_add_term_to_A_for_hnf(tv::term(i));
253         }
254         return hnf_has_var_with_non_integral_value();
255     }
256 
make_hnf_cut()257     lia_move hnf_cutter::make_hnf_cut() {
258         if (!init_terms_for_hnf_cut()) {
259             return lia_move::undef;
260         }
261         lia.settings().stats().m_hnf_cutter_calls++;
262         TRACE("hnf_cut", tout << "settings().stats().m_hnf_cutter_calls = " << lia.settings().stats().m_hnf_cutter_calls << "\n";
263               for (unsigned i : constraints_for_explanation()) {
264                   lra.constraints().display(tout, i);
265               }
266               tout << lra.constraints();
267               );
268 #ifdef Z3DEBUG
269         vector<mpq> x0 = transform_to_local_columns(lra.r_x());
270 #endif
271         lia_move r =  create_cut(lia.m_t, lia.m_k, lia.m_ex, lia.m_upper
272 #ifdef Z3DEBUG
273                                  , x0
274 #endif
275                                  );
276 
277         if (r == lia_move::cut) {
278             TRACE("hnf_cut",
279                   lra.print_term(lia.m_t, tout << "cut:");
280                   tout << " <= " << lia.m_k << std::endl;
281                   for (unsigned i : constraints_for_explanation()) {
282                       lra.constraints().display(tout, i);
283                   }
284                   );
285             lp_assert(lia.current_solution_is_inf_on_cut());
286             lia.settings().stats().m_hnf_cuts++;
287             lia.m_ex->clear();
288             for (unsigned i : constraints_for_explanation()) {
289                 lia.m_ex->push_back(i);
290             }
291         }
292         return r;
293     }
294 
295 }
296