Polytope_Rn.cpp

Go to the documentation of this file.

// politopix allows to make computations on polytopes such as finding vertices, intersecting, Minkowski sums, ...
//     Copyright (C) 2011-2016 : 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)
 
#include <iostream>
#include <sstream>
#include <vector>
#include <math.h>
#include "Rn.h"
#include "Polytope_Rn.h"
#include "IO_Polytope.h"
#include "DoubleDescription_Rn.h"
#include "PolyhedralAlgorithms_Rn.h"
 
 
bool Polytope_Rn::checkEdges() const {
  // Check all edges, built from neighbor generators and topological data, are on frontier and not inside the polytope.
  std::cout << "Check edges.......... ";
  unsigned int RnDIM= this->dimension();
  bool checkEdgesOK = true;
  {for (unsigned int i=0; i<numberOfGenerators(); i++) {
    boost::shared_ptr<Generator_Rn> V1 = _listOfGenerators[i];
    for (unsigned int j=0; j<numberOfGenerators(); j++) {
      boost::shared_ptr<Generator_Rn> V2 = _listOfGenerators[j];
      std::vector< boost::shared_ptr<HalfSpace_Rn> > commonFacets;
      if (V1 != V2 && checkNeighbours(V1, V2, commonFacets) == true) {
        // Build the edge with V1 & V2 and check it is inside all facets.
        std::vector<double> edge(RnDIM);
        for (unsigned int k=0; k<RnDIM; k++)
          edge[k] = V1->getCoordinate(k) - V2->getCoordinate(k);
        unsigned int countFacets = 0;
        double scalarProduct = 0.;
        for (unsigned int l=0; l<commonFacets.size(); l++) {
          for (unsigned int n=0; n<RnDIM; n++)
            scalarProduct = scalarProduct + edge[n]*commonFacets[l]->getCoefficient(n);
          if (scalarProduct < -Rn::getTolerance() || scalarProduct > Rn::getTolerance()) {
            std::cout << "\t### Edge for generators " << i << " and " << j << " not on facet (";
            {for (unsigned int jj=0; jj<RnDIM; jj++) {
               std::cout << commonFacets[l]->getCoefficient(jj) << " ";
               if (jj == RnDIM-1)
                 std::cout << commonFacets[l]->getSideAsText() << " " << commonFacets[l]->getConstant();
            }}
            std::cout << ") " << std::endl;
            checkEdgesOK = false;
          }
          else
            countFacets++;
        }
        // For polytopes only.
        if (countFacets < RnDIM-1) {
          std::cout << "\t### Edge for generators " << i << " and " << j << " has not enough facets" << std::endl;
          checkEdgesOK = false;
        }
      }
    }
  }}
  if (checkEdgesOK == true)
    std::cout << "OK" << std::endl;
  return checkEdgesOK;
}
 
