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