Actual source code: Sieve.hh

  2: #ifndef included_ALE_Sieve_hh
  3: #define included_ALE_Sieve_hh

  5: #include <boost/multi_index_container.hpp>
  6: #include <boost/multi_index/member.hpp>
  7: #include <boost/multi_index/ordered_index.hpp>
  8: #include <boost/multi_index/composite_key.hpp>
  9: #include <iostream>

 11: #ifndef  included_ALE_hh
 12: #include <ALE.hh>
 13: #endif

 15: // ALE extensions

 17: #ifndef  included_ALE_Sifter_hh
 18: #include <Sifter.hh>
 19: #endif


 22: namespace ALE {

 24:   namespace SieveDef {
 25:     //
 26:     // Rec & RecContainer definitions.
 27:     // Rec is intended to denote a graph point record.
 28:     //
 29:     template <typename Point_, typename Marker_>
 30:     struct Rec {
 31:       typedef Point_  point_type;
 32:       typedef Marker_ marker_type;
 33:       template<typename OtherPoint_, typename OtherMarker_ = Marker_>
 34:       struct rebind {
 35:         typedef Rec<OtherPoint_, OtherMarker_> type;
 36:       };
 37:       point_type  point;
 38:       int         degree;
 39:       int         depth;
 40:       int         height;
 41:       marker_type marker;
 42:       // Basic interface
 43:       Rec() : point(point_type()), degree(0), depth(0), height(0), marker(marker_type()) {};
 44:       Rec(const Rec& r) : point(r.point), degree(r.degree), depth(r.depth), height(r.height), marker(r.marker) {}
 45:       Rec(const point_type& p) : point(p), degree(0), depth(0), height(0), marker(marker_type()) {};
 46:       Rec(const point_type& p, const int degree) : point(p), degree(degree), depth(0), height(0), marker(marker_type()) {};
 47:       Rec(const point_type& p, const int degree, const int depth, const int height, const marker_type marker) : point(p), degree(degree), depth(depth), height(height), marker(marker) {};
 48:       // Printing
 49:       friend std::ostream& operator<<(std::ostream& os, const Rec& p) {
 50:         os << "<" << p.point << ", "<< p.degree << ", "<< p.depth << ", "<< p.height << ", "<< p.marker << ">";
 51:         return os;
 52:       };

 54:       struct degreeAdjuster {
 55:         degreeAdjuster(int newDegree) : _newDegree(newDegree) {};
 56:         void operator()(Rec& r) { r.degree = this->_newDegree; }
 57:       private:
 58:         int _newDegree;
 59:       };
 60:     };// class Rec

 62:     template <typename Point_, typename Rec_>
 63:     struct RecContainerTraits {
 64:       typedef Rec_ rec_type;
 65:       typedef typename rec_type::marker_type marker_type;
 66:       // Index tags
 67:       struct pointTag{};
 68:       struct degreeTag{};
 69:       struct markerTag{};
 70:       struct depthMarkerTag{};
 71:       struct heightMarkerTag{};
 72:       // Rec set definition
 73:       typedef ::boost::multi_index::multi_index_container<
 74:         rec_type,
 75:         ::boost::multi_index::indexed_by<
 76:           ::boost::multi_index::ordered_unique<
 77:             ::boost::multi_index::tag<pointTag>, BOOST_MULTI_INDEX_MEMBER(rec_type, typename rec_type::point_type, point)
 78:           >,
 79: //           ::boost::multi_index::ordered_non_unique<
 80: //             ::boost::multi_index::tag<degreeTag>, BOOST_MULTI_INDEX_MEMBER(rec_type, int, degree)
 81: //           >,
 82:           ::boost::multi_index::ordered_non_unique<
 83:             ::boost::multi_index::tag<markerTag>, BOOST_MULTI_INDEX_MEMBER(rec_type, marker_type, marker)
 84:           >,
 85:           ::boost::multi_index::ordered_non_unique<
 86:             ::boost::multi_index::tag<depthMarkerTag>,
 87:             ::boost::multi_index::composite_key<
 88:               rec_type, BOOST_MULTI_INDEX_MEMBER(rec_type,int,depth), BOOST_MULTI_INDEX_MEMBER(rec_type,marker_type,marker)>
 89:           >,
 90:           ::boost::multi_index::ordered_non_unique<
 91:             ::boost::multi_index::tag<heightMarkerTag>,
 92:             ::boost::multi_index::composite_key<
 93:               rec_type, BOOST_MULTI_INDEX_MEMBER(rec_type,int,height), BOOST_MULTI_INDEX_MEMBER(rec_type,marker_type,marker)>
 94:           >
 95:         >,
 96:         ALE_ALLOCATOR<rec_type>
 97:       > set_type;
 98:       //
 99:       // Return types
100:       //

102:      class PointSequence {
103:      public:
104:        typedef ALE::SifterDef::IndexSequenceTraits<typename ::boost::multi_index::index<set_type, pointTag>::type,
105:                                     BOOST_MULTI_INDEX_MEMBER(rec_type, typename rec_type::point_type,point)>
106:        traits;
107:       protected:
108:         const typename traits::index_type& _index;
109:       public:
110: 
111:        // Need to extend the inherited iterator to be able to extract the degree
112:        class iterator : public traits::iterator {
113:        public:
114:          iterator(const typename traits::iterator::itor_type& itor) : traits::iterator(itor) {};
115:          virtual const int& degree() const {return this->_itor->degree;};
116:          virtual const int& marker() const {return this->_itor->marker;};
117:          virtual const int& depth()  const {return this->_itor->depth;};
118:          virtual const int& height() const {return this->_itor->height;};
119:          //void setDegree(const int degree) {this->_index.modify(this->_itor, typename traits::rec_type::degreeAdjuster(degree));};
120:        };
121: 
122:        PointSequence(const PointSequence& seq)            : _index(seq._index) {};
123:        PointSequence(typename traits::index_type& index) : _index(index)     {};
124:        virtual ~PointSequence(){};
125: 
126:        virtual bool empty(){return this->_index.empty();};
127: 
128:        virtual typename traits::index_type::size_type size() {return this->_index.size();};

130:        virtual iterator begin() {
131:          // Retrieve the beginning iterator of the index
132:          return iterator(this->_index.begin());
133:        };
134:        virtual iterator end() {
135:          // Retrieve the ending iterator of the index
136:          // Since the elements in this index are ordered by degree, this amounts to the end() of the index.
137:          return iterator(this->_index.end());
138:        };
139:        virtual bool contains(const typename rec_type::point_type& p) {
140:          // Check whether a given point is in the index
141:          return (this->_index.find(p) != this->_index.end());
142:        };
143:      }; // class PointSequence

145:      template<typename Tag_, typename Value_>
146:      class ValueSequence {
147:      public:
148:        typedef Value_ value_type;
149:        typedef ALE::SifterDef::IndexSequenceTraits<typename ::boost::multi_index::index<set_type, Tag_>::type,
150:                                    BOOST_MULTI_INDEX_MEMBER(rec_type, typename rec_type::point_type,point)>
151:        traits;
152:      protected:
153:        const typename traits::index_type& _index;
154:        const value_type _value;
155:      public:
156:        // Need to extend the inherited iterator to be able to extract the degree
157:        class iterator : public traits::iterator {
158:        public:
159:          iterator(const typename traits::iterator::itor_type& itor) : traits::iterator(itor) {};
160:          virtual const int& degree()  const {return this->_itor->degree;};
161:          virtual const int& marker()  const {return this->_itor->marker;};
162:          virtual const int& depth()   const {return this->_itor->depth;};
163:          virtual const int& height()  const {return this->_itor->height;};
164:        };
165: 
166:        ValueSequence(const PointSequence& seq)           : _index(seq._index), _value(seq._value) {};
167:        ValueSequence(typename traits::index_type& index, const value_type& value) : _index(index), _value(value) {};
168:        virtual ~ValueSequence(){};
169: 
170:        virtual bool empty(){return this->_index.empty();};
171: 
172:        virtual typename traits::index_type::size_type size() {return this->_index.count(this->_value);};

174:        virtual iterator begin() {
175:          return iterator(this->_index.lower_bound(this->_value));
176:        };
177:        virtual iterator end() {
178:          return iterator(this->_index.upper_bound(this->_value));
179:        };
180:      }; // class ValueSequence

182:      template<typename Tag_, typename Value_>
183:      class TwoValueSequence {
184:      public:
185:        typedef Value_ value_type;
186:        typedef ALE::SifterDef::IndexSequenceTraits<typename ::boost::multi_index::index<set_type, Tag_>::type,
187:                                    BOOST_MULTI_INDEX_MEMBER(rec_type, typename rec_type::point_type,point)>
188:        traits;
189:      protected:
190:        const typename traits::index_type& _index;
191:        const value_type _valueA, _valueB;
192:        const bool _useTwoValues;
193:      public:
194:        // Need to extend the inherited iterator to be able to extract the degree
195:        class iterator : public traits::iterator {
196:        public:
197:          iterator(const typename traits::iterator::itor_type& itor) : traits::iterator(itor) {};
198:          virtual const int& degree()  const {return this->_itor->degree;};
199:          virtual const int& marker()  const {return this->_itor->marker;};
200:        };
201: 
202:        TwoValueSequence(const PointSequence& seq)           : _index(seq._index), _valueA(seq._valueA), _valueB(seq._valueB), _useTwoValues(seq._useTwoValues) {};
203:        TwoValueSequence(typename traits::index_type& index, const value_type& valueA) : _index(index), _valueA(valueA), _valueB(value_type()), _useTwoValues(false) {};
204:        TwoValueSequence(typename traits::index_type& index, const value_type& valueA, const value_type& valueB) : _index(index), _valueA(valueA), _valueB(valueB), _useTwoValues(true) {};
205:        virtual ~TwoValueSequence(){};
206: 
207:        virtual bool empty(){return this->_index.empty();};
208: 
209:        virtual typename traits::index_type::size_type size() {
210:          if (this->_useTwoValues) {
211:            return this->_index.count(::boost::make_tuple(this->_valueA,this->_valueB));
212:          } else {
213:            return this->_index.count(::boost::make_tuple(this->_valueA));
214:          }
215:        };

217:        virtual iterator begin() {
218:          if (this->_useTwoValues) {
219:            return iterator(this->_index.lower_bound(::boost::make_tuple(this->_valueA,this->_valueB)));
220:          } else {
221:            return iterator(this->_index.lower_bound(::boost::make_tuple(this->_valueA)));
222:          }
223:        };
224:        virtual iterator end() {
225:          if (this->_useTwoValues) {
226:            return iterator(this->_index.upper_bound(::boost::make_tuple(this->_valueA,this->_valueB)));
227:          } else {
228:            return iterator(this->_index.upper_bound(::boost::make_tuple(this->_valueA)));
229:          }
230:        };
231:      }; // class TwoValueSequence
232:     };// struct RecContainerTraits

234:     template <typename Point_, typename Rec_>
235:     struct RecContainer {
236:       typedef RecContainerTraits<Point_, Rec_> traits;
237:       typedef typename traits::set_type set_type;
238:       template <typename OtherPoint_, typename OtherRec_>
239:       struct rebind {
240:         typedef RecContainer<OtherPoint_, OtherRec_> type;
241:       };
242:       set_type set;

244:       void removePoint(const typename traits::rec_type::point_type& p) {
245:         typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type& index =
246:           ::boost::multi_index::get<typename traits::pointTag>(this->set);
247:         typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type::iterator i = index.find(p);
248:         if (i != index.end()) { // Point exists
249:           index.erase(i);
250:         }
251:       };
252:       void adjustDegree(const typename traits::rec_type::point_type& p, int delta) {
253:         typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type& index =
254:           ::boost::multi_index::get<typename traits::pointTag>(this->set);
255:         typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type::iterator i = index.find(p);
256:         if (i == index.end()) { // No such point exists
257:           if(delta < 0) { // Cannot decrease degree of a non-existent point
258:             ostringstream err;
259:             err << "ERROR: adjustDegree: Non-existent point " << p;
260:             std::cout << err << std::endl;
261:             throw(Exception(err.str().c_str()));
262:           }
263:           else { // We CAN INCREASE the degree of a non-existent point: simply insert a new element with degree == delta
264:             std::pair<typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type::iterator, bool> ii;
265:             typename traits::rec_type r(p,delta);
266:             ii = index.insert(r);
267:             if(ii.second == false) {
268:               ostringstream err;
269:               err << "ERROR: adjustDegree: Failed to insert a rec " << r;
270:               std::cout << err << std::endl;
271:               throw(Exception(err.str().c_str()));
272:             }
273:           }
274:         }
275:         else { // Point exists, so we try to modify its degree
276:           // If the adjustment is zero, there is nothing to do, otherwise ...
277:           if(delta != 0) {
278:             int newDegree = i->degree + delta;
279:             if(newDegree < 0) {
280:               ostringstream ss;
281:               ss << "adjustDegree: Adjustment of " << *i << " by " << delta << " would result in negative degree: " << newDegree;
282:               throw Exception(ss.str().c_str());
283:             }
284:             index.modify(i, typename traits::rec_type::degreeAdjuster(newDegree));
285:           }
286:         }
287:       }; // adjustDegree()
288:     }; // struct RecContainer
289:   };// namespace SieveDef

291:     //
292:     // Sieve:
293:     //   A Sieve is a set of {\emph arrows} connecting {\emph points} of type Point_. Thus we
294:     // could realize a sieve, for instance, as a directed graph. In addition, we will often
295:     // assume an acyclicity constraint, arriving at a DAG. Each arrow may also carry a label,
296:     // or {\emph color} of type Color_, and the interface allows sets of arrows to be filtered
297:     // by color.
298:     //
299:     template <typename Point_, typename Marker_, typename Color_>
300:     class Sieve : public ALE::Sifter<Point_, Point_, Color_, ::boost::multi_index::composite_key_compare<std::less<Point_>, std::less<Color_>, std::less<Point_> >, SieveDef::RecContainer<Point_, SieveDef::Rec<Point_, Marker_> >, SieveDef::RecContainer<Point_, SieveDef::Rec<Point_, Marker_> > > {
301:     public:
302:       typedef Color_  color_type;
303:       typedef Point_  point_type;
304:       typedef Marker_ marker_type;
305:       typedef struct {
306:         typedef ALE::Sifter<Point_, Point_, Color_, ::boost::multi_index::composite_key_compare<std::less<Point_>, std::less<Color_>, std::less<Point_> >, SieveDef::RecContainer<Point_, SieveDef::Rec<Point_, Marker_> >, SieveDef::RecContainer<Point_, SieveDef::Rec<Point_, Marker_> > > baseType;
307:         // Encapsulated container types
308:         typedef typename baseType::traits::arrow_container_type arrow_container_type;
309:         typedef typename baseType::traits::cap_container_type   cap_container_type;
310:         typedef typename baseType::traits::base_container_type  base_container_type;
311:         // Types associated with records held in containers
312:         typedef typename baseType::traits::arrow_type           arrow_type;
313:         typedef typename baseType::traits::source_type          source_type;
314:         typedef typename baseType::traits::sourceRec_type       sourceRec_type;
315:         typedef typename baseType::traits::target_type          target_type;
316:         typedef typename baseType::traits::targetRec_type       targetRec_type;
317:         typedef typename baseType::traits::color_type           color_type;
318:         typedef Point_                                          point_type;
319:         // Convenient tag names
320:         typedef typename baseType::traits::supportInd           supportInd;
321:         typedef typename baseType::traits::coneInd              coneInd;
322:         typedef typename baseType::traits::arrowInd             arrowInd;
323:         typedef typename baseType::traits::baseInd              baseInd;
324:         typedef typename baseType::traits::capInd               capInd;

326:         //
327:         // Return types
328:         //
329:         typedef typename baseType::traits::arrowSequence        arrowSequence;
330:         typedef typename baseType::traits::coneSequence         coneSequence;
331:         typedef typename baseType::traits::supportSequence      supportSequence;
332:         typedef typename baseType::traits::baseSequence         baseSequence;
333:         typedef typename baseType::traits::capSequence          capSequence;
334:         typedef typename base_container_type::traits::template TwoValueSequence<typename base_container_type::traits::depthMarkerTag,int> depthSequence;
335:         typedef typename cap_container_type::traits::template TwoValueSequence<typename cap_container_type::traits::heightMarkerTag,int> heightSequence;
336:         typedef typename cap_container_type::traits::template ValueSequence<typename cap_container_type::traits::markerTag,marker_type> markerSequence;
337:       } traits;
338:       typedef std::set<point_type>    pointSet;
339:       typedef ALE::array<point_type>  pointArray;
340:       typedef std::set<marker_type>   markerSet;
341:       typedef pointSet                coneSet;
342:       typedef pointSet                supportSet;
343:       typedef pointArray              coneArray;
344:       typedef pointArray              supportArray;
345:     public:
346:       Sieve(MPI_Comm comm = PETSC_COMM_SELF, const int& debug = 0) : ALE::Sifter<Point_, Point_, Color_, ::boost::multi_index::composite_key_compare<std::less<Point_>, std::less<Color_>, std::less<Point_> >, SieveDef::RecContainer<Point_, SieveDef::Rec<Point_, Marker_> >, SieveDef::RecContainer<Point_, SieveDef::Rec<Point_, Marker_> > >(comm, debug), doStratify(false), maxDepth(-1), maxHeight(-1), graphDiameter(-1) {
347:         this->_markers = new markerSet();
348:         this->_meetSet = new coneSet();
349:       };
350:       virtual ~Sieve() {};
351:       // Printing
352:       friend std::ostream& operator<<(std::ostream& os, Obj<Sieve<Point_,Marker_,Color_> > s) {
353:         os << *s;
354:         return os;
355:       };
356: 
357:       friend std::ostream& operator<<(std::ostream& os, Sieve<Point_,Marker_,Color_>& s) {
358:         Obj<typename traits::baseSequence> base = s.base();
359:         for(typename traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
360:           Obj<typename traits::coneSequence> cone = s.cone(*b_iter);
361:           os << "Base point " << *b_iter << " with cone:" << std::endl;
362:           for(typename traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
363:             os << "  " << *c_iter << std::endl;
364:           }
365:         }
366:         return os;
367:       };

369:       template<typename ostream_type>
370:       void view(ostream_type& os, const char* label, bool rawData);
371:       void view(const char* label, MPI_Comm comm);

373:       Obj<Sieve> copy() {
374:         Obj<Sieve> s = Sieve(this->comm(), this->debug);
375:         Obj<typename traits::capSequence>  cap  = this->cap();
376:         Obj<typename traits::baseSequence> base = this->base();

378:         for(typename traits::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
379:           s->addCapPoint(*c_iter);
380:         }
381:         for(typename traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
382:           Obj<typename traits::coneSequence> cone = this->cone(*b_iter);

384:           for(typename traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
385:             s->addArrow(*c_iter, *b_iter, c_iter.color());
386:           }
387:         }
388:         s->stratify();
389:         return s;
390:       };
391:       bool hasPoint(const point_type& point) {
392:         if (this->baseContains(point) || this->capContains(point)) return true;
393:         return false;
394:       };
395:     private:
396:       template<class InputSequence> Obj<coneSet> __nCone(Obj<InputSequence>& cone, int n, const Color_& color, bool useColor);
397:       template<class pointSequence> void __nCone(const Obj<pointSequence>& points, int n, const Color_& color, bool useColor, Obj<coneArray> cone, Obj<coneSet> seen);
398:       template<class pointSequence> void __nSupport(const Obj<pointSequence>& points, int n, const Color_& color, bool useColor, Obj<supportArray> cone, Obj<supportSet> seen);
399:     public:
400:       //
401:       // The basic Sieve interface (extensions to Sifter)
402:       //
403:       Obj<coneArray> nCone(const Point_& p, int n);
404:       Obj<coneArray> nCone(const Point_& p, int n, const Color_& color, bool useColor = true);
405:       template<class InputSequence> Obj<coneSet> nCone(const Obj<InputSequence>& points, int n);
406:       template<class InputSequence> Obj<coneSet> nCone(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor = true);

408:       Obj<supportArray> nSupport(const Point_& p, int n);
409:       Obj<supportArray> nSupport(const Point_& p, int n, const Color_& color, bool useColor = true);
410:       template<class InputSequence> Obj<supportSet> nSupport(const Obj<InputSequence>& points, int n);
411:       template<class InputSequence> Obj<supportSet> nSupport(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor = true);
412:     public:
413:       virtual bool checkArrow(const typename traits::arrow_type& a) {
414:         if ((this->_cap.set.find(a.source) == this->_cap.set.end()) && (this->_base.set.find(a.source) == this->_base.set.end())) return false;
415:         if ((this->_cap.set.find(a.target) == this->_cap.set.end()) && (this->_base.set.find(a.target) == this->_base.set.end())) return false;
416:         return true;
417:       };
418:       //
419:       // Iterated versions
420:       //
421:       Obj<coneSet> closure(const Point_& p);

423:       Obj<coneSet> closure(const Point_& p, const Color_& color);

425:       template<class InputSequence>
426:       Obj<coneSet> closure(const Obj<InputSequence>& points);

428:       template<class InputSequence>
429:       Obj<coneSet> closure(const Obj<InputSequence>& points, const Color_& color);

431:       Obj<coneSet> nClosure(const Point_& p, int n);

433:       Obj<coneSet> nClosure(const Point_& p, int n, const Color_& color, bool useColor = true);

435:       template<class InputSequence>
436:       Obj<coneSet> nClosure(const Obj<InputSequence>& points, int n);

438:       template<class InputSequence>
439:       Obj<coneSet> nClosure(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor = true);

441:       Obj<Sieve<Point_,Marker_,Color_> > closureSieve(const Point_& p);

443:       Obj<Sieve<Point_,Marker_,Color_> > closureSieve(const Point_& p, const Color_& color);

445:       template<class InputSequence>
446:       Obj<Sieve<Point_,Marker_,Color_> > closureSieve(const Obj<InputSequence>& points);

448:       template<class InputSequence>
449:       Obj<Sieve<Point_,Marker_,Color_> > closureSieve(const Obj<InputSequence>& points, const Color_& color);

451:       Obj<Sieve<Point_,Marker_,Color_> > nClosureSieve(const Point_& p, int n);

453:       Obj<Sieve<Point_,Marker_,Color_> > nClosureSieve(const Point_& p, int n, const Color_& color, bool useColor = true);

455:       template<class InputSequence>
456:       Obj<Sieve<Point_,Marker_,Color_> > nClosureSieve(const Obj<InputSequence>& points, int n);

458:       template<class InputSequence>
459:       Obj<Sieve<Point_,Marker_,Color_> > nClosureSieve(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor = true);

461:       Obj<supportSet> star(const Point_& p);

463:       Obj<supportSet> star(const Point_& p, const Color_& color);

465:       template<class InputSequence>
466:       Obj<supportSet> star(const Obj<InputSequence>& points);

468:       template<class InputSequence>
469:       Obj<supportSet> star(const Obj<InputSequence>& points, const Color_& color);

471:       Obj<supportSet> nStar(const Point_& p, int n);

473:       Obj<supportSet> nStar(const Point_& p, int n, const Color_& color, bool useColor = true);

475:       template<class InputSequence>
476:       Obj<supportSet> nStar(const Obj<InputSequence>& points, int n);

478:       template<class InputSequence>
479:       Obj<supportSet> nStar(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor = true);

481:     private:
482:       template<class InputSequence>
483:       Obj<coneSet> __nClosure(Obj<InputSequence>& cone, int n, const Color_& color, bool useColor);

485:       template<class InputSequence>
486:       Obj<Sieve<Point_,Marker_,Color_> > __nClosureSieve(Obj<InputSequence>& cone, int n, const Color_& color, bool useColor);

488:       template<class InputSequence>
489:       Obj<supportSet> __nStar(Obj<InputSequence>& support, int n, const Color_& color, bool useColor);

491:     public:
492:       //
493:       // Lattice methods
494:       //
495:       const Obj<coneSet>& meet(const Point_& p, const Point_& q);

497:       const Obj<coneSet>& meet(const Point_& p, const Point_& q, const Color_& color);

499:       template<class InputSequence>
500:       const Obj<coneSet>& meet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1);

502:       template<class InputSequence>
503:       const Obj<coneSet>& meet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, const Color_& color);