void Polytope_Rn::createBoundingSimplex(double M) {
  double MIN =-M;
  double MAX = (2*dimension()-1)*M;
  unsigned int dim = dimension();
  unsigned int currentNumberOfFacets = (unsigned int)dimension()+1;
  unsigned int currentNumberOfGenerators = (unsigned int)dimension()+1;
  //std::cout << "Number of generators = " << currentNumberOfGenerators << std::endl;
 
  // Fill the data structure with blank Generators, 4 in dim 2, 8 in dim 3, ...
  {for (unsigned int vtx_count=0; vtx_count<currentNumberOfGenerators; vtx_count++) {
    boost::shared_ptr<Generator_Rn> VX;
    VX.reset(new Generator_Rn(dim));
    addGenerator(VX);
  }}
  // Create the simplex generators (MIN, MIN, ..., MAX, ..., MIN) where MAX
  // moves from 0 to (n-1) and the last vertex is (MIN, ..., MIN).
  {for (unsigned int coord_count=0; coord_count<dim; coord_count++) {
    {for (unsigned int vtx_count=0; vtx_count<currentNumberOfGenerators; vtx_count++) {
      //vtx_count==currentNumberOfGenerators-1 => this vertex is the corner i.e. (MIN, ..., MIN)
      if (vtx_count==coord_count && vtx_count<currentNumberOfGenerators-1)
        getGenerator(vtx_count)->setCoordinate(coord_count, MAX);
      else
        getGenerator(vtx_count)->setCoordinate(coord_count, MIN);
    }}
  }}
  // Create the simplex array of facets excepted the oblic one, for the moment.
  std::vector< boost::shared_ptr<HalfSpace_Rn> > supportFacets;
  // M + xi >= 0.
  {for (unsigned int facet_count=0; facet_count<currentNumberOfFacets-1; facet_count++) {
    boost::shared_ptr<HalfSpace_Rn> HalfSp(new HalfSpace_Rn(dim));
    {for (unsigned int coord_count=0; coord_count<dim; coord_count++) {
      if (facet_count == coord_count)
        HalfSp->setCoefficient(coord_count, 1.);
      else
        HalfSp->setCoefficient(coord_count, 0.);
    }}
    HalfSp->setConstant(M);
    supportFacets.push_back(HalfSp);
    // The first generator is the corner so it has
    // to contain all the perpendicular half-spaces.
    getGenerator(currentNumberOfGenerators-1)->setFacet(HalfSp);
  }}
  // Oblic half-space (Rn::getDimension()*M, -1., ..., -1.)
  boost::shared_ptr<HalfSpace_Rn> oblicHalfSp(new HalfSpace_Rn(dim));
  {for (unsigned int coord_count=0; coord_count<dim; coord_count++) {
    oblicHalfSp->setCoefficient(coord_count, -1.);
  }}
  oblicHalfSp->setConstant(dimension()*M);
  supportFacets.push_back(oblicHalfSp);
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteH(getListOfHalfSpaces());
  {for (iteH.begin(); iteH.end()!=true; iteH.next()) {
    supportFacets.push_back(iteH.current());
  }}
  _listOfHalfSpaces.clear();
  std::vector< boost::shared_ptr<HalfSpace_Rn> >::const_iterator itF;
  {for (itF=supportFacets.begin(); itF!=supportFacets.end();++itF) {
      _listOfHalfSpaces.push_back(*itF);
  }}
  // Make connections between the generators excepted the corner (already treated).
  {for (unsigned int vtx_count1=0; vtx_count1<currentNumberOfGenerators-1; vtx_count1++) {
    boost::shared_ptr<Generator_Rn> VX1 = getGenerator(vtx_count1);
    VX1->setFacet(oblicHalfSp);
    {for (unsigned int coord_count=0; coord_count<dim; coord_count++) {
      if (VX1->getCoordinate(coord_count) == MIN) {
        VX1->setFacet(_listOfHalfSpaces[coord_count]);
      }
    }}
  }}
  //std::string FileOut("simplex.ptop");
  //boost::shared_ptr<Polytope_Rn> A(this);
  //IO_Polytope::save(FileOut, A);
  //exit(0);
}
 
