Actual source code: Sifter.hh

  1: #ifndef included_ALE_Sifter_hh
  2: #define included_ALE_Sifter_hh

  4: /*
  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: */
 10: #include <iostream>

 12: // ALE extensions

 14: #ifndef  included_ALE_hh
 15: #include <ALE.hh>
 16: #endif


 20: namespace ALE {

 22:   namespace SifterDef {
 23:     // Defines the traits of a sequence representing a subset of a multi_index container Index_.
 24:     // A sequence defines output (input in std terminology) iterators for traversing an Index_ object.
 25:     // Upon dereferencing values are extracted from each result record using a ValueExtractor_ object.
 26:     template <typename Index_, typename ValueExtractor_>
 27:     struct IndexSequenceTraits {
 28:       typedef Index_ index_type;
 29:       class iterator_base {
 30:       public:
 31:         // Standard iterator typedefs
 32:         typedef ValueExtractor_                        extractor_type;
 33:         typedef std::input_iterator_tag                iterator_category;
 34:         typedef typename extractor_type::result_type   value_type;
 35:         typedef int                                    difference_type;
 36:         typedef value_type*                            pointer;
 37:         typedef value_type&                            reference;
 38: 
 39:         // Underlying iterator type
 40:         typedef typename index_type::iterator          itor_type;
 41:       protected:
 42:         // Underlying iterator
 43:         itor_type      _itor;
 44:         // Member extractor
 45:         extractor_type _ex;
 46:       public:
 47:         iterator_base(itor_type itor) {
 48:           this->_itor = itor_type(itor);
 49:         };
 50:         virtual ~iterator_base() {};
 51:         virtual bool              operator==(const iterator_base& iter) const {return this->_itor == iter._itor;};
 52:         virtual bool              operator!=(const iterator_base& iter) const {return this->_itor != iter._itor;};
 53:         // FIX: operator*() should return a const reference, but it won't compile that way, because _ex() returns const value_type
 54:         virtual const value_type  operator*() const {return _ex(*(this->_itor));};
 55:       };// class iterator_base
 56:       class iterator : public iterator_base {
 57:       public:
 58:         // Standard iterator typedefs
 59:         typedef typename iterator_base::iterator_category  iterator_category;
 60:         typedef typename iterator_base::value_type         value_type;
 61:         typedef typename iterator_base::extractor_type     extractor_type;
 62:         typedef typename iterator_base::difference_type    difference_type;
 63:         typedef typename iterator_base::pointer            pointer;
 64:         typedef typename iterator_base::reference          reference;
 65:         // Underlying iterator type
 66:         typedef typename iterator_base::itor_type          itor_type;
 67:       public:
 68:         iterator(const itor_type& itor) : iterator_base(itor) {};
 69:         virtual ~iterator() {};
 70:         //
 71:         virtual iterator   operator++() {++this->_itor; return *this;};
 72:         virtual iterator   operator++(int n) {iterator tmp(this->_itor); ++this->_itor; return tmp;};
 73:       };// class iterator
 74:     }; // struct IndexSequenceTraits
 75: 
 76:     template <typename Index_, typename ValueExtractor_>
 77:     struct ReversibleIndexSequenceTraits {
 78:       typedef IndexSequenceTraits<Index_, ValueExtractor_> base_traits;
 79:       typedef typename base_traits::iterator_base   iterator_base;
 80:       typedef typename base_traits::iterator        iterator;
 81:       typedef typename base_traits::index_type      index_type;

 83:       // reverse_iterator is the reverse of iterator
 84:       class reverse_iterator : public iterator_base {
 85:       public:
 86:         // Standard iterator typedefs
 87:         typedef typename iterator_base::iterator_category  iterator_category;
 88:         typedef typename iterator_base::value_type         value_type;
 89:         typedef typename iterator_base::extractor_type     extractor_type;
 90:         typedef typename iterator_base::difference_type    difference_type;
 91:         typedef typename iterator_base::pointer            pointer;
 92:         typedef typename iterator_base::reference          reference;
 93:         // Underlying iterator type
 94:         typedef typename iterator_base::itor_type          itor_type;
 95:       public:
 96:         reverse_iterator(const itor_type& itor) : iterator_base(itor) {};
 97:         virtual ~reverse_iterator() {};
 98:         //
 99:         virtual reverse_iterator     operator++() {--this->_itor; return *this;};
100:         virtual reverse_iterator     operator++(int n) {reverse_iterator tmp(this->_itor); --this->_itor; return tmp;};
101:       };
102:     }; // class ReversibleIndexSequenceTraits


105:     //
106:     // Rec & RecContainer definitions.
107:     // Rec is intended to denote a graph point record.
108:     //
109:     template <typename Point_>
110:     struct Rec {
111:       typedef Point_ point_type;
112:       template<typename OtherPoint_>
113:       struct rebind {
114:         typedef Rec<OtherPoint_> type;
115:       };
116:       point_type     point;
117:       int            degree;
118:       // Basic interface
119:       Rec() : degree(0){};
120:       Rec(const Rec& r) : point(r.point), degree(r.degree) {}
121:       //Rec(const point_type& p) : point(p), degree(0) {};
122:       Rec(const point_type& p, const int d) : point(p), degree(d) {};
123:       // Printing
124:       friend std::ostream& operator<<(std::ostream& os, const Rec& p) {
125:         os << "<" << p.point << ", "<< p.degree << ">";
126:         return os;
127:       };
128: 
129:       struct degreeAdjuster {
130:         degreeAdjuster(int newDegree) : _newDegree(newDegree) {};
131:         void operator()(Rec& r) { r.degree = this->_newDegree; }
132:       private:
133:         int _newDegree;
134:       };// degreeAdjuster()

136:     };// class Rec

138:     template <typename Point_, typename Rec_>
139:     struct RecContainerTraits {
140:       typedef Rec_ rec_type;
141:       // Index tags
142:       struct pointTag{};
143:       // Rec set definition
144:       typedef ::boost::multi_index::multi_index_container<
145:         rec_type,
146:         ::boost::multi_index::indexed_by<
147:           ::boost::multi_index::ordered_unique<
148:             ::boost::multi_index::tag<pointTag>, BOOST_MULTI_INDEX_MEMBER(rec_type, typename rec_type::point_type, point)
149:           >
150:         >,
151:         ALE_ALLOCATOR<rec_type>
152:       > set_type;
153:       //
154:       // Return types
155:       //

157:      class PointSequence {
158:      public:
159:         typedef IndexSequenceTraits<typename ::boost::multi_index::index<set_type, pointTag>::type,
160:                                     BOOST_MULTI_INDEX_MEMBER(rec_type, typename rec_type::point_type,point)>
161:         traits;
162:       protected:
163:         const typename traits::index_type& _index;
164:       public:
165: 
166:        // Need to extend the inherited iterator to be able to extract the degree
167:        class iterator : public traits::iterator {
168:        public:
169:          iterator(const typename traits::iterator::itor_type& itor) : traits::iterator(itor) {};
170:          virtual const int& degree() const {return this->_itor->degree;};
171:        };
172: 
173:        PointSequence(const PointSequence& seq)            : _index(seq._index) {};
174:        PointSequence(typename traits::index_type& index) : _index(index)     {};
175:        virtual ~PointSequence(){};
176: 
177:        virtual bool empty(){return this->_index.empty();};
178: 
179:        virtual typename traits::index_type::size_type size() {return this->_index.size();};

181:        virtual iterator begin() {
182:          // Retrieve the beginning iterator of the index
183:          return iterator(this->_index.begin());
184:        };
185:        virtual iterator end() {
186:          // Retrieve the ending iterator of the index
187:          // Since the elements in this index are ordered by degree, this amounts to the end() of the index.
188:          return iterator(this->_index.end());
189:        };
190:        virtual bool contains(const typename rec_type::point_type& p) {
191:          // Check whether a given point is in the index
192:          return (this->_index.find(p) != this->_index.end());
193:        }
194:      }; // class PointSequence
195:     };// struct RecContainerTraits


198:     template <typename Point_, typename Rec_>
199:     struct RecContainer {
200:       typedef RecContainerTraits<Point_, Rec_> traits;
201:       typedef typename traits::set_type set_type;
202:       template <typename OtherPoint_, typename OtherRec_>
203:       struct rebind {
204:         typedef RecContainer<OtherPoint_, OtherRec_> type;
205:       };
206:       set_type set;
207:       //
208:       void removePoint(const typename traits::rec_type::point_type& p) {
209:         /*typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type& index = 
210:           ::boost::multi_index::get<typename traits::pointTag>(this->set);
211:         typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type::iterator i = index.find(p);
212:         if (i != index.end()) { // Point exists
213:           i = index.erase(i);
214:         }*/
215:         this->erase(p);
216:       };
217:       //
218:       void adjustDegree(const typename traits::rec_type::point_type& p, int delta) {
219:         typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type& index =
220:           ::boost::multi_index::get<typename traits::pointTag>(this->set);
221:         typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type::iterator i = index.find(p);
222:         if (i == index.end()) { // No such point exists
223:           if(delta < 0) { // Cannot decrease degree of a non-existent point
224:             ostringstream err;
225:             err << "ERROR: adjustDegree: Non-existent point " << p;
226:             std::cout << err << std::endl;
227:             throw(Exception(err.str().c_str()));
228:           }
229:           else { // We CAN INCREASE the degree of a non-existent point: simply insert a new element with degree == delta
230:             std::pair<typename ::boost::multi_index::index<set_type, typename traits::pointTag>::type::iterator, bool> ii;
231:             typename traits::rec_type r(p,delta);
232:             ii = index.insert(r);
233:             if(ii.second == false) {
234:               ostringstream err;
235:               err << "ERROR: adjustDegree: Failed to insert a rec " << r;
236:               std::cout << err << std::endl;
237:               throw(Exception(err.str().c_str()));
238:             }
239:           }
240:         }
241:         else { // Point exists, so we try to modify its degree
242:           // If the adjustment is zero, there is nothing to do, otherwise ...
243:           if(delta != 0) {
244:             int newDegree = i->degree + delta;
245:             if(newDegree < 0) {
246:               ostringstream ss;
247:               ss << "adjustDegree: Adjustment of " << *i << " by " << delta << " would result in negative degree: " << newDegree;
248:               throw Exception(ss.str().c_str());
249:             }
250:             index.modify(i, typename traits::rec_type::degreeAdjuster(newDegree));
251:           }
252:         }
253:       }; // adjustDegree()
254:     }; // struct RecContainer

256:     //
257:     // Arrow & ArrowContainer definitions
258:     //
259:     template<typename Source_, typename Target_, typename Color_>
260:     struct  Arrow { //: public ALE::def::Arrow<Source_, Target_, Color_> {
261:       typedef Arrow   arrow_type;
262:       typedef Source_ source_type;
263:       typedef Target_ target_type;
264:       typedef Color_  color_type;
265:       source_type source;
266:       target_type target;
267:       color_type  color;
268:       Arrow(const source_type& s, const target_type& t, const color_type& c) : source(s), target(t), color(c) {};
269:       // Flipping
270:       template <typename OtherSource_, typename OtherTarget_, typename OtherColor_>
271:       struct rebind {
272:         typedef Arrow<OtherSource_, OtherTarget_, OtherColor_> type;
273:       };
274:       struct flip {
275:         typedef Arrow<target_type, source_type, color_type> type;
276:         type arrow(const arrow_type& a) { return type(a.target, a.source, a.color);};
277:       };

279:       // Printing
280:       friend std::ostream& operator<<(std::ostream& os, const Arrow& a) {
281:         os << a.source << " --(" << a.color << ")--> " << a.target;
282:         return os;
283:       }

285:       // Arrow modifiers
286:       struct sourceChanger {
287:         sourceChanger(const source_type& newSource) : _newSource(newSource) {};
288:         void operator()(arrow_type& a) {a.source = this->_newSource;}
289:       private:
290:         source_type _newSource;
291:       };

293:       struct targetChanger {
294:         targetChanger(const target_type& newTarget) : _newTarget(newTarget) {};
295:         void operator()(arrow_type& a) { a.target = this->_newTarget;}
296:       private:
297:         const target_type _newTarget;
298:       };
299:     };// struct Arrow
300: 

302:     template<typename Source_, typename Target_, typename Color_, typename SupportCompare_>
303:     struct ArrowContainerTraits {
304:     public:
305:       //
306:       // Encapsulated types
307:       //
308:       typedef Arrow<Source_,Target_,Color_>    arrow_type;
309:       typedef typename arrow_type::source_type source_type;
310:       typedef typename arrow_type::target_type target_type;
311:       typedef typename arrow_type::color_type  color_type;
312:       typedef SupportCompare_                  support_compare_type;
313:       // Index tags
314:       struct                                   sourceColorTag{};
315:       struct                                   targetColorTag{};
316:       struct                                   sourceTargetTag{};

318:       // Sequence traits and sequence types
319:       template <typename Index_, typename Key_, typename SubKey_, typename ValueExtractor_>
320:       class ArrowSequence {
321:         // ArrowSequence implements ReversibleIndexSequencTraits with Index_ and ValueExtractor_ types.
322:         // A Key_ object and an optional SubKey_ object are used to extract the index subset.
323:       public:
324:         typedef ReversibleIndexSequenceTraits<Index_, ValueExtractor_>  traits;
325:         //typedef source_type                                             source_type;
326:         //typedef target_type                                             target_type;
327:         //typedef arrow_type                                              arrow_type;
328:         //
329:         typedef Key_                                                    key_type;
330:         typedef SubKey_                                                 subkey_type;
331:       protected:
332:         typename traits::index_type&                                    _index;
333:         key_type                                                  key;
334:         subkey_type                                               subkey;
335:         bool                                                      useSubkey;
336:       public:
337:         // Need to extend the inherited iterators to be able to extract arrow color
338:         class iterator : public traits::iterator {
339:         public:
340:           iterator(const typename traits::iterator::itor_type& itor) : traits::iterator(itor) {};
341:           virtual const source_type& source() const {return this->_itor->source;};
342:           virtual const color_type&  color()  const {return this->_itor->color;};
343:           virtual const target_type& target() const {return this->_itor->target;};
344:           virtual const arrow_type&  arrow()  const {return *(this->_itor);};
345:         };
346:         class reverse_iterator : public traits::reverse_iterator {
347:         public:
348:           reverse_iterator(const typename traits::reverse_iterator::itor_type& itor) : traits::reverse_iterator(itor) {};
349:           virtual const source_type& source() const {return this->_itor->source;};
350:           virtual const color_type&  color()  const {return this->_itor->color;};
351:           virtual const target_type& target() const {return this->_itor->target;};
352:           virtual const arrow_type&  arrow()  const {return *(this->_itor);};
353:         };
354:       public:
355:         //
356:         // Basic ArrowSequence interface
357:         //
358:         ArrowSequence(const ArrowSequence& seq) : _index(seq._index), key(seq.key), subkey(seq.subkey), useSubkey(seq.useSubkey) {};
359:         ArrowSequence(typename traits::index_type& index, const key_type& k) :
360:           _index(index), key(k), subkey(subkey_type()), useSubkey(0) {};
361:         ArrowSequence(typename traits::index_type& index, const key_type& k, const subkey_type& kk) :
362:           _index(index), key(k), subkey(kk), useSubkey(1){};
363:         virtual ~ArrowSequence() {};

365:         void setKey(const key_type& key) {this->key = key;};
366:         void setSubkey(const subkey_type& subkey) {this->subkey = subkey;};
367:         void setUseSubkey(const bool& useSubkey) {this->useSubkey = useSubkey;};
368: 
369:         virtual bool         empty() {return this->_index.empty();};

371:         virtual typename traits::index_type::size_type  size()  {
372:           if (this->useSubkey) {
373:             return this->_index.count(::boost::make_tuple(this->key,this->subkey));
374:           } else {
375:             return this->_index.count(::boost::make_tuple(this->key));
376:           }
377:         };

379:         virtual iterator begin() {
380:           if (this->useSubkey) {
381:             return iterator(this->_index.lower_bound(::boost::make_tuple(this->key,this->subkey)));
382:           } else {
383:             return iterator(this->_index.lower_bound(::boost::make_tuple(this->key)));
384:           }
385:         };
386: 
387:         virtual iterator end() {
388:           if (this->useSubkey) {
389:             return iterator(this->_index.upper_bound(::boost::make_tuple(this->key,this->subkey)));
390:           } else {
391:             return iterator(this->_index.upper_bound(::boost::make_tuple(this->key)));
392:           }
393:         };
394: 
395:         virtual reverse_iterator rbegin() {
396:           if (this->useSubkey) {
397:             return reverse_iterator(--this->_index.upper_bound(::boost::make_tuple(this->key,this->subkey)));
398:           } else {
399:             return reverse_iterator(--this->_index.upper_bound(::boost::make_tuple(this->key)));
400:           }
401:         };
402: 
403:         virtual reverse_iterator rend() {
404:           if (this->useSubkey) {
405:             return reverse_iterator(--this->_index.lower_bound(::boost::make_tuple(this->key,this->subkey)));
406:           } else {
407:             return reverse_iterator(--this->_index.lower_bound(::boost::make_tuple(this->key)));
408:           }
409:         };

411:         template<typename ostream_type>
412:         void view(ostream_type& os, const bool& useColor = false, const char* label = NULL){
413:           if(label != NULL) {
414:             os << "Viewing " << label << " sequence:" << std::endl;
415:           }
416:           os << "[";
417:           for(iterator i = this->begin(); i != this->end(); i++) {
418:             os << " (" << *i;
419:             if(useColor) {
420:               os << "," << i.color();
421:             }
422:             os  << ")";
423:           }
424:           os << " ]" << std::endl;
425:         };
426:       };// class ArrowSequence
427:     };// class ArrowContainerTraits
428: 

430:     // The specialized ArrowContainer types distinguish the cases of unique and multiple colors of arrows on
431:     // for each (source,target) pair (i.e., a single arrow, or multiple arrows between each pair of points).
432:     typedef enum {multiColor, uniColor} ColorMultiplicity;

434:     template<typename Source_, typename Target_, typename Color_, ColorMultiplicity colorMultiplicity, typename SupportCompare_>
435:     struct ArrowContainer {};
436: 
437:     template<typename Source_, typename Target_, typename Color_, typename SupportCompare_>
438:     struct ArrowContainer<Source_, Target_, Color_, multiColor, SupportCompare_> {
439:       // Define container's encapsulated types
440:       typedef ArrowContainerTraits<Source_, Target_, Color_, SupportCompare_>      traits;
441:       // need to def arrow_type locally, since BOOST_MULTI_INDEX_MEMBER barfs when first template parameter starts with 'typename'
442:       typedef typename traits::arrow_type                         arrow_type;
443:       // Container set type
444:       typedef ::boost::multi_index::multi_index_container<
445:         typename traits::arrow_type,
446:         ::boost::multi_index::indexed_by<
447:           ::boost::multi_index::ordered_non_unique<
448:             ::boost::multi_index::tag<typename traits::sourceTargetTag>,
449:             ::boost::multi_index::composite_key<
450:               typename traits::arrow_type,
451:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::source_type, source),
452:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::target_type, target),
453:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::color_type,  color)
454:             >
455:           >,
456:           ::boost::multi_index::ordered_non_unique<
457:             ::boost::multi_index::tag<typename traits::sourceColorTag>,
458:             ::boost::multi_index::composite_key<
459:               typename traits::arrow_type,
460:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::source_type, source),
461:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::color_type,  color),
462:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::target_type, target)
463:             >
464:           >,
465:           ::boost::multi_index::ordered_non_unique<
466:             ::boost::multi_index::tag<typename traits::targetColorTag>,
467:             ::boost::multi_index::composite_key<
468:               typename traits::arrow_type,
469:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::target_type, target),
470:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::color_type,  color),
471:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::source_type, source)
472:             >
473:           >
474:         >,
475:         ALE_ALLOCATOR<typename traits::arrow_type>
476:       > set_type;
477:      // multi-index set of multicolor arrows
478:       set_type set;
479:     }; // class ArrowContainer<multiColor>
480: 
481:     template<typename Source_, typename Target_, typename Color_, typename SupportCompare_>
482:     struct ArrowContainer<Source_, Target_, Color_, uniColor, SupportCompare_> {
483:       // Define container's encapsulated types
484:       typedef ArrowContainerTraits<Source_, Target_, Color_, SupportCompare_> traits;
485:       // need to def arrow_type locally, since BOOST_MULTI_INDEX_MEMBER barfs when first template parameter starts with 'typename'
486:       typedef typename traits::arrow_type                                   arrow_type;

488:       // multi-index set type -- arrow set
489:       typedef ::boost::multi_index::multi_index_container<
490:         typename traits::arrow_type,
491:         ::boost::multi_index::indexed_by<
492:           ::boost::multi_index::ordered_unique<
493:             ::boost::multi_index::tag<typename traits::sourceTargetTag>,
494:             ::boost::multi_index::composite_key<
495:               typename traits::arrow_type,
496:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::source_type, source),
497:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::target_type, target),
498:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::color_type,  color)
499:             >
500:           >,
501:           ::boost::multi_index::ordered_non_unique<
502:             ::boost::multi_index::tag<typename traits::sourceColorTag>,
503:             ::boost::multi_index::composite_key<
504:               typename traits::arrow_type,
505:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::source_type, source),
506:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::color_type,  color),
507:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::target_type, target)
508:             >,
509:             SupportCompare_
510:           >,
511:           ::boost::multi_index::ordered_non_unique<
512:             ::boost::multi_index::tag<typename traits::targetColorTag>,
513:             ::boost::multi_index::composite_key<
514:               typename traits::arrow_type,
515:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::target_type, target),
516:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::color_type,  color),
517:               BOOST_MULTI_INDEX_MEMBER(arrow_type, typename traits::source_type, source)
518:             >
519:           >
520:         >,
521:         ALE_ALLOCATOR<typename traits::arrow_type>
522:       > set_type;
523:       // multi-index set of unicolor arrow records
524:       set_type set;
525:     }; // class ArrowContainer<uniColor>
526:   }; // namespace SifterDef

