00001 /* Polyhedron class implementation: simplify(). 00002 Copyright (C) 2001-2008 Roberto Bagnara <bagnara@cs.unipr.it> 00003 00004 This file is part of the Parma Polyhedra Library (PPL). 00005 00006 The PPL is free software; you can redistribute it and/or modify it 00007 under the terms of the GNU General Public License as published by the 00008 Free Software Foundation; either version 3 of the License, or (at your 00009 option) any later version. 00010 00011 The PPL is distributed in the hope that it will be useful, but WITHOUT 00012 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00013 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 00014 for more details. 00015 00016 You should have received a copy of the GNU General Public License 00017 along with this program; if not, write to the Free Software Foundation, 00018 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA. 00019 00020 For the most up-to-date information see the Parma Polyhedra Library 00021 site: http://www.cs.unipr.it/ppl/ . */ 00022 00023 #include <ppl-config.h> 00024 00025 #include "Linear_Row.defs.hh" 00026 #include "Linear_System.defs.hh" 00027 #include "Bit_Matrix.defs.hh" 00028 #include "Polyhedron.defs.hh" 00029 00030 namespace PPL = Parma_Polyhedra_Library; 00031 00080 PPL::dimension_type 00081 PPL::Polyhedron::simplify(Linear_System& sys, Bit_Matrix& sat) { 00082 // This method is only applied to a well-formed system `sys'. 00083 assert(sys.OK(true)); 00084 assert(sys.num_columns() >= 1); 00085 00086 dimension_type num_rows = sys.num_rows(); 00087 const dimension_type num_columns = sys.num_columns(); 00088 const dimension_type num_cols_sat = sat.num_columns(); 00089 00090 // Looking for the first inequality in `sys'. 00091 dimension_type num_lines_or_equalities = 0; 00092 while (num_lines_or_equalities < num_rows 00093 && sys[num_lines_or_equalities].is_line_or_equality()) 00094 ++num_lines_or_equalities; 00095 00096 // `num_saturators[i]' will contain the number of generators 00097 // that saturate the constraint `sys[i]'. 00098 if (num_rows > simplify_num_saturators_size) { 00099 delete [] simplify_num_saturators_p; 00100 simplify_num_saturators_p = 0; 00101 simplify_num_saturators_size = 0; 00102 size_t new_size = compute_capacity(num_rows); 00103 simplify_num_saturators_p = new dimension_type[new_size]; 00104 simplify_num_saturators_size = new_size; 00105 } 00106 dimension_type* num_saturators = simplify_num_saturators_p; 00107 00108 // Computing the number of saturators for each inequality, 00109 // possibly identifying and swapping those that happen to be 00110 // equalities (see Proposition above). 00111 for (dimension_type i = num_lines_or_equalities; i < num_rows; ++i) { 00112 if (sat[i].empty()) { 00113 // The constraint `sys[i]' is saturated by all the generators. 00114 // Thus, either it is already an equality or it can be transformed 00115 // to an equality (see Proposition above). 00116 sys[i].set_is_line_or_equality(); 00117 // Note: simple normalization already holds. 00118 sys[i].sign_normalize(); 00119 // We also move it just after all the other equalities, 00120 // so that system `sys' keeps its partial sortedness. 00121 if (i != num_lines_or_equalities) { 00122 std::swap(sys[i], sys[num_lines_or_equalities]); 00123 std::swap(sat[i], sat[num_lines_or_equalities]); 00124 std::swap(num_saturators[i], num_saturators[num_lines_or_equalities]); 00125 } 00126 ++num_lines_or_equalities; 00127 // `sys' is no longer sorted. 00128 sys.set_sorted(false); 00129 } 00130 else 00131 // There exists a generator which does not saturate `sys[i]', 00132 // so that `sys[i]' is indeed an inequality. 00133 // We store the number of its saturators. 00134 num_saturators[i] = num_cols_sat - sat[i].count_ones(); 00135 } 00136 00137 // At this point, all the equalities of `sys' (included those 00138 // inequalities that we just transformed to equalities) have 00139 // indexes between 0 and `num_lines_or_equalities' - 1, 00140 // which is the property needed by method gauss(). 00141 // We can simplify the system of equalities, obtaining the rank 00142 // of `sys' as result. 00143 const dimension_type rank = sys.gauss(num_lines_or_equalities); 00144 00145 // Now the irredundant equalities of `sys' have indexes from 0 00146 // to `rank' - 1, whereas the equalities having indexes from `rank' 00147 // to `num_lines_or_equalities' - 1 are all redundant. 00148 // (The inequalities in `sys' have been left untouched.) 00149 // The rows containing equalities are not sorted. 00150 00151 if (rank < num_lines_or_equalities) { 00152 // We identified some redundant equalities. 00153 // Moving them at the bottom of `sys': 00154 // - index `redundant' runs through the redundant equalities 00155 // - index `erasing' identifies the first row that should 00156 // be erased after this loop. 00157 // Note that we exit the loop either because we have moved all 00158 // redundant equalities or because we have moved all the 00159 // inequalities. 00160 for (dimension_type redundant = rank, 00161 erasing = num_rows; 00162 redundant < num_lines_or_equalities 00163 && erasing > num_lines_or_equalities; 00164 ) { 00165 --erasing; 00166 std::swap(sys[redundant], sys[erasing]); 00167 std::swap(sat[redundant], sat[erasing]); 00168 std::swap(num_saturators[redundant], num_saturators[erasing]); 00169 sys.set_sorted(false); 00170 ++redundant; 00171 } 00172 // Adjusting the value of `num_rows' to the number of meaningful 00173 // rows of `sys': `num_lines_or_equalities' - `rank' is the number of 00174 // redundant equalities moved to the bottom of `sys', which are 00175 // no longer meaningful. 00176 num_rows -= num_lines_or_equalities - rank; 00177 // Adjusting the value of `num_lines_or_equalities'. 00178 num_lines_or_equalities = rank; 00179 } 00180 00181 // Now we use the definition of redundancy (given in the Introduction) 00182 // to remove redundant inequalities. 00183 00184 // First we check the saturation rule, which provides a necessary 00185 // condition for an inequality to be irredundant (i.e., it provides 00186 // a sufficient condition for identifying redundant inequalities). 00187 // Let 00188 // num_saturators[i] = num_sat_lines[i] + num_sat_rays_or_points[i]; 00189 // dim_lin_space = num_irred_lines; 00190 // dim_ray_space 00191 // = dim_vector_space - num_irred_equalities - dim_lin_space 00192 // = num_columns - 1 - num_lines_or_equalities - dim_lin_space; 00193 // min_sat_rays_or_points = dim_ray_space. 00194 // 00195 // An inequality saturated by less than `dim_ray_space' _rays/points_ 00196 // is redundant. Thus we have the implication 00197 // 00198 // (num_saturators[i] - num_sat_lines[i] < dim_ray_space) 00199 // ==> 00200 // redundant(sys[i]). 00201 // 00202 // Moreover, since every line saturates all inequalities, we also have 00203 // dim_lin_space = num_sat_lines[i] 00204 // so that we can rewrite the condition above as follows: 00205 // 00206 // (num_saturators[i] < num_columns - num_lines_or_equalities - 1) 00207 // ==> 00208 // redundant(sys[i]). 00209 // 00210 const dimension_type min_saturators 00211 = num_columns - num_lines_or_equalities - 1; 00212 for (dimension_type i = num_lines_or_equalities; i < num_rows; ) { 00213 if (num_saturators[i] < min_saturators) { 00214 // The inequality `sys[i]' is redundant. 00215 --num_rows; 00216 std::swap(sys[i], sys[num_rows]); 00217 std::swap(sat[i], sat[num_rows]); 00218 std::swap(num_saturators[i], num_saturators[num_rows]); 00219 sys.set_sorted(false); 00220 } 00221 else 00222 ++i; 00223 } 00224 00225 // Now we check the independence rule. 00226 for (dimension_type i = num_lines_or_equalities; i < num_rows; ) { 00227 bool redundant = false; 00228 // NOTE: in the inner loop, index `j' runs through _all_ the 00229 // inequalities and we do not test if `sat[i]' is strictly 00230 // contained into `sat[j]'. Experimentation has shown that this 00231 // is faster than having `j' only run through the indexes greater 00232 // than `i' and also doing the test `strict_subset(sat[i], 00233 // sat[k])'. 00234 for (dimension_type j = num_lines_or_equalities; j < num_rows; ) { 00235 if (i == j) 00236 // We want to compare different rows of `sys'. 00237 ++j; 00238 else { 00239 // Let us recall that each generator lies on a facet of the 00240 // polyhedron (see the Introduction). 00241 // Given two constraints `c_1' and `c_2', if there are `m' 00242 // generators lying on the hyper-plane corresponding to `c_1', 00243 // the same `m' generators lie on the hyper-plane 00244 // corresponding to `c_2', too, and there is another one lying 00245 // on the latter but not on the former, then `c_2' is more 00246 // restrictive than `c_1', i.e., `c_1' is redundant. 00247 bool strict_subset; 00248 if (subset_or_equal(sat[j], sat[i], strict_subset)) 00249 if (strict_subset) { 00250 // All the saturators of the inequality `sys[i]' are 00251 // saturators of the inequality `sys[j]' too, 00252 // and there exists at least one saturator of `sys[j]' 00253 // which is not a saturator of `sys[i]'. 00254 // It follows that inequality `sys[i]' is redundant. 00255 redundant = true; 00256 break; 00257 } 00258 else { 00259 // We have `sat[j] == sat[i]'. Hence inequalities 00260 // `sys[i]' and `sys[j]' are saturated by the same set of 00261 // generators. Then we can remove either one of the two 00262 // inequalities: we remove `sys[j]'. 00263 --num_rows; 00264 std::swap(sys[j], sys[num_rows]); 00265 std::swap(sat[j], sat[num_rows]); 00266 std::swap(num_saturators[j], num_saturators[num_rows]); 00267 sys.set_sorted(false); 00268 } 00269 else 00270 // If we reach this point then we know that `sat[i]' does 00271 // not contain (and is different from) `sat[j]', so that 00272 // `sys[i]' is not made redundant by inequality `sys[j]'. 00273 ++j; 00274 } 00275 } 00276 if (redundant) { 00277 // The inequality `sys[i]' is redundant. 00278 --num_rows; 00279 std::swap(sys[i], sys[num_rows]); 00280 std::swap(sat[i], sat[num_rows]); 00281 std::swap(num_saturators[i], num_saturators[num_rows]); 00282 sys.set_sorted(false); 00283 } 00284 else 00285 // The inequality `sys[i]' is not redundant. 00286 ++i; 00287 } 00288 00289 // Here we physically remove the redundant inequalities previously 00290 // moved to the bottom of `sys' and the corresponding `sat' rows. 00291 sys.erase_to_end(num_rows); 00292 sys.unset_pending_rows(); 00293 sat.rows_erase_to_end(num_rows); 00294 // At this point the first `num_lines_or_equalities' rows of 'sys' 00295 // represent the irredundant equalities, while the remaining rows 00296 // (i.e., those having indexes from `num_lines_or_equalities' to 00297 // `num_rows' - 1) represent the irredundant inequalities. 00298 #ifndef NDEBUG 00299 // Check if the flag is set (that of the equalities is already set). 00300 for (dimension_type i = num_lines_or_equalities; i < num_rows; ++i) 00301 assert(sys[i].is_ray_or_point_or_inequality()); 00302 #endif 00303 00304 // Finally, since now the sub-system (of `sys') of the irredundant 00305 // equalities is in triangular form, we back substitute each 00306 // variables with the expression obtained considering the equalities 00307 // starting from the last one. 00308 sys.back_substitute(num_lines_or_equalities); 00309 00310 // The returned value is the number of irredundant equalities i.e., 00311 // the rank of the sub-system of `sys' containing only equalities. 00312 // (See the Introduction for definition of lineality space dimension.) 00313 return num_lines_or_equalities; 00314 }