Polyhedron_minimize_templates.hh

Go to the documentation of this file.

/* 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 {

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;
  }
}


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);
}

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)