00001 /* Polyhedron class implementation: simplify(). 00002 Copyright (C) 2001-2006 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 2 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 <config.h> 00024 00025 #include "Linear_Row.defs.hh" 00026 #include "Linear_System.defs.hh" 00027 #include "Saturation_Row.defs.hh" 00028 #include "Saturation_Matrix.defs.hh" 00029 #include "Polyhedron.defs.hh" 00030 00031 namespace PPL = Parma_Polyhedra_Library; 00032 00081 int 00082 PPL::Polyhedron::simplify(Linear_System& sys, Saturation_Matrix& sat) { 00083 // This method is only applied to a well-formed system `sys'. 00084 assert(sys.OK(true)); 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 static std::vector<dimension_type> num_saturators; 00099 num_saturators.reserve(num_rows); 00100 00101 // Computing the number of saturators for each inequality, 00102 // possibly identifying and swapping those that happen to be 00103 // equalities (see Proposition above). 00104 for (dimension_type i = num_lines_or_equalities; i < num_rows; ++i) { 00105 if (sat[i].empty()) { 00106 // The constraint `sys[i]' is saturated by all the generators. 00107 // Thus, either it is already an equality or it can be transformed 00108 // to an equality (see Proposition above). 00109 sys[i].set_is_line_or_equality(); 00110 // Note: simple normalization already holds. 00111 sys[i].sign_normalize(); 00112 // We also move it just after all the other equalities, 00113 // so that system `sys' keeps its partial sortedness. 00114 if (i != num_lines_or_equalities) { 00115 std::swap(sys[i], sys[num_lines_or_equalities]); 00116 std::swap(sat[i], sat[num_lines_or_equalities]); 00117 std::swap(num_saturators[i], num_saturators[num_lines_or_equalities]); 00118 } 00119 ++num_lines_or_equalities; 00120 // `sys' is no longer sorted. 00121 sys.set_sorted(false); 00122 } 00123 else 00124 // There exists a generator which does not saturate `sys[i]', 00125 // so that `sys[i]' is indeed an inequality. 00126 // We store the number of its saturators. 00127 num_saturators[i] = num_cols_sat - sat[i].count_ones(); 00128 } 00129 00130 // At this point, all the equalities of `sys' (included those 00131 // inequalities that we just transformed to equalities) have 00132 // indexes between 0 and `num_lines_or_equalities' - 1, 00133 // which is the property needed by method gauss(). 00134 // We can simplify the system of equalities, obtaining the rank 00135 // of `sys' as result. 00136 const dimension_type rank = sys.gauss(num_lines_or_equalities); 00137 00138 // Now the irredundant equalities of `sys' have indexes from 0 00139 // to `rank' - 1, whereas the equalities having indexes from `rank' 00140 // to `num_lines_or_equalities' - 1 are all redundant. 00141 // (The inequalities in `sys' have been left untouched.) 00142 // The rows containing equalities are not sorted. 00143 00144 if (rank < num_lines_or_equalities) { 00145 // We identified some redundant equalities. 00146 // Moving them at the bottom of `sys': 00147 // - index `redundant' runs through the redundant equalities 00148 // - index `erasing' identifies the first row that should 00149 // be erased after this loop. 00150 // Note that we exit the loop either because we have moved all 00151 // redundant equalities or because we have moved all the 00152 // inequalities. 00153 for (dimension_type redundant = rank, 00154 erasing = num_rows; 00155 redundant < num_lines_or_equalities 00156 && erasing > num_lines_or_equalities; 00157 ) { 00158 --erasing; 00159 std::swap(sys[redundant], sys[erasing]); 00160 std::swap(sat[redundant], sat[erasing]); 00161 std::swap(num_saturators[redundant], num_saturators[erasing]); 00162 sys.set_sorted(false); 00163 ++redundant; 00164 } 00165 // Adjusting the value of `num_rows' to the number of meaningful 00166 // rows of `sys': `num_lines_or_equalities' - `rank' is the number of 00167 // redundant equalities moved to the bottom of `sys', which are 00168 // no longer meaningful. 00169 num_rows -= num_lines_or_equalities - rank; 00170 // Adjusting the value of `num_lines_or_equalities'. 00171 num_lines_or_equalities = rank; 00172 } 00173 00174 // Now we use the definition of redundancy (given in the Introduction) 00175 // to remove redundant inequalities. 00176 00177 // First we check the saturation rule, which provides a necessary 00178 // condition for an inequality to be irredundant (i.e., it provides 00179 // a sufficient condition for identifying redundant inequalities). 00180 // Let 00181 // num_saturators[i] = num_sat_lines[i] + num_sat_rays_or_points[i]; 00182 // dim_lin_space = num_irred_lines; 00183 // dim_ray_space 00184 // = dim_vector_space - num_irred_equalities - dim_lin_space 00185 // = num_columns - 1 - num_lines_or_equalities - dim_lin_space; 00186 // min_sat_rays_or_points = dim_ray_space. 00187 // 00188 // An inequality saturated by less than `dim_ray_space' _rays/points_ 00189 // is redundant. Thus we have the implication 00190 // 00191 // (num_saturators[i] - num_sat_lines[i] < dim_ray_space) 00192 // ==> 00193 // redundant(sys[i]). 00194 // 00195 // Moreover, since every line saturates all inequalities, we also have 00196 // dim_lin_space = num_sat_lines[i] 00197 // so that we can rewrite the condition above as follows: 00198 // 00199 // (num_saturators[i] < num_columns - num_lines_or_equalities - 1) 00200 // ==> 00201 // redundant(sys[i]). 00202 // 00203 const dimension_type min_saturators 00204 = num_columns - num_lines_or_equalities - 1; 00205 for (dimension_type i = num_lines_or_equalities; i < num_rows; ) { 00206 if (num_saturators[i] < min_saturators) { 00207 // The inequality `sys[i]' is redundant. 00208 --num_rows; 00209 std::swap(sys[i], sys[num_rows]); 00210 std::swap(sat[i], sat[num_rows]); 00211 std::swap(num_saturators[i], num_saturators[num_rows]); 00212 sys.set_sorted(false); 00213 } 00214 else 00215 ++i; 00216 } 00217 00218 // Now we check the independence rule. 00219 for (dimension_type i = num_lines_or_equalities; i < num_rows; ) { 00220 bool redundant = false; 00221 // NOTE: in the inner loop, index `j' runs through _all_ the 00222 // inequalities and we do not test if `sat[i]' is strictly 00223 // contained into `sat[j]'. Experimentation has shown that this 00224 // is faster than having `j' only run through the indexes greater 00225 // than `i' and also doing the test `strict_subset(sat[i], 00226 // sat[k])'. 00227 for (dimension_type j = num_lines_or_equalities; j < num_rows; ) { 00228 if (i == j) 00229 // We want to compare different rows of `sys'. 00230 ++j; 00231 else { 00232 // Let us recall that each generator lies on a facet of the 00233 // polyhedron (see the Introduction). 00234 // Given two constraints `c_1' and `c_2', if there are `m' 00235 // generators lying on the hyper-plane corresponding to `c_1', 00236 // the same `m' generators lie on the hyper-plane 00237 // corresponding to `c_2', too, and there is another one lying 00238 // on the latter but not on the former, then `c_2' is more 00239 // restrictive than `c_1', i.e., `c_1' is redundant. 00240 bool strict_subset; 00241 if (subset_or_equal(sat[j], sat[i], strict_subset)) 00242 if (strict_subset) { 00243 // All the saturators of the inequality `sys[i]' are 00244 // saturators of the inequality `sys[j]' too, 00245 // and there exists at least one saturator of `sys[j]' 00246 // which is not a saturator of `sys[i]'. 00247 // It follows that inequality `sys[i]' is redundant. 00248 redundant = true; 00249 break; 00250 } 00251 else { 00252 // We have `sat[j] == sat[i]'. Hence inequalities 00253 // `sys[i]' and `sys[j]' are saturated by the same set of 00254 // generators. Then we can remove either one of the two 00255 // inequalities: we remove `sys[j]'. 00256 --num_rows; 00257 std::swap(sys[j], sys[num_rows]); 00258 std::swap(sat[j], sat[num_rows]); 00259 std::swap(num_saturators[j], num_saturators[num_rows]); 00260 sys.set_sorted(false); 00261 } 00262 else 00263 // If we reach this point then we know that `sat[i]' does 00264 // not contain (and is different from) `sat[j]', so that 00265 // `sys[i]' is not made redundant by inequality `sys[j]'. 00266 ++j; 00267 } 00268 } 00269 if (redundant) { 00270 // The inequality `sys[i]' is redundant. 00271 --num_rows; 00272 std::swap(sys[i], sys[num_rows]); 00273 std::swap(sat[i], sat[num_rows]); 00274 std::swap(num_saturators[i], num_saturators[num_rows]); 00275 sys.set_sorted(false); 00276 } 00277 else 00278 // The inequality `sys[i]' is not redundant. 00279 ++i; 00280 } 00281 00282 // Here we physically remove the redundant inequalities previously 00283 // moved to the bottom of `sys' and the corresponding `sat' rows. 00284 sys.erase_to_end(num_rows); 00285 sys.unset_pending_rows(); 00286 sat.rows_erase_to_end(num_rows); 00287 // At this point the first `num_lines_or_equalities' rows of 'sys' 00288 // represent the irredundant equalities, while the remaining rows 00289 // (i.e., those having indexes from `num_lines_or_equalities' to 00290 // `num_rows' - 1) represent the irredundant inequalities. 00291 #ifndef NDEBUG 00292 // Check if the flag is set (that of the equalities is already set). 00293 for (dimension_type i = num_lines_or_equalities; i < num_rows; ++i) 00294 assert(sys[i].is_ray_or_point_or_inequality()); 00295 #endif 00296 00297 // Finally, since now the sub-system (of `sys') of the irredundant 00298 // equalities is in triangular form, we back substitute each 00299 // variables with the expression obtained considering the equalities 00300 // starting from the last one. 00301 sys.back_substitute(num_lines_or_equalities); 00302 00303 // The returned value is the number of irredundant equalities i.e., 00304 // the rank of the sub-system of `sys' containing only equalities. 00305 // (See the Introduction for definition of lineality space dimension.) 00306 return num_lines_or_equalities; 00307 }