528:   //
529:   // ASifter (short for Abstract Sifter, structurally a bipartite graph with colored arrows) implements a sequential interface
530:   // similar to that of Sieve, except the source and target points may have different types and iterated operations (e.g., nCone,
531:   // closure) are not available.
532:   //
533: template<typename Source_, typename Target_, typename Color_, SifterDef::ColorMultiplicity colorMultiplicity, typename SupportCompare_ = ::boost::multi_index::composite_key_compare<std::less<Source_>, std::less<Color_>, std::less<Target_> >, typename SourceCtnr_ = SifterDef::RecContainer<Source_, SifterDef::Rec<Source_> >, typename TargetCtnr_ = SifterDef::RecContainer<Target_, SifterDef::Rec<Target_> > >
534:   class ASifter { // class ASifter
535:   public:
536:     typedef struct {
537:       typedef ASifter<Source_, Target_, Color_, colorMultiplicity, SupportCompare_, SourceCtnr_, TargetCtnr_> graph_type;
538:       // Encapsulated container types
539:       typedef SifterDef::ArrowContainer<Source_, Target_, Color_, colorMultiplicity, SupportCompare_> arrow_container_type;
540:       typedef SourceCtnr_                                                            cap_container_type;
541:       typedef TargetCtnr_                                                            base_container_type;
542:       // Types associated with records held in containers
543:       typedef typename arrow_container_type::traits::arrow_type                      arrow_type;
544:       typedef typename arrow_container_type::traits::source_type                     source_type;
545:       typedef typename cap_container_type::traits::rec_type                          sourceRec_type;
546:       typedef typename arrow_container_type::traits::target_type                     target_type;
547:       typedef typename base_container_type::traits::rec_type                         targetRec_type;
548:       typedef typename arrow_container_type::traits::color_type                      color_type;
549:       typedef typename arrow_container_type::traits::support_compare_type            support_compare_type;
550:       // Convenient tag names
551:       typedef typename arrow_container_type::traits::sourceColorTag                  supportInd;
552:       typedef typename arrow_container_type::traits::targetColorTag                  coneInd;
553:       typedef typename arrow_container_type::traits::sourceTargetTag                 arrowInd;
554:       typedef typename base_container_type::traits::pointTag                         baseInd;
555:       typedef typename cap_container_type::traits::pointTag                          capInd;
556:       //
557:       // Return types
558:       //
559:       typedef typename
560:       arrow_container_type::traits::template ArrowSequence<typename ::boost::multi_index::index<typename arrow_container_type::set_type,arrowInd>::type, source_type, target_type, BOOST_MULTI_INDEX_MEMBER(arrow_type, color_type, color)>
561:       arrowSequence;

563:       // FIX: This is a temp fix to include addArrow into the interface; should probably be pushed up to ArrowSequence
564:       struct coneSequence : public arrow_container_type::traits::template ArrowSequence<typename ::boost::multi_index::index<typename arrow_container_type::set_type,coneInd>::type, target_type, color_type, BOOST_MULTI_INDEX_MEMBER(arrow_type, source_type, source)> {
565:       protected:
566:         graph_type& _graph;
567:       public:
568:         typedef typename
569:           arrow_container_type::traits::template ArrowSequence<typename ::boost::multi_index::index<typename arrow_container_type::set_type,coneInd>::type, target_type, color_type, BOOST_MULTI_INDEX_MEMBER(arrow_type, source_type, source)> base_type;
570:         // Encapsulated types
571:         typedef typename base_type::traits traits;
572:         typedef typename base_type::iterator iterator;
573:         typedef typename base_type::reverse_iterator reverse_iterator;
574:         // Basic interface
575:         coneSequence(const coneSequence& seq) : base_type(seq), _graph(seq._graph) {};
576:           coneSequence(graph_type& graph, typename traits::index_type& index, const typename base_type::key_type& k) : base_type(index, k), _graph(graph){};
577:             coneSequence(graph_type& graph, typename traits::index_type& index, const typename base_type::key_type& k, const typename base_type::subkey_type& kk) : base_type(index, k, kk), _graph(graph) {};
578:               virtual ~coneSequence() {};
579: 
580:         // Extended interface
581:         void addArrow(const arrow_type& a) {
582:           // if(a.target != this->key) {
583:           //               throw ALE::Exception("Arrow target mismatch in a coneSequence");
584:           //             }
585:           this->_graph.addArrow(a);
586:         };
587:         void addArrow(const source_type& s, const color_type& c){
588:           this->_graph.addArrow(arrow_type(s,this->key,c));
589:         };
590: 
591:         virtual bool contains(const source_type& s) {
592:           // Check whether a given point is in the index
593:           typename ::boost::multi_index::index<typename ASifter::traits::arrow_container_type::set_type,typename ASifter::traits::arrowInd>::type& index = ::boost::multi_index::get<typename ASifter::traits::arrowInd>(this->_graph._arrows.set);
594:           return (index.find(::boost::make_tuple(s,this->key)) != index.end());
595:         };
596:       };// struct coneSequence
597: 
598:       // FIX: This is a temp fix to include addArrow into the interface; should probably be pushed up to ArrowSequence
599:       struct supportSequence : public arrow_container_type::traits::template ArrowSequence<typename ::boost::multi_index::index<typename arrow_container_type::set_type,supportInd>::type, source_type, color_type, BOOST_MULTI_INDEX_MEMBER(arrow_type, target_type, target)> {
600:       protected:
601:         graph_type& _graph;
602:       public:
603:         typedef typename
604:           arrow_container_type::traits::template ArrowSequence<typename ::boost::multi_index::index<typename arrow_container_type::set_type,supportInd>::type, source_type, color_type, BOOST_MULTI_INDEX_MEMBER(arrow_type, target_type, target)> base_type;
605:         // Encapsulated types
606:         typedef typename base_type::traits traits;
607:         typedef typename base_type::iterator iterator;
608:         typedef typename base_type::reverse_iterator reverse_iterator;
609:         // Basic interface
610:         supportSequence(const supportSequence& seq) : base_type(seq), _graph(seq._graph) {};
611:         supportSequence(graph_type& graph, typename traits::index_type& index, const typename base_type::key_type& k) : base_type(index, k), _graph(graph){};
612:         supportSequence(graph_type& graph, typename traits::index_type& index, const typename base_type::key_type& k, const typename base_type::subkey_type& kk) : base_type(index, k, kk), _graph(graph) {};
613:         virtual ~supportSequence() {};
614: 
615:         // FIX: WARNING: (or a HACK?): we flip the arrow on addition here.
616:         // Fancy interface
617:         void addArrow(const typename arrow_type::flip::type& af) {
618:           this->_graph.addArrow(af.target, af.source, af.color);
619:         };
620:         void addArrow(const target_type& t, const color_type& c){
621:           this->_graph.addArrow(arrow_type(this->key,t,c));
622:         };
623:       };// struct supportSequence

625: 
626:       typedef typename base_container_type::traits::PointSequence baseSequence;
627:       typedef typename cap_container_type::traits::PointSequence  capSequence;
628:       typedef std::set<source_type>   coneSet;
629:       typedef ALE::array<source_type> coneArray;
630:       typedef std::set<target_type>   supportSet;
631:       typedef ALE::array<target_type> supportArray;
632:     } traits;

634:     template <typename OtherSource_, typename OtherTarget_, typename OtherColor_, SifterDef::ColorMultiplicity otherColorMultiplicity,
635:               typename OtherSupportCompare_  = ::boost::multi_index::composite_key_compare<std::less<OtherSource_>, std::less<OtherColor_>, std::less<OtherTarget_> >,
636:               typename OtherSourceCtnr_ = SifterDef::RecContainer<OtherSource_, SifterDef::Rec<OtherSource_> >,
637:               typename OtherTargetCtnr_ = SifterDef::RecContainer<OtherTarget_, SifterDef::Rec<OtherTarget_> > >
638:     struct rebind {
639:       typedef ASifter<OtherSource_, OtherTarget_, OtherColor_, otherColorMultiplicity, OtherSupportCompare_, OtherSourceCtnr_, OtherTargetCtnr_> type;
640:     };

642:   public:
643:     // Debug level
644:     int _debug;
645:     //protected:
646:     typename traits::arrow_container_type _arrows;
647:     typename traits::base_container_type  _base;
648:     typename traits::cap_container_type   _cap;
649:   protected:
650:     MPI_Comm    _comm;
651:     int         _commRank;
652:     int         _commSize;
653:     PetscObject _petscObj;
654:     void __init(MPI_Comm comm) {
655:       static PetscCookie sifterType = -1;
656:       //const char        *id_name = ALE::getClassName<T>();
657:       const char        *id_name = "Sifter";
658:       PetscErrorCode     ierr;

660:       if (sifterType < 0) {
661:         PetscCookieRegister(id_name,&sifterType);CHKERROR(ierr, "Error in MPI_Comm_rank");
662:       }
663:       this->_comm = comm;
664:       MPI_Comm_rank(this->_comm, &this->_commRank); CHKERROR(ierr, "Error in MPI_Comm_rank");
665:       MPI_Comm_size(this->_comm, &this->_commSize); CHKERROR(ierr, "Error in MPI_Comm_rank");
666: #ifdef USE_PETSC_OBJ
667:       PetscObjectCreateGeneric(this->_comm, sifterType, id_name, &this->_petscObj);CHKERROR(ierr, "Error in PetscObjectCreate");
668: #endif
669:       //ALE::restoreClassName<T>(id_name);
670:     };
671:     // We store these sequence objects to avoid creating them each query
672:     Obj<typename traits::coneSequence> _coneSeq;
673:     Obj<typename traits::supportSequence> _supportSeq;
674:   public:
675:     //
676:     // Basic interface
677:     //
678:     ASifter(MPI_Comm comm = PETSC_COMM_SELF, const int& debug = 0) : _debug(debug), _petscObj(NULL) {
679:       __init(comm);
680:       this->_coneSeq    = new typename traits::coneSequence(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), typename traits::target_type());
681:       this->_supportSeq = new typename traits::supportSequence(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), typename traits::source_type());
682:    }
683:     virtual ~ASifter() {
684: #ifdef USE_PETSC_OBJ
685:       if (this->_petscObj) {
687:         PetscObjectDestroy(this->_petscObj); CHKERROR(ierr, "Failed in PetscObjectDestroy");
688:         this->_petscObj = NULL;
689:       }
690: #endif
691:     };
692:     //
693:     // Query methods
694:     //
695:     int         debug()    const {return this->_debug;};
696:     void        setDebug(const int debug) {this->_debug = debug;};
697:     MPI_Comm    comm()     const {return this->_comm;};
698:     int         commSize() const {return this->_commSize;};
699:     int         commRank() const {return this->_commRank;}
700: #ifdef USE_PETSC_OBJ
701:     PetscObject petscObj() const {return this->_petscObj;};
702: #endif

