1
2 /*++
3 Copyright (c) 2015 Microsoft Corporation
4
5 --*/
6
7 #include "math/simplex/sparse_matrix_def.h"
8 #include "math/simplex/simplex.h"
9 #include "math/simplex/simplex_def.h"
10 #include "util/mpq_inf.h"
11 #include "util/vector.h"
12 #include "util/rational.h"
13 #include "util/rlimit.h"
14
15 #define R rational
16 typedef simplex::simplex<simplex::mpz_ext> Simplex;
17 typedef simplex::sparse_matrix<simplex::mpz_ext> sparse_matrix;
18
vec(int i,int j)19 static vector<R> vec(int i, int j) {
20 vector<R> nv;
21 nv.resize(2);
22 nv[0] = R(i);
23 nv[1] = R(j);
24 return nv;
25 }
26
27 // static vector<R> vec(int i, int j, int k) {
28 // vector<R> nv = vec(i, j);
29 // nv.push_back(R(k));
30 // return nv;
31 // }
32
33 // static vector<R> vec(int i, int j, int k, int l) {
34 // vector<R> nv = vec(i, j, k);
35 // nv.push_back(R(l));
36 // return nv;
37 // }
38
39 /// static vector<R> vec(int i, int j, int k, int l, int x) {
40 /// vector<R> nv = vec(i, j, k, l);
41 /// nv.push_back(R(x));
42 /// return nv;
43 /// }
44
45 // static vector<R> vec(int i, int j, int k, int l, int x, int y) {
46 // vector<R> nv = vec(i, j, k, l, x);
47 // nv.push_back(R(y));
48 // return nv;
49 // }
50
51 // static vector<R> vec(int i, int j, int k, int l, int x, int y, int z) {
52 // vector<R> nv = vec(i, j, k, l, x, y);
53 // nv.push_back(R(z));
54 // return nv;
55 // }
56
57
58
add_row(Simplex & S,vector<R> const & _v,R const & _b,bool is_eq=false)59 void add_row(Simplex& S, vector<R> const& _v, R const& _b, bool is_eq = false) {
60 unsynch_mpz_manager m;
61 unsigned_vector vars;
62 vector<R> v(_v);
63 R b(_b);
64 R l(denominator(b));
65 scoped_mpz_vector coeffs(m);
66 for (unsigned i = 0; i < v.size(); ++i) {
67 l = lcm(l, denominator(v[i]));
68 vars.push_back(i);
69 S.ensure_var(i);
70 }
71 b *= l;
72 b.neg();
73 for (unsigned i = 0; i < v.size(); ++i) {
74 v[i] *= l;
75 coeffs.push_back(v[i].to_mpq().numerator());
76 }
77 unsigned nv = S.get_num_vars();
78 vars.push_back(nv);
79 vars.push_back(nv+1);
80 S.ensure_var(nv);
81 S.ensure_var(nv+1);
82 coeffs.push_back(mpz(-1));
83 coeffs.push_back(b.to_mpq().numerator());
84 mpq_inf one(mpq(1),mpq(0));
85 mpq_inf zero(mpq(0),mpq(0));
86 ENSURE(vars.size() == coeffs.size());
87 S.set_lower(nv, zero);
88 if (is_eq) S.set_upper(nv, zero);
89 S.set_lower(nv+1, one);
90 S.set_upper(nv+1, one);
91 S.add_row(nv, coeffs.size(), vars.data(), coeffs.data());
92 }
93
feas(Simplex & S)94 static void feas(Simplex& S) {
95 S.display(std::cout);
96 lbool is_sat = S.make_feasible();
97 std::cout << "feasible: " << is_sat << "\n";
98 S.display(std::cout);
99 }
100
test1()101 static void test1() {
102 reslimit rl;
103 Simplex S(rl);
104 add_row(S, vec(1,0), R(1));
105 add_row(S, vec(0,1), R(1));
106 add_row(S, vec(1,1), R(1));
107 feas(S);
108 }
109
test2()110 static void test2() {
111 reslimit rl; Simplex S(rl);
112 add_row(S, vec(1, 0), R(1));
113 add_row(S, vec(0, 1), R(1));
114 add_row(S, vec(1, 1), R(1), true);
115 feas(S);
116 }
117
test3()118 static void test3() {
119 reslimit rl; Simplex S(rl);
120 add_row(S, vec(-1, 0), R(-1));
121 add_row(S, vec(0, -1), R(-1));
122 add_row(S, vec(1, 1), R(1), true);
123 feas(S);
124 }
125
test4()126 static void test4() {
127 reslimit rl; Simplex S(rl);
128 add_row(S, vec(1, 0), R(1));
129 add_row(S, vec(0, -1), R(-1));
130 add_row(S, vec(1, 1), R(1), true);
131 feas(S);
132 }
133
tst_simplex()134 void tst_simplex() {
135 reslimit rl; Simplex S(rl);
136
137 std::cout << "simplex\n";
138
139 lbool is_sat = S.make_feasible();
140 std::cout << "feasible: " << is_sat << "\n";
141
142 unsynch_mpz_manager m;
143 unsynch_mpq_inf_manager em;
144 scoped_mpz_vector coeffs(m);
145 svector<unsigned> vars;
146 for (unsigned i = 0; i < 5; ++i) {
147 S.ensure_var(i);
148 vars.push_back(i);
149 coeffs.push_back(mpz(i+1));
150 }
151
152 // Simplex::row r = S.add_row(1, coeffs.size(), vars.c_ptr(), coeffs.c_ptr());
153 is_sat = S.make_feasible();
154 std::cout << "feasible: " << is_sat << "\n";
155 S.display(std::cout);
156 _scoped_numeral<unsynch_mpq_inf_manager> num(em);
157 num = std::make_pair(mpq(1), mpq(0));
158 S.set_lower(0, num);
159 S.set_upper(0, num);
160
161 is_sat = S.make_feasible();
162 std::cout << "feasible: " << is_sat << "\n";
163 S.display(std::cout);
164
165 test1();
166 test2();
167 test3();
168 test4();
169 }
170