void Polytope_Rn::createBoundingBox(double M) {
  double MAX = M;
  unsigned int dim = dimension();
  unsigned int currentNumberOfGenerators = (unsigned int)pow(static_cast<double>(2), static_cast<int>(dim));
  //std::cout << "Number of generators = " << currentNumberOfGenerators << std::endl;
 
  // Fill the data structure with blank Generators, 4 in dim 2, 8 in dim 3, ...
  for (unsigned int vtx_count=0; vtx_count<currentNumberOfGenerators; vtx_count++) {
    boost::shared_ptr<Generator_Rn> VX;
    VX.reset(new Generator_Rn(dim));
    addGenerator(VX);
  }
  // Create the cube generators.
  for (unsigned int coord_count=0; coord_count<dim; coord_count++) {
    int counter = 0;
    for (unsigned int vtx_count=0; vtx_count<currentNumberOfGenerators; vtx_count++) {
      if (counter >= pow(static_cast<double>(2), static_cast<int>(coord_count))) {
        MAX = -1. * MAX;
        counter = 0;
      }
      getGenerator(vtx_count)->setCoordinate(coord_count, MAX);
      counter++;
    }
  }
  // Create the cube array of facets.
  std::vector< boost::shared_ptr<HalfSpace_Rn> > supportFacets;
  // xi < M
  for (unsigned int facet_count=0; facet_count<dim; facet_count++) {
    boost::shared_ptr<HalfSpace_Rn> HalfSp(new HalfSpace_Rn(dim));
    for (unsigned int coord_count=0; coord_count<dim; coord_count++) {
      if (facet_count == coord_count)
        HalfSp->setCoefficient(coord_count, -1.);
      else
        HalfSp->setCoefficient(coord_count, 0.);
    }
    HalfSp->setConstant(M);
    supportFacets.push_back(HalfSp);
  }
  // xi > -M
  for (unsigned int facet_count=0; facet_count<dim; facet_count++) {
    boost::shared_ptr<HalfSpace_Rn> HalfSp(new HalfSpace_Rn(dim));
    for (unsigned int coord_count=0; coord_count<dim; coord_count++) {
      if (facet_count == coord_count)
        HalfSp->setCoefficient(coord_count, 1.);
      else
        HalfSp->setCoefficient(coord_count, 0.);
    }
    HalfSp->setConstant(M);
    supportFacets.push_back(HalfSp);
  }
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteH(getListOfHalfSpaces());
  {for (iteH.begin(); iteH.end()!=true; iteH.next()) {
    supportFacets.push_back(iteH.current());
  }}
  _listOfHalfSpaces.clear();
  std::vector< boost::shared_ptr<HalfSpace_Rn> >::const_iterator itF;
  {for (itF=supportFacets.begin(); itF!=supportFacets.end();++itF) {
      _listOfHalfSpaces.push_back(*itF);
  }}
  // Make connections between the generators.
  for (unsigned int vtx_count1=0; vtx_count1<currentNumberOfGenerators; vtx_count1++) {
    boost::shared_ptr<Generator_Rn> VX1 = getGenerator(vtx_count1);
    for (unsigned int coord_count=0; coord_count<dim; coord_count++) {
      if (VX1->getCoordinate(coord_count) == M)
        VX1->setFacet(_listOfHalfSpaces[coord_count]);
      else if (VX1->getCoordinate(coord_count) == -M)
        VX1->setFacet(_listOfHalfSpaces[coord_count+dim]);
    }
  }
}
 
double Polytope_Rn::createTruncatedGenerator(
  const boost::shared_ptr<Generator_Rn_SD>& out,
  const boost::shared_ptr<Generator_Rn_SD>& in,
  boost::shared_ptr<Generator_Rn_SD> newV, double ay, double az, double b) const {
  // Formula for polytopes.
  double t = (-b-az) / (ay-az);
  //for (unsigned int k=0; k<Rn::getDimension(); k++) {
    //double yi = out->getCoordinate(k);
    //double zi =  in->getCoordinate(k);
    //newV->setCoordinate(k, zi+(yi-zi)*t);
  //}
  newV->makeDiff(out, in);
  newV->makeCoefSum(in, newV, 1., t);
  return t;
}
 
bool Polytope_Rn::checkEqualityOfVertices(const boost::shared_ptr<Polytope_Rn>& B, bool printOnScreen) const {
  double RnDIM = Rn::getDimension();
  double TOL2 = Rn::getTolerance()*Rn::getTolerance();
  std::vector<bool> vtxOfA(numberOfGenerators());
  std::vector<bool> vtxOfB(B->numberOfGenerators());
  for (unsigned int i=0; i<vtxOfA.size(); ++i)
    vtxOfA[i] = false;
  for (unsigned int i=0; i<vtxOfB.size(); ++i)
    vtxOfB[i] = false;
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN_A(getListOfGenerators());
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN_B(B->getListOfGenerators());
  {for (iteGN_A.begin(); iteGN_A.end()!=true; iteGN_A.next()) {
    {for (iteGN_B.begin(); iteGN_B.end()!=true; iteGN_B.next()) {
      if (iteGN_A.current()->isEqual2(iteGN_B.current(), RnDIM, TOL2)) {
        vtxOfA[ iteGN_A.currentIteratorNumber() ] = true;
        vtxOfB[ iteGN_B.currentIteratorNumber() ] = true;
      }
    }}
  }}
  bool AinsideB = true, BinsideA = true;
  {for (unsigned int i=0; i<vtxOfA.size(); ++i) {
    if (vtxOfA[i] == false) {
      AinsideB = false;
      break;
    }
  }}
  {for (unsigned int i=0; i<vtxOfB.size(); ++i) {
    if (vtxOfB[i] == false) {
      BinsideA = false;
      break;
    }
  }}
  if (printOnScreen == true) {
    if (AinsideB == true)
      std::cout << "The set of vertices of A in included in the set of vertices of B." << std::endl;
    if (BinsideA == true)
      std::cout << "The set of vertices of B in included in the set of vertices of A." << std::endl;
  }
 
  return (AinsideB && BinsideA);
}
 