704:     // FIX: need const_cap, const_base returning const capSequence etc, but those need to have const_iterators, const_begin etc.
705:     Obj<typename traits::capSequence> cap() {
706:       return typename traits::capSequence(::boost::multi_index::get<typename traits::capInd>(this->_cap.set));
707:     };
708:     Obj<typename traits::baseSequence> base() {
709:       return typename traits::baseSequence(::boost::multi_index::get<typename traits::baseInd>(this->_base.set));
710:     };
711:     bool capContains(const typename traits::source_type& p) {
712:       typename traits::capSequence cap(::boost::multi_index::get<typename traits::capInd>(this->_cap.set));

714:       //for(typename traits::capSequence::iterator c_iter = cap.begin(); c_iter != cap.end(); ++c_iter) {
715:       //}
716:       return cap.contains(p);
717:     };
718:     bool baseContains(const typename traits::target_type& p) {
719:       typename traits::baseSequence base(::boost::multi_index::get<typename traits::baseInd>(this->_base.set));

721:       //for(typename traits::capSequence::iterator c_iter = cap.begin(); c_iter != cap.end(); ++c_iter) {
722:       //}
723:       return base.contains(p);
724:     };
725:     // FIX: should probably have cone and const_cone etc, since arrows can be modified through an iterator (modifyColor).
726:     Obj<typename traits::arrowSequence>
727:     arrows(const typename traits::source_type& s, const typename traits::target_type& t) {
728:       return typename traits::arrowSequence(::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set), s, t);
729:     };
730:     Obj<typename traits::arrowSequence>
731:     arrows(const typename traits::source_type& s) {
732:       return typename traits::arrowSequence(::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set), s);
733:     };
734: #ifdef SLOW
735:     Obj<typename traits::coneSequence>
736:     cone(const typename traits::target_type& p) {
737:       return typename traits::coneSequence(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p);
738:     };
739: #else
740:     const Obj<typename traits::coneSequence>&
741:     cone(const typename traits::target_type& p) {
742:       this->_coneSeq->setKey(p);
743:       this->_coneSeq->setUseSubkey(false);
744:       return this->_coneSeq;
745:     };
746: #endif
747:     template<class InputSequence>
748:     Obj<typename traits::coneSet>
749:     cone(const Obj<InputSequence>& points) {
750:       return this->cone(points, typename traits::color_type(), false);
751:     };
752: #ifdef SLOW
753:     Obj<typename traits::coneSequence>
754:     cone(const typename traits::target_type& p, const typename traits::color_type& color) {
755:       return typename traits::coneSequence(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p, color);
756:     };
757: #else
758:     const Obj<typename traits::coneSequence>&
759:     cone(const typename traits::target_type& p, const typename traits::color_type& color) {
760:       this->_coneSeq->setKey(p);
761:       this->_coneSeq->setSubkey(color);
762:       this->_coneSeq->setUseSubkey(true);
763:       return this->_coneSeq;
764:     };
765: #endif
766:     template<class InputSequence>
767:     Obj<typename traits::coneSet>
768:     cone(const Obj<InputSequence>& points, const typename traits::color_type& color, bool useColor = true) {
769:       Obj<typename traits::coneSet> cone = typename traits::coneSet();
770:       for(typename InputSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
771:         Obj<typename traits::coneSequence> pCone;
772:         if (useColor) {
773:           pCone = this->cone(*p_itor, color);
774:         } else {
775:           pCone = this->cone(*p_itor);
776:         }
777:         cone->insert(pCone->begin(), pCone->end());
778:       }
779:       return cone;
780:     };
781:     template<typename PointCheck>
782:     bool coneContains(const typename traits::target_type& p, const PointCheck& checker) {
783:       typename traits::coneSequence cone(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p);

785:       for(typename traits::coneSequence::iterator c_iter = cone.begin(); c_iter != cone.end(); ++c_iter) {
786:         if (checker(*c_iter, p)) return true;
787:       }
788:       return false;
789:     };
790:     template<typename PointProcess>
791:     void coneApply(const typename traits::target_type& p, PointProcess& processor) {
792:       typename traits::coneSequence cone(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p);

794:       for(typename traits::coneSequence::iterator c_iter = cone.begin(); c_iter != cone.end(); ++c_iter) {
795:         processor(*c_iter, p);
796:       }
797:     };
798:     int getConeSize(const typename traits::target_type& p) {
799:       return this->cone(p)->size();
800:     }
801: #ifdef SLOW
802:     Obj<typename traits::supportSequence>
803:     support(const typename traits::source_type& p) {
804:       return typename traits::supportSequence(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p);
805:     };
806: #else
807:     const Obj<typename traits::supportSequence>&
808:     support(const typename traits::source_type& p) {
809:       this->_supportSeq->setKey(p);
810:       this->_supportSeq->setUseSubkey(false);
811:       return this->_supportSeq;
812:     };
813: #endif
814: #ifdef SLOW
815:     Obj<typename traits::supportSequence>
816:     support(const typename traits::source_type& p, const typename traits::color_type& color) {
817:       return typename traits::supportSequence(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p, color);
818:     };
819: #else
820:     const Obj<typename traits::supportSequence>&
821:     support(const typename traits::source_type& p, const typename traits::color_type& color) {
822:       this->_supportSeq->setKey(p);
823:       this->_supportSeq->setSubkey(color);
824:       this->_supportSeq->setUseSubkey(true);
825:       return this->_supportSeq;
826:     };
827: #endif
828:     template<class InputSequence>
829:     Obj<typename traits::supportSet>
830:     support(const Obj<InputSequence>& sources) {
831:       return this->support(sources, typename traits::color_type(), false);
832:     };
833:     template<class InputSequence>
834:     Obj<typename traits::supportSet>
835:     support(const Obj<InputSequence>& points, const typename traits::color_type& color, bool useColor = true){
836:       Obj<typename traits::supportSet> supp = typename traits::supportSet();
837:       for(typename InputSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
838:         Obj<typename traits::supportSequence> pSupport;
839:         if (useColor) {
840:           pSupport = this->support(*p_itor, color);
841:         } else {
842:           pSupport = this->support(*p_itor);
843:         }
844:         supp->insert(pSupport->begin(), pSupport->end());
845:       }
846:       return supp;
847:     };
848:     template<typename PointCheck>
849:     bool supportContains(const typename traits::source_type& p, const PointCheck& checker) {
850:       typename traits::supportSequence support(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p);

852:       for(typename traits::supportSequence::iterator s_iter = support.begin(); s_iter != support.end(); ++s_iter) {
853:         if (checker(*s_iter, p)) return true;
854:       }
855:       return false;
856:     };
857:     template<typename PointProcess>
858:     void supportApply(const typename traits::source_type& p, PointProcess& processor) {
859:       typename traits::supportSequence support(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p);

861:       for(typename traits::supportSequence::iterator s_iter = support.begin(); s_iter != support.end(); ++s_iter) {
862:         processor(*s_iter, p);
863:       }
864:     };
865:     int getSupportSize(const typename traits::source_type& p) {
866:       return this->support(p)->size();
867:     }

