politopix  5.0.0
DoubleDescription_Rn.h
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) 2012-2022 : Delos Vincent
3 //
4 // This program is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU 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 General Public License for more details.
13 //
14 // You should have received a copy of the GNU 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 #ifndef DOUBLEDESCRIPTION_Rn
22 #define DOUBLEDESCRIPTION_Rn
23 
24 #include <math.h>
25 #include <numeric>
26 #include <iterator>
27 #include <set>
28 #include <map>
29 //#include <boost/unordered_set.hpp>
30 #include "Tracking.h"
31 #include "Neighbours_Rn.h"
32 #include "UpdDoubleDescription_Rn.h"
33 
34 
36 class Segment_Rn {
37 
38 public:
39  Segment_Rn() {_isActive = true;}
40  Segment_Rn(const std::vector< unsigned int >& LC) {_isActive = true;}
41  Segment_Rn(boost::shared_ptr< Generator_Rn_SD>& e1, boost::shared_ptr< Generator_Rn_SD>& e2):_firstEnd(e1),_secondEnd(e2) {_isActive = true;}
42  Segment_Rn(boost::shared_ptr< Generator_Rn_SD>& e1, boost::shared_ptr< Generator_Rn_SD>& e2, const std::vector< unsigned int >& LC):_firstEnd(e1),_secondEnd(e2) {_isActive = true;}
43  boost::shared_ptr< Generator_Rn_SD >& getFirstGenerator() {return _firstEnd;}
44  boost::shared_ptr< Generator_Rn_SD >& getSecondGenerator() {return _secondEnd;}
45  void setFirstGenerator( const boost::shared_ptr< Generator_Rn_SD>& e1) {_firstEnd = e1;}
46  void setSecondGenerator(const boost::shared_ptr< Generator_Rn_SD>& e2) {_secondEnd = e2;}
47  bool getActive() {return _isActive;}
48  void setActive() {_isActive = true;}
49  void setInactive() {_isActive = false;}
50  void dump(std::ostream &ofs) {
51  ofs << "[[ "; (_isActive == true) ? ofs << "t, " : ofs << "f, ";
52  getFirstGenerator()->dump(ofs); std::cout << " , ";
53  getSecondGenerator()->dump(ofs);
54  //std::copy(_commonLinearConstraints.begin(), _commonLinearConstraints.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
55  ofs << " ]]";
56  }
57 
58 protected:
60  //std::vector< unsigned int > _commonLinearConstraints;
62  boost::shared_ptr< Generator_Rn_SD > _firstEnd;
64  boost::shared_ptr< Generator_Rn_SD > _secondEnd;
66  bool _isActive;
67 };
68 
86 template< class POLYHEDRON, class ITERATOR > class DoubleDescription {
87 
88  public:
89  DoubleDescription(POLYHEDRON poly, ITERATOR ite, int truncationStep) {
90  _TOL = Rn::getTolerance();
91  _RnDIM = poly->dimension();
92  std::vector< boost::shared_ptr<Generator_Rn_SD> > listOfGenSD;
93  // Generator_Rn => Generator_Rn_SD
94  // From the current generators make a list of generators dedicated to work in the data structure algorithm.
95  // Make sure the list of half-space numbers is sorted as we want to use algorithms further.
96  poly->getListOfGeneratorsSD(listOfGenSD);
97  unsigned int nbProcessHyp = init(poly, ite, truncationStep);
98  // Do the hard job now.
99  _isEmpty = !compute(poly, ite, truncationStep, listOfGenSD);
100 
101  if (_isEmpty == false) {
102  // Set the pointers from the generators to the half-spaces before dealing with redundancy - as the redundancy
103  // processing algorithm will change the numbers of the half-spaces. Here recycle the old pointers stored in
104  // the polyhedron to build the new list with the created generators.
105 #ifdef DEBUG
106  std::cout << "Entering poly->setListOfGeneratorsSD(listOfGenSD) ..." << std::endl;
107 #endif
108  poly->setListOfGeneratorsSD(listOfGenSD);
109 #ifdef DEBUG
110  std::cout << "Leaving poly->setListOfGeneratorsSD(listOfGenSD) ..." << std::endl;
111  std::cout << std::endl << "// Display the final list of generators." << std::endl;
112  dumpListOfGenerators(poly, listOfGenSD, std::cout);
113 #endif
114 
115  // Include all the generators now.
116  std::vector< unsigned int > HS2Remove;
117  {
118  std::set< unsigned int >::const_iterator it;
119  for (it=_listOfRedundantLinearConstraints.begin(); it!=_listOfRedundantLinearConstraints.end(); ++it)
120  HS2Remove.push_back(*it);
121  }
122 #ifdef DEBUG
123  std::cout << "Before sort() ..." << std::endl;
124  {
125  for (unsigned i=0; i<poly->numberOfHalfSpaces(); ++i) {
126  poly->getHalfSpace(i)->dump(std::cout);
127  std::cout << std::endl;
128  }
129  }
130  std::cout << std::endl;
131 #endif
132  std::sort(HS2Remove.begin(), HS2Remove.end(), std::greater<unsigned int>());
133  {
134  std::vector< unsigned int >::const_iterator it2;
135  for (it2=HS2Remove.begin(); it2!=HS2Remove.end(); ++it2) {
136  //std::cout << "Delete half-space " << *it2;
137  //poly->getHalfSpace(*it2)->dump(std::cout);
138  // Do not forget to begin by the bigger to remove as we could change the ordering.
139  poly->removeHalfSpace(*it2);
140  //std::cout << std::endl;
141  }
142  }
143 #ifdef DEBUG
144  std::cout << "Display HS ..." << std::endl;
145  {
146  for (unsigned i=0; i<poly->numberOfHalfSpaces(); ++i) {
147  poly->getHalfSpace(i)->dump(std::cout);
148  std::cout << std::endl;
149  }
150  }
151  std::cout << std::endl;
152  std::cout << "Before unhookRedundantGenerators() ..." << std::endl;
153  std::cout << std::endl << "// Display the final list of generators." << std::endl;
154  dumpListOfGenerators(poly, listOfGenSD, std::cout, false);
155 #endif
156  //poly->dump(std::cout);
157  // In some fuzzy cases, computing the intersection between the space spanned by two generators and a hyperplane
158  // provides a fuzzy generator somewhere inside the fuzzy zone. The location of such a generator needs to be
159  // refined later on in the process. Meanwhile the current one is likely to lose half-spaces as the probability
160  // that some of them are redundant is high. So remove them as they are not genuine generators.
161  // Case test/DATA1/P38_R6.ptop ./trunk/test_politopix.cpp line 274 revision 409
162  unhookRedundantGenerators(poly);
163 #ifdef DEBUG
164  std::cout << "After unhookRedundantGenerators() ..." << std::endl;
165 #endif
166  }
167  }
168 
171  DoubleDescription(POLYHEDRON poly, ITERATOR ite, int truncationStep, TrackingOperatorToResult& trackerVdesc, TrackingOperatorToResult& trackerHdesc) {
172  _TOL = Rn::getTolerance();
173  _RnDIM = poly->dimension();
174  std::vector< boost::shared_ptr<Generator_Rn_SD> > listOfGenSD;
175  // Generator_Rn => Generator_Rn_SD
176  // From the current generators make a list of generators dedicated to work in the data structure algorithm.
177  // Make sure the list of half-space numbers is sorted as we want to use algorithms further.
178  poly->getListOfGeneratorsSD(listOfGenSD);
179  unsigned int nbProcessHyp = init(poly, ite, truncationStep);
180  // Do the hard job now.
181  _isEmpty = !compute(poly, ite, truncationStep, listOfGenSD);
182 
183  if (_isEmpty == false) {
184  // Store the modifications in the binary operation tracker.
185  // Operator1_Entity[i] = (Generator_Rn_SD::Status, Result_Entity[j])
186  // Operator2_Entity[k] = (Generator_Rn_SD::Status, Result_Entity[l])
187  // Result_Entity[m] = (Generator_Rn_SD::Status, (Operator1_Entity[p], Operator2_Entity[q]))
188  unsigned int nbRes=0;
189  std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGenSD;
190  {
191  for (iteGenSD=listOfGenSD.begin(); iteGenSD!=listOfGenSD.end(); ++iteGenSD) {
192  // In the tracker, the operand entities have been set to Operator_DELETED by default and the
193  // result entities have been set by default to Result_UNCHANGED.
194  Generator_Rn_SD::Status thisStatus = (*iteGenSD)->getStatus();
195  unsigned int nbOp1 = (*iteGenSD)->getGeneratorNumber();
196  if (thisStatus == Generator_Rn_SD::CREATED || thisStatus == Generator_Rn_SD::CREATED_AND_MODIFIED) {
197  trackerVdesc.setResultEntityAsCreated(nbRes);
198  trackerVdesc.setResult_Operator(nbRes,-1);
199  }
200  else if (thisStatus == Generator_Rn_SD::MODIFIED) {
201  trackerVdesc.setResultEntityAsModified(nbRes);
202  trackerVdesc.setResult_Operator(nbRes,nbOp1);
203  trackerVdesc.setOperatorEntityAsModified(nbOp1);
204  trackerVdesc.setOperator_Result(nbOp1,nbRes);
205  }
206  else if (thisStatus == Generator_Rn_SD::UNCHANGED) {
207  trackerVdesc.setResultEntityAsUnchanged(nbRes);
208  trackerVdesc.setResult_Operator(nbRes,nbOp1);
209  trackerVdesc.setOperatorEntityAsUnchanged(nbOp1);
210  trackerVdesc.setOperator_Result(nbOp1,nbRes);
211  }
212  ++nbRes;
213  }
214  }
215  // Set the pointers from the generators to the half-spaces before dealing with redundancy - as the
216  // redundancy processing algorithm will change the numbers of the half-spaces.
217  poly->setListOfGeneratorsSD(listOfGenSD);
218  // Process the half-spaces now.
219  std::vector< unsigned int > HS2Remove;
220  {
221  std::set< unsigned int >::const_iterator it;
222  for (it=_listOfRedundantLinearConstraints.begin(); it!=_listOfRedundantLinearConstraints.end(); ++it)
223  HS2Remove.push_back(*it);
224  }
225  std::sort(HS2Remove.begin(), HS2Remove.end(), std::greater<unsigned int>());
226  {
227  std::vector< unsigned int >::const_iterator it2;
228  for (it2=HS2Remove.begin(); it2!=HS2Remove.end(); ++it2) {
229  //std::cout << "Delete half-space " << *it2;
230  //poly->getHalfSpace(*it2)->dump(std::cout);
231  // Do not forget to begin by the bigger to remove as we could change the ordering.
232  poly->removeHalfSpace(*it2);
233  //std::cout << std::endl;
234  }
235  }
236  //poly->dump(std::cout);
237  // In some fuzzy cases, computing the intersection between the space spanned by two generators and a hyperplane
238  // provides a fuzzy generator somewhere inside the fuzzy zone. The location of such a generator needs to be
239  // refined later on in the process. Meanwhile the current one is likely to lose half-spaces as the probability
240  // that some of them are redundant is high. So remove them as they are not genuine generators.
241  unhookRedundantGenerators(poly);
242  }
243  }
244 
245  virtual ~DoubleDescription() {}
246 
247  unsigned int init(POLYHEDRON poly, ITERATOR ite, int truncationStep) {
248  // For polyhedral cones neigbourhoodCondition() is Rn::getDimension()-1, for polytopes it is Rn::getDimension()-2.
249  _neigbourhoodCondition = poly->neigbourhoodCondition();
250  // Count the number of hyperplanes versus the number of half-spaces.
251  unsigned int nbHyp = 0, nbProcessHyp = 0;
252  {
253  for (ite.begin(); ite.end()!=true; ite.next()) {
254  const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace = ite.current();
255  boost::shared_ptr<Hyperplane_Rn> currentHyp = boost::dynamic_pointer_cast<Hyperplane_Rn>(currentHalfSpace);
256  if (currentHyp) {
257  nbHyp++;
258  if (ite.currentIteratorNumber() < truncationStep)
259  nbProcessHyp++;
260  }
261  double halfSpaceNorm = std::inner_product(currentHalfSpace->begin(), currentHalfSpace->end(), currentHalfSpace->begin(), 0.);
262  // The array of generators attached to each half-space.
263  //std::unordered_set< boost::shared_ptr< Generator_Rn_SD > > genSet;
264  std::set< boost::shared_ptr< Generator_Rn_SD > > genSet;
265  _listOfGeneratorsPerLinearConstraint.push_back(genSet);
266  // Initialize _listOfLinearConstraints and _listOfLinearConstraintNorms.
267  _listOfLinearConstraints.push_back(currentHalfSpace);
268  _listOfLinearConstraintNorms.push_back(sqrt(halfSpaceNorm));
269  boost::shared_ptr<NormedHalfSpace_Rn> currentNormedHalfSpace;
270  currentNormedHalfSpace.reset( new NormedHalfSpace_Rn(*currentHalfSpace) );
271  _listOfNormedLinearConstraints.push_back(currentNormedHalfSpace);
272  }
273  }
274  return nbProcessHyp;
275  }
276 
277  bool getIsEmpty() const {return _isEmpty;}
278 
280  void computeStatesOfAllGenerators(std::vector< boost::shared_ptr<Generator_Rn_SD> >& arrayOfGenOUT, const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace,
281  unsigned int halfspaceNumber, unsigned int& nbGeneratorsIn, unsigned int& nbGeneratorsOut, unsigned int& nbGeneratorsOn) {
282  nbGeneratorsIn = 0;
283  nbGeneratorsOn = 0;
284  nbGeneratorsOut = 0;
285  double halfSpaceNorm = std::inner_product(currentHalfSpace->begin(), currentHalfSpace->end(), currentHalfSpace->begin(), 0.);
286  halfSpaceNorm = sqrt(halfSpaceNorm);
287  //_listOfLinearConstraintNorms.push_back(halfSpaceNorm);
288  //_listOfLinearConstraints.push_back(currentHalfSpace);
289 
290  // First of all reset the state to hs_UNKNOWN in the generators.
291  {
292  std::vector< Segment_Rn >::iterator it_GN;
293  for (it_GN = _allNeighbours.begin(); it_GN != _allNeighbours.end(); ++it_GN) {
294  it_GN->getFirstGenerator()->resetState();
295  it_GN->getSecondGenerator()->resetState();
296  }
297  }
298  // Compute all states of all generators with respect to the half-space.
299 #ifdef DEBUG
300  std::vector< boost::shared_ptr<Generator_Rn_SD> > allCurrentGenerators;
301 #endif
302  {
303  HalfSpace_Rn::State currentState = HalfSpace_Rn::hs_UNKNOWN;
304  std::vector< Segment_Rn >::iterator it_GN;
305  for (it_GN = _allNeighbours.begin(); it_GN != _allNeighbours.end(); ++it_GN) {
306  if (it_GN->getActive() == true) {
307 #ifdef DEBUG
308  if (it_GN->getFirstGenerator()->getState() == HalfSpace_Rn::hs_UNKNOWN)
309  allCurrentGenerators.push_back( it_GN->getFirstGenerator() );
310  if (it_GN->getSecondGenerator()->getState() == HalfSpace_Rn::hs_UNKNOWN)
311  allCurrentGenerators.push_back( it_GN->getSecondGenerator() );
312 #endif
313  if (it_GN->getFirstGenerator()->getState() == HalfSpace_Rn::hs_UNKNOWN) {
314  if (computeGeneratorState(it_GN->getFirstGenerator(), currentHalfSpace, halfspaceNumber, halfSpaceNorm, nbGeneratorsIn, nbGeneratorsOut, nbGeneratorsOn) == HalfSpace_Rn::hs_OUT)
315  arrayOfGenOUT.push_back(it_GN->getFirstGenerator());
316  }
317  if (it_GN->getSecondGenerator()->getState() == HalfSpace_Rn::hs_UNKNOWN) {
318  if (computeGeneratorState(it_GN->getSecondGenerator(), currentHalfSpace, halfspaceNumber, halfSpaceNorm, nbGeneratorsIn, nbGeneratorsOut, nbGeneratorsOn) == HalfSpace_Rn::hs_OUT)
319  arrayOfGenOUT.push_back(it_GN->getSecondGenerator());
320  }
321  }
322  }
323  }
324 #ifdef DEBUG
325  {
326  std::sort( allCurrentGenerators.begin(), allCurrentGenerators.end(), &inferior );
327  std::vector< boost::shared_ptr<Generator_Rn_SD> >::iterator it_GN0;
328  for (it_GN0 = allCurrentGenerators.begin(); it_GN0 != allCurrentGenerators.end(); ++it_GN0) {
329  boost::shared_ptr<Generator_Rn_SD> Gen = *it_GN0;
330  std::cout.precision(15);
331  //std::cout << "# V" << Gen->getGeneratorNumber() << " = [";
332  std::cout << "# G" << Gen->getGeneratorNumber() << " = [";
333  Gen->save(std::cout);
334  std::cout << "]" << std::endl << "{ ";
335  {
336  for (unsigned int ii=0; ii<Gen->numberOfFacets(); ii++) {
337  std::cout << Gen->getFacet(ii);
338  if (ii != Gen->numberOfFacets()-1)
339  std::cout << " ";
340  }
341  }
342  double scalarProduct = Gen->getScalarProduct();
343  std::cout << "}" << std::endl;
344  std::cout << "dotP = " << scalarProduct;
345  std::cout << ", dist = " << (scalarProduct+currentHalfSpace->getConstant()) / halfSpaceNorm;
346  //std::cout << "sp = " << scalarProduct << " ";
347  double distanceToHyperplane = (scalarProduct+currentHalfSpace->getConstant()) / halfSpaceNorm;
348  HalfSpace_Rn::State currentState = HalfSpace_Rn::hs_UNKNOWN;
349  if (distanceToHyperplane > _TOL)
350  currentState = HalfSpace_Rn::hs_IN;
351  else if (distanceToHyperplane < -_TOL)
352  currentState = HalfSpace_Rn::hs_OUT;
353  else
354  currentState = HalfSpace_Rn::hs_ON;
355  std::cout << ", state = " << HalfSpace_Rn::getStateAsText(currentState) << std::endl;
356  }
357  }
358 #endif
359  }
360 
362  static bool inferior(const boost::shared_ptr<Generator_Rn_SD>& HS1, const boost::shared_ptr<Generator_Rn_SD>& HS2) {
363  unsigned int i=0;
364  unsigned int n=HS1->dimension();
365  do {
366  if (HS1->getCoordinate(i) < HS2->getCoordinate(i))
367  return true;
368  else if (HS1->getCoordinate(i) > HS2->getCoordinate(i))
369  return false;
370  else
371  i++;
372  } while (i < n);
373  return true;
374  }
375 
377  HalfSpace_Rn::State computeGeneratorState(const boost::shared_ptr<Generator_Rn_SD>& Gen, const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace, unsigned int halfspaceNumber,
378  double halfSpaceNorm, unsigned int& nbGeneratorsIn, unsigned int& nbGeneratorsOut, unsigned int& nbGeneratorsOn) {
379  double scalarProduct = std::inner_product(Gen->begin(), Gen->end(), currentHalfSpace->begin(), 0.);
380  //boost::numeric::ublas::inner_prod(currentGenerator->vect(), currentHalfSpace->vect());
381 #ifdef DEBUG
382  // std::cout.precision(15);
383  // std::cout << "# V" << Gen->getGeneratorNumber() << " = [";
384  // Gen->save(std::cout);
385  // std::cout << "]" << std::endl << "{ ";
386  // {
387  // for (unsigned int ii=0; ii<Gen->numberOfFacets(); ii++) {
388  // std::cout << Gen->getFacet(ii);
389  // if (ii != Gen->numberOfFacets()-1)
390  // std::cout << " ";
391  // }
392  // }
393  // std::cout << "}" << std::endl;
394  // std::cout << "dotP = " << scalarProduct;
395  // std::cout << ", dist = " << (scalarProduct+currentHalfSpace->getConstant()) / halfSpaceNorm;
396 #endif
397  //std::cout << "sp = " << scalarProduct << " ";
398  double distanceToHyperplane = (scalarProduct+currentHalfSpace->getConstant()) / halfSpaceNorm;
399  if (distanceToHyperplane > _TOL) {
400  nbGeneratorsIn++;
401  Gen->setState(HalfSpace_Rn::hs_IN, scalarProduct);
402  return HalfSpace_Rn::hs_IN;
403  }
404  else if (distanceToHyperplane < -_TOL) {
405  nbGeneratorsOut++;
406  Gen->setState(HalfSpace_Rn::hs_OUT, scalarProduct);
407  return HalfSpace_Rn::hs_OUT;
408  }
409  else {
410  nbGeneratorsOn++;
411  Gen->setState(HalfSpace_Rn::hs_ON, scalarProduct);
412  return HalfSpace_Rn::hs_ON;
413  }
414 #ifdef DEBUG
415  // HalfSpace_Rn::State currentState = HalfSpace_Rn::hs_UNKNOWN;
416  // if (distanceToHyperplane > _TOL)
417  // currentState = HalfSpace_Rn::hs_IN;
418  // else if (distanceToHyperplane < -_TOL)
419  // currentState = HalfSpace_Rn::hs_OUT;
420  // else
421  // currentState = HalfSpace_Rn::hs_ON;
422  // std::cout << ", state = " << HalfSpace_Rn::getStateAsText(currentState) << std::endl;
423 #endif
424  }
425 
427  bool findBestLinearConstraint(std::vector< boost::shared_ptr<HalfSpace_Rn> >::iterator currentLC, std::vector< boost::shared_ptr<HalfSpace_Rn> >::iterator bestLC) {
428  double bestParameter = std::numeric_limits<double>::max();
429  unsigned int hsNb = currentLC - _listOfLinearConstraints.begin();
430  for (; currentLC != _listOfLinearConstraints.end(); ++currentLC, ++hsNb) {
431  double halfSpaceNorm = _listOfLinearConstraintNorms[hsNb];
432  // std::cout << "[0.1]" << std::endl;
433  const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace = _listOfLinearConstraints[hsNb];
434  //std::cout << "[0.2]" << std::endl;
435  boost::shared_ptr<Generator_Rn_SD>& Gen1 = _allNeighbours[ _currentSegmentIndexAllNeighbours ].getFirstGenerator();
436  //std::cout << "[0.3]" << std::endl;
437  boost::shared_ptr<Generator_Rn_SD>& Gen2 = _allNeighbours[ _currentSegmentIndexAllNeighbours ].getSecondGenerator();
438  //std::cout << "[0.4] " << _currentSegmentIndexAllNeighbours << " " << _allNeighbours.size() << std::endl;
439  //(_allNeighbours[ _currentSegmentIndexAllNeighbours ].getActive() == true) ? std::cout << "t" << std::endl : std::cout << "f" << std::endl;
440  //if (currentHalfSpace->getCoefficient(0) == -10) {
441  //}
442  //currentHalfSpace->dump(std::cout);
443  //std::cout << "[0.41]" << std::endl;
444  //Gen1->save(std::cout);
445  //std::cout << "[0.42]" << std::endl;
446  //Gen2->save(std::cout);
447  double scalarProduct1 = std::inner_product(Gen1->begin(), Gen1->end(), currentHalfSpace->begin(), 0.);
448  //std::cout << "[0.5]" << std::endl;
449  double scalarProduct2 = std::inner_product(Gen2->begin(), Gen2->end(), currentHalfSpace->begin(), 0.);
450  //std::cout << "[0.6]" << std::endl;
451  HalfSpace_Rn::State state1 = HalfSpace_Rn::hs_ON, state2 = HalfSpace_Rn::hs_ON;
452  //std::cout << "[0.7]" << std::endl;
453  double distanceToHyperplane1 = (scalarProduct1+currentHalfSpace->getConstant()) / halfSpaceNorm;
454  //std::cout << "[0.8]" << std::endl;
455  double distanceToHyperplane2 = (scalarProduct2+currentHalfSpace->getConstant()) / halfSpaceNorm;
456  //std::cout << "[1]" << std::endl;
457  if (distanceToHyperplane1 > _TOL)
458  state1 = HalfSpace_Rn::hs_IN;
459  else if (distanceToHyperplane1 < -_TOL)
460  state1 = HalfSpace_Rn::hs_OUT;
461  if (distanceToHyperplane2 > _TOL)
462  state2 = HalfSpace_Rn::hs_IN;
463  else if (distanceToHyperplane2 < -_TOL)
464  state2 = HalfSpace_Rn::hs_OUT;
465  //std::cout << "State1 = " << HalfSpace_Rn::getStateAsText(state1) << ", State2 = " << HalfSpace_Rn::getStateAsText(state2) << std::endl;
466  // Make sure we don't have a ON-ON configuration as we need to divide by (scalarProduct1 - scalarProduct2).
467  if (state1 != HalfSpace_Rn::hs_ON || state2 != HalfSpace_Rn::hs_ON) {
468  if (state1 == HalfSpace_Rn::hs_OUT && state2 == HalfSpace_Rn::hs_OUT) {
469  // We know the segment corresponding to currentSegmentIndex in _allNeighbours will be discarded by the truncation process.
470  bestLC = currentLC;
471  //std::cout << "OUT-OUT" << std::endl;
472  //std::cout << "[2] OUT-OUT" << std::endl;
473  return true;
474  }
475  else if (state1 == HalfSpace_Rn::hs_IN && state2 == HalfSpace_Rn::hs_OUT) {
476  // t = (Hcst + scalarProduct1) / (scalarProduct1 - scalarProduct2)
477  // newGen = t.Gen2 + (1.-t).Gen1
478  double param = (currentHalfSpace->getConstant() + scalarProduct1) / (scalarProduct1 - scalarProduct2);
479  //std::cout << "param12 = " << param;
480  if (param < bestParameter) {
481  bestParameter = param;
482  bestLC = currentLC;
483  }
484  //std::cout << "[2] IN-OUT" << std::endl;
485  //std::cout << ", bestParameter12 = " << bestParameter << std::endl;
486  }
487  else if (state1 == HalfSpace_Rn::hs_OUT && state2 == HalfSpace_Rn::hs_IN) {
488  // t = (Hcst + scalarProduct2) / (scalarProduct2 - scalarProduct1)
489  // newGen = t.Gen1 + (1.-t).Gen2
490  double param = (currentHalfSpace->getConstant() + scalarProduct2) / (scalarProduct2 - scalarProduct1);
491  //std::cout << "param21 = " << param;
492  if (param < bestParameter) {
493  bestParameter = param;
494  bestLC = currentLC;
495  }
496  //std::cout << "[2] OUT-IN" << std::endl;
497  //std::cout << ", bestParameter21 = " << bestParameter << std::endl;
498  }
499  }
500  } // for (; currentLC != _listOfLinearConstraints.end(); ++currentLC, ++hsNb) {
501  //std::cout << "[3]" << std::endl;
502  return true;
503  }
504 
506  bool compute(POLYHEDRON poly, ITERATOR iteHS, int truncationStep, std::vector< boost::shared_ptr<Generator_Rn_SD> >& listOfGenSD) {
507  // Store the scalar product result.
508  // Store generators IN, ON, OUT or newly created.
509  std::vector< boost::shared_ptr<Generator_Rn_SD> > GN_IN;
510  std::vector< boost::shared_ptr<Generator_Rn_SD> > GN_OUT;
511  std::vector< boost::shared_ptr<Generator_Rn_SD> > GN_ON;
512  std::vector< boost::shared_ptr<Generator_Rn_SD> > GN_new;
513  std::vector< boost::shared_ptr<Generator_Rn_SD> > listOfGeneratorsOUT2recycle;
514  int halfspaceNumber = truncationStep, generatorNumber = 0;
515  _currentSegmentIndexAllNeighbours = 0;
516  double TOL2 = _TOL * _TOL;
517  std::vector< unsigned int > listOfPairsToRemove;
518 
519 #ifdef DEBUG
520  std::cout << std::endl << "DoubleDescription::compute(" << truncationStep << ")" << std::endl;
521 #endif
522  // Attach the generators to their linear constraints.
523  {
524  std::vector< boost::shared_ptr<Generator_Rn_SD> >::iterator iteGN_1;
525  for (iteGN_1 = listOfGenSD.begin(); iteGN_1 != listOfGenSD.end(); ++iteGN_1) {
526  //boost::shared_ptr< Generator_Rn_SD > tt;
527  //_listOfGeneratorsPerLinearConstraint[ 0 ].push_back(tt);
528  for (unsigned int i = 0; i < (*iteGN_1)->numberOfFacets(); ++i)
529  _listOfGeneratorsPerLinearConstraint[ (*iteGN_1)->getFacet(i) ].insert( *iteGN_1 );
530  // std::vector< std::unordered_set< boost::shared_ptr< Generator_Rn_SD > > > _listOfGeneratorsPerLinearConstraint;
531  // std::vector< boost::shared_ptr<Generator_Rn_SD> >& listOfGenSD
532  }
533  }
534 #ifdef DEBUG
535  std::cout << " // Begin: list of generators per linear constraint //" << std::endl;
536  {
537  std::vector< std::set< boost::shared_ptr< Generator_Rn_SD > > >::iterator listOfGen1 = _listOfGeneratorsPerLinearConstraint.begin();
538  for (; listOfGen1 != _listOfGeneratorsPerLinearConstraint.end(); ++listOfGen1) {
539  std::cout << "LC" << listOfGen1 - _listOfGeneratorsPerLinearConstraint.begin() << ": ";
540  std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator itGN1 = listOfGen1->begin();
541  for (; itGN1 != listOfGen1->end(); ++itGN1) {
542  (*itGN1)->dump(std::cout);
543  std::cout << "; ";
544  }
545  std::cout << std::endl;
546  }
547  }
548  std::cout << " // End: list of generators per linear constraint //" << std::endl;
549 #endif
550  // Compute all neighbours
551  {
552  unsigned int count_IN = 0;
553  std::set< std::pair< unsigned int, unsigned int > > allNgb;
554  std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_1;
555  for (iteGN_1 = listOfGenSD.begin(); iteGN_1 != listOfGenSD.end(); ++iteGN_1, ++count_IN) {
556  unsigned int count_OUT = 0;
557  Neighbours_Rn newGeneratorsTool;
558  std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_2;
559  for (iteGN_2 = listOfGenSD.begin(); iteGN_2 != listOfGenSD.end(); ++iteGN_2, ++count_OUT) {
560  if (iteGN_1 != iteGN_2) {
561  std::vector< unsigned int > commonPFacets;
562  if (checkNeighbours(poly, *iteGN_1, *iteGN_2, commonPFacets) == true)
563  newGeneratorsTool.addGenerator(
564  commonPFacets,
565  count_IN,
566  count_OUT,
567  HalfSpace_Rn::hs_IN_OR_OUT);
568  }
569  }
570  std::set< boost::shared_ptr< Generator_Rn_SD > > setOfNGB;
571  boost::shared_ptr<Generator_Rn_SD> currentGeneratorIn;
572  for (newGeneratorsTool.begin(); newGeneratorsTool.end()!=true; newGeneratorsTool.next()) {
573  unsigned int firstNgbNb = newGeneratorsTool.currentGenInNumber();
574  unsigned int secondNgbNb = newGeneratorsTool.currentGenOutNumber();
575  // Enrich the list of neighbours
576  if ( allNgb.find( std::make_pair(firstNgbNb,secondNgbNb) ) == allNgb.end() &&
577  allNgb.find( std::make_pair(secondNgbNb,firstNgbNb) ) == allNgb.end() ) {
578  _allNeighbours.push_back( Segment_Rn(
579  listOfGenSD[newGeneratorsTool.currentGenInNumber()],
580  listOfGenSD[newGeneratorsTool.currentGenOutNumber()],
581  newGeneratorsTool.commonHalfSpaces()) );
582  allNgb.insert( std::make_pair(firstNgbNb, secondNgbNb) );
583  allNgb.insert( std::make_pair(secondNgbNb, firstNgbNb) );
584  }
585  }
586  //std::cout << "Insertion: ";
587  //(currentGeneratorIn)->save(std::cout);
588  //std::cout << std::endl;
589  }
590  }
591 
592  // We use this algorithm to compute vertices for a single polyhedron or to intersect two polyhedra.
593  // In the last case we truncate the polyhedron A generators net by the the polyhedron B half-spaces
594  // so we need to step forward in the half-space list where B constraints begin.
595  // We also use it in a generic way to intersect two polyhedral cones.
596  //constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS(poly->getListOfHalfSpaces());
597  generatorNumber = poly->numberOfGenerators();
598  {
599  _currentSegmentIndexAllNeighbours = 0;
600  std::vector< boost::shared_ptr<HalfSpace_Rn> >::iterator currentLC;
601  for (currentLC = _listOfLinearConstraints.begin()+truncationStep; currentLC != _listOfLinearConstraints.end(); ++currentLC) {
602  std::vector< boost::shared_ptr<HalfSpace_Rn> >::iterator bestLC = currentLC;
603  //std::cout << "Before" << std::endl;
604  if (_allNeighbours.size() == 0) {
605 #ifdef DEBUG
606  std::cout << "All segments all inactive (1)." << std::endl;
607 #endif
608  _isEmpty = true;
609  poly->reset();
610  return false;
611  }
612  findBestLinearConstraint(currentLC, bestLC);
613  //std::cout << "After" << std::endl;
614  if (currentLC != bestLC) {
615  // Swap the two linear constraints as they are not equal.
616  std::iter_swap(currentLC, bestLC);
617  std::vector< boost::shared_ptr<NormedHalfSpace_Rn> >::iterator bestNLC = _listOfNormedLinearConstraints.begin() + std::distance(_listOfLinearConstraints.begin(), bestLC);
618  std::vector< boost::shared_ptr<NormedHalfSpace_Rn> >::iterator currentNLC = _listOfNormedLinearConstraints.begin() + std::distance(_listOfLinearConstraints.begin(), currentLC);
619  std::iter_swap(currentNLC, bestNLC);
620  //boost::shared_ptr<HalfSpace_Rn> bHS = _listOfLinearConstraints[bestLC - _listOfLinearConstraints.begin()];
621  //_listOfLinearConstraints[currentLC - _listOfLinearConstraints.begin()] = _listOfLinearConstraints[bestLC - _listOfLinearConstraints.begin()];
622  //_listOfLinearConstraints[ bestLC - _listOfLinearConstraints.begin()] = bHS;
623  std::swap(_listOfLinearConstraintNorms[currentLC - _listOfLinearConstraints.begin()], _listOfLinearConstraintNorms[bestLC - _listOfLinearConstraints.begin()]);
624  }
625  //iteHS.setStep(truncationStep);
626  //{
627  //for (iteHS.begin(); iteHS.end()!=true; iteHS.next()) {
628  //const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace = iteHS.current();
629  const boost::shared_ptr<HalfSpace_Rn> currentHalfSpace = *currentLC;
630  boost::shared_ptr<Hyperplane_Rn> currentHyp = boost::dynamic_pointer_cast<Hyperplane_Rn>(currentHalfSpace);
631  bool isHyp = false;
632  if (currentHyp)
633  isHyp = true;
634  // When you add a linear constraint tell whether it is a half-space or a hyperplane.
635 #ifdef DEBUG
636  double halfSpaceNorm = std::inner_product(currentHalfSpace->begin(), currentHalfSpace->end(), currentHalfSpace->begin(), 0.);
637  halfSpaceNorm = sqrt(halfSpaceNorm);
638  std::cout << "Facet number " << halfspaceNumber << std::endl;
639  std::cout << "H" << halfspaceNumber << " = ";
640  currentHalfSpace->dump(std::cout);
641  std::cout << ", norm=" << halfSpaceNorm;
642  std::cout << std::endl;
643  // Show the content of _allNeighbours:
644  //{
645  //std::cout << " // Show the content of _allNeighbours: " << std::endl;
646  //std::vector< Segment_Rn >::iterator it_all1;
647  //for (it_all1 = _allNeighbours.begin(); it_all1 != _allNeighbours.end(); ++it_all1) {
648  //if (it_all1->getActive() == true) {
649  //std::cout << "( ";
650  //it_all1->getFirstGenerator()->dump(std::cout);
651  //std::cout << " , ";
652  //it_all1->getSecondGenerator()->dump(std::cout);
653  //std::cout << " )" << std::endl;
654  //}
655  //}
656  //}
657 #endif
658  // Clean the segment states.
659  {
660  std::vector< Segment_Rn >::iterator pairOfNeighbours3;
661  for (pairOfNeighbours3 = _allNeighbours.begin(); pairOfNeighbours3 != _allNeighbours.end(); ++pairOfNeighbours3) {
662  pairOfNeighbours3->getFirstGenerator()->setState(HalfSpace_Rn::hs_UNKNOWN);
663  pairOfNeighbours3->getSecondGenerator()->setState(HalfSpace_Rn::hs_UNKNOWN);
664  }
665  }
666 
667  // Fill the arrays of generators according to their position with respect to the current half-space.
668  std::vector< boost::shared_ptr<Generator_Rn_SD> > arrayOfGenOUT;
669  unsigned int nbGeneratorsIn = 0, nbGeneratorsOut = 0, nbGeneratorsOn = 0;
670  computeStatesOfAllGenerators(arrayOfGenOUT, currentHalfSpace, halfspaceNumber, nbGeneratorsIn, nbGeneratorsOut, nbGeneratorsOn);
671 
672 #ifdef DEBUG
673  {
674  std::cout << "Begin dump DS" << std::endl;
675  std::vector< Segment_Rn >::iterator pairOfNeighbours2print;
676  for (pairOfNeighbours2print = _allNeighbours.begin(); pairOfNeighbours2print != _allNeighbours.end(); ++pairOfNeighbours2print) {
677  std::cout << pairOfNeighbours2print - _allNeighbours.begin() << ". " ;
678  (pairOfNeighbours2print->getActive() == true) ? std::cout << "t, " : std::cout << "f, ";
679  std::cout << "NGB1= ";
680  pairOfNeighbours2print->getFirstGenerator()->dump(std::cout);
681  std::cout << ", NGB2= ";
682  pairOfNeighbours2print->getSecondGenerator()->dump(std::cout);
683  std::cout << std::endl;
684  }
685  std::cout << " End dump DS" << std::endl;
686  }
687 #endif
688  // Check whether the current constraint, i.e. a half-space or a hyperplane, excludes all Generators.
689  // In the case of a half-space, it used to be: if (GN_IN.empty() == true) {
690  if (currentHalfSpace->testEmptyness(nbGeneratorsIn, nbGeneratorsOut, nbGeneratorsOn)) {
691  // NO INTERSECTION
692 #ifdef DEBUG
693  std::cout << "Truncation is empty or flat." << std::endl;
694 #endif
695  poly->reset();
696  return false;
697  }
698 
699  // Check whether the current constraint is redundant. Only available for polytopes.
700  // For polyhedral cones, replace Rn::getDimension() by (Rn::getDimension()-1).
701  if (nbGeneratorsOut == 0) {// && VX_ON.size() < numberOfGeneratorsPerFacet()) {
702  // REDUNDANT HALFSPACE (be careful when modifying a list under iterator)
703  _listOfRedundantLinearConstraints.insert(halfspaceNumber);
704  //_listOfHalfSpaces.removeCurrentHalfSpace(iteHS.currentIteratorNumber());
705  //_listOfHalfSpaces.erase(HSPos);
706  // When chopping a polyhedral body with a hyperplane, if the set of generators OUT is empty,
707  // only keep the generators ON and ignore the generators IN.
708  // In the case of a half-space, keep things unchanged i.e. listOfGeneratorSD = GN_IN + GN_ON
709  currentHalfSpace->noGeneratorsOUT(listOfGenSD, GN_ON);
710 #ifdef DEBUG
711  std::cout << "VX_OUT= " << nbGeneratorsOut << std::endl;
712  std::cout << "VX_IN = " << nbGeneratorsIn << std::endl;
713  std::cout << "VX_ON = " << nbGeneratorsOn << std::endl;
714  std::cout << "redundantHS = " << _listOfRedundantLinearConstraints.size() << std::endl;
715  std::cout << "============================================" << std::endl;
716 #endif
717  }
718  else if (nbGeneratorsIn != 0) {
719  // Iterate through the list of all pairs of neighbours to do the hard job.
720  std::set< std::pair< unsigned int, unsigned int > > listOfAlreadyCreatedSegments;
721 
722  //std::cout << "LOOP0" << std::endl;
724  // LOOP0: select the generators with ON state //
726  {
727  std::set< boost::shared_ptr< Generator_Rn_SD > > setOfGeneratorsOn;
728  std::vector< Segment_Rn >::iterator pairOfNeighbours0;
729  for (pairOfNeighbours0 = _allNeighbours.begin(); pairOfNeighbours0 != _allNeighbours.end(); ++pairOfNeighbours0) {
730  if (pairOfNeighbours0->getActive()) {
731  HalfSpace_Rn::State firstState = pairOfNeighbours0->getFirstGenerator()->getState();
732  HalfSpace_Rn::State secondState = pairOfNeighbours0->getSecondGenerator()->getState();
733  if (firstState == HalfSpace_Rn::hs_ON) {
734  // ON - ?
735  // Test whether it was already inserted.
736  std::pair< std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator, bool > ret = setOfGeneratorsOn.insert(pairOfNeighbours0->getFirstGenerator());
737  if (ret.second == true) {
738  GN_ON.push_back(pairOfNeighbours0->getFirstGenerator());
739  pairOfNeighbours0->getFirstGenerator()->setFacet(halfspaceNumber);
740  _listOfGeneratorsPerLinearConstraint[halfspaceNumber].insert( pairOfNeighbours0->getFirstGenerator() );
741  // Case: ./trunk/test_politopix.cpp line 676 revision 409
742  // cout << "# intersections 11 & 12 #" << endl;
743  // BOOST_REQUIRE( politopixAPI::computeIntersection(arrayOfPolytopes[0], arrayOfPolytopes[10]) == TEST_OK );
744  Generator_Rn_SD::Status& currentStatus = pairOfNeighbours0->getFirstGenerator()->getStatus();
745  if (currentStatus == Generator_Rn_SD::CREATED)
746  currentStatus = Generator_Rn_SD::CREATED_AND_MODIFIED;
747  else if (currentStatus == Generator_Rn_SD::UNCHANGED)
748  currentStatus = Generator_Rn_SD::MODIFIED;
749  }
750  }
751  if (secondState == HalfSpace_Rn::hs_ON) {
752  // ? - ON
753  // Test whether it was already inserted.
754  std::pair< std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator, bool > ret = setOfGeneratorsOn.insert(pairOfNeighbours0->getSecondGenerator());
755  if (ret.second == true) {
756  GN_ON.push_back(pairOfNeighbours0->getSecondGenerator());
757  pairOfNeighbours0->getSecondGenerator()->setFacet(halfspaceNumber);
758  _listOfGeneratorsPerLinearConstraint[halfspaceNumber].insert( pairOfNeighbours0->getSecondGenerator() );
759  // Case: ./trunk/test_politopix.cpp line 676 revision 409
760  // cout << "# intersections 11 & 12 #" << endl;
761  // BOOST_REQUIRE( politopixAPI::computeIntersection(arrayOfPolytopes[0], arrayOfPolytopes[10]) == TEST_OK );
762  Generator_Rn_SD::Status& currentStatus = pairOfNeighbours0->getSecondGenerator()->getStatus();
763  if (currentStatus == Generator_Rn_SD::CREATED)
764  currentStatus = Generator_Rn_SD::CREATED_AND_MODIFIED;
765  else if (currentStatus == Generator_Rn_SD::UNCHANGED)
766  currentStatus = Generator_Rn_SD::MODIFIED;
767  }
768  }
769  if (firstState == HalfSpace_Rn::hs_ON && secondState == HalfSpace_Rn::hs_ON) {
770  unsigned int genON1 = pairOfNeighbours0->getFirstGenerator()->getGeneratorNumber();
771  unsigned int genON2 = pairOfNeighbours0->getSecondGenerator()->getGeneratorNumber();
772  listOfAlreadyCreatedSegments.insert( std::make_pair(genON1, genON2) );
773  listOfAlreadyCreatedSegments.insert( std::make_pair(genON2, genON1) );
774  }
775  } // if (pairOfNeighbours1->getActive()) {
776  }
777  }
778 
779  //std::cout << "LOOP1" << std::endl;
781  // LOOP1: create the new generators //
783  {
784  std::vector< unsigned int > commonFacets;
785  std::vector< Segment_Rn >::iterator pairOfNeighbours1;
786  for (pairOfNeighbours1 = _allNeighbours.begin(); pairOfNeighbours1 != _allNeighbours.end(); ++pairOfNeighbours1) {
787  if (pairOfNeighbours1->getActive()) {
788  //pairOfNeighbours1->getFirstGenerator()->dump(std::cout);
789  //pairOfNeighbours1->getSecondGenerator()->dump(std::cout);
790  //std::cout << std::endl;
791 
792  HalfSpace_Rn::State firstState = pairOfNeighbours1->getFirstGenerator()->getState();
793  HalfSpace_Rn::State secondState = pairOfNeighbours1->getSecondGenerator()->getState();
794  // Main job.
795  if ((firstState == HalfSpace_Rn::hs_IN && secondState == HalfSpace_Rn::hs_OUT)||
796  (firstState == HalfSpace_Rn::hs_OUT && secondState == HalfSpace_Rn::hs_IN)) {
797  boost::shared_ptr< Generator_Rn_SD > GenA, GenB;
798  if (firstState == HalfSpace_Rn::hs_IN && secondState == HalfSpace_Rn::hs_OUT) {
799  // IN - OUT
800  GenA = pairOfNeighbours1->getFirstGenerator();
801  GenB = pairOfNeighbours1->getSecondGenerator();
802  }
803  else {
804  // OUT - IN
805  GenB = pairOfNeighbours1->getFirstGenerator();
806  GenA = pairOfNeighbours1->getSecondGenerator();
807  }
808  // Create a new Generator as a combination of first and second.
809  double ay = GenB->getScalarProduct();
810  double az = GenA->getScalarProduct();
811  boost::shared_ptr< Generator_Rn_SD > newV;
812  // If an unused generator is referenced in the array listOfGeneratorsOUT2recycle, use it for newV.
813  if (!listOfGeneratorsOUT2recycle.empty()) {
814  newV = listOfGeneratorsOUT2recycle.back();
815  listOfGeneratorsOUT2recycle.pop_back();
816  newV->init(generatorNumber, Generator_Rn_SD::CREATED);
817  }
818  else
819  newV.reset(new Generator_Rn_SD(_RnDIM, generatorNumber, Generator_Rn_SD::CREATED));
820  // We move the generator on the frontier.
821  newV->setState(HalfSpace_Rn::hs_OUT_TO_ON);
822  ++generatorNumber;
823  poly->createTruncatedGenerator(GenB, GenA, newV, ay, az, currentHalfSpace->getConstant());
824  // Get facets computing the intersection on sorted lists of facets.
825  commonFacets.clear();
826  std::set_intersection(GenA->facetsBegin(), GenA->facetsEnd(), GenB->facetsBegin(), GenB->facetsEnd(), std::inserter(commonFacets, commonFacets.end()));
827  newV->setAllFacets(commonFacets);
828  // Test if this generator has already been created and inserted.
829  std::set< unsigned int > setOfFacets;
830  bool alreadyCreated = false;
831  // Check if this new vertex is coincident with an already created vertex with state ON.
832  boost::shared_ptr< Generator_Rn_SD > alreadyExistingGenerator;
833  {
834  std::vector< boost::shared_ptr< Generator_Rn_SD > >::const_iterator iteONGen;
835  for (iteONGen = GN_ON.begin(); iteONGen != GN_ON.end() && alreadyCreated == false; ++iteONGen) {
836  alreadyCreated = (*iteONGen)->isEqual2(newV, _RnDIM, TOL2);
837  if (alreadyCreated == true) {
838  // newV has no future so just transfer its facets to the current equal generator.
839  for (unsigned int i = 0; i < newV->numberOfFacets(); ++i)
840  _listOfGeneratorsPerLinearConstraint[ newV->getFacet(i) ].insert((*iteONGen));
841  // Make sure not to include twice currentHalfSpace.
842  alreadyExistingGenerator = (*iteONGen);
843  newV->exportFacets(setOfFacets);
844  (*iteONGen)->exportFacets(setOfFacets);
845  (*iteONGen)->importFacets(setOfFacets);
846  (*iteONGen)->orderFacets();
847  //(*iteONGen)->setState(HalfSpace_Rn::hs_OUT_TO_ON);
848  if (firstState == HalfSpace_Rn::hs_IN && secondState == HalfSpace_Rn::hs_OUT)
849  pairOfNeighbours1->setSecondGenerator(*iteONGen);
850  else
851  pairOfNeighbours1->setFirstGenerator(*iteONGen);
852  }
853  }
854  }
855  // Check if this new vertex is coincident with a freshly created new vertex.
856  {
857  std::vector< boost::shared_ptr< Generator_Rn_SD > >::const_iterator iteNewGen;
858  for (iteNewGen = GN_new.begin(); iteNewGen != GN_new.end() && alreadyCreated == false; ++iteNewGen) {
859  //(*iteNewGen)->dump(std::cout);
860  //std::cout << std::endl;
861  alreadyCreated = (*iteNewGen)->isEqual2(newV, _RnDIM, TOL2);
862  if (alreadyCreated == true) {
863  // newV has no future so just transfer its facets to the current equal generator.
864  for (unsigned int i = 0; i < newV->numberOfFacets(); ++i)
865  _listOfGeneratorsPerLinearConstraint[ newV->getFacet(i) ].insert((*iteNewGen));
866  newV->setFacet(halfspaceNumber);
867  newV->exportFacets(setOfFacets);
868  (*iteNewGen)->exportFacets(setOfFacets);
869  (*iteNewGen)->importFacets(setOfFacets);
870  (*iteNewGen)->orderFacets();
871  //(*iteNewGen)->setState(HalfSpace_Rn::hs_OUT_TO_ON);
872  alreadyExistingGenerator = (*iteNewGen);
873  if (firstState == HalfSpace_Rn::hs_IN && secondState == HalfSpace_Rn::hs_OUT)
874  pairOfNeighbours1->setSecondGenerator(*iteNewGen);
875  else
876  pairOfNeighbours1->setFirstGenerator(*iteNewGen);
877  // Let the LOOP2 deal with the rest of the data structures.
878  }
879  }
880  }
881 
882  if (alreadyCreated == false) {
883 #ifdef DEBUG
884  std::cout << "New Generator " << newV->getGeneratorNumber() << " = [";
885  newV->dump(std::cout);
886  std::cout << "] GeneratorIN ";
887  (firstState == HalfSpace_Rn::hs_IN) ? std::cout << pairOfNeighbours1->getFirstGenerator()->getGeneratorNumber() : std::cout << pairOfNeighbours1->getSecondGenerator()->getGeneratorNumber();
888  std::cout << " [";
889  (firstState == HalfSpace_Rn::hs_IN) ? pairOfNeighbours1->getFirstGenerator()->dump(std::cout) : pairOfNeighbours1->getSecondGenerator()->dump(std::cout);
890  std::cout << "], VX_IN_sp=" << az << " GeneratorOUT ";
891  (firstState == HalfSpace_Rn::hs_IN) ? std::cout << pairOfNeighbours1->getSecondGenerator()->getGeneratorNumber() : std::cout << pairOfNeighbours1->getFirstGenerator()->getGeneratorNumber();
892  std::cout << " [";
893  (firstState == HalfSpace_Rn::hs_IN) ? pairOfNeighbours1->getSecondGenerator()->dump(std::cout) : pairOfNeighbours1->getFirstGenerator()->dump(std::cout);
894  std::cout << "], VX_OUT_sp=" << ay;
895  std::cout << ", Common half-spaces: { ";
896  std::copy( commonFacets.begin(), commonFacets.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
897  std::cout << "}" << std::endl;
898 #endif
899  std::vector< unsigned int > potentialNewFacets1, newFacets1;
900  //std::set_symmetric_difference(pairOfNeighbours1->getFirstGenerator()->facetsBegin(), pairOfNeighbours1->getFirstGenerator()->facetsEnd(),
901  // pairOfNeighbours1->getSecondGenerator()->facetsBegin(), pairOfNeighbours1->getSecondGenerator()->facetsEnd(), std::inserter(potentialNewFacets1, potentialNewFacets1.end()));
902  std::set_difference(pairOfNeighbours1->getFirstGenerator()->facetsBegin(), pairOfNeighbours1->getFirstGenerator()->facetsEnd(),
903  pairOfNeighbours1->getSecondGenerator()->facetsBegin(), pairOfNeighbours1->getSecondGenerator()->facetsEnd(), std::inserter(potentialNewFacets1, potentialNewFacets1.end()));
904  // Project the new generator on all the facets of the first generator which are not common to both to detect if newV belong to them when we deal with fuzzy cases.
905  bool addFacets = false;
906  unsigned int size1 = potentialNewFacets1.size();
907  for (unsigned int i1 = 0; i1 < size1; ++i1) {
908  double scalarProduct = std::inner_product(newV->begin(), newV->end(), _listOfNormedLinearConstraints[ potentialNewFacets1[i1] ]->begin(), 0.);
909  double distanceToHyperplane = scalarProduct + _listOfNormedLinearConstraints[ potentialNewFacets1[i1] ]->getConstant();
910  //std::cout << " .distanceToHyperplane " << potentialNewFacets1[i1] << " = " << distanceToHyperplane << std::endl;
911  if (-_TOL < distanceToHyperplane && distanceToHyperplane < _TOL) {
912  newFacets1.push_back(potentialNewFacets1[i1]);
913  addFacets = true;
914  }
915  }
916  // Check the other end now.
917  // Project the new generator on all the facets of the second generator which are not common to both to detect if newV belongs to them when we deal with fuzzy cases.
918  std::vector< unsigned int > potentialNewFacets2, newFacets2;
919  std::set_difference(pairOfNeighbours1->getSecondGenerator()->facetsBegin(), pairOfNeighbours1->getSecondGenerator()->facetsEnd(),
920  pairOfNeighbours1->getFirstGenerator()->facetsBegin(), pairOfNeighbours1->getFirstGenerator()->facetsEnd(), std::inserter(potentialNewFacets2, potentialNewFacets2.end()));
921  unsigned int size2 = potentialNewFacets2.size();
922  for (unsigned int i2 = 0; i2 < size2; ++i2) {
923  double scalarProduct = std::inner_product(newV->begin(), newV->end(), _listOfNormedLinearConstraints[ potentialNewFacets2[i2] ]->begin(), 0.);
924  double distanceToHyperplane = scalarProduct + _listOfNormedLinearConstraints[ potentialNewFacets2[i2] ]->getConstant();
925  //std::cout << "..distanceToHyperplane " << potentialNewFacets2[i2] << " = " << distanceToHyperplane << std::endl;
926  if (-_TOL < distanceToHyperplane && distanceToHyperplane < _TOL) {
927  newFacets2.push_back(potentialNewFacets2[i2]);
928  addFacets = true;
929  }
930  }
931 #ifdef DEBUG
932  //std::cout << "potentialNewFacets2 { ";
933  //std::copy( potentialNewFacets2.begin(), potentialNewFacets2.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
934  //std::cout << "}, newFacets2 { ";
935  //std::copy( newFacets2.begin(), newFacets2.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
936  //std::cout << "}" << std::endl;
937  //std::cout << "potentialNewFacets1 { ";
938  //std::copy( potentialNewFacets1.begin(), potentialNewFacets1.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
939  //std::cout << "}, newFacets1 { ";
940  //std::copy( newFacets1.begin(), newFacets1.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
941  //std::cout << "}" << std::endl;
942 #endif
943  if (!newFacets1.empty() && newFacets1.size() == potentialNewFacets1.size() && !newFacets2.empty() && newFacets2.size() == potentialNewFacets2.size()) {
944  // Rare case but it can happen when the edge [pairOfNeighbours1->getFirstGenerator(), pairOfNeighbours1->getSecondGenerator()] is very short and when newV is insde all the facets of both ends.
945 #ifdef DEBUG
946  std::cout << "--Generator " << newV->getGeneratorNumber() << " contains all facets of generator " << pairOfNeighbours1->getFirstGenerator()->getGeneratorNumber();
947  std::cout << " and generator " << pairOfNeighbours1->getSecondGenerator()->getGeneratorNumber() << std::endl;
948 #endif
949  if (pairOfNeighbours1->getFirstGenerator()->getState() == HalfSpace_Rn::hs_IN) {
950  // Erase pairOfNeighbours1->getFirstGenerator() with newV
951  pairOfNeighbours1->getFirstGenerator()->setCoordinates(newV->vect());
952  Generator_Rn_SD::Status& currentStatus = pairOfNeighbours1->getFirstGenerator()->getStatus();
953  if (currentStatus == Generator_Rn_SD::CREATED)
954  currentStatus = Generator_Rn_SD::CREATED_AND_MODIFIED;
955  else if (currentStatus == Generator_Rn_SD::UNCHANGED)
956  currentStatus = Generator_Rn_SD::MODIFIED;
957  // Now recalculate the scalar product with the current half-space in case we need it later.
958  double scalarProduct = std::inner_product(pairOfNeighbours1->getFirstGenerator()->begin(), pairOfNeighbours1->getFirstGenerator()->end(), currentHalfSpace->begin(), 0.);
959  pairOfNeighbours1->getFirstGenerator()->setScalarProduct(scalarProduct);
960  // Now the segment is [NEW, OUT] so we can remove it.
961  pairOfNeighbours1->setInactive();
962  newV = pairOfNeighbours1->getFirstGenerator();
963  // Now transfer all the other facets from the second generator onto the first one.
964 #ifdef DEBUG
965  std::cout << "--Add to the generator " << newV->getGeneratorNumber() << " the following facets:";
966 #endif
967  for (unsigned int i2 = 0; i2 < newFacets2.size(); ++i2) {
968  newV->setFacet( newFacets2[i2] );
969  _listOfGeneratorsPerLinearConstraint[ newFacets2[i2] ].insert(newV);
970 #ifdef DEBUG
971  std::cout << " " << newFacets2[i2];
972 #endif
973  }
974 #ifdef DEBUG
975  std::cout << std::endl;
976 #endif
977  }
978  else if (pairOfNeighbours1->getSecondGenerator()->getState() == HalfSpace_Rn::hs_IN) {
979  // If all the facets of potentialNewFacets1 belong to pairOfNeighbours1->getFirstGenerator(), it means that (facets of pairOfNeighbours1->getFirstGenerator()) \subset (facets of newV)
980  // so erase pairOfNeighbours1->getFirstGenerator() if its state is HalfSpace_Rn::hs_IN to have it substituted by newV (if it's OUT it will be erased anyway).
981  // Erase the generator IN with the new generator as it contains all the facets of the generator IN.
982  pairOfNeighbours1->getSecondGenerator()->setCoordinates(newV->vect());
983  Generator_Rn_SD::Status& currentStatus = pairOfNeighbours1->getSecondGenerator()->getStatus();
984  if (currentStatus == Generator_Rn_SD::CREATED)
985  currentStatus = Generator_Rn_SD::CREATED_AND_MODIFIED;
986  else if (currentStatus == Generator_Rn_SD::UNCHANGED)
987  currentStatus = Generator_Rn_SD::MODIFIED;
988  // Now recalculate the scalar product with the current half-space in case we need it later.
989  double scalarProduct = std::inner_product(pairOfNeighbours1->getSecondGenerator()->begin(), pairOfNeighbours1->getSecondGenerator()->end(), currentHalfSpace->begin(), 0.);
990  pairOfNeighbours1->getSecondGenerator()->setScalarProduct(scalarProduct);
991  // Now the segment is [OUT, NEW] so we can remove it.
992  pairOfNeighbours1->setInactive();
993  newV = pairOfNeighbours1->getSecondGenerator();
994  // Transfer the remaining facets
995 #ifdef DEBUG
996  std::cout << "--Add to the generator " << newV->getGeneratorNumber() << " the following facets:";
997 #endif
998  for (unsigned int i1 = 0; i1 < newFacets1.size(); ++i1) {
999  newV->setFacet( newFacets1[i1] );
1000  _listOfGeneratorsPerLinearConstraint[ newFacets1[i1] ].insert(newV);
1001 #ifdef DEBUG
1002  std::cout << " " << newFacets1[i1];
1003 #endif
1004  }
1005 #ifdef DEBUG
1006  std::cout << std::endl;
1007 #endif
1008  }
1009  }
1010  else if (!newFacets1.empty() && newFacets1.size() == potentialNewFacets1.size()) {
1011  if (pairOfNeighbours1->getFirstGenerator()->getState() == HalfSpace_Rn::hs_IN) {
1012 #ifdef DEBUG
1013  std::cout << " .Change generator IN " << pairOfNeighbours1->getFirstGenerator()->getGeneratorNumber() << " coordinates with generator " << newV->getGeneratorNumber() << std::endl;
1014 #endif
1015  // Erase pairOfNeighbours1->getFirstGenerator() with newV
1016  pairOfNeighbours1->getFirstGenerator()->setCoordinates(newV->vect());
1017  Generator_Rn_SD::Status& currentStatus = pairOfNeighbours1->getFirstGenerator()->getStatus();
1018  if (currentStatus == Generator_Rn_SD::CREATED)
1019  currentStatus = Generator_Rn_SD::CREATED_AND_MODIFIED;
1020  else if (currentStatus == Generator_Rn_SD::UNCHANGED)
1021  currentStatus = Generator_Rn_SD::MODIFIED;
1022  // Now recalculate the scalar product with the current half-space in case we need it later.
1023  double scalarProduct = std::inner_product(pairOfNeighbours1->getFirstGenerator()->begin(), pairOfNeighbours1->getFirstGenerator()->end(), currentHalfSpace->begin(), 0.);
1024  pairOfNeighbours1->getFirstGenerator()->setScalarProduct(scalarProduct);
1025  // Now the segment is [NEW, OUT] so we can remove it.
1026  pairOfNeighbours1->setInactive();
1027  newV = pairOfNeighbours1->getFirstGenerator();
1028  }
1029  else if (pairOfNeighbours1->getFirstGenerator()->getState() == HalfSpace_Rn::hs_OUT) {
1030  // Now the segment is [NEW, IN].
1031  // Just transfer to newV the facets of pairOfNeighbours1->getFirstGenerator(), the OUT generator will be discarded later on.
1032 #ifdef DEBUG
1033  std::cout << " .Change generator OUT " << pairOfNeighbours1->getFirstGenerator()->getGeneratorNumber() << " coordinates with generator " << newV->getGeneratorNumber() << std::endl;
1034  std::cout << " .Add to this generator " << newV->getGeneratorNumber() << " the following facets:";
1035 #endif
1036  for (unsigned int i1 = 0; i1 < newFacets1.size(); ++i1) {
1037  newV->setFacet( newFacets1[i1] );
1038  _listOfGeneratorsPerLinearConstraint[ newFacets1[i1] ].insert(newV);
1039 #ifdef DEBUG
1040  std::cout << " " << newFacets1[i1];
1041 #endif
1042  }
1043 #ifdef DEBUG
1044  std::cout << std::endl;
1045 #endif
1046  }
1047  }
1048  else if (!newFacets2.empty() && newFacets2.size() == potentialNewFacets2.size()) {
1049  if (pairOfNeighbours1->getSecondGenerator()->getState() == HalfSpace_Rn::hs_IN) {
1050 #ifdef DEBUG
1051  std::cout << "..Change generator IN " << pairOfNeighbours1->getSecondGenerator()->getGeneratorNumber() << " coordinates with generator " << newV->getGeneratorNumber() << std::endl;
1052 #endif
1053  // Erase the generator IN with the new generator as it contains all the facets of the generator IN.
1054  pairOfNeighbours1->getSecondGenerator()->setCoordinates(newV->vect());
1055  Generator_Rn_SD::Status& currentStatus = pairOfNeighbours1->getSecondGenerator()->getStatus();
1056  if (currentStatus == Generator_Rn_SD::CREATED)
1057  currentStatus = Generator_Rn_SD::CREATED_AND_MODIFIED;
1058  else if (currentStatus == Generator_Rn_SD::UNCHANGED)
1059  currentStatus = Generator_Rn_SD::MODIFIED;
1060  // Now recalculate the scalar product with the current half-space in case we need it later.
1061  double scalarProduct = std::inner_product(pairOfNeighbours1->getSecondGenerator()->begin(), pairOfNeighbours1->getSecondGenerator()->end(), currentHalfSpace->begin(), 0.);
1062  pairOfNeighbours1->getSecondGenerator()->setScalarProduct(scalarProduct);
1063  // Now the segment is [OUT, NEW] so we can remove it.
1064  pairOfNeighbours1->setInactive();
1065  newV = pairOfNeighbours1->getSecondGenerator();
1066  }
1067  else if (pairOfNeighbours1->getSecondGenerator()->getState() == HalfSpace_Rn::hs_OUT) {
1068  // Now the segment is [IN, NEW].
1069  // Just transfer to newV the facets of pairOfNeighbours1->getSecondGenerator(), the OUT generator will be discarded later on.
1070 #ifdef DEBUG
1071  std::cout << "..Change generator OUT " << pairOfNeighbours1->getSecondGenerator()->getGeneratorNumber() << " coordinates with generator " << newV->getGeneratorNumber() << std::endl;
1072  std::cout << "..Add to this generator " << newV->getGeneratorNumber() << " the following facets:";
1073 #endif
1074  for (unsigned int i2 = 0; i2 < newFacets2.size(); ++i2) {
1075  newV->setFacet( newFacets2[i2] );
1076  _listOfGeneratorsPerLinearConstraint[ newFacets2[i2] ].insert(newV);
1077 #ifdef DEBUG
1078  std::cout << " " << newFacets2[i2];
1079 #endif
1080  }
1081 #ifdef DEBUG
1082  std::cout << std::endl;
1083 #endif
1084  }
1085  }
1086  // Make the connections with the DS.
1087  {
1088  if (!newFacets1.empty() && newFacets1.size() < size1) {
1089  // Just insert all the new facets in the data structure.
1090  size1 = newFacets1.size();
1091 #ifdef DEBUG
1092  std::cout << " .Add to the generator " << newV->getGeneratorNumber() << " the following facets:";
1093 #endif
1094  for (unsigned int i1 = 0; i1 < size1; ++i1) {
1095  newV->setFacet( newFacets1[i1] );
1096  _listOfGeneratorsPerLinearConstraint[ newFacets1[i1] ].insert(newV);
1097 #ifdef DEBUG
1098  std::cout << " " << newFacets1[i1];
1099 #endif
1100  }
1101 #ifdef DEBUG
1102  std::cout << std::endl;
1103 #endif
1104  }
1105  if (!newFacets2.empty() && newFacets2.size() < size2) {
1106  size2 = newFacets2.size();
1107 #ifdef DEBUG
1108  std::cout << "..Add to the generator " << newV->getGeneratorNumber() << " the following facets:";
1109 #endif
1110  for (unsigned int i2 = 0; i2 < size2; ++i2) {
1111  newV->setFacet( newFacets2[i2] );
1112  _listOfGeneratorsPerLinearConstraint[ newFacets2[i2] ].insert(newV);
1113 #ifdef DEBUG
1114  std::cout << " " << newFacets2[i2];
1115 #endif
1116  }
1117 #ifdef DEBUG
1118  std::cout << std::endl;
1119 #endif
1120  }
1121  }
1122  for (unsigned int i = 0; i < newV->numberOfFacets(); ++i)
1123  _listOfGeneratorsPerLinearConstraint[ newV->getFacet(i) ].insert(newV);
1124  // If fuzzy facets were added, we need to reorder the list to make sure it is compliant with std::set_intersection() in the future.
1125  if (addFacets == true)
1126  newV->orderFacets();
1127  // Add the current half-space at this step if we don't fuse it with an ON generator.
1128  newV->setFacet(halfspaceNumber);
1129  _listOfGeneratorsPerLinearConstraint[halfspaceNumber].insert(newV);
1130  //for (unsigned int i = 0; i < newV->numberOfFacets(); ++i)
1131  //_listOfGeneratorsPerLinearConstraint[ newV->getFacet(i) ].insert(newV);
1132  //pairOfNeighbours1->dump(std::cout);
1133 
1134  // Erase the OUT generator with the freshly computed one.
1135  boost::shared_ptr< Generator_Rn_SD > ptr = newV;
1136  if (pairOfNeighbours1->getActive() == true) {
1137  if (firstState == HalfSpace_Rn::hs_IN && secondState == HalfSpace_Rn::hs_OUT) {
1138  // Now that we have destroyed the connections between the OUT generator and its facets, we can discard it.
1139  pairOfNeighbours1->setSecondGenerator(newV);
1140  ptr = pairOfNeighbours1->getSecondGenerator();
1141  }
1142  else {
1143  // Now that we have destroyed the connections between the OUT generator and its facets, we can discard it.
1144  pairOfNeighbours1->setFirstGenerator(newV);
1145  ptr = pairOfNeighbours1->getFirstGenerator();
1146  }
1147  }
1148  //ptr->dump(std::cout);
1149  //std::cout << std::endl;
1150  //pairOfNeighbours1->dump(std::cout);
1151 
1152  GN_new.push_back(ptr);
1153  //pairOfNeighbours1->dump(std::cout);
1154  //newV->dump(std::cout);
1155  //std::cout << std::endl;
1156  //pairOfNeighbours1->dump(std::cout);
1157  //std::cout << std::endl;
1158  } // if (alreadyCreated == false) {
1159 #ifdef DEBUG
1160  else {
1161  std::cout << "No new Generator " << newV->getGeneratorNumber() << " = [";
1162  newV->dump(std::cout);
1163  std::cout << "] GeneratorIN ";
1164  (firstState == HalfSpace_Rn::hs_IN) ? std::cout << pairOfNeighbours1->getFirstGenerator()->getGeneratorNumber() : std::cout << pairOfNeighbours1->getSecondGenerator()->getGeneratorNumber();
1165  std::cout << " [";
1166  (firstState == HalfSpace_Rn::hs_IN) ? pairOfNeighbours1->getFirstGenerator()->dump(std::cout) : pairOfNeighbours1->getSecondGenerator()->dump(std::cout);
1167  std::cout << "], VX_IN_sp=" << az << " GeneratorOUT ";
1168  (firstState == HalfSpace_Rn::hs_IN) ? std::cout << pairOfNeighbours1->getSecondGenerator()->getGeneratorNumber() : std::cout << pairOfNeighbours1->getFirstGenerator()->getGeneratorNumber();
1169  std::cout << " [";
1170  (firstState == HalfSpace_Rn::hs_IN) ? pairOfNeighbours1->getSecondGenerator()->dump(std::cout) : pairOfNeighbours1->getFirstGenerator()->dump(std::cout);
1171  std::cout << "], VX_OUT_sp= " << ay << ". This generator is equal to this one: ";
1172  alreadyExistingGenerator->dump(std::cout);
1173  std::cout << std::endl;
1174  }
1175 #endif
1176  } // if ((firstState == hs_IN && secondState == hs_OUT)|| (firstState == hs_OUT && secondState == hs_IN)) {
1177  } // if (pairOfNeighbours1->getActive()) {
1178 
1179  } // for (pairOfNeighbours1 = _allNeighbours.begin(); pairOfNeighbours1 != _allNeighbours.end(); ++pairOfNeighbours1) {
1180  }
1181 #ifdef DEBUG
1182  //std::cout << "Before LOOP2" << std::endl;
1183  //{
1184  //std::cout << "==========================================================================" << std::endl;
1185  //std::vector< Segment_Rn >::iterator it_all1;
1186  //for (it_all1 = _allNeighbours.begin(); it_all1 != _allNeighbours.end(); ++it_all1) {
1187  //if (it_all1->getActive() == true) {
1188  //std::cout << it_all1 - _allNeighbours.begin() << ". ";
1189  //it_all1->dump(std::cout);
1190  //std::cout << std::endl;
1191  //}
1192  //}
1193  //std::cout << "==========================================================================" << std::endl;
1194  //}
1195 #endif
1196 
1197  //std::cout << "LOOP2" << std::endl;
1199  // LOOP2: reorganize the list of all segments //
1201  {
1202  std::vector< Segment_Rn >::iterator pairOfNeighbours2;
1203  for (pairOfNeighbours2 = _allNeighbours.begin(); pairOfNeighbours2 != _allNeighbours.end(); ++pairOfNeighbours2) {
1204  if (pairOfNeighbours2->getActive()) {
1205  HalfSpace_Rn::State firstState = pairOfNeighbours2->getFirstGenerator()->getState();
1206  HalfSpace_Rn::State secondState = pairOfNeighbours2->getSecondGenerator()->getState();
1207  if (firstState == HalfSpace_Rn::hs_OUT && secondState == HalfSpace_Rn::hs_OUT) {
1208  // OUT - OUT
1209  listOfPairsToRemove.push_back( pairOfNeighbours2 - _allNeighbours.begin() );
1210  _allNeighbours[pairOfNeighbours2 - _allNeighbours.begin()].setInactive();
1211  }
1212  else if (firstState == HalfSpace_Rn::hs_OUT && (secondState == HalfSpace_Rn::hs_ON || secondState == HalfSpace_Rn::hs_OUT_TO_ON)) {
1213  // OUT - ON
1214  // OUT - OUT_TO_ON
1215  listOfPairsToRemove.push_back( pairOfNeighbours2 - _allNeighbours.begin() );
1216  _allNeighbours[pairOfNeighbours2 - _allNeighbours.begin()].setInactive();
1217  }
1218  else if ((firstState == HalfSpace_Rn::hs_ON || firstState == HalfSpace_Rn::hs_OUT_TO_ON) && secondState == HalfSpace_Rn::hs_OUT) {
1219  // ON - OUT
1220  // OUT_TO_ON - OUT
1221  listOfPairsToRemove.push_back( pairOfNeighbours2 - _allNeighbours.begin() );
1222  _allNeighbours[pairOfNeighbours2 - _allNeighbours.begin()].setInactive();
1223  }
1224  else if (firstState == HalfSpace_Rn::hs_ON && secondState == HalfSpace_Rn::hs_IN) {
1225  // ON - IN
1226  // Job already done in LOOP0 when we selected the generators with ON state.
1227  }
1228  else if (firstState == HalfSpace_Rn::hs_IN && secondState == HalfSpace_Rn::hs_ON) {
1229  // IN - ON
1230  // Job already done in LOOP0 when we selected the generators with ON state.
1231  }
1232  else if (firstState == HalfSpace_Rn::hs_ON && secondState == HalfSpace_Rn::hs_ON) {
1233  // ON - ON
1234  // Job already done in LOOP0 when we selected the generators with ON state.
1235  }
1236  else {
1237  // IN - IN
1238  // OUT_TO_ON - ?
1239  // ? - OUT_TO_ON
1240  }
1241  } // if (pairOfNeighbours2->getActive()) {
1242  }
1243  }
1244  {
1245  // Make sure the linear constraints do not reference anymore an OUT generator.
1246  std::vector< boost::shared_ptr<Generator_Rn_SD> >::iterator iteOUT;
1247  for (iteOUT=arrayOfGenOUT.begin(); iteOUT!=arrayOfGenOUT.end(); ++iteOUT) {
1248  std::vector<unsigned int>::const_iterator itFacet;
1249  for (itFacet = (*iteOUT)->facetsBegin(); itFacet != (*iteOUT)->facetsEnd(); ++itFacet)
1250  _listOfGeneratorsPerLinearConstraint[*itFacet].erase(*iteOUT);
1251  }
1252  }
1253  listOfGeneratorsOUT2recycle.insert(listOfGeneratorsOUT2recycle.end(), arrayOfGenOUT.begin(), arrayOfGenOUT.end());
1254  arrayOfGenOUT.clear();
1256  // Make the new connections for the new generators //
1258  // For each new generator build the list of its neighbours on the basis of an algorithm that removes a potential neighbour G if its list
1259  // of linear constraints L is such as \f[ L \subset L' \f] where L' is the list of linear constraints of G'. In this process we don't
1260  // include the already known neighbour of G with state IN as it would be useless (its list does not contain halfspaceNumber).
1261  {
1262  // Segments [NEW,NEW] and [NEW,ON]
1263  {
1264  unsigned int count_new1 = 0, count_new2 = 0;
1265  std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_new1, iteGN_new2, iteGN_ON1;
1266  for (iteGN_new1 = GN_new.begin(); iteGN_new1 != GN_new.end(); ++iteGN_new1, ++count_new1) {
1267  count_new2 = 0;
1268  //(*iteGN_new1)->dump(std::cout);
1269  //std::cout << std::endl;
1270  Neighbours_Rn newGeneratorsTool;
1271  {
1272  for (iteGN_new2 = GN_new.begin(); iteGN_new2 != GN_new.end(); ++iteGN_new2, ++count_new2) {
1273  std::vector< unsigned int > commonPFacets;
1274  //(*iteGN_new2)->dump(std::cout);
1275  //std::cout << std::endl;
1276  if (iteGN_new1 != iteGN_new2 && checkNeighbours(poly, *iteGN_new1, *iteGN_new2, commonPFacets) == true) {
1277  //if (_RnDIM > 3)
1278  newGeneratorsTool.addGenerator(commonPFacets,
1279  count_new1,
1280  count_new2,
1281  HalfSpace_Rn::hs_IN_OR_OUT);
1282  //else
1283  //newGeneratorsTool.addGeneratorWithoutCheck(commonPFacets,
1284  //count_ON1,
1285  //count_ON2,
1286  //HalfSpace_Rn::hs_ON);
1287  }
1288  }
1289  }
1290  //std::cout << " GN_new done " << std::endl;
1291  {
1292  // Do not forget the generators with state ON.
1293  for (iteGN_ON1 = GN_ON.begin(); iteGN_ON1 != GN_ON.end(); ++iteGN_ON1, ++count_new2) {
1294  std::vector< unsigned int > commonPFacets;
1295  if (checkNeighbours(poly, *iteGN_new1, *iteGN_ON1, commonPFacets) == true) {
1296  //if (_RnDIM > 3)
1297  newGeneratorsTool.addGenerator(commonPFacets,
1298  count_new1,
1299  count_new2,
1300  HalfSpace_Rn::hs_IN_OR_OUT);
1301  //else
1302  //newGeneratorsTool.addGeneratorWithoutCheck(
1303  //commonPFacets,
1304  //count_ON,
1305  //count_OUT,
1306  //HalfSpace_Rn::hs_ON);
1307  }
1308  }
1309  }
1310 #ifdef DEBUG
1311  //std::cout << " //Begin dump DS //" << std::endl;
1312  //newGeneratorsTool.dump(std::cout);
1313  //std::cout << " // End dump DS //" << std::endl;
1314 #endif
1315 #ifdef DEBUG
1316  std::cout << " // Get the new segments //" << std::endl;
1317 #endif
1318  for (newGeneratorsTool.begin(); newGeneratorsTool.end() != true; newGeneratorsTool.next()) {
1319  unsigned int gen1 = newGeneratorsTool.currentGenInNumber();
1320  unsigned int gen2 = newGeneratorsTool.currentGenOutNumber();
1321  unsigned int genInNb = GN_new[gen1]->getGeneratorNumber();
1322  unsigned int genOutNb = (gen2 < GN_new.size()) ? GN_new[gen2]->getGeneratorNumber() : GN_ON[gen2 - GN_new.size()]->getGeneratorNumber();
1323  if (listOfAlreadyCreatedSegments.find( std::make_pair(genInNb, genOutNb) ) == listOfAlreadyCreatedSegments.end() &&
1324  listOfAlreadyCreatedSegments.find( std::make_pair(genOutNb, genInNb) ) == listOfAlreadyCreatedSegments.end()) {
1325 #ifdef DEBUG
1326  std::cout << "NGB1 = ";
1327  GN_new[gen1]->dump(std::cout);
1328  std::cout << ", NGB2 = ";
1329  (gen2 < GN_new.size()) ? GN_new[gen2]->dump(std::cout) : GN_ON[gen2 - GN_new.size()]->dump(std::cout);
1330  std::cout << std::endl;
1331 #endif
1332  Segment_Rn newSegment(newGeneratorsTool.commonHalfSpaces());
1333  newSegment.setFirstGenerator( GN_new[gen1] );
1334  // Test in which array we get the second generator.
1335  //std::cout << " // New segments 0 // " << std::endl;
1336  //std::cout << gen1 << " " << gen2 << " " << genInNb << " " << genOutNb << std::endl;
1337  if (gen2 < GN_new.size())
1338  newSegment.setSecondGenerator( GN_new[gen2] );
1339  else
1340  newSegment.setSecondGenerator( GN_ON[gen2 - GN_new.size()] );
1341  // Try to find an empty place to locate the new segment.
1342  //std::cout << " // New segments 1 //" << std::endl;
1343  //newSegment.dump(std::cout);
1344  //std::cout << std::endl;
1345  //std::cout << "_allNeighbours.size() = " << _allNeighbours.size() << ", listOfPairsToRemove.size() = " << listOfPairsToRemove.size() << std::endl;
1346  if (listOfPairsToRemove.empty() == true) {
1347  //std::cout << " // New segments 1.1 a //" << std::endl;
1348  _allNeighbours.push_back(newSegment);
1349  }
1350  else {
1351  //std::cout << " // New segments 1.1 b //" << std::endl;
1352  _allNeighbours[ listOfPairsToRemove[listOfPairsToRemove.size()-1] ] = newSegment;
1353  listOfPairsToRemove.pop_back();
1354  }
1355  //if (listOfPairsToRemove.empty() != false) {
1356  //std::cout << " // New segments 1.1 a //" << std::endl;
1357  //_allNeighbours[listOfPairsToRemove.size()-1] = newSegment;
1358  //listOfPairsToRemove.pop_back();
1359  //}
1360  //else {
1361  //std::cout << " // New segments 1.1 b //" << std::endl;
1362  //_allNeighbours.push_back(newSegment);
1363  //}
1364  //std::cout << " // New segments 1.2 //" << std::endl;
1365  // Mark this segment as done to make sure that when [gen1, gen2] is stored we don't do the same with [gen2, gen1].
1366  listOfAlreadyCreatedSegments.insert( std::make_pair(genInNb, genOutNb) );
1367  listOfAlreadyCreatedSegments.insert( std::make_pair(genOutNb, genInNb) );
1368  //std::cout << " // New segments 2 //" << std::endl;
1369  }
1370  }
1371  }
1372  } // Segments [NEW,NEW] and [NEW,ON]
1373  // Segments [ON,ON] and [ON,NEW]
1374  {
1375  //std::cout << "Size of GN_ON: " << GN_ON.size() << std::endl;
1376  unsigned int count_on1 = 0, count_on2 = 0, count_new1 = 0;
1377  std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_new1, iteGN_ON1, iteGN_ON2;
1378  for (iteGN_ON1 = GN_ON.begin(); iteGN_ON1 != GN_ON.end(); ++iteGN_ON1, ++count_on1) {
1379  count_on2 = 0;
1380  Neighbours_Rn newGeneratorsTool;
1381  {
1382  for (iteGN_ON2 = GN_ON.begin(); iteGN_ON2 != GN_ON.end(); ++iteGN_ON2, ++count_on2) {
1383  std::vector< unsigned int > commonPFacets;
1384  bool toTest = checkNeighbours(poly, *iteGN_ON1, *iteGN_ON2, commonPFacets);
1385  if (iteGN_ON1 != iteGN_ON2 && toTest == true) {
1386 #ifdef DEBUG
1387  //std::cout << "test = " << toTest << ", ngbC = " << poly->neigbourhoodCondition() << ", ";
1388  //std::cout << "commonPFacets.size() = " << commonPFacets.size() << ", ";
1389  //(*iteGN_ON1)->dump(std::cout); std::cout << " <=> ";
1390  //(*iteGN_ON2)->dump(std::cout); std::cout << " -> ";
1391  //std::copy(commonPFacets.begin(), commonPFacets.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
1392  //std::cout << std::endl;
1393 #endif
1394  //if (_RnDIM > 3)
1395  newGeneratorsTool.addGenerator(commonPFacets,
1396  count_on1,
1397  count_on2,
1398  HalfSpace_Rn::hs_IN_OR_OUT);
1399  //else
1400  //newGeneratorsTool.addGeneratorWithoutCheck(commonPFacets,
1401  //count_ON1,
1402  //count_ON2,
1403  //HalfSpace_Rn::hs_ON);
1404  }
1405  }
1406  }
1407  {
1408  // Do not forget the generators with state NEW.
1409  for (iteGN_new1 = GN_new.begin(); iteGN_new1 != GN_new.end(); ++iteGN_new1, ++count_on2) {
1410  std::vector< unsigned int > commonPFacets;
1411  if (checkNeighbours(poly, *iteGN_ON1, *iteGN_new1, commonPFacets) == true) {
1412  //if (_RnDIM > 3)
1413  newGeneratorsTool.addGenerator(commonPFacets,
1414  count_on1,
1415  count_on2,
1416  HalfSpace_Rn::hs_IN_OR_OUT);
1417  //else
1418  //newGeneratorsTool.addGeneratorWithoutCheck(
1419  //commonPFacets,
1420  //count_ON,
1421  //count_OUT,
1422  //HalfSpace_Rn::hs_ON);
1423  }
1424  }
1425  }
1426 #ifdef DEBUG
1427  std::cout << " //Begin dump DS_ //" << std::endl;
1428  newGeneratorsTool.dump(std::cout);
1429  std::cout << " // End dump DS_ //" << std::endl;
1430 #endif
1431  for (newGeneratorsTool.begin(); newGeneratorsTool.end() != true; newGeneratorsTool.next()) {
1432  unsigned int gen1 = newGeneratorsTool.currentGenInNumber();
1433  unsigned int gen2 = newGeneratorsTool.currentGenOutNumber();
1434  unsigned int genInNb = GN_ON[gen1]->getGeneratorNumber();
1435  unsigned int genOutNb = (gen2 < GN_ON.size()) ? GN_ON[gen2]->getGeneratorNumber() : GN_new[gen2 - GN_ON.size()]->getGeneratorNumber();
1436  if (listOfAlreadyCreatedSegments.find( std::make_pair(genInNb, genOutNb) ) == listOfAlreadyCreatedSegments.end() &&
1437  listOfAlreadyCreatedSegments.find( std::make_pair(genOutNb, genInNb) ) == listOfAlreadyCreatedSegments.end()) {
1438 #ifdef DEBUG
1439  std::cout << "NGB1_= ";
1440  GN_ON[gen1]->dump(std::cout);
1441  std::cout << ", NGB2_= ";
1442  (gen2 < GN_ON.size()) ? GN_ON[gen2]->dump(std::cout) : GN_new[gen2 - GN_ON.size()]->dump(std::cout);
1443  std::cout << std::endl;
1444 #endif
1445  Segment_Rn newSegment(newGeneratorsTool.commonHalfSpaces());
1446  newSegment.setFirstGenerator( GN_ON[gen1] );
1447  // Test in which array we get the second generator.
1448  if (gen2 < GN_ON.size())
1449  newSegment.setSecondGenerator( GN_ON[gen2] );
1450  else
1451  newSegment.setSecondGenerator( GN_new[gen2 - GN_ON.size()] );
1452  // Try to find an empty place to locate the new segment.
1453  if (listOfPairsToRemove.empty() == true) {
1454  //std::cout << " // New segments 1.1 a //" << std::endl;
1455  _allNeighbours.push_back(newSegment);
1456  }
1457  else {
1458  //std::cout << " // New segments 1.1 b //" << std::endl;
1459  _allNeighbours[ listOfPairsToRemove[listOfPairsToRemove.size()-1] ] = newSegment;
1460  listOfPairsToRemove.pop_back();
1461  }
1462  //if (listOfPairsToRemove.empty() != false) {
1463  //_allNeighbours[listOfPairsToRemove.size()-1] = newSegment;
1464  //listOfPairsToRemove.pop_back();
1465  //}
1466  //else
1467  //_allNeighbours.push_back(newSegment);
1468  // Mark this segment as done to make sure that when [gen1, gen2] is stored we don't do the same with [gen2, gen1].
1469  listOfAlreadyCreatedSegments.insert( std::make_pair(genInNb, genOutNb) );
1470  listOfAlreadyCreatedSegments.insert( std::make_pair(genOutNb, genInNb) );
1471  }
1472  }
1473  }
1474  } // Segments [ON,NEW] and [ON,ON]
1475  }
1476 
1477 #ifdef DEBUG
1478  {
1479  std::cout << "Begin dump DS" << std::endl;
1480  std::vector< Segment_Rn >::iterator pairOfNeighbours2print;
1481  for (pairOfNeighbours2print = _allNeighbours.begin(); pairOfNeighbours2print != _allNeighbours.end(); ++pairOfNeighbours2print) {
1482  std::cout << pairOfNeighbours2print - _allNeighbours.begin() << ". " ;
1483  (pairOfNeighbours2print->getActive() == true) ? std::cout << "t: " : std::cout << "f: ";
1484  std::cout << "NGB1= ";
1485  pairOfNeighbours2print->getFirstGenerator()->dump(std::cout);
1486  std::cout << ", NGB2= ";
1487  pairOfNeighbours2print->getSecondGenerator()->dump(std::cout);
1488  std::cout << std::endl;
1489  }
1490  std::cout << " End dump DS" << std::endl;
1491  }
1492 #endif
1493 
1495  // Deal with redundancy //
1497 #ifdef DEBUG
1498  std::cout << " // Begin redundancy processing //" << std::endl;
1499  {
1500  std::vector< std::set< boost::shared_ptr< Generator_Rn_SD > > >::iterator listOfGen1 = _listOfGeneratorsPerLinearConstraint.begin();
1501  for (; listOfGen1 != _listOfGeneratorsPerLinearConstraint.end(); ++listOfGen1) {
1502  std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator itGN1 = listOfGen1->begin();
1503  std::vector< unsigned int > gnNbinLC;
1504  for (; itGN1 != listOfGen1->end(); ++itGN1) {
1505  //(*itGN1)->dump(std::cout);
1506  //std::cout << "; ";
1507  gnNbinLC.push_back( (*itGN1)->getGeneratorNumber() );
1508  //std::cout << (*itGN1)->getGeneratorNumber() << " ";
1509  }
1510  if (!gnNbinLC.empty()) {
1511  std::sort(gnNbinLC.begin(), gnNbinLC.end());
1512  std::cout << "LC" << listOfGen1 - _listOfGeneratorsPerLinearConstraint.begin() << ": ";
1513  std::copy(gnNbinLC.begin(), gnNbinLC.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
1514  if (listOfGen1 - _listOfGeneratorsPerLinearConstraint.begin() > truncationStep && listOfGen1 - _listOfGeneratorsPerLinearConstraint.begin() <= halfspaceNumber)
1515  _listOfLinearConstraints[ listOfGen1 - _listOfGeneratorsPerLinearConstraint.begin() ]->dump(std::cout);
1516  std::cout << std::endl;
1517  }
1518  }
1519  }
1520  std::cout << " // End redundancy processing //" << std::endl;
1521 #endif
1522  // Declare a linear constraint redundant if it does not contain the sufficent number of generators.
1523  {
1524  unsigned int countLC = 0;
1525  std::vector< std::set< boost::shared_ptr< Generator_Rn_SD > > >::iterator listOfGen = _listOfGeneratorsPerLinearConstraint.begin();
1526  for (; listOfGen != _listOfGeneratorsPerLinearConstraint.end() && countLC != halfspaceNumber; ++listOfGen, ++countLC) {
1527 #ifdef DEBUG
1528  std::cout << listOfGen - _listOfGeneratorsPerLinearConstraint.begin() << ". ";
1529  {
1530  std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator itGN0 = listOfGen->begin();
1531  if (itGN0 != listOfGen->end()) {
1532  for (; itGN0 != listOfGen->end(); ++itGN0) {
1533  (*itGN0)->dump(std::cout);
1534  std::cout << " ";
1535  }
1536  std::cout << std::endl;
1537  }
1538  }
1539 #endif
1540  if (listOfGen->size() < poly->numberOfGeneratorsPerFacet()) {
1541  // This linear constraint has to be marked redundant.
1542  unsigned int redundantLC = listOfGen - _listOfGeneratorsPerLinearConstraint.begin();
1543  _listOfRedundantLinearConstraints.insert(redundantLC);
1544  // Run the list of the generators belonging to this linear constraint to unhook them.
1545  {
1546 #ifdef DEBUG
1547  std::cout << "Redundant linear constraint " << redundantLC << std::endl;
1548 #endif
1549  std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator itGN = listOfGen->begin();
1550  for (; itGN != listOfGen->end(); ++itGN) {
1551 #ifdef DEBUG
1552  std::cout << " To be removed from ";
1553  (*itGN)->dump(std::cout);
1554  std::cout << std::endl;
1555 #endif
1556  (*itGN)->removeCurrentLinearConstraint(redundantLC);
1557  }
1558  }
1559  // Then delete all traces of this linear constraint.
1560  listOfGen->clear();
1561  } // if (listOfGen->size() < poly->numberOfGeneratorsPerFacet()) {
1562  }
1563  }
1564  // Declare a linear constraint redundant if the list of its generators is included into another linear constraint generator list.
1565  {
1566  std::vector< std::set< boost::shared_ptr< Generator_Rn_SD > > >::iterator listOfGen1 = _listOfGeneratorsPerLinearConstraint.begin();
1567  for (; listOfGen1 != _listOfGeneratorsPerLinearConstraint.end(); ++listOfGen1) {
1568  unsigned int redundantLC1 = listOfGen1 - _listOfGeneratorsPerLinearConstraint.begin();
1569  std::vector< std::set< boost::shared_ptr< Generator_Rn_SD > > >::iterator listOfGen2 = listOfGen1 + 1;
1570  for (; listOfGen2 != _listOfGeneratorsPerLinearConstraint.end(); ++listOfGen2) {
1571  unsigned int redundantLC2 = listOfGen2 - _listOfGeneratorsPerLinearConstraint.begin();
1572  if (!listOfGen1->empty() && !listOfGen2->empty() &&
1573  listOfGen1->size() > listOfGen2->size() &&
1574  std::includes(listOfGen1->begin(), listOfGen1->end(), listOfGen2->begin(), listOfGen2->end())) {
1575  // Run the list of the generators belonging to this linear constraint to unhook them.
1576 #ifdef DEBUG
1577  std::cout << "Redundant linear constraint " << redundantLC2 << " included into linear constraint " << redundantLC1 << std::endl;
1578 #endif
1579  _listOfRedundantLinearConstraints.insert(redundantLC2);
1580  std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator itGN = listOfGen2->begin();
1581  for (; itGN != listOfGen2->end(); ++itGN) {
1582 #ifdef DEBUG
1583  std::cout << " To be removed from ";
1584  (*itGN)->dump(std::cout);
1585  std::cout << std::endl;
1586 #endif
1587  (*itGN)->removeCurrentLinearConstraint(redundantLC2);
1588  }
1589  // Then delete all traces of this linear constraint.
1590  listOfGen2->clear();
1591  }
1592  else if (!listOfGen1->empty() && !listOfGen2->empty() &&
1593  listOfGen1->size() <= listOfGen2->size() &&
1594  std::includes(listOfGen2->begin(), listOfGen2->end(), listOfGen1->begin(), listOfGen1->end())) {
1595  // Run the list of the generators belonging to this linear constraint to unhook them.
1596 #ifdef DEBUG
1597  std::cout << "Redundant linear constraint " << redundantLC1 << " included into linear constraint " << redundantLC2 << std::endl;
1598 #endif
1599  _listOfRedundantLinearConstraints.insert(redundantLC1);
1600  std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator itGN = listOfGen1->begin();
1601  for (; itGN != listOfGen1->end(); ++itGN) {
1602 #ifdef DEBUG
1603  std::cout << " To be removed from ";
1604  (*itGN)->dump(std::cout);
1605  std::cout << std::endl;
1606 #endif
1607  (*itGN)->removeCurrentLinearConstraint(redundantLC1);
1608  }
1609  // Then delete all traces of this linear constraint.
1610  listOfGen1->clear();
1611  }
1612  }
1613  }
1614  }
1615  } // else if (nbGeneratorsIn != 0) {
1616 
1617 #ifdef DEBUG
1618  //std::cout << "VX_new= " << GN_new.size() << std::endl;
1619  //std::cout << "VX_OUT= " << GN_OUT.size() << std::endl;
1620  //std::cout << "VX_IN = " << GN_IN.size() << std::endl;
1621  //std::cout << "VX_ON = " << GN_ON.size() << std::endl;
1622  //std::cout << "redundantHS = " << _listOfRedundantLinearConstraints.size() << std::endl;
1623  //std::cout << "============================================" << std::endl << std::endl;
1624 #endif
1625  // When we perform a flat truncation the set of generators IN
1626  // . is transferred into the set of generators OUT as both of them will be dropped
1627  // . then it is emptied before receiving the generators ON and NEW i.e. the ones we keep
1628  // In the other case of an intersection with a half-space this function does nothing.
1629  //currentHalfSpace->prepareSetOfGenerators(GN_IN, GN_OUT);
1630 
1631  // All the generators with state ON must have a new facet.
1632  //std::vector< boost::shared_ptr<Generator_Rn_SD> >::const_iterator iteGN_ON;
1633  //for (iteGN_ON=GN_ON.begin(); iteGN_ON!=GN_ON.end(); ++iteGN_ON) {
1634  // (*iteGN_ON)->setFacet(halfspaceNumber);
1635  // if ((*iteGN_ON)->getStatus() == Generator_Rn_SD::CREATED)
1636  // (*iteGN_ON)->setStatus(Generator_Rn_SD::CREATED_AND_MODIFIED);
1637  // else if ((*iteGN_ON)->getStatus() == Generator_Rn_SD::UNCHANGED)
1638  // (*iteGN_ON)->setStatus(Generator_Rn_SD::MODIFIED);
1639  // Insert the ON generators in the IN list as this list will be the official one soon.
1640  // GN_IN.push_back(*iteGN_ON);
1641  //}
1642 
1643 #ifdef DEBUG
1644  std::cout << "VX_OUT= " << nbGeneratorsOut << std::endl;
1645  std::cout << "VX_IN = " << nbGeneratorsIn << std::endl;
1646  std::cout << "VX_ON = " << nbGeneratorsOn << std::endl;
1647  std::cout << "VX_new= " << GN_new.size() << std::endl;
1648  std::cout << "redundantHS = " << _listOfRedundantLinearConstraints.size() << std::endl;
1649  std::cout << "============================================" << std::endl << std::endl;
1650 #endif
1651  GN_OUT.clear();
1652  GN_IN.clear();
1653  GN_ON.clear();
1654  GN_new.clear();
1655  halfspaceNumber++;
1656  // Job related to doing the search for the best chopping half-space.
1657  _currentSegmentIndexAllNeighbours++;
1658  while (_allNeighbours[ _currentSegmentIndexAllNeighbours ].getActive() == false || _currentSegmentIndexAllNeighbours >= _allNeighbours.size()) {
1659  _currentSegmentIndexAllNeighbours++;
1660  if (_currentSegmentIndexAllNeighbours >= _allNeighbours.size()) {
1661  _currentSegmentIndexAllNeighbours = 0;
1662  // Safety check to make sure not all the elements in _allNeighbours are inactive.
1663  bool active = false;
1664  for (unsigned int seg = 0; seg < _allNeighbours.size(); ++seg) {
1665  if (_allNeighbours[seg].getActive() == true) {
1666  active = true;
1667  _currentSegmentIndexAllNeighbours = seg;
1668  break;
1669  }
1670  }
1671  if (active == false) {
1672 #ifdef DEBUG
1673  std::cout << "All segments all inactive (2)." << std::endl;
1674 #endif
1675  _isEmpty = true;
1676  poly->reset();
1677  return false;
1678  }
1679  }
1680  } // while (_allNeighbours[ _currentSegmentIndexAllNeighbours ].getActive() == false) {
1681 
1682 #ifdef DEBUG
1683  checkDataStructuresConsistency();
1684 #endif
1685 
1686  } // for (iteHS.begin(); iteHS.end()!=true; iteHS.next()) {
1687  }
1688 
1689 #ifdef DEBUG
1690  std::cout << " // Display redundancy processing //" << std::endl;
1691  {
1692  std::vector< std::set< boost::shared_ptr< Generator_Rn_SD > > >::iterator listOfGen1 = _listOfGeneratorsPerLinearConstraint.begin();
1693  for (; listOfGen1 != _listOfGeneratorsPerLinearConstraint.end(); ++listOfGen1) {
1694  std::cout << "LC" << listOfGen1 - _listOfGeneratorsPerLinearConstraint.begin() << ": ";
1695  std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator itGN1 = listOfGen1->begin();
1696  std::vector< unsigned int > gnNbinLC;
1697  for (; itGN1 != listOfGen1->end(); ++itGN1) {
1698  //(*itGN1)->dump(std::cout);
1699  //std::cout << "; ";
1700  gnNbinLC.push_back( (*itGN1)->getGeneratorNumber() );
1701  //std::cout << (*itGN1)->getGeneratorNumber() << " ";
1702  }
1703  std::sort(gnNbinLC.begin(), gnNbinLC.end());
1704  std::copy(gnNbinLC.begin(), gnNbinLC.end(), std::ostream_iterator< unsigned int >(std::cout, " ") );
1705  if (listOfGen1 - _listOfGeneratorsPerLinearConstraint.begin() > truncationStep && listOfGen1 - _listOfGeneratorsPerLinearConstraint.begin() <= halfspaceNumber)
1706  _listOfLinearConstraints[ listOfGen1 - _listOfGeneratorsPerLinearConstraint.begin() ]->dump(std::cout);
1707  std::cout << std::endl;
1708  }
1709  }
1710 #endif
1711 
1712  // Build the final list of generators without duplications.
1713  {
1714  listOfGenSD.clear();
1715  std::vector< Segment_Rn >::iterator it_all1;
1716  std::set< boost::shared_ptr<Generator_Rn_SD> > setOfUniqueGen;
1717  for (it_all1 = _allNeighbours.begin(); it_all1 != _allNeighbours.end(); ++it_all1) {
1718  if (it_all1->getActive() == true) {
1719  if (setOfUniqueGen.find( it_all1->getFirstGenerator() ) == setOfUniqueGen.end()) {
1720  setOfUniqueGen.insert( it_all1->getFirstGenerator() );
1721  listOfGenSD.push_back( it_all1->getFirstGenerator() );
1722  }
1723  if (setOfUniqueGen.find( it_all1->getSecondGenerator() ) == setOfUniqueGen.end()) {
1724  setOfUniqueGen.insert( it_all1->getSecondGenerator() );
1725  listOfGenSD.push_back( it_all1->getSecondGenerator() );
1726  }
1727  }
1728  }
1729  //dumpListOfGenerators(poly, listOfGenSD, std::cout);
1730  //std::cout << std::endl;
1731  }
1732 
1733  return true;
1734  }
1735 
1736  static void dumpListOfGenerators(POLYHEDRON poly, std::vector< boost::shared_ptr<Generator_Rn_SD> > listOfGenSD, std::ostream &this_ostream, bool displayAll = true) {
1737  std::vector< boost::shared_ptr<Generator_Rn_SD> >::iterator it_all1;
1738  for (it_all1 = listOfGenSD.begin(); it_all1 != listOfGenSD.end(); ++it_all1) {
1739  (*it_all1)->dump(this_ostream);
1740  this_ostream << std::endl;
1741  if (displayAll == true) {
1742  std::vector< unsigned int >::const_iterator it = (*it_all1)->facetsBegin();
1743  for (; it != (*it_all1)->facetsEnd(); ++it) {
1744  poly->getHalfSpace(*it)->dump(std::cout);
1745  this_ostream << " ";
1746  }
1747  this_ostream << std::endl;
1748  }
1749  }
1750  }
1751 
1754  bool checkNeighbours(
1755  POLYHEDRON poly,
1756  const boost::shared_ptr<Generator_Rn_SD>& genIn,
1757  const boost::shared_ptr<Generator_Rn_SD>& genOut,
1758  std::vector< unsigned int >& commonFacets) {
1759  // Compute the intersection on sorted lists of facets.
1760  std::set_intersection(genIn->facetsBegin(), genIn->facetsEnd(), genOut->facetsBegin(), genOut->facetsEnd(), std::inserter(commonFacets, commonFacets.end()));
1761  // For polyhedral cones neigbourhoodCondition() is Rn::getDimension()-1, for polytopes it is Rn::getDimension()-2.
1762  if (commonFacets.size() >= _neigbourhoodCondition)
1763  return true;
1764 
1765  return false;
1766  }
1767 
1768  // Mark as UNCHANGED, CREATED or DELETED the half-spaces going through the double description process.
1769  void markHdescription(TrackingOperatorToResult& trackerHdesc, unsigned int truncationStep, unsigned int numberOfProcessedLinearConstraints) {
1770  trackerHdesc.setNumbersOfEntities(numberOfProcessedLinearConstraints, numberOfProcessedLinearConstraints);
1771  // Before truncation step they are UNCHANGED by default, then they are CREATED by default.
1772  for (unsigned int i=0; i<truncationStep; ++i)
1773  trackerHdesc.setOperatorEntityAsUnchanged(i);
1774  for (unsigned int i=truncationStep; i<numberOfProcessedLinearConstraints; ++i)
1775  trackerHdesc.setResultEntityAsCreated(i);
1776  std::set< unsigned int >::const_iterator iteRED=_listOfRedundantLinearConstraints.begin();
1777  for ( ; iteRED!=_listOfRedundantLinearConstraints.end(); ++iteRED)
1778  trackerHdesc.setOperatorEntityAsDeleted(*iteRED);
1779  }
1780 
1783  void unhookRedundantGenerators(POLYHEDRON poly) {
1784  unsigned int minNbFacetsPerGen = poly->neigbourhoodCondition() + 1;
1785  std::set< boost::shared_ptr<Generator_Rn> > setToRemove;
1786  constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(poly->getListOfGenerators());
1787  {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
1788  //std::cout << "gen->numberOfFacets() = " << iteGN.current()->numberOfFacets();
1789  if (iteGN.current()->numberOfFacets() < minNbFacetsPerGen) {
1790  setToRemove.insert( iteGN.current() );
1791  //std::cout << " removed! ";
1792  }
1793  }}
1794  if (!setToRemove.empty())
1795  poly->removeGenerators(setToRemove);
1796  //std::cout << std::endl;
1797  }
1798 
1800  bool checkFuzziness(
1801  POLYHEDRON poly,
1802  const boost::shared_ptr<Generator_Rn_SD>& gen1,
1803  const boost::shared_ptr<Generator_Rn_SD>& fuzz,
1804  std::vector< unsigned int >& commonFacets) {
1805  // Compute the intersection on sorted lists of facets.
1806  std::set_intersection(gen1->facetsBegin(), gen1->facetsEnd(), fuzz->fuzzyFacetsBegin(), fuzz->fuzzyFacetsEnd(), std::inserter(commonFacets, commonFacets.end()));
1807  //std::cout << "GEN1: ";
1808  //std::copy(gen1->facetsBegin(), gen1->facetsEnd(), std::ostream_iterator< unsigned int >(std::cout, " ") );
1809  //std::cout << "FUZZ: ";
1810  //std::copy(fuzz->fuzzyFacetsBegin(), fuzz->fuzzyFacetsEnd(),std::ostream_iterator< unsigned int >(std::cout, " ") );
1811  //std::cout << "COMM: ";
1812  //std::copy(commonFacets.begin(), commonFacets.end(),std::ostream_iterator< unsigned int >(std::cout, " ") );
1813  //std::cout << std::endl;
1814  // For polyhedral cones neigbourhoodCondition() is Rn::getDimension()-1, for polytopes it is Rn::getDimension()-2.
1815  if (commonFacets.empty())
1816  return false;
1817 
1818  return true;
1819  }
1820 
1821  protected:
1823  bool checkDataStructuresConsistency() {
1824  std::vector< std::set< boost::shared_ptr< Generator_Rn_SD > > > list2CompareOfGeneratorsPerLinearConstraint;
1825  {
1826  unsigned int nbF = _listOfLinearConstraintNorms.size();
1827  for (unsigned int i=0; i<nbF; ++i) {
1828  std::set< boost::shared_ptr< Generator_Rn_SD > > genSet;
1829  list2CompareOfGeneratorsPerLinearConstraint.push_back(genSet);
1830  }
1831  }
1832  std::set< boost::shared_ptr< Generator_Rn_SD > > alreadyProcessedGen;
1833  std::vector< Segment_Rn >::iterator it_GN;
1834  for (it_GN = _allNeighbours.begin(); it_GN != _allNeighbours.end(); ++it_GN) {
1835  if (it_GN->getActive() == true) {
1836  std::pair< std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator, bool > ret;
1837  ret = alreadyProcessedGen.insert(it_GN->getFirstGenerator());
1838  //std::cout << "nb = " << it_GN - _allNeighbours.begin() << std::endl;
1839  if (ret.second == true) {
1840  for (unsigned int ii=0; ii<it_GN->getFirstGenerator()->numberOfFacets(); ii++)
1841  list2CompareOfGeneratorsPerLinearConstraint[ it_GN->getFirstGenerator()->getFacet(ii) ].insert( it_GN->getFirstGenerator() );
1842  }
1843  //std::cout << "nb1 " << std::endl;
1844  ret = alreadyProcessedGen.insert(it_GN->getSecondGenerator());
1845  if (ret.second == true) {
1846  for (unsigned int ii=0; ii<it_GN->getSecondGenerator()->numberOfFacets(); ii++)
1847  list2CompareOfGeneratorsPerLinearConstraint[ it_GN->getSecondGenerator()->getFacet(ii) ].insert( it_GN->getSecondGenerator() );
1848  }
1849  //std::cout << "nb2 " << std::endl;
1850  }
1851  }
1852  std::cout << "======================" << std::endl;
1853  std::cout << "= checkConsistency() =" << std::endl;
1854  std::cout << "======================" << std::endl;
1855  bool OK1 = true;
1856  {
1857  unsigned int nbF = list2CompareOfGeneratorsPerLinearConstraint.size();
1858  for (unsigned int i=0; i<nbF; ++i) {
1859  if (list2CompareOfGeneratorsPerLinearConstraint[i] != _listOfGeneratorsPerLinearConstraint[i]) {
1860  OK1 = false;
1861  std::cout << "Problem found for constraint " << i << ": ";
1862  _listOfLinearConstraints[i]->dump(std::cout);
1863  std::cout << std::endl << "Lists of generators for both data structures: " << std::endl;
1864  {
1865  std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator itGN1 = list2CompareOfGeneratorsPerLinearConstraint[i].begin();
1866  for (; itGN1 != list2CompareOfGeneratorsPerLinearConstraint[i].end(); ++itGN1)
1867  std::cout << (*itGN1)->getGeneratorNumber() << " ";
1868  std::cout << std::endl;
1869  }
1870  {
1871  std::set< boost::shared_ptr< Generator_Rn_SD > >::iterator itGN1 = _listOfGeneratorsPerLinearConstraint[i].begin();
1872  for (; itGN1 != _listOfGeneratorsPerLinearConstraint[i].end(); ++itGN1)
1873  std::cout << (*itGN1)->getGeneratorNumber() << " ";
1874  std::cout << std::endl;
1875  }
1876  }
1877  }
1878  }
1879  (OK1 == true) ? std::cout << "OK1 " : std::cout << "KO1 ";
1880  std::cout << std::endl << std::endl;
1881 
1882  return OK1;
1883  }
1884 
1885  protected:
1886  double _TOL;
1887  unsigned int _RnDIM;
1889  bool _isEmpty;
1890  unsigned int _neigbourhoodCondition;
1892  std::vector< std::set< boost::shared_ptr< Generator_Rn_SD > > > _listOfGeneratorsPerLinearConstraint;
1893  std::set< unsigned int > _listOfRedundantLinearConstraints;
1894  std::vector< boost::shared_ptr<HalfSpace_Rn> > _listOfLinearConstraints;
1895  std::vector< double > _listOfLinearConstraintNorms;
1897  std::vector< boost::shared_ptr<NormedHalfSpace_Rn> > _listOfNormedLinearConstraints;
1899  std::vector< Segment_Rn > _allNeighbours;
1900  unsigned int _currentSegmentIndexAllNeighbours;
1901 
1902 };
1903 
1904 
1905 #endif
HalfSpace_Rn::hs_OUT
@ hs_OUT
Definition: HalfSpace_Rn.h:47
DoubleDescription::computeGeneratorState
HalfSpace_Rn::State computeGeneratorState(const boost::shared_ptr< Generator_Rn_SD > &Gen, const boost::shared_ptr< HalfSpace_Rn > &currentHalfSpace, unsigned int halfspaceNumber, double halfSpaceNorm, unsigned int &nbGeneratorsIn, unsigned int &nbGeneratorsOut, unsigned int &nbGeneratorsOn)
Compute the state hs_ON,hs_IN,hs_OUT,hs_UNKNOWN,hs_IN_OR_OUT for a given generator.
Definition: DoubleDescription_Rn.h:377
HalfSpace_Rn::hs_UNKNOWN
@ hs_UNKNOWN
Definition: HalfSpace_Rn.h:48
DoubleDescription::checkNeighbours
bool checkNeighbours(POLYHEDRON poly, const boost::shared_ptr< Generator_Rn_SD > &genIn, const boost::shared_ptr< Generator_Rn_SD > &genOut, std::vector< unsigned int > &commonFacets)
Definition: DoubleDescription_Rn.h:1754
Segment_Rn::getSecondGenerator
boost::shared_ptr< Generator_Rn_SD > & getSecondGenerator()
Definition: DoubleDescription_Rn.h:44
Neighbours_Rn.h
TrackingOperatorToResult::setResult_Operator
void setResult_Operator(unsigned int nbRes, int nb)
Set the link between the nb-th entity of operator1 and the nbRes-th entity of the result.
Definition: Tracking.h:89
TrackingOperatorToResult::setResultEntityAsUnchanged
void setResultEntityAsUnchanged(unsigned int nb)
Mark as Result_UNCHANGED the nb-th entity before the operation.
Definition: Tracking.h:67
DoubleDescription::DoubleDescription
DoubleDescription(POLYHEDRON poly, ITERATOR ite, int truncationStep, TrackingOperatorToResult &trackerVdesc, TrackingOperatorToResult &trackerHdesc)
Definition: DoubleDescription_Rn.h:171
Neighbours_Rn::end
bool end() const
Iterator function.
Definition: Neighbours_Rn.h:151
HalfSpace_Rn::hs_IN
@ hs_IN
Definition: HalfSpace_Rn.h:46
DoubleDescription::dumpListOfGenerators
static void dumpListOfGenerators(POLYHEDRON poly, std::vector< boost::shared_ptr< Generator_Rn_SD > > listOfGenSD, std::ostream &this_ostream, bool displayAll=true)
Definition: DoubleDescription_Rn.h:1736
Generator_Rn_SD::CREATED_AND_MODIFIED
@ CREATED_AND_MODIFIED
Definition: Generator_Rn.h:293
UpdDoubleDescription_Rn.h
DoubleDescription::inferior
static bool inferior(const boost::shared_ptr< Generator_Rn_SD > &HS1, const boost::shared_ptr< Generator_Rn_SD > &HS2)
Internal function: tell whether a given object is declared inferior to another one.
Definition: DoubleDescription_Rn.h:362
HalfSpace_Rn::hs_IN_OR_OUT
@ hs_IN_OR_OUT
Definition: HalfSpace_Rn.h:49
Generator_Rn_SD::Status
Status
Definition: Generator_Rn.h:285
DoubleDescription::_listOfNormedLinearConstraints
std::vector< boost::shared_ptr< NormedHalfSpace_Rn > > _listOfNormedLinearConstraints
To speed up computations.
Definition: DoubleDescription_Rn.h:1897
TrackingOperatorToResult::setOperatorEntityAsDeleted
void setOperatorEntityAsDeleted(unsigned int nb)
Mark as Operator_DELETED the nb-th entity before the operation.
Definition: Tracking.h:64
Neighbours_Rn::currentGenInNumber
unsigned int currentGenInNumber() const
Iterator function.
Definition: Neighbours_Rn.h:154
DoubleDescription::findBestLinearConstraint
bool findBestLinearConstraint(std::vector< boost::shared_ptr< HalfSpace_Rn > >::iterator currentLC, std::vector< boost::shared_ptr< HalfSpace_Rn > >::iterator bestLC)
Given a list of half-spaces, try to find the one able to cut most of a given segment or to discard it...
Definition: DoubleDescription_Rn.h:427
Segment_Rn
Define a segment whose ends are generators.
Definition: DoubleDescription_Rn.h:36
Generator_Rn_SD::MODIFIED
@ MODIFIED
Definition: Generator_Rn.h:289
DoubleDescription::_listOfRedundantLinearConstraints
std::set< unsigned int > _listOfRedundantLinearConstraints
Definition: DoubleDescription_Rn.h:1893
Neighbours_Rn::next
void next()
Iterator function.
Definition: Neighbours_Rn.h:142
DoubleDescription::checkDataStructuresConsistency
bool checkDataStructuresConsistency()
In debug mode only, check the equality of two data structures.
Definition: DoubleDescription_Rn.h:1823
Neighbours_Rn::begin
void begin()
Iterator function.
Definition: Neighbours_Rn.h:139
TrackingOperatorToResult::setNumbersOfEntities
void setNumbersOfEntities(unsigned int nbEntBefore, unsigned int nbEntAfter)
Definition: Tracking.h:47
Segment_Rn::setSecondGenerator
void setSecondGenerator(const boost::shared_ptr< Generator_Rn_SD > &e2)
Definition: DoubleDescription_Rn.h:46
constIteratorOfListOfGeometricObjects::current
const GEOMETRIC_OBJECT current()
Return the current geometric element.
Definition: GeometricObjectIterator_Rn.h:177
Segment_Rn::Segment_Rn
Segment_Rn(const std::vector< unsigned int > &LC)
Definition: DoubleDescription_Rn.h:40
Neighbours_Rn
Class dedicated to degeneration processing when looking for neighbours. Let A be a polytope of wher...
Definition: Neighbours_Rn.h:59
constIteratorOfListOfGeometricObjects::end
bool end() const
Tell whether we have reached the end of the list.
Definition: GeometricObjectIterator_Rn.h:174
Neighbours_Rn::addGenerator
void addGenerator(const std::vector< unsigned int > &commonFacets, unsigned int numbergenIN, unsigned int numbergenOUT, HalfSpace_Rn::State state)
Tell whether a pseudo neighbor is a genuine one comparing set of half-spaces.
Definition: Neighbours_Rn.h:69
DoubleDescription::_listOfLinearConstraintNorms
std::vector< double > _listOfLinearConstraintNorms
Definition: DoubleDescription_Rn.h:1895
DoubleDescription::unhookRedundantGenerators
void unhookRedundantGenerators(POLYHEDRON poly)
Definition: DoubleDescription_Rn.h:1783
Neighbours_Rn::currentGenOutNumber
unsigned int currentGenOutNumber() const
Iterator function.
Definition: Neighbours_Rn.h:157
TrackingOperatorToResult::setResultEntityAsCreated
void setResultEntityAsCreated(unsigned int nb)
Mark as Result_CREATED the nb-th entity after the operation.
Definition: Tracking.h:76
Segment_Rn::Segment_Rn
Segment_Rn()
Definition: DoubleDescription_Rn.h:39
TrackingOperatorToResult::setOperatorEntityAsUnchanged
void setOperatorEntityAsUnchanged(unsigned int nb)
Mark as Operator_UNCHANGED the nb-th entity before the operation.
Definition: Tracking.h:58
Generator_Rn_SD::CREATED
@ CREATED
Definition: Generator_Rn.h:291
NormedHalfSpace_Rn
Definition: HalfSpace_Rn.h:172
DoubleDescription::markHdescription
void markHdescription(TrackingOperatorToResult &trackerHdesc, unsigned int truncationStep, unsigned int numberOfProcessedLinearConstraints)
Definition: DoubleDescription_Rn.h:1769
DoubleDescription::_RnDIM
unsigned int _RnDIM
Definition: DoubleDescription_Rn.h:1887
HalfSpace_Rn::getStateAsText
static std::string getStateAsText(const HalfSpace_Rn::State &)
Definition: HalfSpace_Rn.cpp:55
DoubleDescription::computeStatesOfAllGenerators
void computeStatesOfAllGenerators(std::vector< boost::shared_ptr< Generator_Rn_SD > > &arrayOfGenOUT, const boost::shared_ptr< HalfSpace_Rn > &currentHalfSpace, unsigned int halfspaceNumber, unsigned int &nbGeneratorsIn, unsigned int &nbGeneratorsOut, unsigned int &nbGeneratorsOn)
For each generator compute its state according to the current half-space.
Definition: DoubleDescription_Rn.h:280
constIteratorOfListOfGeometricObjects
This class is designed to run the list of all geometric objects representing a polytope.
Definition: GeometricObjectIterator_Rn.h:142
DoubleDescription::_TOL
double _TOL
Definition: DoubleDescription_Rn.h:1886
Segment_Rn::setFirstGenerator
void setFirstGenerator(const boost::shared_ptr< Generator_Rn_SD > &e1)
Definition: DoubleDescription_Rn.h:45
DoubleDescription
The algorithm implemented here is an incremental algorithm as mentioned in How Good are Convex Hull ...
Definition: DoubleDescription_Rn.h:86
HalfSpace_Rn::State
State
Definition: HalfSpace_Rn.h:44
DoubleDescription::checkFuzziness
bool checkFuzziness(POLYHEDRON poly, const boost::shared_ptr< Generator_Rn_SD > &gen1, const boost::shared_ptr< Generator_Rn_SD > &fuzz, std::vector< unsigned int > &commonFacets)
Check whether a fuzzy half-space of the generator fuzz can be a half-space of gen1 when we chop the s...
Definition: DoubleDescription_Rn.h:1800
Generator_Rn_SD
A n-coordinates generator for internal data structure. It can be a vertex or an edge whether it is em...
Definition: Generator_Rn.h:280
Segment_Rn::getActive
bool getActive()
Definition: DoubleDescription_Rn.h:47
constIteratorOfListOfGeometricObjects::begin
void begin()
Move the iterator at the beginning of the list.
Definition: GeometricObjectIterator_Rn.h:150
Neighbours_Rn::dump
void dump(std::ostream &ofs) const
Display the content on the stream passed as an argument.
Definition: Neighbours_Rn.h:163
DoubleDescription::_currentSegmentIndexAllNeighbours
unsigned int _currentSegmentIndexAllNeighbours
Definition: DoubleDescription_Rn.h:1900
Tracking.h
DoubleDescription::DoubleDescription
DoubleDescription(POLYHEDRON poly, ITERATOR ite, int truncationStep)
Definition: DoubleDescription_Rn.h:89
Segment_Rn::Segment_Rn
Segment_Rn(boost::shared_ptr< Generator_Rn_SD > &e1, boost::shared_ptr< Generator_Rn_SD > &e2, const std::vector< unsigned int > &LC)
Definition: DoubleDescription_Rn.h:42
Segment_Rn::getFirstGenerator
boost::shared_ptr< Generator_Rn_SD > & getFirstGenerator()
Definition: DoubleDescription_Rn.h:43
DoubleDescription::_isEmpty
bool _isEmpty
Store the current state of the intersection.
Definition: DoubleDescription_Rn.h:1889
DoubleDescription::_listOfLinearConstraints
std::vector< boost::shared_ptr< HalfSpace_Rn > > _listOfLinearConstraints
Definition: DoubleDescription_Rn.h:1894
HalfSpace_Rn::hs_ON
@ hs_ON
Definition: HalfSpace_Rn.h:45
TrackingOperatorToResult::setOperator_Result
void setOperator_Result(unsigned int nb, int nbRes)
Set the link between the nb-th entity of operator1 and the nbRes-th entity of the result.
Definition: Tracking.h:86
HalfSpace_Rn::hs_OUT_TO_ON
@ hs_OUT_TO_ON
Definition: HalfSpace_Rn.h:50
Neighbours_Rn::commonHalfSpaces
const std::vector< unsigned int > & commonHalfSpaces() const
Iterator function.
Definition: Neighbours_Rn.h:160
TrackingOperatorToResult::setOperatorEntityAsModified
void setOperatorEntityAsModified(unsigned int nb)
Mark as Operator_MODIFIED the nb-th entity before the operation.
Definition: Tracking.h:61
DoubleDescription::_listOfGeneratorsPerLinearConstraint
std::vector< std::set< boost::shared_ptr< Generator_Rn_SD > > > _listOfGeneratorsPerLinearConstraint
For each linear constraint, update and store the list of its associated generators.
Definition: DoubleDescription_Rn.h:1892
Segment_Rn::setActive
void setActive()
Definition: DoubleDescription_Rn.h:48
constIteratorOfListOfGeometricObjects::next
void next()
Move the iterator one step forward.
Definition: GeometricObjectIterator_Rn.h:153
Generator_Rn_SD::UNCHANGED
@ UNCHANGED
Definition: Generator_Rn.h:287
Segment_Rn::Segment_Rn
Segment_Rn(boost::shared_ptr< Generator_Rn_SD > &e1, boost::shared_ptr< Generator_Rn_SD > &e2)
Definition: DoubleDescription_Rn.h:41
Segment_Rn::dump
void dump(std::ostream &ofs)
Definition: DoubleDescription_Rn.h:50
TrackingOperatorToResult::setResultEntityAsModified
void setResultEntityAsModified(unsigned int nb)
Mark as Result_MODIFIED the nb-th entity after the operation.
Definition: Tracking.h:70
DoubleDescription::getIsEmpty
bool getIsEmpty() const
Definition: DoubleDescription_Rn.h:277
Rn::getTolerance
static polito_EXPORT double getTolerance()
Give the minimum distance between two points.
Definition: Rn.cpp:31
DoubleDescription::_allNeighbours
std::vector< Segment_Rn > _allNeighbours
The connections between generators { {gen_a, gen_b}, {gen_a, gen_c}, ...}.
Definition: DoubleDescription_Rn.h:1899
DoubleDescription::~DoubleDescription
virtual ~DoubleDescription()
Definition: DoubleDescription_Rn.h:245
DoubleDescription::_neigbourhoodCondition
unsigned int _neigbourhoodCondition
Definition: DoubleDescription_Rn.h:1890
DoubleDescription::compute
bool compute(POLYHEDRON poly, ITERATOR iteHS, int truncationStep, std::vector< boost::shared_ptr< Generator_Rn_SD > > &listOfGenSD)
The main function splitting the polyhedron cone or polytope 1-skeleton with a list of half-spaces.
Definition: DoubleDescription_Rn.h:506
Segment_Rn::setInactive
void setInactive()
Definition: DoubleDescription_Rn.h:49
TrackingOperatorToResult
This class stores static function that dispatch the main geometric values we use.
Definition: Tracking.h:41
DoubleDescription::init
unsigned int init(POLYHEDRON poly, ITERATOR ite, int truncationStep)
Definition: DoubleDescription_Rn.h:247