505:       const Obj<coneSet>& nMeet(const Point_& p, const Point_& q, int n);

507:       const Obj<coneSet>& nMeet(const Point_& p, const Point_& q, int n, const Color_& color, bool useColor = true);

509:       template<class InputSequence>
510:       const Obj<coneSet>& nMeet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n);

512:       template<class InputSequence>
513:       const Obj<coneSet>& nMeet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n,
514:                                 const Color_& color, bool useColor = true);

516:       Obj<supportSet> join(const Point_& p, const Point_& q);

518:       Obj<supportSet> join(const Point_& p, const Point_& q, const Color_& color);

520:       template<class InputSequence>
521:       Obj<supportSet> join(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1);

523:       template<class InputSequence>
524:       Obj<supportSet> join(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, const Color_& color);

526:       Obj<supportSet> nJoin(const Point_& p, const Point_& q, int n);

528:       Obj<supportSet> nJoin(const Point_& p, const Point_& q, int n, const Color_& color, bool useColor = true);

530:       template<class InputSequence>
531:       Obj<supportSet> nJoin(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n);

533:       template<class InputSequence>
534:       Obj<supportSet> nJoin(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n, const Color_& color, bool useColor = true);

536:     public:
537:       Obj<typename traits::depthSequence> roots() {
538:         return this->depthStratum(0);
539:       };
540:       Obj<typename traits::heightSequence> leaves() {
541:         return this->heightStratum(0);
542:       };
543:     private:
544:       bool doStratify;
545:       int  maxDepth, maxHeight, graphDiameter;
546:     public:
547:       //
548:       // Structural queries
549:       //
550:       int depth();
551:       int depth(const point_type& p);
552:       template<typename InputSequence> int depth(const Obj<InputSequence>& points);