869:     template<typename ostream_type>
870:     void view(ostream_type& os, const char* label = NULL, bool rawData = false){
871:       int rank = this->commRank();

873:       if(label != NULL) {
874:         os << "["<<rank<<"]Viewing Sifter '" << label << "':" << std::endl;
875:       }
876:       else {
877:         os << "["<<rank<<"]Viewing a Sifter:" << std::endl;
878:       }
879:       if(!rawData) {
880:         os << "cap --> base:" << std::endl;
881:         Obj<typename traits::capSequence> cap = this->cap();
882:         for(typename traits::capSequence::iterator capi = cap->begin(); capi != cap->end(); capi++) {
883:           const Obj<typename traits::supportSequence>& supp = this->support(*capi);

885:           for(typename traits::supportSequence::iterator suppi = supp->begin(); suppi != supp->end(); suppi++) {
886:             os << *capi << "--(" << suppi.color() << ")-->" << *suppi << std::endl;
887:           }
888:         }
889:         os << "base <-- cap:" << std::endl;
890:         Obj<typename traits::baseSequence> base = this->base();
891:         for(typename traits::baseSequence::iterator basei = base->begin(); basei != base->end(); basei++) {
892:           const Obj<typename traits::coneSequence>& cone = this->cone(*basei);

894:           for(typename traits::coneSequence::iterator conei = cone->begin(); conei != cone->end(); conei++) {
895:             os << *basei <<  "<--(" << conei.color() << ")--" << *conei << std::endl;
896:           }
897:         }
898:         os << "cap --> outdegrees:" << std::endl;
899:         for(typename traits::capSequence::iterator capi = cap->begin(); capi != cap->end(); capi++) {
900:           os << *capi <<  "-->" << capi.degree() << std::endl;
901:         }
902:         os << "base <-- indegrees:" << std::endl;
903:         for(typename traits::baseSequence::iterator basei = base->begin(); basei != base->end(); basei++) {
904:           os << *basei <<  "<--" << basei.degree() << std::endl;
905:         }
906:       }
907:       else {
908:         os << "'raw' arrow set:" << std::endl;
909:         for(typename traits::arrow_container_type::set_type::iterator ai = _arrows.set.begin(); ai != _arrows.set.end(); ai++)
910:         {
911:           typename traits::arrow_type arr = *ai;
912:           os << arr << std::endl;
913:         }
914:         os << "'raw' base set:" << std::endl;
915:         for(typename traits::base_container_type::set_type::iterator bi = _base.set.begin(); bi != _base.set.end(); bi++)
916:         {
917:           typename traits::base_container_type::traits::rec_type bp = *bi;
918:           os << bp << std::endl;
919:         }
920:         os << "'raw' cap set:" << std::endl;
921:         for(typename traits::cap_container_type::set_type::iterator ci = _cap.set.begin(); ci != _cap.set.end(); ci++)
922:         {
923:           typename traits::cap_container_type::traits::rec_type cp = *ci;
924:           os << cp << std::endl;
925:         }
926:       }
927:     };
928:     // A parallel viewer
929:     PetscErrorCode view(const char* label = NULL, bool raw = false){
931:       ostringstream txt;
933:       if(this->_debug) {
934:         std::cout << "viewing a Sifter, comm = " << this->comm() << ", PETSC_COMM_SELF = " << PETSC_COMM_SELF << ", commRank = " << this->commRank() << std::endl;
935:       }
936:       if(label != NULL) {
937:         PetscPrintf(this->comm(), "viewing Sifter: '%s'\n", label);
938:       } else {
939:         PetscPrintf(this->comm(), "viewing a Sifter: \n");
940:       }
941:       if(!raw) {
942:         ostringstream txt;
943:         if(this->commRank() == 0) {
944:           txt << "cap --> base:\n";
945:         }
946:         typename traits::capSequence cap   = this->cap();
947:         typename traits::baseSequence base = this->base();
948:         if(cap.empty()) {
949:           txt << "[" << this->commRank() << "]: empty" << std::endl;
950:         }
951:         for(typename traits::capSequence::iterator capi = cap.begin(); capi != cap.end(); capi++) {
952:           const Obj<typename traits::supportSequence>& supp = this->support(*capi);
953:           for(typename traits::supportSequence::iterator suppi = supp->begin(); suppi != supp->end(); suppi++) {
954:             txt << "[" << this->commRank() << "]: " << *capi << "--(" << suppi.color() << ")-->" << *suppi << std::endl;
955:           }
956:         }
957:         //
958:         PetscSynchronizedPrintf(this->comm(), txt.str().c_str()); CHKERROR(ierr, "Error in PetscSynchronizedFlush");
959:         PetscSynchronizedFlush(this->comm());  CHKERROR(ierr, "Error in PetscSynchronizedFlush");
960: #if 0
961:         //
962:         ostringstream txt1;
963:         if(this->commRank() == 0) {
964:           //txt1 << "cap <point,degree>:\n";
965:           txt1 << "cap:\n";
966:         }
967:         txt1 << "[" << this->commRank() << "]:  [";
968:         for(typename traits::capSequence::iterator capi = cap.begin(); capi != cap.end(); capi++) {
969:           //txt1 << " <" << *capi << "," << capi.degree() << ">";
970:           txt1 << "  " << *capi;
971:         }
972:         txt1 << " ]" << std::endl;
973:         //
974:         PetscSynchronizedPrintf(this->comm(), txt1.str().c_str()); CHKERROR(ierr, "Error in PetscSynchronizedFlush");
975:         PetscSynchronizedFlush(this->comm());  CHKERROR(ierr, "Error in PetscSynchronizedFlush");
976:         //
977:         ostringstream txt2;
978:         if(this->commRank() == 0) {
979:           //txt2 << "base <point,degree>:\n";
980:           txt2 << "base:\n";
981:         }
982:         txt2 << "[" << this->commRank() << "]:  [";
983:         for(typename traits::baseSequence::iterator basei = base.begin(); basei != base.end(); basei++) {
984:           txt2 << "  " << *basei;
985:           //txt2 << " <" << *basei << "," << basei.degree() << ">";
986:         }
987:         txt2 << " ]" << std::endl;
988:         //
989:         PetscSynchronizedPrintf(this->comm(), txt2.str().c_str()); CHKERROR(ierr, "Error in PetscSynchronizedFlush");
990:         PetscSynchronizedFlush(this->comm());  CHKERROR(ierr, "Error in PetscSynchronizedFlush");
991: #endif
992:       }
993:       else { // if(raw)
994:         ostringstream txt;
995:         if(this->commRank() == 0) {
996:           txt << "'raw' arrow set:" << std::endl;
997:         }
998:         for(typename traits::arrow_container_type::set_type::iterator ai = _arrows.set.begin(); ai != _arrows.set.end(); ai++)
999:         {
1000:           typename traits::arrow_type arr = *ai;
1001:           txt << "[" << this->commRank() << "]: " << arr << std::endl;
1002:         }
1003:         PetscSynchronizedPrintf(this->comm(), txt.str().c_str()); CHKERROR(ierr, "Error in PetscSynchronizedFlush");
1004:         PetscSynchronizedFlush(this->comm());  CHKERROR(ierr, "Error in PetscSynchronizedFlush");
1005:         //
1006:         ostringstream txt1;
1007:         if(this->commRank() == 0) {
1008:           txt1 << "'raw' base set:" << std::endl;
1009:         }
1010:         for(typename traits::base_container_type::set_type::iterator bi = _base.set.begin(); bi != _base.set.end(); bi++)
1011:         {
1012:           typename traits::base_container_type::traits::rec_type bp = *bi;
1013:           txt1 << "[" << this->commRank() << "]: " << bp << std::endl;
1014:         }
1015:         PetscSynchronizedPrintf(this->comm(), txt1.str().c_str()); CHKERROR(ierr, "Error in PetscSynchronizedFlush");
1016:         PetscSynchronizedFlush(this->comm());  CHKERROR(ierr, "Error in PetscSynchronizedFlush");
1017:         //
1018:         ostringstream txt2;
1019:         if(this->commRank() == 0) {
1020:           txt2 << "'raw' cap set:" << std::endl;
1021:         }
1022:         for(typename traits::cap_container_type::set_type::iterator ci = _cap.set.begin(); ci != _cap.set.end(); ci++)
1023:         {
1024:           typename traits::cap_container_type::traits::rec_type cp = *ci;
1025:           txt2 << "[" << this->commRank() << "]: " << cp << std::endl;
1026:         }
1027:         PetscSynchronizedPrintf(this->comm(), txt2.str().c_str()); CHKERROR(ierr, "Error in PetscSynchronizedFlush");
1028:         PetscSynchronizedFlush(this->comm());  CHKERROR(ierr, "Error in PetscSynchronizedFlush");
1029:       }// if(raw)
1030: 
1031:       return(0);
1032:     };
1033:   public:
1034:     //
1035:     // Lattice queries
1036:     //
1037:     template<class targetInputSequence>
1038:     Obj<typename traits::coneSequence> meet(const Obj<targetInputSequence>& targets);
1039:     // unimplemented
1040:     template<class targetInputSequence>
1041:     Obj<typename traits::coneSequence> meet(const Obj<targetInputSequence>& targets, const typename traits::color_type& color);
1042:     // unimplemented
1043:     template<class sourceInputSequence>
1044:     Obj<typename traits::coneSequence> join(const Obj<sourceInputSequence>& sources);
1045:     // unimplemented
1046:     template<class sourceInputSequence>
1047:     Obj<typename traits::coneSequence> join(const Obj<sourceInputSequence>& sources, const typename traits::color_type& color);
1048:   public:
1049:     //
1050:     // Structural manipulation
1051:     //
1052:     void clear() {
1053:       this->_arrows.set.clear(); this->_base.set.clear(); this->_cap.set.clear();
1054:     };
1055:     void addBasePoint(const typename traits::target_type t) {
1056:       /* // Increase degree by 0, which won't affect an existing point and will insert a new point, if necessery
1057:          this->_base.adjustDegree(t,0); */
1058:       this->_base.set.insert(typename traits::targetRec_type(t,0));
1059:     };
1060:     void addBasePoint(const typename traits::targetRec_type b) {
1061:       this->_base.set.insert(b);
1062:     };
1063:     void removeBasePoint(const typename traits::target_type t) {
1064:       if (this->_debug) {std::cout << " Removing " << t << " from the base" << std::endl;}
1065:       // Clear the cone and remove the point from _base
1066:       this->clearCone(t);
1067:       this->_base.removePoint(t);
1068:     };
1069:     void addCapPoint(const typename traits::source_type s) {
1070:       /* // Increase degree by 0, which won't affect an existing point and will insert a new point, if necessery
1071:          this->_cap.adjustDegree(s,0); */
1072:       this->_cap.set.insert(typename traits::sourceRec_type(s,0));
1073:     };
1074:     void addCapPoint(const typename traits::sourceRec_type c) {
1075:       this->_cap.set.insert(c);
1076:     };
1077:     void removeCapPoint(const typename traits::source_type s) {
1078:       if (this->_debug) {std::cout << " Removing " << s << " from the cap" << std::endl;}
1079:       // Clear the support and remove the point from _cap
1080:       this->clearSupport(s);
1081:       this->_cap.removePoint(s);
1082:     };
1083:     virtual void addArrow(const typename traits::source_type& p, const typename traits::target_type& q) {
1084:       this->addArrow(p, q, typename traits::color_type());
1085:     };
1086:     virtual void addArrow(const typename traits::source_type& p, const typename traits::target_type& q, const typename traits::color_type& color) {
1087:       this->addArrow(typename traits::arrow_type(p, q, color));
1088:       //std::cout << "Added " << arrow_type(p, q, color);
1089:     };
1090:     virtual bool checkArrow(const typename traits::arrow_type& a) {
1091:       if (this->_cap.set.find(a.source) == this->_cap.set.end()) return false;
1092:       if (this->_base.set.find(a.target) == this->_base.set.end()) return false;
1093:       return true;
1094:     };
1095:     virtual void addArrow(const typename traits::arrow_type& a, bool noNewPoints = false) {
1096:       if (noNewPoints && !this->checkArrow(a)) return;
1097:       this->_arrows.set.insert(a);
1098:       this->addBasePoint(a.target);
1099:       this->addCapPoint(a.source);
1100:     };
1101:     virtual void removeArrow(const typename traits::source_type& p, const typename traits::target_type& q) {
1102:       this->removeArrow(typename traits::arrow_type(p, q, typename traits::color_type()));
1103:     };
1104:     virtual void removeArrow(const typename traits::arrow_type& a) {
1105:       // First, produce an arrow sequence for the given source, target combination.
1106:       typename traits::arrowSequence::traits::index_type& arrowIndex =
1107:         ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);
1108:       typename traits::arrowSequence::traits::index_type::iterator i,ii,j;
1109:       i = arrowIndex.lower_bound(::boost::make_tuple(a.source,a.target));
1110:       ii = arrowIndex.upper_bound(::boost::make_tuple(a.source, a.target));
1111:       if (this->_debug) {
1112:         std::cout << "removeArrow: attempting to remove arrow:" << a << std::endl;
1113:         std::cout << "removeArrow: candidate arrows are:" << std::endl;
1114:       }
1115:       for(j = i; j != ii; j++) {
1116:         if (this->_debug) {
1117:           std::cout << " " << *j;
1118:         }
1119:         // Find the arrow of right color and remove it
1120:         if(j->color == a.color) {
1121:           if (this->_debug) {
1122:             std::cout << std::endl << "removeArrow: found:" << *j << std::endl;
1123:           }
1124:           /* this->_base.adjustDegree(a.target, -1); this->_cap.adjustDegree(a.source,-1); */
1125:           arrowIndex.erase(j);
1126:           break;
1127:         }
1128:       }
1129:     };

