Actual source code: Sifter.hh
petsc-3.4.2 2013-07-02

      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 <sieve/ALE.hh>
     16: #endif

     18: extern PetscErrorCode PetscObjectDestroy_PetscObject(PetscObject*);

     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;

     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

     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:       };

    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:

    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:        };

    173:        PointSequence(const PointSequence& seq)            : _index(seq._index) {};
    174:        PointSequence(typename traits::index_type& index) : _index(index)     {};
    175:        virtual ~PointSequence(){};

    177:        virtual bool empty(){return this->_index.empty();};

    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->set.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

    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;};

    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:         };

    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:         };

    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:         };

    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

    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 {};

    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>

    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() {};

    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:         };

    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

    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() {};

    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

    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 PetscClassId sifterType = -1;
    656:       //const char        *id_name = ALE::getClassName<T>();
    657:       const char        *id_name = "Sifter";
    658:       PetscErrorCode     ierr;

    660:       if (sifterType < 0) {
    661:         PetscClassIdRegister(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
    686:       PetscObjectDestroy(&this->_petscObj);CHKERROR(ierr, "Failed in PetscObjectDestroy");
    687: #endif
    688:     };
    689:     //
    690:     // Query methods
    691:     //
    692:     int         debug()    const {return this->_debug;};
    693:     void        setDebug(const int debug) {this->_debug = debug;};
    694:     MPI_Comm    comm()     const {return this->_comm;};
    695:     int         commSize() const {return this->_commSize;};
    696:     int         commRank() const {return this->_commRank;}
    697: #ifdef USE_PETSC_OBJ
    698:     PetscObject petscObj() const {return this->_petscObj;};
    699: #endif

    701:     // Added to allow optimized versions
    702:     void assemble() {};
    703:     void assemblePoints() {};
    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:     typename traits::capSequence::iterator capBegin() {
    712:       return this->cap()->begin();
    713:     };
    714:     typename traits::capSequence::iterator capEnd() {
    715:       return this->cap()->end();
    716:     };
    717:     typename traits::baseSequence::iterator baseBegin() {
    718:       return this->base()->begin();
    719:     };
    720:     typename traits::baseSequence::iterator baseEnd() {
    721:       return this->base()->end();
    722:     };
    723:     int getBaseSize() {return this->base()->size();};
    724:     void setBaseSize(int size) {};
    725:     int getCapSize() {return this->cap()->size();};
    726:     void setCapSize(int size) {};
    727:     bool capContains(const typename traits::source_type& p) {
    728:       typename traits::capSequence cap(::boost::multi_index::get<typename traits::capInd>(this->_cap.set));

    730:       //for(typename traits::capSequence::iterator c_iter = cap.begin(); c_iter != cap.end(); ++c_iter) {
    731:       //}
    732:       return cap.contains(p);
    733:     };
    734:     bool baseContains(const typename traits::target_type& p) {
    735:       typename traits::baseSequence base(::boost::multi_index::get<typename traits::baseInd>(this->_base.set));

    737:       //for(typename traits::capSequence::iterator c_iter = cap.begin(); c_iter != cap.end(); ++c_iter) {
    738:       //}
    739:       return base.contains(p);
    740:     };
    741:     // FIX: should probably have cone and const_cone etc, since arrows can be modified through an iterator (modifyColor).
    742:     Obj<typename traits::arrowSequence>
    743:     arrows(const typename traits::source_type& s, const typename traits::target_type& t) {
    744:       return typename traits::arrowSequence(::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set), s, t);
    745:     };
    746:     Obj<typename traits::arrowSequence>
    747:     arrows(const typename traits::source_type& s) {
    748:       return typename traits::arrowSequence(::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set), s);
    749:     };
    750: #ifdef SLOW
    751:     Obj<typename traits::coneSequence>
    752:     cone(const typename traits::target_type& p) {
    753:       return typename traits::coneSequence(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p);
    754:     };
    755: #else
    756:     const Obj<typename traits::coneSequence>&
    757:     cone(const typename traits::target_type& p) {
    758:       this->_coneSeq->setKey(p);
    759:       this->_coneSeq->setUseSubkey(false);
    760:       return this->_coneSeq;
    761:     };
    762:     const typename traits::coneSequence::iterator
    763:     coneBegin(const typename traits::target_type& p) {
    764:       return this->cone(p)->begin();
    765:     };
    766:     const typename traits::coneSequence::iterator
    767:     coneEnd(const typename traits::target_type& p) {
    768:       return this->cone(p)->end();
    769:     };
    770:     void setConeSize(const typename traits::target_type& p, int size) {};
    771: #endif
    772:     template<class InputSequence>
    773:     Obj<typename traits::coneSet>
    774:     cone(const Obj<InputSequence>& points) {
    775:       return this->cone(points, typename traits::color_type(), false);
    776:     }
    777: #ifdef SLOW
    778:     Obj<typename traits::coneSequence>
    779:     cone(const typename traits::target_type& p, const typename traits::color_type& color) {
    780:       return typename traits::coneSequence(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p, color);
    781:     }
    782: #else
    783:     const Obj<typename traits::coneSequence>&
    784:     cone(const typename traits::target_type& p, const typename traits::color_type& color) {
    785:       this->_coneSeq->setKey(p);
    786:       this->_coneSeq->setSubkey(color);
    787:       this->_coneSeq->setUseSubkey(true);
    788:       return this->_coneSeq;
    789:     };
    790: #endif
    791:     template<class InputSequence>
    792:     Obj<typename traits::coneSet>
    793:     cone(const Obj<InputSequence>& points, const typename traits::color_type& color, bool useColor = true) {
    794:       Obj<typename traits::coneSet> cone = typename traits::coneSet();
    795:       for(typename InputSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
    796:         Obj<typename traits::coneSequence> pCone;
    797:         if (useColor) {
    798:           pCone = this->cone(*p_itor, color);
    799:         } else {
    800:           pCone = this->cone(*p_itor);
    801:         }
    802:         cone->insert(pCone->begin(), pCone->end());
    803:       }
    804:       return cone;
    805:     }
    806:     template<typename PointCheck>
    807:     bool coneContains(const typename traits::target_type& p, const PointCheck& checker) {
    808:       typename traits::coneSequence cone(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p);

    810:       for(typename traits::coneSequence::iterator c_iter = cone.begin(); c_iter != cone.end(); ++c_iter) {
    811:         if (checker(*c_iter, p)) return true;
    812:       }
    813:       return false;
    814:     }
    815:     template<typename PointProcess>
    816:     void coneApply(const typename traits::target_type& p, PointProcess& processor) {
    817:       typename traits::coneSequence cone(*this, ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set), p);

    819:       for(typename traits::coneSequence::iterator c_iter = cone.begin(); c_iter != cone.end(); ++c_iter) {
    820:         processor(*c_iter, p);
    821:       }
    822:     }
    823:     int getConeSize(const typename traits::target_type& p) {
    824:       return this->cone(p)->size();
    825:     }
    826: #ifdef SLOW
    827:     Obj<typename traits::supportSequence>
    828:     support(const typename traits::source_type& p) {
    829:       return typename traits::supportSequence(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p);
    830:     }
    831: #else
    832:     const Obj<typename traits::supportSequence>&
    833:     support(const typename traits::source_type& p) {
    834:       this->_supportSeq->setKey(p);
    835:       this->_supportSeq->setUseSubkey(false);
    836:       return this->_supportSeq;
    837:     };
    838:     const typename traits::supportSequence::iterator
    839:     supportBegin(const typename traits::source_type& p) {
    840:       return this->support(p)->begin();
    841:     };
    842:     const typename traits::supportSequence::iterator
    843:     supportEnd(const typename traits::source_type& p) {
    844:       return this->support(p)->end();
    845:     };
    846:     const typename traits::supportSequence::iterator
    847:     supportBegin(const typename traits::source_type& p, const typename traits::color_type& color) {
    848:       return this->support(p, color)->begin();
    849:     };
    850:     const typename traits::supportSequence::iterator
    851:     supportEnd(const typename traits::source_type& p, const typename traits::color_type& color) {
    852:       return this->support(p, color)->end();
    853:     };
    854:     void setSupportSize(const typename traits::source_type& p, int size) {};
    855: #endif
    856: #ifdef SLOW
    857:     Obj<typename traits::supportSequence>
    858:     support(const typename traits::source_type& p, const typename traits::color_type& color) {
    859:       return typename traits::supportSequence(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p, color);
    860:     };
    861: #else
    862:     const Obj<typename traits::supportSequence>&
    863:     support(const typename traits::source_type& p, const typename traits::color_type& color) {
    864:       this->_supportSeq->setKey(p);
    865:       this->_supportSeq->setSubkey(color);
    866:       this->_supportSeq->setUseSubkey(true);
    867:       return this->_supportSeq;
    868:     };
    869: #endif
    870:     template<class InputSequence>
    871:     Obj<typename traits::supportSet>
    872:     support(const Obj<InputSequence>& sources) {
    873:       return this->support(sources, typename traits::color_type(), false);
    874:     }
    875:     template<class InputSequence>
    876:     Obj<typename traits::supportSet>
    877:     support(const Obj<InputSequence>& points, const typename traits::color_type& color, bool useColor = true){
    878:       Obj<typename traits::supportSet> supp = typename traits::supportSet();
    879:       for(typename InputSequence::iterator p_itor = points->begin(); p_itor != points->end(); ++p_itor) {
    880:         Obj<typename traits::supportSequence> pSupport;
    881:         if (useColor) {
    882:           pSupport = this->support(*p_itor, color);
    883:         } else {
    884:           pSupport = this->support(*p_itor);
    885:         }
    886:         supp->insert(pSupport->begin(), pSupport->end());
    887:       }
    888:       return supp;
    889:     }
    890:     template<typename PointCheck>
    891:     bool supportContains(const typename traits::source_type& p, const PointCheck& checker) {
    892:       typename traits::supportSequence support(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p);

    894:       for(typename traits::supportSequence::iterator s_iter = support.begin(); s_iter != support.end(); ++s_iter) {
    895:         if (checker(*s_iter, p)) return true;
    896:       }
    897:       return false;
    898:     }
    899:     template<typename PointProcess>
    900:     void supportApply(const typename traits::source_type& p, PointProcess& processor) {
    901:       typename traits::supportSequence support(*this, ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set), p);

    903:       for(typename traits::supportSequence::iterator s_iter = support.begin(); s_iter != support.end(); ++s_iter) {
    904:         processor(*s_iter, p);
    905:       }
    906:     }
    907:     int getSupportSize(const typename traits::source_type& p) {
    908:       return this->support(p)->size();
    909:     }
    910:     int getSupportSize(const typename traits::source_type& p, const typename traits::color_type& color) {
    911:       return this->support(p, color)->size();
    912:     }

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

    918:       if(label != NULL) {
    919:         os << "["<<rank<<"]Viewing Sifter '" << label << "':" << std::endl;
    920:       }
    921:       else {
    922:         os << "["<<rank<<"]Viewing a Sifter:" << std::endl;
    923:       }
    924:       if(!rawData) {
    925:         os << "cap --> base:" << std::endl;
    926:         Obj<typename traits::capSequence> cap = this->cap();
    927:         for(typename traits::capSequence::iterator capi = cap->begin(); capi != cap->end(); capi++) {
    928:           const Obj<typename traits::supportSequence>& supp = this->support(*capi);

    930:           for(typename traits::supportSequence::iterator suppi = supp->begin(); suppi != supp->end(); suppi++) {
    931:             os << *capi << "--(" << suppi.color() << ")-->" << *suppi << std::endl;
    932:           }
    933:         }
    934:         os << "base <-- cap:" << std::endl;
    935:         Obj<typename traits::baseSequence> base = this->base();
    936:         for(typename traits::baseSequence::iterator basei = base->begin(); basei != base->end(); basei++) {
    937:           const Obj<typename traits::coneSequence>& cone = this->cone(*basei);

    939:           for(typename traits::coneSequence::iterator conei = cone->begin(); conei != cone->end(); conei++) {
    940:             os << *basei <<  "<--(" << conei.color() << ")--" << *conei << std::endl;
    941:           }
    942:         }
    943:         os << "cap --> outdegrees:" << std::endl;
    944:         for(typename traits::capSequence::iterator capi = cap->begin(); capi != cap->end(); capi++) {
    945:           os << *capi <<  "-->" << capi.degree() << std::endl;
    946:         }
    947:         os << "base <-- indegrees:" << std::endl;
    948:         for(typename traits::baseSequence::iterator basei = base->begin(); basei != base->end(); basei++) {
    949:           os << *basei <<  "<--" << basei.degree() << std::endl;
    950:         }
    951:       }
    952:       else {
    953:         os << "'raw' arrow set:" << std::endl;
    954:         for(typename traits::arrow_container_type::set_type::iterator ai = _arrows.set.begin(); ai != _arrows.set.end(); ai++)
    955:         {
    956:           typename traits::arrow_type arr = *ai;
    957:           os << arr << std::endl;
    958:         }
    959:         os << "'raw' base set:" << std::endl;
    960:         for(typename traits::base_container_type::set_type::iterator bi = _base.set.begin(); bi != _base.set.end(); bi++)
    961:         {
    962:           typename traits::base_container_type::traits::rec_type bp = *bi;
    963:           os << bp << std::endl;
    964:         }
    965:         os << "'raw' cap set:" << std::endl;
    966:         for(typename traits::cap_container_type::set_type::iterator ci = _cap.set.begin(); ci != _cap.set.end(); ci++)
    967:         {
    968:           typename traits::cap_container_type::traits::rec_type cp = *ci;
    969:           os << cp << std::endl;
    970:         }
    971:       }
    972:     }
    973:     // A parallel viewer
    976:     PetscErrorCode view(const char* label = NULL, bool raw = false){
    978:       ostringstream txt;
    980:       if(this->_debug) {
    981:         std::cout << "viewing a Sifter, comm = " << this->comm() << ", PETSC_COMM_SELF = " << PETSC_COMM_SELF << ", commRank = " << this->commRank() << std::endl;
    982:       }
    983:       if(label != NULL) {
    984:         PetscPrintf(this->comm(), "viewing Sifter: '%s'\n", label);
    985:       } else {
    986:         PetscPrintf(this->comm(), "viewing a Sifter: \n");
    987:       }
    988:       if(!raw) {
    989:         ostringstream txt;
    990:         if(this->commRank() == 0) {
    991:           txt << "cap --> base:\n";
    992:         }
    993:         typename traits::capSequence cap   = this->cap();
    994:         typename traits::baseSequence base = this->base();
    995:         if(cap.empty()) {
    996:           txt << "[" << this->commRank() << "]: empty" << std::endl;
    997:         }
    998:         for(typename traits::capSequence::iterator capi = cap.begin(); capi != cap.end(); capi++) {
    999:           const Obj<typename traits::supportSequence>& supp = this->support(*capi);
    1000:           for(typename traits::supportSequence::iterator suppi = supp->begin(); suppi != supp->end(); suppi++) {
    1001:             txt << "[" << this->commRank() << "]: " << *capi << "--(" << suppi.color() << ")-->" << *suppi << std::endl;
    1002:           }
    1003:         }
    1004:         //
    1005:         PetscSynchronizedPrintf(this->comm(), txt.str().c_str());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1006:         PetscSynchronizedFlush(this->comm());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1007: #if 0
    1008:         //
    1009:         ostringstream txt1;
    1010:         if(this->commRank() == 0) {
    1011:           //txt1 << "cap <point,degree>:\n";
    1012:           txt1 << "cap:\n";
    1013:         }
    1014:         txt1 << "[" << this->commRank() << "]:  [";
    1015:         for(typename traits::capSequence::iterator capi = cap.begin(); capi != cap.end(); capi++) {
    1016:           //txt1 << " <" << *capi << "," << capi.degree() << ">";
    1017:           txt1 << "  " << *capi;
    1018:         }
    1019:         txt1 << " ]" << std::endl;
    1020:         //
    1021:         PetscSynchronizedPrintf(this->comm(), txt1.str().c_str());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1022:         PetscSynchronizedFlush(this->comm());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1023:         //
    1024:         ostringstream txt2;
    1025:         if(this->commRank() == 0) {
    1026:           //txt2 << "base <point,degree>:\n";
    1027:           txt2 << "base:\n";
    1028:         }
    1029:         txt2 << "[" << this->commRank() << "]:  [";
    1030:         for(typename traits::baseSequence::iterator basei = base.begin(); basei != base.end(); basei++) {
    1031:           txt2 << "  " << *basei;
    1032:           //txt2 << " <" << *basei << "," << basei.degree() << ">";
    1033:         }
    1034:         txt2 << " ]" << std::endl;
    1035:         //
    1036:         PetscSynchronizedPrintf(this->comm(), txt2.str().c_str());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1037:         PetscSynchronizedFlush(this->comm());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1038: #endif
    1039:       }
    1040:       else { // if(raw)
    1041:         ostringstream txt;
    1042:         if(this->commRank() == 0) {
    1043:           txt << "'raw' arrow set:" << std::endl;
    1044:         }
    1045:         for(typename traits::arrow_container_type::set_type::iterator ai = _arrows.set.begin(); ai != _arrows.set.end(); ai++)
    1046:         {
    1047:           typename traits::arrow_type arr = *ai;
    1048:           txt << "[" << this->commRank() << "]: " << arr << std::endl;
    1049:         }
    1050:         PetscSynchronizedPrintf(this->comm(), txt.str().c_str());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1051:         PetscSynchronizedFlush(this->comm());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1052:         //
    1053:         ostringstream txt1;
    1054:         if(this->commRank() == 0) {
    1055:           txt1 << "'raw' base set:" << std::endl;
    1056:         }
    1057:         for(typename traits::base_container_type::set_type::iterator bi = _base.set.begin(); bi != _base.set.end(); bi++)
    1058:         {
    1059:           typename traits::base_container_type::traits::rec_type bp = *bi;
    1060:           txt1 << "[" << this->commRank() << "]: " << bp << std::endl;
    1061:         }
    1062:         PetscSynchronizedPrintf(this->comm(), txt1.str().c_str());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1063:         PetscSynchronizedFlush(this->comm());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1064:         //
    1065:         ostringstream txt2;
    1066:         if(this->commRank() == 0) {
    1067:           txt2 << "'raw' cap set:" << std::endl;
    1068:         }
    1069:         for(typename traits::cap_container_type::set_type::iterator ci = _cap.set.begin(); ci != _cap.set.end(); ci++)
    1070:         {
    1071:           typename traits::cap_container_type::traits::rec_type cp = *ci;
    1072:           txt2 << "[" << this->commRank() << "]: " << cp << std::endl;
    1073:         }
    1074:         PetscSynchronizedPrintf(this->comm(), txt2.str().c_str());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1075:         PetscSynchronizedFlush(this->comm());CHKERROR(ierr, "Error in PetscSynchronizedFlush");
    1076:       }// if(raw)

    1078:       return(0);
    1079:     };
    1080:   public:
    1081:     void copy(ASifter *newSifter) {
    1082:       const typename traits::baseSequence::iterator sBegin = this->baseBegin();
    1083:       const typename traits::baseSequence::iterator sEnd   = this->baseEnd();

    1085:       for(typename traits::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
    1086:         const typename traits::coneSequence::iterator pBegin = this->coneBegin(*r_iter);
    1087:         const typename traits::coneSequence::iterator pEnd   = this->coneEnd(*r_iter);

    1089:         for(typename traits::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
    1090:           const Obj<typename traits::supportSequence>&              support  = this->support(*p_iter);
    1091:           const typename traits::supportSequence::iterator supBegin = support->begin();
    1092:           const typename traits::supportSequence::iterator supEnd   = support->end();

    1094:           for(typename traits::supportSequence::iterator s_iter = supBegin; s_iter != supEnd; ++s_iter) {
    1095:             newSifter->addArrow(*p_iter, *s_iter, s_iter.color());
    1096:           }
    1097:         }
    1098:       }
    1099:     };
    1100:     //
    1101:     // Lattice queries
    1102:     //
    1103:     template<class targetInputSequence>
    1104:     Obj<typename traits::coneSequence> meet(const Obj<targetInputSequence>& targets);
    1105:     // unimplemented
    1106:     template<class targetInputSequence>
    1107:     Obj<typename traits::coneSequence> meet(const Obj<targetInputSequence>& targets, const typename traits::color_type& color);
    1108:     // unimplemented
    1109:     template<class sourceInputSequence>
    1110:     Obj<typename traits::coneSequence> join(const Obj<sourceInputSequence>& sources);
    1111:     // unimplemented
    1112:     template<class sourceInputSequence>
    1113:     Obj<typename traits::coneSequence> join(const Obj<sourceInputSequence>& sources, const typename traits::color_type& color);
    1114:   public:
    1115:     //
    1116:     // Structural manipulation
    1117:     //
    1118:     void clear() {
    1119:       this->_arrows.set.clear(); this->_base.set.clear(); this->_cap.set.clear();
    1120:     };
    1121:     void addBasePoint(const typename traits::target_type t) {
    1122:       /* // Increase degree by 0, which won't affect an existing point and will insert a new point, if necessery
    1123:          this->_base.adjustDegree(t,0); */
    1124:       this->_base.set.insert(typename traits::targetRec_type(t,0));
    1125:     };
    1126:     void addBasePoint(const typename traits::targetRec_type b) {
    1127:       this->_base.set.insert(b);
    1128:     };
    1129:     void removeBasePoint(const typename traits::target_type t) {
    1130:       if (this->_debug) {std::cout << " Removing " << t << " from the base" << std::endl;}
    1131:       // Clear the cone and remove the point from _base
    1132:       this->clearCone(t);
    1133:       this->_base.removePoint(t);
    1134:     };
    1135:     void addCapPoint(const typename traits::source_type s) {
    1136:       /* // Increase degree by 0, which won't affect an existing point and will insert a new point, if necessery
    1137:          this->_cap.adjustDegree(s,0); */
    1138:       this->_cap.set.insert(typename traits::sourceRec_type(s,0));
    1139:     };
    1140:     void addCapPoint(const typename traits::sourceRec_type c) {
    1141:       this->_cap.set.insert(c);
    1142:     };
    1143:     void removeCapPoint(const typename traits::source_type s) {
    1144:       if (this->_debug) {std::cout << " Removing " << s << " from the cap" << std::endl;}
    1145:       // Clear the support and remove the point from _cap
    1146:       this->clearSupport(s);
    1147:       this->_cap.removePoint(s);
    1148:     };
    1149:     virtual void addArrow(const typename traits::source_type& p, const typename traits::target_type& q) {
    1150:       this->addArrow(p, q, typename traits::color_type());
    1151:     };
    1152:     virtual void addArrow(const typename traits::source_type& p, const typename traits::target_type& q, const typename traits::color_type& color) {
    1153:       this->addArrow(typename traits::arrow_type(p, q, color));
    1154:       //std::cout << "Added " << arrow_type(p, q, color);
    1155:     };
    1156:     virtual bool checkArrow(const typename traits::arrow_type& a) {
    1157:       if (this->_cap.set.find(a.source) == this->_cap.set.end()) return false;
    1158:       if (this->_base.set.find(a.target) == this->_base.set.end()) return false;
    1159:       return true;
    1160:     };
    1161:     virtual void addArrow(const typename traits::arrow_type& a, bool noNewPoints = false) {
    1162:       if (noNewPoints && !this->checkArrow(a)) return;
    1163:       this->_arrows.set.insert(a);
    1164:       this->addBasePoint(a.target);
    1165:       this->addCapPoint(a.source);
    1166:     };
    1167:     virtual void removeArrow(const typename traits::source_type& p, const typename traits::target_type& q) {
    1168:       this->removeArrow(typename traits::arrow_type(p, q, typename traits::color_type()));
    1169:     };
    1170:     virtual void removeArrow(const typename traits::arrow_type& a) {
    1171:       // First, produce an arrow sequence for the given source, target combination.
    1172:       typename traits::arrowSequence::traits::index_type& arrowIndex =
    1173:         ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);
    1174:       typename traits::arrowSequence::traits::index_type::iterator i,ii,j;
    1175:       i = arrowIndex.lower_bound(::boost::make_tuple(a.source,a.target));
    1176:       ii = arrowIndex.upper_bound(::boost::make_tuple(a.source, a.target));
    1177:       if (this->_debug) {
    1178:         std::cout << "removeArrow: attempting to remove arrow:" << a << std::endl;
    1179:         std::cout << "removeArrow: candidate arrows are:" << std::endl;
    1180:       }
    1181:       for(j = i; j != ii; j++) {
    1182:         if (this->_debug) {
    1183:           std::cout << " " << *j;
    1184:         }
    1185:         // Find the arrow of right color and remove it
    1186:         if(j->color == a.color) {
    1187:           if (this->_debug) {
    1188:             std::cout << std::endl << "removeArrow: found:" << *j << std::endl;
    1189:           }
    1190:           /* this->_base.adjustDegree(a.target, -1); this->_cap.adjustDegree(a.source,-1); */
    1191:           arrowIndex.erase(j);
    1192:           break;
    1193:         }
    1194:       }
    1195:     };

    1197:     void addCone(const typename traits::source_type& source, const typename traits::target_type& target){
    1198:       this->addArrow(source, target);
    1199:     }
    1200:     template<class sourceInputSequence>
    1201:     void addCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target) {
    1202:       this->addCone(sources, target, typename traits::color_type());
    1203:     }
    1204:     void addCone(const typename traits::source_type& source, const typename traits::target_type& target, const typename traits::color_type& color) {
    1205:       this->addArrow(source, target, color);
    1206:     }
    1207:     template<class sourceInputSequence>
    1208:     void
    1209:     addCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target, const typename traits::color_type& color){
    1210:       if (this->_debug > 1) {std::cout << "Adding a cone " << std::endl;}
    1211:       for(typename sourceInputSequence::iterator iter = sources->begin(); iter != sources->end(); ++iter) {
    1212:         if (this->_debug > 1) {std::cout << "Adding arrow from " << *iter << " to " << target << "(" << color << ")" << std::endl;}
    1213:         this->addArrow(*iter, target, color);
    1214:       }
    1215:     }
    1216:     void clearCone(const typename traits::target_type& t) {
    1217:       clearCone(t, typename traits::color_type(), false);
    1218:     };

    1220:     void clearCone(const typename traits::target_type& t, const typename traits::color_type&  color, bool useColor = true) {
    1221:       // Use the cone sequence types to clear the cone
    1222:       typename traits::coneSequence::traits::index_type& coneIndex =
    1223:         ::boost::multi_index::get<typename traits::coneInd>(this->_arrows.set);
    1224:       typename traits::coneSequence::traits::index_type::iterator i, ii, j;
    1225:       if (this->_debug > 20) {
    1226:         std::cout << "clearCone: removing cone over " << t;
    1227:         if(useColor) {
    1228:           std::cout << " with color" << color << std::endl;
    1229:           const Obj<typename traits::coneSequence>& cone = this->cone(t,color);
    1230:           std::cout << "[";
    1231:           for(typename traits::coneSequence::iterator ci = cone->begin(); ci != cone->end(); ci++) {
    1232:             std::cout << "  " << ci.arrow();
    1233:           }
    1234:           std::cout << "]" << std::endl;
    1235:         }
    1236:         else {
    1237:           std::cout << std::endl;
    1238:           const Obj<typename traits::coneSequence>& cone = this->cone(t);
    1239:           std::cout << "[";
    1240:           for(typename traits::coneSequence::iterator ci = cone->begin(); ci != cone->end(); ci++) {
    1241:             std::cout << "  " << ci.arrow();
    1242:           }
    1243:           std::cout << "]" << std::endl;
    1244:         }
    1245:       }
    1246:       if (useColor) {
    1247:         i = coneIndex.lower_bound(::boost::make_tuple(t,color));
    1248:         ii = coneIndex.upper_bound(::boost::make_tuple(t,color));
    1249:       } else {
    1250:         i = coneIndex.lower_bound(::boost::make_tuple(t));
    1251:         ii = coneIndex.upper_bound(::boost::make_tuple(t));
    1252:       }
    1253:       for(j = i; j != ii; j++){
    1254:         // Adjust the degrees before removing the arrow; use a different iterator, since we'll need i,ii to do the arrow erasing.
    1255:         if (this->_debug) {
    1256:           std::cout << "clearCone: adjusting degrees for endpoints of arrow: " << *j << std::endl;
    1257:         }
    1258:         /* this->_cap.adjustDegree(j->source, -1);
    1259:            this->_base.adjustDegree(j->target, -1); */
    1260:       }
    1261:       coneIndex.erase(i,ii);
    1262:     };// clearCone()

    1264:     template<class InputSequence>
    1265:     void
    1266:     restrictBase(const Obj<InputSequence>& points) {
    1267:       typename traits::baseSequence base = this->base();
    1268:       typename std::set<typename traits::target_type> remove;

    1270:       for(typename traits::baseSequence::iterator bi = base.begin(); bi != base.end(); bi++) {
    1271:         // Check whether *bi is in points, if it is NOT, remove it
    1272:         //           if (!points->contains(*bi)) {
    1273:         if (points->find(*bi) == points->end()) {
    1274:           //             this->removeBasePoint(*bi);
    1275:           remove.insert(*bi);
    1276:         }
    1277:       }
    1278:       //FIX
    1279:       for(typename std::set<typename traits::target_type>::iterator r_iter = remove.begin(); r_iter != remove.end(); ++r_iter) {
    1280:         this->removeBasePoint(*r_iter);
    1281:       }
    1282:     }

    1284:     template<class InputSequence>
    1285:     void
    1286:     excludeBase(const Obj<InputSequence>& points) {
    1287:       for(typename InputSequence::iterator pi = points->begin(); pi != points->end(); pi++) {
    1288:         this->removeBasePoint(*pi);
    1289:       }
    1290:     }

    1292:     template<class InputSequence>
    1293:     void
    1294:     restrictCap(const Obj<InputSequence>& points) {
    1295:       typename traits::capSequence cap = this->cap();
    1296:       for(typename traits::capSequence::iterator ci = cap.begin(); ci != cap.end(); ci++) {
    1297:         // Check whether *ci is in points, if it is NOT, remove it
    1298:         if(points->find(*ci) == points->end()) {
    1299:           this->removeCapPoint(*ci);
    1300:         }
    1301:       }
    1302:     }

    1304:     template<class InputSequence>
    1305:     void
    1306:     excludeCap(const Obj<InputSequence>& points) {
    1307:       for(typename InputSequence::iterator pi = points->begin(); pi != points->end(); pi++) {
    1308:         this->removeCapPoint(*pi);
    1309:       }
    1310:     }

    1312:     void clearSupport(const typename traits::source_type& s) {
    1313:       clearSupport(s, typename traits::color_type(), false);
    1314:     };
    1315:     void clearSupport(const typename traits::source_type& s, const typename traits::color_type&  color, bool useColor = true) {
    1316:       // Use the cone sequence types to clear the cone
    1317:       typename
    1318:         traits::supportSequence::traits::index_type& suppIndex = ::boost::multi_index::get<typename traits::supportInd>(this->_arrows.set);
    1319:       typename traits::supportSequence::traits::index_type::iterator i, ii, j;
    1320:       if (useColor) {
    1321:         i = suppIndex.lower_bound(::boost::make_tuple(s,color));
    1322:         ii = suppIndex.upper_bound(::boost::make_tuple(s,color));
    1323:       } else {
    1324:         i = suppIndex.lower_bound(::boost::make_tuple(s));
    1325:         ii = suppIndex.upper_bound(::boost::make_tuple(s));
    1326:       }
    1327:       for(j = i; j != ii; j++){
    1328:         // Adjust the degrees before removing the arrow
    1329:         /* this->_cap.adjustDegree(j->source, -1);
    1330:            this->_base.adjustDegree(j->target, -1); */
    1331:       }
    1332:       suppIndex.erase(i,ii);
    1333:     }
    1334:     void setCone(const typename traits::source_type& source, const typename traits::target_type& target){
    1335:       this->clearCone(target, typename traits::color_type(), false); this->addCone(source, target);
    1336:     }
    1337:     template<class sourceInputSequence>
    1338:     void setCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target) {
    1339:       this->clearCone(target, typename traits::color_type(), false); this->addCone(sources, target, typename traits::color_type());
    1340:     }
    1341:     void setCone(const typename traits::source_type& source, const typename traits::target_type& target, const typename traits::color_type& color) {
    1342:       this->clearCone(target, color, true); this->addCone(source, target, color);
    1343:     }
    1344:     template<class sourceInputSequence>
    1345:     void setCone(const Obj<sourceInputSequence>& sources, const typename traits::target_type& target, const typename traits::color_type& color){
    1346:       this->clearCone(target, color, true); this->addCone(sources, target, color);
    1347:     }
    1348:     template<class targetInputSequence>
    1349:     void addSupport(const typename traits::source_type& source, const Obj<targetInputSequence >& targets) {
    1350:       this->addSupport(source, targets, typename traits::color_type());
    1351:     }
    1352:     template<class targetInputSequence>
    1353:     void addSupport(const typename traits::source_type& source, const Obj<targetInputSequence>& targets, const typename traits::color_type& color) {
    1354:       const typename targetInputSequence::iterator end = targets->end();

    1356:       for(typename targetInputSequence::iterator iter = targets->begin(); iter != end; ++iter) {
    1357:         this->addArrow(source, *iter, color);
    1358:       }
    1359:     }
    1360:     template<typename Sifter_>
    1361:     void add(const Obj<Sifter_>& cbg, bool noNewPoints = false) {
    1362:       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);

    1364:       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) {
    1365:         this->addArrow(*a_iter, noNewPoints);
    1366:       }
    1367:       if (!noNewPoints) {
    1368:         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);

    1370:         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) {
    1371:           this->addBasePoint(*b_iter);
    1372:         }
    1373:         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);

    1375:         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) {
    1376:           this->addCapPoint(*c_iter);
    1377:         }
    1378:       }
    1379:     }
    1380:   }; // class ASifter

    1382:   // A UniSifter aka Sifter
    1383:   template <typename Source_, typename Target_, typename Color_,
    1384:             typename SupportCompare_ = ::boost::multi_index::composite_key_compare<std::less<Source_>, std::less<Color_>, std::less<Target_> >,
    1385:             typename SourceCtnr_ = SifterDef:: RecContainer<Source_, SifterDef::Rec<Source_> >, typename TargetCtnr_= SifterDef::RecContainer<Target_, SifterDef::Rec<Target_> > >
    1386:   class Sifter : public ASifter<Source_, Target_, Color_, SifterDef::uniColor, SupportCompare_, SourceCtnr_, TargetCtnr_> {
    1387:   public:
    1388:       typedef typename ASifter<Source_, Target_, Color_, SifterDef::uniColor, SupportCompare_, SourceCtnr_, TargetCtnr_>::traits       traits;
    1389:     template <typename OtherSource_, typename OtherTarget_, typename OtherColor_,
    1390:               typename OtherSupportCompare_  = ::boost::multi_index::composite_key_compare<std::less<OtherSource_>, std::less<OtherColor_>, std::less<OtherTarget_> >,
    1391:               typename OtherSourceCtnr_ = SifterDef::RecContainer<OtherSource_, SifterDef::Rec<OtherSource_> >,
    1392:               typename OtherTargetCtnr_ = SifterDef::RecContainer<OtherTarget_, SifterDef::Rec<OtherTarget_> >      >
    1393:     struct rebind {
    1394:       typedef Sifter<OtherSource_, OtherTarget_, OtherColor_, OtherSupportCompare_, OtherSourceCtnr_, OtherTargetCtnr_> type;
    1395:     };
    1396:     // Re-export some typedefs expected by CoSifter
    1397:     typedef typename traits::source_type                                            source_type;
    1398:     typedef typename traits::target_type                                            target_type;
    1399:     typedef typename traits::color_type                                             color_type;
    1400:     typedef typename traits::arrow_type                                             Arrow_;
    1401:     typedef typename traits::coneSequence                                           coneSequence;
    1402:     typedef typename traits::supportSequence                                        supportSequence;
    1403:     typedef typename traits::baseSequence                                           baseSequence;
    1404:     typedef typename traits::capSequence                                            capSequence;
    1405:     // Basic interface
    1406:     Sifter(MPI_Comm comm = PETSC_COMM_SELF, const int& debug = 0) :
    1407:       ASifter<Source_, Target_, Color_, SifterDef::uniColor, SupportCompare_, SourceCtnr_, TargetCtnr_>(comm, debug) {};

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

    1413:       const index_type& _index = ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);
    1414: #if 0
    1415:       ::boost::tuple<typename traits::source_type, typename traits::target_type> key = ::boost::make_tuple(s, t);
    1416:       typename index_type::iterator begin = _index.lower_bound(key);
    1417:       if(begin != _index.upper_bound(key)) {
    1418:         return begin->color;
    1419:       }
    1420: #else
    1421:       const typename index_type::iterator begin = _index.find(::boost::make_tuple(s, t));
    1422:       if (begin != _index.end()) {
    1423:         return begin->color;
    1424:       }
    1425: #endif
    1426: //       typename traits::arrowSequence arr(::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set), s, t);
    1427: //       if(arr.begin() != arr.end()) {
    1428: //         return arr.begin().color();
    1429: //       }
    1430:       if (fail) {
    1431:         ostringstream o;
    1432:         o << "Arrow " << s << " --> " << t << " not present";
    1433:         throw ALE::Exception(o.str().c_str());
    1434:       } else {
    1435:         static typename traits::color_type c;
    1436:         return c;
    1437:       }
    1438:     };

    1440:     template<typename ColorChanger>
    1441:     void modifyColor(const typename traits::source_type& s, const typename traits::target_type& t, const ColorChanger& changeColor) {
    1442:       typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type, typename traits::arrowInd>::type& index =
    1443:         ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);
    1444:       typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type, typename traits::arrowInd>::type::iterator i =
    1445:         index.find(::boost::make_tuple(s,t));
    1446:       if (i != index.end()) {
    1447:         index.modify(i, changeColor);
    1448:       } else {
    1449:         typename traits::arrow_type a(s, t, typename traits::color_type());
    1450:         changeColor(a);
    1451:         this->addArrow(a);
    1452:       }
    1453:     }

    1455:     struct ColorSetter {
    1456:       ColorSetter(const typename traits::color_type& color) : _color(color) {};
    1457:       void operator()(typename traits::arrow_type& p) const {
    1458:         p.color = _color;
    1459:       }
    1460:     private:
    1461:       const typename traits::color_type& _color;
    1462:     };

    1464:     void setColor(const typename traits::source_type& s, const typename traits::target_type& t, const typename traits::color_type& color) {
    1465:       ColorSetter colorSetter(color);
    1466:       typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type, typename traits::arrowInd>::type& index =
    1467:         ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);
    1468:       typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type, typename traits::arrowInd>::type::iterator i =
    1469:         index.find(::boost::make_tuple(s,t));
    1470:       if (i != index.end()) {
    1471:         index.modify(i, colorSetter);
    1472:       } else {
    1473:         typename traits::arrow_type a(s, t, color);
    1474:         this->addArrow(a);
    1475:       }
    1476:     };
    1477:     template<typename Labeling, typename AnotherSifter>
    1478:     void relabel(Labeling& relabeling, AnotherSifter& newLabel) {
    1479:       typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type, typename traits::arrowInd>::type& aInd = ::boost::multi_index::get<typename traits::arrowInd>(this->_arrows.set);

    1481:       for(typename ::boost::multi_index::index<typename traits::arrow_container_type::set_type, typename traits::arrowInd>::type::iterator a_iter = aInd.begin(); a_iter != aInd.end(); ++a_iter) {
    1482:         const typename traits::source_type newSource = relabeling.restrictPoint(a_iter->source)[0];
    1483:         const typename traits::target_type newTarget = relabeling.restrictPoint(a_iter->target)[0];

    1485:         newLabel.addArrow(newSource, newTarget);
    1486:       }
    1487:     }
    1488:   };// class Sifter

    1490:   class SifterSerializer {
    1491:   public:
    1492:     template<typename Sifter>
    1493:     static void writeSifter(std::ofstream& fs, Sifter& sifter) {
    1494:       typename Sifter::traits::arrow_container_type::set_type::size_type numArrows;

    1496:       if (sifter.commRank() == 0) {
    1497:         // Write local
    1498:         fs << sifter._arrows.set.size() << std::endl;
    1499:         for(typename Sifter::traits::arrow_container_type::set_type::iterator ai = sifter._arrows.set.begin(); ai != sifter._arrows.set.end(); ai++) {
    1500:           fs << ai->source << " " << ai->target << " " << ai->color << std::endl;
    1501:         }
    1502:         // Receive and write remote
    1503:         for(int p = 1; p < sifter.commSize(); ++p) {
    1504:           PetscInt       size;
    1505:           PetscInt      *arrows;
    1506:           MPI_Status     status;

    1509:           MPI_Recv(&size, 1, MPIU_INT, p, 1, sifter.comm(), &status);CHKERRXX(ierr);
    1510:           numArrows = size;
    1511:           fs << numArrows << std::endl;
    1512:           PetscMalloc(size*3 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
    1513:           MPI_Recv(arrows, size*3, MPIU_INT, p, 1, sifter.comm(), &status);CHKERRXX(ierr);
    1514:           for(PetscInt a = 0; a < size; ++a) {
    1515:             typename Sifter::traits::arrow_type::source_type source = arrows[a*3+0];
    1516:             typename Sifter::traits::arrow_type::target_type target = arrows[a*3+1];
    1517:             typename Sifter::traits::arrow_type::color_type  color  = arrows[a*3+2];

    1519:             fs << source << " " << target << " " << color << std::endl;
    1520:           }
    1521:           PetscFree(arrows);CHKERRXX(ierr);
    1522:         }
    1523:       } else {
    1524:         // Send remote
    1525:         PetscInt       size = sifter._arrows.set.size();
    1526:         PetscInt       a    = 0;
    1527:         PetscInt      *arrows;

    1530:         MPI_Send(&size, 1, MPIU_INT, 0, 1, sifter.comm());CHKERRXX(ierr);
    1531:         // There is no nice way to make a generic MPI type here. Really sucky
    1532:         PetscMalloc(size*3 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
    1533:         for(typename Sifter::traits::arrow_container_type::set_type::iterator ai = sifter._arrows.set.begin(); ai != sifter._arrows.set.end(); ai++, ++a) {
    1534:           arrows[a*3+0] = ai->source;
    1535:           arrows[a*3+1] = ai->target;
    1536:           arrows[a*3+2] = ai->color;
    1537:         }
    1538:         MPI_Send(arrows, size*3, MPIU_INT, 0, 1, sifter.comm());CHKERRXX(ierr);
    1539:         PetscFree(arrows);CHKERRXX(ierr);
    1540:       }
    1541:     }
    1542:     template<typename Sifter>
    1543:     static void loadSifter(std::ifstream& fs, Sifter& sifter) {
    1544:       typedef typename Sifter::traits::arrow_container_type::set_type::size_type size_type;
    1545:       if (sifter.commRank() == 0) {
    1546:         // Load local
    1547:         size_type numArrows;

    1549:         fs >> numArrows;
    1550:         for(size_type a = 0; a < numArrows; ++a) {
    1551:           typename Sifter::traits::arrow_type::source_type source;
    1552:           typename Sifter::traits::arrow_type::target_type target;
    1553:           typename Sifter::traits::arrow_type::color_type  color;

    1555:           fs >> source;
    1556:           fs >> target;
    1557:           fs >> color;
    1558:           sifter.addArrow(typename Sifter::traits::arrow_type(source, target, color));
    1559:         }
    1560:         // Load and send remote
    1561:         for(int p = 1; p < sifter.commSize(); ++p) {
    1562:           PetscInt       size;
    1563:           PetscInt      *arrows;

    1566:           fs >> numArrows;
    1567:           size = numArrows;
    1568:           MPI_Send(&size, 1, MPIU_INT, p, 1, sifter.comm());CHKERRXX(ierr);
    1569:           PetscMalloc(size*3 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
    1570:           for(PetscInt a = 0; a < size; ++a) {
    1571:             typename Sifter::traits::arrow_type::source_type source;
    1572:             typename Sifter::traits::arrow_type::target_type target;
    1573:             typename Sifter::traits::arrow_type::color_type  color;

    1575:             fs >> source;
    1576:             fs >> target;
    1577:             fs >> color;
    1578:             arrows[a*3+0] = source;
    1579:             arrows[a*3+1] = target;
    1580:             arrows[a*3+2] = color;
    1581:           }
    1582:           MPI_Send(arrows, size*3, MPIU_INT, p, 1, sifter.comm());CHKERRXX(ierr);
    1583:           PetscFree(arrows);CHKERRXX(ierr);
    1584:         }
    1585:       } else {
    1586:         // Load remote
    1587:         PetscInt       size;
    1588:         PetscInt      *arrows;
    1589:         MPI_Status     status;

    1592:         MPI_Recv(&size, 1, MPIU_INT, 0, 1, sifter.comm(), &status);CHKERRXX(ierr);
    1593:         PetscMalloc(size*3 * sizeof(PetscInt), &arrows);CHKERRXX(ierr);
    1594:         MPI_Recv(arrows, size*3, MPIU_INT, 0, 1, sifter.comm(), &status);CHKERRXX(ierr);
    1595:         for(PetscInt a = 0; a < size; ++a) {
    1596:           sifter.addArrow(typename Sifter::traits::arrow_type(arrows[a*3+0], arrows[a*3+1], arrows[a*3+2]));
    1597:         }
    1598:         PetscFree(arrows);CHKERRXX(ierr);
    1599:       }
    1600:     }
    1601:   };
    1602: } // namespace ALE

    1604: #endif