554:       int height();
555:       int height(const point_type& p);
556:       template<typename InputSequence> int height(const Obj<InputSequence>& points);

558:       int diameter();
559:       int diameter(const point_type& p);

561:       Obj<typename traits::depthSequence> depthStratum(int d);
562:       Obj<typename traits::depthSequence> depthStratum(int d, marker_type m);

564:       Obj<typename traits::heightSequence> heightStratum(int h);
565:       Obj<typename traits::heightSequence> heightStratum(int h, marker_type m);

567:       Obj<typename traits::markerSequence> markerStratum(marker_type m);
568: 
569:       void setStratification(bool doStratify) {this->doStratify = doStratify;};

571:       bool getStratification() {return this->doStratify;};

573:       void stratify(bool show = false);
574:     protected:
575:       Obj<markerSet> _markers;
576:       Obj<coneSet>   _meetSet;
577:     public:
578:       //
579:       // Structural manipulation
580:       //

582:       struct changeMarker {
583:         changeMarker(int newMarker) : newMarker(newMarker) {};

585:         void operator()(typename traits::base_container_type::traits::rec_type& p) {
586:           p.marker = newMarker;
587:         }
588:       private:
589:         marker_type newMarker;
590:       };

592:       void setMarker(const point_type& p, const marker_type& marker);
593:       template<class InputSequence> void setMarker(const Obj<InputSequence>& points, const marker_type& marker);