boost::shared_ptr<PolyhedralCone_Rn> Polytope_Rn::getPrimalCone(unsigned int u, const std::vector<unsigned int>& allNgb) const {
  if (u >= numberOfGenerators())
    throw std::out_of_range("Polytope_Rn::getPrimalCone() index out of range");
  unsigned int RnDIM=dimension();
  boost::shared_ptr<Generator_Rn> vx1 = getGenerator(u);
  boost::shared_ptr<PolyhedralCone_Rn> primalCone;
  // Split the polytope into its primal polyhedral cones //
  primalCone.reset(new PolyhedralCone_Rn());
  // Insert all half-spaces connected to the current vertex into the primal cone.
  for (unsigned int i=0; i<vx1->numberOfFacets(); i++) {
    primalCone->addHalfSpace(vx1->getFacet(i));
  }
  {
    std::vector<unsigned int>::const_iterator ngb;
    for (ngb=allNgb.begin(); ngb!=allNgb.end(); ngb++) {
      boost::shared_ptr<Generator_Rn> vx2 = getGenerator(*ngb);
      // Build the edge.
      boost::shared_ptr<Generator_Rn> edge12(new Generator_Rn(RnDIM));
      edge12->makeDiff(vx2, vx1);
      // Build the facet list.
      std::vector< boost::shared_ptr<HalfSpace_Rn> > commonFacets;
      // Get facets.
      checkNeighbours(vx1, vx2, commonFacets);
      for (unsigned int j=0; j<commonFacets.size(); ++j) {
        boost::shared_ptr<HalfSpace_Rn> Fj = commonFacets[j];
        edge12->setFacet(Fj);
      }
      // Insert the edge into the primal cone.
      primalCone->addGenerator(edge12);
    }
  }
  return primalCone;
}
 
boost::shared_ptr<Polytope_Rn> Polytope_Rn::sum(const boost::shared_ptr<Polytope_Rn>& pol2sum) {
  boost::shared_ptr<Polytope_Rn> currentSum(new Polytope_Rn());
  boost::shared_ptr<Polytope_Rn> operand1(this);
  MinkowskiSum Ope(operand1, pol2sum, currentSum);
  currentSum->checkTopologyAndGeometry();
  return currentSum;
}
 
boost::shared_ptr<Polytope_Rn> Polytope_Rn::intersect(const boost::shared_ptr<Polytope_Rn>& pol2intersect) {
  boost::shared_ptr<Polytope_Rn> C(new Polytope_Rn()), A(this);
 // Copy the generators of A.
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteHS0(A->getListOfGenerators());
  for (iteHS0.begin(); iteHS0.end()!=true; iteHS0.next()) {
    C->addGenerator(iteHS0.current());
  }
  // Copy the half-spaces of A.
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS1(A->getListOfHalfSpaces());
  for (iteHS1.begin(); iteHS1.end()!=true; iteHS1.next()) {
    C->addHalfSpace(iteHS1.current());
  }
  // Copy the half-spaces of B.
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS2(pol2intersect->getListOfHalfSpaces());
  for (iteHS2.begin(); iteHS2.end()!=true; iteHS2.next()) {
    C->addHalfSpace(iteHS2.current());
  }
  unsigned int truncationStep = A->numberOfHalfSpaces();
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite(C->getListOfHalfSpaces());
  //StrongRedundancyProcessing< boost::shared_ptr<Polytope_Rn> > NRP(C->neigbourhoodCondition());
  // The number of facets for a cube.
  DoubleDescription<
    boost::shared_ptr<Polytope_Rn>,
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > >
    //StrongRedundancyProcessing< boost::shared_ptr<Polytope_Rn> > >
    DD(C, lexmin_ite, truncationStep);
  C->checkTopologyAndGeometry();
  return C;
}
 
boost::shared_ptr<CappedPolytope_Rn> CappedPolytope_Rn::sum(const boost::shared_ptr<CappedPolytope_Rn>& pol2sum, double bb_size) {
  std::set< unsigned int > newCaps;
  boost::shared_ptr<Polytope_Rn> C(new Polytope_Rn());
  boost::shared_ptr<Polytope_Rn> A(this->_polytope);
  PseudoSumWithoutCaps(A, pol2sum->_polytope, C, _operandCaps, pol2sum->_operandCaps, newCaps, bb_size);
  boost::shared_ptr<CappedPolytope_Rn> finalSum(new CappedPolytope_Rn(C, newCaps) );
  return finalSum;
}
 