1131:     void addCone(const typename traits::source_type& source, const typename traits::target_type& target){
1132:       this->addArrow(source, target);
1133:     };
1134:     template<class sourceInputSequence>
1135:     void addCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target) {
1136:       this->addCone(sources, target, typename traits::color_type());
1137:     };
1138:     void addCone(const typename traits::source_type& source, const typename traits::target_type& target, const typename traits::color_type& color) {
1139:       this->addArrow(source, target, color);
1140:     };
1141:     template<class sourceInputSequence>
1142:     void
1143:     addCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target, const typename traits::color_type& color){
1144:       if (this->_debug > 1) {std::cout << "Adding a cone " << std::endl;}
1145:       for(typename sourceInputSequence::iterator iter = sources->begin(); iter != sources->end(); ++iter) {
1146:         if (this->_debug > 1) {std::cout << "Adding arrow from " << *iter << " to " << target << "(" << color << ")" << std::endl;}
1147:         this->addArrow(*iter, target, color);
1148:       }
1149:     };
1150:     void clearCone(const typename traits::target_type& t) {
1151:       clearCone(t, typename traits::color_type(), false);
1152:     };

1154:     void clearCone(const typename traits::target_type& t, const typename traits::color_type&  color, bool useColor = true) {
1155:       // Use the cone sequence types to clear the cone
1156:       typename traits::coneSequence::traits::index_type& coneIndex =
1157:         ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set);
1158:       typename traits::coneSequence::traits::index_type::iterator i, ii, j;
1159:       if (this->_debug > 20) {
1160:         std::cout << "clearCone: removing cone over " << t;
1161:         if(useColor) {
1162:           std::cout << " with color" << color << std::endl;
1163:           const Obj<typename traits::coneSequence>& cone = this->cone(t,color);
1164:           std::cout << "[";
1165:           for(typename traits::coneSequence::iterator ci = cone->begin(); ci != cone->end(); ci++) {
1166:             std::cout << "  " << ci.arrow();
1167:           }
1168:           std::cout << "]" << std::endl;
1169:         }
1170:         else {
1171:           std::cout << std::endl;
1172:           const Obj<typename traits::coneSequence>& cone = this->cone(t);
1173:           std::cout << "[";
1174:           for(typename traits::coneSequence::iterator ci = cone->begin(); ci != cone->end(); ci++) {
1175:             std::cout << "  " << ci.arrow();
1176:           }
1177:           std::cout << "]" << std::endl;
1178:         }
1179:       }
1180:       if (useColor) {
1181:         i = coneIndex.lower_bound(::boost::make_tuple(t,color));
1182:         ii = coneIndex.upper_bound(::boost::make_tuple(t,color));
1183:       } else {
1184:         i = coneIndex.lower_bound(::boost::make_tuple(t));
1185:         ii = coneIndex.upper_bound(::boost::make_tuple(t));
1186:       }
1187:       for(j = i; j != ii; j++){
1188:         // Adjust the degrees before removing the arrow; use a different iterator, since we'll need i,ii to do the arrow erasing.
1189:         if (this->_debug) {
1190:           std::cout << "clearCone: adjusting degrees for endpoints of arrow: " << *j << std::endl;
1191:         }
1192:         /* this->_cap.adjustDegree(j->source, -1);
1193:            this->_base.adjustDegree(j->target, -1); */
1194:       }
1195:       coneIndex.erase(i,ii);
1196:     };// clearCone()