595:       void clearMarkers() {this->_markers.clear();};
596:       Obj<markerSet> markers() {return this->_markers;};
597:     private:
598:       struct changeHeight {
599:         changeHeight(int newHeight) : newHeight(newHeight) {};

601:         void operator()(typename traits::base_container_type::traits::rec_type& p) {
602:           p.height = newHeight;
603:         }
604:       private:
605:         int newHeight;
606:       };
607: 
608:       template<class InputSequence> void __computeClosureHeights(const Obj<InputSequence>& points);

610:       struct changeDepth {
611:         changeDepth(int newDepth) : newDepth(newDepth) {};

613:         void operator()(typename traits::base_container_type::traits::rec_type& p) {
614:           p.depth = newDepth;
615:         }
616:       private:
617:         int newDepth;
618:       };

620:       template<class InputSequence> void __computeStarDepths(const Obj<InputSequence>& points);
621:     };

623:     template <typename Point_, typename Marker_, typename Color_>
624:     Obj<typename Sieve<Point_,Marker_,Color_>::coneArray> Sieve<Point_,Marker_,Color_>::nCone(const Point_& p, int n) {
625:       return this->nCone(p, n, Color_(), false);
626:     };

628:     template <typename Point_, typename Marker_, typename Color_>
629:     Obj<typename Sieve<Point_,Marker_,Color_>::coneArray> Sieve<Point_,Marker_,Color_>::nCone(const Point_& p, int n, const Color_& color, bool useColor) {
630:       Obj<coneArray> cone = new coneArray();
631:       Obj<coneSet>   seen = new coneSet();

633:       this->__nCone(this->cone(p), n-1, color, useColor, cone, seen);
634:       return cone;
635:     };

637:     template <typename Point_, typename Marker_, typename Color_>
638:     template<class pointSequence>
639:     void Sieve<Point_,Marker_,Color_>::__nCone(const Obj<pointSequence>& points, int n, const Color_& color, bool useColor, Obj<coneArray> cone, Obj<coneSet> seen) {
640:       if (n == 0) {
641:         for(typename pointSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
642:           if (seen->find(*p_itor) == seen->end()) {
643:             cone->push_back(*p_itor);
644:             seen->insert(*p_itor);
645:           }
646:         }
647:       } else {
648:         typename pointSequence::iterator end = points->end();
649:         for(typename pointSequence::iterator p_itor = points->begin(); p_itor != end; ++p_itor) {
650:           if (useColor) {
651:             this->__nCone(this->cone(*p_itor, color), n-1, color, useColor, cone, seen);
652:           } else {
653:             this->__nCone(this->cone(*p_itor), n-1, color, useColor, cone, seen);
654:           }
655:         }
656:       }
657:     };

659:     template <typename Point_, typename Marker_, typename Color_>
660:     template<class InputSequence>
661:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::nCone(const Obj<InputSequence>& points, int n) {
662:       return this->nCone(points, n, Color_(), false);
663:     };

665:     template <typename Point_, typename Marker_, typename Color_>
666:     template<class InputSequence>
667:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::nCone(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor ) {
668:       Obj<coneSet> cone = new coneSet();
669:       cone->insert(points->begin(), points->end());
670:       return this->__nCone(cone, n, color, useColor);
671:     };

673:     template <typename Point_, typename Marker_, typename Color_>
674:     template<class InputSequence>
675:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::__nCone(Obj<InputSequence>& cone, int n, const Color_& color, bool useColor) {
676:       Obj<coneSet> base = new coneSet();

678:       for(int i = 0; i < n; ++i) {
679:         Obj<coneSet> tmp = cone; cone = base; base = tmp;
680: 
681:         cone->clear();
682:         for(typename coneSet::iterator b_itor = base->begin(); b_itor != base->end(); ++b_itor) {
683:           Obj<typename traits::coneSequence> pCone;
684: 
685:           if (useColor) {
686:             pCone = this->cone(*b_itor, color);
687:           } else {
688:             pCone = this->cone(*b_itor);
689:           }
690:           cone->insert(pCone->begin(), pCone->end());
691:         }
692:       }
693:       return cone;
694:     };

696:     template <typename Point_, typename Marker_, typename Color_>
697:     Obj<typename Sieve<Point_,Marker_,Color_>::supportArray> Sieve<Point_,Marker_,Color_>::nSupport(const Point_& p, int n) {
698:       return this->nSupport(p, n, Color_(), false);
699:     };

701:     template <typename Point_, typename Marker_, typename Color_>
702:     Obj<typename Sieve<Point_,Marker_,Color_>::supportArray> Sieve<Point_,Marker_,Color_>::nSupport(const Point_& p, int n, const Color_& color, bool useColor) {
703:       Obj<supportArray> cone = new supportArray();
704:       Obj<supportSet>   seen = new supportSet();

706:       this->__nSupport(this->support(p), n-1, color, useColor, cone, seen);
707:       return cone;
708:     };

710:     template <typename Point_, typename Marker_, typename Color_>
711:     template<class pointSequence>
712:     void Sieve<Point_,Marker_,Color_>::__nSupport(const Obj<pointSequence>& points, int n, const Color_& color, bool useColor, Obj<supportArray> support, Obj<supportSet> seen) {
713:       if (n == 0) {
714:         for(typename pointSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
715:           if (seen->find(*p_itor) == seen->end()) {
716:             support->push_back(*p_itor);
717:             seen->insert(*p_itor);
718:           }
719:         }
720:       } else {
721:         typename pointSequence::iterator end = points->end();
722:         for(typename pointSequence::iterator p_itor = points->begin(); p_itor != end; ++p_itor) {
723:           if (useColor) {
724:             this->__nSupport(this->support(*p_itor, color), n-1, color, useColor, support, seen);
725:           } else {
726:             this->__nSupport(this->support(*p_itor), n-1, color, useColor, support, seen);
727:           }
728:         }
729:       }
730:     };

732:     template <typename Point_, typename Marker_, typename Color_>
733:     template<class InputSequence>
734:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nSupport(const Obj<InputSequence>& points, int n) {
735:       return this->nSupport(points, n, Color_(), false);
736:     };

738:     template <typename Point_, typename Marker_, typename Color_>
739:     template<class InputSequence>
740:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nSupport(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor ) {
741:       Obj<supportSet> support = supportSet();
742:       Obj<supportSet> cap = supportSet();
743: 
744:       support->insert(points->begin(), points->end());
745:       for(int i = 0; i < n; ++i) {
746:         Obj<supportSet> tmp = support; support = cap; cap = tmp;
747: 
748:         support->clear();
749:         for(typename supportSet::iterator c_itor = cap->begin(); c_itor != cap->end(); ++c_itor) {
750:           Obj<typename traits::supportSequence> pSupport;
751: 
752:           if (useColor) {
753:             pSupport = this->support(*c_itor, color);
754:           } else {
755:             pSupport = this->support(*c_itor);
756:           }
757:           support->insert(pSupport->begin(), pSupport->end());
758:         }
759:       }
760:       return support;
761:     };
762:     //
763:     // Iterated versions
764:     //
765: 
766:     template <typename Point_, typename Marker_, typename Color_>
767:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::closure(const Point_& p) {
768:       return nClosure(p, this->depth());
769:     };
770: 
771:     template <typename Point_, typename Marker_, typename Color_>
772:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::closure(const Point_& p, const Color_& color) {
773:       return nClosure(p, this->depth(), color);
774:     };

776:     template <typename Point_, typename Marker_, typename Color_>
777:     template<class InputSequence>
778:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::closure(const Obj<InputSequence>& points) {
779:       return nClosure(points, this->depth());
780:     };
781: 
782:     template <typename Point_, typename Marker_, typename Color_>
783:     template<class InputSequence>
784:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::closure(const Obj<InputSequence>& points, const Color_& color) {
785:       return nClosure(points, this->depth(), color);
786:     };
787: 
788:     template <typename Point_, typename Marker_, typename Color_>
789:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::nClosure(const Point_& p, int n) {
790:       return this->nClosure(p, n, Color_(), false);
791:     };

793:     template <typename Point_, typename Marker_, typename Color_>
794:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::nClosure(const Point_& p, int n, const Color_& color, bool useColor ) {
795:       Obj<coneSet> cone = coneSet();
796:       cone->insert(p);
797:       return this->__nClosure(cone, n, color, useColor);
798:     };

800:     template <typename Point_, typename Marker_, typename Color_>
801:     template<class InputSequence>
802:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::nClosure(const Obj<InputSequence>& points, int n) {
803:       return this->nClosure(points, n, Color_(), false);
804:     };

806:     template <typename Point_, typename Marker_, typename Color_>
807:     template<class InputSequence>
808:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::nClosure(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor ) {
809:       Obj<coneSet> cone = coneSet();
810:       cone->insert(points->begin(), points->end());
811:       return this->__nClosure(cone, n, color, useColor);
812:     };

