1 namespace lp { 2 #include "math/lp/lp_utils.h" 3 struct gomory_test { gomory_testgomory_test4 gomory_test( 5 std::function<std::string (unsigned)> name_function_p, 6 std::function<mpq (unsigned)> get_value_p, 7 std::function<bool (unsigned)> at_low_p, 8 std::function<bool (unsigned)> at_upper_p, 9 std::function<impq (unsigned) > lower_bound_p, 10 std::function<impq (unsigned) > upper_bound_p, 11 std::function<unsigned (unsigned) > column_lower_bound_constraint_p, 12 std::function<unsigned (unsigned) > column_upper_bound_constraint_p 13 ) : 14 m_name_function(name_function_p), 15 get_value(get_value_p), 16 at_low(at_low_p), 17 at_upper(at_upper_p), 18 lower_bound(lower_bound_p), 19 upper_bound(upper_bound_p), 20 column_lower_bound_constraint(column_lower_bound_constraint_p), 21 column_upper_bound_constraint(column_upper_bound_constraint_p) 22 {} 23 24 std::function<std::string (unsigned)> m_name_function; 25 std::function<mpq (unsigned)> get_value; 26 std::function<bool (unsigned)> at_low; 27 std::function<bool (unsigned)> at_upper; 28 std::function<impq (unsigned) > lower_bound; 29 std::function<impq (unsigned) > upper_bound; 30 std::function<unsigned (unsigned) > column_lower_bound_constraint; 31 std::function<unsigned (unsigned) > column_upper_bound_constraint; 32 is_realgomory_test33 bool is_real(unsigned) { return false; } // todo: test real case print_rowgomory_test34 void print_row(std::ostream& out, vector<std::pair<mpq, unsigned>> & row ) { 35 bool first = true; 36 for (const auto & it : row) { 37 auto val = it.first; 38 if (first) { 39 first = false; 40 } else { 41 if (numeric_traits<mpq>::is_pos(val)) { 42 out << " + "; 43 } else { 44 out << " - "; 45 val = -val; 46 } 47 } 48 if (val == -numeric_traits<mpq>::one()) 49 out << " - "; 50 else if (val != numeric_traits<mpq>::one()) 51 out << T_to_string(val); 52 53 out << m_name_function(it.second); 54 } 55 } 56 real_case_in_gomory_cutgomory_test57 void real_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term& pol, explanation & expl, const mpq& f_0, const mpq& one_minus_f_0) { 58 TRACE("gomory_cut_detail_real", tout << "real\n";); 59 mpq new_a; 60 if (at_low(x_j)) { 61 if (a.is_pos()) { 62 new_a = a / (1 - f_0); 63 } 64 else { 65 new_a = a / f_0; 66 new_a.neg(); 67 } 68 k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation than 69 // k += lower_bound(x_j).x * new_a; 70 expl.add_pair(column_lower_bound_constraint(x_j), new_a); 71 } 72 else { 73 lp_assert(at_upper(x_j)); 74 if (a.is_pos()) { 75 new_a = a / f_0; 76 new_a.neg(); // the upper terms are inverted. 77 } 78 else { 79 new_a = a / (mpq(1) - f_0); 80 } 81 k.addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a; 82 expl.add_pair(column_upper_bound_constraint(x_j), new_a); 83 } 84 TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << k << "\n";); 85 pol.add_monomial(new_a, x_j); 86 } 87 int_case_in_gomory_cutgomory_test88 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) { 89 lp_assert(is_integer(x_j)); 90 lp_assert(!a.is_int()); 91 lp_assert(f_0 > zero_of_type<mpq>() && f_0 < one_of_type<mpq>()); 92 mpq f_j = fractional_part(a); 93 TRACE("gomory_cut_detail", 94 tout << a << " x_j = " << x_j << ", k = " << k << "\n"; 95 tout << "f_j: " << f_j << "\n"; 96 tout << "f_0: " << f_0 << "\n"; 97 tout << "1 - f_0: " << one_minus_f_0 << "\n"; 98 tout << "at_low(" << x_j << ") = " << at_low(x_j) << std::endl; 99 ); 100 lp_assert (!f_j.is_zero()); 101 mpq new_a; 102 if (at_low(x_j)) { 103 if (f_j <= one_minus_f_0) { 104 new_a = f_j / one_minus_f_0; 105 } 106 else { 107 new_a = (1 - f_j) / f_0; 108 } 109 k.addmul(new_a, lower_bound(x_j).x); 110 expl.add_pair(column_lower_bound_constraint(x_j), new_a); 111 } 112 else { 113 lp_assert(at_upper(x_j)); 114 if (f_j <= f_0) { 115 new_a = f_j / f_0; 116 } 117 else { 118 new_a = (mpq(1) - f_j) / (one_minus_f_0); 119 } 120 new_a.neg(); // the upper terms are inverted 121 k.addmul(new_a, upper_bound(x_j).x); 122 expl.add_pair(column_upper_bound_constraint(x_j), new_a); 123 } 124 TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << k << "\n";); 125 t.add_monomial(new_a, x_j); 126 lcm_den = lcm(lcm_den, denominator(new_a)); 127 } 128 129 report_conflict_from_gomory_cutgomory_test130 void report_conflict_from_gomory_cut(mpq &k) { 131 lp_assert(false); 132 } 133 adjust_term_and_k_for_some_ints_case_gomorygomory_test134 void adjust_term_and_k_for_some_ints_case_gomory(lar_term& t, mpq& k, mpq &lcm_den) { 135 lp_assert(!t.is_empty()); 136 auto pol = t.coeffs_as_vector(); 137 t.clear(); 138 if (pol.size() == 1) { 139 TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;); 140 unsigned v = pol[0].second; 141 lp_assert(is_integer(v)); 142 const mpq& a = pol[0].first; 143 k /= a; 144 if (a.is_pos()) { // we have av >= k 145 if (!k.is_int()) 146 k = ceil(k); 147 // switch size 148 t.add_monomial(- mpq(1), v); 149 k.neg(); 150 } else { 151 if (!k.is_int()) 152 k = floor(k); 153 t.add_monomial(mpq(1), v); 154 } 155 } else { 156 TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;); 157 lcm_den = lcm(lcm_den, denominator(k)); 158 TRACE("gomory_cut_detail", tout << "k: " << k << " lcm_den: " << lcm_den << "\n"; 159 for (unsigned i = 0; i < pol.size(); i++) { 160 tout << pol[i].first << " " << pol[i].second << "\n"; 161 } 162 tout << "k: " << k << "\n";); 163 lp_assert(lcm_den.is_pos()); 164 if (!lcm_den.is_one()) { 165 // normalize coefficients of integer parameters to be integers. 166 for (auto & pi: pol) { 167 pi.first *= lcm_den; 168 SASSERT(!is_integer(pi.second) || pi.first.is_int()); 169 } 170 k *= lcm_den; 171 } 172 TRACE("gomory_cut_detail", tout << "after *lcm\n"; 173 for (unsigned i = 0; i < pol.size(); i++) { 174 tout << pol[i].first << " * v" << pol[i].second << "\n"; 175 } 176 tout << "k: " << k << "\n";); 177 178 // negate everything to return -pol <= -k 179 for (const auto & pi: pol) 180 t.add_monomial(-pi.first, pi.second); 181 k.neg(); 182 } 183 TRACE("gomory_cut_detail", tout << "k = " << k << std::endl;); 184 lp_assert(k.is_int()); 185 } 186 print_termgomory_test187 void print_term(lar_term & t, std::ostream & out) { 188 vector<std::pair<mpq, unsigned>> row; 189 for (auto p : t) 190 row.push_back(std::make_pair(p.coeff(), p.column().index())); 191 print_row(out, row); 192 } 193 mk_gomory_cutgomory_test194 void mk_gomory_cut(lar_term& t, mpq& k, explanation & expl, unsigned inf_col, vector<std::pair<mpq, unsigned>> & row) { 195 enable_trace("gomory_cut"); 196 enable_trace("gomory_cut_detail"); 197 198 TRACE("gomory_cut", 199 tout << "applying cut at:\n"; print_row(tout, row); 200 tout << std::endl << "inf_col = " << inf_col << std::endl; 201 ); 202 203 // gomory will be t >= k 204 k = 1; 205 mpq lcm_den(1); 206 unsigned x_j; 207 mpq a; 208 bool some_int_columns = false; 209 mpq f_0 = fractional_part(get_value(inf_col)); 210 mpq one_min_f_0 = 1 - f_0; 211 for ( auto pp : row) { 212 a = pp.first; 213 x_j = pp.second; 214 if (x_j == inf_col) 215 continue; 216 // make the format compatible with the format used in: Integrating Simplex with DPLL(T) 217 a.neg(); 218 if (is_real(x_j)) 219 real_case_in_gomory_cut(a, x_j, k, t, expl, f_0, one_min_f_0); 220 else { 221 if (a.is_int()) continue; // f_j will be zero and no monomial will be added 222 some_int_columns = true; 223 int_case_in_gomory_cut(a, x_j, k, t, expl, lcm_den, f_0, one_min_f_0); 224 } 225 } 226 227 if (t.is_empty()) 228 return report_conflict_from_gomory_cut(k); 229 if (some_int_columns) 230 adjust_term_and_k_for_some_ints_case_gomory(t, k, lcm_den); 231 232 TRACE("gomory_cut", tout<<"new cut :"; print_term(t, tout); tout << " >= " << k << std::endl;); 233 234 } 235 }; 236 } 237