1 /* Polyhedron class implementation: simplify().
2 Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
3 Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)
4
5 This file is part of the Parma Polyhedra Library (PPL).
6
7 The PPL is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 The PPL is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software Foundation,
19 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
20
21 For the most up-to-date information see the Parma Polyhedra Library
22 site: http://bugseng.com/products/ppl/ . */
23
24 #ifndef PPL_Polyhedron_simplify_templates_hh
25 #define PPL_Polyhedron_simplify_templates_hh 1
26
27 #include "Bit_Matrix_defs.hh"
28 #include "Polyhedron_defs.hh"
29 #include <cstddef>
30 #include <limits>
31
32 namespace Parma_Polyhedra_Library {
33
34 /*!
35 \return
36 The rank of \p sys.
37
38 \param sys
39 The system to simplify: it will be modified;
40
41 \param sat
42 The saturation matrix corresponding to \p sys.
43
44 \p sys may be modified by swapping some of its rows and by possibly
45 removing some of them, if they turn out to be redundant.
46
47 If \p sys is a system of constraints, then the rows of \p sat are
48 indexed by constraints and its columns are indexed by generators;
49 otherwise, if \p sys is a system of generators, then the rows of
50 \p sat are indexed by generators and its columns by constraints.
51
52 Given a system of constraints or a system of generators, this function
53 simplifies it using Gauss' elimination method (to remove redundant
54 equalities/lines), deleting redundant inequalities/rays/points and
55 making back-substitution.
56 The explanation that follows assumes that \p sys is a system of
57 constraints. For the case when \p sys is a system of generators,
58 a similar explanation can be obtain by applying duality.
59
60 The explanation relies on the notion of <EM>redundancy</EM>.
61 (See the Introduction.)
62
63 First we make some observations that can help the reader
64 in understanding the function:
65
66 Proposition: An inequality that is saturated by all the generators
67 can be transformed to an equality.
68
69 In fact, by combining any number of generators that saturate the
70 constraints, we obtain a generator that saturates the constraints too:
71 \f[
72 \langle \vect{c}, \vect{r}_1 \rangle = 0 \land
73 \langle \vect{c}, \vect{r}_2 \rangle = 0
74 \Rightarrow
75 \langle \vect{c}, (\lambda_1 \vect{r}_1 + \lambda_2 \vect{r}_2) \rangle =
76 \lambda_1 \langle \vect{c}, \vect{r}_1 \rangle
77 + \lambda_2 \langle \vect{c}, \vect{r}_2 \rangle
78 = 0,
79 \f]
80 where \f$\lambda_1, \lambda_2\f$ can be any real number.
81 */
82 template <typename Linear_System1>
83 dimension_type
simplify(Linear_System1 & sys,Bit_Matrix & sat)84 Polyhedron::simplify(Linear_System1& sys, Bit_Matrix& sat) {
85 dimension_type num_rows = sys.num_rows();
86 const dimension_type num_cols_sat = sat.num_columns();
87
88 using std::swap;
89
90 // Looking for the first inequality in `sys'.
91 dimension_type num_lines_or_equalities = 0;
92 while (num_lines_or_equalities < num_rows
93 && sys[num_lines_or_equalities].is_line_or_equality()) {
94 ++num_lines_or_equalities;
95 }
96
97 // `num_saturators[i]' will contain the number of generators
98 // that saturate the constraint `sys[i]'.
99 if (num_rows > simplify_num_saturators_size) {
100 delete [] simplify_num_saturators_p;
101 simplify_num_saturators_p = 0;
102 simplify_num_saturators_size = 0;
103 const size_t max_size
104 = std::numeric_limits<size_t>::max() / sizeof(dimension_type);
105 const size_t new_size = compute_capacity(num_rows, max_size);
106 simplify_num_saturators_p = new dimension_type[new_size];
107 simplify_num_saturators_size = new_size;
108 }
109 dimension_type* const num_saturators = simplify_num_saturators_p;
110
111 bool sys_sorted = sys.is_sorted();
112
113 // Computing the number of saturators for each inequality,
114 // possibly identifying and swapping those that happen to be
115 // equalities (see Proposition above).
116 for (dimension_type i = num_lines_or_equalities; i < num_rows; ++i) {
117 if (sat[i].empty()) {
118 // The constraint `sys_rows[i]' is saturated by all the generators.
119 // Thus, either it is already an equality or it can be transformed
120 // to an equality (see Proposition above).
121 sys.sys.rows[i].set_is_line_or_equality();
122 // Note: simple normalization already holds.
123 sys.sys.rows[i].sign_normalize();
124 // We also move it just after all the other equalities,
125 // so that system `sys_rows' keeps its partial sortedness.
126 if (i != num_lines_or_equalities) {
127 sys.sys.rows[i].m_swap(sys.sys.rows[num_lines_or_equalities]);
128 swap(sat[i], sat[num_lines_or_equalities]);
129 swap(num_saturators[i], num_saturators[num_lines_or_equalities]);
130 }
131 ++num_lines_or_equalities;
132 // `sys' is no longer sorted.
133 sys_sorted = false;
134 }
135 else {
136 // There exists a generator which does not saturate `sys[i]',
137 // so that `sys[i]' is indeed an inequality.
138 // We store the number of its saturators.
139 num_saturators[i] = num_cols_sat - sat[i].count_ones();
140 }
141 }
142
143 sys.set_sorted(sys_sorted);
144 PPL_ASSERT(sys.OK());
145
146 // At this point, all the equalities of `sys' (included those
147 // inequalities that we just transformed to equalities) have
148 // indexes between 0 and `num_lines_or_equalities' - 1,
149 // which is the property needed by method gauss().
150 // We can simplify the system of equalities, obtaining the rank
151 // of `sys' as result.
152 const dimension_type rank = sys.gauss(num_lines_or_equalities);
153
154 // Now the irredundant equalities of `sys' have indexes from 0
155 // to `rank' - 1, whereas the equalities having indexes from `rank'
156 // to `num_lines_or_equalities' - 1 are all redundant.
157 // (The inequalities in `sys' have been left untouched.)
158 // The rows containing equalities are not sorted.
159
160 if (rank < num_lines_or_equalities) {
161 // We identified some redundant equalities.
162 // Moving them at the bottom of `sys':
163 // - index `redundant' runs through the redundant equalities
164 // - index `erasing' identifies the first row that should
165 // be erased after this loop.
166 // Note that we exit the loop either because we have removed all
167 // redundant equalities or because we have moved all the
168 // inequalities.
169 for (dimension_type redundant = rank,
170 erasing = num_rows;
171 redundant < num_lines_or_equalities
172 && erasing > num_lines_or_equalities;
173 ) {
174 --erasing;
175 sys.remove_row(redundant);
176 swap(sat[redundant], sat[erasing]);
177 swap(num_saturators[redundant], num_saturators[erasing]);
178 ++redundant;
179 }
180 // Adjusting the value of `num_rows' to the number of meaningful
181 // rows of `sys': `num_lines_or_equalities' - `rank' is the number of
182 // redundant equalities moved to the bottom of `sys', which are
183 // no longer meaningful.
184 num_rows -= num_lines_or_equalities - rank;
185
186 // If the above loop exited because it moved all inequalities, it may not
187 // have removed all the rendundant rows.
188 sys.remove_trailing_rows(sys.num_rows() - num_rows);
189
190 PPL_ASSERT(sys.num_rows() == num_rows);
191
192 sat.remove_trailing_rows(num_lines_or_equalities - rank);
193
194 // Adjusting the value of `num_lines_or_equalities'.
195 num_lines_or_equalities = rank;
196 }
197
198 const dimension_type old_num_rows = sys.num_rows();
199
200 // Now we use the definition of redundancy (given in the Introduction)
201 // to remove redundant inequalities.
202
203 // First we check the saturation rule, which provides a necessary
204 // condition for an inequality to be irredundant (i.e., it provides
205 // a sufficient condition for identifying redundant inequalities).
206 // Let
207 //
208 // num_saturators[i] = num_sat_lines[i] + num_sat_rays_or_points[i],
209 // dim_lin_space = num_irredundant_lines,
210 // dim_ray_space
211 // = dim_vector_space - num_irredundant_equalities - dim_lin_space
212 // = num_columns - 1 - num_lines_or_equalities - dim_lin_space,
213 // min_sat_rays_or_points = dim_ray_space.
214 //
215 // An inequality saturated by less than `dim_ray_space' _rays/points_
216 // is redundant. Thus we have the implication
217 //
218 // (num_saturators[i] - num_sat_lines[i] < dim_ray_space)
219 // ==>
220 // redundant(sys[i]).
221 //
222 // Moreover, since every line saturates all inequalities, we also have
223 // dim_lin_space = num_sat_lines[i]
224 // so that we can rewrite the condition above as follows:
225 //
226 // (num_saturators[i] < num_columns - num_lines_or_equalities - 1)
227 // ==>
228 // redundant(sys[i]).
229 //
230 const dimension_type sys_num_columns
231 = sys.topology() == NECESSARILY_CLOSED ? sys.space_dimension() + 1
232 : sys.space_dimension() + 2;
233 const dimension_type min_saturators
234 = sys_num_columns - num_lines_or_equalities - 1;
235 for (dimension_type i = num_lines_or_equalities; i < num_rows; ) {
236 if (num_saturators[i] < min_saturators) {
237 // The inequality `sys[i]' is redundant.
238 --num_rows;
239 sys.remove_row(i);
240 swap(sat[i], sat[num_rows]);
241 swap(num_saturators[i], num_saturators[num_rows]);
242 }
243 else {
244 ++i;
245 }
246 }
247
248 // Now we check the independence rule.
249 for (dimension_type i = num_lines_or_equalities; i < num_rows; ) {
250 bool redundant = false;
251 // NOTE: in the inner loop, index `j' runs through _all_ the
252 // inequalities and we do not test if `sat[i]' is strictly
253 // contained into `sat[j]'. Experimentation has shown that this
254 // is faster than having `j' only run through the indexes greater
255 // than `i' and also doing the test `strict_subset(sat[i],
256 // sat[k])'.
257 for (dimension_type j = num_lines_or_equalities; j < num_rows; ) {
258 if (i == j) {
259 // We want to compare different rows of `sys'.
260 ++j;
261 }
262 else {
263 // Let us recall that each generator lies on a facet of the
264 // polyhedron (see the Introduction).
265 // Given two constraints `c_1' and `c_2', if there are `m'
266 // generators lying on the hyper-plane corresponding to `c_1',
267 // the same `m' generators lie on the hyper-plane
268 // corresponding to `c_2', too, and there is another one lying
269 // on the latter but not on the former, then `c_2' is more
270 // restrictive than `c_1', i.e., `c_1' is redundant.
271 bool strict_subset;
272 if (subset_or_equal(sat[j], sat[i], strict_subset)) {
273 if (strict_subset) {
274 // All the saturators of the inequality `sys[i]' are
275 // saturators of the inequality `sys[j]' too,
276 // and there exists at least one saturator of `sys[j]'
277 // which is not a saturator of `sys[i]'.
278 // It follows that inequality `sys[i]' is redundant.
279 redundant = true;
280 break;
281 }
282 else {
283 // We have `sat[j] == sat[i]'. Hence inequalities
284 // `sys[i]' and `sys[j]' are saturated by the same set of
285 // generators. Then we can remove either one of the two
286 // inequalities: we remove `sys[j]'.
287 --num_rows;
288 sys.remove_row(j);
289 PPL_ASSERT(sys.num_rows() == num_rows);
290 swap(sat[j], sat[num_rows]);
291 swap(num_saturators[j], num_saturators[num_rows]);
292 }
293 }
294 else {
295 // If we reach this point then we know that `sat[i]' does
296 // not contain (and is different from) `sat[j]', so that
297 // `sys[i]' is not made redundant by inequality `sys[j]'.
298 ++j;
299 }
300 }
301 }
302 if (redundant) {
303 // The inequality `sys[i]' is redundant.
304 --num_rows;
305 sys.remove_row(i);
306 PPL_ASSERT(sys.num_rows() == num_rows);
307 swap(sat[i], sat[num_rows]);
308 swap(num_saturators[i], num_saturators[num_rows]);
309 }
310 else {
311 // The inequality `sys[i]' is not redundant.
312 ++i;
313 }
314 }
315
316 // Here we physically remove the `sat' rows corresponding to the redundant
317 // inequalities previously removed from `sys'.
318 sat.remove_trailing_rows(old_num_rows - num_rows);
319
320 // At this point the first `num_lines_or_equalities' rows of 'sys'
321 // represent the irredundant equalities, while the remaining rows
322 // (i.e., those having indexes from `num_lines_or_equalities' to
323 // `num_rows' - 1) represent the irredundant inequalities.
324 #ifndef NDEBUG
325 // Check if the flag is set (that of the equalities is already set).
326 for (dimension_type i = num_lines_or_equalities; i < num_rows; ++i) {
327 PPL_ASSERT(sys[i].is_ray_or_point_or_inequality());
328 }
329 #endif
330
331 // Finally, since now the sub-system (of `sys') of the irredundant
332 // equalities is in triangular form, we back substitute each
333 // variables with the expression obtained considering the equalities
334 // starting from the last one.
335 sys.back_substitute(num_lines_or_equalities);
336
337 // The returned value is the number of irredundant equalities i.e.,
338 // the rank of the sub-system of `sys' containing only equalities.
339 // (See the Introduction for definition of lineality space dimension.)
340 return num_lines_or_equalities;
341 }
342
343 } // namespace Parma_Polyhedra_Library
344
345 #endif // !defined(PPL_Polyhedron_simplify_templates_hh)
346