814:     template <typename Point_, typename Marker_, typename Color_>
815:     template<class InputSequence>
816:     Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::__nClosure(Obj<InputSequence>& cone, int n, const Color_& color, bool useColor) {
817:       Obj<coneSet> base    = coneSet();
818:       Obj<coneSet> closure = coneSet();
819: 
820:       closure->insert(cone->begin(), cone->end());
821:       for(int i = 0; i < n; ++i) {
822:         Obj<coneSet> tmp = cone; cone = base; base = tmp;
823:         cone->clear();
824:         for(typename coneSet::iterator b_itor = base->begin(); b_itor != base->end(); ++b_itor) {
825:           Obj<typename traits::coneSequence> pCone;
826:           if (useColor) {
827:             pCone = this->cone(*b_itor, color);
828:           } else {
829:             pCone = this->cone(*b_itor);
830:           }
831:           cone->insert(pCone->begin(), pCone->end());
832:           closure->insert(pCone->begin(), pCone->end());
833:         }
834:       }
835:       return closure;
836:     };

838:     template <typename Point_, typename Marker_, typename Color_>
839:     Obj<Sieve<Point_,Marker_,Color_> > Sieve<Point_,Marker_,Color_>::closureSieve(const Point_& p) {
840:       return nClosureSieve(p, this->depth());
841:     };

843:     template <typename Point_, typename Marker_, typename Color_>
844:     Obj<Sieve<Point_,Marker_,Color_> > Sieve<Point_,Marker_,Color_>::closureSieve(const Point_& p, const Color_& color) {
845:       return nClosureSieve(p, this->depth(), color);
846:     };
847: 
848:     template <typename Point_, typename Marker_, typename Color_>
849:     template<class InputSequence>
850:     Obj<Sieve<Point_,Marker_,Color_> > Sieve<Point_,Marker_,Color_>::closureSieve(const Obj<InputSequence>& points) {
851:       return nClosureSieve(points, this->depth());
852:     };
853: 
854:     template <typename Point_, typename Marker_, typename Color_>
855:     template<class InputSequence>
856:     Obj<Sieve<Point_,Marker_,Color_> > Sieve<Point_,Marker_,Color_>::closureSieve(const Obj<InputSequence>& points, const Color_& color) {
857:       return nClosureSieve(points, this->depth(), color);
858:     };
859: 
860:     template <typename Point_, typename Marker_, typename Color_>
861:     Obj<Sieve<Point_,Marker_,Color_> > Sieve<Point_,Marker_,Color_>::nClosureSieve(const Point_& p, int n) {
862:       return this->nClosureSieve(p, n, Color_(), false);
863:     };
864: 
865:     template <typename Point_, typename Marker_, typename Color_>
866:     Obj<Sieve<Point_,Marker_,Color_> > Sieve<Point_,Marker_,Color_>::nClosureSieve(const Point_& p, int n, const Color_& color, bool useColor ) {
867:       Obj<coneSet> cone = coneSet();
868:       cone->insert(p);
869:       return this->__nClosureSieve(cone, n, color, useColor);
870:     };
871: 
872:     template <typename Point_, typename Marker_, typename Color_>
873:     template<class InputSequence>
874:     Obj<Sieve<Point_,Marker_,Color_> > Sieve<Point_,Marker_,Color_>::nClosureSieve(const Obj<InputSequence>& points, int n) {
875:       return this->nClosureSieve(points, n, Color_(), false);
876:     };
877: 
878:     template <typename Point_, typename Marker_, typename Color_>
879:     template<class InputSequence>
880:     Obj<Sieve<Point_,Marker_,Color_> > Sieve<Point_,Marker_,Color_>::nClosureSieve(const Obj<InputSequence>& points, int n,
881:                                                                    const Color_& color, bool useColor ) {
882:       Obj<coneSet> cone = coneSet();
883:       cone->insert(points->begin(), points->end());
884:       return this->__nClosureSieve(cone, n, color, useColor);
885:     };
886: 
887:     template <typename Point_, typename Marker_, typename Color_>
888:     template<class InputSequence>
889:     Obj<Sieve<Point_,Marker_,Color_> > Sieve<Point_,Marker_,Color_>::__nClosureSieve(Obj<InputSequence>& cone,int n,const Color_& color,bool useColor) {
890:       Obj<coneSet> base = coneSet();
891:       Obj<Sieve<Point_,Marker_,Color_> > closure = Sieve<Point_,Marker_,Color_>();
892:       // FIX: need to add initial stuff
893:       for(int i = 0; i < n; ++i) {
894:         Obj<coneSet> tmp = cone; cone = base; base = tmp;
895:         cone->clear();
896:         for(typename coneSet::iterator b_itor = base->begin(); b_itor != base->end(); ++b_itor) {
897:           Obj<typename traits::coneSequence> pCone;
898: 
899:           if (useColor) {
900:             pCone = this->cone(*b_itor, color);
901:           } else {
902:             pCone = this->cone(*b_itor);
903:           }
904:           cone->insert(pCone->begin(), pCone->end());
905:           closure->addCone(pCone, *b_itor);
906:         }
907:       }
908:       return closure;
909:     };
910: 
911:     template <typename Point_, typename Marker_, typename Color_>
912:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::star(const Point_& p) {
913:       return nStar(p, this->height());
914:     };
915: 
916:     template <typename Point_, typename Marker_, typename Color_>
917:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::star(const Point_& p, const Color_& color) {
918:       return nStar(p, this->depth(), color);
919:     };
920: 
921:     template <typename Point_, typename Marker_, typename Color_>
922:     template<class InputSequence>
923:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::star(const Obj<InputSequence>& points) {
924:       return nStar(points, this->height());
925:     };
926: 
927:     template <typename Point_, typename Marker_, typename Color_>
928:     template<class InputSequence>
929:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::star(const Obj<InputSequence>& points, const Color_& color) {
930:       return nStar(points, this->height(), color);
931:     };

933:     template <typename Point_, typename Marker_, typename Color_>
934:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nStar(const Point_& p, int n) {
935:       return this->nStar(p, n, Color_(), false);
936:     };
937: 
938:     template <typename Point_, typename Marker_, typename Color_>
939:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nStar(const Point_& p, int n, const Color_& color, bool useColor ) {
940:       Obj<supportSet> support = supportSet();
941:       support->insert(p);
942:       return this->__nStar(support, n, color, useColor);
943:     };

945:     template <typename Point_, typename Marker_, typename Color_>
946:     template<class InputSequence>
947:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nStar(const Obj<InputSequence>& points, int n) {
948:       return this->nStar(points, n, Color_(), false);
949:     };

951:     template <typename Point_, typename Marker_, typename Color_>
952:     template<class InputSequence>
953:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nStar(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor ) {
954:       Obj<supportSet> support = supportSet();
955:       support->insert(points->begin(), points->end());
956:       return this->__nStar(support, n, color, useColor);
957:     };
958: 
959:     template <typename Point_, typename Marker_, typename Color_>
960:     template<class InputSequence>
961:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::__nStar(Obj<InputSequence>& support, int n, const Color_& color, bool useColor) {
962:       Obj<supportSet> cap = supportSet();
963:       Obj<supportSet> star = supportSet();
964:       star->insert(support->begin(), support->end());
965:       for(int i = 0; i < n; ++i) {
966:         Obj<supportSet> tmp = support; support = cap; cap = tmp;
967:         support->clear();
968:         for(typename supportSet::iterator c_itor = cap->begin(); c_itor != cap->end(); ++c_itor) {
969:           Obj<typename traits::supportSequence> pSupport;
970:           if (useColor) {
971:             pSupport = this->support(*c_itor, color);
972:           } else {
973:             pSupport = this->support(*c_itor);
974:           }
975:           support->insert(pSupport->begin(), pSupport->end());
976:           star->insert(pSupport->begin(), pSupport->end());
977:         }
978:       }
979:       return star;
980:     };

982:     //
983:     // Lattice methods
984:     //

986:     template <typename Point_, typename Marker_, typename Color_>
987:     const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::meet(const Point_& p, const Point_& q) {
988:       return nMeet(p, q, this->depth());
989:     };
990: 
991:     template <typename Point_, typename Marker_, typename Color_>
992:     const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::meet(const Point_& p, const Point_& q, const Color_& color) {
993:       return nMeet(p, q, this->depth(), color);
994:     };

996:     template <typename Point_, typename Marker_, typename Color_>
997:     template<class InputSequence>
998:     const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::meet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1) {
999:       return nMeet(chain0, chain1, this->depth());
1000:     };

1002:     template <typename Point_, typename Marker_, typename Color_>
1003:     template<class InputSequence>
1004:     const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::meet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, const Color_& color) {
1005:         return nMeet(chain0, chain1, this->depth(), color);
1006:     };

1008:     template <typename Point_, typename Marker_, typename Color_>
1009:     const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::nMeet(const Point_& p, const Point_& q, int n) {
1010:       if (n == 1) {
1011:         std::vector<point_type> vecA, vecB;
1012:         const Obj<typename traits::coneSequence>&     coneA  = this->cone(p);
1013:         const typename traits::coneSequence::iterator beginA = coneA->begin();
1014:         const typename traits::coneSequence::iterator endA   = coneA->end();
1015:         const Obj<typename traits::coneSequence>&     coneB  = this->cone(q);
1016:         const typename traits::coneSequence::iterator beginB = coneB->begin();
1017:         const typename traits::coneSequence::iterator endB   = coneB->end();

1019:         vecA.insert(vecA.begin(), beginA, endA);
1020:         std::sort(vecA.begin(), vecA.end());
1021:         vecB.insert(vecB.begin(), beginB, endB);
1022:         std::sort(vecB.begin(), vecB.end());
1023:         this->_meetSet->clear();
1024:         std::set_intersection(vecA.begin(), vecA.end(), vecB.begin(), vecB.end(), std::insert_iterator<typename Sieve<Point_,Marker_,Color_>::coneSet>(*this->_meetSet, this->_meetSet->begin()));
1025:         return this->_meetSet;
1026:       }
1027:       return nMeet(p, q, n, Color_(), false);
1028:     };