1198:     template<class InputSequence>
1199:     void
1200:     restrictBase(const Obj<InputSequence>& points) {
1201:       typename traits::baseSequence base = this->base();
1202:       typename std::set<typename traits::target_type> remove;
1203: 
1204:       for(typename traits::baseSequence::iterator bi = base.begin(); bi != base.end(); bi++) {
1205:         // Check whether *bi is in points, if it is NOT, remove it
1206:         //           if (!points->contains(*bi)) {
1207:         if (points->find(*bi) == points->end()) {
1208:           //             this->removeBasePoint(*bi);
1209:           remove.insert(*bi);
1210:         }
1211:       }
1212:       //FIX
1213:       for(typename std::set<typename traits::target_type>::iterator r_iter = remove.begin(); r_iter != remove.end(); ++r_iter) {
1214:         this->removeBasePoint(*r_iter);
1215:       }
1216:     };

1218:     template<class InputSequence>
1219:     void
1220:     excludeBase(const Obj<InputSequence>& points) {
1221:       for(typename InputSequence::iterator pi = points->begin(); pi != points->end(); pi++) {
1222:         this->removeBasePoint(*pi);
1223:       }
1224:     };

1226:     template<class InputSequence>
1227:     void
1228:     restrictCap(const Obj<InputSequence>& points) {
1229:       typename traits::capSequence cap = this->cap();
1230:       for(typename traits::capSequence::iterator ci = cap.begin(); ci != cap.end(); ci++) {
1231:         // Check whether *ci is in points, if it is NOT, remove it
1232:         if(points->find(*ci) == points->end()) {
1233:           this->removeCapPoint(*ci);
1234:         }
1235:       }
1236:     };

