Actual source code: Sieve.hh
petsc-3.3-p7 2013-05-11
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_Sifter_hh
12: #include <sieve/Sifter.hh>
13: #endif
16: namespace ALE {
18: namespace SieveDef {
19: //
20: // Rec & RecContainer definitions.
21: // Rec is intended to denote a graph point record.
22: //
23: template <typename Point_, typename Marker_>
24: struct Rec {
25: typedef Point_ point_type;
26: typedef Marker_ marker_type;
27: template<typename OtherPoint_, typename OtherMarker_ = Marker_>
28: struct rebind {
29: typedef Rec<OtherPoint_, OtherMarker_> type;
30: };
31: point_type point;
32: int degree;
33: int depth;
34: int height;
35: marker_type marker;
36: // Basic interface
37: Rec() : point(point_type()), degree(0), depth(0), height(0), marker(marker_type()) {};
38: Rec(const Rec& r) : point(r.point), degree(r.degree), depth(r.depth), height(r.height), marker(r.marker) {}
39: Rec(const point_type& p) : point(p), degree(0), depth(0), height(0), marker(marker_type()) {};
40: Rec(const point_type& p, const int degree) : point(p), degree(degree), depth(0), height(0), marker(marker_type()) {};
41: 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) {};
42: // Printing
43: friend std::ostream& operator<<(std::ostream& os, const Rec& p) {
44: os << "<" << p.point << ", "<< p.degree << ", "<< p.depth << ", "<< p.height << ", "<< p.marker << ">";
45: return os;
46: };
48: struct degreeAdjuster {
49: degreeAdjuster(int newDegree) : _newDegree(newDegree) {};
50: void operator()(Rec& r) { r.degree = this->_newDegree; }
51: private:
52: int _newDegree;
53: };
54: };// class Rec
56: template <typename Point_, typename Rec_>
57: struct RecContainerTraits {
58: typedef Rec_ rec_type;
59: typedef typename rec_type::marker_type marker_type;
60: // Index tags
61: struct pointTag{};
62: struct degreeTag{};
63: struct markerTag{};
64: struct depthMarkerTag{};
65: struct heightMarkerTag{};
66: // Rec set definition
67: typedef ::boost::multi_index::multi_index_container<
68: rec_type,
69: ::boost::multi_index::indexed_by<
70: ::boost::multi_index::ordered_unique<
71: ::boost::multi_index::tag<pointTag>, BOOST_MULTI_INDEX_MEMBER(rec_type, typename rec_type::point_type, point)
72: >,
73: // ::boost::multi_index::ordered_non_unique<
74: // ::boost::multi_index::tag<degreeTag>, BOOST_MULTI_INDEX_MEMBER(rec_type, int, degree)
75: // >,
76: ::boost::multi_index::ordered_non_unique<
77: ::boost::multi_index::tag<markerTag>, BOOST_MULTI_INDEX_MEMBER(rec_type, marker_type, marker)
78: >,
79: ::boost::multi_index::ordered_non_unique<
80: ::boost::multi_index::tag<depthMarkerTag>,
81: ::boost::multi_index::composite_key<
82: rec_type, BOOST_MULTI_INDEX_MEMBER(rec_type,int,depth), BOOST_MULTI_INDEX_MEMBER(rec_type,marker_type,marker)>
83: >,
84: ::boost::multi_index::ordered_non_unique<
85: ::boost::multi_index::tag<heightMarkerTag>,
86: ::boost::multi_index::composite_key<
87: rec_type, BOOST_MULTI_INDEX_MEMBER(rec_type,int,height), BOOST_MULTI_INDEX_MEMBER(rec_type,marker_type,marker)>
88: >
89: >,
90: ALE_ALLOCATOR<rec_type>
91: > set_type;
92: //
93: // Return types
94: //
96: class PointSequence {
97: public:
98: typedef ALE::SifterDef::IndexSequenceTraits<typename ::boost::multi_index::index<set_type, pointTag>::type,
99: BOOST_MULTI_INDEX_MEMBER(rec_type, typename rec_type::point_type,point)>
100: traits;
101: protected:
102: const typename traits::index_type& _index;
103: public:
104:
105: // Need to extend the inherited iterator to be able to extract the degree
106: class iterator : public traits::iterator {
107: public:
108: iterator(const typename traits::iterator::itor_type& itor) : traits::iterator(itor) {};
109: virtual const int& degree() const {return this->_itor->degree;};
110: virtual const int& marker() const {return this->_itor->marker;};
111: virtual const int& depth() const {return this->_itor->depth;};
112: virtual const int& height() const {return this->_itor->height;};
113: //void setDegree(const int degree) {this->_index.modify(this->_itor, typename traits::rec_type::degreeAdjuster(degree));};
114: };
115:
116: PointSequence(const PointSequence& seq) : _index(seq._index) {};
117: PointSequence(typename traits::index_type& index) : _index(index) {};
118: virtual ~PointSequence(){};
119:
120: virtual bool empty(){return this->_index.empty();};
121:
122: virtual typename traits::index_type::size_type size() {return this->_index.size();};
124: virtual iterator begin() const {
125: // Retrieve the beginning iterator of the index
126: return iterator(this->_index.begin());
127: };
128: virtual iterator end() {
129: // Retrieve the ending iterator of the index
130: // Since the elements in this index are ordered by degree, this amounts to the end() of the index.
131: return iterator(this->_index.end());
132: };
133: virtual bool contains(const typename rec_type::point_type& p) {
134: // Check whether a given point is in the index
135: return (this->_index.find(p) != this->_index.end());
136: };
137: }; // class PointSequence
139: template<typename Tag_, typename Value_>
140: class ValueSequence {
141: public:
142: typedef Value_ value_type;
143: typedef ALE::SifterDef::IndexSequenceTraits<typename ::boost::multi_index::index<set_type, Tag_>::type,
144: BOOST_MULTI_INDEX_MEMBER(rec_type, typename rec_type::point_type,point)>
145: traits;
146: protected:
147: const typename traits::index_type& _index;
148: const value_type _value;
149: public:
150: // Need to extend the inherited iterator to be able to extract the degree
151: class iterator : public traits::iterator {
152: public:
153: iterator(const typename traits::iterator::itor_type& itor) : traits::iterator(itor) {};
154: virtual const int& degree() const {return this->_itor->degree;};
155: virtual const int& marker() const {return this->_itor->marker;};
156: virtual const int& depth() const {return this->_itor->depth;};
157: virtual const int& height() const {return this->_itor->height;};
158: };
159:
160: ValueSequence(const ValueSequence& seq) : _index(seq._index), _value(seq._value) {};
161: ValueSequence(typename traits::index_type& index, const value_type& value) : _index(index), _value(value) {};
162: virtual ~ValueSequence(){};
163:
164: virtual bool empty(){return this->_index.empty();};
165:
166: virtual typename traits::index_type::size_type size() {return this->_index.count(this->_value);};
168: virtual iterator begin() {
169: return iterator(this->_index.lower_bound(this->_value));
170: };
171: virtual iterator end() {
172: return iterator(this->_index.upper_bound(this->_value));
173: };
174: }; // class ValueSequence
176: template<typename Tag_, typename Value_>
177: class TwoValueSequence {
178: public:
179: typedef Value_ value_type;
180: typedef ALE::SifterDef::IndexSequenceTraits<typename ::boost::multi_index::index<set_type, Tag_>::type,
181: BOOST_MULTI_INDEX_MEMBER(rec_type, typename rec_type::point_type,point)>
182: traits;
183: protected:
184: const typename traits::index_type& _index;
185: const value_type _valueA, _valueB;
186: const bool _useTwoValues;
187: public:
188: // Need to extend the inherited iterator to be able to extract the degree
189: class iterator : public traits::iterator {
190: public:
191: iterator(const typename traits::iterator::itor_type& itor) : traits::iterator(itor) {};
192: virtual const int& degree() const {return this->_itor->degree;};
193: virtual const int& marker() const {return this->_itor->marker;};
194: };
195:
196: TwoValueSequence(const TwoValueSequence& seq) : _index(seq._index), _valueA(seq._valueA), _valueB(seq._valueB), _useTwoValues(seq._useTwoValues) {};
197: TwoValueSequence(typename traits::index_type& index, const value_type& valueA) : _index(index), _valueA(valueA), _valueB(value_type()), _useTwoValues(false) {};
198: TwoValueSequence(typename traits::index_type& index, const value_type& valueA, const value_type& valueB) : _index(index), _valueA(valueA), _valueB(valueB), _useTwoValues(true) {};
199: virtual ~TwoValueSequence(){};
200:
201: virtual bool empty(){return this->_index.empty();};
202:
203: virtual typename traits::index_type::size_type size() {
204: if (this->_useTwoValues) {
205: return this->_index.count(::boost::make_tuple(this->_valueA,this->_valueB));
206: } else {
207: return this->_index.count(::boost::make_tuple(this->_valueA));
208: }
209: };
211: virtual iterator begin() {
212: if (this->_useTwoValues) {
213: return iterator(this->_index.lower_bound(::boost::make_tuple(this->_valueA,this->_valueB)));
214: } else {
215: return iterator(this->_index.lower_bound(::boost::make_tuple(this->_valueA)));
216: }
217: };
218: virtual iterator end() {
219: if (this->_useTwoValues) {
220: return iterator(this->_index.upper_bound(::boost::make_tuple(this->_valueA,this->_valueB)));
221: } else {
222: return iterator(this->_index.upper_bound(::boost::make_tuple(this->_valueA)));
223: }
224: };
225: }; // class TwoValueSequence
226: };// struct RecContainerTraits
228: template <typename Point_, typename Rec_>
229: struct RecContainer {
230: typedef RecContainerTraits<Point_, Rec_> traits;
231: typedef typename traits::set_type set_type;
232: template <typename OtherPoint_, typename OtherRec_>
233: struct rebind {
234: typedef RecContainer<OtherPoint_, OtherRec_> type;
235: };
236: set_type set;
238: void removePoint(const typename traits::rec_type::point_type& p) {
239: /*typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type& index =
240: ::boost::multi_index::get<typename traits::pointTag>(this->set);
241: typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type::iterator i = index.find(p);
242: if (i != index.end()) { // Point exists
243: index.erase(i);
244: }*/
245: this->set.erase(p);
246: };
247: void adjustDegree(const typename traits::rec_type::point_type& p, int delta) {
248: typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type& index =
249: ::boost::multi_index::get<typename traits::pointTag>(this->set);
250: typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type::iterator i = index.find(p);
251: if (i == index.end()) { // No such point exists
252: if(delta < 0) { // Cannot decrease degree of a non-existent point
253: ostringstream err;
254: err << "ERROR: adjustDegree: Non-existent point " << p;
255: std::cout << err << std::endl;
256: throw(Exception(err.str().c_str()));
257: }
258: else { // We CAN INCREASE the degree of a non-existent point: simply insert a new element with degree == delta
259: std::pair<typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type::iterator, bool> ii;
260: typename traits::rec_type r(p,delta);
261: ii = index.insert(r);
262: if(ii.second == false) {
263: ostringstream err;
264: err << "ERROR: adjustDegree: Failed to insert a rec " << r;
265: std::cout << err << std::endl;
266: throw(Exception(err.str().c_str()));
267: }
268: }
269: }
270: else { // Point exists, so we try to modify its degree
271: // If the adjustment is zero, there is nothing to do, otherwise ...
272: if(delta != 0) {
273: int newDegree = i->degree + delta;
274: if(newDegree < 0) {
275: ostringstream ss;
276: ss << "adjustDegree: Adjustment of " << *i << " by " << delta << " would result in negative degree: " << newDegree;
277: throw Exception(ss.str().c_str());
278: }
279: index.modify(i, typename traits::rec_type::degreeAdjuster(newDegree));
280: }
281: }
282: }; // adjustDegree()
283: }; // struct RecContainer
284: };// namespace SieveDef
286: //
287: // Sieve:
288: // A Sieve is a set of {\emph arrows} connecting {\emph points} of type Point_. Thus we
289: // could realize a sieve, for instance, as a directed graph. In addition, we will often
290: // assume an acyclicity constraint, arriving at a DAG. Each arrow may also carry a label,
291: // or {\emph color} of type Color_, and the interface allows sets of arrows to be filtered
292: // by color.
293: //
294: template <typename Point_, typename Marker_, typename Color_>
295: 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_> > > {
296: public:
297: typedef Color_ color_type;
298: typedef Point_ point_type;
299: typedef Marker_ marker_type;
300: typedef struct {
301: 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;
302: // Encapsulated container types
303: typedef typename baseType::traits::arrow_container_type arrow_container_type;
304: typedef typename baseType::traits::cap_container_type cap_container_type;
305: typedef typename baseType::traits::base_container_type base_container_type;
306: // Types associated with records held in containers
307: typedef typename baseType::traits::arrow_type arrow_type;
308: typedef typename baseType::traits::source_type source_type;
309: typedef typename baseType::traits::sourceRec_type sourceRec_type;
310: typedef typename baseType::traits::target_type target_type;
311: typedef typename baseType::traits::targetRec_type targetRec_type;
312: typedef typename baseType::traits::color_type color_type;
313: typedef Point_ point_type;
314: // Convenient tag names
315: typedef typename baseType::traits::supportInd supportInd;
316: typedef typename baseType::traits::coneInd coneInd;
317: typedef typename baseType::traits::arrowInd arrowInd;
318: typedef typename baseType::traits::baseInd baseInd;
319: typedef typename baseType::traits::capInd capInd;
321: //
322: // Return types
323: //
324: typedef typename baseType::traits::arrowSequence arrowSequence;
325: typedef typename baseType::traits::coneSequence coneSequence;
326: typedef typename baseType::traits::supportSequence supportSequence;
327: typedef typename baseType::traits::baseSequence baseSequence;
328: typedef typename baseType::traits::capSequence capSequence;
329: typedef typename base_container_type::traits::template TwoValueSequence<typename base_container_type::traits::depthMarkerTag,int> depthSequence;
330: typedef typename cap_container_type::traits::template TwoValueSequence<typename cap_container_type::traits::heightMarkerTag,int> heightSequence;
331: typedef typename cap_container_type::traits::template ValueSequence<typename cap_container_type::traits::markerTag,marker_type> markerSequence;
332: } traits;
333: typedef std::set<point_type> pointSet;
334: typedef ALE::array<point_type> pointArray;
335: typedef std::set<marker_type> markerSet;
336: typedef pointSet coneSet;
337: typedef pointSet supportSet;
338: typedef pointArray coneArray;
339: typedef pointArray supportArray;
340: public:
341: 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) {
342: this->_markers = new markerSet();
343: this->_meetSet = new coneSet();
344: //std::cout << "["<<this->commRank()<<"]: Creating an ALE::Sieve" << std::endl;
345: };
346: virtual ~Sieve() {
347: //std::cout << "["<<this->commRank()<<"]: Destroying an ALE::Sieve" << std::endl;
348: };
349: // Printing
350: friend std::ostream& operator<<(std::ostream& os, Obj<Sieve<Point_,Marker_,Color_> > s) {
351: os << *s;
352: return os;
353: };
354:
355: friend std::ostream& operator<<(std::ostream& os, Sieve<Point_,Marker_,Color_>& s) {
356: Obj<typename traits::baseSequence> base = s.base();
357: for(typename traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
358: Obj<typename traits::coneSequence> cone = s.cone(*b_iter);
359: os << "Base point " << *b_iter << " with cone:" << std::endl;
360: for(typename traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
361: os << " " << *c_iter << std::endl;
362: }
363: }
364: return os;
365: };
367: template<typename ostream_type>
368: void view(ostream_type& os, const char* label = NULL, bool rawData = false);
369: void view(const char* label = NULL, MPI_Comm comm = MPI_COMM_NULL);
371: Obj<Sieve> copy() {
372: Obj<Sieve> s = Sieve(this->comm(), this->debug);
373: Obj<typename traits::capSequence> cap = this->cap();
374: Obj<typename traits::baseSequence> base = this->base();
376: for(typename traits::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
377: s->addCapPoint(*c_iter);
378: }
379: for(typename traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
380: Obj<typename traits::coneSequence> cone = this->cone(*b_iter);
382: for(typename traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
383: s->addArrow(*c_iter, *b_iter, c_iter.color());
384: }
385: }
386: s->stratify();
387: return s;
388: };
389: bool hasPoint(const point_type& point) {
390: if (this->baseContains(point) || this->capContains(point)) return true;
391: return false;
392: };
393: private:
394: template<class InputSequence> Obj<coneSet> __nCone(Obj<InputSequence>& cone, int n, const Color_& color, bool useColor);
395: template<class pointSequence> void __nCone(const Obj<pointSequence>& points, int n, const Color_& color, bool useColor, Obj<coneArray> cone, Obj<coneSet> seen);
396: template<class pointSequence> void __nSupport(const Obj<pointSequence>& points, int n, const Color_& color, bool useColor, Obj<supportArray> cone, Obj<supportSet> seen);
397: public:
398: //
399: // The basic Sieve interface (extensions to Sifter)
400: //
401: Obj<coneArray> nCone(const Point_& p, int n);
402: Obj<coneArray> nCone(const Point_& p, int n, const Color_& color, bool useColor = true);
403: template<class InputSequence> Obj<coneSet> nCone(const Obj<InputSequence>& points, int n);
404: template<class InputSequence> Obj<coneSet> nCone(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor = true);
406: Obj<supportArray> nSupport(const Point_& p, int n);
407: Obj<supportArray> nSupport(const Point_& p, int n, const Color_& color, bool useColor = true);
408: template<class InputSequence> Obj<supportSet> nSupport(const Obj<InputSequence>& points, int n);
409: template<class InputSequence> Obj<supportSet> nSupport(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor = true);
410: public:
411: virtual bool checkArrow(const typename traits::arrow_type& a) {
412: if ((this->_cap.set.find(a.source) == this->_cap.set.end()) && (this->_base.set.find(a.source) == this->_base.set.end())) return false;
413: if ((this->_cap.set.find(a.target) == this->_cap.set.end()) && (this->_base.set.find(a.target) == this->_base.set.end())) return false;
414: return true;
415: };
416: //
417: // Iterated versions
418: //
419: Obj<supportSet> star(const Point_& p);
421: Obj<supportSet> star(const Point_& p, const Color_& color);
423: template<class InputSequence>
424: Obj<supportSet> star(const Obj<InputSequence>& points);
426: template<class InputSequence>
427: Obj<supportSet> star(const Obj<InputSequence>& points, const Color_& color);
429: Obj<supportSet> nStar(const Point_& p, int n);
431: Obj<supportSet> nStar(const Point_& p, int n, const Color_& color, bool useColor = true);
433: template<class InputSequence>
434: Obj<supportSet> nStar(const Obj<InputSequence>& points, int n);
436: template<class InputSequence>
437: Obj<supportSet> nStar(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor = true);
439: private:
440: template<class InputSequence>
441: Obj<supportSet> __nStar(Obj<InputSequence>& support, int n, const Color_& color, bool useColor);
443: public:
444: //
445: // Lattice methods
446: //
447: const Obj<coneSet>& meet(const Point_& p, const Point_& q);
449: const Obj<coneSet>& meet(const Point_& p, const Point_& q, const Color_& color);
451: template<class InputSequence>
452: const Obj<coneSet>& meet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1);
454: template<class InputSequence>
455: const Obj<coneSet>& meet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, const Color_& color);
457: const Obj<coneSet>& nMeet(const Point_& p, const Point_& q, int n);
459: const Obj<coneSet>& nMeet(const Point_& p, const Point_& q, int n, const Color_& color, bool useColor = true);
461: template<class InputSequence>
462: const Obj<coneSet>& nMeet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n);
464: template<class InputSequence>
465: const Obj<coneSet>& nMeet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n,
466: const Color_& color, bool useColor = true);
468: Obj<supportSet> join(const Point_& p, const Point_& q);
470: Obj<supportSet> join(const Point_& p, const Point_& q, const Color_& color);
472: template<class InputSequence>
473: Obj<supportSet> join(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1);
475: template<class InputSequence>
476: Obj<supportSet> join(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, const Color_& color);
478: template<class InputSequence>
479: Obj<supportSet> nJoin1(const Obj<InputSequence>& chain);
481: Obj<supportSet> nJoin(const Point_& p, const Point_& q, int n);
483: Obj<supportSet> nJoin(const Point_& p, const Point_& q, int n, const Color_& color, bool useColor = true);
485: template<class InputSequence>
486: Obj<supportSet> nJoin(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n);
488: template<class InputSequence>
489: Obj<supportSet> nJoin(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n, const Color_& color, bool useColor = true);
491: template<class InputSequence>
492: Obj<supportSet> nJoin(const Obj<InputSequence>& chain, const int depth);
494: public:
495: Obj<typename traits::depthSequence> roots() {
496: return this->depthStratum(0);
497: };
498: Obj<typename traits::heightSequence> leaves() {
499: return this->heightStratum(0);
500: };
501: private:
502: bool doStratify;
503: int maxDepth, maxHeight, graphDiameter;
504: public:
505: //
506: // Structural queries
507: //
508: int depth();
509: int depth(const point_type& p);
510: template<typename InputSequence> int depth(const Obj<InputSequence>& points);
512: int height();
513: int height(const point_type& p);
514: template<typename InputSequence> int height(const Obj<InputSequence>& points);
516: int diameter();
517: int diameter(const point_type& p);
519: Obj<typename traits::depthSequence> depthStratum(int d) {
520: if (d == 0) {
521: return typename traits::depthSequence(::boost::multi_index::get<typename traits::cap_container_type::traits::depthMarkerTag>(this->_cap.set), d);
522: } else {
523: return typename traits::depthSequence(::boost::multi_index::get<typename traits::base_container_type::traits::depthMarkerTag>(this->_base.set), d);
524: }
525: };
526: Obj<typename traits::depthSequence> depthStratum(int d, marker_type m) {
527: if (d == 0) {
528: return typename traits::depthSequence(::boost::multi_index::get<typename traits::cap_container_type::traits::depthMarkerTag>(this->_cap.set), d, m);
529: } else {
530: return typename traits::depthSequence(::boost::multi_index::get<typename traits::base_container_type::traits::depthMarkerTag>(this->_base.set), d, m);
531: }
532: };
534: Obj<typename traits::heightSequence> heightStratum(int h) {
535: if (h == 0) {
536: return typename traits::heightSequence(::boost::multi_index::get<typename traits::base_container_type::traits::heightMarkerTag>(this->_base.set), h);
537: } else {
538: return typename traits::heightSequence(::boost::multi_index::get<typename traits::cap_container_type::traits::heightMarkerTag>(this->_cap.set), h);
539: }
540: };
541: Obj<typename traits::heightSequence> heightStratum(int h, marker_type m) {
542: if (h == 0) {
543: return typename traits::heightSequence(::boost::multi_index::get<typename traits::base_container_type::traits::heightMarkerTag>(this->_base.set), h, m);
544: } else {
545: return typename traits::heightSequence(::boost::multi_index::get<typename traits::cap_container_type::traits::heightMarkerTag>(this->_cap.set), h, m);
546: }
547: };
549: Obj<typename traits::markerSequence> markerStratum(marker_type m);
550:
551: void setStratification(bool doStratify) {this->doStratify = doStratify;};
553: bool getStratification() {return this->doStratify;};
555: void stratify(bool show = false);
556: protected:
557: Obj<markerSet> _markers;
558: Obj<coneSet> _meetSet;
559: public:
560: //
561: // Structural manipulation
562: //
564: struct changeMarker {
565: changeMarker(int newMarker) : newMarker(newMarker) {};
567: void operator()(typename traits::base_container_type::traits::rec_type& p) {
568: p.marker = newMarker;
569: }
570: private:
571: marker_type newMarker;
572: };
574: void setMarker(const point_type& p, const marker_type& marker);
575: template<class InputSequence> void setMarker(const Obj<InputSequence>& points, const marker_type& marker);
577: void clearMarkers() {this->_markers.clear();};
578: Obj<markerSet> markers() {return this->_markers;};
579: private:
580: struct changeHeight {
581: changeHeight(int newHeight) : newHeight(newHeight) {};
583: void operator()(typename traits::base_container_type::traits::rec_type& p) {
584: p.height = newHeight;
585: }
586: private:
587: int newHeight;
588: };
589:
590: template<class InputSequence> void __computeClosureHeights(const Obj<InputSequence>& points);
592: struct changeDepth {
593: changeDepth(int newDepth) : newDepth(newDepth) {};
595: void operator()(typename traits::base_container_type::traits::rec_type& p) {
596: p.depth = newDepth;
597: }
598: private:
599: int newDepth;
600: };
602: template<class InputSequence> void __computeStarDepths(const Obj<InputSequence>& points);
603: };
605: template <typename Point_, typename Marker_, typename Color_>
606: Obj<typename Sieve<Point_,Marker_,Color_>::coneArray> Sieve<Point_,Marker_,Color_>::nCone(const Point_& p, int n) {
607: return this->nCone(p, n, Color_(), false);
608: };
610: template <typename Point_, typename Marker_, typename Color_>
611: Obj<typename Sieve<Point_,Marker_,Color_>::coneArray> Sieve<Point_,Marker_,Color_>::nCone(const Point_& p, int n, const Color_& color, bool useColor) {
612: Obj<coneArray> cone = new coneArray();
613: Obj<coneSet> seen = new coneSet();
615: this->__nCone(this->cone(p), n-1, color, useColor, cone, seen);
616: return cone;
617: };
619: template <typename Point_, typename Marker_, typename Color_>
620: template<class pointSequence>
621: void Sieve<Point_,Marker_,Color_>::__nCone(const Obj<pointSequence>& points, int n, const Color_& color, bool useColor, Obj<coneArray> cone, Obj<coneSet> seen) {
622: if (n == 0) {
623: for(typename pointSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
624: if (seen->find(*p_itor) == seen->end()) {
625: cone->push_back(*p_itor);
626: seen->insert(*p_itor);
627: }
628: }
629: } else {
630: typename pointSequence::iterator end = points->end();
631: for(typename pointSequence::iterator p_itor = points->begin(); p_itor != end; ++p_itor) {
632: if (useColor) {
633: this->__nCone(this->cone(*p_itor, color), n-1, color, useColor, cone, seen);
634: } else {
635: this->__nCone(this->cone(*p_itor), n-1, color, useColor, cone, seen);
636: }
637: }
638: }
639: };
641: template <typename Point_, typename Marker_, typename Color_>
642: template<class InputSequence>
643: Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::nCone(const Obj<InputSequence>& points, int n) {
644: return this->nCone(points, n, Color_(), false);
645: };
647: template <typename Point_, typename Marker_, typename Color_>
648: template<class InputSequence>
649: Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::nCone(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor ) {
650: Obj<coneSet> cone = new coneSet();
651: cone->insert(points->begin(), points->end());
652: return this->__nCone(cone, n, color, useColor);
653: };
655: template <typename Point_, typename Marker_, typename Color_>
656: template<class InputSequence>
657: Obj<typename Sieve<Point_,Marker_,Color_>::coneSet> Sieve<Point_,Marker_,Color_>::__nCone(Obj<InputSequence>& cone, int n, const Color_& color, bool useColor) {
658: Obj<coneSet> base = new coneSet();
660: for(int i = 0; i < n; ++i) {
661: Obj<coneSet> tmp = cone; cone = base; base = tmp;
662:
663: cone->clear();
664: for(typename coneSet::iterator b_itor = base->begin(); b_itor != base->end(); ++b_itor) {
665: Obj<typename traits::coneSequence> pCone;
666:
667: if (useColor) {
668: pCone = this->cone(*b_itor, color);
669: } else {
670: pCone = this->cone(*b_itor);
671: }
672: cone->insert(pCone->begin(), pCone->end());
673: }
674: }
675: return cone;
676: };
678: template <typename Point_, typename Marker_, typename Color_>
679: Obj<typename Sieve<Point_,Marker_,Color_>::supportArray> Sieve<Point_,Marker_,Color_>::nSupport(const Point_& p, int n) {
680: return this->nSupport(p, n, Color_(), false);
681: };
683: template <typename Point_, typename Marker_, typename Color_>
684: Obj<typename Sieve<Point_,Marker_,Color_>::supportArray> Sieve<Point_,Marker_,Color_>::nSupport(const Point_& p, int n, const Color_& color, bool useColor) {
685: Obj<supportArray> cone = new supportArray();
686: Obj<supportSet> seen = new supportSet();
688: this->__nSupport(this->support(p), n-1, color, useColor, cone, seen);
689: return cone;
690: };
692: template <typename Point_, typename Marker_, typename Color_>
693: template<class pointSequence>
694: void Sieve<Point_,Marker_,Color_>::__nSupport(const Obj<pointSequence>& points, int n, const Color_& color, bool useColor, Obj<supportArray> support, Obj<supportSet> seen) {
695: if (n == 0) {
696: for(typename pointSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
697: if (seen->find(*p_itor) == seen->end()) {
698: support->push_back(*p_itor);
699: seen->insert(*p_itor);
700: }
701: }
702: } else {
703: typename pointSequence::iterator end = points->end();
704: for(typename pointSequence::iterator p_itor = points->begin(); p_itor != end; ++p_itor) {
705: if (useColor) {
706: this->__nSupport(this->support(*p_itor, color), n-1, color, useColor, support, seen);
707: } else {
708: this->__nSupport(this->support(*p_itor), n-1, color, useColor, support, seen);
709: }
710: }
711: }
712: };
714: template <typename Point_, typename Marker_, typename Color_>
715: template<class InputSequence>
716: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nSupport(const Obj<InputSequence>& points, int n) {
717: return this->nSupport(points, n, Color_(), false);
718: };
720: template <typename Point_, typename Marker_, typename Color_>
721: template<class InputSequence>
722: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nSupport(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor ) {
723: Obj<supportSet> support = supportSet();
724: Obj<supportSet> cap = supportSet();
725:
726: support->insert(points->begin(), points->end());
727: for(int i = 0; i < n; ++i) {
728: Obj<supportSet> tmp = support; support = cap; cap = tmp;
729:
730: support->clear();
731: for(typename supportSet::iterator c_itor = cap->begin(); c_itor != cap->end(); ++c_itor) {
732: Obj<typename traits::supportSequence> pSupport;
733:
734: if (useColor) {
735: pSupport = this->support(*c_itor, color);
736: } else {
737: pSupport = this->support(*c_itor);
738: }
739: support->insert(pSupport->begin(), pSupport->end());
740: }
741: }
742: return support;
743: };
744: //
745: // Iterated versions
746: //
747: template <typename Point_, typename Marker_, typename Color_>
748: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::star(const Point_& p) {
749: return nStar(p, this->height());
750: };
751:
752: template <typename Point_, typename Marker_, typename Color_>
753: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::star(const Point_& p, const Color_& color) {
754: return nStar(p, this->depth(), color);
755: };
756:
757: template <typename Point_, typename Marker_, typename Color_>
758: template<class InputSequence>
759: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::star(const Obj<InputSequence>& points) {
760: return nStar(points, this->height());
761: };
762:
763: template <typename Point_, typename Marker_, typename Color_>
764: template<class InputSequence>
765: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::star(const Obj<InputSequence>& points, const Color_& color) {
766: return nStar(points, this->height(), color);
767: };
769: template <typename Point_, typename Marker_, typename Color_>
770: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nStar(const Point_& p, int n) {
771: return this->nStar(p, n, Color_(), false);
772: };
773:
774: template <typename Point_, typename Marker_, typename Color_>
775: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nStar(const Point_& p, int n, const Color_& color, bool useColor ) {
776: Obj<supportSet> support = supportSet();
777: support->insert(p);
778: return this->__nStar(support, n, color, useColor);
779: };
781: template <typename Point_, typename Marker_, typename Color_>
782: template<class InputSequence>
783: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nStar(const Obj<InputSequence>& points, int n) {
784: return this->nStar(points, n, Color_(), false);
785: };
787: template <typename Point_, typename Marker_, typename Color_>
788: template<class InputSequence>
789: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nStar(const Obj<InputSequence>& points, int n, const Color_& color, bool useColor ) {
790: Obj<supportSet> support = supportSet();
791: support->insert(points->begin(), points->end());
792: return this->__nStar(support, n, color, useColor);
793: };
794:
795: template <typename Point_, typename Marker_, typename Color_>
796: template<class InputSequence>
797: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::__nStar(Obj<InputSequence>& support, int n, const Color_& color, bool useColor) {
798: Obj<supportSet> cap = supportSet();
799: Obj<supportSet> star = supportSet();
800: star->insert(support->begin(), support->end());
801: for(int i = 0; i < n; ++i) {
802: Obj<supportSet> tmp = support; support = cap; cap = tmp;
803: support->clear();
804: for(typename supportSet::iterator c_itor = cap->begin(); c_itor != cap->end(); ++c_itor) {
805: Obj<typename traits::supportSequence> pSupport;
806: if (useColor) {
807: pSupport = this->support(*c_itor, color);
808: } else {
809: pSupport = this->support(*c_itor);
810: }
811: support->insert(pSupport->begin(), pSupport->end());
812: star->insert(pSupport->begin(), pSupport->end());
813: }
814: }
815: return star;
816: };
818: //
819: // Lattice methods
820: //
822: template <typename Point_, typename Marker_, typename Color_>
823: const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::meet(const Point_& p, const Point_& q) {
824: return nMeet(p, q, this->depth());
825: };
826:
827: template <typename Point_, typename Marker_, typename Color_>
828: const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::meet(const Point_& p, const Point_& q, const Color_& color) {
829: return nMeet(p, q, this->depth(), color);
830: };
832: template <typename Point_, typename Marker_, typename Color_>
833: template<class InputSequence>
834: const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::meet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1) {
835: return nMeet(chain0, chain1, this->depth());
836: };
838: template <typename Point_, typename Marker_, typename Color_>
839: template<class InputSequence>
840: const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::meet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, const Color_& color) {
841: return nMeet(chain0, chain1, this->depth(), color);
842: };
844: template <typename Point_, typename Marker_, typename Color_>
845: const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::nMeet(const Point_& p, const Point_& q, int n) {
846: if (n == 1) {
847: std::vector<point_type> vecA, vecB;
848: const Obj<typename traits::coneSequence>& coneA = this->cone(p);
849: const typename traits::coneSequence::iterator beginA = coneA->begin();
850: const typename traits::coneSequence::iterator endA = coneA->end();
851: const Obj<typename traits::coneSequence>& coneB = this->cone(q);
852: const typename traits::coneSequence::iterator beginB = coneB->begin();
853: const typename traits::coneSequence::iterator endB = coneB->end();
855: vecA.insert(vecA.begin(), beginA, endA);
856: std::sort(vecA.begin(), vecA.end());
857: vecB.insert(vecB.begin(), beginB, endB);
858: std::sort(vecB.begin(), vecB.end());
859: this->_meetSet->clear();
860: 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()));
861: return this->_meetSet;
862: }
863: return nMeet(p, q, n, Color_(), false);
864: };
866: template <typename Point_, typename Marker_, typename Color_>
867: 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 ) {
868: Obj<coneSet> chain0 = new coneSet();
869: Obj<coneSet> chain1 = new coneSet();
870: chain0->insert(p);
871: chain1->insert(q);
872: return this->nMeet(chain0, chain1, n, color, useColor);
873: };
875: template <typename Point_, typename Marker_, typename Color_>
876: template<class InputSequence>
877: const Obj<typename Sieve<Point_,Marker_,Color_>::coneSet>& Sieve<Point_,Marker_,Color_>::nMeet(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n) {
878: return this->nMeet(chain0, chain1, n, Color_(), false);
879: };
880:
881: template <typename Point_, typename Marker_, typename Color_>
882: template<class InputSequence>
883: 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){
884: // The strategy is to compute the intersection of cones over the chains, remove the intersection
885: // and use the remaining two parts -- two disjoined components of the symmetric difference of cones -- as the new chains.
886: // The intersections at each stage are accumulated and their union is the meet.
887: // 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.
888: Obj<coneSet> cone;
890: this->_meetSet->clear();
891: if((chain0->size() != 0) && (chain1->size() != 0)) {
892: for(int i = 0; i <= n; ++i) {
893: // Compute the intersection of chains and put it in meet at the same time removing it from c and cc
894: std::set<point_type> intersect;
895: //std::set_intersection(chain0->begin(), chain0->end(), chain1->begin(), chain1->end(), std::insert_iterator<coneSet>(meet, meet->begin()));
896: std::set_intersection(chain0->begin(), chain0->end(), chain1->begin(), chain1->end(), std::insert_iterator<std::set<point_type> >(intersect, intersect.begin()));
897: this->_meetSet->insert(intersect.begin(), intersect.end());
898: for(typename std::set<point_type>::iterator i_iter = intersect.begin(); i_iter != intersect.end(); ++i_iter) {
899: chain0->erase(chain0->find(*i_iter));
900: chain1->erase(chain1->find(*i_iter));
901: }
902: // 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.
903: cone = this->cone(chain0);
904: chain0->insert(cone->begin(), cone->end());
905: if(chain0->size() == 0) {
906: break;
907: }
908: cone = this->cone(chain1);
909: chain1->insert(cone->begin(), cone->end());
910: if(chain1->size() == 0) {
911: break;
912: }
913: // If both cones are empty, we should quit
914: }
915: }
916: return this->_meetSet;
917: };
918:
919: template <typename Point_, typename Marker_, typename Color_>
920: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::join(const Point_& p, const Point_& q) {
921: return this->nJoin(p, q, this->depth());
922: };
923:
924: template <typename Point_, typename Marker_, typename Color_>
925: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::join(const Point_& p, const Point_& q, const Color_& color) {
926: return this->nJoin(p, q, this->depth(), color);
927: };
929: template <typename Point_, typename Marker_, typename Color_>
930: template<class InputSequence>
931: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::join(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1) {
932: return this->nJoin(chain0, chain1, this->depth());
933: };
934:
935: template <typename Point_, typename Marker_, typename Color_>
936: template<class InputSequence>
937: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::join(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, const Color_& color) {
938: return this->nJoin(chain0, chain1, this->depth(), color);
939: };
941: // Warning: I think this can be much more efficient by eliminating copies
942: template <typename Point_, typename Marker_, typename Color_>
943: template<class InputSequence>
944: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nJoin1(const Obj<InputSequence>& chain) {
945: Obj<supportSet> join = new supportSet();
946: std::set<point_type> intersectA;
947: std::set<point_type> intersectB;
948: int p = 0;
950: //std::cout << "Doing nJoin1:" << std::endl;
951: for(typename InputSequence::iterator p_iter = chain->begin(); p_iter != chain->end(); ++p_iter) {
952: //std::cout << " point " << *p_iter << std::endl;
953: const Obj<typename traits::supportSequence>& support = this->support(*p_iter);
955: join->insert(support->begin(), support->end());
956: if (p == 0) {
957: intersectB.insert(support->begin(), support->end());
958: p++;
959: } else {
960: std::set_intersection(intersectA.begin(), intersectA.end(), join->begin(), join->end(), std::insert_iterator<std::set<point_type> >(intersectB, intersectB.begin()));
961: }
962: intersectA.clear();
963: intersectA.insert(intersectB.begin(), intersectB.end());
964: intersectB.clear();
965: join->clear();
966: //std::cout << " intersection:" << std::endl;
967: //for(typename std::set<point_type>::iterator i_iter = intersectA.begin(); i_iter != intersectA.end(); ++i_iter) {
968: // std::cout << " " << *i_iter << std::endl;
969: //}
970: }
971: join->insert(intersectA.begin(), intersectA.end());
972: return join;
973: };
975: // Warning: I think this can be much more efficient by eliminating copies
976: template <typename Point_, typename Marker_, typename Color_>
977: template<class InputSequence>
978: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nJoin(const Obj<InputSequence>& chain, const int depth) {
979: Obj<supportSet> join = new supportSet();
980: std::set<point_type> intersectA;
981: std::set<point_type> intersectB;
982: int p = 0;
984: //std::cout << "Doing nJoin1:" << std::endl;
985: for(typename InputSequence::iterator p_iter = chain->begin(); p_iter != chain->end(); ++p_iter) {
986: //std::cout << " point " << *p_iter << std::endl;
987: const Obj<supportArray> support = this->nSupport(*p_iter, depth);
989: join->insert(support->begin(), support->end());
990: if (p == 0) {
991: intersectB.insert(support->begin(), support->end());
992: p++;
993: } else {
994: std::set_intersection(intersectA.begin(), intersectA.end(), join->begin(), join->end(), std::insert_iterator<std::set<point_type> >(intersectB, intersectB.begin()));
995: }
996: intersectA.clear();
997: intersectA.insert(intersectB.begin(), intersectB.end());
998: intersectB.clear();
999: join->clear();
1000: //std::cout << " intersection:" << std::endl;
1001: //for(typename std::set<point_type>::iterator i_iter = intersectA.begin(); i_iter != intersectA.end(); ++i_iter) {
1002: // std::cout << " " << *i_iter << std::endl;
1003: //}
1004: }
1005: join->insert(intersectA.begin(), intersectA.end());
1006: return join;
1007: };
1008:
1009: template <typename Point_, typename Marker_, typename Color_>
1010: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nJoin(const Point_& p, const Point_& q, int n) {
1011: return this->nJoin(p, q, n, Color_(), false);
1012: };
1013:
1014: template <typename Point_, typename Marker_, typename Color_>
1015: 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) {
1016: Obj<supportSet> chain0 = supportSet();
1017: Obj<supportSet> chain1 = supportSet();
1018: chain0->insert(p);
1019: chain1->insert(q);
1020: return this->nJoin(chain0, chain1, n, color, useColor);
1021: };
1023: template <typename Point_, typename Marker_, typename Color_>
1024: template<class InputSequence>
1025: Obj<typename Sieve<Point_,Marker_,Color_>::supportSet> Sieve<Point_,Marker_,Color_>::nJoin(const Obj<InputSequence>& chain0, const Obj<InputSequence>& chain1, int n) {
1026: return this->nJoin(chain0, chain1, n, Color_(), false);
1027: };
1029: template <typename Point_, typename Marker_, typename Color_>
1030: template<class InputSequence>
1031: 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){
1032: // The strategy is to compute the intersection of supports over the chains, remove the intersection
1033: // and use the remaining two parts -- two disjoined components of the symmetric difference of supports -- as the new chains.
1034: // The intersections at each stage are accumulated and their union is the join.
1035: // 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.
1036: Obj<supportSet> join = supportSet();
1037: Obj<supportSet> support;
1038: // std::cout << "Computing nJoin" << std::endl;
1039: // std::cout << " chain 0:" << std::endl;
1040: // for(typename InputSequence::iterator i_iter = chain0->begin(); i_iter != chain0->end(); ++i_iter) {
1041: // std::cout << " " << *i_iter << std::endl;
1042: // }
1043: // std::cout << " chain 1:" << std::endl;
1044: // for(typename InputSequence::iterator i_iter = chain1->begin(); i_iter != chain1->end(); ++i_iter) {
1045: // std::cout << " " << *i_iter << std::endl;
1046: // }
1047: if((chain0->size() != 0) && (chain1->size() != 0)) {
1048: for(int i = 0; i <= n; ++i) {
1049: // std::cout << "Level " << i << std::endl;
1050: // Compute the intersection of chains and put it in meet at the same time removing it from c and cc
1051: std::set<point_type> intersect;
1052: //std::set_intersection(chain0->begin(), chain0->end(), chain1->begin(), chain1->end(), std::insert_iterator<supportSet>(join.obj(), join->begin()));
1053: std::set_intersection(chain0->begin(), chain0->end(), chain1->begin(), chain1->end(), std::insert_iterator<std::set<point_type> >(intersect, intersect.begin()));
1054: join->insert(intersect.begin(), intersect.end());
1055: // std::cout << " Join set:" << std::endl;
1056: // for(typename supportSet::iterator i_iter = join->begin(); i_iter != join->end(); ++i_iter) {
1057: // std::cout << " " << *i_iter << std::endl;
1058: // }
1059: for(typename std::set<point_type>::iterator i_iter = intersect.begin(); i_iter != intersect.end(); ++i_iter) {
1060: chain0->erase(chain0->find(*i_iter));
1061: chain1->erase(chain1->find(*i_iter));
1062: }
1063: // 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.
1064: support = this->support(chain0);
1065: chain0->insert(support->begin(), support->end());
1066: if(chain0->size() == 0) {
1067: break;
1068: }
1069: // std::cout << " chain 0:" << std::endl;
1070: // for(typename InputSequence::iterator i_iter = chain0->begin(); i_iter != chain0->end(); ++i_iter) {
1071: // std::cout << " " << *i_iter << std::endl;
1072: // }
1073: support = this->support(chain1);
1074: chain1->insert(support->begin(), support->end());
1075: if(chain1->size() == 0) {
1076: break;
1077: }
1078: // std::cout << " chain 1:" << std::endl;
1079: // for(typename InputSequence::iterator i_iter = chain1->begin(); i_iter != chain1->end(); ++i_iter) {
1080: // std::cout << " " << *i_iter << std::endl;
1081: // }
1082: // If both supports are empty, we should quit
1083: }
1084: }
1085: return join;
1086: };
1088: template <typename Point_, typename Marker_, typename Color_>
1089: template<typename ostream_type>
1090: void Sieve<Point_,Marker_,Color_>::view(ostream_type& os, const char* label, bool rawData){
1091: if(label != NULL) {
1092: os << "Viewing Sieve '" << label << "':" << std::endl;
1093: }
1094: else {
1095: os << "Viewing a Sieve:" << std::endl;
1096: }
1097: if(!rawData) {
1098: os << "cap --> base:" << std::endl;
1099: Obj<typename traits::capSequence> cap = this->cap();
1100: for(typename traits::capSequence::iterator capi = cap->begin(); capi != cap->end(); ++capi) {
1101: Obj<typename traits::supportSequence> supp = this->support(*capi);
1102: for(typename traits::supportSequence::iterator suppi = supp->begin(); suppi != supp->end(); ++suppi) {
1103: os << *capi << "--(" << suppi.color() << ")-->" << *suppi << std::endl;
1104: }
1105: }
1106: os << "base <-- cap:" << std::endl;
1107: Obj<typename traits::baseSequence> base = this->base();
1108: for(typename traits::baseSequence::iterator basei = base->begin(); basei != base->end(); ++basei) {
1109: Obj<typename traits::coneSequence> cone = this->cone(*basei);
1110: for(typename traits::coneSequence::iterator conei = cone->begin(); conei != cone->end(); ++conei) {
1111: os << *basei << "<--(" << conei.color() << ")--" << *conei << std::endl;
1112: }
1113: }
1114: #if 0
1115: os << "cap --> (outdegree, marker, depth, height):" << std::endl;
1116: for(typename traits::capSequence::iterator capi = cap->begin(); capi != cap->end(); ++capi) {
1117: os << *capi << "-->" << capi.degree() << ", " << capi.marker() << ", " << capi.depth() << ", " << capi.height() << std::endl;
1118: }
1119: os << "base --> (indegree, marker, depth, height):" << std::endl;
1120: for(typename traits::baseSequence::iterator basei = base->begin(); basei != base->end(); ++basei) {
1121: os << *basei << "-->" << basei.degree() << ", " << basei.marker() << ", " << basei.depth() << ", " << basei.height() << std::endl;
1122: }
1123: #endif
1124: }
1125: else {
1126: os << "'raw' arrow set:" << std::endl;
1127: for(typename traits::arrow_container_type::set_type::iterator ai = this->_arrows.set.begin(); ai != this->_arrows.set.end(); ai++)
1128: {
1129: typename traits::arrow_type arr = *ai;
1130: os << arr << std::endl;
1131: }
1132: os << "'raw' base set:" << std::endl;
1133: for(typename traits::base_container_type::set_type::iterator bi = this->_base.set.begin(); bi != this->_base.set.end(); bi++)
1134: {
1135: typename traits::base_container_type::traits::rec_type bp = *bi;
1136: os << bp << std::endl;
1137: }
1138: os << "'raw' cap set:" << std::endl;
1139: for(typename traits::cap_container_type::set_type::iterator ci = this->_cap.set.begin(); ci != this->_cap.set.end(); ci++)
1140: {
1141: typename traits::cap_container_type::traits::rec_type cp = *ci;
1142: os << cp << std::endl;
1143: }
1144: }
1145: };
1146: template <typename Point_, typename Marker_, typename Color_>
1147: void Sieve<Point_,Marker_,Color_>::view(const char* label, MPI_Comm comm) {
1148: ostringstream txt;
1150: if (this->debug()) {
1151: std::cout << "viewing a Sieve, comm = " << this->comm() << ", commRank = " << this->commRank() << std::endl;
1152: }
1153: if(label != NULL) {
1154: if(this->commRank() == 0) {
1155: txt << "viewing Sieve :'" << label << "'" << std::endl;
1156: }
1157: }
1158: else {
1159: if(this->commRank() == 0) {
1160: txt << "viewing a Sieve" << std::endl;
1161: }
1162: }
1163: if(this->commRank() == 0) {
1164: txt << "cap --> base:\n";
1165: }
1166: typename traits::capSequence cap = this->cap();
1167: typename traits::baseSequence base = this->base();
1168: if(cap.empty()) {
1169: txt << "[" << this->commRank() << "]: empty" << std::endl;
1170: }
1171: for(typename traits::capSequence::iterator capi = cap.begin(); capi != cap.end(); capi++) {
1172: const Obj<typename traits::supportSequence>& supp = this->support(*capi);
1173: const typename traits::supportSequence::iterator suppEnd = supp->end();
1175: for(typename traits::supportSequence::iterator suppi = supp->begin(); suppi != suppEnd; suppi++) {
1176: txt << "[" << this->commRank() << "]: " << *capi << "--(" << suppi.color() << ")-->" << *suppi << std::endl;
1177: }
1178: }
1179: PetscSynchronizedPrintf(this->comm(), txt.str().c_str());
1180: PetscSynchronizedFlush(this->comm());
1181: //
1182: ostringstream txt1;
1183: if(this->commRank() == 0) {
1184: txt1 << "base --> cap:\n";
1185: }
1186: if(base.empty()) {
1187: txt1 << "[" << this->commRank() << "]: empty" << std::endl;
1188: }
1189: for(typename traits::baseSequence::iterator basei = base.begin(); basei != base.end(); basei++) {
1190: const Obj<typename traits::coneSequence>& cone = this->cone(*basei);
1191: const typename traits::coneSequence::iterator coneEnd = cone->end();
1193: for(typename traits::coneSequence::iterator conei = cone->begin(); conei != coneEnd; conei++) {
1194: txt1 << "[" << this->commRank() << "]: " << *basei << "<--(" << conei.color() << ")--" << *conei << std::endl;
1195: }
1196: }
1197: //
1198: PetscSynchronizedPrintf(this->comm(), txt1.str().c_str());
1199: PetscSynchronizedFlush(this->comm());
1200: #if 0
1201: //
1202: ostringstream txt2;
1203: if(this->commRank() == 0) {
1204: txt2 << "cap <point, outdegree, marker, depth, height>:\n";
1205: }
1206: txt2 << "[" << this->commRank() << "]: [";
1207: for(typename traits::capSequence::iterator capi = cap.begin(); capi != cap.end(); capi++) {
1208: txt2 << " <" << *capi << ", " << capi.degree() << ", " << capi.marker() << ", " << capi.depth() << ", " << capi.height() << ">";
1209: }
1210: txt2 << " ]" << std::endl;
1211: //
1212: PetscSynchronizedPrintf(this->comm(), txt2.str().c_str());
1213: PetscSynchronizedFlush(this->comm());
1214: //
1215: ostringstream txt3;
1216: if(this->commRank() == 0) {
1217: txt3 << "base <point, indegree, marker, depth, height>:\n";
1218: }
1219: txt3 << "[" << this->commRank() << "]: [";
1220: for(typename traits::baseSequence::iterator basei = base.begin(); basei != base.end(); basei++) {
1221: txt3 << " <" << *basei << "," << basei.degree() << ", " << basei.marker() << ", " << basei.depth() << ", " << basei.height() << ">";
1222: }
1223: txt3 << " ]" << std::endl;
1224: //
1225: PetscSynchronizedPrintf(this->comm(), txt3.str().c_str());
1226: PetscSynchronizedFlush(this->comm());
1227: #endif
1228: };
1229: //
1230: // Structural queries
1231: //
1232: template <typename Point_, typename Marker_, typename Color_>
1233: void Sieve<Point_,Marker_,Color_>::setMarker(const point_type& p, const marker_type& marker) {
1234: 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);
1235: 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);
1237: if (bIndex.find(p) != bIndex.end()) {
1238: bIndex.modify(bIndex.find(p), changeMarker(marker));
1239: }
1240: if (cIndex.find(p) != cIndex.end()) {
1241: cIndex.modify(cIndex.find(p), changeMarker(marker));
1242: }
1243: this->_markers->insert(marker);
1244: };
1246: template <typename Point_, typename Marker_, typename Color_>
1247: template <typename Sequence>
1248: void Sieve<Point_,Marker_,Color_>::setMarker(const Obj<Sequence>& points, const marker_type& marker) {
1249: for(typename Sequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
1250: this->setMarker(*p_iter, marker);
1251: }
1252: };
1254: template <typename Point_, typename Marker_, typename Color_>
1255: int Sieve<Point_,Marker_,Color_>::depth() {
1256: return this->maxDepth;
1257: };
1258: template <typename Point_, typename Marker_, typename Color_>
1259: int Sieve<Point_,Marker_,Color_>::depth(const point_type& p) {
1260: 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);
1261: if (i.find(p) != i.end()) {
1262: return i.find(p)->depth;
1263: }
1264: return 0;
1265: };
1266: template <typename Point_, typename Marker_, typename Color_>
1267: template<typename InputSequence>
1268: int Sieve<Point_,Marker_,Color_>::depth(const Obj<InputSequence>& points) {
1269: 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);
1270: int maxDepth = 0;
1271:
1272: for(typename InputSequence::iterator iter = points->begin(); iter != points->end(); ++iter) {
1273: if (i.find(*iter) != i.end()) {
1274: maxDepth = std::max(maxDepth, i.find(*iter)->depth);
1275: }
1276: }
1277: return maxDepth;
1278: };
1279: template <typename Point_, typename Marker_, typename Color_>
1280: int Sieve<Point_,Marker_,Color_>::height() {
1281: return this->maxHeight;
1282: };
1283: template <typename Point_, typename Marker_, typename Color_>
1284: int Sieve<Point_,Marker_,Color_>::height(const point_type& p) {
1285: 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);
1286: if (i.find(p) != i.end()) {
1287: return i.find(p)->height;
1288: }
1289: return 0;
1290: };
1291: template <typename Point_, typename Marker_, typename Color_>
1292: template<typename InputSequence>
1293: int Sieve<Point_,Marker_,Color_>::height(const Obj<InputSequence>& points) {
1294: 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);
1295: int maxHeight = 0;
1296:
1297: for(typename InputSequence::iterator iter = points->begin(); iter != points->end(); ++iter) {
1298: if (i.find(*iter) != i.end()) {
1299: maxHeight = std::max(maxHeight, i.find(*iter)->height);
1300: }
1301: }
1302: return maxHeight;
1303: };
1305: template <typename Point_, typename Marker_, typename Color_>
1306: int Sieve<Point_,Marker_,Color_>::diameter() {
1307: int globalDiameter;
1308: int MPI_Allreduce(&this->graphDiameter, &globalDiameter, 1, MPI_INT, MPI_MAX, this->comm());
1309: CHKMPIERROR(ierr, ERRORMSG("Error in MPI_Allreduce"));
1310: return globalDiameter;
1311: };
1312: template <typename Point_, typename Marker_, typename Color_>
1313: int Sieve<Point_,Marker_,Color_>::diameter(const point_type& p) {
1314: return this->depth(p) + this->height(p);
1315: };
1317: template <typename Point_, typename Marker_, typename Color_>
1318: template<class InputSequence>
1319: void Sieve<Point_,Marker_,Color_>::__computeClosureHeights(const Obj<InputSequence>& points) {
1320: 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);
1321: 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);
1322: Obj<coneSet> modifiedPoints = coneSet();
1323:
1324: for(typename InputSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
1325: // Compute the max height of the points in the support of p, and add 1
1326: int h0 = this->height(*p_itor);
1327: int h1 = this->height(this->support(*p_itor)) + 1;
1328: if(h1 != h0) {
1329: 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);
1331: index.modify(index.find(*p_itor), changeHeight(h1));
1332: if (bIter != bIndex.end()) {
1333: bIndex.modify(bIter, changeHeight(h1));
1334: }
1335: if (h1 > this->maxHeight) this->maxHeight = h1;
1336: modifiedPoints->insert(*p_itor);
1337: }
1338: }
1339: // FIX: We would like to avoid the copy here with cone()
1340: if(modifiedPoints->size() > 0) {
1341: this->__computeClosureHeights(this->cone(modifiedPoints));
1342: }
1343: };
1344: template <typename Point_, typename Marker_, typename Color_>
1345: template<class InputSequence>
1346: void Sieve<Point_,Marker_,Color_>::__computeStarDepths(const Obj<InputSequence>& points) {
1347: 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);
1348: 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);
1349: Obj<supportSet> modifiedPoints = supportSet();
1350: for(typename InputSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
1351: // Compute the max depth of the points in the support of p, and add 1
1352: int d0 = this->depth(*p_itor);
1353: int d1 = this->depth(this->cone(*p_itor)) + 1;
1354: if(d1 != d0) {
1355: 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);
1357: index.modify(index.find(*p_itor), changeDepth(d1));
1358: if (cIter != cIndex.end()) {
1359: cIndex.modify(cIter, changeDepth(d1));
1360: }
1361: if (d1 > this->maxDepth) this->maxDepth = d1;
1362: modifiedPoints->insert(*p_itor);
1363: }
1364: }
1365: // FIX: We would like to avoid the copy here with cone()
1366: if(modifiedPoints->size() > 0) {
1367: this->__computeStarDepths(this->support(modifiedPoints));
1368: }
1369: };
1372: template <typename Point_, typename Marker_, typename Color_>
1373: void Sieve<Point_,Marker_,Color_>::stratify(bool show) {
1374: ALE_LOG_EVENT_BEGIN;
1375: // FIX: We would like to avoid the copy here with cone() and support()
1376: this->__computeClosureHeights(this->cone(this->leaves()));
1377: this->__computeStarDepths(this->support(this->roots()));
1378:
1379: Obj<typename traits::capSequence> base = this->base();
1381: for(typename traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1382: maxDepth = std::max(maxDepth, b_iter.depth());
1383: //b_iter.setDegree(this->cone(*b_iter)->size());
1384: this->_base.adjustDegree(*b_iter, this->cone(*b_iter)->size());
1385: }
1386: Obj<typename traits::capSequence> cap = this->cap();
1388: for(typename traits::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1389: maxHeight = std::max(maxHeight, c_iter.height());
1390: //c_iter.setDegree(this->support(*c_iter)->size());
1391: this->_cap.adjustDegree(*c_iter, this->support(*c_iter)->size());
1392: }
1393: if (this->debug() || show) {
1394: // const typename ::boost::multi_index::index<StratumSet,point>::type& points = ::boost::multi_index::get<point>(this->strata);
1395: // for(typename ::boost::multi_index::index<StratumSet,point>::type::iterator i = points.begin(); i != points.end(); i++) {
1396: // std::cout << *i << std::endl;
1397: // }
1398: }
1399: ALE_LOG_EVENT_END;
1400: };
1401: //
1402: // Structural manipulation
1403: //
1404:
1405: } // namespace ALE
1407: #endif