1030:     template <typename Point_, typename Marker_, typename Color_>
1031:     const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::nMeet(const Point_& p, const Point_& q, int n, const Color_& color, bool useColor ) {
1032:       Obj<coneSet> chain0 = new coneSet();
1033:       Obj<coneSet> chain1 = new coneSet();
1034:       chain0->insert(p);
1035:       chain1->insert(q);
1036:       return this->nMeet(chain0, chain1, n, color, useColor);
1037:     };

1039:     template <typename Point_, typename Marker_, typename Color_>
1040:     template<class InputSequence>
1041:     const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::nMeet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n) {
1042:       return this->nMeet(chain0, chain1, n, Color_(), false);
1043:     };
1044: 
1045:     template <typename Point_, typename Marker_, typename Color_>
1046:     template<class InputSequence>
1047:     const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::nMeet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1,int n,const Color_& color, bool useColor){
1048:       // The strategy is to compute the intersection of cones over the chains, remove the intersection
1049:       // and use the remaining two parts -- two disjoined components of the symmetric difference of cones -- as the new chains.
1050:       // The intersections at each stage are accumulated and their union is the meet.
1051:       // The iteration stops after n steps in addition to the meet of the initial chains or sooner if at least one of the chains is empty.
1052:       Obj<coneSet> cone;

1054:       this->_meetSet->clear();
1055:       if((chain0->size() != 0) && (chain1->size() != 0)) {
1056:         for(int i = 0; i <= n; ++i) {
1057:           // Compute the intersection of chains and put it in meet at the same time removing it from c and cc
1058:           std::set<point_type> intersect;
1059:           //std::set_intersection(chain0->begin(), chain0->end(), chain1->begin(), chain1->end(), std::insert_iterator<coneSet>(meet, meet->begin()));
1060:           std::set_intersection(chain0->begin(), chain0->end(), chain1->begin(), chain1->end(), std::insert_iterator<std::set<point_type> >(intersect, intersect.begin()));
1061:           this->_meetSet->insert(intersect.begin(), intersect.end());
1062:           for(typename std::set<point_type>::iterator i_iter = intersect.begin(); i_iter != intersect.end(); ++i_iter) {
1063:             chain0->erase(chain0->find(*i_iter));
1064:             chain1->erase(chain1->find(*i_iter));
1065:           }
1066:           // Replace each of the cones with a cone over it, and check if either is empty; if so, return what's in meet at the moment.
1067:           cone = this->cone(chain0);
1068:           chain0->insert(cone->begin(), cone->end());
1069:           if(chain0->size() == 0) {
1070:             break;
1071:           }
1072:           cone = this->cone(chain1);
1073:           chain1->insert(cone->begin(), cone->end());
1074:           if(chain1->size() == 0) {
1075:             break;
1076:           }
1077:           // If both cones are empty, we should quit
1078:         }
1079:       }
1080:       return this->_meetSet;
1081:     };
1082: 
1083:     template <typename Point_, typename Marker_, typename Color_>
1084:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::join(const Point_& p, const Point_& q) {
1085:       return this->nJoin(p, q, this->depth());
1086:     };
1087: 
1088:     template <typename Point_, typename Marker_, typename Color_>
1089:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::join(const Point_& p, const Point_& q, const Color_& color) {
1090:       return this->nJoin(p, q, this->depth(), color);
1091:     };

1093:     template <typename Point_, typename Marker_, typename Color_>
1094:     template<class InputSequence>
1095:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::join(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1) {
1096:       return this->nJoin(chain0, chain1, this->depth());
1097:     };
1098: 
1099:     template <typename Point_, typename Marker_, typename Color_>
1100:     template<class InputSequence>
1101:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::join(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, const Color_& color) {
1102:       return this->nJoin(chain0, chain1, this->depth(), color);
1103:     };
1104: 
1105:     template <typename Point_, typename Marker_, typename Color_>
1106:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nJoin(const Point_& p, const Point_& q, int n) {
1107:       return this->nJoin(p, q, n, Color_(), false);
1108:     };
1109: 
1110:     template <typename Point_, typename Marker_, typename Color_>
1111:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nJoin(const Point_& p, const Point_& q, int n, const Color_& color, bool useColor) {
1112:       Obj<supportSet> chain0 = supportSet();
1113:       Obj<supportSet> chain1 = supportSet();
1114:       chain0->insert(p);
1115:       chain1->insert(q);
1116:       return this->nJoin(chain0, chain1, n, color, useColor);
1117:     };

1119:     template <typename Point_, typename Marker_, typename Color_>
1120:     template<class InputSequence>
1121:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nJoin(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n) {
1122:       return this->nJoin(chain0, chain1, n, Color_(), false);
1123:     };

1125:     template <typename Point_, typename Marker_, typename Color_>
1126:     template<class InputSequence>
1127:     Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nJoin(const Obj<InputSequence>& chain0,const Obj<InputSequence>& chain1,int n,const Color_& color,bool useColor){
1128:       // The strategy is to compute the intersection of supports over the chains, remove the intersection
1129:       // and use the remaining two parts -- two disjoined components of the symmetric difference of supports -- as the new chains.
1130:       // The intersections at each stage are accumulated and their union is the join.
1131:       // The iteration stops after n steps in addition to the join of the initial chains or sooner if at least one of the chains is empty.
1132:       Obj<supportSet> join = supportSet();
1133:       Obj<supportSet> support;
1134:       if((chain0->size() != 0) && (chain1->size() != 0)) {
1135:         for(int i = 0; i <= n; ++i) {
1136:           // Compute the intersection of chains and put it in meet at the same time removing it from c and cc
1137:           std::set<point_type> intersect;
1138:           //std::set_intersection(chain0->begin(), chain0->end(), chain1->begin(), chain1->end(), std::insert_iterator<supportSet>(join.obj(), join->begin()));
1139:           std::set_intersection(chain0->begin(), chain0->end(), chain1->begin(), chain1->end(), std::insert_iterator<std::set<point_type> >(intersect, intersect.begin()));
1140:           join->insert(intersect.begin(), intersect.end());
1141:           for(typename std::set<point_type>::iterator i_iter = intersect.begin(); i_iter != intersect.end(); ++i_iter) {
1142:             chain0->erase(chain0->find(*i_iter));
1143:             chain1->erase(chain1->find(*i_iter));
1144:           }
1145:           // Replace each of the supports with the support over it, and check if either is empty; if so, return what's in join at the moment.
1146:           support = this->support(chain0);
1147:           chain0->insert(support->begin(), support->end());
1148:           if(chain0->size() == 0) {
1149:             break;
1150:           }
1151:           support = this->support(chain1);
1152:           chain1->insert(support->begin(), support->end());
1153:           if(chain1->size() == 0) {
1154:             break;
1155:           }
1156:           // If both supports are empty, we should quit
1157:         }
1158:       }
1159:       return join;
1160:     };