1238:     template<class InputSequence>
1239:     void
1240:     excludeCap(const Obj<InputSequence>& points) {
1241:       for(typename InputSequence::iterator pi = points->begin(); pi != points->end(); pi++) {
1242:         this->removeCapPoint(*pi);
1243:       }
1244:     };

1246:     void clearSupport(const typename traits::source_type& s) {
1247:       clearSupport(s, typename traits::color_type(), false);
1248:     };
1249:     void clearSupport(const typename traits::source_type& s, const typename traits::color_type&  color, bool useColor = true) {
1250:       // Use the cone sequence types to clear the cone
1251:       typename
1252:         traits::supportSequence::traits::index_type& suppIndex = ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set);
1253:       typename traits::supportSequence::traits::index_type::iterator i, ii, j;
1254:       if (useColor) {
1255:         i = suppIndex.lower_bound(::boost::make_tuple(s,color));
1256:         ii = suppIndex.upper_bound(::boost::make_tuple(s,color));
1257:       } else {
1258:         i = suppIndex.lower_bound(::boost::make_tuple(s));
1259:         ii = suppIndex.upper_bound(::boost::make_tuple(s));
1260:       }
1261:       for(j = i; j != ii; j++){
1262:         // Adjust the degrees before removing the arrow
1263:         /* this->_cap.adjustDegree(j->source, -1);
1264:            this->_base.adjustDegree(j->target, -1); */
1265:       }
1266:       suppIndex.erase(i,ii);
1267:     };
1268:     void setCone(const typename traits::source_type& source, const typename traits::target_type& target){
1269:       this->clearCone(target, typename traits::color_type(), false); this->addCone(source, target);
1270:     };
1271:     template<class sourceInputSequence>
1272:     void setCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target) {
1273:       this->clearCone(target, typename traits::color_type(), false); this->addCone(sources, target, typename traits::color_type());
1274:     };
1275:     void setCone(const typename traits::source_type& source, const typename traits::target_type& target, const typename traits::color_type& color) {
1276:       this->clearCone(target, color, true); this->addCone(source, target, color);
1277:     };
1278:     template<class sourceInputSequence>
1279:     void setCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target, const typename traits::color_type& color){
1280:       this->clearCone(target, color, true); this->addCone(sources, target, color);
1281:     };
1282:     template<class targetInputSequence>
1283:     void addSupport(const typename traits::source_type& source, const Obj<targetInputSequence >& targets) {
1284:       this->addSupport(source, targets, typename traits::color_type());
1285:     };
1286:     template<class targetInputSequence>
1287:     void addSupport(const typename traits::source_type& source, const Obj<targetInputSequence>& targets, const typename traits::color_type& color) {
1288:       const typename targetInputSequence::iterator end = targets->end();

1290:       for(typename targetInputSequence::iterator iter = targets->begin(); iter != end; ++iter) {
1291:         this->addArrow(source, *iter, color);
1292:       }
1293:     };
1294:     template<typename Sifter_>
1295:     void add(const Obj<Sifter_>& cbg, bool noNewPoints = false) {
1296:       typename ::boost::multi_index::index<typename Sifter_::traits::arrow_container_type::set_type, typename Sifter_::traits::arrowInd>::type& aInd = ::boost::multi_index::get<typename Sifter_::traits::arrowInd>(cbg->_arrows.set);
1297: 
1298:       for(typename ::boost::multi_index::index<typename Sifter_::traits::arrow_container_type::set_type, typename Sifter_::traits::arrowInd>::type::iterator a_iter = aInd.begin(); a_iter != aInd.end(); ++a_iter) {
1299:         this->addArrow(*a_iter, noNewPoints);
1300:       }
1301:       if (!noNewPoints) {
1302:         typename ::boost::multi_index::index<typename Sifter_::traits::base_container_type::set_type, typename Sifter_::traits::baseInd>::type& bInd = ::boost::multi_index::get<typename Sifter_::traits::baseInd>(this->_base.set);
1303: 
1304:         for(typename ::boost::multi_index::index<typename Sifter_::traits::base_container_type::set_type, typename Sifter_::traits::baseInd>::type::iterator b_iter = bInd.begin(); b_iter != bInd.end(); ++b_iter) {
1305:           this->addBasePoint(*b_iter);
1306:         }
1307:         typename ::boost::multi_index::index<typename Sifter_::traits::cap_container_type::set_type, typename Sifter_::traits::capInd>::type& cInd = ::boost::multi_index::get<typename Sifter_::traits::capInd>(this->_cap.set);
1308: 
1309:         for(typename ::boost::multi_index::index<typename Sifter_::traits::cap_container_type::set_type, typename Sifter_::traits::capInd>::type::iterator c_iter = cInd.begin(); c_iter != cInd.end(); ++c_iter) {
1310:           this->addCapPoint(*c_iter);
1311:         }
1312:       }
1313:     };
1314:   }; // class ASifter

