UpdDoubleDescription_Rn.h

Go to the documentation of this file.

// politopix allows to make computations on polytopes such as finding vertices, intersecting, Minkowski sums, ...
//     Copyright (C) 2021 : Delos Vincent
//
//     This program 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.
//
//     This program 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, see <http://www.gnu.org/licenses/>.
//
// I2M (UMR CNRS 5295 / University of Bordeaux)
 
#ifndef UPDDOUBLEDESCRIPTION_Rn
#define UPDDOUBLEDESCRIPTION_Rn
 
#include <math.h>
#include <numeric>
#include <set>
#include <map>
#include "Tracking.h"
#include "Neighbours_Rn.h"
 
 
 
template< class POLYHEDRON, class ITERATOR, class REDUNDANCY_PROCESSING > class UpdDoubleDescription {
 
 public:
 UpdDoubleDescription(POLYHEDRON poly, ITERATOR ite, REDUNDANCY_PROCESSING redproc, int truncationStep):_redundancyProcessing(redproc) {
    std::vector< boost::shared_ptr<Generator_Rn_SD> > listOfGenSD;
    // Generator_Rn => Generator_Rn_SD
    // From the current generators make a list of generators dedicated to work in the data structure algorithm.
    // Make sure the list of half-space numbers is sorted as we want to use algorithms further.
    poly->getListOfGeneratorsSD(listOfGenSD);
    // Count the number of hyperplanes versus the number of half-spaces.
    unsigned int nbHyp = 0, nbProcessHyp = 0;
    {for (ite.begin(); ite.end()!=true; ite.next()) {
      const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace = ite.current();
      boost::shared_ptr<Hyperplane_Rn> currentHyp = boost::dynamic_pointer_cast<Hyperplane_Rn>(currentHalfSpace);
      if (currentHyp) {
        nbHyp++;
        if (ite.currentIteratorNumber() < truncationStep)
          nbProcessHyp++;
      }
    }}
    _redundancyProcessing.initNumberOfVerticesPerHalfSpace(listOfGenSD, truncationStep, nbProcessHyp); // Deal with redundancy.
    //_redundancyProcessing.dumpListOfRedundantHS( std::cout );
 
    // Do the hard job now.
    _isEmpty = !compute(poly, ite, truncationStep, listOfGenSD);
 
    if (_isEmpty == false) {
      // Set the pointers from the generators to the half-spaces before dealing with redundancy - as the redundancy
      // processing algorithm will change the numbers of the half-spaces. Here recycle the old pointers stored in
      // the polyhedron to build the new list with the created generators.
      poly->setListOfGeneratorsSD(listOfGenSD);
      // Include all the generators now.
      _redundancyProcessing.unhookRedundantHalfSpaces(poly);
      //poly->dump(std::cout);
      // In some fuzzy cases, computing the intersection between the space spanned by two generators and a hyperplane
      // provides a fuzzy generator somewhere inside the fuzzy zone. The location of such a generator needs to be
      // refined later on in the process. Meanwhile the current one is likely to lose half-spaces as the probability
      // that some of them are redundant is high. So remove them as they are not genuine generators.
      _redundancyProcessing.unhookRedundantGenerators(poly);
    }
  }
 
  UpdDoubleDescription(
      POLYHEDRON poly,
      ITERATOR ite,
      REDUNDANCY_PROCESSING redproc,
      int truncationStep,
      TrackingOperatorToResult& trackerVdesc,
      TrackingOperatorToResult& trackerHdesc):_redundancyProcessing(redproc) {
    std::vector< boost::shared_ptr<Generator_Rn_SD> > listOfGenSD;
    // Generator_Rn => Generator_Rn_SD
    // From the current generators make a list of generators dedicated to work in the data structure algorithm.
    // Make sure the list of half-space numbers is sorted as we want to use algorithms further.
    poly->getListOfGeneratorsSD(listOfGenSD);
    // Count the number of hyperplanes versus the number of half-spaces.
    unsigned int nbHyp = 0, nbProcessHyp = 0;
    {for (ite.begin(); ite.end()!=true; ite.next()) {
      const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace = ite.current();
      boost::shared_ptr<Hyperplane_Rn> currentHyp = boost::dynamic_pointer_cast<Hyperplane_Rn>(currentHalfSpace);
      if (currentHyp) {
        nbHyp++;
        if (ite.currentIteratorNumber() < truncationStep)
          nbProcessHyp++;
      }
    }}
    _redundancyProcessing.initNumberOfVerticesPerHalfSpace(listOfGenSD, truncationStep, nbProcessHyp); // Deal with redundancy.
    //_redundancyProcessing.dumpListOfRedundantHS( std::cout );
 
    // Do the hard job now.
    _isEmpty = !compute(poly, ite, truncationStep, listOfGenSD);
 
    if (_isEmpty == false) {
      // Store the modifications in the binary operation tracker.
      // Operator1_Entity[i] = (Generator_Rn_SD::Status, Result_Entity[j])
      // Operator2_Entity[k] = (Generator_Rn_SD::Status, Result_Entity[l])
      //    Result_Entity[m] = (Generator_Rn_SD::Status, (Operator1_Entity[p], Operator2_Entity[q]))
      unsigned int nbRes=0;
      std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGenSD;
      {for (iteGenSD=listOfGenSD.begin(); iteGenSD!=listOfGenSD.end(); ++iteGenSD) {
        // In the tracker, the operand entities have been set to Operator_DELETED by default and the
        // result entities have been set by default to Result_UNCHANGED.
        Generator_Rn_SD::Status thisStatus = (*iteGenSD)->getStatus();
        unsigned int nbOp1 = (*iteGenSD)->getGeneratorNumber();
        if (thisStatus == Generator_Rn_SD::CREATED || thisStatus == Generator_Rn_SD::CREATED_AND_MODIFIED) {
          trackerVdesc.setResultEntityAsCreated(nbRes);
          trackerVdesc.setResult_Operator(nbRes,-1);
        }
        else if (thisStatus == Generator_Rn_SD::MODIFIED) {
          trackerVdesc.setResultEntityAsModified(nbRes);
          trackerVdesc.setResult_Operator(nbRes,nbOp1);
          trackerVdesc.setOperatorEntityAsModified(nbOp1);
          trackerVdesc.setOperator_Result(nbOp1,nbRes);
        }
        else if (thisStatus == Generator_Rn_SD::UNCHANGED) {
          trackerVdesc.setResultEntityAsUnchanged(nbRes);
          trackerVdesc.setResult_Operator(nbRes,nbOp1);
          trackerVdesc.setOperatorEntityAsUnchanged(nbOp1);
          trackerVdesc.setOperator_Result(nbOp1,nbRes);
        }
        ++nbRes;
      }}
      // Set the pointers from the generators to the half-spaces before dealing with redundancy - as the
      // redundancy processing algorithm will change the numbers of the half-spaces.
      poly->setListOfGeneratorsSD(listOfGenSD);
      // Process the half-spaces now.
      _redundancyProcessing.markHdescription(trackerHdesc, truncationStep);
      _redundancyProcessing.unhookRedundantHalfSpaces(poly);
      //poly->dump(std::cout);
      // In some fuzzy cases, computing the intersection between the space spanned by two generators and a hyperplane
      // provides a fuzzy generator somewhere inside the fuzzy zone. The location of such a generator needs to be
      // refined later on in the process. Meanwhile the current one is likely to lose half-spaces as the probability
      // that some of them are redundant is high. So remove them as they are not genuine generators.
      _redundancyProcessing.unhookRedundantGenerators(poly);
    }
  }
 
  virtual ~UpdDoubleDescription() {}
 
  bool getIsEmpty() const {return _isEmpty;}
 
  void computeVertexStates(
    std::vector< boost::shared_ptr<Generator_Rn_SD> >& GN_list,
    const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace,
    std::vector<double>& GN_IN_sp,
    std::vector<double>& GN_OUT_sp,
    std::vector< boost::shared_ptr<Generator_Rn_SD> >& GN_IN,
    std::vector< boost::shared_ptr<Generator_Rn_SD> >& GN_OUT,
    std::vector< boost::shared_ptr<Generator_Rn_SD> >& GN_ON)
  {
    int GeneratorNumber=0;
    double TOL = Rn::getTolerance();
    double halfSpaceNorm =
      std::inner_product(currentHalfSpace->begin(), currentHalfSpace->end(), currentHalfSpace->begin(), 0.);
    halfSpaceNorm = sqrt(halfSpaceNorm);
    std::vector< boost::shared_ptr<Generator_Rn_SD> >::iterator iteGN;
    for (iteGN=GN_list.begin(); iteGN!=GN_list.end(); ++iteGN) {
      //const boost::shared_ptr<Generator_Rn_SD>& currentGenerator = iteGN.current();
      double scalarProduct =
          std::inner_product((*iteGN)->begin(), (*iteGN)->end(), currentHalfSpace->begin(), 0.);
          //boost::numeric::ublas::inner_prod(currentGenerator->vect(), currentHalfSpace->vect());
#ifdef DEBUG
      unsigned int RnDIM=currentHalfSpace->dimension();
      std::cout.precision(15);
      std::cout << "# V" << GeneratorNumber << " = [";
      {for (unsigned int ii=0; ii<RnDIM; ii++) {
        std::cout << (*iteGN)->getCoordinate(ii);
        if (ii != RnDIM-1)
          std::cout << " ";
      }}
      std::cout << "]" << std::endl << "{ ";
      {for (unsigned int ii=0; ii<(*iteGN)->numberOfFacets(); ii++) {
        std::cout << (*iteGN)->getFacet(ii) << " ";
      }}
      std::cout << "}" << std::endl;
      std::cout << "dotP = " << scalarProduct;
      std::cout << ", dist = " << (scalarProduct+currentHalfSpace->getConstant()) / halfSpaceNorm;
#endif
      //std::cout << "sp = " << scalarProduct << " ";
      double distanceToHyperplane = (scalarProduct+currentHalfSpace->getConstant()) / halfSpaceNorm;
      if (distanceToHyperplane > TOL) {
        GN_IN.push_back((*iteGN));
        GN_IN_sp.push_back(scalarProduct);
      }
      else if (distanceToHyperplane < -TOL) {
        GN_OUT.push_back((*iteGN));
        GN_OUT_sp.push_back(scalarProduct);
      }
      else {
        GN_ON.push_back((*iteGN));
      }
#ifdef DEBUG
      HalfSpace_Rn::State currentState = HalfSpace_Rn::hs_UNKNOWN;
      if (distanceToHyperplane > TOL)
        currentState = HalfSpace_Rn::hs_IN;
      else if (distanceToHyperplane < -TOL)
        currentState = HalfSpace_Rn::hs_OUT;
      else
        currentState = HalfSpace_Rn::hs_ON;
      std::cout << ", state = " << HalfSpace_Rn::getStateAsText(currentState) << std::endl;
#endif
      GeneratorNumber++;
    }
  }
 
  bool compute(POLYHEDRON poly, ITERATOR iteHS, int truncationStep, std::vector< boost::shared_ptr<Generator_Rn_SD> >& listOfGenSD) {
    // Store the scalar product result.
    std::vector<double> GN_IN_sp;
    std::vector<double> GN_OUT_sp;
    // Store generators IN, ON, OUT or newly created.
    std::vector< boost::shared_ptr<Generator_Rn_SD> > GN_IN;
    std::vector< boost::shared_ptr<Generator_Rn_SD> > GN_OUT;
    std::vector< boost::shared_ptr<Generator_Rn_SD> > GN_ON;
    std::vector< boost::shared_ptr<Generator_Rn_SD> > GN_new;
    // Remove redundant half spaces.
    std::vector< boost::shared_ptr<HalfSpace_Rn> > redundantHS;
    int halfspaceNumber=truncationStep, generatorNumber=0;
    unsigned int RnDIM=poly->dimension();
    double TOL2=Rn::getTolerance()*Rn::getTolerance();
 
#ifdef DEBUG
    std::cout << std::endl << "UpdDoubleDescription::compute(" << truncationStep << ")" << std::endl;
#endif
    // We use this algorithm to compute vertices for a single polyhedron or to intersect two polyhedra.
    // In the last case we truncate the polyhedron A generators net by the the polyhedron B half-spaces
    // so we need to step forward in the half-space list where B constraints begin.
    //constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS(poly->getListOfHalfSpaces());
    generatorNumber = poly->numberOfGenerators();
    iteHS.setStep(truncationStep);
    {for (iteHS.begin(); iteHS.end()!=true; iteHS.next()) {
      const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace = iteHS.current();
      boost::shared_ptr<Hyperplane_Rn> currentHyp = boost::dynamic_pointer_cast<Hyperplane_Rn>(currentHalfSpace);
      bool isHyp = false;
      if (currentHyp)
        isHyp = true;
      // When you add a linear constraint tell whether it is a half-space or a hyperplane.
      _redundancyProcessing.addHalfSpace(isHyp);
#ifdef DEBUG
      double halfSpaceNorm = std::inner_product(currentHalfSpace->begin(), currentHalfSpace->end(), currentHalfSpace->begin(), 0.);
      halfSpaceNorm = sqrt(halfSpaceNorm);
      std::cout << "Facet number " << halfspaceNumber << std::endl;
      std::cout << "H" << halfspaceNumber << " = ";
      currentHalfSpace->dump(std::cout);
      std::cout << ", norm=" << halfSpaceNorm;
      std::cout << std::endl;
#endif
 
      // Fill the arrays of generators according to their position with respect to the current half-space.
      computeVertexStates(listOfGenSD, currentHalfSpace, GN_IN_sp, GN_OUT_sp, GN_IN, GN_OUT, GN_ON);
 
      // Check whether the current constraint, i.e. a half-space or a hyperplane, excludes all Generators.
      // In the case of a half-space, it used to be: if (GN_IN.empty() == true) { 
      if (currentHalfSpace->testEmptyness(GN_IN, GN_OUT, GN_ON)) {
         // NO INTERSECTION
#ifdef DEBUG
        std::cout << "Truncation is empty or flat." << std::endl;
#endif
        poly->reset();
        return false;
      }
 
      // Check whether the current constraint is redundant. Only available for polytopes.
      // For polyhedral cones, replace Rn::getDimension() by (Rn::getDimension()-1).
      if (GN_OUT.empty() == true) {// && VX_ON.size() < numberOfGeneratorsPerFacet()) {
        // REDUNDANT HALFSPACE (be careful when modifying a list under iterator)
        redundantHS.push_back(currentHalfSpace);
        //_listOfHalfSpaces.removeCurrentHalfSpace(iteHS.currentIteratorNumber());
        //_listOfHalfSpaces.erase(HSPos);
        // When chopping a polyhedral body with a hyperplane, if the set of generators OUT is empty, 
        // only keep the generators ON and ignore the generators IN.
        // In the case of a half-space, keep things unchanged i.e. listOfGeneratorSD = GN_IN + GN_ON
        currentHalfSpace->noGeneratorsOUT(listOfGenSD, GN_ON);
#ifdef DEBUG
        std::cout << "VX_OUT= " << GN_OUT.size() << std::endl;
        std::cout << "VX_IN = " <<  GN_IN.size() << std::endl;
        std::cout << "VX_ON = " <<  GN_ON.size() << std::endl;
        std::cout << "redundantHS = " <<  redundantHS.size() << std::endl;
        std::cout << "============================================" << std::endl;
#endif
      }
      else if (GN_IN.empty() != true) {
        // Build the list of all possible new generators.
        Neighbours_Rn newGeneratorsTool;
        {
          // OUT-IN & OUT-ON
          unsigned int count_OUT=0;
          std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_OUT;
          {for (iteGN_OUT=GN_OUT.begin(); iteGN_OUT!=GN_OUT.end(); ++iteGN_OUT, ++count_OUT) {
            unsigned int count_IN=0;
            std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_IN;
            {for (iteGN_IN=GN_IN.begin(); iteGN_IN!=GN_IN.end(); ++iteGN_IN, ++count_IN) {
              std::vector< unsigned int > commonPFacets;
              if (_redundancyProcessing.checkNeighbours(poly, *iteGN_IN, *iteGN_OUT, commonPFacets) == true) {
                // Uncomment only when computing 2D or 3D non fuzzy cases.
                //if (RnDIM > 3)
                  newGeneratorsTool.addGenerator(
                                commonPFacets,
                                count_IN,
                                count_OUT,
                                HalfSpace_Rn::hs_IN_OR_OUT);
                //else
                  //newGeneratorsTool.addGeneratorWithoutCheck(
                                //commonPFacets,
                                //count_IN,
                                //count_OUT,
                                //HalfSpace_Rn::hs_IN_OR_OUT);
              }
            }}
            unsigned int count_ON=0;
            std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_ON;
            {for (iteGN_ON=GN_ON.begin(); iteGN_ON!=GN_ON.end(); ++iteGN_ON, ++count_ON) {
              std::vector< unsigned int > commonPFacets;
              if (_redundancyProcessing.checkNeighbours(poly, *iteGN_ON, *iteGN_OUT, commonPFacets) == true) {
                //if (RnDIM > 3)
                  newGeneratorsTool.addGenerator(
                                commonPFacets,
                                count_ON,
                                count_OUT,
                                HalfSpace_Rn::hs_ON);
                //else
                  //newGeneratorsTool.addGeneratorWithoutCheck(
                                //commonPFacets,
                                //count_ON,
                                //count_OUT,
                                //HalfSpace_Rn::hs_ON);
              }
            }}
          }}
          // ON-IN & ON-ON
          unsigned int count_ON1=0;
          std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_ON1;
          {for (iteGN_ON1=GN_ON.begin(); iteGN_ON1!=GN_ON.end(); ++iteGN_ON1, ++count_ON1) {
            unsigned int count_IN=0;
            std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_IN;
            {for (iteGN_IN=GN_IN.begin(); iteGN_IN!=GN_IN.end(); ++iteGN_IN, ++count_IN) {
              std::vector< unsigned int > commonPFacets;
              if (_redundancyProcessing.checkNeighbours(poly, *iteGN_IN, *iteGN_ON1, commonPFacets) == true) {
                //if (RnDIM > 3)
                  newGeneratorsTool.addGenerator(
                                commonPFacets,
                                count_IN,
                                count_ON1,
                                HalfSpace_Rn::hs_ON);
                //else
                  //newGeneratorsTool.addGeneratorWithoutCheck(
                                //commonPFacets,
                                //count_IN,
                                //count_ON1,
                                //HalfSpace_Rn::hs_ON);
              }
            }}
            unsigned int count_ON2=0;
            std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_ON2;
            {for (iteGN_ON2=GN_ON.begin(); iteGN_ON2!=GN_ON.end(); ++iteGN_ON2, ++count_ON2) {
              std::vector< unsigned int > commonPFacets;
              if (iteGN_ON1 != iteGN_ON2 &&
                  _redundancyProcessing.checkNeighbours(poly, *iteGN_ON1, *iteGN_ON2, commonPFacets) == true) {
                //if (RnDIM > 3)
                  newGeneratorsTool.addGenerator(
                                commonPFacets,
                                count_ON1,
                                count_ON2,
                                HalfSpace_Rn::hs_ON);
                //else
                  //newGeneratorsTool.addGeneratorWithoutCheck(
                                //commonPFacets,
                                //count_ON1,
                                //count_ON2,
                                //HalfSpace_Rn::hs_ON);
              }
            }}
          }}
        }
        //newGeneratorsTool.dump(std::cout);
        for (newGeneratorsTool.begin(); newGeneratorsTool.end()!=true; newGeneratorsTool.next()) {
#ifdef DEBUG
          std::cout << "NGB1= " << newGeneratorsTool.currentGenInNumber() << ", NGB2= " << newGeneratorsTool.currentGenOutNumber() << std::endl;
#endif
          // Create a new Generator as a combination of currentGeneratorIn and currentGeneratorOut.
          double ay = GN_OUT_sp[newGeneratorsTool.currentGenOutNumber()];
          double az =  GN_IN_sp[newGeneratorsTool.currentGenInNumber()];
          // currentGeneratorIn and currentGeneratorOut are genuine neighbours.
          const boost::shared_ptr<Generator_Rn_SD>& currentGeneratorIn  =  GN_IN[newGeneratorsTool.currentGenInNumber()];
          const boost::shared_ptr<Generator_Rn_SD>& currentGeneratorOut = GN_OUT[newGeneratorsTool.currentGenOutNumber()];
          // Mark the new as CREATED.
          boost::shared_ptr<Generator_Rn_SD> newV(new Generator_Rn_SD(RnDIM, generatorNumber, Generator_Rn_SD::CREATED));
          ++generatorNumber;
          poly->createTruncatedGenerator(currentGeneratorOut, currentGeneratorIn, newV, ay, az, currentHalfSpace->getConstant());
          // Fill in the list of facets for the new Generator.
          std::vector< unsigned int > commonFacets;
          // Get facets.
          _redundancyProcessing.checkNeighbours(poly, currentGeneratorIn, currentGeneratorOut, commonFacets);
          newV->setAllFacets(commonFacets);
          // Do not add the current half-space at this step but wait a little bit more.
 
          // The new Generator is ready to be inserted in the data structure if not equal to
          // another one with state ON.
          bool notEq=true;
          std::set< unsigned int > setOfFacets;
          std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_ON_eq;
          {for (iteGN_ON_eq=GN_ON.begin(); iteGN_ON_eq!=GN_ON.end() && notEq==true; ++iteGN_ON_eq) {
            if ((*iteGN_ON_eq)->isEqual2(newV, RnDIM, TOL2)) {
              notEq = false;
              // newV has no future so just transfer its facets to the current equal generator.
              // Make sure not to include twice currentHalfSpace.
              newV->exportFacets(setOfFacets);
              (*iteGN_ON_eq)->exportFacets(setOfFacets);
              (*iteGN_ON_eq)->importFacets(setOfFacets);
              (*iteGN_ON_eq)->orderFacets();
            }
          }}
          std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_NEW_eq;
          {for (iteGN_NEW_eq=GN_new.begin(); iteGN_NEW_eq!=GN_new.end() && notEq==true; ++iteGN_NEW_eq) {
            if ((*iteGN_NEW_eq)->isEqual2(newV, RnDIM, TOL2)) {
              notEq = false;
              // newV has no future so just transfer its facets to the current equal generator.
              newV->setFacet(halfspaceNumber);
              newV->exportFacets(setOfFacets);
              (*iteGN_NEW_eq)->exportFacets(setOfFacets);
              (*iteGN_NEW_eq)->importFacets(setOfFacets);
              (*iteGN_NEW_eq)->orderFacets();
            }
          }}
          if (notEq == true) {
            newV->setFacet(halfspaceNumber);
            GN_new.push_back(newV);
#ifdef DEBUG
            std::cout << "New Generator = [";
            newV->save(std::cout);
            std::cout << "] GeneratorIN [";
            currentGeneratorIn->save(std::cout);
            std::cout << "], VX_IN_sp=" << az << "  GeneratorOUT [";
            currentGeneratorOut->save(std::cout);
            std::cout << "], VX_OUT_sp=" << ay << std::endl;
#endif
          }
          else {
#ifdef DEBUG
            std::cout << "No new Generator as [";
            newV->save(std::cout);
            std::cout << "] is equal to a previous one." << std::endl;
#endif
          }
        } // for (newGeneratorsTool.begin(); newGeneratorsTool.end()!=true; newGeneratorsTool.next()) {
 
#ifdef DEBUG
        std::cout << "VX_OUT= " << GN_OUT.size() << std::endl;
        std::cout << "VX_IN = " <<  GN_IN.size() << std::endl;
        std::cout << "VX_ON = " <<  GN_ON.size() << std::endl;
        std::cout << "redundantHS = " <<  redundantHS.size() << std::endl;
        std::cout << "============================================" << std::endl;
#endif
        // When we perform a flat truncation the set of generators IN
        // . is transferred into the set of generators OUT as both of them will be dropped
        // . then it is emptied before receiving the generators ON and NEW i.e. the ones we keep
        // In the other case of an intersection with a half-space this function does nothing.
        currentHalfSpace->prepareSetOfGenerators(GN_IN, GN_OUT);
 
        // All the generators with state ON must have a new facet.
        std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_ON;
        for (iteGN_ON=GN_ON.begin(); iteGN_ON!=GN_ON.end(); ++iteGN_ON) {
          (*iteGN_ON)->setFacet(halfspaceNumber);
          if ((*iteGN_ON)->getStatus() == Generator_Rn_SD::CREATED)
            (*iteGN_ON)->setStatus(Generator_Rn_SD::CREATED_AND_MODIFIED);
          else if ((*iteGN_ON)->getStatus() == Generator_Rn_SD::UNCHANGED)
            (*iteGN_ON)->setStatus(Generator_Rn_SD::MODIFIED);
          // Insert the ON generators in the IN list as this list will be the official one soon.
          GN_IN.push_back(*iteGN_ON);
        }
        // The current face must mark all the vertices with state ON. 
        // No need to do more than marking as all vertices ON are now embedded in GN_IN.
        _redundancyProcessing.updateNumberOfVerticesPerHalfSpace(halfspaceNumber, GN_ON); // Deal with redundancy.
 
        std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_new;
        for (iteGN_new=GN_new.begin(); iteGN_new!=GN_new.end(); ++iteGN_new) {
          // Insert the new generators in the IN list as this list will be the official one soon.
          GN_IN.push_back(*iteGN_new);
          _redundancyProcessing.incrementNumberForVerticesForHalfSpace(*iteGN_new); // Deal with redundancy.
        }
        // For the redundancy processing :
        // . all the vertices with state NEW must have their face counter incremented++
        // . all the vertices with state ON  must have their face counter incremented++ (in a global way)
        // . all the vertices with state OUT must have their face counter decremented--
 
        std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_OUT;
        for (iteGN_OUT=GN_OUT.begin(); iteGN_OUT!=GN_OUT.end(); ++iteGN_OUT) {
          _redundancyProcessing.decrementNumberForVerticesForHalfSpace(*iteGN_OUT); // Deal with redundancy.
        }
 
        _redundancyProcessing.updateListOfRedundantHalfSpaces(GN_IN, poly->numberOfGeneratorsPerFacet());  // Deal with redundancy.
        
        // Remove all the generators OUT and replace them with the set IN-ON-new.
        listOfGenSD = GN_IN;
 
      } // else if (GN_IN.empty() != true) {
 
      GN_OUT_sp.clear();
      GN_OUT.clear();
      GN_IN_sp.clear();
      GN_IN.clear();
      GN_ON.clear();
      GN_new.clear();
      halfspaceNumber++;
      //_redundancyProcessing.dumpSD(std::cout);
    }} // for (iteHS.begin(); iteHS.end()!=true; iteHS.next()) {
 
    return true;
  }
 
 
 protected:
  bool _isEmpty;
  REDUNDANCY_PROCESSING _redundancyProcessing;
 
};
 
 
 