1162:     template <typename Point_, typename Marker_, typename Color_>
1163:     template<typename ostream_type>
1164:     void Sieve<Point_,Marker_,Color_>::view(ostream_type& os, const char* label = NULL, bool rawData = false){
1165:       if(label != NULL) {
1166:         os << "Viewing Sieve '" << label << "':" << std::endl;
1167:       }
1168:       else {
1169:         os << "Viewing a Sieve:" << std::endl;
1170:       }
1171:       if(!rawData) {
1172:         os << "cap --> base:" << std::endl;
1173:         Obj<typename traits::capSequence> cap = this->cap();
1174:         for(typename traits::capSequence::iterator capi = cap->begin(); capi != cap->end(); ++capi) {
1175:           Obj<typename traits::supportSequence> supp = this->support(*capi);
1176:           for(typename traits::supportSequence::iterator suppi = supp->begin(); suppi != supp->end(); ++suppi) {
1177:             os << *capi << "--(" << suppi.color() << ")-->" << *suppi << std::endl;
1178:           }
1179:         }
1180:         os << "base <-- cap:" << std::endl;
1181:         Obj<typename traits::baseSequence> base = this->base();
1182:         for(typename traits::baseSequence::iterator basei = base->begin(); basei != base->end(); ++basei) {
1183:           Obj<typename traits::coneSequence> cone = this->cone(*basei);
1184:           for(typename traits::coneSequence::iterator conei = cone->begin(); conei != cone->end(); ++conei) {
1185:             os << *basei <<  "<--(" << conei.color() << ")--" << *conei << std::endl;
1186:           }
1187:         }
1188:         os << "cap --> (outdegree, marker, depth, height):" << std::endl;
1189:         for(typename traits::capSequence::iterator capi = cap->begin(); capi != cap->end(); ++capi) {
1190:           os << *capi <<  "-->" << capi.degree() << ", " << capi.marker() << ", " << capi.depth() << ", " << capi.height() << std::endl;
1191:         }
1192:         os << "base --> (indegree, marker, depth, height):" << std::endl;
1193:         for(typename traits::baseSequence::iterator basei = base->begin(); basei != base->end(); ++basei) {
1194:           os << *basei <<  "-->" << basei.degree() << ", " << basei.marker() << ", " << basei.depth() << ", " << basei.height() << std::endl;
1195:         }
1196:       }
1197:       else {
1198:         os << "'raw' arrow set:" << std::endl;
1199:         for(typename traits::arrow_container_type::set_type::iterator ai = this->_arrows.set.begin(); ai != this->_arrows.set.end(); ai++)
1200:         {
1201:           typename traits::arrow_type arr = *ai;
1202:           os << arr << std::endl;
1203:         }
1204:         os << "'raw' base set:" << std::endl;
1205:         for(typename traits::base_container_type::set_type::iterator bi = this->_base.set.begin(); bi != this->_base.set.end(); bi++)
1206:         {
1207:           typename traits::base_container_type::traits::rec_type bp = *bi;
1208:           os << bp << std::endl;
1209:         }
1210:         os << "'raw' cap set:" << std::endl;
1211:         for(typename traits::cap_container_type::set_type::iterator ci = this->_cap.set.begin(); ci != this->_cap.set.end(); ci++)
1212:         {
1213:           typename traits::cap_container_type::traits::rec_type cp = *ci;
1214:           os << cp << std::endl;
1215:         }
1216:       }
1217:     };
1218:     template <typename Point_, typename Marker_, typename Color_>
1219:     void Sieve<Point_,Marker_,Color_>::view(const char* label = NULL, MPI_Comm comm = MPI_COMM_NULL) {
1220:         ostringstream txt;

1222:         if (this->debug()) {
1223:           std::cout << "viewing a Sieve, comm = " << this->comm() << ", commRank = " << this->commRank() << std::endl;
1224:         }
1225:         if(label != NULL) {
1226:           if(this->commRank() == 0) {
1227:             txt << "viewing Sieve :'" << label << "'" << std::endl;
1228:           }
1229:         }
1230:         else {
1231:           if(this->commRank() == 0) {
1232:             txt << "viewing a Sieve" << std::endl;
1233:           }
1234:         }
1235:         if(this->commRank() == 0) {
1236:           txt << "cap --> base:\n";
1237:         }
1238:         typename traits::capSequence cap   = this->cap();
1239:         typename traits::baseSequence base = this->base();
1240:         if(cap.empty()) {
1241:           txt << "[" << this->commRank() << "]: empty" << std::endl;
1242:         }
1243:         for(typename traits::capSequence::iterator capi = cap.begin(); capi != cap.end(); capi++) {
1244:           const Obj<typename traits::supportSequence>& supp = this->support(*capi);
1245:           const typename traits::supportSequence::iterator suppEnd = supp->end();

1247:           for(typename traits::supportSequence::iterator suppi = supp->begin(); suppi != suppEnd; suppi++) {
1248:             txt << "[" << this->commRank() << "]: " << *capi << "--(" << suppi.color() << ")-->" << *suppi << std::endl;
1249:           }
1250:         }
1251:         PetscSynchronizedPrintf(this->comm(), txt.str().c_str());
1252:         PetscSynchronizedFlush(this->comm());
1253:         //
1254:         ostringstream txt1;
1255:         if(this->commRank() == 0) {
1256:           txt1 << "base --> cap:\n";
1257:         }
1258:         if(base.empty()) {
1259:           txt1 << "[" << this->commRank() << "]: empty" << std::endl;
1260:         }
1261:         for(typename traits::baseSequence::iterator basei = base.begin(); basei != base.end(); basei++) {
1262:           const Obj<typename traits::coneSequence>& cone = this->cone(*basei);
1263:           const typename traits::coneSequence::iterator coneEnd = cone->end();

1265:           for(typename traits::coneSequence::iterator conei = cone->begin(); conei != coneEnd; conei++) {
1266:             txt1 << "[" << this->commRank() << "]: " << *basei << "<--(" << conei.color() << ")--" << *conei << std::endl;
1267:           }
1268:         }
1269:         //
1270:         PetscSynchronizedPrintf(this->comm(), txt1.str().c_str());
1271:         PetscSynchronizedFlush(this->comm());
1272:         //
1273:         ostringstream txt2;
1274:         if(this->commRank() == 0) {
1275:           txt2 << "cap <point, outdegree, marker, depth, height>:\n";
1276:         }
1277:         txt2 << "[" << this->commRank() << "]:  [";
1278:         for(typename traits::capSequence::iterator capi = cap.begin(); capi != cap.end(); capi++) {
1279:           txt2 << " <" << *capi << ", " << capi.degree() << ", " << capi.marker() << ", " << capi.depth() << ", " << capi.height() << ">";
1280:         }
1281:         txt2 << " ]" << std::endl;
1282:         //
1283:         PetscSynchronizedPrintf(this->comm(), txt2.str().c_str());
1284:         PetscSynchronizedFlush(this->comm());
1285:         //
1286:         ostringstream txt3;
1287:         if(this->commRank() == 0) {
1288:           txt3 << "base <point, indegree, marker, depth, height>:\n";
1289:         }
1290:         txt3 << "[" << this->commRank() << "]:  [";
1291:         for(typename traits::baseSequence::iterator basei = base.begin(); basei != base.end(); basei++) {
1292:           txt3 << " <" << *basei << "," << basei.degree() << ", " << basei.marker() << ", " << basei.depth() << ", " << basei.height() << ">";
1293:         }
1294:         txt3 << " ]" << std::endl;
1295:         //
1296:         PetscSynchronizedPrintf(this->comm(), txt3.str().c_str());
1297:         PetscSynchronizedFlush(this->comm());
1298:     };
1299:     //
1300:     // Structural queries
1301:     //
1302:     template <typename Point_, typename Marker_, typename Color_>
1303:     void Sieve<Point_,Marker_,Color_>::setMarker(const point_type& p, const marker_type& marker) {
1304:       typename ::boost::multi_index::index<typename traits::base_container_type::set_type,typename traits::base_container_type::traits::pointTag>::type& bIndex = ::boost::multi_index::get<typename traits::base_container_type::traits::pointTag>(this->_base.set);
1305:       typename ::boost::multi_index::index<typename traits::cap_container_type::set_type,typename traits::cap_container_type::traits::pointTag>::type& cIndex = ::boost::multi_index::get<typename traits::cap_container_type::traits::pointTag>(this->_cap.set);

1307:       if (bIndex.find(p) != bIndex.end()) {
1308:         bIndex.modify(bIndex.find(p), changeMarker(marker));
1309:       }
1310:       if (cIndex.find(p) != cIndex.end()) {
1311:         cIndex.modify(cIndex.find(p), changeMarker(marker));
1312:       }
1313:       this->_markers->insert(marker);
1314:     };

1316:     template <typename Point_, typename Marker_, typename Color_>
1317:     template <typename Sequence>
1318:     void Sieve<Point_,Marker_,Color_>::setMarker(const Obj<Sequence>& points, const marker_type& marker) {
1319:       for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1320:         this->setMarker(*p_iter, marker);
1321:       }
1322:     };

1324:     template <typename Point_, typename Marker_, typename Color_>
1325:     int Sieve<Point_,Marker_,Color_>::depth() {
1326:       return this->maxDepth;
1327:     };
1328:     template <typename Point_, typename Marker_, typename Color_>
1329:     int Sieve<Point_,Marker_,Color_>::depth(const point_type& p) {
1330:       const typename ::boost::multi_index::index<typename traits::base_container_type::set_type,typename traits::cap_container_type::traits::pointTag>::type& i = ::boost::multi_index::get<typename traits::base_container_type::traits::pointTag>(this->_base.set);
1331:       if (i.find(p) != i.end()) {
1332:         return i.find(p)->depth;
1333:       }
1334:       return 0;
1335:     };
1336:     template <typename Point_, typename Marker_, typename Color_>
1337:     template<typename InputSequence>
1338:     int Sieve<Point_,Marker_,Color_>::depth(const Obj<InputSequence>& points) {
1339:       const typename ::boost::multi_index::index<typename traits::base_container_type::set_type,typename traits::cap_container_type::traits::pointTag>::type& i = ::boost::multi_index::get<typename traits::base_container_type::traits::pointTag>(this->_base.set);
1340:       int maxDepth = 0;
1341: 
1342:       for(typename InputSequence::iterator iter = points->begin(); iter != points->end(); ++iter) {
1343:         if (i.find(*iter) != i.end()) {
1344:           maxDepth = std::max(maxDepth, i.find(*iter)->depth);
1345:         }
1346:       }
1347:       return maxDepth;
1348:     };
1349:     template <typename Point_, typename Marker_, typename Color_>
1350:     int Sieve<Point_,Marker_,Color_>::height() {
1351:       return this->maxHeight;
1352:     };
1353:     template <typename Point_, typename Marker_, typename Color_>
1354:     int Sieve<Point_,Marker_,Color_>::height(const point_type& p) {
1355:       const typename ::boost::multi_index::index<typename traits::cap_container_type::set_type,typename traits::cap_container_type::traits::pointTag>::type& i = ::boost::multi_index::get<typename traits::cap_container_type::traits::pointTag>(this->_cap.set);
1356:       if (i.find(p) != i.end()) {
1357:         return i.find(p)->height;
1358:       }
1359:       return 0;
1360:     };
1361:     template <typename Point_, typename Marker_, typename Color_>
1362:     template<typename InputSequence>
1363:     int Sieve<Point_,Marker_,Color_>::height(const Obj<InputSequence>& points) {
1364:       const typename ::boost::multi_index::index<typename traits::cap_container_type::set_type,typename traits::cap_container_type::traits::pointTag>::type& i = ::boost::multi_index::get<typename traits::cap_container_type::traits::pointTag>(this->_cap.set);
1365:       int maxHeight = 0;
1366: 
1367:       for(typename InputSequence::iterator iter = points->begin(); iter != points->end(); ++iter) {
1368:         if (i.find(*iter) != i.end()) {
1369:           maxHeight = std::max(maxHeight, i.find(*iter)->height);
1370:         }
1371:       }
1372:       return maxHeight;
1373:     };

