Polytope_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) 2011-2020 : 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 POLYTOPE_Rn
#define POLYTOPE_Rn
 
#include <boost/shared_ptr.hpp>
#include <stdexcept>
#include <exception>
#include <iostream>
#include <vector>
#include "polito_Export.h"
#include "PolyhedralCone_Rn.h"
#include "Generator_Rn.h"
#include "HalfSpace_Rn.h"
#include "IO_Polytope.h"
 
class polito_EXPORT Polytope_Rn : public PolyhedralCone_Rn {
 
 public:
  Polytope_Rn():PolyhedralCone_Rn() {}
 
  virtual ~Polytope_Rn() {}
 
  static boost::shared_ptr<Polytope_Rn> create() { return boost::shared_ptr<Polytope_Rn>(new Polytope_Rn); }
 
  virtual bool isBounded() const {return true;}
 
  virtual unsigned int neigbourhoodCondition() const {return (dimension()-1);}
 
  virtual unsigned int numberOfGeneratorsPerFacet() const {return (dimension());}
 
  virtual void createBoundingSimplex(double M);
 
  virtual void createBoundingBox(double M);
 
  virtual bool checkEdges() const;
 
  bool checkEqualityOfVertices(const boost::shared_ptr<Polytope_Rn>& B, bool printOnScreen = false) const;
 
  virtual double createTruncatedGenerator(
      const boost::shared_ptr<Generator_Rn_SD>& in,
      const boost::shared_ptr<Generator_Rn_SD>& out,
      boost::shared_ptr<Generator_Rn_SD> newV, double ay, double az, double b=0.) const;
 
  boost::shared_ptr<PolyhedralCone_Rn> getPrimalCone(unsigned int u, const std::vector<unsigned int>& allNgb) const;
 
  boost::shared_ptr<PolyhedralCone_Rn> getPrimalCone(unsigned int i) const;
 
  boost::shared_ptr<Polytope_Rn> sum(const boost::shared_ptr<Polytope_Rn>& pol2sum) ;
 
  boost::shared_ptr<Polytope_Rn> intersect(const boost::shared_ptr<Polytope_Rn>& pol2intersect) ;
 
  bool isIncluded(const boost::shared_ptr<Polytope_Rn>& B) const {
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS_B(B->_listOfHalfSpaces);
    {for (iteHS_B.begin(); iteHS_B.end()!=true; iteHS_B.next()) {
      double halfSpaceNorm = std::inner_product(iteHS_B.current()->begin(), iteHS_B.current()->end(), iteHS_B.current()->begin(), 0.);
      halfSpaceNorm = sqrt(halfSpaceNorm);
      constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN_A(_listOfGenerators);
      {for (iteGN_A.begin(); iteGN_A.end()!=true; iteGN_A.next()) {
        HalfSpace_Rn::State currentState = checkPoint(iteGN_A.current(), iteHS_B.current(), halfSpaceNorm);
        if (currentState == HalfSpace_Rn::hs_OUT)
          return false;
      }}
    }}
    return true;
  }
 
  bool isIncluded(const boost::shared_ptr<Polytope_Rn>& B, double& minmaxDistance) const {
    bool isOut = false;
    double TOL = Rn::getTolerance();
    double maxDist = std::numeric_limits<double>::min();
    double minDist = std::numeric_limits<double>::max();
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS_B(B->_listOfHalfSpaces);
    {for (iteHS_B.begin(); iteHS_B.end()!=true; iteHS_B.next()) {
      double halfSpaceNorm = std::inner_product(iteHS_B.current()->begin(), iteHS_B.current()->end(), iteHS_B.current()->begin(), 0.);
      halfSpaceNorm = sqrt(halfSpaceNorm);
      constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN_A(_listOfGenerators);
      {for (iteGN_A.begin(); iteGN_A.end()!=true; iteGN_A.next()) {
        double scalarProduct = std::inner_product(iteGN_A.current()->begin(), iteGN_A.current()->end(), iteHS_B.current()->begin(), 0.);
        double distanceToHyperplane = (scalarProduct+iteHS_B.current()->getConstant()) / halfSpaceNorm;
        //std::cout << "thisPoint" << i << " = " << thisPoint << std::endl;
        //std::cout << "dist" << i << " = " << distanceToHyperplane << std::endl;
        if (distanceToHyperplane < -TOL) {
          isOut = true;
          // As distanceToHyperplane, in this case, is always negative, use its opposite.
          double opp = -distanceToHyperplane;
          if (opp > maxDist)
            maxDist = opp;
        }
        else {
          // Here distanceToHyperplane \in [-TOL, +inf[ so turn it into a positive number in case it was negative.
          if (distanceToHyperplane < 0.)
            distanceToHyperplane = -distanceToHyperplane;
          if (distanceToHyperplane < minDist)
            minDist = distanceToHyperplane;
        }
      }}
    }}
    // If at least one vertex is out, then exit claiming we have no inclusion.
    if (isOut == true) {
      minmaxDistance = maxDist;
      return false;
    }
    else
      minmaxDistance = minDist;
    return true;
  }
 
  std::string getName() const {return std::string("Polytope");}
 
  std::string getFileExtension() const {return std::string(".ptop");}
 
};
 
 
class polito_EXPORT CappedPolytope_Rn {
 
 public:
  CappedPolytope_Rn(boost::shared_ptr<Polytope_Rn> pol, const std::set< unsigned int >& set) {
    _operandCaps = set;
    _polytope = pol;
  }
 