template< class POLYHEDRON > class NoRedundancyProcessing {
 public:
  NoRedundancyProcessing() {}
 
  virtual ~NoRedundancyProcessing() {}
 
  virtual bool checkNeighbours(
    POLYHEDRON poly,
    const boost::shared_ptr<Generator_Rn>& genIn,
    const boost::shared_ptr<Generator_Rn>& genOut,
    std::vector< boost::shared_ptr<HalfSpace_Rn> >& commonFacets) {
    return poly->checkNeighbours(genIn, genOut, commonFacets);
  }
 
  virtual bool checkNeighbours(
    POLYHEDRON poly,
    const boost::shared_ptr<Generator_Rn_SD>& genIn,
    const boost::shared_ptr<Generator_Rn_SD>& genOut,
    std::vector< unsigned int >& commonFacets) {
    return poly->checkNeighbours(genIn, genOut, commonFacets);
  }
 
  void addHalfSpace(bool isHyperplane) {}
 
  void updateNumberOfVerticesPerHalfSpace(
    unsigned int HS,
    const std::vector< boost::shared_ptr<Generator_Rn_SD> >& GN_ON) {}
 
  virtual bool checkNeighbours(
    POLYHEDRON poly,
    const boost::shared_ptr<Generator_Rn>& genIn,
    const boost::shared_ptr<Generator_Rn>& genOut,
    std::vector< HalfSpace_Rn* >& commonFacets) {
    return poly->checkNeighbours(genIn, genOut, commonFacets);
  }
 
  void initNumberOfVerticesPerHalfSpace(const std::vector< boost::shared_ptr<Generator_Rn_SD> >&, unsigned int nbHS, unsigned int nbHP=0) {}
 
  void updateNumberOfVerticesPerHalfSpace(const boost::shared_ptr<HalfSpace_Rn>& , unsigned int) {}
 
  void updateNumberOfVerticesPerHalfSpace(const boost::shared_ptr<HalfSpace_Rn>& , const std::vector< boost::shared_ptr<Generator_Rn_SD> >&) {}
 
  void incrementNumberForVerticesForHalfSpace(const boost::shared_ptr<Generator_Rn_SD>&) {}
 
  void decrementNumberForVerticesForHalfSpace(const boost::shared_ptr<Generator_Rn_SD>&) {}
 
  void updateListOfRedundantHalfSpaces(const std::vector< boost::shared_ptr<Generator_Rn_SD> >& listOfGen, int numberOfGeneratorsPerFacet) {}
 
  void unhookRedundantHalfSpaces(POLYHEDRON) {}
 
  void unhookRedundantGenerators(POLYHEDRON) {}
 
  void dumpSD(std::ostream &this_stream) {}
 
};
 
 
template< class POLYHEDRON > class StrongRedundancyProcessing {
 
 public:
 StrongRedundancyProcessing(unsigned int ngbCond):_numberOfProcessedLinearConstraints(0),_numberOfProcessedHyperplanes(0) {
    _neigbourhoodCondition = ngbCond;
  }
 
  virtual ~StrongRedundancyProcessing() {}
 
  void fillListOfRedundantHS(
    const POLYHEDRON,
    std::vector< unsigned int >&,
    std::set< unsigned int >& getListOfRedundantHS) {
    getListOfRedundantHS = _listOfRedundantHS;
  }
 
  void initNumberOfVerticesPerHalfSpace(const std::vector< boost::shared_ptr<Generator_Rn_SD> >& LG, unsigned int nbHS, unsigned int nbHP=0) {
    _allGenPerLinearConstraint.clear();
    _isHalfSpace.clear();
    _listOfRedundantHS.clear();
    _numberOfProcessedLinearConstraints = nbHS;
    _numberOfProcessedHyperplanes = nbHP;
 
    // First init to 0 to be able to increment
    {for (unsigned int j=0; j<nbHS; ++j) {
      _isHalfSpace.push_back(true);
    }}
 
    //std::map< boost::shared_ptr<HalfSpace_Rn>, unsigned int >::const_iterator
        //mit(numberOfVertexPerHalfSpace.begin()), mend(numberOfVertexPerHalfSpace.end());
    //for (; mit!=mend; ++mit)
      //std::cout << mit->first << '\t' << mit->second << std::endl;
 
    // Now deal with the set of vertices itself contained by the faces
    // to see if a given set can be included into another one. But first of
    // all set the list of back pointers.
    {for (unsigned int j=0; j<nbHS; ++j) {
      std::set< unsigned int > genSet;
      _allGenPerLinearConstraint.push_back(genSet);
    }}
    std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator constITER_GN;
    for (constITER_GN=LG.begin(); constITER_GN!=LG.end(); ++constITER_GN) {
      {for (unsigned int j=0; j<(*constITER_GN)->numberOfFacets(); j++) {
        unsigned int F = (*constITER_GN)->getFacet(j);
        // GenPerFacet: number(HalfSpace_Rn) => {Generator_Rn_SD_0*, Generator_Rn_SD_1*, Generator_Rn_SD_2*, ...}
        // Fill the data structure storing all the generators for a given facet.
        _allGenPerLinearConstraint[F].insert( (*constITER_GN)->getGeneratorNumber() );
      }}
    }
    //dumpSD(std::cout);
  }
 
  void addHalfSpace(bool isHyperplane) {
    std::set< unsigned int > genSet;
    _allGenPerLinearConstraint.push_back(genSet);
    _numberOfProcessedLinearConstraints++;
    if (isHyperplane) {
      _numberOfProcessedHyperplanes++;
      _isHalfSpace.push_back(false);
    }
    else
      _isHalfSpace.push_back(true);
    //std::cout << "size _allGenPerLinearConstraint = " << _allGenPerLinearConstraint.size() << "size _isHalfSpace = " << _isHalfSpace.size() << "size _numberOfVerticesPerLinearConstraint = " << _numberOfVerticesPerLinearConstraint.size() << std::endl;
  }
 
  void updateNumberOfVerticesPerHalfSpace(
    unsigned int HS,
    const std::vector< boost::shared_ptr<Generator_Rn_SD> >& GN_ON) {
    std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_ON;
    for (iteGN_ON=GN_ON.begin(); iteGN_ON!=GN_ON.end(); ++iteGN_ON)
      _allGenPerLinearConstraint[ HS ].insert( (*iteGN_ON)->getGeneratorNumber() );
  }
 
  void incrementNumberForVerticesForHalfSpace(const boost::shared_ptr<Generator_Rn_SD>& GEN) {
    for (unsigned int l=0; l<GEN->numberOfFacets(); ++l) {
      _allGenPerLinearConstraint[ GEN->getFacet(l) ].insert( GEN->getGeneratorNumber() );
    }
  }
 
  void decrementNumberForVerticesForHalfSpace(const boost::shared_ptr<Generator_Rn_SD>& GEN) {
    for (unsigned int l=0; l<GEN->numberOfFacets(); ++l) {
      _allGenPerLinearConstraint[ GEN->getFacet(l) ].erase( GEN->getGeneratorNumber() );
    }
  }
 
  bool checkNeighbours(
    POLYHEDRON poly,
    const boost::shared_ptr<Generator_Rn_SD>& genIn,
    const boost::shared_ptr<Generator_Rn_SD>& genOut,
    std::vector< unsigned int >& commonFacets) {
    // Compute the intersection on sorted lists of facets.
    std::set_intersection(genIn->facetsBegin(), genIn->facetsEnd(), genOut->facetsBegin(), genOut->facetsEnd(),
      std::inserter(commonFacets, commonFacets.end()));
    // For polyhedral cones neigbourhoodCondition() is Rn::getDimension()-1, for polytopes it is Rn::getDimension()-2.
    if (commonFacets.size() >= _neigbourhoodCondition)
      return true;
 
    return false;
  }
 
  void unhookThisRedundantHalfSpace(std::vector< boost::shared_ptr<Generator_Rn_SD> >& listOfGen, 
                                    unsigned int thisRedundantHS, unsigned int numberOfGeneratorsPerFacetWithoutHyp) {
    unsigned int numberOfGeneratorsPerFacet = numberOfGeneratorsPerFacetWithoutHyp - _numberOfProcessedHyperplanes;
    std::set< unsigned int >& allGenForThisRedundantHS = _allGenPerLinearConstraint[thisRedundantHS];
   {
      std::vector< boost::shared_ptr<Generator_Rn_SD> >::iterator itGN = listOfGen.begin();
      while (itGN != listOfGen.end() || allGenForThisRedundantHS.empty() == false) {
        if (allGenForThisRedundantHS.find((*itGN)->getGeneratorNumber()) != allGenForThisRedundantHS.end()) {
          // This generator has to be modified as it is part of the list of generators that point to the redundant facet.
          // The list of facets is ordered.
          unsigned int indexOfFacet = 0;
          while ((*itGN)->getFacet(indexOfFacet) < thisRedundantHS)
            indexOfFacet++;
          if ((*itGN)->getFacet(indexOfFacet) == thisRedundantHS) {
            (*itGN)->removeFacet(indexOfFacet);
          }
          else {
            std::cerr << "indexOfGen = " << itGN-listOfGen.begin() << std::endl;
            std::cerr << "indexOfFacet = " << indexOfFacet<< std::endl;
            std::cerr << (*itGN)->vect() << std::endl;
            std::cerr << "thisRedundantHS = " << thisRedundantHS << std::endl;
            std::string errorMessage("StrongRedundancyProcessing::unhookThisRedundantHalfSpace() Facet not found.");
            throw std::out_of_range(errorMessage);
          }
          allGenForThisRedundantHS.erase((*itGN)->getGeneratorNumber() );
          // May be the generator that lost a linear constraint does not hold anymore the right number of contstraints?
          if ((*itGN)->numberOfFacets() < numberOfGeneratorsPerFacet) {
            unsigned int genNb = (*itGN)->getGeneratorNumber();
            {
              for (unsigned int idFct=0; idFct<(*itGN)->numberOfFacets(); ++idFct) {
                // Run the list of linear constraints to remove genNb each time it is present.
                _allGenPerLinearConstraint[ (*itGN)->getFacet(idFct) ].erase(genNb);
              }
            }
            // An iterator pointing to the new location of the element that followed the last element erased by the function call. 
            // This is the container end if the operation erased the last element in the sequence.
            itGN = listOfGen.erase(itGN);
          }
          else 
            itGN++;
        }
        else
          itGN++;
      }
    }
 
    if (_allGenPerLinearConstraint[thisRedundantHS].empty() == false) {
      std::cerr << std::endl << "Remaining generators: ";
      std::copy(_allGenPerLinearConstraint[thisRedundantHS].begin(), _allGenPerLinearConstraint[thisRedundantHS].end(), std::ostream_iterator<unsigned int>(std::cerr, " ") );
      std::cerr << std::endl;
      std::string errorMessage("StrongRedundancyProcessing::unhookThisRedundantHalfSpace() Remaining generators.");
      throw std::domain_error(errorMessage);
    }
 
    //std::cout << "Leaving unhookThisRedundantHalfSpace()" << std::endl;
  }
  // number(HalfSpace_Rn1) => {Generator_Rn0*, Generator_Rn1*, Generator_Rn2*, ...}
  // number(HalfSpace_Rn2) => {Generator_Rn3*, Generator_Rn4*, Generator_Rn5*, ...}
  // ...
 
  virtual void updateListOfRedundantHalfSpaces(std::vector< boost::shared_ptr<Generator_Rn_SD> >& listOfGen, int numberOfGeneratorsPerFacetWithoutHyp) {
    // Take into account the case where we have hyperplanes.
    // NB: in the future do not only take into account the number of hyperplanes but also the
    // dimension of the space they create intersecting each other. Here we make the assumption
    // we have a full rank matrix.
    unsigned int numberOfGeneratorsPerFacet = numberOfGeneratorsPerFacetWithoutHyp - _numberOfProcessedHyperplanes;
    //std::cout << "poly->numberOfGeneratorsPerFacet() " << poly->numberOfGeneratorsPerFacet() << std::endl;
    //std::cout << "numberOfGeneratorsPerFacet = " << numberOfGeneratorsPerFacet << ", _numberOfProcessedHyperplanes = " << _numberOfProcessedHyperplanes << std::endl;
    //std::cout << "_numberOfVerticesPerLinearConstraint.size() = " << _numberOfVerticesPerLinearConstraint.size() << std::endl;
    {for (unsigned int i=0; i<_allGenPerLinearConstraint.size(); ++i) {
      // Stop when we reach the current half-space being processed, it was the last inserted.
      //std::cout << "_numberOfVerticesPerLinearConstraint[i] = " << _numberOfVerticesPerLinearConstraint[i] << std::endl;
      if (_allGenPerLinearConstraint[i].size()  < numberOfGeneratorsPerFacet) {
        if (_listOfRedundantHS.insert( i ).second == true) {
#ifdef DEBUG
          std::cout << "Redundant " << i;
          std::cout << std::endl;
#endif
        }
      }
    }}
 
 
    // Now we've checked that half-spaces have the correct minimum number of generators we need to 
    // perform a more complicated task. Run the list to see if a set of generators can be included
    // into another one. In this case, and if we don't deal with flat polytopes or polyhedra,
    // it means that the current half-space is redundant.
    {
      unsigned int i1,i2;
      std::vector< std::set< unsigned int > >::const_iterator it_1, it_2;
      for (it_1=_allGenPerLinearConstraint.begin(), i1=0; it_1!=_allGenPerLinearConstraint.end(); ++it_1, ++i1) {
        for (it_2=_allGenPerLinearConstraint.begin(), i2=0; it_2!=_allGenPerLinearConstraint.end(); ++it_2, ++i2) {
          if (it_1!=it_2 && 
              _isHalfSpace[i1] == true && _isHalfSpace[i2] == true &&
              _listOfRedundantHS.find(i1)==_listOfRedundantHS.end() && _listOfRedundantHS.find(i2)==_listOfRedundantHS.end()) {
            if (it_1->size() > it_2->size()) {
              if (std::includes(it_1->begin(), it_1->end(), it_2->begin(), it_2->end())) {
                // We have it_2 \subset it_1.
                  _listOfRedundantHS.insert(i2);
                  unhookThisRedundantHalfSpace(listOfGen, i2, numberOfGeneratorsPerFacetWithoutHyp);
#ifdef DEBUG
                  std::cout << "New redundant facet: " << i2 << std::endl;
#endif
              }
            }
            else if (it_1->size() < it_2->size()) {
              if (std::includes(it_2->begin(), it_2->end(), it_1->begin(), it_1->end())) {
                // We have it_1 \subset it_2.
                _listOfRedundantHS.insert(i1);
                unhookThisRedundantHalfSpace(listOfGen, i1, numberOfGeneratorsPerFacetWithoutHyp);
#ifdef DEBUG
                std::cout << "New redundant facet: " << i1 << std::endl;
#endif
              }
            }
            else {
              // Unnecessary check for equality
            }
          }
        }
      }
    }
 
  }
 
  void unhookRedundantHalfSpaces(POLYHEDRON poly) {
    //fillListOfRedundantHS( poly );
    // Unhook redundant half-spaces from the vertices.
    unsigned int hsNb;
    if (!_listOfRedundantHS.empty()) {
      constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS(poly->getListOfHalfSpaces());
      {for (iteHS.begin(), hsNb=0; iteHS.end()!=true; iteHS.next(), ++hsNb) {
        if (_listOfRedundantHS.find(hsNb) != _listOfRedundantHS.end()) {
          // We have a genuine redundant half-space.
          constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(poly->getListOfGenerators());
          {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
            boost::shared_ptr<Generator_Rn> gen = iteGN.current();
            {for (int j=gen->numberOfFacets()-1; j>=0; j--) {
              if (iteHS.current() == gen->getFacet(j))
                gen->removeFacet(j);
            }}
          }}
        }
      }}
    }
    // Add to the list of redundant half-spaces the ones which are not
    // referenced at all by vertices from the beginning of the process.
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS2(poly->getListOfHalfSpaces());
    {for (iteHS2.begin(), hsNb=0; iteHS2.end()!=true; iteHS2.next(), ++hsNb) {
      if (_listOfRedundantHS.find(hsNb) == _listOfRedundantHS.end()) {
        //std::cout << "_numberOfVerticesPerLinearConstraint[] = " << _numberOfVerticesPerLinearConstraint[hsNb] << std::endl;
        if (_allGenPerLinearConstraint[hsNb].size() == 0) {
          // The half-space has not been found so it's not in use at all and is redundant.
          _listOfRedundantHS.insert(hsNb);
        }
      }
    }}
    if (!_listOfRedundantHS.empty()) {
#ifdef DEBUG
      std::cout << "Remove " << _listOfRedundantHS.size() << " redundant facets: ";
#endif
      std::vector< unsigned int > HS2Remove;
      std::set< unsigned int >::const_iterator it;
      {for (it=_listOfRedundantHS.begin(); it!=_listOfRedundantHS.end(); ++it) {
#ifdef DEBUG
        std::cout << *it << " ";
#endif
        HS2Remove.push_back(*it);
      }}
      std::sort(HS2Remove.begin(), HS2Remove.end(), std::greater<unsigned int>());
      std::vector< unsigned int >::const_iterator it2;
      {for (it2=HS2Remove.begin(); it2!=HS2Remove.end(); ++it2) {
        // Do not forget to begin by the bigger to remove as we could change the ordering.
        poly->removeHalfSpace(*it2);
      }}
#ifdef DEBUG
      std::cout << std::endl;
#endif
    }
  }
 
  void unhookRedundantGenerators(POLYHEDRON poly) {
    unsigned int minNbFacetsPerGen = poly->neigbourhoodCondition() + 1;
    std::set< boost::shared_ptr<Generator_Rn> > setToRemove;
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(poly->getListOfGenerators());
    {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
      //std::cout << "gen->numberOfFacets() = " << iteGN.current()->numberOfFacets();
      if (iteGN.current()->numberOfFacets() < minNbFacetsPerGen) {
        setToRemove.insert( iteGN.current() );
        //std::cout << " removed! ";
      }
    }}
    if (!setToRemove.empty())
      poly->removeGenerators(setToRemove);
    //std::cout << std::endl;
  }
 
  // Mark as UNCHANGED, CREATED or DELETED the half-spaces going through the double description process.
  void markHdescription(TrackingOperatorToResult& trackerHdesc, unsigned int truncationStep) {
    trackerHdesc.setNumbersOfEntities(_numberOfProcessedLinearConstraints, _numberOfProcessedLinearConstraints);
    // Before truncation step they are UNCHANGED by default, then they are CREATED by default.
    for (unsigned int i=0; i<truncationStep; ++i)
      trackerHdesc.setOperatorEntityAsUnchanged(i);
    for (unsigned int i=truncationStep; i<_numberOfProcessedLinearConstraints; ++i)
      trackerHdesc.setResultEntityAsCreated(i);
    std::set< unsigned int >::const_iterator iteRED=_listOfRedundantHS.begin();
    for ( ; iteRED!=_listOfRedundantHS.end(); ++iteRED)
      trackerHdesc.setOperatorEntityAsDeleted(*iteRED);
  }
 
  void dumpListOfRedundantHS(POLYHEDRON poly, std::ostream &this_stream) {
    unsigned int RnDIM=poly->dimension();
    this_stream << "List of redundant half-spaces :" << std::endl;
    unsigned int hsNb;
    if (!_listOfRedundantHS.empty()) {
      constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS(poly->getListOfHalfSpaces());
      {for (iteHS.begin(), hsNb=0; iteHS.end()!=true; iteHS.next(), ++hsNb) {
        if (_listOfRedundantHS.find(hsNb) != _listOfRedundantHS.end()) {
          boost::shared_ptr<HalfSpace_Rn> currentHalfSpace = iteHS.current();
          this_stream << "H = (" << currentHalfSpace->getConstant() << ", ";
          {for (unsigned int ii=0; ii<RnDIM; ii++) {
            this_stream << currentHalfSpace->getCoefficient(ii);
            if (ii != RnDIM-1)
              this_stream << ", ";
          }}
          this_stream << ")" << std::endl;
        }
      }}
    }
  }
 
  void dumpSD(std::ostream &this_stream) {
    std::vector< std::set< unsigned int > >::const_iterator iteGHS;
    this_stream << "*** number(HalfSpace_Rn0) => {Generator_Rn0*, Generator_Rn1*, Generator_Rn2*, ...}" << std::endl;
    unsigned int counter=0;
    {for (iteGHS=_allGenPerLinearConstraint.begin(); iteGHS!=_allGenPerLinearConstraint.end(); ++iteGHS) {
      this_stream << counter <<  " => { ";
      std::copy(iteGHS->begin(), iteGHS->end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
      this_stream << "}" << std::endl;
      ++counter;
    }}
    //this_stream << "*** number of generators per half-spaces : ";
    //std::copy(_allGenPerLinearConstraint.begin(), _allGenPerLinearConstraint.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
    this_stream << std::endl;
    this_stream << "*** list of redundant half-spaces : ";
    std::copy(_listOfRedundantHS.begin(), _listOfRedundantHS.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
    this_stream << std::endl;
  }
 
 protected:
  // Deal with the number of generators per half-space.
  // number(HalfSpace_Rn0) => {Generator_Rn0*, Generator_Rn1*, Generator_Rn2*, ...}
  std::vector< std::set< unsigned int > > _allGenPerLinearConstraint;
  std::vector< bool > _isHalfSpace;
  std::set< unsigned int > _listOfRedundantHS;
  unsigned int _numberOfProcessedLinearConstraints;
  unsigned int _numberOfProcessedHyperplanes;
  unsigned int _neigbourhoodCondition;
};
 
#endif