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