  CappedPolytope_Rn() {}
 
  virtual ~CappedPolytope_Rn() {}
 
  void setPolytope(boost::shared_ptr<Polytope_Rn> pol) {_polytope = pol;}
 
  boost::shared_ptr<Polytope_Rn> getPolytope() {return _polytope;}
 
  const boost::shared_ptr<PolyhedralCone_Rn>& getPrismaticHPolyhedron() {
    _prismaticPolyhedron.reset(new PolyhedralCone_Rn());
    {for (unsigned int i=0; i<_polytope->numberOfHalfSpaces(); ++i) {
      if (_operandCaps.find(i) == _operandCaps.end())
        _prismaticPolyhedron->addHalfSpace( _polytope->getHalfSpace(i) );
    }}
    return _prismaticPolyhedron;
  }
 
  void resetConstantsOfAllHalfSpaces(double cst) {
    if (_polytope) {
      for (unsigned int i=0; i<_polytope->numberOfHalfSpaces(); ++i) {
        _polytope->getHalfSpace(i)->setConstant(cst);
      }
    }
    if (_prismaticPolyhedron) {
      for (unsigned int i=0; i<_prismaticPolyhedron->numberOfHalfSpaces(); ++i) {
        _prismaticPolyhedron->getHalfSpace(i)->setConstant(cst);
      }
    }
  }
 
  void resetConstantOfGivenHalfSpace(double cst, unsigned int fctNb) {
    if (_polytope) {
      if (fctNb >= _polytope->numberOfHalfSpaces())
        throw std::out_of_range("CappedPolytope_Rn::resetConstantOfGivenHalfSpace() polytope index out of range");
      _polytope->getHalfSpace(fctNb)->setConstant(cst);
    }
    if (_prismaticPolyhedron) {
      if (fctNb >= _prismaticPolyhedron->numberOfHalfSpaces())
        throw std::out_of_range("CappedPolytope_Rn::resetConstantOfGivenHalfSpace() prismaticPolyhedron index out of range");
      _prismaticPolyhedron->getHalfSpace(fctNb)->setConstant(cst);
    }
  }
 
  void setCappedFacets(const std::set< unsigned int >& set) {_operandCaps = set;}
 
  boost::shared_ptr<CappedPolytope_Rn> sum(const boost::shared_ptr<CappedPolytope_Rn>& pol2sum, double bb_size=1000.) ;
 
  void sum(const std::vector< boost::shared_ptr<CappedPolytope_Rn> >& allPolytopes0) {
    std::vector< boost::shared_ptr<CappedPolytope_Rn> > allPolytopes = allPolytopes0;
    {for (unsigned int i=0; i<allPolytopes.size()-1; ++i) {
      boost::shared_ptr<CappedPolytope_Rn> currentSum(new CappedPolytope_Rn());
      currentSum = allPolytopes[i]->sum(allPolytopes[i+1]);
      // Check whether the job is finish or not.
      if (i == allPolytopes.size()-2) {
        _operandCaps = currentSum->_operandCaps;
        _polytope = currentSum->_polytope;
      }
      else
        allPolytopes[i+1] = currentSum;
    }}
  }
 
  void intersect(const std::vector< boost::shared_ptr<CappedPolytope_Rn> >& allPolytopes0) {
    std::vector< boost::shared_ptr<CappedPolytope_Rn> > allPolytopes = allPolytopes0;
    {for (unsigned int i=0; i<allPolytopes.size()-1; ++i) {
      boost::shared_ptr<CappedPolytope_Rn> currentInt(new CappedPolytope_Rn());
      currentInt = allPolytopes[i]->intersect(allPolytopes[i+1]);
      // Check whether the job is finish or not.
      if (i == allPolytopes.size()-2) {
        _operandCaps = currentInt->_operandCaps;
        _polytope = currentInt->_polytope;
      }
      else
        allPolytopes[i+1] = currentInt;
    }}
  }
 
  boost::shared_ptr<CappedPolytope_Rn> intersect(const boost::shared_ptr<CappedPolytope_Rn>& pol2intersect) ;
 
  bool isIncluded(const boost::shared_ptr<CappedPolytope_Rn>& B) const { return _polytope->isIncluded(B->_polytope); }
 
  static bool isIncluded(const boost::shared_ptr<CappedPolytope_Rn>& firstBody, const boost::shared_ptr<CappedPolytope_Rn>& secondBody, double& minmaxDistance) {
    return firstBody->_polytope->isIncluded(secondBody->_polytope, minmaxDistance);
  }
 
  int scalingFactor(double factor);
 
  bool checkTopologyAndGeometry() const { return _polytope->checkTopologyAndGeometry(); }
 
  void save(const std::string& pathA);
 
  void load(const std::string& pathA, double bb_size=10000.);
 
  void dump(std::ostream &current_ostream) const {
    current_ostream << std::endl << "#CAPPED FACETS:" << std::endl;
    std::copy(_operandCaps.begin(), _operandCaps.end(), std::ostream_iterator<unsigned int>(current_ostream, " ") );
    current_ostream << std::endl;
    _polytope->dump(current_ostream);
  }
 
  static std::string getName() {return std::string("CappedPolytope");}
 
  static std::string getFileExtension() {return std::string(".ctop");}
 
 protected:
  std::set< unsigned int > _operandCaps;
  boost::shared_ptr<Polytope_Rn> _polytope;
  boost::shared_ptr<PolyhedralCone_Rn> _prismaticPolyhedron;
 
};
 
#endif // POLYTOPE_Rn