politopix  4.1.0
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
PolyhedralAlgorithms_Rn.cpp
Go to the documentation of this file.
1 // politopix allows to make computations on polytopes such as finding vertices, intersecting, Minkowski sums, ...
2 // Copyright (C) 2011-2016 : Delos Vincent
3 //
4 // This program is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU Lesser General Public License as published by
6 // the Free Software Foundation, either version 3 of the License, or
7 // (at your option) any later version.
8 //
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU Lesser General Public License for more details.
13 //
14 // You should have received a copy of the GNU Lesser General Public License
15 // along with this program. If not, see <http://www.gnu.org/licenses/>.
16 //
19 // I2M (UMR CNRS 5295 / University of Bordeaux)
20 
21 #include <iostream>
22 #include <sstream>
23 #include <vector>
24 #include <math.h>
25 #include <numeric>
26 #include <boost/numeric/ublas/matrix.hpp>
27 #include <boost/numeric/ublas/matrix_proxy.hpp>
28 #include <boost/numeric/ublas/operation.hpp>
29 #include <boost/numeric/ublas/io.hpp>
30 #include <boost/timer.hpp>
31 #include "Rn.h"
32 #include "Polytope_Rn.h"
33 #include "IO_Polytope.h"
34 #include "NormalFan_Rn.h"
35 #include "PolyhedralAlgorithms_Rn.h"
36 
37 //typedef std::vector< unsigned int > ListOfFaces;
38 
39 
40 void FaceEnumeration::Compute(const boost::shared_ptr<Polytope_Rn>& thisPolytope) {
41  FaceEnumeration FE(thisPolytope);
42  FaceEnumeration::ComputeWithFacets(thisPolytope, FE);
43  FaceEnumeration::ComputeWithVertices(thisPolytope, FE);
44 }
45 
46 void FaceEnumeration::Compute(const boost::shared_ptr<Polytope_Rn>& thisPolytope, FaceEnumeration& FaceEnum) {
47  FaceEnumeration::ComputeWithFacets(thisPolytope, FaceEnum);
48  FaceEnumeration::ComputeWithVertices(thisPolytope, FaceEnum);
49 }
50 
51 void FaceEnumeration::ComputeWithFacets(const boost::shared_ptr<Polytope_Rn>& thisPolytope, FaceEnumeration& FaceEnum) {
52  unsigned int RnDIM = Rn::getDimension();
53  FaceEnum._allFacesWithFacets.clear();
54  std::set< ListOfFaces > allPotentialFaces;
55  std::vector< ListOfFaces > thisListOfFacetsPerVertices;
56  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > constITER_GN1(thisPolytope->getListOfGenerators());
57  {for (constITER_GN1.begin(); constITER_GN1.end()!=true; constITER_GN1.next()) {
58  unsigned int NbFacets = constITER_GN1.current()->numberOfFacets();
59  ListOfFaces thisListOfFaces;
60  {for (unsigned int j=0; j<NbFacets; j++) {
61  unsigned int NbFj = thisPolytope->getHalfSpaceNumber( constITER_GN1.current()->getFacet(j) );
62  // GenPerFacet: number(HalfSpace_Rn) => {Generator_Rn_0, Generator_Rn_1, Generator_Rn_2, ...}
63  thisListOfFaces.push_back( NbFj );
64  }}
65  thisListOfFacetsPerVertices.push_back( thisListOfFaces );
66  allPotentialFaces.insert(thisListOfFaces);
67  }}
68  FaceEnum._allFacesWithFacets.push_back( thisListOfFacetsPerVertices );
69  //std::cout << "RnDIM = " << RnDIM << std::endl;
70  {for (unsigned int dimensionOfFace=1; dimensionOfFace<=RnDIM; ++dimensionOfFace) {
71  //std::cout << std::endl << "D=" << dimensionOfFace << std::endl;
72  std::vector< ListOfFaces > kp1_Faces;
73  std::set< ListOfFaces > kp1_FacesSet;
74  std::vector< ListOfFaces >& k_Faces = FaceEnum._allFacesWithFacets[dimensionOfFace-1];
75  std::vector< ListOfFaces >::iterator it1 = k_Faces.begin(), it2 = k_Faces.begin();
76  for (it1 = k_Faces.begin(); it1 != k_Faces.end(); ) {
77  // To know whether we had to move forward.
78  bool it1StepForward = false;
79  for (it2 = it1+1; it2 != k_Faces.end(); ) {
80  std::vector< unsigned int > interFace;
81  std::set_intersection(it1->begin(), it1->end(), it2->begin(), it2->end(), std::inserter(interFace, interFace.end()));
82  // Check whether this face as already been processed.
83  //std::cout << std::endl << "| it1 = { ";
84  //std::copy(it1->begin(), it1->end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
85  //std::cout << "} INTER it2 = { ";
86  //std::copy(it2->begin(), it2->end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
87  //std::cout << "} = { ";
88  //std::copy(interFace.begin(), interFace.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
89  //std::cout << "} ";
90  //std::cout << "inserted ";
91  if (interFace.empty() == false) {
92  if (interFace == *it2) {
93  it2 = k_Faces.erase(it2);
94  //std::cout << "(it2 erased) ";
95  }
96  else
97  ++it2;
98  if (interFace == *it1) {
99  it1 = k_Faces.erase(it1);
100  //std::cout << "(it1 erased) ";
101  it1StepForward = true;
102  if (it1 != k_Faces.end())
103  it2 = it1+1;
104  }
105  //std::cout << "in kp1 ";
106  // Check it was not already in.
107  if (kp1_FacesSet.insert(interFace).second == true)
108  kp1_Faces.push_back(interFace);
109  }
110  else {
111  //std::cout << "empty intersection ";
112  ++it2;
113  }
114  } // for (it2 = k_Faces.begin(); it2 != k_Faces.end(); ) {
115  // Only move forward iif we didn't do it before.
116  if (it1StepForward == false)
117  ++it1;
118  }
119  FaceEnum._allFacesWithFacets.push_back(kp1_Faces);
120  }}
121  //FaceEnum.printFacesWithFacets(std::cout);
122 }
123 
124 void FaceEnumeration::ComputeWithVertices(const boost::shared_ptr<Polytope_Rn>& thisPolytope, FaceEnumeration& FaceEnum) {
125  unsigned int nbHS = thisPolytope->numberOfHalfSpaces();
126  FaceEnum._allFacesWithVertices.clear();
127  FaceEnum._allFacesWithVertices.resize(Rn::getDimension()+1);
128  // number(HalfSpace_Rn0) => {Generator_Rn0*, Generator_Rn1*, Generator_Rn2*, ...}
129  std::vector< std::vector< unsigned int > > allGenPerHS;
130  {for (unsigned int j=0; j<nbHS; ++j) {
131  std::vector< unsigned int > genSet;
132  allGenPerHS.push_back(genSet);
133  }}
134 
135  unsigned int generatorNumber=0;
136  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(thisPolytope->getListOfGenerators());
137  {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
138  {for (unsigned int ii=0; ii<iteGN.current()->numberOfFacets(); ii++) {
139  unsigned int fctNumber = thisPolytope->getHalfSpaceNumber( iteGN.current()->getFacet(ii) );
140  allGenPerHS[fctNumber].push_back( iteGN.currentIteratorNumber() );
141  }}
142  ++generatorNumber;
143  }}
144 
145  unsigned int dim = 0;
146  std::vector< std::vector< ListOfFaces > >::const_iterator iteF;
147  for (iteF=FaceEnum._allFacesWithFacets.begin(); iteF!=FaceEnum._allFacesWithFacets.end(); ++iteF) {
148  //this_ostream << "F" << dim << ": ";
149  std::vector< ListOfFaces >::const_iterator iteFF;
150  for (iteFF=iteF->begin(); iteFF!=iteF->end(); ++iteFF) {
151  ListOfFaces::const_iterator iteFFF;
152  //this_ostream << " { ";
153  iteFFF=iteFF->begin();
154  std::vector< unsigned int > INTER = allGenPerHS[*iteFFF];
155  ++iteFFF;
156  for (; iteFFF!=iteFF->end(); ++iteFFF) {
157  //this_ostream << *iteFFF << " ";
158  std::vector< unsigned int > partial_INTER;
159  std::set_intersection(INTER.begin(), INTER.end(),
160  allGenPerHS[*iteFFF].begin(), allGenPerHS[*iteFFF].end(),
161  std::inserter(partial_INTER, partial_INTER.end()));
162  INTER = partial_INTER;
163  }
164  FaceEnum._allFacesWithVertices[dim].push_back(INTER);
165  //std::copy(INTER.begin(), INTER.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
166  //this_ostream << "} ";
167  }
168  //this_ostream << std::endl;
169  ++dim;
170  }
171 }
172 
173 void FaceEnumeration::printFacesWithFacets(std::ostream &this_ostream) const {
174  unsigned int dim = 0, allSizes = 0;
175  std::vector< std::vector< ListOfFaces > >::const_iterator iteF;
176  for (iteF=_allFacesWithFacets.begin(); iteF!=_allFacesWithFacets.end(); ++iteF) {
177  allSizes = allSizes + iteF->size();
178  }
179  // Add the empty polytope and the total polytope.
180  allSizes = allSizes + 2;
181  this_ostream << "FaceEnumeration::printFacesWithFacets, total number of elements = " << allSizes << std::endl;
182  for (iteF=_allFacesWithFacets.begin(); iteF!=_allFacesWithFacets.end(); ++iteF) {
183  this_ostream << "F" << dim << ": ";
184  std::vector< ListOfFaces >::const_iterator iteFF;
185  for (iteFF=iteF->begin(); iteFF!=iteF->end(); ++iteFF) {
186  ListOfFaces::const_iterator iteFFF;
187  this_ostream << " { ";
188  for (iteFFF=iteFF->begin(); iteFFF!=iteFF->end(); ++iteFFF) {
189  this_ostream << *iteFFF << " ";
190  }
191  this_ostream << "} ";
192  }
193  this_ostream << std::endl;
194  ++dim;
195  }
196 }
197 
198 void FaceEnumeration::printFacesWithVertices(std::ostream &this_ostream) const {
199  unsigned int dim = 0, allSizes = 0;
200  std::vector< std::vector< ListOfFaces > >::const_iterator iteF;
201  for (iteF=_allFacesWithVertices.begin(); iteF!=_allFacesWithVertices.end(); ++iteF) {
202  allSizes = allSizes + iteF->size();
203  }
204  // Add the empty polytope and the total polytope.
205  allSizes = allSizes + 2;
206  this_ostream << "FaceEnumeration::printFacesWithVertices, total number of elements = " << allSizes << std::endl;
207  for (iteF=_allFacesWithVertices.begin(); iteF!=_allFacesWithVertices.end(); ++iteF) {
208  this_ostream << "F" << dim << ": ";
209  std::vector< ListOfFaces >::const_iterator iteFF;
210  for (iteFF=iteF->begin(); iteFF!=iteF->end(); ++iteFF) {
211  ListOfFaces::const_iterator iteFFF;
212  this_ostream << " { ";
213  for (iteFFF=iteFF->begin(); iteFFF!=iteFF->end(); ++iteFFF) {
214  this_ostream << *iteFFF << " ";
215  }
216  this_ostream << "} ";
217  }
218  this_ostream << std::endl;
219  ++dim;
220  }
221 }
222 
223 void FaceEnumeration::printFacesWithVerticesToSage(std::ostream &this_ostream) const {
224  unsigned int dim = 0, allSizes = 0;
225  std::vector< std::vector< ListOfFaces > >::const_iterator iteF;
226  for (iteF=_allFacesWithVertices.begin(); iteF!=_allFacesWithVertices.end(); ++iteF) {
227  allSizes = allSizes + iteF->size();
228  }
229  // Add the empty polytope and the total polytope.
230  allSizes = allSizes + 2;
231  this_ostream << "FaceEnumeration::printFacesWithVerticesToSage, total number of elements = " << allSizes << std::endl;
232  for (iteF=_allFacesWithVertices.begin(); iteF!=_allFacesWithVertices.end(); ++iteF) {
233  this_ostream << "F" << dim << ": " << std::endl;
234  std::vector< ListOfFaces >::const_iterator iteFF;
235  //std::cout << "ok "; //@
236  std::vector< ListOfFaces >::const_iterator stopFF=iteF->end()-1;
237  //std::cout << ", stopFF=" << std::endl; //@
238  for (iteFF=iteF->begin(); iteFF!=stopFF && iteFF!=iteF->end(); ++iteFF) {
239  ListOfFaces::const_iterator iteFFF, stopFFF=iteFF->end()-1;
240  this_ostream << "(";
241  for (iteFFF=iteFF->begin(); iteFFF!=stopFFF; ++iteFFF) {
242  this_ostream << *iteFFF << ", ";
243  }
244  this_ostream << *iteFFF << "), ";
245  }
246  //std::cout << "ok2 " << std::endl; //@
247  // Print the last one.
248  if (iteFF!=iteF->end()) {
249  ListOfFaces::const_iterator iteFFF, stopFFF=iteFF->end()-1;
250  this_ostream << "(";
251  for (iteFFF=iteFF->begin(); iteFFF!=stopFFF; ++iteFFF) {
252  this_ostream << *iteFFF << ", ";
253  }
254  this_ostream << *iteFFF << ") ";
255  }
256  //std::cout << "ok3 " << std::endl; //@
257  this_ostream << std::endl;
258  ++dim;
259  }
260 }
261 
262 void FaceEnumeration::load(const std::string& filename, std::vector< std::vector< ListOfFaces > >& latt)
263  throw (std::ios_base::failure) {
264  std::ifstream file(filename.c_str(), std::ifstream::in);
265  if (!file) {
266  std::string s("Unable to open ");
267  s += filename;
268  s += "\n";
269  throw std::ios_base::failure(s);
270  }
271  FaceEnumeration::load(file, latt);
272  file.close();
273 }
274 
275 void FaceEnumeration::save(const std::string& filename, const std::vector< std::vector< ListOfFaces > >& latt) {
276  std::ofstream file(filename.c_str());
277  if (!file) {
278  std::string s("Unable to open ");
279  s += filename;
280  s += "\n";
281  throw std::ios_base::failure(s);
282  }
283  FaceEnumeration::save(file, latt);
284  file.close();
285 }
286 
287 void FaceEnumeration::load(std::istream& this_stream, std::vector< std::vector< ListOfFaces > >& latt)
288  throw (std::out_of_range){
289  std::string line;
290  // Read the comment line.
291  std::getline(this_stream, line);
292  // Get the 2 integers : space dimension, number of faces.
293  std::getline(this_stream, line);
294  std::istringstream iline(line);
295  unsigned int spaceDimension, totalNumberOfFaces;
296  iline >> spaceDimension;
297  iline >> totalNumberOfFaces;
298  latt.resize(spaceDimension+1);
299  unsigned nbReadFaces = 0;
300  while (nbReadFaces != totalNumberOfFaces) {
301  ListOfFaces thisOne;
302  std::getline(this_stream, line);
303  std::istringstream iline3(line);
304  unsigned int dimensionOfFace, numberOfVtx, val;
305  iline3 >> dimensionOfFace;
306  iline3 >> numberOfVtx;
307  if (dimensionOfFace > spaceDimension) {
308  std::string errorMsg("FaceEnumeration::load() wrong face dimension");
309  throw std::out_of_range(errorMsg);
310  }
311  // dimensionOfFace numberOfVtx a0 a1 ... an
312  for (unsigned int count=0; count<numberOfVtx; ++count) {
313  iline3 >> val;
314  thisOne.push_back(val);
315  }
316  latt[dimensionOfFace].push_back(thisOne);
317  ++nbReadFaces;
318  }
319 }
320 
321 void FaceEnumeration::save(std::ostream& this_stream, const std::vector< std::vector< ListOfFaces > >& latt) {
322  unsigned int allSizes = 0;
323  std::vector< std::vector< ListOfFaces > >::const_iterator iteF;
324  for (iteF=latt.begin(); iteF!=latt.end(); ++iteF) {
325  allSizes = allSizes + iteF->size();
326  }
327  // Don't add the empty polytope and the total polytope.
328  //allSizes = allSizes + 2;
329  //this_ostream << "FaceEnumeration::printFacesWithFacets, total number of elements = " << allSizes << std::endl;
330  unsigned int spaceDimension = Rn::getDimension();
331  unsigned int totalNumberOfFaces = allSizes;
332  this_stream << "# Dimension LatticeSize" << std::endl;
333  this_stream << spaceDimension << " ";
334  this_stream << totalNumberOfFaces << std::endl;
335  unsigned dimReadFaces = 0;
336  for (iteF=latt.begin(); iteF!=latt.end(); ++iteF) {
337  std::vector< ListOfFaces >::const_iterator iteFF;
338  for (iteFF=iteF->begin(); iteFF!=iteF->end(); ++iteFF) {
339  this_stream << dimReadFaces << " " << iteFF->size() << " ";
340  ListOfFaces::const_iterator iteFFF;
341  for (iteFFF=iteFF->begin(); iteFFF!=iteFF->end(); ++iteFFF)
342  this_stream << *iteFFF << " ";
343  this_stream << std::endl;
344  }
345  ++dimReadFaces;
346  }
347 }
348 
349 MinkowskiSum::MinkowskiSum(const std::vector< boost::shared_ptr<Polytope_Rn> >& arrayOfPolytopes, boost::shared_ptr<Polytope_Rn>& C) {
350 
351  boost::shared_ptr<Polytope_Rn> currentSum = *(arrayOfPolytopes.begin());
352  boost::shared_ptr<Polytope_Rn> A = currentSum;
353  std::vector< boost::shared_ptr<Generator_Rn> > listOfGenerators_A, listOfGenerators_B, listOfGenerators_C;
354  std::vector< boost::shared_ptr<PolyhedralCone_Rn> > listOfDualCones_A, listOfDualCones_B, listOfDualCones_C;
355  // Build the list of vertices and dual cones for the first operand.
356  {for (unsigned int i=0; i<A->numberOfGenerators(); ++i) {
357  listOfGenerators_A.push_back(A->getGenerator(i));
358  const boost::shared_ptr<PolyhedralCone_Rn>& cone_a = A->getPrimalCone( A->getGenerator(i) );
359  listOfDualCones_A.push_back( cone_a->computeDualPolyhedralCone() );
360  }}
361  {
362  std::vector< boost::shared_ptr<Polytope_Rn> >::const_iterator itePol;
363  for (itePol = arrayOfPolytopes.begin()+1; itePol != arrayOfPolytopes.end(); ++itePol) {
364  boost::shared_ptr<Polytope_Rn> B = *itePol;
365  {for (unsigned int i=0; i<B->numberOfGenerators(); ++i) {
366  listOfGenerators_B.push_back(B->getGenerator(i));
367  const boost::shared_ptr<PolyhedralCone_Rn>& cone_b = B->getPrimalCone( B->getGenerator(i) );
368  listOfDualCones_B.push_back( cone_b->computeDualPolyhedralCone() );
369  }}
370  _A2C.resize(listOfGenerators_A.size());
371  _B2C.resize(listOfGenerators_B.size());
372  _MinkowskiDecompositionOK.clear();
373  _MinkowskiDecomposition.clear();
374  computeCommonRefinement(listOfGenerators_A, listOfGenerators_B, listOfGenerators_C, listOfDualCones_A, listOfDualCones_B, listOfDualCones_C);
375  listOfDualCones_A.clear();
376  listOfDualCones_A.swap(listOfDualCones_C);
377  listOfGenerators_A.clear();
378  listOfGenerators_A.swap(listOfGenerators_C);
379  listOfGenerators_C.clear();
380  listOfGenerators_B.clear();
381  listOfDualCones_C.clear();
382  listOfDualCones_B.clear();
383  }
384  }
385  // The global common refinement is in A now.
386  MinkowskiSum globalSum;
387  globalSum._sum = C;
388  globalSum._A2C = _A2C;
389  globalSum._B2C = _B2C;
390  globalSum._NF_Cones = listOfDualCones_A;
391  globalSum._NF_Vertices = listOfGenerators_A;
392  globalSum._MinkowskiDecomposition = _MinkowskiDecomposition;
393  globalSum._MinkowskiDecompositionOK = _MinkowskiDecompositionOK;
394  // Build the final sum.
395  globalSum.processNormalFan2();
396 }
397 
398 void MinkowskiSum::computeCommonRefinement(
399  const std::vector< boost::shared_ptr<Generator_Rn> >& listOfGenerators_A,
400  const std::vector< boost::shared_ptr<Generator_Rn> >& listOfGenerators_B,
401  std::vector< boost::shared_ptr<Generator_Rn> >& listOfGenerators_C,
402  const std::vector< boost::shared_ptr<PolyhedralCone_Rn> >& listOfDualCones_A,
403  const std::vector< boost::shared_ptr<PolyhedralCone_Rn> >& listOfDualCones_B,
404  std::vector< boost::shared_ptr<PolyhedralCone_Rn> >& listOfDualCones_C) throw (std::domain_error) {
405  if (listOfGenerators_A.size()==0 || listOfGenerators_B.size()==0)
406  throw std::domain_error("MinkowskiSum::computeCommonRefinement() cannot work without vertices.");
407  unsigned int RnDIM = Rn::getDimension();
408  std::vector< boost::shared_ptr<Generator_Rn> >::const_iterator LVX_A, LVX_B;
409  std::vector< boost::shared_ptr<PolyhedralCone_Rn> >::const_iterator LPC_A, LPC_B;
410  unsigned int minkowskiVertexNumber=0;
411  {for (LPC_B=listOfDualCones_B.begin(), LVX_B=listOfGenerators_B.begin(); LPC_B!=listOfDualCones_B.end(); ++LPC_B, ++LVX_B) {
412  boost::shared_ptr<PolyhedralCone_Rn> Cb = *(LPC_B);
413  {for (LPC_A=listOfDualCones_A.begin(), LVX_A=listOfGenerators_A.begin(); LPC_A!=listOfDualCones_A.end(); ++LPC_A, ++LVX_A) {
414  unsigned int a_i = LPC_A-listOfDualCones_A.begin();
415  unsigned int b_j = LPC_B-listOfDualCones_B.begin();
416  boost::shared_ptr<PolyhedralCone_Rn> Ca = *(LPC_A);
417 #ifdef DEBUG
418  //std::cout << "=== Cone from A number " << LPC_A.currentIteratorNumber() << std::endl;
419  //std::cout << "=== Cone from B number " << LPC_B.currentIteratorNumber() << std::endl;
420  //std::ostringstream stream_A_;
421  //stream_A_ << "EX/ConeA";
422  //stream_A_ << LPC_A.currentIteratorNumber();
423  //stream_A_ << ".pcon";
424  //std::string coneAfile_ = stream_A_.str();
425  // Primal cones.
426  //IO_Polytope::save(coneAfile_, LPC_A.current());
427  //LPC_A.current()->checkTopologyAndGeometry();
428  //std::ostringstream stream_B_;
429  //stream_B_ << "EX/ConeB";
430  //stream_B_ << LPC_B.currentIteratorNumber();
431  //stream_B_ << ".pcon";
432  //std::string coneBfile_ = stream_B_.str();
433  //IO_Polytope::save(coneBfile_, LPC_B.current());
434  //LPC_B.current()->checkTopologyAndGeometry();
435 #endif
436  boost::shared_ptr<PolyhedralCone_Rn> intersectionCone;
437  intersectionCone.reset(new PolyhedralCone_Rn());
438 #ifdef DEBUG
439  //unsigned int i=LPC_A.currentIteratorNumber();
440  //unsigned int j=LPC_B.currentIteratorNumber();
441  //std::cout << "=== Cone from A number " << i << std::endl;
442  //std::cout << "=== Cone from B number " << j << std::endl;
443  //std::ostringstream stream_A;
444  //stream_A << "EX/ConeA";
445  //stream_A << i;
446  //stream_A << ".pcon";
447  //std::string coneAfile = stream_A.str();
448  //IO_Polytope::save(coneAfile, Ca);
449  //Ca->checkTopologyAndGeometry();
450  //std::ostringstream stream_B;
451  //stream_B << "EX/ConeB";
452  //stream_B << j;
453  //stream_B << ".pcon";
454  //std::string coneBfile = stream_B.str();
455  //IO_Polytope::save(coneBfile, Cb);
456  //Cb->checkTopologyAndGeometry();
457 #endif
458  // Fill the data structures.
459  int truncationStep = Ca->numberOfHalfSpaces();
460  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHSA(Ca->getListOfHalfSpaces());
461  for (iteHSA.begin(); iteHSA.end()!=true; iteHSA.next())
462  intersectionCone->addHalfSpace(iteHSA.current());
463  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGNA(Ca->getListOfGenerators());
464  for (iteGNA.begin(); iteGNA.end()!=true; iteGNA.next()) {
465  // Make a deep copy of each generator
466  boost::shared_ptr<Generator_Rn> gn(new Generator_Rn( *(iteGNA.current()) ));
467  intersectionCone->addGenerator(gn);
468  }
469  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHSB(Cb->getListOfHalfSpaces());
470  for (iteHSB.begin(); iteHSB.end()!=true; iteHSB.next())
471  intersectionCone->addHalfSpace(iteHSB.current());
472  // Compute the intersection between the two polyhedral cones.
473  //bool notEmpty = intersectionCone->truncate(truncationStep);
474  //std::cout << iteHS.current()->getConstant() << "\t";
475  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> >
476  lexmin_ite(intersectionCone->getListOfHalfSpaces());
477  //NoRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP;
478  StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP(intersectionCone->neigbourhoodCondition());
479  DoubleDescription<
480  boost::shared_ptr<PolyhedralCone_Rn>,
481  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> >,
482  //NoRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
483  StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
484  DD(intersectionCone, lexmin_ite, NRP, truncationStep);
485  bool notEmpty = !DD.getIsEmpty();
486 #ifdef DEBUG
487  std::cout << "@@@ nbEDG=" << intersectionCone->numberOfGenerators() << "( ";
488 #endif
489  // The second test is necessary when dealing with parallelism.
490  if (notEmpty == true && intersectionCone->numberOfGenerators() >= RnDIM) {
491  listOfDualCones_C.push_back(intersectionCone);
492  boost::shared_ptr<Generator_Rn> VX(new Generator_Rn(RnDIM));
493  VX->makeSum(*LVX_A, *LVX_B);
494 #ifdef DEBUG
495  for (unsigned int k=0; k<RnDIM; k++) {
496  std::cout << (*LVX_A)->getCoordinate(k) << "+" << (*LVX_B)->getCoordinate(k)
497  << "=" << VX->getCoordinate(k) << "\t";
498  }
499  //unsigned int u=LPC_A.currentIteratorNumber();
500  //unsigned int v=LPC_B.currentIteratorNumber();
501  //std::ostringstream stream_C;
502  //stream_C << "EX_2/ConeC";
503  //stream_C << minkowskiVertexNumber;
504  //stream_C << "--";
505  //stream_C << u;
506  //stream_C << "_";
507  //stream_C << v;
508  //stream_C << ".pcon";
509  //std::string coneCfile = stream_C.str();
510  //intersectionCone->checkTopologyAndGeometry();
511  //IO_Polytope::save(coneCfile, intersectionCone);
512 #endif
513  listOfGenerators_C.push_back(VX);
514  // Now store the bijection (a_i, b_j) <-> c_k
515  _MinkowskiDecomposition.push_back( std::make_pair(a_i, b_j) );
516  _MinkowskiDecompositionOK.push_back(true);
517  // Now write the connection between a_i and c_k and also between b_j and c_k.
518  _A2C[a_i].push_back(minkowskiVertexNumber);
519  _B2C[b_j].push_back(minkowskiVertexNumber);
520  ++minkowskiVertexNumber;
521  }
522 #ifdef DEBUG
523  std::cout << ")" << std::endl;
524 #endif
525  }}
526  }}
527 }
528 
529 boost::shared_ptr<Polytope_Rn> MinkowskiSum::compute() throw (std::domain_error) {
530  if (_firstOperand->numberOfGenerators()==0 || _firstOperand->numberOfHalfSpaces()==0 ||
531  _secondOperand->numberOfGenerators()==0 || _secondOperand->numberOfHalfSpaces()==0)
532  throw std::domain_error("MinkowskiSum::compute() needs two double-description polytopes i.e. with both vertices and half-spaces.");
533 
534  //minkowskiVertices(i,j)==1 <=> a_i + b_j is a Minkowski vertex.
535  //MinkowskiDecomposition[k] = (a_i,b_j) <=> c_k = a_i + b_j is a Minkowski vertex.
536  //_neighboursA(a_i,a_j)==1 <=> a_i R a_j
537  //_neighboursB(b_i,b_j)==1 <=> b_u R b_v
538  //_neighboursA.resize(_firstOperand->numberOfGenerators());
539  //_neighboursB.resize(_secondOperand->numberOfGenerators());
540  // _firstOperand->fillNeighbourMatrix(_neighboursA, RnDIM-1);
541  //_secondOperand->fillNeighbourMatrix(_neighboursB, RnDIM-1);
542  //std::cout << "_neighboursA: " << std::endl;
543  //{for (unsigned int i=0; i<_neighboursA.size(); ++i) {
544  //std::cout << i << ": ";
545  //std::copy(_neighboursA[i].begin(), _neighboursA[i].end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
546  //std::cout << std::endl;
547  //}}
548  //std::cout << "_neighboursB: " << std::endl;
549  //{for (unsigned int i=0; i<_neighboursB.size(); ++i) {
550  //std::cout << i << ": ";
551  //std::copy(_neighboursB[i].begin(), _neighboursB[i].end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
552  //std::cout << std::endl;
553  //}}
554 
555  // Give the relation between a vertex of A and all of its associated Minkowski vertices i.e. its polyhedrical cap
556  // _A2C[i] = [ c_u, c_v, ... ] }
557  // } => c_u = a_i + b_j
558  // _B2C[j] = [ c_u, c_w, ... ] }
559  _A2C.resize(_firstOperand->numberOfGenerators());
560  _B2C.resize(_secondOperand->numberOfGenerators());
561  std::vector< boost::shared_ptr<Generator_Rn> > listOfGenerators_A;
562  std::vector< boost::shared_ptr<Generator_Rn> > listOfGenerators_B;
563  std::vector< boost::shared_ptr<PolyhedralCone_Rn> > listOfDualCones_A;
564  std::vector< boost::shared_ptr<PolyhedralCone_Rn> > listOfDualCones_B;
565  {for (unsigned int i=0; i<_firstOperand->numberOfGenerators(); ++i) {
566  listOfGenerators_A.push_back(_firstOperand->getGenerator(i));
567  const boost::shared_ptr<PolyhedralCone_Rn>& cone_a = _firstOperand->getPrimalCone( _firstOperand->getGenerator(i) );
568  listOfDualCones_A.push_back( cone_a->computeDualPolyhedralCone() );
569  }}
570  {for (unsigned int i=0; i<_secondOperand->numberOfGenerators(); ++i) {
571  listOfGenerators_B.push_back(_secondOperand->getGenerator(i));
572  const boost::shared_ptr<PolyhedralCone_Rn>& cone_b = _secondOperand->getPrimalCone( _secondOperand->getGenerator(i) );
573  listOfDualCones_B.push_back( cone_b->computeDualPolyhedralCone() );
574  }}
575 
576  computeCommonRefinement(listOfGenerators_A, listOfGenerators_B, _NF_Vertices, listOfDualCones_A, listOfDualCones_B, _NF_Cones);
577 
578  processNormalFan2();
579 
580  //_sum->checkTopologyAndGeometry();
581  return _sum;
582 }
583 
584 void MinkowskiSum::processNormalFan2()
585 {
586  unsigned int RnDIM=Rn::getDimension();
587  double TOL=Rn::getTolerance();
588 
589  boost::numeric::ublas::matrix<double> matOfVtx(_NF_Vertices.size(), RnDIM);
590  {
591  unsigned int vtxNumber=0;
592  std::vector< boost::shared_ptr<Generator_Rn> >::iterator iteVX2;
593  {for (iteVX2 = _NF_Vertices.begin(); iteVX2!=_NF_Vertices.end(); ++iteVX2) {
594  boost::numeric::ublas::matrix_row< boost::numeric::ublas::matrix<double> > matRow(matOfVtx, vtxNumber);
595  std::copy( (*iteVX2)->begin(), (*iteVX2)->end(), matRow.begin());
596  ++vtxNumber;
597  }}
598  }
599  std::vector< std::vector< unsigned int > > VerticesPerFacet;
600  // To count the total number of edges
601  std::vector< boost::shared_ptr<HalfSpace_Rn> > setOfAllHalfSpaces;
602  std::vector< boost::shared_ptr<Generator_Rn> >::iterator iteVX = _NF_Vertices.begin();
603  std::vector< boost::shared_ptr<PolyhedralCone_Rn> >::iterator itePC = _NF_Cones.begin();
604  {for (; itePC!=_NF_Cones.end(); ++itePC, ++iteVX) {
605  // The big job begins
606  {for (unsigned int u=0; u<(*itePC)->numberOfGenerators(); ++u) {
607  // Build the corresponding half-space.
608  boost::shared_ptr<HalfSpace_Rn> HS;
609  HS.reset(new HalfSpace_Rn(RnDIM));
610  double sum=0.;
611  // The normal vector
612  for (unsigned int coord_count=0; coord_count<RnDIM; coord_count++) {
613  double this_coord = (*itePC)->getGenerator(u)->getCoordinate(coord_count);
614  //HYP->setCoefficient(coord_count, iteEDG.current()->getCoordinate(coord_count));
615  HS->setCoefficient(coord_count, this_coord);
616  sum = sum + this_coord*(*iteVX)->getCoordinate(coord_count);
617  }
618  HS->setConstant(-sum);
619  // Prepare the array of vertices IN/ON.
620  std::vector< unsigned int > vtxStates;
621  double halfSpaceNorm = std::inner_product(HS->begin(), HS->end(), HS->begin(), 0.);
622  halfSpaceNorm = sqrt(halfSpaceNorm);
623  // Now iterate on the list of vertices to classify them.
624  boost::numeric::ublas::vector<double> dist2Hyp = prod(matOfVtx, HS->vect());
625  // Possible future optimization in high dimensions
626  //boost::numeric::ublas::vector<double> dist2Hyp(NF_Vertices.size());
627  //axpy_prod(matOfVtx, HS->vect(), dist2Hyp, true);
628  boost::numeric::ublas::scalar_vector<double> scalsum(_NF_Vertices.size(),sum);
629  dist2Hyp = dist2Hyp - scalsum;
630  dist2Hyp = dist2Hyp / halfSpaceNorm;;
631  {for (unsigned int vtxNumber2=0; vtxNumber2<dist2Hyp.size(); ++vtxNumber2) {
632  if (dist2Hyp(vtxNumber2) > TOL) {
633  // Nothing to do !
634  }
635  //else if (distanceToHyperplane < -TOL) {
636  //std::cerr.precision(15);
637  //std::cerr << "distanceToHyperplane=" << distanceToHyperplane;
638  //std::cerr << std::endl;
639  //HS->dump(std::cerr);
640  //std::cerr << std::endl;
641  //(*iteVX2)->dump(std::cerr);
642  //std::cerr << std::endl;
643  //throw std::domain_error("Polytope_Rn::processNormalFan2() current vertex is out !");
644  //}
645  else {
646  vtxStates.push_back(vtxNumber2);
647  }
648  }}
649  // Now all the states are set for this current half-space.
650 
651  bool inserted=false;
652  {for (unsigned int i=0; i<VerticesPerFacet.size() && inserted==false; ++i) {
653  if (VerticesPerFacet[i].size() == vtxStates.size()) {
654  if (VerticesPerFacet[i] == vtxStates)
655  inserted = true;
656  }
657  else if (VerticesPerFacet[i].size() > vtxStates.size()) {
658  if (std::includes(VerticesPerFacet[i].begin(), VerticesPerFacet[i].end(), vtxStates.begin(), vtxStates.end()) == true) {
659  // No need to do anything as the current set of
660  // generators is included in the set number j.
661  //std::cout << "false1 (size=" << i << ")" << std::endl;
662  inserted = true;
663  }
664  }
665  else {
666  if (std::includes(vtxStates.begin(), vtxStates.end(), VerticesPerFacet[i].begin(), VerticesPerFacet[i].end()) == true) {
667  VerticesPerFacet[i] = vtxStates;
668  setOfAllHalfSpaces[i] = HS;
669  inserted = true;
670  // Do not return yet as we could need to erase another cell in the array.
671  }
672  }
673  }}
674  if (inserted == false) {
675  VerticesPerFacet.push_back(vtxStates);
676  setOfAllHalfSpaces.push_back(HS);
677  //std::cout << "inserted at the end (" << VerticesPerFacet.size() << ")" << std::endl;
678  }
679 
680  }} // {for (unsigned int u=0; u<(*itePC)->numberOfGenerators(); ++u) {
681  }} // {for (; itePC!=NF_Cones.end(); ++itePC, ++iteVX) {
682 
683  // Last job: reconstruct the polytope
684  unsigned int vtxNb=0;
685  iteVX=_NF_Vertices.begin();
686  while (iteVX != _NF_Vertices.end()) {
687  // We insert a vertex only if it has the correct number of half-spaces.
688  // This is due to the fact that some dual cones have the right number of edges
689  // but in fact are "flat" numerically speaking so in the reduction process they
690  // "lose" edges and furthermore do not have at least RnDIM edges.
691  bool foundEnoughHS=false;
692  unsigned int nbHS=0;
693  {for (unsigned int i=0; i<VerticesPerFacet.size() && foundEnoughHS==false; ++i) {
694  {for (unsigned int j=0; j<VerticesPerFacet[i].size(); ++j) {
695  if (VerticesPerFacet[i][j] == vtxNb)
696  ++nbHS;
697  if (nbHS >= RnDIM)
698  foundEnoughHS = true;
699  }}
700  }}
701  if (foundEnoughHS == true) {
702  _sum->addGenerator( (*iteVX) );
703  }
704  else {
705  // A very important step if we want to use the array _MinkowskiDecomposition: some dual cones are flat
706  // so they lose edges during the process and have foundEnoughHS to false. So we need to mark them in
707  // the corresponding array.
708  _MinkowskiDecompositionOK[vtxNb] = false;
709  }
710  ++vtxNb;
711  ++iteVX;
712  }
713  {for (unsigned int i=0; i<VerticesPerFacet.size(); ++i) {
714  {for (unsigned int j=0; j<VerticesPerFacet[i].size(); ++j) {
715  // Give the vertex its corresponding facets.
716  _NF_Vertices[ VerticesPerFacet[i][j] ]->setFacet( setOfAllHalfSpaces[i] );
717  }}
718  _sum->addHalfSpace( setOfAllHalfSpaces[i], false);
719  }}
720 }
721 
722 void MinkowskiSum::processNormalFan1()
723 {
724  unsigned int RnDIM=Rn::getDimension();
725  double TOL=Rn::getTolerance();
726  double TOL2=TOL*TOL;
727  // Give the relation between a vertex of A and all of its associated Minkowski vertices i.e. its polyhedrical cap
728  // _A2C[i] = [ c_u, c_v, ... ] }
729  // } => c_u = a_i + b_j
730  // _B2C[j] = [ c_u, c_w, ... ] }
731  // minkowskiVertices(i,j)==1 <=> a_i + b_j is a Minkowski vertex.
732  // MinkowskiDecomposition[k] = (a_i,b_j) <=> c_k = a_i + b_j is a Minkowski vertex.
733  // _neighboursA(a_i,a_j)==1 <=> a_i R a_j
734  // _neighboursB(b_i,b_j)==1 <=> b_u R b_v
735  unsigned int currentNumber=0;
736  std::vector< boost::shared_ptr<HalfSpace_Rn> > setOfAllHalfSpaces;
737  std::vector< boost::shared_ptr<Generator_Rn> >::iterator iteVX = _NF_Vertices.begin();
738  std::vector< boost::shared_ptr<PolyhedralCone_Rn> >::iterator itePC = _NF_Cones.begin();
739  // Map between an edge and its norm
740  std::map< boost::shared_ptr<Generator_Rn>, double > allEdgesNorms;
741  {for (; itePC!=_NF_Cones.end(); ++itePC, ++iteVX) {
742  // Build a half-space from the dual polyhedral cone generator.
743  // So iterate on the list of its edges.
744 
745  // First of all get the ancestors of the current Minkowski vertex
746  unsigned int a_i = _MinkowskiDecomposition[ currentNumber ].first;
747  unsigned int b_j = _MinkowskiDecomposition[ currentNumber ].second;
748  // Now get all of their facet neighbours
749  // Now inspect all the possible combinations to know whether they provide a Minkowski vertex.
750  std::vector< unsigned int > c_kMinkowskiVertices;
751  std::vector< unsigned int >::const_iterator iteNgb;
752  // Now get all a_i neighbours and find their contributions in terms of Minkowski vertices
753  {for (unsigned int j=0; j<_neighboursA[a_i].size(); ++j) {
754  unsigned int a_i_ngb = _neighboursA[a_i][j];
755  {for (unsigned int k=0; k<_A2C[a_i_ngb].size(); ++k) {
756  // Make sure we will not process the current normal fan with itself.
757  if (_A2C[a_i_ngb][k] != currentNumber)
758  c_kMinkowskiVertices.push_back( _A2C[a_i_ngb][k] );
759  }}
760  }}
761  {for (unsigned int i=0; i<_neighboursB[b_j].size(); ++i) {
762  unsigned int b_j_ngb=_neighboursB[b_j][i];
763  {for (unsigned int k=0; k<_B2C[b_j_ngb].size(); ++k) {
764  if (_B2C[b_j_ngb][k] != currentNumber)
765  c_kMinkowskiVertices.push_back( _B2C[b_j_ngb][k] );
766  }}
767  }}
768  //std::copy(c_kMinkowskiVertices.begin(), c_kMinkowskiVertices.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
769 
770  // The big job begins
771  {for (unsigned int u=0; u<(*itePC)->numberOfGenerators(); ++u) {
772  // Build the corresponding half-space.
773  boost::shared_ptr<HalfSpace_Rn> HS;
774  HS.reset(new HalfSpace_Rn(RnDIM));
775  double sum=0.;
776  // The normal vector
777  for (unsigned int coord_count=0; coord_count<RnDIM; coord_count++) {
778  double this_coord = (*itePC)->getGenerator(u)->getCoordinate(coord_count);
779  //HYP->setCoefficient(coord_count, iteEDG.current()->getCoordinate(coord_count));
780  HS->setCoefficient(coord_count, this_coord);
781  sum = sum + this_coord*(*iteVX)->getCoordinate(coord_count);
782  }
783  HS->setConstant(-sum);
784  const boost::shared_ptr<Generator_Rn>& gn1=(*itePC)->getGenerator(u);
785  double gn1_sum2 = std::inner_product(gn1->begin(), gn1->end(), gn1->begin(), 0.);
786  {for (iteNgb=c_kMinkowskiVertices.begin(); iteNgb!=c_kMinkowskiVertices.end(); ++iteNgb) {
787  // Check all generators of the neighbour cone against the current one.
788  bool check=true;
789  {for (unsigned int v=0; v<_NF_Cones[*iteNgb]->numberOfGenerators() && check==true; ++v) {
790  const boost::shared_ptr<Generator_Rn>& gn2=_NF_Cones[*iteNgb]->getGenerator(v);
791  double scal_prod = std::inner_product(gn1->begin(), gn1->end(), gn2->begin(), 0.);
792  if (scal_prod > 0.) {
793  double coef1 = scal_prod / gn1_sum2;
794  // Check whether the current edge is in the "tube" of the second one, if not switch them.
795  if (gn1->getNormalDistance(gn2, coef1, RnDIM) < TOL2) {
796  // We found a copy of the current edge so remove it from the explored polyhedrical cone.
797  _NF_Cones[*iteNgb]->removeGenerator(v);
798  _NF_Vertices[*iteNgb]->setFacet(HS);
799  check = false;
800  }
801  else if (scal_prod > 0.) {
802  double gn2_sum2 = 0.;
803  // First try to get the norm if it was previously computed
804  std::map< boost::shared_ptr<Generator_Rn>, double >::iterator iteNorm = allEdgesNorms.find( gn2 );
805  if (iteNorm == allEdgesNorms.end())
806  gn2_sum2 = std::inner_product(gn2->begin(), gn2->end(), gn2->begin(), 0.);
807  else
808  gn2_sum2 = iteNorm->second;
809  double coef2 = scal_prod / gn2_sum2;
810  if (gn2->getNormalDistance(gn1, coef2, RnDIM) < TOL2) {
811  _NF_Cones[*iteNgb]->removeGenerator(v);
812  _NF_Vertices[*iteNgb]->setFacet(HS);
813  check = false;
814  if (iteNorm != allEdgesNorms.end())
815  allEdgesNorms.erase(iteNorm);
816  }
817  }
818  }
819  }}
820  }}
821  // Insert the half-space in the polytope list, if not already in,
822  // and insert it into the vertex facet list.
823  (*iteVX)->setFacet(HS);
824  setOfAllHalfSpaces.push_back( HS );
825  // Now remove the current generator as it has become useless.
826  (*itePC)->removeGenerator(u);
827  }} // {for (unsigned int u=0; u<(*itePC)->numberOfGenerators() ; ++u) {
828 
829  // Insert the vertex into the polytope vertex list.
830  _sum->addGenerator( (*iteVX) );
831  ++currentNumber;
832  }} // {for (; itePC!=NF_Cones.end(); ++itePC, ++iteVX) {
833 
834  // Check everything's fine
835  {for (itePC=_NF_Cones.begin(); itePC!=_NF_Cones.end(); ++itePC) {
836  if ((*itePC)->numberOfGenerators() != 0)
837  throw std::domain_error("MinkowskiSum::processNormalFan1() dual cones reduction not operational !");
838  }}
839 
840  // Transfer all half-spaces into the main list.
841  std::vector< boost::shared_ptr<HalfSpace_Rn> >::const_iterator iteHS=setOfAllHalfSpaces.begin();
842  {for (iteHS=setOfAllHalfSpaces.begin(); iteHS!=setOfAllHalfSpaces.end(); ++iteHS) {
843  // Thanks to the previous work we do not need to check.
844  _sum->addHalfSpace( *iteHS, false);
845  }}
846 
847 }
848 
849 void MinkowskiSum::processNormalFan0()
850 {
851  unsigned int RnDIM=Rn::getDimension();
852 
853  std::vector< boost::shared_ptr<Generator_Rn> >::iterator iteVX=_NF_Vertices.begin();
854  std::vector< boost::shared_ptr<PolyhedralCone_Rn> >::iterator itePC=_NF_Cones.begin();
855  for (; itePC!=_NF_Cones.end(); ++itePC, ++iteVX) {
856  // Build a half-space from the dual polyhedral cone generator.
857  // So iterate on the list of its edges.
858 
859  //constIteratorOfListOfGenerators iteEDG(itePC.current()->getListOfGenerators());
860  //for (iteEDG.begin(); iteEDG.end()!=true; iteEDG.next()) {
861  for (unsigned int u=0; u<(*itePC)->numberOfGenerators(); u++) {
862  boost::shared_ptr<HalfSpace_Rn> HS;
863  HS.reset(new HalfSpace_Rn(RnDIM));
864  double sum=0.;
865  // The normal vector
866  for (unsigned int coord_count=0; coord_count<RnDIM; coord_count++) {
867  double this_coord = (*itePC)->getGenerator(u)->getCoordinate(coord_count);
868  //HYP->setCoefficient(coord_count, iteEDG.current()->getCoordinate(coord_count));
869  HS->setCoefficient(coord_count, this_coord);
870  sum = sum + this_coord*(*iteVX)->getCoordinate(coord_count);
871  }
872  HS->setConstant(-sum);
873  // Insert the half-space in the polytope list, if not already in,
874  // and insert it into the vertex facet list.
875  (*iteVX)->setFacet(_sum->addHalfSpace(HS, true));
876  }
877  // Insert the vertex into the polytope vertex list.
878  _sum->addGenerator( (*iteVX) );
879  }
880 
881 }
882 
883 //========================================================================================
884 
886 void PseudoSumWithoutCaps::computeCapHalfSpaces(
887  const std::set< unsigned int >& firstOperandCaps,
888  const std::set< unsigned int >& secondOperandCaps,
889  std::set< unsigned int >& sumCaps) throw (std::domain_error)
890 {
891 #ifdef DEBUG
892  std::cout << "PseudoSumWithoutCaps::computeCapHalfSpaces()" << endl;
893 #endif
894  if (_sum->numberOfGenerators()==0)
895  throw std::domain_error("PseudoSumWithoutCaps::computeCapHalfSpaces() We need to compute the Minkowski sum first.");
896 
897  // For each facet of the sum, store its vertices.
898  unsigned int facetNumber=0;
899  std::vector< std::vector<unsigned int> > listOfVerticesPerFacet(_sum->numberOfHalfSpaces());
900  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS(_sum->getListOfHalfSpaces());
901  {for (iteHS.begin(); iteHS.end()!=true; iteHS.next()) {
902 #ifdef DEBUG
903  std::cout << "facetNumber=" << facetNumber << ":";
904 #endif
905  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(_sum->getListOfGenerators());
906  {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
907  for (unsigned int i=0; i<iteGN.current()->numberOfFacets(); ++i) {
908  if (iteGN.current()->getFacet(i) == iteHS.current()) {
909  listOfVerticesPerFacet[facetNumber].push_back(iteGN.currentIteratorNumber());
910 #ifdef DEBUG
911  std::cout << " " << iteGN.currentIteratorNumber();
912 #endif
913  }
914  }
915  }}
916  ++facetNumber;
917 #ifdef DEBUG
918  std::cout << endl;
919 #endif
920  }}
921 
922  // Reorder the data structures in the father class.
923  std::vector< std::pair< unsigned int, unsigned int > > thisMinkDecomposition;
924  std::vector< boost::shared_ptr<PolyhedralCone_Rn> > theseNF_Cones;
925  std::vector< boost::shared_ptr<Generator_Rn> > theseNF_Vertices;
926  {for (unsigned int i=0; i<_MinkowskiDecompositionOK.size(); ++i) {
927  if (_MinkowskiDecompositionOK[i] == true) {
928  thisMinkDecomposition.push_back( _MinkowskiDecomposition[i] );
929  theseNF_Vertices.push_back( _NF_Vertices[i] );
930  theseNF_Cones.push_back( _NF_Cones[i] );
931  }
932  }}
933  _MinkowskiDecomposition = thisMinkDecomposition;
934  _NF_Vertices = theseNF_Vertices;
935  _NF_Cones = theseNF_Cones;
936 
937  // Now facet by facet, decompose the Minkowski vertices into vertices of A and B.
938  unsigned int sumFacetNb = 0;
939  std::vector< std::vector<unsigned int> >::const_iterator itVtx = listOfVerticesPerFacet.begin();
940  {for (; itVtx!=listOfVerticesPerFacet.end(); ++itVtx) {
941  std::set< unsigned int > listOfVtxA, listOfVtxB;
942  std::vector<unsigned int>::const_iterator MinkVtx = itVtx->begin();
943 #ifdef DEBUG
944  std::cout << endl << "Compute F_A and F_B for facet " << sumFacetNb << endl;
945 #endif
946  {for (; MinkVtx!=itVtx->end(); ++MinkVtx) {
947  unsigned int MinkNumber = *MinkVtx;
948  listOfVtxA.insert( _MinkowskiDecomposition[MinkNumber].first );
949  listOfVtxB.insert( _MinkowskiDecomposition[MinkNumber].second);
950  }}
951 
952  std::set<unsigned int> listOfFacetsOfVtxA;
953  std::set<unsigned int> tmpInterResForA, tmp2InterResForA;
954  std::set<unsigned int>::const_iterator itVtxNbOperand;
955  itVtxNbOperand = listOfVtxA.begin();
956  for (unsigned int ngb_count=0; ngb_count<_firstOperand->getGenerator(*itVtxNbOperand)->numberOfFacets(); ngb_count++) {
957  boost::shared_ptr<HalfSpace_Rn> Fi = _firstOperand->getGenerator(*itVtxNbOperand)->getFacet(ngb_count);
958  listOfFacetsOfVtxA.insert(_firstOperand->getHalfSpaceNumber(Fi));
959  }
960  //_firstOperand->getGenerator(*itVtxNbOperand)->exportFacets(listOfFacetsOfVtxA);
961 #ifdef DEBUG
962  std::cout << "F_A:";
963  std::cout << " V = {";
964  std::cout << " " << *itVtxNbOperand;
965 #endif
966  tmpInterResForA = listOfFacetsOfVtxA;
967  itVtxNbOperand++;
968  {for (; itVtxNbOperand!=listOfVtxA.end(); ++itVtxNbOperand) {
969 #ifdef DEBUG
970  std::cout << " " << *itVtxNbOperand;
971 #endif
972  std::set<unsigned int> curListOfFacetsOfVtxA;
973  for (unsigned int ngb_count=0; ngb_count<_firstOperand->getGenerator(*itVtxNbOperand)->numberOfFacets(); ngb_count++) {
974  boost::shared_ptr<HalfSpace_Rn> Fi = _firstOperand->getGenerator(*itVtxNbOperand)->getFacet(ngb_count);
975  curListOfFacetsOfVtxA.insert(_firstOperand->getHalfSpaceNumber(Fi));
976  }
977  tmp2InterResForA = tmpInterResForA;
978  tmpInterResForA.clear();
979  std::set_intersection(
980  tmp2InterResForA.begin(),
981  tmp2InterResForA.end(),
982  curListOfFacetsOfVtxA.begin(),
983  curListOfFacetsOfVtxA.end(),
984  std::inserter(tmpInterResForA, tmpInterResForA.end()));
985  if (listOfFacetsOfVtxA.empty() == true)
986  throw std::domain_error("PseudoSumWithoutCaps::computeCapHalfSpaces() the i-Face of polytope A has 0 half-space.");
987  }}
988 #ifdef DEBUG
989  std::cout << " }";
990 #endif
991  // At this step the half-spaces support of F_A are computed.
992  std::set<unsigned int> InterResForA;
993  std::set_intersection(
994  tmpInterResForA.begin(),
995  tmpInterResForA.end(),
996  firstOperandCaps.begin(),
997  firstOperandCaps.end(),
998  std::inserter(InterResForA, InterResForA.end()));
999 #ifdef DEBUG
1000  std::cout << " ; H = { ";
1001  std::copy(InterResForA.begin(), InterResForA.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
1002  //std::cout << listOfFacetsOfVtxA.size();
1003  std::cout << "}";
1004  std::cout << endl;
1005 #endif
1006 
1007  if (InterResForA.empty() != true)
1008  sumCaps.insert( sumFacetNb );
1009  else {
1011  // 2nd part //
1013  // Get the facets of all the vertices of the j-face in B and intersect them to extract the support.
1014  std::set<unsigned int> listOfFacetsOfVtxB;
1015  std::set<unsigned int> tmpInterResForB, tmp2InterResForB;
1016  itVtxNbOperand = listOfVtxB.begin();
1017  for (unsigned int ngb_count=0; ngb_count<_secondOperand->getGenerator(*itVtxNbOperand)->numberOfFacets(); ngb_count++) {
1018  boost::shared_ptr<HalfSpace_Rn> Fi = _secondOperand->getGenerator(*itVtxNbOperand)->getFacet(ngb_count);
1019  listOfFacetsOfVtxB.insert(_secondOperand->getHalfSpaceNumber(Fi));
1020  }
1021 #ifdef DEBUG
1022  std::cout << "F_B:";
1023  std::cout << " V = {";
1024  std::cout << " " << *itVtxNbOperand;
1025 #endif
1026  tmpInterResForB = listOfFacetsOfVtxB;
1027  itVtxNbOperand++;
1028  {for (; itVtxNbOperand!=listOfVtxB.end(); ++itVtxNbOperand) {
1029 #ifdef DEBUG
1030  std::cout << " " << *itVtxNbOperand;
1031 #endif
1032  std::set<unsigned int> curListOfFacetsOfVtxB;
1033  for (unsigned int ngb_count=0; ngb_count<_secondOperand->getGenerator(*itVtxNbOperand)->numberOfFacets(); ngb_count++) {
1034  boost::shared_ptr<HalfSpace_Rn> Fi = _secondOperand->getGenerator(*itVtxNbOperand)->getFacet(ngb_count);
1035  curListOfFacetsOfVtxB.insert(_secondOperand->getHalfSpaceNumber(Fi));
1036  }
1037  tmp2InterResForB = tmpInterResForB;
1038  tmpInterResForB.clear();
1039  std::set_intersection(
1040  tmp2InterResForB.begin(),
1041  tmp2InterResForB.end(),
1042  curListOfFacetsOfVtxB.begin(),
1043  curListOfFacetsOfVtxB.end(),
1044  std::inserter(tmpInterResForB, tmpInterResForB.end()));
1045  if (listOfFacetsOfVtxB.empty() == true)
1046  throw std::domain_error("PseudoSumWithoutCaps::computeCapHalfSpaces() the j-Face of polytope B has 0 half-space.");
1047  }}
1048 #ifdef DEBUG
1049  std::cout << " }";
1050 #endif
1051  // Now we have all the geometric supports of F_B, see whether there is a cap half-space inside.
1052  std::set<unsigned int> InterResForB;
1053  std::set_intersection(
1054  tmpInterResForB.begin(),
1055  tmpInterResForB.end(),
1056  secondOperandCaps.begin(),
1057  secondOperandCaps.end(),
1058  std::inserter(InterResForB, InterResForB.end()));
1059 #ifdef DEBUG
1060  std::cout << " ; H = { ";
1061  std::copy(InterResForB.begin(), InterResForB.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
1062  //std::cout << listOfFacetsOfVtxB.size();
1063  std::cout << "}";
1064  std::cout << endl;
1065 #endif
1066 
1067  // If at least one cap half-space is in listOfFacetsOfVtxA or listOfFacetsOfVtxB then the current half-space of the sum is capped too.
1068  if (InterResForB.empty() != true)
1069  sumCaps.insert( sumFacetNb );
1070  }
1071 
1072  sumFacetNb++;
1073  }}
1074 
1075 }
1076 
1077 // Remove the cap half-spaces stored in sets and then truncate again.
1078 boost::shared_ptr<Polytope_Rn> PseudoSumWithoutCaps::rebuildSum(
1079  const std::set<unsigned int>& firstOperandCaps,
1080  const std::set<unsigned int>& secondOperandCaps,
1081  std::set<unsigned int>& newCaps,
1082  double bb_size)
1083 {
1084  if (_sum->numberOfHalfSpaces() == 0 || _sum->numberOfGeneratorsPerFacet() == 0)
1085  throw std::domain_error("PseudoSumWithoutCaps::rebuildSum() _sum is not computed.");
1086 
1087  std::set< unsigned int > sumCaps;
1088  computeCapHalfSpaces(firstOperandCaps, secondOperandCaps, sumCaps);
1089 #ifdef DEBUG
1090  std::cout << "firstOperandCaps=" << firstOperandCaps.size() << std::endl;
1091  std::cout << "secondOperandCaps=" << secondOperandCaps.size() << std::endl;
1092  std::cout << "sumCaps=" << sumCaps.size() << std::endl;
1093 #endif
1094 
1095  // Now remove all cap half-spaces from the current polytope and create another one without them.
1096  boost::shared_ptr<Polytope_Rn> newSum;
1097  newSum.reset(new Polytope_Rn());
1098  newSum->createBoundingBox(bb_size);
1099  std::vector< boost::shared_ptr<HalfSpace_Rn> > tryCapsForNewSum;
1100  // Create the new cap half-spaces.
1101  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS1(newSum->getListOfHalfSpaces());
1102  for (iteHS1.begin(); iteHS1.end()!=true; iteHS1.next()) {
1103  tryCapsForNewSum.push_back(iteHS1.current());
1104  }
1105  // Copy the non cap half-spaces.
1106  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS2(_sum->getListOfHalfSpaces());
1107  for (iteHS2.begin(); iteHS2.end()!=true; iteHS2.next()) {
1108  if (sumCaps.find(iteHS2.currentIteratorNumber()) == sumCaps.end())
1109  newSum->addHalfSpace(iteHS2.current());
1110  }
1111 
1112  // Now the truncation.
1113  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite(newSum->getListOfHalfSpaces());
1114  StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP(newSum->neigbourhoodCondition());
1115  // The number of facets for a cube.
1116  unsigned int truncationStep = 2*Rn::getDimension();
1117  DoubleDescription<
1118  boost::shared_ptr<PolyhedralCone_Rn>,
1119  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> >,
1120  StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
1121  DD(newSum, lexmin_ite, NRP, truncationStep);
1122 
1123  // Build the new list of cap half-spaces for _sum
1124  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS3(newSum->getListOfHalfSpaces());
1125  for (iteHS3.begin(); iteHS3.end()!=true; iteHS3.next()) {
1126  std::vector< boost::shared_ptr<HalfSpace_Rn> >::const_iterator iteCaps = tryCapsForNewSum.begin();
1127  for (; iteCaps!=tryCapsForNewSum.end(); ++iteCaps) {
1128  if (iteHS3.current() == *iteCaps) {
1129  newCaps.insert(iteHS3.currentIteratorNumber());
1130  break;
1131  }
1132  }
1133  }
1134  _sum = newSum;
1135 
1136  return newSum;
1137 }
1138 
1139 //========================================================================================
1140 
1141 PseudoIntersectionWithoutCaps::PseudoIntersectionWithoutCaps(
1142  const boost::shared_ptr<Polytope_Rn>& A,
1143  const boost::shared_ptr<Polytope_Rn>& B,
1144  boost::shared_ptr<Polytope_Rn>& C,
1145  const std::set< unsigned int >& firstOperandCaps,
1146  const std::set< unsigned int >& secondOperandCaps,
1147  std::set< unsigned int >& newCaps,
1148  double bb_size) {
1149 
1150  // Prepare the intersection polytope.
1151  C.reset(new Polytope_Rn());
1152  C->createBoundingBox(bb_size);
1153  std::vector< boost::shared_ptr<HalfSpace_Rn> > RnCaps;
1154  // Get the cube cap half-spaces
1155  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS0(C->getListOfHalfSpaces());
1156  for (iteHS0.begin(); iteHS0.end()!=true; iteHS0.next()) {
1157  RnCaps.push_back(iteHS0.current());
1158  }
1159  // Copy the non cap half-spaces of A.
1160  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS1(A->getListOfHalfSpaces());
1161  for (iteHS1.begin(); iteHS1.end()!=true; iteHS1.next()) {
1162  if (firstOperandCaps.find(iteHS1.currentIteratorNumber()) == firstOperandCaps.end())
1163  C->addHalfSpace(iteHS1.current());
1164  }
1165  // Copy the non cap half-spaces of B.
1166  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS2(B->getListOfHalfSpaces());
1167  for (iteHS2.begin(); iteHS2.end()!=true; iteHS2.next()) {
1168  if (secondOperandCaps.find(iteHS2.currentIteratorNumber()) == secondOperandCaps.end())
1169  C->addHalfSpace(iteHS2.current());
1170  }
1171 
1172  // Now the truncation.
1173  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite(C->getListOfHalfSpaces());
1174  StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP(C->neigbourhoodCondition());
1175  // The number of facets for a cube.
1176  unsigned int truncationStep = 2*Rn::getDimension();
1177  DoubleDescription<
1178  boost::shared_ptr<PolyhedralCone_Rn>,
1179  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> >,
1180  StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
1181  DD(C, lexmin_ite, NRP, truncationStep);
1182 
1183  // Build the new list of cap half-spaces for C
1184  newCaps.clear();
1185  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS3(C->getListOfHalfSpaces());
1186  for (iteHS3.begin(); iteHS3.end()!=true; iteHS3.next()) {
1187  std::vector< boost::shared_ptr<HalfSpace_Rn> >::const_iterator iteCaps = RnCaps.begin();
1188  for (; iteCaps!=RnCaps.end(); ++iteCaps) {
1189  if (iteHS3.current() == *iteCaps) {
1190  newCaps.insert(iteHS3.currentIteratorNumber());
1191  break;
1192  }
1193  }
1194  }
1195 }
1196 
1197 //========================================================================================
1198 
1199 int TopGeomTools::Translate(boost::shared_ptr<Polytope_Rn>& pol, const boost::numeric::ublas::vector<double>& v2t) {
1200  if (pol->numberOfGenerators() != 0) {
1201  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(pol->getListOfGenerators());
1202  {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
1203  const boost::shared_ptr<Generator_Rn>& v2_A = iteGN.current();
1204  boost::numeric::ublas::vector<double> coord = v2_A->vect() + v2t;
1205  v2_A->setCoordinates(coord);
1206  }}
1207  }
1208  if (pol->numberOfHalfSpaces() != 0) {
1209  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS(pol->getListOfHalfSpaces());
1210  {for (iteHS.begin(); iteHS.end()!=true; iteHS.next()) {
1211  const boost::shared_ptr<HalfSpace_Rn>& hs = iteHS.current();
1212  double prod = std::inner_product(hs->begin(), hs->end(), v2t.begin(), 0.);
1213  hs->setConstant( hs->getConstant()-prod );
1214  }}
1215  }
1216  return 0;
1217 }
1218 
1219 int TopGeomTools::GravityCenter(boost::shared_ptr<Polytope_Rn>& pol, boost::numeric::ublas::vector<double>& gravity_center) {
1220  if (pol->numberOfGenerators() != 0) {
1221  for (unsigned int i=0; i<gravity_center.size (); ++i)
1222  gravity_center(i) = 0.;
1223  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(pol->getListOfGenerators());
1224  {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
1225  gravity_center += iteGN.current()->vect();
1226  }}
1227  gravity_center /= pol->numberOfGenerators();
1228  }
1229  else
1230  throw std::domain_error("DoubleDescriptionFromGenerators::Compute() the polytope already does not have generators.");
1231  return 0;
1232 }
1233 
1234 int TopGeomTools::scalingFactor(boost::shared_ptr<Polytope_Rn>& pol, double factor) {
1235  unsigned int dim = Rn::getDimension();
1236  if (fabs(factor) < Rn::getTolerance())
1237  throw std::domain_error("TopGeomTools::scalingFactor() scaling factor too small.");
1238  if (pol->numberOfGenerators() != 0) {
1239  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(pol->getListOfGenerators());
1240  {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
1241  for (unsigned int i=0; i<dim; ++i)
1242  iteGN.current()->setCoordinate(i, iteGN.current()->getCoordinate(i)*factor);
1243  }}
1244  }
1245  if (pol->numberOfHalfSpaces() != 0) {
1246  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS(pol->getListOfHalfSpaces());
1247  for (iteHS.begin(); iteHS.end()!=true; iteHS.next())
1248  iteHS.current()->setConstant( iteHS.current()->getConstant()*factor);
1249  }
1250 
1251  return 0;
1252 }
1253 
1254 int TopGeomTools::PolarPolytope(
1255  const boost::shared_ptr<Polytope_Rn>& original_pol,
1256  boost::shared_ptr<Polytope_Rn>& polar_pol,
1257  bool forceComputation, double bb_size) throw (invalid_argument) {
1258  polar_pol.reset(new Polytope_Rn());
1259  unsigned int dim = original_pol->dimension();
1260  double TOL = Rn::getTolerance();
1261  if (original_pol->numberOfGenerators() != 0) {
1262  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(original_pol->getListOfGenerators());
1263  {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
1264  boost::shared_ptr<HalfSpace_Rn> HS;
1265  HS.reset(new HalfSpace_Rn(dim));
1266  HS->setConstant(1.);
1267  boost::numeric::ublas::vector< double >::const_iterator iteCoord = iteGN.current()->begin();
1268  unsigned int coord_count = 0;
1269  {for ( ; iteCoord != iteGN.current()->end(); ++iteCoord ) {
1270  //std::cout << *iteCoord << " ";
1271  HS->setCoefficient(coord_count, *iteCoord);
1272  ++coord_count;
1273  }}
1274  //std::cout << std::endl;
1275  polar_pol->addHalfSpace(HS);
1276  }}
1277  if (forceComputation == true) {
1278  try {
1279  polar_pol->createBoundingBox(bb_size);
1280  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite(polar_pol->getListOfHalfSpaces());
1281  StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP(polar_pol->neigbourhoodCondition());
1282  // The number of facets for a cube.
1283  unsigned int truncationStep = 2*dim;
1284  DoubleDescription<
1285  boost::shared_ptr<PolyhedralCone_Rn>,
1286  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> >,
1287  StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
1288  DD(polar_pol, lexmin_ite, NRP, truncationStep);
1289  }
1290  catch(invalid_argument& except) {
1291  cerr << "Invalid argument exception in TopGeomTools::PolarPolytope() " << except.what() << endl;
1292  return -1;
1293  }
1294  catch(out_of_range& except) {
1295  cerr << "Out of range exception in TopGeomTools::PolarPolytope() " << except.what() << endl;
1296  return -1;
1297  }
1298  catch(ios_base::failure& except) {
1299  cerr << "In/out exception in TopGeomTools::PolarPolytope() " << except.what() << endl;
1300  return -1;
1301  }
1302  catch(logic_error& except) {
1303  cerr << "Logic error exception in TopGeomTools::PolarPolytope() " << except.what() << endl;
1304  return -1;
1305  }
1306  catch(...) {
1307  cerr << "Unexpected exception caught in TopGeomTools::PolarPolytope() !" << endl;
1308  return -1;
1309  }
1310  }
1311  // Whether there were or not half-spaces, we leave now.
1312  return 0;
1313  }
1314  // At this step we know that generators were not provided.
1315  if (original_pol->numberOfHalfSpaces() != 0) {
1316  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS(original_pol->getListOfHalfSpaces());
1317  {for (iteHS.begin(); iteHS.end()!=true; iteHS.next()) {
1318  if (iteHS.current()->getConstant() < TOL)
1319  throw std::invalid_argument("TopGeomTools::PolarPolytope() the input polytope should contain the origin.");
1320  boost::numeric::ublas::vector<double> copyOfHSCoef(iteHS.current()->vect());
1321  // Scale the hyperplane such as the constant is equal to 1.
1322  copyOfHSCoef = copyOfHSCoef / iteHS.current()->getConstant();
1323  boost::shared_ptr<Generator_Rn> GN;
1324  GN.reset(new Generator_Rn(dim));
1325  GN->setCoordinates(copyOfHSCoef);
1326  //if (original_pol->numberOfGenerators() != 0) {
1327  //constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > polar_iteHS(polar_pol->getListOfHalfSpaces());
1328  //{for (polar_iteHS.begin(); polar_iteHS.end()!=true; polar_iteHS.next()) {
1329  //double halfSpaceNorm = std::inner_product(polar_iteHS.current()->begin(), polar_iteHS.current()->end(), polar_iteHS.current()->begin(), 0.);
1330  //halfSpaceNorm = sqrt(halfSpaceNorm);
1331  //double scalarProduct = std::inner_product(GN->begin(), GN->end(), polar_iteHS.current()->begin(), 0.);
1332  //double distanceToHyperplane = (scalarProduct+polar_iteHS.current()->getConstant()) / halfSpaceNorm;
1333  //if (distanceToHyperplane<TOL && distanceToHyperplane>-TOL)
1334  //GN->setFacet(polar_iteHS.current());
1335  //}}
1336  //}
1337  polar_pol->addGenerator(GN);
1338  }}
1339  }
1340 
1341  return 0;
1342 }
1343 
1344 int TopGeomTools::projectPolytopeOnCanonicalHyperplanes(
1345  const std::set< unsigned int >& listOfHyperplanes,
1346  const boost::shared_ptr<Polytope_Rn>& original_pol,
1347  boost::shared_ptr<Polytope_Rn>& proj_pol) throw (invalid_argument)
1348 {
1349  // Here we really need to work with Rn::getDimension()
1350  unsigned int currentDimension = Rn::getDimension();
1351  if (listOfHyperplanes.empty()==true || listOfHyperplanes.size()>Rn::getDimension())
1352  throw std::invalid_argument("TopGeomTools::projectPolytopeOnCanonicalHyperplanes() wrong list of hyperplanes to project");
1353  unsigned int projectedDimension = listOfHyperplanes.size();
1354 
1355  Rn::setDimension(projectedDimension);
1356  proj_pol.reset(new Polytope_Rn());
1357 
1358  std::vector< boost::shared_ptr<Generator_Rn> > newArrayOfGN;
1359  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(original_pol->getListOfGenerators());
1360  {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
1361  unsigned int j = 0;
1362  boost::shared_ptr<Generator_Rn> VX;
1363  VX.reset(new Generator_Rn(projectedDimension));
1364  for (unsigned int i=1; i<=currentDimension; ++i) {
1365  if (listOfHyperplanes.find(i) != listOfHyperplanes.end()) {
1366  // The direction i is significant so store the corresponding coordinate.
1367  VX->setCoordinate(j, iteGN.current()->getCoordinate(i-1));
1368  ++j;
1369  }
1370  }
1371  bool foundEqual = false;
1372  std::vector< boost::shared_ptr<Generator_Rn> >::const_iterator iteNew = newArrayOfGN.begin();
1373  for ( ; iteNew != newArrayOfGN.end() && foundEqual == false; ++iteNew) {
1374  if (VX->distanceFrom(**iteNew) < Rn::getTolerance())
1375  foundEqual = true;
1376  }
1377  if (foundEqual == false)
1378  newArrayOfGN.push_back( VX );
1379  }}
1380 
1381  std::vector< boost::shared_ptr<Generator_Rn> >::const_iterator iteNew2 = newArrayOfGN.begin();
1382  for ( ; iteNew2 != newArrayOfGN.end(); ++iteNew2)
1383  proj_pol->addGenerator(*iteNew2);
1385  //std::cout << "DUMP" << std::endl;
1386  //proj_pol->dump(std::cout);
1387 
1388  DoubleDescriptionFromGenerators::Compute(proj_pol, 10000.);
1389  proj_pol->checkTopologyAndGeometry();
1390 
1391  Rn::setDimension(currentDimension);
1392  return 0;
1393 }
1394 
1395 int TopGeomTools::extrudeInCanonicalDirections(
1396  const std::set< unsigned int >& originalSpaceDirections,
1397  unsigned int totalSpaceDimension,
1398  const boost::shared_ptr<Polytope_Rn>& original_polytope,
1399  boost::shared_ptr<PolyhedralCone_Rn>& extruded_polyhedron) throw (invalid_argument) {
1400 
1401  unsigned int originalDimension = Rn::getDimension();
1402  if (originalSpaceDirections.empty()==true || originalSpaceDirections.size()>Rn::getDimension())
1403  throw std::invalid_argument("TopGeomTools::extrudeInCanonicalDirections() wrong definition of the original space");
1404 
1405  // Careful as we are going to mix the dimensions of the different spaces.
1406  Rn::setDimension(totalSpaceDimension);
1407  //extruded_polyhedron.reset(new PolyhedralCone_Rn());
1408 
1409  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS(original_polytope->getListOfHalfSpaces());
1410  {for (iteHS.begin(); iteHS.end()!=true; iteHS.next()) {
1411  boost::shared_ptr<HalfSpace_Rn> newHS;
1412  newHS.reset(new HalfSpace_Rn(totalSpaceDimension));
1413  newHS->setConstant(iteHS.current()->getConstant());
1414  unsigned int shift = 0;
1415  for (unsigned int i=1; i<=totalSpaceDimension; ++i) {
1416  // Test whether i belong to the previously referenced directions.
1417  if (originalSpaceDirections.find(i) != originalSpaceDirections.end()) {
1418  newHS->setCoefficient(i-1, iteHS.current()->getCoefficient(i-shift-1));
1419  }
1420  else {
1421  newHS->setCoefficient(i-1, 0.);
1422  ++shift;
1423  }
1424  }
1425  extruded_polyhedron->addHalfSpace(newHS);
1426 
1427  }}
1428 
1429  // Back to the previous space
1430  Rn::setDimension(originalDimension);
1431  return 0;
1432 }
1433 
1434 //========================================================================================
1435 
1436 int DoubleDescriptionFromGenerators::Compute(boost::shared_ptr<Polytope_Rn>& pol, double bb_size)
1437  throw (invalid_argument, out_of_range, ios_base::failure, logic_error) {
1438  unsigned int dim = pol->dimension();
1439  //double TOL = Rn::getTolerance();
1440  if (pol->numberOfHalfSpaces() != 0)
1441  throw std::domain_error("DoubleDescriptionFromGenerators::Compute() the polytope already has half-spaces.");
1442  if (pol->numberOfGenerators() == 0)
1443  throw std::domain_error("DoubleDescriptionFromGenerators::Compute() the polytope does not have generators.");
1444  // Step 1: computes the gravity center
1445  // Step 2: translate the polytope to have gravity centered on the origin
1446  // Step 3: turn the V-description polytope into its H-description polar
1447  // Step 4: get the V-description from the H-description polar
1448  // Step 5: get the polar of the polar
1449  // Step 6: translate it back
1450 
1451  // Step 1: compute the gravity center of the cloud of points.
1452  boost::numeric::ublas::vector<double> gravity_center(dim);
1453  TopGeomTools::GravityCenter(pol, gravity_center);
1454 
1455  // Step 2: translate the cloud to have in centered on the origin.
1456  TopGeomTools::Translate(pol, -gravity_center);
1457 
1458  // Step 3: compute the polar H-description and truncate it.
1459  boost::shared_ptr<Polytope_Rn> polar_pol(new Polytope_Rn());
1460  TopGeomTools::PolarPolytope(pol, polar_pol, true, bb_size);
1461 
1462  // Step 4: compute the polar of the polar.
1463  pol.reset(new Polytope_Rn());
1464  TopGeomTools::PolarPolytope(polar_pol, pol);
1465 
1466  // Step 5: translate it back.
1467  TopGeomTools::Translate(pol, gravity_center);
1468 
1469  return 0;
1470 }
1471 
1472 //========================================================================================
1473 
1474 // set palette defined ( 0 "#A90D0D", 0.25 "#22A90D", 0.50 "#0DA9A9", 0.75 "#FFFFFF", 1.0 "#C0C0C0" )
1475 // .. # The background color is set first, then the border colors, then the X & Y axis colors, then the plotting colors
1476 // set terminal png size 1280,960 #ffffff #ffffff #ffffff #000000
1477 // set output "biarn.png"
1478 // plot "biarn.dat" using 3:4:(($5)/1000.):6 with circles linecolor palette z fill solid border lt 1 notitle, 'biarn.dat' using ($3):($4):1 with labels notitle;void Visualization::gnuplot2D(const boost::shared_ptr<Polytope_Rn>& polygon, const std::string& name, double col, std::ostream& out) throw (std::domain_error) {
1479 void Visualization::gnuplot2D(const boost::shared_ptr<Polytope_Rn>& polygon, const std::string& name, double col, std::ostream& out) throw (std::domain_error) {
1480  if (Rn::getDimension() != 2)
1481  throw std::domain_error("Visualization::gnuplot(std::ostream& out) dimension is not 2 ");
1482  out << "set object " << name << " polygon from ";
1483  boost::shared_ptr<HalfSpace_Rn> oldHS;
1484  boost::shared_ptr<Generator_Rn> startVTX,referenceVTX;
1485  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > LVX_A(polygon->getListOfGenerators());
1486  startVTX = LVX_A.current();
1487  std::set< boost::shared_ptr<Generator_Rn> > alreadyProcessedVTX;
1488  {for (LVX_A.begin(); LVX_A.end()!=true; LVX_A.next()) {
1489  if (alreadyProcessedVTX.find(LVX_A.current()) == alreadyProcessedVTX.end()) {
1490  referenceVTX = LVX_A.current();
1491  alreadyProcessedVTX.insert(referenceVTX);
1492  boost::shared_ptr<HalfSpace_Rn> referenceHS = referenceVTX->getFacet(0);
1493  if (referenceHS == oldHS)
1494  referenceHS = referenceVTX->getFacet(1);
1495  // Write the point
1496  out << LVX_A.current()->getCoordinate(0) << "," << LVX_A.current()->getCoordinate(1) << " to ";
1497  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > LVX_B(polygon->getListOfGenerators());
1498  {for (LVX_B.begin(); LVX_B.end()!=true; LVX_B.next()) {
1499  if (alreadyProcessedVTX.find(LVX_B.current()) == alreadyProcessedVTX.end()) {
1500  if (referenceHS == LVX_B.current()->getFacet(0)) {
1501  referenceHS = referenceVTX->getFacet(1);
1502  // Write the point
1503  out << LVX_B.current()->getCoordinate(0) << "," << LVX_B.current()->getCoordinate(1) << " to ";
1504  referenceVTX = LVX_B.current();
1505  alreadyProcessedVTX.insert(referenceVTX);
1506  }
1507  else if (referenceHS == LVX_B.current()->getFacet(1)) {
1508  referenceHS = referenceVTX->getFacet(0);
1509  // Write the point
1510  out << LVX_B.current()->getCoordinate(0) << "," << LVX_B.current()->getCoordinate(1) << " to ";
1511  referenceVTX = LVX_B.current();
1512  alreadyProcessedVTX.insert(referenceVTX);
1513  }
1514  }
1515  }}
1516  }
1517  }}
1518  out << startVTX->getCoordinate(0) << "," << startVTX->getCoordinate(1);
1519  out << " lc palette " << col << std::endl;
1520  //set object 1 fc rgb '#000000' fillstyle solid lw 0
1521 }
1522 
1523 
TopGeomTools::PolarPolytope
static int PolarPolytope(const boost::shared_ptr< Polytope_Rn > &original_pol, boost::shared_ptr< Polytope_Rn > &polar_pol, bool forceComputation=true, double bb_size=1000.)
Compute the polar polytope.
Definition: PolyhedralAlgorithms_Rn.cpp:1254
MinkowskiSum::_MinkowskiDecompositionOK
std::vector< bool > _MinkowskiDecompositionOK
Tell whether _MinkowskiDecomposition had to be considered or not.
Definition: PolyhedralAlgorithms_Rn.h:257
FaceEnumeration::printFacesWithVerticesToSage
void printFacesWithVerticesToSage(std::ostream &this_ostream) const
Definition: PolyhedralAlgorithms_Rn.cpp:223
Rn::setDimension
static polito_EXPORT void setDimension(unsigned int dim)
Set the dimension for the cartesian space we work in.
Definition: Rn.cpp:27
MinkowskiSum::_neighboursB
std::vector< std::vector< unsigned int > > _neighboursB
Store the neighbours of each vertex of B.
Definition: PolyhedralAlgorithms_Rn.h:248
FaceEnumeration::printFacesWithVertices
void printFacesWithVertices(std::ostream &this_ostream) const
Definition: PolyhedralAlgorithms_Rn.cpp:198
MinkowskiSum::MinkowskiSum
MinkowskiSum()
Sets up the data structure for the computation of the Minkowski sum of two polytopes.
Definition: PolyhedralAlgorithms_Rn.h:145
Polytope_Rn
Model a polytope using its two equivalent definitions : the convex hull and the half-space intersecti...
Definition: Polytope_Rn.h:34
FaceEnumeration::_allFacesWithFacets
std::vector< std::vector< ListOfFaces > > _allFacesWithFacets
Definition: PolyhedralAlgorithms_Rn.h:129
constIteratorOfListOfGeometricObjects::currentIteratorNumber
int currentIteratorNumber() const
Return the current position in the list.
Definition: GeometricObjectIterator_Rn.h:183
MinkowskiSum::processNormalFan1
void processNormalFan1()
Do the final job after having intersected all dual cones. The reduction process uses neighbourhood pr...
Definition: PolyhedralAlgorithms_Rn.cpp:722
HalfSpace_Rn
A half-space whose frontier is a linear (n-1) dimension space. _constant + _coefficients[0].x1 + ... + _coefficients[n-1].xn >= 0.
Definition: HalfSpace_Rn.h:37
PolyhedralCone_Rn
Model a polyhedral cone using its two equivalent definitions : the convex hull and the half-space int...
Definition: PolyhedralCone_Rn.h:45
MinkowskiSum
Compute the Minkowski sum of two polytopes.
Definition: PolyhedralAlgorithms_Rn.h:140
FaceEnumeration::_allFacesWithVertices
std::vector< std::vector< ListOfFaces > > _allFacesWithVertices
Definition: PolyhedralAlgorithms_Rn.h:131
FaceEnumeration::printFacesWithFacets
void printFacesWithFacets(std::ostream &this_ostream) const
Definition: PolyhedralAlgorithms_Rn.cpp:173
Visualization::gnuplot2D
static void gnuplot2D(const boost::shared_ptr< Polytope_Rn > &polygon, const std::string &name, double col, std::ostream &out)
Provide the drwing of polygon under the gnuplot format.
Definition: PolyhedralAlgorithms_Rn.cpp:1479
TopGeomTools::extrudeInCanonicalDirections
static int extrudeInCanonicalDirections(const std::set< unsigned int > &originalSpaceDirections, unsigned int dimensionOfTotalSpace, const boost::shared_ptr< Polytope_Rn > &original_polytope, boost::shared_ptr< PolyhedralCone_Rn > &extruded_polyhedron)
Compute the extrusion of a polytope belonging to a space S1 into a higher dimension space S2 followin...
Definition: PolyhedralAlgorithms_Rn.cpp:1395
constIteratorOfListOfGeometricObjects::current
const GEOMETRIC_OBJECT current()
Return the current geometric element.
Definition: GeometricObjectIterator_Rn.h:174
TopGeomTools::scalingFactor
static int scalingFactor(boost::shared_ptr< Polytope_Rn > &pol, double factor)
Multiply a polytope by a given real number.
Definition: PolyhedralAlgorithms_Rn.cpp:1234
PseudoSumWithoutCaps::rebuildSum
boost::shared_ptr< Polytope_Rn > rebuildSum(const std::set< unsigned int > &firstOperandCaps, const std::set< unsigned int > &secondOperandCaps, std::set< unsigned int > &newCaps, double bb_size=1000.)
Remove the cap half-spaces stored in sets and then truncate again.
Definition: PolyhedralAlgorithms_Rn.cpp:1078
MinkowskiSum::_sum
boost::shared_ptr< Polytope_Rn > _sum
Definition: PolyhedralAlgorithms_Rn.h:232
MinkowskiSum::_A2C
std::vector< std::vector< unsigned int > > _A2C
Store the polyhedrical cap in C of each vertex of A.
Definition: PolyhedralAlgorithms_Rn.h:239
ListOfFaces
std::vector< unsigned int > ListOfFaces
Definition: PolyhedralAlgorithms_Rn.h:40
TopGeomTools::GravityCenter
static int GravityCenter(boost::shared_ptr< Polytope_Rn > &pol, boost::numeric::ublas::vector< double > &gravity_center)
Compute the gravity center of a given polytope.
Definition: PolyhedralAlgorithms_Rn.cpp:1219
MinkowskiSum::_firstOperand
const boost::shared_ptr< Polytope_Rn > _firstOperand
Definition: PolyhedralAlgorithms_Rn.h:230
constIteratorOfListOfGeometricObjects::begin
void begin()
Move the iterator at the beginning of the list.
Definition: GeometricObjectIterator_Rn.h:150
MinkowskiSum::processNormalFan2
void processNormalFan2()
Do the final job after having intersected all dual cones. The reduction process builds half-spaces an...
Definition: PolyhedralAlgorithms_Rn.cpp:584
MinkowskiSum::_neighboursA
std::vector< std::vector< unsigned int > > _neighboursA
Store the neighbours of each vertex of A.
Definition: PolyhedralAlgorithms_Rn.h:245
PolyhedralAlgorithms_Rn.h
NormalFan_Rn.h
FaceEnumeration::ComputeWithFacets
static void ComputeWithFacets(const boost::shared_ptr< Polytope_Rn > &A, FaceEnumeration &FaceEnum)
Definition: PolyhedralAlgorithms_Rn.cpp:51
Rn::getTolerance
static polito_EXPORT double getTolerance()
Give the minimum distance between two points.
Definition: Rn.cpp:31
MinkowskiSum::_secondOperand
const boost::shared_ptr< Polytope_Rn > _secondOperand
Definition: PolyhedralAlgorithms_Rn.h:231
FaceEnumeration
Combinatorial face enumeration for polytopes.
Definition: PolyhedralAlgorithms_Rn.h:44
MinkowskiSum::_B2C
std::vector< std::vector< unsigned int > > _B2C
Store the polyhedrical cap in C of each vertex of B.
Definition: PolyhedralAlgorithms_Rn.h:241
FaceEnumeration::save
static void save(const std::string &filename, const std::vector< std::vector< ListOfFaces > > &latt)
Save the polytope lattice.
Definition: PolyhedralAlgorithms_Rn.cpp:275
Rn::getDimension
static polito_EXPORT unsigned int getDimension()
Return the dimension of the cartesian space we work in.
Definition: Rn.cpp:29
Rn.h
StrongRedundancyProcessing
This class can be more time-consuming than WeakRedundancyProcessing or NoRedundancyProcessing because...
Definition: DoubleDescription_Rn.h:514
Generator_Rn
A n-coordinates generator, which can be a vertex or an edge whether it is contained by a polytope or ...
Definition: Generator_Rn.h:38
FaceEnumeration::ComputeWithVertices
static void ComputeWithVertices(const boost::shared_ptr< Polytope_Rn > &A, FaceEnumeration &FaceEnum)
Definition: PolyhedralAlgorithms_Rn.cpp:124
constIteratorOfListOfGeometricObjects
This class is designed to run the list of all geometric objects representing a polytope.
Definition: GeometricObjectIterator_Rn.h:142
IO_Polytope.h
Polytope_Rn.h
MinkowskiSum::_NF_Cones
std::vector< boost::shared_ptr< PolyhedralCone_Rn > > _NF_Cones
The normal fan polyhedrical cones list.
Definition: PolyhedralAlgorithms_Rn.h:260
MinkowskiSum::computeCommonRefinement
void computeCommonRefinement(const std::vector< boost::shared_ptr< Generator_Rn > > &listOfGenerators_A, const std::vector< boost::shared_ptr< Generator_Rn > > &listOfGenerators_B, std::vector< boost::shared_ptr< Generator_Rn > > &listOfGenerators_C, const std::vector< boost::shared_ptr< PolyhedralCone_Rn > > &listOfDualCones_A, const std::vector< boost::shared_ptr< PolyhedralCone_Rn > > &listOfDualCones_B, std::vector< boost::shared_ptr< PolyhedralCone_Rn > > &listOfDualCones_C)
Just compute the common refinement of the two operands normal fans.
Definition: PolyhedralAlgorithms_Rn.cpp:398
TopGeomTools::projectPolytopeOnCanonicalHyperplanes
static int projectPolytopeOnCanonicalHyperplanes(const std::set< unsigned int > &listOfHyperplanes, const boost::shared_ptr< Polytope_Rn > &original_pol, boost::shared_ptr< Polytope_Rn > &proj_pol)
Compute the projection of a polytope on the intersection of canonical hyperplanes of the shape xi = 0...
Definition: PolyhedralAlgorithms_Rn.cpp:1344
MinkowskiSum::_MinkowskiDecomposition
std::vector< std::pair< unsigned int, unsigned int > > _MinkowskiDecomposition
Store the genitors in A and B of each vertex of C.
Definition: PolyhedralAlgorithms_Rn.h:255
FaceEnumeration::Compute
static void Compute(const boost::shared_ptr< Polytope_Rn > &A)
Definition: PolyhedralAlgorithms_Rn.cpp:40
PseudoSumWithoutCaps::computeCapHalfSpaces
void computeCapHalfSpaces(const std::set< unsigned int > &firstOperandCaps, const std::set< unsigned int > &secondOperandCaps, std::set< unsigned int > &sumCaps)
Return the cap half-spaces of the sum in function of the two operands cap half-spaces.
Definition: PolyhedralAlgorithms_Rn.cpp:886
TopGeomTools::Translate
static int Translate(boost::shared_ptr< Polytope_Rn > &pol, const boost::numeric::ublas::vector< double > &v2t)
Translate a polytope by the given vector.
Definition: PolyhedralAlgorithms_Rn.cpp:1199
PseudoIntersectionWithoutCaps::PseudoIntersectionWithoutCaps
PseudoIntersectionWithoutCaps(const boost::shared_ptr< Polytope_Rn > &A, const boost::shared_ptr< Polytope_Rn > &B, boost::shared_ptr< Polytope_Rn > &C, const std::set< unsigned int > &firstOperandCaps, const std::set< unsigned int > &secondOperandCaps, std::set< unsigned int > &newCaps, double bb_size=1000.)
Remove all cap half-spaces and then compute the intersection of two capped polytopes.
Definition: PolyhedralAlgorithms_Rn.cpp:1141
DoubleDescriptionFromGenerators::Compute
static int Compute(boost::shared_ptr< Polytope_Rn > &pol, double bb_size=1000.)
Use the polarity to get the facets from the generators.
Definition: PolyhedralAlgorithms_Rn.cpp:1436
MinkowskiSum::compute
boost::shared_ptr< Polytope_Rn > compute()
Compute the common refinement of the two operands normal fans and then build the sum.
Definition: PolyhedralAlgorithms_Rn.cpp:529
MinkowskiSum::processNormalFan0
void processNormalFan0()
Do the final job after having intersected all dual cones. The reduction process simply compares all d...
Definition: PolyhedralAlgorithms_Rn.cpp:849
DoubleDescription
The algorithm implemented here is an incremental algorithm as mentioned in How Good are Convex Hull ...
Definition: DoubleDescription_Rn.h:49
FaceEnumeration::load
static void load(const std::string &filename, std::vector< std::vector< ListOfFaces > > &latt)
Load the polytope lattice.
Definition: PolyhedralAlgorithms_Rn.cpp:262
MinkowskiSum::_NF_Vertices
std::vector< boost::shared_ptr< Generator_Rn > > _NF_Vertices
The list of C vertices.
Definition: PolyhedralAlgorithms_Rn.h:262