1 #include "math/simplex/model_based_opt.h"
2 #include "util/util.h"
3 #include "util/uint_set.h"
4 
5 typedef opt::model_based_opt::var var;
6 
add_ineq(opt::model_based_opt & mbo,unsigned x,int a,int k,opt::ineq_type rel)7 static void add_ineq(opt::model_based_opt& mbo, unsigned x, int a, int k, opt::ineq_type rel) {
8     vector<var> vars;
9     vars.push_back(var(x, rational(a)));
10     mbo.add_constraint(vars, rational(k), rel);
11 }
12 
add_ineq(opt::model_based_opt & mbo,unsigned x,int a,unsigned y,int b,int k,opt::ineq_type rel)13 static void add_ineq(opt::model_based_opt& mbo,
14                      unsigned x, int a,
15                      unsigned y, int b, int k,
16                      opt::ineq_type rel) {
17     vector<var> vars;
18     vars.push_back(var(x, rational(a)));
19     vars.push_back(var(y, rational(b)));
20     mbo.add_constraint(vars, rational(k), rel);
21 }
22 
23 #if 0
24 static void add_ineq(opt::model_based_opt& mbo,
25                      unsigned x, int a,
26                      unsigned y, int b,
27                      unsigned z, int c, int k,
28                      opt::ineq_type rel) {
29     vector<var> vars;
30     vars.push_back(var(x, rational(a)));
31     vars.push_back(var(y, rational(b)));
32     vars.push_back(var(z, rational(c)));
33     mbo.add_constraint(vars, rational(k), rel);
34 }
35 #endif
36 
add_random_ineq(opt::model_based_opt & mbo,random_gen & r,svector<int> const & values,unsigned max_vars,unsigned max_coeff)37 static void add_random_ineq(opt::model_based_opt& mbo,
38                             random_gen& r,
39                             svector<int>  const& values,
40                             unsigned max_vars,
41                             unsigned max_coeff) {
42     unsigned num_vars = values.size();
43     uint_set used_vars;
44     vector<var> vars;
45     int value = 0;
46     for (unsigned i = 0; i < max_vars; ++i) {
47         unsigned x = r(num_vars);
48         if (used_vars.contains(x)) {
49             continue;
50         }
51         used_vars.insert(x);
52         int coeff = r(max_coeff + 1);
53         if (coeff == 0) {
54             continue;
55         }
56         unsigned sign = r(2);
57         coeff = sign == 0 ? coeff : -coeff;
58         vars.push_back(var(x, rational(coeff)));
59         value += coeff*values[x];
60     }
61     unsigned abs_value = value < 0 ? - value : value;
62     // value + k <= 0
63     // k <= - value
64     // range for k is 2*|value|
65     // k <= - value - range
66     opt::ineq_type rel = opt::t_le;
67 
68     int coeff = 0;
69     if (r(4) == 0) {
70         rel = opt::t_eq;
71         coeff = -value;
72     }
73     else {
74         if (abs_value > 0) {
75             coeff = -value - r(2*abs_value);
76         }
77         else {
78             coeff = 0;
79         }
80         if (coeff != -value && r(3) == 0) {
81             rel = opt::t_lt;
82         }
83     }
84     mbo.add_constraint(vars, rational(coeff), rel);
85 }
86 
check_random_ineqs(random_gen & r,unsigned num_vars,unsigned max_value,unsigned num_ineqs,unsigned max_vars,unsigned max_coeff)87 static void check_random_ineqs(random_gen& r, unsigned num_vars, unsigned max_value, unsigned num_ineqs, unsigned max_vars, unsigned max_coeff) {
88     opt::model_based_opt mbo;
89     svector<int> values;
90     for (unsigned i = 0; i < num_vars; ++i) {
91         values.push_back(r(max_value + 1));
92         mbo.add_var(rational(values.back()));
93     }
94     for (unsigned i = 0; i < num_ineqs; ++i) {
95         add_random_ineq(mbo, r, values, max_vars, max_coeff);
96     }
97 
98     vector<var> vars;
99     vars.reset();
100     vars.push_back(var(0, rational(2)));
101     vars.push_back(var(1, rational(-2)));
102     mbo.set_objective(vars, rational(0));
103 
104     mbo.display(std::cout);
105     opt::inf_eps value = mbo.maximize();
106     std::cout << "optimal: " << value << "\n";
107     mbo.display(std::cout);
108     for (unsigned i = 0; i < values.size(); ++i) {
109         std::cout << i << ": " << values[i] << " -> " << mbo.get_value(i) << "\n";
110     }
111 }
112 
check_random_ineqs()113 static void check_random_ineqs() {
114     random_gen r(1);
115 
116     for (unsigned i = 0; i < 1009; ++i) {
117         check_random_ineqs(r, 4, 5, 5, 3, 6);
118     }
119 }
120 
121 // test with upper bounds
test1()122 static void test1() {
123     opt::model_based_opt mbo;
124     vector<var> vars;
125     unsigned x = mbo.add_var(rational(2));
126     unsigned y = mbo.add_var(rational(3));
127     unsigned z = mbo.add_var(rational(4));
128     unsigned u = mbo.add_var(rational(5));
129 
130     add_ineq(mbo, x, 1, y, -1, 0, opt::t_le);
131     add_ineq(mbo, x, 1, z, -1, 0, opt::t_le);
132     add_ineq(mbo, y, 1, u, -1, 0, opt::t_le);
133     add_ineq(mbo, z, 1, u, -1, 1, opt::t_le);
134     add_ineq(mbo, u, 1, -6, opt::t_le);
135 
136     vars.reset();
137     vars.push_back(var(x, rational(2)));
138     mbo.set_objective(vars, rational(0));
139 
140     opt::inf_eps value = mbo.maximize();
141     std::cout << value << "\n";
142     std::cout << "x: " << mbo.get_value(x) << "\n";
143     std::cout << "y: " << mbo.get_value(y) << "\n";
144     std::cout << "z: " << mbo.get_value(z) << "\n";
145     std::cout << "u: " << mbo.get_value(u) << "\n";
146 }
147 
148 // test with lower bounds
test2()149 static void test2() {
150     opt::model_based_opt mbo;
151     vector<var> vars;
152     unsigned x = mbo.add_var(rational(5));
153     unsigned y = mbo.add_var(rational(4));
154     unsigned z = mbo.add_var(rational(3));
155     unsigned u = mbo.add_var(rational(2));
156 
157     add_ineq(mbo, x, -1, y, 1, 0, opt::t_le);
158     add_ineq(mbo, x, -1, z, 1, 0, opt::t_le);
159     add_ineq(mbo, y, -1, u, 1, 0, opt::t_le);
160     add_ineq(mbo, z, -1, u, 1, 1, opt::t_le);
161     add_ineq(mbo, u, -1, -6, opt::t_le);
162 
163     vars.reset();
164     vars.push_back(var(x, rational(-2)));
165     mbo.set_objective(vars, rational(0));
166 
167     opt::inf_eps value = mbo.maximize();
168     std::cout << value << "\n";
169 }
170 
171 // test unbounded
test3()172 static void test3() {
173     opt::model_based_opt mbo;
174     vector<var> vars;
175     unsigned x = mbo.add_var(rational(2));
176     unsigned y = mbo.add_var(rational(3));
177     unsigned z = mbo.add_var(rational(4));
178     unsigned u = mbo.add_var(rational(5));
179 
180     add_ineq(mbo, x, 1, y, -1, 0, opt::t_le);
181     add_ineq(mbo, x, 1, z, -1, 0, opt::t_le);
182     add_ineq(mbo, y, 1, u, -1, 0, opt::t_le);
183     add_ineq(mbo, z, 1, u, -1, 1, opt::t_le);
184 
185     vars.reset();
186     vars.push_back(var(x, rational(2)));
187     mbo.set_objective(vars, rational(0));
188 
189     opt::inf_eps value = mbo.maximize();
190     std::cout << value << "\n";
191 
192 }
193 
194 // test strict
test4()195 static void test4() {
196     opt::model_based_opt mbo;
197     vector<var> vars;
198     unsigned x = mbo.add_var(rational(2));
199     unsigned y = mbo.add_var(rational(3));
200     unsigned z = mbo.add_var(rational(4));
201     unsigned u = mbo.add_var(rational(5));
202 
203     add_ineq(mbo, x, 1, y, -1, 0, opt::t_lt);
204     add_ineq(mbo, x, 1, z, -1, 0, opt::t_lt);
205     add_ineq(mbo, y, 1, u, -1, 0, opt::t_le);
206     add_ineq(mbo, z, 1, u, -1, 1, opt::t_le);
207     add_ineq(mbo, u, 1, -6, opt::t_le);
208 
209     vars.reset();
210     vars.push_back(var(x, rational(2)));
211     mbo.set_objective(vars, rational(0));
212 
213     opt::inf_eps value = mbo.maximize();
214     std::cout << value << "\n";
215     std::cout << "x: " << mbo.get_value(x) << "\n";
216     std::cout << "y: " << mbo.get_value(y) << "\n";
217     std::cout << "z: " << mbo.get_value(z) << "\n";
218     std::cout << "u: " << mbo.get_value(u) << "\n";
219 }
220 
display_project(vector<opt::model_based_opt::def> const & defs)221 static void display_project(vector<opt::model_based_opt::def> const& defs) {
222     for (auto const& d : defs) {
223         std::cout << d << "\n";
224     }
225 }
226 
test5()227 static void test5() {
228     opt::model_based_opt mbo;
229     unsigned x = mbo.add_var(rational(2));
230     unsigned y = mbo.add_var(rational(3));
231     unsigned z = mbo.add_var(rational(4));
232     unsigned u = mbo.add_var(rational(5));
233 
234     add_ineq(mbo, x, 1, y, -1, 0, opt::t_le);
235     add_ineq(mbo, x, 1, z, -1, 0, opt::t_le);
236     add_ineq(mbo, y, 1, u, -1, 0, opt::t_le);
237     add_ineq(mbo, z, 1, u, -1, 1, opt::t_le);
238 
239     unsigned vars[2] = { y, z };
240     display_project(mbo.project(1, vars, true));
241     mbo.display(std::cout);
242 
243     display_project(mbo.project(1, vars, true));
244     mbo.display(std::cout);
245 
246     display_project(mbo.project(1, vars+1, true));
247     mbo.display(std::cout);
248 
249     vector<opt::model_based_opt::row> rows;
250     mbo.get_live_rows(rows);
251 }
252 
253 
test6()254 static void test6() {
255     opt::model_based_opt mbo;
256     unsigned x0 = mbo.add_var(rational(1));
257     unsigned x = mbo.add_var(rational(2));
258     unsigned y = mbo.add_var(rational(3));
259     unsigned z = mbo.add_var(rational(4));
260     unsigned u = mbo.add_var(rational(5));
261     unsigned v = mbo.add_var(rational(6));
262     unsigned w = mbo.add_var(rational(6));
263 
264     add_ineq(mbo, x0, 1, y, -1, 0, opt::t_le);
265     add_ineq(mbo, x, 1, y, -1, 0, opt::t_le);
266     add_ineq(mbo, y, 1, u, -1, 0, opt::t_le);
267     add_ineq(mbo, y, 1, z, -1, 1, opt::t_le);
268     add_ineq(mbo, y, 1, v, -1, 1, opt::t_le);
269     add_ineq(mbo, y, 1, w, -1, 1, opt::t_le);
270 
271     mbo.display(std::cout);
272     display_project(mbo.project(1, &y, true));
273     mbo.display(std::cout);
274 }
275 
test7()276 static void test7() {
277     opt::model_based_opt mbo;
278     unsigned x0 = mbo.add_var(rational(2));
279     unsigned x = mbo.add_var(rational(1));
280     unsigned y = mbo.add_var(rational(3));
281     unsigned z = mbo.add_var(rational(4));
282     unsigned u = mbo.add_var(rational(5));
283     unsigned v = mbo.add_var(rational(6));
284     unsigned w = mbo.add_var(rational(6));
285 
286     add_ineq(mbo, x0, 1, y, -1, 0, opt::t_le);
287     add_ineq(mbo, x, 1, y, -1, 0, opt::t_lt);
288     add_ineq(mbo, y, 1, u, -1, 0, opt::t_le);
289     add_ineq(mbo, y, 1, z, -1, 1, opt::t_le);
290     add_ineq(mbo, y, 1, v, -1, 1, opt::t_le);
291     add_ineq(mbo, y, 1, w, -1, 1, opt::t_lt);
292 
293     mbo.display(std::cout);
294     display_project(mbo.project(1, &y, true));
295     mbo.display(std::cout);
296 }
297 
test8()298 static void test8() {
299     opt::model_based_opt mbo;
300     unsigned x0 = mbo.add_var(rational(2));
301     unsigned x = mbo.add_var(rational(1));
302     unsigned y = mbo.add_var(rational(3));
303     unsigned z = mbo.add_var(rational(4));
304     unsigned u = mbo.add_var(rational(5));
305     unsigned v = mbo.add_var(rational(6));
306     //    unsigned w = mbo.add_var(rational(6));
307 
308     add_ineq(mbo, x0, 1, y, -1, 0, opt::t_le);
309     add_ineq(mbo, x, 1, y, -1, 0, opt::t_lt);
310     add_ineq(mbo, y, 1, u, -1, 0, opt::t_le);
311     add_ineq(mbo, y, 1, z, -1, 1, opt::t_le);
312     add_ineq(mbo, y, 1, v, -1, 1, opt::t_le);
313 
314     mbo.display(std::cout);
315     display_project(mbo.project(1, &y, true));
316     mbo.display(std::cout);
317 }
318 
319 
test9()320 static void test9() {
321     opt::model_based_opt mbo;
322     unsigned x0 = mbo.add_var(rational(2), true);
323     unsigned x = mbo.add_var(rational(1), true);
324     unsigned y = mbo.add_var(rational(3), true);
325     unsigned z = mbo.add_var(rational(4), true);
326     unsigned u = mbo.add_var(rational(5), true);
327     unsigned v = mbo.add_var(rational(6), true);
328 
329     add_ineq(mbo, x0, 1, y, -1, 0, opt::t_le);
330     add_ineq(mbo, x, 1, y, -1, 0, opt::t_lt);
331     add_ineq(mbo, y, 1, u, -1, 0, opt::t_le);
332     add_ineq(mbo, y, 1, z, -1, 1, opt::t_le);
333     add_ineq(mbo, y, 1, v, -1, 1, opt::t_le);
334 
335     mbo.display(std::cout);
336     display_project(mbo.project(1, &y, true));
337     mbo.display(std::cout);
338 }
339 
340 
test10()341 static void test10() {
342     opt::model_based_opt mbo;
343     unsigned x0 = mbo.add_var(rational(2), true);
344     unsigned x = mbo.add_var(rational(1), true);
345     unsigned y = mbo.add_var(rational(3), true);
346     unsigned z = mbo.add_var(rational(4), true);
347     unsigned u = mbo.add_var(rational(5), true);
348     unsigned v = mbo.add_var(rational(6), true);
349 
350     add_ineq(mbo, x0, 1, y, -2, 0, opt::t_le);
351     add_ineq(mbo, x, 1, y, -2, 0, opt::t_lt);
352     add_ineq(mbo, y, 3, u, -4, 0, opt::t_le);
353     add_ineq(mbo, y, 3, z, -5, 1, opt::t_le);
354     add_ineq(mbo, y, 3, v, -6, 1, opt::t_le);
355 
356     mbo.display(std::cout);
357     display_project(mbo.project(1, &y, true));
358     mbo.display(std::cout);
359     display_project(mbo.project(1, &x0, true));
360     mbo.display(std::cout);
361 }
362 
test11()363 static void test11() {
364     {
365         opt::model_based_opt mbo;
366         unsigned s = mbo.add_var(rational(2), true);
367         unsigned a = mbo.add_var(rational(1), true);
368         unsigned t = mbo.add_var(rational(3), true);
369 
370         add_ineq(mbo, s, 1, a, -2, 0, opt::t_le); // s - 2a <= 0
371         add_ineq(mbo, a, 2, t, -1, 0, opt::t_le); // 2a - t <= 0
372 
373         mbo.display(std::cout);
374         display_project(mbo.project(1, &a, true));
375         mbo.display(std::cout);
376     }
377 
378     {
379         opt::model_based_opt mbo;
380         unsigned s = mbo.add_var(rational(3), true);
381         unsigned a = mbo.add_var(rational(2), true);
382         unsigned t = mbo.add_var(rational(4), true);
383 
384         add_ineq(mbo, s, 1, a, -2, 0, opt::t_le); // s - 2a <= 0
385         add_ineq(mbo, a, 2, t, -1, 0, opt::t_le); // 2a - t <= 0
386 
387         mbo.display(std::cout);
388         display_project(mbo.project(1, &a, true));
389         mbo.display(std::cout);
390     }
391 
392 }
393 
394 // test with mix of upper and lower bounds
395 
tst_model_based_opt()396 void tst_model_based_opt() {
397     test10();
398     check_random_ineqs();
399     test1();
400     test2();
401     test3();
402     test4();
403     test5();
404     test6();
405     test7();
406     test8();
407     test9();
408     test11();
409 }
410