1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456
|
/* Polyhedron class implementation: minimize() and add_and_minimize().
Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)
This file is part of the Parma Polyhedra Library (PPL).
The PPL is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.
The PPL is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software Foundation,
Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
For the most up-to-date information see the Parma Polyhedra Library
site: http://bugseng.com/products/ppl/ . */
#ifndef PPL_Polyhedron_minimize_templates_hh
#define PPL_Polyhedron_minimize_templates_hh 1
#include "Bit_Matrix_defs.hh"
#include "Polyhedron_defs.hh"
#include <stdexcept>
namespace Parma_Polyhedra_Library {
/*!
\return
<CODE>true</CODE> if the polyhedron is empty, <CODE>false</CODE>
otherwise.
\param con_to_gen
<CODE>true</CODE> if \p source represents the constraints,
<CODE>false</CODE> otherwise;
\param source
The given system, which is not empty;
\param dest
The system to build and minimize;
\param sat
The saturation matrix.
\p dest is not <CODE>const</CODE> because it will be built (and then
modified) during minimize(). Also, \p sat and \p source are
not <CODE>const</CODE> because the former will be built during
\p dest creation and the latter will maybe be sorted and modified by
<CODE>conversion()</CODE> and <CODE>simplify()</CODE>.
\p sat has the generators on its columns and the constraints on its rows
if \p con_to_gen is <CODE>true</CODE>, otherwise it has the generators on
its rows and the constraints on its columns.
Given \p source, this function builds (by means of
<CODE>conversion()</CODE>) \p dest and then simplifies (invoking
<CODE>simplify()</CODE>) \p source, erasing redundant rows.
For the sequel we assume that \p source is the system of constraints
and \p dest is the system of generators.
This will simplify the description of the function; the dual case is
similar.
*/
template <typename Source_Linear_System, typename Dest_Linear_System>
bool
Polyhedron::minimize(const bool con_to_gen,
Source_Linear_System& source,
Dest_Linear_System& dest,
Bit_Matrix& sat) {
typedef typename Dest_Linear_System::row_type dest_row_type;
// Topologies have to agree.
PPL_ASSERT(source.topology() == dest.topology());
// `source' cannot be empty: even if it is an empty constraint system,
// representing the universe polyhedron, homogenization has added
// the positive constraint. It also cannot be an empty generator system,
// since this function is always called starting from a non-empty
// polyhedron.
PPL_ASSERT(!source.has_no_rows());
// Sort the source system, if necessary.
if (!source.is_sorted()) {
source.sort_rows();
}
// Initialization of the system of generators `dest'.
// The algorithm works incrementally and we haven't seen any
// constraint yet: as a consequence, `dest' should describe
// the universe polyhedron of the appropriate dimension.
// To this end, we initialize it to the identity matrix of dimension
// `source.num_columns()': the rows represent the lines corresponding
// to the canonical basis of the vector space.
dimension_type dest_num_rows
= source.topology() == NECESSARILY_CLOSED ? source.space_dimension() + 1
: source.space_dimension() + 2;
dest.clear();
dest.set_space_dimension(source.space_dimension());
// Initialize `dest' to the identity matrix.
for (dimension_type i = 0; i < dest_num_rows; ++i) {
Linear_Expression expr;
expr.set_space_dimension(dest_num_rows - 1);
if (i == 0) {
expr += 1;
}
else {
expr += Variable(i - 1);
}
dest_row_type dest_i(expr, dest_row_type::LINE_OR_EQUALITY, NECESSARILY_CLOSED);
if (dest.topology() == NOT_NECESSARILY_CLOSED) {
dest_i.mark_as_not_necessarily_closed();
}
dest.sys.insert_no_ok(dest_i, Recycle_Input());
}
// The identity matrix `dest' is not sorted (see the sorting rules
// in Constrant.cc and Generator.cc).
dest.set_sorted(false);
// NOTE: the system `dest', as it is now, is not a _legal_ system of
// generators, because in the first row we have a line with a
// non-zero divisor (which should only happen for
// points). However, this is NOT a problem, because `source'
// necessarily contains the positivity constraint (or a
// combination of it with another constraint) which will
// restore things as they should be.
// Building a saturation matrix and initializing it by setting
// all of its elements to zero. This matrix will be modified together
// with `dest' during the conversion.
// NOTE: since we haven't seen any constraint yet, the relevant
// portion of `tmp_sat' is the sub-matrix consisting of
// the first 0 columns: thus the relevant portion correctly
// characterizes the initial saturation information.
Bit_Matrix tmp_sat(dest_num_rows, source.num_rows());
// By invoking the function conversion(), we populate `dest' with
// the generators characterizing the polyhedron described by all
// the constraints in `source'.
// The `start' parameter is zero (we haven't seen any constraint yet)
// and the 5th parameter (representing the number of lines in `dest'),
// by construction, is equal to `dest_num_rows'.
const dimension_type num_lines_or_equalities
= conversion(source, 0U, dest, tmp_sat, dest_num_rows);
// conversion() may have modified the number of rows in `dest'.
dest_num_rows = dest.num_rows();
#ifndef NDEBUG
for (dimension_type i = dest.num_rows(); i-- > 0; ) {
PPL_ASSERT(dest[i].OK());
}
#endif
// Checking if the generators in `dest' represent an empty polyhedron:
// the polyhedron is empty if there are no points
// (because rays, lines and closure points need a supporting point).
// Points can be detected by looking at:
// - the divisor, for necessarily closed polyhedra;
// - the epsilon coordinate, for NNC polyhedra.
dimension_type first_point;
if (dest.is_necessarily_closed()) {
for (first_point = num_lines_or_equalities;
first_point < dest_num_rows;
++first_point) {
if (dest[first_point].expr.inhomogeneous_term() > 0) {
break;
}
}
}
else {
for (first_point = num_lines_or_equalities;
first_point < dest_num_rows;
++first_point) {
if (dest[first_point].expr.get(Variable(dest.space_dimension())) > 0) {
break;
}
}
}
if (first_point == dest_num_rows) {
if (con_to_gen) {
// No point has been found: the polyhedron is empty.
return true;
}
else {
// Here `con_to_gen' is false: `dest' is a system of constraints.
// In this case the condition `first_point == dest_num_rows'
// actually means that all the constraints in `dest' have their
// inhomogeneous term equal to 0.
// This is an ILLEGAL situation, because it implies that
// the constraint system `dest' lacks the positivity constraint
// and no linear combination of the constraints in `dest'
// can reintroduce the positivity constraint.
PPL_UNREACHABLE;
return false;
}
}
else {
// A point has been found: the polyhedron is not empty.
// Now invoking simplify() to remove all the redundant constraints
// from the system `source'.
// Since the saturation matrix `tmp_sat' returned by conversion()
// has rows indexed by generators (the rows of `dest') and columns
// indexed by constraints (the rows of `source'), we have to
// transpose it to obtain the saturation matrix needed by simplify().
sat.transpose_assign(tmp_sat);
simplify(source, sat);
return false;
}
}
/*!
\return
<CODE>true</CODE> if the obtained polyhedron is empty,
<CODE>false</CODE> otherwise.
\param con_to_gen
<CODE>true</CODE> if \p source1 and \p source2 are system of
constraints, <CODE>false</CODE> otherwise;
\param source1
The first element of the given DD pair;
\param dest
The second element of the given DD pair;
\param sat
The saturation matrix that bind \p source1 to \p dest;
\param source2
The new system of generators or constraints.
It is assumed that \p source1 and \p source2 are sorted and have
no pending rows. It is also assumed that \p dest has no pending rows.
On entry, the rows of \p sat are indexed by the rows of \p dest
and its columns are indexed by the rows of \p source1.
On exit, the rows of \p sat are indexed by the rows of \p dest
and its columns are indexed by the rows of the system obtained
by merging \p source1 and \p source2.
Let us suppose we want to add some constraints to a given system of
constraints \p source1. This method, given a minimized double description
pair (\p source1, \p dest) and a system of new constraints \p source2,
modifies \p source1 by adding to it the constraints of \p source2 that
are not in \p source1. Then, by invoking
<CODE>add_and_minimize(bool, Linear_System_Class&, Linear_System_Class&, Bit_Matrix&)</CODE>,
processes the added constraints obtaining a new DD pair.
This method treats also the dual case, i.e., adding new generators to
a previous system of generators. In this case \p source1 contains the
old generators, \p source2 the new ones and \p dest is the system
of constraints in the given minimized DD pair.
Since \p source2 contains the constraints (or the generators) that
will be added to \p source1, it is constant: it will not be modified.
*/
template <typename Source_Linear_System1, typename Source_Linear_System2,
typename Dest_Linear_System>
bool
Polyhedron::add_and_minimize(const bool con_to_gen,
Source_Linear_System1& source1,
Dest_Linear_System& dest,
Bit_Matrix& sat,
const Source_Linear_System2& source2) {
// `source1' and `source2' cannot be empty.
PPL_ASSERT(!source1.has_no_rows() && !source2.has_no_rows());
// `source1' and `source2' must have the same number of columns
// to be merged.
PPL_ASSERT(source1.num_columns() == source2.num_columns());
// `source1' and `source2' are fully sorted.
PPL_ASSERT(source1.is_sorted() && source1.num_pending_rows() == 0);
PPL_ASSERT(source2.is_sorted() && source2.num_pending_rows() == 0);
PPL_ASSERT(dest.num_pending_rows() == 0);
const dimension_type old_source1_num_rows = source1.num_rows();
// `k1' and `k2' run through the rows of `source1' and `source2', resp.
dimension_type k1 = 0;
dimension_type k2 = 0;
dimension_type source2_num_rows = source2.num_rows();
while (k1 < old_source1_num_rows && k2 < source2_num_rows) {
// Add to `source1' the constraints from `source2', as pending rows.
// We exploit the property that initially both `source1' and `source2'
// are sorted and index `k1' only scans the non-pending rows of `source1',
// so that it is not influenced by the pending rows appended to it.
// This way no duplicate (i.e., trivially redundant) constraint
// is introduced in `source1'.
const int cmp = compare(source1[k1], source2[k2]);
if (cmp == 0) {
// We found the same row: there is no need to add `source2[k2]'.
++k2;
// By sortedness, since `k1 < old_source1_num_rows',
// we can increment index `k1' too.
++k1;
}
else if (cmp < 0) {
// By sortedness, we can increment `k1'.
++k1;
}
else {
// Here `cmp > 0'.
// By sortedness, `source2[k2]' cannot be in `source1'.
// We add it as a pending row of `source1' (sortedness unaffected).
source1.add_pending_row(source2[k2]);
// We can increment `k2'.
++k2;
}
}
// Have we scanned all the rows in `source2'?
if (k2 < source2_num_rows) {
// By sortedness, all the rows in `source2' having indexes
// greater than or equal to `k2' were not in `source1'.
// We add them as pending rows of 'source1' (sortedness not affected).
for ( ; k2 < source2_num_rows; ++k2) {
source1.add_pending_row(source2[k2]);
}
}
if (source1.num_pending_rows() == 0) {
// No row was appended to `source1', because all the constraints
// in `source2' were already in `source1'.
// There is nothing left to do ...
return false;
}
return add_and_minimize(con_to_gen, source1, dest, sat);
}
/*!
\return
<CODE>true</CODE> if the obtained polyhedron is empty,
<CODE>false</CODE> otherwise.
\param con_to_gen
<CODE>true</CODE> if \p source is a system of constraints,
<CODE>false</CODE> otherwise;
\param source
The first element of the given DD pair. It also contains the pending
rows to be processed;
\param dest
The second element of the given DD pair. It cannot have pending rows;
\param sat
The saturation matrix that bind the upper part of \p source to \p dest.
On entry, the rows of \p sat are indexed by the rows of \p dest
and its columns are indexed by the non-pending rows of \p source.
On exit, the rows of \p sat are indexed by the rows of \p dest
and its columns are indexed by the rows of \p source.
Let us suppose that \p source is a system of constraints.
This method assumes that the non-pending part of \p source and
system \p dest form a double description pair in minimal form and
will build a new DD pair in minimal form by processing the pending
constraints in \p source. To this end, it will call
<CODE>conversion()</CODE>) and <CODE>simplify</CODE>.
This method treats also the dual case, i.e., processing pending
generators. In this case \p source contains generators and \p dest
is the system of constraints corresponding to the non-pending part
of \p source.
*/
template <typename Source_Linear_System, typename Dest_Linear_System>
bool
Polyhedron::add_and_minimize(const bool con_to_gen,
Source_Linear_System& source,
Dest_Linear_System& dest,
Bit_Matrix& sat) {
PPL_ASSERT(source.num_pending_rows() > 0);
PPL_ASSERT(source.space_dimension() == dest.space_dimension());
PPL_ASSERT(source.is_sorted());
// First, pad the saturation matrix with new columns (of zeroes)
// to accommodate for the pending rows of `source'.
sat.resize(dest.num_rows(), source.num_rows());
// Incrementally compute the new system of generators.
// Parameter `start' is set to the index of the first pending constraint.
const dimension_type num_lines_or_equalities
= conversion(source, source.first_pending_row(),
dest, sat,
dest.num_lines_or_equalities());
// conversion() may have modified the number of rows in `dest'.
const dimension_type dest_num_rows = dest.num_rows();
// Checking if the generators in `dest' represent an empty polyhedron:
// the polyhedron is empty if there are no points
// (because rays, lines and closure points need a supporting point).
// Points can be detected by looking at:
// - the divisor, for necessarily closed polyhedra;
// - the epsilon coordinate, for NNC polyhedra.
dimension_type first_point;
if (dest.is_necessarily_closed()) {
for (first_point = num_lines_or_equalities;
first_point < dest_num_rows;
++first_point) {
if (dest[first_point].expr.inhomogeneous_term() > 0) {
break;
}
}
}
else {
for (first_point = num_lines_or_equalities;
first_point < dest_num_rows;
++first_point) {
if (dest[first_point].expr.get(Variable(dest.space_dimension())) > 0) {
break;
}
}
}
if (first_point == dest_num_rows) {
if (con_to_gen) {
// No point has been found: the polyhedron is empty.
return true;
}
else {
// Here `con_to_gen' is false: `dest' is a system of constraints.
// In this case the condition `first_point == dest_num_rows'
// actually means that all the constraints in `dest' have their
// inhomogeneous term equal to 0.
// This is an ILLEGAL situation, because it implies that
// the constraint system `dest' lacks the positivity constraint
// and no linear combination of the constraints in `dest'
// can reintroduce the positivity constraint.
PPL_UNREACHABLE;
return false;
}
}
else {
// A point has been found: the polyhedron is not empty.
// Now invoking `simplify()' to remove all the redundant constraints
// from the system `source'.
// Since the saturation matrix `sat' returned by `conversion()'
// has rows indexed by generators (the rows of `dest') and columns
// indexed by constraints (the rows of `source'), we have to
// transpose it to obtain the saturation matrix needed by `simplify()'.
sat.transpose();
simplify(source, sat);
// Transposing back.
sat.transpose();
return false;
}
}
} // namespace Parma_Polyhedra_Library
#endif // !defined(PPL_Polyhedron_minimize_templates_hh)
|