boost::shared_ptr<CappedPolytope_Rn> CappedPolytope_Rn::intersect(const boost::shared_ptr<CappedPolytope_Rn>& pol2intersect) {
  std::set< unsigned int > newCaps;
  boost::shared_ptr<Polytope_Rn> C(new Polytope_Rn());
  boost::shared_ptr<Polytope_Rn> A(this->_polytope);
  PseudoIntersectionWithoutCaps(A, pol2intersect->_polytope, C, _operandCaps, pol2intersect->_operandCaps, newCaps);
  boost::shared_ptr<CappedPolytope_Rn> finalSum(new CappedPolytope_Rn(C, newCaps) );
  return finalSum;
}
 
void CappedPolytope_Rn::load(const std::string& pathA, double bb_size) {
  boost::shared_ptr<Polytope_Rn> loadedPolytope(new Polytope_Rn()), cube(new Polytope_Rn());
  IO_Polytope::load(pathA, loadedPolytope);
  _polytope.reset(new Polytope_Rn());
  _operandCaps.clear();
  // Use a trick based on the size of the bounding box, if it is null consider we do not truncate and return now.
  if (bb_size == 0.) {
    // Copy the half-spaces of loadedPolytope.
    {
      constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS1(loadedPolytope->getListOfHalfSpaces());
      for (iteHS1.begin(); iteHS1.end()!=true; iteHS1.next())
        _polytope->addHalfSpace(iteHS1.current());
    }
    return;
  }
  // Create the bounding box centered on the origin.
  _polytope->createBoundingBox(bb_size);
  // Make a copy the half-spaces of the cube, after the truncation we want to know which one remained.
  {
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS1(_polytope->getListOfHalfSpaces());
    for (iteHS1.begin(); iteHS1.end()!=true; iteHS1.next())
      cube->addHalfSpace(iteHS1.current());
  }
  // Copy the half-spaces of loadedPolytope.
  {
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS1(loadedPolytope->getListOfHalfSpaces());
    for (iteHS1.begin(); iteHS1.end()!=true; iteHS1.next())
      _polytope->addHalfSpace(iteHS1.current());
  }
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite(_polytope->getListOfHalfSpaces());
  //StrongRedundancyProcessing< boost::shared_ptr<Polytope_Rn> > NRP(_polytope->neigbourhoodCondition());
  // The number of facets for a cube.
  unsigned int truncationStep = 2*_polytope->dimension();
  DoubleDescription<
    boost::shared_ptr<Polytope_Rn>,
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > >
    //StrongRedundancyProcessing< boost::shared_ptr<Polytope_Rn> > >
      DD(_polytope, lexmin_ite, truncationStep);
  // Now try to find all the cube half-space remaining inside _polytope
  {
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > allCubeHS(cube->getListOfHalfSpaces());
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > allPolHS(_polytope->getListOfHalfSpaces());
    unsigned int polFacetNb = 0;
    for (allPolHS.begin(); allPolHS.end()!=true; allPolHS.next()) {
      bool isEqual = false;
      for (allCubeHS.begin(); allCubeHS.end()!=true; allCubeHS.next()) {
        if (allCubeHS.current() == allPolHS.current())
          isEqual = true;
      }
      if (isEqual == true)
        _operandCaps.insert(polFacetNb);
      polFacetNb++;
    }
  }
}
 
void CappedPolytope_Rn::save(const std::string& pathA) {
  boost::shared_ptr<Polytope_Rn> loadedPolytope(new Polytope_Rn());
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > allPolHS(_polytope->getListOfHalfSpaces());
  for (allPolHS.begin(); allPolHS.end()!=true; allPolHS.next()) {
    if (_operandCaps.find(allPolHS.currentIteratorNumber()) == _operandCaps.end())
      loadedPolytope->addHalfSpace(allPolHS.current());
  }
  IO_Polytope::save(pathA, loadedPolytope);
}
 
int CappedPolytope_Rn::scalingFactor(double factor) {
  return TopGeomTools::scalingFactor(this->_polytope, factor);
}