1316:   // A UniSifter aka Sifter
1317:   template <typename Source_, typename Target_, typename Color_,
1318:             typename SupportCompare_ = ::boost::multi_index::composite_key_compare<std::less<Source_>, std::less<Color_>, std::less<Target_> >,
1319:             typename SourceCtnr_ = SifterDef:: RecContainer<Source_, SifterDef::Rec<Source_> >, typename TargetCtnr_= SifterDef::RecContainer<Target_, SifterDef::Rec<Target_> > >
1320:   class Sifter : public ASifter<Source_, Target_, Color_, SifterDef::uniColor, SupportCompare_, SourceCtnr_, TargetCtnr_> {
1321:   public:
1322:       typedef typename ASifter<Source_, Target_, Color_, SifterDef::uniColor, SupportCompare_, SourceCtnr_, TargetCtnr_>::traits       traits;
1323:     template <typename OtherSource_, typename OtherTarget_, typename OtherColor_,
1324:               typename OtherSupportCompare_  = ::boost::multi_index::composite_key_compare<std::less<OtherSource_>, std::less<OtherColor_>, std::less<OtherTarget_> >,
1325:               typename OtherSourceCtnr_ = SifterDef::RecContainer<OtherSource_, SifterDef::Rec<OtherSource_> >,
1326:               typename OtherTargetCtnr_ = SifterDef::RecContainer<OtherTarget_, SifterDef::Rec<OtherTarget_> >      >
1327:     struct rebind {
1328:       typedef Sifter<OtherSource_, OtherTarget_, OtherColor_, OtherSupportCompare_, OtherSourceCtnr_, OtherTargetCtnr_> type;
1329:     };
1330:     // Re-export some typedefs expected by CoSifter
1331:     typedef typename traits::source_type                                            source_type;
1332:     typedef typename traits::target_type                                            target_type;
1333:     typedef typename traits::arrow_type                                             Arrow_;
1334:     typedef typename traits::coneSequence                                           coneSequence;
1335:     typedef typename traits::supportSequence                                        supportSequence;
1336:     typedef typename traits::baseSequence                                           baseSequence;
1337:     typedef typename traits::capSequence                                            capSequence;
1338:     // Basic interface
1339:     Sifter(MPI_Comm comm = PETSC_COMM_SELF, const int& debug = 0) :
1340:       ASifter<Source_, Target_, Color_, SifterDef::uniColor, SupportCompare_, SourceCtnr_, TargetCtnr_>(comm, debug) {};

1342:     const typename traits::color_type&
1343:     getColor(const typename traits::source_type& s, const typename traits::target_type& t, bool fail = true) {
1344:       typedef typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type,typename traits::arrowInd>::type index_type;

1346:       const index_type& _index = ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);
1347: #if 0
1348:       ::boost::tuple<typename traits::source_type, typename traits::target_type> key = ::boost::make_tuple(s, t);
1349:       typename index_type::iterator begin = _index.lower_bound(key);
1350:       if(begin != _index.upper_bound(key)) {
1351:         return begin->color;
1352:       }
1353: #else
1354:       const typename index_type::iterator begin = _index.find(::boost::make_tuple(s, t));
1355:       if (begin != _index.end()) {
1356:         return begin->color;
1357:       }
1358: #endif
1359: //       typename traits::arrowSequence arr(::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set), s, t);
1360: //       if(arr.begin() != arr.end()) {
1361: //         return arr.begin().color();
1362: //       }
1363:       if (fail) {
1364:         ostringstream o;
1365:         o << "Arrow " << s << " --> " << t << " not present";
1366:         throw ALE::Exception(o.str().c_str());
1367:       } else {
1368:         static typename traits::color_type c;
1369:         return c;
1370:       }
1371:     };

1373:     template<typename ColorChanger>
1374:     void modifyColor(const typename traits::source_type& s, const typename traits::target_type& t, const ColorChanger& changeColor) {
1375:       typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type, typename traits::arrowInd>::type& index =
1376:         ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);
1377:       typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type, typename traits::arrowInd>::type::iterator i =
1378:         index.find(::boost::make_tuple(s,t));
1379:       if (i != index.end()) {
1380:         index.modify(i, changeColor);
1381:       } else {
1382:         typename traits::arrow_type a(s, t, typename traits::color_type());
1383:         changeColor(a);
1384:         this->addArrow(a);
1385:       }
1386:     };

1388:     struct ColorSetter {
1389:       ColorSetter(const typename traits::color_type& color) : _color(color) {};
1390:       void operator()(typename traits::arrow_type& p) const {
1391:         p.color = _color;
1392:       }
1393:     private:
1394:       const typename traits::color_type& _color;
1395:     };

1397:     void setColor(const typename traits::source_type& s, const typename traits::target_type& t, const typename traits::color_type& color) {
1398:       ColorSetter colorSetter(color);
1399:       typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type, typename traits::arrowInd>::type& index =
1400:         ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);
1401:       typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type, typename traits::arrowInd>::type::iterator i =
1402:         index.find(::boost::make_tuple(s,t));
1403:       if (i != index.end()) {
1404:         index.modify(i, colorSetter);
1405:       } else {
1406:         typename traits::arrow_type a(s, t, color);
1407:         this->addArrow(a);
1408:       }
1409:     };
1410:   };// class Sifter

1412: } // namespace ALE

1414: #endif