1375:     template <typename Point_, typename Marker_, typename Color_>
1376:     int Sieve<Point_,Marker_,Color_>::diameter() {
1377:       int globalDiameter;
1378:       int MPI_Allreduce(&this->graphDiameter, &globalDiameter, 1, MPI_INT, MPI_MAX, this->comm());
1379:       CHKMPIERROR(ierr, ERRORMSG("Error in MPI_Allreduce"));
1380:       return globalDiameter;
1381:     };
1382:     template <typename Point_, typename Marker_, typename Color_>
1383:     int Sieve<Point_,Marker_,Color_>::diameter(const point_type& p) {
1384:       return this->depth(p) + this->height(p);
1385:     };

1387:     template <typename Point_, typename Marker_, typename Color_>
1388:     Obj<typename Sieve<Point_,Marker_,Color_>::traits::depthSequence> Sieve<Point_,Marker_,Color_>::depthStratum(int d) {
1389:       if (d == 0) {
1390:         return typename traits::depthSequence(::boost::multi_index::get<typename traits::cap_container_type::traits::depthMarkerTag>(this->_cap.set), d);
1391:       } else {
1392:         return typename traits::depthSequence(::boost::multi_index::get<typename traits::base_container_type::traits::depthMarkerTag>(this->_base.set), d);
1393:       }
1394:     };
1395:     template <typename Point_, typename Marker_, typename Color_>
1396:     Obj<typename Sieve<Point_,Marker_,Color_>::traits::depthSequence> Sieve<Point_,Marker_,Color_>::depthStratum(int d, marker_type m) {
1397:       if (d == 0) {
1398:         return typename traits::depthSequence(::boost::multi_index::get<typename traits::cap_container_type::traits::depthMarkerTag>(this->_cap.set), d, m);
1399:       } else {
1400:         return typename traits::depthSequence(::boost::multi_index::get<typename traits::base_container_type::traits::depthMarkerTag>(this->_base.set), d, m);
1401:       }
1402:     };
1403:     template <typename Point_, typename Marker_, typename Color_>
1404:     Obj<typename Sieve<Point_,Marker_,Color_>::traits::heightSequence> Sieve<Point_,Marker_,Color_>::heightStratum(int h) {
1405:       if (h == 0) {
1406:         return typename traits::heightSequence(::boost::multi_index::get<typename traits::base_container_type::traits::heightMarkerTag>(this->_base.set), h);
1407:       } else {
1408:         return typename traits::heightSequence(::boost::multi_index::get<typename traits::cap_container_type::traits::heightMarkerTag>(this->_cap.set), h);
1409:       }
1410:     };
1411:     template <typename Point_, typename Marker_, typename Color_>
1412:     Obj<typename Sieve<Point_,Marker_,Color_>::traits::heightSequence> Sieve<Point_,Marker_,Color_>::heightStratum(int h, marker_type m) {
1413:       if (h == 0) {
1414:         return typename traits::heightSequence(::boost::multi_index::get<typename traits::base_container_type::traits::heightMarkerTag>(this->_base.set), h, m);
1415:       } else {
1416:         return typename traits::heightSequence(::boost::multi_index::get<typename traits::cap_container_type::traits::heightMarkerTag>(this->_cap.set), h, m);
1417:       }
1418:     };
1419:     template <typename Point_, typename Marker_, typename Color_>
1420:     template<class InputSequence>
1421:     void Sieve<Point_,Marker_,Color_>::__computeClosureHeights(const Obj<InputSequence>& points) {
1422:       typename ::boost::multi_index::index<typename traits::cap_container_type::set_type,typename traits::cap_container_type::traits::pointTag>::type& index = ::boost::multi_index::get<typename traits::cap_container_type::traits::pointTag>(this->_cap.set);
1423:       typename ::boost::multi_index::index<typename traits::base_container_type::set_type,typename traits::base_container_type::traits::pointTag>::type& bIndex = ::boost::multi_index::get<typename traits::base_container_type::traits::pointTag>(this->_base.set);
1424:       Obj<coneSet> modifiedPoints = coneSet();
1425: 
1426:       for(typename InputSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
1427:         // Compute the max height of the points in the support of p, and add 1
1428:         int h0 = this->height(*p_itor);
1429:         int h1 = this->height(this->support(*p_itor)) + 1;
1430:         if(h1 != h0) {
1431:           typename ::boost::multi_index::index<typename traits::base_container_type::set_type,typename traits::base_container_type::traits::pointTag>::type::iterator bIter = bIndex.find(*p_itor);

1433:           index.modify(index.find(*p_itor), changeHeight(h1));
1434:           if (bIter != bIndex.end()) {
1435:             bIndex.modify(bIter, changeHeight(h1));
1436:           }
1437:           if (h1 > this->maxHeight) this->maxHeight = h1;
1438:           modifiedPoints->insert(*p_itor);
1439:         }
1440:       }
1441:       // FIX: We would like to avoid the copy here with cone()
1442:       if(modifiedPoints->size() > 0) {
1443:         this->__computeClosureHeights(this->cone(modifiedPoints));
1444:       }
1445:     };
1446:     template <typename Point_, typename Marker_, typename Color_>
1447:     template<class InputSequence>
1448:     void Sieve<Point_,Marker_,Color_>::__computeStarDepths(const Obj<InputSequence>& points) {
1449:       typename ::boost::multi_index::index<typename traits::base_container_type::set_type,typename traits::base_container_type::traits::pointTag>::type& index = ::boost::multi_index::get<typename traits::base_container_type::traits::pointTag>(this->_base.set);
1450:       typename ::boost::multi_index::index<typename traits::cap_container_type::set_type,typename traits::cap_container_type::traits::pointTag>::type& cIndex = ::boost::multi_index::get<typename traits::cap_container_type::traits::pointTag>(this->_cap.set);
1451:       Obj<supportSet> modifiedPoints = supportSet();
1452:       for(typename InputSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
1453:         // Compute the max depth of the points in the support of p, and add 1
1454:         int d0 = this->depth(*p_itor);
1455:         int d1 = this->depth(this->cone(*p_itor)) + 1;
1456:         if(d1 != d0) {
1457:           typename ::boost::multi_index::index<typename traits::cap_container_type::set_type,typename traits::cap_container_type::traits::pointTag>::type::iterator cIter = cIndex.find(*p_itor);

1459:           index.modify(index.find(*p_itor), changeDepth(d1));
1460:           if (cIter != cIndex.end()) {
1461:             cIndex.modify(cIter, changeDepth(d1));
1462:           }
1463:           if (d1 > this->maxDepth) this->maxDepth = d1;
1464:           modifiedPoints->insert(*p_itor);
1465:         }
1466:       }
1467:       // FIX: We would like to avoid the copy here with cone()
1468:       if(modifiedPoints->size() > 0) {
1469:         this->__computeStarDepths(this->support(modifiedPoints));
1470:       }
1471:     };
1474:     template <typename Point_, typename Marker_, typename Color_>
1475:     void Sieve<Point_,Marker_,Color_>::stratify(bool show) {
1476:       ALE_LOG_EVENT_BEGIN;
1477:       // FIX: We would like to avoid the copy here with cone() and support()
1478:       this->__computeClosureHeights(this->cone(this->leaves()));
1479:       this->__computeStarDepths(this->support(this->roots()));
1480: 
1481:       Obj<typename traits::capSequence> base = this->base();

1483:       for(typename traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1484:         maxDepth = std::max(maxDepth, b_iter.depth());
1485:         //b_iter.setDegree(this->cone(*b_iter)->size());
1486:         this->_base.adjustDegree(*b_iter, this->cone(*b_iter)->size());
1487:       }
1488:       Obj<typename traits::capSequence> cap = this->cap();

1490:       for(typename traits::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1491:         maxHeight = std::max(maxHeight, c_iter.height());
1492:         //c_iter.setDegree(this->support(*c_iter)->size());
1493:         this->_cap.adjustDegree(*c_iter, this->support(*c_iter)->size());
1494:       }
1495:       if (this->debug() || show) {
1496: //         const typename ::boost::multi_index::index<StratumSet,point>::type& points = ::boost::multi_index::get<point>(this->strata);
1497: //         for(typename ::boost::multi_index::index<StratumSet,point>::type::iterator i = points.begin(); i != points.end(); i++) {
1498: //           std::cout << *i << std::endl;
1499: //         }
1500:       }
1501:       ALE_LOG_EVENT_END;
1502:     };
1503:     //
1504:     // Structural manipulation
1505:     //
1506: 
1507: } // namespace ALE

1509: #endif